G4MismatchAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: G4MismatchAna
3 // Plugin Type: analyzer (art v2_13_00)
4 // File: G4MismatchAna_module.cc
5 //
6 // Generated at Wed Mar 4 11:29:05 2020 by Jeremy Hewes using cetskelgen
7 // from cetlib version v3_06_01.
8 ////////////////////////////////////////////////////////////////////////
9 
17 #include "fhiclcpp/ParameterSet.h"
20 
23 #include "Simulation/Particle.h"
25 
26 #include "TTree.h"
27 
28 class G4MismatchAna;
29 
30 
32 public:
33  explicit G4MismatchAna(fhicl::ParameterSet const & p);
34  // The compiler-generated destructor is fine for non-base
35  // classes without bare pointers or other resource use.
36 
37  // Plugins should not be copied or assigned.
38  G4MismatchAna(G4MismatchAna const &) = delete;
39  G4MismatchAna(G4MismatchAna &&) = delete;
40  G4MismatchAna & operator = (G4MismatchAna const &) = delete;
42 
43  // Required functions.
44  void analyze(art::Event const & e) override;
45 
46 private:
47 
48  // Module labels
52 
53  // Information tree
54  TTree* fGENIETree;
55  TTree* fSliceTree;
56 
57  // Tree branches
58  int fRun; // Run number
59  int fSubRun; // SubRun number
60  int fEvent; // Event number
61  int fCycle; // Cycle number
62  int fTruth; // GENIE truth index
63  int fSlice; // Slice index
64 
65  bool fIsMismatched; // Whether GENIE primaries are missing from the G4 record
66  double fNuEnergy; // True neutrino energy
67  int fNuPDG; // True neutrino PDG code
68  bool fIsLeptonMissing; // Whether final state lepton is missing
69  double fLepEnergy; // True final state lepton energy
70  int fLepPDG; // True final state lepton PDG code
71 
72  double fW; // Hadronic invariant mass
73  double fX; // Bjorken X
74  double fY; // Inelasticity
75  double fQSqr; // Momentum transfer
76 
77  int fMissingParticles; // How many GENIE particles are missing from G4 truth
78  double fMissingEnergy; // Total energy of GENIE particles missing from G4 truth
79 
80  double fNumuCosRej; // SR value sel.cosrej.numucontpid2020
81  double fNumuReMID; // SR value sel.remid.pid
82  double fNumuCVN; // SR value sel.cvnloosepreselptp.numuid
83 
84  double fNueCosRejCore; // SR value sel.nuecosrej.cospidcorebdt
85  double fNueCosRejPeri; // SR value sel.nuecosrej.cospidperibdt
86  double fNueCVN; // SR value sel.cvnloosepreselptp.nueid
87 
88  double fNCCVN; // SR value sel.cvnloosepreselptp.ncid
89 
90 };
91 
92 
94  :
95  EDAnalyzer(p),
96  fGenLabel(p.get<std::string>("GenLabel")),
97  fG4Label(p.get<std::string>("G4Label")),
98  fSRLabel(p.get<std::string>("SRLabel"))
99 {
101 
102  // Set up input tree
103  fGENIETree = tfs->make<TTree>("G4MismatchGENIE", "G4MismatchGENIE");
104 
105  fGENIETree->Branch("Run", &fRun);
106  fGENIETree->Branch("SubRun", &fSubRun);
107  fGENIETree->Branch("Event", &fEvent);
108  fGENIETree->Branch("Cycle", &fCycle);
109  fGENIETree->Branch("Truth", &fTruth);
110  fGENIETree->Branch("IsMismatched", &fIsMismatched);
111  fGENIETree->Branch("NuEnergy", &fNuEnergy);
112  fGENIETree->Branch("NuPDG", &fNuPDG);
113  fGENIETree->Branch("IsLeptonMissing", &fIsLeptonMissing);
114  fGENIETree->Branch("LepEnergy", &fLepEnergy);
115  fGENIETree->Branch("LepPDG", &fLepPDG);
116  fGENIETree->Branch("W", &fW);
117  fGENIETree->Branch("X", &fX);
118  fGENIETree->Branch("Y", &fY);
119  fGENIETree->Branch("QSqr", &fQSqr);
120  fGENIETree->Branch("MissingParticles", &fMissingParticles);
121  fGENIETree->Branch("MissingEnergy", &fMissingEnergy);
122 
123  fSliceTree = tfs->make<TTree>("G4MismatchSlice", "G4MismatchSlice");
124 
125  fSliceTree->Branch("Run", &fRun);
126  fSliceTree->Branch("SubRun", &fSubRun);
127  fSliceTree->Branch("Event", &fEvent);
128  fSliceTree->Branch("Cycle", &fCycle);
129  fSliceTree->Branch("Slice", &fSlice);
130  fSliceTree->Branch("IsMismatched", &fIsMismatched);
131  fSliceTree->Branch("NuEnergy", &fNuEnergy);
132  fSliceTree->Branch("NuPDG", &fNuPDG);
133  fSliceTree->Branch("IsLeptonMissing", &fIsLeptonMissing);
134  fSliceTree->Branch("LepEnergy", &fLepEnergy);
135  fSliceTree->Branch("LepPDG", &fLepPDG);
136  fSliceTree->Branch("W", &fW);
137  fSliceTree->Branch("X", &fX);
138  fSliceTree->Branch("Y", &fY);
139  fSliceTree->Branch("QSqr", &fQSqr);
140  fSliceTree->Branch("MissingParticles", &fMissingParticles);
141  fSliceTree->Branch("MissingEnergy", &fMissingEnergy);
142  fSliceTree->Branch("NumuCosRej", &fNumuCosRej);
143  fSliceTree->Branch("NumuReMID", &fNumuReMID);
144  fSliceTree->Branch("NumuCVN", &fNumuCVN);
145  fSliceTree->Branch("NueCosRejCore", &fNueCosRejCore);
146  fSliceTree->Branch("NueCosRejPeri", &fNueCosRejCore);
147  fSliceTree->Branch("NueCVN", &fNueCVN);
148  fSliceTree->Branch("NCCVN", &fNCCVN);
149 }
150 
152 {
153  // Get GENIE information
155  try { e.getByLabel(fGenLabel, truthList); }
156  catch(...) { throw art::Exception(art::errors::ProductNotFound) << "Failed to get MC truth information."; }
157  std::vector<art::Ptr<simb::MCTruth>> truthPtrs;
158  art::fill_ptr_vector(truthPtrs, truthList);
159  // const simb::MCTruth* truth = truthPtrs[0].get();
160 
161  // Get G4 information
163  try { e.getByLabel(fG4Label, g4List); }
164  catch(...) { throw art::Exception(art::errors::ProductNotFound) << "Failed to get simulated particles."; }
165  std::vector<art::Ptr<sim::Particle>> particles;
166  art::fill_ptr_vector(particles, g4List);
167 
168  // Get SR information
170  try { e.getByLabel(fSRLabel, srList); }
171  catch(...) { throw art::Exception(art::errors::ProductNotFound) << "Failed to get standard record."; }
172  std::vector<art::Ptr<caf::StandardRecord>> srPtrs;
173  art::fill_ptr_vector(srPtrs, srList);
174 
175  // Reset tree variables that need resetting
176  fIsMismatched = false;
177  fIsLeptonMissing = true;
178  fMissingParticles = 0;
179  fMissingEnergy = 0;
180  fCycle = -5;
181 
182  // Fill event info
183  fRun = e.id().run();
184  fSubRun = e.id().subRun();
185  fEvent = e.id().event();
186 
187  // Loop over G4 particles to find primary lepton
188  for (size_t i = 0; i < truthPtrs.size(); ++i) {
189 
190  const simb::MCTruth* truth = truthPtrs[i].get();
191 
192  fTruth = i;
193 
194  // Neutrino and lepton truth information
195  fIsMismatched = false;
196  fIsLeptonMissing = true;
197  simb::MCNeutrino nu = truth->GetNeutrino();
198  fNuEnergy = nu.Nu().E();
199  fNuPDG = nu.Nu().PdgCode();
200  fLepEnergy = nu.Lepton().E();
201  fLepPDG = nu.Lepton().PdgCode();
202  fW = nu.W();
203  fX = nu.X();
204  fY = nu.Y();
205  fQSqr = nu.QSqr();
206 
207  for (art::Ptr<sim::Particle> part : particles) {
208  if (part->Mother() == 0 && part->PdgCode() == fLepPDG && fabs(part->E()-fLepEnergy) < 0.01) {
209  fIsLeptonMissing = false; // found the lepton in G4 record
210  break;
211  }
212  }
213 
214  // Construct a list of GENIE primary particles
215  std::vector<std::pair<int, double>> primaries;
216  for (size_t i = 0; i < (size_t)truth->NParticles(); ++i) {
218  if (part.StatusCode() == 1 && part.PdgCode() < 2000000000) {
219  primaries.push_back(std::make_pair(part.PdgCode(), part.E()));
220  }
221  }
222 
223  // Look for all primaries in G4 truth list
224  for (auto primary : primaries) {
225  bool found = false;
226  for (art::Ptr<sim::Particle> part : particles) {
227  // If it's a G4 particle with the same PDG code and within 10 MeV energy,
228  // call them matched. This isn't infallible but the chances of false positives
229  // should be small compared to the rate at which mismatches occur
230  if (part->Mother() == 0 && part->PdgCode() == primary.first
231  && fabs(part->E()-primary.second) < 0.01) {
232  found = true;
233  break; // if we found this primary, don't need to keep looping
234  }
235  } // for G4 particle
236  if (!found) { // if there's a missing particle, then make note of that
237  fIsMismatched = true;
239  fMissingEnergy += primary.second;
240  } // if particle is missing
241  } // for primary
242 
243  // Spit out a warning message if we have a mismatch here
244  if (fIsMismatched) {
246  << "////////////////////////////////////////////" << std::endl
247  << "/ /" << std::endl
248  << "/ GENIE VS GEANT MISMATCH FOUND! /" << std::endl
249  << "/ /" << std::endl
250  << "////////////////////////////////////////////" << std::endl
251  << std::endl;
252  }
253 
254  fGENIETree->Fill();
255 
256  }
257 
258  for (size_t i = 0; i < srPtrs.size(); ++i) {
259 
260  const caf::StandardRecord* sr = srPtrs[i].get();
261 
262  // If we haven't matched any GENIE truth to this slice, just skip it
263  if (sr->mc.nu.size() == 0) continue;
264 
265  // Loop over generators to find the one that matches this slice
266  const simb::MCTruth* tru = nullptr;
267  for (art::Ptr<simb::MCTruth> truth : truthPtrs) {
268  TVector3 genPos = truth->GetNeutrino().Nu().Position().Vect();
269  TVector3 slcPos = sr->mc.nu[0].vtx;
270  TVector3 diff = genPos - slcPos;
271  if (diff.Mag() < 0.01) { // If vertices match within 0.01cm, call it a match
272  tru = truth.get();
273  break; // Found what we were looking for
274  }
275  }
276  if (!tru) throw art::Exception(art::errors::NotFound) << "No valid GENIE truth found for this slice.";
277 
278  fSlice = i;
279  fCycle = sr->hdr.cycle;
280 
282  fNumuReMID = sr->sel.remid.pid;
284 
288 
290 
291 
292  // Neutrino and lepton truth information
293  fIsMismatched = false;
294  fIsLeptonMissing = true;
295  fMissingParticles = 0;
296  fMissingEnergy = 0;
297  simb::MCNeutrino nu = tru->GetNeutrino();
298  fNuEnergy = nu.Nu().E();
299  fNuPDG = nu.Nu().PdgCode();
300  fLepEnergy = nu.Lepton().E();
301  fLepPDG = nu.Lepton().PdgCode();
302  fW = nu.W();
303  fX = nu.X();
304  fY = nu.Y();
305  fQSqr = nu.QSqr();
306 
307  // Loop over G4 particles to find primary lepton
308  for (size_t j = 0; j < sr->mc.nu[0].prim.size(); ++j) {
309  caf::SRTrueParticle p = sr->mc.nu[0].prim[j];
310  if (p.pdg == fLepPDG && fabs(p.p.E-fLepEnergy) < 0.01) {
311  fIsLeptonMissing = false; // found the lepton in G4 record
312  break;
313  }
314  }
315 
316  // Get GENIE primaries for this specific interaction
317  std::vector<std::pair<int, double>> primaries;
318  for (size_t i = 0; i < (size_t)tru->NParticles(); ++i) {
320  if (part.StatusCode() == 1 && part.PdgCode() < 2000000000) {
321  primaries.push_back(std::make_pair(part.PdgCode(), part.E()));
322  }
323  }
324 
325  // Loop over GENIE primaries
326  for (auto primary : primaries) {
327  bool found = false;
328  // Loop over G4 particles in StandardRecord
329  for (size_t j = 0; j < sr->mc.nu[0].prim.size(); ++j) {
330  caf::SRTrueParticle p = sr->mc.nu[0].prim[j];
331  if (primary.first == p.pdg && fabs(primary.second-p.p.E) < 0.01) {
332  found = true;
333  break;
334  }
335  }
336  if (!found) { // If there's a particle missing, keep track of what it is
337  fIsMismatched = true;
339  fMissingEnergy += primary.second;
340  }
341  }
342  fSliceTree->Fill();
343 
344  } // for slice i
345 
346 }
347 
double E(const int i=0) const
Definition: MCParticle.h:232
float ncid
Likelihood Neutral Current.
Definition: SRCVNResult.h:23
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
double QSqr() const
Definition: MCNeutrino.h:157
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
SRHeader hdr
Header branch: run, subrun, etc.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
SRCosRej cosrej
Output from CosRej (Cosmic Rejection)
Definition: SRIDBranch.h:47
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
int StatusCode() const
Definition: MCParticle.h:210
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
SRCVNResult cvnloosepreselptp
Output from CVN - Loose Presel plus PtP cut (many-class PID)
Definition: SRIDBranch.h:54
int NParticles() const
Definition: MCTruth.h:74
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: EventID.h:98
float nueid
Likelihood Charge Current NuE.
Definition: SRCVNResult.h:21
G4MismatchAna(fhicl::ParameterSet const &p)
float cospidperibdt
Nue cosrej PID for the peripheral sample for 2020+.
Definition: SRNueCosRej.h:126
TString part[npart]
Definition: Style.C:32
float pid
PID value output by kNN.
Definition: SRRemid.h:25
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
double W() const
Definition: MCNeutrino.h:154
double Y() const
Definition: MCNeutrino.h:156
SRNueCosRej nuecosrej
Output from NueCosRej (Nue Cosmic Rejection)
Definition: SRIDBranch.h:48
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
float cospidcorebdt
Nue cosrej PID for the core sample for 2020+.
Definition: SRNueCosRej.h:128
G4MismatchAna & operator=(G4MismatchAna const &)=delete
caf::StandardRecord * sr
double X() const
Definition: MCNeutrino.h:155
const double j
Definition: BetheBloch.cxx:29
void analyze(art::Event const &e) override
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
OStream cout
Definition: OStream.cxx:6
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
EventNumber_t event() const
Definition: EventID.h:116
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
The StandardRecord is the primary top-level object in the Common Analysis File trees.
SubRunNumber_t subRun() const
Definition: EventID.h:110
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
SRIDBranch sel
Selector (PID) branch.
float numuid
Likelihood Charge Current NuMu.
Definition: SRCVNResult.h:20
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
Event generator information.
Definition: MCTruth.h:32
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
float numucontpid2020
cosmic rejection PID for contained events; 2020 Analysis
Definition: SRCosRej.h:34
Float_t e
Definition: plot.C:35
Event generator information.
Definition: MCNeutrino.h:18
Definition: fwd.h:28
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
EventID id() const
Definition: Event.h:56
int cycle
MC simulation cycle number.
Definition: SRHeader.h:23
SRLorentzVector p
Momentum 4-vector.
enum BeamMode string