#include "GFMaterialEffects.h"

Inheritance diagram for genf::GFMaterialEffects:

Public Member Functions
void	setEnergyLossBetheBloch (bool opt=true)

void	setNoiseBetheBloch (bool opt=true)

void	setNoiseCoulomb (bool opt=true)

void	setEnergyLossBrems (bool opt=true)

void	setNoiseBrems (bool opt=true)

double	effects (const std::vector< TVector3 > &points, const std::vector< double > &pointPaths, const double &mom, const int &pdg, const bool &doNoise=false, TMatrixT< Double_t > noise=NULL, const TMatrixT< Double_t > jacobian=NULL, const TVector3 directionBefore=NULL, const TVector3 directionAfter=NULL)
	Calculates energy loss in the travelled path, optional calculation of noise matrix. More...

double	stepper (const double &maxDist, const double &posx, const double &posy, const double &posz, const double &dirx, const double &diry, const double &dirz, const double &mom, const int &pdg)
	Returns maximum length so that a specified momentum loss will not be exceeded. More...

double	stepper (const double &maxDist, const TVector3 &pos, const TVector3 &dir, const double &mom, const int &pdg)

Static Public Member Functions
static GFMaterialEffects *	getInstance ()

static void	destruct ()

Private Member Functions
	GFMaterialEffects ()

virtual	~GFMaterialEffects ()

void	getParameters ()
	contains energy loss classes More...

void	calcBeta (double mom)
	sets fbeta, fgamma, fgammasquare; must only be used after calling getParameters() More...

double	energyLossBetheBloch (const double &mom)
	Returns energy loss. More...

void	noiseBetheBloch (const double &mom, TMatrixT< double > *noise) const
	calculation of energy loss straggling More...

void	noiseCoulomb (const double &mom, TMatrixT< double > noise, const TMatrixT< double > jacobian, const TVector3 directionBefore, const TVector3 directionAfter) const
	calculation of multiple scattering More...

double	energyLossBrems (const double &mom) const
	Returns energy loss. More...

void	noiseBrems (const double &mom, TMatrixT< double > *noise) const
	calculation of energy loss straggeling More...

double	MeanExcEnergy_get (int Z)

double	MeanExcEnergy_get (TGeoMaterial *)

Private Attributes
bool	fEnergyLossBetheBloch

bool	fNoiseBetheBloch

bool	fNoiseCoulomb

bool	fEnergyLossBrems

bool	fNoiseBrems

const double	me

double	fstep

double	fbeta

double	fdedx

double	fgamma

double	fgammaSquare

double	fmatDensity

double	fmatZ

double	fmatA

double	fradiationLength

double	fmEE

int	fpdg

double	fcharge

double	fmass

Static Private Attributes
static GFMaterialEffects *	finstance = NULL

Detailed Description

Definition at line 58 of file GFMaterialEffects.h.

Constructor & Destructor Documentation

genf::GFMaterialEffects::GFMaterialEffects ( )

private

Definition at line 55 of file GFMaterialEffects.cxx.

Referenced by getInstance().

                                         :
   fEnergyLossBetheBloch(true), fNoiseBetheBloch(true),
   fNoiseCoulomb(true),
   fEnergyLossBrems(true), fNoiseBrems(true),
   me(0.510998910E-3),
   fstep(0),
   fbeta(0),
   fdedx(0),
   fgamma(0),
   fgammaSquare(0),
   fmatDensity(0),
   fmatZ(0),
   fmatA(0),
   fradiationLength(0),
   fmEE(0),
   fpdg(0),
   fcharge(0),
   fmass(0) {
 }

genf::GFMaterialEffects::~GFMaterialEffects ( )

privatevirtual

Definition at line 39 of file GFMaterialEffects.cxx.

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

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().

                                               {
   fbeta = mom/sqrt(fmass*fmass+mom*mom);
 
   //for numerical stability
   fgammaSquare = 1.-fbeta*fbeta;
   if(fgammaSquare>1.E-10) fgammaSquare = 1./fgammaSquare;
   else fgammaSquare = 1.E10;
   fgamma = sqrt(fgammaSquare);
 }

void genf::GFMaterialEffects::destruct ( )

static

Definition at line 80 of file GFMaterialEffects.cxx.

References finstance.

                                      {
   if(finstance != NULL) {
     delete finstance;
     finstance = NULL;
   }
 }

double genf::GFMaterialEffects::effects	(	const std::vector< TVector3 > &	points,
		const std::vector< double > &	pointPaths,
		const double &	mom,
		const int &	pdg,
		const bool &	doNoise = `false`,
		TMatrixT< Double_t > *	noise = `NULL`,
		const TMatrixT< Double_t > *	jacobian = `NULL`,
		const TVector3 *	directionBefore = `NULL`,
		const TVector3 *	directionAfter = `NULL`
	)

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

Definition at line 87 of file GFMaterialEffects.cxx.

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

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

                                                                        {
 
   //assert(points.size()==pointPaths.size());
   fpdg = pdg;
 
   double momLoss=0.;
                      
   for(unsigned int i=1;i<points.size();++i){
     //    TVector3 p1=points.at(i-1);
     //TVector3 p2=points.at(i);
     //TVector3 dir=p2-p1;
     TVector3 dir=points.at(i)-points.at(i-1);
     double dist=dir.Mag();
     double realPath = pointPaths.at(i);
     
     if (dist > 1.E-8) { // do material effects only if distance is not too small
       dir*=1./dist; //normalize dir
 
       double X(0.);
       /*
       double matDensity, matZ, matA, radiationLength, mEE; 
       double step;
       */
       
       gGeoManager->InitTrack(points.at(i-1).X(),points.at(i-1).Y(),points.at(i-1).Z(),
                                dir.X(),dir.Y(),dir.Z());
 
       while(X<dist){
         
         getParameters();
         //        geoMatManager->getMaterialParameters(matDensity, matZ, matA, radiationLength, mEE); 
 
         //        step = geoMatManager->stepOrNextBoundary(dist-X);
         gGeoManager->FindNextBoundaryAndStep(dist-X);
         fstep = gGeoManager->GetStep();
 
         // Loop over EnergyLoss classes
         if(fmatZ>1.E-3){
           /*
           for(unsigned int j=0;j<fEnergyLoss.size();++j){
             momLoss += realPath/dist*fEnergyLoss.at(j)->energyLoss(step,
                                                                    mom,
                                                                    pdg,              
                                                                    matDensity,
                                                                    matZ,
                                                                    matA,
                                                                    radiationLength,
                                                                    mEE,
                                                                    doNoise,
                                                                    noise,
                                                                    jacobian,
                                                                    directionBefore,
                                                                    directionAfter);
                                                                    }
           */
           calcBeta(mom);
 
           if (fEnergyLossBetheBloch)
             momLoss += realPath/dist * this->energyLossBetheBloch(mom);
           if (doNoise && fEnergyLossBetheBloch && fNoiseBetheBloch)
             this->noiseBetheBloch(mom, noise);
 
           if (/*doNoise &&*/ fNoiseCoulomb) // Force it
             this->noiseCoulomb(mom, noise, jacobian, directionBefore, directionAfter);
 
           if (fEnergyLossBrems)
             momLoss += realPath/dist * this->energyLossBrems(mom);
           if (doNoise && fEnergyLossBrems && fNoiseBrems)
             this->noiseBrems(mom, noise);
 
           
 
         }
 
         X += fstep;
 
       }
     }
   }
   //std::cout << "GFMaterialEffects::effects(): momLoss " << momLoss << std::endl;
   return momLoss;
 }

double genf::GFMaterialEffects::energyLossBetheBloch ( const double & mom )

private

Returns energy loss.

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

Definition at line 289 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

                                                                    {
 
   // calc fdedx, also needed in noiseBetheBloch!
   fdedx = 0.307075*fmatZ/fmatA*fmatDensity/(fbeta*fbeta)*fcharge*fcharge;
   double massRatio = me/fmass;
   // me=0.000511 here is in GeV. So fmEE must come in here in eV to get converted to MeV.
   double argument = fgammaSquare*fbeta*fbeta*me*1.E3*2./((1.E-6*fmEE) * sqrt(1+2*sqrt(fgammaSquare)*massRatio + massRatio*massRatio));
   
   if (fmass==0.0) return(0.0);
   if (argument <= exp(fbeta*fbeta))
     { 
       fdedx = 0.;
       // so-called Anderson-Ziegler domain ... Let's approximate it with a flat
       // 100 MeV/cm, looking at the muon dE/dx curve in the PRD Review, or, ahem, wikipedia.
       // http://pdg.lbl.gov/2011/reviews/rpp2011-rev-passage-particles-matter.pdf
       // But there's a practical problem: must not reduce momentum to zero for tracking in G3,
       // so allow only 50% reduction of kinetic energy.
       //      fdedx = 100.0*1.E-3/fstep;  // in GeV/cm, hence 1.e-3
       //if (fdedx > 0.5*0.5*fmass*fbeta*fbeta/fstep) fdedx = 0.5*0.5*fmass*fbeta*fbeta/fstep;
     }
   else{
     fdedx *= (log(argument)-fbeta*fbeta); // Bethe-Bloch [MeV/cm]
     fdedx *= 1.E-3;  // in GeV/cm, hence 1.e-3
     if (fdedx<0.) fdedx = 0;
   }
 
   double DE = fstep * fdedx; //always positive
   double momLoss = sqrt(mom*mom+2.*sqrt(mom*mom+fmass*fmass)*DE+DE*DE) - mom; //always positive
 
   //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
   if(fabs(momLoss)<1.E-11)momLoss=1.E-11;
 
   //  if(fabs(momLoss)>mom)momLoss=mom; // For going fwd and ending its life. EC.
   return momLoss;
 }

double genf::GFMaterialEffects::energyLossBrems ( const double & mom ) const

private

Returns energy loss.

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

Definition at line 483 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

                                                                     {
 
   if (fabs(fpdg)!=11) return 0; // only for electrons and positrons
 
   #if !defined(BETHE)
     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};
     static const double xi=2.51, beta=0.99, vl=0.00004;
   #endif
   #if defined(BETHE) // no MIGDAL corrections
     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};
     static const double xi=2.10, fbeta=1.00, vl=0.001;
   #endif
 
   double BCUT=10000.; // energy up to which soft bremsstrahlung energy loss is calculated
 
   double THIGH=100., CHIGH=50.;
   double dedxBrems=0.;
 
   if(BCUT>0.){
     double T, kc;
 
     if(BCUT>=mom) BCUT=mom; // confine BCUT to mom
 
     // T=mom,  confined to THIGH
     // kc=BCUT, confined to CHIGH ??
     if(mom>=THIGH) {
       T=THIGH;
       if(BCUT>=THIGH) kc=CHIGH;
       else kc=BCUT;
     }
     else {
       T=mom;
       kc=BCUT;
     }
 
     double E=T+me; // total electron energy
     if(BCUT>T) kc=T;
 
     double X=log(T/me);
     double Y=log(kc/(E*vl));
 
     double XX;
     int    K;
     double S=0., YY=1.;
 
     for (unsigned int I=1; I<=2; ++I) {
       XX=1.;
       for (unsigned int J=1; J<=6; ++J) {
         K=6*I+J-6;
         S=S+C[K]*XX*YY;
         XX=XX*X;
       }
       YY=YY*Y;
     }
 
     for (unsigned int I=3; I<=6; ++I) {
        XX=1.;
        for (unsigned int J=1; J<=6; ++J) {
           K=6*I+J-6;
           if(Y<=0.) S=S+C[K]*XX*YY;
           else      S=S+C[K+24]*XX*YY;
           XX=XX*X;
        }
        YY=YY*Y;
     }
 
     double SS=0.;
     YY=1.;
 
     for (unsigned int I=1; I<=2; ++I) {
       XX=1.;
       for (unsigned int J=1; J<=5; ++J) {
         K=5*I+J+55;
         SS=SS+C[K]*XX*YY;
         XX=XX*X;
       }
       YY=YY*Y;
     }
 
     for (unsigned int I=3; I<=5; ++I) {
       XX=1.;
       for (unsigned int J=1; J<=5; ++J) {
         K=5*I+J+55;
         if(Y<=0.) SS=SS+C[K]*XX*YY;
         else      SS=SS+C[K+15]*XX*YY;
         XX=XX*X;
       }
       YY=YY*Y;
     }
 
     S=S+fmatZ*SS;
 
     if(S>0.){
       double CORR=1.;
       #if !defined(BETHE)
         CORR=1./(1.+0.805485E-10*fmatDensity*fmatZ*E*E/(fmatA*kc*kc)); // MIGDAL correction factor
       #endif
 
       double FAC=fmatZ*(fmatZ+xi)*E*E * pow((kc*CORR/T),beta) / (E+me);
       if(FAC<=0.) return 0.;
       dedxBrems=FAC*S;
 
       double RAT;
 
       if(mom>THIGH) {
         if(BCUT<THIGH) {
           RAT=BCUT/mom;
           S=(1.-0.5*RAT+2.*RAT*RAT/9.);
           RAT=BCUT/T;
           S=S/(1.-0.5*RAT+2.*RAT*RAT/9.);
         }
         else {
           RAT=BCUT/mom;
           S=BCUT*(1.-0.5*RAT+2.*RAT*RAT/9.);
           RAT=kc/T;
           S=S/(kc*(1.-0.5*RAT+2.*RAT*RAT/9.));
         }
         dedxBrems=dedxBrems*S; // GeV barn
       }
 
       dedxBrems = 0.60221367*fmatDensity*dedxBrems/fmatA; // energy loss dE/dx [GeV/cm]
     }
   }
 
   if (dedxBrems<0.) dedxBrems = 0;
 
   double factor=1.; // positron correction factor
 
   if (fpdg==-11){
       static const double AA=7522100., A1=0.415, A3=0.0021, A5=0.00054;
 
       double ETA=0.;
       if(fmatZ>0.) {
         double X=log(AA*mom/fmatZ*fmatZ);
         if(X>-8.) {
           if(X>=+9.) ETA=1.;
           else {
              double W=A1*X+A3*pow(X,3.)+A5*pow(X,5.);
              ETA=0.5+atan(W)/M_PI;
           }
         }
       }
 
       double E0;
 
       if(ETA<0.0001) factor=1.E-10;
       else if (ETA>0.9999) factor=1.;
       else {
          E0=BCUT/mom;
          if(E0>1.) E0=1.;
          if(E0<1.E-8) factor=1.;
          else factor = ETA * ( 1.-pow(1.-E0, 1./ETA) ) / E0;
       }
   }
 
   double DE = fstep * factor*dedxBrems; //always positive
   double momLoss = sqrt(mom*mom+2.*sqrt(mom*mom+fmass*fmass)*DE+DE*DE) - mom; //always positive
 
   return momLoss;
 }

genf::GFMaterialEffects * genf::GFMaterialEffects::getInstance ( )

static

Definition at line 75 of file GFMaterialEffects.cxx.

References finstance, and GFMaterialEffects().

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

                                                           {
   if(finstance == NULL) finstance = new GFMaterialEffects();
   return finstance;
 }

void genf::GFMaterialEffects::getParameters ( )

private

contains energy loss classes

interface to material and geometry

Definition at line 254 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

                                          {
   if (!gGeoManager->GetCurrentVolume()->GetMedium())
     throw GFException(std::string(__func__) + ": no medium", __LINE__, __FILE__).setFatal();
   TGeoMaterial * mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
   fmatDensity      = mat->GetDensity();
   fmatZ            = mat->GetZ();
   fmatA            = mat->GetA();
   fradiationLength = mat->GetRadLen();
   fmEE             = MeanExcEnergy_get(mat);
 
   // You know what? F*ck it. Just force this to be LAr.... is what I *could/will* say here ....
   // See comment in energyLossBetheBloch() for why fmEE is in eV here.
   fmatDensity = 1.40; fmatZ = 18.0; fmatA = 39.95; fradiationLength=13.947; fmEE=188.0;
   
 
   TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fpdg);
   fcharge = part->Charge()/(3.);
   fmass = part->Mass();
 }

double genf::GFMaterialEffects::MeanExcEnergy_get ( int Z )

private

Definition at line 669 of file GFMaterialEffects.cxx.

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

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

                                                     {
   if ((Z < 0) || (Z > MeanExcEnergy_NELEMENTS))
     throw GFException(std::string(__func__) + ": unsupported Z", __LINE__, __FILE__).setFatal();
   return MeanExcEnergy_vals[Z];
 }

double genf::GFMaterialEffects::MeanExcEnergy_get ( TGeoMaterial * mat )

private

Definition at line 675 of file GFMaterialEffects.cxx.

References MeanExcEnergy_get().

                                                                 {
   if(mat->IsMixture()){
     double logMEE = 0.;
     double denom  = 0.;
     TGeoMixture *mix = (TGeoMixture*)mat;
     for(int i=0;i<mix->GetNelements();++i){
       int index = int(floor((mix->GetZmixt())[i]));
       logMEE += 1./(mix->GetAmixt())[i]*(mix->GetWmixt())[i]*(mix->GetZmixt())[i]*log(MeanExcEnergy_get(index));
       denom  += (mix->GetWmixt())[i]*(mix->GetZmixt())[i]*1./(mix->GetAmixt())[i];
     }
     logMEE/=denom;
     return exp(logMEE);
   }
   else{ // not a mixture
     int index = int(floor(mat->GetZ()));
     return MeanExcEnergy_get(index);
   }
 }

void genf::GFMaterialEffects::noiseBetheBloch	(	const double &	mom,
		TMatrixT< double > *	noise
	)		const

private

calculation of energy loss straggling

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

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

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

Definition at line 326 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

                                                                       {
 
 
   // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E);
   double sigma2E = 0.;
   double zeta  = 153.4E3 * fcharge*fcharge/(fbeta*fbeta) * fmatZ/fmatA * fmatDensity * fstep; // eV
   double Emax  = 2.E9*me*fbeta*fbeta*fgammaSquare / (1. + 2.*fgamma*me/fmass + (me/fmass)*(me/fmass) ); // eV
   double kappa = zeta/Emax;
 
   if (kappa > 0.01) { // Vavilov-Gaussian regime
     sigma2E += zeta*Emax*(1.-fbeta*fbeta/2.);  // eV^2
   }
   else { // Urban/Landau approximation
     double alpha = 0.996;
     double sigmaalpha = 15.76;
     // calculate number of collisions Nc
     double I = 16. * pow(fmatZ, 0.9); // eV
     double f2 = 0.;
     if (fmatZ > 2.) f2 = 2./fmatZ;
     double f1 = 1. - f2;
     double e2 = 10.*fmatZ*fmatZ; // eV
     double e1 = pow( (I/pow(e2,f2)), 1./f1);  // eV
 
     double mbbgg2 = 2.E9*fmass*fbeta*fbeta*fgammaSquare; // eV
     double Sigma1 = fdedx*1.0E9 * f1/e1 * (log(mbbgg2 / e1) - fbeta*fbeta) / (log(mbbgg2 / I) - fbeta*fbeta) * 0.6; // 1/cm
     double Sigma2 = fdedx*1.0E9 * f2/e2 * (log(mbbgg2 / e2) - fbeta*fbeta) / (log(mbbgg2 / I) - fbeta*fbeta) * 0.6; // 1/cm
     double Sigma3 = fdedx*1.0E9 * Emax / ( I*(Emax+I)*log((Emax+I)/I) ) * 0.4; // 1/cm
 
     double Nc = (Sigma1 + Sigma2 + Sigma3)*fstep;
 
     if (Nc > 50.) { // truncated Landau distribution
       // calculate sigmaalpha  (see GEANT3 manual W5013)
       double RLAMED = -0.422784 - fbeta*fbeta - log(zeta/Emax);
       double RLAMAX =  0.60715 + 1.1934*RLAMED +(0.67794 + 0.052382*RLAMED)*exp(0.94753+0.74442*RLAMED);
       // from lambda max to sigmaalpha=sigma (empirical polynomial)
       if(RLAMAX <= 1010.) {
          sigmaalpha =  1.975560
                       +9.898841e-02 *RLAMAX
                       -2.828670e-04 *RLAMAX*RLAMAX
                       +5.345406e-07 *pow(RLAMAX,3.)
                       -4.942035e-10 *pow(RLAMAX,4.)
                       +1.729807e-13 *pow(RLAMAX,5.);
       }
       else { sigmaalpha = 1.871887E+01 + 1.296254E-02 *RLAMAX; }
       // alpha=54.6  corresponds to a 0.9996 maximum cut
       if(sigmaalpha > 54.6) sigmaalpha=54.6;
       sigma2E += sigmaalpha*sigmaalpha * zeta*zeta;  // eV^2
     }
     else { // Urban model
       double Ealpha  = I / (1.-(alpha*Emax/(Emax+I)));   // eV
       double meanE32 = I*(Emax+I)/Emax * (Ealpha-I);     // eV^2
       sigma2E += fstep * (Sigma1*e1*e1 + Sigma2*e2*e2 + Sigma3*meanE32); // eV^2
     }
   }
 
   sigma2E*=1.E-18; // eV -> GeV
 
   // update noise matrix
   (*noise)[6][6] += (mom*mom+fmass*fmass)/pow(mom,6.)*sigma2E;
 }

void genf::GFMaterialEffects::noiseBrems	(	const double &	mom,
		TMatrixT< double > *	noise
	)		const

private

calculation of energy loss straggeling

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

Definition at line 645 of file GFMaterialEffects.cxx.

References fmass, fpdg, fradiationLength, and fstep.

Referenced by effects(), and stepper().

                                                                  {
 
   if (fabs(fpdg)!=11) return; // only for electrons and positrons
 
   double LX  = 1.442695*fstep/fradiationLength;
   double S2B = mom*mom * ( 1./pow(3.,LX) - 1./pow(4.,LX) );
   double DEDXB  = pow(fabs(S2B),0.5);
   DEDXB = 1.2E9*DEDXB; //eV
   double sigma2E = DEDXB*DEDXB; //eV^2
   sigma2E*=1.E-18; // eV -> GeV
 
   (*noise)[6][6] += (mom*mom+fmass*fmass)/pow(mom,6.)*sigma2E;
 }

void genf::GFMaterialEffects::noiseCoulomb	(	const double &	mom,
		TMatrixT< double > *	noise,
		const TMatrixT< double > *	jacobian,
		const TVector3 *	directionBefore,
		const TVector3 *	directionAfter
	)		const

private

calculation of multiple scattering

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

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

This is a detailed description of the mathematics involved:
We define a local coordinate system cs' with initial momentum in z-direction:
$\left(\begin{array}{c} x'\\ y'\\ z'\\ a_{x}'\\ a_{y}'\\ a_{z}'\\ \frac{q}{p}'\end{array}\right)=\left(\begin{array}{ccccccc} \cos\psi & \sin\psi & 0\\ -\cos\vartheta\sin\psi & \cos\vartheta\cos\psi & \sin\psi\\ \sin\vartheta\sin\psi & -\sin\vartheta\cos\psi & \cos\vartheta\\ & & & \cos\psi & \sin\psi & 0\\ & & & -\cos\vartheta\sin\psi & \cos\vartheta\cos\psi & \sin\psi\\ & & & \sin\vartheta\sin\psi & -\sin\vartheta\cos\psi & \cos\vartheta\\ & & & & & & 1\end{array}\right)\left(\begin{array}{c} x\\ y\\ z\\ a_{x}\\ a_{y}\\ a_{z}\\ \frac{q}{p}\end{array}\right)=R^{T}\left(\begin{array}{c} x\\ y\\ z\\ a_{x}\\ a_{y}\\ a_{z}\\ \frac{q}{p}\end{array}\right)$

Now the global coordinate system cs is:

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

with the Euler angles

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Definition at line 389 of file GFMaterialEffects.cxx.

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

Referenced by effects(), and stepper().

                                                                           {
 
   // MULTIPLE SCATTERING; calculate sigma^2
   // PANDA report PV/01-07 eq(43); linear in step length
   double sigma2 = 225.E-6/(fbeta*fbeta*mom*mom) * fstep/fradiationLength * fmatZ/(fmatZ+1) * log(159.*pow(fmatZ,-1./3.))/log(287.*pow(fmatZ,-0.5)); // sigma^2 = 225E-6/mom^2 * XX0/fbeta^2 * Z/(Z+1) * ln(159*Z^(-1/3))/ln(287*Z^(-1/2)
 
   // noiseBefore
     TMatrixT<double> noiseBefore(7,7);
 
     // calculate euler angles theta, psi (so that directionBefore' points in z' direction)
     double psi = 0;
     if (fabs((*directionBefore)[1]) < 1E-14) {  // numerical case: arctan(+-inf)=+-PI/2
       if ((*directionBefore)[0] >= 0.) psi = M_PI/2.;
       else psi = 3.*M_PI/2.;
     }
     else {
       if ((*directionBefore)[1] > 0.) psi = M_PI - atan((*directionBefore)[0]/(*directionBefore)[1]);
       else psi = -atan((*directionBefore)[0]/(*directionBefore)[1]);
     }
     // cache sin and cos
     double sintheta = sqrt(1-(*directionBefore)[2]*(*directionBefore)[2]);  // theta = arccos(directionBefore[2])
     double costheta = (*directionBefore)[2];
     double sinpsi = sin(psi);
     double cospsi = cos(psi);
 
     // calculate NoiseBefore Matrix R M R^T
     double noiseBefore34 =  sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseBefore_ij = noiseBefore_ji
     double noiseBefore35 = -sigma2 * costheta * sinpsi * sintheta;
     double noiseBefore45 =  sigma2 * costheta * cospsi * sintheta;
 
     noiseBefore[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
     noiseBefore[4][3] = noiseBefore34;
     noiseBefore[5][3] = noiseBefore35;
 
     noiseBefore[3][4] = noiseBefore34;
     noiseBefore[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
     noiseBefore[5][4] = noiseBefore45;
 
     noiseBefore[3][5] = noiseBefore35;
     noiseBefore[4][5] = noiseBefore45;
     noiseBefore[5][5] = sigma2 * sintheta*sintheta;
 
     TMatrixT<double> jacobianT(7,7);
     jacobianT = (*jacobian);
     jacobianT.T();
 
     noiseBefore = jacobianT*noiseBefore*(*jacobian); //propagate
 
   // noiseAfter
     TMatrixT<double> noiseAfter(7,7);
 
     // calculate euler angles theta, psi (so that A' points in z' direction)
     psi = 0;
     if (fabs((*directionAfter)[1]) < 1E-14) {  // numerical case: arctan(+-inf)=+-PI/2
       if ((*directionAfter)[0] >= 0.) psi = M_PI/2.;
       else psi = 3.*M_PI/2.;
     }
     else {
       if ((*directionAfter)[1] > 0.) psi = M_PI - atan((*directionAfter)[0]/(*directionAfter)[1]);
       else psi = -atan((*directionAfter)[0]/(*directionAfter)[1]);
     }
     // cache sin and cos
     sintheta = sqrt(1-(*directionAfter)[2]*(*directionAfter)[2]);  // theta = arccos(directionAfter[2])
     costheta = (*directionAfter)[2];
     sinpsi = sin(psi);
     cospsi = cos(psi);
 
     // calculate NoiseAfter Matrix R M R^T
     double noiseAfter34 =  sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseAfter_ij = noiseAfter_ji
     double noiseAfter35 = -sigma2 * costheta * sinpsi * sintheta;
     double noiseAfter45 =  sigma2 * costheta * cospsi * sintheta;
 
     noiseAfter[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
     noiseAfter[4][3] = noiseAfter34;
     noiseAfter[5][3] = noiseAfter35;
 
     noiseAfter[3][4] = noiseAfter34;
     noiseAfter[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
     noiseAfter[5][4] = noiseAfter45;
 
     noiseAfter[3][5] = noiseAfter35;
     noiseAfter[4][5] = noiseAfter45;
     noiseAfter[5][5] = sigma2 * sintheta*sintheta;
 
   //calculate mean of noiseBefore and noiseAfter and update noise
     (*noise) += 0.5*noiseBefore + 0.5*noiseAfter;
 
 }

void genf::GFMaterialEffects::setEnergyLossBetheBloch ( bool opt = true )

inline

Definition at line 69 of file GFMaterialEffects.h.

References fEnergyLossBetheBloch.

69 {fEnergyLossBetheBloch=opt;}

genf::GFMaterialEffects::fEnergyLossBetheBloch

bool fEnergyLossBetheBloch

Definition: GFMaterialEffects.h:347

void genf::GFMaterialEffects::setEnergyLossBrems ( bool opt = true )

inline

Definition at line 72 of file GFMaterialEffects.h.

References fEnergyLossBrems.

72 {fEnergyLossBrems=opt;}

genf::GFMaterialEffects::fEnergyLossBrems

bool fEnergyLossBrems

Definition: GFMaterialEffects.h:350

void genf::GFMaterialEffects::setNoiseBetheBloch ( bool opt = true )

inline

Definition at line 70 of file GFMaterialEffects.h.

References fNoiseBetheBloch.

70 {fNoiseBetheBloch=opt;}

genf::GFMaterialEffects::fNoiseBetheBloch

bool fNoiseBetheBloch

Definition: GFMaterialEffects.h:348

void genf::GFMaterialEffects::setNoiseBrems ( bool opt = true )

inline

Definition at line 73 of file GFMaterialEffects.h.

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

73 {fNoiseBrems=opt;}

genf::GFMaterialEffects::fNoiseBrems

bool fNoiseBrems

Definition: GFMaterialEffects.h:351

void genf::GFMaterialEffects::setNoiseCoulomb ( bool opt = true )

inline

Definition at line 71 of file GFMaterialEffects.h.

References fNoiseCoulomb.

71 {fNoiseCoulomb=opt;}

genf::GFMaterialEffects::fNoiseCoulomb

bool fNoiseCoulomb

Definition: GFMaterialEffects.h:349

double genf::GFMaterialEffects::stepper	(	const double &	maxDist,
		const double &	posx,
		const double &	posy,
		const double &	posz,
		const double &	dirx,
		const double &	diry,
		const double &	dirz,
		const double &	mom,
		const int &	pdg
	)

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

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

Definition at line 182 of file GFMaterialEffects.cxx.

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

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

                                                {
 
   static const double maxPloss = .005; // maximum relative momentum loss allowed
 
   gGeoManager->InitTrack(posx,posy,posz,dirx,diry,dirz);
 
   double X(0.);
   double dP = 0.;
   double momLoss = 0.;
   /*
   double matDensity, matZ, matA, radiationLength, mEE;
   double step;
   */
 
   while(X<maxDist){
     
     getParameters();
 
     gGeoManager->FindNextBoundaryAndStep(maxDist-X);
     fstep = gGeoManager->GetStep();
     //
     //    step = geoMatManager->stepOrNextBoundary(maxDist-X);
     
     // Loop over EnergyLoss classes
     
     if(fmatZ>1.E-3){
       /*  
       for(unsigned int j=0;j<fEnergyLoss.size();++j){
         momLoss += fEnergyLoss.at(j)->energyLoss(step,
                                                  mom,
                                                  pdg,              
                                                  matDensity,
                                                  matZ,
                                                  matA,
                                                  radiationLength,
                                                  mEE);  
                                                  }
     */
       calcBeta(mom);
  
       if (fEnergyLossBetheBloch)
         momLoss += this->energyLossBetheBloch(mom);
 
       if (fEnergyLossBrems)
         momLoss += this->energyLossBrems(mom);
       
     }
     
     if(dP + momLoss > mom*maxPloss){
       double fraction = (mom*maxPloss-dP)/momLoss;
       if ((fraction <= 0.) || (fraction >= 1.))
         throw GFException(std::string(__func__) + ": invalid fraction", __LINE__, __FILE__).setFatal();
       dP+=fraction*momLoss;
       X+=fraction*fstep;
       break;
     }
     
     dP += momLoss;
     X += fstep;
   }
 
   return X;
 }

double genf::GFMaterialEffects::stepper	(	const double &	maxDist,
		const TVector3 &	pos,
		const TVector3 &	dir,
		const double &	mom,
		const int &	pdg
	)

inline

Definition at line 98 of file GFMaterialEffects.h.

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

102 {

103 return stepper(maxDist, pos.X(),pos.Y(),pos.Z(),dir.X(),dir.Y(),dir.Z(),mom,pdg);};

dir

TDirectory * dir

Definition: macro.C:5

genf::GFMaterialEffects::stepper

double stepper(const double &maxDist, const double &posx, const double &posy, const double &posz, const double &dirx, const double &diry, const double &dirz, const double &mom, const int &pdg)

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

Definition: GFMaterialEffects.cxx:182