#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", 200, 0.0, 180.0, 200, 0.0, 7.0);
	axes->SetXTitle("Polar angle #theta (degrees)");
	axes->SetYTitle("Total momentum (GeV/c)");
	axes->SetStats(0);
	axes->SetFillColor(kWhite);
	
	// Efficiency
	eff_numerator = (TH2D*)axes->Clone("eff_numerator");
	eff_denominator = (TH2D*)axes->Clone("eff_denominator");
	eff_numerator_wb = (TH2D*)axes->Clone("eff_numerator_wb");
	eff_denominator_wb = (TH2D*)axes->Clone("eff_denominator_wb");
	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();
	momRes = new TH3D("momRes", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 400, -0.4, 0.4);
	momRes->SetFillColor(kWhite);
	trans_momRes = (TH3D*)momRes->Clone("trans_momRes");

	// Theta and phi res
	thetaRes = new TH3D("thetaRes", ""
		, xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax()
		, yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax()
		, 50, -50.0, 50.0);
	thetaRes->SetXTitle(xaxis->GetTitle());
	thetaRes->SetYTitle(yaxis->GetTitle());
	thetaRes->SetZTitle("Resolution (mrad)");
	thetaRes->SetFillColor(kWhite);
	phiRes = (TH3D*)thetaRes->Clone("phiRes");

	// Nhits
	Ncdchits = (TH2D*)axes->Clone("Ncdchits");
	Nfdchits = (TH2D*)axes->Clone("Nfdchits");
	Ncdchits->SetTitle("CDC Nhits per event");
	Nfdchits->SetTitle("FDC Nhits per event");

	Ncdchits_on_track = (TH2D*)axes->Clone("Ncdchits_on_track");
	Nfdchits_on_track = (TH2D*)axes->Clone("Nfdchits_on_track");
	Ncdchits_on_track_wb = (TH2D*)axes->Clone("Ncdchits_on_track_wb");
	Nfdchits_on_track_wb = (TH2D*)axes->Clone("Nfdchits_on_track_wb");

	Ncdchits_on_track->SetTitle("CDC Nhits per track (time-based)");
	Nfdchits_on_track->SetTitle("FDC Nhits per track (time-based)");
	Ncdchits_on_track_wb->SetTitle("CDC Nhits per track (wire-based)");
	Nfdchits_on_track_wb->SetTitle("FDC Nhits per track (wire-based)");

	// Chisq
	chisq_per_Ndof = new TH3D("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);
	chisq_per_Ndof->SetFillColor(kWhite);
	chisq_per_Ndof_wb = (TH3D*)chisq_per_Ndof->Clone("chisq_per_Ndof_wb");
}

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
	GetEntry(entry);
	Nevents_processed++;

	// Older versions of ROOT used the TVector3 objects as members of the
	// selector class. Now it store the separate components forcing us
	// to create the vectors as local variables (yechh!).
	TVector3 pthrown(pthrown_fX, pthrown_fY, pthrown_fZ);
	TVector3 can_p(can_p_fX, can_p_fY, can_p_fZ);
	TVector3 trk_p(trk_p_fX, trk_p_fY, trk_p_fZ);
	TVector3 part_p(part_p_fX, part_p_fY, part_p_fZ);

	// Useful constants
	bool good_event = part_trk_Ndof>0 ? (part_trk_chisq/part_trk_Ndof < 10000.0):kFALSE;
	bool good_event_wb = trk_trk_Ndof>0 ? (trk_trk_chisq/trk_trk_Ndof < 10000.0):kFALSE;
	double p = pthrown.Mag();
	double pt = pthrown.Pt();
	double theta = pthrown.Theta()*TMath::RadToDeg();
	double phi = pthrown.Phi();
	double p_fit = part_p.Mag();
	double pt_fit = part_p.Pt();
	double theta_fit = part_p.Theta()*TMath::RadToDeg();
	double phi_fit = part_p.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);
	eff_denominator_wb->Fill(theta, p);
	prob_vs_theta->Fill(theta, TMath::Prob(part_trk_chisq, part_trk_Ndof));
	prob_vs_theta_wb->Fill(theta, TMath::Prob(trk_trk_chisq, trk_trk_Ndof));

	// MomRes, ThetaRes, PhiRes
	if(good_event){
		if(p!=0.0)momRes->Fill(theta, p, (p_fit-p)/p);
		if(pt!=0.0)trans_momRes->Fill(theta, p, (pt_fit-pt)/pt);
		thetaRes->Fill(theta, p, (theta_fit - theta)*1000.0/TMath::RadToDeg()); // in mrad
		phiRes->Fill(theta, p, (phi_fit - phi)*1000.0); // in mrad
	}
	
	// Nhits
	Ncdchits->Fill(theta, p, Ncdc);
	Nfdchits->Fill(theta, p, Nfdc);
	if(good_event)Ncdchits_on_track->Fill(theta, p, part_Ncdc);
	if(good_event)Nfdchits_on_track->Fill(theta, p, part_Nfdc);
	if(good_event_wb)Ncdchits_on_track_wb->Fill(theta, p, trk_Ncdc);
	if(good_event_wb)Nfdchits_on_track_wb->Fill(theta, p, trk_Nfdc);
	
	// Chisq
	if(part_trk_Ndof>0)chisq_per_Ndof->Fill(theta, p, part_trk_chisq/part_trk_Ndof);
	if(trk_trk_Ndof>0)chisq_per_Ndof_wb->Fill(theta, p, trk_trk_chisq/trk_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_wb);
	eff_vs_p_vs_theta_wb->SetTitle("Wire-based Tracking Efficiency vs. total momentum and #theta angle");

	// MomRes
	cout<<"    momRes:"<<endl;
	momRes->FitSlicesZ();
	TH2D *dp_over_p_amp = (TH2D*)gROOT->FindObject("momRes_0");
	TH2D *dp_over_p_mean = (TH2D*)gROOT->FindObject("momRes_1");
	TH2D *dp_over_p_sigma = (TH2D*)gROOT->FindObject("momRes_2");
	TH2D *dp_over_p_chi2 = (TH2D*)gROOT->FindObject("momRes_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");

	cout<<"    trans_momRes:"<<endl;
	trans_momRes->FitSlicesZ();
	TH2D *dpt_over_pt_amp = (TH2D*)gROOT->FindObject("trans_momRes_0");
	TH2D *dpt_over_pt_mean = (TH2D*)gROOT->FindObject("trans_momRes_1");
	TH2D *dpt_over_pt_sigma = (TH2D*)gROOT->FindObject("trans_momRes_2");
	TH2D *dpt_over_pt_chi2 = (TH2D*)gROOT->FindObject("trans_momRes_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");

	// Theta and phi res
	cout<<"    thetaRes:"<<endl;
	thetaRes->FitSlicesZ();
	TH2D *dtheta_amp = (TH2D*)gROOT->FindObject("thetaRes_0");
	TH2D *dtheta_mean = (TH2D*)gROOT->FindObject("thetaRes_1");
	TH2D *dtheta_sigma = (TH2D*)gROOT->FindObject("thetaRes_2");
	TH2D *dtheta_chi2 = (TH2D*)gROOT->FindObject("thetaRes_chi2");
	dtheta_amp->SetName("dtheta_amp");
	dtheta_mean->SetName("dtheta_mean");
	dtheta_sigma->SetName("dtheta_sigma");
	dtheta_chi2->SetName("dtheta_chi2");

	// PhiRes
	cout<<"    phiRes:"<<endl;
	phiRes->FitSlicesZ();
	TH2D *dphi_amp = (TH2D*)gROOT->FindObject("phiRes_0");
	TH2D *dphi_mean = (TH2D*)gROOT->FindObject("phiRes_1");
	TH2D *dphi_sigma = (TH2D*)gROOT->FindObject("phiRes_2");
	TH2D *dphi_chi2 = (TH2D*)gROOT->FindObject("phiRes_chi2");
	dphi_amp->SetName("dphi_amp");
	dphi_mean->SetName("dphi_mean");
	dphi_sigma->SetName("dphi_sigma");
	dphi_chi2->SetName("dphi_chi2");

	// Nhits
	cout<<"    Nhits:"<<endl;
	Ncdchits->Divide(eff_denominator);
	Nfdchits->Divide(eff_denominator);
	Ncdchits_on_track->Divide(eff_numerator);
	Nfdchits_on_track->Divide(eff_numerator);
	Ncdchits_on_track_wb->Divide(eff_numerator_wb);
	Nfdchits_on_track_wb->Divide(eff_numerator_wb);

}