Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GFlashHomoShowerParameterisation.cc
Go to the documentation of this file.
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 69796 2013-05-15 13:26:12Z gcosmo $
27 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class implementation
31 //
32 // ------- GFlashHomoShowerParameterisation -------
33 //
34 // Authors: E.Barberio & Joanna Weng - 9.11.2004
35 // ------------------------------------------------------------
36 
37 #include <cmath>
38 
39 #include "GFlashHomoShowerParameterisation.hh"
40 #include "GVFlashShowerParameterisation.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh"
44 #include "G4ios.hh"
45 #include "G4Material.hh"
46 #include "G4MaterialTable.hh"
47 
48 GFlashHomoShowerParameterisation::
49 GFlashHomoShowerParameterisation(G4Material * aMat,
50  GVFlashHomoShowerTuning * aPar)
51  : GVFlashShowerParameterisation(),
52  ConstantResolution(0.), NoiseResolution(0.), SamplingResolution(0.),
53  AveLogAlphah(0.), AveLogTmaxh(0.), SigmaLogAlphah(0.), SigmaLogTmaxh(0.),
54  Rhoh(0.), Alphah(0.), Tmaxh(0.), Betah(0.)
55 
56 {
57  if(!aPar) { thePar = new GVFlashHomoShowerTuning; owning = true; }
58  else { thePar = aPar; owning = false; }
59 
60  SetMaterial(aMat);
61  PrintMaterial(aMat);
62 
63  /********************************************/
64  /* Homo Calorimeter */
65  /********************************************/
66  // Longitudinal Coefficients for a homogenious calo
67  // shower max
68  //
69  ParAveT1 = thePar->ParAveT1(); // ln (ln y -0.812)
70  ParAveA1 = thePar->ParAveA1(); // ln a (0.81 + (0.458 + 2.26/Z)ln y)
71  ParAveA2 = thePar->ParAveA2();
72  ParAveA3 = thePar->ParAveA3();
73 
74  // Variance of shower max
75  ParSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-1.4 + 1.26 ln y)**-1
76  ParSigLogT2 = thePar->ParSigLogT2();
77 
78  // variance of 'alpha'
79  //
80  ParSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.58 + 0.86 ln y)**-1
81  ParSigLogA2 = thePar->ParSigLogA2();
82 
83  // correlation alpha%T
84  //
85  ParRho1 = thePar->ParRho1(); // Rho = 0.705 -0.023 ln y
86  ParRho2 = thePar->ParRho2();
87 
88  // Radial Coefficients
89  // r_C (tau)= z_1 +z_2 tau
90  // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2))))
91  //
92  ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E
93  ParRC2 = thePar->ParRC2();
94 
95  ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z
96  ParRC4 = thePar->ParRC4();
97 
98  ParWC1 = thePar->ParWC1();
99  ParWC2 = thePar->ParWC2();
100  ParWC3 = thePar->ParWC3();
101  ParWC4 = thePar->ParWC4();
102  ParWC5 = thePar->ParWC5();
103  ParWC6 = thePar->ParWC6();
104 
105  ParRT1 = thePar->ParRT1();
106  ParRT2 = thePar->ParRT2();
107  ParRT3 = thePar->ParRT3();
108  ParRT4 = thePar->ParRT4();
109  ParRT5 = thePar->ParRT5();
110  ParRT6 = thePar->ParRT6();
111 
112  // Coeff for fluctueted radial profiles for a uniform media
113  //
114  ParSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212)
115  ParSpotT2 = thePar->ParSpotT2();
116 
117  ParSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334)
118  ParSpotA2 = thePar->ParSpotA2();
119 
120  ParSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876
121  ParSpotN2 = thePar->ParSpotN2();
122 
123  // Inits
124 
125  NSpot = 0.00;
126  AlphaNSpot = 0.00;
127  TNSpot = 0.00;
128  BetaNSpot = 0.00;
129 
130  RadiusCore = 0.00;
131  WeightCore = 0.00;
132  RadiusTail = 0.00;
133 
134  G4cout << "/********************************************/ " << G4endl;
135  G4cout << " - GFlashHomoShowerParameterisation::Constructor - " << G4endl;
136  G4cout << "/********************************************/ " << G4endl;
137 }
138 
139 void GFlashHomoShowerParameterisation::SetMaterial(G4Material *mat)
140 {
141  material= mat;
142  Z = GetEffZ(material);
143  A = GetEffA(material);
144  density = material->GetDensity()/(g/cm3);
145  X0 = material->GetRadlen();
146  Ec = 2.66 * std::pow((X0 * Z / A),1.1);
147  G4double Es = 21*MeV;
148  Rm = X0*Es/Ec;
149  // PrintMaterial();
150 }
151 
152 GFlashHomoShowerParameterisation::~GFlashHomoShowerParameterisation()
153 {
154  if(owning) { delete thePar; }
155 }
156 
157 void GFlashHomoShowerParameterisation::
158 GenerateLongitudinalProfile(G4double Energy)
159 {
160  if (material==0)
161  {
162  G4Exception("GFlashHomoShowerParameterisation::GenerateLongitudinalProfile()",
163  "InvalidSetup", FatalException, "No material initialized!");
164  }
165 
166  G4double y = Energy/Ec;
167  ComputeLongitudinalParameters(y);
168  GenerateEnergyProfile(y);
169  GenerateNSpotProfile(y);
170 }
171 
172 void
173 GFlashHomoShowerParameterisation::ComputeLongitudinalParameters(G4double y)
174 {
175  AveLogTmaxh = std::log(ParAveT1 + std::log(y));
176  //ok <ln T hom>
177  AveLogAlphah = std::log(ParAveA1 + (ParAveA2+ParAveA3/Z)*std::log(y));
178  //ok <ln alpha hom>
179 
180  SigmaLogTmaxh = 1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ;
181  //ok sigma (ln T hom)
182  SigmaLogAlphah = 1.00/( ParSigLogA1 + ParSigLogA2*std::log(y));
183  //ok sigma (ln alpha hom)
184  Rhoh = ParRho1+ParRho2*std::log(y); //ok
185 }
186 
187 void GFlashHomoShowerParameterisation::GenerateEnergyProfile(G4double /* y */)
188 {
189  G4double Correlation1h = std::sqrt((1+Rhoh)/2);
190  G4double Correlation2h = std::sqrt((1-Rhoh)/2);
191 
192  G4double Random1 = G4RandGauss::shoot();
193  G4double Random2 = G4RandGauss::shoot();
194 
195  // Parameters for Enenrgy Profile including correaltion and sigmas
196 
197  Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh *
198  (Correlation1h*Random1 + Correlation2h*Random2) );
199  Alphah = std::exp( AveLogAlphah + SigmaLogAlphah *
200  (Correlation1h*Random1 - Correlation2h*Random2) );
201  Betah = (Alphah-1.00)/Tmaxh;
202 }
203 
204 void GFlashHomoShowerParameterisation::GenerateNSpotProfile(const G4double y)
205 {
206  TNSpot = Tmaxh * (ParSpotT1+ParSpotT2*Z); // ok
207  AlphaNSpot = Alphah * (ParSpotA1+ParSpotA2*Z);
208  BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok
209  NSpot = ParSpotN1 * std::log(Z)*std::pow((y*Ec)/GeV,ParSpotN2 ); // ok
210 }
211 
212 G4double GFlashHomoShowerParameterisation::
213 IntegrateEneLongitudinal(G4double LongitudinalStep)
214 {
215  G4double LongitudinalStepInX0 = LongitudinalStep / X0;
216  G4float x1= Betah*LongitudinalStepInX0;
217  G4float x2= Alphah;
218  float x3 = gam(x1,x2);
219  G4double DEne=x3;
220  return DEne;
221 }
222 
223 G4double GFlashHomoShowerParameterisation::
224 IntegrateNspLongitudinal(G4double LongitudinalStep)
225 {
226  G4double LongitudinalStepInX0 = LongitudinalStep / X0;
227  G4float x1 = BetaNSpot*LongitudinalStepInX0;
228  G4float x2 = AlphaNSpot;
229  G4float x3 = gam(x1,x2);
230  G4double DNsp = x3;
231  return DNsp;
232 }
233 
234 
235 G4double GFlashHomoShowerParameterisation::
236 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition)
237 {
238  if(ispot < 1)
239  {
240  // Determine lateral parameters in the middle of the step.
241  // They depend on energy & position along step.
242  //
243  G4double Tau = ComputeTau(LongitudinalPosition);
244  ComputeRadialParameters(Energy,Tau);
245  }
246 
247  G4double Radius;
248  G4double Random1 = G4UniformRand();
249  G4double Random2 = G4UniformRand();
250 
251  if(Random1 <WeightCore) //WeightCore = p < w_i
252  {
253  Radius = Rm * RadiusCore * std::sqrt( Random2/(1. - Random2) );
254  }
255  else
256  {
257  Radius = Rm * RadiusTail * std::sqrt( Random2/(1. - Random2) );
258  }
259  Radius = std::min(Radius,DBL_MAX);
260  return Radius;
261 }
262 
263 G4double GFlashHomoShowerParameterisation::
264 ComputeTau(G4double LongitudinalPosition)
265 {
266  G4double tau = LongitudinalPosition / Tmaxh / X0 //<t> = T* a /(a - 1)
267  * (Alphah-1.00) /Alphah *
268  std::exp(AveLogAlphah)/(std::exp(AveLogAlphah)-1.); //ok
269  return tau;
270 }
271 
272 void GFlashHomoShowerParameterisation::
273 ComputeRadialParameters(G4double Energy, G4double Tau)
274 {
275  G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV) ; //ok
276  G4double z2 = ParRC3+ParRC4*Z ; //ok
277  RadiusCore = z1 + z2 * Tau ; //ok
278 
279  G4double p1 = ParWC1+ParWC2*Z; //ok
280  G4double p2 = ParWC3+ParWC4*Z; //ok
281  G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok
282 
283  WeightCore = p1 * std::exp( (p2-Tau)/p3 - std::exp( (p2-Tau) /p3) ); //ok
284 
285  G4double k1 = ParRT1+ParRT2*Z; // ok
286  G4double k2 = ParRT3; // ok
287  G4double k3 = ParRT4; // ok
288  G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok
289 
290  RadiusTail = k1*(std::exp(k3*(Tau-k2)) +
291  std::exp(k4*(Tau-k2)) ); //ok
292 }
293 
294 G4double GFlashHomoShowerParameterisation::
295 GenerateExponential(const G4double /* Energy */ )
296 {
297  G4double ParExp1 = 9./7.*X0;
298  G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
299  return random;
300 }
Generated on Sat May 25 2013 14:33:21 for Geant4 by   doxygen 1.8.4