File indexing completed on 2024-05-10 02:21:29
0001
0002
0003
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
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
0029 chargeSquare = 1.;
0030
0031 ipotFluct = 0.0001736;
0032 electronDensity = 6.797E+20;
0033 f1Fluct = 0.8571;
0034 f2Fluct = 0.1429;
0035 e1Fluct = 0.000116;
0036 e2Fluct = 0.00196;
0037 e1LogFluct = -9.063;
0038 e2LogFluct = -6.235;
0039 rateFluct = 0.4;
0040 ipotLogFluct = -8.659;
0041 e0 = 1.E-5;
0042
0043
0044 }
0045
0046
0047
0048
0049
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
0060
0061
0062
0063
0064
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
0077
0078
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
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
0129
0130 if (suma > sumalim) {
0131 if ((a1 + a2) > 0.) {
0132 double p1, p2;
0133
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
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
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
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
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
0219 if (p3 <= 2.)
0220 loss += e0 * (1. - 2. * CLHEP::RandFlat::shoot(engine));
0221 }
0222
0223 return loss;
0224 }