CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoJets/JetPlusTracks/src/JetPlusTrackCorrector.cc

Go to the documentation of this file.
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 //  std::cout<<" JetPlusTrackCorrector::correction "<<std::endl;
00108 
00109   theResponseOfChargedWithEff = 0.;
00110   theResponseOfChargedWithoutEff = 0.;
00111   theSumPtWithEff = 0.;
00112   theSumPtWithoutEff = 0.;
00113   theSumEnergyWithEff = 0.;
00114   theSumEnergyWithoutEff = 0.;
00115   
00116   // Corrected 4-momentum for jet
00117   corrected = fJet.p4();
00118   
00119   // Check if jet can be JPT-corrected
00120   if ( !canCorrect(fJet) ) { return 1.; }
00121 
00122 //  std::cout<<" Jet can be corrected "<<std::endl;
00123   
00124   // Match tracks to different particle types
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  // std::cout<<" Tracks are matched "<<std::endl;
00132   
00133   // Debug
00134   if ( verbose_ ) {
00135     edm::LogInfo("JetPlusTrackCorrector")
00136       << "[JetPlusTrackCorrector::" << __func__ << "]"
00137       << " Applying JPT corrections...";
00138   }
00139 
00140   // Pion corrections (both scalar and vectorial)
00141   if ( usePions_ ) { corrected += pionCorrection( fJet.p4(), pions ); }
00142   
00143   // Muon corrections (both scalar and vectorial)
00144   if ( useMuons_ ) { corrected += muonCorrection( fJet.p4(), muons ); }
00145   
00146   // Electrons corrections (both scalar and vectorial)
00147   if ( useElecs_ ) { corrected += elecCorrection( fJet.p4(), elecs ); }
00148 
00149   // Define jet direction using total 3-momentum of tracks (overrides above)
00150   if ( vectorial_ && !vecResponse_ ) { corrected = jetDirFromTracks( corrected, pions, muons, elecs ); }
00151   
00152   // Check if corrected 4-momentum gives negative scale
00153   double scale = checkScale( fJet.p4(), corrected );
00154   
00155   // Debug
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. );// << std::endl
00186     edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
00187   }
00188   
00189 //   std::cout << " mScale= " << scale
00190 //          << " NewResponse " << corrected.energy() 
00191 //          << " Jet energy " << fJet.energy()
00192 //          << " event " << event.id().event() << std::endl;
00193 
00194 /*
00195   std::cout<<" Correction::Corrections for jet "<<  theResponseOfChargedWithEff <<" "<<
00196   theResponseOfChargedWithoutEff <<" "<<
00197   theSumPtWithEff <<" "<<
00198   theSumPtWithoutEff <<" "<<
00199   theSumEnergyWithEff <<" "<<
00200   theSumEnergyWithoutEff <<std::endl;
00201 */
00202 
00203 
00204   // Return energy correction
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, //@@ required by method in derived class
00230                                          jpt::MatchedTracks& pions, 
00231                                          jpt::MatchedTracks& muons, 
00232                                          jpt::MatchedTracks& elecs ) const {
00233   
00234   // Associate tracks to jet at both the Vertex and CaloFace
00235   JetTracks jet_tracks;
00236   bool ok = jetTrackAssociation( fJet, event, setup, jet_tracks ); 
00237  
00238   // std::cout<<"JetPlusTrackCorrector::matchTracks, JetTrackAssociation ok="<<ok<<std::endl;
00239   
00240   if ( !ok ) { return false; }
00241   
00242   // Track collections propagated to Vertex and CaloFace for "pions", muons and electrons
00243   matchTracks( jet_tracks, event, pions, muons, elecs );
00244 
00245   // Debug
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   // Some init
00266   trks.clear();
00267   
00268   // Check if labels are given
00269   if ( !jetTracksAtVertex_.label().empty() && !jetTracksAtCalo_.label().empty() ) { 
00270 
00271     //std::cout<<" Try to associate with splitting "<< std::endl;  
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   // Get Jet-track association at Vertex
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   // std::cout<<" JetPlusTrackCorrector::jtaUsingEventData::here "<<jetSplitMerge_<<std::endl;
00319     
00320   // Retrieve jet-tracks association for given jet
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   // Check if any tracks are associated to jet at vertex
00327   if ( trks.vertex_.empty() ) { return false; }
00328 
00329   // Get Jet-track association at Calo
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   // Retrieve jet-tracks association for given jet
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 //   std::cout<<" JTA:Tracks in vertex "<<trks.vertex_.size()<<" in calo "<<trks.caloFace_.size()<<std::endl;
00352  
00353   // Successful
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   // Some init  
00384   pions.clear(); 
00385   muons.clear(); 
00386   elecs.clear(); 
00387 
00388   // Need vertex for track cleaning
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   // Get RECO muons
00396   edm::Handle<RecoMuons> reco_muons;
00397   bool found_reco_muons = true;
00398   if ( useMuons_ ) { getMuons( event, reco_muons ); }
00399   
00400   // Get RECO electrons and their ids
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   // Check RECO products found
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   // Identify pions/muons/electrons that are "in/in" and "in/out"
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 //     std::cout<<" Track accepted "<<std::endl;
00424       
00425       TrackRefs::iterator it = jet_tracks.caloFace_.end();
00426       bool found = findTrack( jet_tracks, itrk, it );
00427 
00428 //       std::cout<<"JetPlusTrackCorrector::matchTracks:  jet_tracks.caloFace_.size()="<<jet_tracks.caloFace_.size()
00429 //                <<" found = "<<found<<std::endl;
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 //         std::cout<<"JetPlusTrackCorrector::matchTracks: pions.inVertexOutOfCalo_.push_back(*itrk), size="
00444 //                  <<pions.inVertexOutOfCalo_.size() <<std::endl; 
00445                            }
00446       } 
00447     } 
00448   }
00449   
00450   // Identify pions/muons/electrons that are "out/in"
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 //  if ( useTrackQuality_ && !(*itrk)->quality(trackQuality_) ) { return true; }
00529 //  else { return false; }
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   std::cout<<" pionCorrection::Corrections for jet "<<  theResponseOfChargedWithEff <<" "<<
00598   theResponseOfChargedWithoutEff <<" "<<
00599   theSumPtWithEff <<" "<<
00600   theSumPtWithoutEff <<" "<<
00601   theSumEnergyWithEff <<" "<<
00602   theSumEnergyWithoutEff <<std::endl;
00603 */
00604 
00605   // In-cone
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   // Out-of-cone
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   // Out-of-vertex
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; //@@ null 4-momentum
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; //@@ to be implemented
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   // Correction to be applied to jet 4-momentum
00730   P4 corr;
00731   
00732   bool tracks_present = false;
00733   
00734   // Correct using pions in-cone at vertex
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   // Correct using muons in-cone at vertex
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   // Correct using electrons in-cone at vertex
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   // Adjust direction if tracks are present
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   // Correction to be applied to jet 4-momentum
00818   P4 correction;
00819 
00820 /*
00821   std::cout<<" >>>>> Jet "<<jet.pt()<<" "<<jet.eta()<<" "<<jet.phi()<<" Number of tracks "<<tracks.size()<<" Pions? "<<is_pion<<" InVert "<<in_cone_at_vertex<<
00822   " InCalo "<<in_cone_at_calo_face<<std::endl;
00823   
00824   std::cout<<" calculateCorr Start::Corrections for jet "<<  theResponseOfChargedWithEff <<" "<<
00825   theResponseOfChargedWithoutEff <<" "<<
00826   theSumPtWithEff <<" "<<
00827   theSumPtWithoutEff <<" "<<
00828   theSumEnergyWithEff <<" "<<
00829   theSumEnergyWithoutEff <<std::endl;
00830 */
00831 
00832 
00833   // Reset efficiency container
00834   eff.reset();
00835   
00836   double theSumResp = 0;
00837   double theSumPt = 0;
00838   double theSumEnergy = 0;
00839   
00840   
00841   // Iterate through tracks
00842   if ( !tracks.empty() ) {
00843     TrackRefs::iterator itrk = tracks.begin();
00844     TrackRefs::iterator jtrk = tracks.end();
00845     
00846     for ( ; itrk != jtrk; ++itrk ) {
00847 
00848       // Ignore high-pt tracks (only when in-cone and not a mip)
00849        if ( in_cone_at_calo_face && is_pion && (*itrk)->pt() >= 50. ) { continue; }
00850  
00851       // Inner track 4-momentum
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       // Add track momentum (if in-cone at vertex)
00864       if ( in_cone_at_vertex ) { correction += inner; }
00865       
00866       // Find appropriate eta/pt bin for given track
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       // Check bins (not for mips)
00873       if ( is_pion && ( ieta == response_.nEtaBins() || ipt == response_.nPtBins() ) ) { continue; }
00874       
00875       // Outer track 4-momentum 
00876       P4 outer;
00877       if ( in_cone_at_calo_face ) { 
00878         if ( vectorial_ && vecResponse_ ) {
00879           // Build 4-momentum from outer track (SHOULD USE IMPACT POINT?!)
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(); //@@ outerMomentum().eta()
00886             outer_phi = (*itrk)->outerPosition().phi(); //@@ outerMomentum().phi()
00887           }
00888           outer = PtEtaPhiM( outer_pt, outer_eta, outer_phi, mass );
00889           // Check if mip or not
00890           if ( !is_pion ) { outer *= ( outer.energy() > 0. ? mip / outer.energy() : 1. ); } //@@ Scale to mip energy
00891           else { outer *= ( outer.energy() > 0. ? inner.energy() / outer.energy() : 1. ); } //@@ Scale to inner track energy
00892         } else {
00893           // Check if mip or not
00894           if ( !is_pion ) { outer = ( jet.energy() > 0. ? mip / jet.energy() : 1. ) * jet; } //@@ Jet 4-mom scaled by mip energy
00895           else { outer = inner; } //@@ Set to inner track 4-momentum
00896         }      
00897         if ( is_pion ) { outer *= response_.value(ieta,ipt); } //@@ Scale by pion response
00898         correction -= outer; //@@ Subtract 
00899         
00900 // Calculate the sum of responses       
00901         theSumResp += response_.value(ieta,ipt);
00902       }
00903       
00904 // Calculate the sum of pt and energies 
00905         theSumPt += inner.pt();
00906         theSumEnergy += inner.energy(); 
00907       
00908       // Record inner track energy for pion efficiency correction
00909       if ( is_pion ) { eff.addE( ieta, ipt, inner.energy() ); }
00910         
00911       // Debug
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     } // loop through tracks
00927   } // ntracks != 0
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   std::cout<<" calculateCorr End::Corrections for jet "<<  theResponseOfChargedWithEff <<" "<<
00941   theResponseOfChargedWithoutEff <<" "<<
00942   theSumPtWithEff <<" "<<
00943   theSumPtWithoutEff <<" "<<
00944   theSumEnergyWithEff <<" "<<
00945   theSumEnergyWithoutEff <<std::endl;
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   // Total correction to be applied
00958   P4 correction;
00959   
00960   double theSumResp = 0;
00961   double theSumPt = 0;
00962   double theSumEnergy = 0;
00963  
00964   // Iterate through eta/pt bins
00965   for ( uint32_t ieta = 0; ieta < response_.nEtaBins()-1; ++ieta ) {
00966     for ( uint32_t ipt = 0; ipt < response_.nPtBins()-1; ++ipt ) {
00967 
00968       // Check tracks are found in this eta/pt bin
00969       if ( !eff.nTrks(ieta,ipt) ) { continue; }
00970 
00971       for ( uint16_t ii = 0; ii < 2; ++ii  ) {
00972         
00973         // Check which correction should be applied
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         // Calculate correction to be applied   
00980         P4 corr_p4;
00981         if ( vectorial_ && vecResponse_ ) {
00982           double corr_eta = response_.binCenterEta(ieta);
00983           double corr_phi = jet.phi(); //@@ jet phi!
00984           double corr_pt  = response_.binCenterPt(ipt);
00985           corr_p4 = PtEtaPhiM( corr_pt, corr_eta, corr_phi, pionMass_ ); //@@ E^2 = p^2 + m_pion^2, |p| calc'ed from pt bin
00986           corr_p4 *= ( corr_p4.energy() > 0. ? corr / corr_p4.energy() : 1. ); //@@ p4 scaled up by mean energy for bin
00987         } else { 
00988           corr_p4 = ( jet.energy() > 0. ? corr / jet.energy() : 1. ) * jet;
00989         }      
00990 
00991         // Apply correction
00992         if ( ii == 0 )      { correction += corr_p4; theSumPt += corr_p4.pt(); theSumEnergy += corr_p4.energy();} //@@ Add out-of-cone
00993         else if ( ii == 1 ) { correction -= corr_p4; theSumResp += corr_p4.energy();} //@@ Subtract in-cone
00994 
00995       } 
00996 
00997     }
00998   }
00999 
01000       theResponseOfChargedWithEff += theSumResp;
01001       theSumPtWithEff += theSumPt;
01002       theSumEnergyWithEff += theSumEnergy;
01003 /*      
01004       std::cout<<" Efficiency correction End::Corrections for jet "<<  theResponseOfChargedWithEff <<" "<<
01005       theResponseOfChargedWithoutEff <<" "<<
01006       theSumPtWithEff <<" "<<
01007       theSumPtWithoutEff <<" "<<
01008       theSumEnergyWithEff <<" "<<
01009       theSumEnergyWithoutEff <<std::endl;
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; } //@@ Check for null value 
01067     
01068     // DR matching b/w electron and track
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   //std::cout<<" NEW1 Merge/Split schema "<<jetSplitMerge_<<std::endl;
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   //std::cout<<" Size of initial vector "<<tracks.size()<<" "<<fJet.et()<<" "<<fJet.eta()<<" "<<fJet.phi()<<std::endl;
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       //       double dfi = fabs(fJet.phi()-(**it).phi());
01111       //       if(dfi>4.*atan(1.))dfi = 8.*atan(1.)-dfi;
01112       //       double deta = fJet.eta() - (**it).eta();
01113       //       double dR2check = sqrt(dfi*dfi+deta*deta);
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             //std::cout<<" Track belong to another jet also "<<dR2<<" "<<
01131             //(*ii)->et()<<" "<<(*ii)->eta()<<" "<< (*ii)->phi()<<" Track "<<(**it).eta()<<" "<<(**it).phi()<<" "<<scalethis<<" "<<scale<<" "<<flag<<std::endl;
01132             break;
01133           }
01134         }
01135 
01136       //std::cout<<" Track "<<tr<<" "<<flag<<" "<<dR2this<<" "<<dR2check<<" Jet "<<fJet.eta()<<" "<< fJet.phi()<<" Track "<<(**it).eta()<<" "<<(**it).phi()<<std::endl;
01137       if(flag == 1) {tracksthis.push_back (*it);}else{Excl.push_back (*it);}
01138     }
01139 
01140   //std::cout<<" The new size of tracks "<<tracksthis.size()<<" Excludede "<<Excl.size()<<std::endl;
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   //std::cout<<" NEW2" << std::endl;
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   //std::cout<<" Size of initial vector "<<tracks.size()<<" "<<fJet.et()<<" "<<fJet.eta()<<" "<<fJet.phi()<<std::endl;
01162 
01163   for(TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++ )
01164     {
01165 
01166       //std::cout<<" Track at calo surface "
01167       //<<" Track "<<(**it).eta()<<" "<<(**it).phi()<<std::endl;
01168       TrackRefs::iterator itold = find(Excl.begin(),Excl.end(),(*it));
01169       if(itold == Excl.end()) {
01170         tracksthis.push_back (*it);
01171       } 
01172       //else { std::cout<<"Exclude "<<(**it).eta()<<" "<<(**it).phi()<<std::endl; }
01173 
01174     }
01175 
01176   //std::cout<<" Size of calo tracks "<<tracksthis.size()<<std::endl;
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 //std::cout<<"---JetPlusTrackCorrector::correctAA, Jet initial: NewResponse="<<NewResponse <<" et = "<<fJet.et()
01194 //         <<" pt= "<<fJet.pt()<<" eta="<<fJet.eta()<<" phi="<<fJet.phi() <<" ja="<<ja
01195 //        <<" trBgOutOfVertex="<<trBgOutOfVertex.size()<< " trBgOutOfCalo="<<trBgOutOfCalo.size()<<std::endl;
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 // calculate the mean response
01205 //
01206 //     double MeanResponseOfBackgroundCharged = 0.;
01207 //     double MeanEnergyOfBackgroundCharged   = 0.;
01208 
01209 //================= EnergyOfBackgroundCharged ==================>   
01210      for( reco::TrackRefVector::iterator iBgtV = trBgOutOfVertex.begin(); iBgtV != trBgOutOfVertex.end(); iBgtV++)
01211        {
01212        // Temporary solution>>>>>> Remove tracks with pt>50 GeV
01213        //   if( (**iBgtV).pt() >= 50. ) continue;
01214                //response_.value(ieta,ipt);
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 //      // Check bins (not for mips)
01223 //         if ( ieta >= response_.nEtaBins() ) { continue; }
01224 //       if ( ipt >= response_.nPtBins() ) { ipt = response_.nPtBins() - 1; }
01225          
01226          double echarBg=sqrt((**iBgtV).px()*(**iBgtV).px()+(**iBgtV).py()*(**iBgtV).py()+(**iBgtV).pz()*(**iBgtV).pz()+0.14*0.14);
01227          
01228 //       ResponseOfBackgroundCharged += echarBg*response_.value(ieta,ipt)/efficiency_.value(ieta,ipt);
01229          
01230 //       std::cout<<" Efficiency of bg tracks "<<efficiency_.value(ieta,ipt)<<std::endl;
01231          
01232          EnergyOfBackgroundCharged += echarBg/efficiency_.value(ieta,ipt);
01233 
01234          NumberOfBackgroundChargedVertex++;
01235          
01236        } // Energy BG tracks
01237 
01238 //        std::cout<<" JetPlusTrackCorrector.cc, NumberOfBackgroundChargedVertex ="<<NumberOfBackgroundChargedVertex<<std::endl;
01239 
01240 //============= ResponseOfBackgroundCharged =======================>
01241 
01242      for( reco::TrackRefVector::iterator iBgtC = trBgOutOfCalo.begin(); iBgtC != trBgOutOfCalo.end(); iBgtC++)
01243        {
01244        // Temporary solution>>>>>> Remove tracks with pt>50 GeV
01245        //   if( (**iBgtC).pt() >= 50. ) continue;
01246                //response_.value(ieta,ipt);
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       // Check bins (not for mips)
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 //       std::cout<<" Efficiency of bg tracks "<<efficiency_.value(ieta,ipt)<<std::endl;
01263 
01264 
01265          NumberOfBackgroundChargedCalo++;
01266 
01267        } // Response of BG tracks
01268 
01269 //  std::cout<<" JetPlusTrackCorrector.cc, NumberOfBackgroundChargedCalo ="<<NumberOfBackgroundChargedCalo<<std::endl;
01270 //=================================================================>
01271 
01272 //=================================================================>
01273 // Look for in-out tracks
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 //          std::cout<<" Track in out, size= "<<pioninout.size()<<" p= "<<(*it)->p()<<" pt= "<<(*it)->pt()
01283 //                   <<" eta=  "<<(*it)->eta()<<" phi= "<<(*it)->phi()<<std::endl;
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 // Look for in-in tracks
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 //          std::cout<<" Track in in, size= "<<pioninin.size()<<" p= "<<(*it)->p()<<" pt= "<<(*it)->pt()
01301 //                   <<" eta= "<<(*it)->eta()<<" phi= "<<(*it)->phi()<<std::endl;
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 //    double SquareEtaRingWithoutJets = 4*(M_PI*mConeSize - mConeSize*mConeSize);
01314     double SquareEtaRingWithoutJets = ja;
01315     
01316 //   std::cout<<"---SquareEtaRingWithoutJets="<<SquareEtaRingWithoutJets<<std::endl;
01317        
01318     EnergyOfBackgroundCharged       = EnergyOfBackgroundCharged/SquareEtaRingWithoutJets;
01319     ResponseOfBackgroundCharged     = ResponseOfBackgroundCharged/SquareEtaRingWithoutJets;
01320     NumberOfBackgroundChargedVertex = NumberOfBackgroundChargedVertex/SquareEtaRingWithoutJets;
01321     NumberOfBackgroundChargedCalo   = NumberOfBackgroundChargedCalo/SquareEtaRingWithoutJets;
01322 //    NumberOfBackgroundCharged   = NumberOfBackgroundCharged/SquareEtaRingWithoutJets;
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 //    NumberOfBackgroundCharged   = M_PI*mConeSize*mConeSize*NumberOfBackgroundCharged;
01330 
01331 
01332     NewResponse = NewResponse - EnergyOfBackgroundCharged + ResponseOfBackgroundCharged;
01333     
01334 //      std::cout<<"====JetPlusTrackCorrector, Old response=fJet.energy()"<<fJet.energy()
01335 //               <<" EnergyOfBackgroundCharged="<<EnergyOfBackgroundCharged
01336 //               <<" ResponseOfBackgroundCharged="<<ResponseOfBackgroundCharged
01337 //               <<" NewResponse="<<NewResponse<<" NumberOfBackgroundChargedVertex="<<NumberOfBackgroundChargedVertex
01338 //               <<" NumberOfBackgroundChargedCalo="<<NumberOfBackgroundChargedCalo<<std::endl;
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   // Some init
01353   clear();
01354   std::vector<Element> data;
01355 
01356   // Parse file
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   // Populate container
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   // Check
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   // Check
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   // Debug
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 //    edm::LogWarning("JetPlusTrackCorrector") 
01441 //      << "[jpt::Map::" << __func__ << "]"
01442 //      << " Trying to access element " << eta_bin
01443 //      << " of a vector with size " << eta_.size()
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 //    edm::LogWarning("JetPlusTrackCorrector") 
01455 //      << "[jpt::Map::" << __func__ << "]"
01456 //      << " Trying to access element " << pt_bin
01457 //      << " of a vector with size " << pt_.size()
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 //    edm::LogWarning("JetPlusTrackCorrector") 
01470 //      << "[jpt::Map::" << __func__ << "]"
01471 //      << " Trying to access element " << eta_bin+1
01472 //      << " of a vector with size " << eta_.size()
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 //    edm::LogWarning("JetPlusTrackCorrector") 
01485 //      << "[jpt::Map::" << __func__ << "]"
01486 //      << " Trying to access element " << pt_bin+1
01487 //      << " of a vector with size " << pt_.size()
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 ) { //@@ "-1" is bug?
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 ) { //@@ "-1" is bug?
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 //    edm::LogWarning("JetPlusTrackCorrector") 
01520 //      << "[jpt::Map::" << __func__ << "]"
01521 //      << " Trying to access element (" << eta_bin << "," << pt_bin << ")"
01522 //      << " of a vector with size (" << data_.size() << "," << ( data_.empty() ? 0 : data_[0].size() ) << ")"
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 //    edm::LogWarning("JetPlusTrackCorrector") 
01672 //      << "[jpt::Efficiency::" << method << "]"
01673 //      << " Trying to access element (" << eta_bin << "," << pt_bin << ")"
01674 //      << " of a vector with size (" << data_.size() << "," << ( data_.empty() ? 0 : data_[0].size() ) << ")"
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 }