PYTHIA8/pythia8175/examples/main18.cc

   1 // main18.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2013 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 "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 }