src/evb/CompassRun.cpp

@@ -369,7 +369,7 @@ namespace EventBuilder {
 		startIndex = 0;
 		CoincEvent this_event;
 		SlowSort coincidizer(m_params.slowCoincidenceWindow, m_params.channelMapFile);
-		SFPAnalyzer analyzer(m_params);
+		SFPAnalyzer analyzer(m_params.ZT, m_params.AT, m_params.ZP, m_params.AP, m_params.ZE, m_params.AE, m_params.beamEnergy, m_params.spsAngle, m_params.BField);
 	
 		std::vector<TParameter<Double_t>> parvec;
 		parvec.reserve(9);
@@ -459,7 +459,7 @@ namespace EventBuilder {
 		std::vector<CoincEvent> fast_events;
 		SlowSort coincidizer(m_params.slowCoincidenceWindow, m_params.channelMapFile);
 		FastSort speedyCoincidizer(m_params.fastCoincidenceWindowSABRE, m_params.fastCoincidenceWindowIonCh);
-		SFPAnalyzer analyzer(m_params);
+		SFPAnalyzer analyzer(m_params.ZT, m_params.AT, m_params.ZP, m_params.AP, m_params.ZE, m_params.AE, m_params.beamEnergy, m_params.spsAngle, m_params.BField);
 	
 		std::vector<TParameter<Double_t>> parvec;
 		parvec.reserve(9);

src/evb/EVBApp.cpp

@@ -72,7 +72,6 @@ namespace EventBuilder {
 		m_params.BField = data["BField(kG)"].as<double>();
 		m_params.beamEnergy = data["BeamEnergy(MeV)"].as<double>();
 		m_params.spsAngle = data["SPSAngle(deg)"].as<double>();
-		m_params.nudge = data["Nudge(cm)"].as<double>();
 		m_params.runMin = data["MinRun"].as<int>();
 		m_params.runMax = data["MaxRun"].as<int>();
 	
@@ -117,7 +116,6 @@ namespace EventBuilder {
 		yamlStream << YAML::Key << "BField(kG)" << YAML::Value << m_params.BField;
 		yamlStream << YAML::Key << "BeamEnergy(MeV)" << YAML::Value << m_params.beamEnergy;
 		yamlStream << YAML::Key << "SPSAngle(deg)" << YAML::Value << m_params.spsAngle;
-		yamlStream << YAML::Key << "Nudge(cm)" << YAML::Value << m_params.nudge;
 		yamlStream << YAML::Key << "MinRun" << YAML::Value << m_params.runMin;
 		yamlStream << YAML::Key << "MaxRun" << YAML::Value << m_params.runMax;
 		yamlStream << YAML::EndMap;

src/evb/EVBParameters.h

@@ -28,8 +28,6 @@ namespace EventBuilder {
 		double BField = 7.8; //kG
 		double spsAngle = 15.0; //degrees
 		double beamEnergy = 16.0; //MeV
-
-		double nudge = 0.0; //User FP offset factor, cm
 	};
 }
 

src/evb/FP_kinematics.cpp

@@ -39,63 +39,68 @@ namespace EventBuilder {
 
 	//requires (Z,A) for T, P, and E, as well as energy of P,
 	// spectrograph angle of interest, and field value
-	double DeltaZ(int ZT, int AT, int ZP, int AP, int ZE, int AE,
+	double Delta_Z(int ZT, int AT, int ZP, int AP, int ZE, int AE,
-		       double EP, double angle, double B, double nudge) 
+		       double EP, double angle, double B) 
 	{
 	
 		/* CONSTANTS */
-		static constexpr double s_speedOfLight = 2.9979E8; //m/s
+		const double UTOMEV = 931.4940954; //MeV per u;
-		static constexpr double s_qbrho2p = s_speedOfLight * 1.0e-9; //Converts QBrho -> p (kG cm -> MeV)
+		const double MEVTOJ = 1.60218E-13; //J per MeV
+		const double RESTMASS_ELECTRON = 0.000548579909; //amu
+		const double UNIT_CHARGE = 1.602E-19; //Coulombs
+		const double C = 2.9979E8; //m/s
 	
 		/* SESPS-SPECIFIC */
-		static constexpr double s_dispersion = 1.96; //dispersion (x/rho)
+		const double DISP = 1.96; //dispersion (x/rho)
-		static constexpr double s_magnification = 0.39; //magnification in x
+		const double MAG = 0.39; //magnification in x
-		static constexpr double s_deg2rad = M_PI/180.;
+		const double DEGTORAD = M_PI/180.;
-		static constexpr double s_centralTrajAngle = 45.0 * s_deg2rad; //Central trajectory angle to detector plane
 	
-		int ZR = ZT + ZP - ZE;
+		int ZR = ZT + ZP - ZE, AR = AT + AP - AE;
-		int AR = AT + AP - AE;
+		double EE=0; //ejectile energy
-		double ejectileEnergy = 0.0; //ejectile energy
 	
-		angle *= s_deg2rad;
+		double MT=0, MP=0, ME=0, MR=0; //masses (MeV)
 	
-		//reactant masses
+		B /= 10000; //convert to tesla
-		auto masses = MassLookup::GetInstance();
+		angle *= DEGTORAD;
-		double mt = masses.FindMass(ZT, AT);
-		double mp = masses.FindMass(ZP, AP);
-		double me = masses.FindMass(ZE, AE);
-		double mr = masses.FindMass(ZR, AR);
 	
-		if (
+		MT = MASS.FindMass(ZT, AT);
-			MassLookup::IsInvalidMass(mt) ||
+		MP = MASS.FindMass(ZP, AP);
-			MassLookup::IsInvalidMass(mp) ||
+		ME = MASS.FindMass(ZE, AE);
-			MassLookup::IsInvalidMass(me) ||
+		MR = MASS.FindMass(ZR, AR);
-			MassLookup::IsInvalidMass(mr)
+		
-		) 
+		if (MT*MP*ME*MR == 0) 
 		{
 			EVB_WARN("Illegal mass at FP_kinematics::Delta_Z! Returning offset of 0.");
 			return 0;
 		}
 	
-		double Q = mt + mp - me - mr; //Q-value
+		double Q = MT + MP - ME - MR; //Q-value
 		
 		//kinematics a la Iliadis p.590
-		double term1 = std::sqrt(mp*me*EP)/(me + mr)*std::cos(angle);
+		double term1 = sqrt(MP*ME*EP)/(ME + MR)*cos(angle);
-		double term2 = (EP*(mr - mp) + mr*Q)/(me + mr);
+		double term2 = (EP*(MR - MP) + MR*Q)/(ME + MR);
 	
-		ejectileEnergy = term1 + std::sqrt(term1*term1 + term2);
+		EE = term1 + sqrt(term1*term1 + term2);
-		ejectileEnergy *= ejectileEnergy;
+		EE *= EE;
 	
 		//momentum
-		double ejectileP = std::sqrt(ejectileEnergy * (ejectileEnergy + 2.0*me));
+		double PE = sqrt(EE*(EE+2*ME));
-		//calculate rho from B a la B*rho = (proj. momentum)/(proj. charge)
-		double rho = ejectileP / (ZE * B * s_qbrho2p ); //in cm
-		//kinematic factor K
-		double K = (std::sqrt(mp*me*EP / ejectileEnergy) * std::sin(angle)) /
-				   (me + mr - std::sqrt(mp*me*EP / ejectileEnergy) * std::cos(angle));
 	
-		return -1.0*rho*s_dispersion*s_magnification*K * std::cos(s_centralTrajAngle) + nudge; //delta-Z in cm
+		//calculate rho from B a la B*rho = (proj. momentum)/(proj. charge)
+		double rho = (PE*MEVTOJ)/(ZE*UNIT_CHARGE*C*B)*100; //in cm
+	
+		double K;
+	
+		K  = sqrt(MP*ME*EP/EE);
+		K *= sin(angle);
+	
+		double denom = ME + MR - sqrt(MP*ME*EP/EE)*cos(angle);
+	
+		K /= denom;
+		return -1*rho*DISP*MAG*K; //delta-Z in cm
 	
 	}
 	
+	double Wire_Dist() {return 4.28625;} //cm
+
 }

src/evb/FP_kinematics.h

@@ -27,8 +27,6 @@
 
   KGH -- Jul19
 
-  B now in kG -- GWM 2023
-
 */
 
 #ifndef FP_KINEMATICS
@@ -38,10 +36,10 @@ namespace EventBuilder {
 
      //requires (Z,A) for T, P, and E, as well as energy of P,
      // spectrograph angle of interest, and field value
-     double DeltaZ(int ZT, int AT, int ZP, int AP, int ZE, int AE,
+     double Delta_Z(int ZT, int AT, int ZP, int AP, int ZE, int AE,
-     	       double EP, double angle, double B, double nudge = 0); 
+     	       double EP, double angle, double B); 
      
-     static constexpr double WireDist() { return 4.28625; } //cm
+     double Wire_Dist();
 
 }
 

src/evb/MassLookup.cpp

@@ -12,9 +12,6 @@ Written by G.W. McCann Aug. 2020
 
 namespace EventBuilder {
 	
-	//instantiate
-	MassLookup* MassLookup::s_instance = new MassLookup();
-	
 	/*
 	  Read in AMDC mass file, preformated to remove excess info. Here assumes that by default
 	  the file is in a local directory etc/
@@ -31,19 +28,18 @@ namespace EventBuilder {
 		std::ifstream massfile(filepath);
 		if(massfile.is_open()) 
 		{
-			std::string junk, element;
+			int Z,A;
-			NuclearData data;
+			std::string junk, element, key;
-			double atomicMassBig, atomicMassSmall;
+			double atomicMassBig, atomicMassSmall, isotopicMass;
-			uint32_t id;
 			std::getline(massfile,junk);
 			std::getline(massfile,junk);
 			while(massfile>>junk) 
 			{
-				massfile >> data.z >> data.a >> element >> atomicMassBig >> atomicMassSmall;
+				massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall;
-				data.isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - data.z*s_electronMass)*s_u2MeV;
+				isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev;
-				data.isotopicSymbol = fmt::format("{}{}", data.a, element);
+				key = "("+std::to_string(Z)+","+A+")";
-				id = GetIsotopeID(data.z, data.a);
+				massTable[key] = isotopicMass;
-				m_dataMap[id] = data;
+				elementTable[Z] = element;
 			}
 		} 
 		else
@@ -53,28 +49,28 @@ namespace EventBuilder {
 	MassLookup::~MassLookup() {}
 	
 	//Returns nuclear mass in MeV
-	double MassLookup::FindMass(uint32_t Z, uint32_t A) 
+	double MassLookup::FindMass(int Z, int A) 
 	{
-		uint32_t id = GetIsotopeID(Z, A);
+		std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")";
-		auto data = m_dataMap.find(id);
+		auto data = massTable.find(key);
-		if(data == m_dataMap.end()) 
+		if(data == massTable.end()) 
 		{
-			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindMass; returning invalid.",Z,A);
+			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindMass; returning zero.",Z,A);
-			return s_invalidMass;
+			return 0;
 		}
-		return data->second.isotopicMass;
+		return data->second;
 	}
 	
 	//returns element symbol
-	std::string MassLookup::FindSymbol(uint32_t Z, uint32_t A) 
+	std::string MassLookup::FindSymbol(int Z, int A) 
 	{
-		uint32_t id = GetIsotopeID(Z, A);
+		auto data = elementTable.find(Z);
-		auto data = m_dataMap.find(id);
+		if(data == elementTable.end()) 
-		if(data == m_dataMap.end()) 
 		{
-			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindSymbol; returning Invalid.",Z,A);
+			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindSymbol; returning empty string.",Z,A);
-			return s_invalidSymbol;
+			return "";
 		}
-		return data->second.isotopicSymbol;
+		std::string fullsymbol = std::to_string(A) + data->second;
+		return fullsymbol;
 	}
 }

src/evb/MassLookup.h

@@ -7,50 +7,33 @@ of this map (MASS) for use throughout code it is included into.
 
 Written by G.W. McCann Aug. 2020
 
-Updated for better code practices and real singleton usage GWM 2023
-
 */
 #ifndef MASS_LOOKUP_H
 #define MASS_LOOKUP_H
 
 namespace EventBuilder {
 
-	static constexpr uint32_t GetIsotopeID(uint32_t z, uint32_t a)
-	{
-		return z >= a ? z*z + z + a : a*a + z;
-	}
-
-	struct NuclearData
-	{
-		std::string isotopicSymbol = "";
-		double isotopicMass = 0.0;
-		uint32_t z = 0;
-		uint32_t a = 0;
-	};
-
 	class MassLookup 
 	{
-	public:
-		~MassLookup();
-		double FindMass(uint32_t Z, uint32_t A);
-		std::string FindSymbol(uint32_t Z, uint32_t A);
 	
-		static MassLookup& GetInstance() { return *s_instance; }
+	public:
-		static bool IsInvalidMass(double mass) { return mass == s_invalidMass; }
+		MassLookup();
-		static bool IsInvalidSymbol(const std::string& sym) { return sym == s_invalidSymbol; }
+		~MassLookup();
+		double FindMass(int Z, int A);
+		std::string FindSymbol(int Z, int A);
 	
 	private:
-		MassLookup();
+		std::unordered_map<std::string, double> massTable;
-		std::unordered_map<uint32_t, NuclearData> m_dataMap;
+		std::unordered_map<int, std::string> elementTable;
-	
-		static MassLookup* s_instance;
 	
 		//constants
-		static constexpr double s_u2MeV = 931.4940954; // U to MeV
+		static constexpr double u_to_mev = 931.4940954;
-		static constexpr double s_electronMass = 0.000548579909; // electron mass in u
+		static constexpr double electron_mass = 0.000548579909;
-		static constexpr double s_invalidMass = -1.0; //Invalid mass in MeV
+		  
-		static constexpr char s_invalidSymbol[] = "Invalid"; //Invalid isotope symbol
 	};
 	
+	//static instance for use throught program
+	static MassLookup MASS;
+
 }
 #endif

src/evb/SFPAnalyzer.cpp

@@ -14,11 +14,10 @@
 namespace EventBuilder {
 
 	/*Constructor takes in kinematic parameters for generating focal plane weights*/
-	SFPAnalyzer::SFPAnalyzer(const EVBParameters& params) 
+	SFPAnalyzer::SFPAnalyzer(int zt, int at, int zp, int ap, int ze, int ae, double ep,
+								double angle, double b) 
 	{
-		zfp = DeltaZ(params.ZT, params.AT, params.ZP, params.AP, params.ZE, params.AE, params.beamEnergy, 
+		zfp = Delta_Z(zt, at, zp, ap, ze, ae, ep, angle, b*1000.0); //Convert kG to G
-					params.spsAngle, params.BField, params.nudge);
-		EVB_INFO("Focal Plane Z-Offset (nudge + kinematics): {}", zfp);
 		event_address = new CoincEvent();
 		rootObj = new THashTable();
 		GetWeights();
@@ -43,7 +42,7 @@ namespace EventBuilder {
 	 */
 	void SFPAnalyzer::GetWeights() 
 	{
-		w1 = (WireDist()/2.0-zfp)/WireDist();
+		w1 = (Wire_Dist()/2.0-zfp)/Wire_Dist();
 		w2 = 1.0-w1;
 		EVB_INFO("Calculated X-Avg weights of w1={0} and w2={1}",w1,w2);
 	}

src/evb/SFPAnalyzer.h

@@ -12,14 +12,14 @@
 
 #include "DataStructs.h"
 #include "FP_kinematics.h"
-#include "EVBParameters.h"
 
 namespace EventBuilder {
 
 	class SFPAnalyzer
 	{
 	public:
-		SFPAnalyzer(const EVBParameters& params);
+		SFPAnalyzer(int zt, int at, int zp, int ap, int ze, int ae, double ep, double angle,
+		            double b);
 		~SFPAnalyzer();
 		ProcessedEvent GetProcessedEvent(CoincEvent& event);
 		inline void ClearHashTable() { rootObj->Clear(); }

src/guidict/EVBMainFrame.cpp

@@ -145,15 +145,9 @@ EVBMainFrame::EVBMainFrame(const TGWindow* p, UInt_t w, UInt_t h) :
 	fThetaField = new TGNumberEntryField(thetaFrame, Theta, 0, TGNumberEntry::kNESRealFour, TGNumberEntry::kNEANonNegative);
 	thetaFrame->AddFrame(thetalabel, lhints);
 	thetaFrame->AddFrame(fThetaField, fhints);
-	TGHorizontalFrame* nudgeFrame = new TGHorizontalFrame(extraFrame, w*0.175, h*0.15);
-	TGLabel* nudgeLabel = new TGLabel(nudgeFrame, "Nudge (cm):");
-	fNudgeField = new TGNumberEntryField(nudgeFrame, Nudge, 0, TGNumberEntry::kNESRealFour);
-	nudgeFrame->AddFrame(nudgeLabel, lhints);
-	nudgeFrame->AddFrame(fNudgeField, fhints);
 	extraFrame->AddFrame(beamFrame, fhints);
 	extraFrame->AddFrame(bfFrame, fhints);
 	extraFrame->AddFrame(thetaFrame, fhints);
-	extraFrame->AddFrame(nudgeFrame, fhints);
 
 	reactionFrame->AddFrame(targFrame, fhints);
 	reactionFrame->AddFrame(projFrame, fhints);
@@ -405,7 +399,6 @@ bool EVBMainFrame::SetParameters()
 	m_parameters.BField = fBField->GetNumber();
 	m_parameters.beamEnergy = fBKEField->GetNumber();
 	m_parameters.spsAngle = fThetaField->GetNumber();
-	m_parameters.nudge = fNudgeField->GetNumber();
 
 	m_builder.SetParameters(m_parameters);
 	return true;

src/guidict/EVBMainFrame.h

@@ -67,7 +67,6 @@ public:
 		TypeBox,
 		RMin,
 		RMax,
-		Nudge,
 		M_Load_Config,
 		M_Save_Config,
 		M_Exit
@@ -86,7 +85,6 @@ private:
 	TGNumberEntryField *fBField, *fBKEField, *fThetaField;
 	TGNumberEntryField *fSlowWindowField, *fFastICField, *fFastSABREField;
 	TGNumberEntryField *fRMinField, *fRMaxField;
-	TGNumberEntryField *fNudgeField;
 
 	TGHProgressBar* fProgressBar;