CMS 3D CMS Logo

CMSSW_4_4_3_patch1/src/GeneratorInterface/Pythia6Interface/src/PtYDistributor.cc

Go to the documentation of this file.
00001 
00002 #include "GeneratorInterface/Pythia6Interface/interface/PtYDistributor.h"
00003 //#include "FWCore/ParameterSet/interface/FileInPath.h"
00004 #include "FWCore/Utilities/interface/EDMException.h"
00005 #include "TFile.h"
00006 #include "TGraph.h"
00007 
00008 using namespace gen;
00009 
00010 PtYDistributor::PtYDistributor(edm::FileInPath fip, CLHEP::HepRandomEngine& fRandomEngine, 
00011                                double ptmax = 100, double ptmin = 0, 
00012                                double ymax = 10, double ymin = -10, 
00013                                int ptbins = 1000, int ybins = 50) 
00014    : ptmax_(ptmax),ptmin_(ptmin),ymax_(ymax),ymin_(ymin), ptbins_(ptbins), ybins_(ybins)
00015 {  
00016    // edm::FileInPath f1(input);
00017    std::string fDataFile = fip.fullPath();
00018 
00019    std::cout<<" File from "<<fDataFile <<std::endl;
00020    TFile f(fDataFile.c_str(),"READ");
00021    TGraph* yfunc = (TGraph*)f.Get("rapidity"); 
00022    TGraph* ptfunc = (TGraph*)f.Get("pt");
00023 
00024    if( !yfunc ) 
00025       throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
00026          <<"Rapidity distribution could not be found in file "<<fDataFile;
00027 
00028    if( !ptfunc )
00029       throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
00030          <<"Pt distribution could not be found in file "<<fDataFile;
00031 
00032    if(ptbins_ > 100000){
00033       ptbins_ = 100000;
00034    }
00035 
00036    if(ybins_ > 100000){
00037       ybins_ = 100000;
00038    }
00039 
00040    double aProbFunc1[100000];
00041    double aProbFunc2[100000];
00042 
00043    for(int i = 0; i < ybins_; ++i){
00044       double xy = ymin_+i*(ymax_-ymin_)/ybins_;
00045       double yy = yfunc->Eval(xy);   
00046       aProbFunc1[i] = yy;
00047    }
00048 
00049    for(int ip = 0; ip < ptbins_; ++ip){
00050       double xpt = ptmin_+ip*(ptmax_-ptmin_)/ptbins_;
00051       double ypt = ptfunc->Eval(xpt);
00052       aProbFunc2[ip] = ypt;
00053    }
00054 
00055   fYGenerator = new CLHEP::RandGeneral(fRandomEngine,aProbFunc1,ybins_);
00056   fPtGenerator = new CLHEP::RandGeneral(fRandomEngine,aProbFunc2,ptbins_);
00057 
00058   f.Close();
00059   
00060 } // from file
00061 
00062 double
00063 PtYDistributor::fireY(){
00064    return fireY(ymin_,ymax_);
00065 }
00066 
00067 double
00068 PtYDistributor::firePt(){
00069    return firePt(ptmin_,ptmax_);
00070 }
00071 
00072 double
00073 PtYDistributor::fireY(double ymin, double ymax){
00074 
00075    double y = -999;
00076 
00077    if(fYGenerator){
00078       while(y < ymin || y > ymax)
00079          y = ymin_+(ymax_-ymin_)*fYGenerator->fire();
00080    }else{
00081       throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
00082          <<"Random y requested but Random Number Generator for y not Initialized!";
00083    }
00084    return y;
00085 }
00086 
00087 double 
00088 PtYDistributor::firePt(double ptmin, double ptmax){
00089 
00090    double pt = -999;
00091 
00092    if(fPtGenerator){
00093       while(pt < ptmin || pt > ptmax)
00094          pt = ptmin_+(ptmax_-ptmin_)*fPtGenerator->fire();
00095    }else{
00096       throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
00097          <<"Random pt requested but Random Number Generator for pt not Initialized!";
00098    }
00099    return pt;
00100 }