34 {
35
36
39 if (!inputMcEventCollection.isValid()) {
40 ATH_MSG_ERROR(
"Could not find input McEventCollection called " << inputMcEventCollection.name() <<
" in store " << inputMcEventCollection.store() <<
".");
41 return StatusCode::FAILURE;
42 }
43 const HepMC::GenEvent& inputEvent(**(inputMcEventCollection->begin()));
44
45
46 bool inputProblem(false);
47 for (const auto& particle: inputEvent) {
49 if (!
particle->production_vertex()) {
50 ATH_MSG_ERROR(
"Stable particle without a production vertex!! " << particle);
51 inputProblem = true;
52 }
53 }
55 if (!
particle->production_vertex()) {
56 ATH_MSG_ERROR(
"Decyed particle without a production vertex!! " << particle);
57 inputProblem = true;
58 }
60 ATH_MSG_ERROR(
"Decyed particle without an end vertex!! " << particle);
61 inputProblem = true;
62 }
63 }
64 }
65 if (inputProblem) {
67 return StatusCode::FAILURE;
68 }
69#ifdef HEPMC3
72 std::unique_ptr<HepMC::GenEvent> outputEvent = std::make_unique<HepMC::GenEvent>(inputEvent);
73 if (inputEvent.run_info()) outputEvent->set_run_info(std::make_shared<HepMC3::GenRunInfo>(*(inputEvent.run_info().get())));
74 for (;;) {
75 std::vector<HepMC::GenParticlePtr> p_to_remove;
76 std::vector<HepMC::GenVertexPtr> v_to_remove;
77 for (auto& particle: outputEvent->particles()) {
79 p_to_remove.push_back(particle);
80 }
81 }
82 for (auto& particle: p_to_remove) outputEvent->remove_particle(particle);
83 for (auto& vertex: outputEvent->vertices()) {
85 v_to_remove.push_back(vertex);
86 }
87 }
88 for (auto& vertex: v_to_remove) outputEvent->remove_vertex(vertex);
89 if (p_to_remove.empty() && v_to_remove.empty()) break;
90 }
91#else
92
93 std::unique_ptr<HepMC::GenEvent> outputEvent = std::make_unique<HepMC::GenEvent>(inputEvent.signal_process_id(),
94 inputEvent.event_number());
95
96 outputEvent->set_mpi ( inputEvent.mpi() );
97 outputEvent->set_event_scale ( inputEvent.event_scale() );
98 outputEvent->set_alphaQCD ( inputEvent.alphaQCD() );
99 outputEvent->set_alphaQED ( inputEvent.alphaQED() );
100 if ( inputEvent.cross_section() ) {
101 outputEvent->set_cross_section ( *inputEvent.cross_section());
102 }
103 if (inputEvent.heavy_ion()) {
104 outputEvent->set_heavy_ion(*(inputEvent.heavy_ion()));
105 }
106 if (inputEvent.pdf_info()) {
107 outputEvent->set_pdf_info(*(inputEvent.pdf_info()));
108 }
109 outputEvent->define_units( inputEvent.momentum_unit(), inputEvent.length_unit() );
110
111
112
113
114
115
116 std::map<const HepMC::GenVertex*,HepMC::GenVertex*> inputEvtVtxToOutputEvtVtx;
117 HepMC::GenEvent::vertex_const_iterator currentVertexIter(inputEvent.vertices_begin());
118 const HepMC::GenEvent::vertex_const_iterator endOfCurrentListOfVertices(inputEvent.vertices_end());
120 for (; currentVertexIter != endOfCurrentListOfVertices; ++currentVertexIter) {
121 const HepMC::GenVertex *pCurrentVertex(*currentVertexIter);
123 continue;
124 }
125 std::unique_ptr<HepMC::GenVertex> copyOfGenVertex =std::make_unique<HepMC::GenVertex>(pCurrentVertex->position(), pCurrentVertex->id(), pCurrentVertex->weights() );
126 copyOfGenVertex->suggest_barcode(
HepMC::barcode(pCurrentVertex) );
127 inputEvtVtxToOutputEvtVtx[pCurrentVertex] = copyOfGenVertex.get();
128 outputEvent->add_vertex( copyOfGenVertex.release() );
129 }
130
131
132 if ( inputEvent.signal_process_vertex() ) {
133 outputEvent->set_signal_process_vertex( inputEvtVtxToOutputEvtVtx[inputEvent.signal_process_vertex()] );
134 }
135 else {
136 outputEvent->set_signal_process_vertex( nullptr );
137 }
138
139
140
141 HepMC::GenParticle*
beam1{};
142 HepMC::GenParticle*
beam2{};
143 for ( HepMC::GenEvent::particle_const_iterator particleIter = inputEvent.particles_begin();
144 particleIter != inputEvent.particles_end(); ++particleIter ) {
145 const HepMC::GenParticle* currentParticle = *particleIter;
147 continue;
148 }
149 std::unique_ptr<HepMC::GenParticle> copyOfGenParticle = std::make_unique<HepMC::GenParticle>(*currentParticle);
150 const bool isBeamParticle1(currentParticle == inputEvent.beam_particles().first);
151 const bool isBeamParticle2(currentParticle == inputEvent.beam_particles().second);
152
153
154
155 const bool shouldAddProdVertex(currentParticle->production_vertex() && inputEvtVtxToOutputEvtVtx[ currentParticle->production_vertex() ]);
156 const bool shouldAddEndVertex(currentParticle->end_vertex() && inputEvtVtxToOutputEvtVtx[ currentParticle->end_vertex() ]);
157 if ( isBeamParticle1 || isBeamParticle2 || shouldAddProdVertex || shouldAddEndVertex ) {
158 HepMC::GenParticle* particleCopy = copyOfGenParticle.release();
159 if ( isBeamParticle1 ) {
160 beam1 = particleCopy;
161 }
162 if ( isBeamParticle2 ) {
163 beam2 = particleCopy;
164 }
165 if ( shouldAddProdVertex || shouldAddEndVertex ) {
166 if ( shouldAddEndVertex ) {
167 inputEvtVtxToOutputEvtVtx[ currentParticle->end_vertex() ]->
168 add_particle_in(particleCopy);
169 }
170 if ( shouldAddProdVertex ) {
171 inputEvtVtxToOutputEvtVtx[ currentParticle->production_vertex() ]->
172 add_particle_out(particleCopy);
173 }
174 }
175 else {
176 ATH_MSG_WARNING (
"Found GenParticle with no production or end vertex! \n" << *currentParticle);
177 }
178 }
179 }
180 outputEvent->set_beam_particles( beam1, beam2 );
181
182
183 outputEvent->set_random_states( inputEvent.random_states() );
184 outputEvent->weights() = inputEvent.weights();
185#endif
186
187
188 bool outputProblem(false);
189 for (const auto& particle: *(outputEvent.get())) {
191 if (!
particle->production_vertex()) {
192 ATH_MSG_ERROR(
"Stable particle without a production vertex!! " << particle);
193 outputProblem = true;
194 }
196 ATH_MSG_ERROR(
"Stable particle with an end vertex!! " << particle);
197 outputProblem = true;
198 }
199 }
201 if (!
particle->production_vertex()) {
202 ATH_MSG_ERROR(
"Decayed particle without a production vertex!! " << particle);
203 outputProblem = true;
204 }
206 ATH_MSG_ERROR(
"Decayed particle without an end vertex!! " << particle);
207 outputProblem = true;
208 }
209 }
210 }
211 if (outputProblem) {
213 return StatusCode::FAILURE;
214 }
215
217 ATH_CHECK(outputMcEventCollection.record(std::make_unique<McEventCollection>()));
218 outputMcEventCollection->push_back(outputEvent.release());
219 if (!outputMcEventCollection.isValid()) {
220 ATH_MSG_ERROR(
"Could not record output McEventCollection called " << outputMcEventCollection.name() <<
" in store " << outputMcEventCollection.store() <<
".");
221 return StatusCode::FAILURE;
222 }
223
225
226 return StatusCode::SUCCESS;
227
228}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
SG::WriteHandleKey< McEventCollection > m_outputMcEventCollection
SG::ReadHandleKey< McEventCollection > m_inputMcEventCollection
bool is_simulation_vertex(const T &v)
Method to establish if the vertex was created during simulation (TODO migrate to be based on status).
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...
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
bool isDecayed(const T &p)
Identify if the particle decayed.
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
1.14.0