Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-08-06 09:38:22

0001 // -*- C++ -*-
0002 //
0003 // SamplerBase.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_SamplerBase_H
0010 #define ThePEG_SamplerBase_H
0011 // This is the declaration of the SamplerBase class.
0012 
0013 #include "ThePEG/Interface/Interfaced.h"
0014 #include "SamplerBase.fh"
0015 #include "ThePEG/Handlers/StandardEventHandler.fh"
0016 // #include "SamplerBase.xh"
0017 
0018 namespace ThePEG {
0019 
0020 /**
0021  * This is the base class for all phase space sampler classes to be
0022  * used by the EventHandler class to sample the phase space according
0023  * to the cross sections for the processes in the EventHandler. The
0024  * class should be able to sample a unit hyper-cube in arbitrary
0025  * dimensions. The points need not necessarily be sampled with unit
0026  * weight.
0027  *
0028  * The virtual methods to be implemented by concrete sub-classes are
0029  * initialize(), generate() and rejectLast().
0030  *
0031  * @see \ref SamplerBaseInterfaces "The interfaces"
0032  * defined for SamplerBase.
0033  * @see EventHandler
0034  */
0035 class SamplerBase: public Interfaced {
0036 
0037 public:
0038 
0039   /** @name Standard constructors and destructors. */
0040   //@{
0041 
0042   /**
0043    * Constructor
0044    */
0045   SamplerBase()
0046     : Interfaced(), 
0047       theIntegrationList("") {}
0048 
0049   /**
0050    * Destructor.
0051    */
0052   virtual ~SamplerBase();
0053   //@}
0054 
0055 public:
0056 
0057   /**
0058    * Set the event handler for which the function
0059    * StandardEventHandler::dSigDR(const vector<double> &) function
0060    * returns the cross section for the chosen phase space point.
0061    */
0062   void setEventHandler(tStdEHPtr eh) { theEventHandler = eh; }
0063 
0064   /** @name Virtual functions to be overridden by sub-classes. */
0065   //@{
0066   /**
0067    * Initialize the the sampler, possibly doing presampling of the
0068    * phase space.
0069    */
0070   virtual void initialize() = 0;
0071 
0072   /**
0073    * An external hook to prepare the sampler for generating events, e.g. by
0074    * combining grid files from parallel integration runs.
0075    */
0076   virtual void prepare() {}
0077 
0078   /**
0079    * Generarate a new phase space point and return a weight associated
0080    * with it. This weight should preferably be 1.
0081    */
0082   virtual double generate() = 0;
0083 
0084   /**
0085    * Reject the last chosen phase space point.
0086    */
0087   virtual void rejectLast() = 0;
0088 
0089   /**
0090    * Return the last generated phase space point.
0091    */
0092   const vector<double> & lastPoint() const { return theLastPoint; }
0093 
0094   /**
0095    * If the sampler is able to sample several different functions
0096    * separately, this function should return the last chosen
0097    * function. This default version always returns 0.
0098    */
0099   virtual int lastBin() const { return 0; }
0100 
0101   /**
0102    * Return the total integrated cross section determined from the
0103    * Monte Carlo sampling so far.
0104    */
0105   virtual CrossSection integratedXSec() const = 0;
0106 
0107   /**
0108    * Return the error on the total integrated cross section determined
0109    * from the Monte Carlo sampling so far.
0110    */
0111   virtual CrossSection integratedXSecErr() const = 0;
0112 
0113   /**
0114    * Return the reference cross section, a.k.a. maximum weight. When
0115    * not provided directly, this will be determined effectively from
0116    * the sum of weights and sum of weights squared to match up the
0117    * standard definition of a Monte Carlo cross section along with the
0118    * cross section and error quoted.
0119    */
0120   virtual CrossSection maxXSec() const {
0121     if ( sumWeights2() <= 0.0 ) return ZERO;
0122     return integratedXSec()*attempts()/sumWeights();
0123   }
0124 
0125   /**
0126    * Return the number of attempts. When not provided directly, this
0127    * will be determined effectively from the sum of weights and sum of
0128    * weights squared to match up the standard definition of a Monte
0129    * Carlo cross section along with the cross section and error
0130    * quoted.
0131    */
0132   virtual double attempts() const {
0133     CrossSection sigma = integratedXSec();
0134     CrossSection esigma = integratedXSecErr();
0135     double sw = sumWeights(); double sw2 = sumWeights2();
0136     if ( sw2 <= 0.0 ) return 0.0;
0137     return 
0138       sqr(sw)*(sqr(esigma)-sqr(sigma))/(sqr(sw)*sqr(esigma) - sw2*sqr(sigma));
0139   }
0140 
0141   /**
0142    * Return the sum of the weights returned by generate() so far (of
0143    * the events that were not rejeted).
0144    */
0145   virtual double sumWeights() const = 0;
0146 
0147   /**
0148    * Return the sum of the weights squared returned by generate() so
0149    * far (of the events that were not rejeted).
0150    */
0151   virtual double sumWeights2() const = 0;
0152 
0153   /** 
0154    * Return true if this sampler is generating almost unweighted events.
0155    */ 
0156   virtual bool almostUnweighted() const { return false; }
0157   //@}
0158 
0159   /** @name Controlling of run levels and grid handling*/
0160   //@{
0161 
0162   /**
0163    * Set a file containing a list of subprocesses to integrate
0164    */
0165   void integrationList(const string& newIntegrationList) { theIntegrationList = newIntegrationList; }
0166 
0167   /**
0168    * Return a file containing a list of subprocesses to integrate
0169    */
0170   const string& integrationList() const { return theIntegrationList; }
0171 
0172   /**
0173    * Enumerate the possible run levels
0174    */
0175   enum RunLevels {
0176 
0177     UnknownMode = 0,
0178     InitMode,
0179     ReadMode,
0180     BuildMode,
0181     IntegrationMode,
0182     RunMode
0183 
0184   };
0185 
0186   /**
0187    * Return the run level
0188    */
0189   static int runLevel() {
0190     return theRunLevel();
0191   }
0192 
0193   /**
0194    * Set the run level
0195    */
0196   static void setRunLevel(int level) {
0197     theRunLevel() = level;
0198   }
0199 
0200   /**
0201    * Return true, if a setupfile is in use
0202    */
0203   static bool hasSetupFile() {
0204     return theHasSetupFile();
0205   }
0206 
0207   /**
0208    * Indicate that a setupfile is in use.
0209    */
0210   static void setupFileUsed(bool yes = true) {
0211     theHasSetupFile() = yes;
0212   }
0213 
0214   /**
0215    * Return the seed that has been used for this run to disentangle
0216    * grids whihch have been adapted further
0217    */
0218   static long seed() {
0219     return theSeed();
0220   }
0221 
0222   /**
0223    * Set the seed that has been used for this run to disentangle
0224    * grids whihch have been adapted further
0225    */
0226   static void setSeed(long s) {
0227     theSeed() = s;
0228   }
0229 
0230   /**
0231    * Return the number of subprocesses to be integrated per job.
0232    */
0233   static unsigned int integratePerJob() {
0234     return theIntegratePerJob();
0235   }
0236 
0237   /**
0238    * Set the number of subprocesses to be integrated per job.
0239    */
0240   static void setIntegratePerJob(unsigned int s) {
0241     theIntegratePerJob() = s;
0242   }
0243 
0244   /**
0245    * Return the maximum number of integration jobs to be created.
0246    */
0247   static unsigned int integrationJobs() {
0248     return theIntegrationJobs();
0249   }
0250 
0251   /**
0252    * Set the maximum number of integration jobs to be created.
0253    */
0254   static void setIntegrationJobs(unsigned int s) {
0255     theIntegrationJobs() = s;
0256   }
0257 
0258   //@}
0259 
0260 protected:
0261 
0262   /**
0263    * Return the last generated phase space point.
0264    */
0265   vector<double> & lastPoint() { return theLastPoint; }
0266 
0267   /**
0268    * Return the associated event handler.
0269    */
0270   tStdEHPtr eventHandler() const { return theEventHandler; }
0271 
0272 public:
0273 
0274   /** @name Functions used by the persistent I/O system. */
0275   //@{
0276   /**
0277    * Function used to write out object persistently.
0278    * @param os the persistent output stream written to.
0279    */
0280   void persistentOutput(PersistentOStream & os) const;
0281 
0282   /**
0283    * Function used to read in object persistently.
0284    * @param is the persistent input stream read from.
0285    * @param version the version number of the object when written.
0286    */
0287   void persistentInput(PersistentIStream & is, int version);
0288   //@}
0289 
0290   /**
0291    * Standard Init function used to initialize the interfaces.
0292    */
0293   static void Init();
0294 
0295 private:
0296 
0297   /**
0298    * The associated event handler.
0299    */
0300   tStdEHPtr theEventHandler;
0301 
0302   /**
0303    * The last generated phase space point.
0304    */
0305   vector<double> theLastPoint;
0306 
0307   /**
0308    * A file containing a list of subprocesses to integrate
0309    */
0310   string theIntegrationList;
0311 
0312   /**
0313    * The run level
0314    */
0315   static int& theRunLevel() {
0316     static int lvl = UnknownMode;
0317     return lvl;
0318   }
0319 
0320   /**
0321    * True, if a setupfile is in use
0322    */
0323   static bool& theHasSetupFile() {
0324     static bool flag = false;
0325     return flag;
0326   }
0327 
0328   /**
0329    * The seed that has been used for this run to disentangle
0330    * grids whihch have been adapted further
0331    */
0332   static long& theSeed() {
0333     static long s = 0;
0334     return s;
0335   }
0336 
0337   /**
0338    * The number of subprocesses to be integrated per job.
0339    */
0340   static unsigned int& theIntegratePerJob() {
0341     static unsigned int s = 0;
0342     return s;
0343   }
0344 
0345   /**
0346    * The maximum number of integration jobs to be created.
0347    */
0348   static unsigned int& theIntegrationJobs() {
0349     static unsigned int s = 0;
0350     return s;
0351   }
0352 
0353 private:
0354 
0355   /**
0356    * Describe an abstract base class with persistent data.
0357    */
0358   static AbstractClassDescription<SamplerBase> initSamplerBase;
0359 
0360   /**
0361    *  Private and non-existent assignment operator.
0362    */
0363   SamplerBase & operator=(const SamplerBase &) = delete;
0364 
0365 };
0366 
0367 }
0368 
0369 
0370 namespace ThePEG {
0371 
0372 /** @cond TRAITSPECIALIZATIONS */
0373 
0374 /**
0375  * This template specialization informs ThePEG about the base class of
0376  * SamplerBase.
0377  */
0378 template <>
0379 struct BaseClassTrait<SamplerBase,1>: public ClassTraitsType {
0380   /** Typedef of the base class of SamplerBase. */
0381   typedef Interfaced NthBase;
0382 };
0383 
0384 /**
0385  * This template specialization informs ThePEG about the name of the
0386  * SamplerBase class.
0387  */
0388 template <>
0389 struct ClassTraits<SamplerBase>: public ClassTraitsBase<SamplerBase> {
0390   /** Return the class name. */
0391   static string className() { return "ThePEG::SamplerBase"; }
0392 
0393 };
0394 
0395 /** @endcond */
0396 
0397 }
0398 
0399 #endif /* ThePEG_SamplerBase_H */