#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TKey.h"
#include "TString.h"

#include "km3net-dataformat/online/JDAQHeader.hh"

#include "JROOT/JRootFileReader.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JSupport/JMeta.hh"

#include "JCalibrate/JCalibrateK40.hh"

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

namespace {

  double LIVETIME_S       =  10.0;        //!< minimal live time [s]
  double ERROR            =   1.0;        //!< default error [Hz/ns]
}


/**
 * \file
 * Auxiliary program to merge K40 data.
 *
 * \author mkarel, mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  vector<string>     inputFile;
  string             outputFile;
  int                debug;

  try { 

    JProperties properties;

    properties.insert(gmake_property(LIVETIME_S));
    properties.insert(gmake_property(ERROR));

    JParser<> zap("Auxiliary program to merge K40 data.");
    
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    zap['f'] = make_field(inputFile,         "input file (output from JCalibrateK40).");
    zap['o'] = make_field(outputFile,        "output file (input to JFitK40).")   = "mergek40.root";
    zap['d'] = make_field(debug,             "debug.")                            = 1;

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


  gErrorIgnoreLevel = kError;

  TH1* weights_hist = NULL;

  typedef map<TString, TH1*>  map_type;

  map_type zmap;


  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    DEBUG("Processing " << *i << endl) ;

    TFile in(i->c_str(), "read");
    
    // check if file contains weights vector

    if (!in.GetListOfKeys()->Contains(weights_hist_t)) {
      FATAL("Histogram does not contain weights histogram." << endl);
    }

    TH1* weights_hist_i = (TH1*) in.Get(weights_hist_t);

    if (weights_hist == NULL) {

      weights_hist = (TH1*) weights_hist_i->Clone();

    } else {

      // check normalization

      TAxis* x_axis = weights_hist->GetXaxis();

      if (weights_hist_i->GetBinContent(x_axis->FindBin(WS_t)) != 
	  weights_hist  ->GetBinContent(x_axis->FindBin(WS_t))) {
	FATAL("Number of y-bins does not match for this file " << *i << endl);
      }

      if (weights_hist_i->GetBinContent(x_axis->FindBin(WB_t)) != 
	  weights_hist  ->GetBinContent(x_axis->FindBin(WB_t))) {
        FATAL("TMax_ns is not the same for this file " << *i << endl);
      }

      // add normalization

      weights_hist->SetBinContent(x_axis->FindBin(W1_overall_t), 
				  weights_hist  ->GetBinContent(x_axis->FindBin(W1_overall_t)) + 
				  weights_hist_i->GetBinContent(x_axis->FindBin(W1_overall_t)));
    }


    TIter iter(in.GetListOfKeys());

    for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {

      if (TString(key->GetName()).EndsWith(_2S) ||
	  TString(key->GetName()).EndsWith(_1B) ||
	  TString(key->GetName()).EndsWith(_1L)) {

        TH1* h1 = dynamic_cast<TH1*>(key->ReadObj());

	map_type::iterator p = zmap.find(h1->GetName());

	if (p == zmap.end()) {

	  DEBUG("Clone " << h1->GetName() << endl);

	  p = zmap.insert(make_pair(h1->GetName(), (TH1*) h1->Clone())).first;

	} else {

	  DEBUG("Add   " << h1->GetName() << endl);

	  p->second->Add(h1);
	}
      }
    }

    weights_hist->SetDirectory(0);
    
    for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
      i->second->SetDirectory(0);
    }
    
    in.Close();
  }


  if (weights_hist == NULL) {
    FATAL("Missing histogram " << weights_hist_t << endl);
  }


  // write file

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

  const Double_t WS = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WS_t)) ;  // [ns]
  const Double_t WB = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WB_t)) ;  // [ns]

  for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {

    if (i->first.EndsWith(_2S)) {

      TH2D* h2s = (TH2D*) i->second;

      if (h2s->GetEntries() != 0) {

	TH1D* h1b = (TH1D*) zmap.find(TString(i->first).ReplaceAll(_2S, _1B))->second;
	TH1D* h1L = (TH1D*) zmap.find(TString(i->first).ReplaceAll(_2S, _1L))->second;

	//----------------------------------------------------------
	// normalisation and background estimation
	//----------------------------------------------------------

	for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {

	  // normalised background estimation
	  // if livetime is compatible with zero then set bin contents
	  // and erros to zero and default value, respectively

	  if (h1L->GetBinContent(ix) <= LIVETIME_S) {

	    h1b->SetBinContent(ix, 0.0);
	    h1b->SetBinError  (ix, 0.0);

	    for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
	      h2s->SetBinContent(ix, iy, 0.0);
	      h2s->SetBinError  (ix, iy, ERROR);
	    }

	  } else {

	    const double Y = h1b->GetBinContent(ix) / h1L->GetBinContent(ix) / WB;
	    const double W = h1b->GetBinError  (ix) / h1L->GetBinContent(ix) / WB;

	    h1b->SetBinContent(ix, Y);
	    h1b->SetBinError  (ix, W);

	    for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {

	      double y = h2s->GetBinContent(ix, iy) / h1L->GetBinContent(ix) / WS;
	      double w = h2s->GetBinError  (ix, iy) / h1L->GetBinContent(ix) / WS;

	      if (w == 0.0) {

		// for chi2 minimisation; accordingly set error of empty bins

		w = 1.0 / h1L->GetBinContent(ix) / WS;
	      }

	      h2s->SetBinContent(ix, iy, y - Y);
	      h2s->SetBinError  (ix, iy, sqrt(w*w + W*W));
	    }
	  } 
	}

	// reset under and overflows.

	for (int ix = 0; ix <= h2s->GetXaxis()->GetNbins() + 1; ++ix) {
	  h2s->SetBinContent(ix, 0, 0.0);  h2s->SetBinContent(ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
	  h2s->SetBinError  (ix, 0, 0.0);  h2s->SetBinError  (ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
	}

	for (int iy = 0; iy <= h2s->GetYaxis()->GetNbins() + 1; ++iy) {
	  h2s->SetBinContent(0, iy, 0.0);  h2s->SetBinContent(h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
	  h2s->SetBinError  (0, iy, 0.0);  h2s->SetBinError  (h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
	}

	NOTICE("Biased coincidence rate [Hz] " << h2s->GetName() << ' ' << h2s->GetSumOfWeights() * WS << endl);

	h2s->SetName(TString(i->first).ReplaceAll(_2S, _2R));

	h2s->Write();
      }
    } 
  }

  putObject(out, JMeta(argc, argv));

  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    JMeta::copy(i->c_str(), out);

    for (JRootFileReader<JDAQHeader> in(i->c_str()); in.hasNext(); ) {
      putObject(out, *in.next());
    }
  }

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