CMS 3D CMS Logo
edm::HydjetProducer

edm::HydjetProducer Class Reference

#include <GeneratorInterface/HydjetInterface/interface/HydjetProducer.h>

Inheritance diagram for edm::HydjetProducer:

edm::EDProducer edm::ProducerBase edm::ProductRegistryHelper

List of all members.

Public Member Functions

 HydjetProducer (const ParameterSet &)
virtual ~HydjetProducer ()

Private Member Functions

void add_heavy_ion_rec (HepMC::GenEvent *evt)
HepMC::GenParticle * build_hyjet (int index, int barcode)
HepMC::GenVertex * build_hyjet_vertex (int i, int id)
bool call_hyinit (double energy, double a, int ifb, double bmin, double bmax, double bfix, int nh)
bool call_pygive (const std::string &iParm)
void clear ()
bool get_hard_particles (HepMC::GenEvent *evt, std::vector< SubEvent > &subs)
bool get_soft_particles (HepMC::GenEvent *evt, std::vector< SubEvent > &subs)
bool hydjet_init (const ParameterSet &pset)
bool hyjpythia_init (const ParameterSet &pset)
double nuclear_radius () const
virtual void produce (Event &e, const EventSetup &es)

Private Attributes

double abeamtarget_
double bfixed_
double bmax_
double bmin_
int cflag_
double comenergy
bool docollisionalenloss_
 DEFAULT = true.
bool doradiativeenloss_
unsigned int eventNumber_
HepMC::GenEvent * evt
double fracsoftmult_
 DEFAULT = true.
CLHEP::HepRandomEngine * fRandomEngine
double hadfreeztemp_
std::string hymode_
unsigned int maxEventsToPrint_
double maxlongy_
double maxtrany_
int nhard_
int nmultiplicity_
unsigned int nquarkflavor_
int nsoft_
int nsub_
unsigned int pythiaPylistVerbosity_
 number of active quark flavors in qgp DEFAULT=0; allowed values: 0,1,2,3.
double qgpt0_
double qgptau0_
unsigned int shadowingswitch_
double signn_

Detailed Description

Definition at line 35 of file HydjetProducer.h.

Constructor & Destructor Documentation

HydjetProducer::HydjetProducer ( const ParameterSet &  pset  ) 

Definition at line 40 of file HydjetProducer.cc.

References abeamtarget_, bfixed_, bmax_, bmin_, call_hyinit(), cflag_, comenergy, edm::ParameterSet::getUntrackedParameter(), hydjet_init(), hyjpythia_init(), LogDebug, maxEventsToPrint_, nmultiplicity_, nuclear_radius(), and pythiaPylistVerbosity_.

00040                                                        :
00041     EDProducer(), evt(0), 
00042     abeamtarget_(pset.getParameter<double>("aBeamTarget")),
00043     bfixed_(pset.getParameter<double>("bFixed")),
00044     bmax_(pset.getParameter<double>("bMax")),
00045     bmin_(pset.getParameter<double>("bMin")),
00046     cflag_(pset.getParameter<int>("cFlag")),
00047     comenergy(pset.getParameter<double>("comEnergy")),
00048     doradiativeenloss_(pset.getParameter<bool>("doRadiativeEnLoss")),
00049     docollisionalenloss_(pset.getParameter<bool>("doCollisionalEnLoss")),
00050     fracsoftmult_(pset.getParameter<double>("fracSoftMultiplicity")),
00051     hadfreeztemp_(pset.getParameter<double>("hadronFreezoutTemperature")),
00052     hymode_(pset.getParameter<string>("hydjetMode")),
00053     maxEventsToPrint_(pset.getUntrackedParameter<int>("maxEventsToPrint", 1)),
00054     maxlongy_(pset.getParameter<double>("maxLongitudinalRapidity")),
00055     maxtrany_(pset.getParameter<double>("maxTransverseRapidity")),
00056     nsub_(0),
00057     nhard_(0),
00058     nmultiplicity_(pset.getParameter<int>("nMultiplicity")),
00059     nsoft_(0),
00060     nquarkflavor_(pset.getParameter<int>("qgpNumQuarkFlavor")),
00061     pythiaPylistVerbosity_(pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
00062     qgpt0_(pset.getParameter<double>("qgpInitialTemperature")),
00063     qgptau0_(pset.getParameter<double>("qgpProperTimeFormation")),
00064     shadowingswitch_(pset.getParameter<int>("shadowingSwitch")),
00065     signn_(pset.getParameter<double>("sigmaInelNN")),
00066     eventNumber_(0)
00067 {
00068   // Default constructor
00069 
00070   // PYLIST Verbosity Level
00071   // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
00072   pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0);
00073   LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_;
00074 
00075   //Max number of events printed on verbosity level 
00076   maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint",0);
00077   LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_;
00078 
00079   // the input impact parameter (bxx_) is in [fm]; transform in [fm/RA] for hydjet usage
00080   const float ra = nuclear_radius();
00081   LogInfo("RAScaling")<<"Nuclear radius(RA) =  "<<ra;
00082   bmin_     /= ra;
00083   bmax_     /= ra;
00084   bfixed_   /= ra;
00085 
00086   // initialize pythia
00087   hyjpythia_init(pset);
00088 
00089   // hydjet running options
00090   hydjet_init(pset);
00091 
00092   // initialize hydjet
00093   LogInfo("HYDJETinAction") << "##### Calling HYINIT("<<comenergy<<","<<abeamtarget_<<","<<cflag_<<","<<bmin_<<","<<bmax_<<","<<bfixed_<<","<<nmultiplicity_<<") ####";
00094 
00095   call_hyinit(comenergy,abeamtarget_,cflag_,bmin_,bmax_,bfixed_,nmultiplicity_);
00096     
00097   produces<HepMCProduct>();
00098   produces<GenHIEvent>();
00099 }

HydjetProducer::~HydjetProducer (  )  [virtual]

Definition at line 103 of file HydjetProducer.cc.

References clear().

00104 {
00105   // destructor
00106 
00107   call_pystat(1);
00108 
00109   clear();
00110 }

Member Function Documentation

void HydjetProducer::add_heavy_ion_rec ( HepMC::GenEvent *  evt  )  [private]

Definition at line 114 of file HydjetProducer.cc.

References hyfpar, hyjpar, nsub_, and nuclear_radius().

Referenced by produce().

00115 {
00116   // heavy ion record in the final CMSSW Event
00117 
00118   HepMC::HeavyIon* hi = new HepMC::HeavyIon(
00119     nsub_,                               // Ncoll_hard/N of SubEvents
00120     static_cast<int>(hyfpar.npart / 2),  // Npart_proj
00121     static_cast<int>(hyfpar.npart / 2),  // Npart_targ
00122     static_cast<int>(hyfpar.nbcol),      // Ncoll
00123     0,                                   // spectator_neutrons
00124     0,                                   // spectator_protons
00125     0,                                   // N_Nwounded_collisions
00126     0,                                   // Nwounded_N_collisions
00127     0,                                   // Nwounded_Nwounded_collisions
00128     hyfpar.bgen * nuclear_radius(),      // impact_parameter in [fm]
00129     0,                                   // event_plane_angle
00130     0,                                   // eccentricity
00131     hyjpar.sigin                         // sigma_inel_NN
00132   );
00133 
00134   evt->set_heavy_ion(*hi);
00135   delete hi;
00136 
00137 }

HepMC::GenParticle * HydjetProducer::build_hyjet ( int  index,
int  barcode 
) [private]

Definition at line 141 of file HydjetProducer.cc.

References hyjets, and p.

00142 {
00143    // Build particle object corresponding to index in hyjets (soft+hard)
00144 
00145    HepMC::GenParticle* p = new HepMC::GenParticle(
00146                                                   HepMC::FourVector(hyjets.phj[0][index],  // px
00147                                                                     hyjets.phj[1][index],  // py
00148                                                                     hyjets.phj[2][index],  // pz
00149                                                                     hyjets.phj[3][index]), // E
00150                                                   hyjets.khj[1][index],// id
00151                                                   hyjets.khj[0][index] // status
00152                                                   );
00153    p->suggest_barcode(barcode);
00154 
00155    return p;
00156 }

HepMC::GenVertex * HydjetProducer::build_hyjet_vertex ( int  i,
int  id 
) [private]

Definition at line 160 of file HydjetProducer.cc.

References hyjets, t, x, y, and z.

00161 {
00162   // build verteces for the hyjets stored events 
00163 
00164    double x=hyjets.vhj[0][i];
00165    double y=hyjets.vhj[1][i];
00166    double z=hyjets.vhj[2][i];
00167    double t=hyjets.vhj[4][i];
00168 
00169    HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(x,y,z,t),id);
00170    return vertex;
00171 }

bool HydjetProducer::call_hyinit ( double  energy,
double  a,
int  ifb,
double  bmin,
double  bmax,
double  bfix,
int  nh 
) [private]

Definition at line 294 of file HydjetProducer.cc.

References HYINIT, and pydatr.

Referenced by HydjetProducer().

00296 {
00297   // initialize hydjet  
00298 
00299    pydatr.mrpy[2]=1;
00300    HYINIT(energy,a,ifb,bmin,bmax,bfix,nh);
00301     return true;
00302 }

bool HydjetProducer::call_pygive ( const std::string &  iParm  )  [private]

Definition at line 175 of file HydjetProducer.cc.

References pydat1, and PYGIVE.

Referenced by hyjpythia_init().

00176 {
00177   // Set Pythia parameters
00178 
00179   int numWarn = pydat1.mstu[26]; //# warnings
00180   int numErr = pydat1.mstu[22];// # errors
00181   // call the fortran routine pygive with a fortran string
00182   PYGIVE( iParm.c_str(), iParm.length() );  
00183 
00184   // if an error or warning happens it is problem
00185   return pydat1.mstu[26] == numWarn && pydat1.mstu[22] == numErr;   
00186 }

void HydjetProducer::clear ( void   )  [private]

Definition at line 190 of file HydjetProducer.cc.

Referenced by ~HydjetProducer().

00191 {
00192 
00193 }

bool edm::HydjetProducer::get_hard_particles ( HepMC::GenEvent *  evt,
std::vector< SubEvent > &  subs 
) [private]

Referenced by produce().

bool edm::HydjetProducer::get_soft_particles ( HepMC::GenEvent *  evt,
std::vector< SubEvent > &  subs 
) [private]

Referenced by produce().

bool HydjetProducer::hydjet_init ( const ParameterSet &  pset  )  [private]

Definition at line 306 of file HydjetProducer.cc.

References docollisionalenloss_, doradiativeenloss_, fracsoftmult_, fRandomEngine, hadfreeztemp_, hyflow, hyjpar, hymode_, ludatr, maxlongy_, maxtrany_, nquarkflavor_, pyqpar, qgpt0_, qgptau0_, randomEngine, shadowingswitch_, and signn_.

Referenced by HydjetProducer().

00307 {
00308   // set hydjet options
00309 
00310   edm::Service<RandomNumberGenerator> rng;
00311   randomEngine = fRandomEngine = &(rng->getEngine());
00312   uint32_t seed = rng->mySeed();
00313   ludatr.mrlu[0]=seed;
00314 
00315   // hydjet running mode mode
00316   // kHydroOnly --- nhsel=0 jet production off (pure HYDRO event), nhsel=0
00317   // kHydroJets --- nhsle=1 jet production on, jet quenching off (HYDRO+njet*PYTHIA events)
00318   // kHydroQJet --- nhsel=2 jet production & jet quenching on (HYDRO+njet*PYQUEN events)
00319   // kJetsOnly  --- nhsel=3 jet production on, jet quenching off, HYDRO off (njet*PYTHIA events)
00320   // kQJetsOnly --- nhsel=4 jet production & jet quenching on, HYDRO off (njet*PYQUEN events)
00321 
00322   if(hymode_ == "kHydroOnly") hyjpar.nhsel=0;
00323   else if ( hymode_ == "kHydroJets") hyjpar.nhsel=1;
00324   else if ( hymode_ == "kHydroQJets") hyjpar.nhsel=2;
00325   else if ( hymode_ == "kJetsOnly") hyjpar.nhsel=3;
00326   else if ( hymode_ == "kQJetsOnly") hyjpar.nhsel=4;
00327   else  hyjpar.nhsel=2;
00328 
00329   // fraction of soft hydro induced multiplicity 
00330   hyflow.fpart =  fracsoftmult_; 
00331 
00332   // hadron freez-out temperature
00333   hyflow.Tf   = hadfreeztemp_;
00334 
00335   // maximum longitudinal collective rapidity
00336   hyflow.ylfl = maxlongy_;
00337   
00338   // maximum transverse collective rapidity
00339   hyflow.ytfl = maxtrany_;  
00340 
00341   // shadowing on=1, off=0
00342   hyjpar.ishad  = shadowingswitch_;
00343 
00344   // set inelastic nucleon-nucleon cross section
00345   hyjpar.sigin  = signn_;
00346 
00347   // number of active quark flavors in qgp
00348   pyqpar.nfu    = nquarkflavor_;
00349 
00350   // initial temperature of QGP
00351   pyqpar.T0u    = qgpt0_;
00352 
00353   // proper time of QGP formation
00354   pyqpar.tau0u  = qgptau0_;
00355 
00356   // type of medium induced partonic energy loss
00357   if( doradiativeenloss_ && docollisionalenloss_ ){
00358     edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
00359     pyqpar.ienglu = 0; 
00360   } else if ( doradiativeenloss_ ) {
00361     edm::LogInfo("HydjetenLoss") << "##### Only RADIATIVE partonic energy loss ON ####";
00362     pyqpar.ienglu = 1; 
00363   } else if ( docollisionalenloss_ ) {
00364     edm::LogInfo("HydjetEnLoss") << "##### Only COLLISIONAL partonic energy loss ON ####";
00365     pyqpar.ienglu = 2; 
00366   } else {
00367     edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
00368     pyqpar.ienglu = 0; 
00369   }
00370 
00371   return true;
00372 }

bool HydjetProducer::hyjpythia_init ( const ParameterSet &  pset  )  [private]

Definition at line 376 of file HydjetProducer.cc.

References call_pygive(), edm::errors::Configuration, GenMuonPlsPt100GeV_cfg::cout, lat::endl(), edm::ParameterSet::getParameter(), i, and pars.

Referenced by HydjetProducer().

00377 {
00378   //Set PYTHIA parameters
00379 
00380   //random number seed
00381   edm::Service<RandomNumberGenerator> rng;
00382   uint32_t seed = rng->mySeed();
00383   ostringstream sRandomSet;
00384   sRandomSet << "MRPY(1)=" << seed;
00385   call_pygive(sRandomSet.str());
00386 
00387     // Set PYTHIA parameters in a single ParameterSet
00388   ParameterSet pythia_params = pset.getParameter<ParameterSet>("PythiaParameters") ;
00389   // The parameter sets to be read 
00390   vector<string> setNames = pythia_params.getParameter<vector<string> >("parameterSets");
00391 
00392     // Loop over the sets
00393   for ( unsigned i=0; i<setNames.size(); ++i ) {
00394     string mySet = setNames[i];
00395     
00396     // Read the PYTHIA parameters for each set of parameters
00397     vector<string> pars = pythia_params.getParameter<vector<string> >(mySet);
00398     
00399     cout << "----------------------------------------------" << endl;
00400     cout << "Read PYTHIA parameter set " << mySet << endl;
00401     cout << "----------------------------------------------" << endl;
00402     
00403     // Loop over all parameters and stop in case of mistake
00404     for( vector<string>::const_iterator itPar = pars.begin(); itPar != pars.end(); ++itPar ) {
00405       static string sRandomValueSetting("MRPY(1)");
00406       if( 0 == itPar->compare(0,sRandomValueSetting.size(),sRandomValueSetting) ) {
00407         throw edm::Exception(edm::errors::Configuration,"PythiaError")
00408           <<" Attempted to set random number using 'MRPY(1)'. NOT ALLOWED! \n Use RandomNumberGeneratorService to set the random number seed.";
00409       }
00410       if( !call_pygive(*itPar) ) {
00411         throw edm::Exception(edm::errors::Configuration,"PythiaError") 
00412           <<"PYTHIA did not accept \""<<*itPar<<"\"";
00413       }
00414     }
00415   }
00416 
00417   return true;
00418 }

double HydjetProducer::nuclear_radius (  )  const [inline, private]

Definition at line 104 of file HydjetProducer.h.

References funct::pow().

Referenced by add_heavy_ion_rec(), and HydjetProducer().

00105 {
00106   // Return the nuclear radius derived from the 
00107   // beam/target atomic mass number.
00108 
00109   return 1.15 * pow((double)abeamtarget_, 1./3.);
00110 }

void HydjetProducer::produce ( Event &  e,
const EventSetup &  es 
) [private, virtual]

Implements edm::EDProducer.

Definition at line 422 of file HydjetProducer.cc.

References add_heavy_ion_rec(), call_pylist(), edm::EventID::event(), edm::errors::EventCorruption, eventNumber_, except, get_hard_particles(), get_soft_particles(), HYEVNT, hyflow, hyfpar, hyjpar, edm::Event::id(), maxEventsToPrint_, nhard_, edm::Event::put(), pypars, pyqpar, and pythiaPylistVerbosity_.

00423 {
00424   // generate single event
00425 
00426   nsoft_    = 0;
00427   nhard_    = 0;
00428 
00429   edm::LogInfo("HYDJETmode") << "##### HYDJET  nhsel = " << hyjpar.nhsel;
00430   edm::LogInfo("HYDJETfpart") << "##### HYDJET fpart = " << hyflow.fpart;
00431   edm::LogInfo("HYDJETtf") << "##### HYDJET hadron freez-out temp, Tf = " << hyflow.Tf;
00432   edm::LogInfo("HYDJETinTemp") << "##### HYDJET: QGP init temperature, T0 ="<<pyqpar.T0u;
00433   edm::LogInfo("HYDJETinTau") << "##### HYDJET: QGP formation time,tau0 ="<<pyqpar.tau0u;
00434 
00435 
00436   // generate a HYDJET event
00437   int ntry = 0;
00438   while(nsoft_ == 0 && nhard_ == 0){
00439      if(ntry > 100){
00440        edm::LogError("HydjetEmptyEvent") << "##### HYDJET: No Particles generated, Number of tries ="<<ntry;
00441 
00442 // Throw an exception.  Use the EventCorruption exception since it maps onto SkipEvent
00443 // which is what we want to do here.
00444 
00445        std::ostringstream sstr;
00446        sstr << "HydjetProducerProducer: No particles generated after " << ntry << " tries.\n";
00447        edm::Exception except(edm::errors::EventCorruption, sstr.str());
00448        throw except;
00449      } else {
00450        HYEVNT();
00451        nsoft_    = hyfpar.nhyd;
00452        nsub_     = hyjpar.njet;
00453        nhard_    = hyfpar.npyt;
00454        ++ntry;
00455      }
00456   }
00457 
00458   std::vector<SubEvent> subvector;
00459 
00460   // event information
00461   HepMC::GenEvent *evt = new HepMC::GenEvent();
00462 
00463   if(nhard_>0) get_hard_particles(evt,subvector); 
00464   if(nsoft_>0) get_soft_particles(evt,subvector);
00465 
00466   evt->set_signal_process_id(pypars.msti[0]);      // type of the process
00467   evt->set_event_scale(pypars.pari[16]);           // Q^2
00468   ++eventNumber_;
00469   evt->set_event_number(eventNumber_);
00470 
00471   add_heavy_ion_rec(evt);
00472 
00473   auto_ptr<HepMCProduct> bare_product(new HepMCProduct());
00474   bare_product->addHepMCData(evt );
00475   auto_ptr<GenHIEvent> genhi(new GenHIEvent(subvector,hyjpar.nhsel));
00476   e.put(bare_product);
00477   e.put(genhi);
00478 
00479   // print PYLIST info
00480   if (e.id().event() <= maxEventsToPrint_ && pythiaPylistVerbosity_)     
00481      call_pylist(pythiaPylistVerbosity_);      
00482 }

Member Data Documentation

double edm::HydjetProducer::abeamtarget_ [private]

Definition at line 56 of file HydjetProducer.h.

Referenced by HydjetProducer().

double edm::HydjetProducer::bfixed_ [private]

Definition at line 57 of file HydjetProducer.h.

Referenced by HydjetProducer().

double edm::HydjetProducer::bmax_ [private]

Definition at line 58 of file HydjetProducer.h.

Referenced by HydjetProducer().

double edm::HydjetProducer::bmin_ [private]

Definition at line 60 of file HydjetProducer.h.

Referenced by HydjetProducer().

int edm::HydjetProducer::cflag_ [private]

Definition at line 62 of file HydjetProducer.h.

Referenced by HydjetProducer().

double edm::HydjetProducer::comenergy [private]

Definition at line 65 of file HydjetProducer.h.

Referenced by HydjetProducer().

bool edm::HydjetProducer::docollisionalenloss_ [private]

DEFAULT = true.

Definition at line 67 of file HydjetProducer.h.

Referenced by hydjet_init().

bool edm::HydjetProducer::doradiativeenloss_ [private]

Definition at line 66 of file HydjetProducer.h.

Referenced by hydjet_init().

unsigned int edm::HydjetProducer::eventNumber_ [private]

Definition at line 101 of file HydjetProducer.h.

Referenced by produce().

HepMC::GenEvent* edm::HydjetProducer::evt [private]

Definition at line 55 of file HydjetProducer.h.

double edm::HydjetProducer::fracsoftmult_ [private]

DEFAULT = true.

Definition at line 68 of file HydjetProducer.h.

Referenced by hydjet_init().

CLHEP::HepRandomEngine* edm::HydjetProducer::fRandomEngine [private]

Definition at line 100 of file HydjetProducer.h.

Referenced by hydjet_init().

double edm::HydjetProducer::hadfreeztemp_ [private]

Definition at line 74 of file HydjetProducer.h.

Referenced by hydjet_init().

std::string edm::HydjetProducer::hymode_ [private]

Definition at line 76 of file HydjetProducer.h.

Referenced by hydjet_init().

unsigned int edm::HydjetProducer::maxEventsToPrint_ [private]

Definition at line 77 of file HydjetProducer.h.

Referenced by HydjetProducer(), and produce().

double edm::HydjetProducer::maxlongy_ [private]

Definition at line 78 of file HydjetProducer.h.

Referenced by hydjet_init().

double edm::HydjetProducer::maxtrany_ [private]

Definition at line 81 of file HydjetProducer.h.

Referenced by hydjet_init().

int edm::HydjetProducer::nhard_ [private]

Definition at line 85 of file HydjetProducer.h.

Referenced by produce().

int edm::HydjetProducer::nmultiplicity_ [private]

Definition at line 86 of file HydjetProducer.h.

Referenced by HydjetProducer().

unsigned int edm::HydjetProducer::nquarkflavor_ [private]

Definition at line 89 of file HydjetProducer.h.

Referenced by hydjet_init().

int edm::HydjetProducer::nsoft_ [private]

Definition at line 88 of file HydjetProducer.h.

int edm::HydjetProducer::nsub_ [private]

Definition at line 84 of file HydjetProducer.h.

Referenced by add_heavy_ion_rec().

unsigned int edm::HydjetProducer::pythiaPylistVerbosity_ [private]

number of active quark flavors in qgp DEFAULT=0; allowed values: 0,1,2,3.

Definition at line 91 of file HydjetProducer.h.

Referenced by HydjetProducer(), and produce().

double edm::HydjetProducer::qgpt0_ [private]

Definition at line 92 of file HydjetProducer.h.

Referenced by hydjet_init().

double edm::HydjetProducer::qgptau0_ [private]

Definition at line 94 of file HydjetProducer.h.

Referenced by hydjet_init().

unsigned int edm::HydjetProducer::shadowingswitch_ [private]

Definition at line 96 of file HydjetProducer.h.

Referenced by hydjet_init().

double edm::HydjetProducer::signn_ [private]

Definition at line 98 of file HydjetProducer.h.

Referenced by hydjet_init().

The documentation for this class was generated from the following files:
Generated on Tue Jun 9 18:41:14 2009 for CMSSW by  doxygen 1.5.4