3.5.  Hadronic Physics

3.5.1.  Introduction

Optimal exploitation of hadronic final states played a key role in successes of all recent collider experiment in HEP, and the ability to use hadronic final states will continue to be one of the decisive issues during the analysis phase of the LHC experiments. Monte Carlo programs like Geant4 facilitate the use of hadronic final states, and have been developed for many years.

We give an overview of the Object Oriented frameworks for hadronic generators in Geant4, and illustrate the physics models underlying hadronic shower simulation on examples, including the three basic types of modeling; data-driven, parametrisation-driven, and theory-driven modeling, and their possible realisations in the Object Oriented component system of Geant4. We put particular focus on the level of extendibility that can and has been achieved by our Russian dolls approach to Object Oriented design, and the role and importance of the frameworks in a component system.

3.5.2.  Principal Considerations

The purpose of this section is to explain the implementation frameworks used in and provided by Geant4 for hadronic shower simulation as in the 1.1 release of the program. The implementation frameworks follow the Russian dolls approach to implementation framework design. A top-level, very abstracting implementation framework provides the basic interface to the other Geant4 categories, and fulfils the most general use-case for hadronic shower simulation. It is refined for more specific use-cases by implementing a hierarchy of implementation frameworks, each level implementing the common logic of a particular use-case package in a concrete implementation of the interface specification of one framework level above, this way refining the granularity of abstraction and delegation. This defines the Russian dolls architectural pattern. Abstract classes are used as the delegation mechanism [1] .

All framework functional requirements were obtained through use-case analysis. In the following we present for each framework level the compressed use-cases, requirements, designs including the flexibility provided, and illustrate the framework functionality with examples. All design patterns cited are to be read as defined in [ Gamma1995 ].

3.5.3.  Level 1 Framework - processes

There are two principal use-cases of the level 1 framework. A user will want to choose the processes used for his particular simulation run, and a physicist will want to write code for processes of his own and use these together with the rest of the system in a seamless manner.

Requirements

  1. Provide a standard interface to be used by process implementations.
  2. Provide registration mechanisms for processes.

Design and interfaces

Both requirements are implemented in a sub-set of the tracking-physics interface in Geant4}. The class diagram is shown in Figure 3.3.

Level 1 implementation framework of the hadronic category of GEANT4.

Figure 3.3.  Level 1 implementation framework of the hadronic category of GEANT4.


All processes have a common base-class G4VProcess, from which a set of specialised classes are derived. Two of them are used as base classes for hadronic processes at rest and in flight (G4VDiscreteProcess), and for processes like radioactive decay where the same implementation can represent both these extreme cases (G4VRestDiscreteProcess).

Each of these classes declares two types of methods; one for calculating the time to interaction or the physical interaction length, allowing tracking to request the information necessary to decide on the process responsible for final state production, and one to compute the final state. These are pure virtual methods, and have to be implemented in each individual derived class, as enforced by the compiler.

Note on at-rest processes: starting with Geant4 version 9.6 - when the Bertini and Fritiof final-state models have been extended down to zero kinetic energy and used also for simulating the nuclear capture at-rest - the at-rest processes derive from G4HadronicProcess, hence from G4VDiscreteProcess, instead than from G4VRestProcess as in the initial design of at-rest processes. This requires some adaptation a discrete process to handle an at-rest one using top level interface G4VProcess. A different solution, under consideration but not yet implemented, would be instead to have G4HadronicProcess inheriting from G4VRestDiscreteProcess: in this way, G4HadronicProcess, and therefore any theory-driven final-state model, could be deployed for any kind of hadronic process, including capture-at-rest processes and radioactive decays.

Framework functionality

The functionality provided is through the use of process base-class pointers in the tracking-physics interface, and the G4Process-Manager. All functionality is implemented in abstract, and registration of derived process classes with the G4Process-Manager of an individual particle allows for arbitrary combination of both Geant4 provided processes, and user-implemented processes. This registration mechanism is a modification on a Chain of Responsibility. It is outside the scope of the current paper, and its description is available from G4Manual.

3.5.4.  Level 2 Framework - Cross Sections and Models

At the next level of abstraction, only processes that occur for particles in flight are considered. For these, it is easily observed that the sources of cross sections and final state production are rarely the same. Also, different sources will come with different restrictions. The principal use-cases of the framework are addressing these commonalities. A user might want to combine different cross sections and final state or isotope production models as provided by Geant4, and a physicist might want to implement his own model for particular situation, and add cross-section data sets that are relevant for his particular analysis to the system in a seamless manner.

Requirements

  1. Flexible choice of inclusive scattering cross-sections.
  2. Ability to use different data-sets for different parts of the simulation, depending on the conditions at the point of interaction.
  3. Ability to add user-defined data-sets in a seamless manner.
  4. Flexible, unconstrained choice of final state production models.
  5. Ability to use different final state production codes for different parts of the simulation, depending on the conditions at the point of interaction.
  6. Ability to add user-defined final state production models in a seamless manner.
  7. Flexible choice of isotope production models, to run in parasitic mode to any kind of transport models.
  8. Ability to use different isotope production codes for different parts of the simulation, depending on the conditions at the point of interaction.
  9. Ability to add user-defined isotope production models in a seamless manner.

Design and interfaces

The above requirements are implemented in three framework components, one for cross-sections, final state production, and isotope production each. The class diagrams are shown in Figure 3.4 for the cross-section aspects, Figure 3.5 for the final state production aspects, and figure Figure 3.6 for the isotope production aspects.

Level 2 implementation framework of the hadronic category of Geant4; cross-section aspect.

Figure 3.4.  Level 2 implementation framework of the hadronic category of Geant4; cross-section aspect.


Level 2 implementation framework of the hadronic category of Geant4; final state production aspect.

Figure 3.5.  Level 2 implementation framework of the hadronic category of Geant4; final state production aspect.


Level 2 implementation framework of the hadronic category of Geant4; isotope production aspect

Figure 3.6.  Level 2 implementation framework of the hadronic category of Geant4; isotope production aspect


The three parts are integrated in the G4Hadronic-Process class, that serves as base-class for all hadronic processes of particles in flight.

Cross-sections

Each hadronic process is derived from G4Hadronic-Process}, which holds a list of "cross section data sets". The term "data set" is representing an object that encapsulates methods and data for calculating total cross sections for a given process in a certain range of validity. The implementations may take any form. It can be a simple equation as well as sophisticated parameterisations, or evaluated data. All cross section data set classes are derived from the abstract class G4VCrossSection-DataSet}, which declares methods that allow the process inquire, about the applicability of an individual data-set through IsApplicable(const G4DynamicParticle*, const G4Element*), and to delegate the calculation of the actual cross-section value through GetCrossSection(const G4DynamicParticle*, const G4Element*).

Final state production

The G4HadronicInteraction base class is provided for final state generation. It declares a minimal interface of only one pure virtual method: G4VParticleChange* ApplyYourself(const G4Track &, G4Nucleus &)}. G4HadronicProcess provides a registry for final state production models inheriting from G4Hadronic-Interaction. Again, final state production model is meant in very general terms. This can be an implementation of a quark gluon string model [QGSM], a sampling code for ENDF/B data formats [ ENDFForm ], or a parametrisation describing only neutron elastic scattering off Silicon up to 300~MeV.

Isotope production

For isotope production, a base class (G4VIsotope-Production) is provided. It declares a method G4IsoResult * GetIsotope(const G4Track &, const G4Nucleus &) that calculates and returns the isotope production information. Any concrete isotope production model will inherit from this class, and implement the method. Again, the modeling possibilities are not limited, and the implementation of concrete production models is not restricted in any way. By convention, the GetIsotope method returns NULL, if the model is not applicable for the current projectile target combination.

Framework functionality:

Cross Sections

G4HadronicProcess provides registering possibilities for data sets. A default is provided covering all possible conditions to some approximation. The process stores and retrieves the data sets through a data store that acts like a FILO stack (a Chain of Responsibility with a First In Last Out decision strategy). This allows the user to map out the entire parameter space by overlaying cross section data sets to optimise the overall result. Examples are the cross sections for low energy neutron transport. If these are registered last by the user, they will be used whenever low energy neutrons are encountered. In all other conditions the system falls back on the default, or data sets with earlier registration dates. The fact that the registration is done through abstract base classes with no side-effects allows the user to implement and use his own cross sections. Examples are special reaction cross sections of K0-nuclear interactions that might be used for ε/ε ' analysis at LHC to control the systematic error.

Final state production

The G4HadronicProcess class provides a registration service for classes deriving from G4Hadronic-Interaction, and delegates final state production to the applicable model. G4Hadronic-Interactionprovides the functionality needed to define and enforce the applicability of a particular model. Models inheriting from G4Hadronic-Interaction can be restricted in applicability in projectile type and energy, and can be activated/deactivated for individual materials and elements. This allows a user to use final state production models in arbitrary combinations, and to write his own models for final state production. The design is a variant of a Chain of Responsibility. An example would be the likely CMS scenario - the combination of low energy neutron transport with a quantum molecular dynamics [QMD], invariant phase space decay [CHIPS], and fast parametrised models for calorimeter materials, with user defined modeling of interactions of spallation nucleons with the most abundant tracker and calorimeter materials.

Isotope production

The G4HadronicProcess by default calculates the isotope production information from the final state given by the transport model. In addition, it provides a registering mechanism for isotope production models that run in parasitic mode to the transport models and inherit from G4VIsotope-Production. The registering mechanism behaves like a FILO stack, again based on Chain of Responsibility. The models will be asked for isotope production information in inverse order of registration. The first model that returns a non-NULL value will be applied. In addition, the G4Hadronic-Process provides the basic infrastructure for accessing and steering of isotope-production information. It allows to enable and disable the calculation of isotope production information globally, or for individual processes, and to retrieve the isotope production information through the G4IsoParticleChange * GetIsotopeProductionInfo()} method at the end of each step. The G4HadronicProcess is a finite state machine that will ensure the GetIsotope-ProductionInfo returns a non-zero value only at the first call after isotope production occurred. An example of the use of this functionality is the study of activation of a Germanium detector in a high precision, low background experiment.

3.5.5.  Level 3 Framework - Theoretical Models

Level 3 implementation framework of the hadronic category of Geant4; theoretical model aspect.

Figure 3.7.  Level 3 implementation framework of the hadronic category of Geant4; theoretical model aspect.


Geant4 provides at present one implementation framework for theory driven models. The main use-case is that of a user wishing to use theoretical models in general, and to use various intra-nuclear transport or pre-compound models together with models simulating the initial interactions at very high energies.

Requirements

  1. Allow to use or adapt any string-parton or parton transport [VNI],
  2. Allow to adapt event generators, ex. [PYTHIA7], state production in shower simulation.
  3. Allow for combination of the above with any intra-nuclear transport (INT).
  4. Allow stand-alone use of intra-nuclear transport.
  5. Allow for combination of the above with any pre-compound model.
  6. Allow stand-alone use of any pre-compound model.
  7. Allow for use of any evaporation code.
  8. Allow for seamless integration of user defined components for any of the above.

Design and interfaces

To provide the above flexibility, the following abstract base classes have been implemented:

  • G4VHighEnergyGenerator
  • G4VIntranuclearTransportModel
  • G4VPreCompoundModel
  • G4VExcitationHandler

In addition, the class G4TheoFS-Generator is provided to orchestrate interactions between these classes. The class diagram is shown in Figure 3.7.

G4VHighEnergy-Generator serves as base class for parton transport or parton string models, and for Adapters to event generators. This class declares two methods, Scatter, and GetWoundedNucleus.

The base class for INT inherits from G4Hadronic-Interaction, making any concrete implementation usable as a stand-alone model. In doing so, it re-declares the ApplyYourself interface of G4Hadronic-Interaction, and adds a second interface, Propagate, for further propagation after high energy interactions. Propagate takes as arguments a three-dimensional model of a wounded nucleus, and a set of secondaries with energies and positions.

The base class for pre-equilibrium decay models, G4VPre-CompoundModel, inherits from G4Hadronic-Interaction, again making any concrete implementation usable as stand-alone model. It adds a pure virtual DeExcite method for further evolution of the system when intra-nuclear transport assumptions break down. DeExcite takes a G4Fragment, the Geant4 representation of an excited nucleus, as argument.

The base class for evaporation phases, G4VExcitation-Handler, declares an abstract method, BreakItUP(), for compound decay.

Framework functionality

The G4TheoFSGenerator class inherits from G4Hadronic-Interaction, and hence can be registered as a model for final state production with a hadronic process. It allows a concrete implementation of G4VIntranuclear-TransportModel and G4VHighEnergy-Generator to be registered, and delegates initial interactions, and intra-nuclear transport of the corresponding secondaries to the respective classes. The design is a complex variant of a Strategy. The most spectacular application of this pattern is the use of parton-string models for string excitation, quark molecular dynamics for correlated string decay, and quantum molecular dynamics for transport, a combination which promises to result in a coherent description of quark gluon plasma in high energy nucleus-nucleus interactions.

The class G4VIntranuclearTransportModel provides registering mechanisms for concrete implementations of G4VPreCompound-Model, and provides concrete intra-nuclear transports with the possibility of delegating pre-compound decay to these models.

G4VPreCompoundModel provides a registering mechanism for compound decay through the G4VExcitation-Handler interface, and provides concrete implementations with the possibility of delegating the decay of a compound nucleus.

The concrete scenario of G4TheoFS-Generator using a dual parton model and a classical cascade, which in turn uses an exciton pre-compound model that delegates to an evaporation phase, would be the following: G4TheoFS-Generator receives the conditions of the interaction; a primary particle and a nucleus. It asks the dual parton model to perform the initial scatterings, and return the final state, along with the by then damaged nucleus. The nucleus records all information on the damage sustained. G4TheoFS-Generator forwards all information to the classical cascade, that propagates the particles in the already damaged nucleus, keeping track of interactions, further damage to the nucleus, etc.. Once the cascade assumptions break down, the cascade will collect the information of the current state of the hadronic system, like excitation energy and number of excited particles, and interpret it as a pre-compound system. It delegates the decay of this to the exciton model. The exciton model will take the information provided, and calculate transitions and emissions, until the number of excitons in the system equals the mean number of excitons expected in equilibrium for the current excitation energy. Then it hands over to the evaporation phase. The evaporation phase decays the residual nucleus, and the Chain of Command rolls back to G4TheoFS-Generator, accumulating the information produced in the various levels of delegation.

3.5.6.  Level 4 Frameworks - String Parton Models and Intra-Nuclear Cascade

Level 4 implementation framework of the hadronic category of Geant4; parton string aspect.

Figure 3.8.  Level 4 implementation framework of the hadronic category of Geant4; parton string aspect.


Level 4 implementation framework of the hadronic category of Geant4; intra-nuclear transport aspect.

Figure 3.9.  Level 4 implementation framework of the hadronic category of Geant4; intra-nuclear transport aspect.


The use-cases of this level are related to commonalities and detailed choices in string-parton models and cascade models. They are use-cases of an expert user wishing to alter details of a model, or a theoretical physicist, wishing to study details of a particular model.

Requirements

  1. Allow to select string decay algorithm
  2. Allow to select string excitation.
  3. Allow the selection of concrete implementations of three-dimensional models of the nucleus
  4. Allow the selection of concrete implementations of final state and cross sections in intra-nuclear scattering.

Design and interfaces

To fulfil the requirements on string models, two abstract classes are provided, the G4VParton-StringModel and the G4VString-Fragmentation. The base class for parton string models, G4VParton-StringModel, declares the GetStrings() pure virtual method. G4VString-Fragmentation, the pure abstract base class for string fragmentation, declares the interface for string fragmentation.

To fulfill the requirements on intra-nuclear transport, two abstract classes are provided, G4V3DNucleus, and G4VScatterer. At this point in time, the usage of these intra-nuclear transport related classes by concrete codes is not enforced by designs, as the details of the cascade loop are still model dependent, and more experience has to be gathered to achieve standardisation. It is within the responsibility of the implementers of concrete intra-nuclear transport codes to use the abstract interfaces as defined in these classes.

The class diagram is shown in Figure 3.8 for the string parton model aspects, and in Figure 3.9 for the intra-nuclear transport.

Framework functionality

Again variants of Strategy, Bridge and Chain of Responsibility are used. G4VParton-StringModel implements the initialisation of a three-dimensional model of a nucleus, and the logic of scattering. It delegates secondary production to string fragmentation through a G4VString-Fragmentation pointer. It provides a registering service for the concrete string fragmentation, and delegates the string excitation to derived classes. Selection of string excitation is through selection of derived class. Selection of string fragmentation is through registration.

3.5.7.  Level 5 Framework - String De-excitation}

Level 5 implementation framework of the hadronic category of Geant4; string fragmentation aspect.

Figure 3.10.  Level 5 implementation framework of the hadronic category of Geant4; string fragmentation aspect.


The use-case of this level is that of a user or theoretical physicist wishing to understand the systematic effects involved in combining various fragmentation functions with individual types of string fragmentation. Note that this framework level is meeting the current state of the art, making extensions and changes of interfaces in subsequent releases likely.

Requirements

  1. Allow the selection of fragmentation function.

Design and interfaces

A base class for fragmentation functions, G4VFragmentation-Function}, is provided. It declares the GetLightConeZ() interface.

Framework functionality

The design is a basic Strategy. The class diagram is shown in Figure 3.10. At this point in time, the usage of the G4VFragmentation-Function is not enforced by design, but made available from the G4VString-Fragmentation to an implementer of a concrete string decay. G4VString-Fragmentation provides a registering mechanism for the concrete fragmentation function. It delegates the calculation of zf of the hadron to split of the string to the concrete implementation. Standardisation in this area is expected.

3.5.8. Creating Your Own Hadronic Process

For some applications Geant4 might not provide the most appropriate physics implementantion or, in fact, any physics implementation at all. In such cases, it is up to the user to develop the necessary processes, models and cross sections and integrate them into his version of the Geant4 toolkit. The user's process, model or cross section may then be used in a physics list to replace some of those already provided by the toolkit. This modularity requires that user classes be derived from a set of base classes which have been provided to aid integration with the toolkit and to spare the user from consideration of many details not related to the physics at hand.

Processes communicate with Geant4 tracking, telling it where or when an interaction is supposed to occur, and what is supposed to happen at that point. A hadronic process may be implemented directly, or through the use af a framework of classes that modularize physics functionality, make available several utilities and reduce unneccesary code duplication. In the latter, recommended, approach the user must in general develop as many as three classes: a process, a cross section and a model. Instances of the cross section and model classes must then be assigned to the process. In practice it is usually necessary to develop only a model class or a cross section class, since a number of processes are already provided by Geant4. Before writing any code, users should check that Geant4 has not already povided the necessary models, cross sections or processes.

3.5.8.1. Developing a new hadronic model

A hadronic model is repsonsible for the generation of a set of final state four-vectors, given an initial projectile and target. New models should derive from the G4HadronicInteraction base class and at least two methods in this class must be implemented:

virtual G4HadFinalState*
ApplyYourself(const G4HadProjectile& aProjectile, G4Nucleus& targetNucleus)

which is responsible for generating the final state of the interaction including the specification of all particle types and four-momenta, and

virtual G4bool IsApplicable(const G4HadProjectile& aProjectile, G4Nucleus& targetNucleus)

which is responsible for checking that the incident particle type and energy, and the Z and A of the target nucleus, can be adequately handled by the model. G4HadronicInteraction provides a number of utilities to aid in the implementation of these methods.

When implementing ApplyYourself(), the Get() methods of the G4HadProjectile and G4Nucleus classes provide all the initial state information necessary for the generation of the final state. For G4HadProjectile, GetDefinition() provides the particle type, and Get4Momentum() provides the total energy and momentum. For G4Nucleus, GetZ_asInt(), GetN_asInt() and GetA_asInt() provide Z, N and A, while AtomicMass() provides the mass. Additional utility methods are available for both G4HadProjectile and G4Nucleus .

Coordinate systems

The inputs to the model assume that the incident particle (G4HadProjectile) travels along the z axis and interacts with the target (G4Nucleus) which is at rest in the lab frame. Before invoking the ApplyYourself() method, the process rotates the direction of the projectile to be along the z axis and then performs the inverse rotation on the final state particles produced by the model.

The model must perform two additional transformations: into the CM frame, and back out of it after the interaction is complete.

Writing the ApplyYourself() method

It is thus the model developer's responsibility to:

  • boost the projectile and target into the CM frame using the necessary Lorentz transformations,
  • perform all calculations required to generate the final state set of particles,
  • boost the final state particles back to the lab frame with the inverse transformation, and
  • send the final state particles to the process by filling G4HadFinalState and setting its status.

Step 4) is accomplished by using the various Get() and Set() methods provided by the class G4HadFinalState. The developer must also decide whether the original projectile survives the interaction, disappears or is suspended. This is done with the SetStatusChange() method. If the projectile survives, the change in its energy and momentum must be set with the provided methods and it must be flagged as "isAlive". If the particle disappears it must be flagged as "stopAndKill".

Geant4 provides a large number of Lorentz transformation tools which may be used to complete steps 1) and 3).

How step 2) is accomplished is entirely up to the developer. This could be as simple as a look-up table which assigns a final state to an initial state, or as complex as a theoretical high energy generator. Typically, the user will have to provide methods of sampling final state multiplicities, energies and angles using random number generators provided by Geant4. For example, the cosine of the polar angle of an isotropic angular distribution could be sampled as follows:

  G4double cosTheta = 2.*G4UniformRandom() - 1.;

The developer must also see to it that the model conserves energy and momentum. Currently the hadronic framework checks that final states do not exceed reasonably small limits of non-conservation.

Using the hadronic framework

For complex models it is recommended that the user become familiar with the Geant4 hadronic framework. This is covered in detail in the chapter on Extended Functionality in this manual. The framework uses the object-oriented principles of abstraction and re-use to provide a number of services to the developer. The part of the framework used will depend on the type of model. For example, high energy models can take advantage of already-developed string excitation and decay functions and medium energy models can use the intra-nuclear propagation base class and the nuclear de-excitation handler.

Writing the IsApplicable() method

The is a straightforward, but important, method. Most models are quite specific in their range of use and the developer must codify this. It is recommended that this method test for ranges of projectile energy, particle type and target atomic number and weight, and return false when these ranges are exceeded.

3.5.8.2. Developing a new cross section set

New cross section sets should derive from G4VCrossSectionDataSet . This class serves as a container of cross section data and provides a number of access methods that must be implemented by the developer. The essential methods are:

G4double GetElementCrossSection(const G4DynamicParticle*, G4int Z)

which retrieves element-based cross sections,

G4double GetIsoCrossSection(const G4DynamicParticle*, G4int Z, G4int A)

which retrieves isotope-based cross sections,

G4bool IsElementApplicable(const G4DynamicParticle*, G4int Z)

which sets the Z range of the element-based data set,

G4bool IsIsoApplicable(const G4DynamicParticle*, G4int, G4int A)

which sets the Z and A range of the isotope-based data set, and

SetMinKinEnergy(G4double)
SetMaxKinEnergy(G4double)
GetMinKinEnergy()
and GetMaxKinEnergy()

for defining the applicable energy range of the data set.

3.5.8.3. Developing a new hadronic process

As mentioned above, it is preferable to add new physics in terms of a model, and assign the model to an existing process, rather than develop a new, specific process. Under certain circumstances though, a directly implemented process may be necessary. In that case it must derive from G4HadronicProcess and three methods of that class must be implemented:

virtual G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&) ,
virtual G4bool IsApplicable(const G4ParticleDefinition&) , and
G4double GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*).

PostStepDoIt() is responsible for generating the final state of an interaction given the track and step information. It must update the state of the track, flagging it as "Alive", "StopButAlive", "StopAndKill", "KillTrackAndSecondaries", "Suspend", or "PostponeToNextEvent". It is roughly analagous to the ApplyYourself() method in models.

IsApplicable() serves the same purpose in processes as it does in models.

GetMeanFreePath() gets the cross section as a function of particle type, energy, and target material, and converts it to a mean free path, which is in turn passed on to the tracking. This method can be quite simple:

G4double particle = aTrack.GetDynamicParticle();
G4double material = aTrack.GetMaterial();
return factor/theCrossSectionDataStore->GetCrossSection(particle, material);

provided an appropriate cross section data set is already available in the data store.

The above discussion refers to in-flight processes. At-rest hadronic processes do not currently derive from G4HadronicProcess, but from G4VRestProcess or G4VRestDiscreteProcess. As such they do not employ the full hadronic framework and must be implemented directly without models. The methods to be implemented are similar to those in the in-flight case:

virtual G4VParticleChange* AtRestDoIt(const G4Track&, const G4Step&) , and
G4bool IsApplicable(const G4ParticleDefinition&) .

3.5.9. Changing parameters of an existing model: Fritiof (FTF) use-case

3.5.9.1. Total, elastic, inelastic and annihilation cross sections

Parameters of the FTF model are concentrated in G4FTFParameters.cc located in geant4/source/processes/hadronic/models/parton_string/diffraction/src

Input parameters of the code are: type of the projectile particle, mass number and charge of the target nucleus, and the projectile momentum for elementary projectiles, or momentum per nucleon in the case of a projectile nucleus (MeV/c).

A simplified Glauber model is used in the code for calculation of the multiplicity of intra-nuclear collisions. Thus, it is needed to estimate the average total, elastic and inelastic cross sections for interactions of a projectile with target protons and neutrons. In the case of nucleus-nucleus interactions it is needed to know average NN cross sections (averaged over PP, PN, and NN collisions). According to the Glauber approximation, changes of the cross sections during a collision are not considered.

Chips cross sections of elementary interactions is used in the code. A pointer to the cross sections (FTFxsManager) is determined in the lines 178 — 183.

The average cross sections of a projectile interactions with protons and neutrons are determined in the lines 192 — 475. A special calculations for projectile anti-baryons are presented in the lines 236 — 357.

According to a developed approach for antibaryon-nucleon interactions the list of processes shown in Figure 3.11 is considered.

Quark flow diagrams of $\bar pp$-interactions.

Figure 3.11.  Quark flow diagrams of $\bar pp$-interactions.


As usual, quarks and anti-quarks are shown by solid lines. Dashed lines present so-called string junctions. It is assumed that the gluon field in baryons has a non-trivial topology. This heterogeneity is called a "string junction". Quark-gluon strings produced in the reaction are shown by wavy lines.

The diagram of Figure 3.11a represents a process with a string junction annihilation and the creation of three strings. Diagram Figure 3.11b describes quark-antiquark annihilation and string creation between the di-quark and anti-di-quark. Quark-anti-quark and string junction annihilation is shown in Figure 3.11c. Finally, one string is created in the process of Figure 3.11e. Hadrons appear at the fragmentation of the strings in the same way that they appear in $e^+e^-$-annihilation. One can assume that excited strings with complicated gluonic field configurations are created in processes Figure 3.11d and Figure 3.11f. If the collision energy is sufficiently small, glueballs can be formed in the process Figure 3.11f. Mesons with constituent gluons or with hidden baryon number can be created in process Figure 3.11d. Of course the standard FTF processes shown in the bottom of the figure are also allowed. In the simplest approach it is assumed that the energy dependence of the cross sections of these processes vary inversely with a power of s as depicted in Figure 3.11. Here s is the center-of-mass energy squared. This is dictated by the reggeon phenomenology. Calculating the cross sections of binary reactions is a rather complicated procedure (see [ Kaidalov_Bin ]) because there can be interactions in the initial and final states. These interactions reflect also on cross sections of other reactions [ Uzhinsky_GaloyanPbarP ]. The annihilation cross section is given as:
where Xi are yields of the diagrams of Figure 3.11. All cross sections are given in [mb].
Table of the coefficients B, C and D:
¯pp ¯pn ¯np ¯nn ¯Λp ¯Λn ¯Σ-p ¯Σ-n ¯Σ0p ¯Σ0n ¯Σ+p ¯Σ+n
B 5 4 4 5 3 3 2 4 3 3 4 2
C 5 4 4 5 3 3 2 4 3 3 4 2
D 6 4 4 6 3 3 2 2 2 2 2 0
¯Ξ-p ¯Ξ-n ¯Ξ0p ¯Ξ0n ¯Ω-p ¯Ω-n
B 1 2 2 1 0 0
C 1 2 2 1 0 0
D 0 0 0 0 0 0
The coefficients B, C and D are pure combinatorial coefficients calculated on the assumption that the same conditions apply to all quarks and anti-quarks. For example, in ¯pp interactions there are five possibilities to annihilate a quark and an anti-quark, and six possibilities to annihilate two quarks and two anti-quarks. Thus, B=C=5 and D=6. Note that final state particles in the process of Figure 3.11b can coincide with initial state particles. Thus the true elastic cross section is not given by the experimental cross section. At Plab < 40 MeV/c anti-proton-nucleon cross sections are:
All cross sections are given in [mb]. At Plab < 40 MeV/c σb must be equal 0 for ¯pp-interactions because the process ¯pp¯nn is impossible at these energies. It is not so for ¯nn interactions because the process ¯nn¯pp can take place. In the case ¯np interactions, the process "b" will lead to elastic scattering which is not considered. The parameterizations are presented in the lines 236 — 357. Their changes influence on ¯pp channel cross sections. Thus, the changes must be done with caution. To complete the cross section calculations, σFTF responsible for processes presented in bottom row of Figure 3.11 is determined (line 346) as:
where sth is the threshold energy of meson production 2 [MeV]. After that, after line 478, needed quantities are stored in G4FTFParameters class:
  SetTotalCrossSection( Xtotal );
  SetElastisCrossSection( Xelastic );
  SetInelasticCrossSection( Xtotal - Xelastic );

  SetProbabilityOfElasticScatt( Xtotal, Xelastic );
  SetRadiusOfHNinteractions2( Xtotal/pi/10.0 );

  SetProbabilityOfAnnihilation( Xannihilation / (Xtotal - Xelastic) );

  SetSlope( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 );
  SetGamma0( GetSlope()*Xtotal/10.0/2.0/pi );

  // Parameters of elastic scattering
  // Gaussian parametrization of elastic scattering amplitude assumed
  SetAvaragePt2ofElasticScattering( 1.0/( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 )

  G4double Xinel = Xtotal - Xelastic;

3.5.9.2. Cross sections of elementary processes

The Fritiof model[ Fritiof1 ][ Fritiof2 ] assumes that all hadron-hadron interactions are binary reactions,

h1+h2h1'+h2' where h1' and h2' are excited states of the hadrons with discrete or continuous mass spectra (see Figure 3.12). If one of the final hadrons is in its ground state (h1+h2h1+h2' ) the reaction is called "single diffraction dissociation", and if neither hadron is in its ground state it is called a "non-diffractive" interaction.

Non-diffractive and diffractive interactions considered in the Fritiof model.

Figure 3.12.  Non-diffractive and diffractive interactions considered in the Fritiof model.


The excited hadrons are considered as QCD-strings, and the corresponding LUND-string fragmentation model is applied in order to simulate their decays. The key ingredient of the Fritiof model is the sampling of the string masses. In general, the set of final states for the interactions can be represented by Figure 3.13, where samples of possible string masses are shown. There is a point corresponding to elastic scattering, a group of points which represents final states of binary hadron-hadron interactions, lines corresponding to the diffractive interactions, and various intermediate regions. The region populated with the red points is responsible for the non-diffractive interactions. In the model, the mass sampling threshold is set equal to the ground state hadron masses, but in principle the threshold can be lower than these masses. The string masses are sampled in the triangular region restricted by the diagonal line corresponding to the kinematical limit M1+M2 = Ecms where M1 and M2 are the masses of the h1' and h2' hadrons, and also of the threshold lines. If a point is below the string mass threshold, it is shifted to the nearest diffraction line.
Diagram of the final states of hadron-hadron interactions.

Figure 3.13.  Diagram of the final states of hadron-hadron interactions.


Unlike the original Fritiof model, the final state diagram of the current model is complicated, which leads to a non-straightforward mass sampling algorithm. This will be considered below. The original model had no points corresponding to elastic scattering or to the binary final states. As it was known at the time, the mass of an object produced by diffraction dissociation, Mx, for example from the reaction p+p → p + X, is distributed as dMx∕ Mx ∝ dMx2∕ Mx2, so it was natural to assume that the object mass distributions in all inelastic interactions obey the same law. This can be re-written using the light-cone momentum variables, P+ or P-,
where E is the energy of a particle, and pz is its longitudinal momentum along the collision axis. At large energy and positive pz, P-≈ M2+PT2 ∕ 2pz, P+≈ M2+PT2 ∕ 2|pz|. Usually, the transferred transverse momentum, PT, is small and can be neglected. Thus, it was assumed that P- and P+ of a projectile and target associated hadron, respectively, are distributed as
A gaussian distribution was used to sample PT. In the case of hadron-nucleus or nucleus-nucleus interactions it was assumed that the created objects can interact further with other nuclear nucleons and create new objects. Assuming equal masses of the objects, the multiplicity of particles produced in these interactions will be proportional to the number of participating nuclear nucleons, or to the multiplicity of intra-nuclear collisions. Due to this, the multiplicity of particles produced in hadron-nucleus or nucleus-nucleus interactions is larger than that in hadron-hadron ones. The probabilities of multiple intra-nuclear collisions were sampled with the help of a simplified Glauber model. Cascading of secondary particles was not considered. Because the Fermi motion of nuclear nucleons was simulated in a simple manner, the original Fritiof model could not work at Plab< 10—20 GeV/c. It was assumed in the model that the created objects are quark-gluon strings with constituent quarks at their ends originating from the primary colliding hadrons. Thus, the LUND-string fragmentation model was applied for a simulation of the object decays. It was assumed also that the strings with sufficiently large masses have "kinks" — additional radiated gluons. This is very important for a correct reproduction of particle multiplicities in the interactions. All of the above assumptions were reconsidered in the implementation of the Geant4 Fritiof model, and new features were added. These will be presented below.
Quark flow diagrams of π N interactions.

Figure 3.14.  Quark flow diagrams of π N interactions.


The original Fritiof model contains only the pomeron exchange process shown in Figure 3.14d . It would be useful to extend the model by adding the exchange processes shown in Figure 3.11b and Figure 3.11c, and the annihilation process of Figure 3.11a . This could probably be done by introducing a restricted set of mesonic and baryonic resonances and a corresponding set of parameters. This procedure was employed in the binary cascade model of Geant4 (BIC) [ BIC ] and in the Ultra-Relativistic-Quantum-Molecular-Dynamic model (UrQMD) [ UrQMD1 ][ UrQMD2 ]. However, it is complicated to use this solution for the simulation of hadron-nucleus and nucleus-nucleus interactions. The problem is that one has to consider resonance propagation in the nuclear medium and take into account their possible decays which enormously increases computing time. Thus, in the current version of the FTF model only quark exchange processes have been added to account for meson and baryon interactions with nucleons, without considering resonance propagation and decay. This is a reasonable hypothesis at sufficiently high energies. For each projectile hadrons the following probabilities are set up:

  • Probability of quark exchange process without excitation of participants (Figure 3.14b);
  • Probability of quark exchange process with excitation of participants (Figure 3.14c);
  • Probability of projectile diffraction dissociation;
  • Probability of target diffraction dissociation.

All these probabilities have the same functional form:
where y is the projectile rapidity in the target rest frame. For each projectile hadron and for each possible process, SetParams method is called for setting up the parameters A1, B1, A2, B2, A3. For projectile baryons, this is implemented in lines 524 — 542.
  //         Proc#  A1        B1            A2       B2    A3  Atop  Ymin
  SetParams(   0, 13.71,     1.75,          -214.5, 4.25, 0.0, 0.5 , 1.1 );
  SetParams(   1, 25.0 ,     1.0 ,          -50.34, 1.5 , 0.0, 0.0 , 1.4 );
  if ( Xinel > 0.0 ) {
    SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 );
    SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 );
    SetParams( 4, 1.0 ,      0.0 ,           -2.01, 0.5 , 0.0, 0.0 , 1.4 );
  } else {
    SetParams( 2, 0.0,       0.0 ,            0.0 , 0.0 , 0.0, 0.0 , 0.0 );
    SetParams( 3, 0.0,       0.0 ,            0.0 , 0.0 , 0.0, 0.0 , 0.0 );
    SetParams( 4, 0.0,       0.0 ,            0.0 , 0.0 , 0.0, 0.0 , 0.0 );
  }
  if ( AbsProjectileBaryonNumber > 1  ||  NumberOfTargetNucleons > 1 ) {
    SetParams( 2, 0.0,       0.0 ,            0.0 , 0.0 , 0.0, 0.0 , -100.0 );
    //SetParams( 3, 0.0,       0.0 ,            0.0 , 0.0 , 0.0, 0.0 , -100.0 );
  }
Proc# are the quark exchange without participant excitation (0), the quark exchange with participant excitation (1), projectile diffraction (2), target diffraction (3). The probability of the quark exchange with participant excitation has an additional multiplier (Proc#=4). Probability of non-diffractive interactions is Pnd=1-P0 -P1*P4P2-P3. If the target is a nucleus, P2 is set to zero, because it is not clear up to now how to deal with projectile diffraction on nuclei. This is implemented in lines 537 — 542. Due to quark exchange without participant excitation one or two final particles can be in Δ resonance state. The probability of this is set up in the line 543 : SetDeltaProbAtQuarkExchange( 0.0 );. Next lines: 544 — 548, set up a probability of the same quark's exchange between interacting hadrons.
    if ( NumberOfTargetNucleons > 26 ) {
      SetProbOfSameQuarkExchange( 1.0);
    } else {
      SetProbOfSameQuarkExchange( 0.0 );
    }
An exchange with the same quark in baryon-baryon interactions is expected to be suppressed at low energies. It cannot be so at high energies. The energy dependence of the parameter and its influence on final results is not yet well studied. Using the probability of the same quark's exchange, one can obtain the yield to the elastic scattering cross section. Thus, self-consistent change of the parametrization of the elastic cross section can be needed.

3.5.9.3. Parameters of participating hadron excitations

In the lines 549 — 555 the parameters of the excitations are set up.
    SetProjMinDiffMass( 1.16 );     // GeV
    SetProjMinNonDiffMass( 1.16 );  // GeV
    SetTarMinDiffMass( 1.16 );      // GeV
    SetTarMinNonDiffMass( 1.16 );   // GeV
    SetAveragePt2( 0.15 );          // GeV^2
    SetProbLogDistrPrD( 0.3 );
    SetProbLogDistr( 0.3 );
The first line determines the string mass sampling threshold (see Figure 3.12) in diffractive processes. For baryons, it is 1160=mN+mπ + 80 [MeV]. The second line is the analogous threshold for non-diffractive processes. The third and fourth lines set up analogous values for the target nucleon. A change of the parameters will lead to a change of the threshold behaviour of mass distribution in the reaction p+p → p+X. An average transverse momentum in the excitation is set up in the fifth line. Its change will change behaviour of particle distributions at large xF. For pure diffractive processes the average is determined in G4ElasticHNScattering as for the elastic scattering. It is known that the slope of pT distribution in diffraction processes (p+p → p+X, for example) depends on the produced mass of system X. Here an improvement can be introduced. Parameters set up in lines sixth and seventh are very important for a correct description of produced particles multiplicity in non-diffractive interactions. It is assumed in the code that mass distribution of strings produced in non-diffractive interactions has the form:
α can be different for projectile and target hadrons. The sixth line sets up the parameter for a projectile hadron, and the seventh line for a target hadron. A variation of the parameters can have an essential influence on multiplicity of produced particles. The analogous lines are presented for each type of projectile hadron.

3.5.9.4. Nuclear destruction parameters

The nuclear destruction parameters are set up in the lines 672 — 717. For a meson projectile they are:
 if ( ProjectileabsPDGcode < 1000 ) {  // Meson projectile
   SetMaxNumberOfCollisions( Plab, 2.0 );
   SetCofNuclearDestruction( 1.0*
           G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
   SetR2ofNuclearDestruction( 1.5*fermi*fermi );
   SetDofNuclearDestruction( 0.3 );
   SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
                                         ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
   SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
   SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
 }
The second line introduces the maximum number of intra-nuclear collisions used as a correction of the Abramovski-Gribov-Kancheli (AGK) cutting rules (see Geant4 Physics Reference Manual). The FTF model uses the reggeon cascading for simulation of fast nucleon production in hadron-nucleus interactions. As known, the Glauber approximation used in the Fritiof model and in the other string models does not provide enough amount of intra-nuclear collisions for a correct description of nuclear destruction. Additional cascading in nuclei is therefore needed. The usage of standard cascade for secondary particle interactions would lead to a too large multiplicity of produced particles. Thus, in Ref. [ RTIM1 ][ RTIM2 ] it was proposed to use reggeon cascading in the impact parameter space. For its description it is needed to determine the probability to involve a nuclear nucleons into the "cascade". It is obvious that the probability depends on the difference of impact coordinates of new and previous involved nucleons. Looking at the yield of reggeon enhanced diagrams, the functional form of the probability was chosen as:
where i and j are projections of the radii of ith and jth nucleons on the impact parameter plane. The "cascade" is initiated by primary involved nucleons. These nucleons are determined with the help of the Glauber approach. All involved nucleons are ejected from nuclei. The cascade looks like that: a projectile particle interacts with some intra-nuclear nucleons. These nucleons are called "wounded" or "participating" nucleons. The nucleons initiate the cascade. A wounded nucleon can involve a spectator nucleon into the cascade with the probability P. The latter one can involve another nucleon. The second nucleon can involve a third one and so on. This algorithm is implemented in the FTF model. We have tuned Cnd using the HARP-CDP data on proton production in the p+Cu interactions [ HARP_CDP ]. According to our estimations,
y is the projectile rapidity. The value, 2.1 standing in the exponents, corresponds to Plab ∼ 4 GeV/c, which is the middle of the transition region between low and high energies. The formula is implemented in the third and fourth lines. Its parameters, of course, can be changed to improve description of the fast nucleon production in the target fragmentation region. The formula is implemented in the third and fourth lines. Its parameters, of course, can be changed to improve description of the fast nucleon production in the target fragmentation region. The parameters Cnd and Rc are strongly correlated. To avoid this, Rc was fixed (fifth line). Momentum distribution of wounded and involved nucleons of a nucleus is sampled according to the distribution:
where MN is the multiplicity of involved nucleons, x is the light-cone momentum fraction carried by the nucleon relative to light-cone momentum of the nucleus, d is the dispersion parameter responsible for a spread of the momentum distribution, T is the transverse momentum, <PT2> is the average PT squared. d is setting up in the sixth line. Possible value of d can be 0.01 — 0.4. <PT2> is determined in the lines 7, 8. Again, an expression giving asymptotically constant value is used:
<PT2> must be small at low energies, and grows up as energy increase. All coefficients in the expression can be tuned. Their preliminary values were determined from the HARP-CDP experimental data. The last, tenth line sets up the yield of an involved nucleon to the excitation energy of a nuclear residual in MeV. The commented lines 719 — 752 are needed for testing purposes, for easily setting up of the parameters.

3.5.10. Changing parameters of an existing model: Lund use-case

The main parameters of the Geant4 LUND string fragmentation model are in the file G4LundStringFragmentation.cc located in geant4/source/processes/hadronic/models/parton\_string/hadronization/src. The main method of the class is FragmentString which is called by the class G4ExcitedStringDecay. The input parameter of the method is an object type G4ExcitedString. At the beginning of this method it is checked that the string can fragment (line 537, !IsFragmentable(&aString)). If the mass of the string does not allow a two particles decay, then the string is considered as a single hadron out of mass-shell. The minimal mass of the string that can decay is determined in the method SetMinimalStringMass (line 413). If the string can fragment, the main routine starts to work — Loop\_toFragmentString. The loop is working until it is decided to stop fragmentation of the string in the method StopFragmenting (lines 643 — 655). There are two important parameters in the method:
G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
{
  SetMinimalStringMass(string);
  if (string->FourQuarkString())
  {
    return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
  } else {
    return G4UniformRand() < G4Exp(-0.66e-6*(string->Mass()*string->Mass() -
                                             MinimalStringMass*MinimalStringMass));
  }
}
"0.0005" and "0.66e-6". The first one is responsible for the energy dependence of cross sections of the processes: ¯p + p → ¯q¯q + qq → h1 + h2. The second one determines the energy dependence of cross sections for hadronic interactions at low energies. If the fragmentation is not stopped, the code is trying to create a hadron splitting the string (line 1033, G4KineticTrack * Hadron=Splitup(currentString,newString)). The method "G4VLongitudinalStringDecay::Splitup" samples a side of the string decay: left or right string ends. Depending on type of a quark on the decaying end, a quark splitting, or a di-quark splitting is called. After them, a produced hadron momentum is determined in the method SplitEandP. The method QuarkSplitup is implemented in G4VLongitudinalStringDecay class. There can be two possibilities: q → h(q¯q') + q' and q → B(q (qq)') + (¯q¯q)'. A quark-antiquark pair is produced from the vacuum in the first case, and a meson h appears. A di-quark-anti-diquark pair is produced from the vacuum in the second case, and a baryon B appears. The probability of the second case is given by the parameter DiquarkSuppress} inherited from the parent class G4VLongitudinalStringDecay. It is reset in G4LundStringFragmentation.cc (line 69, SetDiquarkSuppression(0.05)). The probability of the first case is 1-DiquarkSuppress. The pairs created in the first case can consist of light quarks (u¯u, d¯d), or strange quarks (s¯s). The probability of the last possibility is StrangeSuppress/(2+StrangeSuppress). The parameter - StrangeSuppress, is set up in G4LundStringFragmentation.cc (line 395, SetStrangenessSuppression(0.46)). By increasing this parameter the strange meson production increases. At creation of the diquark-anti-diquark pair, the flavors of both quarks are sampled independently. Thus, the probability to produce a strange baryon in central region is roughly proportional to StrangeSuppress. The created meson can be a pseudo-scalar, or pseudo-vector one. The probability of a pseudo-scalar meson creation (pspin_meson) is set up in the parent class - G4VLongitudinalStringDecay (see line 82 in G4VLongitudinalStringDecay.cc). The probability of a pseudo-vector meson creation is 1-pspin_meson. The created baryon can belong to triplet or octet multiplets. The probability of the triplet state is given by the parameter pspin_barion determined in the parent class G4VLongitudinalStringDecay (line 85 of G4VLongitudinalStringDecay.cc). It is foreseen the possibility to change pspin_barion in the G4LundStringFragmentation.cc (line 176, //pspin_barion=0.75;).

It is not recommended to change pspin_meson and pspin_barion!

In the method G4LundStringFragmentation::DiQuarkSplitup two possibilities are considered: q1q2 → h(q1 ¯q') + q'q2 and q1q2 → B(q1 q2q') + ¯q'. The probability of the first process is regulated by the parameter DiquarkBreakProb. It is set up in G4LundStringFragmentation.cc (line 67, SetDiquarkBreakProbability(0.05)). This parameter is responsible for the meson production in hadronic interactions at large values |xF|.

The StrangeSuppress parameter is changed for the second process (see line 1487, StrangeSuppress=0.43).

After the creation of a hadron its momentum is determined in the method SplitEandP which has as input the type of the produced hadron and the string, and as input-output the residual string left after the hadron momentum sampling. At the beginning, only the quark content of the residual string is known. If the hadron mass plus a minimal mass of the residual string is larger than the mass of the initial string, then the sampling of the hadron is repeated as described above. If it is not so, an attempt to determine the hadron momentum is undertaken.

First of all, the transverse momentum of the hadron is sampled. The transverse momentum distribution is given by the formula:
If the sum of the transverse mass of the hadron and the transverse mass of the residual string is larger than the initial string mass, then the attempt is repeated. If it is not so, the longitudinal momentum of the hadron is sampled with the help of the GetLightConeZ method (lines 928 — 981). The LUND string fragmentation functions is implemented in GetLightConeZ. The method gives light-cone-momentum fraction of the initial string carried away by the hadron.
a=0.7 [GeV-2] for quark fragmentation, and for a production of a meson by a di-quark (line 941). The di-quark fragmentation function is assumed to be a flat distribution between zmin and zmax (line 978). When the fragmentation of a string is stopped, SplitLast method (lines 658 — 753) is called. Depending of the quark content of the string, Diquark_AntiDiquark_belowThreshold_lastSplitting, or Diquark_AntiDiquark_aboveThreshold_lastSplitting, or Quark_AntiQuark_lastSplitting, or Quark_Diquark_lastSplitting is called. The methods determine the types of the two last produced particles. After that, 4-momentum of the particles are calculated (Sample4Momentum, line 744). If at least one of the two final hadrons is a baryon, an isotropic decay is simulated with a probability ProbIsotropy (line 764),
The power 4 and functional form of the dependence can be tuned. It can reflect on angular distributions of baryons in annihilation processes, or in low energy interactions. If anisotropic decay happens, then the pT of final hadrons is simulated as before. For the sampling of the two final hadrons special arrays are introduced in the lines 81 — 375.


[1] The same can be achieved with template specialisations with slightly improved CPU performance but at the cost of significantly more complex designs and, with present compilers, significantly reduced portability.