thepeg is hosted by Hepforge, IPPP Durham
ThePEG  2.2.1
LesHouchesReader.h
1 // -*- C++ -*-
2 //
3 // LesHouchesReader.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2019 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef THEPEG_LesHouchesReader_H
10 #define THEPEG_LesHouchesReader_H
11 // This is the declaration of the LesHouchesReader class.
12 
13 #include "LesHouches.h"
14 #include "ThePEG/Handlers/HandlerBase.h"
15 #include "ThePEG/Utilities/ObjectIndexer.h"
16 #include "ThePEG/Utilities/Exception.h"
17 #include "ThePEG/Utilities/XSecStat.h"
18 #include "ThePEG/PDF/PartonBinInstance.h"
19 #include "ThePEG/PDF/PartonBin.fh"
20 #include "ThePEG/MatrixElement/ReweightBase.h"
21 #include "LesHouchesEventHandler.fh"
22 #include "LesHouchesReader.fh"
23 #include "ThePEG/Utilities/CFile.h"
24 #include <cstdio>
25 #include <cstring>
26 
27 namespace ThePEG {
28 
77 class LesHouchesReader: public HandlerBase, public LastXCombInfo<> {
78 
82  friend class LesHouchesEventHandler;
83 
88  typedef map<int,XSecStat> StatMap;
89 
94  typedef map<tcPBPair,XCombPtr> XCombMap;
95 
99  typedef vector<ReweightPtr> ReweightVector;
100 
101 public:
102 
110  LesHouchesReader(bool active = false);
111 
116 
120  virtual ~LesHouchesReader();
122 
123 public:
124 
134  virtual void open() = 0;
135 
141  virtual bool doReadEvent() = 0;
142 
146  virtual void close() = 0;
147 
151  virtual vector<string> optWeightsNamesFunc() = 0;
152 
157  vector<string> optionalWeightsNames;
158 
160 
170  virtual void initialize(LesHouchesEventHandler & eh);
171 
184  virtual double getEvent();
185 
191  virtual bool readEvent();
192 
198  virtual void skip(long n);
199 
207 
213 
220  virtual long scan();
221 
226  virtual void initStat();
227 
233  double reweight();
234 
239  virtual void fillEvent();
240 
245  void reset();
246 
254  virtual void setWeightScale(long neve);
255 
257 
260 
266  static size_t eventSize(int N) {
267  return (N + 1)*sizeof(int) + // IDPRUP, ISTUP
268  (7*N + 4)*sizeof(double) + // XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP,
269  // VTIMUP, SPINUP
270  N*sizeof(long) + // IDUP
271  2*N*sizeof(pair<int,int>) + // MOTHUP, ICOLUP
272  sizeof(pair<double,double>) + // XPDWUP.
273  2*sizeof(double); // lastweight and preweight
274  }
275 
282  double eventWeight() const { return hepeup.XWGTUP*lastweight; }
283 
287  const map<string,double>& optionalEventWeights() const { return optionalWeights; }
288 
292  const int& optionalEventnpLO() const { return optionalnpLO; }
293  const int& optionalEventnpNLO() const { return optionalnpNLO; }
294 
303  const PPair & beams() const { return theBeams; }
308  const PPair & incoming() const { return theIncoming; }
313  const PVector & outgoing() const { return theOutgoing; }
318  const PVector & intermediates() const { return theIntermediates; }
325  int maxMultCKKW() const { return theMaxMultCKKW; }
332  int minMultCKKW() const { return theMinMultCKKW; } //@}
333 
340  long NEvents() const { return theNEvents; }
341 
346  long currentPosition() const { return position; }
347 
353  long maxScan() const { return theMaxScan; }
354 
358  bool active() const { return isActive; }
359 
363  bool negativeWeights() const { return heprup.IDWTUP < 0; }
364 
368  const XSecStat & xSecStats() const { return stats; }
369 
373  const StatMap & processStats() const { return statmap; }
374 
379  void select(double weight) {
380  stats.select(weight);
381  statmap[hepeup.IDPRUP].select(weight);
382  }
383 
387  void accept() {
388  stats.accept();
389  statmap[hepeup.IDPRUP].accept();
390  }
391 
395  void reject(double w) {
396  stats.reject(w);
397  statmap[hepeup.IDPRUP].reject(w);
398  }
399 
403  virtual void increaseMaxXSec(CrossSection maxxsec);
404 
409 
414  tCascHdlPtr CKKWHandler() const { return theCKKW; }
415 
420  const PartonPairVec & partonBins() const { return thePartonBins; }
421 
426  const XCombMap & xCombs() const { return theXCombs; }
427 
431  const Cuts & cuts() const { return *theCuts; }
432 
434 
435 protected:
436 
439 
444  string cacheFileName() const { return theCacheFileName; }
445 
450  bool cutEarly() const { return doCutEarly; }
451 
455  CFile cacheFile() const { return theCacheFile;}
456 
460  void openReadCacheFile();
461 
465  void openWriteCacheFile();
466 
470  void closeCacheFile();
471 
475  void cacheEvent() const;
476 
480  bool uncacheEvent();
481 
487  void reopen();
488 
492  template <typename T>
493  static char * mwrite(char * pos, const T & t, size_t n = 1) {
494  std::memcpy(pos, &t, n*sizeof(T));
495  return pos + n*sizeof(T);
496  }
497 
501  template <typename T>
502  static const char * mread(const char * pos, T & t, size_t n = 1) {
503  std::memcpy(&t, pos, n*sizeof(T));
504  return pos + n*sizeof(T);
505  }
506 
508 
517  virtual bool checkPartonBin();
518 
523  virtual void createParticles();
524 
530  virtual tcPBPair createPartonBinInstances();
531 
537  virtual void createBeams();
538 
542  virtual void connectMothers();
544 
545 public:
546 
553  void persistentOutput(PersistentOStream & os) const;
554 
560  void persistentInput(PersistentIStream & is, int version);
562 
566  static void Init();
567 
568 protected:
569 
576  void NEvents(long x) { theNEvents = x; }
577 
582  XCombMap & xCombs() { return theXCombs; }
584 
592  virtual void doinit();
593 
598  virtual void doinitrun();
599 
604  virtual void dofinish() {
605  close();
607  }
608 
613  virtual bool preInitialize() const;
614 
619  virtual void initPDFs();
621 
622 protected:
623 
628 
633 
638 
644  pair<PDFPtr,PDFPtr> inPDF;
645 
649  pair<cPDFPtr,cPDFPtr> outPDF;
650 
655 
660 
665  PartonPairVec thePartonBins;
666 
671  XCombMap theXCombs;
672 
676  CutsPtr theCuts;
677 
683 
688  long position;
689 
693  int reopened;
694 
701 
705  bool scanning;
706 
710  bool isActive;
711 
717 
723 
728 
732  StatMap statmap;
733 
739 
745 
751 
756 
762 
768 
774 
779 
783  ReweightVector reweights;
784 
788  ReweightVector preweights;
789 
793  double preweight;
794 
799 
805 
813 
821 
827  double lastweight;
828 
832  map<string,double> optionalWeights;
833 
834 
840  double maxFactor;
841 
846 
851 
857 
861  vector<double> xSecWeights;
862 
867  map<int,double> maxWeights;
868 
872  bool skipping;
873 
877  unsigned int theMomentumTreatment;
878 
884 
889 
894 
895 private:
896 
898  void setBeamA(long id);
900  long getBeamA() const;
902  void setBeamB(long id);
904  long getBeamB() const;
906  void setEBeamA(Energy e);
908  Energy getEBeamA() const;
910  void setEBeamB(Energy e);
912  Energy getEBeamB() const;
914  void setPDFA(PDFPtr);
916  PDFPtr getPDFA() const;
918  void setPDFB(PDFPtr);
920  PDFPtr getPDFB() const;
921 
922 private:
923 
928 
932  LesHouchesReader & operator=(const LesHouchesReader &) = delete;
933 
934 public:
935 
939  class LesHouchesInconsistencyError: public Exception {};
940 
943  class LesHouchesReopenWarning: public Exception {};
944 
947  class LesHouchesReopenError: public Exception {};
948 
951  class LesHouchesInitError: public InitException {};
954 };
955 
957 ostream & operator<<(ostream & os, const HEPEUP & h);
958 
959 }
960 
961 
962 #include "ThePEG/Utilities/ClassTraits.h"
963 
964 namespace ThePEG {
965 
972 template <>
973 struct BaseClassTrait<LesHouchesReader,1>: public ClassTraitsType {
975  typedef HandlerBase NthBase;
976 };
977 
983 template <>
984 struct ClassTraits<LesHouchesReader>
985  : public ClassTraitsBase<LesHouchesReader> {
989  static string className() { return "ThePEG::LesHouchesReader"; }
995  static string library() { return "LesHouches.so"; }
996 
997 };
998 
1001 }
1002 
1003 #endif /* THEPEG_LesHouchesReader_H */
map< string, double > optionalWeights
The optional weights associated to the last read events.
pair< tcPDPtr, tcPDPtr > tcPDPair
A pair of transient pointers to const ParticleData objects.
Definition: Containers.h:124
long theMaxScan
The maximum number of events to scan to collect information about processes and cross sections...
tCascHdlPtr CKKWHandler() const
Return a possibly null pointer to a CascadeHandler to be used for CKKW-reweighting.
PersistentIStream is used to read persistent objects from a stream where they were previously written...
void select(double weight)
An event of the corresponding class has been attempted.
Definition: XSecStat.h:117
map< int, double > maxWeights
Individual maximum weights for individual (possibly reweighted) processes.
bool active() const
Return true if this reader is active.
long getBeamA() const
Access function for the interface.
vector< double > xSecWeights
Individual scales for different sub-processes if reweighted.
virtual ~LesHouchesReader()
Destructor.
void reset()
Removes the particles created in the last generated event, preparing to produce a new one...
XCombMap & xCombs()
The map of XComb objects indexed by the corresponding PartonBin pair.
static const char * mread(const char *pos, T &t, size_t n=1)
Helper function to read a variable from a memory location.
The LesHouchesEventHandler inherits from the general EventHandler class and administers the reading o...
long maxScan() const
The maximum number of events to scan to collect information about processes and cross sections...
const PPair & beams() const
Return the instances of the beam particles for the current event.
static AbstractClassDescription< LesHouchesReader > initLesHouchesReader
Describe an abstract base class with persistent data.
ClassTraitsType is an empty, non-polymorphic, base class.
Definition: ClassTraits.h:30
string theCacheFileName
Name of file used to cache the events form the reader in a fast-readable form.
PExtrPtr thePartonExtractor
The PartonExtractor object used to construct remnants.
void accept()
An event of the corresponding class has been accepted.
Definition: XSecStat.h:109
bool isActive
True if this is an active reader.
virtual void initialize(LesHouchesEventHandler &eh)
Initialize.
int maxMultCKKW() const
If this reader is to be used (possibly together with others) for CKKW reweighting and veto...
XCombMap theXCombs
The map of XComb objects indexed by the corresponding PartonBin pair.
virtual void open()=0
Open a file or stream with events and read in the run information into the heprup variable...
virtual void fillEvent()
Converts the information in the Les Houches common block variables into a Particle objects...
static char * mwrite(char *pos, const T &t, size_t n=1)
Helper function to write a variable to a memory location.
virtual void skip(long n)
Skip n events.
void setEBeamB(Energy e)
Access function for the interface.
const map< string, double > & optionalEventWeights() const
Return the optional named weights associated to the current event.
long getBeamB() const
Access function for the interface.
const PBIPair & partonBinInstances() const
The pair of PartonBinInstance objects describing the current incoming partons in the event...
tPExtrPtr partonExtractor() const
The PartonExtractor object used to construct remnants.
PersistentOStream is used to write objects persistently to a stream from which they can be read in ag...
const PPair & incoming() const
Return the instances of the incoming particles to the sub process for the current event...
A concreate implementation of ClassDescriptionBase describing an abstract class with persistent data...
LastXCombInfo is a templated class giving easy access to the information in an XComb object...
Definition: LastXCombInfo.h:32
virtual void createBeams()
Create instances of the incoming beams in the event and store them in particleIndex.
unsigned int theMomentumTreatment
Option for the treatment of the momenta supplied.
const StatMap & processStats() const
Collected statistics about the individual processes.
bool doInitPDFs
Should PDFBase objects be constructed from the information in the event file in the initialization...
bool uncacheEvent()
Read an event from the cache file.
void select(double weight)
Select the current event.
TransientRCPtr is a simple wrapper around a bare pointer which can be assigned to and from an RCPtr a...
Definition: RCPtr.h:519
virtual vector< string > optWeightsNamesFunc()=0
return the weight names
const PVector & outgoing() const
Return the instances of the outgoing particles from the sub process for the current event...
bool reweightPDF
Should the event be reweighted by PDFs used by the PartonExtractor?
const XCombMap & xCombs() const
The map of XComb objects indexed by the corresponding PartonBin pair.
virtual void doinit()
Initialize this object after the setup phase before saving an EventGenerator to disk.
bool skipping
Is set to true when getEvent() is called from skip(int).
bool theIncludeSpin
Use the spin information.
virtual void close()=0
Close the file or stream from which events have been read.
int IDWTUP
Master switch indicating how the ME generator envisages the events weights should be interpreted acco...
Definition: LesHouches.h:89
virtual bool checkPartonBin()
Check the existence of a pair of PartonBin objects corresponding to the current event.
LesHouchesReader is an abstract base class to be used for objects which reads event files or streams ...
LesHouchesReader & operator=(const LesHouchesReader &)=delete
Private and non-existent assignment operator.
PBIPair thePartonBinInstances
The pair of PartonBinInstance objects describing the current incoming partons in the event...
virtual bool preInitialize() const
Return true if this object needs to be initialized before all other objects because it needs to extra...
void setBeamB(long id)
Access function for the interface.
This is the main namespace within which all identifiers in ThePEG are declared.
Definition: FactoryBase.h:28
virtual double getEvent()
Calls readEvent() or uncacheEvent() to read information into the LesHouches common block variables...
int theMinMultCKKW
If this reader is to be used (possibly together with others) for CKKW reweighting and veto...
virtual void doinitrun()
Initialize this object.
PVector theOutgoing
The instances of the outgoing particles from the sub process for the current event.
virtual void initPDFs()
Called from doinit() to extract PDFs from the event file and add the corresponding objects to the cur...
bool cutEarly() const
Determines whether to apply cuts to events converting them to ThePEG format.
void reject(double weight=1.0)
Reject the event which was last accepted with accept() or selected with select(double).
Definition: XSecStat.h:142
StatMap statmap
Collect statistics for each individual process.
long currentPosition() const
The number of events produced so far.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
virtual void dofinish()
Finalize this object.
string cacheFileName() const
Name of file used to cache the events form the reader in a fast-readable form.
void reject(double w)
Reject the current event assuming it was previously accepted.
ObjectIndexer< long, Particle > particleIndex
Association between Particles and indices in the current translation.
CrossSection weightScale
The (reweighted) XWGTUP value should be scaled with this cross section when compared to the overestim...
ReweightVector reweights
The reweight objects modifying the weights of this reader.
bool theReOpenAllowed
Option to allow reopening of the file.
virtual tcPBPair createPartonBinInstances()
Using the already created particles create a pair of PartonBinInstance objects corresponding to the i...
HEPRUP heprup
The HEPRUP common block.
virtual void connectMothers()
Go through the mother indices and connect up the Particles.
PVector theIntermediates
The instances of the intermediate particles in the sub process for the current event.
PPair theIncoming
The instances of the incoming particles to the sub process for the current event. ...
Here is the documentation of the CFile class.
Definition: CFile.h:15
void setPDFB(PDFPtr)
Access function for the interface.
int theMaxMultCKKW
If this reader is to be used (possibly together with others) for CKKW reweighting and veto...
LesHouchesReader(bool active=false)
Default constructor.
virtual void initStat()
Take the information corresponding to the HEPRUP common block and initialize the statistics for this ...
int reopened
The number of times this reader has been reopened.
long NEvents() const
The number of events found in this reader.
int optionalnpLO
npLO for LesHouches merging
pair< PPtr, PPtr > PPair
A pair of pointers to Particle objects.
Definition: Containers.h:127
pair< PBIPtr, PBIPtr > PBIPair
A pair of pointers to PartonBinInstance objects.
static size_t eventSize(int N)
Return the size of this event in bytes.
ObjectIndexer< long, ColourLine > colourIndex
Association between ColourLines and colour indices in the current translation.
RCPtr is a reference counted (smart) pointer.
Definition: RCPtr.h:60
virtual void dofinish()
Finalize this object.
map< int, XSecStat > StatMap
Map for accumulating statistics of cross sections per process number.
bool scanning
Flag to tell whether we are in the process of scanning.
bool negativeWeights() const
True if negative weights may be produced.
virtual void createParticles()
Create instances of all particles in the event and store them in particleIndex.
const PartonPairVec & partonBins() const
The pairs of PartonBin objects describing the partons which can be extracted by the PartonExtractor o...
bool useWeightWarnings
Set to true if warnings about possible weight incompatibilities should be issued. ...
vector< string > optionalWeightsNames
vector with the optional weights names
PDFPtr getPDFB() const
Access function for the interface.
HandlerBase is an abstract base class derived from the Interfaced class via the HandlerBaseT class ad...
Definition: HandlerBase.h:151
void setPDFA(PDFPtr)
Access function for the interface.
tcPDPair inData
The ParticleData objects corresponding to the incoming particles.
double eventWeight() const
The current event weight given by XWGTUP times possible reweighting.
Exception is the base class for all exceptions to be used in ThePEG.
Definition: Exception.h:44
double maxFactor
If the maximum cross section of this reader has been increased with increaseMaxXSec(), this is the total factor with which it has been increased.
int optionalnpNLO
npNLO for LesHouches merging
void setBeamA(long id)
Access function for the interface.
The HEPEUP class is a simple container corresponding to the Les Houches accord (hep-ph/0109068) commo...
Definition: LesHouches.h:128
int IDPRUP
The subprocess code for this event (as given in LPRUP).
Definition: LesHouches.h:177
vector< ReweightPtr > ReweightVector
A vector of pointers to ReweightBase objects.
The HEPRUP class is a simple container corresponding to the Les Houches accord (hep-ph/0109068) commo...
Definition: LesHouches.h:26
void setEBeamA(Energy e)
Access function for the interface.
tCascHdlPtr theCKKW
A pointer to a CascadeHandler to be used for CKKW-reweighting.
long position
The number of events produced by this reader so far.
virtual void setWeightScale(long neve)
Possibility for subclasses to recover from non-conformant settings of XMAXUP when an event file has b...
CFile theCacheFile
File stream for the cache.
map< tcPBPair, XCombPtr > XCombMap
Map of XComb objects describing the incoming partons indexed by the corresponding PartonBin pair...
The default concrete implementation of ClassTraitsBase.
Definition: ClassTraits.h:134
This is a templated class which dynamically associates (reference counted) objects to integer indices...
Definition: ObjectIndexer.h:25
PDFPtr getPDFA() const
Access function for the interface.
virtual bool readEvent()
Calls doReadEvent() and performs pre-defined reweightings.
ReweightVector preweights
The preweight objects modifying the weights of this reader.
vector< T > & operator<<(vector< T > &tv, const U &u)
Overload the left shift operator for vector to push_back objects to a vector.
Definition: Containers.h:179
tXCombPtr getXComb()
Get an XComb object.
const XSecStat & xSecStats() const
The collected cross section statistics for this reader.
pair< PDFPtr, PDFPtr > inPDF
The PDFBase objects which has been used for the beam particle when generating the events being read...
PPair theBeams
The instances of the beam particles for the current event.
virtual bool doReadEvent()=0
Read the next event from the file or stream into the corresponding protected variables.
void reopen()
Reopen a reader.
bool doCutEarly
Determines whether to apply cuts to events before converting them to ThePEG format.
pair< cPDFPtr, cPDFPtr > outPDF
The PDFBase object to be used in the subsequent generation.
Energy getEBeamB() const
Access function for the interface.
virtual long scan()
Scan the file or stream to obtain information about cross section weights and particles etc...
void openReadCacheFile()
Open the cache file for reading.
void openWriteCacheFile()
Open the cache file for writing.
double preweight
The factor with which this reader was last pre-weighted.
Energy getEBeamA() const
Access function for the interface.
const int & optionalEventnpLO() const
Return the optional npLO and npNLO.
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
int minMultCKKW() const
If this reader is to be used (possibly together with others) for CKKW reweighting and veto...
void NEvents(long x)
The number of events in this reader.
void accept()
Accept the current event assuming it was previously selcted.
vector< PPtr > PVector
A vector of pointers to Particle objects.
Definition: Containers.h:76
XSecStat stats
Collect statistics for this reader.
tSubProPtr getSubProcess()
Get a SubProcess object corresponding to the information in the Les Houches common block variables...
double XWGTUP
The weight for this event.
Definition: LesHouches.h:182
BaseClassTraits describes the base classes of the templated class.
Definition: ClassTraits.h:156
XSecStat is a concrete helper class used to collect statistics about the cross section for a specific...
Definition: XSecStat.h:36
long theNEvents
The number of events in this reader.
void cacheEvent() const
Write the current event to the cache file.
CFile cacheFile() const
File stream for the cache.
double lastweight
The weight multiplying the last read event due to PDF reweighting, CKKW reweighting or assigned rewei...
static void Init()
Standard Init function used to initialize the interfaces.
const PVector & intermediates() const
Return the instances of the intermediate particles in the sub process for the current event...
PartonPairVec thePartonBins
The pairs of PartonBin objects describing the partons which can be extracted by the PartonExtractor o...
Cuts is a class for implementing kinematical cuts in ThePEG.
Definition: Cuts.h:53
virtual void increaseMaxXSec(CrossSection maxxsec)
Increase the overestimated cross section for this reader.
double reweight()
Reweights the current event using the reweights and preweights vectors.
CutsPtr theCuts
The Cuts object to be used for this reader.
The templated ClassTraitsBase class defines a set of default information about classes used by ThePEG...
Definition: ClassTraits.h:52
const Cuts & cuts() const
The Cuts object to be used for this reader.
HEPEUP hepeup
The HEPEUP common block.
void closeCacheFile()
Close the cache file;.