#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 "JPhysics/JConstants.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"

#include <string>

#include "JSupernova/JKexing2D.hh"

#include "JROOT/JRootToolkit.hh"
#include "JMath/JMathToolkit.hh"

/**
 * \file
 *
 * Auxiliary program to build supernova background from JKexing2D output
 * \author mlincett
 */
int main(int argc, char **argv) {
  using namespace std;
  using namespace JPP;

  string      inputFile;
  string      outputFile;
  string      summaryFile;
  int         windowSize;
  double      fractionThreshold;
  int         debug;
  JRange<int> multiplicity;

  try { 

    JParser<> zap("Auxiliary program to build supernova background from JKexing2D output");
    
    zap['f'] = make_field(inputFile,   "input file (JKexing2D).");
    zap['o'] = make_field(outputFile,  "output file.") = "sn_background.root";
    zap['w'] = make_field(windowSize,  "size of the sliding window to test") = 5;
    zap['M'] = make_field(multiplicity, "final multiplicity range") = JRange<int>(7,11);
    zap['F'] = make_field(fractionThreshold, "minimum fraction of active channels to compute distribution") = 0.99;
    zap['d'] = make_field(debug) = 1; 

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

  typedef JManager<int,    TH1D>  JManager_i1D_t;
  typedef JManager<string, TH1D>  JManager_s1D_t;
  typedef JManager<string, TH2D>  JManager2D_t;

  // input
  
  TFile in(inputFile.c_str(), "exist");

  TParameter<int>* runNumber;

  in.GetObject("RUNNR", runNumber);

  JManager2D_t   MT =   JManager2D_t::Read(in, mul_p   , '%');
  JManager_s1D_t ST = JManager_s1D_t::Read(in, status_p, '%');

  in.Close();

  // output
  
  const int    factoryLimit_peak = 250;
  const int    factoryLimit_runs = 10000;

  // distributions of individual multiplicities

  vector<string> filters = { "RAW", "TAV", "TA1" };

  map<string, JManager_i1D_t> mmap;
  
  for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
    string title = "MD_" + (*f) + "_%";
    mmap[*f] = JManager_i1D_t(new TH1D(title.c_str(), NULL, factoryLimit_peak, 0, factoryLimit_peak)); 
  }

  const int    nx   = 1 + NUMBER_OF_PMTS; 
  const double xmin = -0.5;
  const double xmax = nx - 0.5;

  // distributions of the trigger level

  JManager_s1D_t TD(new TH1D(Form("SNT_[%d,%d]_", multiplicity.first, multiplicity.second) + TString("%"), NULL, factoryLimit_peak, -0.5, -0.5 + factoryLimit_peak));

  double epsilon = 1e-4;

  // bioluminescence benchmarking
  JManager2D_t  BL(new TH2D("BL_%", NULL, 100, epsilon, 1 + epsilon, factoryLimit_peak, -0.5, -0.5 + factoryLimit_peak));
  JManager2D_t  MUL_EFF(new TH2D("MUL_EFF_%", NULL, 100, epsilon, 1 + epsilon, nx, xmin, xmax));

  // history (value vs run number)
  JManager_s1D_t H(new TH1D("H_%", NULL, factoryLimit_runs, 0, factoryLimit_runs));

  // livetime vs fraction of active channels

  TH1D* LT = new TH1D("LIVETIME_ACF", NULL, 100, epsilon, 1 + epsilon);

  // fraction of active PMT

  TGraph* activeChannelFraction = histogramToGraph(*ST["PMT"]);

  const int       nb        = activeChannelFraction->GetN();
  const Double_t* arr_acf_y = activeChannelFraction->GetY();

  const vector<Double_t> vec_acf_y(arr_acf_y, arr_acf_y + nb);

  LT->FillN(nb, arr_acf_y, NULL);

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

    map<string, TGraph*> count_vs_frame;

    // individual multiplicities

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

      // number of events with multiplicity M as a function of the frame index

      const TH1* px = projectHistogram(*MT[*f], M, M, 'x');

      count_vs_frame[*f] = histogramToGraph(*px);

      delete px;

      // distribution of the number of events with multiplicity M
      // selection criterion based on the minimum fraction of active PMTs

      vector<double> select;
      vector<double> threshold(nb, fractionThreshold);

      // stores in 'select' the result of the element-wise comparison between 'vec_acf_y' and 'threshold'

      transform(vec_acf_y.begin(),
		vec_acf_y.end(),
		threshold.begin(),
		back_inserter(select),
		greater_equal<double> {});
      
      mmap[*f][M]->FillN(nb, count_vs_frame[*f]->GetY(), &select[0]);

      // multiplicity spectrum as a function of number of active PMTs

      const vector<double> nb_M(nb, M);

      MUL_EFF[*f]->FillN(nb,
			 arr_acf_y, 
			 &nb_M[0],
			 count_vs_frame[*f]->GetY());

    }
  }

  int RUNNR = runNumber->GetVal();

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

    // distribution of number of events per frame in the selected multiplicity range (A)

    const TH1* px = projectHistogram(*MT[*f], multiplicity.getLowerLimit(), multiplicity.getUpperLimit(), 'x');

    TGraph* events_per_frame = histogramToGraph(*px);

    delete px;

    TD[*f]->FillN(events_per_frame->GetN(), events_per_frame->GetY(), NULL);

    // distribution of number of events per frame in the selected multiplicity range (A)
    // sliding window applied

    vector<double> vec(events_per_frame->GetY(), events_per_frame->GetY() + nb);
    vector<double> ker(windowSize, 1.0);

    vector<double> events_per_frame_sliding = convolve(vec, ker);

    TD[(*f) + "_nTS"]->FillN(events_per_frame_sliding.size(),
			     &events_per_frame_sliding[0],
			     NULL);
    
    // distribution of number of events per frame in the selected multiplicity range (A)
    // as a function of the fraction of active channels
    BL[*f]->FillN(nb,
		  arr_acf_y,
		  events_per_frame->GetY(),
		  NULL, 1);

    // peak as a function of the run number

    px = projectHistogram(*MT[*f], multiplicity.getLowerLimit(), multiplicity.getUpperLimit(), 'x');

    int peak = px->GetMaximum();

    delete px;

    H["PK_" + (*f)]->Fill(RUNNR, peak);

  }

  // history of bioluminescence

  double BI = (ST["PMT"]->GetSumOfWeights() / ST["PMT"]->GetEntries());

  H["BIOLUM"]->Fill(RUNNR, BI);

  // write output file

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

  TD.Write(out);
  
  TDirectory* bm = out.mkdir("BIOLUM");
  TDirectory* hs = out.mkdir("HISTORY");  

  H.Write(*hs);
  BL.Write(*bm);
  
  MUL_EFF.Write(*bm);

  for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
    string dir_name = "MUL_" + (*f); 
    TDirectory* dir = out.mkdir(dir_name.c_str());
    mmap[*f].Write(*dir);
  }

  LT->Write();

  out.Close();

}