| 12 |
|
*/ |
| 13 |
|
// |
| 14 |
|
// Original Author: Yetkin Yilmaz |
| 15 |
< |
// Modified: Frank Ma |
| 15 |
> |
// Modified: Frank Ma, Yen-Jie Lee |
| 16 |
|
// Created: Tue Sep 7 11:38:19 EDT 2010 |
| 17 |
|
// $Id$ |
| 18 |
|
// |
| 19 |
|
// |
| 20 |
|
|
| 21 |
< |
|
| 21 |
> |
#define versionTag "v1" |
| 22 |
|
// system include files |
| 23 |
|
#include <memory> |
| 24 |
|
#include <vector> |
| 56 |
|
|
| 57 |
|
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" |
| 58 |
|
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" |
| 59 |
+ |
#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h" |
| 60 |
|
#include "DataFormats/HcalDetId/interface/HcalDetId.h" |
| 61 |
|
#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h" |
| 62 |
|
#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h" |
| 65 |
|
#include "DataFormats/VertexReco/interface/Vertex.h" |
| 66 |
|
#include "DataFormats/VertexReco/interface/VertexFwd.h" |
| 67 |
|
|
| 68 |
+ |
#include "CalibFormats/HcalObjects/interface/HcalCoderDb.h" |
| 69 |
+ |
#include "CalibFormats/HcalObjects/interface/HcalDbService.h" |
| 70 |
+ |
#include "CalibFormats/HcalObjects/interface/HcalDbRecord.h" |
| 71 |
+ |
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h" |
| 72 |
|
|
| 73 |
|
#include "TNtuple.h" |
| 74 |
|
|
| 75 |
|
using namespace std; |
| 76 |
|
|
| 77 |
< |
#define MAXHITS 1000000 |
| 77 |
> |
#define MAXHITS 100000 |
| 78 |
> |
|
| 79 |
|
|
| 80 |
|
struct MyRecHit{ |
| 81 |
|
int depth[MAXHITS]; |
| 82 |
|
int n; |
| 83 |
|
|
| 84 |
+ |
int ieta[MAXHITS]; |
| 85 |
+ |
int iphi[MAXHITS]; |
| 86 |
+ |
|
| 87 |
|
float e[MAXHITS]; |
| 88 |
|
float et[MAXHITS]; |
| 89 |
|
float eta[MAXHITS]; |
| 90 |
|
float phi[MAXHITS]; |
| 91 |
+ |
float perp[MAXHITS]; |
| 92 |
+ |
float emEt[MAXHITS]; |
| 93 |
+ |
float hadEt[MAXHITS]; |
| 94 |
+ |
|
| 95 |
|
bool isjet[MAXHITS]; |
| 96 |
< |
float etvtx[MAXHITS]; |
| 97 |
< |
float etavtx[MAXHITS]; |
| 96 |
> |
float etVtx[MAXHITS]; |
| 97 |
> |
float etaVtx[MAXHITS]; |
| 98 |
> |
float phiVtx[MAXHITS]; |
| 99 |
> |
float perpVtx[MAXHITS]; |
| 100 |
|
float emEtVtx[MAXHITS]; |
| 101 |
|
float hadEtVtx[MAXHITS]; |
| 102 |
|
|
| 106 |
|
|
| 107 |
|
}; |
| 108 |
|
|
| 109 |
+ |
struct MyZDCRecHit{ |
| 110 |
+ |
int n; |
| 111 |
+ |
float e[18]; |
| 112 |
+ |
int zside[18]; |
| 113 |
+ |
int section [18]; |
| 114 |
+ |
int channel[18]; |
| 115 |
+ |
int saturation[18]; |
| 116 |
+ |
}; |
| 117 |
+ |
|
| 118 |
+ |
struct MyZDCDigi{ |
| 119 |
+ |
int n; |
| 120 |
+ |
float chargefC[10][18]; |
| 121 |
+ |
int adc[10][18]; |
| 122 |
+ |
int zside[18]; |
| 123 |
+ |
int section [18]; |
| 124 |
+ |
int channel[18]; |
| 125 |
+ |
}; |
| 126 |
|
|
| 127 |
|
struct MyBkg{ |
| 128 |
|
int n; |
| 139 |
|
public: |
| 140 |
|
explicit RecHitTreeProducer(const edm::ParameterSet&); |
| 141 |
|
~RecHitTreeProducer(); |
| 142 |
< |
double getEt(const DetId &id, double energy); |
| 143 |
< |
double getEta(const DetId &id); |
| 144 |
< |
double getPhi(const DetId &id); |
| 142 |
> |
double getEt(const DetId &id, double energy, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 143 |
> |
double getEta(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 144 |
> |
double getPhi(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 145 |
> |
double getPerp(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 146 |
> |
|
| 147 |
> |
math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 148 |
> |
double getEt(math::XYZPoint pos, double energy); |
| 149 |
> |
|
| 150 |
|
reco::Vertex::Point getVtx(const edm::Event& ev); |
| 151 |
|
|
| 152 |
|
|
| 183 |
|
MyRecHit hfRecHit; |
| 184 |
|
MyRecHit ebRecHit; |
| 185 |
|
MyRecHit eeRecHit; |
| 186 |
< |
MyRecHit myBC; |
| 187 |
< |
MyRecHit myTowers; |
| 188 |
< |
MyBkg bkg; |
| 186 |
> |
MyRecHit myBC; |
| 187 |
> |
MyRecHit myTowers; |
| 188 |
> |
MyRecHit castorRecHit; |
| 189 |
> |
|
| 190 |
> |
MyZDCRecHit zdcRecHit; |
| 191 |
> |
MyZDCDigi zdcDigi; |
| 192 |
> |
|
| 193 |
> |
MyBkg bkg; |
| 194 |
|
|
| 195 |
|
TNtuple* nt; |
| 196 |
|
TTree* hbheTree; |
| 200 |
|
TTree* bcTree; |
| 201 |
|
TTree* towerTree; |
| 202 |
|
TTree* bkgTree; |
| 203 |
+ |
TTree* castorTree; |
| 204 |
+ |
TTree* zdcRecHitTree; |
| 205 |
+ |
TTree* zdcDigiTree; |
| 206 |
|
|
| 207 |
|
double cone; |
| 208 |
|
double hfTowerThreshold_; |
| 217 |
|
double ebPtMin_; |
| 218 |
|
double eePtMin_; |
| 219 |
|
double towerPtMin_; |
| 220 |
+ |
|
| 221 |
+ |
int nZdcTs_; |
| 222 |
+ |
bool calZDCDigi_; |
| 223 |
|
|
| 224 |
|
edm::Service<TFileService> fs; |
| 225 |
|
const CentralityBins * cbins_; |
| 242 |
|
bool doTowers_; |
| 243 |
|
bool doEcal_; |
| 244 |
|
bool doHcal_; |
| 245 |
+ |
bool doHF_; |
| 246 |
+ |
bool doCastor_; |
| 247 |
+ |
bool doZDCRecHit_; |
| 248 |
+ |
bool doZDCDigi_; |
| 249 |
+ |
|
| 250 |
+ |
bool hasVtx_; |
| 251 |
+ |
bool saveBothVtx_; |
| 252 |
|
|
| 253 |
|
bool doFastJet_; |
| 254 |
|
|
| 255 |
|
bool doEbyEonly_; |
| 256 |
+ |
|
| 257 |
|
}; |
| 258 |
|
|
| 259 |
|
// |
| 286 |
|
doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true); |
| 287 |
|
doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true); |
| 288 |
|
doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true); |
| 289 |
+ |
doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true); |
| 290 |
+ |
doCastor_ = iConfig.getUntrackedParameter<bool>("doCASTOR",true); |
| 291 |
+ |
doZDCRecHit_ = iConfig.getUntrackedParameter<bool>("doZDCRecHit",true); |
| 292 |
+ |
doZDCDigi_ = iConfig.getUntrackedParameter<bool>("doZDCDigi",true); |
| 293 |
+ |
|
| 294 |
+ |
hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",true); |
| 295 |
+ |
saveBothVtx_ = iConfig.getUntrackedParameter<bool>("saveBothVtx",false); |
| 296 |
+ |
|
| 297 |
|
doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true); |
| 298 |
|
FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets")); |
| 299 |
|
doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false); |
| 305 |
|
ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0); |
| 306 |
|
eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.); |
| 307 |
|
towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.); |
| 308 |
+ |
nZdcTs_=iConfig.getUntrackedParameter<int>("nZdcTs",10); |
| 309 |
+ |
calZDCDigi_=iConfig.getUntrackedParameter<bool>("calZDCDigi",true); |
| 310 |
|
} |
| 311 |
|
|
| 312 |
|
|
| 337 |
|
bkg.n = 0; |
| 338 |
|
|
| 339 |
|
// get vertex |
| 340 |
< |
reco::Vertex::Point vtx = getVtx(ev); |
| 340 |
> |
reco::Vertex::Point vtx(0,0,0); |
| 341 |
> |
if (hasVtx_) vtx = getVtx(ev); |
| 342 |
|
|
| 343 |
< |
if(doEcal_){ |
| 344 |
< |
ev.getByLabel(EBSrc_,ebHits); |
| 345 |
< |
ev.getByLabel(EESrc_,eeHits); |
| 343 |
> |
if(doEcal_){ |
| 344 |
> |
ev.getByLabel(EBSrc_,ebHits); |
| 345 |
> |
ev.getByLabel(EESrc_,eeHits); |
| 346 |
|
} |
| 347 |
+ |
|
| 348 |
|
if(doHcal_){ |
| 349 |
< |
ev.getByLabel(HcalRecHitHFSrc_,hfHits); |
| 282 |
< |
ev.getByLabel(HcalRecHitHBHESrc_,hbheHits); |
| 349 |
> |
ev.getByLabel(HcalRecHitHBHESrc_,hbheHits); |
| 350 |
|
} |
| 351 |
+ |
if(doHF_){ |
| 352 |
+ |
ev.getByLabel(HcalRecHitHFSrc_,hfHits); |
| 353 |
+ |
} |
| 354 |
+ |
|
| 355 |
|
if(useJets_) { |
| 356 |
|
ev.getByLabel(JetSrc_,jets); |
| 357 |
|
} |
| 388 |
|
int nHFlongMinus = 0; |
| 389 |
|
int nHFshortMinus = 0; |
| 390 |
|
int nHFtowerMinus = 0; |
| 391 |
< |
|
| 392 |
< |
|
| 322 |
< |
|
| 323 |
< |
if(doHcal_){ |
| 391 |
> |
|
| 392 |
> |
if(doHF_){ |
| 393 |
|
for(unsigned int i = 0; i < hfHits->size(); ++i){ |
| 394 |
|
const HFRecHit & hit= (*hfHits)[i]; |
| 395 |
|
hfRecHit.e[hfRecHit.n] = hit.energy(); |
| 396 |
< |
hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy()); |
| 397 |
< |
hfRecHit.eta[hfRecHit.n] = getEta(hit.id()); |
| 398 |
< |
hfRecHit.phi[hfRecHit.n] = getPhi(hit.id()); |
| 396 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 397 |
> |
|
| 398 |
> |
if(!saveBothVtx_){ |
| 399 |
> |
hfRecHit.et[hfRecHit.n] = getEt(pos,hit.energy()); |
| 400 |
> |
hfRecHit.eta[hfRecHit.n] = pos.eta(); |
| 401 |
> |
hfRecHit.phi[hfRecHit.n] = pos.phi(); |
| 402 |
> |
hfRecHit.perp[hfRecHit.n] = pos.rho(); |
| 403 |
> |
}else{ |
| 404 |
> |
hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy()); |
| 405 |
> |
hfRecHit.eta[hfRecHit.n] = getEta(hit.id()); |
| 406 |
> |
hfRecHit.phi[hfRecHit.n] = getPhi(hit.id()); |
| 407 |
> |
hfRecHit.perp[hfRecHit.n] = getPerp(hit.id()); |
| 408 |
> |
|
| 409 |
> |
hfRecHit.etVtx[hfRecHit.n] = getEt(pos,hit.energy()); |
| 410 |
> |
hfRecHit.etaVtx[hfRecHit.n] = pos.eta(); |
| 411 |
> |
hfRecHit.phiVtx[hfRecHit.n] = pos.phi(); |
| 412 |
> |
hfRecHit.perpVtx[hfRecHit.n] = pos.rho(); |
| 413 |
> |
|
| 414 |
> |
} |
| 415 |
> |
|
| 416 |
|
hfRecHit.isjet[hfRecHit.n] = false; |
| 417 |
|
hfRecHit.depth[hfRecHit.n] = hit.id().depth(); |
| 418 |
|
|
| 419 |
|
if(hit.id().ieta() > 0){ |
| 420 |
< |
if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++; |
| 421 |
< |
if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++; |
| 420 |
> |
if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++; |
| 421 |
> |
if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++; |
| 422 |
|
}else{ |
| 423 |
|
if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++; |
| 424 |
|
if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++; |
| 433 |
|
} |
| 434 |
|
if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++; |
| 435 |
|
} |
| 436 |
< |
if(!doEbyEonly_){ |
| 436 |
> |
|
| 437 |
> |
if(doHcal_ && !doEbyEonly_){ |
| 438 |
|
for(unsigned int i = 0; i < hbheHits->size(); ++i){ |
| 439 |
|
const HBHERecHit & hit= (*hbheHits)[i]; |
| 440 |
|
if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue; |
| 441 |
+ |
|
| 442 |
|
hbheRecHit.e[hbheRecHit.n] = hit.energy(); |
| 443 |
< |
hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy()); |
| 444 |
< |
hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id()); |
| 445 |
< |
hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id()); |
| 446 |
< |
hbheRecHit.isjet[hbheRecHit.n] = false; |
| 443 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 444 |
> |
|
| 445 |
> |
if(!saveBothVtx_){ |
| 446 |
> |
hbheRecHit.et[hbheRecHit.n] = getEt(pos,hit.energy()); |
| 447 |
> |
hbheRecHit.eta[hbheRecHit.n] = pos.eta(); |
| 448 |
> |
hbheRecHit.phi[hbheRecHit.n] = pos.phi(); |
| 449 |
> |
hbheRecHit.perp[hbheRecHit.n] = pos.rho(); |
| 450 |
> |
}else{ |
| 451 |
> |
hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy()); |
| 452 |
> |
hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id()); |
| 453 |
> |
hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id()); |
| 454 |
> |
hbheRecHit.perp[hbheRecHit.n] = getPerp(hit.id()); |
| 455 |
> |
|
| 456 |
> |
hbheRecHit.etVtx[hbheRecHit.n] = getEt(pos,hit.energy()); |
| 457 |
> |
hbheRecHit.etaVtx[hbheRecHit.n] = pos.eta(); |
| 458 |
> |
hbheRecHit.phiVtx[hbheRecHit.n] = pos.phi(); |
| 459 |
> |
hbheRecHit.perpVtx[hbheRecHit.n] = pos.rho(); |
| 460 |
> |
|
| 461 |
> |
} |
| 462 |
> |
|
| 463 |
> |
hbheRecHit.isjet[hbheRecHit.n] = false; |
| 464 |
> |
hbheRecHit.depth[hbheRecHit.n] = hit.id().depth(); |
| 465 |
> |
|
| 466 |
|
if(useJets_){ |
| 467 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 468 |
|
const reco::Jet& jet = (*jets)[j]; |
| 478 |
|
for(unsigned int i = 0; i < ebHits->size(); ++i){ |
| 479 |
|
const EcalRecHit & hit= (*ebHits)[i]; |
| 480 |
|
if (getEt(hit.id(),hit.energy())<ebPtMin_) continue; |
| 481 |
+ |
|
| 482 |
|
ebRecHit.e[ebRecHit.n] = hit.energy(); |
| 483 |
< |
ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy()); |
| 484 |
< |
ebRecHit.eta[ebRecHit.n] = getEta(hit.id()); |
| 485 |
< |
ebRecHit.phi[ebRecHit.n] = getPhi(hit.id()); |
| 483 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 484 |
> |
|
| 485 |
> |
if(!saveBothVtx_){ |
| 486 |
> |
ebRecHit.et[ebRecHit.n] = getEt(pos,hit.energy()); |
| 487 |
> |
ebRecHit.eta[ebRecHit.n] = pos.eta(); |
| 488 |
> |
ebRecHit.phi[ebRecHit.n] = pos.phi(); |
| 489 |
> |
ebRecHit.perp[ebRecHit.n] = pos.rho(); |
| 490 |
> |
}else{ |
| 491 |
> |
ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy()); |
| 492 |
> |
ebRecHit.eta[ebRecHit.n] = getEta(hit.id()); |
| 493 |
> |
ebRecHit.phi[ebRecHit.n] = getPhi(hit.id()); |
| 494 |
> |
ebRecHit.perp[ebRecHit.n] = getPerp(hit.id()); |
| 495 |
> |
ebRecHit.etVtx[ebRecHit.n] = getEt(pos,hit.energy()); |
| 496 |
> |
ebRecHit.etaVtx[ebRecHit.n] = pos.eta(); |
| 497 |
> |
ebRecHit.phiVtx[ebRecHit.n] = pos.phi(); |
| 498 |
> |
ebRecHit.perpVtx[ebRecHit.n] = pos.rho(); |
| 499 |
> |
|
| 500 |
> |
} |
| 501 |
> |
|
| 502 |
|
ebRecHit.isjet[ebRecHit.n] = false; |
| 503 |
|
if(useJets_){ |
| 504 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 513 |
|
for(unsigned int i = 0; i < eeHits->size(); ++i){ |
| 514 |
|
const EcalRecHit & hit= (*eeHits)[i]; |
| 515 |
|
if (getEt(hit.id(),hit.energy())<eePtMin_) continue; |
| 516 |
+ |
|
| 517 |
|
eeRecHit.e[eeRecHit.n] = hit.energy(); |
| 518 |
< |
eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy()); |
| 519 |
< |
eeRecHit.eta[eeRecHit.n] = getEta(hit.id()); |
| 520 |
< |
eeRecHit.phi[eeRecHit.n] = getPhi(hit.id()); |
| 518 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 519 |
> |
|
| 520 |
> |
if(!saveBothVtx_){ |
| 521 |
> |
eeRecHit.et[eeRecHit.n] = getEt(pos,hit.energy()); |
| 522 |
> |
eeRecHit.eta[eeRecHit.n] = pos.eta(); |
| 523 |
> |
eeRecHit.phi[eeRecHit.n] = pos.phi(); |
| 524 |
> |
eeRecHit.perp[eeRecHit.n] = pos.rho(); |
| 525 |
> |
}else{ |
| 526 |
> |
eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy()); |
| 527 |
> |
eeRecHit.eta[eeRecHit.n] = getEta(hit.id()); |
| 528 |
> |
eeRecHit.phi[eeRecHit.n] = getPhi(hit.id()); |
| 529 |
> |
eeRecHit.perp[eeRecHit.n] = getPerp(hit.id()); |
| 530 |
> |
eeRecHit.etVtx[eeRecHit.n] = getEt(pos,hit.energy()); |
| 531 |
> |
eeRecHit.etaVtx[eeRecHit.n] = pos.eta(); |
| 532 |
> |
eeRecHit.phiVtx[eeRecHit.n] = pos.phi(); |
| 533 |
> |
eeRecHit.perpVtx[eeRecHit.n] = pos.rho(); |
| 534 |
> |
|
| 535 |
> |
} |
| 536 |
> |
|
| 537 |
|
eeRecHit.isjet[eeRecHit.n] = false; |
| 538 |
+ |
|
| 539 |
|
if(useJets_){ |
| 540 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 541 |
|
const reco::Jet& jet = (*jets)[j]; |
| 552 |
|
for(unsigned int i = 0; i < towers->size(); ++i){ |
| 553 |
|
const CaloTower & hit= (*towers)[i]; |
| 554 |
|
if (getEt(hit.id(),hit.energy())<towerPtMin_) continue; |
| 555 |
< |
myTowers.e[myTowers.n] = hit.energy(); |
| 556 |
< |
myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy()); |
| 557 |
< |
myTowers.eta[myTowers.n] = getEta(hit.id()); |
| 558 |
< |
myTowers.phi[myTowers.n] = getPhi(hit.id()); |
| 555 |
> |
|
| 556 |
> |
myTowers.et[myTowers.n] = hit.p4(vtx).Et(); |
| 557 |
> |
myTowers.eta[myTowers.n] = hit.p4(vtx).Eta(); |
| 558 |
> |
myTowers.phi[myTowers.n] = hit.p4(vtx).Phi(); |
| 559 |
> |
myTowers.emEt[myTowers.n] = hit.emEt(vtx); |
| 560 |
> |
myTowers.hadEt[myTowers.n] = hit.hadEt(vtx); |
| 561 |
> |
|
| 562 |
> |
if (saveBothVtx_) { |
| 563 |
> |
myTowers.e[myTowers.n] = hit.energy(); |
| 564 |
> |
myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy()); |
| 565 |
> |
myTowers.eta[myTowers.n] = getEta(hit.id()); |
| 566 |
> |
myTowers.phi[myTowers.n] = getPhi(hit.id()); |
| 567 |
> |
myTowers.isjet[myTowers.n] = false; |
| 568 |
> |
myTowers.etVtx[myTowers.n] = hit.p4(vtx).Et(); |
| 569 |
> |
myTowers.etaVtx[myTowers.n] = hit.p4(vtx).Eta(); |
| 570 |
> |
myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx); |
| 571 |
> |
myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx); |
| 572 |
> |
} |
| 573 |
> |
|
| 574 |
|
myTowers.isjet[myTowers.n] = false; |
| 418 |
– |
myTowers.etvtx[myTowers.n] = hit.p4(vtx).Et(); |
| 419 |
– |
myTowers.etavtx[myTowers.n] = hit.p4(vtx).Eta(); |
| 420 |
– |
myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx); |
| 421 |
– |
myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx); |
| 575 |
|
|
| 576 |
|
if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++; |
| 577 |
|
if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++; |
| 611 |
|
} |
| 612 |
|
} |
| 613 |
|
|
| 614 |
+ |
if(doCastor_){ |
| 615 |
+ |
|
| 616 |
+ |
edm::Handle<CastorRecHitCollection> casrechits; |
| 617 |
+ |
try{ ev.getByLabel("castorreco",casrechits); } |
| 618 |
+ |
catch(...) { edm::LogWarning(" CASTOR ") << " Cannot get Castor RecHits " << std::endl; } |
| 619 |
+ |
|
| 620 |
+ |
int nhits = 0; |
| 621 |
+ |
double energyCastor = 0; |
| 622 |
+ |
|
| 623 |
+ |
if(casrechits.failedToGet()!=0 || !casrechits.isValid()) { |
| 624 |
+ |
edm::LogWarning(" CASTOR ") << " Cannot read CastorRecHitCollection" << std::endl; |
| 625 |
+ |
} else { |
| 626 |
+ |
for(size_t i1 = 0; i1 < casrechits->size(); ++i1) { |
| 627 |
+ |
const CastorRecHit & rh = (*casrechits)[i1]; |
| 628 |
+ |
HcalCastorDetId castorid = rh.id(); |
| 629 |
+ |
energyCastor += rh.energy(); |
| 630 |
+ |
if (nhits < 224) { |
| 631 |
+ |
castorRecHit.e[nhits] = rh.energy(); |
| 632 |
+ |
castorRecHit.iphi[nhits] = castorid.sector(); |
| 633 |
+ |
castorRecHit.depth[nhits] = castorid.module(); |
| 634 |
+ |
castorRecHit.phi[nhits] = getPhi(castorid); |
| 635 |
+ |
} |
| 636 |
+ |
|
| 637 |
+ |
nhits++; |
| 638 |
+ |
|
| 639 |
+ |
} // end loop castor rechits |
| 640 |
+ |
} |
| 641 |
+ |
|
| 642 |
+ |
castorRecHit.n = nhits; |
| 643 |
+ |
castorTree->Fill(); |
| 644 |
+ |
} |
| 645 |
+ |
|
| 646 |
+ |
if(doZDCRecHit_){ |
| 647 |
+ |
|
| 648 |
+ |
edm::Handle<ZDCRecHitCollection> zdcrechits; |
| 649 |
+ |
|
| 650 |
+ |
try{ ev.getByLabel("zdcreco",zdcrechits); } |
| 651 |
+ |
catch(...) { edm::LogWarning(" ZDC ") << " Cannot get ZDC RecHits " << std::endl; } |
| 652 |
+ |
|
| 653 |
+ |
int nhits = 0; |
| 654 |
+ |
|
| 655 |
+ |
if (zdcrechits.failedToGet()!=0 || !zdcrechits.isValid()) { |
| 656 |
+ |
edm::LogWarning(" ZDC ") << " Cannot read ZDCRecHitCollection" << std::endl; |
| 657 |
+ |
} else { |
| 658 |
+ |
for(size_t i1 = 0; i1 < zdcrechits->size(); ++i1) { |
| 659 |
+ |
const ZDCRecHit & rh = (*zdcrechits)[i1]; |
| 660 |
+ |
HcalZDCDetId zdcid = rh.id(); |
| 661 |
+ |
if (nhits < 18) { |
| 662 |
+ |
zdcRecHit.e[nhits] = rh.energy(); |
| 663 |
+ |
zdcRecHit.zside[nhits] = zdcid.zside(); |
| 664 |
+ |
zdcRecHit.section[nhits] = zdcid.section(); |
| 665 |
+ |
zdcRecHit.channel[nhits] = zdcid.channel(); |
| 666 |
+ |
zdcRecHit.saturation[nhits] = static_cast<int>( rh.flagField(HcalCaloFlagLabels::ADCSaturationBit) ); |
| 667 |
+ |
} |
| 668 |
+ |
|
| 669 |
+ |
nhits++; |
| 670 |
+ |
|
| 671 |
+ |
} // end loop zdc rechits |
| 672 |
+ |
} |
| 673 |
+ |
|
| 674 |
+ |
zdcRecHit.n = nhits; |
| 675 |
+ |
zdcRecHitTree->Fill(); |
| 676 |
+ |
|
| 677 |
+ |
} |
| 678 |
+ |
|
| 679 |
+ |
if(doZDCDigi_){ |
| 680 |
+ |
|
| 681 |
+ |
edm::Handle<ZDCDigiCollection> zdcdigis; |
| 682 |
+ |
|
| 683 |
+ |
try{ ev.getByLabel("hcalDigis",zdcdigis); } |
| 684 |
+ |
catch(...) { edm::LogWarning(" ZDC ") << " Cannot get ZDC Digis " << std::endl; } |
| 685 |
+ |
|
| 686 |
+ |
int nhits = 0; |
| 687 |
+ |
|
| 688 |
+ |
if (zdcdigis.failedToGet()!=0 || !zdcdigis.isValid()) { |
| 689 |
+ |
edm::LogWarning(" ZDC ") << " Cannot read ZDCDigiCollection" << std::endl; |
| 690 |
+ |
} else { |
| 691 |
+ |
|
| 692 |
+ |
edm::ESHandle<HcalDbService> conditions; |
| 693 |
+ |
iSetup.get<HcalDbRecord>().get(conditions); |
| 694 |
+ |
|
| 695 |
+ |
for(size_t i1 = 0; i1 < zdcdigis->size(); ++i1) { |
| 696 |
+ |
CaloSamples caldigi; |
| 697 |
+ |
const ZDCDataFrame & rh = (*zdcdigis)[i1]; |
| 698 |
+ |
HcalZDCDetId zdcid = rh.id(); |
| 699 |
+ |
|
| 700 |
+ |
if(calZDCDigi_){ |
| 701 |
+ |
const HcalQIEShape* qieshape=conditions->getHcalShape(); |
| 702 |
+ |
const HcalQIECoder* qiecoder=conditions->getHcalCoder(zdcid); |
| 703 |
+ |
HcalCoderDb coder(*qiecoder,*qieshape); |
| 704 |
+ |
coder.adc2fC(rh,caldigi); |
| 705 |
+ |
} |
| 706 |
+ |
|
| 707 |
+ |
if (nhits < 18) { |
| 708 |
+ |
int ts = 0; |
| 709 |
+ |
zdcDigi.zside[nhits] = zdcid.zside(); |
| 710 |
+ |
zdcDigi.section[nhits] = zdcid.section(); |
| 711 |
+ |
zdcDigi.channel[nhits] = zdcid.channel(); |
| 712 |
+ |
|
| 713 |
+ |
for(int j1 = 0; j1 < rh.size(); j1++){ |
| 714 |
+ |
zdcDigi.chargefC[ts][nhits]=calZDCDigi_?caldigi[ts]:rh[ts].nominal_fC(); |
| 715 |
+ |
zdcDigi.adc[ts][nhits]= rh[ts].adc(); |
| 716 |
+ |
ts++; |
| 717 |
+ |
} |
| 718 |
+ |
} |
| 719 |
+ |
nhits++; |
| 720 |
+ |
} // end loop zdc rechits |
| 721 |
+ |
} |
| 722 |
+ |
|
| 723 |
+ |
zdcDigi.n = nhits; |
| 724 |
+ |
zdcDigiTree->Fill(); |
| 725 |
+ |
|
| 726 |
+ |
} |
| 727 |
+ |
|
| 728 |
|
if(!doEbyEonly_){ |
| 729 |
|
towerTree->Fill(); |
| 730 |
|
|
| 749 |
|
RecHitTreeProducer::beginJob() |
| 750 |
|
{ |
| 751 |
|
|
| 752 |
< |
hbheTree = fs->make<TTree>("hbhe",""); |
| 752 |
> |
hbheTree = fs->make<TTree>("hbhe",versionTag); |
| 753 |
|
hbheTree->Branch("n",&hbheRecHit.n,"n/I"); |
| 754 |
|
hbheTree->Branch("e",hbheRecHit.e,"e[n]/F"); |
| 755 |
|
hbheTree->Branch("et",hbheRecHit.et,"et[n]/F"); |
| 756 |
|
hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F"); |
| 757 |
|
hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F"); |
| 758 |
+ |
hbheTree->Branch("perp",hbheRecHit.perp,"perp[n]/F"); |
| 759 |
+ |
|
| 760 |
|
hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O"); |
| 761 |
< |
|
| 762 |
< |
hfTree = fs->make<TTree>("hf",""); |
| 761 |
> |
hbheTree->Branch("depth",hfRecHit.depth,"depth[n]/I"); |
| 762 |
> |
|
| 763 |
> |
hfTree = fs->make<TTree>("hf",versionTag); |
| 764 |
|
hfTree->Branch("n",&hfRecHit.n,"n/I"); |
| 765 |
|
hfTree->Branch("e",hfRecHit.e,"e[n]/F"); |
| 766 |
|
hfTree->Branch("et",hfRecHit.et,"et[n]/F"); |
| 767 |
|
hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F"); |
| 768 |
|
hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F"); |
| 769 |
+ |
hfTree->Branch("perp",hfRecHit.perp,"perp[n]/F"); |
| 770 |
|
hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I"); |
| 771 |
|
hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O"); |
| 772 |
|
|
| 773 |
< |
eeTree = fs->make<TTree>("ee",""); |
| 773 |
> |
eeTree = fs->make<TTree>("ee",versionTag); |
| 774 |
|
eeTree->Branch("n",&eeRecHit.n,"n/I"); |
| 775 |
|
eeTree->Branch("e",eeRecHit.e,"e[n]/F"); |
| 776 |
|
eeTree->Branch("et",eeRecHit.et,"et[n]/F"); |
| 777 |
|
eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F"); |
| 778 |
|
eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F"); |
| 779 |
+ |
eeTree->Branch("perp",eeRecHit.perp,"perp[n]/F"); |
| 780 |
+ |
|
| 781 |
|
eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O"); |
| 782 |
|
|
| 783 |
< |
ebTree = fs->make<TTree>("eb",""); |
| 783 |
> |
ebTree = fs->make<TTree>("eb",versionTag); |
| 784 |
|
ebTree->Branch("n",&ebRecHit.n,"n/I"); |
| 785 |
|
ebTree->Branch("e",ebRecHit.e,"e[n]/F"); |
| 786 |
|
ebTree->Branch("et",ebRecHit.et,"et[n]/F"); |
| 787 |
|
ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F"); |
| 788 |
|
ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F"); |
| 789 |
+ |
ebTree->Branch("perp",ebRecHit.perp,"perp[n]/F"); |
| 790 |
+ |
|
| 791 |
|
ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O"); |
| 792 |
|
|
| 793 |
< |
towerTree = fs->make<TTree>("tower",""); |
| 793 |
> |
towerTree = fs->make<TTree>("tower",versionTag); |
| 794 |
|
towerTree->Branch("n",&myTowers.n,"n/I"); |
| 795 |
|
towerTree->Branch("e",myTowers.e,"e[n]/F"); |
| 796 |
|
towerTree->Branch("et",myTowers.et,"et[n]/F"); |
| 797 |
|
towerTree->Branch("eta",myTowers.eta,"eta[n]/F"); |
| 798 |
|
towerTree->Branch("phi",myTowers.phi,"phi[n]/F"); |
| 799 |
|
towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O"); |
| 800 |
< |
towerTree->Branch("etvtx",myTowers.etvtx,"etvtx[n]/F"); |
| 801 |
< |
towerTree->Branch("etavtx",myTowers.etavtx,"etavtx[n]/F"); |
| 802 |
< |
towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F"); |
| 803 |
< |
towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F"); |
| 800 |
> |
towerTree->Branch("emEt",myTowers.emEt,"emEt[n]/F"); |
| 801 |
> |
towerTree->Branch("hadEt",myTowers.hadEt,"hadEt[n]/F"); |
| 802 |
> |
|
| 803 |
> |
|
| 804 |
> |
if(doCastor_){ |
| 805 |
> |
castorTree = fs->make<TTree>("castor",versionTag); |
| 806 |
> |
castorTree->Branch("n",&castorRecHit.n,"n/I"); |
| 807 |
> |
castorTree->Branch("e",castorRecHit.e,"e[n]/F"); |
| 808 |
> |
castorTree->Branch("iphi",castorRecHit.iphi,"iphi[n]/I"); |
| 809 |
> |
castorTree->Branch("phi",castorRecHit.phi,"phi[n]/F"); |
| 810 |
> |
castorTree->Branch("depth",castorRecHit.depth,"depth[n]/I"); |
| 811 |
> |
} |
| 812 |
> |
|
| 813 |
> |
if(doZDCRecHit_){ |
| 814 |
> |
zdcRecHitTree = fs->make<TTree>("zdcrechit",versionTag); |
| 815 |
> |
zdcRecHitTree->Branch("n",&zdcRecHit.n,"n/I"); |
| 816 |
> |
zdcRecHitTree->Branch("e",zdcRecHit.e,"e[n]/F"); |
| 817 |
> |
zdcRecHitTree->Branch("saturation",zdcRecHit.saturation,"saturation[n]/F"); |
| 818 |
> |
zdcRecHitTree->Branch("zside",zdcRecHit.zside,"zside[n]/I"); |
| 819 |
> |
zdcRecHitTree->Branch("section",zdcRecHit.section,"section[n]/I"); |
| 820 |
> |
zdcRecHitTree->Branch("channel",zdcRecHit.channel,"channel[n]/I"); |
| 821 |
> |
} |
| 822 |
> |
|
| 823 |
> |
if(doZDCDigi_){ |
| 824 |
> |
TString nZdcTsSt=""; |
| 825 |
> |
nZdcTsSt+=nZdcTs_; |
| 826 |
> |
|
| 827 |
> |
zdcDigiTree = fs->make<TTree>("zdcdigi",versionTag); |
| 828 |
> |
zdcDigiTree->Branch("n",&zdcDigi.n,"n/I"); |
| 829 |
> |
zdcDigiTree->Branch("zside",zdcDigi.zside,"zside[n]/I"); |
| 830 |
> |
zdcDigiTree->Branch("section",zdcDigi.section,"section[n]/I"); |
| 831 |
> |
zdcDigiTree->Branch("channel",zdcDigi.channel,"channel[n]/I"); |
| 832 |
> |
|
| 833 |
> |
for( int i=0; i<nZdcTs_;i++){ |
| 834 |
> |
TString adcTsSt("adcTs"), chargefCTsSt("chargefCTs"); |
| 835 |
> |
adcTsSt+=i; chargefCTsSt+=i; |
| 836 |
> |
|
| 837 |
> |
zdcDigiTree->Branch(adcTsSt,zdcDigi.adc[i],adcTsSt+"[n]/I"); |
| 838 |
> |
zdcDigiTree->Branch(chargefCTsSt,zdcDigi.chargefC[i],chargefCTsSt+"[n]/F"); |
| 839 |
> |
} |
| 840 |
> |
} |
| 841 |
|
|
| 842 |
+ |
if (saveBothVtx_) { |
| 843 |
+ |
towerTree->Branch("etVtx",myTowers.etVtx,"etvtx[n]/F"); |
| 844 |
+ |
towerTree->Branch("etaVtx",myTowers.etaVtx,"etavtx[n]/F"); |
| 845 |
+ |
towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F"); |
| 846 |
+ |
towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F"); |
| 847 |
+ |
} |
| 848 |
|
|
| 849 |
|
if(doBasicClusters_){ |
| 850 |
< |
bcTree = fs->make<TTree>("bc",""); |
| 850 |
> |
bcTree = fs->make<TTree>("bc",versionTag); |
| 851 |
|
bcTree->Branch("n",&myBC.n,"n/I"); |
| 852 |
|
bcTree->Branch("e",myBC.e,"e[n]/F"); |
| 853 |
|
bcTree->Branch("et",myBC.et,"et[n]/F"); |
| 860 |
|
} |
| 861 |
|
|
| 862 |
|
if(doFastJet_){ |
| 863 |
< |
bkgTree = fs->make<TTree>("bkg",""); |
| 863 |
> |
bkgTree = fs->make<TTree>("bkg",versionTag); |
| 864 |
|
bkgTree->Branch("n",&bkg.n,"n/I"); |
| 865 |
|
bkgTree->Branch("rho",bkg.rho,"rho[n]/F"); |
| 866 |
|
bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F"); |
| 875 |
|
RecHitTreeProducer::endJob() { |
| 876 |
|
} |
| 877 |
|
|
| 878 |
< |
double RecHitTreeProducer::getEt(const DetId &id, double energy){ |
| 878 |
> |
math::XYZPoint RecHitTreeProducer::getPosition(const DetId &id, reco::Vertex::Point vtx){ |
| 879 |
|
const GlobalPoint& pos=geo->getPosition(id); |
| 880 |
+ |
math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z()); |
| 881 |
+ |
return posV; |
| 882 |
+ |
} |
| 883 |
+ |
|
| 884 |
+ |
double RecHitTreeProducer::getEt(math::XYZPoint pos, double energy){ |
| 885 |
|
double et = energy*sin(pos.theta()); |
| 886 |
|
return et; |
| 887 |
|
} |
| 888 |
|
|
| 889 |
< |
double RecHitTreeProducer::getEta(const DetId &id){ |
| 889 |
> |
double RecHitTreeProducer::getEt(const DetId &id, double energy, reco::Vertex::Point vtx){ |
| 890 |
> |
const GlobalPoint& pos=geo->getPosition(id); |
| 891 |
> |
double et = energy*sin(pos.theta()); |
| 892 |
> |
return et; |
| 893 |
> |
} |
| 894 |
> |
|
| 895 |
> |
double RecHitTreeProducer::getEta(const DetId &id, reco::Vertex::Point vtx){ |
| 896 |
|
const GlobalPoint& pos=geo->getPosition(id); |
| 897 |
|
double et = pos.eta(); |
| 898 |
|
return et; |
| 899 |
|
} |
| 900 |
|
|
| 901 |
< |
double RecHitTreeProducer::getPhi(const DetId &id){ |
| 901 |
> |
double RecHitTreeProducer::getPhi(const DetId &id, reco::Vertex::Point vtx){ |
| 902 |
|
const GlobalPoint& pos=geo->getPosition(id); |
| 903 |
|
double et = pos.phi(); |
| 904 |
|
return et; |
| 905 |
|
} |
| 906 |
|
|
| 907 |
+ |
double RecHitTreeProducer::getPerp(const DetId &id, reco::Vertex::Point vtx){ |
| 908 |
+ |
const GlobalPoint& pos=geo->getPosition(id); |
| 909 |
+ |
double et = pos.perp(); |
| 910 |
+ |
return et; |
| 911 |
+ |
} |
| 912 |
+ |
|
| 913 |
|
reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev) |
| 914 |
|
{ |
| 915 |
|
ev.getByLabel(VtxSrc_,vtxs); |
| 923 |
|
} |
| 924 |
|
|
| 925 |
|
if(nVertex<=0){ |
| 926 |
< |
return reco::Vertex::Point(-999,-999,-999); |
| 926 |
> |
return reco::Vertex::Point(0,0,0); |
| 927 |
|
} |
| 928 |
|
return (*vtxs)[greatestvtx].position(); |
| 929 |
|
} |
