JGizmo/JHistogram2D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <vector>

#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"
#include "TF2.h"

#include "JLang/JToken.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JGizmo/JGizmoToolkit.hh"

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


/**
 * \file
 *
 * Program to create TH1D and fill according given formula.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JToken<';'>                  JToken_t;
  typedef JAbstractHistogram<Double_t> JHistogram_t;

  string            outputFile;
  string            inputFile;
  string            formula;
  vector<JToken_t>  parameters;
  Int_t             numberOfEvents;
  JHistogram_t      X;
  JHistogram_t      Y;
  int               debug;

  try {

    JParser<> zap("Program to create TH1D and fill according given formula.");
    
    zap['o'] = make_field(outputFile);
    zap['f'] = make_field(inputFile)       = "";
    zap['F'] = make_field(formula)         = "";
    zap['@'] = make_field(parameters)      = JPARSER::initialised();
    zap['n'] = make_field(numberOfEvents)  = 0;
    zap['x'] = make_field(X)               = JHistogram_t(100, -1.0, +1.0);
    zap['y'] = make_field(Y)               = JHistogram_t(100, -1.0, +1.0);
    zap['d'] = make_field(debug)           = 1;

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


  if ((formula != "" && inputFile != "") ||
      (formula == "" && inputFile == "")) {
    FATAL("Specify input file or formula." << endl);
  }


  TFile out(outputFile.c_str(), "recreate");

  TH2D h0("h0", NULL,
          X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(),
          Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit());

  if        (formula != "") {

    TF2 f2("f2", formula.c_str());

    for (vector<JToken_t>::const_iterator i = parameters.begin(); i != parameters.end(); ++i) {
      f2.FixParameter(getParameter(*i), getValue(*i));
    }

    if (numberOfEvents > 0) {

      h0.Sumw2();
      h0.FillRandom(f2.GetName(), numberOfEvents);

    } else {

      for (Int_t ix = 1; ix <= h0.GetXaxis()->GetNbins(); ++ix) {
        for (Int_t iy = 1; iy <= h0.GetYaxis()->GetNbins(); ++iy) {
          h0.SetBinContent(ix, iy, f2.Eval(h0.GetXaxis()->GetBinCenter(ix),
                                           h0.GetYaxis()->GetBinCenter(iy)));
        }
      }
    }
  } else if (inputFile != "") {

    h0.Sumw2();

    ifstream in(inputFile.c_str());

    for (double x, y; in >> x >> y; ) {
      h0.Fill(x, y);
    }

    in.close();
  }

  out.Write();
  out.Close();
}