MuondEdxAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MuondEdxAna
3 // Module Type: filter
4 // File: MuondEdxAna_module.cc
5 //
6 // Generated at Tue May 7 16:23:55 2013 by Abigail Waldron using artmod
7 // from cetpkgsupport v1_02_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 // framework includes....
18 #include "fhiclcpp/ParameterSet.h"
19 #include "art_root_io/TFileService.h"
22 #include "Geometry/LiveGeometry.h"
24 
25 // novasoft stuff
27 #include "RecoBase/Track.h"
28 #include "RecoBase/CellHit.h"
29 //#include "MCCheater/BackTracker.h"
30 #include "Calibrator/Calibrator.h"
31 #include "NovaTimingUtilities/TimingUtilities.h"
32 
33 // C++
34 #include <memory>
35 #include <string>
36 #include <iostream>
37 #include <vector>
38 
39 // ROOT
40 #include "TH3D.h"
41 #include "TH2D.h"
42 #include "TProfile.h"
43 #include "TProfile2D.h"
44 #include "TFile.h"
45 #include "TVector3.h"
46 #include "TTree.h"
47 
48 namespace
49 {
50  struct MuondEdxAnaParams
51  {
52  template<class T> using Atom = fhicl::Atom<T>;
53  template<class T> using Table = fhicl::Table<T>;
54  using Comment = fhicl::Comment;
55  using Name = fhicl::Name;
56  Atom<std::string> PCHitsLabel{
57  Name("PCHitsLabel"),
58  Comment("Label of PCHitsList module")
59  };
60  Atom<std::string> stopperLabel{
61  Name("StopperLabel"),
62  Comment("Label of module which selected stopping tracks")
63  };
64  Atom<std::string> qualTreeHitsName{
65  Name("QualTreeHitsName"),
66  Comment("Fill tree with hits of this quality label, options listed in standard_pclist pset")
67  };
68  Atom<std::string> qualXYName{
69  Name("QualXYName"),
70  Comment("Instance label for XY-quality (Tricell) Hits")
71  };
72  Atom<std::string> qualZName{
73  Name("QualZName"),
74  Comment("Instance label for Z-quality (Triplane) Hits")
75  };
76  Atom<std::string> qualAvgName{
77  Name("QualAvgName"),
78  Comment("Instance label for Avg-quality (the rest...) Hits")
79  };
80  Atom<double> endPSystShift{
81  Name("EndPSystShift"),
82  Comment("cm to shift the track endpoint up/down")
83  };
84  };
85 }
86 
87 
88 namespace caldp {
89  class PCHit;
90 }
91 
92 namespace calib {
93  class MuondEdxAna;
94 
95  class MuondEdxAna : public art::EDFilter {
96  public:
98  explicit MuondEdxAna(const Parameters& params);
99  virtual ~MuondEdxAna();
100 
101  bool filter(art::Event & e) override;
102  void beginJob() override;
103  void endJob() override;
104  bool beginRun(art::Run & r) override;
105  bool endRun(art::Run & r) override;
106 
107  protected:
108  MuondEdxAnaParams fParams;
109 
110  double getPECorr(art::Event & e,caldp::PCHit const &pchit);
111  double getPECorrToGeV(art::Event & e, caldp::PCHit const &pchit);
112 
113  void fillRespVar(double PECorr, double PE, double trueE, double RecoGeV, double path, double w, int plane, int cell,
114  int view, double tns, double EventTimeHigh, int run, double dCosX, double dCosY, double dCosZ,
115  int diblock);
116  void fillLowLevel(double PECorr, double PE, double trueE, double RecoGeV, double path, double w, int plane, int cell,
117  int view, int diblock, double x, double trackLength, double FlightLength, double TotalTrackLength,
118  double dCosX, double dCosY, double dCosZ);
119 
122 
123  TH2D* fdEdx;
125  TH2D* fdEdx_minus;
126  TH2D* fdEdx_plus;
127  TH2D* fdEdx_true;
129 
132 
133  TH1D* fPECorr;
134  TH1D* fPE;
135  TH1D* fTrueE;
136  TH1D* fPECorr_X;
137  TH1D* fPE_X;
138  TH1D* fTrueE_X;
139  TH1D* fPECorr_Y;
140  TH1D* fPE_Y;
141  TH1D* fTrueE_Y;
142  TH1D* fRecoGeV;
143 
144  TH1D* fDiblock;
145  TProfile* fPECorrDiblock;
146  TProfile* fMevDiblock;
147 
148  TH1D* fPlane;
149  TH1D* fPlane_X;
150  TH1D* fPlane_Y;
151  TProfile* fPECorrPlane;
152  TProfile* fPECorrPlane_X;
153  TProfile* fPECorrPlane_Y;
154  TProfile* fMevPlane;
155  TProfile* fMevPlane_X;
156  TProfile* fMevPlane_Y;
157 
158  TH1D* fPathLength;
161  TProfile* fPECorrPath;
162  TProfile* fPECorrPath_X;
163  TProfile* fPECorrPath_Y;
164  TProfile* fMevPath;
165  TProfile* fMevPath_X;
166  TProfile* fMevPath_Y;
167 
168  TH1D* fW;
169  TH1D* fW_X;
170  TH1D* fW_Y;
171  TProfile* fPECorrW;
172  TProfile* fPECorrW_X;
173  TProfile* fPECorrW_Y;
174  TProfile* fMevW;
175  TProfile* fMevW_X;
176  TProfile* fMevW_Y;
177 
178  TH2D* fPathW_X;
179  TH2D* fPathW_Y;
180  TProfile2D* fPECorrPathW_X;
181  TProfile2D* fPECorrPathW_Y;
182  TProfile2D* fMeVPathW_X;
183  TProfile2D* fMeVPathW_Y;
184 
185  TH1D* fCell_X;
186  TH1D* fCell_Y;
187  TProfile* fPECorrCell_X;
188  TProfile* fPECorrCell_Y;
189  TProfile* fMevCell_X;
190  TProfile* fMevCell_Y;
191 
198  TH2D* fWvsCell_0;
199  TH2D* fWvsCell_1;
200 
201  TH2D* fdPEdx;
207 
208  TH1D* fStopPlane;
212 
215 
216  TH1D* fTNS;
217  TH1D* fRun;
218  TH1D* fEventTime;
219  TProfile* fPECorrTNS;
220  TProfile* fPECorrRun;
221  TProfile* fPECorrEventTime;
222 
223  // track level plots:
229 
234 
235  TH1D* fDCosX;
236  TH1D* fDCosY;
237  TH1D* fDCosZ;
238  TProfile* fRespDCosX;
239  TProfile* fRespDCosY;
240  TProfile* fRespDCosZ;
241 
242  TH1D* fEvtTime;
243  TProfile* fDetRespEvtTime;
244 
249 
250  TTree* fTree;
251 
252  int diblock;
253  int view;
256  float x;
257  float xCust;
258  float plane;
259  float path;
260  float truepath;
261  float trueW;
262  float trueE;
263  float PE;
264  float truePE;
265  float PECorr;
266  float tns;
267  uint32_t EventTimeHigh;
268  float cell;
269  float w;
270  float PECorrToGeV;
271  float RecoGeV;
273  int pdg;
274  int run;
275  int evt_year;
277  int evt_day;
278  int evt_hour;
279  int evt_min;
280  int evt_sec;
281  int APD;
284  int nhits;
285  int nXYHits;
286  bool badTNS;
287  float dCosX;
288  float dCosY;
289  float dCosZ;
291 
292 
293  // *time_t* of start of Nova epoch, 01-Jan-2010 00:00:00, UTC
294  // This is the value subtracted from the UNIX time_t, timeval or timespec
295  // seconds field.
296  // Using unsigned long long rather than uint64_t to ensure consistent behavior
297  // on OSX and LINUX
298  const uint32_t NOVA_EPOCH = 1262304000;
299  const unsigned long long NOVA_TIME_FACTOR = 64000000;
300 
301  };
302 
303 
304  MuondEdxAna::MuondEdxAna(const Parameters& params)
305  : EDFilter(params)
306  , fParams(params())
307  {
308  produces< TH2D, art::InRun >("dEdx");
309  produces< TH2D, art::InRun >("dEdxMinus");
310  produces< TH2D, art::InRun >("dEdxPlus");
311  produces< TH2D, art::InRun >("dEdxTrue");
312 
313  mf::LogInfo("MuondEdxAna") << "Filling Tree with hits labelled: "
314  << fParams.qualTreeHitsName();
315 
316 
317  }
318 
320  {
321  // Clean up dynamic memory and other resources here.
322  }
323 
325  {
326 
327  run = (int)(e.run());
328 
329  art::Timestamp ts = e.time();
330  EventTimeHigh = ts.timeHigh();
331  const time_t timeSec = ts.timeHigh();
332  struct tm* timeStruct = localtime(&timeSec);
333  evt_year=timeStruct->tm_year+1900;
334  evt_month=timeStruct->tm_mon+1;
335  evt_day=timeStruct->tm_mday;
336  evt_hour=timeStruct->tm_hour+1;
337  evt_min=timeStruct->tm_min;
338  evt_sec=timeStruct->tm_sec;
339 
341  //art::ServiceHandle<cheat::BackTracker> bt;
342 
343  /// Get the reconstructed tracks
345  e.getByLabel(fParams.stopperLabel(), tracks);
346  if( tracks->size() == 0 ) return false;
347 
348 
349  // Annoyingly need all since (XY + Z + Avg) = all hits in track
350  // PCHitTreeGetter gets the main hits to be put in tree
351  art::FindManyP<caldp::PCHit> PCHitTreeGetter(tracks, e,
352  art::InputTag(fParams.PCHitsLabel(),
353  fParams.qualTreeHitsName()));
354  art::FindManyP<caldp::PCHit> PCHitXYGetter( tracks, e,
355  art::InputTag(fParams.PCHitsLabel(),
356  fParams.qualXYName()));
357  art::FindManyP<caldp::PCHit> PCHitAvgGetter( tracks, e,
358  art::InputTag(fParams.PCHitsLabel(),
359  fParams.qualAvgName()));
360  art::FindManyP<caldp::PCHit> PCHitZGetter( tracks, e,
361  art::InputTag(fParams.PCHitsLabel(),
362  fParams.qualZName()));
363  if( !PCHitTreeGetter.isValid() ){
364  mf::LogInfo("MuondEdxAna") << "No hits for tree with label "
365  << fParams.qualTreeHitsName();
366  return false; }
367  if( !PCHitXYGetter.isValid() ){
368  mf::LogInfo("MuondEdxAna") << "No XY hits with label " << fParams.qualXYName();
369  return false; }
370  if( !PCHitAvgGetter.isValid() ){
371  mf::LogInfo("MuondEdxAna") << "No Avg hits with label " << fParams.qualAvgName();
372  return false; }
373  if( !PCHitZGetter.isValid() ){
374  mf::LogInfo("MuondEdxAna") << "No Z hits with label " << fParams.qualZName();
375  return false; }
376 
377 
378 
379  // loop over tracks in event
380  for( unsigned int i=0; i<tracks->size(); ++i ){
381 
382  //track level counters:
383  int windowTricells = 0;
384  int windowAvgHits = 0;
385  int windowZHits = 0;
386  int windowPlanes = 0;
387  std::vector<bool> trackWindowPlanes (1000, false); //vector to prevent double counting planes
388 
389  bool windowXView = false;
390  bool windowYView = false;
391 
392  //make custom track length to check against inbuilt length
393  TVector3 startXYZ = (*tracks)[i].Start();
394  TVector3 stopXYZ = (*tracks)[i].Stop();
395  double trackLength = (startXYZ-stopXYZ).Mag();
396 
397  // needed for track window, but good generally to reject protons
398  if( (*tracks)[i].TotalLength() < 200. ) continue;
399 
401  int stopPlaneID = 10000;
402  int firstMuCPlane = -1;
403  bool isNearDet = false;
404 
405  //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.
407  distanceToLiveDetEdgeZ = fLiveGeom->DistToEdgeZ((*tracks)[i].Stop());
408 
409  //continue if track ends in the muon catcher
411  isNearDet = true;
412  cid = fGeo->CellId(stopXYZ[0], stopXYZ[1], stopXYZ[2], 1.5, 1.5, 6., 1.0);
413  if(cid) stopPlaneID = fGeo->getPlaneID(cid);
414  firstMuCPlane = fGeo->FirstPlaneInMuonCatcher();
415  if(stopPlaneID >= firstMuCPlane) continue;
416  }
417 
419  fDistToFront_Stop -> Fill(livegeom->DistToFront(stopXYZ));
420  fDistToBack_Stop -> Fill(livegeom->DistToBack(stopXYZ));
421 
422  // xy and z conditions are mutually exclusive,
423  // avg hits are just the rest in the track
424  // so xy + z + avg is all the hits in the track
425  std::vector<art::Ptr<caldp::PCHit> > xyhits = PCHitXYGetter.at(i);
426  std::vector<art::Ptr<caldp::PCHit> > avghits = PCHitAvgGetter.at(i);
427  std::vector<art::Ptr<caldp::PCHit> > zhits = PCHitZGetter.at(i);
428  nhits = xyhits.size() + avghits.size() + zhits.size();
429 
430  // hits to loop through and fill tree
431  std::vector<art::Ptr<caldp::PCHit> > hits = PCHitTreeGetter.at(i);
432 
433  // ---------------- Fill vector containing hits per APD ---------------------- //
434  std::vector<int> APDhits (100,0);
435  std::vector< std::vector<int> > DiblockAPDhits (100, std::vector<int> (100,0));
436  badTNS = false; //temporary handle for showery events in old data
437 
438  // tricell hits
439  for( unsigned int n=0; n<xyhits.size(); ++n ){
440  APDhits[xyhits[n]->APD()] += 1;
441  DiblockAPDhits.at(xyhits[n]->Diblock()).at(xyhits[n]->APD()) += 1;
442  if(xyhits[n]->TNS() < 0 || xyhits[n]->TNS() > 1400) badTNS = true;
443  }
444  //avg. hits
445  for( unsigned int p=0; p<avghits.size(); ++p ) {
446  APDhits[avghits[p]->APD()] += 1;
447  DiblockAPDhits.at(avghits[p]->Diblock()).at(avghits[p]->APD()) += 1;
448  if(avghits[p]->TNS() < 0 || avghits[p]->TNS() > 1400) badTNS = true;
449  FlightLength = avghits[p]->FlightLen();
450  TotalTrackLength = (*tracks)[i].TotalLength();
451  x = trackLength - FlightLength;
452  if(x<100 || x>200) continue;
453  windowAvgHits ++; //track level counters
454  }
455  //z hits
456  for( unsigned int q=0; q<zhits.size(); ++q ){
457  APDhits[zhits[q]->APD()] += 1;
458  DiblockAPDhits.at(zhits[q]->Diblock()).at(zhits[q]->APD()) += 1;
459  if(zhits[q]->TNS() < 0 || zhits[q]->TNS() > 1400) badTNS = true;
460  FlightLength = zhits[q]->FlightLen();
461  TotalTrackLength = (*tracks)[i].TotalLength();
462  x = trackLength - FlightLength;
463  if(x<100 || x>200) continue;
464  windowZHits ++; //track level counters
465  }
466  // ---------------- Filled APD vector -------------------------------------------- //
467 
468  //cosine giving initial direction of track
469  dCosX = (*tracks)[i].Dir().X()/(*tracks)[i].Dir().Mag();
470  dCosY = (*tracks)[i].Dir().Y()/(*tracks)[i].Dir().Mag();
471  dCosZ = (*tracks)[i].Dir().Z()/(*tracks)[i].Dir().Mag();
472 
473  for( unsigned int j=0; j<hits.size(); ++j ){
474 
475  if(isNearDet)
476  if(hits[j]->Plane() >= firstMuCPlane) continue;
477 
478  /*
479  int p = hits[j]->Plane();
480  int c = hits[j]->Cell();
481  rb::CellHit chit;
482  chit.SetPlane( p );
483  chit.SetCell ( c );
484  //double tns_calib = hits[j]->TNS() + fCalib->GetTimingOffset(p,c,e.isRealData());
485  double tns_calib = hits[j]->TNS();
486  chit.SetTNS(tns_calib,hits[j]->GoodTime()); // may need to get offset from calibrator...
487  if(!e.isRealData())chit.SetMC(); // assume only filled for MC...
488  const std::vector<sim::FLSHit> flss = bt->HitToFLSHit(chit);
489  int lastpdg=0;
490  bool uniformPdg=true;
491  for(const sim::FLSHit& fls: flss){
492  if(lastpdg!=fls.GetPDG() && lastpdg!=0){
493  std::cout << "multiple pdg in flshit: "
494  << lastpdg << " vs. " << fls.GetPDG() << std::endl;
495  uniformPdg = false;
496  }
497  lastpdg=fls.GetPDG();
498  }
499  if(flss.size()>0 && uniformPdg)
500  pdg = flss[0].GetPDG();
501  else pdg = -1;
502  */
503 
504  //test poisson lambda:
505  //double poissonLambda = hits[j]->PoissonLambda();
506  //std::cout<< "poissonLambda: " << poissonLambda <<std::endl;
507 
508  nXYHits = hits.size();
509  diblock = hits[j]->Diblock();
510  view = (int)hits[j]->View();
511  FlightLength = hits[j]->FlightLen();
512  TotalTrackLength= (*tracks)[i].TotalLength();
513  x = trackLength - FlightLength;
514  plane = hits[j]->Plane();
515  path = hits[j]->Path();
516  truepath = hits[j]->TruePath();
517  trueW = hits[j]->TrueW();
518  trueE = hits[j]->TrueMeV();
519  PE = hits[j]->PE();
520  truePE = hits[j]->TruePE();
521  PECorr = this->getPECorr(e,*hits[j]);
522  PECorrToGeV = this->getPECorrToGeV(e,*hits[j]);
524  tns = hits[j]->TNS();
525  cell = hits[j]->Cell();
526  w = hits[j]->W();
527  APD = hits[j]->APD();
528  sameAPDhitsEvent= APDhits[APD];
529  sameAPDhitDiblockEvent = DiblockAPDhits[diblock][APD];
531 
532  fTree->Fill(); //fill ttree
533 
534  fdPEdx ->Fill(x,PE/path);
535  fdEdx ->Fill(x,PECorr/path);
536  fdEdx_minus->Fill(x-fParams.endPSystShift(),PECorr/path);
537  fdEdx_plus ->Fill(x+fParams.endPSystShift(),PECorr/path);
538  fdEdx_true ->Fill(x,trueE/path);
539 
540  if(w > 0){
543  }
544 
545  if(x<100 || x>200) continue;
546 
547  //track level counters
548  windowTricells ++;
549  if(view == 0) windowXView = true;
550  if(view == 1) windowYView = true;
551  //check if plane already counted
552  if(trackWindowPlanes[plane] == false) windowPlanes ++;
553  trackWindowPlanes[plane] = true;
554 
556  fillLowLevel(PECorr, PE, trueE, RecoGeV, path, w, plane, cell, view, diblock, x, trackLength, FlightLength, TotalTrackLength, dCosX, dCosY, dCosZ);
557 
558  } // loop through hits
559 
560  fWindowTricellsPerTrack -> Fill(windowTricells);
561  fWindowPlanesPerTrack -> Fill(windowPlanes);
562  fWindowHitsPerTrack -> Fill(windowTricells + windowAvgHits + windowZHits);
563  fWinTricellsHitsPerTrack -> Fill(windowTricells,windowTricells + windowAvgHits + windowZHits);
564  fWinTricellsPlanesPerTrack -> Fill(windowTricells,windowPlanes);
565  fWinHitsPlanesPerTrack -> Fill(windowTricells + windowAvgHits + windowZHits, windowPlanes);
566 
567  int nViewsInWindow = 0;
568  if(windowXView == true) nViewsInWindow++;
569  if(windowYView == true) nViewsInWindow++;
570 
571  fWinTricellsViewsPerTrack -> Fill(windowTricells,nViewsInWindow);
572  fWindowViewsPerTrack -> Fill(nViewsInWindow);
573  if(nViewsInWindow == 1){
574  if(windowXView) fWindowOneViewTracks -> Fill(0);
575  if(windowYView) fWindowOneViewTracks -> Fill(1);
576  }
577  }
578 
579  return true;
580  }
581 
582  void MuondEdxAna::fillLowLevel(double PECorr, double PE, double trueE, double RecoGeV, double path, double w, int plane, int cell,
583  int view, int diblock, double x, double trackLength, double FlightLength, double TotalTrackLength,
584  double dCosX, double dCosY, double dCosZ){
585 
586  fDCosX -> Fill(dCosX);
587  fDCosY -> Fill(dCosY);
588  fDCosZ -> Fill(dCosZ);
589  fDiblock -> Fill(diblock);
590  fPECorr -> Fill(PECorr/path);
591  fPE -> Fill(PE/path);
592  fTrueE -> Fill(trueE/path);
593  fRecoGeV -> Fill(1000 * RecoGeV/path);
594 
595 
596  /*
597  if( (
598  fGeo->DetId() == novadaq::cnv::kNEARDET &&
599  (
600  (view==0 && fMinFlatW_ND_x < w && w < fMaxFlatW_ND_x ) ||
601  (view==1 && fMinFlatW_ND_y < w && w < fMaxFlatW_ND_y )
602  )
603  ) || (
604  fGeo->DetId() == novadaq::cnv::kFARDET &&
605  (
606  (view==0 && fMinFlatW_FD_x < w && w < fMaxFlatW_FD_x ) ||
607  (view==1 && fMinFlatW_FD_y < w && w < fMaxFlatW_FD_y )
608  )
609  )
610  ){
611  fMEUreco_flat_w->Fill(PECorr/path);
612  fMEUtrue_flat_w->Fill(trueE/path);
613  }
614  */
615 
616  if(view == 0){
617  fPECorr_X -> Fill(PECorr/path);
618  fPE_X -> Fill(PE/path);
619  fTrueE_X -> Fill(trueE/path);
620  }
621 
622  if(view == 1){
623  fPECorr_Y -> Fill(PECorr/path);
624  fPE_Y -> Fill(PE/path);
625  fTrueE_Y -> Fill(trueE/path);
626  }
627 
628  fTrackEndDist -> Fill(x);
629  fTrackEndDistCustom-> Fill(trackLength - FlightLength);
630  fFlightLength -> Fill(FlightLength);
631  fTrackLength -> Fill(TotalTrackLength);
632  fTrackLengthCustom -> Fill(trackLength);
633 
634  }
635 
636  void MuondEdxAna::fillRespVar(double PECorr, double PE, double trueE, double RecoGeV, double path, double w,
637  int plane, int cell, int view, double tns, double EventTimeHigh, int run, double dCosX,
638  double dCosY, double dCosZ, int diblock){
639 
640  fPathLength -> Fill(path);
641  fPlane -> Fill(plane);
642  fW -> Fill(w);
643  fPECorrPlane -> Fill(plane, PECorr/path);
644  fMevPlane -> Fill(plane, trueE/path);
645  fPECorrW -> Fill(w, PECorr/path);
646  fMevW -> Fill(w, trueE/path);
647  fPECorrPath -> Fill(path,PECorr/path);
648  fMevPath -> Fill(path,trueE/path);
649  fPECorrDiblock -> Fill(diblock, PECorr/path);
650  fMevDiblock -> Fill(diblock, trueE/path);
651  fTNS -> Fill(tns);
652  fPECorrTNS -> Fill(tns,PECorr/path);
653  fRun -> Fill(run);
654  fPECorrRun -> Fill(run,PECorr/path);
655  fEventTime -> Fill(EventTimeHigh);
656  fPECorrEventTime-> Fill(EventTimeHigh,PECorr/path);
657  fRespDCosX -> Fill(dCosX, PECorr/path);
658  fRespDCosY -> Fill(dCosY, PECorr/path);
659  fRespDCosZ -> Fill(dCosZ, PECorr/path);
660 
661  if(view == 0){
662  fPathW_X -> Fill(w, path);
663  fPECorrPathW_X -> Fill(w, path,PECorr/path);
664  fMeVPathW_X -> Fill(w, path,trueE/path);
665  fPathLength_X -> Fill(path);
666  fPECorrPath_X -> Fill(path,PECorr/path);
667  fMevPath_X -> Fill(path,trueE/path);
668  fPlane_X -> Fill(plane);
669  fW_X -> Fill(w);
670  fPECorrPlane_X -> Fill(plane, PECorr/path);
671  fMevPlane_X -> Fill(plane, trueE/path);
672  fPECorrW_X -> Fill(w, PECorr/path);
673  fMevW_X -> Fill(w, trueE/path);
674  fCell_X -> Fill(cell);
675  fPECorrCell_X -> Fill(cell, PECorr/path);
676  fMevCell_X -> Fill(cell, trueE/path);
677  fWvsCell_PECorr_0 -> Fill(cell, w, PECorr/path);
678  fWvsCell_MEV_0 -> Fill(cell, w, trueE/path);
679  fWvsCell_0 -> Fill(cell, w);
680 
681  }
682 
683  if(view == 1){
684  fPathW_Y -> Fill(w, path);
685  fPECorrPathW_Y -> Fill(w, path,PECorr/path);
686  fMeVPathW_Y -> Fill(w, path,trueE/path);
687  fPathLength_Y -> Fill(path);
688  fPECorrPath_Y -> Fill(path,PECorr/path);
689  fMevPath_Y -> Fill(path,trueE/path);
690  fPlane_Y -> Fill(plane);
691  fW_Y -> Fill(w);
692  fPECorrPlane_Y -> Fill(plane, PECorr/path);
693  fMevPlane_Y -> Fill(plane, trueE/path);
694  fPECorrW_Y -> Fill(w, PECorr/path);
695  fMevW_Y -> Fill(w, trueE/path);
696  fCell_Y -> Fill(cell);
697  fPECorrCell_Y -> Fill(cell, PECorr/path);
698  fMevCell_Y -> Fill(cell, trueE/path);
699  fWvsCell_PECorr_1-> Fill(cell, w, PECorr/path);
700  fWvsCell_MEV_1 -> Fill(cell, w, trueE/path);
701  fWvsCell_1 -> Fill(cell, w);
702 
703  }
704 
705  }
706 
708 
709  //Make mock cellhit
710  rb::CellHit cellhit;
711  cellhit.SetPlane(pchit.Plane());
712  cellhit.SetCell(pchit.Cell());
713  cellhit.SetView(pchit.View());
714  cellhit.SetTNS(pchit.TNS(),pchit.GoodTime());
715  if(!e.isRealData())cellhit.SetMC(); // assume only filled for MC...
716  cellhit.SetPE(pchit.PE());
717 
718  double pecorr = fCalib->GetPECorr(cellhit, pchit.W());
719  return pecorr;
720  }
721 
723 
724  rb::CellHit cellhit;
725  if(!e.isRealData()) cellhit.SetMC(); // assume only filled for MC...
726 
727  double pecorrtogev = fCalib->GetPECorrToGeVScale(cellhit);
728  return pecorrtogev;
729  }
730 
731 
733  {
734  // histograms for output
736 
737  fTree = tfs->make<TTree>("fTree","tree to hold muondedx variables");
738  fTree->Branch("nXYHits",&nXYHits);
739  fTree->Branch("nhits",&nhits);
740  fTree->Branch("dCosX",&dCosX);
741  fTree->Branch("dCosY",&dCosY);
742  fTree->Branch("dCosZ",&dCosZ);
743  fTree->Branch("x",&x);
744  fTree->Branch("badTNS",&badTNS);
745  fTree->Branch("run",&run);
746  fTree->Branch("diblock",&diblock);
747  fTree->Branch("view",&view);
748  fTree->Branch("FlightLength",&FlightLength);
749  fTree->Branch("TotalTrackLength",&TotalTrackLength);
750  fTree->Branch("xCust",&xCust);
751  fTree->Branch("distanceToLiveDetEdgeZ",&distanceToLiveDetEdgeZ);
752  fTree->Branch("plane",&plane);
753  fTree->Branch("path",&path);
754  fTree->Branch("truepath",&truepath);
755  fTree->Branch("trueW",&trueW);
756  fTree->Branch("trueE",&trueE);
757  fTree->Branch("truePE",&truePE);
758  fTree->Branch("PE",&PE);
759  fTree->Branch("PECorr",&PECorr);
760  fTree->Branch("tns",&tns);
761  fTree->Branch("cell",&cell);
762  fTree->Branch("w",&w);
763  fTree->Branch("pdg",&pdg);
764  fTree->Branch("PECorrToGeV",&PECorrToGeV);
765  fTree->Branch("RecoGeV",&RecoGeV);
766  fTree->Branch("EventTimeHigh",&EventTimeHigh);
767  fTree->Branch("evt_year",&evt_year);
768  fTree->Branch("evt_month",&evt_month);
769  fTree->Branch("evt_day",&evt_day);
770  fTree->Branch("evt_hour",&evt_hour);
771  fTree->Branch("evt_min",&evt_min);
772  fTree->Branch("evt_sec",&evt_sec);
773  fTree->Branch("APD",&APD);
774  fTree->Branch("sameAPDhitsEvent",&sameAPDhitsEvent);
775  fTree->Branch("sameAPDhitDiblockEvent",&sameAPDhitDiblockEvent);
776  fTree->Branch("fibBrightnessBin",&fibBrightnessBin);
777 
778  fDeltaStartPlane = tfs->make<TH1D>("DeltaStartPlane",";startPlane diff. (x,y);",100,0,100);
779  fDeltaEndPlane = tfs->make<TH1D>("DeltaEndPlane",";endPlane diff. (x,y);",100,0,100);
780  fPlaneAsymmetry = tfs->make<TH1D>("PlaneAsymmetry",";plane asymm. (x,y);",100,0,100);
781  fPlaneAsymmetry2D = tfs->make<TH2D>("PlaneAsymmetry2D",";diff in view-planes(x-y); sum of view-planes (x+y);",100,0,100,1000,0,1000);
782 
783  fEvtTime = tfs->make<TH1D>("EvtTime",";EvtTime (YYYYMMDD);",std::pow(10,6),0,std::pow(10,6));
784  fDetRespEvtTime = tfs->make<TProfile>("DetRespEvtTime",";EvtTime (YYYYMMDD); PECorr/cm",std::pow(10,6),0,std::pow(10,6));
785 
786  fDCosX = tfs->make<TH1D>("DCosX",";dCosX;",200,-1,1);
787  fDCosY = tfs->make<TH1D>("DCosY",";dCosY;",200,-1,1);
788  fDCosZ = tfs->make<TH1D>("DCosZ",";dCosZ;",200,-1,1);
789  fRespDCosX = tfs->make<TProfile>("RespDCosX",";dCosX;",200,-1,1);
790  fRespDCosY = tfs->make<TProfile>("RespDCosY",";dCosY;",200,-1,1);
791  fRespDCosZ = tfs->make<TProfile>("RespDCosZ",";dCosZ;",200,-1,1);
792 
793  fdPEdx = tfs->make<TH2D>("dPEdx",";Distance to track end (cm);PE / cm;",
794  50,0.,500.,200,0.,200.);
795  fdEdx = tfs->make<TH2D>("dEdx",";Distance to track end (cm);PECorr / cm;",
796  50,0.,500.,200,0.,200.);
797  fdEdx_minus = tfs->make<TH2D>("dEdx_minus",";Distance to track end (cm);PECorr / cm;",
798  50,0.,500.,200,0.,200.);
799  fdEdx_plus = tfs->make<TH2D>("dEdx_plus",";Distance to track end (cm);PECorr / cm;",
800  50,0.,500.,200,0.,200.);
801  fdEdx_true = tfs->make<TH2D>("dEdx_true",";Distance to track end (cm);True MeV / cm;",
802  50,0.,500.,200,0.,5.);
803 
804  fdEdx_MuonCatch = tfs->make<TH2D>("dEdx_MuonCatch",";Distance to track end (cm);PECorr / cm;",
805  50,0.,500.,200,0.,200.);
806  fdEdx_true_MuonCatch = tfs->make<TH2D>("dEdx_true_MuonCatch",";Distance to track end (cm);True MeV / cm;",
807  50,0.,500.,200,0.,5.);
808 
809  fdEdx_NearReadout = tfs->make<TH2D>("dEdx_NearReadout",";Distance to track end (cm);PECorr / cm;",
810  50,0.,500.,200,0.,200.);
811  fdEdx_true_NearReadout = tfs->make<TH2D>("dEdx_true_NearReadout",";Distance to track end (cm);True MeV / cm;",
812  50,0.,500.,200,0.,5.);
813 
814 
815  //fMEUreco_flat_w = tfs->make<TH1D>("meureco_flat_w",";PECorr/cm;",1000,0,500);
816  //fMEUtrue_flat_w = tfs->make<TH1D>("meutrue_flat_w",";MeV/cm;",1000,0,25);
817 
818  fPECorr = tfs->make<TH1D>("PECorr",";PECorr/cm;",1000,0,500);
819  fPE = tfs->make<TH1D>("PE",";PE/cm;",1000,0,500);
820  fTrueE = tfs->make<TH1D>("TrueE",";MeV/cm;",1000,0,25);
821  fRecoGeV = tfs->make<TH1D>("RecoGeV",";MeV/cm;",1000,0,25);
822  fPECorr_X = tfs->make<TH1D>("PECorr_X",";PECorr/cm;",1000,0,500);
823  fPE_X = tfs->make<TH1D>("PE_X",";PE/cm;",1000,0,500);
824  fTrueE_X = tfs->make<TH1D>("TrueE_X",";MeV/cm;",1000,0,25);
825  fPECorr_Y = tfs->make<TH1D>("PECorr_Y",";PECorr/cm;",1000,0,500);
826  fPE_Y = tfs->make<TH1D>("PE_Y",";PE/cm;",1000,0,500);
827  fTrueE_Y = tfs->make<TH1D>("TrueE_Y",";MeV/cm;",1000,0,25);
828 
829  fFlightLength = tfs->make<TH1D>("FlightLength",";distance from track start (cm);",1000,0,10000);
830  fTrackLength = tfs->make<TH1D>("TrackLength",";Total track length (cm);",1000,0,10000);
831  fTrackEndDist = tfs->make<TH1D>("TrackEndDist",";Distance to end of track (cm)",1000,0,10000);
832  fTrackLengthCustom = tfs->make<TH1D>("TrackLengthCustom",";Total track length (cm);",1000,0,10000);
833  fTrackEndDistCustom = tfs->make<TH1D>("TrackEndDistCustom",";Distance to end of track (cm)",1000,0,10000);
834  fTrackEndPlane = tfs->make<TH1D>("TrackEndPlane",";plane of track end;",1000,0,1000);
835  fStopPlane = tfs->make<TH1D>("StopPlane",";start and stop z positions;",1000,0,1000);
836  fStopPlane_liveGeomCut = tfs->make<TH1D>("StopPlane_liveGeomCut",";start and stop z positions;",1000,0,1000);
837  fDistToFront_Stop = tfs->make<TH1D>("DistToFront_Stop",";distance to front (cm);",1000,0,1000);
838  fDistToBack_Stop = tfs->make<TH1D>("DistToBack_Stop",";distance to back (cm);",1000,0,1000);
839 
840  fTNS = tfs->make<TH1D>("TNS",";TNS (nova time);",std::pow(10,3),-5*std::pow(10,4),5*std::pow(10,4));
841  fPECorrTNS = tfs->make<TProfile>("PECorrTNS",";TNS (nova time); PECorr/cm;",std::pow(10,3),-5*std::pow(10,4),5*std::pow(10,4));
842  fRun = tfs->make<TH1D>("Run",";Run;",10000,0,100000);
843  fPECorrRun = tfs->make<TProfile>("PECorrRun",";Run; PECorr/cm;",10000,0,100000);
844 
845  fEventTime = tfs->make<TH1D>("EventTime",";date and time (seconds);",std::pow(10,6),0,std::pow(10,10));
846  fPECorrEventTime = tfs->make<TProfile>("PECorrEventTime",";date and time (seconds); PECorr/cm;",std::pow(10,6),0,std::pow(10,10));
847 
848  fPathW_X = tfs->make<TH2D>("PathW_X",";W (cm); Path (cm)",160,-800,800, 100, 2, 12);
849  fPathW_Y = tfs->make<TH2D>("PathW_Y",";W (cm); Path (cm)",160,-800,800, 100, 2, 12);
850  fPECorrPathW_X = tfs->make<TProfile2D>("PECorrPathW_X",";W (cm); Path (cm); PECorr/cm",160,-800,800, 100, 2, 12);
851  fPECorrPathW_Y = tfs->make<TProfile2D>("PECorrPathW_Y",";W (cm); Path (cm); PECorr/cm",160,-800,800, 100, 2, 12);
852  fMeVPathW_X = tfs->make<TProfile2D>("MeVPathW_X",";W (cm); Path (cm); MeV/cm",160,-800,800, 100, 2, 12);
853  fMeVPathW_Y = tfs->make<TProfile2D>("MeVPathW_Y",";W (cm); Path (cm); MeV/cm",160,-800,800, 100, 2, 12);
854 
855  fPathLength = tfs->make<TH1D>("PathLength",";path length (cm);",200,0,20);
856  fPathLength_X = tfs->make<TH1D>("PathLength_X",";path length (cm);",200,0,20);
857  fPathLength_Y = tfs->make<TH1D>("PathLength_Y",";path length (cm);",200,0,20);
858  fPECorrPath = tfs->make<TProfile>("PECorrPath",";path length X-view (cm); PECorr/cm;",200,0,20);
859  fPECorrPath_X = tfs->make<TProfile>("PECorrPath_X",";path length X-view (cm); PECorr/cm;",200,0,20);
860  fPECorrPath_Y = tfs->make<TProfile>("PECorrPath_Y",";path length X-view (cm); PECorr/cm;",200,0,20);
861  fMevPath = tfs->make<TProfile>("MevPath",";path length X-view (cm); Mev/cm;",200,0,20);
862  fMevPath_X = tfs->make<TProfile>("MevPath_X",";path length X-view (cm); Mev/cm;",200,0,20);
863  fMevPath_Y = tfs->make<TProfile>("MevPath_Y",";path length X-view (cm); Mev/cm;",200,0,20);
864 
865  fPlane = tfs->make<TH1D>("Plane",";plane;",1000,0,1000);
866  fPlane_X = tfs->make<TH1D>("Plane_X",";plane;",1000,0,1000);
867  fPlane_Y = tfs->make<TH1D>("Plane_Y",";plane;",1000,0,1000);
868  fPECorrPlane = tfs->make<TProfile>("PECorrPlane",";plane; PECorr;",1000,0,1000);
869  fPECorrPlane_X = tfs->make<TProfile>("PECorrPlane_X",";plane; PECorr;",1000,0,1000);
870  fPECorrPlane_Y = tfs->make<TProfile>("PECorrPlane_Y",";plane; PECorr;",1000,0,1000);
871  fMevPlane = tfs->make<TProfile>("MevPlane",";plane; Mev;",1000,0,1000);
872  fMevPlane_X = tfs->make<TProfile>("MevPlane_X",";plane; Mev;",1000,0,1000);
873  fMevPlane_Y = tfs->make<TProfile>("MevPlane_Y",";plane; Mev;",1000,0,1000);
874 
875  fW = tfs->make<TH1D>("W",";w;",160,-800,800);
876  fW_X = tfs->make<TH1D>("W_X",";w;",160,-800,800);
877  fW_Y = tfs->make<TH1D>("W_Y",";w;",160,-800,800);
878  fPECorrW = tfs->make<TProfile>("PECorrW",";w; PECorr/cm;",160,-800,800);
879  fPECorrW_X = tfs->make<TProfile>("PECorrW_X",";w; PECorr/cm;",160,-800,800);
880  fPECorrW_Y = tfs->make<TProfile>("PECorrW_Y",";w; PECorr/cm;",160,-800,800);
881  fMevW = tfs->make<TProfile>("MevW",";w; Mev/cm;",160,-800,800);
882  fMevW_X = tfs->make<TProfile>("MevW_X",";w; Mev/cm;",160,-800,800);
883  fMevW_Y = tfs->make<TProfile>("MevW_Y",";w; Mev/cm;",160,-800,800);
884 
885  fCell_X = tfs->make<TH1D>("Cell_X",";cell;",1000,0,1000);
886  fCell_Y = tfs->make<TH1D>("Cell_Y",";cell;",1000,0,1000);
887  fPECorrCell_X = tfs->make<TProfile>("PECorrCell_X",";cell; PECorr/cm;",1000,0,1000);
888  fPECorrCell_Y = tfs->make<TProfile>("PECorrCell_Y",";cell; PECorr/cm;",1000,0,1000);
889  fMevCell_X = tfs->make<TProfile>("MevCell_X",";cell; Mev/cm;",1000,0,1000);
890  fMevCell_Y = tfs->make<TProfile>("MevCell_Y",";cell; Mev/cm;",1000,0,1000);
891 
892  //track level histograms
893  fWindowTricellsPerTrack = tfs->make<TH1D>("WindowTricellsPerTrack",";tricells per track in track window;",50,0,50);
894  fWindowHitsPerTrack = tfs->make<TH1D>("WindowHitsPerTrack",";hits per track in track window;",1000,0,1000);
895  fWindowPlanesPerTrack = tfs->make<TH1D>("WindowPlanesPerTrack",";planes per track in track window;",1000,0,1000);
896  fWindowViewsPerTrack = tfs->make<TH1D>("WindowViewsPerTrack",";views per track in track window;",10,0,10);
897  fWindowOneViewTracks = tfs->make<TH1D>("WindowOneViewTracks",";View in window when only one present;",10,0,10);
898 
899  fWinTricellsHitsPerTrack = tfs->make<TH2D>("WinTricellsHitsPerTrack",";tricells per track in track window;hits per track in window;",
900  50,0,50,60,0,60);
901  fWinTricellsPlanesPerTrack = tfs->make<TH2D>("WinTricellsPlanesPerTrack",";tricells per track in track window;planes hit per track in window;",
902  50,0,50,10,0,10);
903  fWinHitsPlanesPerTrack = tfs->make<TH2D>("WinHitsPlanesPerTrack",";hits per track in window;planes per track in window",
904  60,0,60,10,0,10);
905  fWinTricellsViewsPerTrack = tfs->make<TH2D>("WinTricellsViewsPerTrack",";tricells per track in track window;views per track in window;",
906  50,0,50,4,0,4);
907 
908  fWvsCell_0 = tfs->make<TH2D>("WvsCell_0",";Cell;W;",
909  40,0.,400.,160,-800.,800.);
910  fWvsCell_1 = tfs->make<TH2D>("WvsCell_1",";Cell;W;",
911  40,0.,400.,160,-800.,800.);
912  fWvsCell_PECorr_0 = tfs->make<TH3D>("WvsCell_PECorr_0",";Cell;W;",
913  40,0.,400.,160,-800.,800.,100,0.,100.);
914  fWvsCell_PECorr_1 = tfs->make<TH3D>("WvsCell_PECorr_1",";Cell;W;",
915  40,0.,400.,160,-800.,800.,100,0.,100.);
916  fWvsCell_MEV_0 = tfs->make<TH3D>("WvsCell_MEV_0",";Cell;W;",
917  40,0.,400.,160,-800.,800.,100,0.,10.);
918  fWvsCell_MEV_1 = tfs->make<TH3D>("WvsCell_MEV_1",";Cell;W;",
919  40,0.,400.,160,-800.,800.,100,0.,10.);
920 
921  fDiblock = tfs->make<TH1D>("Diblock","diblock",20,0,20);
922  fPECorrDiblock = tfs->make<TProfile>("PECorrDiblock",";diblock; PECorr/cm",20,0,20);
923  fMevDiblock = tfs->make<TProfile>("MevDiblock",";diblock; MeV/cm",20,0,20);
924 
925  ftrackEndPlaneWithinDiblock = tfs->make<TH1D>("trackEndPlaneWithinDiblock", ";plane of track end within diblocks;",100,0,100);
926 
927  }
928 
930  return true;
931  }
932 
934 
936  std::unique_ptr<TH2D> pdEdx(tfs->make<TH2D>(*fdEdx));
937  std::unique_ptr<TH2D> pdEdx_minus(tfs->make<TH2D>(*fdEdx_minus));
938  std::unique_ptr<TH2D> pdEdx_plus(tfs->make<TH2D>(*fdEdx_plus));
939  std::unique_ptr<TH2D> pdEdx_true(tfs->make<TH2D>(*fdEdx_true));
940 
941  r.put(std::move(pdEdx),"dEdx");
942  r.put(std::move(pdEdx_minus),"dEdxMinus");
943  r.put(std::move(pdEdx_plus),"dEdxPlus");
944  r.put(std::move(pdEdx_true),"dEdxTrue");
945  return true;
946  }
947 
949  {
950 
951  }
952 }// end namespace
float TNS() const
Return uncorrected hit time.
Definition: PCHit.h:54
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
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const uint64_t NOVA_TIME_FACTOR
void fillLowLevel(double PECorr, double PE, double trueE, double RecoGeV, double path, double w, int plane, int cell, int view, int diblock, double x, double trackLength, double FlightLength, double TotalTrackLength, double dCosX, double dCosY, double dCosZ)
const char * p
Definition: xmltok.h:285
bool GoodTime() const
Return quality of timing fit for cell.
Definition: PCHit.h:56
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
constexpr T pow(T x)
Definition: pow.h:72
const uint32_t NOVA_EPOCH
int Cell() const
Return cell value.
Definition: PCHit.h:26
double DistToFront(TVector3 vertex)
double DistToBack(TVector3 vertex)
DEFINE_ART_MODULE(TestTMapFile)
bool endRun(art::Run &r) override
"Pre-calibration hit". Common input to calibration procedures
Definition: PCHit.h:16
Definition: Run.h:21
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:44
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
void SetCell(unsigned short cell)
Definition: CellHit.h:52
bool isRealData() const
art::ServiceHandle< calib::Calibrator > fCalib
CDPStorage service.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void hits()
Definition: readHits.C:15
void SetPE(float pe)
Definition: CellHit.h:55
Timestamp time() const
float PE() const
Return PE value.
Definition: PCHit.h:38
double DistToEdgeZ(TVector3 vertex)
void beginJob()
double getPECorrToGeV(art::Event &e, caldp::PCHit const &pchit)
Encapsulate the geometry of one entire detector (near, far, ndos)
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
std::void_t< T > n
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
RunNumber_t run() const
int Plane() const
Return plane value.
Definition: PCHit.h:24
geo::View_t View() const
Return view.
Definition: PCHit.h:36
void fillRespVar(double PECorr, double PE, double trueE, double RecoGeV, double path, double w, int plane, int cell, int view, double tns, double EventTimeHigh, int run, double dCosX, double dCosY, double dCosZ, int diblock)
Definition: run.py:1
void SetMC(bool isMC=true)
Definition: RawDigit.h:106
bool filter(art::Event &e) override
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
void beginJob() override
float Mag() const
TRandom3 r(0)
void endJob() override
double GetPECorr(rb::CellHit const &cellhit, double w)
art::ServiceHandle< geo::Geometry > fGeo
Float_t e
Definition: plot.C:35
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
double getPECorr(art::Event &e, caldp::PCHit const &pchit)
Float_t w
Definition: plot.C:20
bool beginRun(art::Run &r) override
MuondEdxAnaParams fParams
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.