JGizmoToolkit.hh

Go to the documentation of this file.

#ifndef __JGIZMOTOOLKIT__
#define __JGIZMOTOOLKIT__

#include <string>
#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 "TString.h"
#include "TRegexp.h"
#include "TFormula.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
  };


  /**
   * 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.\n
   * The TFile pointer of an already opened file is recovered, else a new file is opened.
   *
   * \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.\n
   * 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 TObject.
   *
   * \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 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.\n
   * The formula may contain names of member methods of the object pointed to.
   * These methods should have no arguments and the return type <tt>Double_t</tt>.
   * 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) {
        
        TIterator* iter = p->GetListOfAllPublicMethods()->MakeIterator();
        
        for (TMethod* method; (method = (TMethod*) iter->Next()) != NULL; ) {

          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) {

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

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

            } else if (pos == index + (Ssiz_t) strlen(method->GetName())) {

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

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

            } else {

              continue;
            }

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

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


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


  /**
   * Make axis logarithmic (e.g. after filling with log10()).
   *
   * \param  axis           axis
   */
  inline void setLogarithm(TAxis* axis)
  {
    if (axis != NULL) {

      const int N = axis->GetNbins();
      Double_t  buffer[N+1];
      
      buffer[0] = TMath::Power(10.0, axis->GetBinLowEdge(1));
      
      for (int i = 1; i <= N; ++i) {
        buffer[i] = TMath::Power(10.0, axis->GetBinLowEdge(i) + axis->GetBinWidth(i));
      }
      
      axis->Set(N, buffer);
    }
  }

    
  /**
   * 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 i = 1; i <= h2.GetXaxis()->GetNbins(); ++i) {
          for (Int_t j = 1; j <= h2.GetYaxis()->GetNbins(); ++j) {
            W += h2.GetBinContent(i,j);
          }
        }
      }

      if (W != 0.0) {

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

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

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

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

        if (normalise) {

          W = 0.0;

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

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

            h2.SetBinContent(i, j, h2.GetBinContent(i,j) / w);
          
            if (use_error) {
              h2.SetBinError(i, j, h2.GetBinError(i,j) / 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>::max();

    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>::max();

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


  /**
   * initialize axis with PMT address labels
   *
   * \param axis           axis
   * \param memo           default 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& pmtAddress = memo[i];

        axis->SetBinLabel(i + 1, TString(pmtAddress.toString()));
      }

    } else {

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

#endif