//--------------------------------------------------------------------------------
// skeleton.cc
//
// This is a skeleton program to read/analyse TRTNtuple files via standalnoe c++
// If you have installed root, you can compile this file as follows:
// g++ skeleton.cc -o skeleton.exe `root-config --cflags --glibs`
//
// version 1.0
// first written on 23 Sep 2022
//--------------------------------------------------------------------------------
 
// C++ headers
#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <iomanip>
using namespace std;
 
// ROOT headers
#include "TApplication.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TF1.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TNtuple.h"
#include "TStyle.h"
 
// MakeFile headers
#include "../Sina/TRT.C"
//#include "../Sina/TRT.h"
 
int main(int argc, char** argv){
  // to start application
  TApplication *myApp = new TApplication("myapp",&argc, argv);
 
  // determine the file type MC=1 or data=0
  bool isSimulation = false;
 
  if(isSimulation) cout << "Running Monte Carlo Simulation File" << endl;
  else             cout << "Running Recorded Data File" << endl;
 
  // output root file where histograms are written
  TFile *fout;
  if(isSimulation) fout = new TFile("histoMC.root","recreate");
  else             fout = new TFile("histoDA.root","recreate");
 
  // book histograms
  TH1F *AvgMu      = new TH1F("AvgMu", ";< #mu >;",                 71, 0.5,71.5);
  TH1F *BeamSpot_x = new TH1F("BeamSpot_x", ";beam spot x(mm);",   100, -5,5);
  TH1F *InDet_pt   = new TH1F("InDet_pt",";InDet P_t (GeV)",       100, 1,10);
  TH1F *DriftTime  = new TH1F("DriftTime",";DriftTime (ns)",        60, 0,60);
  TH1F *Straw_r    = new TH1F("Straw_r",";;Straw r ?",             100, 0,2);
 
  // chains to read tree in data files
  TChain *chSimu = new TChain("TRT","TRT"); // MC
  TChain *chData = new TChain("TRT","TRT"); // Data
 
  // add files to chains
  if(isSimulation) chSimu->Add("../Sina/xTRTNTuple_Zmm.root"); // MC
  else             chData->Add("../Sina/xTRTNTuple_Zmm.root"); // Data
 
  // define TRT sample pointer
  TRT *TRTSample;
  if(isSimulation) TRTSample = new TRT(chSimu); // MC
  else             TRTSample = new TRT(chData); // Data
 
  // divisor for total number tracks to process
  int eventdivisor = 1;
  Long64_t nTracks;
  if(isSimulation) nTracks = chSimu->GetEntries() / eventdivisor;
  else             nTracks = chData->GetEntries() / eventdivisor;
  cout << "Total events desired = " << nTracks << " ( " << 100 / eventdivisor << "% ) "<< endl;
 
  // start of main event loop ------------------------------------------------------
  for(Long64_t trk = 0; trk < nTracks; trk++)
  {
    //get entries in the event (either MC or DA)
    TRTSample->GetEntry(trk);
 
    // this is a scalar variable
    AvgMu->Fill(TRTSample->averageInteractionsPerCrossing);
 
    // this is a TVector3 variable
    BeamSpot_x->Fill( TRTSample->beamSpot->x() );
 
    // this is a TLorentzVector variable
    InDet_pt->Fill( TRTSample->InDetTrackParticles_p4->Pt()/1000.0 );
 
    // this is a vector<float> variable
    for(unsigned int i =0; i<TRTSample->InDetTrackParticles_hit_driftTime->size(); i++){
      double t = TRTSample->InDetTrackParticles_hit_driftTime->at(i);
      DriftTime->Fill(t);
    }
 
    // this is a vector<Vector3> variable
    for(unsigned int i =0; i<TRTSample->InDetTrackParticles_hit_strawPosition->size(); i++){
      double r = TRTSample->InDetTrackParticles_hit_strawPosition->at(i).Pt() / 1000.0;
      Straw_r->Fill(r);
    }
 
  }
  //end of event main loop
  // ------------------------------------------------------------------------
 
  // You can either
  // create a canvas and draw histograms on it
  TCanvas *c1 = new TCanvas("c1","c1",1200,800);
  c1->Divide(2,3);
  c1->cd(1); AvgMu->Draw();
  c1->cd(2); BeamSpot_x->Draw();
  c1->cd(3); InDet_pt->Draw();
  c1->cd(4); DriftTime->Draw();
  c1->cd(5); Straw_r->Draw();
  c1->Print("output.png");
 
  // or write all histograms to an output histogram file defined above
  fout->Write();
  fout->Close();
 
  return 0;
  myApp->Run();
}
// EOF