PYTHIA8/pythia8175/examples/main06.cc

   1 // main06.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 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   pythia.readString("Beams:idA =  11");
  28   pythia.readString("Beams:idB = -11");
  29   double mZ = pythia.particleData.m0(23);
  30   pythia.settings.parm("Beams:eCM", mZ);
  31   pythia.init();
  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
  61     // Find and histogram charged multiplicity.
  62     int nCh = 0;
  63     for (int i = 0; i < pythia.event.size(); ++i)
  64       if (pythia.event[i].isFinal() && pythia.event[i].isCharged()) ++nCh;
  65     nCharge.fill( nCh );
  66
  67     // Find and histogram sphericity.
  68     if (sph.analyze( pythia.event )) {
  69       if (iEvent < 3) sph.list();
  70       spheri.fill( sph.sphericity() );
  71       sAxis.fill( sph.eventAxis(1).pz() );
  72       double e1 = sph.eigenValue(1);
  73       double e2 = sph.eigenValue(2);
  74       double e3 = sph.eigenValue(3);
  75       if (e2 > e1 || e3 > e2) cout << "eigenvalues out of order: "
  76       << e1 << "  " << e2 << "  " << e3 << endl;
  77     }
  78
  79     // Find and histogram linearized sphericity.
  80     if (lin.analyze( pythia.event )) {
  81       if (iEvent < 3) lin.list();
  82       linea.fill( lin.sphericity() );
  83       lAxis.fill( lin.eventAxis(1).pz() );
  84       double e1 = lin.eigenValue(1);
  85       double e2 = lin.eigenValue(2);
  86       double e3 = lin.eigenValue(3);
  87       if (e2 > e1 || e3 > e2) cout << "eigenvalues out of order: "
  88       << e1 << "  " << e2 << "  " << e3 << endl;
  89     }
  90
  91     // Find and histogram thrust.
  92     if (thr.analyze( pythia.event )) {
  93       if (iEvent < 3) thr.list();
  94       thrust.fill( thr.thrust() );
  95       oblateness.fill( thr.oblateness() );
  96       tAxis.fill( thr.eventAxis(1).pz() );
  97       if ( abs(thr.eventAxis(1).pAbs() - 1.) > 1e-8
  98         || abs(thr.eventAxis(2).pAbs() - 1.) > 1e-8
  99         || abs(thr.eventAxis(3).pAbs() - 1.) > 1e-8
 100         || abs(thr.eventAxis(1) * thr.eventAxis(2)) > 1e-8
 101         || abs(thr.eventAxis(1) * thr.eventAxis(3)) > 1e-8
 102         || abs(thr.eventAxis(2) * thr.eventAxis(3)) > 1e-8 ) {
 103         cout << " suspicious Thrust eigenvectors " << endl;
 104         thr.list();
 105       }
 106     }
 107
 108     // Find and histogram cluster jets: Lund, Jade and Durham distance.
 109     if (lund.analyze( pythia.event, 0.01, 0.)) {
 110       if (iEvent < 3) lund.list();
 111       nLund.fill( lund.size() );
 112       for (int j = 0; j < lund.size() - 1; ++j)
 113         eDifLund.fill( lund.p(j).e() - lund.p(j+1).e() );
 114     }
 115     if (jade.analyze( pythia.event, 0.01, 0.)) {
 116       if (iEvent < 3) jade.list();
 117       nJade.fill( jade.size() );
 118       for (int j = 0; j < jade.size() - 1; ++j)
 119         eDifJade.fill( jade.p(j).e() - jade.p(j+1).e() );
 120     }
 121     if (durham.analyze( pythia.event, 0.01, 0.)) {
 122       if (iEvent < 3) durham.list();
 123       nDurham.fill( durham.size() );
 124       for (int j = 0; j < durham.size() - 1; ++j)
 125         eDifDurham.fill( durham.p(j).e() - durham.p(j+1).e() );
 126     }
 127
 128   // End of event loop. Statistics. Output histograms.
 129   }
 130   pythia.stat();
 131   cout << nCharge << spheri << linea << thrust << oblateness
 132        << sAxis << lAxis << tAxis
 133        << nLund << nJade << nDurham
 134        << eDifLund << eDifJade << eDifDurham;
 135
 136   // Done.
 137   return 0;
 138 }