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 "GFlashSamplingShowerParameterisation.hh"
00040 #include "GVFlashShowerParameterisation.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "Randomize.hh"
00043 #include "G4ios.hh"
00044 #include "G4Material.hh"
00045 #include "G4MaterialTable.hh"
00046
00047 GFlashSamplingShowerParameterisation::
00048 GFlashSamplingShowerParameterisation(G4Material* aMat1, G4Material* aMat2,
00049 G4double dd1, G4double dd2,
00050 GFlashSamplingShowerTuning* aPar)
00051 : GVFlashShowerParameterisation(),
00052 ParAveT2(0.), ParSigLogT1(0.), ParSigLogT2(0.),
00053 ParSigLogA1(0.), ParSigLogA2(0.), ParRho1(0.), ParRho2(0.), ParsAveA2(0.),
00054 AveLogAlphah(0.), AveLogTmaxh(0.), SigmaLogAlphah(0.), SigmaLogTmaxh(0.),
00055 Rhoh(0.), Alphah(0.), Tmaxh(0.), Betah(0.), AveLogAlpha(0.), AveLogTmax(0.),
00056 SigmaLogAlpha(0.), SigmaLogTmax(0.), Rho(0.), Alpha(0.), Tmax(0.), Beta(0.)
00057 {
00058 if(!aPar) { thePar = new GFlashSamplingShowerTuning; owning = true; }
00059 else { thePar = aPar; owning = false; }
00060
00061 SetMaterial(aMat1,aMat2 );
00062 d1=dd1;
00063 d2=dd2;
00064
00065
00066
00067
00068 ParAveT1 = thePar->ParAveT1();
00069 ParAveA1 = thePar->ParAveA1();
00070 ParAveA2 = thePar->ParAveA2();
00071 ParAveA3 = thePar->ParAveA3();
00072
00073 ParsAveT1 = thePar->ParsAveT1();
00074 ParsAveT2 = thePar->ParsAveT2();
00075 ParsAveA1 = thePar->ParsAveA1();
00076
00077 ParsSigLogT1 = thePar->ParSigLogT1();
00078 ParsSigLogT2 = thePar->ParSigLogT2();
00079
00080 ParsSigLogA1 = thePar->ParSigLogA1();
00081 ParsSigLogA2 = thePar->ParSigLogA2();
00082
00083 ParsRho1 = thePar->ParRho1();
00084 ParsRho2 = thePar->ParRho2();
00085
00086
00087
00088
00089 ParRC1 = thePar->ParRC1();
00090 ParRC2 = thePar->ParRC2();
00091 ParRC3 = thePar->ParRC3();
00092 ParRC4 = thePar->ParRC4();
00093
00094 ParWC1 = thePar->ParWC1();
00095 ParWC2 = thePar->ParWC2();
00096 ParWC3 = thePar->ParWC3();
00097 ParWC4 = thePar->ParWC4();
00098 ParWC5 = thePar->ParWC5();
00099 ParWC6 = thePar->ParWC6();
00100 ParRT1 = thePar->ParRT1();
00101 ParRT2 = thePar->ParRT2();
00102 ParRT3 = thePar->ParRT3();
00103 ParRT4 = thePar->ParRT4();
00104 ParRT5 = thePar->ParRT5();
00105 ParRT6 = thePar->ParRT6();
00106
00107
00108 ParsRC1= thePar->ParsRC1();
00109 ParsRC2= thePar->ParsRC2();
00110 ParsWC1= thePar->ParsWC1();
00111 ParsWC2= thePar->ParsWC2();
00112 ParsRT1= thePar->ParsRT1();
00113 ParsRT2= thePar->ParsRT2();
00114
00115
00116 ParsSpotT1 = thePar->ParSpotT1();
00117 ParsSpotT2 = thePar->ParSpotT2();
00118 ParsSpotA1 = thePar->ParSpotA1();
00119 ParsSpotA2 = thePar->ParSpotA2();
00120 ParsSpotN1 = thePar->ParSpotN1();
00121 ParsSpotN2 = thePar->ParSpotN2();
00122 SamplingResolution = thePar->SamplingResolution();
00123 ConstantResolution = thePar->ConstantResolution();
00124 NoiseResolution = thePar->NoiseResolution();
00125
00126
00127 NSpot = 0.00;
00128 AlphaNSpot = 0.00;
00129 TNSpot = 0.00;
00130 BetaNSpot = 0.00;
00131 RadiusCore = 0.00;
00132 WeightCore = 0.00;
00133 RadiusTail = 0.00;
00134 ComputeZAX0EFFetc();
00135
00136 G4cout << "/********************************************/ " << G4endl;
00137 G4cout << " - GFlashSamplingShowerParameterisation::Constructor - " << G4endl;
00138 G4cout << "/********************************************/ " << G4endl;
00139 }
00140
00141
00142
00143 GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation()
00144 {
00145 if(owning) { delete thePar; }
00146 }
00147
00148
00149
00150 void GFlashSamplingShowerParameterisation::
00151 SetMaterial(G4Material *mat1, G4Material *mat2)
00152 {
00153 G4double Es = 21*MeV;
00154 material1= mat1;
00155 Z1 = GetEffZ(material1);
00156 A1 = GetEffA(material1);
00157 density1 = material1->GetDensity();
00158 X01 = material1->GetRadlen();
00159 Ec1 = 2.66 * std::pow((X01 * Z1 / A1),1.1);
00160 Rm1 = X01*Es/Ec1;
00161
00162 material2= mat2;
00163 Z2 = GetEffZ(material2);
00164 A2 = GetEffA(material2);
00165 density2 = material2->GetDensity();
00166 X02 = material2->GetRadlen();
00167 Ec2 = 2.66 * std::pow((X02 * Z2 / A2),1.1);
00168 Rm2 = X02*Es/Ec2;
00169
00170 }
00171
00172
00173
00174 void GFlashSamplingShowerParameterisation::ComputeZAX0EFFetc()
00175 {
00176 G4cout << "/************ ComputeZAX0EFFetc ************/" << G4endl;
00177 G4cout << " - GFlashSamplingShowerParameterisation::Material - " << G4endl;
00178
00179 G4double Es = 21*MeV;
00180
00181
00182 G4double denominator = (d1*density1 + d2*density2);
00183 G4double W1 = (d1*density1) / denominator;
00184 G4double W2 = (d2*density2)/denominator;
00185 Zeff = ( W1*Z2 ) + (W2*Z1);
00186 Aeff = ( W1*A1 ) + (W2*A2);
00187 X0eff =(1/ (( W1 / X01) +( W2 / X02)));
00188 Rhoeff = ( (d1 *density1 ) + (d2 * density2 ))/G4double (d2 + d1 );
00189 Rmeff = 1/ ((((W1*Ec1)/ X01) + ((W2* Ec2)/ X02) ) / Es ) ;
00190 Eceff = X0eff *((W1*Ec1)/ X01 + (W2* Ec2)/ X02 );
00191 Fs = X0eff/G4double ((d1/mm )+(d2/mm) );
00192 ehat = (1. / (1+ 0.007*(Z1- Z2)));
00193
00194 G4cout << "W1= " << W1 << G4endl;
00195 G4cout << "W2= " << W2 << G4endl;
00196 G4cout << "effective quantities Zeff = "<<Zeff<< G4endl;
00197 G4cout << "effective quantities Aeff = "<<Aeff<< G4endl;
00198 G4cout << "effective quantities Rhoeff = "<<Rhoeff/g *cm3<<" g/cm3" << G4endl;
00199 G4cout << "effective quantities X0eff = "<<X0eff/cm <<" cm" << G4endl;
00200
00201 X0eff = X0eff * Rhoeff;
00202
00203 G4cout << "effective quantities X0eff = "<<X0eff/g*cm2 <<" g/cm2" << G4endl;
00204 X0eff = X0eff /Rhoeff;
00205 G4cout << "effective quantities RMeff = "<<Rmeff/cm<<" cm" << G4endl;
00206 Rmeff = Rmeff* Rhoeff;
00207 G4cout << "effective quantities RMeff = "<<Rmeff/g *cm2<<" g/cm2" << G4endl;
00208 Rmeff = Rmeff/ Rhoeff;
00209 G4cout << "effective quantities Eceff = "<<Eceff/MeV<< " MeV"<< G4endl;
00210 G4cout << "effective quantities Fs = "<<Fs<<G4endl;
00211 G4cout << "effective quantities ehat = "<<ehat<<G4endl;
00212 G4cout << "/********************************************/ " <<G4endl;
00213 }
00214
00215
00216
00217 void GFlashSamplingShowerParameterisation::
00218 GenerateLongitudinalProfile(G4double Energy)
00219 {
00220 if ((material1==0) || (material2 ==0))
00221 {
00222 G4Exception("GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile()",
00223 "InvalidSetup", FatalException, "No material initialized!");
00224 }
00225 G4double y = Energy/Eceff;
00226 ComputeLongitudinalParameters(y);
00227 GenerateEnergyProfile(y);
00228 GenerateNSpotProfile(y);
00229 }
00230
00231
00232
00233 void
00234 GFlashSamplingShowerParameterisation::ComputeLongitudinalParameters(G4double y)
00235 {
00236 AveLogTmaxh = std::log(std::max(ParAveT1 +std::log(y),0.1));
00237 AveLogAlphah = std::log(std::max(ParAveA1 + (ParAveA2+ParAveA3/Zeff)*std::log(y),.1));
00238
00239 SigmaLogTmaxh = std::min(0.5,1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) );
00240 SigmaLogAlphah = std::min(0.5,1.00/( ParSigLogA1 + ParSigLogA2*std::log(y)));
00241 Rhoh = ParRho1+ParRho2*std::log(y);
00242
00243 AveLogTmax = std::max(0.1,std::log(std::exp(AveLogTmaxh)
00244 + ParsAveT1/Fs + ParsAveT2*(1-ehat)));
00245 AveLogAlpha = std::max(0.1,std::log(std::exp(AveLogAlphah)
00246 + (ParsAveA1/Fs)));
00247
00248 SigmaLogTmax = std::min(0.5,1.00/( ParsSigLogT1
00249 + ParsSigLogT2*std::log(y)) );
00250 SigmaLogAlpha = std::min(0.5,1.00/( ParsSigLogA1
00251 + ParsSigLogA2*std::log(y)));
00252 Rho = ParsRho1+ParsRho2*std::log(y);
00253 }
00254
00255
00256
00257 void GFlashSamplingShowerParameterisation::GenerateEnergyProfile(G4double )
00258 {
00259 G4double Correlation1 = std::sqrt((1+Rho)/2);
00260 G4double Correlation2 = std::sqrt((1-Rho)/2);
00261 G4double Correlation1h = std::sqrt((1+Rhoh)/2);
00262 G4double Correlation2h = std::sqrt((1-Rhoh)/2);
00263 G4double Random1 = G4RandGauss::shoot();
00264 G4double Random2 = G4RandGauss::shoot();
00265
00266 Tmax = std::max(1.,std::exp( AveLogTmax + SigmaLogTmax *
00267 (Correlation1*Random1 + Correlation2*Random2) ));
00268 Alpha = std::max(1.1,std::exp( AveLogAlpha + SigmaLogAlpha *
00269 (Correlation1*Random1 - Correlation2*Random2) ));
00270 Beta = (Alpha-1.00)/Tmax;
00271
00272 Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh *
00273 (Correlation1h*Random1 + Correlation2h*Random2) );
00274 Alphah = std::exp( AveLogAlphah + SigmaLogAlphah *
00275 (Correlation1h*Random1 - Correlation2h*Random2) );
00276 Betah = (Alphah-1.00)/Tmaxh;
00277 }
00278
00279
00280
00281 void GFlashSamplingShowerParameterisation::GenerateNSpotProfile(const G4double y)
00282 {
00283 TNSpot = Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff);
00284 TNSpot = std::max(0.5,Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff));
00285 AlphaNSpot = Alphah * (ParsSpotA1+ParsSpotA2*Zeff);
00286 BetaNSpot = (AlphaNSpot-1.00)/TNSpot;
00287 NSpot = ParsSpotN1 /SamplingResolution * std::pow(y*Eceff/GeV,ParsSpotN2 );
00288 }
00289
00290
00291
00292 G4double
00293 GFlashSamplingShowerParameterisation::
00294 ApplySampling(const G4double DEne, const G4double )
00295 {
00296 G4double DEneFluctuated = DEne;
00297 G4double Resolution = std::pow(SamplingResolution,2);
00298
00299
00300
00301
00302
00303
00304 if(Resolution >0.0 && DEne > 0.00)
00305 {
00306 G4float x1=DEne/Resolution;
00307 G4float x2 = CLHEP::RandGamma::shoot(x1, 1.0)*Resolution;
00308 DEneFluctuated=x2;
00309 }
00310 return DEneFluctuated;
00311 }
00312
00313
00314
00315 G4double GFlashSamplingShowerParameterisation::
00316 IntegrateEneLongitudinal(G4double LongitudinalStep)
00317 {
00318 G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
00319 G4float x1= Betah*LongitudinalStepInX0;
00320 G4float x2= Alphah;
00321 float x3 = gam(x1,x2);
00322 G4double DEne=x3;
00323 return DEne;
00324 }
00325
00326
00327
00328 G4double GFlashSamplingShowerParameterisation::
00329 IntegrateNspLongitudinal(G4double LongitudinalStep)
00330 {
00331 G4double LongitudinalStepInX0 = LongitudinalStep / X0eff;
00332 G4float x1 = BetaNSpot*LongitudinalStepInX0;
00333 G4float x2 = AlphaNSpot;
00334 G4float x3 = gam(x1,x2);
00335 G4double DNsp = x3;
00336 return DNsp;
00337 }
00338
00339
00340
00341 G4double GFlashSamplingShowerParameterisation::
00342 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition)
00343 {
00344 if(ispot < 1)
00345 {
00346
00347
00348
00349 G4double Tau = ComputeTau(LongitudinalPosition);
00350 ComputeRadialParameters(Energy,Tau);
00351 }
00352
00353 G4double Radius;
00354 G4double Random1 = G4UniformRand();
00355 G4double Random2 = G4UniformRand();
00356 if(Random1 <WeightCore)
00357 {
00358 Radius = Rmeff * RadiusCore * std::sqrt( Random2/(1. - Random2) );
00359 }
00360 else
00361 {
00362 Radius = Rmeff * RadiusTail * std::sqrt( Random2/(1. - Random2) );
00363 }
00364 Radius = std::min(Radius,DBL_MAX);
00365 return Radius;
00366 }
00367
00368
00369
00370 G4double
00371 GFlashSamplingShowerParameterisation::
00372 ComputeTau(G4double LongitudinalPosition)
00373 {
00374 G4double tau = LongitudinalPosition / Tmax/ X0eff
00375 * (Alpha-1.00) /Alpha
00376 * std::exp(AveLogAlpha)/(std::exp(AveLogAlpha)-1.);
00377 return tau;
00378 }
00379
00380
00381
00382 void GFlashSamplingShowerParameterisation::
00383 ComputeRadialParameters(G4double Energy, G4double Tau)
00384 {
00385 G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV);
00386 G4double z2 = ParRC3+ParRC4*Zeff;
00387 RadiusCore = z1 + z2 * Tau;
00388 G4double p1 = ParWC1+ParWC2*Zeff;
00389 G4double p2 = ParWC3+ParWC4*Zeff;
00390 G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV);
00391 WeightCore = p1 * std::exp( (p2-Tau)/p3- std::exp( (p2-Tau) /p3) );
00392
00393 G4double k1 = ParRT1+ParRT2*Zeff;
00394 G4double k2 = ParRT3;
00395 G4double k3 = ParRT4;
00396 G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV);
00397
00398 RadiusTail = k1*(std::exp(k3*(Tau-k2))
00399 + std::exp(k4*(Tau-k2)) );
00400
00401
00402
00403 RadiusCore = RadiusCore + ParsRC1*(1-ehat) + ParsRC2/Fs*std::exp(-Tau);
00404 WeightCore = WeightCore + (1-ehat)
00405 * (ParsWC1+ParsWC2/Fs * std::exp(-std::pow((Tau-1.),2)));
00406 RadiusTail = RadiusTail + (1-ehat)* ParsRT1+ ParsRT2/Fs *std::exp(-Tau);
00407 }
00408
00409
00410
00411 G4double GFlashSamplingShowerParameterisation::
00412 GenerateExponential(const G4double )
00413 {
00414 G4double ParExp1 = 9./7.*X0eff;
00415 G4double random = -ParExp1*CLHEP::RandExponential::shoot() ;
00416 return random;
00417 }