ATLAS Offline Software
initMcEventCollection.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 #include "GeneratorObjectsTPCnv/initMcEventCollection.h"
5 
6 // HepMC includes
7 #include "AtlasHepMC/GenEvent.h"
8 #include "AtlasHepMC/GenParticle.h"
9 #include "AtlasHepMC/GenVertex.h"
10 
11 // CLHEP includes
12 #include "CLHEP/Units/SystemOfUnits.h"
13 
14 #include "GeneratorObjectsTPCnv/HepMcParticleLinkCnv_p1.h"
15 #include "StoreGate/WriteHandle.h"
16 #include "GeneratorObjects/McEventCollection.h"
17 #include "TruthUtils/MagicNumbers.h"
18 #include "GaudiKernel/ThreadLocalContext.h"
19 
20 #include "TestTools/initGaudi.h"
21 
22 namespace Athena_test {
23  bool initMcEventCollection(ISvcLocator*& pSvcLoc, std::vector<HepMC::GenParticlePtr>& genPartList, const bool initGaudi)
24  {
25  if (initGaudi && !Athena_test::initGaudi(pSvcLoc)) {
26  std::cerr << "This test can not be run" << std::endl;
27  return false;
28  }
29  // create dummy input McEventCollection with a name that
30  // HepMcParticleLink knows about
31  SG::WriteHandle<McEventCollection> inputTestDataHandle{"TruthEvent"};
32  inputTestDataHandle = std::make_unique<McEventCollection>();
33 
34  // create a dummy EventContext
35  EventContext ctx;
36  ctx.setExtension( Atlas::ExtendedEventContext( SG::CurrentEventStore::store() ) );
37  Gaudi::Hive::setCurrentContext( ctx );
38 
39  // Add a dummy GenEvent
40  const int process_id1(20);
41  const int event_number1(17);
42  inputTestDataHandle->push_back(HepMC::newGenEvent(process_id1, event_number1));
43  HepMC::GenEvent& ge1 = *(inputTestDataHandle->at(0));
44  populateGenEvent(ge1,-11,11,genPartList);
45  populateGenEvent(ge1,-13,13,genPartList);
46  populateGenEvent(ge1,-11,11,genPartList);
47  populateGenEvent(ge1,-13,13,genPartList);
48  populateGenEvent(ge1,-11,11,genPartList);
49  populateGenEvent(ge1,22,22,genPartList);
50  inputTestDataHandle->push_back(new HepMC::GenEvent(ge1));
51  HepMC::GenEvent& ge2 = *(inputTestDataHandle->at(1));
52  const int event_number2(89);
53  ge2.set_event_number(event_number2);
54  populateFilteredGenEvent(ge2,genPartList);
55  return true;
56  }
57 
58  int maximumBarcode(std::vector<HepMC::GenParticlePtr>& genPartList)
59  {
60  int maxBarcode{0};
61  if (genPartList.empty()) { return maxBarcode; }
62  for (HepMC::GenParticlePtr genPart : genPartList)
63  {
64  maxBarcode = std::max(maxBarcode, HepMC::barcode(genPart));
65  }
66  return maxBarcode;
67  }
68 
69  void populateGenEvent(HepMC::GenEvent & ge, int pdgid1, int pdgid2, std::vector<HepMC::GenParticlePtr>& genPartList)
70  {
71  int maxBarcode = maximumBarcode(genPartList);
72  HepMC::FourVector myPos( 0.0, 0.0, 0.0, 0.0);
73  HepMC::GenVertexPtr myVertex = HepMC::newGenVertexPtr( myPos, -1 );
74  HepMC::FourVector fourMomentum1( 0.0, 0.0, 1.0, 1.0*CLHEP::TeV);
75  HepMC::GenParticlePtr inParticle1 = HepMC::newGenParticlePtr(fourMomentum1, pdgid1, 2);
76  myVertex->add_particle_in(inParticle1);
77  HepMC::FourVector fourMomentum2( 0.0, 0.0, -1.0, 1.0*CLHEP::TeV);
78  HepMC::GenParticlePtr inParticle2 = HepMC::newGenParticlePtr(fourMomentum2, pdgid2, 2);
79  myVertex->add_particle_in(inParticle2);
80  HepMC::FourVector fourMomentum3( 0.0, 1.0, 0.0, 1.0*CLHEP::TeV);
81  HepMC::GenParticlePtr inParticle3 = HepMC::newGenParticlePtr(fourMomentum3, pdgid1, 1);
82  myVertex->add_particle_out(inParticle3);
83  genPartList.push_back(inParticle3);
84  HepMC::FourVector fourMomentum4( 0.0, -1.0, 0.0, 1.0*CLHEP::TeV);
85  HepMC::GenParticlePtr inParticle4 = HepMC::newGenParticlePtr(fourMomentum4, pdgid2, 1);
86  myVertex->add_particle_out(inParticle4);
87  genPartList.push_back(inParticle4);
88  ge.add_vertex( myVertex );
89  HepMC::suggest_barcode(inParticle1,maxBarcode+1);
90  HepMC::suggest_barcode(inParticle2,maxBarcode+2);
91  HepMC::suggest_barcode(inParticle3,maxBarcode+3);
92  HepMC::suggest_barcode(inParticle4,maxBarcode+4);
93  HepMC::set_signal_process_vertex(&ge, myVertex );
94  ge.set_beam_particles(inParticle1,inParticle2);
95  }
96 
97  void populateFilteredGenEvent(HepMC::GenEvent & ge, std::vector<HepMC::GenParticlePtr>& genPartVector)
98  {
99  //.......Create new particle (geantino) to link hits from pileup
100  HepMC::GenParticlePtr genPart=HepMC::newGenParticlePtr();
101  genPart->set_pdg_id(999); //Geantino
102  genPart->set_status(1);
103  HepMC::suggest_barcode(genPart, HepMC::SUPPRESSED_PILEUP_BARCODE );
104 
105  HepMC::GenVertexPtr genVertex=HepMC::newGenVertexPtr();
106  genVertex->add_particle_out(genPart);
107  genPartVector.push_back(genPart);
108 
109  //to set geantino vertex as a truth primary vertex
110  HepMC::GenVertexPtr hScatVx = HepMC::barcode_to_vertex(&ge,-3);
111  if (hScatVx!=nullptr) {
112  HepMC::FourVector pmvxpos=hScatVx->position();
113  genVertex->set_position(pmvxpos);
114  //to set geantino kinematic phi=eta=0, E=p=E_hard_scat
115 #ifdef HEPMC3
116  auto itrp =hScatVx->particles_in().cbegin();
117  if (hScatVx->particles_in().size()==2) {
118  HepMC::FourVector mom1=(*itrp)->momentum();
119  HepMC::FourVector mom2=(*(++itrp))->momentum();
120  HepMC::FourVector vxmom;
121  vxmom.setPx(mom1.e()+mom2.e());
122  vxmom.setPy(0.);
123  vxmom.setPz(0.);
124  vxmom.setE(mom1.e()+mom2.e());
125  genPart->set_momentum(vxmom);
126  }
127 #else
128  HepMC::GenVertex::particles_in_const_iterator itrp =hScatVx->particles_in_const_begin();
129  if (hScatVx->particles_in_size()==2) {
130  HepMC::FourVector mom1=(*itrp)->momentum();
131  HepMC::FourVector mom2=(*(++itrp))->momentum();
132  HepMC::FourVector vxmom;
133  vxmom.setPx(mom1.e()+mom2.e());
134  vxmom.setPy(0.);
135  vxmom.setPz(0.);
136  vxmom.setE(mom1.e()+mom2.e());
137  genPart->set_momentum(vxmom);
138  }
139 #endif
140  }
141 
142 #ifdef HEPMC3
143  if (!ge.vertices().empty()) {
144  std::vector<HepMC::GenVertexPtr> vtxvec;
145  for (const auto& vtx: ge.vertices()) {
146  vtxvec.push_back(vtx);
147  ge.remove_vertex(vtx);
148  }
149  vtxvec.clear();
150  }
151 #else
152  if (!ge.vertices_empty()) {
153  HepMC::GenEvent::vertex_iterator itvtx = ge.vertices_begin();
154  while (itvtx != ge.vertices_end()) {
155  HepMC::GenVertexPtr pvtx = *itvtx;
156  ++itvtx;
157  ge.remove_vertex(pvtx);
158  delete pvtx;
159  }
160  }
161 #endif
162 
163  //.....add new vertex with geantino
164  ge.add_vertex(genVertex);
165  HepMC::suggest_barcode(genPart, HepMC::SUPPRESSED_PILEUP_BARCODE );
166  }
167 }
HepMC::GenVertexPtr
HepMC::GenVertex * GenVertexPtr
Definition: GenVertex.h:59
HepMC::suggest_barcode
bool suggest_barcode(T &p, int i)
Definition: GenEvent.h:670
GenEvent.h
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
GenVertex.h
HepMC::GenParticlePtr
GenParticle * GenParticlePtr
Definition: GenParticle.h:37
HepMC::newGenEvent
GenEvent * newGenEvent(const int a, const int b)
Definition: GenEvent.h:624
SG::CurrentEventStore::store
static IProxyDict * store()
Fetch the current store.
python.SystemOfUnits.TeV
int TeV
Definition: SystemOfUnits.py:158
GenParticle.h
Athena_test::populateFilteredGenEvent
void populateFilteredGenEvent(HepMC::GenEvent &ge, std::vector< HepMC::GenParticlePtr > &genPartVector)
Definition: initMcEventCollection.cxx:97
WriteHandle.h
Handle class for recording to StoreGate.
HepMC::SUPPRESSED_PILEUP_BARCODE
constexpr int SUPPRESSED_PILEUP_BARCODE(std::numeric_limits< int32_t >::max())
This barcode is used by objects matched to particles from pile-up interactions in standard MC Product...
HepMC::set_signal_process_vertex
void set_signal_process_vertex(GenEvent *e, T v)
Definition: GenEvent.h:650
McEventCollection.h
HepMC::newGenVertexPtr
GenVertexPtr newGenVertexPtr(const HepMC::FourVector &pos=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), const int i=0)
Definition: GenVertex.h:64
Atlas::ExtendedEventContext
Definition: ExtendedEventContext.h:23
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
Athena_test
functions & macros to test the difference between floats
Definition: InitGaudiGoogleTest.h:31
initMcEventCollection.h
minimal gaudi initialization and record an McEventCollection in StoreGate
initGaudi.h
minimal gaudi initialization for AthenaServices unit testing
Athena_test::initGaudi
bool initGaudi(ISvcLocator *&pSvcLoc)
Minimal Gaudi initialization for unit testing without job options.
Definition: initGaudi.cxx:28
Athena_test::populateGenEvent
void populateGenEvent(HepMC::GenEvent &ge, int pdgid1, int pdgid2, std::vector< HepMC::GenParticlePtr > &genPartVector)
Definition: initMcEventCollection.cxx:69
Athena_test::initMcEventCollection
bool initMcEventCollection(ISvcLocator *&pSvcLoc, std::vector< HepMC::GenParticlePtr > &genPartVector, const bool initGaudi=true)
Definition: initMcEventCollection.cxx:23
HepMcParticleLinkCnv_p1.h
MagicNumbers.h
SG::WriteHandle< McEventCollection >
HepMC::newGenParticlePtr
GenParticlePtr newGenParticlePtr(const HepMC::FourVector &mom=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), int pid=0, int status=0)
Definition: GenParticle.h:39
HepMC::barcode_to_vertex
GenVertex * barcode_to_vertex(const GenEvent *e, int id)
Definition: GenEvent.h:627
Athena_test::maximumBarcode
int maximumBarcode(std::vector< HepMC::GenParticlePtr > &genPartList)
Definition: initMcEventCollection.cxx:58
Generated on Sun Dec 22 2024 21:11:50 for ATLAS Offline Software by   doxygen 1.8.18