From 362e70b0ba4b41647fc3f0672d115b2bd494845d Mon Sep 17 00:00:00 2001 From: vvolkl Date: Wed, 4 Nov 2020 22:46:37 +0100 Subject: [PATCH] Fix HEPMC units (#419) * up * explicitly set hepmc genevent to use GEV/MM * up * up * up * up * up * up * up * adapt init * up * up * up --- .github/workflows/test.yml | 2 +- Generation/Generation/Units.h | 32 ------------------- Generation/options/particleGun.py | 7 ++-- .../src/components/ConstPtParticleGun.cpp | 24 ++++++++++++-- .../src/components/EDMToHepMCConverter.cpp | 18 ++++------- Generation/src/components/GenAlg.cpp | 1 + .../src/components/HepMCToEDMConverter.cpp | 32 +++++++------------ .../components/MomentumRangeParticleGun.cpp | 7 +++- Generation/src/components/PythiaInterface.cpp | 1 - init.sh | 2 +- ..._lcg_97a_FCC_2.sh => init_lcg_97a_FCC_4.sh | 4 ++- 11 files changed, 55 insertions(+), 75 deletions(-) delete mode 100644 Generation/Generation/Units.h rename init_lcg_97a_FCC_2.sh => init_lcg_97a_FCC_4.sh (52%) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index c50d4584c..70e02cb9d 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -8,7 +8,7 @@ jobs: strategy: fail-fast: false matrix: - SETUP: ["init_key4hep.sh", 'init_lcg_97a_FCC_2.sh', 'init_lcg_98.sh'] + SETUP: ["init_key4hep.sh", 'init_lcg_97a_FCC_4.sh', 'init_lcg_98.sh'] steps: - uses: actions/checkout@v2 - name: Install CVMFS diff --git a/Generation/Generation/Units.h b/Generation/Generation/Units.h deleted file mode 100644 index b3b94a4ca..000000000 --- a/Generation/Generation/Units.h +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef GENERATION_UNITS_H -#define GENERATION_UNITS_H - -#include "GaudiKernel/SystemOfUnits.h" -#include "HepMC/Units.h" - -// FIXME: All these should be a constexpr, but in CLHEP / HepMC they are only const - -namespace gen { -/// Default EDM units for the Generation package -namespace edmdefault { -const HepMC::Units::LengthUnit length = HepMC::Units::MM; -const HepMC::Units::MomentumUnit energy = HepMC::Units::GEV; -} -/// Default HepMC units -namespace hepmcdefault { -const HepMC::Units::LengthUnit length = HepMC::Units::MM; -const HepMC::Units::MomentumUnit energy = HepMC::Units::GEV; -} -/// Conversion factors from EDM to HepMC -namespace edm2hepmc { -const double length = HepMC::Units::conversion_factor(edmdefault::length, hepmcdefault::length); -const double energy = HepMC::Units::conversion_factor(edmdefault::energy, hepmcdefault::energy); -} -/// Conversion factors from HepMC to EDM -namespace hepmc2edm { -const double length = HepMC::Units::conversion_factor(hepmcdefault::length, edmdefault::length); -const double energy = HepMC::Units::conversion_factor(hepmcdefault::energy, edmdefault::energy); -} -} - -#endif /*GENERATION_UNITS_H*/ diff --git a/Generation/options/particleGun.py b/Generation/options/particleGun.py index 225537c49..6f8b69d93 100644 --- a/Generation/options/particleGun.py +++ b/Generation/options/particleGun.py @@ -5,7 +5,7 @@ ApplicationMgr( EvtSel='NONE', EvtMax=1, - OutputLevel=VERBOSE, + OutputLevel=INFO, ) from Configurables import FCCDataSvc @@ -15,7 +15,7 @@ from Configurables import ConstPtParticleGun -guntool1 = ConstPtParticleGun("SignalProvider", PdgCodes=[-211]) +guntool1 = ConstPtParticleGun("SignalProvider", PdgCodes=[-211], PtMin=50, PtMax=50) guntool2 = ConstPtParticleGun("PileUpProvider", PdgCodes=[11], writeParticleGunBranches=False) from Configurables import FlatSmearVertex smeartool = FlatSmearVertex() @@ -63,11 +63,12 @@ THistSvc().PrintAll=True THistSvc().AutoSave=True THistSvc().AutoFlush=True -THistSvc().OutputLevel=VERBOSE +THistSvc().OutputLevel=INFO from Configurables import PodioOutput out = PodioOutput("out", filename = "output_particleGun.root") out.outputCommands = ["keep *"] +ApplicationMgr().TopAlg += [out] diff --git a/Generation/src/components/ConstPtParticleGun.cpp b/Generation/src/components/ConstPtParticleGun.cpp index 60748c11c..c01e95666 100644 --- a/Generation/src/components/ConstPtParticleGun.cpp +++ b/Generation/src/components/ConstPtParticleGun.cpp @@ -86,10 +86,11 @@ void ConstPtParticleGun::generateParticle(Gaudi::LorentzVector& momentum, Gaudi: debug() << " -> " << pdgId << endmsg << " P = " << momentum << endmsg; /// additional branches in rootfile if (m_writeParticleGunBranches) { + const double hepmcMomentumConversionFactor = 0.001; float* _particlegun_phi = new float(phi); float* _particlegun_eta = new float(eta); float* _particlegun_costheta = new float (cos(2*atan(exp(eta)))); - float* _particlegun_pt = new float(pt); + float* _particlegun_pt = new float(pt * hepmcMomentumConversionFactor); m_datahandle_particlegun_pt->put(_particlegun_pt); m_datahandle_particlegun_eta->put(_particlegun_eta); m_datahandle_particlegun_costheta->put(_particlegun_costheta); @@ -103,13 +104,30 @@ StatusCode ConstPtParticleGun::getNextEvent(HepMC::GenEvent& theEvent) { // note: pgdid is set in function generateParticle int thePdgId; generateParticle(theFourMomentum, origin, thePdgId); + + + // create HepMC Vertex -- // by calling add_vertex(), the hepmc event is given ownership of the vertex - HepMC::GenVertex* v = new HepMC::GenVertex(HepMC::FourVector(origin.X(), origin.Y(), origin.Z(), origin.T())); + + HepMC::GenVertex* v = new HepMC::GenVertex(HepMC::FourVector( + origin.X(), + origin.Y(), + origin.Z(), + origin.T() + ) + ); // create HepMC particle -- // by calling add_particle_out(), the hepmc vertex is given ownership of the particle + // need to convert from Gaudi Units (MeV) to (GeV) + const double hepmcMomentumConversionFactor = 0.001; HepMC::GenParticle* p = new HepMC::GenParticle( - HepMC::FourVector(theFourMomentum.Px(), theFourMomentum.Py(), theFourMomentum.Pz(), theFourMomentum.E()), + HepMC::FourVector( + theFourMomentum.Px() * hepmcMomentumConversionFactor, + theFourMomentum.Py() * hepmcMomentumConversionFactor, + theFourMomentum.Pz() * hepmcMomentumConversionFactor, + theFourMomentum.E() * hepmcMomentumConversionFactor + ), thePdgId, 1); // hepmc status code for final state particle v->add_particle_out(p); diff --git a/Generation/src/components/EDMToHepMCConverter.cpp b/Generation/src/components/EDMToHepMCConverter.cpp index c038e45c8..05898c06e 100644 --- a/Generation/src/components/EDMToHepMCConverter.cpp +++ b/Generation/src/components/EDMToHepMCConverter.cpp @@ -5,10 +5,8 @@ #include "datamodel/GenVertexCollection.h" #include "datamodel/MCParticleCollection.h" -#include "Generation/Units.h" using HepMC::GenParticle; -using HepMC::Units::conversion_factor; DECLARE_COMPONENT(EDMToHepMCConverter) @@ -25,27 +23,23 @@ StatusCode EDMToHepMCConverter::execute() { const fcc::MCParticleCollection* particles = m_genphandle.get(); // ownership of event given to data service at the end of execute HepMC::GenEvent* event = new HepMC::GenEvent; + event->use_units(HepMC::Units::GEV, HepMC::Units::MM); - // conversion of units to EDM standard units: - // First cover the case that hepMC file is not in expected units and then convert to EDM default - double hepmc2EdmLength = conversion_factor(event->length_unit(), gen::hepmcdefault::length) * gen::hepmc2edm::length; - double hepmc2EdmEnergy = - conversion_factor(event->momentum_unit(), gen::hepmcdefault::energy) * gen::hepmc2edm::energy; for (auto p : *(particles)) { if (p.status() == 1) { // only final state particles GenParticle* pHepMC = - new GenParticle(HepMC::FourVector(p.p4().px, p.p4().py, p.p4().pz, p.p4().mass / hepmc2EdmEnergy), + new GenParticle(HepMC::FourVector(p.p4().px, p.p4().py, p.p4().pz, p.p4().mass), p.pdgId(), p.status()); // hepmc status code for final state particle fcc::ConstGenVertex vStart = p.startVertex(); if (p.startVertex().isAvailable()) { HepMC::GenVertex* v = - new HepMC::GenVertex(HepMC::FourVector(vStart.position().x / hepmc2EdmLength, - vStart.position().y / hepmc2EdmLength, - vStart.position().z / hepmc2EdmLength, - vStart.ctau() / Gaudi::Units::c_light / hepmc2EdmLength)); + new HepMC::GenVertex(HepMC::FourVector(vStart.position().x, + vStart.position().y, + vStart.position().z, + vStart.ctau() / Gaudi::Units::c_light)); v->add_particle_out(pHepMC); event->add_vertex(v); diff --git a/Generation/src/components/GenAlg.cpp b/Generation/src/components/GenAlg.cpp index 6d08846ad..0b63cb7b0 100644 --- a/Generation/src/components/GenAlg.cpp +++ b/Generation/src/components/GenAlg.cpp @@ -29,6 +29,7 @@ StatusCode GenAlg::initialize() { StatusCode GenAlg::execute() { auto theEvent = m_hepmchandle.createAndPut(); + theEvent->use_units(HepMC::Units::GEV, HepMC::Units::MM); const unsigned int numPileUp = m_pileUpTool->numberOfPileUp(); std::vector eventVector; eventVector.reserve(numPileUp + 1); diff --git a/Generation/src/components/HepMCToEDMConverter.cpp b/Generation/src/components/HepMCToEDMConverter.cpp index 5f7dae1da..a28926cb1 100644 --- a/Generation/src/components/HepMCToEDMConverter.cpp +++ b/Generation/src/components/HepMCToEDMConverter.cpp @@ -5,7 +5,6 @@ #include "datamodel/GenVertexCollection.h" #include "datamodel/MCParticleCollection.h" -#include "Generation/Units.h" #include "HepPDT/ParticleID.hh" DECLARE_COMPONENT(HepMCToEDMConverter) @@ -24,13 +23,6 @@ StatusCode HepMCToEDMConverter::execute() { const HepMC::GenEvent* event = m_hepmchandle.get(); fcc::MCParticleCollection* particles = new fcc::MCParticleCollection(); fcc::GenVertexCollection* vertices = new fcc::GenVertexCollection(); - - // conversion of units to EDM standard units: - // First cover the case that hepMC file is not in expected units and then convert to EDM default - double hepmc2EdmLength = - HepMC::Units::conversion_factor(event->length_unit(), gen::hepmcdefault::length) * gen::hepmc2edm::length; - double hepmc2EdmEnergy = - HepMC::Units::conversion_factor(event->momentum_unit(), gen::hepmcdefault::energy) * gen::hepmc2edm::energy; // bookkeeping of particle / vertex relations std::unordered_map hepmcToEdmVertexMap; @@ -51,10 +43,10 @@ StatusCode HepMCToEDMConverter::execute() { auto& p4 = particle.p4(); tmp = (*particle_i)->momentum(); - p4.px = tmp.px() * hepmc2EdmEnergy; - p4.py = tmp.py() * hepmc2EdmEnergy; - p4.pz = tmp.pz() * hepmc2EdmEnergy; - p4.mass = (*particle_i)->generated_mass() * hepmc2EdmEnergy; + p4.px = tmp.px(); + p4.py = tmp.py(); + p4.pz = tmp.pz(); + p4.mass = (*particle_i)->generated_mass(); /// create production vertex, if it has not already been created and logged in the map HepMC::GenVertex* productionVertex = (*particle_i)->production_vertex(); @@ -66,10 +58,10 @@ StatusCode HepMCToEDMConverter::execute() { tmp = productionVertex->position(); auto vertex = vertices->create(); auto& position = vertex.position(); - position.x = tmp.x() * hepmc2EdmLength; - position.y = tmp.y() * hepmc2EdmLength; - position.z = tmp.z() * hepmc2EdmLength; - vertex.ctau(tmp.t() * Gaudi::Units::c_light * hepmc2EdmLength); // is ctau like this? + position.x = tmp.x(); + position.y = tmp.y(); + position.z = tmp.z(); + vertex.ctau(tmp.t() * Gaudi::Units::c_light); // is ctau like this? // add vertex to map for further particles hepmcToEdmVertexMap.insert({productionVertex, vertex}); particle.startVertex(vertex); @@ -86,10 +78,10 @@ StatusCode HepMCToEDMConverter::execute() { tmp = decayVertex->position(); auto vertex = vertices->create(); auto& position = vertex.position(); - position.x = tmp.x() * hepmc2EdmLength; - position.y = tmp.y() * hepmc2EdmLength; - position.z = tmp.z() * hepmc2EdmLength; - vertex.ctau(tmp.t() * Gaudi::Units::c_light * hepmc2EdmLength); // is ctau like this? + position.x = tmp.x(); + position.y = tmp.y(); + position.z = tmp.z(); + vertex.ctau(tmp.t() * Gaudi::Units::c_light); // is ctau like this? // add vertex to map for further particles hepmcToEdmVertexMap.insert({decayVertex, vertex}); particle.endVertex(vertex); diff --git a/Generation/src/components/MomentumRangeParticleGun.cpp b/Generation/src/components/MomentumRangeParticleGun.cpp index 984315a9f..eb5df9b10 100644 --- a/Generation/src/components/MomentumRangeParticleGun.cpp +++ b/Generation/src/components/MomentumRangeParticleGun.cpp @@ -113,8 +113,13 @@ StatusCode MomentumRangeParticleGun::getNextEvent(HepMC::GenEvent& theEvent) { HepMC::GenVertex* v = new HepMC::GenVertex(HepMC::FourVector(origin.X(), origin.Y(), origin.Z(), origin.T())); // create HepMC particle -- // by calling add_particle_out(), the hepmc vertex is given ownership of the particle + const double hepmcMomentumConversionFactor = 0.001; HepMC::GenParticle* p = new HepMC::GenParticle( - HepMC::FourVector(theFourMomentum.Px(), theFourMomentum.Py(), theFourMomentum.Pz(), theFourMomentum.E()), + HepMC::FourVector(theFourMomentum.Px() * hepmcMomentumConversionFactor, + theFourMomentum.Py() * hepmcMomentumConversionFactor, + theFourMomentum.Pz() * hepmcMomentumConversionFactor, + theFourMomentum.E() * hepmcMomentumConversionFactor + ), thePdgId, 1); // hepmc status code for final state particle diff --git a/Generation/src/components/PythiaInterface.cpp b/Generation/src/components/PythiaInterface.cpp index 3f0eeff58..954316a8e 100644 --- a/Generation/src/components/PythiaInterface.cpp +++ b/Generation/src/components/PythiaInterface.cpp @@ -1,5 +1,4 @@ #include "PythiaInterface.h" -#include "Generation/Units.h" #include "GaudiKernel/IIncidentSvc.h" #include "GaudiKernel/Incident.h" diff --git a/init.sh b/init.sh index 215dd7e88..ab54b0bf8 120000 --- a/init.sh +++ b/init.sh @@ -1 +1 @@ -init_lcg_97a_FCC_2.sh \ No newline at end of file +init_lcg_97a_FCC_4.sh \ No newline at end of file diff --git a/init_lcg_97a_FCC_2.sh b/init_lcg_97a_FCC_4.sh similarity index 52% rename from init_lcg_97a_FCC_2.sh rename to init_lcg_97a_FCC_4.sh index 1cae60d3d..7b8d886e7 100644 --- a/init_lcg_97a_FCC_2.sh +++ b/init_lcg_97a_FCC_4.sh @@ -1,9 +1,11 @@ LCGPREFIX=/cvmfs/sft.cern.ch/lcg export BINARY_TAG=x86_64-centos7-gcc8-opt -LCGPATH=$LCGPREFIX/views/LCG_97a_FCC_2/$BINARY_TAG +LCGPATH=$LCGPREFIX/views/LCG_97a_FCC_4/$BINARY_TAG source $LCGPATH/bin/thisdd4hep_only.sh source $LCGPATH/setup.sh +export PATH=/cvmfs/sft.cern.ch/lcg/releases/CMake/3.17.3-75516/x86_64-centos7-gcc8-opt/bin:$PATH +export PATH=/cvmfs/sft.cern.ch/lcg/releases/ninja/1.10.0-3d8f6/x86_64-centos7-gcc8-opt/bin:$PATH export Gaudi_DIR=$(dirname $(readlink -f "$(which gaudirun.py)"))/../ export CMAKE_PREFIX_PATH=$Gaudi_DIR:$CMAKE_PREFIX_PATH