LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PhotonCounterT0Matching_module.cc
Go to the documentation of this file.
1 
32 // Framework includes
37 #include "art_root_io/TFileService.h"
40 #include "fhiclcpp/ParameterSet.h"
41 
42 #include <iostream>
43 #include <memory>
44 
45 // LArSoft
55 
56 // ROOT
57 #include "TH1D.h"
58 #include "TH2D.h"
59 #include "TTree.h"
60 
61 namespace lbne {
63 }
64 
66 public:
68  // The destructor generated by the compiler is fine for classes
69  // without bare pointers or other resource use.
70 
71  // Plugins should not be copied or assigned.
76 
77  // Required functions.
78  void produce(art::Event& e) override;
79 
80  // Selected optional functions.
81  void beginJob() override;
82 
83 private:
84  // Internal functions.....
85  void TrackProp(double TrackStart_X,
86  double TrackEnd_X,
87  double& TrackLength_X,
88  double& TrackCentre_X,
89  double TrackStart_Y,
90  double TrackEnd_Y,
91  double& TrackLength_Y,
92  double& TrackCentre_Y,
93  double TrackStart_Z,
94  double TrackEnd_Z,
95  double& TrackLength_Z,
96  double& TrackCentre_Z,
97  double trkTimeStart,
98  double trkTimeEnd,
99  double& trkTimeLengh,
100  double& trkTimeCentre,
101  double& TrackLength);
102  double DistFromPoint(double StartY,
103  double EndY,
104  double StartZ,
105  double EndZ,
106  double PointY,
107  double PointZ);
108 
109  // Params got from fcl file.......
110  std::string fTrackModuleLabel;
111  std::string fShowerModuleLabel;
112  std::string fHitsModuleLabel;
113  std::string fFlashModuleLabel;
114  std::string fTruthT0ModuleLabel;
116  double fPredictedXPower = 1;
122  double fPEThreshold;
124 
125  // Variables used in module.......
127  double TrackLength_Y, TrackCentre_Y /* , BestTrackCentre_Y */;
128  double TrackLength_Z, TrackCentre_Z /* , BestTrackCentre_Z */;
130 
142 
143  double YZSep, MCTruthT0;
144  // Histograms in TFS branches
145  TTree* fTree;
146  TH2D* hPredX_T;
147  TH2D* hPredX_PE;
148  TH2D* hPredX_T_PE;
155 };
156 
158 {
159  // Call appropriate produces<>() functions here.
160  produces<std::vector<anab::T0>>();
161  produces<art::Assns<recob::Track, anab::T0>>();
162  produces<art::Assns<recob::Shower, anab::T0>>();
163 
164  fTrackModuleLabel = (p.get<std::string>("TrackModuleLabel"));
165  fShowerModuleLabel = (p.get<std::string>("ShowerModuleLabel"));
166  fHitsModuleLabel = (p.get<std::string>("HitsModuleLabel"));
167  fFlashModuleLabel = (p.get<std::string>("FlashModuleLabel"));
168  fTruthT0ModuleLabel = (p.get<std::string>("TruthT0ModuleLabel"));
169 
170  fPredictedXConstant = (p.get<double>("PredictedXConstant"));
171  fPredictedExpConstant = (p.get<double>("PredictedExpConstant"));
172  fPredictedExpGradient = (p.get<double>("PredictedExpGradient"));
173 
174  fDriftWindowSize = (p.get<double>("DriftWindowSize"));
175  fWeightOfDeltaYZ = (p.get<double>("WeightOfDeltaYZ"));
176  fMatchCriteria = (p.get<double>("MatchCriteria"));
177  fPEThreshold = (p.get<double>("PEThreshold"));
178 
179  fVerbosity = (p.get<bool>("Verbose", false));
180 }
181 
183 {
184  // Implementation of optional member function here.
186  fTree = tfs->make<TTree>("PhotonCounterT0Matching", "PhotonCounterT0");
187  fTree->Branch("TrackCentre_X", &BestTrackCentre_X, "TrackCentre_X/D");
188  fTree->Branch("PredictedX", &BestPredictedX, "PredictedX/D");
189  fTree->Branch("TrackTimeCent", &BesttrkTimeCentre, "TimeSepPredX/D");
190  fTree->Branch("FlashTime", &BestFlashTime, "FlashTime/D");
191  fTree->Branch("TimeSep", &BestTimeSep, "TimeSep/D");
192  fTree->Branch("TimeSepPredX", &BestTimeSepPredX, "TimeSepPredX/D");
193  fTree->Branch("minYZSep", &BestminYZSep, "minYZSep/D");
194  fTree->Branch("FitParam", &BestFitParam, "FitParam/D");
195  fTree->Branch("MCTruthT0", &MCTruthT0, "MCTruthT0/D");
196 
197  hPredX_T = tfs->make<TH2D>("hPredX_T",
198  "Predicted X from timing information against reconstructed X; "
199  "Reconstructed X (cm); Predicted X (cm)",
200  30,
201  0,
202  300,
203  30,
204  0,
205  300);
206  hPredX_PE = tfs->make<TH2D>("hPredX_PE",
207  "Predicted X from PE information against reconstructed X; "
208  "Reconstructed X (cm); Predicted X (cm)",
209  30,
210  0,
211  300,
212  30,
213  0,
214  300);
215  hPredX_T_PE = tfs->make<TH2D>("hPredX_T_PE",
216  "Predicted X position from time and PE information; Predicted X "
217  "from timing information (cm); Predicted X from PE information",
218  60,
219  0,
220  300,
221  60,
222  0,
223  300);
224  hdeltaX_deltaYZ = tfs->make<TH2D>(
225  "hdeltaX_deltaYZ",
226  "Difference between X predicted from PE's and T agaisnt distance of flash from track in YZ; "
227  "Difference in X predicted from PE's and T (cm); Distance of flash from track in YZ (cm)",
228  40,
229  0,
230  200,
231  40,
232  0,
233  100);
234  hdeltaYZ_Length = tfs->make<TH2D>("hdeltaYZ_Length",
235  "Distance of flash from track against track length; Distance "
236  "from flash to track (cm); Track length (cm)",
237  20,
238  0,
239  100,
240  60,
241  0,
242  300);
244  tfs->make<TH2D>("hFitParam_Length",
245  "How fit correlates with track length; Fit correlation; Track Length (cm)",
246  50,
247  0,
248  250,
249  30,
250  0,
251  300);
252  hPhotonT0_MCT0 = tfs->make<TH2D>("hPhotonT0_MCT0",
253  "Comparing Photon Counter reconstructed T0 against MCTruth T0; "
254  "Photon Counter T0 (us); MCTruthT0 T0 (us)",
255  1760,
256  -1600,
257  16000,
258  1760,
259  -1600,
260  16000);
261  hT0_diff_full = tfs->make<TH1D>(
262  "hT0_diff_full",
263  "Difference between MCTruth T0 and photon detector T0; Time difference (us); Number",
264  500,
265  -2500,
266  2500);
267  hT0_diff_zoom = tfs->make<TH1D>(
268  "hT0_diff_zoom",
269  "Difference between MCTruth T0 and photon detector T0; Time difference (us); Number",
270  320,
271  -1.6,
272  1.6);
273 }
274 
276 {
277  // Access art services...
279  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
280  auto const detprop =
282 
283  //TrackList handle
284  art::Handle<std::vector<recob::Track>> trackListHandle;
285  std::vector<art::Ptr<recob::Track>> tracklist;
286  if (evt.getByLabel(fTrackModuleLabel, trackListHandle))
287  art::fill_ptr_vector(tracklist, trackListHandle);
288 
289  //ShowerList handle
290  art::Handle<std::vector<recob::Shower>> showerListHandle;
291  std::vector<art::Ptr<recob::Shower>> showerlist;
292  if (evt.getByLabel(fShowerModuleLabel, showerListHandle))
293  art::fill_ptr_vector(showerlist, showerListHandle);
294 
295  //HitList Handle
297  std::vector<art::Ptr<recob::Hit>> hitlist;
298  if (evt.getByLabel(fHitsModuleLabel, hitListHandle)) art::fill_ptr_vector(hitlist, hitListHandle);
299 
300  //FlashList Handle
302  std::vector<art::Ptr<recob::OpFlash>> flashlist;
303  if (evt.getByLabel(fFlashModuleLabel, flashListHandle))
304  art::fill_ptr_vector(flashlist, flashListHandle);
305 
306  // Create anab::T0 objects and make association with recob::Track
307 
308  std::unique_ptr<std::vector<anab::T0>> T0col(new std::vector<anab::T0>);
309  std::unique_ptr<art::Assns<recob::Track, anab::T0>> Trackassn(
311  std::unique_ptr<art::Assns<recob::Shower, anab::T0>> Showerassn(
313 
314  if (trackListHandle.isValid() && flashListHandle.isValid()) {
315  //Access tracks and hits
316  art::FindManyP<recob::Hit> fmtht(trackListHandle, evt, fTrackModuleLabel);
317  art::FindMany<anab::T0> fmtruth(trackListHandle, evt, fTruthT0ModuleLabel);
318 
319  size_t NTracks = tracklist.size();
320  size_t NFlashes = flashlist.size();
321 
322  if (fVerbosity)
323  std::cout << "There were " << NTracks << " tracks and " << NFlashes
324  << " flashes in this event." << std::endl;
325 
326  // Now to access PhotonCounter for each track...
327  for (size_t iTrk = 0; iTrk < NTracks; ++iTrk) {
328  if (fVerbosity) std::cout << "\n New Track " << (int)iTrk << std::endl;
329  // Reset Variables.
332  bool ValidTrack = false;
333 
334  // Work out Properties of the track.
335  recob::Track::Point_t trackStart, trackEnd;
336  std::tie(trackStart, trackEnd) = tracklist[iTrk]->Extent();
337  std::vector<art::Ptr<recob::Hit>> allHits = fmtht.at(iTrk);
338  size_t nHits = allHits.size();
339  trkTimeStart = allHits[nHits - 1]->PeakTime() /
340  clock_data.TPCClock().Frequency(); // Got in ticks, now in us!
341  trkTimeEnd =
342  allHits[0]->PeakTime() / clock_data.TPCClock().Frequency(); // Got in ticks, now in us!
343  TrackProp(trackStart.X(),
344  trackEnd.X(),
347  trackStart.Y(),
348  trackEnd.Y(),
351  trackStart.Z(),
352  trackEnd.Z(),
355  trkTimeStart,
356  trkTimeEnd,
357  trkTimeLengh,
358  trkTimeCentre, // times in us!
359  TrackLength);
360 
361  // Some cout statement about track properties.
362  if (fVerbosity) {
363  std::cout << trackStart.X() << " " << trackEnd.X() << " " << TrackLength_X << " "
364  << TrackCentre_X << "\n"
365  << trackStart.Y() << " " << trackEnd.Y() << " " << TrackLength_Y << " "
366  << TrackCentre_Y << "\n"
367  << trackStart.Z() << " " << trackEnd.Z() << " " << TrackLength_Z << " "
368  << TrackCentre_Z << "\n"
369  << trkTimeStart << " " << trkTimeEnd << " " << trkTimeLengh << " "
370  << trkTimeCentre << std::endl;
371  }
372  // ----- Loop over flashes ------
373  for (size_t iFlash = 0; iFlash < NFlashes; ++iFlash) {
374  //Reset some flash specific quantities
375  YZSep = minYZSep = 9999;
376  FlashTime = TimeSep = 9999;
378  // Check flash could be caused by track...
379  FlashTime = flashlist[iFlash]->Time(); // Got in us!
380  TimeSep = trkTimeCentre - FlashTime; // Time in us!
381  if (TimeSep < 0 || TimeSep > (fDriftWindowSize / clock_data.TPCClock().Frequency()))
382  continue; // Times compared in us!
383 
384  // Check flash has enough PE's to satisfy our threshold
385  if (flashlist[iFlash]->TotalPE() < fPEThreshold) continue;
386 
387  // Work out some quantities for this flash...
388  // PredictedX = ( A / x^n ) + exp ( B + Cx )
389  PredictedX =
390  (fPredictedXConstant / pow(flashlist[iFlash]->TotalPE(), fPredictedXPower)) +
391  (exp(fPredictedExpConstant + (fPredictedExpGradient * flashlist[iFlash]->TotalPE())));
392  TimeSepPredX = TimeSep * detprop.DriftVelocity(); // us * cm/us = cm!
394  // Dependant on each point...
395  for (size_t Point = 1; Point < tracklist[iTrk]->NumberTrajectoryPoints(); ++Point) {
396  auto NewPoint = tracklist[iTrk]->LocationAtPoint(Point);
397  auto PrevPoint = tracklist[iTrk]->LocationAtPoint(Point - 1);
398  YZSep = DistFromPoint(NewPoint.Y(),
399  PrevPoint.Y(),
400  NewPoint.Z(),
401  PrevPoint.Z(),
402  flashlist[iFlash]->YCenter(),
403  flashlist[iFlash]->ZCenter());
404  if (Point == 1) minYZSep = YZSep;
405  if (YZSep < minYZSep) minYZSep = YZSep;
406  }
407 
408  // Determine how well matched this track is......
409  if (fMatchCriteria == 0)
410  FitParam =
411  pow(((DeltaPredX * DeltaPredX) + (minYZSep * minYZSep * fWeightOfDeltaYZ)), 0.5);
412  else if (fMatchCriteria == 1)
413  FitParam = minYZSep;
414  else if (fMatchCriteria == 2)
416 
417  //----FLASH INFO-----
418  if (fVerbosity) {
419  std::cout << "\nFlash " << (int)iFlash << " " << TrackCentre_X << ", " << TimeSepPredX
420  << " - " << PredictedX << " = " << DeltaPredX << ", " << minYZSep << " -> "
421  << FitParam << std::endl;
422  }
423  //----Select best flash------
424  if (FitParam < BestFitParam) {
425  ValidTrack = true;
426  BestFlash = (int)iFlash;
435  BestFlashTime = FlashTime;
437  } // Find best Flash
438  } // Loop over Flashes
439 
440  // ---- Now Make association and fill TTree/Histos with the best matched flash.....
441  if (ValidTrack) {
442 
443  // -- Fill Histos --
450  // ------ Compare Photon Matched to MCTruth Matched -------
451  if (fmtruth.isValid()) {
452  std::vector<const anab::T0*> T0s = fmtruth.at((int)iTrk);
453  for (size_t i = 0; i < T0s.size(); ++i) {
454  MCTruthT0 = T0s[i]->Time() / 1e3; // Got in ns, now in us!!
455  hPhotonT0_MCT0->Fill(BestFlashTime, MCTruthT0);
456  hT0_diff_full->Fill(MCTruthT0 - BestFlashTime);
457  hT0_diff_zoom->Fill(MCTruthT0 - BestFlashTime);
458  }
459  }
460  // -- Fill TTree --
461  fTree->Fill();
462  //Make Association
463  T0col->push_back(anab::T0(
464  BestFlashTime * 1e3, FlashTriggerType, (int)BestFlash, (*T0col).size(), BestFitParam));
465  util::CreateAssn(evt, *T0col, tracklist[iTrk], *Trackassn);
466  } // Valid Track
467  } // Loop over tracks
468  }
469  evt.put(std::move(T0col));
470  evt.put(std::move(Trackassn));
471  evt.put(std::move(Showerassn));
472 
473 } // Produce
474 // ----------------------------------------------------------------------------------------------------------------------------
476  double TrackEnd_X,
477  double& TrackLength_X,
478  double& TrackCentre_X,
479  double TrackStart_Y,
480  double TrackEnd_Y,
481  double& TrackLength_Y,
482  double& TrackCentre_Y,
483  double TrackStart_Z,
484  double TrackEnd_Z,
485  double& TrackLength_Z,
486  double& TrackCentre_Z,
487  double trkTimeStart,
488  double trkTimeEnd,
489  double& trkTimeLengh,
490  double& trkTimeCentre,
491  double& TrackLength)
492 {
494  TrackLength_X = fabs(TrackEnd_X - TrackStart_X);
495  if (TrackStart_X < TrackEnd_X)
496  TrackCentre_X = TrackStart_X + 0.5 * TrackLength_X;
497  else
498  TrackCentre_X = TrackStart_X - 0.5 * TrackLength_X;
499 
500  TrackLength_Y = fabs(TrackEnd_Y - TrackStart_Y);
501  if (TrackStart_Y < TrackEnd_Y)
502  TrackCentre_Y = TrackStart_Y + 0.5 * TrackLength_Y;
503  else
504  TrackCentre_Y = TrackStart_Y - 0.5 * TrackLength_Y;
505 
506  TrackLength_Z = fabs(TrackEnd_Z - TrackStart_Z);
507  if (TrackStart_Z < TrackEnd_Z)
508  TrackCentre_Z = TrackStart_Z + 0.5 * TrackLength_Z;
509  else
510  TrackCentre_Z = TrackStart_Z - 0.5 * TrackLength_Z;
511 
512  trkTimeLengh = trkTimeEnd - trkTimeStart;
513  trkTimeCentre = trkTimeStart + 0.5 * trkTimeLengh;
514 
515  TrackLength = pow(pow((TrackEnd_X - TrackStart_X), 2) + pow((TrackEnd_Y - TrackStart_Y), 2) +
516  pow((TrackEnd_Z - TrackStart_Z), 2),
517  0.5);
518 
519  return;
520 }
521 // ----------------------------------------------------------------------------------------------------------------------------
523  double EndY,
524  double StartZ,
525  double EndZ,
526  double PointY,
527  double PointZ)
528 {
530  double Length = hypot(fabs(EndY - StartY), fabs(EndZ - StartZ));
531  double distance =
532  ((PointZ - StartZ) * (EndY - StartY) - (PointY - StartY) * (EndZ - StartZ)) / Length;
533  return fabs(distance);
534 }
535 // ----------------------------------------------------------------------------------------------------------------------------
float Length(const PFPStruct &pfp)
Definition: PFPUtils.cxx:3269
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Definition: T0.h:16
double DistFromPoint(double StartY, double EndY, double StartZ, double EndZ, double PointY, double PointZ)
bool isValid() const noexcept
Definition: Handle.h:203
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
Provides recob::Track data product.
void TrackProp(double TrackStart_X, double TrackEnd_X, double &TrackLength_X, double &TrackCentre_X, double TrackStart_Y, double TrackEnd_Y, double &TrackLength_Y, double &TrackCentre_Y, double TrackStart_Z, double TrackEnd_Z, double &TrackLength_Z, double &TrackCentre_Z, double trkTimeStart, double trkTimeEnd, double &trkTimeLengh, double &trkTimeCentre, double &TrackLength)
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
PhotonCounterT0Matching(fhicl::ParameterSet const &p)
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:42
tracking::Point_t Point_t
Definition: Track.h:52
PhotonCounterT0Matching & operator=(PhotonCounterT0Matching const &)=delete
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Float_t e
Definition: plot.C:35
art framework interface to geometry description