JDeltaRays.cc File Reference

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "JTools/JConstants.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.

Author

mdejong

Definition in file JDeltaRays.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 39 of file JDeltaRays.cc.

{
  using namespace std;
 
  typedef JTOOLS::JRange<double> JRange_t;
  
  string   outputFile;
  JRange_t T_GeV;
  int      numberOfPoints;
  string   lepton;
  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['T'] = make_field(T_GeV,          "kinetic energy range of electron [GeV]")  = JRange_t();
    zap['n'] = make_field(numberOfPoints, "points for integration")                  = 1000000;
    zap['L'] = make_field(lepton)         = muon, tau, positron, electron;
    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);
  }
 
  using namespace JPP;
 
  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 RATIO =  MASS_ELECTRON / MASS_LEPTON;
 
  
  // Minimal energy of electron for Cherenkov threshold [GeV]
 
  const double EMIN  =  MASS_ELECTRON * INDEX_OF_REFRACTION_WATER / sqrt(INDEX_OF_REFRACTION_WATER*INDEX_OF_REFRACTION_WATER - 1.0);
 
  // Minimal kinetic energy of electron [GeV]
 
  double Tmin  =  sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
  
  NOTICE("Tmin [GeV] " << FIXED(8,6) << Tmin << ' ' << FIXED(8,6) << T_GeV.getLowerLimit() << endl);
 
  if (T_GeV.is_valid()) {
    if (Tmin < T_GeV.getLowerLimit()) {
      Tmin = T_GeV.getLowerLimit();
    }
  }
 
  
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D h0("h0", NULL, 80, log10(MASS_LEPTON), +9.25);
  
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
 
    const double x     = h0.GetBinCenter(i);
    const double E     = pow(10.0,x);         // muon energy [GeV]
 
    const double gamma = E / MASS_LEPTON;
    const double beta  = sqrt(1.0 - 1.0/(gamma*gamma));
 
 
    // Maximal kinetic energy of electron [GeV]
    
    double Tmax  = (lepton == electron ?
                    (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
                    (2.0*MASS_ELECTRON*beta*beta*gamma*gamma)  /  (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
 
    if (T_GeV.is_valid()) {
      if (Tmax > T_GeV.getUpperLimit()) {
        Tmax = T_GeV.getUpperLimit();
      }
    }
 
 
    double value = 0.0;                       // energy of delta-rays
      
    if (Tmax > Tmin) {
    
      const double a =  1.0;                  // Form factor formula
      const double b = -beta*beta/Tmax;       // 
      const double c =  0.5/(E*E);            // F(T)  =  a + b*T + c*T^2
      
      const double xmin = log(Tmin);          // T * 1/T^2 * dT = d log(T)
      const double xmax = log(Tmax);
      const double dx   = (xmax - xmin) / numberOfPoints;
      
      const double W =  dx * 0.5 * K * (Z/A) * (1.0/(beta*beta));
      
      for (double x = xmin; x <= xmax; x += dx) {
        
        const double T = exp(x);
        
        const double F = a + T*(b + T*(c));   // Form factor formula
        const double y = W * F;
        
        value += y;
      }
    }
    
    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)      << value << endl);      
 
    h0.SetBinContent(i, value);
  }
 
 
  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
  
  h0.Fit(&f1, option.c_str(), "same");
 
 
  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 Tmin  =  sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
 
    double Emin  = 1 * MASS_LEPTON;
    double Emax  = 5 * MASS_LEPTON;
 
    for (double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {
 
      const double gamma = E / MASS_LEPTON;
      const double beta  = sqrt(1.0 - 1.0/(gamma*gamma));
      
      const double Tmax  = (lepton == electron ?
                            (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
                            (2.0*MASS_ELECTRON*beta*beta*gamma*gamma)  /  (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
 
      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();
}