Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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=true); //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(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(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 
361  if (m_cm_energy > 0 && std::abs(cmenergy - m_cm_energy) > m_cme_diff) {
362  ATH_MSG_FATAL("Beam particles have incorrect energy: " << m_cm_energy/Gaudi::Units::GeV << " GeV expected, vs. " << cmenergy/Gaudi::Units::GeV << " GeV found");
363  setFilterPassed(false);
364  if (m_doHist){
365  m_h_beamparticle1_Energy->Fill(beams.first->momentum().e()/Gaudi::Units::GeV);
366  m_h_beamparticle2_Energy->Fill(beams.second->momentum().e()/Gaudi::Units::GeV);
367  m_h_cmEnergyDiff->Fill((cmenergy-m_cm_energy)/Gaudi::Units::GeV);
368  }
369  ++m_beamEnergyCheckRate;
370  // Special case: this is so bad that we immediately fail out
371  return StatusCode::FAILURE;
372  }
373  }
374 
375  // Check vertices
376  int vtxDisplacedstatuscode12CheckRateCnt=0;
377  int vtxDisplacedstatuscodenot12CheckRateCnt=0;
378  int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
379 #ifdef HEPMC3
380  for (const auto& vtx: evt->vertices()) {
381 #else
382  for (auto vitr = evt->vertices_begin(); vitr != evt->vertices_end(); ++vitr ) {
383  const HepMC::GenVertex* vtx = *vitr;
384 #endif
385  const HepMC::FourVector pos = vtx->position();
386 
387  // Check for NaNs and infs in vertex position components
388  if ( std::isnan(pos.x()) || std::isinf(pos.x()) ||
389  std::isnan(pos.y()) || std::isinf(pos.y()) ||
390  std::isnan(pos.z()) || std::isinf(pos.z()) ) {
391  ATH_MSG_WARNING("NaN (Not A Number) or inf found in the event record vertex positions");
392 
393  ++m_vtxNANandINFCheckRate;
394  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
395  if (m_vtxNaNTest) {
396  filter_pass = false;
397  }
398  } // Done of checking for nans and infinities
399 
400  // Check for too-far-displaced vertices
401  // Anything which propagates macroscopically should be set stable in evgen for G4 to handle
402  const double dist_trans2 = pos.x()*pos.x() + pos.y()*pos.y(); // in mm2
403  const double dist2 = dist_trans2 + pos.z()*pos.z(); // in mm2
404  const double dist_trans = std::sqrt(dist_trans2); // in mm
405  const double dist = std::sqrt(dist2); // in mm
406  if (dist2 > m_max_dist*m_max_dist) {
407  ATH_MSG_WARNING("Found vertex position displaced by more than " << m_max_dist << "mm: " << dist << "mm");
408  ++vtxDisplacedMoreThan_1m_CheckRateCnt;
409 
410  if (m_vtxDisplacedTest) {
411  filter_pass = false;
412  }
413  }
414  if (dist_trans2 > m_max_dist_trans*m_max_dist_trans) {
415  ATH_MSG_WARNING("Found vertex position displaced by more than " << m_max_dist_trans << "mm in transverse distance: " << dist_trans << "mm");
416 
417 #ifdef HEPMC3
418  for (const auto& part: vtx->particles_in()) {
419 #else
420  for (auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
421  auto part=(*part_it);
422 #endif
423  if (m_dumpEvent){
424  ATH_MSG_WARNING("Outgoing particle : ");
425  HepMC::Print::line(msg( MSG::WARNING ).stream(),part);
426  }
427  ATH_MSG_WARNING("production vertex = " << part->production_vertex()->position().x() << ", " << part->production_vertex()->position().y() << ", " << part->production_vertex()->position().z());
428  ATH_MSG_WARNING("end vertex = " << part->end_vertex()->position().x() << ", " << part->end_vertex()->position().y() << ", " << part->end_vertex()->position().z());
429  if (m_dumpEvent) ATH_MSG_WARNING("parents info: ");
430  if (part->production_vertex()) {
431 #ifdef HEPMC3
432  for(const auto& p_parents: part->production_vertex()->particles_in()) {
433 #else
434  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) {
435  auto p_parents=(*p_parents_it);
436 #endif
437  if (m_dumpEvent){
438  msg(MSG::WARNING) << "\t";
439  HepMC::Print::line( msg( MSG::WARNING ).stream() , p_parents );
440  }
441  }
442  } // Done with fancy print
443 
444  if (part->status() == 1 || part->status() == 2){
445  vtxDisplacedstatuscode12CheckRateCnt += 1;
446  } else {
447  vtxDisplacedstatuscodenot12CheckRateCnt += 1;
448  }
449 
450  if (m_doHist){
451  m_h_energy_dispVtxCheck->Fill(part->momentum().e()/Gaudi::Units::GeV);
452  if (part->momentum().e()/Gaudi::Units::GeV < 10.) {
453  m_h_energy_dispVtxCheck_lt10->Fill(part->momentum().e()/Gaudi::Units::GeV);
454  }
455  m_h_pdgid_dispVtxCheck->Fill(part->pdg_id());
456  m_h_status_dispVtxCheck->Fill(part->status());
457  m_h_px_dispVtxCheck->Fill(part->momentum().px()/Gaudi::Units::GeV);
458  m_h_py_dispVtxCheck->Fill(part->momentum().py()/Gaudi::Units::GeV);
459  m_h_pz_dispVtxCheck->Fill(part->momentum().pz()/Gaudi::Units::GeV);
460  m_h_vx_dispVtxCheck->Fill(part->end_vertex()->position().x());
461  m_h_vy_dispVtxCheck->Fill(part->end_vertex()->position().y());
462  m_h_vz_dispVtxCheck->Fill(part->end_vertex()->position().z());
463  m_h_vxprod_dispVtxCheck->Fill(part->production_vertex()->position().x());
464  m_h_vyprod_dispVtxCheck->Fill(part->production_vertex()->position().y());
465  m_h_vzprod_dispVtxCheck->Fill(part->production_vertex()->position().z());
466  double endvx = part->end_vertex()->position().x();
467  double endvy = part->end_vertex()->position().y();
468  double endvz = part->end_vertex()->position().z();
469  double prodvx = part->production_vertex()->position().x();
470  double prodvy = part->production_vertex()->position().y();
471  double prodvz = part->production_vertex()->position().z();
472  double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
473  double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
474  m_h_vtxend_dispVtxCheck->Fill(enddis);
475  m_h_vtxprod_dispVtxCheck->Fill(proddis);
476  } // End of the filling of histograms for bad vertices
477  } // End of a loop over theparents of the bad vertex
478  } // Found a bad vertex
479  } // Loop over all vertices
480  if (vtxDisplacedstatuscode12CheckRateCnt>0) ++m_vtxDisplacedstatuscode12CheckRate;
481  if (vtxDisplacedstatuscodenot12CheckRateCnt>0) ++m_vtxDisplacedstatuscodenot12CheckRate;
482  if (vtxDisplacedMoreThan_1m_CheckRateCnt>0) ++m_vtxDisplacedMoreThan_1m_CheckRate;
483 
484  // Check particles
485  for (auto pitr: *evt) {
486 
487  // Local loop variables to clean up the check code
488  const HepMC::FourVector pmom = pitr->momentum();
489  const int pstatus = pitr->status();
490  const int ppdgid = pitr->pdg_id();
491  // Check for NaNs and infs in momentum components
492  if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
493  std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
494  std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
495  std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
496  ATH_MSG_WARNING("NaN (Not A Number) or inf found in the event record momenta");
497  ++m_partMomentumNANandINFCheckRate;
498 
499  if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
500  if (m_momNaNTest) {
501  filter_pass = false;
502  }
503  } // End of check for NaNs and infinities
504 
505  // Check for undecayed pi0s
506  if (MC::isStable(pstatus) || MC::isDecayed(pstatus)) {
507  if (ppdgid == 111 && !pitr->end_vertex() ) {
508  unDecPi0.push_back( pitr);
509  ++m_undecayedPi0CheckRate;
510  }
511  } // End of check for undecayed pi0s
512 
513  //check stable particle lifetimes
514  if (MC::isStable(pstatus)) {
515  const HepPDT::ParticleData* pd = particleData(ppdgid);
516  if (pd != NULL) {
517  double plifetime = pd->lifetime()*1e+12; // why lifetime doesn't come in common units???
518  if (plifetime != 0 && plifetime < m_min_tau) { // particles with infinite lifetime get a 0 in the PDT
519  ATH_MSG_WARNING("Stable particle found with lifetime = " << plifetime << "~ns!!");
520  if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
521 
522  ++m_Status1ShortLifetime;
523 
524  if (m_lifeTimeTest) {
525  filter_pass = false;
526  }
527  } // Particle did not have infinite lifetime
528  } // The particle has a data table (so a lifetime)
529  else{
530  int susyPart = 0;
531  std::vector<int>::size_type count = 0;
532  while (susyPart==0 && (count < m_SusyPdgID_tab.size() )){
533  // no warning for SUSY particles from the list susyParticlePdgid.txt
534  if (m_SusyPdgID_tab[count] == std::abs(ppdgid)) {
535  susyPart=1;
536  }
537  count++;
538  } // Look through the SUSY table to see if this one should be counted
539  if (susyPart==0){
540  ATH_MSG_WARNING("Stable particle not found in PDT, no lifetime check done");
541  if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
542  } // It's a SUSY particle -- skip the lifetime check
543  } // The particle has no data table
544  } // Test if the particle is stable
545 
546  //Check that stable particles are known by G4 or they are non-interacting
547  const MC::DecodedPID decodedPID(ppdgid);
548  const int first_dig = decodedPID(0);
549 
550  if (MC::isStable(pstatus) && (!pitr->end_vertex()) && (MC::isSimInteracting(pitr)) && (!MC::isNucleus(ppdgid)) && (first_dig != 9) ) {
551 
552  int known_byG4 = 0;
553  std::vector<int>::size_type count =0;
554 
555  while (known_byG4==0 && count < m_G4pdgID_tab.size()){
556  if(ppdgid == m_G4pdgID_tab[count]) known_byG4=1;
557  count++;
558  }
559  if(known_byG4==0){
560  nonG4_energy += pmom.e();
561  ATH_MSG_WARNING("Interacting particle not known by Geant4 with ID " << ppdgid);
562  }
563  } // End of check that stable particles are known to G4 or are non-interacting
564 
565  // Check for bad PDG IDs
566  if (!MC::isValid(ppdgid)){
567  ATH_MSG_DEBUG("Invalid PDG ID found: " << ppdgid);
568  if (m_unknownPDGIDTest && std::find(m_uknownPDGID_tab.begin(),m_uknownPDGID_tab.end(),ppdgid)==m_uknownPDGID_tab.end()){
569  ATH_MSG_WARNING("Invalid and unmasked PDG ID found: " << ppdgid);
570  filter_pass = false;
571  }
572  } // End of check for invalid PDG IDs
573 
574  // Check for unstables with no end vertex,
575  if (!pitr->end_vertex() && MC::isDecayed(pstatus)) {
576  unstNoEnd.push_back(pitr);
577  ++m_unstableNoEndVtxCheckRate;
578  } // End of check for unstable with no end vertex
579 
580  // Sum final state mom/energy, and note negative energy / tachyonic particles
581  // std::cout << "status " << pstatus << " e " << pmom.e() << " pz " << pmom.pz()<< std::endl;
582  if ( MC::isStable(pstatus) && !pitr->end_vertex() ) {
583  totalPx += pmom.px();
584  totalPy += pmom.py();
585  totalPz += pmom.pz();
586  totalE += pmom.e();
587  if (pmom.e() < 0) {
588  negEnPart.push_back(pitr);
589  ++m_negativeEnergyTachyonicCheckRate;
590  }
591  const double aener = std::abs(pmom.e());
592  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()) ) {
593  tachyons.push_back(pitr);
594  ++m_negativeEnergyTachyonicCheckRate;
595  }
596  } // End of sums for momentum and energy conservation
597 
598  // Decay checks (uses PdgToSearch attr value, for tau by default)
600  int tau_child = 0;
601  if (std::abs(ppdgid) == m_pdg && (MC::isStable(pstatus) || MC::isDecayed(pstatus))) {
602  ++m_TotalTaus;
603  auto vtx = pitr->end_vertex();
604  if (vtx) {
605  double p_energy = 0;
606 #ifdef HEPMC3
607  for (auto desc: HepMC::descendant_particles(vtx)) {
608 #else
609  for (auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
610  auto desc=(*desc_it);
611 #endif
612  if (std::abs(desc->pdg_id()) == m_pdg) tau_child = 1;
613  if ( MC::isStable(desc) ) p_energy += desc->momentum().e();
614  }
615  if (std::abs( p_energy - pmom.e()) > m_energy_diff && !tau_child) {
616  ATH_MSG_WARNING("Energy sum (decay products): "
617  << "Energy (original particle) > " << m_energy_diff << " MeV, "
618  << "Event #" << evt->event_number() << ", "
619  << "The original particle = " << pitr);
620  ++m_decayCheckRate;
621  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
622  }
623  //most taus should not decay immediately
624  const HepMC::FourVector tau_decaypos = vtx->position();
625  const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
626  //tau_child != 1 exclude cases in which a tau is copied to another vertex or emits a photon
627  if ((tau_displacement < 1.e-6) && (tau_child!=1)) ++m_FastDecayedTau;
628  } else {
629  ATH_MSG_WARNING("UNDECAYED PARTICLE WITH PDG_ID = " << m_pdg);
630  ++m_decayCheckRate;
631  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
632  }
633  } // End of checks for specific particle (tau by default)
634 
635  // Check for undisplaced decay daughters from long-lived hadrons
636  if (pitr->end_vertex()) {
637  auto decayvtx = pitr->end_vertex();
638  const HepMC::FourVector decaypos = decayvtx->position();
639  const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
640  if (displacement > 1e-6) {
641  for (auto ip: *decayvtx) {
642  const HepMC::FourVector pos2 = ip->production_vertex()->position();
643  const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
644  if (displacement2 < 1e-6) {
645  ATH_MSG_WARNING("Decay child " << ip << " from " << pitr
646  << " has undisplaced vertex (" << ip->production_vertex()
647  << " @ " << displacement2 << "mm) "
648  << " but parent vertex is displaced (" << decayvtx
649  << " @ " << displacement << "mm)");
650  undisplaceds.push_back(ip);
651  ++m_undisplacedLLHdaughtersCheckRate;
652  } // Check for displacement below 1 um
653  } // Loop over all particles coming from the decay vertex
654  } // Displacement of greater than 1 um
655  } // End of check for undisplaced decay daughters from long-lived hadrons
656 
657  // Check for photons with non-zero masses
659  if (MC::isPhoton(ppdgid) && MC::isStable(pstatus)) {
660  const double mass = pitr->generated_mass();
661  if (std::abs(mass) > 1.0) { // in MeV
662  ATH_MSG_WARNING("Photon with non-zero mass found! Mass: " << mass << " MeV" << pitr);
663  ++m_nonZeroPhotonMassCheckRate;
664  }
665  } // End check for photons with too-large a mass
666 
667  } // End of loop over particles in the event
668 
669  // Energy of interacting particles not known by Geant4
670  if(nonG4_energy > m_nonG4_energy_threshold) {
671  ATH_MSG_WARNING("The energy of interacting particles not known by Geant4 is = " << nonG4_energy << " MeV");
672  if (m_energyG4Test) {
673  filter_pass = false;
674  }
675  ++m_nonG4_energyCheckRate;
676  } // End of check for interacting particles not known by G4
677 
678  // Energy balance
679  double lostE = std::abs(totalE - cmenergy);
680  if (lostE > m_energy_diff) {
681  ATH_MSG_WARNING("ENERGY BALANCE FAILED : E-difference = " << lostE << " MeV");
682 
683  ATH_MSG_WARNING("balance " << totalPx << " " << totalPy << " " << totalPz << " " << totalE);
684 
685  if (m_doHist){
686  m_h_energyImbalance->Fill(lostE/Gaudi::Units::GeV);
687  }
688  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
689  if (m_energyImbalanceTest) {
690  filter_pass = false;
691  }
692  ++m_energyBalanceCheckRate;
693  } // End of energy balance check
694 
695  // Momentum balance
696  if ( std::abs(totalPx) > m_energy_diff || std::abs(totalPy) > m_energy_diff || std::abs(totalPz) > m_energy_diff ) {
697  ATH_MSG_WARNING("MOMENTUM BALANCE FAILED : SumPx = " << totalPx << " SumPy = " << totalPy << " SumPz = " << totalPz << " MeV");
698  if (m_doHist){
699  m_h_momentumImbalance_px->Fill(std::abs(totalPx)/Gaudi::Units::GeV);
700  m_h_momentumImbalance_py->Fill(std::abs(totalPy)/Gaudi::Units::GeV);
701  m_h_momentumImbalance_pz->Fill(std::abs(totalPz)/Gaudi::Units::GeV);
702  }
703  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
704  if (m_momImbalanceTest) {
705  filter_pass = false;
706  }
707  ++m_momentumBalanceCheckRate;
708  } // End of momentum balance check
709 
710  // Negative energy particles
711  if (!negEnPart.empty()) {
712  std::stringstream ss;
713  ss << "NEGATIVE ENERGY PARTICLES FOUND :";
714  for (auto b: negEnPart){
715  ss << " " << b;
716  }
717  ATH_MSG_WARNING(ss.str());
718  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
719  if (m_negativeEnergyTest) {
720  filter_pass = false;
721  }
722  ++m_negativeEnergyCheckRate;
723  } // End of negative energy particle chedk
724 
725  // Tachyons
726  if (!tachyons.empty()) {
727  std::stringstream ss;
728  ss << "PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
729  for (auto b: tachyons){
730  ss << " " << b;
731  }
732  ATH_MSG_WARNING(ss.str());
733  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
734  if (m_tachyonsTest) {
735  filter_pass = false;
736  }
737  ++m_energyBalanceCheckRate;
738  } // End of tachyon check
739 
740  // Unstable particles with no decay vertex
741  if (!unstNoEnd.empty()) {
742  std::stringstream ss;
743  ss << "Unstable particle with no decay vertex found: ";
744  for (auto b: unstNoEnd){
745  ss << " " << b;
746  }
747  ATH_MSG_WARNING(ss.str());
748  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
749  if (m_unstableNoVtxTest) {
750  filter_pass = false;
751  }
752  ++m_unstablePartNoDecayVtxCheckRate;
753  } // End of unstable particle with no decay vertex check
754 
755  // Undecayed pi0
756  if (!unDecPi0.empty()) {
757  std::stringstream ss;
758  ss << "pi0 with no decay vertex found:";
759  for (auto b: unDecPi0){
760  ss << " " << b;
761  }
762  ATH_MSG_WARNING(ss.str());
763  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
764  if (m_pi0NoVtxTest) {
765  filter_pass = false;
766  }
767  ++m_undecayedPi0CheckRate;
768  } // End of undecayed pi0 check
769 
770  // Undisplaced decay daughters of displaced vertices
771  if (!undisplaceds.empty()) {
772  std::stringstream ss;
773  ss << "Undisplaced decay vertices from displaced particle: ";
774  for (auto b: undisplaceds){
775  ss << " " << b;
776  }
777  ATH_MSG_WARNING(ss.str());
778  if (m_dumpEvent) HepMC::Print::content(std::cout,*evt);
779  if (m_undisplacedDaughtersTest) {
780  filter_pass = false;
781  }
782  ++m_undisplacedDecayDaughtersOfDisplacedVtxCheckRate;
783  } // End of undisplaced decay daughter of displaced vertices check
784 
785  } // End of loop over MCEventCollection
786 
787  // End of execution for each event - update filter value
788  if (!filter_pass){
789  setFilterPassed(false);
790  ++m_nFail;
791  } else {
792  ++m_nPass;
793  }
794 
795  // If the efficiency after 100 events is below 10%, there is an important bug going on:
796  // we fail the job immediately so it doesn't run for ever
797  const double tmp_efficiency = double(m_nPass) / double(m_nPass + m_nFail);
798  if ((m_nPass + m_nFail) > 100 && tmp_efficiency < 0.1) {
799  ATH_MSG_FATAL("The efficiency after " << m_nPass + m_nFail << " events is " << tmp_efficiency*100. << "% !!!");
800  return StatusCode::FAILURE;
801  }
802 
803  return StatusCode::SUCCESS;
804 }
805 
806 
807 StatusCode TestHepMC::finalize() {
808 
809  ATH_MSG_INFO("Events passed = " << m_nPass << ", Events Failed = " << m_nFail);
810 
811  ATH_MSG_INFO(" Event rate with invalid Beam Particles = " << m_invalidBeamParticlesCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
812  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)");
813  ATH_MSG_INFO(" Event rate with incorrect beam particle energies = " << m_beamEnergyCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
814  ATH_MSG_INFO(" Event rate with NaN (Not A Number) or inf found in the event record vertex positions = " << m_vtxNANandINFCheckRate*100.0/double(m_nPass + m_nFail) << "%");
815  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 ");
816  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/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
817  ATH_MSG_INFO(" Event rate with vertices displaced more than " << m_max_dist << "~mm = " << m_vtxDisplacedMoreThan_1m_CheckRate*100.0/double(m_nPass + m_nFail) << "%");
818  if (!m_vtxDisplacedTest) ATH_MSG_INFO(" The check for displaced vertices is switched off, so is not included in the final TestHepMC efficiency ");
819  ATH_MSG_INFO(" Event rate with NAN (Not A Number) or inf found in particle momentum values = " << m_partMomentumNANandINFCheckRate*100.0/double(m_nPass + m_nFail) << "%");
820  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 ");
821  ATH_MSG_INFO(" Event rate with undecayed pi0's with status 1 or 2 = " << m_undecayedPi0statuscode12CheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
822  ATH_MSG_INFO(" Event rate with unstable particles with no end vertex = " << m_unstableNoEndVtxCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
823  ATH_MSG_INFO(" Event rate with negative total energy like for tachyonic particles = " << m_negativeEnergyTachyonicCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
824  ATH_MSG_INFO(" Event rate with particles with improper decay properties = " << m_decayCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
825  ATH_MSG_INFO(" Event rate with undisplaced daughters of long lived hadrons = " << m_undisplacedLLHdaughtersCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
826  ATH_MSG_INFO(" Event rate with non zero photon mass = " << m_nonZeroPhotonMassCheckRate*100.0/double(m_nPass + m_nFail) << "% (not included in test efficiency)");
827  ATH_MSG_INFO(" Event rate with no energy balance = " << m_energyBalanceCheckRate*100.0/double(m_nPass + m_nFail) << "%");
828  if (!m_energyImbalanceTest) ATH_MSG_INFO(" The check for energy imbalance is switched off, so is not included in the final TestHepMC efficiency ");
829  ATH_MSG_INFO(" Event rate with no momentum balance = " << m_momentumBalanceCheckRate*100.0/double(m_nPass + m_nFail) << "%");
830  if (!m_momImbalanceTest) ATH_MSG_INFO(" The check for momentum imbalance is switched off, so is not included in the final TestHepMC efficiency ");
831  ATH_MSG_INFO(" Event rate with negative energy particles = " << m_negativeEnergyCheckRate*100.0/double(m_nPass + m_nFail) << "%");
832  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 ");
833  ATH_MSG_INFO(" Event rate with tachyons = " << m_tachyonCheckRate*100.0/double(m_nPass + m_nFail) << "%");
834  if (!m_tachyonsTest) ATH_MSG_INFO(" The check for tachyons is switched off, so is not included in the final TestHepMC efficiency ");
835  ATH_MSG_INFO(" Event rate with stable or unstable particles with no parents = " << m_stableUnstableNoParentCheckRate*100.0/double(m_nPass + m_nFail) << "%");
836  ATH_MSG_INFO(" Event rate with unstable particle with no decay vertex = " << m_unstablePartNoDecayVtxCheckRate*100.0/double(m_nPass + m_nFail) << "%");
837  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 ");
838  ATH_MSG_INFO(" Event rate with undecayed Pi0's = " << m_undecayedPi0CheckRate*100.0/double(m_nPass + m_nFail) << "%");
839  if (!m_pi0NoVtxTest) ATH_MSG_INFO(" The check for undecayed pi0's is switched off, so is not included in the final TestHepMC efficiency ");
840  ATH_MSG_INFO(" Event rate with undisplaced decay daughters of displaced vertices = " << m_undisplacedDecayDaughtersOfDisplacedVtxCheckRate*100.0/double(m_nPass + m_nFail) << "%");
841  if (!m_undisplacedDaughtersTest) ATH_MSG_INFO(" The check for undisplaced daughters is switched off, so is not included in the final TestHepMC efficiency ");
842  ATH_MSG_INFO(" Event rate with particles with status 1 but lifetime < " << m_min_tau << "~ns = " << m_Status1ShortLifetime*100.0/double(m_nPass + m_nFail) << "%");
843  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 ");
844  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/double(m_nPass + m_nFail) << "%");
845  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 ");
846  ATH_MSG_INFO(" Event rate with unknown PDG IDs " << m_unknownPDGIDCheckRate*100.0/double(m_nPass+m_nFail) << "%");
847  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 ");
848 
849  const double tau_fastDrate = double(m_FastDecayedTau) / double(m_TotalTaus);
850  if(tau_fastDrate > m_tau_eff_threshold){
851  ATH_MSG_FATAL("MORE THAN " << 100.*m_tau_eff_threshold << "% OF TAUS DECAYING IMMEDIATELY! " << m_FastDecayedTau << " found, out of: " << m_TotalTaus);
852  return StatusCode::FAILURE;
853  }
854 
855  if (m_noXSECset) {
856  if (m_allowMissingXSec) {
857  ATH_MSG_WARNING(m_noXSECset << " EVENTS WITHOUT CROSS-SECTION!!! Added dummy cross-section instead.");
858  }
859  else {
860  ATH_MSG_FATAL(m_noXSECset << " EVENTS WITHOUT CROSS-SECTION!! Check the setup before production!");
861  return StatusCode::FAILURE;
862  }
863  }
864 
865  const double efficiency = double(m_nPass) / double(m_nPass + m_nFail);
866  ATH_MSG_INFO("Efficiency = " << efficiency * 100 << "%");
867 
868  // Check efficiency, and fail (to kill production jobs) if the pass rate is too low
869  if (efficiency < m_eff_fail_threshold) {
870  ATH_MSG_FATAL("EFFICIENCY ABOVE ERROR THRESHOLD! " << 100*efficiency << "% found, but at least: " << 100*m_eff_fail_threshold << "% required");
871  return StatusCode::FAILURE;
872  } else if (efficiency <= m_eff_warn_threshold) {
873  ATH_MSG_WARNING("EFFICIENCY ABOVE WARNING THRESHOLD! " << 100*efficiency << "% found, but at least: " << 100*m_eff_warn_threshold << "% expected");
874  }
875 
876  return StatusCode::SUCCESS;
877 }
878 
879 
880 #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:17
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
checkFileSG.line
line
Definition: checkFileSG.py:75
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:644
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:710
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
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
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:812
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
AthenaPoolTestWrite.stream
string stream
Definition: AthenaPoolTestWrite.py:12
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
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:210
CaloCondBlobAlgs_fillNoiseFromASCII.desc
desc
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:54
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:807
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
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
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
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:77
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:18
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:379
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:770
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 Tue Apr 8 2025 21:19:53 for ATLAS Offline Software by   doxygen 1.8.18