#include <string.h>
#include <iostream>
#include <vector>
#include "TFile.h"
#include "TLorentzVector.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TRandom3.h"
#include "TVector3.h"
#include <fstream>

using namespace std;
 
//MAIN
int main(int argc, char **argv) {
  char  hId[120];

  //TFile *inSino=new TFile("./medDat_all.root"); //Get the file that has histo
  TFile *inSino=new TFile("./myStage1OutFile.root"); //Get the file that has histos
  TH2D* hSino=(TH2D*)inSino->Get("hSino"); //Get a histogram 
  TH2D* hSinoF=(TH2D*)inSino->Get("hSinoF"); //Get a histogram 
  double sMax = hSino->GetXaxis()->GetXmax();
  double sMin = hSino->GetXaxis()->GetXmin();
  //double thetaMax = hSino->GetYaxis()->GetXmax();
  //double thetaMin = hSino->GetYaxis()->GetXmin();
  ////////////

  int nBinsS     = hSino->GetNbinsX();
  int nBinsTheta = hSino->GetNbinsY();

  TH1D *hSinoSlice[nBinsTheta];
  TH2D *hBP[nBinsTheta];


  for (int iBinTheta = 0; iBinTheta < nBinsTheta; iBinTheta++) {
    sprintf(hId,"hSinoSlice_%d",iBinTheta);
    hSinoSlice[iBinTheta] = (TH1D*) hSino->ProjectionX(hId,iBinTheta+1,iBinTheta+1,"");

    hSinoSlice[iBinTheta]->Smooth(5);

    sprintf(hId,"hBP_%d",iBinTheta);
    hBP[iBinTheta] = new TH2D(hId,"",nBinsS,sMin,sMax,nBinsS,sMin,sMax);
    
  }


  TH2D *hBPall = new TH2D("hBPall","",nBinsS,sMin,sMax,nBinsS,sMin,sMax);
  TH2D *hBPFall = new TH2D("hBPFall","",nBinsS,sMin,sMax,nBinsS,sMin,sMax);

  nBinsS = 4*nBinsS;
  int nBinsX = nBinsS;  int nBinsY = nBinsS;


  for (int iBinTheta = 0; iBinTheta < nBinsTheta; iBinTheta++) {
    for (int iBinY = 0; iBinY < nBinsY; iBinY++) {
      for (int iBinX = 0; iBinX < nBinsX; iBinX++) {
	//double xVal = hBP[iBinTheta]->GetXaxis()->GetBinCenter(iBinX);
	//double yVal = hBP[iBinTheta]->GetYaxis()->GetBinCenter(iBinY);
	double xVal = hBPFall->GetXaxis()->GetBinCenter(iBinX);
	double yVal = hBPFall->GetYaxis()->GetBinCenter(iBinY);
	double thetaDeg = hSino->GetYaxis()->GetBinCenter(iBinTheta);
	double thetaRad = thetaDeg*3.14159/180.0;
	double sVal = xVal*cos(thetaRad) + yVal*sin(thetaRad);
	int sBin = hSino->GetXaxis()->FindBin(sVal);

	double intensity = hSino->GetBinContent(sBin,iBinTheta);
	//hBP[iBinTheta]->SetBinContent(iBinX,iBinY,intensity);
	intensity += hBPall->GetBinContent(iBinX,iBinY);
	hBPall->SetBinContent(iBinX,iBinY,intensity);

	double intensityF = hSinoF->GetBinContent(sBin,iBinTheta);


	//hBP[iBinTheta]->SetBinContent(iBinX,iBinY,intensity);
	intensityF += hBPFall->GetBinContent(iBinX,iBinY);
	hBPFall->SetBinContent(iBinX,iBinY,intensityF);




      }
    }
  }



  //hBPall->RebinX(4);
  //hBPall->RebinY(4);

  //hBPFall->RebinX(4);
  //hBPFall->RebinY(4);

  hBPall->Scale(-1.0);
  hBPFall->Scale(-1.0);

  //hBPall->GetXaxis()->SetRangeUser(-sqrt(2*sMax)*0.95,sqrt(2*sMax)*0.95);
  //hBPall->GetYaxis()->SetRangeUser(-sqrt(2*sMax)*0.95,sqrt(2*sMax)*0.95);
  //hBPFall->GetXaxis()->SetRangeUser(-sqrt(2*sMax)*0.95,sqrt(2*sMax)*0.95);
  //hBPFall->GetYaxis()->SetRangeUser(-sqrt(2*sMax)*0.95,sqrt(2*sMax)*0.95);

  /////////////
  TFile *outFile = new TFile("./myOutFile.root","RECREATE");
  outFile->cd();
  hSino->Write();
  hBPall->Write();
  hBPFall->GetXaxis()->SetRangeUser(-20,20);
  hBPFall->GetYaxis()->SetRangeUser(-20,20);
  hBPFall->Scale(-1);
  hBPFall->SetMaximum(-250);
  hBPFall->SetMinimum(-1000);
  hBPFall->Write();
  for (int iBinTheta = 0; iBinTheta < nBinsTheta; iBinTheta++) {
    //hSinoSlice[iBinTheta]->Write();
    //hBP[iBinTheta]->Write();
  }
  
  outFile->Close();
  return 0;
}

// gStyle->SetPalette(52,0)

////////////
//Rebin by 2 in x and y
//hBPFall->GetXaxis()->SetRangeUser(-20,20)
//hBPFall->GetYaxis()->SetRangeUser(-20,20)
//hBPFall->Scale(-1)
//hBPFall->SetMaximum(-250)
//hBPFall->SetMinimum(-1000)

//gStyle->SetPalette(8,0)
//hBPFall->Draw("colz")