Hadron-Level Standalone

The key trick is to set the flag ProcessLevel:all = off. When pythia.next() is called it then does not try to generate a hard process. Since there are no beams, it is also not possible to perform the normal PartonLevel step. (It is still possible to generate final-state radiation, but this is not automatic. It would have to be done by hand, using the pythia.forceTimeShower(...) method, before pythia.next() is called.) Thus only the HadronLevel methods are called, to take the current content of the event record stored in pythia.event as a starting point for any hadronization and decays that are allowed by the normal parameters of this step. Often the input would consist solely of partons grouped into colour singlets, but also (colour-singlet) particles are allowed.

To set up all the parameters, a pythia.init() call has to be used, without any arguments. In brief, the structure of the main program therefore should be something like

 
  Pythia pythia;                               // Declare generator. 
  Event& event = pythia.event                  // Convenient shorthand. 
  pythia.readString("ProcessLevel:all = off"); // The trick! 
  pythia.init();                               // Initialization. 
  for (int iEvent = 0; iEvent < nEvent; ++iEvent) { 
    // Insert filling of event here! 
    pythia.next();                             // Do the hadron level. 
  } 

Input configuration

With the event.append(...) methods a new entry is added at the bottom of the current record, i.e. the first time it is called entry number 1 is filled, and so on. The append method basically exists in four variants, either without or with history information, and with four-momentum provided either as a Vec4 four-vector or as four individual components:

 
  append( id, status, col, acol, p, m) 
  append( id, status, col, acol, px, py, pz, e, m) 
  append( id, status, mother1, mother2, daughter1, daughter2, col, acol, p, m) 
  append( id, status, mother1, mother2, daughter1, daughter2, col, acol, px, py, pz, e, m) 

All the four methods have two final, optional arguments. The scale one is highly relevant if you want to perform parton showers in addition to hadronization; see the pythia.forceTimeShower() description below. The final pol one denotes polarization, and could be used to perform polarized tau decays.

You can find descriptions of the input variables here. The PDG particle code id and the Les Houches Accord colour col and anticolour acol tags must be set correctly. The four-momentum and mass have to be provided in units of GeV; if you omit the mass it defaults to 0.

Outgoing particles that should hadronize should be given status code 23. Often this is the only status code you need. You could e.g. also fill in incoming partons with -21 and intermediate ones with -22, if you so wish. Usually the choice of status codes is not crucial, so long as you recall that positive numbers correspond to particles that are still around, while negative numbers denote ones that already hadronized or decayed. However, so as not to run into contradictions with the internal PYTHIA checks (when Check:event = on), or with external formats such as HepMC, we do recommend the above codes. When pythia.next() is called the positive-status particles that hadronize/decay get the sign of the status code flipped to negative but the absolute value is retained. The new particles are added with normal PYTHIA status codes.

For normal hadronization/decays in pythia.next() the history encoded in the mother and daughter indices is not used. Therefore the first two append methods, which set all these indices vanishing, should suffice. The subsequent hadronization/decays will still be properly documented.

The exception is when you want to include junctions in your string topology, i.e. have three string pieces meet. Then you must insert in your event record the (decayed) particle that is the reason for the presence of a junction, e.g. a baryon beam remnant from which several valence quarks have been kicked out, or a neutralino that underwent a baryon-number-violating decay. This particle must have as daughters the three partons that together carry the baryon number.

When ProcessLevel:all = off the pythia.next() call applied to a parton-level confguration will hadronize it, without generating any parton showers. Indeed, the point of the framework described here is to be able to feed in complete showered parton topologies for hadronization. (The exception is if you feed in a resonance, see next section.) As an option, however, it is possible to generate a shower before the pythia.next() step by using the pythia.forceTimeShower( int iBeg, int iEnd, double pTmax, int nBranchMax = 0) method. Here iBeg and iEnd give the range of partons that should be allowed to shower, pTmax the maximum pT scale of emissions, and a nonzero nBranchMax a maximum number of allowed branchings. Additionally, a scale has to be set for each parton that should shower, which requires an additional final argument to the append methods above, or alternatively separately with the pythia.event[i].scale(...) method. This scale limits the maximum pT allowed for each parton, in addition to the global pTmax. When not set the scale defaults to 0, meaning no radiation for that parton.

The sample program in main21.cc illustrates how you can work with this facility, both for simple parton configurations and for more complicated ones with junctions, and also how to force a shower.

As an alternative to setting up a topology with the methods above, a Les Houches Event File (LHEF) can also provide the configurations, using the "no-beams" extension. For parsing reasons the <init> and </init> tags need to be present as two separate lines, but there need not be anything between them. If there is, then the beam identities should be picked to be 0. A standard <LesHouchesEvents version="1.0"> line must also be at the top of the file. For the rest only the <event>....</event> blocks need to be present, one for each event. You should select Beams:frameType = 4 and provide the file name in Beams:LHEF, but setting ProcessLevel:all = off here is superfluous since the absence of beams is enough to make this apparent. Needless to say, an externally linked LHAup class works as well as an LHEF, with Beams:frameType = 5.

The event information to store in the LHEF, or provide by the LHAup, is essentially the same as above. The only difference is in status codes: outgoing particles should have 1 instead of 23, and intermediate resonances 2 instead of -22. Incoming partons, if any, are -1 instead of -21.

Extensions to resonance decays

The default when a resonance is encountered is to decay it, let the decay products shower, and finally hadronize the partons. Should a decay sequence already be provided at input, this sequence will be used as input for the showers, which are handled consecutively, followed by hadronization. Thus, a Higgs could be provided alone, or decaying to a pair of W bosons, or the same with the W's decaying further to fermion pairs. Needless to say, correct process-specific angular correlations in decays should not be expected when the process is unspecified.

If you do not want resonances to decay then you can use the ProcessLevel:resonanceDecays = off setting. If instead you want them to decay but not shower, you can use either PartonLevel:FSR = off or PartonLevel:FSRinResonances = off. A warning here is that, generally, it is not a good idea to provide a part of the shower history but not all, e.g. Z^0 → q qbar g: it is not straightforward to avoid doublecouning or other problems within this simpler alternative to a full-scale event generation.

The input configuration has to follow the rules described above, i.e. ProcessLevel:all = off should be set for internal input, but is not necessary for LHEF input. It is possible to combine several resonances, and other coloured or uncoloured particles into the same event. Partonic daughters of resonances would then shower, but other partons not. It may be possible to fool the program, however, since this is not a fully tested core functionality, so don't combine wildly if there is no reason to.

Repeated hadronization or decay

This method truly sticks to the hadron level, and thus cannot handle resonance decays. The real application instead is for repeated hadronization of the same PYTHIA process- and parton-level event. This may for some studies help to save time, given that these two first step are more time-consuming than the hadronization one.

For repeated hadronization you should first generate an event as usual, but with HadronLevel:all = off. This event you can save in a temporary copy, e.g. Event savedEvent = pythia.event. Inside a loop you copy back with pythia.event = savedEvent, and call pythia.forceHadronLevel() to obtain a new hadronization history.

A more limited form of repetition is if you want to decay a given kind of particle repeatedly, without having to generate the rest of the event anew. This could be the case e.g. in B physics applications. Then you can use the pythia.moreDecays() method, which decays all particles in the event record that have not been decayed but should have been done so. The pythia.particleData.mayDecay( id, false/true) method may be used to switch off/on the decays of a particle species id, so that it is not decayed in the pythia.next() call but only inside a loop over a number of tries.

Between each loop the newly produced decay products must be removed and the decayed particle status restored to undecayed. The former is simple, since the new products are appended to the end of the event record: event.saveSize() saves the initial size of the event record, and event.restoreSize() can later be used repeatedly to restore this original size, which means that the new particles at the end are thrown away. The latter is more complicated, and requires the user to identify the positions of all particles of the species and restore a positive status code with event[i].statusPos().

The main15.cc program illustrates both these methods, i.e. either repeated hadronization or repeated decay of PYTHIA events.