Adding an External Process to PYTHIA

This document refers to adding an external process to PYTHIA 6.102 within ATGEN 2.7. From here you can also view commented versions of the routines atpyini and pyupev which can be used to add Drell-Yan W production to PYTHIA.

Adding an external process to PYTHIA is done in two steps:

  1. initialize the process
  2. process each event

Initialize the process

When ATGEN initializes PYTHIA, it must also initialize your external process. So, in atgen.cmz in the pythia subdirectory edit the file "atpyini" and add the following call 

call pyupin(pid, title, max_xsec)

where pid is a process id which is not already in use by PYTHIA, title refers to the name of the process and max_xsec is a double precision value giving the maximum possible cross section for this process in millibarns. 

eg.
call pyupin(195, 'external process 1', 1.333D-04)

In order to determine the maximum possible cross section, it is necessary to first run your process in PYTHIA! Set the maximum to a very low value and run your process. PYTHIA will issue a warning and reset the max every time it calculates an event with a weight too high for the max you have set. Take the last value PYTHIA chooses and enter it in as the maximum in this call. Then run PYTHIA again......

Process an Event

The event processing level is where you provide PYTHIA with all of the physics it needs to create the event. In ATGEN's pythia directory, edit the file "pyupev". This routine is called every event in the following way 

call pyupev(pid,sigev)

where pid is the process id being requested and sigev is the weight (cross section) you will calculate in the routine and return to PYTHIA. The external routine will not need to take any information from PYTHIA, ALL of the physics is controlled from this routine and the information is passed to PYTHIA for further processing. By this I mean that the user specifies the incoming partons, calculates the matrix element that describes the interaction, produces the 4-vectors of the final state particle(s), determines the status of each particle (will it decay further?) and gives this information to PYTHIA without needing to use any internal PYTHIA information.

In order to calculate sigev (the event weight) in the event routine (pyupev), it is necessary to calculate the matrix element corresponding to the Feynman diagram for the process and to calculate the phase space factor for that event. Once sigev has been calculated, PYTHIA uses sigev/max_xsec to determine whether or not to keep the event. PYTHIA generates a random number between 0 and 1 and if it is less than sigev/max_xsec then the event is good, otherwise it calls pyupev again and generates another event.

In addition to calculating sigev the pyupev routine must fill a common block with particle information for the event. For example, PYTHIA must be given the four vector for each particle produced and the decay status of the particle (you can tell PYTHIA not to decay the particle any further).