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

#include <TROOT.h>
#include <TStyle.h>
#include <TH1D.h>
#include <TTree.h>
#include <TCanvas.h>

double r2pdf(double r2)
{
   double z = r2 / 0.1849;
   return  1 / (1 + 3.6 * z + 1.11 * pow(z, 1.9));
}

TH1D *r2dist(TH1D *htemp)
{
   TH1D *hfit = (TH1D*)htemp->Clone("hfit");
   Long64_t nbins = hfit->GetNbinsX();
   for (int n=1; n <= nbins; ++n) {
      double z = exp(hfit->GetXaxis()->GetBinCenter(n));
      double dPdlnz = z / (1 + 3.6 * z + 1.25 * pow(z, 1.8));
      hfit->SetBinContent(n, dPdlnz);
   }
   double height1 = htemp->GetMaximum();
   double height0 = hfit->GetMaximum();
   hfit->Scale(height1 / height0);
   hfit->SetLineColor(kRed);

   gStyle->SetPadLeftMargin(0.15);
   TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1 == 0)
      c1 = new TCanvas("c1", "c1", 40, 20, 600, 500);
   hfit->Draw("same");
   return hfit;
}

TH1D *r2sim()
{
   TTree *h10 = (TTree*)gROOT->FindObject("h10");
   if (h10 == 0) {
      std::cerr << "Error in r2sim: unable to find input tree h10,"
                << " please open geant.root and try again."
                << std::endl;
      return 0;
   }
   TString name("log(x[0]*x[0]+x[1]*x[1])-2*log(0.43)");
   TH1D *hsim = new TH1D("hsim", "r{^2} distribution from simulation",
                         120, -10., 10.);

   gStyle->SetPadLeftMargin(0.15);
   TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1 == 0)
      c1 = new TCanvas("c1", "c1", 40, 20, 600, 500);
   h10->Draw(name + ">>hsim", "");
   return hsim;
}