FLSHit.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief FLSHit class
3 /// \author brebel@fnal.gov
4 /// \date
5 ////////////////////////////////////////////////////////////////////////
6 #include "Simulation/FLSHit.h"
7 #include <cmath>
8 
9 namespace sim{
10 
11  //---------------------------------------------------------------
13  fId (0)
14  ,fPlaneId(0)
15  ,fCellId (0)
16  ,fPDG (0)
17  ,fTrackId(0)
18  { }
19  //__________________________________________________________
21 
22  //__________________________________________________________
23  bool FLSHit::operator==(const FLSHit& rhs) const
24  {
25  if (fId == rhs.fId && GetEntryT() == rhs.GetEntryT()) return true;
26  return false;
27  }
28 
29  //__________________________________________________________
30  bool FLSHit::operator<(const FLSHit& rhs) const
31  {
32  if (fId<rhs.fId) return true;
33  if (fId == rhs.fId)
34  if (GetEntryT() < rhs.GetEntryT()) return true;
35  return false;
36  }
37 
38  //__________________________________________________________
39  void FLSHit::Clear() {
40 
41  // Clear almost everything. We don't want to clear fId, though
42 
43  fPlaneId = 0;
44  fCellId = 0;
45  fTrackId = 0;
46  fPDG = 0;
47 
48  fX.clear();
49  fY.clear();
50  fZ.clear();
51  fT.clear();
52  fEnergy.clear();
53  fVectorEdep.clear();
54  fVectorEdepBirks.clear();
55  fVectorNCerenkov.clear();
56  }
57 
58 //----------------------------------------------------------------------------
59 // Get path length of a particular step
60 float FLSHit::GetPathLength(const int step) const{
61 
62  const double x = fX[step+1] - fX[step];
63  const double y = fY[step+1] - fY[step];
64  const double z = fZ[step+1] - fZ[step];
65 
66  return sqrt(x*x + y*y + z*z);
67 }
68 //----------------------------------------------------------------------------
69 // Get path length of all steps in FLSHit
71 
72  const int nsteps = GetNSteps();
73 
74  double path_length = 0.0;
75  for(int i=0; i<nsteps; ++i){
76  path_length += GetPathLength(i);
77  }
78 
79  return path_length;
80 }
81 
82 //----------------------------------------------------------------------------
83 // Get energy loss of a particular step
84 float FLSHit:: GetEnergyLoss(const int step) const{
85  return fEnergy[step] - fEnergy[step+1];
86 }
87 
88 //----------------------------------------------------------------------------
89 // Get total energy loss in the entire FLSHit
91  return GetEntryEnergy() - GetExitEnergy();
92 }
93 
94 //----------------------------------------------------------------------------
95 float FLSHit::GetXYZAverage(const std::vector<float>& xyz, const int step) const{
96  return (xyz[step] + xyz[step+1]) * 0.5;
97 }
98 
99 //----------------------------------------------------------------------------
100 float FLSHit::GetXYZAverage(const std::vector<float>& xyz) const{
101  const int nsteps = GetNSteps();
102 
103  double xyz_av = 0.0;
104  double weight = 0.0;
105 
106  for(int i=0; i<nsteps; ++i){
107  const double path_length = GetPathLength(i);
108 
109  weight += path_length;
110  xyz_av += path_length * GetXYZAverage(xyz, i);
111  }
112 
113  return xyz_av / weight;
114 }
115 
116  //----------------------------------------------------------------------------
117  // Performs various checks on the size of the data and its content
118  // Has an option to print everything out
119  bool FLSHit::IsFLSHitReasonable(const bool print_errors
120  ,std::ostream& os) const
121  {
122  bool is_everything_ok = true;
123 
124  /// All sizes should be at least 2 (Entry+Exit)
125  if(fX.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fX is empty\n";}}
126  if(fY.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fY is empty\n";}}
127  if(fZ.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fZ is empty\n";}}
128  if(fT.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fT is empty\n";}}
129  if(fEnergy.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fEnergy is empty\n";}}
130  if(fVectorEdep.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fVectorEdep is empty\n";}}
131  if(fVectorEdepBirks.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fVectorEdepBirks is empty\n";}}
132 
133  /// Compare the sizes. Should be identical
134  if(fX.size() != fY.size() ||
135  fX.size() != fZ.size() ||
136  fX.size() != fT.size() ||
137  fX.size() != fEnergy.size()
138  ){
139  is_everything_ok = false;
140  if(print_errors){
141  os<<"Vectors fX, fY, fZ, fT, fEnergy should have the same size, but they have sizes";
142  os<<" "<<fX.size();
143  os<<" "<<fY.size();
144  os<<" "<<fZ.size();
145  os<<" "<<fT.size();
146  os<<" "<<fEnergy.size();
147  os<<" respectively.\n";
148  }// end of printing
149  }// end of checking
150 
151  /// Compare the sizes. Should be identical
152  if(fVectorEdep.size() != fVectorEdepBirks.size()){
153  is_everything_ok = false;
154  if(print_errors){
155  os<<"Vectors fVectorEdep and fVectorEdepBirks, but they have sizes";
156  os<<" "<<fVectorEdep.size();
157  os<<" "<<fVectorEdepBirks.size();
158  os<<" respectively.\n";
159  }// end of printing
160  }// end of checking
161 
162  /// Compare the sizes. fX should have one more entry than fVectorEdep
163  if(fX.size() != fVectorEdep.size() + 1){
164  is_everything_ok = false;
165  if(print_errors){
166  os<<"Vector fX should have one more entry than fVectorEdep, but they have sizes";
167  os<<" "<<fX.size();
168  os<<" "<<fVectorEdep.size();
169  os<<" respectively.\n";
170  }// end of printing
171  }// end of checking
172 
173  /// Time should increase monotonically
175  is_everything_ok = false;
176  if(print_errors){
177  os<<"Vector fT did not increase monotonically\n";
178  }// end of printing
179  }// end of checking
180 
181  /// Energy should decrease monotonically
183  is_everything_ok = false;
184  if(print_errors){
185  os<<"Vector fEnergy did not decrease monotonically\n";
186  }// end of printing
187  }// end of checking
188 
189  /// ID of the FLSHit should be non-zero
190  if(fId == 0){
191  is_everything_ok = false;
192  if(print_errors){
193  os<<"fID is zero\n";
194  }// end of printing
195  }// end of checking fID
196 
197  return is_everything_ok;
198  }// end of IsFLSHitReasonable
199 
200  //----------------------------------------------------------------------
201  // ostream operator.
202  //
203  std::ostream& operator<< (std::ostream & o, const FLSHit & a)
204  {
205  o << std::setiosflags(std::ios::fixed) << std::setprecision(2);
206  o << " CellUniqueId " << std::setw(10) << std::right << a.GetCellUniqueId()
207  << " Entry (x,y,z): (" << a.GetEntryX() << "," << a.GetEntryY() << "," << a.GetEntryZ() << ")"
208  << " Entry E: " << std::setw(5) << std::right << a.GetEntryEnergy()
209  << " Exit (x,y,z): (" << a.GetExitX() << "," << a.GetExitY() << "," << a.GetExitZ() << ")"
210  << " Exit E: " << std::setw(5) << std::right << a.GetExitEnergy()
211  << " Entry Time: " << std::setw(5) << std::right << a.GetEntryT()
212  << " Exit Time: " << std::setw(5) << std::right << a.GetExitT()
213  << " E deposited: " << std::setw(5) << std::right << a.GetEdep()
214  << " E Birks Weight: " << std::setw(5) << std::right << a.GetEdepBirks()
215  << " N Cerenkov Photons " << std::setw(5) << std::right << a.GetNCerenkov()
216  << " Plane,Cell: " << a.GetPlaneID() << "," << a.GetCellID()
217  << " Track ID: " << std::setw(5) << std::right << a.GetTrackID();
218 
219  return o;
220  }
221 
222 }// end namespace sim
223 
224 ////////////////////////////////////////////////////////////////////////
geo::CellUniqueId GetCellUniqueId() const
Unique Cell ID.
Definition: FLSHit.h:41
float GetEntryX() const
Entry point of the particle (position, time and energy)
Definition: FLSHit.h:48
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
std::vector< float > fY
Start or end of the step inside scintillator. Y position, cm.
Definition: FLSHit.h:192
int GetPlaneID() const
Plane ID.
Definition: FLSHit.h:37
float GetEnergyLoss(const int step) const
Get energy loss of a particular step.
Definition: FLSHit.cxx:84
const Var weight
float GetEntryEnergy() const
Definition: FLSHit.h:52
float GetPathLength(const int step) const
Get path length of a particular step.
Definition: FLSHit.cxx:60
int GetCellID() const
Cell ID.
Definition: FLSHit.h:39
float GetExitT() const
Definition: FLSHit.h:57
const int nsteps
std::vector< float > fZ
Start or end of the step inside scintillator. Z position, cm.
Definition: FLSHit.h:193
float GetTotalPathLength() const
Get path length of all steps in FLSHit.
Definition: FLSHit.cxx:70
int GetNSteps() const
Number of steps in FLSHit. Should be equal to GetNPoints() - 1.
Definition: FLSHit.h:63
T sqrt(T number)
Definition: d0nt_math.hpp:156
float GetTotalEnergyLoss() const
Get total energy loss in the entire FLSHit.
Definition: FLSHit.cxx:90
std::vector< float > fVectorEdep
Amount of energy deposited during a Geant4 step (GeV)
Definition: FLSHit.h:197
unsigned short int fCellId
Cell number.
Definition: FLSHit.h:187
float GetNCerenkov() const
Get total N Cerenkov photons.
Definition: FLSHit.h:35
std::vector< double > fT
Start or end of the step inside scintillator. Time, ns.
Definition: FLSHit.h:195
std::vector< float > fVectorNCerenkov
Amount of Cerenkov light emitteded in scintillator.
Definition: FLSHit.h:199
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
float GetExitX() const
Exit point of the particle (position, time and energy)
Definition: FLSHit.h:54
bool IsFLSHitReasonable(const bool print_errors=false, std::ostream &error_printint_stream=std::cout) const
Definition: FLSHit.cxx:119
bool operator<(const FLSHit &rhs) const
Definition: FLSHit.cxx:30
unsigned short int fPlaneId
Plane number.
Definition: FLSHit.h:186
geo::CellUniqueId fId
Unique cell ID.
Definition: FLSHit.h:184
std::ostream & operator<<(std::ostream &o, const FLSHit &a)
Definition: FLSHit.cxx:203
std::vector< float > fVectorEdepBirks
Amount of energy deposited during a Geant4 step ; Birks weighted (GeV)
Definition: FLSHit.h:198
const double a
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:15
z
Definition: test.py:28
float GetEdepBirks() const
Get total Energy with Birks suppression deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:33
int fPDG
PDG code of parent particle.
Definition: FLSHit.h:188
float GetEntryZ() const
Definition: FLSHit.h:50
float GetExitEnergy() const
Definition: FLSHit.h:58
float GetEntryT() const
Definition: FLSHit.h:51
bool operator==(const FLSHit &rhs) const
Definition: FLSHit.cxx:23
int fTrackId
Index number of the Global track producing hit.
Definition: FLSHit.h:189
int GetTrackID() const
Track ID.
Definition: FLSHit.h:45
std::vector< float > fX
Start or end of the step inside scintillator. X position, cm.
Definition: FLSHit.h:191
void Clear()
Clear the FLS hit.
Definition: FLSHit.cxx:39
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:31
float GetExitY() const
Definition: FLSHit.h:55
float GetExitZ() const
Definition: FLSHit.h:56
float GetXYZAverage(const std::vector< float > &xyz, const int step) const
Definition: FLSHit.cxx:95
bool DoesDecreaseMonotonically(const std::vector< T > &input_vector) const
Checks whether the input vector decreases monotonically.
Definition: FLSHit.h:166
float GetEntryY() const
Definition: FLSHit.h:49
std::vector< float > fEnergy
Start or end of the step inside scintillator. Energy, GeV.
Definition: FLSHit.h:194
bool DoesIncreaseMonotonically(const std::vector< T > &input_vector) const
Checks whether the input vector increases monotonically.
Definition: FLSHit.h:149