15_0_X Added Pythia-based Rhadron decays during Geant simulation

PhysicsTools/HepMCCandAlgos/plugins/GenPlusSimParticleProducer.cc

@@ -264,52 +264,88 @@ void GenPlusSimParticleProducer::produce(Event &event, const EventSetup &iSetup)
   }
 
   for (SimTrackContainer::const_iterator isimtrk = simtracks->begin(); isimtrk != simtracks->end(); ++isimtrk) {
-    // Skip PYTHIA tracks.
-    if (isimtrk->genpartIndex() != -1)
-      continue;
+    // Special case for when simTracks need to be added to pythia track filtered by particleType, e.g. a particle with ctau=0.1mm
+    if (isimtrk->genpartIndex() != -1 && !pdgIds_.empty() &&
+        (pdgIds_.find(std::abs(isimtrk->type())) != pdgIds_.end())) {
+      std::vector<int>::const_iterator itIndex;
+      if (barcodesAreSorted) {
+        itIndex = std::lower_bound(genBarcodes->begin(), genBarcodes->end(), isimtrk->genpartIndex());
+      } else {
+        itIndex = std::find(genBarcodes->begin(), genBarcodes->end(), isimtrk->genpartIndex());
+      }
 
-    // Maybe apply the PdgId filter
-    if (!pdgIds_.empty()) {  // if we have a filter on pdg ids
-      if (pdgIds_.find(std::abs(isimtrk->type())) == pdgIds_.end())
-        continue;
+      if ((itIndex != genBarcodes->end()) && (*itIndex == isimtrk->genpartIndex())) {
+        const unsigned int momidx = itIndex - genBarcodes->begin();
+        for (SimTrackContainer::const_iterator idau = simtracksSorted->begin(); idau != simtracksSorted->end();
+             ++idau) {
+          if (idau->noVertex())
+            continue;
+          const SimVertex &dvtx = (*simvertices)[idau->vertIndex()];
+          if (dvtx.noParent())
+            continue;
+          if (dvtx.parentIndex() == static_cast<int>(isimtrk->trackId())) {
+            addGenParticle(*isimtrk,
+                           *idau,
+                           momidx,
+                           *simtracksSorted,
+                           *simvertices,
+                           *candsPtr,
+                           ref,
+                           *newGenBarcodes,
+                           barcodesAreSorted);
+          }
+        }
+      }
     }
 
-    // find simtrack that has a genParticle match to its parent
-    // Look at the production vertex. If there is no vertex, I can do nothing...
-    if (!isimtrk->noVertex()) {
-      // Pick the vertex (isimtrk.vertIndex() is really an index)
-      const SimVertex &vtx = (*simvertices)[isimtrk->vertIndex()];
+    else {
+      // Skip PYTHIA tracks.
+      if (isimtrk->genpartIndex() != -1)
+        continue;
 
-      // Check if the vertex has a parent track (otherwise, we're lost)
-      if (!vtx.noParent()) {
-        // Now note that vtx.parentIndex() is NOT an index, it's a track id, so I have to search for it
-        unsigned int idx = vtx.parentIndex();
-        SimTrackContainer::const_iterator it =
-            std::lower_bound(simtracksSorted->begin(), simtracksSorted->end(), idx, LessById());
-        if ((it != simtracksSorted->end()) && (it->trackId() == idx)) {  //it is the parent sim track
-          if (it->genpartIndex() != -1) {
-            std::vector<int>::const_iterator itIndex;
-            if (barcodesAreSorted) {
-              itIndex = std::lower_bound(genBarcodes->begin(), genBarcodes->end(), it->genpartIndex());
-            } else {
-              itIndex = std::find(genBarcodes->begin(), genBarcodes->end(), it->genpartIndex());
-            }
+      // Maybe apply the PdgId filter
+      if (!pdgIds_.empty()) {  // if we have a filter on pdg ids
+        if (pdgIds_.find(std::abs(isimtrk->type())) == pdgIds_.end())
+          continue;
+      }
 
-            // Check that I found something
-            // I need to check '*itIndex == it->genpartIndex()' because lower_bound just finds the right spot for an item in a sorted list, not the item
-            if ((itIndex != genBarcodes->end()) && (*itIndex == it->genpartIndex())) {
-              // Ok, I'll make the genParticle for st and add it to the new collection updating the map and parent-child relationship
-              // pass the mother and daughter sim tracks and the mother genParticle to method to create the daughter genParticle and recur
-              unsigned int momidx = itIndex - genBarcodes->begin();
-              addGenParticle(*it,
-                             *isimtrk,
-                             momidx,
-                             *simtracksSorted,
-                             *simvertices,
-                             *candsPtr,
-                             ref,
-                             *newGenBarcodes,
-                             barcodesAreSorted);
+      // find simtrack that has a genParticle match to its parent
+      // Look at the production vertex. If there is no vertex, I can do nothing...
+      if (!isimtrk->noVertex()) {
+        // Pick the vertex (isimtrk.vertIndex() is really an index)
+        const SimVertex &vtx = (*simvertices)[isimtrk->vertIndex()];
+
+        // Check if the vertex has a parent track (otherwise, we're lost)
+        if (!vtx.noParent()) {
+          // Now note that vtx.parentIndex() is NOT an index, it's a track id, so I have to search for it
+          unsigned int idx = vtx.parentIndex();
+          SimTrackContainer::const_iterator it =
+              std::lower_bound(simtracksSorted->begin(), simtracksSorted->end(), idx, LessById());
+          if ((it != simtracksSorted->end()) && (it->trackId() == idx)) {  //it is the parent sim track
+            if (it->genpartIndex() != -1) {
+              std::vector<int>::const_iterator itIndex;
+              if (barcodesAreSorted) {
+                itIndex = std::lower_bound(genBarcodes->begin(), genBarcodes->end(), it->genpartIndex());
+              } else {
+                itIndex = std::find(genBarcodes->begin(), genBarcodes->end(), it->genpartIndex());
+              }
+
+              // Check that I found something
+              // I need to check '*itIndex == it->genpartIndex()' because lower_bound just finds the right spot for an item in a sorted list, not the item
+              if ((itIndex != genBarcodes->end()) && (*itIndex == it->genpartIndex())) {
+                // Ok, I'll make the genParticle for st and add it to the new collection updating the map and parent-child relationship
+                // pass the mother and daughter sim tracks and the mother genParticle to method to create the daughter genParticle and recur
+                unsigned int momidx = itIndex - genBarcodes->begin();
+                addGenParticle(*it,
+                               *isimtrk,
+                               momidx,
+                               *simtracksSorted,
+                               *simvertices,
+                               *candsPtr,
+                               ref,
+                               *newGenBarcodes,
+                               barcodesAreSorted);
+              }
             }
           }
         }

SimG4Core/CustomPhysics/BuildFile.xml

@@ -7,6 +7,7 @@
 <use name="SimG4Core/Physics"/>
 <use name="SimG4Core/PhysicsLists"/>
 <use name="SimG4Core/Watcher"/>
+<use name="GeneratorInterface/Pythia8Interface"/>
 <use name="geant4core"/>
 <use name="clhep"/>
 <use name="boost"/>

SimG4Core/CustomPhysics/data/RhadronPythiaDecayerCommands.txt

@@ -0,0 +1,7 @@
+ProcessLevel:all = off
+SUSY:all = on
+RHadrons:allow = off
+1000021:mWidth = 1000.0
+1000006:mWidth = 1000.0
+RHadrons:probGluinoball = 0.1
+PartonLevel:FSR = off

SimG4Core/CustomPhysics/interface/RHadronPythiaDecayer.h

@@ -0,0 +1,59 @@
+#ifndef SimG4Core_CustomPhysics_RHadronPythiaDecayer_H
+#define SimG4Core_CustomPhysics_RHadronPythiaDecayer_H
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "G4Decay.hh"
+#include "G4VExtDecayer.hh"
+#include "G4ThreeVector.hh"
+#include <string>
+#include <vector>
+#include <utility>
+#include <memory>
+
+namespace gen {
+  class P8RndmEngine;
+  typedef std::shared_ptr<P8RndmEngine> P8RndmEnginePtr;
+}  // namespace gen
+
+namespace Pythia8 {
+  class Pythia;
+  class Event;
+  class Rndm;
+}  // namespace Pythia8
+
+class G4DynamicParticle;
+class G4DecayProducts;
+class G4VParticleChange;
+class G4Step;
+class G4Track;
+class RHadronPythiaDecayer : public G4Decay, public G4VExtDecayer {
+public:
+  RHadronPythiaDecayer(edm::ParameterSet const& p);
+  ~RHadronPythiaDecayer() override;
+
+  //What Geant calls to decay the Rhadron
+  G4VParticleChange* DecayIt(const G4Track& aTrack, const G4Step& aStep) override;
+
+  //Tell pythia to decay the Rhadron and return the products in Geant format
+  G4DecayProducts* ImportDecayProducts(const G4Track&) override;
+
+private:
+  // Strip the RHadron down to its constituents in preperation for decaying the gluino or squark
+  void RHadronToConstituents(Pythia8::Event& event);
+
+  std::pair<int, int> fromIdWithSquark(int idRHad) const;
+  std::pair<int, int> fromIdWithGluino(int idRHad, Pythia8::Rndm* rndmPtr) const;
+
+  bool isGluinoRHadron(int pdgId) const;
+
+  //Fill a Pythia8 event with the information from a G4Track
+  void fillParticle(const G4Track&, Pythia8::Event& event) const;
+
+  //Function to decay the RHadron and return products in G4 format
+  void pythiaDecay(const G4Track&, std::vector<G4DynamicParticle*>&);
+
+  std::unique_ptr<Pythia8::Pythia> pythia_;
+  std::vector<G4ThreeVector> secondaryDisplacements_;
+};
+
+#endif

SimG4Core/CustomPhysics/python/CustomPhysics_cfi.py

@@ -11,6 +11,9 @@
     processesDef = cms.FileInPath('SimG4Core/CustomPhysics/data/RhadronProcessList.txt'),
     amplitude = cms.double(100.0), ##
 
+    # Path to pythia arguments for rhadron decay
+    RhadronPythiaDecayerCommandFile = cms.FileInPath('SimG4Core/CustomPhysics/data/RhadronPythiaDecayerCommands.txt'),
+
     # R-hadron physics setup
     rhadronPhysics = cms.bool(True),
     resonant = cms.bool(False),

SimG4Core/CustomPhysics/python/Exotica_HSCP_SIM_cfi.py

@@ -1,4 +1,7 @@
 import FWCore.ParameterSet.Config as cms
+from SimG4Core.CustomPhysics.GenPlusSimParticles_cfi import customizeProduce, customizeKeep
+import os
+import re
 
 def customise(process):
 
@@ -11,6 +14,18 @@ def customise(process):
     process.customPhysicsSetup.particlesDef = PARTICLE_FILE
     process.customPhysicsSetup.reggeModel = USE_REGGE
 
+    # Read in the SLHA file to get the particles definition
+    if hasattr(process, 'generator') and hasattr(process.generator, 'SLHAFileForPythia8'):
+        process.customPhysicsSetup.particlesDef = process.generator.SLHAFileForPythia8.value()
+    elif hasattr(process, 'g4SimHits') and hasattr(process.g4SimHits, 'SLHAFileForPythia8'):
+        process.customPhysicsSetup.particlesDef = process.g4SimHits.SLHAFileForPythia8.value()
+
+    # Passing pythia settings to the Rhadron decayer is optional
+    if hasattr(process, 'generator') and hasattr(process.generator, 'RhadronPythiaDecayerCommandFile'):
+        process.customPhysicsSetup.RhadronPythiaDecayerCommandFile = process.generator.RhadronPythiaDecayerCommandFile.value()
+    elif hasattr(process, 'g4SimHits') and hasattr(process.g4SimHits, 'RhadronPythiaDecayerCommandFile'):
+        process.customPhysicsSetup.RhadronPythiaDecayerCommandFile = process.g4SimHits.RhadronPythiaDecayerCommandFile.value()
+
     if hasattr(process,'g4SimHits'):
         # defined watches
         process.g4SimHits.Watchers = cms.VPSet (
@@ -42,5 +57,49 @@ def customise(process):
             process.customPhysicsSetup
         )	
 
+        # Add customize produce to add genParticlePlusGeant collection. Make sure particleTypes are only the HSCPs
+        if hasattr(process, 'generator') and hasattr(process.generator, 'pdtFile'):
+            process.HepPDTESSource.pdtFileName = process.generator.pdtFile.value()
+        elif hasattr(process, 'g4SimHits') and hasattr(process.g4SimHits, 'pdtFile'):
+            process.HepPDTESSource.pdtFileName = process.g4SimHits.pdtFile.value()
+
+        process = customizeProduce(process)
+        hscp_particle_names = get_hscp_particle_names_from_pdt(FLAVOR, process.HepPDTESSource.pdtFileName.value())
+        process.genParticlePlusGeant.particleTypes = cms.vstring(hscp_particle_names)
+        process.genParticlePlusGeant.filter = cms.vstring("")
+        process = customizeKeep(process)
+
         return (process)
 
+
+def get_hscp_particle_names_from_pdt(flavor, pdt_file_path):
+    flavor_to_pdg_prefix = {
+        "gluino": "~g_",
+        "stop": "~T",
+        "stau": "~tau",
+    }
+    pdg_prefix = flavor_to_pdg_prefix.get(flavor)
+    if not pdg_prefix:
+        return []
+
+    hscp_names = []
+    seen = set()
+
+    if not os.path.isabs(pdt_file_path):
+        scram_arch = os.environ.get('SCRAM_ARCH', '')
+        pdt_file_path = os.path.join(os.environ.get('CMSSW_RELEASE_BASE', ''), 'external', scram_arch, 'data', pdt_file_path)
+
+    with open(pdt_file_path, 'r') as handle:
+        for raw_line in handle:
+            line = raw_line.strip()
+            if not line or line.startswith('#') or line.startswith('//'):
+                continue
+            cols = line.split()
+            if len(cols) < 2:
+                continue
+            name = cols[1]
+            if name.startswith(pdg_prefix) and name not in seen:
+                seen.add(name)
+                hscp_names.append(name)
+
+    return hscp_names