ThePEG  1.8.0
XSecStat.h
1 // -*- C++ -*-
2 //
3 // XSecStat.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_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 
47  theSumWeights (vector<double>(5,0.0)),
48  theSumWeights2(vector<double>(5,0.0)), theLastWeight(0.0) {}
49 
55  : theMaxXSec(xsecmax), theAttempts(0), theAccepted(0), theVetoed(0),
56  theSumWeights (vector<double>(5,0.0)),
57  theSumWeights2(vector<double>(5,0.0)), theLastWeight(0.0) {}
58 
62  XSecStat & operator=(const XSecStat & x) {
66  theVetoed = x.theVetoed;
70  return *this;
71  }
72 
76  XSecStat & operator+=(const XSecStat & x) {
80  theVetoed += x.theVetoed;
81  for(unsigned int ix=0;ix<5;++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() {
94  theSumWeights = theSumWeights2 = vector<double>(5,0.0);
95  theLastWeight = 0.0;
96  }
97 
99 
100 public:
101 
104 
109  void accept() {
110  ++theAccepted;
111  theSumWeights [1] += 1.;
112  theSumWeights2[1] += 1.;
113  }
114 
119  void select(double weight) {
120  ++theAttempts;
121  theSumWeights [0] += weight ;
122  theSumWeights2[0] += sqr(weight);
123  theSumWeights [3] += weight ;
124  theSumWeights2[3] += sqr(weight);
125  theLastWeight = weight;
126  }
127 
131  void reweight(double oldWeight, double newWeight) {
132  theSumWeights [0] += newWeight - oldWeight ;
133  theSumWeights2[0] += sqr(newWeight) - sqr(oldWeight);
134  }
135 
144  void reject(double weight = 1.0) {
145  theSumWeights [1] -= weight ;
146  theSumWeights2[1] -= sqr(weight);
147  theSumWeights [2] += weight ;
148  theSumWeights2[2] += sqr(weight);
150  theSumWeights2[4] += sqr(theLastWeight);
151  ++theVetoed;
152  }
153 
159  CrossSection xSec() const {
160  return attempts() ? maxXSec()*(theSumWeights[0]-theSumWeights[2])/attempts() : maxXSec();
161  }
162 
169  return attempts() ? maxXSec()*sqrt(theSumWeights2[0]+
170  theSumWeights2[2])/attempts() : maxXSec();
171  }
172 
179  return attempts() ? maxXSec()*(theSumWeights[3]-theSumWeights[4])/attempts() : maxXSec();
180  }
181 
188  return attempts() ? maxXSec()*sqrt(theSumWeights2[3]+
189  theSumWeights2[4])/attempts() : maxXSec();
190  }
191 
195  CrossSection maxXSec() const { return theMaxXSec; }
196 
200  double sumWeights() const { return theSumWeights[0] - theSumWeights[2]; }
201 
205  double sumWeights2() const { return theSumWeights2[0] + theSumWeights2[2]; }
206 
210  double sumWeightsNoReweight() const { return theSumWeights[3] - theSumWeights[4]; }
211 
215  double sumWeights2NoReweight() const { return theSumWeights2[3] + theSumWeights2[4]; }
216 
220  long attempts() const { return theAttempts; }
221 
225  long accepted() const { return theAccepted-theVetoed; }
226 
230  void maxXSec(CrossSection x) { theMaxXSec = x; }
232 
233 public:
234 
240  void output(PersistentOStream & os) const;
241 
245  void input(PersistentIStream & is);
247 
248 private:
249 
254 
259 
264 
268  long theVetoed;
269 
273  vector<double> theSumWeights;
274 
278  vector<double> theSumWeights2;
279 
284 
285 };
286 
289 
292 
294 inline XSecStat operator+(const XSecStat & x1, const XSecStat & x2) {
295  XSecStat x = x1;
296  return x += x2;
297 }
298 
299 }
300 
301 #endif /* THEPEG_XSecStat_H */