// This macro will make a ROOT file containing histograms
// of the charged tracking resolutions suitable for use
// with the src/programs/Simulation/plugins/hdparsim
// plugin for doing parameteric simulations

#include "filelist.h"

//--------------------
// mk_hd_charged_res
//--------------------
void mk_hd_res_charged(void)
{
	gROOT->Reset();
	gStyle->SetPalette(1);

	// Create chain of all the input trees
	TChain *trkeff = new TChain("TRACKING/trkeff","trkeff");
	AddFilesToChain(trkeff);
	
	// Open output file to hold histograms
	TFile *f = new TFile("hd_res_charged.root", "RECREATE");
	
	// All of the real work is done in the TSelector derived
	// class in hd_res_charged_selector.cc
	trkeff->Process("hd_res_charged_selector.C");

	// Close resolutions file
	f->Write();
}

#if 0
//--------------------
// MakeEfficiency
//--------------------
void MakeEfficiency(TChain *trkeff, TH2D *axes, TCut &cut1)
{
	TH2D *h = axes->Clone("eff_vs_p_vs_theta");
	h->SetTitle("Tracking Efficiency vs. total momentum and #theta angle");
	cout<<"Making "<<h->GetName()<<endl;

	TH2D *numerator = axes->Clone("eff_numerator");
	TH2D *denominator = axes->Clone("eff_denominator");
	
	cout<<"   Projecting numerator ..."<<endl;
	trkeff->Project("eff_numerator","E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()", cut1);
	cout<<"   Projecting denominator ..."<<endl;
	trkeff->Project("eff_denominator","E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()");

	cout<<"   Finalizing ..."<<endl;
	h->Divide(numerator, denominator);
	h->GetZaxis()->SetRangeUser(0.2, 1.00);
	h->Draw("colz");	
}

//--------------------
// MakeMomRes
//--------------------
void MakeMomRes(TChain *trkeff, TH2D *axes, TCut &cut1)
{
	TH2D *h = axes->Clone("dp_over_p_vs_p_vs_theta");
	h->SetTitle("Relative Momentum resolution vs. total momentum and #theta angle");
	cout<<"Making "<<h->GetName()<<endl;
	
	// Make 3D histo to project onto
	TAxis *xaxis = axes->GetXaxis();
	TAxis *yaxis = axes->GetYaxis();
	TH3D *tmp3d = new TH3D("tmp3d", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 1000, -1.0, 1.0);
	
	cout<<"   Projecting  ..."<<endl;
	trkeff->Project("tmp3d","(E.pfit.Mag()-E.pthrown.Mag())/E.pthrown.Mag():E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()", cut1);
	cout<<"   Fitting ..."<<endl;
	tmp3d->FitSlicesZ();

	cout<<"   Finalizing ..."<<endl;

	// Find histos resulting from fits
	TH2D *tmp3d_0 = (TH2D*)gROOT->FindObject("tmp3d_0");
	TH2D *tmp3d_1 = (TH2D*)gROOT->FindObject("tmp3d_1");
	TH2D *tmp3d_2 = (TH2D*)gROOT->FindObject("tmp3d_2");
	TH2D *tmp3d_chi2 = (TH2D*)gROOT->FindObject("tmp3d_chi2");
	
	// Copy the contents of sigma histo into h
	h->Add(tmp3d_2);
	h->GetZaxis()->SetRangeUser(0.0, 0.15);
	h->Draw("colz");
	
	// Clean up
	delete tmp3d_0;
	delete tmp3d_1;
	delete tmp3d_2;
	delete tmp3d_chi2;
	delete tmp3d;
}

//--------------------
// MakeTransMomRes
//--------------------
void MakeTransMomRes(TChain *trkeff, TH2D *axes, TCut &cut1)
{
	TH2D *h = axes->Clone("dpt_over_pt_vs_p_vs_theta");
	h->SetTitle("Relative Transverse Momentum resolution vs. total momentum and #theta angle");
	cout<<"Making "<<h->GetName()<<endl;
	
	// Make 3D histo to project onto
	TAxis *xaxis = axes->GetXaxis();
	TAxis *yaxis = axes->GetYaxis();
	TH3D *tmp3d = new TH3D("tmp3d", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 1000, -1.0, 1.0);
	
	cout<<"   Projecting  ..."<<endl;
	trkeff->Project("tmp3d","(E.pfit.Pt()-E.pthrown.Pt())/E.pthrown.Pt():E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()", cut1);
	cout<<"   Fitting ..."<<endl;
	tmp3d->FitSlicesZ();

	cout<<"   Finalizing ..."<<endl;

	// Find histos resulting from fits
	TH2D *tmp3d_0 = (TH2D*)gROOT->FindObject("tmp3d_0");
	TH2D *tmp3d_1 = (TH2D*)gROOT->FindObject("tmp3d_1");
	TH2D *tmp3d_2 = (TH2D*)gROOT->FindObject("tmp3d_2");
	TH2D *tmp3d_chi2 = (TH2D*)gROOT->FindObject("tmp3d_chi2");
	
	// Copy the contents of sigma histo into h
	h->Add(tmp3d_2);
	h->GetZaxis()->SetRangeUser(0.0, 0.15);
	h->Draw("colz");
	
	// Clean up
	delete tmp3d_0;
	delete tmp3d_1;
	delete tmp3d_2;
	delete tmp3d_chi2;
	delete tmp3d;
}

//--------------------
// MakeThetaRes
//--------------------
void MakeThetaRes(TChain *trkeff, TH2D *axes, TCut &cut1)
{
	TH2D *h = axes->Clone("dtheta_vs_p_vs_theta");
	h->SetTitle("Polar angular resolution vs. total momentum and #theta angle");
	cout<<"Making "<<h->GetName()<<endl;
	
	// Make 3D histo to project onto
	TAxis *xaxis = axes->GetXaxis();
	TAxis *yaxis = axes->GetYaxis();
	TH3D *tmp3d = new TH3D("tmp3d", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 50, -50.0, 50.0);
	
	cout<<"   Projecting  ..."<<endl;
	trkeff->Project("tmp3d","(E.pfit.Theta()-E.pthrown.Theta())*1000.0:E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()", cut1);
	cout<<"   Fitting ..."<<endl;
	tmp3d->FitSlicesZ();

	cout<<"   Finalizing ..."<<endl;

	// Find histos resulting from fits
	TH2D *tmp3d_0 = (TH2D*)gROOT->FindObject("tmp3d_0");
	TH2D *tmp3d_1 = (TH2D*)gROOT->FindObject("tmp3d_1");
	TH2D *tmp3d_2 = (TH2D*)gROOT->FindObject("tmp3d_2");
	TH2D *tmp3d_chi2 = (TH2D*)gROOT->FindObject("tmp3d_chi2");
	
	// Copy the contents of sigma histo into h
	h->Add(tmp3d_2);
	h->GetZaxis()->SetRangeUser(0.0, 10.0);
	h->Draw("colz");
	
	// Clean up
	delete tmp3d_0;
	delete tmp3d_1;
	delete tmp3d_2;
	delete tmp3d_chi2;
	delete tmp3d;
}

//--------------------
// MakePhiRes
//--------------------
void MakePhiRes(TChain *trkeff, TH2D *axes, TCut &cut1)
{
	TH2D *h = axes->Clone("dphi_vs_p_vs_theta");
	h->SetTitle("Azimuthal angular resolution vs. total momentum and #theta angle");
	cout<<"Making "<<h->GetName()<<endl;
	
	// Make 3D histo to project onto
	TAxis *xaxis = axes->GetXaxis();
	TAxis *yaxis = axes->GetYaxis();
	TH3D *tmp3d = new TH3D("tmp3d", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 50, -50.0, 50.0);
	
	cout<<"   Projecting  ..."<<endl;
	trkeff->Project("tmp3d","(E.pfit.Phi()-E.pthrown.Phi())*1000.0:E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()", cut1);
	cout<<"   Fitting ..."<<endl;
	tmp3d->FitSlicesZ();

	cout<<"   Finalizing ..."<<endl;

	// Find histos resulting from fits
	TH2D *tmp3d_0 = (TH2D*)gROOT->FindObject("tmp3d_0");
	TH2D *tmp3d_1 = (TH2D*)gROOT->FindObject("tmp3d_1");
	TH2D *tmp3d_2 = (TH2D*)gROOT->FindObject("tmp3d_2");
	TH2D *tmp3d_chi2 = (TH2D*)gROOT->FindObject("tmp3d_chi2");
	
	// Copy the contents of sigma histo into h
	h->Add(tmp3d_2);
	h->GetZaxis()->SetRangeUser(0.0, 10.0);
	h->Draw("colz");
	//return;	

	// Clean up
	delete tmp3d_0;
	delete tmp3d_1;
	delete tmp3d_2;
	delete tmp3d_chi2;
	delete tmp3d;
}

//--------------------
// MakeTrkChisqNdof
//--------------------
void MakeTrkChisqNdof(TChain *trkeff, TH2D *axes, TCut &cut1)
{
	
	// Make 3D histo to project onto
	TAxis *xaxis = axes->GetXaxis();
	TAxis *yaxis = axes->GetYaxis();
	TH3D *h = new TH3D("trk_chisq_per_Ndof_vs_p_vs_theta", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 150, 0.0, 15.0);
	h->SetTitle("Tracking #chi^2/Ndof vs. total momentum and #theta angle");
	cout<<"Making "<<h->GetName()<<endl;
	
	cout<<"   Projecting  ..."<<endl;
	trkeff->Project(h->GetName(),"E.trk_chisq/E.trk_Ndof:E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()");

}

//--------------------
// MakeNhits
//--------------------
void MakeNhits(TChain *trkeff, TH2D *axes, TCut &cut1)
{
	TH2D *Ncdchits_vs_p_vs_theta = axes->Clone("Ncdchits_vs_p_vs_theta");
	TH2D *Nfdchits_vs_p_vs_theta = axes->Clone("Nfdchits_vs_p_vs_theta");
	Ncdchits_vs_p_vs_theta->SetTitle("CDC Nhits per track vs. total momentum and #theta angle");
	Nfdchits_vs_p_vs_theta->SetTitle("FDC Nhits per track vs. total momentum and #theta angle");

	cout<<"Making Nhits histos"<<endl;
	
	cout<<"   Projecting  Ncdchits_vs_p_vs_theta ..."<<endl;
	trkeff->Project("Ncdchits_vs_p_vs_theta","E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()", cut1*"E.Ncdc");

	cout<<"   Projecting  Nfdchits_vs_p_vs_theta ..."<<endl;
	trkeff->Project("Nfdchits_vs_p_vs_theta","E.pthrown.Mag():E.pthrown.Theta()*TMath::RadToDeg()", cut1*"E.Nfdc");

	cout<<"Finalizing ..."<<endl;
	TH2D *eff_numerator = gROOT->FindObject("eff_numerator");
	Ncdchits_vs_p_vs_theta->Divide(eff_numerator);
	Nfdchits_vs_p_vs_theta->Divide(eff_numerator);
}
#endif