StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main18.cc
1 // main18.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 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.
10 
11 #include "Pythia8/Pythia.h"
12 
13 using namespace Pythia8;
14 
15 //==========================================================================
16 
17 // The EventFilter class.
18 
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
24 // |eta| < etaMax.
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.
29 
30 // Main methods:
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.
37 
38 class EventFilter {
39 
40 public:
41 
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) {}
47 
48  // Analysis of each new event to find acceptable particles.
49  void filter(Event& event);
50 
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];}
54 
55  // Return pointer or reference to a particle.
56  Particle* particlePtr(int i) {return keptPtrs[i];}
57  Particle& particleRef(int i) {return *keptPtrs[i];}
58 
59  // List kept particles only.
60  void list(ostream& os = cout);
61 
62 private:
63 
64  // Filter properties, set by constructor.
65  int select;
66  double etaMax, pTminCharged, pTminNeutral;
67 
68  // Kept particle indices and pointers, referring to original event.
69  vector<int> keptIndx;
70  vector<Particle*> keptPtrs;
71 
72 };
73 
74 //--------------------------------------------------------------------------
75 
76 // The filter method.
77 
78 void EventFilter::filter(Event& event) {
79 
80  // Reset arrays in preparation for new event.
81  keptIndx.resize(0);
82  keptPtrs.resize(0);
83 
84  // Loop over all particles in the event record.
85  for (int i = 0; i < event.size(); ++i) {
86 
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;
92 
93  // Skip if too large pseudorapidity.
94  if (abs(event[i].eta()) > etaMax) continue;
95 
96  // Skip if too small pT.
97  if (isCharged && event[i].pT() < pTminCharged) continue;
98  else if (!isCharged && event[i].pT() < pTminNeutral) continue;
99 
100  // Add particle to vectors of indices and pointers.
101  keptIndx.push_back( i );
102  keptPtrs.push_back( &event[i] );
103 
104  // End of particle loop. Done.
105  }
106 
107 }
108 
109 //--------------------------------------------------------------------------
110 
111 // The list method: downscaled version of Event::list.
112 
113 void EventFilter::list(ostream& os) {
114 
115  // Header.
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 "
120  << " m \n";
121 
122  // At high energy switch to scientific format for momenta.
123  double eSum = 0.;
124  for (int iKept = 0; iKept < size(); ++iKept) eSum += keptPtrs[iKept]->e();
125  bool useFixed = (eSum < 1e5);
126 
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];
131 
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";
141  }
142 
143  // Listing finished.
144  os << "\n -------- End PYTHIA Event Listing ----------------------------"
145  << "-------------------------------------------------------------------"
146  << endl;
147 }
148 
149 
150 //==========================================================================
151 
152 // Use the EventFilter method to plot some event properties.
153 
154 int main() {
155 
156  // Number of events to generate, to list, to allow aborts.
157  int nEvent = 100;
158  int nList = 1;
159  int nAbort = 3;
160 
161  // Declare generator.
162  Pythia pythia;
163 
164  // Hard QCD events with pThat > 100.
165  pythia.readString("HardQCD:all = on");
166  pythia.readString("PhaseSpace:pTHatMin = 100.");
167 
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");
172 
173  // Initialization for LHC.
174  pythia.init();
175 
176  // Values for filter.
177  int select = 3;
178  double etaMax = 3.;
179  double pTminChg = 1.;
180 
181  // Declare Event Filter according to specification.
182  EventFilter filter( select, etaMax, pTminChg);
183 
184  // Histograms.
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);
188 
189  // Begin event loop.
190  int iAbort = 0;
191  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
192 
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";
197  break;
198  }
199 
200  // Find final charged particles with |eta| < 3 and pT > 1 GeV.
201  filter.filter( pythia.event);
202 
203  // List first few events, both complete and after filtering.
204  if (iEvent < nList) {
205  pythia.info.list();
206  pythia.process.list();
207  pythia.event.list();
208  filter.list();
209  }
210 
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() );
217  }
218 
219  // End of event loop.
220  }
221 
222  // Final statistics.
223  pythia.stat();
224 
225  // Histograms.
226  cout << nCharged << etaCharged << pTCharged;
227 
228  // Done.
229  return 0;
230 }