JGizmoToolkit.hh

Go to the documentation of this file.

#ifndef __JGIZMOTOOLKIT__
#define __JGIZMOTOOLKIT__

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

#include "TError.h"
#include "TFile.h"
#include "TClass.h"
#include "TObject.h"
#include "TKey.h"
#include "TH1.h"
#include "TH2.h"
#include "TGraph.h"
#include "TGraph2D.h"
#include "TMultiGraph.h"
#include "TString.h"
#include "TRegexp.h"
#include "TFormula.h"
#include "TF1.h"
#include "TF2.h"
#include "TIterator.h"
#include "TMethod.h"
#include "TMethodCall.h"
#include "TAxis.h"
#include "TMath.h"

#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 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 TH1.
   *
   * \param  object         pointer to TObject
   * \return                true if TH1 looks like a line; else false
   */
  inline bool isTAttLine(const TObject* object)
  {
    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;

    } else {

      return false;
    }
  }


  /**
   * Get result of given textual formula.
   *
   * The formula may contain names of member methods of the object pointed to.\n
   * These methods could have arguments and the return type should be compatible with <tt>Double_t</tt>.\n
   * 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 ( ; ; ) {

          TIterator* iter   = p->GetListOfAllPublicMethods()->MakeIterator();
          TMethod*   method = NULL;

          for (TMethod* p; (p = (TMethod*) iter->Next()) != NULL; ) {
            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("%15.5e", 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 could have arguments and the return type should be compatible with <tt>Double_t</tt>.\n
   * 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 could have arguments and the return type should be compatible with <tt>Double_t</tt>.\n
   * 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 could have arguments and the return type should be compatible with <tt>Double_t</tt>.\n
   * 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 des not contain a number at given position " << buffer << ' ' << index);
  }


  /**
   * Make histogram axis logarithmic (e.g.\ after filling with log10()).
   *
   * \param  axis           axis
   */
  inline void setLogarithmic(TAxis* axis)
  {
    if (axis != NULL) {
      
      const Int_t    first = axis->GetFirst();
      const Int_t    last  = axis->GetLast();

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

      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)) {
        axis->SetRangeUser(TMath::Power(10.0, xmin), TMath::Power(10.0, xmax));
      }
    }
  }


  /**
   * Make given parameter in formula logarithmic (e.g.\ after filling with log10()).
   *
   * \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 parameter <tt>x</tt> of function logarithmic (e.g.\ after filling with log10()).
   *
   * \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 parameter <tt>x</tt> of function logarithmic (e.g.\ after filling with log10()).
   *
   * \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 parameter <tt>y</tt> of function logarithmic (e.g.\ after filling with log10()).
   *
   * \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 filling with log10()).
   *
   * \param  h1             histogram
   */
  inline void setLogarithmicX(TH1* h1)
  {
    if (h1 != NULL) {

      for (TIter iter(h1->GetListOfFunctions()); TF1* f1 = dynamic_cast<TF1*>(iter.Next()); ) {
        setLogarithmicX(f1);
      }
      
      setLogarithmic(h1->GetXaxis());
    }
  }


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

    if (h2 != NULL) {

      for (TIter iter(h2->GetListOfFunctions()); TF2* f2 = dynamic_cast<TF2*>(iter.Next()); ) {
        setLogarithmicX(f2);
      }

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


  /**
   * Make Y-axis and associated functions of given histogram logarithmic (e.g.\ after filling with log10()).
   *
   * \param  h2             histogram
   */
  inline void setLogarithmicY(TH2* h2)
  {
    if (h2 != NULL) {

      for (TIter iter(h2->GetListOfFunctions()); TF2* f2 = dynamic_cast<TF2*>(iter.Next()); ) {
        setLogarithmicY(f2);
      }

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


  /**
   * Make parameter <tt>x</tt> of graph logarithmic (e.g.\ after filling with log10()).
   *
   * \param  g1             graph
   */
  inline void setLogarithmicX(TGraph* g1)
  {
    if (g1 != NULL) {

      for (TIter iter(g1->GetListOfFunctions()); TF1* f1 = dynamic_cast<TF1*>(iter.Next()); ) {
        setLogarithmicX(f1);
      }

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


  /**
   * Make parameter <tt>x</tt> of graph logarithmic (e.g.\ after filling with log10()).
   *
   * \param  g2             graph
   */
  inline void setLogarithmicX(TGraph2D* g2)
  {
    if (g2 != NULL) {

      for (TIter iter(g2->GetListOfFunctions()); TF2* f2 = dynamic_cast<TF2*>(iter.Next()); ) {
        setLogarithmicX(f2);
      }

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


  /**
   * Make parameter <tt>y</tt> of graph logarithmic (e.g.\ after filling with log10()).
   *
   * \param  g2             graph
   */
  inline void setLogarithmicY(TGraph2D* g2)
  {
    if (g2 != NULL) {

      for (TIter iter(g2->GetListOfFunctions()); TF2* f2 = dynamic_cast<TF2*>(iter.Next()); ) {
        setLogarithmicY(f2);
      }

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


  /**
   * Make parameter <tt>x</tt> of multi graph logarithmic (e.g.\ after filling with log10()).
   *
   * \param  m1             multi graph
   */
  inline void setLogarithmicX(TMultiGraph* m1)
  {
    if (m1 != NULL) {

      for (TIter i1(m1->GetListOfGraphs()); TGraph* g1 = dynamic_cast<TGraph*>(i1()); ) {
        setLogarithmicX(g1);
      }
    }
  }

    
  /**
   * 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);
            }
          }
        }
      }
    }
  }


  /**
   * 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 (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());
    }
  } 
}

#endif