[823] | 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 | // $Id: GFlashSamplingShowerParameterisation.cc,v 1.6 2006/06/29 19:14:20 gunter Exp $ |
---|
[850] | 27 | // GEANT4 tag $Name: HEAD $ |
---|
[823] | 28 | // |
---|
| 29 | // |
---|
| 30 | // ------------------------------------------------------------ |
---|
| 31 | // GEANT 4 class implementation |
---|
| 32 | // |
---|
| 33 | // ------- GFlashSamplingShowerParameterisation ------- |
---|
| 34 | // |
---|
| 35 | // Authors: E.Barberio & Joanna Weng - 11.2005 |
---|
| 36 | // ------------------------------------------------------------ |
---|
| 37 | |
---|
| 38 | #include "GVFlashShowerParameterisation.hh" |
---|
| 39 | #include "GFlashSamplingShowerParameterisation.hh" |
---|
| 40 | #include <cmath> |
---|
| 41 | #include "Randomize.hh" |
---|
| 42 | #include "G4ios.hh" |
---|
| 43 | #include "G4Material.hh" |
---|
| 44 | #include "G4MaterialTable.hh" |
---|
| 45 | |
---|
| 46 | GFlashSamplingShowerParameterisation:: |
---|
| 47 | GFlashSamplingShowerParameterisation(G4Material* aMat1, G4Material* aMat2, |
---|
| 48 | G4double dd1, G4double dd2, |
---|
| 49 | GFlashSamplingShowerTuning* aPar) |
---|
| 50 | : GVFlashShowerParameterisation() |
---|
| 51 | { |
---|
| 52 | if(!aPar) { thePar = new GFlashSamplingShowerTuning; } |
---|
| 53 | else { thePar = aPar; } |
---|
| 54 | |
---|
| 55 | SetMaterial(aMat1,aMat2 ); |
---|
| 56 | d1=dd1; |
---|
| 57 | d2=dd2; |
---|
| 58 | |
---|
| 59 | // Longitudinal Coefficients for a homogenious calo |
---|
| 60 | |
---|
| 61 | // shower max |
---|
| 62 | ParAveT1 = thePar->ParAveT1(); // ln (ln y -0.812) |
---|
| 63 | ParAveA1 = thePar->ParAveA1(); // ln a (0.81 + (0.458 + 2.26/Z)ln y) |
---|
| 64 | ParAveA2 = thePar->ParAveA2(); |
---|
| 65 | ParAveA3 = thePar->ParAveA3(); |
---|
| 66 | // Sampling |
---|
| 67 | ParsAveT1 = thePar->ParsAveT1(); // T_sam = log(exp( log T_hom) + t1*Fs-1 + t2*(1-ehat)); |
---|
| 68 | ParsAveT2 = thePar->ParsAveT2(); |
---|
| 69 | ParsAveA1 = thePar->ParsAveA1(); |
---|
| 70 | // Variance of shower max sampling |
---|
| 71 | ParsSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-2.5 + 1.25 ln y)**-1 |
---|
| 72 | ParsSigLogT2 = thePar->ParSigLogT2(); |
---|
| 73 | // variance of 'alpha' |
---|
| 74 | ParsSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.82 + 0.79 ln y)**-1 |
---|
| 75 | ParsSigLogA2 = thePar->ParSigLogA2(); |
---|
| 76 | // correlation alpha%T |
---|
| 77 | ParsRho1 = thePar->ParRho1(); // Rho = 0.784 -0.023 ln y |
---|
| 78 | ParsRho2 = thePar->ParRho2(); |
---|
| 79 | |
---|
| 80 | // Radial Coefficients |
---|
| 81 | // r_C (tau)= z_1 +z_2 tau |
---|
| 82 | // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2)))) |
---|
| 83 | ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E |
---|
| 84 | ParRC2 = thePar->ParRC2(); |
---|
| 85 | ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z |
---|
| 86 | ParRC4 = thePar->ParRC4(); |
---|
| 87 | |
---|
| 88 | ParWC1 = thePar->ParWC1(); |
---|
| 89 | ParWC2 = thePar->ParWC2(); |
---|
| 90 | ParWC3 = thePar->ParWC3(); |
---|
| 91 | ParWC4 = thePar->ParWC4(); |
---|
| 92 | ParWC5 = thePar->ParWC5(); |
---|
| 93 | ParWC6 = thePar->ParWC6(); |
---|
| 94 | ParRT1 = thePar->ParRT1(); |
---|
| 95 | ParRT2 = thePar->ParRT2(); |
---|
| 96 | ParRT3 = thePar->ParRT3(); |
---|
| 97 | ParRT4 = thePar->ParRT4(); |
---|
| 98 | ParRT5 = thePar->ParRT5(); |
---|
| 99 | ParRT6 = thePar->ParRT6(); |
---|
| 100 | |
---|
| 101 | //additional sampling parameter |
---|
| 102 | ParsRC1= thePar->ParsRC1(); |
---|
| 103 | ParsRC2= thePar->ParsRC2(); |
---|
| 104 | ParsWC1= thePar->ParsWC1(); |
---|
| 105 | ParsWC2= thePar->ParsWC2(); |
---|
| 106 | ParsRT1= thePar->ParsRT1(); |
---|
| 107 | ParsRT2= thePar->ParsRT2(); |
---|
| 108 | |
---|
| 109 | // Coeff for fluctuedted radial profiles for a sampling media |
---|
| 110 | ParsSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212) |
---|
| 111 | ParsSpotT2 = thePar->ParSpotT2(); |
---|
| 112 | ParsSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334) |
---|
| 113 | ParsSpotA2 = thePar->ParSpotA2(); |
---|
| 114 | ParsSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876 |
---|
| 115 | ParsSpotN2 = thePar->ParSpotN2(); |
---|
| 116 | SamplingResolution = thePar->SamplingResolution(); |
---|
| 117 | ConstantResolution = thePar->ConstantResolution(); |
---|
| 118 | NoiseResolution = thePar->NoiseResolution(); |
---|
| 119 | |
---|
| 120 | // Inits |
---|
| 121 | NSpot = 0.00; |
---|
| 122 | AlphaNSpot = 0.00; |
---|
| 123 | TNSpot = 0.00; |
---|
| 124 | BetaNSpot = 0.00; |
---|
| 125 | RadiusCore = 0.00; |
---|
| 126 | WeightCore = 0.00; |
---|
| 127 | RadiusTail = 0.00; |
---|
| 128 | ComputeZAX0EFFetc(); |
---|
| 129 | |
---|
| 130 | G4cout << "/********************************************/ " << G4endl; |
---|
| 131 | G4cout << " - GFlashSamplingShowerParameterisation::Constructor - " << G4endl; |
---|
| 132 | G4cout << "/********************************************/ " << G4endl; |
---|
| 133 | } |
---|
| 134 | |
---|
| 135 | // ------------------------------------------------------------ |
---|
| 136 | |
---|
| 137 | GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation() |
---|
| 138 | {} |
---|
| 139 | |
---|
| 140 | // ------------------------------------------------------------ |
---|
| 141 | |
---|
| 142 | void GFlashSamplingShowerParameterisation:: |
---|
| 143 | SetMaterial(G4Material *mat1, G4Material *mat2) |
---|
| 144 | { |
---|
| 145 | G4double Es = 21*MeV; |
---|
| 146 | material1= mat1; |
---|
| 147 | Z1 = GetEffZ(material1); |
---|
| 148 | A1 = GetEffA(material1); |
---|
| 149 | density1 = material1->GetDensity(); |
---|
| 150 | X01 = material1->GetRadlen(); |
---|
| 151 | Ec1 = 2.66 * std::pow((X01 * Z1 / A1),1.1); |
---|
| 152 | Rm1 = X01*Es/Ec1; |
---|
| 153 | |
---|
| 154 | material2= mat2; |
---|
| 155 | Z2 = GetEffZ(material2); |
---|
| 156 | A2 = GetEffA(material2); |
---|
| 157 | density2 = material2->GetDensity(); |
---|
| 158 | X02 = material2->GetRadlen(); |
---|
| 159 | Ec2 = 2.66 * std::pow((X02 * Z2 / A2),1.1); |
---|
| 160 | Rm2 = X02*Es/Ec2; |
---|
| 161 | // PrintMaterial(); |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | // ------------------------------------------------------------ |
---|
| 165 | |
---|
| 166 | void GFlashSamplingShowerParameterisation::ComputeZAX0EFFetc() |
---|
| 167 | { |
---|
| 168 | G4cout << "/************ ComputeZAX0EFFetc ************/" << G4endl; |
---|
| 169 | G4cout << " - GFlashSamplingShowerParameterisation::Material - " << G4endl; |
---|
| 170 | |
---|
| 171 | G4double Es = 21*MeV; //constant |
---|
| 172 | |
---|
| 173 | // material and geometry parameters for a sampling calorimeter |
---|
| 174 | G4double denominator = (d1*density1 + d2*density2); |
---|
| 175 | G4double W1 = (d1*density1) / denominator; |
---|
| 176 | G4double W2 = (d2*density2)/denominator; |
---|
| 177 | Zeff = ( W1*Z2 ) + (W2*Z1); //X0*Es/Ec; |
---|
| 178 | Aeff = ( W1*A1 ) + (W2*A2); |
---|
| 179 | X0eff =(1/ (( W1 / X01) +( W2 / X02))); |
---|
| 180 | Rhoeff = ( (d1 *density1 ) + (d2 * density2 ))/G4double (d2 + d1 ); |
---|
| 181 | Rmeff = 1/ ((((W1*Ec1)/ X01) + ((W2* Ec2)/ X02) ) / Es ) ; |
---|
| 182 | Eceff = X0eff *((W1*Ec1)/ X01 + (W2* Ec2)/ X02 ); |
---|
| 183 | Fs = X0eff/G4double ((d1/mm )+(d2/mm) ); |
---|
| 184 | ehat = (1. / (1+ 0.007*(Z1- Z2))); |
---|
| 185 | |
---|
| 186 | G4cout << "W1= " << W1 << G4endl; |
---|
| 187 | G4cout << "W2= " << W2 << G4endl; |
---|
| 188 | G4cout << "effective quantities Zeff = "<<Zeff<< G4endl; |
---|
| 189 | G4cout << "effective quantities Aeff = "<<Aeff<< G4endl; |
---|
| 190 | G4cout << "effective quantities Rhoeff = "<<Rhoeff/g *cm3<<" g/cm3" << G4endl; |
---|
| 191 | G4cout << "effective quantities X0eff = "<<X0eff/cm <<" cm" << G4endl; |
---|
| 192 | |
---|
| 193 | X0eff = X0eff * Rhoeff; |
---|
| 194 | |
---|
| 195 | G4cout << "effective quantities X0eff = "<<X0eff/g*cm2 <<" g/cm2" << G4endl; |
---|
| 196 | X0eff = X0eff /Rhoeff; |
---|
| 197 | G4cout << "effective quantities RMeff = "<<Rmeff/cm<<" cm" << G4endl; |
---|
| 198 | Rmeff = Rmeff* Rhoeff; |
---|
| 199 | G4cout << "effective quantities RMeff = "<<Rmeff/g *cm2<<" g/cm2" << G4endl; |
---|
| 200 | Rmeff = Rmeff/ Rhoeff; |
---|
| 201 | G4cout << "effective quantities Eceff = "<<Eceff/MeV<< " MeV"<< G4endl; |
---|
| 202 | G4cout << "effective quantities Fs = "<<Fs<<G4endl; |
---|
| 203 | G4cout << "effective quantities ehat = "<<ehat<<G4endl; |
---|
| 204 | G4cout << "/********************************************/ " <<G4endl; |
---|
| 205 | } |
---|
| 206 | |
---|
| 207 | // ------------------------------------------------------------ |
---|
| 208 | |
---|
| 209 | void GFlashSamplingShowerParameterisation:: |
---|
| 210 | GenerateLongitudinalProfile(G4double Energy) |
---|
| 211 | { |
---|
| 212 | if ((material1==0) || (material2 ==0)) |
---|
| 213 | { |
---|
| 214 | G4Exception("GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile()", |
---|
| 215 | "InvalidSetup", FatalException, "No material initialized!"); |
---|
| 216 | } |
---|
| 217 | G4double y = Energy/Eceff; |
---|
| 218 | ComputeLongitudinalParameters(y); |
---|
| 219 | GenerateEnergyProfile(y); |
---|
| 220 | GenerateNSpotProfile(y); |
---|
| 221 | } |
---|
| 222 | |
---|
| 223 | // ------------------------------------------------------------ |
---|
| 224 | |
---|
| 225 | void |
---|
| 226 | GFlashSamplingShowerParameterisation::ComputeLongitudinalParameters(G4double y) |
---|
| 227 | { |
---|
| 228 | AveLogTmaxh = std::log(std::max(ParAveT1 +std::log(y),0.1)); //ok |
---|
| 229 | AveLogAlphah = std::log(std::max(ParAveA1 + (ParAveA2+ParAveA3/Zeff)*std::log(y),.1)); //ok |
---|
| 230 | //hom |
---|
| 231 | SigmaLogTmaxh = std::min(0.5,1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ); //ok |
---|
| 232 | SigmaLogAlphah = std::min(0.5,1.00/( ParSigLogA1 + ParSigLogA2*std::log(y))); //ok |
---|
| 233 | Rhoh = ParRho1+ParRho2*std::log(y);//ok |
---|
| 234 | // if sampling |
---|
| 235 | AveLogTmax = std::max(0.1,std::log(std::exp(AveLogTmaxh) |
---|
| 236 | + ParsAveT1/Fs + ParsAveT2*(1-ehat))); //ok |
---|
| 237 | AveLogAlpha = std::max(0.1,std::log(std::exp(AveLogAlphah) |
---|
| 238 | + (ParsAveA1/Fs))); //ok |
---|
| 239 | // |
---|
| 240 | SigmaLogTmax = std::min(0.5,1.00/( ParsSigLogT1 |
---|
| 241 | + ParsSigLogT2*std::log(y)) ); //ok |
---|
| 242 | SigmaLogAlpha = std::min(0.5,1.00/( ParsSigLogA1 |
---|
| 243 | + ParsSigLogA2*std::log(y))); //ok |
---|
| 244 | Rho = ParsRho1+ParsRho2*std::log(y); //ok |
---|
| 245 | } |
---|
| 246 | |
---|
| 247 | // ------------------------------------------------------------ |
---|
| 248 | |
---|
| 249 | void GFlashSamplingShowerParameterisation::GenerateEnergyProfile(G4double /* y */) |
---|
| 250 | { |
---|
| 251 | G4double Correlation1 = std::sqrt((1+Rho)/2); |
---|
| 252 | G4double Correlation2 = std::sqrt((1-Rho)/2); |
---|
| 253 | G4double Correlation1h = std::sqrt((1+Rhoh)/2); |
---|
| 254 | G4double Correlation2h = std::sqrt((1-Rhoh)/2); |
---|
| 255 | G4double Random1 = G4RandGauss::shoot(); |
---|
| 256 | G4double Random2 = G4RandGauss::shoot(); |
---|
| 257 | |
---|
| 258 | Tmax = std::max(1.,std::exp( AveLogTmax + SigmaLogTmax * |
---|
| 259 | (Correlation1*Random1 + Correlation2*Random2) )); |
---|
| 260 | Alpha = std::max(1.1,std::exp( AveLogAlpha + SigmaLogAlpha * |
---|
| 261 | (Correlation1*Random1 - Correlation2*Random2) )); |
---|
| 262 | Beta = (Alpha-1.00)/Tmax; |
---|
| 263 | //Parameters for Enenrgy Profile including correaltion and sigmas |
---|
| 264 | Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh * |
---|
| 265 | (Correlation1h*Random1 + Correlation2h*Random2) ); |
---|
| 266 | Alphah = std::exp( AveLogAlphah + SigmaLogAlphah * |
---|
| 267 | (Correlation1h*Random1 - Correlation2h*Random2) ); |
---|
| 268 | Betah = (Alphah-1.00)/Tmaxh; |
---|
| 269 | } |
---|
| 270 | |
---|
| 271 | // ------------------------------------------------------------ |
---|
| 272 | |
---|
| 273 | void GFlashSamplingShowerParameterisation::GenerateNSpotProfile(const G4double y) |
---|
| 274 | { |
---|
| 275 | TNSpot = Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff); //ok. |
---|
| 276 | TNSpot = std::max(0.5,Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff)); |
---|
| 277 | AlphaNSpot = Alphah * (ParsSpotA1+ParsSpotA2*Zeff); |
---|
| 278 | BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok |
---|
| 279 | NSpot = ParsSpotN1 /SamplingResolution * std::pow(y*Eceff/GeV,ParsSpotN2 ); |
---|
| 280 | } |
---|
| 281 | |
---|
| 282 | // ------------------------------------------------------------ |
---|
| 283 | |
---|
| 284 | G4double |
---|
| 285 | GFlashSamplingShowerParameterisation:: |
---|
| 286 | ApplySampling(const G4double DEne, const G4double ) |
---|
| 287 | { |
---|
| 288 | G4double DEneFluctuated = DEne; |
---|
| 289 | G4double Resolution = std::pow(SamplingResolution,2); |
---|
| 290 | |
---|
| 291 | // +pow(NoiseResolution,2)/ //@@@@@@@@ FIXME |
---|
| 292 | // Energy*(1.*MeV)+ |
---|
| 293 | // pow(ConstantResolution,2)* |
---|
| 294 | // Energy/(1.*MeV); |
---|
| 295 | |
---|
| 296 | if(Resolution >0.0 && DEne > 0.00) |
---|
| 297 | { |
---|
| 298 | G4float x1=DEne/Resolution; |
---|
| 299 | G4float x2 = CLHEP::RandGamma::shoot(x1, 1.0)*Resolution; |
---|
| 300 | DEneFluctuated=x2; |
---|
| 301 | } |
---|
| 302 | return DEneFluctuated; |
---|
| 303 | } |
---|
| 304 | |
---|
| 305 | // ------------------------------------------------------------ |
---|
| 306 | |
---|
| 307 | G4double GFlashSamplingShowerParameterisation:: |
---|
| 308 | IntegrateEneLongitudinal(G4double LongitudinalStep) |
---|
| 309 | { |
---|
| 310 | G4double LongitudinalStepInX0 = LongitudinalStep / X0eff; |
---|
| 311 | G4float x1= Betah*LongitudinalStepInX0; |
---|
| 312 | G4float x2= Alphah; |
---|
| 313 | float x3 = gam(x1,x2); |
---|
| 314 | G4double DEne=x3; |
---|
| 315 | return DEne; |
---|
| 316 | } |
---|
| 317 | |
---|
| 318 | // ------------------------------------------------------------ |
---|
| 319 | |
---|
| 320 | G4double GFlashSamplingShowerParameterisation:: |
---|
| 321 | IntegrateNspLongitudinal(G4double LongitudinalStep) |
---|
| 322 | { |
---|
| 323 | G4double LongitudinalStepInX0 = LongitudinalStep / X0eff; |
---|
| 324 | G4float x1 = BetaNSpot*LongitudinalStepInX0; |
---|
| 325 | G4float x2 = AlphaNSpot; |
---|
| 326 | G4float x3 = gam(x1,x2); |
---|
| 327 | G4double DNsp = x3; |
---|
| 328 | return DNsp; |
---|
| 329 | } |
---|
| 330 | |
---|
| 331 | // ------------------------------------------------------------ |
---|
| 332 | |
---|
| 333 | G4double GFlashSamplingShowerParameterisation:: |
---|
| 334 | GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition) |
---|
| 335 | { |
---|
| 336 | if(ispot < 1) |
---|
| 337 | { |
---|
| 338 | // Determine lateral parameters in the middle of the step. |
---|
| 339 | // They depend on energy & position along step |
---|
| 340 | // |
---|
| 341 | G4double Tau = ComputeTau(LongitudinalPosition); |
---|
| 342 | ComputeRadialParameters(Energy,Tau); |
---|
| 343 | } |
---|
| 344 | |
---|
| 345 | G4double Radius; |
---|
| 346 | G4double Random1 = G4UniformRand(); |
---|
| 347 | G4double Random2 = G4UniformRand(); |
---|
| 348 | if(Random1 <WeightCore) //WeightCore = p < w_i |
---|
| 349 | { |
---|
| 350 | Radius = Rmeff * RadiusCore * std::sqrt( Random2/(1. - Random2) ); |
---|
| 351 | } |
---|
| 352 | else |
---|
| 353 | { |
---|
| 354 | Radius = Rmeff * RadiusTail * std::sqrt( Random2/(1. - Random2) ); |
---|
| 355 | } |
---|
| 356 | Radius = std::min(Radius,DBL_MAX); |
---|
| 357 | return Radius; |
---|
| 358 | } |
---|
| 359 | |
---|
| 360 | // ------------------------------------------------------------ |
---|
| 361 | |
---|
| 362 | G4double |
---|
| 363 | GFlashSamplingShowerParameterisation:: |
---|
| 364 | ComputeTau(G4double LongitudinalPosition) |
---|
| 365 | { |
---|
| 366 | G4double tau = LongitudinalPosition / Tmax/ X0eff //<t> = T* a /(a - 1) |
---|
| 367 | * (Alpha-1.00) /Alpha |
---|
| 368 | * std::exp(AveLogAlpha)/(std::exp(AveLogAlpha)-1.); //ok |
---|
| 369 | return tau; |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | // ------------------------------------------------------------ |
---|
| 373 | |
---|
| 374 | void GFlashSamplingShowerParameterisation:: |
---|
| 375 | ComputeRadialParameters(G4double Energy, G4double Tau) |
---|
| 376 | { |
---|
| 377 | G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV); //ok |
---|
| 378 | G4double z2 = ParRC3+ParRC4*Zeff; //ok |
---|
| 379 | RadiusCore = z1 + z2 * Tau; //ok |
---|
| 380 | G4double p1 = ParWC1+ParWC2*Zeff; //ok |
---|
| 381 | G4double p2 = ParWC3+ParWC4*Zeff; //ok |
---|
| 382 | G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok |
---|
| 383 | WeightCore = p1 * std::exp( (p2-Tau)/p3- std::exp( (p2-Tau) /p3) ); //ok |
---|
| 384 | |
---|
| 385 | G4double k1 = ParRT1+ParRT2*Zeff; // ok |
---|
| 386 | G4double k2 = ParRT3; // ok |
---|
| 387 | G4double k3 = ParRT4; // ok |
---|
| 388 | G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok |
---|
| 389 | |
---|
| 390 | RadiusTail = k1*(std::exp(k3*(Tau-k2)) |
---|
| 391 | + std::exp(k4*(Tau-k2)) ); //ok |
---|
| 392 | |
---|
| 393 | // sampling calorimeter |
---|
| 394 | |
---|
| 395 | RadiusCore = RadiusCore + ParsRC1*(1-ehat) + ParsRC2/Fs*std::exp(-Tau); //ok |
---|
| 396 | WeightCore = WeightCore + (1-ehat) |
---|
| 397 | * (ParsWC1+ParsWC2/Fs * std::exp(-std::pow((Tau-1.),2))); //ok |
---|
| 398 | RadiusTail = RadiusTail + (1-ehat)* ParsRT1+ ParsRT2/Fs *std::exp(-Tau); //ok |
---|
| 399 | } |
---|
| 400 | |
---|
| 401 | // ------------------------------------------------------------ |
---|
| 402 | |
---|
| 403 | G4double GFlashSamplingShowerParameterisation:: |
---|
| 404 | GenerateExponential(const G4double /* Energy */ ) |
---|
| 405 | { |
---|
| 406 | G4double ParExp1 = 9./7.*X0eff; |
---|
| 407 | G4double random = -ParExp1*CLHEP::RandExponential::shoot() ; |
---|
| 408 | return random; |
---|
| 409 | } |
---|