| 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: GFlashSamplingShowerParameterisation.cc,v 1.6 2006/06/29 19:14:20 gunter Exp $
>> 27 // GEANT4 tag $Name: geant4-08-03-patch-01 $
26 // 28 //
27 // 29 //
28 // ------------------------------------------- 30 // ------------------------------------------------------------
29 // GEANT 4 class implementation 31 // GEANT 4 class implementation
30 // 32 //
31 // ------- GFlashSamplingShowerParameteri 33 // ------- GFlashSamplingShowerParameterisation -------
32 // 34 //
33 // Authors: E.Barberio & Joanna Weng - 11.2005 35 // Authors: E.Barberio & Joanna Weng - 11.2005
34 // ------------------------------------------- 36 // ------------------------------------------------------------
35 37
36 #include <cmath> <<
37 <<
38 #include "GFlashSamplingShowerParameterisation <<
39 #include "GVFlashShowerParameterisation.hh" 38 #include "GVFlashShowerParameterisation.hh"
40 #include "G4SystemOfUnits.hh" << 39 #include "GFlashSamplingShowerParameterisation.hh"
>> 40 #include <cmath>
41 #include "Randomize.hh" 41 #include "Randomize.hh"
42 #include "G4ios.hh" 42 #include "G4ios.hh"
43 #include "G4Material.hh" 43 #include "G4Material.hh"
44 #include "G4MaterialTable.hh" 44 #include "G4MaterialTable.hh"
45 45
46 GFlashSamplingShowerParameterisation::GFlashSa << 46 GFlashSamplingShowerParameterisation::
47 G4Material* aMat1, G4Material* aMat2, G4doub << 47 GFlashSamplingShowerParameterisation(G4Material* aMat1, G4Material* aMat2,
48 GFlashSamplingShowerTuning* aPar) << 48 G4double dd1, G4double dd2,
49 : GVFlashShowerParameterisation(), << 49 GFlashSamplingShowerTuning* aPar)
50 ParAveT2(0.), << 50 : GVFlashShowerParameterisation()
51 ParSigLogT1(0.), << 51 {
52 ParSigLogT2(0.), << 52 if(!aPar) { thePar = new GFlashSamplingShowerTuning; }
53 ParSigLogA1(0.), << 53 else { thePar = aPar; }
54 ParSigLogA2(0.), << 54
55 ParRho1(0.), << 55 SetMaterial(aMat1,aMat2 );
56 ParRho2(0.), << 56 d1=dd1;
57 ParsAveA2(0.), << 57 d2=dd2;
58 AveLogAlphah(0.), <<
59 AveLogTmaxh(0.), <<
60 SigmaLogAlphah(0.), <<
61 SigmaLogTmaxh(0.), <<
62 Rhoh(0.), <<
63 Alphah(0.), <<
64 Tmaxh(0.), <<
65 Betah(0.), <<
66 AveLogAlpha(0.), <<
67 AveLogTmax(0.), <<
68 SigmaLogAlpha(0.), <<
69 SigmaLogTmax(0.), <<
70 Rho(0.), <<
71 Alpha(0.), <<
72 Tmax(0.), <<
73 Beta(0.) <<
74 { <<
75 if (!aPar) { <<
76 thePar = new GFlashSamplingShowerTuning; <<
77 } <<
78 else { <<
79 thePar = aPar; <<
80 } <<
81 <<
82 SetMaterial(aMat1, aMat2); <<
83 d1 = dd1; <<
84 d2 = dd2; <<
85 58
86 // Longitudinal Coefficients for a homogenio 59 // Longitudinal Coefficients for a homogenious calo
87 60
88 // shower max 61 // shower max
89 ParAveT1 = thePar->ParAveT1(); // ln (ln y << 62 ParAveT1 = thePar->ParAveT1(); // ln (ln y -0.812)
90 ParAveA1 = thePar->ParAveA1(); // ln a (0.8 << 63 ParAveA1 = thePar->ParAveA1(); // ln a (0.81 + (0.458 + 2.26/Z)ln y)
91 ParAveA2 = thePar->ParAveA2(); << 64 ParAveA2 = thePar->ParAveA2();
92 ParAveA3 = thePar->ParAveA3(); << 65 ParAveA3 = thePar->ParAveA3();
93 // Variance of shower max sampling <<
94 ParSigLogT1 = <<
95 thePar <<
96 ->ParsSigLogT1(); // Sigma T1 (-1.4 + 1 <<
97 ParSigLogT2 = thePar->ParsSigLogT2(); // le <<
98 // variance of 'alpha' <<
99 ParSigLogA1 = <<
100 thePar <<
101 ->ParSigLogA1(); // Sigma a (-0.58 + 0. <<
102 ParSigLogA2 = thePar->ParSigLogA2(); // lea <<
103 // correlation alpha%T <<
104 ParRho1 = thePar->ParRho1(); // Rho = 0.705 <<
105 ParRho2 = thePar->ParRho2(); // leaving Par <<
106 // Sampling 66 // Sampling
107 ParsAveT1 = thePar->ParsAveT1(); // T_sam = << 67 ParsAveT1 = thePar->ParsAveT1(); // T_sam = log(exp( log T_hom) + t1*Fs-1 + t2*(1-ehat));
108 ParsAveT2 = thePar->ParsAveT2(); << 68 ParsAveT2 = thePar->ParsAveT2();
109 ParsAveA1 = thePar->ParsAveA1(); << 69 ParsAveA1 = thePar->ParsAveA1();
110 // Variance of shower max sampling << 70 // Variance of shower max sampling
111 ParsSigLogT1 = thePar->ParsSigLogT1(); // S << 71 ParsSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-2.5 + 1.25 ln y)**-1
112 // w << 72 ParsSigLogT2 = thePar->ParSigLogT2();
113 ParsSigLogT2 = thePar->ParsSigLogT2(); <<
114 // variance of 'alpha' 73 // variance of 'alpha'
115 ParsSigLogA1 = thePar->ParsSigLogA1(); // S << 74 ParsSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.82 + 0.79 ln y)**-1
116 // w << 75 ParsSigLogA2 = thePar->ParSigLogA2();
117 ParsSigLogA2 = thePar->ParsSigLogA2(); <<
118 // correlation alpha%T 76 // correlation alpha%T
119 ParsRho1 = << 77 ParsRho1 = thePar->ParRho1(); // Rho = 0.784 -0.023 ln y
120 thePar->ParsRho1(); // Rho = 0.784 -0.023 << 78 ParsRho2 = thePar->ParRho2();
121 ParsRho2 = thePar->ParsRho2(); <<
122 79
123 // Radial Coefficients 80 // Radial Coefficients
124 // r_C (tau)= z_1 +z_2 tau 81 // r_C (tau)= z_1 +z_2 tau
125 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+st 82 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2))))
126 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 << 83 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E
127 ParRC2 = thePar->ParRC2(); << 84 ParRC2 = thePar->ParRC2();
128 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 << 85 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z
129 ParRC4 = thePar->ParRC4(); << 86 ParRC4 = thePar->ParRC4();
130 87
131 ParWC1 = thePar->ParWC1(); 88 ParWC1 = thePar->ParWC1();
132 ParWC2 = thePar->ParWC2(); 89 ParWC2 = thePar->ParWC2();
133 ParWC3 = thePar->ParWC3(); 90 ParWC3 = thePar->ParWC3();
134 ParWC4 = thePar->ParWC4(); 91 ParWC4 = thePar->ParWC4();
135 ParWC5 = thePar->ParWC5(); << 92 ParWC5 = thePar->ParWC5();
136 ParWC6 = thePar->ParWC6(); 93 ParWC6 = thePar->ParWC6();
137 ParRT1 = thePar->ParRT1(); 94 ParRT1 = thePar->ParRT1();
138 ParRT2 = thePar->ParRT2(); 95 ParRT2 = thePar->ParRT2();
139 ParRT3 = thePar->ParRT3(); 96 ParRT3 = thePar->ParRT3();
140 ParRT4 = thePar->ParRT4(); << 97 ParRT4 = thePar->ParRT4();
141 ParRT5 = thePar->ParRT5(); 98 ParRT5 = thePar->ParRT5();
142 ParRT6 = thePar->ParRT6(); 99 ParRT6 = thePar->ParRT6();
143 100
144 // additional sampling parameter << 101 //additional sampling parameter
145 ParsRC1 = thePar->ParsRC1(); << 102 ParsRC1= thePar->ParsRC1();
146 ParsRC2 = thePar->ParsRC2(); << 103 ParsRC2= thePar->ParsRC2();
147 ParsWC1 = thePar->ParsWC1(); << 104 ParsWC1= thePar->ParsWC1();
148 ParsWC2 = thePar->ParsWC2(); << 105 ParsWC2= thePar->ParsWC2();
149 ParsRT1 = thePar->ParsRT1(); << 106 ParsRT1= thePar->ParsRT1();
150 ParsRT2 = thePar->ParsRT2(); << 107 ParsRT2= thePar->ParsRT2();
151 108
152 // Coeff for fluctuedted radial profiles fo 109 // Coeff for fluctuedted radial profiles for a sampling media
153 ParsSpotT1 = thePar->ParSpotT1(); // T_spot << 110 ParsSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212)
154 ParsSpotT2 = thePar->ParSpotT2(); << 111 ParsSpotT2 = thePar->ParSpotT2();
155 ParsSpotA1 = thePar->ParSpotA1(); // a_spot << 112 ParsSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334)
156 ParsSpotA2 = thePar->ParSpotA2(); << 113 ParsSpotA2 = thePar->ParSpotA2();
157 ParsSpotN1 = thePar->ParSpotN1(); // N_Spot << 114 ParsSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876
158 ParsSpotN2 = thePar->ParSpotN2(); << 115 ParsSpotN2 = thePar->ParSpotN2();
159 SamplingResolution = thePar->SamplingResolut << 116 SamplingResolution = thePar->SamplingResolution();
160 ConstantResolution = thePar->ConstantResolut << 117 ConstantResolution = thePar->ConstantResolution();
161 NoiseResolution = thePar->NoiseResolution(); << 118 NoiseResolution = thePar->NoiseResolution();
162 119
163 // Inits 120 // Inits
164 NSpot = 0.00; << 121 NSpot = 0.00;
165 AlphaNSpot = 0.00; << 122 AlphaNSpot = 0.00;
166 TNSpot = 0.00; << 123 TNSpot = 0.00;
167 BetaNSpot = 0.00; << 124 BetaNSpot = 0.00;
168 RadiusCore = 0.00; << 125 RadiusCore = 0.00;
169 WeightCore = 0.00; << 126 WeightCore = 0.00;
170 RadiusTail = 0.00; << 127 RadiusTail = 0.00;
171 ComputeZAX0EFFetc(); 128 ComputeZAX0EFFetc();
172 129
173 G4cout << "/******************************** 130 G4cout << "/********************************************/ " << G4endl;
174 G4cout << " - GFlashSamplingShowerParameter 131 G4cout << " - GFlashSamplingShowerParameterisation::Constructor - " << G4endl;
175 G4cout << "/******************************** << 132 G4cout << "/********************************************/ " << G4endl;
176 } 133 }
177 134
178 // ------------------------------------------- 135 // ------------------------------------------------------------
179 136
180 GFlashSamplingShowerParameterisation::~GFlashS 137 GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation()
181 { << 138 {}
182 delete thePar; <<
183 } <<
184 139
185 // ------------------------------------------- 140 // ------------------------------------------------------------
186 141
187 void GFlashSamplingShowerParameterisation::Set << 142 void GFlashSamplingShowerParameterisation::
>> 143 SetMaterial(G4Material *mat1, G4Material *mat2)
188 { 144 {
189 const G4double Es = 21.2 * MeV; << 145 G4double Es = 21*MeV;
190 material1 = mat1; << 146 material1= mat1;
191 Z1 = GetEffZ(material1); 147 Z1 = GetEffZ(material1);
192 A1 = GetEffA(material1); 148 A1 = GetEffA(material1);
193 density1 = material1->GetDensity(); 149 density1 = material1->GetDensity();
194 X01 = material1->GetRadlen(); << 150 X01 = material1->GetRadlen();
195 Ec1 = 2.66 * std::pow((X01 * Z1 / A1), 1.1); << 151 Ec1 = 2.66 * std::pow((X01 * Z1 / A1),1.1);
196 // Ec1 = 610.0 * MeV / (Z1 + 1.24); << 152 Rm1 = X01*Es/Ec1;
197 Rm1 = X01 * Es / Ec1; <<
198 153
199 material2 = mat2; << 154 material2= mat2;
200 Z2 = GetEffZ(material2); 155 Z2 = GetEffZ(material2);
201 A2 = GetEffA(material2); 156 A2 = GetEffA(material2);
202 density2 = material2->GetDensity(); 157 density2 = material2->GetDensity();
203 X02 = material2->GetRadlen(); << 158 X02 = material2->GetRadlen();
204 Ec2 = 2.66 * std::pow((X02 * Z2 / A2), 1.1); << 159 Ec2 = 2.66 * std::pow((X02 * Z2 / A2),1.1);
205 // Ec2 = 610.0 * MeV / (Z2 + 1.24); << 160 Rm2 = X02*Es/Ec2;
206 Rm2 = X02 * Es / Ec2; << 161 // PrintMaterial();
207 // PrintMaterial(); <<
208 } 162 }
209 163
210 // ------------------------------------------- 164 // ------------------------------------------------------------
211 165
212 void GFlashSamplingShowerParameterisation::Com 166 void GFlashSamplingShowerParameterisation::ComputeZAX0EFFetc()
213 { 167 {
214 G4cout << "/************ ComputeZAX0EFFetc * 168 G4cout << "/************ ComputeZAX0EFFetc ************/" << G4endl;
215 G4cout << " - GFlashSamplingShowerParameter 169 G4cout << " - GFlashSamplingShowerParameterisation::Material - " << G4endl;
216 170
217 const G4double Es = 21.2 * MeV; << 171 G4double Es = 21*MeV; //constant
218 172
219 // material and geometry parameters for a sa 173 // material and geometry parameters for a sampling calorimeter
220 G4double denominator = (d1 * density1 + d2 * << 174 G4double denominator = (d1*density1 + d2*density2);
221 G4double W1 = (d1 * density1) / denominator; << 175 G4double W1 = (d1*density1) / denominator;
222 G4double W2 = (d2 * density2) / denominator; << 176 G4double W2 = (d2*density2)/denominator;
223 Zeff = (W1 * Z1) + (W2 * Z2); // X0*Es/Ec; << 177 Zeff = ( W1*Z2 ) + (W2*Z1); //X0*Es/Ec;
224 Aeff = (W1 * A1) + (W2 * A2); << 178 Aeff = ( W1*A1 ) + (W2*A2);
225 Rhoeff = ((d1 * density1) + (d2 * density2)) << 179 X0eff =(1/ (( W1 / X01) +( W2 / X02)));
226 X0eff = (W1 * Rhoeff) / (X01 * density1) + ( << 180 Rhoeff = ( (d1 *density1 ) + (d2 * density2 ))/G4double (d2 + d1 );
227 X0eff = 1. / X0eff; << 181 Rmeff = 1/ ((((W1*Ec1)/ X01) + ((W2* Ec2)/ X02) ) / Es ) ;
228 Rmeff = 1. / ((((W1 * Ec1) / X01) + ((W2 * E << 182 Eceff = X0eff *((W1*Ec1)/ X01 + (W2* Ec2)/ X02 );
229 Eceff = X0eff * ((W1 * Ec1) / X01 + (W2 * Ec << 183 Fs = X0eff/G4double ((d1/mm )+(d2/mm) );
230 Fs = X0eff / (d1 + d2); // --> was G4double << 184 ehat = (1. / (1+ 0.007*(Z1- Z2)));
231 // mm makes sense.. <<
232 ehat = (1. / (1 + 0.007 * (Z1 - Z2))); <<
233 185
234 G4cout << "W1= " << W1 << G4endl; << 186 G4cout << "W1= " << W1 << G4endl;
235 G4cout << "W2= " << W2 << G4endl; 187 G4cout << "W2= " << W2 << G4endl;
236 G4cout << "effective quantities Zeff = " << << 188 G4cout << "effective quantities Zeff = "<<Zeff<< G4endl;
237 G4cout << "effective quantities Aeff = " << << 189 G4cout << "effective quantities Aeff = "<<Aeff<< G4endl;
238 G4cout << "effective quantities Rhoeff = " < << 190 G4cout << "effective quantities Rhoeff = "<<Rhoeff/g *cm3<<" g/cm3" << G4endl;
239 G4cout << "effective quantities X0eff = " << << 191 G4cout << "effective quantities X0eff = "<<X0eff/cm <<" cm" << G4endl;
240 192
241 X0eff = X0eff * Rhoeff; 193 X0eff = X0eff * Rhoeff;
242 194
243 G4cout << "effective quantities X0eff = " << << 195 G4cout << "effective quantities X0eff = "<<X0eff/g*cm2 <<" g/cm2" << G4endl;
244 X0eff = X0eff / Rhoeff; << 196 X0eff = X0eff /Rhoeff;
245 G4cout << "effective quantities RMeff = " << << 197 G4cout << "effective quantities RMeff = "<<Rmeff/cm<<" cm" << G4endl;
246 Rmeff = Rmeff * Rhoeff; << 198 Rmeff = Rmeff* Rhoeff;
247 G4cout << "effective quantities RMeff = " << << 199 G4cout << "effective quantities RMeff = "<<Rmeff/g *cm2<<" g/cm2" << G4endl;
248 Rmeff = Rmeff / Rhoeff; << 200 Rmeff = Rmeff/ Rhoeff;
249 G4cout << "effective quantities Eceff = " << << 201 G4cout << "effective quantities Eceff = "<<Eceff/MeV<< " MeV"<< G4endl;
250 G4cout << "effective quantities Fs = " << Fs << 202 G4cout << "effective quantities Fs = "<<Fs<<G4endl;
251 G4cout << "effective quantities ehat = " << << 203 G4cout << "effective quantities ehat = "<<ehat<<G4endl;
252 G4cout << "/******************************** << 204 G4cout << "/********************************************/ " <<G4endl;
253 } 205 }
254 206
255 // ------------------------------------------- 207 // ------------------------------------------------------------
256 208
257 void GFlashSamplingShowerParameterisation::Gen << 209 void GFlashSamplingShowerParameterisation::
>> 210 GenerateLongitudinalProfile(G4double Energy)
258 { 211 {
259 if ((material1 == 0) || (material2 == 0)) { << 212 if ((material1==0) || (material2 ==0))
>> 213 {
260 G4Exception("GFlashSamplingShowerParameter 214 G4Exception("GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile()",
261 "InvalidSetup", FatalException 215 "InvalidSetup", FatalException, "No material initialized!");
262 } << 216 }
263 G4double y = Energy / Eceff; << 217 G4double y = Energy/Eceff;
264 ComputeLongitudinalParameters(y); << 218 ComputeLongitudinalParameters(y);
265 GenerateEnergyProfile(y); 219 GenerateEnergyProfile(y);
266 GenerateNSpotProfile(y); 220 GenerateNSpotProfile(y);
267 } 221 }
268 222
269 // ------------------------------------------- 223 // ------------------------------------------------------------
270 224
271 void GFlashSamplingShowerParameterisation::Com << 225 void
>> 226 GFlashSamplingShowerParameterisation::ComputeLongitudinalParameters(G4double y)
272 { 227 {
273 AveLogTmaxh = std::log(std::max(ParAveT1 + s << 228 AveLogTmaxh = std::log(std::max(ParAveT1 +std::log(y),0.1)); //ok
274 AveLogAlphah = << 229 AveLogAlphah = std::log(std::max(ParAveA1 + (ParAveA2+ParAveA3/Zeff)*std::log(y),.1)); //ok
275 std::log(std::max(ParAveA1 + (ParAveA2 + P << 230 //hom
276 // hom << 231 SigmaLogTmaxh = std::min(0.5,1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ); //ok
277 SigmaLogTmaxh = std::min(0.5, 1.00 / (ParSig << 232 SigmaLogAlphah = std::min(0.5,1.00/( ParSigLogA1 + ParSigLogA2*std::log(y))); //ok
278 SigmaLogAlphah = std::min(0.5, 1.00 / (ParSi << 233 Rhoh = ParRho1+ParRho2*std::log(y);//ok
279 Rhoh = ParRho1 + ParRho2 * std::log(y); // << 234 // if sampling
280 // if sampling << 235 AveLogTmax = std::max(0.1,std::log(std::exp(AveLogTmaxh)
281 AveLogTmax = << 236 + ParsAveT1/Fs + ParsAveT2*(1-ehat))); //ok
282 std::max(0.1, std::log(std::exp(AveLogTmax << 237 AveLogAlpha = std::max(0.1,std::log(std::exp(AveLogAlphah)
283 AveLogAlpha = std::max(0.1, std::log(std::ex << 238 + (ParsAveA1/Fs))); //ok
284 // 239 //
285 SigmaLogTmax = std::min(0.5, 1.00 / (ParsSig << 240 SigmaLogTmax = std::min(0.5,1.00/( ParsSigLogT1
286 SigmaLogAlpha = std::min(0.5, 1.00 / (ParsSi << 241 + ParsSigLogT2*std::log(y)) ); //ok
287 Rho = ParsRho1 + ParsRho2 * std::log(y); // << 242 SigmaLogAlpha = std::min(0.5,1.00/( ParsSigLogA1
288 << 243 + ParsSigLogA2*std::log(y))); //ok
289 if (0) { << 244 Rho = ParsRho1+ParsRho2*std::log(y); //ok
290 G4cout << " y = " << y << G4end <<
291 G4cout << " std::log(std::exp(AveLogTmaxh) <<
292 << " std::log(" << std::exp(AveLogT <<
293 << ParsAveT2 * (1 - ehat) << ") = " <<
294 << " std::log(" << std::exp(AveLogT <<
295 << ParsAveT2 << "*" << (1 - ehat) < <<
296 << " std::log(" << std::exp(AveLogT <<
297 << G4endl; <<
298 G4cout << " AveLogTmaxh " << AveLogTmax <<
299 G4cout << " AveLogAlphah " << AveLogAlph <<
300 G4cout << " SigmaLogTmaxh " << SigmaLogTm <<
301 G4cout << " 1.00/( ParSigLogT1 + ParSigLog <<
302 << (ParSigLogT1 + ParSigLogT2 * std <<
303 << ParSigLogT1 << " + " << ParSigLo <<
304 << ParSigLogT1 << " + " << ParSigLo <<
305 G4cout << " SigmaLogAlphah " << SigmaLogAl <<
306 G4cout << " Rhoh " << Rhoh << G4 <<
307 G4cout << " AveLogTmax " << AveLogTmax <<
308 G4cout << " AveLogAlpha " << AveLogAlph <<
309 G4cout << " SigmaLogTmax " << SigmaLogTm <<
310 G4cout << " SigmaLogAlpha " << SigmaLogAl <<
311 G4cout << " Rho " << Rho << G4e <<
312 } <<
313 } 245 }
314 246
315 // ------------------------------------------- 247 // ------------------------------------------------------------
316 248
317 void GFlashSamplingShowerParameterisation::Gen 249 void GFlashSamplingShowerParameterisation::GenerateEnergyProfile(G4double /* y */)
318 { << 250 {
319 G4double Correlation1 = std::sqrt((1 + Rho) << 251 G4double Correlation1 = std::sqrt((1+Rho)/2);
320 G4double Correlation2 = std::sqrt((1 - Rho) << 252 G4double Correlation2 = std::sqrt((1-Rho)/2);
321 G4double Correlation1h = std::sqrt((1 + Rhoh << 253 G4double Correlation1h = std::sqrt((1+Rhoh)/2);
322 G4double Correlation2h = std::sqrt((1 - Rhoh << 254 G4double Correlation2h = std::sqrt((1-Rhoh)/2);
323 G4double Random1 = G4RandGauss::shoot(); 255 G4double Random1 = G4RandGauss::shoot();
324 G4double Random2 = G4RandGauss::shoot(); 256 G4double Random2 = G4RandGauss::shoot();
325 257
326 Tmax = std::max( << 258 Tmax = std::max(1.,std::exp( AveLogTmax + SigmaLogTmax *
327 1., std::exp(AveLogTmax + SigmaLogTmax * ( << 259 (Correlation1*Random1 + Correlation2*Random2) ));
328 Alpha = std::max( << 260 Alpha = std::max(1.1,std::exp( AveLogAlpha + SigmaLogAlpha *
329 1.1, std::exp(AveLogAlpha + SigmaLogAlpha << 261 (Correlation1*Random1 - Correlation2*Random2) ));
330 Beta = (Alpha - 1.00) / Tmax; << 262 Beta = (Alpha-1.00)/Tmax;
331 // Parameters for Enenrgy Profile including << 263 //Parameters for Enenrgy Profile including correaltion and sigmas
332 Tmaxh = << 264 Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh *
333 std::exp(AveLogTmaxh + SigmaLogTmaxh * (Co << 265 (Correlation1h*Random1 + Correlation2h*Random2) );
334 Alphah = << 266 Alphah = std::exp( AveLogAlphah + SigmaLogAlphah *
335 std::exp(AveLogAlphah + SigmaLogAlphah * ( << 267 (Correlation1h*Random1 - Correlation2h*Random2) );
336 Betah = (Alphah - 1.00) / Tmaxh; << 268 Betah = (Alphah-1.00)/Tmaxh;
337 } 269 }
338 270
339 // ------------------------------------------- 271 // ------------------------------------------------------------
340 272
341 void GFlashSamplingShowerParameterisation::Gen 273 void GFlashSamplingShowerParameterisation::GenerateNSpotProfile(const G4double y)
342 { 274 {
343 TNSpot = Tmaxh * (ParsSpotT1 + ParsSpotT2 * << 275 TNSpot = Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff); //ok.
344 TNSpot = std::max(0.5, Tmaxh * (ParsSpotT1 + << 276 TNSpot = std::max(0.5,Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff));
345 AlphaNSpot = Alphah * (ParsSpotA1 + ParsSpot << 277 AlphaNSpot = Alphah * (ParsSpotA1+ParsSpotA2*Zeff);
346 BetaNSpot = (AlphaNSpot - 1.00) / TNSpot; / << 278 BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok
347 NSpot = ParsSpotN1 / SamplingResolution * st << 279 NSpot = ParsSpotN1 /SamplingResolution * std::pow(y*Eceff/GeV,ParsSpotN2 );
348 } 280 }
349 281
350 // ------------------------------------------- 282 // ------------------------------------------------------------
351 283
352 G4double GFlashSamplingShowerParameterisation: << 284 G4double
>> 285 GFlashSamplingShowerParameterisation::
>> 286 ApplySampling(const G4double DEne, const G4double )
353 { 287 {
354 G4double DEneFluctuated = DEne; 288 G4double DEneFluctuated = DEne;
355 G4double Resolution = std::pow(SamplingResol << 289 G4double Resolution = std::pow(SamplingResolution,2);
356 290
357 // +pow(NoiseResolution,2)/ //@@@@@@@ << 291 // +pow(NoiseResolution,2)/ //@@@@@@@@ FIXME
358 // Energy*(1.*MeV)+ 292 // Energy*(1.*MeV)+
359 // pow(ConstantResol 293 // pow(ConstantResolution,2)*
360 // Energy/(1.*MeV); 294 // Energy/(1.*MeV);
361 295
362 if (Resolution > 0.0 && DEne > 0.00) { << 296 if(Resolution >0.0 && DEne > 0.00)
363 // G4float x1 = DEne / Resolution; << 297 {
364 // G4float x2 = G4RandGamma::shoot(x1, 1.0 << 298 G4float x1=DEne/Resolution;
365 // DEneFluctuated = x2; << 299 G4float x2 = CLHEP::RandGamma::shoot(x1, 1.0)*Resolution;
366 G4double x1 = DEne / Resolution; << 300 DEneFluctuated=x2;
367 G4double x2 = 1.0 / Resolution; <<
368 DEneFluctuated = G4RandGamma::shoot(x1, x2 <<
369 } 301 }
370 return DEneFluctuated; 302 return DEneFluctuated;
371 } 303 }
372 304
373 // ------------------------------------------- 305 // ------------------------------------------------------------
374 306
375 G4double GFlashSamplingShowerParameterisation: 307 G4double GFlashSamplingShowerParameterisation::
376 IntegrateEneLongitudinal(G4double Longitudinal 308 IntegrateEneLongitudinal(G4double LongitudinalStep)
377 { 309 {
378 G4double LongitudinalStepInX0 = Longitudinal 310 G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
379 G4float x1= Betah*LongitudinalStepInX0; 311 G4float x1= Betah*LongitudinalStepInX0;
380 G4float x2= Alphah; 312 G4float x2= Alphah;
381 float x3 = gam(x1,x2); 313 float x3 = gam(x1,x2);
382 G4double DEne=x3; 314 G4double DEne=x3;
383 return DEne; 315 return DEne;
384 } 316 }
385 317
386 // ------------------------------------------- 318 // ------------------------------------------------------------
387 319
388 G4double GFlashSamplingShowerParameterisation: << 320 G4double GFlashSamplingShowerParameterisation::
>> 321 IntegrateNspLongitudinal(G4double LongitudinalStep)
389 { 322 {
390 G4double LongitudinalStepInX0 = Longitudinal << 323 G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
391 G4float x1 = BetaNSpot * LongitudinalStepInX << 324 G4float x1 = BetaNSpot*LongitudinalStepInX0;
392 G4float x2 = AlphaNSpot; 325 G4float x2 = AlphaNSpot;
393 G4float x3 = gam(x1, x2); << 326 G4float x3 = gam(x1,x2);
394 G4double DNsp = x3; 327 G4double DNsp = x3;
395 return DNsp; 328 return DNsp;
396 } 329 }
397 330
398 // ------------------------------------------- 331 // ------------------------------------------------------------
399 332
400 G4double GFlashSamplingShowerParameterisation: << 333 G4double GFlashSamplingShowerParameterisation::
401 << 334 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition)
402 { 335 {
403 if (ispot < 1) { << 336 if(ispot < 1)
>> 337 {
404 // Determine lateral parameters in the mid 338 // Determine lateral parameters in the middle of the step.
405 // They depend on energy & position along 339 // They depend on energy & position along step
406 // 340 //
407 G4double Tau = ComputeTau(LongitudinalPosi 341 G4double Tau = ComputeTau(LongitudinalPosition);
408 ComputeRadialParameters(Energy, Tau); << 342 ComputeRadialParameters(Energy,Tau);
409 } 343 }
410 << 344
411 G4double Radius; 345 G4double Radius;
412 G4double Random1 = G4UniformRand(); 346 G4double Random1 = G4UniformRand();
413 G4double Random2 = G4UniformRand(); << 347 G4double Random2 = G4UniformRand();
414 if (Random1 < WeightCore) // WeightCore = p << 348 if(Random1 <WeightCore) //WeightCore = p < w_i
415 { 349 {
416 Radius = Rmeff * RadiusCore * std::sqrt(Ra << 350 Radius = Rmeff * RadiusCore * std::sqrt( Random2/(1. - Random2) );
417 } 351 }
418 else { << 352 else
419 Radius = Rmeff * RadiusTail * std::sqrt(Ra << 353 {
420 } << 354 Radius = Rmeff * RadiusTail * std::sqrt( Random2/(1. - Random2) );
421 Radius = std::min(Radius, DBL_MAX); << 355 }
>> 356 Radius = std::min(Radius,DBL_MAX);
422 return Radius; 357 return Radius;
423 } 358 }
424 359
425 // ------------------------------------------- 360 // ------------------------------------------------------------
426 361
427 G4double GFlashSamplingShowerParameterisation: << 362 G4double
428 { << 363 GFlashSamplingShowerParameterisation::
429 G4double tau = LongitudinalPosition / Tmax / << 364 ComputeTau(G4double LongitudinalPosition)
430 * (Alpha - 1.00) / Alpha * st << 365 {
431 / (std::exp(AveLogAlpha) - 1. << 366 G4double tau = LongitudinalPosition / Tmax/ X0eff //<t> = T* a /(a - 1)
>> 367 * (Alpha-1.00) /Alpha
>> 368 * std::exp(AveLogAlpha)/(std::exp(AveLogAlpha)-1.); //ok
432 return tau; 369 return tau;
433 } 370 }
434 371
435 // ------------------------------------------- 372 // ------------------------------------------------------------
436 373
437 void GFlashSamplingShowerParameterisation::Com << 374 void GFlashSamplingShowerParameterisation::
>> 375 ComputeRadialParameters(G4double Energy, G4double Tau)
438 { 376 {
439 G4double z1 = ParRC1 + ParRC2 * std::log(Ene << 377 G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV); //ok
440 G4double z2 = ParRC3 + ParRC4 * Zeff; // ok << 378 G4double z2 = ParRC3+ParRC4*Zeff; //ok
441 RadiusCore = z1 + z2 * Tau; // ok << 379 RadiusCore = z1 + z2 * Tau; //ok
442 G4double p1 = ParWC1 + ParWC2 * Zeff; // ok << 380 G4double p1 = ParWC1+ParWC2*Zeff; //ok
443 G4double p2 = ParWC3 + ParWC4 * Zeff; // ok << 381 G4double p2 = ParWC3+ParWC4*Zeff; //ok
444 G4double p3 = ParWC5 + ParWC6 * std::log(Ene << 382 G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok
445 WeightCore = p1 * std::exp((p2 - Tau) / p3 - << 383 WeightCore = p1 * std::exp( (p2-Tau)/p3- std::exp( (p2-Tau) /p3) ); //ok
446 << 384
447 G4double k1 = ParRT1 + ParRT2 * Zeff; // ok << 385 G4double k1 = ParRT1+ParRT2*Zeff; // ok
448 G4double k2 = ParRT3; // ok << 386 G4double k2 = ParRT3; // ok
449 G4double k3 = ParRT4; // ok << 387 G4double k3 = ParRT4; // ok
450 G4double k4 = ParRT5 + ParRT6 * std::log(Ene << 388 G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok
451 << 389
452 RadiusTail = k1 * (std::exp(k3 * (Tau - k2)) << 390 RadiusTail = k1*(std::exp(k3*(Tau-k2))
453 << 391 + std::exp(k4*(Tau-k2)) ); //ok
454 // sampling calorimeter << 392
455 << 393 // sampling calorimeter
456 RadiusCore = RadiusCore + ParsRC1 * (1 - eha << 394
457 WeightCore = << 395 RadiusCore = RadiusCore + ParsRC1*(1-ehat) + ParsRC2/Fs*std::exp(-Tau); //ok
458 WeightCore + (1 - ehat) * (ParsWC1 + ParsW << 396 WeightCore = WeightCore + (1-ehat)
459 RadiusTail = RadiusTail + (1 - ehat) * ParsR << 397 * (ParsWC1+ParsWC2/Fs * std::exp(-std::pow((Tau-1.),2))); //ok
>> 398 RadiusTail = RadiusTail + (1-ehat)* ParsRT1+ ParsRT2/Fs *std::exp(-Tau); //ok
460 } 399 }
461 400
462 // ------------------------------------------- 401 // ------------------------------------------------------------
463 402
464 G4double GFlashSamplingShowerParameterisation: 403 G4double GFlashSamplingShowerParameterisation::
465 GenerateExponential(const G4double /* Energy * 404 GenerateExponential(const G4double /* Energy */ )
466 { 405 {
467 G4double ParExp1 = 9./7.*X0eff; 406 G4double ParExp1 = 9./7.*X0eff;
468 G4double random = -ParExp1*G4RandExponentia << 407 G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
469 return random; 408 return random;
470 } 409 }
471 410