#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

#include "km3net-dataformat/definitions/weightlist.hh"

#include "km3net-dataformat/tools/reconstruction.hh"

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

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

#include "JDetector/JDetector.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JDetectorToolkit.hh"

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

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

#include "JTools/JRange.hh"

namespace {

  struct JEvtWeightOption {

    static const char* const unweighted() { return "unweighted";   }
    static const char* const standard()   { return "standardFlux"; }
    static const char* const rescaled()   { return "rescaledFlux"; }
  };
}


/**
 * \file
 *
 * Example program to analyse track fit results from Evt formatted data.
 * \author bofearraigh, bjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  JMultipleFileScanner_t inputFiles;
  
  JLimit_t               numberOfEvents;
  
  string                 detectorFile;
  string                 outputFile;

  JRange<double>         energyRange;
  JRange<double>         coszenithRange;
  JRange<unsigned int>   triggeredHitsRange;
  
  string                 weightOption;
  int                    debug;
  
  try {

    JProperties cuts;

    cuts.insert(gmake_property(energyRange));
    cuts.insert(gmake_property(coszenithRange));    
    cuts.insert(gmake_property(triggeredHitsRange));

    JParser<> zap("Example program to analyse track fit results from Evt formatted data.");

    zap['f'] = make_field(inputFiles);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile)        = "";
    zap['o'] = make_field(outputFile)          = "histogram.root";
    zap['@'] = make_field(cuts)                = JPARSER::initialised();
    zap['W'] = make_field(weightOption)        = 
      JEvtWeightOption::unweighted(),
      JEvtWeightOption::standard(),
      JEvtWeightOption::rescaled();
    zap['d'] = make_field(debug)               = 2;

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

  if (weightOption   != JEvtWeightOption::unweighted() &&
      numberOfEvents != JLimit::max()) {
    FATAL("Cannot apply weighting to limited number of events.");
  }
  
    
  JDetector detector;

  if (detectorFile != "") {
    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
  }

  
  const JPMTRouter router(detector);  

  // Histogram bins

  const Int_t    N_Nhits       =               400;
  const Double_t Nhits_min     =         0.0 - 0.5;
  const Double_t Nhits_max     =      1000.0 - 0.5;
  
  const Int_t    N_Noverlays   =                40;
  const Double_t Noverlays_min =         0.0 - 0.5;
  const Double_t Noverlays_max =      1000.0 - 0.5;

  const Int_t    N_radius      =               201;
  const Double_t radius_min    =      -700.0 - 0.5;
  const Double_t radius_max    =      +700.0 + 0.5;

  const Int_t    N_height      =                18;
  const Double_t height_min    =              20.0;
  const Double_t height_max    =             200.0;

  const Int_t    N_ToTfrac     =                20;
  const Double_t ToTfrac_min   =               0.0;
  const Double_t ToTfrac_max   =               1.0;

  const Int_t    N_dt          =              1500;
  const Double_t dt_min        =      -500.0 - 0.5;
  const Double_t dt_max        =      1000.0 - 0.5;

  const Int_t    N_zenith      =                21;
  const Double_t zenith_min    =              -1.0;
  const Double_t zenith_max    =              +1.0;
  
  const Int_t    N_energy      =                60;
  const Double_t energy_min    =              -2.0;
  const Double_t energy_max    =               4.0;

  const Int_t    N_quality     =               100;
  const Double_t quality_min   =            -200.0;
  const Double_t quality_max   =            +800.0;

  const Int_t    N_beta0       =                60;
  const Double_t beta0_min     =              -2.0;
  const Double_t beta0_max     =              +3.5;
  

  // Set of histograms for data/MC-comparisons

  string ordinateLabel(weightOption == JEvtWeightOption::unweighted() ? 
		       "Number of events" : "Rate [Hz]");

  TFile out(outputFile.c_str(), "recreate");

  TH1D hs  ("hs",   MAKE_CSTRING("; #snapshot hits; "                                  << ordinateLabel),
	    N_Nhits,     Nhits_min,     Nhits_max);
  TH1D hn  ("hn",   MAKE_CSTRING("; #triggered hits; "                                 << ordinateLabel),
	    N_Nhits,     Nhits_min,     Nhits_max);
  TH1D ho  ("ho",   MAKE_CSTRING("; #overlays; "                                       << ordinateLabel),
	    N_Noverlays, Noverlays_min, Noverlays_max);
  
  TH1D hz  ("hz",   MAKE_CSTRING("; cos(#theta); "                                     << ordinateLabel),
	    N_zenith,    zenith_min,    zenith_max);
  TH1D he  ("he",   MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); "                      << ordinateLabel),
	    N_energy,    energy_min,    energy_max);
  TH1D ht  ("ht",   MAKE_CSTRING("; #Delta t [ns]; "                                   << ordinateLabel),
	    N_dt,        dt_min,        dt_max);

  TH1D hZ  ("hZ",   MAKE_CSTRING("; longitudinal ToT-barycenter [m]; "                 << ordinateLabel),
	    N_height,    height_min,    height_max);
  TH1D hT0 ("hT0",  MAKE_CSTRING("; ToT-fraction below earliest triggered hit [ns]; "  << ordinateLabel),
	    N_ToTfrac,   ToTfrac_min,   ToTfrac_max);
  TH1D hT1 ("hT1",  MAKE_CSTRING("; ToT-fraction above earliest triggered hit [ns]; "  << ordinateLabel),
	    N_ToTfrac,   ToTfrac_min,   ToTfrac_max);  

  TH1D hq  ("hq",   MAKE_CSTRING("; quality; "                                         << ordinateLabel),
	    N_quality,   quality_min,   quality_max);
  TH1D hb0 ("hb0",  MAKE_CSTRING("; #beta_{0}; "                                       << ordinateLabel),
	    N_beta0,     beta0_min,     beta0_max);
  
  TH2D hzo ("hzo",  MAKE_CSTRING("; cos(#theta); #overlays; "                          << ordinateLabel),
	    N_zenith,    zenith_min,    zenith_max,
	    N_Noverlays, Noverlays_min, Noverlays_max);
  TH2D hze ("hze",  MAKE_CSTRING("; cos(#theta); log_{10}(E_{reco} / 1 GeV); "         << ordinateLabel),
	    N_zenith,    zenith_min,    zenith_max,
	    N_energy,    energy_min,    energy_max);
  TH2D heo ("heo",  MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #overlays; "           << ordinateLabel),
	    N_energy,    energy_min,    energy_max,
	    N_Noverlays, Noverlays_min, Noverlays_max);
  TH2D hzt ("hzt",  MAKE_CSTRING("; cos(#theta); #Delta t [ns]; "                      << ordinateLabel),
	    N_zenith,    zenith_min,    zenith_max,
	    N_dt,        dt_min,        dt_max);
  TH2D het ("het",  MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; "       << ordinateLabel),
	    N_energy,    energy_min,    energy_max,
	    N_dt,        dt_min,        dt_max);
  TH2D hot ("hot",  MAKE_CSTRING("; #overlays; #Delta t [ns]; "                        << ordinateLabel),
	    N_Noverlays, Noverlays_min, Noverlays_max,
	    N_dt,        dt_min,        dt_max);
  
  TH2D hzq ("hzq",  MAKE_CSTRING("; cos(#theta); quality; "                            << ordinateLabel),
	    N_zenith,    zenith_min,    zenith_max,
	    N_quality,   quality_min,   quality_max);
  TH2D heq ("heq",  MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); quality; "             << ordinateLabel),
	    N_energy,    energy_min,    energy_max,
	    N_quality,   quality_min,   quality_max);
  TH2D hoq ("hoq",  MAKE_CSTRING("; #overlays; quality; "                              << ordinateLabel),
	    N_Noverlays, Noverlays_min, Noverlays_max,
	    N_quality,   quality_min,   quality_max);
  
  TH2D hob0("hob0", MAKE_CSTRING("; #overlays; log_{10}(#beta_{0}); "                  << ordinateLabel),
	    N_Noverlays, Noverlays_min, Noverlays_max,
	    N_beta0,     beta0_min,     beta0_max);
  TH2D heb0("heb0", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
	    N_energy,    energy_min,    energy_max,
	    N_beta0,     beta0_min,     beta0_max);     
  TH2D hzb0("hzb0", MAKE_CSTRING("; cos(#theta); log_{10}(#beta_{0}); "                << ordinateLabel),
	    N_zenith,    zenith_min,    zenith_max,
	    N_beta0,     beta0_min,     beta0_max);
  
  TH2D hxy ("hxy",  MAKE_CSTRING("; x [m]; y [m]; "                                    << ordinateLabel),
	    N_radius,    radius_min,    radius_max,
	    N_radius,    radius_min,    radius_max);
  

  JEvtWeightFileScannerSet<> scanners(inputFiles);
  
  STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
  
  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
    
    scanner->setLimit(numberOfEvents);
    
    while (scanner->hasNext()) {
      
      const Evt*    evt = scanner->next();

      if (has_reconstructed_jppmuon(*evt)) {
      
	DEBUG(get_best_reconstructed_jppmuon(*evt).comment << endl);
	
	const Trk&      best_trk = get_best_reconstructed_jppmuon(*evt);
	const JTrack3E& tb       = getTrack(best_trk);

	const double    logE           = log10(tb.getE());
	const size_t    NsnapshotHits  = evt->hits.size();
	const size_t    NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
					      make_predicate(&Hit::trig, (long long unsigned int) 0, JComparison::gt()));
	
	if (triggeredHitsRange(NtriggeredHits) &&
	    coszenithRange    (tb.getDZ())     &&
	    energyRange       (tb.getE())) { // Apply cuts

	  // Extract event-weight

	  double weight = 0.0;
	  
	  if        (weightOption == JEvtWeightOption::unweighted()) {
	    weight = 1.0;
	  } else if (weightOption == JEvtWeightOption::standard()) {
	    weight = scanner->getWeight(*evt);
	  } else if (weightOption == JEvtWeightOption::rescaled()) {
	    weight = evt->w.at(WEIGHTLIST_RESCALED_EVENT_RATE);
	  } 
      
	  STATUS(RIGHT     (15)    << distance(scanners.begin(), scanner) <<
		 RIGHT     (15)    << scanner->getCounter()               <<
		 SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);

	  // Fill histograms

	  hs .Fill(NsnapshotHits,                          weight);
	  hn .Fill(NtriggeredHits,                         weight);
	  ho .Fill(evt->overlays,                          weight);
	
	  hz .Fill(tb.getDZ(),                             weight);
	  he .Fill(logE,                                   weight);
	  hq .Fill(best_trk.lik,                           weight);

	  hzo.Fill(tb.getDZ(),            evt->overlays,   weight);
	  hze.Fill(tb.getDZ(),            logE,            weight);
	  heo.Fill(logE,                  evt->overlays,   weight);

	  hoq.Fill(evt->overlays,         best_trk.lik,    weight);
	  hzq.Fill(tb.getDZ(),            best_trk.lik,    weight);
	  heq.Fill(logE,                  best_trk.lik,    weight);

	  hxy.Fill(tb.getX(),             tb.getY(),       weight);

	  if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {

	    const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
	  
	    hb0 .Fill(logb0,                               weight);
	    hob0.Fill(evt->overlays,      logb0,           weight);	  
	    heb0.Fill(logE,               logb0,           weight);
	    hzb0.Fill(tb.getDZ(),         logb0,           weight);
	  }

	  double ToTbelow     = 0.0;
	  double ToTabove     = 0.0;
	  double ToTtotal     = 0.0;
	  double ToTweightedZ = 0.0;
	
	  vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
							     make_comparator(&Hit::t));
	
	  for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
	  
	    if (router.hasPMT(hit->pmt_id)) {
	    
	      const JPMT&     pmt  = router.getPMT(hit->pmt_id);
	    
	      const double    t0   = tb.getT(pmt);
	      const double    t1   = getTime(*hit);
	    
	      ht .Fill(t1 - t0,                            weight);
	      hot.Fill(evt->overlays,     t1 - t0,         weight);
	      hzt.Fill(tb.getDZ(),        t1 - t0,         weight);
	      het.Fill(log10(best_trk.E), t1 - t0,         weight);
	      
	      if (hit->trig > 0) {

		ToTtotal     += hit->tot;
		ToTweightedZ += hit->tot * hit->pos.z;
		
		if        (hit->pos.z < firstHit->pos.z)  {
		  ToTbelow   += hit->tot;
		} else if (hit->pos.z >= firstHit->pos.z) {
		  ToTabove   += hit->tot;
		}
	      }
	    }
	  }
	  
	  hT0.Fill(ToTbelow / ToTtotal,                     weight);
	  hT1.Fill(ToTabove / ToTtotal,                     weight);
	  hZ .Fill(ToTweightedZ / ToTtotal,                 weight);
	}
      }
    }
  }
  
  STATUS(endl);

  out.Write();
  out.Close();
}