Skip to content

Commit

Permalink
Merge pull request #38443 from wouf/from-CMSSW_12_5_X_2022-06-20-1100
Browse files Browse the repository at this point in the history 
Hydjet2 bugfixes
  • Loading branch information
@cmsbuild
cmsbuild committed Aug 11, 2022
2 parents f32617f + fbba702 commit 1ad0fb5
Show file tree
Hide file tree
Showing 6 changed files with 509 additions and 93 deletions.
4 changes: 2 additions & 2 deletions GeneratorInterface/Hydjet2Interface/interface/Hydjet2Hadronizer.h
View file
Expand Up @@ -4,7 +4,7 @@
/**
* \class HydjetHadronizer
* \brief Interface to the HYDJET++ (Hydjet2) generator (since core v. 2.4.2), produces HepMC events
* \version 1.0
* \version 1.1
* \author Andrey Belyaev
*/

Expand Down Expand Up @@ -74,7 +74,7 @@ namespace gen {

inline double nuclear_radius() const;

int convertStatusForComponents(int, int);
int convertStatusForComponents(int, int, int);
int convertStatus(int);

InitialParamsHydjet_t fParams;
Expand Down
14 changes: 7 additions & 7 deletions GeneratorInterface/Hydjet2Interface/python/hydjet2DefaultParameters_cff.py
View file
Expand Up @@ -50,7 +50,7 @@
fR = cms.double(13.45), # Maximal transverse radius at thermal freeze-out for central collisions, fR [fm]
fSigmaTau = cms.double(3.5), # Duration of emission at thermal freeze-out for central collisions, fSigmaTau [fm/c]

fWeakDecay = cms.double(0.000000000000001), # Low decay width threshold fWeakDecay[GeV]: width<fWeakDecay decay off, width>=fDecayWidth decay on; can be used to switch off weak decays
fWeakDecay = cms.double(0.00000000000001), # Low decay width threshold fWeakDecay[GeV]: width<fWeakDecay decay off, width>=fDecayWidth decay on; can be used to switch off weak decays

)

Expand All @@ -66,19 +66,19 @@
fMu_th_pip = cms.double(0.), # Chemical potential of pi+ at thermal freeze-out, fMu_th_pip [GeV]

### Maximal longitudinal flow rapidity at thermal freeze-out ###
fYlmax = cms.double(4.5), # Maximal longitudinal flow rapidity at thermal freeze-out, fYlmax
fYlmax = cms.double(3.99), # Maximal longitudinal flow rapidity at thermal freeze-out, fYlmax


fUmax = cms.double(1.35), # Maximal transverse flow rapidity at thermal freeze-out for central collisions, fUmax
fPtmin = cms.double(10.), # Minimal pt of parton-parton scattering in PYTHIA event, fPtmin [GeV/c]
fUmax = cms.double(1.280), # Maximal transverse flow rapidity at thermal freeze-out for central collisions, fUmax
fPtmin = cms.double(9.06), # Minimal pt of parton-parton scattering in PYTHIA event, fPtmin [GeV/c]
fT0 = cms.double(1.1), # Initial QGP temperature for central Pb+Pb collisions in mid-rapidity, fT0 [GeV]; allowed range [0.2,2.0]GeV;

### Volume parameters at thermal freeze-out ###
fTau = cms.double(13.2), # Proper time proper at thermal freeze-out for central collisions, fTau [fm/c]
fR = cms.double(13.9), # Maximal transverse radius at thermal freeze-out for central collisions, fR [fm]
fTau = cms.double(11.5), # Proper time proper at thermal freeze-out for central collisions, fTau [fm/c]
fR = cms.double(16.), # Maximal transverse radius at thermal freeze-out for central collisions, fR [fm]
fSigmaTau = cms.double(2.), # Duration of emission at thermal freeze-out for central collisions, fSigmaTau [fm/c]

fWeakDecay = cms.double(0.000000000000001), # Low decay width threshold fWeakDecay[GeV]: width<fWeakDecay decay off, width>=fDecayWidth decay on; can be used to switch off weak decays
fWeakDecay = cms.double(0.00000000000001), # Low decay width threshold fWeakDecay[GeV]: width<fWeakDecay decay off, width>=fDecayWidth decay on; can be used to switch off weak decays
)

qgpParametersLHC = cms.PSet(
Expand Down
88 changes: 50 additions & 38 deletions GeneratorInterface/Hydjet2Interface/src/Hydjet2Hadronizer.cc
View file
@@ -1,6 +1,6 @@
/**
* \brief Interface to the HYDJET++ (Hydjet2) generator (since core v. 2.4.3), produces HepMC events
* \version 1.2
* \version 1.3
* \author Andrey Belyaev
*/

Expand Down Expand Up @@ -40,32 +40,40 @@ using namespace edm;
using namespace std;
using namespace gen;

int Hydjet2Hadronizer::convertStatusForComponents(int sta, int typ) {
if (sta == 1 && typ == 0)
return 6;
if (sta == 1 && typ == 1)
return 7;
if (sta == 2 && typ == 0)
return 16;
if (sta == 2 && typ == 1)
return 17;

else
return sta;
int Hydjet2Hadronizer::convertStatusForComponents(int sta, int typ, int pySta) {
int st = -1;
if (typ == 0) //soft
st = 2 - sta;
else if (typ == 1)
st = convertStatus(pySta);

if (st == -1)
throw cms::Exception("ConvertStatus") << "Wrong status code!" << endl;

if (separateHydjetComponents_) {
if (st == 1 && typ == 0)
return 6;
if (st == 1 && typ == 1)
return 7;
if (st == 2 && typ == 0)
return 16;
if (st == 2 && typ == 1)
return 17;
}
return st;
}

int Hydjet2Hadronizer::convertStatus(int st) {
if (!separateHydjetComponents_) {
if (st <= 0)
return 0;
if (st <= 10)
return 1;
if (st <= 20)
return 2;
if (st <= 30)
return 3;
}
return st;
if (st <= 0)
return 0;
if (st <= 10)
return 1;
if (st <= 20)
return 2;
if (st <= 30)
return 3;
else
return -1;
}

const std::vector<std::string> Hydjet2Hadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6};
Expand Down Expand Up @@ -408,7 +416,8 @@ bool Hydjet2Hadronizer::get_particles(HepMC::GenEvent *evt) {
isub_l = hj2->GetiJet().at(ihy);
}

if ((hj2->GetFinal().at(ihy)) == 1) //convertStatus(hj2->GetPythiaStatus().at(ihy)) == 1)
if ((convertStatusForComponents(
(hj2->GetFinal()).at(ihy), (hj2->GetType()).at(ihy), (hj2->GetPythiaStatus().at(ihy)))) == 1)
stab++;
LogDebug("Hydjet2_array") << ihy << " MULTin ev.:" << hj2->GetNtot() << " SubEv.#" << hj2->GetiJet().at(ihy)
<< " Part #" << ihy + 1 << ", PDG: " << hj2->GetPdg().at(ihy) << " (st. "
Expand Down Expand Up @@ -471,15 +480,15 @@ HepMC::GenParticle *Hydjet2Hadronizer::build_hyjet2(int index, int barcode) {
double px = px0 * cosphi0_ - py0 * sinphi0_;
double py = py0 * cosphi0_ + px0 * sinphi0_;

HepMC::GenParticle *p =
new HepMC::GenParticle(HepMC::FourVector(px, // px
py, // py
(hj2->GetPz()).at(index), // pz
(hj2->GetE()).at(index)), // E
(hj2->GetPdg()).at(index), // id
convertStatus(convertStatusForComponents((hj2->GetFinal()).at(index),
(hj2->GetType()).at(index))) // status
);
HepMC::GenParticle *p = new HepMC::GenParticle(
HepMC::FourVector(px, // px
py, // py
(hj2->GetPz()).at(index), // pz
(hj2->GetE()).at(index)), // E
(hj2->GetPdg()).at(index), // id
convertStatusForComponents(
(hj2->GetFinal()).at(index), (hj2->GetType()).at(index), (hj2->GetPythiaStatus()).at(index)) // status
);

p->suggest_barcode(barcode);
return p;
Expand All @@ -490,10 +499,13 @@ HepMC::GenVertex *Hydjet2Hadronizer::build_hyjet2_vertex(int i, int id) {
// build verteces for the hyjets stored events
double x0 = (hj2->GetX()).at(i);
double y0 = (hj2->GetY()).at(i);
double x = x0 * cosphi0_ - y0 * sinphi0_;
double y = y0 * cosphi0_ + x0 * sinphi0_;
double z = (hj2->GetZ()).at(i);
double t = (hj2->GetT()).at(i);

// convert to mm (as in PYTHIA6)
const double fm_to_mm = 1e-12;
double x = fm_to_mm * (x0 * cosphi0_ - y0 * sinphi0_);
double y = fm_to_mm * (y0 * cosphi0_ + x0 * sinphi0_);
double z = fm_to_mm * (hj2->GetZ()).at(i);
double t = fm_to_mm * (hj2->GetT()).at(i);

HepMC::GenVertex *vertex = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t), id);

Expand Down

0 comments on commit 1ad0fb5

Please sign in to comment.