#ifndef SHOWERECUTILS
#define SHOWERECUTILS

#include "Trk.hh"
#include "Table.hh"
#include "stringutil.hh"
using stringutil::split;

// --------------------------------------------------
// functions for converting between double[] and Trk.
// in aid of fitting.
//---------------------------------------------------

inline vector<string> parnames()
{
  static vector<string> r;
  if (r.size()==0) r = split(string("x y z theta phi t logE"));
  return r;
}

inline map<uint, std::pair<double, double> > parlimits()
{
  map<uint, std::pair<double, double> > l;
  l[0] = std::make_pair( -1000, 1000 );
  l[1] = std::make_pair( -1000, 1000 );
  l[2] = std::make_pair( -1000, 1000 );
  l[3] = std::make_pair( -0.1, 3.2 );
  l[4] = std::make_pair( -0.1, 6.3 );
  l[5] = std::make_pair( -3000, 3000 );
  l[6] = std::make_pair(  0.0, 10 );
  return l;
}

inline void to_pars( Trk& t, double* r )
{
  r[0] = t.pos.x; r[1] = t.pos.y; r[2] = t.pos.z;
  r[3] = acos( t.dir.z );
  r[4] = atan2( t.dir.y, t.dir.x );
  r[5] = t.t;
  r[6] = log10(t.E);
}

inline void to_trk( const double* r, Trk& t )
{
  t.pos.set( r[0], r[1], r[2] );
  t.dir.set_angles( r[3],r[4] );
  t.t = r[5];
  try{
    t.E = pow(10,r[6]);
  } catch( string except ){
    std::cout << "No energy: " << except << std::endl;
    t.E = -1;
  }
}

struct JppUtil{

public:

  JppUtil(){};
  
  vector<string> parnames()
  {
    static vector<string> r;
    if (r.size()==0) r = split(string("logE x y z theta phi t"));
    return r;
  }

  map<uint, std::pair<double, double> > parlimits()
  {
    map<uint, std::pair<double, double> > l;
    l[0] = std::make_pair(  0.0, 14 );
    l[1] = std::make_pair( -1000, 1000 );
    l[2] = std::make_pair( -1000, 1000 );
    l[3] = std::make_pair( -1000, 1000 );
    l[4] = std::make_pair( -0.1, 3.2 );
    l[5] = std::make_pair( -0.1, 6.3 );
    l[6] = std::make_pair( -3000, 3000 );
    return l;
  }

  void to_pars( Trk& t, double* r )
  {
    r[0] = log10(t.E);
    r[1] = t.pos.x; r[2] = t.pos.y; r[3] = t.pos.z;
    r[4] = acos( t.dir.z );
    if   ( r[4] <    0.0 ) r[4] = -r[4];
    while( r[4] > 2*M_PI ) r[4] -= 2 * M_PI;
    if   ( r[4] >   M_PI ) r[4] = 2*M_PI - r[4];
    r[5] = atan2( t.dir.y, t.dir.x );
    while( r[5] > 2*M_PI ) r[5] -= 2 * M_PI;
    while( r[5] <    0.0 ) r[5] += 2 * M_PI;
    r[6] = t.t;
  }

  void to_trk( const double* r, Trk& t )
  {
    t.pos.set( r[1], r[2], r[3] );
    t.dir.set_angles( r[4],r[5] );
    t.t = r[6];
    try{
      t.E = pow(10,r[0]);
    } catch( string except ){
      std::cout << "No energy: " << except << std::endl;
      t.E = -1;
    }
  }
};

inline void print_error_matrix( Trk& t , ostream& out = cout )
{
  
  vector<string> s = parnames();
  s.insert(s.begin(), " x ");
  Table T( s );
  int N = parnames().size();

  vector<double> scaler = pack( 1.0,1,1, 180/M_PI, 180/M_PI, 1);

  enumerate( j, pn, parnames() )
    {
      T << pn;
      for( int i=0; i<N; i++)
	{
	  double v = t.error_matrix[i+N*j];
	  v*=scaler[i];
	  v*=scaler[j];

	  T << v;
	}
    }
  
  T.print_ascii( out );
  
}

inline vector< Vec > icosahedron()
{
  static vector<Vec> r;
  
  if ( r.size() == 0 ) 
    {
      // see wikipedia : icosahedron
      const double a = ( sqrt(5.0) - 1 ) / 2.0;
      foreach( s, pack(-1,1) ) foreach( f, pack(-a,a) )
        {
          r.push_back( Vec( 0, s, f ) );
          r.push_back( Vec( s, f, 0 ) );
          r.push_back( Vec( f, 0, s ) ); 
        }
    }
  return r;
}




// --- error handler for root to suppress Minuit msgs ----

#include "TSystem.h"

inline void fit_errorhandler (int level, Bool_t abort, const char *location,
		       const char *msg)
{
  if (level < 1001 ) return;
  cout << " There was a root error. level=" << level << ", msg: " << msg << endl;

  if (level < 2001 ) return;
  cout << " stack trace will follow and then we'll exit " << endl;
  cout << endl << endl;
  gSystem->StackTrace();
  exit(1);
  return;
}



#endif