void pion_res(void)
{
	gROOT->Reset();
	//TFile *F1 = new TFile("pion.root");

	TChain chain("TRACKING/trkeff");
	chain.Add("pion1/pion1.root");
	chain.Add("pion2/pion2.root");
	chain.Add("pion3/pion3.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,30,0.,6.0,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,30,0.,6.0);
	for (int i=2;i<=30;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.05);
	  h->SetMinimum(0);

	  switch(i){
	  case 3:
	    h->Clone("h0");
	    break;
	  case 5: 
	    h->Clone("h1");
	   
	    break;  
	  case 8:
	    h->Clone("h2");
	    break;
	  case 13: 
	    h->Clone("h3");
	   
	    break;  
	  case 18:
	    h->Clone("h4");
	    break;
	  case 28: 
	    h->Clone("h5");
	   
	    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 pions");

	h0->Draw("l");
	legend->AddEntry(h0,"0.5 GeV/c","pl");

	/*
	  TH1D *h1=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(6);
	h2->SetLineColor(6);
	h2->Draw("same");
	legend->AddEntry(h2,"1.5 GeV/c","pl");
	
	/*
	TH1D *h3=gROOT->FindObject("h3");
	h3->SetMarkerStyle(20);	
	h3->SetMarkerSize(1.0);
	h3->SetMarkerColor(7);
	h3->SetLineColor(7);
	h3->Draw("same");
	legend->AddEntry(h3,"2.5 GeV/c","pl");
	*/
	
	TH1D *h4=(TH1D *)gROOT->FindObject("h4");
	h4->SetMarkerStyle(20);
	h4->SetMarkerSize(1.0);
	h4->SetMarkerColor(4);
	h4->SetLineColor(4);
	h4->Draw("same");
	legend->AddEntry(h4,"3.5 GeV/c","pl");
	
	TH1D *h5=(TH1D *)gROOT->FindObject("h5");
	h5->SetMarkerStyle(20);
	h5->SetMarkerSize(1.0);
	h5->SetMarkerColor(8);
	h5->SetLineColor(8);
	h5->Draw("same");
	legend->AddEntry(h5,"5.5 GeV/c","pl");
	legend->Draw();
	

	c1->SaveAs("pion_dp_vs_theta.png");	
	c1->SaveAs("pion_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 pions");

	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->GetYaxis()->SetLimits(0.,6.);
	res->Draw("colz");

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

	ofstream ofile("pion_res.txt"); 
	ofile <<"<br/>"<<endl;
	ofile <<"<img src=\"pion_dp_vs_theta.png\">" <<endl;
	ofile <<"<img src=\"pion_dp_vs_p_vs_theta.png\">" <<endl;

	ofile.close();
}