From 0b6319607760db927cb6f794b28badb3eb08fcc2 Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Thu, 9 Jun 2022 11:29:02 -0400 Subject: [PATCH] Added catima for better energy loss calculations --- .gitmodules | 3 +++ CMakeLists.txt | 1 + include/MassLookup.h | 1 + include/SabreEfficiency.h | 2 ++ include/Target.h | 4 ++++ include/Vec4.h | 4 ++-- src/Detectors/AnasenEfficiency.cpp | 1 + src/Detectors/CMakeLists.txt | 1 + src/Mask/CMakeLists.txt | 1 + src/Mask/Target.cpp | 34 ++++++++++++++++++++++++------ src/Plotters/ROOT/RootPlotter.cpp | 15 +++++++++++++ src/vendor/catima | 1 + 12 files changed, 59 insertions(+), 9 deletions(-) create mode 100644 .gitmodules create mode 160000 src/vendor/catima diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..0f578de --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "src/vendor/catima"] + path = src/vendor/catima + url = https://github.com/gwm17/catima.git diff --git a/CMakeLists.txt b/CMakeLists.txt index 060a36e..a507986 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,6 +8,7 @@ set(MASK_INCLUDE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/include) set(CMAKE_CXX_STANDARD 17) +add_subdirectory(src/vendor/catima) add_subdirectory(src/Mask) add_subdirectory(src/MaskApp) add_subdirectory(src/Detectors) diff --git a/include/MassLookup.h b/include/MassLookup.h index 3ca695d..5bba409 100644 --- a/include/MassLookup.h +++ b/include/MassLookup.h @@ -22,6 +22,7 @@ namespace Mask { public: ~MassLookup(); double FindMass(int Z, int A); + double FindMassU(int Z, int A) { return FindMass(Z, A)/u_to_mev; } std::string FindSymbol(int Z, int A); static MassLookup& GetInstance() { diff --git a/include/SabreEfficiency.h b/include/SabreEfficiency.h index ea866c6..03c18cd 100644 --- a/include/SabreEfficiency.h +++ b/include/SabreEfficiency.h @@ -25,6 +25,7 @@ private: Mask::Target deadlayer; Mask::Target sabre_eloss; + Mask::Target degrader; SabreDeadChannelMap dmap; //Sabre constants @@ -41,6 +42,7 @@ private: const double DEG2RAD = M_PI/180.0; static constexpr double DEADLAYER_THIN = 50 * 1e-7 * 2.3296 * 1e6; // ug/cm^2 (50 nm thick * density) static constexpr double SABRE_THICKNESS = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density) + static constexpr double DEGRADER_THICKNESS = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick) const double ENERGY_THRESHOLD = 0.2; //in MeV diff --git a/include/Target.h b/include/Target.h index bcaafd9..83839eb 100644 --- a/include/Target.h +++ b/include/Target.h @@ -17,6 +17,8 @@ Written by G.W. McCann Aug. 2020 #include #include #include "EnergyLoss.h" +#include "catima/gwm_integrators.h" +#include "MassLookup.h" namespace Mask { @@ -39,7 +41,9 @@ namespace Mask { private: EnergyLoss eloss; + catima::Material target_material; double thickness; + double thickness_gcm2; std::vector Z, A, Stoich; }; diff --git a/include/Vec4.h b/include/Vec4.h index 15ca5b3..51eb879 100644 --- a/include/Vec4.h +++ b/include/Vec4.h @@ -15,7 +15,7 @@ namespace Mask { public: Vec4(); Vec4(double px, double py, double pz, double E); - virtual ~Vec4(); + virtual ~Vec4(); void SetVectorCartesian(double px, double py, double pz, double E); void SetVectorSpherical(double theta, double phi, double p, double E); @@ -66,4 +66,4 @@ namespace Mask { } -#endif \ No newline at end of file +#endif diff --git a/src/Detectors/AnasenEfficiency.cpp b/src/Detectors/AnasenEfficiency.cpp index 9a191a4..be71cce 100644 --- a/src/Detectors/AnasenEfficiency.cpp +++ b/src/Detectors/AnasenEfficiency.cpp @@ -1,6 +1,7 @@ #include "AnasenEfficiency.h" #include #include +#include AnasenEfficiency::AnasenEfficiency() : DetectorEfficiency(), det_silicon(si_thickness) diff --git a/src/Detectors/CMakeLists.txt b/src/Detectors/CMakeLists.txt index 46d2b16..1852b72 100644 --- a/src/Detectors/CMakeLists.txt +++ b/src/Detectors/CMakeLists.txt @@ -14,6 +14,7 @@ target_sources(DetectEff PUBLIC target_link_libraries(DetectEff Mask + catima ) set_target_properties(DetectEff PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR}) \ No newline at end of file diff --git a/src/Mask/CMakeLists.txt b/src/Mask/CMakeLists.txt index f808410..cb12a41 100644 --- a/src/Mask/CMakeLists.txt +++ b/src/Mask/CMakeLists.txt @@ -26,4 +26,5 @@ target_sources(Mask PRIVATE Vec4.cpp ) +target_link_libraries(Mask catima) set_target_properties(Mask PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR}) \ No newline at end of file diff --git a/src/Mask/Target.cpp b/src/Mask/Target.cpp index b3a55f7..0a2a908 100644 --- a/src/Mask/Target.cpp +++ b/src/Mask/Target.cpp @@ -11,12 +11,14 @@ Written by G.W. McCann Aug. 2020 */ #include "Target.h" +#include "catima/nucdata.h" namespace Mask { /*Targets must be of known thickness*/ Target::Target(double thick) { thickness = thick; + thickness_gcm2 = thickness*1.0e-6; } Target::~Target() {} @@ -26,8 +28,12 @@ namespace Mask { Z = z; A = a; Stoich = stoich; - + MassLookup& masses = MassLookup::GetInstance(); eloss.SetTargetComponents(Z, A, Stoich); + for(size_t i=0; i M_PI/2.) theta = M_PI - theta; - - return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta))); + catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0); + proj.T = startEnergy/proj.A; + target_material.thickness(thickness_gcm2/fabs(cos(theta))); + return catima::integrate_energyloss(proj, target_material); + //return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta))); } /*Calculates the energy loss for traveling some fraction through the target*/ @@ -56,7 +65,11 @@ namespace Mask { else if (theta > M_PI/2.) theta = M_PI-theta; - return eloss.GetEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta)))); + catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0); + proj.T = finalEnergy/proj.A; + target_material.thickness(thickness_gcm2*percent_depth/fabs(cos(theta))); + return catima::integrate_energyloss(proj, target_material); + //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*/ @@ -66,7 +79,11 @@ namespace Mask { else if (theta > M_PI/2.) theta = M_PI - theta; - return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta))); + catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0); + proj.T = finalEnergy/proj.A; + target_material.thickness(thickness_gcm2/fabs(cos(theta))); + return catima::reverse_integrate_energyloss(proj, target_material); + //return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta))); } /*Calculates the reverse energy loss for traveling some fraction through the target*/ @@ -76,8 +93,11 @@ namespace Mask { 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)))); + catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0); + proj.T = finalEnergy/proj.A; + target_material.thickness(thickness_gcm2*percent_depth/fabs(cos(theta))); + return catima::reverse_integrate_energyloss(proj, target_material); + //return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta)))); } } diff --git a/src/Plotters/ROOT/RootPlotter.cpp b/src/Plotters/ROOT/RootPlotter.cpp index c1176df..e6cc0be 100644 --- a/src/Plotters/ROOT/RootPlotter.cpp +++ b/src/Plotters/ROOT/RootPlotter.cpp @@ -1,5 +1,6 @@ #include "RootPlotter.h" #include +#include #include @@ -57,12 +58,19 @@ void RootPlotter::FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType if(type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) { + TVector3 p1, p2; + p1.SetMagThetaPhi(1.0, data.theta[3], data.phi[3]); + p2.SetMagThetaPhi(1.0, data.theta[4], data.phi[4]); + double theta_resid_break1 = std::acos(p1.Dot(p2)); std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor"; std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)"; MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4); std::string theta_resid_theta_break1_name = "theta_resid_theta_break1_cor"; std::string theta_resid_theta_break1_title = theta_resid_theta_break1_name + ";#theta_{lab} Residual (deg);#theta_{lab} Breakup1 (deg)"; MyFill(theta_resid_theta_break1_name, theta_resid_theta_break1_title, data.theta[3]*rad2deg, data.theta[4]*rad2deg, 4); + std::string ke_break1_theta_rel_name = "ke_break1_theta_rel"; + std::string ke_break1_theta_rel_title = ke_break1_theta_rel_name + ";#theta_{resid-break1};KE_{break1} (MeV)"; + MyFill(ke_break1_theta_rel_name, ke_break1_theta_rel_title, theta_resid_break1*rad2deg, data.KE[4], 4); } if(type == Mask::RxnType::ThreeStepRxn) { @@ -86,9 +94,16 @@ void RootPlotter::FillCorrelationsDetected(const Mask::MaskFileData& data, Mask: } if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4]) { + TVector3 p1, p2; + p1.SetMagThetaPhi(1.0, data.theta[3], data.phi[3]); + p2.SetMagThetaPhi(1.0, data.theta[4], data.phi[4]); + double theta_resid_break1 = std::acos(p1.Dot(p2)); std::string theta_resid_theta_break1_name = "theta_resid_theta_break1_cor_detected"; std::string theta_resid_theta_break1_title = theta_resid_theta_break1_name + ";#theta_{lab} Residual (deg);#theta_{lab} Breakup1 (deg)"; MyFill(theta_resid_theta_break1_name, theta_resid_theta_break1_title, data.theta[3]*rad2deg, data.theta[4]*rad2deg, 4); + std::string ke_break1_theta_rel_name = "ke_break1_theta_rel_detected"; + std::string ke_break1_theta_rel_title = ke_break1_theta_rel_name + ";#theta_{resid-break1};KE_{break1} (MeV)"; + MyFill(ke_break1_theta_rel_name, ke_break1_theta_rel_title, theta_resid_break1*rad2deg, data.KE[4], 4); } if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4] && data.detect_flag[5]) diff --git a/src/vendor/catima b/src/vendor/catima new file mode 160000 index 0000000..629690a --- /dev/null +++ b/src/vendor/catima @@ -0,0 +1 @@ +Subproject commit 629690a6f90ad9d343e02f72d4fad860fce1ed6d