CMSSW Reference Manual

00001 /*
00002  *
00003  * Generates PYQUEN HepMC events
00004  *
00005  * Original Author: Camelia Mironov
00006  * $Id: PyquenProducer.cc,v 1.2 2008/04/21 16:51:48 yilmaz Exp $
00007 */
00008 
00009 #include <iostream>
00010 #include "time.h"
00011 
00012 #include "GeneratorInterface/PyquenInterface/interface/PyquenProducer.h"
00013 #include "GeneratorInterface/PyquenInterface/interface/PyquenWrapper.h"
00014 #include "GeneratorInterface/PyquenInterface/interface/PYR.h"
00015 #include "GeneratorInterface/CommonInterface/interface/PythiaCMS.h"
00016 
00017 #include "SimDataFormats/HepMCProduct/interface/GenInfoProduct.h"
00018 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
00019 #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
00020 #include "FWCore/Framework/interface/Event.h"
00021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00022 #include "FWCore/ServiceRegistry/interface/Service.h"
00023 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00024 #include "CLHEP/Random/RandomEngine.h"
00025 
00026 #include "HepMC/IO_HEPEVT.h"
00027 #include "HepMC/PythiaWrapper.h"
00028 
00029 using namespace edm;
00030 using namespace std;
00031 
00032 HepMC::IO_HEPEVT hepevtio2;
00033 
00034 PyquenProducer :: PyquenProducer(const ParameterSet & pset):
00035 EDProducer(), evt(0), 
00036 abeamtarget_(pset.getParameter<double>("aBeamTarget")),
00037 angularspecselector_(pset.getParameter<int>("angularSpectrumSelector")),
00038 bfixed_(pset.getParameter<double>("bFixed")),
00039 cflag_(pset.getParameter<int>("cFlag")),
00040 comenergy(pset.getParameter<double>("comEnergy")),
00041 doquench_(pset.getParameter<bool>("doQuench")),
00042 doradiativeenloss_(pset.getParameter<bool>("doRadiativeEnLoss")),
00043 docollisionalenloss_(pset.getParameter<bool>("doCollisionalEnLoss")),
00044 nquarkflavor_(pset.getParameter<int>("numQuarkFlavor")),
00045 qgpt0_(pset.getParameter<double>("qgpInitialTemperature")),
00046 qgptau0_(pset.getParameter<double>("qgpProperTimeFormation")),
00047 maxEventsToPrint_(pset.getUntrackedParameter<int>("maxEventsToPrint",1)),
00048 pythiaHepMCVerbosity_(pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false)),
00049 pythiaPylistVerbosity_(pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0)),
00050 eventNumber_(0)
00051 {
00052   // Default constructor
00053 
00054   // Verbosity Level
00055   // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
00056   pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0);
00057   LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_ << endl;
00058   
00059   // HepMC event verbosity Level
00060   pythiaHepMCVerbosity_ = pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false);
00061   LogDebug("HepMCverbosity")  << "Pythia HepMC verbosity = " << pythiaHepMCVerbosity_ << endl; 
00062 
00063   //Max number of events printed on verbosity level 
00064   maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint",0);
00065   LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_ << endl;
00066 
00067   //initialize pythia
00068   pyqpythia_init(pset);
00069 
00070   //initilize pyquen
00071   pyquen_init(pset);
00072 
00073   // Call PYTHIA
00074   call_pyinit("CMS", "p", "p", comenergy);  
00075   
00076   cout<<endl;
00077 
00078   produces<HepMCProduct>();
00079   produces<GenInfoProduct, edm::InRun>();
00080 }
00081 
00082 
00083 //_____________________________________________________________________
00084 PyquenProducer::~PyquenProducer()
00085 {
00086   // distructor
00087 
00088   call_pystat(1);
00089 
00090   clear();
00091 }
00092 
00093 
00094 //_____________________________________________________________________
00095 void PyquenProducer::add_heavy_ion_rec(HepMC::GenEvent *evt)
00096 {
00097   HepMC::HeavyIon *hi = new HepMC::HeavyIon(
00098     -1,                                 // Ncoll_hard
00099     -1,                                 // Npart_proj
00100     -1,                                 // Npart_targ
00101     -1,                                 // Ncoll
00102     -1,                                 // spectator_neutrons
00103     -1,                                 // spectator_protons
00104     -1,                                 // N_Nwounded_collisions
00105     -1,                                 // Nwounded_N_collisions
00106     -1,                                 // Nwounded_Nwounded_collisions
00107     plfpar.bgen,                        // impact_parameter in [fm]
00108     0,                                  // event_plane_angle
00109     0,                                  // eccentricity
00110     -1                                  // sigma_inel_NN
00111   );
00112 
00113   evt->set_heavy_ion(*hi);
00114   delete hi;
00115 }
00116 
00117 
00118 //______________________________________________________________________
00119 bool PyquenProducer::call_pygive(const std::string& iParm ) 
00120 {
00121   // Set Pythia parameters
00122 
00123   int numWarn = pydat1.mstu[26];//# warnings
00124   int numErr = pydat1.mstu[22]; //# errors
00125   // call the fortran routine pygive with a fortran string
00126   PYGIVE( iParm.c_str(), iParm.length() );  
00127 
00128   // if an error or warning happens it is problem
00129   return pydat1.mstu[26] == numWarn && pydat1.mstu[22] == numErr;   
00130 }
00131 
00132 
00133 //____________________________________________________________________
00134 void PyquenProducer::clear()
00135 {
00136 }
00137 
00138 
00139 //_____________________________________________________________________
00140 void PyquenProducer::produce(Event & e, const EventSetup& es)
00141 {
00142   edm::LogInfo("PYQUENabeamtarget") << "##### PYQUEN: beam/target A = "                     << abeamtarget_;
00143   edm::LogInfo("PYQUENcflag")       << "##### PYQUEN: centrality flag cflag_ = "            << cflag_;
00144   edm::LogInfo("PYQUENbfixed")      << "##### PYQUEN: fixed impact parameter bFixed = "     << bfixed_;
00145   edm::LogInfo("PYQUENinNFlav")     << "##### PYQUEN: No active quark flavor nf = "         << pyqpar.nfu;
00146   edm::LogInfo("PYQUENinTemp")      << "##### PYQUEN: Initial temperature of QGP, T0 = "    << pyqpar.T0u;
00147   edm::LogInfo("PYQUENinTau")       << "##### PYQUEN: Proper formation time of QGP, tau0 =" << pyqpar.tau0u;
00148 
00149   // Generate PYQUEN event
00150   // generate single partonic PYTHIA jet event
00151   call_pyevnt();
00152 
00153   // call PYQUEN to apply parton rescattering and energy loss 
00154   // if doQuench=FALSE, it is pure PYTHIA
00155   if( doquench_ ){
00156     PYQUEN(abeamtarget_,cflag_,bfixed_);
00157     edm::LogInfo("PYQUENinAction") << "##### Calling PYQUEN("<<abeamtarget_<<","<<cflag_<<","<<bfixed_<<") ####";
00158   } else {
00159     edm::LogInfo("PYQUENinAction") << "##### Calling PYQUEN: QUENCHING OFF!! This is just PYTHIA !!!! ####";
00160   }
00161 
00162   // call PYTHIA to finish the hadronization
00163   PYEXEC();
00164 
00165   // fill the HEPEVT with the PYJETS event record
00166   call_pyhepc(1);
00167 
00168   // event information
00169   HepMC::GenEvent* evt = hepevtio2.read_next_event();
00170   evt->set_signal_process_id(pypars.msti[0]);      // type of the process
00171   evt->set_event_scale(pypars.pari[16]);           // Q^2
00172   ++eventNumber_;
00173   evt->set_event_number(eventNumber_);
00174 
00175   add_heavy_ion_rec(evt);
00176 
00177   auto_ptr<HepMCProduct> bare_product(new HepMCProduct());
00178   bare_product->addHepMCData(evt );
00179   e.put(bare_product); 
00180 
00181   // verbosity
00182   if( e.id().event() <= maxEventsToPrint_ && ( pythiaPylistVerbosity_ || pythiaHepMCVerbosity_ )) { 
00183     // Prints PYLIST info
00184      if( pythiaPylistVerbosity_ ){
00185        call_pylist(pythiaPylistVerbosity_);
00186      }
00187       
00188      // Prints HepMC event
00189      if( pythiaHepMCVerbosity_ ){
00190         cout << "Event process = " << pypars.msti[0] << endl; 
00191         evt->print(); 
00192      }
00193   }
00194     
00195   return;
00196 }
00197 
00198 
00199 //_____________________________________________________________________
00200 bool PyquenProducer::pyqpythia_init(const ParameterSet & pset)
00201 {
00202   //initialize PYTHIA
00203 
00204   //random number seed
00205   edm::Service<RandomNumberGenerator> rng;
00206   randomEngine = fRandomEngine = &(rng->getEngine());
00207   uint32_t seed = rng->mySeed();
00208   ostringstream sRandomSet;
00209   sRandomSet << "MRPY(1)=" << seed;
00210   call_pygive(sRandomSet.str());
00211 
00212   //Turn Hadronization Off if there is quenching
00213   if(doquench_){
00214      string sHadOff("MSTP(111)=0");
00215      call_pygive(sHadOff);
00216   }
00217 
00218     // Set PYTHIA parameters in a single ParameterSet  
00219   ParameterSet pythia_params = pset.getParameter<ParameterSet>("PythiaParameters") ;
00220   
00221   // The parameter sets to be read
00222   vector<string> setNames = pythia_params.getParameter<vector<string> >("parameterSets");
00223 
00224     // Loop over the sets
00225   for ( unsigned i=0; i<setNames.size(); ++i ) {
00226     string mySet = setNames[i];
00227     
00228     // Read the PYTHIA parameters for each set of parameters
00229     vector<string> pars = pythia_params.getParameter<vector<string> >(mySet);
00230     
00231     cout << "----------------------------------------------" << endl;
00232     cout << "Read PYTHIA parameter set " << mySet << endl;
00233     cout << "----------------------------------------------" << endl;
00234     
00235     // Loop over all parameters and stop in case of mistake
00236     for( vector<string>::const_iterator itPar = pars.begin(); itPar != pars.end(); ++itPar ) {
00237       static string sRandomValueSetting("MRPY(1)");
00238       if( 0 == itPar->compare(0,sRandomValueSetting.size(),sRandomValueSetting) ) {
00239          throw edm::Exception(edm::errors::Configuration,"PythiaError")
00240            << " Attempted to set random number using 'MRPY(1)'. NOT ALLOWED!\n"
00241               " Use RandomNumberGeneratorService to set the random number seed.";
00242       }
00243       if( !call_pygive(*itPar) ) {
00244         throw edm::Exception(edm::errors::Configuration,"PythiaError") 
00245            << "PYTHIA did not accept \""<<*itPar<<"\"";
00246       }
00247     }
00248   }
00249 
00250   return true;
00251 }
00252 
00253 
00254 //_________________________________________________________________
00255 bool PyquenProducer::pyquen_init(const ParameterSet &pset)
00256 {
00257   // PYQUEN initialization
00258 
00259   // angular emitted gluon  spectrum selection 
00260   pyqpar.ianglu = angularspecselector_;
00261 
00262   // type of medium induced partonic energy loss
00263   if( doradiativeenloss_ && docollisionalenloss_ ){
00264     edm::LogInfo("PYQUENinEnLoss") << "##### PYQUEN: Radiative AND Collisional partonic energy loss ON ####";
00265     pyqpar.ienglu = 0; 
00266   } else if ( doradiativeenloss_ ) {
00267     edm::LogInfo("PYQUENinRad") << "##### PYQUEN: Only RADIATIVE partonic energy loss ON ####";
00268     pyqpar.ienglu = 1; 
00269   } else if ( docollisionalenloss_ ) {
00270     edm::LogInfo("PYQUENinColl") << "##### PYQUEN: Only COLLISIONAL partonic energy loss ON ####";
00271     pyqpar.ienglu = 2; 
00272   } else {
00273     edm::LogInfo("PYQUENinEnLoss") << "##### PYQUEN: Radiative AND Collisional partonic energy loss ON ####";
00274     pyqpar.ienglu = 0; 
00275   }
00276 
00277   // number of active quark flavors in qgp
00278   pyqpar.nfu    = nquarkflavor_;
00279 
00280   // initial temperature of QGP
00281   pyqpar.T0u    = qgpt0_;
00282 
00283   // proper time of QGP formation
00284   pyqpar.tau0u  = qgptau0_;
00285 
00286   return true;
00287 }
00288 
00289 
00290 //____________________________________________________________________