EcalDetId/src/EEDetId.cc

0001 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003
0004 #include <algorithm>
0005
0006 const int EEDetId::QuadColLimits[EEDetId::nCols + 1] = {0, 8, 17, 27, 36, 45, 54, 62, 70, 76, 79};
0007
0008 const int EEDetId::iYoffset[EEDetId::nCols + 1] = {0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0};
0009
0010 const unsigned short EEDetId::kxf[] = {
0011     41, 51, 41, 51, 41, 51, 36, 51, 36, 51, 26, 51, 26, 51, 26, 51, 21, 51, 21, 51, 21, 51, 21, 51, 21, 51, 16, 51, 16,
0012     51, 14, 51, 14, 51, 14, 51, 14, 51, 14, 51, 9,  51, 9,  51, 9,  51, 9,  51, 9,  51, 6,  51, 6,  51, 6,  51, 6,  51,
0013     6,  51, 6,  51, 6,  51, 6,  51, 6,  51, 6,  51, 4,  51, 4,  51, 4,  51, 4,  51, 4,  56, 1,  58, 1,  59, 1,  60, 1,
0014     61, 1,  61, 1,  62, 1,  62, 1,  62, 1,  62, 1,  62, 1,  62, 1,  62, 1,  62, 1,  62, 1,  62, 1,  61, 1,  61, 1,  60,
0015     1,  59, 1,  58, 4,  56, 4,  51, 4,  51, 4,  51, 4,  51, 6,  51, 6,  51, 6,  51, 6,  51, 6,  51, 6,  51, 6,  51, 6,
0016     51, 6,  51, 6,  51, 9,  51, 9,  51, 9,  51, 9,  51, 9,  51, 14, 51, 14, 51, 14, 51, 14, 51, 14, 51, 16, 51, 16, 51,
0017     21, 51, 21, 51, 21, 51, 21, 51, 21, 51, 26, 51, 26, 51, 26, 51, 36, 51, 36, 51, 41, 51, 41, 51, 41, 51};
0018
0019 const unsigned short EEDetId::kdi[] = {
0020     0,    10,   20,   30,   40,   50,   60,   75,   90,   105,  120,  145,  170,  195,  220,  245,  270,  300,  330,
0021     360,  390,  420,  450,  480,  510,  540,  570,  605,  640,  675,  710,  747,  784,  821,  858,  895,  932,  969,
0022     1006, 1043, 1080, 1122, 1164, 1206, 1248, 1290, 1332, 1374, 1416, 1458, 1500, 1545, 1590, 1635, 1680, 1725, 1770,
0023     1815, 1860, 1905, 1950, 1995, 2040, 2085, 2130, 2175, 2220, 2265, 2310, 2355, 2400, 2447, 2494, 2541, 2588, 2635,
0024     2682, 2729, 2776, 2818, 2860, 2903, 2946, 2988, 3030, 3071, 3112, 3152, 3192, 3232, 3272, 3311, 3350, 3389, 3428,
0025     3467, 3506, 3545, 3584, 3623, 3662, 3701, 3740, 3779, 3818, 3857, 3896, 3935, 3974, 4013, 4052, 4092, 4132, 4172,
0026     4212, 4253, 4294, 4336, 4378, 4421, 4464, 4506, 4548, 4595, 4642, 4689, 4736, 4783, 4830, 4877, 4924, 4969, 5014,
0027     5059, 5104, 5149, 5194, 5239, 5284, 5329, 5374, 5419, 5464, 5509, 5554, 5599, 5644, 5689, 5734, 5779, 5824, 5866,
0028     5908, 5950, 5992, 6034, 6076, 6118, 6160, 6202, 6244, 6281, 6318, 6355, 6392, 6429, 6466, 6503, 6540, 6577, 6614,
0029     6649, 6684, 6719, 6754, 6784, 6814, 6844, 6874, 6904, 6934, 6964, 6994, 7024, 7054, 7079, 7104, 7129, 7154, 7179,
0030     7204, 7219, 7234, 7249, 7264, 7274, 7284, 7294, 7304, 7314};
0031
0032 EEDetId::EEDetId(int index1, int index2, int iz, int mode) : DetId(Ecal, EcalEndcap) {
0033   int crystal_ix = 0;
0034   int crystal_iy = 0;
0035   if (mode == XYMODE) {
0036     crystal_ix = index1;
0037     crystal_iy = index2;
0038   } else if (mode == SCCRYSTALMODE) {
0039     int SC = index1;
0040     int crystal = index2;
0041     //      std::cout << "iz " << iz << " SC " << index1 << "crystal " << index2  << std::endl;
0042
0043     crystal_ix = iz * ix(SC, crystal);
0044     if (crystal_ix < 0)
0045       crystal_ix++;
0046     crystal_ix += 50;
0047     crystal_iy = iy(SC, crystal);
0048     if (crystal_iy < 0)
0049       crystal_iy++;
0050     crystal_iy += 50;
0051
0052   } else {
0053     throw cms::Exception("InvalidDetId")
0054         << "EEDetId:  Cannot create object.  Unknown mode for (int, int, int) constructor.";
0055   }
0056
0057   if (!validDetId(crystal_ix, crystal_iy, iz)) {
0058     throw cms::Exception("InvalidDetId") << "EEDetId:  Cannot create object.  Indexes out of bounds \n"
0059                                          << "x = " << crystal_ix << " y = " << crystal_iy << " z = " << iz;
0060   }
0061
0062   id_ |= (crystal_iy & 0x7f) | ((crystal_ix & 0x7f) << 7) | ((iz > 0) ? (0x4000) : (0));
0063 }
0064
0065 EEDetId EEDetId::unhashIndex(int hi) {
0066   if (validHashIndex(hi)) {
0067     const int iz(hi < kEEhalf ? -1 : 1);
0068     const uint32_t di(hi % kEEhalf);
0069     const int ii((std::upper_bound(kdi, kdi + (2 * IY_MAX), di) - kdi) - 1);
0070     const int iy(1 + ii / 2);
0071     const int ix(kxf[ii] + di - kdi[ii]);
0072     return EEDetId(ix, iy, iz);
0073   } else {
0074     return EEDetId();
0075   }
0076 }
0077
0078 int EEDetId::ix(int iSC, int iCrys) const {
0079   /*
0080    *  ix() return individual crystal x-coordinate
0081    *
0082    *  Author    : B W Kennedy
0083    *  Version   : 1.00
0084    *  Created   : 21 December 2005
0085    *  Last Mod  : 31 January 2006
0086    *
0087    *  Input     : iSC, iCrys - Supercrystal and crystal ids
0088    */
0089
0090   int nSCinQuadrant = QuadColLimits[nCols];
0091
0092   if (iSC > 4 * nSCinQuadrant || iSC < 1) {
0093     throw new std::exception();
0094   }
0095
0096   //  Map SC number into (x>0,y>0) quadrant.
0097   int iSCmap, iqx, iq;
0098   if (iSC > 3 * nSCinQuadrant) {
0099     iSCmap = iSC - 3 * nSCinQuadrant;
0100     iqx = 1;
0101     iq = 4;
0102   } else if (iSC > 2 * nSCinQuadrant) {
0103     iSCmap = iSC - 2 * nSCinQuadrant;
0104     iqx = -1;
0105     iq = 3;
0106   } else if (iSC > nSCinQuadrant) {
0107     iSCmap = iSC - nSCinQuadrant;
0108     iqx = -1;
0109     iq = 2;
0110   } else {
0111     iSCmap = iSC;
0112     iqx = 1;
0113     iq = 1;
0114   }
0115
0116   // Decide which column the SC is in
0117   int iCol = 0;
0118   while (iSCmap > QuadColLimits[iCol++])
0119     ;
0120   iCol--;
0121
0122   int ixCrys = -1;
0123   if (iq == 1 || iq == 3)
0124     ixCrys = iqx * (5 * (iCol - 1) + (int)(iCrys + 4) / 5);
0125   else if (iq == 2 || iq == 4)
0126     ixCrys = iqx * (5 * (iCol - 1) + (iCrys - 1) % 5 + 1);
0127
0128   // returning a value from 1 to 100
0129
0130   return ixCrys;
0131 }
0132
0133 int EEDetId::iy(int iSC, int iCrys) const {
0134   /*
0135    *  iy() return individual crystal y-coordinate
0136    *
0137    *  Author    : B W Kennedy
0138    *  Version   : 1.00
0139    *  Created   : 21 December 2005
0140    *  Last Mod  : 31 January 2006
0141    *
0142    *  Input     : iSC, iCrys - Supercrystal and crystal ids
0143    */
0144
0145   int nSCinQuadrant = QuadColLimits[nCols];
0146   if (iSC > 4 * nSCinQuadrant || iSC < 1) {
0147     throw new std::exception();
0148   }
0149
0150   //  Map SC number into (x>0,y>0) quadrant
0151   int iSCmap, iqy, iq;
0152   if (iSC > 3 * nSCinQuadrant) {
0153     iSCmap = iSC - 3 * nSCinQuadrant;
0154     iqy = -1;
0155     iq = 4;
0156   } else if (iSC > 2 * nSCinQuadrant) {
0157     iSCmap = iSC - 2 * nSCinQuadrant;
0158     iqy = -1;
0159     iq = 3;
0160   } else if (iSC > nSCinQuadrant) {
0161     iSCmap = iSC - nSCinQuadrant;
0162     iqy = 1;
0163     iq = 2;
0164   } else {
0165     iSCmap = iSC;
0166     iqy = 1;
0167     iq = 1;
0168   }
0169
0170   // Decide which column the SC is in
0171   int iCol = 0;
0172   while (iSCmap > QuadColLimits[iCol++])
0173     ;
0174   iCol--;
0175
0176   int iSCy = iSCmap - QuadColLimits[iCol - 1] + iYoffset[iCol];
0177
0178   int iyCrys = -1;
0179   if (iq == 1 || iq == 3)
0180     iyCrys = iqy * (5 * (iSCy - 1) + (iCrys - 1) % 5 + 1);
0181   else if (iq == 2 || iq == 4)
0182     iyCrys = iqy * (5 * (iSCy - 1) + (int)(iCrys + 4) / 5);
0183   return iyCrys;
0184 }
0185
0186 int EEDetId::ixQuadrantOne() const {
0187   int iQuadrant = iquadrant();
0188   if (iQuadrant == 1 || iQuadrant == 4)
0189     return (ix() - 50);
0190   else if (iQuadrant == 2 || iQuadrant == 3)
0191     return (51 - ix());
0192   //Should never be reached
0193   return -1;
0194 }
0195
0196 int EEDetId::iyQuadrantOne() const {
0197   int iQuadrant = iquadrant();
0198   if (iQuadrant == 1 || iQuadrant == 2)
0199     return (iy() - 50);
0200   else if (iQuadrant == 3 || iQuadrant == 4)
0201     return 51 - iy();
0202   //Should never be reached
0203   return -1;
0204 }
0205
0206 int EEDetId::iquadrant() const {
0207   if (ix() > 50) {
0208     if (iy() > 50)
0209       return 1;
0210     else
0211       return 4;
0212   } else {
0213     if (iy() > 50)
0214       return 2;
0215     else
0216       return 3;
0217   }
0218   //Should never be reached
0219   return -1;
0220 }
0221
0222 int EEDetId::isc() const { return isc(1 + (ix() - 1) / nCrys, 1 + (iy() - 1) / nCrys); }
0223
0224 int EEDetId::isc(int jx, int jy) {
0225   if (0 < jx && 21 > jx && 0 < jy && 21 > jy) {
0226     const int iquad((10 < jx && 10 < jy ? 1 : (11 > jx && 10 < jy ? 2 : (11 > jx && 11 > jy ? 3 : 4))));
0227
0228     const int iCol = (1 == iquad || 4 == iquad ? jx - 10 : 11 - jx);
0229     const int iRow = (1 == iquad || 2 == iquad ? jy - 10 : 11 - jy);
0230
0231     static const int nSCinQuadrant = ISC_MAX / 4;
0232
0233     const int yOff(iYoffset[iCol]);
0234
0235     const int qOff(nSCinQuadrant * (iquad - 1));
0236
0237     const int iscOne(QuadColLimits[iCol - 1] + iRow - yOff);
0238
0239     return (yOff >= iRow ? -1 : (QuadColLimits[iCol] < iscOne ? -2 : iscOne + qOff));
0240   } else {
0241     return -3;  // bad inputs
0242   }
0243 }
0244
0245 int EEDetId::ic() const {
0246   /*
0247     *  Return crystal number from (x,y) coordinates.
0248     *
0249     *  Author    : B W Kennedy
0250     *  Version   : 1.00
0251     *  Created   : 5 May 2006
0252     *  Last Mod  :
0253     *
0254     *  Input     : ix, iy - (x,y) position of crystal
0255     */
0256
0257   /*  Useful constants . */
0258   int iQuadrant = iquadrant();
0259   int icrCol = -1;
0260   int icrRow = -1;
0261
0262   if (iQuadrant == 1 || iQuadrant == 3) {
0263     icrCol = (ixQuadrantOne() - 1) % nCrys;
0264     icrRow = (iyQuadrantOne() - 1) % nCrys;
0265   }
0266
0267   else if (iQuadrant == 2 || iQuadrant == 4) {
0268     icrRow = (ixQuadrantOne() - 1) % nCrys;
0269     icrCol = (iyQuadrantOne() - 1) % nCrys;
0270   }
0271
0272   int icrys = 5 * icrCol + icrRow + 1;
0273
0274   return icrys;
0275 }
0276
0277 bool EEDetId::isNextToBoundary(EEDetId id) { return isNextToDBoundary(id) || isNextToRingBoundary(id); }
0278
0279 bool EEDetId::isNextToDBoundary(EEDetId id) {
0280   // hardcoded values for D boundary
0281   return id.ix() == 50 || id.ix() == 51;
0282 }
0283
0284 bool EEDetId::isNextToRingBoundary(EEDetId id) {
0285   for (int i = -1; i <= 1; ++i) {
0286     for (int j = -1; j <= 1; ++j) {
0287       if (!validDetId(id.ix() + i, id.iy() + j, id.zside())) {
0288         return true;
0289       }
0290     }
0291   }
0292   return false;
0293 }
0294
0295 int EEDetId::iPhiOuterRing() const {
0296   int returnValue(0);
0297   if (isOuterRing()) {
0298     const int ax(abs(ix() - IX_MAX / 2));
0299     const int ay(abs(iy() - IY_MAX / 2));
0300     returnValue = ax + 50 - ay;
0301     if (ay <= 47)
0302       --returnValue;
0303     if (ay <= 45)
0304       --returnValue;
0305     if (ay <= 42)
0306       --returnValue;
0307     if (ay <= 37)
0308       --returnValue;
0309     if (ay <= 35)
0310       --returnValue;
0311     if (ay <= 30)
0312       --returnValue;
0313     if (ay <= 25)
0314       --returnValue;
0315     if (ay <= 15)
0316       --returnValue;
0317     if (ay <= 10)
0318       --returnValue;
0319     const int iq(iquadrant());
0320     if (1 == iq) {
0321       returnValue = 91 - returnValue;
0322     } else {
0323       if (2 == iq) {
0324         returnValue += 90;
0325       } else {
0326         if (3 == iq) {
0327           returnValue = 271 - returnValue;
0328         } else {
0329           returnValue += 270;
0330         }
0331       }
0332     }
0333     returnValue = 1 + (360 + returnValue - 10 - 1) % 360;
0334   }
0335   //   if( positiveZ() ) returnValue += 360 ;
0336   return returnValue;
0337 }
0338
0339 EEDetId EEDetId::idOuterRing(int iPhi, int zEnd) {
0340   iPhi -= 10;  // phi=1 in barrel is at -10deg
0341   while (iPhi < 1)
0342     iPhi += 360;
0343   while (iPhi > 360)
0344     iPhi -= 360;
0345
0346   const int index1(iPhi - 1);
0347   const int quad(index1 / 90);
0348   int indexq(index1 - quad * 90 + 1);
0349   if (0 == quad || 2 == quad)
0350     indexq = 91 - indexq;
0351   const int indexh(indexq > 45 ? 91 - indexq : indexq);
0352   const int axh(
0353       indexh <= 10
0354           ? indexh
0355           : (indexh <= 12
0356                  ? 10
0357                  : (indexh <= 17
0358                         ? indexh - 2
0359                         : (indexh <= 18
0360                                ? 15
0361                                : (indexh <= 28
0362                                       ? indexh - 3
0363                                       : (indexh <= 30
0364                                              ? 25
0365                                              : (indexh <= 35
0366                                                     ? indexh - 5
0367                                                     : (indexh <= 39 ? 30 : (indexh <= 44 ? indexh - 9 : 35)))))))));
0368   const int ayh(
0369       indexh <= 10
0370           ? 50
0371           : (indexh <= 12
0372                  ? 60 - indexh
0373                  : (indexh <= 17
0374                         ? 47
0375                         : (indexh <= 18
0376                                ? 64 - indexh
0377                                : (indexh <= 28
0378                                       ? 45
0379                                       : (indexh <= 30 ? 73 - indexh
0380                                                       : (indexh <= 35 ? 42
0381                                                                       : (indexh <= 39 ? 77 - indexh
0382                                                                                       : (indexh <= 44 ? 37 : 36)))))))));
0383   const int bxh(indexq > 45 ? ayh : axh);
0384   const int byh(indexq > 45 ? axh : ayh);
0385   const int cx((quad == 0 || quad == 3 ? bxh : -bxh + 1) + IX_MAX / 2);
0386   const int cy((quad == 0 || quad == 1 ? byh : -byh + 1) + IY_MAX / 2);
0387
0388   return EEDetId(cx, cy, (zEnd > 0 ? 1 : -1));
0389 }
0390
0391 EEDetId EEDetId::offsetBy(int nrStepsX, int nrStepsY) const {
0392   int newX = ix() + nrStepsX;
0393   int newY = iy() + nrStepsY;
0394
0395   if (validDetId(newX, newY, zside())) {
0396     return EEDetId(newX, newY, zside());
0397   } else {
0398     return EEDetId(0);
0399   }
0400 }
0401
0402 EEDetId EEDetId::switchZSide() const {
0403   int newZSide = -1 * zside();
0404   if (validDetId(ix(), iy(), newZSide)) {
0405     return EEDetId(ix(), iy(), newZSide);
0406   } else {
0407     return EEDetId(0);
0408   }
0409 }
0410
0411 DetId EEDetId::offsetBy(const DetId startId, int nrStepsX, int nrStepsY) {
0412   if (startId.det() == DetId::Ecal && startId.subdetId() == EcalEndcap) {
0413     EEDetId eeStartId(startId);
0414     return eeStartId.offsetBy(nrStepsX, nrStepsY).rawId();
0415   } else {
0416     return DetId(0);
0417   }
0418 }
0419
0420 DetId EEDetId::switchZSide(const DetId startId) {
0421   if (startId.det() == DetId::Ecal && startId.subdetId() == EcalEndcap) {
0422     EEDetId eeStartId(startId);
0423     return eeStartId.switchZSide().rawId();
0424   } else {
0425     return DetId(0);
0426   }
0427 }
0428
0429 bool EEDetId::isOuterRing() const {
0430   const int kx(ix());
0431   const int ky(iy());
0432   const int ax(kx > IX_MAX / 2 ? kx - IX_MAX / 2 : IX_MAX / 2 + 1 - kx);
0433   const int ay(ky > IY_MAX / 2 ? ky - IY_MAX / 2 : IY_MAX / 2 + 1 - ky);
0434   return (isOuterRingXY(ax, ay) || isOuterRingXY(ay, ax));
0435 }
0436
0437 bool EEDetId::isOuterRingXY(int ax, int ay) {
0438   return ((ax <= 10 && ay == 50) || (ax == 10 && ay >= 48) || (ax <= 15 && ax >= 11 && ay == 47) ||
0439           (ax == 15 && ay == 46) || (ax <= 25 && ax >= 16 && ay == 45) || (ax == 25 && ay <= 44 && ay >= 43) ||
0440           (ax <= 30 && ax >= 26 && ay == 42) || (ax == 30 && ay <= 41 && ay >= 38) ||
0441           (ax <= 35 && ax >= 31 && ay == 37) || (ax == 35 && ay == 36));
0442 }
0443
0444 bool EEDetId::slowValidDetId(int crystal_ix, int crystal_iy) {
0445   return  // negative logic!
0446       !((crystal_ix >= 1 && crystal_ix <= 3 && (crystal_iy <= 40 || crystal_iy > 60)) ||
0447         (crystal_ix >= 4 && crystal_ix <= 5 && (crystal_iy <= 35 || crystal_iy > 65)) ||
0448         (crystal_ix >= 6 && crystal_ix <= 8 && (crystal_iy <= 25 || crystal_iy > 75)) ||
0449         (crystal_ix >= 9 && crystal_ix <= 13 && (crystal_iy <= 20 || crystal_iy > 80)) ||
0450         (crystal_ix >= 14 && crystal_ix <= 15 && (crystal_iy <= 15 || crystal_iy > 85)) ||
0451         (crystal_ix >= 16 && crystal_ix <= 20 && (crystal_iy <= 13 || crystal_iy > 87)) ||
0452         (crystal_ix >= 21 && crystal_ix <= 25 && (crystal_iy <= 8 || crystal_iy > 92)) ||
0453         (crystal_ix >= 26 && crystal_ix <= 35 && (crystal_iy <= 5 || crystal_iy > 95)) ||
0454         (crystal_ix >= 36 && crystal_ix <= 39 && (crystal_iy <= 3 || crystal_iy > 97)) ||
0455         (crystal_ix >= 98 && crystal_ix <= 100 && (crystal_iy <= 40 || crystal_iy > 60)) ||
0456         (crystal_ix >= 96 && crystal_ix <= 97 && (crystal_iy <= 35 || crystal_iy > 65)) ||
0457         (crystal_ix >= 93 && crystal_ix <= 95 && (crystal_iy <= 25 || crystal_iy > 75)) ||
0458         (crystal_ix >= 88 && crystal_ix <= 92 && (crystal_iy <= 20 || crystal_iy > 80)) ||
0459         (crystal_ix >= 86 && crystal_ix <= 87 && (crystal_iy <= 15 || crystal_iy > 85)) ||
0460         (crystal_ix >= 81 && crystal_ix <= 85 && (crystal_iy <= 13 || crystal_iy > 87)) ||
0461         (crystal_ix >= 76 && crystal_ix <= 80 && (crystal_iy <= 8 || crystal_iy > 92)) ||
0462         (crystal_ix >= 66 && crystal_ix <= 75 && (crystal_iy <= 5 || crystal_iy > 95)) ||
0463         (crystal_ix >= 62 && crystal_ix <= 65 && (crystal_iy <= 3 || crystal_iy > 97)) ||
0464         ((crystal_ix == 40 || crystal_ix == 61) &&
0465          ((crystal_iy >= 46 && crystal_iy <= 55) || crystal_iy <= 3 || crystal_iy > 97)) ||
0466         ((crystal_ix == 41 || crystal_ix == 60) && crystal_iy >= 44 && crystal_iy <= 57) ||
0467         ((crystal_ix == 42 || crystal_ix == 59) && crystal_iy >= 43 && crystal_iy <= 58) ||
0468         ((crystal_ix == 43 || crystal_ix == 58) && crystal_iy >= 42 && crystal_iy <= 59) ||
0469         ((crystal_ix == 44 || crystal_ix == 45 || crystal_ix == 57 || crystal_ix == 56) && crystal_iy >= 41 &&
0470          crystal_iy <= 60) ||
0471         (crystal_ix >= 46 && crystal_ix <= 55 && crystal_iy >= 40 && crystal_iy <= 61));
0472 }
0473
0474 int EEDetId::distanceX(const EEDetId& a, const EEDetId& b) { return abs(a.ix() - b.ix()); }
0475
0476 int EEDetId::distanceY(const EEDetId& a, const EEDetId& b) { return abs(a.iy() - b.iy()); }
0477
0478 #include <ostream>
0479 std::ostream& operator<<(std::ostream& s, const EEDetId& id) {
0480   return s << "(EE iz " << ((id.zside() > 0) ? ("+ ") : ("- ")) << " ix " << id.ix() << " , iy " << id.iy() << ')';
0481 }