LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
GFEnergyLossBrems.cxx
Go to the documentation of this file.
1 /* Copyright 2008-2009, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "larreco/Genfit/GFEnergyLossBrems.h"
21 #include "larreco/Genfit/GFException.h"
22 #include "math.h"
23 
24 //#define BETHE // don't use MIGDAL correction for electron bremsstrahlung (only Bethe Heitler)
25 
26 genf::GFEnergyLossBrems::~GFEnergyLossBrems() {}
27 
28 double genf::GFEnergyLossBrems::energyLoss(const double& step,
29  const double& mom,
30  const int& pdg,
31  const double& matDensity,
32  const double& matZ,
33  const double& matA,
34  const double& radiationLength,
35  const double& /* meanExcitationEnergy */,
36  const bool& doNoise,
37  TMatrixT<Double_t>* noise,
38  const TMatrixT<Double_t>* /* jacobian */,
39  const TVector3* /* directionBefore */,
40  const TVector3* /* directionAfter */)
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 }
275 
276 //ClassImp(GFEnergyLossBrems)
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
genf::GFEnergyLossBrems::energyLoss
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.
Definition: GFEnergyLossBrems.cxx:28