LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
genf::GFMaterialEffects Class Reference

#include "GFMaterialEffects.h"

Inheritance diagram for genf::GFMaterialEffects:

Public Member Functions

void setEnergyLossBetheBloch (bool opt=true)
 
void setNoiseBetheBloch (bool opt=true)
 
void setNoiseCoulomb (bool opt=true)
 
void setEnergyLossBrems (bool opt=true)
 
void setNoiseBrems (bool opt=true)
 
double effects (const std::vector< TVector3 > &points, const std::vector< double > &pointPaths, const double &mom, const int &pdg, const bool &doNoise=false, TMatrixT< Double_t > *noise=NULL, const TMatrixT< Double_t > *jacobian=NULL, const TVector3 *directionBefore=NULL, const TVector3 *directionAfter=NULL)
 Calculates energy loss in the travelled path, optional calculation of noise matrix. More...
 
double stepper (const double &maxDist, const double &posx, const double &posy, const double &posz, const double &dirx, const double &diry, const double &dirz, const double &mom, const int &pdg)
 Returns maximum length so that a specified momentum loss will not be exceeded. More...
 
double stepper (const double &maxDist, const TVector3 &pos, const TVector3 &dir, const double &mom, const int &pdg)
 

Static Public Member Functions

static GFMaterialEffects * getInstance ()
 
static void destruct ()
 

Private Member Functions

 GFMaterialEffects ()
 
virtual ~GFMaterialEffects ()
 
void getParameters ()
 contains energy loss classes More...
 
void calcBeta (double mom)
 sets fbeta, fgamma, fgammasquare; must only be used after calling getParameters() More...
 
double energyLossBetheBloch (const double &mom)
 Returns energy loss. More...
 
void noiseBetheBloch (const double &mom, TMatrixT< double > *noise) const
 calculation of energy loss straggling More...
 
void noiseCoulomb (const double &mom, TMatrixT< double > *noise, const TMatrixT< double > *jacobian, const TVector3 *directionBefore, const TVector3 *directionAfter) const
 calculation of multiple scattering More...
 
double energyLossBrems (const double &mom) const
 Returns energy loss. More...
 
void noiseBrems (const double &mom, TMatrixT< double > *noise) const
 calculation of energy loss straggeling More...
 
double MeanExcEnergy_get (int Z)
 
double MeanExcEnergy_get (TGeoMaterial *)
 

Private Attributes

bool fEnergyLossBetheBloch
 
bool fNoiseBetheBloch
 
bool fNoiseCoulomb
 
bool fEnergyLossBrems
 
bool fNoiseBrems
 
const double me
 
double fstep
 
double fbeta
 
double fdedx
 
double fgamma
 
double fgammaSquare
 
double fmatDensity
 
double fmatZ
 
double fmatA
 
double fradiationLength
 
double fmEE
 
int fpdg
 
double fcharge
 
double fmass
 

Static Private Attributes

static GFMaterialEffects * finstance = NULL
 

Detailed Description

Definition at line 58 of file GFMaterialEffects.h.

Constructor & Destructor Documentation

genf::GFMaterialEffects::GFMaterialEffects ( )
private

Definition at line 55 of file GFMaterialEffects.cxx.

Referenced by getInstance().

55  :
56  fEnergyLossBetheBloch(true), fNoiseBetheBloch(true),
57  fNoiseCoulomb(true),
58  fEnergyLossBrems(true), fNoiseBrems(true),
59  me(0.510998910E-3),
60  fstep(0),
61  fbeta(0),
62  fdedx(0),
63  fgamma(0),
64  fgammaSquare(0),
65  fmatDensity(0),
66  fmatZ(0),
67  fmatA(0),
68  fradiationLength(0),
69  fmEE(0),
70  fpdg(0),
71  fcharge(0),
72  fmass(0) {
73 }
E
Float_t E
Definition: plot.C:23
genf::GFMaterialEffects::~GFMaterialEffects ( )
privatevirtual

Definition at line 39 of file GFMaterialEffects.cxx.

39  {
40  // for(unsigned int i=0;i<fEnergyLoss.size();++i) delete fEnergyLoss.at(i);
41  //delete geoMatManager;
42 }

Member Function Documentation

void genf::GFMaterialEffects::calcBeta ( double  mom)
private

sets fbeta, fgamma, fgammasquare; must only be used after calling getParameters()

Definition at line 275 of file GFMaterialEffects.cxx.

References E, fbeta, fgamma, fgammaSquare, and fmass.

Referenced by effects(), and stepper().

275  {
276  fbeta = mom/sqrt(fmass*fmass+mom*mom);
277 
278  //for numerical stability
279  fgammaSquare = 1.-fbeta*fbeta;
280  if(fgammaSquare>1.E-10) fgammaSquare = 1./fgammaSquare;
281  else fgammaSquare = 1.E10;
282  fgamma = sqrt(fgammaSquare);
283 }
E
Float_t E
Definition: plot.C:23
void genf::GFMaterialEffects::destruct ( )
static

Definition at line 80 of file GFMaterialEffects.cxx.

References finstance.

80  {
81  if(finstance != NULL) {
82  delete finstance;
83  finstance = NULL;
84  }
85 }
genf::GFMaterialEffects::finstance
static GFMaterialEffects * finstance
Definition: GFMaterialEffects.h:63
double genf::GFMaterialEffects::effects ( const std::vector< TVector3 > &  points,
const std::vector< double > &  pointPaths,
const double &  mom,
const int &  pdg,
const bool &  doNoise = false,
TMatrixT< Double_t > *  noise = NULL,
const TMatrixT< Double_t > *  jacobian = NULL,
const TVector3 *  directionBefore = NULL,
const TVector3 *  directionAfter = NULL 
)

Calculates energy loss in the travelled path, optional calculation of noise matrix.

Definition at line 87 of file GFMaterialEffects.cxx.

References calcBeta(), dir, E, energyLossBetheBloch(), energyLossBrems(), fEnergyLossBetheBloch, fEnergyLossBrems, fmatZ, fNoiseBetheBloch, fNoiseBrems, fNoiseCoulomb, fpdg, fstep, getParameters(), noiseBetheBloch(), noiseBrems(), noiseCoulomb(), and X.

Referenced by genf::RKTrackRep::Extrap(), and setNoiseBrems().

95  {
96 
97  //assert(points.size()==pointPaths.size());
98  fpdg = pdg;
99 
100  double momLoss=0.;
101 
102  for(unsigned int i=1;i<points.size();++i){
103  // TVector3 p1=points.at(i-1);
104  //TVector3 p2=points.at(i);
105  //TVector3 dir=p2-p1;
106  TVector3 dir=points.at(i)-points.at(i-1);
107  double dist=dir.Mag();
108  double realPath = pointPaths.at(i);
109 
110  if (dist > 1.E-8) { // do material effects only if distance is not too small
111  dir*=1./dist; //normalize dir
112 
113  double X(0.);
114  /*
115  double matDensity, matZ, matA, radiationLength, mEE;
116  double step;
117  */
118 
119  gGeoManager->InitTrack(points.at(i-1).X(),points.at(i-1).Y(),points.at(i-1).Z(),
120  dir.X(),dir.Y(),dir.Z());
121 
122  while(X<dist){
123 
124  getParameters();
125  // geoMatManager->getMaterialParameters(matDensity, matZ, matA, radiationLength, mEE);
126 
127  // step = geoMatManager->stepOrNextBoundary(dist-X);
128  gGeoManager->FindNextBoundaryAndStep(dist-X);
129  fstep = gGeoManager->GetStep();
130 
131  // Loop over EnergyLoss classes
132  if(fmatZ>1.E-3){
133  /*
134  for(unsigned int j=0;j<fEnergyLoss.size();++j){
135  momLoss += realPath/dist*fEnergyLoss.at(j)->energyLoss(step,
136  mom,
137  pdg,
138  matDensity,
139  matZ,
140  matA,
141  radiationLength,
142  mEE,
143  doNoise,
144  noise,
145  jacobian,
146  directionBefore,
147  directionAfter);
148  }
149  */
150  calcBeta(mom);
151 
152  if (fEnergyLossBetheBloch)
153  momLoss += realPath/dist * this->energyLossBetheBloch(mom);
154  if (doNoise && fEnergyLossBetheBloch && fNoiseBetheBloch)
155  this->noiseBetheBloch(mom, noise);
156 
157  if (/*doNoise &&*/ fNoiseCoulomb) // Force it
158  this->noiseCoulomb(mom, noise, jacobian, directionBefore, directionAfter);
159 
160  if (fEnergyLossBrems)
161  momLoss += realPath/dist * this->energyLossBrems(mom);
162  if (doNoise && fEnergyLossBrems && fNoiseBrems)
163  this->noiseBrems(mom, noise);
164 
165 
166 
167  }
168 
169  X += fstep;
170 
171  }
172  }
173  }
174  //std::cout << "GFMaterialEffects::effects(): momLoss " << momLoss << std::endl;
175  return momLoss;
176 }
genf::GFMaterialEffects::noiseBrems
void noiseBrems(const double &mom, TMatrixT< double > *noise) const
calculation of energy loss straggeling
Definition: GFMaterialEffects.cxx:645
E
Float_t E
Definition: plot.C:23
genf::GFMaterialEffects::noiseCoulomb
void noiseCoulomb(const double &mom, TMatrixT< double > *noise, const TMatrixT< double > *jacobian, const TVector3 *directionBefore, const TVector3 *directionAfter) const
calculation of multiple scattering
Definition: GFMaterialEffects.cxx:389
genf::GFMaterialEffects::energyLossBrems
double energyLossBrems(const double &mom) const
Returns energy loss.
Definition: GFMaterialEffects.cxx:483
genf::GFMaterialEffects::noiseBetheBloch
void noiseBetheBloch(const double &mom, TMatrixT< double > *noise) const
calculation of energy loss straggling
Definition: GFMaterialEffects.cxx:326
noise
unsigned int noise()
Definition: chem4.cc:265
genf::GFMaterialEffects::calcBeta
void calcBeta(double mom)
sets fbeta, fgamma, fgammasquare; must only be used after calling getParameters() ...
Definition: GFMaterialEffects.cxx:275
dir
TDirectory * dir
Definition: macro.C:5
genf::GFMaterialEffects::getParameters
void getParameters()
contains energy loss classes
Definition: GFMaterialEffects.cxx:254
X
Float_t X
Definition: plot.C:39
genf::GFMaterialEffects::energyLossBetheBloch
double energyLossBetheBloch(const double &mom)
Returns energy loss.
Definition: GFMaterialEffects.cxx:289
double genf::GFMaterialEffects::energyLossBetheBloch ( const double &  mom)
private

Returns energy loss.

Uses Bethe Bloch formula to calculate energy loss. Calcuates and sets fdedx which needed also for noiseBetheBloch. Therefore it is not a const function!

Definition at line 289 of file GFMaterialEffects.cxx.

References E, fbeta, fcharge, fdedx, fgammaSquare, fmass, fmatA, fmatDensity, fmatZ, fmEE, fstep, and me.

Referenced by effects(), and stepper().

289  {
290 
291  // calc fdedx, also needed in noiseBetheBloch!
292  fdedx = 0.307075*fmatZ/fmatA*fmatDensity/(fbeta*fbeta)*fcharge*fcharge;
293  double massRatio = me/fmass;
294  // me=0.000511 here is in GeV. So fmEE must come in here in eV to get converted to MeV.
295  double argument = fgammaSquare*fbeta*fbeta*me*1.E3*2./((1.E-6*fmEE) * sqrt(1+2*sqrt(fgammaSquare)*massRatio + massRatio*massRatio));
296 
297  if (fmass==0.0) return(0.0);
298  if (argument <= exp(fbeta*fbeta))
299  {
300  fdedx = 0.;
301  // so-called Anderson-Ziegler domain ... Let's approximate it with a flat
302  // 100 MeV/cm, looking at the muon dE/dx curve in the PRD Review, or, ahem, wikipedia.
303  // http://pdg.lbl.gov/2011/reviews/rpp2011-rev-passage-particles-matter.pdf
304  // But there's a practical problem: must not reduce momentum to zero for tracking in G3,
305  // so allow only 50% reduction of kinetic energy.
306  // fdedx = 100.0*1.E-3/fstep; // in GeV/cm, hence 1.e-3
307  //if (fdedx > 0.5*0.5*fmass*fbeta*fbeta/fstep) fdedx = 0.5*0.5*fmass*fbeta*fbeta/fstep;
308  }
309  else{
310  fdedx *= (log(argument)-fbeta*fbeta); // Bethe-Bloch [MeV/cm]
311  fdedx *= 1.E-3; // in GeV/cm, hence 1.e-3
312  if (fdedx<0.) fdedx = 0;
313  }
314 
315  double DE = fstep * fdedx; //always positive
316  double momLoss = sqrt(mom*mom+2.*sqrt(mom*mom+fmass*fmass)*DE+DE*DE) - mom; //always positive
317 
318  //in vacuum it can numerically happen that momLoss becomes a small negative number. A cut-off at 0.01 eV for momentum loss seems reasonable
319  if(fabs(momLoss)<1.E-11)momLoss=1.E-11;
320 
321  // if(fabs(momLoss)>mom)momLoss=mom; // For going fwd and ending its life. EC.
322  return momLoss;
323 }
E
Float_t E
Definition: plot.C:23
double genf::GFMaterialEffects::energyLossBrems ( const double &  mom) const
private

Returns energy loss.

Can be called with any pdg, but only calculates energy loss 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.

Definition at line 483 of file GFMaterialEffects.cxx.

References beta, E, fbeta, fmass, fmatA, fmatDensity, fmatZ, fpdg, fstep, K, util::kc, me, X, and Y.

Referenced by effects(), and stepper().

483  {
484 
485  if (fabs(fpdg)!=11) return 0; // only for electrons and positrons
486 
487  #if !defined(BETHE)
488  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};
489  static const double xi=2.51, beta=0.99, vl=0.00004;
490  #endif
491  #if defined(BETHE) // no MIGDAL corrections
492  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};
493  static const double xi=2.10, fbeta=1.00, vl=0.001;
494  #endif
495 
496  double BCUT=10000.; // energy up to which soft bremsstrahlung energy loss is calculated
497 
498  double THIGH=100., CHIGH=50.;
499  double dedxBrems=0.;
500 
501  if(BCUT>0.){
502  double T, kc;
503 
504  if(BCUT>=mom) BCUT=mom; // confine BCUT to mom
505 
506  // T=mom, confined to THIGH
507  // kc=BCUT, confined to CHIGH ??
508  if(mom>=THIGH) {
509  T=THIGH;
510  if(BCUT>=THIGH) kc=CHIGH;
511  else kc=BCUT;
512  }
513  else {
514  T=mom;
515  kc=BCUT;
516  }
517 
518  double E=T+me; // total electron energy
519  if(BCUT>T) kc=T;
520 
521  double X=log(T/me);
522  double Y=log(kc/(E*vl));
523 
524  double XX;
525  int K;
526  double S=0., YY=1.;
527 
528  for (unsigned int I=1; I<=2; ++I) {
529  XX=1.;
530  for (unsigned int J=1; J<=6; ++J) {
531  K=6*I+J-6;
532  S=S+C[K]*XX*YY;
533  XX=XX*X;
534  }
535  YY=YY*Y;
536  }
537 
538  for (unsigned int I=3; I<=6; ++I) {
539  XX=1.;
540  for (unsigned int J=1; J<=6; ++J) {
541  K=6*I+J-6;
542  if(Y<=0.) S=S+C[K]*XX*YY;
543  else S=S+C[K+24]*XX*YY;
544  XX=XX*X;
545  }
546  YY=YY*Y;
547  }
548 
549  double SS=0.;
550  YY=1.;
551 
552  for (unsigned int I=1; I<=2; ++I) {
553  XX=1.;
554  for (unsigned int J=1; J<=5; ++J) {
555  K=5*I+J+55;
556  SS=SS+C[K]*XX*YY;
557  XX=XX*X;
558  }
559  YY=YY*Y;
560  }
561 
562  for (unsigned int I=3; I<=5; ++I) {
563  XX=1.;
564  for (unsigned int J=1; J<=5; ++J) {
565  K=5*I+J+55;
566  if(Y<=0.) SS=SS+C[K]*XX*YY;
567  else SS=SS+C[K+15]*XX*YY;
568  XX=XX*X;
569  }
570  YY=YY*Y;
571  }
572 
573  S=S+fmatZ*SS;
574 
575  if(S>0.){
576  double CORR=1.;
577  #if !defined(BETHE)
578  CORR=1./(1.+0.805485E-10*fmatDensity*fmatZ*E*E/(fmatA*kc*kc)); // MIGDAL correction factor
579  #endif
580 
581  double FAC=fmatZ*(fmatZ+xi)*E*E * pow((kc*CORR/T),beta) / (E+me);
582  if(FAC<=0.) return 0.;
583  dedxBrems=FAC*S;
584 
585  double RAT;
586 
587  if(mom>THIGH) {
588  if(BCUT<THIGH) {
589  RAT=BCUT/mom;
590  S=(1.-0.5*RAT+2.*RAT*RAT/9.);
591  RAT=BCUT/T;
592  S=S/(1.-0.5*RAT+2.*RAT*RAT/9.);
593  }
594  else {
595  RAT=BCUT/mom;
596  S=BCUT*(1.-0.5*RAT+2.*RAT*RAT/9.);
597  RAT=kc/T;
598  S=S/(kc*(1.-0.5*RAT+2.*RAT*RAT/9.));
599  }
600  dedxBrems=dedxBrems*S; // GeV barn
601  }
602 
603  dedxBrems = 0.60221367*fmatDensity*dedxBrems/fmatA; // energy loss dE/dx [GeV/cm]
604  }
605  }
606 
607  if (dedxBrems<0.) dedxBrems = 0;
608 
609  double factor=1.; // positron correction factor
610 
611  if (fpdg==-11){
612  static const double AA=7522100., A1=0.415, A3=0.0021, A5=0.00054;
613 
614  double ETA=0.;
615  if(fmatZ>0.) {
616  double X=log(AA*mom/fmatZ*fmatZ);
617  if(X>-8.) {
618  if(X>=+9.) ETA=1.;
619  else {
620  double W=A1*X+A3*pow(X,3.)+A5*pow(X,5.);
621  ETA=0.5+atan(W)/M_PI;
622  }
623  }
624  }
625 
626  double E0;
627 
628  if(ETA<0.0001) factor=1.E-10;
629  else if (ETA>0.9999) factor=1.;
630  else {
631  E0=BCUT/mom;
632  if(E0>1.) E0=1.;
633  if(E0<1.E-8) factor=1.;
634  else factor = ETA * ( 1.-pow(1.-E0, 1./ETA) ) / E0;
635  }
636  }
637 
638  double DE = fstep * factor*dedxBrems; //always positive
639  double momLoss = sqrt(mom*mom+2.*sqrt(mom*mom+fmass*fmass)*DE+DE*DE) - mom; //always positive
640 
641  return momLoss;
642 }
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
util::kc
constexpr double kc
Speed of light in vacuum in LArSoft units [cm/ns].
Definition: PhysicalConstants.h:59
X
Float_t X
Definition: plot.C:39
genf::GFMaterialEffects * genf::GFMaterialEffects::getInstance ( )
static

Definition at line 75 of file GFMaterialEffects.cxx.

References finstance, and GFMaterialEffects().

Referenced by genf::RKTrackRep::Extrap(), and genf::RKTrackRep::RKutta().

75  {
76  if(finstance == NULL) finstance = new GFMaterialEffects();
77  return finstance;
78 }
genf::GFMaterialEffects::finstance
static GFMaterialEffects * finstance
Definition: GFMaterialEffects.h:63
void genf::GFMaterialEffects::getParameters ( )
private

contains energy loss classes

interface to material and geometry

Definition at line 254 of file GFMaterialEffects.cxx.

References fcharge, fmass, fmatA, fmatDensity, fmatZ, fmEE, fpdg, fradiationLength, mat, MeanExcEnergy_get(), and part.

Referenced by effects(), and stepper().

254  {
255  if (!gGeoManager->GetCurrentVolume()->GetMedium())
256  throw GFException(std::string(__func__) + ": no medium", __LINE__, __FILE__).setFatal();
257  TGeoMaterial * mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
258  fmatDensity = mat->GetDensity();
259  fmatZ = mat->GetZ();
260  fmatA = mat->GetA();
261  fradiationLength = mat->GetRadLen();
262  fmEE = MeanExcEnergy_get(mat);
263 
264  // You know what? F*ck it. Just force this to be LAr.... is what I *could/will* say here ....
265  // See comment in energyLossBetheBloch() for why fmEE is in eV here.
266  fmatDensity = 1.40; fmatZ = 18.0; fmatA = 39.95; fradiationLength=13.947; fmEE=188.0;
267 
268 
269  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fpdg);
270  fcharge = part->Charge()/(3.);
271  fmass = part->Mass();
272 }
part
TString part[npart]
Definition: Style.C:32
mat
Float_t mat
Definition: plot.C:40
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
double genf::GFMaterialEffects::MeanExcEnergy_get ( int  Z)
private

Definition at line 669 of file GFMaterialEffects.cxx.

References MeanExcEnergy_NELEMENTS, MeanExcEnergy_vals, GFException::setFatal(), and Z.

Referenced by getParameters(), MeanExcEnergy_get(), and stepper().

669  {
670  if ((Z < 0) || (Z > MeanExcEnergy_NELEMENTS))
671  throw GFException(std::string(__func__) + ": unsupported Z", __LINE__, __FILE__).setFatal();
672  return MeanExcEnergy_vals[Z];
673 }
MeanExcEnergy_NELEMENTS
const int MeanExcEnergy_NELEMENTS
Definition: GFMaterialEffects.cxx:666
Z
Float_t Z
Definition: plot.C:39
MeanExcEnergy_vals
const float MeanExcEnergy_vals[MeanExcEnergy_NELEMENTS]
Definition: GFMaterialEffects.cxx:667
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
double genf::GFMaterialEffects::MeanExcEnergy_get ( TGeoMaterial *  mat)
private

Definition at line 675 of file GFMaterialEffects.cxx.

References MeanExcEnergy_get().

675  {
676  if(mat->IsMixture()){
677  double logMEE = 0.;
678  double denom = 0.;
679  TGeoMixture *mix = (TGeoMixture*)mat;
680  for(int i=0;i<mix->GetNelements();++i){
681  int index = int(floor((mix->GetZmixt())[i]));
682  logMEE += 1./(mix->GetAmixt())[i]*(mix->GetWmixt())[i]*(mix->GetZmixt())[i]*log(MeanExcEnergy_get(index));
683  denom += (mix->GetWmixt())[i]*(mix->GetZmixt())[i]*1./(mix->GetAmixt())[i];
684  }
685  logMEE/=denom;
686  return exp(logMEE);
687  }
688  else{ // not a mixture
689  int index = int(floor(mat->GetZ()));
690  return MeanExcEnergy_get(index);
691  }
692 }
mat
Float_t mat
Definition: plot.C:40
void genf::GFMaterialEffects::noiseBetheBloch ( const double &  mom,
TMatrixT< double > *  noise 
) const
private

calculation of energy loss straggling

For the energy loss straggeling, different formulas are used for different regions:

  • Vavilov-Gaussian regime
  • Urban/Landau approximation
  • truncated Landau distribution
  • Urban model

Needs fdedx, which is calculated in energyLossBetheBloch, so it has to be calles afterwards!

Definition at line 326 of file GFMaterialEffects.cxx.

References larg4::f1, larg4::f2, fbeta, fcharge, fdedx, fgamma, fgammaSquare, fmass, fmatA, fmatDensity, fmatZ, fstep, and me.

Referenced by effects(), and stepper().

327  {
328 
329 
330  // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E);
331  double sigma2E = 0.;
332  double zeta = 153.4E3 * fcharge*fcharge/(fbeta*fbeta) * fmatZ/fmatA * fmatDensity * fstep; // eV
333  double Emax = 2.E9*me*fbeta*fbeta*fgammaSquare / (1. + 2.*fgamma*me/fmass + (me/fmass)*(me/fmass) ); // eV
334  double kappa = zeta/Emax;
335 
336  if (kappa > 0.01) { // Vavilov-Gaussian regime
337  sigma2E += zeta*Emax*(1.-fbeta*fbeta/2.); // eV^2
338  }
339  else { // Urban/Landau approximation
340  double alpha = 0.996;
341  double sigmaalpha = 15.76;
342  // calculate number of collisions Nc
343  double I = 16. * pow(fmatZ, 0.9); // eV
344  double f2 = 0.;
345  if (fmatZ > 2.) f2 = 2./fmatZ;
346  double f1 = 1. - f2;
347  double e2 = 10.*fmatZ*fmatZ; // eV
348  double e1 = pow( (I/pow(e2,f2)), 1./f1); // eV
349 
350  double mbbgg2 = 2.E9*fmass*fbeta*fbeta*fgammaSquare; // eV
351  double Sigma1 = fdedx*1.0E9 * f1/e1 * (log(mbbgg2 / e1) - fbeta*fbeta) / (log(mbbgg2 / I) - fbeta*fbeta) * 0.6; // 1/cm
352  double Sigma2 = fdedx*1.0E9 * f2/e2 * (log(mbbgg2 / e2) - fbeta*fbeta) / (log(mbbgg2 / I) - fbeta*fbeta) * 0.6; // 1/cm
353  double Sigma3 = fdedx*1.0E9 * Emax / ( I*(Emax+I)*log((Emax+I)/I) ) * 0.4; // 1/cm
354 
355  double Nc = (Sigma1 + Sigma2 + Sigma3)*fstep;
356 
357  if (Nc > 50.) { // truncated Landau distribution
358  // calculate sigmaalpha (see GEANT3 manual W5013)
359  double RLAMED = -0.422784 - fbeta*fbeta - log(zeta/Emax);
360  double RLAMAX = 0.60715 + 1.1934*RLAMED +(0.67794 + 0.052382*RLAMED)*exp(0.94753+0.74442*RLAMED);
361  // from lambda max to sigmaalpha=sigma (empirical polynomial)
362  if(RLAMAX <= 1010.) {
363  sigmaalpha = 1.975560
364  +9.898841e-02 *RLAMAX
365  -2.828670e-04 *RLAMAX*RLAMAX
366  +5.345406e-07 *pow(RLAMAX,3.)
367  -4.942035e-10 *pow(RLAMAX,4.)
368  +1.729807e-13 *pow(RLAMAX,5.);
369  }
370  else { sigmaalpha = 1.871887E+01 + 1.296254E-02 *RLAMAX; }
371  // alpha=54.6 corresponds to a 0.9996 maximum cut
372  if(sigmaalpha > 54.6) sigmaalpha=54.6;
373  sigma2E += sigmaalpha*sigmaalpha * zeta*zeta; // eV^2
374  }
375  else { // Urban model
376  double Ealpha = I / (1.-(alpha*Emax/(Emax+I))); // eV
377  double meanE32 = I*(Emax+I)/Emax * (Ealpha-I); // eV^2
378  sigma2E += fstep * (Sigma1*e1*e1 + Sigma2*e2*e2 + Sigma3*meanE32); // eV^2
379  }
380  }
381 
382  sigma2E*=1.E-18; // eV -> GeV
383 
384  // update noise matrix
385  (*noise)[6][6] += (mom*mom+fmass*fmass)/pow(mom,6.)*sigma2E;
386 }
larg4::f2
SynchrotronAndGN f2
Definition: CustomPhysicsBuiltIns.cxx:11
larg4::f1
EmPhysicsFactory f1
Definition: CustomPhysicsBuiltIns.cxx:10
void genf::GFMaterialEffects::noiseBrems ( const double &  mom,
TMatrixT< double > *  noise 
) const
private

calculation of energy loss straggeling

Can be called with any pdg, but only calculates straggeling for electrons and positrons.

Definition at line 645 of file GFMaterialEffects.cxx.

References fmass, fpdg, fradiationLength, and fstep.

Referenced by effects(), and stepper().

646  {
647 
648  if (fabs(fpdg)!=11) return; // only for electrons and positrons
649 
650  double LX = 1.442695*fstep/fradiationLength;
651  double S2B = mom*mom * ( 1./pow(3.,LX) - 1./pow(4.,LX) );
652  double DEDXB = pow(fabs(S2B),0.5);
653  DEDXB = 1.2E9*DEDXB; //eV
654  double sigma2E = DEDXB*DEDXB; //eV^2
655  sigma2E*=1.E-18; // eV -> GeV
656 
657  (*noise)[6][6] += (mom*mom+fmass*fmass)/pow(mom,6.)*sigma2E;
658 }
void genf::GFMaterialEffects::noiseCoulomb ( const double &  mom,
TMatrixT< double > *  noise,
const TMatrixT< double > *  jacobian,
const TVector3 *  directionBefore,
const TVector3 *  directionAfter 
) const
private

calculation of multiple scattering

With the calculated multiple scattering angle, two noise matrices are calculated:

  • with respect to #directionBefore: #noiseBefore, which is then propagated with the jacobian
  • with respect to #directionAfter: #noiseAfter The mean value of these two matrices is added to the noise matrix noise. This method gives better results than either calculating only noiseBefore or noiseAfter.

    This is a detailed description of the mathematics involved:
    We define a local coordinate system cs' with initial momentum in z-direction:

    \[ \left(\begin{array}{c} x'\\ y'\\ z'\\ a_{x}'\\ a_{y}'\\ a_{z}'\\ \frac{q}{p}'\end{array}\right)=\left(\begin{array}{ccccccc} \cos\psi & \sin\psi & 0\\ -\cos\vartheta\sin\psi & \cos\vartheta\cos\psi & \sin\psi\\ \sin\vartheta\sin\psi & -\sin\vartheta\cos\psi & \cos\vartheta\\ & & & \cos\psi & \sin\psi & 0\\ & & & -\cos\vartheta\sin\psi & \cos\vartheta\cos\psi & \sin\psi\\ & & & \sin\vartheta\sin\psi & -\sin\vartheta\cos\psi & \cos\vartheta\\ & & & & & & 1\end{array}\right)\left(\begin{array}{c} x\\ y\\ z\\ a_{x}\\ a_{y}\\ a_{z}\\ \frac{q}{p}\end{array}\right)=R^{T}\left(\begin{array}{c} x\\ y\\ z\\ a_{x}\\ a_{y}\\ a_{z}\\ \frac{q}{p}\end{array}\right)\]


Now the global coordinate system cs is:

\[ \left(\begin{array}{c} x\\ y\\ z\\ a_{x}\\ a_{y}\\ a_{z}\\ \frac{q}{p}\end{array}\right)=\left(\begin{array}{ccccccc} \cos\psi & -\cos\vartheta\sin\psi & \sin\vartheta\sin\psi\\ \sin\psi & \cos\vartheta\cos\psi & -\sin\vartheta\cos\psi\\ 0 & \sin\psi & \cos\vartheta\\ & & & \cos\psi & -\cos\vartheta\sin\psi & \sin\vartheta\sin\psi\\ & & & \sin\psi & \cos\vartheta\cos\psi & -\sin\vartheta\cos\psi\\ & & & 0 & \sin\psi & \cos\vartheta\\ & & & & & & 1\end{array}\right)\left(\begin{array}{c} x'\\ y'\\ z'\\ a_{x}'\\ a_{y}'\\ a_{z}'\\ \frac{q}{p}'\end{array}\right)=R\left(\begin{array}{c} x'\\ y'\\ z'\\ a_{x}'\\ a_{y}'\\ a_{z}'\\ \frac{q}{p}'\end{array}\right) \]


with the Euler angles

\begin{eqnarray*} \psi & = & \begin{cases} \begin{cases} \frac{\pi}{2} & a_{x} \geq 0 \\ \frac{3\pi}{2} & a_{x} < 0 \end{cases} & a_{y}=0 \mbox{ resp. } |a_{y}|<10^{-14} \\ - \arctan \frac{a_{x}}{a_{y}} & a_{y} < 0 \\ \pi - \arctan \frac{a_{x}}{a_{y}} & a_{y} > 0 \end{cases} \\ \vartheta & = & \arccos a_{z} \end{eqnarray*}


$M$ is the multiple scattering error-matrix in the $\theta$ coordinate system. $\theta_{1/2}=0$ are the multiple scattering angles. There is only an error in direction (which later leads to an error in position, when $N_{before}$ is propagated with $T$).

\[ M=\left(\begin{array}{cc} \sigma^{2} & 0\\ 0 & \sigma^{2}\end{array}\right)\]


This translates into the noise matrix $\overline{M}$ in the local cs' via jacobian J. The mean scattering angle is always 0, this means that $\theta_{1/2}=0$):

\begin{eqnarray*} x' & = & x'\\ y' & = & y'\\ z' & = & z'\\ a_{x}' & = & \sin\theta_{1}\\ a_{y}' & = & \sin\theta_{2}\\ a_{z}' & = & \sqrt{1-\sin^{2}\theta_{1}-\sin^{2}\theta_{2}}\\ \frac{q}{p}' & = & \frac{q}{p}'\end{eqnarray*}

\[ M=\left(\begin{array}{cc} \sigma^{2} & 0\\ 0 & \sigma^{2}\end{array}\right)\]


\[ J=\frac{\partial\left(x',y',z',a_{x}',a_{y}',a_{z}',\frac{q}{p}'\right)}{\partial\left(\theta_{1},\theta_{2}\right)}\]


\[ J=\left(\begin{array}{cc} 0 & 0\\ 0 & 0\\ 0 & 0\\ \cos\theta_{1} & 0\\ 0 & \cos\theta_{2}\\ -\frac{\cos\theta_{1}\sin\theta_{1}}{\sqrt{\cos^{2}\theta_{1}-\sin^{2}\theta_{2}}} & -\frac{\cos\theta_{2}\sin\theta_{2}}{\sqrt{\cos^{2}\theta_{1}-\sin^{2}\theta_{2}}}\\ 0 & 0\end{array}\right) \overset{\theta_{1/2}=0}{=} \left(\begin{array}{cc} 0 & 0\\ 0 & 0\\ 0 & 0\\ 1 & 0\\ 0 & 1\\ 0 & 0\\ 0 & 0\end{array}\right)\]


\[ \overline{M}=J\: M\: J^{T}=\left(\begin{array}{ccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \sigma^{2} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \sigma^{2} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right)\]

The following transformation brings the noise matrix into the global coordinate system cs, resulting in $N$ :

\[ N=R\overline{M}R^{T}=\sigma^{2}\left(\begin{array}{ccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \cos^{2}\psi+\cos^{2}\theta-\cos^{2}\theta\cos^{2}\psi & \cos\psi\sin\psi\sin^{2}\theta & -\cos\theta\sin\psi\sin\theta & 0\\ 0 & 0 & 0 & \cos\psi\sin\psi\sin^{2}\theta & \sin^{2}\psi+\cos^{2}\theta\cos^{2}\psi & \cos\theta\cos\psi\sin\theta & 0\\ 0 & 0 & 0 & -\cos\theta\sin\psi\sin\theta & \cos\theta\cos\psi\sin\theta & \sin^{2}\theta & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right)\]


Now two $N$ are calculated: $N_{before}$ (at the start point, with initial direction #directionBefore) and $N_{after}$ (at the final point, with direction #directionAfter). $N_{before}$ is the propagated with the #jacobian of extrapolation $T$. The new covariance matrix with multiple scattering in global cs is:

\begin{eqnarray*} C_{new} & = C_{old} + 0.5 \cdot T N_{before} T^{T} + 0.5 \cdot N_{after} \end{eqnarray*}



See also: Track following in dense media and inhomogeneous magnetic fields, A.Fontana, P.Genova, L.Lavezzi and A.Rotondi; PANDA Report PV/01-07

Definition at line 389 of file GFMaterialEffects.cxx.

References E, fbeta, fmatZ, fradiationLength, and fstep.

Referenced by effects(), and stepper().

393  {
394 
395  // MULTIPLE SCATTERING; calculate sigma^2
396  // PANDA report PV/01-07 eq(43); linear in step length
397  double sigma2 = 225.E-6/(fbeta*fbeta*mom*mom) * fstep/fradiationLength * fmatZ/(fmatZ+1) * log(159.*pow(fmatZ,-1./3.))/log(287.*pow(fmatZ,-0.5)); // sigma^2 = 225E-6/mom^2 * XX0/fbeta^2 * Z/(Z+1) * ln(159*Z^(-1/3))/ln(287*Z^(-1/2)
398 
399  // noiseBefore
400  TMatrixT<double> noiseBefore(7,7);
401 
402  // calculate euler angles theta, psi (so that directionBefore' points in z' direction)
403  double psi = 0;
404  if (fabs((*directionBefore)[1]) < 1E-14) { // numerical case: arctan(+-inf)=+-PI/2
405  if ((*directionBefore)[0] >= 0.) psi = M_PI/2.;
406  else psi = 3.*M_PI/2.;
407  }
408  else {
409  if ((*directionBefore)[1] > 0.) psi = M_PI - atan((*directionBefore)[0]/(*directionBefore)[1]);
410  else psi = -atan((*directionBefore)[0]/(*directionBefore)[1]);
411  }
412  // cache sin and cos
413  double sintheta = sqrt(1-(*directionBefore)[2]*(*directionBefore)[2]); // theta = arccos(directionBefore[2])
414  double costheta = (*directionBefore)[2];
415  double sinpsi = sin(psi);
416  double cospsi = cos(psi);
417 
418  // calculate NoiseBefore Matrix R M R^T
419  double noiseBefore34 = sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseBefore_ij = noiseBefore_ji
420  double noiseBefore35 = -sigma2 * costheta * sinpsi * sintheta;
421  double noiseBefore45 = sigma2 * costheta * cospsi * sintheta;
422 
423  noiseBefore[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
424  noiseBefore[4][3] = noiseBefore34;
425  noiseBefore[5][3] = noiseBefore35;
426 
427  noiseBefore[3][4] = noiseBefore34;
428  noiseBefore[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
429  noiseBefore[5][4] = noiseBefore45;
430 
431  noiseBefore[3][5] = noiseBefore35;
432  noiseBefore[4][5] = noiseBefore45;
433  noiseBefore[5][5] = sigma2 * sintheta*sintheta;
434 
435  TMatrixT<double> jacobianT(7,7);
436  jacobianT = (*jacobian);
437  jacobianT.T();
438 
439  noiseBefore = jacobianT*noiseBefore*(*jacobian); //propagate
440 
441  // noiseAfter
442  TMatrixT<double> noiseAfter(7,7);
443 
444  // calculate euler angles theta, psi (so that A' points in z' direction)
445  psi = 0;
446  if (fabs((*directionAfter)[1]) < 1E-14) { // numerical case: arctan(+-inf)=+-PI/2
447  if ((*directionAfter)[0] >= 0.) psi = M_PI/2.;
448  else psi = 3.*M_PI/2.;
449  }
450  else {
451  if ((*directionAfter)[1] > 0.) psi = M_PI - atan((*directionAfter)[0]/(*directionAfter)[1]);
452  else psi = -atan((*directionAfter)[0]/(*directionAfter)[1]);
453  }
454  // cache sin and cos
455  sintheta = sqrt(1-(*directionAfter)[2]*(*directionAfter)[2]); // theta = arccos(directionAfter[2])
456  costheta = (*directionAfter)[2];
457  sinpsi = sin(psi);
458  cospsi = cos(psi);
459 
460  // calculate NoiseAfter Matrix R M R^T
461  double noiseAfter34 = sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseAfter_ij = noiseAfter_ji
462  double noiseAfter35 = -sigma2 * costheta * sinpsi * sintheta;
463  double noiseAfter45 = sigma2 * costheta * cospsi * sintheta;
464 
465  noiseAfter[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
466  noiseAfter[4][3] = noiseAfter34;
467  noiseAfter[5][3] = noiseAfter35;
468 
469  noiseAfter[3][4] = noiseAfter34;
470  noiseAfter[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
471  noiseAfter[5][4] = noiseAfter45;
472 
473  noiseAfter[3][5] = noiseAfter35;
474  noiseAfter[4][5] = noiseAfter45;
475  noiseAfter[5][5] = sigma2 * sintheta*sintheta;
476 
477  //calculate mean of noiseBefore and noiseAfter and update noise
478  (*noise) += 0.5*noiseBefore + 0.5*noiseAfter;
479 
480 }
E
Float_t E
Definition: plot.C:23
void genf::GFMaterialEffects::setEnergyLossBetheBloch ( bool  opt = true)
inline

Definition at line 69 of file GFMaterialEffects.h.

References fEnergyLossBetheBloch.

69 {fEnergyLossBetheBloch=opt;}
void genf::GFMaterialEffects::setEnergyLossBrems ( bool  opt = true)
inline

Definition at line 72 of file GFMaterialEffects.h.

References fEnergyLossBrems.

72 {fEnergyLossBrems=opt;}
void genf::GFMaterialEffects::setNoiseBetheBloch ( bool  opt = true)
inline

Definition at line 70 of file GFMaterialEffects.h.

References fNoiseBetheBloch.

70 {fNoiseBetheBloch=opt;}
void genf::GFMaterialEffects::setNoiseBrems ( bool  opt = true)
inline

Definition at line 73 of file GFMaterialEffects.h.

References effects(), fNoiseBrems, noise(), and stepper().

73 {fNoiseBrems=opt;}
void genf::GFMaterialEffects::setNoiseCoulomb ( bool  opt = true)
inline

Definition at line 71 of file GFMaterialEffects.h.

References fNoiseCoulomb.

71 {fNoiseCoulomb=opt;}
double genf::GFMaterialEffects::stepper ( const double &  maxDist,
const double &  posx,
const double &  posy,
const double &  posz,
const double &  dirx,
const double &  diry,
const double &  dirz,
const double &  mom,
const int &  pdg 
)

Returns maximum length so that a specified momentum loss will not be exceeded.

The stepper returns the maximum length that the particle may travel, so that a specified relative momentum loss will not be exceeded.

Definition at line 182 of file GFMaterialEffects.cxx.

References calcBeta(), E, energyLossBetheBloch(), energyLossBrems(), fEnergyLossBetheBloch, fEnergyLossBrems, fmatZ, fstep, getParameters(), and X.

Referenced by genf::RKTrackRep::RKutta(), setNoiseBrems(), and stepper().

190  {
191 
192  static const double maxPloss = .005; // maximum relative momentum loss allowed
193 
194  gGeoManager->InitTrack(posx,posy,posz,dirx,diry,dirz);
195 
196  double X(0.);
197  double dP = 0.;
198  double momLoss = 0.;
199  /*
200  double matDensity, matZ, matA, radiationLength, mEE;
201  double step;
202  */
203 
204  while(X<maxDist){
205 
206  getParameters();
207 
208  gGeoManager->FindNextBoundaryAndStep(maxDist-X);
209  fstep = gGeoManager->GetStep();
210  //
211  // step = geoMatManager->stepOrNextBoundary(maxDist-X);
212 
213  // Loop over EnergyLoss classes
214 
215  if(fmatZ>1.E-3){
216  /*
217  for(unsigned int j=0;j<fEnergyLoss.size();++j){
218  momLoss += fEnergyLoss.at(j)->energyLoss(step,
219  mom,
220  pdg,
221  matDensity,
222  matZ,
223  matA,
224  radiationLength,
225  mEE);
226  }
227  */
228  calcBeta(mom);
229 
230  if (fEnergyLossBetheBloch)
231  momLoss += this->energyLossBetheBloch(mom);
232 
233  if (fEnergyLossBrems)
234  momLoss += this->energyLossBrems(mom);
235 
236  }
237 
238  if(dP + momLoss > mom*maxPloss){
239  double fraction = (mom*maxPloss-dP)/momLoss;
240  if ((fraction <= 0.) || (fraction >= 1.))
241  throw GFException(std::string(__func__) + ": invalid fraction", __LINE__, __FILE__).setFatal();
242  dP+=fraction*momLoss;
243  X+=fraction*fstep;
244  break;
245  }
246 
247  dP += momLoss;
248  X += fstep;
249  }
250 
251  return X;
252 }
E
Float_t E
Definition: plot.C:23
genf::GFMaterialEffects::energyLossBrems
double energyLossBrems(const double &mom) const
Returns energy loss.
Definition: GFMaterialEffects.cxx:483
genf::GFMaterialEffects::calcBeta
void calcBeta(double mom)
sets fbeta, fgamma, fgammasquare; must only be used after calling getParameters() ...
Definition: GFMaterialEffects.cxx:275
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
genf::GFMaterialEffects::getParameters
void getParameters()
contains energy loss classes
Definition: GFMaterialEffects.cxx:254
X
Float_t X
Definition: plot.C:39
genf::GFMaterialEffects::energyLossBetheBloch
double energyLossBetheBloch(const double &mom)
Returns energy loss.
Definition: GFMaterialEffects.cxx:289
double genf::GFMaterialEffects::stepper ( const double &  maxDist,
const TVector3 &  pos,
const TVector3 &  dir,
const double &  mom,
const int &  pdg 
)
inline

Definition at line 98 of file GFMaterialEffects.h.

References calcBeta(), energyLossBetheBloch(), energyLossBrems(), getParameters(), MeanExcEnergy_get(), noise(), noiseBetheBloch(), noiseBrems(), noiseCoulomb(), stepper(), and Z.

102  {
103  return stepper(maxDist, pos.X(),pos.Y(),pos.Z(),dir.X(),dir.Y(),dir.Z(),mom,pdg);};
dir
TDirectory * dir
Definition: macro.C:5
genf::GFMaterialEffects::stepper
double stepper(const double &maxDist, const double &posx, const double &posy, const double &posz, const double &dirx, const double &diry, const double &dirz, const double &mom, const int &pdg)
Returns maximum length so that a specified momentum loss will not be exceeded.
Definition: GFMaterialEffects.cxx:182

Member Data Documentation

double genf::GFMaterialEffects::fbeta
private

Definition at line 358 of file GFMaterialEffects.h.

Referenced by calcBeta(), energyLossBetheBloch(), energyLossBrems(), noiseBetheBloch(), and noiseCoulomb().

double genf::GFMaterialEffects::fcharge
private

Definition at line 370 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), getParameters(), and noiseBetheBloch().

double genf::GFMaterialEffects::fdedx
private

Definition at line 359 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), and noiseBetheBloch().

bool genf::GFMaterialEffects::fEnergyLossBetheBloch
private

Definition at line 347 of file GFMaterialEffects.h.

Referenced by effects(), setEnergyLossBetheBloch(), and stepper().

bool genf::GFMaterialEffects::fEnergyLossBrems
private

Definition at line 350 of file GFMaterialEffects.h.

Referenced by effects(), setEnergyLossBrems(), and stepper().

double genf::GFMaterialEffects::fgamma
private

Definition at line 360 of file GFMaterialEffects.h.

Referenced by calcBeta(), and noiseBetheBloch().

double genf::GFMaterialEffects::fgammaSquare
private

Definition at line 361 of file GFMaterialEffects.h.

Referenced by calcBeta(), energyLossBetheBloch(), and noiseBetheBloch().

genf::GFMaterialEffects * genf::GFMaterialEffects::finstance = NULL
staticprivate

Definition at line 63 of file GFMaterialEffects.h.

Referenced by destruct(), and getInstance().

double genf::GFMaterialEffects::fmass
private

Definition at line 371 of file GFMaterialEffects.h.

Referenced by calcBeta(), energyLossBetheBloch(), energyLossBrems(), getParameters(), noiseBetheBloch(), and noiseBrems().

double genf::GFMaterialEffects::fmatA
private

Definition at line 365 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), energyLossBrems(), getParameters(), and noiseBetheBloch().

double genf::GFMaterialEffects::fmatDensity
private

Definition at line 363 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), energyLossBrems(), getParameters(), and noiseBetheBloch().

double genf::GFMaterialEffects::fmatZ
private

Definition at line 364 of file GFMaterialEffects.h.

Referenced by effects(), energyLossBetheBloch(), energyLossBrems(), getParameters(), noiseBetheBloch(), noiseCoulomb(), and stepper().

double genf::GFMaterialEffects::fmEE
private

Definition at line 367 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), and getParameters().

bool genf::GFMaterialEffects::fNoiseBetheBloch
private

Definition at line 348 of file GFMaterialEffects.h.

Referenced by effects(), and setNoiseBetheBloch().

bool genf::GFMaterialEffects::fNoiseBrems
private

Definition at line 351 of file GFMaterialEffects.h.

Referenced by effects(), and setNoiseBrems().

bool genf::GFMaterialEffects::fNoiseCoulomb
private

Definition at line 349 of file GFMaterialEffects.h.

Referenced by effects(), and setNoiseCoulomb().

int genf::GFMaterialEffects::fpdg
private

Definition at line 369 of file GFMaterialEffects.h.

Referenced by effects(), energyLossBrems(), getParameters(), and noiseBrems().

double genf::GFMaterialEffects::fradiationLength
private

Definition at line 366 of file GFMaterialEffects.h.

Referenced by getParameters(), noiseBrems(), and noiseCoulomb().

double genf::GFMaterialEffects::fstep
private

Definition at line 355 of file GFMaterialEffects.h.

Referenced by effects(), energyLossBetheBloch(), energyLossBrems(), noiseBetheBloch(), noiseBrems(), noiseCoulomb(), and stepper().

const double genf::GFMaterialEffects::me
private

Definition at line 353 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), energyLossBrems(), and noiseBetheBloch().

The documentation for this class was generated from the following files: