#include <fstream>
#include <iostream>
#include <string>

using namespace std;


void fillHistWithFS ( TString A_File, TH2F * A_hist, Int_t A_verboseLevel = 0 )
{
  
  TFile *f = new TFile(  A_File );
  TTree *t = (TTree*)f->Get("cbmsim");
  
  Int_t nevent = t->GetEntries();
  
  FS * ffs;
  t->SetBranchAddress( "FASTSLOW", &ffs);

  if( A_verboseLevel > 0 ) cout << "No Entries= " << nevent << endl;
  
  for (Int_t i=0; i<nevent; i++) {

    t-> GetEntry(i);  

    Int_t fMulti = ffs->multi; 
    if( A_verboseLevel > 0 ) cout << "Event id = " << i << endl;
    if( A_verboseLevel > 0 ) cout << "Multiplicity = " << fMulti << endl;
    for (Int_t pi = 0; pi < fMulti; pi++ )
    {
      Double_t fFast = ffs->fast[pi];
      Double_t fSlow = ffs->slow[pi];
      
      if( A_verboseLevel > 1)
      {
        cout << "  Fast=\t" << fFast << endl;
        cout << "  Slow=\t" << fSlow << endl;
        cout << "  NumTel=\t" <<ffs->numtel[pi] << endl;
        cout << "  De=\t" <<ffs->de[pi] << endl;
        cout << "  Ekin=\t" <<ffs->ekin[pi] << endl;
      }
      A_hist->Fill( fSlow, fFast );
    }
  }
  f->Close();
  delete f;
  return;
}

void read_output_0911 (  ) {
  
  
  TString dir = getenv("VMCWORKDIR");
  TString r3bdir = dir + "/macros";

  TString r3b_geomdir = dir + "/geometry";
  gSystem->Setenv("GEOMPATH",r3b_geomdir.Data());

  TString r3b_confdir = dir + "gconfig";
  gSystem->Setenv("CONFIG_DIR",r3b_confdir.Data());

  TString inFile1 = "./proton/output_p.root";
  TString inFile2 = "./alpha/output_alpha.root";
  
  // In general, the following parts need not be touched
  // ========================================================================

  // ----    Debug option   -------------------------------------------------
  gDebug = 0;
  // ------------------------------------------------------------------------

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  // ----  Load libraries   -------------------------------------------------
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();

  gSystem->Load("libGenVector");

  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libMCStack");
  gSystem->Load("libField");
  gSystem->Load("libGen");

  //----  Load R3B specific libraries ---------------------------------------
  gSystem->Load("libR3Bbase");
  gSystem->Load("libR3BGen");
  gSystem->Load("libR3BPassive");
  gSystem->Load("libR3BData");
  gSystem->Load("libR3BCal");
  gSystem->Load("libR3BCalo");
  gSystem->Load("libR3BDch");
  gSystem->Load("libR3BGfi");
  gSystem->Load("libR3BLand");
  gSystem->Load("libR3BmTof");
  gSystem->Load("libR3BTof");
  gSystem->Load("libR3BATof");
  gSystem->Load("libR3BTra");
  gSystem->Load("libR3BChimera");
  gSystem->Load("libELILuMon");
  gSystem->Load("libData");
  gSystem->Load("libKratta"); 
  gSystem->Load("libKrattaTasks"); 
  
  
  TH2F * hist1 = new TH2F("hist1", "Proton", 1000, 0, 1000,1000, 0, 1000);
  TH2F * hist2 = new TH2F("hist2", "Alpha", 1000, 0, 1000,1000, 0, 1000);
  hist1->SetXTitle("slow");
  hist1->SetYTitle("fast");  
 
 
  fillHistWithFS ( inFile1, hist1 ) ;
  fillHistWithFS ( inFile2, hist2 ) ;
  
  gStyle->SetPalette(1);
  hist1->SetTitle("Fast vs slow for protons and alphas");
  
  hist1->Draw("COLSCAT");
  hist2->Draw("COLSCATSAME");
  
  cout << "END" << endl;

}