JDeltaRays.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "TFitResult.h"

#include "JROOT/JMinimizer.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JPDFToolkit.hh"

#include "JAAnet/JAAnetToolkit.hh"

#include "JTools/JRange.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {
  
  using namespace JPP;
  
  static const char* const electron = "electron";
  static const char* const positron = "positron";
  static const char* const muon     = "muon";
  static const char* const tau      = "tau";


  /**
   * Auxiliary class for 2nd order polynomial form factor.
   */
  struct JFormFactor
  {
    /**
     * Constructor.
     *
     * \param  a   second order polynomial coefficient [GeV^-2]
     * \param  b   first  order polynomial coefficient [GeV^-1]
     * \param  c   zeroth order polynomial coefficient [unit]
     */
    JFormFactor(const double a,
                const double b,
                const double c) :
      a(a), b(b), c(c)
    {}


    /**
     * Get form factor for given delta-ray kinetic energy.
     *
     * \param  T   delta-ray kinetic energy [GeV]
     * \return     form factor              [unit]
     */
    double operator()(const double T) const
    {
      return c + T*(b + T*a);
    }

    
  private:
    double a; // second order polynomial coefficient [GeV^-2]
    double b; // first  order polynomial coefficient [GeV^-1]
    double c; // zeroth order polynomial coefficient [unit]
  };
  

  /**
   * Get number of delta-rays due to delta-rays per unit track length\n
   * for an ionising particle with given energy and given mass.
   *
   * N.B: For this function a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$ is assumed, where\n:
   *      - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
   *      - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
   *      - \f$E\f$ to the ionising particle's energy and
   *      - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
   *
   * \param  E      particle energy                  [GeV]
   * \param  M      particle mass                    [GeV]
   * \param  Tmin   minimum delta-ray kinetic energy [GeV]
   * \param  Tmax   maximum delta-ray kinetic energy [GeV]
   * \param  Z      atomic number                    [unit]
   * \param  A      atomic mass                      [g/mol]
   * \return        delta-ray count                  [m^-1]
   */
  inline double getCount(const double E,
                         const double M,
                         const double Tmin,
                         const double Tmax,
                         const double Z,
                         const double A)
  {      
    if (Tmax > Tmin) {

      const double K = 0.307075;                                         // [MeV mol^-1 cm^2]
    
      const double gamma = E / M;
      const double beta  = sqrt((1.0 - 1.0/gamma) * (1.0 + 1.0/gamma));
      
      const double a = 0.5/(E*E);
      const double b = beta*beta/Tmax;
      const double c = 1.0;

      const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta));              // [MeV g^-1 cm^2]

      const double dT = Tmax - Tmin;
      const double rT = Tmax / Tmin;
      const double mT = Tmax * Tmin;

      return W * (a*dT - b*log(rT) + c*dT/mT) * DENSITY_SEA_WATER * 1.0e2; // [m^-1]
      
    } else {

      return 0.0;
    }
  }


  /**
   * Get number of delta-rays due to delta-rays per unit track length\n
   * for an ionising particle with given energy and given mass.\n\n
   *
   * The template parameter corresponds to a class which contains an `operator()(const double)`\n
   * to compute the form factor corresponding to a given delta-ray kinetic energy.
   *
   * \param  E      particle energy                                     [GeV]
   * \param  M      particle mass                                       [GeV]
   * \param  Tmin   minimum delta-ray kinetic energy                    [GeV]
   * \param  Tmax   maximum delta-ray kinetic energy                    [GeV]
   * \param  Z      atomic number                                       [unit]
   * \param  A      atomic mass                                         [g/mol]
   * \param  F      form factor functor
   * \param  N      number of points for numeric integration
   * \return        delta-ray count                                     [m^-1]
   */
  template<class JFormFactor_t>
  inline double getCount(const double         E,
                         const double         M,
                         const double         Tmin,
                         const double         Tmax,
                         const double         Z,
                         const double         A,
                         const JFormFactor_t& F,
                         const int            N = 1000000)
  {
    using namespace JPP;

    if (Tmax > Tmin) {

      const double K = 0.307075;                            // [MeV mol^-1 cm^2]

      const double gamma = E / M;
      const double beta  = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));

      const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)); // [MeV g^-1 cm^2]

      const double xmin = 1.0/Tmax;
      const double xmax = 1.0/Tmin;
      const double dx   = (xmax - xmin) / ((double) N);
        
      double count = 0.0;

      for (double x = xmin; x <= xmax; x += dx) {

        const double T = 1.0/x;
        const double y = W * F(T) * dx;                     // [MeV g^-1 cm^-2]
          
        count += y;
      }

      return count * DENSITY_SEA_WATER * 1.0e2;             // [m^-1]
        
    } else {

      return 0.0;
    }
  }
}


/**
 * \file
 *
 * Program to determine the energy loss due to visible delta-rays.\n
 * Energy loss due to delta-rays and maximal kinetic energy (Tmax) are taken from reference
 * https://pdg.lbl.gov/2020/reviews/rpp2020-rev-passage-particles-matter.pdf
 * \author mdejong, bjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef  JRange<double> JRange_t;

  string   outputFile;
  JRange_t T_GeV;
  int      numberOfPoints;
  bool     use_numerical;  
  string   lepton;
  double   A;
  double   Z;
  string   option;
  double   precision;
  int      debug;

  try {

    JParser<> zap("Program to determine the energy loss due to visible delta-rays.");

    zap['o'] = make_field(outputFile)     = "delta-rays.root";
    zap['n'] = make_field(numberOfPoints, "points for integration")                  = 1000000;
    zap['N'] = make_field(use_numerical,  "perform numeric integration");
    zap['T'] = make_field(T_GeV,          "kinetic energy range of electron [GeV]")  = JRange_t();
    zap['L'] = make_field(lepton)         = muon, tau, positron, electron;
    zap['A'] = make_field(A,              "atomic mass")                             = 18.0;
    zap['Z'] = make_field(Z,              "atomic number")                           = 10.0;    
    zap['O'] = make_field(option)         = "";
    zap['e'] = make_field(precision)      = 1.0e-6;
    zap['d'] = make_field(debug)          = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  if (option.find('R') == string::npos) { option += 'R'; }
  if (option.find('S') == string::npos) { option += 'S'; }
  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }

  double MASS_LEPTON;

  if      (lepton == muon)
    MASS_LEPTON = MASS_MUON;
  else if (lepton == tau)
    MASS_LEPTON = MASS_TAU;
  else if (lepton == positron)
    MASS_LEPTON = MASS_ELECTRON;
  else if (lepton == electron)
    MASS_LEPTON = MASS_ELECTRON;
  else
    FATAL("Invalid lepton " << lepton << endl);  

  const double Tmin = (T_GeV.is_valid() ?
                       T_GeV.constrain(getDeltaRayTmin()) : getDeltaRayTmin());
  
  TFile out(outputFile.c_str(), "recreate");

  const double error = 1.0e-12;

  const double xmin = log10(MASS_LEPTON);
  const double xmax = 9.25;

  TH1D h0("h0", NULL, 320, xmin, xmax);
  TH1D h1("h1", NULL, 320, xmin, xmax);

  for (int i = 1; i <= h0.GetNbinsX(); ++i) {

    const double x    = h0.GetBinCenter(i);
    const double E    = pow(10.0,x);         // particle energy [GeV]

    double       Tmax = (lepton != electron ?
                         getDeltaRayTmax       (E, MASS_LEPTON) :                        
                         0.5 * getKineticEnergy(E, MASS_ELECTRON));
    
    if (T_GeV.is_valid()) { Tmax = T_GeV.constrain(Tmax); }

    double dEdx = 0.0;
    double dNdx = 0.0;
    
    if (use_numerical) {

      const double gamma = E / MASS_LEPTON;
      const double beta  = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));

      const JFormFactor F(0.5/(E*E), -beta*beta/Tmax, 1.0);

      dEdx = getDeltaRays(E, MASS_LEPTON, Tmin, Tmax, Z, A, F, numberOfPoints) / DENSITY_SEA_WATER * 1e1;  // [MeV g^-1 cm^2]
      dNdx = getCount    (E, MASS_LEPTON, Tmin, Tmax, Z, A, F, numberOfPoints) / DENSITY_SEA_WATER * 1e-2; // [g^-1 cm^2]
      
    } else {

      dEdx = getDeltaRays(E, MASS_LEPTON, Tmin, Tmax, Z, A) / DENSITY_SEA_WATER * 1e1;  // [MeV g^-1 cm^2]
      dNdx = getCount    (E, MASS_LEPTON, Tmin, Tmax, Z, A) / DENSITY_SEA_WATER * 1e-2; // [g^-1 cm^2]
    }
    
    DEBUG("E     [GeV]            " << SCIENTIFIC(8,2) << E    << endl);
    DEBUG("Tmin  [GeV]            " << SCIENTIFIC(8,2) << Tmin << endl);
    DEBUG("Tmax  [GeV]            " << SCIENTIFIC(8,2) << Tmax << endl);
    DEBUG("dE/dx [MeV g^-1 cm^2]  " << FIXED(5,4)      << dEdx << endl);
    DEBUG("dE/dx [g^-1 cm^2]      " << FIXED(5,2)      << dNdx << endl);

    h0.SetBinContent(i, dEdx);
    h0.SetBinError  (i, error);
    h1.SetBinContent(i, dNdx);
  }

  if (use_numerical) {
    
    TF1 f1("f1", "[0] + [1]*x + [2]*x*x + [3]*x*x*x");
      
    if (option.find('W') == string::npos) {
      option += "W";
    }
  
    f1.SetParameter(0,  h0.GetMinimum());
    f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
    f1.SetParameter(2,  0.0);
    f1.SetParameter(3,  0.0);

    f1.SetLineColor(kRed);
  
    // fit
  
    const TFitResultPtr result = h0.Fit(&f1, option.c_str(), "same", xmin, xmax);

    cout << "chi2/NDF " << result->Chi2() << '/' << result->Ndf() << endl;
  
    cout << "\t// " << lepton << endl;
    cout << "\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
  
    for (int i = 0; i != 4; ++i) {
      cout << "\tstatic const double " << (char) ('a' + i) << " = " << SCIENTIFIC(10,3) << f1.GetParameter(i) << ";" << endl;
    }

    double Emin  = 1 * MASS_LEPTON;
    double Emax  = 5 * MASS_LEPTON;

    for (double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {

      const double Tmax  = getDeltaRayTmax(E, MASS_LEPTON);

      if (fabs(Tmax - Tmin) < precision) {
        cout << "\tstatic const double Emin = " << FIXED(7,5) << E << "; // [GeV]" <<  endl;
        break;
      }
      
      if (Tmax > Tmin)
        Emax = E;
      else
        Emin = E;
    }
  }
  
  out.Write();
  out.Close();
}