ATLFieldManager for parameter optimization

FullSimLight/include/ATLFieldManager.hh

+// ATLFieldManager.h
+
+#ifndef ATLFieldManager_h
+#define ATLFieldManager_h 1
+
+#include "G4FieldManager.hh"
+
+#include <string>
+#include <unordered_map>
+#include <memory>  // Add this include for std::unique_ptr
+#include "G4Types.hh"
+#include "G4MagneticField.hh"
+
+#include <vector>
+#include <sstream>
+#include <iterator>
+
+
+//#include "G4bool.hh"   // Include this for G4bool
+#include "G4ThreeVector.hh"
+#include "G4SystemOfUnits.hh"
+
+class G4ChordFinder;
+class G4PropagatorInField;
+class G4Track;
+class G4Region;
+class MagneticField;
+
+
+namespace atl {
+    class Field : public G4MagneticField {
+    public:
+        // Constructor
+        Field(const MagneticField* f, double d);
+
+        // Destructor
+        ~Field() override;
+
+        void GetFieldValue(const G4double* pos, G4double* B) const override;
+
+    private:
+        const MagneticField* theATLMagneticField;
+        double theDelta;
+    };
+}
+
+// Class ParameterSet
+/*
+class ParameterSet {
+public:
+    ParameterSet(const std::unordered_map<std::string, double>& params)
+        : parameters(params) {}
+
+    double getDoubleParameter(const std::string& key) const {
+        auto it = parameters.find(key);
+        return (it != parameters.end()) ? it->second : 0.0;
+    }
+
+private:
+    std::unordered_map<std::string, double> parameters;
+};  */
+
+
+class ParameterSet {
+public:
+    ParameterSet(const std::unordered_map<std::string, double>& params)
+        : parameters(params) {}
+
+    double getDoubleParameter(const std::string& key) const {
+    auto it = parameters.find(key);
+    return (it != parameters.end()) ? it->second : 0.0;
+}
+    std::string getStringParameter(const std::string& key) const {
+        return parameters.count(key) ? std::to_string(parameters.at(key)) : "";
+    }
+
+    template <typename T>
+    T getParameter(const std::string& key) const;
+    
+private:
+    std::unordered_map<std::string, double> parameters;
+};
+
+template <>
+std::vector<std::string> ParameterSet::getParameter<std::vector<std::string>>(const std::string& key) const {
+    std::vector<std::string> result;
+    if (parameters.count(key)) {
+        std::istringstream iss(getStringParameter(key));
+        std::copy(
+            std::istream_iterator<std::string>(iss),
+            std::istream_iterator<std::string>(),
+            std::back_inserter(result)
+        );
+    }
+    return result;
+};
+
+// Define ParameterSet class
+/*
+class ParameterSet {
+public:
+    ParameterSet(const std::unordered_map<std::string, double>& params,
+                  const std::unordered_map<std::string, std::string>& params_str)
+        : parameters(params), parameters_str(params_str) {}
+
+    double getDoubleParameter(const std::string& key) const {
+        auto it = parameters.find(key);
+        return (it != parameters.end()) ? it->second : 0.0;
+    }
+
+    std::string getStringParameter(const std::string& key) const {
+        auto it = parameters_str.find(key);
+        return (it != parameters_str.end()) ? it->second : "";
+    }
+
+private:
+    std::unordered_map<std::string, double> parameters;
+    std::unordered_map<std::string, std::string> parameters_str;
+};  */
+
+//class ATLFieldManager 
+class ATLFieldManager : public G4FieldManager{
+public:
+    ATLFieldManager();
+    ATLFieldManager(const ATLFieldManager&) = delete;
+    ATLFieldManager(ATLFieldManager&&) = default;
+    ATLFieldManager& operator=(const ATLFieldManager&) = delete;
+    ATLFieldManager& operator=(ATLFieldManager&&) = default;
+    virtual ~ATLFieldManager() = default;
+
+    void InitialiseForVolume(const ParameterSet& params,
+                             atl::Field* field,
+                             G4ChordFinder* cfDefault,
+                            // G4ChordFinder* cfMonopole,
+                             const std::string& vol,
+                             const std::string& fieldType,
+                             const std::string& stepperName,
+                             double delta,
+                             G4PropagatorInField*);
+
+    void setMonopoleTracking(G4bool);
+
+private:
+    bool isInsideVacuum(const G4Track*);
+    bool isInsideTracker(const G4Track*);
+    void setDefaultChordFinder();
+    void setChordFinderForTracker();
+    void setChordFinderForBeampipe();
+    void ConfigureForTrack(const G4Track*) ;
+  //  void InitialiseForVolume();
+
+    std::unique_ptr<atl::Field> theField;
+    G4ChordFinder* m_currChordFinder;
+    G4ChordFinder* m_chordFinder;
+    G4ChordFinder* m_chordFinderMonopole;
+    G4PropagatorInField* m_propagator;
+    std::vector<const G4Region*> m_regions;
+
+    double m_dChordDefault;
+    double m_dIntersectionDefault;
+    double m_dOneStepDefault;
+    double m_stepMaxDefault;
+
+    double m_dChordTracker;
+    double m_dIntersectionTracker;
+    double m_dOneStepTracker;
+    double m_stepMaxTracker;
+
+    double m_Rmax2;
+    double m_Zmax;
+    double m_energyThTracker;
+    double m_energyThreshold;
+    
+    double m_dChordBeampipe;
+    double m_dIntersectionBeampipe;
+    double m_dOneStepBeampipe;
+    double m_stepMaxBeampipe;
+
+    bool m_cfTracker;
+    bool m_cfBeampipe;
+};
+
+#endif  // ATLFieldManager_h
+

FullSimLight/src/ATLFieldManager.cc

+
+#include "ATLFieldManager.hh"
+#include "G4FieldManager.hh"
+#include "G4Types.hh"
+
+#include "FSLDetectorConstruction.hh"
+#include "FSLDetectorMessenger.hh"
+#include "RegionConfigurator.hh"
+#include "FullSimLight/MagFieldPlugin.h"
+#include "MassCalculator.hh"
+#include "ClashDetector.hh"
+
+// Geant4 includes
+#include "G4Version.hh"
+#include "globals.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4VisAttributes.hh"
+
+#include "G4UniformMagField.hh"
+#include "G4QuadrupoleMagField.hh"
+#include "G4FieldManager.hh"
+#include "G4TransportationManager.hh"
+#include "G4VPhysicalVolume.hh"
+#include "G4RunManager.hh"
+#include "G4PVPlacement.hh"
+#include "G4ChordFinder.hh"
+#include "G4Mag_UsualEqRhs.hh"
+#include "G4MagIntegratorStepper.hh"
+#include "G4Version.hh"
+#include "G4LogicalVolume.hh"
+#include "G4UnitsTable.hh"
+#include "G4GDMLParser.hh"
+//#include "G4LogicalVolumeStore.hh"
+
+//chord
+//#include "G4MagneticField.hh"
+//#include "G4IntegrationDriver.hh"
+//#include "G4VIntegrationDriver.hh"
+//#include "G4MagIntegratorStepper.hh"
+#include "G4Field.hh"
+#include "G4VIntegrationDriver.hh"
+#include "G4VIntegrationDriver.hh"
+#include "G4FSALIntegrationDriver.hh"
+#include "G4VFSALIntegrationStepper.hh"
+#include "G4RK547FEq1.hh"
+// #include "G4RK547FEq2.hh"
+// #include "G4RK547FEq3.hh"
+#include "G4NystromRK4.hh"
+
+// New FSAL type driver / steppers -----
+#include "G4IntegrationDriver.hh"
+#include "G4InterpolationDriver.hh"
+// #include "G4FSALBogackiShampine45.hh"
+// #include "G4FSALDormandPrince745.hh"
+#include "G4HelixHeum.hh"
+#include "G4BFieldIntegrationDriver.hh"
+
+#include <cassert>
+
+// Geant4 steppers
+#if G4VERSION_NUMBER>=1040
+#include "G4BogackiShampine23.hh"
+#include "G4BogackiShampine45.hh"
+#include "G4DoLoMcPriRK34.hh"
+#include "G4DormandPrince745.hh"
+#include "G4DormandPrinceRK56.hh"
+#include "G4DormandPrinceRK78.hh"
+#include "G4RK547FEq1.hh"
+#include "G4RK547FEq2.hh"
+#include "G4RK547FEq3.hh"
+#include "G4TsitourasRK45.hh"
+#include "G4VIntegrationDriver.hh"
+#include "G4IntegrationDriver.hh"
+#endif
+
+#if G4VERSION_NUMBER>=1060
+#include "G4InterpolationDriver.hh"
+#endif
+
+#include "G4CashKarpRKF45.hh"
+#include "G4ClassicalRK4.hh"
+#include "G4FieldManager.hh"
+#include "G4HelixExplicitEuler.hh"
+#include "G4HelixImplicitEuler.hh"
+//#include "G4HelixBeampipeRunge.hh"
+#include "G4NystromRK4.hh"
+#include "G4RKG3_Stepper.hh"
+
+
+// **** INCLUDES for GeoModel
+#include "GeoModelRead/ReadGeoModel.h"
+#include "GeoModel2G4/ExtParameterisedVolumeBuilder.h"
+#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoPhysVol.h"
+#include "GeoModelKernel/GeoFullPhysVol.h"
+#include "GeoModelKernel/GeoNameTag.h"
+#include "GeoModelKernel/GeoVGeometryPlugin.h"
+#include "GeoModelKernel/GeoGeometryPluginLoader.h"
+#include "GeoModelKernel/GeoVolumeCursor.h"
+
+// For Sensitive Detector plugins:
+#include "FullSimLight/FSLSensitiveDetectorPlugin.h"
+#include "FSLSDPluginLoader.h"
+#include "G4VSensitiveDetector.hh"
+#include "G4LogicalVolumeStore.hh"
+#include "G4LogicalVolume.hh"
+#include "G4SDManager.hh"
+#include "G4UIcommand.hh"
+#include <iomanip>
+//#include "G4ChordFinder.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4MagneticField.hh"
+#include "G4Mag_UsualEqRhs.hh"
+#include "G4MagIntegratorDriver.hh"
+#include "G4DormandPrince745.hh"
+
+// New FSAL type driver / steppers -----
+#include "G4FSALIntegrationDriver.hh"
+#include "G4VFSALIntegrationStepper.hh"
+#include "G4RK547FEq1.hh"
+#include "G4NystromRK4.hh"
+
+// New FSAL type driver / steppers -----
+#include "G4IntegrationDriver.hh"
+#include "G4InterpolationDriver.hh"
+// #include "G4FSALBogackiShampine45.hh"
+// #include "G4FSALDormandPrince745.hh"
+#include "G4HelixHeum.hh"
+#include "G4BFieldIntegrationDriver.hh"
+
+#include "G4PropagatorInField.hh"
+#include "G4UserLimits.hh"
+#include "G4StepLimiter.hh"
+//#include "StepMax.hh"
+#include <cassert>
+#include "HistoManager.hh"
+
+#include "G4RegionStore.hh"
+#include "G4Region.hh"
+#include "G4Track.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4DynamicParticle.hh"
+#include "G4SteppingManager.hh"
+#include "G4Track.hh"
+#include "G4Step.hh"
+#include "G4ThreeVector.hh"
+////////
+#include "CLHEP/Units/GlobalSystemOfUnits.h"
+#include "G4ChordFinder.hh"
+#include "G4MagIntegratorStepper.hh"
+#include "G4PropagatorInField.hh"
+#include "G4Region.hh"
+#include "G4RegionStore.hh"
+#include "G4Track.hh"
+
+#include "G4Types.hh"
+#include "G4ThreeVector.hh"
+#include "G4Region.hh"
+
+ATLFieldManager::ATLFieldManager() 	
+    : G4FieldManager(),  // G4FieldManager is the base class 
+      m_currChordFinder(nullptr),
+      m_chordFinder(nullptr),
+      m_propagator(nullptr),
+      m_dChordDefault(0.25),	//Default
+      m_dIntersectionDefault(0.001),
+      m_dOneStepDefault(0.01),
+      m_stepMaxDefault(1000.), ////limit
+      m_dChordTracker(0.0001),	//Tracker					//tracker* !!?
+      m_dIntersectionTracker(0.001),					//tracker* !!?
+      m_dOneStepTracker(0.001),					//tracker* !!?
+      m_stepMaxTracker(100.),
+
+      m_Rmax2(1.e+4),
+      m_Zmax(3.e+2),
+
+      m_energyThTracker(0.0),					//tracker* !!?
+      m_energyThreshold(0.0),
+
+      m_dChordBeampipe(0.002),	//Beampipe
+      m_dIntersectionBeampipe(0.0001),
+      m_dOneStepBeampipe(0.02),
+      m_stepMaxBeampipe(1000.0),
+      
+      m_cfTracker(false),					//tracker* !!?
+      m_cfBeampipe(false) {							//beampipe*
+       
+  if (m_chordFinder != m_currChordFinder) {
+    delete m_chordFinder;
+  }
+}
+ 
+void ATLFieldManager::InitialiseForVolume(const ParameterSet &params,
+                                          atl::Field* field,
+                                          G4ChordFinder* cfDefault,
+                                      //    G4ChordFinder* cfMonopole,
+                                          const std::string& vol,
+                                          const std::string& fieldType,
+                                          const std::string& stepperName,
+                                          double delta,
+                                          G4PropagatorInField* pf) {
+
+   double minstep = params.getDoubleParameter("Minstep") * CLHEP::mm;		//define or intilize in 1st section with other parameters
+   double minEpsStep = params.getDoubleParameter("MinimumEpsilonStep") * CLHEP::mm;
+   double maxEpsStep = params.getDoubleParameter("MaximumEpsilonStep") * CLHEP::mm;
+   //int maxLC = static_cast<int>(params.getDoubleParameter("MaximumLoopCounts", 1000));
+   int maxLC = static_cast<int>(params.getDoubleParameter("MaximumLoopCounts"));  // Assuming "MaximumLoopCounts" is the parameter name
+
+   // Tracker
+   m_dChordTracker = params.getDoubleParameter("DeltaChord") * CLHEP::mm;					//tracker* !!?
+   m_dIntersectionTracker = params.getDoubleParameter("DeltaIntersection") * CLHEP::mm;		
+   m_dOneStepTracker = params.getDoubleParameter("DeltaOneStep") * CLHEP::mm;					//tracker* !!?
+   m_stepMaxTracker = params.getDoubleParameter("MaxStep") * CLHEP::cm;
+
+   //Beampipe
+   m_dChordBeampipe = params.getDoubleParameter("DeltaChord") * CLHEP::mm;
+   m_dIntersectionBeampipe = params.getDoubleParameter("DeltaIntersection") * CLHEP::mm;
+   m_dOneStepBeampipe = params.getDoubleParameter("DeltaOneStep") * CLHEP::mm;
+   m_stepMaxBeampipe = params.getDoubleParameter("MaxStep") * CLHEP::cm; 
+   //Default
+   m_dChordDefault = params.getDoubleParameter("DeltaChord") * CLHEP::mm;	
+   m_dIntersectionDefault = params.getDoubleParameter("DeltaIntersection") * CLHEP::mm;
+   m_dOneStepDefault = params.getDoubleParameter("DeltaOneStep") * CLHEP::mm;
+   m_stepMaxDefault = params.getDoubleParameter("MaxStep") * CLHEP::cm;
+  
+   m_energyThreshold = params.getDoubleParameter("EnergyThBeampipe") * CLHEP::GeV;	//...check if ....
+   m_energyThTracker = params.getDoubleParameter("EnergyThTracker") * CLHEP::GeV;					//tracker* !!?
+ 
+   double rmax = params.getDoubleParameter("RmaxTracker") * CLHEP::mm;					//tracker* !!?
+   m_Rmax2 = rmax * rmax;
+   m_Zmax = params.getDoubleParameter("ZmaxTracker") * CLHEP::mm;					//tracker* !!?
+ 
+   std::cout << "is it workinggggggggg!" ">\n"
+   //edm::LogVerbatim("")
+    //G4cout
+       << " New ATLFieldManager: LogicalVolume:      <" << vol << ">\n"
+       << "               Stepper                   <" << stepperName << ">\n"	//define
+       << "               Field type                 <" << fieldType << ">\n"	//define
+       << "               Field const delta           " << delta << " mm\n"
+       << "               MaximumLoopCounts           " << maxLC << "\n"
+       << "               MinimumEpsilonStep          " << minEpsStep << "\n"
+       << "               MaximumEpsilonStep          " << maxEpsStep << "\n"
+       << "               MinStep                     " << minstep << " mm\n"
+       << "               MaxStep Default             " << m_stepMaxDefault / CLHEP::cm << " cm\n"
+       << "               MaxStep Tracker             " << m_stepMaxTracker / CLHEP::cm << " cm\n"
+       << "               DeltaChord Default          " << m_dChordDefault << " mm\n"
+       << "               DeltaOneStep Default        " << m_dOneStepDefault << " mm\n"
+       << "               DeltaIntersection Default   " << m_dIntersectionDefault << " mm\n"
+       << "               DeltaInterTracker           " << m_dIntersectionTracker << " mm\n"					//tracker* !!?
+       << "               EnergyThresholdBeampipe       " << m_energyThreshold / CLHEP::MeV << " MeV\n"
+       << "               EnergyThresholdTracker      " << m_energyThTracker / CLHEP::MeV << " MeV\n"					//tracker* !!?
+       << "               DeltaChordBeampipe            " << m_dChordBeampipe << " mm\n"
+       << "               DeltaOneStepBeampipe          " << m_dOneStepBeampipe << " mm\n"
+       << "               DeltaIntersectionBeampipe     " << m_dIntersectionBeampipe << " mm\n"
+       << "               MaxStepInVacuum             " << m_stepMaxBeampipe / CLHEP::cm << " cm";					//beampipe*
+ 
+   // initialisation of chord finders
+   m_chordFinder = cfDefault; //here
+   m_chordFinder->SetDeltaChord(m_dChordDefault); 	//check replacement of m_dChord
+ //  m_chordFinderMonopole = cfmon; //here
+ //  m_chordFinderMonopole->SetDeltaChord(m_dChord);  //check
+ 
+   // initialisation of field manager
+   theField.reset(field);
+   SetDetectorField(field);  //here
+
+   SetMinimumEpsilonStep(minEpsStep);
+   SetMaximumEpsilonStep(maxEpsStep);
+ 
+   // propagater in field
+   m_propagator = pf;
+   pf->SetMaxLoopCount(maxLC);
+   pf->SetMinimumEpsilonStep(minEpsStep);
+   pf->SetMaximumEpsilonStep(maxEpsStep);
+ 
+   // define regions
+   std::vector<std::string> rnames = params.getParameter<std::vector<std::string>>("VacRegions");
+   //std::vector<std::string> rnames = params.getStringParameter("VacRegions");
+
+   if (!rnames.empty()) {
+     std::stringstream ss;
+     std::vector<G4Region *> *rs = G4RegionStore::GetInstance();
+     for (auto &regnam : rnames) {
+       for (auto &reg : *rs) {
+         if (regnam == reg->GetName()) {
+           m_regions.push_back(reg);
+           ss << "  " << regnam;
+         }
+       }
+     }
+  //    std::cout << "Beampipe field integration in G4Regions....is it workinggggggggg!" << std::endl;
+      std::cout << "Beampipe field integration in G4Regions:\n" << ss.str() << "\n";
+//edm::LogVerbatim("SimG4CoreApplication") << "Beampipe field integration in G4Regions:\n" << ss.str() << "\n";// edm is for CMS, modify it
+   }
+ }
+ 
+ void ATLFieldManager::ConfigureForTrack(const G4Track *track) {
+   // run time parameters per track
+   if (track->GetKineticEnergy() > m_energyThTracker && isInsideTracker(track)) {					//tracker* !!?
+     if (!m_cfTracker) {					
+       setChordFinderForTracker();					
+     }
+ 
+   } else if (track->GetKineticEnergy() <= m_energyThreshold || isInsideVacuum(track)) {					//beampipe*
+     if (!m_cfBeampipe) {					//beampipe*
+       setChordFinderForBeampipe();					//beampipe*
+     }
+ 
+   } else if (m_cfTracker || m_cfBeampipe) {				//tracker* !!?		//beampipe*
+     // restore defaults
+     setDefaultChordFinder();
+   }
+ }
+ 
+ bool ATLFieldManager::isInsideVacuum(const G4Track *track) {					//beampipe*
+   if (!m_regions.empty()) {
+     const G4Region *reg = track->GetVolume()->GetLogicalVolume()->GetRegion();
+     for (auto &areg : m_regions) {
+       if (reg == areg) {
+         return true;
+       }
+     }
+   }
+   return false;
+ }
+ 
+ bool ATLFieldManager::isInsideTracker(const G4Track *track) {					//tracker* !!?
+   const G4ThreeVector &pos = track->GetPosition();
+   const double x = pos.x();
+   const double y = pos.y();
+   return (x * x + y * y < m_Rmax2 && std::abs(pos.z()) < m_Zmax);
+ }
+ 
+ void ATLFieldManager::setDefaultChordFinder() {
+   if (m_currChordFinder != m_chordFinder) {
+     m_currChordFinder = m_chordFinder;
+     SetChordFinder(m_currChordFinder);
+   }
+   m_currChordFinder->SetDeltaChord(m_dChordDefault);//monopole
+   SetDeltaOneStep(m_dOneStepDefault);
+   SetDeltaIntersection(m_dIntersectionDefault);
+   m_propagator->SetLargestAcceptableStep(m_stepMaxDefault);
+   m_cfBeampipe = m_cfTracker = false;					//tracker* !!?	//beampipe*
+ }
+ 
+ void ATLFieldManager::setChordFinderForTracker() {		//tracker* !!?
+   if (m_currChordFinder != m_chordFinder) {
+     m_currChordFinder = m_chordFinder;
+     SetChordFinder(m_currChordFinder);
+   }
+   m_currChordFinder->SetDeltaChord(m_dChordTracker);					//tracker* !!?
+   SetDeltaOneStep(m_dOneStepTracker);					//tracker* !!?
+   SetDeltaIntersection(m_dIntersectionTracker);					//tracker* !!?
+   m_propagator->SetLargestAcceptableStep(m_stepMaxTracker);
+   m_cfBeampipe = false;					//beampipe*
+   m_cfTracker = true;					//tracker* !!?
+ }
+ 
+ void ATLFieldManager::setChordFinderForBeampipe() { 		//beampipe*
+   if (m_currChordFinder != m_chordFinder) {
+     m_currChordFinder = m_chordFinder;
+     SetChordFinder(m_currChordFinder);
+   }
+   m_currChordFinder->SetDeltaChord(m_dChordBeampipe);
+   SetDeltaOneStep(m_dOneStepBeampipe);
+   SetDeltaIntersection(m_dIntersectionBeampipe);
+   m_propagator->SetLargestAcceptableStep(m_stepMaxBeampipe);
+   m_cfBeampipe = true;
+   m_cfTracker = false;					//tracker* !!?
+ }