#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;

int main(int argc, char **argv) {


  //TFile *inFile = new TFile("./AB_V1.root"); //Get the file that has histogram
  TFile *inFile = new TFile("./medDat_all.root"); //Get the new file that has histogram
  TH2D *hThetavsS=(TH2D*) inFile->Get("hSino"); //Get the histogram 

  int nBinsTheta = hThetavsS->GetNbinsY();//Get number of angles on Y axis 
  TH1D *hThetaSlice[nBinsTheta];//Array of 1D histogram created from hThetavsS

  char hID[120];

  for (int iBinTheta=0; iBinTheta<nBinsTheta; iBinTheta++){
    sprintf(hID,"hThetaSlice_%d",iBinTheta);
    
   
    hThetaSlice[iBinTheta] = hThetavsS->ProjectionX(hID,iBinTheta+1,iBinTheta+1,"");//Make 1D histogram slice at iBinTheta+1

    //cout<<"thetaValDegree = "<<thetaValDegree<<endl;
  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  
  //TH1D *hTheta45 = hThetavsS->ProjectionX("hTheta45",45,45,"");//Make 1D histogram at 45, want to remove at some point

  double sMax = hThetavsS->GetXaxis()->GetXmax();//Relies on hThetavsS
  double sMin = hThetavsS->GetXaxis()->GetXmin();//Relies on hThetavsS
  int nBinsS = hThetavsS->GetNbinsX();//Relies on hThetavsS

  int nBinsX = nBinsS;
  double xMin = sMin;
  double xMax = sMax;

  int nBinsY = nBinsS;
  double yMin = sMin;
  double yMax = sMax;

  //double thetaDegree = 45;//Adjust theta angle//Hardcoded degree
  //double thetaVal = thetaDegree*3.14159/180;//Hardcoded angle

  TH2D *hXY = new TH2D("hXY","",nBinsX,xMin,xMax,nBinsY,yMin,yMax);

  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  int nBinsXY = hXY->GetNbinsY();//Gets number for 1D histograms
  TH2D *hXYPaint[nBinsXY];//Array of 1D histogram created from hXY that paints the intensity at all angles
  TH2D *hXYPaintIntegrated = new TH2D("hXYPaintIntegrated","Integration Results",nBinsX,xMin,xMax,nBinsY,yMin,yMax);

  char hxyID[120];

  for (int iBinTheta=0; iBinTheta<nBinsTheta; iBinTheta++){
    sprintf(hxyID,"hXYPaint_%d", iBinTheta);
    
    hXYPaint[iBinTheta] = new TH2D(hxyID,"",nBinsX,xMin,xMax,nBinsY,yMin,yMax);
    
    double thetaValDegree = hThetavsS->GetYaxis()->GetBinCenter(iBinTheta);//Obtain angle for a give iBinTheta
    double thetaVal = thetaValDegree*3.14159/180;
    
    for (int iBinX = 1; iBinX <= nBinsX; iBinX++) {
      for (int iBinY = 1; iBinY <= nBinsY; iBinY++) {
	double xVal = hXY->GetXaxis()->GetBinCenter(iBinX);
	double yVal = hXY->GetYaxis()->GetBinCenter(iBinY);
	double sVal = xVal*cos(thetaVal) + yVal*sin(thetaVal);//For all theta values
	
	int sBin = hThetavsS->GetXaxis()->FindBin(sVal);
	double intensityVal = hThetaSlice[iBinTheta]->GetBinContent(sBin);
	
	//if(intensityVal != 0) {//Only prints when intensityVal does not equal zero.
	
	//cout<<"ix,iy: x,y,s = "<<iBinX<<" , "<<iBinY<<" : "<<xVal<<" , "<<yVal<<" ; "<<sVal<<endl;//Prints ix,ix and x,y,s values.
	//cout<<"intensity = "<<intensityVal<<endl;//Prints intensity value for all sVal.
	
	//double intensityValInitial = hXY->GetBinContent(iBinX,iBinY);//Get initial intensity value at each loop of iBinX and iBinY

	double intensityValInitial = hXYPaintIntegrated->GetBinContent(iBinX,iBinY);

	hXYPaint[iBinTheta]->SetBinContent(iBinX,iBinY,intensityVal+intensityValInitial);
	hXYPaintIntegrated->SetBinContent(iBinX,iBinY,intensityVal+intensityValInitial);

	}
      } 
    }

  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  TFile *outFile = new TFile("./myOutFile.root","RECREATE");
  outFile->cd();
  for (int iBinTheta=0; iBinTheta<nBinsTheta; iBinTheta++){
    
    hThetaSlice[iBinTheta]->Write();
  }
  for (int iBinTheta=0; iBinTheta<nBinsTheta; iBinTheta++){
   
    hXYPaint[iBinTheta]->Write();
  }
  hXYPaintIntegrated->Write();
  hXY->Write();
  //hTheta45->Write();
  //hThetaSlice[45]->Write();
  outFile->Close();
  return 0;


  }