CMS 3D CMS Logo

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