21 #include "EvtGenBase/EvtAbsRadCorr.hh"
22 #include "EvtGenBase/EvtDecayBase.hh"
23 #include "EvtGen_i/EvtGenExternal/EvtExternalGenList.hh"
25 #include "EvtGenBase/EvtVector4R.hh"
26 #include "EvtGenBase/EvtParticle.hh"
27 #include "EvtGenBase/EvtParticleFactory.hh"
28 #include "EvtGen/EvtGen.hh"
29 #include "EvtGenBase/EvtRandomEngine.hh"
30 #include "EvtGenBase/EvtDecayTable.hh"
40 #include "GaudiKernel/MsgStream.h"
41 #include "GaudiKernel/ISvcLocator.h"
42 #include "GaudiKernel/DataSvc.h"
43 #include "GaudiKernel/IPartPropSvc.h"
47 #include "HepPID/ParticleName.hh"
50 #include "CLHEP/Random/RandFlat.h"
51 #include "CLHEP/Vector/LorentzVector.h"
125 msg(MSG::INFO) <<
"EvtInclusiveDecay initialize" <<
endmsg;
130 msg(MSG::INFO) <<
"EvtInclusiveDecay selection parameters:" <<
endmsg;
137 msg(MSG::INFO) <<
"User selection parameters:" <<
endmsg;
149 msg(MSG::INFO) <<
"* blackList; = ";
151 msg(MSG::INFO) << (*i) <<
" ";
156 msg(MSG::INFO) <<
"* whiteList = ";
158 msg(MSG::INFO) << (*i) <<
" ";
166 EvtExternalGenList genList(
true,
xmlpath(),
"gamma");
167 EvtAbsRadCorr* radCorrEngine = genList.getPhotosModel();
168 std::list<EvtDecayBase*> extraModels = genList.getListOfModels();
181 return StatusCode::SUCCESS;
186 const EventContext& ctx)
195 const EventContext& ctx)
const
198 rngWrapper->
setSeed(
streamName, ctx.slot(), randomSeedOffset, ctx.eventID().run_number() );
208 ctx.setEventID (EventIDBase (conditionsRun,
209 EventIDBase::UNDEFEVT,
210 EventIDBase::UNDEFNUM,
211 EventIDBase::UNDEFNUM,
223 const EventContext& ctx = Gaudi::Hive::currentContext();
247 HepMC::GenEvent* hepMC = *mcItr;
252 std::set<HepMC::GenVertexPtr> visited;
254 std::set<HepMC::GenParticlePtr> toBeDecayed;
255 for (
auto p: hepMC->particles()) {
256 if ( (!
p->production_vertex()) ||
257 (
p->production_vertex()->particles_in().size() == 0) ) {
260 return StatusCode::FAILURE;
265 msg(MSG::INFO) <<
"Printing HepMC record at " << hepMC <<
" BEFORE running EvtGen:" <<
endmsg;
272 std::set<int> toBeDecayed;
273 for (HepMC::GenEvent::particle_iterator
itp = hepMC->particles_begin();
itp != hepMC->particles_end(); ++
itp) {
275 if ( (!
p->production_vertex()) ||
276 (
p->production_vertex()->particles_in_size() == 0) ) {
279 return StatusCode::FAILURE;
284 msg(MSG::INFO) <<
"Printing HepMC record at " << hepMC <<
" BEFORE running EvtGen:" <<
endmsg;
293 bool eventPassesCuts(
false);
297 for (
auto p: toBeDecayed) {
299 msg(MSG::ERROR ) <<
"Overlapping decay tree for particle" <<
p <<
endmsg;
300 return StatusCode::FAILURE;
307 auto p = hepMC->barcode_to_particle(*itb);
309 msg(MSG::ERROR ) <<
"Overlapping decay tree encountered for barcode " << *itb <<
endmsg;
310 return StatusCode::FAILURE;
320 eventPassesCuts =
true;
329 hepMC->weight(
"nEvtGenDecayAttempts") = loopCounter;
332 hepMC->weights()[
"nEvtGenDecayAttempts"] = loopCounter;
336 msg(MSG::INFO) <<
"Printing HepMC record at " << hepMC <<
" AFTER running EvtGen:" <<
endmsg;
348 return StatusCode::SUCCESS;
355 ATH_MSG_INFO(
"The following particles were checked and didn't have any decay channels:");
357 std::cout << std::endl;
358 std::cout <<
" Particle code Name from HepPDT # Occurences" << std::endl;
359 std::cout <<
"------------------------------------------------------" << std::endl;
363 std::cout << std::setw(14) <<
id
365 << std::setw(20) <<
count
368 std::cout << std::endl;
373 return StatusCode::SUCCESS;
388 std::set<HepMC::GenVertexPtr>& visited,
389 std::set<HepMC::GenParticlePtr>& toBeDecayed) {
393 std::set<HepMC::GenVertexPtr>& visited,
394 std::set<int>& toBeDecayed) {
397 if (!isToBeRemoved) {
400 toBeDecayed.insert(
p);
404 isToBeRemoved =
true;
414 auto v =
p->end_vertex();
417 if (visited.insert(
v).second) {
419 ATH_MSG_WARNING(
"Found particle to be decayed with vertex with >1 incoming mother particles in decay tree");
422 for (
auto itp:
v->particles_out()) {
429 if (visited.insert(
v).second) {
432 ATH_MSG_WARNING(
"Found particle to be decayed with vertex with >1 incoming mother particles in decay tree");
441 return StatusCode::SUCCESS;
450 auto v =
p->end_vertex();
454 hepMC->remove_vertex(
v);
458 std::set<int> vtxBarCodesToDelete;
459 vtxBarCodesToDelete.insert(
v->barcode());
460 for (HepMC::GenVertex::vertex_iterator itv =
v->vertices_begin(HepMC::descendants);
461 itv !=
v->vertices_end(HepMC::descendants);
463 vtxBarCodesToDelete.insert((*itv)->barcode());
465 auto vdel = hepMC->barcode_to_vertex(*itb);
466 hepMC->remove_vertex(vdel);
471 <<
" decay tree with " << vtxBarCodesToDelete.size() <<
" vertices");
504 int id =
part->pdg_id();
505 EvtId evtId=EvtPDL::evtIdFromStdHep(
id);
506 double en =(
part->momentum()).
e()/1000.;
507 double px=(
part->momentum()).
px()/1000.;
508 double py=(
part->momentum()).
py()/1000.;
509 double pz=(
part->momentum()).
pz()/1000.;
511 EvtParticle* evtPart = EvtParticleFactory::particleFactory(evtId,evtP);
514 if(
m_setVMtransversePol && (
id==113 ||
id== 443 ||
id==100443 ||
id==553 ||
id==100553 ||
id==200553) )evtPart->setVectorSpinDensity();
518 double ct_s =
part->production_vertex()->position().t();
519 double x_s =
part->production_vertex()->position().x();
520 double y_s =
part->production_vertex()->position().y();
521 double z_s =
part->production_vertex()->position().z();
523 EvtVector4R treeStart(ct_s,x_s,y_s,z_s);
526 if(evtPart->getNDaug() !=0)
part->set_status(2);
527 evtPart->deleteTree();
533 EvtParticle* evtPart, EvtVector4R treeStart,
double momentumScaleFactor) {
534 if(evtPart->getNDaug()!=0) {
536 double ct=(evtPart->getDaug(0)->get4Pos()).
get(0)+treeStart.get(0);
537 double x=(evtPart->getDaug(0)->get4Pos()).
get(1)+treeStart.get(1);
538 double y=(evtPart->getDaug(0)->get4Pos()).
get(2)+treeStart.get(2);
539 double z=(evtPart->getDaug(0)->get4Pos()).
get(3)+treeStart.get(3);
543 hepMC->add_vertex(end_vtx);
544 end_vtx->add_particle_in(
part);
547 for(
uint it=0;
it<evtPart->getNDaug();
it++) {
548 double e=(evtPart->getDaug(
it)->getP4Lab()).
get(0) * momentumScaleFactor;
549 double px=(evtPart->getDaug(
it)->getP4Lab()).
get(1) * momentumScaleFactor;
550 double py=(evtPart->getDaug(
it)->getP4Lab()).
get(2) * momentumScaleFactor;
551 double pz=(evtPart->getDaug(
it)->getP4Lab()).
get(3) * momentumScaleFactor;
552 int id=EvtPDL::getStdHep(evtPart->getDaug(
it)->getId());
554 if(evtPart->getDaug(
it)->getNDaug() != 0)
status=2;
556 end_vtx->add_particle_out(daughter);
571 int id =
p->pdg_id();
573 auto v =
p->end_vertex();
574 if (
v) nDaughters =
v->particles_out_size();
577 if (
p->status() == 3)
return false;
583 double m2 =
p->momentum().m2();
590 EvtId evtId = EvtPDL::evtIdFromStdHep(
id);
593 if (evtId.getId()>=0)
595 nModes = EvtDecayTable::getInstance()->getNMode(evtId.getAlias());
598 <<
" (status = " <<
p->status()
599 <<
") -- " << nModes <<
" decay modes found");
617 for (
auto itd:
v->particles_out()) {
618 if (std::abs(itd->pdg_id()) == std::abs(
id))
return false;
621 for (HepMC::GenVertex::particle_iterator itd =
v->particles_begin(
HepMC::children);
624 if (std::abs((*itd)->pdg_id()) == std::abs(
id))
return false;
648 int id = std::abs(pId);
671 std::vector<HepMC::GenParticlePtr> *muons =
new std::vector<HepMC::GenParticlePtr>;
673 for (
const auto&
p: *hepMC) {
674 if( std::abs(
p->pdg_id()) == 13 )
678 for (
auto muItr1 = muons->begin(); muItr1 != muons->end(); ++muItr1) {
679 for (
auto muItr2 = muItr1+1; muItr2 != muons->end(); ++muItr2) {
687 double dimuMass =
invMass((*muItr1),(*muItr2));
700 double p1Px = p1->momentum().px();
701 double p1Py = p1->momentum().py();
702 double p1Pz = p1->momentum().pz();
703 double p1E = p1->momentum().e();
704 double p2Px = p2->momentum().px();
705 double p2Py = p2->momentum().py();
706 double p2Pz = p2->momentum().pz();
707 double p2E = p2->momentum().e();
708 double dimuE = p2E + p1E;
709 double dimuPx = p2Px + p1Px;
710 double dimuPy = p2Py + p1Py;
711 double dimuPz = p2Pz + p1Pz;
712 double invMass = std::sqrt(dimuE*dimuE - dimuPx*dimuPx - dimuPy*dimuPy - dimuPz*dimuPz);
724 std::set<HepMC::GenVertexPtr> visited;
725 unsigned int nParticlesFound = 0;
726 unsigned int nTreesFound = 0;
727 for (
auto p: *hepMC) {
728 if ( (!
p->production_vertex()) ||
729 (
p->production_vertex()->particles_in().size() == 0) ) {
731 std::cout <<
"\n Found new partial decay tree:\n" << std::endl;
732 unsigned int nParticlesVisited =
printTree(
p,visited,1,barcodeList);
733 std::cout <<
"\n " << nParticlesVisited <<
" particles in this subtree" << std::endl;
734 nParticlesFound += nParticlesVisited;
737 std::cout <<
"\n Total of " << nParticlesFound <<
" particles found in "
738 << nTreesFound <<
" decay subtrees in HepMC event record\n" << std::endl;
742 std::set<HepMC::GenVertexPtr> visited;
743 unsigned int nParticlesFound = 0;
744 unsigned int nTreesFound = 0;
745 for (HepMC::GenEvent::particle_iterator
itp = hepMC->particles_begin();
itp != hepMC->particles_end(); ++
itp) {
747 if ( (!
p->production_vertex()) ||
748 (
p->production_vertex()->particles_in_size() == 0) ) {
750 std::cout <<
"\n Found new partial decay tree:\n" << std::endl;
751 unsigned int nParticlesVisited =
printTree(
p,visited,1,barcodeList);
752 std::cout <<
"\n " << nParticlesVisited <<
" particles in this subtree" << std::endl;
753 nParticlesFound += nParticlesVisited;
756 std::cout <<
"\n Total of " << nParticlesFound <<
" particles found in "
757 << nTreesFound <<
" decay subtrees in HepMC event record\n" << std::endl;
763 std::set<HepMC::GenVertexPtr>& visited,
int level, std::set<HepMC::GenParticlePtr>* barcodeList) {
764 unsigned int nParticlesVisited = 1;
765 for (
int i=0;
i<
level;
i++) std::cout <<
" ";
767 auto v =
p->end_vertex();
769 if (
v->particles_in().size() > 1)
770 std::cout <<
" [interaction: " <<
v->particles_in().size() <<
" particles, vertex " <<
v <<
"] --> ";
772 std::cout <<
" --> ";
773 if (visited.insert(
v).second) {
774 for (
auto itp:
v->particles_out()) {
777 std::cout << std::endl;
778 for (
auto itp:
v->particles_out()) {
779 if (
itp->end_vertex())
785 std::cout <<
"see above" << std::endl;
787 std::cout <<
" no decay vertex\n" << std::endl;
788 return nParticlesVisited;
792 std::set<HepMC::GenVertexPtr>& visited,
int level, std::set<int>* barcodeList) {
793 unsigned int nParticlesVisited = 1;
794 for (
int i=0;
i<
level;
i++) std::cout <<
" ";
796 auto v =
p->end_vertex();
799 if (
v->particles_in().size() > 1)
800 std::cout <<
" [interaction: " <<
v->particles_in().size() <<
" particles, vertex " <<
v <<
"] --> ";
802 std::cout <<
" --> ";
803 if (visited.insert(
v).second) {
804 for (
auto itp:
v->particles_out()) {
807 std::cout << std::endl;
808 for (
auto itp:
v->particles_out()) {
809 if (
itp->end_vertex())
815 std:: cout <<
"see above" << std::endl;
817 std::cout <<
" no decay vertex\n" << std::endl;
820 if (
v->particles_in_size() > 1)
821 std::cout <<
" [interaction: " <<
v->particles_in_size() <<
" particles, vertex " <<
v <<
"] --> ";
823 std::cout <<
" --> ";
824 if (visited.insert(
v).second) {
830 std::cout << std::endl;
834 if ((*itp)->end_vertex())
840 std:: cout <<
"see above" << std::endl;
842 std::cout <<
" no decay vertex\n" << std::endl;
844 return nParticlesVisited;
850 std::ostringstream buf;
852 if (barcodeList)
for (
const auto& pinl: *barcodeList)
if (pinl&&
p)
if (pinl.get()==
p.get()) inlist=
true;
853 if (statusHighlighting) {
854 if ( ((barcodeList!=0) && (inlist)) ||
866 if (statusHighlighting) {
873 std::ostringstream buf;
874 if (statusHighlighting) {
875 if ( ((barcodeList!=0) && (barcodeList->find(
HepMC::barcode(
p)) != barcodeList->end())) ||
887 if (statusHighlighting) {
904 return CLHEP::RandFlat::shoot(
m_engine);
909 static void local_split( std::vector <std::string>&
tokens,
const std::string&
input,
const char sep) {
911 for (
size_t i = 0;
i <=
input.size();
i++) {
922 char *cmtpath =
getenv(
"CMTPATH");
923 char *cmtconfig =
getenv(
"CMTCONFIG");
925 std::string foundpath =
"";
927 if(cmtpath != 0 && cmtconfig != 0){
929 std::vector<std::string> cmtpaths;
930 local_split(cmtpaths, cmtpath,
':');
932 std::string
installPath =
"/InstallArea/" + std::string(cmtconfig) +
"/share/Pythia8/xmldoc";
934 for(std::vector<std::string>::const_iterator
path = cmtpaths.begin();
935 path != cmtpaths.end() && foundpath ==
""; ++
path){
937 std::ifstream
testFile(testPath.c_str());
938 if(
testFile.good()) foundpath = testPath;