Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-05-10 02:21:19

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: CastorSD.cc
0003 // Date: 02.04
0004 // UpDate: 07.04 - C3TF & C4TF semi-trapezoid added
0005 // Description: Sensitive Detector class for Castor
0006 ///////////////////////////////////////////////////////////////////////////////
0007 
0008 #include "SimG4Core/Notification/interface/TrackInformation.h"
0009 #include "SimG4Core/Notification/interface/TrackInformationExtractor.h"
0010 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0011 
0012 #include "SimG4CMS/Forward/interface/CastorSD.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 
0015 #include "G4SDManager.hh"
0016 #include "G4Step.hh"
0017 #include "G4Track.hh"
0018 #include "G4VProcess.hh"
0019 
0020 #include "G4ios.hh"
0021 #include "G4Cerenkov.hh"
0022 #include "G4LogicalVolumeStore.hh"
0023 
0024 #include <CLHEP/Units/SystemOfUnits.h>
0025 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0026 #include "Randomize.hh"
0027 #include "G4Poisson.hh"
0028 
0029 //#define EDM_ML_DEBUG
0030 
0031 CastorSD::CastorSD(const std::string& name,
0032                    const SensitiveDetectorCatalog& clg,
0033                    edm::ParameterSet const& p,
0034                    const SimTrackManager* manager)
0035     : CaloSD(name, clg, p, manager),
0036       numberingScheme(nullptr),
0037       lvC3EF(nullptr),
0038       lvC3HF(nullptr),
0039       lvC4EF(nullptr),
0040       lvC4HF(nullptr),
0041       lvCAST(nullptr) {
0042   edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD");
0043   useShowerLibrary = m_CastorSD.getParameter<bool>("useShowerLibrary");
0044   energyThresholdSL = m_CastorSD.getParameter<double>("minEnergyInGeVforUsingSLibrary");
0045   energyThresholdSL = energyThresholdSL * CLHEP::GeV;  //  Convert GeV => MeV
0046 
0047   non_compensation_factor = m_CastorSD.getParameter<double>("nonCompensationFactor");
0048 
0049   if (useShowerLibrary) {
0050     showerLibrary = new CastorShowerLibrary(name, p);
0051     setParameterized(true);
0052   }
0053   setNumberingScheme(new CastorNumberingScheme());
0054 
0055   edm::LogVerbatim("ForwardSim") << "********************************************************\n"
0056                                  << "* Constructing a CastorSD  with name " << GetName() << "\n"
0057                                  << "* Using Castor Shower Library: " << useShowerLibrary << "\n"
0058                                  << "********************************************************";
0059 
0060   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0061   for (auto lv : *lvs) {
0062     if (strcmp((lv->GetName()).c_str(), "C3EF") == 0) {
0063       lvC3EF = lv;
0064     }
0065     if (strcmp((lv->GetName()).c_str(), "C3HF") == 0) {
0066       lvC3HF = lv;
0067     }
0068     if (strcmp((lv->GetName()).c_str(), "C4EF") == 0) {
0069       lvC4EF = lv;
0070     }
0071     if (strcmp((lv->GetName()).c_str(), "C4HF") == 0) {
0072       lvC4HF = lv;
0073     }
0074     if (strcmp((lv->GetName()).c_str(), "CAST") == 0) {
0075       lvCAST = lv;
0076     }
0077     if (lvC3EF != nullptr && lvC3HF != nullptr && lvC4EF != nullptr && lvC4HF != nullptr && lvCAST != nullptr) {
0078       break;
0079     }
0080   }
0081   edm::LogVerbatim("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
0082                                  << lvC3EF << " for C3EF; " << lvC3HF << " for C3HF; " << lvC4EF << " for C4EF; "
0083                                  << lvC4HF << " for C4HF; " << lvCAST << " for CAST. ";
0084 }
0085 
0086 //=============================================================================================
0087 
0088 CastorSD::~CastorSD() { delete showerLibrary; }
0089 
0090 //=============================================================================================
0091 
0092 double CastorSD::getEnergyDeposit(const G4Step* aStep) {
0093   double NCherPhot = 0.;
0094 
0095   // Get theTrack
0096   auto const theTrack = aStep->GetTrack();
0097 
0098   // preStepPoint information *********************************************
0099   auto const preStepPoint = aStep->GetPreStepPoint();
0100   auto const currentPV = preStepPoint->GetPhysicalVolume();
0101   auto const currentLV = currentPV->GetLogicalVolume();
0102 
0103 #ifdef EDM_ML_DEBUG
0104   edm::LogVerbatim("ForwardSim") << "CastorSD::getEnergyDeposit:"
0105                                  << "\n CurrentStepNumber , TrackID , ParentID, Particle , VertexPosition ,"
0106                                  << " LogicalVolumeAtVertex , PV, Time"
0107                                  << "\n  TRACKINFO: " << theTrack->GetCurrentStepNumber() << " , "
0108                                  << theTrack->GetTrackID() << " , " << theTrack->GetParentID() << " , "
0109                                  << theTrack->GetDefinition()->GetParticleName() << " , "
0110                                  << theTrack->GetVertexPosition() << " , "
0111                                  << theTrack->GetLogicalVolumeAtVertex()->GetName() << " , " << currentPV->GetName()
0112                                  << " , " << theTrack->GetGlobalTime();
0113 #endif
0114 
0115   // if particle moves from interaction point or "backwards (halo)
0116   const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
0117   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0118   double zint = hitPoint.z();
0119 
0120   // Check if theTrack is a muon (if so, DO NOT use Shower Library)
0121   G4int parCode = theTrack->GetDefinition()->GetPDGEncoding();
0122   bool isHad = G4TrackToParticleID::isStableHadronIon(theTrack);
0123 
0124   double meanNCherPhot = 0;
0125   G4double charge = preStepPoint->GetCharge();
0126   // VI: no Cerenkov light from neutrals
0127   if (0.0 == charge) {
0128     return meanNCherPhot;
0129   }
0130 
0131   G4double beta = preStepPoint->GetBeta();
0132   const double bThreshold = 0.67;
0133   // VI: no Cerenkov light from non-relativistic particles
0134   if (beta < bThreshold) {
0135     return meanNCherPhot;
0136   }
0137 
0138   // remember primary particle hitting the CASTOR detector
0139   TrackInformationExtractor TIextractor;
0140   TrackInformation& trkInfo = TIextractor(theTrack);
0141   if (!trkInfo.hasCastorHit()) {
0142     trkInfo.setCastorHitPID(parCode);
0143   }
0144   G4double stepl = aStep->GetStepLength() / CLHEP::cm;
0145 
0146   //int castorHitPID = trkInfo.hasCastorHit() ? std::abs(trkInfo.getCastorHitPID())
0147   //  : std::abs(parCode);
0148 
0149   // *************************************************
0150   // take into account light collection curve for plate
0151   //      double weight = curve_Castor(nameVolume, preStepPoint);
0152   //      double edep   = aStep->GetTotalEnergyDeposit() * weight;
0153   // *************************************************
0154 
0155   // *************************************************
0156   /*    comments for sensitive volumes:      
0157         C001 ...-... CP06,CBU1 ...-...CALM --- > fibres and bundle 
0158         for first release of CASTOR
0159         CASF  --- > quartz plate  for first and second releases of CASTOR  
0160         GF2Q, GFNQ, GR2Q, GRNQ 
0161         for tests with my own test geometry of HF (on ask of Gavrilov)
0162         C3TF, C4TF - for third release of CASTOR
0163   */
0164 #ifdef EDM_ML_DEBUG
0165   edm::LogVerbatim("ForwardSim") << "CastorSD::getEnergyDeposit: for ID=" << theTrack->GetTrackID()
0166                                  << " LV: " << currentLV->GetName() << " isHad:" << isHad << " pdg=" << parCode
0167                                  << " castorPID=" << trkInfo.getCastorHitPID() << " sl= " << stepl
0168                                  << " Edep= " << aStep->GetTotalEnergyDeposit();
0169 #endif
0170   if (currentLV == lvC3EF || currentLV == lvC4EF || currentLV == lvC3HF || currentLV == lvC4HF) {
0171     const double nMedium = 1.4925;
0172     //     double photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
0173     //     double photEnSpectrDL = 10714.285714;
0174 
0175     const double photEnSpectrDE = 1.24; /* see below   */
0176     /*     E = 2pi*(1./137.)*(eV*cm/370.)/lambda  =     */
0177     /*       = 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     /* default for Castor nameVolume  == "CASF" or (C3TF & C4TF)  */
0183 
0184     const double thFullRefl = 23.; /* 23.dergee */
0185     const double thFullReflRad = thFullRefl * pi / 180.;
0186 
0187     /* default for Castor nameVolume  == "CASF" or (C3TF & C4TF)  */
0188     const double thFibDir = 45.; /* .dergee */
0189     /* for test HF geometry volumes:   
0190        if(nameVolume == "GF2Q" || nameVolume == "GFNQ" ||
0191        nameVolume == "GR2Q" || nameVolume == "GRNQ")
0192        thFibDir = 0.0; // .dergee
0193     */
0194     const double thFibDirRad = thFibDir * pi / 180.;
0195 
0196     // theta of charged particle in LabRF(hit momentum direction):
0197     double costh =
0198         hit_mom.z() / sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
0199     if (zint < 0)
0200       costh = -costh;
0201     double th = acos(std::min(std::max(costh, double(-1.)), double(1.)));
0202 
0203     // just in case (can do bot use):
0204     if (th < 0.)
0205       th += twopi;
0206 
0207     // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
0208     double costhcher = 1. / (nMedium * beta);
0209     double thcher = acos(std::min(std::max(costhcher, double(-1.)), double(1.)));
0210 
0211     // diff thetas of charged part. and quartz direction in LabRF:
0212     double DelFibPart = std::abs(th - thFibDirRad);
0213 
0214     // define real distances:
0215     double d = std::abs(tan(th) - tan(thFibDirRad));
0216 
0217     double a = tan(thFibDirRad) + tan(std::abs(thFibDirRad - thFullReflRad));
0218     double r = tan(th) + tan(std::abs(th - thcher));
0219 
0220     // define losses d_qz in cone of full reflection inside quartz direction
0221     double d_qz;
0222 #ifdef EDM_ML_DEBUG
0223     int variant(0);
0224 #endif
0225     if (DelFibPart > (thFullReflRad + thcher)) {
0226       d_qz = 0.;
0227     } else {
0228       if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
0229         d_qz = 1.;
0230 #ifdef EDM_ML_DEBUG
0231         variant = 1;
0232 #endif
0233       } else {
0234         if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
0235           d_qz = 0.;
0236 #ifdef EDM_ML_DEBUG
0237           variant = 2;
0238 #endif
0239         } else {
0240 #ifdef EDM_ML_DEBUG
0241           variant = 3;
0242 #endif
0243           // use crossed length of circles(cone projection)
0244           // dC1/dC2 :
0245           double arg_arcos = 0.;
0246           double tan_arcos = 2. * a * d;
0247           if (tan_arcos != 0.)
0248             arg_arcos = (r * r - a * a - d * d) / tan_arcos;
0249           arg_arcos = std::abs(arg_arcos);
0250           double th_arcos = acos(std::min(std::max(arg_arcos, -1.), 1.));
0251           d_qz = std::abs(th_arcos / CLHEP::twopi);
0252         }
0253       }
0254     }
0255 #ifdef EDM_ML_DEBUG
0256     edm::LogVerbatim("ForwardSim") << "variant" << variant;
0257 #endif
0258 
0259     // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
0260 
0261     meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta * beta)) * photEnSpectrDE * stepl;
0262 
0263     double scale = isHad ? non_compensation_factor : 1.0;
0264 
0265     int poissNCherPhot = std::max(0, (int)G4Poisson(meanNCherPhot * scale));
0266 
0267     const double effPMTandTransport = 0.19;
0268     const double ReflPower = 0.1;
0269     double proba = d_qz + (1 - d_qz) * ReflPower;
0270     NCherPhot = poissNCherPhot * effPMTandTransport * proba * 0.307;
0271 #ifdef EDM_ML_DEBUG
0272     edm::LogVerbatim("ForwardSim") << " Nph= " << NCherPhot << " Np= " << poissNCherPhot
0273                                    << " eff= " << effPMTandTransport << " pb= " << proba << " Nmean= " << meanNCherPhot
0274                                    << " q=" << charge << " beta=" << beta << " nMedium= " << nMedium << " sl= " << stepl
0275                                    << " Nde=" << photEnSpectrDE;
0276 #endif
0277   }
0278   return NCherPhot;
0279 }
0280 
0281 //=======================================================================================
0282 
0283 uint32_t CastorSD::setDetUnitId(const G4Step* aStep) {
0284   return (numberingScheme == nullptr ? 0 : numberingScheme->getUnitID(aStep));
0285 }
0286 
0287 //=======================================================================================
0288 
0289 void CastorSD::setNumberingScheme(CastorNumberingScheme* scheme) {
0290   if (scheme != nullptr) {
0291     edm::LogVerbatim("ForwardSim") << "CastorSD: updates numbering scheme for " << GetName();
0292     delete numberingScheme;
0293     numberingScheme = scheme;
0294   }
0295 }
0296 
0297 //=======================================================================================
0298 
0299 uint32_t CastorSD::rotateUnitID(uint32_t unitID, const G4Track* track, const CastorShowerEvent& shower) {
0300   // ==============================================================
0301   //
0302   //   o   Exploit Castor phi symmetry to return newUnitID for
0303   //       shower hits based on track phi coordinate
0304   //
0305   // ==============================================================
0306 
0307   // Get 'track' phi:
0308   double trackPhi = track->GetPosition().phi();
0309   if (trackPhi < 0)
0310     trackPhi += 2 * M_PI;
0311   // Get phi from primary that gave rise to SL 'shower':
0312   double showerPhi = shower.getPrimPhi();
0313   if (showerPhi < 0)
0314     showerPhi += 2 * M_PI;
0315   // Delta phi:
0316 
0317   //  Find the OctSector for which 'track' and 'shower' belong
0318   int trackOctSector = (int)(trackPhi / (M_PI / 4));
0319   int showerOctSector = (int)(showerPhi / (M_PI / 4));
0320 
0321   uint32_t newUnitID;
0322   uint32_t sec = ((unitID >> 4) & 0xF);
0323   uint32_t complement = (unitID & 0xFFFFFF0F);
0324 
0325   // Get 'track' z:
0326   double trackZ = track->GetPosition().z();
0327 
0328   int aux;
0329   int dSec = 2 * (trackOctSector - showerOctSector);
0330   if (trackZ > 0)  // Good for revision 1.9 of CastorNumberingScheme
0331   {
0332     int sec1 = sec - dSec;
0333     //    sec -= dSec ;
0334     if (sec1 < 0)
0335       sec1 += 16;
0336     if (sec1 > 15)
0337       sec1 -= 16;
0338     sec = (uint32_t)(sec1);
0339   } else {
0340     if (dSec < 0)
0341       sec += 16;
0342     sec += dSec;
0343     aux = (int)(sec / 16);
0344     sec -= aux * 16;
0345   }
0346   sec = sec << 4;
0347   newUnitID = complement | sec;
0348 
0349 #ifdef EDM_ML_DEBUG
0350   if (newUnitID != unitID) {
0351     edm::LogVerbatim("ForwardSim") << "\n CastorSD::rotateUnitID:  "
0352                                    << "\n     unitID = " << unitID << "\n  newUnitID = " << newUnitID;
0353   }
0354 #endif
0355 
0356   return newUnitID;
0357 }
0358 
0359 //=======================================================================================
0360 
0361 bool CastorSD::getFromLibrary(const G4Step* aStep) {
0362   /////////////////////////////////////////////////////////////////////
0363   //
0364   //   Method to get hits from the Shower Library
0365   //
0366   //   "updateHit" save the Hits to a CaloG4Hit container
0367   //
0368   /////////////////////////////////////////////////////////////////////
0369 
0370   auto const theTrack = aStep->GetTrack();
0371   auto parCode = theTrack->GetDefinition()->GetPDGEncoding();
0372 
0373   // remember primary particle hitting the CASTOR detector
0374   TrackInformationExtractor TIextractor;
0375   TrackInformation& trkInfo = TIextractor(theTrack);
0376   if (!trkInfo.hasCastorHit()) {
0377     trkInfo.setCastorHitPID(parCode);
0378   }
0379 
0380   if (!useShowerLibrary) {
0381     return false;
0382   }
0383 
0384   // preStepPoint information *********************************************
0385 
0386   auto const preStepPoint = aStep->GetPreStepPoint();
0387   auto const currentPV = preStepPoint->GetPhysicalVolume();
0388   auto const currentLV = currentPV->GetLogicalVolume();
0389 
0390 #ifdef EDM_ML_DEBUG
0391   edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: for ID=" << theTrack->GetTrackID()
0392                                  << " parentID= " << theTrack->GetParentID() << " "
0393                                  << theTrack->GetDefinition()->GetParticleName() << " LV: " << currentLV->GetName()
0394                                  << " PV: " << currentPV->GetName() << "\n eta= " << theTrack->GetPosition().eta()
0395                                  << " phi= " << theTrack->GetPosition().phi()
0396                                  << " z(cm)= " << theTrack->GetPosition().z() / CLHEP::cm
0397                                  << " time(ns)= " << theTrack->GetGlobalTime()
0398                                  << " E(GeV)= " << theTrack->GetTotalEnergy() / CLHEP::GeV;
0399 
0400 #endif
0401 
0402   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0403   const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
0404   double zint = hitPoint.z();
0405   double pz = hit_mom.z();
0406 
0407   // Check if theTrack moves backward
0408   bool backward = (pz * zint < 0.) ? true : false;
0409 
0410   // Check that theTrack is above the energy threshold to use Shower Library
0411   bool aboveThreshold = (theTrack->GetKineticEnergy() > energyThresholdSL) ? true : false;
0412 
0413   // Check if theTrack is a muon (if so, DO NOT use Shower Library)
0414   bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
0415   bool isHad = G4TrackToParticleID::isStableHadronIon(theTrack);
0416 
0417   // angle condition
0418   double theta_max = M_PI - 3.1305;  // angle in radians corresponding to -5.2 eta
0419   double R_mom = sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y());
0420   double theta = atan2(R_mom, std::abs(pz));
0421   bool angleok = (theta < theta_max) ? true : false;
0422 
0423   // OkToUse
0424   double R = sqrt(hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y());
0425   bool dot = (zint < -14450. && R < 45.) ? true : false;
0426   bool inRange = (zint < -14700. || R > 193.) ? false : true;
0427   //bool OkToUse = (inRange && !dot) ? true : false;
0428 
0429   bool particleWithinShowerLibrary =
0430       aboveThreshold && (isEM || isHad) && (!backward) && inRange && !dot && angleok && currentLV == lvCAST;
0431 
0432 #ifdef EDM_ML_DEBUG
0433   edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary: ID= " << theTrack->GetTrackID() << " E>E0 "
0434                                  << aboveThreshold << " isEM " << isEM << " isHad " << isHad << " backword " << backward
0435                                  << " Ok " << (inRange && !dot) << " angle " << angleok
0436                                  << " LV: " << currentLV->GetName() << "  " << (currentLV == lvCAST) << " "
0437                                  << particleWithinShowerLibrary << " Edep= " << aStep->GetTotalEnergyDeposit();
0438 #endif
0439 
0440   // Use Castor shower library if energy is above threshold, is not a muon
0441   // and is not moving backward
0442   if (!particleWithinShowerLibrary) {
0443     return false;
0444   }
0445 
0446   // ****    Call method to retrieve hits from the ShowerLibrary   ****
0447   // always kill primary
0448   bool isKilled(true);
0449   CastorShowerEvent hits = showerLibrary->getShowerHits(aStep, isKilled);
0450 
0451   int primaryID = setTrackID(aStep);
0452 
0453   // Reset entry point for new primary
0454   resetForNewPrimary(aStep);
0455 
0456 #ifdef EDM_ML_DEBUG
0457   edm::LogVerbatim("ForwardSim") << "CastorSD::getFromLibrary:  " << hits.getNhit() << " hits for " << GetName()
0458                                  << " from " << theTrack->GetDefinition()->GetParticleName() << " of "
0459                                  << preStepPoint->GetKineticEnergy() / CLHEP::GeV << " GeV and trackID "
0460                                  << theTrack->GetTrackID() << " isHad: " << isHad;
0461 #endif
0462 
0463   // Scale to correct energy
0464   double E_track = preStepPoint->GetTotalEnergy();
0465   double E_SLhit = hits.getPrimE() * CLHEP::GeV;
0466   double scale = E_track / E_SLhit;
0467 
0468   //Non compensation
0469   if (isHad) {
0470     scale *= non_compensation_factor;  // if hadronic extend the scale with the non-compensation factor
0471   }
0472   //  Loop over hits retrieved from the library
0473   for (unsigned int i = 0; i < hits.getNhit(); ++i) {
0474     // Get nPhotoElectrons and set edepositEM / edepositHAD accordingly
0475     double nPhotoElectrons = hits.getNphotons(i) * scale;
0476 
0477     if (isEM) {
0478       edepositEM = nPhotoElectrons;
0479       edepositHAD = 0.;
0480     } else {
0481       edepositEM = 0.;
0482       edepositHAD = nPhotoElectrons;
0483     }
0484 
0485     // Get hit position and time
0486     double time = hits.getTime(i);
0487 
0488     // Get hit detID
0489     unsigned int unitID = hits.getDetID(i);
0490 
0491     // Make the detID "rotation" from one sector to another taking into account the
0492     // sectors of the impinging particle (theTrack) and of the particle that produced
0493     // the 'hits' retrieved from shower library
0494     unsigned int rotatedUnitID = rotateUnitID(unitID, theTrack, hits);
0495     currentID[0].setID(rotatedUnitID, time, primaryID, 0);
0496     processHit(aStep);
0497   }  //  End of loop over hits
0498   return isKilled;
0499 }