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/JDeltaRays.hh"
 
#include "JAAnet/JAAnetToolkit.hh"
 
#include "JTools/JRange.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
  
  static const char* const electron = "electron";
  static const char* const positron = "positron";
  static const char* const muon     = "muon";
  static const char* const tau      = "tau";
}
 
 
/**
 * \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(JDeltaRays::getTmin()) : JDeltaRays::getTmin());
  
  TFile out(outputFile.c_str(), "recreate");
 
  const double xmin = -3.00; //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 ?
                         JDeltaRays::getTmax   (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 JDeltaRays::JFormFactor F(0.5/(E*E), -beta*beta/Tmax, 1.0);
 
      dEdx = JDeltaRays::getEnergyLoss(E, MASS_LEPTON, Tmin, Tmax, Z, A, F, numberOfPoints) * 1.0e3;  // [MeV g^-1 cm^2]
      dNdx = JDeltaRays::getCount     (E, MASS_LEPTON, Tmin, Tmax, Z, A, F, numberOfPoints);          // [g^-1 cm^2]
      
    } else {
 
      dEdx = JDeltaRays::getEnergyLoss(E, MASS_LEPTON, Tmin, Tmax, Z, A) * 1.0e3;                     // [MeV g^-1 cm^2]
      dNdx = JDeltaRays::getCount     (E, MASS_LEPTON, Tmin, Tmax, Z, A);                             // [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);
    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  = JDeltaRays::getTmax(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();
}