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

#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH2.h"
#include "TF2.h"
#include "TProfile2D.h"
#include "TString.h"
#include "TRegexp.h"

#include "JGizmo/JRootObjectID.hh"
#include "JGizmo/JGizmoToolkit.hh"

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


/**
 * \file
 * Auxiliary program for operations on 2D histograms.\n
 * The option <tt>-f</tt> corresponds to <tt>\<file name\>:\<object name\></tt>.
 *
 * Possible operations (option <tt>-u \<operation\></tt>:
 *     - Add\n
 *       ROOT addition of histograms;
 *     - add\n
 *       Addition of histograms by look up of corresponding bin in second histogram;
 *     - Subtract\n
 *       ROOT subtraction of histograms;
 *     - subtract\n
 *       Subtraction of histograms by look up of corresponding bin in second histogram;
 *     - Multiply\n
 *       ROOT multiplication of histograms;
 *     - multiply\n
 *       Multiplication of histograms by look up of corresponding bin in second histogram;
 *     - Divide\n
 *       ROOT division of histograms;
 *     - divide\n
 *       Division of histograms by look up of corresponding bin in second histogram;
 *     - efficiency\n
 *       Like "Divide", but with correction for errors;
 *     - stdev\n
 *       Standard deviation;
 *     - Replace\n
 *       Replace histogram contents by function value;
 *
 * If only one histogram is specified, the operation applies to the histogram and the first associated function.\n
 * Only the quoted ROOT operations and the additional operation Replace are then allowed.\n
 * The option Replace will replace the histogram contents by the function value.
 * 
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  
  vector<JRootObjectID> inputFile;
  string                outputFile;
  string                opera;
  string                name;
  int                   debug;

  try {

    JParser<> zap("Auxiliary program for histogram operations.");

    zap['f'] = make_field(inputFile,  "<input file>:<object name>");
    zap['u'] = make_field(opera)  = 
      JOpera::Add(), 
      JOpera::add(), 
      JOpera::Subtract(), 
      JOpera::subtract(), 
      JOpera::Multiply(), 
      JOpera::multiply(),  
      JOpera::Divide(), 
      JOpera::divide(), 
      JOpera::efficiency(), 
      JOpera::stdev(),
      JOpera::Replace();
    zap['o'] = make_field(outputFile, "ROOT file with histogram(s)")    = "opera.root";
    zap['O'] = make_field(name,       "histogram name:"
			  << "\n\t\"" << JOpera::SAME_AS_OPERATION() << "\" -> same as operation; or"
			  << "\n\t\"" << JOpera::SAME_AS_INPUT()     << "\" -> same as input; else"
			  << "\n\t as specified")                       = JOpera::SAME_AS_OPERATION();
    zap['d'] = make_field(debug)      = 1;

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


  vector<TObject*> listOfObjects;

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

    DEBUG("Input: " << *input << endl);

    TDirectory* dir = getDirectory(*input);

    if (dir == NULL) {
      ERROR("File: " << input->getFullFilename() << " not opened." << endl);
      continue;
    }

    const TRegexp regexp(input->getObjectName());

    TIter iter(dir->GetListOfKeys());

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

      const TString tag(key->GetName());

      DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);

      // option match

      if (tag.Contains(regexp) && isTObject(key)) {

        TObject* p = key->ReadObj();

        if (dynamic_cast<TH2*>(p) == NULL) {
          FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
        }

        listOfObjects.push_back(p);
      }
    }
  }


  TH2* h3 = NULL;

  if        (listOfObjects.size() == 0) {

    FATAL("Number of histograms " << listOfObjects.size() << " = 0." << endl);

  } else if (listOfObjects.size() == 1 && (opera == JOpera::Add()      ||
                                           opera == JOpera::Subtract() ||
                                           opera == JOpera::Multiply() ||
                                           opera == JOpera::Divide()   ||
                                           opera == JOpera::Replace())) {

    // operation between histogram contents and first associated function

    TH2* h1 = dynamic_cast<TH2*>(listOfObjects[0]);
    TF2* f1 = dynamic_cast<TF2*>(h1->GetListOfFunctions()->First());

    if (f1 == NULL) {
      FATAL(h1->GetName() << " has no associated function." << endl);
    }

    NOTICE(h1->GetName() << ' ' << opera << ' ' << f1->GetName() << endl);

    if      (name == JOpera::SAME_AS_OPERATION())
      h3 = (TH2*) h1->Clone(opera.c_str());
    else if (name == JOpera::SAME_AS_INPUT())
      h3 = (TH2*) h1->Clone(h1->GetName());
    else
      h3 = (TH2*) h1->Clone(name.c_str());

    if        (opera == JOpera::Add()) {

      h3->Add     (f1, +1.0);

    } else if (opera == JOpera::Subtract()) {

      h3->Add     (f1, -1.0);

    } else if (opera == JOpera::Multiply()) {

      h3->Multiply(f1);

    } else if (opera == JOpera::Divide()) {

      h3->Divide  (f1);
    
    } else if (opera == JOpera::Replace()) {

      h3->Reset();
      h3->Fill(0.0, 0.0, 0.0);     // ensure one entry for drawing after reset
      h3->Add     (f1);
    }

  } else if (listOfObjects.size() == 2) {

    TH2* h1 = dynamic_cast<TH2*>(listOfObjects[0]);
    TH2* h2 = dynamic_cast<TH2*>(listOfObjects[1]);

    NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);

    if      (name == JOpera::SAME_AS_OPERATION())
      h3 = (TH2*) h1->Clone(opera.c_str());
    else if (name == JOpera::SAME_AS_INPUT())
      h3 = (TH2*) h1->Clone(h1->GetName());
    else
      h3 = (TH2*) h1->Clone(name.c_str());

    if        (opera == JOpera::Add()) {

      h3->Add     (h1, h2, +1.0, +1.0);

    } else if (opera == JOpera::Subtract()) {
	
      h3->Add     (h1, h2, +1.0, -1.0);

    } else if (opera == JOpera::Multiply()) {
	
      h3->Multiply(h1, h2, +1.0, +1.0);

    } else if (opera == JOpera::Divide()) {

      h3->Divide  (h1, h2, +1.0, +1.0);

    } else if (opera == JOpera::add()       ||
	       opera == JOpera::subtract()  ||
	       opera == JOpera::multiply()  ||
	       opera == JOpera::divide()) {

      for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
	for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {

	  const Double_t x  = h1->GetXaxis()->GetBinCenter(i1);
	  const Double_t y  = h1->GetYaxis()->GetBinCenter(i2);
	  
	  const Int_t    j1 = h2->GetXaxis()->FindBin(x);
	  const Int_t    j2 = h2->GetYaxis()->FindBin(y);
	    
	  const Double_t w1 = h1->GetBinContent(i1,i2);
	  const Double_t w2 = h2->GetBinContent(j1,j2);

	  Double_t w3 = 0.0;
	    
	  if        (opera == JOpera::add()) {

	    w3 = w1 + w2;

	  } else if (opera == JOpera::subtract()) {

	    w3 = w1 - w2;

	  } else if (opera == JOpera::multiply()) {

	    w3 = w1 * w2;

	  } else if (opera == JOpera::divide()) {

	    if (w2 == 0.0) {
	      ERROR("Bin " << h2->GetName() << "[" << j1 << "," << j2 << "] empty." << endl);
	      continue;
	    }

	    w3 = w1 / w2;
	  }

	  h3->SetBinContent(i1, i2, w3);
	}
      }

    } else if (opera == JOpera::efficiency()  ||
	       opera == JOpera::stdev()       ||
	       opera == JOpera::sqrt()) {

      for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
	for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {

	  const Double_t y1 = h1->GetBinContent(i1,i2);
	  const Double_t y2 = h2->GetBinContent(i1,i2);

	  Double_t w1 = h1->GetBinError(i1,i2);
	  Double_t w2 = h2->GetBinError(i1,i2);
	    
	  Double_t y3 = 0.0;
	  Double_t w3 = 0.0;
	    
	  if        (opera == JOpera::efficiency()) {
	      
	    if (y2 == 0.0) {
	      ERROR("Bin " << h2->GetName() << "[" << i1 << "," << i2 << "] empty." << endl);
	      continue;
	    }
	      
	    y3 = y1 / y2;
	      
	    if (y1 != 0.0 &&
		y2 != 0.0) {
		
	      w1 /= y1;
	      w2 /= y2;

	      w3  = y3 * fabs(y3 - 1.0) * sqrt(w1*w1 + w2*w2);
	    }
	      
	  } else if (opera == JOpera::stdev()) {

	    w3 = sqrt(w1*w1 + w2*w2);

	    if (w3 != 0.0) {
	      y3 = (y1 - y2) / w3;
	      w3 = 0.0;
	    }

	  } else if (opera == JOpera::sqrt()) {

	    y3 = (y1+y2) * (y1-y2);

	    if (y3 < 0.0)
	      y3 = -sqrt(-y3);
	    else
	      y3 = +sqrt(+y3);

	    w3 = 0.0;
	  }
	
	  h3->SetBinContent(i1, i2, y3);
	  h3->SetBinError  (i1, i2, w3);
	}
      }
    }

  } else if (opera == JOpera::Add()      ||
             opera == JOpera::Multiply()) {

    for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {

      TH2* h1 = dynamic_cast<TH2*>(*i);

      NOTICE(h1->GetName() << ' ' << opera << endl);

      if (h3 == NULL) {

	if      (name == JOpera::SAME_AS_OPERATION())
	  h3 = (TH2*) h1->Clone(opera.c_str());
	else if (name == JOpera::SAME_AS_INPUT())
	  h3 = (TH2*) h1->Clone(h1->GetName());
	else
	  h3 = (TH2*) h1->Clone(name.c_str());

      } else {

        if        (opera == JOpera::Add()) {

          h3->Add     (h3, h1, +1.0, +1.0);

        } else if (opera == JOpera::Multiply()) {

          h3->Multiply(h3, h1, +1.0, +1.0);
        }
      }
    }

  } else {

    FATAL("Incompatible number of histograms " << listOfObjects.size() << " with option " << opera << endl);
  }

  if (h3 != NULL) {

    TFile out(outputFile.c_str(), "recreate");
    
    h3->Write();
    
    out.Write();
    out.Close();
  }
}