//some standard C++ includes
#include <iostream>
#include <stdlib.h>
#include <string>
#include <sstream>
#include <vector>
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
//some ROOT includes
#include "TInterpreter.h"
#include "TROOT.h"
#include "TH1F.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"

//"art" includes (canvas, and gallery)
#include "canvas/Utilities/InputTag.h"
#include "gallery/Event.h"
#include "gallery/ValidHandle.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "canvas/Persistency/Common/FindOne.h"

//"larsoft" object includes
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Track.h"
//convenient for us! let's not bother with art and std namespaces!
using namespace art;
using namespace std;

//I like doing this to not get fooled by underflow/overflow
void ShowUnderOverFlow( TH1* h1){
  h1->SetBinContent(1, h1->GetBinContent(0)+h1->GetBinContent(1));
  h1->SetBinContent(0,0);

  int nbins = h1->GetNbinsX();
  h1->SetBinContent(nbins, h1->GetBinContent(nbins)+h1->GetBinContent(nbins+1));
  h1->SetBinContent(nbins+1,0);
}

void demo_ReadHits() {

  //By default, we hate the stats box! But by default, we forgets to disable it in his ROOT profile stuff...
  gStyle->SetOptStat(0);

    TH1D *myh = new TH1D("myh","pmtrack",1000,0,0.04);
    TH1D *myh2 = new TH1D("myh2","theta",1000,0.0,90.0);
    
    TH2D* h3 = new TH2D("h3","Number of tracks;Energy (GeV);Theta Angle (degrees)",1000,0,0.04,1000,0,90);
    h3->SetMarkerStyle(20);
    h3->SetMarkerSize(0.5);
    h3->SetMarkerColor(2);
    

    
    
  //We specify our files in a list of file names!
  //Note: multiple files allowed. Just separate by comma.
  vector<string> filenames { "/uboone/data/users/abhat/Supernova/single_gen_uboone_20171130T154525_g4_20171130T204628_detsim_20171130T221618_reco1_20171130T222724_reco2.root"};
    
  //We need to specify the "input tag" for our collection of optical hits.
  //This is like the module label, except it can also include process name
  //and an instance label. Format is like this:
  //InputTag mytag{ "module_label","instance_label","process_name"};
  //You can ignore instance label if there isn't one. If multiple processes
  //used the same module label, the most recent one should be used by default.
  //
  //Check the contents of your file by setting up a version of larsoft, and
  //running an event dump:
  //  'lar -c eventdump.fcl -s MyInputFile_1.root -n 1 | grep "std::vector<recob::Hit>" '
  InputTag hit_tag { "gaushit" };
  // InputTag track_tag {"pandoraCosmicKalmanTrack"};
  InputTag track_tag {"pmtrack"};
  InputTag mctruth_tag{"generator"};
  //ok, now for the event loop! Here's how it works.
  //
  //gallery has these built-in iterator things.
  //
  //You declare an event with a list of file names. Then, you
  //move to the next event by using the "next()" function.
  //Do that until you are "atEnd()".
  //
  //In a for loop, that looks like this:

    
    int totalevents=0;
    int totaltracks=0;
    
    
    
    
    for (gallery::Event ev(filenames) ; !ev.atEnd(); ev.next()) {

    //to get run and event info, you use this "eventAuxillary()" object.
    cout << "Processing "
	 << "Run " << ev.eventAuxiliary().run() << ", "
	 << "Event " << ev.eventAuxiliary().event() << endl;

    //Now, we want to get a "valid handle" (which is like a pointer to our collection")
    //We use auto, cause it's annoying to write out the fill type. But it's like
    //vector<recob::Hit>* object.
    auto const& hit_handle = ev.getValidHandle<vector<recob::Hit>>(hit_tag);
    auto const& track_handle = ev.getValidHandle<vector<recob::Track>>(track_tag);
    auto const& mctruth_handle = ev.getValidHandle<vector<simb::MCTruth>>(mctruth_tag);
    //We can now treat this like a pointer, or dereference it to have it be like a vector.
    //I (Wes) for some reason prefer the latter, so I always like to do ...
    auto const& hit_vec(*hit_handle);
    auto const& track_vec(*track_handle);
    auto const& mctruth_vec(*mctruth_handle);

    //For good measure, print out the number of hits
    cout << "\tThere are " << hit_vec.size() << " Hits in this event." << endl;
    cout << "\tThere are " << track_vec.size() << " Tracks in this event." << endl;
    cout << "\tThere are " << mctruth_vec.size() << " MCTruth objects in this event." << endl;
    
      
      //FOR loop to go over the contents of the vectors and print them out
      for (auto const& c : mctruth_vec){
      
          std::cout << c << ' ';
          
          double myE = c.GetParticle(0).E();//Energy of incoming electron
        
          for(Int_t j=1; j<=track_vec.size(); j++)
              myh->Fill(myE);
          
          cout << "Initial Energy of the particle = "<< myE << ' '<<" GeV"<<endl;

      }
      
      //FOR loop to go over the contents of the vectors and print them out
      for (auto const& c : track_vec){
          
          std::cout << c << ' ';
          double mytheta= c.Theta(); //Theta angle of the track
          for(Int_t j=1; j<=track_vec.size(); j++)
              myh2->Fill(mytheta);
          
          cout << "Theta angle of the track = "<< mytheta << ' '<<" degrees"<<endl;

      }

      
      
      
      //We can fill our histogram for number of hits now!!!
    h_hits_per_ev->Fill(hit_vec.size());
    h_tracks_per_ev->Fill(track_vec.size());

    //We can loop over the vector to get hit info too!
    //
    //Don't remember what's in a recob::Hit? Me neither. So, we can look it up.
    //The environment variable $LARDATAOBJ_DIR contains the directory of the
    //lardataobj product we set up. Inside that directory we can find what we need:
    // ' ls $LARDATAOBJ_DIR/source/lardataobj/RecoBase/Hit.h '
    //Note: it's the same directory structure as the include, after you go into
    //the 'source' directory. Look at that file and see what you can access.
    //
    //So, let's fill our histogram for the hit integral/peak time/peak amplitude.
    //We can use a range-based for loop for ease.
   // for( auto const& hit : hit_vec){
   //   h_hit_integral->Fill(hit.Integral());
   //   h_hit_peaktime->Fill(hit.PeakTime());
  //    h_hit_peakamp->Fill(hit.PeakAmplitude());
  }
      totaltracks+=track_vec.size();
      totalevents++;
      
    
  } //end loop over events!
    //TFile *fileout = new TFile("mytest.root","RECREATE");
  //  myh->GetXaxis()->SetTitle("Incoming Energy [GeV]");
  //  myh->GetYaxis()->SetTitle("No. of tracks/MeV");
  //  myh->Draw();
  //  myh->Write();
   // fileout->Write();
   // fileout->Close();
    
    cout<<"\tThere are "<<totalevents<<" events in this root file"<<endl;
    cout<<"\tThere are "<<totaltracks<<" tracks in this root file"<<endl;

  //now, we're in a macro: we can just draw the histogram!
  //Let's make a TCanvas to draw our two histograms side-by-side
  TCanvas* canvas = new TCanvas("canvas","Hit Info!",1500,500);
    canvas->Divide(3); //divides the canvas in three!
   canvas->cd(1); //moves us to the first canvas
    myh->GetXaxis()->SetTitle("Incoming Energy [GeV]");
    myh->GetYaxis()->SetTitle("No. of tracks/MeV");
    myh->Draw();

 // h_hits_per_ev->Draw();
   canvas->cd(2);     //moves us to the second canvas
    myh2->GetXaxis()->SetTitle("Theta Angle [degrees]");
    myh2->GetYaxis()->SetTitle("No. of tracks");
    myh2->Draw();

    
    // canvas->cd(3);     //moves us to the third canvas

  //and ... done!
}