Merge pull request #34485 from civanch/fix_fast_hits_method

Fixed fast sim sampling of ECAL hits
cms-sw · Jul 15, 2021 · 813317b · 813317b
2 parents 9286da1 + e8e72ea
commit 813317b
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 20 deletions.
diff --git a/SimG4CMS/Calo/src/CaloSD.cc b/SimG4CMS/Calo/src/CaloSD.cc
@@ -598,7 +598,7 @@ CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack) {
   storeHit(aHit);
   TrackInformation* trkInfo = cmsTrackInformation(theTrack);
 
-  bool currentlyInsideFineVolume = isItFineCalo(aStep->GetPostStepPoint()->GetTouchable());
+  bool currentlyInsideFineVolume = !doFineCalo_ ? false : isItFineCalo(aStep->GetPostStepPoint()->GetTouchable());
 
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("DoFineCalo") << "Creating new hit " << aHit->getUnitID() << " using "
@@ -617,7 +617,7 @@ CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack) {
 
   // If fine calo is activated for the current volume, perform track/hit
   // saving logic for fineCalo
-  if (doFineCalo_ && currentlyInsideFineVolume) {
+  if (currentlyInsideFineVolume) {
     hitBookkeepingFineCalo(aStep, theTrack, aHit);
   }
   // 'Traditional', non-fine history bookkeeping

diff --git a/SimG4Core/Application/interface/LowEnergyFastSimModel.h b/SimG4Core/Application/interface/LowEnergyFastSimModel.h
@@ -9,6 +9,7 @@
 #include "GFlashHitMaker.hh"
 #include "G4Region.hh"
 #include "G4Types.hh"
+#include "G4ThreeVector.hh"
 
 class TrackingAction;
 
@@ -24,8 +25,10 @@ class LowEnergyFastSimModel : public G4VFastSimulationModel {
   G4double fEmax;
   const G4Envelope* fRegion;
   const TrackingAction* fTrackingAction;
+  G4bool fCheck;
+  G4ThreeVector fTailPos;
   GFlashHitMaker fHitMaker;
-  LowEnergyFastSimParam param;
+  LowEnergyFastSimParam fParam;
 };
 
 #endif
diff --git a/SimG4Core/Application/src/LowEnergyFastSimModel.cc b/SimG4Core/Application/src/LowEnergyFastSimModel.cc
@@ -11,10 +11,14 @@
 #include "G4Region.hh"
 #include "G4PhysicalConstants.hh"
 
-constexpr double twomass = 2 * CLHEP::electron_mass_c2;
+constexpr G4double twomass = 2 * CLHEP::electron_mass_c2;
 
 LowEnergyFastSimModel::LowEnergyFastSimModel(const G4String& name, G4Region* region, const edm::ParameterSet& parSet)
-    : G4VFastSimulationModel(name, region), fRegion(region), fTrackingAction(nullptr) {
+    : G4VFastSimulationModel(name, region),
+      fRegion(region),
+      fTrackingAction(nullptr),
+      fCheck(false),
+      fTailPos(0., 0., 0.) {
   fEmax = parSet.getParameter<double>("LowEnergyGflashEcalEmax") * CLHEP::GeV;
 }
 
@@ -24,32 +28,32 @@ G4bool LowEnergyFastSimModel::IsApplicable(const G4ParticleDefinition& particle)
 
 G4bool LowEnergyFastSimModel::ModelTrigger(const G4FastTrack& fastTrack) {
   const G4Track* track = fastTrack.GetPrimaryTrack();
-  if (nullptr == fTrackingAction) {
-    fTrackingAction = static_cast<const TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
+  if (fCheck) {
+    if (nullptr == fTrackingAction) {
+      fTrackingAction = static_cast<const TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
+    }
+    int pdgMother = std::abs(fTrackingAction->geant4Track()->GetDefinition()->GetPDGEncoding());
+    if (pdgMother == 11 || pdgMother == 22)
+      return false;
   }
-  int pdgMother = std::abs(fTrackingAction->geant4Track()->GetDefinition()->GetPDGEncoding());
-  if (pdgMother == 11 || pdgMother == 22)
-    return false;
   G4double energy = track->GetKineticEnergy();
-  return energy < fEmax && fRegion == fastTrack.GetEnvelope();
+  return (energy < fEmax && fRegion == fastTrack.GetEnvelope());
 }
 
 void LowEnergyFastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
   fastStep.KillPrimaryTrack();
   fastStep.SetPrimaryTrackPathLength(0.0);
-  double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
+  G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
 
   const G4ThreeVector& pos = fastTrack.GetPrimaryTrack()->GetPosition();
 
-  G4double inPointEnergy = param.GetInPointEnergyFraction(energy) * energy;
+  G4double inPointEnergy = fParam.GetInPointEnergyFraction(energy) * energy;
 
   // take into account positron annihilation (not included in in-point)
   if (-11 == fastTrack.GetPrimaryTrack()->GetDefinition()->GetPDGEncoding())
     energy += twomass;
 
   const G4ThreeVector& momDir = fastTrack.GetPrimaryTrack()->GetMomentumDirection();
-  const G4ThreeVector& ortho = momDir.orthogonal();
-  const G4ThreeVector& cross = momDir.cross(ortho);
 
   // in point energy deposition
   GFlashEnergySpot spot;
@@ -59,17 +63,19 @@ void LowEnergyFastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastS
 
   // tail energy deposition
   G4double etail = energy - inPointEnergy;
-  const G4int nspots = int(etail) + 1;
+  const G4int nspots = (G4int)(etail) + 1;
   const G4double tailEnergy = etail / (G4double)nspots;
   for (G4int i = 0; i < nspots; ++i) {
-    const G4double radius = param.GetRadius(energy);
-    const G4double z = param.GetZ();
+    const G4double r = fParam.GetRadius(energy);
+    const G4double z = fParam.GetZ();
 
     const G4double phi = CLHEP::twopi * G4UniformRand();
-    const G4ThreeVector tailPos = pos + z * momDir + radius * std::cos(phi) * ortho + radius * std::sin(phi) * cross;
+    fTailPos.set(r * std::cos(phi), r * std::sin(phi), z);
+    fTailPos.rotateUz(momDir);
+    fTailPos += pos;
 
     spot.SetEnergy(tailEnergy);
-    spot.SetPosition(tailPos);
+    spot.SetPosition(fTailPos);
     fHitMaker.make(&spot, &fastTrack);
   }
 }