#include "TStyle.h"
#include "TList.h"
#include "TH3.h"
#include "TF1.h"
#include "TFile.h"

Double_t transfer_function(Double_t *x, Double_t * /*param*/)
{
   //This function is required to suppress
   //boxes for empty bins - make them transparent.
   if (x)
      return *x > 0 ? 1. : 0.;
   
   return 0.;
}


void Plot3DHistos() {

  TFile* f1 = new TFile("na22.25.root");
  TH3D* histo_or = (TH3D*) f1->Get("MakeVoxels_reconstructedDEMO/EnergyHits");

  // Put shorter ranges for x, y and z, that contain completely the track, 
  // but not much larger, for better visualization
  double xmin = -100.;
  double xmax = 100.;
  double ymin = -100.;
  double ymax = 100.;
  double zmin = 140.;
  double zmax = 280.;
  
  
  double binsize = 10.;
  double nbinsx = (xmax-xmin)/binsize; // here you need to put the size of the voxel, in mm
  double nbinsy = (ymax-ymin)/binsize;
  double nbinsz = (zmax-zmin)/binsize;

  TH3D* histo = new TH3D("H", "H", nbinsx, xmin, xmax, nbinsy, ymin, ymax, nbinsz, zmin, zmax);

  for (int binX=1;binX<histo_or->GetNbinsX(); binX++) {
     for (int binY=1;binY<histo_or->GetNbinsY(); binY++) {
        for (int binZ=1;binZ<histo_or->GetNbinsZ(); binZ++) {
	  if (histo_or->GetBinContent(binX, binY, binZ) != 0) {
	    //std::cout << binX << ", " <<  binY << ", " <<  binZ << ", " << std::endl;
	    histo->Fill(histo_or->GetXaxis()->GetBinCenter(binX), 
			histo_or->GetYaxis()->GetBinCenter(binY),
			histo_or->GetZaxis()->GetBinCenter(binZ), 
			histo_or->GetBinContent(binX, binY, binZ));	      
	  }
	}
     }
  }
  // histo->Draw();
  
  TList * const lof = histo->GetListOfFunctions();
   if (!lof) {
      std::cout<<"List of functions is null\n";
      delete histo;
      return;
   }

   lof->Add(new TF1("TransferFunction", transfer_function, 0., 1000., 0));

  gStyle->SetCanvasPreferGL(kTRUE);
  //gStyle->SetPalette(55);
  
   histo->Draw("glcolz");
  
}