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 JLANG::JToken<';'>           JToken_t;
  typedef JAbstractHistogram<Double_t> JHistogram_t;

  string            outputFile;
  string            formula;
  vector<JToken_t>  parameterValues;
  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(formula);
    zap['='] = make_field(parameterValues) = 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);
  }


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

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


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

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

  h0.Sumw2();

  if (numberOfEvents > 0) {

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


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