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
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
0205
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
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
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
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
0407 const reco::PFJet* pfjet = dynamic_cast<const reco::PFJet*>(&jet);
0408 if (pfjet) {
0409
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
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
0483 for (unsigned i = 0; i < dim; ++i)
0484 buf[i] *= m_factors[i];
0485 }
