1
// main19.cc is a part of the PYTHIA event generator.
2
// Copyright (C) 2012 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.
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.
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.
23
using namespace Pythia8;
25
//==========================================================================
27
// Method to pick a number according to a Poissonian distribution.
29
int poisson(double nAvg, Rndm& rndm) {
31
// Set maximum to avoid overflow.
35
double rPoisson = rndm.flat() * exp(nAvg);
41
// Add to sum and check whether done.
42
for (int i = 0; i < NMAX; ) {
44
if (rSum > rPoisson) return i;
46
// Evaluate next term.
55
//==========================================================================
59
// Number of signal events to generate.
65
// Average number of pileup events per signal event.
66
double nPileupAvg = 2.5;
68
// Average number of beam-gas events per signal ones, on two sides.
69
double nBeamAGasAvg = 0.5;
70
double nBeamBGasAvg = 0.5;
72
// Four generator instances.
75
Pythia pythiaBeamAGas;
76
Pythia pythiaBeamBGas;
78
// One object where all individual events are to be collected.
81
// Switch off automatic event listing.
82
pythiaSignal.readString("Next:numberShowInfo = 0");
83
pythiaSignal.readString("Next:numberShowProcess = 0");
84
pythiaSignal.readString("Next:numberShowEvent = 0");
85
pythiaPileup.readString("Next:numberShowInfo = 0");
86
pythiaPileup.readString("Next:numberShowProcess = 0");
87
pythiaPileup.readString("Next:numberShowEvent = 0");
88
pythiaBeamAGas.readString("Next:numberShowInfo = 0");
89
pythiaBeamAGas.readString("Next:numberShowProcess = 0");
90
pythiaBeamAGas.readString("Next:numberShowEvent = 0");
91
pythiaBeamBGas.readString("Next:numberShowInfo = 0");
92
pythiaBeamBGas.readString("Next:numberShowProcess = 0");
93
pythiaBeamBGas.readString("Next:numberShowEvent = 0");
95
// Initialize generator for signal processes.
96
pythiaSignal.readString("HardQCD:all = on");
97
pythiaSignal.readString("PhaseSpace:pTHatMin = 50.");
98
pythiaSignal.settings.parm("Beams:eCM", 2. * eBeam);
101
// Initialize generator for pileup (background) processes.
102
pythiaPileup.readString("Random:setSeed = on");
103
pythiaPileup.readString("Random:seed = 10000002");
104
pythiaPileup.readString("SoftQCD:all = on");
105
pythiaPileup.settings.parm("Beams:eCM", 2. * eBeam);
108
// Initialize generators for beam A - gas (background) processes.
109
pythiaBeamAGas.readString("Random:setSeed = on");
110
pythiaBeamAGas.readString("Random:seed = 10000003");
111
pythiaBeamAGas.readString("SoftQCD:all = on");
112
pythiaBeamAGas.readString("Beams:frameType = 2");
113
pythiaBeamAGas.settings.parm("Beams:eA", eBeam);
114
pythiaBeamAGas.settings.parm("Beams:eB", 0.);
115
pythiaBeamAGas.init();
117
// Initialize generators for beam B - gas (background) processes.
118
pythiaBeamBGas.readString("Random:setSeed = on");
119
pythiaBeamBGas.readString("Random:seed = 10000004");
120
pythiaBeamBGas.readString("SoftQCD:all = on");
121
pythiaBeamBGas.readString("Beams:frameType = 2");
122
pythiaBeamBGas.settings.parm("Beams:eA", 0.);
123
pythiaBeamBGas.settings.parm("Beams:eB", eBeam);
124
pythiaBeamBGas.init();
126
// Histograms: number of pileups, total charged multiplicity.
127
Hist nPileH("number of pileup events per signal event", 100, -0.5, 99.5);
128
Hist nAGH("number of beam A + gas events per signal event", 100, -0.5, 99.5);
129
Hist nBGH("number of beam B + gas events per signal event", 100, -0.5, 99.5);
130
Hist nChgH("number of charged multiplicity",100, -0.5, 1999.5);
131
Hist sumPZH("total pZ of system",100, -100000., 100000.);
134
for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
136
// Generate a signal event. Copy this event into sumEvent.
137
if (!pythiaSignal.next()) continue;
138
sumEvent = pythiaSignal.event;
140
// Select the number of pileup events to generate.
141
int nPileup = poisson(nPileupAvg, pythiaPileup.rndm);
142
nPileH.fill( nPileup );
144
// Generate a number of pileup events. Add them to sumEvent.
145
for (int iPileup = 0; iPileup < nPileup; ++iPileup) {
147
sumEvent += pythiaPileup.event;
150
// Select the number of beam A + gas events to generate.
151
int nBeamAGas = poisson(nBeamAGasAvg, pythiaBeamAGas.rndm);
152
nAGH.fill( nBeamAGas );
154
// Generate a number of beam A + gas events. Add them to sumEvent.
155
for (int iAG = 0; iAG < nBeamAGas; ++iAG) {
156
pythiaBeamAGas.next();
157
sumEvent += pythiaBeamAGas.event;
160
// Select the number of beam B + gas events to generate.
161
int nBeamBGas = poisson(nBeamBGasAvg, pythiaBeamBGas.rndm);
162
nBGH.fill( nBeamBGas );
164
// Generate a number of beam B + gas events. Add them to sumEvent.
165
for (int iBG = 0; iBG < nBeamBGas; ++iBG) {
166
pythiaBeamBGas.next();
167
sumEvent += pythiaBeamBGas.event;
170
// List first few events.
172
pythiaSignal.info.list();
173
pythiaSignal.process.list();
177
// Find charged multiplicity.
179
for (int i = 0; i < sumEvent.size(); ++i)
180
if (sumEvent[i].isFinal() && sumEvent[i].isCharged()) ++nChg;
183
// Fill net pZ - nonvanishing owing to beam + gas.
184
sumPZH.fill( sumEvent[0].pz() );
189
// Statistics. Histograms.
192
pythiaBeamAGas.stat();
193
pythiaBeamBGas.stat();
194
cout << nPileH << nAGH << nBGH << nChgH << sumPZH;
added
removed
examples/main19.cc
