--- /dev/null
+// main18.cc is a part of the PYTHIA event generator.
+// Copyright (C) 2012 Torbjorn Sjostrand.
+// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
+// Please respect the MCnet Guidelines, see GUIDELINES for details.
+
+// This is a simple test program.
+// It illustrates how to write an event filter.
+// No new functionality is involved - all could be done in the main program
+// - but the division of tasks may be more convenient for recurrent cuts.
+
+#include "Pythia.h"
+
+using namespace Pythia8;
+
+//==========================================================================
+
+// The EventFilter class.
+
+// The constructor takes the following arguments
+// select = 1 : keep only final particles.
+// = 2 : keep only final visible particles (i.e. not neutrinos).
+// = 3 : keep only final charged particles.
+// etaMax (default = 50) : keep only particles with pseudorapidity
+// |eta| < etaMax.
+// pTminCharged (default = 0) : keep a charged particle only if
+// its transverse momentum pT < pTminCharged.
+// pTminNeutral (default = 0) : keep a neutral particle only if
+// its transverse momentum pT < pTminNeutral.
+
+// Main methods:
+// filter( event) takes an event record as input and analyzes it.
+// size() returns the number of particles kept.
+// index(i) returns the index in the full event of the i'th kept particle.
+// particlePtr(i) returns a pointer to the i'th kept particle.
+// particleRef(i) returns a reference to the i'th kept particle.
+// list() gives a listing of the kept particles only.
+
+class EventFilter {
+
+public:
+
+ // Constructor sets properties of filter.
+ EventFilter( int selectIn, double etaMaxIn = 50.,
+ double pTminChargedIn = 0., double pTminNeutralIn = 0.)
+ : select(selectIn), etaMax(etaMaxIn), pTminCharged(pTminChargedIn),
+ pTminNeutral(pTminNeutralIn) {}
+
+ // Analysis of each new event to find acceptable particles.
+ void filter(Event& event);
+
+ // Return size of array, and index of a particle.
+ int size() const {return keptPtrs.size();}
+ int index(int i) const {return keptIndx[i];}
+
+ // Return pointer or reference to a particle.
+ Particle* particlePtr(int i) {return keptPtrs[i];}
+ Particle& particleRef(int i) {return *keptPtrs[i];}
+
+ // List kept particles only.
+ void list(ostream& os = cout);
+
+private:
+
+ // Filter properties, set by constructor.
+ int select;
+ double etaMax, pTminCharged, pTminNeutral;
+
+ // Kept particle indices and pointers, referring to original event.
+ vector<int> keptIndx;
+ vector<Particle*> keptPtrs;
+
+};
+
+//--------------------------------------------------------------------------
+
+// The filter method.
+
+void EventFilter::filter(Event& event) {
+
+ // Reset arrays in preparation for new event.
+ keptIndx.resize(0);
+ keptPtrs.resize(0);
+
+ // Loop over all particles in the event record.
+ for (int i = 0; i < event.size(); ++i) {
+
+ // Skip if particle kind selection criteria not fulfilled.
+ if (!event[i].isFinal()) continue;
+ if (select == 2 && !event[i].isVisible()) continue;
+ bool isCharged = event[i].isCharged();
+ if (select == 3 && !isCharged) continue;
+
+ // Skip if too large pseudorapidity.
+ if (abs(event[i].eta()) > etaMax) continue;
+
+ // Skip if too small pT.
+ if (isCharged && event[i].pT() < pTminCharged) continue;
+ else if (!isCharged && event[i].pT() < pTminNeutral) continue;
+
+ // Add particle to vectors of indices and pointers.
+ keptIndx.push_back( i );
+ keptPtrs.push_back( &event[i] );
+
+ // End of particle loop. Done.
+ }
+
+}
+
+//--------------------------------------------------------------------------
+
+// The list method: downscaled version of Event::list.
+
+void EventFilter::list(ostream& os) {
+
+ // Header.
+ os << "\n -------- PYTHIA Event Listing (filtered) ------------------"
+ << "-----------------------------------------------------------------"
+ << "----\n \n no id name status mothers "
+ << " daughters colours p_x p_y p_z e "
+ << " m \n";
+
+ // At high energy switch to scientific format for momenta.
+ double eSum = 0.;
+ for (int iKept = 0; iKept < size(); ++iKept) eSum += keptPtrs[iKept]->e();
+ bool useFixed = (eSum < 1e5);
+
+ // Listing of kept particles in event.
+ for (int iKept = 0; iKept < size(); ++iKept) {
+ int i = keptIndx[iKept];
+ Particle& pt = *keptPtrs[iKept];
+
+ // Basic line for a particle, always printed.
+ os << setw(6) << i << setw(10) << pt.id() << " " << left
+ << setw(18) << pt.nameWithStatus(18) << right << setw(4)
+ << pt.status() << setw(6) << pt.mother1() << setw(6)
+ << pt.mother2() << setw(6) << pt.daughter1() << setw(6)
+ << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol()
+ << ( (useFixed) ? fixed : scientific ) << setprecision(3)
+ << setw(11) << pt.px() << setw(11) << pt.py() << setw(11)
+ << pt.pz() << setw(11) << pt.e() << setw(11) << pt.m() << "\n";
+ }
+
+ // Listing finished.
+ os << "\n -------- End PYTHIA Event Listing ----------------------------"
+ << "-------------------------------------------------------------------"
+ << endl;
+}
+
+
+//==========================================================================
+
+// Use the EventFilter method to plot some event properties.
+
+int main() {
+
+ // Number of events to generate, to list, to allow aborts.
+ int nEvent = 100;
+ int nList = 1;
+ int nAbort = 3;
+
+ // Declare generator.
+ Pythia pythia;
+
+ // Hard QCD events with pThat > 100.
+ pythia.readString("HardQCD:all = on");
+ pythia.readString("PhaseSpace:pTHatMin = 100.");
+
+ // No automatic event listings - do it manually below.
+ pythia.readString("Next:numberShowInfo = 0");
+ pythia.readString("Next:numberShowProcess = 0");
+ pythia.readString("Next:numberShowEvent = 0");
+
+ // Initialization for LHC.
+ pythia.init();
+
+ // Values for filter.
+ int select = 3;
+ double etaMax = 3.;
+ double pTminChg = 1.;
+
+ // Declare Event Filter according to specification.
+ EventFilter filter( select, etaMax, pTminChg);
+
+ // Histograms.
+ Hist nCharged( "selected charged multiplicity", 100, -0.5, 199.5);
+ Hist etaCharged( "selected charged eta distribution", 100, -5.0, 5.0);
+ Hist pTCharged( "selected charged pT distribution", 100, 0.0, 50.0);
+
+ // Begin event loop.
+ int iAbort = 0;
+ for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
+
+ // Generate events. Quit if too many failures.
+ if (!pythia.next()) {
+ if (++iAbort < nAbort) continue;
+ cout << " Event generation aborted prematurely, owing to error!\n";
+ break;
+ }
+
+ // Find final charged particles with |eta| < 3 and pT > 1 GeV.
+ filter.filter( pythia.event);
+
+ // List first few events, both complete and after filtering.
+ if (iEvent < nList) {
+ pythia.info.list();
+ pythia.process.list();
+ pythia.event.list();
+ filter.list();
+ }
+
+ // Analyze selected particle sample.
+ nCharged.fill( filter.size() );
+ for (int i = 0; i < filter.size(); ++i) {
+ // Use both reference and pointer notation to illustrate freedom.
+ etaCharged.fill( filter.particleRef(i).eta() );
+ pTCharged.fill( filter.particlePtr(i)->pT() );
+ }
+
+ // End of event loop.
+ }
+
+ // Final statistics.
+ pythia.stat();
+
+ // Histograms.
+ cout << nCharged << etaCharged << pTCharged;
+
+ // Done.
+ return 0;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff