1 |
// -*- C++ -*-
|
2 |
//
|
3 |
// Package: AnalysisSkeleton
|
4 |
// Class: AnalysisSkeleton
|
5 |
//
|
6 |
/**\class AnalysisSkeleton AnalysisSkeleton.cc MyAnalysis/AnalysisSkeleton/src/AnalysisSkeleton.cc
|
7 |
|
8 |
Description: [one line class summary]
|
9 |
|
10 |
Implementation:
|
11 |
[Notes on implementation]
|
12 |
*/
|
13 |
//
|
14 |
// Original Author: Christian Sander,,,
|
15 |
// Created: Wed Jun 30 13:55:44 CEST 2010
|
16 |
// $Id: AnalysisSkeleton.cc,v 1.2 2012/04/17 13:25:03 csander Exp $
|
17 |
//
|
18 |
//
|
19 |
/// system include files
|
20 |
#include <memory>
|
21 |
#include <string>
|
22 |
#include <vector>
|
23 |
#include <map>
|
24 |
//#include <cassert>
|
25 |
#include <cmath>
|
26 |
#include <iostream>
|
27 |
|
28 |
// user include files
|
29 |
#include "FWCore/Framework/interface/Frameworkfwd.h"
|
30 |
#include "FWCore/Framework/interface/EDAnalyzer.h"
|
31 |
#include "FWCore/Framework/interface/Event.h"
|
32 |
#include "FWCore/Framework/interface/EventSetup.h"
|
33 |
#include "FWCore/Framework/interface/ESHandle.h"
|
34 |
#include "FWCore/Framework/interface/MakerMacros.h"
|
35 |
#include "FWCore/ParameterSet/interface/ParameterSet.h"
|
36 |
#include "FWCore/ServiceRegistry/interface/Service.h"
|
37 |
#include "CommonTools/UtilAlgos/interface/TFileService.h"
|
38 |
#include "FWCore/Utilities/interface/EDMException.h"
|
39 |
#include "FWCore/MessageLogger/interface/MessageLogger.h"
|
40 |
|
41 |
#include "DataFormats/Math/interface/LorentzVector.h"
|
42 |
#include "DataFormats/Math/interface/deltaR.h"
|
43 |
#include "DataFormats/Math/interface/deltaPhi.h"
|
44 |
#include "DataFormats/JetReco/interface/Jet.h"
|
45 |
#include "DataFormats/JetReco/interface/CaloJet.h"
|
46 |
#include "DataFormats/JetReco/interface/GenJet.h"
|
47 |
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
|
48 |
#include "DataFormats/PatCandidates/interface/Jet.h"
|
49 |
#include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
|
50 |
#include "DataFormats/BTauReco/interface/JetTag.h"
|
51 |
|
52 |
#include "TH1F.h"
|
53 |
#include "TH2F.h"
|
54 |
#include "TMath.h"
|
55 |
#include "TFile.h"
|
56 |
|
57 |
// Ecal
|
58 |
#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
|
59 |
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
|
60 |
#include "DataFormats/DetId/interface/DetId.h"
|
61 |
#include "DataFormats/EcalDetId/interface/EBDetId.h"
|
62 |
#include "DataFormats/EcalDetId/interface/EEDetId.h"
|
63 |
|
64 |
#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
|
65 |
#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
|
66 |
|
67 |
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
|
68 |
#include "CalibCalorimetry/EcalTPGTools/interface/EcalTPGScale.h"
|
69 |
#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
|
70 |
#include "Geometry/Records/interface/IdealGeometryRecord.h"
|
71 |
|
72 |
#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
|
73 |
#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
|
74 |
#include "Geometry/CaloTopology/interface/CaloTopology.h"
|
75 |
|
76 |
// Geometry
|
77 |
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
|
78 |
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
|
79 |
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
|
80 |
#include "Geometry/Records/interface/CaloGeometryRecord.h"
|
81 |
|
82 |
#include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
|
83 |
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
|
84 |
|
85 |
// b-Tagging
|
86 |
#include "DataFormats/BTauReco/interface/JetTag.h"
|
87 |
|
88 |
//
|
89 |
// class declaration
|
90 |
//
|
91 |
using namespace edm;
|
92 |
using namespace std;
|
93 |
using namespace reco;
|
94 |
|
95 |
class AnalysisSkeleton: public edm::EDAnalyzer {
|
96 |
public:
|
97 |
explicit AnalysisSkeleton(const edm::ParameterSet&);
|
98 |
~AnalysisSkeleton();
|
99 |
bool goodJet(const pat::Jet* jet);
|
100 |
|
101 |
private:
|
102 |
virtual void beginJob();
|
103 |
virtual void analyze(const edm::Event&, const edm::EventSetup&);
|
104 |
virtual void endJob();
|
105 |
|
106 |
//
|
107 |
// constants, enums and typedefs
|
108 |
//
|
109 |
|
110 |
//
|
111 |
// static data member definitions
|
112 |
//
|
113 |
// ----------member data ---------------------------
|
114 |
edm::InputTag _jetTag;
|
115 |
edm::InputTag _genJetTag;
|
116 |
edm::InputTag _weightName;
|
117 |
edm::InputTag _EBRecHits;
|
118 |
edm::InputTag _EERecHits;
|
119 |
|
120 |
double _deltaPhiDiJet;
|
121 |
double _absCut3rdJet;
|
122 |
double _relCut3rdJet;
|
123 |
double _deltaRDeadECal;
|
124 |
vector<int> _maskedEcalChannelStatus;
|
125 |
std::string _bTagName;
|
126 |
double _bTagThreshold;
|
127 |
|
128 |
double weight;
|
129 |
|
130 |
// Channel status related
|
131 |
virtual void envSet(const edm::EventSetup&);
|
132 |
edm::ESHandle<EcalChannelStatus> ecalStatus; // these come from EventSetup
|
133 |
edm::ESHandle<CaloGeometry> geometry;
|
134 |
// Store DetId <==> vector<double> (eta, phi, theta)
|
135 |
std::map<DetId, std::vector<double> > EcalAllDeadChannelsValMap;
|
136 |
int getChannelStatusMaps();
|
137 |
bool mapsReady;
|
138 |
EcalTPGScale ecalScale;
|
139 |
|
140 |
|
141 |
TH1F *h_jetpt_1, *h_jetphi_1, *h_jeteta_1;
|
142 |
TH1F *h_jetpt_2, *h_jetphi_2, *h_jeteta_2;
|
143 |
TH1F *h_jetpt_3, *h_jetphi_3, *h_jeteta_3;
|
144 |
|
145 |
TH2F *h_diJetAsym_tot;
|
146 |
TH2F *h_diJetAsym_nob_nodead;
|
147 |
TH2F *h_diJetAsym_nob_dead;
|
148 |
TH2F *h_diJetAsym_b_nodead;
|
149 |
TH2F *h_diJetAsym_b_dead;
|
150 |
|
151 |
};
|
152 |
|
153 |
//
|
154 |
// constructors and destructor
|
155 |
//
|
156 |
AnalysisSkeleton::AnalysisSkeleton(const edm::ParameterSet& iConfig)
|
157 |
|
158 |
{
|
159 |
//now do what ever initialization is needed
|
160 |
_jetTag = iConfig.getParameter<edm::InputTag> ("jetTag");
|
161 |
_weightName = iConfig.getParameter<edm::InputTag> ("weightName");
|
162 |
_EBRecHits = iConfig.getParameter<edm::InputTag> ("EBRecHits");
|
163 |
_EERecHits = iConfig.getParameter<edm::InputTag> ("EERecHits");
|
164 |
_deltaPhiDiJet = iConfig.getParameter<double> ("deltaPhiDiJet");
|
165 |
_absCut3rdJet = iConfig.getParameter<double> ("absCut3rdJet");
|
166 |
_relCut3rdJet = iConfig.getParameter<double> ("relCut3rdJet");
|
167 |
_deltaRDeadECal = iConfig.getParameter<double> ("deltaRDeadECal");
|
168 |
_maskedEcalChannelStatus = iConfig.getParameter<vector<int> > ("maskedEcalChannelStatus");
|
169 |
_bTagName = iConfig.getParameter<string> ("bTagName");
|
170 |
_bTagThreshold = iConfig.getParameter<double> ("bTagThreshold");
|
171 |
}
|
172 |
|
173 |
AnalysisSkeleton::~AnalysisSkeleton() {
|
174 |
|
175 |
// do anything here that needs to be done at destruction time
|
176 |
// (e.g. close files, deallocate resources etc.)
|
177 |
|
178 |
}
|
179 |
|
180 |
//
|
181 |
// member functions
|
182 |
//
|
183 |
|
184 |
// ------------ method called to for each event ------------
|
185 |
void AnalysisSkeleton::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
|
186 |
|
187 |
// Event setup
|
188 |
envSet(iSetup);
|
189 |
|
190 |
getChannelStatusMaps();
|
191 |
|
192 |
//Weight
|
193 |
edm::Handle<double> event_weight;
|
194 |
bool findWeight = iEvent.getByLabel(_weightName, event_weight);
|
195 |
weight = (event_weight.isValid() ? (*event_weight) : 1.0);
|
196 |
|
197 |
//RecoJets
|
198 |
edm::Handle<edm::View<pat::Jet> > Jets_rec;
|
199 |
iEvent.getByLabel(_jetTag, Jets_rec);
|
200 |
|
201 |
//// Fill the histograms
|
202 |
|
203 |
bool pointsToDead = false;
|
204 |
bool hasBTag = false;
|
205 |
const pat::Jet *jet1, *jet2, *jet3;
|
206 |
int NJet = 0;
|
207 |
for (edm::View<pat::Jet>::const_iterator it = Jets_rec->begin(); it != Jets_rec->end(); ++it) {
|
208 |
if (goodJet(&(*it))) {
|
209 |
++NJet;
|
210 |
//// look if one of two leading jets points into direction of masked ECAL cluster
|
211 |
if (NJet < 3){
|
212 |
for (std::map<DetId, std::vector<double> >::iterator kt = EcalAllDeadChannelsValMap.begin(); kt != EcalAllDeadChannelsValMap.end(); ++kt) {
|
213 |
math::PtEtaPhiMLorentzVectorD Evec(100., kt->second.at(0), kt->second.at(1), 0.);
|
214 |
double dR_dead = deltaR(*it, Evec);
|
215 |
if (dR_dead < _deltaRDeadECal) {
|
216 |
pointsToDead = true;
|
217 |
}
|
218 |
}
|
219 |
}
|
220 |
//// look if one of two leading jets has b-tag
|
221 |
if (NJet < 3){
|
222 |
//cout << it->bDiscriminator(_bTagName.c_str()) << endl;
|
223 |
if (it->bDiscriminator(_bTagName.c_str()) > _bTagThreshold) {
|
224 |
hasBTag = true;
|
225 |
}
|
226 |
}
|
227 |
if (NJet == 1){
|
228 |
jet1 = &(*it);
|
229 |
h_jetpt_1->Fill(it->pt());
|
230 |
h_jeteta_1->Fill(it->eta());
|
231 |
h_jetphi_1->Fill(it->phi());
|
232 |
}
|
233 |
if (NJet == 2){
|
234 |
jet2 = &(*it);
|
235 |
h_jetpt_2->Fill(it->pt());
|
236 |
h_jeteta_2->Fill(it->eta());
|
237 |
h_jetphi_2->Fill(it->phi());
|
238 |
}
|
239 |
if (NJet == 3){
|
240 |
jet3 = &(*it);
|
241 |
h_jetpt_3->Fill(it->pt());
|
242 |
h_jeteta_3->Fill(it->eta());
|
243 |
h_jetphi_3->Fill(it->phi());
|
244 |
}
|
245 |
}
|
246 |
}
|
247 |
|
248 |
// cout << "Jet1 (pt, eta, phi): " << jet1->pt() << ", " << jet1->eta() << ", " << jet1->phi() << endl;
|
249 |
// cout << "Jet2 (pt, eta, phi): " << jet2->pt() << ", " << jet2->eta() << ", " << jet2->phi() << endl;
|
250 |
// cout << "Jet3 (pt, eta, phi): " << jet3->pt() << ", " << jet3->eta() << ", " << jet3->phi() << endl;
|
251 |
|
252 |
double ptAve = (jet1->pt() + jet2->pt()) / 2.;
|
253 |
double deltaPhi12 = deltaPhi(jet1->phi(), jet2->phi());
|
254 |
double alpha = (jet3->pt()/ptAve);
|
255 |
double ptAsym = std::abs(jet1->pt() - jet2->pt())/ ptAve;
|
256 |
if (deltaPhi12 > _deltaPhiDiJet && jet3->pt() < _absCut3rdJet && alpha < _relCut3rdJet){
|
257 |
h_diJetAsym_tot->Fill(ptAsym, ptAve);
|
258 |
if (hasBTag) {
|
259 |
if (pointsToDead){
|
260 |
h_diJetAsym_b_dead->Fill(ptAsym, ptAve);
|
261 |
}
|
262 |
else {
|
263 |
h_diJetAsym_b_nodead->Fill(ptAsym, ptAve);
|
264 |
}
|
265 |
} else {
|
266 |
if (pointsToDead){
|
267 |
h_diJetAsym_nob_dead->Fill(ptAsym, ptAve);
|
268 |
}
|
269 |
else {
|
270 |
h_diJetAsym_nob_nodead->Fill(ptAsym, ptAve);
|
271 |
}
|
272 |
}
|
273 |
}
|
274 |
|
275 |
}
|
276 |
|
277 |
// ------------ method called once each job just before starting event loop ------------
|
278 |
void AnalysisSkeleton::beginJob() {
|
279 |
edm::Service < TFileService > fs;
|
280 |
if (!fs) {
|
281 |
throw edm::Exception(edm::errors::Configuration, "TFile Service is not registered in cfg file");
|
282 |
}
|
283 |
h_jetpt_1 = fs->make < TH1F > ("jetpt_1", "Jet1 Pt", 100, 0., 1000.);
|
284 |
h_jetpt_1->Sumw2();
|
285 |
h_jeteta_1 = fs->make < TH1F > ("jeteta_1", "Jet1 Eta", 100, -5., 5.);
|
286 |
h_jeteta_1->Sumw2();
|
287 |
h_jetphi_1 = fs->make < TH1F > ("jetphi_1", "Jet1 Phi", 100, -TMath::Pi(), -TMath::Pi());
|
288 |
h_jetphi_1->Sumw2();
|
289 |
h_jetpt_2 = fs->make < TH1F > ("jetpt_2", "Jet2 Pt", 100, 0., 1000.);
|
290 |
h_jetpt_2->Sumw2();
|
291 |
h_jeteta_2 = fs->make < TH1F > ("jeteta_2", "Jet2 Eta", 100, -5., 5.);
|
292 |
h_jeteta_2->Sumw2();
|
293 |
h_jetphi_2 = fs->make < TH1F > ("jetphi_2", "Jet2 Phi", 100, -TMath::Pi(), -TMath::Pi());
|
294 |
h_jetphi_2->Sumw2();
|
295 |
h_jetpt_3 = fs->make < TH1F > ("jetpt_3", "Jet3 Pt", 100, 0., 1000.);
|
296 |
h_jetpt_3->Sumw2();
|
297 |
h_jeteta_3 = fs->make < TH1F > ("jeteta_3", "Jet3 Eta", 100, -5., 5.);
|
298 |
h_jeteta_3->Sumw2();
|
299 |
h_jetphi_3 = fs->make < TH1F > ("jetphi_3", "Jet3 Phi", 100, -TMath::Pi(), -TMath::Pi());
|
300 |
h_jetphi_3->Sumw2();
|
301 |
|
302 |
h_diJetAsym_tot = fs->make < TH2F > ("diJetAsym_tot", "diJetAsym_tot", 50, 0., 2.5, 20., 0, 2000.);
|
303 |
h_diJetAsym_nob_nodead = fs->make < TH2F > ("diJetAsym_nob_nodead", "diJetAsym_nob_nodead", 50, 0., 2.5, 20., 0, 2000.);
|
304 |
h_diJetAsym_nob_dead = fs->make < TH2F > ("diJetAsym_nob_dead", "diJetAsym_nob_dead", 50, 0., 2.5, 20., 0, 2000.);
|
305 |
h_diJetAsym_b_nodead = fs->make < TH2F > ("diJetAsym_b_nodead", "diJetAsym_b_nodead", 50, 0., 2.5, 20., 0, 2000.);
|
306 |
h_diJetAsym_b_dead = fs->make < TH2F > ("diJetAsym_b_dead", "diJetAsym_b_dead", 50, 0., 2.5, 20., 0, 2000.);
|
307 |
|
308 |
mapsReady = false;
|
309 |
|
310 |
}
|
311 |
|
312 |
// ------------ method called once each job just after ending the event loop ------------
|
313 |
void AnalysisSkeleton::endJob() {
|
314 |
}
|
315 |
|
316 |
bool AnalysisSkeleton::goodJet(const pat::Jet* jet) {
|
317 |
PFJetIDSelectionFunctor jetIDLoosePF(PFJetIDSelectionFunctor::FIRSTDATA, PFJetIDSelectionFunctor::LOOSE);
|
318 |
pat::strbitset ret = jetIDLoosePF.getBitTemplate();
|
319 |
ret.set(false);
|
320 |
bool loose = jetIDLoosePF(*jet, ret);
|
321 |
if (!loose)
|
322 |
return false;
|
323 |
|
324 |
return true;
|
325 |
}
|
326 |
|
327 |
void AnalysisSkeleton::envSet(const edm::EventSetup& iSetup) {
|
328 |
|
329 |
ecalScale.setEventSetup(iSetup);
|
330 |
|
331 |
iSetup.get<EcalChannelStatusRcd> ().get(ecalStatus);
|
332 |
iSetup.get<CaloGeometryRecord> ().get(geometry);
|
333 |
|
334 |
if (!ecalStatus.isValid())
|
335 |
throw "Failed to get ECAL channel status!";
|
336 |
if (!geometry.isValid())
|
337 |
throw "Failed to get the geometry!";
|
338 |
|
339 |
}
|
340 |
|
341 |
int AnalysisSkeleton::getChannelStatusMaps() {
|
342 |
|
343 |
if (mapsReady)
|
344 |
return -1;
|
345 |
|
346 |
EcalAllDeadChannelsValMap.clear();
|
347 |
|
348 |
// Loop over EB ...
|
349 |
for (int ieta = -85; ieta <= 85; ieta++) {
|
350 |
for (int iphi = 0; iphi <= 360; iphi++) {
|
351 |
if (!EBDetId::validDetId(ieta, iphi))
|
352 |
continue;
|
353 |
|
354 |
const EBDetId detid = EBDetId(ieta, iphi, EBDetId::ETAPHIMODE);
|
355 |
EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
|
356 |
// refer https://twiki.cern.ch/twiki/bin/viewauth/CMS/EcalChannelStatus
|
357 |
int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
|
358 |
//std::cout << ieta << " " << iphi << " " << status << std:: endl;
|
359 |
const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
|
360 |
const CaloCellGeometry* cellGeom = subGeom->getGeometry(detid);
|
361 |
double eta = cellGeom->getPosition().eta();
|
362 |
double phi = cellGeom->getPosition().phi();
|
363 |
double theta = cellGeom->getPosition().theta();
|
364 |
|
365 |
for (vector<int>::iterator it = _maskedEcalChannelStatus.begin(); it != _maskedEcalChannelStatus.end(); ++it){
|
366 |
if (status == *it) {
|
367 |
std::vector<double> valVec;
|
368 |
valVec.push_back(eta);
|
369 |
valVec.push_back(phi);
|
370 |
valVec.push_back(theta);
|
371 |
//std::cout << eta << " " << phi << std::endl;
|
372 |
EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
|
373 |
}
|
374 |
}
|
375 |
} // end loop iphi
|
376 |
} // end loop ieta
|
377 |
|
378 |
// Loop over EE detid
|
379 |
for (int ix = 0; ix <= 100; ix++) {
|
380 |
for (int iy = 0; iy <= 100; iy++) {
|
381 |
for (int iz = -1; iz <= 1; iz++) {
|
382 |
if (iz == 0)
|
383 |
continue;
|
384 |
if (!EEDetId::validDetId(ix, iy, iz))
|
385 |
continue;
|
386 |
|
387 |
const EEDetId detid = EEDetId(ix, iy, iz, EEDetId::XYMODE);
|
388 |
EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
|
389 |
int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
|
390 |
|
391 |
const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
|
392 |
const CaloCellGeometry* cellGeom = subGeom->getGeometry(detid);
|
393 |
double eta = cellGeom->getPosition().eta();
|
394 |
double phi = cellGeom->getPosition().phi();
|
395 |
double theta = cellGeom->getPosition().theta();
|
396 |
|
397 |
for (vector<int>::iterator it = _maskedEcalChannelStatus.begin(); it != _maskedEcalChannelStatus.end(); ++it){
|
398 |
if (status == *it) {
|
399 |
std::vector<double> valVec;
|
400 |
valVec.push_back(eta);
|
401 |
valVec.push_back(phi);
|
402 |
valVec.push_back(theta);
|
403 |
//std::cout << eta << " " << phi << std::endl;
|
404 |
EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
|
405 |
}
|
406 |
}
|
407 |
} // end loop iz
|
408 |
} // end loop iy
|
409 |
} // end loop ix
|
410 |
|
411 |
mapsReady = true;
|
412 |
return 1;
|
413 |
}
|
414 |
|
415 |
|
416 |
//define this as a plug-in
|
417 |
DEFINE_FWK_MODULE (AnalysisSkeleton);
|