void test_shower_profile(void)
{
// 
// generate the radial profile of an EM shower using Eq. 27.35 of PDG 2008
//

#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
gStyle->SetOptFit(kTRUE);
// gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 100000;

   // define histograms

   Int_t nbins=100;

   //
   TCanvas *c2 = new TCanvas("c2","c2 Test random",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   // c2->SetGridx();
   // c2->SetGridy();
   // c2->SetLogy();

   // define input masses and kinematics

  Double_t Rpar=2;

   TF1 *Radial = new TF1("Radial",Radial,0,5,1);
   Radial->SetParameter(0,Rpar);
   Double_t ymin = 0;
   Double_t ymax = 300;
   // Radial->GetYaxis()->SetRangeUser(ymin,ymax);
   Radial->DrawCopy();

   //
   TCanvas *c3 = new TCanvas("c3","c3 Test random",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   //c3->SetGridx();
   //c3->SetGridy();
   //c3->SetLogy();

   TH1F *Rdistr = new TH1F("Rdistr","Rdistr",nbins,0,5);

   for (j=0;j<npts;j++) {
      Double_t  rand = Radial->GetRandom();
      Rdistr->Fill(rand);
      }

   Rdistr->SetTitle("");
   // Rdistr->GetXaxis()->SetRangeUser(1.0,2.5);
   // Rdistr->GetYaxis()->SetRangeUser(ymin,ymax);
   Rdistr->GetXaxis()->SetTitleSize(0.05);
   Rdistr->GetYaxis()->SetTitleSize(0.05);
   Rdistr->GetYaxis()->SetTitleOffset(1.5);
   Rdistr->GetYaxis()->SetTitle("events");
   Rdistr->GetXaxis()->SetTitle("Radius (cm)");
   Rdistr->GetXaxis()->SetNdivisions(505);
   Rdistr->SetMarkerColor(2);
   Rdistr->SetMarkerStyle(21);
   Rdistr->Draw();
      

}
Double_t Radial(Double_t *x, Double_t *par) {

   // generate radial EM shower profile distribution
   // See PDG 2008 eq. 27.35
   // Single parameter is Rpar phenomenological function of x/X0 and lnE
   // r is the radial variable

   char string[256];

   Double_t Rpar=par[0];
   Double_t r=x[0];

   sprintf(string,"Rpar=%f, r=%f\n",Rpar, r);
   // printf ("Radial: %s",string);
  
   Double_t den = r*r + Rpar*Rpar;
   Double_t distr = 2*r*Rpar*Rpar/(den*den);

   return distr;

}