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 
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 
648 
654 
659 
665 
670  long position;
671 
675  int reopened;
676 
683 
687  bool scanning;
688 
692  bool isActive;
693 
699 
705 
710 
715 
721 
727 
733 
738 
744 
750 
756 
761 
766 
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 <>
944 struct BaseClassTrait<LesHouchesReader,1>: public ClassTraitsType {
946  typedef HandlerBase NthBase;
947 };
948 
954 template <>
955 struct ClassTraits<LesHouchesReader>
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 */