104 {
105
107
109
110
111
112
113
114
116 const bool doDebug = false;
117 const bool doExtra = false;
118
119
120
121
125 return StatusCode::FAILURE;
126 }
127
128 if (mcEvts->
empty()) {
130 return StatusCode::SUCCESS;
131 }
132
133
137 delete thinnedMcEvts;
138 thinnedMcEvts = nullptr;
139 return StatusCode::FAILURE;
140 }
142
143
144 const HepMC::GenEvent* mcEvt = mcEvts->
front();
145 const auto &wtCont = mcEvt->weights();
146 if (!wtCont.empty()) {
147 } else {
149 }
150 int inEvent = mcEvt->event_number();
152
154
156
157 HepMC::GenEvent* thinEvt = new HepMC::GenEvent(*mcEvt);
158 int nEvent = thinEvt->event_number();
160
161 if (doPrint) {
162 std::cout << "========== BEGIN EVENT BEFORE THINNING ==========" << std::endl;
164 std::cout << "========== END EVENT BEFORE THINNING ==========" << std::endl;
165 }
166
168
170
171
172
173
174
175
176
177
178
179
180
181
182
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
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
201
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
212
213
214 bool isHadVtx = true;
215 bool isHadOut = false;
216 for (const auto& inp: hadv->particles_in() ) {
218 }
219 for (const auto& vp: hadv->particles_out()) {
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
238
239
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) {
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) {
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()) {
263 return StatusCode::SUCCESS;
264 }
265
267
268
270
271#ifdef HEPMC3
272 for (unsigned int iv = 0; iv < hadVertices.size(); ++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
280
281 for (unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
283 auto descendants = HepMC::descendant_particles(ivtx);
284 for (auto pout: descendants) {
286 if (vpar) removePV.push_back(vpPair(vpar, pout));
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
304
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);
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
320
321
322
323
324#ifdef HEPMC3
325 for (auto hadv:thinEvt->vertices()) {
326
327 if (hadv->particles_in().size() == 0 && hadv->particles_out().size() == 0) {
328 removeV.push_back(hadv);
329 deleteV.push_back(hadv);
330 }
331
333 for (auto pin: hadv->particles_in()) {
334 removePV.push_back(vpPair(hadv, pin));
336 }
337 for (auto pout: hadv->particles_out()) {
338 removePV.push_back(vpPair(hadv, 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
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
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);
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();
364 removePV.emplace_back(*hadv, *pout);
366 }
367 removeV.push_back(*hadv);
368 deleteV.push_back(*hadv);
369 }
370#endif
371
372
373
374 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
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) {
388 v->add_particle_out(std::move(p));
389 }
390#ifdef HEPMC3
391 for (unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
393 if (
v->particles_in().size() != 0 ||
v->particles_out().size() != 0) {
395 <<
" with in/out particles " <<
v->particles_in().size() <<
" " <<
v->particles_out().size());
396 }
397 thinEvt->remove_vertex(hadVertices[iv]);
398 }
399
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 }
408 }
409
410
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) {
417 if (*dvItr) delete (*dvItr);
418 }
419
420 deleteP.sort();
421 deleteP.unique();
422 for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
424 if (*dpItr) delete (*dpItr);
425 }
426#endif
427
429
431
432 if (doDebug && doExtra) {
433 std::cout << "========== BEGIN EVENT BEFORE CLUSTER ==========" << std::endl;
435 std::cout << "========== END EVENT BEFORE CLUSTER ==========" << std::endl;
436 }
437
438
439
440
441
442
443
444
445
446
447
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
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
473 for (auto fp: thinEvt->particles() ) {
474 int iCase = 0;
477
479 if (!pvtx) {
481 continue;
482 }
483 if (doDebug)
ATH_MSG_DEBUG(
"Final parton " << pvtx <<
" " << fp);
485
487
488
489 if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 1) {
490
492 if (!pp || pp->parent_event() == nullptr) {
495 continue;
496 }
497
499 if (!ppvtx || ppvtx->parent_event() == nullptr) {
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
517
518
519
520
521 if (pvtx->particles_in().size() == 2 && pvtx->particles_out().size() == 1) {
522
525
526
527 if (std::abs(pp1->momentum().perp()) < 1.e-3) continue;
528 if (std::abs(pp2->momentum().perp()) < 1.e-3) continue;
529
532 if (!ppvtx1 || ppvtx1->parent_event() == nullptr) {
535 continue;
536 }
537 if (!ppvtx2 || ppvtx2->parent_event() == nullptr) {
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
558
559
560
561 if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 2) {
564
565 if (fp == pout1) {
568 continue;
569 }
570 } else if (fp == pout2) {
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
582
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
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
597 if (doDebug)
ATH_MSG_DEBUG(
"Keeping 1->2: parton mass " << m12);
598 continue;
599 }
600
602 if (!ppvtx || ppvtx->parent_event() == nullptr) {
605 continue;
606 }
607
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 }
626
628
629 if (pvtx->particles_in().size() == 1) {
630
632 if (std::abs(pp->pdg_id()) ==
MC::PROTON) iCase = -1;
633 }
634
635
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 }
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
647 for (finp = finpB; finp != finpE; ++finp) {
648 int iCase = 0;
649
650 HepMC::GenParticle*
fp = *finp;
653
654
655 HepMC::GenVertex* pvtx =
fp->production_vertex();
656 if (!pvtx) {
658 continue;
659 }
661
663
665
666
667
668
669 if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 1) {
670
671 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
672 HepMC::GenParticle* pp = *pitr;
676 continue;
677 }
678
679 HepMC::GenVertex* ppvtx = pp->production_vertex();
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);
695 }
696
698
700
701
702
703
704
705
706 if (pvtx->particles_in_size() == 2 && pvtx->particles_out_size() == 1) {
707
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
714
715 if (std::abs(pp1->momentum().perp()) < 1.e-3) continue;
716 if (std::abs(pp2->momentum().perp()) < 1.e-3) continue;
717
718 HepMC::GenVertex* ppvtx1 = pp1->production_vertex();
719 HepMC::GenVertex* ppvtx2 = pp2->production_vertex();
723 continue;
724 }
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) {
744 }
745 }
746
748
750
751
752
753
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
762 if (fp == pout1) {
765 continue;
766 }
767 } else if (fp == pout2) {
770 continue;
771 }
772 } else {
774 continue;
775 }
776 if (fp != pout1) continue;
777
778 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
779 HepMC::GenParticle* pp = *pitr;
780
781
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
787 double m12 = p12.m();
788 if (m12 < 0) {
789 if (fabs(m12) > 10. + 1.0e-5 * p12.e()) {
791 }
792 m12 = 0;
793 }
794 if (doDebug)
ATH_MSG_DEBUG(
"1->2: parton pair mass " << m12);
795
797 if (doDebug)
ATH_MSG_DEBUG(
"Keeping 1->2: parton mass " << m12);
798 continue;
799 }
800
801
802 HepMC::GenVertex* ppvtx = pp->production_vertex();
806 continue;
807 }
808
809
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 }
830
832
834
835
836 if (pvtx->particles_in_size() == 1) {
837
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
844
845 if (iCase == 0) {
847 }
848
849 }
850#endif
851
852
853
854 if (doDebug)
ATH_MSG_DEBUG(
"Actually removing particles " << removePV.size());
855
856 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
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
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) {
874 v->add_particle_out(std::move(p));
875 }
876
877
878 for (
unsigned int i = 0;
i < changePK.size(); ++
i) {
881 pp->set_momentum(changePK[i].second);
882 }
883
884
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 {
897 }
898#endif
899 }
900 if (doDebug)
ATH_MSG_DEBUG(
"End while(moreP) pass " << moreP);
901
902 }
903
904#ifdef HEPMC3
905
906#else
907
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) {
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) {
921 if (*dpItr) delete (*dpItr);
922 }
923
924#endif
926
928
929 if (doDebug && doExtra) {
930 std::cout << "========== BEGIN EVENT BEFORE SOFT ==========" << std::endl;
932 std::cout << "========== END EVENT BEFORE SOFT ==========" << std::endl;
933 }
934
935
936
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()) {
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
957 pNotHad.push_back(fp);
958 int ida = std::abs(
fp->pdg_id());
959 bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
961 pHard.push_back(fp);
962 pHard.insert(pHard.end(),ancestors.begin(),ancestors.end());
963 }
964
965
966
967
968
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
977 pNotHad.sort();
978 pNotHad.unique();
979 pHard.sort();
980 pHard.unique();
981
982
983 if (doDebug) {
984 int nhard = 0;
985 for (auto ph: pHard) {
986 ++nhard;
988 }
989 if (doDebug)
ATH_MSG_DEBUG(
"Hard GenParticles total " << nhard);
990 }
991
992#else
993
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
1020 pNotHad.push_back(fp);
1021 int ida = std::abs(
fp->pdg_id());
1022 bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
1024 pHard.push_back(fp);
1025 for (gpar = gparB; gpar != gparE; ++gpar)
1026 pHard.push_back(*gpar);
1027 }
1028
1029
1030
1031
1032
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
1043 pNotHad.sort();
1044 pNotHad.unique();
1045 pHard.sort();
1046 pHard.unique();
1047
1048
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;
1056 }
1057 if (doDebug)
ATH_MSG_DEBUG(
"Hard GenParticles total " << nhard);
1058 }
1059
1060#endif
1061
1062
1063
1064#ifdef HEPMC3
1065 for (auto p: pNotHad) {
1066
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;
1077 if (pvtx) pvtx->remove_particle_out(p);
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
1094 bool isHard = false;
1095 for (hItr = hItrB; hItr != hItrE; ++hItr) {
1096 if (p == (*hItr)) {
1097 isHard = true;
1098 break;
1099 }
1100 }
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);
1108 }
1109#endif
1110
1112
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) {
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])) {
1144 } else {
1146 }
1147 }
1148
1149 deleteV.sort();
1150 deleteV.unique();
1151 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
1153 if (*dvItr) delete (*dvItr);
1154 }
1155#endif
1156
1158
1160
1161 if (doDebug && doExtra) {
1162 std::cout << "========== BEGIN EVENT BEFORE 1-BODY ==========" << std::endl;
1164 std::cout << "========== END EVENT BEFORE 1-BODY ==========" << std::endl;
1165 }
1166
1167
1168
1169
1170
1171#ifdef HEPMC3
1172 bool moreV1 = true;
1176
1177 while (moreV1) {
1178 moreV1 = false;
1179
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());
1185 if (pin->pdg_id() !=
pout->pdg_id())
continue;
1186
1187 if (pin == beams[0] || pin == beams[1]) continue;
1189 if (!pvtx || pvtx->parent_event() == nullptr) {
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) {
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
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
1232 if (pin == beams[0] || pin == beams[1]) continue;
1233 HepMC::GenVertex* pvtx = pin->production_vertex();
1237 continue;
1238 }
1239
1240 moreV1 = true;
1241 vtx11 = (*v);
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
1262
1263
1264
1265
1266
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();
1280 if (pp->production_vertex()) continue;
1281 double pt = pp->momentum().perp();
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));
1291 }
1292
1293 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
1296 v->remove_particle_in(p);
1297 v->remove_particle_out(std::move(p));
1298 }
1299
1300 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1301 thinEvt->remove_vertex(removeV[i]);
1302 }
1303
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();
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);
1324 }
1325
1326
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
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
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 }
1352
1354
1356
1357 if (doPrint) {
1358 std::cout << "========== BEGIN EVENT AFTER THINNING ==========" << std::endl;
1360 std::cout << "========== END EVENT AFTER THINNING ==========" << std::endl;
1361 }
1362
1365
1367
1369
1371 return StatusCode::SUCCESS;
1372}
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
ServiceHandle< StoreGateSvc > & evtStore()
value_type push_back(value_type pElem)
Add an element to the end of the collection.
const T * front() const
Access the first element in the collection as an rvalue.
bool empty() const noexcept
Returns true if the collection is empty.
virtual HepMC::FourVector vtxOutMom(HepMC::ConstGenVertexPtr)
void line(std::ostream &os, const GenEvent &e)
bool is_simulation_vertex(const T &v)
Method to establish if the vertex was created during simulation (TODO migrate to be based on status).
HepMC::GenVertex * GenVertexPtr
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...
GenParticle * GenParticlePtr
bool isParton(const T &p)
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...
bool isHadron(const T &p)
pout(output, newline=True)
retrieve(aClass, aKey=None)