diff --git a/.gitignore b/.gitignore index 47a980a..965ada5 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ *.o *.so +*.dylib *.sublime-project *.sublime-workspace *.code-workspace 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 3d5db2c..fdec1b1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,4 +7,5 @@ set(SABRERECON_LIBRARY_DIR ${CMAKE_CURRENT_SOURCE_DIR}/lib) set(CMAKE_CXX_STANDARD 17) +add_subdirectory(src/vendor/catima) add_subdirectory(src) diff --git a/input.txt b/input.txt index d71d634..4d09aa1 100644 --- a/input.txt +++ b/input.txt @@ -1,6 +1,6 @@ begin_data - input /Volumes/Wyndle/10B3He_May2022/calibrated/run_53_149_cal.root - output /Volumes/Wyndle/10B3He_May2022/histograms/calHistograms_53_149.root + input /Volumes/Wyndle/10B3He_May2022/calibrated/run_53_143_cal.root + output /Volumes/Wyndle/10B3He_May2022/histograms/calHistograms_53_143_newEloss.root beamKE(MeV) 24.0 end_data begin_reconstructor @@ -8,18 +8,18 @@ begin_reconstructor B 8.759 theta 15.0 begin_fpcal - 78.5352 - 0.038 + 78.5946 + 0.0382 end_fpcal end_focalplane begin_target thickness 74.0 begin_elements - 5 1 + 5 10 1 end_elements end_target end_reconstructor begin_cuts ede_alphas /Volumes/Wyndle/10B3He_May2022/cuts/edeCut_alphas.root scintE cathodeE x1x2 /Volumes/Wyndle/10B3He_May2022/cuts/x1x2Cut.root x1 x2 -end_cuts \ No newline at end of file +end_cuts diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c8ee0dc..5c1828a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -30,6 +30,7 @@ target_sources(SabreRecon PUBLIC ) target_link_libraries(SabreRecon CalDict + catima ${ROOT_LIBRARIES} ) set_target_properties(SabreRecon PROPERTIES diff --git a/src/EnergyLoss/Target.cpp b/src/EnergyLoss/Target.cpp index 9132ac8..c0f5900 100644 --- a/src/EnergyLoss/Target.cpp +++ b/src/EnergyLoss/Target.cpp @@ -23,10 +23,10 @@ namespace SabreRecon { } /*Targets must be of known thickness*/ - Target::Target(const std::vector& z, const std::vector& stoich, double thick) : + Target::Target(const std::vector& a, const std::vector& z, const std::vector& stoich, double thick) : m_isValid(false) { - SetParameters(z, stoich, thick); + SetParameters(a, z, stoich, thick); } Target::~Target() @@ -34,7 +34,7 @@ namespace SabreRecon { } /*Set target elements of given Z, A, S*/ - void Target::SetParameters(const std::vector& z, const std::vector& stoich, double thick) + void Target::SetParameters(const std::vector& a, const std::vector& z, const std::vector& stoich, double thick) { m_params.ZT = z; double denom = 0; @@ -43,6 +43,13 @@ namespace SabreRecon { for(auto& s : stoich) m_params.composition.push_back(s/denom); m_totalThickness = thick; + m_totalThickness_gcm2 = m_totalThickness*1.0e-6; + + auto& masses = MassLookup::GetInstance(); + for(size_t i=0; i #include #include "EnergyLoss.h" +#include "catima/gwm_integrators.h" namespace SabreRecon { @@ -24,10 +25,10 @@ namespace SabreRecon { public: Target(); - Target(const std::vector& z, const std::vector& stoich, double thick); + Target(const std::vector& a, const std::vector& z, const std::vector& stoich, double thick); ~Target(); - void SetParameters(const std::vector& z, const std::vector& stoich, double thick); + void SetParameters(const std::vector& a, const std::vector& z, const std::vector& stoich, double thick); 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); @@ -39,7 +40,10 @@ namespace SabreRecon { private: EnergyLoss::Parameters m_params; + catima::Material m_material; + catima::Projectile m_projectile; double m_totalThickness; + double m_totalThickness_gcm2; bool m_isValid; }; diff --git a/src/Histogrammer.cpp b/src/Histogrammer.cpp index b877337..5066f92 100644 --- a/src/Histogrammer.cpp +++ b/src/Histogrammer.cpp @@ -40,6 +40,7 @@ namespace SabreRecon { double thickness; std::vector targ_z; + std::vector targ_a; std::vector targ_s; std::vector cuts; @@ -105,8 +106,10 @@ namespace SabreRecon { { targ_z.push_back(std::stoi(junk)); input>>junk; + targ_a.push_back(std::stoi(junk)); + input>>junk; targ_s.push_back(std::stoi(junk)); - std::cout<<"e"<Get("CUTG"); protonGate->SetName("protonPunchGate"); @@ -241,11 +245,14 @@ namespace SabreRecon { if(!m_eventPtr->sabre.empty() && m_eventPtr->sabre[0].ringE > s_weakSabreThreshold) { auto& biggestSabre = m_eventPtr->sabre[0]; + recon9B = m_recon.RunFPResidExcitation(m_eventPtr->xavg, m_beamKE, {{5,10},{2,3},{3,4}}); recon5Li = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{2,4}}); recon8Be = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{1,1}}); recon7Be = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{1,2}}); recon14N = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{8,16},{2,3},{2,4},{1,1}}); sabreCoords = m_recon.GetSabreCoordinates(biggestSabre); + b9Coords.SetMagThetaPhi(1.0, recon9B.residThetaLab, recon9B.residPhiLab); + relAngle = std::acos(b9Coords.Dot(sabreCoords)/(sabreCoords.Mag()*b9Coords.Mag())); FillHistogram1D({"xavg_gated_sabre","xavg_gated_sabre;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg); FillHistogram2D({"scintE_cathodeE","scintE_cathodeE;scintE;cathodeE",512,0,4096,512,0,4096}, @@ -262,6 +269,10 @@ namespace SabreRecon { sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE); FillHistogram2D({"xavg_sabreE","xavg_sabreE;xavg; E(MeV)",600,-300.0,300.0,400,0,20.0}, m_eventPtr->xavg, biggestSabre.ringE); + FillHistogram2D({"9Btheta_sabreTheta","9Btheta_sabreTheta;#theta_{9B};#theta_{SABRE}",180,0.0,180.0, + 180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg); + FillHistogram2D({"sabreE_relAngle","sabreE_relAngle;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg, + biggestSabre.ringE); if(m_eventPtr->xavg > -186.0 && m_eventPtr->xavg < -178.0 && (biggestSabre.detID == 2 || biggestSabre.detID == 3)) { @@ -297,9 +308,11 @@ namespace SabreRecon { { FillHistogram1D({"xavg_gated14Ngs", "xavg_gated14Ngs;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg); } - else + + if(!(recon14N.excitation > -0.1 && recon14N.excitation < 0.1) && !(recon7Be.excitation > -0.1 && recon7Be.excitation < 0.15) + && !(recon8Be.excitation > -0.1 && recon8Be.excitation < 0.1) && !(recon5Li.excitation > -1.5 && recon5Li.excitation < 1.5)) { - FillHistogram1D({"xavg_notGated14Ngs", "xavg_notGated14Ngs;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg); + FillHistogram1D({"xavg_notGatedAllChannels", "xavg_notGatedAllChannels;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg); } //Degrader analysis... SABRE detectors 0, 1, 4 are covered with tantalum @@ -309,6 +322,10 @@ namespace SabreRecon { 180,0.0,180.0,400,0.0,20.0}, sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE); FillHistogram2D({"xavg_sabreE_degraded", "xavg_sabreE_degraded;xavg;E(MeV)", 600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE); + FillHistogram2D({"9Btheta_sabreTheta_degraded","9Btheta_sabreTheta_degraded;#theta_{9B};#theta_{SABRE}",180,0.0,180.0, + 180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg); + FillHistogram2D({"sabreE_relAngle_degraded","sabreE_relAngle_degraded;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg, + biggestSabre.ringE); if(biggestSabre.local_wedge != 0 && biggestSabre.local_wedge != 7 && biggestSabre.local_ring != 15 && biggestSabre.local_ring != 0) //Edges might not be degraded right { FillHistogram2D({"sabreE_sabreTheta_degraded_rejectEdge", "sabreE_sabreTheta_degraded;#theta (deg);E(MeV)", @@ -316,7 +333,12 @@ namespace SabreRecon { FillHistogram2D({"xavg_sabreE_degraded_rejectEdge", "xavg_sabreE_degraded;xavg;E(MeV)", 600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE); FillHistogram1D({"xavg_degraded_rejectEdge","xavg_degraded_rejectEdge;xavg",600,-300.0,300.0}, m_eventPtr->xavg); - if(protonGate->IsInside(sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE)) + FillHistogram2D({"9Btheta_sabreTheta_degraded_rejectEdge","9Btheta_sabreTheta_degraded_rejectEdge;#theta_{9B};#theta_{SABRE}",180,0.0,180.0, + 180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg); + FillHistogram2D({"sabreE_relAngle_degraded_rejectEdge","sabreE_relAngle_degraded_rejectEdge;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg, + biggestSabre.ringE); + + if(protonGate->IsInside(relAngle*s_rad2deg, biggestSabre.ringE)) { FillHistogram2D({"xavg_sabreE_degraded_rejectEdge_pGate", "xavg_sabreE_degraded;xavg;E(MeV)", 600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE); diff --git a/src/MassLookup.h b/src/MassLookup.h index 53b11c0..2992a56 100644 --- a/src/MassLookup.h +++ b/src/MassLookup.h @@ -24,6 +24,7 @@ namespace SabreRecon { MassLookup(); ~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); inline static MassLookup& GetInstance() { return *s_instance; } diff --git a/src/Reconstructor.cpp b/src/Reconstructor.cpp index d52e507..8a4dc49 100644 --- a/src/Reconstructor.cpp +++ b/src/Reconstructor.cpp @@ -187,12 +187,15 @@ namespace SabreRecon { auto resid_vec = targ_vec + proj_vec - eject_vec; result.excitation = resid_vec.M() - massResid; + result.residThetaLab = resid_vec.Theta(); + result.residPhiLab = resid_vec.Phi(); auto parent_vec = targ_vec + proj_vec; auto boost = parent_vec.BoostVector(); eject_vec.Boost(-1.0*boost); - result.theta_cm = eject_vec.Theta(); - result.phi_cm = eject_vec.Phi(); + + result.ejectThetaCM = eject_vec.Theta(); + result.ejectPhiCM = eject_vec.Phi(); return result; } @@ -235,8 +238,8 @@ namespace SabreRecon { auto parent_vec = targ_vec + proj_vec; auto boost = parent_vec.BoostVector(); eject_vec.Boost(-1.0*boost); - result.theta_cm = eject_vec.Theta(); - result.phi_cm = eject_vec.Phi(); + result.ejectThetaCM = eject_vec.Theta(); + result.ejectPhiCM = eject_vec.Phi(); return result; } @@ -279,8 +282,8 @@ namespace SabreRecon { result.excitation = decayFrag_vec.M() - massDecayFrag; auto boost = resid_vec.BoostVector(); decayBreak_vec.Boost(-1.0*boost); - result.theta_cm = decayBreak_vec.Theta(); - result.phi_cm = decayBreak_vec.Phi(); + result.ejectThetaCM = decayBreak_vec.Theta(); + result.ejectPhiCM = decayBreak_vec.Phi(); return result; } @@ -323,8 +326,8 @@ namespace SabreRecon { result.excitation = decayFrag_vec.M() - massDecayFrag; auto boost = resid_vec.BoostVector(); decayBreak_vec.Boost(-1.0*boost); - result.theta_cm = decayBreak_vec.Theta(); - result.phi_cm = decayBreak_vec.Phi(); + result.ejectThetaCM = decayBreak_vec.Theta(); + result.ejectPhiCM = decayBreak_vec.Phi(); return result; } diff --git a/src/Reconstructor.h b/src/Reconstructor.h index 395f2ab..fc9781d 100644 --- a/src/Reconstructor.h +++ b/src/Reconstructor.h @@ -14,8 +14,12 @@ namespace SabreRecon { struct ReconResult { double excitation; - double theta_cm; - double phi_cm; + double ejectThetaCM; + double ejectPhiCM; + double residThetaLab; + double residPhiLab; + double residThetaCM; + double residPhiCM; }; struct NucID 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