File indexing completed on 2022-08-27 01:28:06
0001
0002
0003
0004
0005
0006 #include "SimG4CMS/Calo/interface/HFCherenkov.h"
0007
0008 #include "G4Poisson.hh"
0009 #include "G4ParticleDefinition.hh"
0010 #include "G4NavigationHistory.hh"
0011 #include "TMath.h"
0012
0013 #include "Randomize.hh"
0014 #include "G4SystemOfUnits.hh"
0015
0016
0017
0018 HFCherenkov::HFCherenkov(edm::ParameterSet const& m_HF) {
0019 ref_index = m_HF.getParameter<double>("RefIndex");
0020 lambda1 = ((m_HF.getParameter<double>("Lambda1")) / 1.e+7) * CLHEP::cm;
0021 lambda2 = ((m_HF.getParameter<double>("Lambda2")) / 1.e+7) * CLHEP::cm;
0022 aperture = m_HF.getParameter<double>("Aperture");
0023 gain = m_HF.getParameter<double>("Gain");
0024 checkSurvive = m_HF.getParameter<bool>("CheckSurvive");
0025 UseNewPMT = m_HF.getParameter<bool>("UseR7600UPMT");
0026 fibreR = m_HF.getParameter<double>("FibreR") * CLHEP::mm;
0027
0028 aperturetrapped = aperture / ref_index;
0029 sinPsimax = 1 / ref_index;
0030
0031 edm::LogVerbatim("HFShower") << "HFCherenkov:: initialised with ref_index " << ref_index << " lambda1/lambda2 (cm) "
0032 << lambda1 / CLHEP::cm << "|" << lambda2 / CLHEP::cm << " aperture(trapped) " << aperture
0033 << "|" << aperturetrapped << " sinPsimax " << sinPsimax
0034 << " Check photon survival in HF " << checkSurvive << " Gain " << gain << " useNewPMT "
0035 << UseNewPMT << " FibreR " << fibreR;
0036
0037 clearVectors();
0038 }
0039
0040 HFCherenkov::~HFCherenkov() {}
0041
0042 int HFCherenkov::computeNPE(const G4Step* aStep,
0043 const G4ParticleDefinition* pDef,
0044 double pBeta,
0045 double u,
0046 double v,
0047 double w,
0048 double step_length,
0049 double zFiber,
0050 double dose,
0051 int npe_Dose) {
0052 clearVectors();
0053 if (!isApplicable(pDef)) {
0054 return 0;
0055 }
0056 if (pBeta < (1 / ref_index) || step_length < 0.0001) {
0057 #ifdef EDM_ML_DEBUG
0058 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: pBeta " << pBeta << " 1/mu " << (1 / ref_index)
0059 << " step_length " << step_length;
0060 #endif
0061 return 0;
0062 }
0063
0064 double uv = std::sqrt(u * u + v * v);
0065 int nbOfPhotons = computeNbOfPhotons(pBeta, step_length) * aStep->GetTrack()->GetWeight();
0066 #ifdef EDM_ML_DEBUG
0067 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: pBeta " << pBeta << " u/v/w " << u << "/" << v << "/" << w
0068 << " step_length " << step_length << " zFib " << zFiber << " nbOfPhotons "
0069 << nbOfPhotons;
0070 #endif
0071 if (nbOfPhotons < 0) {
0072 return 0;
0073 } else if (nbOfPhotons > 0) {
0074 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0075 const G4TouchableHandle& theTouchable = preStepPoint->GetTouchableHandle();
0076 G4ThreeVector localprepos =
0077 theTouchable->GetHistory()->GetTopTransform().TransformPoint(aStep->GetPreStepPoint()->GetPosition());
0078 G4ThreeVector localpostpos =
0079 theTouchable->GetHistory()->GetTopTransform().TransformPoint(aStep->GetPostStepPoint()->GetPosition());
0080
0081 double length = std::sqrt((localpostpos.x() - localprepos.x()) * (localpostpos.x() - localprepos.x()) +
0082 (localpostpos.y() - localprepos.y()) * (localpostpos.y() - localprepos.y()));
0083 double yemit = std::sqrt(fibreR * fibreR - length * length / 4.);
0084
0085 double u_ph = 0, v_ph = 0, w_ph = 0;
0086 for (int i = 0; i < nbOfPhotons; i++) {
0087 double rand = G4UniformRand();
0088 double theta_C = acos(1. / (pBeta * ref_index));
0089 double phi_C = 2 * M_PI * rand;
0090 double sinTheta = sin(theta_C);
0091 double cosTheta = cos(theta_C);
0092 double cosPhi = cos(phi_C);
0093 double sinPhi = sin(phi_C);
0094
0095 if (uv < 0.001) {
0096 u_ph = sinTheta * cosPhi;
0097 v_ph = sinTheta * sinPhi;
0098 w_ph = cosTheta;
0099 } else {
0100 u_ph = uv * cosTheta + sinTheta * cosPhi * w;
0101 v_ph = sinTheta * sinPhi;
0102 w_ph = w * cosTheta - sinTheta * cosPhi * uv;
0103 }
0104 double r_lambda = G4UniformRand();
0105 double lambda0 = (lambda1 * lambda2) / (lambda2 - r_lambda * (lambda2 - lambda1));
0106 double lambda = (lambda0 / CLHEP::cm) * 1.e+7;
0107 wlini.push_back(lambda);
0108 #ifdef EDM_ML_DEBUG
0109 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: " << i << " lambda " << lambda << " w_ph " << w_ph;
0110 #endif
0111
0112 double xemit = length * (G4UniformRand() - 0.5);
0113 double gam = atan2(yemit, xemit);
0114 double eps = atan2(v_ph, u_ph);
0115 double sinBeta = sin(gam - eps);
0116 double rho = std::sqrt(xemit * xemit + yemit * yemit);
0117 double sinEta = rho / fibreR * sinBeta;
0118 double cosEta = std::sqrt(1. - sinEta * sinEta);
0119 double sinPsi = std::sqrt(1. - w_ph * w_ph);
0120 double cosKsi = cosEta * sinPsi;
0121
0122 #ifdef EDM_ML_DEBUG
0123 if (cosKsi < aperturetrapped && w_ph > 0.) {
0124 edm::LogVerbatim("HFShower") << "HFCherenkov::Trapped photon : " << u_ph << " " << v_ph << " " << w_ph << " "
0125 << xemit << " " << gam << " " << eps << " " << sinBeta << " " << rho << " "
0126 << sinEta << " " << cosEta << " "
0127 << " " << sinPsi << " " << cosKsi;
0128 } else {
0129 edm::LogVerbatim("HFShower") << "HFCherenkov::Rejected photon : " << u_ph << " " << v_ph << " " << w_ph << " "
0130 << xemit << " " << gam << " " << eps << " " << sinBeta << " " << rho << " "
0131 << sinEta << " " << cosEta << " "
0132 << " " << sinPsi << " " << cosKsi;
0133 }
0134 #endif
0135 if (cosKsi < aperturetrapped
0136 && w_ph > 0.
0137 && sinPsi < sinPsimax) {
0138 wltrap.push_back(lambda);
0139 double prob_HF = 1.0;
0140 if (checkSurvive) {
0141 double a0_inv = 0.1234;
0142 double a_inv = a0_inv + 0.14 * pow(dose, 0.30);
0143 double z_meters = zFiber;
0144 prob_HF = exp(-z_meters * a_inv);
0145 }
0146 rand = G4UniformRand();
0147 #ifdef EDM_ML_DEBUG
0148 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: probHF " << prob_HF << " Random " << rand
0149 << " Survive? " << (rand < prob_HF);
0150 #endif
0151 if (rand < prob_HF) {
0152 wlatten.push_back(lambda);
0153 double effHEM = computeHEMEff(lambda);
0154
0155 rand = G4UniformRand();
0156 double phi = 0.5 * M_PI * rand;
0157 double rad_bundle = 19.;
0158 double rad_lg = 25.;
0159 rand = G4UniformRand();
0160 double rad = rad_bundle * std::sqrt(rand);
0161 double rho = rad * sin(phi);
0162 double tlength = 2. * std::sqrt(rad_lg * rad_lg - rho * rho);
0163 double length_lg = 430;
0164 double sin_air = sinPsi * 1.46;
0165 double tang = sin_air / std::sqrt(1. - sin_air * sin_air);
0166 int nbounce = length_lg / tlength * tang + 0.5;
0167 double eff = pow(effHEM, nbounce);
0168 rand = G4UniformRand();
0169 #ifdef EDM_ML_DEBUG
0170 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: w_ph " << w_ph << " effHEM " << effHEM << " eff "
0171 << eff << " Random " << rand << " Survive? " << (rand < eff);
0172 #endif
0173 if (rand < eff) {
0174 wlhem.push_back(lambda);
0175 double qEffic = computeQEff(lambda);
0176 rand = G4UniformRand();
0177 #ifdef EDM_ML_DEBUG
0178 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: qEffic " << qEffic << " Random " << rand
0179 << " Survive? " << (rand < qEffic);
0180 #endif
0181 if (rand < qEffic) {
0182 npe_Dose += 1;
0183 momZ.push_back(w_ph);
0184 wl.push_back(lambda);
0185 wlqeff.push_back(lambda);
0186 }
0187 }
0188 }
0189 }
0190 }
0191 }
0192 int npe = npe_Dose;
0193 #ifdef EDM_ML_DEBUG
0194 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: npe " << npe;
0195 #endif
0196 return npe;
0197 }
0198
0199 int HFCherenkov::computeNPEinPMT(
0200 const G4ParticleDefinition* pDef, double pBeta, double u, double v, double w, double step_length) {
0201 clearVectors();
0202 int npe_ = 0;
0203 if (!isApplicable(pDef)) {
0204 return 0;
0205 }
0206 if (pBeta < (1 / ref_index) || step_length < 0.0001) {
0207 #ifdef EDM_ML_DEBUG
0208 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: pBeta " << pBeta << " 1/mu " << (1 / ref_index)
0209 << " step_length " << step_length;
0210 #endif
0211 return 0;
0212 }
0213
0214 double uv = std::sqrt(u * u + v * v);
0215 int nbOfPhotons = computeNbOfPhotons(pBeta, step_length);
0216 #ifdef EDM_ML_DEBUG
0217 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: pBeta " << pBeta << " u/v/w " << u << "/" << v << "/"
0218 << w << " step_length " << step_length << " nbOfPhotons " << nbOfPhotons;
0219 #endif
0220 if (nbOfPhotons < 0) {
0221 return 0;
0222 } else if (nbOfPhotons > 0) {
0223 double w_ph = 0;
0224 for (int i = 0; i < nbOfPhotons; i++) {
0225 double rand = G4UniformRand();
0226 double theta_C = acos(1. / (pBeta * ref_index));
0227 double phi_C = 2 * M_PI * rand;
0228 double sinTheta = sin(theta_C);
0229 double cosTheta = cos(theta_C);
0230 double cosPhi = cos(phi_C);
0231
0232 if (uv < 0.001) {
0233 w_ph = cosTheta;
0234 } else {
0235 w_ph = w * cosTheta - sinTheta * cosPhi * uv;
0236 }
0237 double r_lambda = G4UniformRand();
0238 double lambda0 = (lambda1 * lambda2) / (lambda2 - r_lambda * (lambda2 - lambda1));
0239 double lambda = (lambda0 / CLHEP::cm) * 1.e+7;
0240 wlini.push_back(lambda);
0241 #ifdef EDM_ML_DEBUG
0242 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: " << i << " lambda " << lambda << " w_ph " << w_ph
0243 << " aperture " << aperture;
0244 #endif
0245 if (w_ph > aperture) {
0246 wltrap.push_back(lambda);
0247 rand = G4UniformRand();
0248 #ifdef EDM_ML_DEBUG
0249 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: Random " << rand << " Survive? " << (rand < 1.);
0250 #endif
0251 if (rand < 1.0) {
0252 wlatten.push_back(lambda);
0253 double qEffic = computeQEff(lambda);
0254 rand = G4UniformRand();
0255 #ifdef EDM_ML_DEBUG
0256 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: qEffic " << qEffic << " Random " << rand
0257 << " Survive? " << (rand < qEffic);
0258 #endif
0259 if (rand < qEffic) {
0260 npe_ += 1;
0261 momZ.push_back(w_ph);
0262 wl.push_back(lambda);
0263 wlqeff.push_back(lambda);
0264 }
0265 }
0266 }
0267 }
0268 }
0269 #ifdef EDM_ML_DEBUG
0270 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: npe " << npe_;
0271 #endif
0272 return npe_;
0273 }
0274
0275 std::vector<double> HFCherenkov::getWLIni() {
0276 std::vector<double> v = wlini;
0277 return v;
0278 }
0279
0280 std::vector<double> HFCherenkov::getWLTrap() {
0281 std::vector<double> v = wltrap;
0282 return v;
0283 }
0284
0285 std::vector<double> HFCherenkov::getWLAtten() {
0286 std::vector<double> v = wlatten;
0287 return v;
0288 }
0289
0290 std::vector<double> HFCherenkov::getWLHEM() {
0291 std::vector<double> v = wlhem;
0292 return v;
0293 }
0294
0295 std::vector<double> HFCherenkov::getWLQEff() {
0296 std::vector<double> v = wlqeff;
0297 return v;
0298 }
0299
0300 std::vector<double> HFCherenkov::getWL() {
0301 std::vector<double> v = wl;
0302 return v;
0303 }
0304
0305 std::vector<double> HFCherenkov::getMom() {
0306 std::vector<double> v = momZ;
0307 return v;
0308 }
0309
0310 int HFCherenkov::computeNbOfPhotons(double beta, G4double stepL) {
0311 double pBeta = beta;
0312 double alpha = 0.0073;
0313 double step_length = stepL;
0314 double theta_C = acos(1. / (pBeta * ref_index));
0315 double lambdaDiff = (1. / lambda1 - 1. / lambda2);
0316 double cherenPhPerLength = 2 * M_PI * alpha * lambdaDiff * CLHEP::cm;
0317 double d_NOfPhotons = cherenPhPerLength * sin(theta_C) * sin(theta_C) * (step_length / CLHEP::cm);
0318 int nbOfPhotons = int(d_NOfPhotons);
0319 #ifdef EDM_ML_DEBUG
0320 edm::LogVerbatim("HFShower") << "HFCherenkov::computeNbOfPhotons: StepLength " << step_length << " theta_C "
0321 << theta_C << " lambdaDiff " << lambdaDiff << " cherenPhPerLength " << cherenPhPerLength
0322 << " Photons " << d_NOfPhotons << " " << nbOfPhotons;
0323 #endif
0324 return nbOfPhotons;
0325 }
0326
0327 double HFCherenkov::computeQEff(double wavelength) {
0328 double qeff(0.);
0329 if (UseNewPMT) {
0330 if (wavelength <= 350) {
0331 qeff = 2.45867 * (TMath::Landau(wavelength, 353.820, 59.1324));
0332 } else if (wavelength > 350 && wavelength < 500) {
0333 qeff = 0.441989 * exp(-pow((wavelength - 358.371), 2) / (2 * pow((138.277), 2)));
0334 } else if (wavelength >= 500 && wavelength < 550) {
0335 qeff = 0.271862 * exp(-pow((wavelength - 491.505), 2) / (2 * pow((47.0418), 2)));
0336 } else if (wavelength >= 550) {
0337 qeff = 0.137297 * exp(-pow((wavelength - 520.260), 2) / (2 * pow((75.5023), 2)));
0338 }
0339 #ifdef EDM_ML_DEBUG
0340 edm::LogVerbatim("HFShower") << "HFCherenkov:: for new PMT : wavelength === " << wavelength << "\tqeff ===\t"
0341 << qeff;
0342 #endif
0343 } else {
0344 double y = (wavelength - 275.) / 180.;
0345 double func = 1. / (1. + 250. * pow((y / 5.), 4));
0346 double qE_R7525 = 0.77 * y * exp(-y) * func;
0347 qeff = qE_R7525;
0348 #ifdef EDM_ML_DEBUG
0349 edm::LogVerbatim("HFShower") << "HFCherenkov::computeQEff: wavelength " << wavelength << " y/func " << y << "/"
0350 << func << " qeff " << qeff;
0351 edm::LogVerbatim("HFShower") << "HFCherenkov:: for old PMT : wavelength === " << wavelength << "; qeff = " << qeff;
0352 #endif
0353 }
0354 return qeff;
0355 }
0356
0357 double HFCherenkov::computeHEMEff(double wavelength) {
0358 double hEMEff = 0;
0359 if (wavelength < 400.) {
0360 hEMEff = 0.;
0361 } else if (wavelength >= 400. && wavelength < 410.) {
0362
0363 hEMEff = (-1322.453 / wavelength) + 4.2056;
0364 } else if (wavelength >= 410.) {
0365 hEMEff = 0.99;
0366 if (wavelength > 415. && wavelength < 445.) {
0367
0368
0369 hEMEff = 0.97;
0370 } else if (wavelength > 550. && wavelength < 600.) {
0371
0372
0373 hEMEff = 0.98;
0374 } else if (wavelength > 565. && wavelength <= 635.) {
0375
0376 hEMEff = (701.7268 / wavelength) - 0.186;
0377 }
0378 }
0379 #ifdef EDM_ML_DEBUG
0380 edm::LogVerbatim("HFShower") << "HFCherenkov::computeHEMEff: wavelength " << wavelength << " hEMEff " << hEMEff;
0381 #endif
0382 return hEMEff;
0383 }
0384
0385 double HFCherenkov::smearNPE(int npe) {
0386 double pe = 0.;
0387 if (npe > 0) {
0388 for (int i = 0; i < npe; ++i) {
0389 double val = G4Poisson(gain);
0390 pe += (val / gain) + 0.001 * G4UniformRand();
0391 }
0392 }
0393 #ifdef EDM_ML_DEBUG
0394 edm::LogVerbatim("HFShower") << "HFCherenkov::smearNPE: npe " << npe << " pe " << pe;
0395 #endif
0396 return pe;
0397 }
0398
0399 void HFCherenkov::clearVectors() {
0400 wl.clear();
0401 wlini.clear();
0402 wltrap.clear();
0403 wlatten.clear();
0404 wlhem.clear();
0405 wlqeff.clear();
0406 momZ.clear();
0407 }
0408
0409 bool HFCherenkov::isApplicable(const G4ParticleDefinition* aParticleType) {
0410 bool tmp = (aParticleType->GetPDGCharge() != 0);
0411 #ifdef EDM_ML_DEBUG
0412 edm::LogVerbatim("HFShower") << "HFCherenkov::isApplicable: aParticleType " << aParticleType->GetParticleName()
0413 << " PDGCharge " << aParticleType->GetPDGCharge() << " Result " << tmp;
0414 #endif
0415 return tmp;
0416 }