Added a nudge parameter

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.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);
+		SFPAnalyzer analyzer(m_params);
 	
 		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.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);
+		SFPAnalyzer analyzer(m_params);
 	
 		std::vector<TParameter<Double_t>> parvec;
 		parvec.reserve(9);

src/evb/EVBApp.cpp

@@ -72,6 +72,7 @@ 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>();
 	
@@ -116,6 +117,7 @@ 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,6 +28,8 @@ 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,68 +39,63 @@ namespace EventBuilder {
 
 	//requires (Z,A) for T, P, and E, as well as energy of P,
 	// spectrograph angle of interest, and field value
-	double Delta_Z(int ZT, int AT, int ZP, int AP, int ZE, int AE,
-		       double EP, double angle, double B) 
+	double DeltaZ(int ZT, int AT, int ZP, int AP, int ZE, int AE,
+		       double EP, double angle, double B, double nudge) 
 	{
 	
 		/* CONSTANTS */
-		const double UTOMEV = 931.4940954; //MeV per u;
-		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
+		static constexpr double s_speedOfLight = 2.9979E8; //m/s
+		static constexpr double s_qbrho2p = s_speedOfLight * 1.0e-9; //Converts QBrho -> p (kG cm -> MeV)
 	
 		/* SESPS-SPECIFIC */
-		const double DISP = 1.96; //dispersion (x/rho)
-		const double MAG = 0.39; //magnification in x
-		const double DEGTORAD = M_PI/180.;
+		static constexpr double s_dispersion = 1.96; //dispersion (x/rho)
+		static constexpr double s_magnification = 0.39; //magnification in x
+		static constexpr double s_deg2rad = M_PI/180.;
+		static constexpr double s_centralTrajAngle = 45.0 * s_deg2rad; //Central trajectory angle to detector plane
 	
-		int ZR = ZT + ZP - ZE, AR = AT + AP - AE;
-		double EE=0; //ejectile energy
+		int ZR = ZT + ZP - ZE;
+		int AR = AT + AP - AE;
+		double ejectileEnergy = 0.0; //ejectile energy
 	
-		double MT=0, MP=0, ME=0, MR=0; //masses (MeV)
+		angle *= s_deg2rad;
 	
-		B /= 10000; //convert to tesla
-		angle *= DEGTORAD;
+		//reactant masses
+		auto masses = MassLookup::GetInstance();
+		double mt = masses.FindMass(ZT, AT);
+		double mp = masses.FindMass(ZP, AP);
+		double me = masses.FindMass(ZE, AE);
+		double mr = masses.FindMass(ZR, AR);
 		
-		MT = MASS.FindMass(ZT, AT);
-		MP = MASS.FindMass(ZP, AP);
-		ME = MASS.FindMass(ZE, AE);
-		MR = MASS.FindMass(ZR, AR);
-		
-		if (MT*MP*ME*MR == 0) 
+		if (
+			MassLookup::IsInvalidMass(mt) ||
+			MassLookup::IsInvalidMass(mp) ||
+			MassLookup::IsInvalidMass(me) ||
+			MassLookup::IsInvalidMass(mr)
+		) 
 		{
 			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 = sqrt(MP*ME*EP)/(ME + MR)*cos(angle);
-		double term2 = (EP*(MR - MP) + MR*Q)/(ME + MR);
+		double term1 = std::sqrt(mp*me*EP)/(me + mr)*std::cos(angle);
+		double term2 = (EP*(mr - mp) + mr*Q)/(me + mr);
 	
-		EE = term1 + sqrt(term1*term1 + term2);
-		EE *= EE;
+		ejectileEnergy = term1 + std::sqrt(term1*term1 + term2);
+		ejectileEnergy *= ejectileEnergy;
 	
 		//momentum
-		double PE = sqrt(EE*(EE+2*ME));
-	
+		double ejectileP = std::sqrt(ejectileEnergy * (ejectileEnergy + 2.0*me));
 		//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 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));
 
-		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
+		return -1.0*rho*s_dispersion*s_magnification*K * std::cos(s_centralTrajAngle) + nudge; //delta-Z in cm
 	
 	}
 	
-	double Wire_Dist() {return 4.28625;} //cm
-
 }

src/evb/FP_kinematics.h

@@ -27,6 +27,8 @@
 
   KGH -- Jul19
 
+  B now in kG -- GWM 2023
+
 */
 
 #ifndef FP_KINEMATICS
@@ -36,10 +38,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 Delta_Z(int ZT, int AT, int ZP, int AP, int ZE, int AE,
-     	       double EP, double angle, double B); 
+     double DeltaZ(int ZT, int AT, int ZP, int AP, int ZE, int AE,
+     	       double EP, double angle, double B, double nudge = 0); 
      
-     double Wire_Dist();
+     static constexpr double WireDist() { return 4.28625; } //cm
 
 }
 

src/evb/MassLookup.cpp

@@ -12,6 +12,9 @@ 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/
@@ -28,18 +31,19 @@ namespace EventBuilder {
 		std::ifstream massfile(filepath);
 		if(massfile.is_open()) 
 		{
-			int Z,A;
-			std::string junk, element, key;
-			double atomicMassBig, atomicMassSmall, isotopicMass;
+			std::string junk, element;
+			NuclearData data;
+			double atomicMassBig, atomicMassSmall;
+			uint32_t id;
 			std::getline(massfile,junk);
 			std::getline(massfile,junk);
 			while(massfile>>junk) 
 			{
-				massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall;
-				isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev;
-				key = "("+std::to_string(Z)+","+A+")";
-				massTable[key] = isotopicMass;
-				elementTable[Z] = element;
+				massfile >> data.z >> data.a >> element >> atomicMassBig >> atomicMassSmall;
+				data.isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - data.z*s_electronMass)*s_u2MeV;
+				data.isotopicSymbol = fmt::format("{}{}", data.a, element);
+				id = GetIsotopeID(data.z, data.a);
+				m_dataMap[id] = data;
 			}
 		} 
 		else
@@ -49,28 +53,28 @@ namespace EventBuilder {
 	MassLookup::~MassLookup() {}
 	
 	//Returns nuclear mass in MeV
-	double MassLookup::FindMass(int Z, int A) 
+	double MassLookup::FindMass(uint32_t Z, uint32_t A) 
 	{
-		std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")";
-		auto data = massTable.find(key);
-		if(data == massTable.end()) 
+		uint32_t id = GetIsotopeID(Z, A);
+		auto data = m_dataMap.find(id);
+		if(data == m_dataMap.end()) 
 		{
-			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindMass; returning zero.",Z,A);
-			return 0;
+			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindMass; returning invalid.",Z,A);
+			return s_invalidMass;
 		}
-		return data->second;
+		return data->second.isotopicMass;
 	}
 	
 	//returns element symbol
-	std::string MassLookup::FindSymbol(int Z, int A) 
+	std::string MassLookup::FindSymbol(uint32_t Z, uint32_t A) 
 	{
-		auto data = elementTable.find(Z);
-		if(data == elementTable.end()) 
+		uint32_t id = GetIsotopeID(Z, A);
+		auto data = m_dataMap.find(id);
+		if(data == m_dataMap.end()) 
 		{
-			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindSymbol; returning empty string.",Z,A);
-			return "";
+			EVB_WARN("Invalid nucleus (Z,A) ({0},{1}) at MassLookup::FindSymbol; returning Invalid.",Z,A);
+			return s_invalidSymbol;
 		}
-		std::string fullsymbol = std::to_string(A) + data->second;
-		return fullsymbol;
+		return data->second.isotopicSymbol;
 	}
 }

src/evb/MassLookup.h

@@ -7,33 +7,50 @@ 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 {
 
-	class MassLookup 
+	static constexpr uint32_t GetIsotopeID(uint32_t z, uint32_t a)
 	{
+		return z >= a ? z*z + z + a : a*a + z;
+	}
 
-	public:
-		MassLookup();
-		~MassLookup();
-		double FindMass(int Z, int A);
-		std::string FindSymbol(int Z, int A);
-	
-	private:
-		std::unordered_map<std::string, double> massTable;
-		std::unordered_map<int, std::string> elementTable;
-	
-		//constants
-		static constexpr double u_to_mev = 931.4940954;
-		static constexpr double electron_mass = 0.000548579909;
-		  
+	struct NuclearData
+	{
+		std::string isotopicSymbol = "";
+		double isotopicMass = 0.0;
+		uint32_t z = 0;
+		uint32_t a = 0;
 	};
 
-	//static instance for use throught program
-	static MassLookup MASS;
+	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; }
+		static bool IsInvalidMass(double mass) { return mass == s_invalidMass; }
+		static bool IsInvalidSymbol(const std::string& sym) { return sym == s_invalidSymbol; }
+	
+	private:
+		MassLookup();
+		std::unordered_map<uint32_t, NuclearData> m_dataMap;
+	
+		static MassLookup* s_instance;
+
+		//constants
+		static constexpr double s_u2MeV = 931.4940954; // U to MeV
+		static constexpr double s_electronMass = 0.000548579909; // electron mass in u
+		static constexpr double s_invalidMass = -1.0; //Invalid mass in MeV
+		static constexpr char s_invalidSymbol[] = "Invalid"; //Invalid isotope symbol
+	};
 	
 }
 #endif

src/evb/SFPAnalyzer.cpp

@@ -14,10 +14,11 @@
 namespace EventBuilder {
 
 	/*Constructor takes in kinematic parameters for generating focal plane weights*/
-	SFPAnalyzer::SFPAnalyzer(int zt, int at, int zp, int ap, int ze, int ae, double ep,
-								double angle, double b) 
+	SFPAnalyzer::SFPAnalyzer(const EVBParameters& params) 
 	{
-		zfp = Delta_Z(zt, at, zp, ap, ze, ae, ep, angle, b*1000.0); //Convert kG to G
+		zfp = DeltaZ(params.ZT, params.AT, params.ZP, params.AP, params.ZE, params.AE, params.beamEnergy, 
+					params.spsAngle, params.BField, params.nudge);
+		EVB_INFO("Focal Plane Z-Offset (nudge + kinematics): {}", zfp);
 		event_address = new CoincEvent();
 		rootObj = new THashTable();
 		GetWeights();
@@ -42,7 +43,7 @@ namespace EventBuilder {
 	 */
 	void SFPAnalyzer::GetWeights() 
 	{
-		w1 = (Wire_Dist()/2.0-zfp)/Wire_Dist();
+		w1 = (WireDist()/2.0-zfp)/WireDist();
 		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(int zt, int at, int zp, int ap, int ze, int ae, double ep, double angle,
-		            double b);
+		SFPAnalyzer(const EVBParameters& params);
 		~SFPAnalyzer();
 		ProcessedEvent GetProcessedEvent(CoincEvent& event);
 		inline void ClearHashTable() { rootObj->Clear(); }

src/guidict/EVBMainFrame.cpp

@@ -145,9 +145,15 @@ 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);
@@ -399,6 +405,7 @@ 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,6 +67,7 @@ public:
 		TypeBox,
 		RMin,
 		RMax,
+		Nudge,
 		M_Load_Config,
 		M_Save_Config,
 		M_Exit
@@ -85,6 +86,7 @@ private:
 	TGNumberEntryField *fBField, *fBKEField, *fThetaField;
 	TGNumberEntryField *fSlowWindowField, *fFastICField, *fFastSABREField;
 	TGNumberEntryField *fRMinField, *fRMaxField;
+	TGNumberEntryField *fNudgeField;
 
 	TGHProgressBar* fProgressBar;