Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:53:15

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 }