MEFinder_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 // Michel Electron Finder package
3 // Dan Pershey
4 // pershey@caltech.edu
5 //
6 // Jul 9 2015.
7 // MEFinder is a Michel E tagging package designed to meet all of NOvA's
8 // Michel reconstruction needs. It is a revamp of MichelEFilter and can
9 // satisfy both our calibration and current analysis requirements for Michels.
10 //
11 ///////////////////////////////////////////////////////////////////////////////
12 
13 // Framework includes
21 #include "cetlib/search_path.h"
22 
23 
24 // NOvASoft includes
25 #include "Geometry/Geometry.h"
28 #include "MEFinder/MEClusters.h"
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/Cluster.h"
31 #include "RecoBase/Track.h"
32 #include "Utilities/AssociationUtil.h"
33 #include "Utilities/func/MathUtil.h"
34 
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "TFile.h"
38 #include "TMath.h"
39 
40 namespace me
41 {
42 
43  // These guys come up a lot...
45  typedef std::vector<PtrCHit> VecCHits;
46  typedef std::vector<me::VecCHits> VecVecCHits;
47 
48 
49  // Holds hists for computing MID, the Michel E LL Identifier
50  struct LikeHists
51  {
54 
57  };
58 
59  // Collection of relevant variables to ME cluster
60  struct MEVars
61  {
62  double distToSlc;
64  // Only valid for TrkMEs
66  };
67 
68  // Holds fhicl parameters
70  {
71  template<class T> using Atom = fhicl::Atom<T>;
72  template<class T> using Table = fhicl::Table<T>;
74  using Name = fhicl::Name;
75 
76 
77  // Grab necessary input from other reco modules
78  Atom<std::string> slicerInput
79  {
80  Name("SlicerInput"),
81  Comment("Input for physics slices")
82  };
83  Atom<std::string> trackInput
84  {
85  Name("TrackInput"),
86  Comment("Input for tracks with which to associate MEs")
87  };
88 
89  // Params for deciding if noise hit is good enough to physics slice
90  Atom<double> maxDist
91  {
92  Name("MaxDist"),
93  Comment("Max dist from phsyics slice to selected noise hit")
94  };
95  Atom<double> minHitADC
96  {
97  Name("MinHitADC"),
98  Comment("Min ADC for a hit to be considered - needed for FD speed ")
99  };
100 
101  // Should we calculate MID values for this job
102  Atom<bool> calibAvailable
103  {
104  Name("CalibAvailable"),
105  Comment("Are there calibration constants available for getting CalE")
106  };
107 
108  // Preselection parameters
109  Atom<double> maxClosestApproach
110  {
111  Name("MaxClosestApproach"),
112  Comment("ME clust must be at least this close to a hit in phys slice")
113  };
114  Atom<double> trkClosestApproach
115  {
116  Name("TrkClosestApproach"),
117  Comment("Dist cut to be considered TrkME")
118  };
119  Atom<unsigned int> maxNCells
120  {
121  Name("MaxNCells"),
122  Comment("Presel Cut")
123  };
124  Atom<unsigned int> minNCells
125  {
126  Name("MinNCells"),
127  Comment("Presel Cut")
128  };
129  Atom<double> maxDeltaT
130  {
131  Name("MaxDeltaT"),
132  Comment("Presel Cut")
133  };
134  Atom<double> minDeltaT
135  {
136  Name("MinDeltaT"),
137  Comment("Presel Cut")
138  };
139  Atom<double> maxCalE
140  {
141  Name("MaxCalE"),
142  Comment("Presel Cut")
143  };
144  Atom<double> minCalE
145  {
146  Name("MinCalE"),
147  Comment("Presel Cut")
148  };
149  Atom<double> minDeltaTReTrigger
150  {
151  Name("MinDeltaTReTrigger"),
152  Comment("See DocDB 12295, apd dead time not well modeled in MC")
153  };
154 
155  // DBSCAN parameters
156  Atom<double> eps
157  {
158  Name("EpsForClustering"),
159  Comment("Dist to search out to in the DBSCAN algorithm")
160  };
161  Atom<double> planeScale
162  {
163  Name("PlaneScale"),
164  Comment("For distances in DBSCAN")
165  };
166  Atom<double> cellScale
167  {
168  Name("CellScale"),
169  Comment("For distances in DBSCAN")
170  };
171  Atom<double> timeScale
172  {
173  Name("TimeScale"),
174  Comment("For distances in DBSCAN")
175  };
176 
177  Atom<std::string> llFileName
178  {
179  Name("LLFileName"),
180  Comment("Contains hist files for calculating Delta LL ME PID, MID")
181  };
182 
183  Atom<bool> matchMEtoSlcByDist
184  {
185  Name("MatchMEtoSlcByDist"),
186  Comment("Match MEs to Slices by distance only")
187  };
188 
189  };
190 
191  class MEFinder : public art::EDProducer
192  {
193  public:
194 
195  // Allows 'nova --print-description' to work
197 
198  explicit MEFinder(const Parameters& params);
199  virtual ~MEFinder();
200 
201  void produce(art::Event& evt);
202  void beginJob();
203 
204  protected:
205 
206  // Following functions work to making a set of ME clusters of rb::CellHits
207 
208  // Sorts through noise and small physics hits to find candidate Michel hits
209  void GetPhysNoiseSlcs(art::Handle<std::vector<rb::Cluster> > slcs,
210  std::vector<rb::Cluster> &noiseSlcs,
211  std::vector<rb::Cluster> &physSlcs);
212  me::VecCHits MECandHits(std::vector<rb::Cluster> noiseSlcs,
213  std::vector<rb::Cluster> physSlcs);
214  void FindNoiseHitRange(me::VecCHits hits, double physTNS,
215  int& firstIdx, int& latsIdx);
216  // Adds noise hit to list of Michel hits if it is near a physics slice
217  void AddMEHit(me::VecCHits &meHits,
218  me::PtrCHit noiseCell,
219  rb::Cluster physSlc);
220  // These two find min distance between a hit and a physics slice
221  double MinHitSlcDist(me::PtrCHit hit, rb::Cluster slc);
222  double HitToHitDistance(const PtrCHit &hit1,
223  const PtrCHit &hit2);
224  // Sorts these ME hits into clusters of rb::CellHits
225  me::VecVecCHits ClusterMEHits(VecCHits);
226  // Holds the distance function between points for clustering
227  double MEDistanceMetric(me::PtrCHit hit1, me::PtrCHit hit2);
228 
229 
230  // Following functions make art SlcME/TrkME products and compute variables
231 
232  me::SlcME MakeSlcME(me::VecCHits michel,
233  art::Handle<std::vector<rb::Cluster> > slcs,
234  int& parentSlcIdx);
235  me::TrkME MakeTrkME(me::SlcME michel,
236  art::Handle<std::vector<rb::Cluster> > slcs,
237  art::Ptr<rb::Track> parentTrack);
238  void BestSlcMatch(me::SlcME michel,
239  me::MEVars& vars,
240  art::Handle<std::vector<rb::Cluster> > slcs,
241  int& parentSlcIdx);
242  bool PassesMEPresel(me::SlcME michel);
243  double GetMID(double deltaT, double distToSlc, double calE, int nCells);
244 
245 
246  // Asks if ME is a TrkME, and fills out relevant extra variables
247 
248  bool IsTrkME(me::SlcME michel,
249  const std::vector< art::Ptr<rb::Track> > trks,
250  art::Ptr<rb::Track>& parentTrk);
251  double TrkMEDist(me::SlcME slcME, TVector3 trkEnd);
252  std::vector<double> MEMeanXYZ(me::VecCHits me, bool&, bool&);
253 
254 
255  // Need a hold of the geometry to calculated cell-cell distances
257 
259 
261  };
262 }
263 
264 namespace me
265 {
267  : fParams(params())
268  {
269  produces<std::vector<me::SlcME> >();
270  produces<art::Assns<me::SlcME, rb::Cluster> >();
271  produces<std::vector<me::TrkME> >();
272  produces<art::Assns<me::TrkME, rb::Track> >();
273  produces<art::Assns<me::TrkME, rb::Cluster> >();
274  }
275 
277 
278  // beginJob pulls out hist files for MID
280  {
281  cet::search_path sp("FW_SEARCH_PATH");
282  std::string LikeFileName;
283  sp.find_file(fParams.llFileName(), LikeFileName);
284  TFile* like = TFile::Open(LikeFileName.c_str());
285  like->cd();
286 
287  fLikeHists.fMEDeltaT_Dist = (TH2F*)gDirectory->Get("MEDeltaT_Dist");
288  fLikeHists.fMECalE_NCells = (TH2F*)gDirectory->Get("MECalE_NCells");
289 
290  fLikeHists.fHFDeltaT_Dist = (TH2F*)gDirectory->Get("HFDeltaT_Dist");
291  fLikeHists.fHFCalE_NCells = (TH2F*)gDirectory->Get("HFCalE_NCells");
292  }
293 
294  // Main produce function. Finds ME's and puts them in ART
296  {
297  // Grab ahold of slc and trk info from the art event, set up Assns
299  evt.getByLabel(fParams.slicerInput(), slcs);
300  art::FindManyP<rb::Track> trkInfo(slcs, evt, fParams.trackInput());
301 
302  std::unique_ptr< std::vector<me::SlcME> >
303  slcMEs(new std::vector<me::SlcME>);
304  std::unique_ptr< art::Assns<me::SlcME, rb::Cluster> >
306 
307  std::unique_ptr< std::vector<me::TrkME> >
308  trkMEs(new std::vector<me::TrkME>);
309  std::unique_ptr< art::Assns<me::TrkME, rb::Track> >
311  std::unique_ptr< art::Assns<me::TrkME, rb::Cluster> >
312  trkMEToSlcAssn(new art::Assns<me::TrkME, rb::Cluster>);
313 
314  // Vector of all CellHits close to any parent slice in time and space
315  std::vector<rb::Cluster> noiseSlcs, physSlcs;
316  GetPhysNoiseSlcs(slcs, noiseSlcs, physSlcs);
317  me::VecCHits meHits = MECandHits(noiseSlcs, physSlcs);
318 
319  // Then sort this vec of CellHits into ME clusters
320  me::VecVecCHits meClusters = ClusterMEHits(meHits);
321 
322  for (unsigned int i = 0; i < meClusters.size(); i++){
323  art::Ptr<rb::Track> parentTrk;
324  int parentSlcIdx;
325  me::SlcME slcME = MakeSlcME(meClusters[i], slcs, parentSlcIdx);
326  if (!PassesMEPresel(slcME)) continue;
327  const std::vector< art::Ptr<rb::Track> > trks = trkInfo.at(parentSlcIdx);
328  if (IsTrkME(slcME, trks, parentTrk)){
329  me::TrkME trkME = MakeTrkME(slcME, slcs, parentTrk);
330  trkMEs->push_back(trkME);
331  util::CreateAssn(*this,evt,*trkMEs,trkME.ParentTrk(),*(trkAssn.get()));
332  util::CreateAssn(*this,evt,*trkMEs,trkME.ParentSlc(),
333  *(trkMEToSlcAssn.get()));
334  }
335  else {
336  slcMEs->push_back(slcME);
337  util::CreateAssn(*this,evt,*slcMEs,slcME.ParentSlc(),*(slcAssn.get()));
338  }
339  }
340  evt.put(std::move(slcMEs));
341  evt.put(std::move(slcAssn));
342  evt.put(std::move(trkMEs));
343  evt.put(std::move(trkAssn));
344  evt.put(std::move(trkMEToSlcAssn));
345  }
346 
347  // Outputs a list of noise cells close enough to physics to be a ME hit
348  void MEFinder::GetPhysNoiseSlcs(art::Handle< std::vector<rb::Cluster> > slcs,
349  std::vector<rb::Cluster> &noiseSlcs,
350  std::vector<rb::Cluster> &physSlcs)
351  {
352  // Label each slice either noise or physics. A noise slice can be a slice
353  // flagged as noise, or slices with only a few hits ( <= 12 )
354  for (unsigned int i = 0; i < slcs->size(); i++){
355  rb::Cluster slice = slcs->at(i);
356  if (slice.IsNoise())
357  noiseSlcs.push_back(slice);
358  else if (slice.NCell() <= fParams.maxNCells())
359  noiseSlcs.push_back(slice);
360  else
361  physSlcs.push_back(slice);
362  }
363  }
364 
365  me::VecCHits MEFinder::MECandHits(std::vector<rb::Cluster> noiseSlcs,
366  std::vector<rb::Cluster> physSlcs){
367  // Loop over all noise hits, then ask each physics slice if the noise
368  // hit is close in time / space, then add the hit to the list of candidate
369  // hits appropriately.
370  me::VecCHits meHits;
371 
372  for (unsigned int nIdx = 0; nIdx < noiseSlcs.size(); nIdx++){
373  me::VecCHits noiseHits;
374  for (unsigned int cIdx = 0; cIdx < noiseSlcs[nIdx].NCell(); cIdx++)
375  if (noiseSlcs[nIdx].Cell(cIdx)->ADC() > fParams.minHitADC())
376  noiseHits.push_back(noiseSlcs[nIdx].Cell(cIdx));
377  rb::SortByTime(noiseHits);
378  for (unsigned int pIdx = 0; pIdx < physSlcs.size(); pIdx++){
379  int fIdx = 0; int lIdx = 0;
380  FindNoiseHitRange(noiseHits, physSlcs[pIdx].MeanTNS(), fIdx, lIdx);
381  for (int cIdx = fIdx; cIdx <= lIdx; cIdx++){
382  AddMEHit(meHits, noiseHits[cIdx], physSlcs[pIdx]);
383  }
384  }
385  }
386  return meHits;
387  }
388 
389  // Calculates the first / last noise hit with a delta T in the allowed
390  // region for MEs for a given slice
391  void MEFinder::FindNoiseHitRange(me::VecCHits nHits, double physTNS,
392  int& fIdx, int& lIdx)
393  {
394  fIdx = 0; lIdx = 0;
395  double fMETime = physTNS+fParams.minDeltaT();
396  double lMETime = physTNS+fParams.maxDeltaT();
397  int nIter = log2(nHits.size())+1;
398  for (int i = 1; i <= nIter; i++){
399  int fmidpoint = fIdx + (nHits.size()-1)/pow(2,i);
400  fIdx = (nHits[fmidpoint]->TNS() > fMETime) ? fIdx : fmidpoint;
401  int lmidpoint = lIdx + nHits.size()/pow(2,i);
402  lIdx = (nHits[lmidpoint]->TNS() > lMETime) ? lIdx : lmidpoint;
403  }
404  fIdx++;
405  }
406 
407  // Adds current noiseCell to list of ME hits if close to physSlc
409  me::PtrCHit noiseCell,
410  rb::Cluster physSlc){
411  double minDist = MinHitSlcDist(noiseCell, physSlc);
412  double deltaT = noiseCell->TNS() - physSlc.MeanTNS();
413  if (noiseCell->ADC() < fParams.minHitADC())
414  return;
415  // Treat Retriggers differently because of APD dead-time (DocDB 12295)
416  // Additionally require deltaT > 0 for tagging retrigger hits to
417  // accomodate negative search time for measuring uncorrelated background
418  double curMinDeltaT = ( minDist == 0 && deltaT > 0 ) ?
420  if (deltaT > curMinDeltaT &&
421  deltaT < fParams.maxDeltaT() &&
422  minDist < fParams.maxDist()){
423  bool addHit = true;
424  for (int i = 0; i < (int)meHits.size(); i++){
425  if (noiseCell->Plane() == meHits[i]->Plane() &&
426  noiseCell->Cell() == meHits[i]->Cell() &&
427  noiseCell->ADC() == meHits[i]->ADC())
428  addHit = false;
429  }
430  if (addHit) meHits.push_back(noiseCell);
431  }
432  return;
433  }
434 
435  // Sorts these hits into clusters
437  me::VecVecCHits meClusters;
438  if (meHits.size() == 0)
439  return meClusters;
440 
441  // This will keep track of wheter hit is already clustered
442  std::vector<bool> isUsed(meHits.size(), false);
443 
444  for (unsigned int i = 0; i < meHits.size(); i++){
445  if (isUsed[i]) continue;
446  me::VecCHits curCluster;
447  curCluster.push_back(meHits[i]);
448  isUsed[i] = true;
449  for (unsigned int j = i; j < meHits.size(); j++){
450  if (isUsed[j]) continue;
451  if (MEDistanceMetric(meHits[i], meHits[j]) < fParams.eps()){
452  isUsed[j] = true;
453  curCluster.push_back(meHits[j]);
454  // Check if any of the hits we just passed over are close to jth idx
455  for (unsigned int k = i; k < j; k++){
456  if (isUsed[k]) continue;
457  if (MEDistanceMetric(meHits[j], meHits[k]) < fParams.eps()){
458  isUsed[k] = true;
459  curCluster.push_back(meHits[k]);
460  }
461  }
462  }
463  }
464  meClusters.push_back(curCluster);
465  }
466  return meClusters;
467  }
468 
470  {
471  double deltaC = (hit1->View() == hit2->View()) ?
472  double(hit1->Cell() - hit2->Cell()) / fParams.cellScale() : 0;
473  double deltaP = double(hit1->Plane() - hit2->Plane());
474  deltaP /= fParams.planeScale();
475  double deltaT = double(hit1->TNS() - hit2->TNS());
476  deltaT /= fParams.timeScale();
477  return sqrt(deltaC*deltaC + deltaP*deltaP + deltaT*deltaT);
478  }
479 
480  // Min dist between a ME hit and a physics slice
482  double minDist = 10000;
483  for (unsigned int i = 0; i < slc.NCell(); i++){
484  double curDist = HitToHitDistance(hit, slc.Cell(i));
485  if (curDist < minDist){
486  minDist = curDist;
487  }
488  }
489  return minDist;
490  }
491 
492  // Calculates dist between two cells
493  double MEFinder::HitToHitDistance(const PtrCHit &hit1, const PtrCHit &hit2){
494  // Only calculate if both hits are in the same view
495  if (hit1->View() != hit2->View()) return 100000.;
496  else{
497  double xyz1[3], xyz2[3];
498 
499  unsigned short plane1_idx = hit1->Plane();
500  unsigned short cell1_idx = hit1->Cell();
501  const geo::PlaneGeo* plane1 = fgeom->Plane(plane1_idx);
502  const geo::CellGeo* cell1 = plane1->Cell(cell1_idx);
503  cell1->GetCenter(xyz1);
504 
505  unsigned short plane2_idx = hit2->Plane();
506  unsigned short cell2_idx = hit2->Cell();
507  const geo::PlaneGeo* plane2 = fgeom->Plane(plane2_idx);
508  const geo::CellGeo* cell2 = plane2->Cell(cell2_idx);
509  cell2->GetCenter(xyz2);
510 
511  // Agnostic to whether hits are in X or Y views
512  return util::pythag(xyz1[0]-xyz2[0], xyz1[1]-xyz2[1], xyz1[2]-xyz2[2]);
513  }
514  }
515 
516 
518  art::Handle<std::vector<rb::Cluster> > slcs,
519  int& parentSlcIdx){
520  me::SlcME slcME;
522 
523  for (unsigned int i = 0; i < michel.size(); i++)
524  slcME.Add(michel[i]);
525  BestSlcMatch(slcME, vars, slcs, parentSlcIdx);
526  slcME.SetDistToSlc (vars.distToSlc);
527  slcME.SetParentSlc (vars.parentSlc);
528  slcME.SetDeltaT (slcME.MeanTNS() - vars.parentSlc->MeanTNS());
529  double mid = -10;
530  if (fParams.calibAvailable()){
531  mid = GetMID(slcME.DeltaT(),
532  slcME.DistToSlc(),
533  slcME.CalorimetricEnergy(),
534  slcME.NCell());
535  }
536  slcME.SetMID(mid);
537  return slcME;
538  }
539 
541  art::Handle<std::vector<rb::Cluster> > slcs,
542  art::Ptr<rb::Track> parentTrk){
543  me::TrkME trkME;
545 
546  for (unsigned int i = 0; i < slcME.NCell(); i++)
547  trkME.Add(slcME.Cell(i));
548 
549  trkME.SetDistToSlc (slcME.DistToSlc());
550  trkME.SetParentSlc (parentTrk);
551  trkME.SetDeltaT (slcME.MeanTNS() - parentTrk->MeanTNS());
552  double mid = -10;
553  if (fParams.calibAvailable()){
554  mid = GetMID(slcME.DeltaT(),
555  slcME.DistToSlc(),
556  slcME.CalorimetricEnergy(),
557  slcME.NCell());
558  }
559  trkME.SetMID(mid);
560  trkME.SetParentTrk(parentTrk);
561  return trkME;
562  }
563 
564  // Takes a ME Cluster, and finds out what physics slice it is most likely
565  // associated with using deltaT and distToSlc info
567  me::MEVars& vars,
568  art::Handle<std::vector<rb::Cluster> > slcs,
569  int& parentSlcIdx){
570  // Grab the MID building histogram
571  TH2F* hAssn = fLikeHists.fMEDeltaT_Dist;
572  // The assnScore is larger for slices that are appropriately near the ME
573  double assnScore = -1;
574  parentSlcIdx = 0;
575  double allSlcMinDist = 1.0e6;
576  for (unsigned int i = 0; i < slcs->size(); i++){
577  rb::Cluster slc = slcs->at(i);
578  if (slc.IsNoise() || slc.NCell() <= fParams.maxNCells())
579  continue;
580  art::Ptr<rb::Cluster> slcPtr(slcs, i);
581 
582  double deltaT = slcME.MeanTNS() - slc.MeanTNS();
583  double minDist = 10000;
584 
585  for (unsigned int j = 0; j < slcME.NCell(); j++){
586  double dist = MinHitSlcDist(slcME.Cell(j), slcs->at(i));
587  if (dist < minDist)
588  minDist = dist;
589  } // Loop over ME hits
590 
591  // Match by distance only
592  if ( fParams.matchMEtoSlcByDist() )
593  {
594  if ( minDist < allSlcMinDist )
595  {
596  parentSlcIdx = i;
597  allSlcMinDist = minDist;
598  vars.parentSlc = slcPtr;
599  vars.distToSlc = minDist;
600  }
601  }
602 
603  // Match by distance and time (standard method)
604  else
605  {
606  double curScore = hAssn->GetBinContent(hAssn->FindBin(deltaT, minDist));
607  if (curScore > assnScore){
608  parentSlcIdx = i;
609  assnScore = curScore;
610  vars.parentSlc = slcPtr;
611  vars.distToSlc = minDist;
612  }
613  }
614 
615  } // Loop over slices
616  }
617 
619  if (slcME.DistToSlc() > fParams.maxClosestApproach())
620  return false;
621  if (slcME.NCell() > fParams.maxNCells())
622  return false;
623  if (slcME.NCell() < fParams.minNCells())
624  return false;
625  if (slcME.DeltaT() > fParams.maxDeltaT())
626  return false;
627  if (slcME.DeltaT() < fParams.minDeltaT())
628  return false;
629  // Only apply calE presel if a calibration is available
630  if (fParams.calibAvailable()){
631  if (slcME.CalorimetricEnergy() > fParams.maxCalE())
632  return false;
633  if (slcME.CalorimetricEnergy() < fParams.minCalE())
634  return false;
635  }
636  return true;
637  }
638 
639  double MEFinder::GetMID(double DeltaT, double Dist, double CalE, int NCells){
640  TH2* meDeltaT_Dist = fLikeHists.fMEDeltaT_Dist;
641  TH2* meCalE_NCells = fLikeHists.fMECalE_NCells;
642 
643  int meDeltaT_Dist_Bin = meDeltaT_Dist->FindBin(DeltaT, Dist);
644  int meCalE_NCells_Bin = meCalE_NCells->FindBin(CalE, NCells);
645  double meLike = meDeltaT_Dist->GetBinContent(meDeltaT_Dist_Bin) *
646  meCalE_NCells->GetBinContent(meCalE_NCells_Bin);
647 
648  TH2* hfDeltaT_Dist = fLikeHists.fHFDeltaT_Dist;
649  TH2* hfCalE_NCells = fLikeHists.fHFCalE_NCells;
650 
651  int hfDeltaT_Dist_Bin = hfDeltaT_Dist->FindBin(DeltaT, Dist);
652  int hfCalE_NCells_Bin = hfCalE_NCells->FindBin(CalE, NCells);
653  double hfLike = hfDeltaT_Dist->GetBinContent(hfDeltaT_Dist_Bin) *
654  hfCalE_NCells->GetBinContent(hfCalE_NCells_Bin);
655 
656  if (meLike == 0 && hfLike == 0) return -10;
657  if (meLike == 0)
658  meLike = 1e-10;
659  if (hfLike == 0)
660  hfLike = 1e-10;
661 
662  return TMath::Log(meLike) - TMath::Log(hfLike);
663 
664  }
665 
667  const std::vector< art::Ptr<rb::Track> > trks,
668  art::Ptr<rb::Track>& parentTrk){
669  for (unsigned int i = 0; i < trks.size(); i++){
670  if (trks[i]->NXCell() == 0 || trks[i]->NYCell() == 0)
671  continue;
672  double dist = TrkMEDist(slcME, trks[i]->Stop());
673  if (dist < fParams.trkClosestApproach()){
674  parentTrk = trks[i];
675  return true;
676  }
677  }
678  return false;
679  }
680 
681  double MEFinder::TrkMEDist(me::SlcME me, TVector3 trkEnd){
682  double deltaX = ( (me.NXCell() == 0) ? 0 : me.MeanX() - trkEnd.X() );
683  double deltaY = ( (me.NYCell() == 0) ? 0 : me.MeanY() - trkEnd.Y() );
684  double deltaZ = me.MeanZ() - trkEnd.Z();
685  return util::pythag( deltaX, deltaY, deltaZ );
686  }
687 }
688 
689 namespace me
690 {
692 }
float TNS() const
Definition: CellHit.h:46
void SetParentSlc(art::Ptr< rb::Cluster >)
Definition: MEClusters.cxx:26
void produce(art::Event &evt)
me::VecCHits MECandHits(std::vector< rb::Cluster > noiseSlcs, std::vector< rb::Cluster > physSlcs)
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void AddMEHit(me::VecCHits &meHits, me::PtrCHit noiseCell, rb::Cluster physSlc)
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
art::ServiceHandle< geo::Geometry > fgeom
Atom< double > cellScale
LikeHists fLikeHists
double DistToSlc() const
Definition: MEClusters.cxx:12
art::Ptr< rb::CellHit > PtrCHit
void SetDistToSlc(double)
Definition: MEClusters.cxx:22
unsigned short Plane() const
Definition: CellHit.h:39
art::Ptr< rb::Cluster > ParentSlc() const
Definition: MEClusters.cxx:15
geo::View_t View() const
Definition: CellHit.h:41
Atom< double > maxCalE
double GetMID(double deltaT, double distToSlc, double calE, int nCells)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Atom< double > maxDist
std::vector< PtrCHit > VecCHits
A collection of associated CellHits.
Definition: Cluster.h:47
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:233
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
virtual ~MEFinder()
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Atom< double > timeScale
std::string find_file(std::string const &filename) const
double TrkMEDist(me::SlcME slcME, TVector3 trkEnd)
Atom< double > planeScale
art::Ptr< rb::Track > ParentTrk() const
Definition: MEClusters.cxx:32
double MEDistanceMetric(me::PtrCHit hit1, me::PtrCHit hit2)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
double MeanX(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:231
double dist
Definition: runWimpSim.h:113
void SetDeltaT(double)
Definition: MEClusters.cxx:20
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Atom< bool > calibAvailable
Atom< double > eps
art::Ptr< rb::Track > parentTrk
unsigned short Cell() const
Definition: CellHit.h:40
Atom< double > maxDeltaT
Atom< double > minCalE
Atom< unsigned int > maxNCells
double DeltaT() const
Definition: MEClusters.cxx:10
void hits()
Definition: readHits.C:15
Definition: NueSkimmer.h:24
double MeanY(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:232
Atom< unsigned int > minNCells
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
void beginJob()
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
void BestSlcMatch(me::SlcME michel, me::MEVars &vars, art::Handle< std::vector< rb::Cluster > > slcs, int &parentSlcIdx)
Atom< std::string > trackInput
bool IsTrkME(me::SlcME michel, const std::vector< art::Ptr< rb::Track > > trks, art::Ptr< rb::Track > &parentTrk)
MEFinderParams fParams
const std::map< std::pair< std::string, std::string >, Variable > vars
const double j
Definition: BetheBloch.cxx:29
bool PassesMEPresel(me::SlcME michel)
std::vector< me::VecCHits > VecVecCHits
void SetMID(double)
Definition: MEClusters.cxx:18
Atom< std::string > llFileName
double distToSlc
Atom< double > minDeltaT
Atom< double > minDeltaTReTrigger
Atom< double > trkClosestApproach
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
Atom< double > maxClosestApproach
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
void SetParentTrk(art::Ptr< rb::Track >)
Definition: MEClusters.cxx:44
MEFinder(const Parameters &params)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
Atom< double > minHitADC
me::SlcME MakeSlcME(me::VecCHits michel, art::Handle< std::vector< rb::Cluster > > slcs, int &parentSlcIdx)
me::TrkME MakeTrkME(me::SlcME michel, art::Handle< std::vector< rb::Cluster > > slcs, art::Ptr< rb::Track > parentTrack)
Atom< std::string > slicerInput
void GetPhysNoiseSlcs(art::Handle< std::vector< rb::Cluster > > slcs, std::vector< rb::Cluster > &noiseSlcs, std::vector< rb::Cluster > &physSlcs)
Atom< bool > matchMEtoSlcByDist
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Float_t e
Definition: plot.C:35
Encapsulate the cell geometry.
Definition: CellGeo.h:25
art::Ptr< rb::Cluster > parentSlc
fvar< T > log2(const fvar< T > &x)
Definition: log2.hpp:19
void FindNoiseHitRange(me::VecCHits hits, double physTNS, int &firstIdx, int &latsIdx)
double MinHitSlcDist(me::PtrCHit hit, rb::Cluster slc)
double HitToHitDistance(const PtrCHit &hit1, const PtrCHit &hit2)
Encapsulate the geometry of one entire detector (near, far, ndos)
me::VecVecCHits ClusterMEHits(VecCHits)