LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CalorimetryAlg.cxx
Go to the documentation of this file.
1 // \file CalorimetryAlg.cxx
3 //
4 // \brief Functions to calculate dE/dx. Based on code in Calorimetry.cxx
5 //
6 // andrzej.szelc@yale.edu
7 //
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  , fModBoxA{config.ModBoxA()}
30  , fModBoxBF("ModBoxB", config.ModBoxBTF1().c_str(), 0, 90)
31  , fBirksA{config.BirksA()}
32  , fBirksKF("BirksK", config.BirksKTF1().c_str(), 0, 90)
33  {
34  if (fLifeTimeForm != 0 and fLifeTimeForm != 1) {
35  throw cet::exception("CalorimetryAlg")
36  << "Unknow CaloLifeTimeForm " << fLifeTimeForm << '\n'
37  << "Must select either '0' for exponential or '1' for exponential + "
38  "constant.\n";
39  }
40 
41  // Set the parameters for the TF1s
42  std::vector<double> modboxb_param;
43  if (!config.ModBoxBParam(modboxb_param)) modboxb_param = {util::kModBoxB};
44 
45  for (unsigned i = 0; i < modboxb_param.size(); i++) {
46  fModBoxBF.SetParameter(i, modboxb_param[i]);
47  }
48 
49  std::vector<double> birksk_param;
50  if (!config.BirksKParam(birksk_param)) birksk_param = {util::kRecombk};
51 
52  for (unsigned i = 0; i < birksk_param.size(); i++) {
53  fBirksKF.SetParameter(i, birksk_param[i]);
54  }
55  }
56 
57  //------------------------------------------------------------------------------------//
58  // Functions to calculate the dEdX based on the AMPLITUDE of the pulse
59  // ----------------------------------------------------------------------------------//
60  double CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
61  detinfo::DetectorPropertiesData const& det_prop,
62  recob::Hit const& hit,
63  double const pitch,
64  double const T0) const
65  {
66  return dEdx_AMP(clock_data,
67  det_prop,
68  hit.PeakAmplitude() / pitch,
69  hit.PeakTime(),
70  hit.WireID().Plane,
71  T0,
72  det_prop.Efield());
73  }
74 
75  // ----------------------------------------------------------------------------------//
76  double CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
77  detinfo::DetectorPropertiesData const& det_prop,
78  double const dQdx,
79  double const time,
80  unsigned int const plane,
81  double const T0) const
82  {
83  double const fADCtoEl = fCalAmpConstants[plane];
84  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
85  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, det_prop.Efield());
86  }
87 
88  //------------------------------------------------------------------------------------//
89  // Functions to calculate the dEdX based on the AMPLITUDE of the pulse with EField
90  // ----------------------------------------------------------------------------------//
91  double CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
92  detinfo::DetectorPropertiesData const& det_prop,
93  recob::Hit const& hit,
94  double const pitch,
95  double const T0,
96  double const EField,
97  double const phi) const
98  {
99  return dEdx_AMP(clock_data,
100  det_prop,
101  hit.PeakAmplitude() / pitch,
102  hit.PeakTime(),
103  hit.WireID().Plane,
104  T0,
105  EField,
106  phi);
107  }
108 
109  // ----------------------------------------------------------------------------------//
110  double CalorimetryAlg::dEdx_AMP(detinfo::DetectorClocksData const& clock_data,
111  detinfo::DetectorPropertiesData const& det_prop,
112  double const dQdx,
113  double const time,
114  unsigned int const plane,
115  double const T0,
116  double const EField,
117  double const phi) const
118  {
119  double const fADCtoEl = fCalAmpConstants[plane];
120  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
121  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, EField, phi);
122  }
123 
124  //------------------------------------------------------------------------------------//
125  // Functions to calculate the dEdX based on the AREA of the pulse
126  // ----------------------------------------------------------------------------------//
127  double CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
128  detinfo::DetectorPropertiesData const& det_prop,
129  recob::Hit const& hit,
130  double const pitch,
131  double const T0) const
132  {
133  return dEdx_AREA(clock_data,
134  det_prop,
135  hit.Integral() / pitch,
136  hit.PeakTime(),
137  hit.WireID().Plane,
138  T0,
139  det_prop.Efield());
140  }
141 
142  // ----------------------------------------------------------------------------------//
143  double CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
144  detinfo::DetectorPropertiesData const& det_prop,
145  double const dQdx,
146  double const time,
147  unsigned int const plane,
148  double const T0) const
149  {
150  double const fADCtoEl = fCalAreaConstants[plane];
151  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
152  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, det_prop.Efield());
153  }
154 
155  //------------------------------------------------------------------------------------//
156  // Functions to calculate the dEdX based on the AREA of the pulse with EField
157  // ----------------------------------------------------------------------------------//
158  double CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
159  detinfo::DetectorPropertiesData const& det_prop,
160  recob::Hit const& hit,
161  double const pitch,
162  double const T0,
163  double const EField,
164  double const phi) const
165  {
166  return dEdx_AREA(clock_data,
167  det_prop,
168  hit.Integral() / pitch,
169  hit.PeakTime(),
170  hit.WireID().Plane,
171  T0,
172  EField,
173  phi);
174  }
175 
176  // ----------------------------------------------------------------------------------//
177  double CalorimetryAlg::dEdx_AREA(detinfo::DetectorClocksData const& clock_data,
178  detinfo::DetectorPropertiesData const& det_prop,
179  double const dQdx,
180  double const time,
181  unsigned int const plane,
182  double const T0,
183  double const EField,
184  double const phi) const
185  {
186  double const fADCtoEl = fCalAreaConstants[plane];
187  double const dQdx_e = dQdx / fADCtoEl; // Conversion from ADC/cm to e/cm
188  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, EField, phi);
189  }
190 
191  // Apply Lifetime and recombination correction.
192  double CalorimetryAlg::dEdx_from_dQdx_e(detinfo::DetectorClocksData const& clock_data,
193  detinfo::DetectorPropertiesData const& det_prop,
194  double dQdx_e,
195  double const time,
196  double const T0) const
197  {
198  return dEdx_from_dQdx_e(clock_data, det_prop, dQdx_e, time, T0, det_prop.Efield());
199  }
200  double CalorimetryAlg::dEdx_from_dQdx_e(detinfo::DetectorClocksData const& clock_data,
201  detinfo::DetectorPropertiesData const& det_prop,
202  double dQdx_e,
203  double const time,
204  double const T0,
205  double const EField,
206  double const phi) const
207  {
208  if (fDoLifeTimeCorrection) {
209  dQdx_e *= LifetimeCorrection(clock_data, det_prop, time, T0); // (dQdx_e in e/cm)
210  }
211 
212  if (fUseModBox) { return ModBoxCorrection(dQdx_e, phi, det_prop.Density(), EField); }
213 
214  return BirksCorrection(dQdx_e, phi, det_prop.Density(), EField);
215  }
216 
217  //------------------------------------------------------------------------------------//
218  // for the time being copying from Calorimetry.cxx - should be decided where
219  // to keep it.
220  // ----------------------------------------------------------------------------------//
221  double calo::CalorimetryAlg::LifetimeCorrection(detinfo::DetectorClocksData const& clock_data,
222  detinfo::DetectorPropertiesData const& det_prop,
223  double const time,
224  double const T0) const
225  {
226  float const t = time - trigger_offset(clock_data);
227  double const timetick = sampling_rate(clock_data) * 1.e-3; // time sample in microsec
228  double const adjusted_time = t * timetick - T0 * 1e-3; // (in microsec)
229 
230  assert(fLifeTimeForm < 2);
231  if (fLifeTimeForm == 0) {
232  // Exponential form
233  double const tau = det_prop.ElectronLifetime();
234  return exp(adjusted_time / tau);
235  }
236 
237  // Exponential+constant form
238  auto const& elifetime_provider =
239  art::ServiceHandle<lariov::ElectronLifetimeService const>()->GetProvider();
240  return elifetime_provider.Lifetime(adjusted_time);
241  }
242 
243  // Reombination corrections
244 
245  // Modified box: allow for a general behavior of beta from phi, the angle between the track
246  // and the electric field
247  double calo::CalorimetryAlg::ModBoxCorrection(double dQdx,
248  double phi,
249  double rho,
250  double E_field) const
251  {
252  // Modified Box model correction has better behavior than the Birks
253  // correction at high values of dQ/dx.
254  constexpr double Wion = 1000. / util::kGeVToElectrons; // 23.6 eV = 1e, Wion in MeV/e
255  double const Beta = fModBoxBF.Eval(phi) / (rho * E_field);
256  double const Alpha = fModBoxA;
257  double const dEdx = (exp(Beta * Wion * dQdx) - Alpha) / Beta;
258 
259  return dEdx;
260  }
261 
262  double calo::CalorimetryAlg::BirksCorrection(double dQdx,
263  double phi,
264  double rho,
265  double E_field) const
266  {
267  // from: S.Amoruso et al., NIM A 523 (2004) 275
268 
269  double A = fBirksA;
270  double K = fBirksKF.Eval(phi); // in KV/cm*(g/cm^2)/MeV
271  constexpr double Wion = 1000. / util::kGeVToElectrons; // 23.6 eV = 1e, Wion in MeV/e
272  K /= rho; // KV/MeV
273  double const dEdx = dQdx / (A / Wion - K / E_field * dQdx); // MeV/cm
274 
275  return dEdx;
276  }
277 
278 } // namespace
calo::CalorimetryAlg::fCalAreaConstants
std::vector< double > const fCalAreaConstants
Definition: CalorimetryAlg.h:186
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:244
recob::Hit::WireID
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
recob::Hit::PeakAmplitude
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:232
detinfo::DetectorPropertiesData::Efield
double Efield(unsigned int planegap=0) const
kV/cm
Definition: DetectorPropertiesData.cc:16
util::kGeVToElectrons
constexpr double kGeVToElectrons
23.6eV per ion pair, 1e9 eV/GeV
Definition: PhysicalConstants.h:55
calo::CalorimetryAlg::dEdx_AMP
double dEdx_AMP(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
Definition: CalorimetryAlg.cxx:60
util::kModBoxB
constexpr double kModBoxB
Modified Box Beta in g/(MeV cm²)*kV/cm.
Definition: PhysicalConstants.h:50
calo::CalorimetryAlg::fDoLifeTimeCorrection
bool const fDoLifeTimeCorrection
Definition: CalorimetryAlg.h:189
detinfo::DetectorPropertiesData::Density
double Density(double temperature=0.) const
Returns argon density at a given temperature.
Definition: DetectorPropertiesData.cc:67
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
calo::CalorimetryAlg::fCalAmpConstants
std::vector< double > const fCalAmpConstants
Definition: CalorimetryAlg.h:185
calo::CalorimetryAlg::BirksCorrection
double BirksCorrection(double dQdx, double phi, double rho, double E_field) const
Definition: CalorimetryAlg.cxx:262
geo_types.h
Definition of data types for geometry description.
tca::dEdx
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
calo::CalorimetryAlg::CalorimetryAlg
CalorimetryAlg(const fhicl::ParameterSet &pset)
Definition: CalorimetryAlg.h:89
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:282
calo::CalorimetryAlg::dEdx_from_dQdx_e
double dEdx_from_dQdx_e(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double dQdx_e, double time, double T0=0) const
Definition: CalorimetryAlg.cxx:192
util::kRecombk
constexpr double kRecombk
Definition: PhysicalConstants.h:36
calo::CalorimetryAlg::ModBoxCorrection
double ModBoxCorrection(double dQdx, double phi, double rho, double E_field) const
Definition: CalorimetryAlg.cxx:247
calo::CalorimetryAlg::dEdx_AREA
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
Definition: CalorimetryAlg.cxx:127
detinfo::trigger_offset
int trigger_offset(DetectorClocksData const &data)
Definition: DetectorClocksData.h:565
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
calo::CalorimetryAlg::LifetimeCorrection
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Definition: CalorimetryAlg.cxx:221
e
Float_t e
Definition: plot.C:35
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:575
calo::CalorimetryAlg::Config::CalAmpConstants
fhicl::Sequence< double > CalAmpConstants
Definition: CalorimetryAlg.h:43
calo
calorimetry
Definition: CaloChecker_module.cc:43
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33