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"
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  : fParams(params())
306  {
307  produces< TH2D, art::InRun >("dEdx");
308  produces< TH2D, art::InRun >("dEdxMinus");
309  produces< TH2D, art::InRun >("dEdxPlus");
310  produces< TH2D, art::InRun >("dEdxTrue");
311 
312  mf::LogInfo("MuondEdxAna") << "Filling Tree with hits labelled: "
313  << fParams.qualTreeHitsName();
314 
315 
316  }
317 
319  {
320  // Clean up dynamic memory and other resources here.
321  }
322 
324  {
325 
326  run = (int)(e.run());
327 
328  art::Timestamp ts = e.time();
329  EventTimeHigh = ts.timeHigh();
330  const time_t timeSec = ts.timeHigh();
331  struct tm* timeStruct = localtime(&timeSec);
332  evt_year=timeStruct->tm_year+1900;
333  evt_month=timeStruct->tm_mon+1;
334  evt_day=timeStruct->tm_mday;
335  evt_hour=timeStruct->tm_hour+1;
336  evt_min=timeStruct->tm_min;
337  evt_sec=timeStruct->tm_sec;
338 
340  //art::ServiceHandle<cheat::BackTracker> bt;
341 
342  /// Get the reconstructed tracks
344  e.getByLabel(fParams.stopperLabel(), tracks);
345  if( tracks->size() == 0 ) return false;
346 
347 
348  // Annoyingly need all since (XY + Z + Avg) = all hits in track
349  // PCHitTreeGetter gets the main hits to be put in tree
350  art::FindManyP<caldp::PCHit> PCHitTreeGetter(tracks, e,
351  art::InputTag(fParams.PCHitsLabel(),
352  fParams.qualTreeHitsName()));
353  art::FindManyP<caldp::PCHit> PCHitXYGetter( tracks, e,
354  art::InputTag(fParams.PCHitsLabel(),
355  fParams.qualXYName()));
356  art::FindManyP<caldp::PCHit> PCHitAvgGetter( tracks, e,
357  art::InputTag(fParams.PCHitsLabel(),
358  fParams.qualAvgName()));
359  art::FindManyP<caldp::PCHit> PCHitZGetter( tracks, e,
360  art::InputTag(fParams.PCHitsLabel(),
361  fParams.qualZName()));
362  if( !PCHitTreeGetter.isValid() ){
363  mf::LogInfo("MuondEdxAna") << "No hits for tree with label "
364  << fParams.qualTreeHitsName();
365  return false; }
366  if( !PCHitXYGetter.isValid() ){
367  mf::LogInfo("MuondEdxAna") << "No XY hits with label " << fParams.qualXYName();
368  return false; }
369  if( !PCHitAvgGetter.isValid() ){
370  mf::LogInfo("MuondEdxAna") << "No Avg hits with label " << fParams.qualAvgName();
371  return false; }
372  if( !PCHitZGetter.isValid() ){
373  mf::LogInfo("MuondEdxAna") << "No Z hits with label " << fParams.qualZName();
374  return false; }
375 
376 
377 
378  // loop over tracks in event
379  for( unsigned int i=0; i<tracks->size(); ++i ){
380 
381  //track level counters:
382  int windowTricells = 0;
383  int windowAvgHits = 0;
384  int windowZHits = 0;
385  int windowPlanes = 0;
386  std::vector<bool> trackWindowPlanes (1000, false); //vector to prevent double counting planes
387 
388  bool windowXView = false;
389  bool windowYView = false;
390 
391  //make custom track length to check against inbuilt length
392  TVector3 startXYZ = (*tracks)[i].Start();
393  TVector3 stopXYZ = (*tracks)[i].Stop();
394  double trackLength = (startXYZ-stopXYZ).Mag();
395 
396  // needed for track window, but good generally to reject protons
397  if( (*tracks)[i].TotalLength() < 200. ) continue;
398 
400  int stopPlaneID = 10000;
401  int firstMuCPlane = -1;
402  bool isNearDet = false;
403 
404  //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.
406  distanceToLiveDetEdgeZ = fLiveGeom->DistToEdgeZ((*tracks)[i].Stop());
407 
408  //continue if track ends in the muon catcher
410  isNearDet = true;
411  cid = fGeo->CellId(stopXYZ[0], stopXYZ[1], stopXYZ[2], 1.5, 1.5, 6., 1.0);
412  if(cid) stopPlaneID = fGeo->getPlaneID(cid);
413  firstMuCPlane = fGeo->FirstPlaneInMuonCatcher();
414  if(stopPlaneID >= firstMuCPlane) continue;
415  }
416 
418  fDistToFront_Stop -> Fill(livegeom->DistToFront(stopXYZ));
419  fDistToBack_Stop -> Fill(livegeom->DistToBack(stopXYZ));
420 
421  // xy and z conditions are mutually exclusive,
422  // avg hits are just the rest in the track
423  // so xy + z + avg is all the hits in the track
424  std::vector<art::Ptr<caldp::PCHit> > xyhits = PCHitXYGetter.at(i);
425  std::vector<art::Ptr<caldp::PCHit> > avghits = PCHitAvgGetter.at(i);
426  std::vector<art::Ptr<caldp::PCHit> > zhits = PCHitZGetter.at(i);
427  nhits = xyhits.size() + avghits.size() + zhits.size();
428 
429  // hits to loop through and fill tree
430  std::vector<art::Ptr<caldp::PCHit> > hits = PCHitTreeGetter.at(i);
431 
432  // ---------------- Fill vector containing hits per APD ---------------------- //
433  std::vector<int> APDhits (100,0);
434  std::vector< std::vector<int> > DiblockAPDhits (100, std::vector<int> (100,0));
435  badTNS = false; //temporary handle for showery events in old data
436 
437  // tricell hits
438  for( unsigned int n=0; n<xyhits.size(); ++n ){
439  APDhits[xyhits[n]->APD()] += 1;
440  DiblockAPDhits.at(xyhits[n]->Diblock()).at(xyhits[n]->APD()) += 1;
441  if(xyhits[n]->TNS() < 0 || xyhits[n]->TNS() > 1400) badTNS = true;
442  }
443  //avg. hits
444  for( unsigned int p=0; p<avghits.size(); ++p ) {
445  APDhits[avghits[p]->APD()] += 1;
446  DiblockAPDhits.at(avghits[p]->Diblock()).at(avghits[p]->APD()) += 1;
447  if(avghits[p]->TNS() < 0 || avghits[p]->TNS() > 1400) badTNS = true;
448  FlightLength = avghits[p]->FlightLen();
449  TotalTrackLength = (*tracks)[i].TotalLength();
450  x = trackLength - FlightLength;
451  if(x<100 || x>200) continue;
452  windowAvgHits ++; //track level counters
453  }
454  //z hits
455  for( unsigned int q=0; q<zhits.size(); ++q ){
456  APDhits[zhits[q]->APD()] += 1;
457  DiblockAPDhits.at(zhits[q]->Diblock()).at(zhits[q]->APD()) += 1;
458  if(zhits[q]->TNS() < 0 || zhits[q]->TNS() > 1400) badTNS = true;
459  FlightLength = zhits[q]->FlightLen();
460  TotalTrackLength = (*tracks)[i].TotalLength();
461  x = trackLength - FlightLength;
462  if(x<100 || x>200) continue;
463  windowZHits ++; //track level counters
464  }
465  // ---------------- Filled APD vector -------------------------------------------- //
466 
467  //cosine giving initial direction of track
468  dCosX = (*tracks)[i].Dir().X()/(*tracks)[i].Dir().Mag();
469  dCosY = (*tracks)[i].Dir().Y()/(*tracks)[i].Dir().Mag();
470  dCosZ = (*tracks)[i].Dir().Z()/(*tracks)[i].Dir().Mag();
471 
472  for( unsigned int j=0; j<hits.size(); ++j ){
473 
474  if(isNearDet)
475  if(hits[j]->Plane() >= firstMuCPlane) continue;
476 
477  /*
478  int p = hits[j]->Plane();
479  int c = hits[j]->Cell();
480  rb::CellHit chit;
481  chit.SetPlane( p );
482  chit.SetCell ( c );
483  //double tns_calib = hits[j]->TNS() + fCalib->GetTimingOffset(p,c,e.isRealData());
484  double tns_calib = hits[j]->TNS();
485  chit.SetTNS(tns_calib,hits[j]->GoodTime()); // may need to get offset from calibrator...
486  if(!e.isRealData())chit.SetMC(); // assume only filled for MC...
487  const std::vector<sim::FLSHit> flss = bt->HitToFLSHit(chit);
488  int lastpdg=0;
489  bool uniformPdg=true;
490  for(const sim::FLSHit& fls: flss){
491  if(lastpdg!=fls.GetPDG() && lastpdg!=0){
492  std::cout << "multiple pdg in flshit: "
493  << lastpdg << " vs. " << fls.GetPDG() << std::endl;
494  uniformPdg = false;
495  }
496  lastpdg=fls.GetPDG();
497  }
498  if(flss.size()>0 && uniformPdg)
499  pdg = flss[0].GetPDG();
500  else pdg = -1;
501  */
502 
503  //test poisson lambda:
504  //double poissonLambda = hits[j]->PoissonLambda();
505  //std::cout<< "poissonLambda: " << poissonLambda <<std::endl;
506 
507  nXYHits = hits.size();
508  diblock = hits[j]->Diblock();
509  view = (int)hits[j]->View();
510  FlightLength = hits[j]->FlightLen();
511  TotalTrackLength= (*tracks)[i].TotalLength();
512  x = trackLength - FlightLength;
513  plane = hits[j]->Plane();
514  path = hits[j]->Path();
515  truepath = hits[j]->TruePath();
516  trueW = hits[j]->TrueW();
517  trueE = hits[j]->TrueMeV();
518  PE = hits[j]->PE();
519  truePE = hits[j]->TruePE();
520  PECorr = this->getPECorr(e,*hits[j]);
521  PECorrToGeV = this->getPECorrToGeV(e,*hits[j]);
523  tns = hits[j]->TNS();
524  cell = hits[j]->Cell();
525  w = hits[j]->W();
526  APD = hits[j]->APD();
527  sameAPDhitsEvent= APDhits[APD];
528  sameAPDhitDiblockEvent = DiblockAPDhits[diblock][APD];
530 
531  fTree->Fill(); //fill ttree
532 
533  fdPEdx ->Fill(x,PE/path);
534  fdEdx ->Fill(x,PECorr/path);
535  fdEdx_minus->Fill(x-fParams.endPSystShift(),PECorr/path);
536  fdEdx_plus ->Fill(x+fParams.endPSystShift(),PECorr/path);
537  fdEdx_true ->Fill(x,trueE/path);
538 
539  if(w > 0){
542  }
543 
544  if(x<100 || x>200) continue;
545 
546  //track level counters
547  windowTricells ++;
548  if(view == 0) windowXView = true;
549  if(view == 1) windowYView = true;
550  //check if plane already counted
551  if(trackWindowPlanes[plane] == false) windowPlanes ++;
552  trackWindowPlanes[plane] = true;
553 
555  fillLowLevel(PECorr, PE, trueE, RecoGeV, path, w, plane, cell, view, diblock, x, trackLength, FlightLength, TotalTrackLength, dCosX, dCosY, dCosZ);
556 
557  } // loop through hits
558 
559  fWindowTricellsPerTrack -> Fill(windowTricells);
560  fWindowPlanesPerTrack -> Fill(windowPlanes);
561  fWindowHitsPerTrack -> Fill(windowTricells + windowAvgHits + windowZHits);
562  fWinTricellsHitsPerTrack -> Fill(windowTricells,windowTricells + windowAvgHits + windowZHits);
563  fWinTricellsPlanesPerTrack -> Fill(windowTricells,windowPlanes);
564  fWinHitsPlanesPerTrack -> Fill(windowTricells + windowAvgHits + windowZHits, windowPlanes);
565 
566  int nViewsInWindow = 0;
567  if(windowXView == true) nViewsInWindow++;
568  if(windowYView == true) nViewsInWindow++;
569 
570  fWinTricellsViewsPerTrack -> Fill(windowTricells,nViewsInWindow);
571  fWindowViewsPerTrack -> Fill(nViewsInWindow);
572  if(nViewsInWindow == 1){
573  if(windowXView) fWindowOneViewTracks -> Fill(0);
574  if(windowYView) fWindowOneViewTracks -> Fill(1);
575  }
576  }
577 
578  return true;
579  }
580 
581  void MuondEdxAna::fillLowLevel(double PECorr, double PE, double trueE, double RecoGeV, double path, double w, int plane, int cell,
582  int view, int diblock, double x, double trackLength, double FlightLength, double TotalTrackLength,
583  double dCosX, double dCosY, double dCosZ){
584 
585  fDCosX -> Fill(dCosX);
586  fDCosY -> Fill(dCosY);
587  fDCosZ -> Fill(dCosZ);
588  fDiblock -> Fill(diblock);
589  fPECorr -> Fill(PECorr/path);
590  fPE -> Fill(PE/path);
591  fTrueE -> Fill(trueE/path);
592  fRecoGeV -> Fill(1000 * RecoGeV/path);
593 
594 
595  /*
596  if( (
597  fGeo->DetId() == novadaq::cnv::kNEARDET &&
598  (
599  (view==0 && fMinFlatW_ND_x < w && w < fMaxFlatW_ND_x ) ||
600  (view==1 && fMinFlatW_ND_y < w && w < fMaxFlatW_ND_y )
601  )
602  ) || (
603  fGeo->DetId() == novadaq::cnv::kFARDET &&
604  (
605  (view==0 && fMinFlatW_FD_x < w && w < fMaxFlatW_FD_x ) ||
606  (view==1 && fMinFlatW_FD_y < w && w < fMaxFlatW_FD_y )
607  )
608  )
609  ){
610  fMEUreco_flat_w->Fill(PECorr/path);
611  fMEUtrue_flat_w->Fill(trueE/path);
612  }
613  */
614 
615  if(view == 0){
616  fPECorr_X -> Fill(PECorr/path);
617  fPE_X -> Fill(PE/path);
618  fTrueE_X -> Fill(trueE/path);
619  }
620 
621  if(view == 1){
622  fPECorr_Y -> Fill(PECorr/path);
623  fPE_Y -> Fill(PE/path);
624  fTrueE_Y -> Fill(trueE/path);
625  }
626 
627  fTrackEndDist -> Fill(x);
628  fTrackEndDistCustom-> Fill(trackLength - FlightLength);
629  fFlightLength -> Fill(FlightLength);
630  fTrackLength -> Fill(TotalTrackLength);
631  fTrackLengthCustom -> Fill(trackLength);
632 
633  }
634 
635  void MuondEdxAna::fillRespVar(double PECorr, double PE, double trueE, double RecoGeV, double path, double w,
636  int plane, int cell, int view, double tns, double EventTimeHigh, int run, double dCosX,
637  double dCosY, double dCosZ, int diblock){
638 
639  fPathLength -> Fill(path);
640  fPlane -> Fill(plane);
641  fW -> Fill(w);
642  fPECorrPlane -> Fill(plane, PECorr/path);
643  fMevPlane -> Fill(plane, trueE/path);
644  fPECorrW -> Fill(w, PECorr/path);
645  fMevW -> Fill(w, trueE/path);
646  fPECorrPath -> Fill(path,PECorr/path);
647  fMevPath -> Fill(path,trueE/path);
648  fPECorrDiblock -> Fill(diblock, PECorr/path);
649  fMevDiblock -> Fill(diblock, trueE/path);
650  fTNS -> Fill(tns);
651  fPECorrTNS -> Fill(tns,PECorr/path);
652  fRun -> Fill(run);
653  fPECorrRun -> Fill(run,PECorr/path);
654  fEventTime -> Fill(EventTimeHigh);
655  fPECorrEventTime-> Fill(EventTimeHigh,PECorr/path);
656  fRespDCosX -> Fill(dCosX, PECorr/path);
657  fRespDCosY -> Fill(dCosY, PECorr/path);
658  fRespDCosZ -> Fill(dCosZ, PECorr/path);
659 
660  if(view == 0){
661  fPathW_X -> Fill(w, path);
662  fPECorrPathW_X -> Fill(w, path,PECorr/path);
663  fMeVPathW_X -> Fill(w, path,trueE/path);
664  fPathLength_X -> Fill(path);
665  fPECorrPath_X -> Fill(path,PECorr/path);
666  fMevPath_X -> Fill(path,trueE/path);
667  fPlane_X -> Fill(plane);
668  fW_X -> Fill(w);
669  fPECorrPlane_X -> Fill(plane, PECorr/path);
670  fMevPlane_X -> Fill(plane, trueE/path);
671  fPECorrW_X -> Fill(w, PECorr/path);
672  fMevW_X -> Fill(w, trueE/path);
673  fCell_X -> Fill(cell);
674  fPECorrCell_X -> Fill(cell, PECorr/path);
675  fMevCell_X -> Fill(cell, trueE/path);
676  fWvsCell_PECorr_0 -> Fill(cell, w, PECorr/path);
677  fWvsCell_MEV_0 -> Fill(cell, w, trueE/path);
678  fWvsCell_0 -> Fill(cell, w);
679 
680  }
681 
682  if(view == 1){
683  fPathW_Y -> Fill(w, path);
684  fPECorrPathW_Y -> Fill(w, path,PECorr/path);
685  fMeVPathW_Y -> Fill(w, path,trueE/path);
686  fPathLength_Y -> Fill(path);
687  fPECorrPath_Y -> Fill(path,PECorr/path);
688  fMevPath_Y -> Fill(path,trueE/path);
689  fPlane_Y -> Fill(plane);
690  fW_Y -> Fill(w);
691  fPECorrPlane_Y -> Fill(plane, PECorr/path);
692  fMevPlane_Y -> Fill(plane, trueE/path);
693  fPECorrW_Y -> Fill(w, PECorr/path);
694  fMevW_Y -> Fill(w, trueE/path);
695  fCell_Y -> Fill(cell);
696  fPECorrCell_Y -> Fill(cell, PECorr/path);
697  fMevCell_Y -> Fill(cell, trueE/path);
698  fWvsCell_PECorr_1-> Fill(cell, w, PECorr/path);
699  fWvsCell_MEV_1 -> Fill(cell, w, trueE/path);
700  fWvsCell_1 -> Fill(cell, w);
701 
702  }
703 
704  }
705 
707 
708  //Make mock cellhit
709  rb::CellHit cellhit;
710  cellhit.SetPlane(pchit.Plane());
711  cellhit.SetCell(pchit.Cell());
712  cellhit.SetView(pchit.View());
713  cellhit.SetTNS(pchit.TNS(),pchit.GoodTime());
714  if(!e.isRealData())cellhit.SetMC(); // assume only filled for MC...
715  cellhit.SetPE(pchit.PE());
716 
717  double pecorr = fCalib->GetPECorr(cellhit, pchit.W());
718  return pecorr;
719  }
720 
722 
723  rb::CellHit cellhit;
724  if(!e.isRealData()) cellhit.SetMC(); // assume only filled for MC...
725 
726  double pecorrtogev = fCalib->GetPECorrToGeVScale(cellhit);
727  return pecorrtogev;
728  }
729 
730 
732  {
733  // histograms for output
735 
736  fTree = tfs->make<TTree>("fTree","tree to hold muondedx variables");
737  fTree->Branch("nXYHits",&nXYHits);
738  fTree->Branch("nhits",&nhits);
739  fTree->Branch("dCosX",&dCosX);
740  fTree->Branch("dCosY",&dCosY);
741  fTree->Branch("dCosZ",&dCosZ);
742  fTree->Branch("x",&x);
743  fTree->Branch("badTNS",&badTNS);
744  fTree->Branch("run",&run);
745  fTree->Branch("diblock",&diblock);
746  fTree->Branch("view",&view);
747  fTree->Branch("FlightLength",&FlightLength);
748  fTree->Branch("TotalTrackLength",&TotalTrackLength);
749  fTree->Branch("xCust",&xCust);
750  fTree->Branch("distanceToLiveDetEdgeZ",&distanceToLiveDetEdgeZ);
751  fTree->Branch("plane",&plane);
752  fTree->Branch("path",&path);
753  fTree->Branch("truepath",&truepath);
754  fTree->Branch("trueW",&trueW);
755  fTree->Branch("trueE",&trueE);
756  fTree->Branch("truePE",&truePE);
757  fTree->Branch("PE",&PE);
758  fTree->Branch("PECorr",&PECorr);
759  fTree->Branch("tns",&tns);
760  fTree->Branch("cell",&cell);
761  fTree->Branch("w",&w);
762  fTree->Branch("pdg",&pdg);
763  fTree->Branch("PECorrToGeV",&PECorrToGeV);
764  fTree->Branch("RecoGeV",&RecoGeV);
765  fTree->Branch("EventTimeHigh",&EventTimeHigh);
766  fTree->Branch("evt_year",&evt_year);
767  fTree->Branch("evt_month",&evt_month);
768  fTree->Branch("evt_day",&evt_day);
769  fTree->Branch("evt_hour",&evt_hour);
770  fTree->Branch("evt_min",&evt_min);
771  fTree->Branch("evt_sec",&evt_sec);
772  fTree->Branch("APD",&APD);
773  fTree->Branch("sameAPDhitsEvent",&sameAPDhitsEvent);
774  fTree->Branch("sameAPDhitDiblockEvent",&sameAPDhitDiblockEvent);
775  fTree->Branch("fibBrightnessBin",&fibBrightnessBin);
776 
777  fDeltaStartPlane = tfs->make<TH1D>("DeltaStartPlane",";startPlane diff. (x,y);",100,0,100);
778  fDeltaEndPlane = tfs->make<TH1D>("DeltaEndPlane",";endPlane diff. (x,y);",100,0,100);
779  fPlaneAsymmetry = tfs->make<TH1D>("PlaneAsymmetry",";plane asymm. (x,y);",100,0,100);
780  fPlaneAsymmetry2D = tfs->make<TH2D>("PlaneAsymmetry2D",";diff in view-planes(x-y); sum of view-planes (x+y);",100,0,100,1000,0,1000);
781 
782  fEvtTime = tfs->make<TH1D>("EvtTime",";EvtTime (YYYYMMDD);",std::pow(10,6),0,std::pow(10,6));
783  fDetRespEvtTime = tfs->make<TProfile>("DetRespEvtTime",";EvtTime (YYYYMMDD); PECorr/cm",std::pow(10,6),0,std::pow(10,6));
784 
785  fDCosX = tfs->make<TH1D>("DCosX",";dCosX;",200,-1,1);
786  fDCosY = tfs->make<TH1D>("DCosY",";dCosY;",200,-1,1);
787  fDCosZ = tfs->make<TH1D>("DCosZ",";dCosZ;",200,-1,1);
788  fRespDCosX = tfs->make<TProfile>("RespDCosX",";dCosX;",200,-1,1);
789  fRespDCosY = tfs->make<TProfile>("RespDCosY",";dCosY;",200,-1,1);
790  fRespDCosZ = tfs->make<TProfile>("RespDCosZ",";dCosZ;",200,-1,1);
791 
792  fdPEdx = tfs->make<TH2D>("dPEdx",";Distance to track end (cm);PE / cm;",
793  50,0.,500.,200,0.,200.);
794  fdEdx = tfs->make<TH2D>("dEdx",";Distance to track end (cm);PECorr / cm;",
795  50,0.,500.,200,0.,200.);
796  fdEdx_minus = tfs->make<TH2D>("dEdx_minus",";Distance to track end (cm);PECorr / cm;",
797  50,0.,500.,200,0.,200.);
798  fdEdx_plus = tfs->make<TH2D>("dEdx_plus",";Distance to track end (cm);PECorr / cm;",
799  50,0.,500.,200,0.,200.);
800  fdEdx_true = tfs->make<TH2D>("dEdx_true",";Distance to track end (cm);True MeV / cm;",
801  50,0.,500.,200,0.,5.);
802 
803  fdEdx_MuonCatch = tfs->make<TH2D>("dEdx_MuonCatch",";Distance to track end (cm);PECorr / cm;",
804  50,0.,500.,200,0.,200.);
805  fdEdx_true_MuonCatch = tfs->make<TH2D>("dEdx_true_MuonCatch",";Distance to track end (cm);True MeV / cm;",
806  50,0.,500.,200,0.,5.);
807 
808  fdEdx_NearReadout = tfs->make<TH2D>("dEdx_NearReadout",";Distance to track end (cm);PECorr / cm;",
809  50,0.,500.,200,0.,200.);
810  fdEdx_true_NearReadout = tfs->make<TH2D>("dEdx_true_NearReadout",";Distance to track end (cm);True MeV / cm;",
811  50,0.,500.,200,0.,5.);
812 
813 
814  //fMEUreco_flat_w = tfs->make<TH1D>("meureco_flat_w",";PECorr/cm;",1000,0,500);
815  //fMEUtrue_flat_w = tfs->make<TH1D>("meutrue_flat_w",";MeV/cm;",1000,0,25);
816 
817  fPECorr = tfs->make<TH1D>("PECorr",";PECorr/cm;",1000,0,500);
818  fPE = tfs->make<TH1D>("PE",";PE/cm;",1000,0,500);
819  fTrueE = tfs->make<TH1D>("TrueE",";MeV/cm;",1000,0,25);
820  fRecoGeV = tfs->make<TH1D>("RecoGeV",";MeV/cm;",1000,0,25);
821  fPECorr_X = tfs->make<TH1D>("PECorr_X",";PECorr/cm;",1000,0,500);
822  fPE_X = tfs->make<TH1D>("PE_X",";PE/cm;",1000,0,500);
823  fTrueE_X = tfs->make<TH1D>("TrueE_X",";MeV/cm;",1000,0,25);
824  fPECorr_Y = tfs->make<TH1D>("PECorr_Y",";PECorr/cm;",1000,0,500);
825  fPE_Y = tfs->make<TH1D>("PE_Y",";PE/cm;",1000,0,500);
826  fTrueE_Y = tfs->make<TH1D>("TrueE_Y",";MeV/cm;",1000,0,25);
827 
828  fFlightLength = tfs->make<TH1D>("FlightLength",";distance from track start (cm);",1000,0,10000);
829  fTrackLength = tfs->make<TH1D>("TrackLength",";Total track length (cm);",1000,0,10000);
830  fTrackEndDist = tfs->make<TH1D>("TrackEndDist",";Distance to end of track (cm)",1000,0,10000);
831  fTrackLengthCustom = tfs->make<TH1D>("TrackLengthCustom",";Total track length (cm);",1000,0,10000);
832  fTrackEndDistCustom = tfs->make<TH1D>("TrackEndDistCustom",";Distance to end of track (cm)",1000,0,10000);
833  fTrackEndPlane = tfs->make<TH1D>("TrackEndPlane",";plane of track end;",1000,0,1000);
834  fStopPlane = tfs->make<TH1D>("StopPlane",";start and stop z positions;",1000,0,1000);
835  fStopPlane_liveGeomCut = tfs->make<TH1D>("StopPlane_liveGeomCut",";start and stop z positions;",1000,0,1000);
836  fDistToFront_Stop = tfs->make<TH1D>("DistToFront_Stop",";distance to front (cm);",1000,0,1000);
837  fDistToBack_Stop = tfs->make<TH1D>("DistToBack_Stop",";distance to back (cm);",1000,0,1000);
838 
839  fTNS = tfs->make<TH1D>("TNS",";TNS (nova time);",std::pow(10,3),-5*std::pow(10,4),5*std::pow(10,4));
840  fPECorrTNS = tfs->make<TProfile>("PECorrTNS",";TNS (nova time); PECorr/cm;",std::pow(10,3),-5*std::pow(10,4),5*std::pow(10,4));
841  fRun = tfs->make<TH1D>("Run",";Run;",10000,0,100000);
842  fPECorrRun = tfs->make<TProfile>("PECorrRun",";Run; PECorr/cm;",10000,0,100000);
843 
844  fEventTime = tfs->make<TH1D>("EventTime",";date and time (seconds);",std::pow(10,6),0,std::pow(10,10));
845  fPECorrEventTime = tfs->make<TProfile>("PECorrEventTime",";date and time (seconds); PECorr/cm;",std::pow(10,6),0,std::pow(10,10));
846 
847  fPathW_X = tfs->make<TH2D>("PathW_X",";W (cm); Path (cm)",160,-800,800, 100, 2, 12);
848  fPathW_Y = tfs->make<TH2D>("PathW_Y",";W (cm); Path (cm)",160,-800,800, 100, 2, 12);
849  fPECorrPathW_X = tfs->make<TProfile2D>("PECorrPathW_X",";W (cm); Path (cm); PECorr/cm",160,-800,800, 100, 2, 12);
850  fPECorrPathW_Y = tfs->make<TProfile2D>("PECorrPathW_Y",";W (cm); Path (cm); PECorr/cm",160,-800,800, 100, 2, 12);
851  fMeVPathW_X = tfs->make<TProfile2D>("MeVPathW_X",";W (cm); Path (cm); MeV/cm",160,-800,800, 100, 2, 12);
852  fMeVPathW_Y = tfs->make<TProfile2D>("MeVPathW_Y",";W (cm); Path (cm); MeV/cm",160,-800,800, 100, 2, 12);
853 
854  fPathLength = tfs->make<TH1D>("PathLength",";path length (cm);",200,0,20);
855  fPathLength_X = tfs->make<TH1D>("PathLength_X",";path length (cm);",200,0,20);
856  fPathLength_Y = tfs->make<TH1D>("PathLength_Y",";path length (cm);",200,0,20);
857  fPECorrPath = tfs->make<TProfile>("PECorrPath",";path length X-view (cm); PECorr/cm;",200,0,20);
858  fPECorrPath_X = tfs->make<TProfile>("PECorrPath_X",";path length X-view (cm); PECorr/cm;",200,0,20);
859  fPECorrPath_Y = tfs->make<TProfile>("PECorrPath_Y",";path length X-view (cm); PECorr/cm;",200,0,20);
860  fMevPath = tfs->make<TProfile>("MevPath",";path length X-view (cm); Mev/cm;",200,0,20);
861  fMevPath_X = tfs->make<TProfile>("MevPath_X",";path length X-view (cm); Mev/cm;",200,0,20);
862  fMevPath_Y = tfs->make<TProfile>("MevPath_Y",";path length X-view (cm); Mev/cm;",200,0,20);
863 
864  fPlane = tfs->make<TH1D>("Plane",";plane;",1000,0,1000);
865  fPlane_X = tfs->make<TH1D>("Plane_X",";plane;",1000,0,1000);
866  fPlane_Y = tfs->make<TH1D>("Plane_Y",";plane;",1000,0,1000);
867  fPECorrPlane = tfs->make<TProfile>("PECorrPlane",";plane; PECorr;",1000,0,1000);
868  fPECorrPlane_X = tfs->make<TProfile>("PECorrPlane_X",";plane; PECorr;",1000,0,1000);
869  fPECorrPlane_Y = tfs->make<TProfile>("PECorrPlane_Y",";plane; PECorr;",1000,0,1000);
870  fMevPlane = tfs->make<TProfile>("MevPlane",";plane; Mev;",1000,0,1000);
871  fMevPlane_X = tfs->make<TProfile>("MevPlane_X",";plane; Mev;",1000,0,1000);
872  fMevPlane_Y = tfs->make<TProfile>("MevPlane_Y",";plane; Mev;",1000,0,1000);
873 
874  fW = tfs->make<TH1D>("W",";w;",160,-800,800);
875  fW_X = tfs->make<TH1D>("W_X",";w;",160,-800,800);
876  fW_Y = tfs->make<TH1D>("W_Y",";w;",160,-800,800);
877  fPECorrW = tfs->make<TProfile>("PECorrW",";w; PECorr/cm;",160,-800,800);
878  fPECorrW_X = tfs->make<TProfile>("PECorrW_X",";w; PECorr/cm;",160,-800,800);
879  fPECorrW_Y = tfs->make<TProfile>("PECorrW_Y",";w; PECorr/cm;",160,-800,800);
880  fMevW = tfs->make<TProfile>("MevW",";w; Mev/cm;",160,-800,800);
881  fMevW_X = tfs->make<TProfile>("MevW_X",";w; Mev/cm;",160,-800,800);
882  fMevW_Y = tfs->make<TProfile>("MevW_Y",";w; Mev/cm;",160,-800,800);
883 
884  fCell_X = tfs->make<TH1D>("Cell_X",";cell;",1000,0,1000);
885  fCell_Y = tfs->make<TH1D>("Cell_Y",";cell;",1000,0,1000);
886  fPECorrCell_X = tfs->make<TProfile>("PECorrCell_X",";cell; PECorr/cm;",1000,0,1000);
887  fPECorrCell_Y = tfs->make<TProfile>("PECorrCell_Y",";cell; PECorr/cm;",1000,0,1000);
888  fMevCell_X = tfs->make<TProfile>("MevCell_X",";cell; Mev/cm;",1000,0,1000);
889  fMevCell_Y = tfs->make<TProfile>("MevCell_Y",";cell; Mev/cm;",1000,0,1000);
890 
891  //track level histograms
892  fWindowTricellsPerTrack = tfs->make<TH1D>("WindowTricellsPerTrack",";tricells per track in track window;",50,0,50);
893  fWindowHitsPerTrack = tfs->make<TH1D>("WindowHitsPerTrack",";hits per track in track window;",1000,0,1000);
894  fWindowPlanesPerTrack = tfs->make<TH1D>("WindowPlanesPerTrack",";planes per track in track window;",1000,0,1000);
895  fWindowViewsPerTrack = tfs->make<TH1D>("WindowViewsPerTrack",";views per track in track window;",10,0,10);
896  fWindowOneViewTracks = tfs->make<TH1D>("WindowOneViewTracks",";View in window when only one present;",10,0,10);
897 
898  fWinTricellsHitsPerTrack = tfs->make<TH2D>("WinTricellsHitsPerTrack",";tricells per track in track window;hits per track in window;",
899  50,0,50,60,0,60);
900  fWinTricellsPlanesPerTrack = tfs->make<TH2D>("WinTricellsPlanesPerTrack",";tricells per track in track window;planes hit per track in window;",
901  50,0,50,10,0,10);
902  fWinHitsPlanesPerTrack = tfs->make<TH2D>("WinHitsPlanesPerTrack",";hits per track in window;planes per track in window",
903  60,0,60,10,0,10);
904  fWinTricellsViewsPerTrack = tfs->make<TH2D>("WinTricellsViewsPerTrack",";tricells per track in track window;views per track in window;",
905  50,0,50,4,0,4);
906 
907  fWvsCell_0 = tfs->make<TH2D>("WvsCell_0",";Cell;W;",
908  40,0.,400.,160,-800.,800.);
909  fWvsCell_1 = tfs->make<TH2D>("WvsCell_1",";Cell;W;",
910  40,0.,400.,160,-800.,800.);
911  fWvsCell_PECorr_0 = tfs->make<TH3D>("WvsCell_PECorr_0",";Cell;W;",
912  40,0.,400.,160,-800.,800.,100,0.,100.);
913  fWvsCell_PECorr_1 = tfs->make<TH3D>("WvsCell_PECorr_1",";Cell;W;",
914  40,0.,400.,160,-800.,800.,100,0.,100.);
915  fWvsCell_MEV_0 = tfs->make<TH3D>("WvsCell_MEV_0",";Cell;W;",
916  40,0.,400.,160,-800.,800.,100,0.,10.);
917  fWvsCell_MEV_1 = tfs->make<TH3D>("WvsCell_MEV_1",";Cell;W;",
918  40,0.,400.,160,-800.,800.,100,0.,10.);
919 
920  fDiblock = tfs->make<TH1D>("Diblock","diblock",20,0,20);
921  fPECorrDiblock = tfs->make<TProfile>("PECorrDiblock",";diblock; PECorr/cm",20,0,20);
922  fMevDiblock = tfs->make<TProfile>("MevDiblock",";diblock; MeV/cm",20,0,20);
923 
924  ftrackEndPlaneWithinDiblock = tfs->make<TH1D>("trackEndPlaneWithinDiblock", ";plane of track end within diblocks;",100,0,100);
925 
926  }
927 
929  return true;
930  }
931 
933 
935  std::unique_ptr<TH2D> pdEdx(tfs->make<TH2D>(*fdEdx));
936  std::unique_ptr<TH2D> pdEdx_minus(tfs->make<TH2D>(*fdEdx_minus));
937  std::unique_ptr<TH2D> pdEdx_plus(tfs->make<TH2D>(*fdEdx_plus));
938  std::unique_ptr<TH2D> pdEdx_true(tfs->make<TH2D>(*fdEdx_true));
939 
940  r.put(std::move(pdEdx),"dEdx");
941  r.put(std::move(pdEdx_minus),"dEdxMinus");
942  r.put(std::move(pdEdx_plus),"dEdxPlus");
943  r.put(std::move(pdEdx_true),"dEdxTrue");
944  return true;
945  }
946 
948  {
949 
950  }
951 }// end namespace
float TNS() const
Return uncorrected hit time.
Definition: PCHit.h:52
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:54
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
constexpr T pow(T x)
Definition: pow.h:75
const uint32_t NOVA_EPOCH
int Cell() const
Return cell value.
Definition: PCHit.h:26
double DistToFront(TVector3 vertex)
double DistToBack(TVector3 vertex)
bool isRealData() const
Definition: Event.h:83
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
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: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
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
void SetCell(unsigned short cell)
Definition: CellHit.h:52
art::ServiceHandle< calib::Calibrator > fCalib
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()
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))
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
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: View.py:1
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
T * make(ARGS...args) const
void beginJob() override
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
void endJob() override
double GetPECorr(rb::CellHit const &cellhit, double w)
art::ServiceHandle< geo::Geometry > fGeo
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
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.