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 > &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 > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleKeyArrayType &)
 specialization for handling Gaudi::Property<SG::VarHandleKeyArray> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleType &)
 specialization for handling Gaudi::Property<SG::VarHandleBase> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T > &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 47 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 44 of file CompactHardTruth.cxx.

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

◆ ~CompactHardTruth()

DerivationFramework::CompactHardTruth::~CompactHardTruth ( )
virtual

Destructor:

Definition at line 66 of file CompactHardTruth.cxx.

66 {}

◆ 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 > &  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 > &  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 > &  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 > &  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 > &  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 110 of file CompactHardTruth.cxx.

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

◆ 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 94 of file CompactHardTruth.cxx.

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

◆ initialize()

StatusCode DerivationFramework::CompactHardTruth::initialize ( )
virtual

Definition at line 70 of file CompactHardTruth.cxx.

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

◆ 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 1386 of file CompactHardTruth.cxx.

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

◆ vtxOutMom()

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

Definition at line 1408 of file CompactHardTruth.cxx.

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

Member Data Documentation

◆ m_dangleFound

int DerivationFramework::CompactHardTruth::m_dangleFound = 0
private

Definition at line 88 of file CompactHardTruth.h.

◆ m_danglePtCut

float DerivationFramework::CompactHardTruth::m_danglePtCut
private

Definition at line 91 of file CompactHardTruth.h.

◆ m_danglePtMax

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

Definition at line 90 of file CompactHardTruth.h.

◆ m_dangleRemoved

int DerivationFramework::CompactHardTruth::m_dangleRemoved = 0
private

Definition at line 89 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 83 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 86 of file CompactHardTruth.h.

◆ m_maxCount

int DerivationFramework::CompactHardTruth::m_maxCount
private

Definition at line 93 of file CompactHardTruth.h.

◆ m_mcEventsName

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

Containers.

Definition at line 79 of file CompactHardTruth.h.

◆ m_missCount

int DerivationFramework::CompactHardTruth::m_missCount = 0
private

Definition at line 84 of file CompactHardTruth.h.

◆ m_partonCut

float DerivationFramework::CompactHardTruth::m_partonCut
private

Definition at line 85 of file CompactHardTruth.h.

◆ m_thinnedMcEventsName

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

Definition at line 80 of file CompactHardTruth.h.

◆ m_thinParticles

int DerivationFramework::CompactHardTruth::m_thinParticles = 0
private

Definition at line 95 of file CompactHardTruth.h.

◆ m_thinVertices

int DerivationFramework::CompactHardTruth::m_thinVertices = 0
private

Definition at line 96 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
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
DerivationFramework::CompactHardTruth::m_maxCount
int m_maxCount
Definition: CompactHardTruth.h:93
DerivationFramework::CompactHardTruth::m_missCount
int m_missCount
Definition: CompactHardTruth.h:84
calibdata.pout
def pout(output, newline=True)
Definition: calibdata.py:130
test_pyathena.px
px
Definition: test_pyathena.py:18
max
#define max(a, b)
Definition: cfImp.cxx:41
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
HepMC::GenParticlePtr
GenParticle * GenParticlePtr
Definition: GenParticle.h:37
McEventCollection
McEventCollection
Definition: GeneratorObjectsTPCnv.cxx:60
skel.it
it
Definition: skel.GENtoEVGEN.py:396
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:611
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:88
SG::VarHandleKeyArray::setOwner
virtual void setOwner(IDataHandleHolder *o)=0
IDTPMcnv.htype
htype
Definition: IDTPMcnv.py:27
DerivationFramework::CompactHardTruth::m_danglePtCut
float m_danglePtCut
Definition: CompactHardTruth.h:91
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:210
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:342
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
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:109
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:348
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:33
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:96
DerivationFramework::CompactHardTruth::m_mcEventsName
std::string m_mcEventsName
Containers.
Definition: CompactHardTruth.h:79
merge_scale_histograms.doc
string doc
Definition: merge_scale_histograms.py:9
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
isHadron
bool isHadron(const T &p)
Definition: AtlasPID.h:218
DerivationFramework::CompactHardTruth::m_partonCut
float m_partonCut
Definition: CompactHardTruth.h:85
DerivationFramework::CompactHardTruth::m_hardCut
float m_hardCut
Definition: CompactHardTruth.h:86
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:90
DerivationFramework::CompactHardTruth::m_evtCount
int m_evtCount
Definition: CompactHardTruth.h:83
SG::VarHandleBase::vhKey
SG::VarHandleKey & vhKey()
Return a non-const reference to the HandleKey.
Definition: StoreGate/src/VarHandleBase.cxx:623
DerivationFramework::CompactHardTruth::m_thinnedMcEventsName
std::string m_thinnedMcEventsName
Definition: CompactHardTruth.h:80
AthAlgorithm::AthAlgorithm
AthAlgorithm()
Default constructor:
python.Bindings.keys
keys
Definition: Control/AthenaPython/python/Bindings.py:798
LHEF::Writer
Pythia8::Writer Writer
Definition: Prophecy4fMerger.cxx:12
isDiquark
bool isDiquark(const T &p)
Definition: AtlasPID.h:217
DerivationFramework::CompactHardTruth::m_dangleRemoved
int m_dangleRemoved
Definition: CompactHardTruth.h:89
AthCommonDataStore::declareGaudiProperty
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>
Definition: AthCommonDataStore.h:156
DerivationFramework::CompactHardTruth::vtxOutMom
virtual HepMC::FourVector vtxOutMom(HepMC::ConstGenVertexPtr)
Definition: CompactHardTruth.cxx:1408
DataVector::empty
bool empty() const noexcept
Returns true if the collection is empty.
query_example.des
des
Definition: query_example.py:9
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:95