// @(#)root/html:$Name:  $:$Id: THtml.cxx,v 1.124 2006/11/24 16:01:02 $
// Author: Mathieu Benoit <mailto:mbenoit@umontreal.ca>
/*************************************************************************
 * . Class CTDR_Detecteur					         *
 *   Written by : Mathieu Benoit                                         *
 *   m.benoit@umontreal.ca						 *
 *************************************************************************/
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string>
#include "TString.h"
#include "CTDR_Detecteur.h"
using namespace std;

const double Pi=M_PI;
//////////////////////////////////////////////////////////////////////////				
// This class represent the weigting potential  in the CZT detector. 	//
// WP Data are read from a  ROOT  file and stored in a TTree. This class//
// pass to the simulation the value of the WP at a position x,y,z, 	//
// on demand. This class is similar to the class Field,			//
// but with some added function.					//				
//////////////////////////////////////////////////////////////////////////

ClassImp(CTDR_Detecteur)


//______________________________________________________________________________
CTDR_Detecteur::CTDR_Detecteur(TFile &Fichier, int NdivX, int NdivY, int NdivZ,double LX,double LY, double LZ,double periode)
{
// Constructor for the weighting potential 
// We retrieve 2 TTree from a file. These data tree contain data of the weigthing potential inside a logical pixel of
// our detector, for the X, and Y channel. 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 weighting potential value
// by calculating the ID associated to the position where we need the value. We then retreive the weighting potential
// 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  his allow to modify the script to use weighting potential 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.  
// 
// See class Field for comments on function uncommented here, they are clones that will eventually be
// merged as one single polymorphic class.                  
                   char hash[100];
                   NDivx=NdivX;
                   NDivy=NdivY;
                   NDivz=NdivZ;
                   Lx=LX;
                   Ly=LY;
                   Lz=LZ;
		   this->periode=periode;
                   TVx=(TTree*)Fichier.Get("Vx");
                   TVy=(TTree*)Fichier.Get("Vy");
                   //sprintf(hash,"int((%i*%i)*(x+%f/2)/(%f/%i))+int((%i)*(y+%f/2)/(%f/%i))",NDivx,NDivx,Lx,Lx,NDivx,NDivy,Ly,Ly,NDivy);
                   //cout << hash << endl;
                   TVx->BuildIndex("I","0");
                   TVx->SetBranchAddress("x",&x_Vx);
                   TVx->SetBranchAddress("y",&y_Vx);
                   TVx->SetBranchAddress("z",&z_Vx);
                   TVx->SetBranchAddress("V",&V_Vx);
                   TVy->BuildIndex("I","0");
                   TVy->SetBranchAddress("x",&x_Vy);
                   TVy->SetBranchAddress("y",&y_Vy);
                   TVy->SetBranchAddress("z",&z_Vy);
                   TVy->SetBranchAddress("V",&V_Vy);
                   //TVy->Scan("x:y:z:V","z==0.0075");


};	



//______________________________________________________________________________
void CTDR_Detecteur::Potentiel(double x , double y ,double z)
{ 
// For the weigting potential , we pay special attention to the values at the surface,
// So we switch top two dimensional interpolation close to the surface, mainly to avoid some bug
// seen when using 3-D interpolation.          


           if(z>(Lz-1e-6)){
                          Interpole2D(x,y);
                          //cout << " 2D" << endl;
                          }
           else if(z<0 | fabs(x)>Lx | fabs(y)>Ly){
                          Vx=0;
                          Vy=0;
                          //cout << " 0" << endl;
                          }
           else{
                Interpole3D(x,y,z);
                //cout << " 3D" << endl;
                };

};


//______________________________________________________________________________
void CTDR_Detecteur::Interpole3D(double x, double y , double z)
{
// See CTDR_Field 
     double V000,V010,V100,V110,V001,V011,V101,V111;
     double t,u,v;
     SetPoints(x,y,z);
     /*Calcul de Vx*/
     Long64_t Adresse=Hash(lowx,lowy,lowz);
     TVx->GetEntryWithIndex(Adresse,0);
     V000=V_Vx;
     //cout << "V000 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,lowz);
     TVx->GetEntryWithIndex(Adresse,0);
     V010=V_Vx;
     //cout <<"V010 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,lowz);
     TVx->GetEntryWithIndex(Adresse,0);
     V100=V_Vx;
     //cout << "V100 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,lowz);
     TVx->GetEntryWithIndex(Adresse,0);
     V110=V_Vx;
     //cout << "V110 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,lowy,highz);
     TVx->GetEntryWithIndex(Adresse,0);
     V001=V_Vx;
     //cout << "V001 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,highz);
     TVx->GetEntryWithIndex(Adresse,0);
     V011=V_Vx;
     //cout <<"V011 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,highz);
     TVx->GetEntryWithIndex(Adresse,0);
     V101=V_Vx;
     //cout << "V101 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,highz);
     TVx->GetEntryWithIndex(Adresse,0);
     V111=V_Vx;
     //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;};
     Vx=V000*(1 - t)*(1 - u)*(1 - v) +V100*t*(1 - u)*(1 - v) +
                V010*(1 - t)* u *(1 - v) + V001*(1 - t)*(1 - u)*v +
                V101*t*(1 - u)*v + V011*(1 - t)*u*v + V110*t*u*(1 - v) + V111*t*u*v;

     /*Calcul de Ey*/

     Adresse=Hash(lowx,lowy,lowz);
     TVy->GetEntryWithIndex(Adresse,0);
     V000=V_Vy;
      //cout <<"V000 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;
     Adresse=Hash(lowx,-highy,lowz);
     TVy->GetEntryWithIndex(Adresse,0);
     V010=V_Vy;
     //cout <<"V010 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;
     Adresse=Hash(highx,lowy,lowz);
     TVy->GetEntryWithIndex(Adresse,0);
     V100=V_Vy;
     // cout <<"V100 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;
     Adresse=Hash(highx,-highy,lowz);
     TVy->GetEntryWithIndex(Adresse,0);
     V110=V_Vy;
     // cout <<"V110 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;
     Adresse=Hash(lowx,lowy,highz);
     TVy->GetEntryWithIndex(Adresse,0);
     V001=V_Vy;
     // cout <<"V001 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;
     Adresse=Hash(lowx,highy,highz);
     TVy->GetEntryWithIndex(Adresse,0);
     V011=V_Vy;
     // cout <<"V011 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;
     Adresse=Hash(highx,lowy,highz);
     TVy->GetEntryWithIndex(Adresse,0);
     V101=V_Vy;
     // cout <<"V101 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;
     Adresse=Hash(highx,highy,highz);
     TVy->GetEntryWithIndex(Adresse,0);
     V111=V_Vy;
     //cout <<"V111 " <<  x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;

     Vy=V000*(1 - t)*(1 - u)*(1 - v) +V100*t*(1 - u)*(1 - v) +
                V010*(1 - t)* u *(1 - v) + V001*(1 - t)*(1 - u)*v +
                V101*t*(1 - u)*v + V011*(1 - t)*u*v + V110*t*u*(1 - v) + V111*t*u*v;
 };

//______________________________________________________________________________ 
void CTDR_Detecteur::Interpole2D(double x, double y)
{
// The 2-D interpolator, works like the 3-D one.
     double t,u;
     double V00,V01,V10,V11;

     SetPoints(x,y,Lz);
     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;};

     Long64_t  Adresse=0;
     Adresse= Hash(lowx,lowy,NDivz);
     TVx->GetEntryWithIndex(Adresse,0);
     V00=V_Vx;

     //cout <<"V00 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;

     Adresse=Hash(lowx,highy,NDivz);
     TVx->GetEntryWithIndex(Adresse,0);
     V01=V_Vx;

     //cout << "V01 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;

     Adresse=Hash(highx,lowy,NDivz);
     TVx->GetEntryWithIndex(Adresse,0);
     V10=V_Vx;

     //cout << "V10 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;

     Adresse=Hash(highx,highy,NDivz);
     TVx->GetEntryWithIndex(Adresse,0);
     V11=V_Vx;

     //cout << "V11 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
    // cout << "(t,u): " << t << ' ' <<  u << endl ; 

     Vx=V00*(1-t)*(1-u) + V10*t*(1-u) + V01*(1-t)*u +V11*t*u;



     Adresse=Hash(lowx,lowy,NDivz);
     TVy->GetEntryWithIndex(Adresse,0);
     V00=V_Vy;

     //cout << "V00 " << x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;

     Adresse=Hash(lowx,highy,NDivz);
     TVy->GetEntryWithIndex(Adresse,0);
     V01=V_Vy;

     //cout << "V01 "  << x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;

     Adresse=Hash(highx,lowy,NDivz);
     TVy->GetEntryWithIndex(Adresse,0);
     V10=V_Vy;

     //cout << "V10 " << x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;

     Adresse=Hash(highx,highy,NDivz);
     TVy->GetEntryWithIndex(Adresse,0);
     V11=V_Vy;

     //cout << "V11 "  << x_Vy <<' ' << y_Vy << ' ' << z_Vy << ' ' << V_Vy << endl;

     Vy=V00*(1-t)*(1-u) + V10*t*(1-u) + V01*(1-t)*u +V11*t*u;
};

//______________________________________________________________________________
void CTDR_Detecteur::Detecte(CTDR_Interaction &inter)
{
// This methods calculate the integration constant for inducted charge integral
// for a whole CTDR_Interaction :
//
// Inducted charge=sum of final weigting potential - sum of initial weighting potentials 

    Vxf=0;
    Vyf=0;
    for(int i=0;i<inter.GetN();i++){
            if(fabs(inter.GetCharge(i).GetX())>Lx/2  | fabs(inter.GetCharge(i).GetY())>Ly/2 )
            {Vxf+=0;
             Vyf+=0;}
            else{ Potentiel(inter.GetCharge(i).GetX(),inter.GetCharge(i).GetY(),inter.GetCharge(i).GetZ());
                  Vxf+=Vx*inter.GetCharge(i).GetQ();
                  Vyf+=Vy*inter.GetCharge(i).GetQ();
            };
            };

   inter.SetQx(Vxf);
   inter.SetQy(Vyf);
};

//______________________________________________________________________________
void CTDR_Detecteur::Detecte(CTDR_Interaction &inter,ofstream &off)
{
// This function compute the inducted charge during simulated and write to a file 
// the inducted charge in X, and Y channel each time the function is evoked

   Vxf=0;
    Vyf=0;
    for(int i=0;i<inter.GetN();i++){
            if(fabs(inter.GetCharge(i).GetX())>Lx/2  | fabs(inter.GetCharge(i).GetY())>Ly/2  )
            {Vxf+=0;
             Vyf+=0;}
            else if(inter.GetCharge(i).GetZ()>=Lz){Potentiel(inter.GetCharge(i).GetX(),inter.GetCharge(i).GetY(),Lz);
                  Vxf+=Vx*inter.GetCharge(i).GetQ();
                  Vyf+=Vy*inter.GetCharge(i).GetQ();
                  }
            else{ Potentiel(inter.GetCharge(i).GetX(),inter.GetCharge(i).GetY(),inter.GetCharge(i).GetZ());
                  Vxf+=Vx*inter.GetCharge(i).GetQ();
                  Vyf+=Vy*inter.GetCharge(i).GetQ();
            };
            };
    Vxf-=inter.GetQx();
    Vyf-=inter.GetQy();
    off << inter.Gett() << " " << Vxf << " " << Vyf << endl;
    cout  <<"(t,Vx,Vy): " << inter.Gett() << " " << Vxf << " " << Vyf << endl;

};

//______________________________________________________________________________

void CTDR_Detecteur::Detecte(CTDR_Interaction &inter,ofstream &off,int N_canaux)
{
// This function compute the inducted charge during simulated and write to a file 
// the inducted charge in X, and Y channel each time the function is evoked

    double *Vxt=new double[2*N_canaux+1];
    double *Vyt=new double[2*N_canaux+1];
    for(int canal=-N_canaux;canal<=N_canaux;canal++)
    {
	    Vyt[canal+N_canaux]=0;
	    Vxt[canal+N_canaux]=0;
    };

    for(int canalx=-N_canaux;canalx<=N_canaux;canalx++)
    {
	    for(int canaly=-N_canaux;canaly<=N_canaux;canaly++)
	    {
		for(int i=0;i<inter.GetN();i++){
			if(fabs(inter.GetCharge(i).GetX()-canalx*periode)>Lx/2  | fabs(inter.GetCharge(i).GetY()-canaly*periode)>Ly/2  )
			{
				Vxt[canalx+N_canaux]+=0;
				Vyt[canaly+N_canaux]+=0;
			}
			else if(inter.GetCharge(i).GetZ()>=Lz){Potentiel(inter.GetCharge(i).GetX()-canalx*periode,inter.GetCharge(i).GetY()-canaly*periode,Lz);
				Vxt[canalx+N_canaux]+=Vx*inter.GetCharge(i).GetQ();
				Vyt[canaly+N_canaux]+=Vy*inter.GetCharge(i).GetQ();
				}
			else{ Potentiel(inter.GetCharge(i).GetX()-canalx*periode,inter.GetCharge(i).GetY()-canaly*periode,inter.GetCharge(i).GetZ());
				Vxt[canalx+N_canaux]+=Vx*inter.GetCharge(i).GetQ();
				Vyt[canaly+N_canaux]+=Vy*inter.GetCharge(i).GetQ();
			};
		};
	    };};

	off << inter.Gett();
 	cout << inter.Gett();
	TString str , str2;

	for(int canal=0;canal<=2*N_canaux;canal++)
	    {
		    str.Form(" %6.4f ",Vxt[canal]);
		    str2.Form(" %6.4f ",Vyt[canal]);
		    off << str <<  str2;
		    cout << str <<  str2;

	    };
	off << endl;
	cout << endl;
	delete Vxt;
	delete Vyt;
};


//______________________________________________________________________________ 
void CTDR_Detecteur::Detectef(CTDR_Interaction &inter,ofstream &off)
{
// This function is used to write the final inducted charge to a file.
// The difference with the previous function is that the time is not written
// I use this function to accumulate in a file final charge values of many events to create scatter plots

    Vxf=0;
    Vyf=0;
    for(int i=0;i<inter.GetN();i++){
            if(fabs(inter.GetCharge(i).GetX())>Lx/2  | fabs(inter.GetCharge(i).GetY())>Ly/2  )
            {  off << inter.GetCharge(i).GetX() << " " << inter.GetCharge(i).GetY() << " " << inter.GetCharge(i).GetZ() << " " << 0 << " " << 0 << ' ' << inter.GetCharge(i).GetQ()<< endl;
               cout << "hors pixel " <<endl;

            }
            else{ if(inter.GetCharge(i).GetZ()>=Lz){Potentiel(inter.GetCharge(i).GetX(),inter.GetCharge(i).GetY(),Lz);
                  off << inter.GetCharge(i).GetX() << " " << inter.GetCharge(i).GetY() << " " << inter.GetCharge(i).GetZ() << " " << Vx << " " << Vy << ' ' << inter.GetCharge(i).GetQ()<< endl;
                  Vxf+=Vx*inter.GetCharge(i).GetQ();
                  Vyf+=Vy*inter.GetCharge(i).GetQ();
                  cout << "collecte " <<endl;
                  }
                  else{Potentiel(inter.GetCharge(i).GetX(),inter.GetCharge(i).GetY(),inter.GetCharge(i).GetZ());
                  off << inter.GetCharge(i).GetX() << " " << inter.GetCharge(i).GetY() << " " << inter.GetCharge(i).GetZ() << " " << Vx << " " << Vy << ' ' << inter.GetCharge(i).GetQ()<< endl;
                  Vxf+=Vx*inter.GetCharge(i).GetQ();
                  Vyf+=Vy*inter.GetCharge(i).GetQ();
                  cout << "piege " <<endl; };};
                  };

    Vxf-=inter.GetQx();
    Vyf-=inter.GetQy();



};

//______________________________________________________________________________
double CTDR_Detecteur::GetVx()
{
// Data acces method for the wighting potential value for X channel
       return Vxf;
};

//______________________________________________________________________________
double CTDR_Detecteur::GetVy()
{
// Data acces method for the wighting potential value for Y channel
       return Vyf;
};

//______________________________________________________________________________
void CTDR_Detecteur::Test(double x , double y ,double z )
{
// A test function, outputs the potential for a given position
     Potentiel(x,y,z);
     cout << "(z,Vx,Vy): " << z << ' ' << Vx << ' ' << Vy << endl;
};