Merge pull request #726 from OpenGATE/spect_config2

Spect FreeFlight cont.

Author: dsarrut

core/opengate_core/opengate_lib/GateAcceptanceAngleTesterManager.cpp renamed to core/opengate_core/opengate_lib/GateAcceptanceAngleManager.cpp

@@ -5,36 +5,48 @@
    See LICENSE.md for further details
    -------------------------------------------------- */
 
-#include "GateAcceptanceAngleTesterManager.h"
+#include "GateAcceptanceAngleManager.h"
 #include "G4RunManager.hh"
-#include "GateAcceptanceAngleTester.h"
+#include "GateAcceptanceAngleSingleVolume.h"
 #include "GateHelpersDict.h"
 
-GateAcceptanceAngleTesterManager::GateAcceptanceAngleTesterManager() {
+GateAcceptanceAngleManager::GateAcceptanceAngleManager() {
   fEnabledFlag = false;
   fNotAcceptedEvents = 0;
   fAALastRunId = -1;
   fPolicy = AASkipEvent;
-  fMaxNotAcceptedEvents = 100000;
+  fMaxNotAcceptedEvents = 10000;
 }
 
-void GateAcceptanceAngleTesterManager::Initialize(py::dict puser_info,
-                                                  bool is_valid_type) {
-  fAcceptanceAngleVolumeNames = DictGetVecStr(puser_info, "volumes");
+GateAcceptanceAngleManager::~GateAcceptanceAngleManager() {}
+
+void GateAcceptanceAngleManager::Initialize(py::dict user_info,
+                                            bool is_valid_type) {
+  // AA is enabled if volumes is not empty and one of the flags is True
+  // intersection_flag or normal_flag
+  fAcceptanceAngleVolumeNames = DictGetVecStr(user_info, "volumes");
   fEnabledFlag = !fAcceptanceAngleVolumeNames.empty();
+
+  bool b2 = DictGetBool(user_info, "intersection_flag");
+  bool b3 = DictGetBool(user_info, "normal_flag");
+
+  fEnabledFlag = fEnabledFlag && (b2 || b3);
+
   if (!fEnabledFlag)
     return;
   // (we cannot use py::dict here as it is lost at the end of the function)
-  fAcceptanceAngleParam = DictToMap(puser_info);
-  auto s = DictGetStr(puser_info, "skip_policy");
+  fAcceptanceAngleParam = DictToMap(user_info);
+  auto s = DictGetStr(user_info, "skip_policy");
+  fMaxNotAcceptedEvents = DictGetInt(user_info, "max_rejection");
+
   fPolicy = AAUndefined;
   if (s == "ZeroEnergy")
     fPolicy = AAZeroEnergy;
   if (s == "SkipEvents")
     fPolicy = AASkipEvent;
   if (fPolicy == AAUndefined) {
     std::ostringstream oss;
-    oss << "Unknown '" << s << "' mode for GateAcceptanceAngleTesterManager. "
+    oss << "Unknown '" << s << "' mode for GateAcceptanceAngleManager. "
         << "Expected: ZeroEnergy or SkipEvents";
     Fatal(oss.str());
   }
@@ -48,13 +60,14 @@ void GateAcceptanceAngleTesterManager::Initialize(py::dict puser_info,
   }
 }
 
-void GateAcceptanceAngleTesterManager::InitializeAcceptanceAngle() {
+void GateAcceptanceAngleManager::InitializeAcceptanceAngle() {
   if (!fEnabledFlag)
     return;
   // Create the testers (only the first time)
   if (fAATesters.empty()) {
     for (const auto &name : fAcceptanceAngleVolumeNames) {
-      auto *t = new GateAcceptanceAngleTester(name, fAcceptanceAngleParam);
+      auto *t =
+          new GateAcceptanceAngleSingleVolume(name, fAcceptanceAngleParam);
       fAATesters.push_back(t);
     }
   }
@@ -68,34 +81,33 @@ void GateAcceptanceAngleTesterManager::InitializeAcceptanceAngle() {
   fEnabledFlag = !fAcceptanceAngleVolumeNames.empty();
 }
 
-unsigned long
-GateAcceptanceAngleTesterManager::GetNumberOfNotAcceptedEvents() const {
+unsigned long GateAcceptanceAngleManager::GetNumberOfNotAcceptedEvents() const {
   return fNotAcceptedEvents;
 }
 
-bool GateAcceptanceAngleTesterManager::TestIfAccept(
+bool GateAcceptanceAngleManager::TestIfAccept(
     const G4ThreeVector &position, const G4ThreeVector &momentum_direction) {
   if (!fEnabledFlag)
     return true;
 
-  // Loop on all volume to check if it at least one volume is accepted
-  for (auto *tester : fAATesters) {
+  // Loop on all the volumes to check if it at least one volume is accepted
+  for (const auto *tester : fAATesters) {
     bool accept = tester->TestIfAccept(position, momentum_direction);
     if (accept)
       return true;
   }
-  fNotAcceptedEvents++;
   if (fNotAcceptedEvents > fMaxNotAcceptedEvents) {
     std::ostringstream oss;
     oss << "Error, in AcceptanceAngleTest: " << fNotAcceptedEvents
         << " trials has been tested without accepted angle ; probably no "
            "possible direction here. Abort. ";
     Fatal(oss.str());
   }
+  fNotAcceptedEvents++;
   return false;
 }
 
-void GateAcceptanceAngleTesterManager::StartAcceptLoop() {
+void GateAcceptanceAngleManager::StartAcceptLoop() {
   if (!fEnabledFlag)
     return;
   fNotAcceptedEvents = 0;

core/opengate_core/opengate_lib/GateAcceptanceAngleTesterManager.h renamed to core/opengate_core/opengate_lib/GateAcceptanceAngleManager.h

@@ -5,19 +5,21 @@
    See LICENSE.md for further details
    -------------------------------------------------- */
 
-#ifndef GateAcceptanceAngleTesterManager_h
-#define GateAcceptanceAngleTesterManager_h
+#ifndef GateAcceptanceAngleManager_h
+#define GateAcceptanceAngleManager_h
 
-#include "GateAcceptanceAngleTester.h"
+#include "GateAcceptanceAngleSingleVolume.h"
 #include "GateHelpers.h"
 
-class GateAcceptanceAngleTesterManager {
+class GateAcceptanceAngleManager {
 public:
-  GateAcceptanceAngleTesterManager();
+  GateAcceptanceAngleManager();
+
+  ~GateAcceptanceAngleManager();
 
   enum AAPolicyType { AAZeroEnergy, AASkipEvent, AAUndefined };
 
-  void Initialize(py::dict puser_info, bool is_iso);
+  void Initialize(py::dict user_info, bool is_valid_type);
 
   void InitializeAcceptanceAngle();
 
@@ -34,12 +36,12 @@ class GateAcceptanceAngleTesterManager {
 
   AAPolicyType fPolicy;
   std::map<std::string, std::string> fAcceptanceAngleParam;
-  std::vector<GateAcceptanceAngleTester *> fAATesters{};
+  std::vector<GateAcceptanceAngleSingleVolume *> fAATesters{};
   std::vector<std::string> fAcceptanceAngleVolumeNames;
   bool fEnabledFlag;
   unsigned long fNotAcceptedEvents;
   unsigned long fMaxNotAcceptedEvents;
   int fAALastRunId;
 };
 
-#endif // GateAcceptanceAngleTesterManager_h
+#endif // GateAcceptanceAngleManager_h

core/opengate_core/opengate_lib/GateAcceptanceAngleSingleVolume.cpp

@@ -0,0 +1,123 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   This software is distributed under the terms
+   of the GNU Lesser General  Public Licence (LGPL)
+   See LICENSE.md for further details
+   -------------------------------------------------- */
+
+#include "GateAcceptanceAngleSingleVolume.h"
+#include "G4LogicalVolumeStore.hh"
+#include "G4Navigator.hh"
+#include "G4PhysicalVolumeStore.hh"
+#include "GateHelpersDict.h"
+#include "GateHelpersImage.h"
+
+GateAcceptanceAngleSingleVolume::GateAcceptanceAngleSingleVolume(
+    const std::string &volume,
+    const std::map<std::string, std::string> &param) {
+  fAcceptanceAngleVolumeName = volume;
+  fAASolid = nullptr;
+  fAANavigator = nullptr;
+  fAARotation = nullptr;
+  fIntersectionFlag = false;
+  fNormalFlag = false;
+  fNormalAngleTolerance = 0;
+  fMinDistanceNormalAngleTolerance = 0;
+  fDistanceDependentAngleToleranceFlag = false;
+  fAngle1 = 0;
+  fAngle2 = 0;
+  fDistance1 = 0;
+  fDistance2 = 0;
+
+  // Retrieve the solid
+  const auto lvs = G4LogicalVolumeStore::GetInstance();
+  const auto lv = lvs->GetVolume(fAcceptanceAngleVolumeName);
+  fAASolid = lv->GetSolid();
+
+  // Init a navigator that will be used to find the transform
+  const auto pvs = G4PhysicalVolumeStore::GetInstance();
+  const auto world = pvs->GetVolume("world");
+  fAANavigator = new G4Navigator();
+  fAANavigator->SetWorldVolume(world);
+
+  // parameters
+  fIntersectionFlag = StrToBool(param.at("intersection_flag"));
+  fNormalFlag = StrToBool(param.at("normal_flag"));
+  fNormalAngleTolerance = StrToDouble(param.at("normal_tolerance"));
+  fMinDistanceNormalAngleTolerance =
+      StrToDouble(param.at("normal_tolerance_min_distance"));
+  fNormalVector = StrToG4ThreeVector(param.at("normal_vector"));
+
+  if (fNormalAngleTolerance <= 0) {
+    Fatal("Normal angle tolerance must be strictly positive while it is " +
+          std::to_string(fNormalAngleTolerance));
+  }
+
+  // DD -> not recommended (too slow), prefer normal_tolerance_min_distance
+  fDistanceDependentAngleToleranceFlag =
+      StrToBool(param.at("distance_dependent_normal_tolerance"));
+  fAngle1 = StrToDouble(param.at("angle1"));
+  fAngle2 = StrToDouble(param.at("angle2"));
+  fDistance1 = StrToDouble(param.at("distance1"));
+  fDistance2 = StrToDouble(param.at("distance2"));
+  if (fDistance1 >= fDistance2) {
+    Fatal(
+        "For 'distance_dependent_normal_tolerance', Distance1 must be strictly "
+        "less than distance2 " +
+        std::to_string(fDistance1) + " " + std::to_string(fDistance2));
+  }
+  a = (1.0 / tan(fAngle1) - 1.0 / tan(fAngle2)) / (fDistance1 - fDistance2);
+  b = 1.0 / tan(fAngle2) - a * fDistance2;
+}
+
+GateAcceptanceAngleSingleVolume::~GateAcceptanceAngleSingleVolume() {
+  delete fAARotation;
+}
+
+void GateAcceptanceAngleSingleVolume::UpdateTransform() {
+
+  if (fAARotation == nullptr)
+    fAARotation = new G4RotationMatrix;
+
+  // Get the transformation
+  G4ThreeVector tr;
+  ComputeTransformationFromWorldToVolume(fAcceptanceAngleVolumeName, tr,
+                                         *fAARotation, true);
+  // It is not fully clear why the AffineTransform need the inverse
+  fAATransform = G4AffineTransform(fAARotation->inverse(), tr);
+}
+
+bool GateAcceptanceAngleSingleVolume::TestIfAccept(
+    const G4ThreeVector &position,
+    const G4ThreeVector &momentum_direction) const {
+  const auto localDirection = (*fAARotation) * (momentum_direction);
+  double dist = 0;
+  if (fIntersectionFlag || fDistanceDependentAngleToleranceFlag) {
+    const auto localPosition = fAATransform.TransformPoint(position);
+    dist = fAASolid->DistanceToIn(localPosition, localDirection);
+  }
+  if (fIntersectionFlag) {
+    if (dist == kInfinity)
+      return false;
+  }
+  if (fNormalFlag || fDistanceDependentAngleToleranceFlag) {
+    const auto angle = fNormalVector.angle(localDirection);
+    if (dist < fMinDistanceNormalAngleTolerance)
+      return true;
+    if (fDistanceDependentAngleToleranceFlag)
+      return DistanceDependentToleranceTest(angle, dist);
+    return angle < fNormalAngleTolerance;
+  }
+  return true;
+}
+
+bool GateAcceptanceAngleSingleVolume::DistanceDependentToleranceTest(
+    // This is very slow, (3x than other methods
+    // We recommend to not use this method
+    const double angle, const double dist) const {
+  const double tol = atan(1.0 / (a * dist + b));
+
+  if (tol < 0)
+    return true;
+  return angle < tol;
+}

core/opengate_core/opengate_lib/GateAcceptanceAngleTester.h renamed to core/opengate_core/opengate_lib/GateAcceptanceAngleSingleVolume.h

@@ -5,34 +5,44 @@
    See LICENSE.md for further details
    -------------------------------------------------- */
 
-#ifndef GateAcceptanceAngleTester_h
-#define GateAcceptanceAngleTester_h
+#ifndef GateAcceptanceAngleSingleVolume_h
+#define GateAcceptanceAngleSingleVolume_h
 
 #include "G4AffineTransform.hh"
 #include "GateHelpers.h"
 
-class GateAcceptanceAngleTester {
+class GateAcceptanceAngleSingleVolume {
 public:
-  GateAcceptanceAngleTester(const std::string &volume,
-                            std::map<std::string, std::string> &param);
+  GateAcceptanceAngleSingleVolume(
+      const std::string &volume,
+      const std::map<std::string, std::string> &param);
 
-  ~GateAcceptanceAngleTester();
+  ~GateAcceptanceAngleSingleVolume();
 
   bool TestIfAccept(const G4ThreeVector &position,
                     const G4ThreeVector &momentum_direction) const;
 
   void UpdateTransform();
+  bool DistanceDependentToleranceTest(double angle, double dist) const;
 
 protected:
   std::string fAcceptanceAngleVolumeName;
   bool fIntersectionFlag;
   bool fNormalFlag;
   double fNormalAngleTolerance;
+  bool fDistanceDependentAngleToleranceFlag;
+  double fMinDistanceNormalAngleTolerance;
+  double fAngle1;
+  double fAngle2;
+  double fDistance1;
+  double fDistance2;
+  double a;
+  double b;
   G4ThreeVector fNormalVector;
   G4AffineTransform fAATransform;
   G4RotationMatrix *fAARotation;
   G4VSolid *fAASolid;
   G4Navigator *fAANavigator;
 };
 
-#endif // GateAcceptanceAngleTester_h
+#endif // GateAcceptanceAngleSingleVolume_h