| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // -------------------------------------------
29 // GEANT 4 class implementation
30 //
31 // ------- GFlashHomoShowerParameterisati
32 //
33 // Authors: E.Barberio & Joanna Weng - 9.11.20
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::GFlashHomoSh
48
49 : GVFlashShowerParameterisation(),
50 ConstantResolution(0.),
51 NoiseResolution(0.),
52 SamplingResolution(0.),
53 AveLogAlphah(0.),
54 AveLogTmaxh(0.),
55 SigmaLogAlphah(0.),
56 SigmaLogTmaxh(0.),
57 Rhoh(0.),
58 Alphah(0.),
59 Tmaxh(0.),
60 Betah(0.)
61
62 {
63 if (!aPar) {
64 thePar = new GVFlashHomoShowerTuning;
65 }
66 else {
67 thePar = aPar;
68 }
69
70 SetMaterial(aMat);
71 PrintMaterial(aMat);
72
73 /*******************************************
74 /* Homo Calorimeter
75 /*******************************************
76 // Longitudinal Coefficients for a homogenio
77 // shower max
78 //
79 ParAveT1 = thePar->ParAveT1(); // ln (ln y
80 ParAveA1 = thePar->ParAveA1(); // ln a (0.8
81 ParAveA2 = thePar->ParAveA2();
82 ParAveA3 = thePar->ParAveA3();
83
84 // Variance of shower max
85 ParSigLogT1 = thePar->ParSigLogT1(); // Sig
86 ParSigLogT2 = thePar->ParSigLogT2();
87
88 // variance of 'alpha'
89 //
90 ParSigLogA1 = thePar->ParSigLogA1(); // Sig
91 ParSigLogA2 = thePar->ParSigLogA2();
92
93 // correlation alpha%T
94 //
95 ParRho1 = thePar->ParRho1(); // Rho = 0.705
96 ParRho2 = thePar->ParRho2();
97
98 // Radial Coefficients
99 // r_C (tau)= z_1 +z_2 tau
100 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+st
101 //
102 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251
103 ParRC2 = thePar->ParRC2();
104
105 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162
106 ParRC4 = thePar->ParRC4();
107
108 ParWC1 = thePar->ParWC1();
109 ParWC2 = thePar->ParWC2();
110 ParWC3 = thePar->ParWC3();
111 ParWC4 = thePar->ParWC4();
112 ParWC5 = thePar->ParWC5();
113 ParWC6 = thePar->ParWC6();
114
115 ParRT1 = thePar->ParRT1();
116 ParRT2 = thePar->ParRT2();
117 ParRT3 = thePar->ParRT3();
118 ParRT4 = thePar->ParRT4();
119 ParRT5 = thePar->ParRT5();
120 ParRT6 = thePar->ParRT6();
121
122 // Coeff for fluctueted radial profiles for
123 //
124 ParSpotT1 = thePar->ParSpotT1(); // T_spot
125 ParSpotT2 = thePar->ParSpotT2();
126
127 ParSpotA1 = thePar->ParSpotA1(); // a_spot=
128 ParSpotA2 = thePar->ParSpotA2();
129
130 ParSpotN1 = thePar->ParSpotN1(); // N_Spot
131 ParSpotN2 = thePar->ParSpotN2();
132
133 // Inits
134
135 NSpot = 0.00;
136 AlphaNSpot = 0.00;
137 TNSpot = 0.00;
138 BetaNSpot = 0.00;
139
140 RadiusCore = 0.00;
141 WeightCore = 0.00;
142 RadiusTail = 0.00;
143
144 G4cout << "/********************************
145 G4cout << " - GFlashHomoShowerParameterisat
146 G4cout << "/********************************
147 }
148
149 void GFlashHomoShowerParameterisation::SetMate
150 {
151 material = mat;
152 Z = GetEffZ(material);
153 A = GetEffA(material);
154 density = material->GetDensity() / (g / cm3)
155 X0 = material->GetRadlen();
156 // O. I. Dovzhenkko and A. A. Pommanskii
157 Ec = 2.66 * std::pow((X0 * Z / A), 1.1);
158 // // Rossi appriximation
159 // Ec = 610.0 * MeV / (Z + 1.24);
160 const G4double Es = 21.2 * MeV;
161 Rm = X0 * Es / Ec;
162 // PrintMaterial();
163 }
164
165 GFlashHomoShowerParameterisation::~GFlashHomoS
166 {
167 delete thePar;
168 }
169
170 void GFlashHomoShowerParameterisation::Generat
171 {
172 if (material == 0) {
173 G4Exception("GFlashHomoShowerParameterisat
174 FatalException, "No material i
175 }
176
177 G4double y = Energy / Ec;
178 ComputeLongitudinalParameters(y);
179 GenerateEnergyProfile(y);
180 GenerateNSpotProfile(y);
181 }
182
183 void GFlashHomoShowerParameterisation::Compute
184 {
185 AveLogTmaxh = std::log(ParAveT1 + std::log(y
186 // ok <ln T hom>
187 AveLogAlphah = std::log(ParAveA1 + (ParAveA2
188 // ok <ln alpha hom>
189
190 SigmaLogTmaxh = 1.00 / (ParSigLogT1 + ParSig
191 // ok sigma (ln T hom)
192 SigmaLogAlphah = 1.00 / (ParSigLogA1 + ParSi
193 // ok sigma (ln alpha hom)
194 Rhoh = ParRho1 + ParRho2 * std::log(y); //
195 }
196
197 void GFlashHomoShowerParameterisation::Generat
198 {
199 G4double Correlation1h = std::sqrt((1 + Rhoh
200 G4double Correlation2h = std::sqrt((1 - Rhoh
201
202 G4double Random1 = G4RandGauss::shoot();
203 G4double Random2 = G4RandGauss::shoot();
204
205 // Parameters for Enenrgy Profile including
206
207 Tmaxh =
208 std::exp(AveLogTmaxh + SigmaLogTmaxh * (Co
209 Alphah =
210 std::exp(AveLogAlphah + SigmaLogAlphah * (
211 Betah = (Alphah - 1.00) / Tmaxh;
212 }
213
214 void GFlashHomoShowerParameterisation::Generat
215 {
216 TNSpot = Tmaxh * (ParSpotT1 + ParSpotT2 * Z)
217 AlphaNSpot = Alphah * (ParSpotA1 + ParSpotA2
218 BetaNSpot = (AlphaNSpot - 1.00) / TNSpot; /
219 NSpot = ParSpotN1 * std::log(Z) * std::pow((
220 }
221
222 G4double GFlashHomoShowerParameterisation::
223 IntegrateEneLongitudinal(G4double Longitudinal
224 {
225 G4double LongitudinalStepInX0 = Longitudinal
226 G4float x1= Betah*LongitudinalStepInX0;
227 G4float x2= Alphah;
228 float x3 = gam(x1,x2);
229 G4double DEne=x3;
230 return DEne;
231 }
232
233 G4double GFlashHomoShowerParameterisation::Int
234 {
235 G4double LongitudinalStepInX0 = Longitudinal
236 G4float x1 = BetaNSpot * LongitudinalStepInX
237 G4float x2 = AlphaNSpot;
238 G4float x3 = gam(x1, x2);
239 G4double DNsp = x3;
240 return DNsp;
241 }
242
243 G4double GFlashHomoShowerParameterisation::Gen
244
245 {
246 if (ispot < 1) {
247 // Determine lateral parameters in the mid
248 // They depend on energy & position along
249 //
250 G4double Tau = ComputeTau(LongitudinalPosi
251 ComputeRadialParameters(Energy, Tau);
252 }
253
254 G4double Radius;
255 G4double Random1 = G4UniformRand();
256 G4double Random2 = G4UniformRand();
257
258 if (Random1 < WeightCore) // WeightCore = p
259 {
260 Radius = Rm * RadiusCore * std::sqrt(Rando
261 }
262 else {
263 Radius = Rm * RadiusTail * std::sqrt(Rando
264 }
265 Radius = std::min(Radius, DBL_MAX);
266 return Radius;
267 }
268
269 G4double GFlashHomoShowerParameterisation::Com
270 {
271 G4double tau = LongitudinalPosition / Tmaxh
272 * (Alphah - 1.00) / Alphah *
273 / (std::exp(AveLogAlphah) - 1
274 return tau;
275 }
276
277 void GFlashHomoShowerParameterisation::Compute
278 {
279 G4double z1 = ParRC1 + ParRC2 * std::log(Ene
280 G4double z2 = ParRC3 + ParRC4 * Z; // ok
281 RadiusCore = z1 + z2 * Tau; // ok
282
283 G4double p1 = ParWC1 + ParWC2 * Z; // ok
284 G4double p2 = ParWC3 + ParWC4 * Z; // ok
285 G4double p3 = ParWC5 + ParWC6 * std::log(Ene
286
287 WeightCore = p1 * std::exp((p2 - Tau) / p3 -
288
289 G4double k1 = ParRT1 + ParRT2 * Z; // ok
290 G4double k2 = ParRT3; // ok
291 G4double k3 = ParRT4; // ok
292 G4double k4 = ParRT5 + ParRT6 * std::log(Ene
293
294 RadiusTail = k1 * (std::exp(k3 * (Tau - k2))
295 }
296
297 G4double GFlashHomoShowerParameterisation::
298 GenerateExponential(const G4double /* Energy *
299 {
300 G4double ParExp1 = 9./7.*X0;
301 G4double random = -ParExp1*G4RandExponentia
302 return random;
303 }
304