#ifndef __JGIZMOTOOLKIT__
#define __JGIZMOTOOLKIT__

#include <string>
#include <sstream>
#include <map>
#include <cmath>
#include <memory>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wall"
#include "TError.h"
#include "TFile.h"
#include "TClass.h"
#include "TObject.h"
#include "TKey.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TAttLine.h"
#include "TAttFill.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include "TGraph2D.h"
#include "TGraph2DErrors.h"
#include "TMultiGraph.h"
#include "TLine.h"
#include "TEllipse.h"
#include "TString.h"
#include "TRegexp.h"
#include "TFormula.h"
#include "TF1.h"
#include "TF2.h"
#include "TList.h"
#include "TIterator.h"
#include "TMethod.h"
#include "TMethodCall.h"
#include "TAxis.h"
#include "TMath.h"
#pragma GCC diagnostic pop

#include "JLang/JException.hh"
#include "JGizmo/JRootObjectID.hh"

#include "JDetector/JModuleAddressMap.hh"


/**
 * \author mdejong
 */

namespace JGIZMO {}
namespace JPP { using namespace JGIZMO; }

/**
 * Auxiliary applications for use of ROOT and more.
 */
namespace JGIZMO {

  using JLANG::JParseError;
  using JLANG::JValueOutOfRange;
  using JLANG::JFunctionalException;
  using JDETECTOR::JModuleAddressMap;


  /**
   * Auxiliary data structure for JOpera1D.cc and JOpera2D.cc applications.
   */
  struct JOpera {
    //
    // Histogram name.
    //
    static const char* const SAME_AS_OPERATION() { return "%"; }   //!< Set name of output histogram to name of operation
    static const char* const SAME_AS_INPUT()     { return "="; }   //!< Set name of output histogram to name of input histogram

    //
    // Histogram operations.
    //
    static const char* const Add()        { return "Add"; }        //!< ROOT TH1::Add
    static const char* const add()        { return "add"; }        //!< Add      contents with lookup bin in second histogram
    static const char* const Subtract()   { return "Subtract"; }   //!< ROOT TH1::Subtract
    static const char* const subtract()   { return "subtract"; }   //!< Subtract contents with lookup bin in second histogram
    static const char* const Multiply()   { return "Multiply"; }   //!< ROOT TH1::Multiply
    static const char* const multiply()   { return "multiply"; }   //!< Multiply contents with lookup bin in second histogram
    static const char* const Divide()     { return "Divide"; }     //!< ROOT TH1::Divide
    static const char* const divide()     { return "divide"; }     //!< Divide   contents with lookup bin in second histogram
    static const char* const efficiency() { return "efficiency"; } //!< Divide   contents and multiply errors with inefficiency
    static const char* const stdev()      { return "stdev"; }      //!< Set contents to standard deviation
    static const char* const sqrt()       { return "sqrt"; }       //!< Set contents to signed difference between squares
    static const char* const Replace()    { return "Replace"; }    //!< Set contents to associated function
  };


  /**
   * Time stamp of earliest UTC time.
   */
  static const char* const TIMESTAMP = "#splitline{}{#splitline{%d:%m:%y}{  %H:%M}}%F1970-01-01 00:00:00";


  /**
   * Get TFile pointer corresponding to give file name.
   *
   * The TFile pointer of an already opened file is recovered, else a new file is opened.\n
   * Note that the closure of the opened files should be done by the caller of this method.
   *
   * \param  file_name      file name
   * \param  option         TFile::Open option
   * \return                pointer to TFile
   */
  inline TFile* getFile(const std::string& file_name, const std::string& option = "exist")
  {
    using namespace std;

    gErrorIgnoreLevel = kError;

    static map<string, TFile*> zmap;
    
    map<string, TFile*>::iterator i = zmap.find(file_name);
    
    if (i == zmap.end() || i->second == NULL || !i->second->IsOpen()) {
      
      TFile* file = TFile::Open(file_name.c_str(), option.c_str());
      
      zmap[file_name] = file;

      return file;

    } else {
    
      return i->second;
    }
  }
 

  /**
   * Get TDirectory pointer.
   *
   * The TFile pointer of an already opened file is recovered, else a new file is opened.
   *
   * \param  id             identifier
   * \return                pointer to TDirectory
   */
  inline TDirectory* getDirectory(const JRootObjectID& id)
  {
    TFile* in = getFile(id.getFilename().c_str(), "exist");

    if (in == NULL || !in->IsOpen()) {
      return NULL;
    }

    if (id.getDirectory() != "")
      return in->GetDirectory(id.getDirectory());
    else
      return in;
  }


  /**
   * Get first TObject with given identifier.
   *
   * \param  id             identifier
   * \return                pointer to TObject (or NULL)
   */
  inline TObject* getObject(const JRootObjectID& id)
  {
    TDirectory* dir = getDirectory(id);

    if (dir != NULL) {

      const TRegexp regexp(id.getObjectName());

      TIter iter(dir->GetListOfKeys());

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

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

	// option match

	if (tag.Index(regexp) != -1) {
	  return key->ReadObj();
	}
      }
    }

    return NULL;
  }


  /**
   * Get drawing option of object.
   *
   * \param  object         pointer to TObject
   * \return                true if object looks like a line; else false
   */
  inline bool isTAttLine(const TObject* object)
  {
    {
      if (dynamic_cast<const TH3*>(object) != NULL) {
	return false;
      }
    }
    {
      if (dynamic_cast<const TH2*>(object) != NULL) {
	return false;
      }
    }
    {
      const TH1* h1 = dynamic_cast<const TH1*>(object);
    
      if (h1 != NULL) {

	if (h1->GetSumw2N()) {
	  for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
	    if (h1->GetBinError(i) != 0.0) {
	      return false;
	    }
	  }
	}

	return true;
      }
    }
    {
      const TGraphErrors* g1 = dynamic_cast<const TGraphErrors*>(object);

      if (g1 != NULL) {

	for (Int_t i = 0; i != g1->GetN(); ++i) {
	  if (g1->GetEY()[i] != 0.0) {
	    return false;
	  }
	}

	return g1->GetN() > 1;
      }
    }
    {
      const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(object);

      if (g1 != NULL) {

	for (Int_t i = 0; i != g1->GetN(); ++i) {
	  if (g1->GetEYhigh()[i] != 0.0 ||
	      g1->GetEYlow() [i] != 0.0) {
	    return false;
	  }
	}

	return g1->GetN() > 1;
      }
    }
    {
      const TGraph* g1 = dynamic_cast<const TGraph*>(object);

      if (g1 != NULL) {
	return g1->GetN() > 1;
      }
    }

    return true;
  }


  /**
   * Get drawing option of object.
   *
   * \param  object         pointer to TObject
   * \return                true if object looks like an area; else false
   */
  inline bool isTAttFill(const TObject* object)
  {
    {
      const TAttFill* t1 = dynamic_cast<const TAttFill*>(object);
    
      if (t1 != NULL) {
	return t1->GetFillColor() != 1 && t1->GetFillStyle() != 0;
      }
    }

    return false;
  }


  /**
   * Get result of given textual formula.
   *
   * The formula may contain names of member methods of the object pointed to.\n
   * These methods should return a value that is compatible with <tt>Double_t</tt> and could have arguments.\n
   * For example:
   * <pre>
   *         getResult("1.0 / GetEntries", TH1*);
   * </pre>
   *
   * \param  text           text
   * \param  object         pointer to object
   * \return                value
   */
  inline Double_t getResult(const TString& text, TObject* object = NULL)
  {
    TString buffer(text);
    
    if (object != NULL) {

      TClass* p = TClass::GetClass(object->ClassName());
      
      if (p != NULL) {
	
	for ( ; ; ) {

	  TMethod* method = NULL;

	  for (std::unique_ptr<TIterator> iter(p->GetListOfAllPublicMethods()->MakeIterator()); TMethod* p = (TMethod*) iter->Next(); ) {
	    if (buffer.Index(p->GetName()) != -1) {
	      if (method == NULL || strlen(p->GetName()) > strlen(method->GetName())) {
		method = p;
	      }
	    }
	  }

	  if (method == NULL) {
	    break;
	  }

	  for (Ssiz_t index; (index = buffer.Index(method->GetName())) != -1; ) {

	    const TRegexp fp(" *([^)]*)");   // function call

	    Ssiz_t len;
	    Ssiz_t pos = buffer.Index(fp, &len, index);
	    
	    Double_t value;
	    
	    if (pos == -1 || pos != index + (Ssiz_t) strlen(method->GetName())) {

	      TMethodCall(p, method->GetName(), NULL).Execute(object, value);

	      len  = strlen(method->GetName());

	    } else {

	      TMethodCall(p, method->GetName(), NULL).Execute(object, TString(buffer(pos + 1, len - 2)), value);

	      len += strlen(method->GetName());
	    }

	    buffer.Replace(index, len, TString::Format("%20.10e", value));
	  }
	}
      }
    }

    return TFormula("/tmp", buffer.Data()).Eval(0.0);
  }


  /**
   * Get result of given textual formula.
   *
   * The formula may contain names of member methods of the object pointed to.\n
   * These methods should return a value that is compatible with <tt>Double_t</tt> and could have arguments.\n
   * For example:
   * <pre>
   *         getResult("1.0 / GetEntries", TH1*);
   * </pre>
   *
   * \param  text           text
   * \param  object         pointer to object
   * \return                value
   */
  inline Double_t getResult(const std::string& text, TObject* object = NULL)
  {
    return getResult(TString(text.c_str()), object);
  }


  /**
   * Get parameter number from text string.
   *
   * The number corresponds to the value <tt>[0-9]*</tt> in the expression <tt>"p[0-9]* = .."</tt>.
   *
   * \param  text           text
   * \return                parameter number
   */
  inline int getParameter(const std::string& text)
  {
    const char* regexp("p[0-9]* *=");

    TString buffer(text.c_str());
    
    buffer = buffer(TRegexp(regexp));
    buffer = buffer(1, buffer.Length() - 2);

    if (!buffer.IsDigit()) {
      THROW(JParseError, "Text is not a number " << text << ' ' << regexp);
    }
  
    return buffer.Atoi();
  }


  /**
   * Get parameter value from text string.
   *
   * The formula may contain names of member methods of the object pointed to.\n
   * These methods should return a value that is compatible with <tt>Double_t</tt> and could have arguments.\n
   * For example:
   * <pre>
   *         getValue("p[..] = 2 * GetMaximum", TH1*);
   * </pre>
   *
   * \param  text           text
   * \param  object         pointer to object
   * \return                value
   */
  inline Double_t getValue(const std::string& text, TObject* object = NULL)
  {
    const char* regexp("=.*");

    TString buffer(text.c_str());
    
    buffer = buffer(TRegexp(regexp));
    buffer = buffer(1, buffer.Length() - 1);

    return getResult(std::string(buffer), object);
  }


  /**
   * Get parameter value from text string.
   *
   * The formula may contain names of member methods of the object pointed to.\n
   * These methods should return a value that is compatible with <tt>Double_t</tt> and could have arguments.\n
   * For example:
   * <pre>
   *         getValue("p[..] = 1.0 2.0 3.0", 1);
   * </pre>
   * will return <tt>2.0</tt>.
   *
   * \param  text           text
   * \param  index          index
   * \return                value
   */
  inline Double_t getValue(const std::string& text, const int index)
  {
    using namespace std;

    const char* regexp("=.*");

    TString buffer(text.c_str());

    buffer = buffer(TRegexp(regexp));
    buffer = buffer(1, buffer.Length() - 1);


    istringstream is((std::string(buffer)));

    Double_t value;

    for (int i = index; is >> value && i > 0; --i) {}

    if (is)
      return value;
    else
      THROW(JParseError, "Text does not contain a number at given position " << buffer << ' ' << index);
  }


  /**
   * Make histogram axis logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  axis           axis
   */
  inline void setLogarithmic(TAxis* axis)
  {
    if (axis != NULL) {
      
      const Int_t N = axis->GetNbins();
      Double_t    buffer[N+1];
      
      buffer[0] = TMath::Power(10.0, axis->GetBinLowEdge(1));
      
      for (Int_t i = 1; i <= N; ++i) {
	buffer[i] = TMath::Power(10.0, axis->GetBinLowEdge(i) + axis->GetBinWidth(i));
      }
      
      axis->Set(N, buffer);
      /*
      if (axis->TestBit(TAxis::kAxisRange)) {

	const Int_t    first = axis->GetFirst();
	const Int_t    last  = axis->GetLast();

	if (first != last) {

	  const Double_t xmin  = axis->GetBinLowEdge(first);
	  const Double_t xmax  = axis->GetBinLowEdge(last)  +  axis->GetBinWidth(last);

	  axis->SetRangeUser(TMath::Power(10.0, xmin), TMath::Power(10.0, xmax));
	}
      }
      */
    }
  }


  /**
   * Make given parameter in formula logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  formula        formula
   * \param  parameter      parameter
   * \return                formula
   */
  inline TString getLogarithmic(const TString& formula, const char parameter)
  {
    const TRegexp regexp[] = {
      TString("^")          + TString(parameter) + TString("[^a-zA-Z_]"),   // parameter at start  of line
      TString("[^a-zA-Z_]") + TString(parameter) + TString("[^a-zA-Z_]"),   // parameter in middle of line
      TString("[^a-zA-Z_]") + TString(parameter) + TString("$")             // parameter at end    of line
    };

    const TString replacement = TString("log10(") + TString(parameter) + TString(")");

    TString buffer(formula);

    if (buffer.Length() == 1 && buffer[0] == parameter) {

      buffer = replacement;

    } else {
 
      for (Ssiz_t pos = 0, i; pos < buffer.Length(); pos += replacement.Length()) {
	if      ((i = buffer.Index(regexp[0], pos)) != -1) { buffer.Replace((pos = i + 0), 1, replacement); }
	else if ((i = buffer.Index(regexp[1], pos)) != -1) { buffer.Replace((pos = i + 1), 1, replacement); }
	else if ((i = buffer.Index(regexp[2], pos)) != -1) { buffer.Replace((pos = i + 1), 1, replacement); }
	else { break; }
      }
    }

    return buffer;
  }


  /**
   * Copy function parameters.
   *
   * \param  from           function
   * \param  to             function
   */
  inline void copy(const TF1& from, TF1& to)
  {
    static_cast<TAttLine&>  (to) = static_cast<const TAttLine&>  (from);
    static_cast<TAttFill&>  (to) = static_cast<const TAttFill&>  (from);
    static_cast<TAttMarker&>(to) = static_cast<const TAttMarker&>(from);

    to.SetParameters(from.GetParameters());

    to.SetNpx(from.GetNpx());
  }


  /**
   * Copy function parameters.
   *
   * \param  from           function
   * \param  to             function
   */
  inline void copy(const TF2& from, TF2& to)
  {
    copy(static_cast<const TF1&>(from), static_cast<TF1&>(to));

    to.SetNpy(from.GetNpy());
  }


  /**
   * Make x-axis of objects in list logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  list           list
   */
  template<class T>
  inline void setLogarithmicX(TList* list);


  /**
   * Make y-axis of objects in list logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  list           list
   */
  template<class T>
  inline void setLogarithmicY(TList* list);


  /**
   * Make x-axis of given function logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  f1             function
   */
  inline void setLogarithmicX(TF1* f1)
  {
    if (f1 != NULL) {

      TF1 fn(f1->GetName(), getLogarithmic(f1->GetExpFormula(), 'x'));

      copy(*f1, fn);

      fn.SetRange(f1->GetXmin(),
		  f1->GetXmax());

      *f1 = fn;

      f1->SetRange(TMath::Power(10.0,f1->GetXmin()),
		   TMath::Power(10.0,f1->GetXmax()));
    }
  }


  /**
   * Make y-axis of given function logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  f1             function
   */
  inline void setLogarithmicY(TF1* f1)
  {
    if (f1 != NULL) {

      TString buffer = f1->GetExpFormula();

      buffer = "pow(10.0, " + buffer + ")";

      TF1 fn(f1->GetName(), buffer);

      copy(*f1, fn);

      *f1 = fn;
    }
  }


  /**
   * Make x-axis of given function logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  f2             function
   */
  inline void setLogarithmicX(TF2* f2)
  {
    if (f2 != NULL) {

      TF2 fn(f2->GetName(), getLogarithmic(f2->GetExpFormula(), 'x'));

      copy(*f2, fn);

      fn.SetRange(f2->GetXmin(),
		  f2->GetYmin(),
		  f2->GetXmax(),
		  f2->GetYmax());

      *f2 = fn;

      f2->SetRange(TMath::Power(10.0,f2->GetXmin()),
		   f2->GetYmin(),
		   TMath::Power(10.0,f2->GetXmax()),
		   f2->GetYmax());
    }
  }


  /**
   * Make y-axis of given function logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  f2             function
   */
  inline void setLogarithmicY(TF2* f2)
  {
    if (f2 != NULL) {

      TF2 fn(f2->GetName(), getLogarithmic(f2->GetExpFormula(), 'y'));

      copy(*f2, fn);

      fn.SetRange(f2->GetXmin(),
		  f2->GetYmin(),
		  f2->GetXmax(),
		  f2->GetYmax());

      *f2 = fn;

      f2->SetRange(f2->GetXmin(),
		   TMath::Power(10.0,f2->GetYmin()),
		   f2->GetXmax(),
		   TMath::Power(10.0,f2->GetYmax()));
    }
  }


  /**
   * Make x-axis and associated functions of given histogram logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  h1             histogram
   */
  inline void setLogarithmicX(TH1* h1)
  {
    if (h1 != NULL) {

      setLogarithmicX<TF1>(h1->GetListOfFunctions());
      
      setLogarithmic(h1->GetXaxis());
    }
  }


  /**
   * Make x-axis and associated functions of given histogram logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  h2             histogram
   */
  inline void setLogarithmicX(TH2* h2)
  {
    using namespace std;

    if (h2 != NULL) {

      setLogarithmicX<TF2>(h2->GetListOfFunctions());

      setLogarithmic(h2->GetXaxis());
    }
  }


  /**
   * Make y-axis and associated functions of given histogram logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  h2             histogram
   */
  inline void setLogarithmicY(TH2* h2)
  {
    if (h2 != NULL) {

      setLogarithmicY<TF2>(h2->GetListOfFunctions());

      setLogarithmic(h2->GetYaxis());
    }
  }


  /**
   * Make x-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g1             graph
   */
  inline void setLogarithmicX(TGraph* g1)
  {
    if (g1 != NULL) {

      setLogarithmicX<TF1>(g1->GetListOfFunctions());

      for (Int_t i = 0; i != g1->GetN(); ++i) {
	g1->GetX()[i] = pow(10.0, g1->GetX()[i]);
      }
    }
  }


  /**
   * Make y-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g1             graph
   */
  inline void setLogarithmicY(TGraph* g1)
  {
    if (g1 != NULL) {

      setLogarithmicY<TF1>(g1->GetListOfFunctions());

      for (Int_t i = 0; i != g1->GetN(); ++i) {
	g1->GetY()[i] = pow(10.0, g1->GetY()[i]);
      }
    }
  }


  /**
   * Make x-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g1             graph
   */
  inline void setLogarithmicX(TGraphErrors* g1)
  {
    if (g1 != NULL) {

      setLogarithmicX<TF1>(g1->GetListOfFunctions());

      for (Int_t i = 0; i != g1->GetN(); ++i) {
	g1->GetEX()[i] = pow(10.0, g1->GetX()[i] + g1->GetEX()[i])  -  pow(10.0, g1->GetX()[i]);
	g1->GetX() [i] = pow(10.0, g1->GetX()[i]);
      }
    }
  }


  /**
   * Make y-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g1             graph
   */
  inline void setLogarithmicY(TGraphErrors* g1)
  {
    if (g1 != NULL) {

      setLogarithmicY<TF1>(g1->GetListOfFunctions());

      for (Int_t i = 0; i != g1->GetN(); ++i) {
	g1->GetEY()[i] = pow(10.0, g1->GetY()[i] + g1->GetEY()[i])  -  pow(10.0, g1->GetY()[i]);
	g1->GetY() [i] = pow(10.0, g1->GetY()[i]);
      }
    }
  }


  /**
   * Make x-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g2             graph
   */
  inline void setLogarithmicX(TGraph2D* g2)
  {
    if (g2 != NULL) {

      setLogarithmicX<TF2>(g2->GetListOfFunctions());

      for (Int_t i = 0; i != g2->GetN(); ++i) {
	g2->GetX()[i] = pow(10.0, g2->GetX()[i]);
      }
    }
  }


  /**
   * Make y-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g2             graph
   */
  inline void setLogarithmicY(TGraph2D* g2)
  {
    if (g2 != NULL) {

      setLogarithmicY<TF2>(g2->GetListOfFunctions());

      for (Int_t i = 0; i != g2->GetN(); ++i) {
	g2->GetY()[i] = pow(10.0, g2->GetY()[i]);
      }
    }
  }


  /**
   * Make x-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g2             graph
   */
  inline void setLogarithmicX(TGraph2DErrors* g2)
  {
    if (g2 != NULL) {

      setLogarithmicX<TF2>(g2->GetListOfFunctions());

      for (Int_t i = 0; i != g2->GetN(); ++i) {
	g2->GetEX()[i] = pow(10.0, g2->GetX()[i] + g2->GetEX()[i])  -  pow(10.0, g2->GetX()[i]);
	g2->GetX() [i] = pow(10.0, g2->GetX()[i]);
      }
    }
  }


  /**
   * Make y-axis and associated functions of given graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  g2             graph
   */
  inline void setLogarithmicY(TGraph2DErrors* g2)
  {
    if (g2 != NULL) {

      setLogarithmicY<TF2>(g2->GetListOfFunctions());

      for (Int_t i = 0; i != g2->GetN(); ++i) {
	g2->GetEY()[i] = pow(10.0, g2->GetY()[i] + g2->GetEY()[i])  -  pow(10.0, g2->GetY()[i]);
	g2->GetY() [i] = pow(10.0, g2->GetY()[i]);
      }
    }
  }


  /**
   * Make x-axis of given multi-graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  gn             multi graph
   */
  inline void setLogarithmicX(TMultiGraph* gn)
  {
    if (gn != NULL) {
      setLogarithmicX<TGraph>(gn->GetListOfGraphs());
    }
  }


  /**
   * Make y-axis of given multi-graph logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  gn             multi graph
   */
  inline void setLogarithmicY(TMultiGraph* gn)
  {
    if (gn != NULL) {
      setLogarithmicY<TGraph>(gn->GetListOfGraphs());
    }
  }


  /**
   * Make x-axis of given line logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  line           line
   */
  inline void setLogarithmicX(TLine* line)
  {
    if (line != NULL) {
      line->SetX1(pow(10.0, line->GetX1()));
      line->SetX2(pow(10.0, line->GetX2()));
    }
  }


  /**
   * Make y-axis of given line logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  line           line
   */
  inline void setLogarithmicY(TLine* line)
  {
    if (line != NULL) {
      line->SetY1(pow(10.0, line->GetY1()));
      line->SetY2(pow(10.0, line->GetY2()));
    }
  }


  /**
   * Make x-axis of given ellipse logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  ellipse        ellipse
   */
  inline void setLogarithmicX(TEllipse* ellipse)
  {
    THROW(JFunctionalException, "Operation setLogarithmicX on TEllipse not allowed.");
  }


  /**
   * Make y-axis of given ellipse logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  ellipse        ellipse
   */
  inline void setLogarithmicY(TEllipse* ellipse)
  {
    THROW(JFunctionalException, "Operation setLogarithmicY on TEllipse not allowed.");
  }


  /**
   * Make x-axis of objects in list logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  list           list
   */
  template<class T>
  inline void setLogarithmicX(TList* list)
  {
    for (TIter i(list); T* p = dynamic_cast<T*>(i.Next()); ) {
      setLogarithmicX(p);
    }
  }


  /**
   * Make y-axis of objects in list logarithmic (e.g.\ after using <tt>log10()</tt>).
   *
   * \param  list           list
   */
  template<class T>
  inline void setLogarithmicY(TList* list)
  {
    for (TIter i(list); T* p = dynamic_cast<T*>(i.Next()); ) {
      setLogarithmicY(p);
    }
  }

    
  /**
   * Convert 1D histogram to PDF.
   *
   * Possible options are:
   *    - N  normalise to histogram contents;
   *    - W  divide by bin width;
   *    - E  convert also bin errors.
   *
   * \param  h1             histogram
   * \param  option         option
   * \param  factor         scaling factor
   */
  inline void convertToPDF(TH1& h1, const std::string& option = "NW", const double factor = 1.0)
  {
    using namespace std;
    
    const bool normalise = (option.find('N') != string::npos || option.find('n') != string::npos);
    const bool bin_width = (option.find('W') != string::npos || option.find('w') != string::npos);
    const bool use_error = (option.find('E') != string::npos || option.find('e') != string::npos);

    Double_t W = 1.0;

    if (normalise) {

      W = 0.0;
      
      for (Int_t i = 1; i <= h1.GetXaxis()->GetNbins(); ++i) {
	W += h1.GetBinContent(i);
      }
    }
    
    if (W != 0.0) {
      
      for (Int_t i = 1; i <= h1.GetXaxis()->GetNbins(); ++i) {

	const Double_t w = W * (bin_width ? h1.GetXaxis()->GetBinWidth(i) : 1.0);
	
	h1.SetBinContent(i, h1.GetBinContent(i) * factor / w);
	
	if (use_error) {
	  h1.SetBinError(i, h1.GetBinError(i) * factor / w);
	}
      }
    }
  }
  
    
  /**
   * Convert 2D histogram to PDF.
   *
   * Possible options are:
   *    - N  normalise to histogram contents;
   *    - X  convert x-axis to PDF;
   *    - Y  convert y-axis to PDF;
   *    - W  divide by bin width;
   *    - E  convert also bin errors.
   *
   * \param  h2             histogram
   * \param  option         option
   * \param  factor         scaling factor
   */
  inline void convertToPDF(TH2& h2, const std::string& option = "NXYW", const double factor = 1.0)
  {
    using namespace std;
    
    const bool normalise = (option.find('N') != string::npos || option.find('n') != string::npos);
    const bool X         = (option.find('X') != string::npos || option.find('x') != string::npos);
    const bool Y         = (option.find('Y') != string::npos || option.find('y') != string::npos);
    const bool bin_width = (option.find('W') != string::npos || option.find('w') != string::npos);
    const bool use_error = (option.find('E') != string::npos || option.find('e') != string::npos);

    Double_t W = 1.0;
    
    if (X && Y) {

      if (normalise) {
      
	W = 0.0;
	
	for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
	  for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
	    W += h2.GetBinContent(ix,iy);
	  }
	}
      }

      if (W != 0.0) {

	for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
	  for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
	    
	    const Double_t w = W * (bin_width ? h2.GetXaxis()->GetBinWidth(ix) * h2.GetYaxis()->GetBinWidth(iy) : 1.0);
	      
	    h2.SetBinContent(ix, iy, h2.GetBinContent(ix,iy) * factor / w);
	  
	    if (use_error) {
	      h2.SetBinError(ix, iy, h2.GetBinError(ix,iy) * factor / w);
	    }
	  }
	}
      }
      
    } else if (X) {
    
      for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {

	if (normalise) {
	  
	  W = 0.0;
	
	  for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
	    W += h2.GetBinContent(ix,iy);
	  }
	}

	if (W != 0.0) {
	  
	  for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
	    
	    const Double_t w = W * (bin_width ? h2.GetXaxis()->GetBinWidth(ix) : 1.0);

	    h2.SetBinContent(ix, iy, h2.GetBinContent(ix,iy) * factor / w);
	  
	    if (use_error) {
	      h2.SetBinError(ix, iy, h2.GetBinError(ix,iy) * factor / w);
	    }
	  }
	}
      }
      
    } else if (Y) {
    
      for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {

	if (normalise) {

	  W = 0.0;

	  for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
	    W += h2.GetBinContent(ix,iy);
	  }
	}

	if (W != 0.0) {
	  
	  for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
	    
	    const Double_t w = W * (bin_width ? h2.GetYaxis()->GetBinWidth(iy) : 1.0);

	    h2.SetBinContent(ix, iy, h2.GetBinContent(ix,iy) / w);
	  
	    if (use_error) {
	      h2.SetBinError(ix, iy, h2.GetBinError(ix,iy) / w);
	    }
	  }
	}
      }
    }
  }


  /**
   * Convert 3D histogram to PDF.
   *
   * Possible options are:
   *    - N  normalise to histogram contents;
   *    - X  convert x-axis to PDF;
   *    - Y  convert y-axis to PDF;
   *    - Z  convert z-axis to PDF;
   *    - W  divide by bin width;
   *    - E  convert also bin errors.
   *
   * \param  h3             histogram
   * \param  option         option
   * \param  factor         scaling factor
   */
  inline void convertToPDF(TH3& h3, const std::string& option = "NXYW", const double factor = 1.0)
  {
    using namespace std;
    
    const bool normalise = (option.find('N') != string::npos || option.find('n') != string::npos);
    const bool X         = (option.find('X') != string::npos || option.find('x') != string::npos);
    const bool Y         = (option.find('Y') != string::npos || option.find('y') != string::npos);
    const bool Z         = (option.find('Z') != string::npos || option.find('z') != string::npos);    
    const bool bin_width = (option.find('W') != string::npos || option.find('w') != string::npos);
    const bool use_error = (option.find('E') != string::npos || option.find('e') != string::npos);

    Double_t W = 1.0;

    if (X && Y && Z) {

      if (normalise) {
      
	W = 0.0;
	
	for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	  for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {	    
	      W += h3.GetBinContent(ix,iy,iz);
	    }
	  }
	}
      }

      if (W != 0.0) {

	for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	  for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {	    
	    
	      const Double_t w = W * (bin_width ? h3.GetXaxis()->GetBinWidth(ix) * h3.GetYaxis()->GetBinWidth(iy) * h3.GetZaxis()->GetBinWidth(iz) : 1.0);
	      
	      h3.SetBinContent(ix, iy, iz, h3.GetBinContent(ix,iy,iz) * factor / w);
	  
	      if (use_error) {
		h3.SetBinError(ix, iy, iz, h3.GetBinError(ix,iy,iz) * factor / w);
	      }
	    }
	  }
	}
      }
      
    } else if (X && Z) {

      for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	
	if (normalise) {
      
	  W = 0.0;
	
	  for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
	      W += h3.GetBinContent(ix,iy,iz);
	    }
	  }
	}

	if (W != 0.0) {

	  for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
	    
	      const Double_t w = W * (bin_width ? h3.GetXaxis()->GetBinWidth(ix) * h3.GetZaxis()->GetBinWidth(iz) : 1.0);
	      
	      h3.SetBinContent(ix, iy, iz, h3.GetBinContent(ix,iy,iz) * factor / w);
	  
	      if (use_error) {
		h3.SetBinError(ix, iy, iz, h3.GetBinError(ix,iy,iz) * factor / w);
	      }
	    }
	  }
	}
      }
      
    } else if (Y && Z) {

      for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	
	if (normalise) {
      
	  W = 0.0;
	
	  for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
	      W += h3.GetBinContent(ix,iy,iz);
	    }
	  }
	}

	if (W != 0.0) {

	  for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
	    
	      const Double_t w = W * (bin_width ? h3.GetYaxis()->GetBinWidth(iy) * h3.GetZaxis()->GetBinWidth(iz) : 1.0);
	      
	      h3.SetBinContent(ix, iy, iz, h3.GetBinContent(ix,iy,iz) * factor / w);
	  
	      if (use_error) {
		h3.SetBinError(ix, iy, iz, h3.GetBinError(ix,iy,iz) * factor / w);
	      }
	    }
	  }
	}
      }
      
    } else if (X && Y) {

      for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {      

	if (normalise) {
      
	  W = 0.0;
	
	  for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	    for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	      W += h3.GetBinContent(ix,iy,iz);
	    }
	  }
	}

	if (W != 0.0) {

	  for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	    for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	    
	      const Double_t w = W * (bin_width ? h3.GetXaxis()->GetBinWidth(ix) * h3.GetYaxis()->GetBinWidth(iy) : 1.0);
	      
	      h3.SetBinContent(ix, iy, iz, h3.GetBinContent(ix,iy,iz) * factor / w);
	  
	      if (use_error) {
		h3.SetBinError(ix, iy, iz, h3.GetBinError(ix,iy,iz) * factor / w);
	      }
	    }
	  }
	}
      }
      
    } else if (X) {
    
      for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {

	  if (normalise) {
	  
	    W = 0.0;
	
	    for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	      W += h3.GetBinContent(ix,iy,iz);
	    }
	  }

	  if (W != 0.0) {
	  
	    for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	    
	      const Double_t w = W * (bin_width ? h3.GetXaxis()->GetBinWidth(ix) : 1.0);

	      h3.SetBinContent(ix, iy, iz, h3.GetBinContent(ix,iy,iz) * factor / w);
	  
	      if (use_error) {
		h3.SetBinError(ix, iy, iz, h3.GetBinError(ix,iy,iz) * factor / w);
	      }
	    }
	  }
	}
      }
      
    } else if (Y) {
    
      for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {	

	  if (normalise) {

	    W = 0.0;

	    for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	      W += h3.GetBinContent(ix,iy,iz);
	    }
	  }

	  if (W != 0.0) {
	  
	    for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	    
	      const Double_t w = W * (bin_width ? h3.GetYaxis()->GetBinWidth(iy) : 1.0);

	      h3.SetBinContent(ix, iy, iz, h3.GetBinContent(ix,iy,iz) / w);
	  
	      if (use_error) {
		h3.SetBinError(ix, iy, iz, h3.GetBinError(ix,iy,iz) / w);
	      }
	    }
	  }
	}
      }
      
    } else if (Z) {

      for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
	for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {	

	  if (normalise) {

	    W = 0.0;

	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
	      W += h3.GetBinContent(ix,iy,iz);
	    }
	  }

	  if (W != 0.0) {
	  
	    for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {
	    
	      const Double_t w = W * (bin_width ? h3.GetZaxis()->GetBinWidth(iz) : 1.0);

	      h3.SetBinContent(ix, iy, iz, h3.GetBinContent(ix,iy,iz) / w);
	  
	      if (use_error) {
		h3.SetBinError(ix, iy, iz, h3.GetBinError(ix,iy,iz) / w);
	      }
	    }
	  }
	}
      }
    }
  }


  /**
   * Set limits of TGraph.
   *
   * \param  g1             graph
   */
  inline void setLimits(TGraph& g1)
  {
    using namespace std;
    
    Double_t ymin = numeric_limits<Double_t>::max();
    Double_t ymax = numeric_limits<Double_t>::lowest();

    for (Int_t i = 0; i != g1.GetN(); ++i) {

      const Double_t y = g1.GetY()[i];

      if (y > ymax) { ymax = y; }
      if (y < ymin) { ymin = y; }
    }

    g1.SetMinimum(ymin);
    g1.SetMaximum(ymax);
  }


  /**
   * Set limits of TGraph2D.
   *
   * \param  g2             graph
   */
  inline void setLimits(TGraph2D& g2)
  {
    using namespace std;
    
    Double_t zmin = numeric_limits<Double_t>::max();
    Double_t zmax = numeric_limits<Double_t>::lowest();

    for (Int_t i = 0; i != g2.GetN(); ++i) {

      const Double_t z = g2.GetZ()[i];

      if (z > zmax) { zmax = z; }
      if (z < zmin) { zmin = z; }
    }

    g2.SetMinimum(zmin);
    g2.SetMaximum(zmax);
  }


  /**
   * Set axis range.
   *
   * \param  xmin           lower limit (I/O)
   * \param  xmax           upper limit (I/O)
   * \param  logx           logarithmic
   */
  inline void setRange(double& xmin,
		       double& xmax,
		       const bool logx)
  {
    if (logx) {
      xmin = log(xmin);
      xmax = log(xmax);
    }

    double dx = (xmax - xmin) * 0.1;

    if (logx || xmin >= dx || xmin < 0.0)
      xmin -= dx;
    else
      xmin  = 0.0;

    xmax += dx;

    if (logx) {
      xmin = exp(xmin);
      xmax = exp(xmax);
    }
  }


  /**
   * Set axis with PMT address labels.
   *
   * This methods sets the labels of the given axis to the sorted values of the PMT ring and position.\n
   * It should normally be called before filling of the corresponding histogram.\n
   * The filling should then be made with the textual representation of the PMT ring and position (i.e.\ JDETECTOR::JPMTPhysicalAddress::toString).
   *
   * Alternatively, the filling can be made with the index of the PMT in the address map of the corresponding module 
   * (i.e.\ JDETECTOR::JModuleAddressMap::getIndex).\n
   * In that case, the method can be called before or after filling of the histogram.
   *
   * \param axis           axis
   * \param memo           module address map
   */
  inline void setAxisLabels(TAxis *axis, const JModuleAddressMap& memo)
  {
    using namespace JPP;

    if (axis->GetNbins() == (int) memo.size()) {

      for (int i = 0; i != axis->GetNbins(); ++i) {

	const JPMTPhysicalAddress& address = memo[i];

	axis->SetBinLabel(i + 1, address.toString().c_str());
      }

    } else {

      THROW(JValueOutOfRange, "Number of bins " << axis->GetNbins() << " != " << memo.size());
    }
  }


  /**
   * Set axis labels with PMT addresses.
   *
   * This methods sets the labels of the given axis to the sorted values of the PMT ring and position.\n
   * It should be called after filling of the corresponding histogram.\n
   * The filling should have been made with the PMT number (e.g.\ KM3NETDAQ::JDAQHit::getPMT).
   *
   * \param h1             histogram
   * \param axis           axis <tt>(x, X, y, Y, z, Z)</tt>
   * \param memo           module address map
   */
  inline void setAxisLabels(TH1& h1, const std::string axis, const JModuleAddressMap& memo)
  {
    using namespace JPP;

    TAxis* pax = NULL;

    if      (axis == "x" || axis == "X") 
      pax = h1.GetXaxis();
    else if (axis == "y" || axis == "Y") 
      pax = h1.GetYaxis();
    else if (axis == "z" || axis == "Z") 
      pax = h1.GetZaxis();
    else
      THROW(JValueOutOfRange, "Invalid axis " << axis);

    if (pax->GetNbins() == (int) memo.size()) {

      for (int i = 0; i != pax->GetNbins(); ++i) {

	const JPMTPhysicalAddress& address = memo.getAddressTranslator(i);

	pax->SetBinLabel(i + 1, address.toString().c_str());
      }

      h1.LabelsOption("a", axis.c_str());       // sort labels

    } else {

      THROW(JValueOutOfRange, "Number of bins " << pax->GetNbins() << " != " << memo.size());
    }
  } 


  /**
   * Check if given key corresponds to a TObject.
   *
   * \param key            ROOT key
   * \return               true if given key corresponds to a TObject; else false
   */
  inline bool isTObject(const TKey* key)
  {
    return (key != NULL && TClass::GetClass(key->GetClassName())->IsTObject());
  }


  /**
   * Get first object of which name matches given reguar expression.
   *
   * \param  ls            pointer to list of objects
   * \param  name          regular expression
   * \return               pointer to object (maybe NULL)
   */
  inline TObject* getObject(TList* ls, const char* const name)
  {
    const TRegexp regexp(name);

    TIter iter(ls);

    for (TObject* p1; (p1 = (TObject*) iter.Next()) != NULL; ) {

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

      if (tag.Contains(regexp)) {
	return p1;
      }
    }

    return NULL;
  }


  /**
   * Get function.
   *
   * \param  h1            histogram
   * \param  fcn           function name
   * \return               pointer to function (maybe NULL)
   */
  inline TF1* getFunction(TH1* h1, const char* const fcn)
  {
    return (TF1*) getObject(h1->GetListOfFunctions(), fcn);
  }


  /**
   * Get function.
   *
   * \param  g1            graph
   * \param  fcn           function name
   * \return               pointer to function (maybe NULL)
   */
  inline TF1* getFunction(TGraph* g1, const char* const fcn)
  {
    return (TF1*) getObject(g1->GetListOfFunctions(), fcn);
  }


  /**
   * Get function.
   *
   * \param  g2            graph
   * \param  fcn           function name
   * \return               pointer to function (maybe NULL)
   */
  inline TF1* getFunction(TGraph2D* g2, const char* const fcn)
  {
    return (TF1*) getObject(g2->GetListOfFunctions(), fcn);
  }
}

#endif