GeneralHardME.h

// -*- C++ -*-
//
// GeneralHardME.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_GeneralHardME_H
#define HERWIG_GeneralHardME_H
//
// This is the declaration of the GeneralHardME class.
//
 
#include "Herwig/MatrixElement/HwMEBase.h"
#include "ThePEG/Utilities/Exception.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Models/General/HPDiagram.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
#include "Herwig/MatrixElement/HardVertex.h"
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
#include "GeneralHardME.fh"
 
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
 
class GeneralHardME: public HwMEBase {
 
public:
 
  typedef vector<HPDiagram>::size_type HPCount;
 
  enum ColourStructure {UNDEFINED,
                        Colour11to11,Colour11to33bar,Colour11to88,
                        Colour33to33,Colour33barto11,Colour33barto33bar,
                        Colour33barto66bar, Colour33barto6bar6,
                        Colour33to61, Colour3bar3barto6bar1,
                        Colour33to16, Colour3bar3barto16bar,
                        Colour38to3bar6, Colour38to63bar,
                        Colour33barto18,Colour33barto81,Colour33barto88,
                        Colour38to13,Colour38to31,
                        Colour38to83,Colour38to38,
                        Colour3bar3barto3bar3bar,
                        Colour3bar8to13bar,Colour3bar8to3bar1,
                        Colour3bar8to83bar,Colour3bar8to3bar8,
                        Colour88to11,Colour88to33bar,
                        Colour88to66bar,Colour88to88,
                        Colour88to18,Colour88to81,
                        Colour33to13bar,Colour33to3bar1,
                        Colour33to83bar,Colour33to3bar8,
                        Colour3bar3barto13,Colour3bar3barto31,
                        Colour3bar3barto83,Colour3bar3barto38,
                        Colour38to3bar3bar,Colour3bar8to33};
 
public:
 
  GeneralHardME();
 
public:
 
  virtual unsigned int orderInAlphaS() const;
 
  virtual unsigned int orderInAlphaEW() const;
 
  virtual double me2() const = 0;
 
  virtual Energy2 scale() const {
    if(scaleChoice_==0) {
      return scaleFactor_*sHat();
    }
    else if(scaleChoice_==1) {
      Energy2 mbar = 0.5*(meMomenta()[2].mass2()+meMomenta()[3].mass2());
      Energy2 t = 0.5*(tHat()-mbar);
      Energy2 u = 0.5*(uHat()-mbar);
      Energy2 s = 0.5*sHat();
      return scaleFactor_*4.*s*t*u/(s*s+t*t+u*u);
    }
    else if(scaleChoice_ ==2) {
      Energy2 scale1 = meMomenta()[2].mass2()+meMomenta()[2].perp2();
      Energy2 scale2 = meMomenta()[3].mass2()+meMomenta()[3].perp2();
      return scaleFactor_*max(scale1,scale2);
    }
    else {
      assert(false);
      return ZERO;
    }
  }
 
  virtual void getDiagrams() const;
 
  virtual Selector<DiagramIndex> 
  diagrams(const DiagramVector & dv) const;
 
  virtual Selector<const ColourLines *>
  colourGeometries(tcDiagPtr diag) const;
 
  void setProcessInfo(const vector<HPDiagram> & process,
                      ColourStructure colour, bool debug, 
                      unsigned int scaleOption,
                      double scaleFactor);
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual void doinitrun();
 
protected:
 
  virtual void debug(double ) const {}  
 
protected:
 
  const vector<HPDiagram> & getProcessInfo() const {
    return diagrams_;
  }
 
  pair<long, long> getIncoming() const {
    return incoming_;
  }
  
  pair<long, long> getOutgoing() const {
    return outgoing_;
  }
  
  const vector<DVector> & getColourFactors() const {
    return colour_;
  }
 
  HPCount numberOfDiags() const {
    return numberOfDiagrams_;
  } 
  
  size_t numberOfFlows() const {
    return numberOfFlows_;
  }
 
  bool debugME() const {
    return debug_;
  }
 
  unsigned int colourFlow() const {return flow_;}
 
  void colourFlow(unsigned int flow) const {flow_=flow;}
 
  unsigned int diagram() const {return diagram_;}
 
  void diagram(unsigned int diag) const {diagram_=diag;}
 
  double selectColourFlow(vector<double> & flow,
                          vector<double> & me,double average) const;
 
  vector<ProductionMatrixElement> & flowME() const {
    return flowME_;
  }
 
  vector<ProductionMatrixElement> & diagramME() const {
    return diagramME_;
  }
 
  ColourStructure colour() const {return colourStructure_;}
 
  ParticleVector hardParticles(tSubProPtr subp) {
    ParticleVector output(4);
    output[0] = subp->incoming().first; 
    output[1] = subp->incoming().second;
    output[2] = subp->outgoing()[0]; 
    output[3] = subp->outgoing()[1];    
    //ensure particle ordering is the same as it was when
    //the diagrams were created
    if( output[0]->id() != getIncoming().first )
      swap(output[0], output[1]);
    if( output[2]->id() != getOutgoing().first )
      swap(output[2], output[3]);
    // return answer
    return output;
  }
 
  void setRescaledMomenta(const ParticleVector & external) {
    cPDVector data(4);
    vector<Lorentz5Momentum> momenta(4);
    for( size_t i = 0; i < 4; ++i ) {
      data[i] = external[i]->dataPtr();
      momenta[i] = external[i]->momentum();
    }
    rescaleMomenta(momenta, data);
  }
 
  void createVertex(ProductionMatrixElement & me,
                    ParticleVector & external) {
    HardVertexPtr hardvertex = new_ptr(HardVertex());
    hardvertex->ME(me);
    for(ParticleVector::size_type i = 0; i < 4; ++i) {
      tSpinPtr spin = external[i]->spinInfo();
      if(i<2) {
        tcPolarizedBeamPDPtr beam = 
          dynamic_ptr_cast<tcPolarizedBeamPDPtr>(external[i]->dataPtr());
        if(beam) spin->rhoMatrix() = beam->rhoMatrix();
      }
      spin->productionVertex(hardvertex);
    }
  }
 
  void initializeMatrixElements(PDT::Spin  in1, PDT::Spin in2,
                                PDT::Spin out1, PDT::Spin out2) {
    flowME().resize(numberOfFlows(),
                    ProductionMatrixElement(in1,in2,out1,out2));
    diagramME().resize(numberOfDiags(),
                       ProductionMatrixElement(in1,in2,out1,out2));
  }
 
private:
 
  GeneralHardME & operator=(const GeneralHardME &) = delete;
 
private:
  
  pair<long, long> incoming_;
  
  pair<long, long> outgoing_;
 
  vector<HPDiagram> diagrams_;
 
  HPCount numberOfDiagrams_;
 
  ColourStructure colourStructure_;
 
  vector<DVector> colour_;
 
  unsigned int numberOfFlows_;
 
  bool debug_;
 
  unsigned int scaleChoice_;
 
  double scaleFactor_;
 
  mutable unsigned int flow_;
 
  mutable unsigned int diagram_;
 
  mutable vector<ProductionMatrixElement> flowME_;
 
  mutable vector<ProductionMatrixElement> diagramME_;
 
};
 
class MEException : public Exception {};
  
}
 
#endif /* HERWIG_GeneralHardME_H */