/*************************************************************
 * 
 * demo_ReadEvent program
 * 
 * This is a simple demonstration of reading a LArSoft file 
 * and printing out the run and event numbers. You can also
 * put the event numbers into a histogram!
 *
 * Original Author: Wesley Ketchum (wketchum@fnal.gov), Aug28, 2016
 * 
 *************************************************************/


//some standard C++ includes
#include <iostream>
#include <stdlib.h>
#include <string>
#include <vector>

//some ROOT includes
#include "TInterpreter.h"
#include "TROOT.h"
#include "TH1F.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TStopwatch.h"
#include "TTree.h"
#include "TBranch.h"
#include "TVector3.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"

#include "lardataobj/RecoBase/Track.h"

//convenient for us! let's not bother with art and std namespaces!
using namespace art;
using namespace std;

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_ReadEvent(){
 
    gStyle->SetOptStat(0);
    
    TFile f_output("MCTracktree_gallery_output.root","RECREATE");
    
    
    
  //Let's make a histogram to store event numbers.
  //I ran this before, so I know my event range. You can adjust this for your file!

  //note, because I'm in my standalone code now, I'm not going to make this a pointer
  //so I can have nice clean memory
  
    
    
    int mypdg;
    double StartX;
    double StartY;
    double StartZ;
    double EndX;
    double EndY;
    double EndZ;
    
    
    double StartX_det;
    double StartY_det;
    double StartZ_det;
    double EndX_det;
    double EndY_det;
    double EndZ_det;
    
    
    
    double tracklength;
    double tracklength_det;
    
    
    
    double muonStartX;
    double muonStartY;
    double muonStartZ;
    double muonEndX;
    double muonEndY;
    double muonEndZ;
    
    
    
    double muonStartX_det;
    double muonStartY_det;
    double muonStartZ_det;
    double muonEndX_det;
    double muonEndY_det;
    double muonEndZ_det;
    
    
    
    double muontracklength;
    double muontracklength_det;
    
    
    
    TTree* MCTracktree = new TTree("MCTracktree","MCTracktree");
    
    MCTracktree->Branch("mypdg",&mypdg,"mypdg/D");
    MCTracktree->Branch("StartX",&StartX,"StartX/D");
    MCTracktree->Branch("StartY",&StartY,"StartY/D");
    MCTracktree->Branch("StartZ",&StartZ,"StartZ/D");
    MCTracktree->Branch("EndX",&EndX,"EndX/D");
    MCTracktree->Branch("EndY",&EndY,"EndY/D");
    MCTracktree->Branch("EndZ",&EndZ,"EndZ/D");
    
    MCTracktree->Branch("mypdg",&mypdg,"mypdg/D");
    MCTracktree->Branch("StartX_det",&StartX_det,"StartX_det/D");
    MCTracktree->Branch("StartY_det",&StartY_det,"StartY_det/D");
    MCTracktree->Branch("StartZ_det",&StartZ_det,"StartZ_det/D");
    MCTracktree->Branch("EndX_det",&EndX_det,"EndX_det/D");
    MCTracktree->Branch("EndY_det",&EndY_det,"EndY_det/D");
    MCTracktree->Branch("EndZ_det",&EndZ_det,"EndZ_det/D");

    
    MCTracktree->Branch("muonStartX",&muonStartX,"muonStartX/D");
    MCTracktree->Branch("muonStartY",&muonStartY,"muonStartY/D");
    MCTracktree->Branch("muonStartZ",&muonStartZ,"muonStartZ/D");
    MCTracktree->Branch("muonEndX",&muonEndX,"muonEndX/D");
    MCTracktree->Branch("muonEndY",&muonEndY,"muonEndY/D");
    MCTracktree->Branch("muonEndZ",&muonEndZ,"muonEndZ/D");

    
    
    MCTracktree->Branch("muonStartX_det",&muonStartX_det,"muonStartX_det/D");
    MCTracktree->Branch("muonStartY_det",&muonStartY_det,"muonStartY_det/D");
    MCTracktree->Branch("muonStartZ_det",&muonStartZ_det,"muonStartZ_det/D");
    MCTracktree->Branch("muonEndX_det",&muonEndX_det,"muonEndX_det/D");
    MCTracktree->Branch("muonEndY_det",&muonEndY_det,"muonEndY_det/D");
    MCTracktree->Branch("muonEndZ_det",&muonEndZ_det,"muonEndZ_det/D");
    
    
    MCTracktree->Branch("muontracklength_det",&muontracklength_det,"muontracklength_det/D");
    MCTracktree->Branch("tracklength_det",&tracklength_det,"tracklength_det/D");
    
    MCTracktree->Branch("muontracklength_det",&muontracklength_det,"muontracklength_det/D");
    MCTracktree->Branch("tracklength_det",&tracklength_det,"tracklength_det/D");
    

    
    
    
    
    
    
    
    
    
    
    
  
  //We specify our files in a list of file names!
  //Note: multiple files allowed. Just separate by comma.
  //  vector<string> filenames { "MyInputFile_1.root" };
  vector<string> filenames { "/pnfs/uboone/mc/uboone/reconstructed/prod_v06_26_06/prodcosmics_corsika_cmc_uboone_intime_reco/reco/prodcosmics_corsika_cmc_uboone_intime_9_20180106T234834_gen2_b68094ab-d48a-42de-84f1-274a79047091_20180125T160523_reco1_20180125T163723_reco2.root","/pnfs/uboone/mc/uboone/reconstructed/prod_v06_26_06/prodcosmics_corsika_cmc_uboone_intime_reco/reco/prodcosmics_corsika_cmc_uboone_intime_0_20171228T135151_gen2_ad665c20-75db-467b-80fc-d945e794f389_20180111T075909_reco1_20180111T082044_reco2.root" };
  
  InputTag mctrack_tag{"mcreco"};
 // InputTag mcstep_tag{"mcreco"};
  //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:
  
  Double_t x1[1000],y1[1000],x2[1000],y2[1000];
    
    
    
    
    
  int counter=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;

    // Track Handle Information in the loop.
    
    auto const& mctrack_handle = ev.getValidHandle<vector<sim::MCTrack>>(mctrack_tag);
      
    auto const& mctrack_vec(*mctrack_handle);
 
    counter++;
      
  //  cout << "\tThere are " <<mctrack_vec.size()<< " MCTrack tracks in this event." << endl;
  
  // h_MCTrack_length->Fill(mctrack_vec.size());
      
     
      
      
      x1[counter]=counter;
      y1[counter]=mctrack_vec.size();
      
      int muontrackcount=0;
    
      
      
      

      
      for (auto const& a : mctrack_vec){
          
           StartX=a.Start().X();
           StartY=a.Start().Y();
           StartZ=a.Start().Z();
           EndX=a.End().X();
           EndY=a.End().Y();
           EndZ=a.End().Z();
           tracklength=sqrt(pow((EndX-StartX),2)+pow((EndY-StartY),2)+pow((EndZ-StartZ),2));

          
          StartX_det=a.at(0).X();
          StartY_det=a.at(0).Y();
          StartZ_det=a.at(0).Z();
          EndX_det=a.at(a.size()-1).X();
          EndY_det=a.at(a.size()-1).Y();
          EndZ_det=a.at(a.size()-1).Z();
          tracklength_det=sqrt(pow((EndX_det-StartX_det),2)+pow((EndY_det-StartY_det),2)+pow((EndZ_det-StartZ_det),2));
          
          
          
          
          
          
          mypdg= a.PdgCode();
       //   cout<<"Start Y is: "<<StartY<<endl;
        //   cout<<"End Y is: "<<EndY<<endl;
          if (mypdg==13 || mypdg==-13){
              
              muontrackcount++;
              
               muonStartX=a.Start().X();
               muonStartY=a.Start().Y();
               muonStartZ=a.Start().Z();
               muonEndX=a.End().X();
               muonEndY=a.End().Y();
               muonEndZ=a.End().Z();
               muontracklength=sqrt(pow((muonEndX-muonStartX),2)+pow((muonEndY-muonStartY),2)+pow((muonEndZ-muonStartZ),2));
          
              muonStartX_det=a.at(0).X();
              muonStartY_det=a.at(0).Y();
              muonStartZ_det=a.at(0).Z();
              muonEndX_det=a.at(a.size()-1).X();
              muonEndY_det=a.at(a.size()-1).Y();
              muonEndZ_det=a.at(a.size()-1).Z();
              muontracklength_det=sqrt(pow((muonEndX_det-muonStartX_det),2)+pow((muonEndY_det-muonStartY_det),2)+pow((muonEndZ_det-muonStartZ_det),2));
          
              
              
              
              
    //   MCTracktree->Fill();
          }
          
          x2[counter]=counter;
          y2[counter]=muontrackcount;
         
          
          
          MCTracktree->Fill();
          
      }
      
      
     
      
      
      
    
  } //end loop over events!
  cout<<"Total Events: "<<counter<<endl;

 // ShowUnderOverFlow(h_MCTracks);
 
    TCanvas* canvas = new TCanvas("canvas","Hit Info!",1500,500);
    canvas->Divide(2);
    canvas->cd(1);
    TGraph* gr1 = new TGraph(counter,x1,y1);
    gr1->SetFillColor(kBlue);
    gr1->SetTitle("All MCTracks");
    gr1->GetXaxis()->SetTitle("Event");
    gr1->GetYaxis()->SetTitle("MCTracks per Event");
    gr1->Draw("AB");
    canvas->cd(2);
    TGraph* gr2 = new TGraph(counter,x2,y2);
    gr2->SetFillColor(kRed);
    gr2->GetXaxis()->SetTitle("Event");
    gr2->GetYaxis()->SetTitle("Muon MCTracks per Event");
    gr2->GetYaxis()->SetRangeUser(0,130);
    gr2->SetTitle("Muon MCTracks");
    gr2->Draw("AB");
  
    

   
 
    
   // canvas->Divide(2);
 // canvas->cd(1);
//  h_MCTracks->SetLineColor(kRed);
//  h_MCTracks->Draw();
//  canvas->cd(2);
 // MCTracks_muon->SetLineColor(kBlue);
 // MCTracks_muon->Draw();
  




  //and ... write to file!
   // TFile f_output("demo_ReadEvent_output.root","RECREATE");
   //h_events.Write();
   //f_output.Close();
 
/*
    TCanvas* canvas = new TCanvas("canvas","Hit Info!",1500,500);
    TPad *pad = new TPad("pad","",0,0,1,1);
    pad->Draw();
    pad->cd();
    
    // draw a frame to define the range
    TH1F *hr = pad->DrawFrame(-0.4,0,1.2,12);
    hr->SetXTitle("X title");
    hr->SetYTitle("Y title");
    pad->GetFrame()->SetBorderSize(12);
    
    // draw first graph
    TGraph* gr1 = new TGraph(counter,x1,y1);
    gr1->SetFillColor(kBlue);
    gr1->SetTitle("All MCTracks");
    gr1->GetXaxis()->SetTitle("Event");
    gr1->GetYaxis()->SetTitle("MCTracks per Event");
    gr1->Draw("AB");
    
    
    // draw second graph
    canvas->cd();
    TPad *overlay = new TPad("overlay","",0,0,1,1);
    overlay->SetFillStyle(4000);
    overlay->SetFillColor(0);
    overlay->SetFrameFillStyle(4000);
    overlay->Draw();
    overlay->cd();
    TGraph* gr2 = new TGraph(counter,x2,y2);
    gr2->SetFillColor(kRed);
    gr2->GetXaxis()->SetTitle("Event");
    gr2->GetYaxis()->SetTitle("Muon MCTracks per Event");
    gr2->GetYaxis()->SetRangeUser(0,130);
    gr2->SetTitle("Muon MCTracks");
    
    TH1F *hframe = overlay->DrawFrame(0,0,100,200);
    hframe->GetXaxis()->SetLabelOffset(99);
    hframe->GetYaxis()->SetLabelOffset(99);
    gr2->Draw("AB");
*/
    
    
    

    f_output.Write();
    f_output.Close();
    
}