Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/parameterisations/gflash/src/GFlashHomoShowerParameterisation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /parameterisations/gflash/src/GFlashHomoShowerParameterisation.cc (Version 11.3.0) and /parameterisations/gflash/src/GFlashHomoShowerParameterisation.cc (Version 10.1.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: GFlashHomoShowerParameterisation.cc 69579 2013-05-08 13:53:57Z gcosmo $
 26 //                                                 27 //
 27 //                                                 28 //
 28 // -------------------------------------------     29 // ------------------------------------------------------------
 29 // GEANT 4 class implementation                    30 // GEANT 4 class implementation
 30 //                                                 31 //
 31 //      ------- GFlashHomoShowerParameterisati     32 //      ------- GFlashHomoShowerParameterisation -------
 32 //                                                 33 //
 33 // Authors: E.Barberio & Joanna Weng - 9.11.20     34 // Authors: E.Barberio & Joanna Weng - 9.11.2004
 34 // -------------------------------------------     35 // ------------------------------------------------------------
 35                                                    36 
 36 #include <cmath>                                   37 #include <cmath>
 37                                                    38 
 38 #include "GFlashHomoShowerParameterisation.hh"     39 #include "GFlashHomoShowerParameterisation.hh"
 39 #include "GVFlashShowerParameterisation.hh"        40 #include "GVFlashShowerParameterisation.hh"
 40 #include "G4PhysicalConstants.hh"                  41 #include "G4PhysicalConstants.hh"
 41 #include "G4SystemOfUnits.hh"                      42 #include "G4SystemOfUnits.hh"
 42 #include "Randomize.hh"                            43 #include "Randomize.hh"
 43 #include "G4ios.hh"                                44 #include "G4ios.hh"
 44 #include "G4Material.hh"                           45 #include "G4Material.hh"
 45 #include "G4MaterialTable.hh"                      46 #include "G4MaterialTable.hh"
 46                                                    47 
 47 GFlashHomoShowerParameterisation::GFlashHomoSh <<  48 GFlashHomoShowerParameterisation::
 48                                                <<  49 GFlashHomoShowerParameterisation(G4Material * aMat,
                                                   >>  50                                  GVFlashHomoShowerTuning * aPar)
 49   : GVFlashShowerParameterisation(),               51   : GVFlashShowerParameterisation(),
 50     ConstantResolution(0.),                    <<  52     ConstantResolution(0.), NoiseResolution(0.), SamplingResolution(0.),
 51     NoiseResolution(0.),                       <<  53     AveLogAlphah(0.), AveLogTmaxh(0.), SigmaLogAlphah(0.), SigmaLogTmaxh(0.),
 52     SamplingResolution(0.),                    <<  54     Rhoh(0.), Alphah(0.), Tmaxh(0.), Betah(0.)
 53     AveLogAlphah(0.),                          <<  55 
 54     AveLogTmaxh(0.),                           <<  56 {  
 55     SigmaLogAlphah(0.),                        <<  57   if(!aPar) { thePar = new GVFlashHomoShowerTuning; owning = true; }
 56     SigmaLogTmaxh(0.),                         <<  58   else      { thePar = aPar; owning = false; }
 57     Rhoh(0.),                                  << 
 58     Alphah(0.),                                << 
 59     Tmaxh(0.),                                 << 
 60     Betah(0.)                                  << 
 61                                                << 
 62 {                                              << 
 63   if (!aPar) {                                 << 
 64     thePar = new GVFlashHomoShowerTuning;      << 
 65   }                                            << 
 66   else {                                       << 
 67     thePar = aPar;                             << 
 68   }                                            << 
 69                                                    59 
 70   SetMaterial(aMat);                               60   SetMaterial(aMat);
 71   PrintMaterial(aMat);                             61   PrintMaterial(aMat);
 72                                                    62 
 73   /*******************************************     63   /********************************************/
 74   /* Homo Calorimeter                              64   /* Homo Calorimeter                         */
 75   /******************************************* <<  65   /********************************************/ 
 76   // Longitudinal Coefficients for a homogenio     66   // Longitudinal Coefficients for a homogenious calo
 77   // shower max                                    67   // shower max
 78   //                                               68   //
 79   ParAveT1 = thePar->ParAveT1();  // ln (ln y  <<  69   ParAveT1    = thePar->ParAveT1();   // ln (ln y -0.812)  
 80   ParAveA1 = thePar->ParAveA1();  // ln a (0.8 <<  70   ParAveA1    = thePar->ParAveA1();   // ln a (0.81 + (0.458 + 2.26/Z)ln y)
 81   ParAveA2 = thePar->ParAveA2();               <<  71   ParAveA2    = thePar->ParAveA2();
 82   ParAveA3 = thePar->ParAveA3();               <<  72   ParAveA3    = thePar->ParAveA3();
 83                                                    73 
 84   // Variance of shower max                        74   // Variance of shower max
 85   ParSigLogT1 = thePar->ParSigLogT1();  // Sig <<  75   ParSigLogT1 = thePar->ParSigLogT1();     // Sigma T1 (-1.4 + 1.26 ln y)**-1 
 86   ParSigLogT2 = thePar->ParSigLogT2();             76   ParSigLogT2 = thePar->ParSigLogT2();
 87                                                    77 
 88   // variance of 'alpha'                           78   // variance of 'alpha'
 89   //                                               79   //
 90   ParSigLogA1 = thePar->ParSigLogA1();  // Sig <<  80   ParSigLogA1 = thePar->ParSigLogA1();    // Sigma a (-0.58 + 0.86 ln y)**-1 
 91   ParSigLogA2 = thePar->ParSigLogA2();             81   ParSigLogA2 = thePar->ParSigLogA2();
 92                                                    82 
 93   // correlation alpha%T                           83   // correlation alpha%T
 94   //                                               84   //
 95   ParRho1 = thePar->ParRho1();  // Rho = 0.705 <<  85   ParRho1     = thePar->ParRho1();   // Rho = 0.705 -0.023 ln y 
 96   ParRho2 = thePar->ParRho2();                 <<  86   ParRho2     = thePar->ParRho2();
 97                                                    87 
 98   // Radial Coefficients                           88   // Radial Coefficients
 99   // r_C (tau)= z_1 +z_2 tau                       89   // r_C (tau)= z_1 +z_2 tau
100   // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+st     90   // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2))))
101   //                                               91   //
102   ParRC1 = thePar->ParRC1();  // z_1 = 0.0251  <<  92   ParRC1 =   thePar->ParRC1();   // z_1 = 0.0251 + 0.00319 ln E
103   ParRC2 = thePar->ParRC2();                   <<  93   ParRC2 =   thePar->ParRC2();
104                                                    94 
105   ParRC3 = thePar->ParRC3();  // z_2 = 0.1162  <<  95   ParRC3 =   thePar->ParRC3();   // z_2 = 0.1162 + - 0.000381 Z
106   ParRC4 = thePar->ParRC4();                   <<  96   ParRC4 =   thePar->ParRC4();
107                                                    97 
108   ParWC1 = thePar->ParWC1();                       98   ParWC1 = thePar->ParWC1();
109   ParWC2 = thePar->ParWC2();                       99   ParWC2 = thePar->ParWC2();
110   ParWC3 = thePar->ParWC3();                      100   ParWC3 = thePar->ParWC3();
111   ParWC4 = thePar->ParWC4();                      101   ParWC4 = thePar->ParWC4();
112   ParWC5 = thePar->ParWC5();                   << 102   ParWC5 = thePar->ParWC5(); 
113   ParWC6 = thePar->ParWC6();                      103   ParWC6 = thePar->ParWC6();
114                                                   104 
115   ParRT1 = thePar->ParRT1();                      105   ParRT1 = thePar->ParRT1();
116   ParRT2 = thePar->ParRT2();                      106   ParRT2 = thePar->ParRT2();
117   ParRT3 = thePar->ParRT3();                      107   ParRT3 = thePar->ParRT3();
118   ParRT4 = thePar->ParRT4();                   << 108   ParRT4 = thePar->ParRT4(); 
119   ParRT5 = thePar->ParRT5();                      109   ParRT5 = thePar->ParRT5();
120   ParRT6 = thePar->ParRT6();                      110   ParRT6 = thePar->ParRT6();
121                                                   111 
122   // Coeff for fluctueted radial  profiles for    112   // Coeff for fluctueted radial  profiles for a uniform media
123   //                                              113   //
124   ParSpotT1 = thePar->ParSpotT1();  // T_spot  << 114   ParSpotT1   = thePar->ParSpotT1();     // T_spot = T_hom =(0.698 + 0.00212)
125   ParSpotT2 = thePar->ParSpotT2();             << 115   ParSpotT2   = thePar->ParSpotT2();
126                                                   116 
127   ParSpotA1 = thePar->ParSpotA1();  // a_spot= << 117   ParSpotA1   = thePar->ParSpotA1();     // a_spot= a_hom (0.639 + 0.00334)
128   ParSpotA2 = thePar->ParSpotA2();             << 118   ParSpotA2   = thePar->ParSpotA2();  
129                                                   119 
130   ParSpotN1 = thePar->ParSpotN1();  // N_Spot  << 120   ParSpotN1   = thePar->ParSpotN1();     // N_Spot 93 * ln(Z) E ** 0.876   
131   ParSpotN2 = thePar->ParSpotN2();             << 121   ParSpotN2   = thePar->ParSpotN2(); 
132                                                   122 
133   // Inits                                        123   // Inits
134                                                   124 
135   NSpot = 0.00;                                << 125   NSpot         = 0.00;
136   AlphaNSpot = 0.00;                           << 126   AlphaNSpot    = 0.00;
137   TNSpot = 0.00;                               << 127   TNSpot        = 0.00;
138   BetaNSpot = 0.00;                            << 128   BetaNSpot     = 0.00;
139                                                << 129 
140   RadiusCore = 0.00;                           << 130   RadiusCore    = 0.00;
141   WeightCore = 0.00;                           << 131   WeightCore    = 0.00;
142   RadiusTail = 0.00;                           << 132   RadiusTail    = 0.00; 
143                                                   133 
144   G4cout << "/********************************    134   G4cout << "/********************************************/ " << G4endl;
145   G4cout << "  - GFlashHomoShowerParameterisat << 135   G4cout << "  - GFlashHomoShowerParameterisation::Constructor -  " << G4endl;  
146   G4cout << "/********************************    136   G4cout << "/********************************************/ " << G4endl;
147 }                                                 137 }
148                                                   138 
149 void GFlashHomoShowerParameterisation::SetMate << 139 void GFlashHomoShowerParameterisation::SetMaterial(G4Material *mat)
150 {                                                 140 {
151   material = mat;                              << 141   material= mat;
152   Z = GetEffZ(material);                          142   Z = GetEffZ(material);
153   A = GetEffA(material);                          143   A = GetEffA(material);
154   density = material->GetDensity() / (g / cm3) << 144   density = material->GetDensity()/(g/cm3);
155   X0 = material->GetRadlen();                  << 145   X0  = material->GetRadlen(); 
156   // O. I. Dovzhenkko and A. A. Pommanskii     << 146   Ec      = 2.66 * std::pow((X0 * Z / A),1.1); 
157   Ec = 2.66 * std::pow((X0 * Z / A), 1.1);     << 147   G4double Es = 21*MeV;
158   // // Rossi appriximation                    << 148   Rm = X0*Es/Ec;
159   // Ec = 610.0 * MeV / (Z + 1.24);            << 149   // PrintMaterial(); 
160   const G4double Es = 21.2 * MeV;              << 
161   Rm = X0 * Es / Ec;                           << 
162   // PrintMaterial();                          << 
163 }                                                 150 }
164                                                   151 
165 GFlashHomoShowerParameterisation::~GFlashHomoS    152 GFlashHomoShowerParameterisation::~GFlashHomoShowerParameterisation()
166 {                                                 153 {
167   delete thePar;                               << 154   if(owning) { delete thePar; }
168 }                                                 155 }
169                                                   156 
170 void GFlashHomoShowerParameterisation::Generat << 157 void GFlashHomoShowerParameterisation::
                                                   >> 158 GenerateLongitudinalProfile(G4double Energy)
171 {                                                 159 {
172   if (material == 0) {                         << 160   if (material==0) 
173     G4Exception("GFlashHomoShowerParameterisat << 161   {
174                 FatalException, "No material i << 162     G4Exception("GFlashHomoShowerParameterisation::GenerateLongitudinalProfile()",
                                                   >> 163                 "InvalidSetup", FatalException, "No material initialized!");
175   }                                               164   }
176                                                << 165   
177   G4double y = Energy / Ec;                    << 166   G4double y = Energy/Ec;
178   ComputeLongitudinalParameters(y);            << 167   ComputeLongitudinalParameters(y);  
179   GenerateEnergyProfile(y);                       168   GenerateEnergyProfile(y);
180   GenerateNSpotProfile(y);                        169   GenerateNSpotProfile(y);
181 }                                                 170 }
182                                                   171 
183 void GFlashHomoShowerParameterisation::Compute << 172 void
                                                   >> 173 GFlashHomoShowerParameterisation::ComputeLongitudinalParameters(G4double y)
184 {                                                 174 {
185   AveLogTmaxh = std::log(ParAveT1 + std::log(y << 175   AveLogTmaxh  = std::log(ParAveT1 + std::log(y));
186   // ok  <ln T hom>                            << 176     //ok  <ln T hom>
187   AveLogAlphah = std::log(ParAveA1 + (ParAveA2 << 177   AveLogAlphah = std::log(ParAveA1 + (ParAveA2+ParAveA3/Z)*std::log(y));
188   // ok  <ln alpha hom>                        << 178     //ok  <ln alpha hom> 
189                                                << 179 
190   SigmaLogTmaxh = 1.00 / (ParSigLogT1 + ParSig << 180   SigmaLogTmaxh  = 1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ;
191   // ok sigma (ln T hom)                       << 181     //ok sigma (ln T hom)
192   SigmaLogAlphah = 1.00 / (ParSigLogA1 + ParSi << 182   SigmaLogAlphah = 1.00/( ParSigLogA1 + ParSigLogA2*std::log(y));
193   // ok sigma (ln alpha hom)                   << 183     //ok sigma (ln alpha hom)
194   Rhoh = ParRho1 + ParRho2 * std::log(y);  //  << 184   Rhoh           = ParRho1+ParRho2*std::log(y);        //ok             
195 }                                                 185 }
196                                                   186 
197 void GFlashHomoShowerParameterisation::Generat    187 void GFlashHomoShowerParameterisation::GenerateEnergyProfile(G4double /* y */)
198 {                                              << 188 { 
199   G4double Correlation1h = std::sqrt((1 + Rhoh << 189   G4double Correlation1h = std::sqrt((1+Rhoh)/2);
200   G4double Correlation2h = std::sqrt((1 - Rhoh << 190   G4double Correlation2h = std::sqrt((1-Rhoh)/2);
201                                                   191 
202   G4double Random1 = G4RandGauss::shoot();        192   G4double Random1 = G4RandGauss::shoot();
203   G4double Random2 = G4RandGauss::shoot();        193   G4double Random2 = G4RandGauss::shoot();
204                                                   194 
205   // Parameters for Enenrgy Profile including  << 195   // Parameters for Enenrgy Profile including correaltion and sigmas 
206                                                << 196  
207   Tmaxh =                                      << 197   Tmaxh  = std::exp( AveLogTmaxh  + SigmaLogTmaxh  *
208     std::exp(AveLogTmaxh + SigmaLogTmaxh * (Co << 198            (Correlation1h*Random1 + Correlation2h*Random2) );
209   Alphah =                                     << 199   Alphah = std::exp( AveLogAlphah + SigmaLogAlphah *
210     std::exp(AveLogAlphah + SigmaLogAlphah * ( << 200            (Correlation1h*Random1 - Correlation2h*Random2) );
211   Betah = (Alphah - 1.00) / Tmaxh;             << 201   Betah  = (Alphah-1.00)/Tmaxh;
212 }                                                 202 }
213                                                   203 
214 void GFlashHomoShowerParameterisation::Generat    204 void GFlashHomoShowerParameterisation::GenerateNSpotProfile(const G4double y)
215 {                                                 205 {
216   TNSpot = Tmaxh * (ParSpotT1 + ParSpotT2 * Z) << 206   TNSpot     = Tmaxh *  (ParSpotT1+ParSpotT2*Z);   // ok
217   AlphaNSpot = Alphah * (ParSpotA1 + ParSpotA2 << 207   AlphaNSpot = Alphah * (ParSpotA1+ParSpotA2*Z);   
218   BetaNSpot = (AlphaNSpot - 1.00) / TNSpot;  / << 208   BetaNSpot  = (AlphaNSpot-1.00)/TNSpot;           // ok 
219   NSpot = ParSpotN1 * std::log(Z) * std::pow(( << 209   NSpot      = ParSpotN1 * std::log(Z)*std::pow((y*Ec)/GeV,ParSpotN2 ); // ok
220 }                                                 210 }
221                                                   211 
222 G4double GFlashHomoShowerParameterisation::       212 G4double GFlashHomoShowerParameterisation::
223 IntegrateEneLongitudinal(G4double Longitudinal    213 IntegrateEneLongitudinal(G4double LongitudinalStep)
224 {                                                 214 {
225   G4double LongitudinalStepInX0 = Longitudinal    215   G4double LongitudinalStepInX0 = LongitudinalStep / X0;
226   G4float x1= Betah*LongitudinalStepInX0;         216   G4float x1= Betah*LongitudinalStepInX0;
227   G4float x2= Alphah;                             217   G4float x2= Alphah;
228   float x3 =  gam(x1,x2);                         218   float x3 =  gam(x1,x2);
229   G4double DEne=x3;                               219   G4double DEne=x3;
230   return DEne;                                    220   return DEne;
231 }                                                 221 }
232                                                   222 
233 G4double GFlashHomoShowerParameterisation::Int << 223 G4double GFlashHomoShowerParameterisation::
                                                   >> 224 IntegrateNspLongitudinal(G4double LongitudinalStep)
234 {                                                 225 {
235   G4double LongitudinalStepInX0 = Longitudinal << 226   G4double LongitudinalStepInX0 = LongitudinalStep / X0; 
236   G4float x1 = BetaNSpot * LongitudinalStepInX << 227   G4float x1 = BetaNSpot*LongitudinalStepInX0;
237   G4float x2 = AlphaNSpot;                        228   G4float x2 = AlphaNSpot;
238   G4float x3 = gam(x1, x2);                    << 229   G4float x3 =  gam(x1,x2);
239   G4double DNsp = x3;                             230   G4double DNsp = x3;
240   return DNsp;                                    231   return DNsp;
241 }                                                 232 }
242                                                   233 
243 G4double GFlashHomoShowerParameterisation::Gen << 234 
244                                                << 235 G4double GFlashHomoShowerParameterisation::
                                                   >> 236 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition)
245 {                                                 237 {
246   if (ispot < 1) {                             << 238   if(ispot < 1) 
                                                   >> 239   {
247     // Determine lateral parameters in the mid    240     // Determine lateral parameters in the middle of the step.
248     // They depend on energy & position along     241     // They depend on energy & position along step.
249     //                                            242     //
250     G4double Tau = ComputeTau(LongitudinalPosi    243     G4double Tau = ComputeTau(LongitudinalPosition);
251     ComputeRadialParameters(Energy, Tau);      << 244     ComputeRadialParameters(Energy,Tau);  
252   }                                               245   }
253                                                   246 
254   G4double Radius;                                247   G4double Radius;
255   G4double Random1 = G4UniformRand();             248   G4double Random1 = G4UniformRand();
256   G4double Random2 = G4UniformRand();          << 249   G4double Random2 = G4UniformRand(); 
257                                                   250 
258   if (Random1 < WeightCore)  // WeightCore = p << 251   if(Random1  <WeightCore) //WeightCore = p < w_i  
259   {                                               252   {
260     Radius = Rm * RadiusCore * std::sqrt(Rando << 253     Radius = Rm * RadiusCore * std::sqrt( Random2/(1. - Random2) );
261   }                                            << 
262   else {                                       << 
263     Radius = Rm * RadiusTail * std::sqrt(Rando << 
264   }                                               254   }
265   Radius = std::min(Radius, DBL_MAX);          << 255   else
                                                   >> 256   {
                                                   >> 257     Radius = Rm * RadiusTail * std::sqrt( Random2/(1. - Random2) );
                                                   >> 258   }   
                                                   >> 259   Radius =  std::min(Radius,DBL_MAX);
266   return Radius;                                  260   return Radius;
267 }                                                 261 }
268                                                   262 
269 G4double GFlashHomoShowerParameterisation::Com << 263 G4double GFlashHomoShowerParameterisation::
                                                   >> 264 ComputeTau(G4double LongitudinalPosition)
270 {                                                 265 {
271   G4double tau = LongitudinalPosition / Tmaxh  << 266   G4double tau = LongitudinalPosition / Tmaxh / X0     //<t> = T* a /(a - 1) 
272                  * (Alphah - 1.00) / Alphah *  << 267   * (Alphah-1.00) /Alphah * 
273                  / (std::exp(AveLogAlphah) - 1 << 268   std::exp(AveLogAlphah)/(std::exp(AveLogAlphah)-1.);  //ok 
274   return tau;                                     269   return tau;
275 }                                                 270 }
276                                                   271 
277 void GFlashHomoShowerParameterisation::Compute << 272 void GFlashHomoShowerParameterisation::
                                                   >> 273 ComputeRadialParameters(G4double Energy, G4double Tau)
278 {                                                 274 {
279   G4double z1 = ParRC1 + ParRC2 * std::log(Ene << 275   G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV)  ;       //ok
280   G4double z2 = ParRC3 + ParRC4 * Z;  // ok    << 276   G4double z2 = ParRC3+ParRC4*Z ;                              //ok
281   RadiusCore = z1 + z2 * Tau;  // ok           << 277   RadiusCore   =  z1 + z2 * Tau  ;                             //ok 
282                                                << 278 
283   G4double p1 = ParWC1 + ParWC2 * Z;  // ok    << 279   G4double p1 = ParWC1+ParWC2*Z;                               //ok
284   G4double p2 = ParWC3 + ParWC4 * Z;  // ok    << 280   G4double p2 = ParWC3+ParWC4*Z;                               //ok
285   G4double p3 = ParWC5 + ParWC6 * std::log(Ene << 281   G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV);            //ok
286                                                << 282 
287   WeightCore = p1 * std::exp((p2 - Tau) / p3 - << 283   WeightCore   =  p1 * std::exp( (p2-Tau)/p3 - std::exp( (p2-Tau) /p3) ); //ok
288                                                << 284 
289   G4double k1 = ParRT1 + ParRT2 * Z;  // ok    << 285   G4double k1 = ParRT1+ParRT2*Z;                   // ok
290   G4double k2 = ParRT3;  // ok                 << 286   G4double k2 = ParRT3;                            // ok
291   G4double k3 = ParRT4;  // ok                 << 287   G4double k3 = ParRT4;                            // ok
292   G4double k4 = ParRT5 + ParRT6 * std::log(Ene << 288   G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV);    // ok
293                                                   289 
294   RadiusTail = k1 * (std::exp(k3 * (Tau - k2)) << 290   RadiusTail   = k1*(std::exp(k3*(Tau-k2))  +
                                                   >> 291   std::exp(k4*(Tau-k2)) );            //ok
295 }                                                 292 }
296                                                   293 
297 G4double GFlashHomoShowerParameterisation::       294 G4double GFlashHomoShowerParameterisation::
298 GenerateExponential(const G4double /* Energy *    295 GenerateExponential(const G4double /* Energy */ )
299 {                                                 296 {
300   G4double ParExp1 =  9./7.*X0;                   297   G4double ParExp1 =  9./7.*X0;
301   G4double random  = -ParExp1*G4RandExponentia    298   G4double random  = -ParExp1*G4RandExponential::shoot() ;
302   return random;                                  299   return random;
303 }                                                 300 }
304                                                   301