From 55c1966bc90d6cc730d87e34885aa682c3b68a1e Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 15 May 2024 08:56:02 -0400 Subject: [PATCH 01/12] For spherical optical boundary, require that the projected hit position in (theta(sphere frame), phi(lab frame)) is within supplied ( [umin,umax],[vmin,vmax] ) --- source/SphericalSurface.cc | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/source/SphericalSurface.cc b/source/SphericalSurface.cc index 1ec8ad3..5ac9743 100644 --- a/source/SphericalSurface.cc +++ b/source/SphericalSurface.cc @@ -21,7 +21,21 @@ 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; + + TVector3 crsSphereFrame = x0 + t*n0 - GetCenter(); + + double theta = crsSphereFrame.Theta(); + double phi = crsSphereFrame.Phi(); + double thetaLab = crs.Theta(); + double phiLab = crs.Phi(); + + // want to check if this is inside the theta/phi acceptance of a mirror + // \theta cut (in epic geometry) is applied in the sphere's reference frame, + // but \phi cut is applied in ePIC reference frame + bool insidecheck = IsInside(theta, philab); + + if(!(insidecheck)) continue; + if ( check_normal && n0.Dot(GetNormal(*crs)) >= 0.0) continue; return true; } //for iq From 65d6833bd3418b970b586d562139287e33eab5f5 Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 15 May 2024 09:10:06 -0400 Subject: [PATCH 02/12] Fix typo, pointer for position... --- source/SphericalSurface.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/SphericalSurface.cc b/source/SphericalSurface.cc index 5ac9743..c6adcb7 100644 --- a/source/SphericalSurface.cc +++ b/source/SphericalSurface.cc @@ -26,13 +26,13 @@ bool SphericalSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVect double theta = crsSphereFrame.Theta(); double phi = crsSphereFrame.Phi(); - double thetaLab = crs.Theta(); - double phiLab = crs.Phi(); + double thetaLab = crs->Theta(); + double phiLab = crs->Phi(); // want to check if this is inside the theta/phi acceptance of a mirror // \theta cut (in epic geometry) is applied in the sphere's reference frame, // but \phi cut is applied in ePIC reference frame - bool insidecheck = IsInside(theta, philab); + bool insidecheck = IsInside(theta, phiLab); if(!(insidecheck)) continue; if ( check_normal && n0.Dot(GetNormal(*crs)) >= 0.0) continue; From cba7a1b7ff53b8f0129f4dca33d97e15d0968d89 Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 15 May 2024 09:49:37 -0400 Subject: [PATCH 03/12] Add check for if the photon has already been reflected in transport --- source/IRT.cc | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/source/IRT.cc b/source/IRT.cc index 3b90b2f..c573a76 100644 --- a/source/IRT.cc +++ b/source/IRT.cc @@ -8,7 +8,9 @@ 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; @@ -17,20 +19,32 @@ bool IRT::Transport(const TVector3 &xfrom, const TVector3 &nfrom) auto surface = boundary->m_Surface; 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; + + // Connor, multi-mirror PR: + // If we will only have cases of one reflection, but with different possible reflecting + // surfaces, this can be generalized by checking if we already had a reflection when + // the propagation to the next mirror surface fails. if (!ok) { - if (transport_in_progress) + // sensor crossing + if (transport_in_progress && !boundary->m_Radiator.GetObject()){ return false; - else + } + // already reflected and checking other mirror boundary + if(transport_in_progress && reflected && !boundary->m_Refractive){ continue; - } //if - transport_in_progress = true; + } + else{ + continue; + } + } //if + transport_in_progress = true; + boundary->m_IncomingDirection = (boundary->m_ImpactPoint - x0).Unit(); TVector3 ns = surface->GetNormal(boundary->m_ImpactPoint); TVector3 na = ns.Cross(boundary->m_IncomingDirection); @@ -55,6 +69,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 @@ -124,8 +139,7 @@ IRTSolution IRT::Solve(const TVector3 &xfrom, const TVector3 &nfrom, const doubl //printf("%10.4f\n", dist); if (dist < m_Precision) { - solution.m_Converged = true; - + solution.m_Converged = true; { double slope = acos(nfrom.Dot(beam)); auto axis = nfrom.Cross(beam).Unit(); From 5321c78edcb87c69c188de6dd0fba707278b388d Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 15 May 2024 16:51:51 -0400 Subject: [PATCH 04/12] Account for \phi perodicity when checking if point is on spherical patch; Make 'IsInside' virtual, and redefine for spherical surface --- include/ParametricSurface.h | 16 ++++++++++++---- source/SphericalSurface.cc | 5 +++-- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/include/ParametricSurface.h b/include/ParametricSurface.h index f1375ba..c74f3e1 100644 --- a/include/ParametricSurface.h +++ b/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 @@ -76,6 +76,14 @@ class SphericalSurface: public ParametricSurface { ParametricSurface(x0, umin, umax, vmin, vmax), m_Concave(true), m_Radius(r0) {}; ~SphericalSurface() {}; + bool IsInside(double u, double v) const { + // 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)); diff --git a/source/SphericalSurface.cc b/source/SphericalSurface.cc index c6adcb7..c681276 100644 --- a/source/SphericalSurface.cc +++ b/source/SphericalSurface.cc @@ -25,9 +25,10 @@ bool SphericalSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVect TVector3 crsSphereFrame = x0 + t*n0 - GetCenter(); double theta = crsSphereFrame.Theta(); - double phi = crsSphereFrame.Phi(); - double thetaLab = crs->Theta(); + //double phi = crsSphereFrame.Phi(); + //double thetaLab = crs->Theta(); double phiLab = crs->Phi(); + if (phiLab < 0) phiLab += 2*M_PI; // want to check if this is inside the theta/phi acceptance of a mirror // \theta cut (in epic geometry) is applied in the sphere's reference frame, From f76ef980d7350394a15220a44b7fb04ddaaf603c Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 3 Jul 2024 10:30:11 -0400 Subject: [PATCH 05/12] Call IsInside for SphericalSurface with only the given crossing point --- source/SphericalSurface.cc | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/source/SphericalSurface.cc b/source/SphericalSurface.cc index c681276..186b7e5 100644 --- a/source/SphericalSurface.cc +++ b/source/SphericalSurface.cc @@ -22,18 +22,8 @@ bool SphericalSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVect if (t < 0.0) continue; *crs = x0 + t*n0; - TVector3 crsSphereFrame = x0 + t*n0 - GetCenter(); - - double theta = crsSphereFrame.Theta(); - //double phi = crsSphereFrame.Phi(); - //double thetaLab = crs->Theta(); - double phiLab = crs->Phi(); - if (phiLab < 0) phiLab += 2*M_PI; - // want to check if this is inside the theta/phi acceptance of a mirror - // \theta cut (in epic geometry) is applied in the sphere's reference frame, - // but \phi cut is applied in ePIC reference frame - bool insidecheck = IsInside(theta, phiLab); + bool insidecheck = IsInside(*crs); if(!(insidecheck)) continue; if ( check_normal && n0.Dot(GetNormal(*crs)) >= 0.0) continue; From 3a11ae43ea0dcc0a0b73a0d3924b0e117a8a3568 Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 3 Jul 2024 10:31:04 -0400 Subject: [PATCH 06/12] Add SphericalSurfaceWithHalfSpace (spherical surface with requirement that crossing points also fall in given half-space) --- include/ParametricSurface.h | 27 +++++++++++++++++++++++-- source/SphericalSurfaceWithHalfSpace.cc | 21 +++++++++++++++++++ 2 files changed, 46 insertions(+), 2 deletions(-) create mode 100644 source/SphericalSurfaceWithHalfSpace.cc diff --git a/include/ParametricSurface.h b/include/ParametricSurface.h index c74f3e1..a2fd7ac 100644 --- a/include/ParametricSurface.h +++ b/include/ParametricSurface.h @@ -75,8 +75,11 @@ 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() {}; - - bool IsInside(double u, double v) const { + 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()); @@ -116,6 +119,26 @@ 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 TVector3 &pNorm, const TVector3 &pPoint, double umin = 0.0, double umax = M_PI, + double vmin = 0.0, double vmax = 2*M_PI): + SphericalSurface(x0, r0, umin, umax, vmin, vmax), m_HSNorm(pNorm), m_HSPoint(pPoint) {}; + ~SphericalSurfaceWithHalfSpace() {}; + + bool IsInside(TVector3 p) const; +private: + // point and normal defining half-space + TVector3 m_HSNorm; + TVector3 m_HSPoint; +#ifndef DISABLE_ROOT_IO + ClassDef(SphericalSurfaceWithHalfSpace, 1); +#endif +}; + class LocalCoordinatesXY: public TObject { public: LocalCoordinatesXY() {}; diff --git a/source/SphericalSurfaceWithHalfSpace.cc b/source/SphericalSurfaceWithHalfSpace.cc new file mode 100644 index 0000000..bc611f3 --- /dev/null +++ b/source/SphericalSurfaceWithHalfSpace.cc @@ -0,0 +1,21 @@ + +#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(); + + bool insideSphere; + if (Vmin() <= Vmax()) { + insideSphere = (u >= Umin() && u <= Umax() && v >= Vmin() && v <= Vmax()); + } else { + insideSphere = ( (v >= Vmin() || v <= Vmax()) && (u >= Umin() && u <= Umax() ) ); + } + + auto insideHalfSpace = (m_HSNorm.Dot(p-m_HSPoint) > 0); + return (insideSphere && insideHalfSpace); +} // SphericalSurfaceWithHalfSpace::IsInside() From ab0ebac2800087d59a98ad3f4d88383b095a290c Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 3 Jul 2024 10:31:51 -0400 Subject: [PATCH 07/12] Add SphericalSurfaceWithHalfSpace to build --- CMakeLists.txt | 1 + include/irtLinkDef.h | 2 ++ 2 files changed, 3 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d5d420..3af8373 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/irtLinkDef.h b/include/irtLinkDef.h index 7e3ba84..4007778 100644 --- a/include/irtLinkDef.h +++ b/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+; From 6cd9c635d0242b6f2e82d2cfbf6a1360386b3fc5 Mon Sep 17 00:00:00 2001 From: cpecar Date: Wed, 3 Jul 2024 10:32:17 -0400 Subject: [PATCH 08/12] Add SphericalMirrorWithHalfSpace --- include/CherenkovMirror.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/CherenkovMirror.h b/include/CherenkovMirror.h index b585626..cb24905 100644 --- a/include/CherenkovMirror.h +++ b/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 TVector3 &pNorm, const TVector3 &pPoint): + CherenkovMirror(solid, material), SphericalSurfaceWithHalfSpace(x0, r0, pNorm, pPoint) {}; + ~SphericalMirrorWithHalfSpace() {}; + +#ifndef DISABLE_ROOT_IO + ClassDef(SphericalMirrorWithHalfSpace, 2); +#endif +}; + class CherenkovMirrorGroup: public TObject { public: CherenkovMirrorGroup() {}; From 782500a968cde5206e2c1e776680d7ea8a25cb71 Mon Sep 17 00:00:00 2001 From: cpecar Date: Fri, 12 Jul 2024 13:20:28 -0400 Subject: [PATCH 09/12] fix phi range for sphere IsInside check --- source/SphericalSurfaceWithHalfSpace.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/source/SphericalSurfaceWithHalfSpace.cc b/source/SphericalSurfaceWithHalfSpace.cc index bc611f3..1ec2317 100644 --- a/source/SphericalSurfaceWithHalfSpace.cc +++ b/source/SphericalSurfaceWithHalfSpace.cc @@ -8,14 +8,15 @@ 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() ) ); - } - + } auto insideHalfSpace = (m_HSNorm.Dot(p-m_HSPoint) > 0); + return (insideSphere && insideHalfSpace); } // SphericalSurfaceWithHalfSpace::IsInside() From 5b38bfc49b7cd69ff020a20370e1814f0f9563a7 Mon Sep 17 00:00:00 2001 From: cpecar Date: Fri, 12 Jul 2024 13:21:51 -0400 Subject: [PATCH 10/12] Add m_Reflective property to OpticalBoundary, to differentiate mirrors and sensor surface --- include/CherenkovDetectorCollection.h | 6 ++--- include/CherenkovPhotonDetector.h | 2 +- include/OpticalBoundary.h | 8 +++---- source/IRT.cc | 33 +++++++++++---------------- 4 files changed, 21 insertions(+), 28 deletions(-) diff --git a/include/CherenkovDetectorCollection.h b/include/CherenkovDetectorCollection.h index 797107f..01a7560 100644 --- a/include/CherenkovDetectorCollection.h +++ b/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); diff --git a/include/CherenkovPhotonDetector.h b/include/CherenkovPhotonDetector.h index 9394078..1fbdb99 100644 --- a/include/CherenkovPhotonDetector.h +++ b/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); diff --git a/include/OpticalBoundary.h b/include/OpticalBoundary.h index bad9022..ca98702 100644 --- a/include/OpticalBoundary.h +++ b/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; //! diff --git a/source/IRT.cc b/source/IRT.cc index c573a76..93e3e4a 100644 --- a/source/IRT.cc +++ b/source/IRT.cc @@ -8,7 +8,6 @@ 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; @@ -18,40 +17,33 @@ bool IRT::Transport(const TVector3 &xfrom, const TVector3 &nfrom) 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); + printf("boundary %d, isOK %d \n", iq, ok); // 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; - - // Connor, multi-mirror PR: - // If we will only have cases of one reflection, but with different possible reflecting - // surfaces, this can be generalized by checking if we already had a reflection when - // the propagation to the next mirror surface fails. if (!ok) { - // sensor crossing - if (transport_in_progress && !boundary->m_Radiator.GetObject()){ + // 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; - } - // already reflected and checking other mirror boundary - if(transport_in_progress && reflected && !boundary->m_Refractive){ - continue; - } - else{ + else continue; - } - } //if - + } //if transport_in_progress = true; - + boundary->m_IncomingDirection = (boundary->m_ImpactPoint - x0).Unit(); TVector3 ns = surface->GetNormal(boundary->m_ImpactPoint); TVector3 na = ns.Cross(boundary->m_IncomingDirection); 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; @@ -139,7 +131,8 @@ IRTSolution IRT::Solve(const TVector3 &xfrom, const TVector3 &nfrom, const doubl //printf("%10.4f\n", dist); if (dist < m_Precision) { - solution.m_Converged = true; + solution.m_Converged = true; + { double slope = acos(nfrom.Dot(beam)); auto axis = nfrom.Cross(beam).Unit(); From 52bd3e11820ef1d738fb702e7f37d98aa20e86e0 Mon Sep 17 00:00:00 2001 From: cpecar Date: Fri, 2 Aug 2024 12:31:50 -0400 Subject: [PATCH 11/12] Supply list of half space cuts to apply, to allow for use of more than 2 mirrors --- include/CherenkovMirror.h | 4 ++-- include/ParametricSurface.h | 18 ++++++++++++------ source/SphericalSurfaceWithHalfSpace.cc | 10 +++++++--- 3 files changed, 21 insertions(+), 11 deletions(-) diff --git a/include/CherenkovMirror.h b/include/CherenkovMirror.h index cb24905..c35bd31 100644 --- a/include/CherenkovMirror.h +++ b/include/CherenkovMirror.h @@ -74,8 +74,8 @@ class SphericalMirror: public CherenkovMirror, public SphericalSurface { class SphericalMirrorWithHalfSpace: public CherenkovMirror, public SphericalSurfaceWithHalfSpace { public: SphericalMirrorWithHalfSpace() {}; - SphericalMirrorWithHalfSpace(G4VSolid *solid, G4Material *material, const TVector3 &x0, double r0, const TVector3 &pNorm, const TVector3 &pPoint): - CherenkovMirror(solid, material), SphericalSurfaceWithHalfSpace(x0, r0, pNorm, pPoint) {}; + SphericalMirrorWithHalfSpace(G4VSolid *solid, G4Material *material, const TVector3 &x0, double r0, const std::vector> halfSpaces): + CherenkovMirror(solid, material), SphericalSurfaceWithHalfSpace(x0, r0, halfSpaces) {}; ~SphericalMirrorWithHalfSpace() {}; #ifndef DISABLE_ROOT_IO diff --git a/include/ParametricSurface.h b/include/ParametricSurface.h index a2fd7ac..8a1ab61 100644 --- a/include/ParametricSurface.h +++ b/include/ParametricSurface.h @@ -124,16 +124,22 @@ class SphericalSurface: public ParametricSurface { class SphericalSurfaceWithHalfSpace: public SphericalSurface { public: SphericalSurfaceWithHalfSpace(): SphericalSurface() {} ; - SphericalSurfaceWithHalfSpace(const TVector3 &x0, double r0, const TVector3 &pNorm, const TVector3 &pPoint, double umin = 0.0, double umax = M_PI, - double vmin = 0.0, double vmax = 2*M_PI): - SphericalSurface(x0, r0, umin, umax, vmin, vmax), m_HSNorm(pNorm), m_HSPoint(pPoint) {}; + SphericalSurfaceWithHalfSpace(const TVector3 &x0, double r0, const std::vector> &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 - TVector3 m_HSNorm; - TVector3 m_HSPoint; + std::vector m_HSNorm; + std::vector m_HSPoint; #ifndef DISABLE_ROOT_IO ClassDef(SphericalSurfaceWithHalfSpace, 1); #endif diff --git a/source/SphericalSurfaceWithHalfSpace.cc b/source/SphericalSurfaceWithHalfSpace.cc index 1ec2317..0464fb1 100644 --- a/source/SphericalSurfaceWithHalfSpace.cc +++ b/source/SphericalSurfaceWithHalfSpace.cc @@ -15,8 +15,12 @@ bool SphericalSurfaceWithHalfSpace::IsInside(TVector3 p) const { insideSphere = (u >= Umin() && u <= Umax() && v >= Vmin() && v <= Vmax()); } else { insideSphere = ( (v >= Vmin() || v <= Vmax()) && (u >= Umin() && u <= Umax() ) ); - } - auto insideHalfSpace = (m_HSNorm.Dot(p-m_HSPoint) > 0); + } - return (insideSphere && insideHalfSpace); + // 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() From 5efe0e7eb5b4bf25b4e7401172275704dab802b7 Mon Sep 17 00:00:00 2001 From: cpecar Date: Fri, 23 Aug 2024 13:04:10 -0400 Subject: [PATCH 12/12] remove unnecessary prints --- source/IRT.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/IRT.cc b/source/IRT.cc index 93e3e4a..683ce39 100644 --- a/source/IRT.cc +++ b/source/IRT.cc @@ -21,7 +21,7 @@ bool IRT::Transport(const TVector3 &xfrom, const TVector3 &nfrom) if(reflected && boundary->m_Reflective) continue; bool ok = surface->GetCrossing(x0, n0, &boundary->m_ImpactPoint); - printf("boundary %d, isOK %d \n", iq, ok); + // 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]);