File indexing completed on 2024-04-06 12:25:27
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 #include "RecoJets/JetAnalyzers/interface/JetToDigiDump.h"
0016 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0017 #include "DataFormats/JetReco/interface/CaloJet.h"
0018
0019
0020
0021
0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0024 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0025 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0026 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0027 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0028 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0029 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0030 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0031 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include <TROOT.h>
0035 #include <TSystem.h>
0036 #include <TFile.h>
0037 #include <TCanvas.h>
0038 #include <cmath>
0039 using namespace edm;
0040 using namespace reco;
0041 using namespace std;
0042
0043 JetToDigiDump::JetToDigiDump(const ParameterSet& cfg)
0044 : DumpLevel(cfg.getParameter<string>("DumpLevel")),
0045 CaloJetAlg(cfg.getParameter<string>("CaloJetAlg")),
0046 DebugLevel(cfg.getParameter<int>("DebugLevel")),
0047 ShowECal(cfg.getParameter<bool>("ShowECal")) {}
0048
0049 void JetToDigiDump::beginJob() {
0050 if (DumpLevel == "Jets") {
0051 cout << "Dump of Jets" << endl;
0052 Dump = 1;
0053 } else if (DumpLevel == "Towers") {
0054 cout << "Dump of Jets and constituent CaloTowers" << endl;
0055 Dump = 2;
0056 } else if (DumpLevel == "RecHits") {
0057 cout << "Dump of Jets, constituent CaloTowers, and constituent RecHits" << endl;
0058 Dump = 3;
0059 } else if (DumpLevel == "Digis") {
0060 cout << "Dump of Jets, constituent CaloTowers, constituent RecHits and associated Digis" << endl;
0061 Dump = 4;
0062 }
0063 cout << "Jet Algorithm being dumped is " << CaloJetAlg << endl;
0064 cout << "Debug level is " << DebugLevel << endl;
0065
0066 evtCount = 0;
0067 }
0068
0069 void JetToDigiDump::analyze(const Event& evt, const EventSetup& es) {
0070 int jetInd;
0071 Handle<CaloJetCollection> caloJets;
0072 Handle<CaloTowerCollection> caloTowers;
0073 Handle<HBHERecHitCollection> HBHERecHits;
0074 Handle<HORecHitCollection> HORecHits;
0075 Handle<HFRecHitCollection> HFRecHits;
0076 Handle<EBRecHitCollection> EBRecHits;
0077 Handle<EERecHitCollection> EERecHits;
0078 Handle<HBHEDigiCollection> HBHEDigis;
0079 Handle<HODigiCollection> HODigis;
0080 Handle<HFDigiCollection> HFDigis;
0081 Handle<EEDigiCollection> EEDigis;
0082 Handle<EBDigiCollection> EBDigis;
0083
0084
0085
0086
0087
0088 if (DebugLevel)
0089 cout << "Getting caloJets" << endl;
0090
0091 evt.getByLabel(CaloJetAlg, caloJets);
0092 if (Dump >= 2)
0093 evt.getByLabel("towerMaker", caloTowers);
0094 if (Dump >= 3) {
0095 if (DebugLevel)
0096 cout << "Getting recHits" << endl;
0097 evt.getByLabel("hbhereco", HBHERecHits);
0098 evt.getByLabel("horeco", HORecHits);
0099 evt.getByLabel("hfreco", HFRecHits);
0100 evt.getByLabel("ecalRecHit", "EcalRecHitsEB", EBRecHits);
0101 evt.getByLabel("ecalRecHit", "EcalRecHitsEE", EERecHits);
0102 if (DebugLevel)
0103 cout << "# of hits gotten - HBHE: " << HBHERecHits->size() << endl;
0104 evt.getByLabel("hcalDigis", HBHEDigis);
0105 evt.getByLabel("hcalDigis", HODigis);
0106 evt.getByLabel("hcalDigis", HFDigis);
0107 if (DebugLevel)
0108 cout << "# of digis gotten - HBHE: " << HBHEDigis->size() << endl;
0109 if (ShowECal) {
0110 evt.getByLabel("ecalDigis", "ebDigis", EBDigis);
0111 evt.getByLabel("ecalDigis", "eeDigis", EEDigis);
0112 }
0113 }
0114
0115 cout << endl << "Evt: " << evtCount << ", Num Jets=" << caloJets->end() - caloJets->begin() << endl;
0116 if (Dump >= 1)
0117 cout << " *********************************************************" << endl;
0118 jetInd = 0;
0119 if (Dump >= 1)
0120 for (CaloJetCollection::const_iterator jet = caloJets->begin(); jet != caloJets->end(); ++jet) {
0121
0122
0123 std::vector<CaloTowerPtr> towers = jet->getCaloConstituents();
0124 int nConstituents = towers.size();
0125 cout << " Jet: " << jetInd << ", eta=" << jet->eta() << ", phi=" << jet->phi() << ", pt=" << jet->pt()
0126 << ",E=" << jet->energy() << ", EB E=" << jet->emEnergyInEB() << " ,HB E=" << jet->hadEnergyInHB()
0127 << ", HO E=" << jet->hadEnergyInHO() << " ,EE E=" << jet->emEnergyInEE() << ", HE E=" << jet->hadEnergyInHE()
0128 << ", HF E=" << jet->hadEnergyInHF() + jet->emEnergyInHF() << ", Num Towers=" << nConstituents << endl;
0129 if (Dump >= 2)
0130 cout << " =====================================================" << endl;
0131 float sumTowerE = 0.0;
0132 if (Dump >= 2)
0133 for (int i = 0; i < nConstituents; i++) {
0134 CaloTowerCollection::const_iterator theTower =
0135 caloTowers->find(towers[i]->id());
0136 if (theTower == caloTowers->end()) {
0137 cerr << "Bug? Can't find the tower" << endl;
0138 return;
0139 }
0140 int ietaTower = towers[i]->id().ieta();
0141 int iphiTower = towers[i]->id().iphi();
0142 sumTowerE += theTower->energy();
0143 size_t numRecHits = theTower->constituentsSize();
0144 cout << " Tower " << i << ": ieta=" << ietaTower << ", eta=" << theTower->eta() << ", iphi=" << iphiTower
0145 << ", phi=" << theTower->phi() << ", energy=" << theTower->energy() << ", EM=" << theTower->emEnergy()
0146 << ", HAD=" << theTower->hadEnergy() << ", HO=" << theTower->outerEnergy()
0147 << ", Num Rec Hits =" << numRecHits << endl;
0148 if (Dump >= 3)
0149 cout << " ------------------------------------------------" << endl;
0150 float sumRecHitE = 0.0;
0151 if (Dump >= 3)
0152 for (size_t j = 0; j < numRecHits; j++) {
0153 DetId RecHitDetID = theTower->constituent(j);
0154 DetId::Detector DetNum = RecHitDetID.det();
0155 if (DetNum == DetId::Hcal) {
0156
0157 HcalDetId HcalID = RecHitDetID;
0158 HcalSubdetector HcalNum = HcalID.subdet();
0159 if (HcalNum == HcalBarrel) {
0160 HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
0161 sumRecHitE += theRecHit->energy();
0162 HBHEDigiCollection::const_iterator theDigis = HBHEDigis->find(HcalID);
0163 cout << " RecHit: " << j << ": HB, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0164 << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0165 << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0166 << ", presamples =" << theDigis->presamples() << endl;
0167
0168
0169
0170
0171
0172 float SumDigiCharge = 0.0;
0173 float EstimatedPedestal = 0.0;
0174 int SamplesToAdd = 4;
0175 if (Dump >= 4)
0176 cout << " ......................................" << endl;
0177 if (Dump >= 4)
0178 for (int k = 0; k < theDigis->size(); k++) {
0179 const HcalQIESample QIE = theDigis->sample(k);
0180 if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0181 SumDigiCharge += QIE.nominal_fC();
0182 if (k < theDigis->presamples() - 1)
0183 EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0184 cout << " Digi: " << k << ", cap ID = " << QIE.capid()
0185 << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0186 }
0187 if (Dump >= 4)
0188 cout << " 4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0189 << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0190 if (Dump >= 4)
0191 cout << " ......................................" << endl;
0192 } else if (HcalNum == HcalEndcap) {
0193 HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
0194 if ((abs(HcalID.ieta()) == 28 || abs(HcalID.ieta()) == 29) && HcalID.depth() == 3) {
0195 sumRecHitE += theRecHit->energy() / 2;
0196 } else {
0197 sumRecHitE += theRecHit->energy();
0198 }
0199 HBHEDigiCollection::const_iterator theDigis = HBHEDigis->find(HcalID);
0200 cout << " RecHit: " << j << ": HE, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0201 << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0202 << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0203 << ", presamples =" << theDigis->presamples() << endl;
0204 float SumDigiCharge = 0.0;
0205 float EstimatedPedestal = 0.0;
0206 int SamplesToAdd = 4;
0207 if (Dump >= 4)
0208 cout << " ......................................" << endl;
0209 if (Dump >= 4)
0210 for (int k = 0; k < theDigis->size(); k++) {
0211 const HcalQIESample QIE = theDigis->sample(k);
0212 if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0213 SumDigiCharge += QIE.nominal_fC();
0214 if (k < theDigis->presamples() - 1)
0215 EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0216 cout << " Digi: " << k << ", cap ID = " << QIE.capid()
0217 << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0218 }
0219 if (Dump >= 4)
0220 cout << " 4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0221 << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0222 if (Dump >= 4)
0223 cout << " ......................................" << endl;
0224 } else if (HcalNum == HcalOuter) {
0225 HORecHitCollection::const_iterator theRecHit = HORecHits->find(HcalID);
0226 sumRecHitE += theRecHit->energy();
0227 HODigiCollection::const_iterator theDigis = HODigis->find(HcalID);
0228 cout << " RecHit: " << j << ": HO, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0229 << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0230 << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0231 << ", presamples =" << theDigis->presamples() << endl;
0232 float SumDigiCharge = 0.0;
0233 float EstimatedPedestal = 0.0;
0234 int SamplesToAdd = 4;
0235 if (Dump >= 4)
0236 cout << " ......................................" << endl;
0237 if (Dump >= 4)
0238 for (int k = 0; k < theDigis->size(); k++) {
0239 const HcalQIESample QIE = theDigis->sample(k);
0240 if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0241 SumDigiCharge += QIE.nominal_fC();
0242 if (k < theDigis->presamples() - 1)
0243 EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0244 cout << " Digi: " << k << ", cap ID = " << QIE.capid()
0245 << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0246 }
0247 if (Dump >= 4)
0248 cout << " 4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0249 << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0250 if (Dump >= 4)
0251 cout << " ......................................" << endl;
0252 } else if (HcalNum == HcalForward) {
0253 HFRecHitCollection::const_iterator theRecHit = HFRecHits->find(HcalID);
0254 sumRecHitE += theRecHit->energy();
0255 HFDigiCollection::const_iterator theDigis = HFDigis->find(HcalID);
0256 cout << " RecHit: " << j << ": HF, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0257 << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0258 << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0259 << ", presamples =" << theDigis->presamples() << endl;
0260 float SumDigiCharge = 0.0;
0261 float EstimatedPedestal = 0.0;
0262 int SamplesToAdd = 1;
0263 if (Dump >= 4)
0264 cout << " ......................................" << endl;
0265 if (Dump >= 4)
0266 for (int k = 0; k < theDigis->size(); k++) {
0267 const HcalQIESample QIE = theDigis->sample(k);
0268 if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0269 SumDigiCharge += QIE.nominal_fC();
0270 if (k < theDigis->presamples() - 1)
0271 EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0272 cout << " Digi: " << k << ", cap ID = " << QIE.capid()
0273 << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0274 }
0275 if (Dump >= 4)
0276 cout << " 1 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0277 << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0278 if (Dump >= 4)
0279 cout << " ......................................" << endl;
0280 }
0281 }
0282 if (ShowECal && DetNum == DetId::Ecal) {
0283 int EcalNum = RecHitDetID.subdetId();
0284 if (EcalNum == 1) {
0285 EBDetId EcalID = RecHitDetID;
0286 EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(EcalID);
0287 EBDigiCollection::const_iterator theDigis = EBDigis->find(EcalID);
0288 sumRecHitE += theRecHit->energy();
0289 cout << " RecHit " << j << ": EB, ieta=" << EcalID.ieta() << ", iphi=" << EcalID.iphi()
0290 << ", SM=" << EcalID.ism() << ", energy=" << theRecHit->energy()
0291 << ", All Digis=" << theDigis->size() << endl;
0292 if (Dump >= 4)
0293 cout << " ......................................" << endl;
0294 if (Dump >= 4)
0295 for (unsigned int k = 0; k < theDigis->size(); k++) {
0296 EBDataFrame frame(*theDigis);
0297 const EcalMGPASample MGPA = frame.sample(k);
0298 cout << " Digi: " << k << ": ADC Sample = " << MGPA.adc()
0299 << ", Gain ID = " << MGPA.gainId() << endl;
0300 }
0301 if (Dump >= 4)
0302 cout << " ......................................" << endl;
0303 } else if (EcalNum == 2) {
0304 EEDetId EcalID = RecHitDetID;
0305 EERecHitCollection::const_iterator theRecHit = EERecHits->find(EcalID);
0306 EEDigiCollection::const_iterator theDigis = EEDigis->find(EcalID);
0307 sumRecHitE += theRecHit->energy();
0308 cout << " RecHit " << j << ": EE, ix=" << EcalID.ix() << ", iy=" << EcalID.iy()
0309 << ", energy=" << theRecHit->energy() << ", All Digis=" << theDigis->size() << endl;
0310 if (Dump >= 4)
0311 cout << " ......................................" << endl;
0312 if (Dump >= 4)
0313 for (unsigned int k = 0; k < theDigis->size(); k++) {
0314 EEDataFrame frame(*theDigis);
0315 const EcalMGPASample MGPA = frame.sample(k);
0316 cout << " Digi: " << k << ": ADC Sample = " << MGPA.adc()
0317 << ", Gain ID = " << MGPA.gainId() << endl;
0318 }
0319 if (Dump >= 4)
0320 cout << " ......................................" << endl;
0321 }
0322 }
0323 }
0324 if (Dump >= 3) {
0325 if (abs(ietaTower) == 28 || abs(ietaTower) == 29) {
0326 cout << " Splitted Sum of RecHit Energies=" << sumRecHitE
0327 << ", CaloTower energy=" << theTower->energy() << endl;
0328 } else {
0329 cout << " Sum of RecHit Energies=" << sumRecHitE << ", CaloTower energy=" << theTower->energy()
0330 << endl;
0331 }
0332 }
0333 if (Dump >= 3)
0334 cout << " ------------------------------------------------" << endl;
0335 }
0336 if (Dump >= 2)
0337 cout << " Sum of tower energies=" << sumTowerE << ", CaloJet energy=" << jet->energy() << endl;
0338 jetInd++;
0339 if (Dump >= 2)
0340 cout << " =====================================================" << endl;
0341 }
0342 evtCount++;
0343 if (Dump >= 1)
0344 cout << " *********************************************************" << endl;
0345 }
0346
0347 void JetToDigiDump::endJob() {}
0348 #include "FWCore/Framework/interface/MakerMacros.h"
0349 DEFINE_FWK_MODULE(JetToDigiDump);