Back to home page

Project CMSSW displayed by LXR
 
 

    

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