Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-02 04:16:31

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     isub = std::floor((hyjets.khj[2][ihy] / 50000));
0345     int hjoffset = isub * 50000;
0346 
0347     if (isub != isub_l) {
0348       sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), isub);
0349       evt->add_vertex(sub_vertices);
0350       if (!evt->signal_process_vertex())
0351         evt->set_signal_process_vertex(sub_vertices);
0352       index[isub] = ihy - 1;
0353       isub_l = isub;
0354     }
0355 
0356     if (convertStatus(hyjets.khj[0][ihy]) == 1)
0357       stab++;
0358     LogDebug("Hydjet_array") << ihy << " MULTin ev.:" << hyjets.nhj << " SubEv.#" << isub << " Part #" << ihy + 1
0359                              << ", PDG: " << hyjets.khj[1][ihy] << " (st. " << convertStatus(hyjets.khj[0][ihy])
0360                              << ") mother=" << hyjets.khj[2][ihy] - (isub * 50000) + index[isub] + 1 << " ("
0361                              << hyjets.khj[2][ihy] << "), childs ("
0362                              << hyjets.khj[3][ihy] - (isub * 50000) + index[isub] + 1 << "-"
0363                              << hyjets.khj[4][ihy] - (isub * 50000) + index[isub] + 1 << "), vtx ("
0364                              << hyjets.vhj[0][ihy] << "," << hyjets.vhj[1][ihy] << "," << hyjets.vhj[2][ihy] << ") "
0365                              << std::endl;
0366 
0367     if (hyjets.khj[2][ihy] == 0) {
0368       primary_particle[ihy] = build_hyjet(ihy, ihy + 1);
0369       if (!sub_vertices)
0370         LogError("Hydjet_array") << "##### HYDJET2: Vertex not initialized!";
0371       else
0372         sub_vertices->add_particle_out(primary_particle[ihy]);
0373       LogDebug("Hydjet_array") << " ---> " << ihy + 1 << std::endl;
0374     } else {
0375       particle[ihy] = build_hyjet(ihy, ihy + 1);
0376       int mid = hyjets.khj[2][ihy] - hjoffset + index[isub];
0377       int mid_t = mid;
0378       while (((mid + 1) < ihy) && (std::abs(hyjets.khj[1][mid]) < 100) &&
0379              (hyjets.khj[3][mid + 1] - hjoffset + index[isub] <= ihy))
0380         mid++;
0381       if (std::abs(hyjets.khj[1][mid]) < 100)
0382         mid = mid_t;
0383 
0384       HepMC::GenParticle* mother = primary_particle.at(mid);
0385       HepMC::GenVertex* prods = build_hyjet_vertex(ihy, isub);
0386 
0387       if (!mother) {
0388         mother = particle[mid];
0389         primary_particle[mid] = mother;
0390       }
0391 
0392       HepMC::GenVertex* prod_vertex = mother->end_vertex();
0393       if (!prod_vertex) {
0394         prod_vertex = prods;
0395         prod_vertex->add_particle_in(mother);
0396         LogDebug("Hydjet_array") << " <--- " << mid + 1 << std::endl;
0397         evt->add_vertex(prod_vertex);
0398         prods = nullptr;
0399       }
0400 
0401       prod_vertex->add_particle_out(particle[ihy]);
0402       LogDebug("Hydjet_array") << " ---" << mid + 1 << "---> " << ihy + 1 << std::endl;
0403       delete prods;
0404     }
0405     ihy++;
0406   }
0407   LogDebug("Hydjet_array") << " MULTin ev.:" << hyjets.nhj << ", last index: " << ihy - 1
0408                            << ", Sub events: " << isub + 1 << ", stable particles: " << stab << std::endl;
0409 
0410   return true;
0411 }
0412 
0413 //______________________________________________________________
0414 bool HydjetHadronizer::call_hyinit(double energy, double a, int ifb, double bmin, double bmax, double bfix, int nh) {
0415   // initialize hydjet
0416 
0417   pydatr.mrpy[2] = 1;
0418   HYINIT(energy, a, ifb, bmin, bmax, bfix, nh);
0419   return true;
0420 }
0421 
0422 //______________________________________________________________
0423 bool HydjetHadronizer::hydjet_init(const ParameterSet& pset) {
0424   // set hydjet options
0425 
0426   // hydjet running mode mode
0427   // kHydroOnly --- nhsel=0 jet production off (pure HYDRO event), nhsel=0
0428   // kHydroJets --- nhsle=1 jet production on, jet quenching off (HYDRO+njet*PYTHIA events)
0429   // kHydroQJet --- nhsel=2 jet production & jet quenching on (HYDRO+njet*PYQUEN events)
0430   // kJetsOnly  --- nhsel=3 jet production on, jet quenching off, HYDRO off (njet*PYTHIA events)
0431   // kQJetsOnly --- nhsel=4 jet production & jet quenching on, HYDRO off (njet*PYQUEN events)
0432 
0433   if (hymode_ == "kHydroOnly")
0434     hyjpar.nhsel = 0;
0435   else if (hymode_ == "kHydroJets")
0436     hyjpar.nhsel = 1;
0437   else if (hymode_ == "kHydroQJets")
0438     hyjpar.nhsel = 2;
0439   else if (hymode_ == "kJetsOnly")
0440     hyjpar.nhsel = 3;
0441   else if (hymode_ == "kQJetsOnly")
0442     hyjpar.nhsel = 4;
0443   else
0444     hyjpar.nhsel = 2;
0445 
0446   // fraction of soft hydro induced multiplicity
0447   hyflow.fpart = fracsoftmult_;
0448 
0449   // hadron freez-out temperature
0450   hyflow.Tf = hadfreeztemp_;
0451 
0452   // maximum longitudinal collective rapidity
0453   hyflow.ylfl = maxlongy_;
0454 
0455   // maximum transverse collective rapidity
0456   hyflow.ytfl = maxtrany_;
0457 
0458   // shadowing on=1, off=0
0459   hyjpar.ishad = shadowingswitch_;
0460 
0461   // set inelastic nucleon-nucleon cross section
0462   hyjpar.sigin = signn_;
0463 
0464   // angular emitted gluon spectrum selection
0465   pyqpar.ianglu = angularspecselector_;
0466 
0467   // number of active quark flavors in qgp
0468   pyqpar.nfu = nquarkflavor_;
0469 
0470   // initial temperature of QGP
0471   pyqpar.T0u = qgpt0_;
0472 
0473   // proper time of QGP formation
0474   pyqpar.tau0u = qgptau0_;
0475 
0476   // type of medium induced partonic energy loss
0477   if (doradiativeenloss_ && docollisionalenloss_) {
0478     edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
0479     pyqpar.ienglu = 0;
0480   } else if (doradiativeenloss_) {
0481     edm::LogInfo("HydjetenLoss") << "##### Only RADIATIVE partonic energy loss ON ####";
0482     pyqpar.ienglu = 1;
0483   } else if (docollisionalenloss_) {
0484     edm::LogInfo("HydjetEnLoss") << "##### Only COLLISIONAL partonic energy loss ON ####";
0485     pyqpar.ienglu = 2;
0486   } else {
0487     edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
0488     pyqpar.ienglu = 0;
0489   }
0490   return true;
0491 }
0492 
0493 //_____________________________________________________________________
0494 
0495 bool HydjetHadronizer::readSettings(int) {
0496   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0497   pythia6Service_->setGeneralParams();
0498 
0499   return true;
0500 }
0501 
0502 //_____________________________________________________________________
0503 
0504 bool HydjetHadronizer::initializeForInternalPartons() {
0505   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0506   // pythia6Service_->setGeneralParams();
0507 
0508   // the input impact parameter (bxx_) is in [fm]; transform in [fm/RA] for hydjet usage
0509   const float ra = nuclear_radius();
0510   LogInfo("RAScaling") << "Nuclear radius(RA) =  " << ra;
0511   bmin_ /= ra;
0512   bmax_ /= ra;
0513   bfixed_ /= ra;
0514 
0515   // hydjet running options
0516   hydjet_init(pset_);
0517   // initialize hydjet
0518   LogInfo("HYDJETinAction") << "##### Calling HYINIT(" << comenergy << "," << abeamtarget_ << "," << cflag_ << ","
0519                             << bmin_ << "," << bmax_ << "," << bfixed_ << "," << nmultiplicity_ << ") ####";
0520   call_hyinit(comenergy, abeamtarget_, cflag_, bmin_, bmax_, bfixed_, nmultiplicity_);
0521   return true;
0522 }
0523 
0524 bool HydjetHadronizer::declareStableParticles(const std::vector<int>& _pdg) {
0525   std::vector<int> pdg = _pdg;
0526   for (size_t i = 0; i < pdg.size(); i++) {
0527     int pyCode = pycomp_(pdg[i]);
0528     std::ostringstream pyCard;
0529     pyCard << "MDCY(" << pyCode << ",1)=0";
0530     std::cout << pyCard.str() << std::endl;
0531     call_pygive(pyCard.str());
0532   }
0533   return true;
0534 }
0535 
0536 //________________________________________________________________
0537 void HydjetHadronizer::rotateEvtPlane() {
0538   const double pi = 3.14159265358979;
0539   phi0_ = 2. * pi * gen::pyr_(nullptr) - pi;
0540   sinphi0_ = sin(phi0_);
0541   cosphi0_ = cos(phi0_);
0542 }
0543 
0544 //________________________________________________________________
0545 bool HydjetHadronizer::hadronize() { return false; }
0546 
0547 bool HydjetHadronizer::decay() { return true; }
0548 
0549 bool HydjetHadronizer::residualDecay() { return true; }
0550 
0551 void HydjetHadronizer::finalizeEvent() {}
0552 
0553 void HydjetHadronizer::statistics() {}
0554 
0555 const char* HydjetHadronizer::classname() const { return "gen::HydjetHadronizer"; }