Pythia8Rivet.h

// Pythia8Rivet.h 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.

#ifndef PYTHIA8RIVET_H
#define PYTHIA8RIVET_H

#include "Pythia8/HIUserHooks.h"
#include "Pythia8Plugins/HepMC2.h"
#include "HepMC/GenEvent.h"
#include "HepMC/IO_GenEvent.h"
#include <set>
#include <vector>
#include <string>
#include "Rivet/Rivet.hh"

namespace Pythia8 {

class Pythia8Rivet {

public:

  Pythia8Rivet(Pythia & pytin, string fname)
    : pythia(&pytin), filename(fname), rivet(0), igBeam(false) {}

  ~Pythia8Rivet() {
    done();
  }

  void addAnalysis(string ana) {
    analyses.insert(ana);
  }

  void addPreload(string prel) {
    preloads.push_back(prel);
  }

  void ignoreBeams(bool flagIn) {
    igBeam = flagIn;
  }

  void addRunName(const string runname) {
    rname = runname;
  }
  void init(const HepMC::GenEvent & gev) {
    if ( rivet ) return;
    rivet = new Rivet::AnalysisHandler(rname);
    rivet->setIgnoreBeams(igBeam);
    Rivet::addAnalysisLibPath(".");
    for(int i = 0, N = preloads.size(); i < N; ++i)
      rivet->readData(preloads[i]);
    for (set<string>::iterator it = analyses.begin();
      it != analyses.end(); ++it) {
      rivet->addAnalysis(*it);
    }
    rivet->init(gev);
  }

  void operator()() {
    this->operator()(pythia->event, -1, &pythia->info, &pythia->settings);
  }

  void operator()(Event & event, int ievnum = -1,
                  const Pythia8::Info* pyinfo = 0,
                  Pythia8::Settings* pyset = 0, bool append = false,
                  HepMC::GenParticle* rootParticle = 0, int iBarcode = -1) {
    HepMC::GenEvent geneve;
    converter.fill_next_event(event, &geneve, ievnum, pyinfo, pyset,
                              append, rootParticle, iBarcode);
    if ( pyinfo && pyinfo->hiInfo ) {
      HepMC::HeavyIon ion;
      ion.set_Ncoll_hard(pyinfo->hiInfo->nCollNDTot());
      ion.set_Ncoll(pyinfo->hiInfo->nAbsProj() +
                    pyinfo->hiInfo->nDiffProj() +
                    pyinfo->hiInfo->nAbsTarg() +
                    pyinfo->hiInfo->nDiffTarg() -
                    pyinfo->hiInfo->nCollND() -
                    pyinfo->hiInfo->nCollDD());
      ion.set_Npart_proj(pyinfo->hiInfo->nAbsProj() +
                         pyinfo->hiInfo->nDiffProj());
      ion.set_Npart_targ(pyinfo->hiInfo->nAbsTarg() +
                         pyinfo->hiInfo->nDiffTarg());
      ion.set_impact_parameter(pyinfo->hiInfo->b());
      geneve.set_heavy_ion(ion);
    }
    if ( !rivet ) init(geneve);
    rivet->analyze(geneve);
  }

  void done() {
    if ( !rivet ) return;
    rivet->finalize();
    rivet->writeData(filename);
    delete rivet;
    rivet = 0;
  }

private:

  Pythia * pythia;

  string filename;

  set<string> analyses;

  vector<string> preloads;

  Rivet::AnalysisHandler * rivet;

  HepMC::Pythia8ToHepMC converter;

  string rname;

  bool igBeam;

};

}

#endif