All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CalorimetryAlg.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file CalorimetryAlg.cxx
3 //
4 // \brief Functions to calculate dE/dx. Based on code in Calorimetry.cxx
5 //
6 // andrzej.szelc@yale.edu
7 //
8 ////////////////////////////////////////////////////////////////////////
9 
10 // LArSoft includes
11 #include "larreco/Calorimetry/CalorimetryAlg.h"
12 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
13 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
14 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
15 #include "lardataalg/DetectorInfo/DetectorProperties.h"
16 #include "lardataobj/RecoBase/Hit.h"
17 #include "larevt/CalibrationDBI/Interface/ElectronLifetimeProvider.h"
18 #include "larevt/CalibrationDBI/Interface/ElectronLifetimeService.h"
19 
20 namespace calo {
21 
22  //--------------------------------------------------------------------
23  CalorimetryAlg::CalorimetryAlg(const Config& config)
24  : fCalAmpConstants{config.CalAmpConstants()}
25  , fCalAreaConstants{config.CalAreaConstants()}
26  , fUseModBox{config.CaloUseModBox()}
27  , fLifeTimeForm{config.CaloLifeTimeForm()}
28  , fDoLifeTimeCorrection{config.CaloDoLifeTimeCorrection()}
29  {
30  if (fLifeTimeForm != 0 and fLifeTimeForm != 1) {
31  throw cet::exception("CalorimetryAlg")
32  << "Unknow CaloLifeTimeForm " << fLifeTimeForm << '\n'
33  << "Must select either '0' for exponential or '1' for exponential + "
34  "constant.\n";
35  }
36  }
37 
38  //------------------------------------------------------------------------------------//
39  // Functions to calculate the dEdX based on the AMPLITUDE of the pulse
40  // ----------------------------------------------------------------------------------//
41  double
42  CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
43  detinfo::DetectorPropertiesData const& det_prop,
44  recob::Hit const& hit,
45  double const pitch,
46  double const T0) const
47  {
48  return dEdx_AMP(
49  clock_data, det_prop, hit.PeakAmplitude() / pitch, hit.PeakTime(), hit.WireID().Plane, T0, det_prop.Efield());
50  }
51 
52  ///\todo The plane argument should really be for a view instead
53  // ----------------------------------------------------------------------------------//
54  double
55  CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
56  detinfo::DetectorPropertiesData const& det_prop,
57  double const dQ,
58  double const time,
59  double const pitch,
60  unsigned int const plane,
61  double const T0) const
62  {
63  double const dQdx = dQ / pitch; // in ADC/cm
64  return dEdx_AMP(clock_data, det_prop, dQdx, time, plane, T0, det_prop.Efield());
65  }
66 
67  // ----------------------------------------------------------------------------------//
68  double
69  CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
70  detinfo::DetectorPropertiesData const& det_prop,
71  double const dQdx,
72  double const time,
73  unsigned int const plane,
74  double const T0) const
75  {
76  double const fADCtoEl = fCalAmpConstants[plane];
77  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
78  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, det_prop.Efield());
79  }
80 
81  //------------------------------------------------------------------------------------//
82  // Functions to calculate the dEdX based on the AMPLITUDE of the pulse with EField
83  // ----------------------------------------------------------------------------------//
84  double
85  CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
86  detinfo::DetectorPropertiesData const& det_prop,
87  recob::Hit const& hit,
88  double const pitch,
89  double const T0,
90  double const EField) const
91  {
92  return dEdx_AMP(
93  clock_data, det_prop, hit.PeakAmplitude() / pitch, hit.PeakTime(), hit.WireID().Plane, T0, EField);
94  }
95 
96  ///\todo The plane argument should really be for a view instead
97  // ----------------------------------------------------------------------------------//
98  double
99  CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
100  detinfo::DetectorPropertiesData const& det_prop,
101  double const dQ,
102  double const time,
103  double const pitch,
104  unsigned int const plane,
105  double const T0,
106  double const EField) const
107  {
108  double const dQdx = dQ / pitch; // in ADC/cm
109  return dEdx_AMP(clock_data, det_prop, dQdx, time, plane, T0, EField);
110  }
111 
112  // ----------------------------------------------------------------------------------//
113  double
114  CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
115  detinfo::DetectorPropertiesData const& det_prop,
116  double const dQdx,
117  double const time,
118  unsigned int const plane,
119  double const T0,
120  double const EField) const
121  {
122  double const fADCtoEl = fCalAmpConstants[plane];
123  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
124  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, EField);
125  }
126 
127  //------------------------------------------------------------------------------------//
128  // Functions to calculate the dEdX based on the AREA of the pulse
129  // ----------------------------------------------------------------------------------//
130  double
131  CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
132  detinfo::DetectorPropertiesData const& det_prop,
133  recob::Hit const& hit,
134  double const pitch,
135  double const T0) const
136  {
137  return dEdx_AREA(
138  clock_data, det_prop, hit.Integral() / pitch, hit.PeakTime(), hit.WireID().Plane, T0, det_prop.Efield());
139  }
140 
141  // ----------------------------------------------------------------------------------//
142  double
143  CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
144  detinfo::DetectorPropertiesData const& det_prop,
145  double const dQ,
146  double const time,
147  double const pitch,
148  unsigned int const plane,
149  double const T0) const
150  {
151  double const dQdx = dQ / pitch; // in ADC/cm
152  return dEdx_AREA(clock_data, det_prop, dQdx, time, plane, T0, det_prop.Efield());
153  }
154 
155  // ----------------------------------------------------------------------------------//
156  double
157  CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
158  detinfo::DetectorPropertiesData const& det_prop,
159  double const dQdx,
160  double const time,
161  unsigned int const plane,
162  double const T0) const
163  {
164  double const fADCtoEl = fCalAreaConstants[plane];
165  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
166  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, det_prop.Efield());
167  }
168 
169  //------------------------------------------------------------------------------------//
170  // Functions to calculate the dEdX based on the AREA of the pulse with EField
171  // ----------------------------------------------------------------------------------//
172  double
173  CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
174  detinfo::DetectorPropertiesData const& det_prop,
175  recob::Hit const& hit,
176  double const pitch,
177  double const T0,
178  double const EField) const
179  {
180  return dEdx_AREA(
181  clock_data, det_prop, hit.Integral() / pitch, hit.PeakTime(), hit.WireID().Plane, T0, EField);
182  }
183 
184  // ----------------------------------------------------------------------------------//
185  double
186  CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
187  detinfo::DetectorPropertiesData const& det_prop,
188  double const dQ,
189  double const time,
190  double const pitch,
191  unsigned int const plane,
192  double const T0,
193  double const EField) const
194  {
195  double const dQdx = dQ / pitch; // in ADC/cm
196  return dEdx_AREA(clock_data, det_prop, dQdx, time, plane, T0, EField);
197  }
198 
199  // ----------------------------------------------------------------------------------//
200  double
201  CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
202  detinfo::DetectorPropertiesData const& det_prop,
203  double const dQdx,
204  double const time,
205  unsigned int const plane,
206  double const T0,
207  double const EField) const
208  {
209  double const fADCtoEl = fCalAreaConstants[plane];
210  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
211  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, EField);
212  }
213 
214  // Apply Lifetime and recombination correction.
215  double
216  CalorimetryAlg::dEdx_from_dQdx_e(detinfo::DetectorClocksData const& clock_data,
217  detinfo::DetectorPropertiesData const& det_prop,
218  double dQdx_e,
219  double const time,
220  double const T0) const
221  {
222  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, det_prop.Efield());
223  }
224  double
225  CalorimetryAlg::dEdx_from_dQdx_e(detinfo::DetectorClocksData const& clock_data,
226  detinfo::DetectorPropertiesData const& det_prop,
227  double dQdx_e,
228  double const time,
229  double const T0,
230  double const EField) const
231  {
232  if (fDoLifeTimeCorrection) {
233  dQdx_e *= LifetimeCorrection(clock_data, det_prop, time, T0); // (dQdx_e in e/cm)
234  }
235 
236  if (fUseModBox) { return det_prop.ModBoxCorrection(dQdx_e, EField); }
237 
238  return det_prop.BirksCorrection(dQdx_e, EField);
239  }
240 
241  //------------------------------------------------------------------------------------//
242  // for the time being copying from Calorimetry.cxx - should be decided where
243  // to keep it.
244  // ----------------------------------------------------------------------------------//
245  double
246  calo::CalorimetryAlg::LifetimeCorrection(detinfo::DetectorClocksData const& clock_data,
247  detinfo::DetectorPropertiesData const& det_prop,
248  double const time,
249  double const T0) const
250  {
251  float const t = time - trigger_offset(clock_data);
252  double const timetick = sampling_rate(clock_data) * 1.e-3; // time sample in microsec
253  double const adjusted_time = t * timetick - T0 * 1e-3; // (in microsec)
254 
255  assert(fLifeTimeForm < 2);
256  if (fLifeTimeForm == 0) {
257  // Exponential form
258  double const tau = det_prop.ElectronLifetime();
259  return exp(adjusted_time / tau);
260  }
261 
262  // Exponential+constant form
263  auto const& elifetime_provider =
264  art::ServiceHandle<lariov::ElectronLifetimeService const>()->GetProvider();
265  return elifetime_provider.Lifetime(adjusted_time);
266  }
267 
268 } // namespace
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
Hit.h
Declaration of signal hit object.
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
hit
process_name hit
Definition: cheaterreco.fcl:51
detinfo::DetectorPropertiesData::ModBoxCorrection
double ModBoxCorrection(double dQdX) const
Definition: DetectorPropertiesData.cc:40
recob::Hit::PeakAmplitude
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
detinfo::DetectorPropertiesData::Efield
double Efield(unsigned int planegap=0) const
kV/cm
Definition: DetectorPropertiesData.cc:17
calo
process_name can override from command line with o or output calo
Definition: pid.fcl:40
and
return match has_match and(match.match_pdg==11 or match.match_pdg==-11)
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
geo_types.h
Definition of data types for geometry description.
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
calo::CalorimetryAlg::CalorimetryAlg
CalorimetryAlg(const fhicl::ParameterSet &pset)
Definition: CalorimetryAlg.h:57
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
detinfo::DetectorPropertiesData::BirksCorrection
double BirksCorrection(double dQdX) const
dQ/dX in electrons/cm, returns dE/dX in MeV/cm.
Definition: DetectorPropertiesData.cc:29
e
do i e
Definition: updateVersionList.sh:68
detinfo::trigger_offset
int trigger_offset(DetectorClocksData const &data)
Definition: DetectorClocksData.h:632
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
calo::CalorimetryAlg::LifetimeCorrection
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Definition: CalorimetryAlg.cxx:246
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
T0
process_name opdaq physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator T0
Definition: gen_protons.fcl:45
Generated by   doxygen 1.8.5