Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2011/10/10 14:43:31 $
00005  *  $Revision: 1.40 $
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 // ***********************************************************
00308 void PFMETAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00309 {
00310   if ( _HighPtJetEventFlag->on() ) _HighPtJetEventFlag->initRun( iRun, iSetup );
00311   if ( _LowPtJetEventFlag ->on() ) _LowPtJetEventFlag ->initRun( iRun, iSetup );
00312   if ( _MinBiasEventFlag  ->on() ) _MinBiasEventFlag  ->initRun( iRun, iSetup );
00313   if ( _HighMETEventFlag  ->on() ) _HighMETEventFlag  ->initRun( iRun, iSetup );
00314   if ( _LowMETEventFlag   ->on() ) _LowMETEventFlag   ->initRun( iRun, iSetup );
00315   if ( _EleEventFlag      ->on() ) _EleEventFlag      ->initRun( iRun, iSetup );
00316   if ( _MuonEventFlag     ->on() ) _MuonEventFlag     ->initRun( iRun, iSetup );
00317 
00318   if (_HighPtJetEventFlag->on() && _HighPtJetEventFlag->expressionsFromDB(_HighPtJetEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00319     highPtJetExpr_ = _HighPtJetEventFlag->expressionsFromDB(_HighPtJetEventFlag->hltDBKey(), iSetup);
00320   if (_LowPtJetEventFlag->on() && _LowPtJetEventFlag->expressionsFromDB(_LowPtJetEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00321     lowPtJetExpr_  = _LowPtJetEventFlag->expressionsFromDB(_LowPtJetEventFlag->hltDBKey(),   iSetup);
00322   if (_HighMETEventFlag->on() && _HighMETEventFlag->expressionsFromDB(_HighMETEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00323     highMETExpr_   = _HighMETEventFlag->expressionsFromDB(_HighMETEventFlag->hltDBKey(),     iSetup);
00324   if (_LowMETEventFlag->on() && _LowMETEventFlag->expressionsFromDB(_LowMETEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00325     lowMETExpr_    = _LowMETEventFlag->expressionsFromDB(_LowMETEventFlag->hltDBKey(),       iSetup);
00326   if (_MuonEventFlag->on() && _MuonEventFlag->expressionsFromDB(_MuonEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00327     muonExpr_      = _MuonEventFlag->expressionsFromDB(_MuonEventFlag->hltDBKey(),           iSetup);
00328   if (_EleEventFlag->on() && _EleEventFlag->expressionsFromDB(_EleEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00329     elecExpr_      = _EleEventFlag->expressionsFromDB(_EleEventFlag->hltDBKey(),             iSetup);
00330   if (_MinBiasEventFlag->on() && _MinBiasEventFlag->expressionsFromDB(_MinBiasEventFlag->hltDBKey(), iSetup)[0] != "CONFIG_ERROR")
00331     minbiasExpr_   = _MinBiasEventFlag->expressionsFromDB(_MinBiasEventFlag->hltDBKey(),     iSetup);
00332 
00333 }
00334 
00335 // ***********************************************************
00336 void PFMETAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup, DQMStore * dbe)
00337 {
00338   
00339   //
00340   //--- Check the time length of the Run from the lumi section plots
00341 
00342   std::string dirName = "JetMET/MET/"+_source+"/";
00343   _dbe->setCurrentFolder(dirName);
00344 
00345   TH1F* tlumisec;
00346 
00347   MonitorElement *meLumiSec = _dbe->get("aaa");
00348   meLumiSec = _dbe->get("JetMET/lumisec");
00349 
00350   int totlsec=0;
00351   double totltime=0.;
00352   if ( meLumiSec->getRootObject() ) {
00353     tlumisec = meLumiSec->getTH1F();
00354     for (int i=0; i<500; i++){
00355       if (tlumisec->GetBinContent(i+1)) totlsec++;
00356     }
00357     totltime = double(totlsec*90); // one lumi sec ~ 90 (sec)
00358   }
00359 
00360   if (totltime==0.) totltime=1.; 
00361 
00362   //
00363   //--- Make the integrated plots with rate (Hz)
00364 
00365   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); ic != _FolderNames.end(); ic++)
00366     {
00367 
00368       std::string DirName;
00369       DirName = dirName+*ic;
00370 
00371       makeRatePlot(DirName,totltime);
00372       if ( _HighPtJetEventFlag->on() ) 
00373         makeRatePlot(DirName+"/"+"triggerName_HighJetPt",totltime);
00374       if ( _LowPtJetEventFlag->on() ) 
00375         makeRatePlot(DirName+"/"+"triggerName_LowJetPt",totltime);
00376       if ( _MinBiasEventFlag->on() ) 
00377         makeRatePlot(DirName+"/"+"triggerName_MinBias",totltime);
00378       if ( _HighMETEventFlag->on() ) 
00379         makeRatePlot(DirName+"/"+"triggerName_HighMET",totltime);
00380       if ( _LowMETEventFlag->on() ) 
00381         makeRatePlot(DirName+"/"+"triggerName_LowMET",totltime);
00382       if ( _EleEventFlag->on() ) 
00383         makeRatePlot(DirName+"/"+"triggerName_Ele",totltime);
00384       if ( _MuonEventFlag->on() ) 
00385         makeRatePlot(DirName+"/"+"triggerName_Muon",totltime);
00386     }
00387 }
00388 
00389 // ***********************************************************
00390 void PFMETAnalyzer::makeRatePlot(std::string DirName, double totltime)
00391 {
00392 
00393   _dbe->setCurrentFolder(DirName);
00394   MonitorElement *mePfMET = _dbe->get(DirName+"/"+"METTask_PfMET");
00395 
00396   TH1F* tPfMET;
00397   TH1F* tPfMETRate;
00398 
00399   if ( mePfMET )
00400     if ( mePfMET->getRootObject() ) {
00401       tPfMET     = mePfMET->getTH1F();
00402       
00403       // Integral plot & convert number of events to rate (hz)
00404       tPfMETRate = (TH1F*) tPfMET->Clone("METTask_PfMETRate");
00405       for (int i = tPfMETRate->GetNbinsX()-1; i>=0; i--){
00406         tPfMETRate->SetBinContent(i+1,tPfMETRate->GetBinContent(i+2)+tPfMET->GetBinContent(i+1));
00407       }
00408       for (int i = 0; i<tPfMETRate->GetNbinsX(); i++){
00409         tPfMETRate->SetBinContent(i+1,tPfMETRate->GetBinContent(i+1)/double(totltime));
00410       }
00411 
00412       tPfMETRate->SetName("METTask_PfMETRate");
00413       tPfMETRate->SetTitle("METTask_PfMETRate");
00414       mePfMETRate      = _dbe->book1D("METTask_PfMETRate",tPfMETRate);
00415       
00416     }
00417 }
00418 
00419 // ***********************************************************
00420 void PFMETAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
00421                             const edm::TriggerResults& triggerResults) {
00422 
00423   if (_verbose) std::cout << "PfMETAnalyzer analyze" << std::endl;
00424 
00425   LogTrace(metname)<<"[PFMETAnalyzer] Analyze PFMET";
00426 
00427   metME->Fill(3);
00428 
00429   // ==========================================================  
00430   // Trigger information 
00431   //
00432   _trig_JetMB=0;
00433   _trig_HighPtJet=0;
00434   _trig_LowPtJet=0;
00435   _trig_MinBias=0;
00436   _trig_HighMET=0;
00437   _trig_LowMET=0;
00438   _trig_Ele=0;
00439   _trig_Muon=0;
00440   _trig_PhysDec=0;
00441   if(&triggerResults) {   
00442 
00444 
00445     //
00446     //
00447     // Check how many HLT triggers are in triggerResults 
00448     int ntrigs = triggerResults.size();
00449     if (_verbose) std::cout << "ntrigs=" << ntrigs << std::endl;
00450     
00451     //
00452     //
00453     // If index=ntrigs, this HLT trigger doesn't exist in the HLT table for this data.
00454     const edm::TriggerNames & triggerNames = iEvent.triggerNames(triggerResults);
00455     
00456     //
00457     //
00458     const unsigned int nTrig(triggerNames.size());
00459     for (unsigned int i=0;i<nTrig;++i)
00460       {
00461         if (triggerNames.triggerName(i).find(highPtJetExpr_[0].substr(0,highPtJetExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00462           _trig_HighPtJet=true;
00463         else if (triggerNames.triggerName(i).find(lowPtJetExpr_[0].substr(0,lowPtJetExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00464           _trig_LowPtJet=true;
00465         else if (triggerNames.triggerName(i).find(highMETExpr_[0].substr(0,highMETExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00466           _trig_HighMET=true;
00467         else if (triggerNames.triggerName(i).find(lowMETExpr_[0].substr(0,lowMETExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00468           _trig_LowMET=true;
00469         else if (triggerNames.triggerName(i).find(muonExpr_[0].substr(0,muonExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00470           _trig_Muon=true;
00471         else if (triggerNames.triggerName(i).find(elecExpr_[0].substr(0,elecExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00472           _trig_Ele=true;
00473         else if (triggerNames.triggerName(i).find(minbiasExpr_[0].substr(0,minbiasExpr_[0].rfind("_v")+2))!=std::string::npos && triggerResults.accept(i))
00474           _trig_MinBias=true;
00475       }
00476 
00477       // count number of requested Jet or MB HLT paths which have fired
00478     for (unsigned int i=0; i!=HLTPathsJetMBByName_.size(); i++) {
00479       unsigned int triggerIndex = triggerNames.triggerIndex(HLTPathsJetMBByName_[i]);
00480       if (triggerIndex<triggerResults.size()) {
00481         if (triggerResults.accept(triggerIndex)) {
00482           _trig_JetMB++;
00483         }
00484       }
00485     }
00486     // for empty input vectors (n==0), take all HLT triggers!
00487     if (HLTPathsJetMBByName_.size()==0) _trig_JetMB=triggerResults.size()-1;
00488 
00489     //
00490     /*
00491       if ( _HighPtJetEventFlag->on() && _HighPtJetEventFlag->accept( iEvent, iSetup) )
00492       _trig_HighPtJet=1;
00493       
00494       if ( _LowPtJetEventFlag->on() && _LowPtJetEventFlag->accept( iEvent, iSetup) )
00495       _trig_LowPtJet=1;
00496       
00497       if ( _MinBiasEventFlag->on() && _MinBiasEventFlag->accept( iEvent, iSetup) )
00498       _trig_MinBias=1;
00499       
00500       if ( _HighMETEventFlag->on() && _HighMETEventFlag->accept( iEvent, iSetup) )
00501       _trig_HighMET=1;
00502       
00503       if ( _LowMETEventFlag->on() && _LowMETEventFlag->accept( iEvent, iSetup) )
00504       _trig_LowMET=1;
00505       
00506       if ( _EleEventFlag->on() && _EleEventFlag->accept( iEvent, iSetup) )
00507       _trig_Ele=1;
00508       
00509       if ( _MuonEventFlag->on() && _MuonEventFlag->accept( iEvent, iSetup) )
00510       _trig_Muon=1;
00511     */
00512     
00513     if (triggerNames.triggerIndex(_hlt_PhysDec)   != triggerNames.size() &&
00514         triggerResults.accept(triggerNames.triggerIndex(_hlt_PhysDec)))   _trig_PhysDec=1;
00515   } else {
00516 
00517     edm::LogInfo("PFMetAnalyzer") << "TriggerResults::HLT not found, "
00518       "automatically select events"; 
00519 
00520     // TriggerResults object not found. Look at all events.    
00521     _trig_JetMB=1;
00522   }
00523 
00524   // ==========================================================
00525   // PfMET information
00526   
00527   // **** Get the MET container  
00528   edm::Handle<reco::PFMETCollection> pfmetcoll;
00529   iEvent.getByLabel(thePfMETCollectionLabel, pfmetcoll);
00530   
00531   if(!pfmetcoll.isValid()) return;
00532 
00533   const PFMETCollection *pfmetcol = pfmetcoll.product();
00534   const PFMET *pfmet;
00535   pfmet = &(pfmetcol->front());
00536     
00537   LogTrace(metname)<<"[PfMETAnalyzer] Call to the PfMET analyzer";
00538 
00539   // ==========================================================
00540   //
00541   edm::Handle<HcalNoiseRBXCollection> HRBXCollection;
00542   iEvent.getByLabel(HcalNoiseRBXCollectionTag,HRBXCollection);
00543   if (!HRBXCollection.isValid()) {
00544     LogDebug("") << "PfMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00545     if (_verbose) std::cout << "PfMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00546   }
00547   
00548   edm::Handle<HcalNoiseSummary> HNoiseSummary;
00549   iEvent.getByLabel(HcalNoiseSummaryTag,HNoiseSummary);
00550   if (!HNoiseSummary.isValid()) {
00551     LogDebug("") << "PfMETAnalyzer: Could not find Hcal NoiseSummary product" << std::endl;
00552     if (_verbose) std::cout << "PfMETAnalyzer: Could not find Hcal NoiseSummary product" << std::endl;
00553   }
00554 
00555   edm::Handle<reco::CaloJetCollection> caloJets;
00556   iEvent.getByLabel(theJetCollectionLabel, caloJets);
00557   if (!caloJets.isValid()) {
00558     LogDebug("") << "PFMETAnalyzer: Could not find jet product" << std::endl;
00559     if (_verbose) std::cout << "PFMETAnalyzer: Could not find jet product" << std::endl;
00560   }
00561 
00562   edm::Handle<edm::View<PFCandidate> > pfCandidates;
00563   iEvent.getByLabel(PFCandidatesTag, pfCandidates);
00564   if (!pfCandidates.isValid()) {
00565     LogDebug("") << "PfMETAnalyzer: Could not find pfcandidates product" << std::endl;
00566     if (_verbose) std::cout << "PfMETAnalyzer: Could not find pfcandidates product" << std::endl;
00567   }
00568 
00569   edm::Handle<reco::PFJetCollection> pfJets;
00570   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
00571   if (!pfJets.isValid()) {
00572     LogDebug("") << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
00573     if (_verbose) std::cout << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
00574   }
00575   // ==========================================================
00576   // PfMET sanity check
00577 
00578   if (_source=="PfMET") validateMET(*pfmet, pfCandidates);
00579   
00580   // ==========================================================
00581   // JetID 
00582 
00583   if (_verbose) std::cout << "JetID starts" << std::endl;
00584   
00585   //
00586   // --- Minimal cuts
00587   //
00588   bool bJetIDMinimal=true;
00589   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00590        cal!=caloJets->end(); ++cal){
00591     jetID->calculate(iEvent, *cal);
00592     if (cal->pt()>10.){
00593       if (fabs(cal->eta())<=2.6 && 
00594           cal->emEnergyFraction()<=0.01) bJetIDMinimal=false;
00595     }
00596   }
00597 
00598   //
00599   // --- Loose cuts, not PF specific for now!
00600   //
00601   bool bJetIDLoose=true;
00602   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00603        cal!=caloJets->end(); ++cal){ 
00604     jetID->calculate(iEvent, *cal);
00605     if (_verbose) std::cout << jetID->n90Hits() << " " 
00606                             << jetID->restrictedEMF() << " "
00607                             << cal->pt() << std::endl;
00608     if (cal->pt()>10.){
00609       //
00610       // for all regions
00611       if (jetID->n90Hits()<2)  bJetIDLoose=false; 
00612       if (jetID->fHPD()>=0.98) bJetIDLoose=false; 
00613       //if (jetID->restrictedEMF()<0.01) bJetIDLoose=false; 
00614       //
00615       // for non-forward
00616       if (fabs(cal->eta())<2.55){
00617         if (cal->emEnergyFraction()<=0.01) bJetIDLoose=false; 
00618       }
00619       // for forward
00620       else {
00621         if (cal->emEnergyFraction()<=-0.9) bJetIDLoose=false; 
00622         if (cal->pt()>80.){
00623           if (cal->emEnergyFraction()>= 1.0) bJetIDLoose=false; 
00624         }
00625       } // forward vs non-forward
00626     }   // pt>10 GeV/c
00627   }     // calor-jets loop
00628  
00629   //
00630   // --- Tight cuts
00631   //
00632   bool bJetIDTight=true;
00633   bJetIDTight=bJetIDLoose;
00634   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00635        cal!=caloJets->end(); ++cal){
00636     jetID->calculate(iEvent, *cal);
00637     if (cal->pt()>25.){
00638       //
00639       // for all regions
00640       if (jetID->fHPD()>=0.95) bJetIDTight=false; 
00641       //
00642       // for 1.0<|eta|<1.75
00643       if (fabs(cal->eta())>=1.00 && fabs(cal->eta())<1.75){
00644         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00645       }
00646       //
00647       // for 1.75<|eta|<2.55
00648       else if (fabs(cal->eta())>=1.75 && fabs(cal->eta())<2.55){
00649         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00650       }
00651       //
00652       // for 2.55<|eta|<3.25
00653       else if (fabs(cal->eta())>=2.55 && fabs(cal->eta())<3.25){
00654         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3) bJetIDTight=false; 
00655         if (cal->pt()>=50. && cal->pt()< 80. && cal->emEnergyFraction()<=-0.2) bJetIDTight=false; 
00656         if (cal->pt()>=80. && cal->pt()<340. && cal->emEnergyFraction()<=-0.1) bJetIDTight=false; 
00657         if (cal->pt()>=340.                  && cal->emEnergyFraction()<=-0.1 
00658             && cal->emEnergyFraction()>=0.95) bJetIDTight=false; 
00659       }
00660       //
00661       // for 3.25<|eta|
00662       else if (fabs(cal->eta())>=3.25){
00663         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3
00664             && cal->emEnergyFraction()>=0.90) bJetIDTight=false; 
00665         if (cal->pt()>=50. && cal->pt()<130. && cal->emEnergyFraction()<=-0.2
00666             && cal->emEnergyFraction()>=0.80) bJetIDTight=false; 
00667         if (cal->pt()>=130.                  && cal->emEnergyFraction()<=-0.1 
00668             && cal->emEnergyFraction()>=0.70) bJetIDTight=false; 
00669       }
00670     }   // pt>10 GeV/c
00671   }     // calor-jets loop
00672   
00673   if (_verbose) std::cout << "JetID ends" << std::endl;
00674 
00675 
00676   // ==========================================================
00677   // HCAL Noise filter
00678   
00679   bool bHcalNoiseFilter      = HNoiseSummary->passLooseNoiseFilter();
00680   bool bHcalNoiseFilterTight = HNoiseSummary->passTightNoiseFilter();
00681 
00682   // ==========================================================
00683   // Get BeamHaloSummary
00684   edm::Handle<BeamHaloSummary> TheBeamHaloSummary ;
00685   iEvent.getByLabel(BeamHaloSummaryTag, TheBeamHaloSummary) ;
00686 
00687   bool bBeamHaloIDTightPass = true;
00688   bool bBeamHaloIDLoosePass = true;
00689 
00690   if(!TheBeamHaloSummary.isValid()) {
00691 
00692     const BeamHaloSummary TheSummary = (*TheBeamHaloSummary.product() );
00693 
00694     if( !TheSummary.EcalLooseHaloId()  && !TheSummary.HcalLooseHaloId() && 
00695         !TheSummary.CSCLooseHaloId()   && !TheSummary.GlobalLooseHaloId() )
00696       bBeamHaloIDLoosePass = false;
00697 
00698     if( !TheSummary.EcalTightHaloId()  && !TheSummary.HcalTightHaloId() && 
00699         !TheSummary.CSCTightHaloId()   && !TheSummary.GlobalTightHaloId() )
00700       bBeamHaloIDTightPass = false;
00701 
00702   }
00703 
00704   // ==========================================================
00705   //Vertex information
00706   
00707   bool bPrimaryVertex = true;
00708   if(_doPVCheck){
00709     bPrimaryVertex = false;
00710     Handle<VertexCollection> vertexHandle;
00711 
00712     iEvent.getByLabel(vertexTag, vertexHandle);
00713 
00714     if (!vertexHandle.isValid()) {
00715       LogDebug("") << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00716       if (_verbose) std::cout << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00717     }
00718     
00719     if ( vertexHandle.isValid() ){
00720       VertexCollection vertexCollection = *(vertexHandle.product());
00721       int vertex_number     = vertexCollection.size();
00722       VertexCollection::const_iterator v = vertexCollection.begin();
00723       for ( ; v != vertexCollection.end(); ++v) {
00724         double vertex_chi2    = v->normalizedChi2();
00725         double vertex_ndof    = v->ndof();
00726         bool   fakeVtx        = v->isFake();
00727         double vertex_Z       = v->z();
00728         
00729         if (  !fakeVtx
00730               && vertex_number>=_nvtx_min
00731               && vertex_ndof   >_vtxndof_min
00732               && vertex_chi2   <_vtxchi2_max
00733               && fabs(vertex_Z)<_vtxz_max )
00734           bPrimaryVertex = true;
00735       }
00736     }
00737   }
00738   // ==========================================================
00739 
00740   edm::Handle< L1GlobalTriggerReadoutRecord > gtReadoutRecord;
00741   iEvent.getByLabel( gtTag, gtReadoutRecord);
00742 
00743   if (!gtReadoutRecord.isValid()) {
00744     LogDebug("") << "CaloMETAnalyzer: Could not find GT readout record" << std::endl;
00745     if (_verbose) std::cout << "CaloMETAnalyzer: Could not find GT readout record product" << std::endl;
00746   }
00747   
00748   bool bTechTriggers    = true;
00749   bool bTechTriggersAND = true;
00750   bool bTechTriggersOR  = false;
00751   bool bTechTriggersNOT = false;
00752 
00753   if (gtReadoutRecord.isValid()) {
00754     const TechnicalTriggerWord&  technicalTriggerWordBeforeMask = gtReadoutRecord->technicalTriggerWord();
00755     
00756     if (_techTrigsAND.size() == 0)
00757       bTechTriggersAND = true;
00758     else
00759       for (unsigned ttr = 0; ttr != _techTrigsAND.size(); ttr++) {
00760         bTechTriggersAND = bTechTriggersAND && technicalTriggerWordBeforeMask.at(_techTrigsAND.at(ttr));
00761       }
00762     
00763     if (_techTrigsAND.size() == 0)
00764       bTechTriggersOR = true;
00765     else
00766       for (unsigned ttr = 0; ttr != _techTrigsOR.size(); ttr++) {
00767         bTechTriggersOR = bTechTriggersOR || technicalTriggerWordBeforeMask.at(_techTrigsOR.at(ttr));
00768       }
00769     if (_techTrigsNOT.size() == 0)
00770       bTechTriggersNOT = false;
00771     else
00772       for (unsigned ttr = 0; ttr != _techTrigsNOT.size(); ttr++) {
00773         bTechTriggersNOT = bTechTriggersNOT || technicalTriggerWordBeforeMask.at(_techTrigsNOT.at(ttr));
00774       }
00775   }
00776   else
00777     {
00778       bTechTriggersAND = true;
00779       bTechTriggersOR  = true;
00780       bTechTriggersNOT = false;
00781     }
00782   
00783   if (_techTrigsAND.size()==0)
00784     bTechTriggersAND = true;
00785   if (_techTrigsOR.size()==0)
00786     bTechTriggersOR = true;
00787   if (_techTrigsNOT.size()==0)
00788     bTechTriggersNOT = false;
00789   
00790   bTechTriggers = bTechTriggersAND && bTechTriggersOR && !bTechTriggersNOT;
00791 
00792   // ==========================================================
00793   // Reconstructed MET Information - fill MonitorElements
00794   
00795   bool bHcalNoise   = bHcalNoiseFilter;
00796   bool bBeamHaloID  = bBeamHaloIDLoosePass;
00797   bool bJetID       = bJetIDMinimal;
00798 
00799   bool bPhysicsDeclared = true;
00800   if(_doHLTPhysicsOn) bPhysicsDeclared =_trig_PhysDec;
00801 
00802   if      (_tightHcalFiltering)     bHcalNoise  = bHcalNoiseFilterTight;
00803   if      (_tightBHFiltering)       bBeamHaloID = bBeamHaloIDTightPass;
00804 
00805   if      (_tightJetIDFiltering==1)  bJetID      = bJetIDMinimal;
00806   else if (_tightJetIDFiltering==2)  bJetID      = bJetIDLoose;
00807   else if (_tightJetIDFiltering==3)  bJetID      = bJetIDTight;
00808   else if (_tightJetIDFiltering==-1) bJetID      = true;
00809 
00810   bool bBasicCleanup = bTechTriggers && bPrimaryVertex && bPhysicsDeclared;
00811   bool bExtraCleanup = bBasicCleanup && bHcalNoise && bJetID && bBeamHaloID;
00812   
00813   std::string DirName = "JetMET/MET/"+_source;
00814   
00815   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); 
00816        ic != _FolderNames.end(); ic++){
00817     if (*ic=="All")                                             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00818     if (DCSFilter->filter(iEvent, iSetup)) {
00819     if (_cleanupSelection){
00820     if (*ic=="BasicCleanup" && bBasicCleanup)                   fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00821     if (*ic=="ExtraCleanup" && bExtraCleanup)                   fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00822     }
00823     if (_allSelection) {
00824       if (*ic=="HcalNoiseFilter"      && bHcalNoiseFilter )       fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00825       if (*ic=="HcalNoiseFilterTight" && bHcalNoiseFilterTight )  fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00826       if (*ic=="JetIDMinimal"         && bJetIDMinimal)           fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00827       if (*ic=="JetIDLoose"           && bJetIDLoose)             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00828       if (*ic=="JetIDTight"           && bJetIDTight)             fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00829       if (*ic=="BeamHaloIDTightPass"  && bBeamHaloIDTightPass)    fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00830       if (*ic=="BeamHaloIDLoosePass"  && bBeamHaloIDLoosePass)    fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00831       if (*ic=="Triggers"             && bTechTriggers)           fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00832       if (*ic=="PV"                   && bPrimaryVertex)          fillMESet(iEvent, DirName+"/"+*ic, *pfmet);
00833     }
00834     } // DCS
00835   }
00836 }
00837 
00838   
00839 // ***********************************************************
00840 void PFMETAnalyzer::validateMET(const reco::PFMET& pfmet, 
00841                                 edm::Handle<edm::View<PFCandidate> > pfCandidates)
00842 {          
00843   double sumEx = 0;
00844   double sumEy = 0;
00845   double sumEt = 0;
00846       
00847   for( unsigned i=0; i<pfCandidates->size(); i++ ) {
00848         
00849     const reco::PFCandidate& cand = (*pfCandidates)[i];
00850         
00851     double E = cand.energy();
00852         
00854     //  if( cand.particleId()==PFCandidate::h_HF || 
00855     //      cand.particleId()==PFCandidate::egamma_HF )
00856     //    E *= hfCalibFactor_;
00857         
00858     double phi = cand.phi();
00859     double cosphi = cos(phi);
00860     double sinphi = sin(phi);
00861         
00862     double theta = cand.theta();
00863     double sintheta = sin(theta);
00864         
00865     double et = E*sintheta;
00866     double ex = et*cosphi;
00867     double ey = et*sinphi;
00868         
00869     sumEx += ex;
00870     sumEy += ey;
00871     sumEt += et;
00872   }
00873       
00874   double Et = sqrt( sumEx*sumEx + sumEy*sumEy);
00875   XYZTLorentzVector missingEt( -sumEx, -sumEy, 0, Et);
00876       
00877   if(_verbose) 
00878     if (sumEt!=pfmet.sumEt() || sumEx!=pfmet.px() || sumEy!=pfmet.py() || missingEt.T()!=pfmet.pt() )   
00879       {
00880         std::cout<<"PFSumEt: " << sumEt         <<", "<<"PFMETBlock: "<<pfmet.pt()<<std::endl;
00881         std::cout<<"PFMET: "   << missingEt.T() <<", "<<"PFMETBlock: "<<pfmet.pt()<<std::endl;
00882         std::cout<<"PFMETx: "  << missingEt.X() <<", "<<"PFMETBlockx: "<<pfmet.pt()<<std::endl;
00883         std::cout<<"PFMETy: "  << missingEt.Y() <<", "<<"PFMETBlocky: "<<pfmet.pt()<<std::endl;
00884       }
00885 }
00886 
00887 // ***********************************************************
00888 void PFMETAnalyzer::fillMESet(const edm::Event& iEvent, std::string DirName, 
00889                               const reco::PFMET& pfmet)
00890 {
00891 
00892   _dbe->setCurrentFolder(DirName);
00893 
00894   bool bLumiSecPlot=false;
00895   if (DirName.find("All")) bLumiSecPlot=true;
00896 
00897   if (_trig_JetMB)
00898     fillMonitorElement(iEvent,DirName,"",pfmet, bLumiSecPlot);
00899   if (_trig_HighPtJet)
00900     fillMonitorElement(iEvent,DirName,"HighPtJet",pfmet,false);
00901   if (_trig_LowPtJet)
00902     fillMonitorElement(iEvent,DirName,"LowPtJet",pfmet,false);
00903   if (_trig_MinBias)
00904     fillMonitorElement(iEvent,DirName,"MinBias",pfmet,false);
00905   if (_trig_HighMET)
00906     fillMonitorElement(iEvent,DirName,"HighMET",pfmet,false);
00907   if (_trig_LowMET)
00908     fillMonitorElement(iEvent,DirName,"LowMET",pfmet,false);
00909   if (_trig_Ele)
00910     fillMonitorElement(iEvent,DirName,"Ele",pfmet,false);
00911   if (_trig_Muon)
00912     fillMonitorElement(iEvent,DirName,"Muon",pfmet,false);
00913 }
00914 
00915 // ***********************************************************
00916 void PFMETAnalyzer::fillMonitorElement(const edm::Event& iEvent, std::string DirName, 
00917                                        std::string TriggerTypeName, 
00918                                        const reco::PFMET& pfmet, bool bLumiSecPlot)
00919 {
00920 
00921   if (TriggerTypeName=="HighPtJet") {
00922     if (!selectHighPtJetEvent(iEvent)) return;
00923   }
00924   else if (TriggerTypeName=="LowPtJet") {
00925     if (!selectLowPtJetEvent(iEvent)) return;
00926   }
00927   else if (TriggerTypeName=="HighMET") {
00928     if (pfmet.pt()<_highPFMETThreshold) return;
00929   }
00930   else if (TriggerTypeName=="LowMET") {
00931     if (pfmet.pt()<_lowPFMETThreshold) return;
00932   }
00933   else if (TriggerTypeName=="Ele") {
00934     if (!selectWElectronEvent(iEvent)) return;
00935   }
00936   else if (TriggerTypeName=="Muon") {
00937     if (!selectWMuonEvent(iEvent)) return;
00938   }
00939   
00940   // Reconstructed MET Information
00941   double pfSumET  = pfmet.sumEt();
00942   double pfMETSig = pfmet.mEtSig();
00943   //double pfEz     = pfmet.e_longitudinal();
00944   double pfMET    = pfmet.pt();
00945   double pfMEx    = pfmet.px();
00946   double pfMEy    = pfmet.py();
00947   double pfMETPhi = pfmet.phi();
00948 
00949   double pfNeutralEMFraction  = pfmet.NeutralEMFraction();
00950   double pfNeutralHadFraction = pfmet.NeutralHadFraction();
00951   double pfChargedEMFraction  = pfmet.ChargedEMFraction();
00952   double pfChargedHadFraction = pfmet.ChargedHadFraction();
00953   double pfMuonFraction       = pfmet.MuonFraction();
00954 
00955   
00956   //
00957   int myLuminosityBlock;
00958   //  myLuminosityBlock = (evtCounter++)/1000;
00959   myLuminosityBlock = iEvent.luminosityBlock();
00960   //
00961 
00962   if (TriggerTypeName!="") DirName = DirName +"/"+TriggerTypeName;
00963 
00964   if (_verbose) std::cout << "_etThreshold = " << _etThreshold << std::endl;
00965   if (pfSumET>_etThreshold){
00966     
00967     mePfMEx    = _dbe->get(DirName+"/"+"METTask_PfMEx");    if (mePfMEx    && mePfMEx->getRootObject())    mePfMEx->Fill(pfMEx);
00968     mePfMEy    = _dbe->get(DirName+"/"+"METTask_PfMEy");    if (mePfMEy    && mePfMEy->getRootObject())    mePfMEy->Fill(pfMEy);
00969     mePfMET    = _dbe->get(DirName+"/"+"METTask_PfMET");    if (mePfMET    && mePfMET->getRootObject())    mePfMET->Fill(pfMET);
00970     mePfMETPhi = _dbe->get(DirName+"/"+"METTask_PfMETPhi"); if (mePfMETPhi && mePfMETPhi->getRootObject()) mePfMETPhi->Fill(pfMETPhi);
00971     mePfSumET  = _dbe->get(DirName+"/"+"METTask_PfSumET");  if (mePfSumET  && mePfSumET->getRootObject())  mePfSumET->Fill(pfSumET);
00972     mePfMETSig = _dbe->get(DirName+"/"+"METTask_PfMETSig"); if (mePfMETSig && mePfMETSig->getRootObject()) mePfMETSig->Fill(pfMETSig);
00973     //mePfEz     = _dbe->get(DirName+"/"+"METTask_PfEz");     if (mePfEz     && mePfEz->getRootObject())     mePfEz->Fill(pfEz);
00974 
00975     mePfMET_logx    = _dbe->get(DirName+"/"+"METTask_PfMET_logx");    if (mePfMET_logx    && mePfMET_logx->getRootObject())    mePfMET_logx->Fill(log10(pfMET));
00976     mePfSumET_logx  = _dbe->get(DirName+"/"+"METTask_PfSumET_logx");  if (mePfSumET_logx  && mePfSumET_logx->getRootObject())  mePfSumET_logx->Fill(log10(pfSumET));
00977 
00978     //mePfMETIonFeedbck = _dbe->get(DirName+"/"+"METTask_PfMETIonFeedbck");  if (mePfMETIonFeedbck && mePfMETIonFeedbck->getRootObject()) mePfMETIonFeedbck->Fill(pfMET);
00979     //mePfMETHPDNoise   = _dbe->get(DirName+"/"+"METTask_PfMETHPDNoise");    if (mePfMETHPDNoise   && mePfMETHPDNoise->getRootObject())   mePfMETHPDNoise->Fill(pfMET);
00980     //mePfMETRBXNoise   = _dbe->get(DirName+"/"+"METTask_PfMETRBXNoise");    if (mePfMETRBXNoise   && mePfMETRBXNoise->getRootObject())   mePfMETRBXNoise->Fill(pfMET);
00981     
00982     mePfNeutralEMFraction = _dbe->get(DirName+"/"+"METTask_PfNeutralEMFraction"); 
00983     if (mePfNeutralEMFraction   && mePfNeutralEMFraction->getRootObject()) mePfNeutralEMFraction->Fill(pfNeutralEMFraction);
00984     mePfNeutralHadFraction = _dbe->get(DirName+"/"+"METTask_PfNeutralHadFraction"); 
00985     if (mePfNeutralHadFraction   && mePfNeutralHadFraction->getRootObject()) mePfNeutralHadFraction->Fill(pfNeutralHadFraction);
00986     mePfChargedEMFraction = _dbe->get(DirName+"/"+"METTask_PfChargedEMFraction"); 
00987     if (mePfChargedEMFraction   && mePfChargedEMFraction->getRootObject()) mePfChargedEMFraction->Fill(pfChargedEMFraction);
00988     mePfChargedHadFraction = _dbe->get(DirName+"/"+"METTask_PfChargedHadFraction"); 
00989     if (mePfChargedHadFraction   && mePfChargedHadFraction->getRootObject()) mePfChargedHadFraction->Fill(pfChargedHadFraction);
00990     mePfMuonFraction = _dbe->get(DirName+"/"+"METTask_PfMuonFraction"); 
00991     if (mePfMuonFraction   && mePfMuonFraction->getRootObject()) mePfMuonFraction->Fill(pfMuonFraction);
00992     
00993     if (_allhist){
00994       if (bLumiSecPlot){
00995         mePfMExLS = _dbe->get(DirName+"/"+"METTask_PfMExLS"); if (mePfMExLS && mePfMExLS->getRootObject()) mePfMExLS->Fill(pfMEx,myLuminosityBlock);
00996         mePfMEyLS = _dbe->get(DirName+"/"+"METTask_PfMEyLS"); if (mePfMEyLS && mePfMEyLS->getRootObject()) mePfMEyLS->Fill(pfMEy,myLuminosityBlock);
00997       }
00998     } // _allhist
00999   } // et threshold cut
01000 }
01001 
01002 // ***********************************************************
01003 bool PFMETAnalyzer::selectHighPtJetEvent(const edm::Event& iEvent){
01004 
01005   bool return_value=false;
01006 
01007   edm::Handle<reco::PFJetCollection> pfJets;
01008   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
01009   if (!pfJets.isValid()) {
01010     LogDebug("") << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
01011     if (_verbose) std::cout << "PFMETAnalyzer: Could not find pfjet product" << std::endl;
01012   }
01013 
01014   for (reco::PFJetCollection::const_iterator pf = pfJets->begin(); 
01015        pf!=pfJets->end(); ++pf){
01016     if (pf->pt()>_highPtPFJetThreshold){
01017       return_value=true;
01018     }
01019   }
01020   
01021   return return_value;
01022 }
01023 
01024 // // ***********************************************************
01025 bool PFMETAnalyzer::selectLowPtJetEvent(const edm::Event& iEvent){
01026 
01027   bool return_value=false;
01028 
01029   edm::Handle<reco::PFJetCollection> pfJets;
01030   iEvent.getByLabel(thePfJetCollectionLabel, pfJets);
01031   if (!pfJets.isValid()) {
01032     LogDebug("") << "PFMETAnalyzer: Could not find jet product" << std::endl;
01033     if (_verbose) std::cout << "PFMETAnalyzer: Could not find jet product" << std::endl;
01034   }
01035 
01036   for (reco::PFJetCollection::const_iterator cal = pfJets->begin(); 
01037        cal!=pfJets->end(); ++cal){
01038     if (cal->pt()>_lowPtPFJetThreshold){
01039       return_value=true;
01040     }
01041   }
01042 
01043   return return_value;
01044 
01045 }
01046 
01047 // ***********************************************************
01048 bool PFMETAnalyzer::selectWElectronEvent(const edm::Event& iEvent){
01049 
01050   bool return_value=true;
01051 
01052   /*
01053     W-electron event selection comes here
01054   */
01055 
01056   return return_value;
01057 
01058 }
01059 
01060 // ***********************************************************
01061 bool PFMETAnalyzer::selectWMuonEvent(const edm::Event& iEvent){
01062 
01063   bool return_value=true;
01064 
01065   /*
01066     W-muon event selection comes here
01067   */
01068 
01069   return return_value;
01070 
01071 }
01072
CMSSW_4_4_3_patch1/src/DQMOffline/JetMET/src/PFMETAnalyzer.cc
