Common/src/SiG4UniversalFluctuation.cc

0001
0002 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0003 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0004
0005 #include <CLHEP/Random/RandFlat.h>
0006 #include <CLHEP/Random/RandGaussQ.h>
0007 #include <CLHEP/Random/RandPoissonQ.h>
0008 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0009 #include <CLHEP/Units/SystemOfUnits.h>
0010 #include "SimTracker/Common/interface/SiG4UniversalFluctuation.h"
0011 #include "vdt/log.h"
0012 #include <cmath>
0013
0014 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0015
0016 using namespace std;
0017
0018 SiG4UniversalFluctuation::SiG4UniversalFluctuation()
0019     : minNumberInteractionsBohr(10.0),
0020       theBohrBeta2(50.0 * CLHEP::keV / proton_mass_c2),
0021       minLoss(10. * CLHEP::eV),
0022       problim(5.e-3),
0023       alim(10.),
0024       nmaxCont1(4.),
0025       nmaxCont2(16.) {
0026   sumalim = -log(problim);
0027
0028   // Add these definitions d.k.
0029   chargeSquare = 1.;  // Assume all particles have charge 1
0030   // Taken from Geant4 printout, HARDWIRED for Silicon.
0031   ipotFluct = 0.0001736;        // material->GetIonisation()->GetMeanExcitationEnergy();
0032   electronDensity = 6.797E+20;  // material->GetElectronDensity();
0033   f1Fluct = 0.8571;             // material->GetIonisation()->GetF1fluct();
0034   f2Fluct = 0.1429;             // material->GetIonisation()->GetF2fluct();
0035   e1Fluct = 0.000116;           // material->GetIonisation()->GetEnergy1fluct();
0036   e2Fluct = 0.00196;            // material->GetIonisation()->GetEnergy2fluct();
0037   e1LogFluct = -9.063;          // material->GetIonisation()->GetLogEnergy1fluct();
0038   e2LogFluct = -6.235;          // material->GetIonisation()->GetLogEnergy2fluct();
0039   rateFluct = 0.4;              // material->GetIonisation()->GetRateionexcfluct();
0040   ipotLogFluct = -8.659;        // material->GetIonisation()->GetLogMeanExcEnergy();
0041   e0 = 1.E-5;                   // material->GetIonisation()->GetEnergy0fluct();
0042
0043   // cout << " init new fluct +++++++++++++++++++++++++++++++++++++++++"<<endl;
0044 }
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047 // The main dedx fluctuation routine.
0048 // Arguments: momentum in MeV/c, mass in MeV, delta ray cut (tmax) in
0049 // MeV, silicon thickness in mm, mean eloss in MeV.
0050
0051 SiG4UniversalFluctuation::~SiG4UniversalFluctuation() {}
0052
0053 double SiG4UniversalFluctuation::SampleFluctuations(const double momentum,
0054                                                     const double mass,
0055                                                     double &tmax,
0056                                                     const double length,
0057                                                     const double meanLoss,
0058                                                     CLHEP::HepRandomEngine *engine) {
0059   // Calculate actual loss from the mean loss.
0060   // The model used to get the fluctuations is essentially the same
0061   // as in Glandz in Geant3 (Cern program library W5013, phys332).
0062   // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
0063
0064   // shortcut for very very small loss (out of validity of the model)
0065   //
0066   if (meanLoss < minLoss)
0067     return meanLoss;
0068
0069   particleMass = mass;
0070   double gam2 = (momentum * momentum) / (particleMass * particleMass) + 1.0;
0071   double beta2 = 1.0 - 1.0 / gam2;
0072   double gam = sqrt(gam2);
0073
0074   double loss(0.), siga(0.);
0075
0076   // Gaussian regime
0077   // for heavy particles only and conditions
0078   // for Gauusian fluct. has been changed
0079   //
0080   if ((particleMass > electron_mass_c2) && (meanLoss >= minNumberInteractionsBohr * tmax)) {
0081     double massrate = electron_mass_c2 / particleMass;
0082     double tmaxkine = 2. * electron_mass_c2 * beta2 * gam2 / (1. + massrate * (2. * gam + massrate));
0083     if (tmaxkine <= 2. * tmax) {
0084       siga = (1.0 / beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length * electronDensity * chargeSquare;
0085       siga = sqrt(siga);
0086       double twomeanLoss = meanLoss + meanLoss;
0087       if (twomeanLoss < siga) {
0088         double x;
0089         do {
0090           loss = twomeanLoss * CLHEP::RandFlat::shoot(engine);
0091           x = (loss - meanLoss) / siga;
0092         } while (1.0 - 0.5 * x * x < CLHEP::RandFlat::shoot(engine));
0093       } else {
0094         do {
0095           loss = CLHEP::RandGaussQ::shoot(engine, meanLoss, siga);
0096         } while (loss < 0. || loss > twomeanLoss);
0097       }
0098       return loss;
0099     }
0100   }
0101
0102   double a1 = 0., a2 = 0., a3 = 0.;
0103   double p3;
0104   double rate = rateFluct;
0105
0106   double w1 = tmax / ipotFluct;
0107   double w2 = vdt::fast_log(2. * electron_mass_c2 * beta2 * gam2) - beta2;
0108
0109   if (w2 > ipotLogFluct) {
0110     double C = meanLoss * (1. - rateFluct) / (w2 - ipotLogFluct);
0111     a1 = C * f1Fluct * (w2 - e1LogFluct) / e1Fluct;
0112     a2 = C * f2Fluct * (w2 - e2LogFluct) / e2Fluct;
0113     if (a2 < 0.) {
0114       a1 = 0.;
0115       a2 = 0.;
0116       rate = 1.;
0117     }
0118   } else {
0119     rate = 1.;
0120   }
0121
0122   // added
0123   if (tmax > ipotFluct) {
0124     a3 = rate * meanLoss * (tmax - ipotFluct) / (ipotFluct * tmax * vdt::fast_log(w1));
0125   }
0126   double suma = a1 + a2 + a3;
0127
0128   // Glandz regime
0129   //
0130   if (suma > sumalim) {
0131     if ((a1 + a2) > 0.) {
0132       double p1, p2;
0133       // excitation type 1
0134       if (a1 > alim) {
0135         siga = sqrt(a1);
0136         p1 = max(0., CLHEP::RandGaussQ::shoot(engine, a1, siga) + 0.5);
0137       } else {
0138         p1 = double(CLHEP::RandPoissonQ::shoot(engine, a1));
0139       }
0140
0141       // excitation type 2
0142       if (a2 > alim) {
0143         siga = sqrt(a2);
0144         p2 = max(0., CLHEP::RandGaussQ::shoot(engine, a2, siga) + 0.5);
0145       } else {
0146         p2 = double(CLHEP::RandPoissonQ::shoot(engine, a2));
0147       }
0148
0149       loss = p1 * e1Fluct + p2 * e2Fluct;
0150
0151       // smearing to avoid unphysical peaks
0152       if (p2 > 0.)
0153         loss += (1. - 2. * CLHEP::RandFlat::shoot(engine)) * e2Fluct;
0154       else if (loss > 0.)
0155         loss += (1. - 2. * CLHEP::RandFlat::shoot(engine)) * e1Fluct;
0156       if (loss < 0.)
0157         loss = 0.0;
0158     }
0159
0160     // ionisation
0161     if (a3 > 0.) {
0162       if (a3 > alim) {
0163         siga = sqrt(a3);
0164         p3 = max(0., CLHEP::RandGaussQ::shoot(engine, a3, siga) + 0.5);
0165       } else {
0166         p3 = double(CLHEP::RandPoissonQ::shoot(engine, a3));
0167       }
0168       double lossc = 0.;
0169       if (p3 > 0) {
0170         double na = 0.;
0171         double alfa = 1.;
0172         if (p3 > nmaxCont2) {
0173           double rfac = p3 / (nmaxCont2 + p3);
0174           double namean = p3 * rfac;
0175           double sa = nmaxCont1 * rfac;
0176           na = CLHEP::RandGaussQ::shoot(engine, namean, sa);
0177           if (na > 0.) {
0178             alfa = w1 * (nmaxCont2 + p3) / (w1 * nmaxCont2 + p3);
0179             double alfa1 = alfa * vdt::fast_log(alfa) / (alfa - 1.);
0180             double ea = na * ipotFluct * alfa1;
0181             double sea = ipotFluct * sqrt(na * (alfa - alfa1 * alfa1));
0182             lossc += CLHEP::RandGaussQ::shoot(engine, ea, sea);
0183           }
0184         }
0185
0186         if (p3 > na) {
0187           w2 = alfa * ipotFluct;
0188           double w = (tmax - w2) / tmax;
0189           int nb = int(p3 - na);
0190           for (int k = 0; k < nb; k++)
0191             lossc += w2 / (1. - w * CLHEP::RandFlat::shoot(engine));
0192         }
0193       }
0194       loss += lossc;
0195     }
0196     return loss;
0197   }
0198
0199   // suma < sumalim;  very small energy loss;
0200   a3 = meanLoss * (tmax - e0) / (tmax * e0 * vdt::fast_log(tmax / e0));
0201   if (a3 > alim) {
0202     siga = sqrt(a3);
0203     p3 = max(0., CLHEP::RandGaussQ::shoot(engine, a3, siga) + 0.5);
0204   } else {
0205     p3 = double(CLHEP::RandPoissonQ::shoot(engine, a3));
0206   }
0207   if (p3 > 0.) {
0208     double w = (tmax - e0) / tmax;
0209     double corrfac = 1.;
0210     if (p3 > nmaxCont2) {
0211       corrfac = p3 / nmaxCont2;
0212       p3 = nmaxCont2;
0213     }
0214     int ip3 = (int)p3;
0215     for (int i = 0; i < ip3; i++)
0216       loss += 1. / (1. - w * CLHEP::RandFlat::shoot(engine));
0217     loss *= e0 * corrfac;
0218     // smearing for losses near to e0
0219     if (p3 <= 2.)
0220       loss += e0 * (1. - 2. * CLHEP::RandFlat::shoot(engine));
0221   }
0222
0223   return loss;
0224 }