CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/DQMOffline/JetMET/src/PFMETAnalyzer.cc

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