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/TrackInformation.h"
0015
0016 #include "G4SDManager.hh"
0017 #include "G4Step.hh"
0018 #include "G4Track.hh"
0019 #include "G4VProcess.hh"
0020 #include "G4ios.hh"
0021 #include "G4Cerenkov.hh"
0022 #include "G4ParticleTable.hh"
0023 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0024 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0025 #include "Randomize.hh"
0026 #include "G4Poisson.hh"
0027
0028 //#define EDM_ML_DEBUG
0029
0030 ZdcSD::ZdcSD(const std::string& name,
0031              const SensitiveDetectorCatalog& clg,
0032              edm::ParameterSet const& p,
0033              const SimTrackManager* manager)
0034     : CaloSD(name, clg, p, manager) {
0035   edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
0036   useShowerLibrary = m_ZdcSD.getParameter<bool>("UseShowerLibrary");
0037   useShowerHits = m_ZdcSD.getParameter<bool>("UseShowerHits");
0038   zdcHitEnergyCut = m_ZdcSD.getParameter<double>("ZdcHitEnergyCut") * GeV;
0039   thFibDir = m_ZdcSD.getParameter<double>("FiberDirection");
0040   verbosity = m_ZdcSD.getParameter<int>("Verbosity");
0041   int verbn = verbosity / 10;
0042   verbosity %= 10;
0043   setNumberingScheme(new ZdcNumberingScheme(verbn));
0044
0045   edm::LogVerbatim("ForwardSim") << "***************************************************\n"
0046                                  << "*                                                 *\n"
0047                                  << "* Constructing a ZdcSD  with name " << name << "   *\n"
0048                                  << "*                                                 *\n"
0049                                  << "***************************************************";
0050
0051   edm::LogVerbatim("ForwardSim") << "\nUse of shower library is set to " << useShowerLibrary
0052                                  << "\nUse of Shower hits method is set to " << useShowerHits;
0053
0054   edm::LogVerbatim("ForwardSim") << "\nEnergy Threshold Cut set to " << zdcHitEnergyCut / GeV << " (GeV)";
0055
0056   if (useShowerLibrary) {
0057     showerLibrary = std::make_unique<ZdcShowerLibrary>(name, p);
0058     setParameterized(true);
0059   } else {
0060     showerLibrary.reset(nullptr);
0061   }
0062 }
0063
0064 void ZdcSD::initRun() { hits.clear(); }
0065
0066 bool ZdcSD::getFromLibrary(const G4Step* aStep) {
0067   bool ok = true;
0068
0069   auto const preStepPoint = aStep->GetPreStepPoint();
0070
0071   double etrack = preStepPoint->GetKineticEnergy();
0072   int primaryID = setTrackID(aStep);
0073
0074   hits.clear();
0075
0076   // Reset entry point for new primary
0077   resetForNewPrimary(aStep);
0078
0079   if (etrack >= zdcHitEnergyCut) {
0080     // create hits only if above threshold
0081
0082 #ifdef EDM_ML_DEBUG
0083     auto const theTrack = aStep->GetTrack();
0084     edm::LogVerbatim("ForwardSim") << "----------------New track------------------------------\n"
0085                                    << "Incident EnergyTrack: " << etrack << " MeV \n"
0086                                    << "Zdc Cut Energy for Hits: " << zdcHitEnergyCut << " MeV \n"
0087                                    << "ZdcSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
0088                                    << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0089                                    << etrack << " MeV\n";
0090 #endif
0091     hits.swap(showerLibrary.get()->getHits(aStep, ok));
0092   }
0093
0094   incidentEnergy = etrack;
0095   entrancePoint = preStepPoint->GetPosition();
0096   for (unsigned int i = 0; i < hits.size(); i++) {
0097     posGlobal = hits[i].position;
0098     entranceLocal = hits[i].entryLocal;
0099     double time = hits[i].time;
0100     unsigned int unitID = hits[i].detID;
0101     edepositHAD = hits[i].DeHad;
0102     edepositEM = hits[i].DeEM;
0103     currentID.setID(unitID, time, primaryID, 0);
0104     processHit(aStep);
0105
0106 #ifdef EDM_ML_DEBUG
0107     edm::LogVerbatim("ForwardSim") << "ZdcSD: Final Hit number:" << i << "-->"
0108                                    << "New HitID: " << currentHit->getUnitID()
0109                                    << " New Hit trackID: " << currentHit->getTrackID()
0110                                    << " New EM Energy: " << currentHit->getEM() / GeV
0111                                    << " New HAD Energy: " << currentHit->getHadr() / GeV
0112                                    << " New HitEntryPoint: " << currentHit->getEntryLocal()
0113                                    << " New IncidentEnergy: " << currentHit->getIncidentEnergy() / GeV
0114                                    << " New HitPosition: " << posGlobal;
0115 #endif
0116   }
0117   return ok;
0118 }
0119
0120 double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
0121   double NCherPhot = 0.;
0122
0123   // preStepPoint information
0124   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0125   G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
0126   std::string nameVolume = ForwardName::getName(currentPV->GetName());
0127
0128   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0129   const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
0130   G4double stepL = aStep->GetStepLength() / cm;
0131   G4double beta = preStepPoint->GetBeta();
0132   G4double charge = preStepPoint->GetCharge();
0133
0134   // theTrack information
0135   G4Track* theTrack = aStep->GetTrack();
0136   G4String particleType = theTrack->GetDefinition()->GetParticleName();
0137   G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0138
0139 #ifdef EDM_ML_DEBUG
0140   const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
0141
0142   // calculations
0143   float costheta =
0144       vert_mom.z() / sqrt(vert_mom.x() * vert_mom.x() + vert_mom.y() * vert_mom.y() + vert_mom.z() * vert_mom.z());
0145   float theta = std::acos(std::min(std::max(costheta, -1.f), 1.f));
0146   float eta = -std::log(std::tan(theta * 0.5f));
0147   float phi = -100.;
0148   if (vert_mom.x() != 0)
0149     phi = std::atan2(vert_mom.y(), vert_mom.x());
0150   if (phi < 0.)
0151     phi += twopi;
0152
0153   // Get the total energy deposit
0154   double stepE = aStep->GetTotalEnergyDeposit();
0155
0156   // postStepPoint information
0157   G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
0158   G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
0159   std::string postnameVolume = ForwardName::getName(postPV->GetName());
0160   edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: \n"
0161                                  << "  preStepPoint: " << nameVolume << "," << stepL << "," << stepE << "," << beta
0162                                  << "," << charge << "\n"
0163                                  << "  postStepPoint: " << postnameVolume << "," << costheta << "," << theta << ","
0164                                  << eta << "," << phi << "," << particleType << " id= " << theTrack->GetTrackID()
0165                                  << " Etot(GeV)= " << theTrack->GetTotalEnergy() / GeV;
0166 #endif
0167   const double bThreshold = 0.67;
0168   if ((beta > bThreshold) && (charge != 0) && (nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber")) {
0169 #ifdef EDM_ML_DEBUG
0170     edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit:  pass ";
0171 #endif
0172     const float nMedium = 1.4925;
0173     // float photEnSpectrDL = 10714.285714;
0174     //       photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
0175
0176     const float photEnSpectrDE = 1.24;
0177     // E = 2pi*(1./137.)*(eV*cm/370.)/lambda = 12.389184*(eV*cm)/lambda
0178     // Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm  = 3.01 eV
0179     // Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm  = 1.77 eV
0180     // delE = Emax - Emin = 1.24 eV
0181
0182     const float effPMTandTransport = 0.15;
0183
0184     // Check these values
0185     const float thFullRefl = 23.;
0186     float thFullReflRad = thFullRefl * pi / 180.;
0187
0188     float thFibDirRad = thFibDir * pi / 180.;
0189
0190     // at which theta the point is located:
0191     //   float th1 = hitPoint.theta();
0192
0193     // theta of charged particle in LabRF(hit momentum direction):
0194     float costh = hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
0195     float th = acos(std::min(std::max(costh, -1.f), 1.f));
0196     // just in case (can do both standard ranges of phi):
0197     if (th < 0.)
0198       th += twopi;
0199
0200     // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
0201     float costhcher = 1. / (nMedium * beta);
0202     float thcher = acos(std::min(std::max(costhcher, -1.f), 1.f));
0203
0204     // diff thetas of charged part. and quartz direction in LabRF:
0205     float DelFibPart = std::abs(th - thFibDirRad);
0206
0207     // define real distances:
0208     float d = std::abs(std::tan(th) - std::tan(thFibDirRad));
0209
0210     float a = std::tan(thFibDirRad) + std::tan(std::abs(thFibDirRad - thFullReflRad));
0211     float r = std::tan(th) + std::tan(std::abs(th - thcher));
0212
0213     // define losses d_qz in cone of full reflection inside quartz direction
0214     float d_qz = -1;
0215 #ifdef EDM_ML_DEBUG
0216     float variant = -1;
0217 #endif
0218     // if (d > (r+a))
0219     if (DelFibPart > (thFullReflRad + thcher)) {
0220 #ifdef EDM_ML_DEBUG
0221       variant = 0.;
0222 #endif
0223       d_qz = 0.;
0224     } else {
0225       // if ((DelFibPart + thcher) < thFullReflRad )  [(d+r) < a]
0226       if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
0227 #ifdef EDM_ML_DEBUG
0228         variant = 1.;
0229 #endif
0230         d_qz = 1.;
0231       } else {
0232         // if ((thcher - DelFibPart ) > thFullReflRad )  [(r-d) > a]
0233         if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
0234 #ifdef EDM_ML_DEBUG
0235           variant = 2.;
0236 #endif
0237           d_qz = 0.;
0238         } else {
0239 #ifdef EDM_ML_DEBUG
0240           variant = 3.;  // d_qz is calculated below
0241 #endif
0242           // use crossed length of circles(cone projection) - dC1/dC2 :
0243           float arg_arcos = 0.;
0244           float tan_arcos = 2. * a * d;
0245           if (tan_arcos != 0.)
0246             arg_arcos = (r * r - a * a - d * d) / tan_arcos;
0247           arg_arcos = std::abs(arg_arcos);
0248           float th_arcos = acos(std::min(std::max(arg_arcos, -1.f), 1.f));
0249           d_qz = th_arcos / twopi;
0250           d_qz = std::abs(d_qz);
0251 #ifdef EDM_ML_DEBUG
0252           edm::LogVerbatim("ForwardSim") << "  d_qz: " << r << "," << a << "," << d << " " << tan_arcos << " "
0253                                          << arg_arcos;
0254           edm::LogVerbatim("ForwardSim") << "," << arg_arcos;
0255           edm::LogVerbatim("ForwardSim") << " " << d_qz;
0256           edm::LogVerbatim("ForwardSim") << " " << th_arcos;
0257           edm::LogVerbatim("ForwardSim") << "," << d_qz;
0258 #endif
0259         }
0260       }
0261     }
0262     double meanNCherPhot = 0.;
0263     int poissNCherPhot = 0;
0264     if (d_qz > 0) {
0265       meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepL;
0266
0267       poissNCherPhot = std::max((int)G4Poisson(meanNCherPhot), 0);
0268       NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
0269     }
0270
0271 #ifdef EDM_ML_DEBUG
0272     edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit:  gED: " << stepE << "," << costh << "," << th << ","
0273                                    << costhcher << "," << thcher << "," << DelFibPart << "," << d << "," << a << ","
0274                                    << r << "," << hitPoint << "," << hit_mom << "," << vert_mom << "," << localPoint
0275                                    << "," << charge << "," << beta << "," << stepL << "," << d_qz << "," << variant
0276                                    << "," << meanNCherPhot << "," << poissNCherPhot << "," << NCherPhot;
0277 #endif
0278     // --constants-----------------
0279     // << "," << photEnSpectrDE
0280     // << "," << nMedium
0281     // << "," << bThreshold
0282     // << "," << thFibDirRad
0283     // << "," << thFullReflRad
0284     // << "," << effPMTandTransport
0285     // --other variables-----------
0286     // << "," << curprocess
0287     // << "," << nameProcess
0288     // << "," << name
0289     // << "," << rad
0290     // << "," << mat
0291
0292   } else {
0293     // determine failure mode: beta, charge, and/or nameVolume
0294     if (beta <= bThreshold)
0295       edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail beta=" << beta;
0296     if (charge == 0)
0297       edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail charge=0";
0298     if (!(nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber"))
0299       edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail nv=" << nameVolume;
0300   }
0301
0302   return NCherPhot;
0303 }
0304
0305 uint32_t ZdcSD::setDetUnitId(const G4Step* aStep) {
0306   return (numberingScheme.get() == nullptr ? 0 : numberingScheme.get()->getUnitID(aStep));
0307 }
0308
0309 void ZdcSD::setNumberingScheme(ZdcNumberingScheme* scheme) {
0310   if (scheme != nullptr) {
0311     edm::LogVerbatim("ForwardSim") << "ZdcSD: updates numbering scheme for " << GetName();
0312     numberingScheme.reset(scheme);
0313   }
0314 }