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  : EDProducer(params)
268  , fParams(params())
269  {
270  produces<std::vector<me::SlcME> >();
271  produces<art::Assns<me::SlcME, rb::Cluster> >();
272  produces<std::vector<me::TrkME> >();
273  produces<art::Assns<me::TrkME, rb::Track> >();
274  produces<art::Assns<me::TrkME, rb::Cluster> >();
275  }
276 
278 
279  // beginJob pulls out hist files for MID
281  {
282  cet::search_path sp("FW_SEARCH_PATH");
283  std::string LikeFileName;
284  sp.find_file(fParams.llFileName(), LikeFileName);
285  TFile* like = TFile::Open(LikeFileName.c_str());
286  like->cd();
287 
288  fLikeHists.fMEDeltaT_Dist = (TH2F*)gDirectory->Get("MEDeltaT_Dist");
289  fLikeHists.fMECalE_NCells = (TH2F*)gDirectory->Get("MECalE_NCells");
290 
291  fLikeHists.fHFDeltaT_Dist = (TH2F*)gDirectory->Get("HFDeltaT_Dist");
292  fLikeHists.fHFCalE_NCells = (TH2F*)gDirectory->Get("HFCalE_NCells");
293  }
294 
295  // Main produce function. Finds ME's and puts them in ART
297  {
298  // Grab ahold of slc and trk info from the art event, set up Assns
300  evt.getByLabel(fParams.slicerInput(), slcs);
301  art::FindManyP<rb::Track> trkInfo(slcs, evt, fParams.trackInput());
302 
303  std::unique_ptr< std::vector<me::SlcME> >
304  slcMEs(new std::vector<me::SlcME>);
305  std::unique_ptr< art::Assns<me::SlcME, rb::Cluster> >
307 
308  std::unique_ptr< std::vector<me::TrkME> >
309  trkMEs(new std::vector<me::TrkME>);
310  std::unique_ptr< art::Assns<me::TrkME, rb::Track> >
312  std::unique_ptr< art::Assns<me::TrkME, rb::Cluster> >
313  trkMEToSlcAssn(new art::Assns<me::TrkME, rb::Cluster>);
314 
315  // Vector of all CellHits close to any parent slice in time and space
316  std::vector<rb::Cluster> noiseSlcs, physSlcs;
317  GetPhysNoiseSlcs(slcs, noiseSlcs, physSlcs);
318  me::VecCHits meHits = MECandHits(noiseSlcs, physSlcs);
319 
320  // Then sort this vec of CellHits into ME clusters
321  me::VecVecCHits meClusters = ClusterMEHits(meHits);
322 
323  for (unsigned int i = 0; i < meClusters.size(); i++){
324  art::Ptr<rb::Track> parentTrk;
325  int parentSlcIdx;
326  me::SlcME slcME = MakeSlcME(meClusters[i], slcs, parentSlcIdx);
327  if (!PassesMEPresel(slcME)) continue;
328  const std::vector< art::Ptr<rb::Track> > trks = trkInfo.at(parentSlcIdx);
329  if (IsTrkME(slcME, trks, parentTrk)){
330  me::TrkME trkME = MakeTrkME(slcME, slcs, parentTrk);
331  trkMEs->push_back(trkME);
332  util::CreateAssn(evt,*trkMEs,trkME.ParentTrk(),*(trkAssn.get()));
333  util::CreateAssn(evt,*trkMEs,trkME.ParentSlc(),
334  *(trkMEToSlcAssn.get()));
335  }
336  else {
337  slcMEs->push_back(slcME);
338  util::CreateAssn(evt,*slcMEs,slcME.ParentSlc(),*(slcAssn.get()));
339  }
340  }
341  evt.put(std::move(slcMEs));
342  evt.put(std::move(slcAssn));
343  evt.put(std::move(trkMEs));
344  evt.put(std::move(trkAssn));
345  evt.put(std::move(trkMEToSlcAssn));
346  }
347 
348  // Outputs a list of noise cells close enough to physics to be a ME hit
349  void MEFinder::GetPhysNoiseSlcs(art::Handle< std::vector<rb::Cluster> > slcs,
350  std::vector<rb::Cluster> &noiseSlcs,
351  std::vector<rb::Cluster> &physSlcs)
352  {
353  // Label each slice either noise or physics. A noise slice can be a slice
354  // flagged as noise, or slices with only a few hits ( <= 12 )
355  for (unsigned int i = 0; i < slcs->size(); i++){
356  rb::Cluster slice = slcs->at(i);
357  if (slice.IsNoise())
358  noiseSlcs.push_back(slice);
359  else if (slice.NCell() <= fParams.maxNCells())
360  noiseSlcs.push_back(slice);
361  else
362  physSlcs.push_back(slice);
363  }
364  }
365 
366  me::VecCHits MEFinder::MECandHits(std::vector<rb::Cluster> noiseSlcs,
367  std::vector<rb::Cluster> physSlcs){
368  // Loop over all noise hits, then ask each physics slice if the noise
369  // hit is close in time / space, then add the hit to the list of candidate
370  // hits appropriately.
371  me::VecCHits meHits;
372 
373  for (unsigned int nIdx = 0; nIdx < noiseSlcs.size(); nIdx++){
374  me::VecCHits noiseHits;
375  for (unsigned int cIdx = 0; cIdx < noiseSlcs[nIdx].NCell(); cIdx++)
376  if (noiseSlcs[nIdx].Cell(cIdx)->ADC() > fParams.minHitADC())
377  noiseHits.push_back(noiseSlcs[nIdx].Cell(cIdx));
378  rb::SortByTime(noiseHits);
379  for (unsigned int pIdx = 0; pIdx < physSlcs.size(); pIdx++){
380  int fIdx = 0; int lIdx = 0;
381  FindNoiseHitRange(noiseHits, physSlcs[pIdx].MeanTNS(), fIdx, lIdx);
382  for (int cIdx = fIdx; cIdx <= lIdx; cIdx++){
383  AddMEHit(meHits, noiseHits[cIdx], physSlcs[pIdx]);
384  }
385  }
386  }
387  return meHits;
388  }
389 
390  // Calculates the first / last noise hit with a delta T in the allowed
391  // region for MEs for a given slice
392  void MEFinder::FindNoiseHitRange(me::VecCHits nHits, double physTNS,
393  int& fIdx, int& lIdx)
394  {
395  fIdx = 0; lIdx = 0;
396  double fMETime = physTNS+fParams.minDeltaT();
397  double lMETime = physTNS+fParams.maxDeltaT();
398  int nIter = log2(nHits.size())+1;
399  for (int i = 1; i <= nIter; i++){
400  int fmidpoint = fIdx + (nHits.size()-1)/pow(2,i);
401  fIdx = (nHits[fmidpoint]->TNS() > fMETime) ? fIdx : fmidpoint;
402  int lmidpoint = lIdx + nHits.size()/pow(2,i);
403  lIdx = (nHits[lmidpoint]->TNS() > lMETime) ? lIdx : lmidpoint;
404  }
405  fIdx++;
406  }
407 
408  // Adds current noiseCell to list of ME hits if close to physSlc
410  me::PtrCHit noiseCell,
411  rb::Cluster physSlc){
412  double minDist = MinHitSlcDist(noiseCell, physSlc);
413  double deltaT = noiseCell->TNS() - physSlc.MeanTNS();
414  if (noiseCell->ADC() < fParams.minHitADC())
415  return;
416  // Treat Retriggers differently because of APD dead-time (DocDB 12295)
417  // Additionally require deltaT > 0 for tagging retrigger hits to
418  // accomodate negative search time for measuring uncorrelated background
419  double curMinDeltaT = ( minDist == 0 && deltaT > 0 ) ?
421  if (deltaT > curMinDeltaT &&
422  deltaT < fParams.maxDeltaT() &&
423  minDist < fParams.maxDist()){
424  bool addHit = true;
425  for (int i = 0; i < (int)meHits.size(); i++){
426  if (noiseCell->Plane() == meHits[i]->Plane() &&
427  noiseCell->Cell() == meHits[i]->Cell() &&
428  noiseCell->ADC() == meHits[i]->ADC())
429  addHit = false;
430  }
431  if (addHit) meHits.push_back(noiseCell);
432  }
433  return;
434  }
435 
436  // Sorts these hits into clusters
438  me::VecVecCHits meClusters;
439  if (meHits.size() == 0)
440  return meClusters;
441 
442  // This will keep track of wheter hit is already clustered
443  std::vector<bool> isUsed(meHits.size(), false);
444 
445  for (unsigned int i = 0; i < meHits.size(); i++){
446  if (isUsed[i]) continue;
447  me::VecCHits curCluster;
448  curCluster.push_back(meHits[i]);
449  isUsed[i] = true;
450  for (unsigned int j = i; j < meHits.size(); j++){
451  if (isUsed[j]) continue;
452  if (MEDistanceMetric(meHits[i], meHits[j]) < fParams.eps()){
453  isUsed[j] = true;
454  curCluster.push_back(meHits[j]);
455  // Check if any of the hits we just passed over are close to jth idx
456  for (unsigned int k = i; k < j; k++){
457  if (isUsed[k]) continue;
458  if (MEDistanceMetric(meHits[j], meHits[k]) < fParams.eps()){
459  isUsed[k] = true;
460  curCluster.push_back(meHits[k]);
461  }
462  }
463  }
464  }
465  meClusters.push_back(curCluster);
466  }
467  return meClusters;
468  }
469 
471  {
472  double deltaC = (hit1->View() == hit2->View()) ?
473  double(hit1->Cell() - hit2->Cell()) / fParams.cellScale() : 0;
474  double deltaP = double(hit1->Plane() - hit2->Plane());
475  deltaP /= fParams.planeScale();
476  double deltaT = double(hit1->TNS() - hit2->TNS());
477  deltaT /= fParams.timeScale();
478  return sqrt(deltaC*deltaC + deltaP*deltaP + deltaT*deltaT);
479  }
480 
481  // Min dist between a ME hit and a physics slice
483  double minDist = 10000;
484  for (unsigned int i = 0; i < slc.NCell(); i++){
485  double curDist = HitToHitDistance(hit, slc.Cell(i));
486  if (curDist < minDist){
487  minDist = curDist;
488  }
489  }
490  return minDist;
491  }
492 
493  // Calculates dist between two cells
494  double MEFinder::HitToHitDistance(const PtrCHit &hit1, const PtrCHit &hit2){
495  // Only calculate if both hits are in the same view
496  if (hit1->View() != hit2->View()) return 100000.;
497  else{
498  double xyz1[3], xyz2[3];
499 
500  unsigned short plane1_idx = hit1->Plane();
501  unsigned short cell1_idx = hit1->Cell();
502  const geo::PlaneGeo* plane1 = fgeom->Plane(plane1_idx);
503  const geo::CellGeo* cell1 = plane1->Cell(cell1_idx);
504  cell1->GetCenter(xyz1);
505 
506  unsigned short plane2_idx = hit2->Plane();
507  unsigned short cell2_idx = hit2->Cell();
508  const geo::PlaneGeo* plane2 = fgeom->Plane(plane2_idx);
509  const geo::CellGeo* cell2 = plane2->Cell(cell2_idx);
510  cell2->GetCenter(xyz2);
511 
512  // Agnostic to whether hits are in X or Y views
513  return util::pythag(xyz1[0]-xyz2[0], xyz1[1]-xyz2[1], xyz1[2]-xyz2[2]);
514  }
515  }
516 
517 
519  art::Handle<std::vector<rb::Cluster> > slcs,
520  int& parentSlcIdx){
521  me::SlcME slcME;
523 
524  for (unsigned int i = 0; i < michel.size(); i++)
525  slcME.Add(michel[i]);
526  BestSlcMatch(slcME, vars, slcs, parentSlcIdx);
527  slcME.SetDistToSlc (vars.distToSlc);
528  slcME.SetParentSlc (vars.parentSlc);
529  slcME.SetDeltaT (slcME.MeanTNS() - vars.parentSlc->MeanTNS());
530  double mid = -10;
531  if (fParams.calibAvailable()){
532  mid = GetMID(slcME.DeltaT(),
533  slcME.DistToSlc(),
534  slcME.CalorimetricEnergy(),
535  slcME.NCell());
536  }
537  slcME.SetMID(mid);
538  return slcME;
539  }
540 
542  art::Handle<std::vector<rb::Cluster> > slcs,
543  art::Ptr<rb::Track> parentTrk){
544  me::TrkME trkME;
546 
547  for (unsigned int i = 0; i < slcME.NCell(); i++)
548  trkME.Add(slcME.Cell(i));
549 
550  trkME.SetDistToSlc (slcME.DistToSlc());
551  trkME.SetParentSlc (parentTrk);
552  trkME.SetDeltaT (slcME.MeanTNS() - parentTrk->MeanTNS());
553  double mid = -10;
554  if (fParams.calibAvailable()){
555  mid = GetMID(slcME.DeltaT(),
556  slcME.DistToSlc(),
557  slcME.CalorimetricEnergy(),
558  slcME.NCell());
559  }
560  trkME.SetMID(mid);
561  trkME.SetParentTrk(parentTrk);
562  return trkME;
563  }
564 
565  // Takes a ME Cluster, and finds out what physics slice it is most likely
566  // associated with using deltaT and distToSlc info
568  me::MEVars& vars,
569  art::Handle<std::vector<rb::Cluster> > slcs,
570  int& parentSlcIdx){
571  // Grab the MID building histogram
572  TH2F* hAssn = fLikeHists.fMEDeltaT_Dist;
573  // The assnScore is larger for slices that are appropriately near the ME
574  double assnScore = -1;
575  parentSlcIdx = 0;
576  double allSlcMinDist = 1.0e6;
577  for (unsigned int i = 0; i < slcs->size(); i++){
578  rb::Cluster slc = slcs->at(i);
579  if (slc.IsNoise() || slc.NCell() <= fParams.maxNCells())
580  continue;
581  art::Ptr<rb::Cluster> slcPtr(slcs, i);
582 
583  double deltaT = slcME.MeanTNS() - slc.MeanTNS();
584  double minDist = 10000;
585 
586  for (unsigned int j = 0; j < slcME.NCell(); j++){
587  double dist = MinHitSlcDist(slcME.Cell(j), slcs->at(i));
588  if (dist < minDist)
589  minDist = dist;
590  } // Loop over ME hits
591 
592  // Match by distance only
593  if ( fParams.matchMEtoSlcByDist() )
594  {
595  if ( minDist < allSlcMinDist )
596  {
597  parentSlcIdx = i;
598  allSlcMinDist = minDist;
599  vars.parentSlc = slcPtr;
600  vars.distToSlc = minDist;
601  }
602  }
603 
604  // Match by distance and time (standard method)
605  else
606  {
607  double curScore = hAssn->GetBinContent(hAssn->FindBin(deltaT, minDist));
608  if (curScore > assnScore){
609  parentSlcIdx = i;
610  assnScore = curScore;
611  vars.parentSlc = slcPtr;
612  vars.distToSlc = minDist;
613  }
614  }
615 
616  } // Loop over slices
617  }
618 
620  if (slcME.DistToSlc() > fParams.maxClosestApproach())
621  return false;
622  if (slcME.NCell() > fParams.maxNCells())
623  return false;
624  if (slcME.NCell() < fParams.minNCells())
625  return false;
626  if (slcME.DeltaT() > fParams.maxDeltaT())
627  return false;
628  if (slcME.DeltaT() < fParams.minDeltaT())
629  return false;
630  // Only apply calE presel if a calibration is available
631  if (fParams.calibAvailable()){
632  if (slcME.CalorimetricEnergy() > fParams.maxCalE())
633  return false;
634  if (slcME.CalorimetricEnergy() < fParams.minCalE())
635  return false;
636  }
637  return true;
638  }
639 
640  double MEFinder::GetMID(double DeltaT, double Dist, double CalE, int NCells){
641  TH2* meDeltaT_Dist = fLikeHists.fMEDeltaT_Dist;
642  TH2* meCalE_NCells = fLikeHists.fMECalE_NCells;
643 
644  int meDeltaT_Dist_Bin = meDeltaT_Dist->FindBin(DeltaT, Dist);
645  int meCalE_NCells_Bin = meCalE_NCells->FindBin(CalE, NCells);
646  double meLike = meDeltaT_Dist->GetBinContent(meDeltaT_Dist_Bin) *
647  meCalE_NCells->GetBinContent(meCalE_NCells_Bin);
648 
649  TH2* hfDeltaT_Dist = fLikeHists.fHFDeltaT_Dist;
650  TH2* hfCalE_NCells = fLikeHists.fHFCalE_NCells;
651 
652  int hfDeltaT_Dist_Bin = hfDeltaT_Dist->FindBin(DeltaT, Dist);
653  int hfCalE_NCells_Bin = hfCalE_NCells->FindBin(CalE, NCells);
654  double hfLike = hfDeltaT_Dist->GetBinContent(hfDeltaT_Dist_Bin) *
655  hfCalE_NCells->GetBinContent(hfCalE_NCells_Bin);
656 
657  if (meLike == 0 && hfLike == 0) return -10;
658  if (meLike == 0)
659  meLike = 1e-10;
660  if (hfLike == 0)
661  hfLike = 1e-10;
662 
663  return TMath::Log(meLike) - TMath::Log(hfLike);
664 
665  }
666 
668  const std::vector< art::Ptr<rb::Track> > trks,
669  art::Ptr<rb::Track>& parentTrk){
670  for (unsigned int i = 0; i < trks.size(); i++){
671  if (trks[i]->NXCell() == 0 || trks[i]->NYCell() == 0)
672  continue;
673  double dist = TrkMEDist(slcME, trks[i]->Stop());
674  if (dist < fParams.trkClosestApproach()){
675  parentTrk = trks[i];
676  return true;
677  }
678  }
679  return false;
680  }
681 
682  double MEFinder::TrkMEDist(me::SlcME me, TVector3 trkEnd){
683  double deltaX = ( (me.NXCell() == 0) ? 0 : me.MeanX() - trkEnd.X() );
684  double deltaY = ( (me.NYCell() == 0) ? 0 : me.MeanY() - trkEnd.Y() );
685  double deltaZ = me.MeanZ() - trkEnd.Z();
686  return util::pythag( deltaX, deltaY, deltaZ );
687  }
688 }
689 
690 namespace me
691 {
693 }
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:72
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
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
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
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
int evt
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
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)
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: structs.h:12
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
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
art::Ptr< rb::Cluster > parentSlc
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)
enum BeamMode string