00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "EvtDecay.h"
00030 #include "ReadME.h"
00031 #include "EvtGenBase/EvtVector4R.hh"
00032 #include "EvtGenBase/EvtParticle.hh"
00033 #include "EvtGenBase/EvtParticleFactory.hh"
00034 #include "EvtGenBase/EvtDecayTag.hh"
00035 #include "EvtGen.hh"
00036 #include "EvtGenBase/EvtRandomEngine.hh"
00037 #include "EvtGenBase/EvtDecayTable.hh"
00038 #include "EvtGenModels/EvtLundCharm.hh"
00039 #include "EvtGenModels/EvtPsi3Sdecay.hh"
00040 #include "EvtGlobalSet.hh"
00041 #include "EvtGenBase/EvtHelSys.hh"
00042 #include "EvtGenModels/EvtCalHelAmp.hh"
00043 #include "EvtGenModels/EvtConExc.hh"
00044
00045
00046 #include "HepMC/GenEvent.h"
00047 #include "HepMC/GenVertex.h"
00048 #include "HepMC/GenParticle.h"
00049 #include "EventModel/EventHeader.h"
00050 #include "GaudiKernel/SmartDataPtr.h"
00051 #include "DataInfoSvc/IDataInfoSvc.h"
00052 #include "DataInfoSvc/DataInfoSvc.h"
00053 #include "GaudiKernel/MsgStream.h"
00054 #include "GaudiKernel/ISvcLocator.h"
00055 #include "GaudiKernel/AlgFactory.h"
00056 #include "GaudiKernel/DataSvc.h"
00057 #include "GaudiKernel/SmartDataPtr.h"
00058 #include "GaudiKernel/IPartPropSvc.h"
00059
00060 #include "BesKernel/IBesRndmGenSvc.h"
00061 #include "GeneratorObject/McGenEvent.h"
00062 #include "CLHEP/Random/Ranlux64Engine.h"
00063 #include "CLHEP/Random/RandFlat.h"
00064
00065 #include <stdlib.h>
00066 #include <string.h>
00067
00068 static string mydecayrec;
00069 static double _amps_max;
00070 int nchr = 0;
00071 int nchr_e = 0;
00072 int nchr_mu= 0;
00073 int nchr_pi= 0;
00074 int nchr_k = 0;
00075 int nchr_p = 0;
00076 int N_gamma=0;
00077 int N_gammaFSR = 0;
00078 int fst[50];
00079 int EvtDecay::m_runNo=0;
00080 double EvtConExc::SetMthr=0;
00081
00082 EvtDecay::EvtDecay(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
00083 {
00084
00085
00086
00087
00088
00089
00090 declareProperty("DecayDecDir", m_DecayDec = "");
00091 declareProperty("PdtTableDir", m_PdtTable = "");
00092 declareProperty("userDecayTableName", userDecFileName = "NONE");
00093 declareProperty("DecayTopology", m_DecayTop = "");
00094 declareProperty("DecayRecord", m_DecayRec = "");
00095 declareProperty("ParentParticle", m_ParentPart = "NONE");
00096 declareProperty("Multiplicity", m_Ncharge = false);
00097 declareProperty("NutupleFile",m_NtupleFile = false);
00098 declareProperty("mDIY",_mDIY= false);
00099 declareProperty("FDPdata",m_SB3File = "");
00100 declareProperty("FDPHisT",m_SB3HT = "");
00101 declareProperty("FDPpart",m_FDPparticle ="");
00102 declareProperty("TruthFile",m_truthFile ="");
00103 declareProperty("TruthPart",m_truthPart ="");
00104 declareProperty("Psi3SopenCharm",m_Psi4040OpenCharm=false);
00105 declareProperty("Psi2openCharm", m_Psi2openCharm=false);
00106 declareProperty("SetMthrForConExc",m_SetMthr=0.0);
00107 declareProperty("statDecays",m_statDecays=false);
00108 declareProperty("TagLundCharmModel", m_tagLundModel=false);
00109 m_mystring.clear();
00110 declareProperty("FileForTrackGen", m_mystring);
00111
00112 m_polarization.clear();
00113 declareProperty("polarization", m_polarization);
00114 declareProperty("eBeamPolarization", m_eBeamPolarization=-1);
00115 m_cluster0.clear();m_wind0.clear();
00116 m_cluster1.clear();m_wind1.clear();
00117 m_cluster2.clear();m_wind2.clear();
00118 declareProperty("cluster0",m_cluster0);
00119 declareProperty("cluster1",m_cluster1);
00120 declareProperty("cluster2",m_cluster2);
00121 declareProperty("masswin0",m_wind0);
00122 declareProperty("masswin1",m_wind1);
00123 declareProperty("masswin2",m_wind2);
00124 declareProperty("OutputP4",m_outputp4="");
00125 declareProperty("ReadMeasuredEcms", m_RdMeasuredEcms = false);
00126 m_ReadTruth.clear();
00127 declareProperty("ReadTruth",m_ReadTruth);
00128
00129
00130 }
00131
00132
00133 StatusCode EvtDecay::initialize(){
00134
00135 MsgStream log(messageService(), name());
00136 if(m_truthFile!=""){truth.open(m_truthFile.c_str());}
00137 log << MSG::INFO << "EvtDecay initialize" << endreq;
00138 static const bool CREATEIFNOTTHERE(true);
00139 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00140 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00141 {
00142 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00143 return RndmStatus;
00144 }
00145
00146 EvtGlobalSet::SV.clear();
00147 for(int i=0;i<m_mystring.size();i++){EvtGlobalSet::SV.push_back(m_mystring[i]);}
00148
00149 EvtGlobalSet::iVV.clear();
00150 EvtGlobalSet::dVV.clear();
00151 std::vector<std::vector<int> >myivv;
00152 std::vector<std::vector<double> >mydvv;
00153
00154 EvtConExc::SetMthr=m_SetMthr;
00155
00156 myivv.push_back(m_cluster0);
00157 myivv.push_back(m_cluster1);
00158 myivv.push_back(m_cluster2);
00159
00160 mydvv.push_back(m_wind0);
00161 mydvv.push_back(m_wind1);
00162 mydvv.push_back(m_wind2);
00163
00164 for(int i=0;i<myivv.size();i++){
00165 std::vector<int> theivv;
00166 for(int j=0;j<myivv[i].size();j++){
00167 theivv.push_back(myivv[i][j]);
00168 }
00169 if(theivv.size()) EvtGlobalSet::iVV.push_back(theivv);
00170 }
00171
00172 for(int i=0;i<mydvv.size();i++){
00173 std::vector<double> thedvv;
00174 for(int j=0;j<mydvv[i].size();j++){
00175 thedvv.push_back(mydvv[i][j]);
00176 }
00177 if(thedvv.size()) EvtGlobalSet::dVV.push_back(thedvv);
00178 }
00179
00180
00181 if(m_eBeamPolarization>1) {std::cout<<"The e-beam polaziation should be in 0 to 1"<<std::endl;abort();}
00182 m_numberEvent=0;
00183 first = true;
00184 m_ampscalflag = true;
00185
00186 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("EVTGEN");
00187 const long s = engine->getSeed();
00188 p_BesRndmGenSvc->setGenseed(s+1);
00189
00190 m_seeds.clear();
00191 m_seeds.push_back(s);
00192 totalChannels = 0;
00193
00194
00195
00196 if ( m_DecayDec == "") {
00197 m_DecayDec=getenv("BESEVTGENROOT");
00198 m_DecayDec +="/share/DECAY.DEC";
00199 }
00200
00201 if ( m_PdtTable == "") {
00202 m_PdtTable =getenv("BESEVTGENROOT");
00203 m_PdtTable +="/share/pdt.table";
00204 }
00205
00206 char catcmd[300];
00207 std::cout<<"===================== user decay card ========================"<<std::endl;
00208 strcpy(catcmd, "cat ");
00209 strcat(catcmd, userDecFileName.c_str());
00210 system(catcmd);
00211
00212 IDataInfoSvc *tmpInfoSvc;
00213 DataInfoSvc* dataInfoSvc;
00214
00215
00216 StatusCode status;
00217 status = service("DataInfoSvc",tmpInfoSvc);
00218 if (status.isSuccess()) {
00219 log << MSG::INFO << "got the DataInfoSvc" << endreq;
00220 dataInfoSvc=dynamic_cast<DataInfoSvc *>(tmpInfoSvc);
00221 dataInfoSvc->setDecayCard(userDecFileName);
00222 } else {
00223 log << MSG::ERROR << "could not get the DataInfoSvc" << endreq;
00224 return StatusCode::FAILURE;
00225 }
00226
00227
00228
00229 m_RandomEngine = new EvtBesRandom(engine);
00230 m_Gen = new EvtGen( m_DecayDec.c_str(), m_PdtTable.c_str(), m_RandomEngine );
00231
00232 if(userDecFileName =="NONE") log << MSG::INFO << "EvtDecay User did not define his Decay table EvtGen will use standart" << endreq;
00233 if(!(userDecFileName =="NONE")) m_Gen->readUDecay(userDecFileName.c_str());
00234
00235 if(!(m_DecayTop==" ")) {
00236 outfile.open(m_DecayTop.c_str());
00237 }
00238
00239
00240
00241
00242 if(m_NtupleFile) {
00243 NTuplePtr nt1(ntupleSvc(),"MYROOTFILE/Trk");
00244 if(nt1) {m_tuple=nt1;}
00245 else {
00246 m_tuple = ntupleSvc()->book ("MYROOTFILE/Trk", CLID_ColumnWiseTuple, "Generator-trk-Ntuple");
00247 if(m_tuple){
00248 status = m_tuple->addItem ("TotNumTrk", TotNumTrk, 0,100);
00249 status = m_tuple->addIndexedItem ("Trk_index", TotNumTrk, m_Trk_index);
00250 status = m_tuple->addIndexedItem ("m_px_trk", TotNumTrk, m_px_trk);
00251 status = m_tuple->addIndexedItem ("m_py_trk", TotNumTrk, m_py_trk);
00252 status = m_tuple->addIndexedItem ("m_pz_trk", TotNumTrk, m_pz_trk);
00253 status = m_tuple->addIndexedItem ("m_en_trk", TotNumTrk, m_en_trk);
00254 status = m_tuple->addIndexedItem ("FST", TotNumTrk, m_fst);
00255
00256 status = m_tuple->addItem ("nchr", m_nchr);
00257 status = m_tuple->addItem ("nchr_e", m_nchr_e);
00258 status = m_tuple->addItem ("nchr_mu", m_nchr_mu);
00259 status = m_tuple->addItem ("nchr_pi", m_nchr_pi);
00260 status = m_tuple->addItem ("nchr_k", m_nchr_k);
00261 status = m_tuple->addItem ("nchr_p", m_nchr_p);
00262 status = m_tuple->addItem ("N_gamma", m_gamma);
00263 status = m_tuple->addItem ("N_gammaFSR", m_gammaFSR);
00264 status = m_tuple->addItem ("Flag1", m_flag1);
00265 } else {
00266 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
00267 return StatusCode::FAILURE;
00268 }
00269 }
00270
00271 NTuplePtr nt2(ntupleSvc(),"MYROOTFILE/mass");
00272 if(nt2) {mass_tuple=nt2;}
00273 else {
00274 mass_tuple = ntupleSvc()->book ("MYROOTFILE/mass", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
00275 if(mass_tuple){
00276 status = mass_tuple->addItem ("m12", m_m12);
00277 status = mass_tuple->addItem ("m13", m_m13);
00278 status = mass_tuple->addItem ("m23", m_m23);
00279 status = mass_tuple->addItem ("m1", m_m1);
00280 status = mass_tuple->addItem ("m2", m_m2);
00281 status = mass_tuple->addItem ("m3", m_m3);
00282 status = mass_tuple->addItem ("costheta1", m_cos1);
00283 status = mass_tuple->addItem ("costheta2", m_cos2);
00284 status = mass_tuple->addItem ("costheta3", m_cos3);
00285 status = mass_tuple->addItem ("ichannel", m_ich);
00286 } else {
00287 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
00288 return StatusCode::FAILURE;
00289 }
00290 }
00291
00292 NTuplePtr nt3(ntupleSvc(),"MYROOTFILE/massGen");
00293 if(nt3) {massgen_tuple=nt3;}
00294 else {
00295 massgen_tuple = ntupleSvc()->book ("MYROOTFILE/massGen", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
00296 if(massgen_tuple){
00297 status = massgen_tuple->addItem ("m12", _m12);
00298 status = massgen_tuple->addItem ("m13", _m13);
00299 status = massgen_tuple->addItem ("m23", _m23);
00300 status = massgen_tuple->addItem ("m1", _m1);
00301 status = massgen_tuple->addItem ("m2", _m2);
00302 status = massgen_tuple->addItem ("m3", _m3);
00303 status = massgen_tuple->addItem ("costheta1", _cos1);
00304 status = massgen_tuple->addItem ("costheta2", _cos2);
00305 status = massgen_tuple->addItem ("costheta3", _cos3);
00306 status = massgen_tuple->addItem ("ichannel", _ich);
00307 } else {
00308 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
00309 return StatusCode::FAILURE;
00310 }
00311 }
00312
00313
00314 }
00315 for(int i=0;i<500;i++){br[i]=0;}
00316 if(!(m_SB3File=="" && m_SB3HT=="")) {
00317 const char *filename, *histitle;
00318 filename=(char*)m_SB3File.c_str();
00319 histitle=(char*)m_SB3HT.c_str();
00320 SuperBody3decay.setFile(filename);
00321 SuperBody3decay.setHTitle(histitle);
00322 SuperBody3decay.init();
00323 }
00324
00325 return StatusCode::SUCCESS;
00326 }
00327
00328 StatusCode EvtDecay::execute()
00329 {
00330
00331 MsgStream log(messageService(), name());
00332
00333
00334 string key = "GEN_EVENT";
00335
00336 SmartDataPtr<McGenEventCol> McEvtColl(eventSvc(), "/Event/Gen");
00337
00338 m_numberEvent += 1;
00339 writeFlag = 0;
00340
00341 if ( McEvtColl == 0 )
00342 {
00343 HepMC::GenEvent* evt=new GenEvent();
00344 evt->set_event_number(m_numberEvent);
00345
00346
00347 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00348 int runNo=eventHeader->runNumber();
00349 int event=eventHeader->eventNumber();
00350 bool newRunFlag=0;
00351 if(runNo != 0 && runNo != m_runNo){m_runNo=runNo;newRunFlag = true;}else{newRunFlag=false;}
00352 if(m_RdMeasuredEcms&& newRunFlag){
00353 runNo = std::abs(runNo);
00354 ReadME theME(runNo);
00355 if(theME.isRunNoValid()){
00356 dbEcms=theME.getEcms();
00357 std::cout<<"Read Ecms= "<<dbEcms<<std::endl;
00358 if(dbEcms!=0){
00359 std::cout<<"INFO: Read the MeasuredEcms"<<std::endl;
00360 }
00361 else{
00362 std::cout<<"ERROR: Can not read the MeasuredEcms, try to turn off the ReadEcmsDB"<<std::endl;
00363 return StatusCode::FAILURE;
00364 }
00365 }
00366 }
00367
00368
00369 callBesEvtGen( evt );
00370 if(!(m_DecayTop=="")) evt->print(outfile);
00371
00372
00373 McGenEventCol *mcColl = new McGenEventCol;
00374 McGenEvent* mcEvent = new McGenEvent(evt);
00375 mcColl->push_back(mcEvent);
00376 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00377 if(sc != SUCCESS) return StatusCode::FAILURE;
00378
00379 } else
00380 {
00381 McGenEventCol::iterator mcItr;
00382 for( mcItr = McEvtColl->begin(); mcItr != McEvtColl->end(); mcItr++ )
00383 {
00384 HepMC::GenEvent* hEvt = (*mcItr)->getGenEvt();
00385
00386
00387 callEvtGen( hEvt );
00388
00389 if(!(m_DecayTop=="")) hEvt->print(outfile);
00390
00391
00392 if(!(m_DecayRec=="")) std::cout<<std::endl;
00393 };
00394 }
00395 if( m_NtupleFile ){ m_tuple->write();}
00396 return StatusCode::SUCCESS;
00397 }
00398
00399
00400
00401
00402
00403
00404 StatusCode EvtDecay::callEvtGen( HepMC::GenEvent* hepMCevt )
00405 {
00406 MsgStream log(messageService(), name());
00407 nchr = 0;
00408 nchr_e = 0;
00409 nchr_mu= 0;
00410 nchr_pi= 0;
00411 nchr_k = 0;
00412 nchr_p = 0;
00413 N_gamma=0;
00414 N_gammaFSR = 0;
00415 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00416
00417 for (int i=0;i<13;i++) fst[i]=0;
00418 HepMC::GenEvent::particle_const_iterator itp;
00419 AllTrk_index=0;
00420 for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp )
00421 {
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442 HepMC::GenParticle* hepMCpart=*itp;
00443 int stat = hepMCpart->status();
00444
00445
00446 if ( stat != 1 ) continue;
00447
00448 loop_to_select_eventsB:
00449 int id = hepMCpart->pdg_id();
00450 parentPDGcode=id;
00451 hepMCpart->set_status(899);
00452 EvtId eid=EvtPDL::evtIdFromStdHep(id);
00453 string parentName=EvtPDL::name(eid);
00454
00455 double en =(hepMCpart->momentum()).e();
00456 double pex=(hepMCpart->momentum()).px();
00457 double pey=(hepMCpart->momentum()).py();
00458 double pez=(hepMCpart->momentum()).pz();
00459
00460 EvtVector4R p_init(en,pex,pey,pez);
00461 parentMass=p_init.mass();
00462 EvtPDL::reSetMass(eid,parentMass);
00463
00464 EvtParticle* part=EvtParticleFactory::particleFactory(eid,p_init);
00465
00466
00467 bool parentFlag=false;
00468 if( writeFlag==0 && part->getSpinType() == EvtSpinType::VECTOR) {
00469 if(m_polarization.size()>0) {
00470 part->setPolarizedSpinDensity(m_polarization[0],m_polarization[1],m_polarization[2]);
00471 }else if(m_eBeamPolarization>0){
00472 part->setPolarizedSpinDensity(1+ m_eBeamPolarization,0,1- m_eBeamPolarization);
00473 } else{
00474 part->setVectorSpinDensity();
00475 }
00476 parentFlag=true;
00477 writeFlag++;
00478 } else {
00479 int id = hepMCpart->pdg_id();
00480 if( id == 22 ) {N_gamma ++;fst[0]++;}
00481 if( id == -22 ){N_gammaFSR ++;}
00482 if( id == 11 ) {fst[1]++;} else if(id == -11){fst[2]++;}
00483 else if(id == 13){fst[3]++;}
00484 else if(id == -13){fst[4]++;}
00485 else if(id == 211){fst[5]++;}
00486 else if(id ==-211){fst[6]++;}
00487 else if(id == 321){fst[7]++;}
00488 else if(id ==-321){fst[8]++;}
00489 else if(id ==2212){fst[9]++;}
00490 else if(id ==-2212){fst[10]++;}
00491 else if(id == 2112){fst[11]++;}
00492 else if(id ==-2112){fst[12]++;}
00493 if( fabs(id) == 11) {nchr_e++;}
00494 if( fabs(id) == 13) {nchr_mu++;}
00495 if( fabs(id) == 211) {nchr_pi++;}
00496 if( fabs(id) == 321) {nchr_k++;}
00497 if( fabs(id) == 2212) {nchr_p++;}
00498
00499
00500 if( m_NtupleFile==true ){
00501 m_Trk_index[AllTrk_index]=id;
00502 m_px_trk[AllTrk_index]=pex;
00503 m_py_trk[AllTrk_index]=pey;
00504 m_pz_trk[AllTrk_index]=pez;
00505 m_en_trk[AllTrk_index]=en ;
00506
00507 AllTrk_index++;
00508 Trk_index[AllTrk_index]=0;
00509 }
00510
00511
00512 }
00513
00514
00515
00516 EvtVector4R fdp_init;
00517 EvtId fdp_ppid;
00518 if(m_FDPparticle!=""){
00519 findPart(part);
00520 fdp_ppid = FDP_id;
00521 fdp_init = FDP_init;
00522 } else{
00523 fdp_ppid = eid;
00524 fdp_init = p_init;
00525 }
00526
00527 if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && first ){
00528 SuperBody3decay_make(fdp_ppid, fdp_init );
00529 }
00530
00531 loop_to_select_eventsA:
00532 m_Gen->generateDecay(part);
00533 if(m_truthFile!="" && m_truthPart!=""){
00534 if(EvtPDL::name(part->getId())==m_truthPart ){
00535 int ndaug = part->getNDaug();
00536 for(int id=0;id<ndaug;id++){
00537 EvtParticle* sonp=part->getDaug(id);
00538 EvtVector4R son=sonp->getP4();
00539 truth<<son.get(0)<<" "<<son.get(1)<<" "<<son.get(2)<<" "<<son.get(3)<<std::endl;
00540 }
00541 }
00542 }
00543 double ratio,rdm;
00544 if(_mDIY && parentFlag && m_ampscalflag ) _amps_max=CalAmpsMax(part);
00545 if(_mDIY && parentFlag) {
00546 rdm=EvtRandom::Flat(0.0, 1.0);
00547 double amps=CalAmpsMDIY( part );
00548 ratio=amps/_amps_max;
00549 }
00550 if(_mDIY && parentFlag && ratio<=rdm){ goto loop_to_select_eventsA;}
00551
00552
00553 bool accept;
00554 if(m_FDPparticle==""){FDP_part=part;}
00555 if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag ){ accept=SuperBody3decay_judge( FDP_part); }
00556 if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && !accept){
00557 part->deleteTree();
00558 writeFlag=0;
00559 goto loop_to_select_eventsB;
00560 }else if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && accept){
00561 countChannel(FDP_part);
00562 }
00563
00564 if((m_Psi2openCharm || m_Psi4040OpenCharm) && parentFlag ){
00565 if(getModel(part) == "OPENCHARM"){
00566 std::cout<<"OPENCHARM model need to set Psi2openCharm and Psi3SopenCharm to be false!"<<std::endl;
00567 abort();
00568 }
00569 EvtPsi3Sdecay opencharm(en, part);
00570 bool theDecay = opencharm.choseDecay();
00571 if(!theDecay ) {
00572 part->deleteTree();
00573 writeFlag=0;
00574 goto loop_to_select_eventsB;
00575 }
00576 }
00577
00578
00579
00580
00581 if(parentFlag){
00582 EvtDecayTag theTag(part);
00583 unsigned int modeTag, multiplicityTag;
00584 modeTag = theTag.getModeTag();
00585 if( getModel(part) == "OPENCHARM"||getModel(part) == "LUNDCHARM" && m_tagLundModel){
00586 multiplicityTag = decayType(part);
00587
00588
00589 } else {
00590 multiplicityTag = theTag.getMultTag() + decayType(part);
00591 }
00592
00593
00594 if(eventHeader == 0) std::cout << "event header unavailable: " << std::endl;
00595 if (eventHeader != 0) {
00596 eventHeader->setFlag1(modeTag);
00597 eventHeader->setFlag2(multiplicityTag);
00598
00599 }
00600 }
00601
00602 if(!(m_DecayRec=="")) mydecayrec=part->writeTreeRec(m_DecayRec);
00603
00604 if(m_statDecays && parentFlag ) countChannel(part);
00605
00606
00607 if ( log.level() <= MSG::DEBUG ) part->printTree() ;
00608 log << MSG::DEBUG << "Converting particles to HepMC" << endreq;
00609
00610 makeHepMC(part, hepMCevt, hepMCpart);
00611 if(part->getNDaug()!=0)
00612 hepMCpart->set_status(999);
00613
00614
00615
00616 if(part->getId()==EvtPDL::getId(m_outputp4)){
00617 int nds=part->getNDaug();
00618 for(int i=0;i<nds;i++){
00619 if(EvtPDL::name(part->getDaug(i)->getId())=="gammaFSR") continue;
00620 EvtVector4R vp1=part->getDaug(i)->getP4Lab();
00621 std::cout<<"vvpp "<<vp1.get(0)<<" "<<vp1.get(1)<<" "<<vp1.get(2)<<" "<<vp1.get(3)<<std::endl;
00622 }
00623 }
00624
00625 if(m_ReadTruth.size())ReadTruth(part,m_ReadTruth);
00626
00627 part->deleteTree();
00628 };
00629
00630 for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp )
00631 {
00632 if (!(*itp)->end_vertex() && ( (*itp)->status()==899 || (*itp)->status()==999 ))
00633 (*itp)->set_status(1);
00634 }
00635
00636 nchr = nchr_pi + nchr_k +nchr_p;
00637 if(m_Ncharge == true ) {std::cout<<"Multiplicity: TotoalCharge_e_mu_pi_K_p_gamma_gammaFSR: "
00638 <<nchr <<" "
00639 <<nchr_e <<" "
00640 << nchr_mu <<" "
00641 << nchr_pi <<" "
00642 << nchr_k <<" "
00643 <<nchr_p <<" "
00644 <<N_gamma <<" "
00645 <<N_gammaFSR<<endl;}
00646 if(m_Ncharge == true ){std::cout<<"FinalStates: "
00647 <<fst[0]<<" "
00648 <<fst[1]<<" "
00649 <<fst[2]<<" "
00650 <<fst[3]<<" "
00651 <<fst[4]<<" "
00652 <<fst[5]<<" "
00653 <<fst[6]<<" "
00654 <<fst[7]<<" "
00655 <<fst[8]<<" "
00656 <<fst[9]<<" "
00657 <<fst[10]<<" "
00658 <<fst[11]<<" "
00659 <<fst[12]<<std::endl;}
00660 if( m_NtupleFile ){
00661
00662 m_nchr = nchr;
00663 m_nchr_e = nchr_e;
00664 m_nchr_mu = nchr_mu;
00665 m_nchr_pi = nchr_pi;
00666 m_nchr_k = nchr_k;
00667 m_nchr_p = nchr_p;
00668 m_gamma=N_gamma;
00669 m_gammaFSR=N_gammaFSR;
00670 TotNumTrk = AllTrk_index;
00671 }
00672
00673
00674 return StatusCode::SUCCESS;
00675 }
00676
00677
00678
00679
00680 StatusCode EvtDecay::callBesEvtGen( HepMC::GenEvent* hepMCevt )
00681 {
00682 MsgStream log(messageService(), name());
00683
00684 nchr = -1;
00685 nchr_e = 0;
00686 nchr_mu= 0;
00687 nchr_pi= 0;
00688 nchr_k = 0;
00689 nchr_p = 0;
00690 N_gamma=0;
00691 N_gammaFSR = 0;
00692 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00693 for (int i=0;i<13;i++) fst[i]=0;
00694 HepMC::GenEvent::particle_const_iterator itp;
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715 loop_to_select_eventsB:
00716 EvtId ppid=EvtPDL::getId(m_ParentPart);
00717 if(m_RdMeasuredEcms) EvtPDL::reSetMass(ppid,dbEcms);
00718 double ppmass=EvtPDL::getMass(ppid);
00719
00720 parentMass=ppmass;
00721 int pid=EvtPDL::getStdHep(ppid);
00722 parentPDGcode=pid;
00723
00724 EvtVector4R p_init(ppmass,0.0,0.0,0.0);
00725
00726 EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init);
00727 HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,ppmass),pid,3);
00728
00729 bool parentFlag = false;
00730
00731 if(writeFlag ==0 && part->getSpinType() == EvtSpinType::VECTOR) {
00732 if(m_polarization.size()>0) {
00733 part->setPolarizedSpinDensity(m_polarization[0],m_polarization[1],m_polarization[2]);
00734 } else if(m_eBeamPolarization>0){
00735 part->setPolarizedSpinDensity(1+ m_eBeamPolarization,0,1- m_eBeamPolarization);
00736 } else{
00737 part->setVectorSpinDensity();
00738 }
00739 parentFlag = true;
00740 writeFlag++;
00741 }
00742
00743
00744 EvtVector4R fdp_init;
00745 EvtId fdp_ppid;
00746 if(m_FDPparticle!=""){
00747 findPart(part);
00748 fdp_ppid = FDP_id;
00749 fdp_init = FDP_init;
00750 } else{
00751 fdp_ppid = ppid;
00752 fdp_init = p_init;
00753 }
00754
00755 if(!(m_SB3File=="" || m_SB3HT=="") && first ){ SuperBody3decay_make( fdp_ppid, fdp_init);}
00756
00757
00758 loop_to_select_events:
00759 m_Gen->generateDecay(part); if(m_numberEvent<=1){ std::cout<<"after m_Gen->decay "<<std::endl; part->printTree();}
00760
00761 double ratio,rdm;
00762 if(_mDIY && m_ampscalflag ) _amps_max=CalAmpsMax(part);
00763
00764 if(_mDIY ) {
00765 for(int myloop=0;myloop<100;myloop++){
00766 m_Gen->generateDecay(part);
00767 double amps=CalAmpsMDIY( part);
00768 ratio=amps/_amps_max;
00769 rdm=EvtRandom::Flat(0.0, 1.0);
00770 if( !isNumber(amps) || !isNumber(_amps_max) || ratio<=0 ) {
00771 part->deleteTree();
00772 part=EvtParticleFactory::particleFactory(ppid,p_init);
00773 continue;
00774 }else if( rdm<=ratio) {
00775 break;
00776 }
00777 }
00778 }
00779
00780 if(m_truthFile!="" && m_truthPart!=""){
00781 if(EvtPDL::name(part->getId())==m_truthPart ){
00782 int ndaug = part->getNDaug();
00783 for(int id=0;id<ndaug;id++){
00784 EvtParticle* sonp=part->getDaug(id);
00785 EvtVector4R son=sonp->getP4();
00786 truth<<son.get(0)<<" "<<son.get(1)<<" "<<son.get(2)<<" "<<son.get(3)<<std::endl;
00787 }
00788 }
00789 }
00790
00791 bool accept;
00792 if(m_FDPparticle==""){FDP_part=part;}
00793 if(!(m_SB3File=="" || m_SB3HT=="")){
00794 accept=SuperBody3decay_judge( FDP_part);
00795 }
00796 if(!(m_SB3File=="" || m_SB3HT=="") && !accept) {
00797 delete hepMCpart;
00798 part->deleteTree();
00799 goto loop_to_select_eventsB;
00800 }else if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && accept){
00801 countChannel(FDP_part);
00802 }
00803
00804 int Nchannel=part->getChannel();
00805
00806 if(m_statDecays && parentFlag) countChannel(part);
00807
00808
00809
00810
00811 if(parentFlag){
00812 EvtDecayTag theTag(part);
00813 unsigned int modeTag, multiplicityTag;
00814 modeTag = theTag.getModeTag();
00815 if( getModel(part) == "OPENCHARM"|| getModel(part) == "LUNDCHARM" && m_tagLundModel ){
00816 multiplicityTag = decayType(part);
00817
00818
00819 } else {
00820 multiplicityTag = theTag.getMultTag() + decayType(part);
00821 }
00822 if(eventHeader == 0) std::cout << "event header unavailable: " << std::endl;
00823 if (eventHeader != 0) {
00824 eventHeader->setFlag1(modeTag);
00825 eventHeader->setFlag2(multiplicityTag);
00826 }
00827 if(m_NtupleFile)m_flag1 = modeTag;
00828
00829 }
00830
00831 if(!(m_DecayRec=="")) {mydecayrec=part->writeTreeRec(m_DecayRec);std::cout<<std::endl;};
00832
00833
00834 if ( log.level() <= MSG::DEBUG ) part->printTree() ;
00835 log << MSG::DEBUG << "Converting particles to HepMC" << endreq;
00836
00837 makeHepMC(part, hepMCevt, hepMCpart);
00838
00839 if(part->getNDaug()!=0) hepMCpart->set_status(999);
00840
00841
00842
00843 AllTrk_index=0;
00844 for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp )
00845 {
00846 if (!(*itp)->end_vertex() && ( (*itp)->status()==899 || (*itp)->status()==999 ))
00847 (*itp)->set_status(1);
00848
00849 HepMC::GenParticle* hepMCpart=*itp;
00850 int stat = hepMCpart->status();
00851 int id = hepMCpart->pdg_id();
00852 if(abs(id)==411 ||abs(id)==413 ) { (*itp)->set_status(999);stat=999;}
00853
00854 if ( stat != 1 ) continue;
00855 if( id == 22 ) {N_gamma ++;fst[0]++;}
00856 if( id == -22 ){N_gammaFSR ++;}
00857 if(id == 11 ) {fst[1]++;}
00858 else if(id == -11) {fst[2]++;}
00859 else if(id == 13 ) {fst[3]++;}
00860 else if(id ==-13) {fst[4]++;}
00861 else if(id==211) {fst[5]++;}
00862 else if(id==-211) {fst[6]++;}
00863 else if(id== 321) {fst[7]++;}
00864 else if(id ==-321) {fst[8]++;}
00865 else if(id ==2212) {fst[9]++;}
00866 else if(id==-2212) {fst[10]++;}
00867 else if(id == 2112){fst[11]++;}
00868 else if(id==-2112) {fst[12]++;}
00869 if( fabs(id) == 11) {nchr_e++;}
00870 if( fabs(id) == 13) {nchr_mu++;}
00871 if( fabs(id) == 211) {nchr_pi++;}
00872 if( fabs(id) == 321) {nchr_k++;}
00873 if( fabs(id) == 2212) {nchr_p++;}
00874
00875
00876 double en =(hepMCpart->momentum()).e();
00877 double pex=(hepMCpart->momentum()).px();
00878 double pey=(hepMCpart->momentum()).py();
00879 double pez=(hepMCpart->momentum()).pz();
00880
00881 if( m_NtupleFile==true && id!=0){
00882 m_Trk_index[AllTrk_index]=id;
00883 m_px_trk[AllTrk_index]=pex;
00884 m_py_trk[AllTrk_index]=pey;
00885 m_pz_trk[AllTrk_index]=pez;
00886 m_en_trk[AllTrk_index]=en ;
00887
00888
00889
00890
00891
00892
00893 AllTrk_index++;
00894 }
00895
00896 }
00897
00898 itp=hepMCevt->particles_begin();
00899 (*itp)->set_status(2);
00900
00901
00902 nchr = nchr_pi + nchr_k +nchr_p;
00903 if(m_Ncharge == true ) {std::cout<<"Multiplicity: TotoalCharge_e_mu_pi_K_p_gamma_gammaFSR: "
00904 <<nchr<<" "
00905 <<nchr_e<<" "
00906 << nchr_mu <<" "
00907 << nchr_pi <<" "
00908 << nchr_k <<" "
00909 <<nchr_p<<" "
00910 <<N_gamma<<" "
00911 <<N_gammaFSR<<endl;}
00912 if(m_Ncharge == true ){std::cout<<"FinalStates: "
00913 <<fst[0]<<" "
00914 <<fst[1]<<" "
00915 <<fst[2]<<" "
00916 <<fst[3]<<" "
00917 <<fst[4]<<" "
00918 <<fst[5]<<" "
00919 <<fst[6]<<" "
00920 <<fst[7]<<" "
00921 <<fst[8]<<" "
00922 <<fst[9]<<" "
00923 <<fst[10]<<" "
00924 <<fst[11]<<" "
00925 <<fst[12]<<std::endl;}
00926
00927
00928 if( m_NtupleFile ){
00929 m_nchr = nchr;
00930 m_nchr_e = nchr_e;
00931 m_nchr_mu = nchr_mu;
00932 m_nchr_pi = nchr_pi;
00933 m_nchr_k = nchr_k;
00934 m_nchr_p = nchr_p;
00935 m_gamma=N_gamma;
00936 m_gammaFSR=N_gammaFSR;
00937 TotNumTrk = AllTrk_index;
00938 }
00939
00940
00941 if(m_ReadTruth.size())ReadTruth(part,m_ReadTruth);
00942
00943
00944
00945 part->deleteTree();
00946 return StatusCode::SUCCESS;
00947 }
00948
00949
00950 StatusCode EvtDecay::finalize()
00951 {
00952
00953 if(EvtCalHelAmp::nevt>0){
00954 double H2=EvtCalHelAmp::_H2;
00955 double H1=EvtCalHelAmp::_H1;
00956 double H0=EvtCalHelAmp::_H0;
00957 double H2err=EvtCalHelAmp::_H2err;
00958 double H1err=EvtCalHelAmp::_H1err;
00959 double H0err=EvtCalHelAmp::_H0err;
00960 int nd = EvtCalHelAmp::nevt;
00961 std::cout<<"Calculated helicity amplitude sqaured are:"<<std::endl;
00962 std::cout<<"|H_{2}|^2=|H_{-2}|^2= "<<H2/nd<<" +/- "<<sqrt(H2err/nd)<<endl;
00963 std::cout<<"|H_{1}|^2=|H_{-1}|^2= "<<H1/nd<<" +/- "<<sqrt(H1err/nd)<<endl;
00964
00965 }
00966 MsgStream log(messageService(), name());
00967 truth.close();
00968 log << MSG::INFO << "EvtDecay finalized" << endreq;
00969 fstream Forfile;
00970 Forfile.open("fort.16",ios::in);
00971 char delcmd[300];
00972 strcpy(delcmd,"rm -f ");
00973 strcat(delcmd,"fort.16");
00974
00975
00976 fstream Forfile2;
00977 Forfile2.open("fort.22",ios::in);
00978 strcpy(delcmd,"rm -f ");
00979 strcat(delcmd,"fort.22");
00980
00981
00982
00983
00984
00985
00986
00987
00988 int totalEvt=0;
00989 if(!(m_SB3File=="" || m_SB3HT=="")){
00990 for(int i=0;i<500;i++){totalEvt=totalEvt+br[i];}
00991 for(int i=0;i<500;i++){
00992 double bi=1.0*br[i]/totalEvt;
00993 std::cout<<"Branching fraction of channel "<<i<<" is: "<<bi<<std::endl;
00994 if(br[i]==0) break;
00995 }
00996 }
00997
00998 if(m_statDecays){
00999 int totalEvt=0;
01000 for(int i=0;i<=totalChannels;i++){totalEvt=totalEvt+br[i];}
01001 std::cout<<"==================Statistical first chain decay ==============================="<<std::endl;
01002 std::cout<<"i-th channel Events Generated Branching fraction generated "<<std::endl;
01003 for(int i=0;i<=totalChannels;i++){
01004 std::cout<<"| "<<i<<" "<<br[i]<<" "<<br[i]*1.0/totalEvt<<std::endl;
01005
01006 }
01007 std::cout<<"--------------------------------------------------------------------------------"<<std::endl;
01008 std::cout<<" sum: "<<totalEvt<<" "<<"1"<<std::endl;
01009 std::cout<<"================== End of statistical first chain decay ========================"<<std::endl;
01010 }
01011
01012 return StatusCode::SUCCESS;
01013 }
01014
01015
01016 StatusCode EvtDecay::makeHepMC(EvtParticle* part, HepMC::GenEvent* hEvt,
01017 HepMC::GenParticle* hPart)
01018 {
01019 MsgStream log(messageService(), name());
01020
01021 if(part->getNDaug()!=0)
01022 {
01023 double ct=(part->getDaug(0)->get4Pos()).get(0);
01024 double x=(part->getDaug(0)->get4Pos()).get(1);
01025 double y=(part->getDaug(0)->get4Pos()).get(2);
01026 double z=(part->getDaug(0)->get4Pos()).get(3);
01027
01028 HepMC::GenVertex* end_vtx = new HepMC::GenVertex(HepLorentzVector(x,y,z,ct));
01029 hEvt->add_vertex( end_vtx );
01030 end_vtx->add_particle_in(hPart);
01031
01032 int ndaug=part->getNDaug();
01033
01034 for(int it=0;it<ndaug;it++)
01035 {
01036
01037 double e=(part->getDaug(it)->getP4Lab()).get(0);
01038 double px=(part->getDaug(it)->getP4Lab()).get(1);
01039 double py=(part->getDaug(it)->getP4Lab()).get(2);
01040 double pz=(part->getDaug(it)->getP4Lab()).get(3);
01041 int id=EvtPDL::getStdHep(part->getDaug(it)->getId());
01042 int status=1;
01043
01044 if(part->getDaug(it)->getNDaug()!=0) status=999;
01045
01046 HepMC::GenParticle* prod_part = new HepMC::GenParticle(HepLorentzVector(px,py,pz,e),id,status);
01047 end_vtx->add_particle_out(prod_part);
01048 makeHepMC(part->getDaug(it),hEvt,prod_part);
01049
01050
01051
01052 if( log.level()<MSG::INFO )
01053 prod_part->print();
01054 }
01055 }
01056
01057 return StatusCode::SUCCESS;
01058 }
01059
01060
01061 void EvtDecay::MeVToGeV (HepMC::GenEvent* evt)
01062 {
01063 for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) {
01064
01065
01066
01067 HepMC::FourVector tmpfv((*p)->momentum().x()/1000.,
01068 (*p)->momentum().y()/1000.,
01069 (*p)->momentum().z()/1000.,
01070 (*p)->momentum().t()/1000.
01071 );
01072
01073 (*p)->set_momentum( tmpfv );
01074 }
01075 }
01076
01077 void EvtDecay::GeVToMeV (HepMC::GenEvent* evt)
01078 {
01079 for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) {
01080
01081
01082 HepMC::FourVector tmpfv((*p)->momentum().x()*1000.,
01083 (*p)->momentum().y()*1000.,
01084 (*p)->momentum().z()*1000.,
01085 (*p)->momentum().t()*1000.
01086 );
01087
01088 (*p)->set_momentum( tmpfv );
01089 }
01090 }
01091
01092
01093 double EvtDecay::CalAmpsMax( EvtParticle* part )
01094 {
01095 double amps_max=0;
01096 for(int i=0;i<100000;i++){
01097 m_Gen->generateDecay(part);
01098 double amps=CalAmpsMDIY(part);
01099 if(amps > amps_max) amps_max=amps;
01100 }
01101 amps_max=amps_max*1.05;
01102 m_ampscalflag = false;
01103 return amps_max;
01104 }
01105
01106
01107
01108 EvtBesRandom::EvtBesRandom(HepRandomEngine* engine)
01109 {
01110 m_engine = engine;
01111 m_engine->showStatus();
01112 }
01113
01114 EvtBesRandom::~EvtBesRandom()
01115 {}
01116
01117 double EvtBesRandom::random()
01118 {
01119 return RandFlat::shoot(m_engine);
01120 }
01121
01122
01123 void EvtDecay::SuperBody3decay_make(EvtId ppid, EvtVector4R p_init ){
01124 double xmass2,ymass2;
01125 bool accept;
01126 EvtVector4R pd1,pd2,pd3;
01127
01128
01129 FinalState_make( ppid, p_init);
01130
01131
01132 for(int i=0;i<100000;i++){
01133 regen:
01134 EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init);
01135 HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,parentMass),parentPDGcode,3);
01136
01137 if(part->getSpinType() == EvtSpinType::VECTOR) {
01138 part->setVectorSpinDensity();
01139 }
01140 m_Gen->generateDecay(part);
01141
01142
01143 FinalState_sort(part);
01144 pd1=son0;
01145 pd2=son1;
01146 pd3=son2;
01147
01148
01149
01150 xmass2=(pd1+pd2).mass2();
01151 ymass2=(pd1+pd3).mass2();
01152 double xlow=SuperBody3decay.getXlow();
01153 double xup=SuperBody3decay.getXup();
01154 double ylow=SuperBody3decay.getYlow();
01155 double yup=SuperBody3decay.getYup();
01156 if(xmass2<xlow || xmass2>xup || ymass2<ylow || ymass2>yup) { part->deleteTree();goto regen;}
01157 SuperBody3decay.HFill(xmass2,ymass2);
01158
01159 if( m_NtupleFile ){
01160 m_ich=part->getChannel();
01161 m_m1=pd1.mass();
01162 m_m2=pd2.mass();
01163 m_m3=pd3.mass();
01164 m_m12=(pd1+pd2).mass();
01165 m_m23=(pd2+pd3).mass();
01166 m_m13=(pd1+pd3).mass();
01167 m_cos1=pd1.get(3)/pd1.d3mag();
01168 m_cos2=pd2.get(3)/pd2.d3mag();
01169 m_cos3=pd3.get(3)/pd3.d3mag();
01170 mass_tuple->write();
01171 }
01172 double m1=pd1.mass(); double m2=pd2.mass();double m3=pd3.mass();
01173 double mj=(pd2+pd3).mass();
01174
01175
01176 }
01177
01178 SuperBody3decay.HReweight();
01179 SuperBody3decay.setZmax();
01180 first=false;
01181 }
01182
01183 bool EvtDecay::SuperBody3decay_judge(EvtParticle* part){
01184 double xmass2,ymass2;
01185 bool accept;
01186 EvtVector4R pd1,pd2,pd3;
01187
01188
01189 FinalState_sort( part);
01190
01191 pd1=son0;
01192 pd2=son1;
01193 pd3=son2;
01194
01195
01196 xmass2=(pd1+pd2).mass2();
01197 ymass2=(pd1+pd3).mass2();
01198 accept=SuperBody3decay.AR(xmass2,ymass2);
01199
01200
01201 if(accept && m_NtupleFile) {
01202 _ich=part->getChannel();
01203 _m1=pd1.mass();
01204 _m2=pd2.mass();
01205 _m3=pd3.mass();
01206 _m12=(pd1+pd2).mass();
01207 _m23=(pd2+pd3).mass();
01208 _m13=(pd1+pd3).mass();
01209 _cos1=pd1.get(3)/pd1.d3mag();
01210 _cos2=pd2.get(3)/pd2.d3mag();
01211 _cos3=pd3.get(3)/pd3.d3mag();
01212 massgen_tuple->write();
01213 }
01214
01215
01216 return accept;
01217 }
01218
01219
01220 void EvtDecay:: FinalState_make(EvtId ppid, EvtVector4R p_init){
01221
01222 multi=1;
01223 identical_flag=true;
01224
01225 for(int i=1;i<10000;i++){
01226 EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init);
01227 HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,parentMass),parentPDGcode,3);
01228
01229 m_Gen->generateDecay(part);
01230 int nson=part->getNDaug();
01231
01232 if(nson == 2) {continue;} else
01233 if(nson ==3){
01234 EvtId xid0,xid1,xid2;
01235 xid0=part->getDaug(0)->getId();
01236 xid1=part->getDaug(1)->getId();
01237 xid2=part->getDaug(2)->getId();
01238
01239
01240 pdg0=EvtPDL::getStdHep(xid0);
01241 pdg1=EvtPDL::getStdHep(xid1);
01242 pdg2=EvtPDL::getStdHep(xid2);
01243
01244 if(pdg0==pdg1 && pdg1==pdg2) {multi=3;}
01245 else if(pdg0==pdg1 ){multi=2;}
01246 else if(pdg0==pdg2 ){multi=2;pdg=pdg1; pdg1=pdg2; pdg2=pdg;}
01247 else if(pdg1==pdg2) {multi=2;pdg=pdg0; pdg0=pdg1; pdg1=pdg2; pdg2=pdg;}
01248 else {multi=1;}
01249 break;
01250 } else{
01251 std::cout<<"No direct three body decay channel found in decay card, I stop run"<<std::endl;
01252 ::abort();
01253 }
01254 }
01255
01256
01257 }
01258
01259 void EvtDecay::FinalState_sort(EvtParticle* part){
01260
01261 EvtVector4R pd0,pd1,pd2;
01262 EvtId xid0,xid1,xid2;
01263 int id0,id1,id2;
01264
01265 int nson=part->getNDaug();
01266 if(nson==2){
01267 pd0=part->getDaug(0)->getP4Lab();
01268 pd1=part->getDaug(1)->getDaug(0)->getP4Lab();
01269 pd2=part->getDaug(1)->getDaug(1)->getP4Lab();
01270
01271 xid0=part->getDaug(0)->getId();
01272 xid1=part->getDaug(1)->getDaug(0)->getId();
01273 xid2=part->getDaug(1)->getDaug(1)->getId();
01274
01275 } else if(nson==3){
01276 pd0=part->getDaug(0)->getP4Lab();
01277 pd1=part->getDaug(1)->getP4Lab();
01278 pd2=part->getDaug(2)->getP4Lab();
01279
01280 xid0=part->getDaug(0)->getId();
01281 xid1=part->getDaug(1)->getId();
01282 xid2=part->getDaug(2)->getId();
01283 }
01284
01285 id0=EvtPDL::getStdHep(xid0);
01286 id1=EvtPDL::getStdHep(xid1);
01287 id2=EvtPDL::getStdHep(xid2);
01288
01289
01290
01291 if(multi==1){
01292 assign_momentum(id0,pd0);
01293 assign_momentum(id1,pd1);
01294 assign_momentum(id2,pd2);
01295 } else if(multi==2){
01296 assign_momentum2(id0,pd0);
01297 assign_momentum2(id1,pd1);
01298 assign_momentum2(id2,pd2);
01299 if(son0.get(0)>son1.get(0)){son=son0;son0=son1;son1=son;}
01300 } else if(multi==3){
01301 double En0=pd0.get(0);
01302 double En1=pd1.get(0);
01303 double En2=pd2.get(0);
01304 if(En0 < En1 && En1 < En2) {son0=pd0;son1=pd1;son2=pd2;}
01305 if(En0 < En2 && En2 < En1) {son0=pd0;son1=pd2;son2=pd1;}
01306 if(En1 < En0 && En0 < En2) {son0=pd1;son1=pd0;son2=pd2;}
01307 if(En1 < En2 && En2 < En0) {son0=pd1;son1=pd2;son2=pd0;}
01308 if(En2 < En0 && En0 < En1) {son0=pd2;son1=pd0;son2=pd1;}
01309 if(En2 < En1 && En1 < En0) {son0=pd2;son1=pd1;son2=pd0;}
01310 }
01311
01312 }
01313
01314
01315 void EvtDecay::assign_momentum(int id0, EvtVector4R pd0){
01316 if(id0==pdg0){son0=pd0;}
01317 else if(id0==pdg1){son1=pd0;}
01318 else if(id0==pdg2){son2=pd0;}
01319 else {std::cout<<"PID= "<<id0<<" not machted, I stop"<<std::endl; ::abort();}
01320 }
01321
01322 void EvtDecay::assign_momentum2(int id0, EvtVector4R pd0){
01323 if(id0==pdg0 && identical_flag){son0=pd0;identical_flag=false;}
01324 else if(id0==pdg1){son1=pd0;identical_flag=true;}
01325 else if(id0==pdg2){son2=pd0;}
01326 else {std::cout<<"PID not machted, I stop"<<std::endl; ::abort();}
01327 }
01328
01329 void EvtDecay::findPart(EvtParticle* part){
01330 FDP_id=EvtPDL::getId(m_FDPparticle);
01331 int FDP_pdg=EvtPDL::getStdHep(FDP_id);
01332
01333 m_Gen->generateDecay(part);
01334 int dn=part->getNDaug();
01335
01336 if(dn!=0){
01337 for(int i=0;i<dn;i++){
01338 EvtParticle* part1=part->getDaug(i);
01339 EvtId sonid=part1->getId();
01340 int son_pdg = EvtPDL::getStdHep(sonid);
01341 if(son_pdg == FDP_pdg){
01342 FDP_init=part1->getP4Lab();
01343 FDP_part=part1;
01344 return;
01345 } else{
01346 findPart(part1);
01347 }
01348 }
01349 }
01350
01351 }
01352
01353 void EvtDecay::countChannel(EvtParticle* part){
01354 int ich=part->getChannel();
01355
01356
01357
01358
01359
01360 br[ich]++;
01361 if(ich>totalChannels) totalChannels = ich;
01362
01363 }
01364
01365 bool EvtDecay::isCharmonium(EvtId xid){
01366 EvtId psi4415 = EvtPDL::getId(std::string("psi(4415)"));
01367 EvtId psi4160 = EvtPDL::getId(std::string("psi(4160)"));
01368 EvtId psi4040 = EvtPDL::getId(std::string("psi(4040)"));
01369 EvtId psi3770 = EvtPDL::getId(std::string("psi(3770)"));
01370 EvtId psiprim = EvtPDL::getId(std::string("psi(2S)"));
01371 EvtId jpsi = EvtPDL::getId(std::string("J/psi"));
01372 EvtId etac = EvtPDL::getId(std::string("eta_c"));
01373 EvtId etac2s = EvtPDL::getId(std::string("eta_c(2S)"));
01374 EvtId hc = EvtPDL::getId(std::string("h_c"));
01375 EvtId chic0 = EvtPDL::getId(std::string("chi_c0"));
01376 EvtId chic1 = EvtPDL::getId(std::string("chi_c1"));
01377 EvtId chic2 = EvtPDL::getId(std::string("chi_c2"));
01378 EvtId chic0p = EvtPDL::getId(std::string("chi_c0p"));
01379 EvtId chic1p = EvtPDL::getId(std::string("chi_c1p"));
01380 EvtId chic2p = EvtPDL::getId(std::string("chi_c1p"));
01381 std::vector<EvtId> Vid; Vid.clear();
01382 Vid.push_back(psi4415);
01383 Vid.push_back(psi4160);
01384 Vid.push_back(psi4040);
01385 Vid.push_back(psi3770);
01386 Vid.push_back(psiprim);
01387 Vid.push_back(jpsi);
01388 Vid.push_back(etac);
01389 Vid.push_back(etac2s);
01390 Vid.push_back(hc);
01391 Vid.push_back(chic0);
01392 Vid.push_back(chic1);
01393 Vid.push_back(chic2);
01394 Vid.push_back(chic0p);
01395 Vid.push_back(chic1p);
01396 Vid.push_back(chic2p);
01397
01398 bool flag=true;
01399 for(int i=0;i<Vid.size();i++){ if(xid == Vid[i]) return flag;}
01400 return false;
01401 }
01402
01403
01404 bool EvtDecay::isCharm(EvtId xid){
01405 EvtId d0 = EvtPDL::getId(std::string("D0"));
01406 EvtId d0bar = EvtPDL::getId(std::string("anti-D0"));
01407 EvtId dp = EvtPDL::getId(std::string("D+"));
01408 EvtId dm = EvtPDL::getId(std::string("D-"));
01409 EvtId d0h = EvtPDL::getId(std::string("D0H"));
01410 EvtId d0l = EvtPDL::getId(std::string("D0L"));
01411 EvtId dstp = EvtPDL::getId(std::string("D*+"));
01412 EvtId dstm = EvtPDL::getId(std::string("D*-"));
01413 EvtId ds0 = EvtPDL::getId(std::string("D*0"));
01414 EvtId ds0bar = EvtPDL::getId(std::string("anti-D*0"));
01415 EvtId dsp = EvtPDL::getId(std::string("D_s+"));
01416 EvtId dsm = EvtPDL::getId(std::string("D_s-"));
01417 EvtId dsstp = EvtPDL::getId(std::string("D_s*+"));
01418 EvtId dsstm = EvtPDL::getId(std::string("D_s*-"));
01419 EvtId ds0stp = EvtPDL::getId(std::string("D_s0*+"));
01420 EvtId ds0stm = EvtPDL::getId(std::string("D_s0*-"));
01421
01422 std::vector<EvtId> Vid; Vid.clear();
01423 Vid.push_back(d0);
01424 Vid.push_back(d0bar);
01425 Vid.push_back(dp);
01426 Vid.push_back(dm);
01427 Vid.push_back(d0h);
01428 Vid.push_back(d0l);
01429 Vid.push_back(dstp);
01430 Vid.push_back(dstm);
01431 Vid.push_back(ds0);
01432 Vid.push_back(ds0bar );
01433 Vid.push_back(dsp );
01434 Vid.push_back(dsm );
01435 Vid.push_back(dsstp );
01436 Vid.push_back(dsstm );
01437 Vid.push_back(ds0stp );
01438 Vid.push_back(ds0stm );
01439
01440 bool flag=true;
01441 for(int i=0;i<Vid.size();i++){ if(xid == Vid[i]) return flag;}
01442 return false;
01443 }
01444
01445 bool EvtDecay::isRadecay(EvtParticle *par){
01446 int ndg = par->getNDaug();
01447 for(int i=0;i<ndg;i++){
01448 EvtId xid = par->getDaug(i)->getId();
01449 if(EvtPDL::getStdHep(xid) == 22){return true;}
01450 }
01451 return false;
01452 }
01453
01454 int EvtDecay::decayType(EvtParticle *par){
01455
01456
01457
01458
01459
01460
01461
01462 int Nchannel=par->getChannel();
01463 int ndg = par->getNDaug();
01464 int nfsr=0;
01465
01466 std::string model = getModel(par,Nchannel);
01467 if( model == "OPENCHARM" || model == "LUNDCHARM" && m_tagLundModel){
01468 int ldm = par->getGeneratorFlag();
01469
01470 return ldm;
01471
01472 }
01473
01474 for(int i=0;i<ndg;i++){
01475 EvtId xid =par->getDaug(i)->getId();
01476 if(EvtPDL::getStdHep(xid) == -22 ){nfsr++;}
01477 }
01478
01479 if( isRadecay(par) ){
01480 for(int i=0;i<ndg;i++){
01481 EvtId xid = par->getDaug(i)->getId();
01482 if( isCharmonium(xid) ) return 1;
01483 }
01484 return 0;
01485 } else{
01486 if(ndg-nfsr ==2 ){
01487 int ndg = par->getNDaug();
01488 EvtId xd1 = par->getDaug(0)->getId();
01489 EvtId xd2 = par->getDaug(1)->getId();
01490 if( isCharm(xd1) && isCharm(xd2) ) {return 2;} else
01491 if( !isCharmonium(xd1) && !isCharmonium(xd2) && !isCharm(xd1) && !isCharm(xd2) ){
01492 return 3;
01493 }
01494 } else{
01495 bool flag = false;
01496 for(int i=0;i<ndg;i++){
01497 EvtId xid = par->getDaug(i)->getId();
01498 if( isCharmonium(xid) ) {flag = true; break;}
01499 if( isCharm(xid) ) {flag = true; break;}
01500 }
01501 if( !flag ){return 3;} else {return 4;}
01502 }
01503 }
01504
01505 }
01506
01507
01508 std::string EvtDecay::getModel(EvtParticle *par, int mode){
01509 EvtDecayBase* thedecay = EvtDecayTable::gettheDecay(par->getId(),mode);
01510 return thedecay->getModelName();
01511 }
01512
01513 std::string EvtDecay::getModel(EvtParticle *par){
01514 int mode = par->getChannel();
01515 EvtDecayBase* thedecay = EvtDecayTable::gettheDecay(par->getId(),mode);
01516 return thedecay->getModelName();
01517 }
01518
01519
01520 void EvtDecay::ReadTruth(EvtParticle* part,std::vector<std::vector<string> > mylist){
01521 if(mylist.size()<2) {std::cout<<" Empty list"<<std::endl;abort();}
01522 EvtVector4R vp1;
01523 for(int k=0;k<mylist[0].size();k++){
01524 if(part->getId()==EvtPDL::getId(mylist[0][k])){
01525 if(part->getDaug(1)->getId()==EvtPDL::getId("vhdr")){part=part->getDaug(1);continue;}
01526 for(int i=1;i<mylist.size();i++){
01527 EvtParticle* mypart=part;
01528 for(int j=0;j<mylist[i].size();j++){
01529 mypart=mypart->getDaug(atoi(mylist[i][j].c_str()));
01530
01531 }
01532
01533 std::cout<<EvtPDL::name(mypart->getId())<<std::endl;
01534 vp1=mypart->getP4Lab();
01535 if( !isNumber(vp1.get(0)) ) {part->printTree(); abort();}
01536 std::cout<<"truth "<<vp1.get(0)<<" "<<vp1.get(1)<<" "<<vp1.get(2)<<" "<<vp1.get(3)<<std::endl;
01537 }
01538 }
01539 }
01540 }
01541
01542 int EvtDecay::isNumber(double d){return (d==d);}
01543
01544 #include "../user/Lenu.inc"