00001 #include "RecoJets/JetPlusTracks/interface/JetPlusTrackCorrector.h"
00002 #include "DataFormats/Common/interface/Handle.h"
00003 #include "DataFormats/Common/interface/View.h"
00004 #include "DataFormats/Math/interface/deltaR.h"
00005 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00006 #include "FWCore/Framework/interface/ESHandle.h"
00007 #include "FWCore/Framework/interface/Event.h"
00008 #include "FWCore/Framework/interface/EventSetup.h"
00009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00011 #include <fstream>
00012 #include <vector>
00013 #include <iomanip>
00014 #include <math.h>
00015
00016 using namespace std;
00017 using namespace jpt;
00018
00019
00020
00021 JetPlusTrackCorrector::JetPlusTrackCorrector( const edm::ParameterSet& pset )
00022 : verbose_( pset.getParameter<bool>("Verbose") ),
00023 vectorial_( pset.getParameter<bool>("VectorialCorrection") ),
00024 vecResponse_( pset.getParameter<bool>("UseResponseInVecCorr") ),
00025 useInConeTracks_( pset.getParameter<bool>("UseInConeTracks") ),
00026 useOutOfConeTracks_( pset.getParameter<bool>("UseOutOfConeTracks") ),
00027 useOutOfVertexTracks_( pset.getParameter<bool>("UseOutOfVertexTracks") ),
00028 usePions_( pset.getParameter<bool>("UsePions") ),
00029 useEff_( pset.getParameter<bool>("UseEfficiency") ),
00030 useMuons_( pset.getParameter<bool>("UseMuons") ),
00031 useElecs_( pset.getParameter<bool>("UseElectrons") ),
00032 useTrackQuality_( pset.getParameter<bool>("UseTrackQuality") ),
00033 jetTracksAtVertex_( pset.getParameter<edm::InputTag>("JetTracksAssociationAtVertex") ),
00034 jetTracksAtCalo_( pset.getParameter<edm::InputTag>("JetTracksAssociationAtCaloFace") ),
00035 jetSplitMerge_( pset.getParameter<int>("JetSplitMerge") ),
00036 srcPVs_(pset.getParameter<edm::InputTag>("srcPVs") ),
00037 ptErrorQuality_( pset.getParameter<double>("PtErrorQuality") ),
00038 dzVertexCut_( pset.getParameter<double>("DzVertexCut") ),
00039 muons_( pset.getParameter<edm::InputTag>("Muons") ),
00040 electrons_( pset.getParameter<edm::InputTag>("Electrons") ),
00041 electronIds_( pset.getParameter<edm::InputTag>("ElectronIds") ),
00042 trackQuality_( reco::TrackBase::qualityByName( pset.getParameter<std::string>("TrackQuality") ) ),
00043 response_( Map( pset.getParameter<std::string>("ResponseMap"), verbose_ ) ),
00044 efficiency_( Map( pset.getParameter<std::string>("EfficiencyMap"), verbose_ ) ),
00045 leakage_( Map( pset.getParameter<std::string>("LeakageMap"), verbose_ ) ),
00046 pionMass_(0.140),
00047 muonMass_(0.105),
00048 elecMass_(0.000511),
00049 maxEta_( pset.getParameter<double>("MaxJetEta") )
00050 {
00051
00052 if ( verbose_ ) {
00053 std::stringstream ss;
00054 ss << "[JetPlusTrackCorrector::" << __func__ << "] Configuration for JPT corrector: " << std::endl
00055 << " Particles" << std::endl
00056 << " UsePions : " << ( usePions_ ? "true" : "false" ) << std::endl
00057 << " UseMuons : " << ( useMuons_ ? "true" : "false" ) << std::endl
00058 << " UseElecs : " << ( useElecs_ ? "true" : "false" ) << std::endl
00059 << " Corrections" << std::endl
00060 << " UseInConeTracks : " << ( useInConeTracks_ ? "true" : "false" ) << std::endl
00061 << " UseOutOfConeTracks : " << ( useOutOfConeTracks_ ? "true" : "false" ) << std::endl
00062 << " UseOutOfVertexTracks : " << ( useOutOfVertexTracks_ ? "true" : "false" ) << std::endl
00063 << " ResponseMap : " << pset.getParameter<std::string>("ResponseMap") << std::endl
00064 << " Efficiency" << std::endl
00065 << " UsePionEfficiency : " << ( useEff_ ? "true" : "false" ) << std::endl
00066 << " EfficiencyMap : " << pset.getParameter<std::string>("EfficiencyMap") << std::endl
00067 << " LeakageMap : " << pset.getParameter<std::string>("LeakageMap") << std::endl
00068 << " Tracks" << std::endl
00069 << " JetTracksAtVertex : " << jetTracksAtVertex_ << std::endl
00070 << " JetTracksAtCalo : " << jetTracksAtCalo_ << std::endl
00071 << " JetSplitMerge : " << jetSplitMerge_ << std::endl
00072 << " UseTrackQuality : " << ( useTrackQuality_ ? "true" : "false" ) << std::endl
00073 << " Collections" << std::endl
00074 << " Muons : " << muons_ << std::endl
00075 << " Electrons : " << electrons_ << std::endl
00076 << " Vectorial" << std::endl
00077 << " UseTracksAndResponse : " << ( ( vectorial_ && vecResponse_ ) ? "true" : "false" ) << std::endl
00078 << " UseTracksOnly : " << ( ( vectorial_ && !vecResponse_ ) ? "true" : "false" );
00079 edm::LogInfo("JetPlusTrackCorrector") << ss.str();
00080 }
00081
00082 if ( !useInConeTracks_ ||
00083 !useOutOfConeTracks_ ||
00084 !useOutOfVertexTracks_ ) {
00085 std::stringstream ss;
00086 ss << "[JetPlusTrackCorrector::" << __func__ << "]"
00087 << " You are using JPT algorithm in a non-standard way!" << std::endl
00088 << " UseInConeTracks : " << ( useInConeTracks_ ? "true" : "false" ) << std::endl
00089 << " UseOutOfConeTracks : " << ( useOutOfConeTracks_ ? "true" : "false" ) << std::endl
00090 << " UseOutOfVertexTracks : " << ( useOutOfVertexTracks_ ? "true" : "false" );
00091 edm::LogWarning("JetPlusTrackCorrector") << ss.str();
00092 }
00093 }
00094
00095
00096
00097 JetPlusTrackCorrector::~JetPlusTrackCorrector() {;}
00098
00099
00100
00101 double JetPlusTrackCorrector::correction( const reco::Jet& fJet, const reco::Jet& fJetcalo,
00102 const edm::Event& event,
00103 const edm::EventSetup& setup,
00104 P4& corrected) const
00105 {
00106
00107
00108
00109 theResponseOfChargedWithEff = 0.;
00110 theResponseOfChargedWithoutEff = 0.;
00111 theSumPtWithEff = 0.;
00112 theSumPtWithoutEff = 0.;
00113 theSumEnergyWithEff = 0.;
00114 theSumEnergyWithoutEff = 0.;
00115
00116
00117 corrected = fJet.p4();
00118
00119
00120 if ( !canCorrect(fJet) ) { return 1.; }
00121
00122
00123
00124
00125 MatchedTracks pions;
00126 MatchedTracks muons;
00127 MatchedTracks elecs;
00128 bool ok = matchTracks( fJetcalo, event, setup, pions, muons, elecs );
00129 if ( !ok ) { return 1.; }
00130
00131
00132
00133
00134 if ( verbose_ ) {
00135 edm::LogInfo("JetPlusTrackCorrector")
00136 << "[JetPlusTrackCorrector::" << __func__ << "]"
00137 << " Applying JPT corrections...";
00138 }
00139
00140
00141 if ( usePions_ ) { corrected += pionCorrection( fJet.p4(), pions ); }
00142
00143
00144 if ( useMuons_ ) { corrected += muonCorrection( fJet.p4(), muons ); }
00145
00146
00147 if ( useElecs_ ) { corrected += elecCorrection( fJet.p4(), elecs ); }
00148
00149
00150 if ( vectorial_ && !vecResponse_ ) { corrected = jetDirFromTracks( corrected, pions, muons, elecs ); }
00151
00152
00153 double scale = checkScale( fJet.p4(), corrected );
00154
00155
00156 if ( verbose_ ) {
00157 std::stringstream ss;
00158 ss << "Total correction:" << std::endl
00159 << std::fixed << std::setprecision(6)
00160 << " Uncorrected (Px,Py,Pz,E) : "
00161 << "(" << fJet.px()
00162 << "," << fJet.py()
00163 << "," << fJet.pz()
00164 << "," << fJet.energy()
00165 << ")" << std::endl
00166 << " Corrected (Px,Py,Pz,E) : "
00167 << "(" << corrected.px()
00168 << "," << corrected.py()
00169 << "," << corrected.pz()
00170 << "," << corrected.energy()
00171 << ")" << std::endl
00172 << " Uncorrected (Pt,Eta,Phi,M) : "
00173 << "(" << fJet.pt()
00174 << "," << fJet.eta()
00175 << "," << fJet.phi()
00176 << "," << fJet.mass()
00177 << ")" << std::endl
00178 << " Corrected (Pt,Eta,Phi,M) : "
00179 << "(" << corrected.pt()
00180 << "," << corrected.eta()
00181 << "," << corrected.phi()
00182 << "," << corrected.mass()
00183 << ")" << std::endl
00184 << " Scalar correction to E : " << scale << std::endl
00185 << " Scalar correction to Et : " << ( fJet.et() > 0. ? corrected.Et() / fJet.et() : 1. );
00186 edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
00187 }
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 return scale;
00206
00207 }
00208
00209
00210
00211 double JetPlusTrackCorrector::correction( const reco::Jet& jet ) const {
00212 edm::LogError("JetPlusTrackCorrector")
00213 << "JetPlusTrackCorrector can be run on entire event only";
00214 return 1.;
00215 }
00216
00217
00218
00219 double JetPlusTrackCorrector::correction( const reco::Particle::LorentzVector& jet ) const {
00220 edm::LogError("JetPlusTrackCorrector")
00221 << "JetPlusTrackCorrector can be run on entire event only";
00222 return 1.;
00223 }
00224
00225
00226
00227 bool JetPlusTrackCorrector::matchTracks( const reco::Jet& fJet,
00228 const edm::Event& event,
00229 const edm::EventSetup& setup,
00230 jpt::MatchedTracks& pions,
00231 jpt::MatchedTracks& muons,
00232 jpt::MatchedTracks& elecs ) const {
00233
00234
00235 JetTracks jet_tracks;
00236 bool ok = jetTrackAssociation( fJet, event, setup, jet_tracks );
00237
00238
00239
00240 if ( !ok ) { return false; }
00241
00242
00243 matchTracks( jet_tracks, event, pions, muons, elecs );
00244
00245
00246 if ( verbose_ ) {
00247 std::stringstream ss;
00248 ss << "Number of tracks:" << std::endl
00249 << " In-cone at Vertex : " << jet_tracks.vertex_.size() << std::endl
00250 << " In-cone at CaloFace : " << jet_tracks.caloFace_.size();
00251 edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
00252 }
00253
00254 return true;
00255
00256 }
00257
00258
00259
00260 bool JetPlusTrackCorrector::jetTrackAssociation( const reco::Jet& fJet,
00261 const edm::Event& event,
00262 const edm::EventSetup& setup,
00263 JetTracks& trks ) const {
00264
00265
00266 trks.clear();
00267
00268
00269 if ( !jetTracksAtVertex_.label().empty() && !jetTracksAtCalo_.label().empty() ) {
00270
00271
00272
00273 return jtaUsingEventData( fJet, event, trks );
00274 } else {
00275
00276 edm::LogWarning("PatJPTCorrector")
00277 << "[JetPlusTrackCorrector::" << __func__ << "]"
00278 << " Empty label for the reco::JetTracksAssociation::Containers"
00279 << std::endl
00280 << " InputTag for JTA \"at vertex\" (label:instance:process) \""
00281 << jetTracksAtVertex_.label() << ":"
00282 << jetTracksAtVertex_.instance() << ":"
00283 << jetTracksAtVertex_.process() << "\""
00284 << std::endl
00285 << " InputTag for JTA \"at calo face\" (label:instance:process) \""
00286 << jetTracksAtCalo_.label() << ":"
00287 << jetTracksAtCalo_.instance() << ":"
00288 << jetTracksAtCalo_.process() << "\"";
00289 return false;
00290 }
00291
00292 }
00293
00294
00295
00296 bool JetPlusTrackCorrector::jtaUsingEventData( const reco::Jet& fJet,
00297 const edm::Event& event,
00298 JetTracks& trks ) const {
00299
00300
00301 edm::Handle<reco::JetTracksAssociation::Container> jetTracksAtVertex;
00302 event.getByLabel( jetTracksAtVertex_, jetTracksAtVertex );
00303
00304 if ( !jetTracksAtVertex.isValid() || jetTracksAtVertex.failedToGet() ) {
00305 if ( verbose_ && edm::isDebugEnabled() ) {
00306 edm::LogWarning("JetPlusTrackCorrector")
00307 << "[JetPlusTrackCorrector::" << __func__ << "]"
00308 << " Invalid handle to reco::JetTracksAssociation::Container (for Vertex)"
00309 << " with InputTag (label:instance:process) \""
00310 << jetTracksAtVertex_.label() << ":"
00311 << jetTracksAtVertex_.instance() << ":"
00312 << jetTracksAtVertex_.process() << "\"";
00313 }
00314 return false;
00315 }
00316
00317
00318
00319
00320
00321 const reco::JetTracksAssociation::Container jtV = *( jetTracksAtVertex.product() );
00322 TrackRefs excluded;
00323 if ( jetSplitMerge_ < 0 ) { trks.vertex_ = reco::JetTracksAssociation::getValue( jtV, fJet ); }
00324 else { rebuildJta( fJet, jtV, trks.vertex_, excluded ); }
00325
00326
00327 if ( trks.vertex_.empty() ) { return false; }
00328
00329
00330 edm::Handle<reco::JetTracksAssociation::Container> jetTracksAtCalo;
00331 event.getByLabel( jetTracksAtCalo_, jetTracksAtCalo );
00332
00333 if ( !jetTracksAtCalo.isValid() || jetTracksAtCalo.failedToGet() ) {
00334 if ( verbose_ && edm::isDebugEnabled() ) {
00335 edm::LogWarning("JetPlusTrackCorrector")
00336 << "[JetPlusTrackCorrector::" << __func__ << "]"
00337 << " Invalid handle to reco::JetTracksAssociation::Container (for CaloFace)"
00338 << " with InputTag (label:instance:process) \""
00339 << jetTracksAtCalo_.label() << ":"
00340 << jetTracksAtCalo_.instance() << ":"
00341 << jetTracksAtCalo_.process() << "\"";
00342 }
00343 return false;
00344 }
00345
00346
00347 const reco::JetTracksAssociation::Container jtC = *( jetTracksAtCalo.product() );
00348 if ( jetSplitMerge_ < 0 ) { trks.caloFace_ = reco::JetTracksAssociation::getValue( jtC, fJet ); }
00349 else { excludeJta( fJet, jtC, trks.caloFace_, excluded ); }
00350
00351
00352
00353
00354 return true;
00355
00356 }
00357
00358
00359
00360 bool JetPlusTrackCorrector::getMuons( const edm::Event& event, edm::Handle<RecoMuons>& reco_muons ) const {
00361 event.getByLabel( muons_, reco_muons );
00362 if ( !reco_muons.isValid() || reco_muons.failedToGet() ) {
00363 edm::LogError("JetPlusTrackCorrector")
00364 << "[JetPlusTrackCorrector::" << __func__ << "]"
00365 << " Invalid handle to reco::Muon collection"
00366 << " with InputTag (label:instance:process) \""
00367 << muons_.label() << ":"
00368 << muons_.instance() << ":"
00369 << muons_.process() << "\"";
00370 return false;
00371 }
00372 return true;
00373 }
00374
00375
00376
00377 void JetPlusTrackCorrector::matchTracks( const JetTracks& jet_tracks,
00378 const edm::Event& event,
00379 MatchedTracks& pions,
00380 MatchedTracks& muons,
00381 MatchedTracks& elecs ) const {
00382
00383
00384 pions.clear();
00385 muons.clear();
00386 elecs.clear();
00387
00388
00389
00390 vertex_=reco::Particle::Point(0,0,0);
00391 edm::Handle<reco::VertexCollection> pvCollection;
00392 event.getByLabel(srcPVs_,pvCollection);
00393 if ( pvCollection.isValid() && pvCollection->size()>0 ) vertex_=pvCollection->begin()->position();
00394
00395
00396 edm::Handle<RecoMuons> reco_muons;
00397 bool found_reco_muons = true;
00398 if ( useMuons_ ) { getMuons( event, reco_muons ); }
00399
00400
00401 edm::Handle<RecoElectrons> reco_elecs;
00402 edm::Handle<RecoElectronIds> reco_elec_ids;
00403 bool found_reco_elecs = true;
00404 if ( useElecs_ ) { getElectrons( event, reco_elecs, reco_elec_ids ); }
00405
00406
00407 if ( !found_reco_muons || !found_reco_elecs ) {
00408 edm::LogError("JetPlusTrackCorrector")
00409 << "[JetPlusTrackCorrector::" << __func__ << "]"
00410 << " Unable to access RECO collections for muons and electrons";
00411 return;
00412 }
00413
00414
00415 {
00416 TrackRefs::const_iterator itrk = jet_tracks.vertex_.begin();
00417 TrackRefs::const_iterator jtrk = jet_tracks.vertex_.end();
00418
00419 for ( ; itrk != jtrk; ++itrk ) {
00420
00421 if ( failTrackQuality(itrk) ) { continue; }
00422
00423
00424
00425 TrackRefs::iterator it = jet_tracks.caloFace_.end();
00426 bool found = findTrack( jet_tracks, itrk, it );
00427
00428
00429
00430
00431 bool is_muon = useMuons_ && matchMuons( itrk, reco_muons );
00432 bool is_ele = useElecs_ && matchElectrons( itrk, reco_elecs, reco_elec_ids );
00433
00434 if ( found ) {
00435 if ( is_muon ) { muons.inVertexInCalo_.push_back(*it); }
00436 else if ( is_ele ) { elecs.inVertexInCalo_.push_back(*it); }
00437 else { pions.inVertexInCalo_.push_back(*it); }
00438 } else {
00439 if ( is_muon ) { muons.inVertexOutOfCalo_.push_back(*itrk); }
00440 else if ( is_ele ) { elecs.inVertexOutOfCalo_.push_back(*itrk); }
00441 else {
00442 pions.inVertexOutOfCalo_.push_back(*itrk);
00443
00444
00445 }
00446 }
00447 }
00448 }
00449
00450
00451 {
00452 TrackRefs::iterator itrk = jet_tracks.caloFace_.begin();
00453 TrackRefs::iterator jtrk = jet_tracks.caloFace_.end();
00454 for ( ; itrk != jtrk; ++itrk ) {
00455
00456 if ( failTrackQuality(itrk) ) { continue; }
00457
00458 if ( !tracksInCalo( pions, muons, elecs ) ) { continue; }
00459
00460 bool found = findTrack( pions, muons, elecs, itrk );
00461
00462 if ( !found ) {
00463
00464 bool is_muon = useMuons_ && matchMuons( itrk, reco_muons );
00465 bool is_ele = useElecs_ && matchElectrons( itrk, reco_elecs, reco_elec_ids );
00466
00467 if ( is_muon ) { muons.outOfVertexInCalo_.push_back(*itrk); }
00468 else if ( is_ele ) { elecs.outOfVertexInCalo_.push_back(*itrk); }
00469 else { pions.outOfVertexInCalo_.push_back(*itrk); }
00470
00471 }
00472 }
00473 }
00474
00475 if ( verbose_ && edm::isDebugEnabled() ) {
00476 std::stringstream ss;
00477 ss << "[JetPlusTrackCorrector::" << __func__ << "] Number of tracks:" << std::endl
00478 << " In-cone at Vertex and in-cone at CaloFace:" << std::endl
00479 << " Pions : " << pions.inVertexInCalo_.size() << std::endl
00480 << " Muons : " << muons.inVertexInCalo_.size() << std::endl
00481 << " Electrons : " << elecs.inVertexInCalo_.size() << std::endl
00482 << " In-cone at Vertex and out-of-cone at CaloFace:" << std::endl
00483 << " Pions : " << pions.inVertexOutOfCalo_.size() << std::endl
00484 << " Muons : " << muons.inVertexOutOfCalo_.size() << std::endl
00485 << " Electrons : " << elecs.inVertexOutOfCalo_.size() << std::endl
00486 << " Out-of-cone at Vertex and in-cone at CaloFace:" << std::endl
00487 << " Pions : " << pions.outOfVertexInCalo_.size() << std::endl
00488 << " Muons : " << muons.outOfVertexInCalo_.size() << std::endl
00489 << " Electrons : " << elecs.outOfVertexInCalo_.size() << std::endl;
00490 LogTrace("JetPlusTrackCorrector") << ss.str();
00491 }
00492
00493 }
00494
00495
00496
00497 bool JetPlusTrackCorrector::getElectrons( const edm::Event& event,
00498 edm::Handle<RecoElectrons>& reco_elecs,
00499 edm::Handle<RecoElectronIds>& reco_elec_ids ) const {
00500 event.getByLabel( electrons_, reco_elecs );
00501 if ( !reco_elecs.isValid() || reco_elecs.failedToGet() ) {
00502 edm::LogError("JetPlusTrackCorrector")
00503 << "[JetPlusTrackCorrector::" << __func__ << "]"
00504 << " Invalid handle to reco::GsfElectron collection"
00505 << " with InputTag (label:instance:process) \""
00506 << electrons_.label() << ":"
00507 << electrons_.instance() << ":"
00508 << electrons_.process() << "\"";
00509 return false;
00510 }
00511 event.getByLabel( electronIds_, reco_elec_ids );
00512 if ( !reco_elec_ids.isValid() || reco_elec_ids.failedToGet() ) {
00513 edm::LogError("JetPlusTrackCorrector")
00514 << "[JetPlusTrackCorrector::" << __func__ << "]"
00515 << " Invalid handle to reco::GsfElectron collection"
00516 << " with InputTag (label:instance:process) \""
00517 << electronIds_.label() << ":"
00518 << electronIds_.instance() << ":"
00519 << electronIds_.process() << "\"";
00520 return false;
00521 }
00522 return true;
00523 }
00524
00525
00526
00527 bool JetPlusTrackCorrector::failTrackQuality( TrackRefs::const_iterator itrk ) const {
00528
00529
00530
00531 bool retcode = false;
00532
00533 if ( useTrackQuality_ && !(*itrk)->quality(trackQuality_) ) {
00534 retcode = true; return retcode;
00535 }
00536 if(((*itrk)->ptError()/(*itrk)->pt()) > ptErrorQuality_) {
00537 retcode = true; return retcode;
00538 }
00539 if(fabs((*itrk)->dz(vertex_)) > dzVertexCut_) {
00540 retcode = true; return retcode;
00541 }
00542
00543 return retcode;
00544 }
00545
00546
00547
00548 bool JetPlusTrackCorrector::findTrack( const JetTracks& jet_tracks,
00549 TrackRefs::const_iterator itrk,
00550 TrackRefs::iterator& it ) const {
00551 it = find( jet_tracks.caloFace_.begin(),
00552 jet_tracks.caloFace_.end(),
00553 *itrk );
00554 if ( it != jet_tracks.caloFace_.end() ) { return true; }
00555 else { return false; }
00556 }
00557
00558
00559
00560 bool JetPlusTrackCorrector::findTrack( const MatchedTracks& pions,
00561 const MatchedTracks& muons,
00562 const MatchedTracks& elecs,
00563 TrackRefs::const_iterator itrk ) const {
00564 TrackRefs::iterator ip = find( pions.inVertexInCalo_.begin(),
00565 pions.inVertexInCalo_.end(),
00566 *itrk );
00567 TrackRefs::iterator im = find( muons.inVertexInCalo_.begin(),
00568 muons.inVertexInCalo_.end(),
00569 *itrk );
00570 TrackRefs::iterator ie = find( elecs.inVertexInCalo_.begin(),
00571 elecs.inVertexInCalo_.end(),
00572 *itrk );
00573 if ( ip == pions.inVertexInCalo_.end() &&
00574 im == muons.inVertexInCalo_.end() &&
00575 ie == elecs.inVertexInCalo_.end() ) { return false; }
00576 else { return true; }
00577 }
00578
00579
00580
00581 bool JetPlusTrackCorrector::tracksInCalo( const MatchedTracks& pions,
00582 const MatchedTracks& muons,
00583 const MatchedTracks& elecs ) const {
00584 if ( !pions.inVertexInCalo_.empty() ||
00585 !muons.inVertexInCalo_.empty() ||
00586 !elecs.inVertexInCalo_.empty() ) { return true; }
00587 else { return false; }
00588 }
00589
00590
00591
00592 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::pionCorrection( const P4& jet,
00593 const MatchedTracks& pions ) const {
00594
00595 P4 corr_pions;
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607 P4 corr_pions_in_cone;
00608 P4 corr_pions_eff_in_cone;
00609 Efficiency in_cone( responseMap(), efficiencyMap(), leakageMap() );
00610
00611 if ( useInConeTracks_ ) {
00612 corr_pions_in_cone = pionCorrection( jet, pions.inVertexInCalo_, in_cone, true, true );
00613 corr_pions += corr_pions_in_cone;
00614 if ( useEff_ ) {
00615 corr_pions_eff_in_cone = pionEfficiency( jet, in_cone, true );
00616 corr_pions += corr_pions_eff_in_cone;
00617 }
00618 }
00619
00620
00621
00622 P4 corr_pions_out_of_cone;
00623 P4 corr_pions_eff_out_of_cone;
00624 Efficiency out_of_cone( responseMap(), efficiencyMap(), leakageMap() );
00625
00626 if ( useOutOfConeTracks_ ) {
00627 corr_pions_out_of_cone = pionCorrection( jet, pions.inVertexOutOfCalo_, out_of_cone, true, false );
00628 corr_pions += corr_pions_out_of_cone;
00629 if ( useEff_ ) {
00630 corr_pions_eff_out_of_cone = pionEfficiency( jet, out_of_cone, false );
00631 corr_pions += corr_pions_eff_out_of_cone;
00632 }
00633 }
00634
00635
00636
00637 P4 corr_pions_out_of_vertex;
00638 Efficiency not_used( responseMap(), efficiencyMap(), leakageMap() );
00639
00640 if ( useOutOfVertexTracks_ ) {
00641 corr_pions_out_of_vertex = pionCorrection( jet, pions.outOfVertexInCalo_, not_used, false, true );
00642 corr_pions += corr_pions_out_of_vertex;
00643 }
00644
00645 if ( verbose_ ) {
00646 std::stringstream ss;
00647 ss << " Pion corrections:" << std::endl
00648 << " In/In : " << "(" << pions.inVertexInCalo_.size() << ") " << corr_pions_in_cone.energy() << std::endl
00649 << " In/Out : " << "(" << pions.inVertexOutOfCalo_.size() << ") " << corr_pions_out_of_cone.energy() << std::endl
00650 << " Out/In : " << "(" << pions.outOfVertexInCalo_.size() << ") " << corr_pions_out_of_vertex.energy() << std::endl;
00651 if ( useEff_ ) {
00652 ss << " Pion efficiency corrections:" << std::endl
00653 << " In/In : " << "(" << pions.inVertexInCalo_.size() << ") " << corr_pions_eff_in_cone.energy() << std::endl
00654 << " In/Out : " << "(" << pions.inVertexOutOfCalo_.size() << ") " << corr_pions_eff_out_of_cone.energy();
00655 }
00656 edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
00657 }
00658
00659 return corr_pions;
00660
00661 }
00662
00663
00664
00665 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::muonCorrection( const P4& jet,
00666 const MatchedTracks& muons ) const {
00667
00668 P4 corr_muons;
00669
00670 P4 corr_muons_in_cone;
00671 P4 corr_muons_out_of_cone;
00672 P4 corr_muons_out_of_vertex;
00673
00674 if ( useInConeTracks_ ) {
00675 corr_muons_in_cone = muonCorrection( jet, muons.inVertexInCalo_, true, true );
00676 corr_muons += corr_muons_in_cone;
00677 }
00678
00679 if ( useOutOfConeTracks_ ) {
00680 corr_muons_out_of_cone = muonCorrection( jet, muons.inVertexOutOfCalo_, true, false );
00681 corr_muons += corr_muons_out_of_cone;
00682 }
00683
00684 if ( useOutOfVertexTracks_ ) {
00685 corr_muons_out_of_vertex = muonCorrection( jet, muons.outOfVertexInCalo_, false, true );
00686 corr_muons += corr_muons_out_of_vertex;
00687 }
00688
00689 if ( verbose_ ) {
00690 std::stringstream ss;
00691 ss << " Muon corrections:" << std::endl
00692 << " In/In : " << "(" << muons.inVertexInCalo_.size() << ") " << corr_muons_in_cone.energy() << std::endl
00693 << " In/Out : " << "(" << muons.inVertexOutOfCalo_.size() << ") " << corr_muons_out_of_cone.energy() << std::endl
00694 << " Out/In : " << "(" << muons.outOfVertexInCalo_.size() << ") " << corr_muons_out_of_vertex.energy();
00695 edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
00696 }
00697
00698 return corr_muons;
00699
00700 }
00701
00702
00703
00704 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::elecCorrection( const P4& jet,
00705 const MatchedTracks& elecs ) const {
00706
00707 P4 null;
00708
00709 if ( verbose_ ) {
00710 std::stringstream ss;
00711 ss << " Electron corrections:" << std::endl
00712 << " In/In : " << "(" << elecs.inVertexInCalo_.size() << ") " << null.energy() << std::endl
00713 << " In/Out : " << "(" << elecs.inVertexOutOfCalo_.size() << ") " << null.energy() << std::endl
00714 << " Out/In : " << "(" << elecs.outOfVertexInCalo_.size() << ") " << null.energy();
00715 edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
00716 }
00717
00718 return null;
00719
00720 }
00721
00722
00723
00724 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::jetDirFromTracks( const P4& corrected,
00725 const MatchedTracks& pions,
00726 const MatchedTracks& muons,
00727 const MatchedTracks& elecs ) const {
00728
00729
00730 P4 corr;
00731
00732 bool tracks_present = false;
00733
00734
00735
00736 if ( !pions.inVertexInCalo_.empty() ) {
00737 TrackRefs::iterator itrk = pions.inVertexInCalo_.begin();
00738 TrackRefs::iterator jtrk = pions.inVertexInCalo_.end();
00739 for ( ; itrk != jtrk; ++itrk ) {
00740 corr += PtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0. );
00741 tracks_present = true;
00742 }
00743 }
00744
00745 if ( !pions.inVertexOutOfCalo_.empty() ) {
00746 TrackRefs::iterator itrk = pions.inVertexOutOfCalo_.begin();
00747 TrackRefs::iterator jtrk = pions.inVertexOutOfCalo_.end();
00748 for ( ; itrk != jtrk; ++itrk ) {
00749 corr += PtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0. );
00750 tracks_present = true;
00751 }
00752 }
00753
00754
00755
00756 if ( !muons.inVertexInCalo_.empty() ) {
00757 TrackRefs::iterator itrk = muons.inVertexInCalo_.begin();
00758 TrackRefs::iterator jtrk = muons.inVertexInCalo_.end();
00759 for ( ; itrk != jtrk; ++itrk ) {
00760 corr += PtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0. );
00761 tracks_present = true;
00762 }
00763 }
00764
00765 if ( !muons.inVertexOutOfCalo_.empty() ) {
00766 TrackRefs::iterator itrk = muons.inVertexOutOfCalo_.begin();
00767 TrackRefs::iterator jtrk = muons.inVertexOutOfCalo_.end();
00768 for ( ; itrk != jtrk; ++itrk ) {
00769 corr += PtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0. );
00770 tracks_present = true;
00771 }
00772 }
00773
00774
00775
00776 if ( !elecs.inVertexInCalo_.empty() ) {
00777 TrackRefs::iterator itrk = elecs.inVertexInCalo_.begin();
00778 TrackRefs::iterator jtrk = elecs.inVertexInCalo_.end();
00779 for ( ; itrk != jtrk; ++itrk ) {
00780 corr += PtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0. );
00781 tracks_present = true;
00782 }
00783 }
00784
00785 if ( !elecs.inVertexOutOfCalo_.empty() ) {
00786 TrackRefs::iterator itrk = elecs.inVertexOutOfCalo_.begin();
00787 TrackRefs::iterator jtrk = elecs.inVertexOutOfCalo_.end();
00788 for ( ; itrk != jtrk; ++itrk ) {
00789 corr += PtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0. );
00790 tracks_present = true;
00791 }
00792 }
00793
00794
00795
00796 if ( !tracks_present ) { corr = corrected; }
00797 else {
00798 corr *= ( corr.P() > 0. ? corrected.P() / corr.P() : 1. );
00799 corr = P4( corr.px(), corr.py(), corr.pz(), corrected.energy() );
00800 }
00801
00802 return corr;
00803
00804 }
00805
00806
00807
00808 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::calculateCorr( const P4& jet,
00809 const TrackRefs& tracks,
00810 jpt::Efficiency& eff,
00811 bool in_cone_at_vertex,
00812 bool in_cone_at_calo_face,
00813 double mass,
00814 bool is_pion,
00815 double mip ) const {
00816
00817
00818 P4 correction;
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834 eff.reset();
00835
00836 double theSumResp = 0;
00837 double theSumPt = 0;
00838 double theSumEnergy = 0;
00839
00840
00841
00842 if ( !tracks.empty() ) {
00843 TrackRefs::iterator itrk = tracks.begin();
00844 TrackRefs::iterator jtrk = tracks.end();
00845
00846 for ( ; itrk != jtrk; ++itrk ) {
00847
00848
00849 if ( in_cone_at_calo_face && is_pion && (*itrk)->pt() >= 50. ) { continue; }
00850
00851
00852 P4 inner;
00853 if ( vectorial_ && vecResponse_ ) {
00854 inner = PtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), mass );
00855 } else {
00856 double energy = sqrt( (*itrk)->px() * (*itrk)->px() +
00857 (*itrk)->py() * (*itrk)->py() +
00858 (*itrk)->pz() * (*itrk)->pz() +
00859 mass * mass );
00860 inner = ( jet.energy() > 0. ? energy / jet.energy() : 1. ) * jet;
00861 }
00862
00863
00864 if ( in_cone_at_vertex ) { correction += inner; }
00865
00866
00867 double eta = fabs( (*itrk)->eta() );
00868 double pt = fabs( (*itrk)->pt() );
00869 uint32_t ieta = response_.etaBin( eta );
00870 uint32_t ipt = response_.ptBin( pt );
00871
00872
00873 if ( is_pion && ( ieta == response_.nEtaBins() || ipt == response_.nPtBins() ) ) { continue; }
00874
00875
00876 P4 outer;
00877 if ( in_cone_at_calo_face ) {
00878 if ( vectorial_ && vecResponse_ ) {
00879
00880 double outer_pt = (*itrk)->pt();
00881 double outer_eta = (*itrk)->eta();
00882 double outer_phi = (*itrk)->phi();
00883 if ( (*itrk)->extra().isNonnull() ) {
00884 outer_pt = (*itrk)->pt();
00885 outer_eta = (*itrk)->outerPosition().eta();
00886 outer_phi = (*itrk)->outerPosition().phi();
00887 }
00888 outer = PtEtaPhiM( outer_pt, outer_eta, outer_phi, mass );
00889
00890 if ( !is_pion ) { outer *= ( outer.energy() > 0. ? mip / outer.energy() : 1. ); }
00891 else { outer *= ( outer.energy() > 0. ? inner.energy() / outer.energy() : 1. ); }
00892 } else {
00893
00894 if ( !is_pion ) { outer = ( jet.energy() > 0. ? mip / jet.energy() : 1. ) * jet; }
00895 else { outer = inner; }
00896 }
00897 if ( is_pion ) { outer *= response_.value(ieta,ipt); }
00898 correction -= outer;
00899
00900
00901 theSumResp += response_.value(ieta,ipt);
00902 }
00903
00904
00905 theSumPt += inner.pt();
00906 theSumEnergy += inner.energy();
00907
00908
00909 if ( is_pion ) { eff.addE( ieta, ipt, inner.energy() ); }
00910
00911
00912 if ( verbose_ && edm::isDebugEnabled() ) {
00913 std::stringstream temp;
00914 temp << " Response[" << ieta << "," << ipt << "]";
00915 std::stringstream ss;
00916 ss << "[JetPlusTrackCorrector::" << __func__ << "]" << std::endl
00917 << " Track eta / pt : " << eta << " / " << pt << std::endl
00918 << temp.str() << std::setw(21-temp.str().size()) << " : "
00919 << response_.value(ieta,ipt) << std::endl
00920 << " Track momentum added : " << inner.energy() << std::endl
00921 << " Response subtracted : " << outer.energy() << std::endl
00922 << " Energy correction : " << correction.energy();
00923 LogDebug("JetPlusTrackCorrector") << ss.str();
00924 }
00925
00926 }
00927 }
00928
00929 if( in_cone_at_vertex ) {
00930
00931 theResponseOfChargedWithEff += theSumResp;
00932 theResponseOfChargedWithoutEff += theSumResp;
00933 theSumPtWithEff += theSumPt;
00934 theSumPtWithoutEff += theSumPt;
00935 theSumEnergyWithEff += theSumEnergy;
00936 theSumEnergyWithoutEff += theSumEnergy;
00937
00938 }
00939
00940
00941
00942
00943
00944
00945
00946
00947 return correction;
00948
00949 }
00950
00951
00952
00953 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::pionEfficiency( const P4& jet,
00954 const Efficiency& eff,
00955 bool in_cone_at_calo_face ) const {
00956
00957
00958 P4 correction;
00959
00960 double theSumResp = 0;
00961 double theSumPt = 0;
00962 double theSumEnergy = 0;
00963
00964
00965 for ( uint32_t ieta = 0; ieta < response_.nEtaBins()-1; ++ieta ) {
00966 for ( uint32_t ipt = 0; ipt < response_.nPtBins()-1; ++ipt ) {
00967
00968
00969 if ( !eff.nTrks(ieta,ipt) ) { continue; }
00970
00971 for ( uint16_t ii = 0; ii < 2; ++ii ) {
00972
00973
00974 double corr = 0.;
00975 if ( ii == 0 ) { corr = eff.outOfConeCorr( ieta, ipt );}
00976 else if ( ii == 1 && in_cone_at_calo_face ) { corr = eff.inConeCorr( ieta, ipt );}
00977 else { continue; }
00978
00979
00980 P4 corr_p4;
00981 if ( vectorial_ && vecResponse_ ) {
00982 double corr_eta = response_.binCenterEta(ieta);
00983 double corr_phi = jet.phi();
00984 double corr_pt = response_.binCenterPt(ipt);
00985 corr_p4 = PtEtaPhiM( corr_pt, corr_eta, corr_phi, pionMass_ );
00986 corr_p4 *= ( corr_p4.energy() > 0. ? corr / corr_p4.energy() : 1. );
00987 } else {
00988 corr_p4 = ( jet.energy() > 0. ? corr / jet.energy() : 1. ) * jet;
00989 }
00990
00991
00992 if ( ii == 0 ) { correction += corr_p4; theSumPt += corr_p4.pt(); theSumEnergy += corr_p4.energy();}
00993 else if ( ii == 1 ) { correction -= corr_p4; theSumResp += corr_p4.energy();}
00994
00995 }
00996
00997 }
00998 }
00999
01000 theResponseOfChargedWithEff += theSumResp;
01001 theSumPtWithEff += theSumPt;
01002 theSumEnergyWithEff += theSumEnergy;
01003
01004
01005
01006
01007
01008
01009
01010
01011 return correction;
01012
01013 }
01014
01015
01016
01017 bool JetPlusTrackCorrector::matchMuons( TrackRefs::const_iterator itrk,
01018 const edm::Handle<RecoMuons>& muons ) const {
01019
01020 if ( muons->empty() ) { return false; }
01021
01022 RecoMuons::const_iterator imuon = muons->begin();
01023 RecoMuons::const_iterator jmuon = muons->end();
01024 for ( ; imuon != jmuon; ++imuon ) {
01025
01026 if ( imuon->innerTrack().isNull() ||
01027 !muon::isGoodMuon(*imuon,muon::TMLastStationTight) ||
01028 imuon->innerTrack()->pt() < 3.0 ) { continue; }
01029
01030 if ( itrk->id() != imuon->innerTrack().id() ) {
01031 edm::LogError("JetPlusTrackCorrector")
01032 << "[JetPlusTrackCorrector::" << __func__ << "]"
01033 << "Product id of the tracks associated to the jet " << itrk->id()
01034 <<" is different from the product id of the inner track used for muons " << imuon->innerTrack().id()
01035 << "!" << std::endl
01036 << "Cannot compare tracks from different collection. Configuration Error!";
01037 return false;
01038 }
01039
01040 if ( *itrk == imuon->innerTrack() ) { return true; }
01041 }
01042
01043 return false;
01044
01045 }
01046
01047
01048
01049 bool JetPlusTrackCorrector::matchElectrons( TrackRefs::const_iterator itrk,
01050 const edm::Handle<RecoElectrons>& elecs,
01051 const edm::Handle<RecoElectronIds>& elec_ids ) const {
01052
01053 if ( elecs->empty() ) { return false; }
01054
01055 double deltaR = 999.;
01056 double deltaRMIN = 999.;
01057
01058 uint32_t electron_index = 0;
01059 RecoElectrons::const_iterator ielec = elecs->begin();
01060 RecoElectrons::const_iterator jelec = elecs->end();
01061 for ( ; ielec != jelec; ++ielec ) {
01062
01063 edm::Ref<RecoElectrons> electron_ref( elecs, electron_index );
01064 electron_index++;
01065
01066 if ( (*elec_ids)[electron_ref] < 1.e-6 ) { continue; }
01067
01068
01069 double deltaphi = fabs( ielec->phi() - (*itrk)->momentum().phi() );
01070 if ( deltaphi > 6.283185308 ) deltaphi -= 6.283185308;
01071 if ( deltaphi > 3.141592654 ) deltaphi = 6.283185308 - deltaphi;
01072 deltaR = abs( sqrt( pow( (ielec->eta() - (*itrk)->momentum().eta()), 2 ) +
01073 pow( deltaphi , 2 ) ) );
01074 if ( deltaR < deltaRMIN ) { deltaRMIN = deltaR; }
01075
01076 }
01077
01078 if ( deltaR < 0.02 ) return true;
01079 else return false;
01080
01081 }
01082
01083
01084
01085 void JetPlusTrackCorrector::rebuildJta( const reco::Jet& fJet,
01086 const JetTracksAssociations& jtV0,
01087 TrackRefs& tracksthis,
01088 TrackRefs& Excl ) const {
01089
01090
01091
01092 tracksthis = reco::JetTracksAssociation::getValue(jtV0,fJet);
01093
01094 if(jetSplitMerge_<0) return;
01095
01096 typedef std::vector<reco::JetBaseRef>::iterator JetBaseRefIterator;
01097 std::vector<reco::JetBaseRef> theJets = reco::JetTracksAssociation::allJets(jtV0);
01098
01099 TrackRefs tracks = tracksthis;
01100 tracksthis.clear();
01101
01102
01103
01104 int tr=0;
01105
01106 for(TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++ )
01107 {
01108
01109 double dR2this = deltaR2( fJet.eta(), fJet.phi(), (**it).eta(), (**it).phi() );
01110
01111
01112
01113
01114
01115 double scalethis = dR2this;
01116 if(jetSplitMerge_ == 0) scalethis = 1./fJet.et();
01117 if(jetSplitMerge_ == 2) scalethis = dR2this/fJet.et();
01118 tr++;
01119 int flag = 1;
01120 for(JetBaseRefIterator ii = theJets.begin(); ii != theJets.end(); ii++)
01121 {
01122 if(&(**ii) == &fJet ) {continue;}
01123 double dR2 = deltaR2( (*ii)->eta(), (*ii)->phi(), (**it).eta(), (**it).phi() );
01124 double scale = dR2;
01125 if(jetSplitMerge_ == 0) scale = 1./fJet.et();
01126 if(jetSplitMerge_ == 2) scale = dR2/fJet.et();
01127 if(scale < scalethis) flag = 0;
01128
01129 if(flag == 0) {
01130
01131
01132 break;
01133 }
01134 }
01135
01136
01137 if(flag == 1) {tracksthis.push_back (*it);}else{Excl.push_back (*it);}
01138 }
01139
01140
01141 return;
01142
01143 }
01144
01145
01146
01147 void JetPlusTrackCorrector::excludeJta( const reco::Jet& fJet,
01148 const JetTracksAssociations& jtV0,
01149 TrackRefs& tracksthis,
01150 const TrackRefs& Excl ) const {
01151
01152
01153
01154 tracksthis = reco::JetTracksAssociation::getValue(jtV0,fJet);
01155 if(Excl.size() == 0) return;
01156 if(jetSplitMerge_<0) return;
01157
01158 TrackRefs tracks = tracksthis;
01159 tracksthis.clear();
01160
01161
01162
01163 for(TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++ )
01164 {
01165
01166
01167
01168 TrackRefs::iterator itold = find(Excl.begin(),Excl.end(),(*it));
01169 if(itold == Excl.end()) {
01170 tracksthis.push_back (*it);
01171 }
01172
01173
01174 }
01175
01176
01177
01178 return;
01179
01180 }
01181
01182
01183
01184 double JetPlusTrackCorrector::correctAA( const reco::Jet& fJet,
01185 const reco::TrackRefVector& trBgOutOfVertex,
01186 double& mConeSize,
01187 const reco::TrackRefVector& pioninin,
01188 const reco::TrackRefVector& pioninout,double ja,
01189 const reco::TrackRefVector& trBgOutOfCalo) const {
01190 double mScale = 1.;
01191 double NewResponse = fJet.energy();
01192
01193
01194
01195
01196
01197 if( trBgOutOfVertex.size() == 0 ) return mScale;
01198 double EnergyOfBackgroundCharged = 0.;
01199 double ResponseOfBackgroundCharged = 0.;
01200 double NumberOfBackgroundChargedVertex = 0.;
01201 double NumberOfBackgroundChargedCalo = 0.;
01202
01203
01204
01205
01206
01207
01208
01209
01210 for( reco::TrackRefVector::iterator iBgtV = trBgOutOfVertex.begin(); iBgtV != trBgOutOfVertex.end(); iBgtV++)
01211 {
01212
01213
01214
01215 double eta = fabs( (**iBgtV).eta() );
01216 double pt = fabs( (**iBgtV).pt() );
01217 uint32_t ieta = response_.etaBin( eta );
01218 uint32_t ipt = response_.ptBin( pt );
01219
01220 if(fabs(fJet.eta() -(**iBgtV).eta()) > mConeSize) continue;
01221
01222
01223
01224
01225
01226 double echarBg=sqrt((**iBgtV).px()*(**iBgtV).px()+(**iBgtV).py()*(**iBgtV).py()+(**iBgtV).pz()*(**iBgtV).pz()+0.14*0.14);
01227
01228
01229
01230
01231
01232 EnergyOfBackgroundCharged += echarBg/efficiency_.value(ieta,ipt);
01233
01234 NumberOfBackgroundChargedVertex++;
01235
01236 }
01237
01238
01239
01240
01241
01242 for( reco::TrackRefVector::iterator iBgtC = trBgOutOfCalo.begin(); iBgtC != trBgOutOfCalo.end(); iBgtC++)
01243 {
01244
01245
01246
01247 double eta = fabs( (**iBgtC).eta() );
01248 double pt = fabs( (**iBgtC).pt() );
01249 uint32_t ieta = response_.etaBin( eta );
01250 uint32_t ipt = response_.ptBin( pt );
01251
01252 if(fabs(fJet.eta() -(**iBgtC).eta()) > mConeSize) continue;
01253
01254
01255 if ( ieta >= response_.nEtaBins() ) { continue; }
01256 if ( ipt >= response_.nPtBins() ) { ipt = response_.nPtBins() - 1; }
01257
01258 double echarBg=sqrt((**iBgtC).px()*(**iBgtC).px()+(**iBgtC).py()*(**iBgtC).py()+(**iBgtC).pz()*(**iBgtC).pz()+0.14*0.14);
01259
01260 ResponseOfBackgroundCharged += echarBg*response_.value(ieta,ipt)/efficiency_.value(ieta,ipt);
01261
01262
01263
01264
01265 NumberOfBackgroundChargedCalo++;
01266
01267 }
01268
01269
01270
01271
01272
01273
01274
01275 double en = 0.;
01276 double px = 0.;
01277 double py = 0.;
01278 double pz = 0.;
01279
01280 for(reco::TrackRefVector::const_iterator it = pioninout.begin(); it != pioninout.end(); it++) {
01281
01282
01283
01284
01285 px+=(*it)->px();
01286 py+=(*it)->py();
01287 pz+=(*it)->pz();
01288 en+=sqrt((*it)->p()*(*it)->p()+0.14*0.14);
01289 }
01290
01291
01292
01293 double en_in = 0.;
01294 double px_in = 0.;
01295 double py_in = 0.;
01296 double pz_in = 0.;
01297
01298 for(reco::TrackRefVector::const_iterator it = pioninin.begin(); it != pioninin.end(); it++) {
01299
01300
01301
01302
01303 px_in+=(*it)->px();
01304 py_in+=(*it)->py();
01305 pz_in+=(*it)->pz();
01306 en_in+=sqrt((*it)->p()*(*it)->p()+0.14*0.14);
01307 }
01308
01309
01310
01311
01312
01313
01314 double SquareEtaRingWithoutJets = ja;
01315
01316
01317
01318 EnergyOfBackgroundCharged = EnergyOfBackgroundCharged/SquareEtaRingWithoutJets;
01319 ResponseOfBackgroundCharged = ResponseOfBackgroundCharged/SquareEtaRingWithoutJets;
01320 NumberOfBackgroundChargedVertex = NumberOfBackgroundChargedVertex/SquareEtaRingWithoutJets;
01321 NumberOfBackgroundChargedCalo = NumberOfBackgroundChargedCalo/SquareEtaRingWithoutJets;
01322
01323
01324
01325 EnergyOfBackgroundCharged = M_PI*mConeSize*mConeSize*EnergyOfBackgroundCharged;
01326 ResponseOfBackgroundCharged = M_PI*mConeSize*mConeSize*ResponseOfBackgroundCharged;
01327 NumberOfBackgroundChargedVertex = M_PI*mConeSize*mConeSize*NumberOfBackgroundChargedVertex;
01328 NumberOfBackgroundChargedCalo = M_PI*mConeSize*mConeSize*NumberOfBackgroundChargedCalo;
01329
01330
01331
01332 NewResponse = NewResponse - EnergyOfBackgroundCharged + ResponseOfBackgroundCharged;
01333
01334
01335
01336
01337
01338
01339
01340 mScale = NewResponse/fJet.energy();
01341 if(mScale <0.) mScale=0.;
01342 return mScale;
01343 }
01344
01345
01346 Map::Map( std::string input, bool verbose )
01347 : eta_(),
01348 pt_(),
01349 data_()
01350 {
01351
01352
01353 clear();
01354 std::vector<Element> data;
01355
01356
01357 std::string file = edm::FileInPath(input).fullPath();
01358 std::ifstream in( file.c_str() );
01359 string line;
01360 uint32_t ieta_old = 0;
01361 while ( std::getline( in, line ) ) {
01362 if ( !line.size() || line[0]=='#' ) { continue; }
01363 std::istringstream ss(line);
01364 Element temp;
01365 ss >> temp.ieta_ >> temp.ipt_ >> temp.eta_ >> temp.pt_ >> temp.val_;
01366 data.push_back(temp);
01367 if ( !ieta_old || temp.ieta_ != ieta_old ) {
01368 if ( eta_.size() < temp.ieta_+1 ) { eta_.resize(temp.ieta_+1,0.); }
01369 eta_[temp.ieta_] = temp.eta_;
01370 ieta_old = temp.ieta_;
01371 }
01372 if ( pt_.size() < temp.ipt_+1 ) { pt_.resize(temp.ipt_+1,0.); }
01373 pt_[temp.ipt_] = temp.pt_;
01374 }
01375
01376
01377 data_.resize( eta_.size(), VDouble( pt_.size(), 0. ) );
01378 std::vector<Element>::const_iterator idata = data.begin();
01379 std::vector<Element>::const_iterator jdata = data.end();
01380 for ( ; idata != jdata; ++idata ) { data_[idata->ieta_][idata->ipt_] = idata->val_; }
01381
01382
01383 if ( data_.empty() || data_[0].empty() ) {
01384 std::stringstream ss;
01385 ss << "[jpt::Map::" << __func__ << "]"
01386 << " Problem parsing map in location \""
01387 << file << "\"! ";
01388 edm::LogError("JetPlusTrackCorrector") << ss.str();
01389 }
01390
01391
01392 if ( eta_.size() != data_.size() ||
01393 pt_.size() != ( data_.empty() ? 0 : data_[0].size() ) ) {
01394 std::stringstream ss;
01395 ss << "[jpt::Map::" << __func__ << "]"
01396 << " Discrepancy b/w number of bins!";
01397 edm::LogError("JetPlusTrackCorrector") << ss.str();
01398 }
01399
01400
01401 if ( verbose && edm::isDebugEnabled() ) {
01402 std::stringstream ss;
01403 ss << "[jpt::Map::" << __func__ << "]"
01404 << " Parsed contents of map at location:" << std::endl
01405 << "\"" << file << "\"" << std::endl;
01406 print(ss);
01407 LogTrace("JetPlusTrackCorrector") << ss.str();
01408 }
01409
01410 }
01411
01412
01413
01414 Map::Map()
01415 : eta_(),
01416 pt_(),
01417 data_()
01418 {
01419 clear();
01420 }
01421
01422
01423
01424 Map::~Map() {
01425 clear();
01426 }
01427
01428
01429
01430 void Map::clear() {
01431 eta_.clear();
01432 pt_.clear();
01433 data_.clear();
01434 }
01435
01436
01437 double Map::eta( uint32_t eta_bin ) const {
01438 if ( !eta_.empty() && eta_bin < eta_.size() ) { return eta_[eta_bin]; }
01439 else {
01440
01441
01442
01443
01444
01445 return eta_[eta_.size()-1];
01446 }
01447 }
01448
01449
01450
01451 double Map::pt( uint32_t pt_bin ) const {
01452 if ( !pt_.empty() && pt_bin < pt_.size() ) { return pt_[pt_bin]; }
01453 else {
01454
01455
01456
01457
01458
01459 return pt_[pt_.size()-1];
01460 }
01461 }
01462
01463
01464
01465 double Map::binCenterEta( uint32_t eta_bin ) const {
01466 if ( !eta_.empty() && eta_bin+1 < eta_.size() ) {
01467 return eta_[eta_bin] + ( eta_[eta_bin+1] - eta_[eta_bin] ) / 2.;
01468 } else {
01469
01470
01471
01472
01473
01474 return eta_[eta_.size()-1];
01475 }
01476 }
01477
01478
01479
01480 double Map::binCenterPt( uint32_t pt_bin ) const {
01481 if ( !pt_.empty() && pt_bin+1 < pt_.size() ) {
01482 return pt_[pt_bin] + ( pt_[pt_bin+1] - pt_[pt_bin] ) / 2.;
01483 } else {
01484
01485
01486
01487
01488
01489 return pt_[pt_.size()-1];
01490 }
01491 }
01492
01493
01494
01495 uint32_t Map::etaBin( double val ) const {
01496 val = fabs( val );
01497 for ( uint32_t ieta = 0; ieta < nEtaBins()-1; ++ieta ) {
01498 if ( val > eta(ieta) && ( ieta+1 == nEtaBins() || val < eta(ieta+1) ) ) { return ieta; }
01499 }
01500 return nEtaBins();
01501 }
01502
01503
01504
01505 uint32_t Map::ptBin( double val ) const {
01506 val = fabs( val );
01507 for ( uint32_t ipt = 0; ipt < nPtBins()-1; ++ipt ) {
01508 if ( val > pt(ipt) && ( (ipt+1) == nPtBins() || val < pt(ipt+1) ) ) { return ipt; }
01509 }
01510 return nPtBins();
01511 }
01512
01513
01514
01515 double Map::value( uint32_t eta_bin, uint32_t pt_bin ) const {
01516 if ( eta_bin < data_.size() &&
01517 pt_bin < ( data_.empty() ? 0 : data_[0].size() ) ) { return data_[eta_bin][pt_bin]; }
01518 else {
01519
01520
01521
01522
01523
01524 return 1.;
01525 }
01526 }
01527
01528
01529
01530 void Map::print( std::stringstream& ss ) const {
01531 ss << " Number of bins in eta : " << data_.size() << std::endl
01532 << " Number of bins in pt : " << ( data_.empty() ? 0 : data_[0].size() ) << std::endl;
01533 VVDouble::const_iterator ieta = data_.begin();
01534 VVDouble::const_iterator jeta = data_.end();
01535 for ( ; ieta != jeta; ++ieta ) {
01536 VDouble::const_iterator ipt = ieta->begin();
01537 VDouble::const_iterator jpt = ieta->end();
01538 for ( ; ipt != jpt; ++ipt ) {
01539 uint32_t eta_bin = static_cast<uint32_t>( ieta - data_.begin() );
01540 uint32_t pt_bin = static_cast<uint32_t>( ipt - ieta->begin() );
01541 ss << " EtaBinNumber: " << eta_bin
01542 << " PtBinNumber: " << pt_bin
01543 << " EtaValue: " << eta_[ eta_bin ]
01544 << " PtValue: " << pt_[ pt_bin ]
01545 << " Value: " << data_[eta_bin][pt_bin]
01546 << std::endl;
01547 }
01548 }
01549 }
01550
01551
01552
01553 MatchedTracks::MatchedTracks()
01554 : inVertexInCalo_(),
01555 outOfVertexInCalo_(),
01556 inVertexOutOfCalo_()
01557 {
01558 clear();
01559 }
01560
01561
01562
01563 MatchedTracks::~MatchedTracks() {
01564 clear();
01565 }
01566
01567
01568
01569 void MatchedTracks::clear() {
01570 inVertexInCalo_.clear();
01571 outOfVertexInCalo_.clear();
01572 inVertexOutOfCalo_.clear();
01573 }
01574
01575
01576
01577 JetTracks::JetTracks()
01578 : vertex_(),
01579 caloFace_()
01580 {
01581 clear();
01582 }
01583
01584
01585
01586 JetTracks::~JetTracks() {
01587 clear();
01588 }
01589
01590
01591
01592 void JetTracks::clear() {
01593 vertex_.clear();
01594 caloFace_.clear();
01595 }
01596
01597
01598
01599 Efficiency::Efficiency( const jpt::Map& response,
01600 const jpt::Map& efficiency,
01601 const jpt::Map& leakage )
01602 : response_(response),
01603 efficiency_(efficiency),
01604 leakage_(leakage)
01605 {
01606 reset();
01607 }
01608
01609
01610
01611 double Efficiency::inConeCorr( uint32_t eta_bin, uint32_t pt_bin ) const {
01612 if ( check(eta_bin,pt_bin,__func__) ) {
01613 return ( outOfConeCorr( eta_bin, pt_bin ) *
01614 leakage_.value( eta_bin, pt_bin ) *
01615 response_.value( eta_bin, pt_bin ) );
01616 } else { return 0.; }
01617 }
01618
01619
01620
01621 double Efficiency::outOfConeCorr( uint32_t eta_bin, uint32_t pt_bin ) const {
01622 if ( check(eta_bin,pt_bin,__func__) ) {
01623 uint16_t ntrks = nTrks( eta_bin, pt_bin );
01624 double mean = meanE( eta_bin, pt_bin );
01625 double eff = ( 1. - efficiency_.value( eta_bin, pt_bin ) ) / efficiency_.value( eta_bin, pt_bin );
01626 if ( !ntrks ) { return 0.; }
01627 return ( ntrks * eff * mean );
01628 } else { return 0.; }
01629 }
01630
01631
01632
01633 uint16_t Efficiency::nTrks( uint32_t eta_bin, uint32_t pt_bin ) const {
01634 if ( check(eta_bin,pt_bin,__func__) ) {
01635 return data_[eta_bin][pt_bin].first;
01636 } else { return 0; }
01637 }
01638
01639
01640
01641 double Efficiency::sumE( uint32_t eta_bin, uint32_t pt_bin ) const {
01642 if ( check(eta_bin,pt_bin,__func__) ) {
01643 return data_[eta_bin][pt_bin].second;
01644 } else { return 0.; }
01645 }
01646
01647
01648
01649 double Efficiency::meanE( uint32_t eta_bin, uint32_t pt_bin ) const {
01650 if ( check(eta_bin,pt_bin,__func__) ) {
01651 Pair tmp = data_[eta_bin][pt_bin];
01652 if ( tmp.first ) { return tmp.second / tmp.first; }
01653 else { return 0.; }
01654 } else { return 0.; }
01655 }
01656
01657
01658
01659 void Efficiency::addE( uint32_t eta_bin, uint32_t pt_bin, double energy ) {
01660 if ( check(eta_bin,pt_bin,__func__) ) {
01661 data_[eta_bin][pt_bin].first++;
01662 data_[eta_bin][pt_bin].second += energy;
01663 }
01664 }
01665
01666
01667
01668 bool Efficiency::check( uint32_t eta_bin, uint32_t pt_bin, std::string method ) const {
01669 if ( eta_bin < data_.size() && pt_bin < ( data_.empty() ? 0 : data_[0].size() ) ) { return true; }
01670 else {
01671
01672
01673
01674
01675
01676 return false;
01677 }
01678 }
01679
01680
01681
01682 void Efficiency::reset() {
01683 data_.clear();
01684 data_.resize( response_.nEtaBins(), VPair( response_.nPtBins(), Pair(0,0.) ) );
01685 }
01686
01687
01688
01689 void Efficiency::print() const {
01690 if ( edm::isDebugEnabled() ) {
01691 std::stringstream ss;
01692 ss << "[jpt::Efficiency::" << __func__ << "]"
01693 << " Contents of maps:" << std::endl;
01694 ss << "Response map: " << std::endl;
01695 response_.print(ss);
01696 ss << "Efficiency map: " << std::endl;
01697 efficiency_.print(ss);
01698 ss << "Leakage map: " << std::endl;
01699 leakage_.print(ss);
01700 LogTrace("JetPlusTrackCorrector") << ss.str();
01701 }
01702 }