#include "TRestGeant4Event.h"
#include "TRestGeant4Metadata.h"
#include "TRestTask.h"
#include <TH3D.h>
#include <TMath.h>

#include <chrono>
using namespace std;
using namespace std::chrono;

Double_t rmin=400, rmax=1600;

Int_t Data(){

  gROOT->SetBatch(kTRUE);
  const Int_t EventTypes=22;        //numero de tipos de particulas generadas
  
  //gamma, n mu, vessel, gas
  //Xe-0 to 10-5-Cu-1 (cambiar el nombre de la carpera)
   Double_t EventNumber[EventTypes]=
  //{279,139,35,17,35,  1.74e-3,5.96e-3,8.86e-3  ,3.92e-3,5.72e-3,2.24e-2,1.56e-2,7.53e-2,1.16e-1,6.89e-3,4.94e-2   ,3.15e-4,3.19e-4,5.08e-4,2.11e-4,1.84e-4,1.62e-3} ;
  {326726,163363,40841,20475,40625,  2.04,6.98,10.4,  5,7,26,18,88,136,8,58,  3.69e-1,3.74e-1,5.96e-1,2.48e-1,2.16e-1,2} ;
  
  // Xe-1-5-Cu-10
  //{1349,674,169,85,168,  8.43e-3,2.88e-2,4.29e-2,   1.89e-2,2.76e-2,1.08e-1,7.56e-2,3.64e-1,5.6e-1,3.33e-2,2.39e-1,  1.52e-2,1.54e-2,2.46e-2,1.02e-2,8.92e-3,7.81e-2};
  //{326726,163363,40841,20475,40625,   2,7,10,  5,7,26,18,88,136,8,58,  4,4,6,2,2,19}        ;
  
  //Xe-1-10-Cu-100
  //{4685,2342,586,294,582,  6.05e-2,2.07e-1,3.08e-1,  6.58e-2,9.6e-2,3.76e-2,2.63e-1,1,2,1.16e-1,8.3e-1,  5.29e-1,5.36e-1,8.54e-1,3.55e-1,3.1e-1,3}                           ;
  //{326726,163363,40841,20475,40625,   4,14,21,  5,7,26,18,88,136,8,58,  37,37,60,25,22,189} ;
  
  //Xe-1-10-St-115
  //{4812,2406,601,302,598,       21,1,2,12,1,3,  6.25e-1,6.33e-1,1,4.19e-1,3.66e-1,3}        ;
  //{326726,163363,40841,20475,40625,         1396,68,124,836,74,186,  42,43,69,28,25,218}    ;
  
  //Xe-1-10-Titan-90
  //{4578,2289,572,287,569,       8,3.59e-7,   4.66e-1,4.71e-1,7.51e-1,3.12e-1,2.73e-1,2}                                         ;
  //{326726,163363,40841,20475,40625,    576,2.56e-5,   33,34,54,22,19,170}                         ;
  
  //He-0-5-Cu-1
  //{40,20,5,3,5,   2.52e-4,8.59e-4,1.28e-3,   5.65e-4,8.25e-4,3.23e-3,2.26e-3,1.09e-2,1.67e-2,9.95e-4,7.13e-3   }                                             ;
  //{326726,163363,40841,20475,40625,  2,7,10, 5,7,26,18,88,136,8,58  }                       ;
  
  //Ne-0-5-Cu-1
  //{8,4,1,0.53,1,   5.25e-5,1.79e-4,2.67e-4,    1.18e-4,1.72e-4,6.73e-4,4.71e-4,2.27e-3,3.49e-3,2.08e-4,1.49e-3    }                                               ;
  //{326726,163363,40841,20475,40625,  2,7,10, 5,7,26,18,88,136,8,58  }                       ;
  
  //Xe-2500-5-Cu-10
  //{1349,674,169,85,6.71e-3,  8.43e-3,2.88e-2,4.29e-2,   3.52e-5,4.83e-16,7.13e-3,8.7e-5,2.11e-2,3.69e-4,3.61e-7,1.17e-2,  3.78e-4,1.04e-6,1.76e-3,2.14e-4,8.39e-5,1.41e-9}                              ;
  //{326726,163363,40841,20475,2,   2,7,10,  8.52e-3,1.17e-13,2,2.11e-2,5,8.94e-2,8.76e-5,2.83,  9.16e-2,2.52e-4,4.26e-1,5.19e-2,2.03e-2,3.43e-7};
  

   Int_t Color[EventTypes]={803,793,632,600,430
                            ,416,419,828
                            ,594,875,616,810,892,619,871,886
                            //,892,594,875,616,810,619
                            //,594,828
                            ,623,920,396,392,832,831};
                            
                            
   Char_t Title[EventTypes][40]={
  "#gamma (1<E<10) MeV","#gamma (10<E<100) MeV","#gamma (0.1<E<100) GeV","Neutrones", "Muones",
   "{}^{238}U","{}^{232}Th","{}^{40}K",
   "{}^{46}Sc","{}^{48}V","{}^{54}Mn","{}^{56}Co","{}^{57}Co","{}^{58}Co","{}^{59}Fe","{}^{60}Co",
   //"{}^{7}Be","{}^{46}Sc","{}^{48}V","{}^{54}Mn","{}^{56}Co","{}^{58}Co",
   //"{}^{46}Sc","{}^{40}K",
   "{}^{3}H","{}^{7}Be","{}^{125}Sb","{}^{121}Te","{}^{123}Te","{}^{127}Xe"
   };
   
  Double_t LimTh1=200;  // en KeV
  Double_t BinsTh1=0.5;
  string File="Xe-10-5-Cu-1"; 
  string FileData=File+"data.txt";
  string FilePdfOrig=File+"OrigMod.pdf";
  string FilePdfRed=File+"RedMod.pdf";
  string FilePdfr2=File+"r2Mod.pdf";
  char NameData[100], NameFile1[100], NameFile2[100], NameFile3[100];
  strcpy(NameData,FileData.c_str());
  strcpy(NameFile1,FilePdfOrig.c_str());
  strcpy(NameFile2,FilePdfr2.c_str());
  strcpy(NameFile3,FilePdfRed.c_str());
  FILE *f;
  f=fopen(NameData,"r");
  
  Int_t Tracks;
  Float_t Energy,rho;
  
  TCanvas *c1 = new TCanvas("c1","c1",700,500);
  TCanvas *c2 = new TCanvas("c2","c2",700,500);
  TCanvas *c3 = new TCanvas("c3","c3",700,500);
  THStack *hs= new THStack("hs","; #it{E} [keV]; Eventos keV^{-1} s^{-1}");
  THStack *hsr2= new THStack("hsr2","; #it{r} [m]; Eventos s^{-1}");
  THStack *hsreduce= new THStack("hsreduce","; #it{E} [keV]; Eventos keV^{-1} s^{-1}");
  
  Int_t Exist, Existred, Existr2;
  std::vector<Int_t> Exist2;
  std::vector<Int_t> Exist2red;
  std::vector<Int_t> Exist2r2;
  
  
  
  
  
  Double_t RhoLimit=1500;
  for (int i=0;i<EventTypes;i++){
  
    string aux=std::to_string(i);
    string FileRoot=File+"/Runwater"+aux+".root";
    /*
    TRestRun* run = new TRestRun();
    run->OpenInputFile(FileRoot);
    TRestGeant4Metadata *g4=(TRestGeant4Metadata*)run->GetMetadataClass("TRestGeant4Metadata");
    Double_t Generated=g4->GetNumberOfEvents();      //lanzados en simulacion
    Double_t Detected=run->GetEntries();             //los almacenados
    Double_t GeneratedReal=EventNumber[i]/13;                     //los recibidos en nuestro tiempo deriva
    Double_t DetectedRealOriginal=GeneratedReal*Detected/Generated ;    //los que deberiamos haber detectado
    Double_t DetectedReal=0;

    TRandom3 *r3 = new TRandom3();  r3->SetSeed(0);  DetectedReal=r3->Poisson(DetectedRealOriginal);
    cout << "\n######Eventos lanzados(sim): "<<Generated<<"    Eventos finales(sim): "<<Detected<<endl;
    cout << "\n######Eventos lanzados: "<<GeneratedReal<<"    Eventos finales: "<<DetectedRealOriginal<<"  Tras Poisson: "<<DetectedReal<<endl;
    */

    TH1* h1 = new TH1D("h1",";#it{E} [keV]; Eventos keV^{-1} s^{-1}",BinsTh1*LimTh1,0,LimTh1);
    TH1* h1r2  = new TH1D("h1r2","#it{r} [m]; Eventos s^{-1}",(rmax-rmin)/10 ,rmin/1000,rmax/1000);
    TH1* h1reduce = new TH1D("h1reduce",";#it{E} [keV]; Eventos keV^{-1} s^{-1}",BinsTh1*LimTh1,0,LimTh1);
    h1->SetStats(0);
    h1->SetMarkerColor(Color[i]);
    h1->SetFillColor(Color[i]);
    h1r2->SetStats(0);
    h1r2->SetMarkerColor(Color[i]);
    h1r2->SetFillColor(Color[i]);
    h1reduce->SetStats(0);
    h1reduce->SetMarkerColor(Color[i]);
    h1reduce->SetFillColor(Color[i]);
    Exist=0;  Existred=0;  Existr2=0;
    
    fscanf(f,"%d",&Tracks);
    for (int n = 0; n<Tracks;n++){
        fscanf(f,"%f %f",&Energy,&rho);
        //cout<<Energy<<"\t"<<rho<<endl;
        h1->Fill(Energy);
        h1r2->Fill(rho);
        Existr2=1;
        if (Energy <= LimTh1 ){
           Exist=1;
            }
        if( rho < (RhoLimit/1000) && Energy <= LimTh1){  //FIDUCIALIZACION
           h1reduce->Fill(Energy);
           Existred=1;
           }
    }
    
    h1->Scale(BinsTh1);
    h1reduce->Scale(BinsTh1);
    
    c1->cd();  h1->Draw("HIST BAR");  hs->Add(h1);
    c2->cd();  h1r2->Draw("HIST BAR");   hsr2->Add(h1r2);
    c3->cd();  h1reduce->Draw("HIST BAR");  hsreduce->Add(h1reduce);
    if (Exist != 0 ){  Exist2.push_back(1);}  else Exist2.push_back(0);
    if (Existred != 0 ){  Exist2red.push_back(1);}  else Exist2red.push_back(0);
    if (Existr2 != 0 ){  Exist2r2.push_back(1);}  else Exist2r2.push_back(0);
    
  

  }
  fclose(f);
  
  ////////////// Add neutrino's spectrum
 
  FILE *fneu;
  fneu=fopen("NeutrinoXenontest.txt",
             //"NeutrinoNeon.md",
             //"NeutrinoHelio.md",
             "r");
  
  Float_t x,y;
  TGraph *g=new TGraph();  g->SetMarkerStyle(23); g->SetMarkerColor(kBlack); g->SetLineColor(kBlack); g->SetLineWidth(2);
  while(!feof(fneu)){
    fscanf(fneu,"%f  %f",&x,&y);
    //cout<<x<<"\t"<<y<<endl;
    g->AddPoint(x,y);
  }
  Exist2.push_back(1);
  c1->cd();
  g->Draw("C SAME");
  
  Exist2red.push_back(1);
  c3->cd();
  g->Draw("C SAME");
  
  /*
  TH1* hn = new TH1D("h1",";#it{E} [keV]; Eventos keV^{-1} s^{-1}",BinsTh1*LimTh1,0,LimTh1);
  hn->SetStats(0);
  hn->SetMarkerColor(kBlack);
  hn->SetFillColor(kBlack);
  Float_t x,y;
  while(!feof(fneu)){
    fscanf(fneu,"%f  %f",&x,&y);
    hn->Fill(x,y);
  }
  c1->cd();
  hn->Draw("C SAME");
  hs->Add(hn);
  */
///////////////////////////////
  
  Exist=0;
  Existred=0;
  Existr2=0;
  for (int i=0; i < Exist2.size(); i++)
      Exist+=Exist2[i];          //exist pasa a ser el numero de eventos no nulos total
  for (int i=0; i< Exist2red.size(); i++)
      Existred+=Exist2red[i];
  for (int i=0; i < Exist2r2.size(); i++)
      Existr2+=Exist2r2[i];   
      
  Double_t aux=0.35+(0.023*(EventTypes-Exist-2));
  Double_t auxred=0.35+(0.023*(EventTypes-Existred)); //0.025
  Double_t auxr2=0.35+(0.023*(EventTypes-Existr2));
  TLegend *legend1 = new TLegend(0.80,aux,0.99,0.9); //para 22
  TLegend *legend1red = new TLegend(0.80,auxred,0.99,0.9); 
  TLegend *legend1r2 = new TLegend(0.80,auxr2,0.99,0.9); //para 22
  
  legend1->SetTextSizePixels(23);
  legend1red->SetTextSizePixels(23);
  legend1r2->SetTextSizePixels(23);
  
  legend1->AddEntry((TObject*)0, "", "");  legend1red->AddEntry((TObject*)0, "", "");  legend1r2->AddEntry((TObject*)0, "", "");
  /*
  cout<<Exist2.size()<<endl;
  if(Exist2.size()>EventTypes){      //Señal de neutrinos en la leyenda
    legend1->AddEntry(g,"Neutrino","l");
    legend1red->AddEntry(g,"Neutrino","l");}
  */
  
  for (int i=0;i <EventTypes;i++){
    if (Exist2[i]==1. ){
      TH1* h1aux = new TH1D("h1aux","; ; ",1,0,2);
      h1aux->SetMarkerColor(Color[i]);
      h1aux->SetFillColor(Color[i]);
      legend1->AddEntry(h1aux,Title[i],"f");
      }
    if (Exist2red[i]==1. ){
      TH1* h1auxred = new TH1D("h1auxred","; ; ",1,0,2);
      h1auxred->SetMarkerColor(Color[i]);
      h1auxred->SetFillColor(Color[i]);
      legend1red->AddEntry(h1auxred,Title[i],"f");
      }
    if (Exist2r2[i]==1. ){
      TH1* h1auxr2 = new TH1D("h1auxr2","; ; ",1,0,2);
      h1auxr2->SetMarkerColor(Color[i]);
      h1auxr2->SetFillColor(Color[i]);
      legend1r2->AddEntry(h1auxr2,Title[i],"f");
      }
  }
  legend1->AddEntry((TObject*)0, "", "");  legend1red->AddEntry((TObject*)0, "", "");  legend1r2->AddEntry((TObject*)0, "", "");
  
  
  c1->cd();
  TH1D *last=(TH1D*)hs->GetStack()->Last();  hs->SetMaximum((1.25)*(last->GetBinContent(last->GetMaximumBin())));  hs->Draw("HIST BAR");
  gPad->RedrawAxis("Y");
  legend1->Draw();
  
  c2->cd();
  TH1D *lastr2=(TH1D*)hsr2->GetStack()->Last();  hsr2->SetMaximum((1.25)*(lastr2->GetBinContent(lastr2->GetMaximumBin())));  hsr2->Draw("HIST BAR");
  gPad->RedrawAxis("Y");
  legend1r2->Draw();
  
  c3->cd();
  TH1D *lastreduce=(TH1D*)hsreduce->GetStack()->Last();  hsreduce->SetMaximum((1.25)*(last->GetBinContent(last->GetMaximumBin())));  hsreduce->Draw("HIST BAR");
  gPad->RedrawAxis("Y");
  legend1red->Draw();
  
  //last->Draw("HIST BAR");
  

  c1->cd();  c1->SaveAs(NameFile1);  c1->Close();
  //c2->cd();  c2->SaveAs(NameFile2);  c2->Close();
  //c3->cd();  c3->SaveAs(NameFile3);  c3->Close();
  
  
  
  
  
   

return 0;
}