Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:21:21

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