JVariance2D.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH1D.h"
#include "TH2.h"
#include "TString.h"
#include "TRegexp.h"

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

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

namespace {
  /**
   * Additional options for formula.
   */
  static const char* const first_t = "first";
}


/**
 * \file
 * Auxiliary program to histogram bin-by-bin deviations of the contents of one (or more) 1D histogram(s).
 * The option <tt>-f</tt> corresponds to <tt><file name>:<object name></tt>.
 *
 * The formula (option <tt>-F <formula></tt>) refers to a ROOT TFormula.
 * The expression may contain member methods of the corresponding object.
 *
 * In case of multiple histograms, the formula is evaluated stand alone (i.e. without ROOT histogram).
 * If formula equals "first", the bin-by-bin deviations of the contents of the histograms are
 * evaluated with respect to the first histogram.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JAbstractHistogram<Double_t> JHistogram_t;
  
  vector<JRootObjectID> inputFile;
  string                outputFile;
  bool                  reuse;
  JHistogram_t          X;
  TString               formula;
  bool                  reverse;
  int                   debug;

  try {

    JParser<> zap("Auxiliary program to histogram bin-by-bin deviations of a set of 2D histograms.");

    zap['f'] = make_field(inputFile,  "<input file>:<object name>");
    zap['o'] = make_field(outputFile, "ROOT file with histogram (possibly I/O, option -r)")       = "variance.root";
    zap['r'] = make_field(reuse,      "reuse histogram from existing output file");
    zap['x'] = make_field(X,          "histogram binning");
    zap['F'] = make_field(formula,    "ROOT TFormula (may contain method names of object)")       = "0";
    zap['R'] = make_field(reverse,    "reverse sign");
    zap['d'] = make_field(debug)          = 1;

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


  vector<TH2*> listOfHistograms;

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

    TH2* h2 = dynamic_cast<TH2*>(getObject(*input));

    if (h2 != NULL) {
      listOfHistograms.push_back(h2);
    }
  }

  if (listOfHistograms.empty()) {
    FATAL("No histograms." << endl);
  }

  if (formula == "") {
    FATAL("Empty formula." << endl);
  }

  TH1* h0 = NULL;

  if (reuse) {
    h0 = dynamic_cast<TH1*>(getObject(JRootObjectID(outputFile, "h0")));
  }

  if (h0 == NULL) {
    h0 = new TH1D("h0", NULL,
                  X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit());
  }


  TH2* h2 = listOfHistograms[0];

  if (listOfHistograms.size() == 1) {

    const double y1 = getResult(formula, h2);

    DEBUG(h2->GetName() << ' ' << formula << ' ' << y1 << endl);

    for (Int_t ix = 1; ix <= h2->GetNbinsX(); ++ix) {
      for (Int_t iy = 1; iy <= h2->GetNbinsY(); ++iy) {

        const Double_t yp = h2->GetBinContent(ix,iy); 

        h0->Fill(reverse ? y1 - yp : yp - y1);
      }
    }

  } else {

    if (formula == first_t) {

      for (Int_t ix = 1; ix <= h2->GetNbinsX(); ++ix) {
        for (Int_t iy = 1; iy <= h2->GetNbinsY(); ++iy) {

          const double y1 = h2->GetBinContent(ix,iy);
        
          for (vector<TH2*>::const_iterator p = listOfHistograms.begin(); ++p != listOfHistograms.end(); ) {
          
            const Double_t yp = (*p)->GetBinContent(ix,iy);
          
            h0->Fill(reverse ? y1 - yp : yp - y1);
          }
        }
      }

    } else {

      const double y1 = getResult(formula);

      for (Int_t ix = 1; ix <= h2->GetNbinsX(); ++ix) {
        for (Int_t iy = 1; iy <= h2->GetNbinsY(); ++iy) {

          for (vector<TH2*>::const_iterator p = listOfHistograms.begin(); p != listOfHistograms.end(); ++p) {
          
            const Double_t yp = (*p)->GetBinContent(ix,iy);

            h0->Fill(reverse ? y1 - yp : yp - y1);
          }
        }
      }
    }
  }


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