//
// EMshowerdep.C - a root macro for estimating the average shower energy
//                 deposition profile (longitudinal) in units of radiation
//                 lengths. The result is returned in Gev, the same units
//                 as used for E0, such that Int[0,inf]{dEdx(x) dx} = E0.
//
// author: richard.t.jones at uconn.edu
// version: november 19, 2014
//

#include <iostream>
#include <math.h>
#include <TF1.h>

Double_t dEdx(Double_t *x, Double_t *par)
{
   double E0_GeV = par[0];
   double Z = par[1];
   double Ecrit_GeV = 0.61 / (Z + 1.24);    // PDG booklet Sect. 27.4.3
   double y = E0_GeV / Ecrit_GeV;
   double xmax = log(y) + 0.5;
   double b = 0.5;
   double a = xmax * b + 1;
std::cout << "a is " << a << std::endl;
   return E0_GeV * b * pow(b * x[0], a - 1) * exp(-b * x[0]) / tgamma(a);
}

void plot_dEdx(double E0_GeV, double Z)
{
   TF1 *f1 = new TF1("dEdx_GeV",dEdx,0,25.,2);
   f1->SetParameters(E0_GeV, Z);
   f1->SetParNames("endpoint energy (GeV)", "average atomic number");
   f1->Draw();
}