Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    PileupJetIdProducer
00004 // Class:      PileupJetIdProducer
00005 // 
00013 //
00014 // Original Author:  Pasquale Musella,40 2-A12,+41227671706,
00015 //         Created:  Wed Apr 18 15:48:47 CEST 2012
00016 // $Id: PileupJetIdProducer.cc,v 1.3 2013/02/27 20:57:37 eulisse Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 
00024 // user include files
00025 #include "FWCore/Framework/interface/Frameworkfwd.h"
00026 #include "FWCore/Framework/interface/EDProducer.h"
00027 #include "DataFormats/JetReco/interface/Jet.h"
00028 
00029 #include "FWCore/Framework/interface/Event.h"
00030 #include "FWCore/Framework/interface/MakerMacros.h"
00031 #include "DataFormats/Common/interface/ValueMap.h"
00032 #include "FWCore/Framework/interface/ESHandle.h"
00033 #include "FWCore/Framework/interface/EventSetup.h"
00034 
00035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00036 
00037 #include "DataFormats/JetReco/interface/PileupJetIdentifier.h"
00038 #include "RecoJets/JetProducers/interface/PileupJetIdAlgo.h"
00039 #include "DataFormats/VertexReco/interface/Vertex.h"
00040 
00041 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
00042 #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
00043 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
00044 
00045 #include "FWCore/ParameterSet/interface/FileInPath.h"
00046 
00047 // ------------------------------------------------------------------------------------------
00048 class PileupJetIdProducer : public edm::EDProducer {
00049 public:
00050         explicit PileupJetIdProducer(const edm::ParameterSet&);
00051         ~PileupJetIdProducer();
00052 
00053         static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
00054 
00055 private:
00056         virtual void produce(edm::Event&, const edm::EventSetup&);
00057       
00058 
00059         void initJetEnergyCorrector(const edm::EventSetup &iSetup, bool isData);
00060 
00061         edm::InputTag jets_, vertexes_, jetids_, rho_;
00062         std::string jec_;
00063         bool runMvas_, produceJetIds_, inputIsCorrected_, applyJec_;
00064         std::vector<std::pair<std::string, PileupJetIdAlgo *> > algos_;
00065         
00066         bool residualsFromTxt_;
00067         edm::FileInPath residualsTxt_;
00068         FactorizedJetCorrector *jecCor_;
00069         std::vector<JetCorrectorParameters> jetCorPars_;
00070 };
00071 
00072 // ------------------------------------------------------------------------------------------
00073 PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig)
00074 {
00075         runMvas_ = iConfig.getParameter<bool>("runMvas");
00076         produceJetIds_ = iConfig.getParameter<bool>("produceJetIds");
00077         jets_ = iConfig.getParameter<edm::InputTag>("jets");
00078         vertexes_ = iConfig.getParameter<edm::InputTag>("vertexes");
00079         jetids_  = iConfig.getParameter<edm::InputTag>("jetids");
00080         inputIsCorrected_ = iConfig.getParameter<bool>("inputIsCorrected");
00081         applyJec_ = iConfig.getParameter<bool>("applyJec");
00082         jec_ =  iConfig.getParameter<std::string>("jec");
00083         rho_ = iConfig.getParameter<edm::InputTag>("rho");
00084         residualsFromTxt_ = iConfig.getParameter<bool>("residualsFromTxt");
00085         residualsTxt_ = iConfig.getParameter<edm::FileInPath>("residualsTxt");
00086         std::vector<edm::ParameterSet> algos = iConfig.getParameter<std::vector<edm::ParameterSet> >("algos");
00087         
00088         jecCor_ = 0;
00089 
00090         if( ! runMvas_ ) assert( algos.size() == 1 );
00091         
00092         if( produceJetIds_ ) {
00093                 produces<edm::ValueMap<StoredPileupJetIdentifier> > ("");
00094         }
00095         for(std::vector<edm::ParameterSet>::iterator it=algos.begin(); it!=algos.end(); ++it) {
00096                 std::string label = it->getParameter<std::string>("label");
00097                 algos_.push_back( std::make_pair(label,new PileupJetIdAlgo(*it)) );
00098                 if( runMvas_ ) {
00099                         produces<edm::ValueMap<float> > (label+"Discriminant");
00100                         produces<edm::ValueMap<int> > (label+"Id");
00101                 }
00102         }
00103 }
00104 
00105 
00106 
00107 // ------------------------------------------------------------------------------------------
00108 PileupJetIdProducer::~PileupJetIdProducer()
00109 {
00110 }
00111 
00112 
00113 // ------------------------------------------------------------------------------------------
00114 void
00115 PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00116 {
00117         using namespace edm;
00118         using namespace std;
00119         using namespace reco;
00120         
00121         // Input jets
00122         Handle<View<Jet> > jetHandle;
00123         iEvent.getByLabel(jets_,jetHandle);
00124         const View<Jet> & jets = *jetHandle;
00125         // vertexes 
00126         Handle<VertexCollection> vertexHandle;
00127         if(  produceJetIds_ ) {
00128                 iEvent.getByLabel(vertexes_, vertexHandle);
00129         }
00130         const VertexCollection & vertexes = *(vertexHandle.product());
00131         // input variables
00132         Handle<ValueMap<StoredPileupJetIdentifier> > vmap;
00133         if( ! produceJetIds_ ) {
00134                 iEvent.getByLabel(jetids_, vmap);
00135         }
00136         // rho
00137         edm::Handle< double > rhoH;
00138         double rho = 0.;
00139         
00140         // products
00141         vector<StoredPileupJetIdentifier> ids; 
00142         map<string, vector<float> > mvas;
00143         map<string, vector<int> > idflags;
00144 
00145         VertexCollection::const_iterator vtx;
00146         if( produceJetIds_ ) {
00147                 // require basic quality cuts on the vertexes
00148                 vtx = vertexes.begin();
00149                 while( vtx != vertexes.end() && ( vtx->isFake() || vtx->ndof() < 4 ) ) {
00150                         ++vtx;
00151                 }
00152                 if( vtx == vertexes.end() ) { vtx = vertexes.begin(); }
00153         }
00154         
00155         // Loop over input jets
00156         for ( unsigned int i=0; i<jets.size(); ++i ) {
00157                 // Pick the first algo to compute the input variables
00158                 vector<pair<string,PileupJetIdAlgo *> >::iterator algoi = algos_.begin();
00159                 PileupJetIdAlgo * ialgo = algoi->second;
00160                 
00161                 const Jet & jet = jets.at(i);
00162                 //const pat::Jet * patjet =  dynamic_cast<const pat::Jet *>(&jet);
00163                 //bool ispat = patjet != 0;
00164                 
00165                 // Get jet energy correction
00166                 float jec = 0.;
00167                 if( applyJec_ ) {
00168                         // If haven't done it get rho from the event
00169                         if( rho == 0. ) {
00170                                 iEvent.getByLabel(rho_,rhoH);
00171                                 rho = *rhoH;
00172                         }
00173                         // jet corrector
00174                         if( jecCor_ == 0 ) {
00175                                 initJetEnergyCorrector( iSetup, iEvent.isRealData() );
00176                         }
00177                         //if( ispat ) {
00178                         //      jecCor_->setJetPt(patjet->correctedJet(0).pt());
00179                         //} else {
00180                         jecCor_->setJetPt(jet.pt());
00181                         //}
00182                         jecCor_->setJetEta(jet.eta());
00183                         jecCor_->setJetA(jet.jetArea());
00184                         jecCor_->setRho(rho);
00185                         jec = jecCor_->getCorrection();
00186                 }
00187                 
00188                 // If it was requested or the input is an uncorrected jet apply the JEC
00189                 bool applyJec = applyJec_ || !inputIsCorrected_;  //( ! ispat && ! inputIsCorrected_ );
00190                 reco::Jet * corrJet = 0;
00191                 if( applyJec ) {
00192                         float scale = jec;
00193                         //if( ispat ) {
00194                         //      corrJet = new pat::Jet(patjet->correctedJet(0)) ;
00195                         //} else {
00196                         corrJet = dynamic_cast<reco::Jet *>( jet.clone() );
00197                         //}
00198                         corrJet->scaleEnergy(scale);
00199                 }
00200                 const reco::Jet * theJet = ( applyJec ? corrJet : &jet );
00201                 
00202                 PileupJetIdentifier puIdentifier;
00203                 if( produceJetIds_ ) {
00204                         // Compute the input variables
00205                         puIdentifier = ialgo->computeIdVariables(theJet, jec,  &(*vtx), vertexes, runMvas_);
00206                         ids.push_back( puIdentifier );
00207                 } else {
00208                         // Or read it from the value map
00209                         puIdentifier = (*vmap)[jets.refAt(i)]; 
00210                         puIdentifier.jetPt(theJet->pt());    // make sure JEC is applied when computing the MVA
00211                         puIdentifier.jetEta(theJet->eta());
00212                         puIdentifier.jetPhi(theJet->phi());
00213                         ialgo->set(puIdentifier); 
00214                         puIdentifier = ialgo->computeMva();
00215                 }
00216                 
00217                 if( runMvas_ ) {
00218                         // Compute the MVA and WP
00219                         mvas[algoi->first].push_back( puIdentifier.mva() );
00220                         idflags[algoi->first].push_back( puIdentifier.idFlag() );
00221                         for( ++algoi; algoi!=algos_.end(); ++algoi) {
00222                                 ialgo = algoi->second;
00223                                 ialgo->set(puIdentifier);
00224                                 PileupJetIdentifier id = ialgo->computeMva();
00225                                 mvas[algoi->first].push_back( id.mva() );
00226                                 idflags[algoi->first].push_back( id.idFlag() );
00227                         }
00228                 }
00229                 
00230                 // cleanup
00231                 if( corrJet ) { delete corrJet; }
00232         }
00233         
00234         // Produce the output value maps
00235         if( runMvas_ ) {
00236                 for(vector<pair<string,PileupJetIdAlgo *> >::iterator ialgo = algos_.begin(); ialgo!=algos_.end(); ++ialgo) {
00237                         // MVA
00238                         vector<float> & mva = mvas[ialgo->first];
00239                         auto_ptr<ValueMap<float> > mvaout(new ValueMap<float>());
00240                         ValueMap<float>::Filler mvafiller(*mvaout);
00241                         mvafiller.insert(jetHandle,mva.begin(),mva.end());
00242                         mvafiller.fill();
00243                         iEvent.put(mvaout,ialgo->first+"Discriminant");
00244                         
00245                         // WP
00246                         vector<int> & idflag = idflags[ialgo->first];
00247                         auto_ptr<ValueMap<int> > idflagout(new ValueMap<int>());
00248                         ValueMap<int>::Filler idflagfiller(*idflagout);
00249                         idflagfiller.insert(jetHandle,idflag.begin(),idflag.end());
00250                         idflagfiller.fill();
00251                         iEvent.put(idflagout,ialgo->first+"Id");
00252                 }
00253         }
00254         // input variables
00255         if( produceJetIds_ ) {
00256                 assert( jetHandle->size() == ids.size() );
00257                 auto_ptr<ValueMap<StoredPileupJetIdentifier> > idsout(new ValueMap<StoredPileupJetIdentifier>());
00258                 ValueMap<StoredPileupJetIdentifier>::Filler idsfiller(*idsout);
00259                 idsfiller.insert(jetHandle,ids.begin(),ids.end());
00260                 idsfiller.fill();
00261                 iEvent.put(idsout);
00262         }
00263 }
00264 
00265 
00266 
00267 // ------------------------------------------------------------------------------------------
00268 void
00269 PileupJetIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
00270         //The following says we do not know what parameters are allowed so do no validation
00271         // Please change this to state exactly what you do use, even if it is no parameters
00272         edm::ParameterSetDescription desc;
00273         desc.setUnknown();
00274         descriptions.addDefault(desc);
00275 }
00276 
00277 
00278 // ------------------------------------------------------------------------------------------
00279 void 
00280 PileupJetIdProducer::initJetEnergyCorrector(const edm::EventSetup &iSetup, bool isData)
00281 {
00282         //jet energy correction levels to apply on raw jet
00283         std::vector<std::string> jecLevels;
00284         jecLevels.push_back("L1FastJet");
00285         jecLevels.push_back("L2Relative");
00286         jecLevels.push_back("L3Absolute");
00287         if(isData && ! residualsFromTxt_ ) jecLevels.push_back("L2L3Residual");
00288 
00289         //check the corrector parameters needed according to the correction levels
00290         edm::ESHandle<JetCorrectorParametersCollection> parameters;
00291         iSetup.get<JetCorrectionsRecord>().get(jec_,parameters);
00292         for(std::vector<std::string>::const_iterator ll = jecLevels.begin(); ll != jecLevels.end(); ++ll)
00293         { 
00294                 const JetCorrectorParameters& ip = (*parameters)[*ll];
00295                 jetCorPars_.push_back(ip); 
00296         } 
00297         if( isData && residualsFromTxt_ ) {
00298                 jetCorPars_.push_back(JetCorrectorParameters(residualsTxt_.fullPath())); 
00299         }
00300         
00301         //instantiate the jet corrector
00302         jecCor_ = new FactorizedJetCorrector(jetCorPars_);
00303 }
00304 //define this as a plug-in
00305 DEFINE_FWK_MODULE(PileupJetIdProducer);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc
