#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);
double hFunTest(int n);
TH1D * convo(TH1D* hIn);
//MAIN
int main(int argc, char **argv) {
  char  hId[120];

  /////CREATE A TEST HISTOGRAM/////
  TH1D *hTest = new TH1D("hTest","",101,-4.5,4.5);
  int nBinsTest = hTest->GetNbinsX();
  for (int iBin = 0; iBin < nBinsTest; iBin++) {
    double xVal = hTest->GetBinCenter(iBin+1);
    hTest->SetBinContent(iBin+1,pow(xVal,2));
  }
  /////Test 1
  TH2D *hTest2d = new TH2D("hTest2d","",nBinsTest,-4.5,4.5,nBinsTest,-4.5,4.5);
  int nBinsTestX1 = hTest2d->GetNbinsX();
  int nBinsTestX2 = hTest2d->GetNbinsY();
  for (int iBinX1 = 1; iBinX1 <= nBinsTestX1; iBinX1++) {
    for (int iBinX2 = 1; iBinX2 <= nBinsTestX2; iBinX2++) {
      hTest2d->SetBinContent(iBinX1,iBinX2,0);
      if (iBinX1 == iBinX2) {
	double hVal = hTest->GetBinContent(iBinX1);
	hTest2d->SetBinContent(iBinX1,iBinX2,hVal);
      }
    }
  }
  /////Test 2
  //  TH2D *hTest2d2 = new TH2D("hTest2d2","",nBinsTest,-4.5,4.5,nBinsTest,-4.5,4.5);
  TH2D *hTest2d2 = (TH2D *) hTest2d->Clone("hTest2d2");
  for (int iBinX2 = 1; iBinX2 <= nBinsTestX2; iBinX2++) {
    double hVal = hTest2d2->GetBinContent(iBinX2,iBinX2);
    for (int iBinX1 = 1; iBinX1 <= nBinsTestX1; iBinX1++) {
      double newVal = hVal*hFunTest(iBinX1-iBinX2);
      hTest2d2->SetBinContent(iBinX1,iBinX2,newVal);
    }
  }
  TH1D *hConvoTest = hTest2d2->ProjectionX("hConvoTest");

  TH1D *hConvo = convo(hTest); 
  hConvo->SetName("hConvo");


  TFile *outFile = new TFile("./myStage1OutFile.root","RECREATE");
  outFile->cd();
  hTest->Write();
  hTest2d->Write();
  hTest2d2->Write();
  hConvoTest->Write();
  hConvo->Write();
  outFile->Close();
  return 0;
}


double hFunTest(int n){
  double hValTest = 0;
  if (n == 0){
    hValTest =1;
  } else {
    if (abs(n)<5) {
      hValTest = 1.0/abs(n);
    }
  }
  return hValTest;
}

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

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;

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