main201

Back to index.

// main201.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 compare Pythia and Vincia on
// inclusive jet rates at the LHC, for a sample with pThat > 100 GeV.
// Also illustrates simple use of OpenMP (if enabled) to run two instances
// of Pythia in parallel, here initialised for Pythia and Vincia shower
// models respectively.

// Authors:
//            Peter Skands

// Keywords:
//            Vincia
//            Dire
//            OpenMP

#include "Pythia8/Pythia.h"
using namespace Pythia8;

int main() {

  // Common parameters used for both runs
  const int nEvent    = 1000;
  const int nListJets = 5;

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

  // Histograms.
  Hist nJetsModel1("Model1 number of jets", 20, -0.5, 19.5);
  Hist eTjetsModel1("Model1 pT for jets", 50, 0., 500.);
  Hist yJetsModel1("Model1 y for jets", 20, -4., 4.);
  Hist phiJetsModel1("Model1 phi for jets", 25, -M_PI, M_PI);
  Hist distJetsModel1("Model1 R distance between jets", 50, 0., 10.);
  Hist nJetsModel2("Model2 number of jets", 20, -0.5, 19.5);
  Hist eTjetsModel2("Model2 pT for jets", 50, 0., 500.);
  Hist yJetsModel2("Model2 y for jets", 20, -4., 4.);
  Hist phiJetsModel2("Model2 phi for jets", 25, -M_PI, M_PI);
  Hist distJetsModel2("Model2 R distance between jets", 50, 0., 10.);
  Hist nJetsRatio("Model2/Model1 number of jets", 20, -0.5, 19.5);
  Hist eTjetsRatio("Model2/Model1 pT for jets", 50, 0., 500.);
  Hist yJetsRatio("Model2/Model1 y for jets", 20, -4., 4.);
  Hist phiJetsRatio("Model2/Model1 phi for jets", 25, -M_PI, M_PI);
  Hist distJetsRatio("Model2/Model1 R distance between jets", 50, 0., 10.);

  // Loop over generators. Use OpenMP parallelisation if enabled.
#ifdef OPENMP
  #pragma omp parallel for
#endif
  for (int iRun=1; iRun<=2; ++iRun) {
    Pythia pythia;
    // Settings common to both runs
    pythia.readString("Beams:eCM = 7000.");
    pythia.readString("HardQCD:all = on");
    pythia.readString("PhaseSpace:pTHatMin = 100.");
    pythia.readString("Next:numberShowInfo = 0");
    pythia.readString("Next:numberShowProcess = 1");
    pythia.readString("Next:numberShowEvent = 1");
    pythia.readString("PartonLevel:MPI = on");
    pythia.readString("HadronLevel:all = on");
    // Settings specific to second run
    if (iRun == 2) {
      // Switch to VINCIA shower model
      pythia.readString("PartonShowers:Model = 2");
    }
    // Initialise generator for this run
    if(!pythia.init()) {continue;}

    // Set histogram pointers
    Hist* nJetsPtr    = &nJetsModel1;
    Hist* eTjetsPtr   = &eTjetsModel1;
    Hist* yJetsPtr    = &yJetsModel1;
    Hist* phiJetsPtr  = &phiJetsModel1;
    Hist* distJetsPtr = &distJetsModel1;
    // Switch to Model2 histograms for second
    if (iRun == 2) {
      nJetsPtr    = &nJetsModel2;
      eTjetsPtr   = &eTjetsModel2;
      yJetsPtr    = &yJetsModel2;
      phiJetsPtr  = &phiJetsModel2;
      distJetsPtr = &distJetsModel2;
    }

    // Set up SlowJet jet finder, with anti-kT clustering, R = 0.7,
    // pT > 20 GeV, |eta| < 4, and pion mass assumed for non-photons
    double etaMax   = 4.;
    double radius   = 0.7;
    double pTjetMin = 20.;
    // Exclude neutrinos (and other invisible) from study.
    int    nSel     = 2;
    SlowJet slowJet( -1, radius, pTjetMin, etaMax, nSel, 1);

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

      // Generate event.
      if (!pythia.next()) {
        cout << " Event generation aborted prematurely, owing to error!\n";
        cout << " Run number was : "<