thepeg is hosted by Hepforge, IPPP Durham
ThePEG  2.1.5
XSecStat.h
1 // -*- C++ -*-
2 //
3 // XSecStat.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2017 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef THEPEG_XSecStat_H
10 #define THEPEG_XSecStat_H
11 //
12 // This is the declaration of the XSecStat class.
13 //
14 
15 #include "ThePEG/Config/ThePEG.h"
16 
17 namespace ThePEG {
18 
36 class XSecStat {
37 
38 public:
39 
43  enum {
44  plainWeights = 0,
45  plainVetoedWeights,
46  reweightedWeights,
47  reweightedVetoedWeights
48  };
49 
57  theSumWeights (),
58  theSumWeights2(), theLastWeight(0.0) {}
59 
64  explicit XSecStat(CrossSection xsecmax)
65  : theMaxXSec(xsecmax), theAttempts(0), theAccepted(0), theVetoed(0),
66  theSumWeights (),
67  theSumWeights2(), theLastWeight(0.0) {}
68 
72  XSecStat & operator=(const XSecStat & x) = default;
73 
77  XSecStat & operator+=(const XSecStat & x) {
80  theVetoed += x.theVetoed;
81  for( unsigned int ix = 0; ix < 4; ++ix ) {
82  theSumWeights [ix] += x.theSumWeights [ix];
83  theSumWeights2[ix] += x.theSumWeights2[ix];
84  }
85  theLastWeight = 0.0;
86  return *this;
87  }
88 
92  void reset() {
95  theLastWeight = 0.0;
96  }
97 
99 
100 public:
101 
104 
109  void accept() {
110  theAccepted += 1;
111  }
112 
117  void select(double weight) {
118  theAttempts += 1;
119  theSumWeights [reweightedWeights] += weight ;
120  theSumWeights2[reweightedWeights] += sqr(weight);
121  theSumWeights [plainWeights] += weight ;
122  theSumWeights2[plainWeights] += sqr(weight);
123  theLastWeight = weight;
124  }
125 
129  void reweight(double oldWeight, double newWeight) {
130  theSumWeights [reweightedWeights] += newWeight - oldWeight ;
131  theSumWeights2[reweightedWeights] += sqr(newWeight) - sqr(oldWeight);
132  }
133 
142  void reject(double weight = 1.0) {
143  theSumWeights [reweightedVetoedWeights] += weight ;
144  theSumWeights2[reweightedVetoedWeights] += sqr(weight);
145  theSumWeights [plainVetoedWeights] += theLastWeight ;
146  theSumWeights2[plainVetoedWeights] += sqr(theLastWeight);
147  theVetoed += 1;
148  }
149 
153  CrossSection maxXSec() const { return theMaxXSec; }
154 
158  double sumWeights() const {
159  return theSumWeights[reweightedWeights] - theSumWeights[reweightedVetoedWeights];
160  }
161 
165  double sumWeights2() const {
166  return theSumWeights2[reweightedWeights] + theSumWeights2[reweightedVetoedWeights];
167  }
168 
172  double sumWeightsNoReweight() const {
173  return theSumWeights[plainWeights] - theSumWeights[plainVetoedWeights];
174  }
175 
179  double sumWeights2NoReweight() const {
180  return theSumWeights2[plainWeights] + theSumWeights2[plainVetoedWeights];
181  }
182 
188  CrossSection xSec(double att = 0) const {
189  double n = (att == 0.0 ? attempts() : att);
190  return n ? maxXSec()*sumWeights()/n : maxXSec();
191  }
192 
198  CrossSection xSecErr(double att = 0) const {
199  double n = (att == 0.0 ? attempts() : att);
200  if ( n < 2 )
201  return maxXSec();
202  double sw = sumWeights(); double sw2 = sumWeights2();
203  return
204  maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
205  }
206 
212  CrossSection xSecNoReweight(double att = 0) const {
213  double n = (att == 0.0 ? attempts() : att);
214  return n ? maxXSec()*sumWeightsNoReweight()/n : maxXSec();
215  }
216 
222  CrossSection xSecErrNoReweight(double att = 0) const {
223  double n = (att == 0.0 ? attempts() : att);
224  if ( n < 2 )
225  return maxXSec();
226  double sw = sumWeightsNoReweight();
227  double sw2 = sumWeights2NoReweight();
228  return
229  maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
230  }
231 
235  double attempts() const { return theAttempts; }
236 
240  double accepted() const { return theAccepted-theVetoed; }
241 
245  double vetoed() const { return theVetoed; }
246 
250  void maxXSec(CrossSection x) { theMaxXSec = x; }
252 
253 public:
254 
260  void output(PersistentOStream & os) const;
261 
265  void input(PersistentIStream & is);
267 
268 private:
269 
274 
278  double theAttempts;
279 
283  double theAccepted;
284 
288  double theVetoed;
289 
293  array<double,4> theSumWeights;
294 
298  array<double,4> theSumWeights2;
299 
304 
305 };
306 
309 
312 
314 inline XSecStat operator+(const XSecStat & x1, const XSecStat & x2) {
315  XSecStat x = x1;
316  return x += x2;
317 }
318 
319 }
320 
321 #endif /* THEPEG_XSecStat_H */
double sumWeights2NoReweight() const
The sum of the squared weights so far, excluding reweighting.
Definition: XSecStat.h:179
double attempts() const
Number of attempts so far.
Definition: XSecStat.h:235
void reset()
Reset the statistics.
Definition: XSecStat.h:92
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:117
double sumWeightsNoReweight() const
The sum of the weights so far, excluding reweighting.
Definition: XSecStat.h:172
void accept()
An event of the corresponding class has been accepted.
Definition: XSecStat.h:109
PersistentOStream is used to write objects persistently to a stream from which they can be read in ag...
XSecStat()
The default constructor.
Definition: XSecStat.h:55
void input(PersistentIStream &is)
Input from a persistent stream.
double theVetoed
Number of events vetoed after being accepted.
Definition: XSecStat.h:288
This is the main namespace within which all identifiers in ThePEG are declared.
Definition: FactoryBase.h:28
This is the main config header file for ThePEG.
double vetoed() const
Number of vetoes so far.
Definition: XSecStat.h:245
CrossSection theMaxXSec
The overestimated cross section.
Definition: XSecStat.h:273
void reject(double weight=1.0)
Reject the event which was last accepted with accept() or selected with select(double).
Definition: XSecStat.h:142
double sumWeights() const
The sum of the weights so far.
Definition: XSecStat.h:158
double accepted() const
Number of accepts so far.
Definition: XSecStat.h:240
double sumWeights2() const
The sum of the squared weights so far.
Definition: XSecStat.h:165
XSecStat & operator=(const XSecStat &x)=default
The assignment operator.
double theLastWeight
The last selected weight, ignoring reweighting.
Definition: XSecStat.h:303
void reweight(double oldWeight, double newWeight)
Reweight a selected and accepted event.
Definition: XSecStat.h:129
void maxXSec(CrossSection x)
Set the overestimated cross section.
Definition: XSecStat.h:250
CrossSection xSec(double att=0) const
The current estimate of the cross section for the corresponding class of events.
Definition: XSecStat.h:188
CrossSection maxXSec() const
The overestimated cross section.
Definition: XSecStat.h:153
CrossSection xSecErr(double att=0) const
The current estimate of the error in the cross section for the corresponding class of events...
Definition: XSecStat.h:198
vector< T > & operator>>(vector< T > &tv, U &u)
Overload the right shift operator for vector to pop objects from a vector.
Definition: Containers.h:192
double theAttempts
Number of attempts so far.
Definition: XSecStat.h:278
array< double, 4 > theSumWeights
The sum of the weights so far.
Definition: XSecStat.h:293
CrossSection xSecNoReweight(double att=0) const
The current estimate of the cross section for the corresponding class of events, excluding reweightin...
Definition: XSecStat.h:212
void output(PersistentOStream &os) const
Output to a persistent stream.
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
This template class allows the compiler to check calculations with physical quantities for dimensiona...
Definition: PhysicalQty.h:81
array< double, 4 > theSumWeights2
The sum of the squared weights so far.
Definition: XSecStat.h:298
double theAccepted
Number of accepted events so far.
Definition: XSecStat.h:283
CrossSection xSecErrNoReweight(double att=0) const
The current estimate of the error in the cross section for the corresponding class of events...
Definition: XSecStat.h:222
XSecStat & operator+=(const XSecStat &x)
Add the contents of another XSecStat.
Definition: XSecStat.h:77
XSecStat(CrossSection xsecmax)
Constructor taking the overestimated cross section, xsecmax, as argument.
Definition: XSecStat.h:64
XSecStat is a concrete helper class used to collect statistics about the cross section for a specific...
Definition: XSecStat.h:36
constexpr ZeroUnit ZERO
ZERO can be used as zero for any unitful quantity.
Definition: PhysicalQty.h:33