LHCbCommon.hh

#ifndef PROJECTIONS_LHCBCOMMON_HH
#define PROJECTIONS_LHCBCOMMON_HH
 
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/PrimaryParticles.hh"
#include "Rivet/Projections/UnstableParticles.hh"
 
 
namespace Rivet {
 
  namespace LHCb {

    class HVRecoFinalState : public FinalState {
    protected:
 
      Rivet::Particles _theHFHadrons;
      Cut _hfcuts;
 
    public:
 
      HVRecoFinalState(const Cut& c = Cuts::abseta > 2.0 && Cuts::abseta < 5.0,
                       const vector<PdgId>& Hadrons = {},
                       const Cut& hadcut = Cuts::OPEN)
          : FinalState(c), _theHFHadrons(), _hfcuts(hadcut) {
        setName("LHCb::HVRecoFinalState");
        // accept only charged heavy-flavour hadrons and prepare Cut for UnstableParticles projection
        Cut cpid = Cuts::OPEN;
        for (const PdgId& id : Hadrons) {
          if (!PID::isHadron(id)) continue;
          if (PID::threeCharge(id) == 0) continue;
          if (!PID::isHeavyFlavour(id)) continue;
          if (cpid == Cuts::OPEN)
            cpid = (Cuts::pid == id);
          else
            cpid = cpid || (Cuts::pid == id);
        };
        if ((hadcut == Cuts::OPEN) && (c != Cuts::OPEN)) _hfcuts = c;
        if (cpid != Cuts::OPEN) _hfcuts = cpid && _hfcuts;
        const UnstableParticles hvhad(_hfcuts);
        this->declare(hvhad, "HVHadrons");
      }
 
      virtual ~HVRecoFinalState() { }
 
      virtual void project(const Event& e) {
        _theHFHadrons.clear();
        _theHFHadrons = apply<UnstableParticles>(e, "HVHadrons").particles();
        if (_theHFHadrons.empty()) { //event needs to be vetoed in analyze if no HV hadrons found
          _theParticles.clear();
          return;
        }
        MSG_DEBUG("Found " << _theHFHadrons.size() << " heavy-flavour hadrons in event");
        FinalState::project(e);
        if (_theParticles.empty()) {
          MSG_DEBUG("No particles in final state.");
          return;
        }
        MSG_DEBUG("Final state contains " << _theParticles.size()
                                          << " particles before adding heavy-flavour hadrons");
        // add heavy flavour hadrons to final state and remove descendants (avoid double-counting)
        _theParticles = _theParticles + _theHFHadrons;
        MSG_DEBUG("Added heavy-flavour hadrons to final state, now " << _theParticles.size()
                                                                     << " particles in final state");
        Particles to_remove = {};
        for (const Particle& p : _theHFHadrons) {
          to_remove = to_remove + p.stableDescendants();
        };
        MSG_DEBUG("Marked " << to_remove.size()
                            << " stable descendants of heavy-flavour hadrons for removal from final state");
        // eliminate descendants based on GenParticlePtr (better match?)
        idiscard(_theParticles, [&](Particle& p) {
          bool found = false;
          for (const Particle& r : to_remove)
            if (p.genParticle() == r.genParticle()) {
              found = true;
              break;
            };
          return found;
        });
        MSG_DEBUG("Removed stable descendants of heavy-flavour hadrons from final state, now "
                  << _theParticles.size() << " particles in final state. Project done!");
        to_remove.clear();
        return;
      }

      Particles hfHadrons() const {
        return _theHFHadrons;
      }

      virtual std::unique_ptr<Rivet::Projection> clone() const {
        return std::unique_ptr<Projection>(new HVRecoFinalState(*this));
      }

      using Projection::operator=;

      virtual CmpState compare(const Projection& p) const {
        if (FinalState::compare(p) == CmpState::NEQ) return CmpState::NEQ;
        const HVRecoFinalState& other = dynamic_cast<const HVRecoFinalState&>(p);
        if (hasProjection("HVHadrons") != other.hasProjection("HVHadrons")) return CmpState::NEQ;
        if (hasProjection("HVHadrons")) {
          const PCmp prevcmp = mkPCmp(other, "HVHadrons");
          if (prevcmp != CmpState::EQ) return CmpState::NEQ;
        }
        const bool hfcutcmp = _hfcuts == other._hfcuts;
        MSG_TRACE(_hfcuts << " VS " << other._hfcuts << " -> EQ == " << std::boolalpha << hfcutcmp);
        if (!hfcutcmp) return CmpState::NEQ;
        // Checks all passed: these FSes are equivalent
        return CmpState::EQ;
      }
    };
 
 
  }
}
 
#endif