PYTHIA8/pythia8145/examples/main09.cc

   1 // main09.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 studies event properties of LEP1 events.
   8
   9 #include "Pythia.h"
  10
  11 using namespace Pythia8;
  12
  13 int main() {
  14
  15   // Generator.
  16   Pythia pythia;
  17
  18   // Allow no substructure in e+- beams: normal for corrected LEP data.
  19   pythia.readString("PDF:lepton = off");
  20   // Process selection.
  21   pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
  22   // Switch off all Z0 decays and then switch back on those to quarks.
  23   pythia.readString("23:onMode = off");
  24   pythia.readString("23:onIfAny = 1 2 3 4 5");
  25
  26   // LEP1 initialization at Z0 mass.
  27   double mZ = pythia.particleData.m0(23);
  28   pythia.init( 11, -11, mZ);
  29
  30   // Check that Z0 decay channels set correctly.
  31   pythia.particleData.listChanged();
  32
  33   // Histograms.
  34   Hist nCharge("charged multiplicity", 100, -0.5, 99.5);
  35   Hist spheri("Sphericity", 100, 0., 1.);
  36   Hist linea("Linearity", 100, 0., 1.);
  37   Hist thrust("thrust", 100, 0.5, 1.);
  38   Hist oblateness("oblateness", 100, 0., 1.);
  39   Hist sAxis("cos(theta_Sphericity)", 100, -1., 1.);
  40   Hist lAxis("cos(theta_Linearity)", 100, -1., 1.);
  41   Hist tAxis("cos(theta_Thrust)", 100, -1., 1.);
  42   Hist nLund("Lund jet multiplicity", 40, -0.5, 39.5);
  43   Hist nJade("Jade jet multiplicity", 40, -0.5, 39.5);
  44   Hist nDurham("Durham jet multiplicity", 40, -0.5, 39.5);
  45   Hist eDifLund("Lund e_i - e_{i+1}", 100, -5.,45.);
  46   Hist eDifJade("Jade e_i - e_{i+1}", 100, -5.,45.);
  47   Hist eDifDurham("Durham e_i - e_{i+1}", 100, -5.,45.);
  48
  49   // Set up Sphericity, "Linearity", Thrust and cluster jet analyses.
  50   Sphericity sph;
  51   Sphericity lin(1.);
  52   Thrust thr;
  53   ClusterJet lund("Lund");
  54   ClusterJet jade("Jade");
  55   ClusterJet durham("Durham");
  56
  57   // Begin event loop. Generate event. Skip if error. List first few.
  58   for (int iEvent = 0; iEvent < 10000; ++iEvent) {
  59     if (!pythia.next()) continue;
  60     if (iEvent < 1) pythia.event.list();
  61
  62     // Find and histogram charged multiplicity.
  63     int nCh = 0;
  64     for (int i = 0; i < pythia.event.size(); ++i)
  65       if (pythia.event[i].isFinal() && pythia.event[i].isCharged()) ++nCh;
  66     nCharge.fill( nCh );
  67
  68     // Find and histogram sphericity.
  69     if (sph.analyze( pythia.event )) {
  70       if (iEvent < 3) sph.list();
  71       spheri.fill( sph.sphericity() );
  72       sAxis.fill( sph.eventAxis(1).pz() );
  73       double e1 = sph.eigenValue(1);
  74       double e2 = sph.eigenValue(2);
  75       double e3 = sph.eigenValue(3);
  76       if (e2 > e1 || e3 > e2) cout << "eigenvalues out of order: "
  77       << e1 << "  " << e2 << "  " << e3 << endl;
  78     }
  79
  80     // Find and histogram linearized sphericity.
  81     if (lin.analyze( pythia.event )) {
  82       if (iEvent < 3) lin.list();
  83       linea.fill( lin.sphericity() );
  84       lAxis.fill( lin.eventAxis(1).pz() );
  85       double e1 = lin.eigenValue(1);
  86       double e2 = lin.eigenValue(2);
  87       double e3 = lin.eigenValue(3);
  88       if (e2 > e1 || e3 > e2) cout << "eigenvalues out of order: "
  89       << e1 << "  " << e2 << "  " << e3 << endl;
  90     }
  91
  92     // Find and histogram thrust.
  93     if (thr.analyze( pythia.event )) {
  94       if (iEvent < 3) thr.list();
  95       thrust.fill( thr.thrust() );
  96       oblateness.fill( thr.oblateness() );
  97       tAxis.fill( thr.eventAxis(1).pz() );
  98       if ( abs(thr.eventAxis(1).pAbs() - 1.) > 1e-8
  99         || abs(thr.eventAxis(2).pAbs() - 1.) > 1e-8
 100         || abs(thr.eventAxis(3).pAbs() - 1.) > 1e-8
 101         || abs(thr.eventAxis(1) * thr.eventAxis(2)) > 1e-8
 102         || abs(thr.eventAxis(1) * thr.eventAxis(3)) > 1e-8
 103         || abs(thr.eventAxis(2) * thr.eventAxis(3)) > 1e-8 ) {
 104         cout << " suspicious Thrust eigenvectors " << endl;
 105         thr.list();
 106       }
 107     }
 108
 109     // Find and histogram cluster jets: Lund, Jade and Durham distance.
 110     if (lund.analyze( pythia.event, 0.01, 0.)) {
 111       if (iEvent < 3) lund.list();
 112       nLund.fill( lund.size() );
 113       for (int j = 0; j < lund.size() - 1; ++j)
 114         eDifLund.fill( lund.p(j).e() - lund.p(j+1).e() );
 115     }
 116     if (jade.analyze( pythia.event, 0.01, 0.)) {
 117       if (iEvent < 3) jade.list();
 118       nJade.fill( jade.size() );
 119       for (int j = 0; j < jade.size() - 1; ++j)
 120         eDifJade.fill( jade.p(j).e() - jade.p(j+1).e() );
 121     }
 122     if (durham.analyze( pythia.event, 0.01, 0.)) {
 123       if (iEvent < 3) durham.list();
 124       nDurham.fill( durham.size() );
 125       for (int j = 0; j < durham.size() - 1; ++j)
 126         eDifDurham.fill( durham.p(j).e() - durham.p(j+1).e() );
 127     }
 128
 129   // End of event loop. Statistics. Output histograms.
 130   }
 131   pythia.statistics();
 132   cout << nCharge << spheri << linea << thrust << oblateness
 133        << sAxis << lAxis << tAxis
 134        << nLund << nJade << nDurham
 135        << eDifLund << eDifJade << eDifDurham;
 136
 137   // Done.
 138   return 0;
 139 }