LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Protected Member Functions | List of all members
genf::GFEnergyLossCoulomb Class Reference

#include "GFEnergyLossCoulomb.h"

Inheritance diagram for genf::GFEnergyLossCoulomb:
genf::GFAbsEnergyLoss

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

Protected Member Functions

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.

References beta, E, genf::GFAbsEnergyLoss::getParticleParameters(), and GFException::setFatal().

40  {
41 
42  double charge, mass;
43  getParticleParameters(pdg, charge, mass);
44 
45  const double beta = mom/sqrt(mass*mass+mom*mom);
46 
47  if (doNoise) {
48  if (!noise) throw GFException(std::string(__func__) + ": no noise given!", __LINE__, __FILE__).setFatal();
49  if (!jacobian) throw GFException(std::string(__func__) + ": no jacobian given!", __LINE__, __FILE__).setFatal();
50  if (!directionBefore) throw GFException(std::string(__func__) + ": no directionBefore given!", __LINE__, __FILE__).setFatal();
51  if (!directionAfter) throw GFException(std::string(__func__) + ": no directionAfter given!", __LINE__, __FILE__).setFatal();
52 
53  // MULTIPLE SCATTERING; calculate sigma^2
54  // PANDA report PV/01-07 eq(43); linear in step length
55  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)
56 
57 
58  // noiseBefore
59  TMatrixT<Double_t> noiseBefore(7,7);
60 
61  // calculate euler angles theta, psi (so that directionBefore' points in z' direction)
62  double psi = 0;
63  if (fabs((*directionBefore)[1]) < 1E-14) { // numerical case: arctan(+-inf)=+-PI/2
64  if ((*directionBefore)[0] >= 0.) psi = M_PI/2.;
65  else psi = 3.*M_PI/2.;
66  }
67  else {
68  if ((*directionBefore)[1] > 0.) psi = M_PI - atan((*directionBefore)[0]/(*directionBefore)[1]);
69  else psi = -atan((*directionBefore)[0]/(*directionBefore)[1]);
70  }
71  // cache sin and cos
72  double sintheta = sqrt(1-(*directionBefore)[2]*(*directionBefore)[2]); // theta = arccos(directionBefore[2])
73  double costheta = (*directionBefore)[2];
74  double sinpsi = sin(psi);
75  double cospsi = cos(psi);
76 
77  // calculate NoiseBefore Matrix R M R^T
78  double noiseBefore34 = sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseBefore_ij = noiseBefore_ji
79  double noiseBefore35 = -sigma2 * costheta * sinpsi * sintheta;
80  double noiseBefore45 = sigma2 * costheta * cospsi * sintheta;
81 
82  noiseBefore[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
83  noiseBefore[4][3] = noiseBefore34;
84  noiseBefore[5][3] = noiseBefore35;
85 
86  noiseBefore[3][4] = noiseBefore34;
87  noiseBefore[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
88  noiseBefore[5][4] = noiseBefore45;
89 
90  noiseBefore[3][5] = noiseBefore35;
91  noiseBefore[4][5] = noiseBefore45;
92  noiseBefore[5][5] = sigma2 * sintheta*sintheta;
93 
94  TMatrixT<Double_t> jacobianT(7,7);
95  jacobianT = (*jacobian);
96  jacobianT.T();
97 
98  noiseBefore = jacobianT*noiseBefore*(*jacobian); //propagate
99 
100  // noiseAfter
101  TMatrixT<Double_t> noiseAfter(7,7);
102 
103  // calculate euler angles theta, psi (so that A' points in z' direction)
104  psi = 0;
105  if (fabs((*directionAfter)[1]) < 1E-14) { // numerical case: arctan(+-inf)=+-PI/2
106  if ((*directionAfter)[0] >= 0.) psi = M_PI/2.;
107  else psi = 3.*M_PI/2.;
108  }
109  else {
110  if ((*directionAfter)[1] > 0.) psi = M_PI - atan((*directionAfter)[0]/(*directionAfter)[1]);
111  else psi = -atan((*directionAfter)[0]/(*directionAfter)[1]);
112  }
113  // cache sin and cos
114  sintheta = sqrt(1-(*directionAfter)[2]*(*directionAfter)[2]); // theta = arccos(directionAfter[2])
115  costheta = (*directionAfter)[2];
116  sinpsi = sin(psi);
117  cospsi = cos(psi);
118 
119  // calculate NoiseAfter Matrix R M R^T
120  double noiseAfter34 = sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseAfter_ij = noiseAfter_ji
121  double noiseAfter35 = -sigma2 * costheta * sinpsi * sintheta;
122  double noiseAfter45 = sigma2 * costheta * cospsi * sintheta;
123 
124  noiseAfter[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
125  noiseAfter[4][3] = noiseAfter34;
126  noiseAfter[5][3] = noiseAfter35;
127 
128  noiseAfter[3][4] = noiseAfter34;
129  noiseAfter[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
130  noiseAfter[5][4] = noiseAfter45;
131 
132  noiseAfter[3][5] = noiseAfter35;
133  noiseAfter[4][5] = noiseAfter45;
134  noiseAfter[5][5] = sigma2 * sintheta*sintheta;
135 
136  //calculate mean of noiseBefore and noiseAfter and update noise
137  (*noise) += 0.5*noiseBefore + 0.5*noiseAfter;
138  }
139 
140  return 0.;
141 }
E
Float_t E
Definition: plot.C:23
beta
Double_t beta
Definition: TG43_relative_dose.C:15
noise
unsigned int noise()
Definition: chem4.cc:265
genf::GFAbsEnergyLoss::getParticleParameters
void getParticleParameters(const int &pdg, double &charge, double &mass)
Gets particle charge and mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:27
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
GFException::setFatal
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:80
double genf::GFAbsEnergyLoss::getParticleMass ( const int &  pdg)
protectedinherited

Returns particle mass (in GeV)

Definition at line 35 of file GFAbsEnergyLoss.cxx.

References part.

Referenced by genf::GFEnergyLossBrems::energyLoss(), and genf::GFEnergyLossBetheBloch::energyLoss().

35  {
36  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(pdg);
37  return part->Mass();
38 }
part
TString part[npart]
Definition: Style.C:32
void genf::GFAbsEnergyLoss::getParticleParameters ( const int &  pdg,
double &  charge,
double &  mass 
)
protectedinherited

Gets particle charge and mass (in GeV)

Definition at line 27 of file GFAbsEnergyLoss.cxx.

References part.

Referenced by genf::GFEnergyLossBrems::energyLoss(), genf::GFEnergyLossBetheBloch::energyLoss(), and energyLoss().

29  {
30  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(pdg);
31  charge = part->Charge()/(3.);
32  mass = part->Mass();
33 }
part
TString part[npart]
Definition: Style.C:32
The documentation for this class was generated from the following files: