Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:12:00 

0001 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0002 
0003 #include "L1Trigger/L1TCalorimeter/interface/CaloStage2Nav.h"
0004 
0005 const l1t::CaloTower l1t::CaloTools::nullTower_;
0006 const l1t::CaloCluster l1t::CaloTools::nullCluster_;
0007 
0008 const float l1t::CaloTools::kGTEtaLSB = 0.0435;
0009 const float l1t::CaloTools::kGTPhiLSB = 0.0435;
0010 const float l1t::CaloTools::kGTEtLSB = 0.5;
0011 
0012 const int64_t l1t::CaloTools::cos_coeff[72] = {
0013     1023,  1019,  1007,  988,  961,  927,  886,  838,  784,  723,  658,  587,  512,  432,  350,  265,  178,   89,
0014     0,     -89,   -178,  -265, -350, -432, -512, -587, -658, -723, -784, -838, -886, -927, -961, -988, -1007, -1019,
0015     -1023, -1019, -1007, -988, -961, -927, -886, -838, -784, -723, -658, -587, -512, -432, -350, -265, -178,  -89,
0016     0,     89,    178,   265,  350,  432,  511,  587,  658,  723,  784,  838,  886,  927,  961,  988,  1007,  1019};
0017 
0018 const int64_t l1t::CaloTools::sin_coeff[72] = {
0019     0,     89,    178,   265,  350,  432,  512,  587,  658,  723,  784,  838,  886,  927,  961,  988,  1007,  1019,
0020     1023,  1019,  1007,  988,  961,  927,  886,  838,  784,  723,  658,  587,  512,  432,  350,  265,  178,   89,
0021     0,     -89,   -178,  -265, -350, -432, -512, -587, -658, -723, -784, -838, -886, -927, -961, -988, -1007, -1019,
0022     -1023, -1019, -1007, -988, -961, -927, -886, -838, -784, -723, -658, -587, -512, -432, -350, -265, -178,  -89};
0023 
0024 // mapping between sums in emulator and data
0025 const int l1t::CaloTools::emul_to_data_sum_index_map[31] = {
0026     9,  1,  19, 8,  0,  18, 10, 4,  6,  14,  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
0027     28, 24, 13, 27, 23, 15, 5,  7,  22, 12,  // 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
0028     3,  21, 11, 2,  20, 17, 30, 26, 16, 29,  // 20, 21, 22, 23, 24, 25, 26, 27, 28, 29
0029     25                                       // 30, 31
0030 };
0031 
0032 bool l1t::CaloTools::insertTower(std::vector<l1t::CaloTower>& towers, const l1t::CaloTower& tower) {
0033   size_t towerIndex = CaloTools::caloTowerHash(tower.hwEta(), tower.hwPhi());
0034   if (towers.size() > towerIndex) {
0035     towers.at(towerIndex) = tower;
0036     return true;
0037   } else
0038     return false;
0039 }
0040 
0041 //currently implemented as a brute force search but this will hopefully change in the future
0042 //with standarising the layout of std::vector<l1t::CaloTower>
0043 const l1t::CaloTower& l1t::CaloTools::getTower(const std::vector<l1t::CaloTower>& towers, int iEta, int iPhi) {
0044   if (abs(iEta) > CaloTools::kHFEnd)
0045     return nullTower_;
0046 
0047   size_t towerIndex = CaloTools::caloTowerHash(iEta, iPhi);
0048   if (towerIndex < towers.size()) {
0049     if (towers[towerIndex].hwEta() != iEta ||
0050         towers[towerIndex].hwPhi() !=
0051             iPhi) {  //it failed, this is bad, but we will not log the error due to policy and silently attempt to do a brute force search instead
0052       //std::cout <<"error, tower "<<towers[towerIndex].hwEta()<<" "<<towers[towerIndex].hwPhi()<<" does not match "<<iEta<<" "<<iPhi<<" index "<<towerIndex<<" nr towrs "<<towers.size()<<std::endl;
0053       for (size_t towerNr = 0; towerNr < towers.size(); towerNr++) {
0054         if (towers[towerNr].hwEta() == iEta && towers[towerNr].hwPhi() == iPhi)
0055           return towers[towerNr];
0056       }
0057     } else
0058       return towers[towerIndex];
0059 
0060   } else {  // in case the vector of towers do not contain all the towers (towerIndex can be > towers.size())
0061     for (size_t towerNr = 0; towerNr < towers.size(); towerNr++) {
0062       if (towers[towerNr].hwEta() == iEta && towers[towerNr].hwPhi() == iPhi)
0063         return towers[towerNr];
0064     }
0065   }
0066 
0067   return nullTower_;
0068 }
0069 
0070 const l1t::CaloCluster& l1t::CaloTools::getCluster(const std::vector<l1t::CaloCluster>& clusters, int iEta, int iPhi) {
0071   for (size_t clusterNr = 0; clusterNr < clusters.size(); clusterNr++) {
0072     if (clusters[clusterNr].hwEta() == iEta && clusters[clusterNr].hwPhi() == iPhi)
0073       return clusters[clusterNr];
0074   }
0075   return nullCluster_;
0076 }
0077 
0078 //this implimentation has not all the necessary info yet, we need to check the exact HF numbering
0079 //(iEta=-28,iPhi=1)=index 0 to (iEta=28,iPhi=72)=index 28*72*2-1
0080 //HF then runs after that so -32,1 = 28*72*2
0081 size_t l1t::CaloTools::caloTowerHash(int iEta, int iPhi) {
0082   if (!isValidIEtaIPhi(iEta, iPhi))
0083     return caloTowerHashMax();
0084   else {
0085     const int absIEta = abs(iEta);
0086     if (absIEta > kHFEnd)
0087       return kNrTowers;
0088     else if (absIEta <= kHBHEEnd) {  //HBHE
0089       int iEtaNoZero = iEta;
0090       if (iEta > 0)
0091         iEtaNoZero--;
0092       return (iEtaNoZero + kHBHEEnd) * kHBHENrPhi + iPhi - 1;
0093     } else {                          //HF
0094       int iEtaIndex = iEta + kHFEnd;  //iEta=-32 is 0
0095       if (iEta > 0)
0096         iEtaIndex = iEta - kHBHEEnd + (kHFEnd - kHBHEEnd) - 1;  //but iEta=29 is 4
0097       return iEtaIndex * kHFNrPhi + iPhi / kHFPhiSeg + kNrHBHETowers;
0098     }
0099   }
0100 }
0101 
0102 size_t l1t::CaloTools::caloTowerHashMax() { return kNrTowers; }
0103 
0104 bool l1t::CaloTools::isValidIEtaIPhi(int iEta, int iPhi) {
0105   size_t absIEta = abs(iEta);
0106   if (iPhi <= 0 || iPhi > kNPhi)
0107     return false;
0108   if (absIEta == 0 || absIEta > kHFEnd)
0109     return false;
0110   //if(absIEta>kHBHEEnd && iPhi%kHFPhiSeg!=1) return false;
0111   return true;
0112 }
0113 
0114 int l1t::CaloTools::calHwEtSum(int iEta,
0115                                int iPhi,
0116                                const std::vector<l1t::CaloTower>& towers,
0117                                int localEtaMin,
0118                                int localEtaMax,
0119                                int localPhiMin,
0120                                int localPhiMax,
0121                                SubDet etMode) {
0122   return calHwEtSum(iEta, iPhi, towers, localEtaMin, localEtaMax, localPhiMin, localPhiMax, kHFEnd, etMode);
0123 }
0124 
0125 int l1t::CaloTools::calHwEtSum(int iEta,
0126                                int iPhi,
0127                                const std::vector<l1t::CaloTower>& towers,
0128                                int localEtaMin,
0129                                int localEtaMax,
0130                                int localPhiMin,
0131                                int localPhiMax,
0132                                int iEtaAbsMax,
0133                                SubDet etMode) {
0134   int hwEtSum = 0;
0135   for (int etaNr = localEtaMin; etaNr <= localEtaMax; etaNr++) {
0136     for (int phiNr = localPhiMin; phiNr <= localPhiMax; phiNr++) {
0137       int towerIEta = l1t::CaloStage2Nav::offsetIEta(iEta, etaNr);
0138       int towerIPhi = l1t::CaloStage2Nav::offsetIPhi(iPhi, phiNr);
0139       if (abs(towerIEta) <= iEtaAbsMax) {
0140         const l1t::CaloTower& tower = getTower(towers, towerIEta, towerIPhi);
0141         if (etMode == ECAL)
0142           hwEtSum += tower.hwEtEm();
0143         else if (etMode == HCAL)
0144           hwEtSum += tower.hwEtHad();
0145         else if (etMode == CALO)
0146           hwEtSum += tower.hwPt();
0147       }
0148     }
0149   }
0150   return hwEtSum;
0151 }
0152 
0153 size_t l1t::CaloTools::calNrTowers(int iEtaMin,
0154                                    int iEtaMax,
0155                                    int iPhiMin,
0156                                    int iPhiMax,
0157                                    const std::vector<l1t::CaloTower>& towers,
0158                                    int minHwEt,
0159                                    int maxHwEt,
0160                                    SubDet etMode) {
0161   size_t nrTowers = 0;
0162   l1t::CaloStage2Nav nav(iEtaMin, iPhiMin);
0163   while (nav.currIEta() <= iEtaMax) {
0164     bool finishPhi = false;
0165     while (!finishPhi) {
0166       const l1t::CaloTower& tower =
0167           l1t::CaloTools::getTower(towers, CaloTools::caloEta(nav.currIEta()), nav.currIPhi());
0168       int towerHwEt = 0;
0169       if (etMode == ECAL)
0170         towerHwEt += tower.hwEtEm();
0171       else if (etMode == HCAL)
0172         towerHwEt += tower.hwEtHad();
0173       else if (etMode == CALO)
0174         towerHwEt += tower.hwPt();
0175       if (towerHwEt >= minHwEt && towerHwEt <= maxHwEt)
0176         nrTowers++;
0177       finishPhi = (nav.currIPhi() == iPhiMax);
0178       nav.north();
0179     }
0180     nav.east();
0181     nav.resetIPhi();
0182   }
0183   return nrTowers;
0184 }
0185 
0186 std::pair<float, float> l1t::CaloTools::towerEtaBounds(int ieta) {
0187   if (ieta == 0)
0188     ieta = 1;
0189   if (ieta > kHFEnd)
0190     ieta = kHFEnd;
0191   if (ieta < (-1 * kHFEnd))
0192     ieta = -1 * kHFEnd;
0193   //const float towerEtas[33] = {0,0.087,0.174,0.261,0.348,0.435,0.522,0.609,0.696,0.783,0.870,0.957,1.044,1.131,1.218,1.305,1.392,1.479,1.566,1.653,1.740,1.830,1.930,2.043,2.172,2.322,2.5,2.650,3.000,3.5,4.0,4.5,5.0};
0194   const float towerEtas[42] = {0,     0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870,
0195                                0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830,
0196                                1.930, 2.043, 2.172, 2.322, 2.5,   2.650, 2.853, 3.139, 3.314, 3.489, 3.664,
0197                                3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191, 5.191};
0198   return std::make_pair(towerEtas[abs(ieta) - 1], towerEtas[abs(ieta)]);
0199 }
0200 
0201 float l1t::CaloTools::towerEta(int ieta) {
0202   std::pair<float, float> bounds = towerEtaBounds(ieta);
0203   float eta = (bounds.second + bounds.first) / 2.;
0204   float sign = ieta > 0 ? 1. : -1.;
0205   return sign * eta;
0206 }
0207 
0208 float l1t::CaloTools::towerPhi(int ieta, int iphi) {
0209   float phi = (float(iphi) - 0.5) * towerPhiSize(ieta);
0210   if (phi > M_PI)
0211     phi = phi - (2 * M_PI);
0212   return phi;
0213 }
0214 
0215 float l1t::CaloTools::towerEtaSize(int ieta) {
0216   std::pair<float, float> bounds = towerEtaBounds(ieta);
0217   float size = (bounds.second - bounds.first);
0218   return size;
0219 }
0220 
0221 float l1t::CaloTools::towerPhiSize(int ieta) { return 2. * M_PI / kNPhi; }
0222 
0223 // convert from calo ieta to internal MP ieta
0224 int l1t::CaloTools::mpEta(int ieta) {
0225   if (ieta > kHFBegin)
0226     return ieta - 1;
0227   else if (ieta < -1 * kHFBegin)
0228     return ieta + 1;
0229   else
0230     return ieta;
0231 }
0232 
0233 // convert from internal MP ieta to calo ieta
0234 int l1t::CaloTools::caloEta(int mpEta) {
0235   if (mpEta >= kHFBegin)
0236     return mpEta + 1;
0237   else if (mpEta <= -1 * kHFBegin)
0238     return mpEta - 1;
0239   else
0240     return mpEta;
0241 }
0242 
0243 // convert calorimeter ieta to RCT region index
0244 int l1t::CaloTools::regionEta(int ieta) {
0245   // outside HF
0246   if (abs(ieta) > kHFEnd)
0247     return (ieta < 0 ? 0 : 21);
0248 
0249   // inside HBHE
0250   if (abs(ieta) <= kHFBegin) {
0251     if (ieta < 0)
0252       return 11 - ceil(double(abs(ieta) / 4.));
0253     else
0254       return ceil(double(abs(ieta) / 4.)) + 10;
0255   }
0256 
0257   // in HF
0258   if (ieta < 0)
0259     return 4 - ceil(double(abs(ieta) - 29) / 4.);
0260   else
0261     return ceil(double(abs(ieta) - 29) / 4.) + 17;
0262 }
0263 
0264 // convert calorimeter ieta to etaBin16 index
0265 int l1t::CaloTools::bin16Eta(int ieta) {
0266   int absIEta = abs(ieta);
0267 
0268   if (absIEta > 0 && absIEta <= 5)
0269     return 0;
0270   else if (absIEta <= 9)
0271     return 1;
0272   else if (absIEta <= 13)
0273     return 2;
0274   else if (absIEta <= 15)
0275     return 3;
0276   else if (absIEta <= 17)
0277     return 4;
0278   else if (absIEta <= 19)
0279     return 5;
0280   else if (absIEta <= 21)
0281     return 6;
0282   else if (absIEta == 22)
0283     return 7;
0284   else if (absIEta == 23)
0285     return 8;
0286   else if (absIEta == 24)
0287     return 9;
0288   else if (absIEta == 25)
0289     return 10;
0290   else if (absIEta == 26)
0291     return 11;
0292   else if (absIEta <= 28)
0293     return 12;
0294   else if (absIEta <= 32)
0295     return 13;
0296   else if (absIEta <= 36)
0297     return 14;
0298   else if (absIEta <= 41)
0299     return 15;
0300   else
0301     return -1;  // error
0302 }
0303 
0304 int l1t::CaloTools::gtEta(int ieta) {
0305   double eta = towerEta(ieta);
0306   return round(eta / kGTEtaLSB);
0307 }
0308 
0309 int l1t::CaloTools::gtPhi(int ieta, int iphi) {
0310   double phi = towerPhi(ieta, iphi);
0311   if (phi < 0)
0312     phi = phi + 2 * M_PI;
0313   return round(phi / kGTPhiLSB);
0314 }
0315 
0316 // this conversion is based on GT input definitions in CMS DN-2014/029
0317 math::PtEtaPhiMLorentzVector l1t::CaloTools::p4Demux(l1t::L1Candidate* cand) {
0318   return math::PtEtaPhiMLorentzVector(
0319       cand->hwPt() * kGTEtLSB + 1.E-6, cand->hwEta() * kGTEtaLSB, cand->hwPhi() * kGTPhiLSB, 0.);
0320 }
0321 
0322 l1t::EGamma l1t::CaloTools::egP4Demux(l1t::EGamma& eg) {
0323   l1t::EGamma tmpEG(p4Demux(&eg), eg.hwPt(), eg.hwEta(), eg.hwPhi(), eg.hwQual(), eg.hwIso());
0324   tmpEG.setTowerIPhi(eg.towerIPhi());
0325   tmpEG.setTowerIEta(eg.towerIEta());
0326   tmpEG.setRawEt(eg.rawEt());
0327   tmpEG.setIsoEt(eg.isoEt());
0328   tmpEG.setFootprintEt(eg.footprintEt());
0329   tmpEG.setNTT(eg.nTT());
0330   tmpEG.setShape(eg.shape());
0331   tmpEG.setTowerHoE(eg.towerHoE());
0332 
0333   return tmpEG;
0334 }
0335 
0336 l1t::Tau l1t::CaloTools::tauP4Demux(l1t::Tau& tau) {
0337   l1t::Tau tmpTau(p4Demux(&tau), tau.hwPt(), tau.hwEta(), tau.hwPhi(), tau.hwQual(), tau.hwIso());
0338   tmpTau.setTowerIPhi(tau.towerIPhi());
0339   tmpTau.setTowerIEta(tau.towerIEta());
0340   tmpTau.setRawEt(tau.rawEt());
0341   tmpTau.setIsoEt(tau.isoEt());
0342   tmpTau.setNTT(tau.nTT());
0343   tmpTau.setHasEM(tau.hasEM());
0344   tmpTau.setIsMerged(tau.isMerged());
0345 
0346   return tmpTau;
0347 }
0348 
0349 l1t::Jet l1t::CaloTools::jetP4Demux(l1t::Jet& jet) {
0350   l1t::Jet tmpJet(p4Demux(&jet), jet.hwPt(), jet.hwEta(), jet.hwPhi(), jet.hwQual());
0351   tmpJet.setTowerIPhi(jet.towerIPhi());
0352   tmpJet.setTowerIEta(jet.towerIEta());
0353   tmpJet.setRawEt(jet.rawEt());
0354   tmpJet.setSeedEt(jet.seedEt());
0355   tmpJet.setPUEt(jet.puEt());
0356   tmpJet.setPUDonutEt(0, jet.puDonutEt(0));
0357   tmpJet.setPUDonutEt(1, jet.puDonutEt(1));
0358   tmpJet.setPUDonutEt(2, jet.puDonutEt(2));
0359   tmpJet.setPUDonutEt(3, jet.puDonutEt(3));
0360 
0361   return tmpJet;
0362 }
0363 
0364 l1t::EtSum l1t::CaloTools::etSumP4Demux(l1t::EtSum& etsum) {
0365   return l1t::EtSum(p4Demux(&etsum), etsum.getType(), etsum.hwPt(), etsum.hwEta(), etsum.hwPhi(), etsum.hwQual());
0366 }
0367 
0368 //
0369 math::PtEtaPhiMLorentzVector l1t::CaloTools::p4MP(l1t::L1Candidate* cand) {
0370   return math::PtEtaPhiMLorentzVector(
0371       cand->hwPt() * 0.5 + 1.E-6, towerEta(cand->hwEta()), towerPhi(cand->hwEta(), cand->hwPhi()), 0.);
0372 }
0373 
0374 l1t::EGamma l1t::CaloTools::egP4MP(l1t::EGamma& eg) {
0375   l1t::EGamma tmpEG(p4MP(&eg), eg.hwPt(), eg.hwEta(), eg.hwPhi(), eg.hwQual(), eg.hwIso());
0376   tmpEG.setTowerIPhi(eg.towerIPhi());
0377   tmpEG.setTowerIEta(eg.towerIEta());
0378   tmpEG.setRawEt(eg.rawEt());
0379   tmpEG.setIsoEt(eg.isoEt());
0380   tmpEG.setFootprintEt(eg.footprintEt());
0381   tmpEG.setNTT(eg.nTT());
0382   tmpEG.setShape(eg.shape());
0383   tmpEG.setTowerHoE(eg.towerHoE());
0384 
0385   return tmpEG;
0386 }
0387 
0388 l1t::Tau l1t::CaloTools::tauP4MP(l1t::Tau& tau) {
0389   l1t::Tau tmpTau(p4MP(&tau), tau.hwPt(), tau.hwEta(), tau.hwPhi(), tau.hwQual(), tau.hwIso());
0390   tmpTau.setTowerIPhi(tau.towerIPhi());
0391   tmpTau.setTowerIEta(tau.towerIEta());
0392   tmpTau.setRawEt(tau.rawEt());
0393   tmpTau.setIsoEt(tau.isoEt());
0394   tmpTau.setNTT(tau.nTT());
0395   tmpTau.setHasEM(tau.hasEM());
0396   tmpTau.setIsMerged(tau.isMerged());
0397 
0398   return tmpTau;
0399 }
0400 
0401 l1t::Jet l1t::CaloTools::jetP4MP(l1t::Jet& jet) {
0402   l1t::Jet tmpJet(p4MP(&jet), jet.hwPt(), jet.hwEta(), jet.hwPhi(), jet.hwQual());
0403   tmpJet.setTowerIPhi(jet.towerIPhi());
0404   tmpJet.setTowerIEta(jet.towerIEta());
0405   tmpJet.setRawEt(jet.rawEt());
0406   tmpJet.setSeedEt(jet.seedEt());
0407   tmpJet.setPUEt(jet.puEt());
0408   tmpJet.setPUDonutEt(0, jet.puDonutEt(0));
0409   tmpJet.setPUDonutEt(1, jet.puDonutEt(1));
0410   tmpJet.setPUDonutEt(2, jet.puDonutEt(2));
0411   tmpJet.setPUDonutEt(3, jet.puDonutEt(3));
0412 
0413   return tmpJet;
0414 }
0415 
0416 l1t::EtSum l1t::CaloTools::etSumP4MP(l1t::EtSum& etsum) {
0417   return l1t::EtSum(p4MP(&etsum), etsum.getType(), etsum.hwPt(), etsum.hwEta(), etsum.hwPhi(), etsum.hwQual());
0418 }
0419 unsigned int l1t::CaloTools::gloriousDivision(uint32_t aNumerator, uint32_t aDenominator) {
0420   static const uint64_t lLut[] = {
0421       0,     16777215, 4194304, 1864135, 1048576, 671089, 466034, 342392, 262144, 207126, 167772, 138655, 116508, 99273,
0422       85598, 74565,    65536,   58053,   51782,   46474,  41943,  38044,  34664,  31715,  29127,  26844,  24818,  23014,
0423       21400, 19949,    18641,   17458,   16384,   15406,  14513,  13696,  12945,  12255,  11619,  11030,  10486,  9980,
0424       9511,  9074,     8666,    8285,    7929,    7595,   7282,   6988,   6711,   6450,   6205,   5973,   5754,   5546,
0425       5350,  5164,     4987,    4820,    4660,    4509,   4365,   4227,   4096,   3971,   3852,   3737,   3628,   3524,
0426       3424,  3328,     3236,    3148,    3064,    2983,   2905,   2830,   2758,   2688,   2621,   2557,   2495,   2435,
0427       2378,  2322,     2268,    2217,    2166,    2118,   2071,   2026,   1982,   1940,   1899,   1859,   1820,   1783,
0428       1747,  1712,     1678,    1645,    1613,    1581,   1551,   1522,   1493,   1465,   1438,   1412,   1387,   1362,
0429       1337,  1314,     1291,    1269,    1247,    1226,   1205,   1185,   1165,   1146,   1127,   1109,   1091,   1074,
0430       1057,  1040,     1024,    1008,    993,     978,    963,    948,    934,    921,    907,    894,    881,    868,
0431       856,   844,      832,     820,     809,     798,    787,    776,    766,    756,    746,    736,    726,    717,
0432       707,   698,      689,     681,     672,     664,    655,    647,    639,    631,    624,    616,    609,    602,
0433       594,   587,      581,     574,     567,     561,    554,    548,    542,    536,    530,    524,    518,    512,
0434       506,   501,      496,     490,     485,     480,    475,    470,    465,    460,    455,    450,    446,    441,
0435       437,   432,      428,     424,     419,     415,    411,    407,    403,    399,    395,    392,    388,    384,
0436       380,   377,      373,     370,     366,     363,    360,    356,    353,    350,    347,    344,    340,    337,
0437       334,   331,      328,     326,     323,     320,    317,    314,    312,    309,    306,    304,    301,    299,
0438       296,   294,      291,     289,     286,     284,    282,    280,    277,    275,    273,    271,    268,    266,
0439       264,   262,      260,     258,     256,     254,    252,    250,    248,    246,    244,    243,    241,    239,
0440       237,   235,      234,     232,     230,     228,    227,    225,    223,    222,    220,    219,    217,    216,
0441       214,   212,      211,     209,     208,     207,    205,    204,    202,    201,    199,    198,    197,    195,
0442       194,   193,      191,     190,     189,     188,    186,    185,    184,    183,    182,    180,    179,    178,
0443       177,   176,      175,     173,     172,     171,    170,    169,    168,    167,    166,    165,    164,    163,
0444       162,   161,      160,     159,     158,     157,    156,    155,    154,    153,    152,    151,    150,    149,
0445       149,   148,      147,     146,     145,     144,    143,    143,    142,    141,    140,    139,    139,    138,
0446       137,   136,      135,     135,     134,     133,    132,    132,    131,    130,    129,    129,    128,    127,
0447       127,   126,      125,     125,     124,     123,    123,    122,    121,    121,    120,    119,    119,    118,
0448       117,   117,      116,     116,     115,     114,    114,    113,    113,    112,    111,    111,    110,    110,
0449       109,   109,      108,     108,     107,     106,    106,    105,    105,    104,    104,    103,    103,    102,
0450       102,   101,      101,     100,     100,     99,     99,     98,     98,     97,     97,     96,     96,     96,
0451       95,    95,       94,      94,      93,      93,     92,     92,     92,     91,     91,     90,     90,     89,
0452       89,    89,       88,      88,      87,      87,     87,     86,     86,     85,     85,     85,     84,     84,
0453       84,    83,       83,      82,      82,      82,     81,     81,     81,     80,     80,     80,     79,     79,
0454       79,    78,       78,      78,      77,      77,     77,     76,     76,     76,     75,     75,     75,     74,
0455       74,    74,       73,      73,      73,      73,     72,     72,     72,     71,     71,     71,     70,     70,
0456       70,    70,       69,      69,      69,      68,     68,     68,     68,     67,     67,     67,     67,     66,
0457       66,    66,       66,      65,      65,      65,     65,     64};
0458 
0459   // Firmware uses 18bit integers - make sure we are in the same range
0460   aNumerator &= 0x3FFFF;
0461   aDenominator &= 0x3FFFF;
0462 
0463   // Shift the denominator to optimise the polynomial expansion
0464   // I limit the shift to half the denominator size in the firmware to save on resources
0465   uint32_t lBitShift(0);
0466   for (; lBitShift != 9; ++lBitShift) {
0467     if (aDenominator & 0x20000)
0468       break;  // There is a 1 in the MSB
0469     aDenominator <<= 1;
0470   }
0471 
0472   // The magical polynomial expansion Voodoo
0473   uint64_t lInverseDenominator(((aDenominator & 0x3FE00) - (aDenominator & 0x001FF)) * (lLut[aDenominator >> 9]));
0474 
0475   // Save on DSPs by throwing away a bunch of LSBs
0476   lInverseDenominator >>= 17;
0477 
0478   // Multiply the numerator by the inverse denominator
0479   uint64_t lResult(aNumerator * lInverseDenominator);
0480 
0481   // Add two bits to the result, to make the Voodoo below work (saves us having an if-else on the shift direction)
0482   lResult <<= 2;
0483 
0484   // Restore the scale by taking into account the bitshift applied above.
0485   // We are now 18 bit left-shifted, so the 18 LSBs are effectively the fractional part...
0486 
0487   uint32_t aFractional = (lResult >>= (9 - lBitShift)) & 0x3FFFF;
0488   // ...and the top 18 bits are the integer part
0489 
0490   // uint32_t aInteger    = ( lResult >>= 18 ) & 0x3FFFF;
0491 
0492   unsigned int result = aFractional >> 10;
0493 
0494   return result;
0495 
0496   // Simples!
0497 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!