1 // main18.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 is a simple test program.
7 // It illustrates how to write an event filter.
8 // No new functionality is involved - all could be done in the main program
9 // - but the division of tasks may be more convenient for recurrent cuts.
13 using namespace Pythia8;
15 //==========================================================================
17 // The EventFilter class.
19 // The constructor takes the following arguments
20 // select = 1 : keep only final particles.
21 // = 2 : keep only final visible particles (i.e. not neutrinos).
22 // = 3 : keep only final charged particles.
23 // etaMax (default = 50) : keep only particles with pseudorapidity
25 // pTminCharged (default = 0) : keep a charged particle only if
26 // its transverse momentum pT < pTminCharged.
27 // pTminNeutral (default = 0) : keep a neutral particle only if
28 // its transverse momentum pT < pTminNeutral.
31 // filter( event) takes an event record as input and analyzes it.
32 // size() returns the number of particles kept.
33 // index(i) returns the index in the full event of the i'th kept particle.
34 // particlePtr(i) returns a pointer to the i'th kept particle.
35 // particleRef(i) returns a reference to the i'th kept particle.
36 // list() gives a listing of the kept particles only.
42 // Constructor sets properties of filter.
43 EventFilter( int selectIn, double etaMaxIn = 50.,
44 double pTminChargedIn = 0., double pTminNeutralIn = 0.)
45 : select(selectIn), etaMax(etaMaxIn), pTminCharged(pTminChargedIn),
46 pTminNeutral(pTminNeutralIn) {}
48 // Analysis of each new event to find acceptable particles.
49 void filter(Event& event);
51 // Return size of array, and index of a particle.
52 int size() const {return keptPtrs.size();}
53 int index(int i) const {return keptIndx[i];}
55 // Return pointer or reference to a particle.
56 Particle* particlePtr(int i) {return keptPtrs[i];}
57 Particle& particleRef(int i) {return *keptPtrs[i];}
59 // List kept particles only.
60 void list(ostream& os = cout);
64 // Filter properties, set by constructor.
66 double etaMax, pTminCharged, pTminNeutral;
68 // Kept particle indices and pointers, referring to original event.
70 vector<Particle*> keptPtrs;
74 //--------------------------------------------------------------------------
78 void EventFilter::filter(Event& event) {
80 // Reset arrays in preparation for new event.
84 // Loop over all particles in the event record.
85 for (int i = 0; i < event.size(); ++i) {
87 // Skip if particle kind selection criteria not fulfilled.
88 if (!event[i].isFinal()) continue;
89 if (select == 2 && !event[i].isVisible()) continue;
90 bool isCharged = event[i].isCharged();
91 if (select == 3 && !isCharged) continue;
93 // Skip if too large pseudorapidity.
94 if (abs(event[i].eta()) > etaMax) continue;
96 // Skip if too small pT.
97 if (isCharged && event[i].pT() < pTminCharged) continue;
98 else if (!isCharged && event[i].pT() < pTminNeutral) continue;
100 // Add particle to vectors of indices and pointers.
101 keptIndx.push_back( i );
102 keptPtrs.push_back( &event[i] );
104 // End of particle loop. Done.
109 //--------------------------------------------------------------------------
111 // The list method: downscaled version of Event::list.
113 void EventFilter::list(ostream& os) {
116 os << "\n -------- PYTHIA Event Listing (filtered) ------------------"
117 << "-----------------------------------------------------------------"
118 << "----\n \n no id name status mothers "
119 << " daughters colours p_x p_y p_z e "
122 // At high energy switch to scientific format for momenta.
124 for (int iKept = 0; iKept < size(); ++iKept) eSum += keptPtrs[iKept]->e();
125 bool useFixed = (eSum < 1e5);
127 // Listing of kept particles in event.
128 for (int iKept = 0; iKept < size(); ++iKept) {
129 int i = keptIndx[iKept];
130 Particle& pt = *keptPtrs[iKept];
132 // Basic line for a particle, always printed.
133 os << setw(6) << i << setw(10) << pt.id() << " " << left
134 << setw(18) << pt.nameWithStatus(18) << right << setw(4)
135 << pt.status() << setw(6) << pt.mother1() << setw(6)
136 << pt.mother2() << setw(6) << pt.daughter1() << setw(6)
137 << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol()
138 << ( (useFixed) ? fixed : scientific ) << setprecision(3)
139 << setw(11) << pt.px() << setw(11) << pt.py() << setw(11)
140 << pt.pz() << setw(11) << pt.e() << setw(11) << pt.m() << "\n";
144 os << "\n -------- End PYTHIA Event Listing ----------------------------"
145 << "-------------------------------------------------------------------"
150 //==========================================================================
152 // Use the EventFilter method to plot some event properties.
156 // Number of events to generate, to list, to allow aborts.
161 // Declare generator.
164 // Hard QCD events with pThat > 100.
165 pythia.readString("HardQCD:all = on");
166 pythia.readString("PhaseSpace:pTHatMin = 100.");
168 // No automatic event listings - do it manually below.
169 pythia.readString("Next:numberShowInfo = 0");
170 pythia.readString("Next:numberShowProcess = 0");
171 pythia.readString("Next:numberShowEvent = 0");
173 // Initialization for LHC.
176 // Values for filter.
179 double pTminChg = 1.;
181 // Declare Event Filter according to specification.
182 EventFilter filter( select, etaMax, pTminChg);
185 Hist nCharged( "selected charged multiplicity", 100, -0.5, 199.5);
186 Hist etaCharged( "selected charged eta distribution", 100, -5.0, 5.0);
187 Hist pTCharged( "selected charged pT distribution", 100, 0.0, 50.0);
191 for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
193 // Generate events. Quit if too many failures.
194 if (!pythia.next()) {
195 if (++iAbort < nAbort) continue;
196 cout << " Event generation aborted prematurely, owing to error!\n";
200 // Find final charged particles with |eta| < 3 and pT > 1 GeV.
201 filter.filter( pythia.event);
203 // List first few events, both complete and after filtering.
204 if (iEvent < nList) {
206 pythia.process.list();
211 // Analyze selected particle sample.
212 nCharged.fill( filter.size() );
213 for (int i = 0; i < filter.size(); ++i) {
214 // Use both reference and pointer notation to illustrate freedom.
215 etaCharged.fill( filter.particleRef(i).eta() );
216 pTCharged.fill( filter.particlePtr(i)->pT() );
219 // End of event loop.
226 cout << nCharged << etaCharged << pTCharged;