L1TTwinMux/src/RPCtoDTTranslator.cc

0001 //-------------------------------------------------
0002 //
0003 //   Class: RPCtoDTTranslator
0004 //
0005 //   RPCtoDTTranslator
0006 //
0007 //
0008 //   Author :
0009 //   G. Flouris               U Ioannina    Mar. 2015
0010 //   modifications: G Karathanasis U Athens
0011 //--------------------------------------------------
0012
0013 #include <iostream>
0014 #include <iomanip>
0015 #include <iterator>
0016 #include <cmath>
0017 #include <map>
0018
0019 #include "L1Trigger/L1TTwinMux/interface/RPCHitCleaner.h"
0020 #include "L1Trigger/L1TTwinMux/interface/RPCtoDTTranslator.h"
0021 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0022 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0023 #include "CondFormats/DataRecord/interface/RPCEMapRcd.h"
0024 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0025
0026 using namespace std;
0027
0028 RPCtoDTTranslator::RPCtoDTTranslator(RPCDigiCollection const& inrpcDigis) : m_rpcDigis{inrpcDigis} {}
0029
0030 namespace {
0031   constexpr int max_rpc_bx = 2;
0032   constexpr int min_rpc_bx = -2;
0033
0034   struct rpc_hit {
0035     int bx;
0036     int station;
0037     int sector;
0038     int wheel;
0039     RPCDetId detid;
0040     int strip;
0041     int roll;
0042     int layer;
0043     rpc_hit(int pbx, int pstation, int psector, int pwheel, RPCDetId pdet, int pstrip, int proll, int player)
0044         : bx(pbx),
0045           station(pstation),
0046           sector(psector),
0047           wheel(pwheel),
0048           detid(pdet),
0049           strip(pstrip),
0050           roll(proll),
0051           layer(player) {}
0052   };
0053
0054   //Need to shift the index so that index 0
0055   // corresponds to min_rpc_bx
0056   class BxToHit {
0057   public:
0058     BxToHit() : m_hits{} {}  //zero initializes
0059
0060     static bool outOfRange(int iBX) { return (iBX > max_rpc_bx or iBX < min_rpc_bx); }
0061
0062     int& operator[](int iBX) { return m_hits[iBX - min_rpc_bx]; }
0063
0064     size_t size() const { return m_hits.size(); }
0065
0066   private:
0067     std::array<int, max_rpc_bx - min_rpc_bx + 1> m_hits;
0068   };
0069 }  // namespace
0070
0071 void RPCtoDTTranslator::run(const RPCGeometry& rpcGeometry) {
0072   std::vector<L1MuDTChambPhDigi> l1ttma_out;
0073   std::vector<L1MuDTChambPhDigi> l1ttma_hits_out;
0074
0075   std::vector<rpc_hit> vrpc_hit_layer1, vrpc_hit_layer2, vrpc_hit_st3, vrpc_hit_st4;
0076
0077   ///Init structues
0078   for (auto chamber = m_rpcDigis.begin(); chamber != m_rpcDigis.end(); ++chamber) {
0079     RPCDetId detid = (*chamber).first;
0080     for (auto digi = (*chamber).second.first; digi != (*chamber).second.second; ++digi) {
0081       if (detid.region() != 0)
0082         continue;  //Region = 0 Barrel
0083       if (BxToHit::outOfRange(digi->bx()))
0084         continue;
0085       if (detid.layer() == 1)
0086         vrpc_hit_layer1.emplace_back(digi->bx(),
0087                                      detid.station(),
0088                                      detid.sector(),
0089                                      detid.ring(),
0090                                      detid,
0091                                      digi->strip(),
0092                                      detid.roll(),
0093                                      detid.layer());
0094       if (detid.station() == 3)
0095         vrpc_hit_st3.emplace_back(digi->bx(),
0096                                   detid.station(),
0097                                   detid.sector(),
0098                                   detid.ring(),
0099                                   detid,
0100                                   digi->strip(),
0101                                   detid.roll(),
0102                                   detid.layer());
0103       if (detid.layer() == 2)
0104         vrpc_hit_layer2.emplace_back(digi->bx(),
0105                                      detid.station(),
0106                                      detid.sector(),
0107                                      detid.ring(),
0108                                      detid,
0109                                      digi->strip(),
0110                                      detid.roll(),
0111                                      detid.layer());
0112       if (detid.station() == 4)
0113         vrpc_hit_st4.emplace_back(digi->bx(),
0114                                   detid.station(),
0115                                   detid.sector(),
0116                                   detid.ring(),
0117                                   detid,
0118                                   digi->strip(),
0119                                   detid.roll(),
0120                                   detid.layer());
0121     }
0122   }  ///for chamber
0123
0124   vector<int> vcluster_size;
0125   int cluster_id = -1;
0126   int itr = 0;
0127   //      int hits[5][4][12][2][5][3][100]= {{{{{{{0}}}}}}};
0128   std::map<RPCHitCleaner::detId_Ext, int> hits;
0129   int cluster_size = 0;
0130   for (auto chamber = m_rpcDigis.begin(); chamber != m_rpcDigis.end(); ++chamber) {
0131     RPCDetId detid = (*chamber).first;
0132     int strip_n1 = -10000;
0133     int bx_n1 = -10000;
0134     if (detid.region() != 0)
0135       continue;  //Region = 0 Barrel
0136     for (auto digi = (*chamber).second.first; digi != (*chamber).second.second; ++digi) {
0137       if (fabs(digi->bx()) > 3)
0138         continue;
0139       //Create cluster ids and store their size
0140       //if((digi->strip()+1!=strip_n1)|| digi->bx()!=bx_n1){
0141       if (abs(digi->strip() - strip_n1) != 1 || digi->bx() != bx_n1) {
0142         if (itr != 0)
0143           vcluster_size.push_back(cluster_size);
0144         cluster_size = 0;
0145         cluster_id++;
0146       }
0147       itr++;
0148       cluster_size++;
0149       ///hit belongs to cluster with clusterid
0150       //hits[(detid.ring()+2)][(detid.station()-1)][(detid.sector()-1)][(detid.layer()-1)][(digi->bx()+2)][detid.roll()-1][digi->strip()]= cluster_id ;
0151       RPCHitCleaner::detId_Ext tmp{detid, digi->bx(), digi->strip()};
0152       hits[tmp] = cluster_id;
0153       ///strip of i-1
0154       strip_n1 = digi->strip();
0155       bx_n1 = digi->bx();
0156     }
0157   }  ///for chamber
0158   vcluster_size.push_back(cluster_size);
0159
0160   for (int wh = -2; wh <= 2; wh++) {
0161     for (int sec = 1; sec <= 12; sec++) {
0162       for (int st = 1; st <= 4; st++) {
0163         int rpcbx = 0;
0164         std::vector<int> delta_phib;
0165         bool found_hits = false;
0166         std::vector<int> rpc2dt_phi, rpc2dt_phib;
0167         ///Loop over all combinations of layer 1 and 2.
0168         int itr1 = 0;
0169         for (unsigned int l1 = 0; l1 < vrpc_hit_layer1.size(); l1++) {
0170           RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer1[l1].detid, vrpc_hit_layer1[l1].bx, vrpc_hit_layer1[l1].strip};
0171           int id = hits[tmp];
0172           int phi1 = radialAngle(vrpc_hit_layer1[l1].detid, rpcGeometry, vrpc_hit_layer1[l1].strip);
0173           if (vcluster_size[id] == 2 && itr1 == 0) {
0174             itr1++;
0175             continue;
0176           }
0177           if (vcluster_size[id] == 2 && itr1 == 1) {
0178             itr1 = 0;
0179             phi1 = phi1 + (radialAngle(vrpc_hit_layer1[l1 - 1].detid, rpcGeometry, vrpc_hit_layer1[l1 - 1].strip));
0180             phi1 /= 2;
0181           }
0182           int itr2 = 0;
0183           for (unsigned int l2 = 0; l2 < vrpc_hit_layer2.size(); l2++) {
0184             if (vrpc_hit_layer1[l1].station != st || vrpc_hit_layer2[l2].station != st)
0185               continue;
0186             if (vrpc_hit_layer1[l1].sector != sec || vrpc_hit_layer2[l2].sector != sec)
0187               continue;
0188             if (vrpc_hit_layer1[l1].wheel != wh || vrpc_hit_layer2[l2].wheel != wh)
0189               continue;
0190             if (vrpc_hit_layer1[l1].bx != vrpc_hit_layer2[l2].bx)
0191               continue;
0192             RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer2[l2].detid, vrpc_hit_layer2[l2].bx, vrpc_hit_layer2[l2].strip};
0193             int id = hits[tmp];
0194
0195             if (vcluster_size[id] == 2 && itr2 == 0) {
0196               itr2++;
0197               continue;
0198             }
0199
0200             //int phi1 = radialAngle(vrpc_hit_layer1[l1].detid, c, vrpc_hit_layer1[l1].strip) ;
0201             int phi2 = radialAngle(vrpc_hit_layer2[l2].detid, rpcGeometry, vrpc_hit_layer2[l2].strip);
0202             if (vcluster_size[id] == 2 && itr2 == 1) {
0203               itr2 = 0;
0204               phi2 = phi2 + (radialAngle(vrpc_hit_layer2[l2 - 1].detid, rpcGeometry, vrpc_hit_layer2[l2 - 1].strip));
0205               phi2 /= 2;
0206             }
0207             int average = ((phi1 + phi2) / 2) << 2;  //10-bit->12-bit
0208             rpc2dt_phi.push_back(average);           //Convert and store to 12-bit
0209             //int xin = localX(vrpc_hit_layer1[l1].detid, c, vrpc_hit_layer1[l1].strip);
0210             //int xout = localX(vrpc_hit_layer2[l2].detid, c, vrpc_hit_layer2[l2].strip);
0211             //cout<<(phi1<<2)<<"   "<<l1<<"   "<<vrpc_hit_layer1[l1].station<<endl;
0212             //cout<<(phi2<<2)<<"   "<<l1<<"   "<<vrpc_hit_layer1[l1].station<<endl;
0213             int xin = localXX((phi1 << 2), 1, vrpc_hit_layer1[l1].station);
0214             int xout = localXX((phi2 << 2), 2, vrpc_hit_layer2[l2].station);
0215             if (vcluster_size[id] == 2 && itr2 == 1) {
0216               int phi1_n1 = radialAngle(vrpc_hit_layer1[l1 - 1].detid, rpcGeometry, vrpc_hit_layer1[l1 - 1].strip);
0217               int phi2_n1 = radialAngle(vrpc_hit_layer2[l2 - 1].detid, rpcGeometry, vrpc_hit_layer2[l2 - 1].strip);
0218               xin += localXX((phi1_n1 << 2), 1, vrpc_hit_layer1[l1].station);
0219               xout += localXX((phi2_n1 << 2), 2, vrpc_hit_layer2[l2].station);
0220               xin /= 2;
0221               xout /= 2;
0222             }
0223             //cout<<">>"<<xin<<"   "<<xout<<endl;
0224             int phi_b = bendingAngle(xin, xout, average);
0225             //cout<<"phib   "<<phi_b<<endl;
0226             rpc2dt_phib.push_back(phi_b);
0227             //delta_phib to find the highest pt primitve
0228             delta_phib.push_back(abs(phi_b));
0229             found_hits = true;
0230             rpcbx = vrpc_hit_layer1[l1].bx;
0231           }
0232         }
0233         if (found_hits) {
0234           //cout<<"found_hits"<<endl;
0235           int min_index = std::distance(delta_phib.begin(), std::min_element(delta_phib.begin(), delta_phib.end())) + 0;
0236           //cout<<min_index<<endl;
0237           l1ttma_out.emplace_back(rpcbx, wh, sec - 1, st, rpc2dt_phi[min_index], rpc2dt_phib[min_index], 3, 0, 0, 2);
0238         }
0239         ///Use ts2tag variable to store N rpchits for the same st/wheel/sec
0240         BxToHit hit;
0241         itr1 = 0;
0242         for (unsigned int l1 = 0; l1 < vrpc_hit_layer1.size(); l1++) {
0243           if (vrpc_hit_layer1[l1].station != st || st > 2 || vrpc_hit_layer1[l1].sector != sec ||
0244               vrpc_hit_layer1[l1].wheel != wh)
0245             continue;
0246           //int id = hits[vrpc_hit_layer1[l1].wheel+2][(vrpc_hit_layer1[l1].station-1)][(vrpc_hit_layer1[l1].sector-1)][(vrpc_hit_layer1[l1].layer-1)][(vrpc_hit_layer1[l1].bx+2)][vrpc_hit_layer1[l1].roll-1][vrpc_hit_layer1[l1].strip];
0247           RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer1[l1].detid, vrpc_hit_layer1[l1].bx, vrpc_hit_layer1[l1].strip};
0248           int id = hits[tmp];
0249           if (vcluster_size[id] == 2 && itr1 == 0) {
0250             itr1++;
0251             continue;
0252           }
0253           int phi2 = radialAngle(vrpc_hit_layer1[l1].detid, rpcGeometry, vrpc_hit_layer1[l1].strip);
0254           phi2 = phi2 << 2;
0255           if (vcluster_size[id] == 2 && itr1 == 1) {
0256             itr1 = 0;
0257             phi2 = phi2 + (radialAngle(vrpc_hit_layer1[l1 - 1].detid, rpcGeometry, vrpc_hit_layer1[l1 - 1].strip) << 2);
0258             phi2 /= 2;
0259           }
0260
0261           l1ttma_hits_out.emplace_back(
0262               vrpc_hit_layer1[l1].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_layer1[l1].bx], 0, 2);
0263           hit[vrpc_hit_layer1[l1].bx]++;
0264         }
0265         itr1 = 0;
0266         for (unsigned int l2 = 0; l2 < vrpc_hit_layer2.size(); l2++) {
0267           if (vrpc_hit_layer2[l2].station != st || st > 2 || vrpc_hit_layer2[l2].sector != sec ||
0268               vrpc_hit_layer2[l2].wheel != wh)
0269             continue;
0270           RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer2[l2].detid, vrpc_hit_layer2[l2].bx, vrpc_hit_layer2[l2].strip};
0271           int id = hits[tmp];
0272           //              int id = hits[vrpc_hit_layer2[l2].wheel+2][(vrpc_hit_layer2[l2].station-1)][(vrpc_hit_layer2[l2].sector-1)][(vrpc_hit_layer2[l2].layer-1)][(vrpc_hit_layer2[l2].bx+2)][vrpc_hit_layer2[l2].roll-1][vrpc_hit_layer2[l2].strip];
0273           if (vcluster_size[id] == 2 && itr1 == 0) {
0274             itr1++;
0275             continue;
0276           }
0277           int phi2 = radialAngle(vrpc_hit_layer2[l2].detid, rpcGeometry, vrpc_hit_layer2[l2].strip);
0278           phi2 = phi2 << 2;
0279           if (vcluster_size[id] == 2 && itr1 == 1) {
0280             itr1 = 0;
0281             phi2 = phi2 + (radialAngle(vrpc_hit_layer2[l2 - 1].detid, rpcGeometry, vrpc_hit_layer2[l2 - 1].strip) << 2);
0282             phi2 /= 2;
0283           }
0284           l1ttma_hits_out.emplace_back(
0285               vrpc_hit_layer2[l2].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_layer2[l2].bx], 0, 2);
0286           hit[vrpc_hit_layer2[l2].bx]++;
0287         }
0288         itr1 = 0;
0289
0290         for (unsigned int l1 = 0; l1 < vrpc_hit_st3.size(); l1++) {
0291           if (st != 3 || vrpc_hit_st3[l1].station != 3 || vrpc_hit_st3[l1].wheel != wh ||
0292               vrpc_hit_st3[l1].sector != sec)
0293             continue;
0294           RPCHitCleaner::detId_Ext tmp{vrpc_hit_st3[l1].detid, vrpc_hit_st3[l1].bx, vrpc_hit_st3[l1].strip};
0295           int id = hits[tmp];
0296           //int id = hits[vrpc_hit_st3[l1].wheel+2][(vrpc_hit_st3[l1].station-1)][(vrpc_hit_st3[l1].sector-1)][(vrpc_hit_st3[l1].layer-1)][(vrpc_hit_st3[l1].bx+2)][vrpc_hit_st3[l1].roll-1][vrpc_hit_st3[l1].strip];
0297           if (vcluster_size[id] == 2 && itr1 == 0) {
0298             itr1++;
0299             continue;
0300           }
0301           int phi2 = radialAngle(vrpc_hit_st3[l1].detid, rpcGeometry, vrpc_hit_st3[l1].strip);
0302           phi2 = phi2 << 2;
0303           if (vcluster_size[id] == 2 && itr1 == 1) {
0304             itr1 = 0;
0305             phi2 = phi2 + (radialAngle(vrpc_hit_st3[l1 - 1].detid, rpcGeometry, vrpc_hit_st3[l1 - 1].strip) << 2);
0306             phi2 /= 2;
0307           }
0308           l1ttma_hits_out.emplace_back(
0309               vrpc_hit_st3[l1].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_st3[l1].bx], 0, 2);
0310           hit[vrpc_hit_st3[l1].bx]++;
0311         }
0312         itr1 = 0;
0313
0314         for (unsigned int l1 = 0; l1 < vrpc_hit_st4.size(); l1++) {
0315           if (st != 4 || vrpc_hit_st4[l1].station != 4 || vrpc_hit_st4[l1].wheel != wh ||
0316               vrpc_hit_st4[l1].sector != sec)
0317             continue;
0318           //int id = hits[vrpc_hit_st4[l1].wheel+2][(vrpc_hit_st4[l1].station-1)][(vrpc_hit_st4[l1].sector-1)][(vrpc_hit_st4[l1].layer-1)][(vrpc_hit_st4[l1].bx+2)][vrpc_hit_st4[l1].roll-1][vrpc_hit_st4[l1].strip];
0319           RPCHitCleaner::detId_Ext tmp{vrpc_hit_st4[l1].detid, vrpc_hit_st4[l1].bx, vrpc_hit_st4[l1].strip};
0320           int id = hits[tmp];
0321           if (vcluster_size[id] == 2 && itr1 == 0) {
0322             itr1++;
0323             continue;
0324           }
0325           int phi2 = radialAngle(vrpc_hit_st4[l1].detid, rpcGeometry, vrpc_hit_st4[l1].strip);
0326           phi2 = phi2 << 2;
0327           if (vcluster_size[id] == 2 && itr1 == 1) {
0328             itr1 = 0;
0329             phi2 = phi2 + (radialAngle(vrpc_hit_st4[l1 - 1].detid, rpcGeometry, vrpc_hit_st4[l1 - 1].strip) << 2);
0330             phi2 /= 2;
0331           }
0332           l1ttma_hits_out.emplace_back(
0333               vrpc_hit_st4[l1].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_st4[l1].bx], 0, 2);
0334           hit[vrpc_hit_st4[l1].bx]++;
0335           //l1ttma_out.push_back(rpc2dt_out);
0336
0337           //break;
0338         }
0339       }
0340     }
0341   }
0342   ///Container to store RPC->DT for RPC only (only in stations 1 and 2 (2 layers->phib))
0343   m_rpcdt_translated.setContainer(l1ttma_out);
0344   ///Container to store RPC->DT for Bx correction
0345   m_rpchitsdt_translated.setContainer(l1ttma_hits_out);
0346 }
0347
0348 ///function - will be replaced by LUTs(?)
0349 int RPCtoDTTranslator::radialAngle(RPCDetId detid, const RPCGeometry& rpcGeometry, int strip) {
0350   int radialAngle;
0351   // from phiGlobal to radialAngle of the primitive in sector sec in [1..12] <- RPC scheme
0352
0353   const RPCRoll* roll = rpcGeometry.roll(detid);
0354   GlobalPoint stripPosition = roll->toGlobal(roll->centreOfStrip(strip));
0355
0356   double globalphi = stripPosition.phi();
0357   int sector = (roll->id()).sector();
0358   if (sector == 1)
0359     radialAngle = int(globalphi * 1024);
0360   else {
0361     if (globalphi >= 0)
0362       radialAngle = int((globalphi - (sector - 1) * Geom::pi() / 6.) * 1024);
0363     else
0364       radialAngle = int((globalphi + (13 - sector) * Geom::pi() / 6.) * 1024);
0365   }
0366   return radialAngle;
0367 }
0368
0369 ///function - will be replaced by LUTs(?)
0370 int RPCtoDTTranslator::localX(RPCDetId detid, const RPCGeometry& rpcGeometry, int strip) {
0371   const RPCRoll* roll = rpcGeometry.roll(detid);
0372
0373   ///Orientaion of RPCs
0374   GlobalPoint p1cmPRG = roll->toGlobal(LocalPoint(1, 0, 0));
0375   GlobalPoint m1cmPRG = roll->toGlobal(LocalPoint(-1, 0, 0));
0376   float phip1cm = p1cmPRG.phi();
0377   float phim1cm = m1cmPRG.phi();
0378   int direction = (phip1cm - phim1cm) / abs(phip1cm - phim1cm);
0379   ///---Orientaion
0380
0381   return direction * roll->centreOfStrip(strip).x();
0382 }
0383
0384 int RPCtoDTTranslator::bendingAngle(int xin, int xout, int phi) {
0385   // use chamber size and max angle in hw units 512
0386   int atanv = (int)(atan((xout - xin) / 34.6) * 512);
0387   if (atanv > 512)
0388     return 512;
0389   int rvalue = atanv - phi / 8;
0390   return rvalue;
0391 }
0392
0393 int RPCtoDTTranslator::localXX(int phi, int layer, int station) {
0394   //R[stat][layer] - radius of rpc station/layer from center of CMS
0395   double R[2][2] = {{410.0, 444.8}, {492.7, 527.3}};
0396   double rvalue = R[station - 1][layer - 1] * tan(phi / 4096.);
0397   return rvalue;
0398 }