#include <ZdcSD.h>

Inheritance diagram for ZdcSD:

Public Member Functions
bool	ProcessHits (G4Step step, G4TouchableHistory tHistory) override

uint32_t	setDetUnitId (const G4Step *step) override

	ZdcSD (const std::string &, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)

	~ZdcSD () override=default

Public Member Functions inherited from CaloSD
	CaloSD (const std::string &aSDname, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, float timeSlice=1., bool ignoreTkID=false, const std::string &newcolname="")

void	clear () override

void	clearHits () override

void	DrawAll () override

void	EndOfEvent (G4HCofThisEvent *eventHC) override

void	fillHits (edm::PCaloHitContainer &, const std::string &) override

void	Initialize (G4HCofThisEvent *HCE) override

bool	isItFineCalo (const G4VTouchable *touch)

void	newCollection (const std::string &name, edm::ParameterSet const &p)

void	PrintAll () override

bool	ProcessHits (G4GFlashSpot aSpot, G4TouchableHistory ) override

void	reset () override

	~CaloSD () override

Public Member Functions inherited from SensitiveCaloDetector
	SensitiveCaloDetector (const std::string &iname, const SensitiveDetectorCatalog &clg, const std::string &newcollname="")

Public Member Functions inherited from SensitiveDetector
void	EndOfEvent (G4HCofThisEvent *eventHC) override

const std::vector< std::string > &	getNames () const

void	Initialize (G4HCofThisEvent *eventHC) override

bool	isCaloSD () const

	SensitiveDetector (const std::string &iname, const SensitiveDetectorCatalog &, bool calo, const std::string &newcollname="")

	~SensitiveDetector () override

Public Member Functions inherited from Observer< const BeginOfRun *>
	Observer ()

void	slotForUpdate (const BeginOfRun * iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const BeginOfEvent *>
	Observer ()

void	slotForUpdate (const BeginOfEvent * iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const BeginOfTrack *>
	Observer ()

void	slotForUpdate (const BeginOfTrack * iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfTrack *>
	Observer ()

void	slotForUpdate (const EndOfTrack * iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfEvent *>
	Observer ()

void	slotForUpdate (const EndOfEvent * iT)

virtual	~Observer ()

Protected Member Functions
double	calculateCherenkovDeposit (const G4Step *)

double	calculateMeanNumberOfPhotons (double, double, double)

double	calculateN2InvIntegral (double)

double	convertEnergyToWavelength (double)

double	evaluateFunction (const std::vector< double > &, const std::vector< double > &, double)

double	generatePhotonEnergy (double, double, double)

double	getEnergyDeposit (const G4Step *) override

bool	getFromLibrary (const G4Step *) override

void	initRun () override

double	linearInterpolation (double, double, double, double, double)

double	photonEnergyDist (double, double, double)

double	pmtEfficiency (double)

int	setTrackID (const G4Step *step) override

Protected Member Functions inherited from CaloSD
bool	checkHit (int k=0)

CaloG4Hit *	createNewHit (const G4Step , const G4Track , int k)

virtual void	endEvent ()

virtual double	EnergyCorrected (const G4Step &step, const G4Track *)

virtual bool	filterHit (CaloG4Hit *, double)

unsigned int	findBoundaryCrossingParent (const G4Track *track, bool markParentAsSaveable=true)

double	getAttenuation (const G4Step *aStep, double birk1, double birk2, double birk3) const

virtual uint16_t	getDepth (const G4Step *)

int	getNumberOfHits (int k=0)

double	getResponseWt (const G4Track *, int k=0)

virtual int	getTrackID (const G4Track *)

bool	hitExists (const G4Step *, int k)

void	ignoreRejection ()

virtual void	initEvent (const BeginOfEvent *)

void	printDetectorLevels (const G4VTouchable *) const

void	processHit (const G4Step *step)

virtual void	processSecondHit (const G4Step , const G4Track )

void	resetForNewPrimary (const G4Step *)

void	setNumberCheckedHits (int val, int k=0)

void	setParameterized (bool val)

G4ThreeVector	setToGlobal (const G4ThreeVector &, const G4VTouchable *) const

G4ThreeVector	setToLocal (const G4ThreeVector &, const G4VTouchable *) const

void	setUseMap (bool val)

std::string	shortreprID (const CaloHitID &ID)

std::string	shortreprID (const CaloG4Hit *hit)

void	update (const BeginOfRun *) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const BeginOfEvent *) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const BeginOfTrack *trk) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const EndOfTrack *trk) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const ::EndOfEvent *) override

void	updateHit (CaloG4Hit *, int k)

Protected Member Functions inherited from SensitiveDetector
TrackInformation *	cmsTrackInformation (const G4Track *aTrack)

Local3DPoint	ConvertToLocal3DPoint (const G4ThreeVector &point) const

Local3DPoint	FinalStepPosition (const G4Step *step, coordinates) const

Local3DPoint	InitialStepPosition (const G4Step *step, coordinates) const

Local3DPoint	LocalPostStepPosition (const G4Step *step) const

Local3DPoint	LocalPreStepPosition (const G4Step *step) const

void	NaNTrap (const G4Step *step) const

void	setNames (const std::vector< std::string > &)

Protected Member Functions inherited from Observer< const EndOfEvent *>
virtual void	update (const EndOfEvent *)=0
	This routine will be called when the appropriate signal arrives. More...

Private Attributes
std::vector< ZdcShowerLibrary::Hit >	hits

std::unique_ptr< ZdcNumberingScheme >	numberingScheme

std::unique_ptr< ZdcShowerLibrary >	showerLibrary

double	thFibDir

bool	useShowerHits

bool	useShowerLibrary

int	verbosity

double	zdcHitEnergyCut

Additional Inherited Members
Protected Types inherited from SensitiveDetector
enum	coordinates { WorldCoordinates, LocalCoordinates }

Static Protected Member Functions inherited from CaloSD
static std::string	printableDecayChain (const std::vector< unsigned int > &decayChain)

Protected Attributes inherited from CaloSD
std::string	collName_ [2]

CaloG4Hit *	currentHit [2]

CaloHitID	currentID [2]

std::string	detName_

float	edepositEM

float	edepositHAD

double	eminHit

double	energyCut

G4ThreeVector	entranceLocal

G4ThreeVector	entrancePoint

bool	forceSave

std::vector< std::string >	hcn_

float	incidentEnergy

double	kmaxIon

double	kmaxNeutron

double	kmaxProton

int	nHC_

G4ThreeVector	posGlobal

CaloHitID	previousID [2]

bool	suppressHeavy

double	tmaxHit

std::vector< int >	useResMap_

Detailed Description

Definition at line 15 of file ZdcSD.h.

Constructor & Destructor Documentation

◆ ZdcSD()

ZdcSD::ZdcSD	(	const std::string &	name,
		const SensitiveDetectorCatalog &	clg,
		edm::ParameterSet const &	p,
		const SimTrackManager *	manager
	)

Definition at line 34 of file ZdcSD.cc.

References edm::ParameterSet::getParameter(), Skims_PA_cff::name, numberingScheme, AlCaHLTBitMon_ParallelJobs::p, CaloSD::setParameterized(), showerLibrary, thFibDir, useShowerHits, useShowerLibrary, verbosity, and zdcHitEnergyCut.

     : CaloSD(name, clg, p, manager) {
   edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
   useShowerLibrary = m_ZdcSD.getParameter<bool>("UseShowerLibrary");
   useShowerHits = m_ZdcSD.getParameter<bool>("UseShowerHits");
   zdcHitEnergyCut = m_ZdcSD.getParameter<double>("ZdcHitEnergyCut") * CLHEP::GeV;
   thFibDir = m_ZdcSD.getParameter<double>("FiberDirection");
   verbosity = m_ZdcSD.getParameter<int>("Verbosity");
   int verbn = verbosity / 10;
   verbosity %= 10;
   numberingScheme = std::make_unique<ZdcNumberingScheme>(verbn);
 
   edm::LogVerbatim("ForwardSim") << "***************************************************\n"
                                  << "*                                                 *\n"
                                  << "* Constructing a ZdcSD  with name " << name << "   *\n"
                                  << "*                                                 *\n"
                                  << "***************************************************";
 
   edm::LogVerbatim("ForwardSim") << "\nUse of shower library is set to " << useShowerLibrary
                                  << "\nUse of Shower hits method is set to " << useShowerHits;
 
   edm::LogVerbatim("ForwardSim") << "\nEnergy Threshold Cut set to " << zdcHitEnergyCut / CLHEP::GeV << " (GeV)";
 
   if (useShowerLibrary) {
     showerLibrary = std::make_unique<ZdcShowerLibrary>(name, p);
     setParameterized(true);
   } else {
     showerLibrary.reset(nullptr);
   }
 }

◆ ~ZdcSD()

ZdcSD::~ZdcSD ( )

overridedefault

Member Function Documentation

◆ calculateCherenkovDeposit()

double ZdcSD::calculateCherenkovDeposit ( const G4Step * aStep )

protected

Definition at line 361 of file ZdcSD.cc.

References a, b, HLT_2023v12_cff::beta, HltBtagPostValidation_cff::c, calculateMeanNumberOfPhotons(), ALCARECOTkAlJpsiMuMu_cff::charge, convertEnergyToWavelength(), funct::cos(), MillePedeFileConverter_cfg::e, EMAX, EMIN, mps_fire::i, SiStripPI::min, dqmiodumpmetadata::n, NAperRINDEX, run3scouting_cff::nPhotons, pmtEfficiency(), findAndChange::post, findAndChange::pre, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, dttmaxenums::R, RINDEX, CaloSD::setToLocal(), funct::sin(), mathSSE::sqrt(), submitPVValidationJobs::t, and V0Monitor_cff::v0.

Referenced by ProcessHits().

                                                            {
   const G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
   G4double charge = pPreStepPoint->GetCharge() / CLHEP::eplus;
   if (charge == 0.0 || aStep->GetStepLength() < 1e-9 * CLHEP::mm)
     return 0.0;
 
   const G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
 
   G4ThreeVector pre = pPreStepPoint->GetPosition();
   G4ThreeVector post = pPostStepPoint->GetPosition();
 
   //Convert step coordinates to local (fiber) coodinates
   const G4ThreeVector localPre = setToLocal(pre, pPreStepPoint->GetTouchable());
   const G4ThreeVector localPost = setToLocal(post, pPreStepPoint->GetTouchable());
 
   // Calculate the unit direction vector in local coordinates
   const G4ThreeVector particleDirection = (localPost - localPre) / (localPost - localPre).mag();
 
   double beta = 0.5 * (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta());
   double stepLength = aStep->GetStepLength() / 1000;  // Geant4 stepLength is in "mm"
 
   int nPhotons;  // Number of Cherenkov photons
 
   nPhotons = G4Poisson(calculateMeanNumberOfPhotons(charge, beta, stepLength));
 
   double totalE = 0.0;
 
   for (int i = 0; i < nPhotons; ++i) {
     // uniform refractive index in PMT range -> uniform energy distribution
     double photonE = EMIN + G4UniformRand() * (EMAX - EMIN);
     // UPDATE: taking into account dispersion relation -> energy distribution
 
     if (G4UniformRand() > pmtEfficiency(convertEnergyToWavelength(photonE)))
       continue;
 
     double omega = G4UniformRand() * twopi;
     double cosTheta = std::min(1.0 / (beta * RINDEX), 1.0);
     double sinTheta = std::sqrt((1. - cosTheta) * (1.0 + cosTheta));
 
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ZdcSD") << "E_gamma: " << photonE << "\t omega: " << omega << "\t thetaC: " << cosTheta;
 #endif
     // Calculate momentum direction w.r.t primary particle (z-direction)
     double px = photonE * sinTheta * std::cos(omega);
     double py = photonE * sinTheta * std::sin(omega);
     double pz = photonE * cosTheta;
     G4ThreeVector photonMomentum(px, py, pz);
 
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ZdcSD") << "pPR = (" << particleDirection.x() << "," << particleDirection.y() << ","
                               << particleDirection.z() << ")";
     edm::LogVerbatim("ZdcSD") << "pCH = (" << px << "," << py << "," << pz << ")";
 #endif
     // Rotate to the fiber reference frame
     photonMomentum.rotateUz(particleDirection);
 
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ZdcSD") << "pLAB = (" << photonMomentum.x() << "," << photonMomentum.y() << ","
                               << photonMomentum.z() << ")";
 #endif
     // Get random position along G4Step
     G4ThreeVector photonPosition = localPre + G4UniformRand() * (localPost - localPre);
 
     // 2D vectors to calculate impact position (x*,y*)
     G4TwoVector r0(photonPosition);
     G4TwoVector v0(photonMomentum);
 
     double R = 0.3; /*mm, fiber radius*/
     double R2 = 0.3 * 0.3;
 
     if (r0.mag() < R && photonMomentum.z() < 0.0) {
       // 2nd order polynomial coefficients
       double a = v0.mag2();
       double b = 2.0 * r0 * v0;
       double c = r0.mag2() - R2;
 
       if (a < 1E-6)
         totalE += 1;  //photonE /*eV*/;
       else {
         // calculate intersection point - solving 2nd order polynomial
         double t = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
         G4ThreeVector n(r0.x() + v0.x() * t, r0.y() + v0.y() * t, 0.0);  // surface normal
         double cosTheta = (n * photonMomentum) / (n.mag() * photonE);    // cosine of incident angle
 
         if (cosTheta >= NAperRINDEX)  // lightguide condition
           totalE += 1;                //photonE /*eV*/;
       }
     }
 
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ZdcSD") << "r = (" << photonPosition.x() << "," << photonPosition.y() << "," << photonPosition.z()
                               << ")" << std::endl;
 #endif
   }
 
 #ifdef EDM_ML_DEBUG
   if (nPhotons > 30) {
     edm::LogVerbatim("ZdcSD") << totalE;
 
     if (totalE > 0)
       edm::LogVerbatim("ZdcSD") << pre.x() << " " << pre.y() << " " << pre.z() << " " << totalE;
   }
 #endif
   return totalE;
 }

◆ calculateMeanNumberOfPhotons()

double ZdcSD::calculateMeanNumberOfPhotons	(	double	charge,
		double	beta,
		double	stepLength
	)

protected

Definition at line 469 of file ZdcSD.cc.

References ALPHA, HLT_2023v12_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, EMAX, EMIN, HBARC, and RINDEX.

Referenced by calculateCherenkovDeposit().

                                                                                         {
   // Return mean number of Cherenkov photons
   return (ALPHA * charge * charge * stepLength) / HBARC * (EMAX - EMIN) * (1.0 - 1.0 / (beta * beta * RINDEX * RINDEX));
 }

◆ calculateN2InvIntegral()

double ZdcSD::calculateN2InvIntegral ( double Emin )

protected

Definition at line 503 of file ZdcSD.cc.

References evaluateFunction().

                                                 {
   // Hardcoded minimum photon energy table (eV)
   const std::vector<double> EMIN_TAB{1.75715, 1.81902, 1.88311, 1.94944, 2.0183,  2.08939, 2.16302, 2.23919,
                                      2.31789, 2.39954, 2.48416, 2.57175, 2.66232, 2.75643, 2.85349, 2.95411,
                                      3.05756, 3.16528, 3.2774,  3.39218, 3.5123,  3.6359,  3.76394, 3.89642,
                                      4.03332, 4.17596, 4.32302, 4.47596, 4.63319};
 
   // Hardcoded integral values
   const std::vector<double> INTEGRAL_TAB{1.39756,  1.36835,  1.33812,  1.30686,  1.27443,  1.24099,  1.20638,  1.17061,
                                          1.1337,   1.09545,  1.05586,  1.01493,  0.972664, 0.928815, 0.883664, 0.836938,
                                          0.788988, 0.739157, 0.687404, 0.634547, 0.579368, 0.522743, 0.464256, 0.40393,
                                          0.341808, 0.27732,  0.211101, 0.142536, 0.0723891};
   return evaluateFunction(EMIN_TAB, INTEGRAL_TAB, Emin);
 }

◆ convertEnergyToWavelength()

double ZdcSD::convertEnergyToWavelength ( double energy )

protected

Definition at line 576 of file ZdcSD.cc.

References hcalRecHitTable_cff::energy.

Referenced by calculateCherenkovDeposit().

576 { return (1240.0 / energy); }

hcalRecHitTable_cff.energy

energy

Definition: hcalRecHitTable_cff.py:13

◆ evaluateFunction()

double ZdcSD::evaluateFunction	(	const std::vector< double > &	X,
		const std::vector< double > &	Y,
		double	x
	)

protected

Definition at line 548 of file ZdcSD.cc.

References mps_fire::i, linearInterpolation(), x, X, and beamSpotPI::Y.

Referenced by calculateN2InvIntegral(), and photonEnergyDist().

                                                                                                {
   for (unsigned int i = 0; i < X.size() - 1; i++) {
     if (x > X[i] && x < X[i + 1]) {
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("ZdcSD") << X[i] << "\t" << Y[i] << "\t" << X[i + 1] << "\t" << Y[i + 1] << "\t" << x << "\t"
                                 << linearInterpolation(X[i], Y[i], X[i + 1], Y[i + 1], x);
 #endif
       return linearInterpolation(X[i], Y[i], X[i + 1], Y[i + 1], x);
     }
   }
 
   if (fabs(X[0] - x) < 1E-8)
     return Y[0];
   else if (fabs(X[X.size() - 1] - x) < 1E-8)
     return Y[X.size() - 1];
   else
     return NAN;
 }

◆ generatePhotonEnergy()

double ZdcSD::generatePhotonEnergy	(	double	charge,
		double	beta,
		double	Emin
	)

protected

Definition at line 490 of file ZdcSD.cc.

References HLT_2023v12_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, EMAX, and photonEnergyDist().

                                                                           {
   double photonE;
 
   // Use rejection method
   do {
     photonE = G4UniformRand() * (EMAX - Emin) + Emin;
   } while (G4UniformRand() > photonEnergyDist(photonE, charge, beta) / photonEnergyDist(EMAX, charge, beta));
 
   return photonE;
 }

◆ getEnergyDeposit()

double ZdcSD::getEnergyDeposit ( const G4Step * aStep )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 178 of file ZdcSD.cc.

References a, funct::abs(), HLT_2023v12_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, ztail::d, PVValHelper::eta, f, ForwardName::getName(), dqm-mbProfile::log, SiStripPI::max, SiStripPI::min, PbPb_ZMuSkimMuonDPG_cff::particleType, phi, pi, alignCSCRings::r, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, funct::tan(), theta(), and thFibDir.

                                                   {
   double NCherPhot = 0.;
 
   // preStepPoint information
   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
 
   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
   const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
   G4double stepL = aStep->GetStepLength() / cm;
   G4double beta = preStepPoint->GetBeta();
   G4double charge = preStepPoint->GetCharge();
   if (charge == 0.0)
     return 0.0;
 
   // theTrack information
   G4Track* theTrack = aStep->GetTrack();
   G4String particleType = theTrack->GetDefinition()->GetParticleName();
   G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
 
 #ifdef EDM_ML_DEBUG
   const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
 
   // calculations
   float costheta =
       vert_mom.z() / sqrt(vert_mom.x() * vert_mom.x() + vert_mom.y() * vert_mom.y() + vert_mom.z() * vert_mom.z());
   float theta = std::acos(std::min(std::max(costheta, -1.f), 1.f));
   float eta = -std::log(std::tan(theta * 0.5f));
   float phi = -100.;
   if (vert_mom.x() != 0)
     phi = std::atan2(vert_mom.y(), vert_mom.x());
   if (phi < 0.)
     phi += twopi;
 
   // Get the total energy deposit
   double stepE = aStep->GetTotalEnergyDeposit();
 
   // postStepPoint information
   G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
   G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
   std::string postnameVolume = ForwardName::getName(postPV->GetName());
   std::string nameVolume = preStepPoint->GetPhysicalVolume()->GetName();
   edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: \n"
                                  << "  preStepPoint: " << nameVolume << "," << stepL << "," << stepE << "," << beta
                                  << "," << charge << "\n"
                                  << "  postStepPoint: " << postnameVolume << "," << costheta << "," << theta << ","
                                  << eta << "," << phi << "," << particleType << " id= " << theTrack->GetTrackID()
                                  << " Etot(GeV)= " << theTrack->GetTotalEnergy() / CLHEP::GeV;
 #endif
   const double bThreshold = 0.67;
   if (beta > bThreshold) {
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit:  pass ";
 #endif
     const float nMedium = 1.4925;
     // float photEnSpectrDL = 10714.285714;
     //       photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
 
     const float photEnSpectrDE = 1.24;
     // E = 2pi*(1./137.)*(eV*cm/370.)/lambda = 12.389184*(eV*cm)/lambda
     // Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm  = 3.01 eV
     // Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm  = 1.77 eV
     // delE = Emax - Emin = 1.24 eV
 
     const float effPMTandTransport = 0.15;
 
     // Check these values
     const float thFullRefl = 23.;
     float thFullReflRad = thFullRefl * pi / 180.;
 
     float thFibDirRad = thFibDir * pi / 180.;
 
     // at which theta the point is located:
     //   float th1 = hitPoint.theta();
 
     // theta of charged particle in LabRF(hit momentum direction):
     float costh = hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
     float th = acos(std::min(std::max(costh, -1.f), 1.f));
     // just in case (can do both standard ranges of phi):
     if (th < 0.)
       th += CLHEP::twopi;
 
     // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
     float costhcher = 1. / (nMedium * beta);
     float thcher = acos(std::min(std::max(costhcher, -1.f), 1.f));
 
     // diff thetas of charged part. and quartz direction in LabRF:
     float DelFibPart = std::abs(th - thFibDirRad);
 
     // define real distances:
     float d = std::abs(std::tan(th) - std::tan(thFibDirRad));
 
     float a = std::tan(thFibDirRad) + std::tan(std::abs(thFibDirRad - thFullReflRad));
     float r = std::tan(th) + std::tan(std::abs(th - thcher));
 
     // define losses d_qz in cone of full reflection inside quartz direction
     float d_qz = -1;
 #ifdef EDM_ML_DEBUG
     float variant = -1;
 #endif
     // if (d > (r+a))
     if (DelFibPart > (thFullReflRad + thcher)) {
 #ifdef EDM_ML_DEBUG
       variant = 0.;
 #endif
       d_qz = 0.;
     } else {
       // if ((DelFibPart + thcher) < thFullReflRad )  [(d+r) < a]
       if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
 #ifdef EDM_ML_DEBUG
         variant = 1.;
 #endif
         d_qz = 1.;
       } else {
         // if ((thcher - DelFibPart ) > thFullReflRad )  [(r-d) > a]
         if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
 #ifdef EDM_ML_DEBUG
           variant = 2.;
 #endif
           d_qz = 0.;
         } else {
 #ifdef EDM_ML_DEBUG
           variant = 3.;  // d_qz is calculated below
 #endif
           // use crossed length of circles(cone projection) - dC1/dC2 :
           float arg_arcos = 0.;
           float tan_arcos = 2. * a * d;
           if (tan_arcos != 0.)
             arg_arcos = (r * r - a * a - d * d) / tan_arcos;
           arg_arcos = std::abs(arg_arcos);
           float th_arcos = acos(std::min(std::max(arg_arcos, -1.f), 1.f));
           d_qz = th_arcos / CLHEP::twopi;
           d_qz = std::abs(d_qz);
 #ifdef EDM_ML_DEBUG
           edm::LogVerbatim("ForwardSim") << "  d_qz: " << r << "," << a << "," << d << " " << tan_arcos << " "
                                          << arg_arcos;
           edm::LogVerbatim("ForwardSim") << "," << arg_arcos;
           edm::LogVerbatim("ForwardSim") << " " << d_qz;
           edm::LogVerbatim("ForwardSim") << " " << th_arcos;
           edm::LogVerbatim("ForwardSim") << "," << d_qz;
 #endif
         }
       }
     }
     double meanNCherPhot = 0.;
     int poissNCherPhot = 0;
     if (d_qz > 0) {
       meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepL;
 
       poissNCherPhot = std::max((int)G4Poisson(meanNCherPhot), 0);
       NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
     }
 
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit:  gED: " << stepE << "," << costh << "," << th << ","
                                    << costhcher << "," << thcher << "," << DelFibPart << "," << d << "," << a << ","
                                    << r << "," << hitPoint << "," << hit_mom << "," << vert_mom << "," << localPoint
                                    << "," << charge << "," << beta << "," << stepL << "," << d_qz << "," << variant
                                    << "," << meanNCherPhot << "," << poissNCherPhot << "," << NCherPhot;
 #endif
 
   } else {
     // determine failure mode: beta, charge, and/or nameVolume
     if (beta <= bThreshold)
       edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail beta=" << beta;
   }
 
   return NCherPhot;
 }

◆ getFromLibrary()

bool ZdcSD::getFromLibrary ( const G4Step * aStep )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 124 of file ZdcSD.cc.

References CaloSD::currentHit, CaloSD::currentID, CaloSD::edepositEM, CaloSD::edepositHAD, CaloSD::entranceLocal, CaloSD::entrancePoint, CaloG4Hit::getEM(), CaloG4Hit::getEntryLocal(), CaloG4Hit::getHadr(), CaloG4Hit::getIncidentEnergy(), CaloG4Hit::getTrackID(), CaloG4Hit::getUnitID(), hits, mps_fire::i, CaloSD::incidentEnergy, convertSQLiteXML::ok, CaloSD::posGlobal, CaloSD::processHit(), CaloSD::resetForNewPrimary(), CaloHitID::setID(), setTrackID(), showerLibrary, hcalRecHitTable_cff::time, and zdcHitEnergyCut.

Referenced by ProcessHits().

                                               {
   bool ok = true;
 
   auto const preStepPoint = aStep->GetPreStepPoint();
 
   double etrack = preStepPoint->GetKineticEnergy();
   int primaryID = setTrackID(aStep);
 
   hits.clear();
 
   // Reset entry point for new primary
   resetForNewPrimary(aStep);
 
   if (etrack >= zdcHitEnergyCut) {
     // create hits only if above threshold
 
 #ifdef EDM_ML_DEBUG
     auto const theTrack = aStep->GetTrack();
     edm::LogVerbatim("ForwardSim") << "----------------New track------------------------------\n"
                                    << "Incident EnergyTrack: " << etrack << " MeV \n"
                                    << "Zdc Cut Energy for Hits: " << zdcHitEnergyCut << " MeV \n"
                                    << "ZdcSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
                                    << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
                                    << etrack << " MeV\n";
 #endif
     hits.swap(showerLibrary.get()->getHits(aStep, ok));
   }
 
   incidentEnergy = etrack;
   entrancePoint = preStepPoint->GetPosition();
   for (unsigned int i = 0; i < hits.size(); i++) {
     posGlobal = hits[i].position;
     entranceLocal = hits[i].entryLocal;
     double time = hits[i].time;
     unsigned int unitID = hits[i].detID;
     edepositHAD = hits[i].DeHad;
     edepositEM = hits[i].DeEM;
     currentID[0].setID(unitID, time, primaryID, 0);
     processHit(aStep);
 
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ForwardSim") << "ZdcSD: Final Hit number:" << i << "-->"
                                    << "New HitID: " << currentHit[0]->getUnitID()
                                    << " New Hit trackID: " << currentHit[0]->getTrackID()
                                    << " New EM Energy: " << currentHit[0]->getEM() / CLHEP::GeV
                                    << " New HAD Energy: " << currentHit[0]->getHadr() / CLHEP::GeV
                                    << " New HitEntryPoint: " << currentHit[0]->getEntryLocal()
                                    << " New IncidentEnergy: " << currentHit[0]->getIncidentEnergy() / CLHEP::GeV
                                    << " New HitPosition: " << posGlobal;
 #endif
   }
   return ok;
 }

◆ initRun()

void ZdcSD::initRun ( )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 68 of file ZdcSD.cc.

References hits, showerLibrary, and useShowerLibrary.

                     {
   if (useShowerLibrary) {
     G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
     showerLibrary->initRun(theParticleTable);
   }
   hits.clear();
 }

◆ linearInterpolation()

double ZdcSD::linearInterpolation	(	double	x1,
		double	y1,
		double	x2,
		double	y2,
		double	z
	)

protected

Definition at line 568 of file ZdcSD.cc.

References testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, testProducerWithPsetDescEmpty_cfi::y1, testProducerWithPsetDescEmpty_cfi::y2, and z.

Referenced by evaluateFunction().

                                                                                       {
   if (x1 < x2)
     return y1 + (z - x1) * (y2 - y1) / (x2 - x1);
   else
     return y2 + (z - x2) * (y1 - y2) / (x1 - x2);
 }

◆ photonEnergyDist()

double ZdcSD::photonEnergyDist	(	double	charge,
		double	beta,
		double	E
	)

protected

Definition at line 475 of file ZdcSD.cc.

References ALPHA, HLT_2023v12_cff::beta, ALCARECOTkAlJpsiMuMu_cff::charge, evaluateFunction(), and HBARC.

Referenced by generatePhotonEnergy().

                                                                    {
   const std::vector<double> ENERGY_TAB{1.75715, 1.81902, 1.88311, 1.94944, 2.0183,  2.08939, 2.16302, 2.23919,
                                        2.31789, 2.39954, 2.48416, 2.57175, 2.66232, 2.75643, 2.85349, 2.95411,
                                        3.05756, 3.16528, 3.2774,  3.39218, 3.5123,  3.6359,  3.76394, 3.89642,
                                        4.03332, 4.17596, 4.32302, 4.47596, 4.63319, 4.79629};
 
   const std::vector<double> RINDEX_TAB{1.45517, 1.45572, 1.45631, 1.45693, 1.45758, 1.45826, 1.45899, 1.45976,
                                        1.46057, 1.46144, 1.46238, 1.46337, 1.46444, 1.46558, 1.4668,  1.46812,
                                        1.46953, 1.47105, 1.4727,  1.47447, 1.4764,  1.47847, 1.48071, 1.48315,
                                        1.48579, 1.48868, 1.49182, 1.49526, 1.499,   1.5031};
   double rIndex = evaluateFunction(ENERGY_TAB, RINDEX_TAB, E);
   return (ALPHA * charge * charge) / HBARC * (1.0 - 1.0 / (beta * beta * rIndex * rIndex));
 }

◆ pmtEfficiency()

double ZdcSD::pmtEfficiency ( double lambda )

protected

Definition at line 518 of file ZdcSD.cc.

References a, b, and mps_fire::i.

Referenced by calculateCherenkovDeposit().

                                          {
   // Hardcoded wavelength values (nm)
   const std::vector<double> LAMBDA_TAB{263.27, 265.98, 268.69, 271.39, 273.20, 275.90, 282.22, 282.22, 293.04,
                                        308.38, 325.52, 346.26, 367.91, 392.27, 417.53, 440.98, 463.53, 484.28,
                                        502.32, 516.75, 528.48, 539.30, 551.93, 564.56, 574.48, 584.41, 595.23,
                                        606.96, 616.88, 625.00, 632.22, 637.63, 642.14, 647.55, 652.96, 656.57,
                                        661.08, 666.49, 669.20, 673.71, 677.32, 680.93, 686.34, 692.65};
 
   // Hardcoded quantum efficiency values
   const std::vector<double> EFF_TAB{2.215,  2.860,  3.659,  4.724,  5.989,  7.734,  9.806,  9.806,  12.322,
                                     15.068, 17.929, 20.570, 22.963, 24.050, 23.847, 22.798, 20.445, 18.003,
                                     15.007, 12.282, 9.869,  7.858,  6.373,  5.121,  4.077,  3.276,  2.562,
                                     2.077,  1.669,  1.305,  1.030,  0.805,  0.629,  0.492,  0.388,  0.303,
                                     0.239,  0.187,  0.144,  0.113,  0.088,  0.069,  0.054,  0.043};
   //double efficiency = evaluateFunction(LAMBDA_TAB,EFF_TAB,lambda);
 
   // Using linear interpolation to calculate efficiency
   for (int i = 0; i < 44 - 1; i++) {
     if (lambda > LAMBDA_TAB[i] && lambda < LAMBDA_TAB[i + 1]) {
       double a = (EFF_TAB[i] - EFF_TAB[i + 1]) / (LAMBDA_TAB[i] - LAMBDA_TAB[i + 1]);
       double b = EFF_TAB[i] - a * LAMBDA_TAB[i];
       return (a * lambda + b) / 100.0;
     }
   }
 
   return 0;
 }

◆ ProcessHits()

bool ZdcSD::ProcessHits	(	G4Step *	step,
		G4TouchableHistory *	tHistory
	)

overridevirtual

Reimplemented from CaloSD.

Definition at line 76 of file ZdcSD.cc.

References calculateCherenkovDeposit(), CaloSD::createNewHit(), CaloSD::currentHit, CaloSD::currentID, hcalRecHitTable_cff::depth, CaloSD::edepositEM, CaloSD::edepositHAD, hcalRecHitTable_cff::energy, CaloSD::getDepth(), getFromLibrary(), CaloSD::getTrackID(), CaloSD::hitExists(), G4TrackToParticleID::isGammaElectronPositron(), findAndChange::pre, setDetUnitId(), CaloHitID::setID(), hcalRecHitTable_cff::time, useShowerHits, and useShowerLibrary.

                                                           {
   if (useShowerLibrary)
     getFromLibrary(aStep);
 
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("ZdcSD") << "ZdcSD::" << GetName() << " ID= " << aStep->GetTrack()->GetTrackID()
                             << " prID= " << aStep->GetTrack()->GetParentID()
                             << " Eprestep= " << aStep->GetPreStepPoint()->GetKineticEnergy()
                             << " step= " << aStep->GetStepLength() << " Edep= " << aStep->GetTotalEnergyDeposit();
 #endif
   if (useShowerHits) {
     // check unitID
     unsigned int unitID = setDetUnitId(aStep);
     if (unitID == 0)
       return false;
 
     auto const theTrack = aStep->GetTrack();
     uint16_t depth = getDepth(aStep);
 
     double time = theTrack->GetGlobalTime() / nanosecond;
     int primaryID = getTrackID(theTrack);
     currentID[0].setID(unitID, time, primaryID, depth);
     double energy = calculateCherenkovDeposit(aStep);
 
     // Russian Roulette
     double wt2 = theTrack->GetWeight();
     if (wt2 > 0.0) {
       energy *= wt2;
     }
 
     if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
       edepositEM = energy;
       edepositHAD = 0;
     } else {
       edepositEM = 0;
       edepositHAD = energy;
     }
     if (!hitExists(aStep, 0) && energy > 0.) {
 #ifdef EDM_ML_DEBUG
       G4ThreeVector pre = aStep->GetPreStepPoint()->GetPosition();
       edm::LogVerbatim("ZdcSD") << pre.x() << " " << pre.y() << " " << pre.z();
 #endif
       currentHit[0] = CaloSD::createNewHit(aStep, theTrack, 0);
     }
   }
   return true;
 }

◆ setDetUnitId()

uint32_t ZdcSD::setDetUnitId ( const G4Step * step )

overridevirtual

Implements CaloSD.

Definition at line 580 of file ZdcSD.cc.

References numberingScheme.

Referenced by ProcessHits().

                                                 {
   return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
 }

◆ setTrackID()

int ZdcSD::setTrackID ( const G4Step * step )

overrideprotectedvirtual

Reimplemented from CaloSD.

Definition at line 584 of file ZdcSD.cc.

References TrackInformation::getIDonCaloSurface(), CaloSD::previousID, and CaloSD::resetForNewPrimary().

Referenced by getFromLibrary().

                                          {
   auto const theTrack = aStep->GetTrack();
   TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
   int primaryID = trkInfo->getIDonCaloSurface();
   if (primaryID == 0) {
 #ifdef EDM_ML_DEBUG
     auto const preStepPoint = aStep->GetPreStepPoint();
     double etrack = preStepPoint->GetKineticEnergy();
     edm::LogVerbatim("ZdcSD") << "ZdcSD: Problem with primaryID **** set by force to TkID **** "
                               << theTrack->GetTrackID() << " E " << etrack;
 #endif
     primaryID = theTrack->GetTrackID();
   }
   if (primaryID != previousID[0].trackID())
     resetForNewPrimary(aStep);
   return primaryID;
 }

Member Data Documentation

◆ hits

std::vector<ZdcShowerLibrary::Hit> ZdcSD::hits

private

Definition at line 49 of file ZdcSD.h.

Referenced by getFromLibrary(), and initRun().

◆ numberingScheme

std::unique_ptr<ZdcNumberingScheme> ZdcSD::numberingScheme

private

Definition at line 47 of file ZdcSD.h.

Referenced by setDetUnitId(), and ZdcSD().

◆ showerLibrary

std::unique_ptr<ZdcShowerLibrary> ZdcSD::showerLibrary

private

Definition at line 46 of file ZdcSD.h.

Referenced by getFromLibrary(), initRun(), and ZdcSD().

◆ thFibDir

double ZdcSD::thFibDir

private

Definition at line 44 of file ZdcSD.h.

Referenced by getEnergyDeposit(), and ZdcSD().

◆ useShowerHits

bool ZdcSD::useShowerHits

private

Definition at line 43 of file ZdcSD.h.

Referenced by ProcessHits(), and ZdcSD().

◆ useShowerLibrary

bool ZdcSD::useShowerLibrary

private

Definition at line 43 of file ZdcSD.h.

Referenced by initRun(), ProcessHits(), and ZdcSD().

◆ verbosity

int ZdcSD::verbosity

private

Definition at line 42 of file ZdcSD.h.

Referenced by ZdcSD().

◆ zdcHitEnergyCut

double ZdcSD::zdcHitEnergyCut

private

Definition at line 45 of file ZdcSD.h.

Referenced by getFromLibrary(), and ZdcSD().

Public Member Functions

Protected Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ZdcSD()

◆ ~ZdcSD()

Member Function Documentation

◆ calculateCherenkovDeposit()

◆ calculateMeanNumberOfPhotons()

◆ calculateN2InvIntegral()

◆ convertEnergyToWavelength()

◆ evaluateFunction()

◆ generatePhotonEnergy()

◆ getEnergyDeposit()

◆ getFromLibrary()

◆ initRun()

◆ linearInterpolation()

◆ photonEnergyDist()

◆ pmtEfficiency()

◆ ProcessHits()

◆ setDetUnitId()

◆ setTrackID()

Member Data Documentation

◆ hits

◆ numberingScheme

◆ showerLibrary

◆ thFibDir

◆ useShowerHits

◆ useShowerLibrary

◆ verbosity

◆ zdcHitEnergyCut