PyquenInterface/src/PyquenHadronizer.cc

0001 #include <iostream>
0002 #include <ctime>
0003
0004 #include "GeneratorInterface/PyquenInterface/interface/PyquenHadronizer.h"
0005
0006 #include "GeneratorInterface/Core/interface/FortranInstance.h"
0007 #include "GeneratorInterface/PyquenInterface/interface/PyquenWrapper.h"
0008 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Declarations.h"
0009 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
0010
0011 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0012
0013 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0014 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include "GeneratorInterface/HiGenCommon/interface/HiGenEvtSelectorFactory.h"
0019
0020 #include "HepMC/IO_HEPEVT.h"
0021 #include "HepMC/PythiaWrapper.h"
0022
0023 using namespace gen;
0024 using namespace edm;
0025 using namespace std;
0026
0027 HepMC::IO_HEPEVT pyquen_hepevtio;
0028
0029 const std::vector<std::string> PyquenHadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6,
0030                                                                        gen::FortranInstance::kFortranInstance};
0031
0032 PyquenHadronizer ::PyquenHadronizer(const ParameterSet& pset, edm::ConsumesCollector&& iC)
0033     : BaseHadronizer(pset),
0034       pset_(pset),
0035       abeamtarget_(pset.getParameter<double>("aBeamTarget")),
0036       angularspecselector_(pset.getParameter<int>("angularSpectrumSelector")),
0037       bmin_(pset.getParameter<double>("bMin")),
0038       bmax_(pset.getParameter<double>("bMax")),
0039       bfixed_(pset.getParameter<double>("bFixed")),
0040       cflag_(pset.getParameter<int>("cFlag")),
0041       comenergy(pset.getParameter<double>("comEnergy")),
0042       doquench_(pset.getParameter<bool>("doQuench")),
0043       doradiativeenloss_(pset.getParameter<bool>("doRadiativeEnLoss")),
0044       docollisionalenloss_(pset.getParameter<bool>("doCollisionalEnLoss")),
0045       doIsospin_(pset.getParameter<bool>("doIsospin")),
0046       protonSide_(pset.getUntrackedParameter<int>("protonSide", 0)),
0047       embedding_(pset.getParameter<int>("embeddingMode")),
0048       evtPlane_(0),
0049       nquarkflavor_(pset.getParameter<int>("qgpNumQuarkFlavor")),
0050       qgpt0_(pset.getParameter<double>("qgpInitialTemperature")),
0051       qgptau0_(pset.getParameter<double>("qgpProperTimeFormation")),
0052       maxEventsToPrint_(pset.getUntrackedParameter<int>("maxEventsToPrint", 1)),
0053       fVertex_(nullptr),
0054       pythiaHepMCVerbosity_(pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false)),
0055       pythiaPylistVerbosity_(pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
0056       pythia6Service_(new Pythia6Service(pset)),
0057       filterType_(pset.getUntrackedParameter<string>("filterType", "None")) {
0058   if (pset.exists("signalVtx"))
0059     signalVtx_ = pset.getUntrackedParameter<std::vector<double> >("signalVtx");
0060
0061   if (signalVtx_.size() == 4) {
0062     if (!fVertex_)
0063       fVertex_ = new HepMC::FourVector();
0064     LogDebug("EventSignalVertex") << "Setting event signal vertex "
0065                                   << " x = " << signalVtx_.at(0) << " y = " << signalVtx_.at(1)
0066                                   << "  z= " << signalVtx_.at(2) << " t = " << signalVtx_.at(3) << endl;
0067     fVertex_->set(signalVtx_.at(0), signalVtx_.at(1), signalVtx_.at(2), signalVtx_.at(3));
0068   }
0069
0070   // Verbosity Level
0071   // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
0072   LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_ << endl;
0073   // HepMC event verbosity Level
0074   pythiaHepMCVerbosity_ = pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false);
0075   LogDebug("HepMCverbosity") << "Pythia HepMC verbosity = " << pythiaHepMCVerbosity_ << endl;
0076
0077   //Max number of events printed on verbosity level
0078   maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint", 0);
0079   LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_ << endl;
0080
0081   if (embedding_ == 1) {
0082     cflag_ = 0;
0083     src_ = iC.consumes<CrossingFrame<edm::HepMCProduct> >(
0084         pset.getUntrackedParameter<edm::InputTag>("backgroundLabel", edm::InputTag("mix", "generatorSmeared")));
0085   }
0086   selector_ = HiGenEvtSelectorFactory::get(filterType_, pset);
0087
0088   int cm = 1, va, vb, vc;
0089   PYQVER(cm, va, vb, vc);
0090   //HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
0091 }
0092
0093 //_____________________________________________________________________
0094 PyquenHadronizer::~PyquenHadronizer() {
0095   // distructor
0096   call_pystat(1);
0097
0098   delete pythia6Service_;
0099 }
0100
0101 //_____________________________________________________________________
0102 void PyquenHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { pythia6Service_->setRandomEngine(v); }
0103
0104 //_____________________________________________________________________
0105 void PyquenHadronizer::add_heavy_ion_rec(HepMC::GenEvent* evt) {
0106   HepMC::HeavyIon* hi = new HepMC::HeavyIon(1,            // Ncoll_hard
0107                                             -1,           // Npart_proj
0108                                             -1,           // Npart_targ
0109                                             1,            // Ncoll
0110                                             -1,           // spectator_neutrons
0111                                             -1,           // spectator_protons
0112                                             -1,           // N_Nwounded_collisions
0113                                             -1,           // Nwounded_N_collisions
0114                                             -1,           // Nwounded_Nwounded_collisions
0115                                             plfpar.bgen,  // impact_parameter in [fm]
0116                                             evtPlane_,    // event_plane_angle
0117                                             0,            // eccentricity
0118                                             -1            // sigma_inel_NN
0119   );
0120
0121   evt->set_heavy_ion(*hi);
0122
0123   delete hi;
0124 }
0125
0126 //_____________________________________________________________________
0127 bool PyquenHadronizer::generatePartonsAndHadronize() {
0128   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0129
0130   // Not possible to retrieve impact paramter and event plane info
0131   // at this part, need to overwrite filter() in
0132   // PyquenGeneratorFilter
0133
0134   if (embedding_ == 1) {
0135     const edm::Event& e = getEDMEvent();
0136     HepMC::GenVertex* genvtx = nullptr;
0137     const HepMC::GenEvent* inev = nullptr;
0138     Handle<CrossingFrame<HepMCProduct> > cf;
0139     e.getByToken(src_, cf);
0140     MixCollection<HepMCProduct> mix(cf.product());
0141     if (mix.size() < 1) {
0142       throw cms::Exception("MatchVtx") << "Mixing has " << mix.size() << " sub-events, should have been at least 1"
0143                                        << endl;
0144     }
0145     const HepMCProduct& bkg = mix.getObject(0);
0146     if (!(bkg.isVtxGenApplied())) {
0147       throw cms::Exception("MatchVtx") << "Input background does not have smeared vertex!" << endl;
0148     } else {
0149       inev = bkg.GetEvent();
0150     }
0151
0152     genvtx = inev->signal_process_vertex();
0153
0154     if (!genvtx)
0155       throw cms::Exception("MatchVtx") << "Input background does not have signal process vertex!" << endl;
0156
0157     double aX, aY, aZ, aT;
0158
0159     aX = genvtx->position().x();
0160     aY = genvtx->position().y();
0161     aZ = genvtx->position().z();
0162     aT = genvtx->position().t();
0163
0164     if (!fVertex_) {
0165       fVertex_ = new HepMC::FourVector();
0166     }
0167
0168     //LogInfo("MatchVtx")
0169     std::cout << " setting vertex "
0170               << " aX " << aX << " aY " << aY << " aZ " << aZ << " aT " << aT << endl;
0171     fVertex_->set(aX, aY, aZ, aT);
0172
0173     const HepMC::HeavyIon* hi = inev->heavy_ion();
0174
0175     if (hi) {
0176       bfixed_ = hi->impact_parameter();
0177       evtPlane_ = hi->event_plane_angle();
0178     } else {
0179       LogWarning("EventEmbedding") << "Background event does not have heavy ion record!";
0180     }
0181   }
0182
0183   // Generate PYQUEN event
0184   // generate single partonic PYTHIA jet event
0185
0186   // Take into account whether it's a nn or pp or pn interaction
0187   if (doIsospin_) {
0188     string projN = "p";
0189     string targN = "p";
0190     if (protonSide_ != 1)
0191       projN = nucleon();
0192     if (protonSide_ != 2)
0193       targN = nucleon();
0194     call_pyinit("CMS", projN.data(), targN.data(), comenergy);
0195   }
0196   call_pyevnt();
0197
0198   // call PYQUEN to apply parton rescattering and energy loss
0199   // if doQuench=FALSE, it is pure PYTHIA
0200   if (doquench_) {
0201     PYQUEN(abeamtarget_, cflag_, bfixed_, bmin_, bmax_);
0202     edm::LogInfo("PYQUENinAction") << "##### Calling PYQUEN(" << abeamtarget_ << "," << cflag_ << "," << bfixed_
0203                                    << ") ####";
0204   } else {
0205     edm::LogInfo("PYQUENinAction") << "##### Calling PYQUEN: QUENCHING OFF!! This is just PYTHIA !!!! ####";
0206   }
0207
0208   // call PYTHIA to finish the hadronization
0209   pyexec_();
0210
0211   // fill the HEPEVT with the PYJETS event record
0212   call_pyhepc(1);
0213
0214   // event information
0215   // pyquen_hepevtio.set_trust_mothers_before_daughters(true);
0216   HepMC::GenEvent* evt = pyquen_hepevtio.read_next_event();
0217
0218   // signal vertex
0219   HepMC::GenVertex* sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0);  // just initialization
0220   if (!evt->signal_process_vertex())
0221     evt->set_signal_process_vertex(sub_vertices);
0222
0223   delete sub_vertices;
0224   evt->set_signal_process_id(pypars.msti[0]);  // type of the process
0225   evt->set_event_scale(pypars.pari[16]);       // Q^2
0226
0227   if (embedding_)
0228     rotateEvtPlane(evt, evtPlane_);
0229   add_heavy_ion_rec(evt);
0230
0231   if (fVertex_) {
0232     // Copy the HepMC::GenEvent
0233     std::unique_ptr<edm::HepMCProduct> HepMCEvt(new edm::HepMCProduct(evt));
0234
0235     HepMCEvt->applyVtxGen(fVertex_);
0236     evt = new HepMC::GenEvent((*HepMCEvt->GetEvent()));
0237   }
0238
0239   //  HepMC::HEPEVT_Wrapper::check_hepevt_consistency();
0240
0241   event().reset(evt);
0242
0243   return true;
0244 }
0245
0246 bool PyquenHadronizer::readSettings(int) {
0247   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0248   pythia6Service_->setGeneralParams();
0249   pythia6Service_->setCSAParams();
0250
0251   //Proton to Nucleon fraction
0252   pfrac_ = 1. / (1.98 + 0.015 * pow(abeamtarget_, 2. / 3));
0253
0254   //initialize pythia
0255   pyqpythia_init(pset_);
0256
0257   //initilize pyquen
0258   pyquen_init(pset_);
0259
0260   return true;
0261 }
0262
0263 bool PyquenHadronizer::initializeForInternalPartons() {
0264   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0265
0266   // Call PYTHIA
0267   call_pyinit("CMS", "p", "p", comenergy);
0268
0269   return true;
0270 }
0271
0272 //_____________________________________________________________________
0273 bool PyquenHadronizer::pyqpythia_init(const ParameterSet& pset) {
0274   //Turn Hadronization Off whether or not there is quenching
0275   // PYEXEC is called later anyway
0276   string sHadOff("MSTP(111)=0");
0277   gen::call_pygive(sHadOff);
0278
0279   return true;
0280 }
0281
0282 //_________________________________________________________________
0283 bool PyquenHadronizer::pyquen_init(const ParameterSet& pset) {
0284   // PYQUEN initialization
0285
0286   // angular emitted gluon  spectrum selection
0287   pyqpar.ianglu = angularspecselector_;
0288
0289   // type of medium induced partonic energy loss
0290   if (doradiativeenloss_ && docollisionalenloss_) {
0291     edm::LogInfo("PYQUENinEnLoss") << "##### PYQUEN: Radiative AND Collisional partonic energy loss ON ####";
0292     pyqpar.ienglu = 0;
0293   } else if (doradiativeenloss_) {
0294     edm::LogInfo("PYQUENinRad") << "##### PYQUEN: Only RADIATIVE partonic energy loss ON ####";
0295     pyqpar.ienglu = 1;
0296   } else if (docollisionalenloss_) {
0297     edm::LogInfo("PYQUENinColl") << "##### PYQUEN: Only COLLISIONAL partonic energy loss ON ####";
0298     pyqpar.ienglu = 2;
0299   } else {
0300     edm::LogInfo("PYQUENinEnLoss") << "##### PYQUEN: Radiative AND Collisional partonic energy loss ON ####";
0301     pyqpar.ienglu = 0;
0302   }
0303
0304   // number of active quark flavors in qgp
0305   pyqpar.nfu = nquarkflavor_;
0306   // initial temperature of QGP
0307   pyqpar.T0u = qgpt0_;
0308   // proper time of QGP formation
0309   pyqpar.tau0u = qgptau0_;
0310
0311   return true;
0312 }
0313
0314 const char* PyquenHadronizer::nucleon() {
0315   int* dummy = nullptr;
0316   double random = gen::pyr_(dummy);
0317   const char* nuc = nullptr;
0318   if (random > pfrac_)
0319     nuc = "n";
0320   else
0321     nuc = "p";
0322
0323   return nuc;
0324 }
0325
0326 void PyquenHadronizer::rotateEvtPlane(HepMC::GenEvent* evt, double angle) {
0327   double sinphi0 = sin(angle);
0328   double cosphi0 = cos(angle);
0329
0330   for (HepMC::GenEvent::vertex_iterator vt = evt->vertices_begin(); vt != evt->vertices_end(); ++vt) {
0331     double x0 = (*vt)->position().x();
0332     double y0 = (*vt)->position().y();
0333     double z = (*vt)->position().z();
0334     double t = (*vt)->position().t();
0335
0336     double x = x0 * cosphi0 - y0 * sinphi0;
0337     double y = y0 * cosphi0 + x0 * sinphi0;
0338
0339     (*vt)->set_position(HepMC::FourVector(x, y, z, t));
0340   }
0341
0342   for (HepMC::GenEvent::particle_iterator vt = evt->particles_begin(); vt != evt->particles_end(); ++vt) {
0343     double x0 = (*vt)->momentum().x();
0344     double y0 = (*vt)->momentum().y();
0345     double z = (*vt)->momentum().z();
0346     double t = (*vt)->momentum().t();
0347
0348     double x = x0 * cosphi0 - y0 * sinphi0;
0349     double y = y0 * cosphi0 + x0 * sinphi0;
0350
0351     (*vt)->set_momentum(HepMC::FourVector(x, y, z, t));
0352   }
0353 }
0354
0355 bool PyquenHadronizer::declareStableParticles(const std::vector<int>& _pdg) {
0356   std::vector<int> pdg = _pdg;
0357   for (size_t i = 0; i < pdg.size(); i++) {
0358     int pyCode = pycomp_(pdg[i]);
0359     std::ostringstream pyCard;
0360     pyCard << "MDCY(" << pyCode << ",1)=0";
0361     std::cout << pyCard.str() << std::endl;
0362     call_pygive(pyCard.str());
0363   }
0364
0365   return true;
0366 }
0367
0368 //____________________________________________________________________
0369
0370 bool PyquenHadronizer::hadronize() { return false; }
0371
0372 bool PyquenHadronizer::decay() { return true; }
0373
0374 bool PyquenHadronizer::residualDecay() { return true; }
0375
0376 void PyquenHadronizer::finalizeEvent() {}
0377
0378 void PyquenHadronizer::statistics() {}
0379
0380 const char* PyquenHadronizer::classname() const { return "gen::PyquenHadronizer"; }