Major overhaul. Fixed bugs in decay angular distribution, now classed. Added ANASEN geometry. Minor bugfixes throughout. Added center of mass switch for main reaction (better randomization where ejectile is not fixed in theta).

etc/12C3Hed_pdecay_L2_fit.txt

@@ -0,0 +1,4 @@
+AngularMomentum: 2
+A0: 0.2121
+A2: 0.2910
+A4: 0.03636

etc/12C3Hed_pdecay_L2_theory.txt

@@ -0,0 +1,4 @@
+AngularMomentum: 2
+A0: 0.5
+A2: 0.544
+A4: 0.362

etc/isotropic_dist.txt

@@ -0,0 +1,2 @@
+AngularMomentum: 0
+A0: 0.5

include/AnasenEfficiency.h

@@ -0,0 +1,60 @@
+#ifndef ANASEN_EFFICIENCY_H
+#define ANASEN_EFFICIENCY_H
+
+#include <string>
+#include <THashTable.h>
+
+#include "StripDetector.h"
+#include "QQQDetector.h"
+
+class AnasenEfficiency {
+public:
+	AnasenEfficiency();
+	~AnasenEfficiency();
+	void CalculateEfficiency(const std::string& filename);
+	inline void SetReactionType(int type) { m_rxn_type = type; };
+	void DrawDetectorSystem(const std::string& filename);
+	double RunConsistencyCheck();
+
+private:
+	void MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val);
+    void MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy);
+	void RunDecay(const std::string& filename);
+	void RunTwoStep(const std::string& filename);
+	void RunThreeStep(const std::string& filename);
+
+	bool IsRing1(double theta, double phi);
+	bool IsRing2(double theta, double phi);
+	bool IsQQQ(double theta, double phi);
+
+	inline bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < epsilon ? true : false; };
+
+	int m_rxn_type;
+	std::vector<StripDetector> m_Ring1, m_Ring2;
+	std::vector<QQQDetector> m_forwardQQQs;
+	std::vector<QQQDetector> m_backwardQQQs;
+
+	/**** ANASEN geometry constants *****/
+	const int n_sx3_per_ring = 12;
+	const int n_qqq = 4;
+	const double sx3_length = 0.075;
+	const double sx3_width = 0.04;
+	const double barrel_gap = 0.013 + 0.049; //0.049 is base gap due to frames
+	const double ring1_z = sx3_length/2.0 + barrel_gap/2.0;
+	//const double ring2_z = -0.124 + sx3_length/2.0 + 0.0245 - barrel_gap/2.0;
+	const double qqq_nom_z = 0.025 + sx3_length + 0.0245 + barrel_gap/2.0;
+	const double qqq_rinner = 0.0501;
+	const double qqq_router = 0.0990;
+	const double qqq_deltaphi = 1.52119;
+	const double qqq_z[4] = {qqq_nom_z, qqq_nom_z - 0.00828, qqq_nom_z, qqq_nom_z};
+	const double qqq_phi[4] = {5.49779, 0.785398, 2.35619, 3.92699};
+	const double ring_rho[12] = {0.0890601, 0.0889871, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354};
+	const double ring_phi[12] = {0.785795, 0.262014, 6.02132, 5.49779, 4.97426, 4.45052, 3.92699, 3.40346, 2.87972, 2.35619, 1.83266, 1.30893};
+	/*************************/
+
+	static constexpr double threshold = 0.2; //MeV
+	static constexpr double epsilon = 1.0e-6;
+	static constexpr double deg2rad = M_PI/180.0;
+};
+
+#endif

include/AngularDistribution.h

@@ -0,0 +1,37 @@
+#ifndef ANGULARDISTRIBUTION_H
+#define ANGULARDISTRIBUTION_H
+
+#include <TRandom3.h>
+#include <string>
+#include <vector>
+
+class AngularDistribution {
+public:
+	AngularDistribution();
+	AngularDistribution(const std::string& file);
+	~AngularDistribution();
+	void ReadDistributionFile(const std::string& file);
+	void AttachRandomNumberGenerator(TRandom3* random) { generator = random; };
+	double GetRandomCosTheta();
+	int GetL() { return L; };
+	double GetBranchingRatio() { return branchingRatio; };
+
+private:
+	bool IsIsotropic();
+	bool IsGeneratorSet() {
+		if(generator) {
+			return true;
+		} else {
+			return false;
+		}
+	}
+
+	TRandom3* generator; //NOT OWNED BY ANGULAR DISTRIBUTION
+
+	double branchingRatio;
+	int L;
+	std::vector<double> constants;
+	bool isoFlag;
+};
+
+#endif

include/KinematicsExceptions.h

@@ -2,6 +2,7 @@
 #define KINEMATICSEXCEPTIONS_H
 
 #include <exception>
+#include <string>
 #include <stdexcept>
 /*
 	ELossException
@@ -12,8 +13,16 @@
 	for locations where this could be thrown
 */
 struct ELossException : public std::exception {
+	ELossException(const std::string& error) {
+		m_error = error;
+	}
+
+	std::string m_error;
+
 	const char* what() const noexcept {
-		return "Failure to calculate particle energy loss. See KinematicsExceptions.h for documentation.";
+		std::string err_str = "Failure to calculate particle energy loss. Reason: ";
+		err_str += m_error + " See KinematicsExceptions.h for documentation.";
+		return err_str.c_str();
 	};
 };
 
@@ -61,6 +70,17 @@ struct QValueException : public std::exception {
 	};
 };
 
+/*
+	EnergyThresholdException
+	This is an exception thrown when the Reaction attempts to calculate a reaction which does not have enough incoming (read: beam kinetic energy) energy
+	to occur even for the 0 degree case. The reaction is kinematically forbidden.
+*/
+struct EnergyThresholdException : public std::exception {
+	const char* what() const noexcept {
+		return "Reaction does not have enough energy to proceed. See KinematicsExceptions.h for documentation.";
+	};
+};
+
 
 
 #endif

include/LegendrePoly.h

@@ -2,5 +2,11 @@
 #define LEGENDREPOLY_H
 
 double P_l(int l, double x);
+double P_l_ROOT(double* x, double* pars);
+double Normed_P_l_sq(int l, double x);
+
+double P_0(double x);
+double P_1(double x);
+double P_2(double x);
 
 #endif

include/MaskFile.h

@@ -0,0 +1,54 @@
+#ifndef MASKFILE_H
+#define MASKFILE_H
+
+#include <string>
+#include <fstream>
+#include <vector>
+
+#include "Nucleus.h"
+
+namespace Mask {
+
+class MaskFile {
+public:
+	enum class FileType {
+		read,
+		write,
+		append,
+		none
+	};
+
+	MaskFile();
+	MaskFile(const std::string& name, MaskFile::FileType type);
+	bool Open(const std::string& name, MaskFile::FileType type);
+	inline bool IsOpen() { return file.is_open(); }
+	void Close();
+	
+	void WriteHeader(std::vector<Nucleus>& data);
+	void WriteData(std::vector<Nucleus>& data);
+	void ReadHeader();
+	std::vector<Nucleus> ReadData();
+
+	
+
+private:
+	void FlushBuffer();
+	void FillBuffer();
+
+	FileType file_type;
+	std::string filename;
+	unsigned int buffer_position;
+	int count;
+
+	std::fstream file;
+
+	std::vector<std::vector<Nucleus>> data_buffer;
+
+	static constexpr int buffersize = 1000;
+	static constexpr int width = 0;
+	static constexpr int precision = 3;
+};
+
+};
+
+#endif

include/QQQDetector.h

@@ -0,0 +1,48 @@
+/*
+	QQQDetector.h
+	Class implementing geometry for QQQDetector where the detector is perpendicular to the beam axis.
+	Detector is first generated centered on the x-axis (phi=0)
+
+	Coordinate convention : +z is downstream, -z is upstream. +y is vertically up in the lab.
+*/
+#ifndef QQQDETECTOR_H
+#define QQQDETECTOR_H
+
+#include <TRandom3.h>
+#include <cmath>
+
+#include "Vec3.h"
+#include "Rotation.h"
+
+class QQQDetector {
+public:
+	QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x=0, double y=0);
+	~QQQDetector();
+	inline Mask::Vec3 GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; };
+	inline Mask::Vec3 GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; };
+	Mask::Vec3 GetTrajectoryCoordinates(double theta, double phi);
+	std::pair<int, int> GetTrajectoryRingWedge(double theta, double phi);
+	Mask::Vec3 GetHitCoordinates(int ringch, int wedgech);
+
+	inline int GetNumberOfRings() { return nrings; };
+	inline int GetNumberOfWedges() { return nwedges; };
+
+private:
+
+	inline bool CheckChannel(int ch) { return (ch >=0 && ch < nrings); };
+	inline bool CheckCorner(int corner) { return (corner >=0 && corner < 4); };
+
+	void CalculateCorners();
+	Mask::Vec3 TransformCoordinates(Mask::Vec3& vector) { return m_ZRot*vector + m_translation; };
+
+	double m_Rinner, m_Router, m_deltaR, m_deltaPhi, m_deltaPhi_per_wedge, m_phiCentral;
+	std::vector<std::vector<Mask::Vec3>> m_ringCoords, m_wedgeCoords;
+	Mask::Vec3 m_translation;
+	Mask::ZRotation m_ZRot;
+
+	static constexpr int nrings = 16;
+	static constexpr int nwedges = 16;
+	static constexpr double deg2rad  = M_PI/180.0;
+};
+
+#endif

include/Reaction.h

@@ -23,6 +23,7 @@ public:
 	void SetNuclei(int zt, int at, int zp, int ap, int ze, int ae);
 	void SetNuclei(const Nucleus* nucs);
 	void SetBeamKE(double bke);
+	void SetEjectileThetaType(int type);
 
 	/*Setters and getters*/
 	inline void SetLayeredTarget(LayeredTarget* targ) { target = targ; };
@@ -46,10 +47,16 @@ public:
 	inline const Nucleus& GetTarget() const { return reactants[0]; };
 	inline const Nucleus& GetEjectile() const { return reactants[2]; };
 	inline const Nucleus& GetResidual() const { return reactants[3]; };
+	inline void ResetTarget() { reactants[0].SetVectorCartesian(0,0,0,0); };
+	inline void ResetProjectile() { reactants[1].SetVectorCartesian(0,0,0,0); };
+	inline void ResetEjectile() { reactants[2].SetVectorCartesian(0,0,0,0); };
+	inline void ResetResidual() { reactants[3].SetVectorCartesian(0,0,0,0); };
 	inline int GetRxnLayer() { return rxnLayer; };
 
 private:
 	void CalculateReaction(); //target + project -> eject + resid
+	void CalculateReactionThetaLab();
+	void CalculateReactionThetaCM();
 	void CalculateDecay(); //target -> light_decay (eject) + heavy_decay(resid)
 
 	Nucleus reactants[4]; //0=target, 1=projectile, 2=ejectile, 3=residual
@@ -57,10 +64,14 @@ private:
 
 	double m_bke, m_theta, m_phi, m_ex;
 
-	int  rxnLayer;
+	int rxnLayer;
+	int m_eject_theta_type; 
 
 	bool decayFlag, nuc_initFlag, resid_elossFlag;
 
+	static constexpr int lab = 0;
+	static constexpr int center_of_mass = 1;
+
 };
 
 };

include/ReactionSystem.h

@@ -10,6 +10,7 @@
 #define REACTIONSYSTEM_H
 
 #include "Reaction.h"
+#include "AngularDistribution.h"
 #include <vector>
 #include <TRandom3.h>
 
@@ -25,15 +26,14 @@ public:
 	void AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness);
 
 	/*Set sampling parameters*/
-	inline void SetRandomGenerator(TRandom3* gen) { generator = gen; gen_set_flag = true; };
+	void SetRandomGenerator(TRandom3* gen);
 	inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); };
+	inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
 	inline void SetTheta1Range(double min, double max) { m_theta1Range = std::make_pair(min*deg2rad, max*deg2rad); };
+	inline void SetPhi1Range(double min, double max) { m_phi1Range = std::make_pair(min*deg2rad, max*deg2rad); };
 	inline void SetExcitationDistro(double mean, double sigma) { m_exDist = std::make_pair(mean, sigma); };
-	inline void SetDecay1AngularMomentum(double l) { L1 = l; };
+	inline void SetDecay1Distribution(const std::string& file) { decay1dist.ReadDistributionFile(file); };
-	inline void SetDecay2AngularMomentum(double l) { L2 = l; };
+	inline void SetDecay2Distribution (const std::string& file) { decay2dist.ReadDistributionFile(file); };
-
-	/*Sampling over angular distribution*/
-	double GetDecayTheta(int L);
 
 	virtual void RunSystem();
 
@@ -42,6 +42,10 @@ public:
 	inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
 	inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
 	inline const char* GetSystemEquation() const { return m_sys_equation.c_str(); };
+	inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
+	inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
+	inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
+	inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
 
 protected:
 	virtual void LinkTarget();
@@ -51,12 +55,13 @@ protected:
 	LayeredTarget target;
 
 	//Sampling information
-	std::pair<double, double> m_beamDist, m_theta1Range, m_exDist;
+	std::pair<double, double> m_beamDist, m_theta1Range, m_phi1Range, m_exDist;
 	TRandom3* generator; //not owned by ReactionSystem
 
+	AngularDistribution decay1dist, decay2dist;
+
 	bool target_set_flag, gen_set_flag;
 	int rxnLayer;
-	int L1, L2;
 	std::string m_sys_equation;
 	static constexpr double deg2rad = M_PI/180.0;
 };

include/SabreEfficiency.h

@@ -4,6 +4,7 @@
 #include "SabreDetector.h"
 #include "Target.h"
 #include "DeadChannelMap.h"
+#include "Kinematics.h"
 #include <THashTable.h>
 
 class SabreEfficiency {
@@ -12,20 +13,23 @@ public:
 	~SabreEfficiency();
 	inline void SetReactionType(int t) { m_rxn_type = t; };
     void SetDeadChannelMap(std::string& filename) { dmap.LoadMapfile(filename); };
-	void CalculateEfficiency(const char* file);
+	void CalculateEfficiency(const std::string& file);
+    void DrawDetectorSystem(const std::string& filename);
+    double RunConsistencyCheck();
 
 private:
-    void MyFill(THashTable* table, const char* name, const char* title, int bins, float min, float max, double val);
+    void MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val);
-    void MyFill(THashTable* table, const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy);
+    void MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy);
-	void Run2Step(const char*);
+	std::pair<bool,double> IsSabre(Mask::NucData* nucleus);
-	void Run3Step(const char*);
+    void Run2Step(const std::string& filename);
-	void RunDecay(const char*);
+	void Run3Step(const std::string& filename);
+	void RunDecay(const std::string& filename);
 
 	int m_rxn_type;
 	std::vector<SabreDetector> detectors;
-    std::vector<double> ringxs, ringys, ringzs;
+    
-    std::vector<double> wedgexs, wedgeys, wedgezs;
 	Target deadlayer;
+    Target sabre_eloss;
     DeadChannelMap dmap;
 
 
@@ -42,6 +46,7 @@ private:
     const double PHI4 = 90.0;
     const double DEG2RAD = M_PI/180.0;
     static constexpr double DEADLAYER_THIN = 50 * 1e-7 * 2.3296 * 1e6; // ug/cm^2 (50 nm thick * density)
+    static constexpr double SABRE_THICKNESS = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density)
 
     const double ENERGY_THRESHOLD = 0.2; //in MeV
 

include/StripDetector.h

@@ -0,0 +1,60 @@
+#ifndef __STRIP_DETECTOR_H
+#define __STRIP_DETECTOR_H
+
+// +z is along beam axis
+// +y is vertically "upward" in the lab frame
+
+//angles must be in radians, but distances can be whatever
+//PROVIDED all input distances are the same
+
+//Front strips from lowest y to highest y
+
+//Back strips from lowest z to highest z
+
+#include <TRandom3.h>
+#include <cmath>
+#include <vector>
+
+#include "Vec3.h"
+#include "Rotation.h"
+
+class StripDetector {
+public:
+  
+	StripDetector(int ns, double len, double wid, double cphi, double cz, double crho);
+	~StripDetector();
+	inline Mask::Vec3 GetFrontStripCoordinates(int stripch, int corner) { return front_strip_coords[stripch][corner]; };
+	inline Mask::Vec3 GetBackStripCoordinates(int stripch, int corner) { return back_strip_coords[stripch][corner]; };
+	inline Mask::Vec3 GetRotatedFrontStripCoordinates(int stripch, int corner) { return rotated_front_strip_coords[stripch][corner]; };
+	inline Mask::Vec3 GetRotatedBackStripCoordinates(int stripch, int corner) { return rotated_back_strip_coords[stripch][corner]; };
+	inline void SetRandomNumberGenerator(TRandom3* random) { m_random = random; };
+	Mask::Vec3 GetHitCoordinates(int front_stripch, double front_strip_ratio);
+	std::pair<int,double> GetChannelRatio(double theta, double phi);
+
+private:
+	void CalculateCorners();
+
+	int num_strips;
+	static constexpr int num_corners = 4;
+
+	double length; //common to all strips, hence total
+	double total_width;
+	double front_strip_width; //assuming equal widths
+	double back_strip_length; //assuming equal widths
+
+	double center_phi; //assuming det centered above x-axis (corresponds to zero phi)
+	double center_z;
+	double center_rho; //perpendicular radius from axis
+
+	std::vector<std::vector<Mask::Vec3>> front_strip_coords, back_strip_coords;
+	std::vector<std::vector<Mask::Vec3>> rotated_front_strip_coords, rotated_back_strip_coords;
+
+	Mask::ZRotation zRot;
+
+	TRandom3* m_random; //Not owned by StripDetector!
+
+	inline bool ValidChannel(int f) { return ((f >= 0 && f < num_strips) ? true : false); };
+	inline bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); };
+};
+
+#endif

input.txt

@@ -1,28 +1,31 @@
 ----------Data Information----------
-OutputFile: /data1/gwm17/10B3He/Feb2021/12C3Hed_13N3546keV_ps_fixedThetaBins_test.root
+OutputFile: /data1/gwm17/TRIUMF_7Bed/simulation/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root
 SaveTree: yes
 SavePlots: yes
 ----------Reaction Information----------
 ReactionType: 2
 Z	A (order is target, projectile, ejectile, break1, break3)
-6	12
-2	3
 1	2
+4	7
 1	1
+2	4
 ----------Target Information----------
 Name: test_targ
 Layers: 1
 ~Layer1
-Thickness(ug/cm^2): 41
+Thickness(ug/cm^2): 50
 Z	A	Stoich
+1	2	2
 6	12	1
 0
 ~
 ----------Sampling Information----------
 NumberOfSamples: 100000
-BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001
+BeamMeanEnergy(MeV): 0.6 BeamEnergySigma(MeV): 0.001
-EjectileThetaMin(deg): 3.0 EjectileThetaMax(deg): 3.0
+EjectileThetaType(0=Lab,1=CM): 1
-ResidualExMean(MeV): 3.546 ResidualExSigma(MeV): 0.02
+EjectileThetaMin(deg): 0.0 EjectileThetaMax(deg): 180.0
-Decay1_AngularMomentum: 0
+EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 360.0
-Decay2_AngularMomentum: 0
+ResidualExMean(MeV): 0.0 ResidualExSigma(MeV): 0.0
+Decay1_DistributionFile: ./etc/isotropic_dist.txt
+Decay2_DistributionFile: ./etc/isotropic_dist.txt
 --------------------------------------

src/AnasenEfficiency.cpp

@@ -0,0 +1,562 @@
+#include "AnasenEfficiency.h"
+#include "Kinematics.h"
+#include <TH2.h>
+#include <TH1.h>
+#include <TGraph2D.h>
+#include <TGraph.h>
+#include <TCanvas.h>
+#include <TParameter.h>
+
+AnasenEfficiency::AnasenEfficiency() {
+	for(int i=0; i<n_sx3_per_ring; i++) {
+		m_Ring1.emplace_back(4, sx3_length, sx3_width, ring_phi[i], ring1_z, ring_rho[i]);
+		m_Ring2.emplace_back(4, sx3_length, sx3_width, ring_phi[i], -1.0*ring1_z, ring_rho[i]);
+	}
+	for(int i=0; i<n_qqq; i++) {
+		m_forwardQQQs.emplace_back(qqq_rinner, qqq_router, qqq_deltaphi, qqq_phi[i], qqq_z[i]);
+		m_backwardQQQs.emplace_back(qqq_rinner, qqq_router, qqq_deltaphi, qqq_phi[i], (-1.0)*qqq_z[i]);
+	}
+}
+
+AnasenEfficiency::~AnasenEfficiency() {}
+
+void AnasenEfficiency::MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val) {
+	TH1F* h = (TH1F*) table->FindObject(name.c_str());
+	if(h) {
+		h->Fill(val);
+	} else {
+		h = new TH1F(name.c_str(), title.c_str(), bins, min, max);
+		h->Fill(val);
+		table->Add(h);
+	}
+}
+
+void AnasenEfficiency::MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy) {
+	TH2F* h = (TH2F*) table->FindObject(name.c_str());
+	if(h) {
+		h->Fill(valx, valy);
+	} else {
+		h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy);
+		h->Fill(valx, valy);
+		table->Add(h);
+	}
+}
+
+void AnasenEfficiency::DrawDetectorSystem(const std::string& filename) {
+	TFile* file = TFile::Open(filename.c_str(), "RECREATE");
+
+	std::vector<double> x, y, z;
+	std::vector<double> cx, cy, cz;
+	Mask::Vec3 coords;
+	for(int i=0; i<n_sx3_per_ring; i++) {
+		for(int j=0; j<4; j++) {
+			for(int k=0; k<4; k++) {
+				coords = m_Ring1[i].GetRotatedFrontStripCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+				coords = m_Ring1[i].GetRotatedBackStripCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+			}
+			coords = m_Ring1[i].GetHitCoordinates(j, 0);
+			cx.push_back(coords.GetX());
+			cy.push_back(coords.GetY());
+			cz.push_back(coords.GetZ());
+		}
+	}
+	for(int i=0; i<n_sx3_per_ring; i++) {
+		for(int j=0; j<4; j++) {
+			for(int k=0; k<4; k++) {
+				coords = m_Ring2[i].GetRotatedFrontStripCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+				coords = m_Ring2[i].GetRotatedBackStripCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+			}
+			coords = m_Ring2[i].GetHitCoordinates(j, 0);
+			cx.push_back(coords.GetX());
+			cy.push_back(coords.GetY());
+			cz.push_back(coords.GetZ());
+		}
+	}
+	for(int i=0; i<n_qqq; i++) {
+		for(int j=0; j<16; j++) {
+			for(int k=0; k<4; k++) {
+				coords = m_forwardQQQs[i].GetRingCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+				coords = m_forwardQQQs[i].GetWedgeCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+			}
+			for(int k=0; k<16; k++) {
+				coords = m_forwardQQQs[i].GetHitCoordinates(j, k);
+				cx.push_back(coords.GetX());
+				cy.push_back(coords.GetY());
+				cz.push_back(coords.GetZ());
+			}
+		}
+	}
+	for(int i=0; i<n_qqq; i++) {
+		for(int j=0; j<16; j++) {
+			for(int k=0; k<4; k++) {
+				coords = m_backwardQQQs[i].GetRingCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+				coords = m_backwardQQQs[i].GetWedgeCoordinates(j, k);
+				x.push_back(coords.GetX());
+				y.push_back(coords.GetY());
+				z.push_back(coords.GetZ());
+			}
+			for(int k=0; k<16; k++) {
+				coords = m_backwardQQQs[i].GetHitCoordinates(j, k);
+				cx.push_back(coords.GetX());
+				cy.push_back(coords.GetY());
+				cz.push_back(coords.GetZ());
+			}
+		}
+	}
+	TGraph2D* graph = new TGraph2D(x.size(), &(x[0]), &(y[0]), &(z[0]));
+	graph->SetName("CornerGraph");
+	graph->SetMarkerStyle(2);
+	graph->SetLineColor(1);
+
+	TGraph2D* graph2 = new TGraph2D(cx.size(), &(cx[0]), &(cy[0]), &(cz[0]));
+	graph2->SetName("CenterGraph");
+	graph2->SetMarkerStyle(2);
+	graph2->SetMarkerColor(4);
+
+	TCanvas* canvas = new TCanvas();
+	canvas->SetName("ANASEN Detector");
+	graph->Draw("A|P");
+	graph2->Draw("same|P");
+	canvas->Write();
+	graph->Write();
+	graph2->Write();
+	file->Close();
+}
+
+double AnasenEfficiency::RunConsistencyCheck() {
+	std::vector<Mask::Vec3> r1_points;
+	std::vector<Mask::Vec3> r2_points;
+	std::vector<Mask::Vec3> fqqq_points;
+	std::vector<Mask::Vec3> bqqq_points;
+	for(int i=0; i<n_sx3_per_ring; i++) {
+		for(int j=0; j<4; j++) {
+			r1_points.push_back(m_Ring1[i].GetHitCoordinates(j, 0));
+		}
+	}
+	for(int i=0; i<n_sx3_per_ring; i++) {
+		for(int j=0; j<4; j++) {
+			r2_points.push_back(m_Ring2[i].GetHitCoordinates(j, 0));
+		}
+	}
+	for(int i=0; i<n_qqq; i++) {
+		for(int j=0; j<16; j++) {
+			for(int k=0; k<16; k++) {
+				fqqq_points.push_back(m_forwardQQQs[i].GetHitCoordinates(j, k));
+			}
+		}
+	}
+	for(int i=0; i<n_qqq; i++) {
+		for(int j=0; j<16; j++) {
+			for(int k=0; k<16; k++) {
+				bqqq_points.push_back(m_backwardQQQs[i].GetHitCoordinates(j, k));
+			}
+		}
+	}
+
+	int npoints = r1_points.size() + r2_points.size() + fqqq_points.size() + bqqq_points.size();
+	int count = 0;
+	Mask::Vec3 coords;
+	for(auto& point : r1_points) {
+		for(auto& sx3 : m_Ring1) {
+			auto result = sx3.GetChannelRatio(point.GetTheta(), point.GetPhi());
+			coords = sx3.GetHitCoordinates(result.first, result.second);
+			if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
+				count++;
+				break;
+			}
+		}
+	}
+	for(auto& point : r2_points) {
+		for(auto& sx3 : m_Ring2) {
+			auto result = sx3.GetChannelRatio(point.GetTheta(), point.GetPhi());
+			coords = sx3.GetHitCoordinates(result.first, result.second);
+			if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
+				count++;
+				break;
+			}
+		}
+	}
+	for(auto& point : fqqq_points) {
+		for(auto& qqq : m_forwardQQQs) {
+			auto result = qqq.GetTrajectoryRingWedge(point.GetTheta(), point.GetPhi());
+			coords = qqq.GetHitCoordinates(result.first, result.second);
+			if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
+				count++;
+				break;
+			}
+		}
+	}
+	for(auto& point : bqqq_points) {
+		for(auto& qqq : m_backwardQQQs) {
+			auto result = qqq.GetTrajectoryRingWedge(point.GetTheta(), point.GetPhi());
+			coords = qqq.GetHitCoordinates(result.first, result.second);
+			if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
+				count++;
+				break;
+			}
+		}
+	}
+
+	double ratio = ((double)count)/((double)npoints);
+
+	return ratio;
+
+}
+
+bool AnasenEfficiency::IsRing1(double theta, double phi) {
+	for(auto& sx3 : m_Ring1) {
+		auto result = sx3.GetChannelRatio(theta, phi);
+		if(result.first != -1) {
+			return true;
+		}
+	}
+
+	return false;
+}
+
+bool AnasenEfficiency::IsRing2(double theta, double phi) {
+	for(auto& sx3 : m_Ring2) {
+		auto result = sx3.GetChannelRatio(theta, phi);
+		if(result.first != -1) {
+			return true;
+		}
+	}
+
+	return false;
+}
+
+bool AnasenEfficiency::IsQQQ(double theta, double phi) {
+	for(auto& qqq : m_forwardQQQs) {
+		auto result = qqq.GetTrajectoryRingWedge(theta, phi);
+		if(result.first != -1) {
+			return true;
+		}
+	}
+
+	
+	for(auto& qqq : m_backwardQQQs) {
+		auto result = qqq.GetTrajectoryRingWedge(theta, phi);
+		if(result.first != -1) {
+			return true;
+		}
+	}
+
+
+
+	return false;
+}
+
+void AnasenEfficiency::CalculateEfficiency(const std::string& file) {
+	std::cout<<"----------ANASEN Efficiency Calculation----------"<<std::endl;
+	std::cout<<"Loading in output from kinematics simulation: "<<file<<std::endl;
+	std::cout<<"Running efficiency calculation..."<<std::endl;
+
+	switch(m_rxn_type) {
+		case Mask::Kinematics::ONESTEP_DECAY:
+		{
+			RunDecay(file);
+			break;
+		}
+		case Mask::Kinematics::TWOSTEP:
+		{
+			RunTwoStep(file);
+			break;
+		}
+		case Mask::Kinematics::THREESTEP:
+		{
+			RunThreeStep(file);
+			break;
+		}
+	}
+	std::cout<<std::endl;
+	std::cout<<"Complete."<<std::endl;
+	std::cout<<"---------------------------------------------"<<std::endl;
+}
+
+void AnasenEfficiency::RunDecay(const std::string& filename) {
+	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
+	TTree* tree = (TTree*) input->Get("DataTree");
+
+	THashTable* table = new THashTable();
+	Mask::NucData* eject = new Mask::NucData();
+	Mask::NucData* resid = new Mask::NucData();
+
+	tree->SetBranchAddress("ejectile", &eject);
+	tree->SetBranchAddress("residual", &resid);
+
+	double nevents = tree->GetEntries();
+
+	//Progress tracking
+	int percent5 = nevents*0.05;
+	int count = 0;
+	int npercent = 0;
+
+	int detected_eject=0;
+	int detected_resid=0;
+	for(int i=0; i<tree->GetEntries(); i++) {
+		if(++count == percent5) {//Show progress every 5%
+			npercent++;
+			count = 0;
+			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
+		}
+
+		tree->GetEntry(i);
+
+		if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) {
+			detected_eject++;
+			MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE);
+		}
+
+		if(resid->KE > threshold && (IsRing1(resid->theta, resid->phi) || IsRing2(resid->theta, resid->phi) || IsQQQ(resid->theta, resid->phi))) {
+			detected_resid++;
+			MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/deg2rad,300,0.0,30.0,resid->KE);
+		}
+
+	}
+
+	double ejecteff = ((double) detected_eject)/nevents;
+	double resideff = ((double) detected_resid)/nevents;
+	TParameter<double> eject_eff("Light Breakup Efficiency", ejecteff);
+	TParameter<double> resid_eff("Heavy Breakup Efficiency", resideff);
+
+	input->cd();
+	eject_eff.Write();
+	resid_eff.Write();
+	table->Write();
+	input->Close();
+}
+
+void AnasenEfficiency::RunTwoStep(const std::string& filename) {
+
+	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
+	TTree* tree = (TTree*) input->Get("DataTree");
+
+	Mask::NucData* eject = new Mask::NucData();
+	Mask::NucData* break1 = new Mask::NucData();
+	Mask::NucData* break2 = new Mask::NucData();
+
+	THashTable* table = new THashTable();
+
+	tree->SetBranchAddress("ejectile", &eject);
+	tree->SetBranchAddress("breakup1", &break1);
+	tree->SetBranchAddress("breakup2", &break2);
+
+	double nevents = tree->GetEntries();
+
+	//Progress tracking
+	int percent5 = nevents*0.05;
+	int count = 0;
+	int npercent = 0;
+
+	double costheta_cm =0;
+
+	bool break1_flag, break2_flag, eject_flag;
+	int b1_count=0, b2_count=0, eject_count=0;
+	int b1b2_count=0, b1e_count=0, b2e_count=0, b1b2e_count=0;
+	for(int i=0; i<tree->GetEntries(); i++) {
+		if(++count == percent5) {//Show progress every 5%
+			npercent++;
+			count = 0;
+			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
+		}
+
+		break1_flag = false;
+		break2_flag = false;
+		eject_flag = false;
+
+		tree->GetEntry(i);
+
+		if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) {
+			eject_count++;
+			eject_flag = true;
+			MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE);
+		}
+
+		if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) {
+			b1_count++;
+			break1_flag = true;
+			costheta_cm = std::cos(break1->theta_cm);
+			MyFill(table,"detected_break1_cm_theta","detected breakup1 CM theta;cos(#theta_{CM})",20,-1.0,1.0,costheta_cm);
+			MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE);
+		}
+
+		if(break2->KE >= threshold && (IsRing1(break2->theta, break2->phi) || IsRing2(break2->theta, break2->phi) || IsQQQ(break2->theta, break2->phi))) {
+			b2_count++;
+			break2_flag = true;
+			MyFill(table, "detected_break2_theta_ke","detected breakup2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/deg2rad,300,0.0,30.0,break2->KE);
+		}
+
+		if(break1_flag && break2_flag && eject_flag) {
+			b1b2e_count++;
+			b1b2_count++;
+			b1e_count++;
+			b2e_count++;
+			MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad);
+			MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE);
+			MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad);
+		} else if(break1_flag && eject_flag) {
+			b1e_count++;
+		} else if(break2_flag && eject_flag) {
+			b2e_count++;
+		} else if(break1_flag && break2_flag) {
+			b1b2_count++;
+			MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad);
+			MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE);
+			MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad);
+		}
+	}
+
+	double eeff = ((double) eject_count)/nevents;
+	double b1eff = ((double) b1_count)/nevents;
+	double b2eff = ((double) b2_count)/nevents;
+	double b1b2eff = b1b2_count/nevents;
+	double b1eeff = b1e_count/nevents;
+	double b2eeff = b2e_count/nevents;
+	double b1b2eeff = b1b2e_count/nevents;
+	TParameter<double> eject_eff("Ejectile Efficiency", eeff);
+	TParameter<double> break1_eff("Breakup 1 Efficiency", b1eff);
+	TParameter<double> break2_eff("Breakup 2 Efficiency", b2eff);
+	TParameter<double> break1break2_eff("Breakup1 Coincident with Breakup2", b1b2eff);
+	TParameter<double> break1eject_eff("Breakup1 Coincident with Ejectile", b1eeff);
+	TParameter<double> break2eject_eff("Breakup2 Coincident with Ejectile", b2eeff);
+	TParameter<double> break1break2eject_eff("All Three Particles Coincident", b1b2eeff);
+
+
+	input->cd();
+	table->Write();
+	eject_eff.Write();
+	break1_eff.Write();
+	break2_eff.Write();
+	break1break2_eff.Write();
+	break1eject_eff.Write();
+	break2eject_eff.Write();
+	break1break2eject_eff.Write();
+
+	input->Close();
+}
+
+void AnasenEfficiency::RunThreeStep(const std::string& filename) {
+	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
+	TTree* tree = (TTree*) input->Get("DataTree");
+
+	THashTable* table = new THashTable();
+
+	Mask::NucData* break1 = new Mask::NucData();
+	Mask::NucData* break3 = new Mask::NucData();
+	Mask::NucData* break4 = new Mask::NucData();
+
+
+	tree->SetBranchAddress("breakup1", &break1);
+	tree->SetBranchAddress("breakup3", &break3);
+	tree->SetBranchAddress("breakup4", &break4);
+
+	double nevents = tree->GetEntries();
+
+	//Progress tracking
+	int percent5 = nevents*0.05;
+	int count = 0;
+	int npercent = 0;
+
+	bool break1_flag, break3_flag, break4_flag;
+	int b1_count=0, b3_count=0, b4_count=0;
+	int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0;
+
+	for(int i=0; i<tree->GetEntries(); i++) {
+		if(++count == percent5) {//Show progress every 5%
+			npercent++;
+			count = 0;
+			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
+		}
+
+		break1_flag = false;
+		break3_flag = false;
+		break4_flag = false;
+
+		tree->GetEntry(i);
+
+		if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) {
+			b1_count++;
+			break1_flag = true;
+			MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE);
+		}
+
+
+		if(break3->KE >= threshold && (IsRing1(break3->theta, break3->phi) || IsRing2(break3->theta, break3->phi) || IsQQQ(break3->theta, break3->phi))) {
+			b3_count++;
+			break3_flag = true;
+			MyFill(table, "detected_break3_theta_ke","detected breakup3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/deg2rad,300,0.0,30.0,break3->KE);
+		}
+
+		if(break4->KE >= threshold && (IsRing1(break4->theta, break4->phi) || IsRing2(break4->theta, break4->phi) || IsQQQ(break4->theta, break4->phi))) {
+			b4_count++;
+			break4_flag = true;
+			MyFill(table, "detected_break4_theta_ke","detected breakup4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/deg2rad,300,0.0,30.0,break4->KE);
+		}
+
+		if(break1_flag && break3_flag && break4_flag) {
+			b1b3b4_count++;
+			b1b3_count++;
+			b1b4_count++;
+			b3b4_count++;
+			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad);
+			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE);
+			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad);
+		} else if(break1_flag && break3_flag) {
+			b1b3_count++;
+		} else if(break1_flag && break4_flag) {
+			b1b4_count++;
+		} else if(break3_flag && break4_flag) {
+			b3b4_count++;
+			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad);
+			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE);
+			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad);
+		}
+	}
+
+	double b1eff = ((double) b1_count)/nevents;
+	double b3eff = ((double) b3_count)/nevents;
+	double b4eff = ((double) b4_count)/nevents;
+	double b1b3eff = b1b3_count/nevents;
+	double b1b4eff = b1b4_count/nevents;
+	double b3b4eff = b3b4_count/nevents;
+	double b1b3b4eff = b1b3b4_count/nevents;
+	TParameter<double> break1_eff("Breakup1 Efficiency", b1eff);
+	TParameter<double> break3_eff("Breakup3 Efficiency", b3eff);
+	TParameter<double> break4_eff("Breakup4 Efficiency", b4eff);
+	TParameter<double> b1b3_eff("Break1 Coincident with Break3", b1b3eff);
+	TParameter<double> b1b4_eff("Break1 Coincident with Break4", b1b4eff);
+	TParameter<double> b3b4_eff("Break3 Coincident with Break4", b3b4eff);
+	TParameter<double> b1b3b4_eff("All Breakups Coincident", b1b3b4eff);
+
+	input->cd();
+	break1_eff.Write();
+	break3_eff.Write();
+	break4_eff.Write();
+	b1b3_eff.Write();
+	b1b4_eff.Write();
+	b3b4_eff.Write();
+	b1b3b4_eff.Write();
+	table->Write();
+	input->Close();
+}

src/AngularDistribution.cpp

@@ -0,0 +1,104 @@
+#include "AngularDistribution.h"
+#include <fstream>
+#include <cmath>
+#include <iostream>
+#include "LegendrePoly.h"
+
+AngularDistribution::AngularDistribution() :
+	generator(nullptr), branchingRatio(1.0), L(0), isoFlag(true)
+{
+}
+
+AngularDistribution::AngularDistribution(const std::string& file) :
+	generator(nullptr), branchingRatio(1.0), L(0), isoFlag(true)
+{
+	ReadDistributionFile(file);
+}
+
+AngularDistribution::~AngularDistribution() {}
+
+void AngularDistribution::ReadDistributionFile(const std::string& file) {
+
+	if(file == "none" || file == "") {
+		L=0;
+		branchingRatio=1.0;
+		constants.clear();
+		constants.push_back(0.5);
+		isoFlag = true;
+		return;
+	}
+
+	std::ifstream input(file);
+	std::string junk;
+	int l;
+	double par;
+
+	if(!input.is_open()) {
+		std::cerr<<"Unable to open distribution file. All values reset to default."<<std::endl;
+		L=0;
+		branchingRatio=1.0;
+		constants.clear();
+		constants.push_back(0.5);
+		isoFlag = true;
+		return;
+	}
+
+	input>>junk>>l;
+	while(input>>junk) {
+		input>>par;
+		constants.push_back(par);
+	}
+	input.close();
+
+	if(constants.size() != ((unsigned int) l+1)) {
+		std::cerr<<"Unexpected number of constants for given angular momentum! Expected "<<l+1<<" and given "<<constants.size()<<std::endl;
+		std::cerr<<"Setting all values to default."<<std::endl;
+		branchingRatio=1.0;
+		constants.clear();
+		constants.push_back(0.5);
+		isoFlag = true;
+		return;
+	}
+
+	//Total branching ratio
+	branchingRatio = constants[0]*2.0;
+	L = l;
+
+	//Renormalize distribution such that total prob is 1.0.
+	//Test branching ratio to see if we "make" a decay particle,
+	//then use re-normalized distribution to pick an angle. 
+	if(constants[0] < 0.5) {
+		double norm = 0.5/constants[0];
+		for(auto& value : constants)
+			value *= norm;
+	}
+
+	isoFlag = false;
+
+}
+
+double AngularDistribution::GetRandomCosTheta() {
+
+	if(!IsGeneratorSet()) {
+		std::cerr<<"Random number generator is not set in AngularDistribution! Returning default value of 0"<<std::endl;
+		return 0.0;
+	}
+
+	if(isoFlag) return generator->Uniform(-1.0, 1.0);
+
+	double test, probability;
+	double costheta;
+
+	test = generator->Uniform(0.0, 1.0);
+	if(test > branchingRatio) return -10;
+
+	do {
+		probability = 0.0;
+		costheta = generator->Uniform(-1.0, 1.0);
+		test = generator->Uniform(0.0, 1.0);
+		for(unsigned int i=0; i<constants.size(); i++)
+			probability += constants[i]*P_l(i*2, costheta);
+	} while(test > probability);
+
+	return costheta;
+}

src/EnergyLoss.cpp

@@ -29,7 +29,7 @@ void EnergyLoss::SetTargetComponents(std::vector<int>& Zt, std::vector<int>& At,
   for(unsigned int i=0; i<Stoich.size(); i++) {
     comp_denom += Stoich[i];
     if(ZT[i] > MAX_Z) {
-      throw ELossException();
+      throw ELossException("Maximum allowed target Z exceeded");
     }
   }
   targ_composition.resize(Stoich.size());
@@ -136,7 +136,7 @@ double EnergyLoss::GetElectronicStoppingPower(double energy) {
   double e_per_u = energy/MP;
   std::vector<double> values;
   if(e_per_u > MAX_H_E_PER_U) {
-    throw ELossException();
+    throw ELossException("Exceeded maximum allowed energy per nucleon");
   } else if (e_per_u > 1000.0) {
     for(auto& z: ZT) {
       values.push_back(Hydrogen_dEdx_High(e_per_u, energy, z));
@@ -150,11 +150,11 @@ double EnergyLoss::GetElectronicStoppingPower(double energy) {
       values.push_back(Hydrogen_dEdx_Low(e_per_u, z));
     }
   } else {
-    throw ELossException();
+    throw ELossException("Negative energy per nucleon");
   }
 
   if(values.size() == 0) {
-    throw ELossException();
+    throw ELossException("Size of value array is 0. Unable to iterate over target components");
   }
 
   if(ZP > 1) { //not hydrogen, need to account for effective charge
@@ -192,6 +192,9 @@ double EnergyLoss::GetNuclearStoppingPower(double energy) {
 
 /*Wrapper function for aquiring total stopping (elec + nuc)*/
 double EnergyLoss::GetTotalStoppingPower(double energy) {
+  if(ZP == 0) {
+    return GetNuclearStoppingPower(energy);
+  }
   return GetElectronicStoppingPower(energy)+GetNuclearStoppingPower(energy);
 }
 

src/Kinematics.cpp

@@ -1,4 +1,5 @@
 #include "Kinematics.h"
+#include "MaskFile.h"
 #include <fstream>
 #include <iostream>
 
@@ -109,25 +110,31 @@ bool Kinematics::LoadConfig(const char* filename) {
 		input>>junk;
 	}
 
-	double par1, par2, L1, L2;
+	double par1, par2;
+	std::string dfile1, dfile2;
 	getline(input, junk);
 	getline(input, junk);
 
 	input>>junk>>m_nsamples;
 	input>>junk>>par1>>junk>>par2;
 	sys->SetBeamDistro(par1, par2);
+	input>>junk>>par1;
+	sys->SetReactionThetaType(par1);
 	input>>junk>>par1>>junk>>par2;
 	sys->SetTheta1Range(par1, par2);
 	input>>junk>>par1>>junk>>par2;
+	sys->SetPhi1Range(par1, par2);
+	input>>junk>>par1>>junk>>par2;
 	sys->SetExcitationDistro(par1, par2);
-	input>>junk>>L1;
+	input>>junk>>dfile1;
-	input>>junk>>L2;
+	input>>junk>>dfile2;
-	sys->SetDecay1AngularMomentum(L1);
+	sys->SetDecay1Distribution(dfile1);
-	sys->SetDecay2AngularMomentum(L2);
+	sys->SetDecay2Distribution(dfile2);
 	sys->SetRandomGenerator(global_generator);
 
 	std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
-	std::cout<<"Decay1 angular momentum: "<<L1<<" Decay2 angular momentum: "<<L2<<std::endl;
+	std::cout<<"Decay1 angular momentum: "<<sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<sys->GetDecay2AngularMomentum()<<std::endl;
+	std::cout<<"Decay1 total branching ratio: "<<sys->GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<<sys->GetDecay2BranchingRatio()<<std::endl;
 	std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
 
 	return true;
@@ -180,6 +187,7 @@ void Kinematics::Run() {
 
 void Kinematics::RunOneStepRxn() {
 
+	
 	TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
 	TTree* tree;
 	NucData targ, proj, eject, residual;
@@ -190,12 +198,15 @@ void Kinematics::RunOneStepRxn() {
 		tree->Branch("ejectile", &eject);
 		tree->Branch("residual", &residual);
 	}
+	
 
 	//For progress tracking
 	int percent5 = 0.05*m_nsamples;
 	int count = 0;
 	int npercent = 0;
 
+	std::vector<Nucleus> data;
+	data.resize(4);
 	for(int i=0; i<m_nsamples; i++) {
 		if(++count == percent5) {//Show update every 5 percent
 			npercent++;
@@ -204,6 +215,7 @@ void Kinematics::RunOneStepRxn() {
 		}
 
 		sys->RunSystem();
+		
 		if(save_tree_flag) {
 			targ = ConvertNucleus(sys->GetTarget());
 			proj = ConvertNucleus(sys->GetProjectile());
@@ -217,8 +229,10 @@ void Kinematics::RunOneStepRxn() {
 			plotman.FillData(sys->GetEjectile());
 			plotman.FillData(sys->GetResidual());
 		}
+
 	}
 
+	
 	output->cd();
 	if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite);
 	if(do_plotter_flag) {
@@ -226,6 +240,7 @@ void Kinematics::RunOneStepRxn() {
 		plotman.ClearTable();
 	}
 	output->Close();
+	
 }
 
 void Kinematics::RunOneStepDecay() {
@@ -264,6 +279,7 @@ void Kinematics::RunOneStepDecay() {
 			plotman.FillData(sys->GetEjectile());
 			plotman.FillData(sys->GetResidual());
 		}
+
 	}
 
 	output->cd();
@@ -282,6 +298,7 @@ void Kinematics::RunTwoStep() {
 		return;
 	}
 
+	
 	TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
 	TTree* tree;
 	NucData targ, proj, eject, residual, break1, break2;
@@ -294,6 +311,8 @@ void Kinematics::RunTwoStep() {
 		tree->Branch("breakup1", &break1);
 		tree->Branch("breakup2", &break2);
 	}
+	
+
 
 	//For progress tracking
 	int percent5 = 0.05*m_nsamples;
@@ -325,8 +344,10 @@ void Kinematics::RunTwoStep() {
 			plotman.FillData(this_sys->GetBreakup1());
 			plotman.FillData(this_sys->GetBreakup2());
 		}
+		
 	}
 
+	
 	output->cd();
 	if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite);
 	if(do_plotter_flag) {
@@ -334,6 +355,7 @@ void Kinematics::RunTwoStep() {
 		plotman.ClearTable();
 	}
 	output->Close();
+	
 }
 
 void Kinematics::RunThreeStep() {

src/LegendrePoly.cpp

@@ -1,4 +1,5 @@
 #include "LegendrePoly.h"
+#include <cmath>
 
 double P_l(int l, double x) {
 	if(l == 0) {
@@ -6,6 +7,26 @@ double P_l(int l, double x) {
 	} else if (l == 1) {
 		return x;
 	} else {
-		return ((2.0*l + 1.0)*x*P_l(l-1, x) - (l-1.0)*P_l(l-2, x))/(double(l));
+		return (2.0*l - 1.0)/l*x*P_l(l-1, x) - (l-1.0)/l*P_l(l-2, x);
 	}
+}
+
+double Normed_P_l_sq(int l, double x) {
+	return (2.0*l+1.0)/2.0*std::pow(P_l(l, x), 2.0);
+}
+
+double P_0(double x) {
+	return 1.0;
+}
+
+double P_1(double x) {
+	return x;
+}
+
+double P_2(double x) {
+	return 0.5*(3.0*x*x -1.0);
+}
+
+double P_l_ROOT(double* x, double* pars) {
+	return P_l(pars[0], x[0]);
 }

src/MaskFile.cpp

@@ -0,0 +1,179 @@
+#include "MaskFile.h"
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+
+namespace Mask {
+
+MaskFile::MaskFile() : 
+	file_type(MaskFile::FileType::none), filename(""), buffer_position(0), count(0), file()
+{
+}
+
+MaskFile::MaskFile(const std::string& name, MaskFile::FileType type) :
+	file_type(type), filename(name), buffer_position(0), count(0), file()
+{
+	Open(filename, type);
+}
+
+bool MaskFile::Open(const std::string& name, MaskFile::FileType type) {
+	if(IsOpen()) {
+		std::cerr<<"Attempted to open file that is already open!"<<std::endl;
+		return false;
+	}
+
+	if(type == FileType::read) {
+		file.open(name, std::ios::in);
+		buffer_position = 0;
+		if(IsOpen()) ReadHeader();
+	} else if (type == FileType::write) {
+		file.open(name, std::ios::out | std::ios::trunc);
+		file<<std::setprecision(precision);
+		file<<std::scientific;
+	} else if (type == FileType::append) {
+		file.open(name, std::ios::out | std::ios::app);
+		file<<std::setprecision(precision);
+		file<<std::scientific;
+	} else {
+		std::cerr<<"Invalid FileType at MaskFile::Open()"<<std::endl;
+		return IsOpen();
+	}
+
+	count = 0;
+
+	return IsOpen();
+}
+
+void MaskFile::Close() {
+	if(data_buffer.size() > 0) {
+		FlushBuffer();
+	}
+	file.close();
+}
+
+void MaskFile::WriteHeader(std::vector<Nucleus>& nuclei) {
+	file<<std::setw(width)<<"ZT"<<","
+		<<std::setw(width)<<"AT"<<","
+		<<std::setw(width)<<"ZP"<<","
+		<<std::setw(width)<<"AP"<<","
+		<<std::setw(width)<<"ZE"<<","
+		<<std::setw(width)<<"AE"<<","
+		<<std::setw(width)<<"ZR"<<","
+		<<std::setw(width)<<"AR"<<","
+		<<std::setw(width)<<"ZB1"<<","
+		<<std::setw(width)<<"AB1"<<","
+		<<std::setw(width)<<"ZB2"<<","
+		<<std::setw(width)<<"AB2"<<","
+		<<std::setw(width)<<"ZB3"<<","
+		<<std::setw(width)<<"AB3"<<","
+		<<std::setw(width)<<"ZB4"<<","
+		<<std::setw(width)<<"AB4"<<","
+		<<std::setw(width)<<"KEP"<<","
+		<<std::endl;
+
+	file<<std::setw(width)<<nuclei[0].GetZ()<<","
+		<<std::setw(width)<<nuclei[0].GetA()<<","
+		<<std::setw(width)<<nuclei[1].GetZ()<<","
+		<<std::setw(width)<<nuclei[1].GetA()<<","
+		<<std::setw(width)<<nuclei[2].GetZ()<<","
+		<<std::setw(width)<<nuclei[2].GetA()<<","
+		<<std::setw(width)<<nuclei[3].GetZ()<<","
+		<<std::setw(width)<<nuclei[3].GetA()<<",";
+	if(nuclei.size() == 6) {
+		file<<std::setw(width)<<nuclei[4].GetZ()<<","
+			<<std::setw(width)<<nuclei[4].GetA()<<","
+			<<std::setw(width)<<nuclei[5].GetZ()<<","
+			<<std::setw(width)<<nuclei[5].GetA()<<",";
+	} else if(nuclei.size() == 8) {
+		file<<std::setw(width)<<nuclei[4].GetZ()<<","
+			<<std::setw(width)<<nuclei[4].GetA()<<","
+			<<std::setw(width)<<nuclei[5].GetZ()<<","
+			<<std::setw(width)<<nuclei[5].GetA()<<","
+			<<std::setw(width)<<nuclei[6].GetZ()<<","
+			<<std::setw(width)<<nuclei[6].GetA()<<","
+			<<std::setw(width)<<nuclei[7].GetZ()<<","
+			<<std::setw(width)<<nuclei[7].GetA()<<",";
+	}
+	file<<std::setw(width)<<nuclei[1].GetKE()<<",";
+	file<<std::endl;
+
+	file<<std::setw(width)<<"Event"<<","
+		<<std::setw(width)<<"EE"<<","<<std::setw(width)<<"KEE"<<","<<std::setw(width)<<"PE"<<","<<std::setw(width)<<"ThetaE"<<","<<std::setw(width)<<"PhiE"<<","<<std::setw(width)<<"ThetaCME"<<","
+		<<std::setw(width)<<"ER"<<","<<std::setw(width)<<"KER"<<","<<std::setw(width)<<"PR"<<","<<std::setw(width)<<"ThetaR"<<","<<std::setw(width)<<"PhiR"<<","<<std::setw(width)<<"ThetaCMR"<<","
+		<<std::setw(width)<<"EB1"<<","<<std::setw(width)<<"KEB1"<<","<<std::setw(width)<<"PB1"<<","<<std::setw(width)<<"ThetaB1"<<","<<std::setw(width)<<"PhiB1"<<","<<std::setw(width)<<"ThetaCMB1"<<","
+		<<std::setw(width)<<"EB2"<<","<<std::setw(width)<<"KEB2"<<","<<std::setw(width)<<"PB2"<<","<<std::setw(width)<<"ThetaB2"<<","<<std::setw(width)<<"PhiB2"<<","<<std::setw(width)<<"ThetaCMB2"<<","
+		<<std::setw(width)<<"EB3"<<","<<std::setw(width)<<"KEB3"<<","<<std::setw(width)<<"PB3"<<","<<std::setw(width)<<"ThetaB3"<<","<<std::setw(width)<<"PhiB3"<<","<<std::setw(width)<<"ThetaCMB3"<<","
+		<<std::setw(width)<<"EB4"<<","<<std::setw(width)<<"KEB4"<<","<<std::setw(width)<<"PB4"<<","<<std::setw(width)<<"ThetaB4"<<","<<std::setw(width)<<"PhiB4"<<","<<std::setw(width)<<"ThetaCMB4"<<","
+		<<std::endl;
+}
+
+void MaskFile::ReadHeader() {
+	std::string junk;
+	std::getline(file, junk);
+}
+
+void MaskFile::WriteData(std::vector<Nucleus>& data) {
+	if(count == 0 && data_buffer.size() == 0) {
+		WriteHeader(data);
+	}
+	
+	data_buffer.push_back(data);
+	if(data_buffer.size() == buffersize) {
+		FlushBuffer();
+	}
+}
+
+void MaskFile::FlushBuffer() {
+	for(auto& event : data_buffer) {
+		file<<std::setw(width)<<count<<",";
+		for(unsigned int i=2; i< event.size(); i++) {
+			Nucleus& nuc = event[i];
+			file<<std::setw(width)<<nuc.GetE()<<","
+				<<std::setw(width)<<nuc.GetKE()<<","
+				<<std::setw(width)<<nuc.GetP()<<","
+				<<std::setw(width)<<nuc.GetTheta()<<","
+				<<std::setw(width)<<nuc.GetPhi()<<","
+				<<std::setw(width)<<nuc.GetThetaCM()<<",";
+		}
+		file<<std::endl;
+		count++;
+	}
+	data_buffer.clear();
+}
+
+std::vector<Nucleus> MaskFile::ReadData() {
+	if(buffer_position == data_buffer.size()) {
+		FillBuffer();
+		buffer_position = 0;
+	}
+	std::vector<Nucleus> data = data_buffer[buffer_position];
+	buffer_position++;
+	
+	return data;
+}
+
+void MaskFile::FillBuffer() {
+	std::vector<Nucleus> data;
+	std::string line;
+	std::stringstream linebuf;
+	std::string junk;
+	int Z, A;
+	double E, p, theta, phi, tcm;
+	data_buffer.clear();
+
+	while(data_buffer.size() <= buffersize) {
+		if(!std::getline(file, line)) break;
+		linebuf.str(line);
+		while(linebuf>>Z) {
+			linebuf>>A;
+			linebuf>>E>>junk>>junk>>p>>theta>>phi>>tcm;
+			data.emplace_back(Z, A);
+			data[data.size()-1].SetVectorSpherical(theta, phi, p, E);
+			data[data.size()-1].SetThetaCM(tcm);
+		}
+		data_buffer.push_back(data);
+		data.clear();
+	}
+}
+
+};

src/Plotter.cpp

@@ -25,12 +25,12 @@ void Plotter::FillData(const Nucleus& nuc) {
 	std::string ex_name = sym + "_ex";
 	std::string ex_title = ex_name + ";E_{ex} (MeV);counts";
 	std::string angdist_name = sym+"_angDist";
-	std::string angdist_title = angdist_name+";#theta_{CM};counts";
+	std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts";
 	
 	MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.GetTheta()*rad2deg, nuc.GetKE(), 2);
 	MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4);
 	MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
-	MyFill(angdist_name.c_str(), angdist_title.c_str(),180,0,180,nuc.GetThetaCM()*rad2deg);
+	MyFill(angdist_name.c_str(), angdist_title.c_str(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
 	
 }
 

src/QQQDetector.cpp

@@ -0,0 +1,155 @@
+#include "QQQDetector.h"
+
+QQQDetector::QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x, double y) :
+	m_Rinner(R_in), m_Router(R_out), m_deltaPhi(deltaPhi), m_phiCentral(phiCentral), m_translation(x,y,z)
+{
+	m_deltaR = (m_Router - m_Rinner)/nrings;
+	m_deltaPhi_per_wedge = m_deltaPhi/nwedges;
+	m_ZRot.SetAngle(m_phiCentral);
+	m_ringCoords.resize(nrings);
+	m_wedgeCoords.resize(nwedges);
+	for(auto& ring : m_ringCoords) {
+		ring.resize(4);
+	}
+	for(auto& wedge : m_wedgeCoords) {
+		wedge.resize(4);
+	}
+
+	CalculateCorners();
+}
+
+QQQDetector::~QQQDetector() {}
+
+void QQQDetector::CalculateCorners() {
+	double x0, x1, x2, x3;
+	double y0, y1, y2, y3;
+	double z0, z1, z2, z3;
+
+	//Generate flat ring corner coordinates
+	for(int i=0; i<nrings; i++) {
+		x0 = (m_Rinner + m_deltaR*(i+1))*std::cos(-m_deltaPhi/2.0);
+		y0 = (m_Rinner + m_deltaR*(i+1))*std::sin(-m_deltaPhi/2.0);
+		z0 = 0.0;
+		m_ringCoords[i][0].SetVectorCartesian(x0, y0, z0);
+
+		x1 = (m_Rinner + m_deltaR*(i))*std::cos(-m_deltaPhi/2.0);
+		y1 = (m_Rinner + m_deltaR*(i))*std::sin(-m_deltaPhi/2.0);
+		z1 = 0.0;
+		m_ringCoords[i][1].SetVectorCartesian(x1, y1, z1);
+
+		x2 = (m_Rinner + m_deltaR*(i))*std::cos(m_deltaPhi/2.0);
+		y2 = (m_Rinner + m_deltaR*(i))*std::sin(m_deltaPhi/2.0);
+		z2 = 0.0;
+		m_ringCoords[i][2].SetVectorCartesian(x2, y2, z2);
+
+		x3 = (m_Rinner + m_deltaR*(i+1))*std::cos(m_deltaPhi/2.0);
+		y3 = (m_Rinner + m_deltaR*(i+1))*std::sin(m_deltaPhi/2.0);
+		z3 = 0.0;
+		m_ringCoords[i][3].SetVectorCartesian(x3, y3, z3);
+	}
+
+	//Generate flat wedge corner coordinates
+	for(int i=0; i<nwedges; i++) {
+		x0 = m_Router * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
+		y0 = m_Router * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
+		z0 = 0.0;
+		m_wedgeCoords[i][0].SetVectorCartesian(x0, y0, z0);
+
+		x1 = m_Rinner * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
+		y1 = m_Rinner * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
+		z1 = 0.0;
+		m_wedgeCoords[i][1].SetVectorCartesian(x1, y1, z1);
+
+		x2 = m_Rinner * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
+		y2 = m_Rinner * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
+		z2 = 0.0;
+		m_wedgeCoords[i][2].SetVectorCartesian(x2, y2, z2);
+
+		x3 = m_Router * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
+		y3 = m_Router * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
+		z3 = 0.0;
+		m_wedgeCoords[i][3].SetVectorCartesian(x3, y3, z3);
+	}
+
+	for(int i=0; i<nrings; i++) {
+		for(int j=0; j<4; j++) {
+			m_ringCoords[i][j] = TransformCoordinates(m_ringCoords[i][j]);
+		}
+	}
+
+	for(int i=0; i<nwedges; i++) {
+		for(int j=0; j<4; j++) {
+			m_wedgeCoords[i][j] = TransformCoordinates(m_wedgeCoords[i][j]);
+		}
+	}
+
+}
+
+Mask::Vec3 QQQDetector::GetTrajectoryCoordinates(double theta, double phi) {
+	double z_to_detector = m_translation.GetZ();
+	double rho_traj = z_to_detector*std::tan(theta);
+	double r_traj = std::sqrt(rho_traj*rho_traj + z_to_detector*z_to_detector);
+	double min_rho, max_rho, min_phi, max_phi;
+
+	Mask::Vec3 result;
+
+	for(auto& ring : m_ringCoords) {
+		min_rho = ring[1].GetRho();
+		max_rho = ring[0].GetRho();
+		if(rho_traj < max_rho && rho_traj > min_rho) {
+			for(auto& wedge : m_wedgeCoords) {
+				min_phi = wedge[0].GetPhi();
+				max_phi = wedge[3].GetPhi();
+				if(phi < min_phi && phi < max_phi) {
+					result.SetVectorSpherical(r_traj, theta, phi);
+					break;
+				}
+			}
+		}
+	}
+	
+	return result;
+
+}
+
+std::pair<int,int> QQQDetector::GetTrajectoryRingWedge(double theta, double phi) {
+	double z_to_detector = m_translation.GetZ();
+	double rho_traj = z_to_detector*std::tan(theta);
+	double min_rho, max_rho, min_phi, max_phi;
+
+
+	for(int r=0; r<nrings; r++) {
+		auto& ring = m_ringCoords[r];
+		min_rho = ring[1].GetRho();
+		max_rho = ring[0].GetRho();
+		if(rho_traj < max_rho && rho_traj > min_rho) {
+			for(int w=0; w<nwedges; w++) {
+				auto& wedge = m_wedgeCoords[w];
+				min_phi = wedge[0].GetPhi();
+				max_phi = wedge[3].GetPhi();
+				if(phi > min_phi && phi < max_phi) {
+					return std::make_pair(r, w);
+				}
+			}
+		}
+	}
+	
+	return std::make_pair(-1, -1);
+
+}
+
+Mask::Vec3 QQQDetector::GetHitCoordinates(int ringch, int wedgech) {
+	if(!CheckChannel(ringch) || !CheckChannel(wedgech)) {
+		return Mask::Vec3();
+	}
+
+	double r_center  = m_Rinner + (0.5+ringch)*m_deltaR;
+	double phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhi_per_wedge;
+	double x = r_center*std::cos(phi_center);
+	double y = r_center*std::sin(phi_center);
+	double z = 0;
+
+	Mask::Vec3 hit(x, y, z);
+
+	return TransformCoordinates(hit);
+}

src/Reaction.cpp

@@ -12,12 +12,12 @@
 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)
+	target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), nuc_initFlag(false), resid_elossFlag(false)
 {
 }
 
 Reaction::Reaction(int zt, int at, int zp, int ap, int ze, int ae) :
-	target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), resid_elossFlag(false)
+	target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), resid_elossFlag(false)
 {
 	SetNuclei(zt, at, zp, ap, ze, ae);
 }
@@ -77,16 +77,28 @@ void Reaction::SetBeamKE(double bke) {
 	m_bke = bke - target->GetProjectileEnergyLoss(reactants[1].GetZ(), reactants[1].GetA(), bke, rxnLayer, 0);
 };
 
+void Reaction::SetEjectileThetaType(int type) {
+	if(decayFlag) return;
+	if(type != center_of_mass && type != lab) return;
+
+	m_eject_theta_type = type;
+}
+
 //Methods given by Iliadis in Nuclear Physics of Stars, Appendix C
-void Reaction::CalculateReaction() {
+//For use with lab frame restricted angles. May not give appropriate disribution for ejectile
-	//Target assumed at rest, with 0 excitation energy
+void Reaction::CalculateReactionThetaLab() {
 	reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass());
 	double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass()));
 	double beam_E = m_bke + reactants[1].GetGroundStateMass();
 	reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E);
 
-
 	double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex);
+
+	double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass());
+	if(m_bke < Ethresh) {
+		throw EnergyThresholdException();
+	}
+
 	double term1 = sqrt(reactants[1].GetGroundStateMass()*reactants[2].GetGroundStateMass()*m_bke)/(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())*cos(m_theta);
 	double term2 = (m_bke*(reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()) + reactants[3].GetGroundStateMass()*Q)/(reactants[3].GetGroundStateMass() + reactants[2].GetGroundStateMass());
 	double sqrt_pos_ejectKE = term1 + sqrt(term1*term1 + term2);
@@ -97,7 +109,6 @@ void Reaction::CalculateReaction() {
 	} else {
 		ejectKE = sqrt_neg_ejectKE*sqrt_neg_ejectKE;
 	}
-
 	double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass()));
 	double ejectE = ejectKE + reactants[2].GetGroundStateMass();
 
@@ -105,12 +116,59 @@ void Reaction::CalculateReaction() {
 
 	reactants[3] = reactants[0] + reactants[1] - reactants[2];
 
-	//energy loss for ejectile (after reaction!)
 	ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, m_theta);
 	ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass()));
 	ejectE = ejectKE + reactants[2].GetGroundStateMass();
 	reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE);
 
+	if(resid_elossFlag) {
+		double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta());
+		double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass()));
+		double residE = residKE + reactants[3].GetInvMass();
+		reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE);
+	}
+}
+
+//Methods from original ANASEN. Gives proper distribution for inverse kinematics.
+void Reaction::CalculateReactionThetaCM() {
+	//Target assumed at rest, with 0 excitation energy
+	reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass());
+	double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass()));
+	double beam_E = m_bke + reactants[1].GetGroundStateMass();
+	reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E);
+
+
+	double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex);
+
+	double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass());
+	if(m_bke < Ethresh) {
+		throw EnergyThresholdException();
+	}
+	
+	auto parent = reactants[0] + reactants[1];
+	double boost2lab[3];
+	double boost2cm[3];
+	const double* boost = parent.GetBoost();
+	for(int i=0; i<3; i++) {
+		boost2lab[i] = boost[i];
+		boost2cm[i] = boost[i]*(-1.0);
+	}
+	parent.ApplyBoost(boost2cm);
+	double ejectE_cm = (std::pow(reactants[2].GetGroundStateMass(), 2.0) - std::pow(reactants[3].GetGroundStateMass() + m_ex, 2.0) + std::pow(parent.GetE(),2.0))/(2.0*parent.GetE());
+	double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - std::pow(reactants[2].GetGroundStateMass(), 2.0));
+	reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm);
+	reactants[2].ApplyBoost(boost2lab);
+	reactants[3] = reactants[0] + reactants[1] - reactants[2];
+
+	double ejectKE = reactants[2].GetKE();
+	double ejectP = reactants[2].GetP();
+	double ejectE = reactants[2].GetE();
+	//energy loss for ejectile (after reaction!)
+	ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, reactants[2].GetTheta());
+	ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass()));
+	ejectE = ejectKE + reactants[2].GetGroundStateMass();
+	reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE);
+
 	//if on, get eloss for residual (after reaction!)
 	if(resid_elossFlag) {
 		double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta());
@@ -118,7 +176,21 @@ void Reaction::CalculateReaction() {
 		double residE = residKE + reactants[3].GetInvMass();
 		reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE);
 	}
+}
 
+void Reaction::CalculateReaction() {
+	switch(m_eject_theta_type) {
+		case center_of_mass:
+		{
+			CalculateReactionThetaCM();
+			break;
+		}
+		case lab:
+		{
+			CalculateReactionThetaLab();
+			break;
+		}
+	}
 }
 
 //Calculate in CM, where decay is isotropic

src/ReactionSystem.cpp

@@ -5,12 +5,12 @@
 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("")
+	m_beamDist(0,0), m_theta1Range(0,0), m_phi1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("")
 {
 }
 
 ReactionSystem::ReactionSystem(std::vector<int>& z, std::vector<int>& a) :
-	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("")
+	m_beamDist(0,0), m_theta1Range(0,0), m_phi1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("")
 {
 	SetNuclei(z, a);
 }
@@ -29,6 +29,13 @@ bool ReactionSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
 	return true;
 }
 
+void ReactionSystem::SetRandomGenerator(TRandom3* gen) {
+	generator = gen;
+	decay1dist.AttachRandomNumberGenerator(gen);
+	decay2dist.AttachRandomNumberGenerator(gen);
+	gen_set_flag = true;
+}
+
 void ReactionSystem::AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness) {
 	target.AddLayer(zt, at, stoich, thickness);
 }
@@ -55,22 +62,6 @@ void ReactionSystem::SetSystemEquation()  {
 	m_sys_equation += step1.GetResidual().GetIsotopicSymbol();
 }
 
-double ReactionSystem::GetDecayTheta(int L) {
-	if(!gen_set_flag) return 0.0;
-
-	double probability = 0.0;
-	double costheta = 0.0;
-	double check = 0.0;
-
-	do {
-		costheta = generator->Uniform(-1.0, 1.0);
-		check = generator->Uniform(0.0, 1.0);
-		probability = std::pow(P_l(L, costheta), 2.0);
-	} while(check > probability);
-
-	return std::acos(costheta);
-}
-
 void ReactionSystem::RunSystem() {
 	if(!gen_set_flag) return;
 	
@@ -80,7 +71,13 @@ void ReactionSystem::RunSystem() {
 	}
 
 	if(step1.IsDecay()) {
-		double rxnTheta = GetDecayTheta(L1);
+		double decaycostheta = decay1dist.GetRandomCosTheta();
+		if(decaycostheta == -10.0) {
+			step1.ResetEjectile();
+			step1.ResetResidual();
+			return;
+		}
+		double rxnTheta = std::acos(decay1dist.GetRandomCosTheta());
 		double rxnPhi = generator->Uniform(0, 2.0*M_PI);
 		step1.SetPolarRxnAngle(rxnTheta);
 		step1.SetAzimRxnAngle(rxnPhi);
@@ -91,7 +88,7 @@ void ReactionSystem::RunSystem() {
 		//Sample parameters
 		double bke = generator->Gaus(m_beamDist.first, m_beamDist.second);
 		double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second)));
-		double rxnPhi = 0;
+		double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second);
 		double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
 	
 		step1.SetBeamKE(bke);

src/SabreEfficiency.cpp

@@ -1,5 +1,4 @@
 #include "SabreEfficiency.h"
-#include "Kinematics.h"
 #include <fstream>
 #include <iostream>
 #include <TFile.h>
@@ -10,9 +9,10 @@
 #include <TCanvas.h>
 #include <TH2.h>
 #include <TH1.h>
+#include <TCutG.h>
 
 SabreEfficiency::SabreEfficiency() : 
-	m_rxn_type(-1), deadlayer(DEADLAYER_THIN)
+	m_rxn_type(-1), deadlayer(DEADLAYER_THIN), sabre_eloss(SABRE_THICKNESS)
 {
 	detectors.reserve(5);
 	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI0*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
@@ -21,59 +21,39 @@ 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);
 
- 	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++) {
- 				coords = detectors[i].GetRingTiltCoords(j, k);
- 				ringxs.push_back(coords.GetX());
- 				ringys.push_back(coords.GetY());
- 				ringzs.push_back(coords.GetZ());
- 			}
- 		}
- 	}
-
- 	for(int i=0; i<5; i++) {
- 		for(int j=0; j<detectors[i].GetNumberOfWedges(); j++) {
- 			for(int k=0; k<4; k++) {
- 				coords = detectors[i].GetWedgeTiltCoords(j, k);
- 				wedgexs.push_back(coords.GetX());
- 				wedgeys.push_back(coords.GetY());
- 				wedgezs.push_back(coords.GetZ());
- 			}
- 		}
- 	}
 	std::vector<int> dead_z = {14};
 	std::vector<int> dead_a = {28};
 	std::vector<int> dead_stoich = {1};
 	deadlayer.SetElements(dead_z, dead_a, dead_stoich);
+	sabre_eloss.SetElements(dead_z, dead_a, dead_stoich);
 }
 
 SabreEfficiency::~SabreEfficiency() {}
 
-void SabreEfficiency::MyFill(THashTable* table, const char* name, const char* title, int bins, float min, float max, double val) {
+void SabreEfficiency::MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val) {
-	TH1F* h = (TH1F*) table->FindObject(name);
+	TH1F* h = (TH1F*) table->FindObject(name.c_str());
 	if(h) {
 		h->Fill(val);
 	} else {
-		h = new TH1F(name, title, bins, min, max);
+		h = new TH1F(name.c_str(), title.c_str(), bins, min, max);
 		h->Fill(val);
 		table->Add(h);
 	}
 }
 
-void SabreEfficiency::MyFill(THashTable* table, const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy) {
+void SabreEfficiency::MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy) {
-	TH2F* h = (TH2F*) table->FindObject(name);
+	TH2F* h = (TH2F*) table->FindObject(name.c_str());
 	if(h) {
 		h->Fill(valx, valy);
 	} else {
-		h = new TH2F(name, title, binsx, minx, maxx, binsy, miny, maxy);
+		h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy);
 		h->Fill(valx, valy);
 		table->Add(h);
 	}
 }
 
-void SabreEfficiency::CalculateEfficiency(const char* file) {
+void SabreEfficiency::CalculateEfficiency(const std::string& file) {
 	std::cout<<"----------SABRE Efficiency Calculation----------"<<std::endl;
 	std::cout<<"Loading in output from kinematics simulation: "<<file<<std::endl;
 	std::cout<<"Running efficiency calculation..."<<std::endl;
@@ -107,90 +87,35 @@ void SabreEfficiency::CalculateEfficiency(const char* file) {
 	std::cout<<"---------------------------------------------"<<std::endl;
 }
 
-void SabreEfficiency::RunDecay(const char* file) {
+void SabreEfficiency::DrawDetectorSystem(const std::string& filename) {
-	TFile* input = TFile::Open(file, "UPDATE");
+	TFile* output = TFile::Open(filename.c_str(), "RECREATE");
-	TTree* tree = (TTree*) input->Get("DataTree");
-
-	Mask::NucData* eject = new Mask::NucData();
-	Mask::NucData* resid = new Mask::NucData();
-
-
-	tree->SetBranchAddress("ejectile", &eject);
-	tree->SetBranchAddress("residual", &resid);
-
-	double nevents = tree->GetEntries();
-	std::vector<double> resid_xs, eject_xs;
-	std::vector<double> resid_ys, eject_ys;
-	std::vector<double> resid_zs, eject_zs;
-
-	//Progress tracking
-	int percent5 = nevents*0.05;
-	int count = 0;
-	int npercent = 0;
 
+	std::vector<double> ringxs, ringys, ringzs;
+    std::vector<double> wedgexs, wedgeys, wedgezs;
 	Mask::Vec3 coords;
-	double thetaIncident, eloss;
+ 	for(int i=0; i<5; i++) {
+ 		for(int j=0; j<detectors[i].GetNumberOfRings(); j++) {
+ 			for(int k=0; k<4; k++) {
+ 				coords = detectors[i].GetRingTiltCoords(j, k);
+ 				ringxs.push_back(coords.GetX());
+ 				ringys.push_back(coords.GetY());
+ 				ringzs.push_back(coords.GetZ());
+ 			}
+ 		}
+ 	}
 
-	for(int i=0; i<tree->GetEntries(); i++) {
+ 	for(int i=0; i<5; i++) {
-		if(++count == percent5) {//Show progress every 5%
+ 		for(int j=0; j<detectors[i].GetNumberOfWedges(); j++) {
-			npercent++;
+ 			for(int k=0; k<4; k++) {
-			count = 0;
+ 				coords = detectors[i].GetWedgeTiltCoords(j, k);
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
+ 				wedgexs.push_back(coords.GetX());
-		}
+ 				wedgeys.push_back(coords.GetY());
+ 				wedgezs.push_back(coords.GetZ());
+ 			}
+ 		}
+ 	}
 
-		tree->GetEntry(i);
+ 	TGraph2D* gr = new TGraph2D(ringxs.size(), &(ringxs[0]), &(ringys[0]), &(ringzs[0]));
-
-		if(eject->KE >= ENERGY_THRESHOLD) {
-			for(int j=0; j<5; j++) {
-				auto& det = detectors[j];
-				auto chan = det.GetTrajectoryRingWedge(eject->theta, eject->phi);
-				if(chan.first != -1 && chan.second != -1) {
-					if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break; //dead channel check
-					coords = det.GetTrajectoryCoordinates(eject->theta, eject->phi);
-					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eloss = deadlayer.getEnergyLossTotal(eject->Z, eject->A, eject->KE, M_PI - thetaIncident);
-					if((eject->KE - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check
-
-					eject_xs.push_back(coords.GetX());
-					eject_ys.push_back(coords.GetY());
-					eject_zs.push_back(coords.GetZ());
-					break;
-				}
-			}
-		}
-
-		if(resid->KE > ENERGY_THRESHOLD) {
-			for(int j=0; j<5; j++) {
-				auto& det = detectors[j];
-				auto chan = det.GetTrajectoryRingWedge(resid->theta, resid->phi);
-				if(chan.first != -1 && chan.second != -1) {
-					if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break;
-					coords = det.GetTrajectoryCoordinates(resid->theta, resid->phi);
-					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eloss = deadlayer.getEnergyLossTotal(resid->Z, resid->A, resid->KE, M_PI - thetaIncident);
-					if((resid->KE - eloss) <= ENERGY_THRESHOLD) break;
-
-					resid_xs.push_back(coords.GetX());
-					resid_ys.push_back(coords.GetY());
-					resid_zs.push_back(coords.GetZ());
-					break;
-				}
-			}
-		}
-
-	}
-
-	double ejecteff = ((double) eject_xs.size())/nevents;
-	double resideff = ((double) resid_xs.size())/nevents;
-	TParameter<double> eject_eff("Light Breakup Efficiency", ejecteff);
-	TParameter<double> resid_eff("Heavy Breakup Efficiency", resideff);
-
-	TGraph2D* gde = new TGraph2D(eject_xs.size(), &(eject_xs[0]), &(eject_ys[0]), &(eject_zs[0]));
-	gde->SetName("detected_eject_points");
-	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(1);
@@ -201,26 +126,127 @@ void SabreEfficiency::RunDecay(const char* file) {
 	gw->SetMarkerStyle(1);
 
 	TCanvas* canvas = new TCanvas();
-	canvas->SetName("detectors_and_particles");
+	canvas->SetName("detector_system");
 	canvas->cd();
 	gr->Draw("AP");
 	gw->Draw("same P");
-	gde->Draw("same P");
+
+	canvas->Write();
+	gr->Write();
+	gw->Write();
+
+	output->Close();
+}
+
+double SabreEfficiency::RunConsistencyCheck() {
+	double theta, phi;
+	double npoints = 5.0*16.0*4.0;
+	int count=0;
+ 	for(int h=0; h<5; h++) {
+ 		for(int j=0; j<16; j++) {
+ 			for(int k=0; k<4; k ++) {
+ 				theta  = detectors[h].GetRingTiltCoords(j, k).GetTheta();
+ 				phi = detectors[h].GetRingTiltCoords(j, k).GetPhi();
+ 				for(int i=0; i<5; i++) {
+ 					auto channels = detectors[i].GetTrajectoryRingWedge(theta, phi);
+ 					if(channels.first != -1) {
+ 						count++;
+ 					}
+ 				}
+ 			}
+ 		}
+ 	}
+
+ 	return ((double)count)/npoints;
+}
+
+/*Returns if detected, as well as total energy deposited in SABRE*/
+std::pair<bool,double> SabreEfficiency::IsSabre(Mask::NucData* nucleus) {
+	if(nucleus->KE <= ENERGY_THRESHOLD) {
+		return std::make_pair(false, 0.0);
+	}
+
+	Mask::Vec3 coords;
+	double thetaIncident, eloss, e_deposited;
+	for(int i=0; i<5; i++) {
+		auto chan = detectors[i].GetTrajectoryRingWedge(nucleus->theta, nucleus->phi);
+		if(chan.first != -1 && chan.second != -1) {
+			if(dmap.IsDead(i, chan.first, 0) || dmap.IsDead(i, chan.second, 1)) break; //dead channel check
+			coords = detectors[i].GetTrajectoryCoordinates(nucleus->theta, nucleus->phi);
+			thetaIncident = std::acos(coords.Dot(detectors[i].GetNormTilted())/(coords.GetR()));
+			eloss = deadlayer.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE, M_PI - thetaIncident);
+			if((nucleus->KE - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check
+			e_deposited = sabre_eloss.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE - eloss, M_PI - thetaIncident);
+			return std::make_pair(true, e_deposited);
+		}
+	}
+
+	return std::make_pair(false,0.0);
+}
+
+void SabreEfficiency::RunDecay(const std::string& filename) {
+	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
+	TTree* tree = (TTree*) input->Get("DataTree");
+
+	THashTable* table = new THashTable();
+
+	Mask::NucData* eject = new Mask::NucData();
+	Mask::NucData* resid = new Mask::NucData();
+
+
+	tree->SetBranchAddress("ejectile", &eject);
+	tree->SetBranchAddress("residual", &resid);
+
+	double nevents = tree->GetEntries();
+
+	//Progress tracking
+	int percent5 = nevents*0.05;
+	int count = 0;
+	int npercent = 0;
+
+	std::pair<bool, double> eject_result, resid_result;
+	int detected_eject=0, detected_resid=0;
+	for(int i=0; i<tree->GetEntries(); i++) {
+		if(++count == percent5) {//Show progress every 5%
+			npercent++;
+			count = 0;
+			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
+		}
+
+		tree->GetEntry(i);
+
+		eject_result = IsSabre(eject);
+		resid_result = IsSabre(resid);
+
+		if(eject_result.first) {
+			detected_eject++;
+			MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject->KE);
+			MyFill(table, "detected_eject_theta_edep","detected ejectile theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject_result.second);
+		}
+
+		if(resid_result.first) {
+			detected_resid++;
+			MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid->KE);
+			MyFill(table, "detected_resid_theta_edep","detected residual theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid_result.second);
+		}
+	}
+
+	double ejecteff = ((double) detected_eject)/nevents;
+	double resideff = ((double) detected_resid)/nevents;
+	TParameter<double> eject_eff("Ejectile Efficiency", ejecteff);
+	TParameter<double> resid_eff("Residual Efficiency", resideff);
 
 	input->cd();
 	eject_eff.Write();
 	resid_eff.Write();
-	gr->Write();
+	table->Write();
-	gw->Write();
-	gde->Write();
-	canvas->Write();
 	input->Close();
 
 }
 
-void SabreEfficiency::Run2Step(const char* file) {
+void SabreEfficiency::Run2Step(const std::string& filename) {
 
-	TFile* input = TFile::Open(file, "UPDATE");
+	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
 	TTree* tree = (TTree*) input->Get("DataTree");
 
 	Mask::NucData* break1 = new Mask::NucData();
@@ -232,20 +258,17 @@ void SabreEfficiency::Run2Step(const char* file) {
 	tree->SetBranchAddress("breakup2", &break2);
 
 	double nevents = tree->GetEntries();
-	std::vector<double> b1_thetas, b2_thetas;
-	std::vector<double> b1_phis, b2_phis;
-	std::vector<double> b1_kes, b2_kes;
 
 	//Progress tracking
 	int percent5 = nevents*0.05;
 	int count = 0;
 	int npercent = 0;
 
-	int cnt_09=0, cnt_08=0, cnt_07=0, cnt_06=0;
 	double costheta_cm =0;
 
-	Mask::Vec3 coords;
+	std::pair<bool, double> break1_result, break2_result;
-	double thetaIncident, eloss;
+	int detected_b1=0, detected_b2=0;
+
 	for(int i=0; i<tree->GetEntries(); i++) {
 		if(++count == percent5) {//Show progress every 5%
 			npercent++;
@@ -255,97 +278,44 @@ void SabreEfficiency::Run2Step(const char* file) {
 
 		tree->GetEntry(i);
 
-		if(break1->KE >= ENERGY_THRESHOLD) {
+		break1_result = IsSabre(break1);
-			for(int j=0; j<5; j++) {
+		break2_result = IsSabre(break2);
-				auto& det = detectors[j];
-				auto chan = det.GetTrajectoryRingWedge(break1->theta, break1->phi);
-				if(chan.first != -1 && chan.second != -1) {
-					if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break;
-					costheta_cm = std::cos(break1->theta_cm);
-					double this_bin;
-					std::string this_name;
-					for(int k=-9; k<=10; k++) {
-						this_bin = 0.1*k;
-						this_name = "counter_histo_bin"+std::to_string(k);
-						if(costheta_cm > this_bin-0.05 && costheta_cm < this_bin+0.05) {
-							MyFill(table, this_name.c_str(), "ctheta",20,-1.0,1.0,costheta_cm);
-							break;
-						} 
-					}
-					/*
-					if(costheta_cm > -0.95 && costheta_cm < -0.85) {
-						cnt_09++;
-					} else if(costheta_cm > -0.85 && costheta_cm < -0.75) {
-						cnt_08++;
-					} else if(costheta_cm > -0.75 && costheta_cm < -0.65) {
-						cnt_07++;
-					} else if(costheta_cm > -0.65 && costheta_cm < -0.55) {
-						cnt_06++;
-					}
-					*/
-					coords = det.GetTrajectoryCoordinates(break1->theta, break1->phi);
-					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eloss = deadlayer.getEnergyLossTotal(break1->Z, break1->A, break1->KE, M_PI - thetaIncident);
-					if((break1->KE - eloss) <= ENERGY_THRESHOLD) break;
-					
 
-					b1_thetas.push_back(break1->theta);
+		if(break1_result.first) {
-					b1_phis.push_back(break1->phi);
+			detected_b1++;
-					b1_kes.push_back(break1->KE);
+			costheta_cm = std::cos(break1->theta_cm);
-					break;
+			MyFill(table,"detected_break1_cm_theta","cos(#theta_{CM})",20,-1.0,1.0,costheta_cm);
-				}
+			MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE);
-			}
+			MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second);
 		}
 
-		if(break2->KE > ENERGY_THRESHOLD) {
+		if(break2_result.first) {
-			for(int j=0; j<5; j++) {
+			detected_b2++;
-				auto& det = detectors[j];
+			MyFill(table, "detected_break2_theta_ke","detected break2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2->KE);
-				auto chan = det.GetTrajectoryRingWedge(break2->theta, break2->phi);
+			MyFill(table, "detected_break2_theta_edep","detected break2 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2_result.second);
-				if(chan.first != -1 && chan.second != -1) {
+		
-					if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break;
-					coords = det.GetTrajectoryCoordinates(break2->theta, break2->phi);
-					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eloss = deadlayer.getEnergyLossTotal(break2->Z, break2->A, break2->KE, M_PI - thetaIncident);
-					if((break2->KE - eloss) <= ENERGY_THRESHOLD) break;
-
-					b2_thetas.push_back(break2->theta);
-					b2_phis.push_back(break2->phi);
-					b2_kes.push_back(break2->KE);
-					break;
-				}
-			}
 		}
 	}
 
-	double b1eff = ((double) b1_thetas.size())/nevents;
+	double b1eff = ((double) detected_b1)/nevents;
-	double b2eff = ((double) b2_thetas.size())/nevents;
+	double b2eff = ((double) detected_b2)/nevents;
-	double b1eff_09 = cnt_09/nevents;
+	TParameter<double> break1_eff("Breakup1 Efficiency", b1eff);
-	double b1eff_08 = cnt_08/nevents;
+	TParameter<double> break2_eff("Breakup2 Efficiency", b2eff);
-	double b1eff_07 = cnt_07/nevents;
-	double b1eff_06 = cnt_06/nevents;
-	TParameter<double> break1_eff("Light Breakup Efficiency", b1eff);
-	TParameter<double> break1_eff_09("Light Breakup Efficiency CosTheta -0.9", b1eff_09);
-	TParameter<double> break1_eff_08("Light Breakup Efficiency CosTheta -0.8", b1eff_08);
-	TParameter<double> break1_eff_07("Light Breakup Efficiency CosTheta -0.7", b1eff_07);
-	TParameter<double> break1_eff_06("Light Breakup Efficiency CosTheta -0.6", b1eff_06);
-	TParameter<double> break2_eff("Heavy Breakup Efficiency", b2eff);
 
 	input->cd();
 	table->Write();
 	break1_eff.Write();
-	break1_eff_09.Write();
-	break1_eff_08.Write();
-	break1_eff_07.Write();
-	break1_eff_06.Write();
 	break2_eff.Write();
 	input->Close();
 
 }
 
-void SabreEfficiency::Run3Step(const char* file) {
+void SabreEfficiency::Run3Step(const std::string& filename) {
-	TFile* input = TFile::Open(file, "UPDATE");
+	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
 	TTree* tree = (TTree*) input->Get("DataTree");
 
+	THashTable* table = new THashTable();
+
 	Mask::NucData* break1 = new Mask::NucData();
 	Mask::NucData* break3 = new Mask::NucData();
 	Mask::NucData* break4 = new Mask::NucData();
@@ -356,23 +326,16 @@ void SabreEfficiency::Run3Step(const char* file) {
 	tree->SetBranchAddress("breakup4", &break4);
 
 	double nevents = tree->GetEntries();
-	std::vector<double> b1_thetas, b3_thetas, b4_thetas;
-	std::vector<double> b1_phis, b3_phis, b4_phis;
-	std::vector<double> b1_kes, b3_kes, b4_kes;
 
 	//Progress tracking
 	int percent5 = nevents*0.05;
 	int count = 0;
 	int npercent = 0;
 
-	bool break1_flag, break3_flag, break4_flag;
+	std::pair<bool, double> break1_result, break3_result, break4_result;
+	int detected_b1=0, detected_b3=0, detected_b4=0;
 	int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0;
 
-	TH2F* b3b4_thetas = new TH2F("b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0,180,180,0,180);
-	TH2F* b3b4_kes = new TH2F("b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0,15,150,0,15);
-	TH2F* b3b4_phis = new TH2F("b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0,360,360,0,360);
-	Mask::Vec3 coords;
-	double thetaIncident, eloss;
 
 	for(int i=0; i<tree->GetEntries(); i++) {
 		if(++count == percent5) {//Show progress every 5%
@@ -381,100 +344,61 @@ void SabreEfficiency::Run3Step(const char* file) {
 			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
 		}
 
-		break1_flag = false;
+		break1_result = IsSabre(break1);
-		break3_flag = false;
+		break3_result = IsSabre(break3);
-		break4_flag = false;
+		break4_result = IsSabre(break4);
+
 
 		tree->GetEntry(i);
 
-		if(break1->KE > ENERGY_THRESHOLD) {
+		if(break1_result.first) {
-			for(int j=0; j<5; j++) {
+			detected_b1++;
-				auto& det = detectors[j];
+			MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE);
-				auto chan = det.GetTrajectoryRingWedge(break1->theta, break1->phi);
+			MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second);
-				if(chan.first != -1 && chan.second != -1) {
-					coords = det.GetTrajectoryCoordinates(break1->theta, break1->phi);
-					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eloss = deadlayer.getEnergyLossTotal(break1->Z, break1->A, break1->KE, M_PI - thetaIncident);
-					if((break1->KE - eloss) <= ENERGY_THRESHOLD) break;
-
-					b1_thetas.push_back(break1->theta);
-					b1_phis.push_back(break1->phi);
-					b1_kes.push_back(break1->KE);
-					break1_flag = true;
-					break;
-				}
-			}
 		}
 
-
+		if(break3_result.first) {
-		if(break3->KE > ENERGY_THRESHOLD) {
+			detected_b3++;
-			for(int j=0; j<5; j++) {
+			MyFill(table, "detected_break3_theta_ke","detected break3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3->KE);
-				auto& det = detectors[j];
+			MyFill(table, "detected_break3_theta_edep","detected break3 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3_result.second);
-				auto chan = det.GetTrajectoryRingWedge(break3->theta, break3->phi);
-				if(chan.first != -1 && chan.second != -1) {
-					coords = det.GetTrajectoryCoordinates(break3->theta, break3->phi);
-					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eloss = deadlayer.getEnergyLossTotal(break3->Z, break3->A, break3->KE, M_PI - thetaIncident);
-					if((break3->KE - eloss) <= ENERGY_THRESHOLD) break;
-
-					b3_thetas.push_back(break3->theta);
-					b3_phis.push_back(break3->phi);
-					b3_kes.push_back(break3->KE);
-					break3_flag = true;
-					break;
-				}
-			}
 		}
 
-		if(break4->KE > ENERGY_THRESHOLD) {
+		if(break4_result.first) {
-			for(int j=0; j<5; j++) {
+			detected_b1++;
-				auto& det = detectors[j];
+			MyFill(table, "detected_break4_theta_ke","detected break4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4->KE);
-				auto chan = det.GetTrajectoryRingWedge(break4->theta, break4->phi);
+			MyFill(table, "detected_break4_theta_edep","detected break4 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4_result.second);
-				if(chan.first != -1 && chan.second != -1) {
-					coords = det.GetTrajectoryCoordinates(break4->theta, break4->phi);
-					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eloss = deadlayer.getEnergyLossTotal(break4->Z, break4->A, break4->KE, M_PI - thetaIncident);
-					if((break4->KE - eloss) <= ENERGY_THRESHOLD) break;
-
-					b4_thetas.push_back(break4->theta);
-					b4_phis.push_back(break4->phi);
-					b4_kes.push_back(break4->KE);
-					break4_flag = true;
-					break;
-				}
-			}
 		}
 
-		if(break1_flag && break3_flag && break4_flag) {
+		if(break1_result.first && break3_result.first && break4_result.first) {
 			b1b3b4_count++;
 			b1b3_count++;
 			b1b4_count++;
 			b3b4_count++;
-			b3b4_thetas->Fill(b3_thetas[b3_thetas.size()-1]/DEG2RAD, b4_thetas[b4_thetas.size()-1]/DEG2RAD);
+			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD);
-			b3b4_kes->Fill(b3_kes[b3_kes.size()-1], b4_kes[b4_kes.size()-1]);
+			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE);
-			b3b4_phis->Fill(b3_phis[b3_phis.size()-1]/DEG2RAD, b4_phis[b4_phis.size()-1]/DEG2RAD);
+			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD);
-		} else if(break1_flag && break3_flag) {
+		} else if(break1_result.first && break3_result.first) {
 			b1b3_count++;
-		} else if(break1_flag && break4_flag) {
+		} else if(break1_result.first && break4_result.first) {
 			b1b4_count++;
-		} else if(break3_flag && break4_flag) {
+		} else if(break3_result.first && break4_result.first) {
 			b3b4_count++;
-			b3b4_thetas->Fill(b3_thetas[b3_thetas.size()-1]/DEG2RAD, b4_thetas[b4_thetas.size()-1]/DEG2RAD);
+			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD);
-			b3b4_kes->Fill(b3_kes[b3_kes.size()-1], b4_kes[b4_kes.size()-1]);
+			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE);
-			b3b4_phis->Fill(b3_phis[b3_phis.size()-1]/DEG2RAD, b4_phis[b4_phis.size()-1]/DEG2RAD);
+			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD);
 		}
 	}
 
-	double b1eff = ((double) b1_thetas.size())/nevents;
+	double b1eff = ((double) detected_b1)/nevents;
-	double b3eff = ((double) b3_thetas.size())/nevents;
+	double b3eff = ((double) detected_b3)/nevents;
-	double b4eff = ((double) b4_thetas.size())/nevents;
+	double b4eff = ((double) detected_b4)/nevents;
 	double b1b3eff = b1b3_count/nevents;
 	double b1b4eff = b1b4_count/nevents;
 	double b3b4eff = b3b4_count/nevents;
 	double b1b3b4eff = b1b3b4_count/nevents;
-	TParameter<double> break1_eff("Light Initial Breakup Efficiency", b1eff);
+	TParameter<double> break1_eff("Breakup1 Efficiency", b1eff);
-	TParameter<double> break3_eff("Light Final Breakup Efficiency", b3eff);
+	TParameter<double> break3_eff("Breakup3 Efficiency", b3eff);
-	TParameter<double> break4_eff("Heavy Final Breakup Efficiency", b4eff);
+	TParameter<double> break4_eff("Breakup4 Efficiency", b4eff);
 	TParameter<double> b1b3_eff("Break1 Coincident with Break3", b1b3eff);
 	TParameter<double> b1b4_eff("Break1 Coincident with Break4", b1b4eff);
 	TParameter<double> b3b4_eff("Break3 Coincident with Break4", b3b4eff);
@@ -488,8 +412,7 @@ void SabreEfficiency::Run3Step(const char* file) {
 	b1b4_eff.Write();
 	b3b4_eff.Write();
 	b1b3b4_eff.Write();
-	b3b4_thetas->Write();
+	table->Write();
-	b3b4_phis->Write();
+
-	b3b4_kes->Write();
 	input->Close();
 }

src/StripDetector.cpp

@@ -0,0 +1,141 @@
+#include "StripDetector.h"
+
+/*
+  Corner layout for each strip in the un-rotated frame
+  0--------------------------1
+  |                          |              y
+  |                          |              ^
+  |                          |              |
+  |                          |              |
+  |                          |              |
+  2--------------------------3    z<--------X
+                                            x
+*/
+
+StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) :
+	m_random(nullptr)
+{
+
+	num_strips = ns;
+	
+	length = std::fabs(len);
+	total_width = std::fabs(wid);
+	front_strip_width = total_width/num_strips;
+	back_strip_length = length/num_strips;
+
+	while (cphi < 0) cphi += 2*M_PI;
+	center_phi = cphi;
+	center_z   = cz;
+	center_rho = std::fabs(crho);
+
+	zRot.SetAngle(center_phi);
+
+	front_strip_coords.resize(num_strips);
+	back_strip_coords.resize(num_strips);
+	rotated_front_strip_coords.resize(num_strips);
+	rotated_back_strip_coords.resize(num_strips);
+	for(int i=0; i<num_strips; i++) {
+		front_strip_coords[i].resize(num_corners);
+		back_strip_coords[i].resize(num_corners);
+		rotated_front_strip_coords[i].resize(num_corners);
+		rotated_back_strip_coords[i].resize(num_corners);
+	}
+	CalculateCorners();
+	
+}
+
+StripDetector::~StripDetector() {}
+
+void StripDetector::CalculateCorners() {
+	double y_min, y_max, z_min, z_max; 
+	for (int s=0; s<num_strips; s++) {
+		y_max = -total_width/2.0 + front_strip_width*(s+1);
+		y_min = -total_width/2.0 + front_strip_width*s;
+		z_max = center_z + length/2.0;
+		z_min = center_z - length/2.0;
+		front_strip_coords[s][0] = Mask::Vec3(center_rho, y_max, z_max);
+		front_strip_coords[s][1] = Mask::Vec3(center_rho, y_max, z_min);
+		front_strip_coords[s][2] = Mask::Vec3(center_rho, y_min, z_max);
+		front_strip_coords[s][3] = Mask::Vec3(center_rho, y_min, z_min);
+
+		z_max = (center_z - length/2.0) + (s+1)*back_strip_length;
+		z_min = (center_z - length/2.0) + (s)*back_strip_length;
+		y_max = total_width/2.0;
+		y_min = -total_width/2.0;
+		back_strip_coords[s][0] = Mask::Vec3(center_rho, y_max, z_max);
+		back_strip_coords[s][1] = Mask::Vec3(center_rho, y_max, z_min);
+		back_strip_coords[s][2] = Mask::Vec3(center_rho, y_min, z_max);
+		back_strip_coords[s][3] = Mask::Vec3(center_rho, y_min, z_min);
+	}
+
+	for(int s=0; s<num_strips; s++) {
+		rotated_front_strip_coords[s][0] = zRot*front_strip_coords[s][0];
+		rotated_front_strip_coords[s][1] = zRot*front_strip_coords[s][1];
+		rotated_front_strip_coords[s][2] = zRot*front_strip_coords[s][2];
+		rotated_front_strip_coords[s][3] = zRot*front_strip_coords[s][3];
+
+		rotated_back_strip_coords[s][0] = zRot*back_strip_coords[s][0];
+		rotated_back_strip_coords[s][1] = zRot*back_strip_coords[s][1];
+		rotated_back_strip_coords[s][2] = zRot*back_strip_coords[s][2];
+		rotated_back_strip_coords[s][3] = zRot*back_strip_coords[s][3];
+	}
+}
+
+Mask::Vec3 StripDetector::GetHitCoordinates(int front_stripch, double front_strip_ratio) {
+
+	if (!ValidChannel(front_stripch) || !ValidRatio(front_strip_ratio)) return Mask::Vec3(0,0,0);
+
+	double y;
+	//If we have a random number generator, randomize y position within pixel. Otherwise take halfway.
+	if(m_random) {
+		y = -total_width/2.0 + (front_stripch + m_random->Uniform(0.0, 1.0))*front_strip_width;
+	} else {
+		y = -total_width/2.0 + (front_stripch+0.5)*front_strip_width;
+	}
+
+	//recall we're still assuming phi=0 det:
+	Mask::Vec3 coords(center_rho, y, front_strip_ratio*(length/2) + center_z);
+  
+	//NOW rotate by appropriate phi
+	return zRot*coords;
+
+}
+
+std::pair<int,double> StripDetector::GetChannelRatio(double theta, double phi) {
+
+	while (phi < 0) phi += 2*M_PI;
+
+	//to make the math easier (and the same for each det), rotate the input phi
+	//BACKWARD by the phi of the det, s.t. we are looking along x-axis
+	phi -= center_phi;
+
+	if (phi > M_PI) phi -= 2*M_PI;
+
+	//then we can check easily whether it even hit the detector in phi
+	double det_max_phi = atan2(total_width/2,center_rho);
+	double det_min_phi = -det_max_phi;
+  
+	if (phi < det_min_phi || phi > det_max_phi) return std::make_pair(-1,0);
+
+	//for theta it's not so simple, so we have to go through the typical plane-intersect method
+	//first thing's first: we have a fixed x for the entire detector plane:
+	double xhit = center_rho;
+	//thus we find the corresponding y and z for that fixed x, given the input theta and phi:
+	double yhit = xhit*tan(phi);
+	double zhit = sqrt(xhit*xhit+yhit*yhit)/tan(theta);
+
+	int front_ch_hit=-1;
+	double ratio=0;
+	for (int s=0; s<num_strips; s++) {
+		if (xhit >= front_strip_coords[s][0].GetX() && xhit <= front_strip_coords[s][0].GetX() && //Check min and max x (constant in flat)
+			yhit >= front_strip_coords[s][2].GetY() && yhit <= front_strip_coords[s][0].GetY() && //Check min and max y
+			zhit >= front_strip_coords[s][1].GetZ() && zhit <= front_strip_coords[s][0].GetZ()) //Check min and max z
+		{
+			front_ch_hit = s;
+			ratio = (zhit-center_z)/(length/2);
+			break;
+		}
+	}
+
+	return std::make_pair(front_ch_hit,ratio);
+}

src/ThreeStepSystem.cpp

@@ -79,10 +79,12 @@ void ThreeStepSystem::RunSystem() {
 	//Sample parameters
 	double bke = generator->Gaus(m_beamDist.first, m_beamDist.second);
 	double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second)));
-	double rxnPhi = 0;
+	double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second);
-	double decay1Theta = GetDecayTheta(L1);
+	double decay1costheta = decay1dist.GetRandomCosTheta();
+	double decay1Theta = std::acos(decay1costheta);
 	double decay1Phi = generator->Uniform(0, M_PI*2.0);
-	double decay2Theta = GetDecayTheta(L2);
+	double decay2costheta = decay2dist.GetRandomCosTheta();
+	double decay2Theta = std::acos(decay2costheta);
 	double decay2Phi = generator->Uniform(0, M_PI*2.0);
 	double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
 
@@ -100,9 +102,22 @@ void ThreeStepSystem::RunSystem() {
 	step1.Calculate();
 
 	step2.SetTarget(step1.GetResidual());
+	if(decay1costheta == -10) {
+		step2.ResetEjectile();
+		step2.ResetResidual();
+		step3.ResetTarget();
+		step3.ResetEjectile();
+		step3.ResetResidual();
+		return;
+	}
 	step2.Calculate();
 
 	step3.SetTarget(step2.GetResidual());
+	if(decay2costheta == -10) {
+		step3.ResetEjectile();
+		step3.ResetResidual();
+		return;
+	}
 	step3.TurnOnResidualEloss();
 	step3.Calculate();
 

src/TwoStepSystem.cpp

@@ -70,8 +70,9 @@ void TwoStepSystem::RunSystem() {
 	//Sample parameters
 	double bke = generator->Gaus(m_beamDist.first, m_beamDist.second);
 	double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second)));
-	double rxnPhi = 0;
+	double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second);
-	double decay1Theta = GetDecayTheta(L1);
+	double decay1costheta = decay1dist.GetRandomCosTheta();
+	double decay1Theta = std::acos(decay1costheta);
 	double decay1Phi = generator->Uniform(0, M_PI*2.0);
 	double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
 
@@ -86,6 +87,11 @@ void TwoStepSystem::RunSystem() {
 	step1.Calculate();
 
 	step2.SetTarget(step1.GetResidual());
+	if(decay1costheta == -10) {
+		step2.ResetEjectile();
+		step2.ResetResidual();
+		return;
+	}
 	step2.TurnOnResidualEloss();
 	step2.Calculate();
 

src/main.cpp

@@ -2,15 +2,19 @@
 #include <string>
 #include "Kinematics.h"
 #include "SabreEfficiency.h"
-#include "SabreDetector.h"
+#include "AnasenEfficiency.h"
 #include "KinematicsExceptions.h"
 
+#include <TGraph2D.h>
+#include <TCanvas.h>
+
 int main(int argc, char** argv) {
 	if(argc<2) {
 		std::cerr<<"Incorrect number of arguments!"<<std::endl;
 		return 1;
 	}
 	
+	
 	Mask::Kinematics calculator;
 	try {
 		if(!calculator.LoadConfig(argv[1])) {
@@ -22,53 +26,26 @@ int main(int argc, char** argv) {
 		std::cerr<<"Terminating process."<<std::endl;
 		return 1;
 	}
+	
+	
+
+	/*
 	SabreEfficiency sabre;
 	std::string mapfile = "./etc/DeadChannels.txt";
 	sabre.SetReactionType(calculator.GetReactionType());
 	sabre.SetDeadChannelMap(mapfile);
 	sabre.CalculateEfficiency(calculator.GetOutputName());
+	//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
+	//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.root");
+	*/
 
-	/*std::vector<SabreDetector> detectors;
-	const double INNER_R = 0.0326;
-    const double OUTER_R = 0.1351;
-    const double TILT = 40.0;
-    const double DIST_2_TARG = -0.1245;
-    const double PHI_COVERAGE = 54.4; //delta phi for each det
-    const double PHI0 = 234.0; //center phi values for each det in array
-    const double PHI1 = 162.0; //# is equal to detID in channel map
-    const double PHI2 = 306.0;
-    const double PHI3 = 18.0;
-    const double PHI4 = 90.0;
-    const double DEG2RAD = M_PI/180.0;
 
-    detectors.reserve(5);
+ 	AnasenEfficiency anasen;
-	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI0*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
+ 	anasen.SetReactionType(calculator.GetReactionType());
-	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI1*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
+ 	anasen.CalculateEfficiency(calculator.GetOutputName());
- 	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI2*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
+ 	//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
- 	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
+ 	//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_test.root");
- 	detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
 
- 	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 ++) {
- 				theta  = detectors[h].GetRingTiltCoords(j, k).GetTheta();
- 				phi = detectors[h].GetRingTiltCoords(j, k).GetPhi();
- 				expected_flat_p = detectors[h].GetRingFlatCoords(j, k).GetPhi();  
- 				for(int i=0; i<5; i++) {
- 					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;
- 					}
- 				}
- 			}
- 		}
- 	}*/
 
 	return 0;