Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:19:17 

0001 #include <cassert>
0002 #include <cfloat>
0003 
0004 #include "JetMETCorrections/FFTJetObjects/interface/FFTGenericScaleCalculator.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "DataFormats/JetReco/interface/PFJet.h"
0007 
0008 #define int_param(varname) m_##varname(ps.getParameter<int>(#varname))
0009 
0010 #define check_param(varname)                                                          \
0011   if ((m_##varname) >= 0) {                                                           \
0012     if ((m_##varname) >= nFactors)                                                    \
0013       throw cms::Exception("FFTJetBadConfig")                                         \
0014           << "In FFTGenericScaleCalculator constructor: "                             \
0015           << "out of range mapping for variable \"" << #varname << "\"" << std::endl; \
0016     mask[(m_##varname)] = 1;                                                          \
0017     ++dim;                                                                            \
0018   }
0019 
0020 static inline double delPhi(const double phi1, const double phi2) {
0021   double dphi = phi1 - phi2;
0022   if (dphi > M_PI)
0023     dphi -= 2.0 * M_PI;
0024   else if (dphi < -M_PI)
0025     dphi += 2.0 * M_PI;
0026   return dphi;
0027 }
0028 
0029 FFTGenericScaleCalculator::FFTGenericScaleCalculator(const edm::ParameterSet& ps)
0030     : m_factors(ps.getParameter<std::vector<double> >("factors")),
0031       m_minLog(ps.getUntrackedParameter<double>("minLog", -800.0)),
0032       int_param(eta),
0033       int_param(phi),
0034       int_param(pt),
0035       int_param(logPt),
0036       int_param(mass),
0037       int_param(logMass),
0038       int_param(energy),
0039       int_param(logEnergy),
0040       int_param(gamma),
0041       int_param(logGamma),
0042       int_param(pileup),
0043       int_param(ncells),
0044       int_param(etSum),
0045       int_param(etaWidth),
0046       int_param(phiWidth),
0047       int_param(averageWidth),
0048       int_param(widthRatio),
0049       int_param(etaPhiCorr),
0050       int_param(fuzziness),
0051       int_param(convergenceDistance),
0052       int_param(recoScale),
0053       int_param(recoScaleRatio),
0054       int_param(membershipFactor),
0055       int_param(magnitude),
0056       int_param(logMagnitude),
0057       int_param(magS1),
0058       int_param(LogMagS1),
0059       int_param(magS2),
0060       int_param(LogMagS2),
0061       int_param(driftSpeed),
0062       int_param(magSpeed),
0063       int_param(lifetime),
0064       int_param(splitTime),
0065       int_param(mergeTime),
0066       int_param(scale),
0067       int_param(logScale),
0068       int_param(nearestNeighborDistance),
0069       int_param(clusterRadius),
0070       int_param(clusterSeparation),
0071       int_param(dRFromJet),
0072       int_param(LaplacianS1),
0073       int_param(LaplacianS2),
0074       int_param(LaplacianS3),
0075       int_param(HessianS2),
0076       int_param(HessianS4),
0077       int_param(HessianS6),
0078       int_param(nConstituents),
0079       int_param(aveConstituentPt),
0080       int_param(logAveConstituentPt),
0081       int_param(constituentPtDistribution),
0082       int_param(constituentEtaPhiSpread),
0083       int_param(chargedHadronEnergyFraction),
0084       int_param(neutralHadronEnergyFraction),
0085       int_param(photonEnergyFraction),
0086       int_param(electronEnergyFraction),
0087       int_param(muonEnergyFraction),
0088       int_param(HFHadronEnergyFraction),
0089       int_param(HFEMEnergyFraction),
0090       int_param(chargedHadronMultiplicity),
0091       int_param(neutralHadronMultiplicity),
0092       int_param(photonMultiplicity),
0093       int_param(electronMultiplicity),
0094       int_param(muonMultiplicity),
0095       int_param(HFHadronMultiplicity),
0096       int_param(HFEMMultiplicity),
0097       int_param(chargedEmEnergyFraction),
0098       int_param(chargedMuEnergyFraction),
0099       int_param(neutralEmEnergyFraction),
0100       int_param(EmEnergyFraction),
0101       int_param(chargedMultiplicity),
0102       int_param(neutralMultiplicity) {
0103   const int nFactors = m_factors.size();
0104   std::vector<int> mask(nFactors, 0);
0105   int dim = 0;
0106 
0107   check_param(eta);
0108   check_param(phi);
0109   check_param(pt);
0110   check_param(logPt);
0111   check_param(mass);
0112   check_param(logMass);
0113   check_param(energy);
0114   check_param(logEnergy);
0115   check_param(gamma);
0116   check_param(logGamma);
0117   check_param(pileup);
0118   check_param(ncells);
0119   check_param(etSum);
0120   check_param(etaWidth);
0121   check_param(phiWidth);
0122   check_param(averageWidth);
0123   check_param(widthRatio);
0124   check_param(etaPhiCorr);
0125   check_param(fuzziness);
0126   check_param(convergenceDistance);
0127   check_param(recoScale);
0128   check_param(recoScaleRatio);
0129   check_param(membershipFactor);
0130   check_param(magnitude);
0131   check_param(logMagnitude);
0132   check_param(magS1);
0133   check_param(LogMagS1);
0134   check_param(magS2);
0135   check_param(LogMagS2);
0136   check_param(driftSpeed);
0137   check_param(magSpeed);
0138   check_param(lifetime);
0139   check_param(splitTime);
0140   check_param(mergeTime);
0141   check_param(scale);
0142   check_param(logScale);
0143   check_param(nearestNeighborDistance);
0144   check_param(clusterRadius);
0145   check_param(clusterSeparation);
0146   check_param(dRFromJet);
0147   check_param(LaplacianS1);
0148   check_param(LaplacianS2);
0149   check_param(LaplacianS3);
0150   check_param(HessianS2);
0151   check_param(HessianS4);
0152   check_param(HessianS6);
0153   check_param(nConstituents);
0154   check_param(aveConstituentPt);
0155   check_param(logAveConstituentPt);
0156   check_param(constituentPtDistribution);
0157   check_param(constituentEtaPhiSpread);
0158   check_param(chargedHadronEnergyFraction);
0159   check_param(neutralHadronEnergyFraction);
0160   check_param(photonEnergyFraction);
0161   check_param(electronEnergyFraction);
0162   check_param(muonEnergyFraction);
0163   check_param(HFHadronEnergyFraction);
0164   check_param(HFEMEnergyFraction);
0165   check_param(chargedHadronMultiplicity);
0166   check_param(neutralHadronMultiplicity);
0167   check_param(photonMultiplicity);
0168   check_param(electronMultiplicity);
0169   check_param(muonMultiplicity);
0170   check_param(HFHadronMultiplicity);
0171   check_param(HFEMMultiplicity);
0172   check_param(chargedEmEnergyFraction);
0173   check_param(chargedMuEnergyFraction);
0174   check_param(neutralEmEnergyFraction);
0175   check_param(EmEnergyFraction);
0176   check_param(chargedMultiplicity);
0177   check_param(neutralMultiplicity);
0178 
0179   if (dim != nFactors)
0180     throw cms::Exception("FFTJetBadConfig")
0181         << "In FFTGenericScaleCalculator constructor: "
0182         << "incompatible number of scaling factors: expected " << dim << ", got " << nFactors << std::endl;
0183   for (int i = 0; i < nFactors; ++i)
0184     if (mask[i] == 0)
0185       throw cms::Exception("FFTJetBadConfig") << "In FFTGenericScaleCalculator constructor: "
0186                                               << "variable number " << i << " is not mapped" << std::endl;
0187 }
0188 
0189 void FFTGenericScaleCalculator::mapFFTJet(const reco::Jet& jet,
0190                                           const reco::FFTJet<float>& fftJet,
0191                                           const math::XYZTLorentzVector& current,
0192                                           double* buf,
0193                                           const unsigned dim) const {
0194   // Verify that the input is reasonable
0195   if (dim != m_factors.size())
0196     throw cms::Exception("FFTJetBadConfig")
0197         << "In FFTGenericScaleCalculator::mapFFTJet: "
0198         << "incompatible table dimensionality: expected " << m_factors.size() << ", got " << dim << std::endl;
0199   if (dim)
0200     assert(buf);
0201   else
0202     return;
0203 
0204   // Go over all variables and map them as configured.
0205   // Variables from the "current" Lorentz vector.
0206   if (m_eta >= 0)
0207     buf[m_eta] = current.eta();
0208 
0209   if (m_phi >= 0)
0210     buf[m_phi] = current.phi();
0211 
0212   if (m_pt >= 0)
0213     buf[m_pt] = current.pt();
0214 
0215   if (m_logPt >= 0)
0216     buf[m_logPt] = f_safeLog(current.pt());
0217 
0218   if (m_mass >= 0)
0219     buf[m_mass] = current.M();
0220 
0221   if (m_logMass >= 0)
0222     buf[m_mass] = f_safeLog(current.M());
0223 
0224   if (m_energy >= 0)
0225     buf[m_energy] = current.e();
0226 
0227   if (m_logEnergy >= 0)
0228     buf[m_energy] = f_safeLog(current.e());
0229 
0230   if (m_gamma >= 0) {
0231     const double m = current.M();
0232     if (m > 0.0)
0233       buf[m_gamma] = current.e() / m;
0234     else
0235       buf[m_gamma] = DBL_MAX;
0236   }
0237 
0238   if (m_logGamma >= 0) {
0239     const double m = current.M();
0240     if (m > 0.0)
0241       buf[m_gamma] = current.e() / m;
0242     else
0243       buf[m_gamma] = DBL_MAX;
0244     buf[m_gamma] = log(buf[m_gamma]);
0245   }
0246 
0247   // Variables from fftJet
0248   if (m_pileup >= 0)
0249     buf[m_pileup] = fftJet.f_pileup().pt();
0250 
0251   if (m_ncells >= 0)
0252     buf[m_ncells] = fftJet.f_ncells();
0253 
0254   if (m_etSum >= 0)
0255     buf[m_etSum] = fftJet.f_etSum();
0256 
0257   if (m_etaWidth >= 0)
0258     buf[m_etaWidth] = fftJet.f_etaWidth();
0259 
0260   if (m_phiWidth >= 0)
0261     buf[m_phiWidth] = fftJet.f_phiWidth();
0262 
0263   if (m_averageWidth >= 0) {
0264     const double etaw = fftJet.f_etaWidth();
0265     const double phiw = fftJet.f_phiWidth();
0266     buf[m_averageWidth] = sqrt(etaw * etaw + phiw * phiw);
0267   }
0268 
0269   if (m_widthRatio >= 0) {
0270     const double etaw = fftJet.f_etaWidth();
0271     const double phiw = fftJet.f_phiWidth();
0272     if (phiw > 0.0)
0273       buf[m_widthRatio] = etaw / phiw;
0274     else
0275       buf[m_widthRatio] = DBL_MAX;
0276   }
0277 
0278   if (m_etaPhiCorr >= 0)
0279     buf[m_etaPhiCorr] = fftJet.f_etaPhiCorr();
0280 
0281   if (m_fuzziness >= 0)
0282     buf[m_fuzziness] = fftJet.f_fuzziness();
0283 
0284   if (m_convergenceDistance >= 0)
0285     buf[m_convergenceDistance] = fftJet.f_convergenceDistance();
0286 
0287   if (m_recoScale >= 0)
0288     buf[m_recoScale] = fftJet.f_recoScale();
0289 
0290   if (m_recoScaleRatio >= 0)
0291     buf[m_recoScaleRatio] = fftJet.f_recoScaleRatio();
0292 
0293   if (m_membershipFactor >= 0)
0294     buf[m_membershipFactor] = fftJet.f_membershipFactor();
0295 
0296   // Get most often used precluster quantities
0297   const reco::PattRecoPeak<float>& preclus = fftJet.f_precluster();
0298   const double scale = preclus.scale();
0299 
0300   if (m_magnitude >= 0)
0301     buf[m_magnitude] = preclus.magnitude();
0302 
0303   if (m_logMagnitude >= 0)
0304     buf[m_logMagnitude] = f_safeLog(preclus.magnitude());
0305 
0306   if (m_magS1 >= 0)
0307     buf[m_magS1] = preclus.magnitude() * scale;
0308 
0309   if (m_LogMagS1 >= 0)
0310     buf[m_LogMagS1] = f_safeLog(preclus.magnitude() * scale);
0311 
0312   if (m_magS2 >= 0)
0313     buf[m_magS2] = preclus.magnitude() * scale * scale;
0314 
0315   if (m_LogMagS2 >= 0)
0316     buf[m_LogMagS2] = f_safeLog(preclus.magnitude() * scale * scale);
0317 
0318   if (m_driftSpeed >= 0)
0319     buf[m_driftSpeed] = preclus.driftSpeed();
0320 
0321   if (m_magSpeed >= 0)
0322     buf[m_magSpeed] = preclus.magSpeed();
0323 
0324   if (m_lifetime >= 0)
0325     buf[m_lifetime] = preclus.lifetime();
0326 
0327   if (m_splitTime >= 0)
0328     buf[m_splitTime] = preclus.splitTime();
0329 
0330   if (m_mergeTime >= 0)
0331     buf[m_mergeTime] = preclus.mergeTime();
0332 
0333   if (m_scale >= 0)
0334     buf[m_scale] = scale;
0335 
0336   if (m_logScale >= 0)
0337     buf[m_logScale] = f_safeLog(scale);
0338 
0339   if (m_nearestNeighborDistance >= 0)
0340     buf[m_nearestNeighborDistance] = preclus.nearestNeighborDistance();
0341 
0342   if (m_clusterRadius >= 0)
0343     buf[m_clusterRadius] = preclus.clusterRadius();
0344 
0345   if (m_clusterSeparation >= 0)
0346     buf[m_clusterSeparation] = preclus.clusterSeparation();
0347 
0348   if (m_dRFromJet >= 0) {
0349     const double deta = preclus.eta() - current.eta();
0350     const double dphi = delPhi(preclus.phi(), current.phi());
0351     buf[m_dRFromJet] = sqrt(deta * deta + dphi * dphi);
0352   }
0353 
0354   if (m_LaplacianS1 >= 0) {
0355     double h[3];
0356     preclus.hessian(h);
0357     buf[m_LaplacianS1] = fabs(h[0] + h[2]) * scale;
0358   }
0359 
0360   if (m_LaplacianS2 >= 0) {
0361     double h[3];
0362     preclus.hessian(h);
0363     buf[m_LaplacianS2] = fabs(h[0] + h[2]) * scale * scale;
0364   }
0365 
0366   if (m_LaplacianS3 >= 0) {
0367     double h[3];
0368     preclus.hessian(h);
0369     buf[m_LaplacianS3] = fabs(h[0] + h[2]) * scale * scale * scale;
0370   }
0371 
0372   if (m_HessianS2 >= 0) {
0373     double h[3];
0374     preclus.hessian(h);
0375     buf[m_HessianS2] = fabs(h[0] * h[2] - h[1] * h[1]) * scale * scale;
0376   }
0377 
0378   if (m_HessianS4 >= 0) {
0379     double h[3];
0380     preclus.hessian(h);
0381     buf[m_HessianS4] = fabs(h[0] * h[2] - h[1] * h[1]) * pow(scale, 4);
0382   }
0383 
0384   if (m_HessianS6 >= 0) {
0385     double h[3];
0386     preclus.hessian(h);
0387     buf[m_HessianS6] = fabs(h[0] * h[2] - h[1] * h[1]) * pow(scale, 6);
0388   }
0389 
0390   // Variables from reco::Jet
0391   if (m_nConstituents >= 0)
0392     buf[m_nConstituents] = jet.nConstituents();
0393 
0394   if (m_aveConstituentPt >= 0)
0395     buf[m_aveConstituentPt] = current.pt() / jet.nConstituents();
0396 
0397   if (m_logAveConstituentPt >= 0)
0398     buf[m_logAveConstituentPt] = f_safeLog(current.pt() / jet.nConstituents());
0399 
0400   if (m_constituentPtDistribution >= 0)
0401     buf[m_constituentPtDistribution] = jet.constituentPtDistribution();
0402 
0403   if (m_constituentEtaPhiSpread >= 0)
0404     buf[m_constituentEtaPhiSpread] = jet.constituentEtaPhiSpread();
0405 
0406   // Variables from reco::PFJet
0407   const reco::PFJet* pfjet = dynamic_cast<const reco::PFJet*>(&jet);
0408   if (pfjet) {
0409     // Particle flow jet
0410     if (m_chargedHadronEnergyFraction >= 0)
0411       buf[m_chargedHadronEnergyFraction] = pfjet->chargedHadronEnergyFraction();
0412 
0413     if (m_neutralHadronEnergyFraction >= 0)
0414       buf[m_neutralHadronEnergyFraction] = pfjet->neutralHadronEnergyFraction();
0415 
0416     if (m_photonEnergyFraction >= 0)
0417       buf[m_photonEnergyFraction] = pfjet->photonEnergyFraction();
0418 
0419     if (m_electronEnergyFraction >= 0)
0420       buf[m_electronEnergyFraction] = pfjet->electronEnergyFraction();
0421 
0422     if (m_muonEnergyFraction >= 0)
0423       buf[m_muonEnergyFraction] = pfjet->muonEnergyFraction();
0424 
0425     if (m_HFHadronEnergyFraction >= 0)
0426       buf[m_HFHadronEnergyFraction] = pfjet->HFHadronEnergyFraction();
0427 
0428     if (m_HFEMEnergyFraction >= 0)
0429       buf[m_HFEMEnergyFraction] = pfjet->HFEMEnergyFraction();
0430 
0431     if (m_chargedHadronMultiplicity >= 0)
0432       buf[m_chargedHadronMultiplicity] = pfjet->chargedHadronMultiplicity();
0433 
0434     if (m_neutralHadronMultiplicity >= 0)
0435       buf[m_neutralHadronMultiplicity] = pfjet->neutralHadronMultiplicity();
0436 
0437     if (m_photonMultiplicity >= 0)
0438       buf[m_photonMultiplicity] = pfjet->photonMultiplicity();
0439 
0440     if (m_electronMultiplicity >= 0)
0441       buf[m_electronMultiplicity] = pfjet->electronMultiplicity();
0442 
0443     if (m_muonMultiplicity >= 0)
0444       buf[m_muonMultiplicity] = pfjet->muonMultiplicity();
0445 
0446     if (m_HFHadronMultiplicity >= 0)
0447       buf[m_HFHadronMultiplicity] = pfjet->HFHadronMultiplicity();
0448 
0449     if (m_HFEMMultiplicity >= 0)
0450       buf[m_HFEMMultiplicity] = pfjet->HFEMMultiplicity();
0451 
0452     if (m_chargedEmEnergyFraction >= 0)
0453       buf[m_chargedEmEnergyFraction] = pfjet->chargedEmEnergyFraction();
0454 
0455     if (m_chargedMuEnergyFraction >= 0)
0456       buf[m_chargedMuEnergyFraction] = pfjet->chargedMuEnergyFraction();
0457 
0458     if (m_neutralEmEnergyFraction >= 0)
0459       buf[m_neutralEmEnergyFraction] = pfjet->neutralEmEnergyFraction();
0460 
0461     if (m_EmEnergyFraction >= 0)
0462       buf[m_EmEnergyFraction] = pfjet->neutralEmEnergyFraction() + pfjet->chargedEmEnergyFraction();
0463 
0464     if (m_chargedMultiplicity >= 0)
0465       buf[m_chargedMultiplicity] = pfjet->chargedMultiplicity();
0466 
0467     if (m_neutralMultiplicity >= 0)
0468       buf[m_neutralMultiplicity] = pfjet->neutralMultiplicity();
0469   } else {
0470     // Not a particle flow jet
0471     if (m_chargedHadronEnergyFraction >= 0 || m_neutralHadronEnergyFraction >= 0 || m_photonEnergyFraction >= 0 ||
0472         m_electronEnergyFraction >= 0 || m_muonEnergyFraction >= 0 || m_HFHadronEnergyFraction >= 0 ||
0473         m_HFEMEnergyFraction >= 0 || m_chargedHadronMultiplicity >= 0 || m_neutralHadronMultiplicity >= 0 ||
0474         m_photonMultiplicity >= 0 || m_electronMultiplicity >= 0 || m_muonMultiplicity >= 0 ||
0475         m_HFHadronMultiplicity >= 0 || m_HFEMMultiplicity >= 0 || m_chargedEmEnergyFraction >= 0 ||
0476         m_chargedMuEnergyFraction >= 0 || m_neutralEmEnergyFraction >= 0 || m_EmEnergyFraction >= 0 ||
0477         m_chargedMultiplicity >= 0 || m_neutralMultiplicity >= 0)
0478       throw cms::Exception("FFTJetBadConfig") << "In FFTGenericScaleCalculator::mapFFTJet: "
0479                                               << "this configuration is valid for particle flow jets only" << std::endl;
0480   }
0481 
0482   // Apply the scaling factors
0483   for (unsigned i = 0; i < dim; ++i)
0484     buf[i] *= m_factors[i];
0485 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!