Merge pull request #22145 from civanch/g4particles_Showers

Geant4 Sensitive Detectors update 4: revised selection of particle type for shower libraries

SimG4CMS/Calo/interface/CaloSD.h

@@ -22,9 +22,6 @@
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
-// To be replaced by something else 
-/* #include "Utilities/Notification/interface/TimerProxy.h" */
-
 #include "G4VPhysicalVolume.hh"
 #include "G4Track.hh"
 #include "G4VGFlashSensitiveDetector.hh"
@@ -134,7 +131,6 @@ class CaloSD : public SensitiveCaloDetector,
   double                          correctT;
   double                          kmaxIon, kmaxNeutron, kmaxProton;
 
-  G4int                           emPDG, epPDG, gammaPDG;
   bool                            forceSave;
 
 private:

SimG4CMS/Calo/interface/HCalSD.h

@@ -97,9 +97,10 @@ class HCalSD : public CaloSD, public Observer<const BeginOfJob *> {
   bool                          testNumber, neutralDensity, testNS_;
   double                        birk1, birk2, birk3, betaThr;
   bool                          useHF, useShowerLibrary, useParam, applyFidCut;
+  bool                          isEM;
   double                        eminHitHB, eminHitHE, eminHitHO, eminHitHF;
   double                        deliveredLumi;
-  G4int                         mumPDG, mupPDG, depth_;
+  G4int                         depth_;
   std::vector<double>           gpar;
   std::vector<int>              hfLevels;
   std::vector<G4String>         hfNames, fibreNames, matNames;

SimG4CMS/Calo/interface/HFShowerLibrary.h

@@ -76,10 +76,6 @@ class HFShowerLibrary {
   double              dphi, rMin, rMax;
   std::vector<double> gpar;
 
-  int                 emPDG, epPDG, gammaPDG;
-  int                 pi0PDG, etaPDG, nuePDG, numuPDG, nutauPDG;
-  int                 anuePDG, anumuPDG, anutauPDG, geantinoPDG;
-
   int                 npe;
   HFShowerPhotonCollection pe;
   HFShowerPhotonCollection* photo;

SimG4CMS/Calo/src/CaloSD.cc

@@ -6,6 +6,7 @@
 #include "SimG4CMS/Calo/interface/CaloSD.h"
 #include "SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h"
 #include "SimG4Core/Notification/interface/TrackInformation.h"
+#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
 #include "SimG4Core/Application/interface/EventAction.h"
 
 #include "G4EventManager.hh"
@@ -109,11 +110,8 @@ bool CaloSD::ProcessHits(G4GFlashSpot* aSpot, G4TouchableHistory*) {
 
   if (aSpot != nullptr) {   
     theTrack = const_cast<G4Track *>(aSpot->GetOriginatorTrack()->GetPrimaryTrack());
-    G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
 
-    if (particleCode == emPDG ||
-        particleCode == epPDG ||
-        particleCode == gammaPDG ) {
+    if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
       edepositEM  = aSpot->GetEnergySpot()->GetEnergy();
       edepositHAD = 0.;
     } else {
@@ -253,10 +251,7 @@ bool CaloSD::getStepInfo(G4Step* aStep) {
 #endif
   }
 
-  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
-  if (particleCode == emPDG ||
-      particleCode == epPDG ||
-      particleCode == gammaPDG ) {
+  if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
     edepositEM  = getEnergyDeposit(aStep);
     edepositHAD = 0.;
   } else {
@@ -459,15 +454,6 @@ double CaloSD::getAttenuation(const G4Step* aStep, double birk1, double birk2, d
 }
 
 void CaloSD::update(const BeginOfRun *) {
-  G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
-  G4String particleName;
-  emPDG = theParticleTable->FindParticle(particleName="e-")->GetPDGEncoding();
-  epPDG = theParticleTable->FindParticle(particleName="e+")->GetPDGEncoding();
-  gammaPDG = theParticleTable->FindParticle(particleName="gamma")->GetPDGEncoding();
-#ifdef DebugLog
-  edm::LogInfo("CaloSim") << "CaloSD: Particle code for e- = " << emPDG
-			  << " for e+ = " << epPDG << " for gamma = " << gammaPDG;
-#endif
   initRun();
   runInit = true;
 } 

SimG4CMS/Calo/src/HCalSD.cc

@@ -8,6 +8,7 @@
 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
 #include "SimG4CMS/Calo/interface/HcalTestNS.h"
 #include "SimG4Core/Notification/interface/TrackInformation.h"
+#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
 #include "DetectorDescription/Core/interface/DDFilter.h"
 #include "DetectorDescription/Core/interface/DDFilteredView.h"
@@ -296,7 +297,7 @@ HCalSD::HCalSD(const std::string& name, const DDCompactView & cpv,
     edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << matNames[i]
                             << " pointer " << materials[i];
 
-  mumPDG = mupPDG = 0;
+  isEM = false;
 
   if (useLayerWt) readWeightFromFile(file);
 
@@ -365,20 +366,20 @@ bool HCalSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
     depth_ = (aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(0))%10;
     const G4LogicalVolume* lv =
       aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
+    isEM = (G4TrackToParticleID::isGammaElectronPositron(aStep->GetTrack())) ? true : false;
     G4String nameVolume = lv->GetName();
     if (isItHF(aStep)) {
-      G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
       double weight(1.0);
       if (m_HFDarkening) {
 	G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
 	double r = hitPoint.perp()/CLHEP::cm;
 	double z = std::abs(hitPoint.z())/CLHEP::cm;
 	double dose_acquired = 0.;
-  if (z>=HFDarkening::lowZLimit && z <= HFDarkening::upperZLimit) {
-    unsigned int hfZLayer = (int)((z - HFDarkening::lowZLimit)/5);
-    if (hfZLayer >= HFDarkening::upperZLimit) hfZLayer = (HFDarkening::upperZLimit-1);
+	if (z>=HFDarkening::lowZLimit && z <= HFDarkening::upperZLimit) {
+	  unsigned int hfZLayer = (int)((z - HFDarkening::lowZLimit)/5);
+	  if (hfZLayer >= HFDarkening::upperZLimit) hfZLayer = (HFDarkening::upperZLimit-1);
 	  float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
-    for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
+	  for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
 	    dose_acquired = m_HFDarkening->dose(i,r);
 	    weight *= m_HFDarkening->degradation(normalized_lumi*dose_acquired);
 	  }
@@ -403,16 +404,18 @@ bool HCalSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
 			    << " hits afterParamS*";
 #endif 
       } else {
-        bool notaMuon = true;
-        if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
-        if (useShowerLibrary && notaMuon) {
+	// shower library is applied only for gamma, e+- or stable hadrons
+	// muons are tracked via HF, hyperons does not interact
+        if (useShowerLibrary && !G4TrackToParticleID::isMuon(aStep->GetTrack())) {
+	  if(isEM || G4TrackToParticleID::isStableHadronIon(aStep->GetTrack())) {
+	    getFromLibrary(aStep, weight);
+	  }
 #ifdef EDM_ML_DEBUG
           LogDebug("HcalSim") << "HCalSD: Starts shower library from " 
                               << nameVolume << " for Track " 
                               << aStep->GetTrack()->GetTrackID() << " ("
                               << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
 #endif
-          getFromLibrary(aStep, weight);
         } else if (isItFibre(lv)) {
 #ifdef EDM_ML_DEBUG
           LogDebug("HcalSim") << "HCalSD: Hit at Fibre in " << nameVolume 
@@ -608,13 +611,6 @@ void HCalSD::update(const BeginOfJob * job) {
 
 void HCalSD::initRun() {
   G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
-  G4String          particleName;
-  mumPDG = theParticleTable->FindParticle(particleName="mu-")->GetPDGEncoding();
-  mupPDG = theParticleTable->FindParticle(particleName="mu+")->GetPDGEncoding();
-#ifdef EDM_ML_DEBUG
-  LogDebug("HcalSim") << "HCalSD: Particle code for mu- = " << mumPDG
-		      << " for mu+ = " << mupPDG;
-#endif
   if (showerLibrary) showerLibrary->initRun(theParticleTable,hcalConstants);
   if (showerParam)   showerParam->initRun(theParticleTable,hcalConstants);
   if (hfshower)      hfshower->initRun(theParticleTable,hcalConstants);
@@ -785,52 +781,50 @@ bool HCalSD::isItinFidVolume (const G4ThreeVector& hitPoint) {
 void HCalSD::getFromLibrary (G4Step* aStep, double weight) {
   preStepPoint  = aStep->GetPreStepPoint(); 
   theTrack = aStep->GetTrack();   
-  int det       = 5;
-  bool ok;
+  int det  = 5;
+  bool ok(false);
 
   std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, ok, weight, false);
 
-  double etrack    = preStepPoint->GetKineticEnergy();
-  int    primaryID = setTrackID(aStep);
+  if(!hits.empty()) {
+    double etrack    = preStepPoint->GetKineticEnergy();
+    int    primaryID = setTrackID(aStep);
 
-  // Reset entry point for new primary
-  posGlobal = preStepPoint->GetPosition();
-  resetForNewPrimary(posGlobal, etrack);
+    // Reset entry point for new primary
+    posGlobal = preStepPoint->GetPosition();
+    resetForNewPrimary(posGlobal, etrack);
 
-  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
-  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
-    edepositEM  = 1.*GeV;
-    edepositHAD = 0.;
-  } else {
-    edepositEM  = 0.;
-    edepositHAD = 1.*GeV;
-  }
+    if (isEM) {
+      edepositEM  = 1.*GeV;
+      edepositHAD = 0.;
+    } else {
+      edepositEM  = 0.;
+      edepositHAD = 1.*GeV;
+    }
 #ifdef EDM_ML_DEBUG
-  edm::LogInfo("HcalSim") << "HCalSD::getFromLibrary " <<hits.size() 
-                          << " hits for " << GetName() << " of " << primaryID 
-                          << " with " << theTrack->GetDefinition()->GetParticleName() 
-                          << " of " << preStepPoint->GetKineticEnergy()/GeV << " GeV";
+    edm::LogInfo("HcalSim") << "HCalSD::getFromLibrary " <<hits.size() 
+			    << " hits for " << GetName() << " of " << primaryID 
+			    << " with " << theTrack->GetDefinition()->GetParticleName() 
+			    << " of " << preStepPoint->GetKineticEnergy()/GeV << " GeV";
 #endif
-  for (unsigned int i=0; i<hits.size(); ++i) {
-    G4ThreeVector hitPoint = hits[i].position;
-    if (isItinFidVolume (hitPoint)) {
-      int depth              = hits[i].depth;
-      double time            = hits[i].time;
-      unsigned int unitID    = setDetUnitId(det, hitPoint, depth);
-      currentID.setID(unitID, time, primaryID, 0);
+    for (unsigned int i=0; i<hits.size(); ++i) {
+      G4ThreeVector hitPoint = hits[i].position;
+      if (isItinFidVolume (hitPoint)) {
+	int depth              = hits[i].depth;
+	double time            = hits[i].time;
+	unsigned int unitID    = setDetUnitId(det, hitPoint, depth);
+	currentID.setID(unitID, time, primaryID, 0);
 #ifdef plotDebug
-      plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
-      bool emType = false;
-      if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
-	emType = true;
-      plotHF(hitPoint,emType);
+	plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
+	plotHF(hitPoint,isEM);
 #endif
 
-      // check if it is in the same unit and timeslice as the previous one
-      if (currentID == previousID) {
-	updateHit(currentHit);
-      } else {
-	if (!checkHit()) currentHit = createNewHit();
+	// check if it is in the same unit and timeslice as the previous one
+	if (currentID == previousID) {
+	  updateHit(currentHit);
+	} else {
+	  if (!checkHit()) currentHit = createNewHit();
+	}
       }
     }
   }
@@ -848,14 +842,12 @@ void HCalSD::getFromLibrary (G4Step* aStep, double weight) {
 void HCalSD::hitForFibre (const G4Step* aStep, double weight) { // if not ParamShower
 
   const G4StepPoint* preStepPoint  = aStep->GetPreStepPoint();
-  const G4Track*     theTrack      = aStep->GetTrack();
   int primaryID = setTrackID(aStep);
 
   int det   = 5;
   std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight);
 
-  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
-  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
+  if (isEM) {
     edepositEM  = 1.*GeV;
     edepositHAD = 0.;
   } else {
@@ -880,10 +872,7 @@ void HCalSD::hitForFibre (const G4Step* aStep, double weight) { // if not ParamS
 	currentID.setID(unitID, time, primaryID, 0);
 #ifdef plotDebug
 	plotProfile(aStep, hitPoint, edepositEM, time, depth);
-	bool emType = false;
-	if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
-	  emType = true;
-	plotHF(hitPoint,emType);
+	plotHF(hitPoint,isEM);
 #endif
 	// check if it is in the same unit and timeslice as the previous one
 	if (currentID == previousID) {

SimG4CMS/Calo/src/HFShowerLibrary.cc

@@ -8,6 +8,7 @@
 #include "DetectorDescription/Core/interface/DDFilter.h"
 #include "DetectorDescription/Core/interface/DDFilteredView.h"
 #include "DetectorDescription/Core/interface/DDValue.h"
+#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
 
 #include "FWCore/Utilities/interface/Exception.h"
 
@@ -119,9 +120,6 @@ HFShowerLibrary::HFShowerLibrary(const std::string & name, const DDCompactView &
 
   fibre = new HFFibre(name, cpv, p);
   photo = new HFShowerPhotonCollection;
-  emPDG = epPDG = gammaPDG = 0;
-  pi0PDG = etaPDG = nuePDG = numuPDG = nutauPDG= 0;
-  anuePDG= anumuPDG = anutauPDG = geantinoPDG = 0;
 }
 
 HFShowerLibrary::~HFShowerLibrary() {
@@ -135,30 +133,6 @@ void HFShowerLibrary::initRun(G4ParticleTable * theParticleTable,
 			      HcalDDDSimConstants* hcons) {
 
   if (fibre) fibre->initRun(hcons);
-
-  G4String parName;
-  emPDG = theParticleTable->FindParticle(parName="e-")->GetPDGEncoding();
-  epPDG = theParticleTable->FindParticle(parName="e+")->GetPDGEncoding();
-  gammaPDG = theParticleTable->FindParticle(parName="gamma")->GetPDGEncoding();
-  pi0PDG = theParticleTable->FindParticle(parName="pi0")->GetPDGEncoding();
-  etaPDG = theParticleTable->FindParticle(parName="eta")->GetPDGEncoding();
-  nuePDG = theParticleTable->FindParticle(parName="nu_e")->GetPDGEncoding();
-  numuPDG = theParticleTable->FindParticle(parName="nu_mu")->GetPDGEncoding();
-  nutauPDG= theParticleTable->FindParticle(parName="nu_tau")->GetPDGEncoding();
-  anuePDG = theParticleTable->FindParticle(parName="anti_nu_e")->GetPDGEncoding();
-  anumuPDG= theParticleTable->FindParticle(parName="anti_nu_mu")->GetPDGEncoding();
-  anutauPDG= theParticleTable->FindParticle(parName="anti_nu_tau")->GetPDGEncoding();
-  geantinoPDG= theParticleTable->FindParticle(parName="geantino")->GetPDGEncoding();
-#ifdef DebugLog
-  edm::LogInfo("HFShower") << "HFShowerLibrary: Particle codes for e- = " 
-			   << emPDG << ", e+ = " << epPDG << ", gamma = " 
-			   << gammaPDG << ", pi0 = " << pi0PDG << ", eta = " 
-			   << etaPDG << ", geantino = " << geantinoPDG 
-			   << "\n        nu_e = " << nuePDG << ", nu_mu = " 
-			   << numuPDG << ", nu_tau = " << nutauPDG 
-			   << ", anti_nu_e = " << anuePDG << ", anti_nu_mu = " 
-			   << anumuPDG << ", anti_nu_tau = " << anutauPDG;
-#endif
 
   //Radius (minimum and maximum)
   std::vector<double> rTable = hcons->getRTableHF();
@@ -198,6 +172,11 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(G4Step * aStep,
   G4String      partType = track->GetDefinition()->GetParticleName();
   int           parCode  = track->GetDefinition()->GetPDGEncoding();
 
+  // VI: for ions use internally pdg code of alpha in order to keep 
+  // consistency with previous simulation
+  if(track->GetDefinition()->IsGeneralIon()) { parCode = 1000020040; }
+
+
 #ifdef DebugLog
   G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
   double zoff   = localPos.z() + 0.5*gpar[1];
@@ -214,7 +193,10 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(G4Step * aStep,
 #endif
 
   double tSlice = (postStepPoint->GetGlobalTime())/nanosecond;
-  double pin    = preStepPoint->GetTotalEnergy();
+
+  // use kinetic energy for protons and ions
+  double pin = (track->GetDefinition()->GetBaryonNumber() > 0) 
+    ? preStepPoint->GetKineticEnergy() : preStepPoint->GetTotalEnergy();
 
   return fillHits(hitPoint,momDir,parCode,pin,ok,weight,tSlice,onlyLong);
 }
@@ -226,12 +208,17 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::fillHits(const G4ThreeVector
 
   std::vector<HFShowerLibrary::Hit> hit;
   ok = false;
-  if (parCode == pi0PDG || parCode == etaPDG || parCode == nuePDG ||
-      parCode == numuPDG || parCode == nutauPDG || parCode == anuePDG ||
-      parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG)
+  bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
+  // shower is built only for gamma, e+- and stable hadrons
+  if (!isEM && !G4TrackToParticleID::isStableHadron(parCode)) { 
     return hit;
+  }
   ok = true;
 
+  // remove low-energy component
+  const double threshold = 50*MeV;
+  if(pin < threshold) { return hit; }
+
   double pz     = momDir.z(); 
   double zint   = hitPoint.z(); 
 
@@ -244,7 +231,7 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::fillHits(const G4ThreeVector
   double ctheta = cos(momDir.theta());
   double stheta = sin(momDir.theta());
 
-  if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
+  if (isEM) {
     if (pin<pmom[nMomBin-1]) {
       interpolate(0, pin);
     } else {

SimG4CMS/Forward/interface/CastorShowerLibrary.h

@@ -70,10 +70,6 @@ class CastorShowerLibrary {
 
   std::vector<double> pmom;
 
-  int                 emPDG, epPDG, gammaPDG;
-  int                 pi0PDG, etaPDG, nuePDG, numuPDG, nutauPDG;
-  int                 anuePDG, anumuPDG, anutauPDG, geantinoPDG;
-  int                 mumPDG, mupPDG;
 // new variables (bins in eta and phi)
   unsigned int        nBinsE, nBinsEta, nBinsPhi;
   unsigned int        nEvtPerBinE, nEvtPerBinEta, nEvtPerBinPhi;