[Rivet] NLOHisto1D

[Andy Buckley]

On 20/05/15 15:36, Leif Lönnblad wrote:
> 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.

This is what the "multi-event" branch does, along with handling of
weight vectors within each (sub)event. Take a look in
/hepforge/hg/rivet/private/multi-event  I'll let Chris say more about
the status of this -- as far as I know it's working for 1D histograms
and just needs a bit of technical work to generalize to wrappers on the
other YODA data types without masses of unmaintainable code duplication
-- i.e. some template or macro smartness, preferably without being too
clever-clever.

> 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).

Yes, this would be nice. As has been said in later threads, I think it
needs to be switchable at the Rivet level, but exactly where the "smooth
filling" code lives is up for debate. I guess it needs to be
runtime-switchable (as opposed to a compile option) and applied
globally, since users need to be able to compare the results of their
Rivet runs on different systems/builds, and analyses won't know if
they'll be receiving counter-events in a given run.

Aaaaand back to paternity-leave radio-silence ;-)

Andy



>   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);
> 
>   }
> 
> _______________________________________________
> Rivet mailing list
> Rivet at projects.hepforge.org
> https://www.hepforge.org/lists/listinfo/rivet


-- 
Dr Andy Buckley, Lecturer / Royal Society University Research Fellow
Particle Physics Expt Group, University of Glasgow