Cluster.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Cluster.h
3 // \brief Cluster of CellHits. Base of other reconstruction classes.
4 // \version $Id: Cluster.h,v 1.22 2012-10-26 22:28:40 bckhouse Exp $
5 // \author Author: Christopher Backhouse - bckhouse@caltech.edu
6 // \date $Date: 2012-10-26 22:28:40 $
7 ////////////////////////////////////////////////////////////////////////
8 #ifndef RBCLUSTER_H
9 #define RBCLUSTER_H
10 
11 #include <iosfwd>
12 #include <vector>
13 
15 
18 
19 #include "TVector3.h"
20 
21 namespace rb
22 {
23  class CellHit;
24  class RecoHit;
25  struct WeightedHit;
26 
28  {
30  kByEnergy ///< In addition to the intrinsic weights
31  };
33 
34  /// \brief A collection of associated CellHits
35  ///
36  /// A Cluster is simply a "bag" of CellHits which have been reconstructed to
37  /// be associated with each other in some way.
38  ///
39  /// A cluster may be created as either 2D or 3D (ie containing hits from one or
40  /// both views).
41  ///
42  /// This class also contains various utility functions for obtaining the bulk
43  /// properties of these hits.
44  ///
45  /// It forms the basis of the reconstruction inheritance hierarchy, with the
46  /// derived classes adding additional information about the hits.
47  class Cluster
48  {
49  public:
50  /// \name 3D cluster constructors
51  //@{
52  Cluster();
53 #ifndef __GCCXML__
54  explicit Cluster(const art::PtrVector<rb::CellHit>& cells, int id=0);
55  //@}
56  /// \name 2D cluster constructors
57  //@{
58  Cluster(geo::View_t view, int id=0);
59  Cluster(geo::View_t view, const art::PtrVector<rb::CellHit>& cells, int id=0);
60  Cluster(geo::View_t view, const std::vector<art::Ptr<rb::CellHit> >& cells, int id=0);
61 #endif // __GCCXML__
62  //@}
63 
64  virtual ~Cluster();
65 
66 #ifndef __GCCXML__
67 
68  // Note: the redundant 'rb::' is not required by C++ but by
69  // CINT. Best to just keep it until CINT catches up.
70  virtual void Add(const art::Ptr<rb::CellHit>& cell, double weight = 1);
71  virtual void Add(const art::PtrVector<rb::CellHit>& cells,
72  const std::vector<double>& weights = std::vector<double>());
73 
74  void SetID(int id) { fID = id; }
75  const int ID() const { return fID; }
76 
77  /// \brief Estimate the unmeasured coordinate of \a chit
78  ///
79  /// The nomenclature for the coordinates is that XYZ are the physical
80  /// coordinates. Whichever of XY is measured by the cell is "V", the
81  /// direction along the cell is "W". The origin of VW space is the same as
82  /// the origin of XY space, the detector centreline.
83  ///
84  /// This method should be overridden in all descendants of \ref Cluster
85  /// which have additional information to make this determination.
86  ///
87  /// For \ref Cluster, this is simply the average coordinate of all hits in
88  /// the view in question.
89  ///
90  /// \param chit Any \ref CellHit, not necessarily owned by this container
91  /// \return Estimated unmeasured coordinate of this hit (cm) assuming
92  /// it to be associated with this reconstructed object.
93  virtual double W(const rb::CellHit* chit) const;
94 
95  virtual bool Is3D() const {return ((fXCell.size() > 0) && (fYCell.size() > 0));}
96  bool Is2D() const {return !Is3D();}
97 
98  /// kXorY for 3D clusters.
99  virtual geo::View_t View() const {return fView;}
100 
101  /// \name Cell counts
102  //@{
103  /// Number of cells in view \a view
104  unsigned int NCell(geo::View_t view) const;
105  /// Number of cells in the x-view
106  unsigned int NXCell() const {return fXCell.size();}
107  /// Number of cells in the y-view
108  unsigned int NYCell() const {return fYCell.size();}
109  /// Number of cells in either view
110  unsigned int NCell() const {return NXCell()+NYCell();}
111  //@}
112 
113  /// \name Cell Accessors
114  //@{
115  /// Get the ith cell from view \a view
116  art::Ptr<rb::CellHit> Cell(geo::View_t view, unsigned int viewIdx) const;
117  /// Get the ith cell in the x-view
118  art::Ptr<rb::CellHit> XCell(unsigned int xIdx) const;
119  /// Get the ith cell in the y-view
120  art::Ptr<rb::CellHit> YCell(unsigned int yIdx) const;
121  /// Get the ith cell from either view
122  art::Ptr<rb::CellHit> Cell(unsigned int globalIdx) const;
123  /// Get all cells from the x-view
124  const art::PtrVector<rb::CellHit>& XCells() const { return fXCell; }
125  /// Get all cells from the x-view
126  const art::PtrVector<rb::CellHit>& YCells() const { return fYCell; }
127  /// Get all cells from both views
129  /// Positions of all the CellHits
130  std::vector<geo::OfflineChan> OfflineChans() const;
131  /// Get all hits from both views, with weights attached
132  std::vector<rb::WeightedHit> WeightedHits() const;
133  /// Weight assigned to the cell
134  double Weight(unsigned int globalIdx) const;
135  /// Weight assigned to the cell
136  double Weight(geo::View_t view, unsigned int viewIdx) const;
137  //@}
138 
139  /// Create a cluster from this one, but with the hits of \a excl removed
140  rb::Cluster Exclude(const rb::Cluster* excl) const;
141 
142  /// \brief Return calibrated hit based on assumed W coordinate
143  ///
144  /// To calibrate a hit you must know its distance from the readout. This is
145  /// a convenience function which uses the \ref W function to estimate
146  /// the position of a \ref CellHit and then uses \ref calib::Calibrator to
147  /// produce a calibrated hit from that information.
148  rb::RecoHit RecoHit(const art::Ptr<rb::CellHit>& chit) const;
149  rb::RecoHit RecoHit(geo::View_t view, unsigned int viewIdx) const;
150  rb::RecoHit RecoHit(unsigned int globalIdx) const;
151 
152  /// \brief Forget about all owned cell hits
153  ///
154  /// Virtual function because derived classes may have more stuff to clear
155  virtual void Clear();
156  /// Remove \a hit from current cluster
157  void RemoveHit(const art::Ptr<rb::CellHit> hit);
158  /// Remove the ith hit from current cluster
159  void RemoveHit(unsigned int globalIdx);
160  /// Declare the cluster to consist of noise hits or not
161  void SetNoise(bool noise) {fNoiseCluster = noise;}
162  /// Is the noise flag set?
163  bool IsNoise() const {return fNoiseCluster;}
164 
165  /// Set weight of the cell at this index
166  void SetWeight(unsigned int globalIdx, double weight);
167  /// Set weight of the cell at this index in this view
168  void SetWeight(geo::View_t view, unsigned int viewIdx, double weight);
169 
170  /// \name Totals
171  //@{
172  /// Sum of the ADC of all the contained hits
173  double TotalADC() const;
174  /// Sum of the PE value of all the contained hits
175  double TotalPE() const;
176 
178  kRecomputeEnergy = 0, ///< Default. Ask Calibrator about each hit, and sum
179  kUsePrecalcEnergy = 1 ///< Use a value computed by a previous call to \ref SavePrecalcTotalGeV
180  };
181 
182  /// \brief Simple sum of the estimated GeV of all the hits
183  ///
184  /// WARNING: This is not a sensible energy estimator by itself. If you just
185  /// want a simple calorimetric estimate look at \ref
186  /// CalorimetricEnergy. Otherwise, a range of more sophisticated/elaborate
187  /// estimators are available elsewhere.
188  double TotalGeV(EEnergyCalcScheme escheme = kRecomputeEnergy) const;
189  /// Sum of all the weights. The effective number of hits
190  double TotalWeight() const;
191  /// \brief Simple estimate of neutrino energy
192  ///
193  /// Uses a simple scale factor to correct for dead material and thresholds,
194  /// could be more advanced, but this is simple enough to use for
195  /// preselection etc.
196  double CalorimetricEnergy(EEnergyCalcScheme escheme = kRecomputeEnergy) const;
197  //@}
198 
199  /// \name Minima
200  /// The "lower" corner of a box containing all the hits
201  //@{
202  TVector3 MinXYZ() const;
203  double MinV(geo::View_t view) const;
204  double MinX() const {return MinXYZ().X();}
205  double MinY() const {return MinXYZ().Y();}
206  double MinZ() const {return MinXYZ().Z();}
207  double MinTNS() const;
208  unsigned int MinPlane(geo::View_t view = geo::kXorY) const;
209  unsigned int MinCell(geo::View_t view) const;
210  //@}
211 
212  /// \name Maxima
213  /// The "upper" corner of a box containing all the hits
214  //@{
215  TVector3 MaxXYZ() const;
216  double MaxV(geo::View_t view) const;
217  double MaxX() const {return MaxXYZ().X();}
218  double MaxY() const {return MaxXYZ().Y();}
219  double MaxZ() const {return MaxXYZ().Z();}
220  double MaxTNS() const;
221  unsigned int MaxPlane(geo::View_t view = geo::kXorY) const;
222  unsigned int MaxCell(geo::View_t view) const;
223  //@}
224 
225  /// \name Means
226  /// The unweighted mean position and time of all the hits in the cluster
227  //@{
228  TVector3 MeanXYZ(rb::AveragingScheme = kDefaultScheme) const;
229  double MeanV(geo::View_t view,
230  rb::AveragingScheme scheme = kDefaultScheme) const;
231  double MeanX(rb::AveragingScheme scheme = kDefaultScheme) const {return MeanXYZ(scheme).X();}
232  double MeanY(rb::AveragingScheme scheme = kDefaultScheme) const {return MeanXYZ(scheme).Y();}
233  double MeanZ(rb::AveragingScheme scheme = kDefaultScheme) const {return MeanXYZ(scheme).Z();}
234  double MeanTNS(rb::AveragingScheme scheme = kDefaultScheme) const;
235 
236  /// Gets the min/max/mean all at once, called by the functions above
237  void MinMaxMeanXYZ(TVector3& lo, TVector3& hi, TVector3& mean,
238  rb::AveragingScheme scheme = kDefaultScheme) const;
239 
240  //@}
241 
242  /// \name Extents
243  /// The size of a box containing all the hits
244  //@{
245  TVector3 ExtentXYZ() const {return MaxXYZ()-MinXYZ();}
246  double ExtentV(geo::View_t view) const {return MaxV(view)-MinV(view);}
247  double ExtentX() const {return ExtentXYZ().X();}
248  double ExtentY() const {return ExtentXYZ().Y();}
249  double ExtentZ() const {return ExtentXYZ().Z();}
250  unsigned int ExtentPlane(geo::View_t view = geo::kXorY) const {return MaxPlane(view)-MinPlane(view)+1;}
251  unsigned int ExtentCell(geo::View_t view) const;
252  double ExtentTNS() const {return MaxTNS()-MinTNS();}
253  //@}
254 
255  /// \name Connectedness measures
256  //@{
257  /// \brief Longest run of adjacent planes with hits
258  ///
259  /// \a view may be kXorY to require hits on every plane no matter the view
260  int MostContiguousPlanes(geo::View_t view) const;
261 
262  /// \brief Longest run of adjacent planes with no hits
263  ///
264  /// \a view may be kXorY to allow a hit in either view to end the run
265  int MostMissingPlanes(geo::View_t view) const;
266 
267  /// \brief Total number of missing planes in cluster.
268  ///
269  /// Put simply, number of planes spanned minus number of hit planes.
270  /// If track is 2D, only planes in that view are counted.
271  int NMissingPlanes(geo::View_t view) const;
272  //@}
273 
274 
275  /// \brief Put the cells in the cluster into a standard order.
276  ///
277  /// Helps ensure that two clusters composed of the same cell hits
278  /// present in exactly the same way.
279  void StandardSort();
280 
282  kInitializeTotalGeV = 0, ///< Default, this is the initial/only call
283  kResetTotalGeV = 1 ///< Altering the already-set value is not an error
284  };
285 
286  /// \brief Store the current result of \ref TotalGeV / \ref CalorimetricEnergy
287  ///
288  /// If you intend to lock in the current calibrations, call this method
289  /// immediately before storing the product in the event. The result can be
290  /// retrieved by passing \ref kUsePrecalcEnergy to \ref TotalGeV or \ref
291  /// CalorimetricEnergy.
292  ///
293  /// Multiple calls are considered to be an error. If you really do need to
294  /// update the energy, pass \ref kResetTotalGeV.
295  void SavePrecalcTotalGeV(ESaveGeVMode savemode);
296 
297  bool operator< (const Cluster& other) const;
298 
299  friend std::ostream& operator << (std::ostream& o, const Cluster& c);
300 
301  protected:
302  /// \name Internal helpers
303  //@{
304 
305 
306  /// Helper. Resizes weights vectors to match cell vectors
307  void EnsureWeightAlloc();
308 
309  private:
310  /// Optimized form of TotalGeV, only valid for actual Clusters
311  double TotalGeVFastClusterOnly() const;
312  //@}
313 #endif // __GCCXML__
314  protected:
315 
316  geo::View_t fView; ///< view this cluster is in
317  art::PtrVector<rb::CellHit> fXCell; ///< collection of x-view cells in cluster
318  art::PtrVector<rb::CellHit> fYCell; ///< collection of y-view cells in cluster
319  std::vector<double> fXWeights; ///< Weights, matching cell indexing
320  std::vector<double> fYWeights; ///< May be empty, means all weights are 1
321  int fID; ///< ID for cluster
322  bool fNoiseCluster; ///< flag for whether this is a noise cluster
323 
324  double fPrecalcTotalGeV; ///< -1 = uninitialized
325  };
326 }
327 
328 #endif // RBCLUSTER_H
329 /////////////////////////////////////////////////////////////////////////////////////////////
double MinV(geo::View_t view) const
Definition: Cluster.cxx:454
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
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
Default, this is the initial/only call.
Definition: Cluster.h:282
double fPrecalcTotalGeV
-1 = uninitialized
Definition: Cluster.h:324
const art::PtrVector< rb::CellHit > & XCells() const
Get all cells from the x-view.
Definition: Cluster.h:124
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
double MinZ() const
Definition: Cluster.h:206
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
bool Is2D() const
Definition: Cluster.h:96
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
double MaxV(geo::View_t view) const
Definition: Cluster.cxx:500
A collection of associated CellHits.
Definition: Cluster.h:47
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:233
double ExtentZ() const
Definition: Cluster.h:249
double TotalGeVFastClusterOnly() const
Optimized form of TotalGeV, only valid for actual Clusters.
Definition: Cluster.cxx:400
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
void SetWeight(unsigned int globalIdx, double weight)
Set weight of the cell at this index.
Definition: Cluster.cxx:327
double ExtentX() const
Definition: Cluster.h:247
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
double MaxX() const
Definition: Cluster.h:217
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
double MeanX(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:231
Use a value computed by a previous call to SavePrecalcTotalGeV.
Definition: Cluster.h:179
double MinTNS() const
Definition: Cluster.cxx:482
TVector3 ExtentXYZ() const
Definition: Cluster.h:245
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 MeanY(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:232
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
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 art::PtrVector< rb::CellHit > & YCells() const
Get all cells from the x-view.
Definition: Cluster.h:126
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
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
Var weights
void SetID(int id)
Definition: Cluster.h:74
double MinY() const
Definition: Cluster.h:205
Perform a "2 point" Hough transform on a collection of hits.
double MaxTNS() const
Definition: Cluster.cxx:528
size_type size() const
Definition: PtrVector.h:308
double ExtentY() const
Definition: Cluster.h:248
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
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
double MaxZ() const
Definition: Cluster.h:219
TVector3 MinXYZ() const
Definition: Cluster.cxx:446
Definition: event.h:1
double ExtentV(geo::View_t view) const
Definition: Cluster.h:246
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
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
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
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
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
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void EnsureWeightAlloc()
Helper. Resizes weights vectors to match cell vectors.
Definition: Cluster.cxx:712
virtual void Clear()
Forget about all owned cell hits.
Definition: Cluster.cxx:279
In addition to the intrinsic weights.
Definition: Cluster.h:30
double ExtentTNS() const
Definition: Cluster.h:252
double MinX() const
Definition: Cluster.h:204
const AveragingScheme kDefaultScheme
Definition: Cluster.h:32
Simple object representing a (plane, cell) pair.
std::vector< double > fYWeights
May be empty, means all weights are 1.
Definition: Cluster.h:320
EEnergyCalcScheme
Definition: Cluster.h:177
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
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
double MaxY() const
Definition: Cluster.h:218