CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/DQMOffline/JetMET/src/TcMETAnalyzer.cc

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2012/05/20 13:12:05 $
00005  *  $Revision: 1.13 $
00006  *  \author A.Apresyan - Caltech
00007  */
00008 
00009 #include "DQMOffline/JetMET/interface/TcMETAnalyzer.h"
00010 #include "DataFormats/Common/interface/Handle.h"
00011 #include "FWCore/Common/interface/TriggerNames.h"
00012 
00013 #include "DataFormats/Math/interface/LorentzVector.h"
00014 
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 
00017 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00018 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
00019 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00020 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00021 
00022 #include "DataFormats/Math/interface/LorentzVector.h"
00023 
00024 #include <string>
00025 using namespace edm;
00026 using namespace reco;
00027 using namespace math;
00028 
00029 // ***********************************************************
00030 TcMETAnalyzer::TcMETAnalyzer(const edm::ParameterSet& pSet) {
00031 
00032   parameters = pSet;
00033 
00034 }
00035 
00036 // ***********************************************************
00037 TcMETAnalyzer::~TcMETAnalyzer() { }
00038 
00039 void TcMETAnalyzer::beginJob(DQMStore * dbe) {
00040 
00041   evtCounter = 0;
00042   metname = "tcMETAnalyzer";
00043 
00044   // trigger information
00045   HLTPathsJetMBByName_ = parameters.getParameter<std::vector<std::string > >("HLTPathsJetMB");
00046 
00047   _hlt_HighPtJet = parameters.getParameter<std::string>("HLT_HighPtJet");
00048   _hlt_LowPtJet  = parameters.getParameter<std::string>("HLT_LowPtJet");
00049   _hlt_HighMET   = parameters.getParameter<std::string>("HLT_HighMET");
00050   //  _hlt_LowMET    = parameters.getParameter<std::string>("HLT_LowMET");
00051   _hlt_Ele       = parameters.getParameter<std::string>("HLT_Ele");
00052   _hlt_Muon      = parameters.getParameter<std::string>("HLT_Muon");
00053 
00054   // TcMET information
00055   theTcMETCollectionLabel       = parameters.getParameter<edm::InputTag>("TcMETCollectionLabel");
00056   _source                       = parameters.getParameter<std::string>("Source");
00057 
00058   // Other data collections
00059   HcalNoiseRBXCollectionTag   = parameters.getParameter<edm::InputTag>("HcalNoiseRBXCollection");
00060   theJetCollectionLabel       = parameters.getParameter<edm::InputTag>("JetCollectionLabel");
00061   HBHENoiseFilterResultTag    = parameters.getParameter<edm::InputTag>("HBHENoiseFilterResultLabel");
00062 
00063   // misc
00064   _verbose     = parameters.getParameter<int>("verbose");
00065   _etThreshold = parameters.getParameter<double>("etThreshold"); // MET threshold
00066   _allhist     = parameters.getParameter<bool>("allHist");       // Full set of monitoring histograms
00067   _allSelection= parameters.getParameter<bool>("allSelection");  // Plot with all sets of event selection
00068 
00069   _highPtTcJetThreshold = parameters.getParameter<double>("HighPtTcJetThreshold"); // High Pt Jet threshold
00070   _lowPtTcJetThreshold = parameters.getParameter<double>("LowPtTcJetThreshold");   // Low Pt Jet threshold
00071   _highTcMETThreshold = parameters.getParameter<double>("HighTcMETThreshold");     // High MET threshold
00072   _lowTcMETThreshold = parameters.getParameter<double>("LowTcMETThreshold");       // Low MET threshold
00073 
00074   //
00075   jetID = new reco::helper::JetIDHelper(parameters.getParameter<ParameterSet>("JetIDParams"));
00076 
00077   // DQStore stuff
00078   LogTrace(metname)<<"[TcMETAnalyzer] Parameters initialization";
00079   std::string DirName = "JetMET/MET/"+_source;
00080   dbe->setCurrentFolder(DirName);
00081 
00082   metME = dbe->book1D("metReco", "metReco", 4, 1, 5);
00083   metME->setBinLabel(2,"TcMET",1);
00084 
00085   _dbe = dbe;
00086 
00087   _FolderNames.push_back("All");
00088   _FolderNames.push_back("Cleanup");
00089   _FolderNames.push_back("HcalNoiseFilter");
00090   _FolderNames.push_back("JetID");
00091   _FolderNames.push_back("JetIDTight");
00092 
00093   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); 
00094        ic != _FolderNames.end(); ic++){
00095     if (*ic=="All")                  bookMESet(DirName+"/"+*ic);
00096     if (*ic=="Cleanup")              bookMESet(DirName+"/"+*ic);
00097     if (_allSelection){
00098     if (*ic=="HcalNoiseFilter")      bookMESet(DirName+"/"+*ic);
00099     if (*ic=="JetID")                bookMESet(DirName+"/"+*ic);
00100     if (*ic=="JetIDTight")           bookMESet(DirName+"/"+*ic);
00101     }
00102   }
00103 }
00104 
00105 // ***********************************************************
00106 void TcMETAnalyzer::endJob() {
00107 
00108   delete jetID;
00109 
00110 }
00111 
00112 // ***********************************************************
00113 void TcMETAnalyzer::bookMESet(std::string DirName)
00114 {
00115 
00116   bool bLumiSecPlot=false;
00117   if (DirName.find("All")!=std::string::npos) bLumiSecPlot=true;
00118 
00119   bookMonitorElement(DirName,bLumiSecPlot);
00120 
00121   if (_hlt_HighPtJet.size()){
00122     bookMonitorElement(DirName+"/"+"HighPtJet",false);
00123     meTriggerName_HighPtJet = _dbe->bookString("triggerName_HighPtJet", _hlt_HighPtJet);
00124   }  
00125 
00126   if (_hlt_LowPtJet.size()){
00127     bookMonitorElement(DirName+"/"+"LowPtJet",false);
00128     meTriggerName_LowPtJet = _dbe->bookString("triggerName_LowPtJet", _hlt_LowPtJet);
00129   }
00130 
00131   if (_hlt_HighMET.size()){
00132     bookMonitorElement(DirName+"/"+"HighMET",false);
00133     meTriggerName_HighMET = _dbe->bookString("triggerName_HighMET", _hlt_HighMET);
00134   }
00135 
00136   //  if (_hlt_LowMET.size()){
00137   //    bookMonitorElement(DirName+"/"+"LowMET",false);
00138   //    meTriggerName_LowMET = _dbe->bookString("triggerName_LowMET", _hlt_LowMET);
00139   //  }
00140 
00141   if (_hlt_Ele.size()){
00142     bookMonitorElement(DirName+"/"+"Ele",false);
00143     meTriggerName_Ele = _dbe->bookString("triggerName_Ele", _hlt_Ele);
00144   }
00145 
00146   if (_hlt_Muon.size()){
00147     bookMonitorElement(DirName+"/"+"Muon",false);
00148     meTriggerName_Muon = _dbe->bookString("triggerName_Muon", _hlt_Muon);
00149   }
00150 
00151 }
00152 
00153 // ***********************************************************
00154 void TcMETAnalyzer::bookMonitorElement(std::string DirName, bool bLumiSecPlot=false)
00155 {
00156 
00157   if (_verbose) std::cout << "booMonitorElement " << DirName << std::endl;
00158   _dbe->setCurrentFolder(DirName);
00159  
00160   meTcMEx    = _dbe->book1D("METTask_TcMEx",    "METTask_TcMEx",    200, -500,  500);
00161   meTcMEy    = _dbe->book1D("METTask_TcMEy",    "METTask_TcMEy",    200, -500,  500);
00162   meTcEz     = _dbe->book1D("METTask_TcEz",     "METTask_TcEz",     200, -500,  500);
00163   meTcMETSig = _dbe->book1D("METTask_TcMETSig", "METTask_TcMETSig",  51,    0,   51);
00164   meTcMET    = _dbe->book1D("METTask_TcMET",    "METTask_TcMET",    200,    0, 1000);
00165   meTcMETPhi = _dbe->book1D("METTask_TcMETPhi", "METTask_TcMETPhi",  60, -3.2,  3.2);
00166   meTcSumET  = _dbe->book1D("METTask_TcSumET",  "METTask_TcSumET",  400,    0, 4000);
00167 
00168   meTcNeutralEMFraction  = _dbe->book1D("METTask_TcNeutralEMFraction", "METTask_TcNeutralEMFraction" ,50,0.,1.);
00169   meTcNeutralHadFraction = _dbe->book1D("METTask_TcNeutralHadFraction","METTask_TcNeutralHadFraction",50,0.,1.);
00170   meTcChargedEMFraction  = _dbe->book1D("METTask_TcChargedEMFraction", "METTask_TcChargedEMFraction" ,50,0.,1.);
00171   meTcChargedHadFraction = _dbe->book1D("METTask_TcChargedHadFraction","METTask_TcChargedHadFraction",50,0.,1.);
00172   meTcMuonFraction       = _dbe->book1D("METTask_TcMuonFraction",      "METTask_TcMuonFraction"      ,50,0.,1.);
00173 
00174   meTcMETIonFeedbck      = _dbe->book1D("METTask_TcMETIonFeedbck", "METTask_TcMETIonFeedbck" ,500,0,1000);
00175   meTcMETHPDNoise        = _dbe->book1D("METTask_TcMETHPDNoise",   "METTask_TcMETHPDNoise"   ,500,0,1000);
00176   meTcMETRBXNoise        = _dbe->book1D("METTask_TcMETRBXNoise",   "METTask_TcMETRBXNoise"   ,500,0,1000);
00177 
00178   if (_allhist){
00179     if (bLumiSecPlot){
00180       meTcMExLS              = _dbe->book2D("METTask_TcMEx_LS","METTask_TcMEx_LS",200,-200,200,50,0.,500.);
00181       meTcMEyLS              = _dbe->book2D("METTask_TcMEy_LS","METTask_TcMEy_LS",200,-200,200,50,0.,500.);
00182     }
00183   }
00184 }
00185 
00186 // ***********************************************************
00187 void TcMETAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00188 {
00189 
00190 }
00191 
00192 // ***********************************************************
00193 void TcMETAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup, DQMStore * dbe)
00194 {
00195   
00196   //
00197   //--- Check the time length of the Run from the lumi section plots
00198 
00199   std::string dirName = "JetMET/MET/"+_source+"/";
00200   _dbe->setCurrentFolder(dirName);
00201 
00202   TH1F* tlumisec;
00203 
00204   MonitorElement *meLumiSec = _dbe->get("aaa");
00205   meLumiSec = _dbe->get("JetMET/lumisec");
00206 
00207   int totlsec=0;
00208   double totltime=0.;
00209   if ( meLumiSec->getRootObject() ) {
00210     tlumisec = meLumiSec->getTH1F();
00211     for (int i=0; i<500; i++){
00212       if (tlumisec->GetBinContent(i+1)) totlsec++;
00213     }
00214     totltime = double(totlsec*90); // one lumi sec ~ 90 (sec)
00215   }
00216 
00217   if (totltime==0.) totltime=1.; 
00218 
00219   //
00220   //--- Make the integrated plots with rate (Hz)
00221 
00222   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); ic != _FolderNames.end(); ic++)
00223     {
00224 
00225       std::string DirName;
00226       DirName = dirName+*ic;
00227 
00228       makeRatePlot(DirName,totltime);
00229       if (_hlt_HighPtJet.size()) makeRatePlot(DirName+"/"+_hlt_HighPtJet,totltime);
00230       if (_hlt_LowPtJet.size())  makeRatePlot(DirName+"/"+_hlt_LowPtJet,totltime);
00231       if (_hlt_HighMET.size())   makeRatePlot(DirName+"/"+_hlt_HighMET,totltime);
00232       //      if (_hlt_LowMET.size())    makeRatePlot(DirName+"/"+_hlt_LowMET,totltime);
00233       if (_hlt_Ele.size())       makeRatePlot(DirName+"/"+_hlt_Ele,totltime);
00234       if (_hlt_Muon.size())      makeRatePlot(DirName+"/"+_hlt_Muon,totltime);
00235 
00236     }
00237 }
00238 
00239 
00240 // ***********************************************************
00241 void TcMETAnalyzer::makeRatePlot(std::string DirName, double totltime)
00242 {
00243 
00244   _dbe->setCurrentFolder(DirName);
00245   MonitorElement *meTcMET = _dbe->get(DirName+"/"+"METTask_TcMET");
00246 
00247   TH1F* tTcMET;
00248   TH1F* tTcMETRate;
00249 
00250   if ( meTcMET )
00251     if ( meTcMET->getRootObject() ) {
00252       tTcMET     = meTcMET->getTH1F();
00253       
00254       // Integral plot & convert number of events to rate (hz)
00255       tTcMETRate = (TH1F*) tTcMET->Clone("METTask_TcMETRate");
00256       for (int i = tTcMETRate->GetNbinsX()-1; i>=0; i--){
00257         tTcMETRate->SetBinContent(i+1,tTcMETRate->GetBinContent(i+2)+tTcMET->GetBinContent(i+1));
00258       }
00259       for (int i = 0; i<tTcMETRate->GetNbinsX(); i++){
00260         tTcMETRate->SetBinContent(i+1,tTcMETRate->GetBinContent(i+1)/double(totltime));
00261       }      
00262 
00263       meTcMETRate      = _dbe->book1D("METTask_TcMETRate",tTcMETRate);
00264       
00265     }
00266 
00267 }
00268 
00269 // ***********************************************************
00270 void TcMETAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
00271                             const edm::TriggerResults& triggerResults) {
00272 
00273   if (_verbose) std::cout << "TcMETAnalyzer analyze" << std::endl;
00274 
00275   LogTrace(metname)<<"[TcMETAnalyzer] Analyze TcMET";
00276 
00277   metME->Fill(2);
00278 
00279   // ==========================================================  
00280   // Trigger information 
00281   //
00282   _trig_JetMB=0;
00283   _trig_HighPtJet=0;
00284   _trig_LowPtJet=0;
00285   _trig_HighMET=0;
00286   //  _trig_LowMET=0;
00287 
00288   if (&triggerResults) {   
00289 
00291 
00292     //
00293     //
00294     // Check how many HLT triggers are in triggerResults 
00295     int ntrigs = triggerResults.size();
00296     if (_verbose) std::cout << "ntrigs=" << ntrigs << std::endl;
00297     
00298     //
00299     //
00300     // If index=ntrigs, this HLT trigger doesn't exist in the HLT table for this data.
00301     const edm::TriggerNames & triggerNames = iEvent.triggerNames(triggerResults);
00302     
00303     //
00304     //
00305     // count number of requested Jet or MB HLT paths which have fired
00306     for (unsigned int i=0; i!=HLTPathsJetMBByName_.size(); i++) {
00307       unsigned int triggerIndex = triggerNames.triggerIndex(HLTPathsJetMBByName_[i]);
00308       if (triggerIndex<triggerResults.size()) {
00309         if (triggerResults.accept(triggerIndex)) {
00310           _trig_JetMB++;
00311         }
00312       }
00313     }
00314     // for empty input vectors (n==0), take all HLT triggers!
00315     if (HLTPathsJetMBByName_.size()==0) _trig_JetMB=triggerResults.size()-1;
00316 
00317     //
00318     if (_verbose) std::cout << "triggerNames size" << " " << triggerNames.size() << std::endl;
00319     if (_verbose) std::cout << _hlt_HighPtJet << " " << triggerNames.triggerIndex(_hlt_HighPtJet) << std::endl;
00320     if (_verbose) std::cout << _hlt_LowPtJet  << " " << triggerNames.triggerIndex(_hlt_LowPtJet)  << std::endl;
00321     if (_verbose) std::cout << _hlt_HighMET   << " " << triggerNames.triggerIndex(_hlt_HighMET)   << std::endl;
00322     //    if (_verbose) std::cout << _hlt_LowMET    << " " << triggerNames.triggerIndex(_hlt_LowMET)    << std::endl;
00323     if (_verbose) std::cout << _hlt_Ele       << " " << triggerNames.triggerIndex(_hlt_Ele)       << std::endl;
00324     if (_verbose) std::cout << _hlt_Muon      << " " << triggerNames.triggerIndex(_hlt_Muon)      << std::endl;
00325 
00326     if (triggerNames.triggerIndex(_hlt_HighPtJet) != triggerNames.size() &&
00327         triggerResults.accept(triggerNames.triggerIndex(_hlt_HighPtJet))) _trig_HighPtJet=1;
00328 
00329     if (triggerNames.triggerIndex(_hlt_LowPtJet)  != triggerNames.size() &&
00330         triggerResults.accept(triggerNames.triggerIndex(_hlt_LowPtJet)))  _trig_LowPtJet=1;
00331 
00332     if (triggerNames.triggerIndex(_hlt_HighMET)   != triggerNames.size() &&
00333         triggerResults.accept(triggerNames.triggerIndex(_hlt_HighMET)))   _trig_HighMET=1;
00334 
00335     //    if (triggerNames.triggerIndex(_hlt_LowMET)    != triggerNames.size() &&
00336     //        triggerResults.accept(triggerNames.triggerIndex(_hlt_LowMET)))    _trig_LowMET=1;
00337 
00338     if (triggerNames.triggerIndex(_hlt_Ele)       != triggerNames.size() &&
00339         triggerResults.accept(triggerNames.triggerIndex(_hlt_Ele)))       _trig_Ele=1;
00340 
00341     if (triggerNames.triggerIndex(_hlt_Muon)      != triggerNames.size() &&
00342         triggerResults.accept(triggerNames.triggerIndex(_hlt_Muon)))      _trig_Muon=1;
00343     
00344   } else {
00345 
00346     edm::LogInfo("TcMetAnalyzer") << "TriggerResults::HLT not found, "
00347       "automatically select events"; 
00348 
00349     // TriggerResults object not found. Look at all events.    
00350     _trig_JetMB=1;
00351   }
00352 
00353   // ==========================================================
00354   // TcMET information
00355   
00356   // **** Get the MET container  
00357   edm::Handle<reco::METCollection> tcmetcoll;
00358   iEvent.getByLabel(theTcMETCollectionLabel, tcmetcoll);
00359   
00360   if(!tcmetcoll.isValid()) return;
00361 
00362   const METCollection *tcmetcol = tcmetcoll.product();
00363   const MET *tcmet;
00364   tcmet = &(tcmetcol->front());
00365     
00366   LogTrace(metname)<<"[TcMETAnalyzer] Call to the TcMET analyzer";
00367 
00368   // ==========================================================
00369   //
00370   edm::Handle<HcalNoiseRBXCollection> HRBXCollection;
00371   iEvent.getByLabel(HcalNoiseRBXCollectionTag,HRBXCollection);
00372   if (!HRBXCollection.isValid()) {
00373     LogDebug("") << "TcMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00374     if (_verbose) std::cout << "TcMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00375   }
00376 
00377 
00378   edm::Handle<bool> HBHENoiseFilterResultHandle;
00379   iEvent.getByLabel(HBHENoiseFilterResultTag, HBHENoiseFilterResultHandle);
00380   bool HBHENoiseFilterResult = *HBHENoiseFilterResultHandle;
00381   if (!HBHENoiseFilterResultHandle.isValid()) {
00382     LogDebug("") << "TcMETAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
00383     if (_verbose) std::cout << "TcMETAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
00384   }
00385 
00386 
00387   edm::Handle<reco::CaloJetCollection> caloJets;
00388   iEvent.getByLabel(theJetCollectionLabel, caloJets);
00389   if (!caloJets.isValid()) {
00390     LogDebug("") << "TcMETAnalyzer: Could not find jet product" << std::endl;
00391     if (_verbose) std::cout << "TcMETAnalyzer: Could not find jet product" << std::endl;
00392   }
00393 
00394   // ==========================================================
00395   // TcMET sanity check
00396 
00397   //   if (_source=="TcMET") validateMET(*tcmet, tcCandidates);
00398   
00399   // ==========================================================
00400   // JetID 
00401 
00402   if (_verbose) std::cout << "JetID starts" << std::endl;
00403   
00404   //
00405   // --- Loose cuts, not Tc specific for now!
00406   //
00407   bool bJetID=true;
00408   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00409        cal!=caloJets->end(); ++cal){ 
00410     jetID->calculate(iEvent, *cal);
00411     if (_verbose) std::cout << jetID->n90Hits() << " " 
00412                             << jetID->restrictedEMF() << " "
00413                             << cal->pt() << std::endl;
00414     if (cal->pt()>10.){
00415       //
00416       // for all regions
00417       if (jetID->n90Hits()<2)  bJetID=false; 
00418       if (jetID->fHPD()>=0.98) bJetID=false; 
00419       //if (jetID->restrictedEMF()<0.01) bJetID=false; 
00420       //
00421       // for non-forward
00422       if (fabs(cal->eta())<2.55){
00423         if (cal->emEnergyFraction()<=0.01) bJetID=false; 
00424       }
00425       // for forward
00426       else {
00427         if (cal->emEnergyFraction()<=-0.9) bJetID=false; 
00428         if (cal->pt()>80.){
00429         if (cal->emEnergyFraction()>= 1.0) bJetID=false; 
00430         }
00431       } // forward vs non-forward
00432     }   // pt>10 GeV/c
00433   }     // calor-jets loop
00434  
00435   //
00436   // --- Tight cuts
00437   //
00438   bool bJetIDTight=true;
00439   bJetIDTight=bJetID;
00440   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00441        cal!=caloJets->end(); ++cal){
00442     jetID->calculate(iEvent, *cal);
00443     if (cal->pt()>25.){
00444       //
00445       // for all regions
00446       if (jetID->fHPD()>=0.95) bJetIDTight=false; 
00447       //
00448       // for 1.0<|eta|<1.75
00449       if (fabs(cal->eta())>=1.00 && fabs(cal->eta())<1.75){
00450         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00451       }
00452       //
00453       // for 1.75<|eta|<2.55
00454       else if (fabs(cal->eta())>=1.75 && fabs(cal->eta())<2.55){
00455         if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false; 
00456       }
00457       //
00458       // for 2.55<|eta|<3.25
00459       else if (fabs(cal->eta())>=2.55 && fabs(cal->eta())<3.25){
00460         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3) bJetIDTight=false; 
00461         if (cal->pt()>=50. && cal->pt()< 80. && cal->emEnergyFraction()<=-0.2) bJetIDTight=false; 
00462         if (cal->pt()>=80. && cal->pt()<340. && cal->emEnergyFraction()<=-0.1) bJetIDTight=false; 
00463         if (cal->pt()>=340.                  && cal->emEnergyFraction()<=-0.1 
00464                                              && cal->emEnergyFraction()>=0.95) bJetIDTight=false; 
00465       }
00466       //
00467       // for 3.25<|eta|
00468       else if (fabs(cal->eta())>=3.25){
00469         if (cal->pt()< 50.                   && cal->emEnergyFraction()<=-0.3
00470                                              && cal->emEnergyFraction()>=0.90) bJetIDTight=false; 
00471         if (cal->pt()>=50. && cal->pt()<130. && cal->emEnergyFraction()<=-0.2
00472                                              && cal->emEnergyFraction()>=0.80) bJetIDTight=false; 
00473         if (cal->pt()>=130.                  && cal->emEnergyFraction()<=-0.1 
00474                                              && cal->emEnergyFraction()>=0.70) bJetIDTight=false; 
00475       }
00476     }   // pt>10 GeV/c
00477   }     // calor-jets loop
00478   
00479   if (_verbose) std::cout << "JetID ends" << std::endl;
00480 
00481 
00482   // ==========================================================
00483   // HCAL Noise filter
00484   
00485   bool bHcalNoiseFilter = HBHENoiseFilterResult;
00486 
00487   // ==========================================================
00488   // Reconstructed MET Information - fill MonitorElements
00489   
00490   std::string DirName = "JetMET/MET/"+_source;
00491   
00492   for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); 
00493        ic != _FolderNames.end(); ic++){
00494     if (*ic=="All")                                   fillMESet(iEvent, DirName+"/"+*ic, *tcmet);
00495     if (*ic=="Cleanup" && bHcalNoiseFilter && bJetID) fillMESet(iEvent, DirName+"/"+*ic, *tcmet);
00496     if (_allSelection) {
00497     if (*ic=="HcalNoiseFilter"      && bHcalNoiseFilter )       fillMESet(iEvent, DirName+"/"+*ic, *tcmet);
00498     if (*ic=="JetID"      && bJetID)                            fillMESet(iEvent, DirName+"/"+*ic, *tcmet);
00499     if (*ic=="JetIDTight" && bJetIDTight)                       fillMESet(iEvent, DirName+"/"+*ic, *tcmet);
00500     }
00501   }
00502 }
00503 
00504 // ***********************************************************
00505 void TcMETAnalyzer::fillMESet(const edm::Event& iEvent, std::string DirName, 
00506                               const reco::MET& tcmet)
00507 {
00508 
00509   _dbe->setCurrentFolder(DirName);
00510 
00511   bool bLumiSecPlot=false;
00512   if (DirName.find("All")) bLumiSecPlot=true;
00513 
00514   if (_trig_JetMB) fillMonitorElement(iEvent,DirName,"",tcmet, bLumiSecPlot);
00515   if (_hlt_HighPtJet.size() && _trig_HighPtJet) fillMonitorElement(iEvent,DirName,"HighPtJet",tcmet,false);
00516   if (_hlt_LowPtJet.size() && _trig_LowPtJet) fillMonitorElement(iEvent,DirName,"LowPtJet",tcmet,false);
00517   if (_hlt_HighMET.size() && _trig_HighMET) fillMonitorElement(iEvent,DirName,"HighMET",tcmet,false);
00518   //  if (_hlt_LowMET.size() && _trig_LowMET) fillMonitorElement(iEvent,DirName,"LowMET",tcmet,false);
00519   if (_hlt_Ele.size() && _trig_Ele) fillMonitorElement(iEvent,DirName,"Ele",tcmet,false);
00520   if (_hlt_Muon.size() && _trig_Muon) fillMonitorElement(iEvent,DirName,"Muon",tcmet,false);
00521 }
00522 
00523 // ***********************************************************
00524 void TcMETAnalyzer::fillMonitorElement(const edm::Event& iEvent, std::string DirName, 
00525                                          std::string TriggerTypeName, 
00526                                          const reco::MET& tcmet, bool bLumiSecPlot)
00527 {
00528 
00529   if (TriggerTypeName=="HighPtJet") {
00530     if (!selectHighPtJetEvent(iEvent)) return;
00531   }
00532   else if (TriggerTypeName=="LowPtJet") {
00533     if (!selectLowPtJetEvent(iEvent)) return;
00534   }
00535   else if (TriggerTypeName=="HighMET") {
00536     if (tcmet.pt()<_highTcMETThreshold) return;
00537   }
00538   //  else if (TriggerTypeName=="LowMET") {
00539   //    if (tcmet.pt()<_lowTcMETThreshold) return;
00540   //  }
00541   else if (TriggerTypeName=="Ele") {
00542     if (!selectWElectronEvent(iEvent)) return;
00543   }
00544   else if (TriggerTypeName=="Muon") {
00545     if (!selectWMuonEvent(iEvent)) return;
00546   }
00547   
00548 // Reconstructed MET Information
00549   double tcSumET  = tcmet.sumEt();
00550   double tcMETSig = tcmet.mEtSig();
00551   double tcEz     = tcmet.e_longitudinal();
00552   double tcMET    = tcmet.pt();
00553   double tcMEx    = tcmet.px();
00554   double tcMEy    = tcmet.py();
00555   double tcMETPhi = tcmet.phi();
00556 
00557   //
00558   int myLuminosityBlock;
00559   //  myLuminosityBlock = (evtCounter++)/1000;
00560   myLuminosityBlock = iEvent.luminosityBlock();
00561   //
00562 
00563   if (TriggerTypeName!="") DirName = DirName +"/"+TriggerTypeName;
00564 
00565   if (_verbose) std::cout << "_etThreshold = " << _etThreshold << std::endl;
00566   if (tcMET>_etThreshold){
00567     
00568     meTcMEx    = _dbe->get(DirName+"/"+"METTask_TcMEx");    if (meTcMEx    && meTcMEx->getRootObject())    meTcMEx->Fill(tcMEx);
00569     meTcMEy    = _dbe->get(DirName+"/"+"METTask_TcMEy");    if (meTcMEy    && meTcMEy->getRootObject())    meTcMEy->Fill(tcMEy);
00570     meTcMET    = _dbe->get(DirName+"/"+"METTask_TcMET");    if (meTcMET    && meTcMET->getRootObject())    meTcMET->Fill(tcMET);
00571     meTcMETPhi = _dbe->get(DirName+"/"+"METTask_TcMETPhi"); if (meTcMETPhi && meTcMETPhi->getRootObject()) meTcMETPhi->Fill(tcMETPhi);
00572     meTcSumET  = _dbe->get(DirName+"/"+"METTask_TcSumET");  if (meTcSumET  && meTcSumET->getRootObject())  meTcSumET->Fill(tcSumET);
00573     meTcMETSig = _dbe->get(DirName+"/"+"METTask_TcMETSig"); if (meTcMETSig && meTcMETSig->getRootObject()) meTcMETSig->Fill(tcMETSig);
00574     meTcEz     = _dbe->get(DirName+"/"+"METTask_TcEz");     if (meTcEz     && meTcEz->getRootObject())     meTcEz->Fill(tcEz);
00575 
00576     meTcMETIonFeedbck = _dbe->get(DirName+"/"+"METTask_TcMETIonFeedbck");  if (meTcMETIonFeedbck && meTcMETIonFeedbck->getRootObject()) meTcMETIonFeedbck->Fill(tcMET);
00577     meTcMETHPDNoise   = _dbe->get(DirName+"/"+"METTask_TcMETHPDNoise");    if (meTcMETHPDNoise   && meTcMETHPDNoise->getRootObject())   meTcMETHPDNoise->Fill(tcMET);
00578     meTcMETRBXNoise   = _dbe->get(DirName+"/"+"METTask_TcMETRBXNoise");    if (meTcMETRBXNoise   && meTcMETRBXNoise->getRootObject())   meTcMETRBXNoise->Fill(tcMET);
00579         
00580     if (_allhist){
00581       if (bLumiSecPlot){
00582         meTcMExLS = _dbe->get(DirName+"/"+"METTask_TcMExLS"); if (meTcMExLS && meTcMExLS->getRootObject()) meTcMExLS->Fill(tcMEx,myLuminosityBlock);
00583         meTcMEyLS = _dbe->get(DirName+"/"+"METTask_TcMEyLS"); if (meTcMEyLS && meTcMEyLS->getRootObject()) meTcMEyLS->Fill(tcMEy,myLuminosityBlock);
00584       }
00585     } // _allhist
00586   } // et threshold cut
00587 }
00588 
00589 
00590 // ***********************************************************
00591 bool TcMETAnalyzer::selectHighPtJetEvent(const edm::Event& iEvent){
00592 
00593   bool return_value=false;
00594 
00595   edm::Handle<reco::CaloJetCollection> caloJets;
00596   iEvent.getByLabel(theJetCollectionLabel, caloJets);
00597   if (!caloJets.isValid()) {
00598     LogDebug("") << "TcMETAnalyzer: Could not find jet product" << std::endl;
00599     if (_verbose) std::cout << "TcMETAnalyzer: Could not find jet product" << std::endl;
00600   }
00601 
00602   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00603        cal!=caloJets->end(); ++cal){
00604     if (cal->pt()>_highPtTcJetThreshold){
00605       return_value=true;
00606     }
00607   }
00608   
00609   return return_value;
00610 }
00611 
00612 // // ***********************************************************
00613 bool TcMETAnalyzer::selectLowPtJetEvent(const edm::Event& iEvent){
00614 
00615   bool return_value=false;
00616 
00617   edm::Handle<reco::CaloJetCollection> caloJets;
00618   iEvent.getByLabel(theJetCollectionLabel, caloJets);
00619   if (!caloJets.isValid()) {
00620     LogDebug("") << "TcMETAnalyzer: Could not find jet product" << std::endl;
00621     if (_verbose) std::cout << "TcMETAnalyzer: Could not find jet product" << std::endl;
00622   }
00623 
00624   for (reco::CaloJetCollection::const_iterator cal = caloJets->begin(); 
00625        cal!=caloJets->end(); ++cal){
00626     if (cal->pt()>_lowPtTcJetThreshold){
00627       return_value=true;
00628     }
00629   }
00630 
00631   return return_value;
00632 
00633 }
00634 
00635 // ***********************************************************
00636 bool TcMETAnalyzer::selectWElectronEvent(const edm::Event& iEvent){
00637 
00638   bool return_value=false;
00639 
00640   /*
00641     W-electron event selection comes here
00642    */
00643 
00644   return return_value;
00645 
00646 }
00647 
00648 // ***********************************************************
00649 bool TcMETAnalyzer::selectWMuonEvent(const edm::Event& iEvent){
00650 
00651   bool return_value=false;
00652 
00653   /*
00654     W-muon event selection comes here
00655    */
00656 
00657   return return_value;
00658 
00659 }