ReMIdValidate_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////
2 /// \brief A module for analyzing the ReMId objects in
3 /// neutrino interactions
4 /// \author Nicholas Raddatz - raddatz@physics.umn.edu
5 /// \date September 2012
6 ////////////////////////////////////////////////////////////////////
7 
8 #include <string>
9 
10 // ROOT includes
11 #include "TVector3.h"
12 #include "TH1F.h"
13 
14 // Framework includes
17 #include "fhiclcpp/ParameterSet.h"
20 
21 // NOvA includes
22 #include "ReMId/ReMId.h"
23 #include "MCCheater/BackTracker.h"
24 #include "Calibrator/Calibrator.h"
25 #include "RecoBase/Track.h"
26 #include "Utilities/AssociationUtil.h"
27 
28 
29 namespace remid {
30  /// A module to analyze remid objects produced by the RecoMuon module
31  class ReMIdValidate : public art::EDAnalyzer {
32  public:
33 
34  explicit ReMIdValidate(fhicl::ParameterSet const& pset);
35  virtual ~ReMIdValidate();
36 
37  void beginJob();
38  void analyze(art::Event const& evt);
39  void endJob();
40  static bool SortByHits(const art::Ptr<rb::Track>& a, const art::Ptr<rb::Track>& b);
41 
42  private:
43  std::string fSliceLabel; ///<Where to find reconstructed slices
44  std::string fTrackLabel; ///<Where to find recondtructed tracks
45  std::string fPidLabel; ///<Where to find the pids
46 
47  TH1F* fpid; ///<Pid value of every track
48  TH1F* fpidnc;
49  TH1F* fpidcc;
50  TH1F* fpidcccut;
51  TH1F* fnscat;
52  TH1F* fndedx;
53  TH1F* fcosscat;
54  TH1F* favgdedx;
55  TH1F* fdedx;
56  TH1F* fpidhigh;
57  TH1F* fpidnchigh;
58  TH1F* fpidcchigh;
59  TH1F* fnscathigh;
60  TH1F* fndedxhigh;
61  TH1F* fcosscathigh;
62  TH1F* favgdedxhigh;
63  TH1F* fdedxhigh;
64  TH1F* fpidmu;
65  TH1F* fnscatmu;
66  TH1F* fndedxmu;
67  TH1F* fcosscatmu;
68  TH1F* favgdedxmu;
69  TH1F* fdedxmu;
70  TH1F* fpidmucont;
71  TH1F* fnscatmucont;
72  TH1F* fndedxmucont;
75  TH1F* fdedxmucont;
76 
77  bool IsContained(art::Ptr<rb::Track> track); ///< Function to determine whether or not a track is contained
78 
79  };
80 }
81 
82 
83 namespace remid{
84 
85  //.......................................................................
86 
88  : EDAnalyzer(pset)
89  {
90  fSliceLabel = (pset.get< std::string >("SliceLabel"));
91  fTrackLabel = (pset.get< std::string >("TrackLabel"));
92  fPidLabel = (pset.get< std::string >("PidLabel" ));
93  }
94 
95  //.......................................................................
96 
98  {
99  }
100 
101  //.......................................................................
102 
104  {
106  fpid = tfs->make<TH1F>("fpid",";pid;NTracks",95,-0.1,1.1);
107  fpidnc = tfs->make<TH1F>("fpidnc",";pid;NTracks",95,-0.1,1.1);
108  fpidcc = tfs->make<TH1F>("fpidcc",";pid;NTracks",95,-0.1,1.1);
109  fpidcccut = tfs->make<TH1F>("fpidcccut",";pid;NTracks",95,-0.1,1.1);
110  fnscat = tfs->make<TH1F>("fnscat",";nscat;NTracks",301,-0.5,301);
111  fndedx = tfs->make<TH1F>("fndedx",";ndedx;NTracks",301,-0.5,301);
112  fcosscat = tfs->make<TH1F>("fcosscat",";cos #theta_{scat};",300,-1.1,1.1);
113  favgdedx = tfs->make<TH1F>("favgdedx",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
114  fdedx = tfs->make<TH1F>("fdedx",";dE/dx (GeV/cm);",500,0.0,0.01);
115  fpidnchigh = tfs->make<TH1F>("fpidnchigh",";pid;NTracks",95,-0.1,1.1);
116  fpidcchigh = tfs->make<TH1F>("fpidcchigh",";pid;NTracks",95,-0.1,1.1);
117  fpidhigh = tfs->make<TH1F>("fpidhigh",";pid;NTracks",95,-0.1,1.1);
118  fnscathigh = tfs->make<TH1F>("fnscathigh",";nscat;NTracks",301,-0.5,301);
119  fndedxhigh = tfs->make<TH1F>("fndedxhigh",";ndedx;NTracks",301,-0.5,301);
120  fcosscathigh = tfs->make<TH1F>("fcosscathigh",";cos #theta_{scat};",300,-1.1,1.1);
121  favgdedxhigh = tfs->make<TH1F>("favgdedxhigh",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
122  fdedxhigh = tfs->make<TH1F>("fdedxhigh",";dE/dx (GeV/cm);",500,0.0,0.01);
123  fpidmu = tfs->make<TH1F>("fpidmu",";pid;NTracks",95,-0.1,1.1);
124  fnscatmu = tfs->make<TH1F>("fnscatmu",";nscat;NTracks",301,-0.5,301);
125  fndedxmu = tfs->make<TH1F>("fndedxmu",";ndedx;NTracks",301,-0.5,301);
126  fcosscatmu = tfs->make<TH1F>("fcosscatmu",";cos #theta_{scat};",300,-1.1,1.1);
127  favgdedxmu = tfs->make<TH1F>("favgdedxmu",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
128  fdedxmu = tfs->make<TH1F>("fdedxmu",";dE/dx (GeV/cm);",500,0.0,0.01);
129  fpidmucont = tfs->make<TH1F>("fpidmucont",";pid;NTracks",95,-0.1,1.1);
130  fnscatmucont = tfs->make<TH1F>("fnscatmucont",";nscat;NTracks",301,-0.5,301);
131  fndedxmucont = tfs->make<TH1F>("fndedxmucont",";ndedx;NTracks",301,-0.5,301);
132  fcosscatmucont = tfs->make<TH1F>("fcosscatmucont",";cos #theta_{scat};",300,-1.1,1.1);
133  favgdedxmucont = tfs->make<TH1F>("favgdedxmucont",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
134  fdedxmucont = tfs->make<TH1F>("fdedxmucont",";dE/dx (GeV/cm);",500,0.0,0.01);
135  }
136 
137  //.......................................................................
138 
140  {
141  // get the slices out of the event
143  evt.getByLabel(fSliceLabel,slicecol);
144 
145  // get the tracks out of the event
147  evt.getByLabel(fTrackLabel,trkcol);
148 
149  art::PtrVector<rb::Cluster> slicelist;
150  art::PtrVector<rb::CellHit> Allevthits; // PtrVector of all of the hits in the event
151 
152  for(unsigned int i = 0; i<slicecol->size();++i){
153  art::Ptr<rb::Cluster>slice(slicecol,i);
154  slicelist.push_back(slice);
155  for(unsigned int ihit = 0; ihit<slice->NCell();++ihit){
156  Allevthits.push_back(slice->Cell(ihit));
157  }
158  }
159 
160  // get assosciations between slices and tracks
161  art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
163 
164  // loop over all slices and find the most effiecient slice for a neutrino interaction
165  double bestEff = -1;
166  unsigned int mostEffSlice = 0; // Index of the most efficient neutrino slice
167  int intType = 6; // This is a index to keep track of what type of interaction this is
168  // get the neutrino interaction type and slice it came from
169  for(unsigned int iSlice = 0;iSlice<slicelist.size();++iSlice){
170  // get the neutrino interaction from this slice
171  std::vector<cheat::NeutrinoEffPur> nuint = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),Allevthits);
172  if(nuint.size()>0 && nuint[0].efficiency>bestEff){
173  bestEff = nuint[0].efficiency;
174  mostEffSlice = iSlice;
175  // get the neutrino
176  const simb::MCNeutrino nu = nuint[0].neutrinoInt->GetNeutrino();
177  if(nu.CCNC() == simb::kNC){ intType = 5; }
178  else{
179  if(nu.Nu().PdgCode() == 12 || nu.Nu().PdgCode() == -12){ intType = 4; }
180  else if(nu.Mode() == 0){intType = 0;}
181  else if(nu.Mode() == 1){intType = 1;}
182  else if(nu.Mode() == 2){intType = 2;}
183  else if(nu.Mode() == 3){intType = 3;}
184  }
185  }
186  }
187 
188  // loop over the slices
189  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
190  if(intType <0){ continue; }
191  if(mostEffSlice != iSlice){ continue; }
192  // get all of the tracks associated to this slice
193  std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
194  if(tracks.size() == 0){ continue; }
195  art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);
196  double bestpid = -1.0;
197  int bestInd = -1;
198  // find the true primary muon track id if a CC interaction
199  std::set<int> muID;
200  if(intType < 4){
201  std::set<int> trackIDs = bt->GetTrackIDList();
202  for(std::set<int>::iterator it= trackIDs.begin(); it != trackIDs.end(); ++it){
203  // get the particle
204  const sim::Particle* part = bt->TrackIDToParticle(*it);
205  // check if this is the muon
206  // std::cout<<"On particle id: "<<*it<<" with pdg: "<<part->PdgCode()<<" and mother: "<<part->Mother()<<std::endl;
207  if((part->PdgCode() == 13 || part->PdgCode() == -13) && part->Mother() == 0){
208  muID.insert(*it);
209  break;
210  }
211  }
212  }
213  double bestmueff = -1.0;
214  int muInd = -1;
215  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
216  art::Ptr<remid::ReMId> pidobj = foids.at(itrack);
217  fpid->Fill(pidobj->Value());
218  if(intType < 4){
219  fpidcc->Fill(pidobj->Value());
220  // find if this track is the most efficiently reconstructed muon
221  double eff = bt->HitCollectionEfficiency(muID,tracks[itrack]->AllCells(),slicelist[iSlice]->AllCells(),geo::kXorY);
222  if(eff > bestmueff){
223  bestmueff = eff;
224  muInd = itrack;
225  }
226  }
227  else{ fpidnc->Fill(pidobj->Value());}
228  fnscat->Fill(pidobj->NScatters());
229  for(unsigned int i = 0; i<pidobj->NScatters(); ++i){
230  fcosscat->Fill(pidobj->CosScat(i));
231  }
232  favgdedx->Fill(pidobj->AvgDedx());
233  fndedx->Fill(pidobj->NDedx());
234  for(unsigned int i = 0; i<pidobj->NDedx(); ++i){
235  fdedx->Fill(pidobj->Dedx(i));
236  }
237  if(pidobj->Value() > bestpid){
238  bestInd = itrack;
239  bestpid = pidobj->Value();
240  }
241 
242  }
243 
244  // make plots for highest pid track
245  if(bestInd <0){ continue; }
246  art::Ptr<remid::ReMId> pidobj = foids.at(bestInd);
247  fpidhigh->Fill(pidobj->Value());
248  if(intType < 4){fpidcchigh->Fill(pidobj->Value()); }
249  else{ fpidnchigh->Fill(pidobj->Value()); }
250  fnscathigh->Fill(pidobj->NScatters());
251  for(unsigned int i = 0; i<pidobj->NScatters(); ++i){
252  fcosscathigh->Fill(pidobj->CosScat(i));
253  }
254  favgdedxhigh->Fill(pidobj->AvgDedx());
255  fndedxhigh->Fill(pidobj->NDedx());
256  for(unsigned int i = 0; i<pidobj->NDedx(); ++i){
257  fdedxhigh->Fill(pidobj->Dedx(i));
258  }
259 
260 
261  // make plots for true mu track by efficiency
262  if(intType >3){ continue; }
263  if(muInd < 0){continue; }
264  art::Ptr<remid::ReMId> mupidobj = foids.at(muInd);
265  fpidmu->Fill(mupidobj->Value());
266  fnscatmu->Fill(mupidobj->NScatters());
267  for(unsigned int i = 0; i<mupidobj->NScatters(); ++i){
268  fcosscatmu->Fill(mupidobj->CosScat(i));
269  }
270  favgdedxmu->Fill(mupidobj->AvgDedx());
271  fndedxmu->Fill(mupidobj->NDedx());
272  for(unsigned int i = 0; i<mupidobj->NDedx(); ++i){
273  fdedxmu->Fill(mupidobj->Dedx(i));
274  }
275 
276  if(IsContained(tracks[muInd])){
277  fpidmucont->Fill(mupidobj->Value());
278  fnscatmucont->Fill(mupidobj->NScatters());
279  for(unsigned int i = 0; i<mupidobj->NScatters(); ++i){
280  fcosscatmucont->Fill(mupidobj->CosScat(i));
281  }
282  favgdedxmucont->Fill(mupidobj->AvgDedx());
283  fndedxmucont->Fill(mupidobj->NDedx());
284  for(unsigned int i = 0; i<mupidobj->NDedx(); ++i){
285  fdedxmucont->Fill(mupidobj->Dedx(i));
286  }
287  }
288 
289  // make plots for contained qe sample
290  if(!IsContained(tracks[bestInd]) && !IsContained(tracks[0])){ continue; }
291  if(tracks[bestInd]->View() != geo::kXorY){ continue; }
292  if(intType != 0){continue; }
293  fpidcccut->Fill(pidobj->Value());
294 
295  } // end loop over slices
296 
297 
298  }
299 
300  //.......................................................................
301 
303  {
304  }
305 
306  //.......................................................................
307 
309  {
310  return a->NCell()>b->NCell();
311  }
312 
313  //.......................................................................
314 
316  {
317  double xmin,xmax,ymin,ymax,zmin,zmax;
318  // // define the boundaries of containment
319  // int concells = 3;
320  // art::ServiceHandle<geo::Geometry> geom;
321  // double detHalfWidth = geom->DetHalfWidth();
322  // double detHalfHeight = geom->DetHalfHeight();
323  // double detLength = geom->DetLength();
324  // // get the boundaries of the containment boundary
325  // const geo::PlaneGeo* pi = geom->Plane(concells-1);
326  // const geo::PlaneGeo* pf = geom->Plane(geom->NPlanes()-concells);
327  // // get zmin and zmax first
328  // const geo::CellGeo* ci = pi->Cell(concells-1);
329  // const geo::CellGeo* cl = pi->Cell(pi->Ncells()-1-concells);
330  // const geo::CellGeo* cf = pf->Cell(0);
331  // // double xyzf[3];
332  // double xyzi[3], xyzf[3], xyzl[3];
333  // ci->GetCenter(xyzi);
334  // cl->GetCenter(xyzl);
335  // cf->GetCenter(xyzf);
336  // std::cout<<"3 cell z center position: "<<xyzf[2]<<std::endl;
337  // double dzi = ci->HalfD();
338  // zmin = xyzi[2]+dzi;
339 
340  // double dzf = cf->HalfD();
341  // zmax = xyzf[2]-dzf;
342 
343  // if(pi->View() == geo::kX){
344  // xmin = xyzi[0]+ci->HalfW();
345  // xmax = xyzl[0]-cl->HalfW();
346  // int nextplane = geom->NextPlaneOtherView(concells-1);
347  // pi = geom->Plane(nextplane);
348  // ci = pi->Cell(concells-1);
349  // cl = pi->Cell(pi->Ncells()-1-concells);
350  // ci->GetCenter(xyzi);
351  // cl->GetCenter(xyzl);
352  // ymin = xyzi[1]+ci->HalfW();
353  // ymax = xyzl[1]-ci->HalfW();
354  // }
355  // else if(pi->View() == geo::kY){
356  // ymin = xyzi[1]+ci->HalfW();
357  // ymax = xyzl[1]-cl->HalfW();
358  // int nextplane = geom->NextPlaneOtherView(concells-1);
359  // pi = geom->Plane(nextplane);
360  // ci = pi->Cell(concells-1);
361  // cl = pi->Cell(pi->Ncells()-1-concells);
362  // ci->GetCenter(xyzi);
363  // cl->GetCenter(xyzl);
364  // xmin = xyzi[0]+ci->HalfW();
365  // xmax = xyzl[0]-ci->HalfW();
366  // }
367 
368  // contianement cuts for a 28 block far detector. Track ends must be 3 cells + 100 cm from detector edge
369  xmin = -750.213;
370  xmax = 750.213;
371  ymin = -750.213;
372  ymax = 750.213;
373  zmin = 19.9156;
374  zmax = 6201.02; // 30 blocks
375  zmax = 5945.29; // 28 blocks
376  // zmax = detLength-zmin;
377  TVector3 start = track->Start();
378  TVector3 stop = track->Stop();
379  if(start.X() > xmin+100.00 && start.X() < xmax-100.00 && start.Y() > ymin+100.00 && start.Y() < ymax-100.00 && start.Z() > zmin+100.00 && start.Z() < zmax-100.00 &&
380  stop.X() > xmin && stop.X() < xmax && stop.Y() > ymin && stop.Y() < ymax && stop.Z() > zmin && stop.Z() < zmax){ return true; }
381  else{ return false; }
382 
383  }
384 
385 
386 } // end of remid namespace
387 
388 namespace remid
389 {
391 }
back track the reconstruction to the simulation
double Dedx(unsigned int i) const
Definition: ReMId.cxx:126
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
unsigned int NScatters() const
Definition: ReMId.cxx:96
std::map< std::string, double > xmax
set< int >::iterator it
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
std::string fTrackLabel
Where to find recondtructed tracks.
int Mother() const
Definition: MCParticle.h:212
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::string fPidLabel
Where to find the pids.
Definition: event.h:19
static bool SortByHits(const art::Ptr< rb::Track > &a, const art::Ptr< rb::Track > &b)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
double AvgDedx() const
Definition: ReMId.cxx:68
double Value() const
Definition: PID.h:22
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...
Double_t ymax
Definition: plot.C:25
TString part[npart]
Definition: Style.C:32
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
A module to analyze remid objects produced by the RecoMuon module.
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
double CosScat(unsigned int i) const
Definition: ReMId.cxx:114
ReMIdValidate(fhicl::ParameterSet const &pset)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Definition: View.py:1
size_type size() const
Definition: PtrVector.h:308
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
unsigned int NDedx() const
Definition: ReMId.cxx:102
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
TH1F * fpid
Pid value of every track.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void analyze(art::Event const &evt)
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fSliceLabel
Where to find reconstructed slices.
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
Double_t ymin
Definition: plot.C:24
A PID for muons.
Definition: FillPIDs.h:10
bool IsContained(art::Ptr< rb::Track > track)
Function to determine whether or not a track is contained.
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149