Skip to content

Conversation

@ctdax
Copy link
Contributor

@ctdax ctdax commented Nov 18, 2025

PR description:

This pull request introduces the ability to decay R-hadrons during simulation using a side instance of Pythia8. A generator is required to decay the underlying SUSY particle to hadronize the outgoing partons. This will allow users to simulate detector interactions of long-lived R-hadrons with proper decays for the first time.

See the talk given here.

How to use it
  • Create a SLHA decay table that includes the decays and widths of the underlying SUSY particle in the R-hadron, e.g. the gluino or the stop.
  • Pass the SLHA decay table to SLHAFile inside of process.generator in the python config file. Optionally, commands can be passed directly to the instance of pythia that handles the R-hadron decay via RhadronPythiaDecayerCommandFile, e.g.
process.generator = cms.EDFilter("Pythia8ConcurrentGeneratorFilter",
    SLHAFile = cms.string('SimG4Core/CustomPhysics/data/TESTDECAY_GLUINO1800_STOP500.txt'),
    RhadronPythiaDecayerCommandFile = cms.untracked.string('SimG4Core/CustomPhysics/data/RhadronPythiaDecayerCommands.txt'),
...
)
  • Make sure the event generator doesn't decay the SUSY particle before it is handed to Geant for simulation. If the event generator is pythia, you can force the SUSY particle to not decay via pythia8CommonSettings in the python config file:
process.generator = cms.EDFilter("Pythia8ConcurrentGeneratorFilter",
    PythiaParameters = cms.PSet(
        pythia8CommonSettings = cms.vstring(
            '1000021:mayDecay = off'
        )
    )
)
Explanation of changes
  • RHadronPythiaDecayer.cc/RHadronPythiaDecayer.h: Utilizes G4Decay::DecayIt() to decay the R-hadron. Because it is defined as an external decayer, when G4Decay::DecayIt() calls ImportDecayProducts() it calls the function inside of RHadronPythiaDecayer. This file strips the incoming R-hadron down to its constituents, decays the underlying SUSY particle, then lets Pythia handle the hadronization. The decay products are returned to G4Decay::DecayIt() for them to be simulated in Geant.
  • RHadronPythiaDecayDataManager.h: Holds the decay information during simulation to be passed to RHDecayTracer.cc during the EDM step.
  • RHDecayTracer.cc/RHDecayTracer.h: Adds the R-hadron decay vertex to the Gen collection. This will allow users to view all of the relevant information regarding the decay, such as the incoming and outgoing particles, after the simulation ends.
  • SimG4Core/CustomPhysics/BuildFile.xml: Pythia8Interface is required to generate the R-hadron decays in Pythia
  • SimG4Core/CustomPhysics/plugins/module.cc: Define the RHDecayTracer framework module
  • SimG4Core/CustomPhysics/python/Exotica_HSCP_SIM_cfi.py: Handles the inputs SLHAFile and RhadronPythiaDecayerCommandFile. Adds the RHDecayTracer module to the process.
  • SimG4Core/CustomPhysics/src/CustomParticleFactory.cc: Forces the lifetime of R-hadrons to match the lifetime of their underlying SUSY particle, which is read in from the SLHA table. The modification to Line298 if ((ToLower(line).find("block") < line.npos) || (line == "#")) { forces the skipping of commented lines in the SLHA text file.
  • SimG4Core/CustomPhysics/src/CustomPhysicsList.cc: Sets the R-hadron decay process to be the external pythia decayer. Lines145-159 is a bug fix. The process expects transportation to happen last, and multiple scattering to happen second to last (see slide 16 of this Geant talk); this was not previously the case. Before the fix, the lifetime of the R-hadrons did not match the expected exponential. The correction to the code fixed the issue. See the two plots below:
  • SimG4Core/CustomPhysics/data/RhadronPythiaDecayerCommands.txt: A new .txt file that contains default commands for the Pythia instance that decays the Rhadrons.
  • SimG4Core/CustomPhysics/python/CustomPhysics_cfi.py: Adds the RhadronPythiaDecayerCommandFile parameter to the config file.

Before the bug fix
buggy_rhadron-lifetime

After the bug fix (the red and the orange overlap)
corrected_rhadron-lifetime

Confirmation of 4-momentum conservation at the decay

A histogram plotting the difference in 4-momentum between the R-hadron and its decay products in 50 events is shown below. These are the 4-momenta in the $p_T$, $\eta$, $\phi$ basis after the decay products are handed back to Geant. All entries are consistent with 0 difference except for some floating-point error.

4pDifferenceAtDecay

PR validation:

All tests below ran successfully.

cd src
scram b distclean 
git cms-checkdeps -a -A
scram b -j 8
scram b runtests use-ibeos
runTheMatrix.py -l limited -i all --ibeos

Before submitting your pull requests, make sure you followed this checklist:

@cmsbuild
Copy link
Contributor

cmsbuild commented Nov 18, 2025

cms-bot internal usage

@cmsbuild
Copy link
Contributor

-code-checks

Logs: https://cmssdt.cern.ch/SDT/code-checks/cms-sw-PR-49414/46868

Code check has found code style and quality issues which could be resolved by applying following patch(s)

@cmsbuild
Copy link
Contributor

@cmsbuild
Copy link
Contributor

A new Pull Request was created by @ctdax for master.

It involves the following packages:

  • SimG4Core/CustomPhysics (simulation)

@civanch, @cmsbuild, @kpedro88, @mdhildreth can you please review it and eventually sign? Thanks.
@bsunanda, @fabiocos, @makortel, @martinamalberti, @rovere, @slomeo this is something you requested to watch as well.
@ftenchini, @mandrenguyen, @sextonkennedy you are the release manager for this.

cms-bot commands are listed here

@kpedro88
Copy link
Contributor

@ctdax thanks for this contribution! It will require some review. Are you available to present this work at Friday's SIM meeting?

@ctdax
Copy link
Contributor Author

ctdax commented Nov 18, 2025

@ctdax thanks for this contribution! It will require some review. Are you available to present this work at Friday's SIM meeting?

Of course! Where can I find the indico for the meeting?

@kpedro88
Copy link
Contributor

I'll send it to you separately

}

void RHadronPythiaDecayer::RHadronToConstituents(Pythia8::Event& event) {
// This code is very similar to Pythia8::RHadrons::decay(). Unfortunately, it is not possible in this scenario to use Pythia8::RHadrons::decay().
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, unfortunately nRHad is a private member in Pythia8::RHadrons. If it were protected, one could make an inherited class and modify its value manually. We could propose to change this in the future.

class HepMCProduct;
}

class RHDecayTracer : public edm::one::EDProducer<edm::one::SharedResources> {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand the purpose of this class. Why is it a one producer? Why is it a producer at all, if it does not produce anything?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have changed the class from a one producer to edm::stream::EDProducer<>. I am also now implementing the change that you suggested during the SIM meeting last Friday. Before making commits however I would like to make sure I understood your suggestion. You believe that we should not modify edmHepMCProduct_generatorSmeared, but rather we should produce a new edmHepMCProduct that includes the R-hadron decay vertex information. Should the new edmHepMCProduct:

  • A) be a copy of edmHepMCProduct_generatorSmeared that also includes the R-hadron decay vertex?
  • B) be a new edmHepMCProduct that only includes the vertex?
  • or C) be something different?

@kpedro88
Copy link
Contributor

@ctdax it depends on how this information will be used.

If you want it to be straightforward downstream to compute GEN-level objects that use the R-hadron decay information (GenJets, GenMET, mini/nanoAOD storage, etc.), then copying and updating the HepMCProduct is the best option. This will require moving around some of the initial GEN-SIM sequences, but that's a smaller attack surface than trying to propagate multiple separate collections through all the downstream steps.

If you only need the R-hadron decay information for debugging or spot tests, then making a new collection that only includes the decay vertex will be much easier.

@ctdax
Copy link
Contributor Author

ctdax commented Dec 2, 2025

@kpedro88 we decided that we would like to compute GEN-level objects using the decay info.

Should the new HepMCProduct be a shallow copy or a deep copy of the original HepMCProduct?

@kpedro88
Copy link
Contributor

kpedro88 commented Dec 2, 2025

It will need to be a deep copy. I am not aware of any way to have a shallow copy with extra particles/information added on top, and even if it existed, using it would likely require invasive changes to many algorithms/producers.

You will need to modify sequences such that the existing HepMCProduct has a different name, and then the new HepMCProduct uses the existing name, such that downstream modules consume it automatically.

@cmsbuild
Copy link
Contributor

Milestone for this pull request has been moved to CMSSW_16_1_X. Please open a backport if it should also go in to CMSSW_16_0_X.

@cmsbuild cmsbuild modified the milestones: CMSSW_16_0_X, CMSSW_16_1_X Dec 18, 2025
@cmsbuild cmsbuild modified the milestones: CMSSW_16_0_X, CMSSW_16_1_X Dec 18, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants