ATLAS Offline Software
Public Member Functions | Protected Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
DerivationFramework::CompactHardTruth Class Reference

#include <CompactHardTruth.h>

Inheritance diagram for DerivationFramework::CompactHardTruth:
Collaboration diagram for DerivationFramework::CompactHardTruth:

Public Member Functions

 CompactHardTruth (const std::string &name, ISvcLocator *pSvcLocator)
 Constructor with parameters: More...
 
virtual ~CompactHardTruth ()
 Destructor: More...
 
virtual StatusCode initialize ()
 
virtual StatusCode execute ()
 
virtual StatusCode finalize ()
 
virtual HepMC::FourVector vtxInMom (HepMC::ConstGenVertexPtr)
 
virtual HepMC::FourVector vtxOutMom (HepMC::ConstGenVertexPtr)
 
virtual StatusCode sysInitialize () override
 Override sysInitialize. More...
 
virtual const DataObjIDColl & extraOutputDeps () const override
 Return the list of extra output dependencies. More...
 
ServiceHandle< StoreGateSvc > & evtStore ()
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc. More...
 
const ServiceHandle< StoreGateSvc > & evtStore () const
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc. More...
 
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc. More...
 
virtual StatusCode sysStart () override
 Handle START transition. More...
 
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles. More...
 
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles. More...
 
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T, V, H > &t)
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, SG::VarHandleKey &hndl, const std::string &doc, const SG::VarHandleKeyType &)
 Declare a new Gaudi property. More...
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, SG::VarHandleBase &hndl, const std::string &doc, const SG::VarHandleType &)
 Declare a new Gaudi property. More...
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, SG::VarHandleKeyArray &hndArr, const std::string &doc, const SG::VarHandleKeyArrayType &)
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, T &property, const std::string &doc, const SG::NotHandleType &)
 Declare a new Gaudi property. More...
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, T &property, const std::string &doc="none")
 Declare a new Gaudi property. More...
 
void updateVHKA (Gaudi::Details::PropertyBase &)
 
MsgStream & msg () const
 
MsgStream & msg (const MSG::Level lvl) const
 
bool msgLvl (const MSG::Level lvl) const
 

Protected Member Functions

void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution More...
 
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. More...
 

Private Types

typedef ServiceHandle< StoreGateSvcStoreGateSvc_t
 

Private Member Functions

 CompactHardTruth ()
 Default constructor: More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyArrayType &)
 specialization for handling Gaudi::Property<SG::VarHandleKeyArray> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleType &)
 specialization for handling Gaudi::Property<SG::VarHandleBase> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &t, const SG::NotHandleType &)
 specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray> More...
 

Private Attributes

std::string m_mcEventsName
 Containers. More...
 
std::string m_thinnedMcEventsName
 
int m_evtCount = 0
 
int m_missCount = 0
 
float m_partonCut
 
float m_hardCut
 
int m_dangleFound = 0
 
int m_dangleRemoved = 0
 
float m_danglePtMax = 0.0F
 
float m_danglePtCut
 
int m_maxCount
 
int m_thinParticles = 0
 
int m_thinVertices = 0
 
DataObjIDColl m_extendedExtraObjects
 
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default) More...
 
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default) More...
 
std::vector< SG::VarHandleKeyArray * > m_vhka
 
bool m_varHandleArraysDeclared
 

Detailed Description

Definition at line 46 of file CompactHardTruth.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

◆ CompactHardTruth() [1/2]

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

Constructor with parameters:

Definition at line 41 of file CompactHardTruth.cxx.

42  : ::AthAlgorithm(name, pSvcLocator)
43  , m_mcEventsName("GEN_AOD")
44  , m_thinnedMcEventsName("GEN_AODTHIN")
45  , m_partonCut(10000.)
46  , m_hardCut(10000.)
47  , m_danglePtCut(0.)
48  , m_maxCount(0) {
49  //
50  // Property declaration
51  //
52  declareProperty("McEvent", m_mcEventsName, "input McEventCollection container name");
53  declareProperty("McEventOut", m_thinnedMcEventsName, "output McEventCollection container name");
54  declareProperty("ShowerMassCut", m_partonCut, "mass cut for thinning parton shower");
55  declareProperty("SoftMtCut", m_hardCut, "mt cut for underlying event showers");
56  declareProperty("DanglePtCut", m_danglePtCut, "maximum pt for dangling partons");
57 
58  declareProperty("MaxCount", m_maxCount, "maximum number of events to print");
59 }

◆ ~CompactHardTruth()

DerivationFramework::CompactHardTruth::~CompactHardTruth ( )
virtual

Destructor:

Definition at line 63 of file CompactHardTruth.cxx.

63 {}

◆ CompactHardTruth() [2/2]

DerivationFramework::CompactHardTruth::CompactHardTruth ( )
private

Default constructor:

Member Function Documentation

◆ declareGaudiProperty() [1/4]

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

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

Definition at line 170 of file AthCommonDataStore.h.

172  {
173  return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
174  hndl.value(),
175  hndl.documentation());
176 
177  }

◆ declareGaudiProperty() [2/4]

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  {
159  return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
160  hndl.value(),
161  hndl.documentation());
162 
163  }

◆ declareGaudiProperty() [3/4]

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

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

Definition at line 184 of file AthCommonDataStore.h.

186  {
187  return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
188  hndl.value(),
189  hndl.documentation());
190  }

◆ declareGaudiProperty() [4/4]

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

specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray>

Definition at line 199 of file AthCommonDataStore.h.

200  {
201  return PBASE::declareProperty(t);
202  }

◆ declareProperty() [1/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty ( const std::string &  name,
SG::VarHandleBase hndl,
const std::string &  doc,
const SG::VarHandleType  
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
hndlObject holding the property value.
docDocumentation string for the property.

This is the version for types that derive from SG::VarHandleBase. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 245 of file AthCommonDataStore.h.

249  {
250  this->declare(hndl.vhKey());
251  hndl.vhKey().setOwner(this);
252 
253  return PBASE::declareProperty(name,hndl,doc);
254  }

◆ declareProperty() [2/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty ( const std::string &  name,
SG::VarHandleKey hndl,
const std::string &  doc,
const SG::VarHandleKeyType  
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
hndlObject holding the property value.
docDocumentation string for the property.

This is the version for types that derive from SG::VarHandleKey. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 221 of file AthCommonDataStore.h.

225  {
226  this->declare(hndl);
227  hndl.setOwner(this);
228 
229  return PBASE::declareProperty(name,hndl,doc);
230  }

◆ declareProperty() [3/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty ( const std::string &  name,
SG::VarHandleKeyArray hndArr,
const std::string &  doc,
const SG::VarHandleKeyArrayType  
)
inlineinherited

Definition at line 259 of file AthCommonDataStore.h.

263  {
264 
265  // std::ostringstream ost;
266  // ost << Algorithm::name() << " VHKA declareProp: " << name
267  // << " size: " << hndArr.keys().size()
268  // << " mode: " << hndArr.mode()
269  // << " vhka size: " << m_vhka.size()
270  // << "\n";
271  // debug() << ost.str() << endmsg;
272 
273  hndArr.setOwner(this);
274  m_vhka.push_back(&hndArr);
275 
276  Gaudi::Details::PropertyBase* p = PBASE::declareProperty(name, hndArr, doc);
277  if (p != 0) {
278  p->declareUpdateHandler(&AthCommonDataStore<PBASE>::updateVHKA, this);
279  } else {
280  ATH_MSG_ERROR("unable to call declareProperty on VarHandleKeyArray "
281  << name);
282  }
283 
284  return p;
285 
286  }

◆ declareProperty() [4/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty ( const std::string &  name,
T &  property,
const std::string &  doc,
const SG::NotHandleType  
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
propertyObject holding the property value.
docDocumentation string for the property.

This is the generic version, for types that do not derive from SG::VarHandleKey. It just forwards to the base class version of declareProperty.

Definition at line 333 of file AthCommonDataStore.h.

337  {
338  return PBASE::declareProperty(name, property, doc);
339  }

◆ declareProperty() [5/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty ( const std::string &  name,
T &  property,
const std::string &  doc = "none" 
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
propertyObject holding the property value.
docDocumentation string for the property.

This dispatches to either the generic declareProperty or the one for VarHandle/Key/KeyArray.

Definition at line 352 of file AthCommonDataStore.h.

355  {
356  typedef typename SG::HandleClassifier<T>::type htype;
357  return declareProperty (name, property, doc, htype());
358  }

◆ declareProperty() [6/6]

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  }

◆ 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.

95 { return m_detStore; }

◆ evtStore() [1/2]

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.

85 { return m_evtStore; }

◆ evtStore() [2/2]

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

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

Definition at line 90 of file AthCommonDataStore.h.

90 { return m_evtStore; }

◆ execute()

StatusCode DerivationFramework::CompactHardTruth::execute ( )
virtual

float eold = pp->momentum().e();

Definition at line 107 of file CompactHardTruth.cxx.

107  {
108 
109  ++m_evtCount;
110  // if( m_evtCount%100 == 0 ){
111  ATH_MSG_INFO("Executing " << name() << " " << m_evtCount);
112  //}
113 
114  // Normally doPrint is used to print the first m_maxCount events
115  // before and after thinning.
116  // doExtra adds extra intermediate event printouts.
117  // doDebug allows debug printout for a range of events.
118  bool doPrint = m_evtCount < m_maxCount;
119  bool doDebug = false;
120  bool doExtra = false;
121  // doDebug = doPrint;
122  // doExtra = doPrint;
123 
124  // Retrieve input data
125  const McEventCollection* mcEvts = nullptr;
126  if (!evtStore()->retrieve(mcEvts, m_mcEventsName).isSuccess() || nullptr == mcEvts) {
127  ATH_MSG_WARNING("could not retrieve mc collection at [" << m_mcEventsName << "]!");
128  return StatusCode::FAILURE;
129  }
130 
131  if (mcEvts->empty()) {
132  ATH_MSG_WARNING("empty McEventCollection at [" << m_mcEventsName << "]");
133  return StatusCode::SUCCESS;
134  }
135 
136  // Create output collection
137  McEventCollection* thinnedMcEvts = new McEventCollection;
138  if (!evtStore()->record(thinnedMcEvts, m_thinnedMcEventsName).isSuccess()) {
139  ATH_MSG_WARNING("Could not record thinned mc collection at [" << m_thinnedMcEventsName << "]!");
140  delete thinnedMcEvts;
141  thinnedMcEvts = nullptr;
142  return StatusCode::FAILURE;
143  }
144  if (evtStore()->setConst(thinnedMcEvts).isFailure()) { ATH_MSG_WARNING("Could not lock the McEventCollection at [" << m_thinnedMcEventsName << "] !!"); }
145 
146  // Signal event is first (only?) event; front() is from DataVector
147  const HepMC::GenEvent* mcEvt = mcEvts->front();
148  auto wtCont = mcEvt->weights();
149  if (!wtCont.empty()) {
150  } else {
151  ATH_MSG_WARNING("Weights not found for mc collection [" << m_mcEventsName << "]");
152  }
153  int inEvent = mcEvt->event_number();
154  if (doDebug) ATH_MSG_DEBUG("FETCHED count/event " << m_evtCount << " " << inEvent);
155 
157  // New event - copy of original
159 
160  HepMC::GenEvent* thinEvt = new HepMC::GenEvent(*mcEvt);
161  int nEvent = thinEvt->event_number();
162  if (doPrint) ATH_MSG_DEBUG("New event number = " << nEvent);
163 
164  if (doPrint) {
165  std::cout << "========== BEGIN EVENT BEFORE THINNING ==========" << std::endl;
166  HepMC::Print::line(std::cout, *thinEvt);
167  std::cout << "========== END EVENT BEFORE THINNING ==========" << std::endl;
168  }
169 
171  // Containers for manipulating particles/vertices
173 
174  // Hadronization vertex must have at least two incoming partons to
175  // conserve color. All outgoing must be hadrons; c.f. Pythia 2->1
176  // color reconnection vertices.
177  // Const iterators => must(?) save changes, then do them.
178  // Always remove particles from vertices, then delete.
179  // removePV: remove particle from vertex
180  // addinPV: add particle in to vertex
181  // addoutPV: add particle out from vertex
182  // removeV: remove vertex from event
183  // deleteP: delete particle after all passes
184  // deleteV: delete vertex after all passes
185  // HepMC ~GenVertex deletes particles, so remove them from ALL vertices
186 
187  typedef std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr> vpPair;
188  std::vector<vpPair> removePV;
189  std::vector<vpPair> addinPV;
190  std::vector<vpPair> addoutPV;
191  std::vector<HepMC::GenVertexPtr> removeV;
192 
193  // Use list with sort/unique to insure no multiple deletion
194  std::list<HepMC::GenParticlePtr> deleteP;
195  std::list<HepMC::GenVertexPtr> deleteV;
200 
202  // Find hadronization vertices
204 
205  if (doDebug) ATH_MSG_DEBUG("Find hadronization vertices");
206 
207  std::vector<HepMC::GenVertexPtr> hadVertices;
208 
209 #ifdef HEPMC3
210  for (auto& hadv: thinEvt->vertices()) {
211  if (!hadv) continue;
212  if (hadv->particles_in().size() < 2) continue;
213  if (hadv->particles_out().size() < 1) continue;
214  // Check hadronization vertex
215  // isHad is true if at least one hadron
216  // q qbar -> pi is allowed, but q qbar -> W... is not
217  bool isHadVtx = true;
218  bool isHadOut = false;
219  for (const auto& inp: hadv->particles_in() ) {
220  if (!(MC::isParton(inp)|| MC::isDiquark(inp)) ) { isHadVtx = false; break;}
221  }
222  for (const auto& vp: hadv->particles_out()) {
223  if (MC::isParton(vp)|| MC::isDiquark(vp)) isHadVtx = false;
224  if (MC::isHadron(vp)) isHadOut = true;
225  }
226  isHadVtx = isHadVtx && isHadOut;
227  if (isHadVtx) hadVertices.push_back(hadv);
228  if (doDebug && isHadVtx) ATH_MSG_VERBOSE("Hadronization vertex " << hadv);
229  }
230 #else
231  HepMC::GenEvent::vertex_iterator hadv = thinEvt->vertices_begin();
232  HepMC::GenEvent::vertex_iterator hadvB = thinEvt->vertices_begin();
233  HepMC::GenEvent::vertex_iterator hadvE = thinEvt->vertices_end();
234 
235  for (; hadv != hadvE; ++hadv) {
236  if (!(*hadv)) continue;
237  if ((*hadv)->particles_in_size() < 2) continue;
238  if ((*hadv)->particles_out_size() < 1) continue;
239 
240  // Check hadronization vertex
241  // isHad is true if at least one hadron
242  // q qbar -> pi is allowed, but q qbar -> W... is not
243  bool isHadVtx = true;
244  bool isHadOut = false;
245  HepMC::GenVertex::particles_in_const_iterator inp = (*hadv)->particles_in_const_begin();
246  HepMC::GenVertex::particles_in_const_iterator inpE = (*hadv)->particles_in_const_end();
247  for (; inp != inpE; ++inp) {
248  if (!(MC::isParton(*inp)|| MC::isDiquark(*inp))) { isHadVtx = false; break;}
249  }
250  HepMC::GenVertex::particles_out_const_iterator vp = (*hadv)->particles_out_const_begin();
251  HepMC::GenVertex::particles_out_const_iterator vpE = (*hadv)->particles_out_const_end();
252  for (; vp != vpE; ++vp) {
253  if (MC::isParton(*vp)|| MC::isDiquark(*vp)) isHadVtx = false;
254  if (MC::isHadron(*vp)) isHadOut = true;
255  }
256  isHadVtx = isHadVtx && isHadOut;
257  if (isHadVtx) hadVertices.push_back(*hadv);
258  if (doDebug && isHadVtx) ATH_MSG_VERBOSE("Hadronization vertex " << *hadv);
259  }
260 #endif
261 
262  if (hadVertices.empty()) {
263  ATH_MSG_WARNING("No hadronization vertices for event " << nEvent);
264  ATH_MSG_WARNING("Exiting without changing event.");
265  thinnedMcEvts->push_back(thinEvt);
266  return StatusCode::SUCCESS;
267  }
268 
270  // Remove all incoming partons from hadronization vertices
271  // Remove and delete all descendants
273 
274 #ifdef HEPMC3
275  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
276  HepMC::GenVertexPtr ivtx = hadVertices[iv];
277  if (doDebug) ATH_MSG_DEBUG("Removing partons from hadVertex " << ivtx);
278  for (auto pin: ivtx->particles_in()) {
279  removePV.push_back(vpPair(ivtx, pin));
280  }
281  }
282  // Remove all descendant particles. Will remove empty vertices later.
283  // Might have parton decays of hadrons - hence delete sort/unique
284  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
285  HepMC::GenVertexPtr ivtx = hadVertices[iv];
286  auto descendants = HepMC::descendant_particles(ivtx);
287  for (auto pout: descendants) {
288  HepMC::GenVertexPtr vpar = pout->production_vertex();
289  if (vpar) removePV.push_back(vpPair(vpar, pout));
290  HepMC::GenVertexPtr vend = pout->end_vertex();
291  if (vend) removePV.push_back(vpPair(vend, pout));
292  deleteP.push_back(pout);
293  }
294  }
295 #else
296  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
297  HepMC::GenVertex* ivtx = hadVertices[iv];
298  if (doDebug) ATH_MSG_DEBUG("Removing partons from hadVertex " << ivtx);
299  HepMC::GenVertex::particles_in_const_iterator pin = ivtx->particles_in_const_begin();
300  HepMC::GenVertex::particles_in_const_iterator pinE = ivtx->particles_in_const_end();
301  for (; pin != pinE; ++pin) {
302  removePV.emplace_back(ivtx, *pin);
303  }
304  }
305 
306  // Remove all descendant particles. Will remove empty vertices later.
307  // Might have parton decays of hadrons - hence delete sort/unique
308  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
309  HepMC::GenVertex* ivtx = hadVertices[iv];
310  HepMC::GenVertex::particle_iterator pout = ivtx->particles_begin(HepMC::descendants);
311  HepMC::GenVertex::particle_iterator poutE = ivtx->particles_end(HepMC::descendants);
312  for (; pout != poutE; ++pout) {
313  HepMC::GenVertex* vpar = (*pout)->production_vertex();
314  if (vpar) removePV.emplace_back(vpar, *pout);
315  HepMC::GenVertex* vend = (*pout)->end_vertex();
316  if (vend) removePV.emplace_back(vend, *pout);
317  deleteP.push_back(*pout);
318  }
319  }
320 #endif
321 
322  // Remove empty vertices
323  // Remove all particles from Geant vertices
324  // Remove and delete Geant vertices and particles
325  // All Geant particles should have Geant parent vertex
326 
327 #ifdef HEPMC3
328  for (auto hadv:thinEvt->vertices()) {
329  // Empth vertices
330  if (hadv->particles_in().size() == 0 && hadv->particles_out().size() == 0) {
331  removeV.push_back(hadv);
332  deleteV.push_back(hadv);
333  }
334  // Geant vertices/particles
335  if (!HepMC::is_simulation_vertex(hadv)) continue;
336  for (auto pin: hadv->particles_in()) {
337  removePV.push_back(vpPair(hadv, pin));
338  if (HepMC::is_simulation_particle(pin)) { deleteP.push_back(pin); }
339  }
340  for (auto pout: hadv->particles_out()) {
341  removePV.push_back(vpPair(hadv, pout));
342  if (HepMC::is_simulation_particle(pout)) { deleteP.push_back(pout); }
343  }
344  removeV.push_back(hadv);
345  deleteV.push_back(hadv);
346  }
347 #else
348  for (hadv = hadvB; hadv != hadvE; ++hadv) {
349 
350  // Empth vertices
351  if ((*hadv)->particles_in_size() == 0 && (*hadv)->particles_out_size() == 0) {
352  removeV.push_back(*hadv);
353  deleteV.push_back(*hadv);
354  }
355 
356  // Geant vertices/particles
357  if (!HepMC::is_simulation_vertex(*hadv)) continue;
358  HepMC::GenVertex::particles_in_const_iterator pin = (*hadv)->particles_in_const_begin();
359  HepMC::GenVertex::particles_in_const_iterator pinE = (*hadv)->particles_in_const_end();
360  for (; pin != pinE; ++pin) {
361  removePV.emplace_back(*hadv, *pin);
362  if (HepMC::is_simulation_particle(*pin)) { deleteP.push_back(*pin); }
363  }
364  HepMC::GenVertex::particles_out_const_iterator pout = (*hadv)->particles_out_const_begin();
365  HepMC::GenVertex::particles_out_const_iterator poutE = (*hadv)->particles_out_const_end();
366  for (; pout != poutE; ++pout) {
367  removePV.emplace_back(*hadv, *pout);
368  if (HepMC::is_simulation_particle(*pout)) { deleteP.push_back(*pout); }
369  }
370  removeV.push_back(*hadv);
371  deleteV.push_back(*hadv);
372  }
373 #endif
374 
375  // Actually implement changes
376 
377  for (unsigned int i = 0; i < removePV.size(); ++i) {
378  HepMC::GenVertexPtr v = removePV[i].first;
379  HepMC::GenParticlePtr p = removePV[i].second;
380 #ifdef HEPMC3
381  v->remove_particle_in(p);
382  v->remove_particle_out(p);
383 #else
384  v->remove_particle(p);
385 #endif
386  }
387 
388  for (unsigned int i = 0; i < addoutPV.size(); ++i) {
389  HepMC::GenVertexPtr v = addoutPV[i].first;
390  HepMC::GenParticlePtr p = addoutPV[i].second;
391  v->add_particle_out(p);
392  }
393 #ifdef HEPMC3
394  for (unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
395  HepMC::GenVertexPtr v = hadVertices[iv];
396  if (v->particles_in().size() != 0 || v->particles_out().size() != 0) {
397  ATH_MSG_WARNING("Removing vertex " << v << " for event " << nEvent
398  << " with in/out particles " << v->particles_in().size() << " " << v->particles_out().size());
399  }
400  thinEvt->remove_vertex(hadVertices[iv]);
401  }
402 //AV: HepMC3 uses smart pointers
403 #else
404 
405  for (unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
406  HepMC::GenVertex* v = hadVertices[iv];
407  if (v->particles_in_size() != 0 || v->particles_out_size() != 0) {
408  ATH_MSG_WARNING("Removing vertex " << HepMC::uniqueID(v) << " for event " << nEvent << " with in/out particles " << v->particles_in_size() << " " << v->particles_out_size());
409  }
410  if (!thinEvt->remove_vertex(hadVertices[iv])) { ATH_MSG_WARNING("Error removing vertex " << HepMC::uniqueID(v) << " for event " << nEvent); }
411  }
412 
413  // Delete removed particles/vertices
414 
415  if (doDebug) ATH_MSG_DEBUG("Deleting hadronization vertices " << deleteV.size());
416  deleteV.sort();
417  deleteV.unique();
418  for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
419  if (doDebug) ATH_MSG_VERBOSE("Deleting vertex " << HepMC::uniqueID(*dvItr));
420  if (*dvItr) delete (*dvItr);
421  }
422 
423  deleteP.sort();
424  deleteP.unique();
425  for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
426  if (doDebug) ATH_MSG_VERBOSE("Deleting particle " << HepMC::uniqueID(*dpItr));
427  if (*dpItr) delete (*dpItr);
428  }
429 #endif
430 
432  // Cluster final partons
434 
435  if (doDebug && doExtra) {
436  std::cout << "========== BEGIN EVENT BEFORE CLUSTER ==========" << std::endl;
437  HepMC::Print::line(std::cout, *thinEvt);
438  std::cout << "========== END EVENT BEFORE CLUSTER ==========" << std::endl;
439  }
440 
441  // Possible cases:
442  // 1->1 branching: Drop earlier 1 and vertex.
443  // 2->1 connection: Drop 1 and vertex; leave 2 free
444  // Note momentum is conserved for these.
445  // Do not split 2->1 from initial 2->1->2 in Herwig.
446  // 1->2 branching: Set p1 = p2a+p2b, drop 2 and vertex.
447  // Add 1 to first vertex.
448 
449  // Preserving color connections required merging vertices. Now no
450  // longer needed.
451 
452  bool moreP = true;
453  using vpPair = std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr>;
454  removePV.clear();
455  addinPV.clear();
456  addoutPV.clear();
457  removeV.clear();
458  deleteP.clear();
459  deleteV.clear();
460 
461  using pkPair = std::pair<HepMC::GenParticlePtr, HepMC::FourVector>;
462  std::vector<pkPair> changePK;
463 
464  if (doDebug) ATH_MSG_DEBUG("Start parton thinning");
465  while (moreP) {
466  if (doDebug) ATH_MSG_DEBUG("New parton pass " << inEvent << " " << thinEvt->particles_size() << " " << thinEvt->vertices_size());
467 
468  moreP = false;
469  removePV.clear();
470  addinPV.clear();
471  addoutPV.clear();
472  removeV.clear();
473  changePK.clear();
474 #ifdef HEPMC3
475  // Find final partons
476  for (auto fp: thinEvt->particles() ) {
477  int iCase = 0;
478  if (!((MC::isParton(fp)|| MC::isDiquark(fp)) && fp->end_vertex() == nullptr)) continue;
479  if (doDebug) ATH_MSG_DEBUG("Starting final parton " << fp);
480  // Production/end vertices
481  HepMC::GenVertexPtr pvtx = fp->production_vertex();
482  if (!pvtx) {
483  ATH_MSG_WARNING("Missing production for final parton " << fp);
484  continue;
485  }
486  if (doDebug) ATH_MSG_DEBUG("Final parton " << pvtx << " " << fp);
488  // Case 1->1
490  // One-particle decay; use final particle
491  // ppvtx -> pp -> pvtx -> fp
492  if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 1) {
493  // Incoming particle to parent vertex
494  HepMC::GenParticlePtr pp = pvtx->particles_in().front();
495  if (!pp || pp->parent_event() == nullptr) {
496  ATH_MSG_DEBUG("1->1: missing pp for fp " << fp);
497  ++m_missCount;
498  continue;
499  }
500  // Its parent vertex
501  HepMC::GenVertexPtr ppvtx = pp->production_vertex();
502  if (!ppvtx || ppvtx->parent_event() == nullptr) {
503  ATH_MSG_DEBUG("1->1: missing ppvtx for fp " << fp);
504  ++m_missCount;
505  continue;
506  }
507  moreP = true;
508  iCase = 1;
509  removePV.push_back(vpPair(ppvtx, pp));
510  removePV.push_back(vpPair(pvtx, pp));
511  deleteP.push_back(pp);
512  removeV.push_back(pvtx);
513  deleteV.push_back(pvtx);
514  addoutPV.push_back(vpPair(ppvtx, fp));
515  if (doDebug) { ATH_MSG_DEBUG("1->1: ppvtx,pp,pvtx,fp,evtx " << ppvtx << " " << pp << " " << pvtx << " " << fp); }
516  }
518  // Case 2->1
520  // Color recombination. Momentum is conserved so just keep 2.
521  // Drop 1 and vertex.
522  // ppvtx1,ppvtx2 -> pp1,pp2 -> pvtx -> fp
523  // Recombination should not affect hard physics!
524  if (pvtx->particles_in().size() == 2 && pvtx->particles_out().size() == 1) {
525  // Incoming particles to parent vertex
526  HepMC::GenParticlePtr pp1 = pvtx->particles_in().front();
527  HepMC::GenParticlePtr pp2 = pvtx->particles_in().back();
528  // Check for 2->1->2 initial state interactions in Herwig++
529  // Initial partons have pt=0, use pt<0.001MeV
530  if (std::abs(pp1->momentum().perp()) < 1.e-3) continue;
531  if (std::abs(pp2->momentum().perp()) < 1.e-3) continue;
532  // Their parent vertices
533  HepMC::GenVertexPtr ppvtx1 = pp1->production_vertex();
534  HepMC::GenVertexPtr ppvtx2 = pp2->production_vertex();
535  if (!ppvtx1 || ppvtx1->parent_event() == nullptr) {
536  ATH_MSG_DEBUG("2->1: missing ppvtx1 for fp " << fp);
537  ++m_missCount;
538  continue;
539  }
540  if (!ppvtx2 || ppvtx2->parent_event() == nullptr) {
541  ATH_MSG_DEBUG("2->1: missing ppvtx2 for fp " << fp);
542  ++m_missCount;
543  continue;
544  }
545  moreP = true;
546  iCase = 2;
547  removePV.push_back(vpPair(pvtx, fp));
548  removePV.push_back(vpPair(pvtx, pp1));
549  removePV.push_back(vpPair(pvtx, pp2));
550  deleteP.push_back(fp);
551  removeV.push_back(pvtx);
552  deleteV.push_back(pvtx);
553  if (doDebug) {
554  ATH_MSG_DEBUG("2->1: ppvtx1,pp1,ppvtx2,pp2,pvtx,fp " << ppvtx1 << " " << pp1
555  << " " << ppvtx2 << " " << pp2 << " " << pvtx << " "<< fp);
556  }
557  }
559  // Case 1->2
561  // Parton branching. Momentum not conserved; 2 momenta correct
562  // Drop only if mass is below cut
563  // ppvtx -> pp -> pvtx -> pout1,pout2/fp
564  if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 2) {
565  HepMC::GenParticlePtr pout1 = pvtx->particles_out().front();
566  HepMC::GenParticlePtr pout2 = pvtx->particles_out().back();
567  // Require two final partons and avoid duplication
568  if (fp == pout1) {
569  if (!((MC::isParton(pout2)|| MC::isDiquark(pout2)) && pout2->end_vertex() == nullptr)) {
570  if (doDebug) ATH_MSG_DEBUG("1->2: not final " << pout2);
571  continue;
572  }
573  } else if (fp == pout2) {
574  if (!((MC::isParton(pout1)|| MC::isDiquark(pout1)) && pout1->end_vertex() == nullptr)) {
575  if (doDebug) ATH_MSG_DEBUG("1->2: not final " << pout1);
576  continue;
577  }
578  } else {
579  ATH_MSG_WARNING("1->2: No match found for branching " << fp << " " << pvtx << " " << pout1 << " " << pout2);
580  continue;
581  }
582  if (fp != pout1) continue;
583  // Incoming particle
584  HepMC::GenParticlePtr pp = pvtx->particles_in().front();
585  // Do not merge initial partons (pt<1MeV or m<-1MeV)
586  if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0) continue;
587  if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0) continue;
588  // Parton pair mass cut
589  HepMC::FourVector p12 = vtxOutMom(pvtx);
590  double m12 = p12.m();
591  if (m12 < 0) {
592  if (std::abs(m12) > 10. + 1.0e-5 * p12.e()) {
593  ATH_MSG_WARNING("Spacelike mass^2 for parton sum " << m12 << " " << pp << " " << pvtx << " " << pout1 << " " << pout2);
594  }
595  m12 = 0;
596  }
597  if (doDebug) ATH_MSG_DEBUG("1->2: parton pair mass " << m12);
598  // If mass > cut, keep pair
599  if (m12 > m_partonCut) {
600  if (doDebug) ATH_MSG_DEBUG("Keeping 1->2: parton mass " << m12);
601  continue;
602  }
603  // Associated vertices
604  HepMC::GenVertexPtr ppvtx = pp->production_vertex();
605  if (!ppvtx || ppvtx->parent_event() == nullptr) {
606  ATH_MSG_DEBUG("1->2: missing ppvtx for fp " << fp);
607  ++m_missCount;
608  continue;
609  }
610  // Merge branching
611  moreP = true;
612  iCase = 3;
613  if (doDebug) ATH_MSG_DEBUG("Merging 1->2: mass " << p12.m());
614  changePK.push_back(pkPair(pp, p12));
615  removePV.push_back(vpPair(pvtx, pp));
616  removePV.push_back(vpPair(pvtx, pout1));
617  removePV.push_back(vpPair(pvtx, pout2));
618  deleteP.push_back(pout1);
619  deleteP.push_back(pout2);
620  removeV.push_back(pvtx);
621  deleteV.push_back(pvtx);
622  if (doDebug) {
623  ATH_MSG_DEBUG("Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << ppvtx << " " << pp << " " << pvtx << " " << pout1 << " "
624  << pout2);
625  ATH_MSG_DEBUG("Merge 1->2: id " << pp->pdg_id() << " " << pout1->pdg_id() << " " << pout2->pdg_id());
626  }
627  } // end 1->2 case
629  // Incoming proton vertex
631  // Do nothing
632  if (pvtx->particles_in().size() == 1) {
633  // Incoming particle to parent vertex
634  HepMC::GenParticlePtr pp = pvtx->particles_in().front();
635  if (std::abs(pp->pdg_id()) == MC::PROTON) iCase = -1;
636  }
637  // Case not found
638  // Need test for 2->2 in underlying event
639  if (iCase == 0) {
640  if (doDebug) ATH_MSG_DEBUG("Case not found " << pvtx << " " << fp << " " << pvtx->particles_in().size() << " " << pvtx->particles_out().size());
641  }
642  } // end final parton loop
643 #else
644 
645  HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
646  HepMC::GenEvent::particle_iterator finpB = thinEvt->particles_begin();
647  HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();
648 
649  // Find final partons
650  for (finp = finpB; finp != finpE; ++finp) {
651  int iCase = 0;
652 
653  HepMC::GenParticle* fp = *finp;
654  if (!((MC::isParton(fp)|| MC::isDiquark(fp)) && fp->end_vertex() == nullptr)) continue;
655  if (doDebug) ATH_MSG_DEBUG("Starting final parton " << HepMC::uniqueID(fp));
656 
657  // Production/end vertices
658  HepMC::GenVertex* pvtx = fp->production_vertex();
659  if (!pvtx) {
660  ATH_MSG_WARNING("Missing production for final parton " << HepMC::uniqueID(fp));
661  continue;
662  }
663  if (doDebug) ATH_MSG_DEBUG("Final parton " << HepMC::uniqueID(pvtx) << " " << HepMC::uniqueID(fp));
664 
666  // Case 1->1
668 
669  // One-particle decay; use final particle
670  // ppvtx -> pp -> pvtx -> fp
671 
672  if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 1) {
673  // Incoming particle to parent vertex
674  HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
675  HepMC::GenParticle* pp = *pitr;
676  if (!pp || HepMC::uniqueID(pp) == 0) {
677  ATH_MSG_DEBUG("1->1: missing pp for fp " << HepMC::uniqueID(fp));
678  ++m_missCount;
679  continue;
680  }
681  // Its parent vertex
682  HepMC::GenVertex* ppvtx = pp->production_vertex();
683  if (!ppvtx || HepMC::uniqueID(ppvtx) == 0) {
684  ATH_MSG_DEBUG("1->1: missing ppvtx for fp " << HepMC::uniqueID(fp));
685  ++m_missCount;
686  continue;
687  }
688  moreP = true;
689  iCase = 1;
690 
691  removePV.emplace_back(ppvtx, pp);
692  removePV.emplace_back(pvtx, pp);
693  deleteP.push_back(pp);
694  removeV.push_back(pvtx);
695  deleteV.push_back(pvtx);
696  addoutPV.emplace_back(ppvtx, fp);
697  if (doDebug) { ATH_MSG_DEBUG("1->1: ppvtx,pp,pvtx,fp,evtx " << HepMC::uniqueID(ppvtx) << " " << HepMC::uniqueID(pp) << " " << HepMC::uniqueID(pvtx) << " " << HepMC::uniqueID(fp)); }
698  }
699 
701  // Case 2->1
703 
704  // Color recombination. Momentum is conserved so just keep 2.
705  // Drop 1 and vertex.
706  // ppvtx1,ppvtx2 -> pp1,pp2 -> pvtx -> fp
707  // Recombination should not affect hard physics!
708 
709  if (pvtx->particles_in_size() == 2 && pvtx->particles_out_size() == 1) {
710  // Incoming particles to parent vertex
711  HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
712  HepMC::GenParticle* pp1 = *pitr;
713  ++pitr;
714  HepMC::GenParticle* pp2 = *pitr;
715 
716  // Check for 2->1->2 initial state interactions in Herwig++
717  // Initial partons have pt=0, use pt<0.001MeV
718  if (std::abs(pp1->momentum().perp()) < 1.e-3) continue;
719  if (std::abs(pp2->momentum().perp()) < 1.e-3) continue;
720  // Their parent vertices
721  HepMC::GenVertex* ppvtx1 = pp1->production_vertex();
722  HepMC::GenVertex* ppvtx2 = pp2->production_vertex();
723  if (!ppvtx1 || HepMC::uniqueID(ppvtx1) == 0) {
724  ATH_MSG_DEBUG("2->1: missing ppvtx1 for fp " << HepMC::uniqueID(fp));
725  ++m_missCount;
726  continue;
727  }
728  if (!ppvtx2 || HepMC::uniqueID(ppvtx2) == 0) {
729  ATH_MSG_DEBUG("2->1: missing ppvtx2 for fp " << HepMC::uniqueID(fp));
730  ++m_missCount;
731  continue;
732  }
733 
734  moreP = true;
735  iCase = 2;
736 
737  removePV.emplace_back(pvtx, fp);
738  removePV.emplace_back(pvtx, pp1);
739  removePV.emplace_back(pvtx, pp2);
740  deleteP.push_back(fp);
741  removeV.push_back(pvtx);
742  deleteV.push_back(pvtx);
743 
744  if (doDebug) {
745  ATH_MSG_DEBUG("2->1: ppvtx1,pp1,ppvtx2,pp2,pvtx,fp " << HepMC::uniqueID(ppvtx1) << " " << HepMC::uniqueID(pp1) << " " << HepMC::uniqueID(ppvtx2) << " " << HepMC::uniqueID(pp2) << " " << HepMC::uniqueID(pvtx) << " "
746  << HepMC::uniqueID(fp));
747  }
748  }
749 
751  // Case 1->2
753 
754  // Parton branching. Momentum not conserved; 2 momenta correct
755  // Drop only if mass is below cut
756  // ppvtx -> pp -> pvtx -> pout1,pout2/fp
757 
758  if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 2) {
759  HepMC::GenVertex::particles_out_const_iterator poutitr = pvtx->particles_out_const_begin();
760  HepMC::GenParticle* pout1 = *poutitr;
761  ++poutitr;
762  HepMC::GenParticle* pout2 = *poutitr;
763 
764  // Require two final partons and avoid duplication
765  if (fp == pout1) {
766  if (!((MC::isParton(pout2)|| MC::isDiquark(pout2)) && pout2->end_vertex() == nullptr)) {
767  if (doDebug) ATH_MSG_DEBUG("1->2: not final " << HepMC::uniqueID(pout2));
768  continue;
769  }
770  } else if (fp == pout2) {
771  if (!((MC::isParton(pout1)|| MC::isDiquark(pout1)) && pout1->end_vertex() == nullptr)) {
772  if (doDebug) ATH_MSG_DEBUG("1->2: not final " << HepMC::uniqueID(pout1));
773  continue;
774  }
775  } else {
776  ATH_MSG_WARNING("1->2: No match found for branching " << HepMC::uniqueID(fp) << " " << HepMC::uniqueID(pvtx) << " " << HepMC::uniqueID(pout1) << " " << HepMC::uniqueID(pout2));
777  continue;
778  }
779  if (fp != pout1) continue;
780  // Incoming particle
781  HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
782  HepMC::GenParticle* pp = *pitr;
783 
784  // Do not merge initial partons (pt<1MeV or m<-1MeV)
785  if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0) continue;
786  if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0) continue;
787 
788  // Parton pair mass cut
789  HepMC::FourVector p12 = vtxOutMom(pvtx);
790  double m12 = p12.m();
791  if (m12 < 0) {
792  if (fabs(m12) > 10. + 1.0e-5 * p12.e()) {
793  ATH_MSG_WARNING("Spacelike mass^2 for parton sum " << m12 << " " << HepMC::uniqueID(pp) << " " << HepMC::uniqueID(pvtx) << " " << HepMC::uniqueID(pout1) << " " << HepMC::uniqueID(pout2));
794  }
795  m12 = 0;
796  }
797  if (doDebug) ATH_MSG_DEBUG("1->2: parton pair mass " << m12);
798  // If mass > cut, keep pair
799  if (m12 > m_partonCut) {
800  if (doDebug) ATH_MSG_DEBUG("Keeping 1->2: parton mass " << m12);
801  continue;
802  }
803 
804  // Associated vertices
805  HepMC::GenVertex* ppvtx = pp->production_vertex();
806  if (!ppvtx || HepMC::uniqueID(ppvtx) == 0) {
807  ATH_MSG_DEBUG("1->2: missing ppvtx for fp " << HepMC::uniqueID(fp));
808  ++m_missCount;
809  continue;
810  }
811 
812  // Merge branching
813  moreP = true;
814  iCase = 3;
815  if (doDebug) ATH_MSG_DEBUG("Merging 1->2: mass " << p12.m());
816 
817  changePK.emplace_back(pp, p12);
818  removePV.emplace_back(pvtx, pp);
819  removePV.emplace_back(pvtx, pout1);
820  removePV.emplace_back(pvtx, pout2);
821 
822  deleteP.push_back(pout1);
823  deleteP.push_back(pout2);
824  removeV.push_back(pvtx);
825  deleteV.push_back(pvtx);
826 
827  if (doDebug) {
828  ATH_MSG_DEBUG("Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << ppvtx << " " << pp << " " << pvtx << " " << pout1 << " "
829  << pout2);
830  ATH_MSG_DEBUG("Merge 1->2: id " << pp->pdg_id() << " " << pout1->pdg_id() << " " << pout2->pdg_id());
831  }
832  } // end 1->2 case
833 
835  // Incoming proton vertex
837 
838  // Do nothing
839  if (pvtx->particles_in_size() == 1) {
840  // Incoming particle to parent vertex
841  HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
842  HepMC::GenParticle* pp = *pitr;
843  if (std::abs(pp->pdg_id()) == MC::PROTON) iCase = -1;
844  }
845 
846  // Case not found
847  // Need test for 2->2 in underlying event
848  if (iCase == 0) {
849  if (doDebug) ATH_MSG_DEBUG("Case not found " << HepMC::uniqueID(pvtx) << " " << HepMC::uniqueID(fp) << " " << pvtx->particles_in_size() << " " << pvtx->particles_out_size());
850  }
851 
852  } // end final parton loop
853 #endif
854 
855  // Actually implement changes -- remove particles from vertices
856  // Parton ends free, so no addinPV
857  if (doDebug) ATH_MSG_DEBUG("Actually removing particles " << removePV.size());
858 
859  for (unsigned int i = 0; i < removePV.size(); ++i) {
860  HepMC::GenVertexPtr v = removePV[i].first;
861  HepMC::GenParticlePtr p = removePV[i].second;
862  if (doDebug) ATH_MSG_VERBOSE("Removing v,p " << v << " " << p);
863 #ifdef HEPMC3
864  v->remove_particle_in(p);
865  v->remove_particle_out(p);
866 #else
867  v->remove_particle(p);
868 #endif
869  }
870 
871  // Actually implement changes -- add particles to vertices
872  if (doDebug) ATH_MSG_DEBUG("Actually add particles in/out " << addinPV.size() << " " << addoutPV.size());
873  for (unsigned int i = 0; i < addoutPV.size(); ++i) {
874  HepMC::GenVertexPtr v = addoutPV[i].first;
875  HepMC::GenParticlePtr p = addoutPV[i].second;
876  if (doDebug) ATH_MSG_VERBOSE("Adding v,p " << v << " " << p);
877  v->add_particle_out(p);
878  }
879 
880  // Actually implement changes -- change momenta
881  for (unsigned int i = 0; i < changePK.size(); ++i) {
882  HepMC::GenParticlePtr pp = changePK[i].first;
884  pp->set_momentum(changePK[i].second);
885  }
886 
887  // Actually implement changes -- remove vertices
888  if (doDebug) ATH_MSG_DEBUG("Actually remove vertices " << removeV.size());
889  for (unsigned int i = 0; i < removeV.size(); ++i) {
890 #ifdef HEPMC3
891  int nv = thinEvt->vertices().size();
892  if (doDebug) { ATH_MSG_VERBOSE("Removing vertex " << removeV[i] << " " << nv << " " << thinEvt->vertices().size()); }
893  thinEvt->remove_vertex(removeV[i]);
894 #else
895  int nv = thinEvt->vertices_size();
896  if (thinEvt->remove_vertex(removeV[i])) {
897  if (doDebug) { ATH_MSG_VERBOSE("Removed vertex " << removeV[i] << " " << nv << " " << thinEvt->vertices_size()); }
898  } else {
899  ATH_MSG_WARNING("Failed to remove vertex " << removeV[i]);
900  }
901 #endif
902  }
903  if (doDebug) ATH_MSG_DEBUG("End while(moreP) pass " << moreP);
904 
905  } // end moreP
906 
907 #ifdef HEPMC3
908 //AV HepMC3 uses smartpointers. This part is not needed.
909 #else
910  // Delete removed particles/vertices
911  if (doDebug) ATH_MSG_DEBUG("Deleting vertices " << deleteV.size());
912  deleteV.sort();
913  deleteV.unique();
914  for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
915  if (doDebug) ATH_MSG_VERBOSE("Deleting vertex " << (*dvItr));
916  if (*dvItr) delete (*dvItr);
917  }
918 
919  if (doDebug) ATH_MSG_DEBUG("Deleting particles " << deleteP.size());
920  deleteP.sort();
921  deleteP.unique();
922  for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
923  if (doDebug) ATH_MSG_VERBOSE("Deleting particle " << (*dpItr));
924  if (*dpItr) delete (*dpItr);
925  }
926 
927 #endif
928  // Strip soft underlying stuff
931 
932  if (doDebug && doExtra) {
933  std::cout << "========== BEGIN EVENT BEFORE SOFT ==========" << std::endl;
934  HepMC::Print::line(std::cout, *thinEvt);
935  std::cout << "========== END EVENT BEFORE SOFT ==========" << std::endl;
936  }
937 
938  // Have deleted all hadronization particles
939  // Find all particles connected to hard process(es) with m_T>10GeV
940  std::list<HepMC::GenParticlePtr> pNotHad;
941  std::list<HepMC::GenParticlePtr> pHard;
942 #ifdef HEPMC3
943  std::vector<HepMC::GenParticlePtr> beams=thinEvt->beams();
944  for (auto fp: thinEvt->particles()) {
945  HepMC::GenVertexPtr pvtx = fp->production_vertex();
946  if (!pvtx) continue;
947 
948  double ep = fp->momentum().e();
949  double pzp = fp->momentum().pz();
950  double mtmax = (ep + pzp) * (ep - pzp);
951  auto ancestors = HepMC::ancestor_particles(fp);
952 
953  for (auto gpar: ancestors) {
954  double e = gpar->momentum().e();
955  double pz = gpar->momentum().pz();
956  mtmax = std::max((e+pz)*(e-pz), mtmax);
957  }
958 
959  // Keep hard particles and all ancestors
960  pNotHad.push_back(fp);
961  int ida = std::abs(fp->pdg_id());
962  bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
963  if (mtmax > m_hardCut * m_hardCut || keepid) {
964  pHard.push_back(fp);
965  pHard.insert(pHard.end(),ancestors.begin(),ancestors.end());
966  }
967 
968  // Also keep all descendants of interesting particles
969  // Include leptons to get photons in Sherpa with no Z parent
970  // All hard descendants would include soft initial radiation
971  // Will remove duplicates with list sort/unique
972  bool keepid2 = ( ida == 6) || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
973  if (keepid2 && fp->end_vertex()) {
974  auto descendants = HepMC::descendant_particles(fp);
975  pHard.insert(pHard.end(),descendants.begin(),descendants.end());
976  }
977  }
978 
979  // Sort/unique lists
980  pNotHad.sort();
981  pNotHad.unique();
982  pHard.sort();
983  pHard.unique();
984 
985  // Dump information
986  if (doDebug) {
987  int nhard = 0;
988  for (auto ph: pHard) {
989  ++nhard;
990  ATH_MSG_DEBUG("Hard GenParticles " << ph );
991  }
992  if (doDebug) ATH_MSG_DEBUG("Hard GenParticles total " << nhard);
993  }
994 
995 #else
996 //AV: This algorithm is terribly slow. For each particle the descendants and ancestors are called.
997  HepMC::GenParticle* beams[2];
998  beams[0] = thinEvt->beam_particles().first;
999  beams[1] = thinEvt->beam_particles().second;
1000 
1001  HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
1002  HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();
1003 
1004  for (; finp != finpE; ++finp) {
1005  HepMC::GenParticle* fp = *finp;
1006  HepMC::GenVertex* pvtx = fp->production_vertex();
1007  if (!pvtx) continue;
1008 
1009  double ep = fp->momentum().e();
1010  double pzp = fp->momentum().pz();
1011  double mtmax = (ep + pzp) * (ep - pzp);
1012  HepMC::GenVertex::particle_iterator gpar = fp->production_vertex()->particles_begin(HepMC::ancestors);
1013  HepMC::GenVertex::particle_iterator gparB = gpar;
1014  HepMC::GenVertex::particle_iterator gparE = fp->production_vertex()->particles_end(HepMC::ancestors);
1015 
1016  for (; gpar != gparE; ++gpar) {
1017  double e = (*gpar)->momentum().e();
1018  double pz = (*gpar)->momentum().pz();
1019  mtmax = std::max((e+pz)*(e-pz), mtmax);
1020  }
1021 
1022  // Keep hard particles and all ancestors
1023  pNotHad.push_back(fp);
1024  int ida = std::abs(fp->pdg_id());
1025  bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
1026  if (mtmax > m_hardCut * m_hardCut || keepid) {
1027  pHard.push_back(fp);
1028  for (gpar = gparB; gpar != gparE; ++gpar)
1029  pHard.push_back(*gpar);
1030  }
1031 
1032  // Also keep all descendants of interesting particles
1033  // Include leptons to get photons in Sherpa with no Z parent
1034  // All hard descendants would include soft initial radiation
1035  // Will remove duplicates with list sort/unique
1036  bool keepid2 = ida == 6 || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
1037  if (keepid2 && fp->end_vertex()) {
1038  HepMC::GenVertex::particle_iterator des = fp->end_vertex()->particles_begin(HepMC::descendants);
1039  HepMC::GenVertex::particle_iterator desE = fp->end_vertex()->particles_end(HepMC::descendants);
1040  for (; des != desE; ++des)
1041  pHard.push_back(*des);
1042  }
1043  }
1044 
1045  // Sort/unique lists
1046  pNotHad.sort();
1047  pNotHad.unique();
1048  pHard.sort();
1049  pHard.unique();
1050 
1051  // Dump information
1052  if (doDebug) {
1053  std::list<HepMC::GenParticle*>::iterator hItr2 = pHard.begin();
1054  std::list<HepMC::GenParticle*>::iterator hItr2E = pHard.end();
1055  int nhard = 0;
1056  for (; hItr2 != hItr2E; ++hItr2) {
1057  ++nhard;
1058  ATH_MSG_DEBUG("Hard GenParticles " << (*hItr2) );
1059  }
1060  if (doDebug) ATH_MSG_DEBUG("Hard GenParticles total " << nhard);
1061  }
1062 
1063 #endif
1064  // Remove non-hadronization, non-hard GenParticles from vertices
1065  // and delete them using lists constructed above.
1066  // Any 1->1 vertices created will be removed below.
1067 #ifdef HEPMC3
1068  for (auto p: pNotHad) {
1069  // Skip hard ones
1070  bool isHard = false;
1071  for (auto h: pHard) {
1072  if (p == h) {
1073  isHard = true;
1074  break;
1075  }
1076  }
1077  if (doDebug) ATH_MSG_DEBUG("Particle bc/isHard " << p << " " << isHard);
1078  if (isHard) continue;
1079  HepMC::GenVertexPtr pvtx = p->production_vertex();
1080  if (pvtx) pvtx->remove_particle_out(p);
1081  HepMC::GenVertexPtr evtx = p->end_vertex();
1082  if (evtx) evtx->remove_particle_in(p);
1083  }
1084 #else
1085 
1086  std::list<HepMC::GenParticle*>::iterator pItr = pNotHad.begin();
1087  std::list<HepMC::GenParticle*>::iterator pItrE = pNotHad.end();
1088 
1089  std::list<HepMC::GenParticle*>::iterator hItr = pHard.begin();
1090  std::list<HepMC::GenParticle*>::iterator hItrB = pHard.begin();
1091  std::list<HepMC::GenParticle*>::iterator hItrE = pHard.end();
1092 
1093  for (; pItr != pItrE; ++pItr) {
1094  HepMC::GenParticle* p = *pItr;
1095 
1096  // Skip hard ones
1097  bool isHard = false;
1098  for (hItr = hItrB; hItr != hItrE; ++hItr) {
1099  if (p == (*hItr)) {
1100  isHard = true;
1101  break;
1102  }
1103  }
1104  if (doDebug) ATH_MSG_DEBUG("Particle bc/isHard " << HepMC::uniqueID(p) << " " << isHard);
1105  if (isHard) continue;
1106  HepMC::GenVertex* pvtx = p->production_vertex();
1107  if (pvtx) pvtx->remove_particle(p);
1108  HepMC::GenVertex* evtx = p->end_vertex();
1109  if (evtx) evtx->remove_particle(p);
1110  delete p;
1111  }
1112 #endif
1113 
1115  // Remove and delete vertices with no remaining particles
1117 
1118  removeV.clear();
1119  deleteV.clear();
1120 
1121 #ifdef HEPMC3
1122  for (auto vtx: thinEvt->vertices()) {
1123  if (vtx->particles_in().size() != 0) continue;
1124  if (vtx->particles_out().size() != 0) continue;
1125  removeV.push_back(vtx);
1126  deleteV.push_back(vtx);
1127  }
1128  if (doDebug) ATH_MSG_DEBUG("Removing/deleting 0-particle vertices " << removeV.size() << " " << deleteV.size());
1129  for (unsigned int i = 0; i < removeV.size(); ++i) {
1130  if (doDebug) ATH_MSG_VERBOSE("Removing vertex " << removeV[i]);
1131  thinEvt->remove_vertex(removeV[i]);
1132  }
1133 #else
1134  HepMC::GenEvent::vertex_iterator vtx = thinEvt->vertices_begin();
1135  HepMC::GenEvent::vertex_iterator vtxE = thinEvt->vertices_end();
1136  for (; vtx != vtxE; ++vtx) {
1137  if ((*vtx)->particles_in_size() != 0) continue;
1138  if ((*vtx)->particles_out_size() != 0) continue;
1139  removeV.push_back(*vtx);
1140  deleteV.push_back(*vtx);
1141  }
1142 
1143  if (doDebug) ATH_MSG_DEBUG("Removing/deleting 0-particle vertices " << removeV.size() << " " << deleteV.size());
1144  for (unsigned int i = 0; i < removeV.size(); ++i) {
1145  if (thinEvt->remove_vertex(removeV[i])) {
1146  if (doDebug) ATH_MSG_VERBOSE("Removed vertex " << removeV[i]);
1147  } else {
1148  ATH_MSG_WARNING("Failed to remove vertex " << removeV[i]);
1149  }
1150  }
1151 
1152  deleteV.sort();
1153  deleteV.unique();
1154  for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
1155  if (doDebug) ATH_MSG_VERBOSE("Deleting vertex " << (*dvItr));
1156  if (*dvItr) delete (*dvItr);
1157  }
1158 #endif
1159 
1161  // Remove remaining 1-1 vertices
1163 
1164  if (doDebug && doExtra) {
1165  std::cout << "========== BEGIN EVENT BEFORE 1-BODY ==========" << std::endl;
1166  HepMC::Print::line(std::cout, *thinEvt);
1167  std::cout << "========== END EVENT BEFORE 1-BODY ==========" << std::endl;
1168  }
1169 
1170  // Not clear how to order sweep, so do it one at a time.
1171  // Always keep child
1172  // pvtx -> pin -> vtx11 -> pout -> evtx
1173 
1174 #ifdef HEPMC3
1175  bool moreV1 = true;
1176  HepMC::GenVertexPtr vtx11;
1179 
1180  while (moreV1) {
1181  moreV1 = false;
1182  // Find next 1->1 vertex
1183  for (auto v: thinEvt->vertices()) {
1184  if (v->particles_in().size() != 1) continue;
1185  if (v->particles_out().size() != 1) continue;
1186  pin = *(v->particles_in().begin());
1187  pout = *(v->particles_out().begin());
1188  if (pin->pdg_id() != pout->pdg_id()) continue;
1189  // Sherpa does 1-body decay of incoming protons :-(
1190  if (pin == beams[0] || pin == beams[1]) continue;
1191  HepMC::GenVertexPtr pvtx = pin->production_vertex();
1192  if (!pvtx || pvtx->parent_event() == nullptr) {
1193  ATH_MSG_DEBUG("1->1: missing pvtx for vertex " << v);
1194  ++m_missCount;
1195  continue;
1196  }
1197  moreV1 = true;
1198  vtx11 = v;
1199  if (doDebug) ATH_MSG_DEBUG("One-body " << pin << " " << vtx11 << " " << pout);
1200  break;
1201  }
1202  if (moreV1) {
1203  HepMC::GenVertexPtr pvtx = pin->production_vertex();
1204  pvtx->remove_particle_in(pin);
1205  pvtx->remove_particle_out(pin);
1206  pvtx->add_particle_out(pout);
1207  vtx11->remove_particle_in(pin);
1208  vtx11->remove_particle_out(pin);
1209  vtx11->remove_particle_in(pout);
1210  vtx11->remove_particle_out(pout);
1211  thinEvt->remove_vertex(vtx11);
1212  if (doDebug) ATH_MSG_DEBUG("One-body new pvtx " << pvtx << " " << pvtx->particles_in().size() << " " << pvtx->particles_out().size());
1213  }
1214  }
1215 #else
1216  bool moreV1 = true;
1217  HepMC::GenVertex* vtx11;
1218  HepMC::GenParticle* pin;
1220 
1221  while (moreV1) {
1222  moreV1 = false;
1223 
1224  HepMC::GenEvent::vertex_iterator v = thinEvt->vertices_begin();
1225  HepMC::GenEvent::vertex_iterator vE = thinEvt->vertices_end();
1226 
1227  // Find next 1->1 vertex
1228  for (; v != vE; ++v) {
1229  if ((*v)->particles_in_size() != 1) continue;
1230  if ((*v)->particles_out_size() != 1) continue;
1231  pin = *((*v)->particles_in_const_begin());
1232  pout = *((*v)->particles_out_const_begin());
1233  if (pin->pdg_id() != pout->pdg_id()) continue;
1234  // Sherpa does 1-body decay of incoming protons :-(
1235  if (pin == beams[0] || pin == beams[1]) continue;
1236  HepMC::GenVertex* pvtx = pin->production_vertex();
1237  if (!pvtx || HepMC::uniqueID(pvtx) == 0) {
1238  ATH_MSG_DEBUG("1->1: missing pvtx for vertex " << HepMC::uniqueID(*v));
1239  ++m_missCount;
1240  continue;
1241  }
1242 
1243  moreV1 = true;
1244  vtx11 = (*v);
1245  if (doDebug) ATH_MSG_DEBUG("One-body " << HepMC::uniqueID(pin) << " " << HepMC::uniqueID(vtx11) << " " << HepMC::uniqueID(pout));
1246  break;
1247  }
1248  if (moreV1) {
1249  HepMC::GenVertex* pvtx = pin->production_vertex();
1250  pvtx->remove_particle(pin);
1251  pvtx->add_particle_out(pout);
1252  vtx11->remove_particle(pin);
1253  vtx11->remove_particle(pout);
1254  thinEvt->remove_vertex(vtx11);
1255  delete pin;
1256  delete vtx11;
1257  if (doDebug) ATH_MSG_DEBUG("One-body new pvtx " << HepMC::uniqueID(pvtx) << " " << pvtx->particles_in_size() << " " << pvtx->particles_out_size());
1258  }
1259  }
1260 
1261 #endif
1262  // Remove dangling particles/vertices
1265 
1266  // GEN_EVENT -> GEN_AOD conversion applies pt cut to partons, breaking
1267  // tree structure. Result is "dangling" partons with 1 -> 0 vertices.
1268  // FIXME!! Meanwhile discard these if pt < m_danglePtCut.
1269 
1270  if (m_danglePtCut > 0) {
1271 
1272  removePV.clear();
1273  removeV.clear();
1274  deleteP.clear();
1275  deleteV.clear();
1276 
1277 #ifdef HEPMC3
1278  for (auto badv: thinEvt->vertices()) {
1279  if (!badv) continue;
1280  if (badv->particles_in().size() != 1 || badv->particles_out().size() != 0) continue;
1281  auto pitr = badv->particles_in().begin();
1282  HepMC::GenParticlePtr pp = *pitr;
1283  if (pp->production_vertex()) continue;
1284  double pt = pp->momentum().perp();
1285  if (pt > m_danglePtMax) m_danglePtMax = pt;
1286  ++m_dangleFound;
1287  if (pt > m_danglePtCut) continue;
1288  if (doDebug) ATH_MSG_DEBUG("1->0: removing pp,badv,pt " << pp << " " << badv << " " << pt);
1289  removePV.push_back(vpPair(badv, pp));
1290  deleteP.push_back(pp);
1291  removeV.push_back(badv);
1292  deleteV.push_back(badv);
1293  ++m_dangleRemoved;
1294  }
1295  // Actually implement changes -- remove particles from vertices
1296  for (unsigned int i = 0; i < removePV.size(); ++i) {
1297  HepMC::GenVertexPtr v = removePV[i].first;
1298  HepMC::GenParticlePtr p = removePV[i].second;
1299  v->remove_particle_in(p);
1300  v->remove_particle_out(p);
1301  }
1302  // Actually implement changes -- remove vertices
1303  for (unsigned int i = 0; i < removeV.size(); ++i) {
1304  thinEvt->remove_vertex(removeV[i]);
1305  }
1306 //AV: HepMC3 uses smart pointers
1307 #else
1308  HepMC::GenEvent::vertex_iterator badv = thinEvt->vertices_begin();
1309  HepMC::GenEvent::vertex_iterator badvE = thinEvt->vertices_end();
1310 
1311  for (; badv != badvE; ++badv) {
1312  if (!(*badv)) continue;
1313  if ((*badv)->particles_in_size() != 1 || (*badv)->particles_out_size() != 0) continue;
1314  HepMC::GenVertex::particles_in_const_iterator pitr = (*badv)->particles_in_const_begin();
1315  HepMC::GenParticle* pp = *pitr;
1316  if (pp->production_vertex()) continue;
1317  double pt = pp->momentum().perp();
1318  if (pt > m_danglePtMax) m_danglePtMax = pt;
1319  ++m_dangleFound;
1320  if (pt > m_danglePtCut) continue;
1321  if (doDebug) ATH_MSG_DEBUG("1->0: removing pp,badv,pt " << pp << " " << *badv << " " << pt);
1322  removePV.emplace_back(*badv, pp);
1323  deleteP.push_back(pp);
1324  removeV.push_back(*badv);
1325  deleteV.push_back(*badv);
1326  ++m_dangleRemoved;
1327  }
1328 
1329  // Actually implement changes -- remove particles from vertices
1330  for (unsigned int i = 0; i < removePV.size(); ++i) {
1331  HepMC::GenVertex* v = removePV[i].first;
1332  HepMC::GenParticle* p = removePV[i].second;
1333  v->remove_particle(p);
1334  }
1335 
1336  // Actually implement changes -- remove vertices
1337  for (unsigned int i = 0; i < removeV.size(); ++i) {
1338  if (!thinEvt->remove_vertex(removeV[i])) { ATH_MSG_WARNING("1->0: Failed to remove vertex " << removeV[i]); }
1339  }
1340 
1341  // Delete removed particles/vertices
1342  deleteV.sort();
1343  deleteV.unique();
1344  for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
1345  if (*dvItr) delete (*dvItr);
1346  }
1347 
1348  deleteP.sort();
1349  deleteP.unique();
1350  for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
1351  if (*dpItr) delete (*dpItr);
1352  }
1353 #endif
1354  } // end m_danglePtCut
1355 
1357  // Done - examine results
1359 
1360  if (doPrint) {
1361  std::cout << "========== BEGIN EVENT AFTER THINNING ==========" << std::endl;
1362  HepMC::Print::line(std::cout, *thinEvt);
1363  std::cout << "========== END EVENT AFTER THINNING ==========" << std::endl;
1364  }
1365 
1366  m_thinParticles += thinEvt->particles_size();
1367  m_thinVertices += thinEvt->vertices_size();
1368 
1370  // Save thinned event in output container
1372 
1373  thinnedMcEvts->push_back(thinEvt);
1374  return StatusCode::SUCCESS;
1375 }

◆ 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 }

◆ finalize()

StatusCode DerivationFramework::CompactHardTruth::finalize ( )
virtual

Definition at line 91 of file CompactHardTruth.cxx.

91  {
92 
93  ATH_MSG_INFO("Finalizing DerivationFramework::CompactHardTruth ");
94  ATH_MSG_INFO("Missing items limiting reclustering " << m_missCount);
95 
96  ATH_MSG_INFO("Dangling partons pt cut: " << m_danglePtCut);
97  ATH_MSG_INFO("Dangling partons found: " << m_dangleFound);
98  ATH_MSG_INFO("Dangling partons removed: " << m_dangleRemoved);
99  ATH_MSG_INFO("Dangling partons max pt: " << m_danglePtMax);
100 
101  ATH_MSG_INFO("CompactHardTruth total particles: " << m_thinParticles);
102  ATH_MSG_INFO("CompactHardTruth total vertices: " << m_thinVertices);
103 
104  return StatusCode::SUCCESS;
105 }

◆ initialize()

StatusCode DerivationFramework::CompactHardTruth::initialize ( )
virtual

Definition at line 67 of file CompactHardTruth.cxx.

67  {
68  ATH_MSG_INFO("Initializing " << name() << "...");
69 
70  m_evtCount = -1;
71  m_missCount = 0;
72  m_dangleFound = 0;
73  m_dangleRemoved = 0;
74  m_danglePtMax = 0;
75 
76  m_thinParticles = 0;
77  m_thinVertices = 0;
78 
79  // Print jobOption inputs
80  ATH_MSG_INFO("-------------------------------------------------");
81  ATH_MSG_INFO("jobOption McEvent " << m_mcEventsName);
82  ATH_MSG_INFO("jobOption McEventOut " << m_thinnedMcEventsName);
83  ATH_MSG_INFO("jobOption ShowerMassCut " << m_partonCut);
84  ATH_MSG_INFO("jobOption SoftMtCut " << m_hardCut);
85  ATH_MSG_INFO("jobOption MaxCount " << m_maxCount);
86  ATH_MSG_INFO("-------------------------------------------------");
87 
88  return StatusCode::SUCCESS;
89 }

◆ 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.

◆ msg() [1/2]

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

Definition at line 24 of file AthCommonMsg.h.

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

◆ msg() [2/2]

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

Definition at line 27 of file AthCommonMsg.h.

27  {
28  return this->msgStream(lvl);
29  }

◆ 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.

◆ 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();
383  PBASE::renounce (h);
384  }

◆ 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  {
365  handlesArray.renounce();
366  }

◆ 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, PyAthena::Alg, and AthHistogramAlgorithm.

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 }

◆ 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.

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

◆ vtxInMom()

HepMC::FourVector DerivationFramework::CompactHardTruth::vtxInMom ( HepMC::ConstGenVertexPtr  v)
virtual

Definition at line 1383 of file CompactHardTruth.cxx.

1383  {
1384 #ifdef HEPMC3
1385  HepMC::FourVector ret(0.0,0.0,0.0,0.0);
1386  for (auto p: v->particles_in()) ret+=p->momentum();
1387  return ret;
1388 #else
1389  double px = 0;
1390  double py = 0;
1391  double pz = 0;
1392  double e = 0;
1393  HepMC::GenVertex::particles_in_const_iterator it = v->particles_in_const_begin();
1394  HepMC::GenVertex::particles_in_const_iterator itE = v->particles_in_const_end();
1395  for (; it != itE; ++it) {
1396  px += (*it)->momentum().px();
1397  py += (*it)->momentum().py();
1398  pz += (*it)->momentum().pz();
1399  e += (*it)->momentum().e();
1400  }
1401  return HepMC::FourVector(px, py, pz, e);
1402 #endif
1403 }

◆ vtxOutMom()

HepMC::FourVector DerivationFramework::CompactHardTruth::vtxOutMom ( HepMC::ConstGenVertexPtr  v)
virtual

Definition at line 1405 of file CompactHardTruth.cxx.

1405  {
1406 #ifdef HEPMC3
1407  HepMC::FourVector ret(0.0,0.0,0.0,0.0);
1408  for (auto p: v->particles_out()) ret+=p->momentum();
1409  return ret;
1410 #else
1411  double px = 0;
1412  double py = 0;
1413  double pz = 0;
1414  double e = 0;
1415  HepMC::GenVertex::particles_out_const_iterator it = v->particles_out_const_begin();
1416  HepMC::GenVertex::particles_out_const_iterator itE = v->particles_out_const_end();
1417  for (; it != itE; ++it) {
1418  px += (*it)->momentum().px();
1419  py += (*it)->momentum().py();
1420  pz += (*it)->momentum().pz();
1421  e += (*it)->momentum().e();
1422  }
1423  return HepMC::FourVector(px, py, pz, e);
1424 #endif
1425 }

Member Data Documentation

◆ m_dangleFound

int DerivationFramework::CompactHardTruth::m_dangleFound = 0
private

Definition at line 87 of file CompactHardTruth.h.

◆ m_danglePtCut

float DerivationFramework::CompactHardTruth::m_danglePtCut
private

Definition at line 90 of file CompactHardTruth.h.

◆ m_danglePtMax

float DerivationFramework::CompactHardTruth::m_danglePtMax = 0.0F
private

Definition at line 89 of file CompactHardTruth.h.

◆ m_dangleRemoved

int DerivationFramework::CompactHardTruth::m_dangleRemoved = 0
private

Definition at line 88 of file CompactHardTruth.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_evtCount

int DerivationFramework::CompactHardTruth::m_evtCount = 0
private

Definition at line 82 of file CompactHardTruth.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_hardCut

float DerivationFramework::CompactHardTruth::m_hardCut
private

Definition at line 85 of file CompactHardTruth.h.

◆ m_maxCount

int DerivationFramework::CompactHardTruth::m_maxCount
private

Definition at line 92 of file CompactHardTruth.h.

◆ m_mcEventsName

std::string DerivationFramework::CompactHardTruth::m_mcEventsName
private

Containers.

Definition at line 78 of file CompactHardTruth.h.

◆ m_missCount

int DerivationFramework::CompactHardTruth::m_missCount = 0
private

Definition at line 83 of file CompactHardTruth.h.

◆ m_partonCut

float DerivationFramework::CompactHardTruth::m_partonCut
private

Definition at line 84 of file CompactHardTruth.h.

◆ m_thinnedMcEventsName

std::string DerivationFramework::CompactHardTruth::m_thinnedMcEventsName
private

Definition at line 79 of file CompactHardTruth.h.

◆ m_thinParticles

int DerivationFramework::CompactHardTruth::m_thinParticles = 0
private

Definition at line 94 of file CompactHardTruth.h.

◆ m_thinVertices

int DerivationFramework::CompactHardTruth::m_thinVertices = 0
private

Definition at line 95 of file CompactHardTruth.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:
HepMC::GenVertexPtr
HepMC::GenVertex * GenVertexPtr
Definition: GenVertex.h:59
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
DerivationFramework::CompactHardTruth::m_maxCount
int m_maxCount
Definition: CompactHardTruth.h:92
DerivationFramework::CompactHardTruth::m_missCount
int m_missCount
Definition: CompactHardTruth.h:83
calibdata.pout
def pout(output, newline=True)
Definition: calibdata.py:129
test_pyathena.px
px
Definition: test_pyathena.py:18
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
HepMC::GenParticlePtr
GenParticle * GenParticlePtr
Definition: GenParticle.h:37
McEventCollection
McEventCollection
Definition: GeneratorObjectsTPCnv.cxx:60
skel.it
it
Definition: skel.GENtoEVGEN.py:407
python.SystemOfUnits.second
float second
Definition: SystemOfUnits.py:135
test_pyathena.pt
pt
Definition: test_pyathena.py:11
AthCommonDataStore< AthCommonMsg< Algorithm > >::m_evtStore
StoreGateSvc_t m_evtStore
Pointer to StoreGate (event store by default)
Definition: AthCommonDataStore.h:390
AthCommonDataStore< AthCommonMsg< Algorithm > >::m_vhka
std::vector< SG::VarHandleKeyArray * > m_vhka
Definition: AthCommonDataStore.h:398
isParton
bool isParton(const T &p)
Definition: AtlasPID.h:1077
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
HepMC::Print::line
void line(std::ostream &os, const GenEvent &e)
Definition: GenEvent.h:676
TruthTest.itE
itE
Definition: TruthTest.py:25
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
DerivationFramework::CompactHardTruth::m_dangleFound
int m_dangleFound
Definition: CompactHardTruth.h:87
SG::VarHandleKeyArray::setOwner
virtual void setOwner(IDataHandleHolder *o)=0
IDTPMcnv.htype
htype
Definition: IDTPMcnv.py:29
DerivationFramework::CompactHardTruth::m_danglePtCut
float m_danglePtCut
Definition: CompactHardTruth.h:90
AthCommonDataStore::declareGaudiProperty
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>
Definition: AthCommonDataStore.h:156
AthCommonDataStore< AthCommonMsg< Algorithm > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
AthCommonDataStore
Definition: AthCommonDataStore.h:52
AthAlgorithm::sysInitialize
virtual StatusCode sysInitialize() override
Override sysInitialize.
Definition: AthAlgorithm.cxx:66
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
AthCommonDataStore< AthCommonMsg< Algorithm > >::outputHandles
virtual std::vector< Gaudi::DataHandle * > outputHandles() const override
Return this algorithm's output handles.
HepMC::is_simulation_particle
bool is_simulation_particle(const T &p)
Method to establish if a particle (or barcode) was created during the simulation (TODO update to be s...
Definition: MagicNumbers.h:354
lumiFormat.i
int i
Definition: lumiFormat.py:85
trigmenu_modify_prescale_json.fp
fp
Definition: trigmenu_modify_prescale_json.py:53
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition: AthCommonDataStore.h:145
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:116
DataVector::front
const T * front() const
Access the first element in the collection as an rvalue.
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
HepMC::is_simulation_vertex
bool is_simulation_vertex(const T &v)
Method to establish if the vertex was created during simulation (TODO migrate to be based on status).
Definition: MagicNumbers.h:360
AthCommonDataStore< AthCommonMsg< Algorithm > >::m_detStore
StoreGateSvc_t m_detStore
Pointer to StoreGate (detector store by default)
Definition: AthCommonDataStore.h:393
McEventCollection
This defines the McEventCollection, which is really just an ObjectVector of McEvent objects.
Definition: McEventCollection.h:32
SG::VarHandleKeyArray::renounce
virtual void renounce()=0
SG::HandleClassifier::type
std::conditional< std::is_base_of< SG::VarHandleKeyArray, T >::value, VarHandleKeyArrayType, type2 >::type type
Definition: HandleClassifier.h:54
Amg::py
@ py
Definition: GeoPrimitives.h:39
DerivationFramework::CompactHardTruth::m_thinVertices
int m_thinVertices
Definition: CompactHardTruth.h:95
DerivationFramework::CompactHardTruth::m_mcEventsName
std::string m_mcEventsName
Containers.
Definition: CompactHardTruth.h:78
merge_scale_histograms.doc
string doc
Definition: merge_scale_histograms.py:9
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
isHadron
bool isHadron(const T &p)
Definition: AtlasPID.h:348
DerivationFramework::CompactHardTruth::m_partonCut
float m_partonCut
Definition: CompactHardTruth.h:84
DerivationFramework::CompactHardTruth::m_hardCut
float m_hardCut
Definition: CompactHardTruth.h:85
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
python.PyAthena.v
v
Definition: PyAthena.py:154
AthAlgorithm::m_extendedExtraObjects
DataObjIDColl m_extendedExtraObjects
Definition: AthAlgorithm.h:79
a
TList * a
Definition: liststreamerinfos.cxx:10
h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
DerivationFramework::CompactHardTruth::m_danglePtMax
float m_danglePtMax
Definition: CompactHardTruth.h:89
DerivationFramework::CompactHardTruth::m_evtCount
int m_evtCount
Definition: CompactHardTruth.h:82
SG::VarHandleBase::vhKey
SG::VarHandleKey & vhKey()
Return a non-const reference to the HandleKey.
Definition: StoreGate/src/VarHandleBase.cxx:629
DerivationFramework::CompactHardTruth::m_thinnedMcEventsName
std::string m_thinnedMcEventsName
Definition: CompactHardTruth.h:79
AthAlgorithm::AthAlgorithm
AthAlgorithm()
Default constructor:
python.Bindings.keys
keys
Definition: Control/AthenaPython/python/Bindings.py:801
LHEF::Writer
Pythia8::Writer Writer
Definition: Prophecy4fMerger.cxx:12
isDiquark
bool isDiquark(const T &p)
PDG rule 4 Diquarks have 4-digit numbers with nq1 >= nq2 and nq3 = 0 APID: states with top quarks are...
Definition: AtlasPID.h:224
DerivationFramework::CompactHardTruth::m_dangleRemoved
int m_dangleRemoved
Definition: CompactHardTruth.h:88
DerivationFramework::CompactHardTruth::vtxOutMom
virtual HepMC::FourVector vtxOutMom(HepMC::ConstGenVertexPtr)
Definition: CompactHardTruth.cxx:1405
DataVector::empty
bool empty() const noexcept
Returns true if the collection is empty.
fitman.k
k
Definition: fitman.py:528
GenParticle
@ GenParticle
Definition: TruthClasses.h:30
ServiceHandle< ICondSvc >
DerivationFramework::CompactHardTruth::m_thinParticles
int m_thinParticles
Definition: CompactHardTruth.h:94