[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: GFlashHomoShowerParameterisation.cc,v 1.5 2006/06/29 19:14:14 gunter Exp $ |
---|
[1058] | 27 | // GEANT4 tag $Name: geant4-09-02-ref-02 $ |
---|
[823] | 28 | // |
---|
| 29 | // |
---|
| 30 | // ------------------------------------------------------------ |
---|
| 31 | // GEANT 4 class implementation |
---|
| 32 | // |
---|
| 33 | // ------- GFlashHomoShowerParameterisation ------- |
---|
| 34 | // |
---|
| 35 | // Authors: E.Barberio & Joanna Weng - 9.11.2004 |
---|
| 36 | // ------------------------------------------------------------ |
---|
| 37 | |
---|
| 38 | #include "GVFlashShowerParameterisation.hh" |
---|
| 39 | #include "GFlashHomoShowerParameterisation.hh" |
---|
| 40 | #include <cmath> |
---|
| 41 | #include "Randomize.hh" |
---|
| 42 | #include "G4ios.hh" |
---|
| 43 | #include "G4Material.hh" |
---|
| 44 | #include "G4MaterialTable.hh" |
---|
| 45 | |
---|
| 46 | GFlashHomoShowerParameterisation:: |
---|
| 47 | GFlashHomoShowerParameterisation(G4Material * aMat, |
---|
| 48 | GVFlashHomoShowerTuning * aPar) |
---|
| 49 | : GVFlashShowerParameterisation() |
---|
| 50 | { |
---|
| 51 | if(!aPar) { thePar = new GVFlashHomoShowerTuning; } |
---|
| 52 | else { thePar = aPar; } |
---|
| 53 | |
---|
| 54 | SetMaterial(aMat); |
---|
| 55 | PrintMaterial(aMat); |
---|
| 56 | |
---|
| 57 | /********************************************/ |
---|
| 58 | /* Homo Calorimeter */ |
---|
| 59 | /********************************************/ |
---|
| 60 | // Longitudinal Coefficients for a homogenious calo |
---|
| 61 | // shower max |
---|
| 62 | // |
---|
| 63 | ParAveT1 = thePar->ParAveT1(); // ln (ln y -0.812) |
---|
| 64 | ParAveA1 = thePar->ParAveA1(); // ln a (0.81 + (0.458 + 2.26/Z)ln y) |
---|
| 65 | ParAveA2 = thePar->ParAveA2(); |
---|
| 66 | ParAveA3 = thePar->ParAveA3(); |
---|
| 67 | |
---|
| 68 | // Variance of shower max |
---|
| 69 | ParSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-1.4 + 1.26 ln y)**-1 |
---|
| 70 | ParSigLogT2 = thePar->ParSigLogT2(); |
---|
| 71 | |
---|
| 72 | // variance of 'alpha' |
---|
| 73 | // |
---|
| 74 | ParSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.58 + 0.86 ln y)**-1 |
---|
| 75 | ParSigLogA2 = thePar->ParSigLogA2(); |
---|
| 76 | |
---|
| 77 | // correlation alpha%T |
---|
| 78 | // |
---|
| 79 | ParRho1 = thePar->ParRho1(); // Rho = 0.705 -0.023 ln y |
---|
| 80 | ParRho2 = thePar->ParRho2(); |
---|
| 81 | |
---|
| 82 | // Radial Coefficients |
---|
| 83 | // r_C (tau)= z_1 +z_2 tau |
---|
| 84 | // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2)))) |
---|
| 85 | // |
---|
| 86 | ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E |
---|
| 87 | ParRC2 = thePar->ParRC2(); |
---|
| 88 | |
---|
| 89 | ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z |
---|
| 90 | ParRC4 = thePar->ParRC4(); |
---|
| 91 | |
---|
| 92 | ParWC1 = thePar->ParWC1(); |
---|
| 93 | ParWC2 = thePar->ParWC2(); |
---|
| 94 | ParWC3 = thePar->ParWC3(); |
---|
| 95 | ParWC4 = thePar->ParWC4(); |
---|
| 96 | ParWC5 = thePar->ParWC5(); |
---|
| 97 | ParWC6 = thePar->ParWC6(); |
---|
| 98 | |
---|
| 99 | ParRT1 = thePar->ParRT1(); |
---|
| 100 | ParRT2 = thePar->ParRT2(); |
---|
| 101 | ParRT3 = thePar->ParRT3(); |
---|
| 102 | ParRT4 = thePar->ParRT4(); |
---|
| 103 | ParRT5 = thePar->ParRT5(); |
---|
| 104 | ParRT6 = thePar->ParRT6(); |
---|
| 105 | |
---|
| 106 | // Coeff for fluctueted radial profiles for a uniform media |
---|
| 107 | // |
---|
| 108 | ParSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212) |
---|
| 109 | ParSpotT2 = thePar->ParSpotT2(); |
---|
| 110 | |
---|
| 111 | ParSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334) |
---|
| 112 | ParSpotA2 = thePar->ParSpotA2(); |
---|
| 113 | |
---|
| 114 | ParSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876 |
---|
| 115 | ParSpotN2 = thePar->ParSpotN2(); |
---|
| 116 | |
---|
| 117 | // Inits |
---|
| 118 | |
---|
| 119 | NSpot = 0.00; |
---|
| 120 | AlphaNSpot = 0.00; |
---|
| 121 | TNSpot = 0.00; |
---|
| 122 | BetaNSpot = 0.00; |
---|
| 123 | |
---|
| 124 | RadiusCore = 0.00; |
---|
| 125 | WeightCore = 0.00; |
---|
| 126 | RadiusTail = 0.00; |
---|
| 127 | |
---|
| 128 | G4cout << "/********************************************/ " << G4endl; |
---|
| 129 | G4cout << " - GFlashHomoShowerParameterisation::Constructor - " << G4endl; |
---|
| 130 | G4cout << "/********************************************/ " << G4endl; |
---|
| 131 | } |
---|
| 132 | |
---|
| 133 | void GFlashHomoShowerParameterisation::SetMaterial(G4Material *mat) |
---|
| 134 | { |
---|
| 135 | material= mat; |
---|
| 136 | Z = GetEffZ(material); |
---|
| 137 | A = GetEffA(material); |
---|
| 138 | density = material->GetDensity()/(g/cm3); |
---|
| 139 | X0 = material->GetRadlen(); |
---|
| 140 | Ec = 2.66 * std::pow((X0 * Z / A),1.1); |
---|
| 141 | G4double Es = 21*MeV; |
---|
| 142 | Rm = X0*Es/Ec; |
---|
| 143 | // PrintMaterial(); |
---|
| 144 | } |
---|
| 145 | |
---|
| 146 | GFlashHomoShowerParameterisation::~GFlashHomoShowerParameterisation() |
---|
| 147 | {} |
---|
| 148 | |
---|
| 149 | void GFlashHomoShowerParameterisation:: |
---|
| 150 | GenerateLongitudinalProfile(G4double Energy) |
---|
| 151 | { |
---|
| 152 | if (material==0) |
---|
| 153 | { |
---|
| 154 | G4Exception("GFlashHomoShowerParameterisation::GenerateLongitudinalProfile()", |
---|
| 155 | "InvalidSetup", FatalException, "No material initialized!"); |
---|
| 156 | } |
---|
| 157 | |
---|
| 158 | G4double y = Energy/Ec; |
---|
| 159 | ComputeLongitudinalParameters(y); |
---|
| 160 | GenerateEnergyProfile(y); |
---|
| 161 | GenerateNSpotProfile(y); |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | void |
---|
| 165 | GFlashHomoShowerParameterisation::ComputeLongitudinalParameters(G4double y) |
---|
| 166 | { |
---|
| 167 | AveLogTmaxh = std::log(ParAveT1 + std::log(y)); |
---|
| 168 | //ok <ln T hom> |
---|
| 169 | AveLogAlphah = std::log(ParAveA1 + (ParAveA2+ParAveA3/Z)*std::log(y)); |
---|
| 170 | //ok <ln alpha hom> |
---|
| 171 | |
---|
| 172 | SigmaLogTmaxh = 1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ; |
---|
| 173 | //ok sigma (ln T hom) |
---|
| 174 | SigmaLogAlphah = 1.00/( ParSigLogA1 + ParSigLogA2*std::log(y)); |
---|
| 175 | //ok sigma (ln alpha hom) |
---|
| 176 | Rhoh = ParRho1+ParRho2*std::log(y); //ok |
---|
| 177 | } |
---|
| 178 | |
---|
| 179 | void GFlashHomoShowerParameterisation::GenerateEnergyProfile(G4double /* y */) |
---|
| 180 | { |
---|
| 181 | G4double Correlation1h = std::sqrt((1+Rhoh)/2); |
---|
| 182 | G4double Correlation2h = std::sqrt((1-Rhoh)/2); |
---|
| 183 | |
---|
| 184 | G4double Random1 = G4RandGauss::shoot(); |
---|
| 185 | G4double Random2 = G4RandGauss::shoot(); |
---|
| 186 | |
---|
| 187 | // Parameters for Enenrgy Profile including correaltion and sigmas |
---|
| 188 | |
---|
| 189 | Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh * |
---|
| 190 | (Correlation1h*Random1 + Correlation2h*Random2) ); |
---|
| 191 | Alphah = std::exp( AveLogAlphah + SigmaLogAlphah * |
---|
| 192 | (Correlation1h*Random1 - Correlation2h*Random2) ); |
---|
| 193 | Betah = (Alphah-1.00)/Tmaxh; |
---|
| 194 | } |
---|
| 195 | |
---|
| 196 | void GFlashHomoShowerParameterisation::GenerateNSpotProfile(const G4double y) |
---|
| 197 | { |
---|
| 198 | TNSpot = Tmaxh * (ParSpotT1+ParSpotT2*Z); // ok |
---|
| 199 | AlphaNSpot = Alphah * (ParSpotA1+ParSpotA2*Z); |
---|
| 200 | BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok |
---|
| 201 | NSpot = ParSpotN1 * std::log(Z)*std::pow((y*Ec)/GeV,ParSpotN2 ); // ok |
---|
| 202 | } |
---|
| 203 | |
---|
| 204 | G4double GFlashHomoShowerParameterisation:: |
---|
| 205 | IntegrateEneLongitudinal(G4double LongitudinalStep) |
---|
| 206 | { |
---|
| 207 | G4double LongitudinalStepInX0 = LongitudinalStep / X0; |
---|
| 208 | G4float x1= Betah*LongitudinalStepInX0; |
---|
| 209 | G4float x2= Alphah; |
---|
| 210 | float x3 = gam(x1,x2); |
---|
| 211 | G4double DEne=x3; |
---|
| 212 | return DEne; |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | G4double GFlashHomoShowerParameterisation:: |
---|
| 216 | IntegrateNspLongitudinal(G4double LongitudinalStep) |
---|
| 217 | { |
---|
| 218 | G4double LongitudinalStepInX0 = LongitudinalStep / X0; |
---|
| 219 | G4float x1 = BetaNSpot*LongitudinalStepInX0; |
---|
| 220 | G4float x2 = AlphaNSpot; |
---|
| 221 | G4float x3 = gam(x1,x2); |
---|
| 222 | G4double DNsp = x3; |
---|
| 223 | return DNsp; |
---|
| 224 | } |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | G4double GFlashHomoShowerParameterisation:: |
---|
| 228 | GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition) |
---|
| 229 | { |
---|
| 230 | if(ispot < 1) |
---|
| 231 | { |
---|
| 232 | // Determine lateral parameters in the middle of the step. |
---|
| 233 | // They depend on energy & position along step. |
---|
| 234 | // |
---|
| 235 | G4double Tau = ComputeTau(LongitudinalPosition); |
---|
| 236 | ComputeRadialParameters(Energy,Tau); |
---|
| 237 | } |
---|
| 238 | |
---|
| 239 | G4double Radius; |
---|
| 240 | G4double Random1 = G4UniformRand(); |
---|
| 241 | G4double Random2 = G4UniformRand(); |
---|
| 242 | |
---|
| 243 | if(Random1 <WeightCore) //WeightCore = p < w_i |
---|
| 244 | { |
---|
| 245 | Radius = Rm * RadiusCore * std::sqrt( Random2/(1. - Random2) ); |
---|
| 246 | } |
---|
| 247 | else |
---|
| 248 | { |
---|
| 249 | Radius = Rm * RadiusTail * std::sqrt( Random2/(1. - Random2) ); |
---|
| 250 | } |
---|
| 251 | Radius = std::min(Radius,DBL_MAX); |
---|
| 252 | return Radius; |
---|
| 253 | } |
---|
| 254 | |
---|
| 255 | G4double GFlashHomoShowerParameterisation:: |
---|
| 256 | ComputeTau(G4double LongitudinalPosition) |
---|
| 257 | { |
---|
| 258 | G4double tau = LongitudinalPosition / Tmaxh / X0 //<t> = T* a /(a - 1) |
---|
| 259 | * (Alphah-1.00) /Alphah * |
---|
| 260 | std::exp(AveLogAlphah)/(std::exp(AveLogAlphah)-1.); //ok |
---|
| 261 | return tau; |
---|
| 262 | } |
---|
| 263 | |
---|
| 264 | void GFlashHomoShowerParameterisation:: |
---|
| 265 | ComputeRadialParameters(G4double Energy, G4double Tau) |
---|
| 266 | { |
---|
| 267 | G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV) ; //ok |
---|
| 268 | G4double z2 = ParRC3+ParRC4*Z ; //ok |
---|
| 269 | RadiusCore = z1 + z2 * Tau ; //ok |
---|
| 270 | |
---|
| 271 | G4double p1 = ParWC1+ParWC2*Z; //ok |
---|
| 272 | G4double p2 = ParWC3+ParWC4*Z; //ok |
---|
| 273 | G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok |
---|
| 274 | |
---|
| 275 | WeightCore = p1 * std::exp( (p2-Tau)/p3 - std::exp( (p2-Tau) /p3) ); //ok |
---|
| 276 | |
---|
| 277 | G4double k1 = ParRT1+ParRT2*Z; // ok |
---|
| 278 | G4double k2 = ParRT3; // ok |
---|
| 279 | G4double k3 = ParRT4; // ok |
---|
| 280 | G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok |
---|
| 281 | |
---|
| 282 | RadiusTail = k1*(std::exp(k3*(Tau-k2)) + |
---|
| 283 | std::exp(k4*(Tau-k2)) ); //ok |
---|
| 284 | } |
---|
| 285 | |
---|
| 286 | G4double GFlashHomoShowerParameterisation:: |
---|
| 287 | GenerateExponential(const G4double /* Energy */ ) |
---|
| 288 | { |
---|
| 289 | G4double ParExp1 = 9./7.*X0; |
---|
| 290 | G4double random = -ParExp1*CLHEP::RandExponential::shoot() ; |
---|
| 291 | return random; |
---|
| 292 | } |
---|