ThePEG/Utilities/XSecStat.h

0001 // -*- C++ -*-
0002 //
0003 // XSecStat.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_XSecStat_H
0010 #define THEPEG_XSecStat_H
0011 //
0012 // This is the declaration of the XSecStat class.
0013 //
0014
0015 #include "ThePEG/Config/ThePEG.h"
0016
0017 namespace ThePEG {
0018
0019 /**
0020  * XSecStat is a concrete helper class used to collect statistics
0021  * about the cross section for a specific process or group of
0022  * processes. It contains an overestimated cross section and
0023  * information about the number of times the process has been used to
0024  * generate an event and how many times this event has been accepted.
0025  *
0026  * An object of this class must initially be given an overestimated
0027  * cross section in the constructor or with the maxXSec(CrossSection)
0028  * function. Each time the corresponding process is selected
0029  * (according to maxXSec()), the select(double) function should be
0030  * called giving the weight with which the event will be accepted as
0031  * argument. If the event is then accepted, the accept() function
0032  * should be called. If an event is later vetoed, the reject()
0033  * function should be called.
0034  *
0035  */
0036 class XSecStat {
0037
0038 public:
0039
0040   /**
0041    * Enumerate the different weight classes
0042    */
0043   enum {
0044     plainWeights = 0,
0045     plainVetoedWeights,
0046     reweightedWeights,
0047     reweightedVetoedWeights
0048   };
0049
0050   /** @name Standard constructors, destructor and assignment operator. */
0051   //@{
0052   /**
0053    * The default constructor.
0054    */
0055   XSecStat()
0056     : theMaxXSec(ZERO), theAttempts(0), theAccepted(0), theVetoed(0),
0057       theSumWeights (),
0058       theSumWeights2(), theLastWeight(0.0) {}
0059
0060   /**
0061    * Constructor taking the overestimated cross section, \a xsecmax,
0062    * as argument.
0063    */
0064   explicit XSecStat(CrossSection xsecmax)
0065     : theMaxXSec(xsecmax), theAttempts(0), theAccepted(0), theVetoed(0),
0066       theSumWeights (),
0067       theSumWeights2(), theLastWeight(0.0) {}
0068
0069   /**
0070    * The assignment operator.
0071    */
0072   XSecStat & operator=(const XSecStat & x) = default;
0073
0074   /**
0075    * Add the contents of another XSecStat.
0076    */
0077   XSecStat & operator+=(const XSecStat & x) {
0078     theAttempts    += x.theAttempts;
0079     theAccepted    += x.theAccepted;
0080     theVetoed      += x.theVetoed;
0081     for( unsigned int ix = 0; ix < 4; ++ix ) {
0082       theSumWeights [ix] += x.theSumWeights [ix];
0083       theSumWeights2[ix] += x.theSumWeights2[ix];
0084     }
0085     theLastWeight = 0.0;
0086     return *this;
0087   }
0088
0089   /**
0090    * Reset the statistics.
0091    */
0092   void reset() {
0093     theAttempts = theAccepted = theVetoed = 0;
0094     theSumWeights = theSumWeights2 = {};
0095     theLastWeight = 0.0;
0096   }
0097
0098   //@}
0099
0100 public:
0101
0102   /** @name Simple access functions */
0103   //@{
0104
0105   /**
0106    * An event of the corresponding class has been accepted. The
0107    * select() method must have been called before.
0108    */
0109   void accept() {
0110     theAccepted += 1;
0111   }
0112
0113   /**
0114    * An event of the corresponding class has been attempted. It will
0115    * subsequently be accepted with the given \a weight.
0116    */
0117   void select(double weight) {
0118     theAttempts += 1;
0119     theSumWeights [reweightedWeights] +=     weight ;
0120     theSumWeights2[reweightedWeights] += sqr(weight);
0121     theSumWeights [plainWeights]      +=     weight ;
0122     theSumWeights2[plainWeights]      += sqr(weight);
0123     theLastWeight = weight;
0124   }
0125
0126   /**
0127    * Reweight a selected and accepted event.
0128    */
0129   void reweight(double oldWeight, double newWeight) {
0130     theSumWeights [reweightedWeights] +=     newWeight  -     oldWeight ;
0131     theSumWeights2[reweightedWeights] += sqr(newWeight) - sqr(oldWeight);
0132   }
0133
0134   /**
0135    * Reject the event which was last accepted with accept() or
0136    * selected with select(double). The \a weight should be set to the
0137    * value, \f$w\f$, used in the previous call to select(double),
0138    * except if the event has been accepted with the probability
0139    * \f$w\f$, in which case \a weight should be set to \f$sign(1,
0140    * w)\f$.
0141    */
0142   void reject(double weight = 1.0) {
0143     theSumWeights [reweightedVetoedWeights] +=            weight ;
0144     theSumWeights2[reweightedVetoedWeights] +=        sqr(weight);
0145     theSumWeights [plainVetoedWeights]      +=     theLastWeight ;
0146     theSumWeights2[plainVetoedWeights]      += sqr(theLastWeight);
0147     theVetoed += 1;
0148   }
0149
0150   /**
0151    * The overestimated cross section.
0152    */
0153   CrossSection maxXSec() const { return theMaxXSec; }
0154
0155   /**
0156    * The sum of the weights so far.
0157    */
0158   double sumWeights() const {
0159     return theSumWeights[reweightedWeights] - theSumWeights[reweightedVetoedWeights];
0160   }
0161
0162   /**
0163    * The sum of the squared weights so far.
0164    */
0165   double sumWeights2() const {
0166     return theSumWeights2[reweightedWeights] + theSumWeights2[reweightedVetoedWeights];
0167   }
0168
0169   /**
0170    * The sum of the weights so far, excluding reweighting.
0171    */
0172   double sumWeightsNoReweight() const {
0173     return theSumWeights[plainWeights] - theSumWeights[plainVetoedWeights];
0174   }
0175
0176   /**
0177    * The sum of the squared weights so far, excluding reweighting.
0178    */
0179   double sumWeights2NoReweight() const {
0180     return theSumWeights2[plainWeights] + theSumWeights2[plainVetoedWeights];
0181   }
0182
0183   /**
0184    * The current estimate of the cross section for the corresponding
0185    * class of events. If no events have been generated, maxXSec() will
0186    * be returned.
0187    */
0188   CrossSection xSec(double att = 0) const {
0189     double n = (att == 0.0 ? attempts() : att);
0190     return n ? maxXSec()*sumWeights()/n : maxXSec();
0191   }
0192
0193   /**
0194    * The current estimate of the error in the cross section for the
0195    * corresponding class of events. If no events have been generated,
0196    * maxXSec() will be returned.
0197    */
0198   CrossSection xSecErr(double att = 0) const {
0199     double n = (att == 0.0 ? attempts() : att);
0200     if ( n < 2 )
0201       return maxXSec();
0202     double sw = sumWeights(); double sw2 = sumWeights2();
0203     return
0204       maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
0205   }
0206
0207   /**
0208    * The current estimate of the cross section for the corresponding
0209    * class of events, excluding reweighting. If no events have been
0210    * generated, maxXSec() will be returned.
0211    */
0212   CrossSection xSecNoReweight(double att = 0) const {
0213     double n = (att == 0.0 ? attempts() : att);
0214     return n ? maxXSec()*sumWeightsNoReweight()/n : maxXSec();
0215   }
0216
0217   /**
0218    * The current estimate of the error in the cross section for the
0219    * corresponding class of events, excluding reweighting. If no
0220    * events have been generated, maxXSec() will be returned.
0221    */
0222   CrossSection xSecErrNoReweight(double att = 0) const {
0223     double n = (att == 0.0 ? attempts() : att);
0224     if ( n < 2 )
0225       return maxXSec();
0226     double sw = sumWeightsNoReweight();
0227     double sw2 = sumWeights2NoReweight();
0228     return
0229       maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
0230   }
0231
0232   /**
0233    * Number of attempts so far.
0234    */
0235   double attempts() const { return theAttempts; }
0236
0237   /**
0238    * Number of accepts so far.
0239    */
0240   double accepted() const { return theAccepted-theVetoed; }
0241
0242   /**
0243    * Number of vetoes so far.
0244    */
0245   double vetoed() const { return theVetoed; }
0246
0247   /**
0248    * Set the overestimated cross section.
0249    */
0250   void maxXSec(CrossSection x) { theMaxXSec = x; }
0251   //@}
0252
0253 public:
0254
0255   /** @name I/O functions */
0256   //@{
0257   /**
0258    * Output to a persistent stream.
0259    */
0260   void output(PersistentOStream & os) const;
0261
0262   /**
0263    * Input from a persistent stream.
0264    */
0265   void input(PersistentIStream & is);
0266   //@}
0267
0268 private:
0269
0270   /**
0271    * The overestimated cross section.
0272    */
0273   CrossSection theMaxXSec;
0274
0275   /**
0276    * Number of attempts so far.
0277    */
0278   double theAttempts;
0279
0280   /**
0281    * Number of accepted events so far.
0282    */
0283   double theAccepted;
0284
0285   /**
0286    * Number of events vetoed after being accepted
0287    */
0288   double theVetoed;
0289
0290   /**
0291    * The sum of the weights so far.
0292    */
0293   array<double,4> theSumWeights;
0294
0295   /**
0296    * The sum of the squared weights so far.
0297    */
0298   array<double,4> theSumWeights2;
0299
0300   /**
0301    * The last selected weight, ignoring reweighting.
0302    */
0303   double theLastWeight;
0304
0305 };
0306
0307 /** Ouptut an XSecStat to a persistent stream. */
0308 PersistentOStream & operator<<(PersistentOStream &, const XSecStat &);
0309
0310 /** Input an XSecStat from a persistent stream. */
0311 PersistentIStream & operator>>(PersistentIStream &, XSecStat &);
0312
0313 /** Add the contents of two XSecStat objects. */
0314 inline XSecStat operator+(const XSecStat & x1, const XSecStat & x2) {
0315   XSecStat x = x1;
0316   return x += x2;
0317 }
0318
0319 }
0320
0321 #endif /* THEPEG_XSecStat_H */