Tracking/dEdx/bethe_bloch.h

0001 #ifndef BETHE_BLOCH_H
0002 #define BETHE_BLOCH_H
0003
0004 #include "TMath.h"
0005
0006 namespace dedx_constants
0007 {
0008   // hadron masses
0009   constexpr double m_pi = 0.1396; // GeV
0010   constexpr double m_K = 0.4937; // GeV
0011   constexpr double m_p = 0.9382; // GeV
0012   constexpr double m_d = 1.876; // GeV
0013
0014   // electron mass [eV]
0015   constexpr double m_e = 511e3;
0016
0017   // TPC gas fractions
0018   constexpr double ar_frac = 0.75;
0019   constexpr double cf4_frac = 0.2;
0020   constexpr double isobutane_frac = 0.05;
0021
0022   // Mean excitation [src: W. Blum, W. Riegler, L. Rolandi, "Particle Detection with Drift Chambers"]
0023   constexpr double ar_I = 188; // eV
0024   constexpr double cf4_I = 115; // eV
0025   constexpr double isobutane_I = 48.3; // eV
0026
0027   // Mean excitation of mixture approximated using Bragg additivity rule
0028   constexpr double sphenix_I = ar_frac*ar_I + cf4_frac*cf4_I + isobutane_frac*isobutane_I;
0029 }
0030
0031 // Bethe-Bloch fit function, vs. betagamma
0032 // A = normalization constant, equal to (ADC conversion)*4pi*n*Z^2*e^4/(m_e*c^2*4pi*epsilon_0^2)
0033 // B = A*(ln(2*m_e/I)-1) - (zero-suppression loss factor)
0034 const double bethe_bloch_new(const double betagamma, const double A, const double B, const double C)
0035 {
0036   const double beta = betagamma/sqrt(1.+betagamma*betagamma);
0037
0038   return A/(beta*beta)*TMath::Log(betagamma) + A/(beta*beta)*B - A - C;
0039 }
0040
0041 const double bethe_bloch_new_2D(const double betagamma, const double A, const double B)
0042 {
0043   const double beta = betagamma/sqrt(1.+betagamma*betagamma);
0044
0045   return A/(beta*beta)*(2.*TMath::Log(2.*dedx_constants::m_e/dedx_constants::sphenix_I * betagamma) - beta*beta) - B;
0046 }
0047
0048 const double bethe_bloch_new_1D(const double betagamma, const double A)
0049 {
0050   const double beta = betagamma/sqrt(1.+betagamma*betagamma);
0051
0052   return A/(beta*beta)*(2.*TMath::Log(2.*dedx_constants::m_e/dedx_constants::sphenix_I * betagamma) - beta*beta);
0053 }
0054
0055 // dE/dx for one gas species, up to normalization
0056 const double bethe_bloch_species(const double betagamma, const double I)
0057 {
0058   const double m_e = 511e3; // eV
0059
0060   const double beta = betagamma/sqrt(1.+betagamma*betagamma);
0061
0062   return 1./(beta*beta)*(TMath::Log(2.*m_e/I*betagamma*betagamma)-beta*beta);
0063 }
0064
0065 // dE/dx for TPC gas mixture, up to normalization
0066 const double bethe_bloch_total(const double betagamma)
0067 {
0068   return dedx_constants::ar_frac * bethe_bloch_species(betagamma,dedx_constants::ar_I) +
0069          dedx_constants::cf4_frac * bethe_bloch_species(betagamma,dedx_constants::cf4_I) +
0070          dedx_constants::isobutane_frac * bethe_bloch_species(betagamma,dedx_constants::isobutane_I);
0071 }
0072
0073 Double_t bethe_bloch_new_wrapper(Double_t* x, Double_t* par)
0074 {
0075   Double_t betagamma = x[0];
0076   Double_t A = par[0];
0077   Double_t B = par[1];
0078   Double_t C = par[2];
0079
0080   return bethe_bloch_new(betagamma,A,B,C);
0081 }
0082
0083 Double_t bethe_bloch_new_2D_wrapper(Double_t* x, Double_t* par)
0084 {
0085   Double_t betagamma = x[0];
0086   Double_t A = par[0];
0087   Double_t B = par[1];
0088
0089   return bethe_bloch_new_2D(betagamma,A,B);
0090 }
0091
0092 Double_t bethe_bloch_new_1D_wrapper(Double_t* x, Double_t* par)
0093 {
0094   Double_t betagamma = x[0];
0095   Double_t A = par[0];
0096
0097   return bethe_bloch_new_1D(betagamma,A);
0098 }
0099
0100 // wrapper function for TF1 constructor, for fitting
0101 Double_t bethe_bloch_wrapper(Double_t* ln_bg, Double_t* par)
0102 {
0103   Double_t betagamma = exp(ln_bg[0]);
0104
0105   Double_t norm = par[0];
0106
0107   return norm * bethe_bloch_total(betagamma);
0108 }
0109
0110 Double_t bethe_bloch_vs_p_wrapper(Double_t* x, Double_t* par)
0111 {
0112   Double_t p = x[0];
0113   Double_t norm = par[0];
0114   Double_t m = par[1];
0115
0116   return norm * bethe_bloch_total(fabs(p)/m);
0117 }
0118
0119 Double_t bethe_bloch_vs_logp_wrapper(Double_t* x, Double_t* par)
0120 {
0121   Double_t p = pow(10.,x[0]);
0122   Double_t norm = par[0];
0123   Double_t m = par[1];
0124
0125   return norm * bethe_bloch_total(fabs(p)/m);
0126 }
0127
0128 Double_t bethe_bloch_vs_p_wrapper_ZS(Double_t* x, Double_t* par)
0129 {
0130   Double_t p = x[0];
0131   Double_t norm = par[0];
0132   Double_t m = par[1];
0133   Double_t ZS_loss = par[2];
0134
0135   return norm * bethe_bloch_total(fabs(p)/m) - ZS_loss;
0136 }
0137
0138 Double_t bethe_bloch_vs_p_wrapper_new(Double_t* x, Double_t* par)
0139 {
0140   Double_t p = x[0];
0141   Double_t A = par[0];
0142   Double_t B = par[1];
0143   Double_t C = par[2];
0144   Double_t m = par[3];
0145
0146   return bethe_bloch_new(fabs(p)/m,A,B,C);
0147 }
0148
0149 Double_t bethe_bloch_vs_p_wrapper_new_2D(Double_t* x, Double_t* par)
0150 {
0151   Double_t p = x[0];
0152   Double_t A = par[0];
0153   Double_t B = par[1];
0154   Double_t m = par[2];
0155
0156   return bethe_bloch_new_2D(fabs(p)/m,A,B);
0157 }
0158
0159 Double_t bethe_bloch_vs_p_wrapper_new_1D(Double_t* x, Double_t* par)
0160 {
0161   Double_t p = x[0];
0162   Double_t A = par[0];
0163   Double_t m = par[1];
0164
0165   return bethe_bloch_new_1D(fabs(p)/m,A);
0166 }
0167
0168 Double_t bethe_bloch_vs_logp_wrapper_ZS(Double_t* x, Double_t* par)
0169 {
0170   Double_t p = pow(10.,x[0]);
0171   Double_t norm = par[0];
0172   Double_t m = par[1];
0173   Double_t ZS_loss = par[2];
0174
0175   return norm * bethe_bloch_total(fabs(p)/m) - ZS_loss;
0176 }
0177
0178 Double_t bethe_bloch_vs_logp_wrapper_new(Double_t* x, Double_t* par)
0179 {
0180   Double_t p = pow(10.,x[0]);
0181   Double_t A = par[0];
0182   Double_t B = par[1];
0183   Double_t C = par[2];
0184   Double_t m = par[3];
0185
0186   return bethe_bloch_new(fabs(p)/m,A,B,C);
0187 }
0188
0189 Double_t bethe_bloch_vs_logp_wrapper_new_1D(Double_t* x, Double_t* par)
0190 {
0191   Double_t p = pow(10.,x[0]);
0192   Double_t A = par[0];
0193   Double_t m = par[1];
0194
0195   return bethe_bloch_new_1D(fabs(p)/m,A);
0196 }
0197
0198 // ratio of dE/dx between two particle species at the same momentum
0199 // (useful for dE/dx peak fits)
0200 const double dedx_ratio(const double p, const double m1, const double m2)
0201 {
0202   const double betagamma1 = fabs(p)/m1;
0203   const double betagamma2 = fabs(p)/m2;
0204
0205   return bethe_bloch_total(betagamma1)/bethe_bloch_total(betagamma2);
0206 }
0207
0208 #endif // BETHE_BLOCH_H