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