JTestChi2.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TProfile.h"

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

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


/**
 * \file
 * Auxiliary program for chi2 test of 1D histograms.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  
  JRootObjectID  inputFileA;
  JRootObjectID  inputFileB;
  string         option;
  double         P;
  int            debug;

  try {

    JParser<> zap("Auxiliary program for chi2 test of 1D histograms.");

    zap['a'] = make_field(inputFileA,  "measurement histogram, e.g: <file name>:<object name>");
    zap['b'] = make_field(inputFileA,  "reference   histogram, e.g: <file name>:<object name>");
    zap['O'] = make_field(option)     = 
      "UU",     "UW",     "WW",
      "UUUF",   "UWUF",   "WWUF",     // with underflows
      "UUOF",   "UWOF",   "WWOF",     // with overflows
      "UUUFOF", "UWUFOF", "WWUFOF";   // with underflows and overflows
    zap['P'] = make_field(P,            "minimal p-value") = 0.0;
    zap['d'] = make_field(debug)                           = 1;

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


  TObject* ha = getObject(inputFileA);
  TObject* hb = getObject(inputFileA);

  if (ha == NULL) { FATAL("No object at " << inputFileA << endl); };
  if (hb == NULL) { FATAL("No object at " << inputFileB << endl); };

  TH1* h1a = NULL;
  TH1* h1b = NULL;

  if (h1a == NULL && dynamic_cast<TProfile*>(ha) != NULL) { h1a = dynamic_cast<TProfile*>(ha)->ProjectionX(); }
  if (h1b == NULL && dynamic_cast<TProfile*>(hb) != NULL) { h1b = dynamic_cast<TProfile*>(hb)->ProjectionX(); }

  if (h1a == NULL && dynamic_cast<TH1*>(ha)      != NULL) { h1a = dynamic_cast<TH1*>(ha); }
  if (h1b == NULL && dynamic_cast<TH1*>(hb)      != NULL) { h1b = dynamic_cast<TH1*>(hb); }

  if (h1a == NULL) { FATAL("Object at " << inputFileA << " is not TH1." << endl); };
  if (h1b == NULL) { FATAL("Object at " << inputFileB << " is not TH1." << endl); };

  const Double_t result = h1a->Chi2Test(h1b, option.c_str());

  NOTICE("Result " << result << endl);

  if (P > 0.0) {
    ASSERT(result >= P);
  }

  return 0;
}