Skip to content

Apply correction for energy (non)conservation in coalescence for c-deuteron #2372

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

Merged
stefanopolitano merged 1 commit into from
Jun 3, 2026
Merged

Apply correction for energy (non)conservation in coalescence for c-deuteron #2372

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 30 additions & 0 deletions MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,9 @@
#include <string>
#include <vector>

#include "Math/Vector3D.h"
#include "Math/Vector4D.h"

R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
#include "MC/config/common/external/generator/CoalescencePythia8.h"

Expand Down Expand Up @@ -240,6 +243,7 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
// we try the coalescence here, if successful we copy particles in the pythia event and we move to the next charm nucleus
isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector<unsigned int>{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1], 10.);
if (isCoalSuccess) {
restoreEnergyConservation(mPythiaGun.event, idxCharmNucleus);
int offset = mPythia.event.size(); // we need to rescale the indices of mothers and daughters, accounting for the particles that are already appended to the event
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
auto part = mPythiaGun.event[iPart];
Expand Down Expand Up @@ -274,6 +278,32 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8


private:

void restoreEnergyConservation(Pythia8::Event& event, int idxCharmNucleus, float targetMassTolerance = 1e-5) {
/// In the coalescence, the energy is not conserved, we rescale the momentum of the charmed nuclei daughters to restore it to avoid changes in the invariant mass of the charmed nucleus

float scale = 1.f;
float invMass{0.f};
while (abs(invMass - mMassCharmNucleus) > targetMassTolerance) {
ROOT::Math::PxPyPzMVector fourVecCharmNucleus;
for (int iDau{event[idxCharmNucleus].daughter1()}; iDau<=event[idxCharmNucleus].daughter2(); ++iDau) {
auto dau = event[iDau];
fourVecCharmNucleus += ROOT::Math::PxPyPzMVector(dau.px() * scale, dau.py() * scale, dau.pz() * scale, dau.m());
}
invMass = fourVecCharmNucleus.M();
scale *= mMassCharmNucleus / invMass;
}

for (int iDau{event[idxCharmNucleus].daughter1()}; iDau<=event[idxCharmNucleus].daughter2(); ++iDau) {
event[iDau].px(event[iDau].px() * scale);
event[iDau].py(event[iDau].py() * scale);
event[iDau].pz(event[iDau].pz() * scale);
}
event[idxCharmNucleus].px(event[idxCharmNucleus].px() * scale);
event[idxCharmNucleus].py(event[idxCharmNucleus].py() * scale);
event[idxCharmNucleus].pz(event[idxCharmNucleus].pz() * scale);
}

// Properties of selection
float mMassCharmNucleus; /// mass of the charmed nucleus
int mPdgCharmNucleus; /// pdg code of the charmed nucleus
Expand Down
Loading