[delphes] save trkcov information (#399)

* [delphes] save trkcov information * use LCG_97a_FCC_2 * use LCG_97a_FCC_2; add missing setup vars * [ci] gaudi xenv workarounds * switch to spack stack * [ci] update for spack stack * Update Sim/SimDelphesInterface/options/PythiaDelphes_config_IDEAtrkCov.py * Update Sim/SimDelphesInterface/src/DelphesSaveChargedParticles.cpp Co-authored-by: clementhelsens <clement.helsens@cern.ch> * Update Sim/SimDelphesInterface/src/DelphesSaveChargedParticles.cpp * fix to compile with podio v0.11 * fix tests Co-authored-by: clementhelsens <clement.helsens@cern.ch>
HEP-FCC · Sep 19, 2020 · d337c0f · d337c0f
1 parent 7e0ce9d
commit d337c0f
Show file tree

Hide file tree

Showing 9 changed files with 233 additions and 71 deletions.
diff --git a/.github/scripts/compile_and_test.sh b/.github/scripts/compile_and_test.sh
@@ -4,7 +4,7 @@ cd /Package
 source init.sh
 mkdir build install
 cd build
-cmake -DCMAKE_INSTALL_PREFIX=../install -DCMAKE_CXX_STANDARD=17  .. && \
+cmake -DCMAKE_INSTALL_PREFIX=../install -DCMAKE_CXX_STANDARD=17   .. && \
 make  -j `getconf _NPROCESSORS_ONLN` && \
 make install && \
 ls -alhR $G4LEDATA && \

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -33,5 +33,5 @@ jobs:
         docker run -it --name CI_container -v ${GITHUB_WORKSPACE}:/Package -v /cvmfs:/cvmfs:shared -d clicdp/cc7-lcg /bin/bash
     - name: Compile and test
       run: |
-        docker exec CI_container /bin/bash -c "ln -s /usr/lib64/liblzma.so.5.2.2 /usr/lib64/liblzma.so;cd ./Package; ./.github/scripts/compile_and_test.sh;"
+        docker exec CI_container /bin/bash -c "yum -y install environment-modules; source /etc/profile.d/modules.sh; ln -s /usr/lib64/liblzma.so.5.2.2 /usr/lib64/liblzma.so; cd ./Package; ./.github/scripts/compile_and_test.sh;"
 
diff --git a/.gitignore b/.gitignore
@@ -10,6 +10,7 @@
 
 *.log
 build*
+install*
 build.*.log
 Makefile
 NMake

diff --git a/Examples/options/geant_pgun_fullsim.py b/Examples/options/geant_pgun_fullsim.py
@@ -18,6 +18,7 @@
 # MomentumRangeParticleGun generates particles of given type(s) within given momentum, phi and theta range
 # FlatSmearVertex smears the vertex with uniform distribution
 guntool = MomentumRangeParticleGun()
+guntool.PdgCodes = [11]
 gen = GenAlg("ParticleGun", SignalProvider=guntool, VertexSmearingTool="FlatSmearVertex")
 gen.hepmc.Path = "hepmc"
 
@@ -32,11 +33,12 @@
 ## DD4hep geometry service
 # Parses the given xml file
 geoservice = GeoSvc("GeoSvc", detectors=['file:Detector/DetFCChhBaseline1/compact/FCChh_DectEmptyMaster.xml',
-                                         'file:Detector/DetFCChhTrackerTkLayout/compact/Tracker.xml',
+                                         #'file:Detector/DetFCChhTrackerTkLayout/compact/Tracker.xml',
                                          'file:Detector/DetFCChhECalInclined/compact/FCChh_ECalBarrel_withCryostat.xml',
                                          'file:Detector/DetFCChhCalDiscs/compact/Endcaps_coneCryo.xml',
-                                         'file:Detector/DetFCChhCalDiscs/compact/Forward_coneCryo.xml',
-                                         'file:Detector/DetFCChhHCalTile/compact/FCChh_HCalBarrel_TileCal.xml'],
+                                         #'file:Detector/DetFCChhCalDiscs/compact/Forward_coneCryo.xml',
+                                         'file:Detector/DetFCChhHCalTile/compact/FCChh_HCalBarrel_TileCal.xml'
+                                         ],
                     OutputLevel = INFO)
 
 from Configurables import SimG4Svc
@@ -51,9 +53,9 @@
 # first, create a tool that saves the tracker hits
 # Name of that tool in GAUDI is "XX/YY" where XX is the tool class name ("SimG4SaveTrackerHits")
 # and YY is the given name ("saveTrackerHits")
-savetrackertool = SimG4SaveTrackerHits("saveTrackerHits", readoutNames = ["TrackerBarrelReadout", "TrackerEndcapReadout"])
-savetrackertool.positionedTrackHits.Path = "positionedHits"
-savetrackertool.trackHits.Path = "hits"
+#savetrackertool = SimG4SaveTrackerHits("saveTrackerHits", readoutNames = ["TrackerBarrelReadout", "TrackerEndcapReadout"])
+#savetrackertool.positionedTrackHits.Path = "positionedHits"
+#savetrackertool.trackHits.Path = "hits"
 # and a tool that saves the calorimeter hits with a name "SimG4SaveCalHits/saveCalHits"
 saveecaltool = SimG4SaveCalHits("saveECalBarrelHits", readoutNames = ["ECalBarrelEta"])
 saveecaltool.positionedCaloHits.Path = "ECalBarrelPositionedHits"
@@ -71,7 +73,7 @@
 particle_converter = SimG4PrimariesFromEdmTool("EdmConverter")
 particle_converter.genParticles.Path = "allGenParticles"
 geantsim = SimG4Alg("SimG4Alg",
-                    outputs = ["SimG4SaveTrackerHits/saveTrackerHits", "SimG4SaveCalHits/saveECalBarrelHits", "SimG4SaveCalHits/saveECalEndcapHits", "SimG4SaveCalHits/saveECalFwdHits", "SimG4SaveCalHits/saveHCalHits"],
+                    outputs = [ "SimG4SaveCalHits/saveECalBarrelHits",  "SimG4SaveCalHits/saveECalEndcapHits",  "SimG4SaveCalHits/saveHCalHits"],
                     eventProvider=particle_converter)
 
 from Configurables import PodioOutput

diff --git a/Sim/SimDelphesInterface/CMakeLists.txt b/Sim/SimDelphesInterface/CMakeLists.txt
@@ -12,6 +12,7 @@ find_package(PODIO)
 find_package(ROOT COMPONENTS Physics EG)
 find_package(HepMC)
 
+
 gaudi_add_module(SimDelphesInterface
                  src/*.cpp
                  INCLUDE_DIRS FWCore ROOT Delphes FCCEDM PODIO HepMC
@@ -26,3 +27,9 @@ include(CTest)
 gaudi_add_test(PythiaDelphes
                WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
                FRAMEWORK options/PythiaDelphes_config.py)
+
+
+gaudi_add_test(PythiaDelphes_IDEAtrkCov
+               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+                COMMAND gaudirun.py  Sim/SimDelphesInterface/options/PythiaDelphes_config_IDEAtrkCov.py
+                )
diff --git a/Sim/SimDelphesInterface/options/PythiaDelphes_config_IDEAtrkCov.py b/Sim/SimDelphesInterface/options/PythiaDelphes_config_IDEAtrkCov.py
@@ -0,0 +1,149 @@
+import sys, os
+from Gaudi.Configuration import *
+
+# Workflow Steering
+from Configurables import ApplicationMgr
+ApplicationMgr().EvtSel = 'NONE'
+ApplicationMgr().EvtMax = 100
+
+## Data event model based on Podio
+from Configurables import FCCDataSvc
+podioEvent = FCCDataSvc("EventDataSvc")
+ApplicationMgr().ExtSvc += [podioEvent]
+ApplicationMgr().OutputLevel = INFO
+
+
+## Pythia generator
+from Configurables import PythiaInterface
+pythia8gentool = PythiaInterface()
+pythia8gentool.Filename = os.path.join(os.environ.get("FCCSWSHAREDIR", ""),"Generation/data/ee_Z_ddbar.cmd")
+
+
+## Write the HepMC::GenEvent to the data service
+from Configurables import GenAlg
+pythia8gen = GenAlg()
+pythia8gen.SignalProvider = pythia8gentool
+pythia8gen.hepmc.Path = "hepmc"
+ApplicationMgr().TopAlg += [pythia8gen]
+
+
+### Reads an HepMC::GenEvent from the data service and writes a collection of EDM Particles
+from Configurables import HepMCToEDMConverter
+hepmc_converter = HepMCToEDMConverter("Converter")
+hepmc_converter.hepmc.Path = "hepmc"
+hepmc_converter.genparticles.Path = "genParticles"
+hepmc_converter.genvertices.Path = "genVertices"
+ApplicationMgr().TopAlg += [hepmc_converter]
+
+
+# Define all output tools that convert the Delphes collections to FCC-EDM:
+
+from Configurables import DelphesSaveChargedParticles
+
+muonSaveTool = DelphesSaveChargedParticles("muons")
+muonSaveTool.delphesArrayName = "MuonFilter/muons"
+muonSaveTool.particles.Path      = "muons"
+muonSaveTool.particles_trkCov.Path      = "muons_trkCov"
+muonSaveTool.mcAssociations.Path = "muonsToMC"
+muonSaveTool.isolationTags.Path  = "muonITags"
+
+eleSaveTool = DelphesSaveChargedParticles("electrons")
+eleSaveTool.delphesArrayName = "ElectronFilter/electrons"
+eleSaveTool.particles.Path      = "electrons"
+eleSaveTool.particles_trkCov.Path      = "electrons_trkCov"
+eleSaveTool.mcAssociations.Path = "electronsToMC"
+eleSaveTool.isolationTags.Path  = "electronITags"
+
+chhadSaveTool = DelphesSaveChargedParticles("efcharged")
+chhadSaveTool.delphesArrayName = "Calorimeter/eflowTracks"
+chhadSaveTool.saveIsolation = False
+chhadSaveTool.particles.Path      = "efcharged"
+chhadSaveTool.particles_trkCov.Path      = "efcharged_trkCov"
+chhadSaveTool.mcAssociations.Path = "efchargedToMC"
+
+
+from Configurables import DelphesSaveNeutralParticles
+
+# Particle-Flow Photons output tool
+pfphotonsSaveTool = DelphesSaveNeutralParticles("efphotons")
+pfphotonsSaveTool.delphesArrayName="Calorimeter/eflowPhotons"
+pfphotonsSaveTool.saveIsolation=False
+pfphotonsSaveTool.particles.Path      = "efphotons"
+pfphotonsSaveTool.mcAssociations.Path = "efphotonsToMC"
+pfphotonsSaveTool.isolationTags.Path  = "efphotonITags"
+
+# Photons output tool
+photonsSaveTool = DelphesSaveNeutralParticles("photons")
+photonsSaveTool.delphesArrayName = "PhotonEfficiency/photons"
+photonsSaveTool.particles.Path      = "photons"
+photonsSaveTool.mcAssociations.Path = "photonsToMC"
+photonsSaveTool.isolationTags.Path  = "photonITags"
+
+# Particle-Flow Neutral Hadrons output tool
+neuthadSaveTool = DelphesSaveNeutralParticles("efneutrals")
+neuthadSaveTool.delphesArrayName = "Calorimeter/eflowNeutralHadrons"
+neuthadSaveTool.saveIsolation = False
+neuthadSaveTool.particles.Path      = "efneutrals"
+neuthadSaveTool.mcAssociations.Path = "efneutralsToMC"
+
+
+from Configurables import DelphesSaveGenJets
+
+genJetSaveTool = DelphesSaveGenJets("genJets")
+genJetSaveTool.delphesArrayName = "GenJetFinder/jets"
+genJetSaveTool.genJets.Path             = "genJets"
+genJetSaveTool.genJetsFlavorTagged.Path = "genJetsFlavor"
+
+
+from Configurables import DelphesSaveJets
+
+jetSaveTool = DelphesSaveJets("jets")
+jetSaveTool.delphesArrayName = "JetEnergyScale/jets"
+jetSaveTool.jets.Path             = "jets"
+jetSaveTool.jetConstituents.Path  = "jetParts"
+jetSaveTool.jetsFlavorTagged.Path = "jetsFlavor"
+jetSaveTool.jetsBTagged.Path      = "bTags"
+jetSaveTool.jetsCTagged.Path      = "cTags"
+jetSaveTool.jetsTauTagged.Path    = "tauTags"
+
+
+from Configurables import DelphesSaveMet
+
+metSaveTool = DelphesSaveMet("met")
+metSaveTool.delphesMETArrayName = "MissingET/momentum"
+metSaveTool.delphesSHTArrayName = "ScalarHT/energy"
+metSaveTool.missingEt.Path = "met"
+
+
+## Delphes simulator -> define objects to be written out
+from Configurables import DelphesSimulation
+delphessim = DelphesSimulation()
+## Define Delphes card
+delphessim.DelphesCard = os.path.join(os.environ.get("DELPHES_DIR", ""), "cards/delphes_card_IDEAtrkCov.tcl")
+delphessim.ROOTOutputFile = ""
+delphessim.ApplyGenFilter = True
+delphessim.outputs = [
+                       "DelphesSaveChargedParticles/muons",
+                       "DelphesSaveChargedParticles/electrons",
+                       "DelphesSaveNeutralParticles/photons",
+                       "DelphesSaveChargedParticles/efcharged",
+                       "DelphesSaveNeutralParticles/efphotons",
+                       "DelphesSaveNeutralParticles/efneutrals",
+                       "DelphesSaveGenJets/genJets",
+                       "DelphesSaveJets/jets",
+                       "DelphesSaveMet/met",
+                     ]
+delphessim.hepmc.Path               = "hepmc"
+delphessim.genParticles.Path        = "skimmedGenParticles"
+delphessim.mcEventWeights.Path      = "mcEventWeights"
+ApplicationMgr().TopAlg += [delphessim]
+
+
+## FCC event-data model output -> define objects to be written out
+from Configurables import PodioOutput
+out = PodioOutput("out")
+out.filename = "FCCDelphesOutput.root"
+out.outputCommands = [
+                       "keep *", 
+                     ]
+ApplicationMgr().TopAlg += [out]
diff --git a/Sim/SimDelphesInterface/src/DelphesSaveChargedParticles.cpp b/Sim/SimDelphesInterface/src/DelphesSaveChargedParticles.cpp
@@ -6,12 +6,41 @@
 
 // datamodel
 #include "datamodel/MCParticleCollection.h"
+#include "datamodel/TrackStateCollection.h"
 #include "datamodel/ParticleCollection.h"
 #include "datamodel/ParticleMCParticleAssociationCollection.h"
 #include "datamodel/TaggedParticleCollection.h"
 // ROOT
 #include "SimDelphesInterface/ParticleStatus.h"
 #include "TObjArray.h"
+#include <TMatrixD.h>
+
+
+
+// Covariance conversion to ACTS format
+TMatrixDSym CovToACTS(TMatrixDSym Cov, double ct /*cot(theta)*/, double C /*half curvature*/, double fB=2. /*magnetic field in T*/)
+{
+  TMatrixDSym cACTS(6); cACTS.Zero();
+  Double_t b = -0.29988*fB / 2.;
+  //
+  // Fill derivative matrix
+  TMatrixD A(5, 5); A.Zero();
+  A(0, 0) = 1000.;    // D-D  conversion to mm
+  A(1, 2) = 1.0;    // phi0-phi0
+  A(2, 4) = 1.0/(TMath::Sqrt(1.0 + ct*ct) * b); // q/p-C
+  A(3, 1) = 1000.;    // z0-z0 conversion to mm
+  A(4, 3) = -1.0 / (1.0 + ct*ct); // theta - cot(theta)
+  A(4, 4) = -C*ct / (b*pow(1.0 + ct*ct,3.0/2.0)); // q/p-cot(theta)
+  //
+  TMatrixDSym Cv = Cov;
+  TMatrixD At(5, 5);
+  At.Transpose(A);
+  Cv.Similarity(At);
+  TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;
+  cACTS(5, 5) = 0.1;  // Currently undefined: set to arbitrary value to avoid crashes
+  //
+  return cACTS;
+}
 
 DECLARE_COMPONENT(DelphesSaveChargedParticles)
 
@@ -21,6 +50,7 @@ DelphesSaveChargedParticles::DelphesSaveChargedParticles(const std::string& aTyp
     : GaudiTool(aType, aName, aParent) {
   declareInterface<IDelphesSaveOutputTool>(this);
   declareProperty("particles", m_particles, "Handle the particles to be saved");
+  declareProperty("particles_trkCov", m_particles_trkCov, "Handle the particles to be saved");
   declareProperty("mcAssociations", m_mcAssociations, "Handle to associate particles with MCParticles");
   declareProperty("isolationTags", m_isolationTaggedParticles, "Handle for isolation tags");
 }
@@ -34,6 +64,7 @@ StatusCode DelphesSaveChargedParticles::finalize() { return GaudiTool::finalize(
 StatusCode DelphesSaveChargedParticles::saveOutput(Delphes& delphes, const fcc::MCParticleCollection& mcParticles) {
   // Create the collections
   auto colParticles = m_particles.createAndPut();
+  auto colParticles_trkCov = m_particles_trkCov.createAndPut();
   auto ascColParticlesToMC = m_mcAssociations.createAndPut();
 
   fcc::TaggedParticleCollection* colITags(nullptr);
@@ -50,6 +81,34 @@ StatusCode DelphesSaveChargedParticles::saveOutput(Delphes& delphes, const fcc::
     auto cand = static_cast<Candidate*>(delphesColl->At(j));
     auto particle = colParticles->create();
 
+
+    auto cov_d = CovToACTS(cand->TrackCovariance, cand->CtgTheta, cand->C);
+    auto t = colParticles_trkCov->create();
+    t.d0(1000* cand->D0);
+    t.z0(1000* cand->DZ);
+    t.phi(cand->Phi);
+   t.theta(TMath::ATan2(1.0, cand->CtgTheta));
+    Double_t b = -0.29988*2.0 / 2.; //mag field
+    t.qOverP(cand->C / (b*TMath::Sqrt(1 + cand->CtgTheta*cand->CtgTheta)));
+    std::array<float,15> t_c;
+    // save upper right triangle, row first
+    t_c[0] = cov_d[0][0];
+    t_c[1] = cov_d[0][1];
+    t_c[2] = cov_d[0][2];
+    t_c[3] = cov_d[0][3];
+    t_c[4] = cov_d[0][4];
+    t_c[5] = cov_d[1][1];
+    t_c[6] = cov_d[1][2];
+    t_c[7] = cov_d[1][3];
+    t_c[8] = cov_d[1][4];
+    t_c[9] = cov_d[2][2];
+    t_c[10] = cov_d[2][3];
+    t_c[11] = cov_d[2][4];
+    t_c[12] = cov_d[3][3];
+    t_c[13] = cov_d[3][4];
+    t_c[14] = cov_d[4][4];
+    t.cov(t_c);
+
     auto& barePart = particle.core();
     barePart.pdgId = cand->PID;
     barePart.status = cand->Status;

diff --git a/Sim/SimDelphesInterface/src/DelphesSaveChargedParticles.h b/Sim/SimDelphesInterface/src/DelphesSaveChargedParticles.h
@@ -8,9 +8,13 @@
 #include "FWCore/DataHandle.h"
 #include "SimDelphesInterface/IDelphesSaveOutputTool.h"
 
+
+
+
 // datamodel
 namespace fcc {
 class ParticleCollection;
+class TrackStateCollection;
 class ParticleMCParticleAssociationCollection;
 class TaggedParticleCollection;
 }
@@ -46,6 +50,7 @@ class DelphesSaveChargedParticles : public GaudiTool, virtual public IDelphesSav
 
 private:
   /// Handle the particles to be saved
+  DataHandle<fcc::TrackStateCollection> m_particles_trkCov{"particles_trkCov", Gaudi::DataHandle::Writer, this};
   DataHandle<fcc::ParticleCollection> m_particles{"particles", Gaudi::DataHandle::Writer, this};
   /// Handle to associate particles with MCParticles
   DataHandle<fcc::ParticleMCParticleAssociationCollection> m_mcAssociations{"mcAssociations", Gaudi::DataHandle::Writer,