BetheBloch.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file BetheBloch.cxx
3 // \brief Bethe Bloch dE/dx from tables or analytic formula
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
8 
9 #include "Utilities/func/MathUtil.h"
10 
11 #include <algorithm>
12 #include <cassert>
13 #include <cstdio>
14 
15 #include "cetlib/search_path.h"
16 
17 #include "TF1.h"
18 #include "TGraph.h"
19 
20 namespace calib
21 {
22  // Generic physical constants:
23 
24  // From PDG
25  const double Me = 510.998910e3; // eV
26  const double Mmu = 105.658367e6; // eV
27  // From PDG "Passage of particles through matter"
28  const double K = 0.307075; // MeV g^-1 cm^2
29  const double j = 0.200;
30 
31 
32  //......................................................................
33  TGraph* IBetheBloch::GetdEdxGraph(double localDensity) const
34  {
35  TGraph* ret = new TGraph;
36  // 1MeV to 100GeV
37  for(double T = 1; T < 1e5; T *= 1.05){
38  ret->SetPoint(ret->GetN(), T, dEdx(T)*localDensity);
39  }
40 
41  return ret;
42  }
43 
44  //......................................................................
45  TGraph* IBetheBloch::GetMPVGraph(double localDensity, double dist) const
46  {
47  TGraph* ret = new TGraph;
48  // 1MeV to 100GeV
49  for(double T = 1; T < 1e5; T *= 1.05){
50  ret->SetPoint(ret->GetN(), T,
51  MPV(T, dist, localDensity)*localDensity/dist);
52  }
53 
54  return ret;
55  }
56 
57 
58  //......................................................................
60  {
61  char junkStr[1024];
62  double junk;
63 
64  // Load the points from the file
66  cet::search_path sp("FW_SEARCH_PATH");
67  switch(mat){
68  case kPolyethylene:
69  // http://pdg.lbl.gov/2010/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_221.dat
70  sp.find_file("Calibration/func/muonloss_221.dat", fname);
71  break;
72  default:
73  assert(0 && "Unknown material");
74  }
75  FILE* f = fopen(fname.c_str(), "r");
76  assert(f);
77 
78  while(!feof(f)){
79  double T, dEdx;
80  int ret = fscanf(f, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
81  &T, &junk, &junk, &junk, &junk, &junk, &junk,
82  &dEdx, &junk, &junk, &junk);
83  if(ret == 11){
84  fT_pts.push_back(T_pt(T, dEdx));
85  }
86  else{
87  // We couldn't parse it properly. Either preamble or a special line.
88  fgets(junkStr, 1024, f);
89  }
90  }
91  fclose(f);
92  }
93 
95  {
96  return pt.T < T;
97  }
98  //......................................................................
99  double BetheBlochTables::dEdx(double T) const
100  {
101  typedef std::vector<T_pt>::const_iterator it_t;
102 
103  // Find the point after T
104  it_t itNext = std::lower_bound(fT_pts.begin(), fT_pts.end(), T, ltTPt);
105  if(itNext == fT_pts.end()) --itNext;
106  if(itNext == fT_pts.begin()) ++itNext;
107  // And the point before
108  it_t itPrev = itNext; --itPrev;
109 
110  // Interpolate
111  const double T0 = itPrev->T; const double y0 = itPrev->dEdx;
112  const double T1 = itNext->T; const double y1 = itNext->dEdx;
113  return ((T1-T)*y0+(T-T0)*y1)/(T1-T0);
114  }
115 
116  //......................................................................
117  double BetheBlochTables::MIP() const
118  {
119  double mip = 1e10;
120  for(unsigned int n = 0; n < fT_pts.size(); ++n){
121  if(fT_pts[n].dEdx < mip) mip = fT_pts[n].dEdx;
122  }
123  return mip;
124  }
125 
126 
127  //......................................................................
129  {
130  switch(mat){
131  case kPolyethylene:
132  // Numbers from the header of http://pdg.lbl.gov/2010/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_221.dat
133  ZA = 0.57034; // Z/A
134  I = 57.4; // eV
135  a = 0.1211;
136  k = 3.4292;
137  x0 = 0.1489;
138  x1 = 2.5296;
139  Cbar = 3.0563;
140  break;
141  default:
142  assert(0 && "Unkown material");
143  }
144  }
145 
146  //......................................................................
147  double BetheBlochAnalytic::dEdx(double T) const
148  {
149  using util::sqr;
150 
151  T *= 1e6; // Convert to eV
152 
153  double beta, gamma;
154  BetaGamma(T, beta, gamma);
155 
156  const double f = 2*Me*sqr(beta*gamma);
157 
158  const double Tmax = f/(1+2*gamma*Me/Mmu+sqr(Me/Mmu));
159 
160  const double delta = Delta(beta*gamma);
161 
162  const double dEdx = K*ZA/sqr(beta)*(log(f*Tmax/(sqr(I)))/2-sqr(beta)-delta/2);
163 
164  return dEdx;
165  }
166 
167  //......................................................................
168  double BetheBlochAnalytic::MPV(double T, double dist,
169  double localDensity) const
170  {
171  using util::sqr;
172 
173  const double x = dist*localDensity;
174 
175  T *= 1e6; // Convert to eV
176 
177  double beta, gamma;
178  BetaGamma(T, beta, gamma);
179 
180  const double f = 2*Me*sqr(beta*gamma);
181 
182  const double delta = Delta(beta*gamma);
183 
184  const double xi = K/2*ZA/sqr(beta)*x;
185  const double mpv = xi*(log(f/I)+log(1e6*xi/I)+j-sqr(beta)-delta);
186 
187  return mpv;
188  }
189 
190  //......................................................................
191  double BetheBlochAnalytic::MIP() const
192  {
193  double mip = 1e20;
194  // Between 100 MeV and 1GeV, with 0.1% accuracy
195  for(double T = 1e2; T < 1e3; T *= 1.001){
196  const double m = dEdx(T);
197  if(m < mip) mip = m;
198  }
199  return mip;
200  }
201 
202  //......................................................................
204  BetaGamma(double T, double& beta, double& gamma) const
205  {
206  gamma = (T+Mmu)/Mmu;
207  beta = sqrt(1-1/(gamma*gamma));
208  }
209 
210  //......................................................................
211  double BetheBlochAnalytic::Delta(double betagamma) const
212  {
213  double delta = 0;
214  const double x = log10(betagamma);
215  if(x > x0){
216  delta += 2*log(10)*x-Cbar;
217  if(x < x1) delta += a*pow(x1-x, k);
218  }
219  return delta;
220  }
221 
222 }
BetheBlochTables(Material mat)
Definition: BetheBloch.cxx:59
virtual double MIP() const
Definition: BetheBloch.cxx:117
Float_t y1[n_points_granero]
Definition: compare.C:5
double delta
Definition: runWimpSim.h:98
virtual double dEdx(double T) const
Definition: BetheBloch.cxx:147
virtual double MPV(double T, double dist, double localDensity) const
Definition: BetheBloch.h:31
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
const double T0
virtual double dEdx(double T) const
Definition: BetheBloch.cxx:99
Double_t beta
virtual double dEdx(double T) const =0
std::string find_file(std::string const &filename) const
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
const double Mmu
Definition: BetheBloch.cxx:26
void BetaGamma(double T, double &beta, double &gamma) const
Relativity.
Definition: BetheBloch.cxx:204
double dist
Definition: runWimpSim.h:113
bool ltTPt(BetheBlochTables::T_pt pt, double T)
Definition: BetheBloch.cxx:94
fclose(fg1)
CDPStorage service.
TGraph * GetMPVGraph(double localDensity, double dist) const
Definition: BetheBloch.cxx:45
virtual double MPV(double T, double dist, double localDensity) const
Definition: BetheBloch.cxx:168
const double a
TGraph * GetdEdxGraph(double localDensity) const
Definition: BetheBloch.cxx:33
Material
Definition: BetheBloch.h:16
const double j
Definition: BetheBloch.cxx:29
T sqr(T x)
Function to perform, and table to cache, timing fits to ADC values.
Definition: ADCShapeFit.cxx:23
const double K
Definition: BetheBloch.cxx:28
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual double MIP() const
Definition: BetheBloch.cxx:191
T log10(T number)
Definition: d0nt_math.hpp:120
assert(nhit_max >=nhit_nbins)
double T
Definition: Xdiff_gwt.C:5
const double Me
Definition: BetheBloch.cxx:25
double Delta(double betagamma) const
Density term.
Definition: BetheBloch.cxx:211
BetheBlochAnalytic(Material mat)
Definition: BetheBloch.cxx:128
Eigen::MatrixXd mat