All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public 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 ()
 
- Public Member Functions inherited from genf::GFAbsEnergyLoss
virtual ~GFAbsEnergyLoss ()
 

Additional Inherited Members

- Protected Member Functions inherited from genf::GFAbsEnergyLoss
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 46 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.

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 }
I
see a below echo or echo I(indirect symbol).'echo" If the symbol is local (non-external)
S
see a below echo S(symbol in a section other than those above)
pdg
var pdg
Definition: selectors.fcl:14
E
process_name E
Definition: rawtorootconvert.fcl:26
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
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:48
lib.dedx.mass
float mass
Definition: dedx.py: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:78
The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5