QCDSmear.cc

#include "QCDSmear.h"
#include "PythiaJetGun.hh"
#include "rk/rk.hh"


namespace stab {
    bool QCDSmear::passesCut(const double prob, const unsigned recoNum,
                                  const unsigned ihard, const unsigned isoft,
                                  const PseudoJet& hard,
                                  const PseudoJet& soft) const
    {
        const double softPt = soft.pt();
#ifdef STREAMVERB
        std::cout << "QCDCut " << prob/alphaS_.alphaS(softPt*softPt) << std::endl;
#endif
        return prob/alphaS_.alphaS(softPt*softPt) >= cutoff_;
    }

    bool QCDSmear::operator()(const double prob, const unsigned recoNum,
                                   const unsigned ihard, const unsigned isoft,
                                   PseudoJet& hard, PseudoJet& soft)
    {
        if (passesCut(prob, recoNum, ihard, isoft, hard, soft))
        {
        const rk::P4 hardP4(geom3::Vector3(hard.px(),hard.py(),hard.pz()),hard.m());
        const rk::P4 softP4(geom3::Vector3(soft.px(),soft.py(),soft.pz()),soft.m());
            // These particles are coming from the common
            // mother with sufficiently high probability.
            // Fluctuate the energy flow.
            const rk::P4& mother = hardP4 + softP4;
            const double rnd0 = PythiaJetGun::pyrandom();
            rk::P4 dau1, dau2;

            if (uniformPhaseSpace_)
            {
                // Redecay according to simple phase space
                const rk::P4 parent(mother.momentum(), mother.m());
                const double rnd1 = PythiaJetGun::pyrandom();
                rk::phaseSpaceDecay(parent, hard.m(), soft.m(),
                                    rnd0, rnd1, &dau1, &dau2);
            }
            else
            {
                // Rotate the decay around the boost axis by
                // a random angle
                const rk::Boost& motherBoost = mother.restBoost();
                const rk::P4& boostedHard = motherBoost*hardP4;
                const rk::P4& boostedSoft = motherBoost*hardP4;
                geom3::Rotation3 rot(motherBoost.direction(), 2.0*M_PI*rnd0);
                dau1 = rk::P4(rot*boostedHard.momentum(), hard.m());
                dau2 = rk::P4(rot*boostedSoft.momentum(), soft.m());

                // Boost back to the lab
                const rk::Boost& labBoost = motherBoost.inverse();
                dau1.boost(labBoost);
                dau2.boost(labBoost);
            }

            if (dau1.pt() < dau2.pt())
            {
                hard = PseudoJet(dau2.px(),dau2.py(),dau2.pz(), dau2.e());
                soft = PseudoJet(dau1.px(),dau1.py(),dau1.pz(), dau1.e());
            }
            else
            {
                hard = PseudoJet(dau1.px(),dau1.py(),dau1.pz(), dau1.e());
                soft = PseudoJet(dau2.px(),dau2.py(),dau2.pz(), dau2.e());
            }
        }
        return true;
    }

} // end stab namespace