ReMIdDedxFD_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////
2 // \file ReMIdDedxFD_module.cc
3 // \brief A module for creating trees for FD cosmics. Based on the code for ND ReMIdDedx_module.cc
4 // ReMId info is not available for tracks outside the spill window. The pieces of code
5 // to extract the dE and dx info are adapted from RecoMuon_module.cc
6 // \version $Id: ReMIdDedxFD_module.cc,v 2017/02/07 sgermani Exp $
7 // \author Stefano Germani - s.germani@ucl.ac.uk
8 ////////////////////////////////////////////////////////////////////
9 #include <string>
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
63  class ReMIdDedxFD : public art::EDAnalyzer {
64  public:
65  explicit ReMIdDedxFD(fhicl::ParameterSet const& pset);
66  ~ReMIdDedxFD();
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
89 
90 
91  double fTotalPOT;
92  bool fFD; ///< Is detector FD?
93 
94  // variables used in filling the kNN training trees
95  float fTrackLength;
96  //float fProtonTrackLength;
97  //float fDedxSep;
98  //float fScatSep;
99  //float fnuE;
100  //float fvisE;
101  //float frecoE;
102  //float fMeasFrac;
103 
104  TTree* DedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
105 
106  // variables used in filling the sample histogram trees
107  float fDedx;
108  int fnDedx;
110  float fDistDedx;
111  //float fCosScat;
112  //float fDistLastScat;
113  //float fScatVar;
114  //int fnScat;
115  //float fDistScat;
116  //int fnPlane;
117  int fpdg;
118  //int fintType;
119  bool fVert;
120  bool fUsed;
121  float fzStart;
122  float fzEnd;
126 
127  // variables added to trees for Mini Prod mismatch studies (S. Germani) ....
129  float fmyDedx;
130  float fDe;
131  int fView;
132  int fPlane;
133  int fCell;
134  float fPE;
135  float fADC;
136  float fPECorr;
137  float fMIP;
138  float fDx;
139  float fxDead;
140  float fxStart;
141  float fyStart;
142  float fxEnd;
143  float fyEnd;
144  float fxHit;
145  float fyHit;
146  float fzHit;
147  float ftHit;
148  float fxDir;
149  float fyDir;
150  float fzDir;
151  float fmcDe;
152  float fmcDeBirks;
153  float fmcDx;
154 
155  float fTrueE;
157 
158  // ......................
159 
164 
166 
167  //novadaq::cnv::DetId fDetID; ///< Detector ID in nova daq convention namespace typedef
168 
169 
170  //unsigned int fdibfirst = 0;
171  //unsigned int fdiblast = 0;
172 
174 
175 
176  };
177 
178 }
179 
180 
181 namespace remid{
182 
183 
184  //.......................................................................
185 
187  EDAnalyzer(pset)
188  {
189  fHitModuleLabel = (pset.get< std::string >("HitModuleLabel"));
190  fSliceLabel = (pset.get< std::string >("SliceLabel"));
191  fTrackLabel = (pset.get< std::string >("TrackLabel"));
192  // fPidLabel = (pset.get< std::string >("PidLabel" ));
193  fGeneratorLabel = pset.get< std::string >("GeneratorInput");
194  fNumuEnergyLabel = pset.get< std::string >("NumuEnergyInput");
195  // fQepidLabel = pset.get< std::string >("QepidInput");
196  }
197 
198  //.......................................................................
199 
201  {
202  }
203 
204  //.......................................................................
205 
207  {
209 
210  DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
211  DedxSample->Branch("dedx", &fDedx, "dedx/F");
212  DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
213  DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
214  DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
215  DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
216  DedxSample->Branch("vert", &fVert, "vert/O");
217  DedxSample->Branch("used", &fUsed, "used/O");
218  DedxSample->Branch("zStart",&fzStart,"zStart/F");
219  DedxSample->Branch("zEnd",&fzEnd,"zEnd/F");
220  DedxSample->Branch("passOffTrackHits",&fPassOffTrackHits,"passOffTrackHits/O");
221  DedxSample->Branch("passOffTrackEnergy",&fPassOffTrackEnergy,"passOffTrackEnergy/O");
222  DedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");
223 
224  DedxSample->Branch("ncellsPerPlane", &fnCellsPerPlane, "ncellsPerPlane/I");
225  DedxSample->Branch("mydedx", &fmyDedx, "mydedx/F");
226  DedxSample->Branch("de", &fDe, "de/F");
227  DedxSample->Branch("view", &fView, "view/I");
228  DedxSample->Branch("cell", &fCell, "cell/I");
229  DedxSample->Branch("plane", &fPlane, "plane/I");
230  DedxSample->Branch("adc", &fADC, "adc/F");
231  DedxSample->Branch("pe", &fPE, "pe/F");
232  DedxSample->Branch("pecorr", &fPECorr, "pecorr/F");
233  DedxSample->Branch("mip", &fMIP, "mip/F");
234  DedxSample->Branch("dx", &fDx, "dx/F");
235  DedxSample->Branch("deadx", &fxDead, "deadx/F");
236  DedxSample->Branch("xStart", &fxStart, "xStart/F");
237  DedxSample->Branch("yStart", &fyStart, "yStart/F");
238  DedxSample->Branch("xEnd", &fxEnd, "xEnd/F");
239  DedxSample->Branch("yEnd", &fyEnd, "yEnd/F");
240  DedxSample->Branch("xHit", &fxHit, "xHit/F");
241  DedxSample->Branch("yHit", &fyHit, "yHit/F");
242  DedxSample->Branch("zHit", &fzHit, "zHit/F");
243  DedxSample->Branch("tHit", &ftHit, "tHit/F");
244  DedxSample->Branch("xDir", &fxDir, "xDir/F");
245  DedxSample->Branch("yDir", &fyDir, "yDir/F");
246  DedxSample->Branch("zDir", &fzDir, "zDir/F");
247  DedxSample->Branch("mcde", &fmcDe, "mcde/F");
248  DedxSample->Branch("mcdeBirks",&fmcDeBirks, "mcdeBirks/F");
249  DedxSample->Branch("mcdx", &fmcDx, "mcdx/F");
250 
251  DedxSample->Branch("pdg", &fpdg, "pdg/I");
252  DedxSample->Branch("trueE", &fTrueE, "trueE/F");
253  DedxSample->Branch("fiberBrightness", &fFiberBrightness, "fiberBrightness/I");
254 
255  fPOT = tfs->make<TH1D>("POT", ";;Total POT",1, 0.0, 1.0);
256  fSlicesPassingSelectionCuts = tfs->make<TH1D>("SlicesPassingSelectionCuts", ";;Slices Passing Selection Cuts",1, 0.0, 1.0);
257  fTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("TrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
258 
259  fTotalPOT = 0.0;
260 
261  }
262 
263  //......................................................................
264 
266 
267  {
268 
269  fFD = false;
270 
271  novadaq::cnv::DetId detID = fgeom->DetId();
272  if (detID == novadaq::cnv::kFARDET) fFD = true;
273 
274  return;
275  }
276 
277 
278  //.......................................................................
279 
281  {
282  //Only written for FD
283  if (!fFD) return;
284 
285  // get a geometry handle **************************************
288  //***************************************************************
289 
290  //Getting event hits
292  evt.getByLabel(fHitModuleLabel, hithdl);
293 
295  for(size_t h = 0; h < hithdl->size(); ++h){
296  art::Ptr<rb::CellHit> hit(hithdl, h);
297  allhits.push_back(hit);
298  }
299 
300  // get the slices out of the event
302  evt.getByLabel(fSliceLabel,slicecol);
303 
304  art::PtrVector<rb::Cluster> slicelist;
305 
306  for(unsigned int i = 0; i<slicecol->size();++i){
307  art::Ptr<rb::Cluster>slice(slicecol,i);
308  slicelist.push_back(slice);
309  }
310 
311  // get assosciations between slices and tracks
312  art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
313 
314  // get assosciations between slices and numu energy
315  art::FindManyP<numue::NumuE> fndmnyNumuEnergy(slicecol, evt, fNumuEnergyLabel);
316 
317  // loop over the slices
318  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
319 
320  // If slice is noise, skip it
321  if(slicelist[iSlice]->IsNoise()){ continue; }
322 
323  // get all of the tracks associated to this slice
324  std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
325  // if no tracks move on
326  if(tracks.size() == 0){ continue; }
327 
328  //Require precisely one track
329  if(tracks.size() != 1){ continue; }
330 
331 
332  //Require Best Track is 3D
333  bool is3D = true;
334 
335  int bestidx = 0;
336  if(!tracks[bestidx]->Is3D()) is3D = false;
337  if (!is3D){ continue; }
338 
339  //Default all pass cuts to false
340  fPassOffTrackHits = false;
341  fPassOffTrackEnergy = false;
342  fPassHadEOnMuonTrack = false;
343 
344  //Determine if tracks come from outside and stop
346 
347  /*
348  double distStartToFront = livegeo->DistToFront( tracks[bestidx]->Start() );
349  double distStartToTop = livegeo->DistToTop( tracks[bestidx]->Start() );
350  double distStartToSide = livegeo->DistToEdgeX( tracks[bestidx]->Start() );
351 
352 
353  //physical starts
354  if (distStartToFront < 0) continue;
355  if (distStartToTop < 0) continue;
356  if (distStartToSide < 0) continue;
357  */
358 
359  // endpoint isn't contained
360  double distEndToEdge = std::min(livegeo->DistToEdgeXY(tracks[bestidx]->Stop()), livegeo->DistToEdgeZ(tracks[bestidx]->Stop()));
361  bool isStopper = true;
362  if (distEndToEdge < 50) isStopper = false;
363 
364  // Process only contained tracks
365  if (!isStopper) continue;
366 
367  //Filling muon track length properties
368  fzStart = tracks[bestidx]->Start().Z();
369  fzEnd = tracks[bestidx]->Stop().Z();
370  fTrackLength = tracks[bestidx]->TotalLength();
371 
372  fxStart = tracks[bestidx]->Start().X();
373  fyStart = tracks[bestidx]->Start().Y();
374  fxEnd = tracks[bestidx]->Stop().X();
375  fyEnd = tracks[bestidx]->Stop().Y();
376 
377  fxDir = tracks[bestidx]->Dir().X();
378  fyDir = tracks[bestidx]->Dir().Y();
379  fzDir = tracks[bestidx]->Dir().Z();
380 
381  double trklen = tracks[bestidx]->TotalLength();
382 
383  unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();
384 
385  // vector holds effective active track length from the end of the track
386  std::vector<double> lengths(ntraj);
387 
388  // loop over trajectory points starting from the end
389  for(int itrj = ntraj-1; itrj >= 0; --itrj){
390  TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
391  if(itrj == (int)ntraj-1) lengths[itrj] = 0.0;
392  else {
393  TVector3 lasttrajpoint = tracks[bestidx]->TrajectoryPoint(itrj+1);
394  lengths[itrj] = (trajpoint-lasttrajpoint).Mag() + lengths[itrj+1];
395  } // else
396  } // for itrj
397 
398  unsigned int minPlane = tracks[bestidx]->MinPlane();
399  unsigned int maxPlane = tracks[bestidx]->MaxPlane();
400  int maxXPlane = -1;
401  int maxYPlane = -1;
402  double hitCosx[tracks[bestidx]->NXCell()];
403  double hitCosy[tracks[bestidx]->NYCell()];
404 
405  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NXCell(); ++ihit){
406  if(tracks[bestidx]->XCell(ihit)->Plane() > maxXPlane){ maxXPlane = tracks[bestidx]->XCell(ihit)->Plane(); }
407  // get direction with this x hit
408  double xyzplane[3];
409  geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
410  double dz = geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->HalfD();
411 
412  double minz = xyzplane[2]-dz;
413  double maxz = xyzplane[2]+dz;
414  unsigned int currtrj = 0;
415  // find the matching trajectory point for this plane;
416  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
417  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
418  } // for itaraj
419  TVector3 direction = tracks[bestidx]->Dir();
420  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
421  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
422  }
423  direction = direction.Unit();
424  hitCosx[ihit] = sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
425  } // for ihit XCell
426 
427  for(unsigned int ihit = 0; ihit<tracks[bestidx]->NYCell(); ++ihit){
428  if(tracks[bestidx]->YCell(ihit)->Plane() > maxYPlane){ maxYPlane = tracks[bestidx]->YCell(ihit)->Plane(); }
429  // get direction with this y hit
430  double xyzplane[3];
431  geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
432  double dz = geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->HalfD();
433  // double dz = HalfD;
434  double minz = xyzplane[2]-dz;
435  double maxz = xyzplane[2]+dz;
436  unsigned int currtrj = 0;
437  // find the matching trajectory point for this plane;
438  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
439  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
440  }
441  TVector3 direction = tracks[bestidx]->Dir();
442  if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
443  direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
444  }
445  direction = direction.Unit();
446  hitCosy[ihit] = sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
447  } // for ihit YCell
448 
449  ///
450 
451  // calculate dedx info for this track
452  std::vector<double> dedxs;
453  std::vector<double> des;
454  std::vector<double> pecorrs;
455  std::vector<double> mips;
456  std::vector<double> dxs;
457  std::vector<double> deads;
458  std::vector<double> xdedxs;
459  std::vector<unsigned int> dedxpls;
460 
461  std::vector<double> pes;
462  std::vector<double> adcs;
463  std::vector<double> xhits;
464  std::vector<double> yhits;
465  std::vector<double> zhits;
466  std::vector<double> thits;
467  std::vector<int> views;
468  std::vector<int> cells;
469  std::vector<int> planes;
470 
471  std::vector<double> mcdes;
472  std::vector<double> mcdebs;
473  std::vector<double> mcdxs;
474 
475  for(unsigned int iPlane = minPlane;iPlane<=maxPlane;++iPlane){
476 
477  double planeGeV = 0;
478  double planePECorr = 0;
479  double planeMIP = 0;
480 
481  int planeView = -1;
482  int plane = -1;
483  int cell = -1;
484  double planeX = 0;
485  double planeY = 0;
486  double planeZ = 0;
487  double planeT = 0;
488  double planeADC = 0;
489  double planePE = 0;
490 
491  double mcPlaneE = 0;
492  double mcPlaneEBirks = 0;
493  double mcPathLength = 0;
494  double totcosinview = 0;
495  double avgCosOtherView = 1.0;
496  bool useOtherView = true;
497  double nhits = 0;
498  double xyzplane[3];
499  geom->Plane(iPlane)->Cell(0)->GetCenter(xyzplane);
500  double dz = geom->Plane(iPlane)->Cell(0)->HalfD();
501  // double dz = HalfD;
502  double minz = xyzplane[2]-dz;
503  double maxz = xyzplane[2]+dz;
504  int matchpt = 0;
505  int minCell = 10000;
506  int maxCell = 0;
507  double nhitsother = 0;
508  int view = 0;
509  int fbIndex = -1;
510  fFiberBrightness = -1;
511  // find the matching trajectory point for this plane;
512  for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
513  if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){matchpt = itraj;}
514  }
515 
516  trklen = lengths[0];
517 
518  // find the track length up to this trajectory point
519  double xlength = trklen - lengths[matchpt];
520 
521 
522  if(minPlane == maxPlane){xlength = trklen/2;}
523  if(tracks[bestidx]->Stop().Y()>tracks[bestidx]->Start().Y()){xlength = (trklen-xlength);}
524 
525  if(geom->Plane(iPlane)->View() == geo::kX){
526  planeView = 0;
527  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
528  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane && tracks[bestidx]->XCell(iHit)->View() != geo::kY){
529  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
530  planeADC+=chit->ADC();
531  planePE +=chit->PE();
532  cell = chit->Cell();
533  plane = chit->Plane();
534  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->XCell(iHit));
535  if(minPlane == maxPlane){
536  rhit = cal->MakeRecoHit(*(tracks[bestidx]->XCell(iHit)), 0.25);
537  }
538  if(rhit.IsCalibrated()){
539  planeGeV+=rhit.GeV();
540  planePECorr+=rhit.PECorr();
541  planeMIP+=rhit.MIP();
542  planeX += rhit.X();
543  planeY += rhit.Y();
544  planeZ += rhit.Z();
545  planeT += rhit.T();
546  totcosinview+=hitCosx[iHit];
547  nhits+=1.0;
548  if(tracks[bestidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->XCell(iHit)->Cell();}
549  if(tracks[bestidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->XCell(iHit)->Cell();}
550 
551  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->XCell(iHit)->Cell() );
552 
553  if(bt->HaveTruthInfo()){
554  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
555  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
556  for(int j = 0; j < (int)flsHits.size(); j++){
557  mcPlaneE += flsHits[j].GetEdep();
558  mcPlaneEBirks += flsHits[j].GetEdepBirks();
559  mcPathLength += flsHits[j].GetTotalPathLength();
560 
561  }//End loop over flshits in hit
562  } // if BackTracker
563 
564  }
565  }
566  } // for iHit
567 
568  // try to get the avgcos from the other view
569  double nhitsbefore = 0;
570  double nhitsafter = 0;
571  double totCosBefore = 0;
572  double totCosAfter = 0;
573  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
574  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane-1){
575  ++nhitsbefore;
576  totCosBefore+=hitCosy[iHit];
577  }
578  else if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane+1){
579  ++nhitsafter;
580  totCosAfter+=hitCosy[iHit];
581  }
582  }
583  if(nhitsbefore > 0 && nhitsafter > 0){
584  avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
585  nhitsother = (nhitsbefore+nhitsafter)/2;
586  }
587  else if(nhitsbefore > 0){
588  avgCosOtherView = totCosBefore/nhitsbefore;
589  nhitsother = nhitsbefore;
590  }
591  else if(nhitsafter > 0){
592  avgCosOtherView = totCosAfter/nhitsafter;
593  nhitsother = nhitsafter;
594  }
595  else{ useOtherView = false; }
596 
597  } else if(geom->Plane(iPlane)->View() == geo::kY){
598  planeView = 1;
599  view = 1;
600  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
601  if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane && tracks[bestidx]->YCell(iHit)->View() != geo::kX
602  && tracks[bestidx]->YCell(iHit)->View() != geo::kXorY){
603  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
604  planeADC+=chit->ADC();
605  planePE +=chit->PE();
606  cell = chit->Cell();
607  plane = chit->Plane();
608 
609  rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->YCell(iHit));
610  if(minPlane == maxPlane){
611  rhit = cal->MakeRecoHit(*(tracks[bestidx]->YCell(iHit)), 0.25);
612  }
613  if(rhit.IsCalibrated()){
614  planeGeV+=rhit.GeV();
615  planePECorr+=rhit.PECorr();
616  planeMIP+=rhit.MIP();
617  planeX += rhit.X();
618  planeY += rhit.Y();
619  planeZ += rhit.Z();
620  planeT += rhit.T();
621  totcosinview+=hitCosy[iHit];
622  nhits+=1.0;
623  if(tracks[bestidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->YCell(iHit)->Cell();}
624  if(tracks[bestidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->YCell(iHit)->Cell();}
625 
626  fbIndex = fb->getBrightnessIndex(iPlane, tracks[bestidx]->YCell(iHit)->Cell() );
627 
628  if(bt->HaveTruthInfo()){
629  const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
630  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
631  for(int j = 0; j < (int)flsHits.size(); j++){
632  mcPlaneE += flsHits[j].GetEdep();
633  mcPlaneEBirks += flsHits[j].GetEdepBirks();
634  mcPathLength += flsHits[j].GetTotalPathLength();
635 
636  }//End loop over flshits in hit
637  } // if BackTracker
638 
639  } // if isCalibrated
640  } //
641  } // for iHit
642 
643  // try to get the avgcos from the other view
644  double nhitsbefore = 0;
645  double nhitsafter = 0;
646  double totCosBefore = 0;
647  double totCosAfter = 0;
648  for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
649  if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane-1){
650  ++nhitsbefore;
651  totCosBefore+=hitCosx[iHit];
652  }
653  else if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane+1){
654  ++nhitsafter;
655  totCosAfter+=hitCosx[iHit];
656  }
657  }
658  if(nhitsbefore > 0 && nhitsafter > 0){
659  avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
660  nhitsother = (nhitsbefore+nhitsafter)/2;
661  }
662  else if(nhitsbefore > 0){
663  avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
664  nhitsother = nhitsbefore;
665  }
666  else if(nhitsafter > 0){
667  avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
668  nhitsother = nhitsbefore;
669  }
670  else{ useOtherView = false; }
671 
672  } // if kY
673  fnCellsPerPlane = nhits;
674  if(nhits >0){
675 
676  if(nhits == 1) {
677  fFiberBrightness = fbIndex;
678  }
679 
680  double avgCos = TMath::Abs(totcosinview/nhits);
681  double xyzp[3];
682  geom->Plane(iPlane)->Cell(minCell)->GetCenter(xyzp);
683  double xyzp1[3];
684  geom->Plane(iPlane)->Cell(maxCell)->GetCenter(xyzp1);
685  double planewidth = 2*geom->Plane(iPlane)->Cell(1)->HalfD();
686  // double planewidth = 2*HalfD;
687  // calculate the amount of dead material to subtract off
688  // first add up live material
689  double liveLength = 0;
690  for(int i = minCell; i<maxCell+1;++i){
691  // get the cell half width
692  const geo::PlaneGeo* p = geom->Plane(iPlane);
693  const geo::CellGeo* c = p->Cell(i);
694  if(i == minCell || i == maxCell){ liveLength+=c->HalfW(); }
695  else{ liveLength+=(2.0*c->HalfW() ); }
696  }
697  // dead material is total length - live length only if min and max cell are not the same
698  double deadLength = 0;
699  if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
700  double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
701  double planePathLength = TMath::Sqrt(planewidth*planewidth + planeheight*planeheight);
702  // std::cout << " ++++ " << iPlane << " " << planePathLength << " " << planewidth << " " << planeheight << std::endl;
703  // in the case of completely vertical track be a bit careful
704  if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }
705  if(useOtherView){
706  double deltaz = planewidth;
707  double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
708  if(avgCos == 0 || avgCos!=avgCos){
709  deltav = liveLength;
710  // if the track looks completely vertical in both views,
711  // than no part of the track length should be in z direction
712  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
713  }
714 
715  // in opposite view dead length
716  double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
717  // in the case that the other view looks completely vertical scale up
718  // the length between cells in this view by the number of hits in the other view
719  // this case should really never happen
720  if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
721  deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
722  }
723  planePathLength = TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
724  }
725  // also the case of completely vertical track
726  if(minPlane == maxPlane){ planePathLength = liveLength; }
727 
728  xdedxs.push_back(xlength);
729 
730  assert(planeGeV/(planePathLength)>=0);
731  dedxs.push_back(planeGeV/(planePathLength));
732  des.push_back(planeGeV);
733  pecorrs.push_back(planePECorr);
734  mips.push_back(planeMIP);
735  dxs.push_back(planePathLength);
736  deads.push_back(deadLength);
737  xhits.push_back(planeX/nhits);
738  yhits.push_back(planeY/nhits);
739  zhits.push_back(planeZ/nhits);
740  thits.push_back(planeT/nhits);
741  views.push_back(planeView);
742  planes.push_back(plane);
743  cells.push_back(cell);
744  pes.push_back(planePE);
745  adcs.push_back(planeADC);
746 
747  mcdes.push_back(mcPlaneE);
748  mcdebs.push_back(mcPlaneEBirks);
749  mcdxs.push_back(mcPathLength);
750 
751  // std::cout << Form(" Dedx %d %2.6f %2.6f %2.6f", iPlane , planeGeV/(planePathLength), planeGeV, planePathLength ) << std::endl;
752  } //nhits > 0
753  } // iPlane
754  //.....
755 
756  //Looking at numu energy things
757  numue::NumuE numuE;
758  if (!fndmnyNumuEnergy.isValid()){ continue; }
759  std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
760  if (energyVec.empty()){ continue; }
761  numuE = *energyVec[0];
762 
763  //Require decent energy object; not too much hadronic energy in muon catcher
764  if (numuE.E() == -5){ continue; }
765  //Not much hadronic energy contamination on muon track
766  if (numuE.HadTrkE() < 0.025){ fPassHadEOnMuonTrack = true; }
767 
768  //Slice energy mostly on one tracks
769  double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy());
770  if (calEDiff < 0.1){ fPassOffTrackEnergy = true; }
771 
772  //Slice hits mostly on one track
773  int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell());
774  if (hitDiff < 4){ fPassOffTrackHits = true; }
775 
776  fpdg = 0;
777 
778  if (bt->HaveTruthInfo()){
779  const std::vector<const sim::Particle*> test_mu = bt->HitsToParticle(tracks[bestidx]->AllCells());
780  if (!test_mu.empty()) {
781  fpdg = test_mu[0]->PdgCode();
782  fTrueE = test_mu[0]->E();
783  }
784  } // if HaveTruthInfo
785 
786  if(ContainedEvent(slicelist[iSlice],tracks[bestidx])){
787  // get all the dedx variables filled
788  fnDedx = des.size() -1; // pid->NDedx()-1;
789  fnDedxUsed = des.size(); // pid->NMeasurementPlanes();
790  // loop over all of the dE/dx measurements
791  //for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
792  for( int idedx = 0; idedx<fnDedx; ++idedx){
793 
794  fDe = des[idedx];
795  fPECorr = pecorrs[idedx];
796  fMIP = mips[idedx];
797  fDx = dxs[idedx];
798  fxDead = deads[idedx];
799  fmyDedx = fDe/fDx;
800 
801  fxHit = xhits[idedx];
802  fyHit = yhits[idedx];
803  fzHit = zhits[idedx];
804  ftHit = thits[idedx];
805  fView = views[idedx];
806  fPlane = planes[idedx];
807  fCell = cells[idedx];
808  fADC = adcs[idedx];
809  fPE = pes[idedx];
810 
811  fmcDe = mcdes[idedx];
812  fmcDeBirks = mcdebs[idedx];
813  fmcDx = mcdxs[idedx];
814 
815  // ............
816  // get the distance from the end of the track
817  // fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[bestidx]->TotalLength();
818  fDistDedx = (1.0 -xdedxs[idedx]/trklen)*tracks[bestidx]->TotalLength();
819  // fill the dedx sample tree
820  DedxSample->Fill();
821  }//Loop over DedxSample
822 
823  }//If passes contaiment (now temporarely always true)
824 
825  } // end loop over slices
826 
827  }
828 
829 
830  //.......................................................................
831 
833  {
834  fPOT->Fill(.5, fTotalPOT);
835  }
836 
837 
838  //.......................................................................
840  {
841 
842  bool contained = false;
843 
844  // check that this is a contained event
845 
846  // 2) how close is the track to the edge of the detector
847  // in planes
848  //SL - NOPE, need to ask slice, not track
849  // unsigned int firstpl = slice->MinPlane();
850  // unsigned int lastpl = slice->MaxPlane();
851  // in projected cells, just define directions here actually get projected cell below
852  TVector3 dirbwk = -trk->Dir();
853  TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();
854 
855  // 3) hits in slice
856  int sliceNHit = slice->NCell();
857 
858  // 4) slice mean time
859  // double slcTime = slice->MeanTNS();
860 
861  // 5) hits on track
862  int trackNHit = trk->NCell();
863 
864  //std::cout<<"In contained event loop!"<<std::endl;
865  // cut the slice number of cells from edge
866  if((sliceNHit > 20)&&(trackNHit > 5)){
867  int fwdcell = flivegeom->ProjectedLiveCellsToEdge(trk->Stop(), dirfwd);
868  // int bwkcell = flivegeom->ProjectedLiveCellsToEdge(trk->Start(),dirbwk);
869 
870  // also make sure the track doesn't start too close to the end of the active region
871  if(fwdcell > 4) contained = true;
872  } // end if track ends projections contained
873 
874  return contained;
875 
876  }// ContainedEvent
877 
878 
879 }
880 
881 
882 namespace remid{
883 
885 
886 }
float T() const
Definition: RecoHit.h:39
size_t NTrajectoryPoints() const
Definition: Track.h:83
back track the reconstruction to the simulation
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
bool fFD
Is detector FD?
art::ServiceHandle< photrans::FiberBrightness > fb
X or Y views.
Definition: PlaneGeo.h:30
unsigned short Plane() const
Definition: CellHit.h:39
double DistToEdgeXY(TVector3 vertex)
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
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
float Z() const
Definition: RecoHit.h:38
double HalfW() const
Definition: CellGeo.cxx:191
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::ServiceHandle< nova::dbi::RunHistoryService > frh
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
Definition: Run.h:31
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
art::ServiceHandle< cheat::BackTracker > bt
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
double DistToEdgeZ(TVector3 vertex)
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
TH1D * fPOT
Histogram of POT passing spill cuts.
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
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
TH1D * fTrueSlicesPassingSelectionCuts
std::string fNumuEnergyLabel
float PE() const
Definition: CellHit.h:42
Collect Geo headers and supply basic geometry functions.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
std::string fSliceLabel
Where to find reconstructed slices.
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
float MIP() const
Definition: RecoHit.cxx:58
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
art::ServiceHandle< geo::LiveGeometry > flivegeom
T * make(ARGS...args) const
float GeV() const
Definition: RecoHit.cxx:69
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
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:324
Definition: event.h:1
float X() const
Definition: RecoHit.h:36
std::string fTrackLabel
Where to find recondtructed tracks.
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
float Mag() const
art::ServiceHandle< geo::Geometry > fgeom
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
void beginRun(const art::Run &run)
ReMIdDedxFD(fhicl::ParameterSet const &pset)
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
std::string fPidLabel
Where to find the pids.
void analyze(art::Event const &evt)
Encapsulate the geometry of one entire detector (near, far, ndos)
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.