#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 prototype
double hFun(int n);

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

  TFile *inSino=new TFile("./medDat_all.root"); //Get the file that has histo
  TH2D* hSino=(TH2D*)inSino->Get("hSino"); //Get the 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();

  TH2D *hSinoF = new TH2D("hSinoF","",nBinsS,sMin,sMax,nBinsTheta,thetaMin,thetaMax);
  //TH1D *hSinoTest = new TH1D("hSinoTest","",nBinsS,sMin,sMax);
  TH1D *hSinoTest = new TH1D("hSinoTest","",9,-4.5,4.5);
  ////////////

  hSinoTest->Fill(-4.0,hFun(-4));
  hSinoTest->Fill(-3.0,hFun(-3));
  hSinoTest->Fill(-2.0,hFun(-2));
  hSinoTest->Fill(-1.0,hFun(-1));
  hSinoTest->Fill(0.0,hFun(0));
  hSinoTest->Fill(1.0,hFun(1));
  hSinoTest->Fill(2.0,hFun(2));
  hSinoTest->Fill(3.0,hFun(3));
  hSinoTest->Fill(4.0,hFun(4));

  //int iBinTheta = 0;
  //int iBinS = 320;
  //for (int iBinN = 0; iBinN < nBinsS; iBinN++) {
  //double thetaDeg = hSino->GetYaxis()->GetBinCenter(iBinTheta);
    //double hVal = hFun(iBinN-iBinS);
    //double hVal = hFun(iBinN-iBinS);
  //double nVal = hSino->GetXaxis()->GetBinCenter(iBinS);
  //hSinoTest->Fill(nVal,hVal);

  //}
  for (int iBinTheta = 0; iBinTheta < nBinsTheta; iBinTheta++) {
    double thetaDeg = hSino->GetYaxis()->GetBinCenter(iBinTheta);
    for (int iBinS = 0; iBinS < nBinsS; iBinS++) {
      double sVal = hSino->GetXaxis()->GetBinCenter(iBinS);
      double intensity = hSino->GetBinContent(iBinS,iBinTheta);
      for (int iBinN = iBinS-nBinsS; iBinN <= iBinS+nBinsS; iBinN++) {
	double nVal = hSinoF->GetXaxis()->GetBinCenter(iBinN);
	double hVal = hFun(iBinN-iBinS);
	hSinoF->Fill(nVal,thetaDeg,hVal*intensity);
      }
    }
  }

  /////////////
  TFile *outFile = new TFile("./myStage1OutFile.root","RECREATE");
  outFile->cd();
  hSino->Write();
  hSinoF->Write();
  hSinoTest->Write();
  outFile->Close();
  return 0;
}


double hFun(int n){
  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;
}