DipoleEventRecord.h

// -*- C++ -*-
//
// DipoleEventRecord.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DipoleEventRecord_H
#define HERWIG_DipoleEventRecord_H
//
// This is the declaration of the DipoleEventRecord class.
//
 
#include "Herwig/Shower/ShowerEventRecord.h"
#include "Herwig/Shower/PerturbativeProcess.h"
#include "ThePEG/PDF/PDF.h"
#include "Dipole.h"
#include "DipoleChain.h"
#include "Herwig/MatrixElement/Matchbox/Utility/DensityOperator.h"
 
#include <tuple>
 
namespace Herwig {
 
using namespace ThePEG;
 
class SubleadingSplittingInfo
  : public DipoleSplittingInfo {
 
public:
 
  SubleadingSplittingInfo()
    : DipoleSplittingInfo() {}
 
  list<DipoleChain>::iterator emitterChain() const { return theEmitterChain; }
 
  list<Dipole>::iterator emitterDipole() const { return theEmitterDipole; }
 
  list<DipoleChain>::iterator spectatorChain() const { return theSpectatorChain; }
 
  list<Dipole>::iterator spectatorDipole() const { return theSpectatorDipole; }
 
  Energy startScale() const { return theStartScale; }
  
  void emitterChain(list<DipoleChain>::iterator it) { theEmitterChain = it; }
 
  void emitterDipole(list<Dipole>::iterator it) { theEmitterDipole = it; }
 
  void spectatorChain(list<DipoleChain>::iterator it) { theSpectatorChain = it; }
 
  void spectatorDipole(list<Dipole>::iterator it) { theSpectatorDipole = it; }
 
  void startScale(Energy s) { theStartScale = s; }
 
private:
 
  list<DipoleChain>::iterator theEmitterChain;
 
  list<Dipole>::iterator theEmitterDipole;
 
  list<DipoleChain>::iterator theSpectatorChain;
 
  list<Dipole>::iterator theSpectatorDipole;
 
  Energy theStartScale;
 
};
 
class DipoleEventRecord : public ShowerEventRecord {
 
public:
 
  DipoleEventRecord() {}
 
  ~DipoleEventRecord() { clear(); }
 
public:
 
  PList& hard() { return theHard; }
 
  const PList& hard() const { return theHard; }
 
  const Lorentz5Momentum& pX() const { return thePX; }
 
  cPDVector& particlesAfter() { return theParticlesAfter; }
  
  const cPDVector& particlesAfter() const { return theParticlesAfter; }
  
  cPDVector& particlesBefore() { return theParticlesBefore; }
  
  const cPDVector& particlesBefore() const { return theParticlesBefore; }
 
  vector<Lorentz5Momentum>& momentaAfter() { return theMomentaAfter; }
  
  const vector<Lorentz5Momentum>& momentaAfter() const { return theMomentaAfter; }
 
  map<PPtr,size_t>& particleIndices() { return theParticleIndices; }
 
  const map<PPtr,size_t>& particleIndices() const { return theParticleIndices; }
 
  DensityOperator& densityOperator() { return theDensityOperator; }
 
  const DensityOperator& densityOperator() const { return theDensityOperator; }
  
  void setSubleadingNc( bool doSub, size_t emissionsLimit ) { 
    doSubleadingNc = doSub;
    continueSubleadingNc = doSub;
    subleadingNcEmissionsLimit = emissionsLimit;
  }
  
  bool getContinueSubleadingNc() const { return continueSubleadingNc; }
 
  void setDensityOperatorEvolution( int scheme, Energy2 cutoff ) {
    densityOperatorEvolution = scheme;
    densityOperatorCutoff = cutoff;
  }
 
  double dipoleKernelForEvolution(size_t em, size_t spec, 
                                  Energy2 pEmitpSpec, Energy2 pEmitpEmis, 
                                  Energy2 pEmispSpec);
 
  void transform(const LorentzRotation& rot);
 
public:
 
  const list<DipoleChain>& chains() const { return theChains; }
 
  list<DipoleChain>& chains() { return theChains; }
 
  const list<DipoleChain>& doneChains() const { return theDoneChains; }
 
  list<DipoleChain>& doneChains() { return theDoneChains; }
 
  bool haveChain() const { return !theChains.empty(); }
 
  DipoleChain& currentChain() { assert(haveChain()); return theChains.front(); }
 
  void popChain();
 
  void popChain(list<DipoleChain>::iterator);
 
  void popChains(const list<list<DipoleChain>::iterator>&);
 
  DipoleIndex 
  mergeIndex(list<Dipole>::iterator firstDipole, const pair<bool,bool>& whichFirst,
             list<Dipole>::iterator secondDipole, const pair<bool,bool>& whichSecond) const;
 
  SubleadingSplittingInfo 
  mergeSplittingInfo(list<DipoleChain>::iterator firstChain, list<Dipole>::iterator firstDipole, 
                     const pair<bool,bool>& whichFirst,
                     list<DipoleChain>::iterator secondChain, list<Dipole>::iterator secondDipole, 
                     const pair<bool,bool>& whichSecond) const;
 
  void getSubleadingSplittings(list<SubleadingSplittingInfo>&);
 
public:
 
  void split(list<Dipole>::iterator dip,
             DipoleSplittingInfo& dsplit,
             pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
             DipoleChain*& firstChain, DipoleChain*& secondChain) {
    split(dip,theChains.begin(),dsplit,childIterators,firstChain,secondChain,false);
  }
 
  void split(list<Dipole>::iterator dip,
             list<DipoleChain>::iterator ch,
             DipoleSplittingInfo& dsplit,
             pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
             DipoleChain*& firstChain, DipoleChain*& secondChain,
             bool colourSpectator = true);
  
  pair<PVector,PVector> tmpsplit(list<Dipole>::iterator dip,
             DipoleSplittingInfo& dsplit,
             pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
             DipoleChain*& firstChain, DipoleChain*& secondChain) {
    return tmpsplit(dip,theChains.begin(),dsplit,childIterators,firstChain,secondChain,false);
  }
  
  pair<PVector,PVector> tmpsplit(list<Dipole>::iterator dip,
             list<DipoleChain>::iterator ch,
             DipoleSplittingInfo& dsplit,
             pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
             DipoleChain*& firstChain, DipoleChain*& secondChain,
             bool colourSpectator = true);
 
 
  void recoil(list<Dipole>::iterator dip,
              list<DipoleChain>::iterator ch,
              DipoleSplittingInfo& dsplit);
 
  void splitSubleading(SubleadingSplittingInfo& dsplit,
                       pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
                       DipoleChain*& firstChain, DipoleChain*& secondChain);
 
  void update(DipoleSplittingInfo& dsplit);
  
  pair<PVector,PVector> tmpupdate(DipoleSplittingInfo& dsplit);
 
  void updateInverse(DipoleSplittingInfo& dsplit);
 
  list<pair<list<Dipole>::iterator,list<DipoleChain>::iterator> >
  inDipoles();
 
  tPPair fillEventRecord(StepPtr step, bool firstInteraction, bool realigned);
 
public:
 
  const map<PPtr,PPtr>& prepare(tSubProPtr subpro,
                                tStdXCombPtr xc,
                                StepPtr step,
                                const pair<PDF,PDF>& pdf,
                                tPPair beam,
                                bool firstInteraction,
                                const set<long>& offShellPartons,
                                bool dipoles = true);
  void slimprepare(tSubProPtr subpro,
                   tStdXCombPtr xc,
                   const pair<PDF,PDF>& pdf,tPPair beam,
                   const set<long>& offShellPartons,
                   bool dipoles = true);
 
  virtual void clear();
 
  void prepareChainsSubleading(const bool decay) {
    static set<long> empty;
    continueSubleadingNc = false;
    PList cordered = colourOrdered(incoming(),outgoing());
    findChains(cordered,empty,decay);
  }
 
public:
 
  void debugLastEvent(ostream&) const;
 
public:
 
  map<PPtr,PerturbativeProcessPtr> & decays() {return theDecays;}
 
  void fillFromDecays(PerturbativeProcessPtr decayProc, vector<PPtr>& original);
 
  void separateDecay(PerturbativeProcessPtr decayProc);
 
  Energy decay(PPtr incoming, bool& powhegEmission);
 
  bool prepareDecay(PerturbativeProcessPtr decayProc,
                                const set<long>& offShellPartons);
 
  void updateDecayMom(PPtr decayParent, PerturbativeProcessPtr decayProc);
 
  void updateDecayChainMom(PPtr decayParent, PerturbativeProcessPtr decayProc);
 
  void updateDecays(PerturbativeProcessPtr decayProc, bool iterate = true);
 
  PerturbativeProcessPtr currentDecay() {return theCurrentDecay;}
 
  void currentDecay(PerturbativeProcessPtr in) {theCurrentDecay=in;}
 
  PPtr nextDecay() {
    if ( !theNextDecays.empty() )
      return theNextDecays.back();
    else
      return PPtr();
  }
  
 
  // SW - Changed from protected to public so that functions can be used in DipoleShowerHandler
public:
 
  void findChains(const PList& ordered, 
                  const set<long>& offShellpartons,
                  const bool decay = false);
 
  PList colourOrdered(PPair & in,PList & out);
 
 
private:
 
  struct getMomentum {
    const Lorentz5Momentum& operator() (PPtr particle) const {
      return particle->momentum();
    }
  };
 
  Lorentz5Momentum thePX;
 
  PList theHard;
 
  map<PPtr,PPtr> theOriginals;
 
  list<DipoleChain> theChains;
 
  list<DipoleChain> theDoneChains;
  
  cPDVector theParticlesAfter;
  
  cPDVector theParticlesBefore;
  
  vector<Lorentz5Momentum> theMomentaAfter;
  
  pair<size_t,pair<size_t,size_t> > theEmitterEmissionIndices;
  
  pair<size_t,size_t> theSpectatorIndices;
  
  DensityOperator theDensityOperator;
 
  bool doSubleadingNc;
  
  bool continueSubleadingNc;
  
  size_t subleadingNcEmissionsLimit;
 
  size_t subEmDone;
 
  
  map<PPtr,size_t> theParticleIndices;
 
  int densityOperatorEvolution;
 
  Energy2 densityOperatorCutoff;
 
  map<std::tuple<size_t,size_t,size_t>,map<size_t,size_t> > theEmissionsMap;
 
private:
 
  map<PPtr,PerturbativeProcessPtr> theDecays;
 
  PerturbativeProcessPtr theCurrentDecay;
 
  PList theNextDecays;
 
};
 
 
}
 
#endif /* HERWIG_DipoleEventRecord_H */