LesHouchesReader.h

// -*- C++ -*-
//
// LesHouchesReader.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// 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_LesHouchesReader_H
#define THEPEG_LesHouchesReader_H
// This is the declaration of the LesHouchesReader class.
 
#include "LesHouches.h"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Utilities/ObjectIndexer.h"
#include "ThePEG/Utilities/Exception.h"
#include "ThePEG/Utilities/XSecStat.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "ThePEG/PDF/PartonBin.fh"
#include "ThePEG/MatrixElement/ReweightBase.h"
#include "LesHouchesEventHandler.fh"
#include "LesHouchesReader.fh"
#include "ThePEG/Utilities/CFile.h"
#include <cstdio>
#include <cstring>
 
namespace ThePEG {
 
class LesHouchesReader: public HandlerBase, public LastXCombInfo<> {
 
  friend class LesHouchesEventHandler;
 
  typedef map<int,XSecStat> StatMap;
 
  typedef map<tcPBPair,XCombPtr> XCombMap;
 
  typedef vector<ReweightPtr> ReweightVector;
 
public:
 
  LesHouchesReader(bool active = false);
 
  LesHouchesReader(const LesHouchesReader &);
 
  virtual ~LesHouchesReader();
 
public:
 
  virtual void open() = 0;
 
  virtual bool doReadEvent() = 0;
 
  virtual void close() = 0;
 
  virtual vector<string> optWeightsNamesFunc() = 0;
 
  vector<string> optionalWeightsNames;
 
 
  virtual void initialize(LesHouchesEventHandler & eh);
 
  virtual double getEvent();
 
  virtual bool readEvent();
 
  virtual void skip(long n);
 
  tXCombPtr getXComb();
 
  tSubProPtr getSubProcess();
 
  virtual long scan();
 
  virtual void initStat();
 
  double reweight();
 
  virtual void fillEvent();
 
  void reset();
 
  virtual void setWeightScale(long neve);
 
 
 
  static size_t eventSize(int N) {
    return (N + 1)*sizeof(int) +       // IDPRUP, ISTUP
      (7*N + 4)*sizeof(double) +       // XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP,
                                       // VTIMUP, SPINUP
      N*sizeof(long) +                 // IDUP
      2*N*sizeof(pair<int,int>) +      // MOTHUP, ICOLUP
      sizeof(pair<double,double>) +    // XPDWUP.
      2*sizeof(double);                // lastweight and preweight
  }
 
  double eventWeight() const { return hepeup.XWGTUP*lastweight; }
 
  const map<string,double>& optionalEventWeights() const { return optionalWeights; }
 
  const long& LHEEventNum() const { return LHEeventnum; }
 
  const int& optionalEventnpLO() const { return optionalnpLO; }
  const int& optionalEventnpNLO() const { return optionalnpNLO; }
  
  const PBIPair & partonBinInstances() const { return thePartonBinInstances; }
  const PPair & beams() const { return theBeams; }
  const PPair & incoming() const { return theIncoming; }
  const PVector & outgoing() const { return theOutgoing; }
  const PVector & intermediates() const { return theIntermediates; }
  int maxMultCKKW() const { return theMaxMultCKKW; }
  int minMultCKKW() const { return theMinMultCKKW; }  //@}
 
  long NEvents() const { return theNEvents; }
 
  long currentPosition() const { return position; }
 
  long maxScan() const { return theMaxScan; }
 
  bool active() const { return isActive; }
 
  bool negativeWeights() const { return heprup.IDWTUP < 0; }
 
  const XSecStat & xSecStats() const { return stats; }
 
  const StatMap & processStats() const { return statmap; }
 
  void select(double weight) {
    stats.select(weight);
    statmap[hepeup.IDPRUP].select(weight);
  }
 
  void accept() {
    stats.accept();
    statmap[hepeup.IDPRUP].accept();
  }
 
  void reject(double w) {
    stats.reject(w);
    statmap[hepeup.IDPRUP].reject(w);
  }
 
  virtual void increaseMaxXSec(CrossSection maxxsec);
 
  tPExtrPtr partonExtractor() const { return thePartonExtractor; }
 
  tCascHdlPtr CKKWHandler() const { return theCKKW; }
 
  const PartonPairVec & partonBins() const { return thePartonBins; }
 
  const XCombMap & xCombs() const { return theXCombs; }
 
  const Cuts & cuts() const { return *theCuts; }
 
 
protected:
 
 
  string cacheFileName() const { return theCacheFileName; }
 
  bool cutEarly() const { return doCutEarly; }
 
  CFile cacheFile() const { return theCacheFile;}
 
  void openReadCacheFile();
 
  void openWriteCacheFile();
 
  void closeCacheFile();
 
  void cacheEvent() const;
 
  bool uncacheEvent();
 
  void reopen();
 
  template <typename T>
  static char * mwrite(char * pos, const T & t, size_t n = 1) {
    std::memcpy(pos, (void *)&t, n*sizeof(T));
    return pos + n*sizeof(T);
  }
 
  template <typename T>
  static const char * mread(const char * pos, T & t, size_t n = 1) {
    std::memcpy((void *)&t, pos, n*sizeof(T));
    return pos + n*sizeof(T);
  }
 
 
  virtual bool checkPartonBin();
 
  virtual void createParticles();
 
  virtual tcPBPair createPartonBinInstances();
 
  virtual void createBeams();
 
  virtual void connectMothers();
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  void NEvents(long x) { theNEvents = x; }
 
  XCombMap & xCombs() { return theXCombs; }  
 
  virtual void doinit();
 
  virtual void doinitrun();
 
  virtual void dofinish() {
    close();
    HandlerBase::dofinish();
  }
 
  virtual bool preInitialize() const;
 
  virtual void initPDFs();
 
protected:
 
  HEPRUP heprup;
 
  HEPEUP hepeup;
 
  tcPDPair inData;
 
  pair<PDFPtr,PDFPtr> inPDF;
 
  pair<cPDFPtr,cPDFPtr> outPDF;
 
  PExtrPtr thePartonExtractor;
 
  tCascHdlPtr theCKKW;
 
  PartonPairVec thePartonBins;
 
  XCombMap theXCombs;
 
  CutsPtr theCuts;
 
  long theNEvents;
 
  long position;
 
  int reopened;
 
  long theMaxScan;
 
  bool scanning;
 
  bool isActive;
 
  string theCacheFileName;
 
  bool doCutEarly;
 
  XSecStat stats;
 
  StatMap statmap;
 
  PBIPair thePartonBinInstances;
 
  ObjectIndexer<long,ColourLine> colourIndex;
 
  ObjectIndexer<long,Particle> particleIndex;
 
  PPair theBeams;
 
  PPair theIncoming;
 
  PVector theOutgoing;
 
  PVector theIntermediates;
 
  CFile theCacheFile;
 
  ReweightVector reweights;
 
  ReweightVector preweights;
 
  double preweight;
 
  bool reweightPDF;
 
  bool doInitPDFs;
 
  int theMaxMultCKKW;
 
  int theMinMultCKKW;
 
  double lastweight;
 
  map<string,double> optionalWeights;
 
  long LHEeventnum;
 
 
  double maxFactor;
 
  int optionalnpLO;
 
  int optionalnpNLO;
 
  CrossSection weightScale;
 
  vector<double> xSecWeights;
 
  map<int,double> maxWeights;
 
  bool skipping;
 
  unsigned int theMomentumTreatment;
 
  bool useWeightWarnings;
 
  bool theReOpenAllowed;
 
  bool theIncludeSpin;
 
private:
 
  void setBeamA(long id);
  long getBeamA() const;
  void setBeamB(long id);
  long getBeamB() const;
  void setEBeamA(Energy e);
  Energy getEBeamA() const;
  void setEBeamB(Energy e);
  Energy getEBeamB() const;
  void setPDFA(PDFPtr);
  PDFPtr getPDFA() const;
  void setPDFB(PDFPtr);
  PDFPtr getPDFB() const;
 
private:
 
  static AbstractClassDescription<LesHouchesReader> initLesHouchesReader;
 
  LesHouchesReader & operator=(const LesHouchesReader &) = delete;
 
public:
 
  class LesHouchesInconsistencyError: public Exception {};
  
  class LesHouchesReopenWarning: public Exception {};
 
  class LesHouchesReopenError: public Exception {};
 
  class LesHouchesInitError: public InitException {};
};
 
ostream & operator<<(ostream & os, const HEPEUP & h);
 
}
 
 
#include "ThePEG/Utilities/ClassTraits.h"
 
namespace ThePEG {
 
template <>
struct BaseClassTrait<LesHouchesReader,1>: public ClassTraitsType {
  typedef HandlerBase NthBase;
};
 
template <>
struct ClassTraits<LesHouchesReader>
  : public ClassTraitsBase<LesHouchesReader> {
  static string className() { return "ThePEG::LesHouchesReader"; }
  static string library() { return "LesHouches.so"; }
 
};
 
}
 
#endif /* THEPEG_LesHouchesReader_H */