23 #include "HepMC/IO_HEPEVT.h"
24 #include "HepMC/PythiaWrapper6_4.h"
25 #include "HepMC/GenEvent.h"
26 #include "HepMC/HeavyIon.h"
27 #include "HepMC/SimpleVector.h"
39 int convertStatus(
int st) {
61 abeamtarget_(pset.getParameter<double>(
"aBeamTarget")),
62 angularspecselector_(pset.getParameter<int>(
"angularSpectrumSelector")),
63 bfixed_(pset.getParameter<double>(
"bFixed")),
64 bmax_(pset.getParameter<double>(
"bMax")),
65 bmin_(pset.getParameter<double>(
"bMin")),
66 cflag_(pset.getParameter<int>(
"cFlag")),
67 embedding_(pset.getParameter<bool>(
"embeddingMode")),
68 comenergy(pset.getParameter<double>(
"comEnergy")),
69 doradiativeenloss_(pset.getParameter<bool>(
"doRadiativeEnLoss")),
70 docollisionalenloss_(pset.getParameter<bool>(
"doCollisionalEnLoss")),
71 fracsoftmult_(pset.getParameter<double>(
"fracSoftMultiplicity")),
72 hadfreeztemp_(pset.getParameter<double>(
"hadronFreezoutTemperature")),
73 hymode_(pset.getParameter<
string>(
"hydjetMode")),
74 maxEventsToPrint_(pset.getUntrackedParameter<int>(
"maxEventsToPrint", 1)),
75 maxlongy_(pset.getParameter<double>(
"maxLongitudinalRapidity")),
76 maxtrany_(pset.getParameter<double>(
"maxTransverseRapidity")),
79 nmultiplicity_(pset.getParameter<int>(
"nMultiplicity")),
81 nquarkflavor_(pset.getParameter<int>(
"qgpNumQuarkFlavor")),
82 pythiaPylistVerbosity_(pset.getUntrackedParameter<int>(
"pythiaPylistVerbosity", 0)),
83 qgpt0_(pset.getParameter<double>(
"qgpInitialTemperature")),
84 qgptau0_(pset.getParameter<double>(
"qgpProperTimeFormation")),
88 rotate_(pset.getParameter<bool>(
"rotateEventPlane")),
89 shadowingswitch_(pset.getParameter<int>(
"shadowingSwitch")),
90 signn_(pset.getParameter<double>(
"sigmaInelNN")),
95 if (pset.
exists(
"signalVtx"))
101 LogDebug(
"EventSignalVertex") <<
"Setting event signal vertex "
122 int cm = 1, va, vb, vc;
124 HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
141 int nproj =
static_cast<int>(npart / 2);
142 int ntarg =
static_cast<int>(npart - nproj);
144 HepMC::HeavyIon* hi =
new HepMC::HeavyIon(
nsub_,
147 static_cast<int>(
hyfpar.nbcol),
159 evt->set_heavy_ion(*hi);
177 convertStatus(
hyjets.khj[0][index]
180 p->suggest_barcode(barcode);
194 HepMC::GenVertex* vertex =
new HepMC::GenVertex(HepMC::FourVector(x, y, z, t),
id);
206 HepMC::GenVertex* genvtx =
nullptr;
211 if (
mix.size() < 1) {
212 throw cms::Exception(
"MatchVtx") <<
"Mixing has " <<
mix.size() <<
" sub-events, should have been at least 1"
217 throw cms::Exception(
"MatchVtx") <<
"Input background does not have smeared vertex!" << endl;
222 genvtx = inev->signal_process_vertex();
225 throw cms::Exception(
"MatchVtx") <<
"Input background does not have signal process vertex!" << endl;
227 double aX, aY, aZ, aT;
229 aX = genvtx->position().x();
230 aY = genvtx->position().y();
231 aZ = genvtx->position().z();
232 aT = genvtx->position().t();
237 LogInfo(
"MatchVtx") <<
" setting vertex "
238 <<
" aX " << aX <<
" aY " << aY <<
" aZ " << aZ <<
" aT " << aT << endl;
241 const HepMC::HeavyIon* hi = inev->heavy_ion();
245 phi0_ = hi->event_plane_angle();
249 LogWarning(
"EventEmbedding") <<
"Background event does not have heavy ion record!";
261 edm::LogInfo(
"HYDJETinTemp") <<
"##### HYDJET: QGP init temperature, T0 =" <<
pyqpar.T0u;
262 edm::LogInfo(
"HYDJETinTau") <<
"##### HYDJET: QGP formation time,tau0 =" <<
pyqpar.tau0u;
267 edm::LogError(
"HydjetEmptyEvent") <<
"##### HYDJET: No Particles generated, Number of tries =" << ntry;
272 std::ostringstream sstr;
273 sstr <<
"HydjetHadronizerProducer: No particles generated after " << ntry <<
" tries.\n";
294 evt->set_signal_process_id(
pypars.msti[0]);
295 evt->set_event_scale(
pypars.pari[16]);
308 HepMC::HEPEVT_Wrapper::check_hepevt_consistency();
309 LogDebug(
"HEPEVT_info") <<
"Ev numb: " << HepMC::HEPEVT_Wrapper::event_number()
310 <<
" Entries number: " << HepMC::HEPEVT_Wrapper::number_entries() <<
" Max. entries "
311 << HepMC::HEPEVT_Wrapper::max_number_entries() << std::endl;
338 vector<HepMC::GenParticle*> primary_particle(
hyjets.nhj);
339 vector<HepMC::GenParticle*> particle(
hyjets.nhj);
341 HepMC::GenVertex* sub_vertices =
new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0);
346 while (ihy <
hyjets.nhj) {
347 isub = std::floor((
hyjets.khj[2][ihy] / 50000));
348 int hjoffset = isub * 50000;
350 if (isub != isub_l) {
351 sub_vertices =
new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), isub);
352 evt->add_vertex(sub_vertices);
353 if (!evt->signal_process_vertex())
354 evt->set_signal_process_vertex(sub_vertices);
355 index[isub] = ihy - 1;
359 if (convertStatus(
hyjets.khj[0][ihy]) == 1)
361 LogDebug(
"Hydjet_array") << ihy <<
" MULTin ev.:" <<
hyjets.nhj <<
" SubEv.#" << isub <<
" Part #" << ihy + 1
362 <<
", PDG: " <<
hyjets.khj[1][ihy] <<
" (st. " << convertStatus(
hyjets.khj[0][ihy])
363 <<
") mother=" <<
hyjets.khj[2][ihy] - (isub * 50000) + index[isub] + 1 <<
" ("
364 <<
hyjets.khj[2][ihy] <<
"), childs ("
365 <<
hyjets.khj[3][ihy] - (isub * 50000) + index[isub] + 1 <<
"-"
366 <<
hyjets.khj[4][ihy] - (isub * 50000) + index[isub] + 1 <<
"), vtx ("
367 <<
hyjets.vhj[0][ihy] <<
"," <<
hyjets.vhj[1][ihy] <<
"," <<
hyjets.vhj[2][ihy] <<
") "
372 sub_vertices->add_particle_out(primary_particle[ihy]);
373 LogDebug(
"Hydjet_array") <<
" ---> " << ihy + 1 << std::endl;
376 int mid =
hyjets.khj[2][ihy] - hjoffset + index[isub];
378 while ((mid < ihy) && (
hyjets.khj[1][mid] < 100) && (
hyjets.khj[3][mid + 1] - hjoffset + index[isub] == ihy))
380 if (
hyjets.khj[1][mid] < 100)
387 mother = particle[mid];
388 primary_particle[mid] = mother;
391 HepMC::GenVertex* prod_vertex = mother->end_vertex();
394 prod_vertex->add_particle_in(mother);
395 LogDebug(
"Hydjet_array") <<
" <--- " << mid + 1 << std::endl;
396 evt->add_vertex(prod_vertex);
400 prod_vertex->add_particle_out(particle[ihy]);
401 LogDebug(
"Hydjet_array") <<
" ---" << mid + 1 <<
"---> " << ihy + 1 << std::endl;
408 LogDebug(
"Hydjet_array") <<
" MULTin ev.:" <<
hyjets.nhj <<
", last index: " << ihy - 1
409 <<
", Sub events: " << isub + 1 <<
", stable particles: " << stab << std::endl;
419 HYINIT(energy, a, ifb, bmin, bmax, bfix, nh);
436 else if (
hymode_ ==
"kHydroJets")
438 else if (
hymode_ ==
"kHydroQJets")
440 else if (
hymode_ ==
"kJetsOnly")
442 else if (
hymode_ ==
"kQJetsOnly")
479 edm::LogInfo(
"HydjetEnLoss") <<
"##### Radiative AND Collisional partonic energy loss ON ####";
482 edm::LogInfo(
"HydjetenLoss") <<
"##### Only RADIATIVE partonic energy loss ON ####";
485 edm::LogInfo(
"HydjetEnLoss") <<
"##### Only COLLISIONAL partonic energy loss ON ####";
488 edm::LogInfo(
"HydjetEnLoss") <<
"##### Radiative AND Collisional partonic energy loss ON ####";
511 LogInfo(
"RAScaling") <<
"Nuclear radius(RA) = " << ra;
526 std::vector<int> pdg = _pdg;
527 for (
size_t i = 0;
i < pdg.size();
i++) {
529 std::ostringstream pyCard;
530 pyCard <<
"MDCY(" << pyCode <<
",1)=0";
539 const double pi = 3.14159265358979;
double bfixed_
fixed impact param (fm); valid only if cflag_=0
T getUntrackedParameter(std::string const &, T const &) const
void doSetRandomEngine(CLHEP::HepRandomEngine *v) override
bool docollisionalenloss_
DEFAULT = true.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
bool hydjet_init(const edm::ParameterSet &pset)
Sin< T >::type sin(const T &t)
bool call_pygive(const std::string &line)
bool exists(std::string const ¶meterName) const
checks if a parameter exists
void add_heavy_ion_rec(HepMC::GenEvent *evt)
Log< level::Error, false > LogError
double comenergy
collision energy
static const std::string kPythia6
bool generatePartonsAndHadronize()
bool embedding_
Switch for embedding mode.
std::string hymode_
Hydjet running mode.
int nsoft_
multiplicity of HYDRO-induced particles in event
if(conf_.getParameter< bool >("UseStripCablingDB"))
Cos< T >::type cos(const T &t)
const char * classname() const
unsigned int maxEventsToPrint_
Events to print if verbosity.
std::unique_ptr< HepMC::GenEvent > & event()
HepMC::GenParticle * build_hyjet(int index, int barcode)
Pythia6Service * pythia6Service_
~HydjetHadronizer() override
int nsub_
number of sub-events
std::vector< double > signalVtx_
Pset double vector to set event signal vertex.
Log< level::Info, false > LogInfo
static const std::string kFortranInstance
int nhard_
multiplicity of PYTHIA(+PYQUEN)-induced particles in event
HepMC::GenVertex * build_hyjet_vertex(int i, int id)
const HepMC::GenEvent * GetEvent() const
T const * product() const
bool isVtxGenApplied() const
double nuclear_radius() const
bool call_hyinit(double energy, double a, int ifb, double bmin, double bmax, double bfix, int nh)
unsigned int pythiaPylistVerbosity_
pythia verbosity; def=1
void setRandomEngine(CLHEP::HepRandomEngine *v)
double abeamtarget_
beam/target atomic mass number
bool initializeForInternalPartons()
edm::EDGetTokenT< CrossingFrame< edm::HepMCProduct > > src_
unsigned int shadowingswitch_
bool doradiativeenloss_
DEFAULT = true.
unsigned int nquarkflavor_
int angularspecselector_
angular emitted gluon spectrum selection
Log< level::Warning, false > LogWarning
bool declareStableParticles(const std::vector< int > &)
bool rotate_
Switch to rotate event plane.
double phi0_
Event plane angle.
HepMC::FourVector * fVertex_
Event signal vertex.
bool get_particles(HepMC::GenEvent *evt)
edm::Event & getEDMEvent() const