#include "TFile.h"
#include "TTree.h"
#include "TBrowser.h"
#include "TH2.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TROOT.h"
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>

// per eseguire questa macro:
// root> .L generateTreePMT.C++
// root> generateTreePMT() 
// NOTA il file di input si deve chiamare inputFile.dat
// e' conveniente creare un link logico del file da processare denominato inputFile.dat (ln -sf xxx.txt inputData.dat)
// Il file di output e' chiamato treePMT.root

// This example illustrates how to make a Tree from variables or arrays
// in a C struct. Use of C structs is strongly discouraged and one should
// use classes instead. However support for C structs is important for
// legacy applications written in C or Fortran.
//    see tree2a.C for the same example using a class instead of a C-struct.
//
// In this example, we are mapping a C struct to one of the Geant3
// common blocks /gctrak/. In the real life, this common will be filled
// by Geant3 at each step and only the Tree Fill function should be called.
// The example emulates the Geant3 step routines.
//
// to run the example, do:
// .x tree2.C   to execute with the CINT interpreter
// .x tree2.C++ to execute with native compiler
//
//  Author: Rene Brun

//const Int_t MAXMEC = 30;
//      PARAMETER (MAXMEC=30)
//      COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(MAXMEC)
//     + ,NAMEC(MAXMEC),NSTEP ,PID,DESTEP,DESTEL,SAFETY,SLENG
//     + ,STEP  ,SNEXT ,SFIELD,TOFG  ,GEKRAT,UPWGHT
typedef struct {
  Int_t         eventNumber;
  Int_t         eventCounter;
  Float_t       ypeak       ;
  Float_t       peakint     ;
  Float_t       risetime    ;
  Float_t       falltime    ;
  Float_t       pulseheight ;
  Float_t       yintegral   ;
} dataPMT_t;



void tree2w()
{
  //create a Tree file tree2.root

  std::cout<<"In Tree2w"<<std::endl;
  //create the file, the Tree and a few branches with
  //a subset of gctrak
  TFile f("treePMT.root","recreate");
  TTree t("t","a Tree for PMT studies");
  dataPMT_t event;


  t.Branch("eventnumber" ,&event.eventNumber,"eventnumber/I");
  t.Branch("eventcounter",&event.eventCounter,"eventcounter/I");
  t.Branch("ypeak"       ,&event.ypeak       ,"ypeak/F");
  t.Branch("peakint"     ,&event.peakint     ,"peakint/F");
  t.Branch("risetime"    ,&event.risetime    ,"risetime/F");
  t.Branch("falltime"    ,&event.falltime    ,"falltime/F");
  t.Branch("pulseheight" ,&event.pulseheight ,"pulseheight/F");
  t.Branch("yintegral"   ,&event.yintegral   ,"yintegral/F");


  //Initialize event variables 
  event.eventNumber =0;
  event.eventCounter=0;
  event.ypeak       =0.;
  event.peakint     =0.;
  event.risetime    =0.;
  event.falltime    =0.;
  event.pulseheight =0.;
  event.yintegral   =0.;


  //std::string fileNameString = "2016-04-22-15-29-38-HVscan_R9420snDD3851_HV800_LA7_45LB60_500nsec_GrIPE.txt";
  std::string fileNameString = "inputFile.dat";
  ifstream inputFile; //abilita la lettura da file inputFile è un comando per caricare il file

   inputFile.open(fileNameString.c_str());
   if (!inputFile.good()) {cout<<"errore nell'apertura del file "<<endl; return;} 
   else {std::cout<<"Input file is good !!!"<<std::endl;}

   std::string aLine;
   Int_t evc = 0;
   Int_t evn = 0; Float_t yp = 0.; Float_t pi = 0; Float_t rt = 0.; Float_t ft = 0.; Float_t ph = 0; Float_t yi = 0.;
   // skip 7 lines at the start (header)
   for (int i=0;i<7;++i) getline(inputFile, aLine);
   std::cout<<"Header is skipped !!!"<<std::endl;
   while (! inputFile.eof() ) 
     {
       std::cout<<" not yeat at the end of file"<<std::endl;
       getline(inputFile, aLine);
       istringstream aLineStream(aLine);
       std::cout<<"READ <"<<aLine<<">"<<std::endl;
       aLineStream >> evn >> yp >> pi >> rt >> ft >> ph >> yi;
       //Float_t Charge = (yi / 50.)*0.5E-9; // Coulomb;  ( V/50 Ohm ) x 0.5ns (bin width)  
       //Charge = Charge/1.6021766208E-19;   // n. of electrons
       //Charge = Charge / 1.E6;             // Me (Millions of electrons)
       Float_t Charge =  -(yi / 50.)*0.5*10000./1.6021766208;    //  Me (Millions of electrons) 
       ++evc;
       std::cout<<"n. of lines read so far is  "<<evc<<std::endl;
       event.eventCounter=evc;
       event.eventNumber =evn;
       event.ypeak       =yp;
       event.peakint     =pi;
       event.risetime    =rt;
       event.falltime    =ft;
       event.pulseheight =ph; //mV 
       event.yintegral   =Charge ; 
       t.Fill();
   }

  //save the Tree header. The file will be automatically closed
  //when going out of the function scope
  t.Write();
}

void tree2r()
{
  //read the Tree generated by tree2w and fill one histogram
  //we are only interested by the destep branch.

  //note that we use "new" to create the TFile and TTree objects !
  //because we want to keep these objects alive when we leave
  //this function.
  TFile *f = new TFile("treePMT_example.root");
  TTree *t = (TTree*)f->Get("t");
  static Float_t charge;
  TBranch *b_yint = t->GetBranch("yintegral");
  b_yint->SetAddress(&charge);

  //create one histogram
  Float_t minCh = 100.;
  Float_t maxCh = 200.;
  Int_t nbins = 100;
  TH1F *hyint   = new TH1F("hyint","Charge [Me]",nbins,minCh,maxCh);

  //read only the destep branch for all entries
  Long64_t nentries = t->GetEntries();
  for (Long64_t i=0;i<nentries;i++) {
    b_yint->GetEntry(i);
    std::cout<<"printing CHARGE [ Me ] = "<<charge<<std::endl;
    hyint->Fill(charge);
  }

  //we do not close the file.
  //We want to keep the generated histograms
  //We fill a 3-d scatter plot with the particle step coordinates
  TCanvas *c1 = new TCanvas("c1","c1",600,800);
  c1->SetFillColor(0);
  hyint->Draw();
  if (gROOT->IsBatch()) return;

  // invoke the x3d viewer
  gPad->GetViewer3D();
}

void generateTreePMT() {
  tree2w();
  tree2r();
}