JDeltaRays.cc File Reference

Program to determine the energy loss due to visible delta-rays. More...

#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Program to determine the energy loss due to visible delta-rays.

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

Definition in file JDeltaRays.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 194 of file JDeltaRays.cc.

{
  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();
}