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 $ |
---|
27 | // GEANT4 tag $Name: $ |
---|
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 | } |
---|