TestTrackIds_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 // Class: TestTrackIds
3 // Module Type: analyzer
4 // File: TestTrackIds_module.cc
5 ///////////////////////////////////////////////////////////////////////////////
6 
7 // ART includes
14 
15 // NOvA includes
16 #include "MCCheater/BackTracker.h"
18 #include "Simulation/FLSHit.h"
19 #include "Simulation/FLSHitList.h"
21 
22 // ROOT includes
23 #include "TH1.h"
24 
25 namespace cheat
26 {
28  {
29  public:
30  explicit TestTrackIds(const fhicl::ParameterSet& p);
31  virtual ~TestTrackIds();
32 
33  void analyze(const art::Event& evt);
34  void beginJob();
35 
36  private:
37  /// Function that actually aborts job if a particular test failed
38  /// For now, equivalent to a single if condition
39  void CrashIfFCLSaysSo() const;
40 
41  bool fCrashUponError; ///< Flag whether to abort if a particular test failed
42 
43  TH1F* fErrorMode; ///< Histogram breaking down how events failed
44 
45  std::map<std::string, double> fErrorModeMap; ///< Map of bin number to a specific error type
46 
47  static bool EssentiallyEqual(double a, double b, double precision = 0.5)
48  { return a <= (b + precision) && a >= (b - precision); }
49 
50  static bool EssentiallyEqualVertex(TLorentzVector a, TLorentzVector b, double precision = 0.5)
51  {
52  return (EssentiallyEqual(a.X(), b.X(), precision) &&
53  EssentiallyEqual(a.Y(), b.Y(), precision) &&
54  EssentiallyEqual(a.Z(), b.Z(), precision));
55  }
56 
57  };
58 
59  //----------------------------------------------------------------------
61  : EDAnalyzer(pset),
62  fCrashUponError(pset.get<bool>("CrashUponError"))
63  {
64  }
65 
66  //----------------------------------------------------------------------
68  {
69  }
70 
71  //----------------------------------------------------------------------
73  {
75 
76  fErrorMode = tfs->make<TH1F>("fErrorMode", "Error Mode;;Events", 13, 0.5, 13.5);
77 
78  // Start at 1 to match ROOT bin convention
79  // To avoid filling errors, find/replace the strings below,
80  // instead of changing it directly and then searching
81  fErrorModeMap["Real Data"] = 1.;
82  fErrorModeMap["No MCTruth List"] = 2.;
83  fErrorModeMap["No PhotonSignals"] = 3.;
84  fErrorModeMap["No FLSHit match to Track"] = 4.;
85  fErrorModeMap["No FLSHit List"] = 5.;
86  fErrorModeMap["No Particle match to FLSHit"] = 6.;
87  fErrorModeMap["Particle exists, not matched to FLSHit"] = 7.;
88  fErrorModeMap["No Mother match to Track"] = 8.;
89  fErrorModeMap["Mother exists, not matched to Track"] = 9.;
90  fErrorModeMap["No Particle at Track"] = 10.;
91  fErrorModeMap["Missing Daughter"] = 11.;
92  fErrorModeMap["Bad Truth Match"] = 12.;
93  fErrorModeMap["Might Be Cosmic"] = 13.;
94 
95  for(const auto& errorModePair : fErrorModeMap) {
96  int bin = (int)errorModePair.second;
97  std::string label = errorModePair.first;
98  fErrorMode->GetXaxis()->SetBinLabel(bin, label.c_str());
99  }
100  }
101 
102  //----------------------------------------------------------------------
104  {
105  // Check that this is MC, stop if it isn't
106  if(evt.isRealData()) {
107  fErrorMode->Fill(fErrorModeMap["Real Data"]);
108 
110 
111  return;
112  }
113 
114  // Get the MC generator information out of the event
115  bool isMCTruthListEmpty = false;
117  try {
118  evt.getByLabel("generator" /*fMCTruthModuleLabel*/, truths);
119  if(truths->empty()) { isMCTruthListEmpty = true; }
120 
121  if(truths->size() > 1) {
122  mf::LogVerbatim("TestTrackIds") << "Event " << evt.event() << " has "
123  << truths->size() << " truth records." << std::endl;
124  }
125  }
126  catch(...) { isMCTruthListEmpty = true; }
127  if(isMCTruthListEmpty) {
128  mf::LogError("TestTrackIds MCTruth") << "Error retrieving MCTruth list." << std::endl;
129  fErrorMode->Fill(fErrorModeMap["No MCTruth List"]);
130 
132 
133  return;
134  }
135 
137 
138  // Invariant tested here is that every photon signal has a matching FLSHit stored
139  bool isPhotonListError = false;
141  try {
142  evt.getByLabel("photrans" /*fPhotonModuleLabel*/, photlist);
143  }
144  catch(...) {
145  isPhotonListError = true;
146  mf::LogWarning("TestTrackIds PhotonSignal") << "Error retrieving PhotonSignal list." << std::endl;
147  fErrorMode->Fill(fErrorModeMap["No PhotonSignals"]);
148  }
149 
150  if(!isPhotonListError) {
151  std::vector<art::Ptr<sim::PhotonSignal> > phots;
152  art::fill_ptr_vector(phots, photlist);
153 
154  // Loop over all photons
155  // phot should be a sim::PhotonSignal
156  for(const auto& phot : phots) {
157  // Get the FLS hits from the photon
158  const std::vector<sim::FLSHit> flslist = bt->CellToFLSHit(phot->Plane(),
159  phot->Cell());
160 
161  // Loop over FLS hits
162  // Check that the FLS hit and photon have the same track; one match is enough
163  bool matched = false;
164  for(const auto& flshit : flslist) {
165  if(flshit.GetTrackID() == phot->TrackId()) {
166  matched = true;
167  break;
168  }
169  }
170 
171  // Log an error in the (unlikely) event the tracks don't match
172  if(!matched) {
173  mf::LogError("TestTrackIds") << "Event " << evt.event()
174  << ": PhotonSignal in plane " << phot->Plane()
175  << ", cell " << phot->Cell()
176  << ", TrackId " << phot->TrackId()
177  << " does not have any matching FLSHits." << std::endl;
178  fErrorMode->Fill(fErrorModeMap["Track has no FLSHits"]);
179 
181  }
182  } // end of loop over PhotonSignals
183  } // end of conditional that there was no error retrieving the Photon List
184 
185  // Now get the sim::Particles from the event
186  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
187 
188  // Used for a sanity checks of pnav
189  const std::set<int> ids = bt->GetTrackIDList();
190 
191  // Invariant here is that every FLSHit's TrackId corresponds to a particle stored in the event
192  bool isFLSHitListError = false;
194  try {
195  evt.getByLabel("geantgen" /*fGeantModuleLabel*/, flshitlist);
196  }
197  catch(...) {
198  isFLSHitListError = true;
199  mf::LogError("TestTrackIds FLSHitList") << "Error retrieving FLSHitList list." << std::endl;
200  fErrorMode->Fill(fErrorModeMap["No FLSHit List"]);
201  }
202 
203  if(!isFLSHitListError) {
204  for(unsigned int i_List = 0, n_List = flshitlist->size(); i_List < n_List; ++i_List) {
205  const std::vector<sim::FLSHit>& hits = flshitlist->at(i_List).fHits;
206 
207  for(const auto& hit : hits){
208  const sim::ParticleNavigator::const_iterator part = pnav.find(hit.GetTrackID());
209  if(part == pnav.end()){
210  mf::LogError("TestTrackIds") << "Event " << evt.event()
211  << ": FLS hit in plane " << hit.GetPlaneID()
212  << ", cell " << hit.GetCellID()
213  << ", trackId " << hit.GetTrackID()
214  << ", PDG code " << hit.GetPDG()
215  << " has no matching sim::Particle" << std::endl;
216  fErrorMode->Fill(fErrorModeMap["No Particle match to FLSHit"]);
217 
218  // Make sure it really isn't there
219  if(ids.count(hit.GetTrackID()) != 0) {
220  mf::LogWarning("TestTrackIds") << "However, the sim::Particle does exist." << std::endl;
221  fErrorMode->Fill(fErrorModeMap["Particle exists, not matched to FLSHit"]);
222  }
223 
225  }
226  } // end of loop over FLSHits
227  } // end of loop over lists in the FLSHitList
228  } // end of conditional that there was no error retrieving the FLSHitList
229 
230  int badMatch = 0;
231  std::vector<NeutrinoWithIndex> allNeutrinos = bt->allNeutrinoInteractions();
232 
233  // Invariant tested here is that every non-primary particle has its mother stored
234  for(sim::ParticleNavigator::const_iterator it = pnav.begin(); it != pnav.end(); ++it) {
235  const int id = it->first;
236 
237  if(pnav.IsPrimary(id)) {
238  // Test that primary muons match vertex location of neutrino mother
239  const sim::Particle* muon = bt->TrackIDToParticle(id);
240  int pdg = muon->PdgCode();
241  int mother = muon->Mother();
242  TLorentzVector muonVert = muon->Position();
243  double muonE = muon->E();
244 
245  const art::Ptr<simb::MCTruth>& motherTruth = bt->ParticleToMCTruth(muon);
246  if(motherTruth->NeutrinoSet()) {
247  TLorentzVector motherVert = motherTruth->GetNeutrino().Nu().Position();
248  double lepE = motherTruth->GetNeutrino().Lepton().E();
249 
250  if(!EssentiallyEqualVertex(muonVert, motherVert)) {
251  // Want to find neutrino in vector
252  int matchedNeutrino = -1;
253  int muonNeutrino = -1;
254 
255  for(unsigned int q = 0; q < allNeutrinos.size(); ++q) {
256  TLorentzVector otherVert = allNeutrinos[q].neutrinoInt->GetNeutrino().Nu().Position();
257 
258  if(EssentiallyEqualVertex(otherVert, motherVert)) {
259  matchedNeutrino = allNeutrinos[q].truthColIndex;
260  } // Matches our neutrino truth
261 
262  if(EssentiallyEqualVertex(otherVert, muonVert)) {
263  muonNeutrino = allNeutrinos[q].truthColIndex;
264  } // Matches our muon vertex
265  } // end of loop over vector of NeutrinoWithIndex
266 
267  // Sometimes primary particles travel away from the vertex somewhat and don't match any truth perfectly ... don't worry
268  // Special cased particles: charged pions, photons, protons, neutrons, K longs, charged kaons
269  // These were choosen because I saw examples of this behaviour. Others could be added.
270  std::set<int> specialCases = {211, -211, 22, 2212, 130, -321, 321, 2112};
271  if(!(specialCases.find(pdg) != specialCases.end() && muonNeutrino == -1)) {
272  ++badMatch;
273  mf::LogError("TestTrackIds") << " PDG: " << pdg << " Mother: " << mother << std::endl
274  << " Particle Vx: " << muon->Vx()
275  << " Vy: " << muon->Vy()
276  << " Vz: " << muon->Vz()
277  << " Ptc E: " << muonE << std::endl
278  << " Neutrino Vx: " << motherTruth->GetNeutrino().Nu().Vx()
279  << " Vy: " << motherTruth->GetNeutrino().Nu().Vy()
280  << " Vz: " << motherTruth->GetNeutrino().Nu().Vz()
281  << " Lep E: " << lepE << std::endl
282  << " Neutrino matched index: " << matchedNeutrino << std::endl
283  << " Ptc. vtx matched index: " << muonNeutrino << std::endl;
284 
285  fErrorMode->Fill(fErrorModeMap["Bad Truth Match"]);
287  } // We are sad that they don't match - not just because particle traveled somewhat
288  } // Vertex doesn't match between primary particle and neutrino truth
289  }
290  else {
291  mf::LogWarning("Might Be Cosmic") << "The mother truth had no neutrino."
292  << "This is a problem if the file is not cosmics." << std::endl;
293 
294  fErrorMode->Fill(fErrorModeMap["Might Be Cosmic"]);
296  }
297 
298  continue;
299  } // Skip the particle if it is a primary
300 
301  const int mom = it->second->Mother();
302  const int pdg = it->second->PdgCode();
303  const std::string proc = it->second->Process();
304 
305  if(pnav.find(mom) == pnav.end()) {
306  mf::LogError("TestTrackIds") << "Event " << evt.event()
307  << ": Particle " << id
308  << ", PDG code " << pdg
309  << ", formed by process " << proc
310  << ", can't find mother, index " << mom << std::endl;
311  fErrorMode->Fill(fErrorModeMap["No Mother match to Track"]);
312 
313  // Make sure it really isn't there
314  if(ids.count(mom) != 0) {
315  mf::LogWarning("TestTrackIds") << "However, the mother does exist." << std::endl;
316  fErrorMode->Fill(fErrorModeMap["Mother exists, not matched to Track"]);
317 
319  }
320  }
321  } // end of loop over particles
322 
323  if(badMatch != 0) {
324  mf::LogError("TestTrackIds") << " Number of bad matches of primary particles to simb truths is: " << badMatch << std::endl;
325  }
326 
327  // Invaraint tested here is that all daughters of a particle should exist
328  std::set<int> trackIDs = bt->GetTrackIDList();
329  for(const auto& trackID : trackIDs) {
330  const sim::Particle* mother = bt->TrackIDToParticle(trackID);
331 
332  if(!mother) {
333  mf::LogError("TestTrackIds") << "Event " << evt.event()
334  << ", can't find particle corresponding to track id: " << trackID << std::endl;
335  fErrorMode->Fill(fErrorModeMap["No Particle at Track"]);
336 
338 
339  continue;
340  }
341 
342  for(int i_daughter = 0, n_daughter = mother->NumberDaughters(); i_daughter < n_daughter; ++i_daughter) {
343  const int id = mother->Daughter(i_daughter);
344 
345  if(!bt->TrackIDToParticle(id)) {
346  mf::LogError("TestTrackIds") << "Event " << evt.event()
347  << ", can't find the daughter of particle: " << trackID
348  << ", and the daughter has track id: " << id << std::endl;
349  fErrorMode->Fill(fErrorModeMap["Missing Daughter"]);
350 
352  }
353  } // end of loop over the daughters of the current particle
354 
355  } // end of loop over all track ids
356 
357  } // end of analyze
358 
359  //----------------------------------------------------------------------
361  {
362  if(fCrashUponError) { abort(); }
363  return;
364  }
365 
367 } // end of namespace cheat
double E(const int i=0) const
Definition: MCParticle.h:232
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
set< int >::iterator it
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
TH1F * fErrorMode
Histogram breaking down how events failed.
int Mother() const
Definition: MCParticle.h:212
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
list_type::const_iterator const_iterator
std::vector< NeutrinoWithIndex > allNeutrinoInteractions() const
Function primarily for CAFMaker - returns a vector of all MCTruth interactions along with their truth...
static bool EssentiallyEqual(double a, double b, double precision=0.5)
TestTrackIds(const fhicl::ParameterSet &p)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
int NumberDaughters() const
Definition: MCParticle.h:216
void CrashIfFCLSaysSo() const
int Daughter(const int i) const
Definition: MCParticle.cxx:112
const char * label
static bool EssentiallyEqualVertex(TLorentzVector a, TLorentzVector b, double precision=0.5)
TString part[npart]
Definition: Style.C:32
void hits()
Definition: readHits.C:15
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
const double a
const art::Ptr< simb::MCTruth > & ParticleToMCTruth(const sim::Particle *p) const
EventNumber_t event() const
Definition: Event.h:67
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
float bin[41]
Definition: plottest35.C:14
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
bool fCrashUponError
Flag whether to abort if a particular test failed.
bool IsPrimary(int trackID) const
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:16
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: event.h:1
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
const hit & b
Definition: hits.cxx:21
const std::set< int > GetTrackIDList() const
Get all G4 track ids present in the event.
Definition: BackTracker.h:750
void analyze(const art::Event &evt)
double Vz(const int i=0) const
Definition: MCParticle.h:222
std::map< std::string, double > fErrorModeMap
Map of bin number to a specific error type.
const std::vector< sim::FLSHit > CellToFLSHit(unsigned int const &plane, unsigned int const &cell) const
Returns the FLSHits contributing the signal in the specified cell. WARNING: Use with extreme caution!...
bool NeutrinoSet() const
Definition: MCTruth.h:77
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Int_t trackID
Definition: plot.C:84
iterator find(const key_type &key)
double Vy(const int i=0) const
Definition: MCParticle.h:221