From 15020355a686604ad9401db4127d84abedffb0fb Mon Sep 17 00:00:00 2001 From: gwm17 Date: Tue, 30 Aug 2022 10:34:16 -0400 Subject: [PATCH 1/6] Fix naming convention on executables: MaskApp->Kinematics, DetectEff->Detectors. Update README --- README.md | 20 ++++++++++---------- src/Detectors/CMakeLists.txt | 12 +++++------- src/Kinematics/CMakeLists.txt | 12 ++++++++++++ src/{MaskApp => Kinematics}/main.cpp | 0 src/MaskApp/CMakeLists.txt | 12 ------------ src/Plotters/CMakeLists.txt | 1 - 6 files changed, 27 insertions(+), 30 deletions(-) create mode 100644 src/Kinematics/CMakeLists.txt rename src/{MaskApp => Kinematics}/main.cpp (100%) delete mode 100644 src/MaskApp/CMakeLists.txt diff --git a/README.md b/README.md index f28a651..cb1950b 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,8 @@ -# MASK: Monte cArlo Simulation of Kinematics -MASK is a Monte Carlo simulation of reaction kinematics for use detector systems at Florida State University. -MASK is capable of simulating multi-step kinematic reaction-decay sequences, storing data in ROOT Trees, after which the kinematic data can be fed to a detector geometry for efficiency testing. Currently geometries for ANASEN and SABRE are included in the code. +# Mask: Monte cArlo Simulation of Kinematics +Mask is a Monte Carlo simulation of reaction kinematics for use detector systems at Florida State University. +Mask is capable of simulating multi-step kinematic reaction-decay sequences, storing data in ROOT Trees, after which the kinematic data can be fed to a detector geometry for efficiency testing. Currently geometries for ANASEN and SABRE are included in the code. -## Building MASK +## Building Mask Dowload the repository from github. CMake is use to build the project; in most environments you can build Mask using the following methods: - `mkdir build` - `cd build` @@ -14,7 +14,7 @@ Executables will be installed to the repostiory's `bin` directory. Libraries wil By default Mask builds for release. To build for debug replace `cmake ..` with `cmake -DCMAKE_BUILD_TYPE=Debug ..`. Mask uses CMake to find the installed ROOT libraries and headers. ## Using the kinematics simulation -By default MASK is capable of simulating reactions of up to three steps. Here is a brief outline of each type: +By default Mask is capable of simulating reactions of up to three steps. Here is a brief outline of each type: 0. A pure decay involves a "target" decaying into an ejectile and residual. It is assumed isotropic by default; any other case will require the modification of the code. 1. A single step reaction involves a target being hit by a projectile and emitting an ejectile and a residual. It can incorporate all of the input file sampling parameters. @@ -23,9 +23,9 @@ By default MASK is capable of simulating reactions of up to three steps. Here is For decays, a specific angular distribution can be given as input as a text file with values of coefficiencts of a Legendre polynomial series. Examples can be found in the `./etc` directory, including an isotropic case. It is assumed that the decays in the center-of-mass frame are isotropic in phi (i.e. m=0). Decay1 corresponds to the first decay, if there are multiple steps, Decay2 to the second. If there are no decays, these parameters are not used (or if only one decay, Decay2_AngularMomentum is not used). The input file requires that the user include target information, which will be used to calculate energy loss for all of the reactants and reaction products. The energy loss through materials is calculated using the `catima` library (found in `src/vendor`), which is a C/C++ interface to the `atima` library (the same energy loss methods used by LISE). The target can contain layers, and each layer can be composed of a compound of elements with a given stoichiometry. If the user wishes to not include energy loss in the kinematics, simply give all target layers a thickness of 0. Note that more layers and more thickness = more time spent calculating energy loss. These energy loss methods are only applicable for solid targets, and should not be applied to gas or liquid targets. Energy loss calculations have a stated uncertainty of approximately five percent. -To run MASK simply do the following from the MASK repository: +To run Mask simply do the following from the Mask repository: -`./bin/MaskApp input.txt` +`./bin/Kinematics input.txt` Input.txt can be replaced by any text file with the correct format. @@ -34,14 +34,14 @@ Detector geometry is encoded using ROOT math libraries in the `src/Detectors` fo To choose which detector scheme is run, modify the main function in `src/Detectors/main.cpp`. The included geometries also have options to do an internal geometry consistency check and print out coordinates for drawing the detector arrays, which can be useful for testing. -To run DetEff use the format +To run Detectors use the format -`./bin/DetEff ` +`./bin/Detectors ` where the detection datafile contains all of the kinematics data as well as information about which particles are detected and the statsfile is a text file containing efficiency statistics. ## Data visualization -All data is saved as ROOT trees of std::vectors of Mask::Nucleus classes. To enable this, a ROOT dictionary is generated and linked into a shared library found in the `lib` directory of the repository. This allows the user to link to the shared library for accessing and analyzing the data generated by MASK. +All data is saved as ROOT trees of std::vectors of Mask::Nucleus classes. To enable this, a ROOT dictionary is generated and linked into a shared library found in the `lib` directory of the repository. This allows the user to link to the shared library for accessing and analyzing the data generated by Mask. Mask also provides a default visualization tool called RootPlot. RootPlot is run as diff --git a/src/Detectors/CMakeLists.txt b/src/Detectors/CMakeLists.txt index 0165619..c3c5232 100644 --- a/src/Detectors/CMakeLists.txt +++ b/src/Detectors/CMakeLists.txt @@ -1,7 +1,7 @@ -add_executable(DetectEff) -target_include_directories(DetectEff PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +add_executable(Detectors) +target_include_directories(Detectors PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) -target_sources(DetectEff PUBLIC +target_sources(Detectors PUBLIC AnasenDeadChannelMap.cpp AnasenDeadChannelMap.h AnasenEfficiency.cpp @@ -20,10 +20,8 @@ target_sources(DetectEff PUBLIC StripDetector.h ) -target_link_libraries(DetectEff - MaskDict +target_link_libraries(Detectors Mask - catima ) -set_target_properties(DetectEff PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR}) \ No newline at end of file +set_target_properties(Detectors PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR}) \ No newline at end of file diff --git a/src/Kinematics/CMakeLists.txt b/src/Kinematics/CMakeLists.txt new file mode 100644 index 0000000..4b54669 --- /dev/null +++ b/src/Kinematics/CMakeLists.txt @@ -0,0 +1,12 @@ +add_executable(Kinematics) +target_include_directories(Kinematics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/..) + +target_sources(Kinematics PUBLIC + main.cpp +) + +target_link_libraries(Kinematics + Mask +) + +set_target_properties(Kinematics PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR}) \ No newline at end of file diff --git a/src/MaskApp/main.cpp b/src/Kinematics/main.cpp similarity index 100% rename from src/MaskApp/main.cpp rename to src/Kinematics/main.cpp diff --git a/src/MaskApp/CMakeLists.txt b/src/MaskApp/CMakeLists.txt deleted file mode 100644 index 8c4e80c..0000000 --- a/src/MaskApp/CMakeLists.txt +++ /dev/null @@ -1,12 +0,0 @@ -add_executable(MaskApp) -target_include_directories(MaskApp PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/..) - -target_sources(MaskApp PUBLIC - main.cpp -) - -target_link_libraries(MaskApp - Mask -) - -set_target_properties(MaskApp PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR}) \ No newline at end of file diff --git a/src/Plotters/CMakeLists.txt b/src/Plotters/CMakeLists.txt index 976effa..dc9d2d0 100644 --- a/src/Plotters/CMakeLists.txt +++ b/src/Plotters/CMakeLists.txt @@ -12,7 +12,6 @@ target_sources(RootPlot PUBLIC ) target_link_libraries(RootPlot - MaskDict Mask ${ROOT_LIBRARIES} ) From 45a80ab76f8fa6a75f1425c28f917fac68788acd Mon Sep 17 00:00:00 2001 From: gwm17 Date: Tue, 30 Aug 2022 14:19:42 -0400 Subject: [PATCH 2/6] Enforce dictionary linking. Update cmake list --- CMakeLists.txt | 2 +- src/Detectors/main.cpp | 6 ++++++ src/Plotters/main.cpp | 6 ++++++ 3 files changed, 13 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 55d5672..ebf7578 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -20,6 +20,6 @@ find_package(ROOT REQUIRED COMPONENTS GenVector) add_subdirectory(src/vendor/catima) add_subdirectory(src/Mask) -add_subdirectory(src/MaskApp) +add_subdirectory(src/Kinematics) add_subdirectory(src/Detectors) add_subdirectory(src/Plotters) \ No newline at end of file diff --git a/src/Detectors/main.cpp b/src/Detectors/main.cpp index b03a378..ba6f8f6 100644 --- a/src/Detectors/main.cpp +++ b/src/Detectors/main.cpp @@ -13,6 +13,12 @@ int main(int argc, char** argv) return 1; } + if(!Mask::EnforceDictionaryLinked()) + { + std::cerr<<"This should be illegal!"< Date: Tue, 30 Aug 2022 15:01:11 -0400 Subject: [PATCH 3/6] Fix up detector geometry, make names clearer for dets and constexpr more of the geometry --- src/Detectors/AnasenEfficiency.cpp | 8 +- src/Detectors/AnasenEfficiency.h | 6 +- src/Detectors/CMakeLists.txt | 4 +- src/Detectors/QQQDetector.cpp | 47 ++++++------ src/Detectors/QQQDetector.h | 12 +-- .../{StripDetector.cpp => SX3Detector.cpp} | 73 ++++++++----------- .../{StripDetector.h => SX3Detector.h} | 25 ++++--- src/Detectors/SabreDetector.cpp | 53 ++++++-------- src/Detectors/SabreDetector.h | 53 ++++++-------- src/Detectors/SabreEfficiency.cpp | 2 +- src/Detectors/SabreEfficiency.h | 3 - 11 files changed, 123 insertions(+), 163 deletions(-) rename src/Detectors/{StripDetector.cpp => SX3Detector.cpp} (70%) rename src/Detectors/{StripDetector.h => SX3Detector.h} (76%) diff --git a/src/Detectors/AnasenEfficiency.cpp b/src/Detectors/AnasenEfficiency.cpp index 5a0b7d8..8c6cd19 100644 --- a/src/Detectors/AnasenEfficiency.cpp +++ b/src/Detectors/AnasenEfficiency.cpp @@ -11,16 +11,16 @@ AnasenEfficiency::AnasenEfficiency() : { for(int i=0; i #include "DetectorEfficiency.h" -#include "StripDetector.h" +#include "SX3Detector.h" #include "QQQDetector.h" #include "Target.h" #include "Nucleus.h" @@ -39,15 +39,11 @@ private: static constexpr int s_nSX3PerBarrel = 12; static constexpr int s_nQQQ = 4; static constexpr double s_sx3Length = 0.075; - static constexpr double s_sx3Width = 0.04; static constexpr double s_barrelGap = 0.0254; static constexpr double s_sx3FrameGap = 0.049; //0.049 is base gap due to frames static constexpr double s_barrel1Z = s_sx3Length/2.0 + s_sx3FrameGap + s_barrelGap/2.0; static constexpr double s_barrel2Z = (-1.0)*(s_barrelGap/2.0 + s_sx3Length/2.0); static constexpr double s_qqqZ = 0.0125 + s_sx3Length + s_sx3FrameGap + s_barrelGap/2.0; - static constexpr double s_qqqInnerR = 0.0501; - static constexpr double s_qqqOuterR = 0.0990; - static constexpr double s_qqqDeltaPhi = 1.52119; static constexpr double s_qqqZList[4] = {s_qqqZ, s_qqqZ, s_qqqZ, s_qqqZ}; static constexpr double s_qqqPhiList[4] = {5.49779, 0.785398, 2.35619, 3.92699}; static constexpr double s_barrelRhoList[12] = {0.0890601, 0.0889871, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, diff --git a/src/Detectors/CMakeLists.txt b/src/Detectors/CMakeLists.txt index c3c5232..3c949eb 100644 --- a/src/Detectors/CMakeLists.txt +++ b/src/Detectors/CMakeLists.txt @@ -16,8 +16,8 @@ target_sources(Detectors PUBLIC SabreDetector.h SabreEfficiency.cpp SabreEfficiency.h - StripDetector.cpp - StripDetector.h + SX3Detector.cpp + SX3Detector.h ) target_link_libraries(Detectors diff --git a/src/Detectors/QQQDetector.cpp b/src/Detectors/QQQDetector.cpp index 11d5155..0f7046a 100644 --- a/src/Detectors/QQQDetector.cpp +++ b/src/Detectors/QQQDetector.cpp @@ -1,11 +1,8 @@ #include "QQQDetector.h" -QQQDetector::QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x, double y) : - m_innerR(R_in), m_outerR(R_out), m_deltaPhi(deltaPhi), m_centralPhi(phiCentral), m_translation(x,y,z), m_norm(0.0,0.0,1.0), - m_uniformFraction(0.0, 1.0), m_isSmearing(false) +QQQDetector::QQQDetector(double phiCentral, double zOffset, double xOffset, double yOffset) : + m_centralPhi(phiCentral), m_translation(xOffset,yOffset,zOffset), m_norm(0.0,0.0,1.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false) { - m_deltaR = (m_outerR - m_innerR)/s_nRings; - m_deltaPhiWedge = m_deltaPhi/s_nWedges; m_zRotation.SetAngle(m_centralPhi); m_ringCoords.resize(s_nRings); m_wedgeCoords.resize(s_nWedges); @@ -28,23 +25,23 @@ void QQQDetector::CalculateCorners() //Generate flat ring corner coordinates for(int i=0; i> m_ringCoords, m_wedgeCoords; @@ -60,6 +55,11 @@ private: static constexpr int s_nRings = 16; static constexpr int s_nWedges = 16; static constexpr double s_deg2rad = M_PI/180.0; + static constexpr double s_innerR = 0.0501; + static constexpr double s_outerR = 0.0990; + static constexpr double s_deltaR = (s_outerR - s_innerR) / s_nRings; + static constexpr double s_deltaPhiTotal = 1.52119; + static constexpr double s_deltaPhi = s_deltaPhiTotal / s_nWedges; }; #endif \ No newline at end of file diff --git a/src/Detectors/StripDetector.cpp b/src/Detectors/SX3Detector.cpp similarity index 70% rename from src/Detectors/StripDetector.cpp rename to src/Detectors/SX3Detector.cpp index 126e90f..47688ae 100644 --- a/src/Detectors/StripDetector.cpp +++ b/src/Detectors/SX3Detector.cpp @@ -1,5 +1,4 @@ -#include "StripDetector.h" -#include +#include "SX3Detector.h" /* Corner layout for each strip in the un-rotated frame @@ -10,36 +9,22 @@ | | | | x 2--------------------------3 z<--------X - | - | - | - y + | + | + | + y */ -StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) : - m_norm(1.0,0.0,0.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false) +StripDetector::StripDetector(double centerPhi, double centerZ, double centerRho) : + m_centerPhi(centerPhi), m_centerZ(centerZ), m_centerRho(centerRho), m_norm(1.0,0.0,0.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false) { - - m_nStrips = ns; - - m_stripLength = std::fabs(len); - m_totalWidth = std::fabs(wid); - m_frontStripWidth = m_totalWidth/m_nStrips; - m_backStripLength = m_stripLength/m_nStrips; - - while (cphi < 0) - cphi += 2*M_PI; - m_centerPhi = cphi; - m_centerZ = cz; - m_centerRho = std::fabs(crho); - m_zRotation.SetAngle(m_centerPhi); - m_frontStripCoords.resize(m_nStrips); - m_backStripCoords.resize(m_nStrips); - m_rotFrontStripCoords.resize(m_nStrips); - m_rotBackStripCoords.resize(m_nStrips); - for(int i=0; i det_max_phi) @@ -138,18 +123,18 @@ StripHit StripDetector::GetChannelRatio(double theta, double phi) double yhit = xhit*tan(phi); double zhit = sqrt(xhit*xhit+yhit*yhit)/tan(theta); - for (int s=0; s=m_frontStripCoords[s][0].X() && xhit <=m_frontStripCoords[s][0].X() && //Check min and max x (constant in flat) yhit >=m_frontStripCoords[s][1].Y() && yhit <=m_frontStripCoords[s][2].Y() && //Check min and max y zhit >=m_frontStripCoords[s][1].Z() && zhit <=m_frontStripCoords[s][0].Z()) //Check min and max z { hit.front_strip_index = s; - hit.front_ratio = (zhit-m_centerZ)/(m_stripLength/2); + hit.front_ratio = (zhit-m_centerZ)/(s_totalLength/2); break; } } - for (int s=0; s= m_backStripCoords[s][0].X() && xhit <= m_backStripCoords[s][0].X() && //Check min and max x (constant in flat) yhit >= m_backStripCoords[s][1].Y() && yhit <= m_backStripCoords[s][2].Y() && //Check min and max y zhit >= m_backStripCoords[s][1].Z() && zhit <= m_backStripCoords[s][0].Z()) //Check min and max z diff --git a/src/Detectors/StripDetector.h b/src/Detectors/SX3Detector.h similarity index 76% rename from src/Detectors/StripDetector.h rename to src/Detectors/SX3Detector.h index cc6463a..965e488 100644 --- a/src/Detectors/StripDetector.h +++ b/src/Detectors/SX3Detector.h @@ -1,5 +1,5 @@ -#ifndef STRIP_DETECTOR_H -#define STRIP_DETECTOR_H +#ifndef SX3_DETECTOR_H +#define SX3_DETECTOR_H // +z is along beam axis // +y is vertically "downward" in the lab frame @@ -30,7 +30,7 @@ class StripDetector { public: - StripDetector(int ns, double len, double wid, double cphi, double cz, double crho); + StripDetector(double centerPhi, double centerZ, double centerRho); ~StripDetector(); const ROOT::Math::XYZPoint& GetFrontStripCoordinates(int stripch, int corner) const { return m_frontStripCoords[stripch][corner]; } const ROOT::Math::XYZPoint& GetBackStripCoordinates(int stripch, int corner) const { return m_backStripCoords[stripch][corner]; } @@ -50,18 +50,10 @@ public: StripHit GetChannelRatio(double theta, double phi); private: - bool ValidChannel(int f) { return ((f >= 0 && f < m_nStrips) ? true : false); }; + bool ValidChannel(int f) { return ((f >= 0 && f < s_nStrips) ? true : false); }; bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); }; void CalculateCorners(); - int m_nStrips; - static constexpr int s_nCorners = 4; - - double m_stripLength; //common to all strips, hence total - double m_totalWidth; - double m_frontStripWidth; //assuming equal widths - double m_backStripLength; //assuming equal widths - double m_centerPhi; //assuming det centered above x-axis (corresponds to zero phi) double m_centerZ; double m_centerRho; //perpendicular radius from axis @@ -77,6 +69,15 @@ private: bool m_isSmearing; + //Units in meters + static constexpr double s_nStrips = 4; //Same for front and back + static constexpr double s_nCorners = 4; + static constexpr double s_totalLength = 0.075; //length of front strips + static constexpr double s_backStripLength = s_totalLength / s_nStrips; + static constexpr double s_totalWidth = 0.04; //width of back strips + static constexpr double s_frontStripWidth = s_totalWidth / s_nStrips; + static constexpr double s_deg2rad = M_PI/180.0; + }; #endif diff --git a/src/Detectors/SabreDetector.cpp b/src/Detectors/SabreDetector.cpp index 9d9da27..cf143ee 100644 --- a/src/Detectors/SabreDetector.cpp +++ b/src/Detectors/SabreDetector.cpp @@ -53,8 +53,7 @@ #include "SabreDetector.h" SabreDetector::SabreDetector() : -m_outerR(0.1351), m_innerR(0.0326), m_deltaPhiFlat(54.4*s_deg2rad), m_centerPhi(0.0), m_tilt(0.0), - m_translation(0.,0.,0.), m_norm(0,0,1.0), m_detectorID(-1) +m_centerPhi(0.0), m_tilt(0.0), m_translation(0.,0.,0.), m_norm(0,0,1.0), m_detectorID(-1) { m_yRotation.SetAngle(m_tilt); m_zRotation.SetAngle(m_centerPhi); @@ -73,17 +72,11 @@ m_outerR(0.1351), m_innerR(0.0326), m_deltaPhiFlat(54.4*s_deg2rad), m_centerPhi( m_tiltWedgeCoords[i].resize(4); } - m_deltaRFlat = m_outerR - m_innerR; - m_deltaRFlatRing = m_deltaRFlat/s_nRings; - m_deltaPhiFlatWedge = m_deltaPhiFlat/s_nWedges; - CalculateCorners(); } -SabreDetector::SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral, - double tiltFromVert, double zdist, double xdist, double ydist) : - m_outerR(Rout), m_innerR(Rin), m_deltaPhiFlat(deltaPhi_flat), m_centerPhi(phiCentral), m_tilt(tiltFromVert), - m_translation(xdist, ydist, zdist), m_norm(0,0,1.0), m_detectorID(detID) +SabreDetector::SabreDetector(int detID, double phiCentral, double tiltFromVert, double zdist, double xdist, double ydist) : + m_centerPhi(phiCentral), m_tilt(tiltFromVert), m_translation(xdist, ydist, zdist), m_norm(0,0,1.0), m_detectorID(detID) { m_yRotation.SetAngle(m_tilt); m_zRotation.SetAngle(m_centerPhi); @@ -104,10 +97,6 @@ SabreDetector::SabreDetector(int detID, double Rin, double Rout, double deltaPhi m_tiltWedgeCoords[i].resize(4); } - m_deltaRFlat = m_outerR - m_innerR; - m_deltaRFlatRing = m_deltaRFlat/s_nRings; - m_deltaPhiFlatWedge = m_deltaPhiFlat/s_nWedges; - CalculateCorners(); } @@ -123,23 +112,23 @@ void SabreDetector::CalculateCorners() //Generate flat ring corner coordinates for(int i=0; i m_innerR && r < m_outerR) || CheckPositionEqual(r, m_innerR) || CheckPositionEqual(r, m_outerR)) + double phi_1 = s_deltaPhiTotal/2.0; + double phi_2 = M_PI*2.0 - s_deltaPhiTotal/2.0; + return (((r > s_innerR && r < s_outerR) || CheckPositionEqual(r, s_innerR) || CheckPositionEqual(r, s_outerR)) && (phi > phi_2 || phi < phi_1 || CheckAngleEqual(phi, phi_1) || CheckAngleEqual(phi, phi_2))); }; @@ -121,57 +118,43 @@ private: */ bool IsRing(double r, int ringch) { - double ringtop = m_innerR + m_deltaRFlatRing*(ringch + 1); - double ringbottom = m_innerR + m_deltaRFlatRing*(ringch); + double ringtop = s_innerR + s_deltaR*(ringch + 1); + double ringbottom = s_innerR + s_deltaR*(ringch); return (r>ringbottom && rwedgebottom && phi> m_flatRingCoords, m_flatWedgeCoords; std::vector> m_tiltRingCoords, m_tiltWedgeCoords; + /*Class constants*/ + static constexpr double s_deg2rad = M_PI/180.0; + static constexpr int s_nRings = 16; + static constexpr int s_nWedges = 8; + static constexpr double s_outerR = 0.1351; + static constexpr double s_innerR = 0.0326; + static constexpr double s_deltaR = (s_outerR - s_innerR) / s_nRings; + static constexpr double s_deltaPhiTotal = 54.4 * s_deg2rad; + static constexpr double s_deltaPhi = (s_deltaPhiTotal / s_nWedges); + /*These are implicitly the width of the spacing between detector active strips*/ + static constexpr double s_positionTol = 0.0001; //0.1 mm position tolerance + static constexpr double s_angularTol = 0.1*M_PI/180.0; // 0.1 degree angular tolerance }; diff --git a/src/Detectors/SabreEfficiency.cpp b/src/Detectors/SabreEfficiency.cpp index 5d0a7b2..0ce0619 100644 --- a/src/Detectors/SabreEfficiency.cpp +++ b/src/Detectors/SabreEfficiency.cpp @@ -13,7 +13,7 @@ SabreEfficiency::SabreEfficiency() : { for(int i=0; i Date: Tue, 30 Aug 2022 15:42:49 -0400 Subject: [PATCH 4/6] Fix typos in SX3 class. Change name of array classes to array rather than efficiency. Update cmake lists. --- .../{AnasenEfficiency.cpp => AnasenArray.cpp} | 28 +++++++++---------- .../{AnasenEfficiency.h => AnasenArray.h} | 21 +++++++------- src/Detectors/CMakeLists.txt | 10 +++---- .../{DetectorEfficiency.h => DetectorArray.h} | 12 ++++---- src/Detectors/SX3Detector.cpp | 12 ++++---- src/Detectors/SX3Detector.h | 10 +++---- .../{SabreEfficiency.cpp => SabreArray.cpp} | 18 ++++++------ .../{SabreEfficiency.h => SabreArray.h} | 12 ++++---- src/Detectors/main.cpp | 8 +++--- 9 files changed, 66 insertions(+), 65 deletions(-) rename src/Detectors/{AnasenEfficiency.cpp => AnasenArray.cpp} (94%) rename src/Detectors/{AnasenEfficiency.h => AnasenArray.h} (81%) rename src/Detectors/{DetectorEfficiency.h => DetectorArray.h} (67%) rename src/Detectors/{SabreEfficiency.cpp => SabreArray.cpp} (94%) rename src/Detectors/{SabreEfficiency.h => SabreArray.h} (89%) diff --git a/src/Detectors/AnasenEfficiency.cpp b/src/Detectors/AnasenArray.cpp similarity index 94% rename from src/Detectors/AnasenEfficiency.cpp rename to src/Detectors/AnasenArray.cpp index 8c6cd19..7e31045 100644 --- a/src/Detectors/AnasenEfficiency.cpp +++ b/src/Detectors/AnasenArray.cpp @@ -1,4 +1,4 @@ -#include "AnasenEfficiency.h" +#include "AnasenArray.h" #include #include #include @@ -6,8 +6,8 @@ #include "TFile.h" #include "TTree.h" -AnasenEfficiency::AnasenEfficiency() : - DetectorEfficiency(), m_detectorEloss({14}, {28}, {1}, s_detectorThickness) +AnasenArray::AnasenArray() : + DetectorArray(), m_detectorEloss({14}, {28}, {1}, s_detectorThickness) { for(int i=0; i r1_points; std::vector r2_points; @@ -230,7 +230,7 @@ double AnasenEfficiency::RunConsistencyCheck() } -DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) +DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus) { DetectorResult observation; double thetaIncident; @@ -255,7 +255,7 @@ DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) return observation; } -DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) +DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus) { DetectorResult observation; double thetaIncident; @@ -280,7 +280,7 @@ DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) return observation; } -DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) +DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus) { DetectorResult observation; double thetaIncident; @@ -324,7 +324,7 @@ DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) return observation; } -DetectorResult AnasenEfficiency::IsAnasen(Mask::Nucleus& nucleus) +DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus) { DetectorResult result; if(nucleus.GetKE() <= s_energyThreshold) @@ -339,7 +339,7 @@ DetectorResult AnasenEfficiency::IsAnasen(Mask::Nucleus& nucleus) return result; } -void AnasenEfficiency::CountCoincidences(const std::vector& data, std::vector& counts) +void AnasenArray::CountCoincidences(const std::vector& data, std::vector& counts) { if (data.size() == 3 && data[1].isDetected && data[2].isDetected) { @@ -417,11 +417,11 @@ void AnasenEfficiency::CountCoincidences(const std::vector& data, } } -void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) { +void AnasenArray::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) { if(!dmap.IsValid()) { - std::cerr<<"Invalid Dead Channel Map at AnasenEfficiency::CalculateEfficiency()! If you want to run with all possible"; + std::cerr<<"Invalid Dead Channel Map at AnasenArray::CalculateEfficiency()! If you want to run with all possible"; std::cerr<<"channels active, simply load an empty file."<Close(); output->Close(); stats.close(); diff --git a/src/Detectors/AnasenEfficiency.h b/src/Detectors/AnasenArray.h similarity index 81% rename from src/Detectors/AnasenEfficiency.h rename to src/Detectors/AnasenArray.h index 8337793..a492d98 100644 --- a/src/Detectors/AnasenEfficiency.h +++ b/src/Detectors/AnasenArray.h @@ -1,23 +1,23 @@ -#ifndef ANASEN_EFFICIENCY_H -#define ANASEN_EFFICIENCY_H +#ifndef ANASEN_ARRAY_H +#define ANASEN_ARRAY_H #include -#include "DetectorEfficiency.h" +#include "DetectorArray.h" #include "SX3Detector.h" #include "QQQDetector.h" #include "Target.h" #include "Nucleus.h" #include "AnasenDeadChannelMap.h" -class AnasenEfficiency : public DetectorEfficiency +class AnasenArray : public DetectorArray { public: - AnasenEfficiency(); - ~AnasenEfficiency(); - void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override; - void DrawDetectorSystem(const std::string& filename) override; - double RunConsistencyCheck() override; + AnasenArray(); + ~AnasenArray(); + virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override; + virtual void DrawDetectorSystem(const std::string& filename) override; + virtual double RunConsistencyCheck() override; inline void SetDeadChannelMap(const std::string& filename) { dmap.LoadMapfile(filename); } private: @@ -27,7 +27,8 @@ private: DetectorResult IsAnasen(Mask::Nucleus& nucleus); void CountCoincidences(const std::vector& data, std::vector& counts); - std::vector m_Ring1, m_Ring2; + std::vector m_Ring1; + std::vector m_Ring2; std::vector m_forwardQQQs; std::vector m_backwardQQQs; diff --git a/src/Detectors/CMakeLists.txt b/src/Detectors/CMakeLists.txt index 3c949eb..12cb098 100644 --- a/src/Detectors/CMakeLists.txt +++ b/src/Detectors/CMakeLists.txt @@ -4,18 +4,18 @@ target_include_directories(Detectors PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_ target_sources(Detectors PUBLIC AnasenDeadChannelMap.cpp AnasenDeadChannelMap.h - AnasenEfficiency.cpp - AnasenEfficiency.h + AnasenArray.cpp + AnasenArray.h main.cpp - DetectorEfficiency.h + DetectorArray.h QQQDetector.cpp QQQDetector.h SabreDeadChannelMap.cpp SabreDeadChannelMap.h SabreDetector.cpp SabreDetector.h - SabreEfficiency.cpp - SabreEfficiency.h + SabreArray.cpp + SabreArray.h SX3Detector.cpp SX3Detector.h ) diff --git a/src/Detectors/DetectorEfficiency.h b/src/Detectors/DetectorArray.h similarity index 67% rename from src/Detectors/DetectorEfficiency.h rename to src/Detectors/DetectorArray.h index 7f3f292..5e7f37c 100644 --- a/src/Detectors/DetectorEfficiency.h +++ b/src/Detectors/DetectorArray.h @@ -1,5 +1,5 @@ -#ifndef DETECTOREFFICIENCY_H -#define DETECTOREFFICIENCY_H +#ifndef DETECTOR_ARRAY_H +#define DETECTOR_ARRAY_H #include #include @@ -14,18 +14,18 @@ struct DetectorResult std::string det_name = ""; }; -class DetectorEfficiency +class DetectorArray { public: - DetectorEfficiency() {}; - virtual ~DetectorEfficiency() {}; + DetectorArray() {}; + virtual ~DetectorArray() {}; virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) = 0; virtual void DrawDetectorSystem(const std::string& filename) = 0; virtual double RunConsistencyCheck() = 0; protected: - inline bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < s_epsilon ? true : false; }; + bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < s_epsilon ? true : false; }; static constexpr double s_epsilon = 1.0e-6; }; diff --git a/src/Detectors/SX3Detector.cpp b/src/Detectors/SX3Detector.cpp index 47688ae..c2b9c91 100644 --- a/src/Detectors/SX3Detector.cpp +++ b/src/Detectors/SX3Detector.cpp @@ -15,7 +15,7 @@ y */ -StripDetector::StripDetector(double centerPhi, double centerZ, double centerRho) : +SX3Detector::SX3Detector(double centerPhi, double centerZ, double centerRho) : m_centerPhi(centerPhi), m_centerZ(centerZ), m_centerRho(centerRho), m_norm(1.0,0.0,0.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false) { m_zRotation.SetAngle(m_centerPhi); @@ -35,9 +35,9 @@ StripDetector::StripDetector(double centerPhi, double centerZ, double centerRho) } -StripDetector::~StripDetector() {} +SX3Detector::~SX3Detector() {} -void StripDetector::CalculateCorners() +void SX3Detector::CalculateCorners() { double y_min, y_max, z_min, z_max; for (int s=0; s= 0 && f < s_nStrips) ? true : false); }; diff --git a/src/Detectors/SabreEfficiency.cpp b/src/Detectors/SabreArray.cpp similarity index 94% rename from src/Detectors/SabreEfficiency.cpp rename to src/Detectors/SabreArray.cpp index 0ce0619..ca5a35e 100644 --- a/src/Detectors/SabreEfficiency.cpp +++ b/src/Detectors/SabreArray.cpp @@ -1,4 +1,4 @@ -#include "SabreEfficiency.h" +#include "SabreArray.h" #include #include #include @@ -7,8 +7,8 @@ #include "TTree.h" -SabreEfficiency::SabreEfficiency() : - DetectorEfficiency(), m_deadlayerEloss({14}, {28}, {1}, s_deadlayerThickness), +SabreArray::SabreArray() : + DetectorArray(), m_deadlayerEloss({14}, {28}, {1}, s_deadlayerThickness), m_detectorEloss({14}, {28}, {1}, s_detectorThickness), m_degraderEloss({73}, {181}, {1}, s_degraderThickness) { @@ -29,9 +29,9 @@ SabreEfficiency::SabreEfficiency() : m_activeDetectors[4] = false; } -SabreEfficiency::~SabreEfficiency() {} +SabreArray::~SabreArray() {} -void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) +void SabreArray::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) { std::cout<<"----------SABRE Efficiency Calculation----------"<& data, std::vector& counts) +void SabreArray::CountCoincidences(const std::vector& data, std::vector& counts) { if (data.size() == 3 && data[1].isDetected && data[2].isDetected) { diff --git a/src/Detectors/SabreEfficiency.h b/src/Detectors/SabreArray.h similarity index 89% rename from src/Detectors/SabreEfficiency.h rename to src/Detectors/SabreArray.h index 7c43fb5..49eb313 100644 --- a/src/Detectors/SabreEfficiency.h +++ b/src/Detectors/SabreArray.h @@ -1,17 +1,17 @@ -#ifndef SABREEFFICIENCY_H -#define SABREEFFICIENCY_H +#ifndef SABRE_ARRAY_H +#define SABRE_ARRAY_H -#include "DetectorEfficiency.h" +#include "DetectorArray.h" #include "SabreDetector.h" #include "Target.h" #include "SabreDeadChannelMap.h" #include "Mask/Nucleus.h" -class SabreEfficiency : public DetectorEfficiency +class SabreArray : public DetectorArray { public: - SabreEfficiency(); - ~SabreEfficiency(); + SabreArray(); + ~SabreArray(); void SetDeadChannelMap(const std::string& filename) { m_deadMap.LoadMapfile(filename); }; void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override; void DrawDetectorSystem(const std::string& filename) override; diff --git a/src/Detectors/main.cpp b/src/Detectors/main.cpp index ba6f8f6..9ae677a 100644 --- a/src/Detectors/main.cpp +++ b/src/Detectors/main.cpp @@ -1,5 +1,5 @@ -#include "SabreEfficiency.h" -#include "AnasenEfficiency.h" +#include "SabreArray.h" +#include "AnasenArray.h" #include "KinematicsExceptions.h" #include #include @@ -24,7 +24,7 @@ int main(int argc, char** argv) std::string statsname = argv[3]; - SabreEfficiency sabre; + SabreArray sabre; std::string mapfile = "./etc/sabreDeadChannels_May2022.txt"; sabre.SetDeadChannelMap(mapfile); sabre.CalculateEfficiency(inputname, outputname, statsname); @@ -33,7 +33,7 @@ int main(int argc, char** argv) /* - AnasenEfficiency anasen; + AnasenArray anasen; std::string mapfile = "./etc/AnasenDeadChannels.txt"; anasen.SetDeadChannelMap(mapfile); anasen.CalculateEfficiency(inputname, outputname, statsname); From 96a9aef6857d74927fd7f7c208867ff5e6bb045c Mon Sep 17 00:00:00 2001 From: gwm17 Date: Wed, 31 Aug 2022 10:05:13 -0400 Subject: [PATCH 5/6] Fix a bug in SABRE where the detectors were tilted backwards. Cleanup some code --- src/Detectors/SabreArray.cpp | 4 +--- src/Detectors/SabreArray.h | 2 +- src/Plotters/RootPlotter.cpp | 5 ----- 3 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/Detectors/SabreArray.cpp b/src/Detectors/SabreArray.cpp index ca5a35e..9390b86 100644 --- a/src/Detectors/SabreArray.cpp +++ b/src/Detectors/SabreArray.cpp @@ -221,8 +221,6 @@ void SabreArray::DrawDetectorSystem(const std::string& filename) } } - output<<"SABRE Geometry File -- Coordinates for Detectors"< Date: Wed, 31 Aug 2022 11:51:47 -0400 Subject: [PATCH 6/6] Fix another bug in the SABRE geometry. Tilt needs to be flipped again in calc (minus implicitly handled in equation). --- src/Detectors/SabreArray.cpp | 1 + src/Detectors/SabreDetector.cpp | 20 +++++++++----------- src/Detectors/main.cpp | 2 +- 3 files changed, 11 insertions(+), 12 deletions(-) diff --git a/src/Detectors/SabreArray.cpp b/src/Detectors/SabreArray.cpp index 9390b86..01a41ae 100644 --- a/src/Detectors/SabreArray.cpp +++ b/src/Detectors/SabreArray.cpp @@ -248,6 +248,7 @@ double SabreArray::RunConsistencyCheck() if(channels.first != -1) { count++; + break; } } } diff --git a/src/Detectors/SabreDetector.cpp b/src/Detectors/SabreDetector.cpp index cf143ee..226a11b 100644 --- a/src/Detectors/SabreDetector.cpp +++ b/src/Detectors/SabreDetector.cpp @@ -194,8 +194,9 @@ ROOT::Math::XYZPoint SabreDetector::GetTrajectoryCoordinates(double theta, doubl if(m_translation.Vect().X() != 0.0 || m_translation.Vect().Y() != 0.0) return ROOT::Math::XYZPoint(); + double tilt = -1.0*m_tilt; //Calculate the *potential* phi in the flat detector - double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_centerPhi) - std::sin(m_centerPhi)*std::cos(phi)); + double phi_numerator = std::cos(tilt)*(std::sin(phi)*std::cos(m_centerPhi) - std::sin(m_centerPhi)*std::cos(phi)); double phi_denominator = std::cos(m_centerPhi)*std::cos(phi) + std::sin(m_centerPhi)*std::sin(phi); double phi_flat = std::atan2(phi_numerator, phi_denominator); if(phi_flat < 0) @@ -203,13 +204,13 @@ ROOT::Math::XYZPoint SabreDetector::GetTrajectoryCoordinates(double theta, doubl //Calculate the *potential* R in the flat detector double r_numerator = m_translation.Vect().Z()*std::cos(phi)*std::sin(theta); - double r_denominator = std::cos(phi_flat)*std::cos(m_centerPhi)*std::cos(m_tilt)*std::cos(theta) - + double r_denominator = std::cos(phi_flat)*std::cos(m_centerPhi)*std::cos(tilt)*std::cos(theta) - std::sin(phi_flat)*std::sin(m_centerPhi)*std::cos(theta) - - std::cos(phi_flat)*std::sin(m_tilt)*std::cos(phi)*std::sin(theta); + std::cos(phi_flat)*std::sin(tilt)*std::cos(phi)*std::sin(theta); double r_flat = r_numerator/r_denominator; //Calculate the distance from the origin to the hit on the detector - double R_to_detector = (r_flat*std::cos(phi_flat)*std::sin(m_tilt) + m_translation.Vect().Z())/std::cos(theta); + double R_to_detector = (r_flat*std::cos(phi_flat)*std::sin(tilt) + m_translation.Vect().Z())/std::cos(theta); double xhit = R_to_detector*std::sin(theta)*std::cos(phi); double yhit = R_to_detector*std::sin(theta)*std::sin(phi); double zhit = R_to_detector*std::cos(theta); @@ -244,8 +245,9 @@ std::pair SabreDetector::GetTrajectoryRingWedge(double theta, double p if(m_translation.Vect().X() != 0.0 || m_translation.Vect().Y() != 0.0) return std::make_pair(-1, -1); + double tilt = -1.0*m_tilt; //Calculate the *potential* phi in the flat detector - double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_centerPhi) - std::sin(m_centerPhi)*std::cos(phi)); + double phi_numerator = std::cos(tilt)*(std::sin(phi)*std::cos(m_centerPhi) - std::sin(m_centerPhi)*std::cos(phi)); double phi_denominator = std::cos(m_centerPhi)*std::cos(phi) + std::sin(m_centerPhi)*std::sin(phi); double phi_flat = std::atan2(phi_numerator, phi_denominator); if(phi_flat < 0) @@ -253,15 +255,11 @@ std::pair SabreDetector::GetTrajectoryRingWedge(double theta, double p //Calculate the *potential* R in the flat detector double r_numerator = m_translation.Vect().Z()*std::cos(phi)*std::sin(theta); - double r_denominator = std::cos(phi_flat)*std::cos(m_centerPhi)*std::cos(m_tilt)*std::cos(theta) - + double r_denominator = std::cos(phi_flat)*std::cos(m_centerPhi)*std::cos(tilt)*std::cos(theta) - std::sin(phi_flat)*std::sin(m_centerPhi)*std::cos(theta) - - std::cos(phi_flat)*std::sin(m_tilt)*std::cos(phi)*std::sin(theta); + std::cos(phi_flat)*std::sin(tilt)*std::cos(phi)*std::sin(theta); double r_flat = r_numerator/r_denominator; - //Calculate the distance from the origin to the hit on the detector - //double R_to_detector = (r_flat*std::cos(phi_flat)*std::sin(m_tilt) + m_translation.GetZ())/std::cos(theta); - - //Check to see if our flat coords fall inside the flat detector if(IsInside(r_flat, phi_flat)) { diff --git a/src/Detectors/main.cpp b/src/Detectors/main.cpp index 9ae677a..8cfb3b0 100644 --- a/src/Detectors/main.cpp +++ b/src/Detectors/main.cpp @@ -28,7 +28,7 @@ int main(int argc, char** argv) std::string mapfile = "./etc/sabreDeadChannels_May2022.txt"; sabre.SetDeadChannelMap(mapfile); sabre.CalculateEfficiency(inputname, outputname, statsname); - //std::cout<<"Running consistency check(1=success): "<