StandardXComb.h

// -*- C++ -*-
//
// StandardXComb.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
// Copyright (C) 2009-2019 Simon Platzer
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_StandardXComb_H
#define ThePEG_StandardXComb_H
// This is the declaration of the StandardXComb class.
 
#include "ThePEG/Config/ThePEG.h"
#include "SubProcessHandler.fh"
#include "ThePEG/PDF/PartonExtractor.fh"
#include "ThePEG/PDF/PartonBin.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "ThePEG/Utilities/VSelector.h"
#include "ThePEG/Utilities/ClassDescription.h"
#include "ThePEG/Utilities/Maths.h"
#include "ThePEG/Utilities/XSecStat.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/Handlers/XComb.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/SubProcessHandler.fh"
#include "StandardXComb.fh"
 
namespace ThePEG {
 
class StandardXComb: public XComb {
 
public:
 
  typedef MEBase::DiagramVector DiagramVector;
 
  typedef MEBase::DiagramIndex DiagramIndex;
 
  friend class MEBase;
 
public:
 
  StandardXComb(Energy newMaxEnergy, const cPDPair & inc,
                tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
                tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
                const PBPair & newPartonBins, tCutsPtr newCuts, tMEPtr newME,
                const DiagramVector & newDiagrams, bool mir,
                tStdXCombPtr newHead = tStdXCombPtr());
 
  StandardXComb(tStdXCombPtr newHead,
                const PBPair & newPartonBins, tMEPtr newME,
                const DiagramVector & newDiagrams);
 
  StandardXComb();
 
  virtual ~StandardXComb();
 
  StandardXComb(tMEPtr me, const tPVector & parts, DiagramIndex i);
 
 
 
  void recreatePartonBinInstances(Energy2 scale);
 
  void refillPartonBinInstances(const double* r);
 
  bool setIncomingPartons(tStdXCombPtr labHead = tStdXCombPtr());
 
  void fill(const PPair& newParticles,
            const PPair& newPartons,
            const vector<Lorentz5Momentum>& newMEMomenta,
            const DVector& newLastRandomNumbers = DVector());
 
 
  tcSubHdlPtr subProcessHandler() const { return theSubProcessHandler; }
 
  tMEPtr matrixElement() const { return theME; }
 
  tStdXCombPtr head() const { return theHead; }
 
  void head(tStdXCombPtr headXC) { theHead = headXC; }
 
  Selector<tStdXCombPtr>& projectors() { return theProjectors; }
 
  const Selector<tStdXCombPtr>& projectors() const { return theProjectors; }
 
  tStdXCombPtr lastProjector() const { return theProjector; }
 
  void lastProjector(tStdXCombPtr pxc) { theProjector = pxc; }
 
 
  virtual bool checkInit();
 
  virtual int nDim() const { return theNDim; }
 
  const pair<int,int>& partonDimensions() const { return partonDims; }
 
  bool willPassCuts();
 
  double cutWeight() const { return theCutWeight; }
 
  virtual void clean();
 
  bool kinematicsGenerated() const  { return theKinematicsGenerated; }
 
  void didGenerateKinematics() { theKinematicsGenerated = true; }
 
  virtual CrossSection dSigDR(const pair<double,double> ll, int nr, const double * r);
 
  CrossSection dSigDR(const double * r);
 
  virtual map<string,double> generateOptionalWeights();
 
  double lastPDFWeight() const { return theLastPDFWeight; }
 
  CrossSection lastCrossSection() const { return theLastCrossSection; }
 
  void checkReshufflingNeeds();
 
  bool needsReshuffling() const { return theNeedsReshuffling; }
 
  void reshuffle(vector<Lorentz5Momentum>&) const;
 
  virtual tSubProPtr construct();
 
  virtual const XSecStat & stats() const { return theStats; }
 
  virtual void select(double weight) { theStats.select(weight); }
 
  virtual void accept() { theStats.accept(); }
 
  void reweight(double oldWeight, double newWeight) {
    theStats.reweight(oldWeight,newWeight);
  }
 
  virtual void reject(double weight = 1.0) { theStats.reject(weight); }
 
  virtual void reset() { theStats.reset(); }
 
  const DiagramVector & diagrams() const { return theDiagrams; }
 
  bool mirror() const { return isMirror; }
 
  const vector<Lorentz5Momentum> & meMomenta() const { return theMEMomenta; }
 
  tcDiagPtr lastDiagram() const { 
    if ( !theExternalDiagram )
      return diagrams()[lastDiagramIndex()];
    return theExternalDiagram;
  }
 
  const cPDVector & mePartonData() const { return theMEPartonData; }
 
  DiagramIndex lastDiagramIndex() const { return theLastDiagramIndex; }
 
  const DVector & meInfo() const { return theMEInfo; }
 
  void meInfo(const DVector & info) { theMEInfo = info; }
 
  const DVector& lastRandomNumbers() const { return theLastRandomNumbers; }
 
  double jacobian() const { return theLastJacobian; }
 
  double lastME2() const { return theLastME2; }
 
  double lastPreweight() const { return theLastPreweight; }
 
  CrossSection lastMECrossSection() const { return theLastMECrossSection; }
 
  double lastMEPDFWeight() const { return theLastMEPDFWeight; }
 
  double lastMECouplings() const { return theLastMECouplings; }
 
  virtual void newSubProcess(bool group = false);
 
  vector<Lorentz5Momentum> & meMomenta() { return theMEMomenta; }
 
  DVector& lastRandomNumbers() { return theLastRandomNumbers; }
 
  cPDVector & mePartonData() { return theMEPartonData; }
 
  void externalDiagram(tcDiagPtr diag) { theExternalDiagram = diag; }
 
  void lastDiagramIndex(DiagramIndex i) { theLastDiagramIndex = i; }
 
  void lastPDFWeight(double w) { theLastPDFWeight = w; }
 
  void lastCrossSection(CrossSection s) { theLastCrossSection = s; }
 
  void jacobian(double j) { theLastJacobian = j; }
 
  void lastME2(double v) { theLastME2 = v; }
 
  void lastPreweight(double w) { theLastPreweight = w; }
 
  void lastMECrossSection(CrossSection v) { theLastMECrossSection = v; }
 
  void lastMEPDFWeight(double v) { theLastMEPDFWeight = v; }
 
  void lastMECouplings(double v) { theLastMECouplings = v; }
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
private:
 
  tSubHdlPtr theSubProcessHandler;
 
  tMEPtr theME;
 
  XSecStat theStats;
 
  DiagramVector theDiagrams;
 
  bool isMirror;
 
  int theNDim;
 
protected:
 
  pair<int,int> partonDims;
 
private:
 
  bool theKinematicsGenerated;
 
  vector<Lorentz5Momentum> theMEMomenta;
 
  cPDVector theMEPartonData;
 
  tcDiagPtr theExternalDiagram;
 
  DiagramIndex theLastDiagramIndex;
 
  DVector theMEInfo;
 
  DVector theLastRandomNumbers;
 
  double theLastPDFWeight;
 
  CrossSection theLastCrossSection;
 
  double theLastJacobian;
 
  double theLastME2;
 
  double theLastPreweight;
 
  CrossSection theLastMECrossSection;
 
  double theLastMEPDFWeight;
 
  double theLastMECouplings;
 
  tStdXCombPtr theHead;
 
  Selector<tStdXCombPtr> theProjectors;
 
  tStdXCombPtr theProjector;
 
  bool checkedCuts;
 
  bool passedCuts;
 
  double theCutWeight;
 
  bool theNeedsReshuffling;
 
  double reshuffleEquation(double, const vector<pair<Energy2,Energy2> >&, Energy2) const;
 
  double solveReshuffleEquation(const vector<pair<Energy2,Energy2> >&, Energy2) const;
 
private:
 
  static ClassDescription<StandardXComb> initStandardXComb;
 
  StandardXComb & operator=(const StandardXComb &) = delete;
 
};
 
template <>
struct BaseClassTrait<StandardXComb,1>: public ClassTraitsType {
  typedef XComb NthBase;
};
 
template <>
struct ClassTraits<StandardXComb>:
    public ClassTraitsBase<StandardXComb> {
  static string className() { return "ThePEG::StandardXComb"; }
};
 
}
 
#endif /* ThePEG_StandardXComb_H */