// Sara Vanini: simple root macro to read and analyse muon radiography .img file
// April, 1st 2008 
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <iostream>
#include "TFile.h"
#include "TMath.h"
#include "TNtuple.h"
#include "TH3F.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TStyle.h"


//customize default plain style e.g.

gROOT->SetStyle("Plain");
gStyle->SetOptStat(111111);
gStyle->SetOptFit(1);
gStyle->SetCanvasPreferGL(true);
gStyle->SetPalette(1);
gStyle->SetPaletteOpacity(1.);

using namespace std;

// define tree structure
 struct vox_struct{
   Int_t n;
   Double_t x;
   Double_t y;
   Double_t z;
   Int_t weight;
 };

// define objects for plots
TH3F * hvoxel; 
TTree * tree;

void imageAnalysis(char *rundir){

  //instantiate objects for plots
  hvoxel = new TH3F("hvoxel","Voxel density",200,0.5,200.5,200,0.5,200.5,200,0.5,200.5); 

  // compose file names
  char hdr_file_name[80];
  sprintf(hdr_file_name,"/data/radmu/Image3D/%s/xx016.hdr",rundir);
  char img_file_name[80];
  sprintf(img_file_name,"/data/radmu/Image3D/%s/xx016.img",rundir);
  cout << "\nAnalysing Header file : " << hdr_file_name << endl;
  cout << "Analysing Image file : " << img_file_name << endl;

  // open .hdr file to read number of voxels: nx, ny, nz
  short int nbuffer[3];  
  FILE * hdrFile = fopen ( hdr_file_name , "rb" );
  if (hdrFile==NULL) 
  {
	cout << "Error opening file... please check!" << endl;
	exit (1);
  }
  fseek( hdrFile , 42 , SEEK_SET );
  fread( &nbuffer,2,3,hdrFile);
  cout << "\nNumber of Voxels retrieved from header file:" << endl;
  cout << "nx = " << nbuffer[0] << " ny = " << nbuffer[1] << " nz = " << nbuffer[2] << endl;
  int nx = nbuffer[0];
  int ny = nbuffer[1];
  int nz = nbuffer[2];
   
  // number of cm/voxel
  int cm_per_vox = 2;

  // open .img file
  FILE * imgFile;
  imgFile = fopen ( img_file_name , "rb" );

  if (imgFile==NULL) 
  {
	cout << "Error opening file... please check!" << endl;
	exit (1);
  }

  Double_t wmin = 1;
  Double_t wmax = 1;

  // read file and store voxels in array (not a matrix for memory-savings!)
  // voxel[1_y1_z1,...nx_y1_z1,1_y2_z1,..nx_y2_z2,..nx_ny_nz ]
  Float_t voxel[2550000];// dimension: nx_max * ny_max * nz_max
  int out = fread (voxel,sizeof(float),nx*ny*nz,imgFile);
  if (out != nx*ny*nz) 
  {
	cout << "File reading error, number of elements read " << out << endl; 
	exit (3);
  }

  // if(!feof(imgFile))
	//cout << "Error: end of file not reached... " << endl;	

  // create a tree for storing data
  vox_struct vox;
  tree = new TTree("T","ROOT tree with density voxels");
  tree->Branch("vox",&vox.n,"n/I:x/D:y:z:weight/I");  

  // retrive data from array, fill tree for later analysis
  int count = 0; 		//voxel counter
  float myDensityCut = 1; 	//cut on density
  for(int zc=0; zc<nz; zc++)
 	for(int yc=0; yc<ny; yc++)
		for(int xc=0; xc<nx; xc++)
		{	
			//cout 	<< "voxel " << xc << ", " << yc << ", " << zc << ": "
			//cout	<< voxel[xc + nx*yc + nx*ny*zc] << " ";//endl; 
			count++;
  			vox.n = count;		
			vox.x = xc;
			vox.y = yc;
			vox.z = zc;
			if(voxel[xc + nx*yc + nx*ny*zc]*100 < myDensityCut)
				vox.weight = 1;
			else
				vox.weight = int(voxel[xc + nx*yc + nx*ny*zc]*100);

			tree->Fill();

			if(wmin > vox.weight)
				wmin = vox.weight;
			if(wmax < vox.weight)
				wmax = vox.weight;		
		}

  // close files
  fclose (hdrFile);
  fclose (imgFile);

  // print min and max
  cout << "Max density " << wmax << endl;
  cout << "Min density " << wmin << endl;

  // loop on tree to put your favourite cuts, to fill histograms, etc...
  // get number of voxels
  Int_t voxnum = tree->GetEntries();
  cout << "\nTotal number of voxels from tree " << voxnum << endl;
  for (Int_t ientry=0 ; ientry < voxnum; ientry ++)
  {
        //read entry from tree
        tree->GetEntry(ientry);
	//dump voxels, just for test
	//if(vox.weight != 1)
        //	cout 	<< "vox number " << vox.n << ", x " << vox.x << ", y " << vox.y 
	//		<< ", z " << vox.z << ", density " << vox.weight << endl;

	// fill your favourite histogram
	//if(vox.x>20 && vox.x<120 && vox.y>20 && vox.y<100 && vox.z>20 && vox.z<70 )	
	hvoxel->Fill(vox.x,vox.y,vox.z,vox.weight);
  }

  // this for color palette
   Double_t Red[3]    = { 0.00, 1.00, 1.00};
   Double_t Green[3]  = { 0.00, 1.00, 0.00};
   Double_t Blue[3]   = { 1.00, 1.00, 0.00};
   Double_t w    = 0.7;
   Double_t v    = (w-wmin)/(wmax-wmin);
   Double_t Length[3] = { 0.00, v, 1.00 };   
   Int_t nb=50;
   TColor::CreateGradientColorTable(3,Length,Red,Green,Blue,nb);
 /*
  TCanvas * c;
  c = new TCanvas("voxel","voxels",800,800);

  // plot result in canvas
  cout << "\n...now plotting..." << endl;
  c->Divide(1,2);
  c->cd(1);  
  //plot tree
  tree->Draw("z:y:x:weight","");
  // plot histo
  c->cd(2); 
  hvoxel->Draw("openGL"); 

  // another Canvas to play with...
  TCanvas * cc = new TCanvas();
  cc->cd();
*/

  hvoxel->Draw("glboxz");
  return;

}