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