ATLAS Offline Software
Loading...
Searching...
No Matches
EvtInclusiveDecay Class Reference

#include <EvtInclusiveDecay.h>

Inheritance diagram for EvtInclusiveDecay:
Collaboration diagram for EvtInclusiveDecay:

Classes

struct  ParticleIdCompare

Public Member Functions

 EvtInclusiveDecay (const std::string &name, ISvcLocator *pSvcLocator)
virtual ~EvtInclusiveDecay ()
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()
std::string xmlpath (void)
virtual StatusCode sysInitialize () override
 Override sysInitialize.
virtual const DataObjIDColl & extraOutputDeps () const override
 Return the list of extra output dependencies.
ServiceHandle< StoreGateSvc > & evtStore ()
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode sysStart () override
 Handle START transition.
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles.
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T, V, H > &t)
void updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream & msg () const
bool msgLvl (const MSG::Level lvl) const
Event collection accessors (const and non-const)
HepMC::GenEvent *event ATLAS_NOT_CONST_THREAD_SAFE ()
 Access the current signal event (first in the McEventCollection)
McEventCollection *events ATLAS_NOT_CONST_THREAD_SAFE ()
 Access the current event's McEventCollection.
const HepMC::GenEvent * event_const () const
 Access the current signal event (const)
const McEventCollectionevents_const () const
 Access the current event's McEventCollection (const)
const McEventCollectionevents_const (const EventContext &ctx) const
Particle data accessors
const ServiceHandle< IPartPropSvc > partPropSvc () const
 Access the particle property service.
const HepPDT::ParticleDataTable & particleTable () const
 Get a particle data table.
const HepPDT::ParticleDataTable & pdt () const
 Shorter alias to get a particle data table.
const HepPDT::ParticleData * particleData (int pid) const
 Access an element in the particle data table.

Protected Member Functions

void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce (T &h)
void extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
 Add StoreName to extra input/output deps as needed.

Protected Attributes

Properties
std::string m_mcEventKey {}
 StoreGate key for the MC event collection (defaults to GEN_EVENT)
BooleanProperty m_mkMcEvent {this, "MakeMcEvent", false, "Create a new MC event collection if it doesn't exist"}
 Flag to determine if a new MC event collection should be made if it doesn't exist.

Private Types

typedef ServiceHandle< StoreGateSvcStoreGateSvc_t

Private Member Functions

Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey>

Private Attributes

DataObjIDColl m_extendedExtraObjects
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default)
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default)
std::vector< SG::VarHandleKeyArray * > m_vhka
bool m_varHandleArraysDeclared

Features for derived classes to use internally

ServiceHandle< IAthRNGSvcm_rndmSvc {this, "RndmSvc", "AthRNGSvc"}
IntegerProperty m_dsid {this, "Dsid", 999999}
IntegerProperty m_randomSeed {this, "RandomSeed", 1234567, "Random seed for the built-in random engine"}
 Seed for random number engine.
McEventCollectionm_mcEvtColl {}
EvtInclusiveAtRndmGenm_evtAtRndmGen {}
EvtGen * m_myEvtGen {}
std::string m_pdtFile
std::string m_decayFile
std::string m_userDecayFile
std::string m_randomStreamName
std::string m_inputKeyName
std::string m_outputKeyName
bool m_readExisting
bool m_prohibitFinalStateDecay
bool m_prohibitReDecay
bool m_prohibitUnDecay
bool m_prohibitRemoveSelfDecay
std::vector< int > m_blackList
std::set< int > m_blackListSet
bool m_allowAllKnownDecays
bool m_allowDefaultBDecays
std::vector< int > m_whiteList
std::set< int > m_whiteListSet
bool m_printHepMCBeforeEvtGen
bool m_printHepMCAfterEvtGen
bool m_printHepMCHighlighted
bool m_printHepMCHighLightTopLevelDecays
bool m_checkDecayTree
bool m_checkDecayChannels
std::map< int, long > m_noDecayChannels
int m_nRepeatedDecays
int m_maxNRepeatedDecays
bool m_applyUserSelection
bool m_userSelRequireOppositeSignedMu
double m_userSelMu1MinPt
double m_userSelMu2MinPt
double m_userSelMu1MaxEta
double m_userSelMu2MaxEta
double m_userSelMinDimuMass
double m_userSelMaxDimuMass
bool m_isfHerwig
bool m_setVMtransversePol
void reseedRandomEngine (const std::string &streamName, const EventContext &ctx)
CLHEP::HepRandomEngine * getRandomEngine (const std::string &streamName, unsigned long int randomSeedOffset, const EventContext &ctx) const
CLHEP::HepRandomEngine * getRandomEngineDuringInitialize (const std::string &streamName, unsigned long int randomSeedOffset, unsigned int conditionsRun=1, unsigned int lbn=1) const
StatusCode traverseDecayTree (HepMC::GenParticlePtr p, bool isToBeRemoved, std::set< HepMC::GenVertexPtr > &visited, std::set< HepMC::GenParticlePtr, ParticleIdCompare > &toBeDecayed)
void removeDecayTree (HepMC::GenEvent *hepMC, HepMC::GenParticlePtr p)
void decayParticle (HepMC::GenEvent *hepMC, HepMC::GenParticlePtr p)
void addEvtGenDecayTree (HepMC::GenEvent *hepMC, HepMC::GenParticlePtr part, EvtParticle *evtPart, EvtVector4R treeStart, double momentumScaleFactor=1.0)
bool isToBeDecayed (HepMC::ConstGenParticlePtr p, bool doCrossChecks)
bool isDefaultB (const int pId) const
bool passesUserSelection (HepMC::GenEvent *hepMC)
double invMass (HepMC::ConstGenParticlePtr p1, HepMC::ConstGenParticlePtr p2)
void printHepMC (HepMC::GenEvent *hepMC, std::set< HepMC::GenParticlePtr, ParticleIdCompare > *particleSet=nullptr)
unsigned int printTree (HepMC::GenParticlePtr p, std::set< HepMC::GenVertexPtr > &visited, int level, std::set< HepMC::GenParticlePtr, ParticleIdCompare > *particleSet=nullptr)
std::string pdgName (HepMC::ConstGenParticlePtr p, bool statusHighlighting=false, std::set< HepMC::GenParticlePtr, ParticleIdCompare > *particleSet=nullptr)

Utility event-mangling functions

Todo
Replace with HepMC units when available
ServiceHandle< IPartPropSvc > m_ppSvc {this, "PartPropSvc", "PartPropSvc"}
 Handle on the particle property service.
SG::ReadHandleKey< McEventCollectionm_mcevents_const { this, "McEventKey", "GEN_EVENT", "StoreGate key of the MC event collection" }
 Const handle to the MC event collection.
void GeVToMeV (HepMC::GenEvent *evt)
 Scale event energies/momenta by x 1000.
void MeVToGeV (HepMC::GenEvent *evt)
 Scale event energies/momenta by x 1/1000.
void cmTomm (HepMC::GenEvent *evt)
 Scale event lengths by x 10.
void mmTocm (HepMC::GenEvent *evt)
 Scale event lengths by x 1/10.

Detailed Description

Definition at line 55 of file EvtInclusiveDecay.h.

Member Typedef Documentation

◆ StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< Algorithm > >::StoreGateSvc_t
privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

◆ EvtInclusiveDecay()

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

Definition at line 57 of file EvtInclusiveDecay.cxx.

57 :
58 GenBase( name, pSvcLocator ),
60
61 // Basic EvtGen configuration: decay and particle definition files, random number stream
62 declareProperty("pdtFile", m_pdtFile = "inclusive.pdt");
63 declareProperty("decayFile", m_decayFile = "2014inclusive.dec");
64 declareProperty("userDecayFile", m_userDecayFile = "");
65 declareProperty("randomStreamName", m_randomStreamName = "EVTGEN");
66 declareProperty("inputKeyName", m_inputKeyName = "GEN_EVENT");
67 declareProperty("outputKeyName",m_outputKeyName = "GEN_EVENT_EVTGEN");
68 declareProperty("readExisting",m_readExisting=false);
69
70 // Selection of particles to be decayed
71 declareProperty("prohibitFinalStateDecay", m_prohibitFinalStateDecay=false);
72 declareProperty("prohibitReDecay", m_prohibitReDecay=false);
73 declareProperty("prohibitUnDecay", m_prohibitUnDecay=true);
74 declareProperty("prohibitRemoveSelfDecay", m_prohibitRemoveSelfDecay=false);
75 declareProperty("blackList",m_blackList);
76 declareProperty("allowAllKnownDecays", m_allowAllKnownDecays=true);
77 declareProperty("allowDefaultBDecays", m_allowDefaultBDecays=true);
78 declareProperty("whiteList",m_whiteList);
79
80 // Level of output
81 declareProperty("printHepMCBeforeEvtGen", m_printHepMCBeforeEvtGen=false);
82 declareProperty("printHepMCAfterEvtGen", m_printHepMCAfterEvtGen=false);
83 declareProperty("printHepMCHighlighted", m_printHepMCHighlighted=true);
84 declareProperty("printHepMCHighLightTopLevelDecays", m_printHepMCHighLightTopLevelDecays=true);
85
86 // Optional checks
87 declareProperty("checkDecayTree", m_checkDecayTree=false);
88 declareProperty("checkDecayChannels", m_checkDecayChannels=false);
89
90 // Repeated decays
91 declareProperty("maxNRepeatedDecays", m_maxNRepeatedDecays=1);
92
93 // User selection
94 declareProperty("applyUserSelection", m_applyUserSelection=false);
95 declareProperty("userSelRequireOppositeSignedMu", m_userSelRequireOppositeSignedMu=true);
96 declareProperty("userSelMu1MinPt", m_userSelMu1MinPt=0.);
97 declareProperty("userSelMu2MinPt", m_userSelMu2MinPt=0.);
98 declareProperty("userSelMu1MaxEta", m_userSelMu1MaxEta=102.5);
99 declareProperty("userSelMu2MaxEta", m_userSelMu2MaxEta=102.5);
100 declareProperty("userSelMinDimuMass", m_userSelMinDimuMass=0.);
101 declareProperty("userSelMaxDimuMass", m_userSelMaxDimuMass=-1.); // set to negative to not apply cut
102 declareProperty("isfHerwig", m_isfHerwig=false);
103 declareProperty("setVMtransversePol", m_setVMtransversePol=false);
104
105 // We have decided to blacklist Tau decays because we are not sure whether the polarization
106 // would be properly passed to EvtGen
107 m_blackList.push_back(15);
108}
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
std::string m_randomStreamName
std::string m_userDecayFile
std::vector< int > m_blackList
std::string m_outputKeyName
std::vector< int > m_whiteList
GenBase(const std::string &name, ISvcLocator *pSvcLocator)
Constructor.
Definition GenBase.cxx:11

◆ ~EvtInclusiveDecay()

EvtInclusiveDecay::~EvtInclusiveDecay ( )
virtual

Definition at line 112 of file EvtInclusiveDecay.cxx.

112 {
113 delete m_myEvtGen;
114 delete m_evtAtRndmGen;
115}
EvtInclusiveAtRndmGen * m_evtAtRndmGen

Member Function Documentation

◆ addEvtGenDecayTree()

void EvtInclusiveDecay::addEvtGenDecayTree ( HepMC::GenEvent * hepMC,
HepMC::GenParticlePtr part,
EvtParticle * evtPart,
EvtVector4R treeStart,
double momentumScaleFactor = 1.0 )
private

Definition at line 477 of file EvtInclusiveDecay.cxx.

478 {
479 if(evtPart->getNDaug()!=0) {
480 // Add decay vertex, starting from production vertex of particle
481 double ct=(evtPart->getDaug(0)->get4Pos()).get(0)+treeStart.get(0);
482 double x=(evtPart->getDaug(0)->get4Pos()).get(1)+treeStart.get(1);
483 double y=(evtPart->getDaug(0)->get4Pos()).get(2)+treeStart.get(2);
484 double z=(evtPart->getDaug(0)->get4Pos()).get(3)+treeStart.get(3);
485
486 HepMC::GenVertexPtr end_vtx = HepMC::newGenVertexPtr(HepMC::FourVector(x,y,z,ct));
487
488 hepMC->add_vertex(end_vtx);
489 end_vtx->add_particle_in(std::move(part));
490
491 // Add decay daughter with their own decay trees
492 for(uint it=0; it<evtPart->getNDaug(); it++) {
493 double e=(evtPart->getDaug(it)->getP4Lab()).get(0) * momentumScaleFactor;
494 double px=(evtPart->getDaug(it)->getP4Lab()).get(1) * momentumScaleFactor;
495 double py=(evtPart->getDaug(it)->getP4Lab()).get(2) * momentumScaleFactor;
496 double pz=(evtPart->getDaug(it)->getP4Lab()).get(3) * momentumScaleFactor;
497 int id=EvtPDL::getStdHep(evtPart->getDaug(it)->getId());
498 int status=1;
499 if(evtPart->getDaug(it)->getNDaug() != 0) status=2;
500 HepMC::GenParticlePtr daughter = HepMC::newGenParticlePtr(HepMC::FourVector(px,py,pz,e),id,status);
501 end_vtx->add_particle_out(daughter);
502 addEvtGenDecayTree(hepMC, std::move(daughter), evtPart->getDaug(it), treeStart, momentumScaleFactor);
503 }
504 }
505}
unsigned int uint
#define y
#define x
#define z
void addEvtGenDecayTree(HepMC::GenEvent *hepMC, HepMC::GenParticlePtr part, EvtParticle *evtPart, EvtVector4R treeStart, double momentumScaleFactor=1.0)
HepMC::GenVertex * GenVertexPtr
Definition GenVertex.h:59
GenVertexPtr newGenVertexPtr(const HepMC::FourVector &pos=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), const int i=0)
Definition GenVertex.h:64
GenParticlePtr newGenParticlePtr(const HepMC::FourVector &mom=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), int pid=0, int status=0)
Definition GenParticle.h:39
GenParticle * GenParticlePtr
Definition GenParticle.h:37
status
Definition merge.py:16

◆ ATLAS_NOT_CONST_THREAD_SAFE() [1/2]

HepMC::GenEvent *event GenBase::ATLAS_NOT_CONST_THREAD_SAFE ( )
inlineinherited

Access the current signal event (first in the McEventCollection)

Note
This function will make a new McEventCollection if there is not already a valid one and MakeMcEvent=True.

Definition at line 76 of file GenBase.h.

76 {
77 if (events()->empty())
78 ATH_MSG_ERROR("McEventCollection is empty during first event access");
79 return *(events()->begin());
80 }
#define ATH_MSG_ERROR(x)
static const Attributes_t empty

◆ ATLAS_NOT_CONST_THREAD_SAFE() [2/2]

McEventCollection *events GenBase::ATLAS_NOT_CONST_THREAD_SAFE ( )
inherited

Access the current event's McEventCollection.

Note
This function will make a new McEventCollection if there is not already a valid one and MakeMcEvent=True.

◆ cmTomm()

void GenBase::cmTomm ( HepMC::GenEvent * evt)
protectedinherited

Scale event lengths by x 10.

Definition at line 78 of file GenBase.cxx.

78 {
79 for (HepMC::GenEvent::vertex_iterator vtx = evt->vertices_begin(); vtx != evt->vertices_end(); ++vtx) {
80 const HepMC::FourVector fv((*vtx)->position().x() * 10,
81 (*vtx)->position().y() * 10,
82 (*vtx)->position().z() * 10,
83 (*vtx)->position().t() * 10);
84 (*vtx)->set_position(fv);
85 }
86}

◆ decayParticle()

void EvtInclusiveDecay::decayParticle ( HepMC::GenEvent * hepMC,
HepMC::GenParticlePtr p )
private

Definition at line 439 of file EvtInclusiveDecay.cxx.

439 {
440 ATH_MSG_DEBUG("Decaying particle " << pdgName(part) << " " << part);
441 if (msgLvl(MSG::VERBOSE)) HepMC::Print::line(std::cout,part);
442
443 // Remove existing decay tree, if any, and flag particle as being decayed by EvtGen
444 removeDecayTree(hepMC,part);
445 part->set_status(0);
446
447 // Create EvtGen version of part and have EvtGen decay it.
448 // Since EvtGen uses GeV, convert particles momentum from MeV to GeV.
449 int id = part->pdg_id();
450 EvtId evtId=EvtPDL::evtIdFromStdHep(id);
451 double en =(part->momentum()).e()/1000.;
452 double px=(part->momentum()).px()/1000.;
453 double py=(part->momentum()).py()/1000.;
454 double pz=(part->momentum()).pz()/1000.;
455 EvtVector4R evtP(en,px,py,pz);
456 EvtParticle* evtPart = EvtParticleFactory::particleFactory(evtId,evtP);
457
458 // set transverse polarization to vector mesons (relevant for coherent production of J/Psi etc in UPC)
459 if(m_setVMtransversePol && (id==113 || id== 443 || id==100443 || id==553 || id==100553 || id==200553) )evtPart->setVectorSpinDensity();
460
461 m_myEvtGen->generateDecay(evtPart);
462 if (msgLvl(MSG::VERBOSE)) evtPart->printTree();
463 double ct_s = part->production_vertex()->position().t();
464 double x_s = part->production_vertex()->position().x();
465 double y_s = part->production_vertex()->position().y();
466 double z_s = part->production_vertex()->position().z();
467
468 EvtVector4R treeStart(ct_s,x_s,y_s,z_s);
469 // Add new decay tree to hepMC, converting back from GeV to MeV.
470 addEvtGenDecayTree(hepMC, part, evtPart, treeStart, 1000.);
471 if(evtPart->getNDaug() !=0) part->set_status(2);
472 evtPart->deleteTree();
473}
#define ATH_MSG_DEBUG(x)
bool msgLvl(const MSG::Level lvl) const
void removeDecayTree(HepMC::GenEvent *hepMC, HepMC::GenParticlePtr p)
std::string pdgName(HepMC::ConstGenParticlePtr p, bool statusHighlighting=false, std::set< HepMC::GenParticlePtr, ParticleIdCompare > *particleSet=nullptr)
void line(std::ostream &os, const GenEvent &e)
Definition GenEvent.h:677

◆ declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< Algorithm > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl,
const SG::VarHandleKeyType &  )
inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

158 {
160 hndl.value(),
161 hndl.documentation());
162
163 }

◆ declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty ( Gaudi::Property< T, V, H > & t)
inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

145 {
146 typedef typename SG::HandleClassifier<T>::type htype;
148 }
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>

◆ detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< Algorithm > >::detStore ( ) const
inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

◆ event_const()

const HepMC::GenEvent * GenBase::event_const ( ) const
inlineinherited

Access the current signal event (const)

Definition at line 83 of file GenBase.h.

83 {
84 if (events_const()->empty())
85 ATH_MSG_ERROR("Const McEventCollection is empty during first event access");
86 return *(events_const()->begin());
87 }
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
const McEventCollection * events_const() const
Access the current event's McEventCollection (const)
Definition GenBase.h:96

◆ events_const() [1/2]

const McEventCollection * GenBase::events_const ( ) const
inlineinherited

Access the current event's McEventCollection (const)

Definition at line 96 of file GenBase.h.

96 {
97 return events_const( getContext() );
98 }

◆ events_const() [2/2]

const McEventCollection * GenBase::events_const ( const EventContext & ctx) const
inlineinherited

Definition at line 99 of file GenBase.h.

99 {
100 SG::ReadHandle<McEventCollection> ret = SG::makeHandle(m_mcevents_const, ctx);
101 if (!ret.isValid())
102 ATH_MSG_ERROR("No McEventCollection found in StoreGate with key " << m_mcevents_const.key());
103 return ret.cptr();
104 }
SG::ReadHandleKey< McEventCollection > m_mcevents_const
Const handle to the MC event collection.
Definition GenBase.h:163
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())

◆ evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< Algorithm > >::evtStore ( )
inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

◆ execute()

StatusCode EvtInclusiveDecay::execute ( )
virtual

Reimplemented from GenBase.

Definition at line 218 of file EvtInclusiveDecay.cxx.

218 {
219 ATH_MSG_DEBUG("EvtInclusiveDecay executing");
220
221 const EventContext& ctx = Gaudi::Hive::currentContext();
223
224 std::string key = m_inputKeyName;
225 // retrieve event from Transient Store (Storegate)
226
227 // Load HepMC info
228 // FIXME should be using Read/WriteHandles here
229 const McEventCollection* oldmcEvtColl{};
230 if(m_readExisting) {
231 CHECK(evtStore()->retrieve(oldmcEvtColl, key));
232 // Fill the new McEventCollection with a copy of the initial HepMC::GenEvent
233 m_mcEvtColl = new McEventCollection(*oldmcEvtColl);
234 }
235 else {CHECK(evtStore()->retrieve(m_mcEvtColl, key));}
236
237 if(m_readExisting) {
238 if(m_outputKeyName!=key) {
240 }
241 }
242
244 for( mcItr = m_mcEvtColl->begin(); mcItr != m_mcEvtColl->end(); ++mcItr ) {
245 HepMC::GenEvent* hepMC = *mcItr;
246
247 // Search HepMC record for particles to be decayed by EvtGen
248 // NOTE: In order to ensure repeatability, we use customized a std::set to obtain
249 // an ordered list of particles to be decayed by EvtGen.
250 std::set<HepMC::GenVertexPtr> visited;
251 std::set<HepMC::GenParticlePtr,ParticleIdCompare> toBeDecayed;
252 for (auto p: *hepMC) {
253 if ( (!p->production_vertex()) ||
254 (p->production_vertex()->particles_in_size() == 0) ) {
255 StatusCode sc = traverseDecayTree(std::move(p),false,visited,toBeDecayed);
256 if (sc.isFailure())
257 return StatusCode::FAILURE;
258 }
259 }
260 // Print HepMC in tree format if desired (before doing anything)
262 msg(MSG::INFO) << "Printing HepMC record at " << hepMC << " BEFORE running EvtGen:" << endmsg;
264 printHepMC(hepMC,&toBeDecayed);
265 else
266 printHepMC(hepMC);
267 }
268
269
270 // Decay selected particles
271 bool eventPassesCuts(false);
272 int loopCounter(0);
273 while( !eventPassesCuts && loopCounter < m_maxNRepeatedDecays ) {
274
275 for (auto p: toBeDecayed) {
276 if (p == 0) {
277#ifdef HEPMC3
278 msg(MSG::ERROR ) << "Overlapping decay tree for particle" << p <<endmsg;
279#else
280 msg(MSG::ERROR ) << "Overlapping decay tree encountered for barcode " << HepMC::barcode(p) << endmsg;
281#endif
282 return StatusCode::FAILURE;
283 }
284 decayParticle(hepMC,std::move(p));
286 }
287
289 eventPassesCuts = passesUserSelection(hepMC);
290 else
291 eventPassesCuts = true;
292
294 loopCounter++;
295 }
296
297 // Store the number of decay attempts in event weights std::map, only if repeated decays enabled
298
299 if(m_maxNRepeatedDecays > 1) {
300#ifdef HEPMC3
301 hepMC->weight("nEvtGenDecayAttempts") = loopCounter;
302#else
303 hepMC->weights()["nEvtGenDecayAttempts"] = loopCounter;
304#endif
305 }
306 // Print HepMC in tree format if desired (after finishing all EvtGen decays)
308 msg(MSG::INFO) << "Printing HepMC record at " << hepMC << " AFTER running EvtGen:" << endmsg;
310 printHepMC(hepMC,&toBeDecayed);
311 else
312 printHepMC(hepMC);
313 }
314 }
315
316 if(m_readExisting && m_outputKeyName==key) {
318 }
319
320 return StatusCode::SUCCESS;
321}
#define endmsg
#define CHECK(...)
Evaluate an expression and check for errors.
static Double_t sc
DataModel_detail::iterator< DataVector > iterator
Definition DataVector.h:842
void printHepMC(HepMC::GenEvent *hepMC, std::set< HepMC::GenParticlePtr, ParticleIdCompare > *particleSet=nullptr)
void reseedRandomEngine(const std::string &streamName, const EventContext &ctx)
StatusCode traverseDecayTree(HepMC::GenParticlePtr p, bool isToBeRemoved, std::set< HepMC::GenVertexPtr > &visited, std::set< HepMC::GenParticlePtr, ParticleIdCompare > &toBeDecayed)
bool passesUserSelection(HepMC::GenEvent *hepMC)
McEventCollection * m_mcEvtColl
void decayParticle(HepMC::GenEvent *hepMC, HepMC::GenParticlePtr p)
::StatusCode StatusCode
StatusCode definition for legacy code.
int barcode(const T *p)
Definition Barcode.h:16
void fillBarcodesAttribute(GenEvent *)
Definition GenEvent.h:627
retrieve(aClass, aKey=None)
Definition PyKernel.py:110
MsgStream & msg
Definition testRead.cxx:32

◆ extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< Algorithm > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps)
protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

◆ extraOutputDeps()

const DataObjIDColl & AthAlgorithm::extraOutputDeps ( ) const
overridevirtualinherited

Return the list of extra output dependencies.

This list is extended to include symlinks implied by inheritance relations.

Definition at line 50 of file AthAlgorithm.cxx.

51{
52 // If we didn't find any symlinks to add, just return the collection
53 // from the base class. Otherwise, return the extended collection.
54 if (!m_extendedExtraObjects.empty()) {
56 }
57 return Algorithm::extraOutputDeps();
58}
DataObjIDColl m_extendedExtraObjects

◆ finalize()

StatusCode EvtInclusiveDecay::finalize ( )

Definition at line 324 of file EvtInclusiveDecay.cxx.

324 {
325
327 ATH_MSG_INFO("The following particles were checked and didn't have any decay channels:");
328 if (msgLvl(MSG::INFO)) {
329 std::cout << std::endl;
330 std::cout << " Particle code Name from HepPDT # Occurences" << std::endl;
331 std::cout << "------------------------------------------------------" << std::endl;
332 for (std::map<int,long>::iterator p = m_noDecayChannels.begin(); p!=m_noDecayChannels.end(); ++p) {
333 int id = p->first;
334 int count = p->second;
335 std::cout << std::setw(14) << id
336 << std::setw(20) << HepPID::particleName(id)
337 << std::setw(20) << count
338 << std::endl;
339 }
340 std::cout << std::endl;
341 }
342 }
343 ATH_MSG_INFO("Total number of repeated decays: " << m_nRepeatedDecays);
344 ATH_MSG_INFO("EvtInclusiveDecay finalized");
345 return StatusCode::SUCCESS;
346}
#define ATH_MSG_INFO(x)
std::map< int, long > m_noDecayChannels
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:146

◆ getRandomEngine()

CLHEP::HepRandomEngine * EvtInclusiveDecay::getRandomEngine ( const std::string & streamName,
unsigned long int randomSeedOffset,
const EventContext & ctx ) const
private

Definition at line 192 of file EvtInclusiveDecay.cxx.

194{
195 ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this, streamName);
196 rngWrapper->setSeed( streamName, ctx.slot(), randomSeedOffset, ctx.eventID().run_number() );
197 return rngWrapper->getEngine(ctx);
198}
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
Definition RNGWrapper.h:169
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
Definition RNGWrapper.h:134
ServiceHandle< IAthRNGSvc > m_rndmSvc

◆ getRandomEngineDuringInitialize()

CLHEP::HepRandomEngine * EvtInclusiveDecay::getRandomEngineDuringInitialize ( const std::string & streamName,
unsigned long int randomSeedOffset,
unsigned int conditionsRun = 1,
unsigned int lbn = 1 ) const
private

Definition at line 201 of file EvtInclusiveDecay.cxx.

202{
203 const size_t slot=0;
204 EventContext ctx;
205 ctx.setSlot( slot );
206 ctx.setEventID (EventIDBase (conditionsRun,
207 EventIDBase::UNDEFEVT, // event
208 EventIDBase::UNDEFNUM, // timestamp
209 EventIDBase::UNDEFNUM, // timestamp ns
210 lbn));
212 Atlas::ExtendedEventContext( evtStore()->hiveProxyDict(),
213 conditionsRun) );
214 return getRandomEngine(streamName, randomSeedOffset, ctx);
215}
CLHEP::HepRandomEngine * getRandomEngine(const std::string &streamName, unsigned long int randomSeedOffset, const EventContext &ctx) const
void setExtendedEventContext(EventContext &ctx, ExtendedEventContext &&ectx)
Move an extended context into a context object.

◆ GeVToMeV()

void GenBase::GeVToMeV ( HepMC::GenEvent * evt)
protectedinherited

Scale event energies/momenta by x 1000.

Todo
Add HepMC units awareness and do it differently when HepMC provides this functionality directly (and reference-based FourVector accessors)

Definition at line 58 of file GenBase.cxx.

58 {
59 for (HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
60 const HepMC::FourVector fv((*p)->momentum().px() * 1000,
61 (*p)->momentum().py() * 1000,
62 (*p)->momentum().pz() * 1000,
63 (*p)->momentum().e() * 1000);
64 (*p)->set_momentum(fv);
65 (*p)->set_generated_mass(1000 * (*p)->generated_mass());
66 }
67}

◆ initialize()

StatusCode EvtInclusiveDecay::initialize ( )
virtual

Reimplemented from GenBase.

Definition at line 119 of file EvtInclusiveDecay.cxx.

119 {
120
122 // Get the random number service
123 CHECK(m_rndmSvc.retrieve());
124
125 msg(MSG::INFO) << "EvtInclusiveDecay initialize" << endmsg;
126 msg(MSG::INFO) << "Particle properties definition file = " << m_pdtFile << endmsg;
127 msg(MSG::INFO) << "Main decay file = " << m_decayFile << endmsg;
128 msg(MSG::INFO) << "User decay file = " << m_userDecayFile << endmsg;
129 msg(MSG::INFO) << "Max number of repeated decays = " << m_maxNRepeatedDecays << endmsg;
130 msg(MSG::INFO) << "EvtInclusiveDecay selection parameters:" << endmsg;
131 msg(MSG::INFO) << "* prohibitFinalStateDecay = " << m_prohibitFinalStateDecay << endmsg;
132 msg(MSG::INFO) << "* prohibitReDecay = " << m_prohibitReDecay << endmsg;
133 msg(MSG::INFO) << "* prohibitUnDecay = " << m_prohibitUnDecay << endmsg;
134 msg(MSG::INFO) << "* prohibitRemoveSelfDecay = " << m_prohibitRemoveSelfDecay << endmsg;
135 msg(MSG::INFO) << "* allowAllKnownDecays = " << m_allowAllKnownDecays << endmsg;
136 msg(MSG::INFO) << "* allowDefaultBDecays = " << m_allowDefaultBDecays << endmsg;
137 msg(MSG::INFO) << "User selection parameters:" << endmsg;
138 msg(MSG::INFO) << "* applyUserSelection = " << m_applyUserSelection << endmsg;
139 msg(MSG::INFO) << "* userSelRequireOppositeSignedMu = " << m_userSelRequireOppositeSignedMu << endmsg;
140 msg(MSG::INFO) << "* userSelMu1MinPt = " << m_userSelMu1MinPt << endmsg;
141 msg(MSG::INFO) << "* userSelMu2MinPt = " << m_userSelMu2MinPt << endmsg;
142 msg(MSG::INFO) << "* userSelMu1MaxEta = " << m_userSelMu1MaxEta << endmsg;
143 msg(MSG::INFO) << "* userSelMu2MaxEta = " << m_userSelMu2MaxEta << endmsg;
144 msg(MSG::INFO) << "* userSelMinDimuMass = " << m_userSelMinDimuMass << endmsg;
145 msg(MSG::INFO) << "* userSelMaxDimuMass = " << m_userSelMaxDimuMass << endmsg;
146
147 // Initialize and print blackList
148 m_blackListSet.insert(m_blackList.begin(),m_blackList.end());
149 msg(MSG::INFO) << "* blackList; = ";
150 for (std::set<int>::iterator i = m_blackListSet.begin(); i!=m_blackListSet.end(); ++i)
151 msg(MSG::INFO) << (*i) << " ";
152 msg(MSG::INFO)<< endmsg;
153
154 // Initialize and print whiteList
155 m_whiteListSet.insert(m_whiteList.begin(),m_whiteList.end());
156 msg(MSG::INFO) << "* whiteList = ";
157 for (std::set<int>::iterator i = m_whiteListSet.begin(); i!=m_whiteListSet.end(); ++i)
158 msg(MSG::INFO) << (*i) << " ";
159 msg(MSG::INFO) << endmsg;
160
161 CLHEP::HepRandomEngine* rndmEngine = getRandomEngineDuringInitialize(m_randomStreamName, m_randomSeed, m_dsid);
162 // Obtain random number generator for EvtGen
163 m_evtAtRndmGen = new EvtInclusiveAtRndmGen(rndmEngine);
164
165 // Create an instance of EvtGen and read particle properties and decay files
166 EvtExternalGenList genList(true,xmlpath(),"gamma");
167 EvtAbsRadCorr* radCorrEngine = genList.getPhotosModel();
168 std::list<EvtDecayBase*> extraModels = genList.getListOfModels();
169
170 // Create the EvtGen generator object
171 // EvtGen myGenerator("decayFile.dec", "evt.pdl", randomEnginePointer,
172 // radCorrEngine, &extraModels);
173
174
175 m_myEvtGen = new EvtGen( m_decayFile.c_str(), m_pdtFile.c_str(), m_evtAtRndmGen, radCorrEngine, &extraModels);
176 if(!m_userDecayFile.empty())
177 m_myEvtGen->readUDecay(m_userDecayFile.c_str());
178
179 return StatusCode::SUCCESS;
180}
#define ATH_CHECK
Evaluate an expression and check for errors.
std::set< int > m_blackListSet
CLHEP::HepRandomEngine * getRandomEngineDuringInitialize(const std::string &streamName, unsigned long int randomSeedOffset, unsigned int conditionsRun=1, unsigned int lbn=1) const
IntegerProperty m_dsid
std::string xmlpath(void)
std::set< int > m_whiteListSet
IntegerProperty m_randomSeed
Seed for random number engine.
virtual StatusCode initialize() override
Definition GenBase.cxx:17

◆ inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< Algorithm > >::inputHandles ( ) const
overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ invMass()

double EvtInclusiveDecay::invMass ( HepMC::ConstGenParticlePtr p1,
HepMC::ConstGenParticlePtr p2 )
private

Definition at line 636 of file EvtInclusiveDecay.cxx.

636 {
637 double p1Px = p1->momentum().px();
638 double p1Py = p1->momentum().py();
639 double p1Pz = p1->momentum().pz();
640 double p1E = p1->momentum().e();
641 double p2Px = p2->momentum().px();
642 double p2Py = p2->momentum().py();
643 double p2Pz = p2->momentum().pz();
644 double p2E = p2->momentum().e();
645 double dimuE = p2E + p1E;
646 double dimuPx = p2Px + p1Px;
647 double dimuPy = p2Py + p1Py;
648 double dimuPz = p2Pz + p1Pz;
649 double invMass = std::sqrt(dimuE*dimuE - dimuPx*dimuPx - dimuPy*dimuPy - dimuPz*dimuPz);
650
651 return invMass;
652}
double invMass(HepMC::ConstGenParticlePtr p1, HepMC::ConstGenParticlePtr p2)

◆ isDefaultB()

bool EvtInclusiveDecay::isDefaultB ( const int pId) const
private

Definition at line 584 of file EvtInclusiveDecay.cxx.

584 {
585 int id = std::abs(pId);
586 if ( id == 511 ||
587 id == 521 ||
588 id == 531 ||
589 id == 541 ||
590 id == 5122 ||
591 id == 5132 ||
592 id == 5232 ||
593 id == 5112 ||
594 id == 5212 ||
595 id == 5222 )
596 return true;
597 else
598 return false;
599}

◆ isToBeDecayed()

bool EvtInclusiveDecay::isToBeDecayed ( HepMC::ConstGenParticlePtr p,
bool doCrossChecks )
private

Definition at line 515 of file EvtInclusiveDecay.cxx.

515 {
516 int id = p->pdg_id();
517 int nDaughters = 0;
518 auto v = p->end_vertex();
519 if (v) nDaughters = v->particles_out_size();
520
521 // Ignore documentation lines
522 if (p->status() == 3) return false;
523 // And any particles that aren't stable or decayed
524 if(!m_isfHerwig && !MC::isPhysical(p)) return false;
525
526 // Particularly for Herwig, try to ignore particles that really should
527 // be flagged as documentation lines
528 double m2 = p->momentum().m2();
529 if (m2 < -1.0E-3) {
530 ATH_MSG_DEBUG("Ignoring particle " << pdgName(std::move(p)) << " with m^2 = " << m2);
531 return false;
532 }
533
534 // Check whether EvtGen has any decay channels defined for this particle
535 EvtId evtId = EvtPDL::evtIdFromStdHep(id);
536 // std::cout << "EVTID: " << evtId.getId() << " alias " << evtId.getAlias() << std::endl;
537 int nModes = 0;
538 if (evtId.getId()>=0)
539 // nModes = EvtDecayTable::getNMode(evtId.getAlias());
540 nModes = EvtDecayTable::getInstance()->getNMode(evtId.getAlias());
541 if (doCrossChecks) {
542 ATH_MSG_VERBOSE("Checking particle " << pdgName(p)
543 << " (status = " << p->status()
544 <<") -- " << nModes << " decay modes found");
545 if (m_checkDecayChannels && nModes==0) {
546 std::map<int,long>::iterator pos = m_noDecayChannels.find(id);
547 if (pos != m_noDecayChannels.end())
548 (pos->second)++;
549 else
551 }
552 }
553
554 // Check prohibit* settings
555 if (m_prohibitFinalStateDecay && MC::isStable(p)) return false;
556 if (m_prohibitReDecay && nDaughters>0) return false;
557 if (m_prohibitUnDecay && nModes==0) return false;
558 if (m_prohibitRemoveSelfDecay && nDaughters>0) {
559 // For now, check only children - this should be sufficient and checking all
560 // descendants would be very expensive.
561 for (auto itd: *v) {
562 if (std::abs(itd->pdg_id()) == std::abs(id)) return false;
563 }
564 }
565
566 // Check blackList
567 if (m_blackListSet.count(std::abs(id))>0) return false;
568
569 // Check allow* settings
570 if (m_allowAllKnownDecays && nModes>0) return true;
571 if (m_allowDefaultBDecays && isDefaultB(id)) return true;
572
573 // Check whiteList
574 if (m_whiteListSet.count(std::abs(id))>0) return true;
575
576 return false; // Default is NOT to decay through EvtGen
577}
#define ATH_MSG_VERBOSE(x)
bool isDefaultB(const int pId) const
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
bool isPhysical(const T &p)
Identify if the particle is physical, i.e. is stable or decayed.

◆ MeVToGeV()

void GenBase::MeVToGeV ( HepMC::GenEvent * evt)
protectedinherited

Scale event energies/momenta by x 1/1000.

Definition at line 68 of file GenBase.cxx.

68 {
69 for (HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
70 const HepMC::FourVector fv((*p)->momentum().px() / 1000,
71 (*p)->momentum().py() / 1000,
72 (*p)->momentum().pz() / 1000,
73 (*p)->momentum().e() / 1000);
74 (*p)->set_momentum(fv);
75 (*p)->set_generated_mass((*p)->generated_mass() / 1000);
76 }
77}

◆ mmTocm()

void GenBase::mmTocm ( HepMC::GenEvent * evt)
protectedinherited

Scale event lengths by x 1/10.

Definition at line 87 of file GenBase.cxx.

87 {
88 for (HepMC::GenEvent::vertex_iterator vtx = evt->vertices_begin(); vtx != evt->vertices_end(); ++vtx) {
89 const HepMC::FourVector fv((*vtx)->position().x() / 10,
90 (*vtx)->position().y() / 10,
91 (*vtx)->position().z() / 10,
92 (*vtx)->position().t() / 10);
93 (*vtx)->set_position(fv);
94 }
95}

◆ msg()

MsgStream & AthCommonMsg< Algorithm >::msg ( ) const
inlineinherited

Definition at line 24 of file AthCommonMsg.h.

24 {
25 return this->msgStream();
26 }

◆ msgLvl()

bool AthCommonMsg< Algorithm >::msgLvl ( const MSG::Level lvl) const
inlineinherited

Definition at line 30 of file AthCommonMsg.h.

30 {
31 return this->msgLevel(lvl);
32 }

◆ outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< Algorithm > >::outputHandles ( ) const
overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ particleData()

const HepPDT::ParticleData * GenBase::particleData ( int pid) const
inlineinherited

Access an element in the particle data table.

Definition at line 126 of file GenBase.h.

126 {
127 return pdt().particle(HepPDT::ParticleID(std::abs(pid)));
128 }
const HepPDT::ParticleDataTable & pdt() const
Shorter alias to get a particle data table.
Definition GenBase.h:123

◆ particleTable()

const HepPDT::ParticleDataTable & GenBase::particleTable ( ) const
inlineinherited

Get a particle data table.

Definition at line 118 of file GenBase.h.

118 {
119 return *(m_ppSvc->PDT());
120 }
ServiceHandle< IPartPropSvc > m_ppSvc
Handle on the particle property service.
Definition GenBase.h:160

◆ partPropSvc()

const ServiceHandle< IPartPropSvc > GenBase::partPropSvc ( ) const
inlineinherited

Access the particle property service.

Definition at line 113 of file GenBase.h.

113 {
114 return m_ppSvc;
115 }

◆ passesUserSelection()

bool EvtInclusiveDecay::passesUserSelection ( HepMC::GenEvent * hepMC)
private

Definition at line 606 of file EvtInclusiveDecay.cxx.

606 {
607 bool passed(false);
608 std::vector<HepMC::GenParticlePtr> *muons = new std::vector<HepMC::GenParticlePtr>;
609
610 for ( const auto& p: *hepMC) {
611 if( std::abs(p->pdg_id()) == 13 )
612 muons->push_back(p);
613 }
614
615 for (auto muItr1 = muons->begin(); muItr1 != muons->end(); ++muItr1) {
616 for (auto muItr2 = muItr1+1; muItr2 != muons->end(); ++muItr2) {
617 if( m_userSelRequireOppositeSignedMu && (*muItr1)->pdg_id() * (*muItr2)->pdg_id() > 0)
618 continue;
619 if( !( (*muItr1)->momentum().perp() > m_userSelMu1MinPt && std::abs((*muItr1)->momentum().pseudoRapidity()) < m_userSelMu1MaxEta &&
620 (*muItr2)->momentum().perp() > m_userSelMu2MinPt && std::abs((*muItr2)->momentum().pseudoRapidity()) < m_userSelMu2MaxEta ) &&
621 !( (*muItr2)->momentum().perp() > m_userSelMu1MinPt && std::abs((*muItr2)->momentum().pseudoRapidity()) < m_userSelMu1MaxEta &&
622 (*muItr1)->momentum().perp() > m_userSelMu2MinPt && std::abs((*muItr1)->momentum().pseudoRapidity()) < m_userSelMu2MaxEta ) )
623 continue;
624 double dimuMass = invMass((*muItr1),(*muItr2));
625 if( !( dimuMass > m_userSelMinDimuMass && (dimuMass < m_userSelMaxDimuMass || m_userSelMaxDimuMass < 0.) ) )
626 continue;
627 passed = true;
628 }
629 }
630
631 delete muons;
632
633 return passed;
634}
bool passed(DecisionID id, const DecisionIDContainer &)
checks if required decision ID is in the set of IDs in the container

◆ pdgName()

std::string EvtInclusiveDecay::pdgName ( HepMC::ConstGenParticlePtr p,
bool statusHighlighting = false,
std::set< HepMC::GenParticlePtr, ParticleIdCompare > * particleSet = nullptr )
private

Definition at line 707 of file EvtInclusiveDecay.cxx.

707 {
708 std::ostringstream buf;
709 bool inlist = false;
710#ifdef HEPMC3
711 if (particleSet) for (const auto& pinl: *particleSet) if (pinl&&p) if (pinl.get() == p.get()) inlist=true;
712#else
713 auto p_nc ATLAS_THREAD_SAFE = const_cast<HepMC::GenParticlePtr> (p);
714 if (particleSet) inlist = (particleSet->find(p_nc) != particleSet->end());
715#endif
716 if (statusHighlighting) {
717 if ( ((particleSet!=0) && (inlist)) ||
718 ((particleSet==0) && isToBeDecayed(p,false)) )
719 buf << "\033[7m"; // reverse
720 if (!MC::isStable(p)) {
721 if (MC::isDecayed(p))
722 buf << "\033[33m"; // yellow
723 else
724 buf << "\033[31m"; // red
725 }
726 }
727 buf << p->pdg_id();
728 buf << "/" << HepPID::particleName(p->pdg_id());
729 if (statusHighlighting) {
730 buf << "\033[0m"; // revert color attributes
731 }
732 return buf.str();
733}
#define ATLAS_THREAD_SAFE
bool isToBeDecayed(HepMC::ConstGenParticlePtr p, bool doCrossChecks)
bool isDecayed(const T &p)
Identify if the particle decayed.

◆ pdt()

const HepPDT::ParticleDataTable & GenBase::pdt ( ) const
inlineinherited

Shorter alias to get a particle data table.

Definition at line 123 of file GenBase.h.

123{ return particleTable(); }
const HepPDT::ParticleDataTable & particleTable() const
Get a particle data table.
Definition GenBase.h:118

◆ printHepMC()

void EvtInclusiveDecay::printHepMC ( HepMC::GenEvent * hepMC,
std::set< HepMC::GenParticlePtr, ParticleIdCompare > * particleSet = nullptr )
private

Definition at line 660 of file EvtInclusiveDecay.cxx.

660 {
661 std::set<HepMC::GenVertexPtr> visited;
662 unsigned int nParticlesFound = 0;
663 unsigned int nTreesFound = 0;
664 for (auto p: *hepMC) {
665 if ( (!p->production_vertex()) ||
666 (p->production_vertex()->particles_in_size() == 0) ) {
667 nTreesFound++;
668 std::cout << "\n Found new partial decay tree:\n" << std::endl;
669 unsigned int nParticlesVisited = printTree(std::move(p),visited,1,particleSet);
670 std::cout << "\n " << nParticlesVisited << " particles in this subtree" << std::endl;
671 nParticlesFound += nParticlesVisited;
672 }
673 }
674 std::cout << "\n Total of " << nParticlesFound << " particles found in "
675 << nTreesFound << " decay subtrees in HepMC event record\n" << std::endl;
676}
unsigned int printTree(HepMC::GenParticlePtr p, std::set< HepMC::GenVertexPtr > &visited, int level, std::set< HepMC::GenParticlePtr, ParticleIdCompare > *particleSet=nullptr)

◆ printTree()

unsigned int EvtInclusiveDecay::printTree ( HepMC::GenParticlePtr p,
std::set< HepMC::GenVertexPtr > & visited,
int level,
std::set< HepMC::GenParticlePtr, ParticleIdCompare > * particleSet = nullptr )
private

Definition at line 678 of file EvtInclusiveDecay.cxx.

678 {
679
680 unsigned int nParticlesVisited = 1;
681 for (int i=0; i<level; i++) std::cout << " ";
682 std::cout << pdgName(p,m_printHepMCHighlighted,particleSet);
683 auto v = p->end_vertex();
684 if (v) {
685 if (v->particles_in_size() > 1)
686 std::cout << " [interaction: " << v->particles_in_size() << " particles, vertex " << v << "] --> ";
687 else
688 std::cout << " --> ";
689 if (visited.insert(v).second) {
690 for (auto itp: *v) {
691 std::cout << pdgName(itp,m_printHepMCHighlighted,particleSet) << " ";
692 }
693 std::cout << std::endl;
694 for (auto itp: *v) {
695 if (itp->end_vertex())
696 nParticlesVisited += printTree(std::move(itp), visited, level+1, particleSet);
697 else
698 nParticlesVisited++;
699 }
700 } else
701 std::cout << "see above" << std::endl;
702 } else
703 std::cout << " no decay vertex\n" << std::endl;
704 return nParticlesVisited;
705}

◆ removeDecayTree()

void EvtInclusiveDecay::removeDecayTree ( HepMC::GenEvent * hepMC,
HepMC::GenParticlePtr p )
private

Definition at line 396 of file EvtInclusiveDecay.cxx.

396 {
397 auto v = p->end_vertex();
398 if (v) {
399#ifdef HEPMC3
400 //This is recursive in HepMC3. But explicit deletion is allowed as well.
401 hepMC->remove_vertex(std::move(v));
402 p->set_status(1); // For now, flag particle as undecayed (stable)
403 ATH_MSG_DEBUG("Removed existing " << pdgName(p) << " " << p );
404#else
405 std::set<int> vtxBarCodesToDelete;
406 vtxBarCodesToDelete.insert(v->barcode());
407 for (HepMC::GenVertex::vertex_iterator itv = v->vertices_begin(HepMC::descendants);
408 itv != v->vertices_end(HepMC::descendants);
409 ++itv)
410 vtxBarCodesToDelete.insert((*itv)->barcode());
411 for (std::set<int>::iterator itb = vtxBarCodesToDelete.begin(); itb != vtxBarCodesToDelete.end(); ++itb) {
412 auto vdel = hepMC->barcode_to_vertex(*itb);
413 hepMC->remove_vertex(vdel);
414 delete vdel;
415 }
416 p->set_status(1); // For now, flag particle as undecayed (stable)
417 ATH_MSG_DEBUG("Removed existing " << pdgName(p) << " (barcode " << p->barcode() << ")" << " decay tree with " << vtxBarCodesToDelete.size() << " vertices");
418#endif
419 }
420}

◆ renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< Algorithm > >::renounce ( T & h)
inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

381 {
382 h.renounce();
384 }
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce(T &h)

◆ renounceArray()

void AthCommonDataStore< AthCommonMsg< Algorithm > >::renounceArray ( SG::VarHandleKeyArray & handlesArray)
inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

364 {
366 }

◆ reseedRandomEngine()

void EvtInclusiveDecay::reseedRandomEngine ( const std::string & streamName,
const EventContext & ctx )
private

Definition at line 183 of file EvtInclusiveDecay.cxx.

185{
186 long seeds[7];
187 ATHRNG::calculateSeedsMC21(seeds, streamName, ctx.eventID().event_number(), m_dsid, m_randomSeed);
188 m_evtAtRndmGen->getEngine()->setSeeds(seeds, 0); // NOT THREAD-SAFE
189}
void calculateSeedsMC21(long *seeds, const std::string &algName, uint64_t ev, uint64_t run, uint64_t offset=0)
Set the random seed using a string (e.g.

◆ sysInitialize()

StatusCode AthAlgorithm::sysInitialize ( )
overridevirtualinherited

Override sysInitialize.

Override sysInitialize from the base class.

Loop through all output handles, and if they're WriteCondHandles, automatically register them and this Algorithm with the CondSvc

Scan through all outputHandles, and if they're WriteCondHandles, register them with the CondSvc

Reimplemented from AthCommonDataStore< AthCommonMsg< Algorithm > >.

Reimplemented in AthAnalysisAlgorithm, AthFilterAlgorithm, AthHistogramAlgorithm, and PyAthena::Alg.

Definition at line 66 of file AthAlgorithm.cxx.

66 {
68
69 if (sc.isFailure()) {
70 return sc;
71 }
72 ServiceHandle<ICondSvc> cs("CondSvc",name());
73 for (auto h : outputHandles()) {
74 if (h->isCondition() && h->mode() == Gaudi::DataHandle::Writer) {
75 // do this inside the loop so we don't create the CondSvc until needed
76 if ( cs.retrieve().isFailure() ) {
77 ATH_MSG_WARNING("no CondSvc found: won't autoreg WriteCondHandles");
78 return StatusCode::SUCCESS;
79 }
80 if (cs->regHandle(this,*h).isFailure()) {
81 sc = StatusCode::FAILURE;
82 ATH_MSG_ERROR("unable to register WriteCondHandle " << h->fullKey()
83 << " with CondSvc");
84 }
85 }
86 }
87 return sc;
88}
#define ATH_MSG_WARNING(x)
virtual StatusCode sysInitialize() override
Override sysInitialize.
AthCommonDataStore(const std::string &name, T... args)
virtual std::vector< Gaudi::DataHandle * > outputHandles() const override

◆ sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< Algorithm > >::sysStart ( )
overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

◆ traverseDecayTree()

StatusCode EvtInclusiveDecay::traverseDecayTree ( HepMC::GenParticlePtr p,
bool isToBeRemoved,
std::set< HepMC::GenVertexPtr > & visited,
std::set< HepMC::GenParticlePtr, ParticleIdCompare > & toBeDecayed )
private

Definition at line 357 of file EvtInclusiveDecay.cxx.

360 {
361 ATH_MSG_VERBOSE("Inspecting: " << pdgName(p) << " " << p);
362 if (!isToBeRemoved) {
363 if (isToBeDecayed(p,true)) {
364 toBeDecayed.insert(p);
365 isToBeRemoved = true;
366 ATH_MSG_VERBOSE("Selected particle for decay: " << pdgName(p) << " " << p );
367
368 // In principle we could stop the recursion here. However, to prevent
369 // pathological cases in certain decay trees (in particular from Herwig),
370 // we continue in order to mark all descendants of this particle
371 // as visited. Thus none of these descendants can be flagged for further
372 // decay, even if it has several mothers.
373 }
374 }
375 auto v = p->end_vertex();
376 if (v) {
377 if (visited.insert(v).second) {
378 if ( isToBeRemoved && (v->particles_in_size()>1) && m_checkDecayTree ) {
379 ATH_MSG_WARNING("Found particle to be decayed with vertex with >1 incoming mother particles in decay tree");
380 ATH_MSG_WARNING( ([&p, &v](){ std::stringstream ss; HepMC::Print::line(ss,p); HepMC::Print::line(ss,v); return ss.str();})());
381 }
382 for (auto itp: *v) {
383 ATH_CHECK(traverseDecayTree(std::move(itp),isToBeRemoved,visited,toBeDecayed) );
384 }
385
386 }
387 }
388 return StatusCode::SUCCESS;
389}
static Double_t ss

◆ updateVHKA()

void AthCommonDataStore< AthCommonMsg< Algorithm > >::updateVHKA ( Gaudi::Details::PropertyBase & )
inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

308 {
309 // debug() << "updateVHKA for property " << p.name() << " " << p.toString()
310 // << " size: " << m_vhka.size() << endmsg;
311 for (auto &a : m_vhka) {
313 for (auto k : keys) {
314 k->setOwner(this);
315 }
316 }
317 }
std::vector< SG::VarHandleKeyArray * > m_vhka

◆ xmlpath()

std::string EvtInclusiveDecay::xmlpath ( void )

Definition at line 749 of file EvtInclusiveDecay.cxx.

749 {
750 return PathResolverFindCalibDirectory( "Pythia8/xmldoc" );
751}
std::string PathResolverFindCalibDirectory(const std::string &logical_file_name)

Member Data Documentation

◆ m_allowAllKnownDecays

bool EvtInclusiveDecay::m_allowAllKnownDecays
private

Definition at line 126 of file EvtInclusiveDecay.h.

◆ m_allowDefaultBDecays

bool EvtInclusiveDecay::m_allowDefaultBDecays
private

Definition at line 127 of file EvtInclusiveDecay.h.

◆ m_applyUserSelection

bool EvtInclusiveDecay::m_applyUserSelection
private

Definition at line 144 of file EvtInclusiveDecay.h.

◆ m_blackList

std::vector<int> EvtInclusiveDecay::m_blackList
private

Definition at line 123 of file EvtInclusiveDecay.h.

◆ m_blackListSet

std::set<int> EvtInclusiveDecay::m_blackListSet
private

Definition at line 124 of file EvtInclusiveDecay.h.

◆ m_checkDecayChannels

bool EvtInclusiveDecay::m_checkDecayChannels
private

Definition at line 137 of file EvtInclusiveDecay.h.

◆ m_checkDecayTree

bool EvtInclusiveDecay::m_checkDecayTree
private

Definition at line 136 of file EvtInclusiveDecay.h.

◆ m_decayFile

std::string EvtInclusiveDecay::m_decayFile
private

Definition at line 112 of file EvtInclusiveDecay.h.

◆ m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Algorithm > >::m_detStore
privateinherited

Pointer to StoreGate (detector store by default)

Definition at line 393 of file AthCommonDataStore.h.

◆ m_dsid

IntegerProperty EvtInclusiveDecay::m_dsid {this, "Dsid", 999999}
private

Definition at line 99 of file EvtInclusiveDecay.h.

99{this, "Dsid", 999999};

◆ m_evtAtRndmGen

EvtInclusiveAtRndmGen* EvtInclusiveDecay::m_evtAtRndmGen {}
private

Definition at line 107 of file EvtInclusiveDecay.h.

107{};

◆ m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Algorithm > >::m_evtStore
privateinherited

Pointer to StoreGate (event store by default)

Definition at line 390 of file AthCommonDataStore.h.

◆ m_extendedExtraObjects

DataObjIDColl AthAlgorithm::m_extendedExtraObjects
privateinherited

Definition at line 79 of file AthAlgorithm.h.

◆ m_inputKeyName

std::string EvtInclusiveDecay::m_inputKeyName
private

Definition at line 115 of file EvtInclusiveDecay.h.

◆ m_isfHerwig

bool EvtInclusiveDecay::m_isfHerwig
private

Definition at line 152 of file EvtInclusiveDecay.h.

◆ m_maxNRepeatedDecays

int EvtInclusiveDecay::m_maxNRepeatedDecays
private

Definition at line 142 of file EvtInclusiveDecay.h.

◆ m_mcEventKey

std::string GenBase::m_mcEventKey {}
protectedinherited

StoreGate key for the MC event collection (defaults to GEN_EVENT)

Definition at line 137 of file GenBase.h.

137{};

◆ m_mcevents_const

SG::ReadHandleKey<McEventCollection> GenBase::m_mcevents_const { this, "McEventKey", "GEN_EVENT", "StoreGate key of the MC event collection" }
privateinherited

Const handle to the MC event collection.

Definition at line 163 of file GenBase.h.

163{ this, "McEventKey", "GEN_EVENT", "StoreGate key of the MC event collection" };

◆ m_mcEvtColl

McEventCollection* EvtInclusiveDecay::m_mcEvtColl {}
private

Definition at line 104 of file EvtInclusiveDecay.h.

104{};

◆ m_mkMcEvent

BooleanProperty GenBase::m_mkMcEvent {this, "MakeMcEvent", false, "Create a new MC event collection if it doesn't exist"}
protectedinherited

Flag to determine if a new MC event collection should be made if it doesn't exist.

Definition at line 139 of file GenBase.h.

139{this, "MakeMcEvent", false, "Create a new MC event collection if it doesn't exist"};

◆ m_myEvtGen

EvtGen* EvtInclusiveDecay::m_myEvtGen {}
private

Definition at line 108 of file EvtInclusiveDecay.h.

108{};

◆ m_noDecayChannels

std::map<int,long> EvtInclusiveDecay::m_noDecayChannels
private

Definition at line 138 of file EvtInclusiveDecay.h.

◆ m_nRepeatedDecays

int EvtInclusiveDecay::m_nRepeatedDecays
private

Definition at line 140 of file EvtInclusiveDecay.h.

◆ m_outputKeyName

std::string EvtInclusiveDecay::m_outputKeyName
private

Definition at line 116 of file EvtInclusiveDecay.h.

◆ m_pdtFile

std::string EvtInclusiveDecay::m_pdtFile
private

Definition at line 111 of file EvtInclusiveDecay.h.

◆ m_ppSvc

ServiceHandle<IPartPropSvc> GenBase::m_ppSvc {this, "PartPropSvc", "PartPropSvc"}
privateinherited

Handle on the particle property service.

Definition at line 160 of file GenBase.h.

160{this, "PartPropSvc", "PartPropSvc"};

◆ m_printHepMCAfterEvtGen

bool EvtInclusiveDecay::m_printHepMCAfterEvtGen
private

Definition at line 132 of file EvtInclusiveDecay.h.

◆ m_printHepMCBeforeEvtGen

bool EvtInclusiveDecay::m_printHepMCBeforeEvtGen
private

Definition at line 131 of file EvtInclusiveDecay.h.

◆ m_printHepMCHighlighted

bool EvtInclusiveDecay::m_printHepMCHighlighted
private

Definition at line 133 of file EvtInclusiveDecay.h.

◆ m_printHepMCHighLightTopLevelDecays

bool EvtInclusiveDecay::m_printHepMCHighLightTopLevelDecays
private

Definition at line 134 of file EvtInclusiveDecay.h.

◆ m_prohibitFinalStateDecay

bool EvtInclusiveDecay::m_prohibitFinalStateDecay
private

Definition at line 119 of file EvtInclusiveDecay.h.

◆ m_prohibitReDecay

bool EvtInclusiveDecay::m_prohibitReDecay
private

Definition at line 120 of file EvtInclusiveDecay.h.

◆ m_prohibitRemoveSelfDecay

bool EvtInclusiveDecay::m_prohibitRemoveSelfDecay
private

Definition at line 122 of file EvtInclusiveDecay.h.

◆ m_prohibitUnDecay

bool EvtInclusiveDecay::m_prohibitUnDecay
private

Definition at line 121 of file EvtInclusiveDecay.h.

◆ m_randomSeed

IntegerProperty EvtInclusiveDecay::m_randomSeed {this, "RandomSeed", 1234567, "Random seed for the built-in random engine"}
private

Seed for random number engine.

Definition at line 102 of file EvtInclusiveDecay.h.

102{this, "RandomSeed", 1234567, "Random seed for the built-in random engine"}; // FIXME make this into an unsigned long int?

◆ m_randomStreamName

std::string EvtInclusiveDecay::m_randomStreamName
private

Definition at line 114 of file EvtInclusiveDecay.h.

◆ m_readExisting

bool EvtInclusiveDecay::m_readExisting
private

Definition at line 118 of file EvtInclusiveDecay.h.

◆ m_rndmSvc

ServiceHandle<IAthRNGSvc> EvtInclusiveDecay::m_rndmSvc {this, "RndmSvc", "AthRNGSvc"}
private

Definition at line 96 of file EvtInclusiveDecay.h.

96{this, "RndmSvc", "AthRNGSvc"};

◆ m_setVMtransversePol

bool EvtInclusiveDecay::m_setVMtransversePol
private

Definition at line 153 of file EvtInclusiveDecay.h.

◆ m_userDecayFile

std::string EvtInclusiveDecay::m_userDecayFile
private

Definition at line 113 of file EvtInclusiveDecay.h.

◆ m_userSelMaxDimuMass

double EvtInclusiveDecay::m_userSelMaxDimuMass
private

Definition at line 151 of file EvtInclusiveDecay.h.

◆ m_userSelMinDimuMass

double EvtInclusiveDecay::m_userSelMinDimuMass
private

Definition at line 150 of file EvtInclusiveDecay.h.

◆ m_userSelMu1MaxEta

double EvtInclusiveDecay::m_userSelMu1MaxEta
private

Definition at line 148 of file EvtInclusiveDecay.h.

◆ m_userSelMu1MinPt

double EvtInclusiveDecay::m_userSelMu1MinPt
private

Definition at line 146 of file EvtInclusiveDecay.h.

◆ m_userSelMu2MaxEta

double EvtInclusiveDecay::m_userSelMu2MaxEta
private

Definition at line 149 of file EvtInclusiveDecay.h.

◆ m_userSelMu2MinPt

double EvtInclusiveDecay::m_userSelMu2MinPt
private

Definition at line 147 of file EvtInclusiveDecay.h.

◆ m_userSelRequireOppositeSignedMu

bool EvtInclusiveDecay::m_userSelRequireOppositeSignedMu
private

Definition at line 145 of file EvtInclusiveDecay.h.

◆ m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< Algorithm > >::m_varHandleArraysDeclared
privateinherited

Definition at line 399 of file AthCommonDataStore.h.

◆ m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< Algorithm > >::m_vhka
privateinherited

Definition at line 398 of file AthCommonDataStore.h.

◆ m_whiteList

std::vector<int> EvtInclusiveDecay::m_whiteList
private

Definition at line 128 of file EvtInclusiveDecay.h.

◆ m_whiteListSet

std::set<int> EvtInclusiveDecay::m_whiteListSet
private

Definition at line 129 of file EvtInclusiveDecay.h.


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