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