[Rivet] NLOHisto1D

[Leif Lönnblad]

On 2015-05-20 13:52, Frank Siegert wrote:
> 1) Do we want to always use this in all analyses, i.e. do we accept
> the (probably minor) overhead from this bookkeeping? I don't know
> whether it would be an option to introduce a switch which dynamically
> enables this functionality within the histogram class.

Since this is determined by the input events, I guess it is impossible 
to say beforehand which analyses should have it or not. So yes, it 
should be used in all analysis. Always? Well, there is currently no way 
of knowing if a HepMC file has counter events (it is just flagged by two 
consecutive events having the same event number) I suggest we do it always.

> 2) Should this live in a Rivet wrapper around YODA::Histo1D or is this
> more a functionality of the histogram class that should be moved into
> YODA? Of course, in the latter case we would not fill using
> Rivet::Event but just event weight and event number.

I think this should live in Rivet. Cf. also my previous comment about 
only allowing filling of histograms in the analyse function, which would 
also require a wrapper around YODA histos.

We have also discussed to use smooth filling of histogram bins, where 
not only one histogram bin is filled, but depending on how close a point 
is to the bib border, also the adjacent bin would be filled (with 
modified weights). This is to get correct statistics when an event is 
fills close the bin border, and the counter event fills close to the 
adjacent one. (I include code for this below).

Cheers,
Leif



   void Histo1D::sfill(double x, double weight) {
     if ( std::isnan(x) ) throw RangeError("X is NaN");
     if ( std::isinf(x) ) throw RangeError("X is Inf");
     // Fill the overall distribution
     _axis.totalDbn().fill(x, weight);

     // We assume that the weight is actually a rectangular
     // distribution of weights around the given x-value. This means
     // that if the x-value is not exactly at the corresponding bin
     // midpoint only a part of the weight is filled in this bin and
     // with a slightly shifted x-value. The rest of the weight is
     // filled in an adjacent bin using an x-value which is the
     // midpoint of the overflowing rectangular distribution into this
     // bin. If the adjacent bin is outside the histogram, the
     // rectangular distribution entered completely in the bin
     // corresponding to the given x-value.

     //    try {
     if ( binIndexAt(x) != -1 ) {
       // First we get bin corresponding to the x-value (if none was
       // found, go directly to the overflow fill).
       Bin & bin = binAt(x);

       // How far from the midpoint are we?
       double dx = x - bin.midpoint();

       // Calculate the part of the weight corresponding to an adjacent
       // bin, and the correspondingly shifted x-value.
       double wadj = weight*abs(dx)/bin.width();
       double xadj = dx > 0? bin.highEdge() + dx*0.5: bin.lowEdge() + 
dx*0.5;

       // and what is left in the current bin.
       double ww = weight - wadj;
       double xx = x - dx*0.5;

       // First fill the adjacent bin
       //      try {
       if ( binIndexAt(xadj) != -1 ) {
	binAt(xadj).fill(xadj, wadj);
	//      } catch (const RangeError& re) {
       } else {
	// If that failed the main bin should be filled with the whole
	// weight distribution.
	ww = weight;
	xx = bin.midpoint();
       }
       bin.fill(xx, ww);
       //    } catch (const RangeError& re) {
     } else {
       // Overflow is handled as in a normal histogram.
       if      (x <  _axis.xMin()) _axis.underflow().fill(x, weight);
       else if (x >= _axis.xMax()) _axis.overflow().fill(x, weight);
     }

     // Lock the axis now that a fill has happened
     _axis._setLock(true);

   }