ATLAS Offline Software
CompactHardTruth.cxx
Go to the documentation of this file.
1 
3 /*
4  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
5 */
6 
7 // CompactHardTruth.cxx
8 // Implementation file for class CompactHardTruth
9 // Author: Frank Paige <paige@bnl.gov>
11 
12 // McParticleTests includes
14 
15 // STL includes
16 #include <list>
17 
18 // FrameWork includes
19 #include "Gaudi/Property.h"
20 #include "GaudiKernel/ITHistSvc.h"
21 #include "GaudiKernel/ServiceHandle.h"
22 
23 #include "AtlasHepMC/GenEvent.h"
24 #include "AtlasHepMC/GenParticle.h"
25 #include "AtlasHepMC/GenVertex.h"
26 #include "AtlasHepMC/Relatives.h"
29 // Needed for FourVector
32 
33 // ROOT includes
34 #include "TH1F.h"
35 
36 namespace DerivationFramework {
37 
39 // Public methods:
41 
42 // Constructors
44 CompactHardTruth::CompactHardTruth(const std::string& name, ISvcLocator* pSvcLocator)
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 }
63 
64 // Destructor
67 
68 // Athena Algorithm's Hooks
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 }
93 
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 }
109 
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 }
1379 
1381 // Const methods:
1383 
1384 // Total cluster FourVectors
1385 
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 }
1407 
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 }
1429 
1430 } // namespace DerivationFramework
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
DerivationFramework::CompactHardTruth::execute
virtual StatusCode execute()
Definition: CompactHardTruth.cxx:110
calibdata.pout
def pout(output, newline=True)
Definition: calibdata.py:130
GenEvent.h
test_pyathena.px
px
Definition: test_pyathena.py:18
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
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
GenVertex.h
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
isParton
bool isParton(const T &p)
Definition: AtlasPID.h:841
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
GenParticle.h
TruthTest.itE
itE
Definition: TruthTest.py:25
DerivationFramework::CompactHardTruth::m_dangleFound
int m_dangleFound
Definition: CompactHardTruth.h:88
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
SimpleVector.h
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:355
McEventCollection.h
lumiFormat.i
int i
Definition: lumiFormat.py:85
trigmenu_modify_prescale_json.fp
fp
Definition: trigmenu_modify_prescale_json.py:53
DerivationFramework::CompactHardTruth::initialize
virtual StatusCode initialize()
Definition: CompactHardTruth.cxx:70
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: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:361
DerivationFramework
THE reconstruction tool.
Definition: ParticleSortingAlg.h:24
DerivationFramework::CompactHardTruth::~CompactHardTruth
virtual ~CompactHardTruth()
Destructor:
Definition: CompactHardTruth.cxx:66
McEventCollection
This defines the McEventCollection, which is really just an ObjectVector of McEvent objects.
Definition: McEventCollection.h:33
AthAlgorithm
Definition: AthAlgorithm.h:47
CompactHardTruth.h
DerivationFramework::CompactHardTruth::finalize
virtual StatusCode finalize()
Definition: CompactHardTruth.cxx:94
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
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
MagicNumbers.h
isHadron
bool isHadron(const T &p)
Definition: AtlasPID.h:300
DerivationFramework::CompactHardTruth::m_partonCut
float m_partonCut
Definition: CompactHardTruth.h:85
DerivationFramework::CompactHardTruth::m_hardCut
float m_hardCut
Definition: CompactHardTruth.h:86
DerivationFramework::CompactHardTruth::vtxInMom
virtual HepMC::FourVector vtxInMom(HepMC::ConstGenVertexPtr)
Definition: CompactHardTruth.cxx:1386
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
h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Relatives.h
DerivationFramework::CompactHardTruth::m_danglePtMax
float m_danglePtMax
Definition: CompactHardTruth.h:90
DerivationFramework::CompactHardTruth::m_evtCount
int m_evtCount
Definition: CompactHardTruth.h:83
DerivationFramework::CompactHardTruth::m_thinnedMcEventsName
std::string m_thinnedMcEventsName
Definition: CompactHardTruth.h:80
DerivationFramework::CompactHardTruth::CompactHardTruth
CompactHardTruth()
Default constructor:
HepMC::ConstGenVertexPtr
const HepMC::GenVertex * ConstGenVertexPtr
Definition: GenVertex.h:60
isDiquark
bool isDiquark(const T &p)
PDG rule 4 Diquarks have 4-digit numbers with nq1 >= nq2 and nq3 = 0 APID: the diquarks with fourth g...
Definition: AtlasPID.h:186
DerivationFramework::CompactHardTruth::m_dangleRemoved
int m_dangleRemoved
Definition: CompactHardTruth.h:89
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.
HepMCHelpers.h
query_example.des
des
Definition: query_example.py:9
GenParticle
@ GenParticle
Definition: TruthClasses.h:30
DerivationFramework::CompactHardTruth::m_thinParticles
int m_thinParticles
Definition: CompactHardTruth.h:95