#include <iostream>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TMath.h"
#include "TStyle.h"
#include "TString.h"
using namespace std;

void recon_tree()
{

   gROOT->SetBatch(kTRUE);
   gStyle->SetOptStat("nemruoi"); // "nemr"

// ROOT data file referencing variables
	int nph, ev, k;
	float y, z, x;
	float y_adj, z_adj;

// Recon ROOT file variables

	int element, tot_hit_count;
	int hc[19];
	float y_seed, z_seed;

   cout << "Here" << endl;

// Setting up histogram array

   TH1F *h[1601][20];
   for(int a = 0; a <= 1600; a++)
     for(int f = 0; f <= 19; f++) {
       h[a][f] = new TH1F(TString::Format("h%d_at_c%d", f, a),
			  TString::Format("hist num: %d, c%d", f, a),
			  100, 0, 400);
       h[a][f]->SetDirectory(0); // hide it from gROOT
     }

// Opening raw data file and referencing variables within it

   TFile *datafile = TFile::Open("back_ref_20mm.root");
   TTree *datatree; datafile->GetObject("dstree", datatree);

   int nentries = datatree->GetEntries();

   datatree->SetBranchAddress("nph", &nph);
	datatree->SetBranchAddress("ev", &nph);
	datatree->SetBranchAddress("y", &y);
	datatree->SetBranchAddress("z", &z);

	float ph_y[1000], ph_z[1000], ph_wl[1000]; // [nph]
	int ph_pid[1000]; // [nph]

	datatree->SetBranchAddress("ph_y", ph_y);
	datatree->SetBranchAddress("ph_z", ph_z);
	datatree->SetBranchAddress("ph_pid", ph_pid);

// Opening ROOT file for writing and setting variables for it

   TFile* recon_file = TFile::Open("recon_tree.root","RECREATE");

	TTree* rtree = new TTree("rtree","Data for PDF reconstruction");

	rtree->Branch("element", &element, "element/I");
	rtree->Branch("y_seed", &y_seed, "y_seed/F");
	rtree->Branch("z_seed", &z_seed, "z_seed/F");
	rtree->Branch("tot_hit_count", &tot_hit_count, "tot_hit_count/I");
	rtree->Branch("hc", hc, "hc[19]/I");

// Histogram for element selection

   TH2F *h_p = new TH2F("h_p", "Element array", 40, -100, 100, 40, -100, 100);

// Iterating over all entries in the data file

	for(int i = 0; i < nentries; i++)
   {

   	fill_n(hc, 19, 0);        // setting all elements of "hit count" array to zero      

   	tot_hit_count = 0;        // setting total hit counts to zero

   	datatree->GetEntry(i);    // accessing current entry of the ROOT data file
	if ((nph < 0) || (nph >= 1000)) { std::cout << "Error : nph = " << nph << std::endl; } // just a precaution

   	y_adj = y*10;             // changing/adjusting y and z values to mm
      z_adj = z*10;

      y_seed = y_adj;
      z_seed = z_adj;

      element = h_p->FindBin(y_adj,z_adj,0);    // determing the element for reconstruction based on position
      if ((element < 0) || (element > 1600)) continue; // just a precaution

   // Iterating over all photons absorbed in photocathode and adding them to the appropriate "hit count" for PMT

      for(int j = 0; j < nph; j++)
      {
        	if(ph_pid[j] <= 28 && ph_pid[j] >= 10)
        	{
            k = ph_pid[j] - 10;

            hc[k]++; // k = [0 ... 18]

         }

      }

   // Determining the total hits from the hit count array for the current event

      for(int l = 0; l <= 18; l++)
      {
      	tot_hit_count += hc[l];
      }

      h[element][0]->Fill(tot_hit_count);    // Filling the total hits histogram for the current element

   // Filling histograms with the current hit count for that PMT

      for(int b = 1; b <= 19; b++)
      {
         h[element][b]->Fill(hc[b - 1]);

      }

   	rtree->Fill();
   }

   rtree->Write();
   h_p->Write();

// Iterating over all canvases and drawing and storing all histograms 

   TCanvas *c = new TCanvas("c", "multipads", 900, 700);
   c->Divide(4, 5);
   for(int m = 0; m <= 1600; m++)
   {
      c->SetName(TString::Format("c%d", m));
      c->SetTitle(TString::Format("multipads %d", m));
      for(int n = 0; n <= 19; n++)
      {
         c->cd(n + 1);
         h[m][n]->Draw();
	 gPad->Modified(); gPad->Update(); // make sure it's really (re)drawn
      }
      c->cd(0);
      // c->Modified(); c->Update(); // make sure it's really (re)drawn
      c->Write();
   }
   delete c; // no longer needed

   delete recon_file; // automatically deletes "rtree" and "h_p", too

   cout << "Done generating recon data ROOT file" << endl;

   delete datafile; // automatically deletes "datatree", too

   // cout << "Begin cleanup" << endl;
   for(int m = 0; m <= 1600; m++) for(int n = 0; n <= 19; n++) delete h[m][n];
   cout << "Done cleanup" << endl;
}