Semi-Internal Processes

Should you actually go ahead, 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 is rather good at handling the phase space of 2 -> 1 and 2 -> 2 processes, is more primitive for 2 -> 3 ones and does not at all address higher multiplicities. This limits the set of processes that you can implement in this framework. The produced particles may be resonances, however, so it is possible to end up with bigger "final" multiplicities through sequential decays, and to include further matrix-element weighting in those decays.

There are two steps involved in implementing a process:
1) writing a new class, where the matrix elements are implemented, including information on incoming and outgoing flavours and colours, and
2) making the process available.
We consider these two aspects in turn. An example where it all comes together is found in main25.cc.

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 plan to delete the process before the natural end of the run, and require some special behaviour at that point. If you call such a destructor you will leave a pointer dangling inside the Pythia object you gave it in to, if that still exists.

void SigmaProcess::initProc()  
is called once during initalization, 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 further described below. 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 incomimg state.
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.
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.).
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 variables exist. 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 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.

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 occured 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.
* qqbarSame: a quark and its antiquark; this is a subset of the above qq 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.
* 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.

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-Wigners. 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::idTchan1()  
int SigmaProcess::idTchan2()  
the 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

      SigmaProcess* mySigma = new MySigma();
      pythia.setSigmaPtr( mySigma); 

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.