CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
PtYDistributor.cc
Go to the documentation of this file.
1 
3 
6 #include "TFile.h"
7 #include "TGraph.h"
8 
9 #include "CLHEP/Random/RandGeneral.h"
10 
11 #include <iostream>
12 #include <string>
13 
14 using namespace gen;
15 
17  double ptmax = 100,
18  double ptmin = 0,
19  double ymax = 10,
20  double ymin = -10,
21  int ptbins = 1000,
22  int ybins = 50)
23  : ptmax_(ptmax), ptmin_(ptmin), ymax_(ymax), ymin_(ymin), ptbins_(ptbins), ybins_(ybins) {
24  // edm::FileInPath f1(input);
25  std::string fDataFile = fip.fullPath();
26 
27  std::cout << " File from " << fDataFile << std::endl;
28  TFile f(fDataFile.c_str(), "READ");
29  TGraph* yfunc = (TGraph*)f.Get("rapidity");
30  TGraph* ptfunc = (TGraph*)f.Get("pt");
31 
32  if (!yfunc)
33  throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
34  << "Rapidity distribution could not be found in file " << fDataFile;
35 
36  if (!ptfunc)
37  throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
38  << "Pt distribution could not be found in file " << fDataFile;
39 
40  if (ptbins_ > 100000) {
41  ptbins_ = 100000;
42  }
43 
44  if (ybins_ > 100000) {
45  ybins_ = 100000;
46  }
47 
48  double aProbFunc1[100000];
49  double aProbFunc2[100000];
50 
51  for (int i = 0; i < ybins_; ++i) {
52  double xy = ymin_ + i * (ymax_ - ymin_) / ybins_;
53  double yy = yfunc->Eval(xy);
54  aProbFunc1[i] = yy;
55  }
56 
57  for (int ip = 0; ip < ptbins_; ++ip) {
58  double xpt = ptmin_ + ip * (ptmax_ - ptmin_) / ptbins_;
59  double ypt = ptfunc->Eval(xpt);
60  aProbFunc2[ip] = ypt;
61  }
62 
63  fYGenerator = new CLHEP::RandGeneral(nullptr, aProbFunc1, ybins_);
64  fPtGenerator = new CLHEP::RandGeneral(nullptr, aProbFunc2, ptbins_);
65 
66  f.Close();
67 
68 } // from file
69 
70 double PtYDistributor::fireY(CLHEP::HepRandomEngine* engine) { return fireY(ymin_, ymax_, engine); }
71 
72 double PtYDistributor::firePt(CLHEP::HepRandomEngine* engine) { return firePt(ptmin_, ptmax_, engine); }
73 
74 double PtYDistributor::fireY(double ymin, double ymax, CLHEP::HepRandomEngine* engine) {
75  double y = -999;
76 
77  if (fYGenerator) {
78  while (y < ymin || y > ymax)
79  y = ymin_ + (ymax_ - ymin_) * fYGenerator->shoot(engine);
80  } else {
81  throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
82  << "Random y requested but Random Number Generator for y not Initialized!";
83  }
84  return y;
85 }
86 
87 double PtYDistributor::firePt(double ptmin, double ptmax, CLHEP::HepRandomEngine* engine) {
88  double pt = -999;
89 
90  if (fPtGenerator) {
91  while (pt < ptmin || pt > ptmax)
92  pt = ptmin_ + (ptmax_ - ptmin_) * fPtGenerator->shoot(engine);
93  } else {
94  throw edm::Exception(edm::errors::NullPointerError, "PtYDistributor")
95  << "Random pt requested but Random Number Generator for pt not Initialized!";
96  }
97  return pt;
98 }
std::string fullPath() const
Definition: FileInPath.cc:161
double firePt(CLHEP::HepRandomEngine *)
double fireY(CLHEP::HepRandomEngine *)
double f[11][100]
CLHEP::RandGeneral * fYGenerator
CLHEP::RandGeneral * fPtGenerator
double ptmin
Definition: HydjetWrapper.h:84