Add Delphes isolation variable and ParticleFlowCandidates to edm4hep …

…output (#107)
key4hep · Apr 21, 2023 · 430f6e4 · 430f6e4
1 parent 5f136a3
commit 430f6e4
Show file tree

Hide file tree

Showing 3 changed files with 58 additions and 11 deletions.
diff --git a/converter/include/k4SimDelphes/DelphesEDM4HepConverter.h b/converter/include/k4SimDelphes/DelphesEDM4HepConverter.h
@@ -42,6 +42,11 @@ namespace k4SimDelphes {
    */
   constexpr std::array<std::string_view, 1> RECO_CLUSTER_OUTPUT = {"Tower"};
 
+  /**
+   * Classes that will be stored as particle flow candidates
+   */
+  constexpr std::array<std::string_view, 1> RECO_CANDIDATES_OUTPUT = {"ParticleFlowCandidate"};
+
   /**
    * Classes that will be stored as TrackerHits
    */
@@ -53,7 +58,7 @@ namespace k4SimDelphes {
   constexpr auto CALORIMETERHIT_OUTPUT_NAME = "CalorimeterHits";
 
   /**
-   * * Eventheader class will be stored only
+   * * Eventheader class will be stored only once
    */
   constexpr auto EVENTHEADER_NAME = "EventHeader";
 
@@ -105,6 +110,7 @@ namespace k4SimDelphes {
 
     void processParticles(const TClonesArray* delphesCollection, std::string const& branch);
     void processTracks(const TClonesArray* delphesCollection, std::string const& branch);
+    void processPFlowCandidates(const TClonesArray* delphesCollection, std::string const& branch);
     void processClusters(const TClonesArray* delphesCollection, std::string const& branch);
     void processJets(const TClonesArray* delphesCollection, std::string const& branch);
     void processPhotons(const TClonesArray* delphesCollection, std::string const& branch) {

diff --git a/converter/src/DelphesEDM4HepConverter.cc b/converter/src/DelphesEDM4HepConverter.cc
@@ -39,8 +39,9 @@ namespace k4SimDelphes {
    * NOTE: not a configuration parameter. this has to be done in this order to
    * ensure that products required by later stages are producd early enough
    */
-  constexpr std::array<std::string_view, 9> PROCESSING_ORDER = {
-      "GenParticle", "Track", "Tower", "Muon", "Electron", "Photon", "Jet", "MissingET", "SclalarHT"};
+  constexpr std::array<std::string_view, 10> PROCESSING_ORDER = {
+      "GenParticle", "Track",     "Tower",    "ParticleFlowCandidate", "Muon", "Electron", "Photon",
+      "Jet",         "MissingET", "SclalarHT"};
 
   template <size_t N>
   void sortBranchesProcessingOrder(std::vector<BranchSettings>&           branches,
@@ -97,6 +98,11 @@ namespace k4SimDelphes {
         m_processFunctions.emplace(branch.name, &DelphesEDM4HepConverter::processTracks);
       }
 
+      if (contains(outputSettings.ReconstructedParticleCollections, branch.name.c_str()) &&
+          contains(RECO_CANDIDATES_OUTPUT, branch.className.c_str())) {
+        m_processFunctions.emplace(branch.name, &DelphesEDM4HepConverter::processPFlowCandidates);
+      }
+
       if (contains(outputSettings.ReconstructedParticleCollections, branch.name.c_str()) &&
           contains(RECO_CLUSTER_OUTPUT, branch.className.c_str())) {
         m_processFunctions.emplace(branch.name, &DelphesEDM4HepConverter::processClusters);
@@ -371,11 +377,30 @@ namespace k4SimDelphes {
     }
   }
 
+  void DelphesEDM4HepConverter::processPFlowCandidates(const TClonesArray* delphesCollection,
+                                                       std::string const&  branch) {
+    auto* candidateCollection = createCollection<edm4hep::ReconstructedParticleCollection>(branch);
+
+    for (auto iCand = 0; iCand < delphesCollection->GetEntries(); ++iCand) {
+      auto* delphesCand = static_cast<ParticleFlowCandidate*>(delphesCollection->At(iCand));
+      auto  candidate   = candidateCollection->create();
+
+      candidate.setCharge(delphesCand->Charge);
+      candidate.setMass(delphesCand->Mass);
+      const auto momentum = delphesCand->P4();
+      candidate.setEnergy(momentum.E());
+      candidate.setMomentum({(float)momentum.Px(), (float)momentum.Py(), (float)momentum.Pz()});
+    }
+  }
+
   template <typename DelphesT>
   void DelphesEDM4HepConverter::fillReferenceCollection(const TClonesArray* delphesCollection,
                                                         std::string const& branch, std::string_view const type) {
     auto* collection = createCollection<edm4hep::ReconstructedParticleCollection>(branch, true);
 
+    //add collection for the isolation variable calculated by delphes:
+    auto* isoIDColl = createCollection<podio::UserDataCollection<float>>(std::string(branch) + "_IsolationVar");
+
     for (auto iCand = 0; iCand < delphesCollection->GetEntries(); ++iCand) {
       auto* delphesCand = static_cast<DelphesT*>(delphesCollection->At(iCand));
 
@@ -393,6 +418,9 @@ namespace k4SimDelphes {
           matchedReco->setCharge(delphesCand->Charge);
         }
 
+        //fill the isolation var
+        isoIDColl->push_back(delphesCand->IsolationVar);
+
       } else {
         std::cerr << "**** WARNING: No matching ReconstructedParticle was found for a Delphes " << type << std::endl;
       }

diff --git a/doc/output_config.md b/doc/output_config.md
@@ -8,7 +8,7 @@ that defines the necessary parameters:
 
 [edm4hep_output_config.tcl](https://github.com/key4hep/k4SimDelphes/blob/main/examples/edm4hep_output_config.tcl). This `tcl` file has to
 be passed to the conversion executable as parameter from the command line. **No
-changes to the delphes card that is used are necessary.** An exemplary call
+changes to the delphes card that is used are necessary.**[^*] An exemplary call
 could look something like this (replacing the `delphes_card.tcl` and the
 `input_file.stdhep` with actual input files):
 
@@ -24,6 +24,8 @@ easily reuse the `tcl` parser from delphes and to introduce a logical
 containment for the parameters. The names of the delphes branches are taken from the `TreeWriter`
 module defintion and used in the `EDM4HepOutput`.
 
+[^*]: As long as all collections requested for the conversion to `edm4hep` are defined as output branches in the `TreeWriter` section of the Delphes card. If not, they need to be added there.
+
 ## Collection conversions
 
 The configuration mainly steers which Delphes collections are still available in
@@ -41,6 +43,9 @@ the assumption that this is a non-overlapping list of particles. It is the users
 responsibility to make sure that this is the case. (See [known
 issues](#known-issues)).
 
+`ParticleFlowCandidate` collections from Delphes can also be stored as separate `ReconstructedParticleCollection`s in the output, but they currently have no associations to the generated particles. Note that they are **not** added to the global `ReconstructedParticleCollection` described above, as that would lead to double counting.
+
+
 ### `GenParticleCollection`
 All Delphes `GenParticle` collections that will be considered and stored as
 `edm4hep::MCParticle`. Each Delphes collection will be put into its own
@@ -84,7 +89,7 @@ The filled collection contains only one element per event.
 
 The parameters **`RecoParticleCollectionName`** and
 **`MCRecoAssociationCollectionName`** control the names of the [global
-reconstructed particle colleciton](#reconstructedparticlecollections) and the
+reconstructed particle collection](#reconstructedparticlecollections) and the
 collection with the `MCRecoParticleAssociation`s that can be used to find the
 `MCParticle`s associated to `ReconstructedParticle`s (and vice versa).
 
@@ -105,26 +110,34 @@ classes. For the conversion the Delphes classes are taken from the `TreeWriter`
 | `Photon`      | `ReconstructedParticle` (subset collection)       |
 | `MissingET`   | `ReconstructedParticle`                           |
 | `ScalarHT`    | `ParticleID`                                      |
+| `ParticleFlowCandidate`       | `ReconstructedParticle`           |
 | n/a           | `MCRecoParticleAssociation`                       |
 
 All Delphes classes that are not listed here are currently not converted.
 
+## EventHeader Collection
+The `EventHeader` collection is used to store information from the Delphes `Event` classes. It always contains one element pre event with member variables `eventNumber` and `weight`. This is currently only implemented for the `DelphesPythia8Reader`.
+
 
-### Known issues
+## Isolation variable
+The isolation variable as calculated by the Delphes isolation module, will be added as a `UserDataCollection` for each collection of `Electron`, `Muon` and `Photon` during the conversion, called `<collection_name>_IsolationVar`.
 
-- [ ] Double counting of Tracks and Clusters. In Delphes it is possible that a
+
+## Known issues
+
+<!-- - [ ] Double counting of Tracks and Clusters. In Delphes it is possible that a
       `Tower` and a `Track` point back to the same generated particle. This is
       currently not checked and each `Track` and `Tower` will be converted into
       an `edm4hep::ReconstructedParticle`. Hence, it is possible to get more
       than one per generated particle. This means that the number of jet
       constituents will be differend between Jets in the Delphes output and in
-      the `edm4hep` output.
+      the `edm4hep` output. -> solved as of April 2023? -->
 
 - [ ] Not all available information is used in the conversion. An incomplete list
       of things that are currently not available in `edm4hep`:
   - [ ] Jet substructure variables (including subjets)
-  - [ ] Isolation variables
-  - [ ] Flavor tag information
-  - [ ] Tau tag information
+
+- [ ] No conversion of `EventHeader` information for readers other than `DelphesPythia8Reader`
+