Reference Manual

Go to the documentation of this file.
00001 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
00002 #include "DataFormats/JetReco/interface/CaloJet.h"
00003 #include "DataFormats/BTauReco/interface/JetTag.h"
00004 #include "DataFormats/JetReco/interface/PFJet.h"
00005 #include "DQM/Physics/src/TopSingleLeptonDQM.h"
00006 #include "DataFormats/Math/interface/deltaR.h"
00007 #include <iostream>
00008 using namespace std;
00009 namespace TopSingleLepton {
00010 
00011   // maximal number of leading jets 
00012   // to be used for top mass estimate
00013   static const unsigned int MAXJETS = 4;
00014   // nominal mass of the W boson to 
00015   // be used for the top mass estimate
00016   static const double WMASS = 80.4;
00017 
00018   MonitorEnsemble::MonitorEnsemble(const char* label, const edm::ParameterSet& cfg) : 
00019     label_(label), elecIso_(0), elecSelect_(0), pvSelect_(0), muonIso_(0), muonSelect_(0), jetIDSelect_(0), includeBTag_(false), lowerEdge_(-1.), upperEdge_(-1.), logged_(0)
00020   {
00021     // sources have to be given; this PSet is not optional
00022     edm::ParameterSet sources=cfg.getParameter<edm::ParameterSet>("sources");
00023     muons_= sources.getParameter<edm::InputTag>("muons");
00024     elecs_= sources.getParameter<edm::InputTag>("elecs");
00025     jets_ = sources.getParameter<edm::InputTag>("jets" );
00026     mets_ = sources.getParameter<std::vector<edm::InputTag> >("mets" );
00027     pvs_ = sources.getParameter<edm::InputTag>("pvs" );
00028     // electronExtras are optional; they may be omitted or 
00029     // empty
00030     if( cfg.existsAs<edm::ParameterSet>("elecExtras") ){
00031       edm::ParameterSet elecExtras=cfg.getParameter<edm::ParameterSet>("elecExtras");
00032       // select is optional; in case it's not found no
00033       // selection will be applied
00034       if( elecExtras.existsAs<std::string>("select") ){
00035         elecSelect_= new StringCutObjectSelector<reco::GsfElectron>(elecExtras.getParameter<std::string>("select"));
00036       }
00037       // isolation is optional; in case it's not found no
00038       // isolation will be applied
00039       if( elecExtras.existsAs<std::string>("isolation") ){
00040         elecIso_= new StringCutObjectSelector<reco::GsfElectron>(elecExtras.getParameter<std::string>("isolation"));
00041       }
00042       // electronId is optional; in case it's not found the 
00043       // InputTag will remain empty
00044       if( elecExtras.existsAs<edm::ParameterSet>("electronId") ){
00045         edm::ParameterSet elecId=elecExtras.getParameter<edm::ParameterSet>("electronId");
00046         electronId_= elecId.getParameter<edm::InputTag>("src");
00047         eidPattern_= elecId.getParameter<int>("pattern");
00048       }
00049     }
00050     // pvExtras are opetional; they may be omitted or empty
00051     if(cfg.existsAs<edm::ParameterSet>("pvExtras")){
00052       edm::ParameterSet pvExtras=cfg.getParameter<edm::ParameterSet>("pvExtras");
00053       // select is optional; in case it's not found no
00054       // selection will be applied
00055       if( pvExtras.existsAs<std::string>("select") ){
00056         pvSelect_= new StringCutObjectSelector<reco::Vertex>(pvExtras.getParameter<std::string>("select"));
00057       }
00058     }
00059     // muonExtras are optional; they may be omitted or empty
00060     if( cfg.existsAs<edm::ParameterSet>("muonExtras") ){
00061       edm::ParameterSet muonExtras=cfg.getParameter<edm::ParameterSet>("muonExtras");
00062       // select is optional; in case it's not found no
00063       // selection will be applied
00064       if( muonExtras.existsAs<std::string>("select") ){
00065         muonSelect_= new StringCutObjectSelector<reco::Muon>(muonExtras.getParameter<std::string>("select"));
00066       }
00067       // isolation is optional; in case it's not found no
00068       // isolation will be applied
00069       if( muonExtras.existsAs<std::string>("isolation") ){
00070         muonIso_= new StringCutObjectSelector<reco::Muon>(muonExtras.getParameter<std::string>("isolation"));
00071       }
00072     }
00073     
00074     // jetExtras are optional; they may be omitted or 
00075     // empty
00076     if( cfg.existsAs<edm::ParameterSet>("jetExtras") ){
00077       edm::ParameterSet jetExtras=cfg.getParameter<edm::ParameterSet>("jetExtras");
00078       // jetCorrector is optional; in case it's not found 
00079       // the InputTag will remain empty
00080       if( jetExtras.existsAs<std::string>("jetCorrector") ){
00081         jetCorrector_= jetExtras.getParameter<std::string>("jetCorrector");
00082       }
00083       // read jetID information if it exists
00084       if(jetExtras.existsAs<edm::ParameterSet>("jetID")){
00085         edm::ParameterSet jetID=jetExtras.getParameter<edm::ParameterSet>("jetID");
00086         jetIDLabel_ =jetID.getParameter<edm::InputTag>("label");
00087         jetIDSelect_= new StringCutObjectSelector<reco::JetID>(jetID.getParameter<std::string>("select"));
00088       }
00089       // select is optional; in case it's not found no
00090       // selection will be applied (only implemented for 
00091       // CaloJets at the moment)
00092       if( jetExtras.existsAs<std::string>("select") ){
00093         jetSelect_= jetExtras.getParameter<std::string>("select");
00094       }
00095       // jetBDiscriminators are optional; in case they are
00096       // not found the InputTag will remain empty; they 
00097       // consist of pairs of edm::JetFlavorAssociation's & 
00098       // corresponding working points
00099       includeBTag_=jetExtras.existsAs<edm::ParameterSet>("jetBTaggers");
00100       if( includeBTag_ ){
00101         edm::ParameterSet btagEff=jetExtras.getParameter<edm::ParameterSet>("jetBTaggers").getParameter<edm::ParameterSet>("trackCountingEff");
00102         btagEff_= btagEff.getParameter<edm::InputTag>("label"); btagEffWP_= btagEff.getParameter<double>("workingPoint");
00103         edm::ParameterSet btagPur=jetExtras.getParameter<edm::ParameterSet>("jetBTaggers").getParameter<edm::ParameterSet>("trackCountingPur");
00104         btagPur_= btagPur.getParameter<edm::InputTag>("label"); btagPurWP_= btagPur.getParameter<double>("workingPoint");
00105         edm::ParameterSet btagVtx=jetExtras.getParameter<edm::ParameterSet>("jetBTaggers").getParameter<edm::ParameterSet>("secondaryVertex" );
00106         btagVtx_= btagVtx.getParameter<edm::InputTag>("label"); btagVtxWP_= btagVtx.getParameter<double>("workingPoint");
00107       }
00108     }
00109 
00110     // triggerExtras are optional; they may be omitted or empty
00111     if( cfg.existsAs<edm::ParameterSet>("triggerExtras") ){
00112       edm::ParameterSet triggerExtras=cfg.getParameter<edm::ParameterSet>("triggerExtras");
00113       triggerTable_=triggerExtras.getParameter<edm::InputTag>("src");
00114       triggerPaths_=triggerExtras.getParameter<std::vector<std::string> >("paths");
00115     }
00116 
00117     // massExtras is optional; in case it's not found no mass
00118     // window cuts are applied for the same flavor monitor
00119     // histograms
00120     if( cfg.existsAs<edm::ParameterSet>("massExtras") ){
00121       edm::ParameterSet massExtras=cfg.getParameter<edm::ParameterSet>("massExtras");
00122       lowerEdge_= massExtras.getParameter<double>("lowerEdge");
00123       upperEdge_= massExtras.getParameter<double>("upperEdge");
00124     }
00125 
00126     // setup the verbosity level for booking histograms;
00127     // per default the verbosity level will be set to 
00128     // STANDARD. This will also be the chosen level in
00129     // the case when the monitoring PSet is not found
00130     verbosity_=STANDARD;
00131     if( cfg.existsAs<edm::ParameterSet>("monitoring") ){
00132       edm::ParameterSet monitoring=cfg.getParameter<edm::ParameterSet>("monitoring");
00133       if(monitoring.getParameter<std::string>("verbosity") == "DEBUG"   )
00134         verbosity_= DEBUG;
00135       if(monitoring.getParameter<std::string>("verbosity") == "VERBOSE" )
00136         verbosity_= VERBOSE;
00137       if(monitoring.getParameter<std::string>("verbosity") == "STANDARD")
00138         verbosity_= STANDARD;
00139     }
00140     // and don't forget to do the histogram booking
00141     book(cfg.getParameter<std::string>("directory"));
00142   }
00143 
00144   void 
00145   MonitorEnsemble::book(std::string directory)
00146   {
00147     //set up the current directory path
00148     std::string current(directory); current+=label_;
00149     store_=edm::Service<DQMStore>().operator->();
00150     store_->setCurrentFolder(current);
00151 
00152     // determine number of bins for trigger monitoring
00153     unsigned int nPaths=triggerPaths_.size();
00154 
00155     // --- [STANDARD] --- //
00156     // number of selected primary vertices
00157     hists_["pvMult_"     ] = store_->book1D("PvMult"     , "N_{pvs}"          ,     100,     0.,    100.);  
00158     // pt of the leading muon
00159     hists_["muonPt_"     ] = store_->book1D("MuonPt"     , "pt(#mu)"          ,     50,     0.,    250.);   
00160     // muon multiplicity before std isolation
00161     hists_["muonMult_"   ] = store_->book1D("MuonMult"   , "N_{All}(#mu)"     ,     10,     0.,     10.);   
00162     // muon multiplicity after  std isolation
00163     hists_["muonMultIso_"] = store_->book1D("MuonMultIso", "N_{Iso}(#mu)"     ,     10,     0.,     10.);   
00164     // pt of the leading electron
00165     hists_["elecPt_"     ] = store_->book1D("ElecPt"     , "pt(e)"            ,     50,     0.,    250.);   
00166     // electron multiplicity before std isolation
00167     hists_["elecMult_"   ] = store_->book1D("ElecMult"   , "N_{All}(e)"       ,     10,     0.,     10.);   
00168     // electron multiplicity after  std isolation
00169     hists_["elecMultIso_"] = store_->book1D("ElecMultIso", "N_{Iso}(e)"       ,     10,     0.,     10.);   
00170     // multiplicity of jets with pt>20 (corrected to L2+L3)
00171     hists_["jetMult_"    ] = store_->book1D("JetMult"    , "N_{30}(jet)"      ,     10,     0.,     10.);   
00172     // trigger efficiency estimates for single lepton triggers
00173     hists_["triggerEff_" ] = store_->book1D("TriggerEff" , "Eff(trigger)"     , nPaths,     0.,  nPaths);
00174     // monitored trigger occupancy for single lepton triggers
00175     hists_["triggerMon_" ] = store_->book1D("TriggerMon" , "Mon(trigger)"     , nPaths,     0.,  nPaths);
00176     // MET (calo)
00177     hists_["metCalo_"    ] = store_->book1D("METCalo"    , "MET_{Calo}"       ,     50,     0.,    200.);   
00178     // W mass estimate
00179     hists_["massW_"      ] = store_->book1D("MassW"      , "M(W)"             ,     60,     0.,    300.);   
00180     // Top mass estimate
00181     hists_["massTop_"    ] = store_->book1D("MassTop"    , "M(Top)"           ,     50,     0.,    500.);   
00182 
00183     // set bin labels for trigger monitoring
00184     triggerBinLabels(std::string("trigger"), triggerPaths_);
00185 
00186     if( verbosity_==STANDARD) return;
00187 
00188     // --- [VERBOSE] --- //
00189     // eta of the leading muon
00190     hists_["muonEta_"    ] = store_->book1D("MuonEta"    , "#eta(#mu)"        ,     30,    -3.,      3.);   
00191     // std isolation variable of the leading muon
00192     hists_["muonRelIso_" ] = store_->book1D("MuonRelIso" , "Iso_{Rel}(#mu)"   ,     50,     0.,      1.);   
00193     // eta of the leading electron
00194     hists_["elecEta_"    ] = store_->book1D("ElecEta"    , "#eta(e)"          ,     30,    -3.,      3.);   
00195     // std isolation variable of the leading electron
00196     hists_["elecRelIso_" ] = store_->book1D("ElecRelIso" , "Iso_{Rel}(e)"     ,     50,     0.,      1.);   
00197     // multiplicity of btagged jets (for track counting high efficiency) with pt(L2L3)>30
00198     hists_["jetMultBEff_"] = store_->book1D("JetMultBEff", "N_{30}(b/eff)"    ,     10,     0.,     10.);   
00199     // btag discriminator for track counting high efficiency for jets with pt(L2L3)>30
00200     hists_["jetBDiscEff_"] = store_->book1D("JetBDiscEff", "Disc_{b/eff}(jet)",     100,     0.,     10.);   
00201     // pt of the 1. leading jet (corrected to L2+L3)
00202     hists_["jet1Pt_"     ] = store_->book1D("Jet1Pt"     , "pt_{L2L3}(jet1)"  ,     60,     0.,    300.);   
00203     // pt of the 2. leading jet (corrected to L2+L3)
00204     hists_["jet2Pt_"     ] = store_->book1D("Jet2Pt"     , "pt_{L2L3}(jet2)"  ,     60,     0.,    300.);   
00205     // pt of the 3. leading jet (corrected to L2+L3)
00206     hists_["jet3Pt_"     ] = store_->book1D("Jet3Pt"     , "pt_{L2L3}(jet3)"  ,     60,     0.,    300.);   
00207     // pt of the 4. leading jet (corrected to L2+L3)
00208     hists_["jet4Pt_"     ] = store_->book1D("Jet4Pt"     , "pt_{L2L3}(jet4)"  ,     60,     0.,    300.);   
00209     // MET (tc)
00210     hists_["metTC_"      ] = store_->book1D("METTC"      , "MET_{TC}"         ,     50,     0.,    200.);   
00211     // MET (pflow)
00212     hists_["metPflow_"   ] = store_->book1D("METPflow"   , "MET_{Pflow}"      ,     50,     0.,    200.);   
00213     // dz for muons (to suppress cosmis)
00214     hists_["muonDelZ_"    ] = store_->book1D("MuonDelZ"  , "d_{z}(#mu)"       ,     50,   -25.,     25.);
00215     // dxy for muons (to suppress cosmics)
00216     hists_["muonDelXY_"   ] = store_->book2D("MuonDelXY" , "d_{xy}(#mu)"      ,     50,   -0.1,     0.1,   50,   -0.1,   0.1);
00217 
00218     // set axes titles for dxy for muons
00219     hists_["muonDelXY_"   ]->setAxisTitle( "x [cm]", 1); hists_["muonDelXY_"   ]->setAxisTitle( "y [cm]", 2);
00220 
00221     if( verbosity_==VERBOSE) return;
00222 
00223     // --- [DEBUG] --- //
00224     // relative muon isolation in tracker for the leading muon
00225     hists_["muonTrkIso_" ] = store_->book1D("MuonTrkIso" , "Iso_{Trk}(#mu)"   ,     50,     0.,      1.);   
00226     // relative muon isolation in ecal+hcal for the leading muon
00227     hists_["muonCalIso_" ] = store_->book1D("MuonCalIso" , "Iso_{Ecal}(#mu)"  ,     50,     0.,      1.);   
00228     // relative electron isolation in tracker for the leading electron
00229     hists_["elecTrkIso_" ] = store_->book1D("ElecTrkIso" , "Iso_{Trk}(e)"     ,     50,     0.,      1.);   
00230     // relative electron isolation in ecal+hcal for the leading electron
00231     hists_["elecCalIso_" ] = store_->book1D("ElecCalIso" , "Iso_{Ecal}(e)"    ,     50,     0.,      1.);   
00232     // multiplicity of btagged jets (for track counting high purity) with pt(L2L3)>30
00233     hists_["jetMultBPur_"] = store_->book1D("JetMultBPur", "N_{30}(b/pur)"    ,     10,     0.,     10.);   
00234     // btag discriminator for track counting high purity
00235     hists_["jetBDiscPur_"] = store_->book1D("JetBDiscPur", "Disc_{b/pur}(Jet)",     100,     0.,     10.);   
00236     // multiplicity of btagged jets (for simple secondary vertex) with pt(L2L3)>30
00237     hists_["jetMultBVtx_"] = store_->book1D("JetMultBVtx", "N_{30}(b/vtx)"    ,     10,     0.,     10.);   
00238     // btag discriminator for simple secondary vertex
00239     hists_["jetBDiscVtx_"] = store_->book1D("JetBDiscVtx", "Disc_{b/vtx}(Jet)",     35,    -1.,      6.);   
00240     // pt of the 1. leading jet (uncorrected)
00241     hists_["jet1PtRaw_"  ] = store_->book1D("Jet1PtRaw"  , "pt_{Raw}(jet1)"   ,     60,     0.,    300.);   
00242     // pt of the 2. leading jet (uncorrected)
00243     hists_["jet2PtRaw_"  ] = store_->book1D("Jet2PtRaw"  , "pt_{Raw}(jet2)"   ,     60,     0.,    300.);   
00244     // pt of the 3. leading jet (uncorrected)
00245     hists_["jet3PtRaw_"  ] = store_->book1D("Jet3PtRaw"  , "pt_{Raw}(jet3)"   ,     60,     0.,    300.);   
00246     // pt of the 4. leading jet (uncorrected)
00247     hists_["jet4PtRaw_"  ] = store_->book1D("Jet4PtRaw"  , "pt_{Raw}(jet4)"   ,     60,     0.,    300.);   
00248     // selected events
00249     hists_["eventLogger_"] = store_->book2D("EventLogger", "Logged Events"    ,      9,     0.,      9.,   10,   0.,   10.);
00250 
00251     // set axes titles for selected events
00252     hists_["eventLogger_"]->getTH1()->SetOption("TEXT");
00253     hists_["eventLogger_"]->setBinLabel( 1 , "Run"             , 1);
00254     hists_["eventLogger_"]->setBinLabel( 2 , "Block"           , 1);
00255     hists_["eventLogger_"]->setBinLabel( 3 , "Event"           , 1);
00256     hists_["eventLogger_"]->setBinLabel( 4 , "pt_{L2L3}(jet1)" , 1);
00257     hists_["eventLogger_"]->setBinLabel( 5 , "pt_{L2L3}(jet2)" , 1);
00258     hists_["eventLogger_"]->setBinLabel( 6 , "pt_{L2L3}(jet3)" , 1);
00259     hists_["eventLogger_"]->setBinLabel( 7 , "pt_{L2L3}(jet4)" , 1);
00260     hists_["eventLogger_"]->setBinLabel( 8 , "M_{W}"           , 1);
00261     hists_["eventLogger_"]->setBinLabel( 9 , "M_{Top}"         , 1);
00262     hists_["eventLogger_"]->setAxisTitle("logged evts"         , 2);
00263     return;
00264   }
00265 
00266   void 
00267   MonitorEnsemble::fill(const edm::Event& event, const edm::EventSetup& setup)
00268   {
00269     // fetch trigger event if configured such 
00270     edm::Handle<edm::TriggerResults> triggerTable;
00271     if(!triggerTable_.label().empty()) {
00272       if( !event.getByLabel(triggerTable_, triggerTable) ) return;
00273     }
00274 
00275     /*
00276     ------------------------------------------------------------
00277 
00278     Primary Vertex Monitoring
00279 
00280     ------------------------------------------------------------
00281     */
00282     // fill monitoring plots for primary verices
00283     edm::Handle<edm::View<reco::Vertex> > pvs;
00284     if( !event.getByLabel(pvs_, pvs) ) return;
00285     unsigned int pvMult = 0;
00286     for(edm::View<reco::Vertex>::const_iterator pv=pvs->begin(); pv!=pvs->end(); ++pv){
00287       if(!pvSelect_ || (*pvSelect_)(*pv))
00288         pvMult++;
00289     }
00290     fill("pvMult_",    pvMult   );
00291 
00292     /* 
00293     ------------------------------------------------------------
00294 
00295     Electron Monitoring
00296 
00297     ------------------------------------------------------------
00298     */
00299 
00300     // fill monitoring plots for electrons
00301     edm::Handle<edm::View<reco::GsfElectron> > elecs;
00302     if( !event.getByLabel(elecs_, elecs) ) return;
00303 
00304     // check availability of electron id
00305     edm::Handle<edm::ValueMap<float> > electronId; 
00306     if(!electronId_.label().empty()) {
00307       if( !event.getByLabel(electronId_, electronId) ) return;
00308     }
00309 
00310     // loop electron collection
00311     unsigned int eMult=0, eMultIso=0;
00312     std::vector<const reco::GsfElectron*> isoElecs;
00313     for(edm::View<reco::GsfElectron>::const_iterator elec=elecs->begin(); elec!=elecs->end(); ++elec){
00314       unsigned int idx = elec-elecs->begin();
00315       // restrict to electrons with good electronId
00316       if( electronId_.label().empty() ? true : ((int)(*electronId)[elecs->refAt(idx)] & eidPattern_) ){
00317         if(!elecSelect_ || (*elecSelect_)(*elec)){
00318           double isolationTrk = elec->pt()/(elec->pt()+elec->dr03TkSumPt());
00319           double isolationCal = elec->pt()/(elec->pt()+elec->dr03EcalRecHitSumEt()+elec->dr03HcalTowerSumEt());
00320           double isolationRel = (elec->dr03TkSumPt()+elec->dr03EcalRecHitSumEt()+elec->dr03HcalTowerSumEt())/elec->pt();
00321           if( eMult==0 ){
00322             // restrict to the leading electron
00323             fill("elecPt_" , elec->pt() );
00324             fill("elecEta_", elec->eta());
00325             fill("elecRelIso_" , isolationRel );
00326             fill("elecTrkIso_" , isolationTrk );
00327             fill("elecCalIso_" , isolationCal );
00328           }
00329           // in addition to the multiplicity counter buffer the iso 
00330           // electron candidates for later overlap check with jets
00331           ++eMult; if(!elecIso_ || (*elecIso_)(*elec)){ isoElecs.push_back(&(*elec)); ++eMultIso;}
00332         }
00333       }
00334     }
00335     fill("elecMult_",    eMult   );
00336     fill("elecMultIso_", eMultIso);
00337     
00338     /* 
00339     ------------------------------------------------------------
00340 
00341     Muon Monitoring
00342 
00343     ------------------------------------------------------------
00344     */
00345 
00346     // fill monitoring plots for muons
00347     unsigned int mMult=0, mMultIso=0;
00348 
00349     edm::Handle<edm::View<reco::Muon> > muons;
00350     if( !event.getByLabel(muons_, muons) ) return;
00351 
00352     for(edm::View<reco::Muon>::const_iterator muon=muons->begin(); muon!=muons->end(); ++muon){
00353       // restrict to globalMuons
00354       if( muon->isGlobalMuon() ){ 
00355         fill("muonDelZ_" , muon->globalTrack()->vz());
00356         fill("muonDelXY_", muon->globalTrack()->vx(), muon->globalTrack()->vy());
00357         // apply preselection
00358         if(!muonSelect_ || (*muonSelect_)(*muon)){
00359           double isolationTrk = muon->pt()/(muon->pt()+muon->isolationR03().sumPt);
00360           double isolationCal = muon->pt()/(muon->pt()+muon->isolationR03().emEt+muon->isolationR03().hadEt);
00361           double isolationRel = (muon->isolationR03().sumPt+muon->isolationR03().emEt+muon->isolationR03().hadEt)/muon->pt();
00362           if( mMult==0 ){
00363             // restrict to leading muon
00364             fill("muonPt_"     , muon->pt() );
00365             fill("muonEta_"    , muon->eta());
00366             fill("muonRelIso_" , isolationRel );
00367             fill("muonTrkIso_" , isolationTrk );
00368             fill("muonCalIso_" , isolationCal );
00369           }
00370            ++mMult; if(!muonIso_ || (*muonIso_)(*muon)) ++mMultIso;
00371         }
00372       }
00373     }
00374     fill("muonMult_",    mMult   );
00375     fill("muonMultIso_", mMultIso);
00376 
00377     /* 
00378     ------------------------------------------------------------
00379 
00380     Jet Monitoring
00381 
00382     ------------------------------------------------------------
00383     */
00384 
00385     // check availability of the btaggers
00386     edm::Handle<reco::JetTagCollection> btagEff, btagPur, btagVtx;
00387     if( includeBTag_ ){ 
00388       if( !event.getByLabel(btagEff_, btagEff) ) return;
00389       if( !event.getByLabel(btagPur_, btagPur) ) return;
00390       if( !event.getByLabel(btagVtx_, btagVtx) ) return;
00391     }
00392     // load jet corrector if configured such
00393     const JetCorrector* corrector=0;
00394     if(!jetCorrector_.empty()){
00395       // check whether a jet correcto is in the event setup or not
00396       if(setup.find( edm::eventsetup::EventSetupRecordKey::makeKey<JetCorrectionsRecord>() )){
00397         corrector = JetCorrector::getJetCorrector(jetCorrector_, setup);
00398       }
00399       else{ 
00400         edm::LogVerbatim( "TopSingleLeptonDQM" ) 
00401           << "\n"
00402           << "------------------------------------------------------------------------------------- \n"
00403           << " No JetCorrectionsRecord available from EventSetup:                                   \n" 
00404           << "  - Jets will not be corrected.                                                       \n"
00405           << "  - If you want to change this add the following lines to your cfg file:              \n"
00406           << "                                                                                      \n"
00407           << "  ## load jet corrections                                                             \n"
00408           << "  process.load(\"JetMETCorrections.Configuration.JetCorrectionServicesAllAlgos_cff\") \n"
00409           << "  process.prefer(\"ak5CaloL2L3\")                                                     \n"
00410           << "                                                                                      \n"
00411           << "------------------------------------------------------------------------------------- \n";
00412       }
00413     }
00414 
00415     // loop jet collection
00416     std::vector<reco::Jet> correctedJets;
00417     unsigned int mult=0, multBEff=0, multBPur=0, multBVtx=0;
00418 
00419     edm::Handle<edm::View<reco::Jet> > jets; 
00420     if( !event.getByLabel(jets_, jets) ) return;
00421 
00422     edm::Handle<reco::JetIDValueMap> jetID; 
00423     if(jetIDSelect_){ 
00424       if( !event.getByLabel(jetIDLabel_, jetID) ) return;
00425     }
00426 
00427     for(edm::View<reco::Jet>::const_iterator jet=jets->begin(); jet!=jets->end(); ++jet){
00428       // check jetID for calo jets
00429       unsigned int idx = jet-jets->begin();
00430       if( jetIDSelect_ && dynamic_cast<const reco::CaloJet*>(jets->refAt(idx).get())){
00431         if(!(*jetIDSelect_)((*jetID)[jets->refAt(idx)])) continue;
00432       }
00433       // chekc additional jet selection for calo, pf and bare reco jets
00434       if(dynamic_cast<const reco::CaloJet*>(&*jet)){
00435         reco::CaloJet sel = dynamic_cast<const reco::CaloJet&>(*jet); sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
00436         StringCutObjectSelector<reco::CaloJet> jetSelect(jetSelect_); if(!jetSelect(sel)){ continue;}
00437       }
00438       else if(dynamic_cast<const reco::PFJet*>(&*jet)){
00439         reco::PFJet sel= dynamic_cast<const reco::PFJet&>(*jet); sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
00440         StringCutObjectSelector<reco::PFJet> jetSelect(jetSelect_); if(!jetSelect(sel)) continue;
00441       } 
00442       else{
00443         reco::Jet sel = *jet; sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
00444         StringCutObjectSelector<reco::Jet> jetSelect(jetSelect_); if(!jetSelect(sel)) continue;
00445       }
00446       // check for overlaps -- comment this to be synchronous with the selection
00447       //bool overlap=false;
00448       //for(std::vector<const reco::GsfElectron*>::const_iterator elec=isoElecs.begin(); elec!=isoElecs.end(); ++elec){
00449       //  if(reco::deltaR((*elec)->eta(), (*elec)->phi(), jet->eta(), jet->phi())<0.4){overlap=true; break;}
00450       //} if(overlap){continue;}
00451 
00452       // prepare jet to fill monitor histograms
00453       reco::Jet monitorJet = *jet; monitorJet.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
00454       correctedJets.push_back(monitorJet);
00455       ++mult; // determine jet multiplicity
00456       if( includeBTag_ ){
00457         // fill b-discriminators
00458         edm::RefToBase<reco::Jet> jetRef = jets->refAt(idx);    
00459         fill("jetBDiscEff_", (*btagEff)[jetRef]); if( (*btagEff)[jetRef]>btagEffWP_ ) ++multBEff; 
00460         fill("jetBDiscPur_", (*btagPur)[jetRef]); if( (*btagPur)[jetRef]>btagPurWP_ ) ++multBPur; 
00461         fill("jetBDiscVtx_", (*btagVtx)[jetRef]); if( (*btagVtx)[jetRef]>btagVtxWP_ ) ++multBVtx; 
00462       }
00463       // fill pt (raw or L2L3) for the leading four jets  
00464       if(idx==0) {fill("jet1Pt_" , monitorJet.pt()); fill("jet1PtRaw_", jet->pt() );}
00465       if(idx==1) {fill("jet2Pt_" , monitorJet.pt()); fill("jet2PtRaw_", jet->pt() );}
00466       if(idx==2) {fill("jet3Pt_" , monitorJet.pt()); fill("jet3PtRaw_", jet->pt() );}
00467       if(idx==3) {fill("jet4Pt_" , monitorJet.pt()); fill("jet4PtRaw_", jet->pt() );}
00468     }
00469     fill("jetMult_"    , mult    );
00470     fill("jetMultBEff_", multBEff);
00471     fill("jetMultBPur_", multBPur);
00472     fill("jetMultBVtx_", multBVtx);
00473     
00474     /* 
00475     ------------------------------------------------------------
00476 
00477     MET Monitoring
00478 
00479     ------------------------------------------------------------
00480     */
00481 
00482     // fill monitoring histograms for met
00483     for(std::vector<edm::InputTag>::const_iterator met_=mets_.begin(); met_!=mets_.end(); ++met_){
00484       edm::Handle<edm::View<reco::MET> > met;
00485       if( !event.getByLabel(*met_, met) ) continue;
00486       if(met->begin()!=met->end()){
00487         unsigned int idx=met_-mets_.begin();
00488         if(idx==0) fill("metCalo_" , met->begin()->et());
00489         if(idx==1) fill("metTC_"   , met->begin()->et());
00490         if(idx==2) fill("metPflow_", met->begin()->et());
00491       }
00492     }
00493 
00494     /* 
00495     ------------------------------------------------------------
00496 
00497     Event Monitoring
00498 
00499     ------------------------------------------------------------
00500     */
00501 
00502     // fill W boson and top mass estimates
00503     Calculate eventKinematics(MAXJETS, WMASS);
00504     double wMass   = eventKinematics.massWBoson  (correctedJets);
00505     double topMass = eventKinematics.massTopQuark(correctedJets);
00506     if(wMass>=0 && topMass>=0) {fill("massW_" , wMass  ); fill("massTop_" , topMass);}
00507     // fill plots for trigger monitoring
00508     if((lowerEdge_==-1. && upperEdge_==-1.) || (lowerEdge_<wMass && wMass<upperEdge_) ){
00509       if(!triggerTable_.label().empty()) fill(event, *triggerTable, "trigger", triggerPaths_);
00510       if(logged_<=hists_.find("eventLogger_")->second->getNbinsY()){
00511         // log runnumber, lumi block, event number & some
00512         // more pysics infomation for interesting events
00513         fill("eventLogger_", 0.5, logged_+0.5, event.eventAuxiliary().run()); 
00514         fill("eventLogger_", 1.5, logged_+0.5, event.eventAuxiliary().luminosityBlock()); 
00515         fill("eventLogger_", 2.5, logged_+0.5, event.eventAuxiliary().event()); 
00516         if(correctedJets.size()>0) fill("eventLogger_", 3.5, logged_+0.5, correctedJets[0].pt()); 
00517         if(correctedJets.size()>1) fill("eventLogger_", 4.5, logged_+0.5, correctedJets[1].pt()); 
00518         if(correctedJets.size()>2) fill("eventLogger_", 5.5, logged_+0.5, correctedJets[2].pt()); 
00519         if(correctedJets.size()>3) fill("eventLogger_", 6.5, logged_+0.5, correctedJets[3].pt()); 
00520         fill("eventLogger_", 7.5, logged_+0.5, wMass  ); 
00521         fill("eventLogger_", 8.5, logged_+0.5, topMass); 
00522         ++logged_;
00523       }
00524     }
00525   }
00526   
00527 }
00528 
00529 
00530 TopSingleLeptonDQM::TopSingleLeptonDQM(const edm::ParameterSet& cfg): triggerTable_(""), vertexSelect_(0), beamspot_(""), beamspotSelect_(0)
00531 {
00532   // configure preselection
00533   edm::ParameterSet presel=cfg.getParameter<edm::ParameterSet>("preselection");
00534   if( presel.existsAs<edm::ParameterSet>("trigger") ){
00535     edm::ParameterSet trigger=presel.getParameter<edm::ParameterSet>("trigger");
00536     triggerTable_=trigger.getParameter<edm::InputTag>("src");
00537     triggerPaths_=trigger.getParameter<std::vector<std::string> >("select");
00538   } 
00539   if( presel.existsAs<edm::ParameterSet>("vertex" ) ){
00540     edm::ParameterSet vertex=presel.getParameter<edm::ParameterSet>("vertex");
00541     vertex_= vertex.getParameter<edm::InputTag>("src");
00542     vertexSelect_= new StringCutObjectSelector<reco::Vertex>(vertex.getParameter<std::string>("select"));
00543   }
00544   if( presel.existsAs<edm::ParameterSet>("beamspot" ) ){
00545     edm::ParameterSet beamspot=presel.getParameter<edm::ParameterSet>("beamspot");
00546     beamspot_= beamspot.getParameter<edm::InputTag>("src");
00547     beamspotSelect_= new StringCutObjectSelector<reco::BeamSpot>(beamspot.getParameter<std::string>("select"));
00548   }
00549 
00550   // conifgure the selection
00551   std::vector<edm::ParameterSet> sel=cfg.getParameter<std::vector<edm::ParameterSet> >("selection");
00552   for(unsigned int i=0; i<sel.size(); ++i){
00553     selectionOrder_.push_back(sel.at(i).getParameter<std::string>("label"));
00554     selection_[selectionStep(selectionOrder_.back())] = std::make_pair(sel.at(i), new TopSingleLepton::MonitorEnsemble(selectionStep(selectionOrder_.back()).c_str(), cfg.getParameter<edm::ParameterSet>("setup")));
00555   }
00556 }
00557 
00558 void 
00559 TopSingleLeptonDQM::analyze(const edm::Event& event, const edm::EventSetup& setup)
00560 { 
00561   if(!triggerTable_.label().empty()){
00562     edm::Handle<edm::TriggerResults> triggerTable;
00563     if( !event.getByLabel(triggerTable_, triggerTable) ) return;
00564     if(!accept(event, *triggerTable, triggerPaths_)) return;
00565   }
00566   //cout<<"trig passed"<<endl;
00567   if(!beamspot_.label().empty()){
00568     edm::Handle<reco::BeamSpot> beamspot;
00569     if( !event.getByLabel(beamspot_, beamspot) ) return;
00570     if(!(*beamspotSelect_)(*beamspot)) return;
00571   }
00572   // apply selection steps
00573   unsigned int passed=0;
00574   for(std::vector<std::string>::const_iterator selIt=selectionOrder_.begin(); selIt!=selectionOrder_.end(); ++selIt){
00575     std::string key = selectionStep(*selIt), type = objectType(*selIt);
00576     if(selection_.find(key)!=selection_.end()){
00577       if(type=="empty"){
00578         selection_[key].second->fill(event, setup);
00579       }
00580       if(type=="muons"){
00581         SelectionStep<reco::Muon> step(selection_[key].first);
00582         if(step.select(event)){ ++passed;
00583           selection_[key].second->fill(event, setup);
00584         } else break;
00585       }
00586       if(type=="elecs"){
00587         SelectionStep<reco::GsfElectron> step(selection_[key].first);
00588         if(step.select(event)){ ++passed;
00589           selection_[key].second->fill(event, setup);
00590         } else break;
00591       }
00592       if(type=="pvs" ){
00593         SelectionStep<reco::Vertex> step(selection_[key].first);
00594         if(step.selectVertex(event)){ ++passed;
00595           selection_[key].second->fill(event, setup);
00596         } else break;
00597       }
00598       if(type=="jets" ){
00599         SelectionStep<reco::Jet> step(selection_[key].first);
00600         if(step.select(event, setup)){ ++passed;
00601           selection_[key].second->fill(event, setup);
00602         } else break;
00603       }
00604       if(type=="jets/pf" ){
00605         SelectionStep<reco::PFJet> step(selection_[key].first);
00606         if(step.select(event, setup)){ ++passed;
00607           selection_[key].second->fill(event, setup);
00608         } else break;
00609       }
00610       if(type=="jets/calo" ){
00611         SelectionStep<reco::CaloJet> step(selection_[key].first);
00612         if(step.select(event, setup)){ ++passed;
00613           selection_[key].second->fill(event, setup);
00614         } else break;
00615       }
00616       if(type=="met" ){
00617         SelectionStep<reco::MET> step(selection_[key].first);
00618         if(step.select(event)){ ++passed;
00619           selection_[key].second->fill(event, setup);
00620         } else break;
00621       }
00622     }
00623   }
00624 }
00625 
00626
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/DQM/Physics/src/TopSingleLeptonDQM.cc
