ReMIdTrain_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////
2 // \file ReMIdTrain_module.cc
3 // \brief A module for creating trees which can be used for kNN training and generating sample histograms for ReMId
4 // \version $Id: ReMIDTrain_module.cc,v 2013/02/21 16:30:29 raddatz Exp $
5 // \author Nicholas Raddatz - raddatz@physics.umn.edu
6 ////////////////////////////////////////////////////////////////////
7 #include <string>
8 
9 // ROOT includes
10 #include "TVector3.h"
11 #include "TTree.h"
12 #include "TMath.h"
13 
14 // Framework includes
15 #include "fhiclcpp/ParameterSet.h"
17 #include "Utilities/AssociationUtil.h"
21 
22 // NOvA includes
23 #include "Geometry/Geometry.h"
24 #include "Geometry/LiveGeometry.h"
26 #include "ReMId/ReMId.h"
27 #include "GeometryObjects/Geo.h"
30 #include "MCCheater/BackTracker.h"
31 #include "RecoBase/Track.h"
32 #include "Simulation/Simulation.h"
33 
34 
35 namespace remid {
36 
37  // A module to analyze ReMId objects for pid training
38  class ReMIdTrain : public art::EDAnalyzer {
39  public:
40 
41  explicit ReMIdTrain(fhicl::ParameterSet const& pset);
42  virtual ~ReMIdTrain();
43 
44  void beginJob();
45  void analyze(art::Event const& evt);
46  void endJob();
47 
48  private:
49  std::string fSliceLabel; ///<Where to find reconstructed slices
50  std::string fTrackLabel; ///<Where to find recondtructed tracks
51  std::string fPidLabel; ///<Where to find the pids
52  bool fCreateTrainTrees; ///<Create kNN training tree?
53  bool fCreateSampleTrees; ///<Create trees for sample histograms?
54 
55  TTree* TreeS; ///< Tree for holding signal input variables for kNN training
56  TTree* TreeB; ///< Tree for holding background input variables for kNN training
57 
58  // variables used in filling the kNN training trees
59  float fTrackLength;
60  float fDedxSep;
61  float fScatSep;
62  float fnuE;
63  float fvisE;
64  float frecoE;
65  float fMeasFrac;
66 
67  TTree* DedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
68  TTree* ScatSample; ///< Tree for holding scattering variables needed for making signal sample histograms
69 
70  // variables used in filling the sample histogram trees
71  float fDedx;
72  int fnDedx;
74  float fDistDedx;
75  float fCosScat;
77  float fScatVar;
78  int fnScat;
79  float fDistScat;
80  int fnPlane;
81  int fpdg;
82  int fintType;
83  bool fVert;
84  bool fUsed;
85 
89 
90  unsigned int fdibfirst = 0;
91  unsigned int fdiblast = 0;
92 
94 
95  static bool SortByHits(const art::Ptr<rb::Track>& a, const art::Ptr<rb::Track>& b);
96 
97 
98  };
99 
100 }
101 
102 
103 namespace remid{
104 
105  //.......................................................................
106 
107  struct Counts
108  {
110  {
111  nHits[0] = nHits[1] = 0;
112  nOnlyHits[0] = nOnlyHits[1] = 0;
113  }
114  int nHits[2]; // Indexed by view
115  int nOnlyHits[2]; // Hit by only one particle
116  };
117 
118  //.......................................................................
119 
121  EDAnalyzer(pset)
122  {
123  fSliceLabel = (pset.get< std::string >("SliceLabel"));
124  fTrackLabel = (pset.get< std::string >("TrackLabel"));
125  fPidLabel = (pset.get< std::string >("PidLabel" ));
126  fCreateTrainTrees = (pset.get< bool >("CreateTrainTrees"));
127  fCreateSampleTrees = (pset.get< bool >("CreateSampleTrees"));
128  }
129 
130  //.......................................................................
131 
133  {
134  }
135 
136  //.......................................................................
137 
139  {
141 
142  if(fCreateTrainTrees){
143  TreeS = tfs->make<TTree>("fTreeS","kNN Input Variables for Signal");
144  TreeS->Branch("trackLength", &fTrackLength,"trackLength/F");
145  TreeS->Branch("dedxSep", &fDedxSep, "dedxSep/F");
146  TreeS->Branch("scatSep", &fScatSep, "scatSep/F");
147  TreeS->Branch("nuE", &fnuE, "nuE/F");
148  TreeS->Branch("intType", &fintType,"intType/I");
149  TreeS->Branch("visE", &fvisE, "visE/F");
150  TreeS->Branch("recoE", &frecoE, "recoE/F");
151  TreeS->Branch("measFrac", &fMeasFrac, "measFrac/F");
152 
153  TreeB = tfs->make<TTree>("fTreeB","kNN Input Variables for Background");
154  TreeB->Branch("trackLength", &fTrackLength,"trackLength/F");
155  TreeB->Branch("dedxSep", &fDedxSep, "dedxSep/F");
156  TreeB->Branch("scatSep", &fScatSep, "scatSep/F");
157  TreeB->Branch("nuE", &fnuE, "nuE/F");
158  TreeB->Branch("intType", &fintType, "intType/I");
159  TreeB->Branch("pdg", &fpdg, "pdg/I");
160  TreeB->Branch("visE", &fvisE, "visE/F");
161  TreeB->Branch("recoE", &frecoE, "recoE/F");
162  TreeB->Branch("measFrac", &fMeasFrac, "measFrac/F");
163 
164 
165  }
166 
167  if(fCreateSampleTrees){
168  DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
169  DedxSample->Branch("dedx", &fDedx, "dedx/F");
170  DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
171  DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
172  DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
173  DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
174  DedxSample->Branch("nPlane", &fnPlane, "nPlane/I");
175  DedxSample->Branch("pdg", &fpdg, "pdg/I");
176  DedxSample->Branch("intType", &fintType, "intType/I");
177  DedxSample->Branch("vert", &fVert, "vert/O");
178  DedxSample->Branch("used", &fUsed, "used/O");
179  DedxSample->Branch("nuE", &fnuE, "nuE/F");
180 
181  ScatSample = tfs->make<TTree>("fScatSample","Variables for Creating Scattering Sample Histograms");
182  ScatSample->Branch("cosScat", &fCosScat, "cosSact/F");
183  ScatSample->Branch("distFromLastScat", &fDistLastScat, "distFromLastScat/F");
184  ScatSample->Branch("scatVar", &fScatVar, "scatVar/F");
185  ScatSample->Branch("nScat", &fnScat, "nScat/I");
186  ScatSample->Branch("distFromEnd", &fDistScat, "distFromEnd/F");
187  ScatSample->Branch("trackLength", &fTrackLength, "trackLength/F");
188  ScatSample->Branch("nPlane", &fnPlane, "nPlane/I");
189  ScatSample->Branch("pdg", &fpdg, "pdg/I");
190  ScatSample->Branch("intType", &fintType, "intType/I");
191 
192  }
193 
194  }
195 
196  //.......................................................................
197 
199  {
200  // get the slices out of the event
202  evt.getByLabel(fSliceLabel,slicecol);
203 
204  // // get the tracks out of the event
205  // art::Handle<std::vector<rb::Track> > trkcol;
206  // evt.getByLabel(fTrackLabel,trkcol);
207 
208  art::PtrVector<rb::Cluster> slicelist;
209  art::PtrVector<rb::CellHit> Allevthits; // PtrVector of all of the hits in the event
210 
211  for(unsigned int i = 0; i<slicecol->size();++i){
212  art::Ptr<rb::Cluster>slice(slicecol,i);
213  slicelist.push_back(slice);
214  for(unsigned int ihit = 0; ihit<slice->NCell();++ihit){
215  Allevthits.push_back(slice->Cell(ihit));
216  }
217  }
218 
219  // get assosciations between slices and tracks
220  art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
221 
222  // get the detector configuration
223  unsigned int dibmask = (frh->GoodDiBlockMask()&frh->GetConfiguration());
224  if(dibmask == ((1<<16)-1)){ dibmask = frh->GetConfiguration(); }
225  if(dibmask == 0){ dibmask = frh->GetConfiguration(); }
226  if(dibmask){
227  int iD;
228  iD=0;
229  while(!((dibmask>>iD)&1)){
230  iD++;
231  fdibfirst=iD+1;
232  }
233  iD=0;
234  while(dibmask>>iD){
235  iD++;
236  fdiblast=iD;
237  }
238  }
239 
240 
242 
243  // loop over all slices and find the most effiecient slice for a neutrino interaction
244  double bestEff = 0;
245  unsigned int mostEffSlice = 0; // Index of the most efficient neutrino slice
246  int nuintType = 6; // This is a index to keep track of what type of interaction this is
247  // get the neutrino interaction type
248  for(unsigned int iSlice = 0;iSlice<slicelist.size();++iSlice){
249  if(slicelist[iSlice]->IsNoise()){ continue; }
250  // get the neutrino interaction from this slice
251  std::vector<cheat::NeutrinoEffPur> nuint = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),Allevthits);
252  if(nuint.size() > 0 && nuint[0].efficiency > bestEff){
253  bestEff = nuint[0].efficiency;
254  mostEffSlice = iSlice;
255  // get the neutrino
256  const simb::MCNeutrino nu = nuint[0].neutrinoInt->GetNeutrino();
257  fnuE = nu.Nu().E();
258  fvisE = nuint[0].energySlice;
259  if(nu.CCNC() == simb::kNC){ nuintType = 5; }
260  else{
261  if(nu.Nu().PdgCode() == 12 || nu.Nu().PdgCode() == -12){nuintType = 4;}
262  else if(nu.Mode() == 0){nuintType = 0;}
263  else if(nu.Mode() == 1){nuintType = 1;}
264  else if(nu.Mode() == 2){nuintType = 2;}
265  else if(nu.Mode() == 3){nuintType = 3;}
266  }
267  }
268  }
269 
270  // loop over the slices
271  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
272  // if we are not on the most efficient slice don't try to do anything
273  if(iSlice != mostEffSlice){ continue; }
274  // also if the efficiency is not greater than 1 or is a noise slice skip it
275  if(bestEff == 0){ continue; }
276  // If slice is noise, skip it
277  if(slicelist[iSlice]->IsNoise()){ continue; }
278  // If don't know what neutrino this slice comes from skip it
279  if(nuintType == 6){ continue; }
280 
281  // get all of the tracks associated to this slice
282  std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
283  // if no tracks move on
284  if(tracks.size() == 0){ continue; }
285 
286  frecoE = slicelist[iSlice]->TotalGeV();
287 
288  // sort the tracks by the number of hits
289  std::sort(tracks.begin(),tracks.end(),SortByHits);
290 
291  //get the trackIDs for the event
292  std::set<int> trackIDs = bt->GetTrackIDList();
293 
294  // don't penalize purity because of michel electron
295  std::set<int> allowedDaughters;
296  allowedDaughters.insert(-11);
297  allowedDaughters.insert(11);
298 
299  // // get the id of the track to use as either signal or background
300  // int trackid = -1;
301  // // get the reco track id of the particle
302  // int recotrackid = -1;
303 
304  std::set<int> primTrkIDs;
305  // Get the set of primary particle track ids for this event
306  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end();++it){
307  if(*it == bt->TrackIDToMotherParticle(*it)->TrackId()){ primTrkIDs.insert(*it); }
308  }
309 
310  // match tracks to these primary particles
311  std::vector<cheat::ParticleEffPur> pep = bt->TracksToParticles(tracks, primTrkIDs, allowedDaughters, slicelist[iSlice]->AllCells());
312 
313  // get the pid that belongs to this recotrack
314  art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);
315 
316  // loop over all of the tracks
317  for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
318 
319  fpdg = pep[itrk].pdg;
320  fintType = nuintType;
321 
322  // check if this track is matched to a primary particle
323  if(fpdg == 0){ continue; }
324 
325  // Be safe, check that efficiency and purity are both greater than 0
326  if(!(pep[itrk].efficiency > 0 || pep[itrk].purity > 0)){ continue; }
327 
328  // get the pid that belongs to this track
329  art::Ptr<remid::ReMId> pid = foids.at(itrk);
330 
331  fnPlane = tracks[itrk]->ExtentPlane();
332 
333  // fill trees with variable for creating sample histograms if wanted
334  if(fCreateSampleTrees){
335 
336  // Create sample histograms from contained tracks only
337 
338  // check if this track is contained by pid
339  if(!(pid->Contained())){ continue; }
340 
341  // any other track cuts?
342  fTrackLength = tracks[itrk]->TotalLength();
343 
344  // cut out strange tracklengths
345  if(tracks[itrk]->TotalLength() > 6500){ continue; }
346 
347  // make sure its a 3d track
348  if(!tracks[itrk]->Is3D()){ continue; }
349 
350  // get all the dedx variables filled
351  fnDedx = pid->NDedx()-1;
353  // loop over all of the dE/dx measurements
354  for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
355  fDedx = pid->Dedx(idedx);
356  fVert = pid->DedxVertex(idedx);
357  fUsed = pid->DedxUsed(idedx);
358  // get the distance from the end of the track
359  fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[itrk]->TotalLength();
360  // fill the dedx sample tree
361  DedxSample->Fill();
362  }
363 
364  // get all the scattering variables filled
365  fnScat = pid->NScatters();
366  for(unsigned int iscat = 0; iscat<pid->NScatters();++iscat){
367  fCosScat = pid->CosScat(iscat);
368  double ang = TMath::ACos(fCosScat);
369  // get the distance from the end of the track and from the last scattering measurement
370  TVector3 trkstart = tracks[itrk]->Start();
371  double dist = 0;
372  double lastdist = 0;
373  for(unsigned int itraj = 0; itraj<tracks[itrk]->NTrajectoryPoints()-1; ++itraj){
374  TVector3 p0 = tracks[itrk]->TrajectoryPoint(itraj);
375  TVector3 p1 = tracks[itrk]->TrajectoryPoint(itraj+1);
376  lastdist = (p1-p0).Mag();
377  dist+=lastdist;
378  if(p1 == pid->ScatLocation(iscat)){break;}
379  }
380  fDistLastScat = lastdist;
381  fDistScat = tracks[itrk]->TotalLength()-dist;
382  fScatVar = ang*ang/lastdist;
383  // fill the scat sample tree
384  ScatSample->Fill();
385  }
386 
387  } // Fill sample histogram trees
388 
389 
390  } // end of loop over reco tracks
391 
392 
393  if(fCreateTrainTrees){
394 
395  fintType = nuintType;
396  if(fintType < 4){
397  // find the best muon track
398  int bestMu = -1;
399  for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
400  if(abs(pep[itrk].pdg) == 13){ bestMu = itrk; }
401  }
402  if(bestMu >= 0){
403  // get the pid that belongs to this track
404  art::Ptr<remid::ReMId> pid = foids.at(bestMu);
405  // check if this pid has any usable measurements of dedx and that the track is 3d
406  if(pid->NMeasurementPlanes() > 0 && tracks[bestMu]->Is3D()){
407  fMeasFrac = ((float)pid->NMeasurementPlanes())/((float)pid->NMeasurementPlanes()+(float)pid->NVertexPlanes());
408  // check if this track is contained by pid
409  // check if this event is contained
410  if(ContainedEvent(slicelist[iSlice],tracks[bestMu])){
411  // any other track cuts?
412  fTrackLength = tracks[bestMu]->TotalLength();
413  // cut out strange tracklengths
414  if(fTrackLength < 6500){
415  fDedxSep = pid->DedxSeparation();
416  fScatSep = pid->ScatSeparation();
417  TreeS->Fill();
418  }
419  } // if contained
420  }
421  } // if muon matched track
422  } // CC event
423  else{
424  // non CC event
425  int bestPi = -1;
426  double longestPi = -1;
427  int longestTrkIdx = -1;
428  double longestTrk = -1;
429  // find the most muon like track
430  // first look for a primary pion match
431  for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
432  double trklen = tracks[itrk]->TotalLength();
433  if(abs(pep[itrk].pdg) == 211 && trklen > longestPi){
434  bestPi = itrk;
435  longestPi = trklen;
436  }
437  if(trklen > longestTrk){
438  longestTrkIdx = itrk;
439  longestTrk = trklen;
440  }
441  } // end of loop over tracks
442 
443  // train against longest pion track if it exist
444  if(bestPi >= 0){
445  // get the pid that belongs to this track
446  art::Ptr<remid::ReMId> pid = foids.at(bestPi);
447  // check if this pid has any usable measurements of dedx and that the track is 3d
448  if(pid->NMeasurementPlanes() > 0 && tracks[bestPi]->Is3D()){
449  fMeasFrac = ((float)pid->NMeasurementPlanes())/((float)pid->NMeasurementPlanes()+(float)pid->NVertexPlanes());
450  // check if this track is contained by pid
451  if(ContainedEvent(slicelist[iSlice],tracks[bestPi])){
452  fTrackLength = longestPi;
453  // cut out strange tracklengths
454  if(fTrackLength < 6500){
455  fpdg = pep[bestPi].pdg;
456  fDedxSep = pid->DedxSeparation();
457  fScatSep = pid->ScatSeparation();
458  TreeB->Fill();
459  }
460  } // if contained
461  }
462  } // if pion exists
463  else if(longestTrkIdx >= 0){
464  // get the pid that belongs to this track
465  art::Ptr<remid::ReMId> pid = foids.at(longestTrkIdx);
466  // check if this pid has any usable measurements of dedx and that the track is 3d
467  if(pid->NMeasurementPlanes() > 0 && tracks[longestTrkIdx]->Is3D()){
468  fMeasFrac = ((float)pid->NMeasurementPlanes())/((float)pid->NMeasurementPlanes()+(float)pid->NVertexPlanes());
469  // check if this track is contained by pid
470  if(ContainedEvent(slicelist[iSlice],tracks[longestTrkIdx])){
471  fTrackLength = longestTrk;
472  // cut out strange tracklengths
473  if(fTrackLength < 6500){
474  fpdg = pep[longestTrkIdx].pdg;
475  fDedxSep = pid->DedxSeparation();
476  fScatSep = pid->ScatSeparation();
477  TreeB->Fill();
478  }
479  } // if contained
480  }
481  } // if long track exists
482 
483  } // NC event
484 
485  } // if fill trees for kNN training
486 
487 
488  } // end loop over slices
489 
490 
491  }
492 
493  //.......................................................................
494 
496  {
497  }
498 
499  //.......................................................................
500 
502  {
503  return a->NCell()>b->NCell();
504  }
505 
506  //.......................................................................
508  {
509 
510  bool contained = false;
511 
512  // check that this is a contained event
513 
514  // get some useful information here
515  // 1) how close is the slice to the edge of the detector
516  unsigned int slcncellsfromedge = 999999999;
517  for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
518  const art::Ptr<rb::CellHit>& chit = slice->Cell(ihit);
519  const int plane = chit->Plane();
520  unsigned int cellsfromedge = std::min((unsigned int)chit->Cell(), fgeom->Plane(plane)->Ncells() - 1 - chit->Cell());
521  if(cellsfromedge < slcncellsfromedge){
522  slcncellsfromedge = cellsfromedge;
523  }
524  }
525  // 2) how close is the track to the edge of the detector
526  // in planes
527  unsigned int firstpl = trk->MinPlane();
528  unsigned int lastpl = trk->MaxPlane();
529  // in projected cells, just define directions here actually get projected cell below
530  TVector3 dirbwk = -trk->Dir();
531  TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();
532 
534  // cut the slice number of cells from edge
535  if(slcncellsfromedge > 1){
536  unsigned int npltofront = firstpl - (fdibfirst-1)*64;
537  unsigned int npltoback = (fdiblast*64) - 1 - lastpl;
538  // check that the start of track is not the front or back of the detector
539  if(npltofront > 1 && npltoback > 1){
540  int fwdcell = flivegeom->ProjectedLiveCellsToEdge(trk->Stop(),dirfwd);
541  int bwkcell = flivegeom->ProjectedLiveCellsToEdge(trk->Start(),dirbwk);
542  if(fwdcell > 10 && bwkcell > 10){
543  contained = true;
544  } // end if track ends projections contained
545  } // end if track end planes contained
546  } // end if slice contained
547  } // end if far detector
548  else if(fgeom->DetId() == novadaq::cnv::kNEARDET){
549  // cut the slice number of cells from edge
550  std::cout<<" slice cells from edge: "<<slcncellsfromedge<<std::endl;
551  if(slcncellsfromedge > 1){
552  // check that first diblock and last diblock give you the full detector
553  std::cout<<"fdibfirst: "<<fdibfirst<<" fdiblast: "<<fdiblast<<std::endl;
554  // if(fdibfirst == 1 && fdiblast == 4){
555  // check that the start of track is not the front or back of the detector
556  std::cout<<"firstpl: "<<firstpl<<" lastpl: "<<lastpl<<std::endl;
557  if(firstpl > 1 && lastpl < 212){
558  int fwdcell = flivegeom->FullNDProjectedCells(trk->Stop(), dirfwd);
559  int bwkcell = flivegeom->FullNDProjectedCells(trk->Start(),dirbwk);
560  // also make sure the track doesn't start too close to the end of the active region
561  std::cout<<"fwdcell: "<<fwdcell<<" bwkcell: "<<bwkcell<<" start Z: "<<trk->Start().Z()<<std::endl;
562  if(fwdcell > 4 && bwkcell > 8 && trk->Start().Z() < 1150){
563  std::cout<<"stop Z: "<<trk->Stop().Z()<<std::endl;
564  // if the track ends in the muon catcher, check that it doesn't clip the corner
565  if(trk->Stop().Z() < 1275){
566  std::cout<<"contained"<<std::endl;
567  contained = true;
568  }
569  else{
570 
571  // find the trajectory point closest in z to the transition plane
572  double mindist = 99999999;
573  unsigned int mintrajid = 0;
574  for(unsigned int itraj = 0; itraj < trk->NTrajectoryPoints()-1; ++itraj){
575  double dist = 1275 - trk->TrajectoryPoint(itraj).Z();
576  if(dist > 0 && dist < mindist){
577  mindist = dist;
578  mintrajid = itraj;
579  }
580  }
581  // get the direction at this trajectory point
582  TVector3 trandir = trk->Dir();
583  if(mintrajid < trk->NTrajectoryPoints() -2){
584  TVector3 trandir = (trk->TrajectoryPoint(mintrajid+1) - trk->TrajectoryPoint(mintrajid)).Unit();
585  }
586  double ypostran = flivegeom->YPositionAtTransition(trk->TrajectoryPoint(mintrajid),trandir);
587  std::cout<<"ypostran: "<<ypostran<<std::endl;
588  if( ypostran < 55){
589  std::cout<<"contained"<<std::endl;
590  contained = true;
591  }
592 
593  } // end else in muon catcher
594  } // end if track ends projections contained
595  } // end if track end planes contained
596  // } // end if expected number of diblocks
597  } // end if slice contained
598  } // end if near detector
599 
600  return contained;
601 
602  }
603 
604 
605 }
606 
607 
608 namespace remid{
609 
611 
612 }
double E(const int i=0) const
Definition: MCParticle.h:232
size_t NTrajectoryPoints() const
Definition: Track.h:83
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
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
set< int >::iterator it
bool fCreateTrainTrees
Create kNN training tree?
bool fCreateSampleTrees
Create trees for sample histograms?
const sim::Particle * TrackIDToMotherParticle(int const &id) const
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...
unsigned short Plane() const
Definition: CellHit.h:39
unsigned int fdibfirst
double ScatSeparation() const
Return the scattering separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:206
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
art::ServiceHandle< geo::Geometry > fgeom
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
std::string fTrackLabel
Where to find recondtructed tracks.
void abs(TH1 *hist)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
int TrackId() const
Definition: MCParticle.h:209
double dist
Definition: runWimpSim.h:113
bool DedxVertex(unsigned int i) const
Definition: ReMId.cxx:138
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
Far Detector at Ash River, MN.
TTree * TreeS
Tree for holding signal input variables for kNN training.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
int NMeasurementPlanes() const
Definition: ReMId.cxx:223
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 ScatLocation(unsigned int i) const
Definition: ReMId.cxx:108
int evt
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
Near Detector in the NuMI cavern.
double CosScat(unsigned int i) const
Definition: ReMId.cxx:114
Collect Geo headers and supply basic geometry functions.
TTree * ScatSample
Tree for holding scattering variables needed for making signal sample histograms. ...
int NVertexPlanes() const
Definition: ReMId.cxx:233
bool DedxUsed(unsigned int i) const
Definition: ReMId.cxx:132
bool Contained() const
Definition: ReMId.cxx:243
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TVector3 Unit() const
size_type size() const
Definition: PtrVector.h:308
ReMIdTrain(fhicl::ParameterSet const &pset)
OStream cout
Definition: OStream.cxx:6
TTree * TreeB
Tree for holding background input variables for kNN training.
unsigned int NDedx() const
Definition: ReMId.cxx:102
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
art::ServiceHandle< geo::LiveGeometry > flivegeom
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string fSliceLabel
Where to find reconstructed slices.
T * make(ARGS...args) const
double DedxLocation(unsigned int i) const
Definition: ReMId.cxx:120
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
std::string fPidLabel
Where to find the pids.
float Mag() const
void analyze(art::Event const &evt)
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
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
art::ServiceHandle< nova::dbi::RunHistoryService > frh
A PID for muons.
Definition: FillPIDs.h:10
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
int GoodDiBlockMask(int subrun=-1, bool reload=false)
Event generator information.
Definition: MCNeutrino.h:18
void efficiency()
Definition: efficiency.C:58
int Mode() const
Definition: MCNeutrino.h:149
std::vector< ParticleEffPur > TracksToParticles(const std::vector< const rb::Track * > &tracks, const std::set< int > &trkIDs, const std::set< int > &allowedDaughters, const std::vector< const rb::CellHit * > &allHits) const
Given a collection of reconstructed tracks, this method finds the best match for each reco track to t...
double DedxSeparation() const
Return the dE/dx separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:199
Encapsulate the geometry of one entire detector (near, far, ndos)
static bool SortByHits(const art::Ptr< rb::Track > &a, const art::Ptr< rb::Track > &b)