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 "TRandom3.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 TH3D 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;
  size_t            numberOfEvents;
  string            title;
  JHistogram_t      X;
  JHistogram_t      Y;
  JHistogram_t      Z;  
  bool              logx;
  bool              logy;
  bool              logz;
  bool              sumw2;
  UInt_t            seed;
  int               debug;
 
  try {
 
    JParser<> zap("Program to create TH3D 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['X'] = make_field(logx);
    zap['Y'] = make_field(logy);
    zap['Z'] = make_field(logz);
    zap['s'] = make_field(sumw2);
    zap['S'] = make_field(seed)            = 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);
  }
 
  gRandom->SetSeed(seed);
 
  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(),
           logx ? pow(10.0, X.getLowerLimit()) : X.getLowerLimit(),
           logx ? pow(10.0, X.getUpperLimit()) : X.getUpperLimit(),
           logy ? pow(10.0, Y.getLowerLimit()) : Y.getLowerLimit(),
           logy ? pow(10.0, Y.getUpperLimit()) : Y.getUpperLimit(),
           logz ? pow(10.0, Z.getLowerLimit()) : Z.getLowerLimit(),
           logz ? pow(10.0, Z.getUpperLimit()) : Z.getUpperLimit());
 
    for (vector<JToken_t>::const_iterator i = parameters.begin(); i != parameters.end(); ++i) {
      f3.FixParameter(getParameter(*i), getValue(*i));
    }
 
    if (numberOfEvents > 0) {
 
      h0.Sumw2();
 
      Double_t x, y, z;
 
      for (size_t i = 0; i != numberOfEvents; ++i) {
 
        f3.GetRandom3(x, y, z, gRandom);
 
        h0.Fill(logx ? log10(x) : x,
                logy ? log10(y) : y,
                logz ? log10(z) : z);
      }
 
    } 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();
}