Pythia8 afterburner

[jncots]

This should fix #87 mentioned in #165 and some other problems releated to issue that some event generators cannot decay some particles.

[jncots]

About execution time:

DecayAfterburner.__call__(event) typically adds the same execution time as it is required for event generation (~ 1 ms), i.e. the time increases twice. For DpmJet:

37% of the time takes decay itself (forceHadronLevel(), _decay on picture), 31 % takes extraction of event stack from Pythia8 (_extract_result on picture), 10% takes filling of Pythia8 stack (_fill), 6% np.copy and np.where, 15 % - copying initial stack to final one.

[jncots]

Problems noticed during testing:

• DpmjetIII193 produces particles with pid = 99999, Dpmjet306 produces particles with pid = 80000, UrQMD34 produces particles with pid = 0. In case of Dpmjet it is not a problem because status != 1, but UrQMD34 produces pid = 0 with status = 1. However, as charge is returned by Pythia8, it is a problem for afterburner, because Pythia8 doesn't know pid with 0, 99999, or 80000. I changed the charge_from_pid to overcome this problem.

• UrQMD34 returns neutron with wrong masses, so afterburner cannot decay them. The same is for QGSJet, but in this case neutrons are returned with proton mass.

[jncots]

I started to integrate DecayAutoburner with models. I changed interface to provide more flexible settings for decaying/stable particles. "theonly_decaying_pids" sets all particles as stable, except the one that in the list. "theonly_stable_pids" sets all particles as unstable, except the one that in the list. The arguments are convenient to set an initial state of stable/unstable in Pythia. "extra_decaying/stable_pids" toggle this state further.

I noticed that no model has expicit lists of stable/unstable particles. Also, there is no common interface to inquire event generators what they consider stable/unstable. Currently, I set:

self._afterburner = DecayAfterburner(
                seed=seed, theonly_stable_pids=long_lived
            )

,i.e. all other particles that is not in long_lived list will decay by DecayAfterburner.

I guess, it is just enough to keep final_state_particles list for each instance and implicitly assume that all other should decay. Then set_stable should modify this list. final_state_particles could be a property which call set_stable if a new list is assigned to it.

The easiest way to switch on DecayAfterburner for a model is to provide a common method switch_on_autoburner, which will call:

self._afterburner = DecayAfterburner(
                seed=seed, theonly_stable_pids=self.final_state_particles
            )

[HDembinski]

I am unhappy with the naming.

1) I would prefer a more serious name than "DecayAfterburner". I am open for suggestions. Since Pythia8 is doing the work, should we call it one these?

• Pythia8Decay
• Pythia8DecayProcessor
• Pythia8DecayHandler

2) The keyword theonly_stable_pids breaks rules regarding underscores. Either the_only_stable_pids or just stable_pids. Keywords shouldn't be too long and don't have to be self-explanatory. You have to write documentation anyway to explain what the keyword does.

Pythia8Decay is primarily a component for us the library developers. It should be used automatically internally by some model that cannot internally decay all particles. It should only be run if the user requests that particles should decay which the model cannot decay by itself. All this should happen without any special user input. The user just requests which particles should be stable/decay with set_stable.

I noticed that no model has expicit lists of stable/unstable particles. Also, there is no common interface to inquire event generators what they consider stable/unstable. Currently, I set:

We already have an interface to set particles stable with set_stable, but as you say we don't have a way to query which particles are currently set as stable in the model.

I like the idea of replacing this interface with a property on the model called final_state_particles which returns a list of pids that are considered final state by the model. Assigning a list to final_state_particles should be allowed and should set all those particles stable and decay all others. We can then deprecate set_stable which is a bit of a misnomer.

[jncots]

I added the functionality discussed above. By default DecayAfterburner is not used. To switch on/off one should explicitly call switch_afterburner(on=True/False).

_final_state_particles is the list containing list of the pdg ids which are considered as final. They can be obtained, but not directly assigned. set_stable is used to do this. The default value is long_lived.

gen = chromo.models.DpmjetIII193(evt_kin)
gen.switch_afterburner()
for event in gen(10):
     ...

# switch off afterburner
gen.switch_afterburner(on=False)

[jncots]

I am unhappy with the naming.

Ok, I like Pythia8DecayHandler. stable_pids is fine.

It should only be run if the user requests that particles should decay which the model cannot decay by itself.

Ideally, yes. Is there any good way to get those particles?

Assigning a list to final_state_particles should be allowed and should set all those particles stable and decay all others.

Basically, it is done in _update_final_state_particles but currently it is private method. I find it more covenient to be a method instead of property setter because it allows add and remove particles from the list via second parameter (stable/unstable).

[HDembinski]

Ideally, yes. Is there any good way to get those particles?

It is model specific. To my knowledge, only QGSJet does not decay anything. I learned from your comments that there are also some issues with DPMJet. Most of the other models do not need this functionality.

I don't like the current design.

• Pythia8DecayHandler should not be called by MCRun. It is not a generic component that all models need. It is a special component that should be used in specific model implementations, which do not implement decays properly, like the class QGSJet in the qgsjet module.
• If you find that there is common code which all such models can share that require external decay handling, then make a new subclass of MCRun called MCRunWithExternalDecay of something, which derives from MCRun. The models like QGSJet then inherit from MCRunWithExternalDecay while all the other models that don't need this functionality derive from MCRun. This is OOP.
• Please add a public property called final_state_particles to MCRun which calls abstract private implementations (e.g. _set_final_state_particles and _get_final_state_particles) that all concrete models implement. The property should return a tuple of particle ids that the model currently includes in the final state. For most models, this is the sequence of long-lived pids that we have in constants by default, unless the user changes something. When assigned to, this property will set which particles are included in the final state by the model. This can be easily implemented by calling the hidden model-specific method _set_stable. It is important to return a tuple instead of a list, to avoid UX issues. A tuple is read-only, which prevents misusage of the interface, e.g. my_model.final_state_particles.append(111) would not change the internal state of the model, contrary to user expectation. By returning a tuple, this cannot happen.
• Remove switch_afterburner. This should not be toggled manually. The model decides on its own if needs the Pythia8DecayHandler based on the requests that the user made to the attribute final_state_particles.

I find it more covenient to be a method instead of property setter because it allows add and remove particles from the list via second parameter (stable/unstable).

Allowing assignment to final_state_particles is pure convenience, but it is Pythonic. We need this propery anyway to have an interface to get the list of particles. It is then natural/pythonic to also support assignment to this property. We don't need another method to add/remove individual particles, we have set_stable for that.

[jncots]

Pythia8DecayHandler should not be called by MCRun. It is not a generic component that all models need.

I checked all models what they decay and what not. It turns out that only EposLHC decays everything that could be passed to it as unstable particles. Models that use Pythia6 (Sibyll, Sophia) and Pythia6 itself cannot decay neutron. Phojet, Dpmjet - just not decay what set as unstable. UrQMD and QGSJet just don't decay particles.

So, I would say we need Pythia8DecayHandler for pretty much everything.

I removed the possibility for switching on/off Pythia8DecayHandler. Pythia8DecayHandler will not be used where Pythia8 is not available (like on Windows). In other cases, the event.pid is checked, and if there are particles which are not in final_state_particles Pythia8DecayHandler will be used.

In all models _decaying_pids first set to unstable, and then final_state_particles are set for stable. In such case we are sure that there are no unstable particles, except that we are set. _decaying_pids are model-specific because different models know different types of particles.

Even by setting stable particles as above, most of the models ignores some particles to decay:

• Phojet and Dpmjet, UrQMD ignores many different types of particles.
• QGSJet is known to not decay any
• Sibyll, Sophia, Pythia6 do not decay neutron.
• Sibyll21 freezes if [-13, 13, 130, 310] are set as decaying

Only EposLHC decays all that is set as decaying.

Thus, currently Pythia8DecayHandler will be called only if the model produces particles which are not in final_state_particles. For a standard long_lived particles (30 psec) Pythia8DecayHandler will be used by QGSJet, Phojet, Dpmjet, and UrQMD. It probably could be used sometimes by Sibyll because I noticed it happens from time to time that stack contains undecayed particles. Sibyll21 will use Pythia8DecayHandler when one of [-13, 13, 130, 310] is unstable.

[jncots]

I notice one unpleasant side effect of using Pythia8 as DecayHandler: it looks like one can set stable/decay only for the pair of particle-antiparticle. I.e. one cannot set pi+ as stable and pi- as unstable.

[jncots]

_validate_decay function is added which checks whether all particles that should decay have been decayed. If not, it uses decay_handler if it is switched on, or suggest to do that if decay_handler is switched off.

[afedynitch]

@jncots, could you pull from main before preparing for merge? I forgot that SIBYLL* has internal decays and it works differently compared to normal SIBYLL. We won't get short lived particles out of it.

At very least it must XFAIL on test_final_state.py.

[jncots]

@jncots, could you pull from main before preparing for merge?

Because the next merge is for #170, I've done it there.

[jncots]

After @afedynitch 's review:

• Renamed long_lived_for to select_long_lived
• Made lists and function related to unstable particles to be private
• Changed decaying to unstable where possible
• Change an order in _set_final_state_particles
• Add comments to run_decay_handler