PtYDistributor.cc

Go to the documentation of this file.

#include "GeneratorInterface/Pythia6Interface/interface/PtYDistributor.h"

#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "TFile.h"
#include "TGraph.h"

#include "CLHEP/Random/RandGeneral.h"

#include <iostream>
#include <string>

using namespace gen;

PtYDistributor::PtYDistributor(const edm::FileInPath& fip,
                               double ptmax = 100,
                               double ptmin = 0,
                               double ymax = 10,
                               double ymin = -10,
                               int ptbins = 1000,
                               int ybins = 50)
    : ptmax_(ptmax), ptmin_(ptmin), ymax_(ymax), ymin_(ymin), ptbins_(ptbins), ybins_(ybins) {
  // edm::FileInPath f1(input);
  std::string fDataFile = fip.fullPath();

  std::cout << " File from " << fDataFile << std::endl;
  TFile f(fDataFile.c_str(), "READ");
  TGraph* yfunc = (TGraph*)f.Get("rapidity");
  TGraph* ptfunc = (TGraph*)f.Get("pt");

  if (!yfunc)
    throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
        << "Rapidity distribution could not be found in file " << fDataFile;

  if (!ptfunc)
    throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
        << "Pt distribution could not be found in file " << fDataFile;

  if (ptbins_ > 100000) {
    ptbins_ = 100000;
  }

  if (ybins_ > 100000) {
    ybins_ = 100000;
  }

  double aProbFunc1[100000];
  double aProbFunc2[100000];

  for (int i = 0; i < ybins_; ++i) {
    double xy = ymin_ + i * (ymax_ - ymin_) / ybins_;
    double yy = yfunc->Eval(xy);
    aProbFunc1[i] = yy;
  }

  for (int ip = 0; ip < ptbins_; ++ip) {
    double xpt = ptmin_ + ip * (ptmax_ - ptmin_) / ptbins_;
    double ypt = ptfunc->Eval(xpt);
    aProbFunc2[ip] = ypt;
  }

  fYGenerator = new CLHEP::RandGeneral(nullptr, aProbFunc1, ybins_);
  fPtGenerator = new CLHEP::RandGeneral(nullptr, aProbFunc2, ptbins_);

  f.Close();

}  // from file

double PtYDistributor::fireY(CLHEP::HepRandomEngine* engine) { return fireY(ymin_, ymax_, engine); }

double PtYDistributor::firePt(CLHEP::HepRandomEngine* engine) { return firePt(ptmin_, ptmax_, engine); }

double PtYDistributor::fireY(double ymin, double ymax, CLHEP::HepRandomEngine* engine) {
  double y = -999;

  if (fYGenerator) {
    while (y < ymin || y > ymax)
      y = ymin_ + (ymax_ - ymin_) * fYGenerator->shoot(engine);
  } else {
    throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
        << "Random y requested but Random Number Generator for y not Initialized!";
  }
  return y;
}

double PtYDistributor::firePt(double ptmin, double ptmax, CLHEP::HepRandomEngine* engine) {
  double pt = -999;

  if (fPtGenerator) {
    while (pt < ptmin || pt > ptmax)
      pt = ptmin_ + (ptmax_ - ptmin_) * fPtGenerator->shoot(engine);
  } else {
    throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
        << "Random pt requested but Random Number Generator for pt not Initialized!";
  }
  return pt;
}