Mcgpj Class Reference

#include <Mcgpj.h>

List of all members.

Public Member Functions

 Mcgpj (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()

Private Attributes

double cmE
double pc
double de
double nt0
double dt
double dp
double at
double td
double am
double cm
double em
double ti
double al
double thm
double thp
double te
double re
double spread
double phase
int proc
int NRad
int IsHardPhoton
int IsNoVacPol
int IsFSR
int Seed
double fM
int fpid [2]
int m_fmode5pi
std::string m_datapath
std::string m_vpolfname
IBesRndmGenSvc * p_BesRndmGenSvc

Detailed Description

Definition at line 20 of file Mcgpj.h.

Constructor & Destructor Documentation

Mcgpj::Mcgpj ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Definition at line 54 of file Mcgpj.cxx.

References al, am, at, cm, cmE, de, dp, dt, em, IsFSR, IsHardPhoton, IsNoVacPol, m_datapath, m_fmode5pi, m_vpolfname, NRad, nt0, pc, proc, re, spread, td, te, thm, thp, and ti.

00054                                                            :Algorithm(name, pSvcLocator){
00055   declareProperty("Data", m_datapath = "${MCGPJROOT}/data"); // Path to data files
00056   declareProperty("VpolFname", m_vpolfname = "vpol.dat"); // Vacuum polarization data file name
00057   declareProperty("CMEnergy", cmE = 3.097); // 2*Ebeam [GeV]
00058   declareProperty("Process", proc = 0); // process
00059   declareProperty("NRad", NRad = 20000);
00060   declareProperty("HardPhoton", IsHardPhoton = 1);
00061   declareProperty("NoVacPol", IsNoVacPol = 0);
00062   declareProperty("FSR", IsFSR = 1);
00063   declareProperty("AcolAngle", pc  = 0.);
00064   declareProperty("DeltaE", de  = -1.); // [GeV]
00065   declareProperty("NTheta0", nt0 = -1.); // [1/sqrt(gamma)]
00066   declareProperty("DeltaTheta", dt  = 0.5); //[rad]
00067   declareProperty("DeltaPhi", dp  = 0.3); // [rad]
00068   declareProperty("AverageTheta", at  = 1.1); //[rad]
00069   declareProperty("ThetaDetector", td  = 1.1 - 0.5/2); //[rad]
00070   declareProperty("AverageMomentum", am  = 0.090); //[GeV/c]
00071   declareProperty("CrossMomentum", cm  = 0.090); //[GeV/c]
00072   declareProperty("MinimumEnergy", em  = 0.050);//[GeV]
00073   declareProperty("ThetaIntermediate", ti  = 0.473);// [rad]
00074   declareProperty("AlphaScale", al  = 1.); 
00075   declareProperty("ThetaMinus", thm = -1.);  //[rad]
00076   declareProperty("ThetaPlus", thp = -1.);  //[rad]
00077   declareProperty("AbsoluteError", te  = 0.05); // absolute error in [nb]
00078   declareProperty("RelativeError", re  = 0.05); // realative error
00079 //  declareProperty("InitialSeed", Seed = 0);  
00080   declareProperty("BeamSpread", spread  = -1);// [MeV]
00081 //  declareProperty("PhaseShift", phase  = -90);// [deg]
00082   declareProperty("Mode5pi", m_fmode5pi  = 0);// 0 - via b1+/- pi-/+, 1 via b10 pi0
00083 }

Member Function Documentation

StatusCode Mcgpj::execute (  ) 

Definition at line 318 of file Mcgpj.cxx.

References cmE, CS, calibUtil::ERROR, fM, fpid, TCrossPart::GenUnWeightedEvent(), TCrossPart::GetNfinal(), TCrossPart::GetParticles(), TCrossPart::GetPid(), gRad, genRecEmupikp::i, Bes_Common::INFO, TRadCor::MakeEvent(), mass, pid, proc, q, and Bes_Common::WARNING.

00318                          { 
00319   MsgStream log(messageService(), name());
00320   log << MSG::INFO << "Mcgpj executing" << endreq;
00321         log<<MSG::WARNING<<"execute start"<<endreq;
00322 
00323   // Fill event information
00324   GenEvent* evt = new GenEvent(1,1);
00325 
00326   GenVertex* prod_vtx = new GenVertex();
00327 
00328   evt->add_vertex( prod_vtx );
00329 
00330   // incoming beam e+
00331   GenParticle* p = 
00332     new GenParticle(HepLorentzVector(0,0,0.5*cmE*1e-3,0.5*cmE*1e-3),-11, 3);
00333   p->suggest_barcode(1);
00334   prod_vtx->add_particle_in(p);
00335 
00336   // incoming beam e-
00337   p = 
00338     new GenParticle(HepLorentzVector(0,0,0.5*cmE*1e-3, 0.5*cmE*1e-3), 11, 3);
00339   p->suggest_barcode(2);
00340   prod_vtx->add_particle_in(p);
00341   
00342   int npart = 2;
00343   if(proc<100){
00344     double mom[4*6];
00345     int np;
00346     gRad->MakeEvent(mom,np);
00347     
00348     // final particle 1
00349     for(int i=0; i<np;i++){
00350       double ptot = mom[i*4+3];
00351       double px = ptot*mom[i*4+0];
00352       double py = ptot*mom[i*4+1];
00353       double pz = ptot*mom[i*4+2];
00354       double mass = 0;    
00355       int pid = 22;
00356       if(i<2){
00357                 pid = fpid[i];
00358         mass = fM;
00359       }
00360       double etot = sqrt(ptot*ptot + mass*mass*1e-6);
00361       p = new GenParticle( HepLorentzVector(px,py,pz,etot), pid, 1); 
00362       p->suggest_barcode(i+3);
00363       prod_vtx->add_particle_out(p);
00364       npart++;
00365     }
00366   } else {
00367     int ipart = CS->GenUnWeightedEvent();
00368     size_t nmax = CS->GetNfinal() + ((ipart==0)?1:2);
00369     for(size_t i=0;i<nmax;i++){
00370       TLorentzVector &q = *(CS->GetParticles()[i]);
00371       int pid = CS->GetPid(i);
00372       p = 
00373         new GenParticle( HepLorentzVector(q.X()*1e-3,q.Y()*1e-3,q.Z()*1e-3,q.T()*1e-3), pid, 1);
00374       p->suggest_barcode(i+3);
00375       prod_vtx->add_particle_out(p);
00376       npart++;
00377     }
00378   }
00379 
00380   if( log.level() < MSG::INFO ){
00381     evt->print();  
00382   }
00383 
00384   // Check if the McCollection already exists
00385   SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00386   if (anMcCol!=0){
00387     // Add event to existing collection
00388         log<<MSG::WARNING<<"add event"<<endreq;
00389     MsgStream log(messageService(), name());
00390     log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00391     McGenEvent* mcEvent = new McGenEvent(evt);
00392     anMcCol->push_back(mcEvent);
00393   } else {
00394     // Create Collection and add  to the transient store
00395         log<<MSG::WARNING<<"create collection"<<endreq;
00396     McGenEventCol *mcColl = new McGenEventCol;
00397     McGenEvent* mcEvent = new McGenEvent(evt);
00398     mcColl->push_back(mcEvent);
00399     StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00400     if (sc != StatusCode::SUCCESS) {
00401       log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00402       delete mcColl;
00403       delete evt;
00404       delete mcEvent;
00405       return StatusCode::FAILURE;
00406     } else {
00407       log << MSG::INFO << "McGenEventCol created and " << npart 
00408           <<" particles stored in McGenEvent" << endreq;
00409     }
00410   }
00411 
00412         log<<MSG::WARNING<<"execute end"<<endreq;
00413   return StatusCode::SUCCESS;
00414 }

StatusCode Mcgpj::finalize (  ) 

Definition at line 416 of file Mcgpj.cxx.

References gConst, gCut, gGlobal, gRad, Bes_Common::INFO, and MatrixElements.

00416                           {
00417   MsgStream log(messageService(), name());
00418 
00419   delete gRad;
00420   delete MatrixElements;
00421   delete gRandom;
00422   delete gConst;
00423   delete gGlobal;
00424   delete gCut;
00425 
00426   log << MSG::INFO << "Mcgpj finalized" << endreq;
00427 
00428   return StatusCode::SUCCESS;
00429 }

StatusCode Mcgpj::initialize (  ) 

Definition at line 85 of file Mcgpj.cxx.

References al, am, at, TKinemCut::AverageTheta(), cm, cmE, cos(), TKinemCut::CosPsi(), CS, de, TKinemCut::DeltaPhi(), TKinemCut::DeltaTheta(), dp, dt, EBeam, em, TKinemCut::EMin(), calibUtil::ERROR, fM, fpid, gConst, gCut, TRadGlobal::Get_RelativeError(), TRadGlobal::Get_TotalError(), IBesRndmGenSvc::GetEngine(), gGlobal, gRad, Bes_Common::INFO, TRadCor::Init(), TKinemCut::Init(), TRadGlobal::Init(), InList, IsFSR, IsHardPhoton, IsNoVacPol, ganga-rec::j, m_datapath, m_fmode5pi, m_vpolfname, TRadCor::MakeCrossSection(), TCrossPart::MakeParts(), MatrixElements, MaxList, NRad, nt0, p_BesRndmGenSvc, TKinemCut::PAverage(), pc, TKinemCut::PCross(), TKinemCut::Print(), TRadGlobal::Print(), TConstants::Print(), proc, re, TRadGlobal::Set_dE_Abs(), TRadGlobal::Set_E(), TRadGlobal::Set_RelativeError(), TRadGlobal::Set_Theta0_Rel(), TRadGlobal::Set_ThetaInt(), TRadGlobal::Set_ThetaMin(), TRadGlobal::Set_TotalError(), TRadGlobal::Set_Type(), TConstants::SetAlphaScale(), TCrossPart::SetBeamSpread(), TRadGlobal::SetDatadir(), TVCrossPart::SetHardPhoton(), TRadCor::SetNEvents(), TVCrossPart::SetNoFSR(), TRadCor::SetPartList(), TRadGlobal::SetVpolFname(), TVCrossPart::SetZeroVP(), spread, td, te, TKinemCut::ThetaMinM(), TKinemCut::ThetaMinP(), thm, thp, and ti.

00085                             {
00086   MsgStream log(messageService(), name());
00087   
00088   log << MSG::INFO << "Mcgpj initialize" << endreq;
00089 
00090   static const bool CREATEIFNOTTHERE(true);                                                            
00091   StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);                 
00092   if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)                                                 
00093   {                                                                                                    
00094     log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;                      
00095     return RndmStatus;                                                                                 
00096   }                                                                                                    
00097   CLHEP::HepRandomEngine* engine  = p_BesRndmGenSvc->GetEngine("Mcgpj");                            
00098   engine->showStatus();
00099 
00100   // GeV -> MeV
00101   cmE *= 1e3;
00102   de *= 1e3;
00103   am *= 1e3;
00104   cm *= 1e3;
00105   em *= 1e3;
00106 
00107   if(gRandom) delete gRandom;
00108   gRandom = new TRandom3();
00109   gRandom->SetSeed(engine->getSeed());
00110   //gRandom->Dump();
00111 
00112   if (proc<100){
00113 
00114     gConst  = new TConstants();
00115     gGlobal = new TRadGlobal();
00116     gCut    = new TKinemCut();
00117 
00118     if((proc%10)==3){
00119       td = 0.025;
00120       dt = 10;
00121       dp = 10;
00122       cm = 0;
00123       am = 0;
00124     }
00125   
00126     const int MaxList = 20;
00127     bool InList[MaxList];
00128     for(int j = 0; j<MaxList; j++) InList[j] = true;
00129 
00130     double EBeam = 0.5*cmE;
00131 
00132     if( ( EBeam < 100 || EBeam > 2500 ) && !IsNoVacPol ){
00133       log<<MSG::ERROR<<"Invalid value of beam energy:"<<EBeam<<endreq;
00134       return StatusCode::FAILURE;
00135     }
00136 
00137     gGlobal->Set_TotalError(te);
00138   
00139     gGlobal->Set_RelativeError(re);
00140 
00141     gGlobal->Set_Type((proc%10));
00142 
00143     gGlobal->Set_E(EBeam);
00144     
00145     gGlobal->Set_ThetaInt(ti);
00146 
00147     if(0>td)
00148       gGlobal->Set_ThetaMin(at-dt/2);
00149     else
00150       gGlobal->Set_ThetaMin(td);
00151 
00152     if(0>de){
00153       if((proc%10)==0) 
00154         gGlobal->Set_dE_Abs(0.015*EBeam);
00155       else 
00156         gGlobal->Set_dE_Abs(0.0003*EBeam);
00157     }
00158     else
00159       gGlobal->Set_dE_Abs(de);
00160 
00161     if(0>nt0){
00162       if(proc>10) 
00163         gGlobal->Set_Theta0_Rel(0.0);
00164       else
00165         gGlobal->Set_Theta0_Rel(1.5);
00166     }else
00167       gGlobal->Set_Theta0_Rel(nt0);
00168 
00169     if(0<thm){
00170       gCut->ThetaMinM(thm);
00171       gGlobal->Set_ThetaMin(thm);
00172     }
00173 
00174     if(0<thp)
00175       gCut->ThetaMinP(thp);
00176 
00177     gCut->CosPsi(cos(pc));
00178 
00179     gCut->DeltaTheta(dt);
00180   
00181     gCut->AverageTheta(at);
00182 
00183     gCut->DeltaPhi(dp);
00184   
00185     gCut->PAverage(am);
00186     
00187     gCut->PCross(cm);
00188 
00189     gCut->EMin(em);
00190   
00191     gConst->SetAlphaScale(al);
00192 
00193     gConst->Print();
00194 
00195     try{
00196       gGlobal->Init();
00197     }catch(std::logic_error lge){
00198       log<<MSG::ERROR<<lge.what()<<endreq;
00199       return StatusCode::FAILURE;
00200     }
00201 
00202     gGlobal->SetDatadir(m_datapath);
00203     gGlobal->SetVpolFname(m_vpolfname);
00204 //    gGlobal->SetIntFname("");
00205 
00206     gGlobal->Print();
00207 
00208     gCut->Init();
00209     gCut->Print();
00210 
00211     log<<MSG::INFO<<"Cross-section statistical precision will be better than "
00212        <<gGlobal->Get_TotalError()<<" nb and "
00213        <<gGlobal->Get_RelativeError()*100<<"%"<<endreq;
00214 
00215     if(!IsHardPhoton)
00216       log<<MSG::INFO<<"Hard photon on big angle is not included!"<<endreq;
00217     
00218     if(IsNoVacPol)
00219       log<<MSG::INFO<<"Vacuum polarization is not included!"<<endreq;
00220 
00221     if(!IsFSR)
00222       log<<MSG::INFO<<"Final state radiation is not included!"<<endreq;
00223 
00224     if(proc>10)
00225       log<<MSG::INFO<<"Alpha order generation only!"<<endreq;
00226 
00227     log<<std::flush;
00228 
00229     if(proc==0){
00230       MatrixElements = new TEPCrossPart();
00231       InList[18] = false;
00232     }
00233     else if(proc==1 || proc==5)
00234       MatrixElements = new TMuCrossPart();
00235     else if(proc==2)
00236       MatrixElements = new TPiCrossPart();
00237     else if(proc==3)
00238       MatrixElements = new TKnCrossPart();
00239     else if(proc==4)
00240       MatrixElements = new TKcCrossPart();
00241     else if(proc==6)
00242       MatrixElements = new TGGCrossPart();
00243     else if(proc==10){
00244       MatrixElements = new TEPCrossPart();
00245       for(int j = 0; j<MaxList; j++) InList[j] = false;
00246       InList[16] = true;
00247       InList[17] = true;
00248       InList[18] = true;
00249     }
00250     else
00251       return StatusCode::FAILURE;
00252 
00253     if(IsNoVacPol)
00254       MatrixElements->SetZeroVP();
00255 
00256     if(!IsFSR)
00257       MatrixElements->SetNoFSR();
00258 
00259     gRad = new TRadCor(MatrixElements);
00260     gRad->SetNEvents(NRad);
00261     gRad->SetPartList(InList);
00262     gRad->Init();
00263 
00264     MatrixElements->SetHardPhoton(IsHardPhoton);
00265     gRad->MakeCrossSection();
00266     if((proc%10)==2)((TPiCrossPart*)MatrixElements)->GetFormFactor()->Print();
00267 
00268     if((proc%10)==0){
00269       fpid[0] = 11; fpid[1] = -11; fM = 0.51099891;
00270     }
00271     if((proc%10)==1){
00272       fpid[0] = 13; fpid[1] = -13;  fM = 105.658367;
00273     }
00274     if((proc%10)==2){
00275       fpid[0] = 211; fpid[1] = -211; fM = 139.57018;
00276     }
00277     if((proc%10)==3){
00278       fpid[0] = 130; fpid[1] = 310; fM = 497.614;
00279     }
00280     if((proc%10)==4){
00281       fpid[0] = 321; fpid[1] = -321; fM = 493.677;
00282     }
00283     if((proc%10)==5){
00284       fpid[0] = 15; fpid[1] = -15;  fM = 1776.84;
00285     }
00286     if((proc%10)==6){
00287       fpid[0] = 22; fpid[1] = 22; fM = 0;
00288     }
00289   }else{
00290     double EBeam = 0.5*cmE;
00291     if(0>de) de = 0.005*EBeam;
00292     
00293     if(0>nt0) nt0 = 1;
00294 
00295     if(proc==100)
00296     {
00297       CS = new T3piCrossPart(EBeam,de,nt0);
00298       if(spread>0) CS->SetBeamSpread(spread);
00299     } 
00300     else if( proc==101)
00301     { 
00302       CS = new T4piCrossPart(EBeam,de,nt0);
00303       if(spread>0) CS->SetBeamSpread(spread);   
00304     }
00305     else if( proc==102)
00306     {
00307       CS = new T5piCrossPart(EBeam,de,nt0,m_fmode5pi);
00308       if(spread>0) CS->SetBeamSpread(spread);
00309     }
00310     CS->MakeParts(te);
00311   }
00312   
00313   log << MSG::INFO << "Mcgpj initialize finished" << endreq;
00314 
00315   return StatusCode::SUCCESS;
00316 }

Member Data Documentation

double Mcgpj::al [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::am [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::at [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::cm [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::cmE [private]

Definition at line 29 of file Mcgpj.h.

Referenced by execute(), initialize(), and Mcgpj().

double Mcgpj::de [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::dp [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::dt [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::em [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::fM [private]

Definition at line 33 of file Mcgpj.h.

Referenced by execute(), and initialize().

int Mcgpj::fpid[2] [private]

Definition at line 34 of file Mcgpj.h.

Referenced by execute(), and initialize().

int Mcgpj::IsFSR [private]

Definition at line 32 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

int Mcgpj::IsHardPhoton [private]

Definition at line 32 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

int Mcgpj::IsNoVacPol [private]

Definition at line 32 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

std::string Mcgpj::m_datapath [private]

Definition at line 36 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

int Mcgpj::m_fmode5pi [private]

Definition at line 35 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

std::string Mcgpj::m_vpolfname [private]

Definition at line 37 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

int Mcgpj::NRad [private]

Definition at line 32 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::nt0 [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

IBesRndmGenSvc* Mcgpj::p_BesRndmGenSvc [private]

Definition at line 39 of file Mcgpj.h.

Referenced by initialize().

double Mcgpj::pc [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::phase [private]

Definition at line 31 of file Mcgpj.h.

int Mcgpj::proc [private]

Definition at line 32 of file Mcgpj.h.

Referenced by execute(), initialize(), and Mcgpj().

double Mcgpj::re [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

int Mcgpj::Seed [private]

Definition at line 32 of file Mcgpj.h.

double Mcgpj::spread [private]

Definition at line 31 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::td [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::te [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::thm [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::thp [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

double Mcgpj::ti [private]

Definition at line 30 of file Mcgpj.h.

Referenced by initialize(), and Mcgpj().

Generated on Tue Nov 29 23:20:06 2016 for BOSS_7.0.2 by  doxygen 1.4.7