Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:35:07

0001 #include "RecoJets/JetProducers/interface/MVAJetPuId.h"
0002 #include "DataFormats/JetReco/interface/PFJet.h"
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "DataFormats/JetReco/interface/Jet.h"
0005 #include "DataFormats/VertexReco/interface/Vertex.h"
0006 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 #include "FWCore/ParameterSet/interface/FileInPath.h"
0009 #include "CommonTools/MVAUtils/interface/TMVAZipReader.h"
0010 #include "TMatrixDSym.h"
0011 #include "TMatrixDSymEigen.h"
0012 #include "DataFormats/JetReco/interface/PileupJetIdentifier.h"
0013 
0014 const float large_val = std::numeric_limits<float>::max();
0015 
0016 MVAJetPuId::MVAJetPuId(const edm::ParameterSet &ps) {
0017   impactParTkThreshod_ = 1.;
0018 
0019   tmvaWeights_ = edm::FileInPath(ps.getParameter<std::string>("tmvaWeights")).fullPath();
0020   tmvaMethod_ = ps.getParameter<std::string>("tmvaMethod");
0021   tmvaVariables_ = ps.getParameter<std::vector<std::string> >("tmvaVariables");
0022   tmvaSpectators_ = ps.getParameter<std::vector<std::string> >("tmvaSpectators");
0023   version_ = ps.getParameter<int>("version");
0024   reader_ = nullptr;
0025   edm::ParameterSet jetConfig = ps.getParameter<edm::ParameterSet>("JetIdParams");
0026   for (int i0 = 0; i0 < NWPs; i0++) {
0027     std::string lCutType = "Tight";
0028     if (i0 == PileupJetIdentifier::kMedium)
0029       lCutType = "Medium";
0030     if (i0 == PileupJetIdentifier::kLoose)
0031       lCutType = "Loose";
0032     for (int i1 = 0; i1 < 1; i1++) {
0033       std::vector<double> pt010 = jetConfig.getParameter<std::vector<double> >(("Pt010_" + lCutType).c_str());
0034       std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_" + lCutType).c_str());
0035       std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_" + lCutType).c_str());
0036       std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_" + lCutType).c_str());
0037       for (int i2 = 0; i2 < NPts; i2++)
0038         mvacut_[i0][0][i2] = pt010[i2];
0039       for (int i2 = 0; i2 < NPts; i2++)
0040         mvacut_[i0][1][i2] = pt1020[i2];
0041       for (int i2 = 0; i2 < NPts; i2++)
0042         mvacut_[i0][2][i2] = pt2030[i2];
0043       for (int i2 = 0; i2 < NPts; i2++)
0044         mvacut_[i0][3][i2] = pt3050[i2];
0045     }
0046   }
0047   setup();
0048 }
0049 
0050 MVAJetPuId::MVAJetPuId(int version,
0051                        const std::string &tmvaWeights,
0052                        const std::string &tmvaMethod,
0053                        Float_t impactParTkThreshod,
0054                        const std::vector<std::string> &tmvaVariables) {
0055   impactParTkThreshod_ = impactParTkThreshod;
0056   tmvaWeights_ = tmvaWeights;
0057   tmvaMethod_ = tmvaMethod;
0058   tmvaVariables_ = tmvaVariables;
0059   version_ = version;
0060 
0061   reader_ = nullptr;
0062 
0063   setup();
0064 }
0065 
0066 void MVAJetPuId::setup() {
0067   initVariables();
0068 
0069   tmvaVariables_.clear();
0070   tmvaVariables_.push_back("rho");
0071   tmvaVariables_.push_back("nParticles");
0072   tmvaVariables_.push_back("nCharged");
0073   tmvaVariables_.push_back("majW");
0074   tmvaVariables_.push_back("minW");
0075   tmvaVariables_.push_back("frac01");
0076   tmvaVariables_.push_back("frac02");
0077   tmvaVariables_.push_back("frac03");
0078   tmvaVariables_.push_back("frac04");
0079   tmvaVariables_.push_back("ptD");
0080   tmvaVariables_.push_back("beta");
0081   tmvaVariables_.push_back("betaStar");
0082   tmvaVariables_.push_back("dR2Mean");
0083   tmvaVariables_.push_back("pull");
0084   tmvaVariables_.push_back("jetR");
0085   tmvaVariables_.push_back("jetRchg");
0086 
0087   tmvaNames_["rho"] = "rho";
0088   tmvaNames_["nParticles"] = "nParticles";
0089   tmvaNames_["nCharged"] = "nCharged";
0090   tmvaNames_["majW"] = "majW";
0091   tmvaNames_["minW"] = "minW";
0092   tmvaNames_["frac01"] = "frac01";
0093   tmvaNames_["frac02"] = "frac02";
0094   tmvaNames_["frac03"] = "frac03";
0095   tmvaNames_["frac04"] = "frac04";
0096   tmvaNames_["ptD"] = "ptD";
0097   tmvaNames_["beta"] = "beta";
0098   tmvaNames_["betaStar"] = "betaStar";
0099   tmvaNames_["dR2Mean"] = "dR2Mean";
0100   tmvaNames_["pull"] = "pull";
0101   tmvaNames_["jetR"] = "jetR";
0102   tmvaNames_["jetRchg"] = "jetRchg";
0103 }
0104 
0105 MVAJetPuId::~MVAJetPuId() {
0106   if (reader_) {
0107     delete reader_;
0108   }
0109 }
0110 
0111 void Assign(const std::vector<float> &vec, float &a, float &b, float &c, float &d) {
0112   size_t sz = vec.size();
0113   a = (sz > 0 ? vec[0] : 0.);
0114   b = (sz > 1 ? vec[1] : 0.);
0115   c = (sz > 2 ? vec[2] : 0.);
0116   d = (sz > 3 ? vec[3] : 0.);
0117 }
0118 
0119 void SetPtEtaPhi(const reco::Candidate &p, float &pt, float &eta, float &phi) {
0120   pt = p.pt();
0121   eta = p.eta();
0122   phi = p.phi();
0123 }
0124 
0125 void MVAJetPuId::bookReader() {
0126   reader_ = new TMVA::Reader("!Color:Silent");
0127   assert(!tmvaMethod_.empty() && !tmvaWeights_.empty());
0128   for (std::vector<std::string>::iterator it = tmvaVariables_.begin(); it != tmvaVariables_.end(); ++it) {
0129     if (tmvaNames_[*it].empty()) {
0130       tmvaNames_[*it] = *it;
0131     }
0132     reader_->AddVariable(*it, variables_[tmvaNames_[*it]].first);
0133   }
0134   for (std::vector<std::string>::iterator it = tmvaSpectators_.begin(); it != tmvaSpectators_.end(); ++it) {
0135     if (tmvaNames_[*it].empty()) {
0136       tmvaNames_[*it] = *it;
0137     }
0138     reader_->AddSpectator(*it, variables_[tmvaNames_[*it]].first);
0139   }
0140   reco::details::loadTMVAWeights(reader_, tmvaMethod_, tmvaWeights_);
0141 }
0142 
0143 void MVAJetPuId::set(const PileupJetIdentifier &id) { internalId_ = id; }
0144 
0145 void MVAJetPuId::runMva() {
0146   if (!reader_) {
0147     bookReader();
0148   }
0149   if (fabs(internalId_.jetEta_) < 5.0)
0150     internalId_.mva_ = reader_->EvaluateMVA(tmvaMethod_.c_str());
0151   if (fabs(internalId_.jetEta_) >= 5.0)
0152     internalId_.mva_ = -2.;
0153   internalId_.idFlag_ = computeIDflag(internalId_.mva_, internalId_.jetPt_, internalId_.jetEta_);
0154 }
0155 
0156 std::pair<int, int> MVAJetPuId::getJetIdKey(float jetPt, float jetEta) {
0157   int ptId = 0;
0158   if (jetPt > 10 && jetPt < 20)
0159     ptId = 1;
0160   if (jetPt >= 20 && jetPt < 30)
0161     ptId = 2;
0162   if (jetPt >= 30)
0163     ptId = 3;
0164 
0165   int etaId = 0;
0166   if (fabs(jetEta) > 2.5 && fabs(jetEta) < 2.75)
0167     etaId = 1;
0168   if (fabs(jetEta) >= 2.75 && fabs(jetEta) < 3.0)
0169     etaId = 2;
0170   if (fabs(jetEta) >= 3.0 && fabs(jetEta) < 5.0)
0171     etaId = 3;
0172   return std::pair<int, int>(ptId, etaId);
0173 }
0174 
0175 int MVAJetPuId::computeIDflag(float mva, float jetPt, float jetEta) {
0176   std::pair<int, int> jetIdKey = getJetIdKey(jetPt, jetEta);
0177   return computeIDflag(mva, jetIdKey.first, jetIdKey.second);
0178 }
0179 
0180 int MVAJetPuId::computeIDflag(float mva, int ptId, int etaId) {
0181   int idFlag(0);
0182   if (mva > mvacut_[PileupJetIdentifier::kTight][ptId][etaId])
0183     idFlag += 1 << PileupJetIdentifier::kTight;
0184   if (mva > mvacut_[PileupJetIdentifier::kMedium][ptId][etaId])
0185     idFlag += 1 << PileupJetIdentifier::kMedium;
0186   if (mva > mvacut_[PileupJetIdentifier::kLoose][ptId][etaId])
0187     idFlag += 1 << PileupJetIdentifier::kLoose;
0188   return idFlag;
0189 }
0190 
0191 PileupJetIdentifier MVAJetPuId::computeMva() {
0192   runMva();
0193   return PileupJetIdentifier(internalId_);
0194 }
0195 
0196 PileupJetIdentifier MVAJetPuId::computeIdVariables(const reco::Jet *jet,
0197                                                    float jec,
0198                                                    const reco::Vertex *vtx,
0199                                                    const reco::VertexCollection &allvtx,
0200                                                    double rho,
0201                                                    bool calculateMva) {
0202   typedef std::vector<reco::PFCandidatePtr> constituents_type;
0203   typedef std::vector<reco::PFCandidatePtr>::iterator constituents_iterator;
0204 
0205   resetVariables();
0206 
0207   const reco::PFJet *pfjet = dynamic_cast<const reco::PFJet *>(jet);
0208 
0209   if (jec < 0.) {
0210     jec = 1.;
0211   }
0212 
0213   constituents_type constituents = pfjet->getPFConstituents();
0214 
0215   reco::PFCandidatePtr lLead, lSecond, lLeadNeut, lLeadEm, lLeadCh, lTrail;
0216   std::vector<float> frac, fracCh, fracEm, fracNeut;
0217   constexpr int ncones = 4;
0218   std::array<float, ncones> cones{{0.1, 0.2, 0.3, 0.4}};
0219   std::array<float *, ncones> coneFracs{
0220       {&internalId_.frac01_, &internalId_.frac02_, &internalId_.frac03_, &internalId_.frac04_}};
0221   TMatrixDSym covMatrix(2);
0222   covMatrix = 0.;
0223 
0224   reco::TrackRef impactTrack;
0225   float jetPt = jet->pt() / jec;  // use uncorrected pt for shape variables
0226   float sumPt = 0., sumPt2 = 0., sumTkPt = 0., sumPtCh = 0, sumPtNe = 0;
0227   float sum_deta = 0;
0228   float sum_dphi = 0;
0229   float sum_deta2 = 0;
0230   float sum_detadphi = 0;
0231   float sum_dphi2 = 0;
0232   SetPtEtaPhi(
0233       *jet, internalId_.jetPt_, internalId_.jetEta_, internalId_.jetPhi_);  // use corrected pt for jet kinematics
0234   internalId_.jetM_ = jet->mass();
0235   internalId_.rho_ = rho;  //allvtx.size();
0236   for (constituents_iterator it = constituents.begin(); it != constituents.end(); ++it) {
0237     reco::PFCandidatePtr &icand = *it;
0238     float candPt = icand->pt();
0239     float candPtFrac = candPt / jetPt;
0240     float candDr = reco::deltaR(**it, *jet);
0241     float candDeta = fabs((*it)->eta() - jet->eta());
0242     float candDphi = reco::deltaPhi(**it, *jet);
0243     float candPtDr = candPt * candDr;
0244     size_t icone = std::lower_bound(&cones[0], &cones[ncones], candDr) - &cones[0];
0245     float weight2 = candPt * candPt;
0246 
0247     if (lLead.isNull() || candPt > lLead->pt()) {
0248       lSecond = lLead;
0249       lLead = icand;
0250     } else if ((lSecond.isNull() || candPt > lSecond->pt()) && (candPt < lLead->pt())) {
0251       lSecond = icand;
0252     }
0253 
0254     //internalId_.dRMean_     += candPtDr;
0255     internalId_.dR2Mean_ += candPtDr * candPtDr;
0256 
0257     internalId_.ptD_ += candPt * candPt;
0258     sumPt += candPt;
0259     sumPt2 += candPt * candPt;
0260     sum_deta += candDeta * weight2;
0261     sum_dphi += candDphi * weight2;
0262     sum_deta2 += candDeta * candDeta * weight2;
0263     sum_detadphi += candDeta * candDphi * weight2;
0264     sum_dphi2 += candDphi * candDphi * weight2;
0265     //Teta         += candPt * candDR * candDeta;
0266     //Tphi         += candPt * candDR * candDphi;
0267 
0268     frac.push_back(candPtFrac);
0269     if (icone < ncones) {
0270       *coneFracs[icone] += candPt;
0271     }
0272 
0273     if (icand->particleId() == reco::PFCandidate::h0) {
0274       if (lLeadNeut.isNull() || candPt > lLeadNeut->pt()) {
0275         lLeadNeut = icand;
0276       }
0277       internalId_.dRMeanNeut_ += candPtDr;
0278       fracNeut.push_back(candPtFrac);
0279       sumPtNe += candPt;
0280     }
0281 
0282     if (icand->particleId() == reco::PFCandidate::gamma) {
0283       if (lLeadEm.isNull() || candPt > lLeadEm->pt()) {
0284         lLeadEm = icand;
0285       }
0286       internalId_.dRMeanEm_ += candPtDr;
0287       fracEm.push_back(candPtFrac);
0288       sumPtNe += candPt;
0289     }
0290 
0291     if (icand->trackRef().isNonnull() && icand->trackRef().isAvailable()) {
0292       if (lLeadCh.isNull() || candPt > lLeadCh->pt()) {
0293         lLeadCh = icand;
0294       }
0295       //internalId_.jetRchg_  += candPtDr;
0296       fracCh.push_back(candPtFrac);
0297       sumPtCh += candPt;
0298     }
0299 
0300     if (icand->trackRef().isNonnull() && icand->trackRef().isAvailable()) {
0301       float tkpt = icand->trackRef()->pt();
0302       sumTkPt += tkpt;
0303       bool inVtx0 =
0304           find(vtx->tracks_begin(), vtx->tracks_end(), reco::TrackBaseRef(icand->trackRef())) != vtx->tracks_end();
0305       bool inAnyOther = false;
0306 
0307       double dZ0 = fabs(icand->trackRef()->dz(vtx->position()));
0308       double dZ = dZ0;
0309       for (reco::VertexCollection::const_iterator vi = allvtx.begin(); vi != allvtx.end(); ++vi) {
0310         const reco::Vertex &iv = *vi;
0311         if (iv.isFake() || iv.ndof() < 4) {
0312           continue;
0313         }
0314 
0315         bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
0316 
0317         if (!isVtx0 && !inAnyOther) {
0318           inAnyOther =
0319               find(iv.tracks_begin(), iv.tracks_end(), reco::TrackBaseRef(icand->trackRef())) != iv.tracks_end();
0320         }
0321 
0322         dZ = std::min(dZ, fabs(icand->trackRef()->dz(iv.position())));
0323       }
0324       if (inVtx0 && !inAnyOther) {
0325         internalId_.betaClassic_ += tkpt;
0326       } else if (!inVtx0 && inAnyOther) {
0327         internalId_.betaStarClassic_ += tkpt;
0328       }
0329 
0330       if (dZ0 < 0.2) {
0331         internalId_.beta_ += tkpt;
0332       } else if (dZ < 0.2) {
0333         internalId_.betaStar_ += tkpt;
0334       }
0335     }
0336 
0337     if (lTrail.isNull() || candPt < lTrail->pt()) {
0338       lTrail = icand;
0339     }
0340   }
0341 
0342   assert(lLead.isNonnull());
0343   if (lSecond.isNull()) {
0344     lSecond = lTrail;
0345   }
0346   if (lLeadNeut.isNull()) {
0347     lLeadNeut = lTrail;
0348   }
0349   if (lLeadEm.isNull()) {
0350     lLeadEm = lTrail;
0351   }
0352   if (lLeadCh.isNull()) {
0353     lLeadCh = lTrail;
0354   }
0355   impactTrack = lLeadCh->trackRef();
0356 
0357   internalId_.nCharged_ = pfjet->chargedMultiplicity();
0358   internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
0359   internalId_.chgEMfrac_ = pfjet->chargedEmEnergy() / jet->energy();
0360   internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() / jet->energy();
0361   internalId_.chgHadrfrac_ = pfjet->chargedHadronEnergy() / jet->energy();
0362   internalId_.neuHadrfrac_ = pfjet->neutralHadronEnergy() / jet->energy();
0363 
0364   if (impactTrack.isNonnull() && impactTrack.isAvailable()) {
0365     internalId_.d0_ = fabs(impactTrack->dxy(vtx->position()));
0366     internalId_.dZ_ = fabs(impactTrack->dz(vtx->position()));
0367   } else {
0368     internalId_.nParticles_ = constituents.size();
0369     SetPtEtaPhi(*lLead, internalId_.leadPt_, internalId_.leadEta_, internalId_.leadPhi_);
0370     SetPtEtaPhi(*lSecond, internalId_.secondPt_, internalId_.secondEta_, internalId_.secondPhi_);
0371     SetPtEtaPhi(*lLeadNeut, internalId_.leadNeutPt_, internalId_.leadNeutEta_, internalId_.leadNeutPhi_);
0372     SetPtEtaPhi(*lLeadEm, internalId_.leadEmPt_, internalId_.leadEmEta_, internalId_.leadEmPhi_);
0373     SetPtEtaPhi(*lLeadCh, internalId_.leadChPt_, internalId_.leadChEta_, internalId_.leadChPhi_);
0374     std::sort(frac.begin(), frac.end(), std::greater<float>());
0375     std::sort(fracCh.begin(), fracCh.end(), std::greater<float>());
0376     std::sort(fracEm.begin(), fracEm.end(), std::greater<float>());
0377     std::sort(fracNeut.begin(), fracNeut.end(), std::greater<float>());
0378     Assign(frac, internalId_.leadFrac_, internalId_.secondFrac_, internalId_.thirdFrac_, internalId_.fourthFrac_);
0379 
0380     //covMatrix(0,0) /= sumPt2;
0381     //covMatrix(0,1) /= sumPt2;
0382     //covMatrix(1,1) /= sumPt2;
0383     //covMatrix(1,0)  = covMatrix(0,1);
0384     //internalId_.etaW_ = sqrt(covMatrix(0,0));
0385     //internalId_.phiW_ = sqrt(covMatrix(1,1));
0386     //internalId_.jetW_ = 0.5*(internalId_.etaW_+internalId_.phiW_);
0387     //TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
0388     //
0389     if (internalId_.majW_ < internalId_.minW_) {
0390       std::swap(internalId_.majW_, internalId_.minW_);
0391     }
0392 
0393     //internalId_.dRLeadCent_ = reco::deltaR(*jet,*lLead);
0394     if (lSecond.isNonnull()) {
0395       internalId_.dRLead2nd_ = reco::deltaR(*jet, *lSecond);
0396     }
0397     internalId_.dRMeanNeut_ /= jetPt;
0398     internalId_.dRMeanEm_ /= jetPt;
0399     //internalId_.jetRchg_   /= jetPt;
0400     internalId_.dR2Mean_ /= sumPt2;
0401     for (size_t ic = 0; ic < ncones; ++ic) {
0402       *coneFracs[ic] /= jetPt;
0403     }
0404 
0405     double ptMean = sumPt / internalId_.nParticles_;
0406     double ptRMS = 0;
0407     for (unsigned int i0 = 0; i0 < frac.size(); i0++) {
0408       ptRMS += (frac[i0] - ptMean) * (frac[i0] - ptMean);
0409     }
0410     ptRMS /= internalId_.nParticles_;
0411     ptRMS = sqrt(ptRMS);
0412     internalId_.jetRchg_ = internalId_.leadChPt_ / sumPt;
0413     internalId_.jetR_ = internalId_.leadPt_ / sumPt;
0414 
0415     internalId_.ptMean_ = ptMean;
0416     internalId_.ptRMS_ = ptRMS / jetPt;
0417     internalId_.pt2A_ = sqrt(internalId_.ptD_ / internalId_.nParticles_) / jetPt;
0418     internalId_.ptD_ = sqrt(internalId_.ptD_) / sumPt;
0419     internalId_.sumPt_ = sumPt;
0420     internalId_.sumChPt_ = sumPtCh;
0421     internalId_.sumNePt_ = sumPtNe;
0422     if (sumPt > 0) {
0423       internalId_.beta_ /= sumPt;
0424       internalId_.betaStar_ /= sumPt;
0425     } else {
0426       assert(internalId_.beta_ == 0. && internalId_.betaStar_ == 0.);
0427     }
0428     float ave_deta = sum_deta / sumPt2;
0429     float ave_dphi = sum_dphi / sumPt2;
0430     float ave_deta2 = sum_deta2 / sumPt2;
0431     float ave_dphi2 = sum_dphi2 / sumPt2;
0432     float a = ave_deta2 - ave_deta * ave_deta;
0433     float b = ave_dphi2 - ave_dphi * ave_dphi;
0434     float c = -(sum_detadphi / sumPt2 - ave_deta * ave_dphi);
0435     float axis1 = 0;
0436     float axis2 = 0;
0437     if ((((a - b) * (a - b) + 4 * c * c)) > 0) {
0438       float delta = sqrt(((a - b) * (a - b) + 4 * c * c));
0439       if (a + b + delta > 0) {
0440         axis1 = sqrt(0.5 * (a + b + delta));
0441       }
0442       if (a + b - delta > 0) {
0443         axis2 = sqrt(0.5 * (a + b - delta));
0444       }
0445     } else {
0446       axis1 = -1;
0447       axis2 = -1;
0448     }
0449     internalId_.majW_ = axis1;  //sqrt(fabs(eigVals(0)));
0450     internalId_.minW_ = axis2;  //sqrt(fabs(eigVals(1)));
0451     //compute Pull
0452 
0453     float ddetaR_sum(0.0), ddphiR_sum(0.0);
0454     for (int i = 0; i < internalId_.nParticles_; ++i) {
0455       reco::PFCandidatePtr part = pfjet->getPFConstituent(i);
0456       float weight = part->pt() * part->pt();
0457       float deta = part->eta() - jet->eta();
0458       float dphi = reco::deltaPhi(*part, *jet);
0459       float ddeta, ddphi, ddR;
0460       ddeta = deta - ave_deta;
0461       ddphi = reco::deltaPhi(dphi, ave_dphi);  //2*atan(tan((dphi - ave_dphi)/2.)) ;
0462       ddR = sqrt(ddeta * ddeta + ddphi * ddphi);
0463       ddetaR_sum += ddR * ddeta * weight;
0464       ddphiR_sum += ddR * ddphi * weight;
0465     }
0466     if (sumPt2 > 0) {
0467       float ddetaR_ave = ddetaR_sum / sumPt2;
0468       float ddphiR_ave = ddphiR_sum / sumPt2;
0469       internalId_.dRMean_ = sqrt(ddetaR_ave * ddetaR_ave + ddphiR_ave * ddphiR_ave);
0470     }
0471   }
0472 
0473   if (calculateMva) {
0474     runMva();
0475   }
0476 
0477   return PileupJetIdentifier(internalId_);
0478 }
0479 
0480 std::string MVAJetPuId::dumpVariables() const {
0481   std::stringstream out;
0482   for (variables_list_t::const_iterator it = variables_.begin(); it != variables_.end(); ++it) {
0483     out << std::setw(15) << it->first << std::setw(3) << "=" << std::setw(5) << *it->second.first << " ("
0484         << std::setw(5) << it->second.second << ")" << std::endl;
0485   }
0486   return out.str();
0487 }
0488 
0489 void MVAJetPuId::resetVariables() {
0490   internalId_.idFlag_ = 0;
0491   for (variables_list_t::iterator it = variables_.begin(); it != variables_.end(); ++it) {
0492     *it->second.first = it->second.second;
0493   }
0494 }
0495 
0496 #define INIT_VARIABLE(NAME, TMVANAME, VAL) \
0497   internalId_.NAME##_ = VAL;               \
0498   variables_[#NAME] = std::make_pair(&internalId_.NAME##_, VAL);
0499 
0500 void MVAJetPuId::initVariables() {
0501   internalId_.idFlag_ = 0;
0502   INIT_VARIABLE(mva, "", -100.);
0503 
0504   INIT_VARIABLE(jetPt, "jetPt", 0.);
0505   INIT_VARIABLE(jetEta, "jetEta", large_val);
0506   INIT_VARIABLE(jetPhi, "", large_val);
0507   INIT_VARIABLE(jetM, "", 0.);
0508   INIT_VARIABLE(nCharged, "nCharged", 0.);
0509   INIT_VARIABLE(nNeutrals, "", 0.);
0510 
0511   INIT_VARIABLE(chgEMfrac, "", 0.);
0512   INIT_VARIABLE(neuEMfrac, "", 0.);
0513   INIT_VARIABLE(chgHadrfrac, "", 0.);
0514   INIT_VARIABLE(neuHadrfrac, "", 0.);
0515 
0516   INIT_VARIABLE(d0, "", -1000.);
0517   INIT_VARIABLE(dZ, "", -1000.);
0518   INIT_VARIABLE(nParticles, "nParticles", 0.);
0519 
0520   INIT_VARIABLE(leadPt, "", 0.);
0521   INIT_VARIABLE(leadEta, "", large_val);
0522   INIT_VARIABLE(leadPhi, "", large_val);
0523   INIT_VARIABLE(secondPt, "", 0.);
0524   INIT_VARIABLE(secondEta, "", large_val);
0525   INIT_VARIABLE(secondPhi, "", large_val);
0526   INIT_VARIABLE(leadNeutPt, "", 0.);
0527   INIT_VARIABLE(leadNeutEta, "", large_val);
0528 
0529   INIT_VARIABLE(jetR, "jetR", 0.);
0530   INIT_VARIABLE(pull, "pull", 0.);
0531   INIT_VARIABLE(jetRchg, "jetRchg", 0.);
0532   INIT_VARIABLE(dR2Mean, "dR2Mean", 0.);
0533 
0534   INIT_VARIABLE(ptD, "ptD", 0.);
0535   INIT_VARIABLE(ptMean, "", 0.);
0536   INIT_VARIABLE(ptRMS, "", 0.);
0537   INIT_VARIABLE(pt2A, "", 0.);
0538   INIT_VARIABLE(ptDCh, "", 0.);
0539   INIT_VARIABLE(ptDNe, "", 0.);
0540   INIT_VARIABLE(sumPt, "", 0.);
0541   INIT_VARIABLE(sumChPt, "", 0.);
0542   INIT_VARIABLE(sumNePt, "", 0.);
0543   INIT_VARIABLE(secondFrac, "", 0.);
0544   INIT_VARIABLE(thirdFrac, "", 0.);
0545   INIT_VARIABLE(fourthFrac, "", 0.);
0546   INIT_VARIABLE(leadChFrac, "", 0.);
0547   INIT_VARIABLE(secondChFrac, "", 0.);
0548   INIT_VARIABLE(thirdChFrac, "", 0.);
0549   INIT_VARIABLE(fourthChFrac, "", 0.);
0550   INIT_VARIABLE(leadNeutFrac, "", 0.);
0551   INIT_VARIABLE(secondNeutFrac, "", 0.);
0552   INIT_VARIABLE(thirdNeutFrac, "", 0.);
0553   INIT_VARIABLE(fourthNeutFrac, "", 0.);
0554   INIT_VARIABLE(leadEmFrac, "", 0.);
0555   INIT_VARIABLE(secondEmFrac, "", 0.);
0556   INIT_VARIABLE(thirdEmFrac, "", 0.);
0557   INIT_VARIABLE(fourthEmFrac, "", 0.);
0558   INIT_VARIABLE(jetW, "", 1.);
0559   INIT_VARIABLE(etaW, "", 1.);
0560   INIT_VARIABLE(phiW, "", 1.);
0561   INIT_VARIABLE(majW, "majW", 1.);
0562   INIT_VARIABLE(minW, "minW", 1.);
0563   INIT_VARIABLE(frac01, "frac01", 0.);
0564   INIT_VARIABLE(frac02, "frac02", 0.);
0565   INIT_VARIABLE(frac03, "frac03", 0.);
0566   INIT_VARIABLE(frac04, "frac04", 0.);
0567 
0568   INIT_VARIABLE(beta, "beta", 0.);
0569   INIT_VARIABLE(betaStar, "betaStar", 0.);
0570   INIT_VARIABLE(betaClassic, "betaClassic", 0.);
0571   INIT_VARIABLE(betaStarClassic, "betaStarClassic", 0.);
0572   INIT_VARIABLE(rho, "rho", 0.);
0573 }
0574 #undef INIT_VARIABLE