JGizmo/JHistogram3D.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH3D.h"
#include "TF3.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 TH2D 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;
  string            title;
  JHistogram_t      X;
  JHistogram_t      Y;
  JHistogram_t      Z;  
  bool              sumw2;
  int               debug;

  try {

    JParser<> zap("Program to create TH2D 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['T'] = make_field(title)           = "h0";
    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['z'] = make_field(Z)               = JHistogram_t(100, -1.0, +1.0);    
    zap['s'] = make_field(sumw2);
    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");

  const string::size_type pos = title.rfind('/');

  if (pos != string::npos) {

    const string dir = title.substr(0, pos);

    out.mkdir(dir.c_str());
    out.cd   (dir.c_str());

    title = title.substr(pos + 1);
  }

  TH3D h0(title.c_str(), NULL,
          X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(),
          Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit(),
          Z.getNumberOfBins(), Z.getLowerLimit(), Z.getUpperLimit());

  if        (formula != "") {

    TF3 f3("f3", formula.c_str());

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

    if (numberOfEvents > 0) {

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

    } else {

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

    if (sumw2) {
      h0.Sumw2();
    }

    ifstream in(inputFile.c_str());

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

    in.close();
  }

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