void plot_piprob(void)
{
// 
// Plot
//

// #include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
// gStyle->SetOptStat(kFALSE);
// gStyle->SetOptFit(kTRUE);
// gStyle->SetOptFit(kFALSE);
gStyle->SetOptFit(11111);
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;

  Float_t xmin=2;
  Float_t xmax=4;

   TF1 *piprob = new TF1("piprob",piprob,xmin,xmax,4);

   Double_t ratio_mupi=10;
   Double_t ntrk=1;
   Double_t frac=1;
   // Double_t frac=0.5;
   Double_t alpha=0.0013;
   piprob->SetParameters(ntrk,ratio_mupi,frac,alpha);
   
   TCanvas *c1 = new TCanvas("c1","c1 Pion Probability",200,10,700,350);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->Divide(2,1);
   c1->SetGridx();
   c1->SetGridy();
   // c1->SetLogy();

   c1->cd(1);
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);
   piprob1=piprob->DrawCopy();

   piprob1->GetXaxis()->SetRangeUser(xmin,xmax);
   // piprob1->GetYaxis()->SetRangeUser(ymin,ymax);
   piprob1->GetXaxis()->SetTitleSize(0.05);
   piprob1->GetYaxis()->SetTitleSize(0.05);
   piprob1->GetXaxis()->SetTitle("Momentum (GeV)");
   piprob1->GetYaxis()->SetTitle("P_{#pi} for one track");
   piprob1->GetYaxis()->SetTitleOffset(1.5);
   piprob1->GetXaxis()->SetNdivisions(505);
   piprob1->SetLineColor(2);
   piprob1->SetTitle("");


   sprintf (string,"#sigma_{#mu}/#sigma_{#pi}=%.1f\n",ratio_mupi);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.8,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();
   sprintf (string,"#alpha_{0}=%.4f\n",alpha);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();


   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);


   // Double_t ratio_mupi=10;
   Double_t ntrk=2;
   Double_t frac=1;
   // Double_t frac=0.5;
   Double_t alpha=0.0013;
   piprob->SetParameters(ntrk,ratio_mupi,frac,alpha);


   piprob2=piprob->DrawCopy();
   piprob2->GetXaxis()->SetRangeUser(xmin,xmax);
   // piprob2->GetYaxis()->SetRangeUser(ymin,ymax);
   piprob2->GetXaxis()->SetTitleSize(0.05);
   piprob2->GetYaxis()->SetTitleSize(0.05);
   piprob2->GetXaxis()->SetTitle("Momentum (GeV)");
   piprob2->GetYaxis()->SetTitle("P_{#pi} for two tracks");
   piprob2->GetYaxis()->SetTitleOffset(1.5);
   piprob2->GetXaxis()->SetNdivisions(505);
   piprob2->SetLineColor(2);
   piprob2->SetTitle("");


   sprintf (string,"#sigma_{#mu}/#sigma_{#pi}=%.1f\n",ratio_mupi);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.8,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();
   sprintf (string,"#alpha_{0}=%.4f\n",alpha);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();


        sprintf(filename,"plot_piprob_c1_frac%.1f.pdf",frac);
        c1->SaveAs(filename);
        sprintf(filename,"plot_piprob_c1_frac%.1f.pdf",frac);
        c1->SaveAs(filename);

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

   // Returns the probability that the event is a 2pi event

   Double_t ntrk=par[0];
   Double_t ratio_mupi=par[1];
   Double_t frac=par[2];
   Double_t alpha=par[3];
   Double_t x1=x[0];

   Double_t ctau=780;
   Double_t ldecay=800;
   Double_t mpi=0.14;
   Double_t ppi = x1;
   Double_t gamma_beta = ppi/mpi;
   Double_t eff = exp(-ldecay/(gamma_beta*ctau));
  
   // eff = 1 - frac*(1-eff);

   char string[256];
   Double_t pi=3.14159;

   // calculation assuming one track is identified

   Double_t P1_given_pi = eff;
   Double_t Ppi = 1;
   Double_t P1 = 1*eff + ratio_mupi*alpha;
   Double_t Ppi_given_1 = P1_given_pi * Ppi / P1;


   sprintf(string,"ppi=%g,  ratio_mupi=%g, eff=%g,P1_given_pi=%g, Ppi_given_1=%g\n",ppi,ratio_mupi,eff,P1_given_pi,Ppi_given_1);
   printf ("func: %s",string);


   // calculation assuming two tracks are identified

   Double_t P2_given_pi = eff;
   Double_t Ppi = Ppi_given_1;
   Double_t P2 = Ppi*eff + (1-Ppi)*alpha;
   Double_t Ppi_given_2 = P2_given_pi * Ppi / P2;

   sprintf(string,"ppi=%g,  ratio_mupi=%g, eff=%g,Ppi=%g,P2_given_pi=%g, Ppi_given_2=%g\n",ppi,ratio_mupi,eff,Ppi,P2_given_pi,Ppi_given_2);
   printf ("func: %s\n",string);

   if (ntrk <=1) return Ppi_given_1;
   if (ntrk <=2) return Ppi_given_2;
}