#define hd_res_charged_selector_cxx
// The class definition in hd_res_charged_selector.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("hd_res_charged_selector.C")
// Root > T->Process("hd_res_charged_selector.C","some options")
// Root > T->Process("hd_res_charged_selector.C+")
//

#include "hd_res_charged_selector.h"
#include <TH2.h>
#include <TStyle.h>


void hd_res_charged_selector::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();
	
	Nevents_processed = 0;

	// Create a 2D histo with the binning and limits that will be
	// used for all histograms. We want all resolutions and efficiencies
	// to be a function of total momentum and theta angle
	TH2D *axes = new TH2D("axes","nada", 300, 0.0, 180.0, 300, 0.0, 7.0);
	axes->SetXTitle("Polar angle #theta (degrees)");
	axes->SetYTitle("Total momentum (GeV/c)");
	axes->SetStats(0);
	
	// Standard cut
	TCut cut1("trk_chisq/trk_Ndof<5.0");

	// Efficiency
	eff_numerator = (TH2D*)axes->Clone("eff_numerator");
	eff_numerator_wb = (TH2D*)axes->Clone("eff_numerator_wb");
	eff_denominator = (TH2D*)axes->Clone("eff_denominator");
	eff_denominator_nhits_cut = (TH2D*)axes->Clone("eff_denominator_nhits_cut");
	prob_vs_theta = new TH2D("prob_vs_theta","Probability vs. #theta", 180, 0.0, 180.0, 1000, 0.0, 1.0);
	prob_vs_theta_wb = (TH2D*)prob_vs_theta->Clone("prob_vs_theta_wb");

	// MomRes
	TAxis *xaxis = axes->GetXaxis();
	TAxis *yaxis = axes->GetYaxis();
	momRes3D = new TH3D("momRes3D", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 400, -0.4, 0.4);
	trans_momRes3D = (TH3D*)momRes3D->Clone("trans_momRes3D");

	// ThetaRes
	thetaRes3D = new TH3D("thetaRes3D", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 50, -50.0, 50.0);

	// PhiRes
	phiRes3D = new TH3D("phiRes3D", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 50, -50.0, 50.0);

	// Nhits
	Ncdchits = (TH2D*)axes->Clone("Ncdchits");
	Nfdchits = (TH2D*)axes->Clone("Nfdchits");
	Ncdchits->SetTitle("CDC Nhits per track vs. total momentum and #theta angle");
	Nfdchits->SetTitle("FDC Nhits per track vs. total momentum and #theta angle");

	// Chisq
	trk_chisq_per_Ndof = new TH3D("trk_chisq_per_Ndof", "Tracking #chi^2/Ndof vs. total momentum and #theta angle"
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 150, 0.0, 15.0);
}

void hd_res_charged_selector::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

}

Bool_t hd_res_charged_selector::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // It can be passed to either hd_res_charged_selector::GetEntry() or TBranch::GetEntry()
   // to read either all or the required parts of the data. When processing
   // keyed objects with PROOF, the object is already loaded and is available
   // via the fObject pointer.
   //
   // This function should contain the "body" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

	// Read this track into memory
	fChain->GetTree()->GetEntry(entry);
	Nevents_processed++;

	// Useful constants
	bool good_event = trk_Ndof>0 ? (trk_chisq/trk_Ndof < 5.0):kFALSE;
	bool good_event_wb = trk_Ndof_wb>0 ? (trk_chisq_wb/trk_Ndof_wb < 5.0):kFALSE;
	double p = pthrown.Mag();
	double pt = pthrown.Pt();
	double theta = pthrown.Theta()*TMath::RadToDeg();
	double phi = pthrown.Phi();
	double p_fit = pfit.Mag();
	double pt_fit = pfit.Pt();
	double theta_fit = pfit.Theta()*TMath::RadToDeg();
	double phi_fit = pfit.Phi();
	
	// Efficiency
	if(good_event)eff_numerator->Fill(theta, p);
	if(good_event_wb)eff_numerator_wb->Fill(theta, p);
	eff_denominator->Fill(theta, p);
	if((Ncdc+Nfdc)>0)eff_denominator_nhits_cut->Fill(theta, p);
	prob_vs_theta->Fill(theta, TMath::Prob(trk_chisq, trk_Ndof));
	prob_vs_theta_wb->Fill(theta, TMath::Prob(trk_chisq_wb, trk_Ndof_wb));

	// MomRes, ThetaRes, PhiRes
	if(good_event){
		momRes3D->Fill(theta, p, (p_fit-p)/p);
		trans_momRes3D->Fill(theta, p, (pt_fit-pt)/pt);
		thetaRes3D->Fill(theta, p, (theta_fit - theta)*1000.0/TMath::RadToDeg()); // in mrad
		phiRes3D->Fill(theta, p, (phi_fit - phi)*1000.0); // in mrad
	}
	
	// Nhits
	if(good_event){
		Ncdchits->Fill(theta, p, Ncdc);
		Nfdchits->Fill(theta, p, Nfdc);
	}
	
	// Chisq
	if(trk_Ndof>0)trk_chisq_per_Ndof->Fill(theta, p, trk_chisq/trk_Ndof);

	// Print ticker
	if(Nevents_processed%2000 == 0){
		cout<<"\r  "<<Nevents_processed/1000<<"k events processed   ";
		cout.flush();
	}

   return kTRUE;
}

void hd_res_charged_selector::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

}

void hd_res_charged_selector::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

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

	// Efficiency
	cout<<"    Efficiency:"<<endl;
	eff_vs_p_vs_theta = (TH2D*)eff_numerator->Clone("eff_vs_p_vs_theta");
	eff_vs_p_vs_theta->Divide(eff_denominator);
	eff_vs_p_vs_theta->SetTitle("Tracking Efficiency vs. total momentum and #theta angle");
	eff_vs_p_vs_theta_wb = (TH2D*)eff_numerator_wb->Clone("eff_vs_p_vs_theta_wb");
	eff_vs_p_vs_theta_wb->Divide(eff_denominator);
	eff_vs_p_vs_theta_wb->SetTitle("Wire-based Tracking Efficiency vs. total momentum and #theta angle");
	eff_vs_p_vs_theta_wb_nhits_cut = (TH2D*)eff_numerator_wb->Clone("eff_vs_p_vs_theta_wb_nhits_cut");
	eff_vs_p_vs_theta_wb_nhits_cut->Divide(eff_denominator_nhits_cut);
	eff_vs_p_vs_theta_wb_nhits_cut->SetTitle("Wire-based Tracking Efficiency vs. total momentum and #theta angle (Nhits>0)");

	// MomRes
	cout<<"    momRes3D:"<<endl;
	momRes3D->FitSlicesZ();
	TH2D *dp_over_p_amp = (TH2D*)gROOT->FindObject("momRes3D_0");
	TH2D *dp_over_p_mean = (TH2D*)gROOT->FindObject("momRes3D_1");
	TH2D *dp_over_p_sigma = (TH2D*)gROOT->FindObject("momRes3D_2");
	TH2D *dp_over_p_chi2 = (TH2D*)gROOT->FindObject("momRes3D_chi2");
	dp_over_p_amp->SetName("dp_over_p_amp");
	dp_over_p_mean->SetName("dp_over_p_mean");
	dp_over_p_sigma->SetName("dp_over_p_sigma");
	dp_over_p_chi2->SetName("dp_over_p_chi2");
	//delete momRes3D;

	cout<<"    trans_momRes3D:"<<endl;
	trans_momRes3D->FitSlicesZ();
	TH2D *dpt_over_pt_amp = (TH2D*)gROOT->FindObject("trans_momRes3D_0");
	TH2D *dpt_over_pt_mean = (TH2D*)gROOT->FindObject("trans_momRes3D_1");
	TH2D *dpt_over_pt_sigma = (TH2D*)gROOT->FindObject("trans_momRes3D_2");
	TH2D *dpt_over_pt_chi2 = (TH2D*)gROOT->FindObject("trans_momRes3D_chi2");
	dpt_over_pt_amp->SetName("dpt_over_pt_amp");
	dpt_over_pt_mean->SetName("dpt_over_pt_mean");
	dpt_over_pt_sigma->SetName("dpt_over_pt_sigma");
	dpt_over_pt_chi2->SetName("dpt_over_pt_chi2");
	//delete trans_momRes3D;

	// ThetaRes
	cout<<"    thetaRes3D:"<<endl;
	thetaRes3D->FitSlicesZ();
	TH2D *dtheta_amp = (TH2D*)gROOT->FindObject("thetaRes3D_0");
	TH2D *dtheta_mean = (TH2D*)gROOT->FindObject("thetaRes3D_1");
	TH2D *dtheta_sigma = (TH2D*)gROOT->FindObject("thetaRes3D_2");
	TH2D *dtheta_chi2 = (TH2D*)gROOT->FindObject("thetaRes3D_chi2");
	dtheta_amp->SetName("dtheta_amp");
	dtheta_mean->SetName("dtheta_mean");
	dtheta_sigma->SetName("dtheta_sigma");
	dtheta_chi2->SetName("dtheta_chi2");
	delete thetaRes3D;

	// PhiRes
	cout<<"    phiRes3D:"<<endl;
	phiRes3D->FitSlicesZ();
	TH2D *dphi_amp = (TH2D*)gROOT->FindObject("phiRes3D_0");
	TH2D *dphi_mean = (TH2D*)gROOT->FindObject("phiRes3D_1");
	TH2D *dphi_sigma = (TH2D*)gROOT->FindObject("phiRes3D_2");
	TH2D *dphi_chi2 = (TH2D*)gROOT->FindObject("phiRes3D_chi2");
	dphi_amp->SetName("dphi_amp");
	dphi_mean->SetName("dphi_mean");
	dphi_sigma->SetName("dphi_sigma");
	dphi_chi2->SetName("dphi_chi2");
	delete phiRes3D;

	// Nhits
	cout<<"    Nhits:"<<endl;
	Ncdchits->Divide(eff_numerator);
	Nfdchits->Divide(eff_numerator);

}