// ROOT includes
#include "TH1.h"
#include "TFile.h"
#include "TChain.h"
#include "TTree.h"
#include "TMath.h"

//C++ includes
#include <iostream>
#include <iomanip>

void analisis(){
  /////////////////////////////////////////////////////////////////////////
  // 
  // Preparacion de muestras, histogramas, ramas, ...
 
  cout << "------------------------------------------" << endl;
  cout << endl;
  
  //Creemos el TChain que nos permite manejar la muestra
  cout << ">> Creando la muestra..." << endl;
  TChain *muestra = new TChain("MVATree");
  //+ An~adimos uno a uno los ficheros  
  muestra->Add("MVA_SingleTop_PU.root");
  
  
  //Averiguemos el numero de sucesos en la muestra
  cout << ">> Calculando el numero de sucesos..." << endl;
  int nEvents = muestra->GetEntries();
  cout << endl;
  cout << "------------------------------------------" << endl;
  cout << nEvents << " sucesos totales (sin peso)" << endl;
  cout << "------------------------------------------" << endl;
  
  
  //Leemos las ramas del ttree
  cout << ">> Preparando las ramas... " << endl;

  // + Channel 
  int Channel;
  muestra->SetBranchAddress("Channel", &Channel);  

  double dPhi;
  muestra->SetBranchAddress("dPhi", &dPhi);  

  
  //Creamos el fichero de salida
  cout << ">> Creando el fichero de salida ..." << endl;
  TString histoFileName("histo.root");
  TFile* histoAnalisis = TFile::Open(histoFileName, "RECREATE");

  // + comprobamos que lo hemos abierto adecuadamente
  if (histoAnalisis->IsOpen())
    cout << "   ... \"" << histoFileName << "\" creado" << endl;
  else {
    cout << "   ... ERROR: \""<< histoFileName 
	 << "\" ya existe. Borralo." << endl;
    return;
  }
  
  
  
  //Creamos los histogramas
  cout << ">> Creando histogramas..." << endl;
  TH1F *hchannel = new TH1F("hChannel", "Channel", 5, -0.5, 5);
  TH1F *hdphi = new TH1F("hdPhi", "#Delta #phi", 100, 0, TMath::Pi());
  
  //
  /////////////////////////////////////////////////////////////////////////
 
  
  /////////////////////////////////////////////////////////////////////////
  //  
  // Bucle a sucesos (Event loop)
  //
  //
  cout << ">> Iterando sobre sucesos..." << endl;
  // + Un truquillo para el output
  unsigned int ndigits=log10(abs(nEvents))+1;
  char* strdelete= new char[ndigits+1];

  for (unsigned int i=0; i<ndigits;i++)
    strdelete[i]='\b';
  strdelete[ndigits]='\0';

  // + Aqui comienza el bucle
  for(int event = 0; event<nEvents; event++){

    //+ Mostrar cada 100 sucesos donde estamos
    if ((event % 100) == 0)
      cout << strdelete << setw(ndigits) << setfill('-') << event << flush;

    //+ Cargamos el suceso en las variables 
    muestra->GetEntry(event);
    
    //+ Llenamos los histogramas
    hchannel->Fill(Channel);
    if (Channel == 0)
      hdphi->Fill(dPhi);
    
  }
  //Final loop output
  cout << strdelete << setw(ndigits) << setfill('-') << event << flush;
  
  //
  /////////////////////////////////////////////////////////////////////////
 
  /////////////////////////////////////////////////////////////////////////
  // Asuntos finales
  cout << endl << ">> Mostrando los histogramas..." << endl;
  //hchannel->Draw();
  hdphi->Draw();


  cout << endl << ">> Guardando los histogramas..." << endl;
  histoAnalisis->Write();
  //histoAnalisis->Close();
  
  
  cout << "------------------------------------------" << endl;
  cout << endl;
  
  //
  /////////////////////////////////////////////////////////////////////////
}