source: trunk/source/parameterisations/gflash/src/GFlashSamplingShowerParameterisation.cc@ 1128

Last change on this file since 1128 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 15.0 KB
Line 
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: geant4-09-02-ref-02 $
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
46GFlashSamplingShowerParameterisation::
47GFlashSamplingShowerParameterisation(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
137GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation()
138{}
139
140// ------------------------------------------------------------
141
142void GFlashSamplingShowerParameterisation::
143SetMaterial(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
166void 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
209void GFlashSamplingShowerParameterisation::
210GenerateLongitudinalProfile(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
225void
226GFlashSamplingShowerParameterisation::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
249void 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
273void 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
284G4double
285GFlashSamplingShowerParameterisation::
286ApplySampling(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
307G4double GFlashSamplingShowerParameterisation::
308IntegrateEneLongitudinal(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
320G4double GFlashSamplingShowerParameterisation::
321IntegrateNspLongitudinal(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
333G4double GFlashSamplingShowerParameterisation::
334GenerateRadius(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
362G4double
363GFlashSamplingShowerParameterisation::
364ComputeTau(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
374void GFlashSamplingShowerParameterisation::
375ComputeRadialParameters(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
403G4double GFlashSamplingShowerParameterisation::
404GenerateExponential(const G4double /* Energy */ )
405{
406 G4double ParExp1 = 9./7.*X0eff;
407 G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
408 return random;
409}
Note: See TracBrowser for help on using the repository browser.