Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/parameterisations/gflash/src/GFlashSamplingShowerParameterisation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /parameterisations/gflash/src/GFlashSamplingShowerParameterisation.cc (Version 11.3.0) and /parameterisations/gflash/src/GFlashSamplingShowerParameterisation.cc (Version 9.4.p1)


  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-09-04-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