void proton_eff(void)
{
	gROOT->Reset();
	//TFile *F1 = new TFile("proton.root");
	TChain chain("TRACKING/trkeff");
	chain.Add("proton1/proton1.root");
	chain.Add("proton2/proton2.root");
	chain.Add("proton3/proton3.root");

	Double_t Red[5]   = { 0.5,1.00,0,0, 0 };
	Double_t Green[5] = { 0.1, 0.3, 1., 0., 0.};
	Double_t Blue[5]  = { 0.1, 0.35, 0.65,0.8,1.0 };	
	Double_t Length[5] = { 0.00, 0.3, 0.5, 0.85, 1.00 };
	Int_t nb=50;
	TColor::CreateGradientColorTable(5,Length,Red,Green,Blue,nb);

	// Create histograms to calculate efficiency
	TH2D *np = new TH2D("np","Numerator", 70,0,140,20,0,6);
	TH2D *dp = new TH2D("dp","Denominator",70,0,140,20,0,6);	
	TH2D *npcut = new TH2D("npcut","numerator",70,0,140,20,0,6);

	chain.Project("np", "F.pthrown.Mag():180./3.14159*F.pthrown.Theta()", "F.trktb.Nfdc>0 || F.trktb.Ncdc>0");
	chain.Project("dp", "F.pthrown.Mag():180/3.14159*F.pthrown.Theta()", "");
	chain.Project("npcut", "F.pthrown.Mag():180./3.14159*F.pthrown.Theta()", "(F.trktb.Nfdc>0 || F.trktb.Ncdc>0) && TMath::Prob(F.trktb.trk_chisq,F.trktb.trk_Ndof)>0.001");
	
	TH2D *effp = (TH2D *)np->Clone("effp");
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTickx();
	c1->SetTicky();
	effp->Divide(dp);
	effp->SetMarkerStyle(2);
	effp->SetStats(0);
	//gStyle->SetPalette(1,0);
	effp->SetContour(nb);
	effp->SetTitle("Proton reconstruction efficiency");
	
	effp->SetXTitle("#theta [degrees]");
	effp->SetYTitle("p [GeV/c]");

	effp->Draw("colz");

	c1->SaveAs("proton_efficiency.png");
	c1->SaveAs("proton_efficiency.pdf");

	TH2D *effpcut = (TH2D *)npcut->Clone("effpcut");
	
	TCanvas *c2 = new TCanvas("c2");
	c2->SetTickx();
	c2->SetTicky();
	effpcut->Divide(dp);
	effpcut->SetMarkerStyle(2);
	effpcut->SetStats(0);


	effpcut->SetContour(nb);
	effpcut->SetTitle("Proton reconstruction efficiency, P(#chi^2)>0.001");

	effpcut->SetXTitle("#theta [degrees]");
	effpcut->SetYTitle("p [GeV/c]");
	effpcut->Draw("colz");

	c2->SaveAs("proton_efficiency_with_chi2_cut.png");
	c2->SaveAs("proton_efficiency_with_chi2_cut.pdf");

	TH1D *prob = new TH1D("prob","prob",1000,0,1);
	chain.Project("prob","TMath::Prob(F.trktb.trk_chisq,F.trktb.trk_Ndof)");
	TCanvas *c3 = new TCanvas("c3");
	c3->SetTickx();
	c3->SetTicky();	
	c3->SetLogy();
	prob->SetTitle("Proton probability distribution");
	prob->SetXTitle("Confidence level for fit");
	prob->Draw();
	c3->SaveAs("proton_prob.png");	
	c3->SaveAs("proton_prob.pdf");



	ofstream ofile("proton_eff.txt");
	ofile<<"<h2>Proton events</h2>" <<endl;
	
	double eff = np->GetEntries()/dp->GetEntries();
	ofile <<"<p>Overall efficiency = "<<eff<<"<br/>"<<endl;
	double effcut = npcut->GetEntries()/dp->GetEntries();
	ofile <<"  Efficiency (Prob>0.001) = "<<effcut<<"<br/>" << endl;

	ofile << "<img src=\"proton_prob.png\"> <br/>" << endl;
	ofile <<"<img src=\"proton_efficiency.png\">" <<endl;
	ofile <<"<img src=\"proton_efficiency_with_chi2_cut.png\">" <<endl;
	
	ofile.close();
}