Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:27:13 

0001 /*
0002  * \file EcalEndcapSimHitsValidation.cc
0003  *
0004  * \author C.Rovelli
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   // verbosity switch
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   // Do nothing if no EndCap data available
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   // look for max
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!