EvtGenDecays.h

// This file is a copy of Pythia8Plugins/EvtGen.h from Pythia 8,
// but slightly modified in basf2 for extended EvtGen models
 
// EvtGen.h is a part of the PYTHIA event generator.
// Copyright (C) 2016 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.
// Author: Philip Ilten.
 
// This file contains an EvtGen interface. HepMC and EvtGen must be enabled.
 
#pragma once
 
#include <generators/utilities/GeneratorConst.h>
#include "Pythia8/Pythia.h"
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtRandomEngine.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include "EvtGenBase/EvtParticleDecayList.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#include "EvtGenExternal/EvtExternalGenList.hh"
 
//==========================================================================
 
class EvtGenRandom : public EvtRandomEngine {
 
public:
 
  explicit EvtGenRandom(Pythia8::Rndm* rndmPtrIn) {rndmPtr = rndmPtrIn;}
 
  double random() {if (rndmPtr) return rndmPtr->flat(); else return -1.0;}
 
  Pythia8::Rndm* rndmPtr;
 
};
 
//==========================================================================
 
class EvtGenDecays {
 
public:
 
  EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, EvtGen* evtGen, bool limit = true);
 
  EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, std::string decayFile, std::string particleDataFile,
               EvtExternalGenList* extPtrIn = 0, EvtAbsRadCorr* fsrPtrIn = 0,
               int mixing = 1, bool xml = false, bool limit = true,
               bool extUse = true, bool fsrUse = true);
 
  ~EvtGenDecays()
  {
    if (evtgen) delete evtgen;
    if (extOwner && extPtr) delete extPtr;
    if (fsrOwner && fsrPtr) delete fsrPtr;
  }
 
  EvtGenDecays(const EvtGenDecays&) = delete;
 
  EvtGenDecays& operator=(const EvtGenDecays&) = delete;
 
  void getDecayLimits(bool limit);
 
  double decay();
 
  void exclude(int id) {excIds.insert(id);}
 
  void updatePythia();
 
  void updateEvtGen();
 
  void readDecayFile(std::string decayFile, bool xml = false)
  {
    evtgen->readUDecay(decayFile.c_str(), xml); updateData();
  }
 
  bool extOwner; 
  bool fsrOwner; 
  EvtExternalGenList* extPtr;      
  EvtAbsRadCorr*      fsrPtr;      
  std::list<EvtDecayBase*> models; 
  struct Signal {
    int status;                 
    EvtId egId;                 
    std::vector<double> bfs;    
    std::vector<int> map;       
    EvtParticleDecayList modes; 
  };
  std::map<int, Signal> signals;
 
  std::string signalSuffix;
 
protected:
 
  void updateData(bool final = false);
 
  void updateEvent(Pythia8::Particle* pyPro, EvtParticle* egPro,
                   std::vector<int>* pySigs = 0, std::vector<EvtParticle*>* egSigs = 0,
                   std::vector<double>* bfs = 0, double* wgt = 0);
 
  bool checkVertex(Pythia8::Particle* pyPro);
 
  bool checkSignal(Pythia8::Particle* pyPro);
 
  bool checkOsc(EvtParticle* egPro);
 
  static constexpr int NTRYDECAY = 1000;
 
  Pythia8::Pythia* pythiaPtr;
 
  EvtGenRandom rndm;
 
  EvtGen* evtgen;
 
  std::set<int> incIds; 
  std::set<int> excIds; 
  bool updated;
 
  std::map<int, Signal>::iterator signal;
 
  double tau0Max;     
  double tauMax;      
  double rMax;        
  double xyMax;       
  double zMax;        
  bool limitTau0;     
  bool limitTau;      
  bool limitRadius;   
  bool limitCylinder; 
  bool limitDecay;    
};
 
//--------------------------------------------------------------------------
 
// The constructor.
 
// The EvtGenDecays object is associated with a single Pythia
// instance. This is to ensure a consistent random number generator
// across the two, as well as any updates to particle data, etc. Note
// that if multiple EvtGenDecays objects exist, that they will modify
// one anothers particle databases due to the design of EvtGen.
 
// This constructor also sets all particles to be decayed by EvtGen as
// stable within Pythia. The parameters within Pythia used to check if
// a particle should be decayed, as described in the "Particle Decays"
// section of the Pythia manual, are set. Note that if the variable
// "limit" is set to "false", then no check will be made before
// decaying a particle with EvtGen.
 
// The constructor is designed to have the exact same form as the
// EvtGen constructor except for these five differences.
// (1) The first variable is the pointer to the Pythia object.
// (2) The third last argument is a flag to limit decays based on the
//     Pythia criteria (based on the particle decay vertex).
// (3) The second last argument is a flag if external models should be
//     passed to EvtGen (default is true).
// (4) The last argument is a flag if an FSR model should be passed
//     to EvtGen (default is true).
// (5) No random engine pointer is passed, as this is obtained from
//     Pythia.
 
//   pythiaPtrIn:      the pointer to the associated Pythia object.
//   decayFile:        the name of the decay file to pass to EvtGen.
//   particleDataFile: the name of the particle data file to pass to EvtGen.
//   extPtrIn:         the optional EvtExternalGenList pointer, this must be
//                     be provided if fsrPtrIn is provided to avoid double
//                     initializations.
//   fsrPtrIn:         the EvtAbsRadCorr pointer to pass to EvtGen.
//   mixing:           the mixing type to pass to EvtGen.
//   xml:              flag to use XML files to pass to EvtGen.
//   limit:            flag to limit particle decays based on Pythia criteria.
//   extUse:           flag to use external models with EvtGen.
//   fsrUse:           flag to use radiative correction engine with EvtGen.
 
 
EvtGenDecays::EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, EvtGen* evtGen, bool limit):
  extOwner(false), fsrOwner(false), extPtr(nullptr), fsrPtr(nullptr),
  signalSuffix("_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(nullptr),
  evtgen(evtGen), updated(false)
{
  getDecayLimits(limit);
}
 
EvtGenDecays::EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, std::string decayFile,
                           std::string particleDataFile, EvtExternalGenList* extPtrIn,
                           EvtAbsRadCorr* fsrPtrIn, int mixing, bool xml, bool limit,
                           bool extUse, bool fsrUse) : extPtr(extPtrIn), fsrPtr(fsrPtrIn),
  signalSuffix("_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(&pythiaPtr->rndm),
  updated(false)
{
 
  // Initialize EvtGen.
  if (!extPtr && fsrPtr) {
    B2ERROR("Error in EvtGenDecays::EvtGenDecays: extPtr is null but fsrPtr is provided");
    return;
  }
  if (extPtr) extOwner = false;
  else {extOwner = true; extPtr = new EvtExternalGenList();}
  if (fsrPtr) fsrOwner = false;
  else {fsrOwner = true; fsrPtr = extPtr->getPhotosModel();}
  models = extPtr->getListOfModels();
  evtgen = new EvtGen(decayFile.c_str(), particleDataFile.c_str(),
                      &rndm, fsrUse ? fsrPtr : 0, extUse ? &models : 0, mixing, xml);
  getDecayLimits(limit);
}
 
void EvtGenDecays::getDecayLimits(bool limit)
{
  // Get the Pythia decay limits.
  if (!pythiaPtr) return;
  limitTau0     = pythiaPtr->settings.flag("ParticleDecays:limitTau0");
  tau0Max       = pythiaPtr->settings.parm("ParticleDecays:tau0Max");
  limitTau      = pythiaPtr->settings.flag("ParticleDecays:limitTau");
  tauMax        = pythiaPtr->settings.parm("ParticleDecays:tauMax");
  limitRadius   = pythiaPtr->settings.flag("ParticleDecays:limitRadius");
  rMax          = pythiaPtr->settings.parm("ParticleDecays:rMax");
  limitCylinder = pythiaPtr->settings.flag("ParticleDecays:limitCylinder");
  xyMax         = pythiaPtr->settings.parm("ParticleDecays:xyMax");
  zMax          = pythiaPtr->settings.parm("ParticleDecays:zMax");
  limitDecay    = limit && (limitTau0 || limitTau ||
                            limitRadius || limitCylinder);
 
}
 
//--------------------------------------------------------------------------
 
// Perform all decays and return the event weight.
 
// All particles in the event record that can be decayed by EvtGen are
// decayed. If a particle is a signal particle, then this is stored in
// a vector of signal particles. A signal particle is only stored if
// its status is the same as the status provided in the signals map. A
// negative status in the signal map indicates that all statuses
// should be accepted. After all signal particles are identified, one
// is randomly chosen and decayed as signal. The remainder are decayed
// normally.
 
// Forcing a signal decay changes the weight of an event from unity,
// and so the relative event weight is returned, given the forced
// signal decay. A weight of 0 indicates no signal in the event, while
// a weight of -1 indicates something is wrong, e.g. either the Pythia
// or EvtGen pointers are not available or the number of tries has
// been exceeded. For the event weight to be valid, one should not
// change the absolute branching fractions in the signal and inclusive
// definitions, but rather just remove the unwanted decay channels
// from the signal decay definition.
 
double EvtGenDecays::decay()
{
 
  // Reset the signal and signal counters.
  if (!pythiaPtr || !evtgen) return -1;
  if (!updated) updateData(true);
 
  // Loop over all particles in the Pythia event.
  Pythia8::Event& event = pythiaPtr->event;
  std::vector<int> pySigs; std::vector<EvtParticle*> egSigs, egPrts;
  std::vector<double> bfs; double wgt(1.);
  for (int iPro = 0; iPro < event.size(); ++iPro) {
 
    // Check particle is final and can be decayed by EvtGen.
    Pythia8::Particle* pyPro = &event[iPro];
    if (!pyPro->isFinal()) continue;
    if (incIds.find(pyPro->id()) == incIds.end()) continue;
    if (pyPro->status() == 93 || pyPro->status() == 94) continue;
 
    // Decay the progenitor with EvtGen.
    EvtParticle* egPro = EvtParticleFactory::particleFactory
                         (EvtPDL::evtIdFromStdHep(pyPro->id()),
                          EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz()));
    egPrts.push_back(egPro);
    egPro->setDiagonalSpinDensity();
    evtgen->generateDecay(egPro);
    pyPro->tau(egPro->getLifetime());
    if (!checkVertex(pyPro)) continue;
 
    // Add oscillations to event record.
    if (checkOsc(egPro))
      updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
 
    // Undo decay if signal (duplicate to stop oscillations).
    else if (checkSignal(pyPro)) {
      pySigs.push_back(pyPro->index());
      egSigs.push_back(egPro);
      bfs.push_back(signal->second.bfs[0]);
      wgt *= 1 - bfs.back();
      egPro->deleteDaughters();
      // The following problematic part of the Pythia's code was removed, see
      // https://gitlab.com/Pythia8/releases/
      // -/blob/pythia8310/include/Pythia8Plugins/EvtGen.h?ref_type=tags#L321-325
 
      // If not signal, add to event record.
    } else updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
  }
  if (pySigs.size() == 0) {
    for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
      egPrts[iPrt]->deleteTree();
    return 0;
  }
 
  // Determine the decays of the signal particles (signal or background).
  std::vector<int> modes; int force(-1);
  for (int iTry = 1; iTry <= NTRYDECAY; ++iTry) {
    int n(0);
 
    modes.clear(); force = pythiaPtr->rndm.pick(bfs); n = 0;
    for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
      if (iSig == force) modes.push_back(0);
      else modes.push_back(pythiaPtr->rndm.flat() > bfs[iSig]);
      if (modes.back() == 0) ++n;
    }
    if (pythiaPtr->rndm.flat() <= 1.0 / n) break;
    if (iTry == NTRYDECAY) {
      wgt = 2.;
      B2WARNING("Warning in EvtGenDecays::decay: maximum number of signal decay attempts exceeded"
                << LogVar("maxTries", NTRYDECAY));
    }
  }
 
  // Decay the signal particles and mark forced decay.
  for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
    EvtParticle* egSig = egSigs[iSig];
    Pythia8::Particle*    pySig = &event[pySigs[iSig]];
    signal = signals.find(pySig->id());
    if (egSig->getNDaug() > 0) egSig = egSig->getDaug(0);
    if (modes[iSig] == 0) egSig->setId(signal->second.egId);
    else {
      signal->second.modes.getDecayModel(egSig);
      egSig->setChannel(signal->second.map[egSig->getChannel()]);
    }
    if (iSig == force)
      pySig->status(pySig->status() == 92 || pySig->status() == 94 ? 96 : 95);
    evtgen->generateDecay(egSig);
    updateEvent(pySig, egSigs[iSig]);
  }
 
  // Delete all EvtGen particles and return weight.
  for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
    egPrts[iPrt]->deleteTree();
  return 1. - wgt;
 
}
 
//--------------------------------------------------------------------------
 
// Update the Pythia particle database from EvtGen.
 
// Note that only the particle spin type, charge type, nominal mass,
// width, minimum mass, maximum mass, and nominal lifetime are
// set. The name string is not set.
 
void EvtGenDecays::updatePythia()
{
  if (!pythiaPtr || !evtgen) return;
  for (int entry = 0; entry < (int)EvtPDL::entries(); ++entry) {
    EvtId egId = EvtPDL::getEntry(entry);
    int   pyId = EvtPDL::getStdHep(egId);
    pythiaPtr->particleData.spinType(pyId, EvtPDL::getSpinType(egId));
    pythiaPtr->particleData.chargeType(pyId, EvtPDL::chg3(egId));
    pythiaPtr->particleData.m0(pyId, EvtPDL::getMass(egId));
    pythiaPtr->particleData.mWidth(pyId, EvtPDL::getWidth(egId));
    pythiaPtr->particleData.mMin(pyId, EvtPDL::getMinMass(egId));
    pythiaPtr->particleData.mMax(pyId, EvtPDL::getMaxMass(egId));
    pythiaPtr->particleData.tau0(pyId, EvtPDL::getctau(egId));
  }
}
 
//--------------------------------------------------------------------------
 
// Update the EvtGen particle database from Pythia.
 
// The particle update database between Pythia and EvtGen is not
// symmetric. Particularly, it is not possible to update the spin
// type, charge, or nominal lifetime in the EvtGen particle database.
 
void EvtGenDecays::updateEvtGen()
{
  if (!pythiaPtr || !evtgen) return;
  int pyId = pythiaPtr->particleData.nextId(0);
  while (pyId != 0) {
    EvtId egId = EvtPDL::evtIdFromStdHep(pyId);
    EvtPDL::reSetMass(egId, pythiaPtr->particleData.m0(pyId));
    EvtPDL::reSetWidth(egId, pythiaPtr->particleData.mWidth(pyId));
    EvtPDL::reSetMassMin(egId, pythiaPtr->particleData.mMin(pyId));
    EvtPDL::reSetMassMax(egId, pythiaPtr->particleData.mMax(pyId));
    pyId = pythiaPtr->particleData.nextId(pyId);
  }
}
 
//--------------------------------------------------------------------------
 
// Update the particles to decay with EvtGen, and the selected signals.
 
// If final is false, then only signals are initialized in the signal
// map. Any particle or alias that ends with signalSuffix is taken as
// a signal particle. If final is true all particle entries in EvtGen
// are checked to see if they should be set stable in Pythia. If an
// EvtGen particle has no decay modes, then Pythia is still allowed to
// decay the particle. Additionally, the signal decay channels are
// turned off for the non-aliased signal particle.
 
void EvtGenDecays::updateData(bool final)
{
 
  // Loop over the EvtGen entries.
  if (!pythiaPtr) return;
  EvtDecayTable* egTable = EvtDecayTable::getInstance();
  if (!egTable) return;
  for (int iEntry = 0; iEntry < (int)EvtPDL::entries(); ++iEntry) {
    EvtId egId = EvtPDL::getEntry(iEntry);
    int   pyId = EvtPDL::getStdHep(egId);
    if (egTable->getNModes(egId) == 0) continue;
    if (excIds.find(pyId) != excIds.end()) continue;
 
    // Stop Pythia from decaying the particle and include in decay set.
    if (final)  {
      incIds.insert(pyId);
      pythiaPtr->particleData.mayDecay(pyId, false);
    }
 
    // Check for signal.
    std::string egName = EvtPDL::name(egId);
    if (egName.size() <= signalSuffix.size() || egName.substr
        (egName.size() - signalSuffix.size()) != signalSuffix) continue;
    signal = signals.find(pyId);
    if (signal == signals.end()) {
      signals[pyId].status = -1;
      signal = signals.find(pyId);
    }
    signal->second.egId  = egId;
    signal->second.bfs   = std::vector<double>(2, 0);
    if (!final) continue;
 
    // Get the signal and background decay modes.
    std::vector<EvtParticleDecayList> egList = egTable->getDecayTable();
    int sigIdx = egId.getAlias();
    int bkgIdx = EvtPDL::evtIdFromStdHep(pyId).getAlias();
    if (sigIdx > (int)egList.size() || bkgIdx > (int)egList.size()) continue;
    EvtParticleDecayList sigModes = egList[sigIdx];
    EvtParticleDecayList bkgModes = egList[bkgIdx];
    EvtParticleDecayList allModes = egList[bkgIdx];
    double sum(0);
 
    // Sum signal branching fractions.
    for (int iMode = 0; iMode < sigModes.getNMode(); ++iMode) {
      EvtDecayBase* mode = sigModes.getDecayModel(iMode);
      if (!mode) continue;
      signal->second.bfs[0] += mode->getBranchingFraction();
      sum += mode->getBranchingFraction();
    }
 
    // Sum remaining background branching fractions.
    for (int iMode = 0; iMode < allModes.getNMode(); ++iMode) {
      EvtDecayBase* mode = allModes.getDecayModel(iMode);
      if (!mode) continue;
      bool match(false);
      for (int jMode = 0; jMode < sigModes.getNMode(); ++jMode)
        if (mode->matchingDecay(*(sigModes.getDecayModel(jMode)))) {
          match = true; break;
        }
      if (match) bkgModes.removeMode(mode);
      else {
        signal->second.map.push_back(iMode);
        signal->second.bfs[1] += mode->getBranchingFraction();
        sum += mode->getBranchingFraction();
      }
    }
    signal->second.modes = bkgModes;
    for (int iBf = 0; iBf < (int)signal->second.bfs.size(); ++iBf)
      signal->second.bfs[iBf] /= sum;
  }
  if (final) updated = true;
 
}
 
//--------------------------------------------------------------------------
 
// Update the Pythia event record with an EvtGen decay tree.
 
// The production vertex of each particle (which can also be obtained
// in EvtGen via EvtParticle::get4Pos()) is set by the decay vertex of
// its mother, which in turn is calculated from the mother's
// lifetime. The status code 93 is used to indicate an external decay,
// while the status code 94 is used to indicate an oscillated external
// decay. If the progenitor has a single daughter with the same ID,
// this daughter is used as the progenitor. This is used to prevent
// double oscillations.
 
// If the arguments after egPro are no NULL and a particle in the
// decay tree is a signal particle, the decay for this particle is
// removed and the particle is stored as a signal candidate in the
// pySigs and egSigs vectors, to be decayed later. However, if any of
// these arguments is NULL then the entire tree is written.
 
void EvtGenDecays::updateEvent(Pythia8::Particle* pyPro, EvtParticle* egPro,
                               std::vector<int>* pySigs, std::vector<EvtParticle*>* egSigs,
                               std::vector<double>* bfs, double* wgt)
{
 
  // Set up the mother vector.
  if (!pythiaPtr) return;
  EvtParticle* egMom = egPro;
  if (egPro->getNDaug() == 1 && egPro->getPDGId() ==
      egPro->getDaug(0)->getPDGId()) egMom = egPro->getDaug(0);
  Pythia8::Event& event = pythiaPtr->event;
  std::vector< std::pair<EvtParticle*, int> >
  moms(1, std::pair<EvtParticle*, int>(egMom, pyPro->index()));
 
  // Loop over the mothers.
  while (moms.size() != 0) {
 
    // Check if particle can decay.
    egMom = moms.back().first;
    int       iMom  = moms.back().second;
    Pythia8::Particle* pyMom = &event[iMom];
    moms.pop_back();
    if (!checkVertex(pyMom)) continue;
    bool osc(checkOsc(egMom));
 
    // Set the children of the mother.
    pyMom->daughters(event.size(), event.size() + egMom->getNDaug() - 1);
    pyMom->statusNeg();
    Pythia8::Vec4 vProd = pyMom->vDec();
    for (int iDtr = 0 ; iDtr < (int)egMom->getNDaug(); ++iDtr) {
      EvtParticle* egDtr = egMom->getDaug(iDtr);
      int          id    = egDtr->getPDGId();
      EvtVector4R  p     = egDtr->getP4Lab();
      int idx = event.append(id, 93, iMom, 0, 0, 0, 0, 0, p.get(1),
                             p.get(2), p.get(3), p.get(0), egDtr->mass());
      Pythia8::Particle* pyDtr = &event.back();
      pyDtr->vProd(vProd);
      pyDtr->tau(egDtr->getLifetime());
      if (pySigs && egSigs && bfs && wgt && checkSignal(pyDtr)) {
        pySigs->push_back(pyDtr->index());
        egSigs->push_back(egDtr);
        bfs->push_back(signal->second.bfs[0]);
        (*wgt) *= 1 - bfs->back();
        egDtr->deleteDaughters();
      }
      if (osc) pyDtr->status(94);
      if (egDtr->getAttribute("FSR"))
        pyDtr->status(Belle2::GeneratorConst::FSR_STATUS_CODE);
      if (egDtr->getNDaug() > 0)
        moms.push_back(std::pair<EvtParticle*, int>(egDtr, idx));
    }
  }
 
}
 
//--------------------------------------------------------------------------
 
// Check if a particle can decay.
 
// Modified slightly from ParticleDecays::checkVertex.
 
bool EvtGenDecays::checkVertex(Pythia8::Particle* pyPro)
{
  using ::Pythia8::pow2;
 
  if (!limitDecay) return true;
  if (limitTau0 && pyPro->tau0() > tau0Max) return false;
  if (limitTau && pyPro->tau() > tauMax) return false;
  if (limitRadius && pow2(pyPro->xDec()) + pow2(pyPro->yDec())
      + pow2(pyPro->zDec()) > pow2(rMax)) return false;
  if (limitCylinder && (pow2(pyPro->xDec()) + pow2(pyPro->yDec())
                        > pow2(xyMax) || abs(pyPro->zDec()) > zMax)) return false;
  return true;
}
 
//--------------------------------------------------------------------------
 
// Check if a particle is signal.
 
bool EvtGenDecays::checkSignal(Pythia8::Particle* pyPro)
{
  signal = signals.find(pyPro->id());
  if (signal != signals.end() && (signal->second.status < 0 ||
                                  signal->second.status == pyPro->status())) return true;
  else return false;
}
 
//--------------------------------------------------------------------------
 
// Check if an EvtGen particle has oscillated.
 
// The criteria defined here for oscillation is a single daughter but
// with a different ID from the mother.
 
bool EvtGenDecays::checkOsc(EvtParticle* egPro)
{
  if (egPro && egPro->getNDaug() == 1 &&
      egPro->getPDGId() != egPro->getDaug(0)->getPDGId()) return true;
  else return false;
}
 
//==========================================================================