#include <GFEnergyLossCoulomb.h>

Inheritance diagram for genf::GFEnergyLossCoulomb:

Public Member Functions
double	energyLoss (const double &step, const double &mom, const int &pdg, const double &matDensity, const double &matZ, const double &matA, const double &radiationLength, const double &meanExcitationEnergy, const bool &doNoise=false, TMatrixT< Double_t > noise=NULL, const TMatrixT< Double_t > jacobian=NULL, const TVector3 directionBefore=NULL, const TVector3 directionAfter=NULL)
	Optional calculation of multiple scattering. More...

virtual	~GFEnergyLossCoulomb ()

Public Member Functions inherited from genf::GFAbsEnergyLoss
virtual	~GFAbsEnergyLoss ()

Additional Inherited Members
Protected Member Functions inherited from genf::GFAbsEnergyLoss
void	getParticleParameters (const int &pdg, double &charge, double &mass)
	Gets particle charge and mass (in GeV) More...

double	getParticleMass (const int &pdg)
	Returns particle mass (in GeV) More...

Detailed Description

Definition at line 43 of file GFEnergyLossCoulomb.h.

Constructor & Destructor Documentation

genf::GFEnergyLossCoulomb::~GFEnergyLossCoulomb ( )

virtual

Definition at line 25 of file GFEnergyLossCoulomb.cxx.

25 {

26 }

Member Function Documentation

double genf::GFEnergyLossCoulomb::energyLoss	(	const double &	step,
		const double &	mom,
		const int &	pdg,
		const double &	matDensity,
		const double &	matZ,
		const double &	matA,
		const double &	radiationLength,
		const double &	meanExcitationEnergy,
		const bool &	doNoise = `false`,
		TMatrixT< Double_t > *	noise = `NULL`,
		const TMatrixT< Double_t > *	jacobian = `NULL`,
		const TVector3 *	directionBefore = `NULL`,
		const TVector3 *	directionAfter = `NULL`
	)

virtual

Optional calculation of multiple scattering.

Only multiple scattering is calculated, the returned value for the energy loss is always 0. 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

Implements genf::GFAbsEnergyLoss.

Definition at line 28 of file GFEnergyLossCoulomb.cxx.

                                                                       {
 
   double charge, mass;
   getParticleParameters(pdg, charge, mass);
 
   const double beta = mom/sqrt(mass*mass+mom*mom);
 
   if (doNoise) {
     if (!noise) throw GFException(std::string(__func__) + ": no noise given!", __LINE__, __FILE__).setFatal();
     if (!jacobian) throw GFException(std::string(__func__) + ": no jacobian given!", __LINE__, __FILE__).setFatal();
     if (!directionBefore) throw GFException(std::string(__func__) + ": no directionBefore given!", __LINE__, __FILE__).setFatal();
     if (!directionAfter) throw GFException(std::string(__func__) + ": no directionAfter given!", __LINE__, __FILE__).setFatal();
 
     // MULTIPLE SCATTERING; calculate sigma^2
     // PANDA report PV/01-07 eq(43); linear in step length
     double sigma2 = 225.E-6/(beta*beta*mom*mom) * step/radiationLength * matZ/(matZ+1) * log(159.*pow(matZ,-1./3.))/log(287./sqrt(matZ)); // sigma^2 = 225E-6/mom^2 * XX0/beta^2 * Z/(Z+1) * ln(159*Z^(-1/3))/ln(287*Z^(-1/2)
 
 
     // noiseBefore
       TMatrixT<Double_t> 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_t> jacobianT(7,7);
       jacobianT = (*jacobian);
       jacobianT.T();
 
       noiseBefore = jacobianT*noiseBefore*(*jacobian); //propagate
 
     // noiseAfter
       TMatrixT<Double_t> 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;
   }
 
   return 0.;
 }

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

larreco/larreco/Genfit/GFEnergyLossCoulomb.h
larreco/larreco/Genfit/GFEnergyLossCoulomb.cxx

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation