CMS 3D CMS Logo

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