XSecStat.h

// -*- C++ -*-
//
// XSecStat.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_XSecStat_H
#define THEPEG_XSecStat_H
//
// This is the declaration of the XSecStat class.
//
 
#include "ThePEG/Config/ThePEG.h"
 
namespace ThePEG {
 
class XSecStat {
 
public:
 
  enum {
    plainWeights = 0,
    plainVetoedWeights,
    reweightedWeights,
    reweightedVetoedWeights
  };
 
  XSecStat() 
    : theMaxXSec(ZERO), theAttempts(0), theAccepted(0), theVetoed(0),
      theSumWeights (),
      theSumWeights2(), theLastWeight(0.0) {}
 
  explicit XSecStat(CrossSection xsecmax) 
    : theMaxXSec(xsecmax), theAttempts(0), theAccepted(0), theVetoed(0),
      theSumWeights (),
      theSumWeights2(), theLastWeight(0.0) {}
 
  XSecStat & operator=(const XSecStat & x) = default;
 
  XSecStat & operator+=(const XSecStat & x) {
    theAttempts    += x.theAttempts;
    theAccepted    += x.theAccepted;
    theVetoed      += x.theVetoed;
    for( unsigned int ix = 0; ix < 4; ++ix ) {
      theSumWeights [ix] += x.theSumWeights [ix];
      theSumWeights2[ix] += x.theSumWeights2[ix];
    }
    theLastWeight = 0.0;
    return *this;
  }
 
  void reset() {
    theAttempts = theAccepted = theVetoed = 0;
    theSumWeights = theSumWeights2 = {};
    theLastWeight = 0.0;
  }
 
 
public:
 
 
  void accept() { 
    theAccepted += 1;
  }
 
  void select(double weight) {
    theAttempts += 1;
    theSumWeights [reweightedWeights] +=     weight ;
    theSumWeights2[reweightedWeights] += sqr(weight);
    theSumWeights [plainWeights]      +=     weight ;
    theSumWeights2[plainWeights]      += sqr(weight);
    theLastWeight = weight;
  }
 
  void reweight(double oldWeight, double newWeight) {
    theSumWeights [reweightedWeights] +=     newWeight  -     oldWeight ;
    theSumWeights2[reweightedWeights] += sqr(newWeight) - sqr(oldWeight);
  }
 
  void reject(double weight = 1.0) {
    theSumWeights [reweightedVetoedWeights] +=            weight ;
    theSumWeights2[reweightedVetoedWeights] +=        sqr(weight);
    theSumWeights [plainVetoedWeights]      +=     theLastWeight ;
    theSumWeights2[plainVetoedWeights]      += sqr(theLastWeight);
    theVetoed += 1;
  }
 
  CrossSection maxXSec() const { return theMaxXSec; }
 
  double sumWeights() const { 
    return theSumWeights[reweightedWeights] - theSumWeights[reweightedVetoedWeights];
  }
 
  double sumWeights2() const {
    return theSumWeights2[reweightedWeights] + theSumWeights2[reweightedVetoedWeights];
  }
 
  double sumWeightsNoReweight() const { 
    return theSumWeights[plainWeights] - theSumWeights[plainVetoedWeights];
  }
 
  double sumWeights2NoReweight() const { 
    return theSumWeights2[plainWeights] + theSumWeights2[plainVetoedWeights];
  }
 
  CrossSection xSec(double att = 0) const {
    double n = (att == 0.0 ? attempts() : att);
    return n ? maxXSec()*sumWeights()/n : maxXSec();
  }
 
  CrossSection xSecErr(double att = 0) const {
    double n = (att == 0.0 ? attempts() : att);
    if ( n < 2 )
      return maxXSec();
    double sw = sumWeights(); double sw2 = sumWeights2();
    return
      maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
  }
 
  CrossSection xSecNoReweight(double att = 0) const {
    double n = (att == 0.0 ? attempts() : att);
    return n ? maxXSec()*sumWeightsNoReweight()/n : maxXSec();
  }
 
  CrossSection xSecErrNoReweight(double att = 0) const {
    double n = (att == 0.0 ? attempts() : att);
    if ( n < 2 )
      return maxXSec();
    double sw = sumWeightsNoReweight(); 
    double sw2 = sumWeights2NoReweight();
    return
      maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
  }
 
  double attempts() const { return theAttempts; }
 
  double accepted() const { return theAccepted-theVetoed; }
 
  double vetoed() const { return theVetoed; }
 
  void maxXSec(CrossSection x) { theMaxXSec = x; }
 
public:
 
  void output(PersistentOStream & os) const;
 
  void input(PersistentIStream & is);
 
private:
 
  CrossSection theMaxXSec;
 
  double theAttempts;
 
  double theAccepted;
 
  double theVetoed;
 
  array<double,4> theSumWeights;
 
  array<double,4> theSumWeights2;
 
  double theLastWeight;
 
};
 
PersistentOStream & operator<<(PersistentOStream &, const XSecStat &);
 
PersistentIStream & operator>>(PersistentIStream &, XSecStat &);
 
inline XSecStat operator+(const XSecStat & x1, const XSecStat & x2) {
  XSecStat x = x1;
  return x += x2;
}
 
}
 
#endif /* THEPEG_XSecStat_H */