MuonTrackHits_module.cc
Go to the documentation of this file.
1 
2 ////////////////////////////////////////////////////////////////////////
3 // Class: MuonTrackHits //
4 // Module Type: filter //
5 // File: MuonTrackHits_module.cc //
6 // //
7 // Created on Thurs Dec 10 2015 by Diana Patricia Mendez //
8 // Based on MuondEdxAna //
9 // //
10 // Added two new TTrees //
11 // *Eliminated the cut to tricells within tracklength < 200cm //
12 // *fTree stores all the hits (hitId: 0=avghit, 1=zhit, 2=tricell) //
13 // *fTree_CalibTrack saves info about all the muon tracks //
14 ////////////////////////////////////////////////////////////////////////
15 
16 // framework includes....
24 #include "fhiclcpp/ParameterSet.h"
28 #include "Geometry/LiveGeometry.h"
30 
31 // novasoft stuff
33 #include "RecoBase/Track.h"
34 #include "RecoBase/CellHit.h"
35 #include "Calibrator/Calibrator.h"
36 #include "NovaTimingUtilities/TimingUtilities.h"
37 
38 // C++
39 #include <memory>
40 #include <string>
41 #include <iostream>
42 #include <vector>
43 
44 // ROOT
45 #include "TH3D.h"
46 #include "TH2D.h"
47 #include "TProfile.h"
48 #include "TProfile2D.h"
49 #include "TFile.h"
50 #include "TVector3.h"
51 #include "TTree.h"
52 
53 namespace caldp {
54  class PCHit;
55 }
56 
57 namespace calib {
58  class MuonTrackHits;
59 
60  class MuonTrackHits : public art::EDFilter {
61  public:
62  explicit MuonTrackHits(fhicl::ParameterSet const & p);
63  virtual ~MuonTrackHits();
64 
65  bool filter(art::Event & e) override;
66  void reconfigure(fhicl::ParameterSet const & p);
67  void beginJob() override;
68  void endJob() override;
69  bool beginRun(art::Run & r) override;
70  bool endRun(art::Run & r) override;
71 
72  private:
73  double getPECorr(art::Event & e,caldp::PCHit const &pchit);
74  double getPECorrToGeV(art::Event & e, caldp::PCHit const &pchit);
75 
76  //keep
77  TH2D* fdEdx;
78  TH2D* fdEdx_minus;
79  TH2D* fdEdx_plus;
80  TH2D* fdEdx_true;
81 
82 
83  TTree* fTree;
85 
86  int diblock;
87  int view;
88  double FlightLength;
90  double x;
91  double xCust;
92  double plane;
93  double path;
94  double trueE;
95  double PE;
96  double PECorr;
97  double tns;
98  uint32_t EventTimeHigh;
99  double cell;
100  double w;
101  double PECorrToGeV;
102  double RecoGeV;
104  int run;
105  int subRun;
106  int event;
107  int evt_year;
109  int evt_day;
110  int evt_hour;
111  int evt_min;
112  int evt_sec;
113  int APD;
116  int nhits;
117  int navghits;
118  int nzhits;
120  int nXYHits;
121  int hitId;
122 
123  bool badTNS;
124 
125  double xStart;
126  double yStart;
127  double zStart;
128  double xEnd;
129  double yEnd;
130  double zEnd;
131 
132  double dCosX;
133  double dCosY;
134  double dCosZ;
135  //double zenith;
136  //double azimuth;
138  int xPlanes;
139  int yPlanes;
141 
147 
148  double fSysShift;
149 
150  // *time_t* of start of Nova epoch, 01-Jan-2010 00:00:00, UTC
151  // This is the value subtracted from the UNIX time_t, timeval or timespec
152  // seconds field.
153  // Using unsigned long long rather than uint64_t to ensure consistent behavior
154  // on OSX and LINUX
155  //const uint32_t NOVA_EPOCH = 1262304000;
156  //const unsigned long long NOVA_TIME_FACTOR = 64000000;
157 
158  };
159 
160 
161  MuonTrackHits::MuonTrackHits(fhicl::ParameterSet const & p)
162  {
163  this->reconfigure(p);
164  produces< TH2D, art::InRun >("dEdx");
165  produces< TH2D, art::InRun >("dEdxMinus");
166  produces< TH2D, art::InRun >("dEdxPlus");
167  produces< TH2D, art::InRun >("dEdxTrue");
168  }
169 
170  MuonTrackHits::~MuonTrackHits()
171  {
172  // Clean up dynamic memory and other resources here.
173  }
174 
175  bool MuonTrackHits::filter(art::Event & e)
176  {
177 
178  run = (int)(e.run());
179  subRun = (int)(e.subRun());
180  event = (int)(e.event());
181 
182 
183  art::Timestamp ts = e.time();
184  EventTimeHigh = ts.timeHigh();
185  const time_t timeSec = ts.timeHigh();
186  struct tm* timeStruct = localtime(&timeSec);
187  evt_year=timeStruct->tm_year+1900;
188  evt_month=timeStruct->tm_mon+1;
189  evt_day=timeStruct->tm_mday;
190  evt_hour=timeStruct->tm_hour+1;
191  evt_min=timeStruct->tm_min;
192  evt_sec=timeStruct->tm_sec;
193 
194  /// Get the reconstructed tracks
196  e.getByLabel(fStopperLabel, tracks);
197  if( tracks->size() == 0 ) return false;
198 
199  art::FindManyP<caldp::PCHit> PCHitGetter(tracks, e, art::InputTag(fPCHitLabel,fQualXYName));
200  art::FindManyP<caldp::PCHit> PCHitAvgGetter(tracks, e, art::InputTag(fPCHitLabel,fQualAvgName));
201  art::FindManyP<caldp::PCHit> PCHitZGetter(tracks, e, art::InputTag(fPCHitLabel,fQualZName));
202  if( !PCHitGetter.isValid() ) return false;
203 
204  // loop over stopping tracks in event
205  for( unsigned int i=0; i<tracks->size(); ++i ){
206 
207  //track level counters:
208  int windowTricells = 0;
209  int windowAvgHits = 0;
210  int windowZHits = 0;
211  int windowPlanes = 0;
212 
213  totalPlanes = 0;
214  xPlanes = 0;
215  yPlanes = 0;
216 
217  std::vector<bool> trackWindowPlanes (1000, false); //vector to prevent double counting planes
218 
219  std::vector<bool> trackPlanes (1000, false); //vector to prevent double counting total number of planes in track
220  std::vector<bool> trackXPlanes (1000, false); //vector to prevent double counting total number of x planes in track
221  std::vector<bool> trackYPlanes (1000, false); //vector to prevent double counting total number of y planes in track
222 
223  bool windowXView = false;
224  bool windowYView = false;
225 
226  //make custom track length to check against inbuilt length
227  TVector3 startXYZ = (*tracks)[i].Start();
228  TVector3 stopXYZ = (*tracks)[i].Stop();
229  double trackLength = (startXYZ-stopXYZ).Mag();
230  xStart = startXYZ.X();
231  yStart = startXYZ.Y();
232  zStart = startXYZ.Z();
233  xEnd = stopXYZ.X();
234  yEnd = stopXYZ.Y();
235  zEnd = stopXYZ.Z();
236  // zenith = startXYZ.Theta();
237  // azimuth = startXYZ.Phi();
238 
239  //cosine giving initial direction of track
240  dCosX = (*tracks)[i].Dir().X()/(*tracks)[i].Dir().Mag();
241  dCosY = (*tracks)[i].Dir().Y()/(*tracks)[i].Dir().Mag();
242  dCosZ = (*tracks)[i].Dir().Z()/(*tracks)[i].Dir().Mag();
243 
244 
245  // needed for track window, but good generally to reject protons
246  // if( (*tracks)[i].TotalLength() < 200. ) continue;
247 
250  int stopPlaneID = 10000;
251  int firstMuCPlane = -1;
252  bool isNearDet = false;
253 
254  //remove tracks that stop outside/ near edge of the active detector. New files won't need this as the selection is in the updated pchitslist_module.
256  distanceToLiveDetEdgeZ = fLiveGeom->DistToEdgeZ((*tracks)[i].Stop());
257 
258  //continue if track ends in the muon catcher
259  if(fGeo->DetId() == novadaq::cnv::kNEARDET){
260  isNearDet = true;
261  cid = fGeo->CellId(stopXYZ[0], stopXYZ[1], stopXYZ[2], 1.5, 1.5, 6., 1.0);
262  if(cid) stopPlaneID = fGeo->getPlaneID(cid);
263  firstMuCPlane = fGeo->FirstPlaneInMuonCatcher();
264  if(stopPlaneID >= firstMuCPlane) continue;
265  }
266 
268 
269  // loop over PCHits associated to stopping track
270  std::vector<art::Ptr<caldp::PCHit> > hits = PCHitGetter.at(i);
271  std::vector<art::Ptr<caldp::PCHit> > avghits = PCHitAvgGetter.at(i);
272  std::vector<art::Ptr<caldp::PCHit> > zhits = PCHitZGetter.at(i);
273 
274  navghits = avghits.size();
275  nzhits = zhits.size();
276  ntricells = hits.size();
277  nhits = hits.size() + avghits.size() + zhits.size();
278 
279  int tplane;
280  int trackView;
281 
282  for(unsigned int k=0; k<hits.size(); ++k){
283 
284  tplane = hits[k]->Plane();
285  trackView = (int)hits[k]->View();
286 
287  if(trackPlanes[tplane] == false) {
288  totalPlanes ++;
289  trackPlanes[tplane] = true;
290  };
291  if(trackView == 0 && trackXPlanes[tplane] == false){
292  xPlanes ++;
293  trackXPlanes[tplane] = true;
294  };
295  if(trackView == 1 && trackYPlanes[tplane] == false){
296  yPlanes ++;
297  trackYPlanes[tplane] = true;
298  };
299 
300  };
301 
302  for(unsigned int k=0; k<avghits.size(); ++k){
303  tplane = avghits[k]->Plane();
304  trackView = (int)avghits[k]->View();
305 
306  if(trackPlanes[tplane] == false){
307  totalPlanes ++;
308  trackPlanes[tplane] = true;
309  };
310  if(trackView == 0 && trackXPlanes[tplane] == false){
311  xPlanes ++;
312  trackXPlanes[tplane] = true;
313  };
314  if(trackView == 1 && trackYPlanes[tplane] == false){
315  yPlanes ++;
316  trackYPlanes[tplane] = true;
317  };
318 
319  };
320 
321  for(unsigned int k=0; k<zhits.size(); ++k){
322  tplane = zhits[k]->Plane();
323  trackView = (int)zhits[k]->View();
324 
325  if(trackPlanes[tplane] == false){
326  totalPlanes ++;
327  trackPlanes[tplane] = true;
328  };
329  if(trackView == 0 && trackXPlanes[tplane] == false){
330  xPlanes ++;
331  trackXPlanes[tplane] = true;
332  };
333  if(trackView == 1 && trackYPlanes[tplane] == false){
334  yPlanes ++;
335  trackYPlanes[tplane] = true;
336  };
337 
338  };
339 
340 
341  /*
342  if(xPlanes==0 || yPlanes==0){
343  std::cout << "***** HEY! *****" << std::endl;
344  std::cout << "run: " << run << ", subRun: " << subRun << ", event: " << event << " xPlanes: " << xPlanes << " yPlanes: " << yPlanes << ", TotalTrackLength: " << TotalTrackLength << std::endl;
345  std::cout << "xStart: " << xStart << ", yStart: " << yStart << ", zStart: " << zStart << std::endl;
346  std::cout << "xEnd: " << xEnd << ", yEnd: " << yEnd << ", zEnd: " << zEnd << std::endl;
347  std::cout << "* * * * * * * * * *" << std::endl;
348  }
349  */
350 
351  int cont;
352  for(unsigned int k=0; k<1000;k++){
353  cont = 0;
354  if(trackPlanes[k] == true){
355  for(unsigned int kk=0;kk<6;kk++){
356  if(trackPlanes[k+kk]==true) cont ++;
357  //else cont = cont;
358  };
359  if (cont>=5) break;
360  };
361  };
362 
363  if(cont>=5) continuity = 1;
364  else continuity = 0;
365 
366 
367  fTree_CalibTracks->Fill();
368 
369  // ---------------- Fill vector containing hits per APD ---------------------- //
370  std::vector<int> APDhits (100,0);
371  std::vector< std::vector<int> > DiblockAPDhits (100, std::vector<int> (100,0));
372  badTNS = false; //temporary handle for showery events in old data
373 
374  //avg. hits
375  if(avghits.size() != 0)
376  for( unsigned int p=0; p<avghits.size(); ++p ) {
377  APDhits[avghits[p]->APD()] += 1;
378  DiblockAPDhits.at(avghits[p]->Diblock()).at(avghits[p]->APD()) += 1;
379  if(avghits[p]->TNS() < 0 || avghits[p]->TNS() > 1400) badTNS = true;
380  FlightLength = avghits[p]->FlightLen();
381  TotalTrackLength = (*tracks)[i].TotalLength();
382  x = trackLength - FlightLength;
383  if(x<100 || x>200) continue;
384  windowAvgHits ++; //track level counters
385  }
386  //z hits
387  if(zhits.size() != 0)
388  for( unsigned int q=0; q<zhits.size(); ++q ){
389  APDhits[zhits[q]->APD()] += 1;
390  DiblockAPDhits.at(zhits[q]->Diblock()).at(zhits[q]->APD()) += 1;
391  if(zhits[q]->TNS() < 0 || zhits[q]->TNS() > 1400) badTNS = true;
392  FlightLength = zhits[q]->FlightLen();
393  TotalTrackLength = (*tracks)[i].TotalLength();
394  x = trackLength - FlightLength;
395  if(x<100 || x>200) continue;
396  windowZHits ++; //track level counters
397  }
398  // tricell hits
399  if(hits.size() != 0)
400  for( unsigned int n=0; n<hits.size(); ++n ){
401  APDhits[hits[n]->APD()] += 1;
402  DiblockAPDhits.at(hits[n]->Diblock()).at(hits[n]->APD()) += 1;
403  if(hits[n]->TNS() < 0 || hits[n]->TNS() > 1400) badTNS = true;
404  FlightLength = hits[n]->FlightLen();
405  TotalTrackLength = (*tracks)[i].TotalLength();
406  x = trackLength - FlightLength;
407  if(x<100 || x>200) continue;
408  windowTricells ++; //track level counters
409  }
410  // ---------------- Filled APD vector -------------------------------------------- //
411 
412 
413  for( unsigned int j=0; j<avghits.size(); ++j) {
414  if(isNearDet)
415  if(avghits[j]->Plane() >= firstMuCPlane) continue;
416 
417  hitId = 0;
418  nXYHits = avghits.size();
419  diblock = avghits[j]->Diblock();
420  view = (int)avghits[j]->View();
421  FlightLength = avghits[j]->FlightLen();
422  TotalTrackLength= (*tracks)[i].TotalLength();
423  x = trackLength - FlightLength;
424  plane = avghits[j]->Plane();
425  path = avghits[j]->Path();
426  trueE = avghits[j]->TrueMeV();
427  PE = avghits[j]->PE();
428  PECorr = this->getPECorr(e,*avghits[j]);
429  PECorrToGeV = this->getPECorrToGeV(e,*avghits[j]);
430  RecoGeV = PECorr * PECorrToGeV;
431  tns = avghits[j]->TNS();
432  cell = avghits[j]->Cell();
433  w = avghits[j]->W();
434  APD = avghits[j]->APD();
435  sameAPDhitsEvent= APDhits[APD];
436  sameAPDhitDiblockEvent= DiblockAPDhits[diblock][APD];
437 
438  fTree->Fill(); //fill ttree
439 
440  }
441 
442 
443  for( unsigned int j=0; j<zhits.size(); ++j) {
444  if(isNearDet)
445  if(zhits[j]->Plane() >= firstMuCPlane) continue;
446 
447  hitId = 1;
448  nXYHits = zhits.size();
449  diblock = zhits[j]->Diblock();
450  view = (int)zhits[j]->View();
451  FlightLength = zhits[j]->FlightLen();
452  TotalTrackLength= (*tracks)[i].TotalLength();
453  x = trackLength - FlightLength;
454  plane = zhits[j]->Plane();
455  path = zhits[j]->Path();
456  trueE = zhits[j]->TrueMeV();
457  PE = zhits[j]->PE();
458  PECorr = this->getPECorr(e,*zhits[j]);
459  PECorrToGeV = this->getPECorrToGeV(e,*zhits[j]);
460  RecoGeV = PECorr * PECorrToGeV;
461  tns = zhits[j]->TNS();
462  cell = zhits[j]->Cell();
463  w = zhits[j]->W();
464  APD = zhits[j]->APD();
465  sameAPDhitsEvent= APDhits[APD];
466  sameAPDhitDiblockEvent= DiblockAPDhits[diblock][APD];
467 
468  fTree->Fill(); //fill ttree
469 
470  }
471 
472 
473  for( unsigned int j=0; j<hits.size(); ++j ) { // loop over number of tricells per track
474 
475  if(isNearDet)
476  if(hits[j]->Plane() >= firstMuCPlane) continue;
477 
478  //test poisson lambda:
479  double poissonLambda = hits[j]->PoissonLambda();
480  std::cout<< "poissonLambda: " << poissonLambda <<std::endl;
481 
482  hitId = 2;
483  nXYHits = hits.size();
484  diblock = hits[j]->Diblock();
485  view = (int)hits[j]->View();
486  FlightLength = hits[j]->FlightLen();
487  TotalTrackLength= (*tracks)[i].TotalLength();
488  x = trackLength - FlightLength;
489  plane = hits[j]->Plane();
490  path = hits[j]->Path();
491  trueE = hits[j]->TrueMeV();
492  PE = hits[j]->PE();
493  PECorr = this->getPECorr(e,*hits[j]);
494  PECorrToGeV = this->getPECorrToGeV(e,*hits[j]);
495  RecoGeV = PECorr * PECorrToGeV;
496  tns = hits[j]->TNS();
497  cell = hits[j]->Cell();
498  w = hits[j]->W();
499  APD = hits[j]->APD();
500  sameAPDhitsEvent= APDhits[APD];
501  sameAPDhitDiblockEvent= DiblockAPDhits[diblock][APD];
502 
503  fTree->Fill(); //fill ttree
504 
505 
506  if(x<100 || x>200) continue;
507 
508  //track level counters
509  if(view == 0) windowXView = true;
510  if(view == 1) windowYView = true;
511  //check if plane already counted
512  if(trackWindowPlanes[plane] == false) windowPlanes ++;
513  trackWindowPlanes[plane] = true;
514 
515 
516  } // end loop over tricells
517 
518 
519  int nViewsInWindow = 0;
520  if(windowXView == true) nViewsInWindow++;
521  if(windowYView == true) nViewsInWindow++;
522 
523 
524 
525  }// loop over tracks
526 
527  return true;
528  }
529 
530 
531  double MuonTrackHits::getPECorr(art::Event & e, caldp::PCHit const &pchit){
532 
533  //Calibrator necessary to get Atten information
535 
536  //Make mock cellhit
537  rb::CellHit cellhit;
538  cellhit.SetPlane(pchit.Plane());
539  cellhit.SetCell(pchit.Cell());
540  cellhit.SetView(pchit.View());
541  cellhit.SetTNS(pchit.TNS(),pchit.GoodTime());
542  if(!e.isRealData())cellhit.SetMC(); // assume only filled for MC...
543  cellhit.SetPE(pchit.PE());
544 
545  double pecorr = calibrator->GetPECorr(cellhit, pchit.W());
546  return pecorr;
547  }
548 
549  double MuonTrackHits::getPECorrToGeV(art::Event & e, caldp::PCHit const &pchit){
550 
551  //Calibrator necessary to get Atten information
553  rb::CellHit cellhit;
554  if(!e.isRealData())cellhit.SetMC(); // assume only filled for MC...
555 
556  double pecorrtogev = calibrator->GetPECorrToGeVScale(cellhit);
557  return pecorrtogev;
558  }
559 
560  void MuonTrackHits::reconfigure(fhicl::ParameterSet const & pset)
561  {
562  fPCHitLabel = pset.get< std::string >("PCHitLabel"); //Label of PCHitList module
563  fQualXYName = pset.get< std::string >("QualXYName"); //Instance label, "XY" quality hits
564  fQualZName = pset.get< std::string >("QualZName"); //Instance label, "Z" quality hits
565  fQualAvgName = pset.get< std::string >("QualAvgName"); //Instance label, "Avg" quality hits
566  fStopperLabel = pset.get<std::string>("StopperLabel");
567  fSysShift = pset.get<double>("SysShift"); // track end point systematic shift
568  }
569 
571  {
572  // histograms for output
574 
575  fTree = tfs->make<TTree>("fTree","tree to hold muon hits within each track variables");
576  fTree->Branch("nXYHits",&nXYHits);
577  fTree->Branch("nhits",&nhits);
578  fTree->Branch("hitId",&hitId);
579  fTree->Branch("xStart",&xStart);
580  fTree->Branch("yStart",&yStart);
581  fTree->Branch("zStart",&zStart);
582  fTree->Branch("dCosX",&dCosX);
583  fTree->Branch("dCosY",&dCosY);
584  fTree->Branch("dCosZ",&dCosZ);
585  // fTree->Branch("zenith",&zenith);
586  // fTree->Branch("azimuth",&azimuth);
587  fTree->Branch("x",&x);
588  fTree->Branch("badTNS",&badTNS);
589  fTree->Branch("run",&run);
590  fTree->Branch("diblock",&diblock);
591  fTree->Branch("view",&view);
592  fTree->Branch("FlightLength",&FlightLength);
593  fTree->Branch("TotalTrackLength",&TotalTrackLength);
594  fTree->Branch("xCust",&xCust);
595  fTree->Branch("distanceToLiveDetEdgeZ",&distanceToLiveDetEdgeZ);
596  fTree->Branch("plane",&plane);
597  fTree->Branch("path",&path);
598  fTree->Branch("trueE",&trueE);
599  fTree->Branch("PE",&PE);
600  fTree->Branch("PECorr",&PECorr);
601  fTree->Branch("tns",&tns);
602  fTree->Branch("cell",&cell);
603  fTree->Branch("w",&w);
604  fTree->Branch("PECorrToGeV",&PECorrToGeV);
605  fTree->Branch("RecoGeV",&RecoGeV);
606  fTree->Branch("EventTimeHigh",&EventTimeHigh);
607  fTree->Branch("evt_year",&evt_year);
608  fTree->Branch("evt_month",&evt_month);
609  fTree->Branch("evt_day",&evt_day);
610  fTree->Branch("evt_hour",&evt_hour);
611  fTree->Branch("evt_min",&evt_min);
612  fTree->Branch("evt_sec",&evt_sec);
613  fTree->Branch("APD",&APD);
614  fTree->Branch("sameAPDhitsEvent",&sameAPDhitsEvent);
615  fTree->Branch("sameAPDhitDiblockEvent",&sameAPDhitDiblockEvent);
616  fTree->Branch("totalPlanes",&totalPlanes);
617  fTree->Branch("xPlanes",&xPlanes);
618  fTree->Branch("yPlanes",&yPlanes);
619  fTree->Branch("continuity",&continuity);
620 
621  fTree_CalibTracks = tfs->make<TTree>("fTree_CalibTracks","tree to hold info about muon tracks");
622  fTree_CalibTracks->Branch("run",&run);
623  fTree_CalibTracks->Branch("subRun",&subRun);
624  fTree_CalibTracks->Branch("event",&event);
625  fTree_CalibTracks->Branch("nhits",&nhits);
626  fTree_CalibTracks->Branch("totalPlanes",&totalPlanes);
627  fTree_CalibTracks->Branch("xPlanes",&xPlanes);
628  fTree_CalibTracks->Branch("yPlanes",&yPlanes);
629  fTree_CalibTracks->Branch("continuity",&continuity);
630  fTree_CalibTracks->Branch("navghits",&navghits);
631  fTree_CalibTracks->Branch("nzhits",&nzhits);
632  fTree_CalibTracks->Branch("ntricells",&ntricells);
633  fTree_CalibTracks->Branch("totalPlanes",&totalPlanes);
634  fTree_CalibTracks->Branch("xStart",&xStart);
635  fTree_CalibTracks->Branch("yStart",&yStart);
636  fTree_CalibTracks->Branch("zStart",&zStart);
637  fTree_CalibTracks->Branch("xEnd",&xEnd);
638  fTree_CalibTracks->Branch("yEnd",&yEnd);
639  fTree_CalibTracks->Branch("zEnd",&zEnd);
640  fTree_CalibTracks->Branch("dCosX",&dCosX);
641  fTree_CalibTracks->Branch("dCosY",&dCosY);
642  fTree_CalibTracks->Branch("dCosZ",&dCosZ);
643  // fTree_CalibTracks->Branch("zenith",&zenith);
644  // fTree_CalibTracks->Branch("azimuth",&azimuth);
645  fTree_CalibTracks->Branch("TotalTrackLength",&TotalTrackLength);
646 
647 
648  fdEdx = tfs->make<TH2D>("dEdx",";Distance to track end (cm);PECorr / cm;",
649  50,0.,500.,200,0.,200.);
650  fdEdx_minus = tfs->make<TH2D>("dEdx_minus",";Distance to track end (cm);PECorr / cm;",
651  50,0.,500.,200,0.,200.);
652  fdEdx_plus = tfs->make<TH2D>("dEdx_plus",";Distance to track end (cm);PECorr / cm;",
653  50,0.,500.,200,0.,200.);
654  fdEdx_true = tfs->make<TH2D>("dEdx_true",";Distance to track end (cm);True MeV / cm;",
655  50,0.,500.,200,0.,5.);
656 
657  }
658 
659  bool MuonTrackHits::beginRun(art::Run &r){
660  return true;
661  }
662 
663  bool MuonTrackHits::endRun(art::Run &r){
664 
666  std::unique_ptr<TH2D> pdEdx(tfs->make<TH2D>(*fdEdx));
667  std::unique_ptr<TH2D> pdEdx_minus(tfs->make<TH2D>(*fdEdx_minus));
668  std::unique_ptr<TH2D> pdEdx_plus(tfs->make<TH2D>(*fdEdx_plus));
669  std::unique_ptr<TH2D> pdEdx_true(tfs->make<TH2D>(*fdEdx_true));
670 
671  r.put(std::move(pdEdx),"dEdx");
672  r.put(std::move(pdEdx_minus),"dEdxMinus");
673  r.put(std::move(pdEdx_plus),"dEdxPlus");
674  r.put(std::move(pdEdx_true),"dEdxTrue");
675  return true;
676  }
677 
678  void MuonTrackHits::endJob()
679  {
680 
681  }
682 
683 }// end namespace
float TNS() const
Return uncorrected hit time.
Definition: PCHit.h:52
SubRunNumber_t subRun() const
Definition: Event.h:72
diblock
print "ROW IS " print row
Definition: geo2elec.py:31
double GetPECorrToGeVScale(rb::CellHit const &cellhit)
void SetTNS(float tns, bool good)
Definition: CellHit.h:56
const char * p
Definition: xmltok.h:285
bool GoodTime() const
Return quality of timing fit for cell.
Definition: PCHit.h:54
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
int Cell() const
Return cell value.
Definition: PCHit.h:26
bool isRealData() const
Definition: Event.h:83
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
"Pre-calibration hit". Common input to calibration procedures
Definition: PCHit.h:16
Definition: Run.h:31
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
void SetView(geo::View_t view)
Definition: CellHit.h:54
void SetPlane(unsigned short plane)
Definition: CellHit.h:53
float W() const
Return W value.
Definition: PCHit.h:42
void SetCell(unsigned short cell)
Definition: CellHit.h:52
CDPStorage service.
void hits()
Definition: readHits.C:15
void SetPE(float pe)
Definition: CellHit.h:55
float PE() const
Return PE value.
Definition: PCHit.h:38
double DistToEdgeZ(TVector3 vertex)
void beginJob()
Encapsulate the geometry of one entire detector (near, far, ndos)
T get(std::string const &key) const
Definition: ParameterSet.h:231
int getPlaneID(const CellUniqueId &id) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Near Detector in the NuMI cavern.
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
const double j
Definition: BetheBloch.cxx:29
EventNumber_t event() const
Definition: Event.h:67
int Plane() const
Return plane value.
Definition: PCHit.h:24
geo::View_t View() const
Return view.
Definition: PCHit.h:36
Definition: View.py:1
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
void SetMC(bool isMC=true)
Definition: RawDigit.h:106
Histograms used by attenuation calibration.
const std::string path
Definition: plot_BEN.C:43
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
float Mag() const
TRandom3 r(0)
Timestamp time() const
Definition: Event.h:61
double GetPECorr(rb::CellHit const &cellhit, double w)
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
Float_t w
Definition: plot.C:20
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
enum BeamMode string