Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:24

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   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 }