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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF2.h"
#include "TMath.h"
#include "TFitResult.h"

#include "JROOT/JMinimizer.hh"

#include "JPhysics/JConstants.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMeta.hh"

#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JFitK40.hh"

#include "JROOT/JManager.hh"

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

const char* livetime_t = "LiveTime";

/**
 * \file
 *
 * Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.
 * \author mlincett
 */

int main(int argc, char **argv) {
  using namespace std;
  using namespace JPP;

  string    inputFile;
  string    outputFile;
  string    summaryFile;
  string    option;
  int       debug;
  bool      fitBackground;
  bool      nominalLiveTime;

  try { 
    JParser<> zap("Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.");
    
    zap['f'] = make_field(inputFile,         "input file");
    zap['o'] = make_field(outputFile,        "output file")                                           = "fitmultiplicity.root";
    zap['O'] = make_field(option,            "fitting option")                                        = "Q";
    zap['B'] = make_field(fitBackground,     "get signal by fit and subtraction from background from the integral");
    zap['L'] = make_field(nominalLiveTime,   "use nominal live time from LiveTime histogram, instead of individual DOMs");
    zap['d'] = make_field(debug) = 1;   

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

  TFile* in = TFile::Open(inputFile.c_str(), "exist");

  if (in == NULL || !in->IsOpen()) {
    FATAL("File: " << inputFile << " not opened." << endl);
  }

  
  //----------------------------------------------------------
  // general histograms
  //----------------------------------------------------------

  const int    nx   = 1 + NUMBER_OF_PMTS; 
  const double xmin = -0.5;
  const double xmax = nx - 0.5;
  
  typedef JManager<TString, TH1D> JManager_t;

  // inclusive
  JManager_t HTF(new TH1D("%F", "",  nx, xmin, xmax));
  // exclusive
  JManager_t HTFE(new TH1D("%FE", "", nx, xmin, xmax));

  HTF->Sumw2(kFALSE);

  //----------------------------------------------------------
  // get overall normalization constants
  //----------------------------------------------------------

  TH1D* hLiveTime = (TH1D*) in->Get(livetime_t);

  if (hLiveTime == NULL) {
    FATAL("File " << inputFile << " does not contain livetime histogram." << endl);
  }
  
  //----------------------------------------------------------
  // initialise output file
  //----------------------------------------------------------

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

  //----------------------------------------------------------
  // loop over bins of livetime histogram
  //----------------------------------------------------------

  TDirectory *fdir = out.mkdir("FIT");
  TDirectory *pdir = out.mkdir("MUL"); 

  // nominal live time
  double nT = hLiveTime->GetBinContent(1);

  for (int bin = 2; bin <= hLiveTime->GetNbinsX(); bin++) { 
 
    TString prefix = hLiveTime->GetXaxis()->GetBinLabel(bin);

    DEBUG("Processing data from " << prefix << endl);

    double T = nominalLiveTime ? nT : hLiveTime->GetBinContent(bin);

    TString h1D_label = prefix + TString("_P" );
    TString h2D_label = prefix + TString("_HT");

    if (T > 0) {
     
      TH2D* h1D = (TH2D*) in->Get(h1D_label);
      TH2D* h2D = (TH2D*) in->Get(h2D_label);

      if (h1D != NULL) {

	h1D->Scale(1.0 / T);

	pdir->cd();

	h1D->Write();
	
      } else {
	DEBUG(h1D_label << " not found" << endl);
      }

      if (h2D != NULL) {

	TDirectory* cdir = fdir->mkdir(h2D_label); cdir->cd();     

	for (int M = 2; M <= NUMBER_OF_PMTS; M++) {

	  int Mbin = h2D->GetXaxis()->FindBin(M);

	  TH1D* hI = h2D->ProjectionY(h2D_label + Form("_%d", M), Mbin, NUMBER_OF_PMTS);

	  hI->Sumw2();

	  if (hI->GetEntries() > 0) {
	    
	    const double W = (hI->GetXaxis()->GetXmax() - hI->GetXaxis()->GetXmin()); 
	    const double N = hI->GetNbinsX();
	    const double V = W / N;
	    const double minRMS = 0.1;

	    double inclusiveCountVal = 0, inclusiveCountErr = 0;
	    
	    if (hI->GetRMS() > minRMS) {

	      hI->Scale(1.0 / V);

	      TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");

	      f1.SetParameter(0, hI->GetMaximum());
	      f1.SetParameter(1, hI->GetMean());
	      f1.SetParameter(2, hI->GetRMS() * 0.25);
	      f1.SetParameter(3, hI->GetMinimum());
	      
	      hI->Fit(&f1, option.c_str(), "same");

	      double sg_v = f1.GetParameter(0);
	      double sg_e = f1.GetParError( 0);
	      
	      double bg_v = f1.GetParameter(3);
	      double bg_e = f1.GetParError( 3);

	      // 2-sigma acceptance criteria for fitted values
	      if (bg_v / bg_e <= 2.0) { bg_v = bg_e = 0; }
	      if (sg_v / sg_e <= 2.0) { sg_v = sg_e = 0; }

	      double integral = hI->GetSumOfWeights();

	      if (fitBackground) {

		inclusiveCountVal = V * (integral - (bg_v * N));
		inclusiveCountErr = V * bg_e * N;

	      } else {

		inclusiveCountErr = sg_e;
		inclusiveCountVal = sg_v;

	      } 	      

	      hI->Scale(1.0 / T);

	      // re-fit after scaling
	      hI->Fit(&f1, option.c_str(), "same");	    
	    
	    }
				       
	    if (inclusiveCountVal <= 0) { inclusiveCountVal = 0; }

	    // --------------------------------
	    // write fitted histogram
	    // --------------------------------

	    hI->Write();

	    // --------------------------------
	    // build inclusive rates histograms
	    // --------------------------------

	    int bin = HTF[h2D_label]->FindBin(M);

	    HTF[h2D_label]->SetBinContent(bin, inclusiveCountVal);

	    double binError = sqrt(inclusiveCountVal + pow(inclusiveCountErr, 2));

	    HTF[h2D_label]->SetBinError(bin, binError);
	    
	    // --------------------------------
	    // build exclusive rates histograms
	    // --------------------------------

	    double exclusiveCountVal = 0;

	    if (M == NUMBER_OF_PMTS) {
	      exclusiveCountVal = inclusiveCountVal;
	    } else {
	      int prevBin = HTF[h2D_label]->GetXaxis()->FindBin(M+1);
	      exclusiveCountVal = inclusiveCountVal - HTF[h2D_label]->GetBinContent(prevBin); 
	    }
	       
	    HTFE[h2D_label]->SetBinContent(HTFE[h2D_label]->FindBin(M), exclusiveCountVal);
	      
	  
	  }
	  
	}

      } else {
	DEBUG(h1D_label << " not found" << endl);
      }	

      // HTF[hLabel]->Sumw2(kTRUE);
      	
      HTFE[h2D_label]->Sumw2();

      HTF[h2D_label ]->Scale(1.0 / T);
      HTFE[h2D_label]->Scale(1.0 / T);
      
    }
  }

  fdir = out.mkdir("HTF" ); fdir->cd();

  HTF.Write(*fdir);

  fdir = out.mkdir("HTFE"); fdir->cd();

  HTFE.Write(*fdir);

  out.Close();

}