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 11.0.p2)


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