CMSSW Reference Manual

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