README updated. Added MASK namespace, commented math classes and SABRE classes. More to follow.

README.md

@@ -1 +1,41 @@
-#Kinematics
+# MASK: Monte cArlo Simulation of Kinematics
+MASK is a Monte Carlo simulation of reaction kinematics intended for use with the Super-Enge Split-pole Spectrograph (SESPS) at Florida State University.
+MASK is capable of simulating multi-step reaction-decay sequences, however in a purely kinematic sense, as it currently has no quantum mechanical input (this
+is planned to be added in the next version). It is also capable of testing detector efficiency; this version contains the methods necessary to simulate the efficiency
+of a reaction into the Silicion Array for Branching Ratio Detectors (SABRE).
+
+## Building MASK
+Download the repository from github using your favorite method. To build simply run
+
+`make`
+
+in the MASK directory, and the executable should be built and found in the bin directory.
+
+## Running MASK
+By default MASK is capable of simulating reactions of up to three steps. Here is a brief outline of each type:
+
+0. A reaction of type 0 is a pure decay. It is assumed isotropic by default; any other case will require the modification of the code.
+1. A reaction of type 1 is a pure reaction. It can incorporate all of the input file sampling parameters.
+2. A reaction of type 2 is a reaction followed by a subsequent decay of the residual nucleus. Again, all sampling is allowed.
+3. A reaction of type 3 is a reaction followed by a subsequent decay of the residual, followed by a decay of one of the products. Again, all sampling is allowed
+
+Note that for type 2 and type 3, the decays are assumed isotropic in the center-of-mass frame of the decay. 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 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.
+
+The default MASK program includes a calculation of SABRE efficiency, whose methods are contained in the SabreEfficiency and SabreDetector classes. This can be disabled by modifying the
+main.cpp file appropriately.
+
+In the input file the user also has the option to select to save the ROOT tree of the simulated data and the default plots. The options are yes or no. Yes saves them, no doesn't.
+
+To run MASK simply do the following from the MASK directory:
+
+`./bin/mask input.txt`
+
+Input.txt can be replaced by any text file with the correct format.
+
+## Requirements
+MASK requires that ROOT is installed for data writting and visualization, as well as for random number generation. Testing has been done only on ROOT 6. Mileage on all other ROOT versions
+will vary.

include/Kinematics.h

@@ -9,6 +9,7 @@
 #include <TFile.h>
 #include <TTree.h>
 
+namespace Mask {
 //For tree
 struct NucData {
 	double KE = -1;
@@ -60,4 +61,6 @@ private:
 	TRandom3* global_generator;
 };
 
+};
+
 #endif

include/LinkDef_Kinematics.h

@@ -1,5 +1,5 @@
 #ifdef __CLING__
 
-#pragma link C++ struct NucData+;
+#pragma link C++ struct Mask::NucData+;
 
 #endif

include/Nucleus.h

@@ -1,10 +1,12 @@
 #ifndef NUCLEUS_H
 #define NUCLEUS_H
 
-#include "G4Vec.h"
+#include "Vec4.h"
 #include <string>
 
-class Nucleus : public G4Vec {
+namespace Mask {
+
+class Nucleus : public Vec4 {
 public:
 	Nucleus();
 	Nucleus(int Z, int A);
@@ -39,4 +41,6 @@ private:
 
 };
 
+};
+
 #endif

include/Plotter.h

@@ -8,6 +8,8 @@
 #include <THashTable.h>
 #include "Nucleus.h"
 
+namespace Mask {
+
 struct GraphData {
 	std::string name;
 	std::string title;
@@ -42,4 +44,6 @@ private:
 
 };
 
+};
+
 #endif

include/Reaction.h

@@ -4,6 +4,8 @@
 #include "Nucleus.h"
 #include "LayeredTarget.h"
 
+namespace Mask {
+
 class Reaction {
 public:
 	Reaction();
@@ -51,4 +53,6 @@ private:
 
 };
 
+};
+
 #endif

include/ReactionSystem.h

@@ -5,6 +5,8 @@
 #include <vector>
 #include <TRandom3.h>
 
+namespace Mask {
+
 class ReactionSystem {
 public:
 	ReactionSystem();
@@ -22,7 +24,7 @@ public:
 	inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
 	inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
 	inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
-	inline const char* GetSystemEquation() { return m_sys_equation.c_str(); };
+	inline const char* GetSystemEquation() const { return m_sys_equation.c_str(); };
 
 protected:
 	virtual void LinkTarget();
@@ -39,4 +41,6 @@ protected:
 	static constexpr double deg2rad = M_PI/180.0;
 };
 
+};
+
 #endif

include/Rotation.h

@@ -1,21 +1,28 @@
-#ifndef GROTATION_H
+/*
-#define GROTATION_H
+	Classes which define rotations about the x, y, and z axes. Using these,
+	any arbitrary orientation can be described. Methods implemented for vector multiplication
+	as well as generating the inverse of the rotation.
+*/
+#ifndef ROTATION_H
+#define ROTATION_H
 
-#include "G3Vec.h"
+#include "Vec3.h"
 
-class GXRotation {
+namespace Mask {
+
+class XRotation {
 public:
-	GXRotation();
+	XRotation();
-	GXRotation(double ang);
+	XRotation(double ang);
-	~GXRotation();
+	~XRotation();
-	G3Vec Rotate(const G3Vec& vector);
+	Vec3 Rotate(const Vec3& vector);
 	inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix(); };
-	inline GXRotation GetInverse() { return GXRotation(-m_angle); };
+	inline XRotation GetInverse() { return XRotation(-m_angle); };
-	inline G3Vec operator*(const G3Vec& vector) {
+	inline Vec3 operator*(const Vec3& vector) {
 		double x = m_matrix[0][0]*vector[0] + m_matrix[0][1]*vector[1] + m_matrix[0][2]*vector[2];
 		double y = m_matrix[1][0]*vector[0] + m_matrix[1][1]*vector[1] + m_matrix[1][2]*vector[2];
 		double z = m_matrix[2][0]*vector[0] + m_matrix[2][1]*vector[1] + m_matrix[2][2]*vector[2];
-		return G3Vec(x, y, z);
+		return Vec3(x, y, z);
 	};
 
 private:
@@ -24,19 +31,19 @@ private:
 	double m_matrix[3][3];
 };
 
-class GYRotation {
+class YRotation {
 public:
-	GYRotation();
+	YRotation();
-	GYRotation(double ang);
+	YRotation(double ang);
-	~GYRotation();
+	~YRotation();
-	G3Vec Rotate(const G3Vec& vector);
+	Vec3 Rotate(const Vec3& vector);
 	inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix(); };
-	inline GYRotation GetInverse() { return GYRotation(-m_angle); };
+	inline YRotation GetInverse() { return YRotation(-m_angle); };
-	inline G3Vec operator*(const G3Vec& vector) {
+	inline Vec3 operator*(const Vec3& vector) {
 		double x = m_matrix[0][0]*vector[0] + m_matrix[0][1]*vector[1] + m_matrix[0][2]*vector[2];
 		double y = m_matrix[1][0]*vector[0] + m_matrix[1][1]*vector[1] + m_matrix[1][2]*vector[2];
 		double z = m_matrix[2][0]*vector[0] + m_matrix[2][1]*vector[1] + m_matrix[2][2]*vector[2];
-		return G3Vec(x, y, z);
+		return Vec3(x, y, z);
 	};
 
 private:
@@ -45,19 +52,19 @@ private:
 	double m_matrix[3][3];
 };
 
-class GZRotation {
+class ZRotation {
 public:
-	GZRotation();
+	ZRotation();
-	GZRotation(double ang);
+	ZRotation(double ang);
-	~GZRotation();
+	~ZRotation();
-	G3Vec Rotate(const G3Vec& vector);
+	Vec3 Rotate(const Vec3& vector);
 	inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix();};
-	inline GZRotation GetInverse() { return GZRotation(-m_angle); };
+	inline ZRotation GetInverse() { return ZRotation(-m_angle); };
-	inline G3Vec operator*(const G3Vec& vector) {
+	inline Vec3 operator*(const Vec3& vector) {
 		double x = m_matrix[0][0]*vector[0] + m_matrix[0][1]*vector[1] + m_matrix[0][2]*vector[2];
 		double y = m_matrix[1][0]*vector[0] + m_matrix[1][1]*vector[1] + m_matrix[1][2]*vector[2];
 		double z = m_matrix[2][0]*vector[0] + m_matrix[2][1]*vector[1] + m_matrix[2][2]*vector[2];
-		return G3Vec(x, y, z);
+		return Vec3(x, y, z);
 	};
 
 private:
@@ -66,4 +73,6 @@ private:
 	double m_matrix[3][3];
 };
 
+};
+
 #endif

include/SabreDetector.h

@@ -1,11 +1,62 @@
+/*
+
+  Class which represents a single MMM detector in the SABRE array at FSU. Origial code by KGH, re-written by
+  GWM.
+
+  Distances in meters, angles in radians.
+
+  The channel arrays have four points, one for each corner. The corners are
+  as follows, as if looking BACK along beam (i.e. from the target's pov):
+
+  0---------------------1
+  |                     |
+  |                     |      x
+  |                     |      <-----
+  |                     |      		|
+  |                     |      		|
+  3---------------------2      		y
+                               (z is hence positive along beam direction) 
+
+  The channel numbers, also as looking back from target pov, are:
+
+  >> rings are 0 -- 15 from inner to outer:
+
+    15 -------------------
+    14 -------------------
+    13 -------------------
+       .
+       .
+       .
+     2 -------------------
+     1 -------------------
+     0 -------------------
+
+  >> wedges are 0 -- 7 moving counterclockwise:
+
+      7 6 ... 1 0
+     | | |   | | |
+     | | |   | | |
+     | | |   | | |
+     | | |   | | |
+     | | |   | | |
+     | | |   | | |
+
+
+  >> Note that the detector starts centered on the x-axis (central phi = 0) untilted, and then is rotated to wherever the frick
+  	 it is supposed to go; phi = 90 is centered on y axis, pointing down towards the bottom of the scattering chamber
+
+  -- GWM, Dec 2020; based on the og code from kgh
+
+*/
+
 #ifndef SABREDETECTOR_H
 #define SABREDETECTOR_H
 
 #include <vector>
 #include <cmath>
 
-#include "G3Vec.h"
+#include "Vec3.h"
-#include "GRotation.h"
+#include "Rotation.h"
 
 class SabreDetector {
 public:
@@ -13,56 +64,108 @@ public:
 	SabreDetector();
 	SabreDetector(double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0);
 	~SabreDetector();
-	inline G3Vec GetRingFlatCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_flat[ch][corner] : G3Vec(); };
-	inline G3Vec GetWedgeFlatCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_flat[ch][corner] : G3Vec(); };
-	inline G3Vec GetRingTiltCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_tilt[ch][corner] : G3Vec(); };
-	inline G3Vec GetWedgeTiltCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_tilt[ch][corner] : G3Vec(); };
-	G3Vec GetTrajectoryCoordinates(double theta, double phi);
-	G3Vec GetHitCoordinates(int ringch, int wedgech);
 
+	/*Return coordinates of the corners of each ring/wedge in SABRE*/
+	inline Mask::Vec3 GetRingFlatCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_flat[ch][corner] : Mask::Vec3(); };
+	inline Mask::Vec3 GetWedgeFlatCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_flat[ch][corner] : Mask::Vec3(); };
+	inline Mask::Vec3 GetRingTiltCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_tilt[ch][corner] : Mask::Vec3(); };
+	inline Mask::Vec3 GetWedgeTiltCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_tilt[ch][corner] : Mask::Vec3(); };
+
+	Mask::Vec3 GetTrajectoryCoordinates(double theta, double phi);
+	std::pair<int, int> GetTrajectoryRingWedge(double theta, double phi);
+	Mask::Vec3 GetHitCoordinates(int ringch, int wedgech);
+
+	/*Basic getters*/
 	inline int GetNumberOfWedges() { return m_nWedges; };
 	inline int GetNumberOfRings() { return m_nRings; };
 	inline double GetInnerRadius() { return m_Rinner; };
 	inline double GetOuterRadius() { return m_Router; };
 	inline double GetPhiCentral() { return m_phiCentral; };
 	inline double GetTiltAngle() { return m_tilt; };
-	inline G3Vec GetTranslation() { return m_translation; };
+	inline Mask::Vec3 GetTranslation() { return m_translation; };
 
 
 private:
 
+	/*Class constants*/
 	static constexpr int m_nRings = 16;
 	static constexpr int m_nWedges = 8;
 	static constexpr double deg2rad = M_PI/180.0;
-	static constexpr double POSITION_TOL = 0.0001;
+	/*These are implicitly the width of the spacing between detector active strips*/
-	static constexpr double ANGULAR_TOL = M_PI/180.0;
+	static constexpr double POSITION_TOL = 0.0001; //0.1 mm position tolerance
+	static constexpr double ANGULAR_TOL = 0.1*M_PI/180.0; // 0.1 degree angular tolerance
 
 	void CalculateCorners();
-	inline G3Vec TransformToTiltedFrame(G3Vec& vector) { return (m_ZRot*(m_YRot*vector)) + m_translation; };
 
+	/*Performs the transformation to the tilted,rotated,translated frame of the SABRE detector*/
+	inline Mask::Vec3 TransformToTiltedFrame(Mask::Vec3& vector) { return (m_ZRot*(m_YRot*vector)) + m_translation; };
+
+	/*Determine if a given channel/corner combo is valid*/
 	inline bool CheckRingChannel(int ch) { return (ch<m_nRings && ch>=0) ? true : false; };
 	inline bool CheckWedgeChannel(int ch) { return (ch<m_nWedges && ch >=0) ? true : false; };
 	inline bool CheckCorner(int corner) { return (corner < 4 && corner >=0) ? true : false; };
 	inline bool CheckRingLocation(int ch, int corner) { return CheckRingChannel(ch) && CheckCorner(corner); };
 	inline bool CheckWedgeLocation(int ch, int corner) { return CheckWedgeChannel(ch) && CheckCorner(corner); };
 
+	/*
+		For all of the calculations, need a limit precision to determine if values are actually equal or not
+		Here the approx. size of the strip spacing is used as the precision.
+	*/
 	inline bool CheckPositionEqual(double val1,double val2) { return fabs(val1-val2) > POSITION_TOL ? false : true; };
 	inline bool CheckAngleEqual(double val1,double val2) { return fabs(val1-val2) > ANGULAR_TOL ? false : true; };
 
+	/*Determine if a hit is within the bulk detector*/
 	inline bool IsInside(double r, double phi) { 
 		double phi_1 = m_deltaPhi_flat/2.0;
 		double phi_2 = M_PI*2.0 - m_deltaPhi_flat/2.0;
 		return (((r > m_Rinner && r < m_Router) || CheckPositionEqual(r, m_Rinner) || CheckPositionEqual(r, m_Router)) && (phi > phi_2 || phi < phi_1 || CheckAngleEqual(phi, phi_1) || CheckAngleEqual(phi, phi_2)));
 	};
 
-	double m_Router, m_Rinner, m_deltaPhi_flat, m_phiCentral, m_tilt;
+	/*
-	G3Vec m_translation;
+		For a given radius/phi are you inside of a given ring/wedge channel,
-	GYRotation m_YRot;
+		or are you on the spacing between these channels
-	GZRotation m_ZRot;
+	*/
-	double m_deltaR_flat, m_deltaR_flat_ring;
+	inline bool IsRing(double r, int ringch) {
+		double ringtop = m_Rinner + m_deltaR_flat_ring*(ringch + 1);
+		double ringbottom = m_Rinner + m_deltaR_flat_ring*(ringch);
+		return (r>ringbottom && r<ringtop); 
+	};
 
-	std::vector<std::vector<G3Vec>> m_ringCoords_flat, m_wedgeCoords_flat;
+	inline bool IsRingTopEdge(double r, int ringch) {
-	std::vector<std::vector<G3Vec>> m_ringCoords_tilt, m_wedgeCoords_tilt;
+		double ringtop = m_Rinner + m_deltaR_flat_ring*(ringch + 1);
+		return CheckPositionEqual(r, ringtop); 
+	};
+
+	inline bool IsRingBottomEdge(double r, int ringch) {
+		double ringbottom = m_Rinner + m_deltaR_flat_ring*(ringch);
+		return CheckPositionEqual(r, ringbottom); 
+	};
+
+	inline bool IsWedge(double phi, int wedgech) {
+		double wedgetop = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech+1);
+		double wedgebottom = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech);
+		return ((phi>wedgebottom && phi<wedgetop));
+	};
+
+	inline bool IsWedgeTopEdge(double phi, int wedgech) {
+		double wedgetop = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech+1);
+		return CheckAngleEqual(phi, wedgetop);
+	}
+
+	inline bool IsWedgeBottomEdge(double phi, int wedgech) {
+		double wedgebottom = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech);
+		return CheckAngleEqual(phi, wedgebottom);
+	}
+
+	/*Class data*/
+	double m_Router, m_Rinner, m_deltaPhi_flat, m_phiCentral, m_tilt;
+	Mask::Vec3 m_translation;
+	Mask::YRotation m_YRot;
+	Mask::ZRotation m_ZRot;
+	double m_deltaR_flat, m_deltaR_flat_ring, m_deltaPhi_flat_wedge;
+
+	std::vector<std::vector<Mask::Vec3>> m_ringCoords_flat, m_wedgeCoords_flat;
+	std::vector<std::vector<Mask::Vec3>> m_ringCoords_tilt, m_wedgeCoords_tilt;
 
 };
 

include/ThreeStepSystem.h

@@ -3,6 +3,8 @@
 
 #include "ReactionSystem.h"
 
+namespace Mask {
+
 class ThreeStepSystem : public ReactionSystem {
 public:
 	ThreeStepSystem();
@@ -24,4 +26,6 @@ protected:
 
 };
 
+};
+
 #endif

include/TwoStepSystem.h

@@ -3,6 +3,8 @@
 
 #include "ReactionSystem.h"
 
+namespace Mask {
+
 class TwoStepSystem : public ReactionSystem {
 public:
 	TwoStepSystem();
@@ -22,4 +24,6 @@ private:
 
 };
 
+};
+
 #endif

include/Vec3.h

@@ -1,13 +1,21 @@
-#ifndef G3VEC_H
+/*
-#define G3VEC_H
+	Class to represent a 3-space vector in both cartesian and spherical coordinates. Can perform vector
+	addition, subtraction, and dot product. 
+
+	--GWM Dec 2020
+*/
+#ifndef VEC3_H
+#define VEC3_H
 
 #include <cmath>
 
-class G3Vec {
+namespace Mask {
+
+class Vec3 {
 public:
-	G3Vec();
+	Vec3();
-	G3Vec(double x, double y, double z);
+	Vec3(double x, double y, double z);
-	~G3Vec();
+	~Vec3();
 
 	void SetVectorCartesian(double x, double y, double z);
 	void SetVectorSpherical(double r, double theta, double phi);
@@ -24,18 +32,19 @@ public:
 	};
 
 	inline const double operator[](int index) const { return index>2 || index<0 ? 0.0 : m_data[index]; };
-	inline G3Vec& operator=(const G3Vec& rhs) { SetVectorCartesian(rhs.GetX(), rhs.GetY(), rhs.GetZ()); return *this; };
+	inline Vec3& operator=(const Vec3& rhs) { SetVectorCartesian(rhs.GetX(), rhs.GetY(), rhs.GetZ()); return *this; };
-	inline G3Vec operator+(const G3Vec& rhs) const { return G3Vec(this->GetX()+rhs.GetX(), this->GetY()+rhs.GetY(), this->GetZ()+rhs.GetZ()); };
+	inline Vec3 operator+(const Vec3& rhs) const { return Vec3(this->GetX()+rhs.GetX(), this->GetY()+rhs.GetY(), this->GetZ()+rhs.GetZ()); };
-	inline G3Vec operator-(const G3Vec& rhs) const { return G3Vec(this->GetX()-rhs.GetX(), this->GetY()-rhs.GetY(), this->GetZ()-rhs.GetZ()); };
+	inline Vec3 operator-(const Vec3& rhs) const { return Vec3(this->GetX()-rhs.GetX(), this->GetY()-rhs.GetY(), this->GetZ()-rhs.GetZ()); };
 
 
-	double Dot(const G3Vec& rhs) const;
+	double Dot(const Vec3& rhs) const;
-	G3Vec Cross(const G3Vec& rhs) const;
+	Vec3 Cross(const Vec3& rhs) const;
 
 
 
 private:
 
+	//Use instead of std::atan2. Better control over values close to x=0
 	inline double Atan2(double y, double x) const {
 		if(x != 0.0) return std::atan2(y, x);
 		else if(y > 0.0) return M_PI/2.0;
@@ -47,4 +56,6 @@ private:
 
 };
 
+};
+
 #endif

include/Vec4.h

@@ -1,15 +1,24 @@
-#ifndef G4VEC_H
+/*
-#define G4VEC_H
+	Class which represents a 4-momentum vector. Can perform vector addition, subtraction, dot product
+	and generate a boost vector to its rest frame as well as apply a boost to itself.
+
+	--GWM Dec 2020.
+*/
+#ifndef VEC4_H
+#define VEC4_H
 
 #include <cmath>
 
-class G4Vec {
+namespace Mask {
+
+class Vec4 {
 public:
-	G4Vec();
+	Vec4();
-	G4Vec(double px, double py, double pz, double E);
+	Vec4(double px, double py, double pz, double E);
-	virtual ~G4Vec();
+	virtual ~Vec4();
 	void SetVectorCartesian(double px, double py, double pz, double E);
 	void SetVectorSpherical(double theta, double phi, double p, double E);
+
 	inline double GetE() const {return m_data[3];};
 	inline double GetPx() const {return m_data[0];};
 	inline double GetPy() const {return m_data[1];};
@@ -22,31 +31,39 @@ public:
 		if(phi<0) phi += 2.0*M_PI;
 		return GetPx() == 0.0 && GetPy() == 0.0 ? 0.0 : phi;
 	};
+
 	inline double GetInvMass() const {return std::sqrt(GetE()*GetE() - GetP()*GetP());};
 	inline double GetKE() const {return GetE() - GetInvMass();};
 	inline const double* GetBoost() const {return &m_boost[0];};
+
 	void ApplyBoost(const double* boost);
 
 	//Only intended for use in looping access!
 	inline const double operator[] (int index) const {return index>3 || index < 0 ? 0.0 : m_data[index];};
 
-	inline G4Vec& operator=(const G4Vec& rhs) {SetVectorCartesian(rhs.GetPx(), rhs.GetPy(), rhs.GetPz(), rhs.GetE()); return *this;};
+	inline Vec4& operator=(const Vec4& rhs) {SetVectorCartesian(rhs.GetPx(), rhs.GetPy(), rhs.GetPz(), rhs.GetE()); return *this;};
-	inline G4Vec operator+(const G4Vec& rhs) const {return G4Vec(rhs.GetPx()+GetPx(), rhs.GetPy()+GetPy(), rhs.GetPz()+GetPz(), rhs.GetE()+GetE());};
+	inline Vec4 operator+(const Vec4& rhs) const {return Vec4(rhs.GetPx()+GetPx(), rhs.GetPy()+GetPy(), rhs.GetPz()+GetPz(), rhs.GetE()+GetE());};
-	inline G4Vec operator-(const G4Vec& rhs) const {return G4Vec(rhs.GetPx()-GetPx(), rhs.GetPy()-GetPy(), rhs.GetPz()-GetPz(), rhs.GetE()-GetE());};
+	inline Vec4 operator-(const Vec4& rhs) const {return Vec4(rhs.GetPx()-GetPx(), rhs.GetPy()-GetPy(), rhs.GetPz()-GetPz(), rhs.GetE()-GetE());};
-	double Dot(const G4Vec& rhs) const;
+	
-	G4Vec Cross(const G4Vec& rhs) const;
+	double Dot(const Vec4& rhs) const;
+	Vec4 Cross(const Vec4& rhs) const;
 
 private:
 	void CalcBoostToCM();
+
+	//use instead of std::atan2. Better controll over x=0
 	inline double Atan2(double y, double x) const { 
 		if(x != 0) return std::atan2(y, x);
 		else if( y > 0 ) return M_PI/2.0;
 		else if( y < 0 ) return -M_PI/2.0;
 		else return 0.0;
 	};
+
 	double m_data[4];
 	double m_boost[3];
 
 };
 
+};
+
 #endif

input.txt

@@ -3,29 +3,28 @@ OutputFile: /media/gordon/b6414c35-ec1f-4fc1-83bc-a6b68ca4325a/gwm17/test_newkin
 SaveTree: yes
 SavePlots: yes
 ----------Reaction Information----------
-ReactionType: 2
+ReactionType: 0
 Z	A (order is target, projectile, ejectile, break1, break3)
-5	10
+5	9
-2	3
+0	0
-2	4
+1	1
-1	2
 ----------Target Information----------
 Name: test_targ
 Layers: 2
 ~Layer1
-Thickness(ug/cm^2): 10
+Thickness(ug/cm^2): 0
 Z	A	Stoich
 6	12	1
 0
 ~
 ~Layer2
-Thickness(ug/cm^2): 80
+Thickness(ug/cm^2): 0
 Z	A	Stoich
-5	10	1
+5	9	1
 0
 ~
 ----------Sampling Information----------
-NumberOfSamples: 10000
+NumberOfSamples: 1000000
 BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001
 EjectileThetaMin(deg): 20.0 EjectileThetaMax(deg): 20.0
 ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.038

src/G3Vec.cpp

@@ -1,36 +0,0 @@
-#include "G3Vec.h"
-
-G3Vec::G3Vec() {
-	m_data[0] = 0.;
-	m_data[1] = 0.;
-	m_data[2] = 0.;
-}
-
-G3Vec::G3Vec(double x, double y, double z) {
-	m_data[0] = x;
-	m_data[1] = y;
-	m_data[2] = z;
-}
-
-G3Vec::~G3Vec() {}
-
-void G3Vec::SetVectorCartesian(double x, double y, double z) {
-	m_data[0] = x;
-	m_data[1] = y;
-	m_data[2] = z;
-}
-
-void G3Vec::SetVectorSpherical(double r, double theta, double phi) {
-	m_data[0] = r*std::cos(phi)*std::sin(theta);
-	m_data[1] = r*std::sin(phi)*std::sin(theta);
-	m_data[2] = r*std::cos(theta);
-}
-
-double G3Vec::Dot(const G3Vec& rhs) const {
-	return GetX()*rhs.GetX() + GetY()*rhs.GetY() + GetZ()*rhs.GetZ();
-}
-
-//Unimplemented
-G3Vec G3Vec::Cross(const G3Vec& rhs) const {
-	return G3Vec(0.,0.,0.);
-}

src/Kinematics.cpp

@@ -2,6 +2,8 @@
 #include <fstream>
 #include <iostream>
 
+namespace Mask {
+
 Kinematics::Kinematics() :
 	sys(nullptr), save_tree_flag(false), do_plotter_flag(false), global_generator(new TRandom3(0))
 {
@@ -392,3 +394,5 @@ void Kinematics::RunThreeStep() {
 	}
 	output->Close();
 }
+
+};

src/LayeredTarget.cpp

@@ -30,7 +30,7 @@ void LayeredTarget::AddLayer(std::vector<int>& Z, std::vector<int>& A, std::vect
 */
 double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle) {
 
-  if(rxnLayer < 0 || rxnLayer > layers.size()) {
+  if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) {
     std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<<std::endl;
     return 0.0;
   }
@@ -57,7 +57,7 @@ double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy
 */
 double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle) {
 
-  if(rxnLayer < 0 || rxnLayer > layers.size()) {
+  if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) {
     std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<<std::endl;
     return 0.0;
   }
@@ -65,7 +65,7 @@ double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy,
   double eloss = 0.0;
   double newEnergy = startEnergy;
   for(unsigned int i=rxnLayer; i<layers.size(); i++) {
-    if(i == rxnLayer) {
+    if(i == ((unsigned int)rxnLayer)) {
       eloss += layers[i].getEnergyLossHalf(ze, ae, newEnergy, angle);
       newEnergy = startEnergy - eloss;
     } else {
@@ -80,7 +80,7 @@ double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy,
 /*ReverseEnergyLoss version of GetEjectileEnergyLoss*/
 double LayeredTarget::GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle) {
 
-  if(rxnLayer < 0 || rxnLayer > layers.size()) {
+  if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) {
     std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<<std::endl;
     return 0.0;
   }

src/Nucleus.cpp

@@ -1,13 +1,15 @@
 #include "Nucleus.h"
 #include "MassLookup.h"
 
+namespace Mask {
+
 Nucleus::Nucleus () :
-	G4Vec(), m_z(0), m_a(0), m_gs_mass(0), m_symbol("")
+	Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_symbol("")
 {
 }
 
 Nucleus::Nucleus(int Z, int A) :
-	G4Vec(), m_z(Z), m_a(A)
+	Vec4(), m_z(Z), m_a(A)
 {
 	m_gs_mass = MASS.FindMass(Z, A);
 	m_symbol = MASS.FindSymbol(Z, A);
@@ -15,7 +17,7 @@ Nucleus::Nucleus(int Z, int A) :
 }
 
 Nucleus::Nucleus(int Z, int A, double px, double py, double pz, double E) :
-	G4Vec(px, py, pz, E), m_z(Z), m_a(A)
+	Vec4(px, py, pz, E), m_z(Z), m_a(A)
 {
 	m_gs_mass = MASS.FindMass(Z, A);
 	m_symbol = MASS.FindSymbol(Z, A);
@@ -33,3 +35,5 @@ bool Nucleus::SetIsotope(int Z, int A) {
 	SetVectorCartesian(0,0,0,m_gs_mass);
 	return true;
 }
+
+};

src/Plotter.cpp

@@ -1,6 +1,8 @@
 #include "Plotter.h"
 #include <iostream>
 
+namespace Mask {
+
 Plotter::Plotter() :
 	table(new THashTable())
 {
@@ -80,3 +82,5 @@ void Plotter::GenerateGraphs() {
 		garbage_collection.push_back(graph);
 	}
 }
+
+};

src/Reaction.cpp

@@ -1,6 +1,8 @@
 #include "Reaction.h"
 #include "KinematicsExceptions.h"
 
+namespace Mask {
+
 Reaction::Reaction() :
 	target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), nuc_initFlag(false), resid_elossFlag(false)
 {
@@ -153,7 +155,7 @@ void Reaction::CalculateDecay() {
 	}
 }
 
-
+};
 
 
 

src/ReactionSystem.cpp

@@ -1,6 +1,7 @@
 #include "ReactionSystem.h"
 #include "KinematicsExceptions.h"
-#include <iostream>
+
+namespace Mask {
 
 ReactionSystem::ReactionSystem() :
 	m_beamDist(0,0), m_theta1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("")
@@ -86,3 +87,5 @@ void ReactionSystem::RunSystem() {
 	}
 
 }
+
+};

src/Rotation.cpp

@@ -1,62 +1,71 @@
-#include "GRotation.h"
+/*
+	Classes which define rotations about the x, y, and z axes. Using these,
+	any arbitrary orientation can be described. Methods implemented for vector multiplication
+	as well as generating the inverse of the rotation.
+*/
+#include "Rotation.h"
 
-GXRotation::GXRotation() :
+namespace Mask {
+
+XRotation::XRotation() :
 m_angle(0)
 {
 	GenerateMatrix();
 }
 
-GXRotation::GXRotation(double angle) :
+XRotation::XRotation(double angle) :
 m_angle(angle)
 {
 	GenerateMatrix();
 }
 
-GXRotation::~GXRotation() {}
+XRotation::~XRotation() {}
 
-void GXRotation::GenerateMatrix() {
+void XRotation::GenerateMatrix() {
 	m_matrix[0][0] = 1.0; m_matrix[0][1] = 0.0; m_matrix[0][2] = 0.0;
 	m_matrix[1][0] = 0.0; m_matrix[1][1] = std::cos(m_angle); m_matrix[1][2] = -std::sin(m_angle);
 	m_matrix[2][0] = 0.0; m_matrix[2][1] = std::sin(m_angle); m_matrix[2][2] = std::cos(m_angle);
 }
 
-GYRotation::GYRotation() :
+YRotation::YRotation() :
 m_angle(0)
 {
 	GenerateMatrix();
 }
 
-GYRotation::GYRotation(double angle) :
+YRotation::YRotation(double angle) :
 m_angle(angle)
 {
 	GenerateMatrix();
 }
 
-GYRotation::~GYRotation() {}
+YRotation::~YRotation() {}
 
-void GYRotation::GenerateMatrix() {
+void YRotation::GenerateMatrix() {
 	m_matrix[0][0] = std::cos(m_angle); m_matrix[0][1] = 0.0; m_matrix[0][2] = -std::sin(m_angle);
 	m_matrix[1][0] = 0.0; m_matrix[1][1] = 1.0; m_matrix[1][2] = 0.0;
 	m_matrix[2][0] = std::sin(m_angle); m_matrix[2][1] = 0.0; m_matrix[2][2] = std::cos(m_angle);
 }
 
 
-GZRotation::GZRotation() :
+ZRotation::ZRotation() :
 m_angle(0)
 {
 	GenerateMatrix();
 }
 
-GZRotation::GZRotation(double angle) :
+ZRotation::ZRotation(double angle) :
 m_angle(angle)
 {
 	GenerateMatrix();
 }
 
-GZRotation::~GZRotation() {}
+ZRotation::~ZRotation() {}
 
-void GZRotation::GenerateMatrix() {
+void ZRotation::GenerateMatrix() {
 	m_matrix[0][0] = std::cos(m_angle); m_matrix[0][1] = -std::sin(m_angle); m_matrix[0][2] = 0.0;
 	m_matrix[1][0] = std::sin(m_angle); m_matrix[1][1] = std::cos(m_angle); m_matrix[1][2] = 0.0;
 	m_matrix[2][0] = 0.0; m_matrix[2][1] = 0.0; m_matrix[2][2] = 1.0;
 }
+
+};

src/SabreDetector.cpp

@@ -1,7 +1,8 @@
-#include "SabreDetector.h"
-#include <iostream>
-
 /*
+
+  Class which represents a single MMM detector in the SABRE array at FSU. Origial code by KGH, re-written by
+  GWM.
+
   Distances in meters, angles in radians.
 
   The channel arrays have four points, one for each corner. The corners are
@@ -48,6 +49,7 @@
 
 */
 
+#include "SabreDetector.h"
 
 SabreDetector::SabreDetector() :
 m_Router(0.1351), m_Rinner(0.0326), m_deltaPhi_flat(54.4*deg2rad), m_phiCentral(0.0), m_tilt(0.0), m_translation(0.,0.,0.)
@@ -55,6 +57,7 @@ m_Router(0.1351), m_Rinner(0.0326), m_deltaPhi_flat(54.4*deg2rad), m_phiCentral(
 	m_YRot.SetAngle(m_tilt);
 	m_ZRot.SetAngle(m_phiCentral);
 
+	//Initialize the coordinate arrays
 	m_ringCoords_flat.resize(m_nRings);
 	m_ringCoords_tilt.resize(m_nRings);
 	m_wedgeCoords_flat.resize(m_nWedges);
@@ -77,10 +80,10 @@ m_Router(0.1351), m_Rinner(0.0326), m_deltaPhi_flat(54.4*deg2rad), m_phiCentral(
 SabreDetector::SabreDetector(double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist, double ydist) :
 m_Router(Rout), m_Rinner(Rin), m_deltaPhi_flat(deltaPhi_flat), m_phiCentral(phiCentral), m_tilt(tiltFromVert), m_translation(xdist, ydist, zdist)
 {
-	std::cout<<m_tilt<<std::endl;
 	m_YRot.SetAngle(m_tilt);
 	m_ZRot.SetAngle(m_phiCentral);
 
+	//Initialize coordinate arrays
 	m_ringCoords_flat.resize(m_nRings);
 	m_ringCoords_tilt.resize(m_nRings);
 	m_wedgeCoords_flat.resize(m_nWedges);
@@ -94,6 +97,10 @@ m_Router(Rout), m_Rinner(Rin), m_deltaPhi_flat(deltaPhi_flat), m_phiCentral(phiC
 		m_wedgeCoords_tilt[i].resize(4);
 	}
 
+	m_deltaR_flat = m_Router - m_Rinner;
+	m_deltaR_flat_ring = m_deltaR_flat/m_nRings;
+	m_deltaPhi_flat_wedge = m_deltaPhi_flat/m_nWedges;
+
 	CalculateCorners();
 }
 
@@ -101,55 +108,52 @@ SabreDetector::~SabreDetector() {}
 
 void SabreDetector::CalculateCorners() {
 
-	double deltaPhi_per_wedge = m_deltaPhi_flat/double(m_nWedges);
-	double deltaR_per_ring = (m_Router - m_Rinner)/double(m_nRings);
-
 	double x0, x1, x2, x3;
 	double y0, y1, y2, y3;
 	double z0, z1, z2, z3;
 
-	//Generate flat rings
+	//Generate flat ring corner coordinates
 	for(int i=0; i<m_nRings; i++) {
-		x0 = (m_Rinner + deltaR_per_ring*(i+1))*std::cos(-m_deltaPhi_flat/2.0);
+		x0 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::cos(-m_deltaPhi_flat/2.0);
-		y0 = (m_Rinner + deltaR_per_ring*(i+1))*std::sin(-m_deltaPhi_flat/2.0);
+		y0 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::sin(-m_deltaPhi_flat/2.0);
 		z0 = 0.0;
 		m_ringCoords_flat[i][0].SetVectorCartesian(x0, y0, z0);
 
-		x1 = (m_Rinner + deltaR_per_ring*(i))*std::cos(-m_deltaPhi_flat/2.0);
+		x1 = (m_Rinner + m_deltaR_flat_ring*(i))*std::cos(-m_deltaPhi_flat/2.0);
-		y1 = (m_Rinner + deltaR_per_ring*(i))*std::sin(-m_deltaPhi_flat/2.0);
+		y1 = (m_Rinner + m_deltaR_flat_ring*(i))*std::sin(-m_deltaPhi_flat/2.0);
 		z1 = 0.0;
 		m_ringCoords_flat[i][1].SetVectorCartesian(x1, y1, z1);
 
-		x2 = (m_Rinner + deltaR_per_ring*(i))*std::cos(m_deltaPhi_flat/2.0);
+		x2 = (m_Rinner + m_deltaR_flat_ring*(i))*std::cos(m_deltaPhi_flat/2.0);
-		y2 = (m_Rinner + deltaR_per_ring*(i))*std::sin(m_deltaPhi_flat/2.0);
+		y2 = (m_Rinner + m_deltaR_flat_ring*(i))*std::sin(m_deltaPhi_flat/2.0);
 		z2 = 0.0;
 		m_ringCoords_flat[i][2].SetVectorCartesian(x2, y2, z2);
 
-		x3 = (m_Rinner + deltaR_per_ring*(i+1))*std::cos(m_deltaPhi_flat/2.0);
+		x3 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::cos(m_deltaPhi_flat/2.0);
-		y3 = (m_Rinner + deltaR_per_ring*(i+1))*std::sin(m_deltaPhi_flat/2.0);
+		y3 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::sin(m_deltaPhi_flat/2.0);
 		z3 = 0.0;
 		m_ringCoords_flat[i][3].SetVectorCartesian(x3, y3, z3);
 	}
 
-	//Generate flat wedges
+	//Generate flat wedge corner coordinates
 	for(int i=0; i<m_nWedges; i++) {
-		x0 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
+		x0 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
-		y0 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
+		y0 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
 		z0 = 0.0;
 		m_wedgeCoords_flat[i][0].SetVectorCartesian(x0, y0, z0);
 
-		x1 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
+		x1 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
-		y1 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
+		y1 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
 		z1 = 0.0;
 		m_wedgeCoords_flat[i][1].SetVectorCartesian(x1, y1, z1);
 
-		x2 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
+		x2 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
-		y2 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
+		y2 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
 		z2 = 0.0;
 		m_wedgeCoords_flat[i][2].SetVectorCartesian(x2, y2, z2);
 
-		x3 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
+		x3 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
-		y3 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
+		y3 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
 		z3 = 0.0;
 		m_wedgeCoords_flat[i][3].SetVectorCartesian(x3, y3, z3);
 	}
@@ -169,10 +173,26 @@ void SabreDetector::CalculateCorners() {
 	}
 }
 
-//Note: float used for calculations due to lack of precision from sin, cos, and tangent functions
+/*
-G3Vec SabreDetector::GetTrajectoryCoordinates(double theta, double phi) {
+	Given a unit vector (R=1, theta, phi) which corresponds to some particle's trajectory,
+	determine whether that particle will intersect with this SABRE detector. If it does calculate
+	the coordinates of the hit. The equation is as follows:
+	
+	Rz(eta)*Ry(psi)*Flat_vector(R', theta'=PI/2, phi') + translation = Tilted_vector(R, theta, phi)
+
+	Where Rz is the ZRotation, Ry is the YRotation, F_vector is the vector of the coordinates in the flat detector frame,
+	and Tilted_vector is the vector of the hit coordinates in the tilted frame. The theta and phi of the the Tilted_vector correspond
+	to the input arguments of the function.
+
+	It checks to deterime whether or not the particle hits within the borders (read: edges) of the SABRE detector, and does not account for
+	the spacing between rings and wedges.
+
+	!NOTE: This currently only applies to a configuration where there is no translation in x & y. The math becomes significantly messier in these cases.
+	Also, don't use tan(). It's behavior near PI/2 makes it basically useless for these.
+*/
+Mask::Vec3 SabreDetector::GetTrajectoryCoordinates(double theta, double phi) {
 	if(m_translation.GetX() != 0.0 || m_translation.GetY() != 0.0) {
-		return G3Vec();
+		return Mask::Vec3();
 	}
 
 	//Calculate the *potential* phi in the flat detector
@@ -195,27 +215,95 @@ G3Vec SabreDetector::GetTrajectoryCoordinates(double theta, double phi) {
 
 	//Check to see if our flat coords fall inside the flat detector
 	if(IsInside(r_flat, phi_flat)) {
-		return G3Vec(xhit, yhit, zhit);
+		return Mask::Vec3(xhit, yhit, zhit);
 	} else {
-		return G3Vec();
+		return Mask::Vec3();
 	}
 }
 
-G3Vec SabreDetector::GetHitCoordinates(int ringch, int wedgech) {
+/*
+	Given a unit vector (R=1, theta, phi) which corresponds to some particle's trajectory,
+	determine whether that particle will intersect with this SABRE detector. If it does determine 
+	which ring and wedge the hit occurs in. The equation is as follows:
+	
+	Rz(eta)*Rx(psi)*Flat_vector(R', theta'=PI/2, phi') + translation = Tilted_vector(R, theta, phi)
+
+	Where Rz is the ZRotation, Rx is the XRotation, F_vector is the vector of the coordinates in the flat detector frame,
+	and Tilted_vector is the vector of the hit coordinates in the tilted frame. The theta and phi of the the Tilted_vector correspond
+	to the input arguments of the function.
+
+	Then using the flat coordinate R' and phi' determine which ring/wedge channels are hit. For precision purposes, the channel is not calculated, but
+	rather found using comparisions. This method accounts for the spacing between rings and wedges.
+
+	!NOTE: This currently only applies to a configuration where there is no translation in x & y. The math becomes significantly messier in these cases.
+	Also, don't use tan(). It's behavior near PI/2 makes it basically useless for these.
+*/
+std::pair<int, int> SabreDetector::GetTrajectoryRingWedge(double theta, double phi) {
+	if(m_translation.GetX() != 0.0 || m_translation.GetY() != 0.0) {
+		return std::make_pair(-1, -1);
+	}
+
+	//Calculate the *potential* phi in the flat detector
+	double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_phiCentral) - std::sin(m_phiCentral)*std::cos(phi));
+	double phi_denominator = std::cos(m_phiCentral)*std::cos(phi) + std::sin(m_phiCentral)*std::sin(phi);
+	double phi_flat = std::atan2(phi_numerator, phi_denominator);
+	if(phi_flat < 0) phi_flat += M_PI*2.0;
+
+	//Calculate the *potential* R in the flat detector
+	double r_numerator = m_translation.GetZ()*std::cos(phi)*std::sin(theta);
+	double r_denominator = std::cos(phi_flat)*std::cos(m_phiCentral)*std::cos(m_tilt)*std::cos(theta) - std::sin(phi_flat)*std::sin(m_phiCentral)*std::cos(theta) - std::cos(phi_flat)*std::sin(m_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)) {
+		int ringch, wedgech;
+		if(phi_flat > M_PI) phi_flat -= 2.0*M_PI; //Need phi in terms of [-deltaPhi_flat/2, deltaPhi_flat/2]
+		for(int i=0; i<m_nRings; i++) {
+			if(IsRingTopEdge(r_flat, i) || IsRingBottomEdge(r_flat, i)) { //If it falls in the interstrip spacing, kill it
+				ringch = -1;
+				break;
+			} else if(IsRing(r_flat, i)) {
+				ringch = i;
+				break;
+			}
+		}
+		for(int i=0; i<m_nWedges; i++) {
+			if(IsWedgeTopEdge(phi_flat, i) || IsWedgeBottomEdge(phi_flat, i)){ //If it falls in the interstrip spacing, kill it
+				wedgech = -1;
+				break;
+			} else if(IsWedge(phi_flat, i)) {
+				wedgech = i;
+				break;
+			}
+		}
+		return std::make_pair(ringch, wedgech);
+	} else {
+		return std::make_pair(-1,-1);
+	}
+}
+
+/*
+	Given a ring/wedge of this SABRE detector, calculate the coordinates of a hit.
+	Currently gives a point in the *center* of the pixel. Better method would be to
+	randomly wiggle the point within the pixel. Method intended for use with data, or
+	to smear out simulated data to mimic real data.
+*/
+Mask::Vec3 SabreDetector::GetHitCoordinates(int ringch, int wedgech) {
 	if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech)) {
-		return G3Vec();
+		return Mask::Vec3();
 	}
 
-	double deltaR_per_ring = (m_Router - m_Rinner)/double(m_nRings);
+	double r_center  = m_Rinner + (0.5+ringch)*m_deltaR_flat_ring;
-	double deltaPhi_per_wedge = (m_deltaPhi_flat)/double(m_nWedges);
+	double phi_center = -m_deltaPhi_flat/2.0 + (0.5+wedgech)*m_deltaPhi_flat_wedge;
-
-	double r_center  = m_Rinner + (0.5+ringch)*deltaR_per_ring;
-	double phi_center = -m_deltaPhi_flat/2.0 + (0.5+wedgech)*deltaPhi_per_wedge;
 	double x = r_center*std::cos(phi_center);
 	double y = r_center*std::sin(phi_center);
 	double z = 0;
 
-	G3Vec hit(x, y, z);
+	Mask::Vec3 hit(x, y, z);
 
 	return TransformToTiltedFrame(hit);
 }

src/SabreEfficiency.cpp

@@ -18,7 +18,7 @@ SabreEfficiency::SabreEfficiency() :
  	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
  	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
 
- 	G3Vec coords;
+ 	Mask::Vec3 coords;
  	for(int i=0; i<5; i++) {
  		for(int j=0; j<detectors[i].GetNumberOfRings(); j++) {
  			for(int k=0; k<4; k++) {
@@ -50,17 +50,17 @@ void SabreEfficiency::CalculateEfficiency(const char* file) {
 	std::cout<<"Running efficiency calculation..."<<std::endl;
 
 	switch(m_rxn_type) {
-		case Kinematics::ONESTEP_DECAY:
+		case Mask::Kinematics::ONESTEP_DECAY:
 		{
 			RunDecay(file);
 			break;
 		}
-		case Kinematics::TWOSTEP:
+		case Mask::Kinematics::TWOSTEP:
 		{
 			Run2Step(file);
 			break;
 		}
-		case Kinematics::THREESTEP:
+		case Mask::Kinematics::THREESTEP:
 		{
 			Run3Step(file);
 			break;
@@ -75,8 +75,8 @@ void SabreEfficiency::RunDecay(const char* file) {
 	TFile* input = TFile::Open(file, "UPDATE");
 	TTree* tree = (TTree*) input->Get("DataTree");
 
-	NucData* eject = new NucData();
+	Mask::NucData* eject = new Mask::NucData();
-	NucData* resid = new NucData();
+	Mask::NucData* resid = new Mask::NucData();
 
 
 	tree->SetBranchAddress("ejectile", &eject);
@@ -92,7 +92,7 @@ void SabreEfficiency::RunDecay(const char* file) {
 	int count = 0;
 	int npercent = 0;
 
-	G3Vec coordinates;
+	Mask::Vec3 coordinates;
 
 	for(int i=0; i<tree->GetEntries(); i++) {
 		if(++count == percent5) {//Show progress every 5%
@@ -105,8 +105,9 @@ void SabreEfficiency::RunDecay(const char* file) {
 
 		if(eject->KE >= ENERGY_THRESHOLD) {
 			for(auto& det : detectors) {
+				auto chan = det.GetTrajectoryRingWedge(eject->theta, eject->phi);
+				if(chan.first != -1 && chan.second != -1) {
 					coordinates = det.GetTrajectoryCoordinates(eject->theta, eject->phi);
-				if(coordinates.GetX() != 0) {
 					eject_xs.push_back(coordinates.GetX());
 					eject_ys.push_back(coordinates.GetY());
 					eject_zs.push_back(coordinates.GetZ());
@@ -135,18 +136,18 @@ void SabreEfficiency::RunDecay(const char* file) {
 
 	TGraph2D* gde = new TGraph2D(eject_xs.size(), &(eject_xs[0]), &(eject_ys[0]), &(eject_zs[0]));
 	gde->SetName("detected_eject_points");
-	gde->SetMarkerStyle(2);
+	gde->SetMarkerStyle(1);
 	gde->SetMarkerColor(2);
 
 	TGraph2D* gr = new TGraph2D(ringxs.size(), &(ringxs[0]), &(ringys[0]), &(ringzs[0]));
 	gr->SetName("ring_detector_edges");
 	gr->SetTitle("SABRE Detector; x(m); y(m); z(m)");
-	gr->SetMarkerStyle(2);
+	gr->SetMarkerStyle(1);
 
 	TGraph2D* gw = new TGraph2D(wedgexs.size(), &(wedgexs[0]), &(wedgeys[0]), &(wedgezs[0]));
 	gw->SetName("wedge_detector_edges");
 	gw->SetTitle("SABRE Detector Wedges; x(m); y(m); z(m)");
-	gw->SetMarkerStyle(2);
+	gw->SetMarkerStyle(1);
 
 	TCanvas* canvas = new TCanvas();
 	canvas->SetName("detectors_and_particles");
@@ -171,8 +172,8 @@ void SabreEfficiency::Run2Step(const char* file) {
 	TFile* input = TFile::Open(file, "UPDATE");
 	TTree* tree = (TTree*) input->Get("DataTree");
 
-	NucData* break1 = new NucData();
+	Mask::NucData* break1 = new Mask::NucData();
-	NucData* break2 = new NucData();
+	Mask::NucData* break2 = new Mask::NucData();
 
 
 	tree->SetBranchAddress("breakup1", &break1);
@@ -237,9 +238,9 @@ void SabreEfficiency::Run3Step(const char* file) {
 	TFile* input = TFile::Open(file, "UPDATE");
 	TTree* tree = (TTree*) input->Get("DataTree");
 
-	NucData* break1 = new NucData();
+	Mask::NucData* break1 = new Mask::NucData();
-	NucData* break3 = new NucData();
+	Mask::NucData* break3 = new Mask::NucData();
-	NucData* break4 = new NucData();
+	Mask::NucData* break4 = new Mask::NucData();
 
 
 	tree->SetBranchAddress("breakup1", &break1);

src/ThreeStepSystem.cpp

@@ -1,6 +1,8 @@
 #include "ThreeStepSystem.h"
 #include "KinematicsExceptions.h"
 
+namespace Mask {
+
 ThreeStepSystem::ThreeStepSystem() :
 	ReactionSystem()
 {
@@ -105,3 +107,5 @@ void ThreeStepSystem::RunSystem() {
 	step3.Calculate();
 
 }
+
+};

src/TwoStepSystem.cpp

@@ -1,6 +1,8 @@
 #include "TwoStepSystem.h"
 #include "KinematicsExceptions.h"
 
+namespace Mask {
+
 TwoStepSystem::TwoStepSystem() :
 	ReactionSystem()
 {
@@ -89,3 +91,5 @@ void TwoStepSystem::RunSystem() {
 
 
 }
+
+};

src/Vec3.cpp

@@ -0,0 +1,46 @@
+/*
+	Class to represent a 3-space vector in both cartesian and spherical coordinates. Can perform vector
+	addition, subtraction, and dot product. 
+
+	--GWM Dec 2020
+*/
+#include "Vec3.h"
+
+namespace Mask {
+
+Vec3::Vec3() {
+	m_data[0] = 0.;
+	m_data[1] = 0.;
+	m_data[2] = 0.;
+}
+
+Vec3::Vec3(double x, double y, double z) {
+	m_data[0] = x;
+	m_data[1] = y;
+	m_data[2] = z;
+}
+
+Vec3::~Vec3() {}
+
+void Vec3::SetVectorCartesian(double x, double y, double z) {
+	m_data[0] = x;
+	m_data[1] = y;
+	m_data[2] = z;
+}
+
+void Vec3::SetVectorSpherical(double r, double theta, double phi) {
+	m_data[0] = r*std::cos(phi)*std::sin(theta);
+	m_data[1] = r*std::sin(phi)*std::sin(theta);
+	m_data[2] = r*std::cos(theta);
+}
+
+double Vec3::Dot(const Vec3& rhs) const {
+	return GetX()*rhs.GetX() + GetY()*rhs.GetY() + GetZ()*rhs.GetZ();
+}
+
+//Unimplemented
+Vec3 Vec3::Cross(const Vec3& rhs) const {
+	return Vec3(0.,0.,0.);
+}
+
+};

src/Vec4.cpp

@@ -1,14 +1,23 @@
-#include "G4Vec.h"
+/*
-//NOTE: uses (-,-,-,+) metric (same as ROOT convention)
+	Class which represents a 4-momentum vector. Can perform vector addition, subtraction, dot product
+	and generate a boost vector to its rest frame as well as apply a boost to itself.
 
-G4Vec::G4Vec() {
+	--GWM Dec 2020.
+	NOTE: uses (-,-,-,+) metric (same as ROOT convention)
+*/
+#include "Vec4.h"
+
+
+namespace Mask {
+
+Vec4::Vec4() {
 	for(auto& val: m_data)
 		val = 0.0;
 	for(auto& val: m_boost)
 		val = 0.0;
 }
 
-G4Vec::G4Vec(double px, double py, double pz, double E) {
+Vec4::Vec4(double px, double py, double pz, double E) {
 	m_data[0] = px;
 	m_data[1] = py;
 	m_data[2] = pz;
@@ -16,9 +25,9 @@ G4Vec::G4Vec(double px, double py, double pz, double E) {
 	CalcBoostToCM();
 }
 
-G4Vec::~G4Vec() {}
+Vec4::~Vec4() {}
 
-void G4Vec::SetVectorCartesian(double px, double py, double pz, double E) {
+void Vec4::SetVectorCartesian(double px, double py, double pz, double E) {
 	m_data[0] = px;
 	m_data[1] = py;
 	m_data[2] = pz;
@@ -27,7 +36,7 @@ void G4Vec::SetVectorCartesian(double px, double py, double pz, double E) {
 	CalcBoostToCM();
 }
 
-void G4Vec::SetVectorSpherical(double theta, double phi, double p, double E) {
+void Vec4::SetVectorSpherical(double theta, double phi, double p, double E) {
 	m_data[0] = p*cos(phi)*sin(theta);
 	m_data[1] = p*sin(phi)*sin(theta);
 	m_data[2] = p*cos(theta);
@@ -35,13 +44,13 @@ void G4Vec::SetVectorSpherical(double theta, double phi, double p, double E) {
 	CalcBoostToCM();
 }
 
-void G4Vec::CalcBoostToCM() {
+void Vec4::CalcBoostToCM() {
 	m_boost[0] = m_data[0]/m_data[3];
 	m_boost[1] = m_data[1]/m_data[3];
 	m_boost[2] = m_data[2]/m_data[3];
 }
 
-void G4Vec::ApplyBoost(const double* beta) {
+void Vec4::ApplyBoost(const double* beta) {
 	double beta2 = beta[0]*beta[0] + beta[1]*beta[1] + beta[2]*beta[2];
 	double gamma  = 1.0/std::sqrt(1.0 - beta2);
 	double bdotp = beta[0]*m_data[0] + beta[1]*m_data[1] + beta[2]*m_data[2];
@@ -53,10 +62,12 @@ void G4Vec::ApplyBoost(const double* beta) {
 		      		   gamma*(GetE() + bdotp));
 }
 
-double G4Vec::Dot(const G4Vec& rhs) const {
+double Vec4::Dot(const Vec4& rhs) const {
 	return GetE()*rhs.GetE() - GetPx()*rhs.GetPx() - GetPy()*rhs.GetPy() - GetPz()*rhs.GetPz();
 }
 
-G4Vec G4Vec::Cross(const G4Vec& rhs) const {
+Vec4 Vec4::Cross(const Vec4& rhs) const {
-	return G4Vec();
+	return Vec4();
 }
+
+};

src/main.cpp

@@ -10,7 +10,7 @@ int main(int argc, char** argv) {
 		return 1;
 	}
 	
-	Kinematics calculator;
+	Mask::Kinematics calculator;
 	try {
 		if(!calculator.LoadConfig(argv[1])) {
 			return 1;
@@ -45,7 +45,7 @@ int main(int argc, char** argv) {
  	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
  	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
 
- 	double theta, phi, expected_flat_t, expected_flat_p;
+ 	double theta, phi, expected_flat_p;
  	for(int h=0; h<5; h++) {
  		for(int j=0; j<16; j++) {
  			for(int k=0; k<4; k ++) {
@@ -53,7 +53,11 @@ int main(int argc, char** argv) {
  				phi = detectors[h].GetRingTiltCoords(j, k).GetPhi();
  				expected_flat_p = detectors[h].GetRingFlatCoords(j, k).GetPhi();  
  				for(int i=0; i<5; i++) {
- 					if(detectors[i].GetTrajectoryCoordinates(theta, phi).GetX() != 0) {
+ 					auto channels = detectors[i].GetTrajectoryRingWedge(theta, phi);
+ 					if(channels.first != -1) {
+ 						std::cout<<"Detected in detector"<<i<<" ring: "<<channels.first<<" wedge: "<<channels.second<<" Expected -- detector: "<<h<<" ring: "<<j;
+ 						if(k == 0 || k == 1) std::cout<<" wedge: 0"<<std::endl;
+ 						else std::cout<<" wedge: 7"<<std::endl;
  						break;
  					} else if(i == 4) {
  						std::cout<<" Not found! detector: "<<h<<" ring: "<<j<<" corner: "<<k<<" theta: "<<theta/DEG2RAD<<" phi: "<<phi/DEG2RAD<<" flat_p: "<<expected_flat_p/DEG2RAD<<std::endl;