Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2012/05/20 17:57:01 $
00005  *  $Revision: 1.51 $
00006  *  \author K. Hatakeyama - Rockefeller University
00007  *          A.Apresyan - Caltech
00008  */
00009 
00010 #include "DQMOffline/JetMET/interface/PFMETAnalyzer.h"
00011 #include "DataFormats/Common/interface/Handle.h"
00012 #include "FWCore/Common/interface/TriggerNames.h"
00013 
00014 #include "DataFormats/JetReco/interface/PFJetCollection.h"
00015 #include "DataFormats/JetReco/interface/PFJet.h"
00016 #include "DataFormats/Math/interface/LorentzVector.h"
00017 
00018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00019 
00020 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00021 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
00022 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00023 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00024 
00025 #include "DataFormats/Math/interface/LorentzVector.h"
00026 
00027 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
00028 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
00029 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00030 
00031 #include <string>
00032 using namespace edm;
00033 using namespace reco;
00034 using namespace math;
00035 
00036 // ***********************************************************
00037 PFMETAnalyzer::PFMETAnalyzer(const edm::ParameterSet& pSet) {
00038 
00039   parameters = pSet;
00040 
00041   edm::ParameterSet highptjetparms = parameters.getParameter<edm::ParameterSet>("highPtJetTrigger");
00042   edm::ParameterSet lowptjetparms  = parameters.getParameter<edm::ParameterSet>("lowPtJetTrigger" );
00043   edm::ParameterSet minbiasparms   = parameters.getParameter<edm::ParameterSet>("minBiasTrigger"  );
00044   edm::ParameterSet highmetparms   = parameters.getParameter<edm::ParameterSet>("highMETTrigger"  );
00045   //  edm::ParameterSet lowmetparms    = parameters.getParameter<edm::ParameterSet>("lowMETTrigger"   );
00046   edm::ParameterSet eleparms       = parameters.getParameter<edm::ParameterSet>("eleTrigger"      );
00047   edm::ParameterSet muonparms      = parameters.getParameter<edm::ParameterSet>("muonTrigger"     );
00048 
00049   //genericTriggerEventFlag_( new GenericTriggerEventFlag( conf_ ) );
00050   _HighPtJetEventFlag = new GenericTriggerEventFlag( highptjetparms );
00051   _LowPtJetEventFlag  = new GenericTriggerEventFlag( lowptjetparms  );
00052   _MinBiasEventFlag   = new GenericTriggerEventFlag( minbiasparms   );
00053   _HighMETEventFlag   = new GenericTriggerEventFlag( highmetparms   );
00054   //  _LowMETEventFlag    = new GenericTriggerEventFlag( lowmetparms    );
00055   _EleEventFlag       = new GenericTriggerEventFlag( eleparms       );
00056   _MuonEventFlag      = new GenericTriggerEventFlag( muonparms      );
00057 
00058   highPtJetExpr_ = highptjetparms.getParameter<std::vector<std::string> >("hltPaths");
00059   lowPtJetExpr_  = lowptjetparms .getParameter<std::vector<std::string> >("hltPaths");
00060   highMETExpr_   = highmetparms  .getParameter<std::vector<std::string> >("hltPaths");
00061   //  lowMETExpr_    = lowmetparms   .getParameter<std::vector<std::string> >("hltPaths");
00062   muonExpr_      = muonparms     .getParameter<std::vector<std::string> >("hltPaths");
00063   elecExpr_      = eleparms      .getParameter<std::vector<std::string> >("hltPaths");
00064   minbiasExpr_   = minbiasparms  .getParameter<std::vector<std::string> >("hltPaths");
00065 
00066 }
00067 
00068 // ***********************************************************
00069 PFMETAnalyzer::~PFMETAnalyzer() { 
00070 
00071   delete _HighPtJetEventFlag;
00072   delete _LowPtJetEventFlag;
00073   delete _MinBiasEventFlag;
00074   delete _HighMETEventFlag;
00075   //  delete _LowMETEventFlag;
00076   delete _EleEventFlag;
00077   delete _MuonEventFlag;
00078 
00079 }
00080 
00081 void PFMETAnalyzer::beginJob(DQMStore * dbe) {
00082 
00083   evtCounter = 0;
00084   metname = "pfMETAnalyzer";
00085 
00086   // trigger information
00087   HLTPathsJetMBByName_ = parameters.getParameter<std::vector<std::string > >("HLTPathsJetMB");
00088 
00089   theCleaningParameters = parameters.getParameter<ParameterSet>("CleaningParameters"),
00090 
00091   //Trigger parameters
00092   gtTag          = theCleaningParameters.getParameter<edm::InputTag>("gtLabel");
00093   _techTrigsAND  = theCleaningParameters.getParameter<std::vector<unsigned > >("techTrigsAND");
00094   _techTrigsOR   = theCleaningParameters.getParameter<std::vector<unsigned > >("techTrigsOR");
00095   _techTrigsNOT  = theCleaningParameters.getParameter<std::vector<unsigned > >("techTrigsNOT");
00096 
00097   _doHLTPhysicsOn = theCleaningParameters.getParameter<bool>("doHLTPhysicsOn");
00098   _hlt_PhysDec    = theCleaningParameters.getParameter<std::string>("HLT_PhysDec");
00099 
00100   _tightBHFiltering    = theCleaningParameters.getParameter<bool>("tightBHFiltering");
00101   _tightJetIDFiltering = theCleaningParameters.getParameter<int>("tightJetIDFiltering");
00102 
00103 
00104   // ==========================================================
00105   //DCS information
00106   // ==========================================================
00107   DCSFilter = new JetMETDQMDCSFilter(parameters.getParameter<ParameterSet>("DCSFilter"));
00108 
00109   //Vertex requirements
00110   _doPVCheck          = theCleaningParameters.getParameter<bool>("doPrimaryVertexCheck");
00111   vertexTag  = theCleaningParameters.getParameter<edm::InputTag>("vertexLabel");
00112 
00113   if (_doPVCheck) {
00114     _nvtx_min        = theCleaningParameters.getParameter<int>("nvtx_min");
00115     _nvtxtrks_min    = theCleaningParameters.getParameter<int>("nvtxtrks_min");
00116     _vtxndof_min     = theCleaningParameters.getParameter<int>("vtxndof_min");
00117     _vtxchi2_max     = theCleaningParameters.getParameter<double>("vtxchi2_max");
00118     _vtxz_max        = theCleaningParameters.getParameter<double>("vtxz_max");
00119   }
00120 
00121 
00122   // PFMET information
00123   thePfMETCollectionLabel       = parameters.getParameter<edm::InputTag>("METCollectionLabel");
00124   thePfJetCollectionLabel       = parameters.getParameter<edm::InputTag>("PfJetCollectionLabel");
00125   _source                       = parameters.getParameter<std::string>("Source");
00126 
00127   // Other data collections
00128   theJetCollectionLabel       = parameters.getParameter<edm::InputTag>("JetCollectionLabel");
00129   PFCandidatesTag             = parameters.getParameter<edm::InputTag>("PFCandidates");
00130   HcalNoiseRBXCollectionTag   = parameters.getParameter<edm::InputTag>("HcalNoiseRBXCollection");
00131   BeamHaloSummaryTag          = parameters.getParameter<edm::InputTag>("BeamHaloSummaryLabel");
00132   HBHENoiseFilterResultTag    = parameters.getParameter<edm::InputTag>("HBHENoiseFilterResultLabel");
00133 
00134   // misc
00135   _verbose     = parameters.getParameter<int>("verbose");
00136   _etThreshold = parameters.getParameter<double>("etThreshold"); // MET threshold
00137   _allhist     = parameters.getParameter<bool>("allHist");       // Full set of monitoring histograms
00138   _allSelection= parameters.getParameter<bool>("allSelection");  // Plot with all sets of event selection
00139   _cleanupSelection= parameters.getParameter<bool>("cleanupSelection");  // Plot with all sets of event selection
00140 
00141   _highPtPFJetThreshold = parameters.getParameter<double>("HighPtJetThreshold");
00142   _lowPtPFJetThreshold  = parameters.getParameter<double>("LowPtJetThreshold");
00143   _highPFMETThreshold   = parameters.getParameter<double>("HighMETThreshold");
00144 
00145   //
00146   jetID = new reco::helper::JetIDHelper(parameters.getParameter<ParameterSet>("JetIDParams"));
00147 
00148   // DQStore stuff
00149   LogTrace(metname)<<"[PFMETAnalyzer] Parameters initialization";
00150   std::string DirName = "JetMET/MET/"+_source;
00151   dbe->setCurrentFolder(DirName);
00152 
00153   metME = dbe->book1D("metReco", "metReco", 4, 1, 5);
00154   metME->setBinLabel(3,"PFMET",1);
00155 
00156   _dbe = dbe;
00157 
00158   _FolderNames.push_back("All");
00159   _FolderNames.push_back("BasicCleanup");
00160   _FolderNames.push_back("ExtraCleanup");
00161   _FolderNames.push_back("HcalNoiseFilter");
00162   _FolderNames.push_back("JetIDMinimal");
00163   _FolderNames.push_back("JetIDLoose");
00164   _FolderNames.push_back("JetIDTight");
00165   _FolderNames.push_back("BeamHaloIDTightPass");
00166   _FolderNames.push_back("BeamHaloIDLoosePass");
00167   _FolderNames.push_back("Triggers");
00168   _FolderNames.push_back("PV");
00169 
00170   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); 
00171        ic != _FolderNames.end(); ic++){
00172     if (*ic=="All")             bookMESet(DirName+"/"+*ic);
00173     if (_cleanupSelection){
00174     if (*ic=="BasicCleanup")    bookMESet(DirName+"/"+*ic);
00175     if (*ic=="ExtraCleanup")    bookMESet(DirName+"/"+*ic);
00176     }
00177     if (_allSelection){
00178       if (*ic=="HcalNoiseFilter")      bookMESet(DirName+"/"+*ic);
00179       if (*ic=="JetIDMinimal")         bookMESet(DirName+"/"+*ic);
00180       if (*ic=="JetIDLoose")           bookMESet(DirName+"/"+*ic);
00181       if (*ic=="JetIDTight")           bookMESet(DirName+"/"+*ic);
00182       if (*ic=="BeamHaloIDTightPass")  bookMESet(DirName+"/"+*ic);
00183       if (*ic=="BeamHaloIDLoosePass")  bookMESet(DirName+"/"+*ic);
00184       if (*ic=="Triggers")             bookMESet(DirName+"/"+*ic);
00185       if (*ic=="PV")                   bookMESet(DirName+"/"+*ic);
00186     }
00187   }
00188 }
00189 
00190 // ***********************************************************
00191 void PFMETAnalyzer::endJob() {
00192 
00193   delete jetID;
00194   delete DCSFilter;
00195 
00196 }
00197 
00198 // ***********************************************************
00199 void PFMETAnalyzer::bookMESet(std::string DirName)
00200 {
00201 
00202   bool bLumiSecPlot=false;
00203   if (DirName.find("All")!=std::string::npos) bLumiSecPlot=true;
00204 
00205   bookMonitorElement(DirName,bLumiSecPlot);
00206 
00207   if ( _HighPtJetEventFlag->on() ) {
00208     bookMonitorElement(DirName+"/"+"HighPtJet",false);
00209     meTriggerName_HighPtJet = _dbe->bookString("triggerName_HighPtJet", highPtJetExpr_[0]);
00210   }  
00211 
00212   if ( _LowPtJetEventFlag->on() ) {
00213     bookMonitorElement(DirName+"/"+"LowPtJet",false);
00214     meTriggerName_LowPtJet = _dbe->bookString("triggerName_LowPtJet", lowPtJetExpr_[0]);
00215   }
00216 
00217   if ( _MinBiasEventFlag->on() ) {
00218     bookMonitorElement(DirName+"/"+"MinBias",false);
00219     meTriggerName_MinBias = _dbe->bookString("triggerName_MinBias", minbiasExpr_[0]);
00220     if (_verbose) std::cout << "_MinBiasEventFlag is on, folder created\n";
00221   }
00222 
00223   if ( _HighMETEventFlag->on() ) {
00224     bookMonitorElement(DirName+"/"+"HighMET",false);
00225     meTriggerName_HighMET = _dbe->bookString("triggerName_HighMET", highMETExpr_[0]);
00226   }
00227 
00228   if ( _EleEventFlag->on() ) {
00229     bookMonitorElement(DirName+"/"+"Ele",false);
00230     meTriggerName_Ele = _dbe->bookString("triggerName_Ele", elecExpr_[0]);
00231     if (_verbose) std::cout << "_EleEventFlag is on, folder created\n";
00232   }
00233 
00234   if ( _MuonEventFlag->on() ) {
00235     bookMonitorElement(DirName+"/"+"Muon",false);
00236     meTriggerName_Muon = _dbe->bookString("triggerName_Muon", muonExpr_[0]);
00237     if (_verbose) std::cout << "_MuonEventFlag is on, folder created\n";
00238   }
00239 }
00240 
00241 
00242 //------------------------------------------------------------------------------
00243 // bookMonitorElement
00244 //------------------------------------------------------------------------------
00245 void PFMETAnalyzer::bookMonitorElement(std::string DirName,
00246                                        bool        bLumiSecPlot = false)
00247 {
00248   _dbe->setCurrentFolder(DirName);
00249 
00250 
00251   mePfMEx        = _dbe->book1D("METTask_PfMEx",        "pfmet.px()",           200, -500,  500); 
00252   mePfMEy        = _dbe->book1D("METTask_PfMEy",        "pfmet.py()",           200, -500,  500); 
00253   mePfMET        = _dbe->book1D("METTask_PfMET",        "pfmet.pt()",           200,    0, 1000); 
00254   mePfSumET      = _dbe->book1D("METTask_PfSumET",      "pfmet.sumEt()",        400,    0, 4000); 
00255   mePfMETSig     = _dbe->book1D("METTask_PfMETSig",     "pfmet.mEtSig()",        51,    0,   51);
00256   mePfMETPhi     = _dbe->book1D("METTask_PfMETPhi",     "pfmet.phi()",           60, -3.2,  3.2);
00257   mePfMET_logx   = _dbe->book1D("METTask_PfMET_logx",   "log10(pfmet.pt())",     40,   -1,    7);
00258   mePfSumET_logx = _dbe->book1D("METTask_PfSumET_logx", "log10(pfmet.sumEt())",  40,   -1,    7);
00259 
00260 
00261   mePhotonEtFraction        = _dbe->book1D("METTask_PfPhotonEtFraction",        "pfmet.photonEtFraction()",         50, 0,    1);
00262   mePhotonEt                = _dbe->book1D("METTask_PfPhotonEt",                "pfmet.photonEt()",                100, 0, 1000);
00263   meNeutralHadronEtFraction = _dbe->book1D("METTask_PfNeutralHadronEtFraction", "pfmet.neutralHadronEtFraction()",  50, 0,    1);
00264   meNeutralHadronEt         = _dbe->book1D("METTask_PfNeutralHadronEt",         "pfmet.neutralHadronEt()",         100, 0, 1000);
00265   meElectronEtFraction      = _dbe->book1D("METTask_PfElectronEtFraction",      "pfmet.electronEtFraction()",       50, 0,    1);
00266   meElectronEt              = _dbe->book1D("METTask_PfElectronEt",              "pfmet.electronEt()",              100, 0, 1000);
00267   meChargedHadronEtFraction = _dbe->book1D("METTask_PfChargedHadronEtFraction", "pfmet.chargedHadronEtFraction()",  50, 0,    1);
00268   meChargedHadronEt         = _dbe->book1D("METTask_PfChargedHadronEt",         "pfmet.chargedHadronEt()",         100, 0, 1000);
00269   meMuonEtFraction          = _dbe->book1D("METTask_PfMuonEtFraction",          "pfmet.muonEtFraction()",           50, 0,    1);
00270   meMuonEt                  = _dbe->book1D("METTask_PfMuonEt",                  "pfmet.muonEt()",                  100, 0, 1000);
00271   meHFHadronEtFraction      = _dbe->book1D("METTask_PfHFHadronEtFraction",      "pfmet.HFHadronEtFraction()",       50, 0,    1);    
00272   meHFHadronEt              = _dbe->book1D("METTask_PfHFHadronEt",              "pfmet.HFHadronEt()",              100, 0, 1000);
00273   meHFEMEtFraction          = _dbe->book1D("METTask_PfHFEMEtFraction",          "pfmet.HFEMEtFraction()",           50, 0,    1);
00274   meHFEMEt                  = _dbe->book1D("METTask_PfHFEMEt",                  "pfmet.HFEMEt()",                  100, 0, 1000);
00275 
00276 
00277   if (_allhist) {
00278     if (bLumiSecPlot) {
00279 
00280       mePfMExLS = _dbe->book2D("METTask_PfMEx_LS", "METTask_PfMEx_LS", 200, -200, 200, 50, 0, 500);
00281       mePfMEyLS = _dbe->book2D("METTask_PfMEy_LS", "METTask_PfMEy_LS", 200, -200, 200, 50, 0, 500);
00282 
00283       mePfMExLS->setAxisTitle("pfmet.px()", 1);
00284       mePfMEyLS->setAxisTitle("pfmet.px()", 1);
00285 
00286       mePfMExLS->setAxisTitle("event.luminosityBlock()", 2);
00287       mePfMEyLS->setAxisTitle("event.luminosityBlock()", 2);
00288     }
00289   }
00290 
00291 
00292   // Book NPV profiles
00293   //----------------------------------------------------------------------------
00294   mePfMEx_profile   = _dbe->bookProfile("METTask_PfMEx_profile",   "pfmet.px()",    nbinsPV, PVlow, PVup, 200, -500,  500);
00295   mePfMEy_profile   = _dbe->bookProfile("METTask_PfMEy_profile",   "pfmet.py()",    nbinsPV, PVlow, PVup, 200, -500,  500); 
00296   mePfMET_profile   = _dbe->bookProfile("METTask_PfMET_profile",   "pfmet.pt()",    nbinsPV, PVlow, PVup, 200,    0, 1000); 
00297   mePfSumET_profile = _dbe->bookProfile("METTask_PfSumET_profile", "pfmet.sumEt()", nbinsPV, PVlow, PVup, 400,    0, 4000); 
00298 
00299   mePhotonEtFraction_profile        = _dbe->bookProfile("METTask_PfPhotonEtFraction_profile",        "pfmet.photonEtFraction()",        nbinsPV, PVlow, PVup,  50, 0,    1);
00300   mePhotonEt_profile                = _dbe->bookProfile("METTask_PfPhotonEt_profile",                "pfmet.photonEt()",                nbinsPV, PVlow, PVup, 100, 0, 1000);
00301   meNeutralHadronEtFraction_profile = _dbe->bookProfile("METTask_PfNeutralHadronEtFraction_profile", "pfmet.neutralHadronEtFraction()", nbinsPV, PVlow, PVup,  50, 0,    1);
00302   meNeutralHadronEt_profile         = _dbe->bookProfile("METTask_PfNeutralHadronEt_profile",         "pfmet.neutralHadronEt()",         nbinsPV, PVlow, PVup, 100, 0, 1000);
00303   meElectronEtFraction_profile      = _dbe->bookProfile("METTask_PfElectronEtFraction_profile",      "pfmet.electronEtFraction()",      nbinsPV, PVlow, PVup,  50, 0,    1);
00304   meElectronEt_profile              = _dbe->bookProfile("METTask_PfElectronEt_profile",              "pfmet.electronEt()",              nbinsPV, PVlow, PVup, 100, 0, 1000);
00305   meChargedHadronEtFraction_profile = _dbe->bookProfile("METTask_PfChargedHadronEtFraction_profile", "pfmet.chargedHadronEtFraction()", nbinsPV, PVlow, PVup,  50, 0,    1);
00306   meChargedHadronEt_profile         = _dbe->bookProfile("METTask_PfChargedHadronEt_profile",         "pfmet.chargedHadronEt()",         nbinsPV, PVlow, PVup, 100, 0, 1000);
00307   meMuonEtFraction_profile          = _dbe->bookProfile("METTask_PfMuonEtFraction_profile",          "pfmet.muonEtFraction()",          nbinsPV, PVlow, PVup,  50, 0,    1);
00308   meMuonEt_profile                  = _dbe->bookProfile("METTask_PfMuonEt_profile",                  "pfmet.muonEt()",                  nbinsPV, PVlow, PVup, 100, 0, 1000);
00309   meHFHadronEtFraction_profile      = _dbe->bookProfile("METTask_PfHFHadronEtFraction_profile",      "pfmet.HFHadronEtFraction()",      nbinsPV, PVlow, PVup,  50, 0,    1);    
00310   meHFHadronEt_profile              = _dbe->bookProfile("METTask_PfHFHadronEt_profile",              "pfmet.HFHadronEt()",              nbinsPV, PVlow, PVup, 100, 0, 1000);
00311   meHFEMEtFraction_profile          = _dbe->bookProfile("METTask_PfHFEMEtFraction_profile",          "pfmet.HFEMEtFraction()",          nbinsPV, PVlow, PVup,  50, 0,    1);
00312   meHFEMEt_profile                  = _dbe->bookProfile("METTask_PfHFEMEt_profile",                  "pfmet.HFEMEt()",                  nbinsPV, PVlow, PVup, 100, 0, 1000);
00313 
00314 
00315   // Set NPV profiles x-axis title
00316   //----------------------------------------------------------------------------
00317   mePfMEx_profile  ->setAxisTitle("nvtx", 1);
00318   mePfMEy_profile  ->setAxisTitle("nvtx", 1);
00319   mePfMET_profile  ->setAxisTitle("nvtx", 1);
00320   mePfSumET_profile->setAxisTitle("nvtx", 1);
00321 
00322   mePhotonEtFraction_profile       ->setAxisTitle("nvtx", 1);
00323   mePhotonEt_profile               ->setAxisTitle("nvtx", 1);
00324   meNeutralHadronEtFraction_profile->setAxisTitle("nvtx", 1);
00325   meNeutralHadronEt_profile        ->setAxisTitle("nvtx", 1);
00326   meElectronEtFraction_profile     ->setAxisTitle("nvtx", 1);
00327   meElectronEt_profile             ->setAxisTitle("nvtx", 1);
00328   meChargedHadronEtFraction_profile->setAxisTitle("nvtx", 1);
00329   meChargedHadronEt_profile        ->setAxisTitle("nvtx", 1);
00330   meMuonEtFraction_profile         ->setAxisTitle("nvtx", 1);
00331   meMuonEt_profile                 ->setAxisTitle("nvtx", 1);
00332   meHFHadronEtFraction_profile     ->setAxisTitle("nvtx", 1);    
00333   meHFHadronEt_profile             ->setAxisTitle("nvtx", 1);
00334   meHFEMEtFraction_profile         ->setAxisTitle("nvtx", 1);
00335   meHFEMEt_profile                 ->setAxisTitle("nvtx", 1);
00336 }
00337 
00338 
00339 //------------------------------------------------------------------------------
00340 // beginRun
00341 //------------------------------------------------------------------------------
00342 void PFMETAnalyzer::beginRun(const edm::Run&        iRun,
00343                              const edm::EventSetup& iSetup)
00344 {
00345   if (_HighPtJetEventFlag->on()) _HighPtJetEventFlag->initRun(iRun, iSetup);
00346   if (_LowPtJetEventFlag ->on()) _LowPtJetEventFlag ->initRun(iRun, iSetup);
00347   if (_MinBiasEventFlag  ->on()) _MinBiasEventFlag  ->initRun(iRun, iSetup);
00348   if (_HighMETEventFlag  ->on()) _HighMETEventFlag  ->initRun(iRun, iSetup);
00349   if (_EleEventFlag      ->on()) _EleEventFlag      ->initRun(iRun, iSetup);
00350   if (_MuonEventFlag     ->on()) _MuonEventFlag     ->initRun(iRun, iSetup);
00351 
00352   if (_HighPtJetEventFlag->on() && _HighPtJetEventFlag->expressionsFromDB(_HighPtJetEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00353     highPtJetExpr_ = _HighPtJetEventFlag->expressionsFromDB(_HighPtJetEventFlag->hltDBKey(), iSetup);
00354 
00355   if (_LowPtJetEventFlag->on() && _LowPtJetEventFlag->expressionsFromDB(_LowPtJetEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00356     lowPtJetExpr_ = _LowPtJetEventFlag->expressionsFromDB(_LowPtJetEventFlag->hltDBKey(), iSetup);
00357 
00358   if (_HighMETEventFlag->on() && _HighMETEventFlag->expressionsFromDB(_HighMETEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00359     highMETExpr_ = _HighMETEventFlag->expressionsFromDB(_HighMETEventFlag->hltDBKey(), iSetup);
00360 
00361   if (_MuonEventFlag->on() && _MuonEventFlag->expressionsFromDB(_MuonEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00362     muonExpr_ = _MuonEventFlag->expressionsFromDB(_MuonEventFlag->hltDBKey(), iSetup);
00363 
00364   if (_EleEventFlag->on() && _EleEventFlag->expressionsFromDB(_EleEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00365     elecExpr_ = _EleEventFlag->expressionsFromDB(_EleEventFlag->hltDBKey(), iSetup);
00366 
00367   if (_MinBiasEventFlag->on() && _MinBiasEventFlag->expressionsFromDB(_MinBiasEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00368     minbiasExpr_ = _MinBiasEventFlag->expressionsFromDB(_MinBiasEventFlag->hltDBKey(), iSetup);
00369 }
00370 
00371 
00372 //------------------------------------------------------------------------------
00373 // endRun
00374 //------------------------------------------------------------------------------
00375 void PFMETAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup, DQMStore * dbe)
00376 {
00377   
00378   //
00379   //--- Check the time length of the Run from the lumi section plots
00380 
00381   std::string dirName = "JetMET/MET/"+_source+"/";
00382   _dbe->setCurrentFolder(dirName);
00383 
00384   TH1F* tlumisec;
00385 
00386   MonitorElement *meLumiSec = _dbe->get("aaa");
00387   meLumiSec = _dbe->get("JetMET/lumisec");
00388 
00389   int totlsec=0;
00390   double totltime=0.;
00391   if ( meLumiSec->getRootObject() ) {
00392     tlumisec = meLumiSec->getTH1F();
00393     for (int i=0; i<500; i++){
00394       if (tlumisec->GetBinContent(i+1)) totlsec++;
00395     }
00396     totltime = double(totlsec*90); // one lumi sec ~ 90 (sec)
00397   }
00398 
00399   if (totltime==0.) totltime=1.; 
00400 
00401   //
00402   //--- Make the integrated plots with rate (Hz)
00403 
00404   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); ic != _FolderNames.end(); ic++)
00405     {
00406 
00407       std::string DirName;
00408       DirName = dirName+*ic;
00409 
00410       makeRatePlot(DirName,totltime);
00411       if ( _HighPtJetEventFlag->on() ) 
00412         makeRatePlot(DirName+"/"+"triggerName_HighJetPt",totltime);
00413       if ( _LowPtJetEventFlag->on() ) 
00414         makeRatePlot(DirName+"/"+"triggerName_LowJetPt",totltime);
00415       if ( _MinBiasEventFlag->on() ) 
00416         makeRatePlot(DirName+"/"+"triggerName_MinBias",totltime);
00417       if ( _HighMETEventFlag->on() ) 
00418         makeRatePlot(DirName+"/"+"triggerName_HighMET",totltime);
00419       if ( _EleEventFlag->on() ) 
00420         makeRatePlot(DirName+"/"+"triggerName_Ele",totltime);
00421       if ( _MuonEventFlag->on() ) 
00422         makeRatePlot(DirName+"/"+"triggerName_Muon",totltime);
00423     }
00424 }
00425 
00426 
00427 // ***********************************************************
00428 void PFMETAnalyzer::makeRatePlot(std::string DirName, double totltime)
00429 {
00430 
00431   _dbe->setCurrentFolder(DirName);
00432   MonitorElement *mePfMET = _dbe->get(DirName+"/"+"METTask_PfMET");
00433 
00434   TH1F* tPfMET;
00435   TH1F* tPfMETRate;
00436 
00437   if ( mePfMET )
00438     if ( mePfMET->getRootObject() ) {
00439       tPfMET     = mePfMET->getTH1F();
00440       
00441       // Integral plot & convert number of events to rate (hz)
00442       tPfMETRate = (TH1F*) tPfMET->Clone("METTask_PfMETRate");
00443       for (int i = tPfMETRate->GetNbinsX()-1; i>=0; i--){
00444         tPfMETRate->SetBinContent(i+1,tPfMETRate->GetBinContent(i+2)+tPfMET->GetBinContent(i+1));
00445       }
00446       for (int i = 0; i<tPfMETRate->GetNbinsX(); i++){
00447         tPfMETRate->SetBinContent(i+1,tPfMETRate->GetBinContent(i+1)/double(totltime));
00448       }
00449 
00450       tPfMETRate->SetName("METTask_PfMETRate");
00451       tPfMETRate->SetTitle("METTask_PfMETRate");
00452       mePfMETRate      = _dbe->book1D("METTask_PfMETRate",tPfMETRate);
00453       
00454     }
00455 }
00456 
00457 // ***********************************************************
00458 void PFMETAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
00459                             const edm::TriggerResults& triggerResults) {
00460 
00461   if (_verbose) std::cout << "PfMETAnalyzer analyze" << std::endl;
00462 
00463   LogTrace(metname)<<"[PFMETAnalyzer] Analyze PFMET";
00464 
00465   metME->Fill(3);
00466 
00467   // ==========================================================  
00468   // Trigger information 
00469   //
00470   _trig_JetMB=0;
00471   _trig_HighPtJet=0;
00472   _trig_LowPtJet=0;
00473   _trig_MinBias=0;
00474   _trig_HighMET=0;
00475   _trig_Ele=0;
00476   _trig_Muon=0;
00477   _trig_PhysDec=0;
00478   if(&triggerResults) {   
00479 
00481 
00482     //
00483     //
00484     // Check how many HLT triggers are in triggerResults 
00485     int ntrigs = triggerResults.size();
00486     if (_verbose) std::cout << "ntrigs=" << ntrigs << std::endl;
00487     
00488     //
00489     //
00490     // If index=ntrigs, this HLT trigger doesn't exist in the HLT table for this data.
00491     const edm::TriggerNames & triggerNames = iEvent.triggerNames(triggerResults);
00492     
00493     //
00494     //
00495     const unsigned int nTrig(triggerNames.size());
00496     for (unsigned int i=0;i<nTrig;++i)
00497       {
00498         if (triggerNames.triggerName(i).find(highPtJetExpr_[0].substr(0,highPtJetExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00499           _trig_HighPtJet=true;
00500         else if (triggerNames.triggerName(i).find(lowPtJetExpr_[0].substr(0,lowPtJetExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00501           _trig_LowPtJet=true;
00502         else if (triggerNames.triggerName(i).find(highMETExpr_[0].substr(0,highMETExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00503           _trig_HighMET=true;
00504         //        else if (triggerNames.triggerName(i).find(lowMETExpr_[0].substr(0,lowMETExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00505         //        _trig_LowMET=true;
00506         else if (triggerNames.triggerName(i).find(muonExpr_[0].substr(0,muonExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00507           _trig_Muon=true;
00508         else if (triggerNames.triggerName(i).find(elecExpr_[0].substr(0,elecExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00509           _trig_Ele=true;
00510         else if (triggerNames.triggerName(i).find(minbiasExpr_[0].substr(0,minbiasExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00511           _trig_MinBias=true;
00512       }
00513 
00514       // count number of requested Jet or MB HLT paths which have fired
00515     for (unsigned int i=0; i!=HLTPathsJetMBByName_.size(); i++) {
00516       unsigned int triggerIndex = triggerNames.triggerIndex(HLTPathsJetMBByName_[i]);
00517       if (triggerIndex<triggerResults.size()) {
00518         if (triggerResults.accept(triggerIndex)) {
00519           _trig_JetMB++;
00520         }
00521       }
00522     }
00523     // for empty input vectors (n==0), take all HLT triggers!
00524     if (HLTPathsJetMBByName_.size()==0) _trig_JetMB=triggerResults.size()-1;
00525 
00526     
00527     if (triggerNames.triggerIndex(_hlt_PhysDec)   != triggerNames.size() &&
00528         triggerResults.accept(triggerNames.triggerIndex(_hlt_PhysDec)))   _trig_PhysDec=1;
00529   } else {
00530 
00531     edm::LogInfo("PFMetAnalyzer") << "TriggerResults::HLT not found, "
00532       "automatically select events"; 
00533 
00534     // TriggerResults object not found. Look at all events.    
00535     _trig_JetMB=1;
00536   }
00537 
00538   // ==========================================================
00539   // PfMET information
00540   
00541   // **** Get the MET container  
00542   edm::Handle<reco::PFMETCollection> pfmetcoll;
00543   iEvent.getByLabel(thePfMETCollectionLabel, pfmetcoll);
00544   
00545   if(!pfmetcoll.isValid()) return;
00546 
00547   const PFMETCollection *pfmetcol = pfmetcoll.product();
00548   const PFMET *pfmet;
00549   pfmet = &(pfmetcol->front());
00550     
00551   LogTrace(metname)<<"[PfMETAnalyzer] Call to the PfMET analyzer";
00552 
00553   // ==========================================================
00554   //
00555   edm::Handle<HcalNoiseRBXCollection> HRBXCollection;
00556   iEvent.getByLabel(HcalNoiseRBXCollectionTag,HRBXCollection);
00557   if (!HRBXCollection.isValid()) {
00558     LogDebug("") << "PfMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00559     if (_verbose) std::cout << "PfMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00560   }
00561 
00562   
00563   edm::Handle<bool> HBHENoiseFilterResultHandle;
00564   iEvent.getByLabel(HBHENoiseFilterResultTag, HBHENoiseFilterResultHandle);
00565   bool HBHENoiseFilterResult = *HBHENoiseFilterResultHandle;
00566   if (!HBHENoiseFilterResultHandle.isValid()) {
00567     LogDebug("") << "PFMETAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
00568     if (_verbose) std::cout << "PFMETAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
00569   }
00570 
00571 
00572   edm::Handle<reco::CaloJetCollection> caloJets;
00573   iEvent.getByLabel(theJetCollectionLabel, caloJets);
00574   if (!caloJets.isValid()) {
00575     LogDebug("") << "PFMETAnalyzer: Could not find jet product" << std::endl;
00576     if (_verbose) std::cout << "PFMETAnalyzer: Could not find jet product" << std::endl;
00577   }
00578 
00579   edm::Handle<edm::View<PFCandidate> > pfCandidates;
00580   iEvent.getByLabel(PFCandidatesTag, pfCandidates);
00581   if (!pfCandidates.isValid()) {
00582     LogDebug("") << "PfMETAnalyzer: Could not find pfcandidates product" << std::endl;
00583     if (_verbose) std::cout << "PfMETAnalyzer: Could not find pfcandidates product" << std::endl;
00584   }
00585 
00586   edm::Handle<reco::PFJetCollection> pfJets;
00587   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
00588   if (!pfJets.isValid()) {
00589     LogDebug("") << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
00590     if (_verbose) std::cout << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
00591   }
00592   // ==========================================================
00593   // PfMET sanity check
00594 
00595   if (_source=="PfMET") validateMET(*pfmet, pfCandidates);
00596   
00597   // ==========================================================
00598   // JetID 
00599 
00600   if (_verbose) std::cout << "JetID starts" << std::endl;
00601   
00602   //
00603   // --- Minimal cuts
00604   //
00605   bool bJetIDMinimal=true;
00606   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00607        cal!=caloJets->end(); ++cal){
00608     jetID->calculate(iEvent, *cal);
00609     if (cal->pt()>10.){
00610       if (fabs(cal->eta())<=2.6 && 
00611           cal->emEnergyFraction()<=0.01) bJetIDMinimal=false;
00612     }
00613   }
00614 
00615   //
00616   // --- Loose cuts, not PF specific for now!
00617   //
00618   bool bJetIDLoose=true;
00619   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00620        cal!=caloJets->end(); ++cal){ 
00621     jetID->calculate(iEvent, *cal);
00622     if (_verbose) std::cout << jetID->n90Hits() << " " 
00623                             << jetID->restrictedEMF() << " "
00624                             << cal->pt() << std::endl;
00625     if (cal->pt()>10.){
00626       //
00627       // for all regions
00628       if (jetID->n90Hits()<2)  bJetIDLoose=false; 
00629       if (jetID->fHPD()>=0.98) bJetIDLoose=false; 
00630       //if (jetID->restrictedEMF()<0.01) bJetIDLoose=false; 
00631       //
00632       // for non-forward
00633       if (fabs(cal->eta())<2.55){
00634         if (cal->emEnergyFraction()<=0.01) bJetIDLoose=false; 
00635       }
00636       // for forward
00637       else {
00638         if (cal->emEnergyFraction()<=-0.9) bJetIDLoose=false; 
00639         if (cal->pt()>80.){
00640           if (cal->emEnergyFraction()>= 1.0) bJetIDLoose=false; 
00641         }
00642       } // forward vs non-forward
00643     }   // pt>10 GeV/c
00644   }     // calor-jets loop
00645  
00646   //
00647   // --- Tight cuts
00648   //
00649   bool bJetIDTight=true;
00650   bJetIDTight=bJetIDLoose;
00651   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00652        cal!=caloJets->end(); ++cal){
00653     jetID->calculate(iEvent, *cal);
00654     if (cal->pt()>25.){
00655       //
00656       // for all regions
00657       if (jetID->fHPD()>=0.95) bJetIDTight=false; 
00658       //
00659       // for 1.0<|eta|<1.75
00660       if (fabs(cal->eta())>=1.00 && fabs(cal->eta())<1.75){
00661         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00662       }
00663       //
00664       // for 1.75<|eta|<2.55
00665       else if (fabs(cal->eta())>=1.75 && fabs(cal->eta())<2.55){
00666         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00667       }
00668       //
00669       // for 2.55<|eta|<3.25
00670       else if (fabs(cal->eta())>=2.55 && fabs(cal->eta())<3.25){
00671         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3) bJetIDTight=false; 
00672         if (cal->pt()>=50. && cal->pt()< 80. && cal->emEnergyFraction()<=-0.2) bJetIDTight=false; 
00673         if (cal->pt()>=80. && cal->pt()<340. && cal->emEnergyFraction()<=-0.1) bJetIDTight=false; 
00674         if (cal->pt()>=340.                  && cal->emEnergyFraction()<=-0.1 
00675             && cal->emEnergyFraction()>=0.95) bJetIDTight=false; 
00676       }
00677       //
00678       // for 3.25<|eta|
00679       else if (fabs(cal->eta())>=3.25){
00680         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3
00681             && cal->emEnergyFraction()>=0.90) bJetIDTight=false; 
00682         if (cal->pt()>=50. && cal->pt()<130. && cal->emEnergyFraction()<=-0.2
00683             && cal->emEnergyFraction()>=0.80) bJetIDTight=false; 
00684         if (cal->pt()>=130.                  && cal->emEnergyFraction()<=-0.1 
00685             && cal->emEnergyFraction()>=0.70) bJetIDTight=false; 
00686       }
00687     }   // pt>10 GeV/c
00688   }     // calor-jets loop
00689   
00690   if (_verbose) std::cout << "JetID ends" << std::endl;
00691 
00692 
00693   // ==========================================================
00694   // HCAL Noise filter
00695   
00696   bool bHcalNoiseFilter = HBHENoiseFilterResult;
00697 
00698   // ==========================================================
00699   // Get BeamHaloSummary
00700   edm::Handle<BeamHaloSummary> TheBeamHaloSummary ;
00701   iEvent.getByLabel(BeamHaloSummaryTag, TheBeamHaloSummary) ;
00702 
00703   bool bBeamHaloIDTightPass = true;
00704   bool bBeamHaloIDLoosePass = true;
00705 
00706   if(!TheBeamHaloSummary.isValid()) {
00707 
00708     const BeamHaloSummary TheSummary = (*TheBeamHaloSummary.product() );
00709 
00710     if( !TheSummary.EcalLooseHaloId()  && !TheSummary.HcalLooseHaloId() && 
00711         !TheSummary.CSCLooseHaloId()   && !TheSummary.GlobalLooseHaloId() )
00712       bBeamHaloIDLoosePass = false;
00713 
00714     if( !TheSummary.EcalTightHaloId()  && !TheSummary.HcalTightHaloId() && 
00715         !TheSummary.CSCTightHaloId()   && !TheSummary.GlobalTightHaloId() )
00716       bBeamHaloIDTightPass = false;
00717 
00718   }
00719 
00720   // ==========================================================
00721   //Vertex information
00722   
00723   _numPV = 0;
00724   bool bPrimaryVertex = true;
00725   if(_doPVCheck){
00726     bPrimaryVertex = false;
00727     Handle<VertexCollection> vertexHandle;
00728 
00729     iEvent.getByLabel(vertexTag, vertexHandle);
00730 
00731     if (!vertexHandle.isValid()) {
00732       LogDebug("") << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00733       if (_verbose) std::cout << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00734     }
00735     
00736     if ( vertexHandle.isValid() ){
00737       VertexCollection vertexCollection = *(vertexHandle.product());
00738       int vertex_number     = vertexCollection.size();
00739       VertexCollection::const_iterator v = vertexCollection.begin();
00740       for ( ; v != vertexCollection.end(); ++v) {
00741         double vertex_chi2    = v->normalizedChi2();
00742         double vertex_ndof    = v->ndof();
00743         bool   fakeVtx        = v->isFake();
00744         double vertex_Z       = v->z();
00745         
00746         if (  !fakeVtx
00747               && vertex_number>=_nvtx_min
00748               && vertex_ndof   >_vtxndof_min
00749               && vertex_chi2   <_vtxchi2_max
00750               && fabs(vertex_Z)<_vtxz_max ) {
00751           bPrimaryVertex = true;
00752           ++_numPV;
00753         }
00754       }
00755     }
00756   }
00757   // ==========================================================
00758 
00759   edm::Handle< L1GlobalTriggerReadoutRecord > gtReadoutRecord;
00760   iEvent.getByLabel( gtTag, gtReadoutRecord);
00761 
00762   if (!gtReadoutRecord.isValid()) {
00763     LogDebug("") << "CaloMETAnalyzer: Could not find GT readout record" << std::endl;
00764     if (_verbose) std::cout << "CaloMETAnalyzer: Could not find GT readout record product" << std::endl;
00765   }
00766   
00767   bool bTechTriggers    = true;
00768   bool bTechTriggersAND = true;
00769   bool bTechTriggersOR  = false;
00770   bool bTechTriggersNOT = false;
00771 
00772   if (gtReadoutRecord.isValid()) {
00773     const TechnicalTriggerWord&  technicalTriggerWordBeforeMask = gtReadoutRecord->technicalTriggerWord();
00774     
00775     if (_techTrigsAND.size() == 0)
00776       bTechTriggersAND = true;
00777     else
00778       for (unsigned ttr = 0; ttr != _techTrigsAND.size(); ttr++) {
00779         bTechTriggersAND = bTechTriggersAND && technicalTriggerWordBeforeMask.at(_techTrigsAND.at(ttr));
00780       }
00781     
00782     if (_techTrigsAND.size() == 0)
00783       bTechTriggersOR = true;
00784     else
00785       for (unsigned ttr = 0; ttr != _techTrigsOR.size(); ttr++) {
00786         bTechTriggersOR = bTechTriggersOR || technicalTriggerWordBeforeMask.at(_techTrigsOR.at(ttr));
00787       }
00788     if (_techTrigsNOT.size() == 0)
00789       bTechTriggersNOT = false;
00790     else
00791       for (unsigned ttr = 0; ttr != _techTrigsNOT.size(); ttr++) {
00792         bTechTriggersNOT = bTechTriggersNOT || technicalTriggerWordBeforeMask.at(_techTrigsNOT.at(ttr));
00793       }
00794   }
00795   else
00796     {
00797       bTechTriggersAND = true;
00798       bTechTriggersOR  = true;
00799       bTechTriggersNOT = false;
00800     }
00801   
00802   if (_techTrigsAND.size()==0)
00803     bTechTriggersAND = true;
00804   if (_techTrigsOR.size()==0)
00805     bTechTriggersOR = true;
00806   if (_techTrigsNOT.size()==0)
00807     bTechTriggersNOT = false;
00808   
00809   bTechTriggers = bTechTriggersAND && bTechTriggersOR && !bTechTriggersNOT;
00810 
00811   // ==========================================================
00812   // Reconstructed MET Information - fill MonitorElements
00813   
00814   bool bHcalNoise   = bHcalNoiseFilter;
00815   bool bBeamHaloID  = bBeamHaloIDLoosePass;
00816   bool bJetID       = bJetIDMinimal;
00817 
00818   bool bPhysicsDeclared = true;
00819   if(_doHLTPhysicsOn) bPhysicsDeclared =_trig_PhysDec;
00820 
00821   if      (_tightBHFiltering)       bBeamHaloID = bBeamHaloIDTightPass;
00822 
00823   if      (_tightJetIDFiltering==1)  bJetID      = bJetIDMinimal;
00824   else if (_tightJetIDFiltering==2)  bJetID      = bJetIDLoose;
00825   else if (_tightJetIDFiltering==3)  bJetID      = bJetIDTight;
00826   else if (_tightJetIDFiltering==-1) bJetID      = true;
00827 
00828   bool bBasicCleanup = bTechTriggers && bPrimaryVertex && bPhysicsDeclared;
00829   bool bExtraCleanup = bBasicCleanup && bHcalNoise && bJetID && bBeamHaloID;
00830   
00831   std::string DirName = "JetMET/MET/"+_source;
00832   
00833   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); 
00834        ic != _FolderNames.end(); ic++){
00835     if (*ic=="All")                                             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00836     if (DCSFilter->filter(iEvent, iSetup)) {
00837     if (_cleanupSelection){
00838     if (*ic=="BasicCleanup" && bBasicCleanup)                   fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00839     if (*ic=="ExtraCleanup" && bExtraCleanup)                   fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00840     }
00841     if (_allSelection) {
00842       if (*ic=="HcalNoiseFilter"      && bHcalNoiseFilter )       fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00843       if (*ic=="JetIDMinimal"         && bJetIDMinimal)           fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00844       if (*ic=="JetIDLoose"           && bJetIDLoose)             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00845       if (*ic=="JetIDTight"           && bJetIDTight)             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00846       if (*ic=="BeamHaloIDTightPass"  && bBeamHaloIDTightPass)    fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00847       if (*ic=="BeamHaloIDLoosePass"  && bBeamHaloIDLoosePass)    fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00848       if (*ic=="Triggers"             && bTechTriggers)           fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00849       if (*ic=="PV"                   && bPrimaryVertex)          fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00850     }
00851     } // DCS
00852   }
00853 }
00854 
00855   
00856 // ***********************************************************
00857 void PFMETAnalyzer::validateMET(const reco::PFMET& pfmet, 
00858                                 edm::Handle<edm::View<PFCandidate> > pfCandidates)
00859 {          
00860   double sumEx = 0;
00861   double sumEy = 0;
00862   double sumEt = 0;
00863       
00864   for( unsigned i=0; i<pfCandidates->size(); i++ ) {
00865         
00866     const reco::PFCandidate& cand = (*pfCandidates)[i];
00867         
00868     double E = cand.energy();
00869         
00871     //  if( cand.particleId()==PFCandidate::h_HF || 
00872     //      cand.particleId()==PFCandidate::egamma_HF )
00873     //    E *= hfCalibFactor_;
00874         
00875     double phi = cand.phi();
00876     double cosphi = cos(phi);
00877     double sinphi = sin(phi);
00878         
00879     double theta = cand.theta();
00880     double sintheta = sin(theta);
00881         
00882     double et = E*sintheta;
00883     double ex = et*cosphi;
00884     double ey = et*sinphi;
00885         
00886     sumEx += ex;
00887     sumEy += ey;
00888     sumEt += et;
00889   }
00890       
00891   double Et = sqrt( sumEx*sumEx + sumEy*sumEy);
00892   XYZTLorentzVector missingEt( -sumEx, -sumEy, 0, Et);
00893       
00894   if(_verbose) 
00895     if (sumEt!=pfmet.sumEt() || sumEx!=pfmet.px() || sumEy!=pfmet.py() || missingEt.T()!=pfmet.pt() )   
00896       {
00897         std::cout<<"PFSumEt: " << sumEt         <<", "<<"PFMETBlock: "<<pfmet.pt()<<std::endl;
00898         std::cout<<"PFMET: "   << missingEt.T() <<", "<<"PFMETBlock: "<<pfmet.pt()<<std::endl;
00899         std::cout<<"PFMETx: "  << missingEt.X() <<", "<<"PFMETBlockx: "<<pfmet.pt()<<std::endl;
00900         std::cout<<"PFMETy: "  << missingEt.Y() <<", "<<"PFMETBlocky: "<<pfmet.pt()<<std::endl;
00901       }
00902 }
00903 
00904 // ***********************************************************
00905 void PFMETAnalyzer::fillMESet(const edm::Event& iEvent, std::string DirName, 
00906                               const reco::PFMET& pfmet)
00907 {
00908 
00909   _dbe->setCurrentFolder(DirName);
00910 
00911   bool bLumiSecPlot=false;
00912   if (DirName.find("All")) bLumiSecPlot=true;
00913 
00914   if (_trig_JetMB)
00915     fillMonitorElement(iEvent,DirName,"",pfmet, bLumiSecPlot);
00916   if (_trig_HighPtJet)
00917     fillMonitorElement(iEvent,DirName,"HighPtJet",pfmet,false);
00918   if (_trig_LowPtJet)
00919     fillMonitorElement(iEvent,DirName,"LowPtJet",pfmet,false);
00920   if (_trig_MinBias)
00921     fillMonitorElement(iEvent,DirName,"MinBias",pfmet,false);
00922   if (_trig_HighMET)
00923     fillMonitorElement(iEvent,DirName,"HighMET",pfmet,false);
00924   if (_trig_Ele)
00925     fillMonitorElement(iEvent,DirName,"Ele",pfmet,false);
00926   if (_trig_Muon)
00927     fillMonitorElement(iEvent,DirName,"Muon",pfmet,false);
00928 }
00929 
00930 
00931 //------------------------------------------------------------------------------
00932 // fillMonitorElement
00933 //------------------------------------------------------------------------------
00934 void PFMETAnalyzer::fillMonitorElement(const edm::Event&  iEvent,
00935                                        std::string        DirName, 
00936                                        std::string        TriggerTypeName, 
00937                                        const reco::PFMET& pfmet,
00938                                        bool               bLumiSecPlot)
00939 {
00940   if (TriggerTypeName == "HighPtJet") {
00941     if (!selectHighPtJetEvent(iEvent)) return;
00942   }
00943   else if (TriggerTypeName == "LowPtJet") {
00944     if (!selectLowPtJetEvent(iEvent)) return;
00945   }
00946   else if (TriggerTypeName == "HighMET") {
00947     if (pfmet.pt()<_highPFMETThreshold) return;
00948   }
00949   else if (TriggerTypeName == "Ele") {
00950     if (!selectWElectronEvent(iEvent)) return;
00951   }
00952   else if (TriggerTypeName == "Muon") {
00953     if (!selectWMuonEvent(iEvent)) return;
00954   }
00955 
00956   if (TriggerTypeName != "") DirName = DirName + "/" + TriggerTypeName;
00957 
00958 
00959   // Reconstructed PFMET information
00960   //----------------------------------------------------------------------------
00961   double pfSumET  = pfmet.sumEt();
00962   double pfMETSig = pfmet.mEtSig();
00963   double pfMET    = pfmet.pt();
00964   double pfMEx    = pfmet.px();
00965   double pfMEy    = pfmet.py();
00966   double pfMETPhi = pfmet.phi();
00967 
00968 
00969   // PFMET getters
00970   //----------------------------------------------------------------------------
00971   double pfPhotonEtFraction        = pfmet.photonEtFraction();
00972   double pfPhotonEt                = pfmet.photonEt();
00973   double pfNeutralHadronEtFraction = pfmet.neutralHadronEtFraction();
00974   double pfNeutralHadronEt         = pfmet.neutralHadronEt();
00975   double pfElectronEtFraction      = pfmet.electronEtFraction();
00976   double pfElectronEt              = pfmet.electronEt();
00977   double pfChargedHadronEtFraction = pfmet.chargedHadronEtFraction();
00978   double pfChargedHadronEt         = pfmet.chargedHadronEt();
00979   double pfMuonEtFraction          = pfmet.muonEtFraction();
00980   double pfMuonEt                  = pfmet.muonEt();
00981   double pfHFHadronEtFraction      = pfmet.HFHadronEtFraction();
00982   double pfHFHadronEt              = pfmet.HFHadronEt();
00983   double pfHFEMEtFraction          = pfmet.HFEMEtFraction();
00984   double pfHFEMEt                  = pfmet.HFEMEt();
00985 
00986 
00987   if (pfSumET > _etThreshold) {
00988     
00989     mePfMEx        = _dbe->get(DirName + "/METTask_PfMEx");
00990     mePfMEy        = _dbe->get(DirName + "/METTask_PfMEy");
00991     mePfMET        = _dbe->get(DirName + "/METTask_PfMET");
00992     mePfMETPhi     = _dbe->get(DirName + "/METTask_PfMETPhi");
00993     mePfSumET      = _dbe->get(DirName + "/METTask_PfSumET");
00994     mePfMETSig     = _dbe->get(DirName + "/METTask_PfMETSig");
00995     mePfMET_logx   = _dbe->get(DirName + "/METTask_PfMET_logx");
00996     mePfSumET_logx = _dbe->get(DirName + "/METTask_PfSumET_logx");
00997 
00998     if (mePfMEx        && mePfMEx       ->getRootObject()) mePfMEx       ->Fill(pfMEx);
00999     if (mePfMEy        && mePfMEy       ->getRootObject()) mePfMEy       ->Fill(pfMEy);
01000     if (mePfMET        && mePfMET       ->getRootObject()) mePfMET       ->Fill(pfMET);
01001     if (mePfMETPhi     && mePfMETPhi    ->getRootObject()) mePfMETPhi    ->Fill(pfMETPhi);
01002     if (mePfSumET      && mePfSumET     ->getRootObject()) mePfSumET     ->Fill(pfSumET);
01003     if (mePfMETSig     && mePfMETSig    ->getRootObject()) mePfMETSig    ->Fill(pfMETSig);
01004     if (mePfMET_logx   && mePfMET_logx  ->getRootObject()) mePfMET_logx  ->Fill(log10(pfMET));
01005     if (mePfSumET_logx && mePfSumET_logx->getRootObject()) mePfSumET_logx->Fill(log10(pfSumET));
01006 
01007 
01008     mePhotonEtFraction        = _dbe->get(DirName + "/METTask_PfPhotonEtFraction");
01009     mePhotonEt                = _dbe->get(DirName + "/METTask_PfPhotonEt");
01010     meNeutralHadronEtFraction = _dbe->get(DirName + "/METTask_PfNeutralHadronEtFraction");
01011     meNeutralHadronEt         = _dbe->get(DirName + "/METTask_PfNeutralHadronEt");
01012     meElectronEtFraction      = _dbe->get(DirName + "/METTask_PfElectronEtFraction");
01013     meElectronEt              = _dbe->get(DirName + "/METTask_PfElectronEt");
01014     meChargedHadronEtFraction = _dbe->get(DirName + "/METTask_PfChargedHadronEtFraction");
01015     meChargedHadronEt         = _dbe->get(DirName + "/METTask_PfChargedHadronEt");
01016     meMuonEtFraction          = _dbe->get(DirName + "/METTask_PfMuonEtFraction");
01017     meMuonEt                  = _dbe->get(DirName + "/METTask_PfMuonEt");
01018     meHFHadronEtFraction      = _dbe->get(DirName + "/METTask_PfHFHadronEtFraction");
01019     meHFHadronEt              = _dbe->get(DirName + "/METTask_PfHFHadronEt");
01020     meHFEMEtFraction          = _dbe->get(DirName + "/METTask_PfHFEMEtFraction");
01021     meHFEMEt                  = _dbe->get(DirName + "/METTask_PfHFEMEt");
01022 
01023     if (mePhotonEtFraction        && mePhotonEtFraction       ->getRootObject()) mePhotonEtFraction       ->Fill(pfPhotonEtFraction);
01024     if (mePhotonEt                && mePhotonEt               ->getRootObject()) mePhotonEt               ->Fill(pfPhotonEt);
01025     if (meNeutralHadronEtFraction && meNeutralHadronEtFraction->getRootObject()) meNeutralHadronEtFraction->Fill(pfNeutralHadronEtFraction);
01026     if (meNeutralHadronEt         && meNeutralHadronEt        ->getRootObject()) meNeutralHadronEt        ->Fill(pfNeutralHadronEt);
01027     if (meElectronEtFraction      && meElectronEtFraction     ->getRootObject()) meElectronEtFraction     ->Fill(pfElectronEtFraction);
01028     if (meElectronEt              && meElectronEt             ->getRootObject()) meElectronEt             ->Fill(pfElectronEt);   
01029     if (meChargedHadronEtFraction && meChargedHadronEtFraction->getRootObject()) meChargedHadronEtFraction->Fill(pfChargedHadronEtFraction);
01030     if (meChargedHadronEt         && meChargedHadronEt        ->getRootObject()) meChargedHadronEt        ->Fill(pfChargedHadronEt);
01031     if (meMuonEtFraction          && meMuonEtFraction         ->getRootObject()) meMuonEtFraction         ->Fill(pfMuonEtFraction);      
01032     if (meMuonEt                  && meMuonEt                 ->getRootObject()) meMuonEt                 ->Fill(pfMuonEt);       
01033     if (meHFHadronEtFraction      && meHFHadronEtFraction     ->getRootObject()) meHFHadronEtFraction     ->Fill(pfHFHadronEtFraction);
01034     if (meHFHadronEt              && meHFHadronEt             ->getRootObject()) meHFHadronEt             ->Fill(pfHFHadronEt);   
01035     if (meHFEMEtFraction          && meHFEMEtFraction         ->getRootObject()) meHFEMEtFraction         ->Fill(pfHFEMEtFraction);
01036     if (meHFEMEt                  && meHFEMEt                 ->getRootObject()) meHFEMEt                 ->Fill(pfHFEMEt);       
01037 
01038 
01039     if (_allhist) {
01040       if (bLumiSecPlot) {
01041 
01042         mePfMExLS = _dbe->get(DirName + "/METTask_PfMExLS");
01043         mePfMEyLS = _dbe->get(DirName + "/METTask_PfMEyLS");
01044 
01045         if (mePfMExLS && mePfMExLS->getRootObject()) mePfMExLS->Fill(pfMEx, iEvent.luminosityBlock());
01046         if (mePfMEyLS && mePfMEyLS->getRootObject()) mePfMEyLS->Fill(pfMEy, iEvent.luminosityBlock());
01047       }
01048     }
01049 
01050 
01051     // Fill NPV profiles
01052     //--------------------------------------------------------------------------
01053     mePfMEx_profile   = _dbe->get(DirName + "/METTask_PfMEx_profile");
01054     mePfMEy_profile   = _dbe->get(DirName + "/METTask_PfMEy_profile");
01055     mePfMET_profile   = _dbe->get(DirName + "/METTask_PfMET_profile");
01056     mePfSumET_profile = _dbe->get(DirName + "/METTask_PfSumET_profile");
01057 
01058     if (mePfMEx_profile   && mePfMEx_profile  ->getRootObject()) mePfMEx_profile  ->Fill(_numPV, pfMEx);
01059     if (mePfMEy_profile   && mePfMEy_profile  ->getRootObject()) mePfMEy_profile  ->Fill(_numPV, pfMEy);
01060     if (mePfMET_profile   && mePfMET_profile  ->getRootObject()) mePfMET_profile  ->Fill(_numPV, pfMET);
01061     if (mePfSumET_profile && mePfSumET_profile->getRootObject()) mePfSumET_profile->Fill(_numPV, pfSumET);
01062 
01063 
01064     mePhotonEtFraction_profile        = _dbe->get(DirName + "/METTask_PfPhotonEtFraction_profile");
01065     mePhotonEt_profile                = _dbe->get(DirName + "/METTask_PfPhotonEt_profile");
01066     meNeutralHadronEtFraction_profile = _dbe->get(DirName + "/METTask_PfNeutralHadronEtFraction_profile");
01067     meNeutralHadronEt_profile         = _dbe->get(DirName + "/METTask_PfNeutralHadronEt_profile");
01068     meElectronEtFraction_profile      = _dbe->get(DirName + "/METTask_PfElectronEtFraction_profile");
01069     meElectronEt_profile              = _dbe->get(DirName + "/METTask_PfElectronEt_profile");
01070     meChargedHadronEtFraction_profile = _dbe->get(DirName + "/METTask_PfChargedHadronEtFraction_profile");
01071     meChargedHadronEt_profile         = _dbe->get(DirName + "/METTask_PfChargedHadronEt_profile");
01072     meMuonEtFraction_profile          = _dbe->get(DirName + "/METTask_PfMuonEtFraction_profile");
01073     meMuonEt_profile                  = _dbe->get(DirName + "/METTask_PfMuonEt_profile");
01074     meHFHadronEtFraction_profile      = _dbe->get(DirName + "/METTask_PfHFHadronEtFraction_profile");
01075     meHFHadronEt_profile              = _dbe->get(DirName + "/METTask_PfHFHadronEt_profile");
01076     meHFEMEtFraction_profile          = _dbe->get(DirName + "/METTask_PfHFEMEtFraction_profile");
01077     meHFEMEt_profile                  = _dbe->get(DirName + "/METTask_PfHFEMEt_profile");
01078 
01079     if (mePhotonEtFraction_profile        && mePhotonEtFraction_profile       ->getRootObject()) mePhotonEtFraction_profile       ->Fill(_numPV, pfPhotonEtFraction);
01080     if (mePhotonEt_profile                && mePhotonEt_profile               ->getRootObject()) mePhotonEt_profile               ->Fill(_numPV, pfPhotonEt);
01081     if (meNeutralHadronEtFraction_profile && meNeutralHadronEtFraction_profile->getRootObject()) meNeutralHadronEtFraction_profile->Fill(_numPV, pfNeutralHadronEtFraction);
01082     if (meNeutralHadronEt_profile         && meNeutralHadronEt_profile        ->getRootObject()) meNeutralHadronEt_profile        ->Fill(_numPV, pfNeutralHadronEt);
01083     if (meElectronEtFraction_profile      && meElectronEtFraction_profile     ->getRootObject()) meElectronEtFraction_profile     ->Fill(_numPV, pfElectronEtFraction);
01084     if (meElectronEt_profile              && meElectronEt_profile             ->getRootObject()) meElectronEt_profile             ->Fill(_numPV, pfElectronEt);   
01085     if (meChargedHadronEtFraction_profile && meChargedHadronEtFraction_profile->getRootObject()) meChargedHadronEtFraction_profile->Fill(_numPV, pfChargedHadronEtFraction);
01086     if (meChargedHadronEt_profile         && meChargedHadronEt_profile        ->getRootObject()) meChargedHadronEt_profile        ->Fill(_numPV, pfChargedHadronEt);
01087     if (meMuonEtFraction_profile          && meMuonEtFraction_profile         ->getRootObject()) meMuonEtFraction_profile         ->Fill(_numPV, pfMuonEtFraction);      
01088     if (meMuonEt_profile                  && meMuonEt_profile                 ->getRootObject()) meMuonEt_profile                 ->Fill(_numPV, pfMuonEt);       
01089     if (meHFHadronEtFraction_profile      && meHFHadronEtFraction_profile     ->getRootObject()) meHFHadronEtFraction_profile     ->Fill(_numPV, pfHFHadronEtFraction);
01090     if (meHFHadronEt_profile              && meHFHadronEt_profile             ->getRootObject()) meHFHadronEt_profile             ->Fill(_numPV, pfHFHadronEt);   
01091     if (meHFEMEtFraction_profile          && meHFEMEtFraction_profile         ->getRootObject()) meHFEMEtFraction_profile         ->Fill(_numPV, pfHFEMEtFraction);
01092     if (meHFEMEt_profile                  && meHFEMEt_profile                 ->getRootObject()) meHFEMEt_profile                 ->Fill(_numPV, pfHFEMEt);       
01093   }
01094 }
01095 
01096 
01097 //------------------------------------------------------------------------------
01098 // selectHighPtJetEvent
01099 //------------------------------------------------------------------------------
01100 bool PFMETAnalyzer::selectHighPtJetEvent(const edm::Event& iEvent){
01101 
01102   bool return_value=false;
01103 
01104   edm::Handle<reco::PFJetCollection> pfJets;
01105   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
01106   if (!pfJets.isValid()) {
01107     LogDebug("") << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
01108     if (_verbose) std::cout << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
01109   }
01110 
01111   for (reco::PFJetCollection::const_iterator pf = pfJets->begin(); 
01112        pf!=pfJets->end(); ++pf){
01113     if (pf->pt()>_highPtPFJetThreshold){
01114       return_value=true;
01115     }
01116   }
01117   
01118   return return_value;
01119 }
01120 
01121 // // ***********************************************************
01122 bool PFMETAnalyzer::selectLowPtJetEvent(const edm::Event& iEvent){
01123 
01124   bool return_value=false;
01125 
01126   edm::Handle<reco::PFJetCollection> pfJets;
01127   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
01128   if (!pfJets.isValid()) {
01129     LogDebug("") << "PFMETAnalyzer: Could not find jet product" << std::endl;
01130     if (_verbose) std::cout << "PFMETAnalyzer: Could not find jet product" << std::endl;
01131   }
01132 
01133   for (reco::PFJetCollection::const_iterator cal = pfJets->begin(); 
01134        cal!=pfJets->end(); ++cal){
01135     if (cal->pt()>_lowPtPFJetThreshold){
01136       return_value=true;
01137     }
01138   }
01139 
01140   return return_value;
01141 
01142 }
01143 
01144 // ***********************************************************
01145 bool PFMETAnalyzer::selectWElectronEvent(const edm::Event& iEvent){
01146 
01147   bool return_value=true;
01148 
01149   /*
01150     W-electron event selection comes here
01151   */
01152 
01153   return return_value;
01154 
01155 }
01156 
01157 // ***********************************************************
01158 bool PFMETAnalyzer::selectWMuonEvent(const edm::Event& iEvent){
01159 
01160   bool return_value=true;
01161 
01162   /*
01163     W-muon event selection comes here
01164   */
01165 
01166   return return_value;
01167 
01168 }
01169
/data/refman/pasoursint/CMSSW_5_3_4/src/DQMOffline/JetMET/src/PFMETAnalyzer.cc
