CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/DQM/Physics/src/TopDiLeptonOfflineDQM.cc

Go to the documentation of this file.
00001 //#include <algorithm>
00002 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
00003 #include "DataFormats/JetReco/interface/CaloJet.h"
00004 #include "DataFormats/BTauReco/interface/JetTag.h"
00005 #include "DQM/Physics/src/TopDiLeptonOfflineDQM.h"
00006 #include "DataFormats/JetReco/interface/PFJet.h"
00007 #include "DQM/Physics/interface/TopDQMHelpers.h"
00008 
00009 namespace TopDiLeptonOffline {
00010 
00011   MonitorEnsemble::MonitorEnsemble(const char* label, const edm::ParameterSet& cfg) : 
00012    label_(label), elecIso_(0), elecSelect_(0), muonIso_(0), muonSelect_(0), jetIDSelect_(0), 
00013    lowerEdge_(-1.), upperEdge_(-1.), elecMuLogged_(0), diMuonLogged_(0), diElecLogged_(0)
00014   {
00015     // sources have to be given; this PSet is not optional
00016     edm::ParameterSet sources=cfg.getParameter<edm::ParameterSet>("sources");
00017     muons_= sources.getParameter<edm::InputTag>("muons");
00018     elecs_= sources.getParameter<edm::InputTag>("elecs");
00019     jets_ = sources.getParameter<edm::InputTag>("jets" );
00020     mets_ = sources.getParameter<std::vector<edm::InputTag> >("mets" );
00021 
00022     // elecExtras are optional; they may be omitted or empty
00023     if( cfg.existsAs<edm::ParameterSet>("elecExtras") ){
00024       edm::ParameterSet elecExtras=cfg.getParameter<edm::ParameterSet>("elecExtras");
00025       // select is optional; in case it's not found no
00026       // selection will be applied
00027       if( elecExtras.existsAs<std::string>("select") ){
00028         elecSelect_= new StringCutObjectSelector<reco::GsfElectron>(elecExtras.getParameter<std::string>("select"));
00029       }
00030       // isolation is optional; in case it's not found no
00031       // isolation will be applied
00032       if( elecExtras.existsAs<std::string>("isolation") ){
00033         elecIso_= new StringCutObjectSelector<reco::GsfElectron>(elecExtras.getParameter<std::string>("isolation"));
00034       }
00035       // electronId is optional; in case it's not found the 
00036       // InputTag will remain empty
00037       if( elecExtras.existsAs<edm::InputTag>("electronId") ){
00038         electronId_= elecExtras.getParameter<edm::InputTag>("electronId");
00039       }
00040     }
00041     // muonExtras are optional; they may be omitted or empty
00042     if( cfg.existsAs<edm::ParameterSet>("muonExtras") ){
00043       edm::ParameterSet muonExtras=cfg.getParameter<edm::ParameterSet>("muonExtras");
00044       // select is optional; in case it's not found no
00045       // selection will be applied
00046       if( muonExtras.existsAs<std::string>("select") ){
00047         muonSelect_= new StringCutObjectSelector<reco::Muon>(muonExtras.getParameter<std::string>("select"));
00048       }
00049       // isolation is optional; in case it's not found no
00050       // isolation will be applied
00051       if( muonExtras.existsAs<std::string>("isolation") ){
00052         muonIso_= new StringCutObjectSelector<reco::Muon>(muonExtras.getParameter<std::string>("isolation"));
00053       }
00054     }
00055     // jetExtras are optional; they may be omitted or empty
00056     if( cfg.existsAs<edm::ParameterSet>("jetExtras") ){
00057       edm::ParameterSet jetExtras=cfg.getParameter<edm::ParameterSet>("jetExtras");
00058       // jetCorrector is optional; in case it's not found 
00059       // the InputTag will remain empty
00060       if( jetExtras.existsAs<std::string>("jetCorrector") ){
00061         jetCorrector_= jetExtras.getParameter<std::string>("jetCorrector");
00062       }
00063       // read jetID information if it exists
00064       if(jetExtras.existsAs<edm::ParameterSet>("jetID")){
00065         edm::ParameterSet jetID=jetExtras.getParameter<edm::ParameterSet>("jetID");
00066         jetIDLabel_ =jetID.getParameter<edm::InputTag>("label");
00067         jetIDSelect_= new StringCutObjectSelector<reco::JetID>(jetID.getParameter<std::string>("select"));
00068       }
00069       // select is optional; in case it's not found no
00070       // selection will be applied (only implemented for 
00071       // CaloJets at the moment)
00072       if( jetExtras.existsAs<std::string>("select") ){
00073         jetSelect_= jetExtras.getParameter<std::string>("select");
00074       }
00075     }
00076     // triggerExtras are optional; they may be omitted or empty
00077     if( cfg.existsAs<edm::ParameterSet>("triggerExtras") ){
00078       edm::ParameterSet triggerExtras=cfg.getParameter<edm::ParameterSet>("triggerExtras");
00079       triggerTable_=triggerExtras.getParameter<edm::InputTag>("src");
00080       elecMuPaths_ =triggerExtras.getParameter<std::vector<std::string> >("pathsELECMU");
00081       diMuonPaths_ =triggerExtras.getParameter<std::vector<std::string> >("pathsDIMUON");
00082     }
00083     // massExtras is optional; in case it's not found no mass
00084     // window cuts are applied for the same flavor monitor
00085     // histograms
00086     if( cfg.existsAs<edm::ParameterSet>("massExtras") ){
00087       edm::ParameterSet massExtras=cfg.getParameter<edm::ParameterSet>("massExtras");
00088       lowerEdge_= massExtras.getParameter<double>("lowerEdge");
00089       upperEdge_= massExtras.getParameter<double>("upperEdge");
00090     }
00091 
00092     // setup the verbosity level for booking histograms;
00093     // per default the verbosity level will be set to 
00094     // STANDARD. This will also be the chosen level in
00095     // the case when the monitoring PSet is not found
00096     verbosity_=STANDARD;
00097     if( cfg.existsAs<edm::ParameterSet>("monitoring") ){
00098       edm::ParameterSet monitoring=cfg.getParameter<edm::ParameterSet>("monitoring");
00099       if(monitoring.getParameter<std::string>("verbosity") == "DEBUG"   )
00100         verbosity_= DEBUG;
00101       if(monitoring.getParameter<std::string>("verbosity") == "VERBOSE" )
00102         verbosity_= VERBOSE;
00103       if(monitoring.getParameter<std::string>("verbosity") == "STANDARD")
00104         verbosity_= STANDARD;
00105     }
00106     // and don't forget to do the histogram booking
00107     book(cfg.getParameter<std::string>("directory"));
00108   }
00109 
00110   void 
00111   MonitorEnsemble::book(std::string directory)
00112   {
00113     //set up the current directory path
00114     std::string current(directory); current+=label_;
00115     store_=edm::Service<DQMStore>().operator->();
00116     store_->setCurrentFolder(current);
00117 
00118     // determine number of bins for trigger monitoring
00119     unsigned int nElecMu=elecMuPaths_.size();
00120     unsigned int nDiMuon=diMuonPaths_.size();
00121 
00122     // --- [STANDARD] --- //
00123     // invariant mass of opposite charge lepton pair (only filled for same flavor)
00124     hists_["invMass_"     ] = store_->book1D("InvMass"     , "M(lep1, lep2)"           ,       80,   0.,     320.);
00125     // invariant mass of opposite charge lepton pair (only filled for same flavor)
00126     hists_["invMassLog_"  ] = store_->book1D("InvMassLog"  , "log_{10}(M(lep1, lep2))" ,       80,   .1,      2.5);
00127     // invariant mass of same charge lepton pair (log10 for low mass region, only filled for same flavor)
00128     hists_["invMassWC_"   ] = store_->book1D("InvMassWC"   , "M_{WC}(L1, L2)"          ,       80,   0.,     320.);
00129     // invariant mass of same charge lepton pair (log10 for low mass region, only filled for same flavor)
00130     hists_["invMassWCLog_"] = store_->book1D("InvMassLogWC", "log_{10}(M_{WC})"        ,       80,   .1,      2.5);
00131     // decay channel [1]: muon/muon, [2]:elec/elec, [3]:elec/muon 
00132     hists_["decayChannel_"] = store_->book1D("DecayChannel", "Decay Channel"           ,        3,    0,        3);
00133     // trigger efficiency estimates for the electron muon channel
00134     hists_["elecMuEff_"   ] = store_->book1D("ElecMuEff"   , "Eff(e/#mu paths)"        ,  nElecMu,   0.,  nElecMu);
00135     // monitored trigger occupancy for the electron muon channel
00136     hists_["elecMuMon_"   ] = store_->book1D("ElecMuMon"   , "Mon(e/#mu paths)"        ,  nElecMu,   0.,  nElecMu);
00137     // trigger efficiency estimates for the di muon channel
00138     hists_["diMuonEff_"   ] = store_->book1D("DiMuonEff"   , "Eff(#mu/#mu paths)"      ,  nDiMuon,   0.,  nDiMuon);
00139     // monitored trigger occupancy for the di muon channel
00140     hists_["diMuonMon_"   ] = store_->book1D("DiMuonMon"   , "Mon(#mu/#mu paths)"      ,  nDiMuon,   0.,  nDiMuon);
00141     // pt of the leading lepton
00142     hists_["lep1Pt_"      ] = store_->book1D("LeptPt"      , "pt(lep1)"                ,       50,   0.,     200.);
00143     // pt of the 2. leading lepton
00144     hists_["lep2Pt_"      ] = store_->book1D("Lep2Pt"      , "pt(lep2)"                ,       50,   0.,     200.);
00145     // multiplicity of jets with pt>30 (corrected to L2+L3)
00146     hists_["jetMult_"     ] = store_->book1D("JetMult"     , "N_{30}(jet)"             ,       10,   0.,      10.); 
00147     // MET (calo)
00148     hists_["metCalo_"     ] = store_->book1D("METCalo"     , "MET_{Calo}"              ,       50,   0.,     200.);
00149 
00150     // set bin labels for trigger monitoring
00151     triggerBinLabels(std::string("elecMu"), elecMuPaths_);
00152     triggerBinLabels(std::string("diMuon"), diMuonPaths_);
00153     // set bin labels for decayChannel_
00154     hists_["decayChannel_"]->setBinLabel( 1, "#mu e"  , 1);
00155     hists_["decayChannel_"]->setBinLabel( 2, "#mu #mu", 1);
00156     hists_["decayChannel_"]->setBinLabel( 3, "e e"    , 1);
00157 
00158     if( verbosity_==STANDARD) return;
00159 
00160     // --- [VERBOSE] --- //
00161     // mean eta of the candidate leptons
00162     hists_["sumEtaL1L2_"  ] = store_->book1D("SumEtaL1L2"  , "<#eta>(lep1, lep2)"      ,       30,  -5.,       5.); 
00163     // deltaEta between the 2 candidate leptons
00164     hists_["dEtaL1L2_"    ] = store_->book1D("DEtaL1L2"    , "#Delta#eta(lep1,lep2)"   ,       30,   0.,       3.);
00165     // deltaPhi between the 2 candidate leptons
00166     hists_["dPhiL1L2_"    ] = store_->book1D("DPhiL1L2"    , "#Delta#phi(lep1,lep2)"   ,       32,   0.,      3.2);
00167     // pt of the candidate electron (depending on the decay channel)
00168     hists_["elecPt_"      ] = store_->book1D("ElecPt"      , "pt(e)"                   ,       50,   0.,     200.);
00169     // relative isolation of the candidate electron (depending on the decay channel)
00170     hists_["elecRelIso_"  ] = store_->book1D("ElecRelIso"  , "Iso_{Rel}(e)"            ,       50,   0.,       1.);
00171     // pt of the canddiate muon (depending on the decay channel)
00172     hists_["muonPt_"      ] = store_->book1D("MuonPt"      , "pt(#mu)"                 ,       50,   0.,     200.);
00173     // relative isolation of the candidate muon (depending on the decay channel)
00174     hists_["muonRelIso_"  ] = store_->book1D("MuonRelIso"  , "Iso_{Rel}(#mu)"          ,       50,   0.,       1.);
00175     // pt of the 1. leading jet (corrected to L2+L3)
00176     hists_["jet1Pt_"      ] = store_->book1D("Jet1Pt"      , "pt_{L2L3}(jet1)"         ,       60,   0.,     300.);   
00177     // pt of the 2. leading jet (corrected to L2+L3)
00178     hists_["jet2Pt_"      ] = store_->book1D("Jet2Pt"      , "pt_{L2L3}(jet2)"         ,       60,   0.,     300.);
00179     // MET (PF)
00180     hists_["metPflow_"    ] = store_->book1D("METPflow"    , "MET_{Pflow}"             ,       50,   0.,     200.);
00181     // MET (TC)
00182     hists_["metTC_"       ] = store_->book1D("METTC"       , "MET_{TC}"                ,       50,   0.,     200.);
00183     // dz for muons (to suppress cosmis)
00184     hists_["muonDelZ_"    ] = store_->book1D("MuonDelZ"    , "d_{z}(#mu)"              ,       50, -25.,      25.);
00185     // dxy for muons (to suppress cosmics)
00186     hists_["muonDelXY_"   ] = store_->book2D("MuonDelXY"   , "d_{xy}(#mu)"                , 50,  -0.1,  0.1, 50, -0.1,  0.1 );
00187     // lepton multiplicity after std isolation
00188     hists_["lepMultIso_"  ] = store_->book2D("LepMultIso"  , "N_{Iso}(e) vs N_{Iso}(#mu)" ,  5,    0.,   5.,  5,   0.,    5.);
00189 
00190     // set axes titles for dxy for muons
00191     hists_["muonDelXY_"   ]->setAxisTitle( "x [cm]", 1); hists_["muonDelXY_"   ]->setAxisTitle( "y [cm]", 2);
00192     // set axes titles for lepton multiplicity after std isolation
00193     hists_["lepMultIso_"  ]->setAxisTitle( "N_{Iso}(#mu)", 1); hists_["lepMultIso_"   ]->setAxisTitle( "N_{Iso}(elec)", 2);
00194 
00195     if( verbosity_==VERBOSE) return;
00196 
00197     // --- [DEBUG] --- //
00198     // electron multiplicity after std isolation
00199     hists_["elecMultIso_" ] = store_->book1D("ElecMultIso" , "N_{Iso}(e)"              ,       10,   0.,      10.);
00200     // muon multiplicity after std isolation
00201     hists_["muonMultIso_" ] = store_->book1D("MuonMultIso" , "N_{Iso}(#mu)"            ,       10,   0.,      10.);
00202     // calo isolation of the candidate muon (depending on the decay channel)
00203     hists_["muonCalIso_"  ] = store_->book1D("MuonCalIso"  , "Iso_{Cal}(#mu)"          ,       50,   0.,       1.);
00204     // track isolation of the candidate muon (depending on the decay channel)
00205     hists_["muonTrkIso_"  ] = store_->book1D("MuonTrkIso"  , "Iso_{Trk}(#mu)"          ,       50,   0.,       1.);
00206     // calo isolation of the candidate electron (depending on the decay channel)
00207     hists_["elecCalIso_"  ] = store_->book1D("ElecCalIso"  , "Iso_{Cal}(e)"            ,       50,   0.,       1.);
00208     // track isolation of the candidate electron (depending on the decay channel)
00209     hists_["elecTrkIso_"  ] = store_->book1D("ElecTrkIso"  , "Iso_{Trk}(e)"            ,       50,   0.,       1.);
00210     // eta of the leading jet
00211     hists_["jet1Eta_"     ] = store_->book1D("Jet1Eta"     , "#eta(jet1)"              ,       30,  -5.,       5.); 
00212     // eta of the 2. leading jet
00213     hists_["jet2Eta_"     ] = store_->book1D("Jet2Eta"     , "#eta(jet2)"              ,       30,  -5.,       5.);
00214     // pt of the 1. leading jet (not corrected)
00215     hists_["jet1PtRaw_"   ] = store_->book1D("Jet1PtRaw"   , "pt_{Raw}(jet1)"          ,       60,   0.,     300.);   
00216     // pt of the 2. leading jet (not corrected)     
00217     hists_["jet2PtRaw_"   ] = store_->book1D("Jet2PtRaw"   , "pt_{Raw}(jet2)"          ,       60,   0.,     300.);
00218     // deltaEta between the 2 leading jets
00219     hists_["dEtaJet1Jet2_"] = store_->book1D("DEtaJet1Jet2", "#Delta#eta(jet1,jet2)"   ,       30,   0.,       3.);
00220     // deltaEta between the lepton and the leading jet
00221     hists_["dEtaJet1Lep1_"] = store_->book1D("DEtaJet1Lep1", "#Delta#eta(jet1,lep1)"   ,       30,   0.,       3.);
00222     // deltaEta between the lepton and MET
00223     hists_["dEtaLep1MET_" ] = store_->book1D("DEtaLep1MET" , "#Delta#eta(lep1,MET)"    ,       30,   0.,       3.);
00224     // deltaEta between leading jet and MET
00225     hists_["dEtaJet1MET_" ] = store_->book1D("DEtaJet1MET" , "#Delta#eta(jet1,MET)"    ,       30,   0.,       3.);
00226     // deltaPhi of 2 leading jets
00227     hists_["dPhiJet1Jet2_"] = store_->book1D("DPhiJet1Jet2", "#Delta#phi(jet1,jet2)"   ,       32,   0.,      3.2);
00228     // deltaPhi of 1. lepton and 1. jet
00229     hists_["dPhiJet1Lep1_"] = store_->book1D("DPhiJet1Lep1", "#Delta#phi(jet1,lep1)"   ,       32,   0.,      3.2);
00230     // deltaPhi of 1. lepton and MET
00231     hists_["dPhiLep1MET_" ] = store_->book1D("DPhiLep1MET" , "#Delta#phi(lep1,MET)"    ,       32,   0.,      3.2);
00232     // deltaPhi of 1. jet and MET
00233     hists_["dPhiJet1MET_" ] = store_->book1D("DPhiJet1MET" , "#Delta#phi(jet1,MET)"    ,       32,   0.,      3.2);
00234     // selected dimuon events
00235     hists_["diMuonLogger_"] = store_->book2D("DiMuonLogger", "Logged DiMuon Events"    ,        8,   0.,       8.,   10,   0.,   10.);
00236     // selected dielec events
00237     hists_["diElecLogger_"] = store_->book2D("DiElecLogger", "Logged DiElec Events"    ,        8,   0.,       8.,   10,   0.,   10.);
00238     // selected elemu events
00239     hists_["elecMuLogger_"] = store_->book2D("ElecMuLogger", "Logged ElecMu Events"    ,        8,   0.,       8.,   10,   0.,   10.);
00240 
00241     // set bin labels for trigger monitoring
00242     loggerBinLabels(std::string("diMuonLogger_")); 
00243     loggerBinLabels(std::string("diElecLogger_")); 
00244     loggerBinLabels(std::string("elecMuLogger_"));
00245     return;
00246   }
00247 
00248   void 
00249   MonitorEnsemble::fill(const edm::Event& event, const edm::EventSetup& setup)
00250   {
00251     // fetch trigger event if configured such 
00252     edm::Handle<edm::TriggerResults> triggerTable;
00253     if(!triggerTable_.label().empty()) {
00254       if( !event.getByLabel(triggerTable_, triggerTable) ) return;
00255     }
00256 
00257     /* 
00258     ------------------------------------------------------------
00259 
00260     Muon Selection
00261 
00262     ------------------------------------------------------------
00263     */
00264 
00265     // buffer isolated muons
00266     std::vector<const reco::Muon*> isoMuons;
00267 
00268     edm::Handle<edm::View<reco::Muon> > muons;
00269     if( !event.getByLabel(muons_, muons) ) return;
00270 
00271     for(edm::View<reco::Muon>::const_iterator muon=muons->begin(); muon!=muons->end(); ++muon){
00272       // restrict to globalMuons
00273       if( muon->isGlobalMuon() ){ 
00274         fill("muonDelZ_" , muon->globalTrack()->vz());
00275         fill("muonDelXY_", muon->globalTrack()->vx(), muon->globalTrack()->vy());
00276         // apply preselection
00277         if(!muonSelect_ || (*muonSelect_)(*muon)){
00278           double isolationTrk = muon->pt()/(muon->pt()+muon->isolationR03().sumPt);
00279           double isolationCal = muon->pt()/(muon->pt()+muon->isolationR03().emEt+muon->isolationR03().hadEt);
00280           double isolationRel = (muon->isolationR03().sumPt+muon->isolationR03().emEt+muon->isolationR03().hadEt)/muon->pt();
00281           fill("muonTrkIso_" , isolationTrk); fill("muonCalIso_" , isolationCal); fill("muonRelIso_" , isolationRel);
00282           if(!muonIso_ || (*muonIso_)(*muon)) isoMuons.push_back(&(*muon));
00283         }
00284       }
00285     }
00286     fill("muonMultIso_", isoMuons.size());
00287 
00288     /* 
00289     ------------------------------------------------------------
00290 
00291     Electron Selection
00292 
00293     ------------------------------------------------------------
00294     */
00295 
00296     // buffer isolated electronss
00297     std::vector<const reco::GsfElectron*> isoElecs;
00298     edm::Handle<edm::ValueMap<float> > electronId; 
00299     if(!electronId_.label().empty()) {
00300       if( !event.getByLabel(electronId_, electronId) ) return;
00301     }
00302 
00303     edm::Handle<edm::View<reco::GsfElectron> > elecs;
00304     if( !event.getByLabel(elecs_, elecs) ) return;
00305 
00306     for(edm::View<reco::GsfElectron>::const_iterator elec=elecs->begin(); elec!=elecs->end(); ++elec){
00307       // restrict to electrons with good electronId
00308       int idx = elec-elecs->begin();
00309       if( electronId_.label().empty() ? true : (*electronId)[elecs->refAt(idx)]>0.99 ){
00310         // apply preselection
00311         if(!elecSelect_ || (*elecSelect_)(*elec)){
00312           double isolationTrk = elec->pt()/(elec->pt()+elec->dr03TkSumPt());
00313           double isolationCal = elec->pt()/(elec->pt()+elec->dr03EcalRecHitSumEt()+elec->dr03HcalTowerSumEt());
00314           double isolationRel = (elec->dr03TkSumPt()+elec->dr03EcalRecHitSumEt()+elec->dr03HcalTowerSumEt())/elec->pt();
00315           fill("elecTrkIso_" , isolationTrk); fill("elecCalIso_" , isolationCal); fill("elecRelIso_" , isolationRel);
00316           if(!elecIso_ || (*elecIso_)(*elec)) isoElecs.push_back(&(*elec));
00317         }
00318       }
00319     }
00320     fill("elecMultIso_", isoElecs.size());
00321 
00322     /* 
00323     ------------------------------------------------------------
00324 
00325     Jet Selection
00326 
00327     ------------------------------------------------------------
00328     */
00329 
00330     const JetCorrector* corrector=0;
00331     if(!jetCorrector_.empty()){
00332       // check whether a jet correcto is in the event setup or not
00333       if(setup.find( edm::eventsetup::EventSetupRecordKey::makeKey<JetCorrectionsRecord>() )){
00334         corrector = JetCorrector::getJetCorrector(jetCorrector_, setup);
00335       }
00336       else{
00337         edm::LogVerbatim( "TopDiLeptonOfflineDQM" ) 
00338           << "\n"
00339           << "------------------------------------------------------------------------------------- \n"
00340           << " No JetCorrectionsRecord available from EventSetup:                                   \n" 
00341           << "  - Jets will not be corrected.                                                       \n"
00342           << "  - If you want to change this add the following lines to your cfg file:              \n"
00343           << "                                                                                      \n"
00344           << "  ## load jet corrections                                                             \n"
00345           << "  process.load(\"JetMETCorrections.Configuration.JetCorrectionServicesAllAlgos_cff\") \n"
00346           << "  process.prefer(\"ak5CaloL2L3\")                                                     \n"
00347           << "                                                                                      \n"
00348           << "------------------------------------------------------------------------------------- \n";
00349       }
00350     }
00351 
00352     unsigned int mult=0;
00353     // buffer leadingJets
00354     std::vector<reco::Jet> leadingJets;
00355     edm::Handle<edm::View<reco::Jet> > jets; 
00356     if( !event.getByLabel(jets_, jets) ) return;
00357 
00358     edm::Handle<reco::JetIDValueMap> jetID;
00359     if(jetIDSelect_){ 
00360       if( !event.getByLabel(jetIDLabel_, jetID) ) return;
00361     }
00362 
00363     for(edm::View<reco::Jet>::const_iterator jet=jets->begin(); jet!=jets->end(); ++jet){
00364       unsigned int idx=jet-jets->begin();
00365       if( jetIDSelect_ && dynamic_cast<const reco::CaloJet*>(jets->refAt(idx).get())){
00366         if(!(*jetIDSelect_)((*jetID)[jets->refAt(idx)])) continue;
00367       }
00368       // chekc additional jet selection for calo, pf and bare reco jets
00369       if(dynamic_cast<const reco::CaloJet*>(&*jet)){
00370         reco::CaloJet sel = dynamic_cast<const reco::CaloJet&>(*jet); sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
00371         StringCutObjectSelector<reco::CaloJet> jetSelect(jetSelect_); if(!jetSelect(sel)){ continue;}
00372       }
00373       else if(dynamic_cast<const reco::PFJet*>(&*jet)){
00374         reco::PFJet sel= dynamic_cast<const reco::PFJet&>(*jet); sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
00375         StringCutObjectSelector<reco::PFJet> jetSelect(jetSelect_); if(!jetSelect(sel)) continue;
00376       } 
00377       else{
00378         reco::Jet sel = *jet; sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
00379         StringCutObjectSelector<reco::Jet> jetSelect(jetSelect_); if(!jetSelect(sel)) continue;
00380       }
00381       // check for overlaps
00382       bool overlap=false;
00383       for(std::vector<const reco::GsfElectron*>::const_iterator elec=isoElecs.begin(); elec!=isoElecs.end(); ++elec){
00384         if(reco::deltaR((*elec)->eta(), (*elec)->phi(), jet->eta(), jet->phi())<0.4){overlap=true; break;}
00385       } if(overlap){continue;}
00386       // prepare jet to fill monitor histograms
00387       reco::Jet monitorJet=*jet; monitorJet.scaleEnergy(corrector ?  corrector->correction(*jet) : 1.);
00388       ++mult; // determine jet multiplicity
00389       if(idx==0) {
00390         leadingJets.push_back(monitorJet);
00391         fill("jet1Pt_"    , monitorJet.pt());
00392         fill("jet1PtRaw_" , jet->pt() );
00393         fill("jet1Eta_"   , jet->eta());
00394       }
00395       if(idx==1) {
00396         leadingJets.push_back(monitorJet);
00397         fill("jet2Pt_"    , monitorJet.pt());
00398         fill("jet2PtRaw_" , jet->pt() );
00399         fill("jet2Eta_"   , jet->eta());
00400       }
00401     }
00402     if(leadingJets.size()>1){
00403       fill("dEtaJet1Jet2_" , leadingJets[0].eta()-leadingJets[1].eta());
00404       fill("dPhiJet1Jet2_" , reco::deltaPhi(leadingJets[0].phi(), leadingJets[1].phi()));
00405       if( !isoMuons.empty() ){
00406         if( isoElecs.empty() || isoMuons[0]->pt()>isoElecs[0]->pt() ){
00407           fill("dEtaJet1Lep1_" , isoMuons[0]->eta()-leadingJets[0].eta());
00408           fill("dPhiJet1Lep1_" , reco::deltaPhi(isoMuons[0]->phi() , leadingJets[0].phi()));
00409         } 
00410       }
00411       if( !isoElecs.empty() ){
00412         if( isoMuons.empty() || isoElecs[0]->pt()>isoMuons[0]->pt() ){
00413           fill("dEtaJet1Lep1_" , isoElecs[0]->eta()-leadingJets[0].eta());
00414           fill("dPhiJet1Lep1_" , reco::deltaPhi(isoElecs[0]->phi() , leadingJets[0].phi()));
00415         }
00416       }
00417     }
00418     fill("jetMult_", mult);
00419     
00420     /* 
00421     ------------------------------------------------------------
00422 
00423     MET Selection
00424 
00425     ------------------------------------------------------------
00426     */
00427 
00428     // buffer for event logging 
00429     reco::MET caloMET;
00430     for(std::vector<edm::InputTag>::const_iterator met_=mets_.begin(); met_!=mets_.end(); ++met_){
00431 
00432       edm::Handle<edm::View<reco::MET> > met;
00433       if( !event.getByLabel(*met_, met) ) continue;
00434 
00435       if(met->begin()!=met->end()){
00436         unsigned int idx=met_-mets_.begin();
00437         if(idx==0){
00438           caloMET=*met->begin(); 
00439           fill("metCalo_", met->begin()->et());
00440           if(!leadingJets.empty()){
00441             fill("dEtaJet1MET_" , leadingJets[0].eta()-met->begin()->eta());
00442             fill("dPhiJet1MET_" , reco::deltaPhi(leadingJets[0].phi(), met->begin()->phi()));
00443           }
00444           if( !isoMuons.empty() ){
00445             if( isoElecs.empty() || isoMuons[0]->pt()>isoElecs[0]->pt() ){
00446               fill("dEtaLep1MET_" , isoMuons[0]->eta()-met->begin()->eta());
00447               fill("dPhiLep1MET_" , reco::deltaPhi(isoMuons[0]->phi(), met->begin()->phi()));
00448             } 
00449           }
00450           if( !isoElecs.empty() ){
00451             if( isoMuons.empty() || isoElecs[0]->pt()>isoMuons[0]->pt() ){
00452               fill("dEtaLep1MET_" , isoElecs[0]->eta()-met->begin()->eta());
00453               fill("dPhiLep1MET_" , reco::deltaPhi(isoElecs[0]->phi(), met->begin()->phi()));
00454             }
00455           }
00456         }
00457         if(idx==1){ fill("metTC_"   , met->begin()->et());}
00458         if(idx==2){ fill("metPflow_", met->begin()->et());}
00459       }
00460     }
00461 
00462 
00463     /* 
00464     ------------------------------------------------------------
00465 
00466     Event Monitoring
00467 
00468     ------------------------------------------------------------
00469     */
00470 
00471     // check number of isolated leptons
00472     fill("lepMultIso_", isoMuons.size(), isoElecs.size());
00473     // ELECMU channel
00474     if( decayChannel(isoMuons, isoElecs) == ELECMU ){
00475       fill("decayChannel_", 0.5);
00476       double mass = (isoElecs[0]->p4()+isoMuons[0]->p4()).mass();
00477       if( (lowerEdge_==-1. && upperEdge_==-1.) || (lowerEdge_<mass && mass<upperEdge_) ){
00478         fill("dEtaL1L2_"  , isoElecs[0]->eta()-isoMuons[0]->eta()); 
00479         fill("sumEtaL1L2_", (isoElecs[0]->eta()+isoMuons[0]->eta())/2); 
00480         fill("dPhiL1L2_"  , reco::deltaPhi(isoElecs[0]->phi(), isoMuons[0]->eta())); 
00481         fill("elecPt_", isoElecs[0]->pt()); fill("muonPt_", isoMuons[0]->pt()); 
00482         fill("lep1Pt_", isoElecs[0]->pt()>isoMuons[0]->pt() ? isoElecs[0]->pt() : isoMuons[0]->pt());
00483         fill("lep2Pt_", isoElecs[0]->pt()>isoMuons[0]->pt() ? isoMuons[0]->pt() : isoElecs[0]->pt());
00484         // fill plots for trigger monitoring
00485         if(!triggerTable_.label().empty()) fill(event, *triggerTable, "elecMu", elecMuPaths_);
00486         if(elecMuLogged_<=hists_.find("elecMuLogger_")->second->getNbinsY()){
00487           // log runnumber, lumi block, event number & some
00488           // more pysics infomation for interesting events
00489           fill("elecMuLogger_", 0.5, elecMuLogged_+0.5, event.eventAuxiliary().run()); 
00490           fill("elecMuLogger_", 1.5, elecMuLogged_+0.5, event.eventAuxiliary().luminosityBlock()); 
00491           fill("elecMuLogger_", 2.5, elecMuLogged_+0.5, event.eventAuxiliary().event()); 
00492           fill("elecMuLogger_", 3.5, elecMuLogged_+0.5, isoMuons[0]->pt()); 
00493           fill("elecMuLogger_", 4.5, elecMuLogged_+0.5, isoElecs[0]->pt()); 
00494           if(leadingJets.size()>0) fill("elecMuLogger_", 5.5, elecMuLogged_+0.5, leadingJets[0].pt()); 
00495           if(leadingJets.size()>1) fill("elecMuLogger_", 6.5, elecMuLogged_+0.5, leadingJets[1].pt()); 
00496           fill("elecMuLogger_", 7.5, elecMuLogged_+0.5, caloMET.et()); 
00497           ++elecMuLogged_; 
00498         }
00499       }
00500     }
00501 
00502     // DIMUON channel
00503     if( decayChannel(isoMuons, isoElecs) == DIMUON ){
00504       fill("decayChannel_", 1.5);
00505       int charge = isoMuons[0]->charge()*isoMuons[1]->charge();
00506       double mass = (isoMuons[0]->p4()+isoMuons[1]->p4()).mass();
00507       fill(charge<0 ? "invMass_"    : "invMassWC_"    , mass       );
00508       fill(charge<0 ? "invMassLog_" : "invMassWCLog_" , log10(mass));
00509       if((lowerEdge_==-1. && upperEdge_==-1.) || (lowerEdge_<mass && mass<upperEdge_) ){
00510         fill("dEtaL1L2"  , isoMuons[0]->eta()-isoMuons[1]->eta() );
00511         fill("sumEtaL1L2", (isoMuons[0]->eta()+isoMuons[1]->eta())/2);
00512         fill("dPhiL1L2"  , reco::deltaPhi(isoMuons[0]->phi(),isoMuons[1]->phi()) );
00513         fill("muonPt_", isoMuons[0]->pt()); fill("muonPt_", isoMuons[1]->pt()); 
00514         fill("lep1Pt_", isoMuons[0]->pt()); fill("lep2Pt_", isoMuons[1]->pt()); 
00515         // fill plots for trigger monitoring
00516         if(!triggerTable_.label().empty()) fill(event, *triggerTable, "diMuon", diMuonPaths_);
00517         if(diMuonLogged_<=hists_.find("diMuonLogger_")->second->getNbinsY()){
00518           // log runnumber, lumi block, event number & some
00519           // more pysics infomation for interesting events
00520           fill("diMuonLogger_", 0.5, diMuonLogged_+0.5, event.eventAuxiliary().run()); 
00521           fill("diMuonLogger_", 1.5, diMuonLogged_+0.5, event.eventAuxiliary().luminosityBlock()); 
00522           fill("diMuonLogger_", 2.5, diMuonLogged_+0.5, event.eventAuxiliary().event()); 
00523           fill("diMuonLogger_", 3.5, diMuonLogged_+0.5, isoMuons[0]->pt()); 
00524           fill("diMuonLogger_", 4.5, diMuonLogged_+0.5, isoMuons[1]->pt()); 
00525           if(leadingJets.size()>0) fill("diMuonLogger_", 5.5, diMuonLogged_+0.5, leadingJets[0].pt()); 
00526           if(leadingJets.size()>1) fill("diMuonLogger_", 6.5, diMuonLogged_+0.5, leadingJets[1].pt()); 
00527           fill("diMuonLogger_", 7.5, diMuonLogged_+0.5, caloMET.et()); 
00528           ++diMuonLogged_; 
00529         }
00530       }
00531     }
00532 
00533     // DIELEC channel
00534     if( decayChannel(isoMuons, isoElecs) == DIELEC ){
00535       int charge = isoElecs[0]->charge()*isoElecs[1]->charge();
00536       double mass = (isoElecs[0]->p4()+isoElecs[1]->p4()).mass();
00537       fill("decayChannel_", 2.5);
00538       fill(charge<0 ? "invMass_"    : "invMassWC_"    , mass       );
00539       fill(charge<0 ? "invMassLog_" : "invMassWCLog_" , log10(mass));
00540       if((lowerEdge_==-1. && upperEdge_==-1.) || (lowerEdge_<mass && mass<upperEdge_) ){
00541         fill("dEtaL1L2"  , isoElecs[0]->eta()-isoElecs[1]->eta() );
00542         fill("sumEtaL1L2", (isoElecs[0]->eta()+isoElecs[1]->eta())/2);
00543         fill("dPhiL1L2"  , reco::deltaPhi(isoElecs[0]->phi(),isoElecs[1]->phi()) );
00544         fill("elecPt_", isoElecs[0]->pt()); fill("elecPt_", isoElecs[1]->pt()); 
00545         fill("lep1Pt_", isoElecs[0]->pt()); fill("lep2Pt_", isoElecs[1]->pt()); 
00546         if(diElecLogged_<=hists_.find("diElecLogger_")->second->getNbinsY()){
00547           // log runnumber, lumi block, event number & some
00548           // more pysics infomation for interesting events
00549           fill("diElecLogger_", 0.5, diElecLogged_+0.5, event.eventAuxiliary().run()); 
00550           fill("diElecLogger_", 1.5, diElecLogged_+0.5, event.eventAuxiliary().luminosityBlock()); 
00551           fill("diElecLogger_", 2.5, diElecLogged_+0.5, event.eventAuxiliary().event()); 
00552           fill("diElecLogger_", 3.5, diElecLogged_+0.5, isoElecs[0]->pt()); 
00553           fill("diElecLogger_", 4.5, diElecLogged_+0.5, isoElecs[1]->pt()); 
00554           if(leadingJets.size()>0) fill("diElecLogger_", 5.5, diElecLogged_+0.5, leadingJets[0].pt()); 
00555           if(leadingJets.size()>1) fill("diElecLogger_", 6.5, diElecLogged_+0.5, leadingJets[1].pt()); 
00556           fill("diElecLogger_", 7.5, diElecLogged_+0.5, caloMET.et()); 
00557           ++diElecLogged_; 
00558         }
00559       }
00560    }
00561   }
00562   
00563 }
00564 
00565 TopDiLeptonOfflineDQM::TopDiLeptonOfflineDQM(const edm::ParameterSet& cfg): triggerTable_(""), vertex_("") 
00566 {
00567   // configure the preselection
00568   edm::ParameterSet presel=cfg.getParameter<edm::ParameterSet>("preselection");
00569   if( presel.existsAs<edm::ParameterSet>("trigger") ){
00570     edm::ParameterSet trigger=presel.getParameter<edm::ParameterSet>("trigger");
00571     triggerTable_=trigger.getParameter<edm::InputTag>("src");
00572     triggerPaths_=trigger.getParameter<std::vector<std::string> >("select");
00573   } 
00574   if( presel.existsAs<edm::ParameterSet>("vertex" ) ){
00575     edm::ParameterSet vertex=presel.getParameter<edm::ParameterSet>("vertex");
00576     vertex_= vertex.getParameter<edm::InputTag>("src");
00577     vertexSelect_= new StringCutObjectSelector<reco::Vertex>(vertex.getParameter<std::string>("select"));
00578   }
00579   if( presel.existsAs<edm::ParameterSet>("beamspot" ) ){
00580     edm::ParameterSet beamspot=presel.getParameter<edm::ParameterSet>("beamspot");
00581     beamspot_= beamspot.getParameter<edm::InputTag>("src");
00582     beamspotSelect_= new StringCutObjectSelector<reco::BeamSpot>(beamspot.getParameter<std::string>("select"));
00583   }
00584 
00585   // conifgure the selection
00586   std::vector<edm::ParameterSet> sel=cfg.getParameter<std::vector<edm::ParameterSet> >("selection");
00587   for(unsigned int i=0; i<sel.size(); ++i){
00588     selectionOrder_.push_back(sel.at(i).getParameter<std::string>("label"));
00589     selection_[selectionStep(selectionOrder_.back())] = std::make_pair(sel.at(i), new TopDiLeptonOffline::MonitorEnsemble(selectionStep(selectionOrder_.back()).c_str(), cfg.getParameter<edm::ParameterSet>("setup")));
00590   }
00591 }
00592 
00593 void 
00594 TopDiLeptonOfflineDQM::analyze(const edm::Event& event, const edm::EventSetup& setup)
00595 { 
00596   if(!triggerTable_.label().empty()){
00597     edm::Handle<edm::TriggerResults> triggerTable;
00598     if( !event.getByLabel(triggerTable_, triggerTable) ) return;
00599     if(!accept(event, *triggerTable, triggerPaths_)) return;
00600   }
00601   if(!vertex_.label().empty()){
00602     edm::Handle<std::vector<reco::Vertex> > vertex;
00603     if( !event.getByLabel(vertex_, vertex) ) return;
00604     if(vertex->empty() || !(*vertexSelect_)(vertex->front())) return;
00605   }
00606   if(!beamspot_.label().empty()){
00607     edm::Handle<reco::BeamSpot> beamspot;
00608     if( !event.getByLabel(beamspot_, beamspot) ) return;
00609     if(!(*beamspotSelect_)(*beamspot)) return;
00610   }
00611   // apply selection steps
00612   for(std::vector<std::string>::const_iterator selIt=selectionOrder_.begin(); selIt!=selectionOrder_.end(); ++selIt){
00613     std::string key = selectionStep(*selIt), type = objectType(*selIt);
00614     if(selection_.find(key)!=selection_.end()){
00615       if(type=="empty"){
00616         selection_[key].second->fill(event, setup);
00617       }
00618       if(type=="muons"){
00619         SelectionStep<reco::Muon> step(selection_[key].first);
00620         if(step.select(event)){
00621           selection_[key].second->fill(event, setup);
00622         } else break;
00623       }
00624       if(type=="elecs"){
00625         SelectionStep<reco::GsfElectron> step(selection_[key].first);
00626         if(step.select(event)){ 
00627           selection_[key].second->fill(event, setup);
00628         } else break;
00629       }
00630       if(type=="jets" ){
00631         SelectionStep<reco::Jet> step(selection_[key].first);
00632         if(step.select(event, setup)){
00633           selection_[key].second->fill(event, setup);
00634         } else break;
00635       }
00636       if(type=="jets/pf" ){
00637         SelectionStep<reco::PFJet> step(selection_[key].first);
00638         if(step.select(event, setup)){
00639           selection_[key].second->fill(event, setup);
00640         } else break;
00641       }
00642       if(type=="jets/calo" ){
00643         SelectionStep<reco::CaloJet> step(selection_[key].first);
00644         if(step.select(event, setup)){
00645           selection_[key].second->fill(event, setup);
00646         } else break;
00647       }
00648       if(type=="met" ){
00649         SelectionStep<reco::MET> step(selection_[key].first);
00650         if(step.select(event)){
00651           selection_[key].second->fill(event, setup);
00652         } else break;
00653       }
00654     }
00655   }
00656 }
00657 
00658