Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id$ 26 // 27 // 27 // 28 // 28 // ------------------------------------------- 29 // ------------------------------------------------------------ 29 // GEANT 4 class implementation 30 // GEANT 4 class implementation 30 // 31 // 31 // ------- GFlashHomoShowerParameterisati 32 // ------- GFlashHomoShowerParameterisation ------- 32 // 33 // 33 // Authors: E.Barberio & Joanna Weng - 9.11.20 34 // Authors: E.Barberio & Joanna Weng - 9.11.2004 34 // ------------------------------------------- 35 // ------------------------------------------------------------ 35 36 36 #include <cmath> 37 #include <cmath> 37 38 38 #include "GFlashHomoShowerParameterisation.hh" 39 #include "GFlashHomoShowerParameterisation.hh" 39 #include "GVFlashShowerParameterisation.hh" 40 #include "GVFlashShowerParameterisation.hh" 40 #include "G4PhysicalConstants.hh" 41 #include "G4PhysicalConstants.hh" 41 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" 42 #include "Randomize.hh" 43 #include "Randomize.hh" 43 #include "G4ios.hh" 44 #include "G4ios.hh" 44 #include "G4Material.hh" 45 #include "G4Material.hh" 45 #include "G4MaterialTable.hh" 46 #include "G4MaterialTable.hh" 46 47 47 GFlashHomoShowerParameterisation::GFlashHomoSh << 48 GFlashHomoShowerParameterisation:: 48 << 49 GFlashHomoShowerParameterisation(G4Material * aMat, >> 50 GVFlashHomoShowerTuning * aPar) 49 : GVFlashShowerParameterisation(), 51 : GVFlashShowerParameterisation(), 50 ConstantResolution(0.), << 52 ConstantResolution(0.), NoiseResolution(0.), SamplingResolution(0.), 51 NoiseResolution(0.), << 53 AveLogAlphah(0.), AveLogTmaxh(0.), SigmaLogAlphah(0.), SigmaLogTmaxh(0.), 52 SamplingResolution(0.), << 54 Rhoh(0.), Alphah(0.), Tmaxh(0.), Betah(0.) 53 AveLogAlphah(0.), << 55 54 AveLogTmaxh(0.), << 56 { 55 SigmaLogAlphah(0.), << 57 if(!aPar) { thePar = new GVFlashHomoShowerTuning; owning = true; } 56 SigmaLogTmaxh(0.), << 58 else { thePar = aPar; owning = false; } 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 59 70 SetMaterial(aMat); 60 SetMaterial(aMat); 71 PrintMaterial(aMat); 61 PrintMaterial(aMat); 72 62 73 /******************************************* 63 /********************************************/ 74 /* Homo Calorimeter 64 /* Homo Calorimeter */ 75 /******************************************* << 65 /********************************************/ 76 // Longitudinal Coefficients for a homogenio 66 // Longitudinal Coefficients for a homogenious calo 77 // shower max 67 // shower max 78 // 68 // 79 ParAveT1 = thePar->ParAveT1(); // ln (ln y << 69 ParAveT1 = thePar->ParAveT1(); // ln (ln y -0.812) 80 ParAveA1 = thePar->ParAveA1(); // ln a (0.8 << 70 ParAveA1 = thePar->ParAveA1(); // ln a (0.81 + (0.458 + 2.26/Z)ln y) 81 ParAveA2 = thePar->ParAveA2(); << 71 ParAveA2 = thePar->ParAveA2(); 82 ParAveA3 = thePar->ParAveA3(); << 72 ParAveA3 = thePar->ParAveA3(); 83 73 84 // Variance of shower max 74 // Variance of shower max 85 ParSigLogT1 = thePar->ParSigLogT1(); // Sig << 75 ParSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-1.4 + 1.26 ln y)**-1 86 ParSigLogT2 = thePar->ParSigLogT2(); 76 ParSigLogT2 = thePar->ParSigLogT2(); 87 77 88 // variance of 'alpha' 78 // variance of 'alpha' 89 // 79 // 90 ParSigLogA1 = thePar->ParSigLogA1(); // Sig << 80 ParSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.58 + 0.86 ln y)**-1 91 ParSigLogA2 = thePar->ParSigLogA2(); 81 ParSigLogA2 = thePar->ParSigLogA2(); 92 82 93 // correlation alpha%T 83 // correlation alpha%T 94 // 84 // 95 ParRho1 = thePar->ParRho1(); // Rho = 0.705 << 85 ParRho1 = thePar->ParRho1(); // Rho = 0.705 -0.023 ln y 96 ParRho2 = thePar->ParRho2(); << 86 ParRho2 = thePar->ParRho2(); 97 87 98 // Radial Coefficients 88 // Radial Coefficients 99 // r_C (tau)= z_1 +z_2 tau 89 // r_C (tau)= z_1 +z_2 tau 100 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+st 90 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2)))) 101 // 91 // 102 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 << 92 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E 103 ParRC2 = thePar->ParRC2(); << 93 ParRC2 = thePar->ParRC2(); 104 94 105 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 << 95 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z 106 ParRC4 = thePar->ParRC4(); << 96 ParRC4 = thePar->ParRC4(); 107 97 108 ParWC1 = thePar->ParWC1(); 98 ParWC1 = thePar->ParWC1(); 109 ParWC2 = thePar->ParWC2(); 99 ParWC2 = thePar->ParWC2(); 110 ParWC3 = thePar->ParWC3(); 100 ParWC3 = thePar->ParWC3(); 111 ParWC4 = thePar->ParWC4(); 101 ParWC4 = thePar->ParWC4(); 112 ParWC5 = thePar->ParWC5(); << 102 ParWC5 = thePar->ParWC5(); 113 ParWC6 = thePar->ParWC6(); 103 ParWC6 = thePar->ParWC6(); 114 104 115 ParRT1 = thePar->ParRT1(); 105 ParRT1 = thePar->ParRT1(); 116 ParRT2 = thePar->ParRT2(); 106 ParRT2 = thePar->ParRT2(); 117 ParRT3 = thePar->ParRT3(); 107 ParRT3 = thePar->ParRT3(); 118 ParRT4 = thePar->ParRT4(); << 108 ParRT4 = thePar->ParRT4(); 119 ParRT5 = thePar->ParRT5(); 109 ParRT5 = thePar->ParRT5(); 120 ParRT6 = thePar->ParRT6(); 110 ParRT6 = thePar->ParRT6(); 121 111 122 // Coeff for fluctueted radial profiles for 112 // Coeff for fluctueted radial profiles for a uniform media 123 // 113 // 124 ParSpotT1 = thePar->ParSpotT1(); // T_spot << 114 ParSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212) 125 ParSpotT2 = thePar->ParSpotT2(); << 115 ParSpotT2 = thePar->ParSpotT2(); 126 116 127 ParSpotA1 = thePar->ParSpotA1(); // a_spot= << 117 ParSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334) 128 ParSpotA2 = thePar->ParSpotA2(); << 118 ParSpotA2 = thePar->ParSpotA2(); 129 119 130 ParSpotN1 = thePar->ParSpotN1(); // N_Spot << 120 ParSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876 131 ParSpotN2 = thePar->ParSpotN2(); << 121 ParSpotN2 = thePar->ParSpotN2(); 132 122 133 // Inits 123 // Inits 134 124 135 NSpot = 0.00; << 125 NSpot = 0.00; 136 AlphaNSpot = 0.00; << 126 AlphaNSpot = 0.00; 137 TNSpot = 0.00; << 127 TNSpot = 0.00; 138 BetaNSpot = 0.00; << 128 BetaNSpot = 0.00; 139 << 129 140 RadiusCore = 0.00; << 130 RadiusCore = 0.00; 141 WeightCore = 0.00; << 131 WeightCore = 0.00; 142 RadiusTail = 0.00; << 132 RadiusTail = 0.00; 143 133 144 G4cout << "/******************************** 134 G4cout << "/********************************************/ " << G4endl; 145 G4cout << " - GFlashHomoShowerParameterisat << 135 G4cout << " - GFlashHomoShowerParameterisation::Constructor - " << G4endl; 146 G4cout << "/******************************** 136 G4cout << "/********************************************/ " << G4endl; 147 } 137 } 148 138 149 void GFlashHomoShowerParameterisation::SetMate << 139 void GFlashHomoShowerParameterisation::SetMaterial(G4Material *mat) 150 { 140 { 151 material = mat; << 141 material= mat; 152 Z = GetEffZ(material); 142 Z = GetEffZ(material); 153 A = GetEffA(material); 143 A = GetEffA(material); 154 density = material->GetDensity() / (g / cm3) << 144 density = material->GetDensity()/(g/cm3); 155 X0 = material->GetRadlen(); << 145 X0 = material->GetRadlen(); 156 // O. I. Dovzhenkko and A. A. Pommanskii << 146 Ec = 2.66 * std::pow((X0 * Z / A),1.1); 157 Ec = 2.66 * std::pow((X0 * Z / A), 1.1); << 147 G4double Es = 21*MeV; 158 // // Rossi appriximation << 148 Rm = X0*Es/Ec; 159 // Ec = 610.0 * MeV / (Z + 1.24); << 149 // PrintMaterial(); 160 const G4double Es = 21.2 * MeV; << 161 Rm = X0 * Es / Ec; << 162 // PrintMaterial(); << 163 } 150 } 164 151 165 GFlashHomoShowerParameterisation::~GFlashHomoS 152 GFlashHomoShowerParameterisation::~GFlashHomoShowerParameterisation() 166 { 153 { 167 delete thePar; << 154 if(owning) { delete thePar; } 168 } 155 } 169 156 170 void GFlashHomoShowerParameterisation::Generat << 157 void GFlashHomoShowerParameterisation:: >> 158 GenerateLongitudinalProfile(G4double Energy) 171 { 159 { 172 if (material == 0) { << 160 if (material==0) 173 G4Exception("GFlashHomoShowerParameterisat << 161 { 174 FatalException, "No material i << 162 G4Exception("GFlashHomoShowerParameterisation::GenerateLongitudinalProfile()", >> 163 "InvalidSetup", FatalException, "No material initialized!"); 175 } 164 } 176 << 165 177 G4double y = Energy / Ec; << 166 G4double y = Energy/Ec; 178 ComputeLongitudinalParameters(y); << 167 ComputeLongitudinalParameters(y); 179 GenerateEnergyProfile(y); 168 GenerateEnergyProfile(y); 180 GenerateNSpotProfile(y); 169 GenerateNSpotProfile(y); 181 } 170 } 182 171 183 void GFlashHomoShowerParameterisation::Compute << 172 void >> 173 GFlashHomoShowerParameterisation::ComputeLongitudinalParameters(G4double y) 184 { 174 { 185 AveLogTmaxh = std::log(ParAveT1 + std::log(y << 175 AveLogTmaxh = std::log(ParAveT1 + std::log(y)); 186 // ok <ln T hom> << 176 //ok <ln T hom> 187 AveLogAlphah = std::log(ParAveA1 + (ParAveA2 << 177 AveLogAlphah = std::log(ParAveA1 + (ParAveA2+ParAveA3/Z)*std::log(y)); 188 // ok <ln alpha hom> << 178 //ok <ln alpha hom> 189 << 179 190 SigmaLogTmaxh = 1.00 / (ParSigLogT1 + ParSig << 180 SigmaLogTmaxh = 1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ; 191 // ok sigma (ln T hom) << 181 //ok sigma (ln T hom) 192 SigmaLogAlphah = 1.00 / (ParSigLogA1 + ParSi << 182 SigmaLogAlphah = 1.00/( ParSigLogA1 + ParSigLogA2*std::log(y)); 193 // ok sigma (ln alpha hom) << 183 //ok sigma (ln alpha hom) 194 Rhoh = ParRho1 + ParRho2 * std::log(y); // << 184 Rhoh = ParRho1+ParRho2*std::log(y); //ok 195 } 185 } 196 186 197 void GFlashHomoShowerParameterisation::Generat 187 void GFlashHomoShowerParameterisation::GenerateEnergyProfile(G4double /* y */) 198 { << 188 { 199 G4double Correlation1h = std::sqrt((1 + Rhoh << 189 G4double Correlation1h = std::sqrt((1+Rhoh)/2); 200 G4double Correlation2h = std::sqrt((1 - Rhoh << 190 G4double Correlation2h = std::sqrt((1-Rhoh)/2); 201 191 202 G4double Random1 = G4RandGauss::shoot(); 192 G4double Random1 = G4RandGauss::shoot(); 203 G4double Random2 = G4RandGauss::shoot(); 193 G4double Random2 = G4RandGauss::shoot(); 204 194 205 // Parameters for Enenrgy Profile including << 195 // Parameters for Enenrgy Profile including correaltion and sigmas 206 << 196 207 Tmaxh = << 197 Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh * 208 std::exp(AveLogTmaxh + SigmaLogTmaxh * (Co << 198 (Correlation1h*Random1 + Correlation2h*Random2) ); 209 Alphah = << 199 Alphah = std::exp( AveLogAlphah + SigmaLogAlphah * 210 std::exp(AveLogAlphah + SigmaLogAlphah * ( << 200 (Correlation1h*Random1 - Correlation2h*Random2) ); 211 Betah = (Alphah - 1.00) / Tmaxh; << 201 Betah = (Alphah-1.00)/Tmaxh; 212 } 202 } 213 203 214 void GFlashHomoShowerParameterisation::Generat 204 void GFlashHomoShowerParameterisation::GenerateNSpotProfile(const G4double y) 215 { 205 { 216 TNSpot = Tmaxh * (ParSpotT1 + ParSpotT2 * Z) << 206 TNSpot = Tmaxh * (ParSpotT1+ParSpotT2*Z); // ok 217 AlphaNSpot = Alphah * (ParSpotA1 + ParSpotA2 << 207 AlphaNSpot = Alphah * (ParSpotA1+ParSpotA2*Z); 218 BetaNSpot = (AlphaNSpot - 1.00) / TNSpot; / << 208 BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok 219 NSpot = ParSpotN1 * std::log(Z) * std::pow(( << 209 NSpot = ParSpotN1 * std::log(Z)*std::pow((y*Ec)/GeV,ParSpotN2 ); // ok 220 } 210 } 221 211 222 G4double GFlashHomoShowerParameterisation:: 212 G4double GFlashHomoShowerParameterisation:: 223 IntegrateEneLongitudinal(G4double Longitudinal 213 IntegrateEneLongitudinal(G4double LongitudinalStep) 224 { 214 { 225 G4double LongitudinalStepInX0 = Longitudinal 215 G4double LongitudinalStepInX0 = LongitudinalStep / X0; 226 G4float x1= Betah*LongitudinalStepInX0; 216 G4float x1= Betah*LongitudinalStepInX0; 227 G4float x2= Alphah; 217 G4float x2= Alphah; 228 float x3 = gam(x1,x2); 218 float x3 = gam(x1,x2); 229 G4double DEne=x3; 219 G4double DEne=x3; 230 return DEne; 220 return DEne; 231 } 221 } 232 222 233 G4double GFlashHomoShowerParameterisation::Int << 223 G4double GFlashHomoShowerParameterisation:: >> 224 IntegrateNspLongitudinal(G4double LongitudinalStep) 234 { 225 { 235 G4double LongitudinalStepInX0 = Longitudinal << 226 G4double LongitudinalStepInX0 = LongitudinalStep / X0; 236 G4float x1 = BetaNSpot * LongitudinalStepInX << 227 G4float x1 = BetaNSpot*LongitudinalStepInX0; 237 G4float x2 = AlphaNSpot; 228 G4float x2 = AlphaNSpot; 238 G4float x3 = gam(x1, x2); << 229 G4float x3 = gam(x1,x2); 239 G4double DNsp = x3; 230 G4double DNsp = x3; 240 return DNsp; 231 return DNsp; 241 } 232 } 242 233 243 G4double GFlashHomoShowerParameterisation::Gen << 234 244 << 235 G4double GFlashHomoShowerParameterisation:: >> 236 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition) 245 { 237 { 246 if (ispot < 1) { << 238 if(ispot < 1) >> 239 { 247 // Determine lateral parameters in the mid 240 // Determine lateral parameters in the middle of the step. 248 // They depend on energy & position along 241 // They depend on energy & position along step. 249 // 242 // 250 G4double Tau = ComputeTau(LongitudinalPosi 243 G4double Tau = ComputeTau(LongitudinalPosition); 251 ComputeRadialParameters(Energy, Tau); << 244 ComputeRadialParameters(Energy,Tau); 252 } 245 } 253 246 254 G4double Radius; 247 G4double Radius; 255 G4double Random1 = G4UniformRand(); 248 G4double Random1 = G4UniformRand(); 256 G4double Random2 = G4UniformRand(); << 249 G4double Random2 = G4UniformRand(); 257 250 258 if (Random1 < WeightCore) // WeightCore = p << 251 if(Random1 <WeightCore) //WeightCore = p < w_i 259 { 252 { 260 Radius = Rm * RadiusCore * std::sqrt(Rando << 253 Radius = Rm * RadiusCore * std::sqrt( Random2/(1. - Random2) ); 261 } << 262 else { << 263 Radius = Rm * RadiusTail * std::sqrt(Rando << 264 } 254 } 265 Radius = std::min(Radius, DBL_MAX); << 255 else >> 256 { >> 257 Radius = Rm * RadiusTail * std::sqrt( Random2/(1. - Random2) ); >> 258 } >> 259 Radius = std::min(Radius,DBL_MAX); 266 return Radius; 260 return Radius; 267 } 261 } 268 262 269 G4double GFlashHomoShowerParameterisation::Com << 263 G4double GFlashHomoShowerParameterisation:: >> 264 ComputeTau(G4double LongitudinalPosition) 270 { 265 { 271 G4double tau = LongitudinalPosition / Tmaxh << 266 G4double tau = LongitudinalPosition / Tmaxh / X0 //<t> = T* a /(a - 1) 272 * (Alphah - 1.00) / Alphah * << 267 * (Alphah-1.00) /Alphah * 273 / (std::exp(AveLogAlphah) - 1 << 268 std::exp(AveLogAlphah)/(std::exp(AveLogAlphah)-1.); //ok 274 return tau; 269 return tau; 275 } 270 } 276 271 277 void GFlashHomoShowerParameterisation::Compute << 272 void GFlashHomoShowerParameterisation:: >> 273 ComputeRadialParameters(G4double Energy, G4double Tau) 278 { 274 { 279 G4double z1 = ParRC1 + ParRC2 * std::log(Ene << 275 G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV) ; //ok 280 G4double z2 = ParRC3 + ParRC4 * Z; // ok << 276 G4double z2 = ParRC3+ParRC4*Z ; //ok 281 RadiusCore = z1 + z2 * Tau; // ok << 277 RadiusCore = z1 + z2 * Tau ; //ok 282 << 278 283 G4double p1 = ParWC1 + ParWC2 * Z; // ok << 279 G4double p1 = ParWC1+ParWC2*Z; //ok 284 G4double p2 = ParWC3 + ParWC4 * Z; // ok << 280 G4double p2 = ParWC3+ParWC4*Z; //ok 285 G4double p3 = ParWC5 + ParWC6 * std::log(Ene << 281 G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok 286 << 282 287 WeightCore = p1 * std::exp((p2 - Tau) / p3 - << 283 WeightCore = p1 * std::exp( (p2-Tau)/p3 - std::exp( (p2-Tau) /p3) ); //ok 288 << 284 289 G4double k1 = ParRT1 + ParRT2 * Z; // ok << 285 G4double k1 = ParRT1+ParRT2*Z; // ok 290 G4double k2 = ParRT3; // ok << 286 G4double k2 = ParRT3; // ok 291 G4double k3 = ParRT4; // ok << 287 G4double k3 = ParRT4; // ok 292 G4double k4 = ParRT5 + ParRT6 * std::log(Ene << 288 G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok 293 289 294 RadiusTail = k1 * (std::exp(k3 * (Tau - k2)) << 290 RadiusTail = k1*(std::exp(k3*(Tau-k2)) + >> 291 std::exp(k4*(Tau-k2)) ); //ok 295 } 292 } 296 293 297 G4double GFlashHomoShowerParameterisation:: 294 G4double GFlashHomoShowerParameterisation:: 298 GenerateExponential(const G4double /* Energy * 295 GenerateExponential(const G4double /* Energy */ ) 299 { 296 { 300 G4double ParExp1 = 9./7.*X0; 297 G4double ParExp1 = 9./7.*X0; 301 G4double random = -ParExp1*G4RandExponentia << 298 G4double random = -ParExp1*CLHEP::RandExponential::shoot() ; 302 return random; 299 return random; 303 } 300 } 304 301