PYTHIA8/pythia8175/htmldoc/SemiInternalProcesses.html

   1 <html>
   2 <head>
   3 <title>Semi-Internal Processes</title>
   4 <link rel="stylesheet" type="text/css" href="pythia.css"/>
   5 <link rel="shortcut icon" href="pythia32.gif"/>
   6 </head>
   7 <body>
   8
   9 <h2>Semi-Internal Processes</h2>
  10
  11 Normally users are expected to implement new processes via the
  12 <a href="LesHouchesAccord.html" target="page">Les Houches Accord</a>. Then
  13 you do all flavour, colour and phase-space selection externally,
  14 before your process-level events are input for further processing
  15 by PYTHIA. However, it is also possible to implement a
  16 new process in exactly the same way as the internal PYTHIA
  17 ones, thus making use of the internal phase space selection machinery
  18 to sample an externally provided cross-section expression.
  19 The MadGraph 5 program [<a href="Bibliography.html" target="page">Alw11</a>] allows you to do exactly that,
  20 i.e. it can be used to generate C++ code that can be linked into
  21 the existing PYTHIA framework, see
  22 <a href="MadGraph5Processes.html" target="page">here</a>.
  23
  24 <p/>
  25 Should you decide to go ahead on your own,
  26 this page gives a brief summary how to do that. If you additionally
  27 want to introduce a new resonance species, with its own internal
  28 width calculations, you will find further instructions
  29 <a href="SemiInternalResonances.html" target="page">here</a>. It is strongly
  30 recommended to shop around for a similar process that has already
  31 been implemented, and to use that existing code as a template.
  32 Look for processes with the same combinations of incoming flavours
  33 and colour flows, rather than the shape of the cross section itself.
  34 With a reasonable such match the task should be of medium difficulty,
  35 without it more demanding.
  36
  37 <p/>
  38 PYTHIA is rather good at handling the phase space of
  39 <i>2 -> 1</i> and <i>2 -> 2</i> processes, is more primitive for
  40 <i>2 -> 3</i> ones and does not at all address higher multiplicities.
  41 This limits the set of processes that you can implement in this
  42 framework. The produced particles may be resonances, however, so it is
  43 possible to end up with bigger "final" multiplicities through sequential
  44 decays, and to include further matrix-element weighting in those decays.
  45
  46 <p/>
  47 There are three steps involved in implementing a process:
  48 <ol>
  49 <li>making use of the PYTHIA-provided kinematics information to
  50 calculate the relevant cross section,</li>
  51 <li>writing a new class,  where the matrix elements are implemented,
  52 including information on incoming and outgoing flavours and colours,
  53 and</li>
  54 <li>making the process available.</li>
  55 </ol>
  56 We consider these aspects in turn. An example where it all comes
  57 together is found in <code>main22.cc</code>.
  58
  59 <h3>The Cross Section Calculation</h3>
  60
  61 The key method for the cross section calculation is
  62 <code>SigmaProcess::sigmaHat()</code>, described below. At the point when
  63 it is called, the kinematics has already been set up, and from these
  64 phase space variables the differential cross section is to be calculated.
  65
  66 <p/>
  67 For a <i>2 -> 1</i> process, the returned value should be
  68 <i>sigmaHat(sHat)</i>, where <code>mH</code> (= <i>mHat</i>),
  69 <code>sH</code> (= <i>sHat</i>) and <code>sH2</code> (= <i>sHat^2</i>)
  70 are available to be used. Incoming partons are massless. Overload the
  71 <code>convertM2()</code> method below if you instead plan to return
  72 <i>|M|^2</i>.
  73
  74 <p/>
  75 For a <i>2 -> 2</i> process, instead <i>d(sigmaHat)/d(tHat)</i>
  76 should be returned, based on provided
  77 <code>mH, sH, sH2, tH, tH2, uH, uH2, m3, s3, m4, s4</code> and
  78 <code>pT2</code> values (<code>s3 = m3*m3</code> etc.). Incoming
  79 partons are massless. Overload the <code>convertM2()</code> method
  80 below if you instead plan to return <i>|M|^2</i>.
  81
  82 <p/>
  83 For a <i>2 -> 3</i> process, instead <i>|M|^2</i> should be
  84 returned, with normalization such that <i>|M|^2 / (2 sHat)</i> integrated
  85 over the three-body phase space gives the cross section. Here no standard
  86 set of Mandelstam-style variables exists. Instead the obvious ones,
  87 <code>mH, sH, m3, s3, m4, s4, m5, s5</code>, are complemented by the
  88 four-vectors <code>p3cm, p4cm, p5cm</code>, from which further invariants
  89 may be calculated. The four-vectors are defined in the CM frame of the
  90 subcollision, with massless incoming partons along the <i>+-z</i> axis.
  91
  92 <p/>
  93 In either case, <i>alpha_s</i> and <i>alpha_em</i> have already
  94 been calculated, and are stored in <code>alpS</code> and <code>alpEM</code>.
  95 Also other standard variables may be used, like
  96 <code>CoupEW::sin2thetaW()</code>, and related flavour-dependent
  97 vector and axial couplings in <code>CoupEW</code> and CKM combinations
  98 in <code>VCKM</code>.
  99
 100 <p/>
 101 In case some of the final-state particles are resonances, their
 102 squared masses have already been selected according to a Breit-Wigner
 103 with a linearly running width <i>Gamma(m) = Gamma(m_0) * m / m_0</i>.
 104 More precisely, the mass spectrum is weighted according to
 105 <i>w_BW(m^2) d(m^2)</i>, where
 106 <br/><i>
 107 w_BW(m^2) = (1/pi) * (m * Gamma(m)) / ( (m^2 - m_0^2)^2 + (m * Gamma(m))^2 ) .
 108 </i><br/>
 109 If you would like to have another expression, the above weights are stored
 110 in <code>runBW3</code>, <code>runBW4</code> and <code>runBW5</code>,
 111 respectively. If you divide out one of these factors, you just remain with
 112 a phase space selection <i>d(m^2)</i> for this particle,
 113 and can multiply on your desired shape factor instead. Unfortunately, the
 114 Monte Carlo efficiency will drop if your new mass distribution differs
 115 dramatically from the input one. Therefore it does make sense to adjust the
 116 database value of the width to be slightly (but not too much) broader
 117 than the distribution you have in mind. Also note that, already by default,
 118 the wings of the Breit-Wigner are oversampled (with a compensating lower
 119 internal weight) by partly sampling like <i>(a + b/m^2 + c/m^4) d(m^2)</i>,
 120 where the last term is only used for <i>gamma^*/Z^0</i>.
 121
 122 <p/>
 123 As alternative to the kinematics variables defined above, also the two
 124 arrays <code>mME[5]</code> and <code>pME[5]</code>, for masses and
 125 four-momenta, respectively, can be used for cross-section calculations.
 126 Here indices 0 and 1 are the two incoming beams, and 2 and onwards the
 127 outgoing particles. Note that this differs by one step from the normal
 128 internal labeling, where slot 0 is left empty. The four-momenta are
 129 defined in the rest frame of the subcollision, with the incoming partons
 130 along the <i>+-z</i> direction. The kinematics need not agree with the
 131 "correct" one stored in the event record, for three reasons.
 132 <br/>1) Gauge invariance forces matrix-element calculations to use
 133 the same masses for incoming as outgoing legs of a particle species,
 134 say <i>b</i> quarks. Therefore the kinematics of the two incoming
 135 partons is recalculated, relative to the normal event record, to put
 136 the partons on the mass shell. (Note that initial masses is a technical
 137 issue, not the correct physics picture: the incoming partons are likely
 138 to be spacelike virtual rather than on the mass shell.)
 139 <br/>2) In principle each fermion flavour has to be treated separately,
 140 owing to a different mass. However, in many cases fermions can be
 141 assumed massless, which speeds up the calculations, and further gains
 142 occur if then different flavours can use the same cross-section
 143 expression. In MadGraph the default is that fermions up to and including
 144 the <i>c</i> quark and the <i>mu</i> lepton are considered massless,
 145 while the <i>b</i> quark and the <i>tau</i> lepton are considered
 146 massive. This can be modified however, and below we provide four flags
 147 that can be used to consider the "borderline" fermions either as
 148 massless or as massive when matrix elements are evaluated, to match the
 149 assumptions made for the matrix elements themselves.
 150 <br/>3) For <i>2 -> 2</i> and <i>2 -> 3</i> processes of massive
 151 identical particles (or antiparticles) in the final state, such as
 152 <i>t tbar</i> or <i>W^+ W^-</i>, the kinematics is here adjusted
 153 so that the two or three particles have the same mass, formed as a
 154 suitable average of the actual Breit-Wigner-distributed masses. This
 155 allows the evaluation of matrix-element expressions that only have
 156 meaning if the two/three have the same mass.
 157 <br/>Thus the mass array <code>mME[5]</code> and the four-momentum array
 158 <code>pME[5]</code> present values both for initial- and final-state
 159 particles based on these mass principles suited for matrix-element input.
 160 Note that these variables therefore differ from the kinematics stored in
 161 the event record proper, where incoming fermions are always massless and
 162 outgoing resonances have independent Breit-Wigner mass distributions.
 163 <br/>The conversion from the normal to the special kinematics is done
 164 by calling the <code>setupForME()</code> method. This you have to do
 165 yourself in the <code>SigmaHat()</code> member of your derived class.
 166 Alternatively it could be done in <code>SigmaKin()</code>, i.e. before
 167 the loop over incoming flavours, but then these would be considered
 168 massless. The identity of final-state particles is obtained from the
 169 <code>id3Mass()</code>, <code>id4Mass()</code> and <code>id5Mass()</code>
 170 methods. Should the conversion to <code>mME[5]</code> and
 171 <code>pME[5]</code> not work, <code>setupForME()</code> will return
 172 <code>false</code>, and then the cross section should be put zero.
 173
 174 <p/><code>flag&nbsp; </code><strong> SigmaProcess:cMassiveME &nbsp;</strong>
 175  (<code>default = <strong>off</strong></code>)<br/>
 176 Let the <i>c</i> quark be massive or not in the kinematics set up for
 177 external matrix-element evaluation.
 178
 179
 180 <p/><code>flag&nbsp; </code><strong> SigmaProcess:bMassiveME &nbsp;</strong>
 181  (<code>default = <strong>on</strong></code>)<br/>
 182 Let the <i>b</i> quark be massive or not in the kinematics set up for
 183 external matrix-element evaluation.
 184
 185
 186 <p/><code>flag&nbsp; </code><strong> SigmaProcess:muMassiveME &nbsp;</strong>
 187  (<code>default = <strong>off</strong></code>)<br/>
 188 Let the <i>mu</i> lepton be massive or not in the kinematics set up for
 189 external matrix-element evaluation.
 190
 191
 192 <p/><code>flag&nbsp; </code><strong> SigmaProcess:tauMassiveME &nbsp;</strong>
 193  (<code>default = <strong>on</strong></code>)<br/>
 194 Let the <i>tau</i> lepton be massive or not in the kinematics set up for
 195 external matrix-element evaluation.
 196
 197
 198
 199 <h3>The Cross Section Class</h3>
 200
 201 The matrix-element information has to be encoded in a new class.
 202 The relevant code could either be put before the main program in the
 203 same file, or be stored separately, e.g. in a matched pair
 204 of <code>.h</code> and <code>.cc</code> files. The latter may be more
 205 convenient, in particular if the cross sections are lengthy, or if you
 206 intend to build up your own little process library, but of course
 207 requires that these additional files are correctly compiled and linked.
 208
 209 <p/>
 210 The class has to be derived either from
 211 <code>Sigma1Process</code>, for <i>2 -> 1</i> processes, from
 212 <code>Sigma2Process</code>, for <i>2 -> 2</i> ones, or from
 213 <code>Sigma3Process</code>, for <i>2 -> 3</i> ones. (The
 214 <code>Sigma0Process</code> class is used for elastic, diffractive
 215 and minimum-bias events, and is not recommended for use beyond that.)
 216 These are in their turn derived from the <code>SigmaProcess</code>
 217 base class.
 218
 219 <p/>
 220 The class can implement a number of methods. Some of these are
 221 compulsory, others strongly recommended, and the rest are to be
 222 used only when the need arises to override the default behaviour.
 223 The methods are:
 224
 225 <p/>
 226 A <b>constructor</b> for the derived class obviously must be available.
 227 Here you are quite free to allow a list of arguments, to set
 228 the parameters of your model, or even to create a set of closely
 229 related but distinct processes. For instance, <i>g g -> Q Qbar</i>,
 230 <i>Q = c</i> or <i>b</i>, is only coded once, and then the
 231 constructor takes the quark code (4 or 5)  as argument,
 232 to allow the proper amount of differentiation.
 233
 234 <p/>
 235 A <b>destructor</b> is only needed if you plan to delete the process
 236 before the natural end of the run, and require some special behaviour
 237 at that point. If you call such a destructor you will leave a pointer
 238 dangling inside the <code>Pythia</code> object you gave it in to,
 239 if that still exists.
 240
 241 <a name="method1"></a>
 242 <p/><strong>void SigmaProcess::initProc() &nbsp;</strong> <br/>
 243 is called once during initialization, and can then be used to set up
 244 parameters, such as masses and couplings, and perform calculations
 245 that need not be repeated for each new event, thereby saving time.
 246 This method needs not be implemented, since in principle all
 247 calculations can be done in <code>sigmaHat</code> below.
 248
 249
 250 <a name="method2"></a>
 251 <p/><strong>void SigmaProcess::sigmaKin() &nbsp;</strong> <br/>
 252 is called once a kinematical configuration has been determined, but
 253 before the two incoming flavours are known. This routine can therefore
 254 be used to perform calculations that otherwise might have to be repeated
 255 over and over again in <code>sigmaHat</code> below. For instance
 256 a flavour-independent cross section calculation for a <i>q g</i>
 257 initial state would be repeated 20 times in <code>sigmaHat</code>,
 258 five times for the five quark flavours allowed in the incoming beams,
 259 times twice to include antiquarks, times twice since the (anti)quark
 260 could be in either of the two beams. You could therefore calculate the
 261 result once only and store it as a private data member of the class.
 262 It is optional whether you want to use this method, however, or put
 263 everything in <code>sigmaHat</code>.
 264
 265
 266 <a name="method3"></a>
 267 <p/><strong>double SigmaProcess::sigmaHat() &nbsp;</strong> <br/>
 268 is the key method for cross section calculations and returns a cross section
 269 value, as described in the previous section. It is called when also a
 270 preliminary set of incoming flavours has been picked, in addition to the
 271 kinematical ones already available for <code>sigmaKin</code>.
 272 Typically <code>sigmaHat</code> is called inside a loop over all allowed
 273 incoming flavour combinations, stored in <code>id1</code> and
 274 <code>id2</code>, with fixed kinematics, as already illustrated above.
 275 The sum over the different flavour combinations provides the total
 276 cross section, while their relative size is used to make a selection of
 277 a specific incoming state.
 278
 279
 280 <a name="method4"></a>
 281 <p/><strong>bool SigmaProcess::setupForME() &nbsp;</strong> <br/>
 282 to be called by the user from inside <code>sigmaHat()</code>
 283 (or possibly <code>sigmaKin()</code>) to setup alternative kinematics
 284 in the <code>mME[5]</code> and <code>pME[5]</code> arrays, better
 285 suited for matrix-element calculations. See the end of the previous
 286 section for a more detailed description. Should the method return
 287 <code>false</code> then the conversion did not work, and
 288 <code>sigmaHat()</code> (or <code>sigmaKin()</code>) should be set to
 289 vanish.
 290
 291
 292 <a name="method5"></a>
 293 <p/><strong>void SigmaProcess::setIdColAcol() &nbsp;</strong> <br/>
 294 is called only once an initial state and a kinematical configuration has
 295 been picked. This routine must set the complete flavour information and
 296 the colour flow of the process. This may involve further random choices,
 297 between different possible final-state flavours or between possible
 298 competing colour flows. Private data members of the class may be used to
 299 retain some information from the previous steps above.
 300 <br/>When this routine is called the two incoming flavours have already
 301 been selected and are available in <code>id1</code> and <code>id2</code>,
 302 whereas the one, two or three outgoing ones either are fixed for a given
 303 process or can be determined from the instate (e.g. whether a <i>W^+</i>
 304 or <i>W^-</i> was produced).  There is also a standard method in
 305 <code>VCKM</code> to pick a final flavour from an initial one with CKM
 306 mixing. Once you have figured out the value of
 307 <code>id3</code> and, the case being, <code>id4</code> and
 308 <code>id5</code>, you store these values permanently by a call
 309 <code>setId( id1, id2, id3, id4, id5)</code>, where the last two may be
 310 omitted if irrelevant.
 311 <br/>Correspondingly, the colours are stored with
 312 <code>setColAcol( col1, acol1, col2, acol2, col3, acol3, col4, acol4,
 313 col5, acol5)</code>, where the final ones may be omitted if irrelevant.
 314 Les Houches style colour tags are used, but starting with number 1
 315 (and later shifted by the currently requested offset). The
 316 input is grouped particle by particle, with the colour index before the
 317 anticolour one. You may need to select colour flow dynamically, depending
 318 on the kinematics, when several distinct possibilities exist. Trivial
 319 operations, like swapping colours and anticolours, can be done with
 320 existing methods.
 321 <br/>When the <code>id3Mass()</code> and <code>id4Mass()</code>
 322 methods have been used, the order of the outgoing particles may be
 323 inconsistent with the way the <i>tHat</i> and <i>uHat</i>
 324 variables have been defined. A typical example would be a process like
 325 <i>q g -> q' W</i> with <i>tHat</i> defined between incoming and
 326 outgoing quark, but where <code>id3Mass() = 24</code> and so the
 327 process is to be stored as <i>q g -> W q'</i>. One should then put
 328 the variable <code>swapTU = true</code> in <code>setIdColAcol()</code>
 329 for each event where the <i>tHat</i> and <i>uHat</i> variables
 330 should be swapped before the event kinematics is reconstructed. This
 331 variable is automatically restored to <code>false</code> for each new
 332 event.
 333
 334
 335 <a name="method6"></a>
 336 <p/><strong>double SigmaProcess::weightDecayFlav( Event& process) &nbsp;</strong> <br/>
 337 is called to allow a reweighting of the simultaneous flavour choices of
 338 resonance decay products. Is currently only used for the
 339 <i>q qbar -> gamma*/Z^0 gamma*/Z^0</i> process, and will likely not
 340 be of interest for you.
 341
 342
 343 <a name="method7"></a>
 344 <p/><strong>double SigmaProcess::weightDecay( Event& process, int iResBeg, int iResEnd) &nbsp;</strong> <br/>
 345 is called when the basic process has one or several resonances, after each
 346 set of related resonances in <code>process[i]</code>,
 347 <code>iResBeg</code> &lt;= <code>i </code> &lt;= <code>iResEnd</code>,
 348 has been allowed to decay. The calculated weight, to be normalized
 349 to the range between 0 and 1, is used to decide whether to accept the
 350 decay(s) or try for a new decay configuration. The base-class version of
 351 this method returns unity, i.e. gives isotropic decays by default.
 352 This method may be called repeatedly for a single event. For instance, in
 353 <i>q qbar -> H^0 Z^0</i> with <i>H^0 -> W^+ W^-</i>, a first call
 354 would be made after the <i>H^0</i> and <i>Z^0</i> decays, and then
 355 depend only on the <i>Z^0</i> decay angles since the <i>H^0</i>
 356 decays isotropically. The second call would be after the <i>W^+ W^-</i>
 357 decays and then involve correlations between the four daughter fermions.
 358
 359
 360 <a name="method8"></a>
 361 <p/><strong>string SigmaProcess::name() &nbsp;</strong> <br/>
 362 returns the name of the process, as you want it to be shown in listings.
 363
 364
 365 <a name="method9"></a>
 366 <p/><strong>int SigmaProcess::code() &nbsp;</strong> <br/>
 367 returns an integer identifier of the process. This has no internal function,
 368 but is only intended as a service for the user to rapidly (and hopefully
 369 uniquely) identify which process occurred in a given event. Numbers below
 370 10000 are reserved for internal PYTHIA use.
 371
 372
 373 <a name="method10"></a>
 374 <p/><strong>string SigmaProcess::inFlux() &nbsp;</strong> <br/>
 375 this string specifies the combinations of incoming partons that are
 376 allowed for the process under consideration, and thereby which incoming
 377 flavours <code>id1</code> and <code>id2</code> the <code>sigmaHat()</code>
 378 calls will be looped over. It is always possible to pick a wider flavour
 379 selection than strictly required and then put to zero cross sections in
 380 the superfluous channels, but of course this may cost some extra execution
 381 time. Currently allowed options are:
 382 <br/>* <code>gg</code>: two gluons.
 383 <br/>* <code>qg</code>: one (anti)quark and one gluon.
 384 <br/>* <code>qq</code>: any combination of two quarks, two antiquarks or
 385 a quark and an antiquark.
 386 <br/>* <code>qqbarSame</code>: a quark and its antiquark;
 387 this is a subset of the above <code>qq</code> option.
 388 <br/>* <code>ff</code>: any combination of two fermions, two antifermions
 389 or a fermion and an antifermion; is the same as <code>qq</code> for
 390 hadron beams but also allows processes to work with lepton beams.
 391 <br/>* <code>ffbarSame</code>: a fermion and its antifermion; is the
 392 same as <code>qqbarSame</code> for hadron beams but also allows processes
 393 to work with lepton beams.
 394 <br/>* <code>ffbarChg</code>: a fermion and an antifermion that combine
 395 to give charge +-1.
 396 <br/>* <code>fgm</code>: a fermion and a photon (gamma).
 397 <br/>* <code>ggm</code>: a gluon and a photon.
 398 <br/>* <code>gmgm</code>: two photons.
 399
 400
 401 <a name="method11"></a>
 402 <p/><strong>bool SigmaProcess::convert2mb() &nbsp;</strong> <br/>
 403 it is assumed that cross sections normally come in dimensions such that
 404 they, when integrated over the relevant phase space, obtain the dimension
 405 GeV^-2, and therefore need to be converted to mb. If the cross section
 406 is already encoded as mb then <code>convert2mb()</code> should be
 407 overloaded to instead return <code>false</code>.
 408
 409
 410 <a name="method12"></a>
 411 <p/><strong>bool SigmaProcess::convertM2() &nbsp;</strong> <br/>
 412 it is assumed that <i>2 -> 1</i> cross sections are encoded as
 413 <i>sigmaHat(sHat)</i>, and <i>2 -> 2</i> ones as
 414 <i>d(sigmaHat)/d(tHat)</i> in the <code>SigmaProcess::sigmaHat()</code>
 415 methods. If <code>convertM2()</code> is overloaded to instead return
 416 <code>true</code> then the return value is instead assumed to be the
 417 squared matrix element <i>|M|^2</i>, and
 418 <code>SigmaProcess::sigmaHatWrap(...)</code> converts to
 419 <i>sigmaHat(sHat)</i> or <i>d(sigmaHat)/d(tHat)</i>, respectively.
 420 This switch has no effect on <i>2 -> 3</i> processes, where
 421 <i>|M|^2</i> is the only allowed input anyway.
 422
 423
 424 <a name="method13"></a>
 425 <p/><strong>int SigmaProcess::id3Mass() &nbsp;</strong> <br/>
 426
 427 <strong>int SigmaProcess::id4Mass() &nbsp;</strong> <br/>
 428
 429 <strong>int SigmaProcess::id5Mass() &nbsp;</strong> <br/>
 430 are the one, two or three final-state flavours, where masses are to be
 431 selected before the matrix elements are evaluated. Only the absolute value
 432 should be given. For massless particles, like gluons and photons, one need
 433 not give anything, i.e. one defaults to 0. The same goes for normal light
 434 quarks, where masses presumably are not implemented in the matrix elements.
 435 Later on, these quarks can still (automatically) obtain constituent masses,
 436 once a <i>u</i>, <i>d</i> or <i>s</i> flavour has been selected.
 437
 438
 439 <a name="method14"></a>
 440 <p/><strong>int SigmaProcess::resonanceA() &nbsp;</strong> <br/>
 441
 442 <strong>int SigmaProcess::resonanceB() &nbsp;</strong> <br/>
 443 are the codes of up to two <i>s</i>-channel resonances contributing to
 444 the matrix elements. These are used by the program to improve the phase-space
 445 selection efficiency, by partly sampling according to the relevant
 446 Breit-Wigner distributions. Massless resonances (the gluon and photon) need not be
 447 specified.
 448
 449
 450 <a name="method15"></a>
 451 <p/><strong>bool SigmaProcess::isSChannel() &nbsp;</strong> <br/>
 452 normally the choice of renormalization and factorization scales in
 453 <i>2 -> 2</i> and <i>2 -> 3</i> processes is based on the assumption
 454 that <i>t</i>- and <i>u</i>-channel exchanges dominates the
 455 cross section. In cases such as <i>f fbar -> gamma* -> f' fbar'</i> a
 456 <i>2 -> 2</i> process actually ought to be given scales as a
 457 <i>2 -> 1</i> one, in the sense that it proceeds entirely through
 458 an <i>s</i>-channel resonance. This can be achieved if you override the
 459 default <code>false</code> to return <code>true</code>. See further the
 460 page on <a href="CouplingsAndScales.html" target="page">couplings and scales</a>.
 461
 462
 463 <a name="method16"></a>
 464 <p/><strong>int SigmaProcess::idSChannel() &nbsp;</strong> <br/>
 465 normally no intermediate state is shown in the event record for
 466 <i>2 -> 2</i> and <i>2 -> 3</i> processes. However, in case
 467 that <code>idSChannel</code> is overloaded to return a nonzero value,
 468 an intermediate particle with that identity code is inserted into the
 469 event record, to make it a <i>2 -> 1 -> 2</i> or <i>2 -> 1 -> 3</i>
 470 process. Thus if both <code>isSChannel</code> and <code>idSChannel</code>
 471 are overloaded, a process will behave and look like it proceeded through
 472 a resonance. The one difference is that the implementation of the
 473 matrix element is not based on the division into a production and a
 474 decay of an intermediate resonance, but is directly describing the
 475 transition from the initial to the final state.
 476
 477
 478 <a name="method17"></a>
 479 <p/><strong>int SigmaProcess::isQCD3body() &nbsp;</strong> <br/>
 480 there are two different 3-body phase-space selection machineries,
 481 of which the non-QCD one is default. If you overload this method
 482 instead the QCD-inspired machinery will be used. The differences
 483 between these two is related to which
 484 <a href="PhaseSpaceCuts.html" target="page">phase space cuts</a>
 485 can be set, and also that the QCD machinery assumes (almost) massless
 486 outgoing partons.
 487
 488
 489 <a name="method18"></a>
 490 <p/><strong>int SigmaProcess::idTchan1() &nbsp;</strong> <br/>
 491
 492 <strong>int SigmaProcess::idTchan2() &nbsp;</strong> <br/>
 493 the non-QCD <i>2 -> 3</i> phase space selection machinery is rather
 494 primitive, as already mentioned. The efficiency can be improved in
 495 processes that proceed though <i>t</i>-channel exchanges, such as
 496 <i>q qbar' -> H^0 q qbar'</i> via <i>Z^0 Z^0</i> fusion, if the identity
 497 of the  <i>t</i>-channel-exchanged particles on the two side of the
 498 event are provided. Only the absolute value is of interest.
 499
 500
 501 <a name="method19"></a>
 502 <p/><strong>double SigmaProcess::tChanFracPow1() &nbsp;</strong> <br/>
 503
 504 <strong>double SigmaProcess::tChanFracPow2() &nbsp;</strong> <br/>
 505 in the above kind of <i>2 -> 3</i> phase-space selection, the
 506 sampling of <i>pT^2</i> is done with one part flat, one part weighted
 507 like <i>1 / (pT^2 + m_R^2)</i> and one part  like
 508 <i>1 / (pT^2 + m_R^2)^2</i>. The above values provide the relative
 509 amount put in the latter two channels, respectively, with the first
 510 obtaining the rest. Thus the sum of <code>tChanFracPow1()</code> and
 511 <code>tChanFracPow2()</code> must be below unity. The final results
 512 should be independent of these numbers, but the Monte Carlo efficiency
 513 may be quite low for a bad choice. Here <i>m_R</i> is the mass of the
 514 exchanged resonance specified by <code>idTchan1()</code> or
 515 <code>idTchan2()</code>. Note that the order of the final-state
 516 listing is important in the above <i>q qbar' -> H^0 q qbar'</i> example,
 517 i.e. the <i>H^0</i> must be returned by <code>id3Mass()</code>,
 518 since it is actually the <i>pT^2</i> of the latter two that are
 519 selected independently, with the first <i>pT</i> then fixed
 520 by transverse-momentum conservation.
 521
 522
 523 <a name="method20"></a>
 524 <p/><strong>bool SigmaProcess::useMirrorWeight() &nbsp;</strong> <br/>
 525 in <i>2 -> 3</i> processes the phase space selection used here
 526 involves a twofold ambiguity basically corresponding to a flipping of
 527 the positions of last two outgoing particles. These are assumed equally
 528 likely by default, <code>false</code>, but for processes proceeding entirely
 529 through <i>t</i>-channel exchange the Monte Carlo efficiency can be
 530 improved by making a preselection based on the relative propagator
 531 weights, <code>true</code>.
 532
 533
 534 <a name="method21"></a>
 535 <p/><strong>int SigmaProcess::gmZmode() &nbsp;</strong> <br/>
 536 allows a possibility to override the global mode
 537 <code><a href="ElectroweakProcesses.html" target="page">WeakZ0:gmZmode</a></code>
 538 for a specific process. The global mode normally is used to switch off
 539 parts of the <i>gamma^*/Z^0</i> propagator for test purposes. The
 540 above local mode is useful for processes where a <i>Z^0</i> really is
 541 that and nothing more, such as <i>q qbar -> H^0 Z^0</i>. The default
 542 value -1 returned by <code>gmZmode()</code> ensures that the global
 543 mode is used, while 0 gives full <i>gamma^*/Z^0</i> interference,
 544 1 <i>gamma^*</i> only and 2 <i>Z^0</i> only.
 545
 546
 547 <h3>Access to a process</h3>
 548
 549 Once you have implemented a class, it is straightforward to make use of
 550 it in a run. Assume you have written a new class <code>MySigma</code>,
 551 which inherits from <code>Sigma1Process</code>, <code>Sigma2Process</code>
 552 or <code>Sigma3Process</code>, which in their turn inherit from
 553 <code>SigmaProcess</code>. You then create an instance of this class
 554 and hand it in to a <code>pythia</code> object with
 555 <pre>
 556       SigmaProcess* mySigma = new MySigma();
 557       pythia.setSigmaPtr( mySigma);
 558 </pre>
 559 If you have several processes you can repeat the procedure any number
 560 of times. When <code>pythia.init(...)</code> is called these processes
 561 are initialized along with any internal processes you may have switched on,
 562 and treated in exactly the same manner. The  <code>pythia.next()</code>
 563 will therefore generate a mix of the different kinds of processes without
 564 distinction. See also the <a href="ProgramFlow.html" target="page">Program Flow</a>
 565 description.
 566
 567 <p/>
 568 If the code should be of good quality and general usefulness, it would
 569 be simple to include it as a permanently available process in the
 570 standard program distribution. The final step of that integration ought to
 571 be left for the PYTHIA authors, but here is a description of what is
 572 required.
 573
 574 <p/>
 575 A flag has to be defined, that allows the process to be switched on;
 576 by default it should always be off. The name of the flag should be
 577 chosen of the type <code>model:process</code>. Here the
 578 <code>model</code> would be related to the general scenario considered,
 579 e.g. <code>Compositeness</code>, while <code>process</code> would
 580 specify instate and outstate, separated by a 2 (= to), e.g.
 581 <code>ug2u*g</code>.
 582 When several processes are implemented and "belong together" it is
 583 also useful to define a <code>model:all</code> switch that affects
 584 all the separate processes.
 585
 586 <p/>
 587 The flags should normally be stored in the <code>ProcessSelection.xml</code>
 588 file or one of its daughters for a specific kind of processes. This is to
 589 make them easily found by users. You could create and use your own
 590 <code>.xml</code> file, so long as you then add that name to the
 591 list of files in the <code>Index.xml</code> file. (If not,
 592 the flags would never be created and the program would not work.)
 593
 594 <p/>
 595 In the <code>ProcessContainer.c</code> file, the
 596 <code>SetupContainers::init()</code> method needs to be expanded to
 597 create instances of the processes switched on. This code is fairly
 598 repetitive, and should be easy to copy and modify from the code
 599 already there. The basic structure is
 600 <br/>(i) check whether a process is requested by the user and, if so,
 601 <br/>(ii) create an instance of the matrix-element class,
 602 <br/>(iii)create a container for the matrix element and its associated
 603 phase-space handling, and
 604 <br>(iv) add the container to the existing process list.
 605
 606 <p/>
 607 Two minor variations are possible. One is that a set of related
 608 processes are lumped inside the the same initial check, i.e. are
 609 switched on all together. The second is that the matrix-element
 610 constructor may take arguments, as specified by you (see above).
 611 If so, the same basic matrix element may be recycled for a set of
 612 related processes, e.g. one for a composite <i>u</i> and one for
 613 a composite <i>d</i>. Obviously these variations may be combined.
 614
 615 </body>
 616 </html>
 617
 618 <!-- Copyright (C) 2013 Torbjorn Sjostrand -->