source: trunk/source/parameterisations/gflash/src/GFlashHomoShowerParameterisation.cc@ 1049

Last change on this file since 1049 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 9.7 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: GFlashHomoShowerParameterisation.cc,v 1.5 2006/06/29 19:14:14 gunter Exp $
27// GEANT4 tag $Name: HEAD $
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
46GFlashHomoShowerParameterisation::
47GFlashHomoShowerParameterisation(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
133void 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
146GFlashHomoShowerParameterisation::~GFlashHomoShowerParameterisation()
147{}
148
149void GFlashHomoShowerParameterisation::
150GenerateLongitudinalProfile(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
164void
165GFlashHomoShowerParameterisation::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
179void 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
196void 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
204G4double GFlashHomoShowerParameterisation::
205IntegrateEneLongitudinal(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
215G4double GFlashHomoShowerParameterisation::
216IntegrateNspLongitudinal(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
227G4double GFlashHomoShowerParameterisation::
228GenerateRadius(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
255G4double GFlashHomoShowerParameterisation::
256ComputeTau(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
264void GFlashHomoShowerParameterisation::
265ComputeRadialParameters(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
286G4double GFlashHomoShowerParameterisation::
287GenerateExponential(const G4double /* Energy */ )
288{
289 G4double ParExp1 = 9./7.*X0;
290 G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
291 return random;
292}
Note: See TracBrowser for help on using the repository browser.