Forward/src/ZdcSD.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: ZdcSD.cc
0003 // Date: 03.01
0004 // Description: Sensitive Detector class for Zdc
0005 // Modifications:
0006 ///////////////////////////////////////////////////////////////////////////////
0007 #include <memory>
0008
0009 #include "SimG4CMS/Forward/interface/ZdcSD.h"
0010 #include "SimG4CMS/Forward/interface/ForwardName.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0014 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0015 #include "SimG4Core/Notification/interface/TrackInformation.h"
0016 #include "SimG4CMS/Forward/interface/ZdcNumberingScheme.h"
0017
0018 #include "G4SDManager.hh"
0019 #include "G4Step.hh"
0020 #include "G4Track.hh"
0021 #include "G4VProcess.hh"
0022 #include "G4ios.hh"
0023 #include "G4Cerenkov.hh"
0024 #include "G4ParticleTable.hh"
0025 #include "G4PhysicalConstants.hh"
0026 #include <CLHEP/Units/SystemOfUnits.h>
0027 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0028 #include "Randomize.hh"
0029 #include "G4Poisson.hh"
0030 #include "G4TwoVector.hh"
0031
0032 //#define EDM_ML_DEBUG
0033
0034 ZdcSD::ZdcSD(const std::string& name,
0035              const SensitiveDetectorCatalog& clg,
0036              edm::ParameterSet const& p,
0037              const SimTrackManager* manager)
0038     : CaloSD(name, clg, p, manager) {
0039   edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
0040   useShowerLibrary = m_ZdcSD.getParameter<bool>("UseShowerLibrary");
0041   useShowerHits = m_ZdcSD.getParameter<bool>("UseShowerHits");
0042   zdcHitEnergyCut = m_ZdcSD.getParameter<double>("ZdcHitEnergyCut") * CLHEP::GeV;
0043   thFibDir = m_ZdcSD.getParameter<double>("FiberDirection");
0044   verbosity = m_ZdcSD.getParameter<int>("Verbosity");
0045   int verbn = verbosity / 10;
0046   verbosity %= 10;
0047   numberingScheme = std::make_unique<ZdcNumberingScheme>(verbn);
0048
0049   edm::LogVerbatim("ForwardSim") << "***************************************************\n"
0050                                  << "*                                                 *\n"
0051                                  << "* Constructing a ZdcSD  with name " << name << "   *\n"
0052                                  << "*                                                 *\n"
0053                                  << "***************************************************";
0054
0055   edm::LogVerbatim("ForwardSim") << "\nUse of shower library is set to " << useShowerLibrary
0056                                  << "\nUse of Shower hits method is set to " << useShowerHits;
0057
0058   edm::LogVerbatim("ForwardSim") << "\nEnergy Threshold Cut set to " << zdcHitEnergyCut / CLHEP::GeV << " (GeV)";
0059
0060   if (useShowerLibrary) {
0061     showerLibrary = std::make_unique<ZdcShowerLibrary>(name, p);
0062     setParameterized(true);
0063   } else {
0064     showerLibrary.reset(nullptr);
0065   }
0066 }
0067
0068 void ZdcSD::initRun() {
0069   if (useShowerLibrary) {
0070     G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
0071     showerLibrary->initRun(theParticleTable);
0072   }
0073   hits.clear();
0074 }
0075
0076 bool ZdcSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0077   if (useShowerLibrary)
0078     getFromLibrary(aStep);
0079
0080 #ifdef EDM_ML_DEBUG
0081   edm::LogVerbatim("ZdcSD") << "ZdcSD::" << GetName() << " ID= " << aStep->GetTrack()->GetTrackID()
0082                             << " prID= " << aStep->GetTrack()->GetParentID()
0083                             << " Eprestep= " << aStep->GetPreStepPoint()->GetKineticEnergy()
0084                             << " step= " << aStep->GetStepLength() << " Edep= " << aStep->GetTotalEnergyDeposit();
0085 #endif
0086   if (useShowerHits) {
0087     // check unitID
0088     unsigned int unitID = setDetUnitId(aStep);
0089     if (unitID == 0)
0090       return false;
0091
0092     auto const theTrack = aStep->GetTrack();
0093     uint16_t depth = getDepth(aStep);
0094
0095     double time = theTrack->GetGlobalTime() / CLHEP::nanosecond;
0096     int primaryID = getTrackID(theTrack);
0097     currentID[0].setID(unitID, time, primaryID, depth);
0098     double energy = calculateCherenkovDeposit(aStep);
0099
0100     // Russian Roulette
0101     double wt2 = theTrack->GetWeight();
0102     if (wt2 > 0.0) {
0103       energy *= wt2;
0104     }
0105
0106     if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0107       edepositEM = energy;
0108       edepositHAD = 0;
0109     } else {
0110       edepositEM = 0;
0111       edepositHAD = energy;
0112     }
0113     if (!hitExists(aStep, 0) && energy > 0.) {
0114 #ifdef EDM_ML_DEBUG
0115       G4ThreeVector pre = aStep->GetPreStepPoint()->GetPosition();
0116       edm::LogVerbatim("ZdcSD") << pre.x() << " " << pre.y() << " " << pre.z();
0117 #endif
0118       currentHit[0] = CaloSD::createNewHit(aStep, theTrack, 0);
0119     }
0120   }
0121   return true;
0122 }
0123
0124 bool ZdcSD::getFromLibrary(const G4Step* aStep) {
0125   bool ok = true;
0126
0127   auto const preStepPoint = aStep->GetPreStepPoint();
0128
0129   double etrack = preStepPoint->GetKineticEnergy();
0130   int primaryID = setTrackID(aStep);
0131
0132   hits.clear();
0133
0134   // Reset entry point for new primary
0135   resetForNewPrimary(aStep);
0136
0137   if (etrack >= zdcHitEnergyCut) {
0138     // create hits only if above threshold
0139
0140 #ifdef EDM_ML_DEBUG
0141     auto const theTrack = aStep->GetTrack();
0142     edm::LogVerbatim("ForwardSim") << "----------------New track------------------------------\n"
0143                                    << "Incident EnergyTrack: " << etrack << " MeV \n"
0144                                    << "Zdc Cut Energy for Hits: " << zdcHitEnergyCut << " MeV \n"
0145                                    << "ZdcSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
0146                                    << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0147                                    << etrack << " MeV\n";
0148 #endif
0149     hits.swap(showerLibrary.get()->getHits(aStep, ok));
0150   }
0151
0152   incidentEnergy = etrack;
0153   entrancePoint = preStepPoint->GetPosition();
0154   for (unsigned int i = 0; i < hits.size(); i++) {
0155     posGlobal = hits[i].position;
0156     entranceLocal = hits[i].entryLocal;
0157     double time = hits[i].time;
0158     unsigned int unitID = hits[i].detID;
0159     edepositHAD = hits[i].DeHad;
0160     edepositEM = hits[i].DeEM;
0161     currentID[0].setID(unitID, time, primaryID, 0);
0162     processHit(aStep);
0163
0164 #ifdef EDM_ML_DEBUG
0165     edm::LogVerbatim("ForwardSim") << "ZdcSD: Final Hit number:" << i << "-->"
0166                                    << "New HitID: " << currentHit[0]->getUnitID()
0167                                    << " New Hit trackID: " << currentHit[0]->getTrackID()
0168                                    << " New EM Energy: " << currentHit[0]->getEM() / CLHEP::GeV
0169                                    << " New HAD Energy: " << currentHit[0]->getHadr() / CLHEP::GeV
0170                                    << " New HitEntryPoint: " << currentHit[0]->getEntryLocal()
0171                                    << " New IncidentEnergy: " << currentHit[0]->getIncidentEnergy() / CLHEP::GeV
0172                                    << " New HitPosition: " << posGlobal;
0173 #endif
0174   }
0175   return ok;
0176 }
0177
0178 double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
0179   double NCherPhot = 0.;
0180
0181   // preStepPoint information
0182   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0183
0184   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0185   const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
0186   G4double stepL = aStep->GetStepLength() / CLHEP::cm;
0187   G4double beta = preStepPoint->GetBeta();
0188   G4double charge = preStepPoint->GetCharge();
0189   if (charge == 0.0)
0190     return 0.0;
0191
0192   // theTrack information
0193   G4Track* theTrack = aStep->GetTrack();
0194   G4String particleType = theTrack->GetDefinition()->GetParticleName();
0195   G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0196
0197 #ifdef EDM_ML_DEBUG
0198   const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
0199
0200   // calculations
0201   float costheta =
0202       vert_mom.z() / sqrt(vert_mom.x() * vert_mom.x() + vert_mom.y() * vert_mom.y() + vert_mom.z() * vert_mom.z());
0203   float theta = std::acos(std::min(std::max(costheta, -1.f), 1.f));
0204   float eta = -std::log(std::tan(theta * 0.5f));
0205   float phi = -100.;
0206   if (vert_mom.x() != 0)
0207     phi = std::atan2(vert_mom.y(), vert_mom.x());
0208   if (phi < 0.)
0209     phi += twopi;
0210
0211   // Get the total energy deposit
0212   double stepE = aStep->GetTotalEnergyDeposit();
0213
0214   // postStepPoint information
0215   G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
0216   G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
0217   std::string postnameVolume = ForwardName::getName(postPV->GetName());
0218   std::string nameVolume = preStepPoint->GetPhysicalVolume()->GetName();
0219   edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: \n"
0220                                  << "  preStepPoint: " << nameVolume << "," << stepL << "," << stepE << "," << beta
0221                                  << "," << charge << "\n"
0222                                  << "  postStepPoint: " << postnameVolume << "," << costheta << "," << theta << ","
0223                                  << eta << "," << phi << "," << particleType << " id= " << theTrack->GetTrackID()
0224                                  << " Etot(GeV)= " << theTrack->GetTotalEnergy() / CLHEP::GeV;
0225 #endif
0226   const double bThreshold = 0.67;
0227   if (beta > bThreshold) {
0228 #ifdef EDM_ML_DEBUG
0229     edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit:  pass ";
0230 #endif
0231     const float nMedium = 1.4925;
0232     // float photEnSpectrDL = 10714.285714;
0233     //       photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
0234
0235     const float photEnSpectrDE = 1.24;
0236     // E = 2pi*(1./137.)*(eV*cm/370.)/lambda = 12.389184*(eV*cm)/lambda
0237     // Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm  = 3.01 eV
0238     // Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm  = 1.77 eV
0239     // delE = Emax - Emin = 1.24 eV
0240
0241     const float effPMTandTransport = 0.15;
0242
0243     // Check these values
0244     const float thFullRefl = 23.;
0245     float thFullReflRad = thFullRefl * pi / 180.;
0246
0247     float thFibDirRad = thFibDir * pi / 180.;
0248
0249     // at which theta the point is located:
0250     //   float th1 = hitPoint.theta();
0251
0252     // theta of charged particle in LabRF(hit momentum direction):
0253     float costh = hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
0254     float th = acos(std::min(std::max(costh, -1.f), 1.f));
0255     // just in case (can do both standard ranges of phi):
0256     if (th < 0.)
0257       th += CLHEP::twopi;
0258
0259     // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
0260     float costhcher = 1. / (nMedium * beta);
0261     float thcher = acos(std::min(std::max(costhcher, -1.f), 1.f));
0262
0263     // diff thetas of charged part. and quartz direction in LabRF:
0264     float DelFibPart = std::abs(th - thFibDirRad);
0265
0266     // define real distances:
0267     float d = std::abs(std::tan(th) - std::tan(thFibDirRad));
0268
0269     float a = std::tan(thFibDirRad) + std::tan(std::abs(thFibDirRad - thFullReflRad));
0270     float r = std::tan(th) + std::tan(std::abs(th - thcher));
0271
0272     // define losses d_qz in cone of full reflection inside quartz direction
0273     float d_qz = -1;
0274 #ifdef EDM_ML_DEBUG
0275     float variant = -1;
0276 #endif
0277     // if (d > (r+a))
0278     if (DelFibPart > (thFullReflRad + thcher)) {
0279 #ifdef EDM_ML_DEBUG
0280       variant = 0.;
0281 #endif
0282       d_qz = 0.;
0283     } else {
0284       // if ((DelFibPart + thcher) < thFullReflRad )  [(d+r) < a]
0285       if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
0286 #ifdef EDM_ML_DEBUG
0287         variant = 1.;
0288 #endif
0289         d_qz = 1.;
0290       } else {
0291         // if ((thcher - DelFibPart ) > thFullReflRad )  [(r-d) > a]
0292         if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
0293 #ifdef EDM_ML_DEBUG
0294           variant = 2.;
0295 #endif
0296           d_qz = 0.;
0297         } else {
0298 #ifdef EDM_ML_DEBUG
0299           variant = 3.;  // d_qz is calculated below
0300 #endif
0301           // use crossed length of circles(cone projection) - dC1/dC2 :
0302           float arg_arcos = 0.;
0303           float tan_arcos = 2. * a * d;
0304           if (tan_arcos != 0.)
0305             arg_arcos = (r * r - a * a - d * d) / tan_arcos;
0306           arg_arcos = std::abs(arg_arcos);
0307           float th_arcos = acos(std::min(std::max(arg_arcos, -1.f), 1.f));
0308           d_qz = th_arcos / CLHEP::twopi;
0309           d_qz = std::abs(d_qz);
0310 #ifdef EDM_ML_DEBUG
0311           edm::LogVerbatim("ForwardSim") << "  d_qz: " << r << "," << a << "," << d << " " << tan_arcos << " "
0312                                          << arg_arcos;
0313           edm::LogVerbatim("ForwardSim") << "," << arg_arcos;
0314           edm::LogVerbatim("ForwardSim") << " " << d_qz;
0315           edm::LogVerbatim("ForwardSim") << " " << th_arcos;
0316           edm::LogVerbatim("ForwardSim") << "," << d_qz;
0317 #endif
0318         }
0319       }
0320     }
0321     double meanNCherPhot = 0.;
0322     int poissNCherPhot = 0;
0323     if (d_qz > 0) {
0324       meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepL;
0325
0326       poissNCherPhot = std::max((int)G4Poisson(meanNCherPhot), 0);
0327       NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
0328     }
0329
0330 #ifdef EDM_ML_DEBUG
0331     edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit:  gED: " << stepE << "," << costh << "," << th << ","
0332                                    << costhcher << "," << thcher << "," << DelFibPart << "," << d << "," << a << ","
0333                                    << r << "," << hitPoint << "," << hit_mom << "," << vert_mom << "," << localPoint
0334                                    << "," << charge << "," << beta << "," << stepL << "," << d_qz << "," << variant
0335                                    << "," << meanNCherPhot << "," << poissNCherPhot << "," << NCherPhot;
0336 #endif
0337
0338   } else {
0339     // determine failure mode: beta, charge, and/or nameVolume
0340     if (beta <= bThreshold)
0341       edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail beta=" << beta;
0342   }
0343
0344   return NCherPhot;
0345 }
0346
0347 ///////////////////////////////////////
0348 // Functions added by Oliver Suranyi //
0349 ///////////////////////////////////////
0350
0351 // Constants as global variables
0352 const double RINDEX = 1.47;
0353 const double NA = 0.22;  // Numerical aperture, characteristic of the specific fiber
0354 const double NAperRINDEX = NA / RINDEX;
0355 const double EMAX = 4.79629 /*eV*/;                                    // Maximum energy of PMT sensitivity range
0356 const double EMIN = 1.75715 /*eV*/;                                    // Minimum energy of PMT sensitivity range
0357 const double ALPHA = /*1/137=*/0.0072973525693;                        // Fine structure constant
0358 const double HBARC = 6.582119514E-16 /*eV*s*/ * 2.99792458E8 /*m/s*/;  // hbar * c
0359
0360 // Calculate the Cherenkov deposit corresponding to a G4Step
0361 double ZdcSD::calculateCherenkovDeposit(const G4Step* aStep) {
0362   const G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
0363   G4double charge = pPreStepPoint->GetCharge() / CLHEP::eplus;
0364   if (charge == 0.0 || aStep->GetStepLength() < 1e-9 * CLHEP::mm)
0365     return 0.0;
0366
0367   const G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
0368
0369   G4ThreeVector pre = pPreStepPoint->GetPosition();
0370   G4ThreeVector post = pPostStepPoint->GetPosition();
0371
0372   //Convert step coordinates to local (fiber) coodinates
0373   const G4ThreeVector localPre = setToLocal(pre, pPreStepPoint->GetTouchable());
0374   const G4ThreeVector localPost = setToLocal(post, pPreStepPoint->GetTouchable());
0375
0376   // Calculate the unit direction vector in local coordinates
0377   const G4ThreeVector particleDirection = (localPost - localPre) / (localPost - localPre).mag();
0378
0379   double beta = 0.5 * (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta());
0380   double stepLength = aStep->GetStepLength() / 1000;  // Geant4 stepLength is in "mm"
0381
0382   int nPhotons;  // Number of Cherenkov photons
0383
0384   nPhotons = G4Poisson(calculateMeanNumberOfPhotons(charge, beta, stepLength));
0385
0386   double totalE = 0.0;
0387
0388   for (int i = 0; i < nPhotons; ++i) {
0389     // uniform refractive index in PMT range -> uniform energy distribution
0390     double photonE = EMIN + G4UniformRand() * (EMAX - EMIN);
0391     // UPDATE: taking into account dispersion relation -> energy distribution
0392
0393     if (G4UniformRand() > pmtEfficiency(convertEnergyToWavelength(photonE)))
0394       continue;
0395
0396     double omega = G4UniformRand() * twopi;
0397     double cosTheta = std::min(1.0 / (beta * RINDEX), 1.0);
0398     double sinTheta = std::sqrt((1. - cosTheta) * (1.0 + cosTheta));
0399
0400 #ifdef EDM_ML_DEBUG
0401     edm::LogVerbatim("ZdcSD") << "E_gamma: " << photonE << "\t omega: " << omega << "\t thetaC: " << cosTheta;
0402 #endif
0403     // Calculate momentum direction w.r.t primary particle (z-direction)
0404     double px = photonE * sinTheta * std::cos(omega);
0405     double py = photonE * sinTheta * std::sin(omega);
0406     double pz = photonE * cosTheta;
0407     G4ThreeVector photonMomentum(px, py, pz);
0408
0409 #ifdef EDM_ML_DEBUG
0410     edm::LogVerbatim("ZdcSD") << "pPR = (" << particleDirection.x() << "," << particleDirection.y() << ","
0411                               << particleDirection.z() << ")";
0412     edm::LogVerbatim("ZdcSD") << "pCH = (" << px << "," << py << "," << pz << ")";
0413 #endif
0414     // Rotate to the fiber reference frame
0415     photonMomentum.rotateUz(particleDirection);
0416
0417 #ifdef EDM_ML_DEBUG
0418     edm::LogVerbatim("ZdcSD") << "pLAB = (" << photonMomentum.x() << "," << photonMomentum.y() << ","
0419                               << photonMomentum.z() << ")";
0420 #endif
0421     // Get random position along G4Step
0422     G4ThreeVector photonPosition = localPre + G4UniformRand() * (localPost - localPre);
0423
0424     // 2D vectors to calculate impact position (x*,y*)
0425     G4TwoVector r0(photonPosition);
0426     G4TwoVector v0(photonMomentum);
0427
0428     double R = 0.3; /*mm, fiber radius*/
0429     double R2 = 0.3 * 0.3;
0430
0431     if (r0.mag() < R && photonMomentum.z() < 0.0) {
0432       // 2nd order polynomial coefficients
0433       double a = v0.mag2();
0434       double b = 2.0 * r0 * v0;
0435       double c = r0.mag2() - R2;
0436
0437       if (a < 1E-6)
0438         totalE += 1;  //photonE /*eV*/;
0439       else {
0440         // calculate intersection point - solving 2nd order polynomial
0441         double t = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
0442         G4ThreeVector n(r0.x() + v0.x() * t, r0.y() + v0.y() * t, 0.0);  // surface normal
0443         double cosTheta = (n * photonMomentum) / (n.mag() * photonE);    // cosine of incident angle
0444
0445         if (cosTheta >= NAperRINDEX)  // lightguide condition
0446           totalE += 1;                //photonE /*eV*/;
0447       }
0448     }
0449
0450 #ifdef EDM_ML_DEBUG
0451     edm::LogVerbatim("ZdcSD") << "r = (" << photonPosition.x() << "," << photonPosition.y() << "," << photonPosition.z()
0452                               << ")" << std::endl;
0453 #endif
0454   }
0455
0456 #ifdef EDM_ML_DEBUG
0457   if (nPhotons > 30) {
0458     edm::LogVerbatim("ZdcSD") << totalE;
0459
0460     if (totalE > 0)
0461       edm::LogVerbatim("ZdcSD") << pre.x() << " " << pre.y() << " " << pre.z() << " " << totalE;
0462   }
0463 #endif
0464   return totalE;
0465 }
0466
0467 // Calculate mean number of Cherenkov photons in the sensitivity range (300-650 nm)
0468 // for a given step length for a particle with given charge and beta
0469 double ZdcSD::calculateMeanNumberOfPhotons(double charge, double beta, double stepLength) {
0470   // Return mean number of Cherenkov photons
0471   return (ALPHA * charge * charge * stepLength) / HBARC * (EMAX - EMIN) * (1.0 - 1.0 / (beta * beta * RINDEX * RINDEX));
0472 }
0473
0474 // Evaluate photon pdf
0475 double ZdcSD::photonEnergyDist(double charge, double beta, double E) {
0476   const std::vector<double> ENERGY_TAB{1.75715, 1.81902, 1.88311, 1.94944, 2.0183,  2.08939, 2.16302, 2.23919,
0477                                        2.31789, 2.39954, 2.48416, 2.57175, 2.66232, 2.75643, 2.85349, 2.95411,
0478                                        3.05756, 3.16528, 3.2774,  3.39218, 3.5123,  3.6359,  3.76394, 3.89642,
0479                                        4.03332, 4.17596, 4.32302, 4.47596, 4.63319, 4.79629};
0480
0481   const std::vector<double> RINDEX_TAB{1.45517, 1.45572, 1.45631, 1.45693, 1.45758, 1.45826, 1.45899, 1.45976,
0482                                        1.46057, 1.46144, 1.46238, 1.46337, 1.46444, 1.46558, 1.4668,  1.46812,
0483                                        1.46953, 1.47105, 1.4727,  1.47447, 1.4764,  1.47847, 1.48071, 1.48315,
0484                                        1.48579, 1.48868, 1.49182, 1.49526, 1.499,   1.5031};
0485   double rIndex = evaluateFunction(ENERGY_TAB, RINDEX_TAB, E);
0486   return (ALPHA * charge * charge) / HBARC * (1.0 - 1.0 / (beta * beta * rIndex * rIndex));
0487 }
0488
0489 // Generate a photon with the given minimum energy accourding to the energy distribution
0490 double ZdcSD::generatePhotonEnergy(double charge, double beta, double Emin) {
0491   double photonE;
0492
0493   // Use rejection method
0494   do {
0495     photonE = G4UniformRand() * (EMAX - Emin) + Emin;
0496   } while (G4UniformRand() > photonEnergyDist(photonE, charge, beta) / photonEnergyDist(EMAX, charge, beta));
0497
0498   return photonE;
0499 }
0500
0501 // Calculate the integral: int_Emin^Emax 1/n(E)^2 dE
0502 // The integral values are tabulated
0503 double ZdcSD::calculateN2InvIntegral(double Emin) {
0504   // Hardcoded minimum photon energy table (eV)
0505   const std::vector<double> EMIN_TAB{1.75715, 1.81902, 1.88311, 1.94944, 2.0183,  2.08939, 2.16302, 2.23919,
0506                                      2.31789, 2.39954, 2.48416, 2.57175, 2.66232, 2.75643, 2.85349, 2.95411,
0507                                      3.05756, 3.16528, 3.2774,  3.39218, 3.5123,  3.6359,  3.76394, 3.89642,
0508                                      4.03332, 4.17596, 4.32302, 4.47596, 4.63319};
0509
0510   // Hardcoded integral values
0511   const std::vector<double> INTEGRAL_TAB{1.39756,  1.36835,  1.33812,  1.30686,  1.27443,  1.24099,  1.20638,  1.17061,
0512                                          1.1337,   1.09545,  1.05586,  1.01493,  0.972664, 0.928815, 0.883664, 0.836938,
0513                                          0.788988, 0.739157, 0.687404, 0.634547, 0.579368, 0.522743, 0.464256, 0.40393,
0514                                          0.341808, 0.27732,  0.211101, 0.142536, 0.0723891};
0515   return evaluateFunction(EMIN_TAB, INTEGRAL_TAB, Emin);
0516 }
0517
0518 double ZdcSD::pmtEfficiency(double lambda) {
0519   // Hardcoded wavelength values (nm)
0520   const std::vector<double> LAMBDA_TAB{263.27, 265.98, 268.69, 271.39, 273.20, 275.90, 282.22, 282.22, 293.04,
0521                                        308.38, 325.52, 346.26, 367.91, 392.27, 417.53, 440.98, 463.53, 484.28,
0522                                        502.32, 516.75, 528.48, 539.30, 551.93, 564.56, 574.48, 584.41, 595.23,
0523                                        606.96, 616.88, 625.00, 632.22, 637.63, 642.14, 647.55, 652.96, 656.57,
0524                                        661.08, 666.49, 669.20, 673.71, 677.32, 680.93, 686.34, 692.65};
0525
0526   // Hardcoded quantum efficiency values
0527   const std::vector<double> EFF_TAB{2.215,  2.860,  3.659,  4.724,  5.989,  7.734,  9.806,  9.806,  12.322,
0528                                     15.068, 17.929, 20.570, 22.963, 24.050, 23.847, 22.798, 20.445, 18.003,
0529                                     15.007, 12.282, 9.869,  7.858,  6.373,  5.121,  4.077,  3.276,  2.562,
0530                                     2.077,  1.669,  1.305,  1.030,  0.805,  0.629,  0.492,  0.388,  0.303,
0531                                     0.239,  0.187,  0.144,  0.113,  0.088,  0.069,  0.054,  0.043};
0532   //double efficiency = evaluateFunction(LAMBDA_TAB,EFF_TAB,lambda);
0533
0534   // Using linear interpolation to calculate efficiency
0535   for (int i = 0; i < 44 - 1; i++) {
0536     if (lambda > LAMBDA_TAB[i] && lambda < LAMBDA_TAB[i + 1]) {
0537       double a = (EFF_TAB[i] - EFF_TAB[i + 1]) / (LAMBDA_TAB[i] - LAMBDA_TAB[i + 1]);
0538       double b = EFF_TAB[i] - a * LAMBDA_TAB[i];
0539       return (a * lambda + b) / 100.0;
0540     }
0541   }
0542
0543   return 0;
0544 }
0545
0546 // Evaluate a function given by set of datapoints
0547 // Linear interpolation is used to calculate function value between datapoints
0548 double ZdcSD::evaluateFunction(const std::vector<double>& X, const std::vector<double>& Y, double x) {
0549   for (unsigned int i = 0; i < X.size() - 1; i++) {
0550     if (x > X[i] && x < X[i + 1]) {
0551 #ifdef EDM_ML_DEBUG
0552       edm::LogVerbatim("ZdcSD") << X[i] << "\t" << Y[i] << "\t" << X[i + 1] << "\t" << Y[i + 1] << "\t" << x << "\t"
0553                                 << linearInterpolation(X[i], Y[i], X[i + 1], Y[i + 1], x);
0554 #endif
0555       return linearInterpolation(X[i], Y[i], X[i + 1], Y[i + 1], x);
0556     }
0557   }
0558
0559   if (fabs(X[0] - x) < 1E-8)
0560     return Y[0];
0561   else if (fabs(X[X.size() - 1] - x) < 1E-8)
0562     return Y[X.size() - 1];
0563   else
0564     return NAN;
0565 }
0566
0567 // Do linear interpolation between two points
0568 double ZdcSD::linearInterpolation(double x1, double y1, double x2, double y2, double z) {
0569   if (x1 < x2)
0570     return y1 + (z - x1) * (y2 - y1) / (x2 - x1);
0571   else
0572     return y2 + (z - x2) * (y1 - y2) / (x1 - x2);
0573 }
0574
0575 // Energy (eV) to wavelength (nm) conversion
0576 double ZdcSD::convertEnergyToWavelength(double energy) { return (1240.0 / energy); }
0577
0578 /////////////////////////////////////////////////////////////////////
0579
0580 uint32_t ZdcSD::setDetUnitId(const G4Step* aStep) {
0581   return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
0582 }
0583
0584 int ZdcSD::setTrackID(const G4Step* aStep) {
0585   auto const theTrack = aStep->GetTrack();
0586   TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
0587   int primaryID = trkInfo->getIDonCaloSurface();
0588   if (primaryID == 0) {
0589 #ifdef EDM_ML_DEBUG
0590     auto const preStepPoint = aStep->GetPreStepPoint();
0591     double etrack = preStepPoint->GetKineticEnergy();
0592     edm::LogVerbatim("ZdcSD") << "ZdcSD: Problem with primaryID **** set by force to TkID **** "
0593                               << theTrack->GetTrackID() << " E " << etrack;
0594 #endif
0595     primaryID = theTrack->GetTrackID();
0596   }
0597   if (primaryID != previousID[0].trackID())
0598     resetForNewPrimary(aStep);
0599   return primaryID;
0600 }