File indexing completed on 2025-01-09 23:33:35
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
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 }
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 }