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"
19 #include "art_root_io/TFileService.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;
147 
148  // variables added to trees for detailed dE/dx studies
150  float fmyDedx;
151  float fDe;
152  int fView;
153  int fPlane;
154  int fCell;
155  float fPE;
156  float fADC;
157  float fPECorr;
158  float fMIP;
159  float fDx;
160  float fxDead;
161  float fxStart;
162  float fxEnd;
163  float fyStart;
164  float fyEnd;
165  float fxHit;
166  float fyHit;
167  float fzHit;
168  float ftHit;
169  float fxDir;
170  float fyDir;
171  float fzDir;
172  float fmcDe;
173  float fmcDeBirks;
174  float fmcDx;
175 
176  float fTrueE;
178 
179  std::vector <int> fMultiCell;
180  std::vector <float> fMultiPE;
181  std::vector <float> fMultiADC;
182  std::vector <float> fMultiPECorr;
183 
184  // ......................
185 
190 
192 
193  //novadaq::cnv::DetId fDetID; ///< Detector ID in nova daq convention namespace typedef
194 
195 
196  //unsigned int fdibfirst = 0;
197  //unsigned int fdiblast = 0;
198 
200 
201  };
202 
203 }
204 
205 
206 namespace remid{
207 
208 
209  //.......................................................................
210 
212  EDAnalyzer(pset)
213  {
214  fHitModuleLabel = (pset.get< std::string >("HitModuleLabel"));
215  fSliceLabel = (pset.get< std::string >("SliceLabel"));
216  fTrackLabel = (pset.get< std::string >("TrackLabel"));
217  fPidLabel = (pset.get< std::string >("PidLabel" ));
218  fGeneratorLabel = pset.get< std::string >("GeneratorInput");
219  fNumiBeamLabel = pset.get< std::string >("NuMIBeamInput");
220  fSpillQualLabel = pset.get< std::string >("SpillQualInput");
221  fNumuEnergyLabel = pset.get< std::string >("NumuEnergyInput");
222  fQepidLabel = pset.get< std::string >("QepidInput");
223  }
224 
225  //.......................................................................
226 
228  {
229 
230  }
231 
232  //.......................................................................
233 
235  {
237 
238  DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
239  DedxSample->Branch("dedx", &fDedx, "dedx/F");
240  DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
241  DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
242  DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
243  DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
244  DedxSample->Branch("vert", &fVert, "vert/O");
245  DedxSample->Branch("used", &fUsed, "used/O");
246  DedxSample->Branch("zStart",&fzStart,"zStart/F");
247  DedxSample->Branch("zEnd",&fzEnd,"zEnd/F");
248  DedxSample->Branch("passOffTrackHits",&fPassOffTrackHits,"passOffTrackHits/O");
249  DedxSample->Branch("passOffTrackEnergy",&fPassOffTrackEnergy,"passOffTrackEnergy/O");
250  DedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");
251  DedxSample->Branch("passRockMuonCut",&fPassRockMuonCut,"passRockMuonCut/O");
252  DedxSample->Branch("passVertX",&fPassVertX,"passVertX/O");
253  DedxSample->Branch("passVertY",&fPassVertY,"passVertY/O");
254  DedxSample->Branch("passVertZ",&fPassVertZ,"passVertZ/O");
255  DedxSample->Branch("passProtonRemid",&fPassProtonRemid,"passProtonRemid/O");
256  DedxSample->Branch("passDedxRatio",&fPassDedxRatio,"passDedxRatio/O");
257  DedxSample->Branch("passEDiffZTest",&fPassEDiffZTest,"passEDiffZTestX/O");
258  DedxSample->Branch("trueNumuCCQE",&fTrueNumuCCQE,"trueNumuCCQE/O");
259 
260  // Branches added wrt standard ReMIdDedx module for detailed dE/dx studies
261  DedxSample->Branch("ncellsPerPlane", &fnCellsPerPlane, "ncellsPerPlane/I");
262  DedxSample->Branch("mydedx", &fmyDedx, "mydedx/F");
263  DedxSample->Branch("de", &fDe, "de/F");
264  DedxSample->Branch("view", &fView, "view/I");
265  DedxSample->Branch("cell", &fCell, "cell/I");
266  DedxSample->Branch("plane", &fPlane, "plane/I");
267  DedxSample->Branch("adc", &fADC, "adc/F");
268  DedxSample->Branch("pe", &fPE, "pe/F");
269  DedxSample->Branch("pecorr", &fPECorr, "pecorr/F");
270  DedxSample->Branch("mip", &fMIP, "mip/F");
271  DedxSample->Branch("dx", &fDx, "dx/F");
272  DedxSample->Branch("deadx", &fxDead, "deadx/F");
273  DedxSample->Branch("xStart", &fxStart, "xStart/F");
274  DedxSample->Branch("xEnd", &fxEnd, "xEnd/F");
275  DedxSample->Branch("yStart", &fyStart, "yStart/F");
276  DedxSample->Branch("yEnd", &fyEnd, "yEnd/F");
277  DedxSample->Branch("xHit", &fxHit, "xHit/F");
278  DedxSample->Branch("yHit", &fyHit, "yHit/F");
279  DedxSample->Branch("zHit", &fzHit, "zHit/F");
280  DedxSample->Branch("tHit", &ftHit, "tHit/F");
281  DedxSample->Branch("xDir", &fxDir, "xDir/F");
282  DedxSample->Branch("yDir", &fyDir, "yDir/F");
283  DedxSample->Branch("zDir", &fzDir, "zDir/F");
284  DedxSample->Branch("mcde", &fmcDe, "mcde/F");
285  DedxSample->Branch("mcdeBirks",&fmcDeBirks, "mcdeBirks/F");
286  DedxSample->Branch("mcdx", &fmcDx, "mcdx/F");
287 
288  DedxSample->Branch("pdg", &fpdg, "pdg/I");
289  DedxSample->Branch("trueE", &fTrueE, "trueE/F");
290 
291  DedxSample->Branch("fiberBrightness", &fFiberBrightness, "fiberBrightness/I");
292 
293  DedxSample->Branch("multiCell", &fMultiCell );
294  DedxSample->Branch("multiADC", &fMultiADC );
295  DedxSample->Branch("multiPE", &fMultiPE );
296  DedxSample->Branch("multiPECorr", &fMultiPECorr);
297 
298 
299  fPOT = tfs->make<TH1D>("POT", ";;Total POT",1, 0.0, 1.0);
300  fEventsPassingSpillCuts = tfs->make<TH1D>("EventsPassingSpillCuts", ";;Events Passing Spill Cuts",1, 0.0, 1.0);
301  fSlicesPassingSelectionCuts = tfs->make<TH1D>("SlicesPassingSelectionCuts", ";;Slices Passing Selection Cuts",1, 0.0, 1.0);
302  fTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("TrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
303  fSuperTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("SuperTrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
304  fTrueProtonEFromTrackLength = tfs->make<TH2D>("TrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
305  fSuperTrueProtonEFromTrackLength = tfs->make<TH2D>("SuperTrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
306 
307  fTotalPOT = 0.0;
308 
309 
310 
311  }
312 
313  //......................................................................
314 
316 
317  {
318  fND = false;
319 
320  novadaq::cnv::DetId detID = fgeom->DetId();
321  if (detID == novadaq::cnv::kNEARDET) fND = true;
322 
323  return;
324  }
325 
326 
327  //.......................................................................
328 
330  {
331  //Only written for ND
332 
333  if (!fND) return;
334 
335  // get a geometry handle **************************************
338  //***************************************************************
339 
340  //Does this event pass spill cuts?
342  if(bt->HaveTruthInfo()){
343  evt.getByLabel(fGeneratorLabel, spillPot);
344  }
345  else{
346  evt.getByLabel(fNumiBeamLabel, spillPot);
347  }
348 
349  if (spillPot.failedToGet()){
350  mf::LogError("Susie") << "Spill Data not found for real data event";
351  return;
352  }
353 
354  //SL -- only cut on beam quality for data
355  if (!bt->HaveTruthInfo()){
356  if (std::abs(spillPot->deltaspilltimensec) > 0.5e9) return;
357  if ((spillPot->spillpot)*1e12 < 2e12) return;
358  if ((spillPot->hornI < -202)||(spillPot->hornI > -198)) return;
359  if ((spillPot->posx < -2.00)||(spillPot->posx > 2.00)) return;
360  if ((spillPot->posy < -2.00)||(spillPot->posy > 2.00)) return;
361  if ((spillPot->widthx < 0.57)||(spillPot->widthx > 1.58)) return;
362  if ((spillPot->widthy < 0.57)||(spillPot->widthy > 1.58)) return;
363  }
364 
365  //Data quality spill cuts
366 
367  // Get the EventQuality information
369  evt.getByLabel(fSpillQualLabel, spillQual);
370 
371  if(spillQual.failedToGet()){
372  mf::LogError("CAFMaker") << "Spill EventQuality not found for real data event";
373  return;
374  }
375 
376  if (spillQual->fracdcm3hits > 0.45) return;
377  if (spillQual->nmissingdcms > 0) return;
378 
379 
380  //Filling spill pot after it passes the cuts
381  float spillpot = spillPot->spillpot;
382  if (!bt->HaveTruthInfo()) spillpot *= 1e12;
383  fTotalPOT += spillpot;
384 
385  fEventsPassingSpillCuts->Fill(0.5);
386 
387  //Getting event hits
389  evt.getByLabel(fHitModuleLabel, hithdl);
390 
392  for(size_t h = 0; h < hithdl->size(); ++h){
393  art::Ptr<rb::CellHit> hit(hithdl, h);
394  allhits.push_back(hit);
395  }
396 
397 
398 
399  // get the slices out of the event
401  evt.getByLabel(fSliceLabel,slicecol);
402 
403 
404  art::PtrVector<rb::Cluster> slicelist;
405 
406  for(unsigned int i = 0; i<slicecol->size();++i){
407  art::Ptr<rb::Cluster>slice(slicecol,i);
408  slicelist.push_back(slice);
409  }
410 
411  // get assosciations between slices and tracks
412  art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
413 
414  // get assosciations between slices and numu energy
415  art::FindManyP<numue::NumuE> fndmnyNumuEnergy(slicecol, evt, fNumuEnergyLabel);
416  // get assosciations between slices and qepid
417  // ecm art::FindManyP<qeef::QePId> fndmnyQepid(slicecol, evt, fQepidLabel);
418 
419  //if (evt.event() == 443183){
420  // std::cout<<"Right event 2!"<<std::endl;
421  //}
422 
423  // loop over the slices
424  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
425 
426  // If slice is noise, skip it
427  if(slicelist[iSlice]->IsNoise()){ continue; }
428 
429  // get all of the tracks associated to this slice
430  std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
431  // if no tracks move on
432  if(tracks.size() == 0){ continue; }
433 
434  //Require precisely one Kalman tracks
435  if(tracks.size() != 1){ continue; }
436 
437  //if ((evt.event() == 443183)&&(iSlice == 2)){
438  // std::cout<<"Right event and slice!"<<std::endl;
439  //}
440 
441  //Require BOTH are 3D
442  bool all3D = true;
443 
444  // get the pid that belongs to this recotrack
445  art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);
446 
447  int bestidx = -1;
448  double bestpid = -1;
449  // loop over all of the tracks and find the most muon like one;
450  for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
451  if(!tracks[itrk]->Is3D()){
452  all3D = false;
453  continue;
454  }
455  // get the pid that belongs to this track
456  art::Ptr<remid::ReMId> pid = foids.at(itrk);
457  if(pid->Value() > bestpid){
458  bestpid = pid->Value();
459  bestidx = itrk;
460  }
461  }//End of loop over tracks
462 
463  if (!all3D){ continue; }
464 
465  //Default all pass cuts to false
466  fPassOffTrackHits = false;
467  fPassOffTrackEnergy = false;
468  fPassHadEOnMuonTrack = false;
469 
470  //ECM I don't think these make sense anymore
471  fPassRockMuonCut = false;
472  fPassVertX = true;
473  fPassVertY = true;
474  fPassVertZ = true;
475  fPassProtonRemid = true;
476  fPassDedxRatio = true;
477  fPassEDiffZTest = true;
478 
479 
480  //Slightly harsh ReMId cut on muon track
481  if(bestidx == -1 || bestpid < 0.8){ continue; }
482 
483  if((bestidx != 0) && (bestidx !=1)){
484  std::cout<<"BUMMER! A weird track id for muon that shouldn't happen!!"<<bestidx<<std::endl;
485  continue;
486  }
487 
488  int protonidx = -1;
489  if (bestidx == 0){
490  protonidx = 1;
491  }
492  if (bestidx == 1){
493  protonidx = 0;
494  }
495 
496  if(protonidx == -1){
497  std::cout<<"BUMMER! A weird track id for proton that shouldn't happen!!"<<std::endl;
498  continue;
499  }
500 
501  //Filling muon track length properties
502  fzStart = tracks[bestidx]->Start().Z();
503  fzEnd = tracks[bestidx]->Stop().Z();
504  fTrackLength = tracks[bestidx]->TotalLength();
505 
506  // Track related Info
507  fxStart = tracks[bestidx]->Start().X();
508  fyStart = tracks[bestidx]->Start().Y();
509  fxEnd = tracks[bestidx]->Stop().X();
510  fyEnd = tracks[bestidx]->Stop().Y();
511 
512  fxDir = tracks[bestidx]->Dir().X();
513  fyDir = tracks[bestidx]->Dir().Y();
514  fzDir = tracks[bestidx]->Dir().Z();
515 
516 
517 //"(abs(xStart)>195||abs(yStart)>195||zStart<10)&&abs(xEnd)<150&&abs(yEnd)<150&&zEnd>50&&zEnd<1200"
518  if( ( abs(fxStart > 195) || abs(fyStart > 195) || abs(fzStart < 10) ) &&
519  abs(fxEnd) < 150 && abs(fyEnd) < 150 && fzEnd >50 && fzEnd < 1275
520 /*&& fzDir > 0.9 */){
521  fPassRockMuonCut = true;
522  }
523 
524 
525  unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();
526  // loop over trajectory points starting from the end
527  for(int itrj = ntraj-1; itrj >= 0; --itrj){
528  TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
529 
530  } // for itrj
531 
532  unsigned int minPlane = tracks[bestidx]->MinPlane();
533  unsigned int maxPlane = tracks[bestidx]->MaxPlane();
534  int maxXPlane = -1;
535  int maxYPlane = -1;
536  double hitCosx[tracks[bestidx]->NXCell()];
537  double hitCosy[tracks[bestidx]->NYCell()];
538 
539  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NXCell(); ++ihit){
540  if(tracks[bestidx]->XCell(ihit)->Plane() > maxXPlane){ maxXPlane = tracks[bestidx]->XCell(ihit)->Plane(); }
541  // get direction with this x hit
542  double xyzplane[3];
543  geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
544  double dz = geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->HalfD();
545 
546  double minz = xyzplane[2]-dz;
547  double maxz = xyzplane[2]+dz;
548  unsigned int currtrj = 0;
549  // find the matching trajectory point for this plane;
550  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
551  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
552  } // for itaraj
553  TVector3 direction = tracks[bestidx]->Dir();
554  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
555  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
556  }
557  direction = direction.Unit();
558  hitCosx[ihit] = sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
559  } // for ihit XCell
560 
561 
562 
563  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NYCell(); ++ihit){
564  if(tracks[bestidx]->YCell(ihit)->Plane() > maxYPlane){ maxYPlane = tracks[bestidx]->YCell(ihit)->Plane(); }
565  // get direction with this y hit
566  double xyzplane[3];
567  geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
568  double dz = geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->HalfD();
569  // double dz = HalfD;
570  double minz = xyzplane[2]-dz;
571  double maxz = xyzplane[2]+dz;
572  unsigned int currtrj = 0;
573  // find the matching trajectory point for this plane;
574  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
575  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
576  }
577  TVector3 direction = tracks[bestidx]->Dir();
578  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
579  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
580  }
581  direction = direction.Unit();
582  hitCosy[ihit] = sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
583  } // for ihit YCell
584 
585  ///
586 
587  // calculate dedx info for this track
588  std::vector<double> dedxs;
589  std::vector<double> des;
590  std::vector<double> pecorrs;
591  std::vector<double> mips;
592  std::vector<double> dxs;
593  std::vector<double> deads;
594  std::vector<double> xdedxs;
595  std::vector<unsigned int> dedxpls;
596 
597  std::vector<double> pes;
598  std::vector<double> adcs;
599  std::vector<double> xhits;
600  std::vector<double> yhits;
601  std::vector<double> zhits;
602  std::vector<double> thits;
603  std::vector<int> views;
604  std::vector<int> cells;
605  std::vector<int> planes;
606 
607  std::vector<double> mcdes;
608  std::vector<double> mcdebs;
609  std::vector<double> mcdxs;
610 
611  std::vector< std::vector<int> > mcells;
612  std::vector< std::vector<float> > madcs;
613  std::vector< std::vector<float> > mpes;
614  std::vector< std::vector<float> > mpecorrs;
615 
616  for(unsigned int iPlane = minPlane;iPlane<=maxPlane;++iPlane){
617  double planeGeV = 0;
618  double planePECorr = 0;
619  double planeMIP = 0;
620 
621  int planeView = -1;
622  int plane = -1;
623  int cell = -1;
624  double planeX = 0;
625  double planeY = 0;
626  double planeZ = 0;
627  double planeT = 0;
628  double planeADC = 0;
629  double planePE = 0;
630 
631  std::vector<int> mCell;
632  std::vector<float> mPE;
633  std::vector<float> mADC;
634  std::vector<float> mPECorr;
635 
636  double mcPlaneE = 0;
637  double mcPlaneEBirks = 0;
638  double mcPathLength = 0;
639  double totcosinview = 0;
640  double avgCosOtherView = 1.0;
641  bool useOtherView = true;
642  double nhits = 0;
643  double xyzplane[3];
644  geom->Plane(iPlane)->Cell(0)->GetCenter(xyzplane);
645  // double dz = geom->Plane(iPlane)->Cell(0)->HalfD();
646  // double minz = xyzplane[2]-dz;
647  // double maxz = xyzplane[2]+dz;
648  // int matchpt = 0; // Not used for now
649  int minCell = 10000;
650  int maxCell = 0;
651  double nhitsother = 0;
652  int view = 0;
653  int fbIndex = -1;
654  fFiberBrightness = -1;
655  // find the matching trajectory point for this plane;
656  /* Not used for now
657  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
658  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){matchpt = itraj;}
659  }
660  */
661 
662  if(geom->Plane(iPlane)->View() == geo::kX){
663  planeView = 0;
664  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
665  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane && tracks[bestidx]->XCell(iHit)->View() != geo::kY){
666  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
667  planeADC+=chit->ADC();
668  planePE +=chit->PE();
669  cell = chit->Cell();
670  plane = chit->Plane();
671 
672  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->XCell(iHit));
673  if(minPlane == maxPlane){
674  rhit = cal->MakeRecoHit(*(tracks[bestidx]->XCell(iHit)), 0.25);
675  }
676  if(rhit.IsCalibrated()){
677  planeGeV+=rhit.GeV();
678  planePECorr+=rhit.PECorr();
679  planeMIP+=rhit.MIP();
680  planeX += rhit.X();
681  planeY += rhit.Y();
682  planeZ += rhit.Z();
683  planeT += rhit.T();
684  totcosinview+=hitCosx[iHit];
685  mCell.push_back( chit->Plane() );
686  mPE.push_back( chit->PE() );
687  mADC.push_back( chit->ADC() );
688  mPECorr.push_back(rhit.PECorr() );
689 
690  nhits+=1.0;
691  if(tracks[bestidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->XCell(iHit)->Cell();}
692  if(tracks[bestidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->XCell(iHit)->Cell();}
693 
694  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->XCell(iHit)->Cell() );
695 
696  if(bt->HaveTruthInfo()){
697  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
698  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
699  for(int j = 0; j < (int)flsHits.size(); j++){
700  mcPlaneE += flsHits[j].GetEdep();
701  mcPlaneEBirks += flsHits[j].GetEdepBirks();
702  mcPathLength += flsHits[j].GetTotalPathLength();
703 
704  }//End loop over flshits in hit
705  } // if BackTracker
706 
707  }
708  }
709  } // for iHit
710 
711  // try to get the avgcos from the other view
712  double nhitsbefore = 0;
713  double nhitsafter = 0;
714  double totCosBefore = 0;
715  double totCosAfter = 0;
716  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
717  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane-1){
718  ++nhitsbefore;
719  totCosBefore+=hitCosy[iHit];
720  }
721  else if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane+1){
722  ++nhitsafter;
723  totCosAfter+=hitCosy[iHit];
724  }
725  }
726  if(nhitsbefore > 0 && nhitsafter > 0){
727  avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
728  nhitsother = (nhitsbefore+nhitsafter)/2;
729  }
730  else if(nhitsbefore > 0){
731  avgCosOtherView = totCosBefore/nhitsbefore;
732  nhitsother = nhitsbefore;
733  }
734  else if(nhitsafter > 0){
735  avgCosOtherView = totCosAfter/nhitsafter;
736  nhitsother = nhitsafter;
737  }
738  else{ useOtherView = false; }
739 
740  } else if(geom->Plane(iPlane)->View() == geo::kY){
741  planeView = 1;
742  view = 1;
743  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
744  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane && tracks[bestidx]->YCell(iHit)->View() != geo::kX
745  && tracks[bestidx]->YCell(iHit)->View() != geo::kXorY){
746  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
747  planeADC+=chit->ADC();
748  planePE +=chit->PE();
749  cell = chit->Cell();
750  plane = chit->Plane();
751 
752  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->YCell(iHit));
753  if(minPlane == maxPlane){
754  rhit = cal->MakeRecoHit(*(tracks[bestidx]->YCell(iHit)), 0.25);
755  }
756  if(rhit.IsCalibrated()){
757  planeGeV+=rhit.GeV();
758  planePECorr+=rhit.PECorr();
759  planeMIP+=rhit.MIP();
760  planeX += rhit.X();
761  planeY += rhit.Y();
762  planeZ += rhit.Z();
763  planeT += rhit.T();
764  totcosinview+=hitCosy[iHit];
765 
766  mCell.push_back( chit->Plane() );
767  mPE.push_back( chit->PE() );
768  mADC.push_back( chit->ADC() );
769  mPECorr.push_back(rhit.PECorr() );
770 
771  nhits+=1.0;
772  if(tracks[bestidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->YCell(iHit)->Cell();}
773  if(tracks[bestidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->YCell(iHit)->Cell();}
774 
775  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->YCell(iHit)->Cell() );
776 
777  if(bt->HaveTruthInfo()){
778  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
779  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
780  for(int j = 0; j < (int)flsHits.size(); j++){
781  mcPlaneE += flsHits[j].GetEdep();
782  mcPlaneEBirks += flsHits[j].GetEdepBirks();
783  mcPathLength += flsHits[j].GetTotalPathLength();
784 
785  }//End loop over flshits in hit
786  } // if BackTracker
787 
788  } // if isCalibrated
789  } //
790  } // for iHit
791 
792  // try to get the avgcos from the other view
793  double nhitsbefore = 0;
794  double nhitsafter = 0;
795  double totCosBefore = 0;
796  double totCosAfter = 0;
797  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
798  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane-1){
799  ++nhitsbefore;
800  totCosBefore+=hitCosx[iHit];
801  }
802  else if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane+1){
803  ++nhitsafter;
804  totCosAfter+=hitCosx[iHit];
805  }
806  }
807  if(nhitsbefore > 0 && nhitsafter > 0){
808  avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
809  nhitsother = (nhitsbefore+nhitsafter)/2;
810  }
811  else if(nhitsbefore > 0){
812  avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
813  nhitsother = nhitsbefore;
814  }
815  else if(nhitsafter > 0){
816  avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
817  nhitsother = nhitsbefore;
818  }
819  else{ useOtherView = false; }
820 
821  } // if kY
822  fnCellsPerPlane = nhits;
823  if(nhits >0){
824 
825  if(nhits == 1) {
826  fFiberBrightness = fbIndex;
827  }
828 
829  double avgCos = TMath::Abs(totcosinview/nhits);
830  double xyzp[3];
831  geom->Plane(iPlane)->Cell(minCell)->GetCenter(xyzp);
832  double xyzp1[3];
833  geom->Plane(iPlane)->Cell(maxCell)->GetCenter(xyzp1);
834  double planewidth = 2*geom->Plane(iPlane)->Cell(1)->HalfD();
835  // double planewidth = 2*HalfD;
836  // calculate the amount of dead material to subtract off
837  // first add up live material
838  double liveLength = 0;
839  for(int i = minCell; i<maxCell+1;++i){
840  // get the cell half width
841  const geo::PlaneGeo* p = geom->Plane(iPlane);
842  const geo::CellGeo* c = p->Cell(i);
843  if(i == minCell || i == maxCell){ liveLength+=c->HalfW(); }
844  else{ liveLength+=(2.0*c->HalfW() ); }
845  }
846  // dead material is total length - live length only if min and max cell are not the same
847  double deadLength = 0;
848  if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
849  double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
850  double planePathLength = TMath::Sqrt(planewidth*planewidth + planeheight*planeheight);
851 
852  // in the case of completely vertical track be a bit careful
853  if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }
854  if(useOtherView){
855  double deltaz = planewidth;
856  double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
857  if(avgCos == 0 || avgCos!=avgCos){
858  deltav = liveLength;
859  // if the track looks completely vertical in both views,
860  // than no part of the track length should be in z direction
861  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
862  }
863 
864  // in opposite view dead length
865  double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
866  // in the case that the other view looks completely vertical scale up
867  // the length between cells in this view by the number of hits in the other view
868  // this case should really never happen
869  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
870  deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
871  }
872  planePathLength = TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
873  }
874  // also the case of completely vertical track
875  if(minPlane == maxPlane){ planePathLength = liveLength; }
876 
877  assert(planeGeV/(planePathLength)>=0);
878  dedxs.push_back(planeGeV/(planePathLength));
879  des.push_back(planeGeV);
880  pecorrs.push_back(planePECorr);
881  mips.push_back(planeMIP);
882  dxs.push_back(planePathLength);
883  deads.push_back(deadLength);
884  xhits.push_back(planeX/nhits);
885  yhits.push_back(planeY/nhits);
886  zhits.push_back(planeZ/nhits);
887  thits.push_back(planeT/nhits);
888  views.push_back(planeView);
889  planes.push_back(plane);
890  cells.push_back(cell);
891  pes.push_back(planePE);
892  adcs.push_back(planeADC);
893 
894  mcdes.push_back(mcPlaneE);
895  mcdebs.push_back(mcPlaneEBirks);
896  mcdxs.push_back(mcPathLength);
897 
898  mcells.push_back(mCell);
899  madcs.push_back(mADC);
900  mpes.push_back(mPE);
901  mpecorrs.push_back(mPECorr);
902 
903  mcdes.push_back(mcPlaneE);
904  mcdebs.push_back(mcPlaneEBirks);
905  mcdxs.push_back(mcPathLength);
906 
907  } //nhits > 0
908  } // iPlane
909 
910  // ..... .......... .......... .......... .......... .......... .....
911 
912  art::Ptr<remid::ReMId> pid = foids.at(bestidx);
913 
914  //Looking at numu energy things
915  numue::NumuE numuE;
916  if (!fndmnyNumuEnergy.isValid()){ continue; }
917  std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
918  if (energyVec.empty()){ continue; }
919  numuE = *energyVec[0];
920 
921  //Require decent energy object; not too much hadronic energy in muon catcher
922  if (numuE.E() == -5){ continue; }
923  if (!(numuE.NDHadCalCat() + numuE.NDHadCalTran() < 0.03)){ continue; }
924  //Not much hadronic energy contamination on muon track
925  if (numuE.HadTrkE() < 0.025){ fPassHadEOnMuonTrack = true; }
926 
927  //Looking at qepid things //ecm
928 /*
929  qeef::QePId qepid;
930  if (!fndmnyQepid.isValid()){ continue; }
931  std::vector<art::Ptr<qeef::QePId> > qVec = fndmnyQepid.at(iSlice);
932  if (qVec.empty()){ continue; }
933  qepid = *qVec[0];
934 
935  //Require ratio of average dE/dx track is good
936  if (qepid.Dedx() > 1.75){ fPassDedxRatio = true; }
937  //Require Ztest of two methods of calculating QE neutrino energy is good
938  if (qepid.EdiffZ() > -1.3){ fPassEDiffZTest = true; }
939 */
940 
941 // ECM need to calculate qepid variables
942 
943 //Require ratio of average dE/dx of two tracks is good
944 
945  if(dedxs.size()>0){
946  double dedx_avg_muon = std::accumulate(dedxs.begin(),dedxs.end(),0.) / dedxs.size();
947 // double dedx_avg_proton = std::accumulate(dedxp.begin(),dedxp.end(),0.) / dedxp.size();
948 // std::cerr << std::accumulate(dedxs.begin(),dedxs.end(),0.) << " size " << dedxs.size() << std::endl;
949 // std::cerr << std::accumulate(dedxp.begin(),dedxp.end(),0.) << " size " << dedxp.size() << std::endl;
950 // std::cerr << "muon" << dedx_avg_muon << " proton " << dedx_avg_proton << " ratio " << dedx_avg_muon /dedx_avg_proton << std::endl;
951 // if(dedx_avg_proton / dedx_avg_muon > 1.75) fPassDedxRatio = true;
952 
953 //ECM THIS MUST BE A MISTAKE
954  if(dedx_avg_muon > 1.75) fPassDedxRatio = true;
955  }
956 
957  //Require Ztest of two methods of calculating QE neutrino energy is good
958  // std::cerr << numuE.AngleQEE() << " " << numuE.TrkQEE() << " " << numuE.AngleQEError() << std::endl;
959  //std::cerr << (numuE.AngleQEE()-numuE.TrkQEE())/numuE.AngleQEError() << std::endl;
960  double ediffz = (numuE.AngleQEE()-numuE.TrkQEE())/numuE.AngleQEError();
961  if (ediffz > -1.3){ fPassEDiffZTest = true; }
962 
963 
964  //Slice energy mostly on one tracks
965  double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy());
966  if (calEDiff < 0.1){ fPassOffTrackEnergy = true; }
967 
968  //Slice hits mostly on two tracks
969  int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell());
970  //if (!(hitDiff < 4)){ continue; }
971  if (hitDiff < 4){ fPassOffTrackHits = true; }
972 
973  //Check if true signal for sim events
974  fTrueNumuCCQE = false;
975  fTrueNumuCCWithProton = false;
977 
978  fpdg = 0;
979 
980  if (bt->HaveTruthInfo()){
981  std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),allhits);
982  if (!funcReturn.empty()){
983  // const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();
984 
985  const std::vector<const sim::Particle*> test_mu = bt->HitsToParticle(tracks[bestidx]->AllCells());
986  if (!test_mu.empty()) {
987  fpdg = test_mu[0]->PdgCode();
988  fTrueE = test_mu[0]->E();
989  // std::cout << "Mu " << fTrueE << " " << fpdg << std::endl;
990  }
991  }
992  }
993 
994 // if(ContainedEvent(slicelist[iSlice],tracks[bestidx])){
995  if(fPassRockMuonCut )
996 // if(1 )
997  { // Containment cut removed for now. Always pass //ecm try to undo
999  fSlicesPassingSelectionCuts->Fill(0.5);
1000 /*
1001  std::cout << evt.run() << " "
1002  << evt.subRun() << " "
1003  << evt.event() << " "
1004  << iSlice
1005  << std::endl;
1006 */
1007  if (fTrueNumuCCWithProton){
1009  }
1010  }
1011  // get all the dedx variables filled
1012  fnDedx = pid->NDedx()-1;
1013  fnDedxUsed = pid->NMeasurementPlanes();
1014  // loop over all of the dE/dx measurements
1015  for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
1016  fDedx = pid->Dedx(idedx);
1017  fVert = pid->DedxVertex(idedx);
1018  fUsed = pid->DedxUsed(idedx);
1019 
1020  fDe = des[idedx];
1021  fPECorr = pecorrs[idedx];
1022  fMIP = mips[idedx];
1023  fDx = dxs[idedx];
1024  fxDead = deads[idedx];
1025  fmyDedx = fDe/fDx;
1026 
1027  fxHit = xhits[idedx];
1028  fyHit = yhits[idedx];
1029  fzHit = zhits[idedx];
1030  ftHit = thits[idedx];
1031  fView = views[idedx];
1032  fPlane = planes[idedx];
1033  fCell = cells[idedx];
1034  fADC = adcs[idedx];
1035  fPE = pes[idedx];
1036 
1037  fMultiCell = mcells[idedx];
1038  fMultiADC = madcs[idedx];
1039  fMultiPE = mpes[idedx];
1040  fMultiPECorr = mpecorrs[idedx];
1041 
1042  fmcDe = mcdes[idedx];
1043  fmcDeBirks = mcdebs[idedx];
1044  fmcDx = mcdxs[idedx];
1045  // get the distance from the end of the track
1046  fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[bestidx]->TotalLength();
1047  // fill the dedx sample tree
1048  DedxSample->Fill();
1049  }//Loop over DedxSample
1050  }//If passes contaiment
1051  } // end loop over slices
1052 
1053  }
1054 
1055 
1056  //.......................................................................
1057 
1059  {
1060  fPOT->Fill(.5, fTotalPOT);
1061  }
1062 
1063 
1064  //.......................................................................
1066  {
1067 
1068  bool contained = false;
1069 
1070  // check that this is a contained event
1071 
1072  // get some useful information here
1073  // 1) how close is the slice to the edge of the detector
1074  unsigned int slcncellsfromedge = 999999999;
1075  for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
1076  const art::Ptr<rb::CellHit>& chit = slice->Cell(ihit);
1077  const int plane = chit->Plane();
1078  unsigned int cellsfromedge = std::min((unsigned int)chit->Cell(), fgeom->Plane(plane)->Ncells() - 1 - chit->Cell());
1079  if(cellsfromedge < slcncellsfromedge){
1080  slcncellsfromedge = cellsfromedge;
1081  }
1082  }
1083  // 2) how close is the track to the edge of the detector
1084  // in planes
1085  //SL - NOPE, need to ask slice, not track
1086  unsigned int firstpl = slice->MinPlane();
1087  unsigned int lastpl = slice->MaxPlane();
1088  // in projected cells, just define directions here actually get projected cell below
1089  TVector3 dirbwk = -trk->Dir();
1090  TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();
1091 
1092  // 3) hits in slice
1093  int sliceNHit = slice->NCell();
1094 
1095  // 4) slice mean time
1096  double slcTime = slice->MeanTNS();
1097 
1098  // 5) hits on track
1099  int trackNHit = trk->NCell();
1100 
1101  //std::cout<<"In contained event loop!"<<std::endl;
1102  // cut the slice number of cells from edge
1103  if((sliceNHit > 20)&&(slcncellsfromedge > 1)&&(slcTime>50000)&&(slcTime<450000)&&(trackNHit > 5)){
1104  if(firstpl > 1 && lastpl < 212){
1105  //std::cout<<"In contained event loop! 1"<<std::endl;
1106  int fwdcell = flivegeom->FullNDProjectedCells(trk->Stop(), dirfwd);
1107  int bwkcell = flivegeom->FullNDProjectedCells(trk->Start(),dirbwk);
1108  // also make sure the track doesn't start too close to the end of the active region
1109  if(fwdcell > 4 && bwkcell > 8 && trk->Start().Z() < 1150){
1110  // if the track ends in the muon catcher, check that it doesn't clip the corner
1111  //std::cout<<"In contained event loop! 2"<<std::endl;
1112  if(trk->Stop().Z() < 1275){
1113  //std::cout<<"In contained event loop! 2a"<<std::endl;
1114  contained = true;
1115  }
1116  else{
1117  //std::cout<<"In contained event loop! 2b"<<std::endl;
1118  /*
1119  // find the trajectory point closest in z to the transition plane
1120  double mindist = 99999999;
1121  unsigned int mintrajid = 0;
1122  for(unsigned int itraj = 0; itraj < trk->NTrajectoryPoints()-1; ++itraj){
1123  double dist = 1275 - trk->TrajectoryPoint(itraj).Z();
1124  if(dist > 0 && dist < mindist){
1125  mindist = dist;
1126  mintrajid = itraj;
1127  }
1128  }
1129  // get the direction at this trajectory point
1130  TVector3 trandir = trk->Dir();
1131  if(mintrajid < trk->NTrajectoryPoints() -2){
1132  TVector3 trandir = (trk->TrajectoryPoint(mintrajid+1) - trk->TrajectoryPoint(mintrajid)).Unit();
1133  }
1134  */
1135  //std::cout<<"In contained event loop! 3"<<std::endl;
1136  //double ypostran = flivegeom->YPositionAtTransition(trk->TrajectoryPoint(mintrajid),trandir);
1137  double ypostranKirk = flivegeom->YPositionAtTransition(trk->Start(),trk->Dir());
1138  if( ypostranKirk < 55){
1139  //std::cout<<"In contained event loop! 4"<<std::endl;
1140  contained = true;
1141  }
1142 
1143  } // end else in muon catcher
1144  } // end if track ends projections contained
1145  } // end if track end planes contained
1146  // } // end if expected number of diblocks
1147  } // end if slice contained
1148 
1149 
1150  return contained;
1151 
1152  }
1153 
1154 
1155 
1156 
1157 }
1158 
1159 
1160 namespace remid{
1161 
1163 
1164 }
float T() const
Definition: RecoHit.h:39
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
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
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 abs(TH1 *hist)
void beginRun(const art::Run &run)
art::ServiceHandle< geo::Geometry > fgeom
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:27
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:21
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
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
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:435
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
double posy
mm
Definition: SpillData.h:35
ReMIdDedxRock(fhicl::ParameterSet const &pset)
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:302
float NDHadCalTran() const
Definition: NumuE.cxx:334
float AngleQEE() const
Definition: NumuE.cxx:259
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
float GeV() const
Definition: RecoHit.cxx:69
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
float AngleQEError() const
Definition: NumuE.cxx:264
double DedxLocation(unsigned int i) const
Definition: ReMId.cxx:120
art::ServiceHandle< geo::LiveGeometry > flivegeom
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:137
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
float TrkQEE() const
Definition: NumuE.cxx:244
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:190
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.
enum BeamMode string