QCDSmearSequence.cc

#include "QCDSmearSequence.h"
#include "mother_probability.hh"

namespace stab {

    QCDSmearSequence::QCDSmearSequence(const AlphaStrong& alphaS)
        : alphaS_(alphaS)
    {
    }

    QCDSmearSequence::QCDSmearSequence(const AlphaStrong& alphaS,
                     const std::vector<QCDSmearSequenceCallback*> cb)
        : alphaS_(alphaS), callbacks_(cb)
    {
        const unsigned nCallbacks = callbacks_.size();
        for (unsigned icall=0; icall<nCallbacks; ++icall)
            assert(callbacks_[icall]);
    }

    void QCDSmearSequence::smear(Event& smearedEvent)
    {
        m_callBackCount = 0;
        
        // copy event into smearedEvent
        const unsigned sz = m_event.size();
        //smearedEvent.clear();
    for ( unsigned i=0; i != sz; i++ )
      smearedEvent[i] = m_event[i];
    

        if (sz > 1U)
        {
            // Calculate the transverse momentum values
            ptValues_.clear();
            for (unsigned i=0; i<sz; ++i)
                ptValues_.push_back(smearedEvent[i].pt());

            // Calculate pairwise probabilities
            probs_.clear();
            const unsigned szm1 = sz - 1U;
            for (unsigned i=0; i<szm1; ++i)
            {
                const PseudoJet& pi = smearedEvent[i];
                for (unsigned j=i+1; j<sz; ++j)
                {
                    const double prob = mother_probability(pi, smearedEvent[j], alphaS_);
                    probs_.push_back(ProbTriple(prob, i, j));
                }
            }
            std::sort(probs_.begin(), probs_.end(), std::greater<ProbTriple>());

            // Go over all pairs calling the callbacks and marking
            // the "softer" of the two particles as "used"
            used_.clear();
            for (unsigned i=0; i<sz; ++i)
                used_.push_back(0);

            const unsigned nCallbacks = callbacks_.size();
            const unsigned long nprobs = probs_.size();
            for (unsigned long iprob=0; iprob<nprobs; ++iprob)
            {
                const ProbTriple& t(probs_[iprob]);
                unsigned soft, hard;

                if (ptValues_[t.second] < ptValues_[t.third])
                {
                    soft = t.second;
                    hard = t.third;
                }
                else
                {
                    soft = t.third;
                    hard = t.second;
                }
                if (!used_[soft])
                {
                    // Call various callbacks
                    bool continueIterating = true;
                    for (unsigned icall=0; icall<nCallbacks && continueIterating; ++icall)
                    {
                        continueIterating = (*callbacks_[icall])(
                            t.first, m_callBackCount, hard,
                            soft, smearedEvent[hard], smearedEvent[soft]);
                        if (continueIterating)
                            ++m_callBackCount;
                    }
                    used_[soft] = 1;
                    if (!continueIterating)
                        break;
                }
            }
        }

    }

} // end stab namespace