#include <iostream>
#include <string>
#include <vector>

#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"

#include "JMath/JRandom.hh"

#include "JAAnet/JAAnetTestkit.hh"
#include "JAAnet/JParticleTypes.hh"

#include "JSirene/JEventShapeVariables.hh"

#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"

#include "Jeep/JeepToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

#include "TRandom3.h"


/**
 * \file
 *
 * Auxiliary program to test event shape variables.
 * \author bjjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<Evt>     inputFile;
  JLimit_t&                     numberOfEvents = inputFile.getLimit();
  double                        epsilon;
  UInt_t                        seed;
  int                           debug;

  try { 

    JParser<> zap("Auxiliary program to test event shape variables.");

    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['e'] = make_field(epsilon)             = 1e-6;
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 3;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
  
  
  gRandom->SetSeed(seed);

  if (debug >= JEEP::debug_t) {
    setLongprint(cout);
  }

  use_root = true;


  // Back-to-back event

  NOTICE("Checking sphericity tensor for single back-to-back event" << endl);

  Evt event0;

  Trk track0    = getRandom<Trk>();
  track0.status = TRK_ST_FINALSTATE;
  track0.type   = TRACK_TYPE_ELECTRON;
  
  Trk track1    = track0;
  track1.type   = TRACK_TYPE_ANTIELECTRON;  
  track1.dir   *= -1;

  event0.mc_trks.push_back(track0);
  event0.mc_trks.push_back(track1);

  const JSphericityTensor S(event0);
  
  const JMatrix3S::eigen_values& lambda = S.getEigenValues();
  
  const double aplanarity = S.getAplanarity();
  const double sphericity = S.getSphericity();
  
  const double thrust = getThrust(event0);

  DEBUG("lambda = " <<
	SCIENTIFIC(10,5) << lambda[0] << endl <<
	SCIENTIFIC(21,5) << lambda[1] << endl <<
	SCIENTIFIC(21,5) << lambda[2] << endl);

  NOTICE("aplanarity = " << SCIENTIFIC(15,5) << aplanarity << endl); 
  NOTICE("sphericity = " << SCIENTIFIC(15,5) << sphericity << endl);
  NOTICE("thrust     = " << SCIENTIFIC(15,5) << thrust << endl);
  
  ASSERT(lambda[0] + lambda[1] + lambda[2] < 1.0 + epsilon);

  ASSERT(aplanarity < epsilon);
  ASSERT(sphericity < epsilon);
  ASSERT(thrust     < 1.0 + epsilon);


  // Check hit inertia tensor
  
  NOTICE("Checking hit inertia tensor for random event" << endl);
  
  Evt event1;

  for (int i = 0; i < 20; ++i) {
    const Hit hit = getRandom<Hit>();
    event1.hits.push_back(hit);
  }

  const JHitInertiaTensor hitInertiaTensor(event1, JPosition3D());

  const double inertiaRatio = hitInertiaTensor.getEigenvalueRatio();

  NOTICE("inertiaRatio = " << SCIENTIFIC(15,5) << inertiaRatio << endl);

  ASSERT(fabs(inertiaRatio) > -epsilon && fabs(inertiaRatio) < 1.0 + epsilon);


  // Events from file

  NOTICE("Checking event shape variables for test file.");

  NOTICE(LEFT(10) << "Event ID" << RIGHT(15) << "aplanarity" << RIGHT(15) << "sphericity" << RIGHT(15) << "thrust" << RIGHT(15) << "H0" << RIGHT(15) << "H1" << endl);

  while (inputFile.hasNext()) {

    const Evt* event = inputFile.next();

    const double E0 = getE0(*event);
    const double E1 = getE1(*event);

    if (fabs((E0 - E1) / E0) > epsilon) {
      WARNING("Event " << inputFile.getCounter() << " does not uphold energy conservation! skip." << endl);
      numberOfEvents.add(1);
      continue; // Skip events for which energy conservation is not upheld
    }

    const JSphericityTensor S(*event, 2.0);

    const JMatrix3S::eigen_values& lambda = S.getEigenValues();

    const double aplanarity = S.getAplanarity();
    const double sphericity = S.getSphericity();
    
    const double thrust = getThrust(*event);

    const JFoxWolframMoments H(*event, getMaximumContainmentVolume(), 4);

    DEBUG("lambda sphericity = " <<
	  SCIENTIFIC(10,5) << lambda[0] << endl <<
	  SCIENTIFIC(21,5) << lambda[1] << endl <<
	  SCIENTIFIC(21,5) << lambda[2] << endl);

    NOTICE(LEFT(10)         << inputFile.getCounter() <<
	   SCIENTIFIC(15,5) << aplanarity             <<
	   SCIENTIFIC(15,5) << sphericity             <<
	   SCIENTIFIC(15,5) << thrust                 <<
	   SCIENTIFIC(15,5) << H.at(0)                <<
	   SCIENTIFIC(15,5) << H.at(1)                <<	   
	   '\r'); DEBUG(endl);
    
    ASSERT(lambda[0] + lambda[1] + lambda[2] < 1.0 + epsilon);

    ASSERT(aplanarity > -epsilon       && aplanarity < 0.5 + epsilon);
    ASSERT(sphericity > -epsilon       && sphericity < 1.0 + epsilon);    
    ASSERT(thrust     >  0.5 - epsilon && thrust     < 1.0 + epsilon);
    
    ASSERT(H.size() == 5);
    ASSERT(fabs(H[0]) > -epsilon);
    ASSERT(fabs(H[1]) > -epsilon);
    ASSERT(fabs(H[2]) > -epsilon);
    ASSERT(fabs(H[3]) > -epsilon);
    ASSERT(fabs(H[4]) > -epsilon);
  }

  return 0;
}