FixHepMC Class Reference

A "fix-up" algorithm to correct weird event records. More...

#include <FixHepMC.h>

Inheritance diagram for FixHepMC:

Collaboration diagram for FixHepMC:

Public Member Functions
	FixHepMC (const std::string &name, ISvcLocator *pSvcLocator)
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 reduction functions
StatusCode	execute ()
StatusCode	finalize ()
Event loop algorithm methods
virtual StatusCode	initialize () override
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 McEventCollection *	events_const () const
	Access the current event's McEventCollection (const)
const McEventCollection *	events_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.

Private Types
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
	specialization for handling Gaudi::Property<SG::VarHandleKey>
Classifiers for identifying particles to be removed
bool	isPID0 (const HepMC::ConstGenParticlePtr &p) const
bool	isNonTransportableInDecayChain (const HepMC::ConstGenParticlePtr &p) const
bool	isSimpleLoop (const HepMC::ConstGenParticlePtr &p) const
bool	fromDecay (const HepMC::ConstGenParticlePtr &p, std::shared_ptr< std::set< int > > &storage) const

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
Config properties
bool	m_killLoops
bool	m_killPDG0
bool	m_cleanDecays
bool	m_purgeUnstableWithoutEndVtx
bool	m_ignoreSemiDisconnected
std::string	m_forced_momentum {""}
std::string	m_forced_length {""}
bool	m_unitsFix
bool	m_setHasCycles
Cleaned-particle counters
Tell HEPMC3 that there are cycles (loops) in these events
long	m_loopKilled
long	m_pdg0Killed
long	m_decayCleaned
long	m_unstablePurged
long	m_totalSeen
long	m_replacedPIDs
std::map< int, int >	m_pidmap
	map of pids to change.
std::map< int, int >	m_replacedpid_counts
	map of counters of replacements.
MC::Loops< HepMC::GenEvent, HepMC::ConstGenParticlePtr, HepMC::ConstGenVertexPtr >	m_looper
	member to detect loops

Properties
ServiceHandle< IPartPropSvc >	m_ppSvc {this, "PartPropSvc", "PartPropSvc"}
	Handle on the particle property service.
SG::ReadHandleKey< McEventCollection >	m_mcevents_const { this, "McEventKey", "GEN_EVENT", "StoreGate key of the MC event collection" }
	Const handle to the MC event collection.
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.

Detailed Description

A "fix-up" algorithm to correct weird event records.

FixHepMC is run in evgen production, after the chain of input, shower, hadronization, and afterburner generators, but before the TestHepMC check algorithm, any analysis algs, and the POOL event stream output. Its job is to "normalise" strange event features to allow greater consistency of MC event record use in ATLAS.

Examples of its features include removing dummy "padding" particles from event records, stripping parton and EW boson intermediates from hadron/tau decay chains, and ensuring that every event has at least one entry in its event weight vector.

Definition at line 31 of file FixHepMC.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

FixHepMC()

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

Definition at line 14 of file FixHepMC.cxx.

  : GenBase(name, pSvcLocator)
  , m_loopKilled(0)
  , m_pdg0Killed(0)
  , m_decayCleaned(0)
  , m_unstablePurged(0)
  , m_totalSeen(0)
  , m_replacedPIDs(0)
{
  declareProperty("KillLoops", m_killLoops = true, "Remove particles in loops?");
  declareProperty("KillPDG0", m_killPDG0 = true, "Remove particles with PDG ID 0?");
  declareProperty("CleanDecays", m_cleanDecays = true, "Clean decay chains from non-propagating particles?");
  declareProperty("PurgeUnstableWithoutEndVtx", m_purgeUnstableWithoutEndVtx = false, "Remove unstable particles without decay vertex?");
  declareProperty("IgnoreSemiDisconnected", m_ignoreSemiDisconnected = false, "Ignore semi-disconnected particles (normal in Sherpa)");
  declareProperty("PIDmap", m_pidmap = std::map<int,int>(), "Map of PDG IDs to replace");
  declareProperty("forced_momentum", m_forced_momentum = "MEV", "Forced momentum unit");
  declareProperty("forced_length", m_forced_length = "MM", "Forced length unit");
  declareProperty("ApplyUnitsFix", m_unitsFix = true, "Attempt to identify momentum units problems and fix them");
  declareProperty("SetHasCycles", m_setHasCycles = false, "Inform HEPMC3 that this event has cycles (loops)");
}

Member Function Documentation

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.

                                                      {
    if (events()->empty())
      ATH_MSG_ERROR("McEventCollection is empty during first event access");
    return *(events()->begin());
  }

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.

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.

  {
    return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                       hndl.value(), 
                                                       hndl.documentation());
 
  }

declareProperty()

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

inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

                                                                     {
    typedef typename SG::HandleClassifier<T>::type htype;
    return AthCommonDataStore<PBASE>::declareGaudiProperty(t, htype());
  }

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.

{ return m_detStore; }

event_const()

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

inlineinherited

Access the current signal event (const)

Definition at line 83 of file GenBase.h.

                                           {
    if (events_const()->empty())
      ATH_MSG_ERROR("Const McEventCollection is empty during first event access");
    return *(events_const()->begin());
  }

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.

                                                {
    return events_const( getContext() );
  }

events_const() [2/2]

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

inlineinherited

Definition at line 99 of file GenBase.h.

                                                                         {
    SG::ReadHandle<McEventCollection> ret = SG::makeHandle(m_mcevents_const, ctx);
    if (!ret.isValid())
      ATH_MSG_ERROR("No McEventCollection found in StoreGate with key " << m_mcevents_const.key());
    return ret.cptr();
  }

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.

{ return m_evtStore; }

execute()

StatusCode FixHepMC::execute ( )

virtual

Case 1: particles without production vertex, except beam particles (status 4)

Case 2: non-final-state particles without end vertex

AV: In case we have 3 particles, we try to add a vertex that correspond to 1->2 and 1->1 splitting. AV: we can try to do that for 4 particles as well.

The condition below will cover 1->1, 1->2 and 2->1 cases

Remove loops inside decay chains AV: Please note that this approach would distort some branching ratios. If some particle would have decay products with bad PDG ids, after the operation below the visible branching ratio of these decays would be zero.

Check the bad particles have prod and end vertices

Check that all particles coming into the decay vertex of bad particle cam from the same production vertex.

Reimplemented from GenBase.

Definition at line 98 of file FixHepMC.cxx.

                             {
  for (McEventCollection::const_iterator ievt = events()->begin(); ievt != events()->end(); ++ievt) {
    // FIXME: const_cast
    HepMC::GenEvent* evt = const_cast<HepMC::GenEvent*>(*ievt);
    m_looper.findLoops(evt,true);
    /* AV: To understand why some algorithms fail one would have to request debug explicitly, but that is the reviwers consensus.*/
    if (!m_looper.loop_particles().empty() || !m_looper.loop_vertices().empty()) {
      ATH_MSG_DEBUG("Found " << m_looper.loop_vertices().size() << " vertices in loops");
      ATH_MSG_DEBUG("Found " << m_looper.loop_particles().size() << " particles in loops");
      ATH_MSG_DEBUG("Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
    }
#ifdef HEPMC3
    auto old_momentum = evt->momentum_unit();
    auto old_length = evt->length_unit();
    evt->set_units(
        m_forced_momentum != "" ? HepMC3::Units::momentum_unit(m_forced_momentum): old_momentum,
        m_forced_length != "" ? HepMC3::Units::length_unit(m_forced_length): old_length);
    if ( m_forced_momentum != "" && m_forced_momentum != HepMC3::Units::name(old_momentum) ) ATH_MSG_WARNING("Updated momentum units " <<  HepMC3::Units::name(old_momentum) << "->" << m_forced_momentum);
    if ( m_forced_length != "" && m_forced_length != HepMC3::Units::name(old_length) ) ATH_MSG_WARNING("Updated length units " <<  HepMC3::Units::name(old_length) << "->" << m_forced_length);
#else
    std::string old_momentum = (evt->momentum_unit() == HepMC::Units::MomentumUnit::MEV ? "MEV" : "GEV" );
    std::string old_length = (evt->length_unit() == HepMC::Units::LengthUnit::MM ? "MM" : "CM" );
    evt->define_units(
        m_forced_momentum != "" ? m_forced_momentum: old_momentum,
        m_forced_length != "" ? m_forced_length: old_length);
    if ( m_forced_momentum != "" && m_forced_momentum != old_momentum ) ATH_MSG_WARNING("Updated momentum units " <<  old_momentum << "->" << m_forced_momentum);
    if ( m_forced_length != "" && m_forced_length != old_length ) ATH_MSG_WARNING("Updated length units " <<  old_length << "->" << m_forced_length);
#endif
 
    if (!m_pidmap.empty()) {
      for (auto ip: *evt) {
        // Skip this particle if (somehow) its pointer is null
         if (!ip) continue;
         auto newpid = m_pidmap.find(ip->pdg_id());
         if (newpid == m_pidmap.end()) continue;
         ip->set_pdg_id(newpid->second);
         // Increment the counter for replaced PDG IDs
         ++m_replacedPIDs;
         m_replacedpid_counts[ip->pdg_id()]++;
      }
    }
 
#ifdef HEPMC3
    if (m_setHasCycles){
      // If asked, tag the event as having cycles and alert the user that this problem exists
      auto cycles = std::make_shared<HepMC3::IntAttribute>(1);
      evt->add_attribute("cycles",cycles);
    }
#endif
 
#ifdef HEPMC3
    // Add a unit entry to the event weight vector if it's currently empty
    if (evt->weights().empty()) {
      ATH_MSG_DEBUG("Adding a unit weight to empty event weight vector");
      evt->weights().push_back(1);
    }
    // Set a (0,0,0) vertex to be the signal vertex if not already set
    if (!HepMC::signal_process_vertex(evt)) {
      const HepMC::FourVector nullpos;
      for (auto  iv: evt->vertices()) {
        if (iv->position() == nullpos) {
          ATH_MSG_DEBUG("Setting representative event position vertex");
          HepMC::set_signal_process_vertex(evt,iv);
          break;
        }
      }
    }
 
    // Catch cases with more than 2 beam particles (16.11.2021)
    std::vector<HepMC::GenParticlePtr> beams_t;
    for (const HepMC::GenParticlePtr& p : evt->beams()) {
      if (p->status() == 4)  beams_t.push_back(p);
    }
    if (beams_t.size() > 2) {
      ATH_MSG_INFO("Invalid number of beam particles " <<  beams_t.size() << ". Will try to fix.");
      std::vector<HepMC::GenParticlePtr> bparttoremove;
      for (const auto& bpart : beams_t) {
        if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
      }
      for (auto bpart: bparttoremove) {
        bpart->production_vertex()->remove_particle_out(bpart);
      }
    }
 
    // Some generators / samples have a mis-match between the units they report and the units used for momentum, usually because we have applied
    // a correction to the momentum without also correcting the units. Here we're going to check for states that look particularly suspicious
    // and apply a correction if required. We will try to identify those issues based on the beam particles.
    if (m_unitsFix){ // Only if requested - allow folks to disable this if they know what they're doing
      double units_problem = -1.;
      // If we have beam particles, let's use them - it's much faster!
      if (beams_t.size() > 0 ){
        for (const HepMC::GenParticlePtr& p : beams_t){
          if (p->momentum().pz() > 1e9) units_problem = p->momentum().pz();
        }
      // if we didn't have beam particles, we'll go through some of the main record
      } else {
        for (const HepMC::GenParticlePtr& p : evt->particles()) {
          if (p && p->momentum().pz() > 1e9) units_problem = p->momentum().pz();
        }
      }
      if (units_problem>0){ // No particles should have momenta above 1 PeV; this must be a units issue
        ATH_MSG_INFO("Apparent units problem; beam particles have z-momentum " << units_problem << " in MeV. Will divide by 1000.");
        MC::MeVToGeV(evt);  //Only scales momenta and masses
      }
    }
 
    // Some heuristics to catch problematic cases
    std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
    for (auto ip : evt->particles()) {
      // Skip this particle if (somehow) its pointer is null
      if (!ip) continue;
      bool particle_to_fix = false;
      int abspid = std::abs(ip->pdg_id());
      auto vProd = ip->production_vertex();
      auto vEnd  = ip->end_vertex();
      if ( (!vProd || vProd->id() == 0) && vEnd && ip->status() != 4) {
        particle_to_fix = true;
        ATH_MSG_DEBUG("Found particle " << ip->pdg_id() << " without production vertex! HepMC status = " << ip->status());
      }
      if (vProd && !vEnd && ip->status() != 1) {
        particle_to_fix = true;
        ATH_MSG_DEBUG("Found particle " << ip->pdg_id() << " without decay vertex! HepMC status = " << ip->status());
      }
      if (particle_to_fix)  semi_disconnected.push_back(ip);
      // Case 3: keep track of loop particles inside decay chains (seen in H7+EvtGen)
      if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
        decay_loop_particles.push_back(std::move(ip));
      }
    }


    if ( !m_ignoreSemiDisconnected && (semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2)) {
      size_t no_endv = 0;
      size_t no_prov = 0;
      HepMC::FourVector sum(0,0,0,0);
      std::set<HepMC::GenVertexPtr> standalone;
      for (const auto& part : semi_disconnected) {
        if (!part->production_vertex() || !part->production_vertex()->id()) {
          no_prov++; sum += part->momentum();
        }
        if (!part->end_vertex()) { 
          no_endv++;  sum -= part->momentum();
        }
        if (part->production_vertex()) standalone.insert(part->production_vertex());
        if (part->end_vertex()) standalone.insert(part->end_vertex());
      }
      ATH_MSG_INFO("Heuristics: found " << semi_disconnected.size() << " semi-disconnected particles. Momentum sum is " << sum << " Standalone vertices " << standalone.size());
      bool standalonevertex = (standalone.size() == 1 && (*standalone.begin())->particles_in().size() + (*standalone.begin())->particles_out().size() == semi_disconnected.size());
      if (! standalonevertex && no_endv && no_prov  && ( no_endv + no_prov  == semi_disconnected.size() )) {
        if (std::abs(sum.px()) < 1e-2  && std::abs(sum.py()) < 1e-2  && std::abs(sum.pz()) < 1e-2 ) {
          ATH_MSG_INFO("Try " << no_endv << "->" << no_prov << " splitting/merging.");
          auto v = HepMC::newGenVertexPtr();
          for (auto part : semi_disconnected) {
            if (!part->production_vertex() || part->production_vertex()->id() == 0) v->add_particle_out(std::move(part));
          }
          for (auto part : semi_disconnected) {
            if (!part->end_vertex()) v->add_particle_in(std::move(part));
          }
          evt->add_vertex(std::move(v));
        }
      }
    }

    for (auto part: decay_loop_particles) {
      auto vend = part->end_vertex();
      auto vprod = part->production_vertex();
      if (!vprod || !vend)  continue;
      bool loop_in_decay = true;
      auto sisters = vend->particles_in();
      for (auto sister: sisters) {
        if (vprod != sister->production_vertex()) loop_in_decay = false;
      }
      if (!loop_in_decay) continue;

      auto daughters = vend->particles_out();
      for (auto p : daughters) vprod->add_particle_out(std::move(p));
      for (auto sister : sisters) { 
        vprod->remove_particle_out(sister); 
        vend->remove_particle_in(sister); 
        evt->remove_particle(std::move(sister));
      }
      evt->remove_vertex(std::move(vend));
 
    }
 
    // Event particle content cleaning -- remove "bad" structures
    std::vector<HepMC::GenParticlePtr> toremove;
    for (auto ip: evt->particles()) {
      // Skip this particle if (somehow) its pointer is null
      if (!ip) continue;
      m_totalSeen += 1;
      // Flag to declare if a particle should be removed
      bool bad_particle = false;
      // Check for loops
      if ( m_killLoops && isSimpleLoop(ip) ) {
        bad_particle = true;
        m_loopKilled += 1;
        ATH_MSG_DEBUG( "Found a looper : " );
        if ( msgLvl( MSG::DEBUG ) ) HepMC::Print::line(ip);
      }
      // Check on PDG ID 0
      if ( m_killPDG0 && isPID0(ip) ) {
        bad_particle = true;
        m_pdg0Killed += 1;
        ATH_MSG_DEBUG( "Found PDG ID 0 : " );
        if ( msgLvl( MSG::DEBUG ) )HepMC::Print::line(ip);
      }
      // Clean decays
      int abs_pdg_id = std::abs(ip->pdg_id());
      bool is_decayed_weak_boson =  ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && ip->end_vertex();
      if ( m_cleanDecays && isNonTransportableInDecayChain(ip) && !is_decayed_weak_boson ) {
        bad_particle = true;
        m_decayCleaned += 1;
        ATH_MSG_DEBUG( "Found a bad particle in a decay chain : " );
        if ( msgLvl( MSG::DEBUG ) ) HepMC::Print::line(ip);
      }
      // Only add to the toremove vector once, even if multiple tests match
      if (bad_particle) toremove.push_back(std::move(ip));
    }
 
    // Properties before cleaning
    const int num_particles_orig = evt->particles().size();
 
    // Do the cleaning
    if (!toremove.empty()) {
      ATH_MSG_DEBUG("Cleaning event record of " << toremove.size() << " bad particles");
      for (auto part: toremove) evt->remove_particle(std::move(part));
    }
 
    if(m_purgeUnstableWithoutEndVtx) {
      int purged=0;
      do {
        purged=0;
        const std::vector <HepMC::GenParticlePtr> allParticles=evt->particles();
        for(auto p : allParticles) {
          HepMC::ConstGenVertexPtr end_v=p->end_vertex();
          if(p->status() == 2 && !end_v) {
            evt->remove_particle(std::move(p));
            ++purged;
            ++m_unstablePurged;
          } 
        }
      }
      while (purged>0);
    }
 
    const int num_particles_filt = evt->particles().size();
 
    if(num_particles_orig!=num_particles_filt) {
      // Write out the change in the number of particles
      ATH_MSG_INFO("Particles filtered: " << num_particles_orig << " -> " << num_particles_filt);
    }
 #else
 
    // Add a unit entry to the event weight vector if it's currently empty
    if (evt->weights().empty()) {
      ATH_MSG_DEBUG("Adding a unit weight to empty event weight vector");
      evt->weights().push_back(1);
    }
 
    // Set a (0,0,0) vertex to be the signal vertex if not already set
    if (evt->signal_process_vertex() == NULL) {
      const HepMC::FourVector nullpos;
      for (HepMC::GenEvent::vertex_const_iterator iv = evt->vertices_begin(); iv != evt->vertices_end(); ++iv) {
        if ((*iv)->position() == nullpos) {
          ATH_MSG_DEBUG("Setting representative event position vertex");
          evt->set_signal_process_vertex(const_cast<HepMC::GenVertex*>(*iv));
          break;
        }
      }
    }
 
 
    // Some heuristics to catch problematic cases
    std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
    for (auto ip : *evt) {
      // Skip this particle if (somehow) its pointer is null
      if (!ip) continue;
      bool particle_to_fix = false;
      int abspid = std::abs(ip->pdg_id());
      auto vProd = ip->production_vertex();
      auto vEnd  = ip->end_vertex();
      if ( (!vProd || vProd->id() == 0) && vEnd && ip->status() != 4) {
        particle_to_fix = true;
        ATH_MSG_DEBUG("Found particle " << ip->pdg_id() << " without production vertex! HepMC status = " << ip->status());
      }
      if (vProd && !vEnd && ip->status() != 1) {
        particle_to_fix = true;
        ATH_MSG_DEBUG("Found particle " << ip->pdg_id() << " without decay vertex! HepMC status = " << ip->status());
      }
      if (particle_to_fix)  semi_disconnected.push_back(ip);
      // Case 3: keep track of loop particles inside decay chains (seen in H7+EvtGen)
      if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
        decay_loop_particles.push_back(ip);
      }
    }

    if (semi_disconnected.size() == 4 ||semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
      size_t no_endv = 0;
      size_t no_prov = 0;
      double vsum[4] = {0,0,0,0};
      for (auto part : semi_disconnected) {
        if (!part->production_vertex() ) { 
          no_prov++; 
          vsum[0] += part->momentum().px(); 
          vsum[1] += part->momentum().py(); 
          vsum[2] += part->momentum().pz(); 
          vsum[3] += part->momentum().e();
        }
        if (!part->end_vertex()) { 
          no_endv++;  
          vsum[0] -= part->momentum().px(); 
          vsum[1] -= part->momentum().py(); 
          vsum[2] -= part->momentum().pz(); 
          vsum[3] -= part->momentum().e();
        }
      }
      HepMC::FourVector sum(vsum[0],vsum[1],vsum[2],vsum[3]);
      ATH_MSG_INFO("Heuristics: found " << semi_disconnected.size() << " semi-disconnected particles. Momentum sum is " << vsum[0] << " " << vsum[1] << " " << vsum[2] << " " << vsum[3]);
      if (no_endv && no_prov  && ( no_endv + no_prov  == semi_disconnected.size() )) {
        if (std::abs(sum.px()) < 1e-2  && std::abs(sum.py()) < 1e-2  && std::abs(sum.pz()) < 1e-2 ) {
          ATH_MSG_INFO("Try " << no_endv << "->" << no_prov << " splitting/merging.");
          auto v = HepMC::newGenVertexPtr();
          for (auto part : semi_disconnected) {
            if (!part->production_vertex())  v->add_particle_out(part);
          }
          for (auto part : semi_disconnected) {
            if (!part->end_vertex()) v->add_particle_in(part);
          }
          evt->add_vertex(v);
        }
      }
    }

    for (auto part: decay_loop_particles) {
      auto vend = part->end_vertex();
      auto vprod = part->production_vertex();
      if (!vend || !vprod) continue;
      bool loop_in_decay = true;
      std::vector<HepMC::GenParticlePtr> sisters;
      for (auto p = vend->particles_begin(HepMC::parents); p!= vend->particles_end(HepMC::parents); ++p) sisters.push_back(*p);
      for (auto sister : sisters) {
        if (vprod != sister->production_vertex()) loop_in_decay = false;
      }
      if (!loop_in_decay) continue;
 
      std::vector<HepMC::GenParticlePtr> daughters;
      for (auto p = vend->particles_begin(HepMC::children); p!= vend->particles_end(HepMC::children); ++p) daughters.push_back(*p);
      for (auto p : daughters) vprod->add_particle_out(p);
      for (auto sister : sisters) { 
        vprod->remove_particle(sister); 
        vend->remove_particle(sister);
      }
      evt->remove_vertex(vend);
 
    }
 
 
    // Event particle content cleaning -- remove "bad" structures
    std::vector<HepMC::GenParticlePtr> toremove; toremove.reserve(10);
    for (HepMC::GenEvent::particle_const_iterator ip = evt->particles_begin(); ip != evt->particles_end(); ++ip) {
      // Skip this particle if (somehow) its pointer is null
      if (*ip == NULL) continue;
      m_totalSeen += 1;
 
      // Flag to declare if a particle should be removed
      bool bad_particle = false;
 
      // Check for loops
      if ( m_killLoops && isSimpleLoop(*ip) ) {
        bad_particle = true;
        m_loopKilled += 1;
        ATH_MSG_DEBUG( "Found a looper : " );
        if ( msgLvl( MSG::DEBUG ) ) (*ip)->print();
      }
 
      // Check on PDG ID 0
      if ( m_killPDG0 && isPID0(*ip) ) {
        bad_particle = true;
        m_pdg0Killed += 1;
        ATH_MSG_DEBUG( "Found PDG ID 0 : " );
        if ( msgLvl( MSG::DEBUG ) ) (*ip)->print();
      }
 
      // Clean decays
      int abs_pdg_id = std::abs((*ip)->pdg_id());
      bool is_decayed_weak_boson =  ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && (*ip)->end_vertex();
      if ( m_cleanDecays && isNonTransportableInDecayChain(*ip) && !is_decayed_weak_boson ) {
        bad_particle = true;
        m_decayCleaned += 1;
        ATH_MSG_DEBUG( "Found a bad particle in a decay chain : " );
        if ( msgLvl( MSG::DEBUG ) ) (*ip)->print();
      }
 
      // Only add to the toremove vector once, even if multiple tests match
      if (bad_particle) toremove.push_back(*ip);
    }
 
    // Properties before cleaning
    const int num_particles_orig = evt->particles_size();
    int num_orphan_vtxs_orig = 0;
    int num_noparent_vtxs_orig = 0;
    int num_nochild_vtxs_orig = 0;
    for (auto v = evt->vertices_begin(); v != evt->vertices_end(); ++v) {
      if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_orig++;
      if ((*v)->particles_in_size()==0) num_noparent_vtxs_orig++;
      if ((*v)->particles_out_size()==0) num_nochild_vtxs_orig++;
    }
 
    // Do the cleaning
    if (!toremove.empty()) {
      ATH_MSG_DEBUG("Cleaning event record of " << toremove.size() << " bad particles");
      // Clean!
      int signal_vertex_bc = evt->signal_process_vertex() ? evt->signal_process_vertex()->barcode() : 0;
      //This is the only place where reduce is used.
      reduce(evt , toremove);
      if (evt->barcode_to_vertex (signal_vertex_bc) == nullptr) {
        evt->set_signal_process_vertex (nullptr);
      }
    }
 
    if(m_purgeUnstableWithoutEndVtx) {
      int purged=0;
      do {
        for (HepMC::GenParticle* p : *evt) {
          HepMC::ConstGenVertexPtr end_v = p->end_vertex();
          if (p->status() == 2 && !end_v) {
            delete p->production_vertex()->remove_particle(p);
            ++purged;
            ++m_unstablePurged;
          }
        }
      }
      while (purged>0);
    }
 
    // Properties after cleaning
    const int num_particles_filt = evt->particles_size();
    if(num_particles_orig!=num_particles_filt) {
      int num_orphan_vtxs_filt = 0;
      int num_noparent_vtxs_filt = 0;
      int num_nochild_vtxs_filt = 0;
      for (auto v = evt->vertices_begin(); v != evt->vertices_end(); ++v) {
        if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_filt++;
        if ((*v)->particles_in_size()==0) num_noparent_vtxs_filt++;
        if ((*v)->particles_out_size()==0) num_nochild_vtxs_filt++;
      }
 
      // Write out the change in the number of particles
      ATH_MSG_INFO("Particles filtered: " << num_particles_orig << " -> " << num_particles_filt);
      // Warn if the numbers of "strange" vertices have changed
      if (num_orphan_vtxs_filt != num_orphan_vtxs_orig)
        ATH_MSG_WARNING("Change in orphaned vertices: " << num_orphan_vtxs_orig << " -> " << num_orphan_vtxs_filt);
      if (num_noparent_vtxs_filt != num_noparent_vtxs_orig)
        ATH_MSG_WARNING("Change in no-parent vertices: " << num_noparent_vtxs_orig << " -> " << num_noparent_vtxs_filt);
      if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
        ATH_MSG_WARNING("Change in no-parent vertices: " << num_nochild_vtxs_orig << " -> " << num_nochild_vtxs_filt);
    }
#endif
 
  } // End of the loop over events in the MC event collection
  return StatusCode::SUCCESS;
}

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.

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

finalize()

StatusCode FixHepMC::finalize

(

)

Definition at line 592 of file FixHepMC.cxx.

                              {
  if (m_killLoops  ) ATH_MSG_INFO( "Removed " <<   m_loopKilled << " of " << m_totalSeen << " particles because of loops." );
  if (m_killPDG0   ) ATH_MSG_INFO( "Removed " <<   m_pdg0Killed << " of " << m_totalSeen << " particles because of PDG ID 0." );
  if (m_cleanDecays) ATH_MSG_INFO( "Removed " << m_decayCleaned << " of " << m_totalSeen << " particles while cleaning decay chains." );
  if(m_purgeUnstableWithoutEndVtx) ATH_MSG_INFO( "Removed " << m_unstablePurged << " of " << m_totalSeen << " unstable particles because they had no decay vertex." );
  if (!m_pidmap.empty()) {
    ATH_MSG_INFO( "Replaced " << m_replacedPIDs << " PIDs of particles." );
    for (const auto& p: m_replacedpid_counts) ATH_MSG_INFO( "Replaced " << p.first << " PIDs " << p.second << "times." );
  }  
  return StatusCode::SUCCESS;
}

fromDecay()

bool FixHepMC::fromDecay ( const HepMC::ConstGenParticlePtr & p, std::shared_ptr< std::set< int > > & storage ) const

private

Definition at line 614 of file FixHepMC.cxx.

                                                                                                      {
      if (!p) return false;
      auto v=p->production_vertex();
      if (!v) return false;
#ifdef HEPMC3
      for ( const auto& anc: v->particles_in()) {
        if (MC::isDecayed(anc) && (MC::isTau(anc->pdg_id()) || MC::isHadron(anc->pdg_id()))) return true;
      }
      if (storage->find(p->id()) != storage->end()) return false;
      storage->insert(p->id());
      for ( const auto& anc: v->particles_in()) {
            if (fromDecay(anc, storage)) return true;
      }
#else
      for (auto  anc=v->particles_in_const_begin(); anc != v->particles_in_const_end(); ++anc) {
        if (MC::isDecayed((*anc)) && (MC::isTau((*anc)->pdg_id()) || MC::isHadron((*anc)->pdg_id()))) return true;
      }
      if (storage->find(p->barcode()) != storage->end()) return false;
      storage->insert(p->barcode());
      for (auto  anc=v->particles_in_const_begin(); anc != v->particles_in_const_end(); ++anc){
        if (fromDecay(*anc, storage)) return true;
      }
#endif
      return false;
}

initialize()

StatusCode GenBase::initialize ( )

overridevirtualinherited

Reimplemented in CopyEventWeight, CountHepMC, DumpMC, EvtInclusiveDecay, FillFilterValues, GenFilter, GenModule, HepMCReadFromFile, PrintHijingPars, PrintMC, TestHepMC, and WriteHepMC.

Definition at line 17 of file GenBase.cxx.

                               {
  ATH_CHECK( m_mcevents_const.initialize() );
  m_mcEventKey = m_mcevents_const.key();
 
  // Get the particle property service
  ATH_CHECK(m_ppSvc.retrieve());
  return StatusCode::SUCCESS;
}

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.

isNonTransportableInDecayChain()

bool FixHepMC::isNonTransportableInDecayChain ( const HepMC::ConstGenParticlePtr & p ) const

private

Definition at line 641 of file FixHepMC.cxx.

                                                                                     {
  auto storage = std::make_shared<std::set<int>>();
  return !MC::isTransportable(p->pdg_id()) && fromDecay(p, storage);
}

isPID0()

bool FixHepMC::isPID0 ( const HepMC::ConstGenParticlePtr & p ) const

private

Definition at line 609 of file FixHepMC.cxx.

                                                             {
  return p->pdg_id() == 0;
}

isSimpleLoop()

bool FixHepMC::isSimpleLoop ( const HepMC::ConstGenParticlePtr & p ) const

private

Definition at line 647 of file FixHepMC.cxx.

                                                                   {
  return (p->production_vertex() == p->end_vertex() && p->end_vertex() != nullptr);
}

msg()

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

inlineinherited

Definition at line 24 of file AthCommonMsg.h.

                                {
    return this->msgStream();
  }

msgLvl()

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

inlineinherited

Definition at line 30 of file AthCommonMsg.h.

                                               {
    return this->msgLevel(lvl);
  }

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.

                                                        {
    return pdt().particle(HepPDT::ParticleID(std::abs(pid)));
  }

particleTable()

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

inlineinherited

Get a particle data table.

Definition at line 118 of file GenBase.h.

                                                       {
    return *(m_ppSvc->PDT());
  }

partPropSvc()

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

inlineinherited

Access the particle property service.

Definition at line 113 of file GenBase.h.

                                                        {
    return m_ppSvc;
  }

pdt()

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

inlineinherited

Shorter alias to get a particle data table.

Definition at line 123 of file GenBase.h.

{ return particleTable(); }

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.

  {
    h.renounce();
    PBASE::renounce (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.

                                                          {
    handlesArray.renounce();
  }

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.

                                       {
  StatusCode sc=AthCommonDataStore<AthCommonMsg<Algorithm>>::sysInitialize();
 
  if (sc.isFailure()) {
    return sc;
  }
  ServiceHandle<ICondSvc> cs("CondSvc",name());
  for (auto h : outputHandles()) {
    if (h->isCondition() && h->mode() == Gaudi::DataHandle::Writer) {
      // do this inside the loop so we don't create the CondSvc until needed
      if ( cs.retrieve().isFailure() ) {
        ATH_MSG_WARNING("no CondSvc found: won't autoreg WriteCondHandles");
        return StatusCode::SUCCESS;
      }
      if (cs->regHandle(this,*h).isFailure()) {
        sc = StatusCode::FAILURE;
        ATH_MSG_ERROR("unable to register WriteCondHandle " << h->fullKey()
                      << " with CondSvc");
      }
    }
  }
  return sc;
}

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.

updateVHKA()

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

inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

                                               {
    // debug() << "updateVHKA for property " << p.name() << " " << p.toString() 
    //         << "  size: " << m_vhka.size() << endmsg;
    for (auto &a : m_vhka) {
      std::vector<SG::VarHandleKey*> keys = a->keys();
      for (auto k : keys) {
        k->setOwner(this);
      }
    }
  }

Member Data Documentation

m_cleanDecays

bool FixHepMC::m_cleanDecays

private

Definition at line 53 of file FixHepMC.h.

m_decayCleaned

long FixHepMC::m_decayCleaned

private

Definition at line 66 of file FixHepMC.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_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_forced_length

std::string FixHepMC::m_forced_length {""}

private

Definition at line 57 of file FixHepMC.h.

{""}; // Force length unit for the event

m_forced_momentum

std::string FixHepMC::m_forced_momentum {""}

private

Definition at line 56 of file FixHepMC.h.

{""}; // Force momentum unit for the event

m_ignoreSemiDisconnected

bool FixHepMC::m_ignoreSemiDisconnected

private

Definition at line 55 of file FixHepMC.h.

m_killLoops

bool FixHepMC::m_killLoops

private

Definition at line 51 of file FixHepMC.h.

m_killPDG0

bool FixHepMC::m_killPDG0

private

Definition at line 52 of file FixHepMC.h.

m_looper

MC::Loops<HepMC::GenEvent,HepMC::ConstGenParticlePtr,HepMC::ConstGenVertexPtr> FixHepMC::m_looper

private

member to detect loops

Definition at line 76 of file FixHepMC.h.

m_loopKilled

long FixHepMC::m_loopKilled

private

Definition at line 64 of file FixHepMC.h.

m_mcEventKey

std::string GenBase::m_mcEventKey {}

protectedinherited

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

Definition at line 136 of file GenBase.h.

{};

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 147 of file GenBase.h.

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

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 138 of file GenBase.h.

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

m_pdg0Killed

long FixHepMC::m_pdg0Killed

private

Definition at line 65 of file FixHepMC.h.

m_pidmap

std::map<int,int> FixHepMC::m_pidmap

private

map of pids to change.

Definition at line 72 of file FixHepMC.h.

m_ppSvc

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

privateinherited

Handle on the particle property service.

Definition at line 144 of file GenBase.h.

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

m_purgeUnstableWithoutEndVtx

bool FixHepMC::m_purgeUnstableWithoutEndVtx

private

Definition at line 54 of file FixHepMC.h.

m_replacedpid_counts

std::map<int,int> FixHepMC::m_replacedpid_counts

private

map of counters of replacements.

Definition at line 73 of file FixHepMC.h.

m_replacedPIDs

long FixHepMC::m_replacedPIDs

private

Definition at line 69 of file FixHepMC.h.

m_setHasCycles

bool FixHepMC::m_setHasCycles

private

Definition at line 59 of file FixHepMC.h.

m_totalSeen

long FixHepMC::m_totalSeen

private

Definition at line 68 of file FixHepMC.h.

m_unitsFix

bool FixHepMC::m_unitsFix

private

Definition at line 58 of file FixHepMC.h.

m_unstablePurged

long FixHepMC::m_unstablePurged

private

Definition at line 67 of file FixHepMC.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.

---

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

• FixHepMC.h
• FixHepMC.cxx