diff --git a/Armory/Kinematics.h b/Armory/Kinematics.h new file mode 100644 index 0000000..768dcdc --- /dev/null +++ b/Armory/Kinematics.h @@ -0,0 +1,229 @@ +#ifndef KINEMATICS_H +#define KINEMATICS_H +#include +#include +#include +#include +#include +const double u_MeV = 931.49410372; //u in MeV + +class Kinematics { +public: + /* + A(d,p)B is used as template, with A being beam, and p being ejectile. + Always, make m3 the thing you detect, and m1 the beam + */ + double m_A, m_d, m_p, m_B; + double E_beam; + +/** + * @{ \name List of funny impossible default values + */ +/** + * \brief Default values used for all the physics values + */ + double Q0=-9999, Qx=-9999; + double P4=-9999, E4=-9999, T4=-9999; //heavy recoil momentum, totE, KE + double P3=-9999, E3=-9999, T3=-9999; //light recoil + double ET=-9999; + double gamma4=-9999, beta4=-9999, theta4=-9999; //theta=heavy-recoil lab angle + double brho=-9999; +/** + * @} + */ + Kinematics(double m1, double m2, double m3, double m4, double ebeam) { + /* + A(d,p)B is used as template, with A being beam, and p being ejectile. + Always, make m3 the thing you detect, m2 the target, and m1 the beam + + ebeam is in MeV/u, all others are in amu + */ + m_A = m1; + m_d = m2; + m_p = m3, m_B = m4, E_beam = ebeam*m_A; + Q0 = (- m_B - m_p + m_d + m_A)*u_MeV; + } + Kinematics() {} + + double getbeamCM_KE() { + return E_beam*(m_d/(m_A+m_d)); + } + + double getEbeam_givenQ(double e3, double Ex, double angle3_deg) { + double angle = angle3_deg*TMath::Pi()/180.0; + double m1 = m_A * u_MeV; + //double m2 = m_d * u_MeV; + double m3 = m_p * u_MeV; + double m4 = m_B * u_MeV; // Total dynamic mass of particle 4 + double Q = Q0-Ex; + + double AA = m1-m4; + double BB = TMath::Sqrt(m1*m3*e3)*TMath::Cos(angle); + double CC = (m3+m4)*e3 - m4*Q; + + double result_sq1 = (BB - TMath::Sqrt(BB*BB-AA*CC))/(AA); + return result_sq1*result_sq1; + } + + void setValues(double m1, double m2, double m3, double m4, double ebeam) { + /* + Can be used to 'live update' say the beam energy in the case of active target detectors. + */ + + m_A = m1; + m_d = m2; + m_p = m3, m_B = m4, E_beam = ebeam*m_A; + Q0 = (- m_B - m_p + m_d + m_A)*u_MeV; + //std::cout << "Q0 MeV: " << Q0 << std::endl; + } + + void setEBeam(double Ebeam) {E_beam = Ebeam;} + + double getBeta4(double t3, double angle3); + double getTheta4(double t3, double angle3); + double getBrho(double t3, double angle3, double charge_state); + double getExc(double t3, double angle3); //t3 is proton energy detected in ORRUBA, angle3 is proton angle in degrees + + double getBeta4_fromvec(double t3, const TVector3 &pos, const TVector3 &origin) { + TVector3 local = pos - origin; //position w.r.t origin + float angle = local.Theta()*180./M_PI; + return getBeta4(t3, angle); + } + + double getExc_fromvec(double t3, const TVector3 &pos, const TVector3 &origin) { + TVector3 local = pos - origin; //position w.r.t origin + float angle = local.Theta()*180./M_PI; + return getExc(t3, angle); + } + + void setValuesFromFile(const std::string& filename) { + (void) filename; + /*std::ifstream in; + in.open(filename); + if(!in) { + std::cerr<< "File not open at " << filename << std::endl; + return; + } + for(std::string line; std::getline(in, line); ) { + if(line.size()!=0 && line[0]=='#') + ; //don't do anything with '#' lines + else { + std::stringstream ss(line); + ss>>m_A>>m_d>>m_p>>m_B>>E_beam; + } + } + in.close();*/ + + } +}; + +//double Kinematics::getQval(double m1, double m2, double m3, double t1, double t3, double angle3) +double Kinematics::getExc(double t3, double angle3) +/* + \brief Follows convention in Marion, 2013: (1 - beam, 2- target, 3-ejectile, 4-recoil) + + m1 is beam, (typically heavy nucleus) + m2 is 'd', (light target) + m3 is 'p', (light ejectile mass) + t1 is beam kinetic energy, + + All calculations are done here, and other wrapper functions written make derived quantities from stuff calculated here. + + \param t3 lab-kinetic-energy in MeV + \param angle3 lab-angle in deg of detected proton (if d,p) or other charged particle + \return Excitation energy of the heavy-recoil nucleus in MeV +*/ +{ + double m1 = m_A; + double m2 = m_d; + double m3 = m_p; + double m4 = m_B; + double t1 = E_beam; + + m1 *= u_MeV; + m2 *= u_MeV; + m3 *= u_MeV; + m4 *= u_MeV; + + double e1 = m1 + t1; + double e3 = m3 + t3; + ET = t1 + m1 + m2; + double p1 = TMath::Sqrt(t1*t1 + 2*m1*t1); + double p3 = TMath::Sqrt(t3*t3 + 2*m3*t3); + double cosTheta = TMath::Cos(angle3*TMath::Pi()/180.); +// return m1+m2-m3-TMath::Sqrt(m1*m1 + m2*m2 + m3*m3 + 2.*m2*e1 - 2.*e3*(e1+m2)+ 2.*p1*p3*cosTheta); + double Q = m1+m2-m3-TMath::Sqrt(m1*m1 + m2*m2 + m3*m3 + 2.*m2*e1 - 2.*e3*(e1+m2)+ 2.*p1*p3*cosTheta); + + Qx = Q; + //Recoil properties just in case it's useful + T4 = ET - e3 - (m1+m2-m3-Q); + P4 = TMath::Sqrt(T4*T4 + 2*m4*T4); + + //this angle will not be affected by eloss + theta4 = (180./M_PI)*TMath::ASin((p3/P4)*TMath::Sin(angle3*M_PI/180.)); + + //recalculate everything other than angle with lowered Kinetic energy + P4 = TMath::Sqrt(T4*T4 + 2*(m4+Q-Q0)*T4); + gamma4 = T4/m4+1.; + beta4 = TMath::Sqrt(1. - 1./(gamma4*gamma4)); + theta4 = (180./M_PI)*TMath::ASin((p3/P4)*TMath::Sin(angle3*M_PI/180.)); + + //if(m2==m3) std::cout << "Q0 elastic: " << Q0 << " " << Q << "\n" ; + + return Q0 - Q;//Q0 = Q + Exc + +} + +double Kinematics::getBeta4(double t3, double angle3) +/* + \brief Follows convention in Marion, 2013: (1 - beam, 2- target, 3-ejectile, 4-recoil) + + m1 is beam, (typically heavy nucleus) + m2 is 'd', (light target) + m3 is 'p', (light ejectile mass) + t1 is beam kinetic energy, + \param t3 lab-kinetic-energy in MeV + \param angle3 lab-angle in deg of detected proton (if d,p) or other charged particle + \return doppler-shift beta (=v/c) of the heavy-recoil nucleus, calls getExc() to fill the value +*/ +{ + getExc(t3, angle3); + return beta4; +} + +double Kinematics::getBrho(double t3, double angle3, double charge_state) { +/* + \brief Follows convention in Marion, 2013: (1 - beam, 2- target, 3-ejectile, 4-recoil) + + m1 is beam, (typically heavy nucleus) + m2 is 'd', (light target) + m3 is 'p', (light ejectile mass) + t1 is beam kinetic energy, + \param t3 lab-kinetic-energy in MeV + \param angle3 lab-angle in deg of detected proton (if d,p) or other charged particle + \param charge_state charge state of the intended nucleus, in units of elementary charge (= +1 for H+, +2 for He2+ etc). + \return b-rho value, generated from P4*3.3359e-3/charge_state where P4 is 4momentum of heavy-recoil calculated from orruba kinematics +*/ + getExc(t3,angle3); + return P4*3.3359e-3/charge_state; +} + +double Kinematics::getTheta4(double t3, double angle3) { +/* + \brief Follows convention in Marion, 2013: (1 - beam, 2- target, 3-ejectile, 4-recoil) + + m1 is beam, (typically heavy nucleus) + m2 is 'd', (light target) + m3 is 'p', (light ejectile mass) + t1 is beam kinetic energy, + \param t3 lab-kinetic-energy in MeV + \param angle3 lab-angle in deg of detected proton (if d,p) or other charged particle + \param charge_state charge state of the intended nucleus, in units of elementary charge (= +1 for H+, +2 for He2+ etc). + \return lab theta value of heavy recoil in degrees +*/ + getExc(t3,angle3); + return theta4; +} + + +#endif diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 1399c90..e8cb372 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -17,6 +17,7 @@ #include "EnergyLoss.h" // energy loss lookup between two positions in a medium #include "HistPlotter.h" #include "AutoHist2D.h" // auto-ranged, auto-binned 2D histograms written alongside tree1 +#include "Kinematics.h" #include #include #include @@ -406,11 +407,12 @@ int main(int argc, char **argv){ tree1->Branch("qqqUp", &qqqUp, "qqqUp/I"); tree1->Branch("qqqBk", &qqqBk, "qqqBk/I"); - double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN//, Ex_recon=NAN + double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN, Ex_recon=NAN ; tree1->Branch("EBeam_Kin", &EBeam_Kin_gs, "EBeam_Kin/D"); tree1->Branch("EBeam_Kin_2.2", &EBeam_Kin_2_2, "EBeam_Kin_2.2/D"); tree1->Branch("EBeam_Kin_3.4", &EBeam_Kin_3_4, "EBeam_Kin_3.4/D"); + tree1->Branch("Ex_recon", &Ex_recon, "Ex_recon/D"); // reconstructed angles from PW track fit, method 1 and 2 double reTheta, rePhi; @@ -670,15 +672,15 @@ int main(int argc, char **argv){ AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)"); AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(thetab * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)"); //tree1->Fill(); - /* + //Kinematics aakin_27Al(26.981538408,4.00260325413,4.0026035413,26.981538408,beam_energy_at_vertex/26.981538408); //m3 is alpha Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,beamEnergy/26.981538408); //m3 is proton //Kinematics apkin_27Al(1.00782503224,4.00260325413,4.00260325413,1.00782503224,beamEnergy/1.00782503224); //m3 is proton - Ex_recon = apkin_27Al.getExc(Tb, thetab); - EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Tb, 0.0, thetab); - EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Tb, 2.2, thetab); - EBeam_Kin_3_4 = apkin_27Al.getEbeam_givenQ(Tb, 3.4, thetab); - //std::cout << EBeam_Kin << std::endl;*/ + Ex_recon = apkin_27Al.getExc(Esx3, thetab); + EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Esx3, 0.0, thetab); + EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Esx3, 2.2, thetab); + EBeam_Kin_3_4 = apkin_27Al.getEbeam_givenQ(Esx3, 3.4, thetab); + //std::cout << EBeam_Kin << std::endl; }else if (qqqID >= 0){