diff --git a/src/evb/EVBParameters.h b/src/evb/EVBParameters.h index 551336c..8344a45 100644 --- a/src/evb/EVBParameters.h +++ b/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 }; } diff --git a/src/evb/FP_kinematics.cpp b/src/evb/FP_kinematics.cpp index 5a6cdfa..546cad6 100644 --- a/src/evb/FP_kinematics.cpp +++ b/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 DeltaZ(int ZT, int AT, int ZP, int AP, int ZE, int AE, double EP, double angle, double B) { /* 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; - - MT = MASS.FindMass(ZT, AT); - MP = MASS.FindMass(ZP, AP); - ME = MASS.FindMass(ZE, AE); - MR = MASS.FindMass(ZR, AR); + //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); - 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 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 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); //delta-Z in cm } - double Wire_Dist() {return 4.28625;} //cm - } diff --git a/src/evb/FP_kinematics.h b/src/evb/FP_kinematics.h index b3b147c..8813745 100644 --- a/src/evb/FP_kinematics.h +++ b/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 DeltaZ(int ZT, int AT, int ZP, int AP, int ZE, int AE, double EP, double angle, double B); - double Wire_Dist(); + static constexpr double WireDist() { return 4.28625; } //cm } diff --git a/src/evb/MassLookup.cpp b/src/evb/MassLookup.cpp index 972f5de..366caa6 100644 --- a/src/evb/MassLookup.cpp +++ b/src/evb/MassLookup.cpp @@ -11,6 +11,9 @@ Written by G.W. McCann Aug. 2020 #include "MassLookup.h" namespace EventBuilder { + + //instantiate + MassLookup* MassLookup::s_instance = new MassLookup(); /* Read in AMDC mass file, preformated to remove excess info. Here assumes that by default @@ -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; } } \ No newline at end of file diff --git a/src/evb/MassLookup.h b/src/evb/MassLookup.h index e6a0dc5..91ba8ae 100644 --- a/src/evb/MassLookup.h +++ b/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 { + 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(); ~MassLookup(); - double FindMass(int Z, int A); - std::string FindSymbol(int Z, int A); + 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: - std::unordered_map massTable; - std::unordered_map elementTable; + MassLookup(); + std::unordered_map m_dataMap; + static MassLookup* s_instance; + //constants - static constexpr double u_to_mev = 931.4940954; - static constexpr double electron_mass = 0.000548579909; - + 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 }; - //static instance for use throught program - static MassLookup MASS; - } #endif diff --git a/src/evb/SFPAnalyzer.cpp b/src/evb/SFPAnalyzer.cpp index 2196563..8cda322 100644 --- a/src/evb/SFPAnalyzer.cpp +++ b/src/evb/SFPAnalyzer.cpp @@ -17,7 +17,7 @@ namespace EventBuilder { SFPAnalyzer::SFPAnalyzer(int zt, int at, int zp, int ap, int ze, int ae, double ep, double angle, double b) { - zfp = Delta_Z(zt, at, zp, ap, ze, ae, ep, angle, b*1000.0); //Convert kG to G + zfp = DeltaZ(zt, at, zp, ap, ze, ae, ep, angle, b); event_address = new CoincEvent(); rootObj = new THashTable(); GetWeights(); @@ -42,7 +42,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); }