PtYDistributor.cc

Go to the documentation of this file.

#include "GeneratorInterface/Pythia6Interface/interface/PtYDistributor.h"

#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "TFile.h"
#include "TGraph.h"

#include "CLHEP/Random/RandGeneral.h"

#include <iostream>
#include <string>

using namespace gen;

PtYDistributor::PtYDistributor(const edm::FileInPath& fip,
                               double ptmax = 100, double ptmin = 0, 
                   double ymax = 10, double ymin = -10, 
                   int ptbins = 1000, int ybins = 50) 
   : ptmax_(ptmax),ptmin_(ptmin),ymax_(ymax),ymin_(ymin), ptbins_(ptbins), ybins_(ybins)
{  
   // edm::FileInPath f1(input);
   std::string fDataFile = fip.fullPath();

   std::cout<<" File from "<<fDataFile <<std::endl;
   TFile f(fDataFile.c_str(),"READ");
   TGraph* yfunc = (TGraph*)f.Get("rapidity"); 
   TGraph* ptfunc = (TGraph*)f.Get("pt");

   if( !yfunc ) 
      throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
         <<"Rapidity distribution could not be found in file "<<fDataFile;

   if( !ptfunc )
      throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
         <<"Pt distribution could not be found in file "<<fDataFile;

   if(ptbins_ > 100000){
      ptbins_ = 100000;
   }

   if(ybins_ > 100000){
      ybins_ = 100000;
   }

   double aProbFunc1[100000];
   double aProbFunc2[100000];

   for(int i = 0; i < ybins_; ++i){
      double xy = ymin_+i*(ymax_-ymin_)/ybins_;
      double yy = yfunc->Eval(xy);   
      aProbFunc1[i] = yy;
   }

   for(int ip = 0; ip < ptbins_; ++ip){
      double xpt = ptmin_+ip*(ptmax_-ptmin_)/ptbins_;
      double ypt = ptfunc->Eval(xpt);
      aProbFunc2[ip] = ypt;
   }

  fYGenerator = new CLHEP::RandGeneral(nullptr,aProbFunc1,ybins_);
  fPtGenerator = new CLHEP::RandGeneral(nullptr,aProbFunc2,ptbins_);

  f.Close();
  
} // from file

double
PtYDistributor::fireY(CLHEP::HepRandomEngine* engine){
  return fireY(ymin_,ymax_, engine);
}

double
PtYDistributor::firePt(CLHEP::HepRandomEngine* engine){
  return firePt(ptmin_,ptmax_, engine);
}

double
PtYDistributor::fireY(double ymin, double ymax, CLHEP::HepRandomEngine* engine){

   double y = -999;

   if(fYGenerator){
      while(y < ymin || y > ymax)
     y = ymin_+(ymax_-ymin_)*fYGenerator->shoot(engine);
   }else{
      throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
     <<"Random y requested but Random Number Generator for y not Initialized!";
   }
   return y;
}

double 
PtYDistributor::firePt(double ptmin, double ptmax, CLHEP::HepRandomEngine* engine){

   double pt = -999;

   if(fPtGenerator){
      while(pt < ptmin || pt > ptmax)
         pt = ptmin_+(ptmax_-ptmin_)*fPtGenerator->shoot(engine);
   }else{
      throw edm::Exception(edm::errors::NullPointerError,"PtYDistributor")
         <<"Random pt requested but Random Number Generator for pt not Initialized!";
   }
   return pt;
}