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 */