#include "TFile.h"
#include "TH1D.h"

#include "PlotPdf.hh"

//Arguments:
//PDF file
//root outfile
//E
//D
//cd
//pmta

PlotPdf::PlotPdf(){
  pdf = DigitalShowerPdf();
  pdf.init();
}

PlotPdf::PlotPdf(const string infile){
  pdf.infile = infile;
  pdf.init();

  //TFile f( infile.c_str() );
  
  //getobj( f, "h3xa" ,    H3 );
  //getobj( f, "h3normx" , G3 );
  
  //H3.Divide(&G3);
  
  //f.Close();
}

void PlotPdf::set_E   (const double E)   { _E    = E;    };
void PlotPdf::set_D   (const double D)   { _D    = D;    };
void PlotPdf::set_cd  (const double cd)  { _cd   = cd;   };
void PlotPdf::set_pmta(const double pmta){ _pmta = pmta; };
void PlotPdf::set_all(const double E,
		      const double D,
		      const double cd,
		      const double pmta){
  set_E(E);
  set_D(D);
  set_cd(cd);
  set_pmta(pmta);
}

void PlotPdf::plot_cd(const string outfile, const double E, const double D, const double pmta){
  
  TFile out( outfile.c_str(), "RECREATE");
  
  TH1D* H1 = new TH1D("h1", "", 200, -1, 1);
  
  const double scaled_E = E*1e-6; //H3 is scaled to 1 PeV
  for(int icd = 0; icd < H1->GetNbinsX(); icd ++){
    const double cd = H1->GetBinCenter(icd);
    H1->SetBinContent(icd, scaled_E*pdf.get_h3_value(pdf.H3, D, cd, pmta));
  }
  out.Write();
}

void PlotPdf::plot_cd(const string outfile){
  plot_cd(outfile, _E, _D, _pmta);
}

void PlotPdf::plot_D(const string outfile, const double E, const double cd, const double pmta){
  
  TFile out( outfile.c_str(), "RECREATE");
  
  TH1D* H1 = new TH1D("h1", "", 200, 0, 1000);
  
  const double scaled_E = E*1e-6; //H3 is scaled to 1 PeV
  for(int iD = 0; iD < H1->GetNbinsX(); iD ++){
    const double D = H1->GetBinCenter(iD);
    H1->SetBinContent(iD, scaled_E*pdf.get_h3_value(pdf.H3, D, cd, pmta));
  }
  out.Write();
}

void PlotPdf::plot_D(const string outfile){
  plot_D(outfile, _E, _cd, _pmta);
}

void PlotPdf::plot_pmta(const string outfile, const double E, const double D, const double cd){
  
  TFile out( outfile.c_str(), "RECREATE");
  
  TH1D* H1 = new TH1D("h1", "", 200, -1, 1);
  
  const double scaled_E = E*1e-6; //H3 is scaled to 1 PeV
  for(int ipmta = 0; ipmta < H1->GetNbinsX(); ipmta ++){
    const double pmta = H1->GetBinCenter(ipmta);
    H1->SetBinContent(ipmta, scaled_E*pdf.get_h3_value(pdf.H3, D, cd, pmta));
  }
  out.Write();
}

void PlotPdf::plot_pmta(const string outfile){
  plot_pmta(outfile, _E, _D, _cd);
}