// @(#)root/html:$Name:  $:$Id: THtml.cxx,v 1.124 2006/11/24 16:01:02 brun Exp $
// Author: Mathieu Benoit <mailto:mbenoit@umontreal.ca> 

/*************************************************************************
 * . Class CTDR_Field					                 *
 *   Written by : Mathieu Benoit                                         *
 *   m.benoit@umontreal.ca						 *
 *************************************************************************/

#include "CTDR_Field.h"
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <string>
using namespace std;

//////////////////////////////////////////////////////////////////////////		        
// This class represent the CTDR_Field in the CZT detector. CTDR_Field 	//
// data are read from a ROOT  file and stored in a TTree.This class pass// 
// to the simulation the value of the CTDR_Field at a position x,y,z,	//
// on demand. This class is similar to the class Detecteur.		//		        
//////////////////////////////////////////////////////////////////////////

ClassImp(CTDR_Field)


//______________________________________________________________________________
CTDR_Field::CTDR_Field(TFile &Fichier, int NdivX, int NdivY, int NdivZ,double LX,double LY, double LZ)
{
// Constructor for the CTDR_Field 
// We retrieve 3 TTree from a file. These data tree contain data of the CTDR_Field inside a logical pixel of
// our detector. The data are obtained from a finite-emelent analysis software, COMSOL multiphysics, and
// we retreive an array of 100*100*100 for nodes inside the model at constant distance. Each data point
// is assigned a unique ID calculated from its position. It allow us to get the memory adress of a CTDR_Field value
// by calculating the ID associated to the position where we need the value. We then retreive the CTDR_Field
// value from the tree using this unique ID pointing to its adress in the tree's memory 
// The data trees are created one time , prior to simulation, using Treereader.c ROOT script .
// This allow to modify the script to use CTDR_Field Data from various sources without touching this code. 
// The Storage of the data in .root files allow to get in memory large amount of data, compressed with
// an efficient method that allow a fast acces to data in the tree. S 
// parameters: (reference to the ROOT file, number of division in X (100), in Y, in Z, size of the model
// in X,Y,Z)

                  char hash[100];// A character chain to store the hash function
                   NDivx=NdivX;
                   NDivy=NdivY;
                   NDivz=NdivZ;
                   Lx=LX;
                   Ly=LY;
                   Lz=LZ;
		   //We retreive the TTree from the File
                   TEx=(TTree*)Fichier.Get("Ex");
                   TEy=(TTree*)Fichier.Get("Ey");
                   TEz=(TTree*)Fichier.Get("Ez");
		   //We build the ID->Memory adress Index
                   TEx->BuildIndex("I","0");
                   TEx->SetBranchAddress("x",&x_Ex);
                   TEx->SetBranchAddress("y",&y_Ex);
                   TEx->SetBranchAddress("z",&z_Ex);
                   TEx->SetBranchAddress("E",&E_Ex);
                   TEy->BuildIndex("I","0");
                   TEy->SetBranchAddress("x",&x_Ey);
                   TEy->SetBranchAddress("y",&y_Ey);
                   TEy->SetBranchAddress("z",&z_Ey);
                   TEy->SetBranchAddress("E",&E_Ey);
                   TEz->BuildIndex("I","0");
                   TEz->SetBranchAddress("x",&x_Ez);
                   TEz->SetBranchAddress("y",&y_Ez);
                   TEz->SetBranchAddress("z",&z_Ez);
                   TEz->SetBranchAddress("E",&E_Ez);
		   //We print some info on the tree to be sure it is well loaded
                   TEx->Print();
                   TEy->Print();
                   TEz->Print();
};

//______________________________________________________________________________
double CTDR_Field::GetEx()
{
//Data access function for the CTDR_Field components
     return Ex;
};

//______________________________________________________________________________
double CTDR_Field::GetEy()
{
//Data access function for the CTDR_Field components
       return Ey;
};     

//______________________________________________________________________________
double CTDR_Field::GetEz() 
{
//Data access function for the CTDR_Field components
       return Ez;
};


//______________________________________________________________________________
void CTDR_Field::SetChamp(double xx, double yy, double zz)
{
//This function Put in the public variable Ex,Ey,Ez, the CTDR_Field value for the point xx,yy,zz
    double x, y, z;
     //The model is assumed to be periodic in X,Y
     if(fabs(xx)>Lx/2|fabs(yy)>Ly/2){
                                   x=fmod(xx+Lx/2,Lx)-Lx/2;
                                   y=fmod(yy+Ly/2,Ly)-Ly/2;
                                   z=zz;
                                   }
     //No CTDR_Field Data over or under the box
     else if(zz>=Lz){
          x=xx;
          y=yy;
          z=Lz;
          }
     else if(z<0){
          x=xx;
          y=yy;
          z=0;
          }
     else {
          x=xx;
          y=yy;
          z=zz;
          };
     //We interpolate the CTDR_Field at xx,yy,zz, from the surrounding points in the data tree
     Interpole(x,y,z);
     };


//______________________________________________________________________________  
void CTDR_Field::Interpole(double x, double y, double z)
{ 
// This function calculate the linear 3D interpolation of the Field at x,y,z, 
//   1. We calculate the ID of the 8 surrounding points to (x,y,z)
//   2. We get the Field value at he adress corresponding to these ID. 
//   3  We compute the interpolation and put the interpolated values in Ex,Ey,Ez  
   	;
     //Variable for the CTDR_Field at the 8 points
     double E000,E010,E100,E110,E001,E011,E101,E111;
     //variable for variable change x,y,z->t,u,v
     double t,u,v;

     /*Ex computation*/
     SetPoints(x,y,z);//We calculate the Id of the 8 surrounding points
     //We retreive the value of the CTDR_Field at the 8 points using their ID
     Long64_t Adresse=Hash(lowx,lowy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E000=E_Ex;
     //cout << "V000 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E010=E_Ex;
     //cout <<"V010 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E100=E_Ex;
     //cout << "V100 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E110=E_Ex;
     //cout << "V110 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,lowy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E001=E_Ex;
     //cout << "V001 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E011=E_Ex;
     //cout <<"V011 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E101=E_Ex;
     //cout << "V101 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E111=E_Ex;
     //cout << "V111 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;

     //We do a variable change to locate ourself in a cube around origin
     if(highx!=lowx){t=(x-Lx*(lowx/NDivx-0.5))/(Lx*(highx/NDivx-0.5)-Lx*(lowx/NDivx-0.5));}
     else{t=0;};
     if(highy!=lowy){u=(y-Ly*(lowy/NDivy-0.5))/(Ly*(highy/NDivy-0.5)-Ly*(lowy/NDivy-0.5));}
     else{u=0;};
     if(highz!=lowz){v=(z-Lz*(lowz/NDivz))/(Lz*(highz/NDivz)-Lz*(lowz/NDivz));}
     else{v=0;};

     //We compute the interpolation
     Ex=E000*(1 - t)*(1 - u)*(1 - v) +E100*t*(1 - u)*(1 - v) +
                E010*(1 - t)* u *(1 - v) + E001*(1 - t)*(1 - u)*v +
                E101*t*(1 - u)*v + E011*(1 - t)*u*v + E110*t*u*(1 - v) + E111*t*u*v;

     /*Ey computation*/

     Adresse=Hash(lowx,lowy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E000=E_Ey;
     //cout << "V000 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E010=E_Ey;
     //cout <<"V010 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E100=E_Ey;
     //cout << "V100 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E110=E_Ey;
     //cout << "V110 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,lowy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E001=E_Ey;
     //cout << "V001 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E011=E_Ey;
     //cout <<"V011 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E101=E_Ey;
     //cout << "V101 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E111=E_Ey;
     //cout << "V111 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     if(highx!=lowx){t=(x-Lx*(lowx/NDivx-0.5))/(Lx*(highx/NDivx-0.5)-Lx*(lowx/NDivx-0.5));}// on calcule une fois seulement
     else{t=0;};
     if(highy!=lowy){u=(y-Ly*(lowy/NDivy-0.5))/(Ly*(highy/NDivy-0.5)-Ly*(lowy/NDivy-0.5));}// on calcule une fois seulementy
     else{u=0;};
     if(highz!=lowz){v=(z-Lz*(lowz/NDivz))/(Lz*(highz/NDivz)-Lz*(lowz/NDivz));}// on calcule une fois seulement
     else{v=0;};
     Ey=E000*(1 - t)*(1 - u)*(1 - v) +E100*t*(1 - u)*(1 - v) +
                E010*(1 - t)* u *(1 - v) + E001*(1 - t)*(1 - u)*v +
                E101*t*(1 - u)*v + E011*(1 - t)*u*v + E110*t*u*(1 - v) + E111*t*u*v;

     /*Ez computation */

     Adresse=Hash(lowx,lowy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E000=E_Ez;
     //cout << "V000 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E010=E_Ez;
     //cout <<"V010 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E100=E_Ez;
     //cout << "V100 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E110=E_Ez;
     //cout << "V110 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,lowy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E001=E_Ez;
     //cout << "V001 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E011=E_Ez;
     //cout <<"V011 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E101=E_Ez;
     //cout << "V101 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E111=E_Ez;
     //cout << "V111 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     if(highx!=lowx){t=(x-Lx*(lowx/NDivx-0.5))/(Lx*(highx/NDivx-0.5)-Lx*(lowx/NDivx-0.5));}// on calcule une fois seulement
     else{t=0;};
     if(highy!=lowy){u=(y-Ly*(lowy/NDivy-0.5))/(Ly*(highy/NDivy-0.5)-Ly*(lowy/NDivy-0.5));}// on calcule une fois seulementy
     else{u=0;};
     if(highz!=lowz){v=(z-Lz*(lowz/NDivz))/(Lz*(highz/NDivz)-Lz*(lowz/NDivz));}// on calcule une fois seulement
     else{v=0;};
     Ez=E000*(1 - t)*(1 - u)*(1 - v) +E100*t*(1 - u)*(1 - v) +
                E010*(1 - t)* u *(1 - v) + E001*(1 - t)*(1 - u)*v +
                E101*t*(1 - u)*v + E011*(1 - t)*u*v + E110*t*u*(1 - v) + E111*t*u*v;
  ;
};




//______________________________________________________________________________
double CTDR_Field::GetLz(){
//Access method for the Model depth in z
       return Lz;
       };