Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/parameterisations/gflash/src/GFlashSamplingShowerParameterisation.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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 // ------------------------------------------------------------
 29 // GEANT 4 class implementation
 30 //
 31 //      ------- GFlashSamplingShowerParameterisation -------
 32 //
 33 // Authors: E.Barberio & Joanna Weng - 11.2005
 34 // ------------------------------------------------------------
 35 
 36 #include <cmath>
 37 
 38 #include "GFlashSamplingShowerParameterisation.hh"
 39 #include "GVFlashShowerParameterisation.hh"
 40 #include "G4SystemOfUnits.hh"
 41 #include "Randomize.hh"
 42 #include "G4ios.hh"
 43 #include "G4Material.hh"
 44 #include "G4MaterialTable.hh"
 45 
 46 GFlashSamplingShowerParameterisation::GFlashSamplingShowerParameterisation(
 47   G4Material* aMat1, G4Material* aMat2, G4double dd1, G4double dd2,
 48   GFlashSamplingShowerTuning* aPar)
 49   : GVFlashShowerParameterisation(),
 50     ParAveT2(0.),
 51     ParSigLogT1(0.),
 52     ParSigLogT2(0.),
 53     ParSigLogA1(0.),
 54     ParSigLogA2(0.),
 55     ParRho1(0.),
 56     ParRho2(0.),
 57     ParsAveA2(0.),
 58     AveLogAlphah(0.),
 59     AveLogTmaxh(0.),
 60     SigmaLogAlphah(0.),
 61     SigmaLogTmaxh(0.),
 62     Rhoh(0.),
 63     Alphah(0.),
 64     Tmaxh(0.),
 65     Betah(0.),
 66     AveLogAlpha(0.),
 67     AveLogTmax(0.),
 68     SigmaLogAlpha(0.),
 69     SigmaLogTmax(0.),
 70     Rho(0.),
 71     Alpha(0.),
 72     Tmax(0.),
 73     Beta(0.)
 74 {
 75   if (!aPar) {
 76     thePar = new GFlashSamplingShowerTuning;
 77   }
 78   else {
 79     thePar = aPar;
 80   }
 81 
 82   SetMaterial(aMat1, aMat2);
 83   d1 = dd1;
 84   d2 = dd2;
 85 
 86   // Longitudinal Coefficients for a homogenious calo
 87 
 88   // shower max
 89   ParAveT1 = thePar->ParAveT1();  // ln (ln y -0.812)
 90   ParAveA1 = thePar->ParAveA1();  // ln a (0.81 + (0.458 + 2.26/Z)ln y)
 91   ParAveA2 = thePar->ParAveA2();
 92   ParAveA3 = thePar->ParAveA3();
 93   // Variance of shower max sampling
 94   ParSigLogT1 =
 95     thePar
 96       ->ParsSigLogT1();  // Sigma T1 (-1.4 + 1.26 ln y)**-1 --> bug : these two lines were missing,
 97   ParSigLogT2 = thePar->ParsSigLogT2();  // leaving ParSigLogT1, ParSigLogT2 as 0.0
 98   // variance of 'alpha'
 99   ParSigLogA1 =
100     thePar
101       ->ParSigLogA1();  // Sigma a (-0.58 + 0.86 ln y)**-1 --> bug : these two lines were missing
102   ParSigLogA2 = thePar->ParSigLogA2();  // leaving ParSigLogA1 ParSigLogAé as 0.0
103   // correlation alpha%T
104   ParRho1 = thePar->ParRho1();  // Rho = 0.705 -0.023 ln y  --> bug : these two lines were missing,
105   ParRho2 = thePar->ParRho2();  // leaving ParRho1 and ParRho2 being 0.0
106   // Sampling
107   ParsAveT1 = thePar->ParsAveT1();  // T_sam = log(exp( log T_hom) + t1*Fs-1 + t2*(1-ehat));
108   ParsAveT2 = thePar->ParsAveT2();
109   ParsAveA1 = thePar->ParsAveA1();
110   // Variance of shower max sampling
111   ParsSigLogT1 = thePar->ParsSigLogT1();  // Sigma T1 (-2.5 + 1.25 ln y)**-1 --> bug ParSigLogT1()
112                                           // was called instead of ParsSigLogT1(); Same for T2.
113   ParsSigLogT2 = thePar->ParsSigLogT2();
114   // variance of 'alpha'
115   ParsSigLogA1 = thePar->ParsSigLogA1();  // Sigma a (-0.82 + 0.79 ln y)**-1 --> bug ParSigLogA1()
116                                           // was called instead of ParsSigLogA1(); Same for A2
117   ParsSigLogA2 = thePar->ParsSigLogA2();
118   // correlation alpha%T
119   ParsRho1 =
120     thePar->ParsRho1();  // Rho = 0.784 -0.023 ln y --> bug was using ParRho1() and ParRho2()
121   ParsRho2 = thePar->ParsRho2();
122 
123   // Radial Coefficients
124   // r_C (tau)= z_1 +z_2 tau
125   // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2))))
126   ParRC1 = thePar->ParRC1();  // z_1 = 0.0251 + 0.00319 ln E
127   ParRC2 = thePar->ParRC2();
128   ParRC3 = thePar->ParRC3();  // z_2 = 0.1162 + - 0.000381 Z
129   ParRC4 = thePar->ParRC4();
130 
131   ParWC1 = thePar->ParWC1();
132   ParWC2 = thePar->ParWC2();
133   ParWC3 = thePar->ParWC3();
134   ParWC4 = thePar->ParWC4();
135   ParWC5 = thePar->ParWC5();
136   ParWC6 = thePar->ParWC6();
137   ParRT1 = thePar->ParRT1();
138   ParRT2 = thePar->ParRT2();
139   ParRT3 = thePar->ParRT3();
140   ParRT4 = thePar->ParRT4();
141   ParRT5 = thePar->ParRT5();
142   ParRT6 = thePar->ParRT6();
143 
144   // additional sampling parameter
145   ParsRC1 = thePar->ParsRC1();
146   ParsRC2 = thePar->ParsRC2();
147   ParsWC1 = thePar->ParsWC1();
148   ParsWC2 = thePar->ParsWC2();
149   ParsRT1 = thePar->ParsRT1();
150   ParsRT2 = thePar->ParsRT2();
151 
152   // Coeff for fluctuedted radial  profiles for a sampling media
153   ParsSpotT1 = thePar->ParSpotT1();  // T_spot = T_hom =(0.698 + 0.00212)
154   ParsSpotT2 = thePar->ParSpotT2();
155   ParsSpotA1 = thePar->ParSpotA1();  // a_spot= a_hom (0.639 + 0.00334)
156   ParsSpotA2 = thePar->ParSpotA2();
157   ParsSpotN1 = thePar->ParSpotN1();  // N_Spot 93 * ln(Z) E ** 0.876
158   ParsSpotN2 = thePar->ParSpotN2();
159   SamplingResolution = thePar->SamplingResolution();
160   ConstantResolution = thePar->ConstantResolution();
161   NoiseResolution = thePar->NoiseResolution();
162 
163   // Inits
164   NSpot = 0.00;
165   AlphaNSpot = 0.00;
166   TNSpot = 0.00;
167   BetaNSpot = 0.00;
168   RadiusCore = 0.00;
169   WeightCore = 0.00;
170   RadiusTail = 0.00;
171   ComputeZAX0EFFetc();
172 
173   G4cout << "/********************************************/ " << G4endl;
174   G4cout << "  - GFlashSamplingShowerParameterisation::Constructor -  " << G4endl;
175   G4cout << "/********************************************/ " << G4endl;
176 }
177 
178 // ------------------------------------------------------------
179 
180 GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation()
181 {
182    delete thePar;
183 }
184 
185 // ------------------------------------------------------------
186 
187 void GFlashSamplingShowerParameterisation::SetMaterial(G4Material* mat1, G4Material* mat2)
188 {
189   const G4double Es = 21.2 * MeV;
190   material1 = mat1;
191   Z1 = GetEffZ(material1);
192   A1 = GetEffA(material1);
193   density1 = material1->GetDensity();
194   X01 = material1->GetRadlen();
195   Ec1 = 2.66 * std::pow((X01 * Z1 / A1), 1.1);
196   // Ec1 = 610.0 * MeV / (Z1 + 1.24);
197   Rm1 = X01 * Es / Ec1;
198 
199   material2 = mat2;
200   Z2 = GetEffZ(material2);
201   A2 = GetEffA(material2);
202   density2 = material2->GetDensity();
203   X02 = material2->GetRadlen();
204   Ec2 = 2.66 * std::pow((X02 * Z2 / A2), 1.1);
205   // Ec2 = 610.0 * MeV / (Z2 + 1.24);
206   Rm2 = X02 * Es / Ec2;
207   // PrintMaterial();
208 }
209 
210 // ------------------------------------------------------------
211 
212 void GFlashSamplingShowerParameterisation::ComputeZAX0EFFetc()
213 {
214   G4cout << "/************ ComputeZAX0EFFetc ************/" << G4endl;
215   G4cout << "  - GFlashSamplingShowerParameterisation::Material -  " << G4endl;
216 
217   const G4double Es = 21.2 * MeV;
218 
219   // material and geometry parameters for a sampling calorimeter
220   G4double denominator = (d1 * density1 + d2 * density2);
221   G4double W1 = (d1 * density1) / denominator;
222   G4double W2 = (d2 * density2) / denominator;
223   Zeff = (W1 * Z1) + (W2 * Z2);  // X0*Es/Ec;
224   Aeff = (W1 * A1) + (W2 * A2);
225   Rhoeff = ((d1 * density1) + (d2 * density2)) / (d1 + d2);  // --> was G4double ( d2  + d1 );
226   X0eff = (W1 * Rhoeff) / (X01 * density1) + (W2 * Rhoeff) / (X02 * density2);
227   X0eff = 1. / X0eff;
228   Rmeff = 1. / ((((W1 * Ec1) / X01) + ((W2 * Ec2) / X02)) / Es);
229   Eceff = X0eff * ((W1 * Ec1) / X01 + (W2 * Ec2) / X02);
230   Fs = X0eff / (d1 + d2);  // --> was G4double ((d1/mm )+(d2/mm) ); Can't understand if dividing by
231                            // mm makes sense... looks weird.
232   ehat = (1. / (1 + 0.007 * (Z1 - Z2)));
233 
234   G4cout << "W1= " << W1 << G4endl;
235   G4cout << "W2= " << W2 << G4endl;
236   G4cout << "effective quantities Zeff = " << Zeff << G4endl;
237   G4cout << "effective quantities Aeff = " << Aeff << G4endl;
238   G4cout << "effective quantities Rhoeff = " << Rhoeff / g * cm3 << " g/cm3" << G4endl;
239   G4cout << "effective quantities X0eff = " << X0eff / cm << " cm" << G4endl;
240 
241   X0eff = X0eff * Rhoeff;
242 
243   G4cout << "effective quantities X0eff = " << X0eff / g * cm2 << " g/cm2" << G4endl;
244   X0eff = X0eff / Rhoeff;
245   G4cout << "effective quantities RMeff = " << Rmeff / cm << "  cm" << G4endl;
246   Rmeff = Rmeff * Rhoeff;
247   G4cout << "effective quantities RMeff = " << Rmeff / g * cm2 << " g/cm2" << G4endl;
248   Rmeff = Rmeff / Rhoeff;
249   G4cout << "effective quantities Eceff = " << Eceff / MeV << " MeV" << G4endl;
250   G4cout << "effective quantities Fs = " << Fs << G4endl;
251   G4cout << "effective quantities ehat = " << ehat << G4endl;
252   G4cout << "/********************************************/ " << G4endl;
253 }
254 
255 // ------------------------------------------------------------
256 
257 void GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile(G4double Energy)
258 {
259   if ((material1 == 0) || (material2 == 0)) {
260     G4Exception("GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile()",
261                 "InvalidSetup", FatalException, "No material initialized!");
262   }
263   G4double y = Energy / Eceff;
264   ComputeLongitudinalParameters(y);
265   GenerateEnergyProfile(y);
266   GenerateNSpotProfile(y);
267 }
268 
269 // ------------------------------------------------------------
270 
271 void GFlashSamplingShowerParameterisation::ComputeLongitudinalParameters(G4double y)
272 {
273   AveLogTmaxh = std::log(std::max(ParAveT1 + std::log(y), 0.1));  // ok
274   AveLogAlphah =
275     std::log(std::max(ParAveA1 + (ParAveA2 + ParAveA3 / Zeff) * std::log(y), 0.1));  // ok
276   // hom
277   SigmaLogTmaxh = std::min(0.5, 1.00 / (ParSigLogT1 + ParSigLogT2 * std::log(y)));  // ok
278   SigmaLogAlphah = std::min(0.5, 1.00 / (ParSigLogA1 + ParSigLogA2 * std::log(y)));  // ok
279   Rhoh = ParRho1 + ParRho2 * std::log(y);  // ok
280   // if sampling
281   AveLogTmax =
282     std::max(0.1, std::log(std::exp(AveLogTmaxh) + ParsAveT1 / Fs + ParsAveT2 * (1 - ehat)));  // ok
283   AveLogAlpha = std::max(0.1, std::log(std::exp(AveLogAlphah) + ParsAveA1 / Fs));  // ok
284   //
285   SigmaLogTmax = std::min(0.5, 1.00 / (ParsSigLogT1 + ParsSigLogT2 * std::log(y)));  // ok
286   SigmaLogAlpha = std::min(0.5, 1.00 / (ParsSigLogA1 + ParsSigLogA2 * std::log(y)));  // ok
287   Rho = ParsRho1 + ParsRho2 * std::log(y);  // ok
288 
289   if (0) {
290     G4cout << " y            = " << y << G4endl;
291     G4cout << " std::log(std::exp(AveLogTmaxh)  + ParsAveT1/Fs + ParsAveT2*(1-ehat)) = "
292            << " std::log(" << std::exp(AveLogTmaxh) << " + " << ParsAveT1 / Fs << " + "
293            << ParsAveT2 * (1 - ehat) << ") = "
294            << " std::log(" << std::exp(AveLogTmaxh) << " + " << ParsAveT1 << "/" << Fs << " + "
295            << ParsAveT2 << "*" << (1 - ehat) << ") = "
296            << " std::log(" << std::exp(AveLogTmaxh) + ParsAveT1 / Fs + ParsAveT2 * (1 - ehat) << ")"
297            << G4endl;
298     G4cout << " AveLogTmaxh    " << AveLogTmaxh << G4endl;
299     G4cout << " AveLogAlphah   " << AveLogAlphah << G4endl;
300     G4cout << " SigmaLogTmaxh  " << SigmaLogTmaxh << G4endl;
301     G4cout << " 1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)  ) = " << 1.00 << "/"
302            << (ParSigLogT1 + ParSigLogT2 * std::log(y)) << " = " << 1.00 << "/" << "("
303            << ParSigLogT1 << " + " << ParSigLogT2 * std::log(y) << " )  = " << 1.00 << "/" << "("
304            << ParSigLogT1 << " + " << ParSigLogT2 << "*" << std::log(y) << " ) " << G4endl;
305     G4cout << " SigmaLogAlphah " << SigmaLogAlphah << G4endl;
306     G4cout << " Rhoh           " << Rhoh << G4endl;
307     G4cout << " AveLogTmax     " << AveLogTmax << G4endl;
308     G4cout << " AveLogAlpha    " << AveLogAlpha << G4endl;
309     G4cout << " SigmaLogTmax   " << SigmaLogTmax << G4endl;
310     G4cout << " SigmaLogAlpha  " << SigmaLogAlpha << G4endl;
311     G4cout << " Rho            " << Rho << G4endl;
312   }
313 }
314 
315 // ------------------------------------------------------------
316 
317 void GFlashSamplingShowerParameterisation::GenerateEnergyProfile(G4double /* y */)
318 {
319   G4double Correlation1 = std::sqrt((1 + Rho) / 2);
320   G4double Correlation2 = std::sqrt((1 - Rho) / 2);
321   G4double Correlation1h = std::sqrt((1 + Rhoh) / 2);
322   G4double Correlation2h = std::sqrt((1 - Rhoh) / 2);
323   G4double Random1 = G4RandGauss::shoot();
324   G4double Random2 = G4RandGauss::shoot();
325 
326   Tmax = std::max(
327     1., std::exp(AveLogTmax + SigmaLogTmax * (Correlation1 * Random1 + Correlation2 * Random2)));
328   Alpha = std::max(
329     1.1, std::exp(AveLogAlpha + SigmaLogAlpha * (Correlation1 * Random1 - Correlation2 * Random2)));
330   Beta = (Alpha - 1.00) / Tmax;
331   // Parameters for Enenrgy Profile including correaltion and sigmas
332   Tmaxh =
333     std::exp(AveLogTmaxh + SigmaLogTmaxh * (Correlation1h * Random1 + Correlation2h * Random2));
334   Alphah =
335     std::exp(AveLogAlphah + SigmaLogAlphah * (Correlation1h * Random1 - Correlation2h * Random2));
336   Betah = (Alphah - 1.00) / Tmaxh;
337 }
338 
339 // ------------------------------------------------------------
340 
341 void GFlashSamplingShowerParameterisation::GenerateNSpotProfile(const G4double y)
342 {
343   TNSpot = Tmaxh * (ParsSpotT1 + ParsSpotT2 * Zeff);  // ok.
344   TNSpot = std::max(0.5, Tmaxh * (ParsSpotT1 + ParsSpotT2 * Zeff));
345   AlphaNSpot = Alphah * (ParsSpotA1 + ParsSpotA2 * Zeff);
346   BetaNSpot = (AlphaNSpot - 1.00) / TNSpot;  // ok
347   NSpot = ParsSpotN1 / SamplingResolution * std::pow(y * Eceff / GeV, ParsSpotN2);
348 }
349 
350 // ------------------------------------------------------------
351 
352 G4double GFlashSamplingShowerParameterisation::ApplySampling(const G4double DEne, const G4double)
353 {
354   G4double DEneFluctuated = DEne;
355   G4double Resolution = std::pow(SamplingResolution, 2);
356 
357   //       +pow(NoiseResolution,2)/  //@@@@@@@@ FIXME
358   //                         Energy*(1.*MeV)+
359   //                         pow(ConstantResolution,2)*
360   //                          Energy/(1.*MeV);
361 
362   if (Resolution > 0.0 && DEne > 0.00) {
363     // G4float x1 = DEne / Resolution;
364     // G4float x2 = G4RandGamma::shoot(x1, 1.0) * Resolution;
365     // DEneFluctuated = x2;
366     G4double x1 = DEne / Resolution;
367     G4double x2 = 1.0 / Resolution;
368     DEneFluctuated = G4RandGamma::shoot(x1, x2);
369   }
370   return DEneFluctuated;
371 }
372 
373 // ------------------------------------------------------------
374 
375 G4double GFlashSamplingShowerParameterisation::
376 IntegrateEneLongitudinal(G4double LongitudinalStep)
377 {
378   G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
379   G4float x1= Betah*LongitudinalStepInX0;
380   G4float x2= Alphah;
381   float x3 =  gam(x1,x2);
382   G4double DEne=x3;
383   return DEne;
384 }
385 
386 // ------------------------------------------------------------
387 
388 G4double GFlashSamplingShowerParameterisation::IntegrateNspLongitudinal(G4double LongitudinalStep)
389 {
390   G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
391   G4float x1 = BetaNSpot * LongitudinalStepInX0;
392   G4float x2 = AlphaNSpot;
393   G4float x3 = gam(x1, x2);
394   G4double DNsp = x3;
395   return DNsp;
396 }
397 
398 // ------------------------------------------------------------
399 
400 G4double GFlashSamplingShowerParameterisation::GenerateRadius(G4int ispot, G4double Energy,
401                                                               G4double LongitudinalPosition)
402 {
403   if (ispot < 1) {
404     // Determine lateral parameters in the middle of the step.
405     // They depend on energy & position along step
406     //
407     G4double Tau = ComputeTau(LongitudinalPosition);
408     ComputeRadialParameters(Energy, Tau);
409   }
410 
411   G4double Radius;
412   G4double Random1 = G4UniformRand();
413   G4double Random2 = G4UniformRand();
414   if (Random1 < WeightCore)  // WeightCore = p < w_i
415   {
416     Radius = Rmeff * RadiusCore * std::sqrt(Random2 / (1. - Random2));
417   }
418   else {
419     Radius = Rmeff * RadiusTail * std::sqrt(Random2 / (1. - Random2));
420   }
421   Radius = std::min(Radius, DBL_MAX);
422   return Radius;
423 }
424 
425 // ------------------------------------------------------------
426 
427 G4double GFlashSamplingShowerParameterisation::ComputeTau(G4double LongitudinalPosition)
428 {
429   G4double tau = LongitudinalPosition / Tmax / X0eff  //<t> = T* a /(a - 1)
430                  * (Alpha - 1.00) / Alpha * std::exp(AveLogAlpha)
431                  / (std::exp(AveLogAlpha) - 1.);  // ok
432   return tau;
433 }
434 
435 // ------------------------------------------------------------
436 
437 void GFlashSamplingShowerParameterisation::ComputeRadialParameters(G4double Energy, G4double Tau)
438 {
439   G4double z1 = ParRC1 + ParRC2 * std::log(Energy / GeV);  // ok
440   G4double z2 = ParRC3 + ParRC4 * Zeff;  // ok
441   RadiusCore = z1 + z2 * Tau;  // ok
442   G4double p1 = ParWC1 + ParWC2 * Zeff;  // ok
443   G4double p2 = ParWC3 + ParWC4 * Zeff;  // ok
444   G4double p3 = ParWC5 + ParWC6 * std::log(Energy / GeV);  // ok
445   WeightCore = p1 * std::exp((p2 - Tau) / p3 - std::exp((p2 - Tau) / p3));  // ok
446 
447   G4double k1 = ParRT1 + ParRT2 * Zeff;  // ok
448   G4double k2 = ParRT3;  // ok
449   G4double k3 = ParRT4;  // ok
450   G4double k4 = ParRT5 + ParRT6 * std::log(Energy / GeV);  // ok
451 
452   RadiusTail = k1 * (std::exp(k3 * (Tau - k2)) + std::exp(k4 * (Tau - k2)));  // ok
453 
454   // sampling calorimeter
455 
456   RadiusCore = RadiusCore + ParsRC1 * (1 - ehat) + ParsRC2 / Fs * std::exp(-Tau);  // ok
457   WeightCore =
458     WeightCore + (1 - ehat) * (ParsWC1 + ParsWC2 / Fs * std::exp(-std::pow((Tau - 1.), 2)));  // ok
459   RadiusTail = RadiusTail + (1 - ehat) * ParsRT1 + ParsRT2 / Fs * std::exp(-Tau);  // ok
460 }
461 
462 // ------------------------------------------------------------
463 
464 G4double GFlashSamplingShowerParameterisation::
465 GenerateExponential(const G4double /* Energy */ )
466 {
467   G4double ParExp1 =  9./7.*X0eff;
468   G4double random  = -ParExp1*G4RandExponential::shoot() ;
469   return random;
470 }
471