EvtDecay Class Reference

#include <EvtDecay.h>

List of all members.

Public Member Functions
	EvtDecay (const string &name, ISvcLocator *pSvcLocator)
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	initialize ()
Public Attributes
DataInfoSvc *	dataInfoSvc
IDataInfoSvc *	tmpInfoSvc
Private Member Functions
void	assign_momentum (int pdx, EvtVector4R pv4)
void	assign_momentum2 (int pdx, EvtVector4R pv4)
double	CalAmpsMax (EvtParticle *part)
double	CalAmpsMDIY (EvtParticle *part)
StatusCode	callBesEvtGen (HepMC::GenEvent *hepMCevt)
StatusCode	callEvtGen (HepMC::GenEvent *hepMCevt)
void	countChannel (EvtParticle *part)
int	decayType (EvtParticle *par)
void	FinalState_make (EvtId ppid, EvtVector4R p_init)
void	FinalState_sort (EvtParticle *part)
void	findPart (EvtParticle *part)
std::string	getModel (EvtParticle *par)
std::string	getModel (EvtParticle *par, int mode)
void	GeVToMeV (HepMC::GenEvent *hepMCevt)
bool	isCharm (EvtId xid)
bool	isCharmonium (EvtId xid)
bool	isRadecay (EvtParticle *par)
StatusCode	makeHepMC (EvtParticle *, HepMC::GenEvent *, HepMC::GenParticle *)
void	MeVToGeV (HepMC::GenEvent *hepMCevt)
bool	SuperBody3decay_judge (EvtParticle *part)
void	SuperBody3decay_make (EvtId ppid, EvtVector4R p_init)
Private Attributes
NTuple::Item< double >	_cos1
NTuple::Item< double >	_cos2
NTuple::Item< double >	_cos3
NTuple::Item< long >	_ich
NTuple::Item< double >	_m1
NTuple::Item< double >	_m12
NTuple::Item< double >	_m13
NTuple::Item< double >	_m2
NTuple::Item< double >	_m23
NTuple::Item< double >	_m3
bool	_mDIY
int	AllTrk_index
int	br [500]
double	en_trk [50]
EvtId	FDP_id
EvtVector4R	FDP_init
EvtParticle *	FDP_part
bool	first
bool	identical_flag
bool	m_ampscalflag
NTuple::Item< double >	m_cos1
NTuple::Item< double >	m_cos2
NTuple::Item< double >	m_cos3
string	m_DecayDec
string	m_DecayRec
string	m_DecayTop
NTuple::Array< double >	m_en_trk
string	m_FDPparticle
NTuple::Array< long >	m_fst
NTuple::Item< long >	m_gamma
NTuple::Item< long >	m_gammaFSR
EvtGen *	m_Gen
NTuple::Item< long >	m_ich
std::vector< int >	m_InSeeds
EvtId	m_KKMCRes
NTuple::Item< double >	m_m1
NTuple::Item< double >	m_m12
NTuple::Item< double >	m_m13
NTuple::Item< double >	m_m2
NTuple::Item< double >	m_m23
NTuple::Item< double >	m_m3
bool	m_Ncharge
NTuple::Item< long >	m_nchr
NTuple::Item< long >	m_nchr_e
NTuple::Item< long >	m_nchr_k
NTuple::Item< long >	m_nchr_mu
NTuple::Item< long >	m_nchr_p
NTuple::Item< long >	m_nchr_pi
bool	m_NtupleFile
int	m_numberEvent
string	m_ParentPart
string	m_PdtTable
bool	m_Psi4040OpenCharm
NTuple::Array< double >	m_px_trk
NTuple::Array< double >	m_py_trk
NTuple::Array< double >	m_pz_trk
EvtBesRandom *	m_RandomEngine
string	m_SB3File
string	m_SB3HT
vector< long int >	m_seeds
bool	m_statDecays
bool	m_tagLundModel
int	m_targetID
NTuple::Array< long >	m_Trk_index
NTuple::Tuple *	m_tuple
NTuple::Tuple *	mass_tuple
NTuple::Tuple *	massgen_tuple
int	multi
std::ofstream	outfile
std::ofstream	outfile2
IBesRndmGenSvc *	p_BesRndmGenSvc
double	parentMass
int	parentPDGcode
int	pdg
int	pdg0
int	pdg1
int	pdg2
double	px_trk [50]
double	py_trk [50]
double	pz_trk [50]
EvtVector4R	son
EvtVector4R	son0
EvtVector4R	son1
EvtVector4R	son2
EvtHis2F	SuperBody3decay
int	totalChannels
NTuple::Item< long >	TotNumTrk
int	Trk_index [50]
string	userDecFileName
int	vbr [500]
int	writeFlag

---

Constructor & Destructor Documentation

EvtDecay::EvtDecay (  const string &  name, ISvcLocator *  pSvcLocator )

:Algorithm( name, pSvcLocator ) { // these can be used to specify alternative locations or filenames // for the EvtGen particle and channel definition files. declareProperty("DecayDecDir", m_DecayDec = ""); declareProperty("PdtTableDir", m_PdtTable = ""); declareProperty("userDecayTableName", userDecFileName = "NONE"); declareProperty("DecayTopology", m_DecayTop = ""); // output decay topology to a file specified by user declareProperty("DecayRecord", m_DecayRec = ""); // output decay topology of one particle specified by user to a file declareProperty("ParentParticle", m_ParentPart = "NONE"); // Mothor particle name in pdt.table for only running BesEvtGen declareProperty("Multiplicity", m_Ncharge = false); // output ncharge number of an event declareProperty("NutupleFile",m_NtupleFile= false); // output Ntuple file declareProperty("mDIY",_mDIY= false); // mDIY model declareProperty("FDPdata",m_SB3File= ""); // Fit Dalitz Plot (FDP) data file name (*.root) declareProperty("FDPHisT",m_SB3HT= ""); // histogram title of Dalitz plot in data file declareProperty("FDPpart",m_FDPparticle=""); //to assign the FDP parent particle name (string) declareProperty("Psi3SopenCharm",m_Psi4040OpenCharm=false); declareProperty("statDecays",m_statDecays=false); declareProperty("TagLundCharmModel", m_tagLundModel=false); }

---

Member Function Documentation

void EvtDecay::assign_momentum (  int  pdx, EvtVector4R  pv4 )  [private]

{ if(id0==pdg0){son0=pd0;} else if(id0==pdg1){son1=pd0;} else if(id0==pdg2){son2=pd0;} else {std::cout<<"PID= "<<id0<<" not machted, I stop"<<std::endl; ::abort();} }

void EvtDecay::assign_momentum2 (  int  pdx, EvtVector4R  pv4 )  [private]

{ // for two identicle particle case if(id0==pdg0 && identical_flag){son0=pd0;identical_flag=false;} else if(id0==pdg1){son1=pd0;identical_flag=true;} else if(id0==pdg2){son2=pd0;} else {std::cout<<"PID not machted, I stop"<<std::endl; ::abort();} }

double EvtDecay::CalAmpsMax (  EvtParticle *  part  )  [private]

{ double amps_max=0; for(int i=0;i<100000;i++){ m_Gen->generateDecay(part); double amps=CalAmpsMDIY(part); if(amps > amps_max) amps_max=amps; } amps_max=amps_max*1.05; m_ampscalflag = false; return amps_max; }

double EvtDecay::CalAmpsMDIY (  EvtParticle *  part  )  [private]

{ //---- read out the momentum for each particle, then using them to calculate the amplitude squared, for example: /*********** using the following decay cards **************** Decay J/psi 1 gamma eta_c PHSP; Enddecay Decay eta_c 1 phi phi PHSP; Enddecay Decay phi 1 K+ K- PHSP; Enddecay End */ //---------------------------------- EvtVector4R pjsi,petac,pphi1,pphi2,pkp1,pkp2; pjsi = part->getP4(); petac= part->getDaug(1)->getP4(); //etac momentum pphi1= part->getDaug(1)->getDaug(0)->getP4(); //phi1 momentum pphi2= part->getDaug(1)->getDaug(1)->getP4(); //phi2 momentum pkp1 = part->getDaug(1)->getDaug(0)->getDaug(0)->getP4(); //K+ from phi1 pkp2 = part->getDaug(1)->getDaug(1)->getDaug(0)->getP4(); //K+ from phi2 EvtHelSys angles0(pjsi,petac); double theta0 = angles0.getHelAng(1); double phi0 = angles0.getHelAng(2); EvtHelSys angles2(pphi1,pkp1); double theta2 = angles2.getHelAng(1); double phi2 = angles2.getHelAng(2); EvtHelSys angles3(pphi2,pkp2); double theta3 = angles3.getHelAng(1); double phi3 = angles3.getHelAng(2); double amps=(3+cos(2*theta0)) * sin(theta2)*sin(theta2) * sin(theta3)*sin(theta3) * sin(phi2+phi3)*sin(phi2+phi3); if(amps <=0){ report(INFO,"EvtGen") << "Amplitude square of modified DIY should be positive, but found to be equal "<<amps<<endl; abort(); } else { return amps; } }

StatusCode EvtDecay::callBesEvtGen (  HepMC::GenEvent *  hepMCevt  )  [private]

{ MsgStream log(messageService(), name()); nchr = -1; nchr_e = 0; nchr_mu= 0; nchr_pi= 0; nchr_k = 0; nchr_p = 0; N_gamma=0; N_gammaFSR = 0; SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); for (int i=0;i<13;i++) fst[i]=0; HepMC::GenEvent::particle_const_iterator itp; // This is the main loop. // We iterate over particles and we look for ones that // should be decayed with EvtGen. // // status == 1 - undecayed Pythia particle (also for particles that are not supposed to decay) // status == 999 - particle already decayed by EvtGen // status == 899 - particle was supposed to decay by EvtGen - but found no model // this may be also intentional - such events are then excluded from // being written into persistent output. // ISTHEP(IHEP): // status code for entry IHEP, with the following meanings: // // = 0 : null entry. // = 1 : an existing entry, which has not decayed or fragmented. This is the main class of entries, which represents the `final state' given by the generator. // = 2 : an entry which has decayed or fragmented and is therefore not appearing in the final state, but is retained for event history information. // = 3 : a documentation line, defined separately from the event history. This could include the two incoming reacting particles, etc. // = 4 - 10 : undefined, but reserved for future standards. // = 11 - 200 : at the disposal of each model builder for constructs specific to his program, but equivalent to a null line in the context of any other program. // = 201 - : at the disposal of users, in particular for event tracking in the detector. loop_to_select_eventsB: EvtId ppid=EvtPDL::getId(m_ParentPart); //parent particle name double ppmass=EvtPDL::getMass(ppid); parentMass=ppmass; int pid=EvtPDL::getStdHep(ppid); parentPDGcode=pid; EvtVector4R p_init(ppmass,0.0,0.0,0.0); EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init); HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,ppmass),pid,3); bool parentFlag = false; // set spin density matrix for e+ e- -> V if(writeFlag ==0 && part->getSpinType() == EvtSpinType::VECTOR) { part->setVectorSpinDensity(); parentFlag = true; writeFlag++; } // for SuperBody3decay generator EvtVector4R fdp_init; EvtId fdp_ppid; if(m_FDPparticle!=""){ findPart(part); fdp_ppid = FDP_id; fdp_init = FDP_init; } else{ fdp_ppid = ppid; fdp_init = p_init; } if(!(m_SB3File=="" || m_SB3HT=="") && first ){ SuperBody3decay_make( fdp_ppid, fdp_init);} // ----------------------------------- loop_to_select_events: m_Gen->generateDecay(part); double ratio,rdm; if(_mDIY && m_ampscalflag ) _amps_max=CalAmpsMax(part); if(_mDIY ) { rdm=EvtRandom::Flat(0.0, 1.0); double amps=CalAmpsMDIY( part ); ratio=amps/_amps_max; } if(_mDIY && ratio<=rdm) goto loop_to_select_events; //--------- For superbody3 bool accept; if(m_FDPparticle==""){FDP_part=part;} if(!(m_SB3File=="" || m_SB3HT=="")){ accept=SuperBody3decay_judge( FDP_part); } if(!(m_SB3File=="" || m_SB3HT=="") && !accept) { delete hepMCpart; part->deleteTree(); goto loop_to_select_eventsB; }else if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && accept){ countChannel(FDP_part); } int Nchannel=part->getChannel(); if(m_statDecays && parentFlag) countChannel(part); //------------- dump the decay tree to the event model // int Nchannel=part->getChannel(); // if(Nchannel>=0) part->dumpTree(); if(parentFlag){ EvtDecayTag theTag(part); unsigned int modeTag, multiplicityTag; modeTag = theTag.getModeTag(); if(getModel(part) == "LUNDCHARM" && m_tagLundModel ){ multiplicityTag = decayType(part);} else { multiplicityTag = theTag.getMultTag() + decayType(part); } if(eventHeader == 0) std::cout << "event header unavailable: " << std::endl; if (eventHeader != 0) { eventHeader->setFlag1(modeTag); eventHeader->setFlag2(multiplicityTag); } } if(!(m_DecayRec=="")) {mydecayrec=part->writeTreeRec(m_DecayRec);std::cout<<std::endl;}; // ----------- pingrg 2008-03-25 --------- if ( log.level() <= MSG::DEBUG ) part->printTree() ; log << MSG::DEBUG << "Converting particles to HepMC" << endreq; makeHepMC(part, hepMCevt, hepMCpart); if(part->getNDaug()!=0) hepMCpart->set_status(999); part->deleteTree(); AllTrk_index=0; for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp ) { if (!(*itp)->end_vertex() && ( (*itp)->status()==899 || (*itp)->status()==999 )) (*itp)->set_status(1); HepMC::GenParticle* hepMCpart=*itp; int stat = hepMCpart->status(); int id = hepMCpart->pdg_id(); if ( stat != 1 ) continue; if( id == 22 ) {N_gamma ++;fst[0]++;} if( id == -22 ){N_gammaFSR ++;} if( id == 11 ) {fst[1]++;} else if(id == -11) {fst[2]++;} else if(id == 13 ) {fst[3]++;} else if(id ==-13) {fst[4]++;} else if(id==211) {fst[5]++;} else if(id==-211) {fst[6]++;} else if(id== 321) {fst[7]++;} else if(id ==-321) {fst[8]++;} else if(id ==2212) {fst[9]++;} else if(id==-2212) {fst[10]++;} else if(id == 2112){fst[11]++;} else if(id==-2112) {fst[12]++;} if( fabs(id) == 11) {nchr_e++;} if( fabs(id) == 13) {nchr_mu++;} if( fabs(id) == 211) {nchr_pi++;} if( fabs(id) == 321) {nchr_k++;} if( fabs(id) == 2212) {nchr_p++;} //for Nuple double en =(hepMCpart->momentum()).e(); double pex=(hepMCpart->momentum()).px(); double pey=(hepMCpart->momentum()).py(); double pez=(hepMCpart->momentum()).pz(); if( m_NtupleFile==true ){ Trk_index[AllTrk_index]=id; px_trk[AllTrk_index]=pex; py_trk[AllTrk_index]=pey; pz_trk[AllTrk_index]=pez; en_trk[AllTrk_index]=en ; /* std::cout<<"trk# " <<AllTrk_index <<" PID:" <<id <<" px: " <<pex <<" py: " <<pey <<" pz: " <<pez <<" E: " <<en<<std::endl;*/ AllTrk_index++; Trk_index[AllTrk_index]=0; } } itp=hepMCevt->particles_begin(); // parent decaying particle status set (*itp)->set_status(2); nchr = nchr_e + nchr_mu + nchr_pi + nchr_k +nchr_p; if(m_Ncharge == true ) {std::cout<<"MultiplicityTotoalCharge_e_mu_pi_K_p_gamma_gammaFSR: " <<nchr<<" " <<nchr_e<<" " << nchr_mu <<" " << nchr_pi <<" " << nchr_k <<" " <<nchr_p<<" " <<N_gamma<<" " <<N_gammaFSR<<endl;} if(m_Ncharge == true ){std::cout<<"FinalStates: " <<fst[0]<<" " <<fst[1]<<" " <<fst[2]<<" " <<fst[3]<<" " <<fst[4]<<" " <<fst[5]<<" " <<fst[6]<<" " <<fst[7]<<" " <<fst[8]<<" " <<fst[9]<<" " <<fst[10]<<" " <<fst[11]<<" " <<fst[12]<<std::endl;} if( m_NtupleFile ){ TotNumTrk=50; m_nchr = nchr; m_nchr_e = nchr_e; m_nchr_mu = nchr_mu; m_nchr_pi = nchr_pi; m_nchr_k = nchr_k; m_nchr_p = nchr_p; m_gamma=N_gamma; m_gammaFSR=N_gammaFSR; for(int i=0;i<=50;i++){ if(Trk_index[i] == 0) break; m_Trk_index[i] = Trk_index[i]; m_fst[i] = fst[i]; m_px_trk[i] = px_trk[i]; m_py_trk[i] = py_trk[i]; m_pz_trk[i] = pz_trk[i]; m_en_trk[i] = en_trk[i]; } m_tuple->write(); } //std::cout<<"Flag= "<<eventHeader->flag1()<<" "<<eventHeader->flag2()<<std::endl; return StatusCode::SUCCESS; }

StatusCode EvtDecay::callEvtGen (  HepMC::GenEvent *  hepMCevt  )  [private]

{ MsgStream log(messageService(), name()); nchr = 0; nchr_e = 0; nchr_mu= 0; nchr_pi= 0; nchr_k = 0; nchr_p = 0; N_gamma=0; N_gammaFSR = 0; SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); for (int i=0;i<13;i++) fst[i]=0; HepMC::GenEvent::particle_const_iterator itp; AllTrk_index=0; for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp ) { // This is the main loop. // We iterate over particles and we look for ones that // should be decayed with EvtGen. // // status == 1 - undecayed Pythia particle (also for particles that are not supposed to decay) // status == 999 - particle already decayed by EvtGen // status == 899 - particle was supposed to decay by EvtGen - but found no model // this may be also intentional - such events are then excluded from // being written into persistent output. // ISTHEP(IHEP): // status code for entry IHEP, with the following meanings: // // = 0 : null entry. // = 1 : an existing entry, which has not decayed or fragmented. This is the main class of entries, which represents the `final state' given by the generator. // = 2 : an entry which has decayed or fragmented and is therefore not appearing in the final state, but is retained for event history information. // = 3 : a documentation line, defined separately from the event history. This could include the two incoming reacting particles, etc. // = 4 - 10 : undefined, but reserved for future standards. // = 11 - 200 : at the disposal of each model builder for constructs specific to his program, but equivalent to a null line in the context of any other program. // = 201 - : at the disposal of users, in particular for event tracking in the detector. HepMC::GenParticle* hepMCpart=*itp; int stat = hepMCpart->status(); if ( stat != 1 ) continue; loop_to_select_eventsB: int id = hepMCpart->pdg_id(); parentPDGcode=id; hepMCpart->set_status(899); EvtId eid=EvtPDL::evtIdFromStdHep(id); string parentName=EvtPDL::name(eid); double en =(hepMCpart->momentum()).e(); double pex=(hepMCpart->momentum()).px(); double pey=(hepMCpart->momentum()).py(); double pez=(hepMCpart->momentum()).pz(); EvtVector4R p_init(en,pex,pey,pez); parentMass=p_init.mass(); EvtParticle* part=EvtParticleFactory::particleFactory(eid,p_init); // set spin density matrix for e+ e- -> V bool parentFlag=false; if( writeFlag==0 && part->getSpinType() == EvtSpinType::VECTOR) { part->setVectorSpinDensity(); parentFlag=true; writeFlag++; } else { int id = hepMCpart->pdg_id(); if( id == 22 ) {N_gamma ++;fst[0]++;} //fst[0]:photon if( id == -22 ){N_gammaFSR ++;} if( id == 11 ) {fst[1]++;} else if(id == -11){fst[2]++;} // fst[1]: electron ; fst[]: positron else if(id == 13){fst[3]++;} // fst[3]: mu- else if(id == -13){fst[4]++;} // fst[4]: mu+ else if(id == 211){fst[5]++;} // fst[5]: pi+ else if(id ==-211){fst[6]++;} // fst[6]: pi- else if(id == 321){fst[7]++;} // fst[7]: K+ else if(id ==-321){fst[8]++;} // fst[8]: K- else if(id ==2212){fst[9]++;} // fst[9]: pronton else if(id ==-2212){fst[10]++;} // fst[10]: anti-proton else if(id == 2112){fst[11]++;} // fst[11]: nucleon else if(id ==-2112){fst[12]++;} // fst[12]: ant-nucleon if( fabs(id) == 11) {nchr_e++;} if( fabs(id) == 13) {nchr_mu++;} if( fabs(id) == 211) {nchr_pi++;} if( fabs(id) == 321) {nchr_k++;} if( fabs(id) == 2212) {nchr_p++;} if( m_NtupleFile==true ){ Trk_index[AllTrk_index]=id; px_trk[AllTrk_index]=pex; py_trk[AllTrk_index]=pey; pz_trk[AllTrk_index]=pez; en_trk[AllTrk_index]=en ; AllTrk_index++; Trk_index[AllTrk_index]=0; } } // for SuperBody3decay generator // EvtVector4R pd1,pd2,pd3; EvtVector4R fdp_init; EvtId fdp_ppid; if(m_FDPparticle!=""){ findPart(part); fdp_ppid = FDP_id; fdp_init = FDP_init; } else{ fdp_ppid = eid; fdp_init = p_init; } if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && first ){ SuperBody3decay_make(fdp_ppid, fdp_init ); } loop_to_select_eventsA: m_Gen->generateDecay(part); double ratio,rdm; if(_mDIY && parentFlag && m_ampscalflag ) _amps_max=CalAmpsMax(part); if(_mDIY && parentFlag) { rdm=EvtRandom::Flat(0.0, 1.0); double amps=CalAmpsMDIY( part ); ratio=amps/_amps_max; } if(_mDIY && parentFlag && ratio<=rdm){ goto loop_to_select_eventsA;} //-------- For superbody3------------------------------ bool accept; if(m_FDPparticle==""){FDP_part=part;} if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag ){ accept=SuperBody3decay_judge( FDP_part); } if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && !accept){ part->deleteTree(); writeFlag=0; goto loop_to_select_eventsB; }else if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && accept){ countChannel(FDP_part); } //-------- for psi4040 open charm inclusive decay if(m_Psi4040OpenCharm && parentFlag ){ EvtPsi3Sdecay opencharm(en, part); bool theDecay = opencharm.choseDecay(); if(!theDecay ) { part->deleteTree(); writeFlag=0; goto loop_to_select_eventsB; } } //------------- dump the decay tree to the event model // if(Nchannel>=0) part->dumpTree(); // if(Nchannel>=0) part->printTree(); if(parentFlag){ EvtDecayTag theTag(part); unsigned int modeTag, multiplicityTag; modeTag = theTag.getModeTag(); if(getModel(part) == "LUNDCHARM" && m_tagLundModel){ multiplicityTag = decayType(part);} else { multiplicityTag = theTag.getMultTag() + decayType(part); } // std::cout<<"StdHep= "<<EvtPDL::getStdHep(part->getId())<<" getChannel "<<part->getChannel()<<" flag1,flag2= "<<modeTag<<" "<<multiplicityTag<<std::endl; if(eventHeader == 0) std::cout << "event header unavailable: " << std::endl; if (eventHeader != 0) { eventHeader->setFlag1(modeTag); eventHeader->setFlag2(multiplicityTag); // std::cout<<"flag1,flag2= "<<modeTag<<" "<<multiplicityTag<<std::endl; } } if(!(m_DecayRec=="")) mydecayrec=part->writeTreeRec(m_DecayRec); //-- decay statistics if(m_statDecays && parentFlag ) countChannel(part); // ----------- pingrg 2008-03-25 --------- if ( log.level() <= MSG::DEBUG ) part->printTree() ; log << MSG::DEBUG << "Converting particles to HepMC" << endreq; makeHepMC(part, hepMCevt, hepMCpart); if(part->getNDaug()!=0) hepMCpart->set_status(999); part->deleteTree(); }; for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp ) { if (!(*itp)->end_vertex() && ( (*itp)->status()==899 || (*itp)->status()==999 )) (*itp)->set_status(1); } nchr = nchr_e + nchr_mu + nchr_pi + nchr_k +nchr_p; if(m_Ncharge == true ) {std::cout<<"MultiplicityTotoalCharge_e_mu_pi_K_p_gamma_gammaFSR: " <<nchr <<" " <<nchr_e <<" " << nchr_mu <<" " << nchr_pi <<" " << nchr_k <<" " <<nchr_p <<" " <<N_gamma <<" " <<N_gammaFSR<<endl;} if(m_Ncharge == true ){std::cout<<"FinalStates: " <<fst[0]<<" " <<fst[1]<<" " <<fst[2]<<" " <<fst[3]<<" " <<fst[4]<<" " <<fst[5]<<" " <<fst[6]<<" " <<fst[7]<<" " <<fst[8]<<" " <<fst[9]<<" " <<fst[10]<<" " <<fst[11]<<" " <<fst[12]<<std::endl;} if( m_NtupleFile ){ TotNumTrk=50; m_nchr = nchr; m_nchr_e = nchr_e; m_nchr_mu = nchr_mu; m_nchr_pi = nchr_pi; m_nchr_k = nchr_k; m_nchr_p = nchr_p; m_gamma=N_gamma; m_gammaFSR=N_gammaFSR; for(int i=0;i<=50;i++){ if(Trk_index[i] == 0) break; m_Trk_index[i] = Trk_index[i]; m_fst[i] = fst[i]; m_px_trk[i] = px_trk[i]; m_py_trk[i] = py_trk[i]; m_pz_trk[i] = pz_trk[i]; m_en_trk[i] = en_trk[i]; } m_tuple->write(); } // std::cout<<"Flag= "<<eventHeader->flag1()<<" "<<eventHeader->flag2()<<std::endl; return StatusCode::SUCCESS; }

void EvtDecay::countChannel (  EvtParticle *  part  )  [private]

{ int ich=part->getChannel(); //std::cout<<"the channel is "<<ich<<std::endl; br[ich]++; if(ich>totalChannels) totalChannels = ich; //std::cout<<"EVENT IN br_i "<<br[ich]<<std::endl; }

int EvtDecay::decayType (  EvtParticle *  par  )  [private]

{ /************************************* * 0: gamma light_hadrons * 1: gamma charmonium * 2: DDbar (D0 D0bar or D+D-) * 3: ligh_hadrons * 4: others *************************************/ int Nchannel=par->getChannel(); int ndg = par->getNDaug(); int nfsr=0; std::string model = getModel(par,Nchannel); if(model == "LUNDCHARM" && m_tagLundModel){ // lundcharm model tag int ldm = par->getGeneratorFlag(); // std::cout<<"LundCharmFlag = "<<ldm<<std::endl; return ldm; } for(int i=0;i<ndg;i++){ //remove the FSR photon EvtId xid =par->getDaug(i)->getId(); if(EvtPDL::getStdHep(xid) == -22 ){nfsr++;} } if( isRadecay(par) ){ for(int i=0;i<ndg;i++){ EvtId xid = par->getDaug(i)->getId(); if( isCharmonium(xid) ) return 1; } return 0; } else{ if(ndg-nfsr ==2 ){ int ndg = par->getNDaug(); EvtId xd1 = par->getDaug(0)->getId(); EvtId xd2 = par->getDaug(1)->getId(); if( isCharm(xd1) && isCharm(xd2) ) {return 2;} else if( !isCharmonium(xd1) && !isCharmonium(xd2) && !isCharm(xd1) && !isCharm(xd2) ){ return 3; } } else{ // ndg>=3 bool flag = false; for(int i=0;i<ndg;i++){ EvtId xid = par->getDaug(i)->getId(); if( isCharmonium(xid) ) {flag = true; break;} if( isCharm(xid) ) {flag = true; break;} } // for_i block if( !flag ){return 3;} else {return 4;} }// if_ndg block } }

StatusCode EvtDecay::execute (   )

{ MsgStream log(messageService(), name()); // log << MSG::INFO << "EvtDecay executing" << endreq; //------------------ string key = "GEN_EVENT"; // retrieve event from Transient Store (Storegate) SmartDataPtr<McGenEventCol> McEvtColl(eventSvc(), "/Event/Gen"); m_numberEvent += 1; writeFlag = 0; //std::cout<<"EvtNumber= "<<m_numberEvent<<std::endl; if ( McEvtColl == 0 ) { HepMC::GenEvent* evt=new GenEvent(); evt->set_event_number(m_numberEvent); callBesEvtGen( evt ); if(!(m_DecayTop=="")) evt->print(outfile); //create a Transient store McGenEventCol *mcColl = new McGenEventCol; McGenEvent* mcEvent = new McGenEvent(evt); mcColl->push_back(mcEvent); StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl); if(sc != SUCCESS) return StatusCode::FAILURE; } else // (McEvtColl != 0) { McGenEventCol::iterator mcItr; for( mcItr = McEvtColl->begin(); mcItr != McEvtColl->end(); mcItr++ ) { HepMC::GenEvent* hEvt = (*mcItr)->getGenEvt(); // MeVToGeV( hEvt ); callEvtGen( hEvt ); if(!(m_DecayTop=="")) hEvt->print(outfile); // GeVToMeV( hEvt ); // if(!(m_DecayRec=="")) outfile2<<std::endl; if(!(m_DecayRec=="")) std::cout<<std::endl; }; } return StatusCode::SUCCESS; }

StatusCode EvtDecay::finalize (   )

{ MsgStream log(messageService(), name()); log << MSG::INFO << "EvtDecay finalized" << endreq; fstream Forfile; Forfile.open("fort.16",ios::in); char delcmd[300]; strcpy(delcmd,"rm -f "); strcat(delcmd,"fort.16"); if(Forfile)system(delcmd); fstream Forfile2; Forfile2.open("fort.22",ios::in); strcpy(delcmd,"rm -f "); strcat(delcmd,"fort.22"); if(Forfile2)system(delcmd); // To get the branching fraction of the specified mode in Lundcharm Model /* EvtLundCharm lundcharmEvt; int nevt=lundcharmEvt.getTotalEvt(); std::cout<<"The total number of the specified mode is "<<nevt<<std::endl; */ int totalEvt=0; if(!(m_SB3File=="" || m_SB3HT=="")){ for(int i=0;i<500;i++){totalEvt=totalEvt+br[i];} for(int i=0;i<500;i++){ double bi=1.0*br[i]/totalEvt; std::cout<<"Branching fraction of channel "<<i<<" is: "<<bi<<std::endl; if(br[i]==0) break; } } if(m_statDecays){ int totalEvt=0; for(int i=0;i<=totalChannels;i++){totalEvt=totalEvt+br[i];} std::cout<<"==================Statistical first chain decay ==============================="<<std::endl; std::cout<<"i-th channel Events Generated Branching fraction generated "<<std::endl; for(int i=0;i<=totalChannels;i++){ std::cout<<"| "<<i<<" "<<br[i]<<" "<<br[i]*1.0/totalEvt<<std::endl; } std::cout<<"--------------------------------------------------------------------------------"<<std::endl; std::cout<<" sum: "<<totalEvt<<" "<<"1"<<std::endl; std::cout<<"================== End of statistical first chain decay ========================"<<std::endl; } return StatusCode::SUCCESS; }

void EvtDecay::FinalState_make (  EvtId  ppid, EvtVector4R  p_init )  [private]

{ // get the final state pdg cond and count the identicle particle multi=1; identical_flag=true; for(int i=1;i<10000;i++){ //sigle out the final state EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init); HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,parentMass),parentPDGcode,3); m_Gen->generateDecay(part); int nson=part->getNDaug(); if(nson == 2) {continue;} else if(nson ==3){ EvtId xid0,xid1,xid2; xid0=part->getDaug(0)->getId(); xid1=part->getDaug(1)->getId(); xid2=part->getDaug(2)->getId(); //-- pdg code pdg0=EvtPDL::getStdHep(xid0); pdg1=EvtPDL::getStdHep(xid1); pdg2=EvtPDL::getStdHep(xid2); if(pdg0==pdg1 && pdg1==pdg2) {multi=3;} else if(pdg0==pdg1 ){multi=2;} // two identical particle list as 0,1 index else if(pdg0==pdg2 ){multi=2;pdg=pdg1; pdg1=pdg2; pdg2=pdg;} else if(pdg1==pdg2) {multi=2;pdg=pdg0; pdg0=pdg1; pdg1=pdg2; pdg2=pdg;} else {multi=1;} break; } else{ std::cout<<"No direct three body decay channel found in decay card, I stop run"<<std::endl; ::abort(); } } // std::cout<<"pdg_i "<<pdg0<<" "<<pdg1<<" "<<pdg2<<std::endl; // std::cout<<"identical particle "<<multi<<std::endl; }

void EvtDecay::FinalState_sort (  EvtParticle *  part  )  [private]

{ // sort the momentum to aon0, son1, son2 EvtVector4R pd0,pd1,pd2; EvtId xid0,xid1,xid2; int id0,id1,id2; int nson=part->getNDaug(); if(nson==2){ pd0=part->getDaug(0)->getP4Lab(); pd1=part->getDaug(1)->getDaug(0)->getP4Lab(); pd2=part->getDaug(1)->getDaug(1)->getP4Lab(); xid0=part->getDaug(0)->getId(); xid1=part->getDaug(1)->getDaug(0)->getId(); xid2=part->getDaug(1)->getDaug(1)->getId(); } else if(nson==3){ pd0=part->getDaug(0)->getP4Lab(); pd1=part->getDaug(1)->getP4Lab(); pd2=part->getDaug(2)->getP4Lab(); xid0=part->getDaug(0)->getId(); xid1=part->getDaug(1)->getId(); xid2=part->getDaug(2)->getId(); } //-- pdg code id0=EvtPDL::getStdHep(xid0); id1=EvtPDL::getStdHep(xid1); id2=EvtPDL::getStdHep(xid2); // std::cout<<"pid before match: "<<id0<<" "<<id1<<" "<<id2<<std::endl; //-- assign sons momentum if(multi==1){ assign_momentum(id0,pd0); assign_momentum(id1,pd1); assign_momentum(id2,pd2); } else if(multi==2){ assign_momentum2(id0,pd0); assign_momentum2(id1,pd1); assign_momentum2(id2,pd2); if(son0.get(0)>son1.get(0)){son=son0;son0=son1;son1=son;} // change into E_0 < E_1 } else if(multi==3){ // sort sons according to energy E_0 < E_1 < E_2 double En0=pd0.get(0); double En1=pd1.get(0); double En2=pd2.get(0); if(En0 < En1 && En1 < En2) {son0=pd0;son1=pd1;son2=pd2;} if(En0 < En2 && En2 < En1) {son0=pd0;son1=pd2;son2=pd1;} if(En1 < En0 && En0 < En2) {son0=pd1;son1=pd0;son2=pd2;} if(En1 < En2 && En2 < En0) {son0=pd1;son1=pd2;son2=pd0;} if(En2 < En0 && En0 < En1) {son0=pd2;son1=pd0;son2=pd1;} if(En2 < En1 && En1 < En0) {son0=pd2;son1=pd1;son2=pd0;} } }

void EvtDecay::findPart (  EvtParticle *  part  )  [private]

{ FDP_id=EvtPDL::getId(m_FDPparticle); int FDP_pdg=EvtPDL::getStdHep(FDP_id); m_Gen->generateDecay(part); int dn=part->getNDaug(); if(dn!=0){ for(int i=0;i<dn;i++){ EvtParticle* part1=part->getDaug(i); EvtId sonid=part1->getId(); int son_pdg = EvtPDL::getStdHep(sonid); if(son_pdg == FDP_pdg){ FDP_init=part1->getP4Lab(); FDP_part=part1; return; } else{ findPart(part1); } } } //if (dn block }

std::string EvtDecay::getModel (  EvtParticle *  par  )  [private]

{ int mode = par->getChannel(); EvtDecayBase* thedecay = EvtDecayTable::gettheDecay(par->getId(),mode); return thedecay->getModelName(); }

std::string EvtDecay::getModel (  EvtParticle *  par, int  mode )  [private]

{ EvtDecayBase* thedecay = EvtDecayTable::gettheDecay(par->getId(),mode); return thedecay->getModelName(); }

void EvtDecay::GeVToMeV (  HepMC::GenEvent *  hepMCevt  )  [private]

{ for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) { // cout << " PDG, BAR " << (*p)->pdg_id() << " " << (*p)->barcode() << endl; (*p)->set_momentum( (*p)->momentum() * 1000. ); } }

StatusCode EvtDecay::initialize (   )

{ MsgStream log(messageService(), name()); log << MSG::INFO << "EvtDecay initialize" << endreq; static const bool CREATEIFNOTTHERE(true); StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE); if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc) { log << MSG::ERROR << " Could not initialize Random Number Service" << endreq; return RndmStatus; } m_numberEvent=0; first = true; m_ampscalflag = true; // Save the random number seeds in the event HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("EVTGEN"); const long s = engine->getSeed(); p_BesRndmGenSvc->setGenseed(s+1); m_seeds.clear(); m_seeds.push_back(s); totalChannels = 0; // open an interface to EvtGen if ( m_DecayDec == "") { //use default DECAY.DEC and pdt.table if not defined by user m_DecayDec=getenv("BESEVTGENROOT"); m_DecayDec +="/share/DECAY.DEC"; } if ( m_PdtTable == "") { m_PdtTable =getenv("BESEVTGENROOT"); m_PdtTable +="/share/pdt.table"; } char catcmd[300]; //output user decay cards to log file std::cout<<"===================== user decay card ========================"<<std::endl; strcpy(catcmd, "cat "); strcat(catcmd, userDecFileName.c_str()); system(catcmd); // write decay cards to the root file: pingrg@2009-09-09 StatusCode status; status = service("DataInfoSvc",tmpInfoSvc); if (status.isSuccess()) { log << MSG::INFO << "got the DataInfoSvc" << endreq; dataInfoSvc=dynamic_cast<DataInfoSvc *>(tmpInfoSvc); dataInfoSvc->setDecayCard(userDecFileName); } else { log << MSG::ERROR << "could not get the DataInfoSvc" << endreq; return StatusCode::FAILURE; } m_RandomEngine = new EvtBesRandom(engine); m_Gen = new EvtGen( m_DecayDec.c_str(), m_PdtTable.c_str(), m_RandomEngine ); if(userDecFileName =="NONE") log << MSG::INFO << "EvtDecay User did not define his Decay table EvtGen will use standart" << endreq; if(!(userDecFileName =="NONE")) m_Gen->readUDecay(userDecFileName.c_str()); if(!(m_DecayTop==" ")) { outfile.open(m_DecayTop.c_str()); } //for output Ntuple file, pingrg-2009-09-07 if(m_NtupleFile) { NTuplePtr nt1(ntupleSvc(),"MYROOTFILE/Trk"); if(nt1) {m_tuple=nt1;} else { m_tuple = ntupleSvc()->book ("MYROOTFILE/Trk", CLID_ColumnWiseTuple, "Generator-trk-Ntuple"); if(m_tuple){ status = m_tuple->addItem ("TotNumTrk", TotNumTrk, 0,100); status = m_tuple->addIndexedItem ("Trk_index", TotNumTrk, m_Trk_index); status = m_tuple->addIndexedItem ("m_px_trk", TotNumTrk, m_px_trk); status = m_tuple->addIndexedItem ("m_py_trk", TotNumTrk, m_py_trk); status = m_tuple->addIndexedItem ("m_pz_trk", TotNumTrk, m_pz_trk); status = m_tuple->addIndexedItem ("m_en_trk", TotNumTrk, m_en_trk); status = m_tuple->addIndexedItem ("FST", TotNumTrk, m_fst); status = m_tuple->addItem ("nchr", m_nchr); status = m_tuple->addItem ("nchr_e", m_nchr_e); status = m_tuple->addItem ("nchr_mu", m_nchr_mu); status = m_tuple->addItem ("nchr_pi", m_nchr_pi); status = m_tuple->addItem ("nchr_k", m_nchr_k); status = m_tuple->addItem ("nchr_p", m_nchr_p); status = m_tuple->addItem ("N_gamma", m_gamma); status = m_tuple->addItem ("N_gammaFSR", m_gammaFSR); } else { log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg; return StatusCode::FAILURE; } } NTuplePtr nt2(ntupleSvc(),"MYROOTFILE/mass"); if(nt2) {mass_tuple=nt2;} else { mass_tuple = ntupleSvc()->book ("MYROOTFILE/mass", CLID_ColumnWiseTuple, "Generator-mass-Ntuple"); if(mass_tuple){ status = mass_tuple->addItem ("m12", m_m12); status = mass_tuple->addItem ("m13", m_m13); status = mass_tuple->addItem ("m23", m_m23); status = mass_tuple->addItem ("m1", m_m1); status = mass_tuple->addItem ("m2", m_m2); status = mass_tuple->addItem ("m3", m_m3); status = mass_tuple->addItem ("costheta1", m_cos1); status = mass_tuple->addItem ("costheta2", m_cos2); status = mass_tuple->addItem ("costheta3", m_cos3); status = mass_tuple->addItem ("ichannel", m_ich); } else { log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg; return StatusCode::FAILURE; } } NTuplePtr nt3(ntupleSvc(),"MYROOTFILE/massGen"); if(nt3) {massgen_tuple=nt3;} else { massgen_tuple = ntupleSvc()->book ("MYROOTFILE/massGen", CLID_ColumnWiseTuple, "Generator-mass-Ntuple"); if(massgen_tuple){ status = massgen_tuple->addItem ("m12", _m12); status = massgen_tuple->addItem ("m13", _m13); status = massgen_tuple->addItem ("m23", _m23); status = massgen_tuple->addItem ("m1", _m1); status = massgen_tuple->addItem ("m2", _m2); status = massgen_tuple->addItem ("m3", _m3); status = massgen_tuple->addItem ("costheta1", _cos1); status = massgen_tuple->addItem ("costheta2", _cos2); status = massgen_tuple->addItem ("costheta3", _cos3); status = massgen_tuple->addItem ("ichannel", _ich); } else { log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg; return StatusCode::FAILURE; } } } for(int i=0;i<500;i++){br[i]=0;} if(!(m_SB3File=="" && m_SB3HT=="")) { const char *filename, *histitle; filename=(char*)m_SB3File.c_str(); histitle=(char*)m_SB3HT.c_str(); SuperBody3decay.setFile(filename); SuperBody3decay.setHTitle(histitle); SuperBody3decay.init(); } return StatusCode::SUCCESS; }

bool EvtDecay::isCharm (  EvtId  xid  )  [private]

{ EvtId d0 = EvtPDL::getId(std::string("D0")); EvtId d0bar = EvtPDL::getId(std::string("anti-D0")); EvtId dp = EvtPDL::getId(std::string("D+")); EvtId dm = EvtPDL::getId(std::string("D-")); EvtId d0h = EvtPDL::getId(std::string("D0H")); EvtId d0l = EvtPDL::getId(std::string("D0L")); EvtId dstp = EvtPDL::getId(std::string("D*+")); EvtId dstm = EvtPDL::getId(std::string("D*-")); EvtId ds0 = EvtPDL::getId(std::string("D*0")); EvtId ds0bar = EvtPDL::getId(std::string("anti-D*0")); EvtId dsp = EvtPDL::getId(std::string("D_s+")); EvtId dsm = EvtPDL::getId(std::string("D_s-")); EvtId dsstp = EvtPDL::getId(std::string("D_s*+")); EvtId dsstm = EvtPDL::getId(std::string("D_s*-")); EvtId ds0stp = EvtPDL::getId(std::string("D_s0*+")); EvtId ds0stm = EvtPDL::getId(std::string("D_s0*-")); std::vector<EvtId> Vid; Vid.clear(); Vid.push_back(d0); Vid.push_back(d0bar); Vid.push_back(dp); Vid.push_back(dm); Vid.push_back(d0h); Vid.push_back(d0l); Vid.push_back(dstp); Vid.push_back(dstm); Vid.push_back(ds0); Vid.push_back(ds0bar ); Vid.push_back(dsp ); Vid.push_back(dsm ); Vid.push_back(dsstp ); Vid.push_back(dsstm ); Vid.push_back(ds0stp ); Vid.push_back(ds0stm ); bool flag=true; for(int i=0;i<Vid.size();i++){ if(xid == Vid[i]) return flag;} return false; }

bool EvtDecay::isCharmonium (  EvtId  xid  )  [private]

{ EvtId psi4415 = EvtPDL::getId(std::string("psi(4415)")); EvtId psi4160 = EvtPDL::getId(std::string("psi(4160)")); EvtId psi4040 = EvtPDL::getId(std::string("psi(4040)")); EvtId psi3770 = EvtPDL::getId(std::string("psi(3770)")); EvtId psiprim = EvtPDL::getId(std::string("psi(2S)")); EvtId jpsi = EvtPDL::getId(std::string("J/psi")); EvtId etac = EvtPDL::getId(std::string("eta_c")); EvtId etac2s = EvtPDL::getId(std::string("eta_c(2S)")); EvtId hc = EvtPDL::getId(std::string("h_c")); EvtId chic0 = EvtPDL::getId(std::string("chi_c0")); EvtId chic1 = EvtPDL::getId(std::string("chi_c1")); EvtId chic2 = EvtPDL::getId(std::string("chi_c2")); EvtId chic0p = EvtPDL::getId(std::string("chi_c0p")); EvtId chic1p = EvtPDL::getId(std::string("chi_c1p")); EvtId chic2p = EvtPDL::getId(std::string("chi_c1p")); std::vector<EvtId> Vid; Vid.clear(); Vid.push_back(psi4415); Vid.push_back(psi4160); Vid.push_back(psi4040); Vid.push_back(psi3770); Vid.push_back(psiprim); Vid.push_back(jpsi); Vid.push_back(etac); Vid.push_back(etac2s); Vid.push_back(hc); Vid.push_back(chic0); Vid.push_back(chic1); Vid.push_back(chic2); Vid.push_back(chic0p); Vid.push_back(chic1p); Vid.push_back(chic2p); bool flag=true; for(int i=0;i<Vid.size();i++){ if(xid == Vid[i]) return flag;} return false; }

bool EvtDecay::isRadecay (  EvtParticle *  par  )  [private]

{ int ndg = par->getNDaug(); for(int i=0;i<ndg;i++){ EvtId xid = par->getDaug(i)->getId(); if(EvtPDL::getStdHep(xid) == 22){return true;} } return false; }

StatusCode EvtDecay::makeHepMC (  EvtParticle *  , HepMC::GenEvent *  , HepMC::GenParticle *  )  [private]

{ MsgStream log(messageService(), name()); if(part->getNDaug()!=0) { double ct=(part->getDaug(0)->get4Pos()).get(0); double x=(part->getDaug(0)->get4Pos()).get(1); double y=(part->getDaug(0)->get4Pos()).get(2); double z=(part->getDaug(0)->get4Pos()).get(3); HepMC::GenVertex* end_vtx = new HepMC::GenVertex(HepLorentzVector(x,y,z,ct)); hEvt->add_vertex( end_vtx ); end_vtx->add_particle_in(hPart); int ndaug=part->getNDaug(); for(int it=0;it<ndaug;it++) { double e=(part->getDaug(it)->getP4Lab()).get(0); double px=(part->getDaug(it)->getP4Lab()).get(1); double py=(part->getDaug(it)->getP4Lab()).get(2); double pz=(part->getDaug(it)->getP4Lab()).get(3); int id=EvtPDL::getStdHep(part->getDaug(it)->getId()); int status=1; if(part->getDaug(it)->getNDaug()!=0) status=999; HepMC::GenParticle* prod_part = new HepMC::GenParticle(HepLorentzVector(px,py,pz,e),id,status); end_vtx->add_particle_out(prod_part); makeHepMC(part->getDaug(it),hEvt,prod_part); if( log.level()<MSG::INFO ) prod_part->print(); } } return StatusCode::SUCCESS; }

void EvtDecay::MeVToGeV (  HepMC::GenEvent *  hepMCevt  )  [private]

{ for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) { // cout << " PDG, BAR " << (*p)->pdg_id() << " " << (*p)->barcode() << endl; (*p)->set_momentum( (*p)->momentum() / 1000. ); } }

bool EvtDecay::SuperBody3decay_judge (  EvtParticle *  part  )  [private]

{ double xmass2,ymass2; bool accept; EvtVector4R pd1,pd2,pd3; FinalState_sort( part); pd1=son0; pd2=son1; pd3=son2; xmass2=(pd1+pd2).mass2(); // Dalitz plot m12^2 ~ m13^2 ymass2=(pd1+pd3).mass2(); accept=SuperBody3decay.AR(xmass2,ymass2); //debugging if(accept && m_NtupleFile) { _ich=part->getChannel(); _m1=pd1.mass(); _m2=pd2.mass(); _m3=pd3.mass(); _m12=(pd1+pd2).mass(); _m23=(pd2+pd3).mass(); _m13=(pd1+pd3).mass(); _cos1=pd1.get(3)/pd1.d3mag(); _cos2=pd2.get(3)/pd2.d3mag(); _cos3=pd3.get(3)/pd3.d3mag(); massgen_tuple->write(); } // std::cout<<"mass 1 2 3 chicj>> "<<_m1<<" "<<_m2<<" "<<_m3<<" "<<_m23<<std::endl; return accept; }

void EvtDecay::SuperBody3decay_make (  EvtId  ppid, EvtVector4R  p_init )  [private]

{ double xmass2,ymass2; bool accept; EvtVector4R pd1,pd2,pd3; //---- find out the pdg codes of final state and count the identical particles FinalState_make( ppid, p_init); //begin to loop with events for(int i=0;i<100000;i++){ regen: EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init); HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,parentMass),parentPDGcode,3); //this line make the decay to select different channel if(part->getSpinType() == EvtSpinType::VECTOR) { part->setVectorSpinDensity(); } m_Gen->generateDecay(part); FinalState_sort(part); pd1=son0; pd2=son1; pd3=son2; // std::cout<<"multi, pdg, pi= "<<multi<<" "<<pdg0<<" "<<pdg1<<" "<<pdg2<<" "<<son0<<son1<<son2<<std::endl; xmass2=(pd1+pd2).mass2(); // Dalitz plot m12^2 ~ m13^2 ymass2=(pd1+pd3).mass2(); double xlow=SuperBody3decay.getXlow(); double xup=SuperBody3decay.getXup(); double ylow=SuperBody3decay.getYlow(); double yup=SuperBody3decay.getYup(); if(xmass2<xlow || xmass2>xup || ymass2<ylow || ymass2>yup) { part->deleteTree();goto regen;} SuperBody3decay.HFill(xmass2,ymass2); if( m_NtupleFile ){ m_ich=part->getChannel(); m_m1=pd1.mass(); m_m2=pd2.mass(); m_m3=pd3.mass(); m_m12=(pd1+pd2).mass(); m_m23=(pd2+pd3).mass(); m_m13=(pd1+pd3).mass(); m_cos1=pd1.get(3)/pd1.d3mag(); m_cos2=pd2.get(3)/pd2.d3mag(); m_cos3=pd3.get(3)/pd3.d3mag(); mass_tuple->write(); } double m1=pd1.mass(); double m2=pd2.mass();double m3=pd3.mass(); double mj=(pd2+pd3).mass(); // std::cout<<"mass 1 2 3 chicj <<=="<<m1<<" "<<m2<<" "<<m3<<" "<<mj<<std::endl; // delete hepMCpart; } SuperBody3decay.HReweight(); //reweighting Dalitz plot SuperBody3decay.setZmax(); // find out the maximum value after reweighting first=false; }

---

Member Data Documentation

NTuple::Item<double> EvtDecay::_cos1 [private]

NTuple::Item<double> EvtDecay::_cos2 [private]

NTuple::Item<double> EvtDecay::_cos3 [private]

NTuple::Item<long> EvtDecay::_ich [private]

NTuple::Item<double> EvtDecay::_m1 [private]

NTuple::Item<double> EvtDecay::_m12 [private]

NTuple::Item<double> EvtDecay::_m13 [private]

NTuple::Item<double> EvtDecay::_m2 [private]

NTuple::Item<double> EvtDecay::_m23 [private]

NTuple::Item<double> EvtDecay::_m3 [private]

bool EvtDecay::_mDIY [private]

int EvtDecay::AllTrk_index [private]

int EvtDecay::br[500] [private]

DataInfoSvc* EvtDecay::dataInfoSvc

double EvtDecay::en_trk[50] [private]

EvtId EvtDecay::FDP_id [private]

EvtVector4R EvtDecay::FDP_init [private]

EvtParticle* EvtDecay::FDP_part [private]

bool EvtDecay::first [private]

bool EvtDecay::identical_flag [private]

bool EvtDecay::m_ampscalflag [private]

NTuple::Item<double> EvtDecay::m_cos1 [private]

NTuple::Item<double> EvtDecay::m_cos2 [private]

NTuple::Item<double> EvtDecay::m_cos3 [private]

string EvtDecay::m_DecayDec [private]

string EvtDecay::m_DecayRec [private]

string EvtDecay::m_DecayTop [private]

NTuple::Array<double> EvtDecay::m_en_trk [private]

string EvtDecay::m_FDPparticle [private]

NTuple::Array<long> EvtDecay::m_fst [private]

NTuple::Item<long> EvtDecay::m_gamma [private]

NTuple::Item<long> EvtDecay::m_gammaFSR [private]

EvtGen* EvtDecay::m_Gen [private]

NTuple::Item<long> EvtDecay::m_ich [private]

std::vector<int> EvtDecay::m_InSeeds [private]

EvtId EvtDecay::m_KKMCRes [private]

NTuple::Item<double> EvtDecay::m_m1 [private]

NTuple::Item<double> EvtDecay::m_m12 [private]

NTuple::Item<double> EvtDecay::m_m13 [private]

NTuple::Item<double> EvtDecay::m_m2 [private]

NTuple::Item<double> EvtDecay::m_m23 [private]

NTuple::Item<double> EvtDecay::m_m3 [private]

bool EvtDecay::m_Ncharge [private]

NTuple::Item<long> EvtDecay::m_nchr [private]

NTuple::Item<long> EvtDecay::m_nchr_e [private]

NTuple::Item<long> EvtDecay::m_nchr_k [private]

NTuple::Item<long> EvtDecay::m_nchr_mu [private]

NTuple::Item<long> EvtDecay::m_nchr_p [private]

NTuple::Item<long> EvtDecay::m_nchr_pi [private]

bool EvtDecay::m_NtupleFile [private]

int EvtDecay::m_numberEvent [private]

string EvtDecay::m_ParentPart [private]

string EvtDecay::m_PdtTable [private]

bool EvtDecay::m_Psi4040OpenCharm [private]

NTuple::Array<double> EvtDecay::m_px_trk [private]

NTuple::Array<double> EvtDecay::m_py_trk [private]

NTuple::Array<double> EvtDecay::m_pz_trk [private]

EvtBesRandom* EvtDecay::m_RandomEngine [private]

string EvtDecay::m_SB3File [private]

string EvtDecay::m_SB3HT [private]

vector<long int> EvtDecay::m_seeds [private]

bool EvtDecay::m_statDecays [private]

bool EvtDecay::m_tagLundModel [private]

int EvtDecay::m_targetID [private]

NTuple::Array<long> EvtDecay::m_Trk_index [private]

NTuple::Tuple* EvtDecay::m_tuple [private]

NTuple::Tuple* EvtDecay::mass_tuple [private]

NTuple::Tuple* EvtDecay::massgen_tuple [private]

int EvtDecay::multi [private]

std::ofstream EvtDecay::outfile [private]

std::ofstream EvtDecay::outfile2 [private]

IBesRndmGenSvc* EvtDecay::p_BesRndmGenSvc [private]

double EvtDecay::parentMass [private]

int EvtDecay::parentPDGcode [private]

int EvtDecay::pdg [private]

int EvtDecay::pdg0 [private]

int EvtDecay::pdg1 [private]

int EvtDecay::pdg2 [private]

double EvtDecay::px_trk[50] [private]

double EvtDecay::py_trk[50] [private]

double EvtDecay::pz_trk[50] [private]

EvtVector4R EvtDecay::son [private]

EvtVector4R EvtDecay::son0 [private]

EvtVector4R EvtDecay::son1 [private]

EvtVector4R EvtDecay::son2 [private]

EvtHis2F EvtDecay::SuperBody3decay [private]

IDataInfoSvc* EvtDecay::tmpInfoSvc

int EvtDecay::totalChannels [private]

NTuple::Item<long> EvtDecay::TotNumTrk [private]

int EvtDecay::Trk_index[50] [private]

string EvtDecay::userDecFileName [private]

int EvtDecay::vbr[500] [private]

int EvtDecay::writeFlag [private]

---

The documentation for this class was generated from the following files:

• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/EvtGen_i/EvtDecay.h
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtDecay.cxx
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/user/UsermDIY.cc

---