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