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 
39 #include "GFlashSamplingShowerParameterisation.hh"
40 #include "GVFlashShowerParameterisation.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "Randomize.hh"
43 #include "G4ios.hh"
44 #include "G4Material.hh"
45 #include "G4MaterialTable.hh"
46 
47 GFlashSamplingShowerParameterisation::
48 GFlashSamplingShowerParameterisation(G4Material* aMat1, G4Material* aMat2,
49  G4double dd1, G4double dd2,
50  GFlashSamplingShowerTuning* aPar)
51  : GVFlashShowerParameterisation(),
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;
134  ComputeZAX0EFFetc();
135 
136  G4cout << "/********************************************/ " << G4endl;
137  G4cout << " - GFlashSamplingShowerParameterisation::Constructor - " << G4endl;
138  G4cout << "/********************************************/ " << G4endl;
139 }
140 
141 // ------------------------------------------------------------
142 
143 GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation()
144 {
145  if(owning) { delete thePar; }
146 }
147 
148 // ------------------------------------------------------------
149 
150 void GFlashSamplingShowerParameterisation::
151 SetMaterial(G4Material *mat1, G4Material *mat2)
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 
174 void GFlashSamplingShowerParameterisation::ComputeZAX0EFFetc()
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 
217 void GFlashSamplingShowerParameterisation::
218 GenerateLongitudinalProfile(G4double Energy)
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
293 GFlashSamplingShowerParameterisation::
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 
315 G4double GFlashSamplingShowerParameterisation::
316 IntegrateEneLongitudinal(G4double LongitudinalStep)
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 
328 G4double GFlashSamplingShowerParameterisation::
329 IntegrateNspLongitudinal(G4double LongitudinalStep)
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 
341 G4double GFlashSamplingShowerParameterisation::
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
371 GFlashSamplingShowerParameterisation::
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 
382 void GFlashSamplingShowerParameterisation::
383 ComputeRadialParameters(G4double Energy, G4double Tau)
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 
411 G4double GFlashSamplingShowerParameterisation::
412 GenerateExponential(const G4double /* Energy */ )
413 {
414  G4double ParExp1 = 9./7.*X0eff;
415  G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
416  return random;
417 }
Generated on Sat May 25 2013 14:33:21 for Geant4 by   doxygen 1.8.4