00001 #include <cassert>
00002 #include <cfloat>
00003
00004 #include "JetMETCorrections/FFTJetObjects/interface/FFTGenericScaleCalculator.h"
00005 #include "FWCore/Utilities/interface/Exception.h"
00006 #include "DataFormats/JetReco/interface/PFJet.h"
00007
00008 #define int_param(varname) m_ ## varname (ps.getParameter< int >( #varname ))
00009
00010 #define check_param(varname) if ((m_ ## varname) >= 0) {\
00011 if ((m_ ## varname) >= nFactors)\
00012 throw cms::Exception("FFTJetBadConfig")\
00013 << "In FFTGenericScaleCalculator constructor: "\
00014 << "out of range mapping for variable \""\
00015 << #varname << "\"" << std::endl;\
00016 mask[(m_ ## varname)] = 1;\
00017 ++dim;\
00018 }
00019
00020
00021 static inline double delPhi(const double phi1, const double phi2)
00022 {
00023 double dphi = phi1 - phi2;
00024 if (dphi > M_PI)
00025 dphi -= 2.0*M_PI;
00026 else if (dphi < -M_PI)
00027 dphi += 2.0*M_PI;
00028 return dphi;
00029 }
00030
00031
00032 FFTGenericScaleCalculator::FFTGenericScaleCalculator(const edm::ParameterSet& ps)
00033 : m_factors(ps.getParameter<std::vector<double> >("factors")),
00034 m_minLog(ps.getUntrackedParameter<double>("minLog", -800.0)),
00035 int_param(eta),
00036 int_param(phi),
00037 int_param(pt),
00038 int_param(logPt),
00039 int_param(mass),
00040 int_param(logMass),
00041 int_param(energy),
00042 int_param(logEnergy),
00043 int_param(gamma),
00044 int_param(logGamma),
00045 int_param(pileup),
00046 int_param(ncells),
00047 int_param(etSum),
00048 int_param(etaWidth),
00049 int_param(phiWidth),
00050 int_param(averageWidth),
00051 int_param(widthRatio),
00052 int_param(etaPhiCorr),
00053 int_param(fuzziness),
00054 int_param(convergenceDistance),
00055 int_param(recoScale),
00056 int_param(recoScaleRatio),
00057 int_param(membershipFactor),
00058 int_param(magnitude),
00059 int_param(logMagnitude),
00060 int_param(magS1),
00061 int_param(LogMagS1),
00062 int_param(magS2),
00063 int_param(LogMagS2),
00064 int_param(driftSpeed),
00065 int_param(magSpeed),
00066 int_param(lifetime),
00067 int_param(scale),
00068 int_param(logScale),
00069 int_param(nearestNeighborDistance),
00070 int_param(clusterRadius),
00071 int_param(clusterSeparation),
00072 int_param(dRFromJet),
00073 int_param(LaplacianS1),
00074 int_param(LaplacianS2),
00075 int_param(LaplacianS3),
00076 int_param(HessianS2),
00077 int_param(HessianS4),
00078 int_param(HessianS6),
00079 int_param(nConstituents),
00080 int_param(aveConstituentPt),
00081 int_param(logAveConstituentPt),
00082 int_param(constituentPtDistribution),
00083 int_param(constituentEtaPhiSpread),
00084 int_param(chargedHadronEnergyFraction),
00085 int_param(neutralHadronEnergyFraction),
00086 int_param(photonEnergyFraction),
00087 int_param(electronEnergyFraction),
00088 int_param(muonEnergyFraction),
00089 int_param(HFHadronEnergyFraction),
00090 int_param(HFEMEnergyFraction),
00091 int_param(chargedHadronMultiplicity),
00092 int_param(neutralHadronMultiplicity),
00093 int_param(photonMultiplicity),
00094 int_param(electronMultiplicity),
00095 int_param(muonMultiplicity),
00096 int_param(HFHadronMultiplicity),
00097 int_param(HFEMMultiplicity),
00098 int_param(chargedEmEnergyFraction),
00099 int_param(chargedMuEnergyFraction),
00100 int_param(neutralEmEnergyFraction),
00101 int_param(EmEnergyFraction),
00102 int_param(chargedMultiplicity),
00103 int_param(neutralMultiplicity)
00104 {
00105 const int nFactors = m_factors.size();
00106 std::vector<int> mask(nFactors, 0);
00107 int dim = 0;
00108
00109 check_param(eta);
00110 check_param(phi);
00111 check_param(pt);
00112 check_param(logPt);
00113 check_param(mass);
00114 check_param(logMass);
00115 check_param(energy);
00116 check_param(logEnergy);
00117 check_param(gamma);
00118 check_param(logGamma);
00119 check_param(pileup);
00120 check_param(ncells);
00121 check_param(etSum);
00122 check_param(etaWidth);
00123 check_param(phiWidth);
00124 check_param(averageWidth);
00125 check_param(widthRatio);
00126 check_param(etaPhiCorr);
00127 check_param(fuzziness);
00128 check_param(convergenceDistance);
00129 check_param(recoScale);
00130 check_param(recoScaleRatio);
00131 check_param(membershipFactor);
00132 check_param(magnitude);
00133 check_param(logMagnitude);
00134 check_param(magS1);
00135 check_param(LogMagS1);
00136 check_param(magS2);
00137 check_param(LogMagS2);
00138 check_param(driftSpeed);
00139 check_param(magSpeed);
00140 check_param(lifetime);
00141 check_param(scale);
00142 check_param(logScale);
00143 check_param(nearestNeighborDistance);
00144 check_param(clusterRadius);
00145 check_param(clusterSeparation);
00146 check_param(dRFromJet);
00147 check_param(LaplacianS1);
00148 check_param(LaplacianS2);
00149 check_param(LaplacianS3);
00150 check_param(HessianS2);
00151 check_param(HessianS4);
00152 check_param(HessianS6);
00153 check_param(nConstituents);
00154 check_param(aveConstituentPt);
00155 check_param(logAveConstituentPt);
00156 check_param(constituentPtDistribution);
00157 check_param(constituentEtaPhiSpread);
00158 check_param(chargedHadronEnergyFraction);
00159 check_param(neutralHadronEnergyFraction);
00160 check_param(photonEnergyFraction);
00161 check_param(electronEnergyFraction);
00162 check_param(muonEnergyFraction);
00163 check_param(HFHadronEnergyFraction);
00164 check_param(HFEMEnergyFraction);
00165 check_param(chargedHadronMultiplicity);
00166 check_param(neutralHadronMultiplicity);
00167 check_param(photonMultiplicity);
00168 check_param(electronMultiplicity);
00169 check_param(muonMultiplicity);
00170 check_param(HFHadronMultiplicity);
00171 check_param(HFEMMultiplicity);
00172 check_param(chargedEmEnergyFraction);
00173 check_param(chargedMuEnergyFraction);
00174 check_param(neutralEmEnergyFraction);
00175 check_param(EmEnergyFraction);
00176 check_param(chargedMultiplicity);
00177 check_param(neutralMultiplicity);
00178
00179 if (dim != nFactors)
00180 throw cms::Exception("FFTJetBadConfig")
00181 << "In FFTGenericScaleCalculator constructor: "
00182 << "incompatible number of scaling factors: expected "
00183 << dim << ", got " << nFactors << std::endl;
00184 for (int i=0; i<nFactors; ++i)
00185 if (mask[i] == 0)
00186 throw cms::Exception("FFTJetBadConfig")
00187 << "In FFTGenericScaleCalculator constructor: "
00188 << "variable number " << i << " is not mapped" << std::endl;
00189 }
00190
00191 void FFTGenericScaleCalculator::mapFFTJet(
00192 const reco::Jet& jet, const reco::FFTJet<float>& fftJet,
00193 const math::XYZTLorentzVector& current,
00194 double* buf, const unsigned dim) const
00195 {
00196
00197 if (dim != m_factors.size())
00198 throw cms::Exception("FFTJetBadConfig")
00199 << "In FFTGenericScaleCalculator::mapFFTJet: "
00200 << "incompatible table dimensionality: expected "
00201 << m_factors.size() << ", got " << dim << std::endl;
00202 if (dim)
00203 assert(buf);
00204 else
00205 return;
00206
00207
00208
00209 if (m_eta >= 0)
00210 buf[m_eta] = current.eta();
00211
00212 if (m_phi >= 0)
00213 buf[m_phi] = current.phi();
00214
00215 if (m_pt >= 0)
00216 buf[m_pt] = current.pt();
00217
00218 if (m_logPt >= 0)
00219 buf[m_logPt] = f_safeLog(current.pt());
00220
00221 if (m_mass >= 0)
00222 buf[m_mass] = current.M();
00223
00224 if (m_logMass >= 0)
00225 buf[m_mass] = f_safeLog(current.M());
00226
00227 if (m_energy >= 0)
00228 buf[m_energy] = current.e();
00229
00230 if (m_logEnergy >= 0)
00231 buf[m_energy] = f_safeLog(current.e());
00232
00233 if (m_gamma >= 0)
00234 {
00235 const double m = current.M();
00236 if (m > 0.0)
00237 buf[m_gamma] = current.e()/m;
00238 else
00239 buf[m_gamma] = DBL_MAX;
00240 }
00241
00242 if (m_logGamma >= 0)
00243 {
00244 const double m = current.M();
00245 if (m > 0.0)
00246 buf[m_gamma] = current.e()/m;
00247 else
00248 buf[m_gamma] = DBL_MAX;
00249 buf[m_gamma] = log(buf[m_gamma]);
00250 }
00251
00252
00253 if (m_pileup >= 0)
00254 buf[m_pileup] = fftJet.f_pileup().pt();
00255
00256 if (m_ncells >= 0)
00257 buf[m_ncells] = fftJet.f_ncells();
00258
00259 if (m_etSum >= 0)
00260 buf[m_etSum] = fftJet.f_etSum();
00261
00262 if (m_etaWidth >= 0)
00263 buf[m_etaWidth] = fftJet.f_etaWidth();
00264
00265 if (m_phiWidth >= 0)
00266 buf[m_phiWidth] = fftJet.f_phiWidth();
00267
00268 if (m_averageWidth >= 0)
00269 {
00270 const double etaw = fftJet.f_etaWidth();
00271 const double phiw = fftJet.f_phiWidth();
00272 buf[m_averageWidth] = sqrt(etaw*etaw + phiw*phiw);
00273 }
00274
00275 if (m_widthRatio >= 0)
00276 {
00277 const double etaw = fftJet.f_etaWidth();
00278 const double phiw = fftJet.f_phiWidth();
00279 if (phiw > 0.0)
00280 buf[m_widthRatio] = etaw/phiw;
00281 else
00282 buf[m_widthRatio] = DBL_MAX;
00283 }
00284
00285 if (m_etaPhiCorr >= 0)
00286 buf[m_etaPhiCorr] = fftJet.f_etaPhiCorr();
00287
00288 if (m_fuzziness >= 0)
00289 buf[m_fuzziness] = fftJet.f_fuzziness();
00290
00291 if (m_convergenceDistance >= 0)
00292 buf[m_convergenceDistance] = fftJet.f_convergenceDistance();
00293
00294 if (m_recoScale >= 0)
00295 buf[m_recoScale] = fftJet.f_recoScale();
00296
00297 if (m_recoScaleRatio >= 0)
00298 buf[m_recoScaleRatio] = fftJet.f_recoScaleRatio();
00299
00300 if (m_membershipFactor >= 0)
00301 buf[m_membershipFactor] = fftJet.f_membershipFactor();
00302
00303
00304 const reco::PattRecoPeak<float>& preclus = fftJet.f_precluster();
00305 const double scale = preclus.scale();
00306
00307 if (m_magnitude >= 0)
00308 buf[m_magnitude] = preclus.magnitude();
00309
00310 if (m_logMagnitude >= 0)
00311 buf[m_logMagnitude] = f_safeLog(preclus.magnitude());
00312
00313 if (m_magS1 >= 0)
00314 buf[m_magS1] = preclus.magnitude()*scale;
00315
00316 if (m_LogMagS1 >= 0)
00317 buf[m_LogMagS1] = f_safeLog(preclus.magnitude()*scale);
00318
00319 if (m_magS2 >= 0)
00320 buf[m_magS2] = preclus.magnitude()*scale*scale;
00321
00322 if (m_LogMagS2 >= 0)
00323 buf[m_LogMagS2] = f_safeLog(preclus.magnitude()*scale*scale);
00324
00325 if (m_driftSpeed >= 0)
00326 buf[m_driftSpeed] = preclus.driftSpeed();
00327
00328 if (m_magSpeed >= 0)
00329 buf[m_magSpeed] = preclus.magSpeed();
00330
00331 if (m_lifetime >= 0)
00332 buf[m_lifetime] = preclus.lifetime();
00333
00334 if (m_scale >= 0)
00335 buf[m_scale] = scale;
00336
00337 if (m_logScale >= 0)
00338 buf[m_logScale] = f_safeLog(scale);
00339
00340 if (m_nearestNeighborDistance >= 0)
00341 buf[m_nearestNeighborDistance] = preclus.nearestNeighborDistance();
00342
00343 if (m_clusterRadius >= 0)
00344 buf[m_clusterRadius] = preclus.clusterRadius();
00345
00346 if (m_clusterSeparation >= 0)
00347 buf[m_clusterSeparation] = preclus.clusterSeparation();
00348
00349 if (m_dRFromJet >= 0)
00350 {
00351 const double deta = preclus.eta() - current.eta();
00352 const double dphi = delPhi(preclus.phi(), current.phi());
00353 buf[m_dRFromJet] = sqrt(deta*deta + dphi*dphi);
00354 }
00355
00356 if (m_LaplacianS1 >= 0)
00357 {
00358 double h[3];
00359 preclus.hessian(h);
00360 buf[m_LaplacianS1] = fabs(h[0] + h[2])*scale;
00361 }
00362
00363 if (m_LaplacianS2 >= 0)
00364 {
00365 double h[3];
00366 preclus.hessian(h);
00367 buf[m_LaplacianS2] = fabs(h[0] + h[2])*scale*scale;
00368 }
00369
00370 if (m_LaplacianS3 >= 0)
00371 {
00372 double h[3];
00373 preclus.hessian(h);
00374 buf[m_LaplacianS3] = fabs(h[0] + h[2])*scale*scale*scale;
00375 }
00376
00377 if (m_HessianS2 >= 0)
00378 {
00379 double h[3];
00380 preclus.hessian(h);
00381 buf[m_HessianS2] = fabs(h[0]*h[2] - h[1]*h[1])*scale*scale;
00382 }
00383
00384 if (m_HessianS4 >= 0)
00385 {
00386 double h[3];
00387 preclus.hessian(h);
00388 buf[m_HessianS4] = fabs(h[0]*h[2] - h[1]*h[1])*pow(scale, 4);
00389 }
00390
00391 if (m_HessianS6 >= 0)
00392 {
00393 double h[3];
00394 preclus.hessian(h);
00395 buf[m_HessianS6] = fabs(h[0]*h[2] - h[1]*h[1])*pow(scale, 6);
00396 }
00397
00398
00399 if (m_nConstituents >= 0)
00400 buf[m_nConstituents] = jet.nConstituents();
00401
00402 if (m_aveConstituentPt >= 0)
00403 buf[m_aveConstituentPt] = current.pt()/jet.nConstituents();
00404
00405 if (m_logAveConstituentPt >= 0)
00406 buf[m_logAveConstituentPt] = f_safeLog(current.pt()/jet.nConstituents());
00407
00408 if (m_constituentPtDistribution >= 0)
00409 buf[m_constituentPtDistribution] = jet.constituentPtDistribution();
00410
00411 if (m_constituentEtaPhiSpread >= 0)
00412 buf[m_constituentEtaPhiSpread] = jet.constituentEtaPhiSpread();
00413
00414
00415 const reco::PFJet* pfjet = dynamic_cast<const reco::PFJet*>(&jet);
00416 if (pfjet)
00417 {
00418
00419 if (m_chargedHadronEnergyFraction >= 0)
00420 buf[m_chargedHadronEnergyFraction] = pfjet->chargedHadronEnergyFraction();
00421
00422 if (m_neutralHadronEnergyFraction >= 0)
00423 buf[m_neutralHadronEnergyFraction] = pfjet->neutralHadronEnergyFraction();
00424
00425 if (m_photonEnergyFraction >= 0)
00426 buf[m_photonEnergyFraction] = pfjet->photonEnergyFraction();
00427
00428 if (m_electronEnergyFraction >= 0)
00429 buf[m_electronEnergyFraction] = pfjet->electronEnergyFraction();
00430
00431 if (m_muonEnergyFraction >= 0)
00432 buf[m_muonEnergyFraction] = pfjet->muonEnergyFraction();
00433
00434 if (m_HFHadronEnergyFraction >= 0)
00435 buf[m_HFHadronEnergyFraction] = pfjet->HFHadronEnergyFraction();
00436
00437 if (m_HFEMEnergyFraction >= 0)
00438 buf[m_HFEMEnergyFraction] = pfjet->HFEMEnergyFraction();
00439
00440 if (m_chargedHadronMultiplicity >= 0)
00441 buf[m_chargedHadronMultiplicity] = pfjet->chargedHadronMultiplicity();
00442
00443 if (m_neutralHadronMultiplicity >= 0)
00444 buf[m_neutralHadronMultiplicity] = pfjet->neutralHadronMultiplicity();
00445
00446 if (m_photonMultiplicity >= 0)
00447 buf[m_photonMultiplicity] = pfjet->photonMultiplicity();
00448
00449 if (m_electronMultiplicity >= 0)
00450 buf[m_electronMultiplicity] = pfjet->electronMultiplicity();
00451
00452 if (m_muonMultiplicity >= 0)
00453 buf[m_muonMultiplicity] = pfjet->muonMultiplicity();
00454
00455 if (m_HFHadronMultiplicity >= 0)
00456 buf[m_HFHadronMultiplicity] = pfjet->HFHadronMultiplicity();
00457
00458 if (m_HFEMMultiplicity >= 0)
00459 buf[m_HFEMMultiplicity] = pfjet->HFEMMultiplicity();
00460
00461 if (m_chargedEmEnergyFraction >= 0)
00462 buf[m_chargedEmEnergyFraction] = pfjet->chargedEmEnergyFraction();
00463
00464 if (m_chargedMuEnergyFraction >= 0)
00465 buf[m_chargedMuEnergyFraction] = pfjet->chargedMuEnergyFraction();
00466
00467 if (m_neutralEmEnergyFraction >= 0)
00468 buf[m_neutralEmEnergyFraction] = pfjet->neutralEmEnergyFraction();
00469
00470 if (m_EmEnergyFraction >= 0)
00471 buf[m_EmEnergyFraction] = pfjet->neutralEmEnergyFraction() +
00472 pfjet->chargedEmEnergyFraction();
00473
00474 if (m_chargedMultiplicity >= 0)
00475 buf[m_chargedMultiplicity] = pfjet->chargedMultiplicity();
00476
00477 if (m_neutralMultiplicity >= 0)
00478 buf[m_neutralMultiplicity] = pfjet->neutralMultiplicity();
00479 }
00480 else
00481 {
00482
00483 if (m_chargedHadronEnergyFraction >= 0 ||
00484 m_neutralHadronEnergyFraction >= 0 ||
00485 m_photonEnergyFraction >= 0 ||
00486 m_electronEnergyFraction >= 0 ||
00487 m_muonEnergyFraction >= 0 ||
00488 m_HFHadronEnergyFraction >= 0 ||
00489 m_HFEMEnergyFraction >= 0 ||
00490 m_chargedHadronMultiplicity >= 0 ||
00491 m_neutralHadronMultiplicity >= 0 ||
00492 m_photonMultiplicity >= 0 ||
00493 m_electronMultiplicity >= 0 ||
00494 m_muonMultiplicity >= 0 ||
00495 m_HFHadronMultiplicity >= 0 ||
00496 m_HFEMMultiplicity >= 0 ||
00497 m_chargedEmEnergyFraction >= 0 ||
00498 m_chargedMuEnergyFraction >= 0 ||
00499 m_neutralEmEnergyFraction >= 0 ||
00500 m_EmEnergyFraction >= 0 ||
00501 m_chargedMultiplicity >= 0 ||
00502 m_neutralMultiplicity >= 0)
00503 throw cms::Exception("FFTJetBadConfig")
00504 << "In FFTGenericScaleCalculator::mapFFTJet: "
00505 << "this configuration is valid for particle flow jets only"
00506 << std::endl;
00507 }
00508
00509
00510 for (unsigned i=0; i<dim; ++i)
00511 buf[i] *= m_factors[i];
00512 }
1.7.3