///////////////////////////////////////////////////////////////////////////////////////
// Functions demonstrating the use of the TrackNav, TrackCursor and TrackNode classes
// All require a ROOT file generated with tracking information ON.
//
// see https://github.com/snoplus/rat/blob/master/example/macros/tracking/tracking.mac
//     https://snopl.us/docs/rat/user_manual/html/working_with_tracks_in_root.html
//
// Particle tracks are divided into a series of backward looking 'Nodes' containing
// the step information. The TrackCursor class controls the navigation through these
// nodes and between tracks.
//
// Secondary particles appear new 'child' tracks at the node that created them.
//
// Primary particles are the children of an inital track with a single node
// representing the event
//
// J dunger - 2016-03-14 Hide TrackTest example from CINT
// J dunger - 2014-10-14 TrackTest_Energy namespace fix
// J dunger - 2014-08-20 <jack.dunger@ox.ac.uk> : New file
///////////////////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include <RAT/DU/DSReader.hh>
#include <RAT/TrackNav.hh>
#include <RAT/TrackCursor.hh>
#include <RAT/TrackNode.hh>
#include <RAT/DB.hh>

#include <string>
#endif
#include <TH1D.h>

//A function to plot the momentum of the primary particles
TH1D* PrimaryMomentum(const std::string& filename){

  TH1D* hist = new TH1D("PrimaryMomentum","",10,0,5);
    // 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 reader(filename);

    //Event Loop
  for(size_t iEv =0; iEv<reader.GetEntryCount(); iEv++){
    const RAT::DS::Entry& ds = reader.GetEntry(iEv);
    RAT::TrackNav nav(&ds);
      RAT::TrackCursor cursor =nav.Cursor(false);  //toggle human readable cursor

      for(size_t iCh = 0; iCh<(size_t)cursor.ChildCount(); iCh++){
       RAT::TrackNode* node = cursor.Child(iCh);
       hist->Fill(node->GetMomentum().Mag());
     }
   }

   return hist;
 }

//A function to plot the final displacement of the 1st Primary Particle
 TH1D* Primary1Displacement(const std::string& filename){

  TH1D* hist = new TH1D("Primary1Displacement","",10,0,100);
  // 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 reader(filename);

    //Event Loop
  for(size_t iEv =0; iEv<reader.GetEntryCount(); iEv++){
    const RAT::DS::Entry& ds = reader.GetEntry(iEv);
    RAT::TrackNav nav(&ds);
    RAT::TrackCursor cursor =nav.Cursor(false);

    if(cursor.ChildCount()){
     cursor.GoChild(0);
     RAT::TrackNode* node = cursor.Here();
     TVector3 init_pos  = node->GetPosition();

     node = cursor.TrackEnd();
     TVector3 fin_pos   = node->GetPosition();

     double disp = (fin_pos - init_pos).Mag();
     hist->Fill(disp);
   }
 }

 return hist;
}

// A function to cout the name of all processes that create secondary particles
void FindSecondaries(const std::string& filename){

  RAT::DU::DSReader reader(filename);

  for(size_t iEv =0; iEv<reader.GetEntryCount(); iEv++){
    const RAT::DS::Entry& ds = reader.GetEntry(iEv);
    RAT::TrackNav nav(&ds);
    RAT::TrackCursor cursor = nav.Cursor(false);

    for(size_t iCh = 0; iCh<(size_t)cursor.ChildCount(); iCh++){
     cursor.GoChild(iCh);
     while(!cursor.IsTrackEnd()){
       RAT::TrackNode* node = cursor.Here();
       if(cursor.ChildCount())
         std::cout << node->GetProcess() << std::endl;
       cursor.GoNext();
     }
     cursor.GoTrackStart();
     cursor.GoParent();

      } //Primary Particle Tracks

    } //event

    return;
  }

/*
In the following example CINT fails to recognise the inheritance of TrackTestEnergy and tries to call the pure virtual method
TrackTest::operator()(TrackNode *n).
Works if compiled with gcc or root ACLiC.
*/

#ifndef __CINT__
// A class and function to plot the inital radial position of all electrons with inital energy >1MeV
  namespace RAT{
    class TrackTest_Energy : public TrackTest {
      double fEnergy;
    public:
      TrackTest_Energy(const double& Energy) : fEnergy(Energy) {}
      virtual bool operator() (TrackNode *n) {return n->GetKineticEnergy() > fEnergy;}

    };
  }
  TH1D* PlotInitalPosECut(const std::string& filename){
    TH1D* hist = new TH1D("PrimaryMomentum","",10,0,5);
    RAT::DU::DSReader reader(filename);
    RAT::TrackTest *t = new RAT::TrackTest_Energy(1);

    //Event Loop
    for(size_t iEv =0; iEv<reader.GetEntryCount(); iEv++){
      const RAT::DS::Entry& ds = reader.GetEntry(iEv);
      RAT::TrackNav nav(&ds);
      RAT::TrackCursor cursor = nav.Cursor(false);  //toggle human readable Cursor
      RAT::TrackNode   *n     = cursor.Here();
      while (n != NULL){
       hist->Fill(n->GetPosition().Mag());
       n = cursor.FindNextTrack(t);
     }
   }
   return hist;
 }
#endif