Merge pull request #38849 from VinInn/HackDormandPrince

Optimized version of G4 Dormand Prince RK stepper
cms-sw · Aug 3, 2022 · be16b47 · be16b47
2 parents f061890 + 1734467
commit be16b47
Show file tree

Hide file tree

Showing 6 changed files with 555 additions and 27 deletions.
diff --git a/SimG4Core/MagneticField/interface/Field.h b/SimG4Core/MagneticField/interface/Field.h
@@ -3,14 +3,18 @@
 
 #include "G4MagneticField.hh"
 
-class MagneticField;
+#include "MagneticField/Engine/interface/MagneticField.h"
+
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+
+#include "CLHEP/Units/SystemOfUnits.h"
 
 namespace sim {
-  class Field : public G4MagneticField {
+  class Field final : public G4MagneticField {
   public:
     Field(const MagneticField *f, double d);
     ~Field() override;
-    void GetFieldValue(const G4double p[4], G4double b[3]) const override;
+    inline void GetFieldValue(const G4double p[], G4double b[3]) const override;
 
   private:
     const MagneticField *theCMSMagneticField;
@@ -20,4 +24,26 @@ namespace sim {
     mutable double oldb[3];
   };
 };  // namespace sim
+
+void sim::Field::GetFieldValue(const G4double xyz[], G4double bfield[3]) const {
+  if (std::abs(oldx[0] - xyz[0]) > theDelta || std::abs(oldx[1] - xyz[1]) > theDelta ||
+      std::abs(oldx[2] - xyz[2]) > theDelta) {
+    constexpr float lunit = (1.0 / CLHEP::cm);
+    GlobalPoint ggg((float)(xyz[0]) * lunit, (float)(xyz[1]) * lunit, (float)(xyz[2]) * lunit);
+    GlobalVector v = theCMSMagneticField->inTesla(ggg);
+
+    constexpr float btesla = CLHEP::tesla;
+    oldb[0] = (v.x() * btesla);
+    oldb[1] = (v.y() * btesla);
+    oldb[2] = (v.z() * btesla);
+    oldx[0] = xyz[0];
+    oldx[1] = xyz[1];
+    oldx[2] = xyz[2];
+  }
+
+  bfield[0] = oldb[0];
+  bfield[1] = oldb[1];
+  bfield[2] = oldb[2];
+}
+
 #endif