Models/General/PrototypeVertex.h

0001 // -*- C++ -*-
0002 //
0003 // PrototypeVertex.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_PrototypeVertex_H
0010 #define HERWIG_PrototypeVertex_H
0011 #include <stack>
0012 #include "ThePEG/Helicity/Vertex/VertexBase.h"
0013 #include "ThePEG/Persistency/PersistentOStream.h"
0014 #include "ThePEG/Persistency/PersistentIStream.h"
0015 #include "ThePEG/Utilities/EnumIO.h"
0016 #include "NBodyDecayConstructorBase.fh"
0017 //
0018 // This is the declaration of the PrototypeVertex class.
0019 //
0020
0021 namespace Herwig {
0022 using namespace ThePEG;
0023 using Helicity::VertexBasePtr;
0024
0025 class PrototypeVertex;
0026 ThePEG_DECLARE_POINTERS(Herwig::PrototypeVertex,PrototypeVertexPtr);
0027
0028 /** Pair of int,double */
0029 typedef pair<unsigned int, double> CFPair;
0030
0031
0032 /**
0033  *  A struct to order the particles in the same way as in the DecayModes
0034  */
0035 struct ParticleOrdering {
0036   /**
0037    *  Operator for the ordering
0038    * @param p1 The first ParticleData object
0039    * @param p2 The second ParticleData object
0040    */
0041   bool operator() (tcPDPtr p1, tcPDPtr p2) const {
0042     return abs(p1->id()) > abs(p2->id()) ||
0043       ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
0044       ( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
0045   }
0046 };
0047
0048 /**
0049  *  A struct to order the particles in the same way as in the DecayMode's
0050  */
0051 struct VertexOrdering {
0052   /**
0053    *  Operator for the ordering
0054    * @param p1 The first ParticleData object
0055    * @param p2 The second ParticleData object
0056    */
0057   bool operator() (const pair< tPDPtr, PrototypeVertexPtr > & p1,
0058            const pair< tPDPtr, PrototypeVertexPtr > & p2) const  {
0059     return  abs(p1.first->id()) > abs(p2.first->id()) ||
0060       ( abs(p1.first->id()) == abs(p2.first->id()) && p1.first->id() > p2.first->id() ) ||
0061       ( p1.first->id() == p2.first->id() && p1.first->fullName() > p2.first->fullName() );
0062   }
0063 };
0064
0065 typedef multiset<pair< tPDPtr, PrototypeVertexPtr >,VertexOrdering > OrderedVertices;
0066
0067 /**
0068  * A set of ParticleData objects ordered as for the DecayMode's
0069  */
0070 typedef multiset<PDPtr,ParticleOrdering> OrderedParticles;
0071
0072 /**
0073  *  Storage of a potenital n-body decay
0074  */
0075 class PrototypeVertex : public Base {
0076
0077 public:
0078
0079   /**
0080    *  Default Constructor
0081    */
0082   PrototypeVertex() : npart(0), possibleOnShell(false) {}
0083
0084   /**
0085    *  Constructor
0086    */
0087   PrototypeVertex(tPDPtr in, OrderedVertices out,
0088           VertexBasePtr v, int n) :
0089     incoming(in), outgoing(out), vertex(v), npart(n),
0090     possibleOnShell(false) {}
0091
0092   /**
0093    *  Incoming particle
0094    */
0095   tPDPtr incoming;
0096
0097   /**
0098    *  Outgoing particles
0099    */
0100   OrderedVertices outgoing;
0101
0102   /**
0103    *  The vertex for the interaction
0104    */
0105   VertexBasePtr vertex;
0106
0107   /**
0108    *  The parent of the vertex
0109    */
0110   tPrototypeVertexPtr parent;
0111
0112   /**
0113    *  Number of particles
0114    */
0115   unsigned int npart;
0116
0117   /**
0118    *  Outgoing particles
0119    */
0120   mutable OrderedParticles outPart;
0121
0122   /**
0123    *  Can have on-shell intermediates
0124    */
0125   bool possibleOnShell;
0126
0127   /**
0128    *  Increment the number of particles
0129    */
0130   void incrementN(int in) {
0131     npart += in;
0132     if(parent) parent->incrementN(in);
0133   }
0134
0135   /**
0136    *  Mass of the incoming particle
0137    */
0138   Energy incomingMass() {
0139     return incoming->mass();
0140   }
0141
0142   /**
0143    *  Total mass of all the outgoing particles
0144    */
0145   Energy outgoingMass() {
0146     Energy mass(ZERO);
0147     for(OrderedVertices::const_iterator it = outgoing.begin();
0148     it!=outgoing.end();++it) {
0149       mass += it->second ?
0150     it->second->outgoingMass() : it->first->mass();
0151     }
0152     return mass;
0153   }
0154
0155   /**
0156    * Total constituent mass of all the outgoing particles
0157    */
0158   Energy outgoingConstituentMass() {
0159     Energy mass(ZERO);
0160     for(OrderedVertices::const_iterator it = outgoing.begin();
0161     it!=outgoing.end();++it) {
0162       mass += it->second ?
0163     it->second->outgoingConstituentMass() : it->first->constituentMass();
0164     }
0165     return mass;
0166   }
0167
0168   /**
0169    * Check the external particles
0170    */
0171   bool checkExternal(bool first=true) {
0172     if(outPart.empty())     setOutgoing();
0173     if(first&&outPart.find(incoming)!=outPart.end()) return false;
0174     bool output = true;
0175     for(OrderedVertices::const_iterator it = outgoing.begin();
0176     it!=outgoing.end();++it) {
0177       if(it->second&& !it->second->checkExternal(false)) output = false;
0178     }
0179     return output;
0180   }
0181
0182   /**
0183    * Set the outgoing particles
0184    */
0185   void setOutgoing() const {
0186     assert(outPart.empty());
0187     for(OrderedVertices::const_iterator it = outgoing.begin();
0188     it!=outgoing.end();++it) {
0189       if(it->second) {
0190     it->second->setOutgoing();
0191     outPart.insert(it->second->outPart.begin(),
0192                it->second->outPart.end());
0193       }
0194       else
0195     outPart.insert(it->first);
0196     }
0197   }
0198
0199   /**
0200    *  Are there potential on-shell intermediates?
0201    */
0202   bool canBeOnShell(unsigned int opt,Energy maxMass,bool first);
0203
0204   /**
0205    *  Check if same external particles
0206    */
0207   bool sameDecay(const PrototypeVertex & x) const;
0208
0209   /**
0210    *  Create a \f$1\to2\f$ prototype
0211    */
0212   static  void createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
0213                 std::stack<PrototypeVertexPtr> & prototypes,
0214                 NBodyDecayConstructorBasePtr decayCon);
0215
0216   /**
0217    *  Expand the prototypes by adding more legs
0218    */
0219   static void expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
0220                    std::stack<PrototypeVertexPtr> & prototypes,
0221                    const set<PDPtr> & excluded,
0222                    NBodyDecayConstructorBasePtr decayCon);
0223
0224   /**
0225    *  Copy the whole structure with a new branching
0226    */
0227   static PrototypeVertexPtr replicateTree(PrototypeVertexPtr parent,
0228                       PrototypeVertexPtr oldChild,
0229                       PrototypeVertexPtr & newChild);
0230
0231 };
0232
0233 /**
0234  * Output to a stream
0235  */
0236 inline ostream & operator<<(ostream & os, const PrototypeVertex & diag) {
0237   os << diag.incoming->PDGName() << " -> ";
0238   bool seq=false;
0239   for(OrderedVertices::const_iterator it = diag.outgoing.begin();
0240       it!=diag.outgoing.end();++it) {
0241     os << it->first->PDGName() << " ";
0242     if(it->second) seq = true;
0243   }
0244   os << " decays via "
0245      << diag.vertex->fullName() << " in a "
0246      << diag.npart << "-body decay\n";
0247   if(!seq) return os;
0248   os << "Followed by\n";
0249   for(OrderedVertices::const_iterator it = diag.outgoing.begin();
0250       it!=diag.outgoing.end();++it) {
0251     if(it->second) os << *it->second;
0252   }
0253   return os;
0254 }
0255
0256 /**
0257  * Test whether two diagrams are identical.
0258  */
0259 inline bool operator==(const PrototypeVertex & x, const PrototypeVertex & y) {
0260   if(x.incoming != y.incoming) return false;
0261   if(x.vertex != y.vertex) return false;
0262   if(x.npart != y.npart) return false;
0263   if(x.outgoing.empty()&&y.outgoing.empty()) return true;
0264   if(x.outgoing.size() != y.outgoing.size()) return false;
0265   OrderedVertices::const_iterator xt = x.outgoing.begin();
0266   OrderedVertices::const_iterator yt = y.outgoing.begin();
0267   while(xt!=x.outgoing.end()) {
0268     if(xt->first != yt->first) return false;
0269     // special for identical particles
0270     OrderedVertices::const_iterator lxt = x.outgoing.lower_bound(*xt);
0271     OrderedVertices::const_iterator uxt = x.outgoing.upper_bound(*xt);
0272     --uxt;
0273     // just one particle
0274     if(lxt==uxt) {
0275       if(xt->second && yt->second) {
0276     if(*(xt->second)==*(yt->second)) {
0277       ++xt;
0278       ++yt;
0279       continue;
0280     }
0281     else return false;
0282       }
0283       else if(xt->second || yt->second)
0284     return false;
0285       ++xt;
0286       ++yt;
0287     }
0288     // identical particles
0289     else {
0290       ++uxt;
0291       OrderedVertices::const_iterator lyt = y.outgoing.lower_bound(*xt);
0292       OrderedVertices::const_iterator uyt = y.outgoing.upper_bound(*xt);
0293       unsigned int nx=0;
0294       for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {++nx;}
0295       unsigned int ny=0;
0296       for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {++ny;}
0297       if(nx!=ny) return false;
0298       vector<bool> matched(ny,false);
0299       for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {
0300     bool found = false;
0301     unsigned int iy=0;
0302     for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {
0303       if(matched[iy]) {
0304         ++iy;
0305         continue;
0306       }
0307       if( (!ixt->second &&!iyt->second)  ||
0308           ( ixt->second&&iyt->second &&
0309         *(ixt->second)==*(iyt->second)) ) {
0310         matched[iy] = true;
0311         found = true;
0312         break;
0313       }
0314       ++iy;
0315     }
0316     if(!found) return false;
0317       }
0318       xt=uxt;
0319       yt=uyt;
0320     }
0321   }
0322   return true;
0323 }
0324
0325 /**
0326  *  A simple vertex for the N-body diagram
0327  */
0328 struct NBVertex {
0329
0330   /**
0331    * Constructor taking a prototype vertex as the arguments
0332    */
0333   NBVertex(PrototypeVertexPtr proto = PrototypeVertexPtr() );
0334
0335   /**
0336    * Incoming particle
0337    */
0338   tPDPtr incoming;
0339
0340   /**
0341    *  Outgoing particles
0342    */
0343   mutable OrderedParticles outgoing;
0344
0345   /**
0346    *  The vertices
0347    */
0348   list<pair<PDPtr,NBVertex> > vertices;
0349
0350   /**
0351    *  The vertex
0352    */
0353   VertexBasePtr vertex;
0354 };
0355
0356 /**
0357  * The NBDiagram struct contains information about a \f$1\to n\f$ decay
0358  * that has been automatically generated.
0359  */
0360   struct NBDiagram : public NBVertex {
0361
0362   /**
0363    * Constructor taking a prototype vertex as the arguments*/
0364   NBDiagram(PrototypeVertexPtr proto=PrototypeVertexPtr());
0365
0366   /**
0367    *  The type of channel
0368    */
0369   vector<unsigned int> channelType;
0370
0371   /** Store colour flow at \f$N_c=3\f$ information */
0372   mutable vector<CFPair> colourFlow;
0373
0374   /** Store colour flow at \f$N_c=\infty\f$ information */
0375   mutable vector<CFPair> largeNcColourFlow;
0376 };
0377
0378 /**
0379  * Output operator to allow the structure to be persistently written
0380  * @param os The output stream
0381  * @param x The NBVertex
0382  */
0383 inline PersistentOStream & operator<<(PersistentOStream & os,
0384                       const NBVertex  & x) {
0385   os << x.incoming << x.outgoing << x.vertices << x.vertex;
0386   return os;
0387 }
0388
0389 /**
0390  * Input operator to allow persistently written data to be read in
0391  * @param is The input stream
0392  * @param x The NBVertex
0393  */
0394 inline PersistentIStream & operator>>(PersistentIStream & is,
0395                       NBVertex & x) {
0396   is >> x.incoming >> x.outgoing >> x.vertices >> x.vertex;
0397   return is;
0398 }
0399
0400 /**
0401  * Output operator to allow the structure to be persistently written
0402  * @param os The output stream
0403  * @param x The NBDiagram
0404  */
0405 inline PersistentOStream & operator<<(PersistentOStream & os,
0406                       const NBDiagram  & x) {
0407   os << x.incoming << x.channelType << x.outgoing << x.vertices << x.vertex
0408      << x.colourFlow << x.largeNcColourFlow;
0409   return os;
0410 }
0411
0412 /**
0413  * Input operator to allow persistently written data to be read in
0414  * @param is The input stream
0415  * @param x The NBDiagram
0416  */
0417 inline PersistentIStream & operator>>(PersistentIStream & is,
0418                       NBDiagram & x) {
0419   is >> x.incoming >> x.channelType >> x.outgoing >> x.vertices >> x.vertex
0420      >> x.colourFlow >> x.largeNcColourFlow;
0421   return is;
0422 }
0423
0424 /**
0425  * Output a NBVertex to a stream
0426  */
0427 inline ostream & operator<<(ostream & os, const NBVertex & vertex) {
0428   os << vertex.incoming->PDGName() << " -> ";
0429   bool seq=false;
0430   for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
0431       it!=vertex.vertices.end();++it) {
0432     os << it->first->PDGName() << " ";
0433     if(it->second.incoming) seq = true;
0434   }
0435   os << "via vertex " << vertex.vertex->fullName() << "\n";
0436   if(!seq) return os;
0437   os << "Followed by\n";
0438   for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
0439       it!=vertex.vertices.end();++it) {
0440     if(it->second.incoming) os << it->second;
0441   }
0442   return os;
0443 }
0444
0445 /**
0446  * Output a NBDiagram to a stream
0447  */
0448 inline ostream & operator<<(ostream & os, const NBDiagram & diag) {
0449   os << diag.incoming->PDGName() << " -> ";
0450   for(OrderedParticles::const_iterator it=diag.outgoing.begin();
0451       it!=diag.outgoing.end();++it) {
0452     os << (**it).PDGName() << " ";
0453   }
0454   os << " has order ";
0455   for(unsigned int ix=0;ix<diag.channelType.size();++ix)
0456     os << diag.channelType[ix] << " ";
0457   os << "\n";
0458   os << "First decay " << diag.incoming->PDGName() << " -> ";
0459   bool seq=false;
0460   for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
0461       it!=diag.vertices.end();++it) {
0462     os << it->first->PDGName() << " ";
0463     if(it->second.incoming) seq = true;
0464   }
0465   os << "via vertex " << diag.vertex->fullName() << "\n";
0466   if(!seq) return os;
0467   os << "Followed by\n";
0468   for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
0469       it!=diag.vertices.end();++it) {
0470     if(it->second.incoming) os << it->second;
0471   }
0472   return os;
0473 }
0474
0475 }
0476
0477 #endif /* HERWIG_PrototypeVertex_H */