#include <TH1D.h>
#include <TH2D.h>
#include <TFile.h>
#include <TMath.h>
#include <TCanvas.h>
// -- RAT headers
#include <RAT/DU/DSReader.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/MCPE.hh>
#include <RAT/DS/MCPMT.hh>
#include <RAT/DU/LightPathCalculator.hh>

#include <iostream>
#include <vector>
#include <string>
#include <bitset>

using std::cout;
using std::endl;
using std::vector;
using std::string;


void draw_history(const char * input_file) {


  TH1D *phhist1 = new TH1D("phhist1","PE history (accumulated)",16,0,16);
  TH1D *phhist2 = new TH1D("phhist2","PE history (accumulated)",16,0,16);
  const char* labels[] = {"Unknown","Cherenkov","Scintillation","Rayleigh","Mie","WLS","WLS_C0","WLS_C1","WLS_C2","WLS_C3","WLS_C4","WLS_C5"};
  // If this is being done on data that does not require remote database connection
  // eg.: a simple simulation with default run number (0)
  // We can disable the remote connections:
  //
  // NOTE: Don't do this if you are using real data!!!
  RAT::DB::Get()->SetAirplaneModeStatus(true);

  RAT::DU::DSReader ds(input_file);

  cout << "Loaded DS with " << ds.GetEntryCount() << " entries" << endl;
  cout << "This file contains " << ds.GetRunCount() << " runs in it" << endl;


  size_t n_evts = ds.GetEntryCount();
  for (size_t i_ds = 0; i_ds < n_evts; ++i_ds) {
    const RAT::DS::Entry& rDS = ds.GetEntry(i_ds);
    if (!rDS.MCExists()) {
      cout << " ERROR: There is no MC data" << endl;
      return;
    }
    const RAT::DS::MC &rMC = rDS.GetMC();
    // -- NFB: For later
    cout << "Entry " << i_ds << " has " << rDS.GetMCEVCount() << " MCEVs" << endl;
    size_t n_mcpmt = rMC.GetMCPMTCount();

    for (size_t i_mcpmt = 0; i_mcpmt < n_mcpmt; ++i_mcpmt) {
      const RAT::DS::MCPMT &rMCPMT = rMC.GetMCPMT(i_mcpmt);
      size_t n_mcpe = rMCPMT.GetMCPECount();
      for (size_t i_mcpe = 0; i_mcpe < n_mcpe; ++i_mcpe) {
        const RAT::DS::MCPE &rMCPE = rMCPMT.GetMCPE(i_mcpe);
        // -- First method of loading and using the photon history
        UShort_t hist = rMCPE.GetHistory();

        cout << "History : " << hist << " hex " << std::hex << hist
             << std::dec << " bin " << std::bitset<16>(hist) << endl;

        for (unsigned int bit = 0; bit < 16; ++bit) {
          if (hist & (0x1 << bit)) {
            phhist1->Fill(labels[bit],1);
          }
        }

        // -- Method 2 of using the photon history...asking specifically for each bit
        for (unsigned int bit = 0; bit < 11; ++bit) {
          RAT::DS::MCPE::PhotonHistory hbit = static_cast<RAT::DS::MCPE::PhotonHistory>(bit);
          if (rMCPE.GetFromHistory(hbit)) {
            phhist2->Fill(labels[bit],1);
          }
        }
      }
    }
  }

  // -- Make a canvas and draw
  TCanvas *c1 = new TCanvas();
  c1->cd();
  c1->Divide(1,2);
  c1->cd(1);
  phhist1->Draw();
  c1->cd(2);
  phhist2->Draw();

}