EnergyLossUpdator.h 8.35 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/

///////////////////////////////////////////////////////////////////
// EnergyLossUpdator.h, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////

#ifndef TRKEXTOOLS_ENERGYLOSSUPDATOR_H
#define TRKEXTOOLS_ENERGYLOSSUPDATOR_H

// Gaudi
#include "AthenaBaseComps/AthAlgTool.h"
// Trk
#include "TrkExInterfaces/IEnergyLossUpdator.h"
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkExUtils/MaterialInteraction.h"
// STL
#include <utility>

namespace Trk {

  class MaterialProperties;
  class EnergyLoss;
  
  /**@class EnergyLossUpdator

    The mean Energy loss per unit length of a particle is calculated as
    
    @f$ (dE/dx) = (dE/dX)_{io} + (dE/dX)_{rad} @f$
    
    Formulas used:
      - <b>For the ionization :</b>
        Bethe-Bloch for minimum ionizing particles (modified for electrons)
        
        @f$ (dE/dX)_{io,mip} = K \cdot \frac{Z}{A} \frac{1/\beta^2} \cdot [\frac{1}{2}]
        \cdot \ln \frac{2 m_e c^2 \beta^2 \gamma^2 T_max}{I^2} - \beta^2] @f$, (1)
        
        respectively
        
        @f$ (dE/dX)_{io,e}  = \frac{1/2} \cdot K \cdot \frac{Z}{A} 
        [2\ln\frac{2m_e c^2}{I} + 3 \ln{\gamma} - 1.95] @f$ (2)
    
        Alternatively, the most probable value instead of the mean value can be taken (steared by boolean in
        the method signature), taken from (4)

         @f$ \_{L}\Delta_p = \xi \left[ln\frac{2mc^2\beta^2\gamma^2}{I} + ln\frac{\xi}{I} - 0.8 + 4.447 \right] @f$

        with @f$ \xi = Z N_a \frac{k}{\beta^2}\codt t @f$

        The with of both distributions is given as the (crude) approximation, 
        that the full width at half maximum of the landau distribution is :  @f$ 4 \xi @f$

        For the transformation into a parameterized sigma the factor @f$ 2\cdot\sqrt{2.\cdot ln(2.) } @f$ can be used


        Sigma can be multiplied by a form factor (given through job options)
    
     - <b>For radiation corrections: </b>

        The Bethe-Heitler parameterisation is in terms of

	     @f$ Z = \frac{E_f}{E_i} @f$ (2)

	     With this parameterisation the expected energy loss

      @f$ \langle Z \rangle = e^{-t} @f$ (2)

	    where @f$ t = \frac{x}{X_0} @f$ or the fraction of a radiation length of material traversed by the particle.

	    This formula can be transposed in terms of @f$ \Delta (E) @f$. The transposition however depends on the direction of propagation.

	    For forwards propagation we get: @f$ \langle \Delta (E) \rangle = E_i ( 1 - \langle Z \rangle ) @f$ (2)

	    For the propagation against momentum we get: @f$ \langle \Delta (E) \rangle = E_i ( \frac{1}{\langle Z \rangle} - 1 ) @f$ (2)

    	The variance of Z is also defined (2): @f$ Var(Z) = e^{-t ln(3) / ln(2)} - e^{-2t} @f$

	    Transposing this one dervies the following relations for 

	    1. Forwards propagation @f$ Var ( \Delta \frac{q}{p} ) = \frac{Var(Z)}{\langle Z \rangle ^2 p^2} @f$

    	2. Backwards propagation @f$ Var ( \Delta \frac{q}{p} ) = \frac{Var(Z)}{p^2} @f$

    	Where p is the global momentum.

        Bethe-Heitler formula
        
        @f$ (dE/dX)_{rad} = \frac{E_{initial}}{X_o}\cdot (\frac{m}{m_e})^2 @f$ (3)
        
      References:
      (1)      Review of Particle Physics, ELSEVIR, ISSN 0370-2639 (2004)
      (2),(3)  D. Stampfer et al., Track fitting with energy loss, Comp. Phys. Comm 79 (1994) 
      (4)      H. Bichsel, Straggling in thin silicon detectors, Rev. Mod. Phys., Vol. 60., No. 3, July 1988
        
             
     @author Andreas.Salzburger@cern.ch, Tom.Atkinson@cern.ch
    */
    
  class EnergyLossUpdator : public AthAlgTool,
                            virtual public IEnergyLossUpdator {
      
    public:
      /** AlgTool like constructor */
      EnergyLossUpdator(const std::string&,const std::string&,const IInterface*);

      /**Virtual destructor*/
      virtual ~EnergyLossUpdator();
       
      /** AlgTool initailize method.*/
      StatusCode initialize();
      /** AlgTool finalize method */
      StatusCode finalize();
        
      /** dEdX calculation when providing MaterialProperties,
        a momentum, a pathlength, and a ParicleHypothesis:
        
        Units: [MeV/mm]
        */
      double dEdX(const MaterialProperties& mat,
                  double p,
                  ParticleHypothesis particle=pion) const;
                  
      /** deltaE calculation
        using dEdX and integrating along pathlength,
        assuming constant dEdX during for the path.
        The sign depends on the given propagation direction 
        
        Units: [MeV]

        mpv steers the most probable energy loss
        */
     EnergyLoss* energyLoss(
			    const MaterialProperties& mat,
			    double p,
			    double pathcorrection,
			    PropDirection dir=alongMomentum,
			    ParticleHypothesis particle=pion,
			    bool mpv = false) const;   
    

143
144
145
146
147
148
149
150
     /** Method to recalculate Eloss values for the fit setting an elossFlag using as an input
         the detailed Eloss information Calorimeter energy, error momentum and momentum error */
     EnergyLoss* updateEnergyLoss(EnergyLoss* eLoss, double caloEnergy, double caloEnergyError, 
                          double pCaloEntry, double momentumError, int & elossFlag) const;

     /** Routine to calculate X0 and Eloss scale factors for the Calorimeter and Muon System */
     void getX0ElossScales(int icalo, double eta, double phi, double & X0Scale, double & ElossScale ) const;

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
      /** Method to return the variance of the change in q/p for the Bethe-Heitler parameterisation */
      double varianceDeltaQoverP(const MaterialProperties&,
                                 double p,
                                 double pathcorrection,
                                 PropDirection direction = alongMomentum,
                                 ParticleHypothesis particleHypothesis = electron ) const;


    private:

      Trk::MaterialInteraction       m_matInt;
      /** dEdX BetheBloch calculation: 
        Units: [MeV]
        */
      double dEdXBetheBloch(const MaterialProperties& mat,
                            double gamma,
                            double beta,
                            ParticleHypothesis particle=pion) const;
                            
      /** dEdX BetheHeitler calculation: 
        Units: [MeV]
        */
      double dEdXBetheHeitler(const MaterialProperties& mat,
                              double initialE,
                              ParticleHypothesis particle=pion) const;
 
      bool                    m_useTrkUtils;               //!< use eloss parametrisation from TrkUtils MaterialInterAction.h
      bool                    m_gaussianVavilovTheory;     //!< include energy loss straggling or not
      bool                    m_useBetheBlochForElectrons; //!< use adopted bethe bloch for electrons 
      mutable double          m_transKaz;                  //!< transport kaz for straggling 
      mutable double          m_transTmax;                 //!< transport Tmax for straggling
      double                  m_stragglingErrorScale;      //!< stragglingErrorScale
      double                  m_mpvScale;                  //!< a scalor that can be introduced for the MPV
      bool                    m_mpvSigmaParametric;        //!< take the (crude) parametric mpv sigma 
      bool                    m_detailedEloss;             //!< provide extended EnergyLoss info 

      static double           s_ka_BetheBloch;          //!< KOverA factor in Bethe-Bloch equation [MeV*cm2/gram]
      static double           s_eulerConstant;          //!< the euler constant for mip calculations

      static double           s_fwhmToSigma;            //!< 1./[2.*sart(2.*ln(2.))] for fwhm -> sigma

      static double           s_mpv_p0;                 //! variable obtained from a fit to landau distribution for sigma(landau)
      static double           s_mpv_p1;                 //! variable obtained from a fit to landau distribution for sigma(landau)
      static double           s_mpv_p2;                 //! variable obtained from a fit to landau distribution for sigma(landau)

      static ParticleMasses   s_particleMasses;         //!< struct of Particle masses                                                        
                                                      
  };


} // end of namespace


#endif // TRKEXTOOLS_ENERGYLOSSUPDATOR_H