File indexing completed on 2023-03-17 11:27:13
0001
0002
0003
0004
0005
0006
0007
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "Validation/EcalHits/interface/EcalEndcapSimHitsValidation.h"
0011 #include <DataFormats/EcalDetId/interface/EEDetId.h>
0012
0013 using namespace cms;
0014 using namespace edm;
0015 using namespace std;
0016
0017 EcalEndcapSimHitsValidation::EcalEndcapSimHitsValidation(const edm::ParameterSet &ps)
0018 : g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
0019 EEHitsCollection(ps.getParameter<std::string>("EEHitsCollection")),
0020 ValidationCollection(ps.getParameter<std::string>("ValidationCollection")) {
0021 EEHitsToken =
0022 consumes<edm::PCaloHitContainer>(edm::InputTag(std::string(g4InfoLabel), std::string(EEHitsCollection)));
0023 ValidationCollectionToken =
0024 consumes<PEcalValidInfo>(edm::InputTag(std::string(g4InfoLabel), std::string(ValidationCollection)));
0025
0026 verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0027
0028 myEntries = 0;
0029 for (int myStep = 0; myStep < 26; myStep++) {
0030 eRLength[myStep] = 0.0;
0031 }
0032 }
0033
0034 void EcalEndcapSimHitsValidation::bookHistograms(DQMStore::IBooker &ib, edm::Run const &, edm::EventSetup const &c) {
0035 ib.setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
0036 ib.setScope(MonitorElementData::Scope::RUN);
0037
0038 std::string histo = "EE+ hits multiplicity";
0039 meEEzpHits_ = ib.book1D(histo, histo, 50, 0., 5000.);
0040
0041 histo = "EE- hits multiplicity";
0042 meEEzmHits_ = ib.book1D(histo, histo, 50, 0., 5000.);
0043
0044 histo = "EE+ crystals multiplicity";
0045 meEEzpCrystals_ = ib.book1D(histo, histo, 200, 0., 2000.);
0046
0047 histo = "EE- crystals multiplicity";
0048 meEEzmCrystals_ = ib.book1D(histo, histo, 200, 0., 2000.);
0049
0050 histo = "EE+ occupancy";
0051 meEEzpOccupancy_ = ib.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
0052
0053 histo = "EE- occupancy";
0054 meEEzmOccupancy_ = ib.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
0055
0056 histo = "EE longitudinal shower profile";
0057 meEELongitudinalShower_ = ib.bookProfile(histo, histo, 26, 0, 26, 100, 0, 20000);
0058
0059 histo = "EE hits energy spectrum";
0060 meEEHitEnergy_ = ib.book1D(histo, histo, 4000, 0., 400.);
0061
0062 histo = "EE hits log10energy spectrum";
0063 meEEhitLog10Energy_ = ib.book1D(histo, histo, 140, -10., 4.);
0064
0065 histo = "EE hits log10energy spectrum vs normalized energy";
0066 meEEhitLog10EnergyNorm_ = ib.bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
0067
0068 histo = "EE hits log10energy spectrum vs normalized energy25";
0069 meEEhitLog10Energy25Norm_ = ib.bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
0070
0071 histo = "EE hits energy spectrum 2";
0072 meEEHitEnergy2_ = ib.book1D(histo, histo, 1000, 0., 0.001);
0073
0074 histo = "EE crystal energy spectrum";
0075 meEEcrystalEnergy_ = ib.book1D(histo, histo, 5000, 0., 50.);
0076
0077 histo = "EE crystal energy spectrum 2";
0078 meEEcrystalEnergy2_ = ib.book1D(histo, histo, 1000, 0., 0.001);
0079
0080 histo = "EE E1";
0081 meEEe1_ = ib.book1D(histo, histo, 400, 0., 400.);
0082
0083 histo = "EE E4";
0084 meEEe4_ = ib.book1D(histo, histo, 400, 0., 400.);
0085
0086 histo = "EE E9";
0087 meEEe9_ = ib.book1D(histo, histo, 400, 0., 400.);
0088
0089 histo = "EE E16";
0090 meEEe16_ = ib.book1D(histo, histo, 400, 0., 400.);
0091
0092 histo = "EE E25";
0093 meEEe25_ = ib.book1D(histo, histo, 400, 0., 400.);
0094
0095 histo = "EE E1oE4";
0096 meEEe1oe4_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0097
0098 histo = "EE E1oE9";
0099 meEEe1oe9_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0100
0101 histo = "EE E4oE9";
0102 meEEe4oe9_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0103
0104 histo = "EE E9oE16";
0105 meEEe9oe16_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0106
0107 histo = "EE E1oE25";
0108 meEEe1oe25_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0109
0110 histo = "EE E9oE25";
0111 meEEe9oe25_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0112
0113 histo = "EE E16oE25";
0114 meEEe16oe25_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0115 }
0116
0117 void EcalEndcapSimHitsValidation::analyze(const edm::Event &e, const edm::EventSetup &c) {
0118 edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
0119
0120 edm::Handle<edm::PCaloHitContainer> EcalHitsEE;
0121 e.getByToken(EEHitsToken, EcalHitsEE);
0122
0123
0124 if (!EcalHitsEE.isValid())
0125 return;
0126
0127 edm::Handle<PEcalValidInfo> MyPEcalValidInfo;
0128 e.getByToken(ValidationCollectionToken, MyPEcalValidInfo);
0129
0130 std::vector<PCaloHit> theEECaloHits;
0131 theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
0132
0133 myEntries++;
0134
0135 std::map<unsigned int, std::vector<PCaloHit *>, std::less<unsigned int>> CaloHitMap;
0136
0137 double EEetzp_ = 0.;
0138 double EEetzm_ = 0.;
0139
0140 double ee1 = 0.0;
0141 double ee4 = 0.0;
0142 double ee9 = 0.0;
0143 double ee16 = 0.0;
0144 double ee25 = 0.0;
0145 std::vector<double> econtr(140, 0.);
0146 std::vector<double> econtr25(140, 0.);
0147
0148 MapType eemap;
0149 MapType eemapzp;
0150 MapType eemapzm;
0151 uint32_t nEEzpHits = 0;
0152 uint32_t nEEzmHits = 0;
0153
0154 for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin(); isim != theEECaloHits.end(); ++isim) {
0155 if (isim->time() > 500.) {
0156 continue;
0157 }
0158
0159 CaloHitMap[isim->id()].push_back(&(*isim));
0160
0161 EEDetId eeid(isim->id());
0162
0163 LogDebug("HitInfo") << " CaloHit " << isim->getName() << "\n"
0164 << " DetID = " << isim->id() << " EEDetId = " << eeid.ix() << " " << eeid.iy() << "\n"
0165 << " Time = " << isim->time() << "\n"
0166 << " Track Id = " << isim->geantTrackId() << "\n"
0167 << " Energy = " << isim->energy();
0168
0169 uint32_t crystid = eeid.rawId();
0170
0171 if (eeid.zside() > 0) {
0172 nEEzpHits++;
0173 EEetzp_ += isim->energy();
0174 eemapzp[crystid] += isim->energy();
0175 meEEzpOccupancy_->Fill(eeid.ix(), eeid.iy());
0176 } else if (eeid.zside() < 0) {
0177 nEEzmHits++;
0178 EEetzm_ += isim->energy();
0179 eemapzm[crystid] += isim->energy();
0180 meEEzmOccupancy_->Fill(eeid.ix(), eeid.iy());
0181 }
0182
0183 meEEHitEnergy_->Fill(isim->energy());
0184 if (isim->energy() > 0) {
0185 meEEhitLog10Energy_->Fill(log10(isim->energy()));
0186 int log10i = int((log10(isim->energy()) + 10.) * 10.);
0187 if (log10i >= 0 && log10i < 140)
0188 econtr[log10i] += isim->energy();
0189 }
0190 meEEHitEnergy2_->Fill(isim->energy());
0191 eemap[crystid] += isim->energy();
0192 }
0193
0194 meEEzpCrystals_->Fill(eemapzp.size());
0195 meEEzmCrystals_->Fill(eemapzm.size());
0196
0197 for (std::map<uint32_t, float, std::less<uint32_t>>::iterator it = eemap.begin(); it != eemap.end(); ++it)
0198 meEEcrystalEnergy_->Fill((*it).second);
0199 for (std::map<uint32_t, float, std::less<uint32_t>>::iterator it = eemap.begin(); it != eemap.end(); ++it)
0200 meEEcrystalEnergy2_->Fill((*it).second);
0201
0202 meEEzpHits_->Fill(nEEzpHits);
0203 meEEzmHits_->Fill(nEEzmHits);
0204
0205 int nEEHits = nEEzmHits + nEEzpHits;
0206 if (nEEHits > 0) {
0207 uint32_t eecenterid = getUnitWithMaxEnergy(eemap);
0208 EEDetId myEEid(eecenterid);
0209 int bx = myEEid.ix();
0210 int by = myEEid.iy();
0211 int bz = myEEid.zside();
0212 ee1 = energyInMatrixEE(1, 1, bx, by, bz, eemap);
0213 meEEe1_->Fill(ee1);
0214 ee9 = energyInMatrixEE(3, 3, bx, by, bz, eemap);
0215 meEEe9_->Fill(ee9);
0216 ee25 = energyInMatrixEE(5, 5, bx, by, bz, eemap);
0217 meEEe25_->Fill(ee25);
0218
0219 std::vector<uint32_t> ids25;
0220 ids25 = getIdsAroundMax(5, 5, bx, by, bz, eemap);
0221
0222 for (unsigned i = 0; i < 25; i++) {
0223 for (unsigned int j = 0; j < CaloHitMap[ids25[i]].size(); j++) {
0224 if (CaloHitMap[ids25[i]][j]->energy() > 0) {
0225 int log10i = int((log10(CaloHitMap[ids25[i]][j]->energy()) + 10.) * 10.);
0226 if (log10i >= 0 && log10i < 140)
0227 econtr25[log10i] += CaloHitMap[ids25[i]][j]->energy();
0228 }
0229 }
0230 }
0231
0232 MapType neweemap;
0233 if (fillEEMatrix(3, 3, bx, by, bz, neweemap, eemap)) {
0234 ee4 = eCluster2x2(neweemap);
0235 meEEe4_->Fill(ee4);
0236 }
0237 if (fillEEMatrix(5, 5, bx, by, bz, neweemap, eemap)) {
0238 ee16 = eCluster4x4(ee9, neweemap);
0239 meEEe16_->Fill(ee16);
0240 }
0241
0242 if (ee4 > 0.1)
0243 meEEe1oe4_->Fill(ee1 / ee4);
0244 if (ee9 > 0.1)
0245 meEEe1oe9_->Fill(ee1 / ee9);
0246 if (ee9 > 0.1)
0247 meEEe4oe9_->Fill(ee4 / ee9);
0248 if (ee16 > 0.1)
0249 meEEe9oe16_->Fill(ee9 / ee16);
0250 if (ee25 > 0.1)
0251 meEEe16oe25_->Fill(ee16 / ee25);
0252 if (ee25 > 0.1)
0253 meEEe1oe25_->Fill(ee1 / ee25);
0254 if (ee25 > 0.1)
0255 meEEe9oe25_->Fill(ee9 / ee25);
0256
0257 if ((EEetzp_ + EEetzm_) != 0) {
0258 for (int i = 0; i < 140; i++) {
0259 meEEhitLog10EnergyNorm_->Fill(-10. + (float(i) + 0.5) / 10., econtr[i] / (EEetzp_ + EEetzm_));
0260 }
0261 }
0262
0263 if (ee25 != 0) {
0264 for (int i = 0; i < 140; i++) {
0265 meEEhitLog10Energy25Norm_->Fill(-10. + (float(i) + 0.5) / 10., econtr25[i] / ee25);
0266 }
0267 }
0268 }
0269
0270 if (MyPEcalValidInfo.isValid()) {
0271 if (MyPEcalValidInfo->ee1x1() > 0.) {
0272 std::vector<float> EX0 = MyPEcalValidInfo->eX0();
0273 meEELongitudinalShower_->Reset();
0274 for (int myStep = 0; myStep < 26; myStep++) {
0275 eRLength[myStep] += EX0[myStep];
0276 meEELongitudinalShower_->Fill(float(myStep), eRLength[myStep] / myEntries);
0277 }
0278 }
0279 }
0280 }
0281
0282 float EcalEndcapSimHitsValidation::energyInMatrixEE(
0283 int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap) {
0284 int ncristals = 0;
0285 float totalEnergy = 0.;
0286
0287 int goBackInX = nCellInX / 2;
0288 int goBackInY = nCellInY / 2;
0289 int startX = centralX - goBackInX;
0290 int startY = centralY - goBackInY;
0291
0292 for (int ix = startX; ix < startX + nCellInX; ix++) {
0293 for (int iy = startY; iy < startY + nCellInY; iy++) {
0294 uint32_t index;
0295
0296 if (EEDetId::validDetId(ix, iy, centralZ)) {
0297 index = EEDetId(ix, iy, centralZ).rawId();
0298 } else {
0299 continue;
0300 }
0301
0302 totalEnergy += themap[index];
0303 ncristals += 1;
0304 }
0305 }
0306
0307 LogDebug("GeomInfo") << nCellInX << " x " << nCellInY << " EE matrix energy = " << totalEnergy << " for " << ncristals
0308 << " crystals";
0309 return totalEnergy;
0310 }
0311
0312 std::vector<uint32_t> EcalEndcapSimHitsValidation::getIdsAroundMax(
0313 int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap) {
0314 int ncristals = 0;
0315 std::vector<uint32_t> ids(nCellInX * nCellInY);
0316
0317 int goBackInX = nCellInX / 2;
0318 int goBackInY = nCellInY / 2;
0319 int startX = centralX - goBackInX;
0320 int startY = centralY - goBackInY;
0321
0322 for (int ix = startX; ix < startX + nCellInX; ix++) {
0323 for (int iy = startY; iy < startY + nCellInY; iy++) {
0324 uint32_t index;
0325
0326 if (EEDetId::validDetId(ix, iy, centralZ)) {
0327 index = EEDetId(ix, iy, centralZ).rawId();
0328 } else {
0329 continue;
0330 }
0331
0332 ids[ncristals] = index;
0333 ncristals += 1;
0334 }
0335 }
0336
0337 return ids;
0338 }
0339
0340 bool EcalEndcapSimHitsValidation::fillEEMatrix(
0341 int nCellInX, int nCellInY, int CentralX, int CentralY, int CentralZ, MapType &fillmap, MapType &themap) {
0342 int goBackInX = nCellInX / 2;
0343 int goBackInY = nCellInY / 2;
0344
0345 int startX = CentralX - goBackInX;
0346 int startY = CentralY - goBackInY;
0347
0348 int i = 0;
0349 for (int ix = startX; ix < startX + nCellInX; ix++) {
0350 for (int iy = startY; iy < startY + nCellInY; iy++) {
0351 uint32_t index;
0352
0353 if (EEDetId::validDetId(ix, iy, CentralZ)) {
0354 index = EEDetId(ix, iy, CentralZ).rawId();
0355 } else {
0356 continue;
0357 }
0358 fillmap[i++] = themap[index];
0359 }
0360 }
0361 uint32_t centerid = getUnitWithMaxEnergy(themap);
0362
0363 if (fillmap[i / 2] == themap[centerid])
0364 return true;
0365 else
0366 return false;
0367 }
0368
0369 float EcalEndcapSimHitsValidation::eCluster2x2(MapType &themap) {
0370 float E22 = 0.;
0371 float e012 = themap[0] + themap[1] + themap[2];
0372 float e036 = themap[0] + themap[3] + themap[6];
0373 float e678 = themap[6] + themap[7] + themap[8];
0374 float e258 = themap[2] + themap[5] + themap[8];
0375
0376 if ((e012 > e678 || e012 == e678) && (e036 > e258 || e036 == e258))
0377 E22 = themap[0] + themap[1] + themap[3] + themap[4];
0378 else if ((e012 > e678 || e012 == e678) && (e036 < e258 || e036 == e258))
0379 E22 = themap[1] + themap[2] + themap[4] + themap[5];
0380 else if ((e012 < e678 || e012 == e678) && (e036 > e258 || e036 == e258))
0381 E22 = themap[3] + themap[4] + themap[6] + themap[7];
0382 else if ((e012 < e678 || e012 == e678) && (e036 < e258 || e036 == e258))
0383 E22 = themap[4] + themap[5] + themap[7] + themap[8];
0384
0385 return E22;
0386 }
0387
0388 float EcalEndcapSimHitsValidation::eCluster4x4(float e33, MapType &themap) {
0389 float E44 = 0.;
0390 float e0_4 = themap[0] + themap[1] + themap[2] + themap[3] + themap[4];
0391 float e0_20 = themap[0] + themap[5] + themap[10] + themap[15] + themap[20];
0392 float e4_24 = themap[4] + themap[9] + themap[14] + themap[19] + themap[24];
0393 float e0_24 = themap[20] + themap[21] + themap[22] + themap[23] + themap[24];
0394
0395 if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
0396 E44 = e33 + themap[0] + themap[1] + themap[2] + themap[3] + themap[5] + themap[10] + themap[15];
0397 else if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
0398 E44 = e33 + themap[1] + themap[2] + themap[3] + themap[4] + themap[9] + themap[14] + themap[19];
0399 else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
0400 E44 = e33 + themap[5] + themap[10] + themap[15] + themap[20] + themap[21] + themap[22] + themap[23];
0401 else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
0402 E44 = e33 + themap[21] + themap[22] + themap[23] + themap[24] + themap[9] + themap[14] + themap[19];
0403
0404 return E44;
0405 }
0406
0407 uint32_t EcalEndcapSimHitsValidation::getUnitWithMaxEnergy(MapType &themap) {
0408
0409 uint32_t unitWithMaxEnergy = 0;
0410 float maxEnergy = 0.;
0411
0412 MapType::iterator iter;
0413 for (iter = themap.begin(); iter != themap.end(); iter++) {
0414 if (maxEnergy < (*iter).second) {
0415 maxEnergy = (*iter).second;
0416 unitWithMaxEnergy = (*iter).first;
0417 }
0418 }
0419
0420 LogDebug("GeomInfo") << " max energy of " << maxEnergy << " GeV in Unit id " << unitWithMaxEnergy;
0421 return unitWithMaxEnergy;
0422 }