QeFinderVal_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////
2 /// \brief A validation module for the QEEventFinder package
3 /// \author Nicholas Raddatz - raddatz@physics.umn.edu
4 /// \date February 2013
5 ////////////////////////////////////////////////////////////////////
6 
7 // ROOT includes
8 #include "TH1.h"
9 
10 // Framework includes
14 #include "fhiclcpp/ParameterSet.h"
16 
17 //NOvA includes
18 #include "QEEventFinder/QePId.h"
19 #include "ReMId/ReMId.h"
20 #include "RecoBase/Track.h"
21 #include "NumuEnergy/NumuE.h"
22 #include "MCCheater/BackTracker.h"
23 #include "Utilities/AssociationUtil.h"
24 
25 
26 namespace qeef{
27  class QeFinderVal : public art::EDAnalyzer{
28  public:
29  explicit QeFinderVal(fhicl::ParameterSet const& pset);
30  virtual ~QeFinderVal();
31 
32  void beginJob();
33  void analyze(art::Event const& evt);
34 
35  private:
36 
42 
43  // make plots of variables for all slices
45 
46  // plots of qepid variables for 1 3d track slices
48 
49  // plots of qepid variables for 2 track and >=1 3d track slices
51 
52  // plots of qepid variables for slices that don't fall into the above to categories
54 
55  // plots of 1 and 2 track pid slcies that were true qe events
58 
59  // plots of 1 and 2 track pid slcies that were true nonqe events
62 
63  };
64 }
65 
66 namespace qeef{
67 
68  //.......................................................................
70  : EDAnalyzer(pset)
71  {
72  fSliceLabel = (pset.get< std::string >("SliceLabel"));
73  fTrackLabel = (pset.get< std::string >("TrackLabel"));
74  fPIDLabel = (pset.get< std::string >("PIDLabel"));
75  fEnergyLabel = (pset.get< std::string >("EnergyLabel"));
76  fQePIdLabel = (pset.get< std::string >("QePIDLabel"));
77 
78  }
79 
80  //.......................................................................
82  {
83  }
84 
85  //.......................................................................
87  {
88 
90 
91  // qepid info for all slices
92  fNtrkAll = tfs->make<TH1F>("NtrkAll",";Ntrk;Slices",5,-1.5,3.5);
93  fModeAll = tfs->make<TH1F>("ModeAll",";Mode;Slices",3,-1.5,1.5);
94  fOffTrkEAll = tfs->make<TH1F>("OffTrkEAll",";OffE;Slices",100,-0.05,1.05);
95  fEdiffAll = tfs->make<TH1F>("EdiffAll",";Ediff;Slices",100,-1.0,1.0);
96  fEdiffZAll = tfs->make<TH1F>("EdiffZAll",";EdiffZ;Slices",100,-15,5);
97  fDedxAll = tfs->make<TH1F>("DedxAll",";Dedx;Slices",100,-1,6);
98  fPdgAll = tfs->make<TH1F>("PdgAll",";Pdg;Slices",16,-0.5,15.5);
99  fPidAll = tfs->make<TH1F>("PidAll",";PID;Slices",100,-5.5,1.5);
100 
101 
102  // qepid info for slices with 1 and only 1 3d track
103  fNtrk1Trk = tfs->make<TH1F>("Ntrk1Trk",";Ntrk;Slices",5,-1.5,3.5);
104  fMode1Trk = tfs->make<TH1F>("Mode1Trk",";Mode;Slices",3,-1.5,1.5);
105  fOffTrkE1Trk = tfs->make<TH1F>("OffTrkE1Trk",";OffE;Slices",100,-0.05,1.05);
106  fEdiff1Trk = tfs->make<TH1F>("Ediff1Trk",";Ediff;Slices",100,-1.0,1.0);
107  fEdiffZ1Trk = tfs->make<TH1F>("EdiffZ1Trk",";EdiffZ;Slices",100,-15,5);
108  fDedx1Trk = tfs->make<TH1F>("Dedx1Trk",";Dedx;Slices",100,-1,6);
109  fPdg1Trk = tfs->make<TH1F>("Pdg1Trk",";Pdg;Slices",16,-0.5,15.5);
110  fPid1Trk = tfs->make<TH1F>("Pid1Trk",";PID;Slices",100,-5.5,1.5);
111 
112  // qepid info for slices with 2 total tracks, with at least one 3d track
113  fNtrk2Trk = tfs->make<TH1F>("Ntrk2Trk",";Ntrk;Slices",5,-1.5,3.5);
114  fMode2Trk = tfs->make<TH1F>("Mode2Trk",";Mode;Slices",3,-1.5,1.5);
115  fOffTrkE2Trk = tfs->make<TH1F>("OffTrkE2Trk",";OffE;Slices",100,-0.05,1.05);
116  fEdiff2Trk = tfs->make<TH1F>("Ediff2Trk",";Ediff;Slices",100,-1.0,1.0);
117  fEdiffZ2Trk = tfs->make<TH1F>("EdiffZ2Trk",";EdiffZ;Slices",100,-15,5);
118  fDedx2Trk = tfs->make<TH1F>("Dedx2Trk",";Dedx;Slices",100,-1,6);
119  fPdg2Trk = tfs->make<TH1F>("Pdg2Trk",";Pdg;Slices",16,-0.5,15.5);
120  fPid2Trk = tfs->make<TH1F>("Pid2Trk",";PID;Slices",100,-5.5,1.5);
121 
122  // qepid info slices that shouldn't have a qepid pid between 0 and 1
123  fNtrk0Trk = tfs->make<TH1F>("Ntrk0Trk",";Ntrk;Slices",5,-1.5,3.5);
124  fMode0Trk = tfs->make<TH1F>("Mode0Trk",";Mode;Slices",3,-1.5,1.5);
125  fOffTrkE0Trk = tfs->make<TH1F>("OffTrkE0Trk",";OffE;Slices",100,-0.05,1.05);
126  fEdiff0Trk = tfs->make<TH1F>("Ediff0Trk",";Ediff;Slices",100,-1.0,1.0);
127  fEdiffZ0Trk = tfs->make<TH1F>("EdiffZ0Trk",";EdiffZ;Slices",100,-15,5);
128  fDedx0Trk = tfs->make<TH1F>("Dedx0Trk",";Dedx;Slices",100,-1,6);
129  fPdg0Trk = tfs->make<TH1F>("Pdg0Trk",";Pdg;Slices",16,-0.5,15.5);
130  fPid0Trk = tfs->make<TH1F>("Pid0Trk",";PID;Slices",100,-5.5,1.5);
131 
132  // make plots for 1 and 2 track true qe events
133  fNtrk1TrkTrueQE = tfs->make<TH1F>("Ntrk1TrkTrueQE",";Ntrk;Slices",5,-1.5,3.5);
134  fMode1TrkTrueQE = tfs->make<TH1F>("Mode1TrkTrueQE",";Mode;Slices",3,-1.5,1.5);
135  fOffTrkE1TrkTrueQE = tfs->make<TH1F>("OffTrkE1TrkTrueQE",";OffE;Slices",100,-0.05,1.05);
136  fEdiff1TrkTrueQE = tfs->make<TH1F>("Ediff1TrkTrueQE",";Ediff;Slices",100,-1.0,1.0);
137  fEdiffZ1TrkTrueQE = tfs->make<TH1F>("EdiffZ1TrkTrueQE",";EdiffZ;Slices",100,-15,5);
138  fDedx1TrkTrueQE = tfs->make<TH1F>("Dedx1TrkTrueQE",";Dedx;Slices",100,-1,6);
139  fPdg1TrkTrueQE = tfs->make<TH1F>("Pdg1TrkTrueQE",";Pdg;Slices",16,-0.5,15.5);
140  fPid1TrkTrueQE = tfs->make<TH1F>("Pid1TrkTrueQE",";PID;Slices",100,-5.5,1.5);
141  fNtrk2TrkTrueQE = tfs->make<TH1F>("Ntrk2TrkTrueQE",";Ntrk;Slices",5,-1.5,3.5);
142  fMode2TrkTrueQE = tfs->make<TH1F>("Mode2TrkTrueQE",";Mode;Slices",3,-1.5,1.5);
143  fOffTrkE2TrkTrueQE = tfs->make<TH1F>("OffTrkE2TrkTrueQE",";OffE;Slices",100,-0.05,1.05);
144  fEdiff2TrkTrueQE = tfs->make<TH1F>("Ediff2TrkTrueQE",";Ediff;Slices",100,-1.0,1.0);
145  fEdiffZ2TrkTrueQE = tfs->make<TH1F>("EdiffZ2TrkTrueQE",";EdiffZ;Slices",100,-15,5);
146  fDedx2TrkTrueQE = tfs->make<TH1F>("Dedx2TrkTrueQE",";Dedx;Slices",100,-1,6);
147  fPdg2TrkTrueQE = tfs->make<TH1F>("Pdg2TrkTrueQE",";Pdg;Slices",16,-0.5,15.5);
148  fPid2TrkTrueQE = tfs->make<TH1F>("Pid2TrkTrueQE",";PID;Slices",100,-5.5,1.5);
149 
150  // make plots for 1 and 2 track true nonqe events
151  fNtrk1TrkTrueNonQE = tfs->make<TH1F>("Ntrk1TrkTrueNonQE",";Ntrk;Slices",5,-1.5,3.5);
152  fMode1TrkTrueNonQE = tfs->make<TH1F>("Mode1TrkTrueNonQE",";Mode;Slices",3,-1.5,1.5);
153  fOffTrkE1TrkTrueNonQE = tfs->make<TH1F>("OffTrkE1TrkTrueNonQE",";OffE;Slices",100,-0.05,1.05);
154  fEdiff1TrkTrueNonQE = tfs->make<TH1F>("Ediff1TrkTrueNonQE",";Ediff;Slices",100,-1.0,1.0);
155  fEdiffZ1TrkTrueNonQE = tfs->make<TH1F>("EdiffZ1TrkTrueNonQE",";EdiffZ;Slices",100,-15,5);
156  fDedx1TrkTrueNonQE = tfs->make<TH1F>("Dedx1TrkTrueNonQE",";Dedx;Slices",100,-1,6);
157  fPdg1TrkTrueNonQE = tfs->make<TH1F>("Pdg1TrkTrueNonQE",";Pdg;Slices",16,-0.5,15.5);
158  fPid1TrkTrueNonQE = tfs->make<TH1F>("Pid1TrkTrueNonQE",";PID;Slices",100,-5.5,1.5);
159  fNtrk2TrkTrueNonQE = tfs->make<TH1F>("Ntrk2TrkTrueNonQE",";Ntrk;Slices",5,-1.5,3.5);
160  fMode2TrkTrueNonQE = tfs->make<TH1F>("Mode2TrkTrueNonQE",";Mode;Slices",3,-1.5,1.5);
161  fOffTrkE2TrkTrueNonQE = tfs->make<TH1F>("OffTrkE2TrkTrueNonQE",";OffE;Slices",100,-0.05,1.05);
162  fEdiff2TrkTrueNonQE = tfs->make<TH1F>("Ediff2TrkTrueNonQE",";Ediff;Slices",100,-1.0,1.0);
163  fEdiffZ2TrkTrueNonQE = tfs->make<TH1F>("EdiffZ2TrkTrueNonQE",";EdiffZ;Slices",100,-15,5);
164  fDedx2TrkTrueNonQE = tfs->make<TH1F>("Dedx2TrkTrueNonQE",";Dedx;Slices",100,-1,6);
165  fPdg2TrkTrueNonQE = tfs->make<TH1F>("Pdg2TrkTrueNonQE",";Pdg;Slices",16,-0.5,15.5);
166  fPid2TrkTrueNonQE = tfs->make<TH1F>("Pid2TrkTrueNonQE",";PID;Slices",100,-5.5,1.5);
167 
168  }
169 
170  //.......................................................................
172  {
173 
174 
175  // get all of slices out of the event
177  evt.getByLabel(fSliceLabel,slices);
178 
179  // make a vector of all cellhits
181  for(size_t isl = 0; isl < slices->size(); ++isl){
182  art::Ptr<rb::Cluster> slice(slices,isl);
183  for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
184  allHits.push_back(slice->Cell(ihit));
185  }
186  }
187 
188  // set up associations between slices and qepids
189  art::FindManyP<qeef::QePId> qeefAss(slices,evt,fQePIdLabel);
190  // set up associations between slices and tracks
191  art::FindManyP<rb::Track> trkAss(slices,evt,fTrackLabel);
192  // set up associations between slices and energy
194 
195  // get a backtracker service handle
197 
198 
199  for(size_t isl = 0; isl < slices->size(); ++isl){
200 
201  // Get the slice
202  art::Ptr<rb::Cluster> slice(slices,isl);
203  // skip if it is noise
204  if(slice->IsNoise()){ continue; }
205 
206  // get the true mctruth matching for this slice
207  if(!bt->HaveTruthInfo()){
208  mf::LogWarning("QeFinderVal") << "attempting to run MC truth check on "
209  << "things without truth, bail";
210  return;
211  }
212 
213  std::vector<cheat::NeutrinoEffPur> nus = bt->SliceToNeutrinoInteractions(slice->AllCells(),allHits);
214  // Is this a true numu CCQE event
215  bool isTrueQE = false;
216  if(nus.size() > 0){
217  if((nus[0].neutrinoInt->GetNeutrino().Nu().PdgCode() == 14 ||
218  nus[0].neutrinoInt->GetNeutrino().Nu().PdgCode() == -14) &&
219  nus[0].neutrinoInt->GetNeutrino().CCNC() == simb::kCC &&
220  nus[0].neutrinoInt->GetNeutrino().Mode() == simb::kQE){ isTrueQE = true; }
221  }
222 
223  // get the pid that are associated with this slice
224  const std::vector<art::Ptr<qeef::QePId> > qepids = qeefAss.at(isl);
225  if(qepids.size() != 1){ std::cout<<"Missing a qe pid object on this slice"<<std::endl; }
226  else{
227  art::Ptr<qeef::QePId> qepid = qepids[0];
228  fNtrkAll->Fill(qepid->Ntrk());
229  fModeAll->Fill(qepid->Mode());
230  fOffTrkEAll->Fill(qepid->OffTrkE());
231  fEdiffAll->Fill(qepid->Ediff());
232  fEdiffZAll->Fill(qepid->EdiffZ());
233  fDedxAll->Fill(qepid->Dedx());
234  fPdgAll->Fill(qepid->Pdg());
235  fPidAll->Fill(qepid->Value());
236 
237  // get all of the tracks associated with this slice
238  const std::vector<art::Ptr<rb::Track> > tracks = trkAss.at(isl);
239 
240  // count how many trks we have
241  int n3d = 0;
242  int ntrks = tracks.size();
243  for(unsigned int itr = 0; itr < tracks.size(); ++itr){
244  if(tracks[itr]->View() == geo::kXorY){ ++n3d; }
245  }
246 
247  // make plots for the individual categories
248  if(n3d == 1 && ntrks ==1){
249  fNtrk1Trk->Fill(qepid->Ntrk());
250  fMode1Trk->Fill(qepid->Mode());
251  fOffTrkE1Trk->Fill(qepid->OffTrkE());
252  fEdiff1Trk->Fill(qepid->Ediff());
253  fEdiffZ1Trk->Fill(qepid->EdiffZ());
254  fDedx1Trk->Fill(qepid->Dedx());
255  fPdg1Trk->Fill(qepid->Pdg());
256  fPid1Trk->Fill(qepid->Value());
257  if(isTrueQE){
258  fNtrk1TrkTrueQE->Fill(qepid->Ntrk());
259  fMode1TrkTrueQE->Fill(qepid->Mode());
260  fOffTrkE1TrkTrueQE->Fill(qepid->OffTrkE());
261  fEdiff1TrkTrueQE->Fill(qepid->Ediff());
262  fEdiffZ1TrkTrueQE->Fill(qepid->EdiffZ());
263  fDedx1TrkTrueQE->Fill(qepid->Dedx());
264  fPdg1TrkTrueQE->Fill(qepid->Pdg());
265  fPid1TrkTrueQE->Fill(qepid->Value());
266  }
267  else{
268  fNtrk1TrkTrueNonQE->Fill(qepid->Ntrk());
269  fMode1TrkTrueNonQE->Fill(qepid->Mode());
270  fOffTrkE1TrkTrueNonQE->Fill(qepid->OffTrkE());
271  fEdiff1TrkTrueNonQE->Fill(qepid->Ediff());
272  fEdiffZ1TrkTrueNonQE->Fill(qepid->EdiffZ());
273  fDedx1TrkTrueNonQE->Fill(qepid->Dedx());
274  fPdg1TrkTrueNonQE->Fill(qepid->Pdg());
275  fPid1TrkTrueNonQE->Fill(qepid->Value());
276  }
277  }
278  else if(n3d >=1 && ntrks ==2){
279  fNtrk2Trk->Fill(qepid->Ntrk());
280  fMode2Trk->Fill(qepid->Mode());
281  fOffTrkE2Trk->Fill(qepid->OffTrkE());
282  fEdiff2Trk->Fill(qepid->Ediff());
283  fEdiffZ2Trk->Fill(qepid->EdiffZ());
284  fDedx2Trk->Fill(qepid->Dedx());
285  fPdg2Trk->Fill(qepid->Pdg());
286  fPid2Trk->Fill(qepid->Value());
287  if(isTrueQE){
288  fNtrk2TrkTrueQE->Fill(qepid->Ntrk());
289  fMode2TrkTrueQE->Fill(qepid->Mode());
290  fOffTrkE2TrkTrueQE->Fill(qepid->OffTrkE());
291  fEdiff2TrkTrueQE->Fill(qepid->Ediff());
292  fEdiffZ2TrkTrueQE->Fill(qepid->EdiffZ());
293  fDedx2TrkTrueQE->Fill(qepid->Dedx());
294  fPdg2TrkTrueQE->Fill(qepid->Pdg());
295  fPid2TrkTrueQE->Fill(qepid->Value());
296  }
297  else{
298  fNtrk2TrkTrueNonQE->Fill(qepid->Ntrk());
299  fMode2TrkTrueNonQE->Fill(qepid->Mode());
300  fOffTrkE2TrkTrueNonQE->Fill(qepid->OffTrkE());
301  fEdiff2TrkTrueNonQE->Fill(qepid->Ediff());
302  fEdiffZ2TrkTrueNonQE->Fill(qepid->EdiffZ());
303  fDedx2TrkTrueNonQE->Fill(qepid->Dedx());
304  fPdg2TrkTrueNonQE->Fill(qepid->Pdg());
305  fPid2TrkTrueNonQE->Fill(qepid->Value());
306  }
307  }
308  else{
309  fNtrk0Trk->Fill(qepid->Ntrk());
310  fMode0Trk->Fill(qepid->Mode());
311  fOffTrkE0Trk->Fill(qepid->OffTrkE());
312  fEdiff0Trk->Fill(qepid->Ediff());
313  fEdiffZ0Trk->Fill(qepid->EdiffZ());
314  fDedx0Trk->Fill(qepid->Dedx());
315  fPdg0Trk->Fill(qepid->Pdg());
316  fPid0Trk->Fill(qepid->Value());
317  }
318 
319  }
320 
321 
322  } // end loop over slices
323 
324  } // end of analyze
325 
326 } // end of qeef namespace
327 
328 namespace qeef
329 {
331 }
double EdiffZ() const
Definition: QePId.cxx:93
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
int Mode() const
Definition: QePId.cxx:69
double Dedx() const
Definition: QePId.cxx:99
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
DEFINE_ART_MODULE(TestTMapFile)
double Value() const
Definition: PID.h:22
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double OffTrkE() const
Definition: QePId.cxx:81
QeFinderVal(fhicl::ParameterSet const &pset)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
void analyze(art::Event const &evt)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Definition: View.py:1
OStream cout
Definition: OStream.cxx:6
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
int Ntrk() const
Definition: QePId.cxx:75
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
int Pdg() const
Definition: PID.h:21
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double Ediff() const
Definition: QePId.cxx:87
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
A module for finding numu CC QE interactions.
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163