ATLAS Offline Software
TestHepMC.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 #ifndef XAOD_ANALYSIS
5 
6 #include "EvgenProdTools/TestHepMC.h"
7 #include "GaudiKernel/DataSvc.h"
8 #include "GaudiKernel/SystemOfUnits.h"
9 #include "TruthUtils/HepMCHelpers.h"
10 #include "PathResolver/PathResolver.h"
11 
12 // For find
13 #include <algorithm>
14 
15 
16 
17 TestHepMC::TestHepMC(const std::string& name, ISvcLocator* pSvcLocator)
18  : GenBase(name, pSvcLocator),
19  m_thistSvc("THistSvc", name)
20 {
21  declareProperty("MaxLoops", m_maxloops = -1); //< Maximal number of particles allowed in the loops. -1 == any number
22  declareProperty("PdgToSearch", m_pdg = 15); //< @todo This test is a bit weirdly specific to taus
23  declareProperty("CmEnergy", m_cm_energy = -1); // in MeV, -1 = get from event
24  declareProperty("MaxTransVtxDisp", m_max_dist_trans = 100.); // mm
25  declareProperty("MaxVtxDisp", m_max_dist = 1000.); // mm;
26  declareProperty("EnergyDifference", m_energy_diff = 1000.); // MeV
27  declareProperty("EnergyDifferenceError", m_max_energy_diff = 100000.); // MeV
28  declareProperty("CmeDifference", m_cme_diff = 1.); // MeV
29  declareProperty("DumpEvent", m_dumpEvent = false);
30  declareProperty("MinTau", m_min_tau = 1/300.); // ns; corresponds to 1mm
31  declareProperty("MaxNonG4Energy", m_nonG4_energy_threshold = 100.); //MeV
32  declareProperty("TauEffThreshold", m_tau_eff_threshold = 0.1); // fraction
33  declareProperty("EffWarnThreshold", m_eff_warn_threshold=0.99); // fraction
34  declareProperty("EffFailThreshold", m_eff_fail_threshold=0.98); // fraction
35  declareProperty("AccuracyMargin", m_accur_margin=0.); //MeV
36 
37  declareProperty("G4ExtraWhiteFile", m_paramFile = "g4_extrawhite.param" );
38  // a list of allowed pdgid which however might not follow the official rules
39  declareProperty("UnknownPDGIDFile", m_unknownPDGIDFile = "pdgid_extras.txt" );
40 
41  declareProperty("NoDecayVertexStatuses", m_vertexStatuses );
42 
43  declareProperty("BeamEnergyTest", m_beamEnergyTest = true); //switching off inactive
44  declareProperty("VtxNaNTest", m_vtxNaNTest = true);
45  declareProperty("VtxDisplacedTest", m_vtxDisplacedTest = true);
46  declareProperty("MomNaNTest", m_momNaNTest = true);
47  declareProperty("LifeTimeTest", m_lifeTimeTest = true);
48  declareProperty("EnergyG4Test", m_energyG4Test = true);
49  declareProperty("EnergyImbalanceTest", m_energyImbalanceTest = true);
50  declareProperty("MomImbalanceTest", m_momImbalanceTest = true);
51  declareProperty("NegativeEnergyTest", m_negativeEnergyTest = true);
52  declareProperty("TachyonsTest", m_tachyonsTest = true);
53  declareProperty("UnstableNoVtxTest", m_unstableNoVtxTest = true);
54  declareProperty("Pi0NoVtxTest", m_pi0NoVtxTest = true);
55  declareProperty("UndisplacedDaughtersTest", m_undisplacedDaughtersTest = true);
56  declareProperty("UknownPDGIDTest", m_unknownPDGIDTest = true);
57  declareProperty("AllowMissingCrossSection", m_allowMissingXSec = false);
58 
59  m_vertexStatuses.push_back( 1 );
60  m_vertexStatuses.push_back( 3 );
61  m_vertexStatuses.push_back( 4 );
62 
63 
64  declareProperty("THistSvc", m_thistSvc);
65 
66  declareProperty("DoHist", m_doHist=false); //histograming yes/no true/false
67 
68  m_nPass = 0;
69  m_nFail = 0;
70 
71  m_TotalTaus = 0;
72  m_FastDecayedTau = 0;
73 
74  // Check counters
75  m_invalidBeamParticlesCheckRate = 0;
76  m_beamParticleswithStatusNotFourCheckRate = 0;
77  m_beamEnergyCheckRate = 0;
78  m_vtxNANandINFCheckRate = 0;
79  m_vtxDisplacedstatuscode12CheckRate = 0;
80  m_vtxDisplacedstatuscodenot12CheckRate = 0;
81  m_partMomentumNANandINFCheckRate = 0;
82  m_undecayedPi0statuscode12CheckRate = 0;
83  m_unstableNoEndVtxCheckRate = 0;
84  m_negativeEnergyTachyonicCheckRate = 0;
85  m_decayCheckRate = 0;
86  m_undisplacedLLHdaughtersCheckRate = 0;
87  m_nonZeroPhotonMassCheckRate = 0;
88  m_energyBalanceCheckRate = 0;
89  m_momentumBalanceCheckRate = 0;
90  m_negativeEnergyCheckRate = 0;
91  m_tachyonCheckRate = 0;
92  m_stableUnstableNoParentCheckRate = 0;
93  m_unstablePartNoDecayVtxCheckRate = 0;
94  m_undecayedPi0CheckRate = 0;
95  m_Status1ShortLifetime = 0;
96  m_vtxDisplacedMoreThan_1m_CheckRate = 0;
97  m_undisplacedDecayDaughtersOfDisplacedVtxCheckRate = 0;
98  m_nonG4_energyCheckRate = 0;
99  m_unknownPDGIDCheckRate = 0;
100  m_noXSECset = 0;
101 
102  m_h_energy_dispVtxCheck = 0;
103  m_h_energy_dispVtxCheck_lt10 = 0;
104  m_h_pdgid_dispVtxCheck = 0;
105  m_h_status_dispVtxCheck = 0;
106  m_h_px_dispVtxCheck = 0;
107  m_h_py_dispVtxCheck = 0;
108  m_h_pz_dispVtxCheck = 0;
109  m_h_vx_dispVtxCheck = 0;
110  m_h_vy_dispVtxCheck = 0;
111  m_h_vz_dispVtxCheck = 0;
112  m_h_vxprod_dispVtxCheck = 0;
113  m_h_vyprod_dispVtxCheck = 0;
114  m_h_vzprod_dispVtxCheck = 0;
115  m_h_vtxprod_dispVtxCheck = 0;
116  m_h_vtxend_dispVtxCheck = 0;
117  m_h_photon_mass = 0;
118  m_h_photon_energy = 0;
119  m_h_photon_e2_p2_e2 = 0;
120  m_h_energyImbalance = 0;
121  m_h_momentumImbalance_px = 0;
122  m_h_momentumImbalance_py = 0;
123  m_h_momentumImbalance_pz = 0;
124  m_h_beamparticle1_Energy = 0;
125  m_h_beamparticle2_Energy = 0;
126  m_h_cmEnergyDiff = 0;
127 }
128 
129 
130 StatusCode TestHepMC::initialize() {
131  CHECK(GenBase::initialize());
132 
133  if (m_doHist){
134  CHECK(m_thistSvc.retrieve());
135 
136  m_h_energy_dispVtxCheck = new TH1F("h_energy_dispVtxCheck", "h_energy_dispVtxCheck", 2000, 0., 2000.);
137  m_h_energy_dispVtxCheck_lt10 = new TH1F("h_energy_dispVtxCheck_lt10", "h_energy_dispVtxCheck_lt10", 1000, 0., 10.);
138  m_h_pdgid_dispVtxCheck = new TH1F("h_pdgid_dispVtxCheck", "h_pdgid_dispVtxCheck", 10000, 0., 10000.);
139  m_h_status_dispVtxCheck = new TH1F("h_status_dispVtxCheck", "h_status_dispVtxCheck", 10000, 0., 10000.);
140  m_h_px_dispVtxCheck = new TH1F("h_px_dispVtxCheck", "h_px_dispVtxCheck", 4000, -2000., 2000.);
141  m_h_py_dispVtxCheck = new TH1F("h_py_dispVtxCheck", "h_py_dispVtxCheck", 4000, -2000., 2000.);
142  m_h_pz_dispVtxCheck = new TH1F("h_pz_dispVtxCheck", "h_pz_dispVtxCheck", 4000, -2000., 2000.);
143  m_h_vx_dispVtxCheck = new TH1F("h_vx_dispVtxCheck", "h_vx_dispVtxCheck", 40000, -200., 200);
144  m_h_vy_dispVtxCheck = new TH1F("h_vy_dispVtxCheck", "h_vy_dispVtxCheck", 40000, -200., 200);
145  m_h_vz_dispVtxCheck = new TH1F("h_vz_dispVtxCheck", "h_vz_dispVtxCheck", 40000, -200., 200);
146  m_h_vxprod_dispVtxCheck = new TH1F("h_vxprod_dispVtxCheck", "h_vxprod_dispVtxCheck", 40000, -200., 200.);
147  m_h_vyprod_dispVtxCheck = new TH1F("h_vyprod_dispVtxCheck", "h_vyprod_dispVtxCheck", 40000, -200., 200.);
148  m_h_vzprod_dispVtxCheck = new TH1F("h_vzprod_dispVtxCheck", "h_vzprod_dispVtxCheck", 40000, -200., 200.);
149  m_h_vtxprod_dispVtxCheck = new TH1F("h_vtxprod_dispVtxCheck", "h_vtxprod_dispVtxCheck", 20000, 0., 200.);
150  m_h_vtxend_dispVtxCheck = new TH1F("h_vtxend_dispVtxCheck", "h_vtxend_dispVtxCheck", 20000, 0., 200.);
151  m_h_photon_mass = new TH1F("h_photon_mass", "h_photon_mass", 20000, -10000., 10000);
152  m_h_photon_energy = new TH1F("h_photon_energy", "h_photon_energy", 20000, -10000., 10000);
153  m_h_photon_e2_p2_e2 = new TH1F("h_photon_e2_p2_e2", "h_photon_e2_p2_e2", 20000, -10., 10);
154  m_h_energyImbalance = new TH1F("h_energyImbalance", "h_energyImbalance", 2000, 0., 2000.);
155  m_h_momentumImbalance_px = new TH1F("h_momentumImbalance_px", "h_momentumImbalance_px", 2000,0., 2000.);
156  m_h_momentumImbalance_py = new TH1F("h_momentumImbalance_py", "h_momentumImbalance_py", 2000,0., 2000.);
157  m_h_momentumImbalance_pz = new TH1F("h_momentumImbalance_pz", "h_momentumImbalance_pz", 2000,0., 2000.);
158  m_h_beamparticle1_Energy = new TH1F("h_beamparticle1_Energy", "h_beamparticle1_Energy", 14000,0., 14000.);
159  m_h_beamparticle2_Energy = new TH1F("h_beamparticle2_Energy", "h_beamparticle2_Energy", 14000,0., 14000.);
160  m_h_cmEnergyDiff = new TH1F("h_cmEnergyDiff", "h_cmEnergyDiff", 8000, -4000., 4000.);
161 
162  CHECK(m_thistSvc->regHist("/TestHepMCname/h_energy_dispVtxCheck", m_h_energy_dispVtxCheck));
163  CHECK(m_thistSvc->regHist("/TestHepMCname/h_energy_dispVtxCheck_lt10", m_h_energy_dispVtxCheck_lt10));
164  CHECK(m_thistSvc->regHist("/TestHepMCname/h_pdgid_dispVtxCheck", m_h_pdgid_dispVtxCheck));
165  CHECK(m_thistSvc->regHist("/TestHepMCname/h_status_dispVtxCheck", m_h_status_dispVtxCheck));
166  CHECK(m_thistSvc->regHist("/TestHepMCname/h_px_dispVtxCheck", m_h_px_dispVtxCheck));
167  CHECK(m_thistSvc->regHist("/TestHepMCname/h_py_dispVtxCheck", m_h_py_dispVtxCheck));
168  CHECK(m_thistSvc->regHist("/TestHepMCname/h_pz_dispVtxCheck", m_h_pz_dispVtxCheck));
169  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vx_dispVtxCheck", m_h_vx_dispVtxCheck));
170  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vy_dispVtxCheck", m_h_vy_dispVtxCheck));
171  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vz_dispVtxCheck", m_h_vz_dispVtxCheck));
172  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vxprod_dispVtxCheck", m_h_vxprod_dispVtxCheck));
173  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vyprod_dispVtxCheck", m_h_vyprod_dispVtxCheck));
174  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vzprod_dispVtxCheck", m_h_vzprod_dispVtxCheck));
175  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vtxprod_dispVtxCheck", m_h_vtxprod_dispVtxCheck));
176  CHECK(m_thistSvc->regHist("/TestHepMCname/h_vtxend_dispVtxCheck", m_h_vtxend_dispVtxCheck));
177  CHECK(m_thistSvc->regHist("/TestHepMCname/h_photon_mass", m_h_photon_mass));
178  CHECK(m_thistSvc->regHist("/TestHepMCname/h_photon_energy", m_h_photon_energy));
179  CHECK(m_thistSvc->regHist("/TestHepMCname/h_photon_e2_p2_e2", m_h_photon_e2_p2_e2));
180  CHECK(m_thistSvc->regHist("/TestHepMCname/h_energyImbalance", m_h_energyImbalance));
181  CHECK(m_thistSvc->regHist("/TestHepMCname/h_momentumImbalance_px", m_h_momentumImbalance_px));
182  CHECK(m_thistSvc->regHist("/TestHepMCname/h_momentumImbalance_py", m_h_momentumImbalance_py));
183  CHECK(m_thistSvc->regHist("/TestHepMCname/h_momentumImbalance_pz", m_h_momentumImbalance_pz));
184  CHECK(m_thistSvc->regHist("/TestHepMCname/h_beamparticle1_Energy", m_h_beamparticle1_Energy));
185  CHECK(m_thistSvc->regHist("/TestHepMCname/h_beamparticle2_Energy", m_h_beamparticle2_Energy));
186  CHECK(m_thistSvc->regHist("/TestHepMCname/h_cmEnergyDiff", m_h_cmEnergyDiff));
187 
188  ATH_MSG_INFO("No decay vertex - is ignored for particles with status (list):" );
189  for ( unsigned int i = 0; i < m_vertexStatuses.size(); i++ ) ATH_MSG_INFO(" : " << m_vertexStatuses.at(i) );
190  ATH_MSG_INFO("Vertex statuses finished");
191 
192  } // End of histogramming setup
193 
194  // open the files and read G4particle_acceptlist.txt
195 
196  const std::string& fileLocation = PathResolverFindDataFile ( "G4particle_acceptlist.txt" );
197  std::ifstream G4file;
198  G4file.open(fileLocation);
199  std::string line;
200  int G4pdgID;
201  if (!G4file.fail()){
202  while(std::getline(G4file,line)){
203  std::stringstream ss(line);
204  ss >> G4pdgID;
205  m_G4pdgID_tab.push_back(G4pdgID);
206  }
207  G4file.close();
208  }
209  else {
210  ATH_MSG_WARNING("Failed to open G4particle_acceptlist.txt, checking that all particles are known by Genat4 cannot be performed");
211  }
212 
213  // Open the param file (G4 accept list)
214  G4file.open(m_paramFile.c_str());
215  if (!G4file.fail()){
216  ATH_MSG_INFO("extra accept list for G4 found " << m_paramFile.c_str());
217  while(std::getline(G4file,line)){
218  std::stringstream ss(line);
219  ss >> G4pdgID;
220  m_G4pdgID_tab.push_back(G4pdgID);
221  }
222  G4file.close();
223  }
224  else {
225  ATH_MSG_INFO("extra accept list for G4 not provided ");
226  }
227 
228  // Open the files and read susyParticlePdgid.txt
229  std::ifstream susyFile;
230  susyFile.open("susyParticlePdgid.txt");
231  int susyPdgID;
232  if (!susyFile.fail()){
233  while(getline(susyFile,line)){
234  std::stringstream ss1(line);
235  ss1 >> susyPdgID;
236  m_SusyPdgID_tab.push_back(susyPdgID);
237  }
238  susyFile.close();
239  }
240  else{
241  ATH_MSG_WARNING("Failed to open susyParticlePdgid.txt, listing particles not present in PDTTable");
242  }
243 
244  // Open the file of extra PDG IDs that don't need to obey the rules
245  std::ifstream pdgFile;
246  pdgFile.open(m_unknownPDGIDFile.c_str());
247  int pdgID;
248  if (!pdgFile.fail()){
249  ATH_MSG_INFO("extra accept list for PDG IDs found " << m_unknownPDGIDFile.c_str());
250  while(std::getline(pdgFile,line)){
251  std::stringstream ss(line);
252  ss >> pdgID;
253  m_uknownPDGID_tab.push_back(pdgID);
254  }
255  pdgFile.close();
256  }
257  else {
258  ATH_MSG_INFO("extra accept list for PDG IDs not provided");
259  }
260 
261  // Print Efficiency warning and error thresholds
262  ATH_MSG_INFO("EffWarnThreshold = " << m_eff_warn_threshold * 100 << " %");
263  ATH_MSG_INFO("EffFailThreshold = " << m_eff_fail_threshold * 100 << " %");
264 
265  return StatusCode::SUCCESS;
266 }
267 
268 
269 StatusCode TestHepMC::execute() {
270 
271  // Holder for filter outcome; allows us to check all filters on each event and diagnose multiple problems at once
272  bool filter_pass = true;
273 
274  // Loop over all events in McEventCollection
275  for(const HepMC::GenEvent* evt : *events_const()) {
276  double totalPx = 0;
277  double totalPy = 0;
278  double totalPz = 0;
279  double totalE = 0;
280  double nonG4_energy = 0;
281  std::vector<HepMC::ConstGenParticlePtr> negEnPart;
282  std::vector<HepMC::ConstGenParticlePtr> tachyons;
283  std::vector<HepMC::ConstGenParticlePtr> unstNoEnd;
284  std::vector<HepMC::ConstGenParticlePtr> unDecPi0;
285  std::vector<HepMC::ConstGenParticlePtr> undisplaceds;
286 
287  m_looper.findLoops(evt,true);
288  if (!m_looper.loop_particles().empty() || !m_looper.loop_vertices().empty()) {
289  ATH_MSG_DEBUG("Found " << m_looper.loop_vertices().size() << " vertices in loops");
290  ATH_MSG_DEBUG("Found " << m_looper.loop_particles().size() << " particles in loops");
291  ATH_MSG_DEBUG("Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
292  if (m_maxloops > 0 && m_looper.loop_particles().size() > static_cast<std::size_t>(m_maxloops) ) filter_pass = false;
293  }
294 
295 #ifdef HEPMC3
296  const auto xsec = evt->cross_section();
297  if (!xsec) {
298  ATH_MSG_WARNING("WATCH OUT: event is missing the generator cross-section!");
299  ++m_noXSECset;
300  if (m_allowMissingXSec) {
301  ATH_MSG_WARNING("-> Adding a dummy cross-section for debugging purposes.");
302  // for debugging purposes only -> set a dummy cross-section to make TestHepMC happy
303  std::shared_ptr<HepMC3::GenCrossSection> dummy_xsec = std::make_shared<HepMC3::GenCrossSection>();
304  dummy_xsec->set_cross_section(1.0,0.0);
305  HepMC::GenEvent* evt_nonconst = const_cast<HepMC::GenEvent*>(evt);
306  evt_nonconst->set_cross_section(std::move(dummy_xsec));
307  }
308  else {
309  ATH_MSG_WARNING("-> Will report this as failure.");
310  }
311  }
312 
313  // Check beams and work out per-event beam energy
314  std::vector<std::shared_ptr<const HepMC3::GenParticle>> beams_t;
315  for (auto p : evt->beams()) { if (p->status() == 4) beams_t.push_back(std::move(p)); }
316  std::pair<std::shared_ptr<const HepMC3::GenParticle>,std::shared_ptr<const HepMC3::GenParticle>> beams;
317  if (beams_t.size() == 2) {
318  beams.first=beams_t.at(0);
319  beams.second=beams_t.at(1);
320  } else {
321  ATH_MSG_WARNING("Invalid number of beam particles " << beams_t.size() << " this generator interface should be fixed");
323  for (const auto& part: beams_t) HepMC3::Print::line(part);
324  }
325 #else
326  auto beams = evt->beam_particles();
327 #endif
328  double cmenergy = m_cm_energy;
329  if (!HepMC::valid_beam_particles(evt)) {
330  ATH_MSG_WARNING("Invalid beam particle pointers -- this generator interface should be fixed");
331  if (cmenergy < 0) ATH_MSG_WARNING("Invalid expected beam energy: " << cmenergy << " MeV");
332  ++m_invalidBeamParticlesCheckRate;
333  } else {
334  if (!MC::isBeam(beams.first) || !MC::isBeam(beams.second)) {
335  ATH_MSG_WARNING("Beam particles have incorrectly set status -- this generator interface should be fixed");
336  ++m_beamParticleswithStatusNotFourCheckRate;
337  }
338  const double sumE = beams.first->momentum().e() + beams.second->momentum().e();
339  const double sumP = beams.first->momentum().pz() + beams.second->momentum().pz();
340  cmenergy = std::sqrt(sumE*sumE - sumP*sumP);
341 
342  if(beams.first->pdg_id() == MC::OXYGEN && beams.second->pdg_id() == MC::OXYGEN){//OO collisions
343  cmenergy /= MC::baryonNumber(MC::OXYGEN); // divided by the total number of nucleons per nucleus
344  }
345  if(beams.first->pdg_id() == MC::LEAD && beams.second->pdg_id() == MC::LEAD){//PbPb collisions
346  cmenergy /= MC::baryonNumber(MC::LEAD); // divided by the total number of nucleons per nucleus
347  }
348  if(beams.first->pdg_id() == MC::OXYGEN && beams.second->pdg_id() == MC::PROTON){//Op collisions
349  cmenergy = -2.0*beams.second->momentum().pz()*std::sqrt(static_cast<double>(MC::numberOfProtons(MC::OXYGEN))/MC::baryonNumber(MC::OXYGEN));
350  }
351  if(beams.first->pdg_id() == MC::PROTON && beams.second->pdg_id() == MC::OXYGEN){//pO collisions
352  cmenergy = 2.0*beams.first->momentum().pz()*std::sqrt(static_cast<double>(MC::numberOfProtons(MC::OXYGEN))/MC::baryonNumber(MC::OXYGEN));
353  }
354  if(beams.first->pdg_id() == MC::LEAD && beams.second->pdg_id() == MC::PROTON){//Pbp collisions
355  cmenergy = -2.0*beams.second->momentum().pz()*std::sqrt(static_cast<double>(MC::numberOfProtons(MC::LEAD))/MC::baryonNumber(MC::LEAD));
356  }
357  if(beams.first->pdg_id() == MC::PROTON && beams.second->pdg_id() == MC::LEAD){//pPb collisions
358  cmenergy = 2.0*beams.first->momentum().pz()*std::sqrt(static_cast<double>(MC::numberOfProtons(MC::LEAD))/MC::baryonNumber(MC::LEAD));
359  }
360  if(beams.first->pdg_id() == MC::OXYGEN && beams.second->pdg_id() == MC::HELIUM){//OHe collisions
361  cmenergy /= std::sqrt(static_cast<double>(MC::baryonNumber(MC::OXYGEN)*MC::baryonNumber(MC::HELIUM)));
362  }
363  if(beams.first->pdg_id() == MC::HELIUM && beams.second->pdg_id() == MC::OXYGEN){//HeO collisions
364  cmenergy /= std::sqrt(static_cast<double>(MC::baryonNumber(MC::OXYGEN)*MC::baryonNumber(MC::HELIUM)));
365  }
366 
367  if (m_cm_energy > 0 && std::abs(cmenergy - m_cm_energy) > m_cme_diff) {
368  ATH_MSG_FATAL("Beam particles have incorrect energy: " << m_cm_energy/Gaudi::Units::GeV << " GeV expected, vs. " << cmenergy/Gaudi::Units::GeV << " GeV found");
369  setFilterPassed(false);
370  if (m_doHist){
371  m_h_beamparticle1_Energy->Fill(beams.first->momentum().e()/Gaudi::Units::GeV);
372  m_h_beamparticle2_Energy->Fill(beams.second->momentum().e()/Gaudi::Units::GeV);
373  m_h_cmEnergyDiff->Fill((cmenergy-m_cm_energy)/Gaudi::Units::GeV);
374  }
375  ++m_beamEnergyCheckRate;
376  // Special case: this is so bad that we immediately fail out
377  return StatusCode::FAILURE;
378  }
379  }
380 
381  // Check vertices
382  int vtxDisplacedstatuscode12CheckRateCnt=0;
383  int vtxDisplacedstatuscodenot12CheckRateCnt=0;
384  int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
385 #ifdef HEPMC3
386  for (const auto& vtx: evt->vertices()) {
387 #else
388  for (auto vitr = evt->vertices_begin(); vitr != evt->vertices_end(); ++vitr ) {
389  const HepMC::GenVertex* vtx = *vitr;
390 #endif
391  const HepMC::FourVector pos = vtx->position();
392 
393  // Check for NaNs and infs in vertex position components
394  if ( std::isnan(pos.x()) || std::isinf(pos.x()) ||
395  std::isnan(pos.y()) || std::isinf(pos.y()) ||
396  std::isnan(pos.z()) || std::isinf(pos.z()) ) {
397  ATH_MSG_WARNING("NaN (Not A Number) or inf found in the event record vertex positions");
398 
399  ++m_vtxNANandINFCheckRate;
400  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
401  if (m_vtxNaNTest) {
402  filter_pass = false;
403  }
404  } // Done of checking for nans and infinities
405 
406  // Check for too-far-displaced vertices
407  // Anything which propagates macroscopically should be set stable in evgen for G4 to handle
408  const double dist_trans2 = pos.x()*pos.x() + pos.y()*pos.y(); // in mm2
409  const double dist2 = dist_trans2 + pos.z()*pos.z(); // in mm2
410  const double dist_trans = std::sqrt(dist_trans2); // in mm
411  const double dist = std::sqrt(dist2); // in mm
412  if (dist2 > m_max_dist*m_max_dist) {
413  ATH_MSG_WARNING("Found vertex position displaced by more than " << m_max_dist << "mm: " << dist << "mm");
414  ++vtxDisplacedMoreThan_1m_CheckRateCnt;
415 
416  if (m_vtxDisplacedTest) {
417  filter_pass = false;
418  }
419  }
420  if (dist_trans2 > m_max_dist_trans*m_max_dist_trans) {
421  ATH_MSG_WARNING("Found vertex position displaced by more than " << m_max_dist_trans << "mm in transverse distance: " << dist_trans << "mm");
422 
423 #ifdef HEPMC3
424  for (const auto& part: vtx->particles_in()) {
425 #else
426  for (auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
427  auto part=(*part_it);
428 #endif
429  if (m_dumpEvent){
430  ATH_MSG_WARNING("Outgoing particle : ");
431  HepMC::Print::line(msg( MSG::WARNING ).stream(),part);
432  }
433  ATH_MSG_WARNING("production vertex = " << part->production_vertex()->position().x() << ", " << part->production_vertex()->position().y() << ", " << part->production_vertex()->position().z());
434  ATH_MSG_WARNING("end vertex = " << part->end_vertex()->position().x() << ", " << part->end_vertex()->position().y() << ", " << part->end_vertex()->position().z());
435  if (m_dumpEvent) ATH_MSG_WARNING("parents info: ");
436  if (part->production_vertex()) {
437 #ifdef HEPMC3
438  for(const auto& p_parents: part->production_vertex()->particles_in()) {
439 #else
440  for(auto p_parents_it = part->production_vertex()->particles_in_const_begin(); p_parents_it != part->production_vertex()->particles_in_const_end(); ++p_parents_it) {
441  auto p_parents=(*p_parents_it);
442 #endif
443  if (m_dumpEvent){
444  msg(MSG::WARNING) << "\t";
445  HepMC::Print::line( msg( MSG::WARNING ).stream() , p_parents );
446  }
447  }
448  } // Done with fancy print
449 
450  if (part->status() == 1 || part->status() == 2){
451  vtxDisplacedstatuscode12CheckRateCnt += 1;
452  } else {
453  vtxDisplacedstatuscodenot12CheckRateCnt += 1;
454  }
455 
456  if (m_doHist){
457  m_h_energy_dispVtxCheck->Fill(part->momentum().e()/Gaudi::Units::GeV);
458  if (part->momentum().e()/Gaudi::Units::GeV < 10.) {
459  m_h_energy_dispVtxCheck_lt10->Fill(part->momentum().e()/Gaudi::Units::GeV);
460  }
461  m_h_pdgid_dispVtxCheck->Fill(part->pdg_id());
462  m_h_status_dispVtxCheck->Fill(part->status());
463  m_h_px_dispVtxCheck->Fill(part->momentum().px()/Gaudi::Units::GeV);
464  m_h_py_dispVtxCheck->Fill(part->momentum().py()/Gaudi::Units::GeV);
465  m_h_pz_dispVtxCheck->Fill(part->momentum().pz()/Gaudi::Units::GeV);
466  m_h_vx_dispVtxCheck->Fill(part->end_vertex()->position().x());
467  m_h_vy_dispVtxCheck->Fill(part->end_vertex()->position().y());
468  m_h_vz_dispVtxCheck->Fill(part->end_vertex()->position().z());
469  m_h_vxprod_dispVtxCheck->Fill(part->production_vertex()->position().x());
470  m_h_vyprod_dispVtxCheck->Fill(part->production_vertex()->position().y());
471  m_h_vzprod_dispVtxCheck->Fill(part->production_vertex()->position().z());
472  double endvx = part->end_vertex()->position().x();
473  double endvy = part->end_vertex()->position().y();
474  double endvz = part->end_vertex()->position().z();
475  double prodvx = part->production_vertex()->position().x();
476  double prodvy = part->production_vertex()->position().y();
477  double prodvz = part->production_vertex()->position().z();
478  double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
479  double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
480  m_h_vtxend_dispVtxCheck->Fill(enddis);
481  m_h_vtxprod_dispVtxCheck->Fill(proddis);
482  } // End of the filling of histograms for bad vertices
483  } // End of a loop over theparents of the bad vertex
484  } // Found a bad vertex
485  } // Loop over all vertices
486  if (vtxDisplacedstatuscode12CheckRateCnt>0) ++m_vtxDisplacedstatuscode12CheckRate;
487  if (vtxDisplacedstatuscodenot12CheckRateCnt>0) ++m_vtxDisplacedstatuscodenot12CheckRate;
488  if (vtxDisplacedMoreThan_1m_CheckRateCnt>0) ++m_vtxDisplacedMoreThan_1m_CheckRate;
489 
490  // Check particles
491  for (auto pitr: *evt) {
492 
493  // Local loop variables to clean up the check code
494  const HepMC::FourVector pmom = pitr->momentum();
495  const int pstatus = pitr->status();
496  const int ppdgid = pitr->pdg_id();
497  // Check for NaNs and infs in momentum components
498  if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
499  std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
500  std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
501  std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
502  ATH_MSG_WARNING("NaN (Not A Number) or inf found in the event record momenta");
503  ++m_partMomentumNANandINFCheckRate;
504 
505  if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
506  if (m_momNaNTest) {
507  filter_pass = false;
508  }
509  } // End of check for NaNs and infinities
510 
511  // Check for undecayed pi0s
512  if (MC::isStable(pstatus) || MC::isDecayed(pstatus)) {
513  if (ppdgid == 111 && !pitr->end_vertex() ) {
514  unDecPi0.push_back( pitr);
515  ++m_undecayedPi0CheckRate;
516  }
517  } // End of check for undecayed pi0s
518 
519  //check stable particle lifetimes
520  if (MC::isStable(pstatus)) {
521  const HepPDT::ParticleData* pd = particleData(ppdgid);
522  if (pd != NULL) {
523  double plifetime = pd->lifetime()*1e+12; // why lifetime doesn't come in common units???
524  if (plifetime != 0 && plifetime < m_min_tau) { // particles with infinite lifetime get a 0 in the PDT
525  ATH_MSG_WARNING("Stable particle found with lifetime = " << plifetime << "~ns!!");
526  if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
527 
528  ++m_Status1ShortLifetime;
529 
530  if (m_lifeTimeTest) {
531  filter_pass = false;
532  }
533  } // Particle did not have infinite lifetime
534  } // The particle has a data table (so a lifetime)
535  else{
536  int susyPart = 0;
537  std::vector<int>::size_type count = 0;
538  while (susyPart==0 && (count < m_SusyPdgID_tab.size() )){
539  // no warning for SUSY particles from the list susyParticlePdgid.txt
540  if (m_SusyPdgID_tab[count] == std::abs(ppdgid)) {
541  susyPart=1;
542  }
543  count++;
544  } // Look through the SUSY table to see if this one should be counted
545  if (susyPart==0){
546  ATH_MSG_WARNING("Stable particle not found in PDT, no lifetime check done");
547  if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
548  } // It's a SUSY particle -- skip the lifetime check
549  } // The particle has no data table
550  } // Test if the particle is stable
551 
552  //Check that stable particles are known by G4 or they are non-interacting
553  const MC::DecodedPID decodedPID(ppdgid);
554  const int first_dig = decodedPID(0);
555 
556  if (MC::isStable(pstatus) && (!pitr->end_vertex()) && (MC::isSimInteracting(pitr)) && (!MC::isNucleus(ppdgid)) && (first_dig != 9) ) {
557 
558  int known_byG4 = 0;
559  std::vector<int>::size_type count =0;
560 
561  while (known_byG4==0 && count < m_G4pdgID_tab.size()){
562  if(ppdgid == m_G4pdgID_tab[count]) known_byG4=1;
563  count++;
564  }
565  if(known_byG4==0){
566  nonG4_energy += pmom.e();
567  ATH_MSG_WARNING("Interacting particle not known by Geant4 with ID " << ppdgid);
568  }
569  } // End of check that stable particles are known to G4 or are non-interacting
570 
571  // Check for bad PDG IDs
572  if (!MC::isValid(ppdgid)){
573  ATH_MSG_DEBUG("Invalid PDG ID found: " << ppdgid);
574  if (m_unknownPDGIDTest && std::find(m_uknownPDGID_tab.begin(),m_uknownPDGID_tab.end(),ppdgid)==m_uknownPDGID_tab.end()){
575  ATH_MSG_WARNING("Invalid and unmasked PDG ID found: " << ppdgid);
576  filter_pass = false;
577  }
578  } // End of check for invalid PDG IDs
579 
580  // Check for unstables with no end vertex,
581  if (!pitr->end_vertex() && MC::isDecayed(pstatus)) {
582  unstNoEnd.push_back(pitr);
583  ++m_unstableNoEndVtxCheckRate;
584  } // End of check for unstable with no end vertex
585 
586  // Sum final state mom/energy, and note negative energy / tachyonic particles
587  // std::cout << "status " << pstatus << " e " << pmom.e() << " pz " << pmom.pz()<< std::endl;
588  if ( MC::isStable(pstatus) && !pitr->end_vertex() ) {
589  totalPx += pmom.px();
590  totalPy += pmom.py();
591  totalPz += pmom.pz();
592  totalE += pmom.e();
593  if (pmom.e() < 0) {
594  negEnPart.push_back(pitr);
595  ++m_negativeEnergyTachyonicCheckRate;
596  }
597  const double aener = std::abs(pmom.e());
598  if ( aener+m_accur_margin < std::abs(pmom.px()) || aener+m_accur_margin < std::abs(pmom.py()) || aener+m_accur_margin < std::abs(pmom.pz()) ) {
599  tachyons.push_back(pitr);
600  ++m_negativeEnergyTachyonicCheckRate;
601  }
602  } // End of sums for momentum and energy conservation
603 
604  // Decay checks (uses PdgToSearch attr value, for tau by default)
606  int tau_child = 0;
607  if (std::abs(ppdgid) == m_pdg && (MC::isStable(pstatus) || MC::isDecayed(pstatus))) {
608  ++m_TotalTaus;
609  auto vtx = pitr->end_vertex();
610  if (vtx) {
611  double p_energy = 0;
612 #ifdef HEPMC3
613  for (auto desc: HepMC::descendant_particles(vtx)) {
614 #else
615  for (auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
616  auto desc=(*desc_it);
617 #endif
618  if (std::abs(desc->pdg_id()) == m_pdg) tau_child = 1;
619  if ( MC::isStable(desc) ) p_energy += desc->momentum().e();
620  }
621  if (std::abs( p_energy - pmom.e()) > m_energy_diff && !tau_child) {
622  ATH_MSG_WARNING("Energy sum (decay products): "
623  << "Energy (original particle) > " << m_energy_diff << " MeV, "
624  << "Event #" << evt->event_number() << ", "
625  << "The original particle = " << pitr);
626  ++m_decayCheckRate;
627  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
628  }
629  //most taus should not decay immediately
630  const HepMC::FourVector tau_decaypos = vtx->position();
631  const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
632  //tau_child != 1 exclude cases in which a tau is copied to another vertex or emits a photon
633  if ((tau_displacement < 1.e-6) && (tau_child!=1)) ++m_FastDecayedTau;
634  } else {
635  ATH_MSG_WARNING("UNDECAYED PARTICLE WITH PDG_ID = " << m_pdg);
636  ++m_decayCheckRate;
637  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
638  }
639  } // End of checks for specific particle (tau by default)
640 
641  // Check for undisplaced decay daughters from long-lived hadrons
642  if (pitr->end_vertex()) {
643  auto decayvtx = pitr->end_vertex();
644  const HepMC::FourVector decaypos = decayvtx->position();
645  const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
646  if (displacement > 1e-6) {
647  for (auto ip: *decayvtx) {
648  const HepMC::FourVector pos2 = ip->production_vertex()->position();
649  const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
650  if (displacement2 < 1e-6) {
651  ATH_MSG_WARNING("Decay child " << ip << " from " << pitr
652  << " has undisplaced vertex (" << ip->production_vertex()
653  << " @ " << displacement2 << "mm) "
654  << " but parent vertex is displaced (" << decayvtx
655  << " @ " << displacement << "mm)");
656  undisplaceds.push_back(std::move(ip));
657  ++m_undisplacedLLHdaughtersCheckRate;
658  } // Check for displacement below 1 um
659  } // Loop over all particles coming from the decay vertex
660  } // Displacement of greater than 1 um
661  } // End of check for undisplaced decay daughters from long-lived hadrons
662 
663  // Check for photons with non-zero masses
665  if (MC::isPhoton(ppdgid) && MC::isStable(pstatus)) {
666  const double mass = pitr->generated_mass();
667  if (std::abs(mass) > 1.0) { // in MeV
668  ATH_MSG_WARNING("Photon with non-zero mass found! Mass: " << mass << " MeV" << pitr);
669  ++m_nonZeroPhotonMassCheckRate;
670  }
671  } // End check for photons with too-large a mass
672 
673  } // End of loop over particles in the event
674 
675  // Energy of interacting particles not known by Geant4
676  if(nonG4_energy > m_nonG4_energy_threshold) {
677  ATH_MSG_WARNING("The energy of interacting particles not known by Geant4 is = " << nonG4_energy << " MeV");
678  if (m_energyG4Test) {
679  filter_pass = false;
680  }
681  ++m_nonG4_energyCheckRate;
682  } // End of check for interacting particles not known by G4
683 
684  // Energy balance
685  double lostE = std::abs(totalE - cmenergy);
686  if (lostE > m_energy_diff) {
687  ATH_MSG_WARNING("ENERGY BALANCE FAILED : E-difference = " << lostE << " MeV");
688 
689  ATH_MSG_WARNING("balance " << totalPx << " " << totalPy << " " << totalPz << " " << totalE);
690 
691  if (m_doHist){
692  m_h_energyImbalance->Fill(lostE/Gaudi::Units::GeV);
693  }
694  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
695  if (m_energyImbalanceTest) {
696  filter_pass = false;
697  }
698  ++m_energyBalanceCheckRate;
699  } // End of energy balance check
700 
701  // Momentum balance
702  if ( std::abs(totalPx) > m_energy_diff || std::abs(totalPy) > m_energy_diff || std::abs(totalPz) > m_energy_diff ) {
703  ATH_MSG_WARNING("MOMENTUM BALANCE FAILED : SumPx = " << totalPx << " SumPy = " << totalPy << " SumPz = " << totalPz << " MeV");
704  if (m_doHist){
705  m_h_momentumImbalance_px->Fill(std::abs(totalPx)/Gaudi::Units::GeV);
706  m_h_momentumImbalance_py->Fill(std::abs(totalPy)/Gaudi::Units::GeV);
707  m_h_momentumImbalance_pz->Fill(std::abs(totalPz)/Gaudi::Units::GeV);
708  }
709  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
710  if (m_momImbalanceTest) {
711  filter_pass = false;
712  }
713  ++m_momentumBalanceCheckRate;
714  } // End of momentum balance check
715 
716  // Negative energy particles
717  if (!negEnPart.empty()) {
718  std::stringstream ss;
719  ss << "NEGATIVE ENERGY PARTICLES FOUND :";
720  for (const auto &b: negEnPart){
721  ss << " " << b;
722  }
723  ATH_MSG_WARNING(ss.str());
724  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
725  if (m_negativeEnergyTest) {
726  filter_pass = false;
727  }
728  ++m_negativeEnergyCheckRate;
729  } // End of negative energy particle chedk
730 
731  // Tachyons
732  if (!tachyons.empty()) {
733  std::stringstream ss;
734  ss << "PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
735  for (auto b: tachyons){
736  ss << " " << b;
737  }
738  ATH_MSG_WARNING(ss.str());
739  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
740  if (m_tachyonsTest) {
741  filter_pass = false;
742  }
743  ++m_energyBalanceCheckRate;
744  } // End of tachyon check
745 
746  // Unstable particles with no decay vertex
747  if (!unstNoEnd.empty()) {
748  std::stringstream ss;
749  ss << "Unstable particle with no decay vertex found: ";
750  for (auto b: unstNoEnd){
751  ss << " " << b;
752  }
753  ATH_MSG_WARNING(ss.str());
754  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
755  if (m_unstableNoVtxTest) {
756  filter_pass = false;
757  }
758  ++m_unstablePartNoDecayVtxCheckRate;
759  } // End of unstable particle with no decay vertex check
760 
761  // Undecayed pi0
762  if (!unDecPi0.empty()) {
763  std::stringstream ss;
764  ss << "pi0 with no decay vertex found:";
765  for (auto b: unDecPi0){
766  ss << " " << b;
767  }
768  ATH_MSG_WARNING(ss.str());
769  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
770  if (m_pi0NoVtxTest) {
771  filter_pass = false;
772  }
773  ++m_undecayedPi0CheckRate;
774  } // End of undecayed pi0 check
775 
776  // Undisplaced decay daughters of displaced vertices
777  if (!undisplaceds.empty()) {
778  std::stringstream ss{"Undisplaced decay vertices from displaced particle: "};
779  for (HepMC::ConstGenParticlePtr b: undisplaceds){
780  // coverity[COPY_INSTEAD_OF_MOVE]
781  ss << " " << b;
782  }
783  ATH_MSG_WARNING(ss.str());
784  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
785  if (m_undisplacedDaughtersTest) {
786  filter_pass = false;
787  }
788  ++m_undisplacedDecayDaughtersOfDisplacedVtxCheckRate;
789  } // End of undisplaced decay daughter of displaced vertices check
790 
791  } // End of loop over MCEventCollection
792 
793  // End of execution for each event - update filter value
794  if (!filter_pass){
795  setFilterPassed(false);
796  ++m_nFail;
797  } else {
798  ++m_nPass;
799  }
800 
801  // If the efficiency after 100 events is below 10%, there is an important bug going on:
802  // we fail the job immediately so it doesn't run for ever
803  const double tmp_efficiency = double(m_nPass) / double(m_nPass + m_nFail);
804  if ((m_nPass + m_nFail) > 100 && tmp_efficiency < 0.1) {
805  ATH_MSG_FATAL("The efficiency after " << m_nPass + m_nFail << " events is " << tmp_efficiency*100. << "% !!!");
806  return StatusCode::FAILURE;
807  }
808 
809  return StatusCode::SUCCESS;
810 }
811 
812 
813 StatusCode TestHepMC::finalize() {
814 
815  ATH_MSG_INFO("Events passed = " << m_nPass << ", Events Failed = " << m_nFail);
816  // Keep a denominator for all the fractions, to ensure we don't get a bunch of nans.
817  // If nothing passed or failed, all the other counters should be zero; this is just avoiding FPEs etc.
818  double denom = m_nPass + m_nFail > 0. ? m_nPass + m_nFail : 1.;
819 
820  ATH_MSG_INFO(" Event rate with invalid Beam Particles = " << m_invalidBeamParticlesCheckRate*100.0/denom << "% (not included in test efficiency)");
821  ATH_MSG_INFO(" Event rate with beam particles and status not equal to 4 = " << m_beamParticleswithStatusNotFourCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
822  ATH_MSG_INFO(" Event rate with incorrect beam particle energies = " << m_beamEnergyCheckRate*100.0/denom << "% (not included in test efficiency)");
823  ATH_MSG_INFO(" Event rate with NaN (Not A Number) or inf found in the event record vertex positions = " << m_vtxNANandINFCheckRate*100.0/denom << "%");
824  if (!m_vtxNaNTest) ATH_MSG_INFO(" The check for NaN or inf in vtx. record is switched off, so is not included in the final TestHepMC efficiency ");
825  ATH_MSG_INFO(" Event rate with vertices displaced more than " << m_max_dist_trans << "~mm in transverse direction for particles with status code other than 1 and 2 = " << m_vtxDisplacedstatuscodenot12CheckRate*100.0/denom << "% (not included in test efficiency)");
826  ATH_MSG_INFO(" Event rate with vertices displaced more than " << m_max_dist << "~mm = " << m_vtxDisplacedMoreThan_1m_CheckRate*100.0/denom << "%");
827  if (!m_vtxDisplacedTest) ATH_MSG_INFO(" The check for displaced vertices is switched off, so is not included in the final TestHepMC efficiency ");
828  ATH_MSG_INFO(" Event rate with NAN (Not A Number) or inf found in particle momentum values = " << m_partMomentumNANandINFCheckRate*100.0/denom << "%");
829  if (!m_momNaNTest) ATH_MSG_INFO(" The check for NaN/inf in momentum record is switched off, so is not included in the final TestHepMC efficiency ");
830  ATH_MSG_INFO(" Event rate with undecayed pi0's with status 1 or 2 = " << m_undecayedPi0statuscode12CheckRate*100.0/denom << "% (not included in test efficiency)");
831  ATH_MSG_INFO(" Event rate with unstable particles with no end vertex = " << m_unstableNoEndVtxCheckRate*100.0/denom << "% (not included in test efficiency)");
832  ATH_MSG_INFO(" Event rate with negative total energy like for tachyonic particles = " << m_negativeEnergyTachyonicCheckRate*100.0/denom << "% (not included in test efficiency)");
833  ATH_MSG_INFO(" Event rate with particles with improper decay properties = " << m_decayCheckRate*100.0/denom << "% (not included in test efficiency)");
834  ATH_MSG_INFO(" Event rate with undisplaced daughters of long lived hadrons = " << m_undisplacedLLHdaughtersCheckRate*100.0/denom << "% (not included in test efficiency)");
835  ATH_MSG_INFO(" Event rate with non zero photon mass = " << m_nonZeroPhotonMassCheckRate*100.0/denom << "% (not included in test efficiency)");
836  ATH_MSG_INFO(" Event rate with no energy balance = " << m_energyBalanceCheckRate*100.0/denom << "%");
837  if (!m_energyImbalanceTest) ATH_MSG_INFO(" The check for energy imbalance is switched off, so is not included in the final TestHepMC efficiency ");
838  ATH_MSG_INFO(" Event rate with no momentum balance = " << m_momentumBalanceCheckRate*100.0/denom << "%");
839  if (!m_momImbalanceTest) ATH_MSG_INFO(" The check for momentum imbalance is switched off, so is not included in the final TestHepMC efficiency ");
840  ATH_MSG_INFO(" Event rate with negative energy particles = " << m_negativeEnergyCheckRate*100.0/denom << "%");
841  if (!m_negativeEnergyTest) ATH_MSG_INFO(" The check for particles with negative energy is switched off, so is not included in the final TestHepMC efficiency ");
842  ATH_MSG_INFO(" Event rate with tachyons = " << m_tachyonCheckRate*100.0/denom << "%");
843  if (!m_tachyonsTest) ATH_MSG_INFO(" The check for tachyons is switched off, so is not included in the final TestHepMC efficiency ");
844  ATH_MSG_INFO(" Event rate with stable or unstable particles with no parents = " << m_stableUnstableNoParentCheckRate*100.0/denom << "%");
845  ATH_MSG_INFO(" Event rate with unstable particle with no decay vertex = " << m_unstablePartNoDecayVtxCheckRate*100.0/denom << "%");
846  if (!m_unstableNoVtxTest) ATH_MSG_INFO(" The check for unstable part. without end vertex is switched off, so is not included in the final TestHepMC efficiency ");
847  ATH_MSG_INFO(" Event rate with undecayed Pi0's = " << m_undecayedPi0CheckRate*100.0/denom << "%");
848  if (!m_pi0NoVtxTest) ATH_MSG_INFO(" The check for undecayed pi0's is switched off, so is not included in the final TestHepMC efficiency ");
849  ATH_MSG_INFO(" Event rate with undisplaced decay daughters of displaced vertices = " << m_undisplacedDecayDaughtersOfDisplacedVtxCheckRate*100.0/denom << "%");
850  if (!m_undisplacedDaughtersTest) ATH_MSG_INFO(" The check for undisplaced daughters is switched off, so is not included in the final TestHepMC efficiency ");
851  ATH_MSG_INFO(" Event rate with particles with status 1 but lifetime < " << m_min_tau << "~ns = " << m_Status1ShortLifetime*100.0/denom << "%");
852  if (!m_lifeTimeTest) ATH_MSG_INFO(" The check for status 1 particles with too short lifetime is switched off, so is not included in the final TestHepMC efficiency ");
853  ATH_MSG_INFO(" Event rate with energy sum of interacting particles non known by Geant4 above " << m_nonG4_energy_threshold << " MeV = " << m_nonG4_energyCheckRate*100.0/denom << "%");
854  if (!m_energyG4Test) ATH_MSG_INFO(" The check for energy not known by G4 is switched off, so is not included in the final TestHepMC efficiency ");
855  ATH_MSG_INFO(" Event rate with unknown PDG IDs = " << m_unknownPDGIDCheckRate*100.0/denom << "%");
856  if (!m_unknownPDGIDTest) ATH_MSG_INFO(" The check for unknown PDG IDs is sitched off, so it is not included in the final TestHepMC efficiency ");
857 
858  const double tau_fastDrate = double(m_FastDecayedTau) / double(m_TotalTaus);
859  if(tau_fastDrate > m_tau_eff_threshold){
860  ATH_MSG_FATAL("MORE THAN " << 100.*m_tau_eff_threshold << "% OF TAUS DECAYING IMMEDIATELY! " << m_FastDecayedTau << " found, out of: " << m_TotalTaus);
861  return StatusCode::FAILURE;
862  }
863 
864  if (m_noXSECset) {
865  if (m_allowMissingXSec) {
866  ATH_MSG_WARNING(m_noXSECset << " EVENTS WITHOUT CROSS-SECTION!!! Added dummy cross-section instead.");
867  }
868  else {
869  ATH_MSG_FATAL(m_noXSECset << " EVENTS WITHOUT CROSS-SECTION!! Check the setup before production!");
870  return StatusCode::FAILURE;
871  }
872  }
873 
874  const double efficiency = double(m_nPass) / double(m_nPass + m_nFail);
875  ATH_MSG_INFO("Efficiency = " << efficiency * 100 << "%");
876 
877  // Check efficiency, and fail (to kill production jobs) if the pass rate is too low
878  if (efficiency < m_eff_fail_threshold) {
879  ATH_MSG_FATAL("EFFICIENCY ABOVE ERROR THRESHOLD! " << 100*efficiency << "% found, but at least: " << 100*m_eff_fail_threshold << "% required");
880  return StatusCode::FAILURE;
881  } else if (efficiency <= m_eff_warn_threshold) {
882  ATH_MSG_WARNING("EFFICIENCY ABOVE WARNING THRESHOLD! " << 100*efficiency << "% found, but at least: " << 100*m_eff_warn_threshold << "% expected");
883  }
884 
885  return StatusCode::SUCCESS;
886 }
887 
888 
889 #endif
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
TestHepMC::m_cm_energy
double m_cm_energy
Definition: TestHepMC.h:45
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
TestHepMC::m_h_status_dispVtxCheck
TH1F * m_h_status_dispVtxCheck
Definition: TestHepMC.h:104
TestHepMC::m_tachyonCheckRate
int m_tachyonCheckRate
Definition: TestHepMC.h:82
TestHepMC::m_h_vx_dispVtxCheck
TH1F * m_h_vx_dispVtxCheck
Definition: TestHepMC.h:108
TestHepMC::m_h_py_dispVtxCheck
TH1F * m_h_py_dispVtxCheck
Definition: TestHepMC.h:106
TestHepMC::m_doHist
bool m_doHist
Definition: TestHepMC.h:51
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:18
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
TestHepMC::m_unstableNoVtxTest
bool m_unstableNoVtxTest
Definition: TestHepMC.h:53
isNucleus
bool isNucleus(const T &p)
PDG rule 16 Nuclear codes are given as 10-digit numbers ±10LZZZAAAI.
Definition: AtlasPID.h:697
HepMC::Print::content
void content(std::ostream &os, const GenEvent &e)
Definition: GenEvent.h:678
TestHepMC::m_TotalTaus
int m_TotalTaus
Definition: TestHepMC.h:62
TestHepMC::m_energyBalanceCheckRate
int m_energyBalanceCheckRate
Definition: TestHepMC.h:79
TestHepMC::m_unstablePartNoDecayVtxCheckRate
int m_unstablePartNoDecayVtxCheckRate
Definition: TestHepMC.h:84
PowhegControl_ttHplus_NLO.ss
ss
Definition: PowhegControl_ttHplus_NLO.py:83
TestHepMC::m_eff_fail_threshold
double m_eff_fail_threshold
Definition: TestHepMC.h:49
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
TestHepMC::m_beamParticleswithStatusNotFourCheckRate
int m_beamParticleswithStatusNotFourCheckRate
Definition: TestHepMC.h:66
TestHepMC::m_negativeEnergyTest
bool m_negativeEnergyTest
Definition: TestHepMC.h:53
baryonNumber
double baryonNumber(const T &p)
Definition: AtlasPID.h:768
GenBase::events_const
const McEventCollection * events_const() const
Access the current event's McEventCollection (const)
Definition: GenBase.h:96
TestHepMC::m_noXSECset
int m_noXSECset
Definition: TestHepMC.h:60
TestHepMC.h
TestHepMC::m_negativeEnergyCheckRate
int m_negativeEnergyCheckRate
Definition: TestHepMC.h:81
TestHepMC::m_tachyonsTest
bool m_tachyonsTest
Definition: TestHepMC.h:53
TestHepMC::m_maxloops
int m_maxloops
Definition: TestHepMC.h:43
TestHepMC::m_h_vtxprod_dispVtxCheck
TH1F * m_h_vtxprod_dispVtxCheck
Definition: TestHepMC.h:115
TestHepMC::m_FastDecayedTau
int m_FastDecayedTau
Definition: TestHepMC.h:63
TestHepMC::m_thistSvc
ServiceHandle< ITHistSvc > m_thistSvc
Definition: TestHepMC.h:99
TestHepMC::m_h_momentumImbalance_pz
TH1F * m_h_momentumImbalance_pz
Definition: TestHepMC.h:124
TestHepMC::m_vtxDisplacedstatuscodenot12CheckRate
int m_vtxDisplacedstatuscodenot12CheckRate
Definition: TestHepMC.h:70
TestHepMC::m_vtxDisplacedstatuscode12CheckRate
int m_vtxDisplacedstatuscode12CheckRate
Definition: TestHepMC.h:69
TestHepMC::m_unknownPDGIDFile
std::string m_unknownPDGIDFile
Definition: TestHepMC.h:92
makeDTCalibBlob_pickPhase.pd
pd
Definition: makeDTCalibBlob_pickPhase.py:342
LArG4FSStartPointFilter.evt
evt
Definition: LArG4FSStartPointFilter.py:42
TestHepMC::m_energy_diff
double m_energy_diff
Definition: TestHepMC.h:46
isValid
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
Definition: AtlasPID.h:872
TestHepMC::m_unknownPDGIDCheckRate
int m_unknownPDGIDCheckRate
Definition: TestHepMC.h:89
TestHepMC::m_momentumBalanceCheckRate
int m_momentumBalanceCheckRate
Definition: TestHepMC.h:80
TestHepMC::m_pdg
int m_pdg
Definition: TestHepMC.h:44
HepMC::Print::line
void line(std::ostream &os, const GenEvent &e)
Definition: GenEvent.h:676
TestHepMC::m_G4pdgID_tab
std::vector< int > m_G4pdgID_tab
Definition: TestHepMC.h:94
TestHepMC::m_h_px_dispVtxCheck
TH1F * m_h_px_dispVtxCheck
Definition: TestHepMC.h:105
TestHepMC::m_uknownPDGID_tab
std::vector< int > m_uknownPDGID_tab
Definition: TestHepMC.h:96
dq_defect_bulk_create_defects.line
line
Definition: dq_defect_bulk_create_defects.py:27
AthenaPoolTestWrite.stream
string stream
Definition: AthenaPoolTestWrite.py:12
XMLtoHeader.count
count
Definition: XMLtoHeader.py:84
TestHepMC::m_eff_warn_threshold
double m_eff_warn_threshold
Definition: TestHepMC.h:49
TestHepMC::m_accur_margin
double m_accur_margin
Definition: TestHepMC.h:50
TestHepMC::m_dumpEvent
bool m_dumpEvent
Definition: TestHepMC.h:47
TestHepMC::m_nonG4_energy_threshold
double m_nonG4_energy_threshold
Definition: TestHepMC.h:48
TestHepMC::m_h_pz_dispVtxCheck
TH1F * m_h_pz_dispVtxCheck
Definition: TestHepMC.h:107
TestHepMC::m_momImbalanceTest
bool m_momImbalanceTest
Definition: TestHepMC.h:53
TestHepMC::execute
StatusCode execute()
Definition: TestHepMC.cxx:269
TestHepMC::m_max_dist
double m_max_dist
Definition: TestHepMC.h:48
TestHepMC::m_allowMissingXSec
bool m_allowMissingXSec
Definition: TestHepMC.h:47
TestHepMC::m_stableUnstableNoParentCheckRate
int m_stableUnstableNoParentCheckRate
Definition: TestHepMC.h:83
TestHepMC::m_unstableNoEndVtxCheckRate
int m_unstableNoEndVtxCheckRate
Definition: TestHepMC.h:74
TestHepMC::m_SusyPdgID_tab
std::vector< int > m_SusyPdgID_tab
Definition: TestHepMC.h:95
efficiency
void efficiency(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:128
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
CaloCondBlobAlgs_fillNoiseFromASCII.desc
desc
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:53
GenBase
Base class for common behaviour of MC truth algorithms.
Definition: GenBase.h:47
TestHepMC::m_vertexStatuses
std::vector< int > m_vertexStatuses
Definition: TestHepMC.h:56
TestHepMC::m_undecayedPi0statuscode12CheckRate
int m_undecayedPi0statuscode12CheckRate
Definition: TestHepMC.h:73
TestHepMC::m_vtxNANandINFCheckRate
int m_vtxNANandINFCheckRate
Definition: TestHepMC.h:68
TestHepMC::m_h_photon_energy
TH1F * m_h_photon_energy
Definition: TestHepMC.h:118
lumiFormat.i
int i
Definition: lumiFormat.py:85
TestHepMC::finalize
StatusCode finalize()
Definition: TestHepMC.cxx:813
TestHepMC::m_partMomentumNANandINFCheckRate
int m_partMomentumNANandINFCheckRate
Definition: TestHepMC.h:72
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
TestHepMC::m_h_beamparticle2_Energy
TH1F * m_h_beamparticle2_Energy
Definition: TestHepMC.h:127
TestHepMC::m_h_vyprod_dispVtxCheck
TH1F * m_h_vyprod_dispVtxCheck
Definition: TestHepMC.h:112
TestHepMC::m_undecayedPi0CheckRate
int m_undecayedPi0CheckRate
Definition: TestHepMC.h:85
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition: AthCommonDataStore.h:145
TestHepMC::m_h_momentumImbalance_py
TH1F * m_h_momentumImbalance_py
Definition: TestHepMC.h:123
TestHepMC::m_vtxNaNTest
bool m_vtxNaNTest
Definition: TestHepMC.h:52
TestHepMC::m_nFail
int m_nFail
Definition: TestHepMC.h:59
TestHepMC::m_vtxDisplacedMoreThan_1m_CheckRate
int m_vtxDisplacedMoreThan_1m_CheckRate
Definition: TestHepMC.h:71
find_tgc_unfilled_channelids.ip
ip
Definition: find_tgc_unfilled_channelids.py:3
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
TestHepMC::m_h_vtxend_dispVtxCheck
TH1F * m_h_vtxend_dispVtxCheck
Definition: TestHepMC.h:114
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
TestHepMC::m_momNaNTest
bool m_momNaNTest
Definition: TestHepMC.h:52
TestHepMC::m_pi0NoVtxTest
bool m_pi0NoVtxTest
Definition: TestHepMC.h:54
TestHepMC::m_energyG4Test
bool m_energyG4Test
Definition: TestHepMC.h:52
TestHepMC::m_paramFile
std::string m_paramFile
Definition: TestHepMC.h:91
TestHepMC::m_h_pdgid_dispVtxCheck
TH1F * m_h_pdgid_dispVtxCheck
Definition: TestHepMC.h:103
compute_lumi.denom
denom
Definition: compute_lumi.py:76
TestHepMC::m_h_photon_e2_p2_e2
TH1F * m_h_photon_e2_p2_e2
Definition: TestHepMC.h:119
TestHepMC::m_beamEnergyTest
bool m_beamEnergyTest
Definition: TestHepMC.h:52
TestHepMC::m_h_vzprod_dispVtxCheck
TH1F * m_h_vzprod_dispVtxCheck
Definition: TestHepMC.h:113
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition: GenParticle.h:38
PathResolver.h
TestHepMC::m_h_vy_dispVtxCheck
TH1F * m_h_vy_dispVtxCheck
Definition: TestHepMC.h:109
MC::Loops::loop_particles
const std::vector< Prt > & loop_particles() const
Definition: Loops.h:25
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
TestHepMC::m_beamEnergyCheckRate
int m_beamEnergyCheckRate
Definition: TestHepMC.h:67
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
TestHepMC::m_energyImbalanceTest
bool m_energyImbalanceTest
Definition: TestHepMC.h:53
TestHepMC::m_nonG4_energyCheckRate
int m_nonG4_energyCheckRate
Definition: TestHepMC.h:88
TestHepMC::m_nonZeroPhotonMassCheckRate
int m_nonZeroPhotonMassCheckRate
Definition: TestHepMC.h:78
TestHepMC::m_min_tau
double m_min_tau
Definition: TestHepMC.h:48
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
TestHepMC::m_looper
MC::Loops< HepMC::GenEvent, HepMC::ConstGenParticlePtr, HepMC::ConstGenVertexPtr > m_looper
member to detect loops
Definition: TestHepMC.h:130
TestHepMC::m_invalidBeamParticlesCheckRate
int m_invalidBeamParticlesCheckRate
Definition: TestHepMC.h:65
PathResolverFindDataFile
std::string PathResolverFindDataFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:278
MC::Loops::loop_vertices
const std::vector< Vtx > & loop_vertices() const
Definition: Loops.h:28
TestHepMC::m_h_vxprod_dispVtxCheck
TH1F * m_h_vxprod_dispVtxCheck
Definition: TestHepMC.h:111
TestHepMC::TestHepMC
TestHepMC(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TestHepMC.cxx:17
MC::isStable
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
Definition: HepMCHelpers.h:45
TestHepMC::m_nPass
int m_nPass
Definition: TestHepMC.h:58
TestHepMC::m_tau_eff_threshold
double m_tau_eff_threshold
Definition: TestHepMC.h:49
TestHepMC::m_cme_diff
double m_cme_diff
Definition: TestHepMC.h:45
TestHepMC::m_unknownPDGIDTest
bool m_unknownPDGIDTest
Definition: TestHepMC.h:54
TestHepMC::m_h_energy_dispVtxCheck
TH1F * m_h_energy_dispVtxCheck
Definition: TestHepMC.h:101
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
HepMC::valid_beam_particles
bool valid_beam_particles(const GenEvent *e)
Definition: GenEvent.h:681
numberOfProtons
int numberOfProtons(const T &p)
Definition: AtlasPID.h:828
MC::isDecayed
bool isDecayed(const T &p)
Identify if the particle decayed.
Definition: HepMCHelpers.h:42
xAOD::EgammaHelpers::isPhoton
bool isPhoton(const xAOD::Egamma *eg)
is the object a photon
Definition: EgammaxAODHelpers.cxx:21
GetAllXsec.xsec
xsec
Definition: GetAllXsec.py:85
MC::isBeam
bool isBeam(const T &p)
Identify if the particle is beam particle.
Definition: HepMCHelpers.h:39
TestHepMC::m_max_energy_diff
double m_max_energy_diff
Definition: TestHepMC.h:46
AthCommonMsg< Algorithm >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
MC::isSimInteracting
bool isSimInteracting(const T &p)
Identify if the particle could interact with the detector during the simulation, e....
Definition: HepMCHelpers.h:60
TestHepMC::initialize
StatusCode initialize()
Definition: TestHepMC.cxx:130
TestHepMC::m_h_vz_dispVtxCheck
TH1F * m_h_vz_dispVtxCheck
Definition: TestHepMC.h:110
TestHepMC::m_Status1ShortLifetime
int m_Status1ShortLifetime
Definition: TestHepMC.h:86
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
MC::Loops::findLoops
int findLoops(const Evt *evt, bool force)
Definition: Loops.h:31
TestHepMC::m_undisplacedDecayDaughtersOfDisplacedVtxCheckRate
int m_undisplacedDecayDaughtersOfDisplacedVtxCheckRate
Definition: TestHepMC.h:87
TestHepMC::m_max_dist_trans
double m_max_dist_trans
Definition: TestHepMC.h:48
TestHepMC::m_h_cmEnergyDiff
TH1F * m_h_cmEnergyDiff
Definition: TestHepMC.h:128
TestHepMC::m_undisplacedLLHdaughtersCheckRate
int m_undisplacedLLHdaughtersCheckRate
Definition: TestHepMC.h:77
TestHepMC::m_h_energyImbalance
TH1F * m_h_energyImbalance
Definition: TestHepMC.h:121
TestHepMC::m_h_beamparticle1_Energy
TH1F * m_h_beamparticle1_Energy
Definition: TestHepMC.h:126
TestHepMC::m_vtxDisplacedTest
bool m_vtxDisplacedTest
Definition: TestHepMC.h:52
TestHepMC::m_negativeEnergyTachyonicCheckRate
int m_negativeEnergyTachyonicCheckRate
Definition: TestHepMC.h:75
GenBase::particleData
const HepPDT::ParticleData * particleData(int pid) const
Access an element in the particle data table.
Definition: GenBase.h:126
TestHepMC::m_h_energy_dispVtxCheck_lt10
TH1F * m_h_energy_dispVtxCheck_lt10
Definition: TestHepMC.h:102
GenBase::initialize
virtual StatusCode initialize() override
Definition: GenBase.cxx:17
TestHepMC::m_decayCheckRate
int m_decayCheckRate
Definition: TestHepMC.h:76
TestHepMC::m_undisplacedDaughtersTest
bool m_undisplacedDaughtersTest
Definition: TestHepMC.h:54
HepMCHelpers.h
TestHepMC::m_lifeTimeTest
bool m_lifeTimeTest
Definition: TestHepMC.h:52
TestHepMC::m_h_photon_mass
TH1F * m_h_photon_mass
Definition: TestHepMC.h:117
TestHepMC::m_h_momentumImbalance_px
TH1F * m_h_momentumImbalance_px
Definition: TestHepMC.h:122
Generated on Sun Sep 28 2025 21:20:10 for ATLAS Offline Software by   doxygen 1.8.18