TransverseFuzzinessCalculator.cc

#include <cmath>
#include <cassert>

#include "TransverseFuzzinessCalculator.h"

namespace stab {
    double TransverseFuzzinessCalculator::calculate(
        const Event &initialParticles, 
        const Event &finalJets,
        const double *weights,
        const unsigned nJets, const unsigned nParts, 
        const bool lastFinalJetIsUnclustered,
        double *jetFuzziness, const unsigned lenJetFuzziness) const
    {
        assert(weights);
        assert(nJets);
        assert(nParts);
        assert(initialParticles.size() >= nParts);
        assert(finalJets.size() >= nJets);
        assert(jetFuzziness);
        assert(lenJetFuzziness >= nJets);

        if (ptStorage_.size() < nParts)
            ptStorage_.resize(nParts);
        double* ptValues = &ptStorage_[0];
        long double ptSum = 0.0L;
        for (unsigned i=0; i<nParts; ++i)
        {
            const double pt = initialParticles[i].pt();
            ptValues[i] = pt;
            ptSum += pt;
        }
        assert(ptSum > 0.0L);

        double eventFuzzinesContribution;
        long double fuzzySum = 0.0L;
        for (unsigned i=0; i<nJets; ++i)
        {
            jetFuzziness[i] = 0.0;
            eventFuzzinesContribution = 0.0;
            singleJetFuzziness(ptValues, weights + i*nParts, nParts,
                               jetFuzziness+i, &eventFuzzinesContribution);
            if ( addInQuadrature_ ) fuzzySum += eventFuzzinesContribution*eventFuzzinesContribution;
            else fuzzySum += eventFuzzinesContribution;
        }

        if (divideByNjets_) 
            fuzzySum /= nJets;
    if ( addInQuadrature_)
       fuzzySum = std::sqrt(fuzzySum);
        if (divideByPt_)
            fuzzySum /= ptSum;

        return fuzzySum;
    }

    void TransverseFuzzinessCalculator::singleJetFuzziness(
        const double *ptValues, const double *weights, const unsigned len,
        double* thisJetFuzzines, double* eventFuzzinesContribution) const
    {
        long double numSum = 0.0L, denomSum = 0.0L;

    double sumW = 0.0;

        for (unsigned i=0; i<len; ++i)
        {
            const double w = weights[i];
        sumW += w;
            if (w > 0.0)
            {
        assert(ptValues[i]>0.0L);
                assert(w <= 1.0);
                const double weightedPt = w*ptValues[i];
                denomSum += weightedPt;
                numSum += weightedPt*(1.0 - w)*ptValues[i];
        assert(denomSum > 0.0L);
            }
        }
    if ( sumW == 0. ) { 
      *eventFuzzinesContribution=0.;
      return ;
    }
    assert( sumW > 0. );
        assert(denomSum > 0.0L);

        const double eUncert = sqrt(static_cast<double>(numSum));
        *thisJetFuzzines = eUncert/static_cast<double>(denomSum);

        if (divideByPt_)
            *eventFuzzinesContribution = eUncert;
        else
            *eventFuzzinesContribution = *thisJetFuzzines;
    }
}