NumuEAna_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 // Class: NumuEAna
3 // Module Type: analyzer
4 // File: NumuEAna_module.cc
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TTree.h"
10 
16 
17 #include "MCCheater/BackTracker.h"
19 #include "CosRej/CosRejObj.h"
20 #include "RecoBase/Track.h"
21 #include "RecoBase/Cluster.h"
22 #include "RecoBase/CellHit.h"
23 #include "RecoBase/Energy.h"
24 #include "Simulation/Simulation.h"
25 #include "Simulation/Particle.h"
26 #include "Geometry/Geometry.h"
29 #include "Utilities/AssociationUtil.h"
30 #include "Utilities/func/MathUtil.h"
31 #include "ReMId/ReMId.h"
32 #include "QEEventFinder/QePId.h"
33 #include "NumuEnergy/NumuE.h"
34 #include "NovaDAQConventions/DAQConventions.h"
35 
36 class TH1D;
37 class TH2D;
38 class TTree;
39 
40 namespace numue {
41 
42  class NumuEAna : public art::EDAnalyzer {
43  public:
44  explicit NumuEAna(fhicl::ParameterSet const &p);
45  virtual ~NumuEAna();
46 
47  virtual void analyze(art::Event const &e);
48 
49  virtual void reconfigure(fhicl::ParameterSet const & p);
50  virtual void beginJob();
51  virtual void beginRun(art::Run const &r);
53 
54  private:
55 
56  std::string fHitModuleLabel; ///< lable for module making the hits
57  std::string fTrackModuleLabel; ///< label for module making the tracks
58  std::string fCosmicTrackModuleLabel; ///< label for module making cosmic tracks
59  std::string fSlicerModuleLabel; ///< label for module making the slices
60  std::string fEnergyModuleLabel; ///< label for module writing energy objects
61  std::string fPIDModuleLabel; ///< label for module writing pid objects
62  std::string fCosRejLabel; ///< label for module writing cosmic rejection
63  std::string fQepidLabel; ///< label for module writing qe pid objects
64 
66 
67  TTree* fNumuETree; ///< n-tuple used for further analyses (pure population)
68  TH1D* fCalCCEnergy; ///< histogram of calorimetric neutrino energy for general CC population
69  TH1D* fTrkQEEnergy; ///< histogram of muon track length neutrino energy for QE super pure population
70  TH1D* fTrkNonQEEnergy; ///< histogram of muon track length neutrino energy for NonQE super pure population
71  TH1D* fAngleQEEnergy; ///< histogram of neutrino energy for QE population using angle formula
72 
73  TH1D* fAngleQEError; ///< histogram of error on the neutrino energy for QE population using angle formula
74 
75  TH1D* fCalCCEnergyDiff; ///< histogram of the diff. between calorimetric neutrino E and true neutrino E (for pure pop.)
76  TH1D* fTrkQEEnergyDiff; ///< histogram of the diff. between muon track length neutrino E and true neutrino E (for pure QE pop.)
77  TH1D* fTrkNonQEEnergyDiff; ///< histogram of the diff. between muon track length neutrino E and true neutrino E (for pure NonQE pop.)
78  TH1D* fAngleQEEnergyDiff; ///< histogram of the diff. between neutrino E using angle formula and true neutrino E (for pure QE pop.)
79 
80  TH2D* fCalCCEnergy2D; ///< histogram of calorimetric neutrino E and true neutrino E (for pure pop.)
81  TH2D* fTrkQEEnergy2D; ///< histogram of muon track length neutrino E and true neutrino E (for pure QE pop.)
82  TH2D* fTrkNonQEEnergy2D; ///< histogram of muon track length neutrino E and true neutrino E (for pure NonQE pop.)
83  TH2D* fAngleQEEnergy2D; ///< histogram of neutrino E using angle formula and true neutrino E (for pure QE pop.)
84 
85  TH1D* fHadronCalEnergy; ///< histogram of hadronic calorimetric energy NOT on the muon track
86  TH1D* fHadronTrkEnergy; ///< histogram of hadronic calorimetric energy on the muon track
87 
88  TH1D* fNDTrkLenActive; ///< histogram of ND muon track length in the active region
89  TH1D* fNDTrkLenCatcher; ///< histogram of ND muon track length in the muon catcher region
90  TH1D* fNDTrkCalActive; ///< histogram of ND muon calorimetric energy in the active region
91  TH1D* fNDTrkCalTransition; ///< histogram of ND muon calorimetric energy in the transition region
92  TH1D* fNDTrkCalCatcher; ///< histogram of ND muon calorimetric energy in the muon catcher region
93  TH1D* fNDHadCalActive; ///< histogram of ND hadron calorimetric energy in the active region
94  TH1D* fNDHadCalTransition; ///< histogram of ND hadron calorimetric energy in the transition region
95  TH1D* fNDHadCalCatcher; ///< histogram of ND hadron calorimetric energy in the muon catcher region
96 
97  // Main energy estimators
98  float CalCCEnergy ;
99  float QEEnergy ;
100  float NonQEEnergy ;
101  float CCEnergy ;
102  float AngleQEEnergy ;
103  float AngleQEError ;
104  float HadCalEnergy ;
105  float HadTrkEnergy ;
106  float ndTrkLenAct ;
107  float ndTrkLenCat ;
108  float ndTrkCalAct ;
109  float ndTrkCalTran ;
110  float ndTrkCalCat ;
111  float ndHadCalAct ;
112  float ndHadCalTran ;
113  float ndHadCalCat ;
114  float ndHadTrkAct ;
115  float ndHadTrkTran ;
116  float ndHadTrkCat ;
117  float ndTrkTranX ;
118  float ndTrkTranY ;
119 
120  bool fND; ///< is detector ND?
121  bool fFD; ///< is detector FD?
122 
123  bool isCC; ///< is the interaction CC?
124  bool TrueQEInt; ///< is the interaction QE?
125  bool trackMatchesMuon; ///< track is muon-like
126  bool contained; ///< does it pass numu containment cuts?
127 
128  int Run;
129  int SubRun;
130  int Evt;
131  int SubEvt;
132  int goodMuon;
134  int PDGCode;
135  int nKalman;
136  int nCosmic;
137  int qepidntrk; ///< Number of tracks used by QePId
138 
139  // All of these have to do with containment
144  unsigned int mincellsfromedge;
147 
148  double NuPurity; ///< Neutrino selection purity
149  double NuEfficiency; ///< Neutrino selection efficiency
150  double trackLength; ///< Reconstructed tracklength
151  double catcherE; ///< Energy before entering the muon catcher (ND)
152  double bestPID; ///< Highest RemID
153  double qepidPID; ///< PID value of output kNN from QePId
154  double cosrejanglekal; ///< Cosine of angle of best ReMId Kalman Track from CosRej
155  double cosrejnumucontpid; ///< Cosmic rejection PID for contained events from CosRej
158  double start_x_reco;
159  double start_y_reco;
160  double start_z_reco;
161  double end_x_reco;
162  double end_y_reco;
163  double end_z_reco;
164  double TrueMuonE; ///< True energy of the muon
165  double RecoMuonE; ///< Reco energy of the muon
166  double TrueNuE; ///< True energy of the neutrino
167  double TrueTrackLength; ///< True tracklength of the muon
168  double fMaxTrkLen; ///< length of detector from corner to corner
169  };
170 
171  //-------------------------------------------------------------------
173  EDAnalyzer(p),
174  fND(false),
175  fFD(false),
176  fMaxTrkLen(-1)
177  {
178  this->reconfigure(p);
179  }
180 
181  //-------------------------------------------------------------------
183  {
184  }
185 
186  //-------------------------------------------------------------------
188  {
189  if(!(fND || fFD)){
190  mf::LogWarning("NumuEAna") << "attempting to run on detector that is not "
191  << "nd or fd, bail";
192  return;
193  }
194 
196 
197  Run = e.run();
198  SubRun = e.subRun();
199  Evt = e.event();
200 
201 
202  // get all hits in the event
204  e.getByLabel(fHitModuleLabel, hithdl);
205 
206 
208  for(size_t h = 0; h < hithdl->size(); ++h){
209  art::Ptr<rb::CellHit> hit(hithdl, h);
210  allhits.push_back(hit);
211  }
212 
213  // get the tracks out of the eventpdg
215  e.getByLabel(fTrackModuleLabel,trkcol);
216 
217  //Load the slicer list from the event
219  e.getByLabel(fSlicerModuleLabel,slicecol);
220  if(slicecol->empty()){
221  mf::LogWarning ("No Slices")<<"No Slices in the input file";
222  return;
223  }
224 
225  art::PtrVector<rb::Cluster> slicelist;
226  for(unsigned int i = 0; i<slicecol->size();++i){
227  art::Ptr<rb::Cluster> clust(slicecol,i);
228  slicelist.push_back(clust);
229  }
230 
232  if(!bt->HaveTruthInfo()){
233  mf::LogWarning("NumuEAna") << "attempting to run MC truth check on "
234  << "things without truth, bail";
235  return;
236  }
237 
238  //Associating slices to tracks and energies
239  art::FindManyP<rb::Track> sliceToTracks(slicecol, e, fTrackModuleLabel);
240  art::FindManyP<rb::Track> sliceToCosmicTracks(slicecol, e, fCosmicTrackModuleLabel);
241  art::FindOneP<cosrej::CosRejObj> sliceToCosRej(slicecol, e, fCosRejLabel);
242  art::FindOneP<NumuE> sliceToEnergy(slicecol, e, fEnergyModuleLabel);
243  art::FindOneP<qeef::QePId> sliceToQepid(slicecol, e, fQepidLabel);
244 
245  //Loop over slices
246  for(size_t s = 0; s < slicelist.size(); ++s){
247 
248  SubEvt = s;
249 
250  const art::Ptr<NumuE> sliceEnergy = sliceToEnergy.at(s);
251  if (!sliceEnergy) continue;
252 
253  const art::Ptr<cosrej::CosRejObj> sliceCosrej = sliceToCosRej.at(s);
254  if (!sliceCosrej) continue;
255 
256  cosrejanglekal = sliceCosrej->AngleKal();
257  cosrejnumucontpid = sliceCosrej->ConCosPID();
258 
259  const art::Ptr<qeef::QePId> sliceQepid = sliceToQepid.at(s);
260  qepidPID = -5;
261  qepidntrk = -5;
262  if (sliceQepid){
263  qepidPID = sliceQepid->Value();
264  qepidntrk = sliceQepid->Ntrk();
265  }
266 
267  // A great exercise of copy and paste starts here
268  unsigned int dibmask = (rh->GoodDiBlockMask()&rh->GetConfiguration());
269  if(dibmask == ((1<<16)-1)) {
270  dibmask = rh->GetConfiguration();
271  }
272  if(dibmask == 0) {
273  dibmask = rh->GetConfiguration();
274  }
275 
276  unsigned int dibfirst = 0;
277  unsigned int diblast = 0;
278  if (dibmask) {
279  int iD;
280  iD=0; while (!((dibmask>>iD)&1)) iD++; dibfirst=iD+1;
281  iD=0; while (dibmask>>iD) iD++; diblast=iD;
282  }
283 
284  unsigned int cellsfromedge = 99999999;
285  mincellsfromedge = 99999999;
286 
287  nslicehit = slicelist.at(s)->NCell();
288  nslicecontplanes = slicelist.at(s)->MostContiguousPlanes(geo::kXorY);
289 
290  //Loop over slice in question here:
291  for(unsigned int hitIdx = 0; hitIdx < slicelist.at(s)->NCell(); ++hitIdx)
292  {
293  const art::Ptr<rb::CellHit>& chit = slicelist.at(s)->Cell(hitIdx);
294  const int planeNum = chit->Plane();
295  cellsfromedge = std::min((unsigned int)chit->Cell(), geom->Plane(planeNum)->Ncells() - 1 - chit->Cell());
296 
297  if(cellsfromedge < mincellsfromedge)
298  mincellsfromedge = cellsfromedge;
299 
300  } //END OF SLICE LOOP HERE
301 
302  // And ends here (lo and behold)
303 
304  //Pulling out all info from the energy object and plotting it - these plots are without any cuts
305  CalCCEnergy = sliceEnergy->CalCCE();
306  QEEnergy = sliceEnergy->TrkQEE();
307  NonQEEnergy = sliceEnergy->TrkNonQEE();
308  CCEnergy = sliceEnergy->TrkCCE();
309  AngleQEEnergy = sliceEnergy->AngleQEE();
310  AngleQEError = sliceEnergy->AngleQEError();
311  HadCalEnergy = sliceEnergy->HadCalE();
312  HadTrkEnergy = sliceEnergy->HadTrkE();
313  ndTrkLenAct = sliceEnergy->NDTrkLenAct();
314  ndTrkLenCat = sliceEnergy->NDTrkLenCat();
315  ndTrkCalAct = sliceEnergy->NDTrkCalAct();
316  ndTrkCalTran = sliceEnergy->NDTrkCalTran();
317  ndTrkCalCat = sliceEnergy->NDTrkCalCat();
318  ndHadCalAct = sliceEnergy->NDHadCalAct();
319  ndHadCalTran = sliceEnergy->NDHadCalTran();
320  ndHadCalCat = sliceEnergy->NDHadCalCat();
321  ndHadTrkAct = sliceEnergy->NDHadTrkAct();
322  ndHadTrkTran = sliceEnergy->NDHadTrkTran();
323  ndHadTrkCat = sliceEnergy->NDHadTrkCat();
324  ndTrkTranX = sliceEnergy->NDTrkTranX();
325  ndTrkTranY = sliceEnergy->NDTrkTranY();
326 
327  fCalCCEnergy->Fill(CalCCEnergy);
328  fTrkQEEnergy->Fill(QEEnergy);
342 
343  //Now start trying to make super pure samples with truth to fill resolution plots for neutrino energy
345  sliceHits = slicelist[s]->AllCells();
346 
347  std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(sliceHits,allhits);
348  if (funcReturn.empty()) continue;
349 
350  NuPurity = funcReturn[0].purity;
351  NuEfficiency = funcReturn[0].efficiency;
352 
353  //Getting CC neutrinos
354  const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();
355  isCC = (test_neutrino.CCNC()==simb::kCC);
356  PDGCode = test_neutrino.Nu().PdgCode();
357 
358  //Checking if QE or not
359  TrueQEInt = false;
360  if ((test_neutrino.InteractionType()==simb::kCCQE)) TrueQEInt = true;
361 
362  //Determining information about the muon - does it have enough hits to make this a decent interaction?
363  std::vector<const sim::Particle*> particleList = bt->MCTruthToParticles(funcReturn[0].neutrinoInt);
364  std::set<int> neutrinoTrackIDs;
365  for (unsigned int i = 0; i < particleList.size(); ++i){
366  neutrinoTrackIDs.insert(particleList[i]->TrackId());
367  }
368 
369  std::set<int>::iterator itr = neutrinoTrackIDs.begin();
370 
371  goodMuon = -1;
372  bestTrack = -1;
373  trackLength = 0.0;
374  catcherE = 0.0;
375  TrueTrackLength = 0.0;
376  TrueMuonE = 0.0;
377 
378  while( itr != neutrinoTrackIDs.end() ){
379 
380  int id_int = *itr;
381  const sim::Particle* particle = bt->TrackIDToParticle(id_int);
382  int PDG = particle->PdgCode();
383 
384  if ((PDG != 13)&&(PDG != -13)){
385  ++itr;
386  continue;
387  }//End of loop over non-muons
388 
389  //Trying to find main muon for the event
390  const sim::Particle* mother = bt->TrackIDToMotherParticle(id_int);
391  if (mother->TrackId() != id_int){
392  ++itr;
393  continue;
394  }//End of loop over muons who don't come from the neutrino
395 
396  bool passCuts = bt->PassMuonCuts(id_int,sliceHits);
397 
398  if (!passCuts){
399  ++itr;
400  continue;
401  }//End of loop over muons who don't pass cuts
402 
403  goodMuon = id_int;
405  TrueMuonE = particle->E();
406 
407  double closeToCentre = 50.0;
408 
409  std::vector<sim::FLSHit> flshits = bt->ParticleToFLSHit(particle->TrackId());
410 
411  for(unsigned int q = 0; q < flshits.size(); ++q){
412  if((flshits[q].GetPlaneID() == 191) && (std::abs(flshits[q].GetPDG())==13)){
413  if(fabs(flshits[q].GetEntryY()) < closeToCentre){
414  catcherE = flshits[q].GetEntryEnergy() + 0.105658;
415  closeToCentre = fabs(flshits[q].GetEntryY());
416  }
417  if(fabs(flshits[q].GetExitY()) < closeToCentre){
418  catcherE = flshits[q].GetExitEnergy() + 0.105658;
419  closeToCentre = fabs(flshits[q].GetExitY());
420  }
421  } // End of loop over plane 191
422  } // End of loop over flshits
423 
424  break;
425  }//End of while loop over neutrino particles
426 
427  if (goodMuon == -1) continue;
428 
429  //Now we know we have an interaction of the correct type, passes muon cuts. We now loop over tracks.
430 
431 
432  // Get all of the tracks in the slice
433  const std::vector< art::Ptr<rb::Track> > sliceTracks = sliceToTracks.at(s);
434  if(sliceTracks.empty()) continue;
435 
436  // Get all of the cosmic tracks in the slice
437  const std::vector< art::Ptr<rb::Track> > sliceCosmicTracks = sliceToCosmicTracks.at(s);
438  nCosmic = sliceCosmicTracks.size();
439 
440  bestPID = 0.0;
441  RecoMuonE = -5.0;
442 
443  art::FindOneP<remid::ReMId> trackRemid(sliceTracks, e, fPIDModuleLabel);
444  art::FindOneP<rb::Energy> trackEnergy(sliceTracks, e, fEnergyModuleLabel);
445 
446  // Get Reco values for the best track
448 
449  if(bestTrack == -1) continue;
450  if(bestTrack == 999) continue;
451 
452  trackLength = sliceTracks[bestTrack]->TotalLength();
453  bestPID = trackRemid.at(bestTrack)->Value();
454 
455  if(trackEnergy.isValid()){
456  RecoMuonE = trackEnergy.at(bestTrack)->E();
457  }
458 
459  //Checking if reco muon track matches truth muon well enough
460  trackMatchesMuon = false;
461  std::set<int> muonTrackID;
462  muonTrackID.insert(goodMuon);
463  art::PtrVector< rb::CellHit > bestTrackHits;
464  bestTrackHits = sliceTracks[bestTrack]->AllCells();
465  trackMuonEffHits = bt->HitCollectionEfficiency(muonTrackID, bestTrackHits, allhits, geo::kXorY, 0, false);
466  trackMuonPurHits = bt->HitCollectionPurity(muonTrackID, bestTrackHits, 0, 0, false);
467 
468  if ((trackMuonEffHits > 0.8)&&(trackMuonPurHits > 0.8)) trackMatchesMuon = true;
469 
470  //Determining if track is contained enough
471  contained = false;
472 
473  firstplane = -1;
474  lastplane = -1;
475  nplanestofront = -1;
476  nplanestoback = -1;
477 
478  if (sliceTracks[bestTrack]->NCell()>0){
479 
480  firstplane = sliceTracks[bestTrack]->MinPlane();
481  lastplane = sliceTracks[bestTrack]->MaxPlane();
482  if(fFD)
483  {
484  nplanestofront = firstplane - (dibfirst-1)*64;
485  nplanestoback = (diblast*64) - 1 - lastplane;
486  }
487  if(fND)
488  {
490  nplanestoback = 214 - lastplane;
491  }
492  start_x_reco = sliceTracks[bestTrack]->Start().X();
493  start_y_reco = sliceTracks[bestTrack]->Start().Y();
494  start_z_reco = sliceTracks[bestTrack]->Start().Z();
495 
496  end_x_reco = sliceTracks[bestTrack]->Stop().X();
497  end_y_reco = sliceTracks[bestTrack]->Stop().Y();
498  end_z_reco = sliceTracks[bestTrack]->Stop().Z();
499 
500  if (fFD){
501  if(mincellsfromedge > 1){
502  if (sliceCosrej->CosFwdCell() > 0 && sliceCosrej->CosBakCell() > 0){
503  if (sliceCosrej->KalFwdCell() > 10 && sliceCosrej->KalBakCell() > 10){
504  if (nplanestofront > 1 && nplanestoback > 1){
505  contained = true;
506  }
507  }
508  }
509  }
510  }//End of FD containment logic chain
511 
512  if (fND){
513  if(mincellsfromedge > 1){
514  if(sliceCosrej->KalFwdCellND() > 4 && sliceCosrej->KalBakCellND() > 8){
515  if(firstplane > 1 && lastplane < 212){
516  if(start_z_reco < 1150 && (end_z_reco < 1275 || sliceCosrej->KalYPosAtTrans() < 55)){
517  if( (ndHadCalTran + ndHadCalCat) < 0.03){
518  contained = true;
519  }
520  }
521  }
522  }
523  }
524  }//End of ND containemnt logic chain
525  }//End of loop over best tracks to find start/stop
526 
527  //Filling remaining histograms with relevant super pure populations
528  TrueNuE = test_neutrino.Nu().E();
529 
532 
533  if (TrueQEInt){
538  }//End of pure QE slices
539  else{
542  }//End of pure NonQE slices
543 
544  nKalman = sliceTracks.size();
545 
546  fNumuETree->Fill(); // Fill the tree for the pure sample
547  }//End of loop over slices
548  return;
549  }//End of analyze
550 
551  //-------------------------------------------------------------------
553  {
554  fHitModuleLabel = p.get< std::string >("HitModuleLabel" );
555  fTrackModuleLabel = p.get< std::string >("TrackModuleLabel" );
556  fCosmicTrackModuleLabel = p.get< std::string >("CosmicTrackModuleLabel");
557  fSlicerModuleLabel = p.get< std::string >("SlicerModuleLabel" );
558  fEnergyModuleLabel = p.get< std::string >("EnergyModuleLabel" );
559  fPIDModuleLabel = p.get< std::string >("PIDModuleLabel" );
560  fCosRejLabel = p.get< std::string >("CosRejLabel" );
561  fQepidLabel = p.get< std::string >("QepidLabel" );
562  }
563 
564  //-------------------------------------------------------------------
566  {
568 
569  fNumuETree = tfs->make<TTree>("NumuETree","Numu analysis Energy tree"); // Output tree
570  fNumuETree->Branch("CalCCEnergy",&CalCCEnergy);
571  fNumuETree->Branch("QEEnergy",&QEEnergy);
572  fNumuETree->Branch("NonQEEnergy",&NonQEEnergy);
573  fNumuETree->Branch("CCEnergy",&CCEnergy);
574  fNumuETree->Branch("AngleQEEnergy",&AngleQEEnergy);
575  fNumuETree->Branch("AngleQEError",&AngleQEError);
576  fNumuETree->Branch("HadCalEnergy",&HadCalEnergy);
577  fNumuETree->Branch("HadTrkEnergy",&HadTrkEnergy);
578  fNumuETree->Branch("ndTrkLenAct",&ndTrkLenAct);
579  fNumuETree->Branch("ndTrkLenCat",&ndTrkLenCat);
580  fNumuETree->Branch("ndTrkCalAct",&ndTrkCalAct);
581  fNumuETree->Branch("ndTrkCalTran",&ndTrkCalTran);
582  fNumuETree->Branch("ndTrkCalCat",&ndTrkCalCat);
583  fNumuETree->Branch("ndHadCalAct",&ndHadCalAct);
584  fNumuETree->Branch("ndHadCalTran",&ndHadCalTran);
585  fNumuETree->Branch("ndHadCalCat",&ndHadCalCat);
586  fNumuETree->Branch("ndHadTrkAct",&ndHadTrkAct);
587  fNumuETree->Branch("ndHadTrkTran",&ndHadTrkTran);
588  fNumuETree->Branch("ndHadTrkCat",&ndHadTrkCat);
589  fNumuETree->Branch("ndTrkTranX",&ndTrkTranX);
590  fNumuETree->Branch("ndTrkTranY",&ndTrkTranY);
591  fNumuETree->Branch("ND",&fND);
592  fNumuETree->Branch("FD",&fFD);
593  // Booleans
594  fNumuETree->Branch("isCC",&isCC);
595  fNumuETree->Branch("TrueQEInt",&TrueQEInt);
596  fNumuETree->Branch("trackMatchesMuon",&trackMatchesMuon);
597  fNumuETree->Branch("contained",&contained);
598  // Integers
599  fNumuETree->Branch("Run",&Run);
600  fNumuETree->Branch("SubRun",&SubRun);
601  fNumuETree->Branch("Evt",&Evt);
602  fNumuETree->Branch("SubEvt",&SubEvt);
603  fNumuETree->Branch("PDGCode",&PDGCode);
604  fNumuETree->Branch("goodMuon",&goodMuon);
605  fNumuETree->Branch("bestTrack",&bestTrack);
606  fNumuETree->Branch("nKalman",&nKalman);
607  fNumuETree->Branch("nCosmic",&nCosmic);
608  // Containment integers
609  fNumuETree->Branch("firstplane",&firstplane);
610  fNumuETree->Branch("lastplane",&lastplane);
611  fNumuETree->Branch("nplanestofront",&nplanestofront);
612  fNumuETree->Branch("nplanestoback",&nplanestoback);
613  fNumuETree->Branch("mincellsfromedge",&mincellsfromedge);
614  fNumuETree->Branch("nslicehit",&nslicehit);
615  fNumuETree->Branch("nslicecontplanes",&nslicecontplanes);
616  // Doubles
617  fNumuETree->Branch("NuPurity",&NuPurity);
618  fNumuETree->Branch("NuEfficiency",&NuEfficiency);
619  fNumuETree->Branch("trackLength",&trackLength);
620  fNumuETree->Branch("catcherE",&catcherE);
621  fNumuETree->Branch("bestPID",&bestPID);
622  fNumuETree->Branch("qepidPID",&qepidPID);
623  fNumuETree->Branch("qepidntrk",&qepidntrk);
624  fNumuETree->Branch("cosrejanglekal",&cosrejanglekal);
625  fNumuETree->Branch("cosrejnumucontpid",&cosrejnumucontpid);
626  fNumuETree->Branch("trackMuonEffHits",&trackMuonEffHits);
627  fNumuETree->Branch("trackMuonPurHits",&trackMuonPurHits);
628  fNumuETree->Branch("start_x_reco",&start_x_reco);
629  fNumuETree->Branch("start_y_reco",&start_y_reco);
630  fNumuETree->Branch("start_z_reco",&start_z_reco);
631  fNumuETree->Branch("end_x_reco",&end_x_reco);
632  fNumuETree->Branch("end_y_reco",&end_y_reco);
633  fNumuETree->Branch("end_z_reco",&end_z_reco);
634  fNumuETree->Branch("TrueNuE",&TrueNuE);
635  fNumuETree->Branch("TrueMuonE",&TrueMuonE);
636  fNumuETree->Branch("RecoMuonE",&RecoMuonE);
637  fNumuETree->Branch("TrueTrackLength",&TrueTrackLength);
638 
639  // And finally the histograms
640  fCalCCEnergy = tfs->make<TH1D>("CalCCEnergy", ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
641  fTrkQEEnergy = tfs->make<TH1D>("TrkQEEnergy", ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
642  fTrkNonQEEnergy = tfs->make<TH1D>("TrkNonQEEnergy", ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
643  fAngleQEEnergy = tfs->make<TH1D>("AngleQEEnergy", ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
644 
645  fAngleQEError = tfs->make<TH1D>("AngleQEError", ";Error on Neutrino Energy;Slices",200, -0.1, 5.0);
646 
647  fCalCCEnergyDiff = tfs->make<TH1D>("CalCCEnergyDiff", ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
648  fTrkQEEnergyDiff = tfs->make<TH1D>("TrkQEEnergyDiff", ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
649  fTrkNonQEEnergyDiff = tfs->make<TH1D>("TrkNonQEEnergyDiff", ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
650  fAngleQEEnergyDiff = tfs->make<TH1D>("AngleQEEnergyDiff", ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
651 
652  fCalCCEnergy2D = tfs->make<TH2D>("CalCCEnergy2D", ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
653  fTrkQEEnergy2D = tfs->make<TH2D>("TrkQEEnergy2D", ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
654  fTrkNonQEEnergy2D = tfs->make<TH2D>("TrkNonQEEnergy2D", ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
655  fAngleQEEnergy2D = tfs->make<TH2D>("AngleQEEnergy2D", ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
656 
657  fHadronCalEnergy = tfs->make<TH1D>("HadronCalEnergy", ";Reco Hadron Energy (GeV);Slices",200, -0.1, 10);
658  fHadronTrkEnergy = tfs->make<TH1D>("HadronTrkEnergy", ";Track Hadron Energy (GeV);Slices",200, -0.1, 5);
659 
660  fNDTrkLenActive = tfs->make<TH1D>("NDTrkLenActive", ";Track Length (cm) in Active Region;Slices",200, 0., 1400);
661  fNDTrkLenCatcher = tfs->make<TH1D>("NDTrkLenCatcher", ";Track Length (cm) in Catcher Region;Slices",200, 0., 800);
662  fNDTrkCalActive = tfs->make<TH1D>("NDTrkCalActive", ";Track Energy (GeV) in Active Region;Slices",200, 0., 10);
663  fNDTrkCalTransition = tfs->make<TH1D>("NDTrkCalTransition", ";Track Energy (GeV) in Tran. Plane;Slices",200, 0., 0.5);
664  fNDTrkCalCatcher = tfs->make<TH1D>("NDTrkCalCatcher", ";Track Energy (GeV) in Catcher Region;Slices",200, 0., 2);
665  fNDHadCalActive = tfs->make<TH1D>("NDHadCalActive", ";Hadron Energy (GeV) in Active Region;Slices",200, 0., 10);
666  fNDHadCalTransition = tfs->make<TH1D>("NDHadCalTransition", ";Hadron Energy (GeV) in Tran. Plane;Slices",200, 0., 0.5);
667  fNDHadCalCatcher = tfs->make<TH1D>("NDHadCalCatcher", ";Hadron Energy (GeV) in Catcher Region;Slices",200, 0., 2);
668 
669  }//End of begin Run
670 
671  //-------------------------------------------------------------------
673  {
675  double detHalfWidth = geom->DetHalfWidth();
676  double detHalfHeight = geom->DetHalfHeight();
677  double detLength = geom->DetLength();
678 
679  fMaxTrkLen = util::pythag(detHalfWidth*2, detHalfHeight*2, detLength);
680 
681  fND = false;
682  fFD = false;
683 
684  novadaq::cnv::DetId detID = geom->DetId();
685  if (detID == novadaq::cnv::kNEARDET) fND = true;
686  if (detID == novadaq::cnv::kFARDET) fFD = true;
687 
688  }//End of beginRun
689 
691  {
692  const unsigned int N = p->NumberTrajectoryPoints();
693 
694  if(N < 2){
695  mf::LogWarning("RecoCheckAna") << "Particle has no trajectory points";
696  return 0;
697  }
698 
700 
701  double dist = 0;
702 
703  for(unsigned int n = 0; n < N-1; ++n){
704  TVector3 p0 = p->Position(n).Vect();
705  TVector3 p1 = p->Position(n+1).Vect();
706  if(p0.Z() < 0 || p1.Z() < 0 ||
707  p0.Z() > geom->DetLength() || p1.Z() > geom->DetLength() ||
708  fabs(p0.X()) > geom->DetHalfWidth() || fabs(p1.X()) > geom->DetHalfWidth() ||
709  fabs(p0.Y()) > geom->DetHalfHeight() || fabs(p1.Y()) > geom->DetHalfHeight()) continue;
710 
711  if(view == geo::kX){
712  p0.SetY(0);
713  p1.SetY(0);
714  }
715  if(view == geo::kY){
716  p0.SetX(0);
717  p1.SetX(0);
718  }
719  dist += (p1-p0).Mag();
720  }
721 
722  return dist;
723 } // end of TotalLengthInDetector
724 
725  //----------------------------------------------------------
727 
728 }//End of namespace
double E(const int i=0) const
Definition: MCParticle.h:232
std::string fPIDModuleLabel
label for module writing pid objects
std::string fCosmicTrackModuleLabel
label for module making cosmic tracks
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
SubRunNumber_t subRun() const
Definition: Event.h:72
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
int CCNC() const
Definition: MCNeutrino.h:148
TH2D * fCalCCEnergy2D
histogram of calorimetric neutrino E and true neutrino E (for pure pop.)
Energy estimators for CC events.
Definition: FillEnergies.h:7
TH1D * fNDHadCalCatcher
histogram of ND hadron calorimetric energy in the muon catcher region
TTree * fNumuETree
n-tuple used for further analyses (pure population)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH1D * fHadronCalEnergy
histogram of hadronic calorimetric energy NOT on the muon track
bool fND
is detector ND?
const sim::Particle * TrackIDToMotherParticle(int const &id) const
bool contained
does it pass numu containment cuts?
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
TH1D * fNDTrkLenActive
histogram of ND muon track length in the active region
std::string fEnergyModuleLabel
label for module writing energy objects
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
virtual void analyze(art::Event const &e)
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
int qepidntrk
Number of tracks used by QePId.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::string fTrackModuleLabel
label for module making the tracks
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
art::ServiceHandle< nova::dbi::RunHistoryService > rh
TH1D * fHadronTrkEnergy
histogram of hadronic calorimetric energy on the muon track
double DetLength() const
double TrueMuonE
True energy of the muon.
TH1D * fNDTrkCalCatcher
histogram of ND muon calorimetric energy in the muon catcher region
int KalBakCellND() const
Definition: CosRejObj.cxx:670
TH1D * fCalCCEnergyDiff
histogram of the diff. between calorimetric neutrino E and true neutrino E (for pure pop...
bool trackMatchesMuon
track is muon-like
double cosrejnumucontpid
Cosmic rejection PID for contained events from CosRej.
charged current quasi-elastic
Definition: MCNeutrino.h:96
const PlaneGeo * Plane(unsigned int i) const
TH1D * fAngleQEError
histogram of error on the neutrino energy for QE population using angle formula
double RecoMuonE
Reco energy of the muon.
DEFINE_ART_MODULE(TestTMapFile)
TH1D * fNDTrkCalTransition
histogram of ND muon calorimetric energy in the transition region
int KalFwdCell() const
Definition: CosRejObj.cxx:652
double bestPID
Highest RemID.
TH1D * fTrkNonQEEnergy
histogram of muon track length neutrino energy for NonQE super pure population
float abs(float number)
Definition: d0nt_math.hpp:39
double Value() const
Definition: PID.h:22
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::string fHitModuleLabel
lable for module making the hits
Definition: Run.h:31
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
int TrackId() const
Definition: MCParticle.h:209
double dist
Definition: runWimpSim.h:113
double catcherE
Energy before entering the muon catcher (ND)
std::string fSlicerModuleLabel
label for module making the slices
NumuEAna(fhicl::ParameterSet const &p)
unsigned short Cell() const
Definition: CellHit.h:40
TH1D * fNDHadCalActive
histogram of ND hadron calorimetric energy in the active region
const XML_Char * s
Definition: expat.h:262
int InteractionType() const
Definition: MCNeutrino.h:150
TH1D * fNDTrkLenCatcher
histogram of ND muon track length in the muon catcher region
TString part[npart]
Definition: Style.C:32
bool fFD
is detector FD?
Far Detector at Ash River, MN.
bool TrueQEInt
is the interaction QE?
TH2D * fTrkQEEnergy2D
histogram of muon track length neutrino E and true neutrino E (for pure QE pop.)
double cosrejanglekal
Cosine of angle of best ReMId Kalman Track from CosRej.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double TrueNuE
True energy of the neutrino.
TH1D * fTrkQEEnergy
histogram of muon track length neutrino energy for QE super pure population
std::string fQepidLabel
label for module writing qe pid objects
double TrueTrackLength
True tracklength of the muon.
virtual void reconfigure(fhicl::ParameterSet const &p)
unsigned int mincellsfromedge
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual void beginJob()
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
float ConCosPID() const
Definition: CosRejObj.cxx:465
int KalBakCell() const
Definition: CosRejObj.cxx:664
Near Detector in the NuMI cavern.
TH1D * fAngleQEEnergy
histogram of neutrino energy for QE population using angle formula
int CosBakCell() const
Definition: CosRejObj.cxx:640
int CosFwdCell() const
Definition: CosRejObj.cxx:628
TH1D * fTrkQEEnergyDiff
histogram of the diff. between muon track length neutrino E and true neutrino E (for pure QE pop...
EventNumber_t event() const
Definition: Event.h:67
reference at(size_type n)
Definition: PtrVector.h:365
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
size_type size() const
Definition: PtrVector.h:308
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
TH2D * fTrkNonQEEnergy2D
histogram of muon track length neutrino E and true neutrino E (for pure NonQE pop.)
TH1D * fNDTrkCalActive
histogram of ND muon calorimetric energy in the active region
int Ntrk() const
Definition: QePId.cxx:75
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double NuPurity
Neutrino selection purity.
bool isCC
is the interaction CC?
double DetHalfWidth() const
double qepidPID
PID value of output kNN from QePId.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
double fMaxTrkLen
length of detector from corner to corner
Definition: event.h:1
TH1D * fNDHadCalTransition
histogram of ND hadron calorimetric energy in the transition region
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
float Mag() const
double trackLength
Reconstructed tracklength.
TRandom3 r(0)
int KalFwdCellND() const
Definition: CosRejObj.cxx:658
particleList
Definition: demo.py:41
TH1D * fAngleQEEnergyDiff
histogram of the diff. between neutrino E using angle formula and true neutrino E (for pure QE pop...
bool PassMuonCuts(int trackID, art::PtrVector< rb::CellHit > const &hits) const
Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function) cont...
T min(const caf::Proxy< T > &a, T b)
double NuEfficiency
Neutrino selection efficiency.
int GoodDiBlockMask(int subrun=-1, bool reload=false)
std::string fCosRejLabel
label for module writing cosmic rejection
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
Event generator information.
Definition: MCNeutrino.h:18
double TotalLengthInDetector(const sim::Particle *part, geo::View_t view)
virtual void beginRun(art::Run const &r)
TH1D * fTrkNonQEEnergyDiff
histogram of the diff. between muon track length neutrino E and true neutrino E (for pure NonQE pop...
Encapsulate the geometry of one entire detector (near, far, ndos)
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249
TH1D * fCalCCEnergy
histogram of calorimetric neutrino energy for general CC population
TH2D * fAngleQEEnergy2D
histogram of neutrino E using angle formula and true neutrino E (for pure QE pop.) ...
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
float AngleKal() const
Definition: CosRejObj.cxx:453