From 51b64e0c9785bc0597f9ce470dbc9ad001ee1785 Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Thu, 23 Sep 2021 11:42:14 -0400 Subject: [PATCH] Fixed energy loss bugs, including case of backwards travel through a non-uniform target. Modified premake file to be more universally compatible for ROOT. Requires user specification of ROOT directories in the premake file. --- etc/convertedAnasenCoords.txt | 30 ++++ include/AnasenEfficiency.h | 23 ++- include/LayeredTarget.h | 3 + include/QQQDetector.h | 23 ++- include/StripDetector.h | 23 ++- include/Target.h | 8 +- input.txt | 13 +- premake5.lua | 47 +++-- src/Detectors/AnasenEfficiency.cpp | 260 +++++++++++++++++++++++++-- src/Detectors/DetectorEfficiency.cpp | 9 +- src/Detectors/QQQDetector.cpp | 18 +- src/Detectors/SabreEfficiency.cpp | 4 +- src/Detectors/StripDetector.cpp | 41 +++-- src/LayeredTarget.cpp | 68 ++++--- src/Target.cpp | 38 ++-- src/TwoStepSystem.cpp | 4 +- 16 files changed, 472 insertions(+), 140 deletions(-) create mode 100644 etc/convertedAnasenCoords.txt diff --git a/etc/convertedAnasenCoords.txt b/etc/convertedAnasenCoords.txt new file mode 100644 index 0000000..703ce91 --- /dev/null +++ b/etc/convertedAnasenCoords.txt @@ -0,0 +1,30 @@ +For SX3's: det id | mid rho | mid z | mid phi | name +For QQQ's: det id | mid phi | z | name +0 5.49779 0 QQQ_0 +1 0.785398 0.8028 QQQ_1 +2 2.35619 0 QQQ_2 +3 3.92699 0 QQQ_3 +6 8.90601 6.25 0.785795 R1a_A +7 8.89871 6.25 0.262014 R1a_B +8 8.90354 6.25 6.02132 R1a_C +9 8.90247 6.25 5.49779 R1a_D +10 8.90354 6.25 4.97426 R1a_E +11 8.90354 6.25 4.45052 R1a_F +12 8.90247 6.25 3.92699 R1b_A +13 8.90354 6.25 3.40346 R1b_B +14 8.90354 6.25 2.87972 R1b_C +15 8.90247 6.25 2.35619 R1b_D +16 8.90354 6.25 1.83266 R1b_E +17 8.90354 6.25 1.30893 R1b_F +18 8.90601 18.65 0.785795 R2a_A +19 8.89871 18.65 0.262014 R2a_B +20 8.90354 18.65 6.02132 R2a_C +21 8.90247 18.65 5.49779 R2a_D +22 8.90354 18.65 4.97426 R2a_E +23 8.90354 18.65 4.45052 R2a_F +24 8.90247 18.65 3.92699 R2b_A +25 8.90354 18.65 3.40346 R2b_B +26 8.90354 18.65 2.87972 R2b_C +27 8.90247 18.65 2.35619 R2b_D +28 8.90354 18.65 1.83266 R2b_E +29 8.90354 18.65 1.30893 R2b_F diff --git a/include/AnasenEfficiency.h b/include/AnasenEfficiency.h index fbe1ec8..d59f85f 100644 --- a/include/AnasenEfficiency.h +++ b/include/AnasenEfficiency.h @@ -6,6 +6,16 @@ #include "DetectorEfficiency.h" #include "StripDetector.h" #include "QQQDetector.h" +#include "Target.h" +#include "Nucleus.h" +#include "MaskFile.h" + +struct DetectorResult { + bool detectFlag = false; + Mask::Vec3 direction; + double energy_deposited = 0.0; + std::string det_name = ""; +}; class AnasenEfficiency : public DetectorEfficiency { public: @@ -16,14 +26,18 @@ public: double RunConsistencyCheck() override; private: - bool IsRing1(double theta, double phi); - bool IsRing2(double theta, double phi); - bool IsQQQ(double theta, double phi); + DetectorResult IsRing1(Mask::Nucleus& nucleus); + DetectorResult IsRing2(Mask::Nucleus& nucleus); + DetectorResult IsQQQ(Mask::Nucleus& nucleus); + DetectorResult IsAnasen(Mask::Nucleus& nucleus); + void CountCoincidences(const Mask::MaskFileData& data, std::vector& counts, int rxn_type); std::vector m_Ring1, m_Ring2; std::vector m_forwardQQQs; std::vector m_backwardQQQs; + Mask::Target det_silicon; + /**** ANASEN geometry constants *****/ const int n_sx3_per_ring = 12; const int n_qqq = 4; @@ -42,8 +56,9 @@ private: const double ring_phi[12] = {0.785795, 0.262014, 6.02132, 5.49779, 4.97426, 4.45052, 3.92699, 3.40346, 2.87972, 2.35619, 1.83266, 1.30893}; /*************************/ - static constexpr double threshold = 0.2; //MeV + static constexpr double threshold = 0.6; //MeV static constexpr double deg2rad = M_PI/180.0; + static constexpr double si_thickness = 1000 * 1e-4 * 2.3926 * 1e6; //thickness in um -> eff thickness in ug/cm^2 for detector }; #endif \ No newline at end of file diff --git a/include/LayeredTarget.h b/include/LayeredTarget.h index 16206e4..8087fbb 100644 --- a/include/LayeredTarget.h +++ b/include/LayeredTarget.h @@ -15,6 +15,8 @@ Written by G.W. McCann Aug. 2020 #include #include +#include +#include "RandomGenerator.h" #include "Target.h" namespace Mask { @@ -37,6 +39,7 @@ namespace Mask { private: std::vector layers; std::string name; + std::uniform_real_distribution m_fractional_depth_dist; }; } diff --git a/include/QQQDetector.h b/include/QQQDetector.h index 9ee6d90..204d8d0 100644 --- a/include/QQQDetector.h +++ b/include/QQQDetector.h @@ -10,36 +10,45 @@ #include #include +#include #include "Vec3.h" #include "Rotation.h" +#include "RandomGenerator.h" class QQQDetector { public: QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x=0, double y=0); ~QQQDetector(); - inline Mask::Vec3 GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; }; - inline Mask::Vec3 GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; }; + inline Mask::Vec3 GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; } + inline Mask::Vec3 GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; } + inline Mask::Vec3 GetNorm() { return m_norm; } Mask::Vec3 GetTrajectoryCoordinates(double theta, double phi); std::pair GetTrajectoryRingWedge(double theta, double phi); Mask::Vec3 GetHitCoordinates(int ringch, int wedgech); + inline void TurnOnRandomizedCoordinates() { rndmFlag = true; } + inline void TurnOffRandomizedCoordinates() { rndmFlag = false; } - inline int GetNumberOfRings() { return nrings; }; - inline int GetNumberOfWedges() { return nwedges; }; + inline int GetNumberOfRings() { return nrings; } + inline int GetNumberOfWedges() { return nwedges; } private: - inline bool CheckChannel(int ch) { return (ch >=0 && ch < nrings); }; - inline bool CheckCorner(int corner) { return (corner >=0 && corner < 4); }; + inline bool CheckChannel(int ch) { return (ch >=0 && ch < nrings); } + inline bool CheckCorner(int corner) { return (corner >=0 && corner < 4); } void CalculateCorners(); - Mask::Vec3 TransformCoordinates(Mask::Vec3& vector) { return m_ZRot*vector + m_translation; }; + Mask::Vec3 TransformCoordinates(Mask::Vec3& vector) { return m_ZRot*vector + m_translation; } double m_Rinner, m_Router, m_deltaR, m_deltaPhi, m_deltaPhi_per_wedge, m_phiCentral; std::vector> m_ringCoords, m_wedgeCoords; Mask::Vec3 m_translation; + Mask::Vec3 m_norm; Mask::ZRotation m_ZRot; + std::uniform_real_distribution m_uniform_fraction; + bool rndmFlag; + static constexpr int nrings = 16; static constexpr int nwedges = 16; static constexpr double deg2rad = M_PI/180.0; diff --git a/include/StripDetector.h b/include/StripDetector.h index 7479cf0..c397db0 100644 --- a/include/StripDetector.h +++ b/include/StripDetector.h @@ -17,21 +17,28 @@ #include "Vec3.h" #include "Rotation.h" +#include "RandomGenerator.h" class StripDetector { public: StripDetector(int ns, double len, double wid, double cphi, double cz, double crho); ~StripDetector(); - inline Mask::Vec3 GetFrontStripCoordinates(int stripch, int corner) { return front_strip_coords[stripch][corner]; }; - inline Mask::Vec3 GetBackStripCoordinates(int stripch, int corner) { return back_strip_coords[stripch][corner]; }; - inline Mask::Vec3 GetRotatedFrontStripCoordinates(int stripch, int corner) { return rotated_front_strip_coords[stripch][corner]; }; - inline Mask::Vec3 GetRotatedBackStripCoordinates(int stripch, int corner) { return rotated_back_strip_coords[stripch][corner]; }; - inline void SetRandomNumberGenerator(std::mt19937* random) { m_random = random; }; + inline Mask::Vec3 GetFrontStripCoordinates(int stripch, int corner) { return front_strip_coords[stripch][corner]; } + inline Mask::Vec3 GetBackStripCoordinates(int stripch, int corner) { return back_strip_coords[stripch][corner]; } + inline Mask::Vec3 GetRotatedFrontStripCoordinates(int stripch, int corner) { return rotated_front_strip_coords[stripch][corner]; } + inline Mask::Vec3 GetRotatedBackStripCoordinates(int stripch, int corner) { return rotated_back_strip_coords[stripch][corner]; } + inline Mask::Vec3 GetNormRotated() { return zRot*m_norm_unrot; } + + inline void TurnOnRandomizedCoordinates() { rndmFlag = true; } + inline void TurnOffRandomizedCoordinates() { rndmFlag = false; } + Mask::Vec3 GetHitCoordinates(int front_stripch, double front_strip_ratio); std::pair GetChannelRatio(double theta, double phi); private: + inline bool ValidChannel(int f) { return ((f >= 0 && f < num_strips) ? true : false); }; + inline bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); }; void CalculateCorners(); int num_strips; @@ -49,14 +56,14 @@ private: std::vector> front_strip_coords, back_strip_coords; std::vector> rotated_front_strip_coords, rotated_back_strip_coords; + Mask::Vec3 m_norm_unrot; + Mask::ZRotation zRot; std::uniform_real_distribution m_uniform_fraction; - std::mt19937* m_random; //Not owned by StripDetector! + bool rndmFlag; - inline bool ValidChannel(int f) { return ((f >= 0 && f < num_strips) ? true : false); }; - inline bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); }; }; #endif diff --git a/include/Target.h b/include/Target.h index 93c977b..bcaafd9 100644 --- a/include/Target.h +++ b/include/Target.h @@ -27,10 +27,10 @@ namespace Mask { ~Target(); void SetElements(std::vector& z, std::vector& a, std::vector& stoich); bool ContainsElement(int z, int a); - double getEnergyLossTotal(int zp, int ap, double startEnergy, double angle); - double getEnergyLossHalf(int zp, int ap, double startEnergy, double angle); - double getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle); - double getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double angle); + double GetEnergyLossTotal(int zp, int ap, double startEnergy, double angle); + double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle); + double GetEnergyLossFractionalDepth(int zp, int ap, double startEnergy, double angle, double percent_depth); + double GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double angle, double percent_depth); inline const double& GetThickness() { return thickness; } inline int GetNumberOfElements() { return Z.size(); } inline int GetElementZ(int index) { return Z[index]; } diff --git a/input.txt b/input.txt index 85bca36..97ca403 100644 --- a/input.txt +++ b/input.txt @@ -1,14 +1,13 @@ ----------Data Information---------- -OutputFile: 7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.mask +OutputFile: /data1/gwm17/mask_tests/7Be12C_870keV_beam_50CD2.mask SaveTree: yes SavePlots: yes ----------Reaction Information---------- -ReactionType: 2 +ReactionType: 1 Z A (order is target, projectile, ejectile, break1, break3(if pure decay is target, ejectile)) -1 2 +6 12 +4 7 4 7 -1 1 -2 4 ----------Target Information---------- Name: test_targ Layers: 1 @@ -20,8 +19,8 @@ Z A Stoich 0 ~ ----------Sampling Information---------- -NumberOfSamples: 10000 -BeamMeanEnergy(MeV): 0.6 BeamEnergySigma(MeV): 0.0 +NumberOfSamples: 100000 +BeamMeanEnergy(MeV): 0.87 BeamEnergySigma(MeV): 0.0 EjectileThetaType(0=Lab,1=CM): 1 EjectileThetaMin(deg): 0.0 EjectileThetaMax(deg): 180.0 EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 360.0 diff --git a/premake5.lua b/premake5.lua index 1eb4b76..cd4c287 100644 --- a/premake5.lua +++ b/premake5.lua @@ -43,36 +43,32 @@ project "RootPlot" "include/*.h" } - filter "system:windows" - systemversion "latest" - --ROOT cannot be located using the config tools, so we must query for the ROOTSYS env variable - --Have to use an if statement to hide this (@penguin for example doesn't have a ROOTSYS) - if os.host() == windows then - rootpath = os.getenv("ROOTSYS") + --User specified path to ROOT CERN libraries-- + ROOTIncludepath = "/usr/include/root" + ROOTLibpath = "/usr/lib64/root" - includedirs { - "include", - "./", - rootpath .. "include" - } + includedirs { + "include" + } - links { - rootpath .. "lib/**.lib" - } - end + sysincludedirs { + ROOTIncludepath + } + + libdirs { + ROOTLibpath + } + + links { + "Gui", "Core", "Imt", "RIO", "Net", "Hist", + "Graf", "Graf3d", "Gpad", "ROOTDataFrame", "ROOTVecOps", + "Tree", "TreePlayer", "Rint", "Postscript", "Matrix", + "Physics", "MathCore", "Thread", "MultiProc", "m", "dl" + } filter "system:macosx or linux" - includedirs { - "include", - "./" - } - - buildoptions { - "`root-config --cflags`" - } - linkoptions { - "`root-config --glibs`" + "-pthread" } filter "configurations:Debug" @@ -98,6 +94,7 @@ project "DetectEff" "src/Rotation.cpp", "src/Target.cpp", "src/EnergyLoss.cpp", + "src/RandomGenerator.cpp", "include/*.h" } diff --git a/src/Detectors/AnasenEfficiency.cpp b/src/Detectors/AnasenEfficiency.cpp index 117644a..560dfe7 100644 --- a/src/Detectors/AnasenEfficiency.cpp +++ b/src/Detectors/AnasenEfficiency.cpp @@ -1,20 +1,29 @@ #include "AnasenEfficiency.h" #include "Kinematics.h" -#include "MaskFile.h" #include #include AnasenEfficiency::AnasenEfficiency() : - DetectorEfficiency() + DetectorEfficiency(), det_silicon(si_thickness) { for(int i=0; i det_z = {14}; + std::vector det_a = {28}; + std::vector det_stoich = {1}; + det_silicon.SetElements(det_z, det_a, det_stoich); + } AnasenEfficiency::~AnasenEfficiency() {} @@ -196,47 +205,199 @@ double AnasenEfficiency::RunConsistencyCheck() { } -bool AnasenEfficiency::IsRing1(double theta, double phi) { +DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) { + + DetectorResult observation; + //Mask::Vec3 coords; + double thetaIncident, eloss, e_dep; for(auto& sx3 : m_Ring1) { - auto result = sx3.GetChannelRatio(theta, phi); + auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - return true; + //coords = sx3.GetHitCoordinates(result.first, result.second); + observation.detectFlag = true; + observation.direction = sx3.GetHitCoordinates(result.first, result.second); + thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR()); + if(thetaIncident > M_PI/2.0) + thetaIncident = M_PI - thetaIncident; + + //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + observation.det_name = "R1"; + //return std::make_pair(true, e_dep); + return observation; } } - return false; + //return std::make_pair(false, 0.0); + return observation; } -bool AnasenEfficiency::IsRing2(double theta, double phi) { +DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) { + + DetectorResult observation; + //Mask::Vec3 coords; + double thetaIncident, eloss, e_dep; for(auto& sx3 : m_Ring2) { - auto result = sx3.GetChannelRatio(theta, phi); + auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - return true; + //coords = sx3.GetHitCoordinates(result.first, result.second); + observation.detectFlag = true; + observation.direction = sx3.GetHitCoordinates(result.first, result.second); + thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR()); + if(thetaIncident > M_PI/2.0) + thetaIncident = M_PI - thetaIncident; + + //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + observation.det_name = "R2"; + //return std::make_pair(true, e_dep); + return observation; } } - return false; + return observation; } -bool AnasenEfficiency::IsQQQ(double theta, double phi) { +DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) { + + DetectorResult observation; + //Mask::Vec3 coords; + double thetaIncident, eloss, e_dep; + for(auto& qqq : m_forwardQQQs) { - auto result = qqq.GetTrajectoryRingWedge(theta, phi); + auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - return true; + //coords = qqq.GetHitCoordinates(result.first, result.second); + observation.detectFlag = true; + observation.direction = qqq.GetHitCoordinates(result.first, result.second); + thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR()); + if(thetaIncident > M_PI/2.0) + thetaIncident = M_PI - thetaIncident; + + //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + //return std::make_pair(true, e_dep); + observation.det_name = "FQQQ"; + return observation; } } for(auto& qqq : m_backwardQQQs) { - auto result = qqq.GetTrajectoryRingWedge(theta, phi); + auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - return true; + //coords = qqq.GetHitCoordinates(result.first, result.second); + observation.detectFlag = true; + observation.direction = qqq.GetHitCoordinates(result.first, result.second); + thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR()); + if(thetaIncident > M_PI/2.0) + thetaIncident = M_PI - thetaIncident; + + //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); + //return std::make_pair(true, e_dep); + observation.det_name = "BQQQ"; + return observation; } } + //return std::make_pair(false, 0.0); + return observation; +} +DetectorResult AnasenEfficiency::IsAnasen(Mask::Nucleus& nucleus) { + DetectorResult result; + if(nucleus.GetKE() <= threshold) + return result; - return false; + if(!result.detectFlag) { + result = IsRing1(nucleus); + } + if(!result.detectFlag) { + result = IsRing2(nucleus); + } + if(!result.detectFlag) { + result = IsQQQ(nucleus); + } + + return result; +} + +void AnasenEfficiency::CountCoincidences(const Mask::MaskFileData& data, std::vector& counts, int rxn_type) { + if (rxn_type == 0 && data.detect_flag[1] && data.detect_flag[2]) + { + counts[0]++; + } + else if (rxn_type == 1 && data.detect_flag[2] && data.detect_flag[3]) + { + counts[0]++; + } + else if(rxn_type == 2) + { + if(data.detect_flag[2] && data.detect_flag[4]) + { + counts[0]++; + } + if(data.detect_flag[2] && data.detect_flag[5]) + { + counts[1]++; + } + if(data.detect_flag[4] && data.detect_flag[5]) + { + counts[2]++; + } + if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[5]) + { + counts[3]++; + } + } + else if(rxn_type == 3) + { + if(data.detect_flag[2] && data.detect_flag[4]) + { + counts[0]++; + } + if(data.detect_flag[2] && data.detect_flag[6]) + { + counts[1]++; + } + if(data.detect_flag[2] && data.detect_flag[7]) + { + counts[2]++; + } + if(data.detect_flag[4] && data.detect_flag[6]) + { + counts[3]++; + } + if(data.detect_flag[4] && data.detect_flag[7]) + { + counts[4]++; + } + if(data.detect_flag[6] && data.detect_flag[7]) + { + counts[5]++; + } + if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6]) + { + counts[6]++; + } + if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[7]) + { + counts[7]++; + } + if(data.detect_flag[2] && data.detect_flag[6] && data.detect_flag[7]) + { + counts[8]++; + } + if(data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7]) + { + counts[9]++; + } + if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7]) + { + counts[10]++; + } + } } void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) { @@ -257,18 +418,23 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s Mask::MaskFileData data; std::vector counts; + std::vector coinc_counts; switch(header.rxn_type) { case 0: counts.resize(3, 0); + coinc_counts.resize(1, 0); break; case 1: counts.resize(4, 0); + coinc_counts.resize(1, 0); break; case 2: counts.resize(6, 0); + coinc_counts.resize(4, 0); break; case 3: counts.resize(8, 0); + coinc_counts.resize(11, 0); break; default: { @@ -284,6 +450,8 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s int count = 0; int npercent = 0; + Mask::Nucleus nucleus; + int index=0; while(true) { if(++count == percent5) {//Show progress every 5% npercent++; @@ -296,15 +464,25 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s break; for(unsigned int i=0; i= threshold && (IsRing1(data.theta[i], data.phi[i]) || IsRing2(data.theta[i], data.phi[i]) || IsQQQ(data.theta[i], data.phi[i]))) { + nucleus.SetIsotope(data.Z[i], data.A[i]); + nucleus.SetVectorSpherical(data.theta[i], data.phi[i], data.p[i], data.E[i]); + auto result = IsAnasen(nucleus); + if(result.detectFlag) { data.detect_flag[i] = true; + data.KE[i] = result.energy_deposited; + data.theta[i] = result.direction.GetTheta(); + data.phi[i] = result.direction.GetPhi(); counts[i]++; } else if(data.detect_flag[i] == true) { data.detect_flag[i] = false; } } + CountCoincidences(data, coinc_counts, header.rxn_type); + output.WriteData(data); + + index++; } input.Close(); @@ -316,6 +494,54 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s stats< #include @@ -23,11 +24,15 @@ int main(int argc, char** argv) { //sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt"); */ - + try { AnasenEfficiency anasen; anasen.CalculateEfficiency(inputname, outputname, statsname); //std::cout<<"Running consistency check(1=success): "< QQQDetector::GetTrajectoryRingWedge(double theta, double phi) } Mask::Vec3 QQQDetector::GetHitCoordinates(int ringch, int wedgech) { - if(!CheckChannel(ringch) || !CheckChannel(wedgech)) { + if(!CheckChannel(ringch) || !CheckChannel(wedgech)) return Mask::Vec3(); - } - double r_center = m_Rinner + (0.5+ringch)*m_deltaR; - double phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhi_per_wedge; + double r_center, phi_center; + if(rndmFlag) + { + r_center = m_Rinner + (m_uniform_fraction(Mask::RandomGenerator::GetInstance().GetGenerator())+ringch)*m_deltaR; + phi_center = -m_deltaPhi/2.0 + (m_uniform_fraction(Mask::RandomGenerator::GetInstance().GetGenerator())+wedgech)*m_deltaPhi_per_wedge; + } + else + { + r_center = m_Rinner + (0.5+ringch)*m_deltaR; + phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhi_per_wedge; + } double x = r_center*std::cos(phi_center); double y = r_center*std::sin(phi_center); double z = 0; diff --git a/src/Detectors/SabreEfficiency.cpp b/src/Detectors/SabreEfficiency.cpp index 6c9ff3d..411988c 100644 --- a/src/Detectors/SabreEfficiency.cpp +++ b/src/Detectors/SabreEfficiency.cpp @@ -197,9 +197,9 @@ std::pair SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) { if(dmap.IsDead(i, chan.first, 0) || dmap.IsDead(i, chan.second, 1)) break; //dead channel check coords = detectors[i].GetTrajectoryCoordinates(nucleus.GetTheta(), nucleus.GetPhi()); thetaIncident = std::acos(coords.Dot(detectors[i].GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), M_PI - thetaIncident); + eloss = deadlayer.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), M_PI - thetaIncident); if((nucleus.GetKE() - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check - e_deposited = sabre_eloss.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE() - eloss, M_PI - thetaIncident); + e_deposited = sabre_eloss.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE() - eloss, M_PI - thetaIncident); return std::make_pair(true, e_deposited); } } diff --git a/src/Detectors/StripDetector.cpp b/src/Detectors/StripDetector.cpp index b08a346..1644d89 100644 --- a/src/Detectors/StripDetector.cpp +++ b/src/Detectors/StripDetector.cpp @@ -1,19 +1,23 @@ #include "StripDetector.h" +#include /* Corner layout for each strip in the un-rotated frame 0--------------------------1 - | | y - | | ^ - | | | - | | | - | | | + | | + | | + | | + | | + | | x 2--------------------------3 z<--------X - x + | + | + | + y */ StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) : - m_random(nullptr), m_uniform_fraction(0.0, 1.0) + m_norm_unrot(1.0,0.0,0.0), m_uniform_fraction(0.0, 1.0), rndmFlag(false) { num_strips = ns; @@ -53,19 +57,19 @@ void StripDetector::CalculateCorners() { y_min = -total_width/2.0 + front_strip_width*s; z_max = center_z + length/2.0; z_min = center_z - length/2.0; - front_strip_coords[s][0] = Mask::Vec3(center_rho, y_max, z_max); - front_strip_coords[s][1] = Mask::Vec3(center_rho, y_max, z_min); - front_strip_coords[s][2] = Mask::Vec3(center_rho, y_min, z_max); - front_strip_coords[s][3] = Mask::Vec3(center_rho, y_min, z_min); + front_strip_coords[s][2] = Mask::Vec3(center_rho, y_max, z_max); + front_strip_coords[s][3] = Mask::Vec3(center_rho, y_max, z_min); + front_strip_coords[s][0] = Mask::Vec3(center_rho, y_min, z_max); + front_strip_coords[s][1] = Mask::Vec3(center_rho, y_min, z_min); z_max = (center_z - length/2.0) + (s+1)*back_strip_length; z_min = (center_z - length/2.0) + (s)*back_strip_length; y_max = total_width/2.0; y_min = -total_width/2.0; - back_strip_coords[s][0] = Mask::Vec3(center_rho, y_max, z_max); - back_strip_coords[s][1] = Mask::Vec3(center_rho, y_max, z_min); - back_strip_coords[s][2] = Mask::Vec3(center_rho, y_min, z_max); - back_strip_coords[s][3] = Mask::Vec3(center_rho, y_min, z_min); + back_strip_coords[s][2] = Mask::Vec3(center_rho, y_max, z_max); + back_strip_coords[s][3] = Mask::Vec3(center_rho, y_max, z_min); + back_strip_coords[s][0] = Mask::Vec3(center_rho, y_min, z_max); + back_strip_coords[s][1] = Mask::Vec3(center_rho, y_min, z_min); } for(int s=0; s StripDetector::GetChannelRatio(double theta, double phi) { double ratio=0; for (int s=0; s= front_strip_coords[s][0].GetX() && xhit <= front_strip_coords[s][0].GetX() && //Check min and max x (constant in flat) - yhit >= front_strip_coords[s][2].GetY() && yhit <= front_strip_coords[s][0].GetY() && //Check min and max y + yhit >= front_strip_coords[s][1].GetY() && yhit <= front_strip_coords[s][2].GetY() && //Check min and max y zhit >= front_strip_coords[s][1].GetZ() && zhit <= front_strip_coords[s][0].GetZ()) //Check min and max z { front_ch_hit = s; diff --git a/src/LayeredTarget.cpp b/src/LayeredTarget.cpp index 268f2c3..cd74e6f 100644 --- a/src/LayeredTarget.cpp +++ b/src/LayeredTarget.cpp @@ -16,7 +16,7 @@ Written by G.W. McCann Aug. 2020 namespace Mask { LayeredTarget::LayeredTarget() : - name("") + name(""), m_fractional_depth_dist(0.0, 1.0) { } @@ -31,7 +31,7 @@ namespace Mask { /* Here projectile refers to the incoming reactant particle (i.e. the beam) - Calculates energy loss assuming that the reaction occurs in the middle of the target + Calculates energy loss assuming that the reaction occurs in the middle of the target layer Note that the layer order can matter! */ double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle) { @@ -43,12 +43,14 @@ namespace Mask { double eloss = 0.0; double newEnergy = startEnergy; + double frac; for(int i=0; i<=rxnLayer; i++) { if(i == rxnLayer) { - eloss += layers[i].getEnergyLossHalf(zp, ap, newEnergy, angle); + frac = m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator()); + eloss += layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, frac); newEnergy = startEnergy - eloss; } else { - eloss += layers[i].getEnergyLossTotal(zp, ap, newEnergy, angle); + eloss += layers[i].GetEnergyLossTotal(zp, ap, newEnergy, angle); newEnergy = startEnergy-eloss; } } @@ -69,14 +71,27 @@ namespace Mask { } double eloss = 0.0; - double newEnergy = startEnergy; - for(unsigned int i=rxnLayer; i=0; i--) { + if(i == ((unsigned int)rxnLayer)) { + eloss += layers[i].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator())); + newEnergy = startEnergy - eloss; + } else { + eloss += layers[i].GetEnergyLossTotal(ze, ae, newEnergy, angle); + newEnergy = startEnergy - eloss; + } } } @@ -92,14 +107,27 @@ namespace Mask { } double eloss = 0.0; - double newEnergy = startEnergy; - for(int i=(layers.size()-1); i>=rxnLayer; i--) { - if(i == rxnLayer) { - eloss += layers[i].getReverseEnergyLossHalf(ze, ae, newEnergy, angle); - newEnergy = startEnergy + eloss; - } else { - eloss += layers[i].getReverseEnergyLossTotal(ze, ae, newEnergy, angle); - newEnergy = startEnergy + eloss; + if(angle < M_PI/2.0) { + double newEnergy = startEnergy; + for(int i=(layers.size()-1); i>=rxnLayer; i--) { + if(i == rxnLayer) { + eloss += layers[i].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator())); + newEnergy = startEnergy + eloss; + } else { + eloss += layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle); + newEnergy = startEnergy + eloss; + } + } + } else { + double newEnergy = startEnergy; + for(int i=0; i <= rxnLayer; i++) { + if(i == rxnLayer) { + eloss += layers[i].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator())); + newEnergy = startEnergy + eloss; + } else { + eloss += layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle); + newEnergy = startEnergy + eloss; + } } } diff --git a/src/Target.cpp b/src/Target.cpp index 83e3599..b3a55f7 100644 --- a/src/Target.cpp +++ b/src/Target.cpp @@ -39,7 +39,7 @@ namespace Mask { } /*Calculates energy loss for travelling all the way through the target*/ - double Target::getEnergyLossTotal(int zp, int ap, double startEnergy, double theta) { + double Target::GetEnergyLossTotal(int zp, int ap, double startEnergy, double theta) { if(theta == M_PI/2.) return startEnergy; else if (theta > M_PI/2.) @@ -47,19 +47,20 @@ namespace Mask { return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta))); } - - /*Calculates energy loss for travelling halfway through the target*/ - double Target::getEnergyLossHalf(int zp, int ap, double startEnergy, double theta) { - if(theta == M_PI/2.) - return startEnergy; - else if (theta > M_PI/2.) - theta = M_PI - theta; - return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/(2.0*fabs(cos(theta)))); + /*Calculates the energy loss for traveling some fraction through the target*/ + double Target::GetEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double theta, double percent_depth) + { + if(theta == M_PI/2.) + return finalEnergy; + else if (theta > M_PI/2.) + theta = M_PI-theta; + + return eloss.GetEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta)))); } /*Calculates reverse energy loss for travelling all the way through the target*/ - double Target::getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta) { + double Target::GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta) { if(theta == M_PI/2.) return finalEnergy; else if (theta > M_PI/2.) @@ -67,15 +68,16 @@ namespace Mask { return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta))); } - - /*Calculates reverse energy loss for travelling half way through the target*/ - double Target::getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double theta) { - if(theta == M_PI/2.) - return finalEnergy; - else if (theta > M_PI/2.) - theta = M_PI - theta; - return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta)))); + /*Calculates the reverse energy loss for traveling some fraction through the target*/ + double Target::GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double theta, double percent_depth) + { + if(theta == M_PI/2.) + return finalEnergy; + else if (theta > M_PI/2.) + theta = M_PI-theta; + + return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta)))); } } diff --git a/src/TwoStepSystem.cpp b/src/TwoStepSystem.cpp index 55530e5..4c3fd27 100644 --- a/src/TwoStepSystem.cpp +++ b/src/TwoStepSystem.cpp @@ -73,7 +73,7 @@ namespace Mask { double decay1costheta = decay1dist.GetRandomCosTheta(); double decay1Theta = std::acos(decay1costheta); double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator()); - double residEx = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator()); + double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator()); step1.SetBeamKE(bke); step1.SetPolarRxnAngle(rxnTheta); @@ -97,4 +97,4 @@ namespace Mask { } -} \ No newline at end of file +}