Tools/src/EcalRingCalibrationTools.cc

0001 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
0002 #include "DataFormats/DetId/interface/DetId.h"
0003
0004 //Includes need to read from geometry
0005 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <iostream>
0012
0013 //by default is not initialized, gets initialized at first call
0014 std::atomic<bool> EcalRingCalibrationTools::isInitializedFromGeometry_(false);
0015 short EcalRingCalibrationTools::endcapRingIndex_[EEDetId::IX_MAX][EEDetId::IY_MAX];
0016 std::once_flag EcalRingCalibrationTools::once_;
0017
0018 constexpr short EcalRingCalibrationTools::N_RING_TOTAL;
0019 constexpr short EcalRingCalibrationTools::N_RING_BARREL;
0020 constexpr short EcalRingCalibrationTools::N_RING_ENDCAP;
0021
0022 constexpr short EcalRingCalibrationTools::N_MODULES_BARREL;
0023
0024 short EcalRingCalibrationTools::getRingIndex(DetId id) {
0025   if (id.det() != DetId::Ecal)
0026     return -1;
0027
0028   if (id.subdetId() == EcalBarrel) {
0029     if (EBDetId(id).ieta() < 0)
0030       return EBDetId(id).ieta() + 85;
0031     else
0032       return EBDetId(id).ieta() + 84;
0033   }
0034   if (id.subdetId() == EcalEndcap) {
0035     //needed only for the EE, it can be replaced at some point with something smarter
0036     if (not isInitializedFromGeometry_)
0037       throw std::logic_error(
0038           "EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
0039
0040     EEDetId eid(id);
0041     short endcapRingIndex = endcapRingIndex_[eid.ix() - 1][eid.iy() - 1] + N_RING_BARREL;
0042     if (eid.zside() == 1)
0043       endcapRingIndex += N_RING_ENDCAP / 2;
0044     return endcapRingIndex;
0045   }
0046   return -1;
0047 }
0048
0049 short EcalRingCalibrationTools::getModuleIndex(DetId id) {
0050   if (id.det() != DetId::Ecal)
0051     return -1;
0052
0053   if (id.subdetId() == EcalBarrel) {
0054     short module = 4 * (EBDetId(id).ism() - 1) + EBDetId(id).im() - 1;
0055
0056     //      std::cout<<"SM construction # : "<<EBDetId(id).ism() <<" SM installation # : "<<  installationSMNumber[ EBDetId(id).ism() -1 ]<< "Xtal is Module :"<<module<< std::endl;
0057     return module;
0058   }
0059   if (id.subdetId() == EcalEndcap) {
0060     return -1;
0061   }
0062
0063   return -1;
0064 }
0065
0066 std::vector<DetId> EcalRingCalibrationTools::getDetIdsInRing(short etaIndex) {
0067   std::vector<DetId> ringIds;
0068   if (etaIndex < 0)
0069     return ringIds;
0070
0071   if (etaIndex < N_RING_BARREL) {
0072     int k = 0;
0073     if (etaIndex < 85)
0074       k = -85 + etaIndex;
0075     else
0076       k = etaIndex - 84;
0077
0078     for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
0079
0080       if (EBDetId::validDetId(k, iphi))
0081         ringIds.push_back(EBDetId(k, iphi));
0082   }
0083
0084   else if (etaIndex < N_RING_TOTAL) {
0085     //needed only for the EE, it can be replaced at some point maybe with something smarter
0086     if (not isInitializedFromGeometry_)
0087       throw std::logic_error(
0088           "EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
0089
0090     int zside = (etaIndex < N_RING_BARREL + (N_RING_ENDCAP / 2)) ? -1 : 1;
0091     short eeEtaIndex = (etaIndex - N_RING_BARREL) % (N_RING_ENDCAP / 2);
0092
0093     for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
0094       for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
0095         if (endcapRingIndex_[ix][iy] == eeEtaIndex)
0096           ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
0097   }
0098
0099   return ringIds;
0100 }
0101
0102 std::vector<DetId> EcalRingCalibrationTools::getDetIdsInECAL() {
0103   std::vector<DetId> ringIds;
0104
0105   for (int ieta = -EBDetId::MAX_IETA; ieta <= EBDetId::MAX_IETA; ++ieta)
0106     for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
0107       if (EBDetId::validDetId(ieta, iphi))
0108         ringIds.push_back(EBDetId(ieta, iphi));
0109
0110   //needed only for the EE, it can be replaced at some point maybe with something smarter
0111   if (not isInitializedFromGeometry_)
0112     throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
0113
0114   for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
0115     for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
0116       for (int zside = -1; zside < 2; zside += 2)
0117         if (EEDetId::validDetId(ix + 1, iy + 1, zside))
0118           ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
0119
0120   //  std::cout<<" [EcalRingCalibrationTools::getDetIdsInECAL()] DetId.size() is "<<ringIds.size()<<std::endl;
0121
0122   return ringIds;
0123 }
0124
0125 std::vector<DetId> EcalRingCalibrationTools::getDetIdsInModule(short moduleIndex) {
0126   std::vector<DetId> ringIds;
0127   if (moduleIndex < 0)
0128     return ringIds;
0129
0130   short moduleBound[5] = {1, 26, 46, 66, 86};
0131   if (moduleIndex < EcalRingCalibrationTools::N_MODULES_BARREL) {
0132     ////////////////////////////nuovo ciclo
0133     short sm, moduleInSm, zsm;
0134
0135     short minModuleiphi, maxModuleiphi, minModuleieta = 360, maxModuleieta = 0;
0136
0137     //  if(moduleIndex%4 != 0 )
0138     sm = moduleIndex / 4 + 1;
0139     //  else
0140     //sm = moduleIndex/4;//i.e. module 8 belongs to sm=3, not sm=3
0141
0142     //if(moduleIndex%4 != 0 )
0143     moduleInSm = moduleIndex % 4;
0144     //else
0145     //moduleInSm = 4;//moduleInSm is [1,2,3,4]
0146
0147     if (moduleIndex > 71)
0148       zsm = -1;
0149     else
0150       zsm = 1;
0151
0152     minModuleiphi = ((sm - 1) % 18 + 1) * 20 - 19;
0153     maxModuleiphi = ((sm - 1) % 18 + 1) * 20;
0154
0155     if (zsm == 1) {
0156       minModuleieta = moduleBound[moduleInSm];
0157       maxModuleieta = moduleBound[moduleInSm + 1] - 1;
0158     } else if (zsm == -1) {
0159       minModuleieta = -moduleBound[moduleInSm + 1] + 1;
0160       maxModuleieta = -moduleBound[moduleInSm];
0161     }
0162     ////////////////////////////nuovo ciclo
0163
0164     std::cout << "Called moduleIndex " << moduleIndex << std::endl;
0165     std::cout << "minModuleieta " << minModuleieta << " maxModuleieta " << maxModuleieta << " minModuleiphi "
0166               << minModuleiphi << " maxModuleiphi " << maxModuleiphi << std::endl;
0167
0168     for (int ieta = minModuleieta; ieta <= maxModuleieta; ++ieta) {
0169       for (int iphi = minModuleiphi; iphi <= maxModuleiphi; ++iphi) {
0170         ringIds.push_back(EBDetId(ieta, iphi));
0171
0172         //        std::cout<<"Putting Xtal with ieta: "<<ieta<<" iphi "<<iphi<<" of SM "<<sm<<" into Module "<<moduleIndex<<std::endl;
0173
0174       }  //close loop on phi
0175     }    //close loop on eta
0176   }      //close if ( moduleInstallationNumber < 144)
0177
0178   return ringIds;
0179 }
0180
0181 void EcalRingCalibrationTools::setCaloGeometry(const CaloGeometry* geometry) {
0182   std::call_once(once_, EcalRingCalibrationTools::initializeFromGeometry, geometry);
0183 }
0184
0185 void EcalRingCalibrationTools::initializeFromGeometry(CaloGeometry const* geometry) {
0186   if (not geometry)
0187     throw std::invalid_argument("EcalRingCalibrationTools::initializeFromGeometry called with a nullptr argument");
0188
0189   float cellPosEta[EEDetId::IX_MAX][EEDetId::IY_MAX];
0190   for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
0191     for (int iy = 0; iy < EEDetId::IY_MAX; ++iy) {
0192       cellPosEta[ix][iy] = -1.;
0193       endcapRingIndex_[ix][iy] = -9;
0194     }
0195
0196   CaloSubdetectorGeometry const* endcapGeometry = geometry->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0197   if (not endcapGeometry)
0198     throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry not found");
0199
0200   std::vector<DetId> const& endcapCells = geometry->getValidDetIds(DetId::Ecal, EcalEndcap);
0201
0202   for (std::vector<DetId>::const_iterator endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); ++endcapIt) {
0203     EEDetId ee(*endcapIt);
0204     if (ee.zside() == -1)
0205       continue;  //Just using +side to fill absEta x,y map
0206     auto cellGeometry = endcapGeometry->getGeometry(*endcapIt);
0207     int ics = ee.ix() - 1;
0208     int ips = ee.iy() - 1;
0209     cellPosEta[ics][ips] = fabs(cellGeometry->getPosition().eta());
0210     //std::cout<<"EE Xtal, |eta| is "<<fabs(cellGeometry->getPosition().eta())<<std::endl;
0211   }
0212
0213   float eta_ring[N_RING_ENDCAP / 2];
0214   for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring)
0215     eta_ring[ring] = cellPosEta[ring][50];
0216
0217   double etaBoundary[N_RING_ENDCAP / 2 + 1];
0218   etaBoundary[0] = 1.47;
0219   etaBoundary[N_RING_ENDCAP / 2] = 4.0;
0220
0221   for (int ring = 1; ring < N_RING_ENDCAP / 2; ++ring)
0222     etaBoundary[ring] = (eta_ring[ring] + eta_ring[ring - 1]) / 2.;
0223
0224   for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring) {
0225     // std::cout<<"***********************EE ring: "<<ring<<" eta "<<(etaBoundary[ring] + etaBoundary[ring+1])/2.<<std::endl;
0226     for (int ix = 0; ix < EEDetId::IX_MAX; ix++)
0227       for (int iy = 0; iy < EEDetId::IY_MAX; iy++)
0228         if (cellPosEta[ix][iy] > etaBoundary[ring] && cellPosEta[ix][iy] < etaBoundary[ring + 1]) {
0229           endcapRingIndex_[ix][iy] = ring;
0230           //std::cout<<"endcapRing_["<<ix+1<<"]["<<iy+1<<"] = "<<ring<<";"<<std::endl;
0231         }
0232   }
0233
0234   std::vector<DetId> const& barrelCells = geometry->getValidDetIds(DetId::Ecal, EcalBarrel);
0235
0236   for (std::vector<DetId>::const_iterator barrelIt = barrelCells.begin(); barrelIt != barrelCells.end(); ++barrelIt) {
0237     EBDetId eb(*barrelIt);
0238   }
0239
0240   //EB
0241
0242   isInitializedFromGeometry_ = true;
0243 }