Merge pull request #2 from gwm17/devel

Devel

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)

README.md

@@ -1,8 +1,8 @@
-# MASK: Monte cArlo Simulation of Kinematics
+# 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 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 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 <kinematics_datafile> <new_detection_datafile> <new_detection_statsfile>`
+`./bin/Detectors <kinematics_datafile> <new_detection_datafile> <new_detection_statsfile>`
 
 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
 

src/Detectors/AnasenArray.cpp

@@ -1,4 +1,4 @@
-#include "AnasenEfficiency.h"
+#include "AnasenArray.h"
 #include <fstream>
 #include <iomanip>
 #include <iostream>
@@ -6,29 +6,29 @@
 #include "TFile.h"
 #include "TTree.h"
 
-AnasenEfficiency::AnasenEfficiency() :
+AnasenArray::AnasenArray() :
-	DetectorEfficiency(), m_detectorEloss({14}, {28}, {1}, s_detectorThickness)
+	DetectorArray(), m_detectorEloss({14}, {28}, {1}, s_detectorThickness)
 {
 	for(int i=0; i<s_nSX3PerBarrel; i++)
 	{
-		m_Ring1.emplace_back(4, s_sx3Length, s_sx3Width, s_barrelPhiList[i], s_barrel1Z, s_barrelRhoList[i]);
+		m_Ring1.emplace_back(s_barrelPhiList[i], s_barrel1Z, s_barrelRhoList[i]);
 		m_Ring1[i].SetPixelSmearing(true);
-		m_Ring2.emplace_back(4, s_sx3Length, s_sx3Width, s_barrelPhiList[i], s_barrel2Z, s_barrelRhoList[i]);
+		m_Ring2.emplace_back(s_barrelPhiList[i], s_barrel2Z, s_barrelRhoList[i]);
 		m_Ring2[i].SetPixelSmearing(true);
 	}
 	for(int i=0; i<s_nQQQ; i++)
 	{
-		m_forwardQQQs.emplace_back(s_qqqInnerR, s_qqqOuterR, s_qqqDeltaPhi, s_qqqPhiList[i], s_qqqZList[i]);
+		m_forwardQQQs.emplace_back(s_qqqPhiList[i], s_qqqZList[i]);
 		m_forwardQQQs[i].SetSmearing(true);
-		m_backwardQQQs.emplace_back(s_qqqInnerR, s_qqqOuterR, s_qqqDeltaPhi, s_qqqPhiList[i], (-1.0)*s_qqqZList[i]);
+		m_backwardQQQs.emplace_back(s_qqqPhiList[i], (-1.0)*s_qqqZList[i]);
 		m_backwardQQQs[i].SetSmearing(true);
 	}
 }
 
-AnasenEfficiency::~AnasenEfficiency() {}
+AnasenArray::~AnasenArray() {}
 
 
-void AnasenEfficiency::DrawDetectorSystem(const std::string& filename)
+void AnasenArray::DrawDetectorSystem(const std::string& filename)
 {
 	std::ofstream output(filename);
 
@@ -135,7 +135,7 @@ void AnasenEfficiency::DrawDetectorSystem(const std::string& filename)
 	output.close();
 }
 
-double AnasenEfficiency::RunConsistencyCheck()
+double AnasenArray::RunConsistencyCheck()
 {
 	std::vector<ROOT::Math::XYZPoint> r1_points;
 	std::vector<ROOT::Math::XYZPoint> 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<Mask::Nucleus>& data, std::vector<int>& counts)
+void AnasenArray::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
 {
 	if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
 	{
@@ -417,11 +417,11 @@ void AnasenEfficiency::CountCoincidences(const std::vector<Mask::Nucleus>& 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."<<std::endl;
 		return;
 	}
@@ -460,7 +460,7 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s
 		case 8: coinc_counts.resize(11, 0); break;
 		default:
 		{
-			std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
+			std::cerr<<"Bad reaction type at AnasenArray::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
 			input->Close();
 			output->Close();
 			stats.close();

src/Detectors/AnasenArray.h

@@ -1,23 +1,23 @@
-#ifndef ANASEN_EFFICIENCY_H
+#ifndef ANASEN_ARRAY_H
-#define ANASEN_EFFICIENCY_H
+#define ANASEN_ARRAY_H
 
 #include <string>
 
-#include "DetectorEfficiency.h"
+#include "DetectorArray.h"
-#include "StripDetector.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();
+	AnasenArray();
-	~AnasenEfficiency();
+	~AnasenArray();
-	void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
+	virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
-	void DrawDetectorSystem(const std::string& filename) override;
+	virtual void DrawDetectorSystem(const std::string& filename) override;
-	double RunConsistencyCheck() 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<Mask::Nucleus>& data, std::vector<int>& counts);
 
-	std::vector<StripDetector> m_Ring1, m_Ring2;
+	std::vector<SX3Detector> m_Ring1;
+	std::vector<SX3Detector> m_Ring2;
 	std::vector<QQQDetector> m_forwardQQQs;
 	std::vector<QQQDetector> 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,

src/Detectors/CMakeLists.txt

@@ -1,29 +1,27 @@
-add_executable(DetectEff)
+add_executable(Detectors)
-target_include_directories(DetectEff PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..)
+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
+    AnasenArray.cpp
-    AnasenEfficiency.h
+    AnasenArray.h
     main.cpp
-    DetectorEfficiency.h
+    DetectorArray.h
     QQQDetector.cpp
     QQQDetector.h
     SabreDeadChannelMap.cpp
     SabreDeadChannelMap.h
     SabreDetector.cpp
     SabreDetector.h
-    SabreEfficiency.cpp
+    SabreArray.cpp
-    SabreEfficiency.h
+    SabreArray.h
-    StripDetector.cpp
+    SX3Detector.cpp
-    StripDetector.h
+    SX3Detector.h
 )
 
-target_link_libraries(DetectEff
+target_link_libraries(Detectors
-    MaskDict
     Mask
-    catima
 )
 
-set_target_properties(DetectEff PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR})
+set_target_properties(Detectors PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR})

src/Detectors/DetectorArray.h

@@ -1,5 +1,5 @@
-#ifndef DETECTOREFFICIENCY_H
+#ifndef DETECTOR_ARRAY_H
-#define DETECTOREFFICIENCY_H
+#define DETECTOR_ARRAY_H
 
 #include <string>
 #include <cmath>
@@ -14,18 +14,18 @@ struct DetectorResult
 	std::string det_name = "";
 };
 
-class DetectorEfficiency
+class DetectorArray
 {
 public:
-	DetectorEfficiency() {};
+	DetectorArray() {};
-	virtual ~DetectorEfficiency() {};
+	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;
 };

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) :
+QQQDetector::QQQDetector(double phiCentral, double zOffset, double xOffset, double yOffset) :
-	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_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_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<s_nRings; i++)
 	{
-		x0 = (m_innerR + m_deltaR*(i+1))*std::cos(-m_deltaPhi/2.0);
+		x0 = (s_innerR + s_deltaR*(i+1))*std::cos(-s_deltaPhiTotal/2.0);
-		y0 = (m_innerR + m_deltaR*(i+1))*std::sin(-m_deltaPhi/2.0);
+		y0 = (s_innerR + s_deltaR*(i+1))*std::sin(-s_deltaPhiTotal/2.0);
 		z0 = 0.0;
 		m_ringCoords[i][0].SetXYZ(x0, y0, z0);
 
-		x1 = (m_innerR + m_deltaR*(i))*std::cos(-m_deltaPhi/2.0);
+		x1 = (s_innerR + s_deltaR*(i))*std::cos(-s_deltaPhiTotal/2.0);
-		y1 = (m_innerR + m_deltaR*(i))*std::sin(-m_deltaPhi/2.0);
+		y1 = (s_innerR + s_deltaR*(i))*std::sin(-s_deltaPhiTotal/2.0);
 		z1 = 0.0;
 		m_ringCoords[i][1].SetXYZ(x1, y1, z1);
 
-		x2 = (m_innerR + m_deltaR*(i))*std::cos(m_deltaPhi/2.0);
+		x2 = (s_innerR + s_deltaR*(i))*std::cos(s_deltaPhiTotal/2.0);
-		y2 = (m_innerR + m_deltaR*(i))*std::sin(m_deltaPhi/2.0);
+		y2 = (s_innerR + s_deltaR*(i))*std::sin(s_deltaPhiTotal/2.0);
 		z2 = 0.0;
 		m_ringCoords[i][2].SetXYZ(x2, y2, z2);
 
-		x3 = (m_innerR + m_deltaR*(i+1))*std::cos(m_deltaPhi/2.0);
+		x3 = (s_innerR + s_deltaR*(i+1))*std::cos(s_deltaPhiTotal/2.0);
-		y3 = (m_innerR + m_deltaR*(i+1))*std::sin(m_deltaPhi/2.0);
+		y3 = (s_innerR + s_deltaR*(i+1))*std::sin(s_deltaPhiTotal/2.0);
 		z3 = 0.0;
 		m_ringCoords[i][3].SetXYZ(x3, y3, z3);
 	}
@@ -52,23 +49,23 @@ void QQQDetector::CalculateCorners()
 	//Generate flat wedge corner coordinates
 	for(int i=0; i<s_nWedges; i++)
 	{
-		x0 = m_outerR * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
+		x0 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
-		y0 = m_outerR * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
+		y0 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
 		z0 = 0.0;
 		m_wedgeCoords[i][0].SetXYZ(x0, y0, z0);
 
-		x1 = m_innerR * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
+		x1 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
-		y1 = m_innerR * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
+		y1 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
 		z1 = 0.0;
 		m_wedgeCoords[i][1].SetXYZ(x1, y1, z1);
 
-		x2 = m_innerR * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
+		x2 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
-		y2 = m_innerR * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
+		y2 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
 		z2 = 0.0;
 		m_wedgeCoords[i][2].SetXYZ(x2, y2, z2);
 
-		x3 = m_outerR * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
+		x3 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
-		y3 = m_outerR * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
+		y3 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
 		z3 = 0.0;
 		m_wedgeCoords[i][3].SetXYZ(x3, y3, z3);
 	}
@@ -156,13 +153,13 @@ ROOT::Math::XYZPoint QQQDetector::GetHitCoordinates(int ringch, int wedgech)
 	double r_center, phi_center;
 	if(m_isSmearing)
 	{
-		r_center  = m_innerR + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+ringch)*m_deltaR;
+		r_center  = s_innerR + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+ringch)*s_deltaR;
-		phi_center = -m_deltaPhi/2.0 + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+wedgech)*m_deltaPhiWedge;
+		phi_center = -s_deltaPhiTotal/2.0 + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+wedgech)*s_deltaPhi;
 	}
 	else
 	{
-		r_center  = m_innerR + (0.5+ringch)*m_deltaR;
+		r_center  = s_innerR + (0.5+ringch)*s_deltaR;
-		phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhiWedge;
+		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);

src/Detectors/QQQDetector.h

@@ -20,7 +20,7 @@
 class QQQDetector
 {
 public:
-	QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x=0, double y=0);
+	QQQDetector(double phiCentral, double zOffset, double xOffset=0, double yOffset=0);
 	~QQQDetector();
 	const ROOT::Math::XYZPoint& GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; }
 	const ROOT::Math::XYZPoint& GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; }
@@ -42,11 +42,6 @@ private:
 	void CalculateCorners();
 	ROOT::Math::XYZPoint TransformCoordinates(ROOT::Math::XYZPoint& vector) { return m_translation * (m_zRotation * vector) ; }
 
-	double m_innerR;
-	double m_outerR;
-	double m_deltaR;
-	double m_deltaPhi;
-	double m_deltaPhiWedge;
 	double m_centralPhi;
 
 	std::vector<std::vector<ROOT::Math::XYZPoint>> 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

src/Detectors/SX3Detector.cpp

@@ -1,5 +1,4 @@
-#include "StripDetector.h"
+#include "SX3Detector.h"
-#include <iostream>
 
 /*
   Corner layout for each strip in the un-rotated frame
@@ -16,30 +15,16 @@
 											y
 */
 
-StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) :
+SX3Detector::SX3Detector(double centerPhi, double centerZ, double centerRho) :
-	m_norm(1.0,0.0,0.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false)
+	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_frontStripCoords.resize(s_nStrips);
-	m_backStripCoords.resize(m_nStrips);
+	m_backStripCoords.resize(s_nStrips);
-	m_rotFrontStripCoords.resize(m_nStrips);
+	m_rotFrontStripCoords.resize(s_nStrips);
-	m_rotBackStripCoords.resize(m_nStrips);
+	m_rotBackStripCoords.resize(s_nStrips);
-	for(int i=0; i<m_nStrips; i++)
+	for(int i=0; i<s_nStrips; i++)
 	{
 		m_frontStripCoords[i].resize(s_nCorners);
 		m_backStripCoords[i].resize(s_nCorners);
@@ -50,33 +35,33 @@ StripDetector::StripDetector(int ns, double len, double wid, double cphi, double
 	
 }
 
-StripDetector::~StripDetector() {}
+SX3Detector::~SX3Detector() {}
 
-void StripDetector::CalculateCorners()
+void SX3Detector::CalculateCorners()
 {
 	double y_min, y_max, z_min, z_max; 
-	for (int s=0; s<m_nStrips; s++)
+	for (int s=0; s<s_nStrips; s++)
 	{
-		y_max = m_totalWidth/2.0 - m_frontStripWidth*s;
+		y_max = s_totalWidth/2.0 - s_frontStripWidth*s;
-		y_min = m_totalWidth/2.0 - m_frontStripWidth*(s+1);
+		y_min = s_totalWidth/2.0 - s_frontStripWidth*(s+1);
-		z_max = m_centerZ + m_stripLength/2.0;
+		z_max = m_centerZ + s_totalLength/2.0;
-		z_min = m_centerZ - m_stripLength/2.0;
+		z_min = m_centerZ - s_totalLength/2.0;
 		m_frontStripCoords[s][2].SetXYZ(m_centerRho, y_max, z_max);
 		m_frontStripCoords[s][3].SetXYZ(m_centerRho, y_max, z_min);
 		m_frontStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max);
 		m_frontStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min);
 
-		z_max = (m_centerZ - m_stripLength/2.0) + (s+1)*m_backStripLength;
+		z_max = (m_centerZ - s_totalLength/2.0) + (s+1)*s_backStripLength;
-		z_min = (m_centerZ - m_stripLength/2.0) + (s)*m_backStripLength;
+		z_min = (m_centerZ - s_totalLength/2.0) + (s)*s_backStripLength;
-		y_max = m_totalWidth/2.0;
+		y_max = s_totalWidth/2.0;
-		y_min = -m_totalWidth/2.0;
+		y_min = -s_totalWidth/2.0;
 		m_backStripCoords[s][2].SetXYZ(m_centerRho, y_max, z_max);
 		m_backStripCoords[s][3].SetXYZ(m_centerRho, y_max, z_min);
 		m_backStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max);
 		m_backStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min);
 	}
 
-	for(int s=0; s<m_nStrips; s++)
+	for(int s=0; s<s_nStrips; s++)
 	{
 		m_rotFrontStripCoords[s][0] = m_zRotation*m_frontStripCoords[s][0];
 		m_rotFrontStripCoords[s][1] = m_zRotation*m_frontStripCoords[s][1];
@@ -90,7 +75,7 @@ void StripDetector::CalculateCorners()
 	}
 }
 
-ROOT::Math::XYZPoint StripDetector::GetHitCoordinates(int front_stripch, double front_strip_ratio)
+ROOT::Math::XYZPoint SX3Detector::GetHitCoordinates(int front_stripch, double front_strip_ratio)
 {
 
 	if (!ValidChannel(front_stripch) || !ValidRatio(front_strip_ratio))
@@ -98,22 +83,22 @@ ROOT::Math::XYZPoint StripDetector::GetHitCoordinates(int front_stripch, double
 
 	double y;
 	if(m_isSmearing)
-		y = -m_totalWidth/2.0 + (front_stripch + m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator()))*m_frontStripWidth;
+		y = -s_totalWidth/2.0 + (front_stripch + m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator()))*s_frontStripWidth;
 	else
-		y = -m_totalWidth/2.0 + (front_stripch+0.5)*m_frontStripWidth;
+		y = -s_totalWidth/2.0 + (front_stripch+0.5)*s_frontStripWidth;
 
 	//recall we're still assuming phi=0 det:
-	ROOT::Math::XYZPoint coords(m_centerRho, y, front_strip_ratio*(m_stripLength/2) + m_centerZ);
+	ROOT::Math::XYZPoint coords(m_centerRho, y, front_strip_ratio*(s_totalLength/2) + m_centerZ);
   
 	//NOW rotate by appropriate phi
 	return m_zRotation*coords;
 
 }
 
-StripHit StripDetector::GetChannelRatio(double theta, double phi)
+SX3Hit SX3Detector::GetChannelRatio(double theta, double phi)
 {
 
-	StripHit hit;
+	SX3Hit hit;
 	while (phi < 0)
 		phi += 2*M_PI;
 
@@ -125,7 +110,7 @@ StripHit StripDetector::GetChannelRatio(double theta, double phi)
 		phi -= 2*M_PI;
 
 	//then we can check easily whether it even hit the detector in phi
-	double det_max_phi = std::atan2(m_totalWidth/2,m_centerRho);
+	double det_max_phi = std::atan2(s_totalWidth/2,m_centerRho);
 	double det_min_phi = -det_max_phi;
   
 	if (phi < det_min_phi || phi > 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_nStrips; s++) {
+	for (int s=0; s<s_nStrips; s++) {
 		if (xhit >=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_nStrips; s++) {
+	for (int s=0; s<s_nStrips; s++) {
 		if (xhit >= 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

src/Detectors/SX3Detector.h

@@ -1,5 +1,5 @@
-#ifndef STRIP_DETECTOR_H
+#ifndef SX3_DETECTOR_H
-#define STRIP_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);
+	SX3Detector(double centerPhi, double centerZ, double centerRho);
-	~StripDetector();
+	~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

src/Detectors/SabreArray.cpp

@@ -1,4 +1,4 @@
-#include "SabreEfficiency.h"
+#include "SabreArray.h"
 #include <fstream>
 #include <iostream>
 #include <iomanip>
@@ -7,13 +7,13 @@
 #include "TTree.h"
 
 
-SabreEfficiency::SabreEfficiency() : 
+SabreArray::SabreArray() : 
-	DetectorEfficiency(), m_deadlayerEloss({14}, {28}, {1}, s_deadlayerThickness), 
+	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<s_nDets; i++)
-		m_detectors.emplace_back(i, s_innerR, s_outerR, s_deltaPhi*s_deg2rad, s_centerPhiList[i]*s_deg2rad, s_tilt*s_deg2rad, s_zOffset);
+		m_detectors.emplace_back(i, s_centerPhiList[i]*s_deg2rad, s_tilt*s_deg2rad, s_zOffset);
 	//Only 0, 1, 4 valid in degrader land
 	m_degradedDetectors[0] = true;
 	m_degradedDetectors[1] = true;
@@ -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----------"<<std::endl;
 	std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
@@ -186,7 +186,7 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
 	std::cout<<"---------------------------------------------"<<std::endl;
 }
 
-void SabreEfficiency::DrawDetectorSystem(const std::string& filename)
+void SabreArray::DrawDetectorSystem(const std::string& filename)
 {
 	std::ofstream output(filename);
 
@@ -221,8 +221,6 @@ void SabreEfficiency::DrawDetectorSystem(const std::string& filename)
  		}
  	}
 
- 	output<<"SABRE Geometry File -- Coordinates for Detectors"<<std::endl;
-	output<<"Edges: x y z"<<std::endl;
 	for(unsigned int i=0; i<ringxs.size(); i++)
 		output<<ringxs[i]<<" "<<ringys[i]<<" "<<ringzs[i]<<std::endl;
 	for(unsigned int i=0; i<wedgexs.size(); i++)
@@ -231,7 +229,7 @@ void SabreEfficiency::DrawDetectorSystem(const std::string& filename)
  	output.close();
 }
 
-double SabreEfficiency::RunConsistencyCheck()
+double SabreArray::RunConsistencyCheck()
 {
 	double theta, phi;
 	double npoints = 5.0*16.0*4.0;
@@ -250,6 +248,7 @@ double SabreEfficiency::RunConsistencyCheck()
  					if(channels.first != -1)
 					{
  						count++;
+						break;
  					}
  				}
  			}
@@ -260,7 +259,7 @@ double SabreEfficiency::RunConsistencyCheck()
 }
 
 /*Returns if detected, as well as total energy deposited in SABRE*/
-DetectorResult SabreEfficiency::IsSabre(Mask::Nucleus& nucleus)
+DetectorResult SabreArray::IsSabre(Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	if(nucleus.GetKE() <= s_energyThreshold)
@@ -279,7 +278,6 @@ DetectorResult SabreEfficiency::IsSabre(Mask::Nucleus& nucleus)
 		if(m_deadMap.IsDead(detector.GetDetectorID(), channel.first, 0) || m_deadMap.IsDead(detector.GetDetectorID(), channel.second, 1))
 				break; //dead channel check
 
-		observation.detectFlag = true;
 		observation.direction = detector.GetTrajectoryCoordinates(nucleus.vec4.Theta(), nucleus.vec4.Phi());
 		thetaIncident = std::acos(observation.direction.Dot(detector.GetNormTilted())/(observation.direction.R()));
 
@@ -293,6 +291,7 @@ DetectorResult SabreEfficiency::IsSabre(Mask::Nucleus& nucleus)
 
 		observation.det_name = "SABRE"+std::to_string(detector.GetDetectorID());
 		observation.energy_deposited = m_detectorEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, ke, M_PI - thetaIncident);
+		observation.detectFlag = true;
 		return observation;
 	}
 
@@ -300,7 +299,7 @@ DetectorResult SabreEfficiency::IsSabre(Mask::Nucleus& nucleus)
 	return observation;
 }
 
-void SabreEfficiency::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
+void SabreArray::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
 {
 	if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
 	{

src/Detectors/SabreArray.h

@@ -1,17 +1,17 @@
-#ifndef SABREEFFICIENCY_H
+#ifndef SABRE_ARRAY_H
-#define SABREEFFICIENCY_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();
+	SabreArray();
-	~SabreEfficiency();
+	~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_tilt = -40.0;
-    static constexpr double s_outerR = 0.1351;
-    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;

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_centerPhi(0.0), m_tilt(0.0), m_translation(0.,0.,0.), m_norm(0,0,1.0), m_detectorID(-1)
-	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, 
+SabreDetector::SabreDetector(int detID, double phiCentral, double tiltFromVert, double zdist, double xdist, double ydist) :
-							 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_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)
 {
 	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<s_nRings; i++)
 	{
-		x0 = (m_innerR + m_deltaRFlatRing*(i+1))*std::cos(-m_deltaPhiFlat/2.0);
+		x0 = (s_innerR + s_deltaR*(i+1))*std::cos(-s_deltaPhiTotal/2.0);
-		y0 = (m_innerR + m_deltaRFlatRing*(i+1))*std::sin(-m_deltaPhiFlat/2.0);
+		y0 = (s_innerR + s_deltaR*(i+1))*std::sin(-s_deltaPhiTotal/2.0);
 		z0 = 0.0;
 		m_flatRingCoords[i][0].SetXYZ(x0, y0, z0);
 
-		x1 = (m_innerR + m_deltaRFlatRing*(i))*std::cos(-m_deltaPhiFlat/2.0);
+		x1 = (s_innerR + s_deltaR*(i))*std::cos(-s_deltaPhiTotal/2.0);
-		y1 = (m_innerR + m_deltaRFlatRing*(i))*std::sin(-m_deltaPhiFlat/2.0);
+		y1 = (s_innerR + s_deltaR*(i))*std::sin(-s_deltaPhiTotal/2.0);
 		z1 = 0.0;
 		m_flatRingCoords[i][1].SetXYZ(x1, y1, z1);
 
-		x2 = (m_innerR + m_deltaRFlatRing*(i))*std::cos(m_deltaPhiFlat/2.0);
+		x2 = (s_innerR + s_deltaR*(i))*std::cos(s_deltaPhiTotal/2.0);
-		y2 = (m_innerR + m_deltaRFlatRing*(i))*std::sin(m_deltaPhiFlat/2.0);
+		y2 = (s_innerR + s_deltaR*(i))*std::sin(s_deltaPhiTotal/2.0);
 		z2 = 0.0;
 		m_flatRingCoords[i][2].SetXYZ(x2, y2, z2);
 
-		x3 = (m_innerR + m_deltaRFlatRing*(i+1))*std::cos(m_deltaPhiFlat/2.0);
+		x3 = (s_innerR + s_deltaR*(i+1))*std::cos(s_deltaPhiTotal/2.0);
-		y3 = (m_innerR + m_deltaRFlatRing*(i+1))*std::sin(m_deltaPhiFlat/2.0);
+		y3 = (s_innerR + s_deltaR*(i+1))*std::sin(s_deltaPhiTotal/2.0);
 		z3 = 0.0;
 		m_flatRingCoords[i][3].SetXYZ(x3, y3, z3);
 	}
@@ -147,23 +136,23 @@ void SabreDetector::CalculateCorners()
 	//Generate flat wedge corner coordinates
 	for(int i=0; i<s_nWedges; i++)
 	{
-		x0 = m_outerR * std::cos(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
+		x0 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
-		y0 = m_outerR * std::sin(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
+		y0 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
 		z0 = 0.0;
 		m_flatWedgeCoords[i][0].SetXYZ(x0, y0, z0);
 
-		x1 = m_innerR * std::cos(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
+		x1 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
-		y1 = m_innerR * std::sin(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
+		y1 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
 		z1 = 0.0;
 		m_flatWedgeCoords[i][1].SetXYZ(x1, y1, z1);
 
-		x2 = m_innerR * std::cos(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
+		x2 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
-		y2 = m_innerR * std::sin(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
+		y2 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
 		z2 = 0.0;
 		m_flatWedgeCoords[i][2].SetXYZ(x2, y2, z2);
 
-		x3 = m_outerR * std::cos(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
+		x3 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
-		y3 = m_outerR * std::sin(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
+		y3 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
 		z3 = 0.0;
 		m_flatWedgeCoords[i][3].SetXYZ(x3, y3, z3);
 	}
@@ -205,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)
@@ -214,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);
@@ -255,8 +245,9 @@ std::pair<int, int> 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<int, int> 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 r_center  = s_innerR + (0.5+ringch)*s_deltaR;
-	double phi_center = -m_deltaPhiFlat/2.0 + (0.5+wedgech)*m_deltaPhiFlatWedge;
+	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;

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_1 = s_deltaPhiTotal/2.0;
-		double phi_2 = M_PI*2.0 - m_deltaPhiFlat/2.0;
+		double phi_2 = M_PI*2.0 - s_deltaPhiTotal/2.0;
-		return (((r > m_innerR && r < m_outerR) || CheckPositionEqual(r, m_innerR) || CheckPositionEqual(r, m_outerR)) 
+		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 ringtop = s_innerR + s_deltaR*(ringch + 1);
-		double ringbottom = m_innerR + m_deltaRFlatRing*(ringch);
+		double ringbottom = s_innerR + s_deltaR*(ringch);
 		return (r>ringbottom && r<ringtop); 
 	};
 
 	inline bool IsRingTopEdge(double r, int ringch)
 	{
-		double ringtop = m_innerR + m_deltaRFlatRing*(ringch + 1);
+		double ringtop = s_innerR + s_deltaR*(ringch + 1);
 		return CheckPositionEqual(r, ringtop); 
 	};
 
 	inline bool IsRingBottomEdge(double r, int ringch)
 	{
-		double ringbottom = m_innerR + m_deltaRFlatRing*(ringch);
+		double ringbottom = s_innerR + s_deltaR*(ringch);
 		return CheckPositionEqual(r, ringbottom); 
 	};
 
 	inline bool IsWedge(double phi, int wedgech)
 	{
-		double wedgetop = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech+1);
+		double wedgetop = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech+1);
-		double wedgebottom = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech);
+		double wedgebottom = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech);
 		return ((phi>wedgebottom && phi<wedgetop));
 	};
 
 	inline bool IsWedgeTopEdge(double phi, int wedgech)
 	{
-		double wedgetop = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech+1);
+		double wedgetop = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech+1);
 		return CheckAngleEqual(phi, wedgetop);
 	}
 
 	inline bool IsWedgeBottomEdge(double phi, int wedgech)
 	{
-		double wedgebottom = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech);
+		double wedgebottom = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech);
 		return CheckAngleEqual(phi, wedgebottom);
 	}
 
-	/*Class constants*/
-	static constexpr int s_nRings = 16;
-	static constexpr int s_nWedges = 8;
-	static constexpr double s_deg2rad = M_PI/180.0;
-	/*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
-
 	/*Class data*/
-	double m_outerR;
-	double m_innerR;
-	double m_deltaRFlat;
-	double m_deltaRFlatRing;
-	double m_deltaPhiFlat;
-	double m_deltaPhiFlatWedge;
 	double m_centerPhi;
 	double m_tilt;
 	ROOT::Math::Translation3D m_translation;
@@ -183,6 +166,18 @@ private:
 	std::vector<std::vector<ROOT::Math::XYZPoint>> m_flatRingCoords, m_flatWedgeCoords;
 	std::vector<std::vector<ROOT::Math::XYZPoint>> 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
 };
 
 

src/Detectors/main.cpp

@@ -1,5 +1,5 @@
-#include "SabreEfficiency.h"
+#include "SabreArray.h"
-#include "AnasenEfficiency.h"
+#include "AnasenArray.h"
 #include "KinematicsExceptions.h"
 #include <iostream>
 #include <string>
@@ -13,21 +13,27 @@ int main(int argc, char** argv)
 		return 1;
 	}
 
+	if(!Mask::EnforceDictionaryLinked())
+	{
+		std::cerr<<"This should be illegal!"<<std::endl;
+		return 1;
+	}
+
 	std::string inputname = argv[1];
 	std::string outputname = argv[2];
 	std::string statsname = argv[3];
 
 	
-	SabreEfficiency sabre;
+	SabreArray sabre;
 	std::string mapfile = "./etc/sabreDeadChannels_May2022.txt";
 	sabre.SetDeadChannelMap(mapfile);
 	sabre.CalculateEfficiency(inputname, outputname, statsname);
-	//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
+	//std::cout<<"Running consistency check(0=success): "<<sabre.RunConsistencyCheck()<<std::endl;
 	//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt");
 	
 
 	/*
- 	AnasenEfficiency anasen;
+ 	AnasenArray anasen;
  	std::string mapfile = "./etc/AnasenDeadChannels.txt";
  	anasen.SetDeadChannelMap(mapfile);
  	anasen.CalculateEfficiency(inputname, outputname, statsname);

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})

src/MaskApp/CMakeLists.txt

@@ -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})

src/Plotters/CMakeLists.txt

@@ -12,7 +12,6 @@ target_sources(RootPlot PUBLIC
 )
 
 target_link_libraries(RootPlot
-    MaskDict
     Mask
     ${ROOT_LIBRARIES}
 )

src/Plotters/RootPlotter.cpp

@@ -13,11 +13,6 @@ static double FullPhi(double phi)
 RootPlotter::RootPlotter()
 {
 	TH1::AddDirectory(kFALSE);
-	//Enforce dictionary linking
-	if(Mask::EnforceDictionaryLinked())
-	{
-		std::cout<<"Dictionary Linked"<<std::endl;
-	}
 }
 
 RootPlotter::~RootPlotter() {}

src/Plotters/main.cpp

@@ -9,6 +9,12 @@ int main(int argc, char** argv)
         return 1;
     }
 
+    if(!Mask::EnforceDictionaryLinked())
+    {
+        std::cerr<<"This should be illegal!"<<std::endl;
+		return 1;
+    }
+
     RootPlotter plotter;
     plotter.Run(argv[1], argv[2]);
 }