Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-10 02:59:06 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File:  HFCherenkov.cc
0003 // Description: Generate Cherenkov photons for HF
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 //#define EDM_ML_DEBUG
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       //photon momentum
0096       if (uv < 0.001) {  // aligned with z-axis
0097         u_ph = sinTheta * cosPhi;
0098         v_ph = sinTheta * sinPhi;
0099         w_ph = cosTheta;
0100       } else {  // general case
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;  // lambda is in nm
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  // photon is trapped inside fiber
0137           && w_ph > 0.              // and moves to PMT
0138           && sinPsi < sinPsimax) {  // and is not reflected at fiber end
0139         wltrap.push_back(lambda);
0140         double prob_HF = 1.0;
0141         if (checkSurvive) {
0142           double a0_inv = 0.1234;                          //meter^-1
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) {  // survived in HF
0153           wlatten.push_back(lambda);
0154           double effHEM = computeHEMEff(lambda);
0155           // compute number of bounces in air guide
0156           rand = G4UniformRand();
0157           double phi = 0.5 * M_PI * rand;
0158           double rad_bundle = 19.;  // bundle radius
0159           double rad_lg = 25.;      //  light guide radius
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) {  // survived HEM
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) {  // made photoelectron
0183               npe_Dose += 1;
0184               momZ.push_back(w_ph);
0185               wl.push_back(lambda);
0186               wlqeff.push_back(lambda);
0187             }  // made pe
0188           }  // passed HEM
0189         }  // passed fiber
0190       }  // end of  if(w_ph < w_aperture), trapped inside fiber
0191     }  // end of ++NbOfPhotons
0192   }  // end of if(NbOfPhotons)}
0193   int npe = npe_Dose;  // Nb of photoelectrons
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       //photon momentum
0233       if (uv < 0.001) {  // aligned with z-axis
0234         w_ph = cosTheta;
0235       } else {  // general case
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;  // lambda is in nm
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) {  // phton trapped inside PMT glass
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) {  // survived all the times and sent to photo-cathode
0253           wlatten.push_back(lambda);
0254           double qEffic = computeQEff(lambda);  //Quantum efficiency of the PMT
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) {  // made photoelectron
0261             npe_ += 1;
0262             momZ.push_back(w_ph);
0263             wl.push_back(lambda);
0264             wlqeff.push_back(lambda);
0265           }  // made pe
0266         }  // accepted all Cherenkov photons
0267       }  // end of  if(w_ph < w_aperture), trapped inside glass
0268     }  // end of ++NbOfPhotons
0269   }  // end of if(NbOfPhotons)}
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     //hEMEff = .99 * (wavelength - 400.) / 10.;
0364     hEMEff = (-1322.453 / wavelength) + 4.2056;
0365   } else if (wavelength >= 410.) {
0366     hEMEff = 0.99;
0367     if (wavelength > 415. && wavelength < 445.) {
0368       //abs(wavelength - 430.) < 15.
0369       //hEMEff = 0.95;
0370       hEMEff = 0.97;
0371     } else if (wavelength > 550. && wavelength < 600.) {
0372       // abs(wavelength - 575.) < 25.)
0373       //hEMEff = 0.96;
0374       hEMEff = 0.98;
0375     } else if (wavelength > 565. && wavelength <= 635.) {  // added later
0376       // abs(wavelength - 600.) < 35.)
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!