PYTHIA8/pythia8145/examples/main21.cc

   1 // main21.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // This is a simple test program.
   7 // It illustrates how to feed in a single particle
   8 // or a toy parton-level configurations.
   9
  10 #include "Pythia.h"
  11 using namespace Pythia8;
  12
  13 //==========================================================================
  14
  15 // Single-particle gun. The particle must be a colour singlet.
  16 // Input: flavour, energy, direction (theta, phi).
  17 // If theta < 0 then random choice over solid angle.
  18
  19 void fillParticle(int id, double ee, double thetaIn, double phiIn,
  20                   Event& event, ParticleData& pdt, Rndm& rndm) {
  21
  22   // Reset event record to allow for new event.
  23   event.reset();
  24
  25   // Select particle mass; where relevant according to Breit-Wigner.
  26   double mm = pdt.mass(id);
  27   double pp = sqrtpos(ee*ee - mm*mm);
  28
  29   // Angles as input or uniform in solid angle.
  30   double cThe, sThe, phi;
  31   if (thetaIn >= 0.) {
  32     cThe = cos(thetaIn);
  33     sThe = sin(thetaIn);
  34     phi  = phiIn;
  35   } else {
  36     cThe = 2. * rndm.flat() - 1.;
  37     sThe = sqrtpos(1. - cThe * cThe);
  38     phi = 2. * M_PI * rndm.flat();
  39   }
  40
  41   // Store the particle in the event record.
  42   event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi),
  43     pp * cThe, ee, mm);
  44 }
  45
  46 //==========================================================================
  47
  48 // Simple method to do the filling of partons into the event record.
  49
  50 void fillPartons(int type, double ee, Event& event, ParticleData& pdt,
  51   Rndm& rndm) {
  52
  53   // Reset event record to allow for new event.
  54   event.reset();
  55
  56   // Information on a q qbar system, to be hadronized.
  57   if (type == 1) {
  58     int    id = 2;
  59     double mm = pdt.m0(id);
  60     double pp = sqrtpos(ee*ee - mm*mm);
  61     event.append(  id, 23, 101,   0, 0., 0.,  pp, ee, mm);
  62     event.append( -id, 23,   0, 101, 0., 0., -pp, ee, mm);
  63
  64   // Information on a g g system, to be hadronized.
  65   } else if (type == 2) {
  66     event.append( 21, 23, 101, 102, 0., 0.,  ee, ee);
  67     event.append( 21, 23, 102, 101, 0., 0., -ee, ee);
  68
  69   // Information on a g g g system, to be hadronized.
  70   } else if (type == 3) {
  71     event.append( 21, 23, 101, 102,        0., 0.,        ee, ee);
  72     event.append( 21, 23, 102, 103,  0.8 * ee, 0., -0.6 * ee, ee);
  73     event.append( 21, 23, 103, 101, -0.8 * ee, 0., -0.6 * ee, ee);
  74
  75   // Information on a q q q junction system, to be hadronized.
  76   } else if (type == 4 || type == 5) {
  77
  78     // Need a colour singlet mother parton to define junction origin.
  79     event.append( 1000022, -21, 0, 0, 2, 4, 0, 0,
  80                   0., 0., 1.01 * ee, 1.01 * ee);
  81
  82     // The three endpoint q q q; the minimal system.
  83     double rt75 = sqrt(0.75);
  84     event.append( 2, 23, 1, 0, 0, 0, 101, 0,
  85                           0., 0., 1.01 * ee, 1.01 * ee);
  86     event.append( 2, 23, 1, 0, 0, 0, 102, 0,
  87                    rt75 * ee, 0., -0.5 * ee,        ee );
  88     event.append( 1, 23, 1, 0, 0, 0, 103, 0,
  89                   -rt75 * ee, 0., -0.5 * ee,        ee );
  90
  91     // Define the qqq configuration as starting point for adding gluons.
  92     if (type == 5) {
  93       int colNow[4] = {0, 101, 102, 103};
  94       Vec4 pQ[4];
  95       pQ[1] = Vec4(0., 0., 1., 0.);
  96       pQ[2] = Vec4( rt75, 0., -0.5, 0.);
  97       pQ[3] = Vec4(-rt75, 0., -0.5, 0.);
  98
  99       // Minimal cos(q-g opening angle), allows more or less nasty events.
 100       double cosThetaMin =0.;
 101
 102       // Add a few gluons (almost) at random.
 103       for (int nglu = 0; nglu < 5; ++nglu) {
 104         int iq = 1 + int( 2.99999 * rndm.flat() );
 105         double px, py, pz, e, prod;
 106         do {
 107           e =  ee * rndm.flat();
 108           double cThe = 2. * rndm.flat() - 1.;
 109           double phi = 2. * M_PI * rndm.flat();
 110           px = e * sqrt(1. - cThe*cThe) * cos(phi);
 111           py = e * sqrt(1. - cThe*cThe) * sin(phi);
 112           pz = e * cThe;
 113           prod = ( px * pQ[iq].px() + py * pQ[iq].py() + pz * pQ[iq].pz() )
 114             / e;
 115         } while (prod < cosThetaMin);
 116         int colNew = 104 + nglu;
 117         event.append( 21, 23, 1, 0, 0, 0, colNew, colNow[iq],
 118           px, py, pz, e, 0.);
 119         colNow[iq] = colNew;
 120       }
 121       // Update daughter range of mother.
 122       event[1].daughters(1, event.size() - 1);
 123
 124     }
 125
 126   // Information on a q q qbar qbar dijunction system, to be hadronized.
 127   } else if (type >= 6) {
 128
 129     // The two fictitious beam remnant particles; needed for junctions.
 130     event.append( 2212, -12, 0, 0, 3, 5, 0, 0, 0., 0., ee, ee, 0.);
 131     event.append(-2212, -12, 0, 0, 6, 8, 0, 0, 0., 0., ee, ee, 0.);
 132
 133     // Opening angle between "diquark" legs.
 134     double theta = 0.2;
 135     double cThe = cos(theta);
 136     double sThe = sin(theta);
 137
 138     // Set one colour depending on whether more gluons or not.
 139     int acol = (type == 6) ? 103 : 106;
 140
 141     // The four endpoint q q qbar qbar; the minimal system.
 142     // Two additional fictitious partons to make up original beams.
 143     event.append(  2,   23, 1, 0, 0, 0, 101, 0,
 144                   ee * sThe, 0.,  ee * cThe, ee, 0.);
 145     event.append(  1,   23, 1, 0, 0, 0, 102, 0,
 146                  -ee * sThe, 0.,  ee * cThe, ee, 0.);
 147     event.append(  2, -21, 1, 0, 0, 0, 103, 0,
 148                          0., 0.,  ee       , ee, 0.);
 149     event.append( -2,   23, 2, 0, 0, 0, 0, 104,
 150                   ee * sThe, 0., -ee * cThe, ee, 0.);
 151     event.append( -1,   23, 2, 0, 0, 0, 0, 105,
 152                  -ee * sThe, 0., -ee * cThe, ee, 0.);
 153     event.append( -2, -21, 2, 0, 0, 0, 0, acol,
 154                          0., 0., -ee       , ee, 0.);
 155
 156     // Add extra gluons on string between junctions.
 157     if (type == 7) {
 158       event.append( 21, 23, 5, 8, 0, 0, 103, 106, 0., ee, 0., ee, 0.);
 159     } else if (type == 8) {
 160       event.append( 21, 23, 5, 8, 0, 0, 103, 108, 0., ee, 0., ee, 0.);
 161       event.append( 21, 23, 5, 8, 0, 0, 108, 106, 0.,-ee, 0., ee, 0.);
 162     } else if (type == 9) {
 163       event.append( 21, 23, 5, 8, 0, 0, 103, 107, 0., ee, 0., ee, 0.);
 164       event.append( 21, 23, 5, 8, 0, 0, 107, 108, ee, 0., 0., ee, 0.);
 165       event.append( 21, 23, 5, 8, 0, 0, 108, 106, 0.,-ee, 0., ee, 0.);
 166     } else if (type == 10) {
 167       event.append( 21, 23, 5, 8, 0, 0, 103, 107, 0., ee, 0., ee, 0.);
 168       event.append( 21, 23, 5, 8, 0, 0, 107, 108, ee, 0., 0., ee, 0.);
 169       event.append( 21, 23, 5, 8, 0, 0, 108, 109, 0.,-ee, 0., ee, 0.);
 170       event.append( 21, 23, 5, 8, 0, 0, 109, 106,-ee, 0., 0., ee, 0.);
 171     }
 172
 173   // No more cases: done.
 174   }
 175 }
 176
 177 //==========================================================================
 178
 179 int main() {
 180
 181   // Pick kind of events to generate:
 182   // 0 = single-particle gun.
 183   // 1 = q qbar.
 184   // 2 = g g.
 185   // 3 = g g g.
 186   // 4 = minimal q q q junction topology.
 187   // 5 = q q q junction topology with gluons on the strings.
 188   // 6 = q q qbar qbar dijunction topology, no gluons.
 189   // 7 - 10 : ditto, but with 1 - 4 gluons on string between junctions.
 190   int type = 0;
 191
 192   // Set particle species and energy for single-particle gun.
 193   int    idGun = 15;
 194   double eeGun = 20.;
 195
 196   // Set typical energy per parton.
 197   double ee = 20.0;
 198
 199   // Set number of events to generate and to list.
 200   int nEvent = 10000;
 201   int nList = 3;
 202
 203   // Generator; shorthand for event and particleData.
 204   Pythia pythia;
 205   Event& event      = pythia.event;
 206   ParticleData& pdt = pythia.particleData;
 207
 208   // Key requirement: switch off ProcessLevel, and thereby also PartonLevel.
 209   pythia.readString("ProcessLevel:all = off");
 210
 211   // Optionally switch off decays.
 212   //pythia.readString("HadronLevel:Decay = off");
 213
 214   // Provide printout of initial information.
 215   pythia.settings.listChanged();
 216
 217   // Initialize.
 218   pythia.init();
 219
 220   // Book histograms.
 221   Hist epCons("deviation from energy-momentum conservation",100,0.,1e-4);
 222   Hist chgCons("deviation from charge conservation",57,-9.5,9.5);
 223   Hist nFinal("final particle multiplicity",100,-0.5,99.5);
 224   Hist dnparticledp("dn/dp for particles",100,0.,ee);
 225   Hist status85("multiplicity status code 85",50,-0.5,49.5);
 226   Hist status86("multiplicity status code 86",50,-0.5,49.5);
 227   Hist status83("multiplicity status code 83",50,-0.5,49.5);
 228   Hist status84("multiplicity status code 84",50,-0.5,49.5);
 229   Hist dndtheta("particle flow in event plane",100,-M_PI,M_PI);
 230   Hist dedtheta("energy flow in event plane",100,-M_PI,M_PI);
 231   Hist dpartondtheta("parton flow in event plane",100,-M_PI,M_PI);
 232   Hist dndyAnti("dn/dy primaries antijunction",100, -10., 10.);
 233   Hist dndyJun("dn/dy primaries junction",100, -10., 10.);
 234   Hist dndySum("dn/dy all primaries",100, -10., 10.);
 235
 236   // Begin of event loop.
 237   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
 238     if (iEvent%(max(1,nEvent/20)) == 0) cout << " Now begin event "
 239       << iEvent << endl;
 240
 241     // Set up single particle, with random direction in solid angle.
 242     if (type == 0) fillParticle( idGun, eeGun, -1., 0., event, pdt, pythia.rndm);
 243
 244     // Set up parton-level configuration.
 245     else fillPartons( type, ee, event, pdt, pythia.rndm);
 246
 247     // Generate events. Quit if failure.
 248     if (!pythia.next()) {
 249       cout << " Event generation aborted prematurely, owing to error!\n";
 250       break;
 251     }
 252
 253     // List first few events.
 254     if (iEvent < nList) {
 255       event.list();
 256       // Also list junctions.
 257       event.listJunctions();
 258     }
 259
 260     // Initialize statistics.
 261     Vec4 pSum = - event[0].p();
 262     double chargeSum = 0.;
 263     if (type == 0) chargeSum = -event[1].charge();
 264     if (type == 4 || type == 5) chargeSum = -1;
 265     int nFin = 0;
 266     int n85 = 0;
 267     int n86 = 0;
 268     int n83 = 0;
 269     int n84 = 0;
 270
 271     // Loop over all particles.
 272     for (int i = 0; i < event.size(); ++i) {
 273       int status = event[i].statusAbs();
 274
 275       // Find any unrecognized particle codes.
 276       int id = event[i].id();
 277       if (id == 0 || !pdt.isParticle(id))
 278         cout << " Error! Unknown code id = " << id << "\n";
 279
 280       // Find particles with E-p mismatch.
 281       double eCalc = event[i].eCalc();
 282       if (abs(eCalc/event[i].e() - 1.) > 1e-6) cout << " e mismatch, i = "
 283         << i << " e_nominal = " << event[i].e() << " e-from-p = "
 284         << eCalc << " m-from-e " << event[i].mCalc() << "\n";
 285
 286       // Parton flow in event plane.
 287       if (status == 71 || status == 72) {
 288         double thetaXZ = event[i].thetaXZ();
 289         dpartondtheta.fill(thetaXZ);
 290       }
 291
 292       // Origin of primary hadrons.
 293       if (status == 85) ++n85;
 294       if (status == 86) ++n86;
 295       if (status == 83) ++n83;
 296       if (status == 84) ++n84;
 297
 298       // Flow of primary hadrons in the event plane.
 299       if (status > 80 && status < 90) {
 300         double eAbs = event[i].e();
 301         if (eAbs < 0.) {cout << " e < 0 line " << i; event.list();}
 302         double thetaXZ = event[i].thetaXZ();
 303         dndtheta.fill(thetaXZ);
 304         dedtheta.fill(thetaXZ, eAbs);
 305
 306         // Rapidity distribution of primary hadrons.
 307         double y = event[i].y();
 308         dndySum.fill(y);
 309         if (type >= 6) {
 310           int motherId = event[event[i].mother1()].id();
 311           if (motherId > 0 ) dndyJun.fill(event[i].y());
 312           else dndyAnti.fill(event[i].y());
 313         }
 314       }
 315
 316       // Study final-state particles.
 317       if (event[i].isFinal()) {
 318         pSum += event[i].p();
 319         chargeSum += event[i].charge();
 320         nFin++;
 321         double pAbs = event[i].pAbs();
 322         dnparticledp.fill(pAbs);
 323       }
 324     }
 325
 326     // Fill histograms once for each event.
 327     double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
 328       + abs(pSum.pz());
 329     epCons.fill(epDev);
 330     chgCons.fill(chargeSum);
 331     nFinal.fill(nFin);
 332     status85.fill(n85);
 333     status86.fill(n86);
 334     status83.fill(n83);
 335     status84.fill(n84);
 336     if (epDev > 1e-3  || abs(chargeSum) > 0.1) event.list();
 337
 338   // End of event loop.
 339   }
 340
 341   // Print statistics, histograms and done.
 342   pythia.statistics();
 343   cout << epCons << chgCons << nFinal << dnparticledp
 344        << dndtheta << dedtheta << dpartondtheta << dndySum;
 345   if (type >= 4) cout << status85 << status86 << status83
 346        << status84;
 347   if (type >= 6) cout << dndyJun << dndyAnti;
 348
 349   // Done.
 350   return 0;
 351 }