00001 #include "RecoJets/JetProducers/interface/PileupJetIdAlgo.h"
00002
00003 #include "DataFormats/JetReco/interface/PFJet.h"
00004 #include "DataFormats/JetReco/interface/Jet.h"
00005 #include "DataFormats/VertexReco/interface/Vertex.h"
00006 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
00007 #include "DataFormats/Math/interface/deltaR.h"
00008 #include "FWCore/ParameterSet/interface/FileInPath.h"
00009 #include "TMatrixDSym.h"
00010 #include "TMatrixDSymEigen.h"
00011
00012
00013 const float large_val = std::numeric_limits<float>::max();
00014
00015
00016 PileupJetIdAlgo::PileupJetIdAlgo(const edm::ParameterSet & ps)
00017 {
00018 impactParTkThreshod_ = 1.;
00019 cutBased_ = false;
00020
00021 cutBased_ = ps.getParameter<bool>("cutBased");
00022 if(!cutBased_)
00023 {
00024 tmvaWeights_ = edm::FileInPath(ps.getParameter<std::string>("tmvaWeights")).fullPath();
00025 tmvaMethod_ = ps.getParameter<std::string>("tmvaMethod");
00026 tmvaVariables_ = ps.getParameter<std::vector<std::string> >("tmvaVariables");
00027 tmvaSpectators_ = ps.getParameter<std::vector<std::string> >("tmvaSpectators");
00028 version_ = ps.getParameter<int>("version");
00029 }
00030 reader_ = 0;
00031 edm::ParameterSet jetConfig = ps.getParameter<edm::ParameterSet>("JetIdParams");
00032 for(int i0 = 0; i0 < 3; i0++) {
00033 std::string lCutType = "Tight";
00034 if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
00035 if(i0 == PileupJetIdentifier::kLoose) lCutType = "Loose";
00036 int nCut = 1;
00037 if(cutBased_) nCut++;
00038 for(int i1 = 0; i1 < nCut; i1++) {
00039 std::string lFullCutType = lCutType;
00040 if(cutBased_ && i1 == 0) lFullCutType = "BetaStar"+ lCutType;
00041 if(cutBased_ && i1 == 1) lFullCutType = "RMS" + lCutType;
00042 std::vector<double> pt010 = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
00043 std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
00044 std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
00045 std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
00046 if(!cutBased_) {
00047 for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][0][i2] = pt010 [i2];
00048 for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][1][i2] = pt1020[i2];
00049 for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][2][i2] = pt2030[i2];
00050 for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][3][i2] = pt3050[i2];
00051 }
00052 if(cutBased_ && i1 == 0) {
00053 for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][0][i2] = pt010 [i2];
00054 for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][1][i2] = pt1020[i2];
00055 for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][2][i2] = pt2030[i2];
00056 for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][3][i2] = pt3050[i2];
00057 }
00058 if(cutBased_ && i1 == 1) {
00059 for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][0][i2] = pt010 [i2];
00060 for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][1][i2] = pt1020[i2];
00061 for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][2][i2] = pt2030[i2];
00062 for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][3][i2] = pt3050[i2];
00063 }
00064 }
00065 }
00066 setup();
00067 }
00068
00069
00070
00071 PileupJetIdAlgo::PileupJetIdAlgo(int version,
00072 const std::string & tmvaWeights,
00073 const std::string & tmvaMethod,
00074 Float_t impactParTkThreshod,
00075 const std::vector<std::string> & tmvaVariables
00076 )
00077 {
00078 impactParTkThreshod_ = impactParTkThreshod;
00079 tmvaWeights_ = tmvaWeights;
00080 tmvaMethod_ = tmvaMethod;
00081 tmvaVariables_ = tmvaVariables;
00082 version_ = version;
00083
00084 reader_ = 0;
00085
00086 setup();
00087 }
00088
00089
00090 void PileupJetIdAlgo::setup()
00091 {
00092 initVariables();
00093
00094 if( version_ == PHILv0 ) {
00095 tmvaVariables_.clear();
00096 tmvaVariables_.push_back( "jspt_1" );
00097 tmvaVariables_.push_back( "jseta_1" );
00098 tmvaVariables_.push_back( "jsphi_1" );
00099 tmvaVariables_.push_back( "jd0_1" );
00100 tmvaVariables_.push_back( "jdZ_1" );
00101 tmvaVariables_.push_back( "jm_1" );
00102 tmvaVariables_.push_back( "npart_1" );
00103 tmvaVariables_.push_back( "lpt_1" );
00104 tmvaVariables_.push_back( "leta_1" );
00105 tmvaVariables_.push_back( "lphi_1" );
00106 tmvaVariables_.push_back( "spt_1" );
00107 tmvaVariables_.push_back( "seta_1" );
00108 tmvaVariables_.push_back( "sphi_1" );
00109 tmvaVariables_.push_back( "lnept_1" );
00110 tmvaVariables_.push_back( "lneeta_1");
00111 tmvaVariables_.push_back( "lnephi_1");
00112 tmvaVariables_.push_back( "lempt_1" );
00113 tmvaVariables_.push_back( "lemeta_1");
00114 tmvaVariables_.push_back( "lemphi_1");
00115 tmvaVariables_.push_back( "lchpt_1" );
00116
00117 tmvaVariables_.push_back( "lchphi_1");
00118 tmvaVariables_.push_back( "lLfr_1" );
00119 tmvaVariables_.push_back( "drlc_1" );
00120 tmvaVariables_.push_back( "drls_1" );
00121 tmvaVariables_.push_back( "drm_1" );
00122 tmvaVariables_.push_back( "drmne_1" );
00123 tmvaVariables_.push_back( "drem_1" );
00124 tmvaVariables_.push_back( "drch_1" );
00125
00126 tmvaNames_["jspt_1"] = "jetPt";
00127 tmvaNames_["jseta_1"] = "jetEta";
00128 tmvaNames_["jsphi_1"] = "jetPhi";
00129 tmvaNames_["jm_1"] = "jetM";
00130 tmvaNames_["jd0_1"] = "d0";
00131 tmvaNames_["jdZ_1"] = "dZ";
00132 tmvaNames_["npart_1"] = "nParticles";
00133
00134 tmvaNames_["lpt_1"] = "leadPt";
00135 tmvaNames_["leta_1"] = "leadEta";
00136 tmvaNames_["lphi_1"] = "leadPhi";
00137 tmvaNames_["spt_1"] = "secondPt";
00138 tmvaNames_["seta_1"] = "secondEta";
00139 tmvaNames_["sphi_1"] = "secondPhi";
00140 tmvaNames_["lnept_1"] = "leadNeutPt";
00141 tmvaNames_["lneeta_1"] = "leadNeutEta";
00142 tmvaNames_["lnephi_1"] = "leadNeutPhi";
00143 tmvaNames_["lempt_1"] = "leadEmPt";
00144 tmvaNames_["lemeta_1"] = "leadEmEta";
00145 tmvaNames_["lemphi_1"] = "leadEmPhi";
00146 tmvaNames_["lchpt_1"] = "leadChPt";
00147 tmvaNames_["lcheta_1"] = "leadChEta";
00148 tmvaNames_["lchphi_1"] = "leadChPhi";
00149 tmvaNames_["lLfr_1"] = "leadFrac";
00150
00151 tmvaNames_["drlc_1"] = "dRLeadCent";
00152 tmvaNames_["drls_1"] = "dRLead2nd";
00153 tmvaNames_["drm_1"] = "dRMean";
00154 tmvaNames_["drmne_1"] = "dRMeanNeut";
00155 tmvaNames_["drem_1"] = "dRMeanEm";
00156 tmvaNames_["drch_1"] = "dRMeanCh";
00157
00158 } else if( ! cutBased_ ){
00159 assert( tmvaMethod_.empty() || (! tmvaVariables_.empty() && version_ == USER) );
00160 }
00161 }
00162
00163
00164 PileupJetIdAlgo::~PileupJetIdAlgo()
00165 {
00166 if( reader_ ) {
00167 delete reader_;
00168 }
00169 }
00170
00171
00172 void assign(const std::vector<float> & vec, float & a, float & b, float & c, float & d )
00173 {
00174 size_t sz = vec.size();
00175 a = ( sz > 0 ? vec[0] : 0. );
00176 b = ( sz > 1 ? vec[1] : 0. );
00177 c = ( sz > 2 ? vec[2] : 0. );
00178 d = ( sz > 3 ? vec[3] : 0. );
00179 }
00180
00181 void setPtEtaPhi(const reco::Candidate & p, float & pt, float & eta, float &phi )
00182 {
00183 pt = p.pt();
00184 eta = p.eta();
00185 phi = p.phi();
00186 }
00187
00188
00189 void PileupJetIdAlgo::bookReader()
00190 {
00191 reader_ = new TMVA::Reader("!Color:!Silent");
00192 assert( ! tmvaMethod_.empty() && ! tmvaWeights_.empty() );
00193 for(std::vector<std::string>::iterator it=tmvaVariables_.begin(); it!=tmvaVariables_.end(); ++it) {
00194 if( tmvaNames_[*it].empty() ) {
00195 tmvaNames_[*it] = *it;
00196 }
00197 reader_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
00198 }
00199 for(std::vector<std::string>::iterator it=tmvaSpectators_.begin(); it!=tmvaSpectators_.end(); ++it) {
00200 if( tmvaNames_[*it].empty() ) {
00201 tmvaNames_[*it] = *it;
00202 }
00203 reader_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
00204 }
00205 reader_->BookMVA(tmvaMethod_.c_str(), tmvaWeights_.c_str());
00206 }
00207
00208
00209 void PileupJetIdAlgo::set(const PileupJetIdentifier & id)
00210 {
00211 internalId_ = id;
00212 }
00213
00214
00215 void PileupJetIdAlgo::runMva()
00216 {
00217 if( cutBased_ ) {
00218 internalId_.idFlag_ = computeCutIDflag(internalId_.betaStarClassic_,internalId_.dR2Mean_,internalId_.nvtx_,internalId_.jetPt_,internalId_.jetEta_);
00219 } else {
00220 if( ! reader_ ) { bookReader(); std::cerr << "Reader booked" << std::endl; }
00221 if(fabs(internalId_.jetEta_) < 5.0) internalId_.mva_ = reader_->EvaluateMVA( tmvaMethod_.c_str() );
00222 if(fabs(internalId_.jetEta_) >= 5.0) internalId_.mva_ = -2.;
00223 internalId_.idFlag_ = computeIDflag(internalId_.mva_,internalId_.jetPt_,internalId_.jetEta_);
00224 }
00225 }
00226
00227
00228 std::pair<int,int> PileupJetIdAlgo::getJetIdKey(float jetPt, float jetEta)
00229 {
00230 int ptId = 0;
00231 if(jetPt > 10 && jetPt < 20) ptId = 1;
00232 if(jetPt > 20 && jetPt < 30) ptId = 2;
00233 if(jetPt > 30 ) ptId = 3;
00234
00235 int etaId = 0;
00236 if(fabs(jetEta) > 2.5 && fabs(jetEta) < 2.75) etaId = 1;
00237 if(fabs(jetEta) > 2.75 && fabs(jetEta) < 3.0 ) etaId = 2;
00238 if(fabs(jetEta) > 3.0 && fabs(jetEta) < 5.0 ) etaId = 3;
00239
00240 return std::pair<int,int>(ptId,etaId);
00241 }
00242
00243 int PileupJetIdAlgo::computeCutIDflag(float betaStarClassic,float dR2Mean,float nvtx, float jetPt, float jetEta)
00244 {
00245 std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
00246 float betaStarModified = betaStarClassic/log(nvtx-0.64);
00247 int idFlag(0);
00248 if(betaStarModified < betaStarCut_[PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second] &&
00249 dR2Mean < rmsCut_ [PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second]
00250 ) idFlag += 1 << PileupJetIdentifier::kTight;
00251
00252 if(betaStarModified < betaStarCut_[PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second] &&
00253 dR2Mean < rmsCut_ [PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second]
00254 ) idFlag += 1 << PileupJetIdentifier::kMedium;
00255
00256 if(betaStarModified < betaStarCut_[PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second] &&
00257 dR2Mean < rmsCut_ [PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second]
00258 ) idFlag += 1 << PileupJetIdentifier::kLoose;
00259 return idFlag;
00260 }
00261
00262 int PileupJetIdAlgo::computeIDflag(float mva, float jetPt, float jetEta)
00263 {
00264 std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
00265 return computeIDflag(mva,jetIdKey.first,jetIdKey.second);
00266 }
00267
00268
00269 int PileupJetIdAlgo::computeIDflag(float mva,int ptId,int etaId)
00270 {
00271 int idFlag(0);
00272 if(mva > mvacut_[PileupJetIdentifier::kTight ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kTight;
00273 if(mva > mvacut_[PileupJetIdentifier::kMedium][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kMedium;
00274 if(mva > mvacut_[PileupJetIdentifier::kLoose ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kLoose;
00275 return idFlag;
00276 }
00277
00278
00279
00280 PileupJetIdentifier PileupJetIdAlgo::computeMva()
00281 {
00282 runMva();
00283 return PileupJetIdentifier(internalId_);
00284 }
00285
00286
00287 PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, float jec, const reco::Vertex * vtx,
00288 const reco::VertexCollection & allvtx,
00289 bool calculateMva)
00290 {
00291 static int printWarning = 10;
00292 typedef std::vector <reco::PFCandidatePtr> constituents_type;
00293 typedef std::vector <reco::PFCandidatePtr>::iterator constituents_iterator;
00294
00295
00296 resetVariables();
00297
00298
00299
00300 const reco::PFJet * pfjet = dynamic_cast<const reco::PFJet *>(jet);
00301
00302
00303
00304
00305 if( jec < 0. ) {
00306 jec = 1.;
00307 }
00308
00309 constituents_type constituents = pfjet->getPFConstituents();
00310
00311 reco::PFCandidatePtr lLead, lSecond, lLeadNeut, lLeadEm, lLeadCh, lTrail;
00312 std::vector<float> frac, fracCh, fracEm, fracNeut;
00313 float cones[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
00314 size_t ncones = sizeof(cones)/sizeof(float);
00315 float * coneFracs[] = { &internalId_.frac01_, &internalId_.frac02_, &internalId_.frac03_, &internalId_.frac04_,
00316 &internalId_.frac05_, &internalId_.frac06_, &internalId_.frac07_ };
00317 float * coneEmFracs[] = { &internalId_.emFrac01_, &internalId_.emFrac02_, &internalId_.emFrac03_, &internalId_.emFrac04_,
00318 &internalId_.emFrac05_, &internalId_.emFrac06_, &internalId_.emFrac07_ };
00319 float * coneNeutFracs[] = { &internalId_.neutFrac01_, &internalId_.neutFrac02_, &internalId_.neutFrac03_, &internalId_.neutFrac04_,
00320 &internalId_.neutFrac05_, &internalId_.neutFrac06_, &internalId_.neutFrac07_ };
00321 float * coneChFracs[] = { &internalId_.chFrac01_, &internalId_.chFrac02_, &internalId_.chFrac03_, &internalId_.chFrac04_,
00322 &internalId_.chFrac05_, &internalId_.chFrac06_, &internalId_.chFrac07_ };
00323 TMatrixDSym covMatrix(2); covMatrix = 0.;
00324
00325 reco::TrackRef impactTrack;
00326 float jetPt = jet->pt() / jec;
00327 float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
00328 setPtEtaPhi(*jet,internalId_.jetPt_,internalId_.jetEta_,internalId_.jetPhi_);
00329 internalId_.jetM_ = jet->mass();
00330 internalId_.nvtx_ = allvtx.size();
00331
00332 for(constituents_iterator it=constituents.begin(); it!=constituents.end(); ++it) {
00333 reco::PFCandidatePtr & icand = *it;
00334 float candPt = icand->pt();
00335 float candPtFrac = candPt/jetPt;
00336 float candDr = reco::deltaR(**it,*jet);
00337 float candDeta = fabs( (*it)->eta() - jet->eta() );
00338 float candDphi = reco::deltaPhi(**it,*jet);
00339 float candPtDr = candPt * candDr;
00340 size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
00341
00342
00343 if( lLead.isNull() || candPt > lLead->pt() ) {
00344 lSecond = lLead;
00345 lLead = icand;
00346 } else if( (lSecond.isNull() || candPt > lSecond->pt()) && (candPt < lLead->pt()) ) {
00347 lSecond = icand;
00348 }
00349
00350
00351 internalId_.dRMean_ += candPtDr;
00352 internalId_.dR2Mean_ += candPtDr*candPtDr;
00353 covMatrix(0,0) += candPt*candPt*candDeta*candDeta;
00354 covMatrix(0,1) += candPt*candPt*candDeta*candDphi;
00355 covMatrix(1,1) += candPt*candPt*candDphi*candDphi;
00356 internalId_.ptD_ += candPt*candPt;
00357 sumPt += candPt;
00358 sumPt2 += candPt*candPt;
00359
00360
00361 frac.push_back(candPtFrac);
00362 if( icone < ncones ) { *coneFracs[icone] += candPt; }
00363
00364
00365 if( icand->particleId() == reco::PFCandidate::h0 ) {
00366 if (lLeadNeut.isNull() || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
00367 internalId_.dRMeanNeut_ += candPtDr;
00368 fracNeut.push_back(candPtFrac);
00369 if( icone < ncones ) { *coneNeutFracs[icone] += candPt; }
00370 internalId_.ptDNe_ += candPt*candPt;
00371 sumPtNe += candPt;
00372 }
00373
00374 if( icand->particleId() == reco::PFCandidate::gamma ) {
00375 if(lLeadEm.isNull() || candPt > lLeadEm->pt()) { lLeadEm = icand; }
00376 internalId_.dRMeanEm_ += candPtDr;
00377 fracEm.push_back(candPtFrac);
00378 if( icone < ncones ) { *coneEmFracs[icone] += candPt; }
00379 internalId_.ptDNe_ += candPt*candPt;
00380 sumPtNe += candPt;
00381 }
00382
00383 if( icand->trackRef().isNonnull() && icand->trackRef().isAvailable() ) {
00384 if (lLeadCh.isNull() || candPt > lLeadCh->pt()) { lLeadCh = icand; }
00385 internalId_.dRMeanCh_ += candPtDr;
00386 internalId_.ptDCh_ += candPt*candPt;
00387 fracCh.push_back(candPtFrac);
00388 if( icone < ncones ) { *coneChFracs[icone] += candPt; }
00389 sumPtCh += candPt;
00390 }
00391
00392 if( icand->trackRef().isNonnull() && icand->trackRef().isAvailable() ) {
00393 try {
00394 float tkpt = icand->trackRef()->pt();
00395 sumTkPt += tkpt;
00396
00397 bool inVtx0 = find( vtx->tracks_begin(), vtx->tracks_end(), reco::TrackBaseRef(icand->trackRef()) ) != vtx->tracks_end();
00398 bool inAnyOther = false;
00399
00400 double dZ0 = fabs(icand->trackRef()->dz(vtx->position()));
00401 double dZ = dZ0;
00402 for(reco::VertexCollection::const_iterator vi=allvtx.begin(); vi!=allvtx.end(); ++vi ) {
00403 const reco::Vertex & iv = *vi;
00404 if( iv.isFake() || iv.ndof() < 4 ) { continue; }
00405
00406 bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
00407
00408 if( ! isVtx0 && ! inAnyOther ) {
00409 inAnyOther = find( iv.tracks_begin(), iv.tracks_end(), reco::TrackBaseRef(icand->trackRef()) ) != iv.tracks_end();
00410 }
00411
00412 dZ = std::min(dZ,fabs(icand->trackRef()->dz(iv.position())));
00413 }
00414
00415 if( inVtx0 && ! inAnyOther ) {
00416 internalId_.betaClassic_ += tkpt;
00417 } else if( ! inVtx0 && inAnyOther ) {
00418 internalId_.betaStarClassic_ += tkpt;
00419 }
00420
00421 if( dZ0 < 0.2 ) {
00422 internalId_.beta_ += tkpt;
00423 } else if( dZ < 0.2 ) {
00424 internalId_.betaStar_ += tkpt;
00425 }
00426 } catch (cms::Exception &e) {
00427 if(printWarning-- > 0) { std::cerr << e << std::endl; }
00428 }
00429 }
00430
00431 if( lTrail.isNull() || candPt < lTrail->pt() ) {
00432 lTrail = icand;
00433 }
00434 }
00435
00436
00437 assert( lLead.isNonnull() );
00438
00439 if ( lSecond.isNull() ) { lSecond = lTrail; }
00440 if ( lLeadNeut.isNull() ) { lLeadNeut = lTrail; }
00441 if ( lLeadEm.isNull() ) { lLeadEm = lTrail; }
00442 if ( lLeadCh.isNull() ) { lLeadCh = lTrail; }
00443 impactTrack = lLeadCh->trackRef();
00444
00445 internalId_.nCharged_ = pfjet->chargedMultiplicity();
00446 internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
00447 internalId_.chgEMfrac_ = pfjet->chargedEmEnergy() /jet->energy();
00448 internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() /jet->energy();
00449 internalId_.chgHadrfrac_ = pfjet->chargedHadronEnergy()/jet->energy();
00450 internalId_.neuHadrfrac_ = pfjet->neutralHadronEnergy()/jet->energy();
00451
00452 if( impactTrack.isNonnull() && impactTrack.isAvailable() ) {
00453 internalId_.d0_ = fabs(impactTrack->dxy(vtx->position()));
00454 internalId_.dZ_ = fabs(impactTrack->dz(vtx->position()));
00455 } else {
00456 if(printWarning-- > 0) { std::cerr << "WARNING : did not find any valid track reference attached to the jet " << std::endl; }
00457 }
00458 internalId_.nParticles_ = constituents.size();
00459
00460 setPtEtaPhi(*lLead,internalId_.leadPt_,internalId_.leadEta_,internalId_.leadPhi_);
00461 setPtEtaPhi(*lSecond,internalId_.secondPt_,internalId_.secondEta_,internalId_.secondPhi_);
00462 setPtEtaPhi(*lLeadNeut,internalId_.leadNeutPt_,internalId_.leadNeutEta_,internalId_.leadNeutPhi_);
00463 setPtEtaPhi(*lLeadEm,internalId_.leadEmPt_,internalId_.leadEmEta_,internalId_.leadEmPhi_);
00464 setPtEtaPhi(*lLeadCh,internalId_.leadChPt_,internalId_.leadChEta_,internalId_.leadChPhi_);
00465
00466 std::sort(frac.begin(),frac.end(),std::greater<float>());
00467 std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
00468 std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
00469 std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
00470 assign(frac, internalId_.leadFrac_, internalId_.secondFrac_, internalId_.thirdFrac_, internalId_.fourthFrac_);
00471 assign(fracCh, internalId_.leadChFrac_, internalId_.secondChFrac_, internalId_.thirdChFrac_, internalId_.fourthChFrac_);
00472 assign(fracEm, internalId_.leadEmFrac_, internalId_.secondEmFrac_, internalId_.thirdEmFrac_, internalId_.fourthEmFrac_);
00473 assign(fracNeut,internalId_.leadNeutFrac_,internalId_.secondNeutFrac_,internalId_.thirdNeutFrac_,internalId_.fourthNeutFrac_);
00474
00475 covMatrix(0,0) /= sumPt2;
00476 covMatrix(0,1) /= sumPt2;
00477 covMatrix(1,1) /= sumPt2;
00478 covMatrix(1,0) = covMatrix(0,1);
00479 internalId_.etaW_ = sqrt(covMatrix(0,0));
00480 internalId_.phiW_ = sqrt(covMatrix(1,1));
00481 internalId_.jetW_ = 0.5*(internalId_.etaW_+internalId_.phiW_);
00482 TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
00483 internalId_.majW_ = sqrt(fabs(eigVals(0)));
00484 internalId_.minW_ = sqrt(fabs(eigVals(1)));
00485 if( internalId_.majW_ < internalId_.minW_ ) { std::swap(internalId_.majW_,internalId_.minW_); }
00486
00487 internalId_.dRLeadCent_ = reco::deltaR(*jet,*lLead);
00488 if( lSecond.isNonnull() ) { internalId_.dRLead2nd_ = reco::deltaR(*jet,*lSecond); }
00489 internalId_.dRMean_ /= jetPt;
00490 internalId_.dRMeanNeut_ /= jetPt;
00491 internalId_.dRMeanEm_ /= jetPt;
00492 internalId_.dRMeanCh_ /= jetPt;
00493 internalId_.dR2Mean_ /= sumPt2;
00494
00495 for(size_t ic=0; ic<ncones; ++ic){
00496 *coneFracs[ic] /= jetPt;
00497 *coneEmFracs[ic] /= jetPt;
00498 *coneNeutFracs[ic] /= jetPt;
00499 *coneChFracs[ic] /= jetPt;
00500 }
00501
00502 double ptMean = sumPt/internalId_.nParticles_;
00503 double ptRMS = 0;
00504 for(unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
00505 ptRMS/=internalId_.nParticles_;
00506 ptRMS=sqrt(ptRMS);
00507
00508 internalId_.ptMean_ = ptMean;
00509 internalId_.ptRMS_ = ptRMS/jetPt;
00510 internalId_.pt2A_ = sqrt( internalId_.ptD_ /internalId_.nParticles_)/jetPt;
00511 internalId_.ptD_ = sqrt( internalId_.ptD_) / sumPt;
00512 internalId_.ptDCh_ = sqrt( internalId_.ptDCh_) / sumPtCh;
00513 internalId_.ptDNe_ = sqrt( internalId_.ptDNe_) / sumPtNe;
00514 internalId_.sumPt_ = sumPt;
00515 internalId_.sumChPt_ = sumPtCh;
00516 internalId_.sumNePt_ = sumPtNe;
00517
00518 if( sumTkPt != 0. ) {
00519 internalId_.beta_ /= sumTkPt;
00520 internalId_.betaStar_ /= sumTkPt;
00521 internalId_.betaClassic_ /= sumTkPt;
00522 internalId_.betaStarClassic_ /= sumTkPt;
00523 } else {
00524 assert( internalId_.beta_ == 0. && internalId_.betaStar_ == 0.&& internalId_.betaClassic_ == 0. && internalId_.betaStarClassic_ == 0. );
00525 }
00526
00527 if( calculateMva ) {
00528 runMva();
00529 }
00530
00531 return PileupJetIdentifier(internalId_);
00532 }
00533
00534
00535
00536
00537 std::string PileupJetIdAlgo::dumpVariables() const
00538 {
00539 std::stringstream out;
00540 for(variables_list_t::const_iterator it=variables_.begin();
00541 it!=variables_.end(); ++it ) {
00542 out << std::setw(15) << it->first << std::setw(3) << "="
00543 << std::setw(5) << *it->second.first
00544 << " (" << std::setw(5) << it->second.second << ")" << std::endl;
00545 }
00546 return out.str();
00547 }
00548
00549
00550 void PileupJetIdAlgo::resetVariables()
00551 {
00552 internalId_.idFlag_ = 0;
00553 for(variables_list_t::iterator it=variables_.begin();
00554 it!=variables_.end(); ++it ) {
00555 *it->second.first = it->second.second;
00556 }
00557 }
00558
00559
00560 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
00561 internalId_.NAME ## _ = VAL; \
00562 variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
00563
00564
00565 void PileupJetIdAlgo::initVariables()
00566 {
00567 internalId_.idFlag_ = 0;
00568 INIT_VARIABLE(mva , "", -100.);
00569 INIT_VARIABLE(jetPt , "jspt_1", 0.);
00570 INIT_VARIABLE(jetEta , "jseta_1", large_val);
00571 INIT_VARIABLE(jetPhi , "jsphi_1", large_val);
00572 INIT_VARIABLE(jetM , "jm_1", 0.);
00573
00574 INIT_VARIABLE(nCharged , "", 0.);
00575 INIT_VARIABLE(nNeutrals , "", 0.);
00576
00577 INIT_VARIABLE(chgEMfrac , "", 0.);
00578 INIT_VARIABLE(neuEMfrac , "", 0.);
00579 INIT_VARIABLE(chgHadrfrac, "", 0.);
00580 INIT_VARIABLE(neuHadrfrac, "", 0.);
00581
00582 INIT_VARIABLE(d0 , "jd0_1" , -1000.);
00583 INIT_VARIABLE(dZ , "jdZ_1" , -1000.);
00584 INIT_VARIABLE(nParticles , "npart_1" , 0.);
00585
00586 INIT_VARIABLE(leadPt , "lpt_1" , 0.);
00587 INIT_VARIABLE(leadEta , "leta_1" , large_val);
00588 INIT_VARIABLE(leadPhi , "lphi_1" , large_val);
00589 INIT_VARIABLE(secondPt , "spt_1" , 0.);
00590 INIT_VARIABLE(secondEta , "seta_1" , large_val);
00591 INIT_VARIABLE(secondPhi , "sphi_1" , large_val);
00592 INIT_VARIABLE(leadNeutPt , "lnept_1" , 0.);
00593 INIT_VARIABLE(leadNeutEta, "lneeta_1" , large_val);
00594 INIT_VARIABLE(leadNeutPhi, "lnephi_1" , large_val);
00595 INIT_VARIABLE(leadEmPt , "lempt_1" , 0.);
00596 INIT_VARIABLE(leadEmEta , "lemeta_1" , large_val);
00597 INIT_VARIABLE(leadEmPhi , "lemphi_1" , large_val);
00598 INIT_VARIABLE(leadChPt , "lchpt_1" , 0.);
00599 INIT_VARIABLE(leadChEta , "lcheta_1" , large_val);
00600 INIT_VARIABLE(leadChPhi , "lchphi_1" , large_val);
00601 INIT_VARIABLE(leadFrac , "lLfr_1" , 0.);
00602
00603 INIT_VARIABLE(dRLeadCent , "drlc_1" , 0.);
00604 INIT_VARIABLE(dRLead2nd , "drls_1" , 0.);
00605 INIT_VARIABLE(dRMean , "drm_1" , 0.);
00606 INIT_VARIABLE(dRMeanNeut , "drmne_1" , 0.);
00607 INIT_VARIABLE(dRMeanEm , "drem_1" , 0.);
00608 INIT_VARIABLE(dRMeanCh , "drch_1" , 0.);
00609 INIT_VARIABLE(dR2Mean , "" , 0.);
00610
00611 INIT_VARIABLE(ptD , "", 0.);
00612 INIT_VARIABLE(ptMean , "", 0.);
00613 INIT_VARIABLE(ptRMS , "", 0.);
00614 INIT_VARIABLE(pt2A , "", 0.);
00615 INIT_VARIABLE(ptDCh , "", 0.);
00616 INIT_VARIABLE(ptDNe , "", 0.);
00617 INIT_VARIABLE(sumPt , "", 0.);
00618 INIT_VARIABLE(sumChPt , "", 0.);
00619 INIT_VARIABLE(sumNePt , "", 0.);
00620
00621 INIT_VARIABLE(secondFrac ,"" ,0.);
00622 INIT_VARIABLE(thirdFrac ,"" ,0.);
00623 INIT_VARIABLE(fourthFrac ,"" ,0.);
00624
00625 INIT_VARIABLE(leadChFrac ,"" ,0.);
00626 INIT_VARIABLE(secondChFrac ,"" ,0.);
00627 INIT_VARIABLE(thirdChFrac ,"" ,0.);
00628 INIT_VARIABLE(fourthChFrac ,"" ,0.);
00629
00630 INIT_VARIABLE(leadNeutFrac ,"" ,0.);
00631 INIT_VARIABLE(secondNeutFrac ,"" ,0.);
00632 INIT_VARIABLE(thirdNeutFrac ,"" ,0.);
00633 INIT_VARIABLE(fourthNeutFrac ,"" ,0.);
00634
00635 INIT_VARIABLE(leadEmFrac ,"" ,0.);
00636 INIT_VARIABLE(secondEmFrac ,"" ,0.);
00637 INIT_VARIABLE(thirdEmFrac ,"" ,0.);
00638 INIT_VARIABLE(fourthEmFrac ,"" ,0.);
00639
00640 INIT_VARIABLE(jetW ,"" ,1.);
00641 INIT_VARIABLE(etaW ,"" ,1.);
00642 INIT_VARIABLE(phiW ,"" ,1.);
00643
00644 INIT_VARIABLE(majW ,"" ,1.);
00645 INIT_VARIABLE(minW ,"" ,1.);
00646
00647 INIT_VARIABLE(frac01 ,"" ,0.);
00648 INIT_VARIABLE(frac02 ,"" ,0.);
00649 INIT_VARIABLE(frac03 ,"" ,0.);
00650 INIT_VARIABLE(frac04 ,"" ,0.);
00651 INIT_VARIABLE(frac05 ,"" ,0.);
00652 INIT_VARIABLE(frac06 ,"" ,0.);
00653 INIT_VARIABLE(frac07 ,"" ,0.);
00654
00655 INIT_VARIABLE(chFrac01 ,"" ,0.);
00656 INIT_VARIABLE(chFrac02 ,"" ,0.);
00657 INIT_VARIABLE(chFrac03 ,"" ,0.);
00658 INIT_VARIABLE(chFrac04 ,"" ,0.);
00659 INIT_VARIABLE(chFrac05 ,"" ,0.);
00660 INIT_VARIABLE(chFrac06 ,"" ,0.);
00661 INIT_VARIABLE(chFrac07 ,"" ,0.);
00662
00663 INIT_VARIABLE(neutFrac01 ,"" ,0.);
00664 INIT_VARIABLE(neutFrac02 ,"" ,0.);
00665 INIT_VARIABLE(neutFrac03 ,"" ,0.);
00666 INIT_VARIABLE(neutFrac04 ,"" ,0.);
00667 INIT_VARIABLE(neutFrac05 ,"" ,0.);
00668 INIT_VARIABLE(neutFrac06 ,"" ,0.);
00669 INIT_VARIABLE(neutFrac07 ,"" ,0.);
00670
00671 INIT_VARIABLE(emFrac01 ,"" ,0.);
00672 INIT_VARIABLE(emFrac02 ,"" ,0.);
00673 INIT_VARIABLE(emFrac03 ,"" ,0.);
00674 INIT_VARIABLE(emFrac04 ,"" ,0.);
00675 INIT_VARIABLE(emFrac05 ,"" ,0.);
00676 INIT_VARIABLE(emFrac06 ,"" ,0.);
00677 INIT_VARIABLE(emFrac07 ,"" ,0.);
00678
00679 INIT_VARIABLE(beta ,"" ,0.);
00680 INIT_VARIABLE(betaStar ,"" ,0.);
00681 INIT_VARIABLE(betaClassic ,"" ,0.);
00682 INIT_VARIABLE(betaStarClassic ,"" ,0.);
00683
00684 INIT_VARIABLE(nvtx ,"" ,0.);
00685
00686 }
00687
00688 #undef INIT_VARIABLE
1.7.3