LArSoft  v09_90_00
Liquid Argon Software toolkit - https://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 42 of file GFEnergyLossCoulomb.h.

Constructor & Destructor Documentation

genf::GFEnergyLossCoulomb::~GFEnergyLossCoulomb ( )
virtual

Definition at line 25 of file GFEnergyLossCoulomb.cxx.

25 {}

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 27 of file GFEnergyLossCoulomb.cxx.

References 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)
49  throw GFException(std::string(__func__) + ": no noise given!", __LINE__, __FILE__).setFatal();
50  if (!jacobian)
51  throw GFException(std::string(__func__) + ": no jacobian given!", __LINE__, __FILE__)
52  .setFatal();
53  if (!directionBefore)
54  throw GFException(std::string(__func__) + ": no directionBefore given!", __LINE__, __FILE__)
55  .setFatal();
56  if (!directionAfter)
57  throw GFException(std::string(__func__) + ": no directionAfter given!", __LINE__, __FILE__)
58  .setFatal();
59 
60  // MULTIPLE SCATTERING; calculate sigma^2
61  // PANDA report PV/01-07 eq(43); linear in step length
62  double sigma2 =
63  225.E-6 / (beta * beta * mom * mom) * step / radiationLength * matZ / (matZ + 1) *
64  log(159. * pow(matZ, -1. / 3.)) /
65  log(
66  287. /
67  sqrt(
68  matZ)); // sigma^2 = 225E-6/mom^2 * XX0/beta^2 * Z/(Z+1) * ln(159*Z^(-1/3))/ln(287*Z^(-1/2)
69 
70  // noiseBefore
71  TMatrixT<Double_t> noiseBefore(7, 7);
72 
73  // calculate euler angles theta, psi (so that directionBefore' points in z' direction)
74  double psi = 0;
75  if (fabs((*directionBefore)[1]) < 1E-14) { // numerical case: arctan(+-inf)=+-PI/2
76  if ((*directionBefore)[0] >= 0.)
77  psi = M_PI / 2.;
78  else
79  psi = 3. * M_PI / 2.;
80  }
81  else {
82  if ((*directionBefore)[1] > 0.)
83  psi = M_PI - atan((*directionBefore)[0] / (*directionBefore)[1]);
84  else
85  psi = -atan((*directionBefore)[0] / (*directionBefore)[1]);
86  }
87  // cache sin and cos
88  double sintheta =
89  sqrt(1 - (*directionBefore)[2] * (*directionBefore)[2]); // theta = arccos(directionBefore[2])
90  double costheta = (*directionBefore)[2];
91  double sinpsi = sin(psi);
92  double cospsi = cos(psi);
93 
94  // calculate NoiseBefore Matrix R M R^T
95  double noiseBefore34 =
96  sigma2 * cospsi * sinpsi * sintheta * sintheta; // noiseBefore_ij = noiseBefore_ji
97  double noiseBefore35 = -sigma2 * costheta * sinpsi * sintheta;
98  double noiseBefore45 = sigma2 * costheta * cospsi * sintheta;
99 
100  noiseBefore[3][3] =
101  sigma2 * (cospsi * cospsi + costheta * costheta - costheta * costheta * cospsi * cospsi);
102  noiseBefore[4][3] = noiseBefore34;
103  noiseBefore[5][3] = noiseBefore35;
104 
105  noiseBefore[3][4] = noiseBefore34;
106  noiseBefore[4][4] = sigma2 * (sinpsi * sinpsi + costheta * costheta * cospsi * cospsi);
107  noiseBefore[5][4] = noiseBefore45;
108 
109  noiseBefore[3][5] = noiseBefore35;
110  noiseBefore[4][5] = noiseBefore45;
111  noiseBefore[5][5] = sigma2 * sintheta * sintheta;
112 
113  TMatrixT<Double_t> jacobianT(7, 7);
114  jacobianT = (*jacobian);
115  jacobianT.T();
116 
117  noiseBefore = jacobianT * noiseBefore * (*jacobian); //propagate
118 
119  // noiseAfter
120  TMatrixT<Double_t> noiseAfter(7, 7);
121 
122  // calculate euler angles theta, psi (so that A' points in z' direction)
123  psi = 0;
124  if (fabs((*directionAfter)[1]) < 1E-14) { // numerical case: arctan(+-inf)=+-PI/2
125  if ((*directionAfter)[0] >= 0.)
126  psi = M_PI / 2.;
127  else
128  psi = 3. * M_PI / 2.;
129  }
130  else {
131  if ((*directionAfter)[1] > 0.)
132  psi = M_PI - atan((*directionAfter)[0] / (*directionAfter)[1]);
133  else
134  psi = -atan((*directionAfter)[0] / (*directionAfter)[1]);
135  }
136  // cache sin and cos
137  sintheta =
138  sqrt(1 - (*directionAfter)[2] * (*directionAfter)[2]); // theta = arccos(directionAfter[2])
139  costheta = (*directionAfter)[2];
140  sinpsi = sin(psi);
141  cospsi = cos(psi);
142 
143  // calculate NoiseAfter Matrix R M R^T
144  double noiseAfter34 =
145  sigma2 * cospsi * sinpsi * sintheta * sintheta; // noiseAfter_ij = noiseAfter_ji
146  double noiseAfter35 = -sigma2 * costheta * sinpsi * sintheta;
147  double noiseAfter45 = sigma2 * costheta * cospsi * sintheta;
148 
149  noiseAfter[3][3] =
150  sigma2 * (cospsi * cospsi + costheta * costheta - costheta * costheta * cospsi * cospsi);
151  noiseAfter[4][3] = noiseAfter34;
152  noiseAfter[5][3] = noiseAfter35;
153 
154  noiseAfter[3][4] = noiseAfter34;
155  noiseAfter[4][4] = sigma2 * (sinpsi * sinpsi + costheta * costheta * cospsi * cospsi);
156  noiseAfter[5][4] = noiseAfter45;
157 
158  noiseAfter[3][5] = noiseAfter35;
159  noiseAfter[4][5] = noiseAfter45;
160  noiseAfter[5][5] = sigma2 * sintheta * sintheta;
161 
162  //calculate mean of noiseBefore and noiseAfter and update noise
163  (*noise) += 0.5 * noiseBefore + 0.5 * noiseAfter;
164  }
165 
166  return 0.;
167 }
noise
unsigned int noise()
Definition: chem4.cc:261
E
Float_t E
Definition: plot.C:20
genf::GFAbsEnergyLoss::getParticleParameters
void getParticleParameters(const int &pdg, double &charge, double &mass)
Gets particle charge and mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:26
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
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:75
double genf::GFAbsEnergyLoss::getParticleMass ( const int &  pdg)
protectedinherited

Returns particle mass (in GeV)

Definition at line 33 of file GFAbsEnergyLoss.cxx.

References part.

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

34 {
35  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg);
36  return part->Mass();
37 }
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 26 of file GFAbsEnergyLoss.cxx.

References part.

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

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