// StateManager.h

   #ifndef __ScoreKeeper_StateManager_h__
#define __ScoreKeeper_StateManager_h__

#include "Track.h"
#include "StackManager.h"
#include "Callbacks.h"

#ifdef __ENABLE_STARMC__
  #include "TGiant3.h"
  #include "St_geant_Maker/St_geant_Maker.h"
  #include   "StarMCCallbacks.h"
#endif

//#define __ENABLE_G4__
#ifdef  __ENABLE_G4__
  #include "G4Event.hh"
  #include "G4Track.hh"
  #include "G4EventManager.hh"
  #include   "G4Callbacks.h"
#endif

#include <iostream>

// Use compensated sum for accumulated quantities.  The second element of accumulated
// arrays will be the compensation term
//   https://en.wikipedia.org/wiki/Kahan_summation_algorithm

namespace ScoreKeeper {

  struct State_t {
    // event
    int    idevent;
    // track identity
    int    pdgid;
    // track enumeration, to be defined by stacking mechanism
    int    idtrack;
    int    idparent; 
    // state flags
    int    flags; 
    // current state of track on geometry
    int    nstep;
    double step;
    double astep[2]; // accumulated step size
    double x, y, z, px, py, pz;
    double tof; // time of flight (from prod vtx to here)
    double etotal; // total energy
    double atof[2]; // time of flight since collion (accumlated over parent tracks)
    // physical response of detector to passage of particle
    double destep;
    double adestep[2]; // accumulated destep

  };

  template<typename MC>
  class StateManager { 
  public:

    typedef MC MonteCarlo_t;
    
    // Update the current tracking state and return a reference
    // to the state data
    State_t& Update() {

      // Previous state
      _prev = _state;

      // Obtain state flags
      _state.flags    = getStateFlags ( _mc );

      _state.idevent  = getIdEvent    ( _mc );
      // _state.idtrack  = getIdTrack    ( _mc );
      _state.pdgid    = getTrackPdgId ( _mc );
      // _state.idparent = getIdParent   ( _mc );
      _state.nstep    = getStepNumber ( _mc );
      _state.step     = getStepSize   ( _mc );

      _state.x         = getx( _mc );
      _state.y         = gety( _mc );
      _state.z         = getz( _mc );
      _state.tof       = gett( _mc );
      _state.px        = getpx( _mc );
      _state.py        = getpy( _mc );
      _state.pz        = getpz( _mc );
      _state.etotal    = getTotalEnergy( _mc );
      _state.destep    = getEnergyDeposit( _mc );

      kahanSum( _state.astep,   _state.step );
      kahanSum( _state.atof,    _state.tof );
      kahanSum( _state.adestep, _state.destep );

      return _state;

    }

    // Signal a new event, clear the state and increment event id.
    // (Event ID may be overwritten by underlying MC)
    State_t& NewEvent() {
      State_t old = _state;
      State_t s = {};
      _state = s;
      _state.idevent = old.idevent + 1;
      return _state;
    }

    // Signal a new track, clear all and increment track id
    // (Track ID may be overwritten by underlying MC)
    State_t& NewTrack() {
      State_t old = _state;
      State_t s = {};
      _state = s;
      _state.idtrack = old.idtrack + 1;
      return _state;
    };

    State_t& Current(){ return _state; }
    State_t& Previous(){ return _prev; }
    
  private:
  protected:
    
    State_t      _state; // current state
    State_t      _prev;  // previous state

    MonteCarlo_t _mc;    // MC type
    
    // Performs kahan compensated summation
    void kahanSum( double* x, const double in ) {
      double& sum = x[0];
      double& c   = x[1];
      double  y   = in - c;  // On first call, c = 0
      double  t   = sum + y; // Alas, sum is big, y small, so low-order digits of y are lost.
      c = (t - sum) - y;     // (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
      sum = t;
    };
    
    

  };

  std::ostream& operator << (std::ostream &out, const State_t& s )
  {
    out << "-- mc state --------------------------------------------------------" << std::endl
	<< "event = " << s.idevent << std::endl
	<< "track = " << s.idtrack << std::endl
	<< "pdgid = " << s.pdgid   << std::endl
	<< "flags = " << s.flags   << std::endl
	<< "astep = " << s.astep[0] << std::endl
	<< "atof  = " << s.atof[0] << std::endl
	<< "adest = " << s.adestep[0] << std::endl;
    return out;
  }

};

#endif



// G4Callbacks.h



#ifndef __ScoreKeeper_Geant4Callbacks_h__
#define __ScoreKeeper_Geant4Callbacks_h__

#include "Callbacks.h"
#include <iostream>

namespace ScoreKeeper {

#ifdef __ENABLEG4__
  template<>   int getStateFlags<Geant4>( Geant4 ){ return 0; }
  template<> 
  int getIdEvent   <Geant4>( Geant4 ){ 
    int id = 0;
    auto* mgr = G4EventManager::GetEventManager();
    auto* evt = mgr->GetConstCurrentEvent();
    id = evt->GetEventID();
    return id; 
  } 
#endif 

};

#endif