From f11ba0f7baac9d00da47c6ac405cc67e8bccb539 Mon Sep 17 00:00:00 2001
From: Rui Zhang <rui.zhang@cern.ch>
Date: Thu, 24 Oct 2024 10:38:51 +0200
Subject: [PATCH] [Simulation] EMEC simulation AF3 in G4

Add Sim.FastCalo.doEMECFCS flag and remove AbsEta and Ekin flags. They can be set via
--postExec 'sim:cfg.getPublicTool(FastSimulationMasterTool).FastSimulations[FastCaloSim].EkinMax=9999;cfg.getPublicTool(FastSimulationMasterTool).FastSimulations[FastCaloSim].EkinMin=99;'
---
 .../G4AtlasTools/src/G4CaloTransportTool.cxx  | 21 ++++---
 .../G4FastSimulation/CMakeLists.txt           |  2 +-
 .../python/G4FastSimulationConfig.py          | 20 +++++++
 .../G4FastSimulation/src/FastCaloSim.cxx      | 59 ++++++++++++++-----
 .../G4FastSimulation/src/FastCaloSim.h        | 18 ++++++
 .../G4FastSimulation/src/FastCaloSimTool.cxx  |  3 +-
 .../G4FastSimulation/src/FastCaloSimTool.h    | 11 ++++
 .../SimulationConfig/python/SimConfigFlags.py |  3 +
 .../python/SimulationMetadata.py              |  4 ++
 9 files changed, 116 insertions(+), 25 deletions(-)

diff --git a/Simulation/G4Atlas/G4AtlasTools/src/G4CaloTransportTool.cxx b/Simulation/G4Atlas/G4AtlasTools/src/G4CaloTransportTool.cxx
index 2c82131ff92..b87b0b6c62a 100644
--- a/Simulation/G4Atlas/G4AtlasTools/src/G4CaloTransportTool.cxx
+++ b/Simulation/G4Atlas/G4AtlasTools/src/G4CaloTransportTool.cxx
@@ -169,13 +169,20 @@ std::vector<G4FieldTrack> G4CaloTransportTool::transport(
     // Fill the output vector with the updated track
     outputStepVector.push_back(tmpFieldTrack);
     // Get the name of the volume in which the particle is located
-    std::string volName =
-        navigator
-            ->LocateGlobalPointAndSetup(tmpFieldTrack.GetPosition(), nullptr)
-            ->GetName();
-    // We stop the track navigation once we have reached the provided volume
-    if (volName.find(m_transportLimitVolume) != std::string::npos)
-      break;
+    auto volumn = navigator->LocateGlobalPointAndSetup(tmpFieldTrack.GetPosition(), nullptr);
+    if (volumn != nullptr) {
+      std::string volName = volumn->GetName();
+      // We stop the track navigation once we have reached the provided volume
+      if (volName.find(m_transportLimitVolume) != std::string::npos)
+        break;
+    } else {
+      G4ExceptionDescription description;
+      description <<"at step " << iStep << "/" << m_maxSteps.value() << " with position " << tmpFieldTrack.GetPosition() << " and momentum " << tmpFieldTrack.GetMomentum();
+      G4Exception("G4CaloTransportTool::transport",
+            "Cannot get LocateGlobalPointAndSetup",
+            JustWarning,
+            description);
+    }
   }
 
   return outputStepVector;
diff --git a/Simulation/G4Utilities/G4FastSimulation/CMakeLists.txt b/Simulation/G4Utilities/G4FastSimulation/CMakeLists.txt
index 591cd68c47e..9f96eda1475 100644
--- a/Simulation/G4Utilities/G4FastSimulation/CMakeLists.txt
+++ b/Simulation/G4Utilities/G4FastSimulation/CMakeLists.txt
@@ -18,7 +18,7 @@ endif()
 # FastCaloSim G4VFastSimulationModel should only be used in Athena builds
 set ( extra_src )
 if( NOT SIMULATIONBASE )
-   set( extra_src src/FastCaloSim.cxx  src/FastCaloSimTool.cxx  )
+   set( extra_src src/FastCaloSim.cxx src/FastCaloSimTool.cxx  )
 endif()
 
 # Component(s) in the package:
diff --git a/Simulation/G4Utilities/G4FastSimulation/python/G4FastSimulationConfig.py b/Simulation/G4Utilities/G4FastSimulation/python/G4FastSimulationConfig.py
index a549b0e3761..e8b0da521f2 100644
--- a/Simulation/G4Utilities/G4FastSimulation/python/G4FastSimulationConfig.py
+++ b/Simulation/G4Utilities/G4FastSimulation/python/G4FastSimulationConfig.py
@@ -40,5 +40,25 @@ def FastCaloSimCfg(flags, **kwargs):
     from G4AtlasTools.G4AtlasToolsConfig import G4CaloTransportToolCfg
     kwargs.setdefault("G4CaloTransportTool", result.addPublicTool(result.popToolsAndMerge(G4CaloTransportToolCfg(flags))))
     
+    # Config FastCaloSim
+    kwargs.setdefault('doEMECFCS', flags.Sim.FastCalo.doEMECFCS)
+
+    if flags.Sim.FastCalo.doEMECFCS:  # AF3 in EMEC and G4 in rest
+        kwargs.setdefault('doPhotons', True)
+        kwargs.setdefault('doElectrons', True)
+        kwargs.setdefault('doHadrons', False)
+        kwargs.setdefault('AbsEtaMin', 1.5)
+        kwargs.setdefault('AbsEtaMax', 3.2)
+        kwargs.setdefault('EkinMin', 0)
+        kwargs.setdefault('EkinMax', 8192)
+    else: # These are set to AF3 configuration
+        kwargs.setdefault('doPhotons', True)
+        kwargs.setdefault('doElectrons', True)
+        kwargs.setdefault('doHadrons', True)
+        kwargs.setdefault('AbsEtaMin', 0)
+        kwargs.setdefault('AbsEtaMax', 10)
+        kwargs.setdefault('EkinMin', 0)
+        kwargs.setdefault('EkinMax', float('inf'))
+
     result.setPrivateTools(CompFactory.FastCaloSimTool(name="FastCaloSim", **kwargs))
     return result
diff --git a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.cxx b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.cxx
index 6c9feecc69d..942d624d411 100644
--- a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.cxx
+++ b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.cxx
@@ -40,6 +40,14 @@ FastCaloSim::FastCaloSim(const std::string& name,
                          const ServiceHandle<ISF::IFastCaloSimParamSvc>& FastCaloSimSvc,
                          const Gaudi::Property<std::string>& CaloCellContainerSDName,
                          const Gaudi::Property<bool>& doG4Transport,
+                         const Gaudi::Property<bool>& doPhotons,
+                         const Gaudi::Property<bool>& doElectrons,
+                         const Gaudi::Property<bool>& doHadrons,
+                         const Gaudi::Property<float>& AbsEtaMin,
+                         const Gaudi::Property<float>& AbsEtaMax,
+                         const Gaudi::Property<float>& EkinMin,
+                         const Gaudi::Property<float>& EkinMax,
+                         const Gaudi::Property<bool>& doEMECFCS,
                          FastCaloSimTool * FastCaloSimTool)
 
 : G4VFastSimulationModel(name),
@@ -50,6 +58,14 @@ FastCaloSim::FastCaloSim(const std::string& name,
   m_FastCaloSimSvc(FastCaloSimSvc),
   m_CaloCellContainerSDName(CaloCellContainerSDName),
   m_doG4Transport(doG4Transport),
+  m_doPhotons(doPhotons),
+  m_doElectrons(doElectrons),
+  m_doHadrons(doHadrons),
+  m_AbsEtaMin(AbsEtaMin),
+  m_AbsEtaMax(AbsEtaMax),
+  m_EkinMin(EkinMin),
+  m_EkinMax(EkinMax),
+  m_doEMECFCS(doEMECFCS),
   m_FastCaloSimTool(FastCaloSimTool)
 {
 }
@@ -79,17 +95,17 @@ G4bool FastCaloSim::IsApplicable(const G4ParticleDefinition& particleType)
   bool isHadron   = MC::isHadron(particleType.GetPDGEncoding());
 
   // FastCaloSim is applicable if it is photon, electron, positron or any hadron
-  bool IsApplicable = isPhoton || isElectron || isPositron || isHadron;
+  bool isApplicable = (isPhoton && m_doPhotons) || (isElectron && m_doElectrons) || (isPositron && m_doElectrons) || (isHadron && m_doHadrons);
 
   #ifdef FCS_DEBUG
     const std::string pName = particleType.GetParticleName();
     G4cout<< "[FastCaloSim::IsApplicable] Got " << pName <<G4endl;
-    if(IsApplicable) G4cout<<"[FastCaloSim::IsApplicable] APPLICABLE"<<G4endl;
+    if(isApplicable) G4cout<<"[FastCaloSim::IsApplicable] APPLICABLE"<<G4endl;
     else G4cout<<"[FastCaloSim::IsApplicable] NOT APPLICABLE"<<G4endl;
   #endif
 
 
-  return IsApplicable;
+  return isApplicable;
 }
 
 G4bool FastCaloSim::ModelTrigger(const G4FastTrack& fastTrack)
@@ -109,6 +125,31 @@ G4bool FastCaloSim::ModelTrigger(const G4FastTrack& fastTrack)
                                     <<G4endl;
   #endif
 
+  // Get particle definition 
+  const G4ParticleDefinition * G4Particle = fastTrack.GetPrimaryTrack() -> GetDefinition();
+  // Get particle kinetic energy
+  const float Ekin = fastTrack.GetPrimaryTrack() -> GetKineticEnergy();
+  // Get particle position eta
+  const float eta_pos = (fastTrack.GetPrimaryTrack() -> GetPosition()).eta();
+  
+  // Check particle type
+  bool isPhoton    = G4Particle == G4Gamma::Definition();
+  bool isElectron  = G4Particle == G4Electron::Definition();
+  bool isPositron  = G4Particle == G4Positron::Definition();
+  bool isPionPlus  = G4Particle == G4PionPlus::Definition();
+  bool isPionMinus = G4Particle == G4PionMinus::Definition();
+
+
+  // Check if there is a configuration for this PID
+  bool withinEtaRange = (std::abs(eta_pos) > m_AbsEtaMin) && (std::abs(eta_pos) < m_AbsEtaMax);
+  bool withinEkinRange = (Ekin > m_EkinMin) && (Ekin < m_EkinMax);
+  
+  if (!(withinEtaRange && withinEkinRange)) {
+    #ifdef FCS_DEBUG
+      G4cout<<"[FastCaloSim::ModelTrigger] Model not triggered"<<G4endl;
+    #endif
+    return false;
+  }
 
   // Simulate particles below 50 keV with Geant4 to have same config as ISF implementation 
   if (fastTrack.GetPrimaryTrack() -> GetKineticEnergy() < 0.05) {
@@ -130,18 +171,6 @@ G4bool FastCaloSim::ModelTrigger(const G4FastTrack& fastTrack)
   float minEkinPions = 200;
   float minEkinOtherHadrons = 400;
 
-  // Get particle definition 
-  const G4ParticleDefinition * G4Particle = fastTrack.GetPrimaryTrack() -> GetDefinition();
-  // Get particle kinetic energy
-  const float Ekin = fastTrack.GetPrimaryTrack() -> GetKineticEnergy();
-  
-  // Check particle type
-  bool isPhoton    = G4Particle == G4Gamma::Definition();
-  bool isElectron  = G4Particle == G4Electron::Definition();
-  bool isPositron  = G4Particle == G4Positron::Definition();
-  bool isPionPlus  = G4Particle == G4PionPlus::Definition();
-  bool isPionMinus = G4Particle == G4PionMinus::Definition();
-
   // Pass all photons, electrons and positrons to FastCaloSim
   if (isPhoton || isElectron || isPositron){
     #ifdef FCS_DEBUG
diff --git a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.h b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.h
index 238df3a4189..63b3201136d 100644
--- a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.h
+++ b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSim.h
@@ -40,6 +40,14 @@ class FastCaloSim: public G4VFastSimulationModel
               const ServiceHandle<ISF::IFastCaloSimParamSvc>& FastCaloSimSvc,
               const Gaudi::Property<std::string>& CaloCellContainerSDName,
               const Gaudi::Property<bool>& doG4Transport,
+              const Gaudi::Property<bool>& doPhotons,
+              const Gaudi::Property<bool>& doElectrons,
+              const Gaudi::Property<bool>& doHadrons,
+              const Gaudi::Property<float>& EtaLow,
+              const Gaudi::Property<float>& EtaHigh,
+              const Gaudi::Property<float>& EkinLow,
+              const Gaudi::Property<float>& EkinHigh,
+              const Gaudi::Property<bool>& doEMECFCS,
               FastCaloSimTool * FastCaloSimTool);
   ~FastCaloSim() {}
 
@@ -81,6 +89,16 @@ class FastCaloSim: public G4VFastSimulationModel
   // Boolean flag to enable Geant4 transportation
   Gaudi::Property<bool> m_doG4Transport;
 
+  // Boundaries to enable AF3 transportation
+  Gaudi::Property<bool> m_doPhotons;
+  Gaudi::Property<bool> m_doElectrons;
+  Gaudi::Property<bool> m_doHadrons;
+  Gaudi::Property<float> m_AbsEtaMin;
+  Gaudi::Property<float> m_AbsEtaMax;
+  Gaudi::Property<float> m_EkinMin;
+  Gaudi::Property<float> m_EkinMax;
+  Gaudi::Property<float> m_doEMECFCS;
+
   // Fast simulation FastCaloSimTool 
   FastCaloSimTool * m_FastCaloSimTool;
 };
diff --git a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.cxx b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.cxx
index e99e3f198b4..138f5dabd01 100644
--- a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.cxx
+++ b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.cxx
@@ -59,7 +59,6 @@ G4VFastSimulationModel* FastCaloSimTool::makeFastSimModel()
 {
   ATH_MSG_DEBUG("Initializing Fast Sim Model");
 
-
   // Create the FastCaloSim fast simulation model
-  return new FastCaloSim(name(), m_rndmGenSvc, m_randomEngineName,  m_FastCaloSimCaloTransportation, m_FastCaloSimCaloExtrapolation, m_G4CaloTransportTool, m_FastCaloSimSvc, m_CaloCellContainerSDName, m_doG4Transport, this);
+  return new FastCaloSim(name(), m_rndmGenSvc, m_randomEngineName,  m_FastCaloSimCaloTransportation, m_FastCaloSimCaloExtrapolation, m_G4CaloTransportTool, m_FastCaloSimSvc, m_CaloCellContainerSDName, m_doG4Transport, m_doPhotons, m_doElectrons, m_doHadrons, m_AbsEtaMin, m_AbsEtaMax, m_EkinMin, m_EkinMax, m_doEMECFCS, this);
 }
diff --git a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.h b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.h
index cc30c0cfc1a..93dd2748b1f 100644
--- a/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.h
+++ b/Simulation/G4Utilities/G4FastSimulation/src/FastCaloSimTool.h
@@ -56,6 +56,17 @@ protected:
   Gaudi::Property<std::string> m_CaloCellContainerSDName{this, "CaloCellContainerSDName", "", "Name of the associated CaloCellContainerSD"};
   // Flag to enable G4 transportation
   Gaudi::Property<bool> m_doG4Transport{this, "doG4Transport", false, "Flag to enable G4 transportation"};
+  /// @brief Optional flags that allow to further limit the usage of fast simulation beyond the default configuration
+  /// Note: For fast simulation jobs, default values should be used. They are only useful in special cases
+  /// where you want to further limit the usage of fast sim, e.g. if we want to replace certain parts of full simulation
+  Gaudi::Property<bool> m_doPhotons{this, "doPhotons", true, "Flag to enable FCS simulation for photons"};
+  Gaudi::Property<bool> m_doElectrons{this, "doElectrons", true, "Flag to enable FCS simulation for electrons and positrons"};
+  Gaudi::Property<bool> m_doHadrons{this, "doHadrons", true, "Flag to enable FCS simulation for pions and other hadrons"};
+  Gaudi::Property<float> m_AbsEtaMin{this, "AbsEtaMin", 0, "Abs(Eta) lower bound for FastCaloSim"};
+  Gaudi::Property<float> m_AbsEtaMax{this, "AbsEtaMax", 10, "Abs(Eta) upper bound for FastCaloSim"};
+  Gaudi::Property<float> m_EkinMin{this, "EkinMin", 0, "Kinetic energy lower bound for FastCaloSim"};
+  Gaudi::Property<float> m_EkinMax{this, "EkinMax", std::numeric_limits<float>::max(), "Kinetic energy upper bound for FastCaloSim"};
+  Gaudi::Property<bool> m_doEMECFCS{this, "doEMECFCS", false, "Run FCS in EMEC region while G4 in the rest region"};
 };
 
 #endif //G4FASTSIMULATION_FASTCALOSIMTOOL_H
diff --git a/Simulation/SimulationConfig/python/SimConfigFlags.py b/Simulation/SimulationConfig/python/SimConfigFlags.py
index 6ceaecd6f98..3086b3d98c9 100644
--- a/Simulation/SimulationConfig/python/SimConfigFlags.py
+++ b/Simulation/SimulationConfig/python/SimConfigFlags.py
@@ -238,6 +238,9 @@ def createSimConfigFlags():
 
     scf.addFlag("Sim.FastShower.InputCollection", "TruthEvent") # StoreGate collection name of modified TruthEvent for legacy FastCaloSim use
 
+    # Config FastCaloSim scheme
+    scf.addFlag("Sim.FastCalo.doEMECFCS", False)
+
     # FastChain
     # Setting the BCID for Out-of-Time PU events, list of int
     scf.addFlag("Sim.FastChain.BCID", [1])
diff --git a/Simulation/SimulationConfig/python/SimulationMetadata.py b/Simulation/SimulationConfig/python/SimulationMetadata.py
index d5e988b4d74..e2033c3def0 100644
--- a/Simulation/SimulationConfig/python/SimulationMetadata.py
+++ b/Simulation/SimulationConfig/python/SimulationMetadata.py
@@ -35,6 +35,10 @@ def fillAtlasMetadata(flags, dbFiller):
             if "SimplifiedGeoPath" in flag and not flags.Sim.SimplifiedGeoPath:
                 # This flag is only written to metadata in case the FastCaloSim simplified geometry path is set
                 continue
+            if "FastCalo.doEMECFCS" in flag and not flags.Sim.FastCalo.doEMECFCS:
+                # This flag is only written to metadata in case FCS is used in EMEC region (1.5 < AbsEta < 3.2, Ekin < 8 GeV) is set
+                continue
+
             key = flag.split(".")[-1] #use final part of flag as the key
             value = flags._get(flag)
             if isinstance(value, FlagEnum):
-- 
GitLab