00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "RecoMET/METProducers/interface/HcalNoiseInfoProducer.h"
00011 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00012 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00013 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00014 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00015 #include "DataFormats/TrackReco/interface/Track.h"
00016 #include "FWCore/Framework/interface/ESHandle.h"
00017 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
00018 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
00019 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
00020 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
00021 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
00022 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
00023 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
00024 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00025 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00026
00027 using namespace reco;
00028
00029
00030
00031
00032
00033 HcalNoiseInfoProducer::HcalNoiseInfoProducer(const edm::ParameterSet& iConfig) : algo_(iConfig)
00034 {
00035
00036 fillDigis_ = iConfig.getParameter<bool>("fillDigis");
00037 fillRecHits_ = iConfig.getParameter<bool>("fillRecHits");
00038 fillCaloTowers_ = iConfig.getParameter<bool>("fillCaloTowers");
00039 fillTracks_ = iConfig.getParameter<bool>("fillTracks");
00040
00041 maxProblemRBXs_ = iConfig.getParameter<int>("maxProblemRBXs");
00042
00043 maxCaloTowerIEta_ = iConfig.getParameter<int>("maxCaloTowerIEta");
00044 maxTrackEta_ = iConfig.getParameter<double>("maxTrackEta");
00045 minTrackPt_ = iConfig.getParameter<double>("minTrackPt");
00046
00047 digiCollName_ = iConfig.getParameter<std::string>("digiCollName");
00048 recHitCollName_ = iConfig.getParameter<std::string>("recHitCollName");
00049 caloTowerCollName_ = iConfig.getParameter<std::string>("caloTowerCollName");
00050 trackCollName_ = iConfig.getParameter<std::string>("trackCollName");
00051
00052 minRecHitE_ = iConfig.getParameter<double>("minRecHitE");
00053 minLowHitE_ = iConfig.getParameter<double>("minLowHitE");
00054 minHighHitE_ = iConfig.getParameter<double>("minHighHitE");
00055
00056 HcalAcceptSeverityLevel_ = iConfig.getParameter<uint32_t>("HcalAcceptSeverityLevel");
00057 if (iConfig.exists("HcalRecHitFlagsToBeExcluded"))
00058 HcalRecHitFlagsToBeExcluded_ = iConfig.getParameter<std::vector<int> >("HcalRecHitFlagsToBeExcluded");
00059 else{
00060 edm::LogWarning("MisConfiguration")<<"the module is missing the parameter HcalAcceptSeverityLevel. created empty.";
00061 HcalRecHitFlagsToBeExcluded_.resize(0);
00062 }
00063
00064 TS4TS5EnergyThreshold_ = iConfig.getParameter<double>("TS4TS5EnergyThreshold");
00065
00066 std::vector<double> TS4TS5UpperThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperThreshold");
00067 std::vector<double> TS4TS5UpperCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperCut");
00068 std::vector<double> TS4TS5LowerThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerThreshold");
00069 std::vector<double> TS4TS5LowerCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerCut");
00070
00071 for(int i = 0; i < (int)TS4TS5UpperThresholdTemp.size() && i < (int)TS4TS5UpperCutTemp.size(); i++)
00072 TS4TS5UpperCut_.push_back(std::pair<double, double>(TS4TS5UpperThresholdTemp[i], TS4TS5UpperCutTemp[i]));
00073 sort(TS4TS5UpperCut_.begin(), TS4TS5UpperCut_.end());
00074
00075 for(int i = 0; i < (int)TS4TS5LowerThresholdTemp.size() && i < (int)TS4TS5LowerCutTemp.size(); i++)
00076 TS4TS5LowerCut_.push_back(std::pair<double, double>(TS4TS5LowerThresholdTemp[i], TS4TS5LowerCutTemp[i]));
00077 sort(TS4TS5LowerCut_.begin(), TS4TS5LowerCut_.end());
00078
00079
00080 if(fillDigis_ && !fillRecHits_) {
00081 fillRecHits_=true;
00082 edm::LogWarning("HCalNoiseInfoProducer") << " forcing fillRecHits to be true if fillDigis is true.\n";
00083 }
00084
00085 const float adc2fCTemp[128]={-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,
00086 13.5,15.,17.,19.,21.,23.,25.,27.,29.5,32.5,35.5,38.5,42.,46.,50.,54.5,59.5,
00087 64.5,59.5,64.5,69.5,74.5,79.5,84.5,89.5,94.5,99.5,104.5,109.5,114.5,119.5,
00088 124.5,129.5,137.,147.,157.,167.,177.,187.,197.,209.5,224.5,239.5,254.5,272.,
00089 292.,312.,334.5,359.5,384.5,359.5,384.5,409.5,434.5,459.5,484.5,509.5,534.5,
00090 559.5,584.5,609.5,634.5,659.5,684.5,709.5,747.,797.,847.,897.,947.,997.,
00091 1047.,1109.5,1184.5,1259.5,1334.5,1422.,1522.,1622.,1734.5,1859.5,1984.5,
00092 1859.5,1984.5,2109.5,2234.5,2359.5,2484.5,2609.5,2734.5,2859.5,2984.5,
00093 3109.5,3234.5,3359.5,3484.5,3609.5,3797.,4047.,4297.,4547.,4797.,5047.,
00094 5297.,5609.5,5984.5,6359.5,6734.5,7172.,7672.,8172.,8734.5,9359.5,9984.5};
00095 for(int i = 0; i < 128; i++)
00096 adc2fC[i] = adc2fCTemp[i];
00097
00098
00099 produces<HcalNoiseRBXCollection>();
00100
00101 produces<HcalNoiseSummary>();
00102 }
00103
00104
00105 HcalNoiseInfoProducer::~HcalNoiseInfoProducer()
00106 {
00107 }
00108
00109
00110
00111
00112
00113
00114
00115 void
00116 HcalNoiseInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00117 {
00118
00119 std::auto_ptr<HcalNoiseRBXCollection> result1(new HcalNoiseRBXCollection);
00120 std::auto_ptr<HcalNoiseSummary> result2(new HcalNoiseSummary);
00121
00122
00123 HcalNoiseRBXArray rbxarray;
00124 HcalNoiseSummary &summary=*result2;
00125
00126
00127
00128 if(fillRecHits_) fillrechits(iEvent, iSetup, rbxarray, summary);
00129 if(fillDigis_) filldigis(iEvent, iSetup, rbxarray, summary);
00130 if(fillCaloTowers_) fillcalotwrs(iEvent, iSetup, rbxarray, summary);
00131 if(fillTracks_) filltracks(iEvent, iSetup, summary);
00132
00133
00134
00135
00136
00137
00138 HcalNoiseRBXArray::iterator maxit=rbxarray.begin();
00139 double maxenergy=-999;
00140 bool maxwritten=false;
00141 for(HcalNoiseRBXArray::iterator rit = rbxarray.begin(); rit!=rbxarray.end(); ++rit) {
00142 HcalNoiseRBX &rbx=(*rit);
00143 CommonHcalNoiseRBXData data(rbx, minRecHitE_, minLowHitE_, minHighHitE_, TS4TS5EnergyThreshold_,
00144 TS4TS5UpperCut_, TS4TS5LowerCut_);
00145
00146
00147 if(data.energy()>maxenergy) {
00148 maxenergy=data.energy();
00149 maxit=rit;
00150 maxwritten=false;
00151 }
00152
00153
00154 bool writerbx = algo_.isProblematic(data) || !algo_.passLooseNoiseFilter(data) ||
00155 !algo_.passTightNoiseFilter(data) || !algo_.passHighLevelNoiseFilter(data);
00156
00157
00158 fillOtherSummaryVariables(summary, data);
00159
00160 if(writerbx) {
00161 summary.nproblemRBXs_++;
00162 if(summary.nproblemRBXs_<=maxProblemRBXs_) {
00163 result1->push_back(rbx);
00164 if(maxit==rit) maxwritten=true;
00165 }
00166 }
00167 }
00168
00169
00170 if(!maxwritten) {
00171 HcalNoiseRBX &rbx=(*maxit);
00172
00173
00174 result1->push_back(rbx);
00175 }
00176
00177
00178 iEvent.put(result1);
00179 iEvent.put(result2);
00180
00181 return;
00182 }
00183
00184
00185 void
00186 HcalNoiseInfoProducer::beginJob()
00187 {
00188 return;
00189 }
00190
00191
00192 void
00193 HcalNoiseInfoProducer::endJob()
00194 {
00195 return;
00196 }
00197
00198
00199
00200
00201 void
00202 HcalNoiseInfoProducer::beginRun(edm::Run&, const edm::EventSetup&)
00203 {
00204 return;
00205 }
00206
00207
00208 void
00209 HcalNoiseInfoProducer::endRun(edm::Run&, const edm::EventSetup&)
00210 {
00211 return;
00212 }
00213
00214
00215 void
00216 HcalNoiseInfoProducer::fillOtherSummaryVariables(HcalNoiseSummary& summary, const CommonHcalNoiseRBXData& data) const
00217 {
00218
00219 if(algo_.passRatioThreshold(data) && data.validRatio()) {
00220 if(data.ratio()<summary.minE2Over10TS()) {
00221 summary.mine2ts_ = data.e2ts();
00222 summary.mine10ts_ = data.e10ts(); }
00223 if(data.ratio()>summary.maxE2Over10TS()) {
00224 summary.maxe2ts_ = data.e2ts();
00225 summary.maxe10ts_ = data.e10ts();
00226 }
00227 }
00228
00229
00230 if(algo_.passZerosThreshold(data)) {
00231 if(data.numZeros()>summary.maxZeros()) {
00232 summary.maxzeros_ = data.numZeros();
00233 }
00234 }
00235
00236
00237 if(data.numHPDHits() > summary.maxHPDHits()) {
00238 summary.maxhpdhits_ = data.numHPDHits();
00239 }
00240 if(data.numRBXHits() > summary.maxRBXHits()) {
00241 summary.maxrbxhits_ = data.numRBXHits();
00242 }
00243 if(data.numHPDNoOtherHits() > summary.maxHPDNoOtherHits()) {
00244 summary.maxhpdhitsnoother_ = data.numHPDNoOtherHits();
00245 }
00246
00247
00248 if(data.PassTS4TS5() == false)
00249 summary.hasBadRBXTS4TS5_ = true;
00250
00251
00252 if(data.minLowEHitTime()<summary.min10GeVHitTime()) {
00253 summary.min10_ = data.minLowEHitTime();
00254 }
00255 if(data.maxLowEHitTime()>summary.max10GeVHitTime()) {
00256 summary.max10_ = data.maxLowEHitTime();
00257 }
00258 summary.rms10_ += data.lowEHitTimeSqrd();
00259 summary.cnthit10_ += data.numLowEHits();
00260 if(data.minHighEHitTime()<summary.min25GeVHitTime()) {
00261 summary.min25_ = data.minHighEHitTime();
00262 }
00263 if(data.maxHighEHitTime()>summary.max25GeVHitTime()) {
00264 summary.max25_ = data.maxHighEHitTime();
00265 }
00266 summary.rms25_ += data.highEHitTimeSqrd();
00267 summary.cnthit25_ += data.numHighEHits();
00268
00269
00270 if(algo_.passEMFThreshold(data)) {
00271 if(summary.minHPDEMF() > data.HPDEMF()) {
00272 summary.minhpdemf_ = data.HPDEMF();
00273 }
00274 if(summary.minRBXEMF() > data.RBXEMF()) {
00275 summary.minrbxemf_ = data.RBXEMF();
00276 }
00277 }
00278
00279
00280 if(!algo_.passLooseRatio(data)) summary.filterstatus_ |= 0x1;
00281 if(!algo_.passLooseHits(data)) summary.filterstatus_ |= 0x2;
00282 if(!algo_.passLooseZeros(data)) summary.filterstatus_ |= 0x4;
00283 if(!algo_.passLooseTiming(data)) summary.filterstatus_ |= 0x8;
00284
00285 if(!algo_.passTightRatio(data)) summary.filterstatus_ |= 0x100;
00286 if(!algo_.passTightHits(data)) summary.filterstatus_ |= 0x200;
00287 if(!algo_.passTightZeros(data)) summary.filterstatus_ |= 0x400;
00288 if(!algo_.passTightTiming(data)) summary.filterstatus_ |= 0x800;
00289
00290 if(!algo_.passHighLevelNoiseFilter(data)) summary.filterstatus_ |= 0x10000;
00291
00292
00293 JoinCaloTowerRefVectorsWithoutDuplicates join;
00294 if(!algo_.passLooseNoiseFilter(data))
00295 join(summary.loosenoisetwrs_, data.rbxTowers());
00296 if(!algo_.passTightNoiseFilter(data))
00297 join(summary.tightnoisetwrs_, data.rbxTowers());
00298 if(!algo_.passHighLevelNoiseFilter(data))
00299 join(summary.hlnoisetwrs_, data.rbxTowers());
00300
00301 return;
00302 }
00303
00304
00305
00306 void
00307 HcalNoiseInfoProducer::filldigis(edm::Event& iEvent, const edm::EventSetup& iSetup, HcalNoiseRBXArray& array, HcalNoiseSummary& summary)
00308 {
00309
00310 TotalCalibCharge = 0;
00311 int NcalibTS45=0;
00312 int NcalibTS45gt15=0;
00313 double chargecalibTS45=0;
00314 double chargecalibgt15TS45=0;
00315
00316
00317 edm::ESHandle<HcalDbService> conditions;
00318 iSetup.get<HcalDbRecord>().get(conditions);
00319 edm::ESHandle<HcalChannelQuality> qualhandle;
00320 iSetup.get<HcalChannelQualityRcd>().get(qualhandle);
00321 const HcalChannelQuality* myqual = qualhandle.product();
00322 edm::ESHandle<HcalSeverityLevelComputer> mycomputer;
00323 iSetup.get<HcalSeverityLevelComputerRcd>().get(mycomputer);
00324 const HcalSeverityLevelComputer* mySeverity = mycomputer.product();
00325
00326
00327 edm::Handle<HBHEDigiCollection> handle;
00328 iEvent.getByLabel(digiCollName_, handle);
00329 if(!handle.isValid()) {
00330 throw edm::Exception(edm::errors::ProductNotFound) << " could not find HBHEDigiCollection named " << digiCollName_ << "\n.";
00331 return;
00332 }
00333
00334
00335 for(HBHEDigiCollection::const_iterator it=handle->begin(); it!=handle->end(); ++it) {
00336 const HBHEDataFrame &digi=(*it);
00337 HcalDetId cell = digi.id();
00338 DetId detcell=(DetId)cell;
00339
00340
00341 const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00342 if(mySeverity->dropChannel(mydigistatus->getValue())) continue;
00343 if(digi.zsMarkAndPass()) continue;
00344
00345 if ((mydigistatus->getValue() & (1 <<HcalChannelStatus::HcalCellExcludeFromHBHENoiseSummary))==1) continue;
00346
00347
00348 const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00349 const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00350 const HcalQIEShape* shape = conditions->getHcalShape(channelCoder);
00351 HcalCoderDb coder (*channelCoder, *shape);
00352
00353
00354 HcalNoiseRBX &rbx=(*array.findRBX(digi));
00355 HcalNoiseHPD &hpd=(*array.findHPD(digi));
00356
00357
00358
00359 bool isBig=false, isBig5=false, isRBX=false;
00360 int counter=0;
00361 edm::RefVector<HBHERecHitCollection> &rechits=hpd.rechits_;
00362 for(edm::RefVector<HBHERecHitCollection>::const_iterator rit=rechits.begin();
00363 rit!=rechits.end(); ++rit, ++counter) {
00364 if((*rit)->id() == digi.id()) {
00365 if(counter==0) isBig=isBig5=true;
00366 if(counter<5) isBig5=true;
00367 isRBX=true;
00368 }
00369 }
00370
00371
00372 int totalzeros=0;
00373 CaloSamples tool;
00374 coder.adc2fC(digi,tool);
00375 for(int ts=0; ts<tool.size(); ++ts) {
00376
00377
00378 if(digi[ts].adc()==0) {
00379 ++hpd.totalZeros_;
00380 ++totalzeros;
00381 }
00382
00383
00384 double corrfc = tool[ts]-calibrations.pedestal(digi[ts].capid());
00385
00386
00387 if(isBig) hpd.bigCharge_[ts]+=corrfc;
00388 if(isBig5) hpd.big5Charge_[ts]+=corrfc;
00389 if(isRBX) rbx.allCharge_[ts]+=corrfc;
00390 }
00391
00392
00393 if(totalzeros>hpd.maxZeros_)
00394 hpd.maxZeros_=totalzeros;
00395 }
00396
00397
00398 edm::Handle<HcalCalibDigiCollection> hCalib;
00399 iEvent.getByLabel("hcalDigis", hCalib);
00400
00401
00402 if(hCalib.isValid() == true)
00403 {
00404 for(HcalCalibDigiCollection::const_iterator digi = hCalib->begin(); digi != hCalib->end(); digi++)
00405 {
00406 if(digi->id().hcalSubdet() == 0)
00407 continue;
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 for(int i = 0; i < (int)digi->size(); i++)
00432 TotalCalibCharge = TotalCalibCharge + adc2fC[digi->sample(i).adc()&0xff];
00433
00434
00435 HcalCalibDetId myid=(HcalCalibDetId)digi->id();
00436 if ( myid.calibFlavor()==HcalCalibDetId::HOCrosstalk)
00437 continue;
00438 if (digi->size()>5)
00439 {
00440 double sumCharge=adc2fC[digi->sample(4).adc()&0xff]+adc2fC[digi->sample(5).adc()&0xff];
00441
00442 ++NcalibTS45;
00443 chargecalibTS45+=sumCharge;
00444 if (sumCharge>15)
00445 {
00446 ++NcalibTS45gt15;
00447 chargecalibgt15TS45+=sumCharge;
00448 }
00449 }
00450
00451 }
00452 }
00453
00454 summary.calibCharge_ = TotalCalibCharge;
00455 summary.calibCountTS45_=NcalibTS45;
00456 summary.calibCountgt15TS45_=NcalibTS45gt15;
00457 summary.calibChargeTS45_=chargecalibTS45;
00458 summary.calibChargegt15TS45_=chargecalibgt15TS45;
00459 return;
00460 }
00461
00462
00463 void
00464 HcalNoiseInfoProducer::fillrechits(edm::Event& iEvent, const edm::EventSetup& iSetup, HcalNoiseRBXArray& array, HcalNoiseSummary& summary) const
00465 {
00466
00467 edm::ESHandle<HcalChannelQuality> hcalChStatus;
00468 iSetup.get<HcalChannelQualityRcd>().get( hcalChStatus );
00469 const HcalChannelQuality* dbHcalChStatus = hcalChStatus.product();
00470
00471
00472 edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
00473 iSetup.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
00474 const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product();
00475
00476
00477 edm::ESHandle<CaloGeometry> pG;
00478 iSetup.get<CaloGeometryRecord>().get(pG);
00479 const CaloGeometry* geo = pG.product();
00480
00481
00482 edm::Handle<HBHERecHitCollection> handle;
00483 iEvent.getByLabel(recHitCollName_, handle);
00484 if(!handle.isValid()) {
00485 throw edm::Exception(edm::errors::ProductNotFound)
00486 << " could not find HBHERecHitCollection named " << recHitCollName_ << "\n.";
00487 return;
00488 }
00489
00490 summary.rechitCount_ = 0;
00491 summary.rechitCount15_ = 0;
00492 summary.rechitEnergy_ = 0;
00493 summary.rechitEnergy15_ = 0;
00494
00495 summary.hitsInLaserRegion_=0;
00496 summary.hitsInNonLaserRegion_=0;
00497 summary.energyInLaserRegion_=0;
00498 summary.energyInNonLaserRegion_=0;
00499
00500
00501
00502
00503 for(HBHERecHitCollection::const_iterator it=handle->begin(); it!=handle->end(); ++it) {
00504 const HBHERecHit &rechit=(*it);
00505
00506
00507 const DetId id = rechit.detid();
00508 uint32_t recHitFlag = rechit.flags();
00509 uint32_t isolbitset = (1 << HcalCaloFlagLabels::HBHEIsolatedNoise);
00510 uint32_t flatbitset = (1 << HcalCaloFlagLabels::HBHEFlatNoise);
00511 uint32_t spikebitset = (1 << HcalCaloFlagLabels::HBHESpikeNoise);
00512 uint32_t trianglebitset = (1 << HcalCaloFlagLabels::HBHETriangleNoise);
00513 uint32_t ts4ts5bitset = (1 << HcalCaloFlagLabels::HBHETS4TS5Noise);
00514 for(unsigned int i=0; i<HcalRecHitFlagsToBeExcluded_.size(); i++) {
00515 uint32_t bitset = (1 << HcalRecHitFlagsToBeExcluded_[i]);
00516 recHitFlag = (recHitFlag & bitset) ? recHitFlag-bitset : recHitFlag;
00517 }
00518 const uint32_t dbStatusFlag = dbHcalChStatus->getValues(id)->getValue();
00519
00520
00521 if ((dbStatusFlag & (1 <<HcalChannelStatus::HcalCellExcludeFromHBHENoiseSummary))==1) continue;
00522
00523 int severityLevel = hcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
00524 bool isRecovered = hcalSevLvlComputer->recoveredRecHit(id, recHitFlag);
00525 if(severityLevel!=0 && !isRecovered && severityLevel>static_cast<int>(HcalAcceptSeverityLevel_)) continue;
00526
00527
00528 summary.rechitCount_ = summary.rechitCount_ + 1;
00529 summary.rechitEnergy_ = summary.rechitEnergy_ + rechit.energy();
00530 if ((dbStatusFlag & (1 <<HcalChannelStatus::HcalBadLaserSignal))==1)
00531 {
00532 ++summary.hitsInNonLaserRegion_;
00533 summary.energyInNonLaserRegion_+=rechit.energy();
00534 }
00535 else
00536 {
00537 ++summary.hitsInLaserRegion_;
00538 summary.energyInLaserRegion_+=rechit.energy();
00539 }
00540
00541 if(rechit.energy() > 1.5)
00542 {
00543 summary.rechitCount15_ = summary.rechitCount15_ + 1;
00544 summary.rechitEnergy15_ = summary.rechitEnergy15_ + rechit.energy();
00545 }
00546
00547
00548 if(rechit.flags() & isolbitset) {
00549 summary.nisolnoise_++;
00550 summary.isolnoisee_ += rechit.energy();
00551 GlobalPoint gp = geo->getPosition(rechit.id());
00552 double et = rechit.energy()*gp.perp()/gp.mag();
00553 summary.isolnoiseet_ += et;
00554 }
00555
00556 if(rechit.flags() & flatbitset) {
00557 summary.nflatnoise_++;
00558 summary.flatnoisee_ += rechit.energy();
00559 GlobalPoint gp = geo->getPosition(rechit.id());
00560 double et = rechit.energy()*gp.perp()/gp.mag();
00561 summary.flatnoiseet_ += et;
00562 }
00563
00564 if(rechit.flags() & spikebitset) {
00565 summary.nspikenoise_++;
00566 summary.spikenoisee_ += rechit.energy();
00567 GlobalPoint gp = geo->getPosition(rechit.id());
00568 double et = rechit.energy()*gp.perp()/gp.mag();
00569 summary.spikenoiseet_ += et;
00570 }
00571
00572 if(rechit.flags() & trianglebitset) {
00573 summary.ntrianglenoise_++;
00574 summary.trianglenoisee_ += rechit.energy();
00575 GlobalPoint gp = geo->getPosition(rechit.id());
00576 double et = rechit.energy()*gp.perp()/gp.mag();
00577 summary.trianglenoiseet_ += et;
00578 }
00579
00580 if(rechit.flags() & ts4ts5bitset) {
00581 if ((dbStatusFlag & (1 <<HcalChannelStatus::HcalCellExcludeFromHBHENoiseSummaryR45))==0)
00582 {
00583 summary.nts4ts5noise_++;
00584 summary.ts4ts5noisee_ += rechit.energy();
00585 GlobalPoint gp = geo->getPosition(rechit.id());
00586 double et = rechit.energy()*gp.perp()/gp.mag();
00587 summary.ts4ts5noiseet_ += et;
00588 }
00589 }
00590
00591
00592 HcalNoiseHPD& hpd=(*array.findHPD(rechit));
00593
00594
00595 edm::Ref<HBHERecHitCollection> myRef(handle, it-handle->begin());
00596
00597
00598 hpd.refrechitset_.insert(myRef);
00599
00600 }
00601
00602
00603 for(HcalNoiseRBXArray::iterator rbxit=array.begin(); rbxit!=array.end(); ++rbxit) {
00604 for(std::vector<HcalNoiseHPD>::iterator hpdit=rbxit->hpds_.begin(); hpdit!=rbxit->hpds_.end(); ++hpdit) {
00605
00606
00607 for(std::set<edm::Ref<HBHERecHitCollection>, RefHBHERecHitEnergyComparison>::const_iterator
00608 it=hpdit->refrechitset_.begin(); it!=hpdit->refrechitset_.end(); ++it) {
00609 hpdit->rechits_.push_back(*it);
00610 }
00611 }
00612 }
00613
00614
00615 return;
00616 }
00617
00618
00619 void
00620 HcalNoiseInfoProducer::fillcalotwrs(edm::Event& iEvent, const edm::EventSetup& iSetup, HcalNoiseRBXArray& array, HcalNoiseSummary& summary) const
00621 {
00622
00623 edm::Handle<CaloTowerCollection> handle;
00624 iEvent.getByLabel(caloTowerCollName_, handle);
00625 if(!handle.isValid()) {
00626 throw edm::Exception(edm::errors::ProductNotFound)
00627 << " could not find CaloTowerCollection named " << caloTowerCollName_ << "\n.";
00628 return;
00629 }
00630
00631 summary.emenergy_ = summary.hadenergy_ = 0.0;
00632
00633
00634 for(CaloTowerCollection::const_iterator it = handle->begin(); it!=handle->end(); ++it) {
00635 const CaloTower& twr=(*it);
00636
00637
00638 edm::Ref<CaloTowerCollection> myRef(handle, it-handle->begin());
00639
00640
00641 std::vector<std::vector<HcalNoiseHPD>::iterator> hpditervec;
00642 array.findHPD(twr, hpditervec);
00643
00644
00645 for(std::vector<std::vector<HcalNoiseHPD>::iterator>::iterator it=hpditervec.begin();
00646 it!=hpditervec.end(); ++it)
00647 (*it)->calotowers_.push_back(myRef);
00648
00649
00650 if(twr.ietaAbs()>maxCaloTowerIEta_) {
00651 summary.emenergy_ += twr.emEnergy();
00652 summary.hadenergy_ += twr.hadEnergy();
00653 }
00654 }
00655
00656 return;
00657 }
00658
00659
00660 void
00661 HcalNoiseInfoProducer::filltracks(edm::Event& iEvent, const edm::EventSetup& iSetup, HcalNoiseSummary& summary) const
00662 {
00663 edm::Handle<reco::TrackCollection> handle;
00664 iEvent.getByLabel(trackCollName_, handle);
00665
00666
00667 if(!handle.isValid()) {
00668
00669
00670 return;
00671 }
00672
00673 summary.trackenergy_=0.0;
00674 for(reco::TrackCollection::const_iterator iTrack = handle->begin(); iTrack!=handle->end(); ++iTrack) {
00675 reco::Track trk=*iTrack;
00676 if(trk.pt()<minTrackPt_ || fabs(trk.eta())>maxTrackEta_) continue;
00677
00678 summary.trackenergy_ += trk.p();
00679 }
00680
00681 return;
00682 }
00683
00684
00685 DEFINE_FWK_MODULE(HcalNoiseInfoProducer);