12 #include "boost/lexical_cast.hpp"
28 #include "HepMC/PythiaWrapper6_4.h"
29 #include "HepMC/GenEvent.h"
30 #include "HepMC/HeavyIon.h"
31 #include "HepMC/SimpleVector.h"
58 abeamtarget_(pset.getParameter<double>(
"aBeamTarget")),
59 bfixed_(pset.getParameter<double>(
"bFixed")),
60 bmax_(pset.getParameter<double>(
"bMax")),
61 bmin_(pset.getParameter<double>(
"bMin")),
62 cflag_(pset.getParameter<int>(
"cFlag")),
63 embedding_(pset.getParameter<bool>(
"embeddingMode")),
64 comenergy(pset.getParameter<double>(
"comEnergy")),
65 doradiativeenloss_(pset.getParameter<bool>(
"doRadiativeEnLoss")),
66 docollisionalenloss_(pset.getParameter<bool>(
"doCollisionalEnLoss")),
67 fracsoftmult_(pset.getParameter<double>(
"fracSoftMultiplicity")),
68 hadfreeztemp_(pset.getParameter<double>(
"hadronFreezoutTemperature")),
69 hymode_(pset.getParameter<
string>(
"hydjetMode")),
70 maxEventsToPrint_(pset.getUntrackedParameter<int>(
"maxEventsToPrint", 1)),
71 maxlongy_(pset.getParameter<double>(
"maxLongitudinalRapidity")),
72 maxtrany_(pset.getParameter<double>(
"maxTransverseRapidity")),
75 nmultiplicity_(pset.getParameter<int>(
"nMultiplicity")),
77 nquarkflavor_(pset.getParameter<int>(
"qgpNumQuarkFlavor")),
78 pythiaPylistVerbosity_(pset.getUntrackedParameter<int>(
"pythiaPylistVerbosity", 0)),
79 qgpt0_(pset.getParameter<double>(
"qgpInitialTemperature")),
80 qgptau0_(pset.getParameter<double>(
"qgpProperTimeFormation")),
84 rotate_(pset.getParameter<bool>(
"rotateEventPlane")),
85 shadowingswitch_(pset.getParameter<int>(
"shadowingSwitch")),
86 signn_(pset.getParameter<double>(
"sigmaInelNN")),
119 int nproj =
static_cast<int>(npart / 2);
120 int ntarg =
static_cast<int>(npart - nproj);
122 HepMC::HeavyIon* hi =
new HepMC::HeavyIon(
126 static_cast<int>(
hyfpar.nbcol),
138 evt->set_heavy_ion(*hi);
163 p->suggest_barcode(barcode);
179 HepMC::GenVertex* vertex =
new HepMC::GenVertex(HepMC::FourVector(x,y,z,t),
id);
195 const HepMC::GenEvent * inev = input->GetEvent();
196 const HepMC::HeavyIon* hi = inev->heavy_ion();
198 bfixed_ = hi->impact_parameter();
199 phi0_ = hi->event_plane_angle();
203 LogWarning(
"EventEmbedding")<<
"Background event does not have heavy ion record!";
213 edm::LogInfo(
"HYDJETinTemp") <<
"##### HYDJET: QGP init temperature, T0 ="<<
pyqpar.T0u;
214 edm::LogInfo(
"HYDJETinTau") <<
"##### HYDJET: QGP formation time,tau0 ="<<
pyqpar.tau0u;
220 edm::LogError(
"HydjetEmptyEvent") <<
"##### HYDJET: No Particles generated, Number of tries ="<<ntry;
225 std::ostringstream sstr;
226 sstr <<
"HydjetHadronizerProducer: No particles generated after " << ntry <<
" tries.\n";
241 HepMC::GenEvent *
evt =
new HepMC::GenEvent();
245 evt->set_signal_process_id(
pypars.msti[0]);
246 evt->set_event_scale(
pypars.pari[16]);
270 vector<HepMC::GenVertex*> sub_vertices(
nsub_);
273 for(
int isub=0;isub<
nsub_;isub++){
274 LogDebug(
"SubEvent") <<
"Sub Event ID : "<<isub;
276 int sub_up = (isub+1)*50000;
277 vector<HepMC::GenParticle*> particles;
278 vector<int> mother_ids;
279 vector<HepMC::GenVertex*> prods;
281 sub_vertices[isub] =
new HepMC::GenVertex(HepMC::FourVector(0,0,0,0),isub);
282 evt->add_vertex(sub_vertices[isub]);
283 if(!evt->signal_process_vertex()) evt->set_signal_process_vertex(sub_vertices[isub]);
288 mother_ids.push_back(
hyjets.khj[2][ihy]);
297 LogDebug(
"Hydjet")<<
"Number of particles in vector "<<particles.size();
299 for (
unsigned int i = 0;
i<particles.size();
i++) {
304 int mid = mother_ids[
i]-isub*50000-1;
306 LogDebug(
"DecayChain")<<
"Mother's ID "<<mid;
307 LogDebug(
"DecayChain")<<
"Particle's PDG ID "<<part->pdg_id();
310 sub_vertices[isub]->add_particle_out(part);
316 LogDebug(
"DecayChain")<<
"Mother's PDG ID "<<mother->pdg_id();
318 HepMC::GenVertex* prod_vertex = mother->end_vertex();
320 prod_vertex = prods[
i];
321 prod_vertex->add_particle_in(mother);
322 evt->add_vertex(prod_vertex);
325 prod_vertex->add_particle_out(part);
329 for (
unsigned int i = 0;
i<prods.size();
i++) {
330 if(prods[
i])
delete prods[
i];
339 double bmax,
double bfix,
int nh)
344 HYINIT(energy,a,ifb,bmin,bmax,bfix,nh);
397 edm::LogInfo(
"HydjetEnLoss") <<
"##### Radiative AND Collisional partonic energy loss ON ####";
400 edm::LogInfo(
"HydjetenLoss") <<
"##### Only RADIATIVE partonic energy loss ON ####";
403 edm::LogInfo(
"HydjetEnLoss") <<
"##### Only COLLISIONAL partonic energy loss ON ####";
406 edm::LogInfo(
"HydjetEnLoss") <<
"##### Radiative AND Collisional partonic energy loss ON ####";
432 LogInfo(
"RAScaling")<<
"Nuclear radius(RA) = "<<ra;
448 std::vector<int> pdg = _pdg;
449 for (
size_t i=0;
i < pdg.size();
i++ ) {
451 std::ostringstream pyCard ;
452 pyCard <<
"MDCY(" << pyCode <<
",1)=0";
462 const double pi = 3.14159265358979;
495 return "gen::HydjetHadronizer";
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
bool docollisionalenloss_
DEFAULT = true.
bool hydjet_init(const edm::ParameterSet &pset)
Sin< T >::type sin(const T &t)
bool call_pygive(const std::string &line)
void add_heavy_ion_rec(HepMC::GenEvent *evt)
std::auto_ptr< HepMC::GenEvent > & event()
bool generatePartonsAndHadronize()
Cos< T >::type cos(const T &t)
const char * classname() const
virtual ~HydjetHadronizer()
unsigned int maxEventsToPrint_
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
HepMC::GenParticle * build_hyjet(int index, int barcode)
Pythia6Service * pythia6Service_
double fracsoftmult_
DEFAULT = true.
HepMC::GenVertex * build_hyjet_vertex(int i, int id)
double nuclear_radius() const
bool call_hyinit(double energy, double a, int ifb, double bmin, double bmax, double bfix, int nh)
unsigned int pythiaPylistVerbosity_
bool initializeForInternalPartons()
unsigned int shadowingswitch_
unsigned int nquarkflavor_
bool declareStableParticles(const std::vector< int > &)
bool get_particles(HepMC::GenEvent *evt)
edm::Event & getEDMEvent() const