Pythia8Rivet.h

// Pythia8Rivet.h is a part of the PYTHIA event generator.
// Copyright (C) 2018 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 {

using namespace std;

class Pythia8Rivet {

public:

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

  ~Pythia8Rivet() {
    done();
  }

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

  void addRunName(const string runname) {
    rname = runname;
  }
  void init(const HepMC::GenEvent & gev) {
    if ( rivet ) return;
    rivet = new Rivet::AnalysisHandler(rname);
    Rivet::addAnalysisLibPath(".");
    vector<string> anals(analyses.begin(), analyses.end());
    rivet->addAnalyses(anals);
    rivet->init(gev);
  }

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

  void operator()(Event & event, int ievnum = -1, 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;

  Rivet::AnalysisHandler * rivet;

  HepMC::Pythia8ToHepMC converter;

  string rname;

};

}

#endif