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

@@ -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

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]
 
 
 

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);

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);

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<HepMC::GenEvent> eventVector;
   eventVector.reserve(numPileUp + 1);

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<const HepMC::GenVertex*, fcc::GenVertex> 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);

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
 

Generation/src/components/PythiaInterface.cpp

@@ -1,5 +1,4 @@
 #include "PythiaInterface.h"
-#include "Generation/Units.h"
 
 #include "GaudiKernel/IIncidentSvc.h"
 #include "GaudiKernel/Incident.h"

init.sh

@@ -1 +1 @@
-init_lcg_97a_FCC_2.sh
+init_lcg_97a_FCC_4.sh

init_lcg_97a_FCC_2.sh → 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