00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <cmath>
00038
00039 #include "GFlashHomoShowerParameterisation.hh"
00040 #include "GVFlashShowerParameterisation.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "Randomize.hh"
00044 #include "G4ios.hh"
00045 #include "G4Material.hh"
00046 #include "G4MaterialTable.hh"
00047
00048 GFlashHomoShowerParameterisation::
00049 GFlashHomoShowerParameterisation(G4Material * aMat,
00050 GVFlashHomoShowerTuning * aPar)
00051 : GVFlashShowerParameterisation(),
00052 ConstantResolution(0.), NoiseResolution(0.), SamplingResolution(0.),
00053 AveLogAlphah(0.), AveLogTmaxh(0.), SigmaLogAlphah(0.), SigmaLogTmaxh(0.),
00054 Rhoh(0.), Alphah(0.), Tmaxh(0.), Betah(0.)
00055
00056 {
00057 if(!aPar) { thePar = new GVFlashHomoShowerTuning; owning = true; }
00058 else { thePar = aPar; owning = false; }
00059
00060 SetMaterial(aMat);
00061 PrintMaterial(aMat);
00062
00063
00064
00065
00066
00067
00068
00069 ParAveT1 = thePar->ParAveT1();
00070 ParAveA1 = thePar->ParAveA1();
00071 ParAveA2 = thePar->ParAveA2();
00072 ParAveA3 = thePar->ParAveA3();
00073
00074
00075 ParSigLogT1 = thePar->ParSigLogT1();
00076 ParSigLogT2 = thePar->ParSigLogT2();
00077
00078
00079
00080 ParSigLogA1 = thePar->ParSigLogA1();
00081 ParSigLogA2 = thePar->ParSigLogA2();
00082
00083
00084
00085 ParRho1 = thePar->ParRho1();
00086 ParRho2 = thePar->ParRho2();
00087
00088
00089
00090
00091
00092 ParRC1 = thePar->ParRC1();
00093 ParRC2 = thePar->ParRC2();
00094
00095 ParRC3 = thePar->ParRC3();
00096 ParRC4 = thePar->ParRC4();
00097
00098 ParWC1 = thePar->ParWC1();
00099 ParWC2 = thePar->ParWC2();
00100 ParWC3 = thePar->ParWC3();
00101 ParWC4 = thePar->ParWC4();
00102 ParWC5 = thePar->ParWC5();
00103 ParWC6 = thePar->ParWC6();
00104
00105 ParRT1 = thePar->ParRT1();
00106 ParRT2 = thePar->ParRT2();
00107 ParRT3 = thePar->ParRT3();
00108 ParRT4 = thePar->ParRT4();
00109 ParRT5 = thePar->ParRT5();
00110 ParRT6 = thePar->ParRT6();
00111
00112
00113
00114 ParSpotT1 = thePar->ParSpotT1();
00115 ParSpotT2 = thePar->ParSpotT2();
00116
00117 ParSpotA1 = thePar->ParSpotA1();
00118 ParSpotA2 = thePar->ParSpotA2();
00119
00120 ParSpotN1 = thePar->ParSpotN1();
00121 ParSpotN2 = thePar->ParSpotN2();
00122
00123
00124
00125 NSpot = 0.00;
00126 AlphaNSpot = 0.00;
00127 TNSpot = 0.00;
00128 BetaNSpot = 0.00;
00129
00130 RadiusCore = 0.00;
00131 WeightCore = 0.00;
00132 RadiusTail = 0.00;
00133
00134 G4cout << "/********************************************/ " << G4endl;
00135 G4cout << " - GFlashHomoShowerParameterisation::Constructor - " << G4endl;
00136 G4cout << "/********************************************/ " << G4endl;
00137 }
00138
00139 void GFlashHomoShowerParameterisation::SetMaterial(G4Material *mat)
00140 {
00141 material= mat;
00142 Z = GetEffZ(material);
00143 A = GetEffA(material);
00144 density = material->GetDensity()/(g/cm3);
00145 X0 = material->GetRadlen();
00146 Ec = 2.66 * std::pow((X0 * Z / A),1.1);
00147 G4double Es = 21*MeV;
00148 Rm = X0*Es/Ec;
00149
00150 }
00151
00152 GFlashHomoShowerParameterisation::~GFlashHomoShowerParameterisation()
00153 {
00154 if(owning) { delete thePar; }
00155 }
00156
00157 void GFlashHomoShowerParameterisation::
00158 GenerateLongitudinalProfile(G4double Energy)
00159 {
00160 if (material==0)
00161 {
00162 G4Exception("GFlashHomoShowerParameterisation::GenerateLongitudinalProfile()",
00163 "InvalidSetup", FatalException, "No material initialized!");
00164 }
00165
00166 G4double y = Energy/Ec;
00167 ComputeLongitudinalParameters(y);
00168 GenerateEnergyProfile(y);
00169 GenerateNSpotProfile(y);
00170 }
00171
00172 void
00173 GFlashHomoShowerParameterisation::ComputeLongitudinalParameters(G4double y)
00174 {
00175 AveLogTmaxh = std::log(ParAveT1 + std::log(y));
00176
00177 AveLogAlphah = std::log(ParAveA1 + (ParAveA2+ParAveA3/Z)*std::log(y));
00178
00179
00180 SigmaLogTmaxh = 1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ;
00181
00182 SigmaLogAlphah = 1.00/( ParSigLogA1 + ParSigLogA2*std::log(y));
00183
00184 Rhoh = ParRho1+ParRho2*std::log(y);
00185 }
00186
00187 void GFlashHomoShowerParameterisation::GenerateEnergyProfile(G4double )
00188 {
00189 G4double Correlation1h = std::sqrt((1+Rhoh)/2);
00190 G4double Correlation2h = std::sqrt((1-Rhoh)/2);
00191
00192 G4double Random1 = G4RandGauss::shoot();
00193 G4double Random2 = G4RandGauss::shoot();
00194
00195
00196
00197 Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh *
00198 (Correlation1h*Random1 + Correlation2h*Random2) );
00199 Alphah = std::exp( AveLogAlphah + SigmaLogAlphah *
00200 (Correlation1h*Random1 - Correlation2h*Random2) );
00201 Betah = (Alphah-1.00)/Tmaxh;
00202 }
00203
00204 void GFlashHomoShowerParameterisation::GenerateNSpotProfile(const G4double y)
00205 {
00206 TNSpot = Tmaxh * (ParSpotT1+ParSpotT2*Z);
00207 AlphaNSpot = Alphah * (ParSpotA1+ParSpotA2*Z);
00208 BetaNSpot = (AlphaNSpot-1.00)/TNSpot;
00209 NSpot = ParSpotN1 * std::log(Z)*std::pow((y*Ec)/GeV,ParSpotN2 );
00210 }
00211
00212 G4double GFlashHomoShowerParameterisation::
00213 IntegrateEneLongitudinal(G4double LongitudinalStep)
00214 {
00215 G4double LongitudinalStepInX0 = LongitudinalStep / X0;
00216 G4float x1= Betah*LongitudinalStepInX0;
00217 G4float x2= Alphah;
00218 float x3 = gam(x1,x2);
00219 G4double DEne=x3;
00220 return DEne;
00221 }
00222
00223 G4double GFlashHomoShowerParameterisation::
00224 IntegrateNspLongitudinal(G4double LongitudinalStep)
00225 {
00226 G4double LongitudinalStepInX0 = LongitudinalStep / X0;
00227 G4float x1 = BetaNSpot*LongitudinalStepInX0;
00228 G4float x2 = AlphaNSpot;
00229 G4float x3 = gam(x1,x2);
00230 G4double DNsp = x3;
00231 return DNsp;
00232 }
00233
00234
00235 G4double GFlashHomoShowerParameterisation::
00236 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition)
00237 {
00238 if(ispot < 1)
00239 {
00240
00241
00242
00243 G4double Tau = ComputeTau(LongitudinalPosition);
00244 ComputeRadialParameters(Energy,Tau);
00245 }
00246
00247 G4double Radius;
00248 G4double Random1 = G4UniformRand();
00249 G4double Random2 = G4UniformRand();
00250
00251 if(Random1 <WeightCore)
00252 {
00253 Radius = Rm * RadiusCore * std::sqrt( Random2/(1. - Random2) );
00254 }
00255 else
00256 {
00257 Radius = Rm * RadiusTail * std::sqrt( Random2/(1. - Random2) );
00258 }
00259 Radius = std::min(Radius,DBL_MAX);
00260 return Radius;
00261 }
00262
00263 G4double GFlashHomoShowerParameterisation::
00264 ComputeTau(G4double LongitudinalPosition)
00265 {
00266 G4double tau = LongitudinalPosition / Tmaxh / X0
00267 * (Alphah-1.00) /Alphah *
00268 std::exp(AveLogAlphah)/(std::exp(AveLogAlphah)-1.);
00269 return tau;
00270 }
00271
00272 void GFlashHomoShowerParameterisation::
00273 ComputeRadialParameters(G4double Energy, G4double Tau)
00274 {
00275 G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV) ;
00276 G4double z2 = ParRC3+ParRC4*Z ;
00277 RadiusCore = z1 + z2 * Tau ;
00278
00279 G4double p1 = ParWC1+ParWC2*Z;
00280 G4double p2 = ParWC3+ParWC4*Z;
00281 G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV);
00282
00283 WeightCore = p1 * std::exp( (p2-Tau)/p3 - std::exp( (p2-Tau) /p3) );
00284
00285 G4double k1 = ParRT1+ParRT2*Z;
00286 G4double k2 = ParRT3;
00287 G4double k3 = ParRT4;
00288 G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV);
00289
00290 RadiusTail = k1*(std::exp(k3*(Tau-k2)) +
00291 std::exp(k4*(Tau-k2)) );
00292 }
00293
00294 G4double GFlashHomoShowerParameterisation::
00295 GenerateExponential(const G4double )
00296 {
00297 G4double ParExp1 = 9./7.*X0;
00298 G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
00299 return random;
00300 }