ATLAS Offline Software
TruthClosureCheck.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "TruthClosureCheck.h"
6 //
7 #include "AtlasHepMC/GenEvent.h"
8 #include "AtlasHepMC/GenVertex.h"
9 #include "GeneratorObjects/McEventCollection.h"
10 #include "TruthUtils/HepMCHelpers.h"
11 
12 
13 #include "StoreGate/ReadHandle.h"
14 
15 #include "CLHEP/Vector/LorentzVector.h"
16 #include "CLHEP/Units/SystemOfUnits.h"
17 //
18 #include "CLHEP/Geometry/Point3D.h"
19 #include "GeoPrimitives/GeoPrimitives.h"
20 #include <climits>
21 
22 TruthClosureCheck::TruthClosureCheck(const std::string& name, ISvcLocator* pSvcLocator):
23  AthAlgorithm(name, pSvcLocator)
24 {
25  declareProperty("OriginalMcEventCollection" , m_originalMcEventCollection = "BeamTruthEvent");
26  declareProperty("ResetMcEventCollection" , m_resetMcEventCollection = "NewTruthEvent");
27  declareProperty("CompareMomenta", m_compareMomenta);
28  declareProperty("PostSimulation", m_postSimulation);
29 }
30 
31 
32 //----------------------------------------------------
33 StatusCode TruthClosureCheck::initialize() {
34  //----------------------------------------------------
35 
36  ATH_CHECK(m_originalMcEventCollection.initialize());
37  ATH_CHECK(m_resetMcEventCollection.initialize());
38  return StatusCode::SUCCESS;
39 
40 }
41 
42 
43 StatusCode TruthClosureCheck::sanityCheck(const HepMC::GenEvent& event) const {
44  //Sanity check
45  bool resetProblem(false);
46  for (const auto& particle: event) {
47  if (MC::isStable(particle)) {
48  if (!particle->production_vertex()) {
49  ATH_MSG_ERROR("Stable particle without a production vertex!! " << particle);
50  resetProblem = true;
51  }
52  if (!m_postSimulation && particle->end_vertex()) {
53  ATH_MSG_ERROR("Stable particle with an end vertex!! " << particle);
54  resetProblem = true;
55  }
56  }
57  else if (MC::isDecayed(particle)) {
58  if (!particle->production_vertex()) {
59  ATH_MSG_ERROR("Decayed particle without a production vertex!! " << particle);
60  resetProblem = true;
61  }
62  if (!particle->end_vertex()) {
63  ATH_MSG_ERROR("Decayed particle without an end vertex!! " << particle);
64  resetProblem = true;
65  }
66  }
67  }
68  if (resetProblem) {
69  return StatusCode::FAILURE;
70  }
71 
72  return StatusCode::SUCCESS;
73 }
74 
75 #ifdef HEPMC3
76 void TruthClosureCheck::printGenVertex(const HepMC::ConstGenVertexPtr& origVertex,
77  const HepMC::ConstGenVertexPtr& resetVertex) const
78 {
79  ATH_MSG_INFO("----------------------------------");
80  ATH_MSG_INFO("Original Vertex:");
81  ATH_MSG_INFO( origVertex );
82  ATH_MSG_INFO("Particles In:");
83  for (const auto& originalPartIn: origVertex->particles_in()) ATH_MSG_INFO( originalPartIn );
84  ATH_MSG_INFO("Particles Out:");
85  for (const auto& originalPartOut: origVertex->particles_out()) ATH_MSG_INFO( originalPartOut );
86  ATH_MSG_INFO("----------------------------------");
87  ATH_MSG_INFO("Reset Vertex:");
88  ATH_MSG_INFO( resetVertex );
89  ATH_MSG_INFO("Particles In:");
90  for (const auto& resetPartIn: resetVertex->particles_in()) ATH_MSG_INFO( resetPartIn );
91  ATH_MSG_INFO("Particles Out:");
92  for (const auto& resetPartOut: resetVertex->particles_out()) ATH_MSG_INFO( resetPartOut );
93  ATH_MSG_INFO("----------------------------------");
94  return;
95 }
96 #else
97 
98 void TruthClosureCheck::printGenVertex(const HepMC::GenVertex& origVertex,
99  const HepMC::GenVertex& resetVertex) const
100 {
101  ATH_MSG_INFO("----------------------------------");
102  ATH_MSG_INFO("Original Vertex:");
103  ATH_MSG_INFO( origVertex );
104  ATH_MSG_INFO("Particles In:");
105  HepMC::GenVertex::particles_in_const_iterator originalPartInIter(origVertex.particles_in_const_begin());
106  const HepMC::GenVertex::particles_in_const_iterator endOfOriginalListOfParticlesIn(origVertex.particles_in_const_end());
107  while( originalPartInIter!=endOfOriginalListOfParticlesIn ) {
108  ATH_MSG_INFO( **originalPartInIter );
109  ++originalPartInIter;
110  }
111  ATH_MSG_INFO("Particles Out:");
112  HepMC::GenVertex::particles_out_const_iterator originalPartOutIter(origVertex.particles_out_const_begin());
113  const HepMC::GenVertex::particles_out_const_iterator endOfOriginalListOfParticlesOut(origVertex.particles_out_const_end());
114  while( originalPartOutIter!=endOfOriginalListOfParticlesOut) {
115  ATH_MSG_INFO( **originalPartOutIter );
116  ++originalPartOutIter;
117  }
118  ATH_MSG_INFO("----------------------------------");
119  ATH_MSG_INFO("Reset Vertex:");
120  ATH_MSG_INFO( resetVertex );
121  ATH_MSG_INFO("Particles In:");
122  HepMC::GenVertex::particles_in_const_iterator resetPartInIter(resetVertex.particles_in_const_begin());
123  const HepMC::GenVertex::particles_in_const_iterator endOfResetListOfParticlesIn(resetVertex.particles_in_const_end());
124  while( resetPartInIter!=endOfResetListOfParticlesIn ) {
125  ATH_MSG_INFO( **resetPartInIter );
126  ++resetPartInIter;
127  }
128  ATH_MSG_INFO("Particles Out:");
129  HepMC::GenVertex::particles_out_const_iterator resetPartOutIter(resetVertex.particles_out_const_begin());
130  const HepMC::GenVertex::particles_out_const_iterator endOfResetListOfParticlesOut(resetVertex.particles_out_const_end());
131  while( resetPartOutIter!=endOfResetListOfParticlesOut) {
132  ATH_MSG_INFO( **resetPartOutIter );
133  ++resetPartOutIter;
134  }
135  ATH_MSG_INFO("----------------------------------");
136  return;
137 }
138 #endif
139 
140 #ifdef HEPMC3
141 StatusCode TruthClosureCheck::compareGenVertex(const HepMC::ConstGenVertexPtr& origVertex,
142  const HepMC::ConstGenVertexPtr& resetVertex) const
143 {
144  if (!origVertex && !resetVertex) return StatusCode::SUCCESS;
145  if (!origVertex || !resetVertex) return StatusCode::FAILURE;
146  bool pass{true};
147 
148  if (HepMC::barcode(origVertex) != HepMC::barcode(resetVertex)) { // FIXME barcode-based
149  ATH_MSG_ERROR ("vertex barcode differs! Original: "<<HepMC::barcode(origVertex)<<", Reset: "<<HepMC::barcode(resetVertex));
150  pass = false;
151  }
152  if (origVertex->position() != resetVertex->position()) {
153  const HepMC::FourVector& oP = origVertex->position();
154  const HepMC::FourVector& rP = resetVertex->position();
155  ATH_MSG_ERROR("vertex position differs! Original: ("<<oP.x()<<","<<oP.y()<<","<<oP.z()<<","<<oP.t()<<"), Reset: ("<<rP.x()<<","<<rP.y()<<","<<rP.z()<<","<<rP.t()<<")");
156  pass = false;
157  }
158  if (origVertex->particles_in().size() != resetVertex->particles_in().size() ) {
159  ATH_MSG_ERROR ("particles_in_size differs! Original: "<<origVertex->particles_in().size()<<", Reset: "<<resetVertex->particles_in().size());
160  pass = false;
161  }
162  if (origVertex->particles_out().size() != resetVertex->particles_out().size() ) {
163  ATH_MSG_ERROR ("particles_out_size differs! Original: "<<origVertex->particles_out().size()<<", Reset: "<<resetVertex->particles_out().size());
164  pass = false;
165  }
166  std::vector<HepMC::ConstGenParticlePtr> OriginalListOfParticlesIn = origVertex->particles_in();
167  std::vector<HepMC::ConstGenParticlePtr> ResetListOfParticlesIn = resetVertex->particles_in();
168  for ( std::vector<HepMC::ConstGenParticlePtr>::iterator originalPartInIter(OriginalListOfParticlesIn.begin()),resetPartInIter(ResetListOfParticlesIn.begin());
169  originalPartInIter != OriginalListOfParticlesIn.end() && resetPartInIter != ResetListOfParticlesIn.end();
170  ++originalPartInIter, ++resetPartInIter
171  ) {
172  if (compareGenParticle(*originalPartInIter,*resetPartInIter).isFailure()) {
173  ATH_MSG_ERROR ( "input particle properties differ!" );
174  pass &= false;
175  }
176  ATH_MSG_VERBOSE ( "particles match!" );
177  }
178 
179 
180  std::vector<HepMC::ConstGenParticlePtr> OriginalListOfParticlesOut = origVertex->particles_out();
181  std::vector<HepMC::ConstGenParticlePtr> ResetListOfParticlesOut = resetVertex->particles_out();
182 //AV: please remember that the best quantities to compare particles are physical quantities.
183  std::sort(OriginalListOfParticlesOut.begin(), OriginalListOfParticlesOut.end(), [](auto& a, auto& b) -> bool{return a->momentum().pz() > b->momentum().pz(); });
184  std::sort(ResetListOfParticlesOut.begin(), ResetListOfParticlesOut.end(), [](auto& a, auto& b) -> bool{return a->momentum().pz() > b->momentum().pz(); });
185 
186  for ( std::vector<HepMC::ConstGenParticlePtr>::iterator originalPartOutIter(OriginalListOfParticlesOut.begin()), resetPartOutIter(ResetListOfParticlesOut.begin());
187  originalPartOutIter != OriginalListOfParticlesOut.end() && resetPartOutIter != ResetListOfParticlesOut.end();
188  ++originalPartOutIter, ++resetPartOutIter
189  ) {
190  if (compareGenParticle(*originalPartOutIter,*resetPartOutIter).isFailure()) {
191  ATH_MSG_ERROR ( "output particle properties differ!" );
192  pass &= false;
193  }
194  ATH_MSG_VERBOSE ( "particles match!" );
195  }
196  if(!pass) { return StatusCode::FAILURE; }
197  return StatusCode::SUCCESS;
198 }
199 
200 
201 #else
202 StatusCode TruthClosureCheck::compareGenVertex(const HepMC::GenVertex& origVertex,
203  const HepMC::GenVertex& resetVertex) const
204 {
205  bool pass{true};
206 
207  if (HepMC::barcode(origVertex) != HepMC::barcode(resetVertex)) { // FIXME HepMC::barcode
208  ATH_MSG_ERROR ("vertex barcode differs! Original: " << HepMC::barcode(origVertex) << ", Reset: " << HepMC::barcode(resetVertex));
209  pass = false;
210  }
211  if (origVertex.position() != resetVertex.position()) {
212  const HepMC::FourVector& oP = origVertex.position();
213  const HepMC::FourVector& rP = resetVertex.position();
214  ATH_MSG_ERROR("vertex position differs! Original: ("<<oP.x()<<","<<oP.y()<<","<<oP.z()<<","<<oP.t()<<"), Reset: ("<<rP.x()<<","<<rP.y()<<","<<rP.z()<<","<<rP.t()<<")");
215  pass = false;
216  }
217  if (origVertex.particles_in_size() != resetVertex.particles_in_size() ) {
218  ATH_MSG_ERROR ("particles_in_size differs! Original: "<<origVertex.particles_in_size()<<", Reset: "<<resetVertex.particles_in_size());
219  pass = false;
220  }
221  if (origVertex.particles_out_size() != resetVertex.particles_out_size() ) {
222  ATH_MSG_ERROR ("particles_out_size differs! Original: "<<origVertex.particles_out_size()<<", Reset: "<<resetVertex.particles_out_size());
223  pass = false;
224  }
225  HepMC::GenVertex::particles_in_const_iterator originalPartInIter(origVertex.particles_in_const_begin());
226  const HepMC::GenVertex::particles_in_const_iterator endOfOriginalListOfParticlesIn(origVertex.particles_in_const_end());
227  HepMC::GenVertex::particles_in_const_iterator resetPartInIter(resetVertex.particles_in_const_begin());
228  const HepMC::GenVertex::particles_in_const_iterator endOfResetListOfParticlesIn(resetVertex.particles_in_const_end());
229  while( originalPartInIter!=endOfOriginalListOfParticlesIn &&
230  resetPartInIter!=endOfResetListOfParticlesIn ) {
231  if (compareGenParticle(**originalPartInIter,**resetPartInIter).isFailure()) {
232  ATH_MSG_ERROR ( "input particle properties differ!" );
233  pass &= false;
234  }
235  ATH_MSG_VERBOSE ( "particles match!" );
236  ++resetPartInIter;
237  ++originalPartInIter;
238  }
239 
240  HepMC::GenVertex::particles_out_const_iterator originalPartOutIter(origVertex.particles_out_const_begin());
241  const HepMC::GenVertex::particles_out_const_iterator endOfOriginalListOfParticlesOut(origVertex.particles_out_const_end());
242  // ordering of particles may differ in each case - sigh..
243  const HepMC::GenVertex::particles_out_const_iterator endOfResetListOfParticlesOut(resetVertex.particles_out_const_end());
244  while( originalPartOutIter!=endOfOriginalListOfParticlesOut) {
245  const int barcodeOrig{HepMC::barcode(*originalPartOutIter)}; // FIXME barcode-based
246  HepMC::GenVertex::particles_out_const_iterator resetPartOutIter(resetVertex.particles_out_const_begin());
247  HepMC::GenVertex::particles_in_const_iterator matchingResetParticleIter{endOfResetListOfParticlesOut};
248  while(resetPartOutIter!=endOfResetListOfParticlesOut) {
249  if ( barcodeOrig == HepMC::barcode(*resetPartOutIter) ) {
250  matchingResetParticleIter = resetPartOutIter;
251  break;
252  }
253  ++resetPartOutIter;
254  }
255  if (matchingResetParticleIter==endOfResetListOfParticlesOut ||
256  compareGenParticle(**originalPartOutIter,**matchingResetParticleIter).isFailure()) {
257  ATH_MSG_ERROR ( "output particle properties differ!" );
258  pass &= false;
259  }
260  ATH_MSG_VERBOSE ( "particles match!" );
261  ++originalPartOutIter;
262  }
263 
264  if(!pass) { return StatusCode::FAILURE; }
265  return StatusCode::SUCCESS;
266 }
267 #endif
268 
269 
270 StatusCode TruthClosureCheck::compareMomenta(const HepMC::FourVector& origMomenta,
271  const HepMC::FourVector& resetMomenta) const
272 {
273  bool pass{true};
274  if (m_compareMomenta) {
275  if (m_momentaLimit<std::abs(origMomenta.px()-resetMomenta.px())) {
276  pass &= false;
277  }
278  if (m_momentaLimit<std::abs(origMomenta.py()-resetMomenta.py())) {
279  pass &= false;
280  }
281  if (m_momentaLimit<std::abs(origMomenta.pz()-resetMomenta.pz())) {
282  pass &= false;
283  }
284  if (m_momentaLimit<std::abs(origMomenta.e()-resetMomenta.e())) {
285  pass &= false;
286  }
287  if(origMomenta != resetMomenta) {
288  ATH_MSG_ERROR ("Exact momenta agreement check failed.");
289  pass &= false;
290  }
291  }
292  if(!pass) { return StatusCode::FAILURE; }
293  return StatusCode::SUCCESS;
294 }
295 
296 #ifdef HEPMC3
297 StatusCode TruthClosureCheck::compareGenParticle(const HepMC::ConstGenParticlePtr& origParticle,
298  const HepMC::ConstGenParticlePtr& resetParticle) const
299 {
300  if (!origParticle && !resetParticle) return StatusCode::SUCCESS;
301  if (!origParticle || !resetParticle) return StatusCode::FAILURE;
302 
303  bool pass{true};
304  if (HepMC::barcode(origParticle) != HepMC::barcode(resetParticle)) { // FIXME barcode-based
305  ATH_MSG_ERROR ("particle barcode differs! Original: "<<HepMC::barcode(origParticle)<<", Reset: "<<HepMC::barcode(resetParticle));
306  pass &= false;
307  }
308  if (origParticle->status() != resetParticle->status()) {
309  ATH_MSG_ERROR ("particle status differs! Original: "<<origParticle->status()<<", Reset: "<<resetParticle->status());
310  pass &= false;
311  }
312  if (origParticle->pdg_id() != resetParticle->pdg_id()) {
313  ATH_MSG_ERROR ("particle pdg_id differs! Original: "<<origParticle->pdg_id()<<", Reset: "<<resetParticle->pdg_id());
314  pass &= false;
315  }
316  if (compareMomenta(origParticle->momentum(), resetParticle->momentum()).isFailure()) {
317  const HepMC::FourVector& oM = origParticle->momentum();
318  const HepMC::FourVector& rM = resetParticle->momentum();
319  ATH_MSG_DEBUG("particle momentum differs! Original: ("<<oM.x()<<","<<oM.y()<<","<<oM.z()<<","<<oM.t()<<"), Reset: ("<<rM.x()<<","<<rM.y()<<","<<rM.z()<<","<<rM.t()<<")");
320  pass &= false;
321  }
322  if(!pass) { return StatusCode::FAILURE; }
323  return StatusCode::SUCCESS;
324 }
325 
326 #else
327 StatusCode TruthClosureCheck::compareGenParticle(const HepMC::GenParticle& origParticle,
328  const HepMC::GenParticle& resetParticle) const
329 {
330  bool pass{true};
331  if (HepMC::barcode(origParticle) != HepMC::barcode(resetParticle)) { // FIXME barcode-based
332  ATH_MSG_ERROR ("particle barcode differs! Original: "<<HepMC::barcode(origParticle)<<", Reset: "<<HepMC::barcode(resetParticle));
333  pass &= false;
334  }
335  if (origParticle.status() != resetParticle.status()) {
336  ATH_MSG_ERROR ("particle status differs! Original: "<<origParticle.status()<<", Reset: "<<resetParticle.status());
337  pass &= false;
338  }
339  if (origParticle.pdg_id() != resetParticle.pdg_id()) {
340  ATH_MSG_ERROR ("particle pdg_id differs! Original: "<<origParticle.pdg_id()<<", Reset: "<<resetParticle.pdg_id());
341  pass &= false;
342  }
343  if (compareMomenta(origParticle.momentum(), resetParticle.momentum()).isFailure()) {
344  const HepMC::FourVector& oM = origParticle.momentum();
345  const HepMC::FourVector& rM = resetParticle.momentum();
346  ATH_MSG_DEBUG("particle momentum differs! Original: ("<<oM.x()<<","<<oM.y()<<","<<oM.z()<<","<<oM.t()<<"), Reset: ("<<rM.x()<<","<<rM.y()<<","<<rM.z()<<","<<rM.t()<<")");
347  pass &= false;
348  }
349  if(!pass) { return StatusCode::FAILURE; }
350  return StatusCode::SUCCESS;
351 }
352 #endif
353 
354 #ifdef HEPMC3
355 //-------------------------------------------------
356 StatusCode TruthClosureCheck::execute() {
357 //-------------------------------------------------
358 
359  ATH_MSG_DEBUG( " execute..... " );
360  SG::ReadHandle<McEventCollection> originalMcEventCollection(m_originalMcEventCollection);
361  if (!originalMcEventCollection.isValid()) {
362  ATH_MSG_ERROR("Could not find original McEventCollection called " << originalMcEventCollection.name() << " in store " << originalMcEventCollection.store() << ".");
363  return StatusCode::FAILURE;
364  }
365  const HepMC::GenEvent& originalEvent(**(originalMcEventCollection->begin()));
366 
367  if (sanityCheck(originalEvent).isFailure()) {
368  ATH_MSG_FATAL("Problems in original GenEvent - bailing out.");
369  return StatusCode::FAILURE;
370  }
371 
372  SG::ReadHandle<McEventCollection> resetMcEventCollection(m_resetMcEventCollection);
373  if (!resetMcEventCollection.isValid()) {
374  ATH_MSG_ERROR("Could not find reset McEventCollection called " << resetMcEventCollection.name() << " in store " << resetMcEventCollection.store() << ".");
375  return StatusCode::FAILURE;
376  }
377  const HepMC::GenEvent& resetEvent(**(resetMcEventCollection->begin()));
378 
379  if (sanityCheck(resetEvent).isFailure()) {
380  ATH_MSG_FATAL("Problems in reset GenEvent - bailing out.");
381  return StatusCode::FAILURE;
382  }
383 
384  if (originalEvent.event_number() != resetEvent.event_number() ) {
385  ATH_MSG_ERROR ("event_number differs! Original: "<<originalEvent.event_number()<<", Reset: "<<resetEvent.event_number());
386  }
387 
388  if (HepMC::signal_process_id(originalEvent) != HepMC::signal_process_id(resetEvent) ) {
389  ATH_MSG_ERROR ("signal_process_id differs! Original: "<<HepMC::signal_process_id(originalEvent)<<", Reset: "<<HepMC::signal_process_id(resetEvent));
390  }
392  std::vector<HepMC::ConstGenParticlePtr> OriginalBeams = originalEvent.beams();
393  std::vector<HepMC::ConstGenParticlePtr> ResetBeams = resetEvent.beams();
394 
395  for ( std::vector<HepMC::ConstGenParticlePtr>::iterator originalBeamIter(OriginalBeams.begin()), resetBeamIter(ResetBeams.begin());
396  originalBeamIter != OriginalBeams.end() && resetBeamIter != ResetBeams.end();
397  ++originalBeamIter, ++resetBeamIter
398  ) {
399  if (compareGenParticle(*originalBeamIter,*resetBeamIter).isFailure()) {
400  ATH_MSG_ERROR ( "Beam particle properties differ!" );
401  }
403  if (originalEvent.particles().size() != resetEvent.particles().size() ) {
404  ATH_MSG_ERROR ("particles_size differs! Original: "<<originalEvent.particles().size()<<", Reset: "<<resetEvent.particles().size());
405  }
406 
407  if (originalEvent.vertices().size() != resetEvent.vertices().size() ) {
408  ATH_MSG_ERROR ("vertices_size differs! Original: "<<originalEvent.vertices().size()<<", Reset: "<<resetEvent.vertices().size());
409  }
410  }
411 
412  //loop over vertices in Background GenEvent
413  std::vector<HepMC::ConstGenVertexPtr> OriginalListOfVertices = originalEvent.vertices();
414  std::vector<HepMC::ConstGenVertexPtr> ResetListOfVertices = resetEvent.vertices();
415 
416  for( std::vector<HepMC::ConstGenVertexPtr>::iterator origVertexIter(OriginalListOfVertices.begin()),resetVertexIter(ResetListOfVertices.begin());
417  origVertexIter != OriginalListOfVertices.end() && resetVertexIter != ResetListOfVertices.end();
418  ++origVertexIter, ++resetVertexIter
419  ) {
420  if (compareGenVertex(*origVertexIter,*resetVertexIter).isFailure()) {
421  ATH_MSG_ERROR ( "vertices differ!" );
422  printGenVertex(*origVertexIter,*resetVertexIter);
423  }
424  }
425  ATH_MSG_INFO( "Completed Truth reset closure check ..... " );
426 
427  return StatusCode::SUCCESS;
428 
429 }
430 
431 #else
432 
433 
434 //-------------------------------------------------
435 StatusCode TruthClosureCheck::execute() {
436  //-------------------------------------------------
437 
438  ATH_MSG_DEBUG( " execute..... " );
439  SG::ReadHandle<McEventCollection> originalMcEventCollection(m_originalMcEventCollection);
440  if (!originalMcEventCollection.isValid()) {
441  ATH_MSG_ERROR("Could not find original McEventCollection called " << originalMcEventCollection.name() << " in store " << originalMcEventCollection.store() << ".");
442  return StatusCode::FAILURE;
443  }
444  const HepMC::GenEvent& originalEvent(**(originalMcEventCollection->begin()));
445 
446  if (sanityCheck(originalEvent).isFailure()) {
447  ATH_MSG_FATAL("Problems in original GenEvent - bailing out.");
448  return StatusCode::FAILURE;
449  }
450 
451  SG::ReadHandle<McEventCollection> resetMcEventCollection(m_resetMcEventCollection);
452  if (!resetMcEventCollection.isValid()) {
453  ATH_MSG_ERROR("Could not find reset McEventCollection called " << resetMcEventCollection.name() << " in store " << resetMcEventCollection.store() << ".");
454  return StatusCode::FAILURE;
455  }
456  const HepMC::GenEvent& resetEvent(**(resetMcEventCollection->begin()));
457 
458  if (sanityCheck(resetEvent).isFailure()) {
459  ATH_MSG_FATAL("Problems in reset GenEvent - bailing out.");
460  return StatusCode::FAILURE;
461  }
462 
463  if (originalEvent.event_number() != resetEvent.event_number() ) {
464  ATH_MSG_ERROR ("event_number differs! Original: "<<originalEvent.event_number()<<", Reset: "<<resetEvent.event_number());
465  }
466 
467  if (originalEvent.signal_process_id() != resetEvent.signal_process_id() ) {
468  ATH_MSG_ERROR ("signal_process_id differs! Original: "<<originalEvent.signal_process_id()<<", Reset: "<<resetEvent.signal_process_id());
469  }
470 
471  if (originalEvent.valid_beam_particles() != resetEvent.valid_beam_particles() ) {
472  ATH_MSG_ERROR ("valid_beam_particles differs! Original: "<<originalEvent.valid_beam_particles()<<", Reset: "<<resetEvent.valid_beam_particles());
473  }
474  else if (originalEvent.valid_beam_particles() && resetEvent.valid_beam_particles()) {
475  std::pair<HepMC::GenParticle*,HepMC::GenParticle*> originalBP = originalEvent.beam_particles();
476  std::pair<HepMC::GenParticle*,HepMC::GenParticle*> resetBP = resetEvent.beam_particles();
477  if ( ( !originalBP.first && resetBP.first ) ||
478  ( originalBP.first && !resetBP.first ) ||
479  ( originalBP.first && resetBP.first && compareGenParticle(*(originalBP.first), *(resetBP.first)).isFailure() ) ) {
480  ATH_MSG_ERROR ( "First beam particle does not match" );
481  }
482  if ( ( !originalBP.second && resetBP.second ) ||
483  ( originalBP.second && !resetBP.second ) ||
484  ( originalBP.second && resetBP.second && compareGenParticle(*(originalBP.second), *(resetBP.second)).isFailure() ) ) {
485  ATH_MSG_ERROR ( "Second beam particle does not match" );
486  }
487 
488  if (originalEvent.particles_size() != resetEvent.particles_size() ) {
489  ATH_MSG_ERROR ("particles_size differs! Original: "<<originalEvent.particles_size()<<", Reset: "<<resetEvent.particles_size());
490  }
491 
492  if (originalEvent.vertices_size() != resetEvent.vertices_size() ) {
493  ATH_MSG_ERROR ("vertices_size differs! Original: "<<originalEvent.vertices_size()<<", Reset: "<<resetEvent.vertices_size());
494  }
495  }
496 
497  //loop over vertices in Background GenEvent
498  HepMC::GenEvent::vertex_const_iterator origVertexIter(originalEvent.vertices_begin());
499  const HepMC::GenEvent::vertex_const_iterator endOfOriginalListOfVertices(originalEvent.vertices_end());
500  HepMC::GenEvent::vertex_const_iterator resetVertexIter(resetEvent.vertices_begin());
501  const HepMC::GenEvent::vertex_const_iterator endOfResetListOfVertices(resetEvent.vertices_end());
502  while( origVertexIter!=endOfOriginalListOfVertices &&
503  resetVertexIter!=endOfResetListOfVertices) {
504  if (compareGenVertex(**origVertexIter,**resetVertexIter).isFailure()) {
505  ATH_MSG_ERROR ( "vertices differ!" );
506  printGenVertex(**origVertexIter,**resetVertexIter);
507  }
508  ++origVertexIter;
509  ++resetVertexIter;
510  }
511  ATH_MSG_INFO( "Completed Truth reset closure check ..... " );
512 
513  return StatusCode::SUCCESS;
514 
515 }
516 #endif
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
TruthClosureCheck::initialize
virtual StatusCode initialize() override final
Definition: TruthClosureCheck.cxx:33
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
TruthClosureCheck::m_postSimulation
bool m_postSimulation
Definition: TruthClosureCheck.h:47
TruthClosureCheck::m_originalMcEventCollection
SG::ReadHandleKey< McEventCollection > m_originalMcEventCollection
Definition: TruthClosureCheck.h:43
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
TruthClosureCheck::m_compareMomenta
bool m_compareMomenta
Definition: TruthClosureCheck.h:45
GenEvent.h
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
SG::ReadHandle< McEventCollection >
TruthClosureCheck::sanityCheck
StatusCode sanityCheck(const HepMC::GenEvent &event) const
Definition: TruthClosureCheck.cxx:43
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
SG::VarHandleBase::name
const std::string & name() const
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleBase.cxx:75
HepMC::signal_process_id
int signal_process_id(const GenEvent &e)
Definition: GenEvent.h:635
GenVertex.h
TruthClosureCheck::compareGenParticle
StatusCode compareGenParticle(const HepMC::GenParticle &origParticle, const HepMC::GenParticle &resetParticle) const
Definition: TruthClosureCheck.cxx:327
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
TruthClosureCheck::compareMomenta
StatusCode compareMomenta(const HepMC::FourVector &origMomenta, const HepMC::FourVector &resetMomenta) const
Definition: TruthClosureCheck.cxx:270
GeoPrimitives.h
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
event
POOL::TEvent event(POOL::TEvent::kClassAccess)
McEventCollection.h
TruthClosureCheck::TruthClosureCheck
TruthClosureCheck(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TruthClosureCheck.cxx:22
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
SG::VarHandleBase::store
std::string store() const
Return the name of the store holding the object we are proxying.
Definition: StoreGate/src/VarHandleBase.cxx:376
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TruthClosureCheck::execute
virtual StatusCode execute() override final
Definition: TruthClosureCheck.cxx:435
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
AthAlgorithm
Definition: AthAlgorithm.h:47
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition: GenParticle.h:38
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
TruthClosureCheck::m_momentaLimit
double m_momentaLimit
Definition: TruthClosureCheck.h:46
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
TruthClosureCheck::printGenVertex
void printGenVertex(const HepMC::GenVertex &origVertex, const HepMC::GenVertex &resetVertex) const
Definition: TruthClosureCheck.cxx:98
TruthClosureCheck::compareGenVertex
StatusCode compareGenVertex(const HepMC::GenVertex &origVertex, const HepMC::GenVertex &resetVertex) const
Definition: TruthClosureCheck.cxx:202
MC::isStable
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
Definition: HepMCHelpers.h:45
a
TList * a
Definition: liststreamerinfos.cxx:10
std::sort
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:623
MC::isDecayed
bool isDecayed(const T &p)
Identify if the particle decayed.
Definition: HepMCHelpers.h:42
HepMC::ConstGenVertexPtr
const HepMC::GenVertex * ConstGenVertexPtr
Definition: GenVertex.h:60
TruthClosureCheck::m_resetMcEventCollection
SG::ReadHandleKey< McEventCollection > m_resetMcEventCollection
Definition: TruthClosureCheck.h:44
ReadHandle.h
Handle class for reading from StoreGate.
TruthClosureCheck.h
HepMCHelpers.h
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
GenParticle
@ GenParticle
Definition: TruthClasses.h:30
Generated on Thu Nov 7 2024 21:31:05 for ATLAS Offline Software by   doxygen 1.8.18