main200

Back to index.

// main200.cc is a part of the PYTHIA event generator.
// Copyright (C) 2020 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This is a simple test program to run Z decays at LEP I, with some basic
// event shapes, spectra, and multiplicity counts.

// Authors:
//            Peter Skands

// Keywords:
//            Vincia
//            Dire
//            Electron‑positron

// Include Pythia8 header(s) and namespace.
#include "Pythia8/Pythia.h"
using namespace Pythia8;

// Main Program
int main() {

  //************************************************************************
  // Define Pythia 8 generator

  Pythia pythia;

  // Read user settings from file
  pythia.readFile("main200.cmnd");

  //************************************************************************

  // Shorthands
  Event& event = pythia.event;
  Settings& settings = pythia.settings;

  // Extract settings to be used in the main program.
  int    nEvent     = settings.mode("Main:numberOfEvents");
  int    nAbort     = settings.mode("Main:timesAllowErrors");
  bool   vinciaOn   = settings.mode("PartonShowers:model") == 2;
  bool   helicityOn = vinciaOn && settings.flag("Vincia:helicityShower");
  int    iEventPri  = -1;
  double eCM        = settings.parm("Beams:eCM");

  //************************************************************************

  // Initialize
  if(!pythia.init()) { return EXIT_FAILURE; }

  //************************************************************************

  // Define a few PYTHIA utilities
  Thrust Thr(1);
  Sphericity SphLin(1,2);

  //************************************************************************

  // Define PYTHIA histograms
  // Rootlink: to use ROOT instead, HIST -> TH1F ("tag",...)

  // Number of quarks, partons, and charged particles
  double nQsum      = 0.0;
  double nPsum      = 0.0;
  double nCsum      = 0.0;
  double nStrangeSum = 0.0;
  double nCharmSum  = 0.0;
  double nBottomSum = 0.0;
  double nBaryonSum = 0.0;
  double nGammaSum  = 0.0;
  Hist histNQuarks("nQuarks", 50, -0.5, 49.5);
  Hist histNPartons("nPartons", 100, -0.5, 99.5);
  Hist histNCharged("nCharged", 100, -1.0, 199.0);
  Hist histNGamma("nPhotons", 50, -0.5, 49.5);
  Hist histXStrange("Ln(x) for strange quarks", 25, -5.0, 0.0);
  Hist histXCharm("Ln(x) for charm quarks", 25, -5.0, 0.0);
  Hist histXBottom("Ln(x) for bottom quarks", 25, -5.0, 0.0);
  Hist histMSS("Invariant mass of s-sbar pairs",50,0.,10.);
  Hist histMCC("Invariant mass of c-cbar pairs",50,0.,25.);
  Hist histMBB("Invariant mass of b-bbar pairs",50,0.,50.);
  Hist histX1("Ln(x) for 1st branching (QCD)",100,-7.5,0.0);
  Hist histX1Gamma("Ln(x) for 1st branching (QED)",100,-7.5,0.0);
  Hist histX2Gluon("Ln(x) for 2nd branching (gluons)",100,-7.5,0.0);
  Hist histX2Quark("Ln(x) for 2nd branching (quarks)",100,-7.5,0.0);

  // Thrust, C, and D parameters
  int nBinsShapes=100;
  Hist histT("1-T",nBinsShapes, 0.0, 0.5);
  Hist histC("C",nBinsShapes,0.0,1.0);
  Hist histD("D",nBinsShapes,0.0,1.0);
  double wHistT=nBinsShapes/0.5;
  double wHistCD=nBinsShapes/1.0;

  //************************************************************************

  // EVENT GENERATION LOOP
  // Generation, event-by-event printout, analysis, and storage

  // (Optionally) re-initialize random number engine before event loop
  const int initseed = settings.mode("Random:seed");
  pythia.rndm.init(initseed);


  // Counter for negative-weight events
  double weight=1.0;
  double sumWeights = 0.0;

  // Begin event loop
  int iAbort = 0;
  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {

    bool aborted = !pythia.next();
    if(aborted){
      event.list();
      if (++iAbort < nAbort){
        continue;
      }
      cout << " Event generation aborted prematurely, owing to error!\n";
      cout<< "Event number was : "< iStrange;
    vector iCharm;
    vector iBottom;
    for (int i=1;i 80) nc++;
      // Find last parton-level partons
      int iDau1 = event[i].daughter1();
      if (iDau1 == 0 || abs(event[iDau1].status()) > 80) {
        // Count up partons and quarks
        if (event[i].isQuark() || event[i].isGluon()) np++;
        if (event[i].id() == 22) nGam++;
        if (event[i].isQuark()) nq++;
        if (event[i].idAbs() == 3) {
          nStrange++;
          histXStrange.fill(log(2*event[i].pAbs()/eCM),weight);
          for (int j=0; j<(int)iStrange.size(); ++j) {
            double m2qq = m2(event[i],event[iStrange[j]]);
            histMSS.fill(sqrt(m2qq),weight);
          }
          iStrange.push_back(i);
        }
        if (event[i].idAbs() == 4) {
          nCharm++;
          histXCharm.fill(log(2*event[i].pAbs()/eCM),weight);
          for (int j=0; j<(int)iCharm.size(); ++j) {
            double m2qq = m2(event[i],event[iCharm[j]]);
            histMCC.fill(sqrt(m2qq),weight);
          }
          iCharm.push_back(i);
        }
        if (event[i].idAbs() == 5) {
          nBottom++;
          histXBottom.fill(log(2*event[i].pAbs()/eCM),weight);
          for (int j=0; j<(int)iBottom.size(); ++j) {
            double m2qq = m2(event[i],event[iBottom[j]]);
            histMBB.fill(sqrt(m2qq),weight);
          }
          iBottom.push_back(i);
        }
      }
      if (event[i].isHadron() && event[i].isFinal()) {
        int idAbs = abs(event[i].id());
        if (idAbs > 1000 && idAbs < 10000) nBaryonSum++;
      }
    }
    histNQuarks.fill(0.5*nq,weight);
    histNPartons.fill(1.0*np,weight);
    histNCharged.fill(1.0*nc,weight);
    histNGamma.fill(1.0*nGam,weight);
    nQsum += 0.5*nq * weight;
    nPsum += np * weight;
    nCsum += nc * weight;
    nStrangeSum += (nStrange - nStrangeBorn) * weight;
    nCharmSum   += (nCharm - nCharmBorn) * weight;
    nBottomSum  += (nBottom - nBottomBorn) * weight;
    nGammaSum   += nGam * weight;

    // Histogram thrust
    Thr.analyze( event );
    histT.fill(1.0-Thr.thrust(),wHistT*weight);

    // Histogram Linear Sphericity values
    if (np >= 2.0) {
      SphLin.analyze( event );
      double evC=3*(SphLin.eigenValue(1)*SphLin.eigenValue(2)
                    + SphLin.eigenValue(2)*SphLin.eigenValue(3)
                    + SphLin.eigenValue(3)*SphLin.eigenValue(1));
      double evD=27*SphLin.eigenValue(1)*SphLin.eigenValue(2)
      *SphLin.eigenValue(3);
      histC.fill(evC,wHistCD*weight);
      histD.fill(evD,wHistCD*weight);
    }

  }

  //************************************************************************

  // POST-RUN FINALIZATION
  // Normalization, Statistics, Output

  //Normalize histograms to effective number of positive-weight events.
  double normFac=1.0/sumWeights;
  histT          *= normFac;
  histC          *= normFac;
  histD          *= normFac;
  histNQuarks    *= normFac;
  histNPartons   *= normFac;
  histNCharged   *= normFac;
  histNGamma     *= normFac;
  histXStrange   *= normFac;
  histXCharm     *= normFac;
  histXBottom    *= normFac;
  histX1         *= normFac;
  histX1Gamma    *= normFac;
  histX2Gluon    *= normFac;
  histX2Quark    *= normFac;

  // Print a few histograms
  cout< 0) cout< = "<SS>      = "<CC>      = "<BB>      = "<    = "<    = "<    = "<    = "<