LArSoft  v09_90_00
Liquid Argon Software toolkit - https://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 50 of file GFMaterialEffects.h.

Constructor & Destructor Documentation

genf::GFMaterialEffects::GFMaterialEffects ( )
private

Definition at line 55 of file GFMaterialEffects.cxx.

Referenced by getInstance().

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

Definition at line 38 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().

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

Definition at line 83 of file GFMaterialEffects.cxx.

References finstance.

84 {
85  if (finstance != NULL) {
86  delete finstance;
87  finstance = NULL;
88  }
89 }
genf::GFMaterialEffects::finstance
static GFMaterialEffects * finstance
Definition: GFMaterialEffects.h:54
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 91 of file GFMaterialEffects.cxx.

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

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

100 {
101 
102  //assert(points.size()==pointPaths.size());
103  fpdg = pdg;
104 
105  double momLoss = 0.;
106 
107  for (unsigned int i = 1; i < points.size(); ++i) {
108  // TVector3 p1=points.at(i-1);
109  //TVector3 p2=points.at(i);
110  //TVector3 dir=p2-p1;
111  TVector3 dir = points.at(i) - points.at(i - 1);
112  double dist = dir.Mag();
113  double realPath = pointPaths.at(i);
114 
115  if (dist > 1.E-8) { // do material effects only if distance is not too small
116  dir *= 1. / dist; //normalize dir
117 
118  double X(0.);
119  /*
120  double matDensity, matZ, matA, radiationLength, mEE;
121  double step;
122  */
123 
124  gGeoManager->InitTrack(points.at(i - 1).X(),
125  points.at(i - 1).Y(),
126  points.at(i - 1).Z(),
127  dir.X(),
128  dir.Y(),
129  dir.Z());
130 
131  while (X < dist) {
132 
133  getParameters();
134  // geoMatManager->getMaterialParameters(matDensity, matZ, matA, radiationLength, mEE);
135 
136  // step = geoMatManager->stepOrNextBoundary(dist-X);
137  gGeoManager->FindNextBoundaryAndStep(dist - X);
138  fstep = gGeoManager->GetStep();
139 
140  // Loop over EnergyLoss classes
141  if (fmatZ > 1.E-3) {
142  /*
143  for(unsigned int j=0;j<fEnergyLoss.size();++j){
144  momLoss += realPath/dist*fEnergyLoss.at(j)->energyLoss(step,
145  mom,
146  pdg,
147  matDensity,
148  matZ,
149  matA,
150  radiationLength,
151  mEE,
152  doNoise,
153  noise,
154  jacobian,
155  directionBefore,
156  directionAfter);
157  }
158  */
159  calcBeta(mom);
160 
161  if (fEnergyLossBetheBloch) momLoss += realPath / dist * this->energyLossBetheBloch(mom);
162  if (doNoise && fEnergyLossBetheBloch && fNoiseBetheBloch)
163  this->noiseBetheBloch(mom, noise);
164 
165  if (/*doNoise &&*/ fNoiseCoulomb) // Force it
166  this->noiseCoulomb(mom, noise, jacobian, directionBefore, directionAfter);
167 
168  if (fEnergyLossBrems) momLoss += realPath / dist * this->energyLossBrems(mom);
169  if (doNoise && fEnergyLossBrems && fNoiseBrems) this->noiseBrems(mom, noise);
170  }
171 
172  X += fstep;
173  }
174  }
175  }
176  //std::cout << "GFMaterialEffects::effects(): momLoss " << momLoss << std::endl;
177  return momLoss;
178 }
genf::GFMaterialEffects::noiseBrems
void noiseBrems(const double &mom, TMatrixT< double > *noise) const
calculation of energy loss straggeling
Definition: GFMaterialEffects.cxx:717
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:394
genf::GFMaterialEffects::energyLossBrems
double energyLossBrems(const double &mom) const
Returns energy loss.
Definition: GFMaterialEffects.cxx:508
genf::GFMaterialEffects::noiseBetheBloch
void noiseBetheBloch(const double &mom, TMatrixT< double > *noise) const
calculation of energy loss straggling
Definition: GFMaterialEffects.cxx:329
noise
unsigned int noise()
Definition: chem4.cc:261
E
Float_t E
Definition: plot.C:20
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
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:537
genf::GFMaterialEffects::getParameters
void getParameters()
contains energy loss classes
Definition: GFMaterialEffects.cxx:251
X
Float_t X
Definition: plot.C:37
genf::GFMaterialEffects::energyLossBetheBloch
double energyLossBetheBloch(const double &mom)
Returns energy loss.
Definition: GFMaterialEffects.cxx:290
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 290 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

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

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

Referenced by effects(), and stepper().

509 {
510 
511  if (fabs(fpdg) != 11) return 0; // only for electrons and positrons
512 
513 #if !defined(BETHE)
514  static const double C[101] = {
515  0.0, -0.960613E-01, 0.631029E-01, -0.142819E-01, 0.150437E-02, -0.733286E-04,
516  0.131404E-05, 0.859343E-01, -0.529023E-01, 0.131899E-01, -0.159201E-02, 0.926958E-04,
517  -0.208439E-05, -0.684096E+01, 0.370364E+01, -0.786752E+00, 0.822670E-01, -0.424710E-02,
518  0.867980E-04, -0.200856E+01, 0.129573E+01, -0.306533E+00, 0.343682E-01, -0.185931E-02,
519  0.392432E-04, 0.127538E+01, -0.515705E+00, 0.820644E-01, -0.641997E-02, 0.245913E-03,
520  -0.365789E-05, 0.115792E+00, -0.463143E-01, 0.725442E-02, -0.556266E-03, 0.208049E-04,
521  -0.300895E-06, -0.271082E-01, 0.173949E-01, -0.452531E-02, 0.569405E-03, -0.344856E-04,
522  0.803964E-06, 0.419855E-02, -0.277188E-02, 0.737658E-03, -0.939463E-04, 0.569748E-05,
523  -0.131737E-06, -0.318752E-03, 0.215144E-03, -0.579787E-04, 0.737972E-05, -0.441485E-06,
524  0.994726E-08, 0.938233E-05, -0.651642E-05, 0.177303E-05, -0.224680E-06, 0.132080E-07,
525  -0.288593E-09, -0.245667E-03, 0.833406E-04, -0.129217E-04, 0.915099E-06, -0.247179E-07,
526  0.147696E-03, -0.498793E-04, 0.402375E-05, 0.989281E-07, -0.133378E-07, -0.737702E-02,
527  0.333057E-02, -0.553141E-03, 0.402464E-04, -0.107977E-05, -0.641533E-02, 0.290113E-02,
528  -0.477641E-03, 0.342008E-04, -0.900582E-06, 0.574303E-05, 0.908521E-04, -0.256900E-04,
529  0.239921E-05, -0.741271E-07, -0.341260E-04, 0.971711E-05, -0.172031E-06, -0.119455E-06,
530  0.704166E-08, 0.341740E-05, -0.775867E-06, -0.653231E-07, 0.225605E-07, -0.114860E-08,
531  -0.119391E-06, 0.194885E-07, 0.588959E-08, -0.127589E-08, 0.608247E-10};
532  static const double xi = 2.51, beta = 0.99, vl = 0.00004;
533 #endif
534 #if defined(BETHE) // no MIGDAL corrections
535  static const double C[101] = {
536  0.0, 0.834459E-02, 0.443979E-02, -0.101420E-02, 0.963240E-04, -0.409769E-05,
537  0.642589E-07, 0.464473E-02, -0.290378E-02, 0.547457E-03, -0.426949E-04, 0.137760E-05,
538  -0.131050E-07, -0.547866E-02, 0.156218E-02, -0.167352E-03, 0.101026E-04, -0.427518E-06,
539  0.949555E-08, -0.406862E-02, 0.208317E-02, -0.374766E-03, 0.317610E-04, -0.130533E-05,
540  0.211051E-07, 0.158941E-02, -0.385362E-03, 0.315564E-04, -0.734968E-06, -0.230387E-07,
541  0.971174E-09, 0.467219E-03, -0.154047E-03, 0.202400E-04, -0.132438E-05, 0.431474E-07,
542  -0.559750E-09, -0.220958E-02, 0.100698E-02, -0.596464E-04, -0.124653E-04, 0.142999E-05,
543  -0.394378E-07, 0.477447E-03, -0.184952E-03, -0.152614E-04, 0.848418E-05, -0.736136E-06,
544  0.190192E-07, -0.552930E-04, 0.209858E-04, 0.290001E-05, -0.133254E-05, 0.116971E-06,
545  -0.309716E-08, 0.212117E-05, -0.103884E-05, -0.110912E-06, 0.655143E-07, -0.613013E-08,
546  0.169207E-09, 0.301125E-04, -0.461920E-04, 0.871485E-05, -0.622331E-06, 0.151800E-07,
547  -0.478023E-04, 0.247530E-04, -0.381763E-05, 0.232819E-06, -0.494487E-08, -0.336230E-04,
548  0.223822E-04, -0.384583E-05, 0.252867E-06, -0.572599E-08, 0.105335E-04, -0.567074E-06,
549  -0.216564E-06, 0.237268E-07, -0.658131E-09, 0.282025E-05, -0.671965E-06, 0.565858E-07,
550  -0.193843E-08, 0.211839E-10, 0.157544E-04, -0.304104E-05, -0.624410E-06, 0.120124E-06,
551  -0.457445E-08, -0.188222E-05, -0.407118E-06, 0.375106E-06, -0.466881E-07, 0.158312E-08,
552  0.945037E-07, 0.564718E-07, -0.319231E-07, 0.371926E-08, -0.123111E-09};
553  static const double xi = 2.10, fbeta = 1.00, vl = 0.001;
554 #endif
555 
556  double BCUT = 10000.; // energy up to which soft bremsstrahlung energy loss is calculated
557 
558  double THIGH = 100., CHIGH = 50.;
559  double dedxBrems = 0.;
560 
561  if (BCUT > 0.) {
562  double T, kc;
563 
564  if (BCUT >= mom) BCUT = mom; // confine BCUT to mom
565 
566  // T=mom, confined to THIGH
567  // kc=BCUT, confined to CHIGH ??
568  if (mom >= THIGH) {
569  T = THIGH;
570  if (BCUT >= THIGH)
571  kc = CHIGH;
572  else
573  kc = BCUT;
574  }
575  else {
576  T = mom;
577  kc = BCUT;
578  }
579 
580  double E = T + me; // total electron energy
581  if (BCUT > T) kc = T;
582 
583  double X = log(T / me);
584  double Y = log(kc / (E * vl));
585 
586  double XX;
587  int K;
588  double S = 0., YY = 1.;
589 
590  for (unsigned int I = 1; I <= 2; ++I) {
591  XX = 1.;
592  for (unsigned int J = 1; J <= 6; ++J) {
593  K = 6 * I + J - 6;
594  S = S + C[K] * XX * YY;
595  XX = XX * X;
596  }
597  YY = YY * Y;
598  }
599 
600  for (unsigned int I = 3; I <= 6; ++I) {
601  XX = 1.;
602  for (unsigned int J = 1; J <= 6; ++J) {
603  K = 6 * I + J - 6;
604  if (Y <= 0.)
605  S = S + C[K] * XX * YY;
606  else
607  S = S + C[K + 24] * XX * YY;
608  XX = XX * X;
609  }
610  YY = YY * Y;
611  }
612 
613  double SS = 0.;
614  YY = 1.;
615 
616  for (unsigned int I = 1; I <= 2; ++I) {
617  XX = 1.;
618  for (unsigned int J = 1; J <= 5; ++J) {
619  K = 5 * I + J + 55;
620  SS = SS + C[K] * XX * YY;
621  XX = XX * X;
622  }
623  YY = YY * Y;
624  }
625 
626  for (unsigned int I = 3; I <= 5; ++I) {
627  XX = 1.;
628  for (unsigned int J = 1; J <= 5; ++J) {
629  K = 5 * I + J + 55;
630  if (Y <= 0.)
631  SS = SS + C[K] * XX * YY;
632  else
633  SS = SS + C[K + 15] * XX * YY;
634  XX = XX * X;
635  }
636  YY = YY * Y;
637  }
638 
639  S = S + fmatZ * SS;
640 
641  if (S > 0.) {
642  double CORR = 1.;
643 #if !defined(BETHE)
644  CORR = 1. / (1. + 0.805485E-10 * fmatDensity * fmatZ * E * E /
645  (fmatA * kc * kc)); // MIGDAL correction factor
646 #endif
647 
648  double FAC = fmatZ * (fmatZ + xi) * E * E * pow((kc * CORR / T), beta) / (E + me);
649  if (FAC <= 0.) return 0.;
650  dedxBrems = FAC * S;
651 
652  double RAT;
653 
654  if (mom > THIGH) {
655  if (BCUT < THIGH) {
656  RAT = BCUT / mom;
657  S = (1. - 0.5 * RAT + 2. * RAT * RAT / 9.);
658  RAT = BCUT / T;
659  S = S / (1. - 0.5 * RAT + 2. * RAT * RAT / 9.);
660  }
661  else {
662  RAT = BCUT / mom;
663  S = BCUT * (1. - 0.5 * RAT + 2. * RAT * RAT / 9.);
664  RAT = kc / T;
665  S = S / (kc * (1. - 0.5 * RAT + 2. * RAT * RAT / 9.));
666  }
667  dedxBrems = dedxBrems * S; // GeV barn
668  }
669 
670  dedxBrems = 0.60221367 * fmatDensity * dedxBrems / fmatA; // energy loss dE/dx [GeV/cm]
671  }
672  }
673 
674  if (dedxBrems < 0.) dedxBrems = 0;
675 
676  double factor = 1.; // positron correction factor
677 
678  if (fpdg == -11) {
679  static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054;
680 
681  double ETA = 0.;
682  if (fmatZ > 0.) {
683  double X = log(AA * mom / fmatZ * fmatZ);
684  if (X > -8.) {
685  if (X >= +9.)
686  ETA = 1.;
687  else {
688  double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.);
689  ETA = 0.5 + atan(W) / M_PI;
690  }
691  }
692  }
693 
694  double E0;
695 
696  if (ETA < 0.0001)
697  factor = 1.E-10;
698  else if (ETA > 0.9999)
699  factor = 1.;
700  else {
701  E0 = BCUT / mom;
702  if (E0 > 1.) E0 = 1.;
703  if (E0 < 1.E-8)
704  factor = 1.;
705  else
706  factor = ETA * (1. - pow(1. - E0, 1. / ETA)) / E0;
707  }
708  }
709 
710  double DE = fstep * factor * dedxBrems; //always positive
711  double momLoss =
712  sqrt(mom * mom + 2. * sqrt(mom * mom + fmass * fmass) * DE + DE * DE) - mom; //always positive
713 
714  return momLoss;
715 }
Y
Float_t Y
Definition: plot.C:37
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
X
Float_t X
Definition: plot.C:37
genf::GFMaterialEffects * genf::GFMaterialEffects::getInstance ( )
static

Definition at line 77 of file GFMaterialEffects.cxx.

References finstance, and GFMaterialEffects().

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

78 {
79  if (finstance == NULL) finstance = new GFMaterialEffects();
80  return finstance;
81 }
genf::GFMaterialEffects::finstance
static GFMaterialEffects * finstance
Definition: GFMaterialEffects.h:54
void genf::GFMaterialEffects::getParameters ( )
private

contains energy loss classes

interface to material and geometry

Definition at line 251 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

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

Definition at line 747 of file GFMaterialEffects.cxx.

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

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

748 {
749  if ((Z < 0) || (Z > MeanExcEnergy_NELEMENTS))
750  throw GFException(std::string(__func__) + ": unsupported Z", __LINE__, __FILE__).setFatal();
751  return MeanExcEnergy_vals[Z];
752 }
MeanExcEnergy_NELEMENTS
const int MeanExcEnergy_NELEMENTS
Definition: GFMaterialEffects.cxx:737
Z
Float_t Z
Definition: plot.C:37
MeanExcEnergy_vals
const float MeanExcEnergy_vals[MeanExcEnergy_NELEMENTS]
Definition: GFMaterialEffects.cxx:738
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
double genf::GFMaterialEffects::MeanExcEnergy_get ( TGeoMaterial *  mat)
private

Definition at line 754 of file GFMaterialEffects.cxx.

References MeanExcEnergy_get().

755 {
756  if (mat->IsMixture()) {
757  double logMEE = 0.;
758  double denom = 0.;
759  TGeoMixture* mix = (TGeoMixture*)mat;
760  for (int i = 0; i < mix->GetNelements(); ++i) {
761  int index = int(floor((mix->GetZmixt())[i]));
762  logMEE += 1. / (mix->GetAmixt())[i] * (mix->GetWmixt())[i] * (mix->GetZmixt())[i] *
763  log(MeanExcEnergy_get(index));
764  denom += (mix->GetWmixt())[i] * (mix->GetZmixt())[i] * 1. / (mix->GetAmixt())[i];
765  }
766  logMEE /= denom;
767  return exp(logMEE);
768  }
769  else { // not a mixture
770  int index = int(floor(mat->GetZ()));
771  return MeanExcEnergy_get(index);
772  }
773 }
mat
Float_t mat
Definition: plot.C:38
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 329 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

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

References fmass, fpdg, fradiationLength, and fstep.

Referenced by effects(), and stepper().

718 {
719 
720  if (fabs(fpdg) != 11) return; // only for electrons and positrons
721 
722  double LX = 1.442695 * fstep / fradiationLength;
723  double S2B = mom * mom * (1. / pow(3., LX) - 1. / pow(4., LX));
724  double DEDXB = pow(fabs(S2B), 0.5);
725  DEDXB = 1.2E9 * DEDXB; //eV
726  double sigma2E = DEDXB * DEDXB; //eV^2
727  sigma2E *= 1.E-18; // eV -> GeV
728 
729  (*noise)[6][6] += (mom * mom + fmass * fmass) / pow(mom, 6.) * sigma2E;
730 }
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 394 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

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

Definition at line 60 of file GFMaterialEffects.h.

References fEnergyLossBetheBloch.

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

Definition at line 63 of file GFMaterialEffects.h.

References fEnergyLossBrems.

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

Definition at line 61 of file GFMaterialEffects.h.

References fNoiseBetheBloch.

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

Definition at line 64 of file GFMaterialEffects.h.

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

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

Definition at line 62 of file GFMaterialEffects.h.

References fNoiseCoulomb.

62 { 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 180 of file GFMaterialEffects.cxx.

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

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

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

Definition at line 89 of file GFMaterialEffects.h.

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

94  {
95  return stepper(maxDist, pos.X(), pos.Y(), pos.Z(), dir.X(), dir.Y(), dir.Z(), mom, pdg);
96  };
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:180

Member Data Documentation

double genf::GFMaterialEffects::fbeta
private

Definition at line 349 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fcharge
private

Definition at line 361 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fdedx
private

Definition at line 350 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), and noiseBetheBloch().

bool genf::GFMaterialEffects::fEnergyLossBetheBloch
private

Definition at line 338 of file GFMaterialEffects.h.

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

bool genf::GFMaterialEffects::fEnergyLossBrems
private

Definition at line 341 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fgamma
private

Definition at line 351 of file GFMaterialEffects.h.

Referenced by calcBeta(), and noiseBetheBloch().

double genf::GFMaterialEffects::fgammaSquare
private

Definition at line 352 of file GFMaterialEffects.h.

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

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

Definition at line 54 of file GFMaterialEffects.h.

Referenced by destruct(), and getInstance().

double genf::GFMaterialEffects::fmass
private

Definition at line 362 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fmatA
private

Definition at line 356 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fmatDensity
private

Definition at line 354 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fmatZ
private

Definition at line 355 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fmEE
private

Definition at line 358 of file GFMaterialEffects.h.

Referenced by energyLossBetheBloch(), and getParameters().

bool genf::GFMaterialEffects::fNoiseBetheBloch
private

Definition at line 339 of file GFMaterialEffects.h.

Referenced by effects(), and setNoiseBetheBloch().

bool genf::GFMaterialEffects::fNoiseBrems
private

Definition at line 342 of file GFMaterialEffects.h.

Referenced by effects(), and setNoiseBrems().

bool genf::GFMaterialEffects::fNoiseCoulomb
private

Definition at line 340 of file GFMaterialEffects.h.

Referenced by effects(), and setNoiseCoulomb().

int genf::GFMaterialEffects::fpdg
private

Definition at line 360 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fradiationLength
private

Definition at line 357 of file GFMaterialEffects.h.

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

double genf::GFMaterialEffects::fstep
private

Definition at line 346 of file GFMaterialEffects.h.

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

const double genf::GFMaterialEffects::me
private

Definition at line 344 of file GFMaterialEffects.h.

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

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