#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <string>
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TFile.h"
using namespace std;


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)));
   } 
   return rms;
}


void mc() {
   evenement Event[1300];

   ifstream montecarlo("select_output.txt", ios::in);

   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;
         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;
         }
      }
   }

   double rmsxmoy[13], rmsymoy[13], rmszmoy[13];
   TCanvas *canx,   *cany,   *canz;
   TCanvas *can1,   *can2,   *can3;
   TH1F    *histox, *histoy, *histoz;

   for (int i=0; i<13; i++) {
      rmsxmoy[i] = 0;
      rmsymoy[i] = 0;
      rmszmoy[i] = 0;
      histox = new TH1F(Form("histox%d",i), Form("histox%d",i), 20, 0, 2.0);
      histoy = new TH1F(Form("histoy%d",i), Form("histoy%d",i), 20, 0, 2.0);
      histoz = new TH1F(Form("histoz%d",i), Form("histoz%d",i), 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] = histox->GetMean();
      rmsymoy[i] = histoy->GetMean();
      rmszmoy[i] = histoz->GetMean();  
      canx = new TCanvas(Form("rmsx = %ikeV", Energies[i]));  
      histox->Draw("");
      cany = new TCanvas(Form("rmsy = %ikeV", Energies[i]));  
      histoy->Draw("");
      canz = new TCanvas(Form("rmsz = %ikeV", Energies[i]));  
      histoz->Draw("");
   }

   TGraph *gr  = new TGraph(13, En, rmsxmoy);
   TGraph *gr2 = new TGraph(13, En, rmsymoy);	
   TGraph *gr3 = new TGraph(13, En, rmszmoy);

   can1 = 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");

   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");

   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");
	
   montecarlo.close();
   return;
}