ParticleIDAlg.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file ParticleIDAlg.h
3 /// \brief Calculate dE/dx and log likelihoods for different particle hypotheses.
4 /// This information will be of use for identification of any particle that
5 /// uses a Shower as a starting point.
6 ///
7 
8 ///
9 /// Upon instantiation the histograms are loaded. Ideally one would reuse an instance of this
10 /// class to save reloading of the histograms.
11 ///
12 /// \author smith@physics.umn.edu
13 ////////////////////////////////////////////////////////////////////////
14 
15 
16 #ifndef ParticleIDAlg_h
17 #define ParticleIDAlg_h
18 #include <vector>
19 #include "fhiclcpp/ParameterSet.h"
23 #include "Geometry/Geometry.h"
24 #include "ShowerLID/NuEEnergyAlg.h"
27 #include "Geometry/Geometry.h"
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/WeightedHit.h"
31 #include "RecoBase/Cluster.h"
32 #include "RecoBase/Track.h"
33 #include "RecoBase/Shower.h"
34 #include "Utilities/func/MathUtil.h"
35 #include "TGeoManager.h"
36 #include "TGeoNode.h"
37 
38 
39 
40 namespace rb {
41  class Cluster;
42  class Track;
43  class CellHit;
44  struct WeightedHit;
45 }
46 
47 
48 namespace slid {
49 
50  class ParticleIDAlg {
51  public:
52  explicit ParticleIDAlg(fhicl::ParameterSet const& pset,
54  NuEEnergyAlg* NuEEnergyAlgIn,
55  bool loadLibs = true);
56 
57  ~ParticleIDAlg();
58 
59  /// \brief Reinitialize parameters read from FCL file
60  void reconfigure(const fhicl::ParameterSet & pset,
61  novadaq::cnv::DetId detId);
62 
63  /// \brief Calculate transverse dE/dx log-likelihood for specific particle hypothesis
64  double DedxTransLL(const slid::DedxParticleType partHyp,
65  const rb::Shower* vShower,
66  const std::vector< const rb::Shower* > showercol,
67  const art::Event& evt);
68 
69 
70  /// \brief Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis
71  double DedxLongLL(const slid::DedxParticleType partHyp,
72  const rb::Shower* vShower,
73  const std::vector< const rb::Shower* > showercol,
74  const art::Event& evt);
75 
76  /// \brief Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis
77  /// assuming that the particle direction if opposite to the reconstructed direction
78  double DedxInverseLongLL(const slid::DedxParticleType partHyp,
79  const rb::Shower* vShower,
80  const std::vector< const rb::Shower* > showercol,
81  const art::Event& evt);
82 
83  /// \brief Set rb::Prong to be analyzed. This must be set before any calculations can be done.
84  void SetShower(const rb::Shower* vShower,
85  const std::vector< const rb::Shower* > showercol,
86  const art::Event& evt);
87 
88 
89 
90  ///\brief Calculate the gap between the vertex and the start of the shower
91  double GetGapVertexToShowerStart(const rb::Shower* vShower, TVector3 evtVtx,
92  const art::Event& evt);
93 
94  ///\brief Calculate the mass of candidate and any other shower in the slice that is closest to the pi0 mass.
95  double Pi0Mass(const std::vector< const rb::Shower* > showerCollection,
96  unsigned int iShower, const art::Event& evt);
97 
98  ///\brief Index in showerCollection of the shower for which the reconstructed
99  /// invariant mass was closest to pi0 mass for the shower at index iShower.
100  int Pi0SecondPhoton(const std::vector< const rb::Shower* > showerCollection,
101  unsigned int iShower, const art::Event& evt);
102 
103  ///\brief Calculate the shower radius
104  double Radius(const std::vector<const rb::Shower*> showerCollection,
105  unsigned int iShower, const art::Event& evt);
106 
107  ///\brief return longitudinal dedx for a specified plane index
108  double PlaneLongDedx(unsigned int pIdx);
109 
110  ///\brief return transverse dedx for specified transverse plane
111  double PlaneTransDedx(unsigned int tpIdx);
112 
113  ///\brief return transverse dedx for specified transverse plane
114  double CellPlaneDedx(int tpIdx, unsigned int pIdx);
115 
116  ///\brief return shower radius for a plane index
117  double PlaneRadius(unsigned int pIdx);
118 
119  ///\brief Return point of intersection between shower core and plane
120  TVector3 PlaneHitXYZ(unsigned int pIdx);
121 
122  ///\brief Return point of intersection cell between shower core and plane
123  double PlaneHitCell(unsigned int pIdx);
124  ///\brief Return the gloabel plane index for an index in the shower coordinate
125  unsigned int PlaneToGlobalPlane(unsigned int pIdx);
126 
127  ///\brief Return the energy weighted mean shower positon for the plane in cm
128  double PlaneCentroid(unsigned int pIdx);
129 
130  ///\brief Return the asymmetry of the shower
131  double ShowerAsym() const {return fAsym;}
132 
133  ///\brief Return the inertia calculation
134  TVector3 InertiaXYZ() const {return fIxyz;}
135 
136  /// \brief Get distance of closest approach between shower and vertex
137  double PointDoca(TVector3 vtx, TVector3 start, TVector3 stop);
138 
139  /// \brief Number of MIP planes within the first 20 planes of the start of the shower
140  int GetMIPPlanes(const rb::Shower* shower,
141  const std::vector< const rb::Shower* > showercol,
142  const art::Event& evt);
143 
144  /// Map of the longitudinal ll by paricle type
145  std::map<int, float> fPartLongLL;
146 
147  /// Map of the transverse ll by particle type
148  std::map<int, float> fPartTransLL;
149 
150  private:
151 
152  /// \brief Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector.
153  double CalcPlaneLongitudinalDedx(unsigned int plane);
154 
155  /// \brief Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index.
156  double CalcPlaneTransverseDedx(unsigned int iTransversePlane);
157 
158 
159  /// \brief Calculate the transverse dE/dx in a given tranverse cell and plane. Input is the transverse slice index, including plus and minus to the shower core, and plane number.
160  double CalcCellPlaneTransverseDedx(int iTransversePlane, unsigned int iPlane);
161 
162  /// \brief Return a map off all hits in a given longitudinal plane indexed by the longitudinalplane
163  /// for shower already read into object by SetShower
164  std::map< int, std::map< int, int > > CalcPlaneHits(const rb::Shower* vShower );
165 
166  /// \brief Return a map of mean cell posisiton for each plane
167  /// for shower already read into object by SetShower
168  std::map< int, double > CalcPlaneCentroid( const rb::Shower* vShower);
169 
170  /// \brief Sum the energy in a longitudinal plane
171  double CalcEnergyInLongitudinalPlane(int iPlane);
172 
173  /// \brief Determine detector XY region
174  void CalcDetectorXYRegionIndex();
175 
176  /// \brief Get longitudinal dE/dx probability for a given longitudinal plane, energy bin, and particle hypothesis
177  double CalcPlaneLongDedxProb(DedxParticleType particleType, int igev, unsigned int iPlane, unsigned int plane);
178 
179  /// \brief Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis
180  double CalcPlaneTransDedxProb(DedxParticleType particleType, int igev, unsigned int iTransversePlane);
181 
182  /// \brief From Jianming's code, still working out changes.
183  double CalcInterPlaneDedxProb(DedxParticleType type, unsigned int iLongPlane, unsigned int plane);
184 
185  /// \brief From Jianming's code, still working out changes.
186  double CalcInterCellTransDedxProb(const DedxParticleType type, unsigned int iTransPlane);
187 
188  /// \brief Load all of the histograms from files into arrays of TH1D* objects. The details of handling the
189  /// histograms are taken care of within object of the helper class slid::DedxDistribution.
190  ///
191  bool LoadDedxHistogramFiles(const fhicl::ParameterSet & pset, novadaq::cnv::DetId geom);
192 
193  /// \brief Calculate path length of shower through this plane.
194  double CalcTrkHitPath(unsigned int plane);
195 
196  /// \brief Calculate energy bin indices and energy values from the histograms
197  /// corresponding energy of shower to be used to interpolate between bins in energy
198  void CalculateEnergyBinIndex();
199 
200  /// \brief Calculate the asymmetry and inertia of the shower
201  void CalcAsymIneria();
202 
203  /// \brief Helper function for interia and asymmetry calculation
204  TVector3 GetCellNodePos(unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2);
205 
206  /// data/MC flag
207  //bool isData;
208 
209  /// event ID
211 
212  /// Closest mass to pi0 mass formed by this shower and any other shower in the slice.
214 
215  /// dE/dx distributions (histogram files and data) for each detector type and particle category
216  std::vector<slid::DedxDistribution> fDedxDistributions;
217 
218  /// rb::Shower object for which particle ID is to be calculated. This is cached.
220  std::vector< const rb::Shower* > fShowerCol;
221 
222  /// Shower energy
224 
225  /// geometry
227 
228  /// Nue energy algorithms
230 
231 
232  /// Map of longitudinal dedx by plane. The first element is
233  /// plane number and the second is the dedx for that plane.
234  /// This is used for lazy caching of the dedx plane values
235  /// in order to save recalculating for the same shower.
236  std::map< int, double> fDedxLongByPlane;
237 
238  /// Map of transverse dedx by plane. The first element is
239  /// plane number and the second is the dedx for that plane.
240  /// This is used for lazy caching of the dedx plane values
241  /// in order to save recalculating for the same shower.
242  std::map< int, double> fDedxTransByPlane;
243 
244 
245  /// Map where the first element is
246  /// plane number and the second is a map with the key being the
247  /// cell number in the plane and the second value is the global
248  /// cell index in a shower
249  std::map< int, std::map< int, int > > fPlaneHits;
250 
251  /// Map where the key is the index of the plane in the shower and the second is the global plane
252  /// This map saves looping in some of the utility functions
253  std::map<unsigned int, unsigned int> fMapIdxToGlobalPlane;
254 
255  /// Map where the first index is
256  /// plane number and the second is the energy weighted center position for that plane
257  std::map< int, double > fPlaneCentroid;
258 
259  /// Map of transverse energy for each transverse plane
260  std::map<int, double > fTransEnergy;
261 
262  /// Map of path length for each transverse plane
263  std::map<int, double > fTransPathLength;
264 
265  /// Index of the other photon shower used for pi0 mass calculation
266  /// for a given shower
268 
269  /// Detector sub-region bin. Calculated once and cached for a given shower/event combo.
270  int iReg;
271 
272  /// Shower start position
273  TVector3 pos;
274 
275  /// Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.
276  int igev0;
277  int igev1;
278 
279  /// Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event combo.
280  double e0;
281  double e1;
282 
283  /// shower moment of ineria
284  TVector3 fIxyz;
285  /// shower asymmetry
286  double fAsym;
287 
288  /// Upper bound of MIP Range
290 
291  /// Lower bound of MIP Range
292  double fMipRangeLow;
293 
294  /// Detector ID
296 
297  };
298 }
299 #endif // ParticleIDAlg_h
300 
301 ////////////////////////////////////////////////////////////////////////
302 
303 
304 
novadaq::cnv::DetId fDetId
Detector ID.
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
std::map< int, double > fTransEnergy
Map of transverse energy for each transverse plane.
double fClosestPi0Mass
Closest mass to pi0 mass formed by this shower and any other shower in the slice. ...
TVector3 InertiaXYZ() const
Return the inertia calculation.
This is a helper class for ParticleIDAlg that provides a tidy structure in which to hold the dE/dx hi...
int iReg
Detector sub-region bin. Calculated once and cached for a given shower/event combo.
double e0
Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event ...
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
std::vector< slid::DedxDistribution > fDedxDistributions
dE/dx distributions (histogram files and data) for each detector type and particle category ...
double fShowerEnergyGev
Shower energy.
double fMipRangeLow
Lower bound of MIP Range.
double Radius(int A, double Ro=constants::kNucRo)
DedxParticleType
An enum used to give allowed particle types a visible name in the code. Note that for electron...
std::vector< const rb::Shower * > fShowerCol
double ShowerAsym() const
Return the asymmetry of the shower.
int evt
TVector3 pos
Shower start position.
art::ServiceHandle< geo::Geometry > fGeom
geometry
std::map< int, float > fPartTransLL
Map of the transverse ll by particle type.
Perform a "2 point" Hough transform on a collection of hits.
double fAsym
shower asymmetry
art::EventID eventID
data/MC flag
std::map< int, double > fTransPathLength
Map of path length for each transverse plane.
std::map< int, double > fDedxTransByPlane
TVector3 fIxyz
shower moment of ineria
Calculates deposited and corrected energy of the electron shower and of electron flavoured neutrino...
A rb::Prong with a length.
Definition: Shower.h:18
void geom(int which=0)
Definition: geom.C:163
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
std::map< int, double > fPlaneCentroid
std::map< int, std::map< int, int > > fPlaneHits
std::map< int, float > fPartLongLL
Map of the longitudinal ll by paricle type.
int igev0
Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.
Encapsulate the geometry of one entire detector (near, far, ndos)
std::map< int, double > fDedxLongByPlane
double fMipRangeHigh
Upper bound of MIP Range.