Belle II Software  release-08-01-10
Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
MaterialEffects Class Reference

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

#include <MaterialEffects.h>

Collaboration diagram for MaterialEffects:
Collaboration graph
[legend]

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 setMagCharge (double magCharge)
 
void setMscModel (const std::string &modelName)
 Select the multiple scattering model that will be used during track fit. More...
 
double effects (const std::vector< 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, Material &currentMaterial, StepLimits &limits, bool varField=true)
 Returns maximum length so that a specified momentum loss will not be exceeded. More...
 
void setDebugLvl (unsigned int lvl=1)
 
void drawdEdx (int pdg=11)
 

Static Public Member Functions

static MaterialEffects * getInstance ()
 
static void destruct ()
 

Private Member Functions

void getParticleParameters ()
 sets charge_, mass_
 
void getMomGammaBeta (double Energy, double &mom, double &gammaSquare, double &gamma, double &betaSquare) const
 
double momentumLoss (double stepSign, double mom, bool linear)
 Returns momentum loss. More...
 
double dEdx (double Energy)
 Calculate dEdx for a given energy.
 
double dEdxBetheBloch (double betaSquare, double gamma, double gammasquare) const
 Uses Bethe Bloch formula to calculate dEdx.
 
void noiseBetheBloch (M7x7 &noise, double mom, double betaSquare, double gamma, double gammaSquare) const
 calculation of energy loss straggeling More...
 
void noiseCoulomb (M7x7 &noise, const M1x3 &direction, double momSquare, double betaSquare) const
 calculation of multiple scattering More...
 
double dEdxBrems (double mom) const
 Returns dEdx. More...
 
void noiseBrems (M7x7 &noise, double momSquare, double betaSquare) const
 calculation of energy loss straggeling More...
 

Private Attributes

bool noEffects_
 
bool energyLossBetheBloch_
 
bool noiseBetheBloch_
 
bool noiseCoulomb_
 
bool energyLossBrems_
 
bool noiseBrems_
 
bool ignoreBoundariesBetweenEqualMaterials_
 
const double me_
 
double stepSize_
 
double dEdx_
 
double E_
 
double matDensity_
 
double matZ_
 
double matA_
 
double radiationLength_
 
double mEE_
 
int pdg_
 
double charge_
 
double mag_charge_
 
double mass_
 
int mscModelCode_
 
AbsMaterialInterface * materialInterface_
 depending on this number a specific msc model is chosen in the noiseCoulomb function.
 
unsigned int debugLvl_
 

Static Private Attributes

static MaterialEffects * instance_ = 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 50 of file MaterialEffects.h.

Member Function Documentation

◆ dEdxBrems()

double dEdxBrems ( double  mom) const
private

Returns dEdx.

Can be called with any pdg, but only calculates dEdx for electrons and positrons (otherwise returns 0). Uses a gaussian approximation (Bethe-Heitler formula with Migdal corrections). For positrons, dEdx is weighed with a correction factor.

Definition at line 628 of file MaterialEffects.cc.

629 {
630 
631  // Code ported from GEANT 3 (gbrele.F)
632 
633  if (abs(pdg_) != 11) return 0; // only for electrons and positrons
634 
635 #if !defined(BETHE)
636  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};
637  static const double xi = 2.51, beta = 0.99, vl = 0.00004;
638 #endif
639 #if defined(BETHE) // no MIGDAL corrections
640  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};
641  static const double xi = 2.10, beta = 1.00, vl = 0.001;
642 #endif
643 
644  double BCUT = 10000.; // energy up to which soft bremsstrahlung energy loss is calculated
645 
646  static const double THIGH = 100., CHIGH = 50.;
647  double dedxBrems = 0.;
648 
649  if (BCUT > 0.) {
650  double T, kc;
651 
652  if (BCUT > mom) BCUT = mom; // confine BCUT to mom_
653 
654  // T=mom_, confined to THIGH
655  // kc=BCUT, confined to CHIGH ??
656  if (mom > THIGH) {
657  T = THIGH;
658  if (BCUT >= THIGH)
659  kc = CHIGH;
660  else
661  kc = BCUT;
662  } else {
663  T = mom;
664  kc = BCUT;
665  }
666 
667  double E = T + me_; // total electron energy
668  if (BCUT > T)
669  kc = T;
670 
671  double X = log(T / me_);
672  double Y = log(kc / (E * vl));
673 
674  double XX;
675  int K;
676  double S = 0., YY = 1.;
677 
678  for (unsigned int I = 1; I <= 2; ++I) {
679  XX = 1.;
680  for (unsigned int J = 1; J <= 6; ++J) {
681  K = 6 * I + J - 6;
682  S += C[K] * XX * YY;
683  XX *= X;
684  }
685  YY *= Y;
686  }
687 
688  for (unsigned int I = 3; I <= 6; ++I) {
689  XX = 1.;
690  for (unsigned int J = 1; J <= 6; ++J) {
691  K = 6 * I + J - 6;
692  if (Y > 0.)
693  K += 24;
694  S += C[K] * XX * YY;
695  XX *= X;
696  }
697  YY *= Y;
698  }
699 
700  double SS = 0.;
701  YY = 1.;
702 
703  for (unsigned int I = 1; I <= 2; ++I) {
704  XX = 1.;
705  for (unsigned int J = 1; J <= 5; ++J) {
706  K = 5 * I + J + 55;
707  SS += C[K] * XX * YY;
708  XX *= X;
709  }
710  YY *= Y;
711  }
712 
713  for (unsigned int I = 3; I <= 5; ++I) {
714  XX = 1.;
715  for (unsigned int J = 1; J <= 5; ++J) {
716  K = 5 * I + J + 55;
717  if (Y > 0.)
718  K += 15;
719  SS += C[K] * XX * YY;
720  XX *= X;
721  }
722  YY *= Y;
723  }
724 
725  S += matZ_ * SS;
726 
727  if (S > 0.) {
728  double CORR = 1.;
729 #if !defined(BETHE)
730  CORR = 1. / (1. + 0.805485E-10 * matDensity_ * matZ_ * E * E / (matA_ * kc * kc)); // MIGDAL correction factor
731 #endif
732 
733  // We use exp(beta * log(...) here because pow(..., beta) is
734  // REALLY slow and we don't need ultimate numerical precision
735  // for this approximation.
736  double FAC = matZ_ * (matZ_ + xi) * E * E / (E + me_);
737  if (beta == 1.) // That is the #ifdef BETHE case
738  FAC *= kc * CORR / T;
739  else
740  FAC *= exp(beta * log(kc * CORR / T));
741  if (FAC <= 0.)
742  return 0.;
743  dedxBrems = FAC * S;
744 
745 
746  if (mom >= THIGH) {
747  double RAT;
748  if (BCUT < THIGH) {
749  RAT = BCUT / mom;
750  S = (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
751  RAT = BCUT / T;
752  S /= 1. - 0.5 * RAT + 2.*RAT * RAT / 9.;
753  } else {
754  RAT = BCUT / mom;
755  S = BCUT * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
756  RAT = kc / T;
757  S /= kc * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
758  }
759  dedxBrems *= S; // GeV barn
760  }
761 
762  dedxBrems *= 0.60221367 * matDensity_ / matA_; // energy loss dE/dx [GeV/cm]
763  }
764  }
765 
766  if (dedxBrems < 0.)
767  dedxBrems = 0;
768 
769  double factor = 1.; // positron correction factor
770 
771  if (pdg_ == -11) {
772  static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054;
773 
774  double ETA = 0.;
775  if (matZ_ > 0.) {
776  double X = log(AA * mom / (matZ_ * matZ_));
777  if (X > -8.) {
778  if (X >= +9.) ETA = 1.;
779  else {
780  double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.);
781  ETA = 0.5 + atan(W) / M_PI;
782  }
783  }
784  }
785 
786  if (ETA < 0.0001)
787  factor = 1.E-10;
788  else if (ETA > 0.9999)
789  factor = 1.;
790  else {
791  double E0 = BCUT / mom;
792  if (E0 > 1.)
793  E0 = 1.;
794  if (E0 < 1.E-8)
795  factor = 1.;
796  else
797  factor = ETA * (1. - pow(1. - E0, 1. / ETA)) / E0;
798  }
799  }
800 
801  return factor * dedxBrems; //always positive
802 }
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36
K
#define K(x)
macro autogenerated by FFTW
Definition: DiscreteCosineTransform_31points.cc:37
Belle2::atan
double atan(double a)
atan for double
Definition: beamHelpers.h:34

◆ momentumLoss()

double momentumLoss ( double  stepSign,
double  mom,
bool  linear 
)
private

Returns momentum loss.

Also sets dEdx_ and E_.

Definition at line 392 of file MaterialEffects.cc.

◆ noiseBetheBloch()

void noiseBetheBloch ( M7x7 &  noise,
double  mom,
double  betaSquare,
double  gamma,
double  gammaSquare 
) 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 momentumLoss, so it has to be called afterwards!

Definition at line 497 of file MaterialEffects.cc.

◆ noiseBrems()

void noiseBrems ( M7x7 &  noise,
double  momSquare,
double  betaSquare 
) const
private

calculation of energy loss straggeling

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

Definition at line 805 of file MaterialEffects.cc.

◆ noiseCoulomb()

void noiseCoulomb ( M7x7 &  noise,
const M1x3 &  direction,
double  momSquare,
double  betaSquare 
) 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 557 of file MaterialEffects.cc.

◆ setMscModel()

void 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 98 of file MaterialEffects.cc.

◆ stepper()

void stepper ( const RKTrackRep *  rep,
M1x7 &  state7,
const double &  mom,
double &  relMomLoss,
const int &  pdg,
Material &  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 219 of file MaterialEffects.cc.

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