ReMIdDedxStudies_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////
2 // \file ReMIdDedxStudies_module.cc
3 // \brief A module for studying dE/dx along tracks for QE like events.
4 // Based on ReMIdDedx_module.cc with extended info in the output Tree.
5 // \version $Id: ReMIdDedxStudies_module.cc,v 2017/06/20
6 // \author Stefano Germani - s.germani@ucl.ac.uk
7 ////////////////////////////////////////////////////////////////////
8 #include <string>
9 
10 
11 // ROOT includes
12 #include "TVector3.h"
13 #include "TTree.h"
14 #include "TH1D.h"
15 #include "TH2D.h"
16 #include "TMath.h"
17 
18 // Framework includes
19 #include "fhiclcpp/ParameterSet.h"
21 #include "Utilities/AssociationUtil.h"
26 
27 // NOvA includes
28 //....................................
29 #include "Calibrator/Calibrator.h"
31 //....................................
32 #include "Geometry/Geometry.h"
33 #include "Geometry/LiveGeometry.h"
35 #include "ReMId/ReMId.h"
36 #include "NumuEnergy/NumuE.h"
37 #include "QEEventFinder/QePId.h"
38 #include "GeometryObjects/Geo.h"
41 #include "MCCheater/BackTracker.h"
42 #include "RecoBase/Track.h"
43 //................................
44 #include "RecoBase/Cluster.h"
45 #include "RecoBase/CellHit.h"
46 #include "RecoBase/RecoHit.h"
47 
48 #include "MCCheater/SimHit.h"
49 
50 //.................................
51 
52 #include "Simulation/Simulation.h"
54 #include "SummaryData/POTSum.h"
55 #include "SummaryData/SpillData.h"
56 #include "MCCheater/BackTracker.h"
57 
58 //....
59 
60 namespace remid {
61 
62  // A module to analyze ReMId objects for pid training
64  public:
65  explicit ReMIdDedxStudies(fhicl::ParameterSet const& pset);
67 
68  void analyze(art::Event const& evt);
69 
70  void beginJob();
71  void beginRun(const art::Run& run);
72  void endJob();
73  //virtual void beginSubRun(art::SubRun& sr);
74 
75  private:
77  std::string fSliceLabel; ///<Where to find reconstructed slices
78  std::string fTrackLabel; ///<Where to find recondtructed tracks
79  std::string fPidLabel; ///<Where to find the pids
85 
86  TH1D* fPOT; ///<Histogram of POT passing spill cuts
93 
94  double fTotalPOT;
95  bool fND; ///< Is detector ND?
96 
97  // variables used in filling the kNN training trees
98  float fTrackLength;
100  //float fDedxSep;
101  //float fScatSep;
102  //float fnuE;
103  //float fvisE;
104  //float frecoE;
105  //float fMeasFrac;
106 
107  TTree* DedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
108  TTree* ProtonDedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
109 
110  // variables used in filling the sample histogram trees
111  float fDedx;
112  int fnDedx;
114  float fDistDedx;
115  float fProtonDedx;
119  //float fCosScat;
120  //float fDistLastScat;
121  //float fScatVar;
122  //int fnScat;
123  //float fDistScat;
124  //int fnPlane;
125  int fpdg;
126  //int fintType;
127  bool fVert;
128  bool fUsed;
129  float fzStart;
130  float fzEnd;
134  float fProtonzEnd;
147 
148  // variables added to trees for detailed studies
150  float fnewDedx;
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  float fTrueE;
177 
178 
181  float fProtonDe;
185  float fProtonADC;
186  float fProtonPE;
188  float fProtonMIP;
189  float fProtonDx;
192  float fProtonxEnd;
194  float fProtonyEnd;
195  float fProtonxHit;
196  float fProtonyHit;
197  float fProtontHit;
198  float fProtonzHit;
199  float fProtonxDir;
200  float fProtonyDir;
201  float fProtonzDir;
202  float fProtonmcDe;
204  float fProtonmcDx;
205 
209  // ......................
210 
215 
217 
218  //novadaq::cnv::DetId fDetID; ///< Detector ID in nova daq convention namespace typedef
219 
220 
221  //unsigned int fdibfirst = 0;
222  //unsigned int fdiblast = 0;
223 
225 
226  };
227 
228 }
229 
230 
231 namespace remid{
232 
233 
234  //.......................................................................
235 
237  EDAnalyzer(pset)
238  {
239  fHitModuleLabel = (pset.get< std::string >("HitModuleLabel"));
240  fSliceLabel = (pset.get< std::string >("SliceLabel"));
241  fTrackLabel = (pset.get< std::string >("TrackLabel"));
242  fPidLabel = (pset.get< std::string >("PidLabel" ));
243  fGeneratorLabel = pset.get< std::string >("GeneratorInput");
244  fNumiBeamLabel = pset.get< std::string >("NuMIBeamInput");
245  fSpillQualLabel = pset.get< std::string >("SpillQualInput");
246  fNumuEnergyLabel = pset.get< std::string >("NumuEnergyInput");
247  fQepidLabel = pset.get< std::string >("QepidInput");
248  }
249 
250  //.......................................................................
251 
253  {
254 
255  }
256 
257  //.......................................................................
258 
260  {
262 
263  // Muon Tree
264  DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
265  DedxSample->Branch("dedx", &fDedx, "dedx/F");
266  DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
267  DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
268  DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
269  DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
270  DedxSample->Branch("vert", &fVert, "vert/O");
271  DedxSample->Branch("used", &fUsed, "used/O");
272  DedxSample->Branch("zStart", &fzStart, "zStart/F");
273  DedxSample->Branch("zEnd", &fzEnd, "zEnd/F");
274  DedxSample->Branch("passOffTrackHits", &fPassOffTrackHits, "passOffTrackHits/O");
275  DedxSample->Branch("passOffTrackEnergy", &fPassOffTrackEnergy, "passOffTrackEnergy/O");
276  DedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");
277  DedxSample->Branch("passVertX", &fPassVertX, "passVertX/O");
278  DedxSample->Branch("passVertY", &fPassVertY, "passVertY/O");
279  DedxSample->Branch("passVertZ", &fPassVertZ, "passVertZ/O");
280  DedxSample->Branch("passProtonRemid", &fPassProtonRemid, "passProtonRemid/O");
281  DedxSample->Branch("passDedxRatio", &fPassDedxRatio, "passDedxRatio/O");
282  DedxSample->Branch("passEDiffZTest", &fPassEDiffZTest, "passEDiffZTestX/O");
283  DedxSample->Branch("trueNumuCCQE", &fTrueNumuCCQE, "trueNumuCCQE/O");
284 
285  // Added Branches wrt standard ReMIdDedx Module
286  DedxSample->Branch("ncellsPerPlane", &fnCellsPerPlane, "ncellsPerPlane/I");
287  DedxSample->Branch("newdedx", &fnewDedx, "newdedx/F");
288  DedxSample->Branch("de", &fDe, "de/F");
289  DedxSample->Branch("view", &fView, "view/I");
290  DedxSample->Branch("cell", &fCell, "cell/I");
291  DedxSample->Branch("plane", &fPlane, "plane/I");
292  DedxSample->Branch("adc", &fADC, "adc/F");
293  DedxSample->Branch("pe", &fPE, "pe/F");
294  DedxSample->Branch("pecorr", &fPECorr, "pecorr/F");
295  DedxSample->Branch("mip", &fMIP, "mip/F");
296  DedxSample->Branch("dx", &fDx, "dx/F");
297  DedxSample->Branch("deadx", &fxDead, "deadx/F");
298  DedxSample->Branch("xStart", &fxStart, "xStart/F");
299  DedxSample->Branch("xEnd", &fxEnd, "xEnd/F");
300  DedxSample->Branch("yStart", &fyStart, "yStart/F");
301  DedxSample->Branch("yEnd", &fyEnd, "yEnd/F");
302  DedxSample->Branch("xHit", &fxHit, "xHit/F");
303  DedxSample->Branch("yHit", &fyHit, "yHit/F");
304  DedxSample->Branch("zHit", &fzHit, "zHit/F");
305  DedxSample->Branch("tHit", &ftHit, "tHit/F");
306  DedxSample->Branch("xDir", &fxDir, "xDir/F");
307  DedxSample->Branch("yDir", &fyDir, "yDir/F");
308  DedxSample->Branch("zDir", &fzDir, "zDir/F");
309  DedxSample->Branch("mcde", &fmcDe, "mcde/F");
310  DedxSample->Branch("mcdeBirks",&fmcDeBirks, "mcdeBirks/F");
311  DedxSample->Branch("mcdx", &fmcDx, "mcdx/F");
312 
313  DedxSample->Branch("pdg", &fpdg, "pdg/I");
314  DedxSample->Branch("trueE", &fTrueE, "trueE/F");
315  DedxSample->Branch("fiberBrightness", &fFiberBrightness, "fiberBrightness/I");
316 
317 
318  // Proton Tree
319  ProtonDedxSample = tfs->make<TTree>("fProtonDedxSample","Variables for Creating Proton Dedx Sample Histograms");
320  ProtonDedxSample->Branch("Protondedx", &fProtonDedx, "Protondedx/F");
321  ProtonDedxSample->Branch("ProtonnDedx", &fProtonnDedx, "ProtonnDedx/I");
322  ProtonDedxSample->Branch("ProtonnDedxUsed", &fProtonnDedxUsed, "ProtonnDedxUsed/I");
323  ProtonDedxSample->Branch("ProtondistFromEnd", &fProtonDistDedx, "ProtondistFromEnd/F");
324  ProtonDedxSample->Branch("ProtontrackLength", &fProtonTrackLength, "ProtontrackLength/F");
325  ProtonDedxSample->Branch("Protonvert", &fProtonVert, "Protonvert/O");
326  ProtonDedxSample->Branch("Protonused", &fProtonUsed, "Protonused/O");
327  ProtonDedxSample->Branch("ProtonzStart",&fProtonzStart,"ProtonzStart/F");
328  ProtonDedxSample->Branch("ProtonzEnd",&fProtonzEnd,"ProtonzEnd/F");
329  ProtonDedxSample->Branch("passOffTrackHits",&fPassOffTrackHits,"passOffTrackHits/O");
330  ProtonDedxSample->Branch("passOffTrackEnergy",&fPassOffTrackEnergy,"passOffTrackEnergy/O");
331  ProtonDedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");
332  ProtonDedxSample->Branch("passVertX",&fPassVertX,"passVertX/O");
333  ProtonDedxSample->Branch("passVertY",&fPassVertY,"passVertY/O");
334  ProtonDedxSample->Branch("passVertZ",&fPassVertZ,"passVertZ/O");
335  ProtonDedxSample->Branch("passProtonRemid",&fPassProtonRemid,"passProtonRemid/O");
336  ProtonDedxSample->Branch("passDedxRatio",&fPassDedxRatio,"passDedxRatio/O");
337  ProtonDedxSample->Branch("passEDiffZTest",&fPassEDiffZTest,"passEDiffZTestX/O");
338  ProtonDedxSample->Branch("trueNumuCCQE",&fTrueNumuCCQE,"trueNumuCCQE/O");
339 
340  // Added Branches wrt standard ReMIdDedx Module
341  ProtonDedxSample->Branch("ProtonncellsPerPlane", &fProtonnCellsPerPlane, "ProtonncellsPerPlane/I");
342  ProtonDedxSample->Branch("Protonde", &fProtonDe, "Protonde/F");
343  ProtonDedxSample->Branch("Protonview", &fProtonView, "Protonview/I");
344  ProtonDedxSample->Branch("Protoncell", &fProtonCell, "Protoncell/I");
345  ProtonDedxSample->Branch("Protonplane", &fProtonPlane, "Protonplane/I");
346  ProtonDedxSample->Branch("Protonadc", &fProtonADC, "Protonadc/F");
347  ProtonDedxSample->Branch("Protonpe", &fProtonPE, "Protonpe/F");
348  ProtonDedxSample->Branch("Protonpecorr", &fProtonPECorr, "Protonpecorr/F");
349  ProtonDedxSample->Branch("Protonmip", &fProtonMIP, "Protonmip/F");
350  ProtonDedxSample->Branch("Protondx", &fProtonDx, "Protondx/F");
351  ProtonDedxSample->Branch("Protondeadx", &fProtonxDead, "Protondeadx/F");
352  ProtonDedxSample->Branch("Protonnewdedx", &fProtonnewDedx, "Protonnewdedx/F");
353  ProtonDedxSample->Branch("ProtonxStart", &fProtonxStart, "ProtonxStart/F");
354  ProtonDedxSample->Branch("ProtonxEnd", &fProtonxEnd, "ProtonxEnd/F");
355  ProtonDedxSample->Branch("ProtonyStart", &fProtonyStart, "ProtonyStart/F");
356  ProtonDedxSample->Branch("ProtonyEnd", &fProtonyEnd, "ProtonyEnd/F");
357  ProtonDedxSample->Branch("ProtonxHit", &fProtonxHit, "ProtonxHit/F");
358  ProtonDedxSample->Branch("ProtonyHit", &fProtonyHit, "ProtonyHit/F");
359  ProtonDedxSample->Branch("ProtonzHit", &fProtonzHit, "ProtonzHit/F");
360  ProtonDedxSample->Branch("ProtontHit", &fProtontHit, "ProtontHit/F");
361  ProtonDedxSample->Branch("ProtonxDir", &fProtonxDir, "ProtonxDir/F");
362  ProtonDedxSample->Branch("ProtonyDir", &fProtonyDir, "ProtonyDir/F");
363  ProtonDedxSample->Branch("ProtonzDir", &fProtonzDir, "ProtonzDir/F");
364  ProtonDedxSample->Branch("Protonmcde", &fProtonmcDe, "Protonmcde/F");
365  ProtonDedxSample->Branch("ProtonmcdeBirks",&fProtonmcDeBirks, "ProtonmcdeBirks/F");
366  ProtonDedxSample->Branch("Protonmcdx", &fProtonmcDx, "Protonmcdx/F");
367 
368  ProtonDedxSample->Branch("ProtonTrueE", &fProtonTrueE, "ProtonTrueE/F");
369  ProtonDedxSample->Branch("Protonpdg", &fProtonpdg, "Protonpdg/I");
370  ProtonDedxSample->Branch("ProtonfiberBrightness", &fProtonFiberBrightness, "ProtonfiberBrightness/I");
371  ProtonDedxSample->Branch("trueNumuCCWithProton", &fTrueNumuCCWithProton, "trueNumuCCWithProton/O");
372  ProtonDedxSample->Branch("superTrueNumuCCWithProton", &fSuperTrueNumuCCWithProton, "superTrueNumuCCWithProton/O");
373 
374 
375  fPOT = tfs->make<TH1D>("POT", ";;Total POT",1, 0.0, 1.0);
376  fEventsPassingSpillCuts = tfs->make<TH1D>("EventsPassingSpillCuts", ";;Events Passing Spill Cuts",1, 0.0, 1.0);
377  fSlicesPassingSelectionCuts = tfs->make<TH1D>("SlicesPassingSelectionCuts", ";;Slices Passing Selection Cuts",1, 0.0, 1.0);
378  fTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("TrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
379  fSuperTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("SuperTrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
380  fTrueProtonEFromTrackLength = tfs->make<TH2D>("TrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
381  fSuperTrueProtonEFromTrackLength = tfs->make<TH2D>("SuperTrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
382 
383  fTotalPOT = 0.0;
384 
385 
386 
387  }
388 
389  //......................................................................
390 
392 
393  {
394  fND = false;
395 
396  novadaq::cnv::DetId detID = fgeom->DetId();
397  if (detID == novadaq::cnv::kNEARDET) fND = true;
398 
399  return;
400  }
401 
402 
403  //.......................................................................
404 
406  {
407  //Only written for ND
408 
409  if (!fND) return;
410 
411  // get a geometry handle **************************************
414  //***************************************************************
415 
416  //Does this event pass spill cuts?
418  if(bt->HaveTruthInfo()){
419  evt.getByLabel(fGeneratorLabel, spillPot);
420  }
421  else{
422  evt.getByLabel(fNumiBeamLabel, spillPot);
423  }
424 
425  if (spillPot.failedToGet()){
426  mf::LogError("Susie") << "Spill Data not found for real data event";
427  return;
428  }
429 
430  //SL -- only cut on beam quality for data
431  if (!bt->HaveTruthInfo()){
432  if (std::abs(spillPot->deltaspilltimensec) > 0.5e9) return;
433  if ((spillPot->spillpot)*1e12 < 2e12) return;
434  if ((spillPot->hornI < -202)||(spillPot->hornI > -198)) return;
435  if ((spillPot->posx < -2.00)||(spillPot->posx > 2.00)) return;
436  if ((spillPot->posy < -2.00)||(spillPot->posy > 2.00)) return;
437  if ((spillPot->widthx < 0.57)||(spillPot->widthx > 1.58)) return;
438  if ((spillPot->widthy < 0.57)||(spillPot->widthy > 1.58)) return;
439  }
440 
441  //Data quality spill cuts
442 
443  // Get the EventQuality information
445  evt.getByLabel(fSpillQualLabel, spillQual);
446 
447  if(spillQual.failedToGet()){
448  mf::LogError("CAFMaker") << "Spill EventQuality not found for real data event";
449  return;
450  }
451 
452  if (spillQual->fracdcm3hits > 0.45) return;
453  if (spillQual->nmissingdcms > 0) return;
454 
455 
456  //Filling spill pot after it passes the cuts
457  float spillpot = spillPot->spillpot;
458  if (!bt->HaveTruthInfo()) spillpot *= 1e12;
459  fTotalPOT += spillpot;
460 
461  fEventsPassingSpillCuts->Fill(0.5);
462 
463  //Getting event hits
465  evt.getByLabel(fHitModuleLabel, hithdl);
466 
468  for(size_t h = 0; h < hithdl->size(); ++h){
469  art::Ptr<rb::CellHit> hit(hithdl, h);
470  allhits.push_back(hit);
471  }
472 
473 
474 
475  // get the slices out of the event
477  evt.getByLabel(fSliceLabel,slicecol);
478 
479 
480  art::PtrVector<rb::Cluster> slicelist;
481 
482  for(unsigned int i = 0; i<slicecol->size();++i){
483  art::Ptr<rb::Cluster>slice(slicecol,i);
484  slicelist.push_back(slice);
485  }
486 
487  // get assosciations between slices and tracks
488  art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
489 
490  // get assosciations between slices and numu energy
491  art::FindManyP<numue::NumuE> fndmnyNumuEnergy(slicecol, evt, fNumuEnergyLabel);
492  // get assosciations between slices and qepid
493  art::FindManyP<qeef::QePId> fndmnyQepid(slicecol, evt, fQepidLabel);
494 
495  // loop over the slices
496  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
497 
498  // If slice is noise, skip it
499  if(slicelist[iSlice]->IsNoise()){ continue; }
500 
501  // get all of the tracks associated to this slice
502  std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
503  // if no tracks move on
504  if(tracks.size() == 0){ continue; }
505  //Require precisely two Kalman tracks
506  if(tracks.size() != 2){ continue; }
507 
508  //Require BOTH are 3D
509  bool all3D = true;
510 
511  // get the pid that belongs to this recotrack
512  art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);
513  if(!foids.isValid()) continue;
514  if(foids.at(0)->Value() >1 || foids.at(0)->Value()<0) continue;
515  if(foids.at(1)->Value() >1 || foids.at(1)->Value()<0) continue;
516 
517  int bestidx = -1;
518  double bestpid = -1;
519  // loop over all of the tracks and find the most muon like one;
520  for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
521  if(!tracks[itrk]->Is3D()){
522  all3D = false;
523  continue;
524  }
525  // get the pid that belongs to this track
526  art::Ptr<remid::ReMId> pid = foids.at(itrk);
527  if(pid->Value() > bestpid){
528  bestpid = pid->Value();
529  bestidx = itrk;
530  }
531  }//End of loop over tracks
532 
533  if (!all3D){ continue; }
534 
535  //Default all pass cuts to false
536  fPassOffTrackHits = false;
537  fPassOffTrackEnergy = false;
538  fPassHadEOnMuonTrack = false;
539  fPassVertX = false;
540  fPassVertY = false;
541  fPassVertZ = false;
542  fPassProtonRemid = false;
543  fPassDedxRatio = false;
544  fPassEDiffZTest = false;
545 
546 
547  //Slightly harsh ReMId cut on muon track
548  if(bestidx == -1 || bestpid < 0.8){ continue; }
549 
550  if((bestidx != 0) && (bestidx !=1)){
551  std::cout<<"BUMMER! A weird track id for muon that shouldn't happen!!"<<bestidx<<std::endl;
552  continue;
553  }
554 
555  int protonidx = -1;
556  if (bestidx == 0){
557  protonidx = 1;
558  }
559  if (bestidx == 1){
560  protonidx = 0;
561  }
562 
563  if(protonidx == -1){
564  std::cout<<"BUMMER! A weird track id for proton that shouldn't happen!!"<<std::endl;
565  continue;
566  }
567 
568  //Filling muon track length properties
569  fzStart = tracks[bestidx]->Start().Z();
570  fzEnd = tracks[bestidx]->Stop().Z();
571  fTrackLength = tracks[bestidx]->TotalLength();
572 
573  //Track related Info
574  fxStart = tracks[bestidx]->Start().X();
575  fyStart = tracks[bestidx]->Start().Y();
576  fxEnd = tracks[bestidx]->Stop().X();
577  fyEnd = tracks[bestidx]->Stop().Y();
578 
579  fxDir = tracks[bestidx]->Dir().X();
580  fyDir = tracks[bestidx]->Dir().Y();
581  fzDir = tracks[bestidx]->Dir().Z();
582 
583  unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();
584  // loop over trajectory points starting from the end
585  for(int itrj = ntraj-1; itrj >= 0; --itrj){
586  TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
587 
588  } // for itrj
589 
590  unsigned int minPlane = tracks[bestidx]->MinPlane();
591  unsigned int maxPlane = tracks[bestidx]->MaxPlane();
592  int maxXPlane = -1;
593  int maxYPlane = -1;
594  double hitCosx[tracks[bestidx]->NXCell()];
595  double hitCosy[tracks[bestidx]->NYCell()];
596 
597  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NXCell(); ++ihit){
598  if(tracks[bestidx]->XCell(ihit)->Plane() > maxXPlane){ maxXPlane = tracks[bestidx]->XCell(ihit)->Plane(); }
599  // get direction with this x hit
600  double xyzplane[3];
601  geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
602  double dz = geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->HalfD();
603 
604  double minz = xyzplane[2]-dz;
605  double maxz = xyzplane[2]+dz;
606  unsigned int currtrj = 0;
607  // find the matching trajectory point for this plane;
608  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
609  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
610  } // for itaraj
611  TVector3 direction = tracks[bestidx]->Dir();
612  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
613  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
614  }
615  direction = direction.Unit();
616  hitCosx[ihit] = sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
617  } // for ihit XCell
618 
619 
620 
621  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NYCell(); ++ihit){
622  if(tracks[bestidx]->YCell(ihit)->Plane() > maxYPlane){ maxYPlane = tracks[bestidx]->YCell(ihit)->Plane(); }
623  // get direction with this y hit
624  double xyzplane[3];
625  geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
626  double dz = geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->HalfD();
627  double minz = xyzplane[2]-dz;
628  double maxz = xyzplane[2]+dz;
629  unsigned int currtrj = 0;
630  // find the matching trajectory point for this plane;
631  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
632  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
633  }
634  TVector3 direction = tracks[bestidx]->Dir();
635  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
636  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
637  }
638  direction = direction.Unit();
639  hitCosy[ihit] = sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
640  } // for ihit YCell
641 
642  ///
643 
644  // calculate dedx info for this track
645  std::vector<double> dedxs;
646  std::vector<double> des;
647  std::vector<double> pecorrs;
648  std::vector<double> mips;
649  std::vector<double> dxs;
650  std::vector<double> deads;
651  std::vector<double> xdedxs;
652  std::vector<unsigned int> dedxpls;
653 
654  std::vector<double> pes;
655  std::vector<double> adcs;
656  std::vector<double> xhits;
657  std::vector<double> yhits;
658  std::vector<double> zhits;
659  std::vector<double> thits;
660  std::vector<int> views;
661  std::vector<int> cells;
662  std::vector<int> planes;
663 
664  std::vector<double> mcdes;
665  std::vector<double> mcdebs;
666  std::vector<double> mcdxs;
667 
668  for(unsigned int iPlane = minPlane;iPlane<=maxPlane;++iPlane){
669  double planeGeV = 0;
670  double planePECorr = 0;
671  double planeMIP = 0;
672 
673  int planeView = -1;
674  int plane = -1;
675  int cell = -1;
676  double planeX = 0;
677  double planeY = 0;
678  double planeZ = 0;
679  double planeT = 0;
680  double planeADC = 0;
681  double planePE = 0;
682 
683  double mcPlaneE = 0;
684  double mcPlaneEBirks = 0;
685  double mcPathLength = 0;
686  double totcosinview = 0;
687  double avgCosOtherView = 1.0;
688  bool useOtherView = true;
689  double nhits = 0;
690  double xyzplane[3];
691  geom->Plane(iPlane)->Cell(0)->GetCenter(xyzplane);
692  // double dz = geom->Plane(iPlane)->Cell(0)->HalfD();
693  // double minz = xyzplane[2]-dz;
694  // double maxz = xyzplane[2]+dz;
695  // int matchpt = 0; // Not used for now
696  int minCell = 10000;
697  int maxCell = 0;
698  double nhitsother = 0;
699  int view = 0;
700  int fbIndex = -1;
701  fFiberBrightness = -1;
702  // find the matching trajectory point for this plane;
703  /* Not Used for now
704  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
705  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){matchpt = itraj;}
706  }
707  */
708 
709  if(geom->Plane(iPlane)->View() == geo::kX){
710  planeView = 0;
711  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
712  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane && tracks[bestidx]->XCell(iHit)->View() != geo::kY){
713  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
714  planeADC+=chit->ADC();
715  planePE +=chit->PE();
716  cell = chit->Cell();
717  plane = chit->Plane();
718 
719  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->XCell(iHit));
720  if(minPlane == maxPlane){
721  rhit = cal->MakeRecoHit(*(tracks[bestidx]->XCell(iHit)), 0.25);
722  }
723  if(rhit.IsCalibrated()){
724  planeGeV+=rhit.GeV();
725  planePECorr+=rhit.PECorr();
726  planeMIP+=rhit.MIP();
727  planeX += rhit.X();
728  planeY += rhit.Y();
729  planeZ += rhit.Z();
730  planeT += rhit.T();
731  totcosinview+=hitCosx[iHit];
732  nhits+=1.0;
733  if(tracks[bestidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->XCell(iHit)->Cell();}
734  if(tracks[bestidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->XCell(iHit)->Cell();}
735 
736  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->XCell(iHit)->Cell() );
737 
738  if(bt->HaveTruthInfo()){
739  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
740  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
741  for(int j = 0; j < (int)flsHits.size(); j++){
742  mcPlaneE += flsHits[j].GetEdep();
743  mcPlaneEBirks += flsHits[j].GetEdepBirks();
744  mcPathLength += flsHits[j].GetTotalPathLength();
745 
746  }//End loop over flshits in hit
747  } // if BackTracker
748 
749  }
750  }
751  } // for iHit
752 
753  // try to get the avgcos from the other view
754  double nhitsbefore = 0;
755  double nhitsafter = 0;
756  double totCosBefore = 0;
757  double totCosAfter = 0;
758  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
759  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane-1){
760  ++nhitsbefore;
761  totCosBefore+=hitCosy[iHit];
762  }
763  else if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane+1){
764  ++nhitsafter;
765  totCosAfter+=hitCosy[iHit];
766  }
767  }
768  if(nhitsbefore > 0 && nhitsafter > 0){
769  avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
770  nhitsother = (nhitsbefore+nhitsafter)/2;
771  }
772  else if(nhitsbefore > 0){
773  avgCosOtherView = totCosBefore/nhitsbefore;
774  nhitsother = nhitsbefore;
775  }
776  else if(nhitsafter > 0){
777  avgCosOtherView = totCosAfter/nhitsafter;
778  nhitsother = nhitsafter;
779  }
780  else{ useOtherView = false; }
781 
782  } else if(geom->Plane(iPlane)->View() == geo::kY){
783  planeView = 1;
784  view = 1;
785  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
786  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane && tracks[bestidx]->YCell(iHit)->View() != geo::kX
787  && tracks[bestidx]->YCell(iHit)->View() != geo::kXorY){
788  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
789  planeADC+=chit->ADC();
790  planePE +=chit->PE();
791  cell = chit->Cell();
792  plane = chit->Plane();
793 
794  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->YCell(iHit));
795  if(minPlane == maxPlane){
796  rhit = cal->MakeRecoHit(*(tracks[bestidx]->YCell(iHit)), 0.25);
797  }
798  if(rhit.IsCalibrated()){
799  planeGeV+=rhit.GeV();
800  planePECorr+=rhit.PECorr();
801  planeMIP+=rhit.MIP();
802  planeX += rhit.X();
803  planeY += rhit.Y();
804  planeZ += rhit.Z();
805  planeT += rhit.T();
806  totcosinview+=hitCosy[iHit];
807  nhits+=1.0;
808  if(tracks[bestidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->YCell(iHit)->Cell();}
809  if(tracks[bestidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->YCell(iHit)->Cell();}
810 
811  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->YCell(iHit)->Cell() );
812 
813  if(bt->HaveTruthInfo()){
814  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
815  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
816  for(int j = 0; j < (int)flsHits.size(); j++){
817  mcPlaneE += flsHits[j].GetEdep();
818  mcPlaneEBirks += flsHits[j].GetEdepBirks();
819  mcPathLength += flsHits[j].GetTotalPathLength();
820 
821  }//End loop over flshits in hit
822  } // if BackTracker
823 
824  } // if isCalibrated
825  } //
826  } // for iHit
827 
828  // try to get the avgcos from the other view
829  double nhitsbefore = 0;
830  double nhitsafter = 0;
831  double totCosBefore = 0;
832  double totCosAfter = 0;
833  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
834  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane-1){
835  ++nhitsbefore;
836  totCosBefore+=hitCosx[iHit];
837  }
838  else if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane+1){
839  ++nhitsafter;
840  totCosAfter+=hitCosx[iHit];
841  }
842  }
843  if(nhitsbefore > 0 && nhitsafter > 0){
844  avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
845  nhitsother = (nhitsbefore+nhitsafter)/2;
846  }
847  else if(nhitsbefore > 0){
848  avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
849  nhitsother = nhitsbefore;
850  }
851  else if(nhitsafter > 0){
852  avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
853  nhitsother = nhitsbefore;
854  }
855  else{ useOtherView = false; }
856 
857  } // if kY
858  fnCellsPerPlane = nhits;
859  if(nhits >0){
860 
861  if(nhits == 1) {
862  fFiberBrightness = fbIndex;
863  }
864 
865  double avgCos = TMath::Abs(totcosinview/nhits);
866  double xyzp[3];
867  geom->Plane(iPlane)->Cell(minCell)->GetCenter(xyzp);
868  double xyzp1[3];
869  geom->Plane(iPlane)->Cell(maxCell)->GetCenter(xyzp1);
870  double planewidth = 2*geom->Plane(iPlane)->Cell(1)->HalfD();
871  // double planewidth = 2*HalfD;
872  // calculate the amount of dead material to subtract off
873  // first add up live material
874  double liveLength = 0;
875  for(int i = minCell; i<maxCell+1;++i){
876  // get the cell half width
877  const geo::PlaneGeo* p = geom->Plane(iPlane);
878  const geo::CellGeo* c = p->Cell(i);
879  if(i == minCell || i == maxCell){ liveLength+=c->HalfW(); }
880  else{ liveLength+=(2.0*c->HalfW() ); }
881  }
882  // dead material is total length - live length only if min and max cell are not the same
883  double deadLength = 0;
884  if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
885  double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
886  double planePathLength = TMath::Sqrt(planewidth*planewidth + planeheight*planeheight);
887  // std::cout << " ++++ " << iPlane << " " << planePathLength << " " << planewidth << " " << planeheight << std::endl;
888  // in the case of completely vertical track be a bit careful
889  if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }
890  if(useOtherView){
891  double deltaz = planewidth;
892  double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
893  if(avgCos == 0 || avgCos!=avgCos){
894  deltav = liveLength;
895  // if the track looks completely vertical in both views,
896  // than no part of the track length should be in z direction
897  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
898  }
899 
900  // in opposite view dead length
901  double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
902  // in the case that the other view looks completely vertical scale up
903  // the length between cells in this view by the number of hits in the other view
904  // this case should really never happen
905  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
906  deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
907  }
908  planePathLength = TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
909  }
910  // also the case of completely vertical track
911  if(minPlane == maxPlane){ planePathLength = liveLength; }
912 
913  // std::cout << " ==> " << iPlane << " " << planeGeV << " " << planePathLength << std::endl;
914  assert(planeGeV/(planePathLength)>=0);
915  dedxs.push_back(planeGeV/(planePathLength));
916  des.push_back(planeGeV);
917  pecorrs.push_back(planePECorr);
918  mips.push_back(planeMIP);
919  dxs.push_back(planePathLength);
920  deads.push_back(deadLength);
921  xhits.push_back(planeX/nhits);
922  yhits.push_back(planeY/nhits);
923  zhits.push_back(planeZ/nhits);
924  thits.push_back(planeT/nhits);
925  views.push_back(planeView);
926  planes.push_back(plane);
927  cells.push_back(cell);
928  pes.push_back(planePE);
929  adcs.push_back(planeADC);
930 
931  mcdes.push_back(mcPlaneE);
932  mcdebs.push_back(mcPlaneEBirks);
933  mcdxs.push_back(mcPathLength);
934 
935  } //nhits > 0
936  } // iPlane
937  //.....
938 
940 
941  //Filling proton track length properties
942  fProtonzStart = tracks[protonidx]->Start().Z();
943  fProtonzEnd = tracks[protonidx]->Stop().Z();
944  fProtonTrackLength = tracks[protonidx]->TotalLength();
945 
946  //Proton Track related Info
947  fProtonxStart = tracks[protonidx]->Start().X();
948  fProtonyStart = tracks[protonidx]->Start().Y();
949  fProtonxEnd = tracks[protonidx]->Stop().X();
950  fProtonyEnd = tracks[protonidx]->Stop().Y();
951 
952  fProtonxDir = tracks[protonidx]->Dir().X();
953  fProtonyDir = tracks[protonidx]->Dir().Y();
954  fProtonzDir = tracks[protonidx]->Dir().Z();
955 
956 
957  unsigned int minPlanep = tracks[protonidx]->MinPlane();
958  unsigned int maxPlanep = tracks[protonidx]->MaxPlane();
959  int maxXPlanep = -1;
960  int maxYPlanep = -1;
961  double hitCosxp[tracks[protonidx]->NXCell()];
962  double hitCosyp[tracks[protonidx]->NYCell()];
963 
964  for(unsigned int ihit = 0; ihit<tracks[protonidx]->NXCell(); ++ihit){
965  if(tracks[protonidx]->XCell(ihit)->Plane() > maxXPlanep){ maxXPlanep = tracks[protonidx]->XCell(ihit)->Plane(); }
966  // get direction with this x hit
967  double xyzplane[3];
968  geom->Plane(tracks[protonidx]->XCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
969  double dz = geom->Plane(tracks[protonidx]->XCell(ihit)->Plane())->Cell(0)->HalfD();
970 
971  double minz = xyzplane[2]-dz;
972  double maxz = xyzplane[2]+dz;
973  unsigned int currtrj = 0;
974  // find the matching trajectory point for this plane;
975  for(unsigned int itraj = 0;itraj<tracks[protonidx]->NTrajectoryPoints();++itraj){
976  if(tracks[protonidx]->TrajectoryPoint(itraj).Z()>minz && tracks[protonidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
977  } // for itaraj
978  TVector3 direction = tracks[protonidx]->Dir();
979  if(currtrj>0 && currtrj<tracks[protonidx]->NTrajectoryPoints()-2){
980  direction = tracks[protonidx]->TrajectoryPoint(currtrj) - tracks[protonidx]->TrajectoryPoint(currtrj-1);
981  }
982  direction = direction.Unit();
983  hitCosxp[ihit] = sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
984  } // for ihit XCell
985 
986 
987 
988  for(unsigned int ihit = 0; ihit<tracks[protonidx]->NYCell(); ++ihit){
989  if(tracks[protonidx]->YCell(ihit)->Plane() > maxYPlanep){ maxYPlanep = tracks[protonidx]->YCell(ihit)->Plane(); }
990  // get direction with this y hit
991  double xyzplane[3];
992  geom->Plane(tracks[protonidx]->YCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
993  double dz = geom->Plane(tracks[protonidx]->YCell(ihit)->Plane())->Cell(0)->HalfD();
994  double minz = xyzplane[2]-dz;
995  double maxz = xyzplane[2]+dz;
996  unsigned int currtrj = 0;
997  // find the matching trajectory point for this plane;
998  for(unsigned int itraj = 0;itraj<tracks[protonidx]->NTrajectoryPoints();++itraj){
999  if(tracks[protonidx]->TrajectoryPoint(itraj).Z()>minz && tracks[protonidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
1000  }
1001  TVector3 direction = tracks[protonidx]->Dir();
1002  if(currtrj>0 && currtrj<tracks[protonidx]->NTrajectoryPoints()-2){
1003  direction = tracks[protonidx]->TrajectoryPoint(currtrj) - tracks[protonidx]->TrajectoryPoint(currtrj-1);
1004  }
1005  direction = direction.Unit();
1006  hitCosyp[ihit] = sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
1007  } // for ihit YCell
1008 
1009  ///
1010 
1011  // calculate dedx info for proton track
1012  std::vector<double> dedxp;
1013  std::vector<double> dep;
1014  std::vector<double> pecorrp;
1015  std::vector<double> mipp;
1016  std::vector<double> dxp;
1017  std::vector<double> deadp;
1018  std::vector<double> xdedxp;
1019  std::vector<unsigned int> dedxplp;
1020  std::vector<double> xhitp;
1021  std::vector<double> yhitp;
1022  std::vector<double> zhitp;
1023  std::vector<double> thitp;
1024  std::vector<int> viewp;
1025  std::vector<int> planep;
1026  std::vector<int> cellp;
1027  std::vector<double> adcp;
1028  std::vector<double> pep;
1029  std::vector<double> mcdep;
1030  std::vector<double> mcdebp;
1031  std::vector<double> mcdxp;
1032 
1033  for(unsigned int iPlane = minPlanep;iPlane<=maxPlanep;++iPlane){
1034  double planeGeV = 0;
1035  double mcPlaneE = 0;
1036  double mcPlaneEBirks = 0;
1037  double mcPathLength = 0;
1038  double planePECorr = 0;
1039  double planeMIP = 0;
1040  double planeX = 0;
1041  double planeY = 0;
1042  double planeZ = 0;
1043  double planeT = 0;
1044  int planeView = -1;
1045  int plane = -1;
1046  int cell = -1;
1047  double planeADC = 0;
1048  double planePE = 0;
1049 
1050  double totcosinview = 0;
1051  double avgCosOtherView = 1.0;
1052  bool useOtherView = true;
1053  double nhits = 0;
1054  double xyzplane[3];
1055  geom->Plane(iPlane)->Cell(0)->GetCenter(xyzplane);
1056  // double dz = geom->Plane(iPlane)->Cell(0)->HalfD();
1057  // double minz = xyzplane[2]-dz;
1058  // double maxz = xyzplane[2]+dz;
1059  // int matchpt = 0; // Not Used for now
1060  int minCell = 10000;
1061  int maxCell = 0;
1062  double nhitsother = 0;
1063  int view = 0;
1064  int fbIndex = -1;
1066 
1067 
1068  // find the matching trajectory point for this plane;
1069  /* Not Used for now
1070  for(unsigned int itraj = 0;itraj<tracks[protonidx]->NTrajectoryPoints();++itraj){
1071  if(tracks[protonidx]->TrajectoryPoint(itraj).Z()>minz && tracks[protonidx]->TrajectoryPoint(itraj).Z()<maxz){matchpt = itraj;}
1072  }
1073  */
1074 
1075  if(geom->Plane(iPlane)->View() == geo::kX){
1076  planeView = 0;
1077  for(unsigned int iHit = 0; iHit<tracks[protonidx]->NXCell();++iHit){
1078  if(tracks[protonidx]->XCell(iHit)->Plane() == iPlane && tracks[protonidx]->XCell(iHit)->View() != geo::kY){
1079  const art::Ptr<rb::CellHit>& chit = tracks[protonidx]->XCell(iHit);
1080  planeADC+=chit->ADC();
1081  planePE +=chit->PE();
1082  cell =chit->Cell();
1083  plane =chit->Plane();
1084 
1085  rb::RecoHit rhit = tracks[protonidx]->RecoHit(tracks[protonidx]->XCell(iHit));
1086  if(minPlanep == maxPlanep){
1087  rhit = cal->MakeRecoHit(*(tracks[protonidx]->XCell(iHit)), 0.25);
1088  }
1089  if(rhit.IsCalibrated()){
1090  planeGeV+=rhit.GeV();
1091  planePECorr+=rhit.PECorr();
1092  planeMIP+=rhit.MIP();
1093  planeX += rhit.X();
1094  planeY += rhit.Y();
1095  planeZ += rhit.Z();
1096  planeT += rhit.T();
1097  totcosinview+=hitCosxp[iHit];
1098  nhits+=1.0;
1099  if(tracks[protonidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[protonidx]->XCell(iHit)->Cell();}
1100  if(tracks[protonidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[protonidx]->XCell(iHit)->Cell();}
1101 
1102  fbIndex = fb->getBrightnessIndex(iPlane, tracks[protonidx]->XCell(iHit)->Cell() );
1103  if(bt->HaveTruthInfo()){
1104  const art::Ptr<rb::CellHit>& chit = tracks[protonidx]->XCell(iHit);
1105  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
1106  for(int j = 0; j < (int)flsHits.size(); j++){
1107  mcPlaneE += flsHits[j].GetEdep();
1108  mcPlaneEBirks += flsHits[j].GetEdepBirks();
1109  mcPathLength += flsHits[j].GetTotalPathLength();
1110  }//End loop over flshits in hit
1111  } // if BackTracker
1112 
1113  }
1114  }
1115  } // for iHit
1116 
1117  // try to get the avgcos from the other view
1118  double nhitsbefore = 0;
1119  double nhitsafter = 0;
1120  double totCosBefore = 0;
1121  double totCosAfter = 0;
1122  for(unsigned int iHit = 0; iHit<tracks[protonidx]->NYCell();++iHit){
1123  if(tracks[protonidx]->YCell(iHit)->Plane() == iPlane-1){
1124  ++nhitsbefore;
1125  totCosBefore+=hitCosyp[iHit];
1126  }
1127  else if(tracks[protonidx]->YCell(iHit)->Plane() == iPlane+1){
1128  ++nhitsafter;
1129  totCosAfter+=hitCosyp[iHit];
1130  }
1131  }
1132  if(nhitsbefore > 0 && nhitsafter > 0){
1133  avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
1134  nhitsother = (nhitsbefore+nhitsafter)/2;
1135  }
1136  else if(nhitsbefore > 0){
1137  avgCosOtherView = totCosBefore/nhitsbefore;
1138  nhitsother = nhitsbefore;
1139  }
1140  else if(nhitsafter > 0){
1141  avgCosOtherView = totCosAfter/nhitsafter;
1142  nhitsother = nhitsafter;
1143  }
1144  else{ useOtherView = false; }
1145 
1146  } else if(geom->Plane(iPlane)->View() == geo::kY){
1147  planeView = 1;
1148  view = 1;
1149  for(unsigned int iHit = 0; iHit<tracks[protonidx]->NYCell();++iHit){
1150  if(tracks[protonidx]->YCell(iHit)->Plane() == iPlane && tracks[protonidx]->YCell(iHit)->View() != geo::kX
1151  && tracks[protonidx]->YCell(iHit)->View() != geo::kXorY){
1152  const art::Ptr<rb::CellHit>& chit = tracks[protonidx]->YCell(iHit);
1153  planeADC+=chit->ADC();
1154  planePE +=chit->PE();
1155  cell =chit->Cell();
1156  plane =chit->Plane();
1157 
1158  rb::RecoHit rhit = tracks[protonidx]->RecoHit(tracks[protonidx]->YCell(iHit));
1159  if(minPlanep == maxPlanep){
1160  rhit = cal->MakeRecoHit(*(tracks[protonidx]->YCell(iHit)), 0.25);
1161  }
1162  if(rhit.IsCalibrated()){
1163  planeGeV+=rhit.GeV();
1164  planePECorr+=rhit.PECorr();
1165  planeMIP+=rhit.MIP();
1166  planeX += rhit.X();
1167  planeY += rhit.Y();
1168  planeZ += rhit.Z();
1169  planeT += rhit.T();
1170  totcosinview+=hitCosyp[iHit];
1171  nhits+=1.0;
1172  if(tracks[protonidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[protonidx]->YCell(iHit)->Cell();}
1173  if(tracks[protonidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[protonidx]->YCell(iHit)->Cell();}
1174 
1175  fbIndex = fb->getBrightnessIndex(iPlane, tracks[protonidx]->YCell(iHit)->Cell() );
1176  if(bt->HaveTruthInfo()){
1177  const art::Ptr<rb::CellHit>& chit = tracks[protonidx]->YCell(iHit);
1178  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
1179  for(int j = 0; j < (int)flsHits.size(); j++){
1180  mcPlaneE += flsHits[j].GetEdep();
1181  mcPlaneEBirks += flsHits[j].GetEdepBirks();
1182  mcPathLength += flsHits[j].GetTotalPathLength();
1183  }//End loop over flshits in hit
1184  // std::cout << mcCellE << " " << planeGeV << std::endl;
1185  } // if BackTracker
1186 
1187  } // if isCalibrated
1188  } //
1189  } // for iHit
1190 
1191  // try to get the avgcos from the other view
1192  double nhitsbefore = 0;
1193  double nhitsafter = 0;
1194  double totCosBefore = 0;
1195  double totCosAfter = 0;
1196  for(unsigned int iHit = 0; iHit<tracks[protonidx]->NXCell();++iHit){
1197  if(tracks[protonidx]->XCell(iHit)->Plane() == iPlane-1){
1198  ++nhitsbefore;
1199  totCosBefore+=hitCosxp[iHit];
1200  }
1201  else if(tracks[protonidx]->XCell(iHit)->Plane() == iPlane+1){
1202  ++nhitsafter;
1203  totCosAfter+=hitCosxp[iHit];
1204  }
1205  }
1206  if(nhitsbefore > 0 && nhitsafter > 0){
1207  avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
1208  nhitsother = (nhitsbefore+nhitsafter)/2;
1209  }
1210  else if(nhitsbefore > 0){
1211  avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
1212  nhitsother = nhitsbefore;
1213  }
1214  else if(nhitsafter > 0){
1215  avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
1216  nhitsother = nhitsbefore;
1217  }
1218  else{ useOtherView = false; }
1219 
1220  } // if kY
1221  fProtonnCellsPerPlane = nhits;
1222  if(nhits >0){
1223 
1224  if( nhits == 1 ) fProtonFiberBrightness = fbIndex;
1225 
1226  double avgCos = TMath::Abs(totcosinview/nhits);
1227  double xyzp[3];
1228  geom->Plane(iPlane)->Cell(minCell)->GetCenter(xyzp);
1229  double xyzp1[3];
1230  geom->Plane(iPlane)->Cell(maxCell)->GetCenter(xyzp1);
1231  double planewidth = 2*geom->Plane(iPlane)->Cell(1)->HalfD();
1232 
1233  // calculate the amount of dead material to subtract off
1234  // first add up live material
1235  double liveLength = 0;
1236  for(int i = minCell; i<maxCell+1;++i){
1237  // get the cell half width
1238  const geo::PlaneGeo* p = geom->Plane(iPlane);
1239  const geo::CellGeo* c = p->Cell(i);
1240  if(i == minCell || i == maxCell){ liveLength+=c->HalfW(); }
1241  else{ liveLength+=(2.0*c->HalfW()); }
1242  }
1243  // dead material is total length - live length only if min and max cell are not the same
1244  double deadLength = 0;
1245  if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
1246  double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
1247  double planePathLength = TMath::Sqrt(planewidth*planewidth + planeheight*planeheight);
1248  // std::cout << " ++++ " << iPlane << " " << planePathLength << " " << planewidth << " " << planeheight << std::endl;
1249  // in the case of completely vertical track be a bit careful
1250  if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }
1251  if(useOtherView){
1252  double deltaz = planewidth;
1253  double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
1254  if(avgCos == 0 || avgCos!=avgCos){
1255  deltav = liveLength;
1256  // if the track looks completely vertical in both views,
1257  // than no part of the track length should be in z direction
1258  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
1259  }
1260 
1261  // in opposite view dead length
1262  double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
1263  // in the case that the other view looks completely vertical scale up
1264  // the length between cells in this view by the number of hits in the other view
1265  // this case should really never happen
1266  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
1267  deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
1268  }
1269  planePathLength = TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
1270  }
1271  // also the case of completely vertical track
1272  if(minPlanep == maxPlanep){ planePathLength = liveLength; }
1273 
1274  // std::cout << " ==> " << iPlane << " " << planeGeV << " " << planePathLength << std::endl;
1275  assert(planeGeV/(planePathLength)>=0);
1276  dedxp.push_back(planeGeV/(planePathLength));
1277  dep.push_back(planeGeV);
1278  pecorrp.push_back(planePECorr);
1279  mipp.push_back(planeMIP);
1280  dxp.push_back(planePathLength);
1281  deadp.push_back(deadLength);
1282  xhitp.push_back(planeX/nhits);
1283  yhitp.push_back(planeY/nhits);
1284  zhitp.push_back(planeZ/nhits);
1285  thitp.push_back(planeT/nhits);
1286  viewp.push_back(planeView);
1287  planep.push_back(plane);
1288  cellp.push_back(cell);
1289  adcp.push_back(planeADC);
1290  pep.push_back(planePE);
1291 
1292  mcdep.push_back(mcPlaneE);
1293  mcdebp.push_back(mcPlaneEBirks);
1294  mcdxp.push_back(mcPathLength);
1295  } //nhits > 0
1296  } // iPlane
1297 
1298 
1299  // ..... .....
1300  art::Ptr<remid::ReMId> pid = foids.at(bestidx);
1301  art::Ptr<remid::ReMId> protonpid = foids.at(protonidx);
1302 
1303  //Looking at numu energy things
1304  numue::NumuE numuE;
1305  if (!fndmnyNumuEnergy.isValid()){ continue; }
1306  std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
1307  if (energyVec.empty()){ continue; }
1308  numuE = *energyVec[0];
1309 
1310  //Require decent energy object; not too much hadronic energy in muon catcher
1311  if (numuE.E() == -5){ continue; }
1312  if (!(numuE.NDHadCalCat() + numuE.NDHadCalTran() < 0.03)){ continue; }
1313  //Not much hadronic energy contamination on muon track
1314  if (numuE.HadTrkE() < 0.025){ fPassHadEOnMuonTrack = true; }
1315 
1316  //Looking at qepid things
1317  qeef::QePId qepid;
1318  if (!fndmnyQepid.isValid()){ continue; }
1319  std::vector<art::Ptr<qeef::QePId> > qVec = fndmnyQepid.at(iSlice);
1320  if (qVec.empty()){ continue; }
1321  qepid = *qVec[0];
1322 
1323  //Require ratio of average dE/dx of two tracks is good
1324  if (qepid.Dedx() > 1.75){ fPassDedxRatio = true; }
1325  //Require Ztest of two methods of calculating QE neutrino energy is good
1326  if (qepid.EdiffZ() > -1.3){ fPassEDiffZTest = true; }
1327 
1328 
1329  //Slice energy mostly on two tracks
1330  double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy()+tracks[protonidx]->CalorimetricEnergy());
1331  if (calEDiff < 0.1){ fPassOffTrackEnergy = true; }
1332 
1333  //Slice hits mostly on two tracks
1334  int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell()+tracks[protonidx]->NCell());
1335  if (hitDiff < 4){ fPassOffTrackHits = true; }
1336 
1337  //Two tracks start near each other
1338  double diffX = fabs(tracks[bestidx]->Start().X()-tracks[protonidx]->Start().X());
1339  double diffY = fabs(tracks[bestidx]->Start().Y()-tracks[protonidx]->Start().Y());
1340  double diffZ = fabs(tracks[bestidx]->Start().Z()-tracks[protonidx]->Start().Z());
1341  if (diffX < 6.0){ fPassVertX = true; }
1342  if (diffY < 6.0){ fPassVertY = true; }
1343  if (diffZ < 10.0){ fPassVertZ = true; }
1344 
1345  //Require proton to have a low ReMId value
1346  if (protonpid->Value() < 0.15){ fPassProtonRemid = true; }
1347 
1348  double trueProtonE = -1.0;
1349 
1350  //Check if true signal for sim events
1351  fTrueNumuCCQE = false;
1352  fTrueNumuCCWithProton = false;
1354 
1355  fpdg = 0;
1356  fProtonpdg = 0;
1357 
1358 
1359  if (bt->HaveTruthInfo()){
1360  std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),allhits);
1361  if (!funcReturn.empty()){
1362  const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();
1363 
1364  const std::vector<const sim::Particle*> test_mu = bt->HitsToParticle(tracks[bestidx]->AllCells());
1365  if (!test_mu.empty()) {
1366  fpdg = test_mu[0]->PdgCode();
1367  fTrueE = test_mu[0]->E();
1368  // std::cout << "Mu " << fTrueE << " " << fpdg << std::endl;
1369  }
1370 
1371  const std::vector<const sim::Particle*> test_particle = bt->HitsToParticle(tracks[protonidx]->AllCells());
1372  if (!test_particle.empty()) {
1373  fProtonpdg = test_particle[0]->PdgCode();
1374  fProtonTrueE = test_particle[0]->E();
1375  }
1376 
1377  if((test_neutrino.CCNC()==simb::kCC) && (std::abs(test_neutrino.Nu().PdgCode()) == 14) && (test_neutrino.InteractionType()==simb::kCCQE)){
1378  fTrueNumuCCQE = true;
1379  //Find if matching proton, true E of proton
1380  // const std::vector<const sim::Particle*> test_particle = bt->HitsToParticle(tracks[protonidx]->AllCells());
1381  if (!test_particle.empty()){
1382  int pdg = fProtonpdg;
1383  //std::cout<<" with pdg: "<<pdg<<std::endl;
1384  if ((pdg == 2212)||(pdg == -2212)){
1385  fTrueNumuCCWithProton = true;
1386  // trueProtonE = test_particle[0]->E();
1387  trueProtonE = fProtonTrueE;
1388  //Test to see if matches to 80% efficiency and purity
1389  std::set<int> protonTrackID;
1390  protonTrackID.insert(test_particle[0]->TrackId());
1391  double trackProtonEff = bt->HitCollectionEfficiency(protonTrackID, tracks[protonidx]->AllCells(), allhits, geo::kXorY, 0, false);
1392  double trackProtonPur = bt->HitCollectionPurity(protonTrackID, tracks[protonidx]->AllCells(), 0, 0, false);
1393  if ((trackProtonEff > 0.8)&&(trackProtonPur > 0.8)){
1395  }
1396  }
1397  }
1398  }
1399  }
1400  }
1401 
1402  if(ContainedEvent(slicelist[iSlice],tracks[bestidx])){
1403 
1405  fSlicesPassingSelectionCuts->Fill(0.5);
1406  if (fTrueNumuCCWithProton){
1409  }
1413  }
1414  }
1415  // get all the dedx variables filled
1416  fnDedx = pid->NDedx()-1;
1417  fnDedxUsed = pid->NMeasurementPlanes();
1418  // loop over all of the dE/dx measurements
1419  for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
1420  fDedx = pid->Dedx(idedx);
1421  fVert = pid->DedxVertex(idedx);
1422  fUsed = pid->DedxUsed(idedx);
1423 
1424  fDe = des[idedx];
1425  fPECorr = pecorrs[idedx];
1426  fMIP = mips[idedx];
1427  fDx = dxs[idedx];
1428  fxDead = deads[idedx];
1429  fnewDedx = fDe/fDx;
1430 
1431  fxHit = xhits[idedx];
1432  fyHit = yhits[idedx];
1433  fzHit = zhits[idedx];
1434  ftHit = thits[idedx];
1435  fView = views[idedx];
1436  fPlane = planes[idedx];
1437  fCell = cells[idedx];
1438  fADC = adcs[idedx];
1439  fPE = pes[idedx];
1440 
1441  fmcDe = mcdes[idedx];
1442  fmcDeBirks = mcdebs[idedx];
1443  fmcDx = mcdxs[idedx];
1444 
1445  // get the distance from the end of the track
1446  fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[bestidx]->TotalLength();
1447  // fill the dedx sample tree
1448  DedxSample->Fill();
1449  }//Loop over DedxSample
1450 
1451  // get all the proton dedx variables filled
1452  fProtonnDedx = protonpid->NDedx()-1;
1453  fProtonnDedxUsed = protonpid->NMeasurementPlanes();
1454  // loop over all of the dE/dx measurements
1455  for(unsigned int sidedx = 0; sidedx<protonpid->NDedx(); ++sidedx){
1456  fProtonDedx = protonpid->Dedx(sidedx);
1457  fProtonVert = protonpid->DedxVertex(sidedx);
1458  fProtonUsed = protonpid->DedxUsed(sidedx);
1459  // get the distance from the end of the track
1460  fProtonDistDedx = (1.0 - protonpid->DedxLocation(sidedx))*tracks[protonidx]->TotalLength();
1461 
1462  fProtonDe = dep[sidedx];
1463  fProtonPECorr = pecorrp[sidedx];
1464  fProtonMIP = mipp[sidedx];
1465  fProtonDx = dxp[sidedx];
1466  fProtonxDead = deadp[sidedx];
1468  fProtonxHit = xhitp[sidedx];
1469  fProtonyHit = yhitp[sidedx];
1470  fProtonzHit = zhitp[sidedx];
1471  fProtontHit = thitp[sidedx];
1472  fProtonView = viewp[sidedx];
1473  fProtonPlane = planep[sidedx];
1474  fProtonCell = cellp[sidedx];
1475  fProtonADC = adcp[sidedx];
1476  fProtonPE = pep[sidedx];
1477 
1478  fProtonmcDe = mcdep[sidedx];
1479  fProtonmcDeBirks = mcdebp[sidedx];
1480  fProtonmcDx = mcdxp[sidedx];
1481 
1482  // ........
1483  // fill the dedx sample tree
1484  //std::cout<<"now we fill!"<<std::endl;
1485  ProtonDedxSample->Fill();
1486  }//Loop over DedxSample
1487 
1488  }//If passes contaiment
1489 
1490  } // end loop over slices
1491 
1492  }
1493 
1494 
1495  //.......................................................................
1496 
1498  {
1499  fPOT->Fill(.5, fTotalPOT);
1500  }
1501 
1502 
1503  //.......................................................................
1505  {
1506 
1507  bool contained = false;
1508 
1509  // check that this is a contained event
1510 
1511  // get some useful information here
1512  // 1) how close is the slice to the edge of the detector
1513  unsigned int slcncellsfromedge = 999999999;
1514  for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
1515  const art::Ptr<rb::CellHit>& chit = slice->Cell(ihit);
1516  const int plane = chit->Plane();
1517  unsigned int cellsfromedge = std::min((unsigned int)chit->Cell(), fgeom->Plane(plane)->Ncells() - 1 - chit->Cell());
1518  if(cellsfromedge < slcncellsfromedge){
1519  slcncellsfromedge = cellsfromedge;
1520  }
1521  }
1522  // 2) how close is the track to the edge of the detector
1523  // in planes
1524  //SL - NOPE, need to ask slice, not track
1525  unsigned int firstpl = slice->MinPlane();
1526  unsigned int lastpl = slice->MaxPlane();
1527  // in projected cells, just define directions here actually get projected cell below
1528  TVector3 dirbwk = -trk->Dir();
1529  TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();
1530 
1531  // 3) hits in slice
1532  int sliceNHit = slice->NCell();
1533 
1534  // 4) slice mean time
1535  double slcTime = slice->MeanTNS();
1536 
1537  // 5) hits on track
1538  int trackNHit = trk->NCell();
1539 
1540  //std::cout<<"In contained event loop!"<<std::endl;
1541  // cut the slice number of cells from edge
1542  if((sliceNHit > 20)&&(slcncellsfromedge > 1)&&(slcTime>50000)&&(slcTime<450000)&&(trackNHit > 5)){
1543  if(firstpl > 1 && lastpl < 212){
1544  //std::cout<<"In contained event loop! 1"<<std::endl;
1545  int fwdcell = flivegeom->FullNDProjectedCells(trk->Stop(), dirfwd);
1546  int bwkcell = flivegeom->FullNDProjectedCells(trk->Start(),dirbwk);
1547  // also make sure the track doesn't start too close to the end of the active region
1548  if(fwdcell > 4 && bwkcell > 8 && trk->Start().Z() < 1150){
1549  // if the track ends in the muon catcher, check that it doesn't clip the corner
1550  //std::cout<<"In contained event loop! 2"<<std::endl;
1551  if(trk->Stop().Z() < 1275){
1552  //std::cout<<"In contained event loop! 2a"<<std::endl;
1553  contained = true;
1554  }
1555  else{
1556  //std::cout<<"In contained event loop! 2b"<<std::endl;
1557  /*
1558  // find the trajectory point closest in z to the transition plane
1559  double mindist = 99999999;
1560  unsigned int mintrajid = 0;
1561  for(unsigned int itraj = 0; itraj < trk->NTrajectoryPoints()-1; ++itraj){
1562  double dist = 1275 - trk->TrajectoryPoint(itraj).Z();
1563  if(dist > 0 && dist < mindist){
1564  mindist = dist;
1565  mintrajid = itraj;
1566  }
1567  }
1568  // get the direction at this trajectory point
1569  TVector3 trandir = trk->Dir();
1570  if(mintrajid < trk->NTrajectoryPoints() -2){
1571  TVector3 trandir = (trk->TrajectoryPoint(mintrajid+1) - trk->TrajectoryPoint(mintrajid)).Unit();
1572  }
1573  */
1574  //std::cout<<"In contained event loop! 3"<<std::endl;
1575  //double ypostran = flivegeom->YPositionAtTransition(trk->TrajectoryPoint(mintrajid),trandir);
1576  double ypostranKirk = flivegeom->YPositionAtTransition(trk->Start(),trk->Dir());
1577  if( ypostranKirk < 55){
1578  //std::cout<<"In contained event loop! 4"<<std::endl;
1579  contained = true;
1580  }
1581 
1582  } // end else in muon catcher
1583  } // end if track ends projections contained
1584  } // end if track end planes contained
1585  // } // end if expected number of diblocks
1586  } // end if slice contained
1587 
1588 
1589  return contained;
1590 
1591  }
1592 
1593 
1594 
1595 
1596 }
1597 
1598 
1599 namespace remid{
1600 
1602 
1603 }
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
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
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
float E() const
Definition: Energy.cxx:27
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
art::ServiceHandle< geo::LiveGeometry > flivegeom
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
std::string fPidLabel
Where to find the pids.
art::ServiceHandle< nova::dbi::RunHistoryService > frh
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
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
TTree * ProtonDedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
charged current quasi-elastic
Definition: MCNeutrino.h:96
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
int nmissingdcms
of missing DCMs
Definition: EventQuality.h:23
float abs(float number)
Definition: d0nt_math.hpp:39
art::ServiceHandle< cheat::BackTracker > bt
double Value() const
Definition: PID.h:22
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:31
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
Float_t Y
Definition: plot.C:38
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
Float_t Z
Definition: plot.C:38
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
int InteractionType() const
Definition: MCNeutrino.h:150
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)
std::string fSliceLabel
Where to find reconstructed slices.
float PE() const
Definition: CellHit.h:42
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
art::ServiceHandle< geo::Geometry > fgeom
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
const double j
Definition: BetheBloch.cxx:29
double hornI
kA
Definition: SpillData.h:27
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double posy
mm
Definition: SpillData.h:35
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
float NDHadCalTran() const
Definition: NumuE.cxx:334
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
float GeV() const
Definition: RecoHit.cxx:69
void analyze(art::Event const &evt)
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
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
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
Definition: structs.h:12
float X() const
Definition: RecoHit.h:36
ReMIdDedxStudies(fhicl::ParameterSet const &pset)
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
void beginRun(const art::Run &run)
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
TH1D * fPOT
Histogram of POT passing spill cuts.
std::string fTrackLabel
Where to find recondtructed tracks.
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
art::ServiceHandle< photrans::FiberBrightness > fb
Encapsulate the cell geometry.
Definition: CellGeo.h:25
Event generator information.
Definition: MCNeutrino.h:18
Float_t X
Definition: plot.C:38
double posx
mm
Definition: SpillData.h:34
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
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.