#ifndef AATRACKINGACTIONHH
#define AATRACKINGACTIONHH

#include "Sim.hh"

#ifndef __CLING__
#include "G4UserSteppingAction.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "G4UserTrackingAction.hh"
#include "G4UserStackingAction.hh"
#include "G4Track.hh"
#endif

#include "CLHEP/Units/PhysicalConstants.h"
using namespace CLHEP;





inline Vec toVec( const G4ThreeVector& v )
{
  using CLHEP::m;
  return Vec( v.x()/m, v.y()/m, v.z()/m );
}

struct AAStackingAction : public G4UserStackingAction
{


  G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track *t)
  {
    const double beta = t->GetVelocity() / CLHEP::c_light;

    int typ = t->GetParticleDefinition()->GetPDGEncoding();

    if (typ == 11 || typ == -11) // electron/ position
    {
      if (beta < 0.7)
      {
        //cout << " kill sub-cherenkov electron " << endl;
        return fKill;
      }
    }

    if (typ == 22)
    {
      if (t->GetTotalEnergy() / GeV < 0.00051099891 * 2) // can't pair produce
      {
        return fKill;
      }
    }

    return fUrgent;
  }
};

struct AATrackingAction : public G4UserTrackingAction
{
  vector<Trk>* output_tracks;
  vector<Trk>* simulated_tracks;

  double energy_threshold; // for tracks to be stored 

  Trk trk;  

  AATrackingAction() {}

  
  virtual void PreUserTrackingAction(const G4Track* t)  	
  {
    TIMESCOPE

    trk.id  = t -> GetTrackID ();
    trk.t   = t -> GetGlobalTime ();    
    trk.E   = t -> GetTotalEnergy () / GeV;
    trk.pos.set(  t -> GetPosition().x() / m ,  t -> GetPosition().y() / m,  t -> GetPosition().z() / m );

    // if ( g_sim -> verbose_tracks)
    //   {	
	  //     get_track_pool().add_track( trk.id, trk );
    //   }
  }

  virtual void PostUserTrackingAction(const G4Track *t)
  {

    // check that we got the track we expected

    if (trk.id != t->GetTrackID())
      fatal("wrong track id ");

    if (g_sim->record_tracks)
    {
      trk.mother_id = t->GetParentID();
      if (trk.mother_id == 0)
        trk.mother_id = -1;
      trk.type = t->GetParticleDefinition()->GetPDGEncoding();
      trk.status = t->GetTrackStatus();

      Vec endpos = toVec(t->GetPosition());
      trk.dir = (endpos - trk.pos).normalize();
      trk.len = t->GetTrackLength() / m;
      g_sim->tracks.push_back(trk);
    }

    // if (g_sim->verbose_tracks)
    // {
    //   get_track_pool().print(t->GetTrackID());
    //   print("after : ", trk.__str__());
    // }
    return;

    //---------------------------------------------------------------
    // find out if the track is one of the primary (input) particles
    //---------------------------------------------------------------

    bool is_input_track = false;

    const G4Event *evt = G4RunManager::GetRunManager()->GetCurrentEvent();

    for (int i = 0; i < evt->GetNumberOfPrimaryVertex(); i++)
    {
      if (t->GetTrackID() == evt->GetPrimaryVertex(i)->GetPrimary(0)->GetTrackID())
      {
        is_input_track = true;
      }
    }

    // skip non-input track that are below threshold
    const double mass = t->GetParticleDefinition()->GetPDGMass() / GeV;

    // do a cut on *kinetic* energy -- otherwise we end up with a bunch of protons and neutrons
    if (!is_input_track && ((trk.E - mass) < energy_threshold))
      return;

    Vec p( t->GetPosition().x() / m, t->GetPosition().y() / m, t->GetPosition().z() / m);

    // trk.rec_type = is_input_track;
    trk.len = (p - trk.pos).len();
    trk.dir = (p - trk.pos) / trk.len;
    trk.type = t->GetParticleDefinition()->GetPDGEncoding();
    trk.id = 1;

    static bool dbg = true;
    if (dbg)
    {
      cout << " post tracking particle " << t->GetParticleDefinition()->GetParticleName()
           << " pdg_code=" << t->GetParticleDefinition()->GetPDGEncoding()
           << endl;
    }
  }
};
#endif