MatrixElement/General/GeneralHardME.h

0001 // -*- C++ -*-
0002 //
0003 // GeneralHardME.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_GeneralHardME_H
0010 #define HERWIG_GeneralHardME_H
0011 //
0012 // This is the declaration of the GeneralHardME class.
0013 //
0014
0015 #include "Herwig/MatrixElement/HwMEBase.h"
0016 #include "ThePEG/Utilities/Exception.h"
0017 #include "ThePEG/Persistency/PersistentOStream.h"
0018 #include "ThePEG/Persistency/PersistentIStream.h"
0019 #include "Herwig/Models/General/HPDiagram.h"
0020 #include "Herwig/MatrixElement/ProductionMatrixElement.h"
0021 #include "Herwig/MatrixElement/HardVertex.h"
0022 #include "ThePEG/EventRecord/SpinInfo.h"
0023 #include "ThePEG/PDF/PolarizedBeamParticleData.h"
0024 #include "GeneralHardME.fh"
0025
0026 namespace Herwig {
0027 using namespace ThePEG;
0028 using Helicity::VertexBasePtr;
0029
0030 /**
0031  * This defines the GeneralHardME class that is designed to serve as a
0032  * base class for matrix elements of specific spin structures when those
0033  * structures are created by a a general model, i.e. a SUSY production
0034  * ME. It stores a vector of diagram structures that contain the required
0035  * to calculate the matrix element.
0036  *
0037  * @see HwMEBase
0038  */
0039
0040 class GeneralHardME: public HwMEBase {
0041
0042 public:
0043
0044   /**
0045    * Convenient typedef for size_type of HPDiagram vector
0046    */
0047   typedef vector<HPDiagram>::size_type HPCount;
0048
0049   /**
0050    *  Enum for the possible colour structures
0051    */
0052   enum ColourStructure {UNDEFINED,
0053                         Colour11to11,Colour11to33bar,Colour11to88,
0054                         Colour33to33,Colour33barto11,Colour33barto33bar,
0055                         Colour33barto66bar, Colour33barto6bar6,
0056             Colour33to61, Colour3bar3barto6bar1,
0057                         Colour33to16, Colour3bar3barto16bar,
0058             Colour38to3bar6, Colour38to63bar,
0059                         Colour33barto18,Colour33barto81,Colour33barto88,
0060                         Colour38to13,Colour38to31,
0061                         Colour38to83,Colour38to38,
0062                         Colour3bar3barto3bar3bar,
0063                         Colour3bar8to13bar,Colour3bar8to3bar1,
0064                         Colour3bar8to83bar,Colour3bar8to3bar8,
0065                         Colour88to11,Colour88to33bar,
0066                         Colour88to66bar,Colour88to88,
0067                         Colour88to18,Colour88to81,
0068             Colour33to13bar,Colour33to3bar1,
0069             Colour33to83bar,Colour33to3bar8,
0070             Colour3bar3barto13,Colour3bar3barto31,
0071             Colour3bar3barto83,Colour3bar3barto38,
0072             Colour38to3bar3bar,Colour3bar8to33};
0073
0074 public:
0075
0076   /**
0077    * The default constructor.
0078    */
0079   GeneralHardME();
0080
0081 public:
0082
0083   /** @name Virtual functions required by the MEBase class. */
0084   //@{
0085   /**
0086    * Return the order in \f$\alpha_S\f$ in which this matrix
0087    * element is given.
0088    */
0089   virtual unsigned int orderInAlphaS() const;
0090
0091   /**
0092    * Return the order in \f$\alpha_{EW}\f$ in which this matrix
0093    * element is given.
0094    */
0095   virtual unsigned int orderInAlphaEW() const;
0096
0097   /**
0098    * The matrix element for the kinematical configuration
0099    * previously provided by the last call to setKinematics(), suitably
0100    * scaled by sHat() to give a dimension-less number.
0101    * @return the matrix element scaled with sHat() to give a
0102    * dimensionless number.
0103    */
0104   virtual double me2() const = 0;
0105
0106   /**
0107    * Return the scale associated with the last set phase space point.
0108    */
0109   virtual Energy2 scale() const {
0110     if(scaleChoice_==0) {
0111       return scaleFactor_*sHat();
0112     }
0113     else if(scaleChoice_==1) {
0114       Energy2 mbar = 0.5*(meMomenta()[2].mass2()+meMomenta()[3].mass2());
0115       Energy2 t = 0.5*(tHat()-mbar);
0116       Energy2 u = 0.5*(uHat()-mbar);
0117       Energy2 s = 0.5*sHat();
0118       return scaleFactor_*4.*s*t*u/(s*s+t*t+u*u);
0119     }
0120     else if(scaleChoice_ ==2) {
0121       Energy2 scale1 = meMomenta()[2].mass2()+meMomenta()[2].perp2();
0122       Energy2 scale2 = meMomenta()[3].mass2()+meMomenta()[3].perp2();
0123       return scaleFactor_*max(scale1,scale2);
0124     }
0125     else {
0126       assert(false);
0127       return ZERO;
0128     }
0129   }
0130
0131   /**
0132    * Add all possible diagrams with the add() function.
0133    */
0134   virtual void getDiagrams() const;
0135
0136   /**
0137    * Get diagram selector. With the information previously supplied with the
0138    * setKinematics method, a derived class may optionally
0139    * override this method to weight the given diagrams with their
0140    * (although certainly not physical) relative probabilities.
0141    * @param dv the diagrams to be weighted.
0142    * @return a Selector relating the given diagrams to their weights.
0143    */
0144   virtual Selector<DiagramIndex>
0145   diagrams(const DiagramVector & dv) const;
0146
0147   /**
0148    * Return a Selector with possible colour geometries for the selected
0149    * diagram weighted by their relative probabilities.
0150    * @param diag the diagram chosen.
0151    * @return the possible colour geometries weighted by their
0152    * relative probabilities.
0153    */
0154   virtual Selector<const ColourLines *>
0155   colourGeometries(tcDiagPtr diag) const;
0156   //@}
0157
0158   /**
0159    * Set the diagrams and matrix of colour factors.
0160    * @param process vector of MEDiagram with information that
0161    * will allow the diagrams to be created in the specific matrix element
0162    * @param colour The colour structure for the process
0163    * @param debug Whether to compare the numerical answer to an analytical
0164    * formula (This is only stored for certain processes. It is intended
0165    * for quick checks of the matrix elements).
0166    * @param scaleOption The option of what scale to use
0167    * @param scaleFactor The prefactor for the scale
0168    */
0169   void setProcessInfo(const vector<HPDiagram> & process,
0170               ColourStructure colour, bool debug,
0171               unsigned int scaleOption,
0172               double scaleFactor);
0173
0174 public:
0175
0176   /** @name Functions used by the persistent I/O system. */
0177   //@{
0178   /**
0179    * Function used to write out object persistently.
0180    * @param os the persistent output stream written to.
0181    */
0182   void persistentOutput(PersistentOStream & os) const;
0183
0184   /**
0185    * Function used to read in object persistently.
0186    * @param is the persistent input stream read from.
0187    * @param version the version number of the object when written.
0188    */
0189   void persistentInput(PersistentIStream & is, int version);
0190   //@}
0191
0192   /**
0193    * The standard Init function used to initialize the interfaces.
0194    * Called exactly once for each class by the class description system
0195    * before the main function starts or
0196    * when this class is dynamically loaded.
0197    */
0198   static void Init();
0199
0200 protected:
0201
0202   /** @name Standard Interfaced functions. */
0203   //@{
0204   /**
0205    * Initialize this object. Called in the run phase just before
0206    * a run begins.
0207    */
0208   virtual void doinitrun();
0209   //@}
0210
0211 protected:
0212
0213   /**
0214    * A debugging function to test the value of me2 against an
0215    * analytic function. This is to be overidden in an inheriting class.
0216    */
0217   virtual void debug(double ) const {}
0218
0219 protected:
0220
0221   /**
0222    * Access the HPDiagrams that store the required information
0223    * to create the diagrams
0224    */
0225   const vector<HPDiagram> & getProcessInfo() const {
0226     return diagrams_;
0227   }
0228
0229   /**
0230    * Return the incoming pair
0231    * @return Pair of particle ids for the incoming particles
0232    */
0233   pair<long, long> getIncoming() const {
0234     return incoming_;
0235   }
0236
0237   /**
0238    * Return the outgoing pair
0239    * @return Pair of particle ids for the outgoing particles
0240    */
0241   pair<long, long> getOutgoing() const {
0242     return outgoing_;
0243   }
0244
0245   /**
0246    * Return the matrix of colour factors
0247    */
0248   const vector<DVector> & getColourFactors() const {
0249     return colour_;
0250   }
0251
0252   /**
0253    * Get the number of diagrams in this process
0254    */
0255   HPCount numberOfDiags() const {
0256     return numberOfDiagrams_;
0257   }
0258
0259   /**
0260    * Access number of colour flows
0261    */
0262   size_t numberOfFlows() const {
0263     return numberOfFlows_;
0264   }
0265
0266   /**
0267    * Whether to print the debug information
0268    */
0269   bool debugME() const {
0270     return debug_;
0271   }
0272
0273   /**
0274    *  Set/Get Info on the selected diagram and colour flow
0275    */
0276   //@{
0277   /**
0278    * Colour flow
0279    */
0280   unsigned int colourFlow() const {return flow_;}
0281
0282   /**
0283    * Colour flow
0284    */
0285   void colourFlow(unsigned int flow) const {flow_=flow;}
0286
0287   /**
0288    * Diagram
0289    */
0290   unsigned int diagram() const {return diagram_;}
0291
0292   /**
0293    * Diagram
0294    */
0295   void diagram(unsigned int diag) const {diagram_=diag;}
0296   //@}
0297
0298   /**
0299    *  Calculate weight and select colour flow
0300    */
0301   double selectColourFlow(vector<double> & flow,
0302               vector<double> & me,double average) const;
0303
0304   /**
0305    *  Access to the colour flow matrix element
0306    */
0307   vector<ProductionMatrixElement> & flowME() const {
0308     return flowME_;
0309   }
0310
0311   /**
0312    *  Access to the diagram matrix element
0313    */
0314   vector<ProductionMatrixElement> & diagramME() const {
0315     return diagramME_;
0316   }
0317
0318   /**
0319    *  Access to the colour structure
0320    */
0321   ColourStructure colour() const {return colourStructure_;}
0322
0323   /**
0324    *  Extract the paricles from the subprocess
0325    */
0326   ParticleVector hardParticles(tSubProPtr subp) {
0327     ParticleVector output(4);
0328     output[0] = subp->incoming().first;
0329     output[1] = subp->incoming().second;
0330     output[2] = subp->outgoing()[0];
0331     output[3] = subp->outgoing()[1];
0332     //ensure particle ordering is the same as it was when
0333     //the diagrams were created
0334     if( output[0]->id() != getIncoming().first )
0335       swap(output[0], output[1]);
0336     if( output[2]->id() != getOutgoing().first )
0337       swap(output[2], output[3]);
0338     // return answer
0339     return output;
0340   }
0341
0342   /**
0343    *  Set the rescaled momenta
0344    */
0345   void setRescaledMomenta(const ParticleVector & external) {
0346     cPDVector data(4);
0347     vector<Lorentz5Momentum> momenta(4);
0348     for( size_t i = 0; i < 4; ++i ) {
0349       data[i] = external[i]->dataPtr();
0350       momenta[i] = external[i]->momentum();
0351     }
0352     rescaleMomenta(momenta, data);
0353   }
0354
0355   /**
0356    *  Create the vertes
0357    */
0358   void createVertex(ProductionMatrixElement & me,
0359             ParticleVector & external) {
0360     HardVertexPtr hardvertex = new_ptr(HardVertex());
0361     hardvertex->ME(me);
0362     for(ParticleVector::size_type i = 0; i < 4; ++i) {
0363       tSpinPtr spin = external[i]->spinInfo();
0364       if(i<2) {
0365     tcPolarizedBeamPDPtr beam =
0366       dynamic_ptr_cast<tcPolarizedBeamPDPtr>(external[i]->dataPtr());
0367     if(beam) spin->rhoMatrix() = beam->rhoMatrix();
0368       }
0369       spin->productionVertex(hardvertex);
0370     }
0371   }
0372
0373   /**
0374    *  Initialize the storage of the helicity matrix elements
0375    */
0376   void initializeMatrixElements(PDT::Spin  in1, PDT::Spin in2,
0377                 PDT::Spin out1, PDT::Spin out2) {
0378     flowME().resize(numberOfFlows(),
0379             ProductionMatrixElement(in1,in2,out1,out2));
0380     diagramME().resize(numberOfDiags(),
0381                ProductionMatrixElement(in1,in2,out1,out2));
0382   }
0383
0384 private:
0385
0386   /**
0387    * The assignment operator is private and must never be called.
0388    * In fact, it should not even be implemented.
0389    */
0390   GeneralHardME & operator=(const GeneralHardME &) = delete;
0391
0392 private:
0393
0394   /**
0395    *  External particles
0396    */
0397   //@{
0398   /**
0399    * Store incoming particles
0400    */
0401   pair<long, long> incoming_;
0402
0403   /**
0404    * Store the outgoing particles
0405    */
0406   pair<long, long> outgoing_;
0407   //@}
0408
0409   /**
0410    *  Diagrams
0411    */
0412   //@{
0413   /**
0414    * Store all diagrams as a vector of structures
0415    */
0416   vector<HPDiagram> diagrams_;
0417
0418   /**
0419    * Store the number of diagrams for fast retrieval
0420    */
0421   HPCount numberOfDiagrams_;
0422   //@}
0423
0424   /**
0425    *  Colour information
0426    */
0427   //@{
0428   /**
0429    *  The colour structure
0430    */
0431   ColourStructure colourStructure_;
0432
0433   /**
0434    * Store colour factors for ME calc.
0435    */
0436   vector<DVector> colour_;
0437
0438   /**
0439    * The number of colourflows.
0440    */
0441   unsigned int numberOfFlows_;
0442   //@}
0443
0444   /**
0445    * Whether to test the value of me2 against the analytical function
0446    */
0447   bool debug_;
0448
0449   /**
0450    * The scale chocie
0451    */
0452   unsigned int scaleChoice_;
0453
0454   /**
0455    *  The scale factor
0456    */
0457   double scaleFactor_;
0458
0459   /**
0460    *  Info on the selected diagram and colour flow
0461    */
0462   //@{
0463   /**
0464    * Colour flow
0465    */
0466   mutable unsigned int flow_;
0467
0468   /**
0469    * Diagram
0470    */
0471   mutable unsigned int diagram_;
0472   //@}
0473
0474   /**
0475    *  Storage of the matrix elements
0476    */
0477   //@{
0478   /**
0479    *  Matrix elements for the different colour flows
0480    */
0481   mutable vector<ProductionMatrixElement> flowME_;
0482
0483   /**
0484    *  Matrix elements for the different Feynman diagrams
0485    */
0486   mutable vector<ProductionMatrixElement> diagramME_;
0487   //@}
0488
0489 };
0490
0491 /** Exception class to indicate a problem has occurred with setting
0492     up to matrix element.*/
0493 class MEException : public Exception {};
0494
0495 }
0496
0497 #endif /* HERWIG_GeneralHardME_H */