PYTHIA8/pythia8145/examples/main19.cc

   1 // main19.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 program runs four instances of Pythia simultaneously,
   7 // one for signal events, one for pileup background ones, and two
   8 // For beam-gas background ones. Note that Pythia does not do nuclear
   9 // effects, so beam-gas is represented by "fixed-target" pp collisions.
  10 // The = and += overloaded operators are used to join several
  11 // event records into one, but should be used with caution.
  12
  13 // Note that each instance of Pythia is independent of any other,
  14 // but with two important points to remember.
  15 // 1) By default all generate the same random number sequence,
  16 //    which has to be corrected if they are to generate the same
  17 //    physics, like the two beam-gas ones below.
  18 // 2) Interfaces to external Fortran programs are "by definition" static.
  19 //    Thus it is not a good idea to use LHAPDF to set different PDF's
  20 //    in different instances.
  21
  22 #include "Pythia.h"
  23 using namespace Pythia8;
  24
  25 //==========================================================================
  26
  27 // Method to pick a number according to a Poissonian distribution.
  28
  29 int poisson(double nAvg, Rndm& rndm) {
  30
  31   // Set maximum to avoid overflow.
  32   const int NMAX = 100;
  33
  34   // Random number.
  35   double rPoisson = rndm.flat() * exp(nAvg);
  36
  37   // Initialize.
  38   double rSum  = 0.;
  39   double rTerm = 1.;
  40
  41   // Add to sum and check whether done.
  42   for (int i = 0; i < NMAX; ) {
  43     rSum += rTerm;
  44     if (rSum > rPoisson) return i;
  45
  46     // Evaluate next term.
  47     ++i;
  48     rTerm *= nAvg / i;
  49   }
  50
  51   // Emergency return.
  52   return NMAX;
  53 }
  54
  55 //==========================================================================
  56
  57 int main() {
  58
  59   // Number of signal events to generate.
  60   int nEvent = 100;
  61
  62   // Beam Energy.
  63   double eBeam = 7000.;
  64
  65   // Average number of pileup events per signal event.
  66   double nPileupAvg = 2.5;
  67
  68   // Average number of beam-gas events per signal ones, on two sides.
  69   double nBeamAGasAvg = 0.5;
  70   double nBeamBGasAvg = 0.5;
  71
  72   // Four generator instances.
  73   Pythia pythiaSignal;
  74   Pythia pythiaPileup;
  75   Pythia pythiaBeamAGas;
  76   Pythia pythiaBeamBGas;
  77
  78   // One object where all individual events are to be collected.
  79   Event sumEvent;
  80
  81   // Initialize generator for signal processes.
  82   pythiaSignal.readString("HardQCD:all = on");
  83   pythiaSignal.readString("PhaseSpace:pTHatMin = 50.");
  84   pythiaSignal.init( 2212, 2212, 2. * eBeam);
  85
  86   // Initialize generator for pileup (background) processes.
  87   pythiaPileup.readString("Random:setSeed = on");
  88   pythiaPileup.readString("Random:seed = 10000002");
  89   pythiaPileup.readString("SoftQCD:all = on");
  90   pythiaPileup.init( 2212, 2212, 2. * eBeam);
  91
  92   // Initialize generators for beam A - gas (background) processes.
  93   pythiaBeamAGas.readString("Random:setSeed = on");
  94   pythiaBeamAGas.readString("Random:seed = 10000003");
  95   pythiaBeamAGas.readString("SoftQCD:all = on");
  96   pythiaBeamAGas.init( 2212, 2212, eBeam, 0.);
  97
  98   // Initialize generators for beam B - gas (background) processes.
  99   pythiaBeamBGas.readString("Random:setSeed = on");
 100   pythiaBeamBGas.readString("Random:seed = 10000004");
 101   pythiaBeamBGas.readString("SoftQCD:all = on");
 102   pythiaBeamBGas.init( 2212, 2212, 0., eBeam);
 103
 104   // Histograms: number of pileups, total charged multiplicity.
 105   Hist nPileH("number of pileup events per signal event", 100, -0.5, 99.5);
 106   Hist nAGH("number of beam A + gas events per signal event", 100, -0.5, 99.5);
 107   Hist nBGH("number of beam B + gas events per signal event", 100, -0.5, 99.5);
 108   Hist nChgH("number of charged multiplicity",100, -0.5, 1999.5);
 109   Hist sumPZH("total pZ of system",100, -100000., 100000.);
 110
 111   // Loop over events.
 112   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
 113
 114     // Generate a signal event. Copy this event into sumEvent.
 115     if (!pythiaSignal.next()) continue;
 116     sumEvent = pythiaSignal.event;
 117
 118     // Select the number of pileup events to generate.
 119     int nPileup = poisson(nPileupAvg, pythiaPileup.rndm);
 120     nPileH.fill( nPileup );
 121
 122     // Generate a number of pileup events. Add them to sumEvent.
 123     for (int iPileup = 0; iPileup < nPileup; ++iPileup) {
 124       pythiaPileup.next();
 125       sumEvent += pythiaPileup.event;
 126     }
 127
 128     // Select the number of beam A + gas events to generate.
 129     int nBeamAGas = poisson(nBeamAGasAvg, pythiaBeamAGas.rndm);
 130     nAGH.fill( nBeamAGas );
 131
 132     // Generate a number of beam A + gas events. Add them to sumEvent.
 133     for (int iAG = 0; iAG < nBeamAGas; ++iAG) {
 134       pythiaBeamAGas.next();
 135       sumEvent += pythiaBeamAGas.event;
 136     }
 137
 138     // Select the number of beam B + gas events to generate.
 139     int nBeamBGas = poisson(nBeamBGasAvg, pythiaBeamBGas.rndm);
 140     nBGH.fill( nBeamBGas );
 141
 142     // Generate a number of beam B + gas events. Add them to sumEvent.
 143     for (int iBG = 0; iBG < nBeamBGas; ++iBG) {
 144       pythiaBeamBGas.next();
 145       sumEvent += pythiaBeamBGas.event;
 146     }
 147
 148     // List first few events.
 149     if (iEvent < 1) {
 150       pythiaSignal.info.list();
 151       pythiaSignal.process.list();
 152       sumEvent.list();
 153     }
 154
 155     // Find charged multiplicity.
 156     int nChg = 0;
 157     for (int i = 0; i < sumEvent.size(); ++i)
 158       if (sumEvent[i].isFinal() && sumEvent[i].isCharged()) ++nChg;
 159     nChgH.fill( nChg );
 160
 161     // Fill net pZ - nonvanishing owing to beam + gas.
 162     sumPZH.fill( sumEvent[0].pz() );
 163
 164   // End of event loop
 165   }
 166
 167   // Statistics. Histograms.
 168   pythiaSignal.statistics();
 169   pythiaPileup.statistics();
 170   pythiaBeamAGas.statistics();
 171   pythiaBeamBGas.statistics();
 172   cout << nPileH << nAGH << nBGH << nChgH << sumPZH;
 173
 174   return 0;
 175 }