Sampling/exsample/generator.h

0001 // -*- C++ -*-
0002 //
0003 // generator.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
0004 //
0005 // Copyright (C) 2008-2019 Simon Platzer -- simon.plaetzer@desy.de, The Herwig Collaboration
0006 //
0007 // ExSample is licenced under version 3 of the GPL, see COPYING for details.
0008 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0009 //
0010 //
0011 #ifndef EXSAMPLE_generator_h_included
0012 #define EXSAMPLE_generator_h_included
0013
0014 #include "cell.h"
0015 #include "selectors.h"
0016 #include "statistics.h"
0017 #include "binary_tree.h"
0018
0019 namespace exsample {
0020
0021   /// \brief Exception thrown, if the generator has just changed its
0022   /// state. The attempt of generating an event should be repeated.
0023   struct generator_update{};
0024
0025   /// \brief A generator for plain sampling and integrating
0026   template<class Function, class Random>
0027   class generator {
0028
0029   public:
0030
0031     /// default constructor
0032     generator()
0033       : function_(0), statistics_(), check_events_(0),
0034     adaption_info_(), root_cell_(),
0035     rnd_gen_(), did_split_(false), initialized_(false),
0036     last_cell_(), last_point_(), last_value_(0.),
0037     compensating_(false) {}
0038
0039   public:
0040
0041     /// initialize this generator
0042     template<class SlaveStatistics>
0043     void initialize(SlaveStatistics&);
0044
0045     /// generate an event, returning
0046     /// the sign of the weight
0047     template<class SlaveStatistics>
0048     double generate(SlaveStatistics&);
0049
0050     /// return the last sampled phase space point
0051     const std::vector<double>& last_point() const { return last_point_; }
0052
0053     /// indicate that the last generated point has been rejected
0054     void reject() {
0055       statistics_.reject(last_value_);
0056       last_cell_->info().reject();
0057     }
0058
0059     /// finalize this generator
0060     void finalize() {
0061       statistics_.reset();
0062     }
0063
0064   public:
0065
0066     /// return true, if this generator has been initialized
0067     bool initialized() const { return initialized_; }
0068
0069     /// return true, if at least one split has been performed
0070     bool did_split() const { return did_split_; }
0071
0072     /// access the function
0073     Function& function() { return *function_; }
0074
0075     /// set the function
0076     void function(Function * f) { function_ = f; }
0077
0078     /// return the statistics object
0079     const statistics& stats() const { return statistics_; }
0080
0081     /// return the sampled volume
0082     double volume() const {
0083       return exsample::volume(adaption_info_.lower_left, adaption_info_.upper_right);
0084     }
0085
0086     /// return the integral
0087     double integral() const {
0088       return volume() * statistics_.average_weight();
0089     }
0090
0091     /// return the error on the integral
0092     double integral_uncertainty() const {
0093       return volume() * std::sqrt(statistics_.average_weight_variance());
0094     }
0095
0096     /// return the integral
0097     double current_integral() const {
0098       return volume() * statistics_.current().first;
0099     }
0100
0101     /// return the error on the integral
0102     double current_integral_uncertainty() const {
0103       return volume() * std::sqrt(statistics_.current().second);
0104     }
0105
0106     /// return the variance of the integral estimate
0107     double integral_variance() const {
0108       return sqr(volume()) * statistics_.average_weight_variance();
0109     }
0110
0111     /// access the adaption_info object
0112     adaption_info& sampling_parameters() { return adaption_info_; }
0113
0114     /// return true, if still compensating
0115     bool compensating() const { return compensating_; }
0116
0117   public:
0118
0119     /// put to ostream
0120     template<class OStream>
0121     void put(OStream& os) const;
0122
0123     /// get from istream
0124     template<class IStream>
0125     void get(IStream& is);
0126
0127   private:
0128
0129     /// check for and possibly split
0130     /// the last selected cell
0131     bool split();
0132
0133     /// compensate the last selected cell indicating the
0134     /// value and position of the new overestimate
0135     void compensate();
0136
0137     /// function to be sampled
0138     Function * function_;
0139
0140     /// the global statistics object
0141     statistics statistics_;
0142
0143     /// the number of events after which
0144     /// a cell is checked for splits
0145     unsigned long check_events_;
0146
0147     /// the adaption info object
0148     adaption_info adaption_info_;
0149
0150     /// the root cell
0151     binary_tree<cell> root_cell_;
0152
0153     /// the random number generator to be used
0154     rnd_generator<Random> rnd_gen_;
0155
0156     /// wether a split has already been performed
0157     bool did_split_;
0158
0159     /// wether this generator has been initialized
0160     bool initialized_;
0161
0162     /// the last selected cell
0163     binary_tree<cell>::iterator last_cell_;
0164
0165     /// the last sampled phasespace point
0166     std::vector<double> last_point_;
0167
0168     /// the last function value
0169     double last_value_;
0170
0171     /// wether or not we are compensating
0172     bool compensating_;
0173
0174   };
0175
0176 }
0177
0178 #include "generator.icc"
0179
0180 #endif // EXSAMPLE_generator_h_included
0181