// purpose: project 2d hist to y axis and find and fit peak position
//          hist2d_4: energy gen vs. energy rec BCAL  bins 7-55
//          hist2d_5: energy gen vs. energy rec FCAL  bins 7-55
//          hist2d_6: angle gen vs. angle rec BCAL
//          hist2d_7: angle gen vs. angle rec FCAL

using namespace std;

#include "TH1F.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH2F.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TF1.h"
#include "TGraphErrors.h"

#define DEBUG 0

Int_t pandf(Double_t *, Double_t *, Double_t *, Double_t *, TH2D **);
Int_t getmax(TH1D* , Double_t*);


Int_t pandf(Double_t *fresBE, Double_t *fresFE, Double_t *fresBA, Double_t *fresFA, TH2D **h) {
  
  TH1D *hbcalE[64];
  TH1D *hfcalE[64];
  TH1D *hbcalA[240];
  TH1D *hfcalA[240];
  
  
  for (Int_t i=0;i<64;i++){
    char str1[128];
    sprintf(str1,"BcalE_proj%d",i);
    hbcalE[i] = h[0]->ProjectionX(str1,i+1,i+1,"o");
    sprintf(str1,"FcalE_proj%d",i);
    hfcalE[i] = h[1]->ProjectionX(str1,i+1,i+1,"o");
  }
  for (Int_t i=0;i<240;i++){
    char str1[128];
    sprintf(str1,"BcalA_proj%d",i);
    hbcalA[i] = h[2]->ProjectionX(str1,i+1,i+1,"o");
    sprintf(str1,"FcalA_proj%d",i);
    hfcalA[i] = h[3]->ProjectionX(str1,i+1,i+1,"o");
  }

  Double_t *p;
  for (Int_t i=0;i<64;i++){
    p =  fresBE + i*4;
    if (hbcalE[i]->GetEntries()>700){
      getmax(hbcalE[i], p);
    } else{
      p[0] = 0.;
      p[1] = 0.;
      p[2] = 0.;
      p[3] = 0.;
    }
    p =  fresFE + i*4;
    if (hfcalE[i]->GetEntries()>700){
      getmax(hfcalE[i], p);
    } else{
      p[0] = 0.;
      p[1] = 0.;
      p[2] = 0.;
      p[3] = 0.;
    }
  }
  
  for (Int_t i=0;i<240;i++){
    p = fresBA + i*4;
    if (hbcalA[i]->GetEntries()>500){
      getmax(hbcalA[i], p);
    } else{
      p[0] = 0.;
      p[1] = 0.;
      p[2] = 0.;
      p[3] = 0.;
    }
    p = fresFA + i*4;
    if (hfcalA[i]->GetEntries()>500){
      getmax(hfcalA[i], p);
    } else{
      p[0] = 0.;
      p[1] = 0.;
      p[2] = 0.;
      p[3] = 0.;
    }
  }
  
  return 0;
}

Int_t getmax(TH1D* h, Double_t *fres){
  
  Int_t bin = h->GetMaximumBin();
  Double_t max = h->GetBinLowEdge(bin);
  Double_t min = max * 0.7;
  max = max * 1.2;
  if (max>h->GetBinLowEdge(h->GetNbinsX())){
    max = h->GetBinLowEdge(h->GetNbinsX());
  }
  h->Fit("gaus","Q","R",min,max);
  TF1 *gf = h->GetFunction("gaus");
  if (!gf){
    fres[0] = 0;
    fres[1] = 0;
    fres[3] = 0;
    fres[4] = 0;
    return 0;
  }
  
  Double_t pos = gf->GetParameter(1);
  Double_t sig = gf->GetParameter(2);
  
  max = pos + 2.*sig;
  min = pos - 2.*sig;
  h->Fit("gaus","Q","R",min,max);
  gf = h->GetFunction("gaus");
  pos = gf->GetParameter(1);
  sig = gf->GetParameter(2);
  
  fres[0] = pos;
  fres[1] = sig;
  fres[2] =  gf->GetParError(1);
  fres[3] =  gf->GetParError(2);
  
  return 1;
  
}

Int_t main(){

  Double_t fresBE[265];
  Double_t fresFE[265];
  Double_t fresBA[960];
  Double_t fresFA[960];


  TH2F *hBCAL_Gen = new TH2F("hBCAL_Gen","Energy vs Angle gen BCAL", 240, 5., 25., 64, 0., 3.2);
  hBCAL_Gen->GetXaxis()->SetTitle("Angle generated [deg]");
  hBCAL_Gen->GetYaxis()->SetTitle("Energy generated [GeV]");
  
  TH2F *hFCAL_Gen = new TH2F("hFCAL_Gen","Energy vs Angle gen FCAL", 240, 5., 25., 64, 0., 3.2);
  hFCAL_Gen->GetXaxis()->SetTitle("Angle generated [deg]");
  hFCAL_Gen->GetYaxis()->SetTitle("Energy generated [GeV]");
  
  TH2F *hBCAL_Rec = new TH2F("hBCAL_Rec","Energy vs Angle rec BCAL", 240, 5., 25., 64, 0., 3.2);
  hBCAL_Rec->GetXaxis()->SetTitle("Angle generated [deg]");
  hBCAL_Rec->GetYaxis()->SetTitle("Energy generated [GeV]");
  
  TH2F *hFCAL_Rec = new TH2F("hFCAL_Rec","Energy vs Angle rec FCAL", 240, 5., 25., 64, 0., 3.2);
  hFCAL_Rec->GetXaxis()->SetTitle("Angle generated [deg]");
  hFCAL_Rec->GetYaxis()->SetTitle("Energy generated [GeV]");
  

  Int_t EventNum;
  Int_t NGen;    // number generated photons
  Int_t NRecB;   // number of reconstructed photons in BCAL
  Int_t NRecF;   // number of reconstructed photons in FCAL
  Float_t GenE[20]; // energated energy
  Float_t GenA[20]; // generatged polar angle
  Float_t BcalE[20];  // reconstucted energy in BCAL
  Float_t BcalA[20];  // reconstruced polar angle in BCAL
  Float_t FcalE[20];  // reconstruced energy in FCAL
  Float_t FcalA[20];  // reconstruced polar angle in FCAL
  
  char inf[128];
  cout<<"File Name: ";
  cin>>inf;

  char outf[128];
  char temps[128];
  char *strp;
  strp = strstr(inf,"_");
  sprintf(temps,"%s",strp);
  sprintf(outf,"hists%s",temps);

  TFile *INF = new TFile(inf,"READ");
  if (!INF){
    cout<<"Error open file "<<inf<<endl;
    return 0;
  }
  TTree *EventTree = (TTree*)INF->Get("EventTree");

  TH2D *h[4];
  h[0] = (TH2D*)INF->Get("hist2d_4"); // BCAL Energy gen vs. Energy rec
  h[1] = (TH2D*)INF->Get("hist2d_5"); // FCAL Energy gen vs. Energy rec
  h[2] = (TH2D*)INF->Get("hist2d_6"); // BCAL Angle gen vs. Angle rec
  h[3] = (TH2D*)INF->Get("hist2d_7"); // FCAL Angle gen vs. Angle rec


  // determine for each generated Energy and Angle bin
  // the reconstructed mean Energy and Angle and width
  // by a gausian fit to the distribution as given by
  // a projection of the 4 histograms above to the horizontal
  // axis for each bin in the vertical axis (E/A generated).

  Int_t res1 = pandf(fresBE, fresFE, fresBA, fresFA, h);

  if (DEBUG){
    for (Int_t k=0;k<64;k++){
      for (Int_t j=0;j<4;j++)
	cout<<fresBE[k*4+j]<<"  ";
      
      cout<<"         ";
      for (Int_t j=0;j<4;j++)
	cout<<fresFE[k*4+j]<<"  ";
      
      cout<<endl;
    }
  }

  EventTree->SetBranchAddress("EventNum", &EventNum);
  EventTree->SetBranchAddress("NGen", &NGen);
  EventTree->SetBranchAddress("GenE",  GenE); 
  EventTree->SetBranchAddress("GenA",  GenA); 
  
  EventTree->SetBranchAddress("NRecB",    &NRecB); 
  EventTree->SetBranchAddress("BcalE",    BcalE); 
  EventTree->SetBranchAddress("BcalA",    BcalA); 
  
  EventTree->SetBranchAddress("NRecF",    &NRecF); 
  EventTree->SetBranchAddress("FcalE",    FcalE); 
  EventTree->SetBranchAddress("FcalA",    FcalA); 
  
  Long64_t nentries = EventTree->GetEntries();
  cout<<"Number of Entries: "<<nentries<<endl;
  
    for (Long64_t i=0; i<nentries;i++) {
  //for (Long64_t i=0; i<1800000;i++) {

    EventTree->GetEntry(i);
    
    hBCAL_Gen->Fill(GenA[0], GenE[0]);  
    hFCAL_Gen->Fill(GenA[0], GenE[0]);  
    

    //cout<< i<<" "<<NGen<<"   "<<GenE[0]<<"  "<<GenA[0]<<" : "<<NRecB<<"  "<<NRecF<<endl;

    Int_t binE = (Int_t)(GenE[0]/0.05);
    Int_t binA = (Int_t)((GenA[0]-5.)/0.08333333333333333);
    
    for (Int_t j=0;j<NRecB;j++){
      Double_t pa = fresBA[binA*4 + 1];
      Double_t pe = fresBE[binE*4 + 1];
      Float_t da = TMath::Abs( BcalA[j] - fresBA[binA*4 + 0] );
      Float_t de = TMath::Abs( BcalE[j] - fresBE[binE*4 + 0] ); 

      if ((da<3.*pa) && (de<3.*pe))
	hBCAL_Rec->Fill(GenA[0], GenE[0]);  
    }
    
    for (Int_t j=0;j<NRecF;j++){
      Double_t pa = fresFA[binA*4+1];
      Double_t pe = fresFE[binE*4+1];
      Float_t da = TMath::Abs( FcalA[j] - fresFA[binA*4 + 0] );
      Float_t de = TMath::Abs( FcalE[j] - fresFE[binE*4 + 0] );
      
      if ((da<3.*pa) && (de<3.*pe)){
	hFCAL_Rec->Fill(GenA[0], GenE[0]);  
      }
    }
    
  } // end event loop
  

  Int_t Ne = 64;
  Int_t Na = 240;
  
  TH1D *ratioAFcal[240];
  TH1D *ratioABcal[240];
  TH1D *h1 = hFCAL_Rec->ProjectionY("proj1",1,1,"eo"); 
  TH1D *h2 = hFCAL_Gen->ProjectionY("proj2",1,1,"eo");
  for (Int_t j=0;j<Na-1;j++){
    char str0[128];
    sprintf(str0,"ratioAFcal%d",j);
    char str1[128];
    Float_t a = j*20./240. + 20./480. + 5.;
    sprintf(str1,"Fcal Angle bin %d = %5.2f deg",j,a);
    h1->Reset();
    h2->Reset();
    h1 = hFCAL_Rec->ProjectionY("proj1",j+1,j+1,"eo"); 
    h2 = hFCAL_Gen->ProjectionY("proj2",j+1,j+1,"eo");

    ratioAFcal[j] = new TH1D(str0, str1, 64, 0., 3.2);
    ratioAFcal[j]->Sumw2();
    ratioAFcal[j]->GetXaxis()->SetTitle("Energy gen. [GeV]");
    ratioAFcal[j]->GetYaxis()->SetTitle("Efficiency");
    ratioAFcal[j]->Divide(h1, h2);
    
    sprintf(str0,"ratioABcal%d",j);
    a = j*20./240. + 20./480. + 5.;
    sprintf(str1,"Bcal Angle bin %d = %5.2f deg",j,a);
    ratioABcal[j] = new TH1D(str0, str1, 64, 0., 3.2);
    ratioABcal[j]->GetXaxis()->SetTitle("Energy gen. [GeV]");
    ratioABcal[j]->GetYaxis()->SetTitle("Efficiency");
    ratioABcal[j]->Sumw2();
    h1->Reset();
    h2->Reset();
    h1 = hBCAL_Rec->ProjectionY("proj1",j+1,j+1,"eo");
    h2 = hBCAL_Gen->ProjectionY("proj2",j+1,j+1,"eo");
    ratioABcal[j]->Divide(h1, h2);
   
  }
  
  TH1D *ratioEFcal[Ne];
  TH1D *ratioEBcal[Ne];
  TH1D *h3 = hFCAL_Rec->ProjectionX("proj3",1,1,"eo"); 
  TH1D *h4 = hFCAL_Gen->ProjectionX("proj4",1,1,"eo");
  for (Int_t j=0;j<Ne-1;j++){
    char str0[128];
    sprintf(str0,"ratioEFcal%d",j);
    char str1[128];
    Float_t a = j*3.2/64. + 3.2/128.;
    sprintf(str1,"Fcal Energy bin %d = %5.2f GeV",j,a);
    ratioEFcal[j] = new TH1D(str0, str1, 240, 5., 25.);
    ratioEFcal[j]->GetXaxis()->SetTitle("Angle gen. [deg]");
    ratioEFcal[j]->GetYaxis()->SetTitle("Efficiency");
    ratioEFcal[j]->Sumw2();

    h3->Reset();
    h4->Reset();
    h3 = hFCAL_Rec->ProjectionX("proj3",j+1,j+1,"eo"); 
    h4 = hFCAL_Gen->ProjectionX("proj4",j+1,j+1,"eo");
    ratioEFcal[j]->Divide(h3, h4);
    
    sprintf(str0,"ratioEBcal%d",j);
    a = j*3.2/64. + 3.2/128.;
    sprintf(str1,"Bcal Energy bin %d = %5.2f GeV",j,a);
    ratioEBcal[j] = new TH1D(str0, str1, 240, 5., 25.);
    ratioEBcal[j]->GetXaxis()->SetTitle("Angle gen. [deg]");
    ratioEBcal[j]->GetYaxis()->SetTitle("Efficiency");
    ratioEBcal[j]->Sumw2();
    h3->Reset();
    h4->Reset();
    h3 = hBCAL_Rec->ProjectionX("proj3",j+1,j+1,"eo"); 
    h4 = hBCAL_Gen->ProjectionX("proj4",j+1,j+1,"eo");    
    ratioEBcal[j]->Divide(h3, h4);
    
  }


  Double_t bina[Na];
  Double_t yafcal[Na];
  Double_t dyafcal[Na];
  Double_t yabcal[Na];
  Double_t dyabcal[Na];
  Double_t safcal[Na];
  Double_t dsafcal[Na];
  Double_t sabcal[Na];
  Double_t dsabcal[Na];
  for (Int_t k=0;k<Na;k++){
    bina[k] = k*20./240. + 20./480. + 5.; 
    yafcal[k] = fresFA[k*4];
    dyafcal[k] = fresFA[k*4+2];
    yabcal[k] = fresBA[k*4];
    dyabcal[k] = fresBA[k*4+2];
    safcal[k] = fresFA[k*4+1];
    dsafcal[k] = fresFA[k*4+3];
    sabcal[k] = fresBA[k*4+1];
    dsabcal[k] = fresBA[k*4+3];
  }
  TGraphErrors *Afcal = new TGraphErrors(Na, bina, yafcal, NULL, dyafcal);
  Afcal->SetTitle("FCAL Reconstructed angle vs. generated angle");
  Afcal->GetXaxis()->SetTitle("Generated angle [deg]");
  Afcal->GetXaxis()->SetTitle("Reconstructed angle Peak pos. [deg]");
  TGraphErrors *Abcal = new TGraphErrors(Na, bina, yabcal, NULL, dyabcal);
  Abcal->SetTitle("BCAL Reconstructed angle vs. generated angle");
  Abcal->GetXaxis()->SetTitle("Generated angle [deg]");
  Abcal->GetXaxis()->SetTitle("Reconstructed angle Peak pos. [deg]");

  TGraphErrors *Asfcal = new TGraphErrors(Na, bina, safcal, NULL, dsafcal);
  Asfcal->SetTitle("FCAL #sigma Reconstructed angle vs. generated angle");
  Asfcal->GetXaxis()->SetTitle("Generated angle [deg]");
  Asfcal->GetXaxis()->SetTitle("#sigma reconstructed. angle [deg]");
  TGraphErrors *Asbcal = new TGraphErrors(Na, bina, sabcal, NULL, dsabcal);
  Asbcal->SetTitle("BCAL #sigma Reconstructed angle vs. generated angle");
  Asbcal->GetXaxis()->SetTitle("Generated angle [deg]");
  Asbcal->GetXaxis()->SetTitle("#sigma reconstructed. angle [deg]");


  Double_t bine[Ne];
  Double_t yefcal[Ne];
  Double_t dyefcal[Ne];
  Double_t yebcal[Ne];
  Double_t dyebcal[Ne];
  Double_t sefcal[Ne];
  Double_t dsefcal[Ne];
  Double_t sebcal[Ne];
  Double_t dsebcal[Ne];
  for (Int_t k=0;k<Ne;k++){
    bine[k] = k*3.2/64. + 3.2/128.;
    yefcal[k] = fresFE[k*4];
    dyefcal[k] = fresFE[k*4+2];
    yebcal[k] = fresBE[k*4];
    dyebcal[k] = fresBE[k*4+2];
    sefcal[k] = fresFE[k*4+1];
    dsefcal[k] = fresFE[k*4+3];
    sebcal[k] = fresBE[k*4+1];
    dsebcal[k] = fresBE[k*4+3];
  }
  TGraphErrors *Efcal = new TGraphErrors(Ne, bine, yefcal, NULL, dyefcal);
  Efcal->SetTitle("FCAL Reconstructed energy vs. generated energy [GeV]");
  Efcal->GetXaxis()->SetTitle("Generated energy [GeV]");
  Efcal->GetXaxis()->SetTitle("Reconstructed energy Peak pos. [GeV]");
  TGraphErrors *Ebcal = new TGraphErrors(Ne, bine, yebcal, NULL, dyebcal);
  Ebcal->SetTitle("BCAL Reconstructed energy vs. generated energy [GeV]");
  Ebcal->GetXaxis()->SetTitle("Generated energy [GeV]");
  Ebcal->GetXaxis()->SetTitle("Reconstructed energy Peak pos. [GeV]");

  TGraphErrors *Esfcal = new TGraphErrors(Ne, bine, sefcal, NULL, dsefcal);
  Esfcal->SetTitle("FCAL #sigma Reconstructed energy vs. generated energy [GeV]");
  Esfcal->GetXaxis()->SetTitle("Generated energy [GeV]");
  Esfcal->GetXaxis()->SetTitle("#sigma Reconstructed energy [GeV]");
  TGraphErrors *Esbcal = new TGraphErrors(Ne, bine, sebcal, NULL, dsebcal);
  Esbcal->SetTitle("BCAL #sigma Reconstructed energy vs. generated energy [GeV]");
  Esbcal->GetXaxis()->SetTitle("Generated energy [GeV]");
  Esbcal->GetXaxis()->SetTitle("#sigma Reconstructed energy [GeV]");
  
  
  TFile *OF = new TFile(outf, "RECREATE");
  OF->cd();
  hFCAL_Gen->Write();
  hBCAL_Gen->Write();
  hFCAL_Rec->Write();
  hBCAL_Rec->Write();

  Afcal->Write("Afcal");
  Asfcal->Write("Asfcal");
  Abcal->Write("Abcal");
  Asbcal->Write("Asbcal");

  Efcal->Write("Efcal");
  Esfcal->Write("Esfcal");
  Ebcal->Write("Ebcal");
  Esbcal->Write("Esbcal");

  for (Int_t j=0;j<Ne-1;j++){
    if (ratioEBcal[j]->GetEntries())
      ratioEBcal[j]->Write();
    
    if (ratioEFcal[j]->GetEntries())
      ratioEFcal[j]->Write();
  }
  
  for (Int_t j=0;j<Na-1;j++){
    if (ratioABcal[j]->GetEntries())
      ratioABcal[j]->Write();
    
    if (ratioAFcal[j]->GetEntries())
      ratioAFcal[j]->Write();
  }
  



  OF->Close();


}