include/G4HepEm/G4HepEmElectronEnergyLossFluctuation.icc

0001
0002 #include "G4HepEmElectronEnergyLossFluctuation.hh"
0003
0004
0005 #include "G4HepEmRandomEngine.hh"
0006
0007 #include "G4HepEmConstants.hh"
0008 #include "G4HepEmMath.hh"
0009
0010
0011 double G4HepEmElectronEnergyLossFluctuation::SampleEnergyLossFLuctuation(double /*ekin*/, double tcut, double /*tmax*/,
0012        double excEner, double /*stepLength*/, double meanELoss, G4HepEmRandomEngine* rnge) {
0013   const double scaling  = G4HepEmMin(1. + 5.E-4/tcut, 1.5);
0014   const double meanLoss = meanELoss/scaling;
0015
0016   const double kFluctParRate     = 0.56;
0017   const double kFluctParE0       = 1.E-5; // 10 eV
0018   const double kFluctParNMaxCont = 8.;
0019
0020   const double w1 = tcut/kFluctParE0;
0021   double a3 = meanLoss*(tcut - kFluctParE0)/(kFluctParE0*tcut*G4HepEmLog(w1));
0022   double a1 = 0.;
0023   double e1 = excEner;
0024   double eloss = 0.0;
0025   // 1. excittaion part
0026   if (tcut > excEner) {
0027     const double a1Tmp = meanLoss*(1. - kFluctParRate)/excEner;
0028     // NOTE: this corresponds to G4UniversalFluctuation as in G4-v11.p01
0029     const double kFluctParA0 = 42.;
0030     const double kFluctParFw =  4.;
0031     const double dum0  = a1Tmp < kFluctParA0
0032                          ? .1 + (kFluctParFw - .1)*std::sqrt(a1Tmp/kFluctParA0)
0033                          : kFluctParFw;
0034     // NOTE: this corresponds to G4UniversalFluctuation as in G4-v11.00
0035     // const double kFluctParA0 = 15.;
0036     // const double kFluctParFw =  5.;
0037     // const double dum0  = a1Tmp < kFluctParA0
0038     //                      ? kFluctParFw*a1Tmp/kFluctParA0
0039     //                      : kFluctParFw;
0040     a1  = a1Tmp/dum0;
0041     e1 *= dum0;
0042     a3 *= kFluctParRate;
0043     //
0044     // add excition (a1 > 0)
0045     if (a1 > kFluctParNMaxCont) {
0046       // Gaussian
0047       const double emean = a1*e1;
0048       const double sig2e = emean*e1;
0049       eloss = SampleGaussianLoss(emean, sig2e, rnge);
0050     } else {
0051       // small number --> sampling from Poisson
0052       const int p = rnge->Poisson(a1);
0053       eloss = p > 0 ? ((p + 1) - 2.*rnge->flat())*e1 : 0.;
0054     }
0055   }
0056   //
0057   // 2. ionisation part
0058   if (a3 > 0.) {
0059     double   p3 = a3;
0060     double alfa = 1.;
0061     if (a3 > kFluctParNMaxCont) {
0062       alfa = w1*(kFluctParNMaxCont + a3)/(w1*kFluctParNMaxCont + a3);
0063       const double alfa1  = alfa*G4HepEmLog(alfa)/(alfa - 1.);
0064       const double namean = a3*w1*(alfa - 1.)/((w1 - 1.)*alfa);
0065       const double emean  = namean*kFluctParE0*alfa1;
0066       const double sig2e  = kFluctParE0*kFluctParE0*namean*(alfa - alfa1*alfa1);
0067       eloss += SampleGaussianLoss(emean, sig2e, rnge);
0068       p3 = a3 - namean;
0069     }
0070     //
0071     const double w3 = alfa*kFluctParE0;
0072     if (tcut > w3) {
0073       const double w = (tcut - w3)/tcut;
0074       const int  nnb = rnge->Poisson(p3);
0075       if (nnb > 0) {
0076         const int kBlockSize = 8;
0077         const int nBlocks    = nnb/kBlockSize;
0078         //
0079         double rndm[kBlockSize];
0080         for (int ib=0; ib<nBlocks; ++ib) {
0081           rnge->flatArray(kBlockSize, rndm);
0082           for (int i=0; i<kBlockSize; ++i) {
0083             eloss += w3/(1.-w*rndm[i]);
0084           }
0085         }
0086         const int nTail = nnb - nBlocks*kBlockSize;
0087         rnge->flatArray(nTail, rndm);
0088         for (int i=0; i<nTail; ++i) {
0089           eloss += w3/(1.-w*rndm[i]);
0090         }
0091       }
0092     }
0093   }
0094   //
0095   // deliver result
0096   return eloss*scaling;
0097 }
0098
0099
0100 double G4HepEmElectronEnergyLossFluctuation::SampleGaussianLoss(double meane, double sig2e, G4HepEmRandomEngine* rnge) {
0101   const double twom = 2.*meane;
0102   if (meane*meane < 0.0625*sig2e) {
0103     return twom*rnge->flat();
0104   }
0105   const double sig = std::sqrt(sig2e);
0106   double eloss;
0107   do {
0108     eloss = rnge->Gauss(meane, sig);
0109   } while (eloss < 0. || eloss > twom);
0110   return eloss;
0111 }