SND@LHC Software
Loading...
Searching...
No Matches
genfit::MaterialEffects Class Reference

Stepper and energy loss/noise matrix calculation. More...

#include <MaterialEffects.h>

Collaboration diagram for genfit::MaterialEffects:

Public Member Functions

void init (AbsMaterialInterface *matIfc)
 set the material interface here. Material interface classes must be derived from AbsMaterialInterface.
 
bool isInitialized ()
 
void setNoEffects (bool opt=true)
 
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)
 
void ignoreBoundariesBetweenEqualMaterials (bool opt=true)
 
void setMscModel (const std::string &modelName)
 Select the multiple scattering model that will be used during track fit.
 
double effects (const std::vector< genfit::RKStep > &steps, int materialsFXStart, int materialsFXStop, const double &mom, const int &pdg, M7x7 *noise=nullptr)
 Calculates energy loss in the traveled path, optional calculation of noise matrix.
 
void stepper (const RKTrackRep *rep, M1x7 &state7, const double &mom, double &relMomLoss, const int &pdg, MaterialProperties &currentMaterial, StepLimits &limits, bool varField=true)
 Returns maximum length so that a specified momentum loss will not be exceeded.
 
void init (AbsMaterialInterface *matIfc)
 set the material interface here. Material interface classes must be derived from AbsMaterialInterface.
 
bool isInitialized ()
 
void setNoEffects (bool opt=true)
 
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)
 
void ignoreBoundariesBetweenEqualMaterials (bool opt=true)
 
void setMscModel (const std::string &modelName)
 Select the multiple scattering model that will be used during track fit.
 
double effects (const std::vector< genfit::RKStep > &steps, int materialsFXStart, int materialsFXStop, const double &mom, const int &pdg, M7x7 *noise=nullptr)
 Calculates energy loss in the traveled path, optional calculation of noise matrix.
 
void stepper (const RKTrackRep *rep, M1x7 &state7, const double &mom, double &relMomLoss, const int &pdg, MaterialProperties &currentMaterial, StepLimits &limits, bool varField=true)
 Returns maximum length so that a specified momentum loss will not be exceeded.
 

Static Public Member Functions

static MaterialEffectsgetInstance ()
 
static void destruct ()
 
static MaterialEffectsgetInstance ()
 
static void destruct ()
 

Private Member Functions

 MaterialEffects ()
 
virtual ~MaterialEffects ()
 
void getParticleParameters (double mom)
 sets charge_, mass_ and calculates beta_, gamma_, fgammasquare;
 
double energyLossBetheBloch ()
 Returns energy loss.
 
void noiseBetheBloch (M7x7 &noise) const
 calculation of energy loss straggeling
 
void noiseCoulomb (M7x7 &noise, const M1x3 &direction) const
 calculation of multiple scattering
 
double energyLossBrems () const
 Returns energy loss.
 
void noiseBrems (M7x7 &noise) const
 calculation of energy loss straggeling
 
 MaterialEffects ()
 
virtual ~MaterialEffects ()
 
void getParticleParameters (double mom)
 sets charge_, mass_ and calculates beta_, gamma_, fgammasquare;
 
double energyLossBetheBloch ()
 Returns energy loss.
 
void noiseBetheBloch (M7x7 &noise) const
 calculation of energy loss straggeling
 
void noiseCoulomb (M7x7 &noise, const M1x3 &direction) const
 calculation of multiple scattering
 
double energyLossBrems () const
 Returns energy loss.
 
void noiseBrems (M7x7 &noise) const
 calculation of energy loss straggeling
 

Private Attributes

bool noEffects_
 
bool energyLossBetheBloch_
 
bool noiseBetheBloch_
 
bool noiseCoulomb_
 
bool energyLossBrems_
 
bool noiseBrems_
 
bool ignoreBoundariesBetweenEqualMaterials_
 
const double me_
 
double stepSize_
 
double mom_
 
double beta_
 
double dEdx_
 
double gamma_
 
double gammaSquare_
 
double matDensity_
 
double matZ_
 
double matA_
 
double radiationLength_
 
double mEE_
 
int pdg_
 
int charge_
 
double mass_
 
int mscModelCode_
 
AbsMaterialInterfacematerialInterface_
 depending on this number a specific msc model is chosen in the noiseCoulomb function.
 

Static Private Attributes

static MaterialEffectsinstance_ = nullptr
 

Detailed Description

Stepper and energy loss/noise matrix calculation.

Author
Christian Höppner (Technische Universität München, original author)
Sebastian Neubert (Technische Universität München, original author)
Johannes Rauch (Technische Universität München, author)

It provides functionality to limit the stepsize of an extrapolation in order not to exceed a specified maximum momentum loss. After propagation, the energy loss for the given length and (optionally) the noise matrix can be calculated. You have to set which energy-loss and noise mechanisms you want to use. At the moment, per default all energy loss and noise options are ON.

Definition at line 51 of file MaterialEffects.h.

Constructor & Destructor Documentation

◆ MaterialEffects() [1/2]

genfit::MaterialEffects::MaterialEffects ( )
private

Definition at line 41 of file MaterialEffects.cc.

41 :
42 noEffects_(false),
44 noiseCoulomb_(true),
45 energyLossBrems_(true), noiseBrems_(true),
47 me_(0.510998910E-3),
48 stepSize_(0),
49 mom_(0),
50 beta_(0),
51 dEdx_(0),
52 gamma_(0),
53 gammaSquare_(0),
54 matDensity_(0),
55 matZ_(0),
56 matA_(0),
58 mEE_(0),
59 pdg_(0),
60 charge_(0),
61 mass_(0),
63 materialInterface_(nullptr)
64{
65}
AbsMaterialInterface * materialInterface_
depending on this number a specific msc model is chosen in the noiseCoulomb function.

◆ ~MaterialEffects() [1/2]

genfit::MaterialEffects::~MaterialEffects ( )
privatevirtual

Definition at line 67 of file MaterialEffects.cc.

68{
69 if (materialInterface_ != nullptr) delete materialInterface_;
70}

◆ MaterialEffects() [2/2]

genfit::MaterialEffects::MaterialEffects ( )
private

◆ ~MaterialEffects() [2/2]

virtual genfit::MaterialEffects::~MaterialEffects ( )
privatevirtual

Member Function Documentation

◆ destruct() [1/2]

void genfit::MaterialEffects::destruct ( )
static

Definition at line 78 of file MaterialEffects.cc.

79{
80 if (instance_ != nullptr) {
81 delete instance_;
82 instance_ = nullptr;
83 }
84}
static MaterialEffects * instance_

◆ destruct() [2/2]

static void genfit::MaterialEffects::destruct ( )
static

◆ effects() [1/2]

double genfit::MaterialEffects::effects ( const std::vector< genfit::RKStep > &  steps,
int  materialsFXStart,
int  materialsFXStop,
const double &  mom,
const int &  pdg,
M7x7 noise = nullptr 
)

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

◆ effects() [2/2]

double genfit::MaterialEffects::effects ( const std::vector< genfit::RKStep > &  steps,
int  materialsFXStart,
int  materialsFXStop,
const double &  mom,
const int &  pdg,
M7x7 noise = nullptr 
)

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

◆ energyLossBetheBloch() [1/2]

double genfit::MaterialEffects::energyLossBetheBloch ( )
private

Returns energy loss.

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

Definition at line 340 of file MaterialEffects.cc.

341{
342
343 // calc dEdx_, also needed in noiseBetheBloch!
344 dEdx_ = 0.307075 * matZ_ / matA_ * matDensity_ / (beta_ * beta_) * charge_ * charge_;
345 double massRatio = me_ / mass_;
346 double argument = gammaSquare_ * beta_ * beta_ * me_ * 1.E3 * 2. / ((1.E-6 * mEE_) * sqrt(1 + 2 * gamma_ * massRatio + massRatio * massRatio));
347 dEdx_ *= log(argument) - beta_ * beta_; // Bethe-Bloch [MeV/cm]
348 dEdx_ *= 1.E-3; // in GeV/cm, hence 1.e-3
349 if (dEdx_ < 0.) dEdx_ = 0;
350
351 double dE = fabs(stepSize_) * dEdx_; //always positive
352 double momLoss = sqrt(mom_ * mom_ + 2.*sqrt(mom_ * mom_ + mass_ * mass_) * dE + dE * dE) - mom_; //always positive
353
354 //in vacuum it can numerically happen that momLoss becomes a small negative number.
355 if (momLoss < 0.) return 0.;
356 return momLoss;
357}

◆ energyLossBetheBloch() [2/2]

double genfit::MaterialEffects::energyLossBetheBloch ( )
private

Returns energy loss.

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

◆ energyLossBrems() [1/2]

double genfit::MaterialEffects::energyLossBrems ( ) 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 492 of file MaterialEffects.cc.

493{
494
495 // Code ported from GEANT 3
496
497 if (abs(pdg_) != 11) return 0; // only for electrons and positrons
498
499#if !defined(BETHE)
500 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};
501 static const double xi = 2.51, beta = 0.99, vl = 0.00004;
502#endif
503#if defined(BETHE) // no MIGDAL corrections
504 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};
505 static const double xi = 2.10, beta = 1.00, vl = 0.001;
506#endif
507
508 double BCUT = 10000.; // energy up to which soft bremsstrahlung energy loss is calculated
509
510 static const double THIGH = 100., CHIGH = 50.;
511 double dedxBrems = 0.;
512
513 if (BCUT > 0.) {
514 double T, kc;
515
516 if (BCUT >= mom_) BCUT = mom_; // confine BCUT to mom_
517
518 // T=mom_, confined to THIGH
519 // kc=BCUT, confined to CHIGH ??
520 if (mom_ >= THIGH) {
521 T = THIGH;
522 if (BCUT >= THIGH) kc = CHIGH;
523 else kc = BCUT;
524 } else {
525 T = mom_;
526 kc = BCUT;
527 }
528
529 double E = T + me_; // total electron energy
530 if (BCUT > T) kc = T;
531
532 double X = log(T / me_);
533 double Y = log(kc / (E * vl));
534
535 double XX;
536 int K;
537 double S = 0., YY = 1.;
538
539 for (unsigned int I = 1; I <= 2; ++I) {
540 XX = 1.;
541 for (unsigned int J = 1; J <= 6; ++J) {
542 K = 6 * I + J - 6;
543 S = S + C[K] * XX * YY;
544 XX = XX * X;
545 }
546 YY = YY * Y;
547 }
548
549 for (unsigned int I = 3; I <= 6; ++I) {
550 XX = 1.;
551 for (unsigned int J = 1; J <= 6; ++J) {
552 K = 6 * I + J - 6;
553 if (Y <= 0.) S = S + C[K] * XX * YY;
554 else S = S + C[K + 24] * XX * YY;
555 XX = XX * X;
556 }
557 YY = YY * Y;
558 }
559
560 double SS = 0.;
561 YY = 1.;
562
563 for (unsigned int I = 1; I <= 2; ++I) {
564 XX = 1.;
565 for (unsigned int J = 1; J <= 5; ++J) {
566 K = 5 * I + J + 55;
567 SS = SS + C[K] * XX * YY;
568 XX = XX * X;
569 }
570 YY = YY * Y;
571 }
572
573 for (unsigned int I = 3; I <= 5; ++I) {
574 XX = 1.;
575 for (unsigned int J = 1; J <= 5; ++J) {
576 K = 5 * I + J + 55;
577 if (Y <= 0.) SS = SS + C[K] * XX * YY;
578 else SS = SS + C[K + 15] * XX * YY;
579 XX = XX * X;
580 }
581 YY = YY * Y;
582 }
583
584 S = S + matZ_ * SS;
585
586 if (S > 0.) {
587 double CORR = 1.;
588#if !defined(BETHE)
589 CORR = 1. / (1. + 0.805485E-10 * matDensity_ * matZ_ * E * E / (matA_ * kc * kc)); // MIGDAL correction factor
590#endif
591
592 // We use exp(beta * log(...) here because pow(..., beta) is
593 // REALLY slow and we don't need ultimate numerical precision
594 // for this approximation.
595 double FAC = matZ_ * (matZ_ + xi) * E * E / (E + me_);
596 if (beta == 1.) { // That is the #ifdef BETHE case
597 FAC *= kc * CORR / T;
598 } else {
599 FAC *= exp(beta * log(kc * CORR / T));
600 }
601 if (FAC <= 0.) return 0.;
602 dedxBrems = FAC * S;
603
604
605 if (mom_ > THIGH) {
606 double RAT;
607 if (BCUT < THIGH) {
608 RAT = BCUT / mom_;
609 S = (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
610 RAT = BCUT / T;
611 S = S / (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
612 } else {
613 RAT = BCUT / mom_;
614 S = BCUT * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
615 RAT = kc / T;
616 S = S / (kc * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.));
617 }
618 dedxBrems = dedxBrems * S; // GeV barn
619 }
620
621 dedxBrems = 0.60221367 * matDensity_ * dedxBrems / matA_; // energy loss dE/dx [GeV/cm]
622 }
623 }
624
625 if (dedxBrems < 0.) dedxBrems = 0;
626
627 double factor = 1.; // positron correction factor
628
629 if (pdg_ == -11) {
630 static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054;
631
632 double ETA = 0.;
633 if (matZ_ > 0.) {
634 double X = log(AA * mom_ / (matZ_ * matZ_));
635 if (X > -8.) {
636 if (X >= +9.) ETA = 1.;
637 else {
638 double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.);
639 ETA = 0.5 + atan(W) / M_PI;
640 }
641 }
642 }
643
644 if (ETA < 0.0001) factor = 1.E-10;
645 else if (ETA > 0.9999) factor = 1.;
646 else {
647 double E0 = BCUT / mom_;
648 if (E0 > 1.) E0 = 1.;
649 if (E0 < 1.E-8) factor = 1.;
650 else factor = ETA * (1. - pow(1. - E0, 1. / ETA)) / E0;
651 }
652 }
653
654 double dE = fabs(stepSize_) * factor * dedxBrems; //always positive
655 double momLoss = sqrt(mom_ * mom_ + 2.*sqrt(mom_ * mom_ + mass_ * mass_) * dE + dE * dE) - mom_; //always positive
656
657 return momLoss;
658}
dict S
Definition MufiCTR.py:12

◆ energyLossBrems() [2/2]

double genfit::MaterialEffects::energyLossBrems ( ) 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.

◆ getInstance() [1/2]

MaterialEffects * genfit::MaterialEffects::getInstance ( )
static

Definition at line 72 of file MaterialEffects.cc.

73{
74 if (instance_ == nullptr) instance_ = new MaterialEffects();
75 return instance_;
76}

◆ getInstance() [2/2]

static MaterialEffects * genfit::MaterialEffects::getInstance ( )
static

◆ getParticleParameters() [1/2]

void genfit::MaterialEffects::getParticleParameters ( double  mom)
private

sets charge_, mass_ and calculates beta_, gamma_, fgammasquare;

Definition at line 324 of file MaterialEffects.cc.

325{
326 TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg_);
327 mom_ = mom;
328 charge_ = int(part->Charge() / 3.); // We only ever use the square
329 mass_ = part->Mass(); // GeV
330
331 beta_ = 1 / hypot(mass_ / mom, 1);
332 gammaSquare_ = 1 + mom*mom / mass_ / mass_;
333 gamma_ = hypot(mom / mass_, 1);
334}

◆ getParticleParameters() [2/2]

void genfit::MaterialEffects::getParticleParameters ( double  mom)
private

sets charge_, mass_ and calculates beta_, gamma_, fgammasquare;

◆ ignoreBoundariesBetweenEqualMaterials() [1/2]

void genfit::MaterialEffects::ignoreBoundariesBetweenEqualMaterials ( bool  opt = true)
inline

◆ ignoreBoundariesBetweenEqualMaterials() [2/2]

void genfit::MaterialEffects::ignoreBoundariesBetweenEqualMaterials ( bool  opt = true)
inline

Definition at line 77 of file MaterialEffects.h.

◆ init() [1/2]

void genfit::MaterialEffects::init ( AbsMaterialInterface matIfc)

set the material interface here. Material interface classes must be derived from AbsMaterialInterface.

Definition at line 86 of file MaterialEffects.cc.

87{
88 if (materialInterface_ != nullptr) {
89 std::string msg("MaterialEffects::initMaterialInterface(): Already initialized! ");
90 std::runtime_error err(msg);
91 }
92 materialInterface_ = matIfc;
93}

◆ init() [2/2]

void genfit::MaterialEffects::init ( AbsMaterialInterface matIfc)

set the material interface here. Material interface classes must be derived from AbsMaterialInterface.

◆ isInitialized() [1/2]

bool genfit::MaterialEffects::isInitialized ( )
inline

Definition at line 68 of file MaterialEffects.h.

68{ return materialInterface_ != nullptr; }

◆ isInitialized() [2/2]

bool genfit::MaterialEffects::isInitialized ( )
inline

Definition at line 68 of file MaterialEffects.h.

68{ return materialInterface_ != nullptr; }

◆ noiseBetheBloch() [1/2]

void genfit::MaterialEffects::noiseBetheBloch ( M7x7 noise) const
private

calculation of energy loss straggeling

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

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

Needs dEdx_, which is calculated in energyLossBetheBloch, so it has to be called afterwards!

Definition at line 360 of file MaterialEffects.cc.

361{
362
363 // Code ported from GEANT 3
364
365 // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E);
366 double sigma2E = 0.;
367 double zeta = 153.4E3 * charge_ * charge_ / (beta_ * beta_) * matZ_ / matA_ * matDensity_ * fabs(stepSize_); // eV
368 double Emax = 2.E9 * me_ * beta_ * beta_ * gammaSquare_ / (1. + 2.*gamma_ * me_ / mass_ + (me_ / mass_) * (me_ / mass_)); // eV
369 double kappa = zeta / Emax;
370
371 if (kappa > 0.01) { // Vavilov-Gaussian regime
372 sigma2E += zeta * Emax * (1. - beta_ * beta_ / 2.); // eV^2
373 } else { // Urban/Landau approximation
374 // calculate number of collisions Nc
375 double I = 16. * pow(matZ_, 0.9); // eV
376 double f2 = 0.;
377 if (matZ_ > 2.) f2 = 2. / matZ_;
378 double f1 = 1. - f2;
379 double e2 = 10.*matZ_ * matZ_; // eV
380 double e1 = pow((I / pow(e2, f2)), 1. / f1); // eV
381
382 double mbbgg2 = 2.E9 * mass_ * beta_ * beta_ * gammaSquare_; // eV
383 double Sigma1 = dEdx_ * 1.0E9 * f1 / e1 * (log(mbbgg2 / e1) - beta_ * beta_) / (log(mbbgg2 / I) - beta_ * beta_) * 0.6; // 1/cm
384 double Sigma2 = dEdx_ * 1.0E9 * f2 / e2 * (log(mbbgg2 / e2) - beta_ * beta_) / (log(mbbgg2 / I) - beta_ * beta_) * 0.6; // 1/cm
385 double Sigma3 = dEdx_ * 1.0E9 * Emax / (I * (Emax + I) * log((Emax + I) / I)) * 0.4; // 1/cm
386
387 double Nc = (Sigma1 + Sigma2 + Sigma3) * fabs(stepSize_);
388
389 if (Nc > 50.) { // truncated Landau distribution
390 double sigmaalpha = 15.76;
391 // calculate sigmaalpha (see GEANT3 manual W5013)
392 double RLAMED = -0.422784 - beta_ * beta_ - log(zeta / Emax);
393 double RLAMAX = 0.60715 + 1.1934 * RLAMED + (0.67794 + 0.052382 * RLAMED) * exp(0.94753 + 0.74442 * RLAMED);
394 // from lambda max to sigmaalpha=sigma (empirical polynomial)
395 if (RLAMAX <= 1010.) {
396 sigmaalpha = 1.975560
397 + 9.898841e-02 * RLAMAX
398 - 2.828670e-04 * RLAMAX * RLAMAX
399 + 5.345406e-07 * pow(RLAMAX, 3.)
400 - 4.942035e-10 * pow(RLAMAX, 4.)
401 + 1.729807e-13 * pow(RLAMAX, 5.);
402 } else { sigmaalpha = 1.871887E+01 + 1.296254E-02 * RLAMAX; }
403 // alpha=54.6 corresponds to a 0.9996 maximum cut
404 if (sigmaalpha > 54.6) sigmaalpha = 54.6;
405 sigma2E += sigmaalpha * sigmaalpha * zeta * zeta; // eV^2
406 } else { // Urban model
407 static const double alpha = 0.996;
408 double Ealpha = I / (1. - (alpha * Emax / (Emax + I))); // eV
409 double meanE32 = I * (Emax + I) / Emax * (Ealpha - I); // eV^2
410 sigma2E += fabs(stepSize_) * (Sigma1 * e1 * e1 + Sigma2 * e2 * e2 + Sigma3 * meanE32); // eV^2
411 }
412 }
413
414 sigma2E *= 1.E-18; // eV -> GeV
415
416 // update noise matrix, using linear error propagation from E to q/p
417 noise[6 * 7 + 6] += charge_*charge_/beta_/beta_ / pow(mom_, 4) * sigma2E;
418}

◆ noiseBetheBloch() [2/2]

void genfit::MaterialEffects::noiseBetheBloch ( M7x7 noise) const
private

calculation of energy loss straggeling

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

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

Needs dEdx_, which is calculated in energyLossBetheBloch, so it has to be called afterwards!

◆ noiseBrems() [1/2]

void genfit::MaterialEffects::noiseBrems ( M7x7 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 661 of file MaterialEffects.cc.

662{
663
664 // Code ported from GEANT 3 and simplified
665 // this formula assumes p >> m and therefore p^2 + m^2 = p^2
666 // the factor 1.44 is not in the original Behta Heitler model. It seems to be some empirical correction copied over from some other project
667
668 if (abs(pdg_) != 11) return; // only for electrons and positrons
669
670 double minusXOverLn2 = -1.442695 * fabs(stepSize_) / radiationLength_;
671 double sigma2 = 1.44*(pow(3., minusXOverLn2) - pow(4., minusXOverLn2)) / (mom_*mom_);
672 //XXX std::cout << "breams sigma: " << sigma2E << std::endl;
673 //assert(sigma2 >= 0.0);
674 sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0);
675 noise[6 * 7 + 6] += charge_*charge_/beta_/beta_ / pow(mom_, 4) * sigma2;
676
677}

◆ noiseBrems() [2/2]

void genfit::MaterialEffects::noiseBrems ( M7x7 noise) const
private

calculation of energy loss straggeling

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

◆ noiseCoulomb() [1/2]

void genfit::MaterialEffects::noiseCoulomb ( M7x7 noise,
const M1x3 direction 
) const
private

calculation of multiple scattering

This function first calcuates a MSC variance based on the current material and step length 2 different formulas for the MSC variance are implemeted. One can select the formula via "setMscModel". With the MSC variance and the current direction of the track a full 7D noise matrix is calculated. This noise matrix is the additional noise at the end of fStep in the 7D globa cooridnate system taking even the (co)variances of the position coordinates into account.

Definition at line 421 of file MaterialEffects.cc.

423{
424
425 // MULTIPLE SCATTERING; calculate sigma^2
426 double sigma2 = 0;
427 assert(mscModelCode_ == 0 || mscModelCode_ == 1);
428 const double step = fabs(stepSize_);
429 const double step2 = step * step;
430 if (mscModelCode_ == 0) {// PANDA report PV/01-07 eq(43); linear in step length
431 sigma2 = 225.E-6 * charge_ * charge_ / (beta_ * beta_ * mom_ * mom_) * step / radiationLength_ * matZ_ / (matZ_ + 1) * log(159.*pow(matZ_, -1. / 3.)) / log(287.*pow(matZ_, -0.5)); // sigma^2 = 225E-6*z^2/mom^2 * XX0/beta_^2 * Z/(Z+1) * ln(159*Z^(-1/3))/ln(287*Z^(-1/2)
432
433 } else if (mscModelCode_ == 1) { //Highland not linear in step length formula taken from PDG book 2011 edition
434 double stepOverRadLength = step / radiationLength_;
435 double logCor = (1 + 0.038 * log(stepOverRadLength));
436 sigma2 = 0.0136 * 0.0136 * charge_ * charge_ / (beta_ * beta_ * mom_ * mom_) * stepOverRadLength * logCor * logCor;
437 }
438 //assert(sigma2 >= 0.0);
439 sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0);
440 //XXX std::cout << "MaterialEffects::noiseCoulomb the MSC variance is " << sigma2 << std::endl;
441
442 double noiseAfter[7 * 7]; // will hold the new MSC noise to cause by the current stepSize_ length
443 memset(noiseAfter, 0x00, 7 * 7 * sizeof(double));
444
445 const double *a = direction;
446 // This calculates the MSC angular spread in the 7D global
447 // coordinate system. See PDG 2010, Sec. 27.3 for formulae.
448 noiseAfter[0 * 7 + 0] = sigma2 * step2 / 3.0 * (1 - a[0]*a[0]);
449 noiseAfter[1 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[1];
450 noiseAfter[2 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[2];
451 noiseAfter[3 * 7 + 0] = sigma2 * step * 0.5 * (1 - a[0]*a[0]);
452 noiseAfter[4 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1];
453 noiseAfter[5 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1];
454 noiseAfter[0 * 7 + 1] = noiseAfter[1 * 7 + 0];
455 noiseAfter[1 * 7 + 1] = sigma2 * step2 / 3.0 * (1 - a[1]*a[1]);
456 noiseAfter[2 * 7 + 1] = -sigma2 * step2 / 3.0 * a[1]*a[2];
457 noiseAfter[3 * 7 + 1] = noiseAfter[4 * 7 + 0]; // Cov(x,a_y) = Cov(y,a_x)
458 noiseAfter[4 * 7 + 1] = sigma2 * step * 0.5 * (1 - a[1] * a[1]);
459 noiseAfter[5 * 7 + 1] = -sigma2 * step * 0.5 * a[1]*a[2];
460 noiseAfter[0 * 7 + 2] = noiseAfter[2 * 7 + 0];
461 noiseAfter[1 * 7 + 2] = noiseAfter[2 * 7 + 1];
462 noiseAfter[2 * 7 + 2] = sigma2 * step2 / 3.0 * (1 - a[2]*a[2]);
463 noiseAfter[3 * 7 + 2] = noiseAfter[5 * 7 + 0]; // Cov(z,a_x) = Cov(x,a_z)
464 noiseAfter[4 * 7 + 2] = noiseAfter[5 * 7 + 1]; // Cov(y,a_z) = Cov(z,a_y)
465 noiseAfter[5 * 7 + 2] = sigma2 * step * 0.5 * (1 - a[2]*a[2]);
466 noiseAfter[0 * 7 + 3] = noiseAfter[3 * 7 + 0];
467 noiseAfter[1 * 7 + 3] = noiseAfter[3 * 7 + 1];
468 noiseAfter[2 * 7 + 3] = noiseAfter[3 * 7 + 2];
469 noiseAfter[3 * 7 + 3] = sigma2 * (1 - a[0]*a[0]);
470 noiseAfter[4 * 7 + 3] = -sigma2 * a[0]*a[1];
471 noiseAfter[5 * 7 + 3] = -sigma2 * a[0]*a[2];
472 noiseAfter[0 * 7 + 4] = noiseAfter[4 * 7 + 0];
473 noiseAfter[1 * 7 + 4] = noiseAfter[4 * 7 + 1];
474 noiseAfter[2 * 7 + 4] = noiseAfter[4 * 7 + 2];
475 noiseAfter[3 * 7 + 4] = noiseAfter[4 * 7 + 3];
476 noiseAfter[4 * 7 + 4] = sigma2 * (1 - a[1]*a[1]);
477 noiseAfter[5 * 7 + 4] = -sigma2 * a[1]*a[2];
478 noiseAfter[0 * 7 + 5] = noiseAfter[5 * 7 + 0];
479 noiseAfter[1 * 7 + 5] = noiseAfter[5 * 7 + 1];
480 noiseAfter[2 * 7 + 5] = noiseAfter[5 * 7 + 2];
481 noiseAfter[3 * 7 + 5] = noiseAfter[5 * 7 + 3];
482 noiseAfter[4 * 7 + 5] = noiseAfter[5 * 7 + 4];
483 noiseAfter[5 * 7 + 5] = sigma2 * (1 - a[2]*a[2]);
484// std::cout << "new noise\n";
485// RKTools::printDim(noiseAfter, 7,7);
486 for (unsigned int i = 0; i < 7 * 7; ++i) {
487 noise[i] += noiseAfter[i];
488 }
489}
int i
Definition ShipAna.py:86

◆ noiseCoulomb() [2/2]

void genfit::MaterialEffects::noiseCoulomb ( M7x7 noise,
const M1x3 direction 
) const
private

calculation of multiple scattering

This function first calcuates a MSC variance based on the current material and step length 2 different formulas for the MSC variance are implemeted. One can select the formula via "setMscModel". With the MSC variance and the current direction of the track a full 7D noise matrix is calculated. This noise matrix is the additional noise at the end of fStep in the 7D globa cooridnate system taking even the (co)variances of the position coordinates into account.

◆ setEnergyLossBetheBloch() [1/2]

void genfit::MaterialEffects::setEnergyLossBetheBloch ( bool  opt = true)
inline

Definition at line 72 of file MaterialEffects.h.

◆ setEnergyLossBetheBloch() [2/2]

void genfit::MaterialEffects::setEnergyLossBetheBloch ( bool  opt = true)
inline

Definition at line 72 of file MaterialEffects.h.

◆ setEnergyLossBrems() [1/2]

void genfit::MaterialEffects::setEnergyLossBrems ( bool  opt = true)
inline

Definition at line 75 of file MaterialEffects.h.

◆ setEnergyLossBrems() [2/2]

void genfit::MaterialEffects::setEnergyLossBrems ( bool  opt = true)
inline

Definition at line 75 of file MaterialEffects.h.

◆ setMscModel() [1/2]

void genfit::MaterialEffects::setMscModel ( const std::string &  modelName)

Select the multiple scattering model that will be used during track fit.

At the moment two model are available GEANE and Highland. GEANE is the model was was present in Genfit first. Note that using this function has no effect if setNoiseCoulomb(false) is set.

Definition at line 97 of file MaterialEffects.cc.

98{
99 if (modelName == std::string("GEANE")) {
100 mscModelCode_ = 0;
101 } else if (modelName == std::string("Highland")) {
102 mscModelCode_ = 1;
103 } else {// throw exception
104 std::string errorMsg = std::string("There is no MSC model called \"") + modelName + "\". Maybe it is not implemented or you misspelled the model name";
105 Exception exc(errorMsg, __LINE__, __FILE__);
106 exc.setFatal();
107 std::cerr << exc.what();
108 throw exc;
109 }
110}

◆ setMscModel() [2/2]

void genfit::MaterialEffects::setMscModel ( const std::string &  modelName)

Select the multiple scattering model that will be used during track fit.

At the moment two model are available GEANE and Highland. GEANE is the model was was present in Genfit first. Note that using this function has no effect if setNoiseCoulomb(false) is set.

◆ setNoEffects() [1/2]

void genfit::MaterialEffects::setNoEffects ( bool  opt = true)
inline

Definition at line 70 of file MaterialEffects.h.

70{noEffects_ = opt;}

◆ setNoEffects() [2/2]

void genfit::MaterialEffects::setNoEffects ( bool  opt = true)
inline

Definition at line 70 of file MaterialEffects.h.

70{noEffects_ = opt;}

◆ setNoiseBetheBloch() [1/2]

void genfit::MaterialEffects::setNoiseBetheBloch ( bool  opt = true)
inline

Definition at line 73 of file MaterialEffects.h.

◆ setNoiseBetheBloch() [2/2]

void genfit::MaterialEffects::setNoiseBetheBloch ( bool  opt = true)
inline

Definition at line 73 of file MaterialEffects.h.

◆ setNoiseBrems() [1/2]

void genfit::MaterialEffects::setNoiseBrems ( bool  opt = true)
inline

Definition at line 76 of file MaterialEffects.h.

76{noiseBrems_ = opt; noEffects_ = false;}

◆ setNoiseBrems() [2/2]

void genfit::MaterialEffects::setNoiseBrems ( bool  opt = true)
inline

Definition at line 76 of file MaterialEffects.h.

76{noiseBrems_ = opt; noEffects_ = false;}

◆ setNoiseCoulomb() [1/2]

void genfit::MaterialEffects::setNoiseCoulomb ( bool  opt = true)
inline

Definition at line 74 of file MaterialEffects.h.

74{noiseCoulomb_ = opt; noEffects_ = false;}

◆ setNoiseCoulomb() [2/2]

void genfit::MaterialEffects::setNoiseCoulomb ( bool  opt = true)
inline

Definition at line 74 of file MaterialEffects.h.

74{noiseCoulomb_ = opt; noEffects_ = false;}

◆ stepper() [1/2]

void genfit::MaterialEffects::stepper ( const RKTrackRep rep,
M1x7 state7,
const double &  mom,
double &  relMomLoss,
const int &  pdg,
MaterialProperties currentMaterial,
StepLimits limits,
bool  varField = true 
)

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, or the next material boundary is reached. The material crossed are stored together with their stepsizes.

Definition at line 196 of file MaterialEffects.cc.

204{
205
206 static const double maxRelMomLoss = .005; // maximum relative momentum loss allowed
207 static const double minStep = 1.E-4; // 1 µm
208
209 // Trivial cases
210
211 if (noEffects_)
212 return;
213
214 if (materialInterface_ == nullptr) {
215 std::string msg("MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!");
216 std::runtime_error err(msg);
217 throw err;
218 }
219
220 if (relMomLoss > maxRelMomLoss) {
221 limits.setLimit(stp_momLoss, 0);
222 return;
223 }
224
225
226 double sMax = limits.getLowestLimitSignedVal(); // signed
227
228 if (fabs(sMax) < minStep)
229 return;
230
231
232
233 pdg_ = pdg;
235
236
237 // make minStep
238 state7[0] += limits.getStepSign() * minStep * state7[3];
239 state7[1] += limits.getStepSign() * minStep * state7[4];
240 state7[2] += limits.getStepSign() * minStep * state7[5];
241
242 materialInterface_->initTrack(state7[0], state7[1], state7[2],
243 limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]);
244
246 currentMaterial.setMaterialProperties(matDensity_, matZ_, matA_, radiationLength_, mEE_);
247
248
249 #ifdef DEBUG
250 std::cout << " currentMaterial "; currentMaterial.Print();
251 #endif
252
253 // limit due to momloss
254 double relMomLossPer_cm(0);
255 stepSize_ = 1; // set stepsize for momLoss calculation
256
257 if (matZ_ > 1.E-3) { // don't calculate energy loss for vacuum
258 if (energyLossBetheBloch_) relMomLossPer_cm += this->energyLossBetheBloch() / mom;
259 if (energyLossBrems_) relMomLossPer_cm += this->energyLossBrems() / mom;
260 }
261
262 double maxStepMomLoss = (maxRelMomLoss - relMomLoss) / relMomLossPer_cm;
263 limits.setLimit(stp_momLoss, maxStepMomLoss);
264
265 #ifdef DEBUG
266 std::cout << " momLoss exceeded after a step of " << maxStepMomLoss << "\n";
267 #endif
268
269
270 // now look for boundaries
271 sMax = limits.getLowestLimitSignedVal();
272
273 stepSize_ = limits.getStepSign() * minStep;
274 MaterialProperties materialAfter;
275 M1x3 SA;
276 double boundaryStep(sMax);
277
278 for (unsigned int i=0; i<100; ++i) {
279 #ifdef DEBUG
280 std::cout << " find next boundary\n";
281 #endif
282 double step = materialInterface_->findNextBoundary(rep, state7, boundaryStep, varField);
283 stepSize_ += step;
284 boundaryStep -= step;
285
286 #ifdef DEBUG
287 std::cout << " made a step of " << step << "\n";
288 #endif
289
291 break;
292
293 if (fabs(stepSize_) >= fabs(sMax))
294 break;
295
296 // propagate with found step to boundary
297 rep->RKPropagate(state7, NULL, SA, step, varField);
298
299 // make minStep to cross boundary
300 state7[0] += limits.getStepSign() * minStep * state7[3];
301 state7[1] += limits.getStepSign() * minStep * state7[4];
302 state7[2] += limits.getStepSign() * minStep * state7[5];
303
304 materialInterface_->initTrack(state7[0], state7[1], state7[2],
305 limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]);
306
308
309 #ifdef DEBUG
310 std::cout << " material after step "; materialAfter.Print();
311 #endif
312
313 if (materialAfter != currentMaterial)
314 break;
315 }
316
317 limits.setLimit(stp_boundary, stepSize_);
318
319
320 relMomLoss += relMomLossPer_cm * limits.getLowestLimitVal();
321}
virtual void getMaterialParameters(double &density, double &Z, double &A, double &radiationLength, double &mEE)=0
Get material parameters in current material.
virtual double findNextBoundary(const RKTrackRep *rep, const M1x7 &state7, double sMax, bool varField=true)=0
Make a step until maxStep or the next boundary is reached.
virtual bool initTrack(double posX, double posY, double posZ, double dirX, double dirY, double dirZ)=0
Initialize the navigator at given position and with given direction. Return true if volume changed.
void getParticleParameters(double mom)
sets charge_, mass_ and calculates beta_, gamma_, fgammasquare;
double energyLossBetheBloch()
Returns energy loss.
double energyLossBrems() const
Returns energy loss.
@ stp_boundary
Definition StepLimits.h:44
@ stp_momLoss
Definition StepLimits.h:39
double M1x3[1 *3]
Definition RKTools.h:33

◆ stepper() [2/2]

void genfit::MaterialEffects::stepper ( const RKTrackRep rep,
M1x7 state7,
const double &  mom,
double &  relMomLoss,
const int &  pdg,
MaterialProperties currentMaterial,
StepLimits limits,
bool  varField = true 
)

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, or the next material boundary is reached. The material crossed are stored together with their stepsizes.

Member Data Documentation

◆ beta_

double genfit::MaterialEffects::beta_
private

Definition at line 175 of file MaterialEffects.h.

◆ charge_

int genfit::MaterialEffects::charge_
private

Definition at line 187 of file MaterialEffects.h.

◆ dEdx_

double genfit::MaterialEffects::dEdx_
private

Definition at line 176 of file MaterialEffects.h.

◆ energyLossBetheBloch_

bool genfit::MaterialEffects::energyLossBetheBloch_
private

Definition at line 161 of file MaterialEffects.h.

◆ energyLossBrems_

bool genfit::MaterialEffects::energyLossBrems_
private

Definition at line 164 of file MaterialEffects.h.

◆ gamma_

double genfit::MaterialEffects::gamma_
private

Definition at line 177 of file MaterialEffects.h.

◆ gammaSquare_

double genfit::MaterialEffects::gammaSquare_
private

Definition at line 178 of file MaterialEffects.h.

◆ ignoreBoundariesBetweenEqualMaterials_

bool genfit::MaterialEffects::ignoreBoundariesBetweenEqualMaterials_
private

Definition at line 167 of file MaterialEffects.h.

◆ instance_

MaterialEffects * genfit::MaterialEffects::instance_ = nullptr
staticprivate

Definition at line 58 of file MaterialEffects.h.

◆ mass_

double genfit::MaterialEffects::mass_
private

Definition at line 188 of file MaterialEffects.h.

◆ matA_

double genfit::MaterialEffects::matA_
private

Definition at line 182 of file MaterialEffects.h.

◆ matDensity_

double genfit::MaterialEffects::matDensity_
private

Definition at line 180 of file MaterialEffects.h.

◆ materialInterface_

AbsMaterialInterface * genfit::MaterialEffects::materialInterface_
private

depending on this number a specific msc model is chosen in the noiseCoulomb function.

Definition at line 192 of file MaterialEffects.h.

◆ matZ_

double genfit::MaterialEffects::matZ_
private

Definition at line 181 of file MaterialEffects.h.

◆ me_

const double genfit::MaterialEffects::me_
private

Definition at line 169 of file MaterialEffects.h.

◆ mEE_

double genfit::MaterialEffects::mEE_
private

Definition at line 184 of file MaterialEffects.h.

◆ mom_

double genfit::MaterialEffects::mom_
private

Definition at line 174 of file MaterialEffects.h.

◆ mscModelCode_

int genfit::MaterialEffects::mscModelCode_
private

Definition at line 190 of file MaterialEffects.h.

◆ noEffects_

bool genfit::MaterialEffects::noEffects_
private

Definition at line 159 of file MaterialEffects.h.

◆ noiseBetheBloch_

bool genfit::MaterialEffects::noiseBetheBloch_
private

Definition at line 162 of file MaterialEffects.h.

◆ noiseBrems_

bool genfit::MaterialEffects::noiseBrems_
private

Definition at line 165 of file MaterialEffects.h.

◆ noiseCoulomb_

bool genfit::MaterialEffects::noiseCoulomb_
private

Definition at line 163 of file MaterialEffects.h.

◆ pdg_

int genfit::MaterialEffects::pdg_
private

Definition at line 186 of file MaterialEffects.h.

◆ radiationLength_

double genfit::MaterialEffects::radiationLength_
private

Definition at line 183 of file MaterialEffects.h.

◆ stepSize_

double genfit::MaterialEffects::stepSize_
private

Definition at line 171 of file MaterialEffects.h.


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