00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "DQMOffline/JetMET/interface/CaloMETAnalyzer.h"
00011 #include "DataFormats/Common/interface/Handle.h"
00012
00013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00014 #include "FWCore/Common/interface/TriggerNames.h"
00015
00016 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00017 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
00018 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00019 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00020
00021 #include "DataFormats/Math/interface/LorentzVector.h"
00022
00023 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
00024 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
00025 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00026
00027 #include "TLorentzVector.h"
00028
00029 #include <string>
00030 using namespace edm;
00031
00032
00033 CaloMETAnalyzer::CaloMETAnalyzer(const edm::ParameterSet& pSet) {
00034
00035 parameters = pSet;
00036
00037 edm::ParameterSet highptjetparms = parameters.getParameter<edm::ParameterSet>("highPtJetTrigger");
00038 edm::ParameterSet lowptjetparms = parameters.getParameter<edm::ParameterSet>("lowPtJetTrigger" );
00039 edm::ParameterSet minbiasparms = parameters.getParameter<edm::ParameterSet>("minBiasTrigger" );
00040 edm::ParameterSet highmetparms = parameters.getParameter<edm::ParameterSet>("highMETTrigger" );
00041 edm::ParameterSet lowmetparms = parameters.getParameter<edm::ParameterSet>("lowMETTrigger" );
00042 edm::ParameterSet eleparms = parameters.getParameter<edm::ParameterSet>("eleTrigger" );
00043 edm::ParameterSet muonparms = parameters.getParameter<edm::ParameterSet>("muonTrigger" );
00044
00045 _hlt_HighPtJet = highptjetparms.getParameter<std::string>("hltDBKey");
00046 _hlt_LowPtJet = lowptjetparms .getParameter<std::string>("hltDBKey");
00047 _hlt_MinBias = minbiasparms .getParameter<std::string>("hltDBKey");
00048 _hlt_HighMET = highmetparms .getParameter<std::string>("hltDBKey");
00049 _hlt_LowMET = lowmetparms .getParameter<std::string>("hltDBKey");
00050 _hlt_Ele = eleparms .getParameter<std::string>("hltDBKey");
00051 _hlt_Muon = muonparms .getParameter<std::string>("hltDBKey");
00052
00053
00054
00055 _HighPtJetEventFlag = new GenericTriggerEventFlag( highptjetparms );
00056 _LowPtJetEventFlag = new GenericTriggerEventFlag( lowptjetparms );
00057 _MinBiasEventFlag = new GenericTriggerEventFlag( minbiasparms );
00058 _HighMETEventFlag = new GenericTriggerEventFlag( highmetparms );
00059 _LowMETEventFlag = new GenericTriggerEventFlag( lowmetparms );
00060 _EleEventFlag = new GenericTriggerEventFlag( eleparms );
00061 _MuonEventFlag = new GenericTriggerEventFlag( muonparms );
00062
00063 }
00064
00065
00066 CaloMETAnalyzer::~CaloMETAnalyzer() {
00067
00068 delete _HighPtJetEventFlag;
00069 delete _LowPtJetEventFlag;
00070 delete _MinBiasEventFlag;
00071 delete _HighMETEventFlag;
00072 delete _LowMETEventFlag;
00073 delete _EleEventFlag;
00074 delete _MuonEventFlag;
00075
00076 }
00077
00078
00079 void CaloMETAnalyzer::beginJob(DQMStore * dbe) {
00080
00081 evtCounter = 0;
00082 metname = "caloMETAnalyzer";
00083
00084
00085 HLTPathsJetMBByName_ = parameters.getParameter<std::vector<std::string > >("HLTPathsJetMB");
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095 theCleaningParameters = parameters.getParameter<ParameterSet>("CleaningParameters"),
00096
00097
00098 gtTag = theCleaningParameters.getParameter<edm::InputTag>("gtLabel");
00099 _techTrigsAND = theCleaningParameters.getParameter<std::vector<unsigned > >("techTrigsAND");
00100 _techTrigsOR = theCleaningParameters.getParameter<std::vector<unsigned > >("techTrigsOR");
00101 _techTrigsNOT = theCleaningParameters.getParameter<std::vector<unsigned > >("techTrigsNOT");
00102
00103 _doHLTPhysicsOn = theCleaningParameters.getParameter<bool>("doHLTPhysicsOn");
00104 _hlt_PhysDec = theCleaningParameters.getParameter<std::string>("HLT_PhysDec");
00105
00106 _tightBHFiltering = theCleaningParameters.getParameter<bool>("tightBHFiltering");
00107 _tightJetIDFiltering = theCleaningParameters.getParameter<int>("tightJetIDFiltering");
00108 _tightHcalFiltering = theCleaningParameters.getParameter<bool>("tightHcalFiltering");
00109
00110
00111
00112
00113 DCSFilter = new JetMETDQMDCSFilter(parameters.getParameter<ParameterSet>("DCSFilter"));
00114
00115
00116 _doPVCheck = theCleaningParameters.getParameter<bool>("doPrimaryVertexCheck");
00117 vertexTag = theCleaningParameters.getParameter<edm::InputTag>("vertexLabel");
00118
00119 if (_doPVCheck) {
00120 _nvtx_min = theCleaningParameters.getParameter<int>("nvtx_min");
00121 _nvtxtrks_min = theCleaningParameters.getParameter<int>("nvtxtrks_min");
00122 _vtxndof_min = theCleaningParameters.getParameter<int>("vtxndof_min");
00123 _vtxchi2_max = theCleaningParameters.getParameter<double>("vtxchi2_max");
00124 _vtxz_max = theCleaningParameters.getParameter<double>("vtxz_max");
00125 }
00126
00127
00128
00129 theCaloMETCollectionLabel = parameters.getParameter<edm::InputTag>("METCollectionLabel");
00130 _source = parameters.getParameter<std::string>("Source");
00131
00132 if (theCaloMETCollectionLabel.label() == "corMetGlobalMuons" ) {
00133 inputBeamSpotLabel = parameters.getParameter<edm::InputTag>("InputBeamSpotLabel");
00134 }
00135
00136
00137 theCaloTowersLabel = parameters.getParameter<edm::InputTag>("CaloTowersLabel");
00138 theJetCollectionLabel = parameters.getParameter<edm::InputTag>("JetCollectionLabel");
00139 HcalNoiseRBXCollectionTag = parameters.getParameter<edm::InputTag>("HcalNoiseRBXCollection");
00140 HcalNoiseSummaryTag = parameters.getParameter<edm::InputTag>("HcalNoiseSummary");
00141 BeamHaloSummaryTag = parameters.getParameter<edm::InputTag>("BeamHaloSummaryLabel");
00142
00143
00144 _verbose = parameters.getParameter<int>("verbose");
00145 _print = parameters.getParameter<int>("printOut");
00146 _etThreshold = parameters.getParameter<double>("etThreshold");
00147 _allhist = parameters.getParameter<bool>("allHist");
00148 _allSelection= parameters.getParameter<bool>("allSelection");
00149 _cleanupSelection= parameters.getParameter<bool>("cleanupSelection");
00150
00151 _highPtJetThreshold = parameters.getParameter<double>("HighPtJetThreshold");
00152 _lowPtJetThreshold = parameters.getParameter<double>("LowPtJetThreshold");
00153 _highMETThreshold = parameters.getParameter<double>("HighMETThreshold");
00154 _lowMETThreshold = parameters.getParameter<double>("LowMETThreshold");
00155
00156
00157 jetID = new reco::helper::JetIDHelper(parameters.getParameter<ParameterSet>("JetIDParams"));
00158
00159
00160 LogTrace(metname)<<"[CaloMETAnalyzer] Parameters initialization";
00161 std::string DirName = "JetMET/MET/"+_source;
00162 dbe->setCurrentFolder(DirName);
00163
00164 hmetME = dbe->book1D("metReco", "metReco", 4, 1, 5);
00165 hmetME->setBinLabel(1,"CaloMET",1);
00166
00167 _dbe = dbe;
00168
00169 _FolderNames.push_back("All");
00170 _FolderNames.push_back("BasicCleanup");
00171 _FolderNames.push_back("ExtraCleanup");
00172 _FolderNames.push_back("HcalNoiseFilter");
00173 _FolderNames.push_back("HcalNoiseFilterTight");
00174 _FolderNames.push_back("JetIDMinimal");
00175 _FolderNames.push_back("JetIDLoose");
00176 _FolderNames.push_back("JetIDTight");
00177 _FolderNames.push_back("BeamHaloIDTightPass");
00178 _FolderNames.push_back("BeamHaloIDLoosePass");
00179 _FolderNames.push_back("Triggers");
00180 _FolderNames.push_back("PV");
00181
00182 for (std::vector<std::string>::const_iterator ic = _FolderNames.begin();
00183 ic != _FolderNames.end(); ic++){
00184 if (*ic=="All") bookMESet(DirName+"/"+*ic);
00185 if (_cleanupSelection){
00186 if (*ic=="BasicCleanup") bookMESet(DirName+"/"+*ic);
00187 if (*ic=="ExtraCleanup") bookMESet(DirName+"/"+*ic);
00188 }
00189 if (_allSelection){
00190 if (*ic=="HcalNoiseFilter") bookMESet(DirName+"/"+*ic);
00191 if (*ic=="HcalNoiseFilterTight") bookMESet(DirName+"/"+*ic);
00192 if (*ic=="JetIDMinimal") bookMESet(DirName+"/"+*ic);
00193 if (*ic=="JetIDLoose") bookMESet(DirName+"/"+*ic);
00194 if (*ic=="JetIDTight") bookMESet(DirName+"/"+*ic);
00195 if (*ic=="BeamHaloIDTightPass") bookMESet(DirName+"/"+*ic);
00196 if (*ic=="BeamHaloIDLoosePass") bookMESet(DirName+"/"+*ic);
00197 if (*ic=="Triggers") bookMESet(DirName+"/"+*ic);
00198 if (*ic=="PV") bookMESet(DirName+"/"+*ic);
00199 }
00200 }
00201
00202 }
00203
00204
00205 void CaloMETAnalyzer::endJob() {
00206
00207 delete jetID;
00208 delete DCSFilter;
00209
00210 }
00211
00212
00213 void CaloMETAnalyzer::bookMESet(std::string DirName)
00214 {
00215
00216 bool bLumiSecPlot=false;
00217 if (DirName.find("All")!=std::string::npos) bLumiSecPlot=true;
00218
00219 bookMonitorElement(DirName,bLumiSecPlot);
00220
00221 if ( _HighPtJetEventFlag->on() ) {
00222 bookMonitorElement(DirName+"/"+"HighPtJet",false);
00223 hTriggerName_HighPtJet = _dbe->bookString("triggerName_HighPtJet", _hlt_HighPtJet);
00224 }
00225
00226 if ( _LowPtJetEventFlag->on() ) {
00227 bookMonitorElement(DirName+"/"+"LowPtJet",false);
00228 hTriggerName_LowPtJet = _dbe->bookString("triggerName_LowPtJet", _hlt_LowPtJet);
00229 }
00230
00231 if ( _MinBiasEventFlag->on() ) {
00232 bookMonitorElement(DirName+"/"+"MinBias",false);
00233 hTriggerName_MinBias = _dbe->bookString("triggerName_MinBias", _hlt_MinBias);
00234 }
00235
00236 if ( _HighMETEventFlag->on() ) {
00237 bookMonitorElement(DirName+"/"+"HighMET",false);
00238 hTriggerName_HighMET = _dbe->bookString("triggerName_HighMET", _hlt_HighMET);
00239 }
00240
00241 if ( _LowMETEventFlag->on() ) {
00242 bookMonitorElement(DirName+"/"+"LowMET",false);
00243 hTriggerName_LowMET = _dbe->bookString("triggerName_LowMET", _hlt_LowMET);
00244 }
00245
00246 if ( _EleEventFlag->on() ) {
00247 bookMonitorElement(DirName+"/"+"Ele",false);
00248 hTriggerName_Ele = _dbe->bookString("triggerName_Ele", _hlt_Ele);
00249 }
00250
00251 if ( _MuonEventFlag->on() ) {
00252 bookMonitorElement(DirName+"/"+"Muon",false);
00253 hTriggerName_Muon = _dbe->bookString("triggerName_Muon", _hlt_Muon);
00254 }
00255 }
00256
00257
00258 void CaloMETAnalyzer::bookMonitorElement(std::string DirName, bool bLumiSecPlot=false)
00259 {
00260
00261 if (_verbose) std::cout << "bookMonitorElement " << DirName << std::endl;
00262 _dbe->setCurrentFolder(DirName);
00263
00264 hNevents = _dbe->book1D("METTask_Nevents", "METTask_Nevents" ,1,0,1);
00265 hCaloMEx = _dbe->book1D("METTask_CaloMEx", "METTask_CaloMEx" ,500,-500,500);
00266 hCaloMEx->setAxisTitle("MEx [GeV]",1);
00267 hCaloMEy = _dbe->book1D("METTask_CaloMEy", "METTask_CaloMEy" ,500,-500,500);
00268 hCaloMEy->setAxisTitle("MEy [GeV]",1);
00269 hCaloEz = _dbe->book1D("METTask_CaloEz", "METTask_CaloEz" ,500,-500,500);
00270 hCaloEz->setAxisTitle("Ez [GeV]",1);
00271 hCaloMETSig = _dbe->book1D("METTask_CaloMETSig","METTask_CaloMETSig",51,0,51);
00272 hCaloMETSig->setAxisTitle("METSig",1);
00273 hCaloMET = _dbe->book1D("METTask_CaloMET", "METTask_CaloMET" ,500,0,1000);
00274 hCaloMET->setAxisTitle("MET [GeV]",1);
00275
00276
00277 hCaloMETPhi = _dbe->book1D("METTask_CaloMETPhi","METTask_CaloMETPhi",80,-TMath::Pi(),TMath::Pi());
00278 hCaloMETPhi->setAxisTitle("METPhi [rad]",1);
00279 hCaloSumET = _dbe->book1D("METTask_CaloSumET", "METTask_CaloSumET" ,500,0,2000);
00280 hCaloSumET->setAxisTitle("SumET [GeV]",1);
00281
00282 hCaloMET_logx = _dbe->book1D("METTask_CaloMET_logx", "METTask_CaloMET_logx" ,40,-1.,7.);
00283 hCaloMET_logx->setAxisTitle("log(MET) [GeV]",1);
00284 hCaloSumET_logx = _dbe->book1D("METTask_CaloSumET_logx", "METTask_CaloSumET_logx" ,40,-1.,7.);
00285 hCaloSumET_logx->setAxisTitle("log(SumET) [GeV]",1);
00286
00287 hCaloMETIonFeedbck = _dbe->book1D("METTask_CaloMETIonFeedbck", "METTask_CaloMETIonFeedbck" ,500,0,1000);
00288 hCaloMETIonFeedbck->setAxisTitle("MET [GeV]",1);
00289 hCaloMETHPDNoise = _dbe->book1D("METTask_CaloMETHPDNoise", "METTask_CaloMETHPDNoise" ,500,0,1000);
00290 hCaloMETHPDNoise->setAxisTitle("MET [GeV]",1);
00291 hCaloMETRBXNoise = _dbe->book1D("METTask_CaloMETRBXNoise", "METTask_CaloMETRBXNoise" ,500,0,1000);
00292 hCaloMETRBXNoise->setAxisTitle("MET [GeV]",1);
00293
00294 hCaloMETPhi002 = _dbe->book1D("METTask_CaloMETPhi002","METTask_CaloMETPhi002",72,-TMath::Pi(),TMath::Pi());
00295 hCaloMETPhi002->setAxisTitle("METPhi [rad] (MET>2 GeV)",1);
00296 hCaloMETPhi010 = _dbe->book1D("METTask_CaloMETPhi010","METTask_CaloMETPhi010",72,-TMath::Pi(),TMath::Pi());
00297 hCaloMETPhi010->setAxisTitle("METPhi [rad] (MET>10 GeV)",1);
00298 hCaloMETPhi020 = _dbe->book1D("METTask_CaloMETPhi020","METTask_CaloMETPhi020",72,-TMath::Pi(),TMath::Pi());
00299 hCaloMETPhi020->setAxisTitle("METPhi [rad] (MET>20 GeV)",1);
00300
00301 if (_allhist){
00302 if (bLumiSecPlot){
00303 hCaloMExLS = _dbe->book2D("METTask_CaloMEx_LS","METTask_CaloMEx_LS",200,-200,200,50,0.,500.);
00304 hCaloMExLS->setAxisTitle("MEx [GeV]",1);
00305 hCaloMExLS->setAxisTitle("Lumi Section",2);
00306 hCaloMEyLS = _dbe->book2D("METTask_CaloMEy_LS","METTask_CaloMEy_LS",200,-200,200,50,0.,500.);
00307 hCaloMEyLS->setAxisTitle("MEy [GeV]",1);
00308 hCaloMEyLS->setAxisTitle("Lumi Section",2);
00309 }
00310
00311 hCaloMaxEtInEmTowers = _dbe->book1D("METTask_CaloMaxEtInEmTowers", "METTask_CaloMaxEtInEmTowers" ,100,0,2000);
00312 hCaloMaxEtInEmTowers->setAxisTitle("Et(Max) in EM Tower [GeV]",1);
00313 hCaloMaxEtInHadTowers = _dbe->book1D("METTask_CaloMaxEtInHadTowers", "METTask_CaloMaxEtInHadTowers" ,100,0,2000);
00314 hCaloMaxEtInHadTowers->setAxisTitle("Et(Max) in Had Tower [GeV]",1);
00315 hCaloEtFractionHadronic = _dbe->book1D("METTask_CaloEtFractionHadronic","METTask_CaloEtFractionHadronic",100,0,1);
00316 hCaloEtFractionHadronic->setAxisTitle("Hadronic Et Fraction",1);
00317 hCaloEmEtFraction = _dbe->book1D("METTask_CaloEmEtFraction", "METTask_CaloEmEtFraction" ,100,0,1);
00318 hCaloEmEtFraction->setAxisTitle("EM Et Fraction",1);
00319
00320 hCaloEmEtFraction002 = _dbe->book1D("METTask_CaloEmEtFraction002", "METTask_CaloEmEtFraction002" ,100,0,1);
00321 hCaloEmEtFraction002->setAxisTitle("EM Et Fraction (MET>2 GeV)",1);
00322 hCaloEmEtFraction010 = _dbe->book1D("METTask_CaloEmEtFraction010", "METTask_CaloEmEtFraction010" ,100,0,1);
00323 hCaloEmEtFraction010->setAxisTitle("EM Et Fraction (MET>10 GeV)",1);
00324 hCaloEmEtFraction020 = _dbe->book1D("METTask_CaloEmEtFraction020", "METTask_CaloEmEtFraction020" ,100,0,1);
00325 hCaloEmEtFraction020->setAxisTitle("EM Et Fraction (MET>20 GeV)",1);
00326
00327 hCaloHadEtInHB = _dbe->book1D("METTask_CaloHadEtInHB","METTask_CaloHadEtInHB",100,0,2000);
00328 hCaloHadEtInHB->setAxisTitle("Had Et [GeV]",1);
00329 hCaloHadEtInHO = _dbe->book1D("METTask_CaloHadEtInHO","METTask_CaloHadEtInHO",100,0,2000);
00330 hCaloHadEtInHO->setAxisTitle("Had Et [GeV]",1);
00331 hCaloHadEtInHE = _dbe->book1D("METTask_CaloHadEtInHE","METTask_CaloHadEtInHE",100,0,2000);
00332 hCaloHadEtInHE->setAxisTitle("Had Et [GeV]",1);
00333 hCaloHadEtInHF = _dbe->book1D("METTask_CaloHadEtInHF","METTask_CaloHadEtInHF",100,0,2000);
00334 hCaloHadEtInHF->setAxisTitle("Had Et [GeV]",1);
00335 hCaloEmEtInHF = _dbe->book1D("METTask_CaloEmEtInHF" ,"METTask_CaloEmEtInHF" ,100,0,2000);
00336 hCaloEmEtInHF->setAxisTitle("EM Et [GeV]",1);
00337 hCaloEmEtInEE = _dbe->book1D("METTask_CaloEmEtInEE" ,"METTask_CaloEmEtInEE" ,100,0,2000);
00338 hCaloEmEtInEE->setAxisTitle("EM Et [GeV]",1);
00339 hCaloEmEtInEB = _dbe->book1D("METTask_CaloEmEtInEB" ,"METTask_CaloEmEtInEB" ,100,0,2000);
00340 hCaloEmEtInEB->setAxisTitle("EM Et [GeV]",1);
00341
00342 hCaloEmMEx= _dbe->book1D("METTask_CaloEmMEx","METTask_CaloEmMEx",500,-500,500);
00343 hCaloEmMEx->setAxisTitle("EM MEx [GeV]",1);
00344 hCaloEmMEy= _dbe->book1D("METTask_CaloEmMEy","METTask_CaloEmMEy",500,-500,500);
00345 hCaloEmMEy->setAxisTitle("EM MEy [GeV]",1);
00346 hCaloEmEz= _dbe->book1D("METTask_CaloEmEz","METTask_CaloEmEz",500,-500,500);
00347 hCaloEmEz->setAxisTitle("EM Ez [GeV]",1);
00348 hCaloEmMET= _dbe->book1D("METTask_CaloEmMET","METTask_CaloEmMET",500,0,1000);
00349 hCaloEmMET->setAxisTitle("EM MET [GeV]",1);
00350 hCaloEmMETPhi= _dbe->book1D("METTask_CaloEmMETPhi","METTask_CaloEmMETPhi",80,-TMath::Pi(),TMath::Pi());
00351 hCaloEmMETPhi->setAxisTitle("EM METPhi [rad]",1);
00352 hCaloEmSumET= _dbe->book1D("METTask_CaloEmSumET","METTask_CaloEmSumET",500,0,2000);
00353 hCaloEmSumET->setAxisTitle("EM SumET [GeV]",1);
00354
00355 hCaloHaMEx= _dbe->book1D("METTask_CaloHaMEx","METTask_CaloHaMEx",500,-500,500);
00356 hCaloHaMEx->setAxisTitle("HA MEx [GeV]",1);
00357 hCaloHaMEy= _dbe->book1D("METTask_CaloHaMEy","METTask_CaloHaMEy",500,-500,500);
00358 hCaloHaMEy->setAxisTitle("HA MEy [GeV]",1);
00359 hCaloHaEz= _dbe->book1D("METTask_CaloHaEz","METTask_CaloHaEz",500,-500,500);
00360 hCaloHaEz->setAxisTitle("HA Ez [GeV]",1);
00361 hCaloHaMET= _dbe->book1D("METTask_CaloHaMET","METTask_CaloHaMET",500,0,1000);
00362 hCaloHaMET->setAxisTitle("HA MET [GeV]",1);
00363 hCaloHaMETPhi= _dbe->book1D("METTask_CaloHaMETPhi","METTask_CaloHaMETPhi",80,-TMath::Pi(),TMath::Pi());
00364 hCaloHaMETPhi->setAxisTitle("HA METPhi [rad]",1);
00365 hCaloHaSumET= _dbe->book1D("METTask_CaloHaSumET","METTask_CaloHaSumET",500,0,2000);
00366 hCaloHaSumET->setAxisTitle("HA SumET [GeV]",1);
00367 }
00368
00369 if (theCaloMETCollectionLabel.label() == "corMetGlobalMuons" ) {
00370 hCalomuPt = _dbe->book1D("METTask_CalomuonPt", "METTask_CalomuonPt", 50, 0, 500);
00371 hCalomuEta = _dbe->book1D("METTask_CalomuonEta", "METTask_CalomuonEta", 60, -3.0, 3.0);
00372 hCalomuNhits = _dbe->book1D("METTask_CalomuonNhits", "METTask_CalomuonNhits", 50, 0, 50);
00373 hCalomuChi2 = _dbe->book1D("METTask_CalomuonNormalizedChi2", "METTask_CalomuonNormalizedChi2", 20, 0, 20);
00374 hCalomuD0 = _dbe->book1D("METTask_CalomuonD0", "METTask_CalomuonD0", 50, -1, 1);
00375 hCaloMExCorrection = _dbe->book1D("METTask_CaloMExCorrection", "METTask_CaloMExCorrection", 100, -500.0,500.0);
00376 hCaloMEyCorrection = _dbe->book1D("METTask_CaloMEyCorrection", "METTask_CaloMEyCorrection", 100, -500.0,500.0);
00377 hCaloMuonCorrectionFlag = _dbe->book1D("METTask_CaloCorrectionFlag","METTask_CaloCorrectionFlag", 5, -0.5, 4.5);
00378 }
00379
00380 }
00381
00382
00383 void CaloMETAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00384 {
00385 if ( _HighPtJetEventFlag->on() ) _HighPtJetEventFlag->initRun( iRun, iSetup );
00386 if ( _LowPtJetEventFlag ->on() ) _LowPtJetEventFlag ->initRun( iRun, iSetup );
00387 if ( _MinBiasEventFlag ->on() ) _MinBiasEventFlag ->initRun( iRun, iSetup );
00388 if ( _HighMETEventFlag ->on() ) _HighMETEventFlag ->initRun( iRun, iSetup );
00389 if ( _LowMETEventFlag ->on() ) _LowMETEventFlag ->initRun( iRun, iSetup );
00390 if ( _EleEventFlag ->on() ) _EleEventFlag ->initRun( iRun, iSetup );
00391 if ( _MuonEventFlag ->on() ) _MuonEventFlag ->initRun( iRun, iSetup );
00392
00393 }
00394
00395
00396 void CaloMETAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup, DQMStore * dbe)
00397 {
00398
00399
00400
00401
00402 std::string dirName = "JetMET/MET/"+_source+"/";
00403 _dbe->setCurrentFolder(dirName);
00404
00405 TH1F* tlumisec;
00406
00407 MonitorElement *meLumiSec = _dbe->get("aaa");
00408 meLumiSec = _dbe->get("JetMET/lumisec");
00409
00410 int totlsec=0;
00411 double totltime=0.;
00412 if ( meLumiSec->getRootObject() ) {
00413 tlumisec = meLumiSec->getTH1F();
00414 for (int i=0; i<500; i++){
00415 if (tlumisec->GetBinContent(i+1)) totlsec++;
00416 }
00417 totltime = double(totlsec*90);
00418 }
00419
00420 if (totltime==0.) totltime=1.;
00421
00422
00423
00424
00425 for (std::vector<std::string>::const_iterator ic = _FolderNames.begin(); ic != _FolderNames.end(); ic++)
00426 {
00427
00428 std::string DirName;
00429 DirName = dirName+*ic;
00430
00431 makeRatePlot(DirName,totltime);
00432 if ( _HighPtJetEventFlag->on() )
00433 makeRatePlot(DirName+"/"+"triggerName_HighJetPt",totltime);
00434 if ( _LowPtJetEventFlag->on() )
00435 makeRatePlot(DirName+"/"+"triggerName_LowJetPt",totltime);
00436 if ( _MinBiasEventFlag->on() )
00437 makeRatePlot(DirName+"/"+"triggerName_MinBias",totltime);
00438 if ( _HighMETEventFlag->on() )
00439 makeRatePlot(DirName+"/"+"triggerName_HighMET",totltime);
00440 if ( _LowMETEventFlag->on() )
00441 makeRatePlot(DirName+"/"+"triggerName_LowMET",totltime);
00442 if ( _EleEventFlag->on() )
00443 makeRatePlot(DirName+"/"+"triggerName_Ele",totltime);
00444 if ( _MuonEventFlag->on() )
00445 makeRatePlot(DirName+"/"+"triggerName_Muon",totltime);
00446 }
00447
00448 }
00449
00450
00451 void CaloMETAnalyzer::makeRatePlot(std::string DirName, double totltime)
00452 {
00453
00454 _dbe->setCurrentFolder(DirName);
00455 MonitorElement *meCaloMET = _dbe->get(DirName+"/"+"METTask_CaloMET");
00456
00457 TH1F* tCaloMET;
00458 TH1F* tCaloMETRate;
00459
00460 if ( meCaloMET )
00461 if ( meCaloMET->getRootObject() ) {
00462 tCaloMET = meCaloMET->getTH1F();
00463
00464
00465 tCaloMETRate = (TH1F*) tCaloMET->Clone("METTask_CaloMETRate");
00466 for (int i = tCaloMETRate->GetNbinsX()-1; i>=0; i--){
00467 tCaloMETRate->SetBinContent(i+1,tCaloMETRate->GetBinContent(i+2)+tCaloMET->GetBinContent(i+1));
00468 }
00469 for (int i = 0; i<tCaloMETRate->GetNbinsX(); i++){
00470 tCaloMETRate->SetBinContent(i+1,tCaloMETRate->GetBinContent(i+1)/double(totltime));
00471 }
00472
00473 tCaloMETRate->SetName("METTask_CaloMETRate");
00474 tCaloMETRate->SetTitle("METTask_CaloMETRate");
00475 hCaloMETRate = _dbe->book1D("METTask_CaloMETRate",tCaloMETRate);
00476 hCaloMETRate->setAxisTitle("MET Threshold [GeV]",1);
00477 }
00478 }
00479
00480
00481
00482 void CaloMETAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup,
00483 const edm::TriggerResults& triggerResults) {
00484
00485 if (_verbose) std::cout << "CaloMETAnalyzer analyze" << std::endl;
00486
00487 std::string DirName = "JetMET/MET/"+_source;
00488
00489 if (_print){
00490 std::cout << " " << std::endl;
00491 std::cout << "Event = " << iEvent.id().event() << std::endl;
00492 }
00493
00494 LogTrace(metname)<<"[CaloMETAnalyzer] Analyze CaloMET";
00495
00496 hmetME->Fill(1);
00497
00498
00499
00500
00501 _trig_JetMB=0;
00502 _trig_HighPtJet=0;
00503 _trig_LowPtJet=0;
00504 _trig_MinBias=0;
00505 _trig_HighMET=0;
00506 _trig_LowMET=0;
00507 _trig_Ele=0;
00508 _trig_Muon=0;
00509 _trig_PhysDec=0;
00510 if(&triggerResults) {
00511
00513
00514
00515
00516
00517 int ntrigs = triggerResults.size();
00518 if (_verbose) std::cout << "ntrigs=" << ntrigs << std::endl;
00519
00520
00521
00522
00523 const edm::TriggerNames & triggerNames = iEvent.triggerNames(triggerResults);
00524
00525
00526
00527
00528 for (unsigned int i=0; i!=HLTPathsJetMBByName_.size(); i++) {
00529 unsigned int triggerIndex = triggerNames.triggerIndex(HLTPathsJetMBByName_[i]);
00530 if (triggerIndex<triggerResults.size()) {
00531 if (triggerResults.accept(triggerIndex)) {
00532 _trig_JetMB++;
00533 }
00534 }
00535 }
00536
00537 if (HLTPathsJetMBByName_.size()==0) _trig_JetMB=triggerResults.size()-1;
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550 if ( _HighPtJetEventFlag->on() && _HighPtJetEventFlag->accept( iEvent, iSetup ) )
00551 _trig_HighPtJet=1;
00552
00553 if ( _LowPtJetEventFlag->on() && _LowPtJetEventFlag->accept( iEvent, iSetup ) )
00554 _trig_LowPtJet=1;
00555
00556 if ( _MinBiasEventFlag->on() && _MinBiasEventFlag->accept( iEvent, iSetup ) )
00557 _trig_MinBias=1;
00558
00559 if ( _HighMETEventFlag->on() && _HighMETEventFlag->accept( iEvent, iSetup ) )
00560 _trig_HighMET=1;
00561
00562 if ( _LowMETEventFlag->on() && _LowMETEventFlag->accept( iEvent, iSetup ) )
00563 _trig_LowMET=1;
00564
00565 if ( _EleEventFlag->on() && _EleEventFlag->accept( iEvent, iSetup ) )
00566 _trig_Ele=1;
00567
00568 if ( _MuonEventFlag->on() && _MuonEventFlag->accept( iEvent, iSetup ) )
00569 _trig_Muon=1;
00570
00571 if (triggerNames.triggerIndex(_hlt_PhysDec) != triggerNames.size() &&
00572 triggerResults.accept(triggerNames.triggerIndex(_hlt_PhysDec))) _trig_PhysDec=1;
00573
00574 } else {
00575
00576 edm::LogInfo("CaloMetAnalyzer") << "TriggerResults::HLT not found, "
00577 "automatically select events";
00578
00579
00580 _trig_JetMB=1;
00581
00582 }
00583
00584
00585
00586
00587
00588 edm::Handle<reco::CaloMETCollection> calometcoll;
00589 iEvent.getByLabel(theCaloMETCollectionLabel, calometcoll);
00590
00591 if(!calometcoll.isValid()) {
00592 std::cout<<"Unable to find MET results for CaloMET collection "<<theCaloMETCollectionLabel<<std::endl;
00593 return;
00594 }
00595
00596 const reco::CaloMETCollection *calometcol = calometcoll.product();
00597 const reco::CaloMET *calomet;
00598 calomet = &(calometcol->front());
00599
00600 LogTrace(metname)<<"[CaloMETAnalyzer] Call to the CaloMET analyzer";
00601
00602
00603 if (theCaloMETCollectionLabel.label() == "corMetGlobalMuons" ) {
00604
00605 iEvent.getByLabel("muonMETValueMapProducer" , "muCorrData", corMetGlobalMuons_ValueMap_Handle);
00606 iEvent.getByLabel("muons", muon_h);
00607 iEvent.getByLabel(inputBeamSpotLabel, beamSpot_h);
00608
00609 if(!beamSpot_h.isValid()) edm::LogInfo("OutputInfo") << "falied to retrieve beam spot data require by MET Task";
00610
00611 bspot = ( beamSpot_h.isValid() ) ? beamSpot_h->position() : math::XYZPoint(0, 0, 0);
00612
00613 }
00614
00615
00616
00617
00618 edm::Handle<reco::HcalNoiseRBXCollection> HRBXCollection;
00619 iEvent.getByLabel(HcalNoiseRBXCollectionTag,HRBXCollection);
00620 if (!HRBXCollection.isValid()) {
00621 LogDebug("") << "CaloMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00622 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find HcalNoiseRBX Collection" << std::endl;
00623 }
00624
00625 edm::Handle<HcalNoiseSummary> HNoiseSummary;
00626 iEvent.getByLabel(HcalNoiseSummaryTag,HNoiseSummary);
00627 if (!HNoiseSummary.isValid()) {
00628 LogDebug("") << "CaloMETAnalyzer: Could not find Hcal NoiseSummary product" << std::endl;
00629 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find Hcal NoiseSummary product" << std::endl;
00630 }
00631
00632 edm::Handle<reco::CaloJetCollection> caloJets;
00633 iEvent.getByLabel(theJetCollectionLabel, caloJets);
00634 if (!caloJets.isValid()) {
00635 LogDebug("") << "CaloMETAnalyzer: Could not find jet product" << std::endl;
00636 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find jet product" << std::endl;
00637 }
00638
00639 edm::Handle<edm::View<reco::Candidate> > towers;
00640 iEvent.getByLabel(theCaloTowersLabel, towers);
00641 if (!towers.isValid()) {
00642 LogDebug("") << "CaloMETAnalyzer: Could not find caltower product" << std::endl;
00643 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find caltower product" << std::endl;
00644 }
00645
00646
00647
00648
00649 if (_source=="CaloMET") validateMET(*calomet,towers);
00650
00651
00652
00653 if (_allhist) computeEmHaMET(towers);
00654
00655
00656
00657
00658 if (_verbose) std::cout << "JetID starts" << std::endl;
00659
00660
00661
00662
00663 bool bJetIDMinimal=true;
00664 int nj=0;
00665 for (reco::CaloJetCollection::const_iterator cal = caloJets->begin();
00666 cal!=caloJets->end(); ++cal){
00667 jetID->calculate(iEvent, *cal);
00668 if (_print && nj<=1) std::cout << "Jet pT = " << cal->pt() << " (GeV) "
00669 << " eta = " << cal->eta() << " "
00670 << " phi = " << cal->phi() << " "
00671 << " emf = " << cal->emEnergyFraction() << std::endl;
00672 nj++;
00673 if (cal->pt()>10.){
00674 if (fabs(cal->eta())<=2.6 &&
00675 cal->emEnergyFraction()<=0.01) bJetIDMinimal=false;
00676 }
00677 }
00678
00679
00680
00681
00682 bool bJetIDLoose=true;
00683 for (reco::CaloJetCollection::const_iterator cal = caloJets->begin();
00684 cal!=caloJets->end(); ++cal){
00685 jetID->calculate(iEvent, *cal);
00686 if (_verbose) std::cout << jetID->n90Hits() << " "
00687 << jetID->restrictedEMF() << " "
00688 << cal->pt() << std::endl;
00689 if (cal->pt()>10.){
00690
00691
00692 if (jetID->n90Hits()<2) bJetIDLoose=false;
00693 if (jetID->fHPD()>=0.98) bJetIDLoose=false;
00694
00695
00696 if (fabs(cal->eta())<2.55){
00697 if (cal->emEnergyFraction()<=0.01) bJetIDLoose=false;
00698 }
00699
00700 else {
00701 if (cal->emEnergyFraction()<=-0.9) bJetIDLoose=false;
00702 if (cal->pt()>80.){
00703 if (cal->emEnergyFraction()>= 1.0) bJetIDLoose=false;
00704 }
00705 }
00706 }
00707 }
00708
00709
00710
00711
00712 bool bJetIDTight=true;
00713 bJetIDTight=bJetIDLoose;
00714 for (reco::CaloJetCollection::const_iterator cal = caloJets->begin();
00715 cal!=caloJets->end(); ++cal){
00716 jetID->calculate(iEvent, *cal);
00717 if (cal->pt()>25.){
00718
00719
00720 if (jetID->fHPD()>=0.95) bJetIDTight=false;
00721
00722
00723 if (fabs(cal->eta())>=1.00 && fabs(cal->eta())<1.75){
00724 if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false;
00725 }
00726
00727
00728 else if (fabs(cal->eta())>=1.75 && fabs(cal->eta())<2.55){
00729 if (cal->pt()>80. && cal->emEnergyFraction()>=1.) bJetIDTight=false;
00730 }
00731
00732
00733 else if (fabs(cal->eta())>=2.55 && fabs(cal->eta())<3.25){
00734 if (cal->pt()< 50. && cal->emEnergyFraction()<=-0.3) bJetIDTight=false;
00735 if (cal->pt()>=50. && cal->pt()< 80. && cal->emEnergyFraction()<=-0.2) bJetIDTight=false;
00736 if (cal->pt()>=80. && cal->pt()<340. && cal->emEnergyFraction()<=-0.1) bJetIDTight=false;
00737 if (cal->pt()>=340. && cal->emEnergyFraction()<=-0.1
00738 && cal->emEnergyFraction()>=0.95) bJetIDTight=false;
00739 }
00740
00741
00742 else if (fabs(cal->eta())>=3.25){
00743 if (cal->pt()< 50. && cal->emEnergyFraction()<=-0.3
00744 && cal->emEnergyFraction()>=0.90) bJetIDTight=false;
00745 if (cal->pt()>=50. && cal->pt()<130. && cal->emEnergyFraction()<=-0.2
00746 && cal->emEnergyFraction()>=0.80) bJetIDTight=false;
00747 if (cal->pt()>=130. && cal->emEnergyFraction()<=-0.1
00748 && cal->emEnergyFraction()>=0.70) bJetIDTight=false;
00749 }
00750 }
00751 }
00752
00753 if (_verbose) std::cout << "JetID ends" << std::endl;
00754
00755
00756
00757
00758 bool bHcalNoiseFilter = HNoiseSummary->passLooseNoiseFilter();
00759 bool bHcalNoiseFilterTight = HNoiseSummary->passTightNoiseFilter();
00760
00761 if (_verbose) std::cout << "HcalNoiseFilter Summary ends" << std::endl;
00762
00763
00764
00765 edm::Handle<reco::BeamHaloSummary> TheBeamHaloSummary ;
00766 iEvent.getByLabel(BeamHaloSummaryTag, TheBeamHaloSummary) ;
00767
00768 bool bBeamHaloIDTightPass = true;
00769 bool bBeamHaloIDLoosePass = true;
00770
00771 if(TheBeamHaloSummary.isValid()) {
00772
00773 const reco::BeamHaloSummary TheSummary = (*TheBeamHaloSummary.product() );
00774
00775
00776
00777
00778
00779
00780 if( TheSummary.EcalLooseHaloId() || TheSummary.HcalLooseHaloId() ||
00781 TheSummary.CSCLooseHaloId() || TheSummary.GlobalLooseHaloId() )
00782 bBeamHaloIDLoosePass = false;
00783
00784 if( TheSummary.EcalTightHaloId() || TheSummary.HcalTightHaloId() ||
00785 TheSummary.CSCTightHaloId() || TheSummary.GlobalTightHaloId() )
00786 bBeamHaloIDTightPass = false;
00787
00788 }
00789
00790 if (_verbose) std::cout << "BeamHaloSummary ends" << std::endl;
00791
00792
00793
00794
00795 bool bPrimaryVertex = true;
00796 if(_doPVCheck){
00797 bPrimaryVertex = false;
00798 Handle<reco::VertexCollection> vertexHandle;
00799
00800 iEvent.getByLabel(vertexTag, vertexHandle);
00801
00802 if (!vertexHandle.isValid()) {
00803 LogDebug("") << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00804 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find vertex collection" << std::endl;
00805 }
00806
00807 if ( vertexHandle.isValid() ){
00808 reco::VertexCollection vertexCollection = *(vertexHandle.product());
00809 int vertex_number = vertexCollection.size();
00810 reco::VertexCollection::const_iterator v = vertexCollection.begin();
00811 double vertex_chi2 = v->normalizedChi2();
00812 double vertex_ndof = v->ndof();
00813 bool fakeVtx = v->isFake();
00814 double vertex_Z = v->z();
00815
00816 if ( !fakeVtx
00817 && vertex_number>=_nvtx_min
00818 && vertex_ndof >_vtxndof_min
00819 && vertex_chi2 <_vtxchi2_max
00820 && fabs(vertex_Z)<_vtxz_max )
00821 bPrimaryVertex = true;
00822 }
00823 }
00824
00825
00826 edm::Handle< L1GlobalTriggerReadoutRecord > gtReadoutRecord;
00827 iEvent.getByLabel( gtTag, gtReadoutRecord);
00828
00829 if (!gtReadoutRecord.isValid()) {
00830 LogDebug("") << "CaloMETAnalyzer: Could not find GT readout record" << std::endl;
00831 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find GT readout record product" << std::endl;
00832 }
00833
00834 bool bTechTriggers = true;
00835 bool bTechTriggersAND = true;
00836 bool bTechTriggersOR = false;
00837 bool bTechTriggersNOT = false;
00838
00839 if (gtReadoutRecord.isValid()) {
00840 const TechnicalTriggerWord& technicalTriggerWordBeforeMask = gtReadoutRecord->technicalTriggerWord();
00841
00842 if (_techTrigsAND.size() == 0)
00843 bTechTriggersAND = true;
00844 else
00845 for (unsigned ttr = 0; ttr != _techTrigsAND.size(); ttr++) {
00846 bTechTriggersAND = bTechTriggersAND && technicalTriggerWordBeforeMask.at(_techTrigsAND.at(ttr));
00847 }
00848
00849 if (_techTrigsAND.size() == 0)
00850 bTechTriggersOR = true;
00851 else
00852 for (unsigned ttr = 0; ttr != _techTrigsOR.size(); ttr++) {
00853 bTechTriggersOR = bTechTriggersOR || technicalTriggerWordBeforeMask.at(_techTrigsOR.at(ttr));
00854 }
00855 if (_techTrigsNOT.size() == 0)
00856 bTechTriggersNOT = false;
00857 else
00858 for (unsigned ttr = 0; ttr != _techTrigsNOT.size(); ttr++) {
00859 bTechTriggersNOT = bTechTriggersNOT || technicalTriggerWordBeforeMask.at(_techTrigsNOT.at(ttr));
00860 }
00861 }
00862 else
00863 {
00864 bTechTriggersAND = true;
00865 bTechTriggersOR = true;
00866 bTechTriggersNOT = false;
00867 }
00868
00869 if (_techTrigsAND.size()==0)
00870 bTechTriggersAND = true;
00871 if (_techTrigsOR.size()==0)
00872 bTechTriggersOR = true;
00873 if (_techTrigsNOT.size()==0)
00874 bTechTriggersNOT = false;
00875
00876 bTechTriggers = bTechTriggersAND && bTechTriggersOR && !bTechTriggersNOT;
00877
00878
00879
00880
00881 bool bHcalNoise = bHcalNoiseFilter;
00882 bool bBeamHaloID = bBeamHaloIDLoosePass;
00883 bool bJetID = bJetIDMinimal;
00884
00885 bool bPhysicsDeclared = true;
00886 if(_doHLTPhysicsOn) bPhysicsDeclared =_trig_PhysDec;
00887
00888 if (_tightHcalFiltering) bHcalNoise = bHcalNoiseFilterTight;
00889 if (_tightBHFiltering) bBeamHaloID = bBeamHaloIDTightPass;
00890
00891 if (_tightJetIDFiltering==1) bJetID = bJetIDMinimal;
00892 else if (_tightJetIDFiltering==2) bJetID = bJetIDLoose;
00893 else if (_tightJetIDFiltering==3) bJetID = bJetIDTight;
00894 else if (_tightJetIDFiltering==-1) bJetID = true;
00895
00896 bool bBasicCleanup = bTechTriggers && bPrimaryVertex && bPhysicsDeclared;
00897 bool bExtraCleanup = bBasicCleanup && bHcalNoise && bJetID && bBeamHaloID;
00898
00899
00900
00901 for (std::vector<std::string>::const_iterator ic = _FolderNames.begin();
00902 ic != _FolderNames.end(); ic++){
00903 if (*ic=="All") fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00904 if (DCSFilter->filter(iEvent, iSetup)) {
00905 if (_cleanupSelection){
00906 if (*ic=="BasicCleanup" && bBasicCleanup) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00907 if (*ic=="ExtraCleanup" && bExtraCleanup) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00908 }
00909 if (_allSelection) {
00910 if (*ic=="HcalNoiseFilter" && bHcalNoiseFilter ) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00911 if (*ic=="HcalNoiseFilterTight" && bHcalNoiseFilterTight ) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00912 if (*ic=="JetIDMinimal" && bJetIDMinimal) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00913 if (*ic=="JetIDLoose" && bJetIDLoose) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00914 if (*ic=="JetIDTight" && bJetIDTight) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00915 if (*ic=="BeamHaloIDTightPass" && bBeamHaloIDTightPass) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00916 if (*ic=="BeamHaloIDLoosePass" && bBeamHaloIDLoosePass) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00917 if (*ic=="Triggers" && bTechTriggers) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00918 if (*ic=="PV" && bPrimaryVertex) fillMESet(iEvent, DirName+"/"+*ic, *calomet);
00919 }
00920 }
00921 }
00922 }
00923
00924
00925 void CaloMETAnalyzer::computeEmHaMET(edm::Handle<edm::View<reco::Candidate> > towers)
00926 {
00927
00928 edm::View<reco::Candidate>::const_iterator towerCand = towers->begin();
00929
00930 double sum_em_et = 0.0;
00931 double sum_em_ex = 0.0;
00932 double sum_em_ey = 0.0;
00933 double sum_em_ez = 0.0;
00934
00935 double sum_ha_et = 0.0;
00936 double sum_ha_ex = 0.0;
00937 double sum_ha_ey = 0.0;
00938 double sum_ha_ez = 0.0;
00939
00940 for ( ; towerCand != towers->end(); towerCand++)
00941 {
00942 const reco::Candidate* candidate = &(*towerCand);
00943 if (candidate)
00944 {
00945 const CaloTower* calotower = dynamic_cast<const CaloTower*> (candidate);
00946 if (calotower){
00947 double Tower_ET = calotower->et();
00948 if (Tower_ET>0.3) {
00949
00950 double phi = candidate->phi();
00951 double theta = candidate->theta();
00952
00953 double e_em = calotower->emEnergy();
00954 double e_ha = calotower->hadEnergy();
00955 double et_em = e_em*sin(theta);
00956 double et_ha = e_ha*sin(theta);
00957
00958 sum_em_ez += e_em*cos(theta);
00959 sum_em_et += et_em;
00960 sum_em_ex += et_em*cos(phi);
00961 sum_em_ey += et_em*sin(phi);
00962
00963 sum_ha_ez += e_ha*cos(theta);
00964 sum_ha_et += et_ha;
00965 sum_ha_ex += et_ha*cos(phi);
00966 sum_ha_ey += et_ha*sin(phi);
00967
00968 }
00969 }
00970 }
00971 }
00972
00973
00974 _EmMEx = -sum_em_ex;
00975 _EmMEy = -sum_em_ey;
00976 _EmMET = pow(_EmMEx*_EmMEx+_EmMEy*_EmMEy,0.5);
00977 _EmCaloEz = sum_em_ez;
00978 _EmSumEt = sum_em_et;
00979 _EmMetPhi = atan2( _EmMEy, _EmMEx );
00980
00981 _HaMEx = -sum_ha_ex;
00982 _HaMEy = -sum_ha_ey;
00983 _HaMET = pow(_HaMEx*_HaMEx+_HaMEy*_HaMEy,0.5);
00984 _HaCaloEz = sum_ha_ez;
00985 _HaSumEt = sum_ha_et;
00986 _HaMetPhi = atan2( _HaMEy, _HaMEx );
00987
00988 }
00989
00990 void CaloMETAnalyzer::validateMET(const reco::CaloMET& calomet,
00991 edm::Handle<edm::View<reco::Candidate> > towers)
00992 {
00993
00994 edm::View<reco::Candidate>::const_iterator towerCand = towers->begin();
00995
00996 double sum_et = 0.0;
00997 double sum_ex = 0.0;
00998 double sum_ey = 0.0;
00999 double sum_ez = 0.0;
01000
01001 for ( ; towerCand != towers->end(); towerCand++)
01002 {
01003 const reco::Candidate* candidate = &(*towerCand);
01004 if (candidate)
01005 {
01006 const CaloTower* calotower = dynamic_cast<const CaloTower*> (candidate);
01007 if (calotower){
01008 double Tower_ET = calotower->et();
01009 if (Tower_ET>0.3) {
01010
01011 double phi = candidate->phi();
01012 double theta = candidate->theta();
01013 double e = candidate->energy();
01014 double et = e*sin(theta);
01015 sum_ez += e*cos(theta);
01016 sum_et += et;
01017 sum_ex += et*cos(phi);
01018 sum_ey += et*sin(phi);
01019
01020 }
01021 }
01022 }
01023 }
01024
01025 double Mex = -sum_ex;
01026 double Mey = -sum_ey;
01027
01028 double Met = sqrt( sum_ex*sum_ex + sum_ey*sum_ey );
01029 double Sumet = sum_et;
01030
01031
01032 if (_verbose){
01033 if (Sumet!=calomet.sumEt() || Mex!=calomet.px() || Mey!=calomet.py() || Met!=calomet.pt() ){
01034 std::cout << _source << std::endl;
01035 std::cout << "SUMET" << Sumet << " METBlock" << calomet.sumEt() << std::endl;
01036 std::cout << "MEX" << Mex << " METBlock" << calomet.px() << std::endl;
01037 std::cout << "MEY" << Mey << " METBlock" << calomet.py() << std::endl;
01038 std::cout << "MET" << Met << " METBlock" << calomet.pt() << std::endl;
01039 }
01040 }
01041
01042 if (_print){
01043 std::cout << "SUMET = " << calomet.sumEt() << " (GeV) "
01044 << "MEX" << calomet.px() << " (GeV) "
01045 << "MEY" << calomet.py() << " (GeV) "
01046 << "MET" << calomet.pt() << " (GeV) " << std::endl;
01047 }
01048
01049 }
01050
01051
01052 void CaloMETAnalyzer::fillMESet(const edm::Event& iEvent, std::string DirName,
01053 const reco::CaloMET& calomet)
01054 {
01055
01056 _dbe->setCurrentFolder(DirName);
01057
01058 bool bLumiSecPlot=false;
01059 if (DirName.find("All")) bLumiSecPlot=true;
01060
01061 if (_trig_JetMB) fillMonitorElement(iEvent,DirName,"",calomet, bLumiSecPlot);
01062
01063 if (_trig_HighPtJet)
01064 fillMonitorElement(iEvent,DirName,"HighPtJet",calomet,false);
01065
01066 if (_trig_LowPtJet)
01067 fillMonitorElement(iEvent,DirName,"LowPtJet",calomet,false);
01068
01069 if (_trig_MinBias)
01070 fillMonitorElement(iEvent,DirName,"MinBias",calomet,false);
01071
01072 if (_trig_HighMET)
01073 fillMonitorElement(iEvent,DirName,"HighMET",calomet,false);
01074
01075 if (_trig_LowMET)
01076 fillMonitorElement(iEvent,DirName,"LowMET",calomet,false);
01077
01078 if (_trig_Ele)
01079 fillMonitorElement(iEvent,DirName,"Ele",calomet,false);
01080
01081 if (_trig_Muon)
01082 fillMonitorElement(iEvent,DirName,"Muon",calomet,false);
01083
01084 }
01085
01086
01087 void CaloMETAnalyzer::fillMonitorElement(const edm::Event& iEvent, std::string DirName,
01088 std::string TriggerTypeName,
01089 const reco::CaloMET& calomet, bool bLumiSecPlot)
01090 {
01091
01092 if (TriggerTypeName=="HighPtJet") {
01093 if (!selectHighPtJetEvent(iEvent)) return;
01094 }
01095 else if (TriggerTypeName=="LowPtJet") {
01096 if (!selectLowPtJetEvent(iEvent)) return;
01097 }
01098 else if (TriggerTypeName=="HighMET") {
01099 if (calomet.pt()<_highMETThreshold) return;
01100 }
01101 else if (TriggerTypeName=="LowMET") {
01102 if (calomet.pt()<_lowMETThreshold) return;
01103 }
01104 else if (TriggerTypeName=="Ele") {
01105 if (!selectWElectronEvent(iEvent)) return;
01106 }
01107 else if (TriggerTypeName=="Muon") {
01108 if (!selectWMuonEvent(iEvent)) return;
01109 }
01110
01111 double caloSumET = calomet.sumEt();
01112 double caloMETSig = calomet.mEtSig();
01113 double caloEz = calomet.e_longitudinal();
01114 double caloMET = calomet.pt();
01115 double caloMEx = calomet.px();
01116 double caloMEy = calomet.py();
01117 double caloMETPhi = calomet.phi();
01118
01119 if (_verbose) std::cout << _source << " " << caloMET << std::endl;
01120
01121 double caloEtFractionHadronic = calomet.etFractionHadronic();
01122 double caloEmEtFraction = calomet.emEtFraction();
01123
01124 double caloMaxEtInEMTowers = calomet.maxEtInEmTowers();
01125 double caloMaxEtInHadTowers = calomet.maxEtInHadTowers();
01126
01127 double caloHadEtInHB = calomet.hadEtInHB();
01128 double caloHadEtInHO = calomet.hadEtInHO();
01129 double caloHadEtInHE = calomet.hadEtInHE();
01130 double caloHadEtInHF = calomet.hadEtInHF();
01131 double caloEmEtInEB = calomet.emEtInEB();
01132 double caloEmEtInEE = calomet.emEtInEE();
01133 double caloEmEtInHF = calomet.emEtInHF();
01134
01135
01136 int myLuminosityBlock;
01137
01138 myLuminosityBlock = iEvent.luminosityBlock();
01139
01140
01141 if (TriggerTypeName!="") DirName = DirName +"/"+TriggerTypeName;
01142
01143 if (_verbose) std::cout << "_etThreshold = " << _etThreshold << std::endl;
01144 if (caloSumET>_etThreshold){
01145
01146 hCaloMEx = _dbe->get(DirName+"/"+"METTask_CaloMEx"); if (hCaloMEx && hCaloMEx->getRootObject() ) hCaloMEx->Fill(caloMEx);
01147 hCaloMEy = _dbe->get(DirName+"/"+"METTask_CaloMEy"); if (hCaloMEy && hCaloMEy->getRootObject() ) hCaloMEy->Fill(caloMEy);
01148 hCaloMET = _dbe->get(DirName+"/"+"METTask_CaloMET"); if (hCaloMET && hCaloMET->getRootObject() ) hCaloMET->Fill(caloMET);
01149 hCaloMETPhi = _dbe->get(DirName+"/"+"METTask_CaloMETPhi"); if (hCaloMETPhi && hCaloMETPhi->getRootObject() ) hCaloMETPhi->Fill(caloMETPhi);
01150 hCaloSumET = _dbe->get(DirName+"/"+"METTask_CaloSumET"); if (hCaloSumET && hCaloSumET->getRootObject() ) hCaloSumET->Fill(caloSumET);
01151 hCaloMETSig = _dbe->get(DirName+"/"+"METTask_CaloMETSig"); if (hCaloMETSig && hCaloMETSig->getRootObject() ) hCaloMETSig->Fill(caloMETSig);
01152 hCaloEz = _dbe->get(DirName+"/"+"METTask_CaloEz"); if (hCaloEz && hCaloEz->getRootObject() ) hCaloEz->Fill(caloEz);
01153
01154 hCaloMET_logx = _dbe->get(DirName+"/"+"METTask_CaloMET_logx"); if (hCaloMET_logx && hCaloMET_logx->getRootObject() ) hCaloMET_logx->Fill(log10(caloMET));
01155 hCaloSumET_logx = _dbe->get(DirName+"/"+"METTask_CaloSumET_logx"); if (hCaloSumET_logx && hCaloSumET_logx->getRootObject() ) hCaloSumET_logx->Fill(log10(caloSumET));
01156
01157 hCaloMETIonFeedbck = _dbe->get(DirName+"/"+"METTask_CaloMETIonFeedbck"); if (hCaloMETIonFeedbck && hCaloMETIonFeedbck->getRootObject() ) hCaloMETIonFeedbck->Fill(caloMET);
01158 hCaloMETHPDNoise = _dbe->get(DirName+"/"+"METTask_CaloMETHPDNoise"); if (hCaloMETHPDNoise && hCaloMETHPDNoise->getRootObject() ) hCaloMETHPDNoise->Fill(caloMET);
01159
01160 hCaloMETPhi002 = _dbe->get(DirName+"/"+"METTask_CaloMETPhi002"); if (caloMET> 2. && hCaloMETPhi002 && hCaloMETPhi002->getRootObject()) { hCaloMETPhi002->Fill(caloMETPhi);}
01161 hCaloMETPhi010 = _dbe->get(DirName+"/"+"METTask_CaloMETPhi010"); if (caloMET> 10. && hCaloMETPhi010 && hCaloMETPhi010->getRootObject()) { hCaloMETPhi010->Fill(caloMETPhi);}
01162 hCaloMETPhi020 = _dbe->get(DirName+"/"+"METTask_CaloMETPhi020"); if (caloMET> 20. && hCaloMETPhi020 && hCaloMETPhi020->getRootObject()) { hCaloMETPhi020->Fill(caloMETPhi);}
01163
01164 if (_allhist){
01165 if (bLumiSecPlot){
01166 hCaloMExLS = _dbe->get(DirName+"/"+"METTask_CaloMExLS"); if (hCaloMExLS && hCaloMExLS->getRootObject()) hCaloMExLS->Fill(caloMEx,myLuminosityBlock);
01167 hCaloMEyLS = _dbe->get(DirName+"/"+"METTask_CaloMEyLS"); if (hCaloMEyLS && hCaloMEyLS->getRootObject()) hCaloMEyLS->Fill(caloMEy,myLuminosityBlock);
01168 }
01169
01170 hCaloEtFractionHadronic = _dbe->get(DirName+"/"+"METTask_CaloEtFractionHadronic"); if (hCaloEtFractionHadronic && hCaloEtFractionHadronic->getRootObject()) hCaloEtFractionHadronic->Fill(caloEtFractionHadronic);
01171 hCaloEmEtFraction = _dbe->get(DirName+"/"+"METTask_CaloEmEtFraction"); if (hCaloEmEtFraction && hCaloEmEtFraction->getRootObject()) hCaloEmEtFraction->Fill(caloEmEtFraction);
01172
01173 hCaloEmEtFraction002 = _dbe->get(DirName+"/"+"METTask_CaloEmEtFraction002"); if (caloMET> 2. && hCaloEmEtFraction002 && hCaloEmEtFraction002->getRootObject()) hCaloEmEtFraction002->Fill(caloEmEtFraction);
01174 hCaloEmEtFraction010 = _dbe->get(DirName+"/"+"METTask_CaloEmEtFraction010"); if (caloMET> 10. && hCaloEmEtFraction010 && hCaloEmEtFraction010->getRootObject()) hCaloEmEtFraction010->Fill(caloEmEtFraction);
01175 hCaloEmEtFraction020 = _dbe->get(DirName+"/"+"METTask_CaloEmEtFraction020"); if (caloMET> 20. && hCaloEmEtFraction020 && hCaloEmEtFraction020->getRootObject()) hCaloEmEtFraction020->Fill(caloEmEtFraction);
01176
01177 hCaloMaxEtInEmTowers = _dbe->get(DirName+"/"+"METTask_CaloMaxEtInEmTowers"); if (hCaloMaxEtInEmTowers && hCaloMaxEtInEmTowers->getRootObject()) hCaloMaxEtInEmTowers->Fill(caloMaxEtInEMTowers);
01178 hCaloMaxEtInHadTowers = _dbe->get(DirName+"/"+"METTask_CaloMaxEtInHadTowers"); if (hCaloMaxEtInHadTowers && hCaloMaxEtInHadTowers->getRootObject()) hCaloMaxEtInHadTowers->Fill(caloMaxEtInHadTowers);
01179
01180 hCaloHadEtInHB = _dbe->get(DirName+"/"+"METTask_CaloHadEtInHB"); if (hCaloHadEtInHB && hCaloHadEtInHB->getRootObject()) hCaloHadEtInHB->Fill(caloHadEtInHB);
01181 hCaloHadEtInHO = _dbe->get(DirName+"/"+"METTask_CaloHadEtInHO"); if (hCaloHadEtInHO && hCaloHadEtInHO->getRootObject()) hCaloHadEtInHO->Fill(caloHadEtInHO);
01182 hCaloHadEtInHE = _dbe->get(DirName+"/"+"METTask_CaloHadEtInHE"); if (hCaloHadEtInHE && hCaloHadEtInHE->getRootObject()) hCaloHadEtInHE->Fill(caloHadEtInHE);
01183 hCaloHadEtInHF = _dbe->get(DirName+"/"+"METTask_CaloHadEtInHF"); if (hCaloHadEtInHF && hCaloHadEtInHF->getRootObject()) hCaloHadEtInHF->Fill(caloHadEtInHF);
01184 hCaloEmEtInEB = _dbe->get(DirName+"/"+"METTask_CaloEmEtInEB"); if (hCaloEmEtInEB && hCaloEmEtInEB->getRootObject()) hCaloEmEtInEB->Fill(caloEmEtInEB);
01185 hCaloEmEtInEE = _dbe->get(DirName+"/"+"METTask_CaloEmEtInEE"); if (hCaloEmEtInEE && hCaloEmEtInEE->getRootObject()) hCaloEmEtInEE->Fill(caloEmEtInEE);
01186 hCaloEmEtInHF = _dbe->get(DirName+"/"+"METTask_CaloEmEtInHF"); if (hCaloEmEtInHF && hCaloEmEtInHF->getRootObject()) hCaloEmEtInHF->Fill(caloEmEtInHF);
01187
01188 hCaloEmMEx = _dbe->get(DirName+"/"+"METTask_CaloEmMEx"); if (hCaloEmMEx && hCaloEmMEx->getRootObject()) hCaloEmMEx->Fill(_EmMEx);
01189 hCaloEmMEy = _dbe->get(DirName+"/"+"METTask_CaloEmMEy"); if (hCaloEmMEy && hCaloEmMEy->getRootObject()) hCaloEmMEy->Fill(_EmMEy);
01190 hCaloEmEz = _dbe->get(DirName+"/"+"METTask_CaloEmEz"); if (hCaloEmEz && hCaloEmEz->getRootObject()) hCaloEmEz->Fill(_EmCaloEz);
01191 hCaloEmMET = _dbe->get(DirName+"/"+"METTask_CaloEmMET"); if (hCaloEmMET && hCaloEmMET->getRootObject()) hCaloEmMET->Fill(_EmMET);
01192 hCaloEmMETPhi = _dbe->get(DirName+"/"+"METTask_CaloEmMETPhi"); if (hCaloEmMETPhi && hCaloEmMETPhi->getRootObject()) hCaloEmMETPhi->Fill(_EmMetPhi);
01193 hCaloEmSumET = _dbe->get(DirName+"/"+"METTask_CaloEmSumET"); if (hCaloEmSumET && hCaloEmSumET->getRootObject()) hCaloEmSumET->Fill(_EmSumEt);
01194
01195 hCaloHaMEx = _dbe->get(DirName+"/"+"METTask_CaloHaMEx"); if (hCaloHaMEx && hCaloHaMEx->getRootObject()) hCaloHaMEx->Fill(_HaMEx);
01196 hCaloHaMEy = _dbe->get(DirName+"/"+"METTask_CaloHaMEy"); if (hCaloHaMEy && hCaloHaMEy->getRootObject()) hCaloHaMEy->Fill(_HaMEy);
01197 hCaloHaEz = _dbe->get(DirName+"/"+"METTask_CaloHaEz"); if (hCaloHaEz && hCaloHaEz->getRootObject()) hCaloHaEz->Fill(_HaCaloEz);
01198 hCaloHaMET = _dbe->get(DirName+"/"+"METTask_CaloHaMET"); if (hCaloHaMET && hCaloHaMET->getRootObject()) hCaloHaMET->Fill(_HaMET);
01199 hCaloHaMETPhi = _dbe->get(DirName+"/"+"METTask_CaloHaMETPhi"); if (hCaloHaMETPhi && hCaloHaMETPhi->getRootObject()) hCaloHaMETPhi->Fill(_HaMetPhi);
01200 hCaloHaSumET = _dbe->get(DirName+"/"+"METTask_CaloHaSumET"); if (hCaloHaSumET && hCaloHaSumET->getRootObject()) hCaloHaSumET->Fill(_HaSumEt);
01201
01202 }
01203
01204 if (theCaloMETCollectionLabel.label() == "corMetGlobalMuons" ) {
01205
01206 for( reco::MuonCollection::const_iterator muonit = muon_h->begin(); muonit != muon_h->end(); muonit++ ) {
01207 const reco::TrackRef siTrack = muonit->innerTrack();
01208 hCalomuPt = _dbe->get(DirName+"/"+"METTask_CalomuPt"); if (hCalomuPt && hCalomuPt->getRootObject()) hCalomuPt->Fill( muonit->p4().pt() );
01209 hCalomuEta = _dbe->get(DirName+"/"+"METTask_CalomuEta"); if (hCalomuEta && hCalomuEta->getRootObject()) hCalomuEta->Fill( muonit->p4().eta() );
01210 hCalomuNhits = _dbe->get(DirName+"/"+"METTask_CalomuNhits"); if (hCalomuNhits && hCalomuNhits->getRootObject()) hCalomuNhits->Fill( siTrack.isNonnull() ? siTrack->numberOfValidHits() : -999 );
01211 hCalomuChi2 = _dbe->get(DirName+"/"+"METTask_CalomuChi2"); if (hCalomuChi2 && hCalomuChi2->getRootObject()) hCalomuChi2->Fill( siTrack.isNonnull() ? siTrack->chi2()/siTrack->ndof() : -999 );
01212 double d0 = siTrack.isNonnull() ? -1 * siTrack->dxy( bspot) : -999;
01213 hCalomuD0 = _dbe->get(DirName+"/"+"METTask_CalomuD0"); if (hCalomuD0 && hCalomuD0->getRootObject()) hCalomuD0->Fill( d0 );
01214 }
01215
01216 const unsigned int nMuons = muon_h->size();
01217 for( unsigned int mus = 0; mus < nMuons; mus++ ) {
01218 reco::MuonRef muref( muon_h, mus);
01219 reco::MuonMETCorrectionData muCorrData = (*corMetGlobalMuons_ValueMap_Handle)[muref];
01220 hCaloMExCorrection = _dbe->get(DirName+"/"+"METTask_CaloMExCorrection"); if (hCaloMExCorrection && hCaloMExCorrection->getRootObject()) hCaloMExCorrection-> Fill(muCorrData.corrY());
01221 hCaloMEyCorrection = _dbe->get(DirName+"/"+"METTask_CaloMEyCorrection"); if (hCaloMEyCorrection && hCaloMEyCorrection->getRootObject()) hCaloMEyCorrection-> Fill(muCorrData.corrX());
01222 hCaloMuonCorrectionFlag = _dbe->get(DirName+"/"+"METTask_CaloMuonCorrectionFlag"); if (hCaloMuonCorrectionFlag && hCaloMuonCorrectionFlag->getRootObject()) hCaloMuonCorrectionFlag-> Fill(muCorrData.type());
01223 }
01224 }
01225 }
01226
01227 }
01228
01229
01230 bool CaloMETAnalyzer::selectHighPtJetEvent(const edm::Event& iEvent){
01231
01232 bool return_value=false;
01233
01234 edm::Handle<reco::CaloJetCollection> caloJets;
01235 iEvent.getByLabel(theJetCollectionLabel, caloJets);
01236 if (!caloJets.isValid()) {
01237 LogDebug("") << "CaloMETAnalyzer: Could not find jet product" << std::endl;
01238 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find jet product" << std::endl;
01239 }
01240
01241 for (reco::CaloJetCollection::const_iterator cal = caloJets->begin();
01242 cal!=caloJets->end(); ++cal){
01243 if (cal->pt()>_highPtJetThreshold){
01244 return_value=true;
01245 }
01246 }
01247
01248 return return_value;
01249
01250 }
01251
01252
01253 bool CaloMETAnalyzer::selectLowPtJetEvent(const edm::Event& iEvent){
01254
01255 bool return_value=false;
01256
01257 edm::Handle<reco::CaloJetCollection> caloJets;
01258 iEvent.getByLabel(theJetCollectionLabel, caloJets);
01259 if (!caloJets.isValid()) {
01260 LogDebug("") << "CaloMETAnalyzer: Could not find jet product" << std::endl;
01261 if (_verbose) std::cout << "CaloMETAnalyzer: Could not find jet product" << std::endl;
01262 }
01263
01264 for (reco::CaloJetCollection::const_iterator cal = caloJets->begin();
01265 cal!=caloJets->end(); ++cal){
01266 if (cal->pt()>_lowPtJetThreshold){
01267 return_value=true;
01268 }
01269 }
01270
01271 return return_value;
01272
01273 }
01274
01275
01276 bool CaloMETAnalyzer::selectWElectronEvent(const edm::Event& iEvent){
01277
01278 bool return_value=true;
01279
01280
01281
01282
01283
01284 return return_value;
01285
01286 }
01287
01288
01289 bool CaloMETAnalyzer::selectWMuonEvent(const edm::Event& iEvent){
01290
01291 bool return_value=true;
01292
01293
01294
01295
01296
01297 return return_value;
01298
01299 }
01300