Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GFlashSamplingShowerParameterisation.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: GFlashSamplingShowerParameterisation.cc 69796 2013-05-15 13:26:12Z gcosmo $
27 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class implementation
31 //
32 // ------- GFlashSamplingShowerParameterisation -------
33 //
34 // Authors: E.Barberio & Joanna Weng - 11.2005
35 // ------------------------------------------------------------
36 
37 #include <cmath>
38 
41 #include "G4SystemOfUnits.hh"
42 #include "Randomize.hh"
43 #include "G4ios.hh"
44 #include "G4Material.hh"
45 #include "G4MaterialTable.hh"
46 
49  G4double dd1, G4double dd2,
52  ParAveT2(0.), ParSigLogT1(0.), ParSigLogT2(0.),
53  ParSigLogA1(0.), ParSigLogA2(0.), ParRho1(0.), ParRho2(0.), ParsAveA2(0.),
54  AveLogAlphah(0.), AveLogTmaxh(0.), SigmaLogAlphah(0.), SigmaLogTmaxh(0.),
55  Rhoh(0.), Alphah(0.), Tmaxh(0.), Betah(0.), AveLogAlpha(0.), AveLogTmax(0.),
56  SigmaLogAlpha(0.), SigmaLogTmax(0.), Rho(0.), Alpha(0.), Tmax(0.), Beta(0.)
57 {
58  if(!aPar) { thePar = new GFlashSamplingShowerTuning; owning = true; }
59  else { thePar = aPar; owning = false; }
60 
61  SetMaterial(aMat1,aMat2 );
62  d1=dd1;
63  d2=dd2;
64 
65  // Longitudinal Coefficients for a homogenious calo
66 
67  // shower max
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  // Sampling
73  ParsAveT1 = thePar->ParsAveT1(); // T_sam = log(exp( log T_hom) + t1*Fs-1 + t2*(1-ehat));
74  ParsAveT2 = thePar->ParsAveT2();
75  ParsAveA1 = thePar->ParsAveA1();
76  // Variance of shower max sampling
77  ParsSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-2.5 + 1.25 ln y)**-1
78  ParsSigLogT2 = thePar->ParSigLogT2();
79  // variance of 'alpha'
80  ParsSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.82 + 0.79 ln y)**-1
81  ParsSigLogA2 = thePar->ParSigLogA2();
82  // correlation alpha%T
83  ParsRho1 = thePar->ParRho1(); // Rho = 0.784 -0.023 ln y
84  ParsRho2 = thePar->ParRho2();
85 
86  // Radial Coefficients
87  // r_C (tau)= z_1 +z_2 tau
88  // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2))))
89  ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E
90  ParRC2 = thePar->ParRC2();
91  ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z
92  ParRC4 = thePar->ParRC4();
93 
94  ParWC1 = thePar->ParWC1();
95  ParWC2 = thePar->ParWC2();
96  ParWC3 = thePar->ParWC3();
97  ParWC4 = thePar->ParWC4();
98  ParWC5 = thePar->ParWC5();
99  ParWC6 = thePar->ParWC6();
100  ParRT1 = thePar->ParRT1();
101  ParRT2 = thePar->ParRT2();
102  ParRT3 = thePar->ParRT3();
103  ParRT4 = thePar->ParRT4();
104  ParRT5 = thePar->ParRT5();
105  ParRT6 = thePar->ParRT6();
106 
107  //additional sampling parameter
108  ParsRC1= thePar->ParsRC1();
109  ParsRC2= thePar->ParsRC2();
110  ParsWC1= thePar->ParsWC1();
111  ParsWC2= thePar->ParsWC2();
112  ParsRT1= thePar->ParsRT1();
113  ParsRT2= thePar->ParsRT2();
114 
115  // Coeff for fluctuedted radial profiles for a sampling media
116  ParsSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212)
117  ParsSpotT2 = thePar->ParSpotT2();
118  ParsSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334)
119  ParsSpotA2 = thePar->ParSpotA2();
120  ParsSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876
121  ParsSpotN2 = thePar->ParSpotN2();
122  SamplingResolution = thePar->SamplingResolution();
123  ConstantResolution = thePar->ConstantResolution();
124  NoiseResolution = thePar->NoiseResolution();
125 
126  // Inits
127  NSpot = 0.00;
128  AlphaNSpot = 0.00;
129  TNSpot = 0.00;
130  BetaNSpot = 0.00;
131  RadiusCore = 0.00;
132  WeightCore = 0.00;
133  RadiusTail = 0.00;
135 
136  G4cout << "/********************************************/ " << G4endl;
137  G4cout << " - GFlashSamplingShowerParameterisation::Constructor - " << G4endl;
138  G4cout << "/********************************************/ " << G4endl;
139 }
140 
141 // ------------------------------------------------------------
142 
144 {
145  if(owning) { delete thePar; }
146 }
147 
148 // ------------------------------------------------------------
149 
152 {
153  G4double Es = 21*MeV;
154  material1= mat1;
155  Z1 = GetEffZ(material1);
156  A1 = GetEffA(material1);
157  density1 = material1->GetDensity();
158  X01 = material1->GetRadlen();
159  Ec1 = 2.66 * std::pow((X01 * Z1 / A1),1.1);
160  Rm1 = X01*Es/Ec1;
161 
162  material2= mat2;
163  Z2 = GetEffZ(material2);
164  A2 = GetEffA(material2);
165  density2 = material2->GetDensity();
166  X02 = material2->GetRadlen();
167  Ec2 = 2.66 * std::pow((X02 * Z2 / A2),1.1);
168  Rm2 = X02*Es/Ec2;
169  // PrintMaterial();
170 }
171 
172 // ------------------------------------------------------------
173 
175 {
176  G4cout << "/************ ComputeZAX0EFFetc ************/" << G4endl;
177  G4cout << " - GFlashSamplingShowerParameterisation::Material - " << G4endl;
178 
179  G4double Es = 21*MeV; //constant
180 
181  // material and geometry parameters for a sampling calorimeter
182  G4double denominator = (d1*density1 + d2*density2);
183  G4double W1 = (d1*density1) / denominator;
184  G4double W2 = (d2*density2)/denominator;
185  Zeff = ( W1*Z2 ) + (W2*Z1); //X0*Es/Ec;
186  Aeff = ( W1*A1 ) + (W2*A2);
187  X0eff =(1/ (( W1 / X01) +( W2 / X02)));
188  Rhoeff = ( (d1 *density1 ) + (d2 * density2 ))/G4double (d2 + d1 );
189  Rmeff = 1/ ((((W1*Ec1)/ X01) + ((W2* Ec2)/ X02) ) / Es ) ;
190  Eceff = X0eff *((W1*Ec1)/ X01 + (W2* Ec2)/ X02 );
191  Fs = X0eff/G4double ((d1/mm )+(d2/mm) );
192  ehat = (1. / (1+ 0.007*(Z1- Z2)));
193 
194  G4cout << "W1= " << W1 << G4endl;
195  G4cout << "W2= " << W2 << G4endl;
196  G4cout << "effective quantities Zeff = "<<Zeff<< G4endl;
197  G4cout << "effective quantities Aeff = "<<Aeff<< G4endl;
198  G4cout << "effective quantities Rhoeff = "<<Rhoeff/g *cm3<<" g/cm3" << G4endl;
199  G4cout << "effective quantities X0eff = "<<X0eff/cm <<" cm" << G4endl;
200 
201  X0eff = X0eff * Rhoeff;
202 
203  G4cout << "effective quantities X0eff = "<<X0eff/g*cm2 <<" g/cm2" << G4endl;
204  X0eff = X0eff /Rhoeff;
205  G4cout << "effective quantities RMeff = "<<Rmeff/cm<<" cm" << G4endl;
206  Rmeff = Rmeff* Rhoeff;
207  G4cout << "effective quantities RMeff = "<<Rmeff/g *cm2<<" g/cm2" << G4endl;
208  Rmeff = Rmeff/ Rhoeff;
209  G4cout << "effective quantities Eceff = "<<Eceff/MeV<< " MeV"<< G4endl;
210  G4cout << "effective quantities Fs = "<<Fs<<G4endl;
211  G4cout << "effective quantities ehat = "<<ehat<<G4endl;
212  G4cout << "/********************************************/ " <<G4endl;
213 }
214 
215 // ------------------------------------------------------------
216 
219 {
220  if ((material1==0) || (material2 ==0))
221  {
222  G4Exception("GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile()",
223  "InvalidSetup", FatalException, "No material initialized!");
224  }
225  G4double y = Energy/Eceff;
226  ComputeLongitudinalParameters(y);
227  GenerateEnergyProfile(y);
228  GenerateNSpotProfile(y);
229 }
230 
231 // ------------------------------------------------------------
232 
233 void
234 GFlashSamplingShowerParameterisation::ComputeLongitudinalParameters(G4double y)
235 {
236  AveLogTmaxh = std::log(std::max(ParAveT1 +std::log(y),0.1)); //ok
237  AveLogAlphah = std::log(std::max(ParAveA1 + (ParAveA2+ParAveA3/Zeff)*std::log(y),.1)); //ok
238  //hom
239  SigmaLogTmaxh = std::min(0.5,1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ); //ok
240  SigmaLogAlphah = std::min(0.5,1.00/( ParSigLogA1 + ParSigLogA2*std::log(y))); //ok
241  Rhoh = ParRho1+ParRho2*std::log(y);//ok
242  // if sampling
243  AveLogTmax = std::max(0.1,std::log(std::exp(AveLogTmaxh)
244  + ParsAveT1/Fs + ParsAveT2*(1-ehat))); //ok
245  AveLogAlpha = std::max(0.1,std::log(std::exp(AveLogAlphah)
246  + (ParsAveA1/Fs))); //ok
247  //
248  SigmaLogTmax = std::min(0.5,1.00/( ParsSigLogT1
249  + ParsSigLogT2*std::log(y)) ); //ok
250  SigmaLogAlpha = std::min(0.5,1.00/( ParsSigLogA1
251  + ParsSigLogA2*std::log(y))); //ok
252  Rho = ParsRho1+ParsRho2*std::log(y); //ok
253 }
254 
255 // ------------------------------------------------------------
256 
257 void GFlashSamplingShowerParameterisation::GenerateEnergyProfile(G4double /* y */)
258 {
259  G4double Correlation1 = std::sqrt((1+Rho)/2);
260  G4double Correlation2 = std::sqrt((1-Rho)/2);
261  G4double Correlation1h = std::sqrt((1+Rhoh)/2);
262  G4double Correlation2h = std::sqrt((1-Rhoh)/2);
263  G4double Random1 = G4RandGauss::shoot();
264  G4double Random2 = G4RandGauss::shoot();
265 
266  Tmax = std::max(1.,std::exp( AveLogTmax + SigmaLogTmax *
267  (Correlation1*Random1 + Correlation2*Random2) ));
268  Alpha = std::max(1.1,std::exp( AveLogAlpha + SigmaLogAlpha *
269  (Correlation1*Random1 - Correlation2*Random2) ));
270  Beta = (Alpha-1.00)/Tmax;
271  //Parameters for Enenrgy Profile including correaltion and sigmas
272  Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh *
273  (Correlation1h*Random1 + Correlation2h*Random2) );
274  Alphah = std::exp( AveLogAlphah + SigmaLogAlphah *
275  (Correlation1h*Random1 - Correlation2h*Random2) );
276  Betah = (Alphah-1.00)/Tmaxh;
277 }
278 
279 // ------------------------------------------------------------
280 
281 void GFlashSamplingShowerParameterisation::GenerateNSpotProfile(const G4double y)
282 {
283  TNSpot = Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff); //ok.
284  TNSpot = std::max(0.5,Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff));
285  AlphaNSpot = Alphah * (ParsSpotA1+ParsSpotA2*Zeff);
286  BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok
287  NSpot = ParsSpotN1 /SamplingResolution * std::pow(y*Eceff/GeV,ParsSpotN2 );
288 }
289 
290 // ------------------------------------------------------------
291 
292 G4double
294 ApplySampling(const G4double DEne, const G4double )
295 {
296  G4double DEneFluctuated = DEne;
297  G4double Resolution = std::pow(SamplingResolution,2);
298 
299  // +pow(NoiseResolution,2)/ //@@@@@@@@ FIXME
300  // Energy*(1.*MeV)+
301  // pow(ConstantResolution,2)*
302  // Energy/(1.*MeV);
303 
304  if(Resolution >0.0 && DEne > 0.00)
305  {
306  G4float x1=DEne/Resolution;
307  G4float x2 = CLHEP::RandGamma::shoot(x1, 1.0)*Resolution;
308  DEneFluctuated=x2;
309  }
310  return DEneFluctuated;
311 }
312 
313 // ------------------------------------------------------------
314 
317 {
318  G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
319  G4float x1= Betah*LongitudinalStepInX0;
320  G4float x2= Alphah;
321  float x3 = gam(x1,x2);
322  G4double DEne=x3;
323  return DEne;
324 }
325 
326 // ------------------------------------------------------------
327 
330 {
331  G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
332  G4float x1 = BetaNSpot*LongitudinalStepInX0;
333  G4float x2 = AlphaNSpot;
334  G4float x3 = gam(x1,x2);
335  G4double DNsp = x3;
336  return DNsp;
337 }
338 
339 // ------------------------------------------------------------
340 
342 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition)
343 {
344  if(ispot < 1)
345  {
346  // Determine lateral parameters in the middle of the step.
347  // They depend on energy & position along step
348  //
349  G4double Tau = ComputeTau(LongitudinalPosition);
350  ComputeRadialParameters(Energy,Tau);
351  }
352 
353  G4double Radius;
354  G4double Random1 = G4UniformRand();
355  G4double Random2 = G4UniformRand();
356  if(Random1 <WeightCore) //WeightCore = p < w_i
357  {
358  Radius = Rmeff * RadiusCore * std::sqrt( Random2/(1. - Random2) );
359  }
360  else
361  {
362  Radius = Rmeff * RadiusTail * std::sqrt( Random2/(1. - Random2) );
363  }
364  Radius = std::min(Radius,DBL_MAX);
365  return Radius;
366 }
367 
368 // ------------------------------------------------------------
369 
370 G4double
372 ComputeTau(G4double LongitudinalPosition)
373 {
374  G4double tau = LongitudinalPosition / Tmax/ X0eff //<t> = T* a /(a - 1)
375  * (Alpha-1.00) /Alpha
376  * std::exp(AveLogAlpha)/(std::exp(AveLogAlpha)-1.); //ok
377  return tau;
378 }
379 
380 // ------------------------------------------------------------
381 
384 {
385  G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV); //ok
386  G4double z2 = ParRC3+ParRC4*Zeff; //ok
387  RadiusCore = z1 + z2 * Tau; //ok
388  G4double p1 = ParWC1+ParWC2*Zeff; //ok
389  G4double p2 = ParWC3+ParWC4*Zeff; //ok
390  G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok
391  WeightCore = p1 * std::exp( (p2-Tau)/p3- std::exp( (p2-Tau) /p3) ); //ok
392 
393  G4double k1 = ParRT1+ParRT2*Zeff; // ok
394  G4double k2 = ParRT3; // ok
395  G4double k3 = ParRT4; // ok
396  G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok
397 
398  RadiusTail = k1*(std::exp(k3*(Tau-k2))
399  + std::exp(k4*(Tau-k2)) ); //ok
400 
401  // sampling calorimeter
402 
403  RadiusCore = RadiusCore + ParsRC1*(1-ehat) + ParsRC2/Fs*std::exp(-Tau); //ok
404  WeightCore = WeightCore + (1-ehat)
405  * (ParsWC1+ParsWC2/Fs * std::exp(-std::pow((Tau-1.),2))); //ok
406  RadiusTail = RadiusTail + (1-ehat)* ParsRT1+ ParsRT2/Fs *std::exp(-Tau); //ok
407 }
408 
409 // ------------------------------------------------------------
410 
412 GenerateExponential(const G4double /* Energy */ )
413 {
414  G4double ParExp1 = 9./7.*X0eff;
415  G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
416  return random;
417 }