source: trunk/source/parameterisations/gflash/src/GFlashSamplingShowerParameterisation.cc @ 1202

Last change on this file since 1202 was 1058, checked in by garnier, 15 years ago

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