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/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/AnasenEfficiency.cpp b/src/Detectors/AnasenArray.cpp similarity index 91% rename from src/Detectors/AnasenEfficiency.cpp rename to src/Detectors/AnasenArray.cpp index 5a0b7d8..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,29 +6,29 @@ #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 74% rename from src/Detectors/AnasenEfficiency.h rename to src/Detectors/AnasenArray.h index 9c2b7a6..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 "StripDetector.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; @@ -39,15 +40,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 0165619..12cb098 100644 --- a/src/Detectors/CMakeLists.txt +++ b/src/Detectors/CMakeLists.txt @@ -1,29 +1,27 @@ -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 - 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 - StripDetector.cpp - StripDetector.h + SabreArray.cpp + SabreArray.h + SX3Detector.cpp + SX3Detector.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/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/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 65% rename from src/Detectors/StripDetector.cpp rename to src/Detectors/SX3Detector.cpp index 126e90f..c2b9c91 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) +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_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 72% rename from src/Detectors/StripDetector.h rename to src/Detectors/SX3Detector.h index cc6463a..79a89b3 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 @@ -19,19 +19,19 @@ #include "Math/RotationZ.h" #include "Mask/RandomGenerator.h" -struct StripHit +struct SX3Hit { int front_strip_index=-1; int back_strip_index=-1; double front_ratio=0.0; }; -class StripDetector +class SX3Detector { public: - StripDetector(int ns, double len, double wid, double cphi, double cz, double crho); - ~StripDetector(); + SX3Detector(double centerPhi, double centerZ, double centerRho); + ~SX3Detector(); 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]; } const ROOT::Math::XYZPoint& GetRotatedFrontStripCoordinates(int stripch, int corner) const @@ -47,21 +47,13 @@ public: void SetPixelSmearing(bool isSmearing) { m_isSmearing = isSmearing; } ROOT::Math::XYZPoint GetHitCoordinates(int front_stripch, double front_strip_ratio); - StripHit GetChannelRatio(double theta, double phi); + SX3Hit 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/SabreEfficiency.cpp b/src/Detectors/SabreArray.cpp similarity index 92% rename from src/Detectors/SabreEfficiency.cpp rename to src/Detectors/SabreArray.cpp index 5d0a7b2..01a41ae 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,13 +7,13 @@ #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) { for(int i=0; i& 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 79% rename from src/Detectors/SabreEfficiency.h rename to src/Detectors/SabreArray.h index 193f0f9..2ae54bc 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; @@ -32,11 +32,8 @@ private: bool m_degradedDetectors[5]; //Sabre constants - static constexpr double s_innerR = 0.0326; - static constexpr double s_outerR = 0.1351; - static constexpr double s_tilt = 40.0; + static constexpr double s_tilt = -40.0; static constexpr double s_zOffset = -0.1245; - static constexpr double s_deltaPhi = 54.4; //delta phi for each det static constexpr int s_nDets = 5; static constexpr double s_centerPhiList[s_nDets] = {306.0, 18.0, 234.0, 162.0, 90.0}; static constexpr double s_deg2rad = M_PI/180.0; diff --git a/src/Detectors/SabreDetector.cpp b/src/Detectors/SabreDetector.cpp index 9d9da27..226a11b 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 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) @@ -264,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)) { @@ -322,8 +309,8 @@ ROOT::Math::XYZPoint SabreDetector::GetHitCoordinates(int ringch, int wedgech) if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech)) return ROOT::Math::XYZPoint(); - double r_center = m_innerR + (0.5+ringch)*m_deltaRFlatRing; - double phi_center = -m_deltaPhiFlat/2.0 + (0.5+wedgech)*m_deltaPhiFlatWedge; + double r_center = s_innerR + (0.5+ringch)*s_deltaR; + double phi_center = -s_deltaPhiTotal/2.0 + (0.5+wedgech)*s_deltaPhi; double x = r_center*std::cos(phi_center); double y = r_center*std::sin(phi_center); double z = 0; diff --git a/src/Detectors/SabreDetector.h b/src/Detectors/SabreDetector.h index b9b1d23..b0cde5c 100644 --- a/src/Detectors/SabreDetector.h +++ b/src/Detectors/SabreDetector.h @@ -65,7 +65,7 @@ class SabreDetector { public: SabreDetector(); - SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0); + SabreDetector(int detID, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0); ~SabreDetector(); /*Return coordinates of the corners of each ring/wedge in SABRE*/ @@ -84,9 +84,6 @@ public: int GetDetectorID() { return m_detectorID; } private: - - - void CalculateCorners(); /*Performs the transformation to the tilted,rotated,translated frame of the SABRE detector*/ @@ -109,9 +106,9 @@ private: /*Determine if a hit is within the bulk detector*/ bool IsInside(double r, double phi) { - double phi_1 = m_deltaPhiFlat/2.0; - double phi_2 = M_PI*2.0 - m_deltaPhiFlat/2.0; - return (((r > 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/main.cpp b/src/Detectors/main.cpp index b03a378..8cfb3b0 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 @@ -13,21 +13,27 @@ int main(int argc, char** argv) return 1; } + if(!Mask::EnforceDictionaryLinked()) + { + std::cerr<<"This should be illegal!"<