ReMIdDedxRock_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////
2 // \file ReMIdDedxRock_module.cc
3 // \brief A module for creating trees which can be used to study dE/dx for Rock Muons
4 // Adapted from ReMIdDedx_module.cc
5 // \version $Id: ReMIdDedxRock_module.cc,v 2017/06/20 sgermani Exp$
6 // \author Stefano Germani - s.germani@ucl.ac.uk
7 ////////////////////////////////////////////////////////////////////
8 #include <string>
9 
10 // ROOT includes
11 #include "TVector3.h"
12 #include "TTree.h"
13 #include "TH1D.h"
14 #include "TH2D.h"
15 #include "TMath.h"
16 
17 // Framework includes
18 #include "fhiclcpp/ParameterSet.h"
20 #include "Utilities/AssociationUtil.h"
25 
26 // NOvA includes
27 //....................................
28 #include "Calibrator/Calibrator.h"
30 //....................................
31 #include "Geometry/Geometry.h"
32 #include "Geometry/LiveGeometry.h"
34 #include "ReMId/ReMId.h"
35 #include "NumuEnergy/NumuE.h"
36 #include "QEEventFinder/QePId.h"
37 #include "GeometryObjects/Geo.h"
40 #include "MCCheater/BackTracker.h"
41 #include "RecoBase/Track.h"
42 //................................
43 #include "RecoBase/Cluster.h"
44 #include "RecoBase/CellHit.h"
45 #include "RecoBase/RecoHit.h"
46 
47 #include "MCCheater/SimHit.h"
48 
49 //.................................
50 
51 #include "Simulation/Simulation.h"
53 #include "SummaryData/POTSum.h"
54 #include "SummaryData/SpillData.h"
55 #include "MCCheater/BackTracker.h"
56 
57 //....
58 
59 namespace remid {
60 
61  // A module to analyze ReMId objects for pid training
62  class ReMIdDedxRock : public art::EDAnalyzer {
63  public:
64  explicit ReMIdDedxRock(fhicl::ParameterSet const& pset);
66 
67  void analyze(art::Event const& evt);
68 
69  void beginJob();
70  void beginRun(const art::Run& run);
71  void endJob();
72  //virtual void beginSubRun(art::SubRun& sr);
73 
74  private:
76  std::string fSliceLabel; ///<Where to find reconstructed slices
77  std::string fTrackLabel; ///<Where to find recondtructed tracks
78  std::string fPidLabel; ///<Where to find the pids
84 
85  TH1D* fPOT; ///<Histogram of POT passing spill cuts
92 
93  double fTotalPOT;
94  bool fND; ///< Is detector ND?
95 
96  // variables used in filling the kNN training trees
97  float fTrackLength;
98  //float fProtonTrackLength;
99  //float fDedxSep;
100  //float fScatSep;
101  //float fnuE;
102  //float fvisE;
103  //float frecoE;
104  //float fMeasFrac;
105 
106  TTree* DedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
107  //TTree* ProtonDedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
108 
109  // variables used in filling the sample histogram trees
110  float fDedx;
111  int fnDedx;
113  float fDistDedx;
114  //float fProtonDedx;
115  //int fProtonnDedx;
116  //int fProtonnDedxUsed;
117  //float fProtonDistDedx;
118  //float fCosScat;
119  //float fDistLastScat;
120  //float fScatVar;
121  //int fnScat;
122  //float fDistScat;
123  //int fnPlane;
124  int fpdg;
125  //int fintType;
126  bool fVert;
127  bool fUsed;
128  float fzStart;
129  float fzEnd;
130  //bool fProtonVert;
131  //bool fProtonUsed;
132  //float fProtonzStart;
133  //float fProtonzEnd;
146 
147  // variables added to trees for detailed dE/dx studies
149  float fmyDedx;
150  float fDe;
151  int fView;
152  int fPlane;
153  int fCell;
154  float fPE;
155  float fADC;
156  float fPECorr;
157  float fMIP;
158  float fDx;
159  float fxDead;
160  float fxStart;
161  float fxEnd;
162  float fyStart;
163  float fyEnd;
164  float fxHit;
165  float fyHit;
166  float fzHit;
167  float ftHit;
168  float fxDir;
169  float fyDir;
170  float fzDir;
171  float fmcDe;
172  float fmcDeBirks;
173  float fmcDx;
174 
175  float fTrueE;
177 
178  std::vector <int> fMultiCell;
179  std::vector <float> fMultiPE;
180  std::vector <float> fMultiADC;
181  std::vector <float> fMultiPECorr;
182 
183  // ......................
184 
189 
191 
192  //novadaq::cnv::DetId fDetID; ///< Detector ID in nova daq convention namespace typedef
193 
194 
195  //unsigned int fdibfirst = 0;
196  //unsigned int fdiblast = 0;
197 
199 
200  };
201 
202 }
203 
204 
205 namespace remid{
206 
207 
208  //.......................................................................
209 
211  EDAnalyzer(pset)
212  {
213  fHitModuleLabel = (pset.get< std::string >("HitModuleLabel"));
214  fSliceLabel = (pset.get< std::string >("SliceLabel"));
215  fTrackLabel = (pset.get< std::string >("TrackLabel"));
216  fPidLabel = (pset.get< std::string >("PidLabel" ));
217  fGeneratorLabel = pset.get< std::string >("GeneratorInput");
218  fNumiBeamLabel = pset.get< std::string >("NuMIBeamInput");
219  fSpillQualLabel = pset.get< std::string >("SpillQualInput");
220  fNumuEnergyLabel = pset.get< std::string >("NumuEnergyInput");
221  fQepidLabel = pset.get< std::string >("QepidInput");
222  }
223 
224  //.......................................................................
225 
227  {
228 
229  }
230 
231  //.......................................................................
232 
234  {
236 
237  DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
238  DedxSample->Branch("dedx", &fDedx, "dedx/F");
239  DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
240  DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
241  DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
242  DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
243  DedxSample->Branch("vert", &fVert, "vert/O");
244  DedxSample->Branch("used", &fUsed, "used/O");
245  DedxSample->Branch("zStart",&fzStart,"zStart/F");
246  DedxSample->Branch("zEnd",&fzEnd,"zEnd/F");
247  DedxSample->Branch("passOffTrackHits",&fPassOffTrackHits,"passOffTrackHits/O");
248  DedxSample->Branch("passOffTrackEnergy",&fPassOffTrackEnergy,"passOffTrackEnergy/O");
249  DedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");
250  DedxSample->Branch("passVertX",&fPassVertX,"passVertX/O");
251  DedxSample->Branch("passVertY",&fPassVertY,"passVertY/O");
252  DedxSample->Branch("passVertZ",&fPassVertZ,"passVertZ/O");
253  DedxSample->Branch("passProtonRemid",&fPassProtonRemid,"passProtonRemid/O");
254  DedxSample->Branch("passDedxRatio",&fPassDedxRatio,"passDedxRatio/O");
255  DedxSample->Branch("passEDiffZTest",&fPassEDiffZTest,"passEDiffZTestX/O");
256  DedxSample->Branch("trueNumuCCQE",&fTrueNumuCCQE,"trueNumuCCQE/O");
257 
258  // Branches added wrt standard ReMIdDedx module for detailed dE/dx studies
259  DedxSample->Branch("ncellsPerPlane", &fnCellsPerPlane, "ncellsPerPlane/I");
260  DedxSample->Branch("mydedx", &fmyDedx, "mydedx/F");
261  DedxSample->Branch("de", &fDe, "de/F");
262  DedxSample->Branch("view", &fView, "view/I");
263  DedxSample->Branch("cell", &fCell, "cell/I");
264  DedxSample->Branch("plane", &fPlane, "plane/I");
265  DedxSample->Branch("adc", &fADC, "adc/F");
266  DedxSample->Branch("pe", &fPE, "pe/F");
267  DedxSample->Branch("pecorr", &fPECorr, "pecorr/F");
268  DedxSample->Branch("mip", &fMIP, "mip/F");
269  DedxSample->Branch("dx", &fDx, "dx/F");
270  DedxSample->Branch("deadx", &fxDead, "deadx/F");
271  DedxSample->Branch("xStart", &fxStart, "xStart/F");
272  DedxSample->Branch("xEnd", &fxEnd, "xEnd/F");
273  DedxSample->Branch("yStart", &fyStart, "yStart/F");
274  DedxSample->Branch("yEnd", &fyEnd, "yEnd/F");
275  DedxSample->Branch("xHit", &fxHit, "xHit/F");
276  DedxSample->Branch("yHit", &fyHit, "yHit/F");
277  DedxSample->Branch("zHit", &fzHit, "zHit/F");
278  DedxSample->Branch("tHit", &ftHit, "tHit/F");
279  DedxSample->Branch("xDir", &fxDir, "xDir/F");
280  DedxSample->Branch("yDir", &fyDir, "yDir/F");
281  DedxSample->Branch("zDir", &fzDir, "zDir/F");
282  DedxSample->Branch("mcde", &fmcDe, "mcde/F");
283  DedxSample->Branch("mcdeBirks",&fmcDeBirks, "mcdeBirks/F");
284  DedxSample->Branch("mcdx", &fmcDx, "mcdx/F");
285 
286  DedxSample->Branch("pdg", &fpdg, "pdg/I");
287  DedxSample->Branch("trueE", &fTrueE, "trueE/F");
288 
289  DedxSample->Branch("fiberBrightness", &fFiberBrightness, "fiberBrightness/I");
290 
291  DedxSample->Branch("multiCell", &fMultiCell );
292  DedxSample->Branch("multiADC", &fMultiADC );
293  DedxSample->Branch("multiPE", &fMultiPE );
294  DedxSample->Branch("multiPECorr", &fMultiPECorr);
295 
296 
297  fPOT = tfs->make<TH1D>("POT", ";;Total POT",1, 0.0, 1.0);
298  fEventsPassingSpillCuts = tfs->make<TH1D>("EventsPassingSpillCuts", ";;Events Passing Spill Cuts",1, 0.0, 1.0);
299  fSlicesPassingSelectionCuts = tfs->make<TH1D>("SlicesPassingSelectionCuts", ";;Slices Passing Selection Cuts",1, 0.0, 1.0);
300  fTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("TrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
301  fSuperTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("SuperTrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
302  fTrueProtonEFromTrackLength = tfs->make<TH2D>("TrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
303  fSuperTrueProtonEFromTrackLength = tfs->make<TH2D>("SuperTrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
304 
305  fTotalPOT = 0.0;
306 
307 
308 
309  }
310 
311  //......................................................................
312 
314 
315  {
316  fND = false;
317 
318  novadaq::cnv::DetId detID = fgeom->DetId();
319  if (detID == novadaq::cnv::kNEARDET) fND = true;
320 
321  return;
322  }
323 
324 
325  //.......................................................................
326 
328  {
329  //Only written for ND
330 
331  if (!fND) return;
332 
333  // get a geometry handle **************************************
336  //***************************************************************
337 
338  //Does this event pass spill cuts?
340  if(bt->HaveTruthInfo()){
341  evt.getByLabel(fGeneratorLabel, spillPot);
342  }
343  else{
344  evt.getByLabel(fNumiBeamLabel, spillPot);
345  }
346 
347  if (spillPot.failedToGet()){
348  mf::LogError("Susie") << "Spill Data not found for real data event";
349  return;
350  }
351 
352  //SL -- only cut on beam quality for data
353  if (!bt->HaveTruthInfo()){
354  if (std::abs(spillPot->deltaspilltimensec) > 0.5e9) return;
355  if ((spillPot->spillpot)*1e12 < 2e12) return;
356  if ((spillPot->hornI < -202)||(spillPot->hornI > -198)) return;
357  if ((spillPot->posx < -2.00)||(spillPot->posx > 2.00)) return;
358  if ((spillPot->posy < -2.00)||(spillPot->posy > 2.00)) return;
359  if ((spillPot->widthx < 0.57)||(spillPot->widthx > 1.58)) return;
360  if ((spillPot->widthy < 0.57)||(spillPot->widthy > 1.58)) return;
361  }
362 
363  //Data quality spill cuts
364 
365  // Get the EventQuality information
367  evt.getByLabel(fSpillQualLabel, spillQual);
368 
369  if(spillQual.failedToGet()){
370  mf::LogError("CAFMaker") << "Spill EventQuality not found for real data event";
371  return;
372  }
373 
374  if (spillQual->fracdcm3hits > 0.45) return;
375  if (spillQual->nmissingdcms > 0) return;
376 
377 
378  //Filling spill pot after it passes the cuts
379  float spillpot = spillPot->spillpot;
380  if (!bt->HaveTruthInfo()) spillpot *= 1e12;
381  fTotalPOT += spillpot;
382 
383  fEventsPassingSpillCuts->Fill(0.5);
384 
385  //Getting event hits
387  evt.getByLabel(fHitModuleLabel, hithdl);
388 
390  for(size_t h = 0; h < hithdl->size(); ++h){
391  art::Ptr<rb::CellHit> hit(hithdl, h);
392  allhits.push_back(hit);
393  }
394 
395 
396 
397  // get the slices out of the event
399  evt.getByLabel(fSliceLabel,slicecol);
400 
401 
402  art::PtrVector<rb::Cluster> slicelist;
403 
404  for(unsigned int i = 0; i<slicecol->size();++i){
405  art::Ptr<rb::Cluster>slice(slicecol,i);
406  slicelist.push_back(slice);
407  }
408 
409  // get assosciations between slices and tracks
410  art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
411 
412  // get assosciations between slices and numu energy
413  art::FindManyP<numue::NumuE> fndmnyNumuEnergy(slicecol, evt, fNumuEnergyLabel);
414  // get assosciations between slices and qepid
415  art::FindManyP<qeef::QePId> fndmnyQepid(slicecol, evt, fQepidLabel);
416 
417  //if (evt.event() == 443183){
418  // std::cout<<"Right event 2!"<<std::endl;
419  //}
420 
421  // loop over the slices
422  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
423 
424  // If slice is noise, skip it
425  if(slicelist[iSlice]->IsNoise()){ continue; }
426 
427  // get all of the tracks associated to this slice
428  std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
429  // if no tracks move on
430  if(tracks.size() == 0){ continue; }
431 
432  //Require precisely one Kalman tracks
433  if(tracks.size() != 1){ continue; }
434 
435  //if ((evt.event() == 443183)&&(iSlice == 2)){
436  // std::cout<<"Right event and slice!"<<std::endl;
437  //}
438 
439  //Require BOTH are 3D
440  bool all3D = true;
441 
442  // get the pid that belongs to this recotrack
443  art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);
444 
445  int bestidx = -1;
446  double bestpid = -1;
447  // loop over all of the tracks and find the most muon like one;
448  for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
449  if(!tracks[itrk]->Is3D()){
450  all3D = false;
451  continue;
452  }
453  // get the pid that belongs to this track
454  art::Ptr<remid::ReMId> pid = foids.at(itrk);
455  if(pid->Value() > bestpid){
456  bestpid = pid->Value();
457  bestidx = itrk;
458  }
459  }//End of loop over tracks
460 
461  if (!all3D){ continue; }
462 
463  //Default all pass cuts to false
464  fPassOffTrackHits = false;
465  fPassOffTrackEnergy = false;
466  fPassHadEOnMuonTrack = false;
467  fPassVertX = false;
468  fPassVertY = false;
469  fPassVertZ = false;
470  fPassProtonRemid = false;
471  fPassDedxRatio = false;
472  fPassEDiffZTest = false;
473 
474 
475  //Slightly harsh ReMId cut on muon track
476  if(bestidx == -1 || bestpid < 0.8){ continue; }
477 
478  if((bestidx != 0) && (bestidx !=1)){
479  std::cout<<"BUMMER! A weird track id for muon that shouldn't happen!!"<<bestidx<<std::endl;
480  continue;
481  }
482 
483  int protonidx = -1;
484  if (bestidx == 0){
485  protonidx = 1;
486  }
487  if (bestidx == 1){
488  protonidx = 0;
489  }
490 
491  if(protonidx == -1){
492  std::cout<<"BUMMER! A weird track id for proton that shouldn't happen!!"<<std::endl;
493  continue;
494  }
495 
496  //Filling muon track length properties
497  fzStart = tracks[bestidx]->Start().Z();
498  fzEnd = tracks[bestidx]->Stop().Z();
499  fTrackLength = tracks[bestidx]->TotalLength();
500 
501  // Track related Info
502  fxStart = tracks[bestidx]->Start().X();
503  fyStart = tracks[bestidx]->Start().Y();
504  fxEnd = tracks[bestidx]->Stop().X();
505  fyEnd = tracks[bestidx]->Stop().Y();
506 
507  fxDir = tracks[bestidx]->Dir().X();
508  fyDir = tracks[bestidx]->Dir().Y();
509  fzDir = tracks[bestidx]->Dir().Z();
510 
511  unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();
512  // loop over trajectory points starting from the end
513  for(int itrj = ntraj-1; itrj >= 0; --itrj){
514  TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
515 
516  } // for itrj
517 
518  unsigned int minPlane = tracks[bestidx]->MinPlane();
519  unsigned int maxPlane = tracks[bestidx]->MaxPlane();
520  int maxXPlane = -1;
521  int maxYPlane = -1;
522  double hitCosx[tracks[bestidx]->NXCell()];
523  double hitCosy[tracks[bestidx]->NYCell()];
524 
525  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NXCell(); ++ihit){
526  if(tracks[bestidx]->XCell(ihit)->Plane() > maxXPlane){ maxXPlane = tracks[bestidx]->XCell(ihit)->Plane(); }
527  // get direction with this x hit
528  double xyzplane[3];
529  geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
530  double dz = geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->HalfD();
531 
532  double minz = xyzplane[2]-dz;
533  double maxz = xyzplane[2]+dz;
534  unsigned int currtrj = 0;
535  // find the matching trajectory point for this plane;
536  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
537  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
538  } // for itaraj
539  TVector3 direction = tracks[bestidx]->Dir();
540  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
541  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
542  }
543  direction = direction.Unit();
544  hitCosx[ihit] = sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
545  } // for ihit XCell
546 
547 
548 
549  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NYCell(); ++ihit){
550  if(tracks[bestidx]->YCell(ihit)->Plane() > maxYPlane){ maxYPlane = tracks[bestidx]->YCell(ihit)->Plane(); }
551  // get direction with this y hit
552  double xyzplane[3];
553  geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
554  double dz = geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->HalfD();
555  // double dz = HalfD;
556  double minz = xyzplane[2]-dz;
557  double maxz = xyzplane[2]+dz;
558  unsigned int currtrj = 0;
559  // find the matching trajectory point for this plane;
560  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
561  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
562  }
563  TVector3 direction = tracks[bestidx]->Dir();
564  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
565  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
566  }
567  direction = direction.Unit();
568  hitCosy[ihit] = sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
569  } // for ihit YCell
570 
571  ///
572 
573  // calculate dedx info for this track
574  std::vector<double> dedxs;
575  std::vector<double> des;
576  std::vector<double> pecorrs;
577  std::vector<double> mips;
578  std::vector<double> dxs;
579  std::vector<double> deads;
580  std::vector<double> xdedxs;
581  std::vector<unsigned int> dedxpls;
582 
583  std::vector<double> pes;
584  std::vector<double> adcs;
585  std::vector<double> xhits;
586  std::vector<double> yhits;
587  std::vector<double> zhits;
588  std::vector<double> thits;
589  std::vector<int> views;
590  std::vector<int> cells;
591  std::vector<int> planes;
592 
593  std::vector<double> mcdes;
594  std::vector<double> mcdebs;
595  std::vector<double> mcdxs;
596 
597  std::vector< std::vector<int> > mcells;
598  std::vector< std::vector<float> > madcs;
599  std::vector< std::vector<float> > mpes;
600  std::vector< std::vector<float> > mpecorrs;
601 
602  for(unsigned int iPlane = minPlane;iPlane<=maxPlane;++iPlane){
603  double planeGeV = 0;
604  double planePECorr = 0;
605  double planeMIP = 0;
606 
607  int planeView = -1;
608  int plane = -1;
609  int cell = -1;
610  double planeX = 0;
611  double planeY = 0;
612  double planeZ = 0;
613  double planeT = 0;
614  double planeADC = 0;
615  double planePE = 0;
616 
617  std::vector<int> mCell;
618  std::vector<float> mPE;
619  std::vector<float> mADC;
620  std::vector<float> mPECorr;
621 
622  double mcPlaneE = 0;
623  double mcPlaneEBirks = 0;
624  double mcPathLength = 0;
625  double totcosinview = 0;
626  double avgCosOtherView = 1.0;
627  bool useOtherView = true;
628  double nhits = 0;
629  double xyzplane[3];
630  geom->Plane(iPlane)->Cell(0)->GetCenter(xyzplane);
631  // double dz = geom->Plane(iPlane)->Cell(0)->HalfD();
632  // double minz = xyzplane[2]-dz;
633  // double maxz = xyzplane[2]+dz;
634  // int matchpt = 0; // Not used for now
635  int minCell = 10000;
636  int maxCell = 0;
637  double nhitsother = 0;
638  int view = 0;
639  int fbIndex = -1;
640  fFiberBrightness = -1;
641  // find the matching trajectory point for this plane;
642  /* Not used for now
643  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
644  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){matchpt = itraj;}
645  }
646  */
647 
648  if(geom->Plane(iPlane)->View() == geo::kX){
649  planeView = 0;
650  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
651  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane && tracks[bestidx]->XCell(iHit)->View() != geo::kY){
652  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
653  planeADC+=chit->ADC();
654  planePE +=chit->PE();
655  cell = chit->Cell();
656  plane = chit->Plane();
657 
658  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->XCell(iHit));
659  if(minPlane == maxPlane){
660  rhit = cal->MakeRecoHit(*(tracks[bestidx]->XCell(iHit)), 0.25);
661  }
662  if(rhit.IsCalibrated()){
663  planeGeV+=rhit.GeV();
664  planePECorr+=rhit.PECorr();
665  planeMIP+=rhit.MIP();
666  planeX += rhit.X();
667  planeY += rhit.Y();
668  planeZ += rhit.Z();
669  planeT += rhit.T();
670  totcosinview+=hitCosx[iHit];
671  mCell.push_back( chit->Plane() );
672  mPE.push_back( chit->PE() );
673  mADC.push_back( chit->ADC() );
674  mPECorr.push_back(rhit.PECorr() );
675 
676  nhits+=1.0;
677  if(tracks[bestidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->XCell(iHit)->Cell();}
678  if(tracks[bestidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->XCell(iHit)->Cell();}
679 
680  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->XCell(iHit)->Cell() );
681 
682  if(bt->HaveTruthInfo()){
683  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
684  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
685  for(int j = 0; j < (int)flsHits.size(); j++){
686  mcPlaneE += flsHits[j].GetEdep();
687  mcPlaneEBirks += flsHits[j].GetEdepBirks();
688  mcPathLength += flsHits[j].GetTotalPathLength();
689 
690  }//End loop over flshits in hit
691  } // if BackTracker
692 
693  }
694  }
695  } // for iHit
696 
697  // try to get the avgcos from the other view
698  double nhitsbefore = 0;
699  double nhitsafter = 0;
700  double totCosBefore = 0;
701  double totCosAfter = 0;
702  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
703  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane-1){
704  ++nhitsbefore;
705  totCosBefore+=hitCosy[iHit];
706  }
707  else if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane+1){
708  ++nhitsafter;
709  totCosAfter+=hitCosy[iHit];
710  }
711  }
712  if(nhitsbefore > 0 && nhitsafter > 0){
713  avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
714  nhitsother = (nhitsbefore+nhitsafter)/2;
715  }
716  else if(nhitsbefore > 0){
717  avgCosOtherView = totCosBefore/nhitsbefore;
718  nhitsother = nhitsbefore;
719  }
720  else if(nhitsafter > 0){
721  avgCosOtherView = totCosAfter/nhitsafter;
722  nhitsother = nhitsafter;
723  }
724  else{ useOtherView = false; }
725 
726  } else if(geom->Plane(iPlane)->View() == geo::kY){
727  planeView = 1;
728  view = 1;
729  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
730  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane && tracks[bestidx]->YCell(iHit)->View() != geo::kX
731  && tracks[bestidx]->YCell(iHit)->View() != geo::kXorY){
732  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
733  planeADC+=chit->ADC();
734  planePE +=chit->PE();
735  cell = chit->Cell();
736  plane = chit->Plane();
737 
738  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->YCell(iHit));
739  if(minPlane == maxPlane){
740  rhit = cal->MakeRecoHit(*(tracks[bestidx]->YCell(iHit)), 0.25);
741  }
742  if(rhit.IsCalibrated()){
743  planeGeV+=rhit.GeV();
744  planePECorr+=rhit.PECorr();
745  planeMIP+=rhit.MIP();
746  planeX += rhit.X();
747  planeY += rhit.Y();
748  planeZ += rhit.Z();
749  planeT += rhit.T();
750  totcosinview+=hitCosy[iHit];
751 
752  mCell.push_back( chit->Plane() );
753  mPE.push_back( chit->PE() );
754  mADC.push_back( chit->ADC() );
755  mPECorr.push_back(rhit.PECorr() );
756 
757  nhits+=1.0;
758  if(tracks[bestidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->YCell(iHit)->Cell();}
759  if(tracks[bestidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->YCell(iHit)->Cell();}
760 
761  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->YCell(iHit)->Cell() );
762 
763  if(bt->HaveTruthInfo()){
764  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
765  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
766  for(int j = 0; j < (int)flsHits.size(); j++){
767  mcPlaneE += flsHits[j].GetEdep();
768  mcPlaneEBirks += flsHits[j].GetEdepBirks();
769  mcPathLength += flsHits[j].GetTotalPathLength();
770 
771  }//End loop over flshits in hit
772  } // if BackTracker
773 
774  } // if isCalibrated
775  } //
776  } // for iHit
777 
778  // try to get the avgcos from the other view
779  double nhitsbefore = 0;
780  double nhitsafter = 0;
781  double totCosBefore = 0;
782  double totCosAfter = 0;
783  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
784  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane-1){
785  ++nhitsbefore;
786  totCosBefore+=hitCosx[iHit];
787  }
788  else if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane+1){
789  ++nhitsafter;
790  totCosAfter+=hitCosx[iHit];
791  }
792  }
793  if(nhitsbefore > 0 && nhitsafter > 0){
794  avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
795  nhitsother = (nhitsbefore+nhitsafter)/2;
796  }
797  else if(nhitsbefore > 0){
798  avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
799  nhitsother = nhitsbefore;
800  }
801  else if(nhitsafter > 0){
802  avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
803  nhitsother = nhitsbefore;
804  }
805  else{ useOtherView = false; }
806 
807  } // if kY
808  fnCellsPerPlane = nhits;
809  if(nhits >0){
810 
811  if(nhits == 1) {
812  fFiberBrightness = fbIndex;
813  }
814 
815  double avgCos = TMath::Abs(totcosinview/nhits);
816  double xyzp[3];
817  geom->Plane(iPlane)->Cell(minCell)->GetCenter(xyzp);
818  double xyzp1[3];
819  geom->Plane(iPlane)->Cell(maxCell)->GetCenter(xyzp1);
820  double planewidth = 2*geom->Plane(iPlane)->Cell(1)->HalfD();
821  // double planewidth = 2*HalfD;
822  // calculate the amount of dead material to subtract off
823  // first add up live material
824  double liveLength = 0;
825  for(int i = minCell; i<maxCell+1;++i){
826  // get the cell half width
827  const geo::PlaneGeo* p = geom->Plane(iPlane);
828  const geo::CellGeo* c = p->Cell(i);
829  if(i == minCell || i == maxCell){ liveLength+=c->HalfW(); }
830  else{ liveLength+=(2.0*c->HalfW() ); }
831  }
832  // dead material is total length - live length only if min and max cell are not the same
833  double deadLength = 0;
834  if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
835  double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
836  double planePathLength = TMath::Sqrt(planewidth*planewidth + planeheight*planeheight);
837 
838  // in the case of completely vertical track be a bit careful
839  if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }
840  if(useOtherView){
841  double deltaz = planewidth;
842  double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
843  if(avgCos == 0 || avgCos!=avgCos){
844  deltav = liveLength;
845  // if the track looks completely vertical in both views,
846  // than no part of the track length should be in z direction
847  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
848  }
849 
850  // in opposite view dead length
851  double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
852  // in the case that the other view looks completely vertical scale up
853  // the length between cells in this view by the number of hits in the other view
854  // this case should really never happen
855  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
856  deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
857  }
858  planePathLength = TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
859  }
860  // also the case of completely vertical track
861  if(minPlane == maxPlane){ planePathLength = liveLength; }
862 
863  assert(planeGeV/(planePathLength)>=0);
864  dedxs.push_back(planeGeV/(planePathLength));
865  des.push_back(planeGeV);
866  pecorrs.push_back(planePECorr);
867  mips.push_back(planeMIP);
868  dxs.push_back(planePathLength);
869  deads.push_back(deadLength);
870  xhits.push_back(planeX/nhits);
871  yhits.push_back(planeY/nhits);
872  zhits.push_back(planeZ/nhits);
873  thits.push_back(planeT/nhits);
874  views.push_back(planeView);
875  planes.push_back(plane);
876  cells.push_back(cell);
877  pes.push_back(planePE);
878  adcs.push_back(planeADC);
879 
880  mcdes.push_back(mcPlaneE);
881  mcdebs.push_back(mcPlaneEBirks);
882  mcdxs.push_back(mcPathLength);
883 
884  mcells.push_back(mCell);
885  madcs.push_back(mADC);
886  mpes.push_back(mPE);
887  mpecorrs.push_back(mPECorr);
888 
889  mcdes.push_back(mcPlaneE);
890  mcdebs.push_back(mcPlaneEBirks);
891  mcdxs.push_back(mcPathLength);
892 
893  } //nhits > 0
894  } // iPlane
895 
896  // ..... .......... .......... .......... .......... .......... .....
897 
898  art::Ptr<remid::ReMId> pid = foids.at(bestidx);
899 
900  //Looking at numu energy things
901  numue::NumuE numuE;
902  if (!fndmnyNumuEnergy.isValid()){ continue; }
903  std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
904  if (energyVec.empty()){ continue; }
905  numuE = *energyVec[0];
906 
907  //Require decent energy object; not too much hadronic energy in muon catcher
908  if (numuE.E() == -5){ continue; }
909  if (!(numuE.NDHadCalCat() + numuE.NDHadCalTran() < 0.03)){ continue; }
910  //Not much hadronic energy contamination on muon track
911  if (numuE.HadTrkE() < 0.025){ fPassHadEOnMuonTrack = true; }
912 
913  //Looking at qepid things
914  qeef::QePId qepid;
915  if (!fndmnyQepid.isValid()){ continue; }
916  std::vector<art::Ptr<qeef::QePId> > qVec = fndmnyQepid.at(iSlice);
917  if (qVec.empty()){ continue; }
918  qepid = *qVec[0];
919 
920  //Require ratio of average dE/dx track is good
921  if (qepid.Dedx() > 1.75){ fPassDedxRatio = true; }
922  //Require Ztest of two methods of calculating QE neutrino energy is good
923  if (qepid.EdiffZ() > -1.3){ fPassEDiffZTest = true; }
924 
925 
926  //Slice energy mostly on one tracks
927  double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy());
928  if (calEDiff < 0.1){ fPassOffTrackEnergy = true; }
929 
930  //Slice hits mostly on two tracks
931  int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell());
932  //if (!(hitDiff < 4)){ continue; }
933  if (hitDiff < 4){ fPassOffTrackHits = true; }
934 
935  //Check if true signal for sim events
936  fTrueNumuCCQE = false;
937  fTrueNumuCCWithProton = false;
939 
940  fpdg = 0;
941 
942  if (bt->HaveTruthInfo()){
943  std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),allhits);
944  if (!funcReturn.empty()){
945  // const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();
946 
947  const std::vector<const sim::Particle*> test_mu = bt->HitsToParticle(tracks[bestidx]->AllCells());
948  if (!test_mu.empty()) {
949  fpdg = test_mu[0]->PdgCode();
950  fTrueE = test_mu[0]->E();
951  // std::cout << "Mu " << fTrueE << " " << fpdg << std::endl;
952  }
953  }
954  }
955 
956  // if(ContainedEvent(slicelist[iSlice],tracks[bestidx])){
957  if(true ){ // Containment cut removed for now. Always pass
959  fSlicesPassingSelectionCuts->Fill(0.5);
962  }
963  }
964  // get all the dedx variables filled
965  fnDedx = pid->NDedx()-1;
967  // loop over all of the dE/dx measurements
968  for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
969  fDedx = pid->Dedx(idedx);
970  fVert = pid->DedxVertex(idedx);
971  fUsed = pid->DedxUsed(idedx);
972 
973  fDe = des[idedx];
974  fPECorr = pecorrs[idedx];
975  fMIP = mips[idedx];
976  fDx = dxs[idedx];
977  fxDead = deads[idedx];
978  fmyDedx = fDe/fDx;
979 
980  fxHit = xhits[idedx];
981  fyHit = yhits[idedx];
982  fzHit = zhits[idedx];
983  ftHit = thits[idedx];
984  fView = views[idedx];
985  fPlane = planes[idedx];
986  fCell = cells[idedx];
987  fADC = adcs[idedx];
988  fPE = pes[idedx];
989 
990  fMultiCell = mcells[idedx];
991  fMultiADC = madcs[idedx];
992  fMultiPE = mpes[idedx];
993  fMultiPECorr = mpecorrs[idedx];
994 
995  fmcDe = mcdes[idedx];
996  fmcDeBirks = mcdebs[idedx];
997  fmcDx = mcdxs[idedx];
998  // get the distance from the end of the track
999  fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[bestidx]->TotalLength();
1000  // fill the dedx sample tree
1001  DedxSample->Fill();
1002  }//Loop over DedxSample
1003  }//If passes contaiment
1004  } // end loop over slices
1005 
1006  }
1007 
1008 
1009  //.......................................................................
1010 
1012  {
1013  fPOT->Fill(.5, fTotalPOT);
1014  }
1015 
1016 
1017  //.......................................................................
1019  {
1020 
1021  bool contained = false;
1022 
1023  // check that this is a contained event
1024 
1025  // get some useful information here
1026  // 1) how close is the slice to the edge of the detector
1027  unsigned int slcncellsfromedge = 999999999;
1028  for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
1029  const art::Ptr<rb::CellHit>& chit = slice->Cell(ihit);
1030  const int plane = chit->Plane();
1031  unsigned int cellsfromedge = std::min((unsigned int)chit->Cell(), fgeom->Plane(plane)->Ncells() - 1 - chit->Cell());
1032  if(cellsfromedge < slcncellsfromedge){
1033  slcncellsfromedge = cellsfromedge;
1034  }
1035  }
1036  // 2) how close is the track to the edge of the detector
1037  // in planes
1038  //SL - NOPE, need to ask slice, not track
1039  unsigned int firstpl = slice->MinPlane();
1040  unsigned int lastpl = slice->MaxPlane();
1041  // in projected cells, just define directions here actually get projected cell below
1042  TVector3 dirbwk = -trk->Dir();
1043  TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();
1044 
1045  // 3) hits in slice
1046  int sliceNHit = slice->NCell();
1047 
1048  // 4) slice mean time
1049  double slcTime = slice->MeanTNS();
1050 
1051  // 5) hits on track
1052  int trackNHit = trk->NCell();
1053 
1054  //std::cout<<"In contained event loop!"<<std::endl;
1055  // cut the slice number of cells from edge
1056  if((sliceNHit > 20)&&(slcncellsfromedge > 1)&&(slcTime>50000)&&(slcTime<450000)&&(trackNHit > 5)){
1057  if(firstpl > 1 && lastpl < 212){
1058  //std::cout<<"In contained event loop! 1"<<std::endl;
1059  int fwdcell = flivegeom->FullNDProjectedCells(trk->Stop(), dirfwd);
1060  int bwkcell = flivegeom->FullNDProjectedCells(trk->Start(),dirbwk);
1061  // also make sure the track doesn't start too close to the end of the active region
1062  if(fwdcell > 4 && bwkcell > 8 && trk->Start().Z() < 1150){
1063  // if the track ends in the muon catcher, check that it doesn't clip the corner
1064  //std::cout<<"In contained event loop! 2"<<std::endl;
1065  if(trk->Stop().Z() < 1275){
1066  //std::cout<<"In contained event loop! 2a"<<std::endl;
1067  contained = true;
1068  }
1069  else{
1070  //std::cout<<"In contained event loop! 2b"<<std::endl;
1071  /*
1072  // find the trajectory point closest in z to the transition plane
1073  double mindist = 99999999;
1074  unsigned int mintrajid = 0;
1075  for(unsigned int itraj = 0; itraj < trk->NTrajectoryPoints()-1; ++itraj){
1076  double dist = 1275 - trk->TrajectoryPoint(itraj).Z();
1077  if(dist > 0 && dist < mindist){
1078  mindist = dist;
1079  mintrajid = itraj;
1080  }
1081  }
1082  // get the direction at this trajectory point
1083  TVector3 trandir = trk->Dir();
1084  if(mintrajid < trk->NTrajectoryPoints() -2){
1085  TVector3 trandir = (trk->TrajectoryPoint(mintrajid+1) - trk->TrajectoryPoint(mintrajid)).Unit();
1086  }
1087  */
1088  //std::cout<<"In contained event loop! 3"<<std::endl;
1089  //double ypostran = flivegeom->YPositionAtTransition(trk->TrajectoryPoint(mintrajid),trandir);
1090  double ypostranKirk = flivegeom->YPositionAtTransition(trk->Start(),trk->Dir());
1091  if( ypostranKirk < 55){
1092  //std::cout<<"In contained event loop! 4"<<std::endl;
1093  contained = true;
1094  }
1095 
1096  } // end else in muon catcher
1097  } // end if track ends projections contained
1098  } // end if track end planes contained
1099  // } // end if expected number of diblocks
1100  } // end if slice contained
1101 
1102 
1103  return contained;
1104 
1105  }
1106 
1107 
1108 
1109 
1110 }
1111 
1112 
1113 namespace remid{
1114 
1116 
1117 }
float T() const
Definition: RecoHit.h:39
size_t NTrajectoryPoints() const
Definition: Track.h:83
double EdiffZ() const
Definition: QePId.cxx:93
back track the reconstruction to the simulation
double Dedx(unsigned int i) const
Definition: ReMId.cxx:126
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
float E() const
Definition: Energy.cxx:27
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
unsigned short Plane() const
Definition: CellHit.h:39
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
float Z() const
Definition: RecoHit.h:38
double HalfW() const
Definition: CellGeo.cxx:191
art::ServiceHandle< photrans::FiberBrightness > fb
void beginRun(const art::Run &run)
art::ServiceHandle< geo::Geometry > fgeom
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
art::ServiceHandle< cheat::BackTracker > bt
virtual TVector3 Start() const
Definition: Prong.h:73
int nmissingdcms
of missing DCMs
Definition: EventQuality.h:23
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
Definition: Run.h:31
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
bool DedxVertex(unsigned int i) const
Definition: ReMId.cxx:138
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
signed long long int deltaspilltimensec
Definition: SpillData.h:24
double fracdcm3hits
fraction of DCM3 hits in horizontal modules
Definition: EventQuality.h:24
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
int NMeasurementPlanes() const
Definition: ReMId.cxx:223
float HadTrkE() const
Definition: NumuE.cxx:284
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
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
double widthx
mm
Definition: SpillData.h:36
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
int evt
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
float PE() const
Definition: CellHit.h:42
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
bool DedxUsed(unsigned int i) const
Definition: ReMId.cxx:132
std::string fPidLabel
Where to find the pids.
const double j
Definition: BetheBloch.cxx:29
double hornI
kA
Definition: SpillData.h:27
std::vector< int > fMultiCell
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double posy
mm
Definition: SpillData.h:35
ReMIdDedxRock(fhicl::ParameterSet const &pset)
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
float NDHadCalTran() const
Definition: NumuE.cxx:334
OStream cout
Definition: OStream.cxx:6
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
float MIP() const
Definition: RecoHit.cxx:58
std::vector< float > fMultiADC
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
float GeV() const
Definition: RecoHit.cxx:69
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
double DedxLocation(unsigned int i) const
Definition: ReMId.cxx:120
art::ServiceHandle< geo::LiveGeometry > flivegeom
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
Definition: structs.h:12
std::string fSliceLabel
Where to find reconstructed slices.
float X() const
Definition: RecoHit.h:36
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
std::string fTrackLabel
Where to find recondtructed tracks.
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
std::vector< float > fMultiPECorr
void analyze(art::Event const &evt)
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
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
float Y() const
Definition: RecoHit.h:37
float PECorr() const
Definition: RecoHit.cxx:47
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Encapsulate the cell geometry.
Definition: CellGeo.h:25
TH1D * fPOT
Histogram of POT passing spill cuts.
double posx
mm
Definition: SpillData.h:34
bool fND
Is detector ND?
float NDHadCalCat() const
Definition: NumuE.cxx:339
double widthy
mm
Definition: SpillData.h:37
double spillpot
POT for spill normalized by 10^12.
Definition: SpillData.h:26
std::vector< float > fMultiPE
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.