Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // 28 // ------------------------------------------- 29 // GEANT 4 class implementation 30 // 31 // ------- GFlashSamplingShowerParameteri 32 // 33 // Authors: E.Barberio & Joanna Weng - 11.2005 34 // ------------------------------------------- 35 36 #include <cmath> 37 38 #include "GFlashSamplingShowerParameterisation 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::GFlashSa 47 G4Material* aMat1, G4Material* aMat2, G4doub 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 homogenio 87 88 // shower max 89 ParAveT1 = thePar->ParAveT1(); // ln (ln y 90 ParAveA1 = thePar->ParAveA1(); // ln a (0.8 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 97 ParSigLogT2 = thePar->ParsSigLogT2(); // le 98 // variance of 'alpha' 99 ParSigLogA1 = 100 thePar 101 ->ParSigLogA1(); // Sigma a (-0.58 + 0. 102 ParSigLogA2 = thePar->ParSigLogA2(); // lea 103 // correlation alpha%T 104 ParRho1 = thePar->ParRho1(); // Rho = 0.705 105 ParRho2 = thePar->ParRho2(); // leaving Par 106 // Sampling 107 ParsAveT1 = thePar->ParsAveT1(); // T_sam = 108 ParsAveT2 = thePar->ParsAveT2(); 109 ParsAveA1 = thePar->ParsAveA1(); 110 // Variance of shower max sampling 111 ParsSigLogT1 = thePar->ParsSigLogT1(); // S 112 // w 113 ParsSigLogT2 = thePar->ParsSigLogT2(); 114 // variance of 'alpha' 115 ParsSigLogA1 = thePar->ParsSigLogA1(); // S 116 // w 117 ParsSigLogA2 = thePar->ParsSigLogA2(); 118 // correlation alpha%T 119 ParsRho1 = 120 thePar->ParsRho1(); // Rho = 0.784 -0.023 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 ))+st 126 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 127 ParRC2 = thePar->ParRC2(); 128 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 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 fo 153 ParsSpotT1 = thePar->ParSpotT1(); // T_spot 154 ParsSpotT2 = thePar->ParSpotT2(); 155 ParsSpotA1 = thePar->ParSpotA1(); // a_spot 156 ParsSpotA2 = thePar->ParSpotA2(); 157 ParsSpotN1 = thePar->ParSpotN1(); // N_Spot 158 ParsSpotN2 = thePar->ParSpotN2(); 159 SamplingResolution = thePar->SamplingResolut 160 ConstantResolution = thePar->ConstantResolut 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 << "/******************************** 174 G4cout << " - GFlashSamplingShowerParameter 175 G4cout << "/******************************** 176 } 177 178 // ------------------------------------------- 179 180 GFlashSamplingShowerParameterisation::~GFlashS 181 { 182 delete thePar; 183 } 184 185 // ------------------------------------------- 186 187 void GFlashSamplingShowerParameterisation::Set 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::Com 213 { 214 G4cout << "/************ ComputeZAX0EFFetc * 215 G4cout << " - GFlashSamplingShowerParameter 216 217 const G4double Es = 21.2 * MeV; 218 219 // material and geometry parameters for a sa 220 G4double denominator = (d1 * density1 + d2 * 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)) 226 X0eff = (W1 * Rhoeff) / (X01 * density1) + ( 227 X0eff = 1. / X0eff; 228 Rmeff = 1. / ((((W1 * Ec1) / X01) + ((W2 * E 229 Eceff = X0eff * ((W1 * Ec1) / X01 + (W2 * Ec 230 Fs = X0eff / (d1 + d2); // --> was G4double 231 // mm makes sense.. 232 ehat = (1. / (1 + 0.007 * (Z1 - Z2))); 233 234 G4cout << "W1= " << W1 << G4endl; 235 G4cout << "W2= " << W2 << G4endl; 236 G4cout << "effective quantities Zeff = " << 237 G4cout << "effective quantities Aeff = " << 238 G4cout << "effective quantities Rhoeff = " < 239 G4cout << "effective quantities X0eff = " << 240 241 X0eff = X0eff * Rhoeff; 242 243 G4cout << "effective quantities X0eff = " << 244 X0eff = X0eff / Rhoeff; 245 G4cout << "effective quantities RMeff = " << 246 Rmeff = Rmeff * Rhoeff; 247 G4cout << "effective quantities RMeff = " << 248 Rmeff = Rmeff / Rhoeff; 249 G4cout << "effective quantities Eceff = " << 250 G4cout << "effective quantities Fs = " << Fs 251 G4cout << "effective quantities ehat = " << 252 G4cout << "/******************************** 253 } 254 255 // ------------------------------------------- 256 257 void GFlashSamplingShowerParameterisation::Gen 258 { 259 if ((material1 == 0) || (material2 == 0)) { 260 G4Exception("GFlashSamplingShowerParameter 261 "InvalidSetup", FatalException 262 } 263 G4double y = Energy / Eceff; 264 ComputeLongitudinalParameters(y); 265 GenerateEnergyProfile(y); 266 GenerateNSpotProfile(y); 267 } 268 269 // ------------------------------------------- 270 271 void GFlashSamplingShowerParameterisation::Com 272 { 273 AveLogTmaxh = std::log(std::max(ParAveT1 + s 274 AveLogAlphah = 275 std::log(std::max(ParAveA1 + (ParAveA2 + P 276 // hom 277 SigmaLogTmaxh = std::min(0.5, 1.00 / (ParSig 278 SigmaLogAlphah = std::min(0.5, 1.00 / (ParSi 279 Rhoh = ParRho1 + ParRho2 * std::log(y); // 280 // if sampling 281 AveLogTmax = 282 std::max(0.1, std::log(std::exp(AveLogTmax 283 AveLogAlpha = std::max(0.1, std::log(std::ex 284 // 285 SigmaLogTmax = std::min(0.5, 1.00 / (ParsSig 286 SigmaLogAlpha = std::min(0.5, 1.00 / (ParsSi 287 Rho = ParsRho1 + ParsRho2 * std::log(y); // 288 289 if (0) { 290 G4cout << " y = " << y << G4end 291 G4cout << " std::log(std::exp(AveLogTmaxh) 292 << " std::log(" << std::exp(AveLogT 293 << ParsAveT2 * (1 - ehat) << ") = " 294 << " std::log(" << std::exp(AveLogT 295 << ParsAveT2 << "*" << (1 - ehat) < 296 << " std::log(" << std::exp(AveLogT 297 << G4endl; 298 G4cout << " AveLogTmaxh " << AveLogTmax 299 G4cout << " AveLogAlphah " << AveLogAlph 300 G4cout << " SigmaLogTmaxh " << SigmaLogTm 301 G4cout << " 1.00/( ParSigLogT1 + ParSigLog 302 << (ParSigLogT1 + ParSigLogT2 * std 303 << ParSigLogT1 << " + " << ParSigLo 304 << ParSigLogT1 << " + " << ParSigLo 305 G4cout << " SigmaLogAlphah " << SigmaLogAl 306 G4cout << " Rhoh " << Rhoh << G4 307 G4cout << " AveLogTmax " << AveLogTmax 308 G4cout << " AveLogAlpha " << AveLogAlph 309 G4cout << " SigmaLogTmax " << SigmaLogTm 310 G4cout << " SigmaLogAlpha " << SigmaLogAl 311 G4cout << " Rho " << Rho << G4e 312 } 313 } 314 315 // ------------------------------------------- 316 317 void GFlashSamplingShowerParameterisation::Gen 318 { 319 G4double Correlation1 = std::sqrt((1 + Rho) 320 G4double Correlation2 = std::sqrt((1 - Rho) 321 G4double Correlation1h = std::sqrt((1 + Rhoh 322 G4double Correlation2h = std::sqrt((1 - Rhoh 323 G4double Random1 = G4RandGauss::shoot(); 324 G4double Random2 = G4RandGauss::shoot(); 325 326 Tmax = std::max( 327 1., std::exp(AveLogTmax + SigmaLogTmax * ( 328 Alpha = std::max( 329 1.1, std::exp(AveLogAlpha + SigmaLogAlpha 330 Beta = (Alpha - 1.00) / Tmax; 331 // Parameters for Enenrgy Profile including 332 Tmaxh = 333 std::exp(AveLogTmaxh + SigmaLogTmaxh * (Co 334 Alphah = 335 std::exp(AveLogAlphah + SigmaLogAlphah * ( 336 Betah = (Alphah - 1.00) / Tmaxh; 337 } 338 339 // ------------------------------------------- 340 341 void GFlashSamplingShowerParameterisation::Gen 342 { 343 TNSpot = Tmaxh * (ParsSpotT1 + ParsSpotT2 * 344 TNSpot = std::max(0.5, Tmaxh * (ParsSpotT1 + 345 AlphaNSpot = Alphah * (ParsSpotA1 + ParsSpot 346 BetaNSpot = (AlphaNSpot - 1.00) / TNSpot; / 347 NSpot = ParsSpotN1 / SamplingResolution * st 348 } 349 350 // ------------------------------------------- 351 352 G4double GFlashSamplingShowerParameterisation: 353 { 354 G4double DEneFluctuated = DEne; 355 G4double Resolution = std::pow(SamplingResol 356 357 // +pow(NoiseResolution,2)/ //@@@@@@@ 358 // Energy*(1.*MeV)+ 359 // pow(ConstantResol 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 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 Longitudinal 377 { 378 G4double LongitudinalStepInX0 = Longitudinal 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: 389 { 390 G4double LongitudinalStepInX0 = Longitudinal 391 G4float x1 = BetaNSpot * LongitudinalStepInX 392 G4float x2 = AlphaNSpot; 393 G4float x3 = gam(x1, x2); 394 G4double DNsp = x3; 395 return DNsp; 396 } 397 398 // ------------------------------------------- 399 400 G4double GFlashSamplingShowerParameterisation: 401 402 { 403 if (ispot < 1) { 404 // Determine lateral parameters in the mid 405 // They depend on energy & position along 406 // 407 G4double Tau = ComputeTau(LongitudinalPosi 408 ComputeRadialParameters(Energy, Tau); 409 } 410 411 G4double Radius; 412 G4double Random1 = G4UniformRand(); 413 G4double Random2 = G4UniformRand(); 414 if (Random1 < WeightCore) // WeightCore = p 415 { 416 Radius = Rmeff * RadiusCore * std::sqrt(Ra 417 } 418 else { 419 Radius = Rmeff * RadiusTail * std::sqrt(Ra 420 } 421 Radius = std::min(Radius, DBL_MAX); 422 return Radius; 423 } 424 425 // ------------------------------------------- 426 427 G4double GFlashSamplingShowerParameterisation: 428 { 429 G4double tau = LongitudinalPosition / Tmax / 430 * (Alpha - 1.00) / Alpha * st 431 / (std::exp(AveLogAlpha) - 1. 432 return tau; 433 } 434 435 // ------------------------------------------- 436 437 void GFlashSamplingShowerParameterisation::Com 438 { 439 G4double z1 = ParRC1 + ParRC2 * std::log(Ene 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(Ene 445 WeightCore = p1 * std::exp((p2 - Tau) / p3 - 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(Ene 451 452 RadiusTail = k1 * (std::exp(k3 * (Tau - k2)) 453 454 // sampling calorimeter 455 456 RadiusCore = RadiusCore + ParsRC1 * (1 - eha 457 WeightCore = 458 WeightCore + (1 - ehat) * (ParsWC1 + ParsW 459 RadiusTail = RadiusTail + (1 - ehat) * ParsR 460 } 461 462 // ------------------------------------------- 463 464 G4double GFlashSamplingShowerParameterisation: 465 GenerateExponential(const G4double /* Energy * 466 { 467 G4double ParExp1 = 9./7.*X0eff; 468 G4double random = -ParExp1*G4RandExponentia 469 return random; 470 } 471