Cluster.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Cluster.cxx
3 // \brief Cluster of CellHits. Base of other reconstruction classes.
4 // \version $Id: Cluster.cxx,v 1.31 2012-10-26 18:35:05 bckhouse Exp $
5 // \author Christopher Backhouse - bckhouse@caltech.edu
6 // \date $Date: 2012-10-26 18:35:05 $
7 ////////////////////////////////////////////////////////////////////////
8 #include <cfloat>
9 #include <iomanip>
10 
13 
14 #include "RecoBase/Cluster.h"
15 #include "Calibrator/Calibrator.h"
16 #include "Geometry/Geometry.h"
18 #include "RecoBase/CellHit.h"
19 #include "RecoBase/RecoHit.h"
20 #include "RecoBase/WeightedHit.h"
21 
22 namespace rb
23 {
24  //------------------------------------------------------------
26  : fView(geo::kXorY)
27  , fID(0)
28  , fNoiseCluster(false)
29  , fPrecalcTotalGeV(-1)
30  {
31  }
32 
33  //------------------------------------------------------------
35  : fView(geo::kXorY)
36  , fID(id)
37  , fNoiseCluster(false)
38  , fPrecalcTotalGeV(-1)
39  {
40  for(unsigned int i = 0; i < cells.size(); ++i) Add(cells[i]);
41  }
42 
43  //------------------------------------------------------------
45  : fView(view)
46  , fID(id)
47  , fNoiseCluster(false)
48  , fPrecalcTotalGeV(-1)
49  {
50  }
51 
52  //------------------------------------------------------------
54  : fView(view)
55  , fID(id)
56  , fNoiseCluster(false)
57  , fPrecalcTotalGeV(-1)
58  {
59  for(unsigned int i = 0; i < cells.size(); ++i) Add(cells[i]);
60 
61  assert((fView == geo::kX && fYCell.empty()) ||
62  (fView == geo::kY && fXCell.empty()));
63  }
64 
65  //------------------------------------------------------------
66  Cluster::Cluster(geo::View_t view, const std::vector<art::Ptr<CellHit> >& cells, int id)
67  : fView(view)
68  , fID(id)
69  , fNoiseCluster(false)
70  , fPrecalcTotalGeV(-1)
71  {
72  for(unsigned int i = 0; i < cells.size(); ++i) Add(cells[i]);
73 
74  assert((fView == geo::kX && fYCell.empty()) ||
75  (fView == geo::kY && fXCell.empty()));
76  }
77 
78  //------------------------------------------------------------
80  {
81  }
82 
83  //------------------------------------------------------------
85  {
86  fPrecalcTotalGeV = -1; // Invalidate any cached value
87 
88  assert(weight >= 0 && weight <= 1);
89 
90  // Are we optimizing (omitting) the weights storage?
91  const bool optWeights = (weight == 1 && fXWeights.empty() && fYWeights.empty());
92  if(!optWeights) EnsureWeightAlloc();
93 
94  if(cell->View() == geo::kX){
95  assert(fView != geo::kY);
96  fXCell.push_back(cell);
97  if(!optWeights) fXWeights.push_back(weight);
98  }
99  else{
100  assert(fView != geo::kX);
101  fYCell.push_back(cell);
102  if(!optWeights) fYWeights.push_back(weight);
103  }
104  }
105 
106  //------------------------------------------------------------
108  const std::vector<double>& weights)
109  {
110  const unsigned int N = cells.size();
111 
112  assert(weights.empty() || weights.size() == N);
113 
114  for(unsigned int n = 0; n < N; ++n){
115  assert(weights.empty() || (weights[n] >= 0 && weights[n] <= 1));
116  Add(cells[n], weights.empty() ? 1 : weights[n]);
117  }
118  }
119 
120  //------------------------------------------------------------
121  double Cluster::W(const rb::CellHit* hit) const
122  {
123  geo::View_t otherView = geo::kX;
124  if(hit->View() == geo::kX) otherView = geo::kY;
125 
126  // TB horizontal modules (Y-view) have readout on opposite end as ND/FD
128  if(geom->DetId() == novadaq::cnv::kTESTBEAM && hit->View() == geo::kY) return -MeanV(otherView);
129 
130  return MeanV(otherView);
131  }
132 
133  //------------------------------------------------------------
134  unsigned int Cluster::NCell(geo::View_t view) const
135  {
136  switch(view){
137  case geo::kX: return NXCell();
138  case geo::kY: return NYCell();
139  case geo::kXorY: return NCell();
140  default: assert(0 && "Unknown view");
141  }
142  }
143 
144  //------------------------------------------------------------
146  unsigned int viewIdx) const
147  {
148  switch(view){
149  case geo::kX: return XCell(viewIdx);
150  case geo::kY: return YCell(viewIdx);
151  case geo::kXorY: return Cell(viewIdx);
152  default: assert(0 && "Unknown view");
153  }
154  }
155 
156  //------------------------------------------------------------
157  art::Ptr<rb::CellHit> Cluster::XCell(unsigned int xIdx) const
158  {
159  assert(xIdx < fXCell.size());
160 
161  return fXCell[xIdx];
162  }
163 
164  //------------------------------------------------------------
165  art::Ptr<rb::CellHit> Cluster::YCell(unsigned int yIdx) const
166  {
167  assert(yIdx < fYCell.size());
168 
169  return fYCell[yIdx];
170  }
171 
172  //------------------------------------------------------------
173  art::Ptr<rb::CellHit> Cluster::Cell(unsigned int globalIdx) const
174  {
175  if(globalIdx < NXCell()) return XCell(globalIdx);
176  return YCell(globalIdx-NXCell());
177  }
178 
179  //------------------------------------------------------------
181  {
183  for(size_t c = 0; c < fXCell.size(); ++c) all.push_back(fXCell[c]);
184  for(size_t c = 0; c < fYCell.size(); ++c) all.push_back(fYCell[c]);
185 
186  return all;
187  }
188 
189  //------------------------------------------------------------
190  std::vector<geo::OfflineChan> Cluster::OfflineChans() const
191  {
192  std::vector<geo::OfflineChan> ret;
193  ret.reserve(NCell());
194  for(unsigned int i = 0; i < NCell(); ++i)
195  ret.push_back(Cell(i)->OfflineChan());
196  return ret;
197  }
198 
199  //------------------------------------------------------------
200  std::vector<rb::WeightedHit> Cluster::WeightedHits() const
201  {
202  std::vector<rb::WeightedHit> all;
203  for(unsigned int cellIdx = 0; cellIdx < NCell(); ++cellIdx)
204  all.push_back(rb::WeightedHit(Cell(cellIdx), Weight(cellIdx)));
205  return all;
206  }
207 
208  //------------------------------------------------------------
209  double Cluster::Weight(unsigned int globalIdx) const
210  {
211  if(fXWeights.empty() && fYWeights.empty()) return 1;
212 
213  if(globalIdx < NXCell()) return Weight(geo::kX, globalIdx);
214  return Weight(geo::kY, globalIdx-NXCell());
215  }
216 
217  //------------------------------------------------------------
218  double Cluster::Weight(geo::View_t view, unsigned int viewIdx) const
219  {
220  if(fXWeights.empty() && fYWeights.empty()) return 1;
221 
222  if(view == geo::kX){
223  assert(viewIdx < fXWeights.size());
224  return fXWeights[viewIdx];
225  }
226  else{
227  assert(viewIdx < fYWeights.size());
228  return fYWeights[viewIdx];
229  }
230  }
231 
232  //------------------------------------------------------------
234  {
236  ret.SetNoise(fNoiseCluster);
237 
238  const unsigned int I = NCell();
239  const unsigned int J = excl->NCell();
240  for(unsigned int i = 0; i < I; ++i){
241  const art::Ptr<rb::CellHit> chit = Cell(i);
242  bool clash = false;
243  for(unsigned int j = 0; j < J; ++j){
244  if(excl->Cell(j) == chit){
245  clash = true;
246  break;
247  }
248  } // end for j
249  if(!clash) ret.Add(chit, Weight(i));
250  }
251 
252  if(ret.NCell() == 0)
253  mf::LogWarning("Cluster") << "Exclude created empty cluster";
254 
255  return ret;
256  }
257 
258  //------------------------------------------------------------
260  {
262  rb::RecoHit rb(cal->MakeRecoHit(*chit, this->W(chit.get())));
263  return rb;
264  }
265 
266  //------------------------------------------------------------
267  rb::RecoHit Cluster::RecoHit(geo::View_t view, unsigned int viewIdx) const
268  {
269  return RecoHit(Cell(view, viewIdx));
270  }
271 
272  //------------------------------------------------------------
273  rb::RecoHit Cluster::RecoHit(unsigned int globalIdx) const
274  {
275  return RecoHit(Cell(globalIdx));
276  }
277 
278  //------------------------------------------------------------
280  {
281  fPrecalcTotalGeV = -1; // Invalidate any cached value
282 
283  fXCell.clear();
284  fYCell.clear();
285  fXWeights.clear();
286  fYWeights.clear();
287  }
288 
289  //-----------------------------------------------------------
291  {
292  fPrecalcTotalGeV = -1; // Invalidate any cached value
293 
294  if(hit->View() == geo::kX){
295  // try to find the hit in the xcells
296  for(unsigned int i = 0;i<NXCell();++i){
297  if(fXCell[i] == hit){
299  if(!fXWeights.empty()) fXWeights.erase(fXWeights.begin()+i);
300  break;
301  }
302  }
303  }
304  else if(hit->View() == geo::kY){
305  // try to find the hit in the ycells
306  for(unsigned int i = 0;i<NYCell();++i){
307  if(fYCell[i] == hit){
309  if(!fYWeights.empty()) fYWeights.erase(fYWeights.begin()+i);
310  break;
311  }
312  }
313  }
314  }
315 
316  //-----------------------------------------------------------
317  void Cluster::RemoveHit(unsigned int globalIdx)
318  {
319  fPrecalcTotalGeV = -1; // Invalidate any cached value
320 
321  assert(globalIdx < NCell());
322  if(globalIdx < NXCell()) fXCell.erase(fXCell.begin()+globalIdx);
323  else fYCell.erase(fYCell.begin()+globalIdx-NXCell());
324  }
325 
326  //------------------------------------------------------------
327  void Cluster::SetWeight(unsigned int globalIdx, double weight)
328  {
329  fPrecalcTotalGeV = -1; // Invalidate any cached value
330 
331  if(globalIdx < NXCell())
332  SetWeight(geo::kX, globalIdx, weight);
333  else
334  SetWeight(geo::kY, globalIdx-NXCell(), weight);
335  }
336 
337  //------------------------------------------------------------
338  void Cluster::SetWeight(geo::View_t view, unsigned int viewIdx,
339  double weight)
340  {
341  fPrecalcTotalGeV = -1; // Invalidate any cached value
342 
343  assert(weight >= 0 && weight <= 1);
344 
345  // Are we optimizing (omitting) the weights storage?
346  const bool optWeights = (weight == 1 && fXWeights.empty() && fYWeights.empty());
347  if(!optWeights) EnsureWeightAlloc();
348 
349  if(view == geo::kX){
350  assert(viewIdx < fXCell.size());
351  if(!optWeights) fXWeights[viewIdx] = weight;
352  }
353  else{
354  assert(viewIdx < fYCell.size());
355  if(!optWeights) fYWeights[viewIdx] = weight;
356  }
357  }
358 
359  //------------------------------------------------------------
360  double Cluster::TotalADC() const
361  {
362  double ret = 0;
363  for(unsigned int i = 0; i < NCell(); ++i)
364  ret += Cell(i)->ADC()*Weight(i);
365  return ret;
366  }
367 
368  //------------------------------------------------------------
369  double Cluster::TotalPE() const
370  {
371  double ret = 0;
372  for(unsigned int i = 0; i < NCell(); ++i)
373  ret += Cell(i)->PE()*Weight(i);
374  return ret;
375  }
376 
377  //------------------------------------------------------------
378  double Cluster::TotalGeV(EEnergyCalcScheme escheme) const
379  {
380  if(escheme == kUsePrecalcEnergy){
381  if(fPrecalcTotalGeV >= 0) return fPrecalcTotalGeV;
382 
383  mf::LogWarning("rb::Cluster::TotalGeV()") << "Passed kUsePrecalcEnergy but no energy was cached in the Cluster. Falling back to recomputing." << std::endl;
384  }
385 
386  // Optimization that's only safe for literal Clusters
387  if(typeid(*this) == typeid(rb::Cluster)) return TotalGeVFastClusterOnly();
388 
389  // For derived classes, fall through to the generic implementation
390 
391  double ret = 0;
392  for(unsigned int i = 0; i < NCell(); ++i){
393  const rb::RecoHit rhit = RecoHit(i);
394  if(rhit.IsCalibrated()) ret += rhit.GeV()*Weight(i);
395  }
396  return ret;
397  }
398 
399  //------------------------------------------------------------
401  {
402  // This is only safe for objects that return the same W for all hits in the
403  // same view. We only know that's the case for Cluster. It's definitely not
404  // correct for Prong and Track.
405  assert(typeid(*this) == typeid(rb::Cluster));
406 
408 
409  double gev = 0;
410 
411  // We just use the mean position for the W estimate in all cases.
412  // Precalculate once up front, instead of NCell() times as we would
413  // otherwise.
414  const TVector3 mean = MeanXYZ();
415 
416  // X and Y cells have the mean of the opposite view as their W
417  for(unsigned int i = 0; i < NXCell(); ++i){
418  const rb::RecoHit rhit(cal->MakeRecoHit(*(XCell(i)), mean.Y()));
419  if(rhit.IsCalibrated()) gev += rhit.GeV()*Weight(geo::kX, i);
420  }
421  for(unsigned int i = 0; i < NYCell(); ++i){
422  const rb::RecoHit rhit(cal->MakeRecoHit(*(YCell(i)), mean.X()));
423  if(rhit.IsCalibrated()) gev += rhit.GeV()*Weight(geo::kY, i);
424  }
425 
426  return gev;
427  }
428 
429  //------------------------------------------------------------
430  double Cluster::TotalWeight() const
431  {
432  double ret = 0;
433  for(unsigned int cellIdx = 0; cellIdx < NCell(); ++cellIdx)
434  ret += Weight(cellIdx);
435  return ret;
436  }
437 
438  //------------------------------------------------------------
440  {
442  return cal->GetGeVToCalorimetricScale()*TotalGeV(escheme);
443  }
444 
445  //------------------------------------------------------------
446  TVector3 Cluster::MinXYZ() const
447  {
448  TVector3 lo, junk;
449  MinMaxMeanXYZ(lo, junk, junk);
450  return lo;
451  }
452 
453  //------------------------------------------------------------
455  {
456  assert(view == geo::kX || view == geo::kY);
457  if(view == geo::kX) return MinXYZ().X();
458  return MinXYZ().Y();
459  }
460 
461  //------------------------------------------------------------
462  unsigned int Cluster::MinPlane(geo::View_t view) const
463  {
464  assert(NCell(view) > 0);
465  unsigned int ret = UINT_MAX;
466  for(unsigned int i = 0; i < NCell(view); ++i)
467  if(Cell(view, i)->Plane() < ret) ret = Cell(view, i)->Plane();
468  return ret;
469  }
470 
471  //------------------------------------------------------------
472  unsigned int Cluster::MinCell(geo::View_t view) const
473  {
474  assert(NCell(view) > 0);
475  unsigned int ret = UINT_MAX;
476  for(unsigned int i = 0; i < NCell(view); ++i)
477  if(Cell(view, i)->Cell() < ret) ret = Cell(view, i)->Cell();
478  return ret;
479  }
480 
481  //------------------------------------------------------------
482  double Cluster::MinTNS() const
483  {
484  assert(NCell() > 0);
485  double ret = DBL_MAX;
486  for(unsigned int i = 0; i < NCell(); ++i)
487  if(Cell(i)->TNS() < ret) ret = Cell(i)->TNS();
488  return ret;
489  }
490 
491  //------------------------------------------------------------
492  TVector3 Cluster::MaxXYZ() const
493  {
494  TVector3 hi, junk;
495  MinMaxMeanXYZ(junk, hi, junk);
496  return hi;
497  }
498 
499  //------------------------------------------------------------
501  {
502  assert(view == geo::kX || view == geo::kY);
503  if(view == geo::kX) return MaxXYZ().X();
504  return MaxXYZ().Y();
505  }
506 
507  //------------------------------------------------------------
508  unsigned int Cluster::MaxPlane(geo::View_t view) const
509  {
510  assert(NCell(view) > 0);
511  unsigned int ret = 0;
512  for(unsigned int i = 0; i < NCell(view); ++i)
513  if(Cell(view, i)->Plane() > ret) ret = Cell(view, i)->Plane();
514  return ret;
515  }
516 
517  //------------------------------------------------------------
518  unsigned int Cluster::MaxCell(geo::View_t view) const
519  {
520  assert(NCell(view) > 0);
521  unsigned int ret = 0;
522  for(unsigned int i = 0; i < NCell(view); ++i)
523  if(Cell(view, i)->Cell() > ret) ret = Cell(view, i)->Cell();
524  return ret;
525  }
526 
527  //------------------------------------------------------------
528  double Cluster::MaxTNS() const
529  {
530  assert(NCell() > 0);
531  double ret = -DBL_MAX;
532  for(unsigned int i = 0; i < NCell(); ++i)
533  if(Cell(i)->TNS() > ret) ret = Cell(i)->TNS();
534  return ret;
535  }
536 
537  //------------------------------------------------------------
538  TVector3 Cluster::MeanXYZ(rb::AveragingScheme scheme) const
539  {
540  TVector3 mean, junk;
541  MinMaxMeanXYZ(junk, junk, mean, scheme);
542  return mean;
543  }
544 
545  //------------------------------------------------------------
547  {
548  assert(view == geo::kX || view == geo::kY);
549  if(view == geo::kX) return MeanXYZ(scheme).X();
550  return MeanXYZ(scheme).Y();
551  }
552 
553  //------------------------------------------------------------
555  {
556  assert(NCell() > 0);
557  double ret = 0;
558  double totW = 0;
559  for(unsigned int i = 0; i < NCell(); ++i){
560  double w = Weight(i);
561  if(scheme == kByEnergy) w *= Cell(i)->PE();
562 
563  ret += Cell(i)->TNS()*w;
564  totW += w;
565  }
566  return ret/totW;
567  }
568 
569  //------------------------------------------------------------
571  {
572  return MaxCell(view)-MinCell(view)+1;
573  }
574 
575  //------------------------------------------------------------
576  bool Cluster::operator<(const Cluster& other) const
577  {
578  return this->MeanZ() < other.MeanZ();
579  }
580 
581  //------------------------------------------------------------
582  void Cluster::MinMaxMeanXYZ(TVector3& lo, TVector3& hi, TVector3& mean,
583  rb::AveragingScheme scheme) const
584  {
585  // If this becomes a bottleneck we can cache the results
586 
587  /// \todo Won't be quite right in the FD when cells aren't at right angles?
588 
589  assert(NCell() > 0);
590 
592 
593  double Wx = 0;
594  double Wy = 0;
595  for(int n = 0; n < 3; ++n){
596  lo(n) = DBL_MAX;
597  hi(n) = -DBL_MAX;
598  mean(n) = 0;
599  }
600 
601  for(unsigned int i = 0; i < NCell(); ++i){
602  const art::Ptr<rb::CellHit>& h = Cell(i);
603  double xyz[3];
604  geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
605  const double x = xyz[0];
606  const double y = xyz[1];
607  const double z = xyz[2];
608  const geo::View_t view = geom->Plane(h->Plane())->View();
609  double weight = Weight(i);
610  if(scheme == kByEnergy) weight *= h->PE();
611 
612  if(view == geo::kX){
613  if(x < lo.X()) lo.SetX(x);
614  if(x > hi.X()) hi.SetX(x);
615  mean.SetX(mean.X()+x*weight);
616  Wx += weight;
617  }
618  if(view == geo::kY){
619  if(y < lo.Y()) lo.SetY(y);
620  if(y > hi.Y()) hi.SetY(y);
621  mean.SetY(mean.Y()+y*weight);
622  Wy += weight;
623  }
624  if(z < lo.Z()) lo.SetZ(z);
625  if(z > hi.Z()) hi.SetZ(z);
626  mean.SetZ(mean.Z()+z*weight);
627  }
628 
629  if(Wx) mean.SetX(mean.X()/Wx);
630  if(Wy) mean.SetY(mean.Y()/Wy);
631  if(Wx || Wy) mean.SetZ(mean.Z()/(Wx+Wy));
632  }
633 
634  //------------------------------------------------------------
636  {
637  std::set<int> planes;
638 
639  for(unsigned int i = 0; i < NCell(view); ++i)
640  planes.insert(Cell(view, i)->Plane());
641 
642  if(planes.empty()) return 0;
643 
644  // In 3D expect a hit every plane, in 2D only every other plane.
645  //
646  // This isn't right if the planes aren't alternating in views, but
647  // thankfully we stopped building detectors like that.
648  const int expGap = (view == geo::kXorY) ? 1 : 2;
649 
650  int prevPlane = -999;
651  int curRun = 1;
652  int longestRun = 0;
653  for(auto plane: planes){
654  if(plane - prevPlane == expGap){
655  ++curRun;
656  longestRun = std::max(longestRun, curRun);
657  }
658  else{
659  curRun = 1;
660  }
661  prevPlane = plane;
662  } // end for plane
663 
664  return longestRun;
665  }
666 
667  //------------------------------------------------------------
669  {
670  std::set<int> planes;
671 
672  for(unsigned int i = 0; i < NCell(view); ++i)
673  planes.insert(Cell(view, i)->Plane());
674 
675  if(planes.empty()) return 0;
676 
677  int maxDiff = 1;
678  int prevPlane = *planes.begin();
679  for(int plane: planes){
680  maxDiff = std::max(maxDiff, plane-prevPlane);
681  prevPlane = plane;
682  }
683 
684  // One plane difference as expected in 3D isn't a gap
685  if(view == geo::kXorY) return maxDiff-1;
686 
687  // Two plane difference in 2D isn't a gap. Beyond that, every two planes
688  // difference only indicates one missing plane of hits.
689  return (maxDiff-2)/2;
690  }
691 
692  //------------------------------------------------------------
694  {
695  if(view == geo::kXorY) view = View();
696 
697  const int stride = (view == geo::kXorY) ? 1 : 2;
698 
699  std::set<int> planes;
700 
701  for(unsigned int i = 0; i < NCell(view); ++i)
702  planes.insert(Cell(view, i)->Plane());
703 
704  if(planes.empty()) return 0;
705 
706  return (*planes.rbegin() - *planes.begin())/stride - (int)planes.size() + 1;
707 
708  }
709 
710 
711  //------------------------------------------------------------
713  {
714  if(fXWeights.size() != fXCell.size()) fXWeights.resize(fXCell.size(), 1);
715  if(fYWeights.size() != fYCell.size()) fYWeights.resize(fYCell.size(), 1);
716  }
717 
718  //------------------------------------------------------------
720  const art::Ptr<rb::CellHit>& b)
721  {
722  if (a->fTDC[0]!=b->fTDC[0]) return a->fTDC[0]<b->fTDC[0];
723  if (a->Plane()!=b->Plane()) return a->Plane()<b->Plane();
724  if (a->Cell() !=b->Cell()) return a->Cell() <b->Cell();
725  if (a->fADC[0]!=b->fADC[0]) return a->fADC[0]<b->fADC[0];
726  if (a->ID() !=b->ID()) return a->ID() <b->ID();
727  return false; // Should never get here...
728  }
729 
730  //------------------------------------------------------------
732  {
733  //
734  // Unweighted clusters are an easy sorting job
735  //
736  if (fXWeights.size()==0 && fYWeights.size()==0) {
737  sort(fXCell.begin(), fXCell.end(), standard_compare);
738  sort(fYCell.begin(), fYCell.end(), standard_compare);
739  return;
740  }
741  //
742  // \todo Sorting weighted hits requires more complicated code
743  //
744  std::cerr << __FILE__ << ":" << __LINE__
745  << " Sort for weighted hits not implemented"
746  << std::endl;
747  abort();
748  }
749 
750  //------------------------------------------------------------
752  {
753  const double gev = TotalGeV(kRecomputeEnergy);
754  if(fPrecalcTotalGeV < 0 || savemode == kResetTotalGeV){
755  fPrecalcTotalGeV = gev;
756  }
757  else{
758  if(fabs(gev-fPrecalcTotalGeV) > .01){
759  mf::LogError("rb::Cluster") << "SaveTotalGeV() already called. Old value was " << fPrecalcTotalGeV << ", now " << gev
760  << ". If you intended to reset the stored value, you must pass kResetTotalGeV. Aborting." << std::endl;
761  abort();
762  }
763  else{
764  mf::LogWarning("rb::Cluster") << "Called SaveTotalGeV() twice for the same cluster, which is inefficient." << std::endl;
765  }
766  }
767  }
768 
769  //------------------------------------------------------------
770  std::ostream& operator<< (std::ostream& o, const Cluster& c)
771  {
772  o << std::setiosflags(std::ios::fixed) << std::setprecision(2);
773  o << "View = " << std::setw(1) << std::right << c.View()
774  << " ID = " << std::setw(1) << std::right << c.ID()
775  << " #Cells = " << std::setw(3) << std::right << c.NCell()
776  << " Is3D = " << std::setw(1) << std::right << c.Is3D()
777  << " TotalADC = " << std::setw(5) << std::right << c.TotalADC();
778 
779  return o;
780  }
781 
782 } // end namespace
783 ////////////////////////////////////////////////////////////////////////
double MinV(geo::View_t view) const
Definition: Cluster.cxx:454
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
float TNS() const
Definition: CellHit.h:46
T max(const caf::Proxy< T > &a, T b)
double TotalPE() const
Sum of the PE value of all the contained hits.
Definition: Cluster.cxx:369
std::vector< geo::OfflineChan > OfflineChans() const
Positions of all the CellHits.
Definition: Cluster.cxx:190
TSpline3 lo("lo", xlo, ylo, 12,"0")
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
void StandardSort()
Put the cells in the cluster into a standard order.
Definition: Cluster.cxx:731
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
double fPrecalcTotalGeV
-1 = uninitialized
Definition: Cluster.h:324
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void SetNoise(bool noise)
Declare the cluster to consist of noise hits or not.
Definition: Cluster.h:161
const Var weight
void SavePrecalcTotalGeV(ESaveGeVMode savemode)
Store the current result of TotalGeV / CalorimetricEnergy.
Definition: Cluster.cxx:751
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
iterator begin()
Definition: PtrVector.h:223
Definition: geo.h:1
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double MaxV(geo::View_t view) const
Definition: Cluster.cxx:500
iterator erase(iterator position)
Definition: PtrVector.h:508
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
double GetGeVToCalorimetricScale() const
For use by RecoBase classes.
A collection of associated CellHits.
Definition: Cluster.h:47
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:233
double TotalGeVFastClusterOnly() const
Optimized form of TotalGeV, only valid for actual Clusters.
Definition: Cluster.cxx:400
OStream cerr
Definition: OStream.cxx:7
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
void RemoveHit(const art::Ptr< rb::CellHit > hit)
Remove hit from current cluster.
Definition: Cluster.cxx:290
double MeanV(geo::View_t view, rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:546
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void SetWeight(unsigned int globalIdx, double weight)
Set weight of the cell at this index.
Definition: Cluster.cxx:327
TVector3 MeanXYZ(rb::AveragingScheme=kDefaultScheme) const
Definition: Cluster.cxx:538
rb::Cluster Exclude(const rb::Cluster *excl) const
Create a cluster from this one, but with the hits of excl removed.
Definition: Cluster.cxx:233
art::PtrVector< rb::CellHit > fYCell
collection of y-view cells in cluster
Definition: Cluster.h:318
bool standard_compare(const art::Ptr< rb::CellHit > &a, const art::Ptr< rb::CellHit > &b)
Definition: Cluster.cxx:719
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
unsigned int ExtentCell(geo::View_t view) const
Definition: Cluster.cxx:570
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
unsigned short Cell() const
Definition: CellHit.h:40
Use a value computed by a previous call to SavePrecalcTotalGeV.
Definition: Cluster.h:179
double MinTNS() const
Definition: Cluster.cxx:482
art::PtrVector< rb::CellHit > fXCell
collection of x-view cells in cluster
Definition: Cluster.h:317
Altering the already-set value is not an error.
Definition: Cluster.h:283
AveragingScheme
Definition: Cluster.h:27
int MostMissingPlanes(geo::View_t view) const
Longest run of adjacent planes with no hits.
Definition: Cluster.cxx:668
double TotalADC() const
Sum of the ADC of all the contained hits.
Definition: Cluster.cxx:360
TSpline3 hi("hi", xhi, yhi, 18,"0")
bool operator<(const Cluster &other) const
Definition: Cluster.cxx:576
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
unsigned int NCell() const
Number of cells in either view.
Definition: Cluster.h:110
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
const double a
std::vector< rb::WeightedHit > WeightedHits() const
Get all hits from both views, with weights attached.
Definition: Cluster.cxx:200
geo::View_t fView
view this cluster is in
Definition: Cluster.h:316
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
iterator end()
Definition: PtrVector.h:237
float PE() const
Definition: CellHit.h:42
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
bool fNoiseCluster
flag for whether this is a noise cluster
Definition: Cluster.h:322
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
Var weights
const double j
Definition: BetheBloch.cxx:29
std::vector< int32_t > fTDC
TDC(-like) time value. Event time is subtracted. Vector structure is a historical artifact...
Definition: RawDigit.h:50
Perform a "2 point" Hough transform on a collection of hits.
z
Definition: test.py:28
bool empty() const
Definition: PtrVector.h:336
double MaxTNS() const
Definition: Cluster.cxx:528
size_type size() const
Definition: PtrVector.h:308
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
int NMissingPlanes(geo::View_t view) const
Total number of missing planes in cluster.
Definition: Cluster.cxx:693
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
std::vector< int16_t > fADC
list of ADC(-like) charge values
Definition: RawDigit.h:49
virtual bool Is3D() const
Definition: Cluster.h:95
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
std::vector< double > fXWeights
Weights, matching cell indexing.
Definition: Cluster.h:319
const int ID() const
Definition: CellHit.h:47
float GeV() const
Definition: RecoHit.cxx:69
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
T const * get() const
Definition: Ptr.h:321
TVector3 MinXYZ() const
Definition: Cluster.cxx:446
Definition: structs.h:12
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
int MostContiguousPlanes(geo::View_t view) const
Longest run of adjacent planes with hits.
Definition: Cluster.cxx:635
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
double TotalWeight() const
Sum of all the weights. The effective number of hits.
Definition: Cluster.cxx:430
int fID
ID for cluster.
Definition: Cluster.h:321
const hit & b
Definition: hits.cxx:21
Simple hit+weight pair, returned from rb::Cluster::WeightedHits.
Definition: WeightedHit.h:18
assert(nhit_max >=nhit_nbins)
virtual ~Cluster()
Definition: Cluster.cxx:79
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
const int ID() const
Definition: Cluster.h:75
void MinMaxMeanXYZ(TVector3 &lo, TVector3 &hi, TVector3 &mean, rb::AveragingScheme scheme=kDefaultScheme) const
Gets the min/max/mean all at once, called by the functions above.
Definition: Cluster.cxx:582
void EnsureWeightAlloc()
Helper. Resizes weights vectors to match cell vectors.
Definition: Cluster.cxx:712
Helper for AttenCurve.
Definition: Path.h:10
virtual void Clear()
Forget about all owned cell hits.
Definition: Cluster.cxx:279
void clear()
Definition: PtrVector.h:537
Float_t w
Definition: plot.C:20
In addition to the intrinsic weights.
Definition: Cluster.h:30
std::vector< double > fYWeights
May be empty, means all weights are 1.
Definition: Cluster.h:320
Definition: fwd.h:28
EEnergyCalcScheme
Definition: Cluster.h:177
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
Encapsulate the geometry of one entire detector (near, far, ndos)
Default. Ask Calibrator about each hit, and sum.
Definition: Cluster.h:178
friend std::ostream & operator<<(std::ostream &o, const Cluster &c)
Definition: Cluster.cxx:770