42 {
43
44
47 if (!inputMcEventCollection.isValid()) {
48 ATH_MSG_ERROR(
"Could not find input McEventCollection called " << inputMcEventCollection.name() <<
" in store " << inputMcEventCollection.store() <<
".");
49 return StatusCode::FAILURE;
50 }
51 const HepMC::GenEvent& inputEvent(**(inputMcEventCollection->begin()));
52
53
54 bool inputProblem(false);
55 for (const auto& particle: inputEvent) {
57 if (!
particle->production_vertex()) {
58 ATH_MSG_ERROR(
"Stable particle without a production vertex!! " << particle);
59 inputProblem = true;
60 }
61 }
63 if (!
particle->production_vertex()) {
64 ATH_MSG_ERROR(
"Decyed particle without a production vertex!! " << particle);
65 inputProblem = true;
66 }
68 ATH_MSG_ERROR(
"Decyed particle without an end vertex!! " << particle);
69 inputProblem = true;
70 }
71 }
72 }
73 if (inputProblem) {
75 return StatusCode::FAILURE;
76 }
77#ifdef HEPMC3
80 std::unique_ptr<HepMC::GenEvent> outputEvent = std::make_unique<HepMC::GenEvent>(inputEvent);
81 if (inputEvent.run_info()) outputEvent->set_run_info(std::make_shared<HepMC3::GenRunInfo>(*(inputEvent.run_info().get())));
82 for (;;) {
83 std::vector<HepMC::GenParticlePtr> p_to_remove;
84 std::vector<HepMC::GenVertexPtr> v_to_remove;
85 for (auto& particle: outputEvent->particles()) {
87 p_to_remove.push_back(particle);
88 }
89 }
90 for (auto& particle: p_to_remove) outputEvent->remove_particle(particle);
91 for (auto& vertex: outputEvent->vertices()) {
93 v_to_remove.push_back(vertex);
94 }
95 }
96 for (auto& vertex: v_to_remove) outputEvent->remove_vertex(vertex);
97 if (p_to_remove.empty() && v_to_remove.empty()) break;
98 }
99#else
100
101 std::unique_ptr<HepMC::GenEvent> outputEvent = std::make_unique<HepMC::GenEvent>(inputEvent.signal_process_id(),
102 inputEvent.event_number());
103
104 outputEvent->set_mpi ( inputEvent.mpi() );
105 outputEvent->set_event_scale ( inputEvent.event_scale() );
106 outputEvent->set_alphaQCD ( inputEvent.alphaQCD() );
107 outputEvent->set_alphaQED ( inputEvent.alphaQED() );
108 if ( inputEvent.cross_section() ) {
109 outputEvent->set_cross_section ( *inputEvent.cross_section());
110 }
111 if (inputEvent.heavy_ion()) {
112 outputEvent->set_heavy_ion(*(inputEvent.heavy_ion()));
113 }
114 if (inputEvent.pdf_info()) {
115 outputEvent->set_pdf_info(*(inputEvent.pdf_info()));
116 }
117 outputEvent->define_units( inputEvent.momentum_unit(), inputEvent.length_unit() );
118
119
120
121
122
123
124 std::map<const HepMC::GenVertex*,HepMC::GenVertex*> inputEvtVtxToOutputEvtVtx;
125 HepMC::GenEvent::vertex_const_iterator currentVertexIter(inputEvent.vertices_begin());
126 const HepMC::GenEvent::vertex_const_iterator endOfCurrentListOfVertices(inputEvent.vertices_end());
128 for (; currentVertexIter != endOfCurrentListOfVertices; ++currentVertexIter) {
129 const HepMC::GenVertex *pCurrentVertex(*currentVertexIter);
131 continue;
132 }
133 std::unique_ptr<HepMC::GenVertex> copyOfGenVertex =std::make_unique<HepMC::GenVertex>(pCurrentVertex->position(), pCurrentVertex->id(), pCurrentVertex->weights() );
134 copyOfGenVertex->suggest_barcode(
HepMC::barcode(pCurrentVertex) );
135 inputEvtVtxToOutputEvtVtx[pCurrentVertex] = copyOfGenVertex.get();
136 outputEvent->add_vertex( copyOfGenVertex.release() );
137 }
138
139
140 if ( inputEvent.signal_process_vertex() ) {
141 outputEvent->set_signal_process_vertex( inputEvtVtxToOutputEvtVtx[inputEvent.signal_process_vertex()] );
142 }
143 else {
144 outputEvent->set_signal_process_vertex( nullptr );
145 }
146
147
148
149 HepMC::GenParticle*
beam1{};
150 HepMC::GenParticle*
beam2{};
151 for ( HepMC::GenEvent::particle_const_iterator particleIter = inputEvent.particles_begin();
152 particleIter != inputEvent.particles_end(); ++particleIter ) {
153 const HepMC::GenParticle* currentParticle = *particleIter;
155 continue;
156 }
157 std::unique_ptr<HepMC::GenParticle> copyOfGenParticle = std::make_unique<HepMC::GenParticle>(*currentParticle);
158 const bool isBeamParticle1(currentParticle == inputEvent.beam_particles().first);
159 const bool isBeamParticle2(currentParticle == inputEvent.beam_particles().second);
160
161
162
163 const bool shouldAddProdVertex(currentParticle->production_vertex() && inputEvtVtxToOutputEvtVtx[ currentParticle->production_vertex() ]);
164 const bool shouldAddEndVertex(currentParticle->end_vertex() && inputEvtVtxToOutputEvtVtx[ currentParticle->end_vertex() ]);
165 if ( isBeamParticle1 || isBeamParticle2 || shouldAddProdVertex || shouldAddEndVertex ) {
166 HepMC::GenParticle* particleCopy = copyOfGenParticle.release();
167 if ( isBeamParticle1 ) {
168 beam1 = particleCopy;
169 }
170 if ( isBeamParticle2 ) {
171 beam2 = particleCopy;
172 }
173 if ( shouldAddProdVertex || shouldAddEndVertex ) {
174 if ( shouldAddEndVertex ) {
175 inputEvtVtxToOutputEvtVtx[ currentParticle->end_vertex() ]->
176 add_particle_in(particleCopy);
177 }
178 if ( shouldAddProdVertex ) {
179 inputEvtVtxToOutputEvtVtx[ currentParticle->production_vertex() ]->
180 add_particle_out(particleCopy);
181 }
182 }
183 else {
184 ATH_MSG_WARNING (
"Found GenParticle with no production or end vertex! \n" << *currentParticle);
185 }
186 }
187 }
188 outputEvent->set_beam_particles( beam1, beam2 );
189
190
191 outputEvent->set_random_states( inputEvent.random_states() );
192 outputEvent->weights() = inputEvent.weights();
193#endif
194
195
196 bool outputProblem(false);
197 for (const auto& particle: *(outputEvent.get())) {
199 if (!
particle->production_vertex()) {
200 ATH_MSG_ERROR(
"Stable particle without a production vertex!! " << particle);
201 outputProblem = true;
202 }
204 ATH_MSG_ERROR(
"Stable particle with an end vertex!! " << particle);
205 outputProblem = true;
206 }
207 }
209 if (!
particle->production_vertex()) {
210 ATH_MSG_ERROR(
"Decayed particle without a production vertex!! " << particle);
211 outputProblem = true;
212 }
214 ATH_MSG_ERROR(
"Decayed particle without an end vertex!! " << particle);
215 outputProblem = true;
216 }
217 }
218 }
219 if (outputProblem) {
221 return StatusCode::FAILURE;
222 }
223
225 ATH_CHECK(outputMcEventCollection.record(std::make_unique<McEventCollection>()));
226 outputMcEventCollection->push_back(outputEvent.release());
227 if (!outputMcEventCollection.isValid()) {
228 ATH_MSG_ERROR(
"Could not record output McEventCollection called " << outputMcEventCollection.name() <<
" in store " << outputMcEventCollection.store() <<
".");
229 return StatusCode::FAILURE;
230 }
231
233
234 return StatusCode::SUCCESS;
235
236}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
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