RPCTrigger/src/RPCStripsRing.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     RPCConeBuilder
0004 // Class  :     RPCStripsRing
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author:  Tomasz Fruboes
0010 //         Created:  Tue Feb 26 15:13:10 CET 2008
0011 //
0012
0013 // system include files
0014
0015 // user include files
0016 //#include "L1TriggerConfig/RPCConeBuilder/interface/RPCStripsRing.h"
0017 #include "L1Trigger/RPCTrigger/interface/RPCStripsRing.h"
0018 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
0019 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0020 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022
0023 RPCStripsRing::RPCStripsRing()
0024     : m_hwPlane(-1),
0025       m_etaPartition(99),
0026       m_region(-2),
0027       m_isReferenceRing(false),
0028       m_didVirtuals(false),
0029       m_didFiltering(false) {}
0030
0031 RPCStripsRing::RPCStripsRing(const RPCRoll* roll, std::shared_ptr<L1RPCConeBuilder::TConMap> cmap)
0032     : m_didVirtuals(false), m_didFiltering(false), m_connectionsMap(cmap) {
0033   RPCDetId detId = roll->id();
0034   RPCGeomServ grs(detId);
0035
0036   m_etaPartition = grs.eta_partition();
0037   m_hwPlane = calculateHwPlane(roll);
0038
0039   m_isReferenceRing = false;
0040
0041   m_region = detId.region();
0042
0043   int ring = detId.ring();
0044
0045   if (m_region == 0 && std::abs(ring) < 2 && m_hwPlane == 2)  // for barell wheel -1,0,1 refplane is hwPlane=2
0046     m_isReferenceRing = true;
0047   else if (m_region == 0 && std::abs(ring) == 2 && m_hwPlane == 6)  // for barell wheel -2,2 refplane is hwPlane=6
0048     m_isReferenceRing = true;
0049   else if (m_region != 0 && m_hwPlane == 2)  // for endcaps
0050     m_isReferenceRing = true;
0051
0052   if (getRingId() == 2008 || getRingId() == 2108)  //exception: endcaps;hwplane 2;farest roll from beam
0053     m_isReferenceRing = false;
0054
0055   addRoll(roll);
0056 }
0057
0058 void RPCStripsRing::addRoll(const RPCRoll* roll) {
0059   //  RPCDetId detId = roll->id();
0060
0061   if (getRingId() != getRingId(roll)) {
0062     throw cms::Exception("RPCInternal") << "RPCStripsRing::addRoll ringsIds dont match \n";
0063   }
0064
0065   //iterate over the strips of this roll
0066   for (int i = 1; i <= roll->nstrips(); i++) {
0067     LocalPoint lStripCentre = roll->centreOfStrip(i);
0068     GlobalPoint gStripCentre = roll->toGlobal(lStripCentre);
0069     float phiRaw = gStripCentre.phi();
0070
0071     TStrip newStrip(roll->id().rawId(), i);
0072     (*this)[phiRaw] = newStrip;
0073   }
0074 }
0075
0076 int RPCStripsRing::getRingId(int etaPart, int hwPlane) {
0077   int sign = 1;  // positive
0078   if (etaPart < 0) {
0079     sign = 0;
0080   }
0081
0082   return 1000 * (hwPlane) +        //1...6
0083          100 * (sign) +            //
0084          1 * (std::abs(etaPart));  //-17...17
0085 }
0086
0087 int RPCStripsRing::getRingId() { return getRingId(m_etaPartition, m_hwPlane); }
0088
0089 int RPCStripsRing::getRingId(const RPCRoll* roll) {
0090   RPCDetId detId = roll->id();
0091   RPCGeomServ grs(detId);
0092   int etaPartition = grs.eta_partition();
0093   int hwPlane = calculateHwPlane(roll);
0094
0095   return getRingId(etaPartition, hwPlane);
0096 }
0097
0098 //  hwPlane is  station number for endcaps
0099 //  for barrell numbering goes 1 5 2 6 3 4 (first number means plane closest to the beam)
0100 int RPCStripsRing::calculateHwPlane(const RPCRoll* roll) {
0101   int hwPlane = -1;
0102   RPCDetId detId = roll->id();
0103   int station = detId.station();
0104   int layer = detId.layer();
0105   int region = detId.region();
0106
0107   if (region != 0) {  // endcaps
0108     hwPlane = station;
0109   }
0110   // Now comes the barell
0111   else if (station > 2) {
0112     hwPlane = station;
0113   } else if (station == 1 && layer == 1) {
0114     hwPlane = 1;
0115   } else if (station == 1 && layer == 2) {
0116     hwPlane = 5;
0117   } else if (station == 2 && layer == 1) {
0118     hwPlane = 2;
0119   } else if (station == 2 && layer == 2) {
0120     hwPlane = 6;
0121   }
0122
0123   /*if (hwPlane < 1)
0124     std::cout << "prb: " << hwPlane << " "
0125         << region << " "
0126         << station << " "
0127   << layer << std::endl;*/
0128   if (hwPlane < 0) {
0129     throw cms::Exception("RPCInternal") << "Calculated negative hwplane \n";
0130   }
0131
0132   return hwPlane;
0133 }
0134
0135 void RPCStripsRing::filterOverlapingChambers() {
0136   if (m_didFiltering)
0137     return;
0138   m_didFiltering = true;
0139
0140   if (m_region != 0 || m_hwPlane != 4)
0141     return;
0142
0143   typedef std::map<uint32_t, int> TDetId2StripNo;
0144   TDetId2StripNo det2stripNo;
0145
0146   // Note: we begin in middle of first chamber (ch1), we have to handle that
0147   int ch1EndStrips = 0;  // no of strips on the end of the map (first=last chamber of map)
0148
0149   // How many strips has each chamber?
0150   RPCStripsRing::iterator it = this->begin();
0151   uint32_t ch1Det = it->second.m_detRawId;
0152   for (; it != this->end(); ++it) {
0153     if (det2stripNo.find(it->second.m_detRawId) == det2stripNo.end()) {
0154       det2stripNo[it->second.m_detRawId] = 1;  // Add new chamber to a map, set strip cnt to 1
0155     } else {
0156       ++det2stripNo[it->second.m_detRawId];  // Increase strip count of a chamber
0157     }
0158
0159     if (det2stripNo.size() != 1 && ch1Det == it->second.m_detRawId) {
0160       ++ch1EndStrips;
0161     }
0162   }
0163
0164   det2stripNo[ch1Det] -= ch1EndStrips;
0165
0166   // std::cout << ch1BegStrips << " " << ch1EndStrips << std::endl;
0167
0168   //TDetId2StripNo::iterator itIds = det2stripNo.begin();
0169   //for(;itIds!=det2stripNo.end();++itIds){
0170   //    std::cout << itIds->first << " " << itIds->second << std::endl;
0171   //  }
0172
0173   it = this->begin();
0174   uint32_t lastDet = it->second.m_detRawId;
0175   while (it != this->end()) {
0176     if (det2stripNo[it->second.m_detRawId] < 0) {
0177       throw cms::Exception("RPCInternal") << " RPCStripsRing::filterOverlapingChambers() - no strips left \n";
0178     }
0179     if (it->second.m_detRawId == lastDet) {
0180       --det2stripNo[lastDet];
0181       ++it;
0182     } else if (det2stripNo[lastDet] == 0) {  // no more strips left in lastDet, proceed to new det
0183
0184       if (lastDet == ch1Det) {
0185         det2stripNo[ch1Det] += ch1EndStrips;
0186       }
0187
0188       lastDet = it->second.m_detRawId;
0189       --det2stripNo[lastDet];
0190       ++it;
0191     } else {  // there are still strips in last det, delete current strip
0192       --det2stripNo[it->second.m_detRawId];
0193       RPCStripsRing::iterator itErase = it;
0194       ++it;
0195       //std::cout << "Removing strip " <<  it->second.m_detRawId << " " << (int)it->second.m_strip << std::endl;
0196       this->erase(itErase);
0197     }
0198   }
0199 }
0200
0201 void RPCStripsRing::fillWithVirtualStrips() {
0202   if (m_didVirtuals)
0203     return;
0204   m_didVirtuals = true;
0205
0206   const float pi = 3.141592654;
0207   double dphi = 2.0 * pi / 1152;  // defines angular granulation of strips.
0208
0209   RPCStripsRing stripsToInsert;
0210
0211   float delta = 0;
0212   int stripsToAdd = 0;
0213
0214   RPCStripsRing::iterator it = this->begin();
0215   RPCStripsRing::iterator itLast = this->begin();
0216   for (; it != this->end(); ++it) {
0217     /*std::cout << it->first << " "
0218         << it->second.m_detRawId << " "
0219         << (int)it->second.m_strip << std::endl;
0220     */
0221
0222     delta = it->first - itLast->first;
0223     if (it == itLast ||                                        // skip first loop iteration
0224         itLast->second.m_detRawId == it->second.m_detRawId ||  // insert strips between two chambers only
0225         delta < 0) {
0226       itLast = it;
0227       continue;
0228     }
0229
0230     stripsToAdd = (int)std::floor(delta / dphi) - 1;
0231     //std::cout << delta << " " << stripsToAdd << std::endl;
0232
0233     if (isReferenceRing() && m_hwPlane == 6)
0234       ++stripsToAdd;
0235
0236     for (int i = 0; i < stripsToAdd; ++i) {
0237       stripsToInsert[itLast->first + dphi * (i + 1)] = TStrip();
0238     }
0239
0240     itLast = it;
0241   }
0242   // TODO: check delta between first and last strip in map
0243
0244   this->insert(stripsToInsert.begin(), stripsToInsert.end());
0245 }
0246 void RPCStripsRing::createRefConnections(TOtherConnStructVec& otherRings, int logplane, int logplaneSize) {
0247   //*
0248   /*std::cout << "RefCon for " << getRingId()
0249        << " (" << getEtaPartition()<<  ")"
0250        << " tower: " << getTowerForRefRing()
0251        << " ; connected: "
0252        << otherRings.size()
0253        << std::endl
0254        << std::endl;
0255   //*/
0256
0257   // XXX - TODO: warning on wrong logplaneSize
0258
0259   if (!this->isReferenceRing()) {
0260     throw cms::Exception("RPCInternal") << " RPCStripsRing::createRefConnections "
0261                                         << " called for non-reference ring \n";
0262   }
0263
0264   /*
0265    if (logplaneSize!=8) {
0266      throw cms::Exception("RPCInternal") << " RPCStripsRing::createRefConnections "
0267          << " called for lpSize " << logplaneSize << " \n";
0268
0269    }*/
0270   const float pi = 3.141592654;
0271   const float offset = (5. / 360.) * 2 * pi;  // XXX
0272
0273   //find first reference strip of first PAC (the strip with phi ~= 5deg)
0274   RPCStripsRing::iterator starEndIt = this->begin();
0275   while ((++starEndIt)->first < offset)
0276     ;
0277
0278   RPCStripsRing::iterator it = starEndIt;
0279   //--starEndIt;
0280
0281   float angle = 0;
0282   int curPACno = -1;
0283   int curStripNo = 0;
0284   int curBegStripNo = 0;
0285
0286   bool firstIter = true;
0287
0288   while (it != starEndIt || firstIter) {  // iterate over strips
0289
0290     firstIter = false;
0291     // New PAC
0292     if (curStripNo % logplaneSize == 0) {
0293       ++curPACno;
0294       curBegStripNo = curStripNo;
0295       RPCStripsRing::iterator plus8 = it;
0296       bool skipOccured = false;
0297       for (int i = 0; i < 7; ++i) {
0298         ++plus8;
0299         if (plus8 == this->end()) {
0300           plus8 = this->begin();
0301           skipOccured = true;
0302         }
0303       }
0304
0305       // calculate angle
0306       float phi = it->first;
0307       float phiP8 = plus8->first;
0308       if (skipOccured) {
0309         // phiP8 is negative
0310         // phi is positive
0311         // xcheck
0312         if (phi * phiP8 > 0) {
0313           throw cms::Exception("RPCInternal") << " RPCStripsRing::createRefConnections phi/phi8 error \n";
0314         }
0315         angle = (2 * pi + phiP8 + phi) / 2;
0316         if (angle > pi) {  // should land on positive side
0317           angle -= 2 * pi;
0318         }
0319
0320         if (std::abs(angle) > pi) {
0321           throw cms::Exception("RPCInternal") << " RPCStripsRing::createRefConnections "
0322                                               << " problem with angle calc \n";
0323         }
0324       } else {
0325         angle = (phiP8 + phi) / 2;
0326       }
0327       //std::cout << curPACno << " " << phiP8 << " " << phi << " "  << angle << std::endl;
0328
0329       TOtherConnStructVec::iterator itOt = otherRings.begin();
0330       for (; itOt != otherRings.end(); ++itOt) {
0331         itOt->m_it->second.createOtherConnections(
0332             getTowerForRefRing(), curPACno, itOt->m_logplane, itOt->m_logplaneSize, angle);
0333       }
0334     }
0335
0336     if (!it->second.isVirtual()) {
0337       L1RPCConeBuilder::TStripCon newCon;
0338       newCon.m_tower = getTowerForRefRing();
0339       newCon.m_PAC = curPACno;
0340       newCon.m_logplane = logplane;
0341       newCon.m_logstrip = curStripNo - curBegStripNo;
0342       //std::cout << " Adding con for " << it->second.m_detRawId << std::endl;
0343       (*m_connectionsMap)[it->second.m_detRawId][it->second.m_strip].push_back(newCon);
0344       //std::cout << " Adding ref connection " << std::endl;
0345     }
0346     ++curStripNo;
0347     ++it;
0348     if (it == this->end()) {
0349       it = this->begin();
0350     }
0351
0352   }  // iteration over strips ends
0353
0354   //std::cout << " refcon: " << curPACno << " PACs" << std::endl;
0355   //std::cout << "After refCon: " << m_connectionsMap.size() << std::endl;
0356 }
0357
0358 void RPCStripsRing::createOtherConnections(int tower, int PACno, int logplane, int logplaneSize, float angle) {
0359   //std::cout << "    OtherCon for " << getRingId() << std::endl;
0360
0361   if (this->isReferenceRing()) {
0362     throw cms::Exception("RPCInternal") << " RPCStripsRing::createOtherConnections "
0363                                         << " called for reference ring \n";
0364   }
0365
0366   RPCStripsRing::const_iterator it = this->lower_bound(angle);
0367
0368   if (it == this->end())
0369     it = this->begin();
0370
0371   for (int i = 0; i < logplaneSize / 2; i++) {
0372     if (it == this->begin())
0373       it = this->end();  // (m_stripPhiMap.end()--) is ok.
0374
0375     --it;
0376   }
0377
0378   for (int i = 0; i < logplaneSize; i++) {
0379     if (!it->second.isVirtual()) {
0380       L1RPCConeBuilder::TStripCon newCon;
0381       newCon.m_tower = tower;
0382       newCon.m_PAC = PACno;
0383       newCon.m_logplane = logplane;
0384       newCon.m_logstrip = i;
0385       (*m_connectionsMap)[it->second.m_detRawId][it->second.m_strip].push_back(newCon);
0386       //std::cout << " Adding other connection " << std::endl;
0387     }
0388
0389     ++it;
0390     if (it == this->end())
0391       it = this->begin();
0392   }
0393 }
0394
0395 // Defines to which tower this ring (only ref ring) belongs
0396 int RPCStripsRing::getTowerForRefRing() {
0397   int ret = 0;
0398
0399   if (!this->isReferenceRing()) {
0400     throw cms::Exception("RPCInternal") << " RPCStripsRing::getTowerForRefRing() "
0401                                         << " called for non reference ring \n";
0402   }
0403
0404   int etaAbs = std::abs(getEtaPartition());
0405   if (etaAbs < 8) {
0406     ret = getEtaPartition();
0407   } else if (etaAbs > 8) {
0408     int sign = (getEtaPartition() > 0 ? 1 : -1);
0409     ret = getEtaPartition() - sign;
0410   } else {
0411     throw cms::Exception("RPCInternal") << " RPCStripsRing::getTowerForRefRing() "
0412                                         << " called for etaPartition 8 \n";
0413   }
0414
0415   return ret;
0416 }
0417 /*
0418       struct TStripCon{
0419         signed char m_tower;
0420         unsigned char m_PAC;
0421         unsigned char m_logplane;
0422         unsigned char m_logstrip;
0423       };
0424       typedef std::vector<TStripCon> TStripConVec;
0425       typedef std::map<unsigned char, TStripConVec> TStrip2ConVec;
0426       typedef std::map<uint32_t, TStrip2ConVec> TConMap;
0427
0428       // compressed connections
0429       struct TCompressedCon{
0430         signed char m_tower;
0431         unsigned char m_PAC;
0432         signed char m_offset;
0433         signed char m_mul;
0434       };
0435       typedef std::vector<TCompressedCon> TCompressedConVec;
0436       typedef std::map<uint32_t, TCompressedConVec> TCompressedConMap;
0437
0438 */
0439
0440 void RPCStripsRing::compressConnections() {
0441   L1RPCConeBuilder::TConMap::iterator itChamber = m_connectionsMap->begin();
0442
0443   auto uncompressedConsLeft = std::make_shared<L1RPCConeBuilder::TConMap>();
0444
0445   m_compressedConnectionMap = std::make_shared<L1RPCConeBuilder::TCompressedConMap>();
0446
0447   int compressedCons = 0, uncompressedConsBefore = 0, uncompressedConsAfter = 0;
0448
0449   //   int offsetMin =0, offsetMax =0;
0450
0451   for (; itChamber != m_connectionsMap->end(); ++itChamber) {
0452     uint32_t detId = itChamber->first;
0453
0454     for (L1RPCConeBuilder::TStrip2ConVec::iterator itStrip = itChamber->second.begin();
0455          itStrip != itChamber->second.end();
0456          ++itStrip) {
0457       // Iterate over strip Connections
0458       for (L1RPCConeBuilder::TStripConVec::iterator itConn = itStrip->second.begin(); itConn != itStrip->second.end();
0459            ++itConn) {
0460         // Check if this connection isn't allready present in the compressed map
0461         ++uncompressedConsBefore;
0462         bool alreadyDone = false;
0463         if (m_compressedConnectionMap->find(detId) != m_compressedConnectionMap->end()) {
0464           // iterate over the vec, check element by element
0465           for (L1RPCConeBuilder::TCompressedConVec::iterator itCompConn = (*m_compressedConnectionMap)[detId].begin();
0466                itCompConn != (*m_compressedConnectionMap)[detId].end();
0467                ++itCompConn) {
0468             if (itCompConn->m_tower == itConn->m_tower && itCompConn->m_PAC == itConn->m_PAC &&
0469                 itCompConn->m_logplane == itConn->m_logplane)  // connection allready compressed
0470             {
0471               alreadyDone = true;
0472
0473               int logStrip = itCompConn->m_mul * itStrip->first + itCompConn->m_offset;
0474               if (logStrip != itConn->m_logstrip) {
0475                 //copy the problematic connection to the "safe" map
0476                 (*uncompressedConsLeft)[detId][itStrip->first].push_back(*itConn);
0477                 ++uncompressedConsAfter;
0478                 edm::LogWarning("RPCTriggerConfig")
0479                     << " Compression failed for det " << detId << " strip " << (int)itStrip->first << " . Got "
0480                     << (int)logStrip << " expected " << (int)itConn->m_logstrip << std::endl;
0481               } else {
0482                 itCompConn->addStrip(itStrip->first);
0483               }
0484             }
0485           }  // compressed connection iteration end
0486         }
0487         //if (detId==637569977) std::cout << " Buld cons for strip " << (int)itStrip->first << std::endl;
0488
0489         if (!alreadyDone) {
0490           // find another strip contributing to the same PAC,tower,logplane
0491           L1RPCConeBuilder::TStrip2ConVec::iterator itStripOther = itStrip;
0492           ++itStripOther;
0493           bool otherStripFound = false;
0494           signed char mul = 1;
0495           for (; itStripOther != itChamber->second.end() && !otherStripFound; ++itStripOther) {
0496             for (L1RPCConeBuilder::TStripConVec::iterator itConnOther = itStripOther->second.begin();
0497                  itConnOther != itStripOther->second.end();
0498                  ++itConnOther) {
0499               if (itConnOther->m_tower == itConn->m_tower && itConnOther->m_PAC == itConn->m_PAC &&
0500                   itConnOther->m_logplane == itConn->m_logplane)  // connection to same PAC,logplane
0501               {
0502                 otherStripFound = true;
0503                 if ((itStripOther->first - itStrip->first) * (itConnOther->m_logstrip - itConn->m_logstrip) < 0) {
0504                   mul = -1;
0505                 }
0506                 break;
0507               }
0508             }  // otherConnections iter ends
0509           }    // otherStrip iter ends
0510
0511           /*
0512           if (itConn->m_tower==3 && itConn->m_PAC==73 && itConn->m_logplane==4 && detId==637569977){
0513             std::cout << " Buld cons for strip " << (int)itStrip->first;
0514             if (otherStripFound)
0515               std::cout << " other strip " << itStrip->first;
0516             else
0517               std::cout << " no other strip ";
0518
0519             std::cout << std::endl;
0520
0521         }*/
0522
0523           L1RPCConeBuilder::TCompressedCon nCompConn;
0524           nCompConn.m_tower = itConn->m_tower;
0525           nCompConn.m_PAC = itConn->m_PAC;
0526           nCompConn.m_logplane = itConn->m_logplane;
0527           nCompConn.m_mul = mul;
0528           nCompConn.m_offset = itConn->m_logstrip - mul * (signed short)(itStrip->first);
0529           nCompConn.addStrip(itStrip->first);
0530
0531           if (otherStripFound) {
0532           } else {
0533             //  uncompressedConsLeft[detId][itStrip->first].push_back(*itConn);
0534             //  ++uncompressedConsAfter;
0535           }
0536           (*m_compressedConnectionMap)[detId].push_back(nCompConn);
0537           ++compressedCons;
0538
0539         }  // if(!allreadyDone)
0540       }    // iterate on connections
0541     }      // iterate on strips
0542   }        // iterate on chambers
0543
0544   // 159 -87
0545   //std::cout << offsetMax << " TT " << offsetMin << std::endl;
0546
0547   edm::LogInfo("RPCTriggerConfig") << " Compressed: " << compressedCons << " "
0548                                    << sizeof(L1RPCConeBuilder::TCompressedCon)
0549                                    << " Uncompressed before: " << uncompressedConsBefore << " "
0550                                    << sizeof(L1RPCConeBuilder::TStripCon)
0551                                    << " Uncompressed after: " << uncompressedConsAfter << " "
0552                                    << sizeof(L1RPCConeBuilder::TStripCon);
0553   m_connectionsMap = uncompressedConsLeft;
0554 }