EnergyLossUpdator.cxx 20.3 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
/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/

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

// Trk include
#include "TrkExTools/EnergyLossUpdator.h"
#include "TrkGeometry/MaterialProperties.h"
#include "TrkMaterialOnTrack/EnergyLoss.h"
// Gaudi
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SystemOfUnits.h"

#include <math.h>

// static particle masses
Trk::ParticleMasses Trk::EnergyLossUpdator::s_particleMasses;
// statics doubles 
double Trk::EnergyLossUpdator::s_ka_BetheBloch = 30.7075 * Gaudi::Units::MeV;
double Trk::EnergyLossUpdator::s_eulerConstant = 0.577215665; // as given by google.com ;-)

// from full with at half maximum to sigma for a gaussian
double Trk::EnergyLossUpdator::s_fwhmToSigma   = 0.424;  //   1./(2.*sqrt(2.*log(2.)));

// mpv to sigma fit values
double Trk::EnergyLossUpdator::s_mpv_p0        = 4.57270e-02;
double Trk::EnergyLossUpdator::s_mpv_p1        = 8.11761e-03;
double Trk::EnergyLossUpdator::s_mpv_p2        = -4.85133e-01;

// constructor
Trk::EnergyLossUpdator::EnergyLossUpdator(const std::string& t, const std::string& n, const IInterface* p) :
  AthAlgTool(t,n,p),
  m_useTrkUtils(true),
  m_gaussianVavilovTheory(false),
  m_useBetheBlochForElectrons(true),
  m_transKaz(0.),
  m_transTmax(0.),
  m_stragglingErrorScale(1.),
  m_mpvScale(0.98),
  m_mpvSigmaParametric(false),
44
  m_detailedEloss(true)
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
{
   declareInterface<Trk::IEnergyLossUpdator>(this);
   // scale from outside
   declareProperty("UseTrkUtils"                , m_useTrkUtils);           
   declareProperty("BetheBlochKAFactor"         , s_ka_BetheBloch);
   // some additional setup
   declareProperty("UseGaussVavilovTheory"      , m_gaussianVavilovTheory);
   declareProperty("UseBetheBlochForElectrons"  , m_useBetheBlochForElectrons);
   // scalor for the straggling
   declareProperty("ScaleStragglingError"       , m_stragglingErrorScale);
   // scalor for the most probable value
   declareProperty("UseParametricMpvError"      , m_mpvSigmaParametric);
   declareProperty("ScaleMostProbableValue"     , m_mpvScale);
   declareProperty("DetailedEloss"              , m_detailedEloss);
}

// destructor
Trk::EnergyLossUpdator::~EnergyLossUpdator()
{}

// Athena standard methods
// initialize
StatusCode Trk::EnergyLossUpdator::initialize()
{
    ATH_MSG_INFO( "initialize()" );    
    return StatusCode::SUCCESS;
}

// finalize
StatusCode Trk::EnergyLossUpdator::finalize()
{
    ATH_MSG_INFO( "finalize() successful" );
    return StatusCode::SUCCESS;
}

// public interface method
double Trk::EnergyLossUpdator::dEdX(const MaterialProperties& mat,
                                    double p,
                                    ParticleHypothesis particle) const
{
  if (particle == Trk::undefined || particle == Trk::nonInteracting) return 0.; 
   
  // preparation of kinetic constants
  double m     = s_particleMasses.mass[particle];
  double E     = sqrt(p*p+m*m);
  double beta  = p/E;
  double gamma = E/m;

  m_transKaz  = 0.;
  m_transTmax = 0.;
  // add ionization and radiation
  double dEdX = dEdXBetheBloch(mat, beta, gamma, particle) + dEdXBetheHeitler(mat, E, particle);

  // add e+e- pair production and photonuclear effect for muons at energies above 8 GeV
  if ((particle == Trk::muon) && (E > 8000.)) {
    if (E < 1.e6) {
      dEdX += - 0.5345/mat.x0() + 6.803e-5*E/mat.x0() + 2.278e-11*E*E/mat.x0() - 9.899e-18*E*E*E/mat.x0(); //E below 1 TeV
    } else {
      dEdX += - 2.986/mat.x0() + 9.253e-5*E/mat.x0(); //E above 1 TeV
    }
  }

  return dEdX;
}

   
// public interface method
Trk::EnergyLoss* Trk::EnergyLossUpdator::energyLoss(const MaterialProperties& mat,
                                                    double p,
                                                    double pathcorrection,
                                                    PropDirection dir,
                                                    ParticleHypothesis particle,
                                                    bool mpvSwitch) const
{ 
 
  if (particle == Trk::undefined ) {
    ATH_MSG_WARNING( "undefined ParticleHypothesis, energy loss calculation cancelled" );
    return 0;
  }

  double deltaE       = 0.;
  double sigmaDeltaE  = 0.;
  m_transKaz  = 0.;
  m_transTmax = 0.;

  // preparation 
131
  double sign = (dir==Trk::oppositeMomentum) ? -1. : 1.;
132
133
134
135
136
137
138
139
140

  double pathLength = pathcorrection * mat.thicknessInX0()*mat.x0();  

  double sigIoni = 0.;
  double sigRad = 0.;
  double kazL = 0.;
  double meanIoni  = m_matInt.dEdl_ionization(p, &(mat.material()), particle, sigIoni, kazL);
  double meanRad  = m_matInt.dEdl_radiation(p, &(mat.material()), particle, sigRad);

141
142
  meanIoni = sign*pathLength*meanIoni;  
  meanRad = sign*pathLength*meanRad;  
143
144
145
146
147
148
149
150
151
152
  sigIoni = pathLength*sigIoni;  
  sigRad  = pathLength*sigRad;  
  kazL    = pathLength*kazL;
//
// include pathlength dependence of Landau ionization 
//
  sigIoni = sigIoni - kazL*log(pathLength);  

  deltaE = meanIoni + meanRad;
  sigmaDeltaE = sqrt(sigIoni*sigIoni+sigRad*sigRad);
153
  ATH_MSG_DEBUG( " Energy loss updator deltaE " << deltaE << " meanIoni " << meanIoni << " meanRad " << meanRad << " sigIoni " << sigIoni << " sigRad " << sigRad << " sign " << sign  << " pathLength " << pathLength );   
154
155

  Trk::EnergyLoss* eloss = !m_detailedEloss ? new Trk::EnergyLoss(deltaE,sigmaDeltaE):
156
     new Trk::EnergyLoss(deltaE, sigmaDeltaE, sigmaDeltaE, sigmaDeltaE, meanIoni, sigIoni, meanRad, sigRad, pathLength );
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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
  if(m_useTrkUtils) return eloss;

// Code below will not be used if the parameterization of TrkUtils is used
  
  double m     = s_particleMasses.mass[particle]; 
  double mfrac = s_particleMasses.mass[Trk::electron]/m;
  double mfrac2 = mfrac*mfrac;
  double E     = sqrt(p*p+m*m);
  // relativistic properties
  double beta  = p/E;
  double gamma = E/m;
  
  // mean radiation loss + sigma
  // Follow CompPhys Comm 79 (1994) P157-164
  // Brem parameterization in terms of path-length in number of radiation lengths
  double dInX0  = pathcorrection*mat.thicknessInX0();
  double meanZ = exp(-mfrac2*dInX0);

  double deltaE_rad =  (dir == Trk::alongMomentum )? sign * E * (1. - meanZ) : sign * E * (1. / meanZ - 1.);

  // add e+e- pair production and photonuclear effect for muons at energies above 8 GeV
  if ((particle == Trk::muon) && (E > 8000.)) {
    double deltaEpair = 0.;
    if (E < 1.e6) {
      deltaEpair = (- 0.5345 + 6.803e-5*E + 2.278e-11*E*E - 9.899e-18*E*E*E)*dInX0; //E below 1 TeV
   } else {
      deltaEpair = (- 2.986 + 9.253e-5*E)*dInX0; //E above 1 TeV
    }
    deltaE_rad += sign*deltaEpair;
  }

//  The sigma radiation loss of CompPhys Comm 79 (1994) P157-164 has PROBLEMS at low Eloss values 

//  double varZ  = exp( -1. * dInX0 * log(3.) / log(2.) ) - exp(-2. * dInX0);
//  double varianceQoverP =  ( dir == Trk::alongMomentum )? 1. / (meanZ * meanZ * p * p) * varZ : varZ / (p * p);
//  double sigmaDeltaE_rad = beta*beta*p*p*sqrt( varianceQoverP );

// Use a simple scale factor applied to the Eloss from Radiation obtained from G4 simulation in stead

  double Radiation_FWHM = m_stragglingErrorScale*4*0.85*fabs(deltaE_rad)/3.59524;
  double sigmaDeltaE_rad = s_fwhmToSigma*Radiation_FWHM;

  // ionization 
  double deltaE_ioni = 0.;
  double sigmaDeltaE_ioni = 0.;
  if (particle != Trk::electron || m_useBetheBlochForElectrons) {

    if (!mpvSwitch) {
      // use the dEdX function 
      deltaE_ioni = sign*pathLength*dEdXBetheBloch(mat, beta, gamma, particle);

      // the different straggling functions
      if (m_gaussianVavilovTheory)     {
	// use the Gaussian approximation for the Vavilov Theory
	double sigmaE2 =  m_transKaz * pathLength * m_transTmax*(1.-beta*beta/2.);
	sigmaDeltaE_ioni = sqrt(sigmaE2);
      } else 

//      Take FWHM maximum of Landau and convert to Gaussian  
//      For the FWHM of the Landau Bichsel/PDG is used: FWHM = 4 xi = 4 m_transKaz * pathLength
	sigmaDeltaE_ioni = 4.*s_fwhmToSigma* m_transKaz * pathLength; 
      
    } else {

      double eta2   = beta*gamma;   
      eta2 *= eta2;   
      // most probable  
      double kaz = 0.5*s_ka_BetheBloch*mat.zOverAtimesRho();
      double xi = kaz * pathLength;
      // get the ionisation potential
      double iPot = 16.e-6 * std::pow( mat.averageZ(),0.9);
      // density effect, only valid for high energies (gamma > 10 -> p > 1GeV for muons)
      double delta = 0.;
      /* ST replace with STEP-like coding
	 if (eta2 > 2.) {
	 // high energy density effect
	 double eplasma = 28.816e-6 * sqrt(1000.*0.5);
	 delta = 2.*log(eplasma/iPot) + log(eta2) - 0.5;
	 }
      */
      if (gamma>10.) {
	double eplasma = 28.816e-6 * sqrt(1000.*mat.zOverAtimesRho());
	delta = 2.*log(eplasma/iPot) + log(eta2) - 1.;
      }	 
      // calculate the most probable value of the Landau distribution
      double mpv = m_mpvScale * xi/(beta*beta) * (log(m*eta2*kaz/(iPot*iPot)) + 0.2-beta*beta-delta);//12.325);

      // sigma 
      // (1) - obtained from fitted function  [  return par0 + par1*std::pow(x,par2); ]
      // (2) - following Bichsel: fwhm (full with at half maximum) = 4 *xi, w = 2.35 fwhm (for gaussian)  
      // use the crude approximation for the width 4*chi
      sigmaDeltaE_ioni = !m_mpvSigmaParametric ? s_fwhmToSigma*4.*xi : mpv*(s_mpv_p0 + s_mpv_p1*std::pow(mpv,s_mpv_p2));
  
      // get the fitted sigma        
      deltaE_ioni = sign* mpv;

    }
    // scale the error
    sigmaDeltaE_ioni *= m_stragglingErrorScale; 
  }

  deltaE = deltaE_ioni + deltaE_rad;

  sigmaDeltaE = sqrt( sigmaDeltaE_rad*sigmaDeltaE_rad + sigmaDeltaE_ioni*sigmaDeltaE_ioni);

262
  return ( m_detailedEloss ? new EnergyLoss(deltaE, sigmaDeltaE, sigmaDeltaE, sigmaDeltaE,
263
					    (mpvSwitch? deltaE_ioni :0.9*deltaE_ioni), sigmaDeltaE_ioni,
264
					    deltaE_rad, sigmaDeltaE_rad, pathLength ):
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
	                     new EnergyLoss(deltaE, sigmaDeltaE) );

}

double Trk::EnergyLossUpdator::dEdXBetheBloch(const MaterialProperties& mat,
                                              double beta,
                                              double gamma,
                                              ParticleHypothesis particle) const
{

  if (particle == Trk::undefined || particle == Trk::nonInteracting ) return 0.;

  if (mat.averageZ()==0. || mat.zOverAtimesRho()==0. ) {
    ATH_MSG_ERROR("empty material properties pass to the EnergyLossUpdator:Z,zOAtr:"<<mat.averageZ()<<","<<mat.zOverAtimesRho());
    return 0.;
  }

  // 16 eV * Z**0.9 - bring to MeV
  double iPot = 16.e-6 * std::pow(mat.averageZ(),0.9);
  // and the electron mass in MeV
  double me    = s_particleMasses.mass[Trk::electron];
  double m     = s_particleMasses.mass[particle];
  double eta2  = beta*gamma;
 
  // K/A*Z = 0.5 * 30.7075MeV/(g/mm2) * Z/A * rho[g/mm3]  / scale to mm by this
  double kaz = 0.5*s_ka_BetheBloch*mat.zOverAtimesRho();
     
  if (particle != Trk::electron){

    // density effect, only valid for high energies (gamma > 10 -> p > 1GeV for muons)
    double delta = 0.;

    /* ST replace with STEP-like coding  
    // high energy density effect --- will be ramped up linearly
    double eplasma = 28.816e-6 * sqrt(1000.*0.5);
    delta = 2.*log(eplasma/iPot) + log(eta2) - 0.5;
    if (eta2 < 100.){
    delta *= (eta2-3.)/97.; 
    }
    */
      
    eta2 *= eta2; 

    if (gamma>10.) {
      double eplasma = 28.816e-6 * sqrt(1000.*mat.zOverAtimesRho());
      delta = 2.*log(eplasma/iPot) + log(eta2) - 1.;
    }
    
    // mass fraction
    double mfrac = me/m;
    // tmax - cut off energy
    double tMax = 2.*eta2*me/(1.+2.*gamma*mfrac+mfrac*mfrac);
    // divide by beta^2 for non-electrons
    kaz /= beta*beta;
    // store the transport variables
    m_transKaz  = kaz;
    m_transTmax = tMax;
    // return
    return  kaz*(log(2.*me*eta2*tMax/(iPot*iPot)) - 2.*(beta*beta) - delta);
  
  }
  m_transKaz = kaz;
  // for electrons use slightly different BetheBloch adaption
  // see Stampfer, et al, "Track Fitting With Energy Loss", Comp. Pyhs. Comm. 79 (1994), 157-164
  return kaz*(2.*log(2.*me/iPot)+3.*log(gamma) - 1.95);
  
}

double Trk::EnergyLossUpdator::dEdXBetheHeitler(const MaterialProperties& mat,
                                                double initialE, 
                                                ParticleHypothesis particle) const
{

  if (particle == Trk::undefined || particle == Trk::nonInteracting ) return 0.;

  double mfrac = (s_particleMasses.mass[Trk::electron]/s_particleMasses.mass[particle]); mfrac *= mfrac;
  
  return initialE/mat.x0()*mfrac;
}

345
346
// public interface method
Trk::EnergyLoss* Trk::EnergyLossUpdator::updateEnergyLoss(Trk::EnergyLoss* eLoss, double caloEnergy, double caloEnergyError, double pCaloEntry, double momentumError, int & elossFlag) const {
347
348

//
349
350
351
352
353
// Input: the detailed EnergyLoss object in the Calorimeter that contains the different Eloss terms 
// and their uncertainties; caloEnergy and error; and the muon momentumError (all in MeV)
//
// For use in the MuonSystem
// Input: caloEnergy = 0. caloEnergyError = 0. and pCaloEntry = pMuonEntry momentum at MuonEntry
354
355
//
// Output: an updated Energy loss values deltaE() 
356
357
358
359
360
361
//         that can be used in the track fit and corresponds to the Most Probable EnergyLoss value
//         taking into account the ionization, radiation and 
//         smearing due to the errors including the momentumError (in MeV)
//
//         elossFlag = false if Calorimeter Energy is NOT stored (and later fitted) on the Eloss object
//                   = true  Calorimeter Energy is stored and will be fitted  
362
363
364
365
//
//  deltaE is used in the final fit
//

366
367
  elossFlag = 0;

368
369
370
371
372
373
374
  int isign = 1;
  if(eLoss->deltaE()<0) isign = -1;
 
  double deltaE = eLoss->deltaE();
  double sigmaDeltaE = eLoss->sigmaDeltaE();
// Detailed Eloss
  double deltaE_ioni = eLoss->meanIoni(); 
375
376
  double sigmaDeltaE_ioni = 0.45*eLoss->sigmaIoni(); // sigma Landau
  double deltaE_rad = eLoss->meanRad();  
377
  double sigmaDeltaE_rad = eLoss->sigmaRad(); // rms and mean of steep exponential
378
  double depth = eLoss->length();
379
380
381
382
383

  double sigmaPlusDeltaE = eLoss->sigmaPlusDeltaE();
  double sigmaMinusDeltaE = eLoss->sigmaMinusDeltaE();

  double MOP = deltaE_ioni - isign*3.59524*sigmaDeltaE_ioni;
384
385

//  std::cout << " update Energyloss old deltaE " << deltaE << " sigmaPlusDeltaE " << sigmaPlusDeltaE << " sigmaMinusDeltaE " << sigmaMinusDeltaE << std::endl;
386
387
388
//
// MOP shift due to ionization and radiation
//   
389
  double MOPshift = isign*50*10000./pCaloEntry + isign*0.75*sqrt(sigmaDeltaE_ioni*sigmaDeltaE_rad);
390
391
392
//
// define sigmas for Landau convoluted with exponential
//
393
394
  double fracErad = sigmaDeltaE_rad + fabs(deltaE_rad)*pCaloEntry/(800000.+pCaloEntry);
  double sigmaL = sigmaDeltaE_ioni + fracErad/3.59524;
395
396
397
398
399
400
401
402
403
404
405
  double sigmaMinus = 1.02*sigmaDeltaE_ioni + 0.08*sigmaDeltaE_rad; 
  double sigmaPlus = 4.65*sigmaDeltaE_ioni + 1.16*sigmaDeltaE_rad;

//
// Case where the measured Calorimeter energy is not available
//
  
  if(caloEnergyError<=0) {
//
// Shift of MOP due to momentum resolution 
//
406
407
    double scale_xc = 2.3;
    double xc = scale_xc*0.87388*momentumError/(3.59524*sigmaL);
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
    double correction = (1.747*xc*xc + 0.97*0.938*xc*xc*xc)/(1+4.346*xc+5.371*xc*xc+0.938*xc*xc*xc); // correction ranges from 0 to 0.97
    double MOPreso = isign*3.59524*sigmaL*correction;

    deltaE = MOP + MOPshift + MOPreso;
    sigmaMinusDeltaE = sigmaMinus;
    sigmaPlusDeltaE  = sigmaPlus;
    sigmaDeltaE = sqrt(0.5*sigmaMinusDeltaE*sigmaMinusDeltaE+0.5*sigmaPlusDeltaE*sigmaPlusDeltaE);

  } else {

    double sigmaPlusTot = sqrt(sigmaPlus*sigmaPlus+caloEnergyError*caloEnergyError);

    if(caloEnergy > fabs(MOP + MOPshift) + 2*sigmaPlusTot) {
//
// Use measured Calorimeter energy
//
      deltaE = isign*caloEnergy;
      sigmaMinusDeltaE = caloEnergyError;  
      sigmaPlusDeltaE = caloEnergyError;
      sigmaDeltaE = caloEnergyError;
428
      elossFlag = 1;
429
430
431

    } else {

432
// Use MOP after corrections
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457

//
// MOPCalo is correction to MOP for Calorimeter energy cut
//
      double xe = caloEnergyError/sigmaL;
      double MOPCalo = -isign*sigmaL*0.1*(1+xe)/(1+0.077*xe+0.077*xe*xe); // ranges from -0.1 sigmaL to 0.
//
// Shift of MOP due to momentum resolution smearing
//
      double errorScale = sqrt(1./4.65/4.65 + caloEnergyError*caloEnergyError/(caloEnergyError*caloEnergyError+sigmaPlus*sigmaPlus)); // ranges from 1/4.65 to 1 
      double xc = 0.87388*momentumError/(3.59524*sigmaL*errorScale);
      double correction = (1.747*xc*xc + 0.97*0.938*xc*xc*xc)/(1+4.346*xc+5.371*xc*xc+0.938*xc*xc*xc); // correction ranges from 0 to 0.97
      double MOPreso = isign*3.59524*sigmaL*errorScale*correction;
// 
// Use MOP after correction for radiation (MOPshift) and Calo cut (MOPCalo) and momentum resolution smearing (MOPreso) 
//
      deltaE = MOP + MOPshift + MOPCalo + MOPreso;
      sigmaMinusDeltaE = sigmaMinus;
      sigmaPlusDeltaE = sigmaPlus*errorScale; 
      sigmaDeltaE = sqrt(0.5*sigmaMinusDeltaE*sigmaMinusDeltaE+0.5*sigmaPlusDeltaE*sigmaPlusDeltaE);

   }

  }

458
459
460
461
462
//  std::cout << " update Energyloss new deltaE " << deltaE << " sigmaPlusDeltaE " << sigmaPlusDeltaE << " sigmaMinusDeltaE " << sigmaMinusDeltaE << std::endl;
 
  return new EnergyLoss(deltaE, sigmaDeltaE, sigmaMinusDeltaE, sigmaPlusDeltaE, deltaE_ioni, sigmaDeltaE_ioni, deltaE_rad, sigmaDeltaE_rad, depth );

}  
463

464
465
// public interface method
  void Trk::EnergyLossUpdator::getX0ElossScales(int icalo, double eta, double phi, double & X0Scale, double & ElossScale ) const {
466

467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
 //
 // for Calorimeter icalo = 1 
 //     Muon System icalo = 0
 // convention eta, phi is at Calorimeter Exit (or Muon Entry)
 //            eta and phi are from the position (not direction)
 //
 // input   X0 and ElossScale = 1
 // output  updated X0Scale and ElossScale
 //

  double X0CaloGirder[4] = {-1.02877e-01, -2.74322e-02, 8.12989e-02,9.73551e-01};

  double X0CaloScale[60] = {1.03178 ,1.03178 ,1.03178 ,1.03178 ,1.01512 ,1.01456 ,1.01266 ,1.02299 ,1.01437 ,1.01561 ,1.01731 ,1.02502 ,1.01171 ,0.987469 ,0.979134 ,1.01018 ,1.01813 ,0.986147 ,0.973686 ,0.983771 ,0.960334 ,0.967629 ,0.982935 ,0.993162 ,1.00241 ,1.00012 ,0.9814 ,0.937632 ,0.927918 ,0.93491 ,0.934184 ,0.921438 ,0.933788 ,0.954208 ,0.973684 ,0.977765 ,1.00697 ,1.01779 ,1.02017 ,1.01604 ,1.02121 ,1.02344 ,1.02888 ,1.02278 ,1.03096 ,1.02409 ,1.03149 ,1.03797 ,1.04254 ,1.04943 ,1.02433 ,1.02379 ,1.02177 ,1.02271 ,1.02011 ,1.02414 ,1.03649 ,1.03739 ,1.03739 ,1.03739};

  double ElossCaloScale[30] = { 1.03504 ,1.03504 ,1.03504 ,1.02631 ,1.04258 ,1.03506 ,1.05307 ,1.02399 ,1.04237 ,1.04368 ,1.04043 ,1.05899 ,1.07933 ,1.08604 ,1.08984 ,1.03564 ,1.04158 ,1.05983 ,1.06291 ,1.06853 ,1.0674 ,1.05427 ,1.06466 ,1.06274 ,1.06141 ,1.06314 ,1.06868 ,1.07242 ,1.07242 ,1.07242};
482
  
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
  double X0MuonScale[60] = {-0.0320612 ,-0.0320612 ,-0.0320612 ,-0.0320612 ,-0.0693796 ,-0.0389677 ,-0.0860891 ,-0.124606 ,-0.0882329 ,-0.100014 ,-0.0790912 ,-0.0745538 ,-0.099088 ,-0.0933711 ,-0.0618782 ,-0.0619762 ,-0.0658361 ,-0.109704 ,-0.129547 ,-0.143364 ,-0.0774768 ,-0.0739859 ,-0.0417835 ,-0.022119 ,0.00308797 ,0.0197657 ,-0.0137871 ,-0.036848 ,-0.0643794 ,-0.0514949 ,-0.0317105 ,0.016539 ,0.0308435 ,-0.00056883 ,-0.00756813 ,-0.00760612 ,-0.0234571 ,-0.0980915 ,-0.101175 ,-0.102354 ,-0.0920337 ,-0.100337 ,-0.0887628 ,0.0660931 ,0.228999 ,0.260675 ,0.266301 ,0.267907 ,0.281668 ,0.194433 ,0.132954 ,0.20707 ,0.220466 ,0.20936 ,0.191441 ,0.191441 ,0.191441 ,0.191441 ,0.191441 ,0.191441};


  int i60 = fabs(eta)*20.;
  if(i60<0)  i60 = 0;
  if(i60>59) i60 = 59; 

  if(icalo==1) {
//
// Girder parametrization
//
    double x = phi+ 3.1416-3.1416/32.*int((3.1416+phi)/(3.1416/32.));
    double scale = 0.;
    double pi = acos(-1.);
    if(x>pi/64.) x = pi/32.-x;

     if(x<0.005) {
       scale = X0CaloGirder[0]*(1-x/0.005)+X0CaloGirder[1] + X0CaloGirder[3];
     } else if(x<0.017) {
       scale = X0CaloGirder[1] + X0CaloGirder[3];
     } else if(x<0.028) {
       scale = X0CaloGirder[2] + X0CaloGirder[3];
     } else {
       scale = X0CaloGirder[3];
     }

     if(fabs(eta)>1.3) scale = 1.;
//
// eta dependence of X0 
//
    scale *= X0CaloScale[i60];
    X0Scale = scale;
//
// eta dependence of Eloss
//    
    int i30 = fabs(eta)*10.;
    if(i30<0)   i30 = 0;
    if(i30>29)  i30 = 29;

    double nfactor = 0.987363/1.05471;

    ElossScale = nfactor*ElossCaloScale[i30]*scale;

  } else {
//
// Muon system 
//
// eta dependence of X0 
//
    double scale = 1. + X0MuonScale[i60];
//
// Muon scale is now 1 with MuonTrackingGeometry and TrkDetDescrGeoModelCnv fixes
//
    scale = 1.0;
    X0Scale = scale;
    ElossScale = 0.93*scale;
539
  }
540
541
542
543

}