HydjetInterface/src/HydjetHadronizer.cc

0001 /**
0002    \brief Interface to the HYDJET generator (since core v. 1.9.1), produces HepMC events
0003    \version 2.2
0004    \authors Camelia Mironov, Andrey Belyaev
0005 */
0006
0007 #include <iostream>
0008 #include <cmath>
0009
0010 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/Run.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Utilities/interface/EDMException.h"
0016
0017 #include "GeneratorInterface/Core/interface/FortranInstance.h"
0018 #include "GeneratorInterface/HydjetInterface/interface/HydjetHadronizer.h"
0019 #include "GeneratorInterface/HydjetInterface/interface/HydjetWrapper.h"
0020 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Declarations.h"
0021 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
0022
0023 #include "HepMC/IO_HEPEVT.h"
0024 #include "HepMC/PythiaWrapper6_4.h"
0025 #include "HepMC/GenEvent.h"
0026 #include "HepMC/HeavyIon.h"
0027 #include "HepMC/SimpleVector.h"
0028
0029 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0030 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0031 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0032 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
0033
0034 using namespace edm;
0035 using namespace std;
0036 using namespace gen;
0037
0038 int HydjetHadronizer::convertStatus(int st) {
0039   if (st <= 0)
0040     return 0;
0041   if (st <= 10)
0042     return 1;
0043   if (st <= 20)
0044     return 2;
0045   if (st <= 30)
0046     return 3;
0047   else
0048     return st;
0049 }
0050
0051 const std::vector<std::string> HydjetHadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6,
0052                                                                        gen::FortranInstance::kFortranInstance};
0053
0054 //_____________________________________________________________________
0055 HydjetHadronizer::HydjetHadronizer(const ParameterSet& pset, edm::ConsumesCollector&& iC)
0056     : BaseHadronizer(pset),
0057       evt(nullptr),
0058       pset_(pset),
0059       abeamtarget_(pset.getParameter<double>("aBeamTarget")),
0060       angularspecselector_(pset.getParameter<int>("angularSpectrumSelector")),
0061       bfixed_(pset.getParameter<double>("bFixed")),
0062       bmax_(pset.getParameter<double>("bMax")),
0063       bmin_(pset.getParameter<double>("bMin")),
0064       cflag_(pset.getParameter<int>("cFlag")),
0065       embedding_(pset.getParameter<int>("embeddingMode")),
0066       comenergy(pset.getParameter<double>("comEnergy")),
0067       doradiativeenloss_(pset.getParameter<bool>("doRadiativeEnLoss")),
0068       docollisionalenloss_(pset.getParameter<bool>("doCollisionalEnLoss")),
0069       fracsoftmult_(pset.getParameter<double>("fracSoftMultiplicity")),
0070       hadfreeztemp_(pset.getParameter<double>("hadronFreezoutTemperature")),
0071       hymode_(pset.getParameter<string>("hydjetMode")),
0072       maxEventsToPrint_(pset.getUntrackedParameter<int>("maxEventsToPrint", 1)),
0073       maxlongy_(pset.getParameter<double>("maxLongitudinalRapidity")),
0074       maxtrany_(pset.getParameter<double>("maxTransverseRapidity")),
0075       nsub_(0),
0076       nhard_(0),
0077       nmultiplicity_(pset.getParameter<int>("nMultiplicity")),
0078       nsoft_(0),
0079       nquarkflavor_(pset.getParameter<int>("qgpNumQuarkFlavor")),
0080       pythiaPylistVerbosity_(pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
0081       qgpt0_(pset.getParameter<double>("qgpInitialTemperature")),
0082       qgptau0_(pset.getParameter<double>("qgpProperTimeFormation")),
0083       phi0_(0.),
0084       sinphi0_(0.),
0085       cosphi0_(1.),
0086       rotate_(pset.getParameter<bool>("rotateEventPlane")),
0087       shadowingswitch_(pset.getParameter<int>("shadowingSwitch")),
0088       signn_(pset.getParameter<double>("sigmaInelNN")),
0089       fVertex_(nullptr),
0090       pythia6Service_(new Pythia6Service(pset)) {
0091   // Default constructor
0092
0093   if (pset.exists("signalVtx"))
0094     signalVtx_ = pset.getUntrackedParameter<std::vector<double> >("signalVtx");
0095
0096   if (signalVtx_.size() == 4) {
0097     if (!fVertex_)
0098       fVertex_ = new HepMC::FourVector();
0099     LogDebug("EventSignalVertex") << "Setting event signal vertex "
0100                                   << " x = " << signalVtx_.at(0) << " y = " << signalVtx_.at(1)
0101                                   << "  z= " << signalVtx_.at(2) << " t = " << signalVtx_.at(3) << endl;
0102     fVertex_->set(signalVtx_.at(0), signalVtx_.at(1), signalVtx_.at(2), signalVtx_.at(3));
0103   }
0104
0105   // PYLIST Verbosity Level
0106   // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
0107   pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0);
0108   LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_;
0109
0110   //Max number of events printed on verbosity level
0111   maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint", 0);
0112   LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_;
0113
0114   if (embedding_) {
0115     cflag_ = 0;
0116     src_ = iC.consumes<CrossingFrame<edm::HepMCProduct> >(
0117         pset.getUntrackedParameter<edm::InputTag>("backgroundLabel", edm::InputTag("mix", "generatorSmeared")));
0118   }
0119
0120   int cm = 1, va, vb, vc;
0121   HYJVER(cm, va, vb, vc);
0122   HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
0123 }
0124
0125 //_____________________________________________________________________
0126 HydjetHadronizer::~HydjetHadronizer() {
0127   // destructor
0128   call_pystat(1);
0129   delete pythia6Service_;
0130 }
0131
0132 //_____________________________________________________________________
0133 void HydjetHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { pythia6Service_->setRandomEngine(v); }
0134
0135 //_____________________________________________________________________
0136 void HydjetHadronizer::add_heavy_ion_rec(HepMC::GenEvent* evt) {
0137   // heavy ion record in the final CMSSW Event
0138   double npart = hyfpar.npart;
0139   int nproj = static_cast<int>(npart / 2);
0140   int ntarg = static_cast<int>(npart - nproj);
0141
0142   HepMC::HeavyIon* hi = new HepMC::HeavyIon(nsub_,                           // Ncoll_hard/N of SubEvents
0143                                             nproj,                           // Npart_proj
0144                                             ntarg,                           // Npart_targ
0145                                             static_cast<int>(hyfpar.nbcol),  // Ncoll
0146                                             0,                               // spectator_neutrons
0147                                             0,                               // spectator_protons
0148                                             0,                               // N_Nwounded_collisions
0149                                             0,                               // Nwounded_N_collisions
0150                                             0,                               // Nwounded_Nwounded_collisions
0151                                             hyfpar.bgen * nuclear_radius(),  // impact_parameter in [fm]
0152                                             phi0_,                           // event_plane_angle
0153                                             0,  //hypsi3.psi3,                                   // eccentricity
0154                                             hyjpar.sigin  // sigma_inel_NN
0155   );
0156
0157   evt->set_heavy_ion(*hi);
0158   delete hi;
0159 }
0160
0161 //___________________________________________________________________
0162 HepMC::GenParticle* HydjetHadronizer::build_hyjet(int index, int barcode) {
0163   // Build particle object corresponding to index in hyjets (soft+hard)
0164   double x0 = hyjets.phj[0][index];
0165   double y0 = hyjets.phj[1][index];
0166
0167   double x = x0 * cosphi0_ - y0 * sinphi0_;
0168   double y = y0 * cosphi0_ + x0 * sinphi0_;
0169
0170   HepMC::GenParticle* p = new HepMC::GenParticle(HepMC::FourVector(x,                      // px
0171                                                                    y,                      // py
0172                                                                    hyjets.phj[2][index],   // pz
0173                                                                    hyjets.phj[3][index]),  // E
0174                                                  hyjets.khj[1][index],                     // id
0175                                                  convertStatus(hyjets.khj[0][index]        // status
0176                                                                ));
0177
0178   p->suggest_barcode(barcode);
0179   return p;
0180 }
0181
0182 //___________________________________________________________________
0183 HepMC::GenVertex* HydjetHadronizer::build_hyjet_vertex(int i, int id) {
0184   // build verteces for the hyjets stored events
0185   double x0 = hyjets.vhj[0][i];
0186   double y0 = hyjets.vhj[1][i];
0187   double x = x0 * cosphi0_ - y0 * sinphi0_;
0188   double y = y0 * cosphi0_ + x0 * sinphi0_;
0189   double z = hyjets.vhj[2][i];
0190   double t = hyjets.vhj[4][i];
0191
0192   HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t), id);
0193   return vertex;
0194 }
0195
0196 //___________________________________________________________________
0197
0198 bool HydjetHadronizer::generatePartonsAndHadronize() {
0199   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0200
0201   // generate single event
0202   if (embedding_) {
0203     const edm::Event& e = getEDMEvent();
0204     HepMC::GenVertex* genvtx = nullptr;
0205     const HepMC::GenEvent* inev = nullptr;
0206     Handle<CrossingFrame<HepMCProduct> > cf;
0207     e.getByToken(src_, cf);
0208     MixCollection<HepMCProduct> mix(cf.product());
0209     if (mix.size() < 1) {
0210       throw cms::Exception("MatchVtx") << "Mixing has " << mix.size() << " sub-events, should have been at least 1"
0211                                        << endl;
0212     }
0213     const HepMCProduct& bkg = mix.getObject(0);
0214     if (!(bkg.isVtxGenApplied())) {
0215       throw cms::Exception("MatchVtx") << "Input background does not have smeared vertex!" << endl;
0216     } else {
0217       inev = bkg.GetEvent();
0218     }
0219
0220     genvtx = inev->signal_process_vertex();
0221
0222     if (!genvtx)
0223       throw cms::Exception("MatchVtx") << "Input background does not have signal process vertex!" << endl;
0224
0225     double aX, aY, aZ, aT;
0226
0227     aX = genvtx->position().x();
0228     aY = genvtx->position().y();
0229     aZ = genvtx->position().z();
0230     aT = genvtx->position().t();
0231
0232     if (!fVertex_) {
0233       fVertex_ = new HepMC::FourVector();
0234     }
0235     LogInfo("MatchVtx") << " setting vertex "
0236                         << " aX " << aX << " aY " << aY << " aZ " << aZ << " aT " << aT << endl;
0237     fVertex_->set(aX, aY, aZ, aT);
0238
0239     const HepMC::HeavyIon* hi = inev->heavy_ion();
0240
0241     if (hi) {
0242       bfixed_ = (hi->impact_parameter()) / nuclear_radius();
0243       phi0_ = hi->event_plane_angle();
0244       sinphi0_ = sin(phi0_);
0245       cosphi0_ = cos(phi0_);
0246     } else {
0247       LogWarning("EventEmbedding") << "Background event does not have heavy ion record!";
0248     }
0249
0250   } else if (rotate_)
0251     rotateEvtPlane();
0252
0253   nsoft_ = 0;
0254   nhard_ = 0;
0255
0256   edm::LogInfo("HYDJETmode") << "##### HYDJET  nhsel = " << hyjpar.nhsel;
0257   edm::LogInfo("HYDJETfpart") << "##### HYDJET fpart = " << hyflow.fpart;
0258   edm::LogInfo("HYDJETtf") << "##### HYDJET hadron freez-out temp, Tf = " << hyflow.Tf;
0259   edm::LogInfo("HYDJETinTemp") << "##### HYDJET: QGP init temperature, T0 =" << pyqpar.T0u;
0260   edm::LogInfo("HYDJETinTau") << "##### HYDJET: QGP formation time,tau0 =" << pyqpar.tau0u;
0261
0262   int ntry = 0;
0263   while (nsoft_ == 0 && nhard_ == 0) {
0264     if (ntry > 100) {
0265       edm::LogError("HydjetEmptyEvent") << "##### HYDJET: No Particles generated, Number of tries =" << ntry;
0266
0267       // Throw an exception.  Use the EventCorruption exception since it maps onto SkipEvent
0268       // which is what we want to do here.
0269
0270       std::ostringstream sstr;
0271       sstr << "HydjetHadronizerProducer: No particles generated after " << ntry << " tries.\n";
0272       edm::Exception except(edm::errors::EventCorruption, sstr.str());
0273       throw except;
0274     } else {
0275       HYEVNT(bfixed_);
0276       nsoft_ = hyfpar.nhyd;
0277       nsub_ = hyjpar.njet;
0278       nhard_ = hyfpar.npyt;
0279       ++ntry;
0280     }
0281   }
0282
0283   if (hyjpar.nhsel < 3)
0284     nsub_++;
0285
0286   // event information
0287   std::unique_ptr<HepMC::GenEvent> evt = std::make_unique<HepMC::GenEvent>();
0288   std::unique_ptr<edm::HepMCProduct> HepMCEvt = std::make_unique<edm::HepMCProduct>();
0289
0290   if (nhard_ > 0 || nsoft_ > 0)
0291     get_particles(evt.get());
0292
0293   evt->set_signal_process_id(pypars.msti[0]);  // type of the process
0294   evt->set_event_scale(pypars.pari[16]);       // Q^2
0295   add_heavy_ion_rec(evt.get());
0296
0297   if (fVertex_) {
0298     // generate new vertex & apply the shift
0299     // Copy the HepMC::GenEvent
0300     HepMCEvt = std::make_unique<edm::HepMCProduct>(evt.get());
0301     HepMCEvt->applyVtxGen(fVertex_);
0302     evt = std::make_unique<HepMC::GenEvent>(*HepMCEvt->GetEvent());
0303   }
0304
0305   HepMC::HEPEVT_Wrapper::check_hepevt_consistency();
0306   LogDebug("HEPEVT_info") << "Ev numb: " << HepMC::HEPEVT_Wrapper::event_number()
0307                           << " Entries number: " << HepMC::HEPEVT_Wrapper::number_entries() << " Max. entries "
0308                           << HepMC::HEPEVT_Wrapper::max_number_entries() << std::endl;
0309
0310   event() = std::move(evt);
0311   return true;
0312 }
0313
0314 //_____________________________________________________________________
0315 bool HydjetHadronizer::get_particles(HepMC::GenEvent* evt) {
0316   // Hard particles. The first nhard_ lines from hyjets array.
0317   // Pythia/Pyquen sub-events (sub-collisions) for a given event
0318   // Return T/F if success/failure
0319   // Create particles from lujet entries, assign them into vertices and
0320   // put the vertices in the GenEvent, for each SubEvent
0321   // The SubEvent information is kept by storing indeces of main vertices
0322   // of subevents as a vector in GenHIEvent.
0323
0324   LogDebug("Hydjet") << " Number of sub events " << nsub_;
0325   LogDebug("Hydjet") << " Number of hard events " << hyjpar.njet;
0326   LogDebug("Hydjet") << " Number of hard particles " << nhard_;
0327   LogDebug("Hydjet") << " Number of soft particles " << nsoft_;
0328   LogDebug("Hydjet") << " nhard_ + nsoft_ = " << nhard_ + nsoft_ << " hyjets.nhj = " << hyjets.nhj << endl;
0329
0330   int ihy = 0;
0331   int isub = -1;
0332   int isub_l = -1;
0333   int stab = 0;
0334
0335   vector<HepMC::GenParticle*> primary_particle(hyjets.nhj);
0336   vector<HepMC::GenParticle*> particle(hyjets.nhj);
0337
0338   HepMC::GenVertex* sub_vertices = nullptr;  // just initialization
0339
0340   // contain the last index in for each subevent
0341   vector<int> index(nsub_);
0342
0343   while (ihy < hyjets.nhj) {
0344     constexpr int kMaxMultiplicity = 200000;
0345     isub = std::floor((hyjets.khj[2][ihy] / kMaxMultiplicity));
0346     int hjoffset = isub * kMaxMultiplicity;
0347
0348     if (isub != isub_l) {
0349       sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), isub);
0350       evt->add_vertex(sub_vertices);
0351       if (!evt->signal_process_vertex())
0352         evt->set_signal_process_vertex(sub_vertices);
0353       if (isub >= static_cast<int>(index.size())) {
0354         index.resize(isub + 1);
0355       }
0356       index[isub] = ihy - 1;
0357       isub_l = isub;
0358     }
0359
0360     if (convertStatus(hyjets.khj[0][ihy]) == 1)
0361       stab++;
0362     LogDebug("Hydjet_array") << ihy << " MULTin ev.:" << hyjets.nhj << " SubEv.#" << isub << " Part #" << ihy + 1
0363                              << ", PDG: " << hyjets.khj[1][ihy] << " (st. " << convertStatus(hyjets.khj[0][ihy])
0364                              << ") mother=" << hyjets.khj[2][ihy] - (isub * kMaxMultiplicity) + index[isub] + 1 << " ("
0365                              << hyjets.khj[2][ihy] << "), childs ("
0366                              << hyjets.khj[3][ihy] - (isub * kMaxMultiplicity) + index[isub] + 1 << "-"
0367                              << hyjets.khj[4][ihy] - (isub * kMaxMultiplicity) + index[isub] + 1 << "), vtx ("
0368                              << hyjets.vhj[0][ihy] << "," << hyjets.vhj[1][ihy] << "," << hyjets.vhj[2][ihy] << ") "
0369                              << std::endl;
0370
0371     if (hyjets.khj[2][ihy] == 0) {
0372       primary_particle[ihy] = build_hyjet(ihy, ihy + 1);
0373       if (!sub_vertices)
0374         LogError("Hydjet_array") << "##### HYDJET2: Vertex not initialized!";
0375       else
0376         sub_vertices->add_particle_out(primary_particle[ihy]);
0377       LogDebug("Hydjet_array") << " ---> " << ihy + 1 << std::endl;
0378     } else {
0379       particle[ihy] = build_hyjet(ihy, ihy + 1);
0380       int mid = hyjets.khj[2][ihy] - hjoffset + index[isub];
0381       if (mid > ihy) {
0382         throw cms::Exception("BadHydjetParent") << "Parent ID " << mid << " is greater than child id " << ihy
0383                                                 << " this is not allowed.\n When this happened in the past it was "
0384                                                    "caused the sub event multiplicity being > "
0385                                                 << kMaxMultiplicity << " which caused an offset overflow.";
0386       }
0387       int mid_t = mid;
0388       while (((mid + 1) < ihy) && (std::abs(hyjets.khj[1][mid]) < 100) &&
0389              (hyjets.khj[3][mid + 1] - hjoffset + index[isub] <= ihy))
0390         mid++;
0391       if (std::abs(hyjets.khj[1][mid]) < 100)
0392         mid = mid_t;
0393
0394       HepMC::GenParticle* mother = primary_particle.at(mid);
0395       HepMC::GenVertex* prods = build_hyjet_vertex(ihy, isub);
0396
0397       if (!mother) {
0398         mother = particle[mid];
0399         primary_particle[mid] = mother;
0400       }
0401
0402       HepMC::GenVertex* prod_vertex = mother->end_vertex();
0403       if (!prod_vertex) {
0404         prod_vertex = prods;
0405         prod_vertex->add_particle_in(mother);
0406         LogDebug("Hydjet_array") << " <--- " << mid + 1 << std::endl;
0407         evt->add_vertex(prod_vertex);
0408         prods = nullptr;
0409       }
0410
0411       prod_vertex->add_particle_out(particle[ihy]);
0412       LogDebug("Hydjet_array") << " ---" << mid + 1 << "---> " << ihy + 1 << std::endl;
0413       delete prods;
0414     }
0415     ihy++;
0416   }
0417   LogDebug("Hydjet_array") << " MULTin ev.:" << hyjets.nhj << ", last index: " << ihy - 1
0418                            << ", Sub events: " << isub + 1 << ", stable particles: " << stab << std::endl;
0419
0420   return true;
0421 }
0422
0423 //______________________________________________________________
0424 bool HydjetHadronizer::call_hyinit(double energy, double a, int ifb, double bmin, double bmax, double bfix, int nh) {
0425   // initialize hydjet
0426
0427   pydatr.mrpy[2] = 1;
0428   HYINIT(energy, a, ifb, bmin, bmax, bfix, nh);
0429   return true;
0430 }
0431
0432 //______________________________________________________________
0433 bool HydjetHadronizer::hydjet_init(const ParameterSet& pset) {
0434   // set hydjet options
0435
0436   // hydjet running mode mode
0437   // kHydroOnly --- nhsel=0 jet production off (pure HYDRO event), nhsel=0
0438   // kHydroJets --- nhsle=1 jet production on, jet quenching off (HYDRO+njet*PYTHIA events)
0439   // kHydroQJet --- nhsel=2 jet production & jet quenching on (HYDRO+njet*PYQUEN events)
0440   // kJetsOnly  --- nhsel=3 jet production on, jet quenching off, HYDRO off (njet*PYTHIA events)
0441   // kQJetsOnly --- nhsel=4 jet production & jet quenching on, HYDRO off (njet*PYQUEN events)
0442
0443   if (hymode_ == "kHydroOnly")
0444     hyjpar.nhsel = 0;
0445   else if (hymode_ == "kHydroJets")
0446     hyjpar.nhsel = 1;
0447   else if (hymode_ == "kHydroQJets")
0448     hyjpar.nhsel = 2;
0449   else if (hymode_ == "kJetsOnly")
0450     hyjpar.nhsel = 3;
0451   else if (hymode_ == "kQJetsOnly")
0452     hyjpar.nhsel = 4;
0453   else
0454     hyjpar.nhsel = 2;
0455
0456   // fraction of soft hydro induced multiplicity
0457   hyflow.fpart = fracsoftmult_;
0458
0459   // hadron freez-out temperature
0460   hyflow.Tf = hadfreeztemp_;
0461
0462   // maximum longitudinal collective rapidity
0463   hyflow.ylfl = maxlongy_;
0464
0465   // maximum transverse collective rapidity
0466   hyflow.ytfl = maxtrany_;
0467
0468   // shadowing on=1, off=0
0469   hyjpar.ishad = shadowingswitch_;
0470
0471   // set inelastic nucleon-nucleon cross section
0472   hyjpar.sigin = signn_;
0473
0474   // angular emitted gluon spectrum selection
0475   pyqpar.ianglu = angularspecselector_;
0476
0477   // number of active quark flavors in qgp
0478   pyqpar.nfu = nquarkflavor_;
0479
0480   // initial temperature of QGP
0481   pyqpar.T0u = qgpt0_;
0482
0483   // proper time of QGP formation
0484   pyqpar.tau0u = qgptau0_;
0485
0486   // type of medium induced partonic energy loss
0487   if (doradiativeenloss_ && docollisionalenloss_) {
0488     edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
0489     pyqpar.ienglu = 0;
0490   } else if (doradiativeenloss_) {
0491     edm::LogInfo("HydjetenLoss") << "##### Only RADIATIVE partonic energy loss ON ####";
0492     pyqpar.ienglu = 1;
0493   } else if (docollisionalenloss_) {
0494     edm::LogInfo("HydjetEnLoss") << "##### Only COLLISIONAL partonic energy loss ON ####";
0495     pyqpar.ienglu = 2;
0496   } else {
0497     edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
0498     pyqpar.ienglu = 0;
0499   }
0500   return true;
0501 }
0502
0503 //_____________________________________________________________________
0504
0505 bool HydjetHadronizer::readSettings(int) {
0506   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0507   pythia6Service_->setGeneralParams();
0508
0509   return true;
0510 }
0511
0512 //_____________________________________________________________________
0513
0514 bool HydjetHadronizer::initializeForInternalPartons() {
0515   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0516   // pythia6Service_->setGeneralParams();
0517
0518   // the input impact parameter (bxx_) is in [fm]; transform in [fm/RA] for hydjet usage
0519   const float ra = nuclear_radius();
0520   LogInfo("RAScaling") << "Nuclear radius(RA) =  " << ra;
0521   bmin_ /= ra;
0522   bmax_ /= ra;
0523   bfixed_ /= ra;
0524
0525   // hydjet running options
0526   hydjet_init(pset_);
0527   // initialize hydjet
0528   LogInfo("HYDJETinAction") << "##### Calling HYINIT(" << comenergy << "," << abeamtarget_ << "," << cflag_ << ","
0529                             << bmin_ << "," << bmax_ << "," << bfixed_ << "," << nmultiplicity_ << ") ####";
0530   call_hyinit(comenergy, abeamtarget_, cflag_, bmin_, bmax_, bfixed_, nmultiplicity_);
0531   return true;
0532 }
0533
0534 bool HydjetHadronizer::declareStableParticles(const std::vector<int>& _pdg) {
0535   std::vector<int> pdg = _pdg;
0536   for (size_t i = 0; i < pdg.size(); i++) {
0537     int pyCode = pycomp_(pdg[i]);
0538     std::ostringstream pyCard;
0539     pyCard << "MDCY(" << pyCode << ",1)=0";
0540     std::cout << pyCard.str() << std::endl;
0541     call_pygive(pyCard.str());
0542   }
0543   return true;
0544 }
0545
0546 //________________________________________________________________
0547 void HydjetHadronizer::rotateEvtPlane() {
0548   const double pi = 3.14159265358979;
0549   phi0_ = 2. * pi * gen::pyr_(nullptr) - pi;
0550   sinphi0_ = sin(phi0_);
0551   cosphi0_ = cos(phi0_);
0552 }
0553
0554 //________________________________________________________________
0555 bool HydjetHadronizer::hadronize() { return false; }
0556
0557 bool HydjetHadronizer::decay() { return true; }
0558
0559 bool HydjetHadronizer::residualDecay() { return true; }
0560
0561 void HydjetHadronizer::finalizeEvent() {}
0562
0563 void HydjetHadronizer::statistics() {}
0564
0565 const char* HydjetHadronizer::classname() const { return "gen::HydjetHadronizer"; }