#ifndef ANTDSTREADERINC
#define ANTDSTREADERINC

#include "EventFile.hh"
#include "io_util.hh"

#define private public
#include "AntDST.h"
#undef private


inline void generate_antdst_project()
{
  string projectdir = "antdst";
  /* you may have to just copy-paste these two lines into root, to boodstrap the process */
  auto f = TFile::Open("~/MC_032189_anumu_a_CC_reco.root");
  f->MakeProject( projectdir.c_str() ,"*","recreate+"); // make library, but do not load
}



// def beta( trk ) :

// 	if trk.error_matrix.size() < 25 : return -1
// 	sigma_theta2 = trk.error_matrix[ 5*3+3]
// 	sigma_phi2   = trk.error_matrix[ 5*4+4]
// 	if sigma_theta2 > 0 and sigma_phi2 > 0 :
// 		return (sigma_theta2 + sigma_phi2)**0.5
// 	else:
// 	 print ("negative sigma2")

// 	return -1


inline double compute_beta( const vector<double>& error_matrix )
{
  if (error_matrix.size() < 25 ) return -1;

  const double sigma_theta2 = error_matrix[ 5*3+3 ];
  const double sigma_phi2   = error_matrix[ 5*4+4 ];
  if ( sigma_theta2 > 0 && sigma_phi2 > 0 ) 
    return sqrt( sigma_theta2 + sigma_phi2 );
  
  //print ("negative sigma2 in beta computation.");
  return -2;
}


inline void read( Trk& t, const Particle& p , bool read_energy = true, int level =3 )
{
  t.type = p.fParticleType;
  //t.E    = p.fEnergy; // that would be too easy.

  if ( read_energy )
  {
    if ( p.fAntEnergy.size() != 1 ) {
      print ("warning: strange size of fAntEnergy for mc particle");
    }
    if ( p.fAntEnergy.size() > 0 ) t.E = p.fAntEnergy[0].fEnergy;
  }

  t.pos.set( p.fPosition.X(), p.fPosition.Y(), p.fPosition.Z() );
  t.dir.set( p.fAxis.X(), p.fAxis.Y(), p.fAxis.Z() );

  // not sure if the following are good for anything, but hey...
  if ( level > 1 ) {
    t.setusr("RightAscension",     p.fRightAscension);
    t.setusr("Declination",        p.fDeclination);
  }

 

  if ( level > 2 ) {
    t.setusr("GalacticLong",       p.fGalacticLong);
    t.setusr("GalacticLat",        p.fGalacticLat);
    t.setusr("RealGalacticLat",    p.fRealGalacticLat);
    t.setusr("RealGalacticLong",   p.fRealGalacticLong);
    t.setusr("RealDeclination",    p.fRealDeclination);
    t.setusr("RealRightAscension", p.fRealRightAscension);
    t.setusr("Zenith",             p.fZenith);
    t.setusr("Azimuth",            p.fAzimuth);
  }

}


inline void read_energy( Trk& t, const RecParticle& rp , int level = 3  )
{
  for(auto e: rp.fAntEnergy ) {

    //print ( "energy", e.fErecStrategy , t.rec_type );

    if ( e.fErecStrategy == 1 ) /// dEdX (CEA: A. Romeyer, G. Vannoni, K. Payet)
    {
      t.E   = e.fEnergy; 
      t.len = e.fParameters["TrackLength"];
      
      if (level < 2 ) {
        t.setusr("dEdX_dEdX", e.fParameters["dEdX"]  );
      }

      else {
        t.setusr("EdEdX",  e.fEnergy ); 
        foreach_map( key,val, e.fParameters ) {  
          t.setusr("dEdX_"+key, val );
        }
      }
    }

    if ( e.fErecStrategy == 2 && level > 2 )  // NeuralNetwork (ECAP: J. Schnabel)
    {
      t.setusr("EANN",  e.fEnergy );
    }

    if ( e.fErecStrategy == eTantraShowerEnergy) {

      if ( t.rec_type != eShowerTantraFit ) {
        print ("? found tantra E reco for other type of track", t.rec_type );
      }
      else {
        t.E = e.fEnergy;
      }
    }

  }
}


inline void read( Trk& t, const RecParticle& rp , int level = 3 ) 
{
  read( t, (Particle&) rp , false , level );

  t.error_matrix.clear();
  foreach( r, rp.fErrorMatrix ) foreach( x, r ) t.error_matrix.push_back(x);

  if ( level < 2 ) { // just store beta, not full error matrix.
    const double beta = compute_beta( t.error_matrix );
    t.error_matrix.resize(1);
    t.error_matrix[0] = beta;
  }

  // energy -- it's complicated.
  read_energy( t, rp , level );

  t.fitinf.push_back(rp.fNUsedLines   );
  t.fitinf.push_back(rp.fNUsedHits    );
  t.fitinf.push_back(rp.fTotalUsedAmp );
}


inline void read(Trk& t, const RecEvent& rev , int level = 3 )
{
  t.rec_stages.resize(1);
  t.rec_type       = rev.fStrategy;
  t.rec_stages[0]  = rev.fRecStep; 
  t.lik            = rev.fRecQuality; 
 
  if ( rev.fRecParticles.size() != 1 ) {

    if ( rev.fRecParticles.size() != 0 ) // seems to be normal?
    {
      dump(rev);
      print ("antdst warning: strange number of recparticles in recevent", rev.fRecParticles.size() );
      if (rev.fRecParticles.size() == 0 ) return; // invalid read.
    }
    t = Trk();
    return;
  } 

  t.fitinf.clear();
  t.fitinf.push_back( rev.fNTotalHits );     //0 
  t.fitinf.push_back( rev.fNTotalLines );   
  t.fitinf.push_back( rev.fTotalAmplitude );   
  t.fitinf.push_back( rev.fNTriggeredHits );   
  t.fitinf.push_back( rev.fNTriggeredLines );   
  t.fitinf.push_back( rev.fTriggerAmplitude );    //6
  //nb: 3 more items added to fitinf in followring read.
  read( t, rev.fRecParticles[0], level );

  foreach( p , rev.fAdditionalParameters ) {
    t.fitinf.push_back(p);
  }

}

inline void read(Hit& h, const AntDSTHit& a )
{
  h.id    = a.fHitId;
  h.dom_id= a.fOMId;
  h.pos.set( a.fOMPosition.x(), a.fOMPosition.y(), a.fOMPosition.z() );
  h.a     = a.fCharge;
  h.t     = a.fTime;
  h.type  = a.fHitType;
}

inline Trk& getat( vector<Trk>& v , unsigned i ) {


  if ( i+1 > v.size() ) v.resize( i+1 );
  return v[i];
}


inline void read( Evt& e, const AntDST& a, int level = 3 )
{
  e.clear(); 
  e.w.resize(4);
  e.w[0] = a.fMCEvent.fIDPrim;      // abusing this
  e.w[1] = a.fMCEvent.fWeight_W2;
  e.w[2] = a.fMCEvent.fWeight_W3;
  e.w[3] = a.fMCEvent.fNGenEvents;  // maybe a good idea.

  foreach( p , a.fMCEvent.fParticles ) {

    int typ =  abs( p.fParticleType );
    bool nu = typ == 12 || typ == 14 || typ == 16;

    if ( nu || level > 2 ) {
      Trk t;
      read( t, p , level, level );
      e.mc_trks.push_back( t );
    }
  }
  
  // recevents

  e.trks.clear();

  for (unsigned i=0; i< a.fRecEvents.size(); i++) 
  {
    if ( level < 2 ) 
    {
      int typ = a.fRecEvents[i].fStrategy;
      if ( typ == eAAFit )          read( getat( e.trks, 0 ) , a.fRecEvents[i] ,level );
      if ( typ == eShowerTantraFit) read( getat( e.trks, 1 ) , a.fRecEvents[i] ,level );
      if ( typ == eShowerDusjFit)   
      { 
          e.setusr("dusjQ", double( a.fRecEvents[i].fRecQuality) );
      }
    }
    else 
    {
      read( getat( e.trks, i) , a.fRecEvents[i] , level );
    }
  }

  e.mc_run_id       = a.fHasMC? 1:0;
  e.run_id          = a.fRunId;
  e.frame_index     = a.fFrameIndex;
  e.trigger_counter = a.fTriggerCounter;
  e.index           = a.fFrameTarget;
  e.id              = a.fEventId; // may go wrong int->longlong
  //                = a.fJD;
  e.t               = a.fEventTime;
  //                  a.fRunStart;
  //usr             = a.fUserData;

  // fUserCommentsMap;

  if (level > 2)
  {
    for (auto p : a.fUserCommentsMap)
    {
      e.comment += p.first + ":" + p.second + ", ";
    }
  }

  // fTriggerType
  for(auto& t : a.fTriggerType ) {
    if ( t > 63 ) print("too large trigger-type in antdst");
    e.trigger_mask &= (1<<unsigned(t));
  }


  if (level>2)
    {
      // fHits
      Hit hit;
      e.hits.clear();
      foreach_map( k,v, a.fHits) {
        (void) k;
        read( hit, v );
        e.hits.push_back(hit);
      }
    }
}

struct AntDstEventFile : public RootEventFile
{

  /*! you can set the level to 1 or 2 to read less information than the default */
  static int level( int l = -1  ) 
  { 
    static int lvl = 3;
    if ( l == -1 ) return lvl;
    return lvl = l;
  }

  AntDST* aevt;
  
  string treename()   { return "AntData"; }
  string branchname() { return "Event.";}
  
  void setbranchaddress( TBranch* b ) { b->SetAddress( &aevt);}
  void readevt(Evt& evt )             { read( evt, *aevt, level() );    }
  bool read_header(Head& h);           // in cc

  AntDstEventFile(): aevt(0) {};

  virtual string info() { return "AntDstEventFile"; }
  virtual AntDstEventFile* clone() { return new AntDstEventFile(*this); }
};


/*! helper function to access the AntDST inside the reader -- in case you need
    access to the original antdst object */
inline AntDST* access_antdst_event( EventFile* f )
{

  AntDstEventFile* a = dynamic_cast<AntDstEventFile*>( f->impl );
  if ( !a ) {
    print ("failed to get antdst from Eventfile.");
    return nullptr;
  }

  return a->aevt;
}


inline void test_antdst()
{
  /* Load the library that is created by generate_project() */

  gInterpreter->Load("antdst/antdst.so");

  
  /* Register the antdsteventfile with the EventFile */

  AntDstEventFile::level(2);

  auto* p =  new AntDstEventFile();  
  EventFile::_registered_filetypes().push_back( p );

  auto f = EventFile( "/pbs/home/h/heijboer//MC_032189_anumu_a_CC_reco.root" ) ;

  // from now on, use f as any aanet input file
  
  for( auto evt : f)
    {
      print ( evt.frame_index );
    }

}

#endif