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.2 2012/09/20 15:43:30 veelken 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 beginJob() ;
00057         virtual void produce(edm::Event&, const edm::EventSetup&);
00058         virtual void endJob() ;
00059       
00060         virtual void beginRun(edm::Run&, edm::EventSetup const&);
00061         virtual void endRun(edm::Run&, edm::EventSetup const&);
00062         virtual void beginLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const&);
00063         virtual void endLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const&);
00064 
00065         void initJetEnergyCorrector(const edm::EventSetup &iSetup, bool isData);
00066 
00067         edm::InputTag jets_, vertexes_, jetids_, rho_;
00068         std::string jec_;
00069         bool runMvas_, produceJetIds_, inputIsCorrected_, applyJec_;
00070         std::vector<std::pair<std::string, PileupJetIdAlgo *> > algos_;
00071         
00072         bool residualsFromTxt_;
00073         edm::FileInPath residualsTxt_;
00074         FactorizedJetCorrector *jecCor_;
00075         std::vector<JetCorrectorParameters> jetCorPars_;
00076 };
00077 
00078 // ------------------------------------------------------------------------------------------
00079 PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig)
00080 {
00081         runMvas_ = iConfig.getParameter<bool>("runMvas");
00082         produceJetIds_ = iConfig.getParameter<bool>("produceJetIds");
00083         jets_ = iConfig.getParameter<edm::InputTag>("jets");
00084         vertexes_ = iConfig.getParameter<edm::InputTag>("vertexes");
00085         jetids_  = iConfig.getParameter<edm::InputTag>("jetids");
00086         inputIsCorrected_ = iConfig.getParameter<bool>("inputIsCorrected");
00087         applyJec_ = iConfig.getParameter<bool>("applyJec");
00088         jec_ =  iConfig.getParameter<std::string>("jec");
00089         rho_ = iConfig.getParameter<edm::InputTag>("rho");
00090         residualsFromTxt_ = iConfig.getParameter<bool>("residualsFromTxt");
00091         residualsTxt_ = iConfig.getParameter<edm::FileInPath>("residualsTxt");
00092         std::vector<edm::ParameterSet> algos = iConfig.getParameter<std::vector<edm::ParameterSet> >("algos");
00093         
00094         jecCor_ = 0;
00095 
00096         if( ! runMvas_ ) assert( algos.size() == 1 );
00097         
00098         if( produceJetIds_ ) {
00099                 produces<edm::ValueMap<StoredPileupJetIdentifier> > ("");
00100         }
00101         for(std::vector<edm::ParameterSet>::iterator it=algos.begin(); it!=algos.end(); ++it) {
00102                 std::string label = it->getParameter<std::string>("label");
00103                 algos_.push_back( std::make_pair(label,new PileupJetIdAlgo(*it)) );
00104                 if( runMvas_ ) {
00105                         produces<edm::ValueMap<float> > (label+"Discriminant");
00106                         produces<edm::ValueMap<int> > (label+"Id");
00107                 }
00108         }
00109 }
00110 
00111 
00112 
00113 // ------------------------------------------------------------------------------------------
00114 PileupJetIdProducer::~PileupJetIdProducer()
00115 {
00116 }
00117 
00118 
00119 // ------------------------------------------------------------------------------------------
00120 void
00121 PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00122 {
00123         using namespace edm;
00124         using namespace std;
00125         using namespace reco;
00126         
00127         // Input jets
00128         Handle<View<Jet> > jetHandle;
00129         iEvent.getByLabel(jets_,jetHandle);
00130         const View<Jet> & jets = *jetHandle;
00131         // vertexes 
00132         Handle<VertexCollection> vertexHandle;
00133         if(  produceJetIds_ ) {
00134                 iEvent.getByLabel(vertexes_, vertexHandle);
00135         }
00136         const VertexCollection & vertexes = *(vertexHandle.product());
00137         // input variables
00138         Handle<ValueMap<StoredPileupJetIdentifier> > vmap;
00139         if( ! produceJetIds_ ) {
00140                 iEvent.getByLabel(jetids_, vmap);
00141         }
00142         // rho
00143         edm::Handle< double > rhoH;
00144         double rho = 0.;
00145         
00146         // products
00147         vector<StoredPileupJetIdentifier> ids; 
00148         map<string, vector<float> > mvas;
00149         map<string, vector<int> > idflags;
00150 
00151         VertexCollection::const_iterator vtx;
00152         if( produceJetIds_ ) {
00153                 // require basic quality cuts on the vertexes
00154                 vtx = vertexes.begin();
00155                 while( vtx != vertexes.end() && ( vtx->isFake() || vtx->ndof() < 4 ) ) {
00156                         ++vtx;
00157                 }
00158                 if( vtx == vertexes.end() ) { vtx = vertexes.begin(); }
00159         }
00160         
00161         // Loop over input jets
00162         for ( unsigned int i=0; i<jets.size(); ++i ) {
00163                 // Pick the first algo to compute the input variables
00164                 vector<pair<string,PileupJetIdAlgo *> >::iterator algoi = algos_.begin();
00165                 PileupJetIdAlgo * ialgo = algoi->second;
00166                 
00167                 const Jet & jet = jets.at(i);
00168                 //const pat::Jet * patjet =  dynamic_cast<const pat::Jet *>(&jet);
00169                 //bool ispat = patjet != 0;
00170                 
00171                 // Get jet energy correction
00172                 float jec = 0.;
00173                 if( applyJec_ ) {
00174                         // If haven't done it get rho from the event
00175                         if( rho == 0. ) {
00176                                 iEvent.getByLabel(rho_,rhoH);
00177                                 rho = *rhoH;
00178                         }
00179                         // jet corrector
00180                         if( jecCor_ == 0 ) {
00181                                 initJetEnergyCorrector( iSetup, iEvent.isRealData() );
00182                         }
00183                         //if( ispat ) {
00184                         //      jecCor_->setJetPt(patjet->correctedJet(0).pt());
00185                         //} else {
00186                         jecCor_->setJetPt(jet.pt());
00187                         //}
00188                         jecCor_->setJetEta(jet.eta());
00189                         jecCor_->setJetA(jet.jetArea());
00190                         jecCor_->setRho(rho);
00191                         jec = jecCor_->getCorrection();
00192                 }
00193                 
00194                 // If it was requested or the input is an uncorrected jet apply the JEC
00195                 bool applyJec = applyJec_ || !inputIsCorrected_;  //( ! ispat && ! inputIsCorrected_ );
00196                 reco::Jet * corrJet = 0;
00197                 if( applyJec ) {
00198                         float scale = jec;
00199                         //if( ispat ) {
00200                         //      corrJet = new pat::Jet(patjet->correctedJet(0)) ;
00201                         //} else {
00202                         corrJet = dynamic_cast<reco::Jet *>( jet.clone() );
00203                         //}
00204                         corrJet->scaleEnergy(scale);
00205                 }
00206                 const reco::Jet * theJet = ( applyJec ? corrJet : &jet );
00207                 
00208                 PileupJetIdentifier puIdentifier;
00209                 if( produceJetIds_ ) {
00210                         // Compute the input variables
00211                         puIdentifier = ialgo->computeIdVariables(theJet, jec,  &(*vtx), vertexes, runMvas_);
00212                         ids.push_back( puIdentifier );
00213                 } else {
00214                         // Or read it from the value map
00215                         puIdentifier = (*vmap)[jets.refAt(i)]; 
00216                         puIdentifier.jetPt(theJet->pt());    // make sure JEC is applied when computing the MVA
00217                         puIdentifier.jetEta(theJet->eta());
00218                         puIdentifier.jetPhi(theJet->phi());
00219                         ialgo->set(puIdentifier); 
00220                         puIdentifier = ialgo->computeMva();
00221                 }
00222                 
00223                 if( runMvas_ ) {
00224                         // Compute the MVA and WP
00225                         mvas[algoi->first].push_back( puIdentifier.mva() );
00226                         idflags[algoi->first].push_back( puIdentifier.idFlag() );
00227                         for( ++algoi; algoi!=algos_.end(); ++algoi) {
00228                                 ialgo = algoi->second;
00229                                 ialgo->set(puIdentifier);
00230                                 PileupJetIdentifier id = ialgo->computeMva();
00231                                 mvas[algoi->first].push_back( id.mva() );
00232                                 idflags[algoi->first].push_back( id.idFlag() );
00233                         }
00234                 }
00235                 
00236                 // cleanup
00237                 if( corrJet ) { delete corrJet; }
00238         }
00239         
00240         // Produce the output value maps
00241         if( runMvas_ ) {
00242                 for(vector<pair<string,PileupJetIdAlgo *> >::iterator ialgo = algos_.begin(); ialgo!=algos_.end(); ++ialgo) {
00243                         // MVA
00244                         vector<float> & mva = mvas[ialgo->first];
00245                         auto_ptr<ValueMap<float> > mvaout(new ValueMap<float>());
00246                         ValueMap<float>::Filler mvafiller(*mvaout);
00247                         mvafiller.insert(jetHandle,mva.begin(),mva.end());
00248                         mvafiller.fill();
00249                         iEvent.put(mvaout,ialgo->first+"Discriminant");
00250                         
00251                         // WP
00252                         vector<int> & idflag = idflags[ialgo->first];
00253                         auto_ptr<ValueMap<int> > idflagout(new ValueMap<int>());
00254                         ValueMap<int>::Filler idflagfiller(*idflagout);
00255                         idflagfiller.insert(jetHandle,idflag.begin(),idflag.end());
00256                         idflagfiller.fill();
00257                         iEvent.put(idflagout,ialgo->first+"Id");
00258                 }
00259         }
00260         // input variables
00261         if( produceJetIds_ ) {
00262                 assert( jetHandle->size() == ids.size() );
00263                 auto_ptr<ValueMap<StoredPileupJetIdentifier> > idsout(new ValueMap<StoredPileupJetIdentifier>());
00264                 ValueMap<StoredPileupJetIdentifier>::Filler idsfiller(*idsout);
00265                 idsfiller.insert(jetHandle,ids.begin(),ids.end());
00266                 idsfiller.fill();
00267                 iEvent.put(idsout);
00268         }
00269 }
00270 
00271 // ------------------------------------------------------------------------------------------
00272 void 
00273 PileupJetIdProducer::beginJob()
00274 {
00275 }
00276 
00277 // ------------------------------------------------------------------------------------------
00278 void 
00279 PileupJetIdProducer::endJob() {
00280 }
00281 
00282 // ------------------------------------------------------------------------------------------
00283 void 
00284 PileupJetIdProducer::beginRun(edm::Run&, edm::EventSetup const&)
00285 {
00286 }
00287 
00288 // ------------------------------------------------------------------------------------------
00289 void 
00290 PileupJetIdProducer::endRun(edm::Run&, edm::EventSetup const&)
00291 {
00292 }
00293 
00294 // ------------------------------------------------------------------------------------------
00295 void 
00296 PileupJetIdProducer::beginLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const&)
00297 {
00298 }
00299 
00300 // ------------------------------------------------------------------------------------------
00301 void 
00302 PileupJetIdProducer::endLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const&)
00303 {
00304 }
00305 
00306 // ------------------------------------------------------------------------------------------
00307 void
00308 PileupJetIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
00309         //The following says we do not know what parameters are allowed so do no validation
00310         // Please change this to state exactly what you do use, even if it is no parameters
00311         edm::ParameterSetDescription desc;
00312         desc.setUnknown();
00313         descriptions.addDefault(desc);
00314 }
00315 
00316 
00317 // ------------------------------------------------------------------------------------------
00318 void 
00319 PileupJetIdProducer::initJetEnergyCorrector(const edm::EventSetup &iSetup, bool isData)
00320 {
00321         //jet energy correction levels to apply on raw jet
00322         std::vector<std::string> jecLevels;
00323         jecLevels.push_back("L1FastJet");
00324         jecLevels.push_back("L2Relative");
00325         jecLevels.push_back("L3Absolute");
00326         if(isData && ! residualsFromTxt_ ) jecLevels.push_back("L2L3Residual");
00327 
00328         //check the corrector parameters needed according to the correction levels
00329         edm::ESHandle<JetCorrectorParametersCollection> parameters;
00330         iSetup.get<JetCorrectionsRecord>().get(jec_,parameters);
00331         for(std::vector<std::string>::const_iterator ll = jecLevels.begin(); ll != jecLevels.end(); ++ll)
00332         { 
00333                 const JetCorrectorParameters& ip = (*parameters)[*ll];
00334                 jetCorPars_.push_back(ip); 
00335         } 
00336         if( isData && residualsFromTxt_ ) {
00337                 jetCorPars_.push_back(JetCorrectorParameters(residualsTxt_.fullPath())); 
00338         }
00339         
00340         //instantiate the jet corrector
00341         jecCor_ = new FactorizedJetCorrector(jetCorPars_);
00342 }
00343 //define this as a plug-in
00344 DEFINE_FWK_MODULE(PileupJetIdProducer);
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2_patch1/src/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc
