BetheBloch.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file   BetheBloch.cxx
// \brief  Bethe Bloch dE/dx from tables or analytic formula
// \author Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "Calibration/func/BetheBloch.h"

#include "Utilities/func/MathUtil.h"

#include <algorithm>
#include <cassert>
#include <cstdio>

#include "cetlib/search_path.h"

#include "TF1.h"
#include "TGraph.h"

namespace calib
{
  // Generic physical constants:

  // From PDG
  const double Me  = 510.998910e3; // eV
  const double Mmu = 105.658367e6; // eV
  // From PDG "Passage of particles through matter"
  const double K = 0.307075; // MeV g^-1 cm^2
  const double j = 0.200;


  //......................................................................
  TGraph* IBetheBloch::GetdEdxGraph(double localDensity) const
  {
    TGraph* ret = new TGraph;
    // 1MeV to 100GeV
    for(double T = 1; T < 1e5; T *= 1.05){
      ret->SetPoint(ret->GetN(), T, dEdx(T)*localDensity);
    }

    return ret;
  }

  //......................................................................
  TGraph* IBetheBloch::GetMPVGraph(double localDensity, double dist) const
  {
    TGraph* ret = new TGraph;
    // 1MeV to 100GeV
    for(double T = 1; T < 1e5; T *= 1.05){
      ret->SetPoint(ret->GetN(), T,
                    MPV(T, dist, localDensity)*localDensity/dist);
    }

    return ret;
  }


  //......................................................................
  BetheBlochTables::BetheBlochTables(Material mat)
  {
    char junkStr[1024];
    double junk;

    // Load the points from the file
    std::string fname;
    cet::search_path sp("FW_SEARCH_PATH");
    switch(mat){
    case kPolyethylene:
      // http://pdg.lbl.gov/2010/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_221.dat
      sp.find_file("Calibration/func/muonloss_221.dat", fname);
      break;
    default:
      assert(0 && "Unknown material");
    }
    FILE* f = fopen(fname.c_str(), "r");
    assert(f);

    while(!feof(f)){
      double T, dEdx;
      int ret = fscanf(f, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
                       &T, &junk, &junk, &junk, &junk, &junk, &junk,
                       &dEdx, &junk, &junk, &junk);
      if(ret == 11){
        fT_pts.push_back(T_pt(T, dEdx));
      }
      else{
        // We couldn't parse it properly. Either preamble or a special line.
        fgets(junkStr, 1024, f);
      }
    }
    fclose(f);
  }

  bool ltTPt(BetheBlochTables::T_pt pt, double T)
  {
    return pt.T < T;
  }
  //......................................................................
  double BetheBlochTables::dEdx(double T) const
  {
    typedef std::vector<T_pt>::const_iterator it_t;

    // Find the point after T
    it_t itNext = std::lower_bound(fT_pts.begin(), fT_pts.end(), T, ltTPt);
    if(itNext == fT_pts.end()) --itNext;
    if(itNext == fT_pts.begin()) ++itNext;
    // And the point before
    it_t itPrev = itNext; --itPrev;

    // Interpolate
    const double T0 = itPrev->T; const double y0 = itPrev->dEdx;
    const double T1 = itNext->T; const double y1 = itNext->dEdx;
    return ((T1-T)*y0+(T-T0)*y1)/(T1-T0);
  }

  //......................................................................
  double BetheBlochTables::MIP() const
  {
    double mip = 1e10;
    for(unsigned int n = 0; n < fT_pts.size(); ++n){
      if(fT_pts[n].dEdx < mip) mip = fT_pts[n].dEdx;
    }
    return mip;
  }


  //......................................................................
  BetheBlochAnalytic::BetheBlochAnalytic(Material mat)
  {
    switch(mat){
    case kPolyethylene:
      // Numbers from the header of http://pdg.lbl.gov/2010/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_221.dat
      ZA   = 0.57034; // Z/A
      I    = 57.4; // eV
      a    = 0.1211;
      k    = 3.4292;
      x0   = 0.1489;
      x1   = 2.5296;
      Cbar = 3.0563;
      break;
    default:
      assert(0 && "Unkown material");
    }
  }

  //......................................................................
  double BetheBlochAnalytic::dEdx(double T) const
  {
    using util::sqr;

    T *= 1e6; // Convert to eV

    double beta, gamma;
    BetaGamma(T, beta, gamma);

    const double f = 2*Me*sqr(beta*gamma);

    const double Tmax = f/(1+2*gamma*Me/Mmu+sqr(Me/Mmu));

    const double delta = Delta(beta*gamma);

    const double dEdx = K*ZA/sqr(beta)*(log(f*Tmax/(sqr(I)))/2-sqr(beta)-delta/2);

    return dEdx;
  }

  //......................................................................
  double BetheBlochAnalytic::MPV(double T, double dist,
                                 double localDensity) const
  {
    using util::sqr;

    const double x = dist*localDensity;

    T *= 1e6; // Convert to eV

    double beta, gamma;
    BetaGamma(T, beta, gamma);

    const double f = 2*Me*sqr(beta*gamma);

    const double delta = Delta(beta*gamma);

    const double xi = K/2*ZA/sqr(beta)*x;
    const double mpv = xi*(log(f/I)+log(1e6*xi/I)+j-sqr(beta)-delta);

    return mpv;
  }

  //......................................................................
  double BetheBlochAnalytic::MIP() const
  {
    double mip = 1e20;
    // Between 100 MeV and 1GeV, with 0.1% accuracy
    for(double T = 1e2; T < 1e3; T *= 1.001){
      const double m = dEdx(T);
      if(m < mip) mip = m;
    }
    return mip;
  }

  //......................................................................
  void BetheBlochAnalytic::
  BetaGamma(double T, double& beta, double& gamma) const
  {
    gamma = (T+Mmu)/Mmu;
    beta = sqrt(1-1/(gamma*gamma));
  }

  //......................................................................
  double BetheBlochAnalytic::Delta(double betagamma) const
  {
    double delta = 0;
    const double x = log10(betagamma);
    if(x > x0){
      delta += 2*log(10)*x-Cbar;
      if(x < x1) delta += a*pow(x1-x, k);
    }
    return delta;
  }

}