LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Protected Member Functions | List of all members
genf::GFEnergyLossBrems Class Reference

#include "GFEnergyLossBrems.h"

Inheritance diagram for genf::GFEnergyLossBrems:
genf::GFAbsEnergyLoss

Public Member Functions

double energyLoss (const double &step, const double &mom, const int &pdg, const double &matDensity, const double &matZ, const double &matA, const double &radiationLength, const double &meanExcitationEnergy, const bool &doNoise=false, TMatrixT< Double_t > *noise=NULL, const TMatrixT< Double_t > *jacobian=NULL, const TVector3 *directionBefore=NULL, const TVector3 *directionAfter=NULL)
 Returns energy loss, optional calculation of energy loss straggeling. More...
 
virtual ~GFEnergyLossBrems ()
 

Protected Member Functions

void getParticleParameters (const int &pdg, double &charge, double &mass)
 Gets particle charge and mass (in GeV) More...
 
double getParticleMass (const int &pdg)
 Returns particle mass (in GeV) More...
 

Detailed Description

Definition at line 45 of file GFEnergyLossBrems.h.

Constructor & Destructor Documentation

genf::GFEnergyLossBrems::~GFEnergyLossBrems ( )
virtual

Definition at line 26 of file GFEnergyLossBrems.cxx.

26 {}

Member Function Documentation

double genf::GFEnergyLossBrems::energyLoss ( const double &  step,
const double &  mom,
const int &  pdg,
const double &  matDensity,
const double &  matZ,
const double &  matA,
const double &  radiationLength,
const double &  meanExcitationEnergy,
const bool &  doNoise = false,
TMatrixT< Double_t > *  noise = NULL,
const TMatrixT< Double_t > *  jacobian = NULL,
const TVector3 *  directionBefore = NULL,
const TVector3 *  directionAfter = NULL 
)
virtual

Returns energy loss, optional calculation of energy loss straggeling.

Can be called with any pdg, but only calculates energy loss and straggeling for electrons and positrons (otherwise returns 0). Uses a gaussian approximation (Bethe-Heitler formula with Migdal corrections). For positrons the energy loss is weighed with a correction factor.

Implements genf::GFAbsEnergyLoss.

Definition at line 28 of file GFEnergyLossBrems.cxx.

References E, genf::GFAbsEnergyLoss::getParticleMass(), genf::GFAbsEnergyLoss::getParticleParameters(), util::kc, GFException::setFatal(), X, and Y.

41 {
42 
43  if (fabs(pdg == 11)) { // only for electrons and positrons
44 #if !defined(BETHE)
45  static const double C[101] = {
46  0.0, -0.960613E-01, 0.631029E-01, -0.142819E-01, 0.150437E-02, -0.733286E-04,
47  0.131404E-05, 0.859343E-01, -0.529023E-01, 0.131899E-01, -0.159201E-02, 0.926958E-04,
48  -0.208439E-05, -0.684096E+01, 0.370364E+01, -0.786752E+00, 0.822670E-01, -0.424710E-02,
49  0.867980E-04, -0.200856E+01, 0.129573E+01, -0.306533E+00, 0.343682E-01, -0.185931E-02,
50  0.392432E-04, 0.127538E+01, -0.515705E+00, 0.820644E-01, -0.641997E-02, 0.245913E-03,
51  -0.365789E-05, 0.115792E+00, -0.463143E-01, 0.725442E-02, -0.556266E-03, 0.208049E-04,
52  -0.300895E-06, -0.271082E-01, 0.173949E-01, -0.452531E-02, 0.569405E-03, -0.344856E-04,
53  0.803964E-06, 0.419855E-02, -0.277188E-02, 0.737658E-03, -0.939463E-04, 0.569748E-05,
54  -0.131737E-06, -0.318752E-03, 0.215144E-03, -0.579787E-04, 0.737972E-05, -0.441485E-06,
55  0.994726E-08, 0.938233E-05, -0.651642E-05, 0.177303E-05, -0.224680E-06, 0.132080E-07,
56  -0.288593E-09, -0.245667E-03, 0.833406E-04, -0.129217E-04, 0.915099E-06, -0.247179E-07,
57  0.147696E-03, -0.498793E-04, 0.402375E-05, 0.989281E-07, -0.133378E-07, -0.737702E-02,
58  0.333057E-02, -0.553141E-03, 0.402464E-04, -0.107977E-05, -0.641533E-02, 0.290113E-02,
59  -0.477641E-03, 0.342008E-04, -0.900582E-06, 0.574303E-05, 0.908521E-04, -0.256900E-04,
60  0.239921E-05, -0.741271E-07, -0.341260E-04, 0.971711E-05, -0.172031E-06, -0.119455E-06,
61  0.704166E-08, 0.341740E-05, -0.775867E-06, -0.653231E-07, 0.225605E-07, -0.114860E-08,
62  -0.119391E-06, 0.194885E-07, 0.588959E-08, -0.127589E-08, 0.608247E-10};
63  static const double xi = 2.51, beta = 0.99, vl = 0.00004;
64 #endif
65 #if defined(BETHE) // no MIGDAL corrections
66  static const double C[101] = {
67  0.0, 0.834459E-02, 0.443979E-02, -0.101420E-02, 0.963240E-04, -0.409769E-05,
68  0.642589E-07, 0.464473E-02, -0.290378E-02, 0.547457E-03, -0.426949E-04, 0.137760E-05,
69  -0.131050E-07, -0.547866E-02, 0.156218E-02, -0.167352E-03, 0.101026E-04, -0.427518E-06,
70  0.949555E-08, -0.406862E-02, 0.208317E-02, -0.374766E-03, 0.317610E-04, -0.130533E-05,
71  0.211051E-07, 0.158941E-02, -0.385362E-03, 0.315564E-04, -0.734968E-06, -0.230387E-07,
72  0.971174E-09, 0.467219E-03, -0.154047E-03, 0.202400E-04, -0.132438E-05, 0.431474E-07,
73  -0.559750E-09, -0.220958E-02, 0.100698E-02, -0.596464E-04, -0.124653E-04, 0.142999E-05,
74  -0.394378E-07, 0.477447E-03, -0.184952E-03, -0.152614E-04, 0.848418E-05, -0.736136E-06,
75  0.190192E-07, -0.552930E-04, 0.209858E-04, 0.290001E-05, -0.133254E-05, 0.116971E-06,
76  -0.309716E-08, 0.212117E-05, -0.103884E-05, -0.110912E-06, 0.655143E-07, -0.613013E-08,
77  0.169207E-09, 0.301125E-04, -0.461920E-04, 0.871485E-05, -0.622331E-06, 0.151800E-07,
78  -0.478023E-04, 0.247530E-04, -0.381763E-05, 0.232819E-06, -0.494487E-08, -0.336230E-04,
79  0.223822E-04, -0.384583E-05, 0.252867E-06, -0.572599E-08, 0.105335E-04, -0.567074E-06,
80  -0.216564E-06, 0.237268E-07, -0.658131E-09, 0.282025E-05, -0.671965E-06, 0.565858E-07,
81  -0.193843E-08, 0.211839E-10, 0.157544E-04, -0.304104E-05, -0.624410E-06, 0.120124E-06,
82  -0.457445E-08, -0.188222E-05, -0.407118E-06, 0.375106E-06, -0.466881E-07, 0.158312E-08,
83  0.945037E-07, 0.564718E-07, -0.319231E-07, 0.371926E-08, -0.123111E-09};
84  static const double xi = 2.10, beta = 1.00, vl = 0.001;
85 #endif
86 
87  double BCUT = 10000.; // energy up to which soft bremsstrahlung energy loss is calculated
88 
89  static const double me = getParticleMass(11); // electron mass (GeV) 0.000519...
90 
91  double charge, mass;
92  getParticleParameters(pdg, charge, mass);
93 
94  double THIGH = 100., CHIGH = 50.;
95 
96  double dedxBrems = 0.;
97 
98  if (BCUT > 0.) {
99 
100  double T, kc;
101 
102  if (BCUT >= mom) BCUT = mom; // confine BCUT to mom
103 
104  // T=mom, confined to THIGH
105  // kc=BCUT, confined to CHIGH ??
106  if (mom >= THIGH) {
107  T = THIGH;
108  if (BCUT >= THIGH)
109  kc = CHIGH;
110  else
111  kc = BCUT;
112  }
113  else {
114  T = mom;
115  kc = BCUT;
116  }
117 
118  double E = T + me; // total electron energy
119  if (BCUT > T) kc = T;
120 
121  double X = log(T / me);
122  double Y = log(kc / (E * vl));
123 
124  double XX;
125  int K;
126  double S = 0., YY = 1.;
127 
128  for (unsigned int I = 1; I <= 2; ++I) {
129  XX = 1.;
130  for (unsigned int J = 1; J <= 6; ++J) {
131  K = 6 * I + J - 6;
132  S = S + C[K] * XX * YY;
133  XX = XX * X;
134  }
135  YY = YY * Y;
136  }
137 
138  for (unsigned int I = 3; I <= 6; ++I) {
139  XX = 1.;
140  for (unsigned int J = 1; J <= 6; ++J) {
141  K = 6 * I + J - 6;
142  if (Y <= 0.)
143  S = S + C[K] * XX * YY;
144  else
145  S = S + C[K + 24] * XX * YY;
146  XX = XX * X;
147  }
148  YY = YY * Y;
149  }
150 
151  double SS = 0.;
152  YY = 1.;
153 
154  for (unsigned int I = 1; I <= 2; ++I) {
155  XX = 1.;
156  for (unsigned int J = 1; J <= 5; ++J) {
157  K = 5 * I + J + 55;
158  SS = SS + C[K] * XX * YY;
159  XX = XX * X;
160  }
161  YY = YY * Y;
162  }
163 
164  for (unsigned int I = 3; I <= 5; ++I) {
165  XX = 1.;
166  for (unsigned int J = 1; J <= 5; ++J) {
167  K = 5 * I + J + 55;
168  if (Y <= 0.)
169  SS = SS + C[K] * XX * YY;
170  else
171  SS = SS + C[K + 15] * XX * YY;
172  XX = XX * X;
173  }
174  YY = YY * Y;
175  }
176 
177  S = S + matZ * SS;
178 
179  if (S > 0.) {
180  double CORR = 1.;
181 #if !defined(CERNLIB_BETHE)
182  CORR = 1. / (1. + 0.805485E-10 * matDensity * matZ * E * E /
183  (matA * kc * kc)); // MIGDAL correction factor
184 #endif
185 
186  double FAC = matZ * (matZ + xi) * E * E * pow((kc * CORR / T), beta) / (E + me);
187  if (FAC <= 0.) return 0.;
188  dedxBrems = FAC * S;
189 
190  double RAT;
191 
192  if (mom > THIGH) {
193  if (BCUT < THIGH) {
194  RAT = BCUT / mom;
195  S = (1. - 0.5 * RAT + 2. * RAT * RAT / 9.);
196  RAT = BCUT / T;
197  S = S / (1. - 0.5 * RAT + 2. * RAT * RAT / 9.);
198  }
199  else {
200  RAT = BCUT / mom;
201  S = BCUT * (1. - 0.5 * RAT + 2. * RAT * RAT / 9.);
202  RAT = kc / T;
203  S = S / (kc * (1. - 0.5 * RAT + 2. * RAT * RAT / 9.));
204  }
205  dedxBrems = dedxBrems * S; // GeV barn
206  }
207 
208  dedxBrems = 0.60221367 * matDensity * dedxBrems / matA; // energy loss dE/dx [GeV/cm]
209  }
210  }
211 
212  if (dedxBrems < 0.)
213  throw GFException("genf::GFEnergyLossBrems::energyLoss(): negative dE/dx", __LINE__, __FILE__)
214  .setFatal();
215 
216  double factor = 1.; // positron correction factor
217 
218  if (pdg == -11) {
219  static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054;
220 
221  double ETA = 0.;
222  if (matZ > 0.) {
223  double X = log(AA * mom / matZ * matZ);
224  if (X > -8.) {
225  if (X >= +9.)
226  ETA = 1.;
227  else {
228  double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.);
229  ETA = 0.5 + atan(W) / M_PI;
230  }
231  }
232  }
233 
234  double E0;
235 
236  if (ETA < 0.0001)
237  factor = 1.E-10;
238  else if (ETA > 0.9999)
239  factor = 1.;
240  else {
241  E0 = BCUT / mom;
242  if (E0 > 1.) E0 = 1.;
243  if (E0 < 1.E-8)
244  factor = 1.;
245  else
246  factor = ETA * (1. - pow(1. - E0, 1. / ETA)) / E0;
247  }
248  }
249 
250  double DE = step * factor * dedxBrems; //always positive
251  double momLoss =
252  sqrt(mom * mom + 2. * sqrt(mom * mom + mass * mass) * DE + DE * DE) - mom; //always positive
253 
254  if (doNoise) {
255  if (!noise)
256  throw GFException(
257  "genf::GFEnergyLossBrems::energyLoss(): no noise matrix!", __LINE__, __FILE__)
258  .setFatal();
259 
260  double LX = 1.442695 * step / radiationLength;
261  double S2B = mom * mom * (1. / pow(3., LX) - 1. / pow(4., LX));
262  double DEDXB = pow(fabs(S2B), 0.5);
263  DEDXB = 1.2E9 * DEDXB; //eV
264  double sigma2E = DEDXB * DEDXB; //eV^2
265  sigma2E *= 1.E-18; // eV -> GeV
266 
267  (*noise)[6][6] += (mom * mom + mass * mass) / pow(mom, 6.) * sigma2E;
268  }
269 
270  return momLoss;
271  }
272  else
273  return 0.; // if not an electron/positron
274 }
Y
Float_t Y
Definition: plot.C:37
noise
unsigned int noise()
Definition: chem4.cc:261
genf::GFAbsEnergyLoss::getParticleMass
double getParticleMass(const int &pdg)
Returns particle mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:33
util::kc
constexpr double kc
Speed of light in vacuum in LArSoft units [cm/ns].
Definition: PhysicalConstants.h:59
E
Float_t E
Definition: plot.C:20
genf::GFAbsEnergyLoss::getParticleParameters
void getParticleParameters(const int &pdg, double &charge, double &mass)
Gets particle charge and mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:26
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
GFException::setFatal
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
X
Float_t X
Definition: plot.C:37
double genf::GFAbsEnergyLoss::getParticleMass ( const int &  pdg)
protectedinherited

Returns particle mass (in GeV)

Definition at line 33 of file GFAbsEnergyLoss.cxx.

References part.

Referenced by energyLoss(), and genf::GFEnergyLossBetheBloch::energyLoss().

34 {
35  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg);
36  return part->Mass();
37 }
part
TString part[npart]
Definition: Style.C:32
void genf::GFAbsEnergyLoss::getParticleParameters ( const int &  pdg,
double &  charge,
double &  mass 
)
protectedinherited

Gets particle charge and mass (in GeV)

Definition at line 26 of file GFAbsEnergyLoss.cxx.

References part.

Referenced by energyLoss(), genf::GFEnergyLossBetheBloch::energyLoss(), and genf::GFEnergyLossCoulomb::energyLoss().

27 {
28  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg);
29  charge = part->Charge() / (3.);
30  mass = part->Mass();
31 }
part
TString part[npart]
Definition: Style.C:32
The documentation for this class was generated from the following files: