#include <HFCherenkov.h>

Public Member Functions
void	clearVectors ()

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

int	computeNPEinPMT (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 15 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::LogInfo("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 41 of file HFCherenkov.cc.

41 {}

Member Function Documentation

void HFCherenkov::clearVectors ( )

Definition at line 444 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 399 of file HFCherenkov.cc.

References LogDebug.

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 DebugLog
   LogDebug("HFShower") << "HFCherenkov::computeHEMEff: wavelength "
                << wavelength << " hEMEff " << hEMEff;
 #endif
   return hEMEff;
 }

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

private

Definition at line 344 of file HFCherenkov.cc.

References alpha, beta, lambda1, lambda2, LogDebug, 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 DebugLog
   LogDebug("HFShower") << "HFCherenkov::computeNbOfPhotons: StepLength " 
                << step_length << " theta_C " << theta_C 
                << " lambdaDiff " << lambdaDiff
                << " cherenPhPerLength " << cherenPhPerLength 
                << " Photons " << d_NOfPhotons << " " << nbOfPhotons;
 #endif
   return nbOfPhotons;
 }

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

Definition at line 79 of file HFCherenkov.cc.

References aperture, aperturetrapped, checkSurvive, clearVectors(), computeHEMEff(), computeNbOfPhotons(), computeQEff(), funct::cos(), create_public_lumi_plots::exp, fibreR, gam, i, isApplicable(), lambda1, lambda2, LogDebug, 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 DebugLog
     LogDebug("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 DebugLog
   LogDebug("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();
     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 DebugLog
       LogDebug("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 DebugLog
       if (cosKsi < aperturetrapped && w_ph>0.) {
     LogDebug("HFShower") << "HFCherenkov::Trapped photon : " << u_ph << " "
                  << v_ph << " " << w_ph << " " << xemit << " "
                  << gam << " " << eps << " " << sinBeta << " "
                  << rho << " " << sinEta << " " << cosEta << " "
                  << " " << sinPsi << " " << cosKsi;
       } else {
     LogDebug("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 DebugLog
     LogDebug("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 DebugLog
       LogDebug("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 DebugLog
         LogDebug("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 DebugLog
   LogDebug("HFShower") << "HFCherenkov::computeNPE: npe " << npe;
 #endif
   return npe;
 }

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

Definition at line 226 of file HFCherenkov.cc.

References aperture, clearVectors(), computeNbOfPhotons(), computeQEff(), funct::cos(), i, isApplicable(), lambda1, lambda2, LogDebug, 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 DebugLog
     LogDebug("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 DebugLog
   LogDebug("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 DebugLog
       LogDebug("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 DebugLog
     LogDebug("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);
       rand = G4UniformRand();
       double qEffic = computeQEff(lambda);//Quantum efficiency of the PMT
       rand = G4UniformRand();
 #ifdef DebugLog
       LogDebug("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 DebugLog
   LogDebug("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 43 of file HFCherenkov.cc.

References apertureTrap, computeNbOfPhotons(), funct::cos(), 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 364 of file HFCherenkov.cc.

References create_public_lumi_plots::exp, LogDebug, funct::pow(), UseNewPMT, and detailsBasic3DVector::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 DebugLog
     LogDebug("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 DebugLog
     LogDebug("HFShower") << "HFCherenkov:: for old PMT : wavelength === "
              << wavelength << "; qeff = " << qeff;
 #endif
   }
 
 #ifdef DebugLog
   LogDebug("HFShower") << "HFCherenkov::computeQEff: wavelength " << wavelength
                << " y/func " << y << "/" << func << " qeff " << qeff;
 #endif
   return qeff;
 }

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

Definition at line 339 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 334 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 319 of file HFCherenkov.cc.

References findQualityFiles::v, and wlatten.

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

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

Definition at line 324 of file HFCherenkov.cc.

References findQualityFiles::v, and wlhem.

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

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

Definition at line 309 of file HFCherenkov.cc.

References findQualityFiles::v, and wlini.

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

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

Definition at line 329 of file HFCherenkov.cc.

References findQualityFiles::v, and wlqeff.

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

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

Definition at line 314 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 455 of file HFCherenkov.cc.

References LogDebug, and tmp.

Referenced by computeNPE(), and computeNPEinPMT().

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

double HFCherenkov::smearNPE ( G4int Npe )

Definition at line 429 of file HFCherenkov.cc.

References gain, i, and LogDebug.

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