JPythia.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

#include "JSirene/JPythia.hh"
#include "JGizmo/JManager.hh"

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


/**
 * \file
 * Auxiliary program histogram plot EM-equivalent energy for hadrons.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string   outputFile;
  double   epsilon;
  bool     debug;

  try {

    JParser<> zap("Auxiliary program histogram plot EM-equivalent energy for hadrons."\
                  "\nNote that the results have a small offset to effectively be visualised.");

    zap['o'] = make_field(outputFile)  = "pythia.root";
    zap['e'] = make_field(epsilon)     = 0.004;
    zap['d'] = make_field(debug);
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }



  JGIZMO::JManager<int, TH1D> zmap(new TH1D("E[%]",  NULL, 500, -0.8, 9.0));

  const int buffer[] = {
    TRACK_TYPE_PHOTON,
    TRACK_TYPE_NEUTRAL_PION,
    TRACK_TYPE_NEUTRAL_ANTIPION,
    TRACK_TYPE_ELECTRON,
    TRACK_TYPE_ANTIELECTRON,
    TRACK_TYPE_CHARGED_PION_PLUS,
    TRACK_TYPE_CHARGED_PION_MINUS,
    -1 };

  map<int, double> Y;

  for (const int *type = buffer; *type != -1; ++type) {
    Y[*type] = epsilon * (type - buffer);
  }
  
  for (int i = 1; i <= zmap->GetNbinsX(); ++i) {

    const double x = zmap->GetBinCenter(i);
    const double E = pow(10.0, x);

    for (const int *type = buffer; *type != -1; ++type) {
      zmap[*type]->SetBinContent(i, pythia(*type, E)/E + Y[*type]);
    }
  }

  zmap.Write(outputFile.c_str());
}