#include "TFile.h"
#include "TTree.h"
#include "TH2.h"
#include "TF1.h"

void offsetcor(string rootfilename, string outfilename)
{

  ofstream outfile(outfilename.c_str(), ofstream::out);

  TFile *rootFile = new TFile(rootfilename.c_str(),"OPEN");
  TTree* tree;
  rootFile->GetObject("t",tree);
  tree->Print();

  TFile *newrootFile = new TFile("offsetcor.root","RECREATE");
  
  TH2D* th2_beamSignal_b0_ic1[30];
  TH2D* th2_beamSignal_b1_ic1[30];
  TH2D* th2_beamSignal_b2_ic1[30];
  TH2D* th2_beamSignal_b3_ic1[30];

  TProfile * prof_beamSignal_b0_ic1[30];
  TProfile * prof_beamSignal_b1_ic1[30];
  TProfile * prof_beamSignal_b2_ic1[30];
  TProfile * prof_beamSignal_b3_ic1[30];

  char histname[30] = "";
  char fitname[30] = "";

  TCanvas*  c1;
  TF1 * fit_beamSignal_b0_ic1[30];
  TF1 * fit_beamSignal_b1_ic1[30];
  TF1 * fit_beamSignal_b2_ic1[30];
  TF1 * fit_beamSignal_b3_ic1[30];

			       
  for (int iii=0; iii<=29; iii++){
    ///b0
    sprintf(histname,"th2_beamSignal_b0_ic1_%d",iii);
    sprintf(fitname,"fit_beamSignal_b0_ic1_%d",iii);
    th2_beamSignal_b0_ic1[iii] = new TH2D(histname,histname,100,0,500,100,0,50000);
    tree->Project(histname,"beamSignal_b0:ic1_avg",Form("energy_index==%i",1+iii*10),"coloz");
    prof_beamSignal_b0_ic1[iii] = th2_beamSignal_b0_ic1[iii]->ProfileX();
    fit_beamSignal_b0_ic1[iii] = new TF1(fitname,"pol1",10,150);
    prof_beamSignal_b0_ic1[iii]->Fit(fit_beamSignal_b0_ic1[iii],"R");

    ///b1
    sprintf(histname,"th2_beamSignal_b1_ic1_%d",iii);
    sprintf(fitname,"fit_beamSignal_b1_ic1_%d",iii);
    th2_beamSignal_b1_ic1[iii] = new TH2D(histname,histname,100,0,500,100,0,50000);
    tree->Project(histname,"beamSignal_b1:ic1_avg",Form("energy_index==%i",1+iii*10),"coloz");
    prof_beamSignal_b1_ic1[iii] = th2_beamSignal_b1_ic1[iii]->ProfileX();
    fit_beamSignal_b1_ic1[iii] = new TF1(fitname,"pol1",10,150);
    prof_beamSignal_b1_ic1[iii]->Fit(fit_beamSignal_b1_ic1[iii],"R");

        ///b2
    sprintf(histname,"th2_beamSignal_b2_ic1_%d",iii);
    sprintf(fitname,"fit_beamSignal_b2_ic1_%d",iii);
    th2_beamSignal_b2_ic1[iii] = new TH2D(histname,histname,100,0,500,100,0,50000);
    tree->Project(histname,"beamSignal_b2:ic1_avg",Form("energy_index==%i",1+iii*10),"coloz");
    prof_beamSignal_b2_ic1[iii] = th2_beamSignal_b2_ic1[iii]->ProfileX();
    fit_beamSignal_b2_ic1[iii] = new TF1(fitname,"pol1",10,150);
    prof_beamSignal_b2_ic1[iii]->Fit(fit_beamSignal_b2_ic1[iii],"R");

    ///b3
    sprintf(histname,"th2_beamSignal_b3_ic1_%d",iii);
    sprintf(fitname,"fit_beamSignal_b3_ic1_%d",iii);
    th2_beamSignal_b3_ic1[iii] = new TH2D(histname,histname,100,0,500,100,0,50000);
    tree->Project(histname,"beamSignal_b3:ic1_avg",Form("energy_index==%i",1+iii*10),"coloz");
    prof_beamSignal_b3_ic1[iii] = th2_beamSignal_b3_ic1[iii]->ProfileX();
    fit_beamSignal_b3_ic1[iii] = new TF1(fitname,"pol1",10,150);
    prof_beamSignal_b3_ic1[iii]->Fit(fit_beamSignal_b3_ic1[iii],"R");

    
    outfile << fit_beamSignal_b0_ic1[iii]->GetParameter(0) << " "  << fit_beamSignal_b1_ic1[iii]->GetParameter(0)  << " "  << fit_beamSignal_b2_ic1[iii]->GetParameter(0)  << " "  << fit_beamSignal_b3_ic1[iii]->GetParameter(0) << std::endl;
 }
 
  newrootFile->Write();
  newrootFile->Close();
  outfile.close();
  rootFile->Close();
}