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

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

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

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

#include "JGeometry2D/JVector2D.hh"
#include "JGeometry2D/JCircle2D.hh"

#include "JGeometry3D/JCylinder3D.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JGeane.hh"

#include "JSirene/JVisibleEnergyToolkit.hh"


/**
 * \file
 * 
 * Program to test visible energy toolkit.
 * \author bjung
 */
int main(int argc, char** argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JAAnetTypes_t> inputFile;
  double                              epsilon;
  int                                 debug;

  try {

    JParser<> zap("Program to test visible energy toolkit");

    zap['f'] = make_field(inputFile);
    zap['e'] = make_field(epsilon) = 1e-4;
    zap['d'] = make_field(debug)   = 2;

    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  unsigned int n = 0;
  unsigned int N = 0;
  
  const Head& header = getHeader(inputFile);

  const JCylinder3D can = get<JCylinder3D>(header);
    
  STATUS(LEFT(15) << "event number" << RIGHT(15) << "E0 [GeV]" << RIGHT(15) << "Evis [GeV]" << endl);
  
  for (JMultipleFileScanner<Evt>& in = inputFile; in.hasNext(); ) {

    const Evt& evt = *(in.next());

    const double E0   = getE0(evt);
    const double Evis = getVisibleEnergy(evt, can);

    STATUS(RIGHT(15) << in.getCounter() << FIXED(15,3) << E0 << FIXED(15,3) << Evis << '\r'); DEBUG(endl);

    if (Evis > 0.0 && Evis > E0) {
      
      STATUS(endl); WARNING("Evis > E0 for event " << in.getCounter() << endl);

      DEBUG(RIGHT(15) << "track PDG ID" << RIGHT(37) << "position" << RIGHT(35) << "direction" << RIGHT(15) << "E [GeV]" << RIGHT(15) << "Evis [GeV]" << endl);      

      for (vector<Trk>::const_iterator i = evt.mc_trks.cbegin(); i != evt.mc_trks.cend(); ++i) {
	DEBUG(RIGHT(15) << i->type << RIGHT(20) << i->pos << RIGHT(20) << i->dir << RIGHT(15) << i->E << RIGHT(15) << getVisibleEnergy(*i, can) << endl);
      }
      
      ++N;
    }

    ++n;
  }

  const double fraction = N / ((double) n);
  
  ASSERT(fraction < epsilon);
}