LesHouchesReader.h

// -*- C++ -*-
//
// LesHouchesReader.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2007 Leif Lonnblad
//
// ThePEG is licenced under version 2 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 <cstdio>
#include <cstring>

namespace ThePEG {

class LesHouchesReader: public HandlerBase, public LastXCombInfo<> {

  friend class LesHouchesEventHandler;

  typedef FILE * CFile;

  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 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 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();

  bool compressedCache() const;

  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, &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(&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;

  double maxFactor;

  CrossSection weightScale;

  vector<double> xSecWeights;

  map<int,double> maxWeights;

  bool skipping;

  unsigned int theMomentumTreatment;

  bool useWeightWarnings;

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 &);

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

---