LArSoft  v07_13_02
Liquid Argon Software toolkit - http://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  {
27 }

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 29 of file GFEnergyLossBrems.cxx.

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

Returns particle mass (in GeV)

Definition at line 35 of file GFAbsEnergyLoss.cxx.

References part.

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

35  {
36  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(pdg);
37  return part->Mass();
38 }
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 27 of file GFAbsEnergyLoss.cxx.

References part.

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

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