PrototypeVertex.h

// -*- C++ -*-
//
// PrototypeVertex.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_PrototypeVertex_H
#define HERWIG_PrototypeVertex_H
#include <stack>
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "NBodyDecayConstructorBase.fh"
//
// This is the declaration of the PrototypeVertex class.
//
 
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
 
class PrototypeVertex;
ThePEG_DECLARE_POINTERS(Herwig::PrototypeVertex,PrototypeVertexPtr);
  
typedef pair<unsigned int, double> CFPair;
 
 
struct ParticleOrdering {
  bool operator() (tcPDPtr p1, tcPDPtr p2) const {
    return abs(p1->id()) > abs(p2->id()) ||
      ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
      ( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
  }
};
 
struct VertexOrdering {
  bool operator() (const pair< tPDPtr, PrototypeVertexPtr > & p1,
                   const pair< tPDPtr, PrototypeVertexPtr > & p2) const  {
    return  abs(p1.first->id()) > abs(p2.first->id()) ||
      ( abs(p1.first->id()) == abs(p2.first->id()) && p1.first->id() > p2.first->id() ) ||
      ( p1.first->id() == p2.first->id() && p1.first->fullName() > p2.first->fullName() );
  }
};
  
typedef multiset<pair< tPDPtr, PrototypeVertexPtr >,VertexOrdering > OrderedVertices;
 
typedef multiset<PDPtr,ParticleOrdering> OrderedParticles;
 
class PrototypeVertex : public Base {
  
public:
 
  PrototypeVertex() : npart(0), possibleOnShell(false) {}
 
  PrototypeVertex(tPDPtr in, OrderedVertices out,
                  VertexBasePtr v, int n) :
    incoming(in), outgoing(out), vertex(v), npart(n),
    possibleOnShell(false) {}
 
  tPDPtr incoming;
 
  OrderedVertices outgoing;
 
  VertexBasePtr vertex;
 
  tPrototypeVertexPtr parent;
  
  unsigned int npart;
 
  mutable OrderedParticles outPart;
 
  bool possibleOnShell;
 
  void incrementN(int in) {
    npart += in;
    if(parent) parent->incrementN(in);
  }
 
  Energy incomingMass() {
    return incoming->mass();
  }
 
  Energy outgoingMass() {
    Energy mass(ZERO);
    for(OrderedVertices::const_iterator it = outgoing.begin();
        it!=outgoing.end();++it) {
      mass += it->second ? 
        it->second->outgoingMass() : it->first->mass();
    }
    return mass;
  }
 
  Energy outgoingConstituentMass() {
    Energy mass(ZERO);
    for(OrderedVertices::const_iterator it = outgoing.begin();
        it!=outgoing.end();++it) {
      mass += it->second ? 
        it->second->outgoingConstituentMass() : it->first->constituentMass();
    }
    return mass;
  }
 
  bool checkExternal(bool first=true) {
    if(outPart.empty())     setOutgoing();
    if(first&&outPart.find(incoming)!=outPart.end()) return false;
    bool output = true;
    for(OrderedVertices::const_iterator it = outgoing.begin();
        it!=outgoing.end();++it) {
      if(it->second&& !it->second->checkExternal(false)) output = false;
    }
    return output;
  }
  
  void setOutgoing() const {
    assert(outPart.empty());
    for(OrderedVertices::const_iterator it = outgoing.begin();
        it!=outgoing.end();++it) {
      if(it->second) {
        it->second->setOutgoing();
        outPart.insert(it->second->outPart.begin(),
                       it->second->outPart.end());
      }
      else
        outPart.insert(it->first);
    }
  }
 
  bool canBeOnShell(unsigned int opt,Energy maxMass,bool first);
 
  bool sameDecay(const PrototypeVertex & x) const;
 
  static  void createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
                                std::stack<PrototypeVertexPtr> & prototypes,
                                NBodyDecayConstructorBasePtr decayCon);
 
  static void expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
                               std::stack<PrototypeVertexPtr> & prototypes,
                               const set<PDPtr> & excluded,
                               NBodyDecayConstructorBasePtr decayCon);
 
  static PrototypeVertexPtr replicateTree(PrototypeVertexPtr parent,
                                          PrototypeVertexPtr oldChild,
                                          PrototypeVertexPtr & newChild);
 
};
 
inline ostream & operator<<(ostream & os, const PrototypeVertex & diag) {
  os << diag.incoming->PDGName() << " -> ";
  bool seq=false;
  for(OrderedVertices::const_iterator it = diag.outgoing.begin();
      it!=diag.outgoing.end();++it) {
    os << it->first->PDGName() << " ";
    if(it->second) seq = true;
  }
  os << " decays via "
     << diag.vertex->fullName() << " in a " 
     << diag.npart << "-body decay\n";
  if(!seq) return os;
  os << "Followed by\n";
  for(OrderedVertices::const_iterator it = diag.outgoing.begin();
      it!=diag.outgoing.end();++it) {
    if(it->second) os << *it->second;
  }
  return os;
}
 
inline bool operator==(const PrototypeVertex & x, const PrototypeVertex & y) {
  if(x.incoming != y.incoming) return false;
  if(x.vertex != y.vertex) return false;
  if(x.npart != y.npart) return false;
  if(x.outgoing.empty()&&y.outgoing.empty()) return true;
  if(x.outgoing.size() != y.outgoing.size()) return false;
  OrderedVertices::const_iterator xt = x.outgoing.begin();
  OrderedVertices::const_iterator yt = y.outgoing.begin();
  while(xt!=x.outgoing.end()) {
    if(xt->first != yt->first) return false;
    // special for identical particles
    OrderedVertices::const_iterator lxt = x.outgoing.lower_bound(*xt);
    OrderedVertices::const_iterator uxt = x.outgoing.upper_bound(*xt);
    --uxt;
    // just one particle
    if(lxt==uxt) {
      if(xt->second && yt->second) {
        if(*(xt->second)==*(yt->second)) {
          ++xt;
          ++yt;
          continue;
        }
        else return false;
      }
      else if(xt->second || yt->second)
        return false;
      ++xt;
      ++yt;
    }
    // identical particles
    else {
      ++uxt;
      OrderedVertices::const_iterator lyt = y.outgoing.lower_bound(*xt);
      OrderedVertices::const_iterator uyt = y.outgoing.upper_bound(*xt);
      unsigned int nx=0;
      for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {++nx;}
      unsigned int ny=0;
      for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {++ny;}
      if(nx!=ny) return false;
      vector<bool> matched(ny,false);
      for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {
        bool found = false;
        unsigned int iy=0;
        for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {
          if(matched[iy]) {
            ++iy;
            continue;
          }
          if( (!ixt->second &&!iyt->second)  ||
              ( ixt->second&&iyt->second &&
                *(ixt->second)==*(iyt->second)) ) {
            matched[iy] = true;
            found = true;
            break;
          }
          ++iy;
        }
        if(!found) return false;
      }
      xt=uxt;
      yt=uyt;
    }
  }
  return true;
}
 
struct NBVertex {
 
  NBVertex(PrototypeVertexPtr proto = PrototypeVertexPtr() );
 
  tPDPtr incoming;
 
  mutable OrderedParticles outgoing;
 
  list<pair<PDPtr,NBVertex> > vertices;
 
  VertexBasePtr vertex;
};
 
  struct NBDiagram : public NBVertex {
 
  NBDiagram(PrototypeVertexPtr proto=PrototypeVertexPtr());
 
  vector<unsigned int> channelType;
 
  mutable vector<CFPair> colourFlow;
  
  mutable vector<CFPair> largeNcColourFlow;
};
 
inline PersistentOStream & operator<<(PersistentOStream & os, 
                                      const NBVertex  & x) {
  os << x.incoming << x.outgoing << x.vertices << x.vertex;
  return os;
}
  
inline PersistentIStream & operator>>(PersistentIStream & is,
                                      NBVertex & x) {
  is >> x.incoming >> x.outgoing >> x.vertices >> x.vertex;
  return is;
}
 
inline PersistentOStream & operator<<(PersistentOStream & os, 
                                      const NBDiagram  & x) {
  os << x.incoming << x.channelType << x.outgoing << x.vertices << x.vertex
     << x.colourFlow << x.largeNcColourFlow;
  return os;
}
  
inline PersistentIStream & operator>>(PersistentIStream & is,
                                      NBDiagram & x) {
  is >> x.incoming >> x.channelType >> x.outgoing >> x.vertices >> x.vertex
     >> x.colourFlow >> x.largeNcColourFlow;
  return is;
}
 
inline ostream & operator<<(ostream & os, const NBVertex & vertex) {
  os << vertex.incoming->PDGName() << " -> ";
  bool seq=false;
  for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
      it!=vertex.vertices.end();++it) {
    os << it->first->PDGName() << " ";
    if(it->second.incoming) seq = true;
  }
  os << "via vertex " << vertex.vertex->fullName() << "\n";
  if(!seq) return os;
  os << "Followed by\n";
  for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
      it!=vertex.vertices.end();++it) {
    if(it->second.incoming) os << it->second;
  }
  return os;
}
 
inline ostream & operator<<(ostream & os, const NBDiagram & diag) {
  os << diag.incoming->PDGName() << " -> ";
  for(OrderedParticles::const_iterator it=diag.outgoing.begin();
      it!=diag.outgoing.end();++it) {
    os << (**it).PDGName() << " ";
  }
  os << " has order ";
  for(unsigned int ix=0;ix<diag.channelType.size();++ix)
    os << diag.channelType[ix] << " ";
  os << "\n";
  os << "First decay " << diag.incoming->PDGName() << " -> ";
  bool seq=false;
  for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
      it!=diag.vertices.end();++it) {
    os << it->first->PDGName() << " ";
    if(it->second.incoming) seq = true;
  }
  os << "via vertex " << diag.vertex->fullName() << "\n";
  if(!seq) return os;
  os << "Followed by\n";
  for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
      it!=diag.vertices.end();++it) {
    if(it->second.incoming) os << it->second;
  }
  return os;
}
 
}
 
#endif /* HERWIG_PrototypeVertex_H */