Utilities/Statistics/EventContribution.h

0001 // -*- C++ -*-
0002 //
0003 // EventContribution.hpp is a part of myStatistics
0004 // Copyright (C) 2012-2019 Simon Platzer, The Herwig Collaboration
0005 //
0006 // myStatistics is licenced under version 3 of the GPL, see COPYING for details.
0007 //
0008 #ifndef MYSTATISTICS_EventContribution_hpp_included
0009 #define MYSTATISTICS_EventContribution_hpp_included
0010
0011 #include <cfloat>
0012
0013 namespace Statistics {
0014
0015   /**
0016    * \brief A pointlike or boxlike eventContribution; serves to define the EventContribution concept
0017    * \author Simon Platzer
0018    */
0019   class EventContribution {
0020
0021   public:
0022
0023     /**
0024      * Construct an eventContribution with given width
0025      */
0026     EventContribution(double newCentralValue,
0027               double newWeight,
0028               double newWidth = 0.0)
0029       : theCentralValue(newCentralValue),
0030     theSupport(newCentralValue - newWidth/2.,
0031            newCentralValue + newWidth/2.),
0032     theWeight(newWeight) {}
0033
0034   public:
0035
0036     /**
0037      * Return the central value
0038      */
0039     double centralValue() const { return theCentralValue; }
0040
0041     /**
0042      * Return the support
0043      */
0044     const std::pair<double,double>& support() const { return theSupport; }
0045
0046     /**
0047      * Return the normalized overlap with an interval
0048      */
0049     double overlap(const std::pair<double,double>& interval) const {
0050       return calculateOverlap(interval,support(),support().second-support().first);
0051     }
0052
0053     /**
0054      * Return the eventContribution weight
0055      */
0056     double weight() const { return theWeight; }
0057
0058   public:
0059
0060     /**
0061      * Remap the central value and support given a periodicity
0062      * interval; if the support exceeds the periodicity it is ajusted
0063      * to the peridocity interval.
0064      */
0065     void periodic(const std::pair<double,double>& periodicity) {
0066       double delta = periodicity.second - periodicity.first;
0067       if ( support().second - support().first > delta ) {
0068     theSupport.first = centralValue() - delta/2.;
0069     theSupport.second = centralValue() + delta/2.;
0070       }
0071       double shift = 0.;
0072       if ( centralValue() >= periodicity.second ) {
0073     while ( centralValue() >= periodicity.second ) {
0074       shift -= delta;
0075       theCentralValue -= delta;
0076     }
0077       } else if ( centralValue() < periodicity.first ) {
0078     while ( centralValue() <= periodicity.first ) {
0079       shift += delta;
0080       theCentralValue += delta;
0081     }
0082       }
0083       theSupport.first += shift;
0084       theSupport.second += shift;
0085       if ( theSupport.first < periodicity.first )
0086     theSupport.first += delta;
0087       if ( theSupport.second > periodicity.second )
0088     theSupport.second -= delta;
0089     }
0090
0091     /**
0092      * Adjust to lower boundary
0093      */
0094     void noUnderflow(double lower) {
0095       if ( support().first < lower ) {
0096     double shift = lower - support().first;
0097     theSupport.first += shift;
0098     theSupport.second += shift;
0099     theCentralValue += shift;
0100       }
0101     }
0102
0103     /**
0104      * Adjust to upper boundary
0105      */
0106     void noOverflow(double upper) {
0107       if ( support().second > upper ) {
0108     double shift = upper - support().second - DBL_EPSILON;
0109     theSupport.first += shift;
0110     theSupport.second += shift;
0111     theCentralValue += shift;
0112       }
0113     }
0114
0115     /**
0116      * Return the normalized overlap with an interval, assuming a
0117      * periodic quantity
0118      */
0119     double overlap(const std::pair<double,double>& interval,
0120            const std::pair<double,double>& periodicity) const {
0121       if ( support().first < support().second )
0122     return calculateOverlap(interval,support(),support().second-support().first);
0123       double norm =
0124     support().second - periodicity.first +
0125     periodicity.second - support().first;
0126       return
0127     calculateOverlap(interval,std::make_pair(periodicity.first,support().second),norm) +
0128     calculateOverlap(interval,std::make_pair(support().first,periodicity.second),norm);
0129     }
0130
0131   private:
0132
0133     /**
0134      * Calculate the normalized overlap with an interval
0135      */
0136     double calculateOverlap(const std::pair<double,double>& interval,
0137                 const std::pair<double,double>& newSupport,
0138                 double norm) const {
0139       if ( newSupport.first == newSupport.second ) {
0140     if ( newSupport.first >= interval.first &&
0141          newSupport.second < interval.second )
0142       return 1.0;
0143     return 0.0;
0144       }
0145       if ( newSupport.first <= interval.first &&
0146        newSupport.second <= interval.first )
0147     return 0.0;
0148       if ( newSupport.first >= interval.second &&
0149        newSupport.second >= interval.second )
0150     return 0.0;
0151       double lower = std::max(newSupport.first,interval.first);
0152       double upper = std::min(newSupport.second,interval.second);
0153       return std::max(0.0,(upper-lower)/norm);
0154     }
0155
0156     /**
0157      * The central value
0158      */
0159     double theCentralValue;
0160
0161     /**
0162      * The support
0163      */
0164     std::pair<double,double> theSupport;
0165
0166     /**
0167      * The eventContribution weight
0168      */
0169     double theWeight;
0170
0171   };
0172
0173 }
0174
0175 #endif // MYSTATISTICS_EventContribution_hpp_included