Added several key new features. Dead channels are now implemented in a text file in etc. Decays can now have an angular distribution based on inputed l values. New dependancy: gsl

README.md

@@ -19,11 +19,7 @@ By default MASK is capable of simulating reactions of up to three steps. Here is
 2. A reaction of type 2 is a reaction followed by a subsequent decay of the residual nucleus. Again, all sampling is allowed.
 3. A reaction of type 3 is a reaction followed by a subsequent decay of the residual, followed by a decay of one of the products. Again, all sampling is allowed
 
-Note that for type 2 and type 3, the decays are assumed isotropic in the center-of-mass frame of the decay. The input file requires that the user include target information,
+For decays, a specific angular momentum L can be assumed. These are given in the input file as Decay1_AngularMomentum, and Decay2_AngularMomentum. This essentially modifies the center-of-mass angular distribution (as well as the lab frame). It is assumed that the decays in the center-of-mass frame are isotropic in phi (i.e. m=0). Decay1 corresponds to the first decay, if there are multiple steps, Decay2 to the second. If there are no decays, these parameters are not used (or if only one decay, Decay2_AngularMomentum is not used). The input file requires that the user include target information, which will be used to calculate energy loss for all of the reactants and reaction products. The target can contain layers, and each layer can be composed of a compound of elements with a given stoichiometry. If the user wishes to not include energy loss in the kinematics, simply give all target layers a thickness of 0. Note that more layers and more thickness = more time spent calculating energy loss. These energy loss methods are only applicable for solid targets, and should not be applied to gas or liquid targets. Energy loss calculations have a stated uncertainty of approximately five percent.
-which will be used to calculate energy loss for all of the reactants and reaction products. The target can contain layers, and each layer can be composed of a compound of elements
-with a given stoichiometry. If the user wishes to not include energy loss in the kinematics, simply give all target layers a thickness of 0. Note that more layers and more thickness = more
-time spent calculating energy loss. These energy loss methods are only applicable for solid targets, and should not be applied to gas or liquid targets. Energy loss calculations have a
-stated uncertainty of approximately five percent.
 
 The default MASK program includes a calculation of SABRE efficiency, whose methods are contained in the SabreEfficiency and SabreDetector classes. This can be disabled by modifying the
 main.cpp file appropriately.
@@ -37,5 +33,4 @@ To run MASK simply do the following from the MASK directory:
 Input.txt can be replaced by any text file with the correct format.
 
 ## Requirements
-MASK requires that ROOT is installed for data writting and visualization, as well as for random number generation. Testing has been done only on ROOT 6. Mileage on all other ROOT versions
+MASK requires that ROOT is installed for data writting and visualization, as well as for random number generation. It also requires gsl to calculate Legendre Polynomials. Testing has been done only on ROOT 6. Mileage on all other ROOT versions will vary.
-will vary.

etc/DeadChannels.txt

@@ -0,0 +1,4 @@
+DetectorID RING/WEDGE Channel
+4 WEDGE 4
+3 RING 2
+0 RING 13

include/DeadChannelMap.h

@@ -0,0 +1,33 @@
+#ifndef DEADCHANNELMAP_H
+#define DEADCHANNELMAP_H
+
+#include <string>
+#include <iostream>
+#include <unordered_map>
+
+class DeadChannelMap {
+public:
+	DeadChannelMap();
+	DeadChannelMap(std::string& name);
+	~DeadChannelMap();
+	void LoadMapfile(std::string& name);
+	bool IsDead(int det, int channel, int ringwedgeFlag);
+	bool IsValid() { return validFlag; };
+
+private:
+	std::unordered_map<int, bool> dcMap;
+	bool validFlag;
+	static constexpr int RING = 0;
+	static constexpr int WEDGE = 1;
+
+	/*
+	Channel identifier calculated like detector*24 + 16*(RING or WEDGE) + channel
+	Example:
+	Detector 1 ring 15
+	1*24 + 16*0 + 15 = 39
+	Detector 1 wedge 0
+	1*24 + 16*1 +0= 40
+	*/
+};
+
+#endif

include/Kinematics.h

@@ -17,6 +17,7 @@ struct NucData {
 	double Ex = -1;
 	double p = -1;
 	double theta = -1;
+	double theta_cm = -1;
 	double phi = -1;
 	int Z = -1;
 	int A = -1;

include/Nucleus.h

@@ -21,15 +21,18 @@ public:
 	Nucleus(int Z, int A, double px, double py, double pz, double E);
 	~Nucleus();
 	bool SetIsotope(int Z, int A);
+	inline void SetThetaCM(double tcm) { m_theta_cm = tcm; }; //save theta in rxn CM frame
 	inline int GetZ() const { return m_z; };
 	inline int GetA() const { return m_a; };
 	inline double GetExcitationEnergy() const { return GetInvMass() - m_gs_mass; };
 	inline double GetGroundStateMass() const { return m_gs_mass; };
 	inline const char* GetIsotopicSymbol() const { return m_symbol.c_str(); };
+	inline double GetThetaCM() const { return m_theta_cm; };
 
 	inline Nucleus& operator=(const Nucleus& rhs) {
 		SetIsotope(rhs.GetZ(), rhs.GetA());
 		SetVectorCartesian(rhs.GetPx(), rhs.GetPy(), rhs.GetPz(), rhs.GetE());
+		SetThetaCM(rhs.GetThetaCM());
 		return *this;
 	};
 
@@ -46,6 +49,7 @@ public:
 private:
 	int m_z, m_a;
 	double m_gs_mass;
+	double m_theta_cm;
 	std::string m_symbol;
 
 };

include/ReactionSystem.h

@@ -29,6 +29,11 @@ public:
 	inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); };
 	inline void SetTheta1Range(double min, double max) { m_theta1Range = 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 SetDecay2AngularMomentum(double l) { L2 = l; };
+
+	/*Sampling over angular distribution*/
+	double GetDecayTheta(int L);
 
 	virtual void RunSystem();
 
@@ -51,6 +56,7 @@ protected:
 
 	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

@@ -3,12 +3,14 @@
 
 #include "SabreDetector.h"
 #include "Target.h"
+#include "DeadChannelMap.h"
 
 class SabreEfficiency {
 public:
 	SabreEfficiency();
 	~SabreEfficiency();
 	inline void SetReactionType(int t) { m_rxn_type = t; };
+    void SetDeadChannelMap(std::string& filename) { dmap.LoadMapfile(filename); };
 	void CalculateEfficiency(const char* file);
 
 private:
@@ -21,13 +23,13 @@ private:
     std::vector<double> ringxs, ringys, ringzs;
     std::vector<double> wedgexs, wedgeys, wedgezs;
 	Target deadlayer;
+    DeadChannelMap dmap;
 
 
 	//Sabre constants
 	const double INNER_R = 0.0326;
     const double OUTER_R = 0.1351;
     const double TILT = 40.0;
-    //const double DIST_2_TARG = 0.14549;
     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

input.txt

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

makefile

@@ -1,8 +1,8 @@
 CC=g++
 ROOTGEN=rootcint
-CFLAGS=-std=c++11 -g -Wall `root-config --cflags`
+CFLAGS=-g -Wall `root-config --cflags`
 CPPFLAGS=-I ./include
-LDFLAGS=`root-config --glibs`
+LDFLAGS=`root-config --glibs` -lgsl
 
 ROOTINCLDIR=./
 INCLDIR=./include

src/DeadChannelMap.cpp

@@ -0,0 +1,68 @@
+#include "DeadChannelMap.h"
+#include <fstream>
+
+DeadChannelMap::DeadChannelMap() :
+	validFlag(false)
+{
+	int maxchan = 5*24+1*16+7;
+	for(int i=0; i<maxchan; i++) {
+		dcMap[i] = false;
+	}
+}
+
+DeadChannelMap::DeadChannelMap(std::string& name) :
+	validFlag(false)
+{
+	LoadMapfile(name);
+}
+
+DeadChannelMap::~DeadChannelMap() {}
+
+void DeadChannelMap::LoadMapfile(std::string& name) {
+	std::ifstream input(name);
+	if(!input.is_open()) {
+		std::cerr<<"Unable to load dead channels, file: "<<name<<" could not be opened"<<std::endl;
+		validFlag = false;
+		return;
+	}
+
+	
+
+	std::string junk, rw;
+	int detID, channel;
+	int this_channel;
+
+	std::getline(input, junk);
+	while(input>>detID) {
+		input>>rw>>channel;
+		if(rw == "RING") {
+			this_channel = detID*24+RING*16+channel;
+			dcMap[this_channel] = true;
+		} else if(rw == "WEDGE") {
+			this_channel = detID*24+WEDGE*16+channel;
+			dcMap[this_channel] = true;
+		} else {
+			std::cerr<<"Invalid ring/wedge designation at DeadChannelMap"<<std::endl;
+		}
+	}
+
+	validFlag = true;
+}
+
+bool DeadChannelMap::IsDead(int detID, int channel, int ringwedgeFlag) {
+	if(!IsValid() || (ringwedgeFlag != 0 && ringwedgeFlag != 1)) {
+		std::cerr<<"Error at DeadChannelMap IsDead(), bad input parameters"<<std::endl;
+		return false;
+	}
+
+	int this_channel = detID*24+ringwedgeFlag*16+channel;
+
+	auto iter = dcMap.find(this_channel);
+	if(iter != dcMap.end()) {
+		return iter->second;
+	} else {
+		std::cerr<<"Error at DeadChannelMap IsDead(), invalid channel: "<<this_channel<<std::endl;
+		std::cerr<<"detID: "<<detID<<" ringwedgeFlag: "<<ringwedgeFlag<<" channel: "<<channel<<std::endl;
+		return false;
+	}
+}

src/Kinematics.cpp

@@ -109,7 +109,7 @@ bool Kinematics::LoadConfig(const char* filename) {
 		input>>junk;
 	}
 
-	double par1, par2;
+	double par1, par2, L1, L2;
 	getline(input, junk);
 	getline(input, junk);
 
@@ -120,9 +120,14 @@ bool Kinematics::LoadConfig(const char* filename) {
 	sys->SetTheta1Range(par1, par2);
 	input>>junk>>par1>>junk>>par2;
 	sys->SetExcitationDistro(par1, par2);
+	input>>junk>>L1;
+	input>>junk>>L2;
+	sys->SetDecay1AngularMomentum(L1);
+	sys->SetDecay2AngularMomentum(L2);
 	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<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
 
 	return true;
@@ -136,6 +141,7 @@ NucData Kinematics::ConvertNucleus(const Nucleus& nuc) {
 	datum.KE = nuc.GetKE();
 	datum.p = nuc.GetP();
 	datum.theta = nuc.GetTheta();
+	datum.theta_cm = nuc.GetThetaCM();
 	datum.phi = nuc.GetPhi();
 	datum.Ex = nuc.GetExcitationEnergy();
 	datum.Z = nuc.GetZ();

src/Nucleus.cpp

@@ -12,12 +12,12 @@
 namespace Mask {
 
 Nucleus::Nucleus () :
-	Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_symbol("")
+	Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_theta_cm(0), m_symbol("")
 {
 }
 
 Nucleus::Nucleus(int Z, int A) :
-	Vec4(), m_z(Z), m_a(A)
+	Vec4(), m_z(Z), m_a(A), m_theta_cm(0)
 {
 	m_gs_mass = MASS.FindMass(Z, A);
 	m_symbol = MASS.FindSymbol(Z, A);

src/Plotter.cpp

@@ -24,10 +24,13 @@ void Plotter::FillData(const Nucleus& nuc) {
 	std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
 	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";
 	
 	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);
 	
 }
 

src/Reaction.cpp

@@ -143,6 +143,7 @@ void Reaction::CalculateDecay() {
 	double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - reactants[2].GetGroundStateMass()*reactants[2].GetGroundStateMass());
 
 	reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm);
+	reactants[2].SetThetaCM(m_theta);
 
 	reactants[0].ApplyBoost(boost2lab);
 	reactants[2].ApplyBoost(boost2lab);

src/ReactionSystem.cpp

@@ -1,5 +1,6 @@
 #include "ReactionSystem.h"
 #include "KinematicsExceptions.h"
+#include <gsl/gsl_sf_legendre.h>
 
 namespace Mask {
 
@@ -54,6 +55,22 @@ 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(gsl_sf_legendre_Pl(L, costheta), 2.0);
+	} while(check > probability);
+
+	return std::acos(costheta);
+}
+
 void ReactionSystem::RunSystem() {
 	if(!gen_set_flag) return;
 	
@@ -63,7 +80,7 @@ void ReactionSystem::RunSystem() {
 	}
 
 	if(step1.IsDecay()) {
-		double rxnTheta = acos(generator->Uniform(-1, 1));
+		double rxnTheta = GetDecayTheta(L1);
 		double rxnPhi = generator->Uniform(0, 2.0*M_PI);
 		step1.SetPolarRxnAngle(rxnTheta);
 		step1.SetAzimRxnAngle(rxnPhi);

src/SabreEfficiency.cpp

@@ -8,6 +8,8 @@
 #include <TGraph.h>
 #include <TGraph2D.h>
 #include <TCanvas.h>
+#include <TH2.h>
+#include <TH1.h>
 
 SabreEfficiency::SabreEfficiency() : 
 	m_rxn_type(-1), deadlayer(DEADLAYER_THIN)
@@ -54,6 +56,13 @@ void SabreEfficiency::CalculateEfficiency(const char* file) {
 	std::cout<<"Loading in output from kinematics simulation: "<<file<<std::endl;
 	std::cout<<"Running efficiency calculation..."<<std::endl;
 
+	if(!dmap.IsValid()) {
+		std::cerr<<"Unable to run SABRE Efficiency without a dead channel map."<<std::endl;
+		std::cerr<<"If you have no dead channels, simply make a file that's empty"<<std::endl;
+		std::cerr<<"Exiting."<<std::endl;
+		std::cout<<"---------------------------------------------"<<std::endl;
+	}
+
 	switch(m_rxn_type) {
 		case Mask::Kinematics::ONESTEP_DECAY:
 		{
@@ -97,7 +106,8 @@ void SabreEfficiency::RunDecay(const char* file) {
 	int count = 0;
 	int npercent = 0;
 
-	Mask::Vec3 coordinates;
+	Mask::Vec3 coords;
+	double thetaIncident, eloss;
 
 	for(int i=0; i<tree->GetEntries(); i++) {
 		if(++count == percent5) {//Show progress every 5%
@@ -109,24 +119,38 @@ void SabreEfficiency::RunDecay(const char* file) {
 		tree->GetEntry(i);
 
 		if(eject->KE >= ENERGY_THRESHOLD) {
-			for(auto& det : detectors) {
+			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) {
-					coordinates = det.GetTrajectoryCoordinates(eject->theta, eject->phi);
+					if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break; //dead channel check
-					eject_xs.push_back(coordinates.GetX());
+					coords = det.GetTrajectoryCoordinates(eject->theta, eject->phi);
-					eject_ys.push_back(coordinates.GetY());
+					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					eject_zs.push_back(coordinates.GetZ());
+					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(auto& det : detectors) {
+			for(int j=0; j<5; j++) {
-				if(det.GetTrajectoryCoordinates(resid->theta, resid->phi).GetX() != 0) {
+				auto& det = detectors[j];
-					resid_xs.push_back(coordinates.GetX());
+				auto chan = det.GetTrajectoryRingWedge(resid->theta, resid->phi);
-					resid_ys.push_back(coordinates.GetY());
+				if(chan.first != -1 && chan.second != -1) {
-					resid_zs.push_back(coordinates.GetZ());
+					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;
 				}
 			}
@@ -189,15 +213,16 @@ void SabreEfficiency::Run2Step(const char* file) {
 	std::vector<double> b1_phis, b2_phis;
 	std::vector<double> b1_kes, b2_kes;
 
-	double avg_ke_per_pixel[640] = {0};
-	int hits_per_pixel[640] = {0};
-
 	//Progress tracking
 	int percent5 = nevents*0.05;
 	int count = 0;
 	int npercent = 0;
 
+	int cnt_03=0, cnt_47=0, cnt_811=0, cnt_1215=0;
+	double avg_theta03=0, avg_theta47=0, avg_theta811=0, avg_theta1215=0;
+
 	Mask::Vec3 coords;
+	double thetaIncident, eloss;
 	for(int i=0; i<tree->GetEntries(); i++) {
 		if(++count == percent5) {//Show progress every 5%
 			npercent++;
@@ -212,14 +237,27 @@ void SabreEfficiency::Run2Step(const char* file) {
 				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;
+
+					if(chan.first >= 0 && chan.first <4) {
+						cnt_03++;
+						avg_theta03 += break1->theta_cm;
+					} else if(chan.first >= 4 && chan.first <8) {
+						cnt_47++;
+						avg_theta47 += break1->theta_cm;
+					} else if(chan.first >= 8 && chan.first <12) {
+						cnt_811++;
+						avg_theta811 += break1->theta_cm;
+					} else if(chan.first >= 12 && chan.first <16) {
+						cnt_1215++;
+						avg_theta1215 += break1->theta_cm;
+					}
+
 					coords = det.GetTrajectoryCoordinates(break1->theta, break1->phi);
-					double thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
+					thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR()));
-					double eloss = deadlayer.getEnergyLossTotal(break1->Z, break1->A, break1->KE, M_PI - thetaIncident);
+					eloss = deadlayer.getEnergyLossTotal(break1->Z, break1->A, break1->KE, M_PI - thetaIncident);
 					if((break1->KE - eloss) <= ENERGY_THRESHOLD) break;
 					
-					int pixel = (chan.first + 16*chan.second) + 128*j; //calc pixel
-					avg_ke_per_pixel[pixel] += (break1->KE - eloss);
-					hits_per_pixel[pixel]++;
 
 					b1_thetas.push_back(break1->theta);
 					b1_phis.push_back(break1->phi);
@@ -230,8 +268,16 @@ void SabreEfficiency::Run2Step(const char* file) {
 		}
 
 		if(break2->KE > ENERGY_THRESHOLD) {
-			for(auto& det : detectors) {
+			for(int j=0; j<5; j++) {
-				if(det.GetTrajectoryCoordinates(break2->theta, break2->phi).GetX() != 0) {
+				auto& det = detectors[j];
+				auto chan = det.GetTrajectoryRingWedge(break2->theta, break2->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(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);
@@ -239,24 +285,39 @@ void SabreEfficiency::Run2Step(const char* file) {
 				}
 			}
 		}
-
 	}
 
 	double b1eff = ((double) b1_thetas.size())/nevents;
 	double b2eff = ((double) b2_thetas.size())/nevents;
+	double b1eff_ring03 = cnt_03/nevents;
+	double b1eff_ring47 = cnt_47/nevents;
+	double b1eff_ring811 = cnt_811/nevents;
+	double b1eff_ring1215 = cnt_1215/nevents;
+	avg_theta03 /= cnt_03;
+	avg_theta47 /= cnt_47;
+	avg_theta811 /= cnt_811;
+	avg_theta1215 /= cnt_1215;
 	TParameter<double> break1_eff("Light Breakup Efficiency", b1eff);
+	TParameter<double> break1_eff_ring03("Light Breakup Efficiency Rings03", b1eff_ring03);
+	TParameter<double> break1_eff_ring47("Light Breakup Efficiency Rings47", b1eff_ring47);
+	TParameter<double> break1_eff_ring811("Light Breakup Efficiency Rings811", b1eff_ring811);
+	TParameter<double> break1_eff_ring1215("Light Breakup Efficiency Rings1215", b1eff_ring1215);
+	TParameter<double> avg_theta_cm_ring03("Avg. ThetaCM Rings03", avg_theta03);
+	TParameter<double> avg_theta_cm_ring47("Avg. ThetaCM Rings47", avg_theta47);
+	TParameter<double> avg_theta_cm_ring811("Avg. ThetaCM Rings811", avg_theta811);
+	TParameter<double> avg_theta_cm_ring1215("Avg. ThetaCM Rings1215", avg_theta1215);
 	TParameter<double> break2_eff("Heavy Breakup Efficiency", b2eff);
 
-	std::ofstream output("/data1/gwm17/test_dead_pixels.txt");
-	output<<"Average particle kinetic energy (MeV) per pixel (pixel = (ringch + wedgech*16) + 128*detID)"<<std::endl;
-	for(int i=0; i<640; i++) {
-		if(hits_per_pixel[i] == 0) output<<i<<" "<<0.0<<std::endl;
-		else output<<i<<" "<<((double) (avg_ke_per_pixel[i]/hits_per_pixel[i]))<<std::endl;
-	}
-	output.close();
-
 	input->cd();
 	break1_eff.Write();
+	break1_eff_ring03.Write();
+	break1_eff_ring47.Write();
+	break1_eff_ring811.Write();
+	break1_eff_ring1215.Write();
+	avg_theta_cm_ring03.Write();
+	avg_theta_cm_ring47.Write();
+	avg_theta_cm_ring811.Write();
+	avg_theta_cm_ring1215.Write();
 	break2_eff.Write();
 	input->Close();
 
@@ -285,6 +346,15 @@ void SabreEfficiency::Run3Step(const char* file) {
 	int count = 0;
 	int npercent = 0;
 
+	bool break1_flag, break3_flag, break4_flag;
+	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%
 			npercent++;
@@ -292,14 +362,26 @@ void SabreEfficiency::Run3Step(const char* file) {
 			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
 		}
 
+		break1_flag = false;
+		break3_flag = false;
+		break4_flag = false;
+
 		tree->GetEntry(i);
 
 		if(break1->KE > ENERGY_THRESHOLD) {
-			for(auto& det : detectors) {
+			for(int j=0; j<5; j++) {
-				if(det.GetTrajectoryCoordinates(break1->theta, break1->phi).GetX() != 0) {
+				auto& det = detectors[j];
+				auto chan = det.GetTrajectoryRingWedge(break1->theta, break1->phi);
+				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;
 				}
 			}
@@ -307,38 +389,88 @@ void SabreEfficiency::Run3Step(const char* file) {
 
 
 		if(break3->KE > ENERGY_THRESHOLD) {
-			for(auto& det : detectors) {
+			for(int j=0; j<5; j++) {
-				if(det.GetTrajectoryCoordinates(break3->theta, break3->phi).GetX() != 0) {
+				auto& det = detectors[j];
+				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) {
-			for(auto& det : detectors) {
+			for(int j=0; j<5; j++) {
-				if(det.GetTrajectoryCoordinates(break4->theta, break4->phi).GetX() != 0) {
+				auto& det = detectors[j];
+				auto chan = det.GetTrajectoryRingWedge(break4->theta, break4->phi);
+				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) {
+			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);
+			b3b4_kes->Fill(b3_kes[b3_kes.size()-1], b4_kes[b4_kes.size()-1]);
+			b3b4_phis->Fill(b3_phis[b3_phis.size()-1]/DEG2RAD, b4_phis[b4_phis.size()-1]/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++;
+			b3b4_thetas->Fill(b3_thetas[b3_thetas.size()-1]/DEG2RAD, b4_thetas[b4_thetas.size()-1]/DEG2RAD);
+			b3b4_kes->Fill(b3_kes[b3_kes.size()-1], b4_kes[b4_kes.size()-1]);
+			b3b4_phis->Fill(b3_phis[b3_phis.size()-1]/DEG2RAD, b4_phis[b4_phis.size()-1]/DEG2RAD);
+		}
 	}
 
 	double b1eff = ((double) b1_thetas.size())/nevents;
 	double b3eff = ((double) b3_thetas.size())/nevents;
 	double b4eff = ((double) b4_thetas.size())/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> break3_eff("Light Final Breakup Efficiency", b3eff);
 	TParameter<double> break4_eff("Heavy Final Breakup 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();
+	b3b4_thetas->Write();
+	b3b4_phis->Write();
+	b3b4_kes->Write();
 	input->Close();
 }

src/ThreeStepSystem.cpp

@@ -80,9 +80,9 @@ void ThreeStepSystem::RunSystem() {
 	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 decay1Theta = acos(generator->Uniform(-1, 1));
+	double decay1Theta = GetDecayTheta(L1);
 	double decay1Phi = generator->Uniform(0, M_PI*2.0);
-	double decay2Theta = acos(generator->Uniform(-1,1));
+	double decay2Theta = GetDecayTheta(L2);
 	double decay2Phi = generator->Uniform(0, M_PI*2.0);
 	double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
 

src/TwoStepSystem.cpp

@@ -71,7 +71,7 @@ void TwoStepSystem::RunSystem() {
 	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 decay1Theta = acos(generator->Uniform(-1, 1));
+	double decay1Theta = GetDecayTheta(L1);
 	double decay1Phi = generator->Uniform(0, M_PI*2.0);
 	double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
 

src/main.cpp

@@ -1,4 +1,5 @@
 #include <iostream>
+#include <string>
 #include "Kinematics.h"
 #include "SabreEfficiency.h"
 #include "SabreDetector.h"
@@ -22,7 +23,9 @@ int main(int argc, char** argv) {
 		return 1;
 	}
 	SabreEfficiency sabre;
+	std::string mapfile = "./etc/DeadChannels.txt";
 	sabre.SetReactionType(calculator.GetReactionType());
+	sabre.SetDeadChannelMap(mapfile);
 	sabre.CalculateEfficiency(calculator.GetOutputName());
 
 	/*std::vector<SabreDetector> detectors;