#include <HFCherenkov.h>

Public Member Functions
void	clearVectors ()

int	computeNPE (const G4Step step, const G4ParticleDefinition pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)

int	computeNPEinPMT (const G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length)

int	computeNPhTrapped (double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)

std::vector< double >	getMom ()

std::vector< double >	getWL ()

std::vector< double >	getWLAtten ()

std::vector< double >	getWLHEM ()

std::vector< double >	getWLIni ()

std::vector< double >	getWLQEff ()

std::vector< double >	getWLTrap ()

	HFCherenkov (edm::ParameterSet const &p)

double	smearNPE (G4int Npe)

virtual	~HFCherenkov ()

Private Member Functions
double	computeHEMEff (double wavelength)

int	computeNbOfPhotons (double pBeta, double step_length)

double	computeQEff (double wavelength)

bool	isApplicable (const G4ParticleDefinition *aParticleType)

Private Attributes
double	aperture

double	apertureTrap

double	aperturetrapped

bool	checkSurvive

double	fibreR

double	gain

double	lambda1

double	lambda2

std::vector< double >	momZ

G4ThreeVector	phMom

double	ref_index

double	sinPsimax

bool	UseNewPMT

std::vector< double >	wl

std::vector< double >	wlatten

std::vector< double >	wlhem

std::vector< double >	wlini

std::vector< double >	wlqeff

std::vector< double >	wltrap

Detailed Description

Definition at line 19 of file HFCherenkov.h.

Constructor & Destructor Documentation

HFCherenkov::HFCherenkov ( edm::ParameterSet const & p )

Definition at line 18 of file HFCherenkov.cc.

References aperture, apertureTrap, aperturetrapped, checkSurvive, clearVectors(), funct::cos(), fibreR, gain, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), lambda1, lambda2, funct::pow(), ref_index, sinPsimax, and UseNewPMT.

                                                     {
   ref_index = m_HF.getParameter<double>("RefIndex");
   lambda1 = ((m_HF.getParameter<double>("Lambda1")) / pow(double(10), 7)) * cm;
   lambda2 = ((m_HF.getParameter<double>("Lambda2")) / pow(double(10), 7)) * cm;
   aperture = cos(asin(m_HF.getParameter<double>("Aperture")));
   apertureTrap = cos(asin(m_HF.getParameter<double>("ApertureTrapped")));
   aperturetrapped = m_HF.getParameter<double>("CosApertureTrapped");
   gain = m_HF.getParameter<double>("Gain");
   checkSurvive = m_HF.getParameter<bool>("CheckSurvive");
   UseNewPMT = m_HF.getParameter<bool>("UseR7600UPMT");
   sinPsimax = m_HF.getUntrackedParameter<double>("SinPsiMax", 0.5);
   fibreR = m_HF.getUntrackedParameter<double>("FibreR", 0.3) * mm;
 
   edm::LogVerbatim("HFShower") << "HFCherenkov:: initialised with ref_index " << ref_index << " lambda1/lambda2 (cm) "
                                << lambda1 / cm << "|" << lambda2 / cm << " aperture(total/trapped) " << aperture << "|"
                                << apertureTrap << "|" << aperturetrapped << " Check photon survival in HF "
                                << checkSurvive << " Gain " << gain << " useNewPMT " << UseNewPMT << " FibreR "
                                << fibreR;
 
   clearVectors();
 }

HFCherenkov::~HFCherenkov ( )

virtual

Definition at line 40 of file HFCherenkov.cc.

40 {}

Member Function Documentation

void HFCherenkov::clearVectors ( )

Definition at line 421 of file HFCherenkov.cc.

References momZ, wl, wlatten, wlhem, wlini, wlqeff, and wltrap.

Referenced by computeNPE(), computeNPEinPMT(), and HFCherenkov().

                                {
   wl.clear();
   wlini.clear();
   wltrap.clear();
   wlatten.clear();
   wlhem.clear();
   wlqeff.clear();
   momZ.clear();
 }

double HFCherenkov::computeHEMEff ( double wavelength )

private

Definition at line 379 of file HFCherenkov.cc.

Referenced by computeNPE().

                                                    {
   double hEMEff = 0;
   if (wavelength < 400.) {
     hEMEff = 0.;
   } else if (wavelength >= 400. && wavelength < 410.) {
     //hEMEff = .99 * (wavelength - 400.) / 10.;
     hEMEff = (-1322.453 / wavelength) + 4.2056;
   } else if (wavelength >= 410.) {
     hEMEff = 0.99;
     if (wavelength > 415. && wavelength < 445.) {
       //abs(wavelength - 430.) < 15.
       //hEMEff = 0.95;
       hEMEff = 0.97;
     } else if (wavelength > 550. && wavelength < 600.) {
       // abs(wavelength - 575.) < 25.)
       //hEMEff = 0.96;
       hEMEff = 0.98;
     } else if (wavelength > 565. && wavelength <= 635.) {  // added later
       // abs(wavelength - 600.) < 35.)
       hEMEff = (701.7268 / wavelength) - 0.186;
     }
   }
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::computeHEMEff: wavelength " << wavelength << " hEMEff " << hEMEff;
 #endif
   return hEMEff;
 }

int HFCherenkov::computeNbOfPhotons	(	double	pBeta,
		double	step_length
	)

private

Definition at line 332 of file HFCherenkov.cc.

References alpha, pfBoostedDoubleSVAK8TagInfos_cfi::beta, createfilelist::int, lambda1, lambda2, M_PI, ref_index, and funct::sin().

Referenced by computeNPE(), computeNPEinPMT(), and computeNPhTrapped().

                                                                {
   double pBeta = beta;
   double alpha = 0.0073;
   double step_length = stepL;
   double theta_C = acos(1. / (pBeta * ref_index));
   double lambdaDiff = (1. / lambda1 - 1. / lambda2);
   double cherenPhPerLength = 2 * M_PI * alpha * lambdaDiff * cm;
   double d_NOfPhotons = cherenPhPerLength * sin(theta_C) * sin(theta_C) * (step_length / cm);
   int nbOfPhotons = int(d_NOfPhotons);
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::computeNbOfPhotons: StepLength " << step_length << " theta_C "
                                << theta_C << " lambdaDiff " << lambdaDiff << " cherenPhPerLength " << cherenPhPerLength
                                << " Photons " << d_NOfPhotons << " " << nbOfPhotons;
 #endif
   return nbOfPhotons;
 }

int HFCherenkov::computeNPE	(	const G4Step *	step,
		const G4ParticleDefinition *	pDef,
		double	pBeta,
		double	u,
		double	v,
		double	w,
		double	step_length,
		double	zFiber,
		double	Dose,
		int	Npe_Dose
	)

Definition at line 77 of file HFCherenkov.cc.

References aperture, aperturetrapped, checkSurvive, clearVectors(), computeHEMEff(), computeNbOfPhotons(), computeQEff(), funct::cos(), JetChargeProducer_cfi::exp, fibreR, mps_fire::i, isApplicable(), lambda1, lambda2, M_PI, momZ, funct::pow(), rand(), ref_index, rho, funct::sin(), sinPsimax, mathSSE::sqrt(), w, wl, wlatten, wlhem, wlini, wlqeff, and wltrap.

Referenced by HFShower::getHits().

                                           {
   clearVectors();
   if (!isApplicable(pDef)) {
     return 0;
   }
   if (pBeta < (1 / ref_index) || step_length < 0.0001) {
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: pBeta " << pBeta << " 1/mu " << (1 / ref_index)
                                  << " step_length " << step_length;
 #endif
     return 0;
   }
 
   double uv = sqrt(u * u + v * v);
   int nbOfPhotons = computeNbOfPhotons(pBeta, step_length) * aStep->GetTrack()->GetWeight();
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: pBeta " << pBeta << " u/v/w " << u << "/" << v << "/" << w
                                << " step_length " << step_length << " zFib " << zFiber << " nbOfPhotons "
                                << nbOfPhotons;
 #endif
   if (nbOfPhotons < 0) {
     return 0;
   } else if (nbOfPhotons > 0) {
     G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
     const G4TouchableHandle& theTouchable = preStepPoint->GetTouchableHandle();
     G4ThreeVector localprepos =
         theTouchable->GetHistory()->GetTopTransform().TransformPoint(aStep->GetPreStepPoint()->GetPosition());
     G4ThreeVector localpostpos =
         theTouchable->GetHistory()->GetTopTransform().TransformPoint(aStep->GetPostStepPoint()->GetPosition());
 
     double length = sqrt((localpostpos.x() - localprepos.x()) * (localpostpos.x() - localprepos.x()) +
                          (localpostpos.y() - localprepos.y()) * (localpostpos.y() - localprepos.y()));
     double yemit = std::sqrt(fibreR * fibreR - length * length / 4.);
 
     double u_ph = 0, v_ph = 0, w_ph = 0;
     for (int i = 0; i < nbOfPhotons; i++) {
       double rand = G4UniformRand();
       double theta_C = acos(1. / (pBeta * ref_index));
       double phi_C = 2 * M_PI * rand;
       double sinTheta = sin(theta_C);
       double cosTheta = cos(theta_C);
       double cosPhi = cos(phi_C);
       double sinPhi = sin(phi_C);
       //photon momentum
       if (uv < 0.001) {  // aligned with z-axis
         u_ph = sinTheta * cosPhi;
         v_ph = sinTheta * sinPhi;
         w_ph = cosTheta;
       } else {  // general case
         u_ph = uv * cosTheta + sinTheta * cosPhi * w;
         v_ph = sinTheta * sinPhi;
         w_ph = w * cosTheta - sinTheta * cosPhi * uv;
       }
       double r_lambda = G4UniformRand();
       double lambda0 = (lambda1 * lambda2) / (lambda2 - r_lambda * (lambda2 - lambda1));
       double lambda = (lambda0 / cm) * pow(double(10), 7);  // lambda is in nm
       wlini.push_back(lambda);
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: " << i << " lambda " << lambda << " w_ph " << w_ph
                                    << " aperture " << aperture;
 #endif
       // --------------
       double xemit = length * (G4UniformRand() - 0.5);
       double gam = atan2(yemit, xemit);
       double eps = atan2(v_ph, u_ph);
       double sinBeta = sin(gam - eps);
       double rho = sqrt(xemit * xemit + yemit * yemit);
       double sinEta = rho / fibreR * sinBeta;
       double cosEta = sqrt(1. - sinEta * sinEta);
       double sinPsi = sqrt(1. - w_ph * w_ph);
       double cosKsi = cosEta * sinPsi;
 #ifdef EDM_ML_DEBUG
       if (cosKsi < aperturetrapped && w_ph > 0.) {
         edm::LogVerbatim("HFShower") << "HFCherenkov::Trapped photon : " << u_ph << " " << v_ph << " " << w_ph << " "
                                      << xemit << " " << gam << " " << eps << " " << sinBeta << " " << rho << " "
                                      << sinEta << " " << cosEta << " "
                                      << " " << sinPsi << " " << cosKsi;
       } else {
         edm::LogVerbatim("HFShower") << "HFCherenkov::Rejected photon : " << u_ph << " " << v_ph << " " << w_ph << " "
                                      << xemit << " " << gam << " " << eps << " " << sinBeta << " " << rho << " "
                                      << sinEta << " " << cosEta << " "
                                      << " " << sinPsi << " " << cosKsi;
       }
 #endif
       if (cosKsi < aperturetrapped  // photon trapped inside fiber
           && w_ph > 0.              // and moves to PMT
           && sinPsi < sinPsimax) {  // and is not reflected at fiber end
         wltrap.push_back(lambda);
         double prob_HF = 1.0;                 //photon survived in HF
         double a0_inv = 0.1234;               //meter^-1
         double prob_MX = exp(-0.5 * a0_inv);  //light mixer
         if (checkSurvive) {
           double a_inv = a0_inv + 0.14 * pow(dose, 0.30);
           double z_meters = zFiber;
           prob_HF = exp(-z_meters * a_inv);  //photon survived in HF
         }
         rand = G4UniformRand();
 #ifdef EDM_ML_DEBUG
         edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: probHF " << prob_HF << " prob_MX " << prob_MX
                                      << " Random " << rand << " Survive? " << (rand < (prob_HF * prob_MX));
 #endif
         if (rand < (prob_HF * prob_MX)) {  // survived and sent to light mixer
           wlatten.push_back(lambda);
           rand = G4UniformRand();
           double effHEM = computeHEMEff(lambda);
 #ifdef EDM_ML_DEBUG
           edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: w_ph " << w_ph << " effHEM " << effHEM << " Random "
                                        << rand << " Survive? " << (w_ph > 0.997 || (rand < effHEM));
 #endif
           if (w_ph > 0.997 || (rand < effHEM)) {  // survived HEM
             wlhem.push_back(lambda);
             double qEffic = computeQEff(lambda);
             rand = G4UniformRand();
 #ifdef EDM_ML_DEBUG
             edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: qEffic " << qEffic << " Random " << rand
                                          << " Survive? " << (rand < qEffic);
 #endif
             if (rand < qEffic) {  // made photoelectron
               npe_Dose += 1;
               momZ.push_back(w_ph);
               wl.push_back(lambda);
               wlqeff.push_back(lambda);
             }          // made pe
           }            // passed HEM
         }              // passed fiber
       }                // end of  if(w_ph < w_aperture), trapped inside fiber
     }                  // end of ++NbOfPhotons
   }                    // end of if(NbOfPhotons)}
   int npe = npe_Dose;  // Nb of photoelectrons
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: npe " << npe;
 #endif
   return npe;
 }

int HFCherenkov::computeNPEinPMT	(	const G4ParticleDefinition *	pDef,
		double	pBeta,
		double	u,
		double	v,
		double	w,
		double	step_length
	)

Definition at line 221 of file HFCherenkov.cc.

References aperture, clearVectors(), computeNbOfPhotons(), computeQEff(), funct::cos(), mps_fire::i, isApplicable(), lambda1, lambda2, M_PI, momZ, funct::pow(), rand(), ref_index, funct::sin(), mathSSE::sqrt(), wl, wlatten, wlini, wlqeff, and wltrap.

Referenced by HFShowerFibreBundle::getHits(), and HFShowerPMT::getHits().

                                                                                                       {
   clearVectors();
   int npe_ = 0;
   if (!isApplicable(pDef)) {
     return 0;
   }
   if (pBeta < (1 / ref_index) || step_length < 0.0001) {
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: pBeta " << pBeta << " 1/mu " << (1 / ref_index)
                                  << " step_length " << step_length;
 #endif
     return 0;
   }
 
   double uv = sqrt(u * u + v * v);
   int nbOfPhotons = computeNbOfPhotons(pBeta, step_length);
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: pBeta " << pBeta << " u/v/w " << u << "/" << v << "/"
                                << w << " step_length " << step_length << " nbOfPhotons " << nbOfPhotons;
 #endif
   if (nbOfPhotons < 0) {
     return 0;
   } else if (nbOfPhotons > 0) {
     double w_ph = 0;
     for (int i = 0; i < nbOfPhotons; i++) {
       double rand = G4UniformRand();
       double theta_C = acos(1. / (pBeta * ref_index));
       double phi_C = 2 * M_PI * rand;
       double sinTheta = sin(theta_C);
       double cosTheta = cos(theta_C);
       double cosPhi = cos(phi_C);
       //photon momentum
       if (uv < 0.001) {  // aligned with z-axis
         w_ph = cosTheta;
       } else {  // general case
         w_ph = w * cosTheta - sinTheta * cosPhi * uv;
       }
       double r_lambda = G4UniformRand();
       double lambda0 = (lambda1 * lambda2) / (lambda2 - r_lambda * (lambda2 - lambda1));
       double lambda = (lambda0 / cm) * pow(double(10), 7);  // lambda is in nm
       wlini.push_back(lambda);
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: " << i << " lambda " << lambda << " w_ph " << w_ph
                                    << " aperture " << aperture;
 #endif
       if (w_ph > aperture) {  // phton trapped inside PMT glass
         wltrap.push_back(lambda);
         rand = G4UniformRand();
 #ifdef EDM_ML_DEBUG
         edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: Random " << rand << " Survive? " << (rand < 1.);
 #endif
         if (rand < 1.0) {  // survived all the times and sent to photo-cathode
           wlatten.push_back(lambda);
           double qEffic = computeQEff(lambda);  //Quantum efficiency of the PMT
           rand = G4UniformRand();
 #ifdef EDM_ML_DEBUG
           edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: qEffic " << qEffic << " Random " << rand
                                        << " Survive? " << (rand < qEffic);
 #endif
           if (rand < qEffic) {  // made photoelectron
             npe_ += 1;
             momZ.push_back(w_ph);
             wl.push_back(lambda);
             wlqeff.push_back(lambda);
           }  // made pe
         }    // accepted all Cherenkov photons
       }      // end of  if(w_ph < w_aperture), trapped inside glass
     }        // end of ++NbOfPhotons
   }          // end of if(NbOfPhotons)}
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: npe " << npe_;
 #endif
   return npe_;
 }

int HFCherenkov::computeNPhTrapped	(	double	pBeta,
		double	u,
		double	v,
		double	w,
		double	step_length,
		double	zFiber,
		double	Dose,
		int	Npe_Dose
	)

Definition at line 42 of file HFCherenkov.cc.

References apertureTrap, computeNbOfPhotons(), funct::cos(), mps_fire::i, M_PI, rand(), ref_index, funct::sin(), and mathSSE::sqrt().

                                                                                                               {
   if (pBeta < (1 / ref_index) || step_length < 0.0001) {
     return 0;
   }
 
   double uv = sqrt(u * u + v * v);
   int nbOfPhotons = computeNbOfPhotons(pBeta, step_length);
 
   if (nbOfPhotons < 0) {
     return 0;
   } else if (nbOfPhotons > 0) {
     double w_ph = 0;
     for (int i = 0; i < nbOfPhotons; i++) {
       double rand = G4UniformRand();
       double theta_C = acos(1. / (pBeta * ref_index));
       double phi_C = 2 * M_PI * rand;
       double sinTheta = sin(theta_C);
       double cosTheta = cos(theta_C);
       double cosPhi = cos(phi_C);
       //photon momentum
       if (uv < 0.001) {  // aligned with z-axis
         w_ph = cosTheta;
       } else {  // general case
         w_ph = w * cosTheta - sinTheta * cosPhi * uv;
       }
       if (w_ph > apertureTrap) {  // phton trapped inside fiber
         npe_Dose += 1;
       }
     }
   }
   int n_photons = npe_Dose;
   return n_photons;
 }

double HFCherenkov::computeQEff ( double wavelength )

private

Definition at line 349 of file HFCherenkov.cc.

References JetChargeProducer_cfi::exp, patCandidates_cff::func, funct::pow(), UseNewPMT, and y.

Referenced by computeNPE(), and computeNPEinPMT().

                                                  {
   double qeff(0.);
   if (UseNewPMT) {
     if (wavelength <= 350) {
       qeff = 2.45867 * (TMath::Landau(wavelength, 353.820, 59.1324));
     } else if (wavelength > 350 && wavelength < 500) {
       qeff = 0.441989 * exp(-pow((wavelength - 358.371), 2) / (2 * pow((138.277), 2)));
     } else if (wavelength >= 500 && wavelength < 550) {
       qeff = 0.271862 * exp(-pow((wavelength - 491.505), 2) / (2 * pow((47.0418), 2)));
     } else if (wavelength >= 550) {
       qeff = 0.137297 * exp(-pow((wavelength - 520.260), 2) / (2 * pow((75.5023), 2)));
     }
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("HFShower") << "HFCherenkov:: for new PMT : wavelength === " << wavelength << "\tqeff  ===\t"
                                  << qeff;
 #endif
   } else {
     double y = (wavelength - 275.) / 180.;
     double func = 1. / (1. + 250. * pow((y / 5.), 4));
     double qE_R7525 = 0.77 * y * exp(-y) * func;
     qeff = qE_R7525;
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("HFShower") << "HFCherenkov::computeQEff: wavelength " << wavelength << " y/func " << y << "/"
                                  << func << " qeff " << qeff;
     edm::LogVerbatim("HFShower") << "HFCherenkov:: for old PMT : wavelength === " << wavelength << "; qeff = " << qeff;
 #endif
   }
   return qeff;
 }

std::vector< double > HFCherenkov::getMom ( )

Definition at line 327 of file HFCherenkov.cc.

References momZ, and findQualityFiles::v.

Referenced by HFShower::getHits().

                                       {
   std::vector<double> v = momZ;
   return v;
 }

std::vector< double > HFCherenkov::getWL ( )

Definition at line 322 of file HFCherenkov.cc.

References findQualityFiles::v, and wl.

Referenced by HFShower::getHits().

                                      {
   std::vector<double> v = wl;
   return v;
 }

std::vector< double > HFCherenkov::getWLAtten ( )

Definition at line 307 of file HFCherenkov.cc.

References findQualityFiles::v, and wlatten.

                                           {
   std::vector<double> v = wlatten;
   return v;
 }

std::vector< double > HFCherenkov::getWLHEM ( )

Definition at line 312 of file HFCherenkov.cc.

References findQualityFiles::v, and wlhem.

                                         {
   std::vector<double> v = wlhem;
   return v;
 }

std::vector< double > HFCherenkov::getWLIni ( )

Definition at line 297 of file HFCherenkov.cc.

References findQualityFiles::v, and wlini.

                                         {
   std::vector<double> v = wlini;
   return v;
 }

std::vector< double > HFCherenkov::getWLQEff ( )

Definition at line 317 of file HFCherenkov.cc.

References findQualityFiles::v, and wlqeff.

                                          {
   std::vector<double> v = wlqeff;
   return v;
 }

std::vector< double > HFCherenkov::getWLTrap ( )

Definition at line 302 of file HFCherenkov.cc.

References findQualityFiles::v, and wltrap.

                                          {
   std::vector<double> v = wltrap;
   return v;
 }

bool HFCherenkov::isApplicable ( const G4ParticleDefinition * aParticleType )

private

Definition at line 431 of file HFCherenkov.cc.

References tmp.

Referenced by computeNPE(), and computeNPEinPMT().

                                                                         {
   bool tmp = (aParticleType->GetPDGCharge() != 0);
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::isApplicable: aParticleType " << aParticleType->GetParticleName()
                                << " PDGCharge " << aParticleType->GetPDGCharge() << " Result " << tmp;
 #endif
   return tmp;
 }

double HFCherenkov::smearNPE ( G4int Npe )

Definition at line 407 of file HFCherenkov.cc.

References gain, mps_fire::i, and heppy_batch::val.

                                     {
   double pe = 0.;
   if (npe > 0) {
     for (int i = 0; i < npe; ++i) {
       double val = G4Poisson(gain);
       pe += (val / gain) + 0.001 * G4UniformRand();
     }
   }
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HFShower") << "HFCherenkov::smearNPE: npe " << npe << " pe " << pe;
 #endif
   return pe;
 }