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 // ------- GFlashHomoShowerParameterisati 32 // 33 // Authors: E.Barberio & Joanna Weng - 9.11.20 34 // ------------------------------------------- 35 36 #include <cmath> 37 38 #include "GFlashHomoShowerParameterisation.hh" 39 #include "GVFlashShowerParameterisation.hh" 40 #include "G4PhysicalConstants.hh" 41 #include "G4SystemOfUnits.hh" 42 #include "Randomize.hh" 43 #include "G4ios.hh" 44 #include "G4Material.hh" 45 #include "G4MaterialTable.hh" 46 47 GFlashHomoShowerParameterisation::GFlashHomoSh 48 49 : GVFlashShowerParameterisation(), 50 ConstantResolution(0.), 51 NoiseResolution(0.), 52 SamplingResolution(0.), 53 AveLogAlphah(0.), 54 AveLogTmaxh(0.), 55 SigmaLogAlphah(0.), 56 SigmaLogTmaxh(0.), 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 70 SetMaterial(aMat); 71 PrintMaterial(aMat); 72 73 /******************************************* 74 /* Homo Calorimeter 75 /******************************************* 76 // Longitudinal Coefficients for a homogenio 77 // shower max 78 // 79 ParAveT1 = thePar->ParAveT1(); // ln (ln y 80 ParAveA1 = thePar->ParAveA1(); // ln a (0.8 81 ParAveA2 = thePar->ParAveA2(); 82 ParAveA3 = thePar->ParAveA3(); 83 84 // Variance of shower max 85 ParSigLogT1 = thePar->ParSigLogT1(); // Sig 86 ParSigLogT2 = thePar->ParSigLogT2(); 87 88 // variance of 'alpha' 89 // 90 ParSigLogA1 = thePar->ParSigLogA1(); // Sig 91 ParSigLogA2 = thePar->ParSigLogA2(); 92 93 // correlation alpha%T 94 // 95 ParRho1 = thePar->ParRho1(); // Rho = 0.705 96 ParRho2 = thePar->ParRho2(); 97 98 // Radial Coefficients 99 // r_C (tau)= z_1 +z_2 tau 100 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+st 101 // 102 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 103 ParRC2 = thePar->ParRC2(); 104 105 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 106 ParRC4 = thePar->ParRC4(); 107 108 ParWC1 = thePar->ParWC1(); 109 ParWC2 = thePar->ParWC2(); 110 ParWC3 = thePar->ParWC3(); 111 ParWC4 = thePar->ParWC4(); 112 ParWC5 = thePar->ParWC5(); 113 ParWC6 = thePar->ParWC6(); 114 115 ParRT1 = thePar->ParRT1(); 116 ParRT2 = thePar->ParRT2(); 117 ParRT3 = thePar->ParRT3(); 118 ParRT4 = thePar->ParRT4(); 119 ParRT5 = thePar->ParRT5(); 120 ParRT6 = thePar->ParRT6(); 121 122 // Coeff for fluctueted radial profiles for 123 // 124 ParSpotT1 = thePar->ParSpotT1(); // T_spot 125 ParSpotT2 = thePar->ParSpotT2(); 126 127 ParSpotA1 = thePar->ParSpotA1(); // a_spot= 128 ParSpotA2 = thePar->ParSpotA2(); 129 130 ParSpotN1 = thePar->ParSpotN1(); // N_Spot 131 ParSpotN2 = thePar->ParSpotN2(); 132 133 // Inits 134 135 NSpot = 0.00; 136 AlphaNSpot = 0.00; 137 TNSpot = 0.00; 138 BetaNSpot = 0.00; 139 140 RadiusCore = 0.00; 141 WeightCore = 0.00; 142 RadiusTail = 0.00; 143 144 G4cout << "/******************************** 145 G4cout << " - GFlashHomoShowerParameterisat 146 G4cout << "/******************************** 147 } 148 149 void GFlashHomoShowerParameterisation::SetMate 150 { 151 material = mat; 152 Z = GetEffZ(material); 153 A = GetEffA(material); 154 density = material->GetDensity() / (g / cm3) 155 X0 = material->GetRadlen(); 156 // O. I. Dovzhenkko and A. A. Pommanskii 157 Ec = 2.66 * std::pow((X0 * Z / A), 1.1); 158 // // Rossi appriximation 159 // Ec = 610.0 * MeV / (Z + 1.24); 160 const G4double Es = 21.2 * MeV; 161 Rm = X0 * Es / Ec; 162 // PrintMaterial(); 163 } 164 165 GFlashHomoShowerParameterisation::~GFlashHomoS 166 { 167 delete thePar; 168 } 169 170 void GFlashHomoShowerParameterisation::Generat 171 { 172 if (material == 0) { 173 G4Exception("GFlashHomoShowerParameterisat 174 FatalException, "No material i 175 } 176 177 G4double y = Energy / Ec; 178 ComputeLongitudinalParameters(y); 179 GenerateEnergyProfile(y); 180 GenerateNSpotProfile(y); 181 } 182 183 void GFlashHomoShowerParameterisation::Compute 184 { 185 AveLogTmaxh = std::log(ParAveT1 + std::log(y 186 // ok <ln T hom> 187 AveLogAlphah = std::log(ParAveA1 + (ParAveA2 188 // ok <ln alpha hom> 189 190 SigmaLogTmaxh = 1.00 / (ParSigLogT1 + ParSig 191 // ok sigma (ln T hom) 192 SigmaLogAlphah = 1.00 / (ParSigLogA1 + ParSi 193 // ok sigma (ln alpha hom) 194 Rhoh = ParRho1 + ParRho2 * std::log(y); // 195 } 196 197 void GFlashHomoShowerParameterisation::Generat 198 { 199 G4double Correlation1h = std::sqrt((1 + Rhoh 200 G4double Correlation2h = std::sqrt((1 - Rhoh 201 202 G4double Random1 = G4RandGauss::shoot(); 203 G4double Random2 = G4RandGauss::shoot(); 204 205 // Parameters for Enenrgy Profile including 206 207 Tmaxh = 208 std::exp(AveLogTmaxh + SigmaLogTmaxh * (Co 209 Alphah = 210 std::exp(AveLogAlphah + SigmaLogAlphah * ( 211 Betah = (Alphah - 1.00) / Tmaxh; 212 } 213 214 void GFlashHomoShowerParameterisation::Generat 215 { 216 TNSpot = Tmaxh * (ParSpotT1 + ParSpotT2 * Z) 217 AlphaNSpot = Alphah * (ParSpotA1 + ParSpotA2 218 BetaNSpot = (AlphaNSpot - 1.00) / TNSpot; / 219 NSpot = ParSpotN1 * std::log(Z) * std::pow(( 220 } 221 222 G4double GFlashHomoShowerParameterisation:: 223 IntegrateEneLongitudinal(G4double Longitudinal 224 { 225 G4double LongitudinalStepInX0 = Longitudinal 226 G4float x1= Betah*LongitudinalStepInX0; 227 G4float x2= Alphah; 228 float x3 = gam(x1,x2); 229 G4double DEne=x3; 230 return DEne; 231 } 232 233 G4double GFlashHomoShowerParameterisation::Int 234 { 235 G4double LongitudinalStepInX0 = Longitudinal 236 G4float x1 = BetaNSpot * LongitudinalStepInX 237 G4float x2 = AlphaNSpot; 238 G4float x3 = gam(x1, x2); 239 G4double DNsp = x3; 240 return DNsp; 241 } 242 243 G4double GFlashHomoShowerParameterisation::Gen 244 245 { 246 if (ispot < 1) { 247 // Determine lateral parameters in the mid 248 // They depend on energy & position along 249 // 250 G4double Tau = ComputeTau(LongitudinalPosi 251 ComputeRadialParameters(Energy, Tau); 252 } 253 254 G4double Radius; 255 G4double Random1 = G4UniformRand(); 256 G4double Random2 = G4UniformRand(); 257 258 if (Random1 < WeightCore) // WeightCore = p 259 { 260 Radius = Rm * RadiusCore * std::sqrt(Rando 261 } 262 else { 263 Radius = Rm * RadiusTail * std::sqrt(Rando 264 } 265 Radius = std::min(Radius, DBL_MAX); 266 return Radius; 267 } 268 269 G4double GFlashHomoShowerParameterisation::Com 270 { 271 G4double tau = LongitudinalPosition / Tmaxh 272 * (Alphah - 1.00) / Alphah * 273 / (std::exp(AveLogAlphah) - 1 274 return tau; 275 } 276 277 void GFlashHomoShowerParameterisation::Compute 278 { 279 G4double z1 = ParRC1 + ParRC2 * std::log(Ene 280 G4double z2 = ParRC3 + ParRC4 * Z; // ok 281 RadiusCore = z1 + z2 * Tau; // ok 282 283 G4double p1 = ParWC1 + ParWC2 * Z; // ok 284 G4double p2 = ParWC3 + ParWC4 * Z; // ok 285 G4double p3 = ParWC5 + ParWC6 * std::log(Ene 286 287 WeightCore = p1 * std::exp((p2 - Tau) / p3 - 288 289 G4double k1 = ParRT1 + ParRT2 * Z; // ok 290 G4double k2 = ParRT3; // ok 291 G4double k3 = ParRT4; // ok 292 G4double k4 = ParRT5 + ParRT6 * std::log(Ene 293 294 RadiusTail = k1 * (std::exp(k3 * (Tau - k2)) 295 } 296 297 G4double GFlashHomoShowerParameterisation:: 298 GenerateExponential(const G4double /* Energy * 299 { 300 G4double ParExp1 = 9./7.*X0; 301 G4double random = -ParExp1*G4RandExponentia 302 return random; 303 } 304