00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "Mcgpj/Mcgpj.h"
00014
00015
00016 #include "HepMC/GenEvent.h"
00017 #include "HepMC/GenVertex.h"
00018 #include "HepMC/GenParticle.h"
00019
00020 #include "CLHEP/Vector/LorentzVector.h"
00021 #include "BesKernel/IBesRndmGenSvc.h"
00022
00023 #include "GaudiKernel/SmartDataPtr.h"
00024
00025 #include "GeneratorObject/McGenEvent.h"
00026
00027 #include "TRandom3.h"
00028 #include "TRadCor.h"
00029 #include "TEPCrossPart.h"
00030 #include "TMuCrossPart.h"
00031 #include "TPiCrossPart.h"
00032 #include "TKnCrossPart.h"
00033 #include "TKcCrossPart.h"
00034 #include "TGGCrossPart.h"
00035 #include "TRadGlobal.h"
00036 #include "TKinemCut.h"
00037
00038 #include "TPhoton_o.h"
00039 #include "TDFun_o.h"
00040 #include "T5piCrossPart.h"
00041 #include "T4piCrossPart.h"
00042 #include "T3piCrossPart.h"
00043
00044 using namespace CLHEP;
00045
00046 TRadCor *gRad;
00047 TRadGlobal *gGlobal;
00048 TKinemCut *gCut;
00049 TConstants *gConst;
00050 TVCrossPart *MatrixElements;
00051 TCrossPart *CS;
00052
00053
00054 Mcgpj::Mcgpj(const std::string& name, ISvcLocator* pSvcLocator):Algorithm(name, pSvcLocator){
00055 declareProperty("Data", m_datapath = "${MCGPJROOT}/data");
00056 declareProperty("VpolFname", m_vpolfname = "vpol.dat");
00057 declareProperty("CMEnergy", cmE = 3.097);
00058 declareProperty("Process", proc = 0);
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.);
00065 declareProperty("NTheta0", nt0 = -1.);
00066 declareProperty("DeltaTheta", dt = 0.5);
00067 declareProperty("DeltaPhi", dp = 0.3);
00068 declareProperty("AverageTheta", at = 1.1);
00069 declareProperty("ThetaDetector", td = 1.1 - 0.5/2);
00070 declareProperty("AverageMomentum", am = 0.090);
00071 declareProperty("CrossMomentum", cm = 0.090);
00072 declareProperty("MinimumEnergy", em = 0.050);
00073 declareProperty("ThetaIntermediate", ti = 0.473);
00074 declareProperty("AlphaScale", al = 1.);
00075 declareProperty("ThetaMinus", thm = -1.);
00076 declareProperty("ThetaPlus", thp = -1.);
00077 declareProperty("AbsoluteError", te = 0.05);
00078 declareProperty("RelativeError", re = 0.05);
00079
00080 declareProperty("BeamSpread", spread = -1);
00081
00082 declareProperty("Mode5pi", m_fmode5pi = 0);
00083 }
00084
00085 StatusCode Mcgpj::initialize(){
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
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
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
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 }
00317
00318 StatusCode Mcgpj::execute(){
00319 MsgStream log(messageService(), name());
00320 log << MSG::INFO << "Mcgpj executing" << endreq;
00321 log<<MSG::WARNING<<"execute start"<<endreq;
00322
00323
00324 GenEvent* evt = new GenEvent(1,1);
00325
00326 GenVertex* prod_vtx = new GenVertex();
00327
00328 evt->add_vertex( prod_vtx );
00329
00330
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
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
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
00385 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00386 if (anMcCol!=0){
00387
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
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 }
00415
00416 StatusCode Mcgpj::finalize(){
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 }