Shower and ntuplizer

main42.cc

// /cvmfs/sft.cern.ch/lcg/releases/gcc/12.1.0-57c96/x86_64-centos7/bin/g++ main42.cc -o main42 -I/cvmfs/sft.cern.ch/lcg/releases/MCGenerators/pythia8/310-4b0bd/x86_64-centos7-gcc12-opt/include -O2 -std=c++11 -pedantic -W -Wall -Wshadow -fPIC -pthread  -L/cvmfs/sft.cern.ch/lcg/releases/MCGenerators/pythia8/310-4b0bd/x86_64-centos7-gcc12-opt/lib -Wl,-rpath,/cvmfs/sft.cern.ch/lcg/releases/MCGenerators/pythia8/310-4b0bd/x86_64-centos7-gcc12-opt/lib -lpythia8 -ldl -I/cvmfs/sft.cern.ch/lcg/releases/fastjet/3.4.1-5af57/x86_64-centos7-gcc12-opt/include/ -L/cvmfs/sft.cern.ch/lcg/releases/fastjet/3.4.1-5af57/x86_64-centos7-gcc12-opt/lib/ -Wl,-rpath,/cvmfs/sft.cern.ch/lcg/releases/fastjet/3.4.1-5af57/x86_64-centos7-gcc12-opt/lib/ -lfastjet `root-config --cflags --glibs`

// main42.cc is a part of the PYTHIA event generator.
// Copyright (C) 2023 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: Mikhail Kirsanov <Mikhail.Kirsanov@cern.ch>

// Keywords: basic usage; command file; command line option; hepmc; LHE file

// This program illustrates how a file with HepMC events
// can be generated by Pythia8.
// Input and output files are specified on the command line, e.g. like
// ./main42.exe main42.cmnd hepmcout42.dat > out
// The main program contains no analysis; this is intended to happen later.
// It therefore "never" has to be recompiled to handle different tasks.

// WARNING: typically one needs 25 MB/100 events at the LHC.
// Therefore large event samples may be impractical.

#include "Pythia8/Pythia.h"
//
// #ifndef HEPMC2
// #include "Pythia8Plugins/HepMC3.h"
// #else
// #include "Pythia8Plugins/HepMC2.h"
// #endif

#include "TFile.h"
#include "TTree.h"

#include <fastjet/PseudoJet.hh>
#include <fastjet/ClusterSequence.hh>
#include <fastjet/Selector.hh>

using namespace Pythia8;
using namespace fastjet;


const Int_t MAX_PARTICLES = 3000;
const Int_t MAX_SLIMMED_PARTICLES = 30;
const Int_t MAX_LHE = 20;



int main(int argc, char* argv[]) {

  // Check that correct number of command-line arguments
  if (argc != 3) {
    cerr << " Unexpected number of command-line arguments. \n You are"
         << " expected to provide one input and one output file name. \n"
         << " Program stopped! " << endl;
    return 1;
  }

  double R_jet = 0.4;
  JetDefinition jet_def(antikt_algorithm, R_jet);
  Selector select_akt = SelectorAbsEtaMax(5.0) && SelectorPtMin(15.);

  vector<PseudoJet> hadrons;

  TFile f("tree.root","recreate");
  TTree Events("Events", "Events");
  Float_t LHEPart_pt[MAX_LHE], LHEPart_eta[MAX_LHE], LHEPart_phi[MAX_LHE], LHEPart_mass[MAX_LHE];
  Int_t LHEPart_pdgId[MAX_LHE], LHEPart_status[MAX_LHE];
  Int_t nLHEPart;

  Events.Branch("nLHEPart", &nLHEPart, "nLHEPart/I");
  Events.Branch("LHEPart_pt", LHEPart_pt, "[nLHEPart]/F");
  Events.Branch("LHEPart_eta", LHEPart_eta, "[nLHEPart]/F");
  Events.Branch("LHEPart_phi", LHEPart_phi, "[nLHEPart]/F");
  Events.Branch("LHEPart_mass", LHEPart_mass, "[nLHEPart]/F");
  Events.Branch("LHEPart_pdgId", LHEPart_pdgId, "[nLHEPart]/I");
  Events.Branch("LHEPart_status", LHEPart_status, "[nLHEPart]/I");

  Float_t GenJet_pt[MAX_SLIMMED_PARTICLES], GenJet_eta[MAX_SLIMMED_PARTICLES], GenJet_phi[MAX_SLIMMED_PARTICLES], GenJet_mass[MAX_SLIMMED_PARTICLES];
  Int_t nGenJet;

  Events.Branch("nGenJet", &nGenJet, "nGenJet/I");
  Events.Branch("GenJet_pt", GenJet_pt, "[nGenJet]/F");
  Events.Branch("GenJet_eta", GenJet_eta, "[nGenJet]/F");
  Events.Branch("GenJet_phi", GenJet_phi, "[nGenJet]/F");
  Events.Branch("GenJet_mass", GenJet_mass, "[nGenJet]/F");

  Float_t GenPart_pt[MAX_PARTICLES], GenPart_eta[MAX_PARTICLES], GenPart_phi[MAX_PARTICLES], GenPart_mass[MAX_PARTICLES];
  Int_t GenPart_pdgId[MAX_PARTICLES];
  Int_t nGenPart;

  Events.Branch("nGenPart", &nGenPart, "nGenPart/I");
  Events.Branch("GenPart_pt", GenPart_pt, "[nGenPart]/F");
  Events.Branch("GenPart_eta", GenPart_eta, "[nGenPart]/F");
  Events.Branch("GenPart_phi", GenPart_phi, "[nGenPart]/F");
  Events.Branch("GenPart_mass", GenPart_mass, "[nGenPart]/F");
  Events.Branch("GenPart_pdgId", GenPart_pdgId, "[nGenPart]/F");

  Int_t Map_genJetIdx[MAX_PARTICLES], Map_genPartIdx[MAX_PARTICLES];
  Int_t nMap;

  Events.Branch("nMap", &nMap, "nMap/I");
  Events.Branch("Map_genJetIdx", Map_genJetIdx, "[nMap]/I");
  Events.Branch("Map_genPartIdx", Map_genPartIdx, "[nMap]/I");

  // Check that the provided input name corresponds to an existing file.
  ifstream is(argv[1]);
  if (!is) {
    cerr << " Command-line file " << argv[1] << " was not found. \n"
         << " Program stopped! " << endl;
    return 1;
  }

  // Confirm that external files will be used for input and output.
  cout << "\n >>> PYTHIA settings will be read from file " << argv[1]
       << " <<< \n >>> HepMC events will be written to file "
       << argv[2] << " <<< \n" << endl;

  // Interface for conversion from Pythia8::Event to HepMC event.
  // Specify file where HepMC events will be stored.
  // Pythia8ToHepMC toHepMC(argv[2]);

  // Generator.
  Pythia pythia;

  // Read in commands from external file.
  pythia.readFile(argv[1]);

  // Extract settings to be used in the main program.
  int    nEvent    = pythia.mode("Main:numberOfEvents");
  int    nAbort    = pythia.mode("Main:timesAllowErrors");

  // Initialization.
  pythia.init();

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

    // Generate event.
    if (!pythia.next()) {

      // If failure because reached end of file then exit event loop.
      if (pythia.info.atEndOfFile()) {
        cout << " Aborted since reached end of Les Houches Event File\n";
        break;
      }

      // First few failures write off as "acceptable" errors, then quit.
      if (++iAbort < nAbort) continue;
      cout << " Event generation aborted prematurely, owing to error!\n";
      break;
    }

    // Construct new empty HepMC event, fill it and write it out.
    // toHepMC.writeNextEvent( pythia );
      int pid, status;
      double px, py, pz, e, mass;
      // double x, y, z, t;

      hadrons.clear () ; 

      nLHEPart = 0;
      nGenPart = 0;
      nGenJet = 0;
      nMap = 0;

      for(uint i = 1; i < pythia.event.size(); ++i)
      {
        Pythia8::Particle &particle = pythia.event[i];
      status = particle.status();
      if ((abs(status) >= 21) && (abs(status) <= 29)){
        // std::cout << particle.index() << " " << status << " " << particle.id() << "\n";
        LHEPart_pt[nLHEPart] = particle.pT();
        LHEPart_eta[nLHEPart] = particle.eta();
        LHEPart_phi[nLHEPart] = particle.phi();
        LHEPart_mass[nLHEPart] = particle.m();
        LHEPart_pdgId[nLHEPart] = particle.id();
        LHEPart_status[nLHEPart] = abs(status)-21;
        nLHEPart++;
      }

        status = particle.statusHepMC();
        if (particle.statusHepMC()!=1 || !particle.isFinal()){
          continue;
        }
      // if (particle.isHadron()){
      if (particle.isVisible()){
        GenPart_pt[nGenPart] = particle.pT();
        GenPart_eta[nGenPart] = particle.eta();
        GenPart_phi[nGenPart] = particle.phi();
        GenPart_mass[nGenPart] = particle.m();
        GenPart_pdgId[nGenPart] = particle.id();

        px = particle.px();
        py = particle.py();
        pz = particle.pz();
        e = particle.e();
        PseudoJet tmp (px, py, pz, e);
        tmp.set_user_index(nGenPart);
        hadrons.push_back(tmp);

        nGenPart++;
      }
    }

    ClusterSequence cluster (hadrons, jet_def) ;
    vector<PseudoJet> antikT_jets = sorted_by_pt (select_akt (cluster.inclusive_jets ())) ;
    // std::cout << "    +-> clustered jet number: " << antikT_jets.size () << "\n" ;
    // // print the jets
    // std::cout <<   "        pt eta phi mass" << endl;
    for (unsigned i = 0; i < antikT_jets.size(); i++) {
      GenJet_pt[nGenJet] = antikT_jets[i].pt();
      GenJet_eta[nGenJet] = antikT_jets[i].eta();
      GenJet_phi[nGenJet] = antikT_jets[i].phi();
      GenJet_mass[nGenJet] = antikT_jets[i].m();
      vector<PseudoJet> constituents = antikT_jets[i].constituents();
      for (unsigned j = 0; j < constituents.size(); j++) {
        // cout  << nGenJet << " " constituents[j].user_index()
        //      << endl;
        Map_genJetIdx[nMap] = nGenJet;
        Map_genPartIdx[nMap] = constituents[j].user_index();
        nMap ++;
      }
      nGenJet++;
  }



    Events.Fill();
  // End of event loop. Statistics.
  }
  Events.Write();
  pythia.stat();

  // Done.
  return 0;
}

main42.cmnd

! File: main42.cmnd
! This file contains commands to be read in for a Pythia8 run.
! Lines not beginning with a letter or digit are comments.
! Names are case-insensitive  -  but spellings-sensitive!
! The changes here are illustrative, not always physics-motivated.

! 1) Settings that will be used in a main program.
Main:numberOfEvents = 10         ! number of events to generate
Main:timesAllowErrors = 3          ! abort run after this many flawed events

! 2) Settings related to output in init(), next() and stat().
Init:showChangedSettings = on      ! list changed settings
Init:showAllSettings = off         ! list all settings
Init:showChangedParticleData = on  ! list changed particle data
Init:showAllParticleData = off     ! list all particle data
Next:numberCount = 1000            ! print message every n events
Next:numberShowLHA = 1             ! print LHA information n times
Next:numberShowInfo = 1            ! print event information n times
Next:numberShowProcess = 1         ! print process record n times
Next:numberShowEvent = 1           ! print event record n times
Stat:showPartonLevel = on          ! additional statistics on MPI

! 3) Beam parameter settings. Values below agree with default ones.
Beams:idA = 2212                   ! first beam, p = 2212, pbar = -2212
Beams:idB = 2212                   ! second beam, p = 2212, pbar = -2212
Beams:eCM = 14000.                 ! CM energy of collision

! 4) PDF settings. Default is to use internal PDFs
! some pdf sets examples: cteq61.LHpdf cteq61.LHgrid MRST2004nlo.LHgrid
#PDF:pSet = LHAPDF5:MRST2001lo.LHgrid
! Allow extrapolation of PDF's beyond x and Q2 boundaries, at own risk.
! Default behaviour is to freeze PDF's at boundaries.
#PDF:extrapolate = on

! 5a) Pick processes and kinematics cuts.
! Colour singlet charmonium production of J/psi and chi_c.
#Charmonium:gg2ccbar(3S1)[3S1(1)]g    = on,off
#Charmonium:gg2ccbar(3PJ)[3PJ(1)]g    = on,on,on
#Charmonium:qg2ccbar(3PJ)[3PJ(1)]q    = on,on,on
#Charmonium:qqbar2ccbar(3PJ)[3PJ(1)]g = on,on,on
#PhaseSpace:pTHatMin = 20.          ! minimum pT of hard process

! 5b) Alternative beam and process selection from a Les Houches Event File.
! NOTE: to use this option, comment out the lines in section 5a above
! and uncomment the ones below. Section 3 is ignored for frameType = 4.
Beams:frameType = 4                ! read info from a LHEF
Beams:LHEF = events.lhe             ! the LHEF to read from

! 6) Other settings. Can be expanded as desired.
! Note: may overwrite some of the values above, so watch out.
#Tune:pp = 6                        ! use Tune 4Cx
#ParticleDecays:limitTau0 = on      ! set long-lived particle stable ...
#ParticleDecays:tau0Max = 10        ! ... if c*tau0 > 10 mm