gwm17/SPS_SABRE_EventBuilder
1
0
Fork 0
mirror of https://github.com/sesps/SPS_SABRE_EventBuilder.git synced 2024-11-22 10:08:50 -05:00
Code Issues Projects Releases Wiki Activity

Cleanup kinematics code and mass code. Add in correction for z-offset being taken from central trajectory.

Browse Source
This commit is contained in:
Gordon McCann 2023-05-19 16:03:30 -04:00
parent ca0a7066df
commit f053847fb9
6 changed files with 98 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/evb/EVBParameters.h
View File

@ -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
};
}

73
src/evb/FP_kinematics.cpp
View File

@ -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;
//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); //delta-Z in cm
}
double Wire_Dist() {return 4.28625;} //cm
}

6
src/evb/FP_kinematics.h
View File

@ -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
}

48
src/evb/MassLookup.cpp
View File

@ -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;
}
}

51
src/evb/MassLookup.h
View File

@ -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

4
src/evb/SFPAnalyzer.cpp
View File

@ -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);
}