void proton_res()
{
	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");


	// Create histograms to calculate resolution
	TH3D *dp_vs_p_vs_theta = new TH3D("dp_vs_p_vs_theta","Normalized momentum residuals vs p vs theta",75,0,130,10,0.,2.,100,-0.5,0.5);
	chain.Project("dp_vs_p_vs_theta", "(sqrt(F.trktb.p.fX*F.trktb.p.fX+F.trktb.p.fY*F.trktb.p.fY+F.trktb.p.fZ*F.trktb.p.fZ)-F.pthrown.Mag())/F.pthrown.Mag():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&&F.pthrown.Theta()>0.017");
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTickx();
	c1->SetTicky();
	gStyle->SetOptStat(0);


	TLegend *legend=new TLegend(0.65,0.7,0.85,0.95);
	TH2D *res=new TH2D("res","Momentum resolution",75,0,130,10,0.,2.);
	for (int i=2;i<=10;i++){
	  dp_vs_p_vs_theta->GetYaxis()->SetRange(i,i);
	  dp_vs_p_vs_theta->Project3D("zx"); 
	  TH2D *hzx=(TH2D *)gROOT->FindObject("dp_vs_p_vs_theta_zx");
	  hzx->FitSlicesY();
	  TH1D *h=(TH1D *)gROOT->FindObject("dp_vs_p_vs_theta_zx_2");
	  h->SetMaximum(0.1);
	  h->SetMinimum(0);

	  switch(i){
	  case 3:
	    h->Clone("h0");
	    break;
	  case 5: 
	    h->Clone("h1");
	   
	    break;  
	  case 8:
	    h->Clone("h2");
	    break;
	  default:
	    break;
	  }

	  for (int j=1;j<=75;j++){
	    res->SetBinContent(j,i,h->GetBinContent(j));
	  }
	  
	}	
	gPad->SetRightMargin(0.05);

	TH1D *h0=(TH1D *)gROOT->FindObject("h0");
	h0->SetMarkerStyle(20);
	h0->SetMarkerSize(1.0);
	h0->SetMarkerColor(1);
	h0->SetLineColor(1);
	h0->SetXTitle("#theta [degrees]");
	h0->SetYTitle("#sigma(#deltap/p)");
	h0->GetYaxis()->SetTitleOffset(1.25);
	h0->SetTitle("Momentum resolution for protons");
	h0->Draw("l");
	legend->AddEntry(h0,"0.5 GeV/c","pl");

     
	TH1D *h1=(TH1D *)gROOT->FindObject("h1");
	h1->SetMarkerStyle(20);
	h1->SetMarkerSize(1.0);
	h1->SetMarkerColor(2);
	h1->SetLineColor(2);
	h1->Draw("same");
	legend->AddEntry(h1,"0.9 GeV/c","pl");
	
	TH1D *h2=(TH1D *)gROOT->FindObject("h2");
	h2->SetMarkerStyle(20);	
	h2->SetMarkerSize(1.0);
	h2->SetMarkerColor(4);
	h2->SetLineColor(4);
	h2->Draw("same");
	legend->AddEntry(h2,"1.5 GeV/c","pl");
	legend->Draw();
	
	c1->SaveAs("proton_dp_vs_theta.png");	
	c1->SaveAs("proton_dp_vs_theta.pdf");

	TCanvas *c2 = new TCanvas("c2");
	c2->SetTickx();
	c2->SetTicky();
	gStyle->SetPalette(1,0);	
	gPad->SetLeftMargin(0.1);
	gPad->SetRightMargin(0.15);

	res->SetTitle("Momentum resolution for protons");
	res->SetZTitle("#sigma(#deltap/p)");
	res->SetXTitle("#theta [degrees]");
	res->SetYTitle("p [GeV/c]");
	res->SetMaximum(0.1);	
	res->SetContour(50);
	res->GetZaxis()->SetTitleOffset(1.25);
	res->Draw("colz");	

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

	ofstream ofile("proton_res.txt");
	ofile <<"<br/>"<<endl;
	ofile <<"<img src=\"proton_dp_vs_theta.png\">" <<endl;
	ofile <<"<img src=\"proton_dp_vs_p_vs_theta.png\">" <<endl;
	ofile <<"</html>" <<endl;
	
	ofile.close();

}