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