fixed merge conflicts

include/EnergyLoss.h

@@ -46,6 +46,7 @@ class EnergyLoss {
 
     //constants for calculations
     static constexpr double MAX_FRACTIONAL_STEP = 0.001;
+    static constexpr double MAX_DEPTH = 50;
     static constexpr double MAX_H_E_PER_U = 100000.0;
     static constexpr double AVOGADRO = 0.60221367; //N_A times 10^(-24) for converting
     static constexpr double MEV2U = 1.0/931.4940954;

include/Nucleus.h

@@ -1,3 +1,11 @@
+/*
+	Nucleus.h
+	Nucleus is a derived class of Vec4. A nucleus is the kinematics is essentially a 4 vector with the
+	additional properties of the number of total nucleons (A), the number of protons (Z), a ground state mass,
+	an exctitation energy, and an isotopic symbol.
+
+	--GWM Jan 2021
+*/
 #ifndef NUCLEUS_H
 #define NUCLEUS_H
 
@@ -25,7 +33,8 @@ public:
 		return *this;
 	};
 
-	inline Nucleus operator+(const Nucleus& daughter) { 
+	//Conservation of nucleons and momentum
+	inline Nucleus operator+(const Nucleus& daughter) {
 		return Nucleus(GetZ()+daughter.GetZ(), GetA()+daughter.GetA(), GetPx()+daughter.GetPx(), GetPy()+daughter.GetPy(), GetPz()+daughter.GetPz(), GetE()+daughter.GetE()); 
 	};
 	inline Nucleus operator-(const Nucleus& daughter) {

include/Reaction.h

@@ -1,3 +1,11 @@
+/*
+	Reaction.h
+	Reaction is a class which implements either a decay or scattering reaction. As such it requires either
+	3 (decay) or 4 (scattering) nuclei to perform any calcualtions. I also links together the target, which provides
+	energy loss calculations, with the kinematics. Note that Reaction does not own the LayeredTarget.
+
+	--GWM Jan. 2021
+*/
 #ifndef REACTION_H
 #define REACTION_H
 
@@ -15,6 +23,8 @@ public:
 	void SetNuclei(int zt, int at, int zp, int ap, int ze, int ae);
 	void SetNuclei(const Nucleus* nucs);
 	void SetBeamKE(double bke);
+
+	/*Setters and getters*/
 	inline void SetLayeredTarget(LayeredTarget* targ) { target = targ; };
 	inline void SetPolarRxnAngle(double theta) { m_theta = theta; };
 	inline void SetAzimRxnAngle(double phi) { m_phi = phi; };

include/ReactionSystem.h

@@ -1,3 +1,11 @@
+/*
+	ReactionSystem.h
+	ReactionSystem is the base class from which all kinematics calculations should inherit. It contains all members and functions
+	to perform a single step reaction/decay, which is the basic building block of all subsequent types of reactions in MASK. More
+	complicated systems (see TwoStepSystem and ThreeStepSystem) add further data members and override the virtual functions.
+
+	--GWM Jan. 2021
+*/
 #ifndef REACTIONSYSTEM_H
 #define REACTIONSYSTEM_H
 
@@ -12,12 +20,16 @@ public:
 	ReactionSystem();
 	ReactionSystem(std::vector<int>& z, std::vector<int>& a);
 	virtual ~ReactionSystem();
+
 	virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a);
 	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; };
 	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); };
+
 	virtual void RunSystem();
 
 	inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
@@ -32,6 +44,8 @@ protected:
 	
 	Reaction step1;
 	LayeredTarget target;
+
+	//Sampling information
 	std::pair<double, double> m_beamDist, m_theta1Range, m_exDist;
 	TRandom3* generator; //not owned by ReactionSystem
 

input.txt

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

input.txt.orig

@@ -0,0 +1,46 @@
+----------Data Information----------
+OutputFile: /data1/gwm17/test_dead.root
+SaveTree: yes
+SavePlots: yes
+----------Reaction Information----------
+ReactionType: 2
+Z	A (order is target, projectile, ejectile, break1, break3)
+<<<<<<< HEAD
+6	12
+2	3
+1	2
+1	1
+=======
+5	10
+2	3
+2	4
+1	2
+>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7
+----------Target Information----------
+Name: test_targ
+Layers: 1
+~Layer1
+<<<<<<< HEAD
+Thickness(ug/cm^2): 40
+=======
+Thickness(ug/cm^2): 10
+>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7
+Z	A	Stoich
+6	12	1
+0
+~
+<<<<<<< HEAD
+=======
+~Layer2
+Thickness(ug/cm^2): 80
+Z	A	Stoich
+5	10	1
+0
+~
+>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7
+----------Sampling Information----------
+NumberOfSamples: 1000000
+BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001
+EjectileThetaMin(deg): 3.0 EjectileThetaMax(deg): 3.0
+ResidualExMean(MeV): 2.364 ResidualExSigma(MeV): 0.0317
+--------------------------------------

src/EnergyLoss.cpp

@@ -55,6 +55,8 @@ double EnergyLoss::GetEnergyLoss(int zp, int ap, double e_initial, double thickn
   double e_step = GetTotalStoppingPower(e_final)*x_step/1000.0; //initial step in e
   double e_threshold = 0.05*e_initial;
 
+  int depth=0;
+
   
   if(thickness == 0.0) return 0;
   else if(e_initial == 0.0) return 0;
@@ -62,16 +64,20 @@ double EnergyLoss::GetEnergyLoss(int zp, int ap, double e_initial, double thickn
   bool go = true;
   while(go) {
     //If intial guess of step size is too large, shrink until in range
-    if(e_step/e_final > MAX_FRACTIONAL_STEP /*&& e_step >= E_PRECISION_LIMIT*/) {
+    if(e_step/e_final > MAX_FRACTIONAL_STEP && depth < MAX_DEPTH) {
+      depth++;
       x_step *= 0.5;
       e_step = GetTotalStoppingPower(e_final)*x_step/1000.0;
     } else if((x_step + x_traversed) >= thickness) { //last chunk
       go = false;
       x_step = thickness - x_traversed; //get valid portion of last chunk
       e_final -= GetTotalStoppingPower(e_final)*x_step/1000.0;
+      if(depth > 20)std::cout<<"depth: "<<depth<<std::endl;
       if(e_final <= e_threshold) {
         return e_initial;
       }
+    } else if(depth == MAX_DEPTH) {
+      return e_initial;
     } else {
       x_traversed += x_step;
       e_step = GetTotalStoppingPower(e_final)*x_step/1000.0;

src/Nucleus.cpp

@@ -1,3 +1,11 @@
+/*
+	Nucleus.cpp
+	Nucleus is a derived class of Vec4. A nucleus is the kinematics is essentially a 4 vector with the
+	additional properties of the number of total nucleons (A), the number of protons (Z), a ground state mass,
+	an exctitation energy, and an isotopic symbol.
+
+	--GWM Jan 2021
+*/
 #include "Nucleus.h"
 #include "MassLookup.h"
 
@@ -13,7 +21,7 @@ Nucleus::Nucleus(int Z, int A) :
 {
 	m_gs_mass = MASS.FindMass(Z, A);
 	m_symbol = MASS.FindSymbol(Z, A);
-	SetVectorCartesian(0,0,0,m_gs_mass);
+	SetVectorCartesian(0,0,0,m_gs_mass); //by defualt a nucleus has mass given by the g.s.
 }
 
 Nucleus::Nucleus(int Z, int A, double px, double py, double pz, double E) :

src/Reaction.cpp

@@ -1,3 +1,11 @@
+/*
+	Reaction.cpp
+	Reaction is a class which implements either a decay or scattering reaction. As such it requires either
+	3 (decay) or 4 (scattering) nuclei to perform any calcualtions. I also links together the target, which provides
+	energy loss calculations, with the kinematics. Note that Reaction does not own the LayeredTarget.
+
+	--GWM Jan. 2021
+*/
 #include "Reaction.h"
 #include "KinematicsExceptions.h"
 
@@ -69,6 +77,7 @@ void Reaction::SetBeamKE(double bke) {
 	m_bke = bke - target->GetProjectileEnergyLoss(reactants[1].GetZ(), reactants[1].GetA(), bke, rxnLayer, 0);
 };
 
+//Methods given by Iliadis in Nuclear Physics of Stars, Appendix C
 void Reaction::CalculateReaction() {
 	//Target assumed at rest, with 0 excitation energy
 	reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass());
@@ -112,6 +121,7 @@ void Reaction::CalculateReaction() {
 
 }
 
+//Calculate in CM, where decay is isotropic
 void Reaction::CalculateDecay() {
 
 	double Q = reactants[0].GetInvMass() - reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass();

src/Target.cpp

@@ -40,28 +40,28 @@ bool Target::ContainsElement(int z, int a) {
 /*Calculates energy loss for travelling all the way through the target*/
 double Target::getEnergyLossTotal(int zp, int ap, double startEnergy, double theta) {
   if(theta == PI/2.) return startEnergy;
-  else if(theta > PI/2.) theta = PI - theta;
+  else if (theta > PI/2.) theta = PI - theta;
   return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta)));
 }
 
 /*Calculates energy loss for travelling halfway through the target*/
 double Target::getEnergyLossHalf(int zp, int ap, double startEnergy, double theta) {
   if(theta == PI/2.) return startEnergy;
-  else if(theta > PI/2.) theta = PI - theta;
+  else if (theta > PI/2.) theta = PI - theta;
   return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/(2.0*fabs(cos(theta))));
 }
 
 /*Calculates reverse energy loss for travelling all the way through the target*/
 double Target::getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta) {
   if(theta == PI/2.) return finalEnergy;
-  else if(theta > PI/2.) theta = PI - theta;
+  else if (theta > PI/2.) theta = PI - theta;
   return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta)));
 }
 
 /*Calculates reverse energy loss for travelling half way through the target*/
 double Target::getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double theta) {
   if(theta == PI/2.) return finalEnergy;
-  else if(theta > PI/2.) theta = PI - theta;
+  else if (theta > PI/2.) theta = PI - theta;
   return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta))));
 }
 

src/Target.cpp.orig

@@ -0,0 +1,104 @@
+/*
+
+Target.cpp
+A basic target unit for use in the SPANCRedux environment. A target
+is defined as a single compound with elements Z,A of a given stoichiometry 
+Holds an energy loss class
+
+Based on code by D.W. Visser written at Yale for the original SPANC
+
+Written by G.W. McCann Aug. 2020
+
+*/
+#include "Target.h"
+
+/*Targets must be of known thickness*/
+Target::Target(double thick) {
+  thickness = thick;
+}
+
+Target::~Target() {
+}
+
+/*Set target elements of given Z, A, S*/
+void Target::SetElements(std::vector<int>& z, std::vector<int>& a, std::vector<int>& stoich) {
+  Z = z;
+  A = a;
+  Stoich = stoich;
+  
+  eloss.SetTargetComponents(Z, A, Stoich);
+}
+
+/*Element verification*/
+bool Target::ContainsElement(int z, int a) {
+  for(unsigned int i=0; i<Z.size(); i++) {
+    if( z == Z[i] && a == A[i]) return true;
+  }
+  return false;
+}
+
+/*Calculates energy loss for travelling all the way through the target*/
+double Target::getEnergyLossTotal(int zp, int ap, double startEnergy, double theta) {
+  if(theta == PI/2.) return startEnergy;
+<<<<<<< HEAD
+  else if(theta > PI/2.) theta = PI - theta;
+=======
+  else if (theta > PI/2.)  theta = PI - theta;
+>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7
+  return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta)));
+}
+
+/*Calculates energy loss for travelling halfway through the target*/
+double Target::getEnergyLossHalf(int zp, int ap, double startEnergy, double theta) {
+  if(theta == PI/2.) return startEnergy;
+<<<<<<< HEAD
+  else if(theta > PI/2.) theta = PI - theta;
+=======
+  else if (theta > PI/2.)  theta = PI - theta;
+>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7
+  return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/(2.0*fabs(cos(theta))));
+}
+
+/*Calculates reverse energy loss for travelling all the way through the target*/
+double Target::getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta) {
+  if(theta == PI/2.) return finalEnergy;
+<<<<<<< HEAD
+  else if(theta > PI/2.) theta = PI - theta;
+=======
+  else if (theta > PI/2.)  theta = PI - theta;
+>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7
+  return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta)));
+}
+
+/*Calculates reverse energy loss for travelling half way through the target*/
+double Target::getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double theta) {
+  if(theta == PI/2.) return finalEnergy;
+<<<<<<< HEAD
+  else if(theta > PI/2.) theta = PI - theta;
+=======
+  else if (theta > PI/2.)  theta = PI - theta;
+>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7
+  return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta))));
+}
+
+/*Getter functions*/
+
+double& Target::GetThickness() {
+  return thickness;
+}
+
+int Target::GetNumberOfElements() {
+  return Z.size();
+}
+
+int Target::GetElementZ(int index) {
+  return Z[index];
+}
+
+int Target::GetElementA(int index) {
+  return A[index];
+}
+
+int Target::GetElementStoich(int index) {
+  return Stoich[index];
+}