CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/DQM/Physics/src/TopSingleLeptonDQM.cc

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