#include "StandardLabels.C"

#include "filelist.h"

void pmean_vs_theta_compare(void)
{
	gROOT->Reset();

	TH1D *h1 = GetResHisto("with_muls_no_eloss", "h1", kRed);
	TH1D *h2 = GetResHisto("with_muls_no_eloss_spoiled", "h2", kBlack);
	TH1D *h3 = GetResHisto("no_muls_no_eloss", "h3", kMagenta);
	TH1D *h4 = GetResHisto("no_muls_no_eloss_spoiled", "h4", kBlue);

	TCanvas *c1 = new TCanvas("c1");
	c1->SetFillColor(kWhite);
	c1->SetTicks();
	c1->SetGrid();

	TH2D *axes = new TH2D("axes","Momentum Systematic Error vs. #theta", 100, 0.0, 40.0, 100, -0.03, 0.04);
	axes->SetStats(0);
	axes->SetYTitle("Total Momentum systematic #deltap/p");
	axes->SetXTitle("#theta (degrees)");
	axes->Draw();
	
	h1->Draw("same");
	h2->Draw("same");
	h3->Draw("same");
	h4->Draw("same");
	
	TLegend *leg = new TLegend(0.50, 0.55, 0.85, 0.89);
	leg->SetFillColor(kWhite);
	leg->AddEntry(h1, "Nominal field with MULS");
	leg->AddEntry(h2, "Spoiled field with MULS");
	leg->AddEntry(h3, "Nominal field w/o MULS");
	leg->AddEntry(h4, "Spoiled field w/o MULS");
	leg->Draw();
	
	AddStandardLabels(axes);
	
	c1->SaveAs("pmean_vs_theta_compare.gif");
	c1->SaveAs("pmean_vs_theta_compare.pdf");
}

TH1D* GetResHisto(string nam, const char *hname, int color)
{
	TChain *track = new TChain("TRACKING/track", "track");
	
	if(nam=="with_muls_no_eloss")
		AddFilesToChain_with_muls_no_eloss(track);
	else if(nam=="with_muls_no_eloss_spoiled")
		AddFilesToChain_with_muls_no_eloss_spoiled(track);
	else if(nam=="no_muls_no_eloss")
		AddFilesToChain_no_muls_no_eloss(track);
	else if(nam=="no_muls_no_eloss_spoiled")
		AddFilesToChain_no_muls_no_eloss_spoiled(track);
	else
		return NULL;

	TH2D *dp_vs_theta = new TH2D("dp_vs_theta", "", 35, 2.0, 40.0, 100, -0.1, 0.1);
	track->Project("dp_vs_theta", "(track.pfit.Mag()-track.pthrown.Mag())/track.pthrown.Mag():track.pthrown.Theta()*57.3");

	dp_vs_theta->FitSlicesY();
	
	TH1D *dp_vs_theta_2 = (TH1D*)gROOT->FindObject("dp_vs_theta_1");
	TH1D *h = dp_vs_theta_2->Clone(hname);
	h->SetLineColor(color);
	h->SetMarkerColor(color);
	h->SetLineWidth(2);
	
	return h;
}