#define BCAL_shower_info_cxx
// The class definition in BCAL_shower_info.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("BCAL_shower_info.C")
// Root > T->Process("BCAL_shower_info.C","some options")
// Root > T->Process("BCAL_shower_info.C+")
//

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

//const double energies[BCAL_shower_info::n_energies] = {0.02, 0.03, 0.04, 0.05, 0.07, 0.1, 0.150, 0.2, 0.5, 1.0, 2.0};
const double energies[BCAL_shower_info::n_energies] = {0.03, 0.05, 0.07, 0.1, 0.150, 0.2, 0.5, 1.0, 2.0};

const double inner_rad = 64.3;

const double bcal_length = 390.0;
const double bcal_center = 212.0;
const double target_z = 65.0;

const double z_min = bcal_center - bcal_length/2.0;
const double z_max = bcal_center + bcal_length/2.0;

void BCAL_shower_info::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();

  //initialize interesting histos first so they will show up first in TBrowser

  //this set of constants if you want to see all energies
  //const int n_E_hist_bins = 200;
  //const double E_hist_min = 0.0;
  //const double E_hist_max = 2.05;

  //this set of constants if you only want to see 50 MeV and higher.
  //const int n_E_hist_bins = 200;
  const int n_E_hist_bins = 100;
  const double E_hist_min = 0.000;
  const double E_hist_max = 2.05;

  n_showers_mean_hist = new TH2F("n_showers_mean_hist","n_showers_mean_hist",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);

  gteq_one_hist = new TH2F("gteq_one_hist","Efficiency (>= 1 shower reconstructed)",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);
  gteq_one_hist->GetXaxis()->SetTitle("z (cm)");
  gteq_one_hist->GetYaxis()->SetTitle("E (GeV)");
  gteq_one_hist->SetStats(0);

  excess_showers_hist = new TH2F("excess_showers_hist","Excess (>1) showers",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);
  excess_showers_hist->GetXaxis()->SetTitle("z (cm)");
  excess_showers_hist->GetYaxis()->SetTitle("E (GeV)");
  excess_showers_hist->SetStats(0);
//   excess_showers_hist->SetLabelSize(0.05,"y");
//   excess_showers_hist->GetYaxis()->SetTitleSize(0.05);
//   excess_showers_hist->SetLabelSize(0.05,"x");
//   excess_showers_hist->GetXaxis()->SetTitleSize(0.05);
//   excess_showers_hist->SetTitleSize(0.5);

  gteq_one_no_conversions_hist = new TH2F("gteq_one_no_conversions_hist","Efficiency (>= 1 shower reconstructed)",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);
  gteq_one_no_conversions_hist->GetXaxis()->SetTitle("z (cm)");
  gteq_one_no_conversions_hist->GetYaxis()->SetTitle("E (GeV)");
  gteq_one_no_conversions_hist->SetStats(0);

  excess_showers_no_conversions_hist = new TH2F("excess_showers_no_conversions_hist","Excess (>1) showers",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);
  excess_showers_no_conversions_hist->GetXaxis()->SetTitle("z (cm)");
  excess_showers_no_conversions_hist->GetYaxis()->SetTitle("E (GeV)");
  excess_showers_no_conversions_hist->SetStats(0);

  conversions_hist = new TH2F("conversions_hist","conversions_hist",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);

  match_hist = new TH2F("match_hist","match_hist",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);
  match_hist->SetStats(0);
  match_hist->SetTitle("BCAL single-shower reconstruction efficiency");
  match_hist->GetXaxis()->SetTitle("z (cm)");
  match_hist->GetYaxis()->SetTitle("E (GeV)");
  match_no_conversions_hist = new TH2F("match_no_conversions_hist","match_no_conversions_hist",n_z_bins,z_min,z_max,n_E_hist_bins,E_hist_min,E_hist_max);

  for (int i=0; i<n_energies; i++) {
    for (int j=0; j<n_z_bins; j++) {
      stringstream histo_name;
      histo_name << "n_showers_hist_" << i << "_" << j;
      n_showers_hist[i][j] = new TH1I(histo_name.str().c_str(),histo_name.str().c_str(),7,0,7);

      histo_name.str("");
      histo_name.clear();
      histo_name << "z_error_hist_" << i << "_" << j;
      z_error_hist[i][j] = new TH1F(histo_name.str().c_str(),histo_name.str().c_str(),200,-60,60);

      histo_name.str("");
      histo_name.clear();
      histo_name << "E_match_hist_" << i << "_" << j;
      double lower_limit = 0;
      double upper_limit = 2.0*energies[i];
      E_match_hist[i][j] = new TH1F(histo_name.str().c_str(),histo_name.str().c_str(),200,lower_limit,upper_limit);

      histo_name.str("");
      histo_name.clear();
      histo_name << "ratio_raw_thrown_" << i << "_" << j;
      ratio_raw_thrown[i][j] = new TH1F(histo_name.str().c_str(),histo_name.str().c_str(),50,0.0,2.0);

      count_eq_one[i][j]=0;
      count_gteq_one[i][j]=0;
      count_eq_one_no_conversions[i][j]=0;
      count_gteq_one_no_conversions[i][j]=0;
      count_conversions[i][j]=0;
      count_total[i][j]=0;
      count_match[i][j]=0;
      count_match_no_conversions[i][j]=0;
    }
  }
}

void BCAL_shower_info::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 BCAL_shower_info::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 BCAL_shower_info::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.

  GetEntry(entry);

  int E_bin = -1;

  for (int i=0; i<n_energies; i++) {
    double tol = 1e-5;
    if (fabs(E_thrown-energies[i])<tol) E_bin=i;
  }

  if (E_bin==-1) return kTRUE;

  int z_bin = n_z_bins*(z_entry_thrown-z_min)/bcal_length;
  if (z_bin < 0 || z_bin >= n_z_bins) return kTRUE;

  n_showers_hist[E_bin][z_bin]->Fill(nShowers);

  if (match) {
    z_error_hist[E_bin][z_bin]->Fill(matched_z_entry-z_entry_thrown);
    E_match_hist[E_bin][z_bin]->Fill(matched_E);
    ratio_raw_thrown[E_bin][z_bin]->Fill(matched_E_raw/E_thrown);
  }

  count_total[E_bin][z_bin]++;
  if (nShowers==1) count_eq_one[E_bin][z_bin]++;
  if (nShowers>=1) count_gteq_one[E_bin][z_bin]++;
  if (nShowers==1 && no_conversions) count_eq_one_no_conversions[E_bin][z_bin]++;
  if (nShowers>=1 && no_conversions) count_gteq_one_no_conversions[E_bin][z_bin]++;
  if (!no_conversions) count_conversions[E_bin][z_bin]++;
  if (match) count_match[E_bin][z_bin]++;
  if (match && no_conversions) count_match_no_conversions[E_bin][z_bin]++;

  return kTRUE;
}

void BCAL_shower_info::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 BCAL_shower_info::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.

  //efficiency is #matches/#events per bin
  //this will be used for constructing an efficiency function
  double efficiency[n_energies][n_z_bins];

  for (int i=0; i<n_energies; i++) {
    for (int j=0; j<n_z_bins; j++) {
      double z_val = (j+0.5)*bcal_length/n_z_bins + z_min;
      double E_val = energies[i];

      double IU_mean = n_showers_hist[i][j]->GetMean();
      n_showers_mean_hist->Fill(z_val,E_val,IU_mean);

      double gteq_one_pct = (double)count_gteq_one[i][j]/count_total[i][j];
      gteq_one_hist->Fill(z_val,E_val,gteq_one_pct);

      int count_no_conversions = count_total[i][j] - count_conversions[i][j];
      double gteq_one_pct_no_conversions = (double)count_gteq_one_no_conversions[i][j]/count_no_conversions;
      gteq_one_no_conversions_hist->Fill(z_val,E_val,gteq_one_pct_no_conversions);

      double excess_pct = 1.0 - (double)count_eq_one[i][j]/count_gteq_one[i][j];
      excess_showers_hist->Fill(z_val,E_val,excess_pct);

      double excess_pct_no_conversions = 1.0 - (double)count_eq_one_no_conversions[i][j]/count_gteq_one_no_conversions[i][j];
      excess_showers_no_conversions_hist->Fill(z_val,E_val,excess_pct_no_conversions);

      double conversion_pct=(double)count_conversions[i][j]/count_total[i][j];
      conversions_hist->Fill(z_val,E_val,conversion_pct);

      double match_pct=(double)count_match[i][j]/count_total[i][j];
      efficiency[i][j] = match_pct;
      match_hist->Fill(z_val,E_val,match_pct);

      double match_no_conv_pct = (double)count_match_no_conversions[i][j]/count_no_conversions;
      match_no_conversions_hist->Fill(z_val,E_val,match_no_conv_pct);

    }
  }

  //print the efficiency array
  cout << "const int n_energies = " << n_energies << ";" << endl;
  cout << "const int n_z_bins = " << n_z_bins << ";" << endl;
  cout << endl << "const double efficiency[n_energies][n_z_bins] = {" << endl;
  for (int i=0; i<n_energies; i++) {
    cout << "{";
    for (int j=0; j<n_z_bins; j++) {
      cout << efficiency[i][j];
      if (j != n_z_bins-1) cout << ", ";
    }
    cout << "}";
    if (i != n_energies-1) cout << ",";
    cout << endl;
  }
  cout << "};" << endl;
  
  cout << endl << "const double z_low[n_z_bins+1] = {";
  for (int i=0; i <= n_z_bins; i++) {
    double z_val = i*bcal_length/n_z_bins + z_min;
    cout << z_val;
    if (i != n_z_bins) cout << ", ";
  }
  cout << "};" << endl;

  cout << "const double energies[n_energies] = {";
  for (int i=0; i<n_energies; i++) {
    cout << energies[i];
    if (i != n_energies-1) cout << ", ";
  }
  cout << "};" << endl;

  //below is energy calibration stuff
  /*
  TGraphErrors *ratio_vs_E[n_z_bins];
  TF1 *func1[n_z_bins];
  double scale[n_z_bins];
  double nonlin[n_z_bins];
  double z[n_z_bins];
  for (int i=0; i<n_z_bins; i++) {
    z[i] = (i+0.5)*(z_max-z_min)/n_z_bins + z_min;

    double ratio[n_energies];
    double zeroes[n_energies];
    double errors[n_energies];
    for (int j=0; j<n_energies; j++) {
      ratio[j] = ratio_raw_thrown[j][i]->GetBinCenter(ratio_raw_thrown[j][i]->GetMaximumBin());
      //Take the error on the mean as the appropriate error of the ratio here
      errors[j] = ratio_raw_thrown[j][i]->GetMeanError();

      zeroes[j]=0;
    }

    int exclusions = 4; //We must exclude the first two energies (20 MeV and 30 MeV) since there is a 30 MeV threshold in clusterization.
    stringstream graphname;
    graphname << "ratio_vs_E_" << i;
    ratio_vs_E[i] = new TGraphErrors(n_energies-exclusions, energies+exclusions, ratio+exclusions, zeroes+exclusions, errors+exclusions);
    ratio_vs_E[i]->SetName(graphname.str().c_str());
    ratio_vs_E[i]->SetMarkerStyle(20);

    stringstream funcname;
    funcname << "func1_" << i;
    func1[i] = new TF1(funcname.str().c_str(),"[0]*pow(x,[1])");
    ratio_vs_E[i]->Fit(func1[i]);
    scale[i] = func1[i]->GetParameter(0);
    nonlin[i] = func1[i]->GetParameter(1);

    ratio_vs_E[i]->Write();
  }
  TCanvas *scale_canvas = new TCanvas("scale_canvas");
  TGraph *scale_graph = new TGraph(n_z_bins,z,scale);
  scale_graph->SetMarkerStyle(20);
  scale_graph->Draw("AP");

  TCanvas *nonlin_canvas = new TCanvas("nonlin_canvas");
  TGraph *nonlin_graph = new TGraph(n_z_bins,z,nonlin);
  nonlin_graph->SetMarkerStyle(20);
  nonlin_graph->Draw("AP");*/
}