CompactHardTruth.cxx 53.9 KB
Newer Older
1
2
3
///////////////////////// -*- C++ -*- /////////////////////////////

/*
4
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
5
6
*/

7
// CompactHardTruth.cxx
8
9
// Implementation file for class CompactHardTruth
// Author: Frank Paige <paige@bnl.gov>
10
///////////////////////////////////////////////////////////////////
11
12
13
14
15
16
17
18

// McParticleTests includes
#include "DerivationFrameworkMCTruth/CompactHardTruth.h"

// STL includes
#include <list>

// FrameWork includes
Charles Leggett's avatar
Charles Leggett committed
19
#include "Gaudi/Property.h"
20
21
22
#include "GaudiKernel/ITHistSvc.h"
#include "GaudiKernel/ServiceHandle.h"

23
24
25
#include "AtlasHepMC/GenEvent.h"
#include "AtlasHepMC/GenParticle.h"
#include "AtlasHepMC/GenVertex.h"
26
#include "AtlasHepMC/Relatives.h"
27
#include "GeneratorObjects/McEventCollection.h"
28
// Needed for FourVector
29
#include "AtlasHepMC/SimpleVector.h"
30
31
32
33
34
35

// ROOT includes
#include "TH1F.h"

namespace DerivationFramework {

36
37
38
///////////////////////////////////////////////////////////////////
// Public methods:
///////////////////////////////////////////////////////////////////
39
40
41

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

  declareProperty("MaxCount", m_maxCount, "maximum number of events to print");
60
61
62
63
}

// Destructor
///////////////
64
CompactHardTruth::~CompactHardTruth() {}
65
66
67

// Athena Algorithm's Hooks
////////////////////////////
68
69
StatusCode CompactHardTruth::initialize() {
  ATH_MSG_INFO("Initializing " << name() << "...");
70
71
72
73
74
75
76
77
78
79
80

  m_evtCount = -1;
  m_missCount = 0;
  m_dangleFound = 0;
  m_dangleRemoved = 0;
  m_danglePtMax = 0;

  m_thinParticles = 0;
  m_thinVertices = 0;

  // Print jobOption inputs
81
82
83
84
85
86
87
  ATH_MSG_INFO("-------------------------------------------------");
  ATH_MSG_INFO("jobOption McEvent            " << m_mcEventsName);
  ATH_MSG_INFO("jobOption McEventOut         " << m_thinnedMcEventsName);
  ATH_MSG_INFO("jobOption ShowerMassCut      " << m_partonCut);
  ATH_MSG_INFO("jobOption SoftMtCut          " << m_hardCut);
  ATH_MSG_INFO("jobOption MaxCount           " << m_maxCount);
  ATH_MSG_INFO("-------------------------------------------------");
88
89
90
91

  return StatusCode::SUCCESS;
}

92
StatusCode CompactHardTruth::finalize() {
93
94

  ATH_MSG_INFO("Finalizing DerivationFramework::CompactHardTruth ");
95
  ATH_MSG_INFO("Missing items limiting reclustering " << m_missCount);
96

97
98
99
100
  ATH_MSG_INFO("Dangling partons pt cut:  " << m_danglePtCut);
  ATH_MSG_INFO("Dangling partons found:   " << m_dangleFound);
  ATH_MSG_INFO("Dangling partons removed: " << m_dangleRemoved);
  ATH_MSG_INFO("Dangling partons max pt:  " << m_danglePtMax);
101

102
103
  ATH_MSG_INFO("CompactHardTruth total particles:  " << m_thinParticles);
  ATH_MSG_INFO("CompactHardTruth total vertices:   " << m_thinVertices);
104
105
106
107

  return StatusCode::SUCCESS;
}

108
StatusCode CompactHardTruth::execute() {
109
110

  ++m_evtCount;
111
112
  // if( m_evtCount%100 == 0 ){
  ATH_MSG_INFO("Executing " << name() << " " << m_evtCount);
113
114
115
  //}

  // Normally doPrint is used to print the first m_maxCount events
116
  // before and after thinning.
117
118
119
120
121
  // doExtra adds extra intermediate event printouts.
  // doDebug allows debug printout for a range of events.
  bool doPrint = m_evtCount < m_maxCount;
  bool doDebug = false;
  bool doExtra = false;
122
123
  // doDebug = doPrint;
  // doExtra = doPrint;
124
125

  // Retrieve input data
126
127
  const McEventCollection* mcEvts = nullptr;
  if (!evtStore()->retrieve(mcEvts, m_mcEventsName).isSuccess() || nullptr == mcEvts) {
128
    ATH_MSG_WARNING("could not retrieve mc collection at [" << m_mcEventsName << "]!");
129
130
131
132
    return StatusCode::FAILURE;
  }

  if (mcEvts->empty()) {
133
    ATH_MSG_WARNING("empty McEventCollection at [" << m_mcEventsName << "]");
134
135
136
137
138
    return StatusCode::SUCCESS;
  }

  // Create output collection
  McEventCollection* thinnedMcEvts = new McEventCollection;
139
140
141
  if (!evtStore()->record(thinnedMcEvts, m_thinnedMcEventsName).isSuccess()) {
    ATH_MSG_WARNING("Could not record thinned mc collection at [" << m_thinnedMcEventsName << "]!");
    delete thinnedMcEvts;
142
    thinnedMcEvts = nullptr;
143
144
    return StatusCode::FAILURE;
  }
145
  if (evtStore()->setConst(thinnedMcEvts).isFailure()) { ATH_MSG_WARNING("Could not lock the McEventCollection at [" << m_thinnedMcEventsName << "] !!"); }
146
147
148

  // Signal event is first (only?) event; front() is from DataVector
  const HepMC::GenEvent* mcEvt = mcEvts->front();
149
  auto wtCont = mcEvt->weights();
150
  if (!wtCont.empty()) {
151
  } else {
152
    ATH_MSG_WARNING("Weights not found for mc collection [" << m_mcEventsName << "]");
153
154
  }
  int inEvent = mcEvt->event_number();
155
  if (doDebug) ATH_MSG_DEBUG("FETCHED count/event " << m_evtCount << " " << inEvent);
156
157
158
159
160
161
162

  ///////////////////////////////
  // New event - copy of original
  ///////////////////////////////

  HepMC::GenEvent* thinEvt = new HepMC::GenEvent(*mcEvt);
  int nEvent = thinEvt->event_number();
163
  if (doPrint) ATH_MSG_DEBUG("New event number = " << nEvent);
164

165
166
  if (doPrint) {
    std::cout << "========== BEGIN EVENT BEFORE THINNING ==========" << std::endl;
167
    HepMC::Print::line(std::cout, *thinEvt);
168
    std::cout << "========== END EVENT BEFORE THINNING ==========" << std::endl;
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
  }

  /////////////////////////////////////////////////
  // Containers for manipulating particles/vertices
  /////////////////////////////////////////////////

  // Hadronization vertex must have at least two incoming partons to
  // conserve color. All outgoing must be hadrons; c.f. Pythia 2->1
  // color reconnection vertices.
  // Const iterators => must(?) save changes, then do them.
  // Always remove particles from vertices, then delete.
  // removePV:  remove particle from vertex
  // addinPV:   add particle in to vertex
  // addoutPV:  add particle out from vertex
  // removeV:   remove vertex from event
  // deleteP:   delete particle after all passes
  // deleteV:   delete vertex after all passes
  // HepMC ~GenVertex deletes particles, so remove them from ALL vertices

188
  typedef std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr> vpPair;
189
190
191
  std::vector<vpPair> removePV;
  std::vector<vpPair> addinPV;
  std::vector<vpPair> addoutPV;
192
  std::vector<HepMC::GenVertexPtr> removeV;
193
194

  // Use list with sort/unique to insure no multiple deletion
195
196
197
198
199
200
  std::list<HepMC::GenParticlePtr> deleteP;
  std::list<HepMC::GenVertexPtr> deleteV;
  std::list<HepMC::GenParticlePtr>::iterator dpItr;
  std::list<HepMC::GenParticlePtr>::iterator dpItrE;
  std::list<HepMC::GenVertexPtr>::iterator dvItr;
  std::list<HepMC::GenVertexPtr>::iterator dvItrE;
201
202
203
204
205

  //////////////////////////////
  // Find hadronization vertices
  //////////////////////////////

206
  if (doDebug) ATH_MSG_DEBUG("Find hadronization vertices");
207

208
  std::vector<HepMC::GenVertexPtr> hadVertices;
209

210
211
212
213
214
215
216
217
218
219
220
#ifdef HEPMC3
  for (auto hadv:thinEvt->vertices()) {
    if (!hadv) continue;
    if (hadv->particles_in().size() < 2) continue;
    if (hadv->particles_out().size() < 1) continue;
    // Check hadronization vertex
    // isHad is true if at least one hadron
    // q qbar -> pi is allowed, but q qbar -> W... is not
    bool isHadVtx = true;
    bool isHadOut = false;
    for (auto inp:  hadv->particles_in() ) {
221
      if (!isParton(inp)) { isHadVtx = false; break;}
222
223
224
225
226
227
228
229
230
231
    }
    for (auto vp:  hadv->particles_out()) {
      if (isParton(vp)) isHadVtx = false;
      if (isHadron(vp)) isHadOut = true;
    }
    isHadVtx = isHadVtx && isHadOut;
    if (isHadVtx) hadVertices.push_back(hadv);
    if (doDebug && isHadVtx) ATH_MSG_VERBOSE("Hadronization vertex " << HepMC::barcode(hadv));
  }
#else
232
233
234
235
  HepMC::GenEvent::vertex_iterator hadv = thinEvt->vertices_begin();
  HepMC::GenEvent::vertex_iterator hadvB = thinEvt->vertices_begin();
  HepMC::GenEvent::vertex_iterator hadvE = thinEvt->vertices_end();

236
237
238
239
  for (; hadv != hadvE; ++hadv) {
    if (!(*hadv)) continue;
    if ((*hadv)->particles_in_size() < 2) continue;
    if ((*hadv)->particles_out_size() < 1) continue;
240
241
242
243
244
245

    // Check hadronization vertex
    // isHad is true if at least one hadron
    // q qbar -> pi is allowed, but q qbar -> W... is not
    bool isHadVtx = true;
    bool isHadOut = false;
246
247
248
    HepMC::GenVertex::particles_in_const_iterator inp = (*hadv)->particles_in_const_begin();
    HepMC::GenVertex::particles_in_const_iterator inpE = (*hadv)->particles_in_const_end();
    for (; inp != inpE; ++inp) {
249
      if (!isParton(*inp)) { isHadVtx = false; break;}
250
    }
251
252
253
254
255
    HepMC::GenVertex::particles_out_const_iterator vp = (*hadv)->particles_out_const_begin();
    HepMC::GenVertex::particles_out_const_iterator vpE = (*hadv)->particles_out_const_end();
    for (; vp != vpE; ++vp) {
      if (isParton(*vp)) isHadVtx = false;
      if (isHadron(*vp)) isHadOut = true;
256
257
    }
    isHadVtx = isHadVtx && isHadOut;
258
259
    if (isHadVtx) hadVertices.push_back(*hadv);
    if (doDebug && isHadVtx) ATH_MSG_VERBOSE("Hadronization vertex " << HepMC::barcode(*hadv));
260
  }
261
#endif
262

263
  if (hadVertices.empty()) {
264
    ATH_MSG_WARNING("No hadronization vertices for event " << nEvent);
265
266
267
268
269
270
271
272
273
274
    ATH_MSG_WARNING("Exiting without changing event.");
    thinnedMcEvts->push_back(thinEvt);
    return StatusCode::SUCCESS;
  }

  //////////////////////////////////////////////////////////
  // Remove all incoming partons from hadronization vertices
  // Remove and delete all descendants
  //////////////////////////////////////////////////////////

275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
#ifdef HEPMC3
  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
    HepMC::GenVertexPtr ivtx = hadVertices[iv];
    if (doDebug) ATH_MSG_DEBUG("Removing partons from hadVertex " << HepMC::barcode(ivtx));
    for (auto  pin:  ivtx->particles_in()) {
      removePV.push_back(vpPair(ivtx, pin));
    }
  }
  // Remove all descendant particles. Will remove empty vertices later.
  // Might have parton decays of hadrons - hence delete sort/unique
  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
    HepMC::GenVertexPtr ivtx = hadVertices[iv]; 
    auto descendants = HepMC::descendant_particles(ivtx);
     for (auto  pout:  descendants) {
      HepMC::GenVertexPtr vpar = pout->production_vertex();
      if (vpar) removePV.push_back(vpPair(vpar, pout));
      HepMC::GenVertexPtr vend = pout->end_vertex();
      if (vend) removePV.push_back(vpPair(vend, pout));
      deleteP.push_back(pout);
    }
  }
#else
297
  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
298
    HepMC::GenVertex* ivtx = hadVertices[iv];
299
300
301
302
    if (doDebug) ATH_MSG_DEBUG("Removing partons from hadVertex " << HepMC::barcode(ivtx));
    HepMC::GenVertex::particles_in_const_iterator pin = ivtx->particles_in_const_begin();
    HepMC::GenVertex::particles_in_const_iterator pinE = ivtx->particles_in_const_end();
    for (; pin != pinE; ++pin) {
303
      removePV.emplace_back(ivtx, *pin);
304
305
306
307
308
    }
  }

  // Remove all descendant particles. Will remove empty vertices later.
  // Might have parton decays of hadrons - hence delete sort/unique
309
  for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
310
    HepMC::GenVertex* ivtx = hadVertices[iv];
311
312
313
    HepMC::GenVertex::particle_iterator pout = ivtx->particles_begin(HepMC::descendants);
    HepMC::GenVertex::particle_iterator poutE = ivtx->particles_end(HepMC::descendants);
    for (; pout != poutE; ++pout) {
314
      HepMC::GenVertex* vpar = (*pout)->production_vertex();
315
      if (vpar) removePV.emplace_back(vpar, *pout);
316
      HepMC::GenVertex* vend = (*pout)->end_vertex();
317
      if (vend) removePV.emplace_back(vend, *pout);
318
319
320
      deleteP.push_back(*pout);
    }
  }
321
#endif
322
323
324
325
326
327
328
329

  // Remove empty vertices
  // Remove all particles from Geant vertices
  // Remove and delete Geant vertices and particles
  // All Geant particles should have Geant parent vertex

  static const int cutG4 = 200000;

330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
#ifdef HEPMC3
  for (auto hadv:thinEvt->vertices()) {
    // Empth vertices
    if (hadv->particles_in().size() == 0 && hadv->particles_out().size() == 0) {
      removeV.push_back(hadv);
      deleteV.push_back(hadv);
    }
    // Geant vertices/particles
    if (HepMC::barcode(hadv) > -cutG4) continue;
    for (auto  pin: hadv->particles_in()) {
      removePV.push_back(vpPair(hadv, pin));
      if (HepMC::barcode(pin) > cutG4) { deleteP.push_back(pin); }
    }
    for (auto  pout: hadv->particles_out()) {
      removePV.push_back(vpPair(hadv, pout));
      if ( HepMC::barcode(pout) > cutG4) { deleteP.push_back(pout); }
    }
    removeV.push_back(hadv);
    deleteV.push_back(hadv);
  }
#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
356
357
358
359
      removeV.push_back(*hadv);
      deleteV.push_back(*hadv);
    }

    // Geant vertices/particles
360
361
362
363
    if (HepMC::barcode(*hadv) > -cutG4) continue;
    HepMC::GenVertex::particles_in_const_iterator pin = (*hadv)->particles_in_const_begin();
    HepMC::GenVertex::particles_in_const_iterator pinE = (*hadv)->particles_in_const_end();
    for (; pin != pinE; ++pin) {
364
      removePV.emplace_back(*hadv, *pin);
365
      if ((*pin)->barcode() > cutG4) { deleteP.push_back(*pin); }
366
    }
367
368
369
    HepMC::GenVertex::particles_out_const_iterator pout = (*hadv)->particles_out_const_begin();
    HepMC::GenVertex::particles_out_const_iterator poutE = (*hadv)->particles_out_const_end();
    for (; pout != poutE; ++pout) {
370
      removePV.emplace_back(*hadv, *pout);
371
      if ((*pout)->barcode() > cutG4) { deleteP.push_back(*pout); }
372
373
374
375
    }
    removeV.push_back(*hadv);
    deleteV.push_back(*hadv);
  }
376
#endif
377
378
379

  // Actually implement changes

380
  for (unsigned int i = 0; i < removePV.size(); ++i) {
381
382
383
384
385
386
    HepMC::GenVertexPtr v = removePV[i].first;
    HepMC::GenParticlePtr p = removePV[i].second;
#ifdef HEPMC3
    v->remove_particle_in(p);
    v->remove_particle_out(p);
#else
387
    v->remove_particle(p);
388
#endif
389
390
  }

391
  for (unsigned int i = 0; i < addoutPV.size(); ++i) {
392
393
    HepMC::GenVertexPtr v = addoutPV[i].first;
    HepMC::GenParticlePtr p = addoutPV[i].second;
394
395
    v->add_particle_out(p);
  }
396
397
398
399
400
401
402
403
404
405
#ifdef HEPMC3
  for (unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
    HepMC::GenVertexPtr v = hadVertices[iv];
    if (v->particles_in().size() != 0 || v->particles_out().size() != 0) {
      ATH_MSG_WARNING("Removing vertex " << HepMC::barcode(v) << " for event " << nEvent << " with in/out particles " << v->particles_in().size() << " " << v->particles_out().size());
    }
    thinEvt->remove_vertex(hadVertices[iv]);
  }
//AV: HepMC3 uses smart pointers
#else
406

407
  for (unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
408
    HepMC::GenVertex* v = hadVertices[iv];
409
410
    if (v->particles_in_size() != 0 || v->particles_out_size() != 0) {
      ATH_MSG_WARNING("Removing vertex " << v->barcode() << " for event " << nEvent << " with in/out particles " << v->particles_in_size() << " " << v->particles_out_size());
411
    }
412
    if (!thinEvt->remove_vertex(hadVertices[iv])) { ATH_MSG_WARNING("Error removing vertex " << v->barcode() << " for event " << nEvent); }
413
414
415
416
  }

  // Delete removed particles/vertices

417
  if (doDebug) ATH_MSG_DEBUG("Deleting hadronization vertices " << deleteV.size());
418
419
  deleteV.sort();
  deleteV.unique();
420
421
422
  for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
    if (doDebug) ATH_MSG_VERBOSE("Deleting vertex " << (*dvItr)->barcode());
    if (*dvItr) delete (*dvItr);
423
424
425
426
  }

  deleteP.sort();
  deleteP.unique();
427
428
429
  for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
    if (doDebug) ATH_MSG_VERBOSE("Deleting particle " << (*dpItr)->barcode());
    if (*dpItr) delete (*dpItr);
430
  }
431
#endif
432
433
434
435
436

  ////////////////////////
  // Cluster final partons
  ////////////////////////

437
438
  if (doDebug && doExtra) {
    std::cout << "========== BEGIN EVENT BEFORE CLUSTER ==========" << std::endl;
439
    HepMC::Print::line(std::cout, *thinEvt);
440
    std::cout << "========== END EVENT BEFORE CLUSTER ==========" << std::endl;
441
442
443
444
445
446
447
448
449
450
451
452
453
454
  }

  // Possible cases:
  // 1->1 branching:  Drop earlier 1 and vertex.
  // 2->1 connection: Drop 1 and vertex; leave 2 free
  //                  Note momentum is conserved for these.
  //                  Do not split 2->1 from initial 2->1->2 in Herwig.
  // 1->2 branching:  Set p1 = p2a+p2b, drop 2 and vertex.
  //                  Add 1 to first vertex.

  // Preserving color connections required merging vertices. Now no
  // longer needed.

  bool moreP = true;
455
  using vpPair = std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr>;
456
457
458
459
460
461
462
  removePV.clear();
  addinPV.clear();
  addoutPV.clear();
  removeV.clear();
  deleteP.clear();
  deleteV.clear();

463
  using pkPair = std::pair<HepMC::GenParticlePtr, HepMC::FourVector>;
464
465
  std::vector<pkPair> changePK;

466
467
468
469
  if (doDebug) ATH_MSG_DEBUG("Start parton thinning");
  while (moreP) {
#ifdef HEPMC3
    if (doDebug) ATH_MSG_DEBUG("New parton pass " << inEvent << " " << thinEvt->particles().size() << " " << thinEvt->vertices().size());
470
#else
471
    if (doDebug) ATH_MSG_DEBUG("New parton pass " << inEvent << " " << thinEvt->particles_size() << " " << thinEvt->vertices_size());
472
#endif
473
474
475
476
477
478
479

    moreP = false;
    removePV.clear();
    addinPV.clear();
    addoutPV.clear();
    removeV.clear();
    changePK.clear();
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
#ifdef HEPMC3
    // Find final partons
    for (auto fp: thinEvt->particles() ) {
      int iCase = 0;
      if (!isFinalParton(fp)) continue;
      if (doDebug) ATH_MSG_DEBUG("Starting final parton " << HepMC::barcode(fp));
      // Production/end vertices
      HepMC::GenVertexPtr pvtx = fp->production_vertex();
      if (!pvtx) {
        ATH_MSG_WARNING("Missing production for final parton " << HepMC::barcode(fp));
        continue;
      }
      if (doDebug) ATH_MSG_DEBUG("Final parton " << HepMC::barcode(pvtx) << " " << HepMC::barcode(fp));
      ////////////
      // Case 1->1
      ////////////
      // One-particle decay; use final particle
      // ppvtx -> pp -> pvtx -> fp
      if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 1) {
        // Incoming particle to parent vertex
        auto pitr = pvtx->particles_in().begin();
        HepMC::GenParticlePtr pp = *pitr;
        if (!pp || HepMC::barcode(pp) == 0) {
          ATH_MSG_DEBUG("1->1: missing pp for fp " << HepMC::barcode(fp));
          ++m_missCount;
          continue;
        }
        // Its parent vertex
        HepMC::GenVertexPtr ppvtx = pp->production_vertex();
        if (!ppvtx || HepMC::barcode(ppvtx) == 0) {
          ATH_MSG_DEBUG("1->1: missing ppvtx for fp " << HepMC::barcode(fp));
          ++m_missCount;
          continue;
        }
        moreP = true;
        iCase = 1;
        removePV.push_back(vpPair(ppvtx, pp));
        removePV.push_back(vpPair(pvtx, pp));
        deleteP.push_back(pp);
        removeV.push_back(pvtx);
        deleteV.push_back(pvtx);
        addoutPV.push_back(vpPair(ppvtx, fp));
        if (doDebug) { ATH_MSG_DEBUG("1->1: ppvtx,pp,pvtx,fp,evtx " << HepMC::barcode(ppvtx) << " " << HepMC::barcode(pp) << " " << HepMC::barcode(pvtx) << " " << HepMC::barcode(fp)); }
      }
      ////////////
      // Case 2->1
      ////////////
      // Color recombination. Momentum is conserved so just keep 2.
      // Drop 1 and vertex.
      // ppvtx1,ppvtx2 -> pp1,pp2 -> pvtx -> fp
      // Recombination should not affect hard physics!
      if (pvtx->particles_in().size() == 2 && pvtx->particles_out().size() == 1) {
        // Incoming particles to parent vertex
        auto pitr = pvtx->particles_in().begin();
        HepMC::GenParticlePtr pp1 = *pitr;
        ++pitr;
        HepMC::GenParticlePtr pp2 = *pitr;
        // Check for 2->1->2 initial state interactions in Herwig++
        // Initial partons have pt=0, use pt<0.001MeV
        if (std::abs(pp1->momentum().perp()) < 1.e-3) continue;
        if (std::abs(pp2->momentum().perp()) < 1.e-3) continue;
        // Their parent vertices
        HepMC::GenVertexPtr ppvtx1 = pp1->production_vertex();
        HepMC::GenVertexPtr ppvtx2 = pp2->production_vertex();
        if (!ppvtx1 || HepMC::barcode(ppvtx1) == 0) {
          ATH_MSG_DEBUG("2->1: missing ppvtx1 for fp " << HepMC::barcode(fp));
          ++m_missCount;
          continue;
        }
        if (!ppvtx2 || HepMC::barcode(ppvtx2) == 0) {
          ATH_MSG_DEBUG("2->1: missing ppvtx2 for fp " << HepMC::barcode(fp));
          ++m_missCount;
          continue;
        }
        moreP = true;
        iCase = 2;
        removePV.push_back(vpPair(pvtx, fp));
        removePV.push_back(vpPair(pvtx, pp1));
        removePV.push_back(vpPair(pvtx, pp2));
        deleteP.push_back(fp);
        removeV.push_back(pvtx);
        deleteV.push_back(pvtx);
        if (doDebug) {
          ATH_MSG_DEBUG("2->1: ppvtx1,pp1,ppvtx2,pp2,pvtx,fp " << HepMC::barcode(ppvtx1) << " " << HepMC::barcode(pp1) 
          << " " << HepMC::barcode(ppvtx2) << " " << HepMC::barcode(pp2) << " " << HepMC::barcode(pvtx) << " "<< HepMC::barcode(fp));
        }
      }
      ////////////
      // Case 1->2
      ////////////
      // Parton branching. Momentum not conserved; 2 momenta correct
      // Drop only if mass is below cut
      // ppvtx -> pp -> pvtx -> pout1,pout2/fp
      if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 2) {
        auto poutitr = pvtx->particles_out().begin();
        HepMC::GenParticlePtr pout1 = *poutitr;
        ++poutitr;
        HepMC::GenParticlePtr pout2 = *poutitr;
        // Require two final partons and avoid duplication
        if (fp == pout1) {
          if (!isFinalParton(pout2)) {
            if (doDebug) ATH_MSG_DEBUG("1->2: not final " << HepMC::barcode(pout2));
            continue;
          }
        } else if (fp == pout2) {
          if (!isFinalParton(pout1)) {
            if (doDebug) ATH_MSG_DEBUG("1->2: not final " << HepMC::barcode(pout1));
            continue;
          }
        } else {
          ATH_MSG_WARNING("1->2: No match found for branching " << HepMC::barcode(fp) << " " << HepMC::barcode(pvtx) << " " << HepMC::barcode(pout1) << " " << HepMC::barcode(pout2));
          continue;
        }
        if (fp != pout1) continue;
        // Incoming particle
        auto pitr = pvtx->particles_in().begin();
        HepMC::GenParticlePtr pp = *pitr;
        // Do not merge initial partons (pt<1MeV or m<-1MeV)
        if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0) continue;
        if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0) continue;
        // Parton pair mass cut
        HepMC::FourVector p12 = vtxOutMom(pvtx);
        double m12 = p12.m();
        if (m12 < 0) {
          if (std::abs(m12) > 10. + 1.0e-5 * p12.e()) {
            ATH_MSG_WARNING("Spacelike mass^2 for parton sum " << m12 << " " << HepMC::barcode(pp) << " " << HepMC::barcode(pvtx) << " " << HepMC::barcode(pout1) << " " << HepMC::barcode(pout2));
          }
          m12 = 0;
        }
        if (doDebug) ATH_MSG_DEBUG("1->2: parton pair mass " << m12);
        // If mass > cut, keep pair
        if (m12 > m_partonCut) {
          if (doDebug) ATH_MSG_DEBUG("Keeping 1->2: parton mass " << m12);
          continue;
        }
        // Associated vertices
        HepMC::GenVertexPtr ppvtx = pp->production_vertex();
        if (!ppvtx || HepMC::barcode(ppvtx) == 0) {
          ATH_MSG_DEBUG("1->2: missing ppvtx for fp " << HepMC::barcode(fp));
          ++m_missCount;
          continue;
        }
        // Merge branching
        moreP = true;
        iCase = 3;
        if (doDebug) ATH_MSG_DEBUG("Merging 1->2: mass " << p12.m());
        changePK.push_back(pkPair(pp, p12));
        removePV.push_back(vpPair(pvtx, pp));
        removePV.push_back(vpPair(pvtx, pout1));
        removePV.push_back(vpPair(pvtx, pout2));
        deleteP.push_back(pout1);
        deleteP.push_back(pout2);
        removeV.push_back(pvtx);
        deleteV.push_back(pvtx);
        if (doDebug) {
          ATH_MSG_DEBUG("Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << HepMC::barcode(ppvtx) << " " << HepMC::barcode(pp) << " " << HepMC::barcode(pvtx) << " " << HepMC::barcode(pout1) << " "
                                                                      << HepMC::barcode(pout2));
          ATH_MSG_DEBUG("Merge 1->2: id " << pp->pdg_id() << " " << pout1->pdg_id() << " " << pout2->pdg_id());
        }
      } // end 1->2 case
      /////////////////////////
      // Incoming proton vertex
      /////////////////////////
      // Do nothing
      if (pvtx->particles_in().size() == 1) {
        // Incoming particle to parent vertex
        auto pitr = pvtx->particles_in().begin();
        HepMC::GenParticlePtr pp = *pitr;
        if (std::abs(pp->pdg_id()) == 2212) iCase = -1;
      }
      // Case not found
      // Need test for 2->2 in underlying event
      if (iCase == 0) {
        if (doDebug) ATH_MSG_DEBUG("Case not found " << HepMC::barcode(pvtx) << " " << HepMC::barcode(fp) << " " << pvtx->particles_in().size() << " " << pvtx->particles_out().size());
      }
    } // end final parton loop
#else
657
658
659
660
661
662

    HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
    HepMC::GenEvent::particle_iterator finpB = thinEvt->particles_begin();
    HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();

    // Find final partons
663
    for (finp = finpB; finp != finpE; ++finp) {
664
665
666
      int iCase = 0;

      HepMC::GenParticle* fp = *finp;
667
668
      if (!isFinalParton(fp)) continue;
      if (doDebug) ATH_MSG_DEBUG("Starting final parton " << fp->barcode());
669
670
671

      // Production/end vertices
      HepMC::GenVertex* pvtx = fp->production_vertex();
672
673
      if (!pvtx) {
        ATH_MSG_WARNING("Missing production for final parton " << fp->barcode());
674
675
        continue;
      }
676
      if (doDebug) ATH_MSG_DEBUG("Final parton " << pvtx->barcode() << " " << fp->barcode());
677
678
679
680
681
682
683
684

      ////////////
      // Case 1->1
      ////////////

      // One-particle decay; use final particle
      // ppvtx -> pp -> pvtx -> fp

685
      if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 1) {
686
        // Incoming particle to parent vertex
687
        HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
688
        HepMC::GenParticle* pp = *pitr;
689
690
        if (!pp || pp->barcode() == 0) {
          ATH_MSG_DEBUG("1->1: missing pp for fp " << fp->barcode());
691
692
693
694
          ++m_missCount;
          continue;
        }
        // Its parent vertex
695
696
697
        HepMC::GenVertex* ppvtx = pp->production_vertex();
        if (!ppvtx || ppvtx->barcode() == 0) {
          ATH_MSG_DEBUG("1->1: missing ppvtx for fp " << fp->barcode());
698
          ++m_missCount;
699
          continue;
700
701
702
703
        }
        moreP = true;
        iCase = 1;

704
705
        removePV.emplace_back(ppvtx, pp);
        removePV.emplace_back(pvtx, pp);
706
707
708
        deleteP.push_back(pp);
        removeV.push_back(pvtx);
        deleteV.push_back(pvtx);
709
        addoutPV.emplace_back(ppvtx, fp);
710
        if (doDebug) { ATH_MSG_DEBUG("1->1: ppvtx,pp,pvtx,fp,evtx " << ppvtx->barcode() << " " << pp->barcode() << " " << pvtx->barcode() << " " << fp->barcode()); }
711
712
713
714
715
716
717
      }

      ////////////
      // Case 2->1
      ////////////

      // Color recombination. Momentum is conserved so just keep 2.
718
      // Drop 1 and vertex.
719
720
721
      // ppvtx1,ppvtx2 -> pp1,pp2 -> pvtx -> fp
      // Recombination should not affect hard physics!

722
      if (pvtx->particles_in_size() == 2 && pvtx->particles_out_size() == 1) {
723
        // Incoming particles to parent vertex
724
        HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
725
726
727
728
729
730
        HepMC::GenParticle* pp1 = *pitr;
        ++pitr;
        HepMC::GenParticle* pp2 = *pitr;

        // Check for 2->1->2 initial state interactions in Herwig++
        // Initial partons have pt=0, use pt<0.001MeV
731
732
        if (std::abs(pp1->momentum().perp()) < 1.e-3) continue;
        if (std::abs(pp2->momentum().perp()) < 1.e-3) continue;
733
        // Their parent vertices
734
735
736
737
        HepMC::GenVertex* ppvtx1 = pp1->production_vertex();
        HepMC::GenVertex* ppvtx2 = pp2->production_vertex();
        if (!ppvtx1 || ppvtx1->barcode() == 0) {
          ATH_MSG_DEBUG("2->1: missing ppvtx1 for fp " << fp->barcode());
738
739
740
          ++m_missCount;
          continue;
        }
741
742
        if (!ppvtx2 || ppvtx2->barcode() == 0) {
          ATH_MSG_DEBUG("2->1: missing ppvtx2 for fp " << fp->barcode());
743
744
745
          ++m_missCount;
          continue;
        }
746

747
748
749
        moreP = true;
        iCase = 2;

750
751
752
        removePV.emplace_back(pvtx, fp);
        removePV.emplace_back(pvtx, pp1);
        removePV.emplace_back(pvtx, pp2);
753
754
755
756
        deleteP.push_back(fp);
        removeV.push_back(pvtx);
        deleteV.push_back(pvtx);

757
758
759
        if (doDebug) {
          ATH_MSG_DEBUG("2->1: ppvtx1,pp1,ppvtx2,pp2,pvtx,fp " << ppvtx1->barcode() << " " << pp1->barcode() << " " << ppvtx2->barcode() << " " << pp2->barcode() << " " << pvtx->barcode() << " "
                                                               << fp->barcode());
760
761
762
763
764
765
766
767
768
769
770
        }
      }

      ////////////
      // Case 1->2
      ////////////

      // Parton branching. Momentum not conserved; 2 momenta correct
      // Drop only if mass is below cut
      // ppvtx -> pp -> pvtx -> pout1,pout2/fp

771
772
      if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 2) {
        HepMC::GenVertex::particles_out_const_iterator poutitr = pvtx->particles_out_const_begin();
773
774
775
776
777
        HepMC::GenParticle* pout1 = *poutitr;
        ++poutitr;
        HepMC::GenParticle* pout2 = *poutitr;

        // Require two final partons and avoid duplication
778
779
780
        if (fp == pout1) {
          if (!isFinalParton(pout2)) {
            if (doDebug) ATH_MSG_DEBUG("1->2: not final " << pout2->barcode());
781
782
            continue;
          }
783
784
785
        } else if (fp == pout2) {
          if (!isFinalParton(pout1)) {
            if (doDebug) ATH_MSG_DEBUG("1->2: not final " << pout1->barcode());
786
787
788
            continue;
          }
        } else {
789
          ATH_MSG_WARNING("1->2: No match found for branching " << fp->barcode() << " " << pvtx->barcode() << " " << pout1->barcode() << " " << pout2->barcode());
790
791
          continue;
        }
792
793
794
        if (fp != pout1) continue;
        // Incoming particle
        HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
795
796
797
        HepMC::GenParticle* pp = *pitr;

        // Do not merge initial partons (pt<1MeV or m<-1MeV)
798
799
        if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0) continue;
        if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0) continue;
800
801
802
803

        // Parton pair mass cut
        HepMC::FourVector p12 = vtxOutMom(pvtx);
        double m12 = p12.m();
804
805
806
        if (m12 < 0) {
          if (fabs(m12) > 10. + 1.0e-5 * p12.e()) {
            ATH_MSG_WARNING("Spacelike mass^2 for parton sum " << m12 << " " << pp->barcode() << " " << pvtx->barcode() << " " << pout1->barcode() << " " << pout2->barcode());
807
808
809
          }
          m12 = 0;
        }
810
        if (doDebug) ATH_MSG_DEBUG("1->2: parton pair mass " << m12);
811
        // If mass > cut, keep pair
812
813
        if (m12 > m_partonCut) {
          if (doDebug) ATH_MSG_DEBUG("Keeping 1->2: parton mass " << m12);
814
815
816
817
818
          continue;
        }

        // Associated vertices
        HepMC::GenVertex* ppvtx = pp->production_vertex();
819
820
        if (!ppvtx || ppvtx->barcode() == 0) {
          ATH_MSG_DEBUG("1->2: missing ppvtx for fp " << fp->barcode());
821
822
823
824
825
826
827
          ++m_missCount;
          continue;
        }

        // Merge branching
        moreP = true;
        iCase = 3;
828
        if (doDebug) ATH_MSG_DEBUG("Merging 1->2: mass " << p12.m());
829

830
831
832
833
        changePK.emplace_back(pp, p12);
        removePV.emplace_back(pvtx, pp);
        removePV.emplace_back(pvtx, pout1);
        removePV.emplace_back(pvtx, pout2);
834
835
836
837
838
839

        deleteP.push_back(pout1);
        deleteP.push_back(pout2);
        removeV.push_back(pvtx);
        deleteV.push_back(pvtx);

840
841
842
843
        if (doDebug) {
          ATH_MSG_DEBUG("Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << HepMC::barcode(ppvtx) << " " << HepMC::barcode(pp) << " " << HepMC::barcode(pvtx) << " " << HepMC::barcode(pout1) << " "
                                                                      << HepMC::barcode(pout2));
          ATH_MSG_DEBUG("Merge 1->2: id " << pp->pdg_id() << " " << pout1->pdg_id() << " " << pout2->pdg_id());
844
        }
845
      } // end 1->2 case
846
847
848
849
850
851

      /////////////////////////
      // Incoming proton vertex
      /////////////////////////

      // Do nothing
852
      if (pvtx->particles_in_size() == 1) {
853
        // Incoming particle to parent vertex
854
        HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
855
        HepMC::GenParticle* pp = *pitr;
856
        if (std::abs(pp->pdg_id()) == 2212) iCase = -1;
857
858
859
860
      }

      // Case not found
      // Need test for 2->2 in underlying event
861
862
      if (iCase == 0) {
        if (doDebug) ATH_MSG_DEBUG("Case not found " << pvtx->barcode() << " " << fp->barcode() << " " << pvtx->particles_in_size() << " " << pvtx->particles_out_size());
863
864
      }

865
    } // end final parton loop
866
#endif
867
868
869

    // Actually implement changes -- remove particles from vertices
    // Parton ends free, so no addinPV
870
    if (doDebug) ATH_MSG_DEBUG("Actually removing particles " << removePV.size());
871

872
    for (unsigned int i = 0; i < removePV.size(); ++i) {
873
874
875
876
877
878
879
      HepMC::GenVertexPtr v = removePV[i].first;
      HepMC::GenParticlePtr p = removePV[i].second;
      if (doDebug) ATH_MSG_VERBOSE("Removing v,p " << HepMC::barcode(v) << " " << HepMC::barcode(p));
#ifdef HEPMC3
      v->remove_particle_in(p);
      v->remove_particle_out(p);
#else      
880
      v->remove_particle(p);
881
#endif
882
883
884
    }

    // Actually implement changes -- add particles to vertices
885
886
    if (doDebug) ATH_MSG_DEBUG("Actually add particles in/out " << addinPV.size() << " " << addoutPV.size());
    for (unsigned int i = 0; i < addoutPV.size(); ++i) {
887
888
889
      HepMC::GenVertexPtr v = addoutPV[i].first;
      HepMC::GenParticlePtr p = addoutPV[i].second;
      if (doDebug) ATH_MSG_VERBOSE("Adding v,p " << HepMC::barcode(v) << " " << HepMC::barcode(p));
890
891
892
893
      v->add_particle_out(p);
    }

    // Actually implement changes -- change momenta
894
    for (unsigned int i = 0; i < changePK.size(); ++i) {
895
      HepMC::GenParticlePtr pp = changePK[i].first;
896
      //! float eold = pp->momentum().e();
897
898
899
900
      pp->set_momentum(changePK[i].second);
    }

    // Actually implement changes -- remove vertices
901
902
    if (doDebug) ATH_MSG_DEBUG("Actually remove vertices " << removeV.size());
    for (unsigned int i = 0; i < removeV.size(); ++i) {
903
904
905
906
907
#ifdef HEPMC3
      int nv = thinEvt->vertices().size();
      if (doDebug) { ATH_MSG_VERBOSE("Removing vertex " << HepMC::barcode(removeV[i]) << " " << nv << " " << thinEvt->vertices().size()); }
      thinEvt->remove_vertex(removeV[i]);
#else
908
      int nv = thinEvt->vertices_size();
909
910
      if (thinEvt->remove_vertex(removeV[i])) {
        if (doDebug) { ATH_MSG_VERBOSE("Removed vertex " << removeV[i]->barcode() << " " << nv << " " << thinEvt->vertices_size()); }
911
      } else {
912
        ATH_MSG_WARNING("Failed to remove vertex " << removeV[i]->barcode());
913
      }
914
#endif
915
    }
916
    if (doDebug) ATH_MSG_DEBUG("End while(moreP) pass " << moreP);
917

918
  } // end moreP
919

920
921
922
#ifdef HEPMC3
//AV HepMC3 uses smartpointers. This part is not needed.
#else
923
  // Delete removed particles/vertices
924
  if (doDebug) ATH_MSG_DEBUG("Deleting vertices " << deleteV.size());
925
926
  deleteV.sort();
  deleteV.unique();
927
  for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
928
    if (doDebug) ATH_MSG_VERBOSE("Deleting vertex " << HepMC::barcode((*dvItr)));
929
    if (*dvItr) delete (*dvItr);
930
931
  }

932
  if (doDebug) ATH_MSG_DEBUG("Deleting particles " << deleteP.size());
933
934
  deleteP.sort();
  deleteP.unique();
935
  for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
936
    if (doDebug) ATH_MSG_VERBOSE("Deleting particle " << HepMC::barcode((*dpItr)));
937
    if (*dpItr) delete (*dpItr);
938
939
  }

940
#endif
941
942
943
944
  //////////////////////////////
  // Strip soft underlying stuff
  //////////////////////////////

945
946
  if (doDebug && doExtra) {
    std::cout << "========== BEGIN EVENT BEFORE SOFT ==========" << std::endl;
947
    HepMC::Print::line(std::cout, *thinEvt);
948
    std::cout << "========== END EVENT BEFORE SOFT ==========" << std::endl;
949
950
  }

951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
  // Have deleted all hadronization particles
  // Find all particles connected to hard process(es) with m_T>10GeV
  std::list<HepMC::GenParticlePtr> pNotHad;
  std::list<HepMC::GenParticlePtr> pHard;
#ifdef HEPMC3
  std::vector<HepMC::GenParticlePtr> beams=thinEvt->beams();
  for (auto fp: thinEvt->particles()) {
     HepMC::GenVertexPtr pvtx = fp->production_vertex();
    if (!pvtx) continue;

    double ep = fp->momentum().e();
    double pzp = fp->momentum().pz();
    double mtmax = (ep + pzp) * (ep - pzp);
    auto ancestors = HepMC::ancestor_particles(fp);
  
    for (auto gpar: ancestors) {
      double e = gpar->momentum().e();
      double pz = gpar->momentum().pz();
969
      mtmax = std::max((e+pz)*(e-pz), mtmax);
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
    }

    // Keep hard particles and all ancestors
    pNotHad.push_back(fp);
    int ida = std::abs(fp->pdg_id());
    bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
    if (mtmax > m_hardCut * m_hardCut || keepid) {
      pHard.push_back(fp);
      pHard.insert(pHard.end(),ancestors.begin(),ancestors.end());
    }    

    // Also keep all descendants of interesting particles
    // Include leptons to get photons in Sherpa with no Z parent
    // All hard descendants would include soft initial radiation
    // Will remove duplicates with list sort/unique
    bool keepid2 = ( ida == 6) || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
    if (keepid2 && fp->end_vertex()) {
      auto descendants = HepMC::descendant_particles(fp);
      pHard.insert(pHard.end(),descendants.begin(),descendants.end());
    }
  }

  // Sort/unique lists
  pNotHad.sort();
  pNotHad.unique();
  pHard.sort();
  pHard.unique();

  // Dump information
  if (doDebug) {
     int nhard = 0;
    for (auto ph: pHard) {
      ++nhard;
      ATH_MSG_DEBUG("Hard GenParticles " << HepMC::barcode(ph) << " " << ph->pdg_id() << " " << ph->momentum().perp() / 1000. << " " << ph->momentum().pz() / 1000.);
    }
    if (doDebug) ATH_MSG_DEBUG("Hard GenParticles total " << nhard);
  }

#else
//AV: This algorithm is terribly slow. For each particle the descendants and ancestors are called.
1010
1011
1012
  HepMC::GenParticle* beams[2];
  beams[0] = thinEvt->beam_particles().first;
  beams[1] = thinEvt->beam_particles().second;
1013

1014
1015
1016
  HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
  HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();

1017
1018
1019
1020
  for (; finp != finpE; ++finp) {
    HepMC::GenParticle* fp = *finp;
    HepMC::GenVertex* pvtx = fp->production_vertex();
    if (!pvtx) continue;
1021
1022
1023

    double ep = fp->momentum().e();
    double pzp = fp->momentum().pz();
1024
1025
    double mtmax = (ep + pzp) * (ep - pzp);
    HepMC::GenVertex::particle_iterator gpar = fp->production_vertex()->particles_begin(HepMC::ancestors);
1026
    HepMC::GenVertex::particle_iterator gparB = gpar;
1027
    HepMC::GenVertex::particle_iterator gparE = fp->production_vertex()->particles_end(HepMC::ancestors);
1028

1029
    for (; gpar != gparE; ++gpar) {
1030
1031
      double e = (*gpar)->momentum().e();
      double pz = (*gpar)->momentum().pz();
1032
      mtmax = std::max((e+pz)*(e-pz), mtmax);
1033
1034
1035
1036
    }

    // Keep hard particles and all ancestors
    pNotHad.push_back(fp);
1037
    int ida = std::abs(fp->pdg_id());
1038
1039
    bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
    if (mtmax > m_hardCut * m_hardCut || keepid) {
1040
      pHard.push_back(fp);
1041
1042
      for (gpar = gparB; gpar != gparE; ++gpar)
        pHard.push_back(*gpar);
1043
1044
1045
1046
1047
1048
    }

    // Also keep all descendants of interesting particles
    // Include leptons to get photons in Sherpa with no Z parent
    // All hard descendants would include soft initial radiation
    // Will remove duplicates with list sort/unique
1049
1050
1051
1052
1053
1054
    bool keepid2 = ida == 6 || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
    if (keepid2 && fp->end_vertex()) {
      HepMC::GenVertex::particle_iterator des = fp->end_vertex()->particles_begin(HepMC::descendants);
      HepMC::GenVertex::particle_iterator desE = fp->end_vertex()->particles_end(HepMC::descendants);
      for (; des != desE; ++des)
        pHard.push_back(*des);
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
    }
  }

  // Sort/unique lists
  pNotHad.sort();
  pNotHad.unique();
  pHard.sort();
  pHard.unique();

  // Dump information
1065
  if (doDebug) {
1066
1067
1068
    std::list<HepMC::GenParticle*>::iterator hItr2 = pHard.begin();
    std::list<HepMC::GenParticle*>::iterator hItr2E = pHard.end();
    int nhard = 0;
1069
    for (; hItr2 != hItr2E; ++hItr2) {
1070
      ++nhard;
1071
      ATH_MSG_DEBUG("Hard GenParticles " << (*hItr2)->barcode() << " " << (*hItr2)->pdg_id() << " " << (*hItr2)->momentum().perp() / 1000. << " " << (*hItr2)->momentum().pz() / 1000.);
1072
    }
1073
    if (doDebug) ATH_MSG_DEBUG("Hard GenParticles total " << nhard);
1074
1075
  }

1076
#endif
1077
1078
1079
  // Remove non-hadronization, non-hard GenParticles from vertices
  // and delete them using lists constructed above.
  // Any 1->1 vertices created will be removed below.
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
#ifdef HEPMC3
  for (auto p: pNotHad) {
    // Skip hard ones
    bool isHard = false;
    for (auto h: pHard) {
      if (p == h) {
        isHard = true;
        break;
      }
    }
    if (doDebug) ATH_MSG_DEBUG("Particle bc/isHard " << HepMC::barcode(p) << " " << isHard);
    if (isHard) continue;
    HepMC::GenVertexPtr pvtx = p->production_vertex();
    if (pvtx) pvtx->remove_particle_out(p);
    HepMC::GenVertexPtr evtx = p->end_vertex();
    if (evtx) evtx->remove_particle_in(p);
   }
#else
1098
1099
1100
1101
1102
1103
1104
1105

  std::list<HepMC::GenParticle*>::iterator pItr = pNotHad.begin();
  std::list<HepMC::GenParticle*>::iterator pItrE = pNotHad.end();

  std::list<HepMC::GenParticle*>::iterator hItr = pHard.begin();
  std::list<HepMC::GenParticle*>::iterator hItrB = pHard.begin();
  std::list<HepMC::GenParticle*>::iterator hItrE = pHard.end();

1106
  for (; pItr != pItrE; ++pItr) {
1107
1108
1109
1110
    HepMC::GenParticle* p = *pItr;

    // Skip hard ones
    bool isHard = false;
1111
1112
    for (hItr = hItrB; hItr != hItrE; ++hItr) {
      if (p == (*hItr)) {
1113
1114
1115
1116
        isHard = true;
        break;
      }
    }
1117
1118
    if (doDebug) ATH_MSG_DEBUG("Particle bc/isHard " << p->barcode() << " " << isHard);
    if (isHard) continue;
1119
    HepMC::GenVertex* pvtx = p->production_vertex();
1120
    if (pvtx) pvtx->remove_particle(p);
1121
    HepMC::GenVertex* evtx = p->end_vertex();
1122
    if (evtx) evtx->remove_particle(p);
1123
1124
    delete p;
  }
1125
#endif
1126
1127
1128
1129
1130
1131
1132
1133

  /////////////////////////////////////////////////////////
  // Remove and delete vertices with no remaining particles
  /////////////////////////////////////////////////////////

  removeV.clear();
  deleteV.clear();

1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
#ifdef HEPMC3
  for (auto vtx: thinEvt->vertices()) {
    if (vtx->particles_in().size() != 0) continue;
    if (vtx->particles_out().size() != 0) continue;
    removeV.push_back(vtx);
    deleteV.push_back(vtx);
  }
  if (doDebug) ATH_MSG_DEBUG("Removing/deleting 0-particle vertices " << removeV.size() << " " << deleteV.size());
  for (unsigned int i = 0; i < removeV.size(); ++i) {
    if (doDebug) ATH_MSG_VERBOSE("Removing vertex " << HepMC::barcode(removeV[i]));
    thinEvt->remove_vertex(removeV[i]);
  }
#else
1147
1148
  HepMC::GenEvent::vertex_iterator vtx = thinEvt->vertices_begin();
  HepMC::GenEvent::vertex_iterator vtxE = thinEvt->vertices_end();
1149
1150
1151
  for (; vtx != vtxE; ++vtx) {
    if ((*vtx)->particles_in_size() != 0) continue;
    if ((*vtx)->particles_out_size() != 0) continue;
1152
1153
1154
1155
    removeV.push_back(*vtx);
    deleteV.push_back(*vtx);
  }

1156
1157
1158
1159
  if (doDebug) ATH_MSG_DEBUG("Removing/deleting 0-particle vertices " << removeV.size() << " " << deleteV.size());
  for (unsigned int i = 0; i < removeV.size(); ++i) {
    if (thinEvt->remove_vertex(removeV[i])) {
      if (doDebug) ATH_MSG_VERBOSE("Removed vertex " << removeV[i]->barcode());
1160
    } else {
1161
      ATH_MSG_WARNING("Failed to remove vertex " << removeV[i]->barcode());
1162
1163
1164
1165
1166
    }
  }

  deleteV.sort();
  deleteV.unique();
1167
1168
1169
  for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
    if (doDebug) ATH_MSG_VERBOSE("Deleting vertex " << (*dvItr)->barcode());
    if (*dvItr) delete (*dvItr);
1170
  }
1171
#endif
1172
1173
1174
1175
1176

  ////////////////////////////////
  // Remove remaining 1-1 vertices
  ////////////////////////////////

1177
1178
  if (doDebug && doExtra) {
    std::cout << "========== BEGIN EVENT BEFORE 1-BODY ==========" << std::endl;
1179
    HepMC::Print::line(std::cout, *thinEvt);
1180
    std::cout << "========== END EVENT BEFORE 1-BODY ==========" << std::endl;
1181
1182
1183
1184
1185
1186
  }

  // Not clear how to order sweep, so do it one at a time.
  // Always keep child
  // pvtx -> pin -> vtx11 -> pout -> evtx

1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
#ifdef HEPMC3
  bool moreV1 = true;
  HepMC::GenVertexPtr vtx11;
  HepMC::</