AbsSmearJet.cc

#include <cmath>
#include <cassert>

#include "AbsSmearJet.hh"
#include "PythiaJetGun.hh"

#define CODE_GAMMA 22

using namespace PythiaJetGun;

namespace stab {
    double AbsSmearJet::calorimeterEtUncertainty(
        const double S, const double N, const double C, const double Et)
    {
        assert(Et >= 0.0);
        return sqrt(Et*S*S + N*N + Et*Et*C*C);
    }

    PseudoJet AbsSmearJet::smear4Mom(const PseudoJet& p4, const double sigmaPt,
                                  const double sigmaEtaPhi,
                                  const bool makePtPositive)
    {
        const double originalPt = p4.pt();
        double pt = originalPt + gaussRandom(0.0, sigmaPt);
        if (pt < 0.0)
            pt = 0.0;
        while (makePtPositive && pt <= 0.0)
            pt = originalPt + gaussRandom(0.0, sigmaPt);
        const double eta = p4.eta() + gaussRandom(0.0, sigmaEtaPhi);
        const double phi = p4.phi() + gaussRandom(0.0, sigmaEtaPhi);
        const double pE = p4.E();
    const double m = p4.m();
    const double px = pt*cos(phi);
    const double py = pt*sin(phi);
    const double pz = pt*sinh(eta);
    const double E = sqrt(m*m+pt*pt+pz*pz);
        PseudoJet smearedPart(px, py, pz, E);
        //fastjet::PseudoJet::UserInfoBase info =  p4.user_info<StabilityUserInfo>();
    //smearedPart.set_user_info( &info );
    return smearedPart;
    }

    void AbsSmearJet::smear(Event & smearedEvent) 
    {
        const unsigned nPart = m_event.size();

        smearedEvent.reserve(nPart);

        for (unsigned i=0; i<nPart; ++i)
        {
            const PseudoJet & part = m_event[i];
        //const StabilityUserInfo & info = part.user_info<StabilityUserInfo>();
        const StabilityUserInfo & info = m_info[i];
        const int code = info.pdgId();
        const double charge = info.charge();
            if (code == CODE_GAMMA)
                processPhoton(part, info, &smearedEvent);
            else if (charge == 0.0)
                // These are typically neutrons (code 2112),
                // antineutrons (code -2112), and K_L (code 130)
                processNeutralHadron(part,info, &smearedEvent);
            else
                processTrack(part, info, &smearedEvent);
        }
    assert( smearedEvent.size() );
    }
}