ThePEG  1.8.0
LesHouchesReader.h
1 // -*- C++ -*-
2 //
3 // LesHouchesReader.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2011 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 2 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;
148 
158  virtual void initialize(LesHouchesEventHandler & eh);
159 
172  virtual double getEvent();
173 
179  virtual bool readEvent();
180 
186  virtual void skip(long n);
187 
195 
201 
208  virtual long scan();
209 
214  virtual void initStat();
215 
221  double reweight();
222 
227  virtual void fillEvent();
228 
233  void reset();
234 
242  virtual void setWeightScale(long neve);
243 
245 
248 
254  static size_t eventSize(int N) {
255  return (N + 1)*sizeof(int) + // IDPRUP, ISTUP
256  (7*N + 4)*sizeof(double) + // XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP,
257  // VTIMUP, SPINUP
258  N*sizeof(long) + // IDUP
259  2*N*sizeof(pair<int,int>) + // MOTHUP, ICOLUP
260  sizeof(pair<double,double>) + // XPDWUP.
261  2*sizeof(double); // lastweight and preweight
262  }
263 
270  double eventWeight() const { return hepeup.XWGTUP*lastweight; }
271 
275  const map<string,double>& optionalEventWeights() const { return optionalWeights; }
276 
285  const PPair & beams() const { return theBeams; }
290  const PPair & incoming() const { return theIncoming; }
295  const PVector & outgoing() const { return theOutgoing; }
300  const PVector & intermediates() const { return theIntermediates; }
307  int maxMultCKKW() const { return theMaxMultCKKW; }
314  int minMultCKKW() const { return theMinMultCKKW; } //@}
315 
322  long NEvents() const { return theNEvents; }
323 
328  long currentPosition() const { return position; }
329 
335  long maxScan() const { return theMaxScan; }
336 
340  bool active() const { return isActive; }
341 
345  bool negativeWeights() const { return heprup.IDWTUP < 0; }
346 
350  const XSecStat & xSecStats() const { return stats; }
351 
355  const StatMap & processStats() const { return statmap; }
356 
361  void select(double weight) {
362  stats.select(weight);
363  statmap[hepeup.IDPRUP].select(weight);
364  }
365 
369  void accept() {
370  stats.accept();
371  statmap[hepeup.IDPRUP].accept();
372  }
373 
377  void reject(double w) {
378  stats.reject(w);
379  statmap[hepeup.IDPRUP].reject(w);
380  }
381 
385  virtual void increaseMaxXSec(CrossSection maxxsec);
386 
391 
396  tCascHdlPtr CKKWHandler() const { return theCKKW; }
397 
402  const PartonPairVec & partonBins() const { return thePartonBins; }
403 
408  const XCombMap & xCombs() const { return theXCombs; }
409 
413  const Cuts & cuts() const { return *theCuts; }
414 
416 
417 protected:
418 
421 
426  string cacheFileName() const { return theCacheFileName; }
427 
432  bool cutEarly() const { return doCutEarly; }
433 
437  CFile cacheFile() const { return theCacheFile;}
438 
442  void openReadCacheFile();
443 
447  void openWriteCacheFile();
448 
452  void closeCacheFile();
453 
457  void cacheEvent() const;
458 
462  bool uncacheEvent();
463 
469  void reopen();
470 
474  template <typename T>
475  static char * mwrite(char * pos, const T & t, size_t n = 1) {
476  std::memcpy(pos, &t, n*sizeof(T));
477  return pos + n*sizeof(T);
478  }
479 
483  template <typename T>
484  static const char * mread(const char * pos, T & t, size_t n = 1) {
485  std::memcpy(&t, pos, n*sizeof(T));
486  return pos + n*sizeof(T);
487  }
488 
490 
499  virtual bool checkPartonBin();
500 
505  virtual void createParticles();
506 
512  virtual tcPBPair createPartonBinInstances();
513 
519  virtual void createBeams();
520 
524  virtual void connectMothers();
526 
527 public:
528 
535  void persistentOutput(PersistentOStream & os) const;
536 
542  void persistentInput(PersistentIStream & is, int version);
544 
548  static void Init();
549 
550 protected:
551 
558  void NEvents(long x) { theNEvents = x; }
559 
564  XCombMap & xCombs() { return theXCombs; }
566 
574  virtual void doinit();
575 
580  virtual void doinitrun();
581 
586  virtual void dofinish() {
587  close();
589  }
590 
595  virtual bool preInitialize() const;
596 
601  virtual void initPDFs();
603 
604 protected:
605 
610 
615 
620 
626  pair<PDFPtr,PDFPtr> inPDF;
627 
631  pair<cPDFPtr,cPDFPtr> outPDF;
632 
637 
642 
647  PartonPairVec thePartonBins;
648 
653  XCombMap theXCombs;
654 
658  CutsPtr theCuts;
659 
665 
670  long position;
671 
675  int reopened;
676 
683 
687  bool scanning;
688 
692  bool isActive;
693 
699 
705 
710 
714  StatMap statmap;
715 
721 
727 
733 
738 
744 
750 
756 
761 
765  ReweightVector reweights;
766 
770  ReweightVector preweights;
771 
775  double preweight;
776 
781 
787 
795 
803 
809  double lastweight;
810 
814  map<string,double> optionalWeights;
815 
821  double maxFactor;
822 
828 
832  vector<double> xSecWeights;
833 
838  map<int,double> maxWeights;
839 
843  bool skipping;
844 
848  unsigned int theMomentumTreatment;
849 
855 
860 
865 
866 private:
867 
869  void setBeamA(long id);
871  long getBeamA() const;
873  void setBeamB(long id);
875  long getBeamB() const;
877  void setEBeamA(Energy e);
879  Energy getEBeamA() const;
881  void setEBeamB(Energy e);
883  Energy getEBeamB() const;
885  void setPDFA(PDFPtr);
887  PDFPtr getPDFA() const;
889  void setPDFB(PDFPtr);
891  PDFPtr getPDFB() const;
892 
893 private:
894 
899 
904 
905 public:
906 
910  class LesHouchesInconsistencyError: public Exception {};
911 
914  class LesHouchesReopenWarning: public Exception {};
915 
918  class LesHouchesReopenError: public Exception {};
919 
922  class LesHouchesInitError: public InitException {};
925 };
926 
928 ostream & operator<<(ostream & os, const HEPEUP & h);
929 
930 }
931 
932 
933 #include "ThePEG/Utilities/ClassTraits.h"
934 
935 namespace ThePEG {
936 
943 template <>
946  typedef HandlerBase NthBase;
947 };
948 
954 template <>
956  : public ClassTraitsBase<LesHouchesReader> {
960  static string className() { return "ThePEG::LesHouchesReader"; }
966  static string library() { return "LesHouches.so"; }
967 
968 };
969 
972 }
973 
974 #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...
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:119
map< int, double > maxWeights
Individual maximum weights for individual (possibly reweighted) processes.
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.
virtual bool preInitialize() const
Return true if this object needs to be initialized before all other objects because it needs to extra...
long maxScan() const
The maximum number of events to scan to collect information about processes and cross sections...
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...
static AbstractClassDescription< LesHouchesReader > initLesHouchesReader
Describe an abstract base class with persistent data.
const Cuts & cuts() const
The Cuts object to be used for this reader.
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.
long currentPosition() const
The number of events produced so far.
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.
tCascHdlPtr CKKWHandler() const
Return a possibly null pointer to a CascadeHandler to be used for CKKW-reweighting.
virtual void initialize(LesHouchesEventHandler &eh)
Initialize.
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.
string cacheFileName() const
Name of file used to cache the events form the reader in a fast-readable form.
int maxMultCKKW() const
If this reader is to be used (possibly together with others) for CKKW reweighting and veto...
PersistentOStream is used to write objects persistently to a stream from which they can be read in ag...
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
Energy getEBeamA() const
Access function for the interface.
virtual void createBeams()
Create instances of the incoming beams in the event and store them in particleIndex.
const XSecStat & xSecStats() const
The collected cross section statistics for this reader.
unsigned int theMomentumTreatment
Option for the treatment of the momenta supplied.
bool doInitPDFs
Should PDFBase objects be constructed from the information in the event file in the initialization...
const map< string, double > & optionalEventWeights() const
Return the optional named weights associated to the current event.
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:509
long NEvents() const
The number of events found in this reader.
bool active() const
Return true if this reader is active.
bool reweightPDF
Should the event be reweighted by PDFs used by the PartonExtractor?
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 ...
PBIPair thePartonBinInstances
The pair of PartonBinInstance objects describing the current incoming partons in the event...
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.
bool negativeWeights() const
True if negative weights may be produced.
PDFPtr getPDFB() const
Access function for the interface.
PVector theOutgoing
The instances of the outgoing particles from the sub process for the current event.
const PartonPairVec & partonBins() const
The pairs of PartonBin objects describing the partons which can be extracted by the PartonExtractor o...
const PBIPair & partonBinInstances() const
The pair of PartonBinInstance objects describing the current incoming partons in the 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:144
StatMap statmap
Collect statistics for each individual process.
const PVector & intermediates() const
Return the instances of the intermediate particles in the sub process for the current event...
virtual void dofinish()
Finalize this object.
Energy getEBeamB() const
Access function for the interface.
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.
long getBeamA() const
Access function for the interface.
QTY< 0, 1, 0 >::Type Energy
Energy.
Definition: Unitsystem.h:34
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.
const XCombMap & xCombs() const
The map of XComb objects indexed by the corresponding PartonBin pair.
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. ...
double eventWeight() const
The current event weight given by XWGTUP times possible reweighting.
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.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
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
LesHouchesReader & operator=(const LesHouchesReader &)
Private and non-existent assignment operator.
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.
PDFPtr getPDFA() const
Access function for the interface.
virtual void createParticles()
Create instances of all particles in the event and store them in particleIndex.
bool useWeightWarnings
Set to true if warnings about possible weight incompatibilities should be issued. ...
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.
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.
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
const StatMap & processStats() const
Collected statistics about the individual processes.
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.
long getBeamB() const
Access function for the interface.
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
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.
pair< PDFPtr, PDFPtr > inPDF
The PDFBase objects which has been used for the beam particle when generating the events being read...
int minMultCKKW() const
If this reader is to be used (possibly together with others) for CKKW reweighting and veto...
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.
This template class allows the compiler to check calculations with physical quantities for dimensiona...
Definition: PhysicalQty.h:81
void reopen()
Reopen a reader.
bool doCutEarly
Determines whether to apply cuts to events before converting them to ThePEG format.
const PVector & outgoing() const
Return the instances of the outgoing particles from the sub process for the current event...
pair< cPDFPtr, cPDFPtr > outPDF
The PDFBase object to be used in the subsequent generation.
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.
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
void NEvents(long x)
The number of events in this reader.
CFile cacheFile() const
File stream for the cache.
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.
double lastweight
The weight multiplying the last read event due to PDF reweighting, CKKW reweighting or assigned rewei...
tPExtrPtr partonExtractor() const
The PartonExtractor object used to construct remnants.
const PPair & incoming() const
Return the instances of the incoming particles to the sub process for the current event...
void cacheEvent() const
Write the current event to the cache file.
const PPair & beams() const
Return the instances of the beam particles for the current event.
static void Init()
Standard Init function used to initialize the interfaces.
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:52
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
HEPEUP hepeup
The HEPEUP common block.
void closeCacheFile()
Close the cache file;.