//Inclusion des librairies utiles au programme
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <string>
#include "/home/gcm/mbenoit/root/include/TH1F.h"
#include "/home/gcm/mbenoit/root/include/TH2F.h"
#include "/home/gcm/mbenoit/root/include/TH3F.h"
#include "/home/gcm/mbenoit/root/include/TGraph.h"
#include "/home/gcm/mbenoit/root/include/TCanvas.h"
#include "/home/gcm/mbenoit/root/include/TFile.h"
using namespace std;

//gROOT->Reset();

typedef struct {double x; double y; double z; double energied;} interaction ;
typedef struct { int Nevent; int Ninteraction; interaction site[210];double energie;} evenement;
int Energies[13]={10, 20, 30, 50, 70, 100, 150, 200, 300, 400, 500, 700, 1000};
double En[13]={10, 20, 30, 50, 70, 100, 150, 200, 300, 400, 500, 700, 1000};
int Nphotons=100;


double rmsevent(int i, int j,int xyz,evenement Event[]){
  double rms=0, carre=0,moy=0,Et=0; 
 if(xyz==0){for(int k=0;k<Event[10*i+j].Ninteraction;k++){
    moy+=Event[10*i+j].site[k].x*(Event[10*i+j].site[k].energied);
    carre+=pow(Event[10*i+j].site[k].x,2)*(Event[10*i+j].site[k].energied);
    Et+=Event[10*i+j].site[k].energied;
  }; 
 rms=sqrt(fabs(carre/Et-pow(moy/Et,2)));
 };
if(xyz==1){for(int k=0;k<Event[10*i+j].Ninteraction;k++){
    moy+=Event[10*i+j].site[k].y*(Event[10*i+j].site[k].energied);
    carre+=pow(Event[10*i+j].site[k].y,2)*(Event[10*i+j].site[k].energied);
    Et+=Event[10*i+j].site[k].energied;
  }; 
rms=sqrt(fabs(carre/Et-pow(moy/Et,2)));
 };
if(xyz==2){for(int k=0;k<Event[10*i+j].Ninteraction;k++){
    moy+=Event[10*i+j].site[k].z*(Event[10*i+j].site[k].energied);
    carre+=pow(Event[10*i+j].site[k].z,2)*(Event[10*i+j].site[k].energied);
    Et+=Event[10*i+j].site[k].energied;
  };
  rms=sqrt(fabs(carre/Et-pow(moy/Et,2)));
}; 
   
  
/*   
     if(xyz==0){for(int k=0;k<Event[10*i+j].Ninteraction;k++){
    moy+=Event[10*i+j].site[k].x/Event[10*i+j].Ninteraction;
    carre+=pow(Event[10*i+j].site[k].x,2)/Event[10*i+j].Ninteraction;
  }; 
 rms=sqrt(fabs(carre-pow(moy,2)));
 };
if(xyz==1){for(int k=0;k<Event[10*i+j].Ninteraction;k++){
    moy+=Event[10*i+j].site[k].y/Event[10*i+j].Ninteraction;
    carre+=pow(Event[10*i+j].site[k].y,2)/Event[10*i+j].Ninteraction;
  }; 
rms=sqrt(fabs(carre-pow(moy,2)));
 };
if(xyz==2){for(int k=0;k<Event[10*i+j].Ninteraction;k++){
    moy+=Event[10*i+j].site[k].z/Event[10*i+j].Ninteraction;
    carre+=pow(Event[10*i+j].site[k].z,2)/Event[10*i+j].Ninteraction;
  }; 
rms=sqrt(fabs(carre-pow(moy,2)));
 };   */ 
 return rms;
};
 
    
  

int mc(){

char nom[100];

ifstream montecarlo("/home/gcm/mbenoit/Monte-carlo_CZT/Montecarlo", ios::in);




//TFile fi("Histo.root","new");

evenement Event[1300];

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

	for(int j=0; j<Nphotons;j++){
		      Event[10*i+j].energie=double(Energies[i]);
		      montecarlo >> Event[10*i+j].Nevent;
		      montecarlo >> Event[10*i+j].Ninteraction;
		      //cout << Event[10*i+j].Nevent << ' ' << Event[10*i+j].Ninteraction << endl ; 

			for(int k=0;k<Event[10*i+j].Ninteraction;k++){
				
				montecarlo >>Event[10*i+j].site[k].x;
				montecarlo >>Event[10*i+j].site[k].y;
				montecarlo >>Event[10*i+j].site[k].z;
				montecarlo >>Event[10*i+j].site[k].energied;
				//cout << Event[10*i+j].site[k].x << ' ' << Event[10*i+j].site[k].y << ' ' <<  Event[10*i+j].site[k].z <<endl;

			};
	
		};
	
	};
double rmsxmoy[13],rmsymoy[13],rmszmoy[13];
	
cout << "Lecture du fichier termine. appuyer sur une touche " << endl;	
 
 	for(int i=0;i<13;i++){
	 rmsxmoy[i]=0;
	 rmsymoy[i]=0;
	 rmszmoy[i]=0;

	 TH1F *histox = new TH1F("histo","h",20,0,2.0);
	 TH1F *histoy = new TH1F("histo","h",20,0,2.0);
	 TH1F *histoz = new TH1F("histoz","h",20,0,2.0);
	 histox->SetBit(TH1::kCanRebin);
	 histoy->SetBit(TH1::kCanRebin);
	 histoz->SetBit(TH1::kCanRebin);
	  for(int j=0; j<Nphotons;j++){
	   histox->Fill(rmsevent(i,j,0,Event));
	   histoy->Fill(rmsevent(i,j,1,Event));
	   histoz->Fill(rmsevent(i,j,2,Event));	   
 	  rmsxmoy[i]+=rmsevent(i,j,0,Event)/Nphotons;
	  rmsymoy[i]+=rmsevent(i,j,1,Event)/Nphotons;
	  rmszmoy[i]+=rmsevent(i,j,2,Event)/Nphotons; 
	  //cout <<  Event[10*i+j].Nevent << endl ;
	  //cout << sqrt(pow( rmsevent(i,j,0,Event) ,2 ) +  pow( rmsevent(i,j,1,Event) ,2 )) << ' ' << rmsevent(i,j,2,Event) << endl ;
	  };
/*  	  rmsxmoy[i]=histox->GetMean();
	  rmsymoy[i]=histoy->GetMean();
	  rmszmoy[i]=histoz->GetMean();  */
   	
	  //cout << rmsxymoy[i]<< endl;
	  //cout << rmszmoy[i]<< endl;
	  
	};
TGraph *gr= new TGraph(13,En,rmsxmoy);
TGraph *gr2= new TGraph(13,En,rmsymoy);	
TGraph *gr3= new TGraph(13,En,rmszmoy);

TCanvas *can = new TCanvas("rms_x_moy");
gr->SetTitle("rayon rms x initial moyen pour l'interaction d'un gamma dans le CZT");
gr->GetXaxis()->SetTitle("Energie(keV)");
gr->GetYaxis()->SetTitle("rayon rms (mm)");
gr->Draw("AC");

TCanvas *can2 = new TCanvas("rms_y_moy");
gr2->SetTitle("rayon rms y initial moyen pour l'interaction d'un gamma dans le CZT");
gr2->GetXaxis()->SetTitle("Energie(keV)");
gr2->GetYaxis()->SetTitle("rayon rms (mm)");
gr2->Draw("AC");

TCanvas *can3 = new TCanvas("rms_z_moy");
gr3->SetTitle("rayon rms z initial moyen pour l'interaction d'un gamma dans le CZT");
gr3->GetXaxis()->SetTitle("Energie(keV)");
gr3->GetYaxis()->SetTitle("rayon rms (mm)");
gr3->Draw("AC");
	
/* for(int i=0;i<13;i++){
	double rms[13];
	sprintf(nom,"E=%ikeV",Energies[i]);
	TCanvas *can = new TCanvas(nom);
	sprintf(nom,"Spectre dans le CZT a E=%i keV",Energies[i]);
	TH1F *spectre = new TH1F("spectre",nom,100,0.99995*double(Energies[i]),1.00005*double(Energies[i]));	
	for(int j=0; j<Nphotons;j++){
		double E=0;
		for(int k=0;k<Event[10*i+j].Ninteraction;k++){
			E+=Event[10*i+j].site[k].energied;
		};
		spectre->Fill(E);
	};
	rms[i]=1e3*spectre->GetRMS(1);
	can->Draw();
	spectre->Draw("");
};
double En[13]=={10, 20, 30, 50, 70, 100, 150, 200, 300, 400, 500, 700, 1000};

TCanvas *can = new TCanvas("rms_raies");
TGraph *graph = new TGraph(13,En,rms);
graph->GetXaxis()->SetTitle("Energie(keV)");
graph->GetYaxis()->SetTitle("Energie rms (eV)");
graph->SetTitle("energie RMS pour un spectre delta de gamma dans le CZT");

graph->Draw("A*"); */			
		

montecarlo.close();
return 0;
};