From f4846ccab854af0e44b8806d9280a20662161b8d Mon Sep 17 00:00:00 2001 From: swenzel Date: Tue, 25 Mar 2025 13:02:28 +0100 Subject: [PATCH] GeneratorHybrid: improve unit treatment So far, units are treated solely in the Generator::addTrack function. This works well for fundamental generators. However, the hybrid generator is a meta generator potentially consisting of a collection of underlying generators that may have completely different units. This may currently lead to wrong generator output, in certain cases. This commit fixes these bugs and introduces unit handling within GeneratorHybrid. --- Generators/include/Generators/Generator.h | 4 ++ .../include/Generators/GeneratorHybrid.h | 1 - Generators/src/GeneratorHybrid.cxx | 60 +++++++++++++++---- 3 files changed, 52 insertions(+), 13 deletions(-) diff --git a/Generators/include/Generators/Generator.h b/Generators/include/Generators/Generator.h index 640cc80931862..bd35a00793e2d 100644 --- a/Generators/include/Generators/Generator.h +++ b/Generators/include/Generators/Generator.h @@ -78,9 +78,13 @@ class Generator : public FairGenerator /** setters **/ void setMomentumUnit(double val) { mMomentumUnit = val; }; + double getMomentumUnit() const { return mMomentumUnit; } void setEnergyUnit(double val) { mEnergyUnit = val; }; + double getEnergyUnit() const { return mEnergyUnit; } void setPositionUnit(double val) { mPositionUnit = val; }; + double getPositionUnit() const { return mPositionUnit; } void setTimeUnit(double val) { mTimeUnit = val; }; + double getTimeUnit() const { return mTimeUnit; } void setBoost(Double_t val) { mBoost = val; }; void setTriggerMode(ETriggerMode_t val) { mTriggerMode = val; }; void addTrigger(Trigger trigger) { mTriggers.push_back(trigger); }; diff --git a/Generators/include/Generators/GeneratorHybrid.h b/Generators/include/Generators/GeneratorHybrid.h index 955240d6a28fa..b92437b02d874 100644 --- a/Generators/include/Generators/GeneratorHybrid.h +++ b/Generators/include/Generators/GeneratorHybrid.h @@ -54,7 +54,6 @@ class GeneratorHybrid : public Generator { public: - GeneratorHybrid() = default; GeneratorHybrid(const std::string& inputgens); ~GeneratorHybrid(); diff --git a/Generators/src/GeneratorHybrid.cxx b/Generators/src/GeneratorHybrid.cxx index 729d69527c384..83a694703c259 100644 --- a/Generators/src/GeneratorHybrid.cxx +++ b/Generators/src/GeneratorHybrid.cxx @@ -25,6 +25,12 @@ namespace eventgen GeneratorHybrid::GeneratorHybrid(const std::string& inputgens) { + // This generator has trivial unit conversions + setTimeUnit(1.); + setPositionUnit(1.); + setMomentumUnit(1.); + setEnergyUnit(1.); + if (!parseJSON(inputgens)) { LOG(fatal) << "Failed to parse JSON configuration from input generators"; exit(1); @@ -382,6 +388,27 @@ bool GeneratorHybrid::importParticles() } } } + + auto unit_transformer = [](auto& p, auto pos_unit, auto time_unit, auto en_unit, auto mom_unit) { + p.SetMomentum(p.Px() * mom_unit, p.Py() * mom_unit, p.Pz() * mom_unit, p.Energy() * en_unit); + p.SetProductionVertex(p.Vx() * pos_unit, p.Vy() * pos_unit, p.Vz() * pos_unit, p.T() * time_unit); + }; + + auto index_transformer = [](auto& p, int offset) { + for (int i = 0; i < 2; ++i) { + if (p.GetMother(i) != -1) { + const auto newindex = p.GetMother(i) + offset; + p.SetMother(i, newindex); + } + } + if (p.GetNDaughters() > 0) { + for (int i = 0; i < 2; ++i) { + const auto newindex = p.GetDaughter(i) + offset; + p.SetDaughter(i, newindex); + } + } + }; + // Clear particles and event header mParticles.clear(); mMCEventHeader.clearInfo(); @@ -391,23 +418,20 @@ bool GeneratorHybrid::importParticles() LOG(info) << "Importing particles for task " << subIndex; auto subParticles = gens[subIndex]->getParticles(); + auto time_unit = gens[subIndex]->getTimeUnit(); + auto pos_unit = gens[subIndex]->getPositionUnit(); + auto mom_unit = gens[subIndex]->getMomentumUnit(); + auto energy_unit = gens[subIndex]->getEnergyUnit(); + // The particles carry mother and daughter indices, which are relative // to the sub-generator. We need to adjust these indices to reflect that particles // are now embedded into a cocktail. auto offset = mParticles.size(); for (auto& p : subParticles) { - for (int i = 0; i < 2; ++i) { - if (p.GetMother(i) != -1) { - const auto newindex = p.GetMother(i) + offset; - p.SetMother(i, newindex); - } - } - if (p.GetNDaughters() > 0) { - for (int i = 0; i < 2; ++i) { - const auto newindex = p.GetDaughter(i) + offset; - p.SetDaughter(i, newindex); - } - } + // apply the mother-daugher index transformation + index_transformer(p, offset); + // apply unit transformation of sub-generator + unit_transformer(p, pos_unit, time_unit, energy_unit, mom_unit); } mParticles.insert(mParticles.end(), subParticles.begin(), subParticles.end()); @@ -420,6 +444,18 @@ bool GeneratorHybrid::importParticles() LOG(info) << "Importing particles for task " << genIndex; // at this moment the mIndex-th generator is ready to be used mParticles = gens[genIndex]->getParticles(); + + auto time_unit = gens[genIndex]->getTimeUnit(); + auto pos_unit = gens[genIndex]->getPositionUnit(); + auto mom_unit = gens[genIndex]->getMomentumUnit(); + auto energy_unit = gens[genIndex]->getEnergyUnit(); + + // transform units to units of the hybrid generator + for (auto& p : mParticles) { + // apply unit transformation + unit_transformer(p, pos_unit, time_unit, energy_unit, mom_unit); + } + // fetch the event Header information from the underlying generator gens[genIndex]->updateHeader(&mMCEventHeader); mInputTaskQueue.push(genIndex);