main442

// main442.cc is a part of the PYTHIA event generator.
// Copyright (C) 2024 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.

// Authors:
//            Christian Bierlich

// Keywords:
//            Hadronization
//            Colour reconnection
//            Rope hadronization

// This program provides a demonstration of the flavour rope model supplied
// in the Rope Hadronization framework. It produces four histograms
// showing the ratio of respectively K^0_s, Lambda_0, Cascade and Omega^-
// to pions as function of event multiplicity.
// No kind of Levy-Tsallis fitting, triggering or similar is done
// is this simple demonstration analysis. It should therefore not
// be taken as anything but a proof of concept.

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

//==========================================================================

int main() {

  // Generator. Process selection. LHC initialization.
  Pythia pythia;
  pythia.readString("Beams:eCM = 7000.");
  pythia.readString("SoftQCD:nonDiffractive = on");

  // Enabling flavour ropes, setting model parameters.
  // The model is still untuned. These parameter values
  // are chosen for illustrative purposes.
  pythia.readString("Ropewalk:RopeHadronization = on");
  pythia.readString("Ropewalk:doShoving = off");
  pythia.readString("Ropewalk:doFlavour = on");
  pythia.readString("Ropewalk:r0 = 0.5");
  pythia.readString("Ropewalk:m0 = 0.2");
  pythia.readString("Ropewalk:beta = 0.1");

  // Enabling setting of vertex information.
  pythia.readString("PartonVertex:setVertex = on");

  // Prevent unstable particles from decaying.
  pythia.readString("ParticleDecays:limitTau0 = on");
  pythia.readString("ParticleDecays:tau0Max = 10");

  // If Pythia fails to initialize, exit with error.
  if (!pythia.init()) return 1;

  // Histograms.
  Hist pion("pions (mult)", 50, 10., 135.);
  Hist kaon("kaons (mult)", 50, 10., 135.);
  Hist lambda("lambdas (mult)", 50, 10., 135.);
  Hist xi("xi (mult)", 50, 10., 135.);
  Hist omega("omega (mult)", 50, 10., 135.);

  // Note: High statistics is needed to fill the high multiplicity end
  // of the histograms, especially for Omega.
  const int nEvent = 4000;

  // Begin event loop. Generate event. Skip if error. List first one.
  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
    if (!pythia.next()) continue;

    // Counters for particle species.
    int nCharged = 0, nPions = 0, nKaons = 0;
    int nLambdas = 0, nXis = 0, nOmegas = 0;
    // Event short notation.
    Event& event = pythia.event;
    for (int i = 0; i < event.size(); ++i){
        Particle& p = event[i];
        // Apply simple, particle level, cuts.
        if (p.isFinal() && abs(p.eta()) < 2.5 && p.pT() > 0.1 ) {
          if(p.isCharged()) ++nCharged;
          int absid = abs(p.id());
          if(absid == 211) ++nPions;
          else if(absid == 310) ++nKaons;
          else if(absid == 3122) ++nLambdas;
          else if(absid == 3312) ++nXis;
          else if(absid == 3334) ++nOmegas;
        }
     }

     // Discard events with event multiplicity less than 10.
     if(nCharged < 10) continue;

     // Fill histograms.
     pion.fill( double(nCharged), double(nPions) );
     kaon.fill( double(nCharged), double(nKaons) );
     lambda.fill( double(nCharged), double(nLambdas) );
     xi.fill( double(nCharged), double(nXis) );
     omega.fill( double(nCharged), double(nOmegas) );

  // End of event loop.
  }

  // Construct ratio histograms.
  Hist kp = kaon / pion;
  kp.title("kaon / pion (multiplicity)");
  Hist lp = lambda / pion;
  lp.title("lambda / pion (multiplicity)");
  Hist xp = xi / pion;
  xp.title("xi / pion (multiplicity)");
  Hist op = omega / pion;
  op.title("omega / pion (multiplicity)");

  // Statistics. Histograms. Done.
  pythia.stat();
  cout << kp << lp << xp << op;
  return 0;
}