Pythia6Interface/src/PtYDistributor.cc

0001
0002 #include "GeneratorInterface/Pythia6Interface/interface/PtYDistributor.h"
0003
0004 #include "FWCore/ParameterSet/interface/FileInPath.h"
0005 #include "FWCore/Utilities/interface/EDMException.h"
0006 #include "TFile.h"
0007 #include "TGraph.h"
0008
0009 #include "CLHEP/Random/RandGeneral.h"
0010
0011 #include <iostream>
0012 #include <string>
0013
0014 using namespace gen;
0015
0016 PtYDistributor::PtYDistributor(const edm::FileInPath& fip,
0017                                double ptmax = 100,
0018                                double ptmin = 0,
0019                                double ymax = 10,
0020                                double ymin = -10,
0021                                int ptbins = 1000,
0022                                int ybins = 50)
0023     : ptmax_(ptmax), ptmin_(ptmin), ymax_(ymax), ymin_(ymin), ptbins_(ptbins), ybins_(ybins) {
0024   // edm::FileInPath f1(input);
0025   const std::string& fDataFile = fip.fullPath();
0026
0027   std::cout << " File from " << fDataFile << std::endl;
0028   TFile f(fDataFile.c_str(), "READ");
0029   TGraph* yfunc = (TGraph*)f.Get("rapidity");
0030   TGraph* ptfunc = (TGraph*)f.Get("pt");
0031
0032   if (!yfunc)
0033     throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
0034         << "Rapidity distribution could not be found in file " << fDataFile;
0035
0036   if (!ptfunc)
0037     throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
0038         << "Pt distribution could not be found in file " << fDataFile;
0039
0040   if (ptbins_ > 100000) {
0041     ptbins_ = 100000;
0042   }
0043
0044   if (ybins_ > 100000) {
0045     ybins_ = 100000;
0046   }
0047
0048   double aProbFunc1[100000];
0049   double aProbFunc2[100000];
0050
0051   for (int i = 0; i < ybins_; ++i) {
0052     double xy = ymin_ + i * (ymax_ - ymin_) / ybins_;
0053     double yy = yfunc->Eval(xy);
0054     aProbFunc1[i] = yy;
0055   }
0056
0057   for (int ip = 0; ip < ptbins_; ++ip) {
0058     double xpt = ptmin_ + ip * (ptmax_ - ptmin_) / ptbins_;
0059     double ypt = ptfunc->Eval(xpt);
0060     aProbFunc2[ip] = ypt;
0061   }
0062
0063   fYGenerator = new CLHEP::RandGeneral(nullptr, aProbFunc1, ybins_);
0064   fPtGenerator = new CLHEP::RandGeneral(nullptr, aProbFunc2, ptbins_);
0065
0066   f.Close();
0067
0068 }  // from file
0069
0070 double PtYDistributor::fireY(CLHEP::HepRandomEngine* engine) { return fireY(ymin_, ymax_, engine); }
0071
0072 double PtYDistributor::firePt(CLHEP::HepRandomEngine* engine) { return firePt(ptmin_, ptmax_, engine); }
0073
0074 double PtYDistributor::fireY(double ymin, double ymax, CLHEP::HepRandomEngine* engine) {
0075   double y = -999;
0076
0077   if (fYGenerator) {
0078     while (y < ymin || y > ymax)
0079       y = ymin_ + (ymax_ - ymin_) * fYGenerator->shoot(engine);
0080   } else {
0081     throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
0082         << "Random y requested but Random Number Generator for y not Initialized!";
0083   }
0084   return y;
0085 }
0086
0087 double PtYDistributor::firePt(double ptmin, double ptmax, CLHEP::HepRandomEngine* engine) {
0088   double pt = -999;
0089
0090   if (fPtGenerator) {
0091     while (pt < ptmin || pt > ptmax)
0092       pt = ptmin_ + (ptmax_ - ptmin_) * fPtGenerator->shoot(engine);
0093   } else {
0094     throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
0095         << "Random pt requested but Random Number Generator for pt not Initialized!";
0096   }
0097   return pt;
0098 }