Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2012/04/12 15:42:58 $
00005  *  $Revision: 1.44 $
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   _tightHcalFiltering   = theCleaningParameters.getParameter<bool>("tightHcalFiltering");
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   HcalNoiseSummaryTag         = parameters.getParameter<edm::InputTag>("HcalNoiseSummary");
00132   BeamHaloSummaryTag          = parameters.getParameter<edm::InputTag>("BeamHaloSummaryLabel");
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"); // High Pt Jet threshold
00142   _lowPtPFJetThreshold  = parameters.getParameter<double>("LowPtJetThreshold");   // Low Pt Jet threshold
00143   _highPFMETThreshold   = parameters.getParameter<double>("HighMETThreshold");     // High MET threshold
00144   _lowPFMETThreshold    = parameters.getParameter<double>("LowMETThreshold");       // Low MET threshold
00145 
00146   //
00147   jetID = new reco::helper::JetIDHelper(parameters.getParameter<ParameterSet>("JetIDParams"));
00148 
00149   // DQStore stuff
00150   LogTrace(metname)<<"[PFMETAnalyzer] Parameters initialization";
00151   std::string DirName = "JetMET/MET/"+_source;
00152   dbe->setCurrentFolder(DirName);
00153 
00154   metME = dbe->book1D("metReco", "metReco", 4, 1, 5);
00155   metME->setBinLabel(3,"PFMET",1);
00156 
00157   _dbe = dbe;
00158 
00159   _FolderNames.push_back("All");
00160   _FolderNames.push_back("BasicCleanup");
00161   _FolderNames.push_back("ExtraCleanup");
00162   _FolderNames.push_back("HcalNoiseFilter");
00163   _FolderNames.push_back("HcalNoiseFilterTight");
00164   _FolderNames.push_back("JetIDMinimal");
00165   _FolderNames.push_back("JetIDLoose");
00166   _FolderNames.push_back("JetIDTight");
00167   _FolderNames.push_back("BeamHaloIDTightPass");
00168   _FolderNames.push_back("BeamHaloIDLoosePass");
00169   _FolderNames.push_back("Triggers");
00170   _FolderNames.push_back("PV");
00171 
00172   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); 
00173        ic != _FolderNames.end(); ic++){
00174     if (*ic=="All")             bookMESet(DirName+"/"+*ic);
00175     if (_cleanupSelection){
00176     if (*ic=="BasicCleanup")    bookMESet(DirName+"/"+*ic);
00177     if (*ic=="ExtraCleanup")    bookMESet(DirName+"/"+*ic);
00178     }
00179     if (_allSelection){
00180       if (*ic=="HcalNoiseFilter")      bookMESet(DirName+"/"+*ic);
00181       if (*ic=="HcalNoiseFilterTight") bookMESet(DirName+"/"+*ic);
00182       if (*ic=="JetIDMinimal")         bookMESet(DirName+"/"+*ic);
00183       if (*ic=="JetIDLoose")           bookMESet(DirName+"/"+*ic);
00184       if (*ic=="JetIDTight")           bookMESet(DirName+"/"+*ic);
00185       if (*ic=="BeamHaloIDTightPass")  bookMESet(DirName+"/"+*ic);
00186       if (*ic=="BeamHaloIDLoosePass")  bookMESet(DirName+"/"+*ic);
00187       if (*ic=="Triggers")             bookMESet(DirName+"/"+*ic);
00188       if (*ic=="PV")                   bookMESet(DirName+"/"+*ic);
00189     }
00190   }
00191 }
00192 
00193 // ***********************************************************
00194 void PFMETAnalyzer::endJob() {
00195 
00196   delete jetID;
00197   delete DCSFilter;
00198 
00199 }
00200 
00201 // ***********************************************************
00202 void PFMETAnalyzer::bookMESet(std::string DirName)
00203 {
00204 
00205   bool bLumiSecPlot=false;
00206   if (DirName.find("All")!=std::string::npos) bLumiSecPlot=true;
00207 
00208   bookMonitorElement(DirName,bLumiSecPlot);
00209 
00210   if ( _HighPtJetEventFlag->on() ) {
00211     bookMonitorElement(DirName+"/"+"HighPtJet",false);
00212     meTriggerName_HighPtJet = _dbe->bookString("triggerName_HighPtJet", highPtJetExpr_[0]);
00213   }  
00214 
00215   if ( _LowPtJetEventFlag->on() ) {
00216     bookMonitorElement(DirName+"/"+"LowPtJet",false);
00217     meTriggerName_LowPtJet = _dbe->bookString("triggerName_LowPtJet", lowPtJetExpr_[0]);
00218   }
00219 
00220   if ( _MinBiasEventFlag->on() ) {
00221     bookMonitorElement(DirName+"/"+"MinBias",false);
00222     meTriggerName_MinBias = _dbe->bookString("triggerName_MinBias", minbiasExpr_[0]);
00223     if (_verbose) std::cout << "_MinBiasEventFlag is on, folder created\n";
00224   }
00225 
00226   if ( _HighMETEventFlag->on() ) {
00227     bookMonitorElement(DirName+"/"+"HighMET",false);
00228     meTriggerName_HighMET = _dbe->bookString("triggerName_HighMET", highMETExpr_[0]);
00229   }
00230 
00231   if ( _LowMETEventFlag->on() ) {
00232     bookMonitorElement(DirName+"/"+"LowMET",false);
00233     meTriggerName_LowMET = _dbe->bookString("triggerName_LowMET", lowMETExpr_[0]);
00234   }
00235 
00236   if ( _EleEventFlag->on() ) {
00237     bookMonitorElement(DirName+"/"+"Ele",false);
00238     meTriggerName_Ele = _dbe->bookString("triggerName_Ele", elecExpr_[0]);
00239     if (_verbose) std::cout << "_EleEventFlag is on, folder created\n";
00240   }
00241 
00242   if ( _MuonEventFlag->on() ) {
00243     bookMonitorElement(DirName+"/"+"Muon",false);
00244     meTriggerName_Muon = _dbe->bookString("triggerName_Muon", muonExpr_[0]);
00245     if (_verbose) std::cout << "_MuonEventFlag is on, folder created\n";
00246   }
00247 }
00248 
00249 // ***********************************************************
00250 void PFMETAnalyzer::bookMonitorElement(std::string DirName, bool bLumiSecPlot=false)
00251 {
00252 
00253   if (_verbose) std::cout << "booMonitorElement " << DirName << std::endl;
00254   _dbe->setCurrentFolder(DirName);
00255  
00256   //meNevents              = _dbe->book1D("METTask_Nevents", "METTask_Nevents"   ,1,0,1);
00257   mePfMEx                = _dbe->book1D("METTask_PfMEx",   "METTask_PfMEx"   ,200,-500,500); 
00258   mePfMEx->setAxisTitle("MEx [GeV]",1);
00259   mePfMEy                = _dbe->book1D("METTask_PfMEy",   "METTask_PfMEy"   ,200,-500,500); 
00260   mePfMEy->setAxisTitle("MEy [GeV]",1);
00261   //mePfEz                 = _dbe->book1D("METTask_PfEz",    "METTask_PfEz"    ,500,-500,500);
00262   //mePfEz->setAxisTitle("MEz [GeV]",1);
00263   mePfMETSig             = _dbe->book1D("METTask_PfMETSig","METTask_PfMETSig",51,0,51);
00264   mePfMETSig->setAxisTitle("METSig",1);
00265   mePfMET                = _dbe->book1D("METTask_PfMET",   "METTask_PfMET"   ,200,0,1000); 
00266   mePfMET->setAxisTitle("MET [GeV]",1);
00267   mePfMETPhi             = _dbe->book1D("METTask_PfMETPhi","METTask_PfMETPhi",60,-TMath::Pi(),TMath::Pi());
00268   mePfMETPhi->setAxisTitle("METPhi [rad]",1);
00269   mePfSumET              = _dbe->book1D("METTask_PfSumET", "METTask_PfSumET" ,400,0,2000); 
00270   mePfSumET->setAxisTitle("SumET [GeV]",1);
00271 
00272   mePfMET_logx           = _dbe->book1D("METTask_PfMET_logx",   "METTask_PfMET_logx"   ,40,-1.,7.);
00273   mePfMET_logx->setAxisTitle("log(MET) [GeV]",1);
00274   mePfSumET_logx         = _dbe->book1D("METTask_PfSumET_logx", "METTask_PfSumET_logx" ,40,-1.,7.);
00275   mePfSumET_logx->setAxisTitle("log(SumET) [GeV]",1);
00276 
00277   mePfNeutralEMFraction  = _dbe->book1D("METTask_PfNeutralEMFraction", "METTask_PfNeutralEMFraction" ,50,0.,1.);
00278   mePfNeutralEMFraction->setAxisTitle("Pf Neutral EM Fraction",1);
00279   mePfNeutralHadFraction = _dbe->book1D("METTask_PfNeutralHadFraction","METTask_PfNeutralHadFraction",50,0.,1.);
00280   mePfNeutralHadFraction->setAxisTitle("Pf Neutral Had Fraction",1);
00281   mePfChargedEMFraction  = _dbe->book1D("METTask_PfChargedEMFraction", "METTask_PfChargedEMFraction" ,50,0.,1.);
00282   mePfChargedEMFraction->setAxisTitle("Pf Charged EM Fraction",1);
00283   mePfChargedHadFraction = _dbe->book1D("METTask_PfChargedHadFraction","METTask_PfChargedHadFraction",50,0.,1.);
00284   mePfChargedHadFraction->setAxisTitle("Pf Charged Had Fraction",1);
00285   mePfMuonFraction       = _dbe->book1D("METTask_PfMuonFraction",      "METTask_PfMuonFraction"      ,50,0.,1.);
00286   mePfMuonFraction->setAxisTitle("Pf Muon Fraction",1);
00287 
00288   //mePfMETIonFeedbck      = _dbe->book1D("METTask_PfMETIonFeedbck", "METTask_PfMETIonFeedbck" ,500,0,1000);
00289   //mePfMETIonFeedbck->setAxisTitle("MET [GeV]",1);
00290   //mePfMETHPDNoise        = _dbe->book1D("METTask_PfMETHPDNoise",   "METTask_PfMETHPDNoise"   ,500,0,1000);
00291   //mePfMETHPDNoise->setAxisTitle("MET [GeV]",1);
00292   //mePfMETRBXNoise        = _dbe->book1D("METTask_PfMETRBXNoise",   "METTask_PfMETRBXNoise"   ,500,0,1000);
00293   //mePfMETRBXNoise->setAxisTitle("MET [GeV]",1);
00294 
00295   if (_allhist){
00296     if (bLumiSecPlot){
00297       mePfMExLS              = _dbe->book2D("METTask_PfMEx_LS","METTask_PfMEx_LS",200,-200,200,50,0.,500.);
00298       mePfMExLS->setAxisTitle("MEx [GeV]",1);
00299       mePfMExLS->setAxisTitle("Lumi Section",2);
00300       mePfMEyLS              = _dbe->book2D("METTask_PfMEy_LS","METTask_PfMEy_LS",200,-200,200,50,0.,500.);
00301       mePfMEyLS->setAxisTitle("MEy [GeV]",1);
00302       mePfMEyLS->setAxisTitle("Lumi Section",2);
00303     }
00304   }
00305 
00306 
00307   // Book NPV profiles
00308   //----------------------------------------------------------------------------
00309   mePfMEx_profile   = _dbe->bookProfile("METTask_PfMEx_profile",   "MEx [GeV]",   nbinsPV, PVlow, PVup, 100, -500,  500);
00310   mePfMEy_profile   = _dbe->bookProfile("METTask_PfMEy_profile",   "MEy [GeV]",   nbinsPV, PVlow, PVup, 100, -500,  500); 
00311   mePfMET_profile   = _dbe->bookProfile("METTask_PfMET_profile",   "MET [GeV]",   nbinsPV, PVlow, PVup, 100,    0, 1000); 
00312   mePfSumET_profile = _dbe->bookProfile("METTask_PfSumET_profile", "SumET [GeV]", nbinsPV, PVlow, PVup, 200,    0, 2000); 
00313 
00314   mePfNeutralEMFraction_profile  = _dbe->bookProfile("METTask_PfNeutralEMFraction_profile",  "PF neutral EM fraction",  nbinsPV, PVlow, PVup, 50, 0, 1);
00315   mePfNeutralHadFraction_profile = _dbe->bookProfile("METTask_PfNeutralHadFraction_profile", "PF neutral HAD fraction", nbinsPV, PVlow, PVup, 50, 0, 1);
00316   mePfChargedEMFraction_profile  = _dbe->bookProfile("METTask_PfChargedEMFraction_profile",  "PF charged EM fraction",  nbinsPV, PVlow, PVup, 50, 0, 1);
00317   mePfChargedHadFraction_profile = _dbe->bookProfile("METTask_PfChargedHadFraction_profile", "PF charged HAD fraction", nbinsPV, PVlow, PVup, 50, 0, 1);
00318   mePfMuonFraction_profile       = _dbe->bookProfile("METTask_PfMuonFraction_profile",       "PF muon fraction",        nbinsPV, PVlow, PVup, 50, 0, 1);
00319 
00320 
00321   // Set NPV profiles x-axis title
00322   //----------------------------------------------------------------------------
00323   mePfMEx_profile  ->setAxisTitle("nvtx",1);
00324   mePfMEy_profile  ->setAxisTitle("nvtx",1);
00325   mePfMET_profile  ->setAxisTitle("nvtx",1);
00326   mePfSumET_profile->setAxisTitle("nvtx",1);
00327 
00328   mePfNeutralEMFraction_profile ->setAxisTitle("nvtx",1);
00329   mePfNeutralHadFraction_profile->setAxisTitle("nvtx",1);
00330   mePfChargedEMFraction_profile ->setAxisTitle("nvtx",1);
00331   mePfChargedHadFraction_profile->setAxisTitle("nvtx",1);
00332   mePfMuonFraction_profile      ->setAxisTitle("nvtx",1);
00333 }
00334 
00335 
00336 // ***********************************************************
00337 void PFMETAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00338 {
00339   if ( _HighPtJetEventFlag->on() ) _HighPtJetEventFlag->initRun( iRun, iSetup );
00340   if ( _LowPtJetEventFlag ->on() ) _LowPtJetEventFlag ->initRun( iRun, iSetup );
00341   if ( _MinBiasEventFlag  ->on() ) _MinBiasEventFlag  ->initRun( iRun, iSetup );
00342   if ( _HighMETEventFlag  ->on() ) _HighMETEventFlag  ->initRun( iRun, iSetup );
00343   if ( _LowMETEventFlag   ->on() ) _LowMETEventFlag   ->initRun( iRun, iSetup );
00344   if ( _EleEventFlag      ->on() ) _EleEventFlag      ->initRun( iRun, iSetup );
00345   if ( _MuonEventFlag     ->on() ) _MuonEventFlag     ->initRun( iRun, iSetup );
00346 
00347   if (_HighPtJetEventFlag->on() && _HighPtJetEventFlag->expressionsFromDB(_HighPtJetEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00348     highPtJetExpr_ = _HighPtJetEventFlag->expressionsFromDB(_HighPtJetEventFlag->hltDBKey(), iSetup);
00349   if (_LowPtJetEventFlag->on() && _LowPtJetEventFlag->expressionsFromDB(_LowPtJetEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00350     lowPtJetExpr_  = _LowPtJetEventFlag->expressionsFromDB(_LowPtJetEventFlag->hltDBKey(),   iSetup);
00351   if (_HighMETEventFlag->on() && _HighMETEventFlag->expressionsFromDB(_HighMETEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00352     highMETExpr_   = _HighMETEventFlag->expressionsFromDB(_HighMETEventFlag->hltDBKey(),     iSetup);
00353   if (_LowMETEventFlag->on() && _LowMETEventFlag->expressionsFromDB(_LowMETEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00354     lowMETExpr_    = _LowMETEventFlag->expressionsFromDB(_LowMETEventFlag->hltDBKey(),       iSetup);
00355   if (_MuonEventFlag->on() && _MuonEventFlag->expressionsFromDB(_MuonEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00356     muonExpr_      = _MuonEventFlag->expressionsFromDB(_MuonEventFlag->hltDBKey(),           iSetup);
00357   if (_EleEventFlag->on() && _EleEventFlag->expressionsFromDB(_EleEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00358     elecExpr_      = _EleEventFlag->expressionsFromDB(_EleEventFlag->hltDBKey(),             iSetup);
00359   if (_MinBiasEventFlag->on() && _MinBiasEventFlag->expressionsFromDB(_MinBiasEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00360     minbiasExpr_   = _MinBiasEventFlag->expressionsFromDB(_MinBiasEventFlag->hltDBKey(),     iSetup);
00361 
00362 }
00363 
00364 // ***********************************************************
00365 void PFMETAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup, DQMStore * dbe)
00366 {
00367   
00368   //
00369   //--- Check the time length of the Run from the lumi section plots
00370 
00371   std::string dirName = "JetMET/MET/"+_source+"/";
00372   _dbe->setCurrentFolder(dirName);
00373 
00374   TH1F* tlumisec;
00375 
00376   MonitorElement *meLumiSec = _dbe->get("aaa");
00377   meLumiSec = _dbe->get("JetMET/lumisec");
00378 
00379   int totlsec=0;
00380   double totltime=0.;
00381   if ( meLumiSec->getRootObject() ) {
00382     tlumisec = meLumiSec->getTH1F();
00383     for (int i=0; i<500; i++){
00384       if (tlumisec->GetBinContent(i+1)) totlsec++;
00385     }
00386     totltime = double(totlsec*90); // one lumi sec ~ 90 (sec)
00387   }
00388 
00389   if (totltime==0.) totltime=1.; 
00390 
00391   //
00392   //--- Make the integrated plots with rate (Hz)
00393 
00394   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); ic != _FolderNames.end(); ic++)
00395     {
00396 
00397       std::string DirName;
00398       DirName = dirName+*ic;
00399 
00400       makeRatePlot(DirName,totltime);
00401       if ( _HighPtJetEventFlag->on() ) 
00402         makeRatePlot(DirName+"/"+"triggerName_HighJetPt",totltime);
00403       if ( _LowPtJetEventFlag->on() ) 
00404         makeRatePlot(DirName+"/"+"triggerName_LowJetPt",totltime);
00405       if ( _MinBiasEventFlag->on() ) 
00406         makeRatePlot(DirName+"/"+"triggerName_MinBias",totltime);
00407       if ( _HighMETEventFlag->on() ) 
00408         makeRatePlot(DirName+"/"+"triggerName_HighMET",totltime);
00409       if ( _LowMETEventFlag->on() ) 
00410         makeRatePlot(DirName+"/"+"triggerName_LowMET",totltime);
00411       if ( _EleEventFlag->on() ) 
00412         makeRatePlot(DirName+"/"+"triggerName_Ele",totltime);
00413       if ( _MuonEventFlag->on() ) 
00414         makeRatePlot(DirName+"/"+"triggerName_Muon",totltime);
00415     }
00416 }
00417 
00418 // ***********************************************************
00419 void PFMETAnalyzer::makeRatePlot(std::string DirName, double totltime)
00420 {
00421 
00422   _dbe->setCurrentFolder(DirName);
00423   MonitorElement *mePfMET = _dbe->get(DirName+"/"+"METTask_PfMET");
00424 
00425   TH1F* tPfMET;
00426   TH1F* tPfMETRate;
00427 
00428   if ( mePfMET )
00429     if ( mePfMET->getRootObject() ) {
00430       tPfMET     = mePfMET->getTH1F();
00431       
00432       // Integral plot & convert number of events to rate (hz)
00433       tPfMETRate = (TH1F*) tPfMET->Clone("METTask_PfMETRate");
00434       for (int i = tPfMETRate->GetNbinsX()-1; i>=0; i--){
00435         tPfMETRate->SetBinContent(i+1,tPfMETRate->GetBinContent(i+2)+tPfMET->GetBinContent(i+1));
00436       }
00437       for (int i = 0; i<tPfMETRate->GetNbinsX(); i++){
00438         tPfMETRate->SetBinContent(i+1,tPfMETRate->GetBinContent(i+1)/double(totltime));
00439       }
00440 
00441       tPfMETRate->SetName("METTask_PfMETRate");
00442       tPfMETRate->SetTitle("METTask_PfMETRate");
00443       mePfMETRate      = _dbe->book1D("METTask_PfMETRate",tPfMETRate);
00444       
00445     }
00446 }
00447 
00448 // ***********************************************************
00449 void PFMETAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
00450                             const edm::TriggerResults& triggerResults) {
00451 
00452   if (_verbose) std::cout << "PfMETAnalyzer analyze" << std::endl;
00453 
00454   LogTrace(metname)<<"[PFMETAnalyzer] Analyze PFMET";
00455 
00456   metME->Fill(3);
00457 
00458   // ==========================================================  
00459   // Trigger information 
00460   //
00461   _trig_JetMB=0;
00462   _trig_HighPtJet=0;
00463   _trig_LowPtJet=0;
00464   _trig_MinBias=0;
00465   _trig_HighMET=0;
00466   _trig_LowMET=0;
00467   _trig_Ele=0;
00468   _trig_Muon=0;
00469   _trig_PhysDec=0;
00470   if(&triggerResults) {   
00471 
00473 
00474     //
00475     //
00476     // Check how many HLT triggers are in triggerResults 
00477     int ntrigs = triggerResults.size();
00478     if (_verbose) std::cout << "ntrigs=" << ntrigs << std::endl;
00479     
00480     //
00481     //
00482     // If index=ntrigs, this HLT trigger doesn't exist in the HLT table for this data.
00483     const edm::TriggerNames & triggerNames = iEvent.triggerNames(triggerResults);
00484     
00485     //
00486     //
00487     const unsigned int nTrig(triggerNames.size());
00488     for (unsigned int i=0;i<nTrig;++i)
00489       {
00490         if (triggerNames.triggerName(i).find(highPtJetExpr_[0].substr(0,highPtJetExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00491           _trig_HighPtJet=true;
00492         else if (triggerNames.triggerName(i).find(lowPtJetExpr_[0].substr(0,lowPtJetExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00493           _trig_LowPtJet=true;
00494         else if (triggerNames.triggerName(i).find(highMETExpr_[0].substr(0,highMETExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00495           _trig_HighMET=true;
00496         else if (triggerNames.triggerName(i).find(lowMETExpr_[0].substr(0,lowMETExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00497           _trig_LowMET=true;
00498         else if (triggerNames.triggerName(i).find(muonExpr_[0].substr(0,muonExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00499           _trig_Muon=true;
00500         else if (triggerNames.triggerName(i).find(elecExpr_[0].substr(0,elecExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00501           _trig_Ele=true;
00502         else if (triggerNames.triggerName(i).find(minbiasExpr_[0].substr(0,minbiasExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00503           _trig_MinBias=true;
00504       }
00505 
00506       // count number of requested Jet or MB HLT paths which have fired
00507     for (unsigned int i=0; i!=HLTPathsJetMBByName_.size(); i++) {
00508       unsigned int triggerIndex = triggerNames.triggerIndex(HLTPathsJetMBByName_[i]);
00509       if (triggerIndex<triggerResults.size()) {
00510         if (triggerResults.accept(triggerIndex)) {
00511           _trig_JetMB++;
00512         }
00513       }
00514     }
00515     // for empty input vectors (n==0), take all HLT triggers!
00516     if (HLTPathsJetMBByName_.size()==0) _trig_JetMB=triggerResults.size()-1;
00517 
00518     //
00519     /*
00520       if ( _HighPtJetEventFlag->on() && _HighPtJetEventFlag->accept( iEvent, iSetup) )
00521       _trig_HighPtJet=1;
00522       
00523       if ( _LowPtJetEventFlag->on() && _LowPtJetEventFlag->accept( iEvent, iSetup) )
00524       _trig_LowPtJet=1;
00525       
00526       if ( _MinBiasEventFlag->on() && _MinBiasEventFlag->accept( iEvent, iSetup) )
00527       _trig_MinBias=1;
00528       
00529       if ( _HighMETEventFlag->on() && _HighMETEventFlag->accept( iEvent, iSetup) )
00530       _trig_HighMET=1;
00531       
00532       if ( _LowMETEventFlag->on() && _LowMETEventFlag->accept( iEvent, iSetup) )
00533       _trig_LowMET=1;
00534       
00535       if ( _EleEventFlag->on() && _EleEventFlag->accept( iEvent, iSetup) )
00536       _trig_Ele=1;
00537       
00538       if ( _MuonEventFlag->on() && _MuonEventFlag->accept( iEvent, iSetup) )
00539       _trig_Muon=1;
00540     */
00541     
00542     if (triggerNames.triggerIndex(_hlt_PhysDec)   != triggerNames.size() &&
00543         triggerResults.accept(triggerNames.triggerIndex(_hlt_PhysDec)))   _trig_PhysDec=1;
00544   } else {
00545 
00546     edm::LogInfo("PFMetAnalyzer") << "TriggerResults::HLT not found, "
00547       "automatically select events"; 
00548 
00549     // TriggerResults object not found. Look at all events.    
00550     _trig_JetMB=1;
00551   }
00552 
00553   // ==========================================================
00554   // PfMET information
00555   
00556   // **** Get the MET container  
00557   edm::Handle<reco::PFMETCollection> pfmetcoll;
00558   iEvent.getByLabel(thePfMETCollectionLabel, pfmetcoll);
00559   
00560   if(!pfmetcoll.isValid()) return;
00561 
00562   const PFMETCollection *pfmetcol = pfmetcoll.product();
00563   const PFMET *pfmet;
00564   pfmet = &(pfmetcol->front());
00565     
00566   LogTrace(metname)<<"[PfMETAnalyzer] Call to the PfMET analyzer";
00567 
00568   // ==========================================================
00569   //
00570   edm::Handle<HcalNoiseRBXCollection> HRBXCollection;
00571   iEvent.getByLabel(HcalNoiseRBXCollectionTag,HRBXCollection);
00572   if (!HRBXCollection.isValid()) {
00573     LogDebug("") << "PfMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00574     if (_verbose) std::cout << "PfMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00575   }
00576   
00577   edm::Handle<HcalNoiseSummary> HNoiseSummary;
00578   iEvent.getByLabel(HcalNoiseSummaryTag,HNoiseSummary);
00579   if (!HNoiseSummary.isValid()) {
00580     LogDebug("") << "PfMETAnalyzer: Could not find Hcal NoiseSummary product" << std::endl;
00581     if (_verbose) std::cout << "PfMETAnalyzer: Could not find Hcal NoiseSummary product" << std::endl;
00582   }
00583 
00584   edm::Handle<reco::CaloJetCollection> caloJets;
00585   iEvent.getByLabel(theJetCollectionLabel, caloJets);
00586   if (!caloJets.isValid()) {
00587     LogDebug("") << "PFMETAnalyzer: Could not find jet product" << std::endl;
00588     if (_verbose) std::cout << "PFMETAnalyzer: Could not find jet product" << std::endl;
00589   }
00590 
00591   edm::Handle<edm::View<PFCandidate> > pfCandidates;
00592   iEvent.getByLabel(PFCandidatesTag, pfCandidates);
00593   if (!pfCandidates.isValid()) {
00594     LogDebug("") << "PfMETAnalyzer: Could not find pfcandidates product" << std::endl;
00595     if (_verbose) std::cout << "PfMETAnalyzer: Could not find pfcandidates product" << std::endl;
00596   }
00597 
00598   edm::Handle<reco::PFJetCollection> pfJets;
00599   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
00600   if (!pfJets.isValid()) {
00601     LogDebug("") << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
00602     if (_verbose) std::cout << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
00603   }
00604   // ==========================================================
00605   // PfMET sanity check
00606 
00607   if (_source=="PfMET") validateMET(*pfmet, pfCandidates);
00608   
00609   // ==========================================================
00610   // JetID 
00611 
00612   if (_verbose) std::cout << "JetID starts" << std::endl;
00613   
00614   //
00615   // --- Minimal cuts
00616   //
00617   bool bJetIDMinimal=true;
00618   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00619        cal!=caloJets->end(); ++cal){
00620     jetID->calculate(iEvent, *cal);
00621     if (cal->pt()>10.){
00622       if (fabs(cal->eta())<=2.6 && 
00623           cal->emEnergyFraction()<=0.01) bJetIDMinimal=false;
00624     }
00625   }
00626 
00627   //
00628   // --- Loose cuts, not PF specific for now!
00629   //
00630   bool bJetIDLoose=true;
00631   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00632        cal!=caloJets->end(); ++cal){ 
00633     jetID->calculate(iEvent, *cal);
00634     if (_verbose) std::cout << jetID->n90Hits() << " " 
00635                             << jetID->restrictedEMF() << " "
00636                             << cal->pt() << std::endl;
00637     if (cal->pt()>10.){
00638       //
00639       // for all regions
00640       if (jetID->n90Hits()<2)  bJetIDLoose=false; 
00641       if (jetID->fHPD()>=0.98) bJetIDLoose=false; 
00642       //if (jetID->restrictedEMF()<0.01) bJetIDLoose=false; 
00643       //
00644       // for non-forward
00645       if (fabs(cal->eta())<2.55){
00646         if (cal->emEnergyFraction()<=0.01) bJetIDLoose=false; 
00647       }
00648       // for forward
00649       else {
00650         if (cal->emEnergyFraction()<=-0.9) bJetIDLoose=false; 
00651         if (cal->pt()>80.){
00652           if (cal->emEnergyFraction()>= 1.0) bJetIDLoose=false; 
00653         }
00654       } // forward vs non-forward
00655     }   // pt>10 GeV/c
00656   }     // calor-jets loop
00657  
00658   //
00659   // --- Tight cuts
00660   //
00661   bool bJetIDTight=true;
00662   bJetIDTight=bJetIDLoose;
00663   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00664        cal!=caloJets->end(); ++cal){
00665     jetID->calculate(iEvent, *cal);
00666     if (cal->pt()>25.){
00667       //
00668       // for all regions
00669       if (jetID->fHPD()>=0.95) bJetIDTight=false; 
00670       //
00671       // for 1.0<|eta|<1.75
00672       if (fabs(cal->eta())>=1.00 && fabs(cal->eta())<1.75){
00673         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00674       }
00675       //
00676       // for 1.75<|eta|<2.55
00677       else if (fabs(cal->eta())>=1.75 && fabs(cal->eta())<2.55){
00678         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00679       }
00680       //
00681       // for 2.55<|eta|<3.25
00682       else if (fabs(cal->eta())>=2.55 && fabs(cal->eta())<3.25){
00683         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3) bJetIDTight=false; 
00684         if (cal->pt()>=50. && cal->pt()< 80. && cal->emEnergyFraction()<=-0.2) bJetIDTight=false; 
00685         if (cal->pt()>=80. && cal->pt()<340. && cal->emEnergyFraction()<=-0.1) bJetIDTight=false; 
00686         if (cal->pt()>=340.                  && cal->emEnergyFraction()<=-0.1 
00687             && cal->emEnergyFraction()>=0.95) bJetIDTight=false; 
00688       }
00689       //
00690       // for 3.25<|eta|
00691       else if (fabs(cal->eta())>=3.25){
00692         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3
00693             && cal->emEnergyFraction()>=0.90) bJetIDTight=false; 
00694         if (cal->pt()>=50. && cal->pt()<130. && cal->emEnergyFraction()<=-0.2
00695             && cal->emEnergyFraction()>=0.80) bJetIDTight=false; 
00696         if (cal->pt()>=130.                  && cal->emEnergyFraction()<=-0.1 
00697             && cal->emEnergyFraction()>=0.70) bJetIDTight=false; 
00698       }
00699     }   // pt>10 GeV/c
00700   }     // calor-jets loop
00701   
00702   if (_verbose) std::cout << "JetID ends" << std::endl;
00703 
00704 
00705   // ==========================================================
00706   // HCAL Noise filter
00707   
00708   bool bHcalNoiseFilter      = HNoiseSummary->passLooseNoiseFilter();
00709   bool bHcalNoiseFilterTight = HNoiseSummary->passTightNoiseFilter();
00710 
00711   // ==========================================================
00712   // Get BeamHaloSummary
00713   edm::Handle<BeamHaloSummary> TheBeamHaloSummary ;
00714   iEvent.getByLabel(BeamHaloSummaryTag, TheBeamHaloSummary) ;
00715 
00716   bool bBeamHaloIDTightPass = true;
00717   bool bBeamHaloIDLoosePass = true;
00718 
00719   if(!TheBeamHaloSummary.isValid()) {
00720 
00721     const BeamHaloSummary TheSummary = (*TheBeamHaloSummary.product() );
00722 
00723     if( !TheSummary.EcalLooseHaloId()  && !TheSummary.HcalLooseHaloId() && 
00724         !TheSummary.CSCLooseHaloId()   && !TheSummary.GlobalLooseHaloId() )
00725       bBeamHaloIDLoosePass = false;
00726 
00727     if( !TheSummary.EcalTightHaloId()  && !TheSummary.HcalTightHaloId() && 
00728         !TheSummary.CSCTightHaloId()   && !TheSummary.GlobalTightHaloId() )
00729       bBeamHaloIDTightPass = false;
00730 
00731   }
00732 
00733   // ==========================================================
00734   //Vertex information
00735   
00736   _numPV = 0;
00737   bool bPrimaryVertex = true;
00738   if(_doPVCheck){
00739     bPrimaryVertex = false;
00740     Handle<VertexCollection> vertexHandle;
00741 
00742     iEvent.getByLabel(vertexTag, vertexHandle);
00743 
00744     if (!vertexHandle.isValid()) {
00745       LogDebug("") << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00746       if (_verbose) std::cout << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00747     }
00748     
00749     if ( vertexHandle.isValid() ){
00750       VertexCollection vertexCollection = *(vertexHandle.product());
00751       int vertex_number     = vertexCollection.size();
00752       VertexCollection::const_iterator v = vertexCollection.begin();
00753       for ( ; v != vertexCollection.end(); ++v) {
00754         double vertex_chi2    = v->normalizedChi2();
00755         double vertex_ndof    = v->ndof();
00756         bool   fakeVtx        = v->isFake();
00757         double vertex_Z       = v->z();
00758         
00759         if (  !fakeVtx
00760               && vertex_number>=_nvtx_min
00761               && vertex_ndof   >_vtxndof_min
00762               && vertex_chi2   <_vtxchi2_max
00763               && fabs(vertex_Z)<_vtxz_max ) {
00764           bPrimaryVertex = true;
00765           ++_numPV;
00766         }
00767       }
00768     }
00769   }
00770   // ==========================================================
00771 
00772   edm::Handle< L1GlobalTriggerReadoutRecord > gtReadoutRecord;
00773   iEvent.getByLabel( gtTag, gtReadoutRecord);
00774 
00775   if (!gtReadoutRecord.isValid()) {
00776     LogDebug("") << "CaloMETAnalyzer: Could not find GT readout record" << std::endl;
00777     if (_verbose) std::cout << "CaloMETAnalyzer: Could not find GT readout record product" << std::endl;
00778   }
00779   
00780   bool bTechTriggers    = true;
00781   bool bTechTriggersAND = true;
00782   bool bTechTriggersOR  = false;
00783   bool bTechTriggersNOT = false;
00784 
00785   if (gtReadoutRecord.isValid()) {
00786     const TechnicalTriggerWord&  technicalTriggerWordBeforeMask = gtReadoutRecord->technicalTriggerWord();
00787     
00788     if (_techTrigsAND.size() == 0)
00789       bTechTriggersAND = true;
00790     else
00791       for (unsigned ttr = 0; ttr != _techTrigsAND.size(); ttr++) {
00792         bTechTriggersAND = bTechTriggersAND && technicalTriggerWordBeforeMask.at(_techTrigsAND.at(ttr));
00793       }
00794     
00795     if (_techTrigsAND.size() == 0)
00796       bTechTriggersOR = true;
00797     else
00798       for (unsigned ttr = 0; ttr != _techTrigsOR.size(); ttr++) {
00799         bTechTriggersOR = bTechTriggersOR || technicalTriggerWordBeforeMask.at(_techTrigsOR.at(ttr));
00800       }
00801     if (_techTrigsNOT.size() == 0)
00802       bTechTriggersNOT = false;
00803     else
00804       for (unsigned ttr = 0; ttr != _techTrigsNOT.size(); ttr++) {
00805         bTechTriggersNOT = bTechTriggersNOT || technicalTriggerWordBeforeMask.at(_techTrigsNOT.at(ttr));
00806       }
00807   }
00808   else
00809     {
00810       bTechTriggersAND = true;
00811       bTechTriggersOR  = true;
00812       bTechTriggersNOT = false;
00813     }
00814   
00815   if (_techTrigsAND.size()==0)
00816     bTechTriggersAND = true;
00817   if (_techTrigsOR.size()==0)
00818     bTechTriggersOR = true;
00819   if (_techTrigsNOT.size()==0)
00820     bTechTriggersNOT = false;
00821   
00822   bTechTriggers = bTechTriggersAND && bTechTriggersOR && !bTechTriggersNOT;
00823 
00824   // ==========================================================
00825   // Reconstructed MET Information - fill MonitorElements
00826   
00827   bool bHcalNoise   = bHcalNoiseFilter;
00828   bool bBeamHaloID  = bBeamHaloIDLoosePass;
00829   bool bJetID       = bJetIDMinimal;
00830 
00831   bool bPhysicsDeclared = true;
00832   if(_doHLTPhysicsOn) bPhysicsDeclared =_trig_PhysDec;
00833 
00834   if      (_tightHcalFiltering)     bHcalNoise  = bHcalNoiseFilterTight;
00835   if      (_tightBHFiltering)       bBeamHaloID = bBeamHaloIDTightPass;
00836 
00837   if      (_tightJetIDFiltering==1)  bJetID      = bJetIDMinimal;
00838   else if (_tightJetIDFiltering==2)  bJetID      = bJetIDLoose;
00839   else if (_tightJetIDFiltering==3)  bJetID      = bJetIDTight;
00840   else if (_tightJetIDFiltering==-1) bJetID      = true;
00841 
00842   bool bBasicCleanup = bTechTriggers && bPrimaryVertex && bPhysicsDeclared;
00843   bool bExtraCleanup = bBasicCleanup && bHcalNoise && bJetID && bBeamHaloID;
00844   
00845   std::string DirName = "JetMET/MET/"+_source;
00846   
00847   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); 
00848        ic != _FolderNames.end(); ic++){
00849     if (*ic=="All")                                             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00850     if (DCSFilter->filter(iEvent, iSetup)) {
00851     if (_cleanupSelection){
00852     if (*ic=="BasicCleanup" && bBasicCleanup)                   fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00853     if (*ic=="ExtraCleanup" && bExtraCleanup)                   fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00854     }
00855     if (_allSelection) {
00856       if (*ic=="HcalNoiseFilter"      && bHcalNoiseFilter )       fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00857       if (*ic=="HcalNoiseFilterTight" && bHcalNoiseFilterTight )  fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00858       if (*ic=="JetIDMinimal"         && bJetIDMinimal)           fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00859       if (*ic=="JetIDLoose"           && bJetIDLoose)             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00860       if (*ic=="JetIDTight"           && bJetIDTight)             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00861       if (*ic=="BeamHaloIDTightPass"  && bBeamHaloIDTightPass)    fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00862       if (*ic=="BeamHaloIDLoosePass"  && bBeamHaloIDLoosePass)    fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00863       if (*ic=="Triggers"             && bTechTriggers)           fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00864       if (*ic=="PV"                   && bPrimaryVertex)          fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00865     }
00866     } // DCS
00867   }
00868 }
00869 
00870   
00871 // ***********************************************************
00872 void PFMETAnalyzer::validateMET(const reco::PFMET& pfmet, 
00873                                 edm::Handle<edm::View<PFCandidate> > pfCandidates)
00874 {          
00875   double sumEx = 0;
00876   double sumEy = 0;
00877   double sumEt = 0;
00878       
00879   for( unsigned i=0; i<pfCandidates->size(); i++ ) {
00880         
00881     const reco::PFCandidate& cand = (*pfCandidates)[i];
00882         
00883     double E = cand.energy();
00884         
00886     //  if( cand.particleId()==PFCandidate::h_HF || 
00887     //      cand.particleId()==PFCandidate::egamma_HF )
00888     //    E *= hfCalibFactor_;
00889         
00890     double phi = cand.phi();
00891     double cosphi = cos(phi);
00892     double sinphi = sin(phi);
00893         
00894     double theta = cand.theta();
00895     double sintheta = sin(theta);
00896         
00897     double et = E*sintheta;
00898     double ex = et*cosphi;
00899     double ey = et*sinphi;
00900         
00901     sumEx += ex;
00902     sumEy += ey;
00903     sumEt += et;
00904   }
00905       
00906   double Et = sqrt( sumEx*sumEx + sumEy*sumEy);
00907   XYZTLorentzVector missingEt( -sumEx, -sumEy, 0, Et);
00908       
00909   if(_verbose) 
00910     if (sumEt!=pfmet.sumEt() || sumEx!=pfmet.px() || sumEy!=pfmet.py() || missingEt.T()!=pfmet.pt() )   
00911       {
00912         std::cout<<"PFSumEt: " << sumEt         <<", "<<"PFMETBlock: "<<pfmet.pt()<<std::endl;
00913         std::cout<<"PFMET: "   << missingEt.T() <<", "<<"PFMETBlock: "<<pfmet.pt()<<std::endl;
00914         std::cout<<"PFMETx: "  << missingEt.X() <<", "<<"PFMETBlockx: "<<pfmet.pt()<<std::endl;
00915         std::cout<<"PFMETy: "  << missingEt.Y() <<", "<<"PFMETBlocky: "<<pfmet.pt()<<std::endl;
00916       }
00917 }
00918 
00919 // ***********************************************************
00920 void PFMETAnalyzer::fillMESet(const edm::Event& iEvent, std::string DirName, 
00921                               const reco::PFMET& pfmet)
00922 {
00923 
00924   _dbe->setCurrentFolder(DirName);
00925 
00926   bool bLumiSecPlot=false;
00927   if (DirName.find("All")) bLumiSecPlot=true;
00928 
00929   if (_trig_JetMB)
00930     fillMonitorElement(iEvent,DirName,"",pfmet, bLumiSecPlot);
00931   if (_trig_HighPtJet)
00932     fillMonitorElement(iEvent,DirName,"HighPtJet",pfmet,false);
00933   if (_trig_LowPtJet)
00934     fillMonitorElement(iEvent,DirName,"LowPtJet",pfmet,false);
00935   if (_trig_MinBias)
00936     fillMonitorElement(iEvent,DirName,"MinBias",pfmet,false);
00937   if (_trig_HighMET)
00938     fillMonitorElement(iEvent,DirName,"HighMET",pfmet,false);
00939   if (_trig_LowMET)
00940     fillMonitorElement(iEvent,DirName,"LowMET",pfmet,false);
00941   if (_trig_Ele)
00942     fillMonitorElement(iEvent,DirName,"Ele",pfmet,false);
00943   if (_trig_Muon)
00944     fillMonitorElement(iEvent,DirName,"Muon",pfmet,false);
00945 }
00946 
00947 // ***********************************************************
00948 void PFMETAnalyzer::fillMonitorElement(const edm::Event& iEvent, std::string DirName, 
00949                                        std::string TriggerTypeName, 
00950                                        const reco::PFMET& pfmet, bool bLumiSecPlot)
00951 {
00952 
00953   if (TriggerTypeName=="HighPtJet") {
00954     if (!selectHighPtJetEvent(iEvent)) return;
00955   }
00956   else if (TriggerTypeName=="LowPtJet") {
00957     if (!selectLowPtJetEvent(iEvent)) return;
00958   }
00959   else if (TriggerTypeName=="HighMET") {
00960     if (pfmet.pt()<_highPFMETThreshold) return;
00961   }
00962   else if (TriggerTypeName=="LowMET") {
00963     if (pfmet.pt()<_lowPFMETThreshold) return;
00964   }
00965   else if (TriggerTypeName=="Ele") {
00966     if (!selectWElectronEvent(iEvent)) return;
00967   }
00968   else if (TriggerTypeName=="Muon") {
00969     if (!selectWMuonEvent(iEvent)) return;
00970   }
00971   
00972   // Reconstructed MET Information
00973   double pfSumET  = pfmet.sumEt();
00974   double pfMETSig = pfmet.mEtSig();
00975   //double pfEz     = pfmet.e_longitudinal();
00976   double pfMET    = pfmet.pt();
00977   double pfMEx    = pfmet.px();
00978   double pfMEy    = pfmet.py();
00979   double pfMETPhi = pfmet.phi();
00980 
00981   double pfNeutralEMFraction  = pfmet.NeutralEMFraction();
00982   double pfNeutralHadFraction = pfmet.NeutralHadFraction();
00983   double pfChargedEMFraction  = pfmet.ChargedEMFraction();
00984   double pfChargedHadFraction = pfmet.ChargedHadFraction();
00985   double pfMuonFraction       = pfmet.MuonFraction();
00986 
00987   
00988   //
00989   int myLuminosityBlock;
00990   //  myLuminosityBlock = (evtCounter++)/1000;
00991   myLuminosityBlock = iEvent.luminosityBlock();
00992   //
00993 
00994   if (TriggerTypeName!="") DirName = DirName +"/"+TriggerTypeName;
00995 
00996   if (_verbose) std::cout << "_etThreshold = " << _etThreshold << std::endl;
00997   if (pfSumET>_etThreshold){
00998     
00999     mePfMEx    = _dbe->get(DirName+"/"+"METTask_PfMEx");    if (mePfMEx    && mePfMEx->getRootObject())    mePfMEx->Fill(pfMEx);
01000     mePfMEy    = _dbe->get(DirName+"/"+"METTask_PfMEy");    if (mePfMEy    && mePfMEy->getRootObject())    mePfMEy->Fill(pfMEy);
01001     mePfMET    = _dbe->get(DirName+"/"+"METTask_PfMET");    if (mePfMET    && mePfMET->getRootObject())    mePfMET->Fill(pfMET);
01002     mePfMETPhi = _dbe->get(DirName+"/"+"METTask_PfMETPhi"); if (mePfMETPhi && mePfMETPhi->getRootObject()) mePfMETPhi->Fill(pfMETPhi);
01003     mePfSumET  = _dbe->get(DirName+"/"+"METTask_PfSumET");  if (mePfSumET  && mePfSumET->getRootObject())  mePfSumET->Fill(pfSumET);
01004     mePfMETSig = _dbe->get(DirName+"/"+"METTask_PfMETSig"); if (mePfMETSig && mePfMETSig->getRootObject()) mePfMETSig->Fill(pfMETSig);
01005     //mePfEz     = _dbe->get(DirName+"/"+"METTask_PfEz");     if (mePfEz     && mePfEz->getRootObject())     mePfEz->Fill(pfEz);
01006 
01007     mePfMET_logx    = _dbe->get(DirName+"/"+"METTask_PfMET_logx");    if (mePfMET_logx    && mePfMET_logx->getRootObject())    mePfMET_logx->Fill(log10(pfMET));
01008     mePfSumET_logx  = _dbe->get(DirName+"/"+"METTask_PfSumET_logx");  if (mePfSumET_logx  && mePfSumET_logx->getRootObject())  mePfSumET_logx->Fill(log10(pfSumET));
01009 
01010     //mePfMETIonFeedbck = _dbe->get(DirName+"/"+"METTask_PfMETIonFeedbck");  if (mePfMETIonFeedbck && mePfMETIonFeedbck->getRootObject()) mePfMETIonFeedbck->Fill(pfMET);
01011     //mePfMETHPDNoise   = _dbe->get(DirName+"/"+"METTask_PfMETHPDNoise");    if (mePfMETHPDNoise   && mePfMETHPDNoise->getRootObject())   mePfMETHPDNoise->Fill(pfMET);
01012     //mePfMETRBXNoise   = _dbe->get(DirName+"/"+"METTask_PfMETRBXNoise");    if (mePfMETRBXNoise   && mePfMETRBXNoise->getRootObject())   mePfMETRBXNoise->Fill(pfMET);
01013     
01014     mePfNeutralEMFraction = _dbe->get(DirName+"/"+"METTask_PfNeutralEMFraction"); 
01015     if (mePfNeutralEMFraction   && mePfNeutralEMFraction->getRootObject()) mePfNeutralEMFraction->Fill(pfNeutralEMFraction);
01016     mePfNeutralHadFraction = _dbe->get(DirName+"/"+"METTask_PfNeutralHadFraction"); 
01017     if (mePfNeutralHadFraction   && mePfNeutralHadFraction->getRootObject()) mePfNeutralHadFraction->Fill(pfNeutralHadFraction);
01018     mePfChargedEMFraction = _dbe->get(DirName+"/"+"METTask_PfChargedEMFraction"); 
01019     if (mePfChargedEMFraction   && mePfChargedEMFraction->getRootObject()) mePfChargedEMFraction->Fill(pfChargedEMFraction);
01020     mePfChargedHadFraction = _dbe->get(DirName+"/"+"METTask_PfChargedHadFraction"); 
01021     if (mePfChargedHadFraction   && mePfChargedHadFraction->getRootObject()) mePfChargedHadFraction->Fill(pfChargedHadFraction);
01022     mePfMuonFraction = _dbe->get(DirName+"/"+"METTask_PfMuonFraction"); 
01023     if (mePfMuonFraction   && mePfMuonFraction->getRootObject()) mePfMuonFraction->Fill(pfMuonFraction);
01024     
01025     if (_allhist){
01026       if (bLumiSecPlot){
01027         mePfMExLS = _dbe->get(DirName+"/"+"METTask_PfMExLS"); if (mePfMExLS && mePfMExLS->getRootObject()) mePfMExLS->Fill(pfMEx,myLuminosityBlock);
01028         mePfMEyLS = _dbe->get(DirName+"/"+"METTask_PfMEyLS"); if (mePfMEyLS && mePfMEyLS->getRootObject()) mePfMEyLS->Fill(pfMEy,myLuminosityBlock);
01029       }
01030     } // _allhist
01031 
01032 
01033     // Fill NPV profiles
01034     //--------------------------------------------------------------------------
01035     mePfMEx_profile   = _dbe->get(DirName + "/METTask_PfMEx_profile");
01036     mePfMEy_profile   = _dbe->get(DirName + "/METTask_PfMEy_profile");
01037     mePfMET_profile   = _dbe->get(DirName + "/METTask_PfMET_profile");
01038     mePfSumET_profile = _dbe->get(DirName + "/METTask_PfSumET_profile");
01039 
01040     if (mePfMEx_profile   && mePfMEx_profile  ->getRootObject()) mePfMEx_profile  ->Fill(_numPV, pfMEx);
01041     if (mePfMEy_profile   && mePfMEy_profile  ->getRootObject()) mePfMEy_profile  ->Fill(_numPV, pfMEy);
01042     if (mePfMET_profile   && mePfMET_profile  ->getRootObject()) mePfMET_profile  ->Fill(_numPV, pfMET);
01043     if (mePfSumET_profile && mePfSumET_profile->getRootObject()) mePfSumET_profile->Fill(_numPV, pfSumET);
01044 
01045     mePfNeutralEMFraction_profile  = _dbe->get(DirName + "/METTask_PfNeutralEMFraction_profile");
01046     mePfNeutralHadFraction_profile = _dbe->get(DirName + "/METTask_PfNeutralHadFraction_profile");
01047     mePfChargedEMFraction_profile  = _dbe->get(DirName + "/METTask_PfChargedEMFraction_profile");
01048     mePfChargedHadFraction_profile = _dbe->get(DirName + "/METTask_PfChargedHadFraction_profile");
01049     mePfMuonFraction_profile       = _dbe->get(DirName + "/METTask_PfMuonFraction_profile");
01050 
01051     if (mePfNeutralEMFraction_profile  && mePfNeutralEMFraction_profile ->getRootObject()) mePfNeutralEMFraction_profile ->Fill(_numPV, pfNeutralEMFraction);
01052     if (mePfNeutralHadFraction_profile && mePfNeutralHadFraction_profile->getRootObject()) mePfNeutralHadFraction_profile->Fill(_numPV, pfNeutralHadFraction);
01053     if (mePfChargedEMFraction_profile  && mePfChargedEMFraction_profile ->getRootObject()) mePfChargedEMFraction_profile ->Fill(_numPV, pfChargedEMFraction);
01054     if (mePfChargedHadFraction_profile && mePfChargedHadFraction_profile->getRootObject()) mePfChargedHadFraction_profile->Fill(_numPV, pfChargedHadFraction);
01055     if (mePfMuonFraction_profile       && mePfMuonFraction_profile      ->getRootObject()) mePfMuonFraction_profile      ->Fill(_numPV, pfMuonFraction);
01056 
01057   } // ET threshold cut
01058 }
01059 
01060 // ***********************************************************
01061 bool PFMETAnalyzer::selectHighPtJetEvent(const edm::Event& iEvent){
01062 
01063   bool return_value=false;
01064 
01065   edm::Handle<reco::PFJetCollection> pfJets;
01066   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
01067   if (!pfJets.isValid()) {
01068     LogDebug("") << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
01069     if (_verbose) std::cout << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
01070   }
01071 
01072   for (reco::PFJetCollection::const_iterator pf = pfJets->begin(); 
01073        pf!=pfJets->end(); ++pf){
01074     if (pf->pt()>_highPtPFJetThreshold){
01075       return_value=true;
01076     }
01077   }
01078   
01079   return return_value;
01080 }
01081 
01082 // // ***********************************************************
01083 bool PFMETAnalyzer::selectLowPtJetEvent(const edm::Event& iEvent){
01084 
01085   bool return_value=false;
01086 
01087   edm::Handle<reco::PFJetCollection> pfJets;
01088   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
01089   if (!pfJets.isValid()) {
01090     LogDebug("") << "PFMETAnalyzer: Could not find jet product" << std::endl;
01091     if (_verbose) std::cout << "PFMETAnalyzer: Could not find jet product" << std::endl;
01092   }
01093 
01094   for (reco::PFJetCollection::const_iterator cal = pfJets->begin(); 
01095        cal!=pfJets->end(); ++cal){
01096     if (cal->pt()>_lowPtPFJetThreshold){
01097       return_value=true;
01098     }
01099   }
01100 
01101   return return_value;
01102 
01103 }
01104 
01105 // ***********************************************************
01106 bool PFMETAnalyzer::selectWElectronEvent(const edm::Event& iEvent){
01107 
01108   bool return_value=true;
01109 
01110   /*
01111     W-electron event selection comes here
01112   */
01113 
01114   return return_value;
01115 
01116 }
01117 
01118 // ***********************************************************
01119 bool PFMETAnalyzer::selectWMuonEvent(const edm::Event& iEvent){
01120 
01121   bool return_value=true;
01122 
01123   /*
01124     W-muon event selection comes here
01125   */
01126 
01127   return return_value;
01128 
01129 }
01130
/data/refman/pasoursint/CMSSW_5_3_0/src/DQMOffline/JetMET/src/PFMETAnalyzer.cc
