{
#include <math.h>;
c1 = new TCanvas("dE/dX","Graph with error bars",200,100,700,500);
c1->SetGrid();
//- Bethe-Bloch formula------------
double NA= 6.022e23;  // [1/mol],=1/g
double re = 2.818e-15; // [m]
double mec2 = 0.511;  // MeV
double z = 1.;
double Z = 82.;
double A = 207.2;
double beta;
double I = 16.*pow(Z,0.9);   // [eV]
double pi = 3.141593; 
//
for (int k=0; k<4; k++){
double mass = 0.105; //  muon
if(k==1) mass = 0.140;
if(k==2) mass = 0.49;
if(k==3) mass = 0.838;

double E[80], p[80], BB[80];
for (int i=0; i<70; i++){
   p[i] = 0.05*pow(10., i/20.);
   E[i] = sqrt(p[i]*p[i] + mass*mass);
   double beta = p[i]/E[i];
   BB[i] = 4.*pi*NA*re*re*mec2*z*z*Z/A/beta/beta*1e4*
    (log(2.*mec2*1.e6*beta*beta/I/(1.-beta*beta))-beta*beta);
//--- dimension=1/g*m*m*MeV=MeV/(g/m/m)=1.e4*MeV/(g/cm/cm)--------------
//----order estimate-----
//------4*3*6e23*3e-15*3e-15*0.5*82/207*1e4=128*e(23-30+4)=0.128 MeV/(g/cm/cm)
    cout<<p[i]<<" "<<E[i]<<" "<<BB[i]<<endl;
}
//
hBB = new TGraph(70, p, BB);
gPad->SetLogx();
hBB->GetXaxis()->SetTitle("particle momentum (GeV/c)");
hBB->GetYaxis()->SetTitle("dE/dX (MeV g^{-1}cm^{2}))");
hBB ->SetMaximum(2.);
hBB->SetMinimum(1.);
hBB->SetLineColor(k+1);
if(k==0) hBB->Draw("AL");
else hBB->Draw("SAME L");
}
//
//-TO BE replaced with Bethe-Bloch formula------------
c1->Print("Bethe-Bloch.png");
}