| 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$ |
| 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" |
| 63 |
|
#include "DataFormats/DetId/interface/DetId.h" |
| 64 |
|
|
| 65 |
+ |
#include "DataFormats/VertexReco/interface/Vertex.h" |
| 66 |
+ |
#include "DataFormats/VertexReco/interface/VertexFwd.h" |
| 67 |
+ |
|
| 68 |
|
|
| 69 |
|
#include "TNtuple.h" |
| 70 |
|
|
| 71 |
|
using namespace std; |
| 72 |
|
|
| 73 |
< |
#define MAXHITS 1000000 |
| 73 |
> |
#define MAXHITS 100000 |
| 74 |
> |
|
| 75 |
|
|
| 76 |
|
struct MyRecHit{ |
| 77 |
|
int depth[MAXHITS]; |
| 78 |
|
int n; |
| 79 |
|
|
| 80 |
+ |
int ieta[MAXHITS]; |
| 81 |
+ |
int iphi[MAXHITS]; |
| 82 |
+ |
|
| 83 |
|
float e[MAXHITS]; |
| 84 |
|
float et[MAXHITS]; |
| 85 |
|
float eta[MAXHITS]; |
| 86 |
|
float phi[MAXHITS]; |
| 87 |
+ |
float perp[MAXHITS]; |
| 88 |
+ |
float emEt[MAXHITS]; |
| 89 |
+ |
float hadEt[MAXHITS]; |
| 90 |
+ |
|
| 91 |
|
bool isjet[MAXHITS]; |
| 92 |
+ |
float etVtx[MAXHITS]; |
| 93 |
+ |
float etaVtx[MAXHITS]; |
| 94 |
+ |
float phiVtx[MAXHITS]; |
| 95 |
+ |
float perpVtx[MAXHITS]; |
| 96 |
+ |
float emEtVtx[MAXHITS]; |
| 97 |
+ |
float hadEtVtx[MAXHITS]; |
| 98 |
|
|
| 99 |
|
float jtpt; |
| 100 |
|
float jteta; |
| 102 |
|
|
| 103 |
|
}; |
| 104 |
|
|
| 105 |
+ |
struct MyZDCRecHit{ |
| 106 |
+ |
int n; |
| 107 |
+ |
float e[18]; |
| 108 |
+ |
int zside[18]; |
| 109 |
+ |
int section [18]; |
| 110 |
+ |
int channel[18]; |
| 111 |
+ |
int saturation[18]; |
| 112 |
+ |
}; |
| 113 |
+ |
|
| 114 |
+ |
struct MyZDCDigi{ |
| 115 |
+ |
int n; |
| 116 |
+ |
int nts; |
| 117 |
+ |
float chargefC[18][10]; |
| 118 |
+ |
int adc[18][10]; |
| 119 |
+ |
int zside[18][10]; |
| 120 |
+ |
int section [18][10]; |
| 121 |
+ |
int channel[18][10]; |
| 122 |
+ |
int ts[18][10]; |
| 123 |
+ |
}; |
| 124 |
|
|
| 125 |
|
struct MyBkg{ |
| 126 |
|
int n; |
| 137 |
|
public: |
| 138 |
|
explicit RecHitTreeProducer(const edm::ParameterSet&); |
| 139 |
|
~RecHitTreeProducer(); |
| 140 |
< |
double getEt(const DetId &id, double energy); |
| 141 |
< |
double getEta(const DetId &id); |
| 142 |
< |
double getPhi(const DetId &id); |
| 140 |
> |
double getEt(const DetId &id, double energy, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 141 |
> |
double getEta(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 142 |
> |
double getPhi(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 143 |
> |
double getPerp(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 144 |
> |
|
| 145 |
> |
math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); |
| 146 |
> |
double getEt(math::XYZPoint pos, double energy); |
| 147 |
> |
|
| 148 |
> |
reco::Vertex::Point getVtx(const edm::Event& ev); |
| 149 |
> |
|
| 150 |
|
|
| 151 |
|
|
| 152 |
|
private: |
| 168 |
|
|
| 169 |
|
edm::Handle<reco::BasicClusterCollection> bClusters; |
| 170 |
|
edm::Handle<CaloTowerCollection> towers; |
| 171 |
< |
|
| 171 |
> |
edm::Handle<reco::VertexCollection> vtxs; |
| 172 |
|
|
| 173 |
|
typedef vector<EcalRecHit>::const_iterator EcalIterator; |
| 174 |
|
|
| 181 |
|
MyRecHit hfRecHit; |
| 182 |
|
MyRecHit ebRecHit; |
| 183 |
|
MyRecHit eeRecHit; |
| 184 |
< |
MyRecHit myBC; |
| 185 |
< |
MyRecHit myTowers; |
| 186 |
< |
MyBkg bkg; |
| 184 |
> |
MyRecHit myBC; |
| 185 |
> |
MyRecHit myTowers; |
| 186 |
> |
MyRecHit castorRecHit; |
| 187 |
> |
|
| 188 |
> |
MyZDCRecHit zdcRecHit; |
| 189 |
> |
MyZDCDigi zdcDigi; |
| 190 |
> |
|
| 191 |
> |
MyBkg bkg; |
| 192 |
|
|
| 193 |
|
TNtuple* nt; |
| 194 |
|
TTree* hbheTree; |
| 198 |
|
TTree* bcTree; |
| 199 |
|
TTree* towerTree; |
| 200 |
|
TTree* bkgTree; |
| 201 |
+ |
TTree* castorTree; |
| 202 |
+ |
TTree* zdcRecHitTree; |
| 203 |
+ |
TTree* zdcDigiTree; |
| 204 |
|
|
| 205 |
|
double cone; |
| 206 |
|
double hfTowerThreshold_; |
| 209 |
|
double hbheThreshold_; |
| 210 |
|
double ebThreshold_; |
| 211 |
|
double eeThreshold_; |
| 212 |
< |
|
| 212 |
> |
|
| 213 |
> |
double hbhePtMin_; |
| 214 |
> |
double hfPtMin_; |
| 215 |
> |
double ebPtMin_; |
| 216 |
> |
double eePtMin_; |
| 217 |
> |
double towerPtMin_; |
| 218 |
> |
|
| 219 |
|
edm::Service<TFileService> fs; |
| 220 |
|
const CentralityBins * cbins_; |
| 221 |
|
const CaloGeometry *geo; |
| 226 |
|
edm::InputTag EESrc_; |
| 227 |
|
edm::InputTag BCSrc_; |
| 228 |
|
edm::InputTag TowerSrc_; |
| 229 |
+ |
edm::InputTag VtxSrc_; |
| 230 |
|
|
| 231 |
|
edm::InputTag JetSrc_; |
| 232 |
|
|
| 237 |
|
bool doTowers_; |
| 238 |
|
bool doEcal_; |
| 239 |
|
bool doHcal_; |
| 240 |
+ |
bool doHF_; |
| 241 |
+ |
bool doCastor_; |
| 242 |
+ |
bool doZDCRecHit_; |
| 243 |
+ |
bool doZDCDigi_; |
| 244 |
+ |
|
| 245 |
+ |
bool hasVtx_; |
| 246 |
+ |
bool saveBothVtx_; |
| 247 |
|
|
| 248 |
|
bool doFastJet_; |
| 249 |
|
|
| 250 |
|
bool doEbyEonly_; |
| 251 |
+ |
|
| 252 |
|
}; |
| 253 |
|
|
| 254 |
|
// |
| 274 |
|
HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco")); |
| 275 |
|
BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO")); |
| 276 |
|
TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker")); |
| 277 |
< |
JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets")); |
| 277 |
> |
VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex")); |
| 278 |
> |
JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets")); |
| 279 |
|
useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true); |
| 280 |
|
doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false); |
| 281 |
|
doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true); |
| 282 |
|
doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true); |
| 283 |
|
doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true); |
| 284 |
+ |
doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true); |
| 285 |
+ |
doCastor_ = iConfig.getUntrackedParameter<bool>("doCASTOR",true); |
| 286 |
+ |
doZDCRecHit_ = iConfig.getUntrackedParameter<bool>("doZDCRecHit",true); |
| 287 |
+ |
doZDCDigi_ = iConfig.getUntrackedParameter<bool>("doZDCDigi",true); |
| 288 |
+ |
|
| 289 |
+ |
hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",true); |
| 290 |
+ |
saveBothVtx_ = iConfig.getUntrackedParameter<bool>("saveBothVtx",false); |
| 291 |
+ |
|
| 292 |
|
doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true); |
| 293 |
|
FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets")); |
| 294 |
|
doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false); |
| 295 |
|
hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.); |
| 296 |
< |
hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",3.); |
| 297 |
< |
hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",3.); |
| 296 |
> |
hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5); |
| 297 |
> |
hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85); |
| 298 |
> |
hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0); |
| 299 |
> |
hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0); |
| 300 |
> |
ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0); |
| 301 |
> |
eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.); |
| 302 |
> |
towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.); |
| 303 |
> |
|
| 304 |
> |
|
| 305 |
|
} |
| 306 |
|
|
| 307 |
|
|
| 331 |
|
myTowers.n = 0; |
| 332 |
|
bkg.n = 0; |
| 333 |
|
|
| 334 |
< |
if(doEcal_){ |
| 335 |
< |
ev.getByLabel(EBSrc_,ebHits); |
| 336 |
< |
ev.getByLabel(EESrc_,eeHits); |
| 334 |
> |
// get vertex |
| 335 |
> |
reco::Vertex::Point vtx(0,0,0); |
| 336 |
> |
if (hasVtx_) vtx = getVtx(ev); |
| 337 |
> |
|
| 338 |
> |
if(doEcal_){ |
| 339 |
> |
ev.getByLabel(EBSrc_,ebHits); |
| 340 |
> |
ev.getByLabel(EESrc_,eeHits); |
| 341 |
|
} |
| 342 |
+ |
|
| 343 |
|
if(doHcal_){ |
| 344 |
< |
ev.getByLabel(HcalRecHitHFSrc_,hfHits); |
| 345 |
< |
ev.getByLabel(HcalRecHitHBHESrc_,hbheHits); |
| 344 |
> |
ev.getByLabel(HcalRecHitHBHESrc_,hbheHits); |
| 345 |
> |
} |
| 346 |
> |
if(doHF_){ |
| 347 |
> |
ev.getByLabel(HcalRecHitHFSrc_,hfHits); |
| 348 |
|
} |
| 349 |
+ |
|
| 350 |
|
if(useJets_) { |
| 351 |
|
ev.getByLabel(JetSrc_,jets); |
| 352 |
|
} |
| 377 |
|
geo = pGeo.product(); |
| 378 |
|
} |
| 379 |
|
|
| 380 |
< |
int nHFlong = 0; |
| 381 |
< |
int nHFshort = 0; |
| 382 |
< |
int nHFtower = 0; |
| 383 |
< |
|
| 384 |
< |
if(doHcal_){ |
| 380 |
> |
int nHFlongPlus = 0; |
| 381 |
> |
int nHFshortPlus = 0; |
| 382 |
> |
int nHFtowerPlus = 0; |
| 383 |
> |
int nHFlongMinus = 0; |
| 384 |
> |
int nHFshortMinus = 0; |
| 385 |
> |
int nHFtowerMinus = 0; |
| 386 |
> |
|
| 387 |
> |
if(doHF_){ |
| 388 |
|
for(unsigned int i = 0; i < hfHits->size(); ++i){ |
| 389 |
|
const HFRecHit & hit= (*hfHits)[i]; |
| 390 |
|
hfRecHit.e[hfRecHit.n] = hit.energy(); |
| 391 |
< |
hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy()); |
| 392 |
< |
hfRecHit.eta[hfRecHit.n] = getEta(hit.id()); |
| 393 |
< |
hfRecHit.phi[hfRecHit.n] = getPhi(hit.id()); |
| 391 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 392 |
> |
|
| 393 |
> |
if(!saveBothVtx_){ |
| 394 |
> |
hfRecHit.et[hfRecHit.n] = getEt(pos,hit.energy()); |
| 395 |
> |
hfRecHit.eta[hfRecHit.n] = pos.eta(); |
| 396 |
> |
hfRecHit.phi[hfRecHit.n] = pos.phi(); |
| 397 |
> |
hfRecHit.perp[hfRecHit.n] = pos.rho(); |
| 398 |
> |
}else{ |
| 399 |
> |
hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy()); |
| 400 |
> |
hfRecHit.eta[hfRecHit.n] = getEta(hit.id()); |
| 401 |
> |
hfRecHit.phi[hfRecHit.n] = getPhi(hit.id()); |
| 402 |
> |
hfRecHit.perp[hfRecHit.n] = getPerp(hit.id()); |
| 403 |
> |
|
| 404 |
> |
hfRecHit.etVtx[hfRecHit.n] = getEt(pos,hit.energy()); |
| 405 |
> |
hfRecHit.etaVtx[hfRecHit.n] = pos.eta(); |
| 406 |
> |
hfRecHit.phiVtx[hfRecHit.n] = pos.phi(); |
| 407 |
> |
hfRecHit.perpVtx[hfRecHit.n] = pos.rho(); |
| 408 |
> |
|
| 409 |
> |
} |
| 410 |
> |
|
| 411 |
|
hfRecHit.isjet[hfRecHit.n] = false; |
| 412 |
|
hfRecHit.depth[hfRecHit.n] = hit.id().depth(); |
| 413 |
< |
if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshort++; |
| 414 |
< |
if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlong++; |
| 413 |
> |
|
| 414 |
> |
if(hit.id().ieta() > 0){ |
| 415 |
> |
if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++; |
| 416 |
> |
if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++; |
| 417 |
> |
}else{ |
| 418 |
> |
if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++; |
| 419 |
> |
if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++; |
| 420 |
> |
} |
| 421 |
|
|
| 422 |
|
if(useJets_){ |
| 423 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 426 |
|
if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; } |
| 427 |
|
} |
| 428 |
|
} |
| 429 |
< |
hfRecHit.n++; |
| 429 |
> |
if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++; |
| 430 |
|
} |
| 431 |
< |
if(!doEbyEonly_){ |
| 431 |
> |
|
| 432 |
> |
if(doHcal_ && !doEbyEonly_){ |
| 433 |
|
for(unsigned int i = 0; i < hbheHits->size(); ++i){ |
| 434 |
|
const HBHERecHit & hit= (*hbheHits)[i]; |
| 435 |
+ |
if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue; |
| 436 |
+ |
|
| 437 |
|
hbheRecHit.e[hbheRecHit.n] = hit.energy(); |
| 438 |
< |
hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy()); |
| 439 |
< |
hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id()); |
| 440 |
< |
hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id()); |
| 441 |
< |
hbheRecHit.isjet[hbheRecHit.n] = false; |
| 438 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 439 |
> |
|
| 440 |
> |
if(!saveBothVtx_){ |
| 441 |
> |
hbheRecHit.et[hbheRecHit.n] = getEt(pos,hit.energy()); |
| 442 |
> |
hbheRecHit.eta[hbheRecHit.n] = pos.eta(); |
| 443 |
> |
hbheRecHit.phi[hbheRecHit.n] = pos.phi(); |
| 444 |
> |
hbheRecHit.perp[hbheRecHit.n] = pos.rho(); |
| 445 |
> |
}else{ |
| 446 |
> |
hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy()); |
| 447 |
> |
hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id()); |
| 448 |
> |
hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id()); |
| 449 |
> |
hbheRecHit.perp[hbheRecHit.n] = getPerp(hit.id()); |
| 450 |
> |
|
| 451 |
> |
hbheRecHit.etVtx[hbheRecHit.n] = getEt(pos,hit.energy()); |
| 452 |
> |
hbheRecHit.etaVtx[hbheRecHit.n] = pos.eta(); |
| 453 |
> |
hbheRecHit.phiVtx[hbheRecHit.n] = pos.phi(); |
| 454 |
> |
hbheRecHit.perpVtx[hbheRecHit.n] = pos.rho(); |
| 455 |
> |
|
| 456 |
> |
} |
| 457 |
> |
|
| 458 |
> |
hbheRecHit.isjet[hbheRecHit.n] = false; |
| 459 |
> |
hbheRecHit.depth[hbheRecHit.n] = hit.id().depth(); |
| 460 |
> |
|
| 461 |
|
if(useJets_){ |
| 462 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 463 |
|
const reco::Jet& jet = (*jets)[j]; |
| 472 |
|
if(doEcal_ && !doEbyEonly_){ |
| 473 |
|
for(unsigned int i = 0; i < ebHits->size(); ++i){ |
| 474 |
|
const EcalRecHit & hit= (*ebHits)[i]; |
| 475 |
+ |
if (getEt(hit.id(),hit.energy())<ebPtMin_) continue; |
| 476 |
+ |
|
| 477 |
|
ebRecHit.e[ebRecHit.n] = hit.energy(); |
| 478 |
< |
ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy()); |
| 479 |
< |
ebRecHit.eta[ebRecHit.n] = getEta(hit.id()); |
| 480 |
< |
ebRecHit.phi[ebRecHit.n] = getPhi(hit.id()); |
| 478 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 479 |
> |
|
| 480 |
> |
if(!saveBothVtx_){ |
| 481 |
> |
ebRecHit.et[ebRecHit.n] = getEt(pos,hit.energy()); |
| 482 |
> |
ebRecHit.eta[ebRecHit.n] = pos.eta(); |
| 483 |
> |
ebRecHit.phi[ebRecHit.n] = pos.phi(); |
| 484 |
> |
ebRecHit.perp[ebRecHit.n] = pos.rho(); |
| 485 |
> |
}else{ |
| 486 |
> |
ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy()); |
| 487 |
> |
ebRecHit.eta[ebRecHit.n] = getEta(hit.id()); |
| 488 |
> |
ebRecHit.phi[ebRecHit.n] = getPhi(hit.id()); |
| 489 |
> |
ebRecHit.perp[ebRecHit.n] = getPerp(hit.id()); |
| 490 |
> |
ebRecHit.etVtx[ebRecHit.n] = getEt(pos,hit.energy()); |
| 491 |
> |
ebRecHit.etaVtx[ebRecHit.n] = pos.eta(); |
| 492 |
> |
ebRecHit.phiVtx[ebRecHit.n] = pos.phi(); |
| 493 |
> |
ebRecHit.perpVtx[ebRecHit.n] = pos.rho(); |
| 494 |
> |
|
| 495 |
> |
} |
| 496 |
> |
|
| 497 |
|
ebRecHit.isjet[ebRecHit.n] = false; |
| 498 |
|
if(useJets_){ |
| 499 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 507 |
|
|
| 508 |
|
for(unsigned int i = 0; i < eeHits->size(); ++i){ |
| 509 |
|
const EcalRecHit & hit= (*eeHits)[i]; |
| 510 |
+ |
if (getEt(hit.id(),hit.energy())<eePtMin_) continue; |
| 511 |
+ |
|
| 512 |
|
eeRecHit.e[eeRecHit.n] = hit.energy(); |
| 513 |
< |
eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy()); |
| 514 |
< |
eeRecHit.eta[eeRecHit.n] = getEta(hit.id()); |
| 515 |
< |
eeRecHit.phi[eeRecHit.n] = getPhi(hit.id()); |
| 513 |
> |
math::XYZPoint pos = getPosition(hit.id(),vtx); |
| 514 |
> |
|
| 515 |
> |
if(!saveBothVtx_){ |
| 516 |
> |
eeRecHit.et[eeRecHit.n] = getEt(pos,hit.energy()); |
| 517 |
> |
eeRecHit.eta[eeRecHit.n] = pos.eta(); |
| 518 |
> |
eeRecHit.phi[eeRecHit.n] = pos.phi(); |
| 519 |
> |
eeRecHit.perp[eeRecHit.n] = pos.rho(); |
| 520 |
> |
}else{ |
| 521 |
> |
eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy()); |
| 522 |
> |
eeRecHit.eta[eeRecHit.n] = getEta(hit.id()); |
| 523 |
> |
eeRecHit.phi[eeRecHit.n] = getPhi(hit.id()); |
| 524 |
> |
eeRecHit.perp[eeRecHit.n] = getPerp(hit.id()); |
| 525 |
> |
eeRecHit.etVtx[eeRecHit.n] = getEt(pos,hit.energy()); |
| 526 |
> |
eeRecHit.etaVtx[eeRecHit.n] = pos.eta(); |
| 527 |
> |
eeRecHit.phiVtx[eeRecHit.n] = pos.phi(); |
| 528 |
> |
eeRecHit.perpVtx[eeRecHit.n] = pos.rho(); |
| 529 |
> |
|
| 530 |
> |
} |
| 531 |
> |
|
| 532 |
|
eeRecHit.isjet[eeRecHit.n] = false; |
| 533 |
+ |
|
| 534 |
|
if(useJets_){ |
| 535 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 536 |
|
const reco::Jet& jet = (*jets)[j]; |
| 546 |
|
|
| 547 |
|
for(unsigned int i = 0; i < towers->size(); ++i){ |
| 548 |
|
const CaloTower & hit= (*towers)[i]; |
| 549 |
< |
myTowers.e[myTowers.n] = hit.energy(); |
| 550 |
< |
myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy()); |
| 551 |
< |
myTowers.eta[myTowers.n] = getEta(hit.id()); |
| 552 |
< |
myTowers.phi[myTowers.n] = getPhi(hit.id()); |
| 549 |
> |
if (getEt(hit.id(),hit.energy())<towerPtMin_) continue; |
| 550 |
> |
|
| 551 |
> |
myTowers.et[myTowers.n] = hit.p4(vtx).Et(); |
| 552 |
> |
myTowers.eta[myTowers.n] = hit.p4(vtx).Eta(); |
| 553 |
> |
myTowers.phi[myTowers.n] = hit.p4(vtx).Phi(); |
| 554 |
> |
myTowers.emEt[myTowers.n] = hit.emEt(vtx); |
| 555 |
> |
myTowers.hadEt[myTowers.n] = hit.hadEt(vtx); |
| 556 |
> |
|
| 557 |
> |
if (saveBothVtx_) { |
| 558 |
> |
myTowers.e[myTowers.n] = hit.energy(); |
| 559 |
> |
myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy()); |
| 560 |
> |
myTowers.eta[myTowers.n] = getEta(hit.id()); |
| 561 |
> |
myTowers.phi[myTowers.n] = getPhi(hit.id()); |
| 562 |
> |
myTowers.isjet[myTowers.n] = false; |
| 563 |
> |
myTowers.etVtx[myTowers.n] = hit.p4(vtx).Et(); |
| 564 |
> |
myTowers.etaVtx[myTowers.n] = hit.p4(vtx).Eta(); |
| 565 |
> |
myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx); |
| 566 |
> |
myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx); |
| 567 |
> |
} |
| 568 |
> |
|
| 569 |
|
myTowers.isjet[myTowers.n] = false; |
| 570 |
|
|
| 571 |
< |
if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtower++; |
| 571 |
> |
if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++; |
| 572 |
> |
if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++; |
| 573 |
|
|
| 574 |
|
if(useJets_){ |
| 575 |
|
for(unsigned int j = 0 ; j < jets->size(); ++j){ |
| 606 |
|
} |
| 607 |
|
} |
| 608 |
|
|
| 609 |
+ |
if(doCastor_){ |
| 610 |
+ |
|
| 611 |
+ |
edm::Handle<CastorRecHitCollection> casrechits; |
| 612 |
+ |
try{ ev.getByLabel("castorreco",casrechits); } |
| 613 |
+ |
catch(...) { edm::LogWarning(" CASTOR ") << " Cannot get Castor RecHits " << std::endl; } |
| 614 |
+ |
|
| 615 |
+ |
int nhits = 0; |
| 616 |
+ |
double energyCastor = 0; |
| 617 |
+ |
|
| 618 |
+ |
if(casrechits.failedToGet()!=0 || !casrechits.isValid()) { |
| 619 |
+ |
edm::LogWarning(" CASTOR ") << " Cannot read CastorRecHitCollection" << std::endl; |
| 620 |
+ |
} else { |
| 621 |
+ |
for(size_t i1 = 0; i1 < casrechits->size(); ++i1) { |
| 622 |
+ |
const CastorRecHit & rh = (*casrechits)[i1]; |
| 623 |
+ |
HcalCastorDetId castorid = rh.id(); |
| 624 |
+ |
energyCastor += rh.energy(); |
| 625 |
+ |
if (nhits < 224) { |
| 626 |
+ |
castorRecHit.e[nhits] = rh.energy(); |
| 627 |
+ |
castorRecHit.iphi[nhits] = castorid.sector(); |
| 628 |
+ |
castorRecHit.depth[nhits] = castorid.module(); |
| 629 |
+ |
castorRecHit.phi[nhits] = getPhi(castorid); |
| 630 |
+ |
} |
| 631 |
+ |
|
| 632 |
+ |
nhits++; |
| 633 |
+ |
|
| 634 |
+ |
} // end loop castor rechits |
| 635 |
+ |
} |
| 636 |
+ |
|
| 637 |
+ |
castorRecHit.n = nhits; |
| 638 |
+ |
castorTree->Fill(); |
| 639 |
+ |
} |
| 640 |
+ |
|
| 641 |
+ |
if(doZDCRecHit_){ |
| 642 |
+ |
|
| 643 |
+ |
edm::Handle<ZDCRecHitCollection> zdcrechits; |
| 644 |
+ |
|
| 645 |
+ |
try{ ev.getByLabel("zdcreco",zdcrechits); } |
| 646 |
+ |
catch(...) { edm::LogWarning(" ZDC ") << " Cannot get ZDC RecHits " << std::endl; } |
| 647 |
+ |
|
| 648 |
+ |
int nhits = 0; |
| 649 |
+ |
|
| 650 |
+ |
if (zdcrechits.failedToGet()!=0 || !zdcrechits.isValid()) { |
| 651 |
+ |
edm::LogWarning(" ZDC ") << " Cannot read ZDCRecHitCollection" << std::endl; |
| 652 |
+ |
} else { |
| 653 |
+ |
for(size_t i1 = 0; i1 < zdcrechits->size(); ++i1) { |
| 654 |
+ |
const ZDCRecHit & rh = (*zdcrechits)[i1]; |
| 655 |
+ |
HcalZDCDetId zdcid = rh.id(); |
| 656 |
+ |
if (nhits < 18) { |
| 657 |
+ |
zdcRecHit.e[nhits] = rh.energy(); |
| 658 |
+ |
zdcRecHit.zside[nhits] = zdcid.zside(); |
| 659 |
+ |
zdcRecHit.section[nhits] = zdcid.section(); |
| 660 |
+ |
zdcRecHit.channel[nhits] = zdcid.channel(); |
| 661 |
+ |
} |
| 662 |
+ |
|
| 663 |
+ |
nhits++; |
| 664 |
+ |
|
| 665 |
+ |
} // end loop zdc rechits |
| 666 |
+ |
} |
| 667 |
+ |
|
| 668 |
+ |
zdcRecHit.n = nhits; |
| 669 |
+ |
zdcRecHitTree->Fill(); |
| 670 |
+ |
|
| 671 |
+ |
} |
| 672 |
+ |
|
| 673 |
+ |
if(doZDCDigi_){ |
| 674 |
+ |
|
| 675 |
+ |
edm::Handle<ZDCDigiCollection> zdcdigis; |
| 676 |
+ |
|
| 677 |
+ |
try{ ev.getByLabel("hcalDigis",zdcdigis); } |
| 678 |
+ |
catch(...) { edm::LogWarning(" ZDC ") << " Cannot get ZDC Digis " << std::endl; } |
| 679 |
+ |
|
| 680 |
+ |
int nhits = 0; |
| 681 |
+ |
|
| 682 |
+ |
if (zdcdigis.failedToGet()!=0 || !zdcdigis.isValid()) { |
| 683 |
+ |
edm::LogWarning(" ZDC ") << " Cannot read ZDCDigiCollection" << std::endl; |
| 684 |
+ |
} else { |
| 685 |
+ |
for(size_t i1 = 0; i1 < zdcdigis->size(); ++i1) { |
| 686 |
+ |
|
| 687 |
+ |
const ZDCDataFrame & rh = (*zdcdigis)[i1]; |
| 688 |
+ |
HcalZDCDetId zdcid = rh.id(); |
| 689 |
+ |
|
| 690 |
+ |
if (nhits < 18) { |
| 691 |
+ |
|
| 692 |
+ |
int ts = 0; |
| 693 |
+ |
|
| 694 |
+ |
for(int j1 = 0; j1 < rh.size(); j1++){ |
| 695 |
+ |
|
| 696 |
+ |
zdcDigi.chargefC[nhits][ts]= rh[ts].nominal_fC(); |
| 697 |
+ |
zdcDigi.adc[nhits][ts]= rh[ts].adc(); |
| 698 |
+ |
zdcDigi.zside[nhits][ts] = zdcid.zside(); |
| 699 |
+ |
zdcDigi.section[nhits][ts] = zdcid.section(); |
| 700 |
+ |
zdcDigi.channel[nhits][ts] = zdcid.channel(); |
| 701 |
+ |
zdcDigi.ts[nhits][ts] = ts; |
| 702 |
+ |
|
| 703 |
+ |
ts++; |
| 704 |
+ |
} |
| 705 |
+ |
|
| 706 |
+ |
zdcDigi.nts=ts; |
| 707 |
+ |
|
| 708 |
+ |
} |
| 709 |
+ |
|
| 710 |
+ |
nhits++; |
| 711 |
+ |
|
| 712 |
+ |
} // end loop zdc rechits |
| 713 |
+ |
} |
| 714 |
+ |
|
| 715 |
+ |
zdcDigi.n = nhits; |
| 716 |
+ |
zdcDigiTree->Fill(); |
| 717 |
+ |
|
| 718 |
+ |
} |
| 719 |
+ |
|
| 720 |
|
if(!doEbyEonly_){ |
| 721 |
|
towerTree->Fill(); |
| 722 |
|
|
| 731 |
|
} |
| 732 |
|
} |
| 733 |
|
|
| 734 |
< |
nt->Fill(nHFtower,nHFlong,nHFshort); |
| 734 |
> |
nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus); |
| 735 |
|
|
| 736 |
|
} |
| 737 |
|
|
| 741 |
|
RecHitTreeProducer::beginJob() |
| 742 |
|
{ |
| 743 |
|
|
| 744 |
< |
hbheTree = fs->make<TTree>("hbhe",""); |
| 744 |
> |
hbheTree = fs->make<TTree>("hbhe",versionTag); |
| 745 |
|
hbheTree->Branch("n",&hbheRecHit.n,"n/I"); |
| 746 |
|
hbheTree->Branch("e",hbheRecHit.e,"e[n]/F"); |
| 747 |
|
hbheTree->Branch("et",hbheRecHit.et,"et[n]/F"); |
| 748 |
|
hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F"); |
| 749 |
|
hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F"); |
| 750 |
+ |
hbheTree->Branch("perp",hbheRecHit.perp,"perp[n]/F"); |
| 751 |
+ |
|
| 752 |
|
hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O"); |
| 753 |
< |
|
| 754 |
< |
hfTree = fs->make<TTree>("hf",""); |
| 753 |
> |
hbheTree->Branch("depth",hfRecHit.depth,"depth[n]/I"); |
| 754 |
> |
|
| 755 |
> |
hfTree = fs->make<TTree>("hf",versionTag); |
| 756 |
|
hfTree->Branch("n",&hfRecHit.n,"n/I"); |
| 757 |
|
hfTree->Branch("e",hfRecHit.e,"e[n]/F"); |
| 758 |
|
hfTree->Branch("et",hfRecHit.et,"et[n]/F"); |
| 759 |
|
hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F"); |
| 760 |
|
hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F"); |
| 761 |
+ |
hfTree->Branch("perp",hfRecHit.perp,"perp[n]/F"); |
| 762 |
|
hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I"); |
| 763 |
|
hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O"); |
| 764 |
|
|
| 765 |
< |
eeTree = fs->make<TTree>("ee",""); |
| 765 |
> |
eeTree = fs->make<TTree>("ee",versionTag); |
| 766 |
|
eeTree->Branch("n",&eeRecHit.n,"n/I"); |
| 767 |
|
eeTree->Branch("e",eeRecHit.e,"e[n]/F"); |
| 768 |
|
eeTree->Branch("et",eeRecHit.et,"et[n]/F"); |
| 769 |
|
eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F"); |
| 770 |
|
eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F"); |
| 771 |
+ |
eeTree->Branch("perp",eeRecHit.perp,"perp[n]/F"); |
| 772 |
+ |
|
| 773 |
|
eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O"); |
| 774 |
|
|
| 775 |
< |
ebTree = fs->make<TTree>("eb",""); |
| 775 |
> |
ebTree = fs->make<TTree>("eb",versionTag); |
| 776 |
|
ebTree->Branch("n",&ebRecHit.n,"n/I"); |
| 777 |
|
ebTree->Branch("e",ebRecHit.e,"e[n]/F"); |
| 778 |
|
ebTree->Branch("et",ebRecHit.et,"et[n]/F"); |
| 779 |
|
ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F"); |
| 780 |
|
ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F"); |
| 781 |
+ |
ebTree->Branch("perp",ebRecHit.perp,"perp[n]/F"); |
| 782 |
+ |
|
| 783 |
|
ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O"); |
| 784 |
|
|
| 785 |
< |
towerTree = fs->make<TTree>("tower",""); |
| 785 |
> |
towerTree = fs->make<TTree>("tower",versionTag); |
| 786 |
|
towerTree->Branch("n",&myTowers.n,"n/I"); |
| 787 |
|
towerTree->Branch("e",myTowers.e,"e[n]/F"); |
| 788 |
|
towerTree->Branch("et",myTowers.et,"et[n]/F"); |
| 789 |
|
towerTree->Branch("eta",myTowers.eta,"eta[n]/F"); |
| 790 |
|
towerTree->Branch("phi",myTowers.phi,"phi[n]/F"); |
| 791 |
|
towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O"); |
| 792 |
+ |
towerTree->Branch("emEt",myTowers.emEt,"emEt[n]/F"); |
| 793 |
+ |
towerTree->Branch("hadEt",myTowers.hadEt,"hadEt[n]/F"); |
| 794 |
+ |
|
| 795 |
+ |
|
| 796 |
+ |
if(doCastor_){ |
| 797 |
+ |
castorTree = fs->make<TTree>("castor",versionTag); |
| 798 |
+ |
castorTree->Branch("n",&castorRecHit.n,"n/I"); |
| 799 |
+ |
castorTree->Branch("e",castorRecHit.e,"e[n]/F"); |
| 800 |
+ |
castorTree->Branch("iphi",castorRecHit.iphi,"iphi[n]/I"); |
| 801 |
+ |
castorTree->Branch("phi",castorRecHit.phi,"phi[n]/F"); |
| 802 |
+ |
castorTree->Branch("depth",castorRecHit.depth,"depth[n]/I"); |
| 803 |
+ |
} |
| 804 |
|
|
| 805 |
+ |
if(doZDCRecHit_){ |
| 806 |
+ |
zdcRecHitTree = fs->make<TTree>("zdcrechit",versionTag); |
| 807 |
+ |
zdcRecHitTree->Branch("n",&zdcRecHit.n,"n/I"); |
| 808 |
+ |
zdcRecHitTree->Branch("e",zdcRecHit.e,"e[n]/F"); |
| 809 |
+ |
zdcRecHitTree->Branch("zside",zdcRecHit.zside,"zside[n]/I"); |
| 810 |
+ |
zdcRecHitTree->Branch("section",zdcRecHit.section,"section[n]/I"); |
| 811 |
+ |
zdcRecHitTree->Branch("channel",zdcRecHit.channel,"channel[n]/I"); |
| 812 |
+ |
} |
| 813 |
+ |
|
| 814 |
+ |
if(doZDCDigi_){ |
| 815 |
+ |
zdcDigiTree = fs->make<TTree>("zdcrechit",versionTag); |
| 816 |
+ |
zdcDigiTree->Branch("n",&zdcDigi.n,"n/I"); |
| 817 |
+ |
zdcDigiTree->Branch("nts",&zdcDigi.nts,"nts/I"); |
| 818 |
+ |
zdcDigiTree->Branch("zside",zdcDigi.zside,"zside[n][nts]/I"); |
| 819 |
+ |
zdcDigiTree->Branch("section",zdcDigi.section,"section[n][nts]/I"); |
| 820 |
+ |
zdcDigiTree->Branch("channel",zdcDigi.channel,"channel[n][nts]/I"); |
| 821 |
+ |
zdcDigiTree->Branch("adc",zdcDigi.adc,"adc[n][nts]/I"); |
| 822 |
+ |
zdcDigiTree->Branch("chargefC",zdcDigi.chargefC,"chargefC[n][nts]/F"); |
| 823 |
+ |
} |
| 824 |
+ |
|
| 825 |
+ |
if (saveBothVtx_) { |
| 826 |
+ |
towerTree->Branch("etVtx",myTowers.etVtx,"etvtx[n]/F"); |
| 827 |
+ |
towerTree->Branch("etaVtx",myTowers.etaVtx,"etavtx[n]/F"); |
| 828 |
+ |
towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F"); |
| 829 |
+ |
towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F"); |
| 830 |
+ |
} |
| 831 |
|
|
| 832 |
|
if(doBasicClusters_){ |
| 833 |
< |
bcTree = fs->make<TTree>("bc",""); |
| 833 |
> |
bcTree = fs->make<TTree>("bc",versionTag); |
| 834 |
|
bcTree->Branch("n",&myBC.n,"n/I"); |
| 835 |
|
bcTree->Branch("e",myBC.e,"e[n]/F"); |
| 836 |
|
bcTree->Branch("et",myBC.et,"et[n]/F"); |
| 843 |
|
} |
| 844 |
|
|
| 845 |
|
if(doFastJet_){ |
| 846 |
< |
bkgTree = fs->make<TTree>("bkg",""); |
| 846 |
> |
bkgTree = fs->make<TTree>("bkg",versionTag); |
| 847 |
|
bkgTree->Branch("n",&bkg.n,"n/I"); |
| 848 |
|
bkgTree->Branch("rho",bkg.rho,"rho[n]/F"); |
| 849 |
|
bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F"); |
| 850 |
|
} |
| 851 |
|
|
| 852 |
< |
nt = fs->make<TNtuple>("ntEvent","","nHF:nHFlong:nHFshort"); |
| 852 |
> |
nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus"); |
| 853 |
|
|
| 854 |
|
} |
| 855 |
|
|
| 858 |
|
RecHitTreeProducer::endJob() { |
| 859 |
|
} |
| 860 |
|
|
| 861 |
< |
double RecHitTreeProducer::getEt(const DetId &id, double energy){ |
| 861 |
> |
math::XYZPoint RecHitTreeProducer::getPosition(const DetId &id, reco::Vertex::Point vtx){ |
| 862 |
> |
const GlobalPoint& pos=geo->getPosition(id); |
| 863 |
> |
math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z()); |
| 864 |
> |
return posV; |
| 865 |
> |
} |
| 866 |
> |
|
| 867 |
> |
double RecHitTreeProducer::getEt(math::XYZPoint pos, double energy){ |
| 868 |
> |
double et = energy*sin(pos.theta()); |
| 869 |
> |
return et; |
| 870 |
> |
} |
| 871 |
> |
|
| 872 |
> |
double RecHitTreeProducer::getEt(const DetId &id, double energy, reco::Vertex::Point vtx){ |
| 873 |
|
const GlobalPoint& pos=geo->getPosition(id); |
| 874 |
|
double et = energy*sin(pos.theta()); |
| 875 |
|
return et; |
| 876 |
|
} |
| 877 |
|
|
| 878 |
< |
double RecHitTreeProducer::getEta(const DetId &id){ |
| 878 |
> |
double RecHitTreeProducer::getEta(const DetId &id, reco::Vertex::Point vtx){ |
| 879 |
|
const GlobalPoint& pos=geo->getPosition(id); |
| 880 |
|
double et = pos.eta(); |
| 881 |
|
return et; |
| 882 |
|
} |
| 883 |
|
|
| 884 |
< |
double RecHitTreeProducer::getPhi(const DetId &id){ |
| 884 |
> |
double RecHitTreeProducer::getPhi(const DetId &id, reco::Vertex::Point vtx){ |
| 885 |
|
const GlobalPoint& pos=geo->getPosition(id); |
| 886 |
|
double et = pos.phi(); |
| 887 |
|
return et; |
| 888 |
|
} |
| 889 |
|
|
| 890 |
+ |
double RecHitTreeProducer::getPerp(const DetId &id, reco::Vertex::Point vtx){ |
| 891 |
+ |
const GlobalPoint& pos=geo->getPosition(id); |
| 892 |
+ |
double et = pos.perp(); |
| 893 |
+ |
return et; |
| 894 |
+ |
} |
| 895 |
|
|
| 896 |
+ |
reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev) |
| 897 |
+ |
{ |
| 898 |
+ |
ev.getByLabel(VtxSrc_,vtxs); |
| 899 |
+ |
int greatestvtx = 0; |
| 900 |
+ |
int nVertex = vtxs->size(); |
| 901 |
+ |
|
| 902 |
+ |
for (unsigned int i = 0 ; i< vtxs->size(); ++i){ |
| 903 |
+ |
unsigned int daughter = (*vtxs)[i].tracksSize(); |
| 904 |
+ |
if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i; |
| 905 |
+ |
//cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl; |
| 906 |
+ |
} |
| 907 |
+ |
|
| 908 |
+ |
if(nVertex<=0){ |
| 909 |
+ |
return reco::Vertex::Point(0,0,0); |
| 910 |
+ |
} |
| 911 |
+ |
return (*vtxs)[greatestvtx].position(); |
| 912 |
+ |
} |
| 913 |
|
|
| 914 |
|
//define this as a plug-in |
| 915 |
|
DEFINE_FWK_MODULE(RecHitTreeProducer); |
