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

//Function prototypes //ASDF
TH1D* convo(TH1D *hIn); //ASDF
double hFun(int n); //ASDF

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


  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

  TH1D *hTmp;//ASDF

  char hID[120];

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

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

  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);
	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();
  outFile->Close();
  return 0;
}

double hFun(int n){    //CONVO UPDATE: added this line of code down to return hVal;
  double hVal = 0;
  int nAbs = abs(n);
  if (nAbs%2 == 0) hVal = 0;
  if (nAbs%2 == 1){
    hVal = -1.0*pow(n*3.14159,-2);
  }
  if (nAbs==0) hVal = 1.0/4.0;
  return hVal;
}
    
//CONVO UPDATE: removed } here, was getting error

//ASDF
TH1D* convo(TH1D *hIn){
  TH1D *hOut = (TH1D *) hIn->Clone("hOut");
  hOut->Reset();
  int nBins = hIn->GetNbinsX();
  for (int iBinX2 = 1; iBinX2 <= nBins; iBinX2++) {
    double hVal = hIn->GetBinContent(iBinX2);
    for (int iBinX1 = 1; iBinX1 <= nBins; iBinX1++) {
      double hValX1 = hOut->GetBinContent(iBinX1);
      //double newVal = hVal*hFunTest(iBinX1-iBinX2) + hValX1; //ASDF
      double newVal = hVal*hFun(iBinX1-iBinX2) + hValX1; //ASDF

      hOut->SetBinContent(iBinX1,newVal);
    }
  }
  return hOut;
}