dEdxTable.cxx
Go to the documentation of this file.
1 ///
2 /// \file dEdxTables.cxx
3 /// \brief Implement energy loss tables using GEANT4
4 /// \author messier@indiana.edu
5 ///
6 /// These tables were calculated using novasoft on June 16, 2015 using GEANT4:
7 ///
8 /// Geant4 version Name: geant4-09-06-patch-04a (04-February-2015)
9 /// Copyright : Geant4 Collaboration
10 /// Reference : NIM A 506 (2003), 250-303
11 /// WWW : http://cern.ch/geant4
12 ///
13 /// In particular, they come from G4EmCalculator::ComputeTotalDEDX()
14 /// for muons in various NOvA materials. This is implemented into the
15 /// g4nova/PrintTables functions.
16 ///
18 #include <string>
19 #include <cmath>
20 #include <iostream>
21 #include <algorithm>
22 using namespace bpfit;
23 //
24 // Muon mass used by Groom, Mokhov, Striganov to compute their tables
25 //
26 static const double gsMmu = 105.6583568E-3;
27 static const double log_gsMmu = log(gsMmu);
28 
29 //......................................................................
30 
31 dEdxTable::dEdxTable(const char* mat) :
32  fSz(0),
33  fT(0),
34  fRho(0)
35 {
36  static const std::string scint ("Scintillator");
37  static const std::string pvc ("PVC");
38  static const std::string glue ("Glue");
39  static const std::string steel ("Steel");
40  static const std::string vac ("Vacuum");
41 
42  if (scint==mat) this->InitScintillator();
43  else if (pvc ==mat) this->InitPVC();
44  else if (steel==mat) this->InitSteel();
45  else if (glue ==mat) this->InitGlue();
46  else if (vac ==mat) this->InitVac();
47  else {
48  static unsigned int nwarn = 10;
49  if (nwarn>0) {
50  --nwarn;
51  std::cerr << __FILE__ << ":" << __LINE__
52  << " Have no table for material " << mat
53  << " Defaulting to Vacuum.";
54  if (nwarn==0) {
55  std::cerr << " Last warning.";
56  }
58  this->InitVac();
59  }
60  }
61 }
62 
63 //......................................................................
64 
65 double dEdxTable::Interp(const double bg,
66  const double logT[],
67  const double logdedx[],
68  const double slope[],
69  unsigned int imax) const
70 {
71  //
72  // Compute equivalent muon kinetic energy. Ie. a muon with
73  // the same beta*gamma
74  //
75  const double x = log_gsMmu + log(bg);
76 
77  // Very very often consecutive values of x are close enough to
78  // be between the same array values, so cache i1 and avoid calling
79  // std::upper_bound most of the time.
80  static unsigned int i1=1;
81  if(i1 >= 1 && i1 < imax && logT[i1-1] < x && logT[i1] >= x)/* no change */;
82  else if (x<=logT[0]) i1=1;
83  else if (x>=logT[imax-1]) i1=imax-1;
84  else i1 = std::upper_bound(logT, logT+imax, x) - logT;
85 
86  const double x0 = logT[i1-1];
87  const double y0 = logdedx[i1-1];
88 
89  const double y = y0 + (x-x0)*slope[i1];
90 
91  return exp(y);
92 }
93 
94 //......................................................................
95 //
96 // Initialize tables for PVC
97 //
99 {
100  fSz = 90;
101  static const double T[] = {
102  1e-05, 1.26219e-05, 1.59312e-05, 2.01081e-05, 2.53802e-05, 3.20345e-05, 4.04335e-05, 5.10347e-05, 6.44153e-05, 8.13041e-05,
103  0.000102621, 0.000129527, 0.000163487, 0.000206351, 0.000260454, 0.000328741, 0.000414933, 0.000523722, 0.000661035, 0.00083435,
104  0.00105311, 0.00132922, 0.00167772, 0.00211759, 0.0026728, 0.00337357, 0.00425808, 0.00537449, 0.00678361, 0.00856218,
105  0.0108071, 0.0136405, 0.0172169, 0.0217309, 0.0274285, 0.0346199, 0.0436968, 0.0551535, 0.069614, 0.0878659,
106  0.110903, 0.13998, 0.176681, 0.223005, 0.281474, 0.355273, 0.44842, 0.56599, 0.714385, 0.901688,
107  1.1381, 1.43649, 1.81312, 2.2885, 2.88851, 3.64584, 4.60173, 5.80825, 7.33109, 9.2532,
108  11.6793, 14.7414, 18.6064, 23.4848, 29.6422, 37.414, 47.2234, 59.6048, 75.2323, 94.9572,
109  119.854, 151.278, 190.941, 241.003, 304.191, 383.946, 484.611, 611.67, 772.041, 974.46,
110  1229.95, 1552.43, 1959.45, 2473.2, 3121.63, 3940.09, 4973.12, 6277.01, 7922.76, 10000
111  };
112  static const double dEdx[] = {
113  433.154, 434.907, 426.57, 410.17, 388.417, 361.638, 330.62, 297.715, 264.919, 233.809,
114  204.698, 178.096, 154.171, 132.923, 114.048, 97.2441, 82.3936, 69.5649, 58.4527, 48.9138,
115  40.8327, 34.0107, 28.2895, 23.4964, 19.5062, 16.1859, 13.4461, 11.1866, 9.32543, 7.79674,
116  6.54546, 5.52379, 4.69278, 4.01968, 3.47735, 3.04335, 2.69917, 2.42961, 2.22215, 2.06648,
117  1.94627, 1.85994, 1.80319, 1.7705, 1.75712, 1.75895, 1.77251, 1.79485, 1.82357, 1.85671,
118  1.89416, 1.93357, 1.97392, 2.01459, 2.0552, 2.09546, 2.13529, 2.17464, 2.21357, 2.25228,
119  2.29107, 2.33036, 2.37076, 2.41308, 2.45836, 2.50798, 2.5637, 2.62775, 2.70301, 2.79306,
120  2.90248, 3.03702, 3.2039, 3.41229, 3.67378, 4.00307, 4.41873, 4.94427, 5.60942, 6.45176,
121  7.51878, 8.87046, 10.5827, 12.7509, 15.4958, 18.9701, 23.3661, 28.9267, 35.9586, 44.849
122  };
123 
124  fRho = 1.49;
125 
126  static double * logT = new double[fSz], * logdEdx = new double[fSz], *slope = new double[fSz];
127  for(unsigned int i = 0; i < fSz; i++){
128  logT[i] = log(T[i]);
129  logdEdx[i] = log(fRho*dEdx[i]*1e-3);
130  if(i == 0) slope[i] = 0;
131  else slope[i] = (logdEdx[i] - logdEdx[i-1])/(logT[i] - logT[i-1]);
132  }
133 
134  flogT = logT;
135  flogdEdx = logdEdx;
136  fSlope = slope;
137 
138  fT = T;
139 }
140 
141 //......................................................................
142 ///
143 /// Initialize tables for steel
144 ///
146 {
147  fSz = 90;
148  static const double T[] = {
149  1e-05, 1.26219e-05, 1.59312e-05, 2.01081e-05, 2.53802e-05, 3.20345e-05, 4.04335e-05, 5.10347e-05, 6.44153e-05, 8.13041e-05,
150  0.000102621, 0.000129527, 0.000163487, 0.000206351, 0.000260454, 0.000328741, 0.000414933, 0.000523722, 0.000661035, 0.00083435,
151  0.00105311, 0.00132922, 0.00167772, 0.00211759, 0.0026728, 0.00337357, 0.00425808, 0.00537449, 0.00678361, 0.00856218,
152  0.0108071, 0.0136405, 0.0172169, 0.0217309, 0.0274285, 0.0346199, 0.0436968, 0.0551535, 0.069614, 0.0878659,
153  0.110903, 0.13998, 0.176681, 0.223005, 0.281474, 0.355273, 0.44842, 0.56599, 0.714385, 0.901688,
154  1.1381, 1.43649, 1.81312, 2.2885, 2.88851, 3.64584, 4.60173, 5.80825, 7.33109, 9.2532,
155  11.6793, 14.7414, 18.6064, 23.4848, 29.6422, 37.414, 47.2234, 59.6048, 75.2323, 94.9572,
156  119.854, 151.278, 190.941, 241.003, 304.191, 383.946, 484.611, 611.67, 772.041, 974.46,
157  1229.95, 1552.43, 1959.45, 2473.2, 3121.63, 3940.09, 4973.12, 6277.01, 7922.76, 10000
158  };
159  static const double dEdx[] = {
160  169.215, 175.49, 178.929, 179.971, 178.078, 173.061, 165.115, 154.965, 143.604, 131.176,
161  118.637, 106.228, 94.1998, 82.437, 73.2226, 63.9816, 55.4551, 47.6246, 40.7865, 34.7465,
162  29.3611, 24.5292, 20.8287, 17.5669, 14.7612, 12.3854, 10.3773, 8.70296, 7.30404, 6.13792,
163  5.17274, 4.3853, 3.74771, 3.22531, 2.80237, 2.46185, 2.191, 1.97863, 1.81539, 1.69346,
164  1.60099, 1.53594, 1.49475, 1.47312, 1.46727, 1.47392, 1.49021, 1.51375, 1.54253, 1.5749,
165  1.61147, 1.64925, 1.68779, 1.72662, 1.76546, 1.80417, 1.84278, 1.88139, 1.92029, 1.95989,
166  2.00082, 2.04391, 2.09027, 2.14139, 2.19917, 2.2661, 2.34535, 2.44098, 2.55816, 2.70343,
167  2.88513, 3.11374, 3.40254, 3.76834, 4.23242, 4.82172, 5.57032, 6.5213, 7.72916, 9.26274,
168  11.209, 13.6779, 16.8082, 20.7745, 25.7982, 32.1583, 40.207, 50.3888, 63.2649, 79.5436
169  };
170 
171  fRho = 7.87;
172 
173  static double * logT = new double[fSz], * logdEdx = new double[fSz], *slope = new double[fSz];
174  for(unsigned int i = 0; i < fSz; i++){
175  logT[i] = log(T[i]);
176  logdEdx[i] = log(fRho*dEdx[i]*1e-3);
177  if(i == 0) slope[i] = 0;
178  else slope[i] = (logdEdx[i] - logdEdx[i-1])/(logT[i] - logT[i-1]);
179  }
180 
181  flogT = logT;
182  flogdEdx = logdEdx;
183  fSlope = slope;
184 
185  fT = T;
186 }
187 
188 //......................................................................
189 //
190 // Vaccuum is a valid GEANT material. Set things very close to
191 // zero. Keep the tables small for sake of efficiency
192 //
194 {
195  fSz = 2;
196  static const double T[] = { 10.0E-3, 800.0 };
197  static const double dEdx[] = { 1.0E-6, 1.0E-6};
198 
199  fRho = 1.0E-6;
200 
201  static double * logT = new double[fSz], * logdEdx = new double[fSz], * slope = new double[fSz];
202  for(unsigned int i = 0; i < fSz; i++){
203  logT[i] = log(T[i]);
204  logdEdx[i] = log(fRho*dEdx[i]*1e-3);
205  if(i == 0) slope[i] = 0;
206  else slope[i] = (logdEdx[i] - logdEdx[i-1])/(logT[i] - logT[i-1]);
207  }
208 
209  flogT = logT;
210  flogdEdx = logdEdx;
211  fSlope = slope;
212 
213  fT = T;
214 }
215 
216 //......................................................................
217 //
218 // Tables for our scintillator
219 //
221 {
222  fSz = 90;
223  static const double T[] = {
224  1e-05, 1.26219e-05, 1.59312e-05, 2.01081e-05, 2.53802e-05, 3.20345e-05, 4.04335e-05, 5.10347e-05, 6.44153e-05, 8.13041e-05,
225  0.000102621, 0.000129527, 0.000163487, 0.000206351, 0.000260454, 0.000328741, 0.000414933, 0.000523722, 0.000661035, 0.00083435,
226  0.00105311, 0.00132922, 0.00167772, 0.00211759, 0.0026728, 0.00337357, 0.00425808, 0.00537449, 0.00678361, 0.00856218,
227  0.0108071, 0.0136405, 0.0172169, 0.0217309, 0.0274285, 0.0346199, 0.0436968, 0.0551535, 0.069614, 0.0878659,
228  0.110903, 0.13998, 0.176681, 0.223005, 0.281474, 0.355273, 0.44842, 0.56599, 0.714385, 0.901688,
229  1.1381, 1.43649, 1.81312, 2.2885, 2.88851, 3.64584, 4.60173, 5.80825, 7.33109, 9.2532,
230  11.6793, 14.7414, 18.6064, 23.4848, 29.6422, 37.414, 47.2234, 59.6048, 75.2323, 94.9572,
231  119.854, 151.278, 190.941, 241.003, 304.191, 383.946, 484.611, 611.67, 772.041, 974.46,
232  1229.95, 1552.43, 1959.45, 2473.2, 3121.63, 3940.09, 4973.12, 6277.01, 7922.76, 10000
233  };
234  static const double dEdx[] = {
235  662.297, 660.884, 643.851, 613.531, 576.982, 532.716, 482.603, 430.261, 378.82, 330.916,
236  286.679, 246.807, 211.44, 180.382, 152.86, 128.949, 108.129, 90.3863, 75.3084, 62.5557,
237  51.8743, 42.9576, 35.5366, 29.3797, 24.2861, 20.0831, 16.6232, 13.7815, 11.4525, 9.5482,
238  7.99474, 6.73087, 5.70578, 4.87751, 4.21146, 3.67924, 3.25753, 2.92727, 2.67283, 2.4814,
239  2.33543, 2.23014, 2.15962, 2.11726, 2.09735, 2.09496, 2.10592, 2.12671, 2.15449, 2.18698,
240  2.22365, 2.26246, 2.30202, 2.3417, 2.38109, 2.41997, 2.45825, 2.49593, 2.53309, 2.5699,
241  2.60664, 2.64354, 2.68087, 2.71902, 2.75851, 2.80006, 2.8446, 2.89335, 2.94783, 3.01003,
242  3.08242, 3.16815, 3.2712, 3.39662, 3.55079, 3.74178, 3.97982, 4.27786, 4.65228, 5.12381,
243  5.71864, 6.46986, 7.41943, 8.61997, 10.1383, 12.0585, 14.487, 17.558, 21.4409, 26.3498
244  };
245 
246  fRho = 0.859;
247 
248  static double * logT = new double[fSz], * logdEdx = new double[fSz], * slope = new double[fSz];
249  for(unsigned int i = 0; i < fSz; i++){
250  logT[i] = log(T[i]);
251  logdEdx[i] = log(fRho*dEdx[i]*1e-3);
252  if(i == 0) slope[i] = 0;
253  else slope[i] = (logdEdx[i] - logdEdx[i-1])/(logT[i] - logT[i-1]);
254  }
255 
256  flogT = logT;
257  flogdEdx = logdEdx;
258  fSlope = slope;
259 
260  fT = T;
261 }
262 
263 //......................................................................
264 //
265 // Tables for glue
266 ///
268 {
269  fSz = 90;
270  static const double T[] = {
271  1e-05, 1.26219e-05, 1.59312e-05, 2.01081e-05, 2.53802e-05, 3.20345e-05, 4.04335e-05, 5.10347e-05, 6.44153e-05, 8.13041e-05,
272  0.000102621, 0.000129527, 0.000163487, 0.000206351, 0.000260454, 0.000328741, 0.000414933, 0.000523722, 0.000661035, 0.00083435,
273  0.00105311, 0.00132922, 0.00167772, 0.00211759, 0.0026728, 0.00337357, 0.00425808, 0.00537449, 0.00678361, 0.00856218,
274  0.0108071, 0.0136405, 0.0172169, 0.0217309, 0.0274285, 0.0346199, 0.0436968, 0.0551535, 0.069614, 0.0878659,
275  0.110903, 0.13998, 0.176681, 0.223005, 0.281474, 0.355273, 0.44842, 0.56599, 0.714385, 0.901688,
276  1.1381, 1.43649, 1.81312, 2.2885, 2.88851, 3.64584, 4.60173, 5.80825, 7.33109, 9.2532,
277  11.6793, 14.7414, 18.6064, 23.4848, 29.6422, 37.414, 47.2234, 59.6048, 75.2323, 94.9572,
278  119.854, 151.278, 190.941, 241.003, 304.191, 383.946, 484.611, 611.67, 772.041, 974.46,
279  1229.95, 1552.43, 1959.45, 2473.2, 3121.63, 3940.09, 4973.12, 6277.01, 7922.76, 10000
280  };
281  static const double dEdx[] = {
282  552.504, 554.049, 543.558, 522.08, 493.313, 458.006, 417.645, 374.441, 331.402, 290.778,
283  252.868, 218.495, 187.771, 160.669, 136.674, 115.537, 97.1704, 81.5452, 68.139, 56.7364,
284  47.1509, 39.1196, 32.413, 26.8358, 22.2121, 18.389, 15.2373, 12.6449, 10.5176, 8.77596,
285  7.35375, 6.19557, 5.25544, 4.49529, 3.88368, 3.39474, 3.00723, 2.70376, 2.47003, 2.29432,
286  2.16365, 2.07074, 2.00943, 1.9738, 1.95865, 1.9595, 1.97251, 1.99444, 2.02264, 2.055,
287  2.09112, 2.12897, 2.16732, 2.2056, 2.24345, 2.2807, 2.31731, 2.35331, 2.38884, 2.42412,
288  2.45947, 2.4952, 2.53163, 2.56921, 2.60856, 2.65051, 2.69613, 2.74682, 2.80437, 2.87104,
289  2.9497, 3.044, 3.15854, 3.29917, 3.47328, 3.69021, 3.96182, 4.30309, 4.73299, 5.27551,
290  5.96095, 6.82761, 7.92401, 9.31105, 11.0659, 13.2861, 16.0945, 19.6464, 24.1377, 29.816
291  };
292 
293  fRho = 1.34;
294 
295  static double * logT = new double[fSz], * logdEdx = new double[fSz], * slope = new double[fSz];
296  for(unsigned int i = 0; i < fSz; i++){
297  logT[i] = log(T[i]);
298  logdEdx[i] = log(fRho*dEdx[i]*1e-3);
299  if(i == 0) slope[i] = 0;
300  else slope[i] = (logdEdx[i] - logdEdx[i-1])/(logT[i] - logT[i-1]);
301  }
302 
303  flogT = logT;
304  flogdEdx = logdEdx;
305  fSlope = slope;
306 
307  fT = T;
308 }
309 
310 //......................................................................
311 
312 double dEdxTable::operator()(double T, double m) const
313 {
314  if (fT==0) return 0.0;
315  return this->Interp(T/m, flogT, flogdEdx, fSlope, fSz);
316 }
317 
318 ////////////////////////////////////////////////////////////////////////
const double * flogT
Log of muon kinetic energy.
Definition: dEdxTable.h:48
const double * fT
Muon kinetic energy.
Definition: dEdxTable.h:47
void InitScintillator()
Definition: dEdxTable.cxx:220
Definition: scint.py:1
OStream cerr
Definition: OStream.cxx:7
Tables of mean energy loss in NOvA materials.
Definition: pvc.py:1
double operator()(double T, double m) const
Definition: dEdxTable.cxx:312
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const double * flogdEdx
Log of muon energy loss.
Definition: dEdxTable.h:49
static const double gsMmu
Definition: dEdxTable.cxx:26
Definition: glue.py:1
static const double log_gsMmu
Definition: dEdxTable.cxx:27
dEdxTable(const char *mat)
Build dEdx table for named material.
Definition: dEdxTable.cxx:31
double Interp(const double bg, const double logT[], const double dedx[], const double slope[], unsigned int imax) const
Definition: dEdxTable.cxx:65
double T
Definition: Xdiff_gwt.C:5
const double * fSlope
slope of the logs of energy vs. energy loss
Definition: dEdxTable.h:50
Float_t e
Definition: plot.C:35
unsigned int fSz
Table size.
Definition: dEdxTable.h:46
Eigen::MatrixXd mat
double fRho
Material density.
Definition: dEdxTable.h:51
enum BeamMode string