ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/csander/MyAnalysis/AnalysisSkeleton/src/AnalysisSkeleton.cc
Revision: 1.3
Committed: Tue Sep 11 11:31:27 2012 UTC (12 years, 7 months ago) by csander
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +311 -214 lines
Error occurred while calculating annotation data.
Log Message:
Updated to 52X for dead ECAL studies

File Contents

# Content
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);