FFTGenericScaleCalculator.cc

Go to the documentation of this file.

#include <cassert>
#include <cfloat>

#include "JetMETCorrections/FFTJetObjects/interface/FFTGenericScaleCalculator.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "DataFormats/JetReco/interface/PFJet.h"

#define int_param(varname) m_##varname(ps.getParameter<int>(#varname))

#define check_param(varname)                                                          \
  if ((m_##varname) >= 0) {                                                           \
    if ((m_##varname) >= nFactors)                                                    \
      throw cms::Exception("FFTJetBadConfig")                                         \
          << "In FFTGenericScaleCalculator constructor: "                             \
          << "out of range mapping for variable \"" << #varname << "\"" << std::endl; \
    mask[(m_##varname)] = 1;                                                          \
    ++dim;                                                                            \
  }

static inline double delPhi(const double phi1, const double phi2) {
  double dphi = phi1 - phi2;
  if (dphi > M_PI)
    dphi -= 2.0 * M_PI;
  else if (dphi < -M_PI)
    dphi += 2.0 * M_PI;
  return dphi;
}

FFTGenericScaleCalculator::FFTGenericScaleCalculator(const edm::ParameterSet& ps)
    : m_factors(ps.getParameter<std::vector<double> >("factors")),
      m_minLog(ps.getUntrackedParameter<double>("minLog", -800.0)),
      int_param(eta),
      int_param(phi),
      int_param(pt),
      int_param(logPt),
      int_param(mass),
      int_param(logMass),
      int_param(energy),
      int_param(logEnergy),
      int_param(gamma),
      int_param(logGamma),
      int_param(pileup),
      int_param(ncells),
      int_param(etSum),
      int_param(etaWidth),
      int_param(phiWidth),
      int_param(averageWidth),
      int_param(widthRatio),
      int_param(etaPhiCorr),
      int_param(fuzziness),
      int_param(convergenceDistance),
      int_param(recoScale),
      int_param(recoScaleRatio),
      int_param(membershipFactor),
      int_param(magnitude),
      int_param(logMagnitude),
      int_param(magS1),
      int_param(LogMagS1),
      int_param(magS2),
      int_param(LogMagS2),
      int_param(driftSpeed),
      int_param(magSpeed),
      int_param(lifetime),
      int_param(splitTime),
      int_param(mergeTime),
      int_param(scale),
      int_param(logScale),
      int_param(nearestNeighborDistance),
      int_param(clusterRadius),
      int_param(clusterSeparation),
      int_param(dRFromJet),
      int_param(LaplacianS1),
      int_param(LaplacianS2),
      int_param(LaplacianS3),
      int_param(HessianS2),
      int_param(HessianS4),
      int_param(HessianS6),
      int_param(nConstituents),
      int_param(aveConstituentPt),
      int_param(logAveConstituentPt),
      int_param(constituentPtDistribution),
      int_param(constituentEtaPhiSpread),
      int_param(chargedHadronEnergyFraction),
      int_param(neutralHadronEnergyFraction),
      int_param(photonEnergyFraction),
      int_param(electronEnergyFraction),
      int_param(muonEnergyFraction),
      int_param(HFHadronEnergyFraction),
      int_param(HFEMEnergyFraction),
      int_param(chargedHadronMultiplicity),
      int_param(neutralHadronMultiplicity),
      int_param(photonMultiplicity),
      int_param(electronMultiplicity),
      int_param(muonMultiplicity),
      int_param(HFHadronMultiplicity),
      int_param(HFEMMultiplicity),
      int_param(chargedEmEnergyFraction),
      int_param(chargedMuEnergyFraction),
      int_param(neutralEmEnergyFraction),
      int_param(EmEnergyFraction),
      int_param(chargedMultiplicity),
      int_param(neutralMultiplicity) {
  const int nFactors = m_factors.size();
  std::vector<int> mask(nFactors, 0);
  int dim = 0;

  check_param(eta);
  check_param(phi);
  check_param(pt);
  check_param(logPt);
  check_param(mass);
  check_param(logMass);
  check_param(energy);
  check_param(logEnergy);
  check_param(gamma);
  check_param(logGamma);
  check_param(pileup);
  check_param(ncells);
  check_param(etSum);
  check_param(etaWidth);
  check_param(phiWidth);
  check_param(averageWidth);
  check_param(widthRatio);
  check_param(etaPhiCorr);
  check_param(fuzziness);
  check_param(convergenceDistance);
  check_param(recoScale);
  check_param(recoScaleRatio);
  check_param(membershipFactor);
  check_param(magnitude);
  check_param(logMagnitude);
  check_param(magS1);
  check_param(LogMagS1);
  check_param(magS2);
  check_param(LogMagS2);
  check_param(driftSpeed);
  check_param(magSpeed);
  check_param(lifetime);
  check_param(splitTime);
  check_param(mergeTime);
  check_param(scale);
  check_param(logScale);
  check_param(nearestNeighborDistance);
  check_param(clusterRadius);
  check_param(clusterSeparation);
  check_param(dRFromJet);
  check_param(LaplacianS1);
  check_param(LaplacianS2);
  check_param(LaplacianS3);
  check_param(HessianS2);
  check_param(HessianS4);
  check_param(HessianS6);
  check_param(nConstituents);
  check_param(aveConstituentPt);
  check_param(logAveConstituentPt);
  check_param(constituentPtDistribution);
  check_param(constituentEtaPhiSpread);
  check_param(chargedHadronEnergyFraction);
  check_param(neutralHadronEnergyFraction);
  check_param(photonEnergyFraction);
  check_param(electronEnergyFraction);
  check_param(muonEnergyFraction);
  check_param(HFHadronEnergyFraction);
  check_param(HFEMEnergyFraction);
  check_param(chargedHadronMultiplicity);
  check_param(neutralHadronMultiplicity);
  check_param(photonMultiplicity);
  check_param(electronMultiplicity);
  check_param(muonMultiplicity);
  check_param(HFHadronMultiplicity);
  check_param(HFEMMultiplicity);
  check_param(chargedEmEnergyFraction);
  check_param(chargedMuEnergyFraction);
  check_param(neutralEmEnergyFraction);
  check_param(EmEnergyFraction);
  check_param(chargedMultiplicity);
  check_param(neutralMultiplicity);

  if (dim != nFactors)
    throw cms::Exception("FFTJetBadConfig")
        << "In FFTGenericScaleCalculator constructor: "
        << "incompatible number of scaling factors: expected " << dim << ", got " << nFactors << std::endl;
  for (int i = 0; i < nFactors; ++i)
    if (mask[i] == 0)
      throw cms::Exception("FFTJetBadConfig") << "In FFTGenericScaleCalculator constructor: "
                                              << "variable number " << i << " is not mapped" << std::endl;
}

void FFTGenericScaleCalculator::mapFFTJet(const reco::Jet& jet,
                                          const reco::FFTJet<float>& fftJet,
                                          const math::XYZTLorentzVector& current,
                                          double* buf,
                                          const unsigned dim) const {
  // Verify that the input is reasonable
  if (dim != m_factors.size())
    throw cms::Exception("FFTJetBadConfig")
        << "In FFTGenericScaleCalculator::mapFFTJet: "
        << "incompatible table dimensionality: expected " << m_factors.size() << ", got " << dim << std::endl;
  if (dim)
    assert(buf);
  else
    return;

  // Go over all variables and map them as configured.
  // Variables from the "current" Lorentz vector.
  if (m_eta >= 0)
    buf[m_eta] = current.eta();

  if (m_phi >= 0)
    buf[m_phi] = current.phi();

  if (m_pt >= 0)
    buf[m_pt] = current.pt();

  if (m_logPt >= 0)
    buf[m_logPt] = f_safeLog(current.pt());

  if (m_mass >= 0)
    buf[m_mass] = current.M();

  if (m_logMass >= 0)
    buf[m_mass] = f_safeLog(current.M());

  if (m_energy >= 0)
    buf[m_energy] = current.e();

  if (m_logEnergy >= 0)
    buf[m_energy] = f_safeLog(current.e());

  if (m_gamma >= 0) {
    const double m = current.M();
    if (m > 0.0)
      buf[m_gamma] = current.e() / m;
    else
      buf[m_gamma] = DBL_MAX;
  }

  if (m_logGamma >= 0) {
    const double m = current.M();
    if (m > 0.0)
      buf[m_gamma] = current.e() / m;
    else
      buf[m_gamma] = DBL_MAX;
    buf[m_gamma] = log(buf[m_gamma]);
  }

  // Variables from fftJet
  if (m_pileup >= 0)
    buf[m_pileup] = fftJet.f_pileup().pt();

  if (m_ncells >= 0)
    buf[m_ncells] = fftJet.f_ncells();

  if (m_etSum >= 0)
    buf[m_etSum] = fftJet.f_etSum();

  if (m_etaWidth >= 0)
    buf[m_etaWidth] = fftJet.f_etaWidth();

  if (m_phiWidth >= 0)
    buf[m_phiWidth] = fftJet.f_phiWidth();

  if (m_averageWidth >= 0) {
    const double etaw = fftJet.f_etaWidth();
    const double phiw = fftJet.f_phiWidth();
    buf[m_averageWidth] = sqrt(etaw * etaw + phiw * phiw);
  }

  if (m_widthRatio >= 0) {
    const double etaw = fftJet.f_etaWidth();
    const double phiw = fftJet.f_phiWidth();
    if (phiw > 0.0)
      buf[m_widthRatio] = etaw / phiw;
    else
      buf[m_widthRatio] = DBL_MAX;
  }

  if (m_etaPhiCorr >= 0)
    buf[m_etaPhiCorr] = fftJet.f_etaPhiCorr();

  if (m_fuzziness >= 0)
    buf[m_fuzziness] = fftJet.f_fuzziness();

  if (m_convergenceDistance >= 0)
    buf[m_convergenceDistance] = fftJet.f_convergenceDistance();

  if (m_recoScale >= 0)
    buf[m_recoScale] = fftJet.f_recoScale();

  if (m_recoScaleRatio >= 0)
    buf[m_recoScaleRatio] = fftJet.f_recoScaleRatio();

  if (m_membershipFactor >= 0)
    buf[m_membershipFactor] = fftJet.f_membershipFactor();

  // Get most often used precluster quantities
  const reco::PattRecoPeak<float>& preclus = fftJet.f_precluster();
  const double scale = preclus.scale();

  if (m_magnitude >= 0)
    buf[m_magnitude] = preclus.magnitude();

  if (m_logMagnitude >= 0)
    buf[m_logMagnitude] = f_safeLog(preclus.magnitude());

  if (m_magS1 >= 0)
    buf[m_magS1] = preclus.magnitude() * scale;

  if (m_LogMagS1 >= 0)
    buf[m_LogMagS1] = f_safeLog(preclus.magnitude() * scale);

  if (m_magS2 >= 0)
    buf[m_magS2] = preclus.magnitude() * scale * scale;

  if (m_LogMagS2 >= 0)
    buf[m_LogMagS2] = f_safeLog(preclus.magnitude() * scale * scale);

  if (m_driftSpeed >= 0)
    buf[m_driftSpeed] = preclus.driftSpeed();

  if (m_magSpeed >= 0)
    buf[m_magSpeed] = preclus.magSpeed();

  if (m_lifetime >= 0)
    buf[m_lifetime] = preclus.lifetime();

  if (m_splitTime >= 0)
    buf[m_splitTime] = preclus.splitTime();

  if (m_mergeTime >= 0)
    buf[m_mergeTime] = preclus.mergeTime();

  if (m_scale >= 0)
    buf[m_scale] = scale;

  if (m_logScale >= 0)
    buf[m_logScale] = f_safeLog(scale);

  if (m_nearestNeighborDistance >= 0)
    buf[m_nearestNeighborDistance] = preclus.nearestNeighborDistance();

  if (m_clusterRadius >= 0)
    buf[m_clusterRadius] = preclus.clusterRadius();

  if (m_clusterSeparation >= 0)
    buf[m_clusterSeparation] = preclus.clusterSeparation();

  if (m_dRFromJet >= 0) {
    const double deta = preclus.eta() - current.eta();
    const double dphi = delPhi(preclus.phi(), current.phi());
    buf[m_dRFromJet] = sqrt(deta * deta + dphi * dphi);
  }

  if (m_LaplacianS1 >= 0) {
    double h[3];
    preclus.hessian(h);
    buf[m_LaplacianS1] = fabs(h[0] + h[2]) * scale;
  }

  if (m_LaplacianS2 >= 0) {
    double h[3];
    preclus.hessian(h);
    buf[m_LaplacianS2] = fabs(h[0] + h[2]) * scale * scale;
  }

  if (m_LaplacianS3 >= 0) {
    double h[3];
    preclus.hessian(h);
    buf[m_LaplacianS3] = fabs(h[0] + h[2]) * scale * scale * scale;
  }

  if (m_HessianS2 >= 0) {
    double h[3];
    preclus.hessian(h);
    buf[m_HessianS2] = fabs(h[0] * h[2] - h[1] * h[1]) * scale * scale;
  }

  if (m_HessianS4 >= 0) {
    double h[3];
    preclus.hessian(h);
    buf[m_HessianS4] = fabs(h[0] * h[2] - h[1] * h[1]) * pow(scale, 4);
  }

  if (m_HessianS6 >= 0) {
    double h[3];
    preclus.hessian(h);
    buf[m_HessianS6] = fabs(h[0] * h[2] - h[1] * h[1]) * pow(scale, 6);
  }

  // Variables from reco::Jet
  if (m_nConstituents >= 0)
    buf[m_nConstituents] = jet.nConstituents();

  if (m_aveConstituentPt >= 0)
    buf[m_aveConstituentPt] = current.pt() / jet.nConstituents();

  if (m_logAveConstituentPt >= 0)
    buf[m_logAveConstituentPt] = f_safeLog(current.pt() / jet.nConstituents());

  if (m_constituentPtDistribution >= 0)
    buf[m_constituentPtDistribution] = jet.constituentPtDistribution();

  if (m_constituentEtaPhiSpread >= 0)
    buf[m_constituentEtaPhiSpread] = jet.constituentEtaPhiSpread();

  // Variables from reco::PFJet
  const reco::PFJet* pfjet = dynamic_cast<const reco::PFJet*>(&jet);
  if (pfjet) {
    // Particle flow jet
    if (m_chargedHadronEnergyFraction >= 0)
      buf[m_chargedHadronEnergyFraction] = pfjet->chargedHadronEnergyFraction();

    if (m_neutralHadronEnergyFraction >= 0)
      buf[m_neutralHadronEnergyFraction] = pfjet->neutralHadronEnergyFraction();

    if (m_photonEnergyFraction >= 0)
      buf[m_photonEnergyFraction] = pfjet->photonEnergyFraction();

    if (m_electronEnergyFraction >= 0)
      buf[m_electronEnergyFraction] = pfjet->electronEnergyFraction();

    if (m_muonEnergyFraction >= 0)
      buf[m_muonEnergyFraction] = pfjet->muonEnergyFraction();

    if (m_HFHadronEnergyFraction >= 0)
      buf[m_HFHadronEnergyFraction] = pfjet->HFHadronEnergyFraction();

    if (m_HFEMEnergyFraction >= 0)
      buf[m_HFEMEnergyFraction] = pfjet->HFEMEnergyFraction();

    if (m_chargedHadronMultiplicity >= 0)
      buf[m_chargedHadronMultiplicity] = pfjet->chargedHadronMultiplicity();

    if (m_neutralHadronMultiplicity >= 0)
      buf[m_neutralHadronMultiplicity] = pfjet->neutralHadronMultiplicity();

    if (m_photonMultiplicity >= 0)
      buf[m_photonMultiplicity] = pfjet->photonMultiplicity();

    if (m_electronMultiplicity >= 0)
      buf[m_electronMultiplicity] = pfjet->electronMultiplicity();

    if (m_muonMultiplicity >= 0)
      buf[m_muonMultiplicity] = pfjet->muonMultiplicity();

    if (m_HFHadronMultiplicity >= 0)
      buf[m_HFHadronMultiplicity] = pfjet->HFHadronMultiplicity();

    if (m_HFEMMultiplicity >= 0)
      buf[m_HFEMMultiplicity] = pfjet->HFEMMultiplicity();

    if (m_chargedEmEnergyFraction >= 0)
      buf[m_chargedEmEnergyFraction] = pfjet->chargedEmEnergyFraction();

    if (m_chargedMuEnergyFraction >= 0)
      buf[m_chargedMuEnergyFraction] = pfjet->chargedMuEnergyFraction();

    if (m_neutralEmEnergyFraction >= 0)
      buf[m_neutralEmEnergyFraction] = pfjet->neutralEmEnergyFraction();

    if (m_EmEnergyFraction >= 0)
      buf[m_EmEnergyFraction] = pfjet->neutralEmEnergyFraction() + pfjet->chargedEmEnergyFraction();

    if (m_chargedMultiplicity >= 0)
      buf[m_chargedMultiplicity] = pfjet->chargedMultiplicity();

    if (m_neutralMultiplicity >= 0)
      buf[m_neutralMultiplicity] = pfjet->neutralMultiplicity();
  } else {
    // Not a particle flow jet
    if (m_chargedHadronEnergyFraction >= 0 || m_neutralHadronEnergyFraction >= 0 || m_photonEnergyFraction >= 0 ||
        m_electronEnergyFraction >= 0 || m_muonEnergyFraction >= 0 || m_HFHadronEnergyFraction >= 0 ||
        m_HFEMEnergyFraction >= 0 || m_chargedHadronMultiplicity >= 0 || m_neutralHadronMultiplicity >= 0 ||
        m_photonMultiplicity >= 0 || m_electronMultiplicity >= 0 || m_muonMultiplicity >= 0 ||
        m_HFHadronMultiplicity >= 0 || m_HFEMMultiplicity >= 0 || m_chargedEmEnergyFraction >= 0 ||
        m_chargedMuEnergyFraction >= 0 || m_neutralEmEnergyFraction >= 0 || m_EmEnergyFraction >= 0 ||
        m_chargedMultiplicity >= 0 || m_neutralMultiplicity >= 0)
      throw cms::Exception("FFTJetBadConfig") << "In FFTGenericScaleCalculator::mapFFTJet: "
                                              << "this configuration is valid for particle flow jets only" << std::endl;
  }

  // Apply the scaling factors
  for (unsigned i = 0; i < dim; ++i)
    buf[i] *= m_factors[i];
}