diff --git a/README.md b/README.md index c5a84a60..b662017c 100644 --- a/README.md +++ b/README.md @@ -65,6 +65,7 @@ It is necessary to enable the `REST_MPFR` cmake option in order to benefit from ### Publications This repository makes use of the following published codes: + - K. Altenmuller et al, REST-for-Physics, a ROOT-based framework for event oriented data analysis and combined Monte Carlo response, [Computer Physics Communications 273, April 2022, 108281](https://doi.org/10.1016/j.cpc.2021.108281). - S.Hoof, J.Jaeckel, T.J.Lennert, Quantifying uncertainties in the solar axion flux and their impact on determining axion model parameters, [JCAP09(2021)006](https://doi.org/10.1088/1475-7516/2021/09/006). - T.Kittelmann, E.Klinkby, E.B.Knudsen, P.Willendrup, X.X.Cai, K.Kanaki, Monte Carlo Particle Lists: MCPL, Computer Physics Communications 218 (2017) 17–42](http://dx.doi.org/10.17632/cby92vsv5g.1). diff --git a/inc/TRestAxionTrueWolterOptics.h b/inc/TRestAxionTrueWolterOptics.h index b8ba4ba0..5dece507 100644 --- a/inc/TRestAxionTrueWolterOptics.h +++ b/inc/TRestAxionTrueWolterOptics.h @@ -35,9 +35,15 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics { /// Entrance radius R1 in mm. See schematic figure. std::vector fR1; //! + /// Radius R2 in mm. See schematic figure. + std::vector fR2; //! + /// Radius R3 in mm. See schematic figure. std::vector fR3; //! + /// Radius R4 in mm. See schematic figure. + std::vector fR4; //! + /// Radius R5 in mm. See schematic figure. std::vector fR5; //! @@ -47,6 +53,11 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics { /// Mirror thickness in mm. See schematic figure. std::vector fThickness; //! + /// Distance between mirror stacks in mm. See schematic figure. + /// This is here calculated using the functions from + /// https://backend.orbit.dtu.dk/ws/portalfiles/portal/122353510/phdthesis_for_DTU_orbit.pdf. + std::vector fXSep; //! + /// The spider structure to be used as an optical opaque mask (common to all planes) TRestSpiderMask* fSpiderMask = nullptr; //< @@ -83,12 +94,24 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics { return r; } + /// It returns a vector with the values of R2 + std::vector GetR2() { + std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 1); + return r; + } + /// It returns a vector with the values of R3 std::vector GetR3() { std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 2); return r; } + /// It returns a vector with the values of R4 + std::vector GetR4() { + std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 3); + return r; + } + /// It returns a vector with the values of R5 std::vector GetR5() { std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 4); @@ -110,13 +133,13 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics { /// It returns the entrance Z-position defined by the optical axis. Double_t GetEntrancePositionZ() override { - if (fCosAlpha.size() > 0) return -fMirrorLength * fCosAlpha[0]; + if (fCosAlpha.size() > 0) return -fMirrorLength * fCosAlpha[0] - 0.5 * fXSep[0]; return 0; } /// It returns the exit Z-position defined by the optical axis Double_t GetExitPositionZ() override { - if (fCosAlpha_3.size() > 0) return fMirrorLength * fCosAlpha_3[0]; + if (fCosAlpha_3.size() > 0) return fMirrorLength * fCosAlpha_3[0] + 0.5 * fXSep[0]; return 0; } diff --git a/inc/TRestAxionWolterOptics.h b/inc/TRestAxionWolterOptics.h index eb1e837a..103f29b9 100644 --- a/inc/TRestAxionWolterOptics.h +++ b/inc/TRestAxionWolterOptics.h @@ -35,9 +35,15 @@ class TRestAxionWolterOptics : public TRestAxionOptics { /// Entrance radius R1 in mm. See schematic figure. std::vector fR1; //! + /// Radius R2 in mm. See schematic figure. + std::vector fR2; //! + /// Radius R3 in mm. See schematic figure. std::vector fR3; //! + /// Radius R4 in mm. See schematic figure. + std::vector fR4; //! + /// Radius R5 in mm. See schematic figure. std::vector fR5; //! @@ -47,6 +53,11 @@ class TRestAxionWolterOptics : public TRestAxionOptics { /// Mirror thickness in mm. See schematic figure. std::vector fThickness; //! + /// Distance between mirror stacks in mm. See schematic figure. + /// This is here calculated using the functions from + /// https://backend.orbit.dtu.dk/ws/portalfiles/portal/122353510/phdthesis_for_DTU_orbit.pdf. + std::vector fXSep; //! + /// The spider structure to be used as an optical opaque mask (common to all planes) TRestSpiderMask* fSpiderMask = nullptr; //< @@ -83,12 +94,24 @@ class TRestAxionWolterOptics : public TRestAxionOptics { return r; } + /// It returns a vector with the values of R2 + std::vector GetR2() { + std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 1); + return r; + } + /// It returns a vector with the values of R3 std::vector GetR3() { std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 2); return r; } + /// It returns a vector with the values of R4 + std::vector GetR4() { + std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 3); + return r; + } + /// It returns a vector with the values of R5 std::vector GetR5() { std::vector r = TRestTools::GetColumnFromTable(fOpticsData, 4); @@ -110,13 +133,13 @@ class TRestAxionWolterOptics : public TRestAxionOptics { /// It returns the entrance Z-position defined by the optical axis. Double_t GetEntrancePositionZ() override { - if (fCosAlpha.size() > 0) return -fMirrorLength * fCosAlpha[0]; + if (fCosAlpha.size() > 0) return -fMirrorLength * fCosAlpha[0] - 0.5 * fXSep[0]; return 0; } /// It returns the exit Z-position defined by the optical axis Double_t GetExitPositionZ() override { - if (fCosAlpha_3.size() > 0) return fMirrorLength * fCosAlpha_3[0]; + if (fCosAlpha_3.size() > 0) return fMirrorLength * fCosAlpha_3[0] + 0.5 * fXSep[0]; return 0; } diff --git a/src/TRestAxionTrueWolterOptics.cxx b/src/TRestAxionTrueWolterOptics.cxx index 592c5ba0..3daa32e0 100644 --- a/src/TRestAxionTrueWolterOptics.cxx +++ b/src/TRestAxionTrueWolterOptics.cxx @@ -169,11 +169,22 @@ void TRestAxionTrueWolterOptics::Initialize() { SetLibraryVersion(LIBRARY_VERSION); fR1 = GetR1(); + fR2 = GetR2(); fR3 = GetR3(); + fR4 = GetR4(); fR5 = GetR5(); fAlpha = GetAlpha(); fThickness = GetThickness(); + fXSep.clear(); + for (unsigned int n = 0; n < fAlpha.size(); n++) + if (fR2[n] == fR3[n]) { + fXSep.push_back(0); + } else { + fXSep.push_back(2 * (fR1[n] - fR3[n] - fMirrorLength * TMath::Sin(fAlpha[n])) / + TMath::Tan(fAlpha[n])); + } + if (fAlpha.size() == 0) return; fCosAlpha.clear(); @@ -275,10 +286,11 @@ Int_t TRestAxionTrueWolterOptics::FirstMirrorReflection(const TVector3& pos, con TVector3 vertex(0, 0, fFrontVertex[mirror]); //// Reflection on first mirror - fFirstInteractionPosition = REST_Physics::GetParabolicVectorIntersection( - pos, dir, fAlpha[mirror], fR3[mirror], fMirrorLength); // should add this: TVector3(0, 0, -1), vertex + fFirstInteractionPosition = + REST_Physics::GetParabolicVectorIntersection(pos, dir, fAlpha[mirror], fR3[mirror]); - if (fFirstInteractionPosition.Z() < GetEntrancePositionZ() || fFirstInteractionPosition.Z() > 0) { + if (fFirstInteractionPosition.Z() < GetEntrancePositionZ() || + fFirstInteractionPosition.Z() > -(0.5 * fXSep[mirror])) { RESTDebug << "TRestAxionTrueWolterOptics::FirstMirrorReflection. No interaction!" << RESTendl; fFirstInteractionPosition = REST_Physics::MoveByDistance(pos, dir, fMirrorLength / 2.); fMiddleDirection = fEntranceDirection; @@ -320,13 +332,14 @@ Int_t TRestAxionTrueWolterOptics::SecondMirrorReflection(const TVector3& pos, co TVector3 vertex(0, 0, fBackVertex[mirror]); Double_t focal = fR3[mirror] / TMath::Tan(4 * fAlpha[mirror]); + Double_t beta = 3 * fAlpha[mirror]; - //// Reflection on first mirror + //// Reflection on second mirror fSecondInteractionPosition = - REST_Physics::GetHyperbolicVectorIntersection(pos, dir, fAlpha[mirror], fR3[mirror], fMirrorLength, - focal); // should add this: TVector3(0, 0, -1), vertex, + REST_Physics::GetHyperbolicVectorIntersection(pos, dir, beta, fR3[mirror], focal); - if (fSecondInteractionPosition.Z() > GetExitPositionZ() || fSecondInteractionPosition.Z() < 0) { + if (fSecondInteractionPosition.Z() > GetExitPositionZ() || + fSecondInteractionPosition.Z() < (0.5 * fXSep[mirror])) { RESTDebug << "TRestAxionTrueWolterOptics::SecondMirrorReflection. No interaction!" << RESTendl; fSecondInteractionPosition = REST_Physics::MoveByDistance(pos, dir, fMirrorLength / 2.); fExitDirection = fMiddleDirection; @@ -335,7 +348,7 @@ Int_t TRestAxionTrueWolterOptics::SecondMirrorReflection(const TVector3& pos, co } TVector3 hyperNormal = - REST_Physics::GetHyperbolicNormal(fSecondInteractionPosition, fAlpha[mirror], fR3[mirror], focal); + REST_Physics::GetHyperbolicNormal(fSecondInteractionPosition, beta, fR3[mirror], focal); fExitDirection = GetVectorReflection(fMiddleDirection, hyperNormal); diff --git a/src/TRestAxionWolterOptics.cxx b/src/TRestAxionWolterOptics.cxx index faac2445..27e84d81 100644 --- a/src/TRestAxionWolterOptics.cxx +++ b/src/TRestAxionWolterOptics.cxx @@ -169,11 +169,22 @@ void TRestAxionWolterOptics::Initialize() { SetLibraryVersion(LIBRARY_VERSION); fR1 = GetR1(); + fR2 = GetR2(); fR3 = GetR3(); + fR4 = GetR4(); fR5 = GetR5(); fAlpha = GetAlpha(); fThickness = GetThickness(); + fXSep.clear(); + for (unsigned int n = 0; n < fAlpha.size(); n++) + if (fR2[n] == fR3[n]) { + fXSep.push_back(0); + } else { + fXSep.push_back(2 * (fR1[n] - fR3[n] - fMirrorLength * TMath::Sin(fAlpha[n])) / + TMath::Tan(fAlpha[n])); + } + if (fAlpha.size() == 0) return; fCosAlpha.clear(); @@ -278,7 +289,8 @@ Int_t TRestAxionWolterOptics::FirstMirrorReflection(const TVector3& pos, const T fFirstInteractionPosition = pos + dir * REST_Physics::GetConeVectorIntersection(pos, dir, TVector3(0, 0, -1), vertex, cosA); - if (fFirstInteractionPosition.Z() < GetEntrancePositionZ() || fFirstInteractionPosition.Z() > 0) { + if (fFirstInteractionPosition.Z() < GetEntrancePositionZ() || + fFirstInteractionPosition.Z() > -(0.5 * fXSep[mirror])) { RESTDebug << "TRestAxionWolterOptics::FirstMirrorReflection. No interaction!" << RESTendl; fFirstInteractionPosition = REST_Physics::MoveByDistance(pos, dir, fMirrorLength / 2.); fMiddleDirection = fEntranceDirection; @@ -319,11 +331,12 @@ Int_t TRestAxionWolterOptics::SecondMirrorReflection(const TVector3& pos, const TVector3 vertex(0, 0, fBackVertex[mirror]); Double_t cosA = fCosAlpha_3[mirror]; - //// Reflection on first mirror + //// Reflection on second mirror fSecondInteractionPosition = pos + dir * REST_Physics::GetConeVectorIntersection(pos, dir, TVector3(0, 0, -1), vertex, cosA); - if (fSecondInteractionPosition.Z() > GetExitPositionZ() || fSecondInteractionPosition.Z() < 0) { + if (fSecondInteractionPosition.Z() > GetExitPositionZ() || + fSecondInteractionPosition.Z() < (0.5 * fXSep[mirror])) { RESTDebug << "TRestAxionWolterOptics::SecondMirrorReflection. No interaction!" << RESTendl; fSecondInteractionPosition = REST_Physics::MoveByDistance(pos, dir, fMirrorLength / 2.); fExitDirection = fMiddleDirection;