Semi-Internal Processes

1. The Cross Section Calculation
2. The Cross Section Class
3. Access to a process
4. Implementing an external phase-space generator

Should you decide to go ahead on your own, this page gives a brief summary how to do that. If you additionally want to introduce a new resonance species, with its own internal width calculations, you will find further instructions here. It is strongly recommended to shop around for a similar process that has already been implemented, and to use that existing code as a template. Look for processes with the same combinations of incoming flavours and colour flows, rather than the shape of the cross section itself. With a reasonable such match the task should be of medium difficulty, without it more demanding.

PYTHIA's internal phase-space generators are rather good at handling the phase space of 2 → 1 and 2 → 2 processes, are more primitive for 2 → 3 ones and do not at all address higher multiplicities. An option is therefore also provided for external phase-space generators to be used, which must then be encapsulated to inherit from PYTHIA's PhaseSpace base class (or one of its derivatives). The set of processes that can be implemented in this framework is therefore in principle unlimited, though the user must supply external phase-space generators for non-trivial 2 → 3 processes and all higher 2 → n multiplicities. Note, however, that the produced particles may be resonances, so it is possible to end up with bigger "final" multiplicities through sequential decays, also with the internal phase-space generators, and to include further matrix-element weighting in those decays.

For processes using PYTHIA's internal phase-space generators, there are three steps involved in implementing a process:

1. making use of the PYTHIA-provided kinematics information to calculate the relevant cross section,
2. writing a new class, where the matrix elements are implemented, including information on incoming and outgoing flavours and colours, and
3. making the process available.

For processes for which an external phase-space generator will be used, step 1 above changes to writing a new class, where the phase-space generator is implemented, and making use of that to calculate the relevant cross section. There are no example programs illustrating how to do this yet, but the methodology is described below, under "Implementing an external phase-space generator".

The Cross Section Calculation

For a 2 → 1 process, the returned value should be sigmaHat(sHat), where mH (= mHat), sH (= sHat) and sH2 (= sHat^2) are available to be used. Incoming partons are massless. Overload the convertM2() method below if you instead plan to return |M|^2.

For a 2 → 2 process, instead d(sigmaHat)/d(tHat) should be returned, based on provided mH, sH, sH2, tH, tH2, uH, uH2, m3, s3, m4, s4 and pT2 values (s3 = m3*m3 etc.). Incoming partons are massless. Overload the convertM2() method below if you instead plan to return |M|^2.

For a 2 → 3 process, instead |M|^2 should be returned, with normalization such that |M|^2 / (2 sHat) integrated over the three-body phase space gives the cross section. Here no standard set of Mandelstam-style variables exists. Instead the obvious ones, mH, sH, m3, s3, m4, s4, m5, s5, are complemented by the four-vectors p3cm, p4cm, p5cm, from which further invariants may be calculated. The four-vectors are defined in the CM frame of the subcollision, with massless incoming partons along the +-z axis.

In either case, alpha_s and alpha_em have already been calculated, and are stored in alpS and alpEM. Also other standard variables may be used, like CoupEW::sin2thetaW(), and related flavour-dependent vector and axial couplings in CoupEW and CKM combinations in VCKM.

In case some of the final-state particles are resonances, their squared masses have already been selected according to a Breit-Wigner with a linearly running width Gamma(m) = Gamma(m_0) * m / m_0. More precisely, the mass spectrum is weighted according to w_BW(m^2) d(m^2), where
w_BW(m^2) = (1/pi) * (m * Gamma(m)) / ( (m^2 - m_0^2)^2 + (m * Gamma(m))^2 ) .
If you would like to have another expression, the above weights are stored in runBW3, runBW4 and runBW5, respectively. If you divide out one of these factors, you just remain with a phase space selection d(m^2) for this particle, and can multiply on your desired shape factor instead. Unfortunately, the Monte Carlo efficiency will drop if your new mass distribution differs dramatically from the input one. Therefore it does make sense to adjust the database value of the width to be slightly (but not too much) broader than the distribution you have in mind. Also note that, already by default, the wings of the Breit-Wigner are oversampled (with a compensating lower internal weight) by partly sampling like (a + b/m^2 + c/m^4) d(m^2), where the last term is only used for gamma^*/Z^0.

As alternative to the kinematics variables defined above, also the two arrays mME[5] and pME[5], for masses and four-momenta, respectively, can be used for cross-section calculations. Here indices 0 and 1 are the two incoming beams, and 2 and onwards the outgoing particles. Note that this differs by one step from the normal internal labeling, where slot 0 is left empty. The four-momenta are defined in the rest frame of the subcollision, with the incoming partons along the +-z direction. The kinematics need not agree with the "correct" one stored in the event record, for three reasons.
1) Gauge invariance forces matrix-element calculations to use the same masses for incoming as outgoing legs of a particle species, say b quarks. Therefore the kinematics of the two incoming partons is recalculated, relative to the normal event record, to put the partons on the mass shell. (Note that initial masses is a technical issue, not the correct physics picture: the incoming partons are likely to be spacelike virtual rather than on the mass shell.)
2) In principle each fermion flavour has to be treated separately, owing to a different mass. However, in many cases fermions can be assumed massless, which speeds up the calculations, and further gains occur if then different flavours can use the same cross-section expression. In MadGraph the default is that fermions up to and including the c quark and the mu lepton are considered massless, while the b quark and the tau lepton are considered massive. This can be modified however, and below we provide four flags that can be used to consider the "borderline" fermions either as massless or as massive when matrix elements are evaluated, to match the assumptions made for the matrix elements themselves.
3) For 2 → 2 and 2 → 3 processes of massive identical particles (or antiparticles) in the final state, such as t tbar or W^+ W^-, the kinematics is here adjusted so that the two or three particles have the same mass, formed as a suitable average of the actual Breit-Wigner-distributed masses. This allows the evaluation of matrix-element expressions that only have meaning if the two/three have the same mass.
Thus the mass array mME[5] and the four-momentum array pME[5] present values both for initial- and final-state particles based on these mass principles suited for matrix-element input. Note that these variables therefore differ from the kinematics stored in the event record proper, where incoming fermions are always massless and outgoing resonances have independent Breit-Wigner mass distributions.
The conversion from the normal to the special kinematics is done by calling the setupForME() method. This you have to do yourself in the SigmaHat() member of your derived class. Alternatively it could be done in SigmaKin(), i.e. before the loop over incoming flavours, but then these would be considered massless. The identity of final-state particles is obtained from the id3Mass(), id4Mass() and id5Mass() methods. Should the conversion to mME[5] and pME[5] not work, setupForME() will return false, and then the cross section should be put zero.

flag  SigmaProcess:cMassiveME   (default = off)
Let the c quark be massive or not in the kinematics set up for external matrix-element evaluation.

flag  SigmaProcess:bMassiveME   (default = on)
Let the b quark be massive or not in the kinematics set up for external matrix-element evaluation.

flag  SigmaProcess:muMassiveME   (default = off)
Let the mu lepton be massive or not in the kinematics set up for external matrix-element evaluation.

flag  SigmaProcess:tauMassiveME   (default = on)
Let the tau lepton be massive or not in the kinematics set up for external matrix-element evaluation.

The Cross Section Class

The class has to be derived either from Sigma1Process, for 2 → 1 processes, from Sigma2Process, for 2 → 2 ones, or from Sigma3Process, for 2 → 3 ones. (The Sigma0Process class is used for elastic, diffractive and minimum-bias events, and is not recommended for use beyond that.) These are in their turn derived from the SigmaProcess base class.

The class can implement a number of methods. Some of these are compulsory, others strongly recommended, and the rest are to be used only when the need arises to override the default behaviour. The methods are:

A constructor for the derived class obviously must be available. Here you are quite free to allow a list of arguments, to set the parameters of your model, or even to create a set of closely related but distinct processes. For instance, g g → Q Qbar, Q = c or b, is only coded once, and then the constructor takes the quark code (4 or 5) as argument, to allow the proper amount of differentiation.

A destructor is only needed if you require some special behaviour.

void SigmaProcess::initProc()  
is called once during initialization, and can then be used to set up parameters, such as masses and couplings, and perform calculations that need not be repeated for each new event, thereby saving time. This method needs not be implemented, since in principle all calculations can be done in sigmaHat below.

void SigmaProcess::sigmaKin()  
is called once a kinematical configuration has been determined, but before the two incoming flavours are known. This routine can therefore be used to perform calculations that otherwise might have to be repeated over and over again in sigmaHat below. For instance a flavour-independent cross section calculation for a q g initial state would be repeated 20 times in sigmaHat, five times for the five quark flavours allowed in the incoming beams, times twice to include antiquarks, times twice since the (anti)quark could be in either of the two beams. You could therefore calculate the result once only and store it as a private data member of the class. It is optional whether you want to use this method, however, or put everything in sigmaHat.

double SigmaProcess::sigmaHat()  
is the key method for cross section calculations and returns a cross section value, as described in the previous section. It is called when also a preliminary set of incoming flavours has been picked, in addition to the kinematical ones already available for sigmaKin. Typically sigmaHat is called inside a loop over all allowed incoming flavour combinations, stored in id1 and id2, with fixed kinematics, as already illustrated above. The sum over the different flavour combinations provides the total cross section, while their relative size is used to make a selection of a specific incoming state.

bool SigmaProcess::setupForME()  
to be called by the user from inside sigmaHat() (or possibly sigmaKin()) to setup alternative kinematics in the mME[5] and pME[5] arrays, better suited for matrix-element calculations. See the end of the previous section for a more detailed description. Should the method return false then the conversion did not work, and sigmaHat() (or sigmaKin()) should be set to vanish.

void SigmaProcess::setIdColAcol()  
is called only once an initial state and a kinematical configuration has been picked. This routine must set the complete flavour information and the colour flow of the process. This may involve further random choices, between different possible final-state flavours or between possible competing colour flows. Private data members of the class may be used to retain some information from the previous steps above.
When this routine is called the two incoming flavours have already been selected and are available in id1 and id2, whereas the one, two or three outgoing ones either are fixed for a given process or can be determined from the instate (e.g. whether a W^+ or W^- was produced). There is also a standard method in VCKM to pick a final flavour from an initial one with CKM mixing. Once you have figured out the value of id3 and, the case being, id4 and id5, you store these values permanently by a call setId( id1, id2, id3, id4, id5), where the last two may be omitted if irrelevant.
Correspondingly, the colours are stored with setColAcol( col1, acol1, col2, acol2, col3, acol3, col4, acol4, col5, acol5), where the final ones may be omitted if irrelevant. Les Houches style colour tags are used, but starting with number 1 (and later shifted by the currently requested offset). The input is grouped particle by particle, with the colour index before the anticolour one. You may need to select colour flow dynamically, depending on the kinematics, when several distinct possibilities exist. Trivial operations, like swapping colours and anticolours, can be done with existing methods.
When the id3Mass() and id4Mass() methods have been used, the order of the outgoing particles may be inconsistent with the way the tHat and uHat variables have been defined. A typical example would be a process like q g → q' W with tHat defined between incoming and outgoing quark, but where id3Mass() = 24 and so the process is to be stored as q g → W q'. One should then put the variable swapTU = true in setIdColAcol() for each event where the tHat and uHat variables should be swapped before the event kinematics is reconstructed. This variable is automatically restored to false for each new event.

double SigmaProcess::weightDecayFlav( Event& process)  
is called to allow a reweighting of the simultaneous flavour choices of resonance decay products. Is currently only used for the q qbar → gamma*/Z^0 gamma*/Z^0 process, and will likely not be of interest for you.

double SigmaProcess::weightDecay( Event& process, int iResBeg, int iResEnd)  
is called when the basic process has one or several resonances, after each set of related resonances in process[i], iResBeg <= i <= iResEnd, has been allowed to decay. The calculated weight, to be normalized to the range between 0 and 1, is used to decide whether to accept the decay(s) or try for a new decay configuration. The base-class version of this method returns unity, i.e. gives isotropic decays by default. This method may be called repeatedly for a single event. For instance, in q qbar → H^0 Z^0 with H^0 → W^+ W^-, a first call would be made after the H^0 and Z^0 decays, and then depend only on the Z^0 decay angles since the H^0 decays isotropically. The second call would be after the W^+ W^- decays and then involve correlations between the four daughter fermions.

string SigmaProcess::name()  
returns the name of the process, as you want it to be shown in listings.

int SigmaProcess::code()  
returns an integer identifier of the process. This has no internal function, but is only intended as a service for the user to rapidly (and hopefully uniquely) identify which process occurred in a given event. Numbers below 10000 are reserved for internal PYTHIA use.

string SigmaProcess::inFlux()  
this string specifies the combinations of incoming partons that are allowed for the process under consideration, and thereby which incoming flavours id1 and id2 the sigmaHat() calls will be looped over. It is always possible to pick a wider flavour selection than strictly required and then put to zero cross sections in the superfluous channels, but of course this may cost some extra execution time. Currently allowed options are:
* gg: two gluons.
* qg: one (anti)quark and one gluon.
* qq: any combination of two quarks, two antiquarks or a quark and an antiquark.
* qqbar: any combination of a quark and an antiquark; a subset of the qq option.
* qqbarSame: a quark and its antiquark; a subset of the qqbar option.
* ff: any combination of two fermions, two antifermions or a fermion and an antifermion; is the same as qq for hadron beams but also allows processes to work with lepton beams.
* ffbar: any combination of a fermion and an antifermion; is the same as qqbar for hadron beams but also allows processes to work with lepton beams.
* ffbarSame: a fermion and its antifermion; is the same as qqbarSame for hadron beams but also allows processes to work with lepton beams.
* ffbarChg: a fermion and an antifermion that combine to give charge +-1.
* fgm: a fermion and a photon (gamma).
* ggm: a gluon and a photon.
* gmgm: two photons.

bool SigmaProcess::convert2mb()  
it is assumed that cross sections normally come in dimensions such that they, when integrated over the relevant phase space, obtain the dimension GeV^-2, and therefore need to be converted to mb. If the cross section is already encoded as mb then convert2mb() should be overloaded to instead return false.

bool SigmaProcess::convertM2()  
it is assumed that 2 → 1 cross sections are encoded as sigmaHat(sHat), and 2 → 2 ones as d(sigmaHat)/d(tHat) in the SigmaProcess::sigmaHat() methods. If convertM2() is overloaded to instead return true then the return value is instead assumed to be the squared matrix element |M|^2, and SigmaProcess::sigmaHatWrap(...) converts to sigmaHat(sHat) or d(sigmaHat)/d(tHat), respectively. This switch has no effect on 2 → 3 processes, where |M|^2 is the only allowed input anyway.

int SigmaProcess::id3Mass()  
int SigmaProcess::id4Mass()  
int SigmaProcess::id5Mass()  
are the one, two or three final-state flavours, where masses are to be selected before the matrix elements are evaluated. Only the absolute value should be given. For massless particles, like gluons and photons, one need not give anything, i.e. one defaults to 0. The same goes for normal light quarks, where masses presumably are not implemented in the matrix elements. Later on, these quarks can still (automatically) obtain constituent masses, once a u, d or s flavour has been selected.

int SigmaProcess::resonanceA()  
int SigmaProcess::resonanceB()  
are the codes of up to two s-channel resonances contributing to the matrix elements. These are used by the program to improve the phase-space selection efficiency, by partly sampling according to the relevant Breit-Wigner distributions. Massless resonances (the gluon and photon) need not be specified.

bool SigmaProcess::isSChannel()  
normally the choice of renormalization and factorization scales in 2 → 2 and 2 → 3 processes is based on the assumption that t- and u-channel exchanges dominates the cross section. In cases such as f fbar → gamma* → f' fbar' a 2 → 2 process actually ought to be given scales as a 2 → 1 one, in the sense that it proceeds entirely through an s-channel resonance. This can be achieved if you override the default false to return true. See further the page on couplings and scales.

int SigmaProcess::idSChannel()  
normally no intermediate state is shown in the event record for 2 → 2 and 2 → 3 processes. However, in case that idSChannel is overloaded to return a nonzero value, an intermediate particle with that identity code is inserted into the event record, to make it a 2 → 1 → 2 or 2 → 1 → 3 process. Thus if both isSChannel and idSChannel are overloaded, a process will behave and look like it proceeded through a resonance. The one difference is that the implementation of the matrix element is not based on the division into a production and a decay of an intermediate resonance, but is directly describing the transition from the initial to the final state.

int SigmaProcess::isQCD3body()  
there are two different 3-body phase-space selection machineries, of which the non-QCD one is default. If you overload this method instead the QCD-inspired machinery will be used. The differences between these two is related to which phase space cuts can be set, and also that the QCD machinery assumes (almost) massless outgoing partons.

int SigmaProcess::idTchan1()  
int SigmaProcess::idTchan2()  
the non-QCD 2 → 3 phase space selection machinery is rather primitive, as already mentioned. The efficiency can be improved in processes that proceed though t-channel exchanges, such as q qbar' → H^0 q qbar' via Z^0 Z^0 fusion, if the identity of the t-channel-exchanged particles on the two side of the event are provided. Only the absolute value is of interest.

double SigmaProcess::tChanFracPow1()  
double SigmaProcess::tChanFracPow2()  
in the above kind of 2 → 3 phase-space selection, the sampling of pT^2 is done with one part flat, one part weighted like 1 / (pT^2 + m_R^2) and one part like 1 / (pT^2 + m_R^2)^2. The above values provide the relative amount put in the latter two channels, respectively, with the first obtaining the rest. Thus the sum of tChanFracPow1() and tChanFracPow2() must be below unity. The final results should be independent of these numbers, but the Monte Carlo efficiency may be quite low for a bad choice. Here m_R is the mass of the exchanged resonance specified by idTchan1() or idTchan2(). Note that the order of the final-state listing is important in the above q qbar' → H^0 q qbar' example, i.e. the H^0 must be returned by id3Mass(), since it is actually the pT^2 of the latter two that are selected independently, with the first pT then fixed by transverse-momentum conservation.

bool SigmaProcess::useMirrorWeight()  
in 2 → 3 processes the phase space selection used here involves a twofold ambiguity basically corresponding to a flipping of the positions of last two outgoing particles. These are assumed equally likely by default, false, but for processes proceeding entirely through t-channel exchange the Monte Carlo efficiency can be improved by making a preselection based on the relative propagator weights, true.

int SigmaProcess::gmZmode()  
allows a possibility to override the global mode WeakZ0:gmZmode for a specific process. The global mode normally is used to switch off parts of the gamma^*/Z^0 propagator for test purposes. The above local mode is useful for processes where a Z^0 really is that and nothing more, such as q qbar → H^0 Z^0. The default value -1 returned by gmZmode() ensures that the global mode is used, while 0 gives full gamma^*/Z^0 interference, 1 gamma^* only and 2 Z^0 only.

Access to a process

 
      SigmaProcessPtr mySigma = make_shared<MySigma()>; 
      pythia.setSigmaPtr( mySigma); 

 
      pythia.setSigmaPtr(make_shared<MySigma()>, 
                         make_shared<MyPhaseSpaceGenerator()>); 

If the code should be of good quality and general usefulness, it would be simple to include it as a permanently available process in the standard program distribution. The final step of that integration ought to be left for the PYTHIA authors, but here is a description of what is required.

A flag has to be defined, that allows the process to be switched on; by default it should always be off. The name of the flag should be chosen of the type model:process. Here the model would be related to the general scenario considered, e.g. Compositeness, while process would specify instate and outstate, separated by a 2 (= to), e.g. ug2u*g. When several processes are implemented and "belong together" it is also useful to define a model:all switch that affects all the separate processes.

The flags should normally be stored in the ProcessSelection.xml file or one of its daughters for a specific kind of processes. This is to make them easily found by users. You could create and use your own .xml file, so long as you then add that name to the list of files in the Index.xml file. (If not, the flags would never be created and the program would not work.)

In the ProcessContainer.c file, the SetupContainers::init() method needs to be expanded to create instances of the processes switched on. This code is fairly repetitive, and should be easy to copy and modify from the code already there. The basic structure is
(i) check whether a process is requested by the user and, if so,
(ii) create an instance of the matrix-element class,
(iii)create a container for the matrix element and its associated phase-space handling, and
(iv) add the container to the existing process list.

Two minor variations are possible. One is that a set of related processes are lumped inside the the same initial check, i.e. are switched on all together. The second is that the matrix-element constructor may take arguments, as specified by you (see above). If so, the same basic matrix element may be recycled for a set of related processes, e.g. one for a composite u and one for a composite d. Obviously these variations may be combined.

Implementing an external phase-space generator

 
      // Determine how phase space should be sampled. 
      virtual bool setupSampling(); 
      // Select a trial event kinematics. 
      virtual bool trialKin(bool inEvent = true, bool repeatSame = false); 
      // Construct final (accepted) event kinematics. 
      virtual bool finalKin(); 

 
      // Inform whether beam particles are resolved or scatter directly. 
      virtual bool isResolved(); 

In the setupSampling() step the main point is to determine the upper estimate of the cross section integrated over the allowed phase space regions, and this should be stored in sigmaMx. The ratio between the correct cross section and its upper estimate is a measure of the phase-space selection efficiency, and the purpose of this step is to optimize the sampling accordingly. To this end any convenient set of phase-space variables may be chosen. The x1H and x2H varables should be used to denote the incoming parton momentum fractions, however, to be used in PDF evaluations.

In the trialKin() intermediate step the same set of internal variables can be used, and fed into the SigmaProcess code to evaluate the cross section in the given phase space point, multiplied by the integrated cross section. This value is to be stored in sigmaNw, and the ratio sigmaNw/sigmaMx will be used to determine whether the trial event is accepted or not.

In the finalKin() step the output is more standardized. The key values are the ones stored in the mH[] and pH[] arrays, the former for masses and the latter for four-momenta. Here the first two slots represent the two incoming partons and the subsequent ones up to ten outgoing particles. Other particle properties, like the number of final-state particles, their identities and colours, and more, are defined by the SigmaProcess class.

A tailor-made 2 → 3 generator could be defined, e.g., by starting from the code for PYTHIA's internal PhaseSpace2to3tauycyl base class, which provides a specific representation of 3-parton phase space, used for generic 2 → 3 processes in PYTHIA. The virtual functions described could then be redefined to generate a different sampling of 3-parton phase space. One example of this is provided by the existing PhaseSpace2to3yyycyl class, which PYTHIA uses for massless QCD processes. Note the interplay between the phase-space variables, generated and saved here, and how they are used by the matrix-element codes. For general processes, the user can define samplings in terms of their own phase-space parametrizations, as long as the corresponding matrix elements use the same variables to evaluate the cross-section expressions.