Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:06 

0001 /*
0002  * \file EcalBarrelSimHitsValidation.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/EcalBarrelSimHitsValidation.h"
0011 #include <DataFormats/EcalDetId/interface/EBDetId.h>
0012 
0013 using namespace cms;
0014 using namespace edm;
0015 using namespace std;
0016 
0017 EcalBarrelSimHitsValidation::EcalBarrelSimHitsValidation(const edm::ParameterSet &ps)
0018     : g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
0019       EBHitsCollection(ps.getParameter<std::string>("EBHitsCollection")),
0020       ValidationCollection(ps.getParameter<std::string>("ValidationCollection")) {
0021   EBHitsToken =
0022       consumes<edm::PCaloHitContainer>(edm::InputTag(std::string(g4InfoLabel), std::string(EBHitsCollection)));
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 EcalBarrelSimHitsValidation::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 = "EB hits multiplicity";
0039   menEBHits_ = ib.book1D(histo, histo, 50, 0., 5000.);
0040 
0041   histo = "EB crystals multiplicity";
0042   menEBCrystals_ = ib.book1D(histo, histo, 200, 0., 2000.);
0043 
0044   histo = "EB occupancy";
0045   meEBOccupancy_ = ib.book2D(histo, histo, 360, 0., 360., 170, -85., 85.);
0046 
0047   histo = "EB longitudinal shower profile";
0048   meEBLongitudinalShower_ = ib.bookProfile(histo, histo, 26, 0., 26., 100, 0., 20000.);
0049 
0050   histo = "EB hits energy spectrum";
0051   meEBhitEnergy_ = ib.book1D(histo, histo, 4000, 0., 400.);
0052 
0053   histo = "EB hits log10energy spectrum";
0054   meEBhitLog10Energy_ = ib.book1D(histo, histo, 140, -10., 4.);
0055 
0056   histo = "EB hits log10energy spectrum vs normalized energy";
0057   meEBhitLog10EnergyNorm_ = ib.bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
0058 
0059   histo = "EB hits log10energy spectrum vs normalized energy25";
0060   meEBhitLog10Energy25Norm_ = ib.bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
0061 
0062   histo = "EB hits energy spectrum 2";
0063   meEBhitEnergy2_ = ib.book1D(histo, histo, 1000, 0., 0.001);
0064 
0065   histo = "EB crystal energy spectrum";
0066   meEBcrystalEnergy_ = ib.book1D(histo, histo, 5000, 0., 50.);
0067 
0068   histo = "EB crystal energy spectrum 2";
0069   meEBcrystalEnergy2_ = ib.book1D(histo, histo, 1000, 0., 0.001);
0070 
0071   histo = "EB E1";
0072   meEBe1_ = ib.book1D(histo, histo, 400, 0., 400.);
0073 
0074   histo = "EB E4";
0075   meEBe4_ = ib.book1D(histo, histo, 400, 0., 400.);
0076 
0077   histo = "EB E9";
0078   meEBe9_ = ib.book1D(histo, histo, 400, 0., 400.);
0079 
0080   histo = "EB E16";
0081   meEBe16_ = ib.book1D(histo, histo, 400, 0., 400.);
0082 
0083   histo = "EB E25";
0084   meEBe25_ = ib.book1D(histo, histo, 400, 0., 400.);
0085 
0086   histo = "EB E1oE4";
0087   meEBe1oe4_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0088 
0089   histo = "EB E1oE9";
0090   meEBe1oe9_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0091 
0092   histo = "EB E4oE9";
0093   meEBe4oe9_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0094 
0095   histo = "EB E9oE16";
0096   meEBe9oe16_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0097 
0098   histo = "EB E1oE25";
0099   meEBe1oe25_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0100 
0101   histo = "EB E9oE25";
0102   meEBe9oe25_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0103 
0104   histo = "EB E16oE25";
0105   meEBe16oe25_ = ib.book1D(histo, histo, 100, 0.4, 1.1);
0106 }
0107 
0108 void EcalBarrelSimHitsValidation::analyze(const edm::Event &e, const edm::EventSetup &c) {
0109   edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
0110 
0111   edm::Handle<edm::PCaloHitContainer> EcalHitsEB;
0112   e.getByToken(EBHitsToken, EcalHitsEB);
0113 
0114   // Do nothing if no Barrel data available
0115   if (!EcalHitsEB.isValid())
0116     return;
0117 
0118   edm::Handle<PEcalValidInfo> MyPEcalValidInfo;
0119   e.getByToken(ValidationCollectionToken, MyPEcalValidInfo);
0120 
0121   std::vector<PCaloHit> theEBCaloHits;
0122   theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
0123 
0124   myEntries++;
0125 
0126   double EBEnergy_ = 0.;
0127   std::map<unsigned int, std::vector<PCaloHit *>, std::less<unsigned int>> CaloHitMap;
0128 
0129   double eb1 = 0.0;
0130   double eb4 = 0.0;
0131   double eb9 = 0.0;
0132   double eb16 = 0.0;
0133   double eb25 = 0.0;
0134   std::vector<double> econtr(140, 0.);
0135   std::vector<double> econtr25(140, 0.);
0136 
0137   MapType ebmap;
0138   uint32_t nEBHits = 0;
0139 
0140   for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin(); isim != theEBCaloHits.end(); ++isim) {
0141     if (isim->time() > 500.) {
0142       continue;
0143     }
0144 
0145     CaloHitMap[isim->id()].push_back(&(*isim));
0146 
0147     EBDetId ebid(isim->id());
0148 
0149     LogDebug("HitInfo") << " CaloHit " << isim->getName() << "\n"
0150                         << " DetID = " << isim->id() << " EBDetId = " << ebid.ieta() << " " << ebid.iphi() << "\n"
0151                         << " Time = " << isim->time() << "\n"
0152                         << " Track Id = " << isim->geantTrackId() << "\n"
0153                         << " Energy = " << isim->energy();
0154 
0155     meEBOccupancy_->Fill(ebid.iphi(), ebid.ieta());
0156 
0157     uint32_t crystid = ebid.rawId();
0158     ebmap[crystid] += isim->energy();
0159 
0160     EBEnergy_ += isim->energy();
0161     nEBHits++;
0162     meEBhitEnergy_->Fill(isim->energy());
0163     if (isim->energy() > 0) {
0164       meEBhitLog10Energy_->Fill(log10(isim->energy()));
0165       int log10i = int((log10(isim->energy()) + 10.) * 10.);
0166       if (log10i >= 0 && log10i < 140)
0167         econtr[log10i] += isim->energy();
0168     }
0169     meEBhitEnergy2_->Fill(isim->energy());
0170   }
0171 
0172   menEBCrystals_->Fill(ebmap.size());
0173   if (meEBcrystalEnergy_) {
0174     for (std::map<uint32_t, float, std::less<uint32_t>>::iterator it = ebmap.begin(); it != ebmap.end(); ++it)
0175       meEBcrystalEnergy_->Fill((*it).second);
0176   }
0177   if (meEBcrystalEnergy2_) {
0178     for (std::map<uint32_t, float, std::less<uint32_t>>::iterator it = ebmap.begin(); it != ebmap.end(); ++it)
0179       meEBcrystalEnergy2_->Fill((*it).second);
0180   }
0181 
0182   menEBHits_->Fill(nEBHits);
0183 
0184   if (nEBHits > 0) {
0185     uint32_t ebcenterid = getUnitWithMaxEnergy(ebmap);
0186     EBDetId myEBid(ebcenterid);
0187     int bx = myEBid.ietaAbs();
0188     int by = myEBid.iphi();
0189     int bz = myEBid.zside();
0190     eb1 = energyInMatrixEB(1, 1, bx, by, bz, ebmap);
0191     meEBe1_->Fill(eb1);
0192     eb9 = energyInMatrixEB(3, 3, bx, by, bz, ebmap);
0193     meEBe9_->Fill(eb9);
0194     eb25 = energyInMatrixEB(5, 5, bx, by, bz, ebmap);
0195     meEBe25_->Fill(eb25);
0196 
0197     std::vector<uint32_t> ids25;
0198     ids25 = getIdsAroundMax(5, 5, bx, by, bz, ebmap);
0199 
0200     for (unsigned i = 0; i < 25; i++) {
0201       for (unsigned int j = 0; j < CaloHitMap[ids25[i]].size(); j++) {
0202         if (CaloHitMap[ids25[i]][j]->energy() > 0) {
0203           int log10i = int((log10(CaloHitMap[ids25[i]][j]->energy()) + 10.) * 10.);
0204           if (log10i >= 0 && log10i < 140)
0205             econtr25[log10i] += CaloHitMap[ids25[i]][j]->energy();
0206         }
0207       }
0208     }
0209 
0210     MapType newebmap;
0211     if (fillEBMatrix(3, 3, bx, by, bz, newebmap, ebmap)) {
0212       eb4 = eCluster2x2(newebmap);
0213       meEBe4_->Fill(eb4);
0214     }
0215     if (fillEBMatrix(5, 5, bx, by, bz, newebmap, ebmap)) {
0216       eb16 = eCluster4x4(eb9, newebmap);
0217       meEBe16_->Fill(eb16);
0218     }
0219 
0220     if (eb4 > 0.1)
0221       meEBe1oe4_->Fill(eb1 / eb4);
0222     if (eb9 > 0.1)
0223       meEBe1oe9_->Fill(eb1 / eb9);
0224     if (eb9 > 0.1)
0225       meEBe4oe9_->Fill(eb4 / eb9);
0226     if (eb16 > 0.1)
0227       meEBe9oe16_->Fill(eb9 / eb16);
0228     if (eb25 > 0.1)
0229       meEBe1oe25_->Fill(eb1 / eb25);
0230     if (eb25 > 0.1)
0231       meEBe9oe25_->Fill(eb9 / eb25);
0232     if (eb25 > 0.1)
0233       meEBe16oe25_->Fill(eb16 / eb25);
0234 
0235     if (EBEnergy_ != 0) {
0236       for (int i = 0; i < 140; i++) {
0237         meEBhitLog10EnergyNorm_->Fill(-10. + (float(i) + 0.5) / 10., econtr[i] / EBEnergy_);
0238       }
0239     }
0240 
0241     if (eb25 != 0) {
0242       for (int i = 0; i < 140; i++) {
0243         meEBhitLog10Energy25Norm_->Fill(-10. + (float(i) + 0.5) / 10., econtr25[i] / eb25);
0244       }
0245     }
0246   }
0247 
0248   if (MyPEcalValidInfo.isValid()) {
0249     if (MyPEcalValidInfo->eb1x1() > 0.) {
0250       std::vector<float> BX0 = MyPEcalValidInfo->bX0();
0251       meEBLongitudinalShower_->Reset();
0252       for (int myStep = 0; myStep < 26; myStep++) {
0253         eRLength[myStep] += BX0[myStep];
0254         meEBLongitudinalShower_->Fill(float(myStep), eRLength[myStep] / myEntries);
0255       }
0256     }
0257   }
0258 }
0259 
0260 float EcalBarrelSimHitsValidation::energyInMatrixEB(
0261     int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType &themap) {
0262   int ncristals = 0;
0263   float totalEnergy = 0.;
0264 
0265   int goBackInEta = nCellInEta / 2;
0266   int goBackInPhi = nCellInPhi / 2;
0267   int startEta = centralZ * centralEta - goBackInEta;
0268   int startPhi = centralPhi - goBackInPhi;
0269 
0270   for (int ieta = startEta; ieta < startEta + nCellInEta; ieta++) {
0271     for (int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++) {
0272       uint32_t index;
0273       if (abs(ieta) > 85 || abs(ieta) < 1) {
0274         continue;
0275       }
0276       if (iphi < 1) {
0277         index = EBDetId(ieta, iphi + 360).rawId();
0278       } else if (iphi > 360) {
0279         index = EBDetId(ieta, iphi - 360).rawId();
0280       } else {
0281         index = EBDetId(ieta, iphi).rawId();
0282       }
0283 
0284       totalEnergy += themap[index];
0285       ncristals += 1;
0286     }
0287   }
0288 
0289   LogDebug("GeomInfo") << nCellInEta << " x " << nCellInPhi << " EB matrix energy = " << totalEnergy << " for "
0290                        << ncristals << " crystals";
0291   return totalEnergy;
0292 }
0293 
0294 std::vector<uint32_t> EcalBarrelSimHitsValidation::getIdsAroundMax(
0295     int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType &themap) {
0296   int ncristals = 0;
0297   std::vector<uint32_t> ids(nCellInEta * nCellInPhi);
0298 
0299   int goBackInEta = nCellInEta / 2;
0300   int goBackInPhi = nCellInPhi / 2;
0301   int startEta = centralZ * centralEta - goBackInEta;
0302   int startPhi = centralPhi - goBackInPhi;
0303 
0304   for (int ieta = startEta; ieta < startEta + nCellInEta; ieta++) {
0305     for (int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++) {
0306       uint32_t index;
0307       if (abs(ieta) > 85 || abs(ieta) < 1) {
0308         continue;
0309       }
0310       if (iphi < 1) {
0311         index = EBDetId(ieta, iphi + 360).rawId();
0312       } else if (iphi > 360) {
0313         index = EBDetId(ieta, iphi - 360).rawId();
0314       } else {
0315         index = EBDetId(ieta, iphi).rawId();
0316       }
0317       ids[ncristals] = index;
0318       ncristals += 1;
0319     }
0320   }
0321 
0322   return ids;
0323 }
0324 
0325 bool EcalBarrelSimHitsValidation::fillEBMatrix(
0326     int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &fillmap, MapType &themap) {
0327   int goBackInEta = nCellInEta / 2;
0328   int goBackInPhi = nCellInPhi / 2;
0329 
0330   int startEta = CentralZ * CentralEta - goBackInEta;
0331   int startPhi = CentralPhi - goBackInPhi;
0332 
0333   int i = 0;
0334   for (int ieta = startEta; ieta < startEta + nCellInEta; ieta++) {
0335     for (int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++) {
0336       uint32_t index;
0337       if (abs(ieta) > 85 || abs(ieta) < 1) {
0338         continue;
0339       }
0340       if (iphi < 1) {
0341         index = EBDetId(ieta, iphi + 360).rawId();
0342       } else if (iphi > 360) {
0343         index = EBDetId(ieta, iphi - 360).rawId();
0344       } else {
0345         index = EBDetId(ieta, iphi).rawId();
0346       }
0347       fillmap[i++] = themap[index];
0348     }
0349   }
0350 
0351   uint32_t ebcenterid = getUnitWithMaxEnergy(themap);
0352 
0353   if (fillmap[i / 2] == themap[ebcenterid])
0354     return true;
0355   else
0356     return false;
0357 }
0358 
0359 float EcalBarrelSimHitsValidation::eCluster2x2(MapType &themap) {
0360   float E22 = 0.;
0361   float e012 = themap[0] + themap[1] + themap[2];
0362   float e036 = themap[0] + themap[3] + themap[6];
0363   float e678 = themap[6] + themap[7] + themap[8];
0364   float e258 = themap[2] + themap[5] + themap[8];
0365 
0366   if ((e012 > e678 || e012 == e678) && (e036 > e258 || e036 == e258))
0367     E22 = themap[0] + themap[1] + themap[3] + themap[4];
0368   else if ((e012 > e678 || e012 == e678) && (e036 < e258 || e036 == e258))
0369     E22 = themap[1] + themap[2] + themap[4] + themap[5];
0370   else if ((e012 < e678 || e012 == e678) && (e036 > e258 || e036 == e258))
0371     E22 = themap[3] + themap[4] + themap[6] + themap[7];
0372   else if ((e012 < e678 || e012 == e678) && (e036 < e258 || e036 == e258))
0373     E22 = themap[4] + themap[5] + themap[7] + themap[8];
0374 
0375   return E22;
0376 }
0377 
0378 float EcalBarrelSimHitsValidation::eCluster4x4(float e33, MapType &themap) {
0379   float E44 = 0.;
0380   float e0_4 = themap[0] + themap[1] + themap[2] + themap[3] + themap[4];
0381   float e0_20 = themap[0] + themap[5] + themap[10] + themap[15] + themap[20];
0382   float e4_24 = themap[4] + themap[9] + themap[14] + themap[19] + themap[24];
0383   float e0_24 = themap[20] + themap[21] + themap[22] + themap[23] + themap[24];
0384 
0385   if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
0386     E44 = e33 + themap[0] + themap[1] + themap[2] + themap[3] + themap[5] + themap[10] + themap[15];
0387   else if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
0388     E44 = e33 + themap[1] + themap[2] + themap[3] + themap[4] + themap[9] + themap[14] + themap[19];
0389   else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
0390     E44 = e33 + themap[5] + themap[10] + themap[15] + themap[20] + themap[21] + themap[22] + themap[23];
0391   else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
0392     E44 = e33 + themap[21] + themap[22] + themap[23] + themap[24] + themap[9] + themap[14] + themap[19];
0393 
0394   return E44;
0395 }
0396 
0397 uint32_t EcalBarrelSimHitsValidation::getUnitWithMaxEnergy(MapType &themap) {
0398   // look for max
0399   uint32_t unitWithMaxEnergy = 0;
0400   float maxEnergy = 0.;
0401 
0402   MapType::iterator iter;
0403   for (iter = themap.begin(); iter != themap.end(); iter++) {
0404     if (maxEnergy < (*iter).second) {
0405       maxEnergy = (*iter).second;
0406       unitWithMaxEnergy = (*iter).first;
0407     }
0408   }
0409 
0410   LogDebug("GeomInfo") << " max energy of " << maxEnergy << " GeV in Unit id " << unitWithMaxEnergy;
0411   return unitWithMaxEnergy;
0412 }
[ 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!