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