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: GFlashSamplingShowerParameterisation.cc,v 1.6 2006/06/29 19:14:20 gunter Exp $ |
---|
27 | // GEANT4 tag $Name: HEAD $ |
---|
28 | // |
---|
29 | // |
---|
30 | // ------------------------------------------------------------ |
---|
31 | // GEANT 4 class implementation |
---|
32 | // |
---|
33 | // ------- GFlashSamplingShowerParameterisation ------- |
---|
34 | // |
---|
35 | // Authors: E.Barberio & Joanna Weng - 11.2005 |
---|
36 | // ------------------------------------------------------------ |
---|
37 | |
---|
38 | #include "GVFlashShowerParameterisation.hh" |
---|
39 | #include "GFlashSamplingShowerParameterisation.hh" |
---|
40 | #include <cmath> |
---|
41 | #include "Randomize.hh" |
---|
42 | #include "G4ios.hh" |
---|
43 | #include "G4Material.hh" |
---|
44 | #include "G4MaterialTable.hh" |
---|
45 | |
---|
46 | GFlashSamplingShowerParameterisation:: |
---|
47 | GFlashSamplingShowerParameterisation(G4Material* aMat1, G4Material* aMat2, |
---|
48 | G4double dd1, G4double dd2, |
---|
49 | GFlashSamplingShowerTuning* aPar) |
---|
50 | : GVFlashShowerParameterisation() |
---|
51 | { |
---|
52 | if(!aPar) { thePar = new GFlashSamplingShowerTuning; } |
---|
53 | else { thePar = aPar; } |
---|
54 | |
---|
55 | SetMaterial(aMat1,aMat2 ); |
---|
56 | d1=dd1; |
---|
57 | d2=dd2; |
---|
58 | |
---|
59 | // Longitudinal Coefficients for a homogenious calo |
---|
60 | |
---|
61 | // shower max |
---|
62 | ParAveT1 = thePar->ParAveT1(); // ln (ln y -0.812) |
---|
63 | ParAveA1 = thePar->ParAveA1(); // ln a (0.81 + (0.458 + 2.26/Z)ln y) |
---|
64 | ParAveA2 = thePar->ParAveA2(); |
---|
65 | ParAveA3 = thePar->ParAveA3(); |
---|
66 | // Sampling |
---|
67 | ParsAveT1 = thePar->ParsAveT1(); // T_sam = log(exp( log T_hom) + t1*Fs-1 + t2*(1-ehat)); |
---|
68 | ParsAveT2 = thePar->ParsAveT2(); |
---|
69 | ParsAveA1 = thePar->ParsAveA1(); |
---|
70 | // Variance of shower max sampling |
---|
71 | ParsSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-2.5 + 1.25 ln y)**-1 |
---|
72 | ParsSigLogT2 = thePar->ParSigLogT2(); |
---|
73 | // variance of 'alpha' |
---|
74 | ParsSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.82 + 0.79 ln y)**-1 |
---|
75 | ParsSigLogA2 = thePar->ParSigLogA2(); |
---|
76 | // correlation alpha%T |
---|
77 | ParsRho1 = thePar->ParRho1(); // Rho = 0.784 -0.023 ln y |
---|
78 | ParsRho2 = thePar->ParRho2(); |
---|
79 | |
---|
80 | // Radial Coefficients |
---|
81 | // r_C (tau)= z_1 +z_2 tau |
---|
82 | // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2)))) |
---|
83 | ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E |
---|
84 | ParRC2 = thePar->ParRC2(); |
---|
85 | ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z |
---|
86 | ParRC4 = thePar->ParRC4(); |
---|
87 | |
---|
88 | ParWC1 = thePar->ParWC1(); |
---|
89 | ParWC2 = thePar->ParWC2(); |
---|
90 | ParWC3 = thePar->ParWC3(); |
---|
91 | ParWC4 = thePar->ParWC4(); |
---|
92 | ParWC5 = thePar->ParWC5(); |
---|
93 | ParWC6 = thePar->ParWC6(); |
---|
94 | ParRT1 = thePar->ParRT1(); |
---|
95 | ParRT2 = thePar->ParRT2(); |
---|
96 | ParRT3 = thePar->ParRT3(); |
---|
97 | ParRT4 = thePar->ParRT4(); |
---|
98 | ParRT5 = thePar->ParRT5(); |
---|
99 | ParRT6 = thePar->ParRT6(); |
---|
100 | |
---|
101 | //additional sampling parameter |
---|
102 | ParsRC1= thePar->ParsRC1(); |
---|
103 | ParsRC2= thePar->ParsRC2(); |
---|
104 | ParsWC1= thePar->ParsWC1(); |
---|
105 | ParsWC2= thePar->ParsWC2(); |
---|
106 | ParsRT1= thePar->ParsRT1(); |
---|
107 | ParsRT2= thePar->ParsRT2(); |
---|
108 | |
---|
109 | // Coeff for fluctuedted radial profiles for a sampling media |
---|
110 | ParsSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212) |
---|
111 | ParsSpotT2 = thePar->ParSpotT2(); |
---|
112 | ParsSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334) |
---|
113 | ParsSpotA2 = thePar->ParSpotA2(); |
---|
114 | ParsSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876 |
---|
115 | ParsSpotN2 = thePar->ParSpotN2(); |
---|
116 | SamplingResolution = thePar->SamplingResolution(); |
---|
117 | ConstantResolution = thePar->ConstantResolution(); |
---|
118 | NoiseResolution = thePar->NoiseResolution(); |
---|
119 | |
---|
120 | // Inits |
---|
121 | NSpot = 0.00; |
---|
122 | AlphaNSpot = 0.00; |
---|
123 | TNSpot = 0.00; |
---|
124 | BetaNSpot = 0.00; |
---|
125 | RadiusCore = 0.00; |
---|
126 | WeightCore = 0.00; |
---|
127 | RadiusTail = 0.00; |
---|
128 | ComputeZAX0EFFetc(); |
---|
129 | |
---|
130 | G4cout << "/********************************************/ " << G4endl; |
---|
131 | G4cout << " - GFlashSamplingShowerParameterisation::Constructor - " << G4endl; |
---|
132 | G4cout << "/********************************************/ " << G4endl; |
---|
133 | } |
---|
134 | |
---|
135 | // ------------------------------------------------------------ |
---|
136 | |
---|
137 | GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation() |
---|
138 | {} |
---|
139 | |
---|
140 | // ------------------------------------------------------------ |
---|
141 | |
---|
142 | void GFlashSamplingShowerParameterisation:: |
---|
143 | SetMaterial(G4Material *mat1, G4Material *mat2) |
---|
144 | { |
---|
145 | G4double Es = 21*MeV; |
---|
146 | material1= mat1; |
---|
147 | Z1 = GetEffZ(material1); |
---|
148 | A1 = GetEffA(material1); |
---|
149 | density1 = material1->GetDensity(); |
---|
150 | X01 = material1->GetRadlen(); |
---|
151 | Ec1 = 2.66 * std::pow((X01 * Z1 / A1),1.1); |
---|
152 | Rm1 = X01*Es/Ec1; |
---|
153 | |
---|
154 | material2= mat2; |
---|
155 | Z2 = GetEffZ(material2); |
---|
156 | A2 = GetEffA(material2); |
---|
157 | density2 = material2->GetDensity(); |
---|
158 | X02 = material2->GetRadlen(); |
---|
159 | Ec2 = 2.66 * std::pow((X02 * Z2 / A2),1.1); |
---|
160 | Rm2 = X02*Es/Ec2; |
---|
161 | // PrintMaterial(); |
---|
162 | } |
---|
163 | |
---|
164 | // ------------------------------------------------------------ |
---|
165 | |
---|
166 | void GFlashSamplingShowerParameterisation::ComputeZAX0EFFetc() |
---|
167 | { |
---|
168 | G4cout << "/************ ComputeZAX0EFFetc ************/" << G4endl; |
---|
169 | G4cout << " - GFlashSamplingShowerParameterisation::Material - " << G4endl; |
---|
170 | |
---|
171 | G4double Es = 21*MeV; //constant |
---|
172 | |
---|
173 | // material and geometry parameters for a sampling calorimeter |
---|
174 | G4double denominator = (d1*density1 + d2*density2); |
---|
175 | G4double W1 = (d1*density1) / denominator; |
---|
176 | G4double W2 = (d2*density2)/denominator; |
---|
177 | Zeff = ( W1*Z2 ) + (W2*Z1); //X0*Es/Ec; |
---|
178 | Aeff = ( W1*A1 ) + (W2*A2); |
---|
179 | X0eff =(1/ (( W1 / X01) +( W2 / X02))); |
---|
180 | Rhoeff = ( (d1 *density1 ) + (d2 * density2 ))/G4double (d2 + d1 ); |
---|
181 | Rmeff = 1/ ((((W1*Ec1)/ X01) + ((W2* Ec2)/ X02) ) / Es ) ; |
---|
182 | Eceff = X0eff *((W1*Ec1)/ X01 + (W2* Ec2)/ X02 ); |
---|
183 | Fs = X0eff/G4double ((d1/mm )+(d2/mm) ); |
---|
184 | ehat = (1. / (1+ 0.007*(Z1- Z2))); |
---|
185 | |
---|
186 | G4cout << "W1= " << W1 << G4endl; |
---|
187 | G4cout << "W2= " << W2 << G4endl; |
---|
188 | G4cout << "effective quantities Zeff = "<<Zeff<< G4endl; |
---|
189 | G4cout << "effective quantities Aeff = "<<Aeff<< G4endl; |
---|
190 | G4cout << "effective quantities Rhoeff = "<<Rhoeff/g *cm3<<" g/cm3" << G4endl; |
---|
191 | G4cout << "effective quantities X0eff = "<<X0eff/cm <<" cm" << G4endl; |
---|
192 | |
---|
193 | X0eff = X0eff * Rhoeff; |
---|
194 | |
---|
195 | G4cout << "effective quantities X0eff = "<<X0eff/g*cm2 <<" g/cm2" << G4endl; |
---|
196 | X0eff = X0eff /Rhoeff; |
---|
197 | G4cout << "effective quantities RMeff = "<<Rmeff/cm<<" cm" << G4endl; |
---|
198 | Rmeff = Rmeff* Rhoeff; |
---|
199 | G4cout << "effective quantities RMeff = "<<Rmeff/g *cm2<<" g/cm2" << G4endl; |
---|
200 | Rmeff = Rmeff/ Rhoeff; |
---|
201 | G4cout << "effective quantities Eceff = "<<Eceff/MeV<< " MeV"<< G4endl; |
---|
202 | G4cout << "effective quantities Fs = "<<Fs<<G4endl; |
---|
203 | G4cout << "effective quantities ehat = "<<ehat<<G4endl; |
---|
204 | G4cout << "/********************************************/ " <<G4endl; |
---|
205 | } |
---|
206 | |
---|
207 | // ------------------------------------------------------------ |
---|
208 | |
---|
209 | void GFlashSamplingShowerParameterisation:: |
---|
210 | GenerateLongitudinalProfile(G4double Energy) |
---|
211 | { |
---|
212 | if ((material1==0) || (material2 ==0)) |
---|
213 | { |
---|
214 | G4Exception("GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile()", |
---|
215 | "InvalidSetup", FatalException, "No material initialized!"); |
---|
216 | } |
---|
217 | G4double y = Energy/Eceff; |
---|
218 | ComputeLongitudinalParameters(y); |
---|
219 | GenerateEnergyProfile(y); |
---|
220 | GenerateNSpotProfile(y); |
---|
221 | } |
---|
222 | |
---|
223 | // ------------------------------------------------------------ |
---|
224 | |
---|
225 | void |
---|
226 | GFlashSamplingShowerParameterisation::ComputeLongitudinalParameters(G4double y) |
---|
227 | { |
---|
228 | AveLogTmaxh = std::log(std::max(ParAveT1 +std::log(y),0.1)); //ok |
---|
229 | AveLogAlphah = std::log(std::max(ParAveA1 + (ParAveA2+ParAveA3/Zeff)*std::log(y),.1)); //ok |
---|
230 | //hom |
---|
231 | SigmaLogTmaxh = std::min(0.5,1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ); //ok |
---|
232 | SigmaLogAlphah = std::min(0.5,1.00/( ParSigLogA1 + ParSigLogA2*std::log(y))); //ok |
---|
233 | Rhoh = ParRho1+ParRho2*std::log(y);//ok |
---|
234 | // if sampling |
---|
235 | AveLogTmax = std::max(0.1,std::log(std::exp(AveLogTmaxh) |
---|
236 | + ParsAveT1/Fs + ParsAveT2*(1-ehat))); //ok |
---|
237 | AveLogAlpha = std::max(0.1,std::log(std::exp(AveLogAlphah) |
---|
238 | + (ParsAveA1/Fs))); //ok |
---|
239 | // |
---|
240 | SigmaLogTmax = std::min(0.5,1.00/( ParsSigLogT1 |
---|
241 | + ParsSigLogT2*std::log(y)) ); //ok |
---|
242 | SigmaLogAlpha = std::min(0.5,1.00/( ParsSigLogA1 |
---|
243 | + ParsSigLogA2*std::log(y))); //ok |
---|
244 | Rho = ParsRho1+ParsRho2*std::log(y); //ok |
---|
245 | } |
---|
246 | |
---|
247 | // ------------------------------------------------------------ |
---|
248 | |
---|
249 | void GFlashSamplingShowerParameterisation::GenerateEnergyProfile(G4double /* y */) |
---|
250 | { |
---|
251 | G4double Correlation1 = std::sqrt((1+Rho)/2); |
---|
252 | G4double Correlation2 = std::sqrt((1-Rho)/2); |
---|
253 | G4double Correlation1h = std::sqrt((1+Rhoh)/2); |
---|
254 | G4double Correlation2h = std::sqrt((1-Rhoh)/2); |
---|
255 | G4double Random1 = G4RandGauss::shoot(); |
---|
256 | G4double Random2 = G4RandGauss::shoot(); |
---|
257 | |
---|
258 | Tmax = std::max(1.,std::exp( AveLogTmax + SigmaLogTmax * |
---|
259 | (Correlation1*Random1 + Correlation2*Random2) )); |
---|
260 | Alpha = std::max(1.1,std::exp( AveLogAlpha + SigmaLogAlpha * |
---|
261 | (Correlation1*Random1 - Correlation2*Random2) )); |
---|
262 | Beta = (Alpha-1.00)/Tmax; |
---|
263 | //Parameters for Enenrgy Profile including correaltion and sigmas |
---|
264 | Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh * |
---|
265 | (Correlation1h*Random1 + Correlation2h*Random2) ); |
---|
266 | Alphah = std::exp( AveLogAlphah + SigmaLogAlphah * |
---|
267 | (Correlation1h*Random1 - Correlation2h*Random2) ); |
---|
268 | Betah = (Alphah-1.00)/Tmaxh; |
---|
269 | } |
---|
270 | |
---|
271 | // ------------------------------------------------------------ |
---|
272 | |
---|
273 | void GFlashSamplingShowerParameterisation::GenerateNSpotProfile(const G4double y) |
---|
274 | { |
---|
275 | TNSpot = Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff); //ok. |
---|
276 | TNSpot = std::max(0.5,Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff)); |
---|
277 | AlphaNSpot = Alphah * (ParsSpotA1+ParsSpotA2*Zeff); |
---|
278 | BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok |
---|
279 | NSpot = ParsSpotN1 /SamplingResolution * std::pow(y*Eceff/GeV,ParsSpotN2 ); |
---|
280 | } |
---|
281 | |
---|
282 | // ------------------------------------------------------------ |
---|
283 | |
---|
284 | G4double |
---|
285 | GFlashSamplingShowerParameterisation:: |
---|
286 | ApplySampling(const G4double DEne, const G4double ) |
---|
287 | { |
---|
288 | G4double DEneFluctuated = DEne; |
---|
289 | G4double Resolution = std::pow(SamplingResolution,2); |
---|
290 | |
---|
291 | // +pow(NoiseResolution,2)/ //@@@@@@@@ FIXME |
---|
292 | // Energy*(1.*MeV)+ |
---|
293 | // pow(ConstantResolution,2)* |
---|
294 | // Energy/(1.*MeV); |
---|
295 | |
---|
296 | if(Resolution >0.0 && DEne > 0.00) |
---|
297 | { |
---|
298 | G4float x1=DEne/Resolution; |
---|
299 | G4float x2 = CLHEP::RandGamma::shoot(x1, 1.0)*Resolution; |
---|
300 | DEneFluctuated=x2; |
---|
301 | } |
---|
302 | return DEneFluctuated; |
---|
303 | } |
---|
304 | |
---|
305 | // ------------------------------------------------------------ |
---|
306 | |
---|
307 | G4double GFlashSamplingShowerParameterisation:: |
---|
308 | IntegrateEneLongitudinal(G4double LongitudinalStep) |
---|
309 | { |
---|
310 | G4double LongitudinalStepInX0 = LongitudinalStep / X0eff; |
---|
311 | G4float x1= Betah*LongitudinalStepInX0; |
---|
312 | G4float x2= Alphah; |
---|
313 | float x3 = gam(x1,x2); |
---|
314 | G4double DEne=x3; |
---|
315 | return DEne; |
---|
316 | } |
---|
317 | |
---|
318 | // ------------------------------------------------------------ |
---|
319 | |
---|
320 | G4double GFlashSamplingShowerParameterisation:: |
---|
321 | IntegrateNspLongitudinal(G4double LongitudinalStep) |
---|
322 | { |
---|
323 | G4double LongitudinalStepInX0 = LongitudinalStep / X0eff; |
---|
324 | G4float x1 = BetaNSpot*LongitudinalStepInX0; |
---|
325 | G4float x2 = AlphaNSpot; |
---|
326 | G4float x3 = gam(x1,x2); |
---|
327 | G4double DNsp = x3; |
---|
328 | return DNsp; |
---|
329 | } |
---|
330 | |
---|
331 | // ------------------------------------------------------------ |
---|
332 | |
---|
333 | G4double GFlashSamplingShowerParameterisation:: |
---|
334 | GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition) |
---|
335 | { |
---|
336 | if(ispot < 1) |
---|
337 | { |
---|
338 | // Determine lateral parameters in the middle of the step. |
---|
339 | // They depend on energy & position along step |
---|
340 | // |
---|
341 | G4double Tau = ComputeTau(LongitudinalPosition); |
---|
342 | ComputeRadialParameters(Energy,Tau); |
---|
343 | } |
---|
344 | |
---|
345 | G4double Radius; |
---|
346 | G4double Random1 = G4UniformRand(); |
---|
347 | G4double Random2 = G4UniformRand(); |
---|
348 | if(Random1 <WeightCore) //WeightCore = p < w_i |
---|
349 | { |
---|
350 | Radius = Rmeff * RadiusCore * std::sqrt( Random2/(1. - Random2) ); |
---|
351 | } |
---|
352 | else |
---|
353 | { |
---|
354 | Radius = Rmeff * RadiusTail * std::sqrt( Random2/(1. - Random2) ); |
---|
355 | } |
---|
356 | Radius = std::min(Radius,DBL_MAX); |
---|
357 | return Radius; |
---|
358 | } |
---|
359 | |
---|
360 | // ------------------------------------------------------------ |
---|
361 | |
---|
362 | G4double |
---|
363 | GFlashSamplingShowerParameterisation:: |
---|
364 | ComputeTau(G4double LongitudinalPosition) |
---|
365 | { |
---|
366 | G4double tau = LongitudinalPosition / Tmax/ X0eff //<t> = T* a /(a - 1) |
---|
367 | * (Alpha-1.00) /Alpha |
---|
368 | * std::exp(AveLogAlpha)/(std::exp(AveLogAlpha)-1.); //ok |
---|
369 | return tau; |
---|
370 | } |
---|
371 | |
---|
372 | // ------------------------------------------------------------ |
---|
373 | |
---|
374 | void GFlashSamplingShowerParameterisation:: |
---|
375 | ComputeRadialParameters(G4double Energy, G4double Tau) |
---|
376 | { |
---|
377 | G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV); //ok |
---|
378 | G4double z2 = ParRC3+ParRC4*Zeff; //ok |
---|
379 | RadiusCore = z1 + z2 * Tau; //ok |
---|
380 | G4double p1 = ParWC1+ParWC2*Zeff; //ok |
---|
381 | G4double p2 = ParWC3+ParWC4*Zeff; //ok |
---|
382 | G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok |
---|
383 | WeightCore = p1 * std::exp( (p2-Tau)/p3- std::exp( (p2-Tau) /p3) ); //ok |
---|
384 | |
---|
385 | G4double k1 = ParRT1+ParRT2*Zeff; // ok |
---|
386 | G4double k2 = ParRT3; // ok |
---|
387 | G4double k3 = ParRT4; // ok |
---|
388 | G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok |
---|
389 | |
---|
390 | RadiusTail = k1*(std::exp(k3*(Tau-k2)) |
---|
391 | + std::exp(k4*(Tau-k2)) ); //ok |
---|
392 | |
---|
393 | // sampling calorimeter |
---|
394 | |
---|
395 | RadiusCore = RadiusCore + ParsRC1*(1-ehat) + ParsRC2/Fs*std::exp(-Tau); //ok |
---|
396 | WeightCore = WeightCore + (1-ehat) |
---|
397 | * (ParsWC1+ParsWC2/Fs * std::exp(-std::pow((Tau-1.),2))); //ok |
---|
398 | RadiusTail = RadiusTail + (1-ehat)* ParsRT1+ ParsRT2/Fs *std::exp(-Tau); //ok |
---|
399 | } |
---|
400 | |
---|
401 | // ------------------------------------------------------------ |
---|
402 | |
---|
403 | G4double GFlashSamplingShowerParameterisation:: |
---|
404 | GenerateExponential(const G4double /* Energy */ ) |
---|
405 | { |
---|
406 | G4double ParExp1 = 9./7.*X0eff; |
---|
407 | G4double random = -ParExp1*CLHEP::RandExponential::shoot() ; |
---|
408 | return random; |
---|
409 | } |
---|