Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Hydjet2 bugfixes #38443

Merged
merged 3 commits into from Aug 11, 2022
Merged

Hydjet2 bugfixes #38443

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b310a57
Hydjet2 bugfixes for status codes and units
Jun 20, 2022
a3b5590
Hydjet2 bugfixes for status codes and units
Jun 20, 2022
fbba702
perrota's comments on units conversation
Jun 21, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
4 changes: 2 additions & 2 deletions GeneratorInterface/Hydjet2Interface/interface/Hydjet2Hadronizer.h
View file
Open in desktop
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
Open in desktop
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
87 changes: 49 additions & 38 deletions GeneratorInterface/Hydjet2Interface/src/Hydjet2Hadronizer.cc
View file
Open in desktop
@@ -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,12 @@ 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)
double x = 0.000000000001 * (x0 * cosphi0_ - y0 * sinphi0_);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is less error prone, and more intuitive, if you define the conversion in a contexpr parameter

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How it is possible? We do not know the x0 and y0 at compile time...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

e.g.

  double constexpr am_to_mm = 0.000000000001;  // or even 1e-12, which is more intuitive
  double x =  am_to_mm * (x0 * cosphi0_ - y0 * sinphi0_);

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I see. But what is the reason of constexpr here? Why not just:

const double fm_to_mm = 1e-12; 
double x =  fm_to_mm * (x0 * cosphi0_ - y0 * sinphi0_);

?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I see. But what is the reason of constexpr here? Why not just:

const double fm_to_mm = 1e-12; 
double x =  fm_to_mm * (x0 * cosphi0_ - y0 * sinphi0_);

?

Well, constexpr is initialized at compile time, while const it is at runtime. Not a big difference here, let it as you updated it.

double y = 0.000000000001 * (y0 * cosphi0_ + x0 * sinphi0_);
double z = 0.000000000001 * (hj2->GetZ()).at(i);
double t = 0.000000000001 * (hj2->GetT()).at(i);

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

Expand Down