ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/RecHitTreeProducer.cc
Revision: 1.20
Committed: Wed Jun 6 15:53:41 2012 UTC (12 years, 11 months ago) by yjlee
Content type: text/plain
Branch: MAIN
CVS Tags: HiForest_V02_34, HiForest_V02_33, HiForest_V02_32, HiForest_V02_31, HiForest_V02_30, HiForest_V02_27, HiForest_V02_26, QM_2012, HiForest_V02_25, HiForest_V02_24, HiForest_V02_23, HiForest_V02_22, HiForest_V02_21, HiForest_V02_20, HiForest_V02_19
Branch point for: branch_44x
Changes since 1.19: +2 -2 lines
Log Message: 
Lower the array size

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: RecHitTreeProducer
4 // Class: RecHitTreeProducer
5 //
6 /**\class RecHitTreeProducer RecHitTreeProducer.cc CmsHi/RecHitTreeProducer/src/RecHitTreeProducer.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Yetkin Yilmaz
15 // Modified: Frank Ma, Yen-Jie Lee
16 // Created: Tue Sep 7 11:38:19 EDT 2010
17 // $Id: RecHitTreeProducer.cc,v 1.19 2011/11/30 21:26:08 yilmaz Exp $
18 //
19 //
20
21 #define versionTag "v1"
22 // system include files
23 #include <memory>
24 #include <vector>
25 #include <iostream>
26
27 // user include files
28 #include "FWCore/Framework/interface/Frameworkfwd.h"
29 #include "FWCore/Framework/interface/EDAnalyzer.h"
30
31 #include "FWCore/Framework/interface/Event.h"
32 #include "FWCore/Framework/interface/MakerMacros.h"
33 #include "FWCore/Framework/interface/ESHandle.h"
34 #include "FWCore/ParameterSet/interface/ParameterSet.h"
35
36 #include "DataFormats/DetId/interface/DetId.h"
37 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
38 #include "Geometry/Records/interface/IdealGeometryRecord.h"
39 #include "Geometry/Records/interface/CaloGeometryRecord.h"
40 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
41
42 #include "CommonTools/UtilAlgos/interface/TFileService.h"
43 #include "FWCore/ServiceRegistry/interface/Service.h"
44
45 #include "DataFormats/Math/interface/deltaR.h"
46
47 #include "DataFormats/Common/interface/Handle.h"
48 #include "DataFormats/FWLite/interface/ChainEvent.h"
49 #include "FWCore/Utilities/interface/InputTag.h"
50 #include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
51 #include "DataFormats/CaloTowers/interface/CaloTower.h"
52 #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
53 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
54 #include "DataFormats/PatCandidates/interface/Jet.h"
55 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
56
57 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
58 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
59 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
60 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
61 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
62 #include "DataFormats/DetId/interface/DetId.h"
63
64 #include "DataFormats/VertexReco/interface/Vertex.h"
65 #include "DataFormats/VertexReco/interface/VertexFwd.h"
66
67
68 #include "TNtuple.h"
69
70 using namespace std;
71
72 #define MAXHITS 100000
73
74 struct MyRecHit{
75 int depth[MAXHITS];
76 int n;
77
78 float e[MAXHITS];
79 float et[MAXHITS];
80 float eta[MAXHITS];
81 float phi[MAXHITS];
82 float perp[MAXHITS];
83 float emEt[MAXHITS];
84 float hadEt[MAXHITS];
85
86 bool isjet[MAXHITS];
87 float etVtx[MAXHITS];
88 float etaVtx[MAXHITS];
89 float phiVtx[MAXHITS];
90 float perpVtx[MAXHITS];
91 float emEtVtx[MAXHITS];
92 float hadEtVtx[MAXHITS];
93
94 float jtpt;
95 float jteta;
96 float jtphi;
97
98 };
99
100
101 struct MyBkg{
102 int n;
103 float rho[50];
104 float sigma[50];
105 };
106
107
108 //
109 // class declaration
110 //
111
112 class RecHitTreeProducer : public edm::EDAnalyzer {
113 public:
114 explicit RecHitTreeProducer(const edm::ParameterSet&);
115 ~RecHitTreeProducer();
116 double getEt(const DetId &id, double energy, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
117 double getEta(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
118 double getPhi(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
119 double getPerp(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
120
121 math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
122 double getEt(math::XYZPoint pos, double energy);
123
124 reco::Vertex::Point getVtx(const edm::Event& ev);
125
126
127
128 private:
129 virtual void beginJob() ;
130 virtual void analyze(const edm::Event&, const edm::EventSetup&);
131 virtual void endJob() ;
132
133 // ----------member data ---------------------------
134
135 edm::Handle<reco::Centrality> cent;
136 edm::Handle<vector<double> > ktRhos;
137 edm::Handle<vector<double> > akRhos;
138
139 edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
140 edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
141
142 edm::Handle<HFRecHitCollection> hfHits;
143 edm::Handle<HBHERecHitCollection> hbheHits;
144
145 edm::Handle<reco::BasicClusterCollection> bClusters;
146 edm::Handle<CaloTowerCollection> towers;
147 edm::Handle<reco::VertexCollection> vtxs;
148
149 typedef vector<EcalRecHit>::const_iterator EcalIterator;
150
151 edm::Handle<reco::CaloJetCollection> jets;
152
153 edm::Handle<std::vector<double> > rhos;
154 edm::Handle<std::vector<double> > sigmas;
155
156 MyRecHit hbheRecHit;
157 MyRecHit hfRecHit;
158 MyRecHit ebRecHit;
159 MyRecHit eeRecHit;
160 MyRecHit myBC;
161 MyRecHit myTowers;
162 MyBkg bkg;
163
164 TNtuple* nt;
165 TTree* hbheTree;
166 TTree* hfTree;
167 TTree* ebTree;
168 TTree* eeTree;
169 TTree* bcTree;
170 TTree* towerTree;
171 TTree* bkgTree;
172
173 double cone;
174 double hfTowerThreshold_;
175 double hfLongThreshold_;
176 double hfShortThreshold_;
177 double hbheThreshold_;
178 double ebThreshold_;
179 double eeThreshold_;
180
181 double hbhePtMin_;
182 double hfPtMin_;
183 double ebPtMin_;
184 double eePtMin_;
185 double towerPtMin_;
186
187 edm::Service<TFileService> fs;
188 const CentralityBins * cbins_;
189 const CaloGeometry *geo;
190
191 edm::InputTag HcalRecHitHFSrc_;
192 edm::InputTag HcalRecHitHBHESrc_;
193 edm::InputTag EBSrc_;
194 edm::InputTag EESrc_;
195 edm::InputTag BCSrc_;
196 edm::InputTag TowerSrc_;
197 edm::InputTag VtxSrc_;
198
199 edm::InputTag JetSrc_;
200
201 edm::InputTag FastJetTag_;
202
203 bool useJets_;
204 bool doBasicClusters_;
205 bool doTowers_;
206 bool doEcal_;
207 bool doHcal_;
208 bool doHF_;
209 bool hasVtx_;
210 bool saveBothVtx_;
211
212 bool doFastJet_;
213
214 bool doEbyEonly_;
215
216 };
217
218 //
219 // constants, enums and typedefs
220 //
221
222 //
223 // static data member definitions
224 //
225
226 //
227 // constructors and destructor
228 //
229 RecHitTreeProducer::RecHitTreeProducer(const edm::ParameterSet& iConfig) :
230 cbins_(0),
231 geo(0),
232 cone(0.5)
233 {
234 //now do what ever initialization is needed
235 EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
236 EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
237 HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
238 HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
239 BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
240 TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
241 VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
242 JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
243 useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
244 doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
245 doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
246 doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
247 doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
248 doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true);
249 hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",true);
250 saveBothVtx_ = iConfig.getUntrackedParameter<bool>("saveBothVtx",false);
251
252 doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
253 FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
254 doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
255 hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
256 hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
257 hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
258 hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
259 hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
260 ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
261 eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
262 towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
263
264
265 }
266
267
268 RecHitTreeProducer::~RecHitTreeProducer()
269 {
270
271 // do anything here that needs to be done at desctruction time
272 // (e.g. close files, deallocate resources etc.)
273
274 }
275
276
277 //
278 // member functions
279 //
280
281 // ------------ method called to for each event ------------
282 void
283 RecHitTreeProducer::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
284 {
285
286 hfRecHit.n = 0;
287 hbheRecHit.n = 0;
288 ebRecHit.n = 0;
289 eeRecHit.n = 0;
290 myBC.n = 0;
291 myTowers.n = 0;
292 bkg.n = 0;
293
294 // get vertex
295 reco::Vertex::Point vtx(0,0,0);
296 if (hasVtx_) vtx = getVtx(ev);
297
298 if(doEcal_){
299 ev.getByLabel(EBSrc_,ebHits);
300 ev.getByLabel(EESrc_,eeHits);
301 }
302
303 if(doHcal_){
304 ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
305 }
306 if(doHF_){
307 ev.getByLabel(HcalRecHitHFSrc_,hfHits);
308 }
309
310 if(useJets_) {
311 ev.getByLabel(JetSrc_,jets);
312 }
313
314 if(doBasicClusters_){
315 ev.getByLabel(BCSrc_,bClusters);
316 }
317
318 if(doTowers_){
319 ev.getByLabel(TowerSrc_,towers);
320 }
321
322 if(doFastJet_){
323 ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
324 ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
325 bkg.n = rhos->size();
326 for(unsigned int i = 0; i < rhos->size(); ++i){
327 bkg.rho[i] = (*rhos)[i];
328 bkg.sigma[i] = (*sigmas)[i];
329 }
330 }
331
332 if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
333
334 if(!geo){
335 edm::ESHandle<CaloGeometry> pGeo;
336 iSetup.get<CaloGeometryRecord>().get(pGeo);
337 geo = pGeo.product();
338 }
339
340 int nHFlongPlus = 0;
341 int nHFshortPlus = 0;
342 int nHFtowerPlus = 0;
343 int nHFlongMinus = 0;
344 int nHFshortMinus = 0;
345 int nHFtowerMinus = 0;
346
347 if(doHF_){
348 for(unsigned int i = 0; i < hfHits->size(); ++i){
349 const HFRecHit & hit= (*hfHits)[i];
350 hfRecHit.e[hfRecHit.n] = hit.energy();
351 math::XYZPoint pos = getPosition(hit.id(),vtx);
352
353 if(!saveBothVtx_){
354 hfRecHit.et[hfRecHit.n] = getEt(pos,hit.energy());
355 hfRecHit.eta[hfRecHit.n] = pos.eta();
356 hfRecHit.phi[hfRecHit.n] = pos.phi();
357 hfRecHit.perp[hfRecHit.n] = pos.rho();
358 }else{
359 hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
360 hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
361 hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
362 hfRecHit.perp[hfRecHit.n] = getPerp(hit.id());
363
364 hfRecHit.etVtx[hfRecHit.n] = getEt(pos,hit.energy());
365 hfRecHit.etaVtx[hfRecHit.n] = pos.eta();
366 hfRecHit.phiVtx[hfRecHit.n] = pos.phi();
367 hfRecHit.perpVtx[hfRecHit.n] = pos.rho();
368
369 }
370
371 hfRecHit.isjet[hfRecHit.n] = false;
372 hfRecHit.depth[hfRecHit.n] = hit.id().depth();
373
374 if(hit.id().ieta() > 0){
375 if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
376 if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
377 }else{
378 if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
379 if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
380 }
381
382 if(useJets_){
383 for(unsigned int j = 0 ; j < jets->size(); ++j){
384 const reco::Jet& jet = (*jets)[j];
385 double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
386 if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
387 }
388 }
389 if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
390 }
391
392 if(doHcal_ && !doEbyEonly_){
393 for(unsigned int i = 0; i < hbheHits->size(); ++i){
394 const HBHERecHit & hit= (*hbheHits)[i];
395 if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
396
397 hbheRecHit.e[hbheRecHit.n] = hit.energy();
398 math::XYZPoint pos = getPosition(hit.id(),vtx);
399
400 if(!saveBothVtx_){
401 hbheRecHit.et[hbheRecHit.n] = getEt(pos,hit.energy());
402 hbheRecHit.eta[hbheRecHit.n] = pos.eta();
403 hbheRecHit.phi[hbheRecHit.n] = pos.phi();
404 hbheRecHit.perp[hbheRecHit.n] = pos.rho();
405 }else{
406 hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
407 hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
408 hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
409 hbheRecHit.perp[hbheRecHit.n] = getPerp(hit.id());
410
411 hbheRecHit.etVtx[hbheRecHit.n] = getEt(pos,hit.energy());
412 hbheRecHit.etaVtx[hbheRecHit.n] = pos.eta();
413 hbheRecHit.phiVtx[hbheRecHit.n] = pos.phi();
414 hbheRecHit.perpVtx[hbheRecHit.n] = pos.rho();
415
416 }
417
418 hbheRecHit.isjet[hbheRecHit.n] = false;
419 hbheRecHit.depth[hbheRecHit.n] = hit.id().depth();
420
421 if(useJets_){
422 for(unsigned int j = 0 ; j < jets->size(); ++j){
423 const reco::Jet& jet = (*jets)[j];
424 double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
425 if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
426 }
427 }
428 hbheRecHit.n++;
429 }
430 }
431 }
432 if(doEcal_ && !doEbyEonly_){
433 for(unsigned int i = 0; i < ebHits->size(); ++i){
434 const EcalRecHit & hit= (*ebHits)[i];
435 if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
436
437 ebRecHit.e[ebRecHit.n] = hit.energy();
438 math::XYZPoint pos = getPosition(hit.id(),vtx);
439
440 if(!saveBothVtx_){
441 ebRecHit.et[ebRecHit.n] = getEt(pos,hit.energy());
442 ebRecHit.eta[ebRecHit.n] = pos.eta();
443 ebRecHit.phi[ebRecHit.n] = pos.phi();
444 ebRecHit.perp[ebRecHit.n] = pos.rho();
445 }else{
446 ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
447 ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
448 ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
449 ebRecHit.perp[ebRecHit.n] = getPerp(hit.id());
450 ebRecHit.etVtx[ebRecHit.n] = getEt(pos,hit.energy());
451 ebRecHit.etaVtx[ebRecHit.n] = pos.eta();
452 ebRecHit.phiVtx[ebRecHit.n] = pos.phi();
453 ebRecHit.perpVtx[ebRecHit.n] = pos.rho();
454
455 }
456
457 ebRecHit.isjet[ebRecHit.n] = false;
458 if(useJets_){
459 for(unsigned int j = 0 ; j < jets->size(); ++j){
460 const reco::Jet& jet = (*jets)[j];
461 double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
462 if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
463 }
464 }
465 ebRecHit.n++;
466 }
467
468 for(unsigned int i = 0; i < eeHits->size(); ++i){
469 const EcalRecHit & hit= (*eeHits)[i];
470 if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
471
472 eeRecHit.e[eeRecHit.n] = hit.energy();
473 math::XYZPoint pos = getPosition(hit.id(),vtx);
474
475 if(!saveBothVtx_){
476 eeRecHit.et[eeRecHit.n] = getEt(pos,hit.energy());
477 eeRecHit.eta[eeRecHit.n] = pos.eta();
478 eeRecHit.phi[eeRecHit.n] = pos.phi();
479 eeRecHit.perp[eeRecHit.n] = pos.rho();
480 }else{
481 eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
482 eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
483 eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
484 eeRecHit.perp[eeRecHit.n] = getPerp(hit.id());
485 eeRecHit.etVtx[eeRecHit.n] = getEt(pos,hit.energy());
486 eeRecHit.etaVtx[eeRecHit.n] = pos.eta();
487 eeRecHit.phiVtx[eeRecHit.n] = pos.phi();
488 eeRecHit.perpVtx[eeRecHit.n] = pos.rho();
489
490 }
491
492 eeRecHit.isjet[eeRecHit.n] = false;
493
494 if(useJets_){
495 for(unsigned int j = 0 ; j < jets->size(); ++j){
496 const reco::Jet& jet = (*jets)[j];
497 double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
498 if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
499 }
500 }
501 eeRecHit.n++;
502 }
503 }
504
505 if(doTowers_){
506
507 for(unsigned int i = 0; i < towers->size(); ++i){
508 const CaloTower & hit= (*towers)[i];
509 if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
510
511 myTowers.et[myTowers.n] = hit.p4(vtx).Et();
512 myTowers.eta[myTowers.n] = hit.p4(vtx).Eta();
513 myTowers.phi[myTowers.n] = hit.p4(vtx).Phi();
514 myTowers.emEt[myTowers.n] = hit.emEt(vtx);
515 myTowers.hadEt[myTowers.n] = hit.hadEt(vtx);
516
517 if (saveBothVtx_) {
518 myTowers.e[myTowers.n] = hit.energy();
519 myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
520 myTowers.eta[myTowers.n] = getEta(hit.id());
521 myTowers.phi[myTowers.n] = getPhi(hit.id());
522 myTowers.isjet[myTowers.n] = false;
523 myTowers.etVtx[myTowers.n] = hit.p4(vtx).Et();
524 myTowers.etaVtx[myTowers.n] = hit.p4(vtx).Eta();
525 myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
526 myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
527 }
528
529 myTowers.isjet[myTowers.n] = false;
530
531 if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
532 if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
533
534 if(useJets_){
535 for(unsigned int j = 0 ; j < jets->size(); ++j){
536 const reco::Jet& jet = (*jets)[j];
537 double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
538 if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
539 }
540 }
541 myTowers.n++;
542 }
543
544 }
545
546 if(doBasicClusters_ && !doEbyEonly_){
547 for(unsigned int j = 0 ; j < jets->size(); ++j){
548 const reco::Jet& jet = (*jets)[j];
549 myBC.n = 0;
550 myBC.jtpt = jet.pt();
551 myBC.jteta = jet.eta();
552 myBC.jtphi = jet.phi();
553
554 for(unsigned int i = 0; i < bClusters->size(); ++i){
555 const reco::BasicCluster & bc= (*bClusters)[i];
556 double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
557 if(dr < cone){
558 myBC.e[myBC.n] = bc.energy();
559 myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
560 myBC.eta[myBC.n] = bc.eta();
561 myBC.phi[myBC.n] = bc.phi();
562 myBC.n++;
563 }
564 }
565 bcTree->Fill();
566 }
567 }
568
569 if(!doEbyEonly_){
570 towerTree->Fill();
571
572 eeTree->Fill();
573 ebTree->Fill();
574
575 hbheTree->Fill();
576 hfTree->Fill();
577
578 if (doFastJet_) {
579 bkgTree->Fill();
580 }
581 }
582
583 nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
584
585 }
586
587
588 // ------------ method called once each job just before starting event loop ------------
589 void
590 RecHitTreeProducer::beginJob()
591 {
592
593 hbheTree = fs->make<TTree>("hbhe",versionTag);
594 hbheTree->Branch("n",&hbheRecHit.n,"n/I");
595 hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
596 hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
597 hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
598 hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
599 hbheTree->Branch("perp",hbheRecHit.perp,"perp[n]/F");
600
601 hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
602 hbheTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
603
604 hfTree = fs->make<TTree>("hf",versionTag);
605 hfTree->Branch("n",&hfRecHit.n,"n/I");
606 hfTree->Branch("e",hfRecHit.e,"e[n]/F");
607 hfTree->Branch("et",hfRecHit.et,"et[n]/F");
608 hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
609 hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
610 hfTree->Branch("perp",hfRecHit.perp,"perp[n]/F");
611 hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
612 hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
613
614 eeTree = fs->make<TTree>("ee",versionTag);
615 eeTree->Branch("n",&eeRecHit.n,"n/I");
616 eeTree->Branch("e",eeRecHit.e,"e[n]/F");
617 eeTree->Branch("et",eeRecHit.et,"et[n]/F");
618 eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
619 eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
620 eeTree->Branch("perp",eeRecHit.perp,"perp[n]/F");
621
622 eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
623
624 ebTree = fs->make<TTree>("eb",versionTag);
625 ebTree->Branch("n",&ebRecHit.n,"n/I");
626 ebTree->Branch("e",ebRecHit.e,"e[n]/F");
627 ebTree->Branch("et",ebRecHit.et,"et[n]/F");
628 ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
629 ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
630 ebTree->Branch("perp",ebRecHit.perp,"perp[n]/F");
631
632 ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
633
634 towerTree = fs->make<TTree>("tower",versionTag);
635 towerTree->Branch("n",&myTowers.n,"n/I");
636 towerTree->Branch("e",myTowers.e,"e[n]/F");
637 towerTree->Branch("et",myTowers.et,"et[n]/F");
638 towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
639 towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
640 towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
641 towerTree->Branch("emEt",myTowers.emEt,"emEt[n]/F");
642 towerTree->Branch("hadEt",myTowers.hadEt,"hadEt[n]/F");
643
644 if (saveBothVtx_) {
645 towerTree->Branch("etVtx",myTowers.etVtx,"etvtx[n]/F");
646 towerTree->Branch("etaVtx",myTowers.etaVtx,"etavtx[n]/F");
647 towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
648 towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
649 }
650
651
652 if(doBasicClusters_){
653 bcTree = fs->make<TTree>("bc",versionTag);
654 bcTree->Branch("n",&myBC.n,"n/I");
655 bcTree->Branch("e",myBC.e,"e[n]/F");
656 bcTree->Branch("et",myBC.et,"et[n]/F");
657 bcTree->Branch("eta",myBC.eta,"eta[n]/F");
658 bcTree->Branch("phi",myBC.phi,"phi[n]/F");
659 bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
660 bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
661 bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
662 // bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
663 }
664
665 if(doFastJet_){
666 bkgTree = fs->make<TTree>("bkg",versionTag);
667 bkgTree->Branch("n",&bkg.n,"n/I");
668 bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
669 bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
670 }
671
672 nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
673
674 }
675
676 // ------------ method called once each job just after ending the event loop ------------
677 void
678 RecHitTreeProducer::endJob() {
679 }
680
681 math::XYZPoint RecHitTreeProducer::getPosition(const DetId &id, reco::Vertex::Point vtx){
682 const GlobalPoint& pos=geo->getPosition(id);
683 math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
684 return posV;
685 }
686
687 double RecHitTreeProducer::getEt(math::XYZPoint pos, double energy){
688 double et = energy*sin(pos.theta());
689 return et;
690 }
691
692 double RecHitTreeProducer::getEt(const DetId &id, double energy, reco::Vertex::Point vtx){
693 const GlobalPoint& pos=geo->getPosition(id);
694 double et = energy*sin(pos.theta());
695 return et;
696 }
697
698 double RecHitTreeProducer::getEta(const DetId &id, reco::Vertex::Point vtx){
699 const GlobalPoint& pos=geo->getPosition(id);
700 double et = pos.eta();
701 return et;
702 }
703
704 double RecHitTreeProducer::getPhi(const DetId &id, reco::Vertex::Point vtx){
705 const GlobalPoint& pos=geo->getPosition(id);
706 double et = pos.phi();
707 return et;
708 }
709
710 double RecHitTreeProducer::getPerp(const DetId &id, reco::Vertex::Point vtx){
711 const GlobalPoint& pos=geo->getPosition(id);
712 double et = pos.perp();
713 return et;
714 }
715
716 reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
717 {
718 ev.getByLabel(VtxSrc_,vtxs);
719 int greatestvtx = 0;
720 int nVertex = vtxs->size();
721
722 for (unsigned int i = 0 ; i< vtxs->size(); ++i){
723 unsigned int daughter = (*vtxs)[i].tracksSize();
724 if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
725 //cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
726 }
727
728 if(nVertex<=0){
729 return reco::Vertex::Point(0,0,0);
730 }
731 return (*vtxs)[greatestvtx].position();
732 }
733
734 //define this as a plug-in
735 DEFINE_FWK_MODULE(RecHitTreeProducer);