NSDPhysicsCalcs.h

Go to the documentation of this file.

#ifndef _NSD_PHYSICS_CALCS_
#define _NSD_PHYSICS_CALCS_
 
#include <cmath>
#include <vector>
#include <numeric>
 
namespace NSDPhysicsCalcs
{
  //physics constants used in calculations
  constexpr double CcmPerNs = 29.979245800;
  constexpr double MnMeV = 939.56542052;
  constexpr double electronRestEnergy = 0.510998950; // MeV
 
class ComptonCalcs
{
public:
  // Computes and returns recoil gamma ray energy from Compton
  //    scattering kinematics
  inline static double recoilGammaEnergy(double a_incidentEnergy,
                                         double a_scatteringAngle);
 
  // Computes and returns recoil electron energy from Compton
  //    scattering kinematics
  inline static double recoilElectronEnergy(double a_incidentEnergy,
                                            double a_scatteringAngle);
 
  // Computes and returns scattering angle (degrees) of recoil gamma ray 
  //  from electron energy deposited in detector volume assuming electron is  
  //  fully stopped
  inline static double scatteringAngleFromED(double a_incidentEnergy, // MeV
                                             double a_electronEnergyDep); // MeV
 
};
 
 
class NTOFCalcs
{
public:
  NTOFCalcs();
  inline static double gT(double a_L
                         );
  inline static double Nen(double a_L,
                           double a_t
                           );
  inline static double NTOF(double a_L,
                            double a_nEn
                            );
  inline static std::vector<double> GenEnBins(double a_L,
                                              double a_maxT,
                                              double a_minT,
                                              double a_dt 
                                              );
};
 
class GenExtraCalcs
{
public:
  // GenExtraCalcs();
  inline static double levelThreshold(double A,
                                  double a_Ex
                                 );
 
  inline static double scintSigma(double a_E,
                         double a_A,
                         double a_B,
                         double a_C
                         ); 
  
  inline static double biExponential(double a_t,
                            double a_tau1,
                            double a_tau2
                            );
 
  inline static double triExponential(double a_t,
                            double a_tau1,
                            double a_tau2,
                            double a_tau3,
                            double a_lambda
                            );
  //Added by Ananya on 2/15/24 for CLYC WF fitting
  inline static double ErfEnvExponential(double a_t,
                            double a_tau1,
                            double a_tau2
                            );
  
  inline static double cloverResol(const std::vector<double>& a_coeffs,
                       double               a_val
                      );
  
  inline static double polyEval(const std::vector<double>& a_coeffs,
                       double               a_val
                      );
  inline static double getNormIntegral(double a_mu,
                                       double a_sigma,
                                       double a_start,
                                       double a_stop
                                      );
 
  inline static std::vector<double> getNormVector(std::vector<double> a_pntA,
                                                  std::vector<double> a_pntB,
                                                  std::vector<double> a_pntC
                                                  );
 
  inline static std::vector<double> getVector(double a_x1,
                                              double a_y1,
                                              double a_z1,
                                              double a_x2,
                                              double a_y2,
                                              double a_z2);
  inline static std::vector<double> getVector(std::vector<double> a_v1, 
                                              std::vector<double> a_v2);
  inline static double magnitude(std::vector<double> a_v);
 
  inline static double calcAngle(std::vector<double> a_v1,
                                 std::vector<double> a_v2);
  inline static double calcAngle(double a_x1,
                                 double a_y1,
                                 double a_z1,
                                 double a_x2,
                                 double a_y2,
                                 double a_z2
                                 );
  inline static double getGaussArea(double  a_norm,
                                    double  a_sig,
                                    double  a_binWidth,
                                    double  a_dNorm,
                                    double  a_dSig,
                                    double* a_dArea
                                    );
 
  inline static unsigned long upper_power_of_two(unsigned long v);
 
}; 
//~~~~~~~~~~~~~~~~~~implementation for inlines~~~~~~~~~~~~~~~~~~~~~~~~
// McGuire, 2023
 
// Calculates recoil gamma ray energy according to compton scattering kinematics
double ComptonCalcs::recoilGammaEnergy(double a_incidentEnergy, // MeV
                                       double a_scatteringAngle) // radians
{
  double A = 1.0;
  A *= a_incidentEnergy;
  double B = 1.0;
  B += (a_incidentEnergy/electronRestEnergy) * (1 - cos(a_scatteringAngle));
  double recoilGammaEnergy = A / B; // MeV
  return recoilGammaEnergy;
}
 
// Calculates recoil electron energy according to compton scattering kinematics
//    and conservation of energy
double ComptonCalcs::recoilElectronEnergy(double a_incidentEnergy, // MeV
                                          double a_scatteringAngle) // radians
{
  double recoilGammaRayEnergy = recoilGammaEnergy(a_incidentEnergy, 
                                               a_scatteringAngle); // MeV
  double recoilElectronEnergy = a_incidentEnergy - recoilGammaRayEnergy; // MeV
  return recoilElectronEnergy;
}
 
// Computes and returns scattering angle (degrees) of recoil gamma ray from 
// electron energy deposited in detector volume assuming electron is fully 
// stopped
// "this is a sligtly dangerous way to calculate the scatter angle
// what if the electron left the volume before depositing it's full energy?" -JB
double ComptonCalcs::scatteringAngleFromED(double a_incidentEnergy, // MeV
                                           double a_electronEnergyDep) // MeV
{
  double A = 1.0;
  double B = 1.0;
  B *= ((a_incidentEnergy / (a_incidentEnergy - a_electronEnergyDep)) - 1);
  B *= electronRestEnergy/a_incidentEnergy;
  A -= B;
  double theta = acos(A) * (180. / M_PI);; // degrees
  return theta; // degrees
}
 
 
double NTOFCalcs::gT(double a_L)
{
  return a_L/CcmPerNs;
}
double NTOFCalcs::Nen(double a_L,
                             double a_t
                            )
{
  double gTime = gT(a_L);
  if(a_t>=gTime)
  {
    double beta = (a_L*a_L)/(a_t*a_t*CcmPerNs*CcmPerNs);
    double gamma = 1.0/sqrt(1-beta);
    return MnMeV*(gamma -1.0);
  }
  //implicit else due to return if the flight time is faster than a gamma
  //this isn't a valid neutron time of flight return -100 for the energy
  return -100;
}
double NTOFCalcs::NTOF(double a_L,
                              double a_nEn
                              )
{
  double numerator = a_L * (a_nEn+MnMeV);
  double denominator = CcmPerNs * sqrt(a_nEn*a_nEn+2.0*MnMeV*a_nEn);
  return numerator/denominator;
}
std::vector<double> NTOFCalcs::GenEnBins(double a_L,
                             double a_maxT,
                             double a_minT,
                             double a_dt 
                          )
{
  std::vector<double> res;
  for(double t = a_maxT;t>a_minT;t-=a_dt)
  {
    res.push_back(Nen(a_L,t));
  }
  res.push_back(Nen(a_L,a_minT));
  return res;
}
 
 
double GenExtraCalcs::levelThreshold(double A,
                            double a_Ex
                            )
{
  return a_Ex*(A+1)/A ;
}
 
double GenExtraCalcs::cloverResol(const std::vector<double>& a_coeffs,
                double               a_val
               )
{
  return (a_coeffs[0]*a_val+a_coeffs[1])/(1+0.221*log(a_val)-1.07) \
                            + (0.221*log(a_val)-1.07)/(1+0.221*log(a_val)-1.07)*(a_coeffs[2]*a_val+2*a_coeffs[1]);
}
 
double GenExtraCalcs::polyEval(const std::vector<double>& a_coeffs,
                double               a_val
               )
{
  int a_coeffsSize = a_coeffs.size();
  if(a_coeffsSize>0)
  {
    double result = a_coeffs[0];
    // for(auto& coeff : a_coeffs)
    for(int a_i=1; a_i<a_coeffsSize; a_i++)
    {
      result = result*a_val + a_coeffs[a_i];
    }
    return result;
  }
  else
  {
    return 0;
  }
}
double GenExtraCalcs::scintSigma(double a_E,
                         double a_A,
                         double a_B,
                         double a_C
                         )
{
  return a_A+a_B*sqrt(a_E)+a_C*a_E;
} 
 
double GenExtraCalcs::biExponential(double a_t,
                            double a_tau1,
                            double a_tau2
                            )
{
  if(a_t>0)
  {
  return std::exp(-1.0 * a_t / a_tau2) * (1. - std::exp(-1.0 * a_t / a_tau1)) / 
                                          a_tau2 / a_tau2 * (a_tau1 + a_tau2);
  }
  return 0;
}
 
double GenExtraCalcs::ErfEnvExponential(double a_t,
                            double a_tau1,
                            double a_tau2
                            )
{
  if(a_t>0)
  {
  return std::exp(-1.0 * a_t / a_tau2) * std::erf(a_t / a_tau1);
  }
  return 0;
}
 
double GenExtraCalcs::triExponential(double a_t,
                            double a_tau1,
                            double a_tau2,
                            double a_tau3,
                            double a_lambda
                            )
{
  if(a_t>0)
  {
  return (a_lambda*std::exp(-1.0 * a_t / a_tau3)+std::exp(-1.0 * a_t / a_tau2)) * 
         (1. - std::exp(-1.0 * a_t / a_tau1)) *(a_tau2+a_tau3)*(a_tau1*a_tau2
          + a_tau1*a_tau3 + a_tau2*a_tau3)/((a_tau2*a_tau3)*(a_tau2*a_tau3));
  }
  return 0;
}
 
double GenExtraCalcs::getNormIntegral(double a_mu,
                                     double a_sigma,
                                     double a_start,
                                     double a_stop
                                    )
{
  double erfDenom = a_sigma*sqrt(2);
  double upperErfArg = (a_stop-a_mu)/erfDenom;
  double lowerErfArg = (a_start-a_mu)/erfDenom;
  return 0.5*(erf(upperErfArg)-erf(lowerErfArg));
}
 
inline std::vector<double> GenExtraCalcs::
                                  getNormVector(std::vector<double> a_pntA,
                                                std::vector<double> a_pntB,
                                                std::vector<double> a_pntC
                                                )
{
  //vector from a to b
  std::vector<double> vecAB;
  vecAB.push_back(a_pntB[0]-a_pntA[0]);
  vecAB.push_back(a_pntB[1]-a_pntA[1]);
  vecAB.push_back(a_pntB[2]-a_pntA[2]);
  std::vector<double> vecAC;
  //vector from a to c
  vecAC.push_back(a_pntC[0]-a_pntA[0]);
  vecAC.push_back(a_pntC[1]-a_pntA[1]);
  vecAC.push_back(a_pntC[2]-a_pntA[2]);
  //orthogonal vector is the cross product of the two resultant vectors
  std::vector<double> orthVector;
  orthVector.push_back(vecAB[1]*vecAC[2]-vecAB[2]*vecAC[1]);
  orthVector.push_back(vecAB[2]*vecAC[0]-vecAB[0]*vecAC[2]);
  orthVector.push_back(vecAB[0]*vecAC[1]-vecAB[1]*vecAC[0]);
  double vecNorm = sqrt(orthVector[0]*orthVector[0]+
                     orthVector[1]*orthVector[1]+
                     orthVector[2]*orthVector[2]
                     );
  orthVector[0]/=vecNorm;
  orthVector[1]/=vecNorm;
  orthVector[2]/=vecNorm;
  return orthVector;
}
 
inline std::vector<double> GenExtraCalcs::getVector(double a_x1,
                                                    double a_y1,
                                                    double a_z1,
                                                    double a_x2,
                                                    double a_y2,
                                                    double a_z2)
{
  double x = a_x2 - a_x1;
  double y = a_y2 - a_y1;
  double z = a_z2 - a_z1;
  std::vector<double> rtnVec = {x , y, z};
  return rtnVec;
}
 
inline std::vector<double> GenExtraCalcs::getVector(std::vector<double> a_v1,
                                                    std::vector<double> a_v2)
{
  return getVector(a_v1[0], a_v1[1], a_v1[2], a_v2[0], a_v2[1], a_v2[2]);
}
 
// WARING: Assumes vector is length 3
inline double GenExtraCalcs::magnitude(std::vector<double> a_v)
{
  return (std::sqrt(a_v[0]*a_v[0] + a_v[1]*a_v[1] + a_v[2]*a_v[2]));
}
 
inline  double GenExtraCalcs::calcAngle(std::vector<double> a_v1,
                                        std::vector<double> a_v2)
{
  //this needs to be completed
  double numerator = std::inner_product(a_v1.begin(), 
                                        a_v1.end(), a_v2.begin(), 0.);
  // another method below
  //double numnerator = a_v1[0]*a_v2[0] + a_v1[1]*a_v2[1] + a_v1[2]*a_v2[2];
  double denominator = magnitude(a_v1) * magnitude(a_v2);
  double cosTheta = numerator / denominator;
  double theta = acos(cosTheta); // radians
  //prefer radians as builtins all use radians
  return theta;
}
inline double GenExtraCalcs::calcAngle(double a_x1,
                                       double a_y1,
                                       double a_z1,
                                       double a_x2,
                                       double a_y2,
                                       double a_z2
                                       )
{
  double dot = a_x1*a_x2+a_y1*a_y2+a_z1*a_z2;
  double mag1 = sqrt(a_x1*a_x1+a_y1*a_y1+a_z1*a_z1);
  double mag2 = sqrt(a_x2*a_x2+a_y2*a_y2+a_z2*a_z2);
  return acos(dot/(mag1*mag2)); // radians
}
//this functon returns the area under a fitted gaussian and it's associated 
  //area when fitting binned data with a gaussian distribution
inline double GenExtraCalcs::
  getGaussArea(double  a_norm,
               double  a_sig,
               double  a_binWidth,
               double  a_dNorm,
               double  a_dSig,
               double* a_dArea
               )
{
  double area = a_norm*a_sig*sqrt(2.*M_PI)/a_binWidth;
  if(a_dArea!=nullptr)
  {
    double percDSig = a_dSig/a_sig;
    double percDNorm = a_dNorm/a_norm;
    *a_dArea = area * sqrt(percDNorm*percDNorm + percDSig*percDSig);
  }
  return area;
}
 
unsigned long GenExtraCalcs::upper_power_of_two(unsigned long v)
{
    v--;
    v |= v >> 1;
    v |= v >> 2;
    v |= v >> 4;
    v |= v >> 8;
    v |= v >> 16;
    v++;
    return v;
}
}//end of namespace
 
#endif