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 */