Compatibility with dual/multi-mirror dRICH

CMakeLists.txt

@@ -55,6 +55,7 @@ list(FILTER HEADERS EXCLUDE REGEX "LinkDef\\.h$")
 # sources
 set(IRT_SRC
   ${PROJECT_SOURCE_DIR}/source/SphericalSurface.cc
+  ${PROJECT_SOURCE_DIR}/source/SphericalSurfaceWithHalfSpace.cc
   ${PROJECT_SOURCE_DIR}/source/FlatSurface.cc
   ${PROJECT_SOURCE_DIR}/source/IRT.cc
   ${PROJECT_SOURCE_DIR}/source/ChargedParticle.cc

include/CherenkovDetectorCollection.h

@@ -55,13 +55,13 @@ class CherenkovDetectorCollection: public BitMask {
     {
       auto boundary = surface->_Clone(0.0, TVector3(0,0,1));
       boundary->Shift(( thickness/2)*surface->GetNormal());
-      det->AddOpticalBoundary(path, new OpticalBoundary(radiator,                  boundary, true));
+      det->AddOpticalBoundary(path, new OpticalBoundary(radiator,                  boundary, true, false));
       radiator->m_Borders[path].first = boundary;
     }
     {
       auto boundary = surface->_Clone(0.0, TVector3(0,0,1));
       boundary->Shift((-thickness/2)*surface->GetNormal());
-      det->AddOpticalBoundary(path, new OpticalBoundary(det->GetContainerVolume(), boundary, true));
+      det->AddOpticalBoundary(path, new OpticalBoundary(det->GetContainerVolume(), boundary, true, false));
 
       radiator->m_Borders[path].second = boundary;
       // This will most likely be a temporary assignment;
@@ -106,7 +106,7 @@ class CherenkovDetectorCollection: public BitMask {
     // This is most likely a temporary assignment;
     radiator->m_Borders[path].first = surface;
 
-    det->AddOpticalBoundary(path, new OpticalBoundary(FindRadiator(lv), surface, true));
+    det->AddOpticalBoundary(path, new OpticalBoundary(FindRadiator(lv), surface, true, false));
     //det->SetContainerVolume(lv);
     det->SetContainerVolume(radiator);
 

include/CherenkovMirror.h

@@ -71,6 +71,18 @@ class SphericalMirror: public CherenkovMirror, public SphericalSurface {
 #endif
 };
 
+class SphericalMirrorWithHalfSpace: public CherenkovMirror, public SphericalSurfaceWithHalfSpace {
+ public:
+  SphericalMirrorWithHalfSpace() {};
+  SphericalMirrorWithHalfSpace(G4VSolid *solid, G4Material *material, const TVector3 &x0, double r0, const std::vector<std::pair<TVector3,TVector3>> halfSpaces): 
+    CherenkovMirror(solid, material), SphericalSurfaceWithHalfSpace(x0, r0, halfSpaces) {};
+  ~SphericalMirrorWithHalfSpace() {};
+
+#ifndef DISABLE_ROOT_IO
+  ClassDef(SphericalMirrorWithHalfSpace, 2);
+#endif
+};
+
 class CherenkovMirrorGroup: public TObject {
  public:
   CherenkovMirrorGroup() {};

include/CherenkovPhotonDetector.h

@@ -31,7 +31,7 @@ class CherenkovPhotonDetector: public G4Object {
   double GetGeometricEfficiency( void ) const { return m_GeometricEfficiency; };
 
   void AddItselfToOpticalBoundaries(IRT *irt, ParametricSurface *surface) /*const*/ {
-    auto boundary = new OpticalBoundary(0, surface, true);
+    auto boundary = new OpticalBoundary(0, surface, true, false);
     irt->AddOpticalBoundary(boundary);
 
     m_OpticalBoundaryStorage.push_back(boundary);

include/OpticalBoundary.h

@@ -12,10 +12,10 @@ class OpticalBoundary: public TObject {
   friend class IRT;
 
  public:
- OpticalBoundary(): m_Radiator(0), m_Surface(0), m_Refractive(true) {};
+  OpticalBoundary(): m_Radiator(0), m_Surface(0), m_Refractive(true), m_Reflective(false) {};
   // The "0" as a radiator ptr implicitly says "there is no photon life beyond this boundary in IRT";
- OpticalBoundary(/*const*/ CherenkovRadiator *radiator, const ParametricSurface *surface, bool refractive): 
-  m_Radiator(radiator), m_Surface(surface), m_Refractive(refractive) {};
+  OpticalBoundary(/*const*/ CherenkovRadiator *radiator, const ParametricSurface *surface, bool refractive, bool reflective): 
+    m_Radiator(radiator), m_Surface(surface), m_Refractive(refractive), m_Reflective(reflective) {};
   ~OpticalBoundary() {};
 
   CherenkovRadiator *GetRadiator( void ) const {
@@ -32,7 +32,7 @@ class OpticalBoundary: public TObject {
 
   // Boundary surface; either refractive or reflective; 
   const ParametricSurface *m_Surface;
-  bool m_Refractive;
+  bool m_Refractive, m_Reflective;
 
   // Working variables; FIXME: not multithreading friendly;
   static thread_local TVector3 m_ImpactPoint, m_IncomingDirection, m_OutgoingDirection; //!

include/ParametricSurface.h

@@ -13,19 +13,19 @@ class ParametricSurface: public TObject {
  ParametricSurface(const TVector3 &x0, double umin, double umax, double vmin, double vmax):
   m_Center(x0), m_Umin(umin), m_Umax(umax), m_Vmin(vmin), m_Vmax(vmax) {};
   ~ParametricSurface() {};
-
+  
   // 2D parameter ranges (call them U&V);
   double Umin( void ) const { return m_Umin; };
   double Umax( void ) const { return m_Umax; };
   double Vmin( void ) const { return m_Vmin; };
   double Vmax( void ) const { return m_Vmax; };
-  bool IsInside(double u, double v) const {
+
+  virtual bool IsInside(double u, double v) const {
     return (u >= m_Umin && u <= m_Umax && v >= m_Vmin && v <= m_Vmax);
-  };
+  };  
   void SetUVranges(double umin, double umax, double vmin, double vmax) {
     m_Umin = umin; m_Umax = umax; m_Vmin = vmin; m_Vmax = vmax;
   };
-
   virtual TVector3 GetSpacePoint(double u, double v) const = 0;
   // There is no check that the point actually belongs to the surface;
   // it is assumed that GEANT stepping was done correctly, so the point 
@@ -75,7 +75,18 @@ class SphericalSurface: public ParametricSurface {
 		  double vmin = 0.0, double vmax = 2*M_PI): 
   ParametricSurface(x0, umin, umax, vmin, vmax), m_Concave(true), m_Radius(r0) {};
   ~SphericalSurface() {};
-
+  using ParametricSurface::IsInside;
+  virtual bool IsInside(TVector3 p) const {
+    double u = p.Theta();
+    double v = p.Phi();
+    if (v < 0) v += 2*M_PI;
+    // override, accounting for \phi periodicity
+    if (Vmin() <= Vmax()) {
+      return (u >= Umin() && u <= Umax() && v >= Vmin() && v <= Vmax());
+    } else {
+      return ( (v >= Vmin() || v <= Vmax()) && (u >= Umin() && u <= Umax() ) );
+    }    
+  }
   // FIXME: no range check?; 
   TVector3 GetSpacePoint(double theta, double phi) const {
     TVector3 nn(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta));
@@ -108,6 +119,32 @@ class SphericalSurface: public ParametricSurface {
 #endif
 };
 
+// want the same behavior as a SphericalSurface, but also checking
+// that any crossing point is also passing the given half-space condition
+class SphericalSurfaceWithHalfSpace: public SphericalSurface {
+public:
+  SphericalSurfaceWithHalfSpace(): SphericalSurface() {} ;
+  SphericalSurfaceWithHalfSpace(const TVector3 &x0, double r0, const std::vector<std::pair<TVector3, TVector3>> &halfSpaces,
+				double umin = 0.0, double umax = M_PI, 
+				double vmin = 0.0, double vmax = 2*M_PI): 
+    SphericalSurface(x0, r0, umin, umax, vmin, vmax){
+    for (const auto &halfSpace : halfSpaces) {
+      m_HSNorm.push_back(halfSpace.first);
+      m_HSPoint.push_back(halfSpace.second);
+    }
+  }
+  ~SphericalSurfaceWithHalfSpace() {};
+
+  bool IsInside(TVector3 p) const;
+private:
+  // point and normal defining half-space
+  std::vector<TVector3> m_HSNorm;
+  std::vector<TVector3> m_HSPoint;
+#ifndef DISABLE_ROOT_IO
+  ClassDef(SphericalSurfaceWithHalfSpace, 1);
+#endif
+};
+
 class LocalCoordinatesXY: public TObject {
  public:
   LocalCoordinatesXY() {};

include/irtLinkDef.h

@@ -30,6 +30,7 @@
 
 #pragma link C++ class ParametricSurface+;
 #pragma link C++ class SphericalSurface+;
+#pragma link C++ class SphericalSurfaceWithHalfSpace+;
 #pragma link C++ class LocalCoordinatesXY+;
 #pragma link C++ class FlatSurface+;
 
@@ -41,6 +42,7 @@
 #pragma link C++ class CherenkovMirrorGroup+;
 #pragma link C++ class FlatMirror+;
 #pragma link C++ class SphericalMirror+;
+#pragma link C++ class SphericalMirrorWithHalfSpace+;
 
 #pragma link C++ class CherenkovPhotonDetector+;
 

source/IRT.cc

@@ -9,22 +9,28 @@ thread_local TVector3 OpticalBoundary::m_OutgoingDirection;
 bool IRT::Transport(const TVector3 &xfrom, const TVector3 &nfrom)
 {
   bool transport_in_progress = false;
+  bool reflected = false;
   TVector3 x0 = xfrom, n0 = nfrom;
   // Just go through the optical boundaries, and calculate either reflection 
   // or refraction on that particular surface;
   for(unsigned iq=0; iq<_m_OpticalBoundaries.size(); iq++) {
     auto boundary = GetOpticalBoundary(iq), prev = iq ? GetOpticalBoundary(iq-1) : 0;
     auto surface = boundary->m_Surface;
 
+    // if already reflected, skip following mirrors
+    if(reflected && boundary->m_Reflective) continue;
+
     bool ok = surface->GetCrossing(x0, n0, &boundary->m_ImpactPoint);
-
+    
     // The logic here is that the first few boundaries may be irrelenat for this 
     // emission point (say for the gas case the emission point is beyond the aerogel-gas
     // refractive boundary; then just skip to the next one without changing [x0,n0]); 
     // however if one valid boundary was handled already, this must be a no-crossing case; 
     // then return false immediately;
     if (!ok) {
-      if (transport_in_progress) 
+      // if a mirror was missed, can still hit the next mirror.
+      // if missed a non-mirror boundary, return false
+      if (transport_in_progress && !boundary->m_Reflective) 
 	return false;
       else
 	continue;
@@ -37,7 +43,7 @@ bool IRT::Transport(const TVector3 &xfrom, const TVector3 &nfrom)
 
     boundary->m_OutgoingDirection = boundary->m_IncomingDirection;
     // Must be the sensor dump; FIXME:: do this check better;
-    if (!boundary->m_Radiator.GetObject()) return true;
+    if (!boundary->m_Radiator.GetObject() && !boundary->m_Reflective) return true;
 
     if (boundary->m_Refractive) {
       // Will not be able to determine the refractive index;
@@ -55,6 +61,7 @@ bool IRT::Transport(const TVector3 &xfrom, const TVector3 &nfrom)
       } //if
     } else {
       // Reflection;
+      reflected = true;
       boundary->m_OutgoingDirection.Rotate(M_PI - 2*acos(ns.Dot(boundary->m_IncomingDirection)), na);
     } //if
 

source/SphericalSurface.cc

@@ -21,7 +21,12 @@ bool SphericalSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVect
 
       if (t < 0.0) continue;
       *crs = x0 + t*n0;
-      if (check_normal && n0.Dot(GetNormal(*crs)) >= 0.0) continue;
+
+      // want to check if this is inside the theta/phi acceptance of a mirror
+      bool insidecheck = IsInside(*crs);
+
+      if(!(insidecheck)) continue;
+      if ( check_normal && n0.Dot(GetNormal(*crs)) >= 0.0) continue;
 
       return true;
     } //for iq

source/SphericalSurfaceWithHalfSpace.cc

@@ -0,0 +1,26 @@
+
+#include "ParametricSurface.h"
+// -------------------------------------------------------------------------------------
+
+
+// check if inside sphere and half space
+bool SphericalSurfaceWithHalfSpace::IsInside(TVector3 p) const {
+  // lab frame theta, phi
+  double u = p.Theta();
+  double v = p.Phi();
+  if (v < 0) v += 2*M_PI;
+
+  bool insideSphere;
+  if (Vmin() <= Vmax()) {
+    insideSphere = (u >= Umin() && u <= Umax() && v >= Vmin() && v <= Vmax());
+  } else {
+    insideSphere = ( (v >= Vmin() || v <= Vmax()) && (u >= Umin() && u <= Umax() ) );
+  }
+
+  // check if outside of any given half spaces
+  for(int i = 0; i < m_HSNorm.size(); i++){
+    if (m_HSNorm[i].Dot(p-m_HSPoint[i]) < 0) return false;
+  }
+
+  return insideSphere;
+} // SphericalSurfaceWithHalfSpace::IsInside()