#include <TMath.h>
#include "TFile.h"
#include "TH1.h"
#include "TString.h"
#include "TGraph.h"
#include "TRandom3.h"
#include "TCanvas.h"
#include "TDirectory.h"
#include <TSystem.h>
using namespace std;


void Visualise_Fluctuation3(){
const Int_t nbfiles = 10;//100; //240 //number of files
const Int_t nLayers = 16; //number of layers
const Int_t nMicrons= 300; //number of slices
const Int_t nParticles=100000; // nb of particles

const Int_t Emaxdep= 4000;
  Int_t xrange= 4600;
  Int_t yrange=4600;

Double_t edep;
Double_t eventID;
Double_t Detector;
Double_t step;

TH1::AddDirectory(kFALSE);

vector<string> ion= {"alpha", "ion_6_12_6","ion_7_14_7","ion_14_28_14","ion_26_56_26"};
//{"proton","alpha", "He3", "ion_3_7_3", "ion_4_9_4", "ion_5_11_5", "ion_6_12_6", "ion_7_14_7", "ion_8_16_8", "ion_9_19_9", "ion_10_20_10", "ion_11_23_11", "ion_12_24_12" , "ion_13_27_13", "ion_14_28_14", "ion_26_56_26"};
Int_t nspecies = ion.size(); // number of particles to plot

vector<Int_t> ener= {1,3,4,5,6,7,8,9,10,11,12,13,14,15,20,22,24,26,28,30,50,70,100};
//{1, 2, 3, 5, 8, 10, 12, 15, 18, 20, 25, 30, 40, 50, 60, 80, 100}
Int_t nener = ener.size(); // number of particles to plot
vector<Int_t> angle= {0};
//{0, 1, 2, 5, 8, 10, 15, 20, 30, 40, 50, 70};

Int_t nangle = angle.size(); // number of particles to plot

//Int_t nI=4;Int_t nE=3;Int_t nA=0; //7 12
TString title = TString("Emax Distribution ")+ Form("%d", nParticles);
TH2D *hEmax = new TH2D("Emax distrib",title,nLayers*nMicrons/2,0,nLayers*nMicrons,Emaxdep/30,0,Emaxdep);


cout<<nener<<endl;
double_t  En;
Int_t range;
Int_t realRange;
Int_t c =0;

for (Int_t nI = 1; nI < 2; ++nI)
{
  for (Int_t nE = 2; nE < 3; ++nE) //14  //21
  {

  TString current_read_file_name = TString("processed/config300_")+ ion[nI] + TString("_")+ Form("%d", ener[nE])+ TString("MeV_nucl_") + Form("%d", angle[0]) +TString("_deg_emp.root");
  cout<< "File "<< current_read_file_name << " opened" << endl;
  TFile *myfile= new TFile(current_read_file_name,"READ");
  gDirectory->pwd();
  myfile->cd("Particles");

  double_t EMParticle[nParticles]={0};
  Int_t RangeParticle[nParticles]={0};
  TIter next(gDirectory->GetListOfKeys());
  TKey *key;
  c=0;
  while ((key = (TKey*)next())&& (c<nParticles)) {

    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1F *myh1 = (TH1F*)key->ReadObj();
    TString name=myh1->GetName();
    TString current_particle=myh1->GetName();
    Int_t first= current_particle.First('_');
    current_particle.Remove(0,first+1);
    first= current_particle.First('_');
    current_particle.Remove(first,current_particle.Length());
    Int_t a = current_particle.Atoi();
    Int_t last= name.Last('_');
    name.Remove(0,last+1);
    Int_t b=name.Atoi(); //layer value
    for (Int_t j = 0; j < myh1->GetNbinsX()+1; ++j)
    {
      range=myh1->GetBin(j) * myh1->GetBinWidth(j);
      realRange=range + 300*(b-1);
      En=myh1->GetBinContent(myh1->GetBin(j));
      if (En>0)
        {
          if (En>EMParticle[a]){
            EMParticle[a]=En;
            RangeParticle[a]=realRange;
          }
        }

    }
    c=c+1;
    delete myh1;
  }

  for (Int_t k = 0; k < nParticles; ++k)
  {
    if(k<nParticles){
      hEmax->Fill(RangeParticle[k],EMParticle[k]);
    }
  }


  myfile->Close();
  delete myfile;
  cout<<"Files closed"<<endl;
  }
}

ProcInfo_t meminfo;
gSystem->GetProcInfo(&meminfo);
cout<< "Mem used= "<<meminfo.fMemResident<<"kB"<<endl;;


  TCanvas * c2 = new TCanvas("c2","c2", 1000, 900);
  hEmax->GetXaxis()->SetTitle("Range [um]");
  hEmax->GetYaxis()->SetTitle("Eloss [kev/um]");
  hEmax->GetXaxis()->SetRange(0,xrange);
  hEmax->GetYaxis()->SetRange(0,yrange);
  hEmax->Draw("colz");
  gStyle->SetOptStat(1111);
/*
TFile *myfile= new TFile("test.root","RECREATE");
myfile->cd();
hEmax->Write();
myfile->Close();
delete myfile;*/


}