Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include <iostream>
0002 #include <algorithm>
0003 #include <strstream>
0004 
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ParameterSet/interface/FileInPath.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 
0010 #include "CondFormats/L1TObjects/interface/L1TMuonOverlapParams.h"
0011 
0012 #include "L1Trigger/L1TMuonOverlap/interface/OMTFProcessor.h"
0013 #include "L1Trigger/L1TMuonOverlap/interface/GoldenPattern.h"
0014 #include "L1Trigger/L1TMuonOverlap/interface/OMTFinput.h"
0015 #include "L1Trigger/L1TMuonOverlap/interface/OMTFResult.h"
0016 
0017 #include "L1Trigger/RPCTrigger/interface/RPCConst.h"
0018 
0019 #include "SimDataFormats/Track/interface/SimTrack.h"
0020 ///////////////////////////////////////////////
0021 ///////////////////////////////////////////////
0022 OMTFProcessor::~OMTFProcessor() {
0023   for (auto it : theGPs)
0024     delete it.second;
0025 }
0026 ///////////////////////////////////////////////
0027 ///////////////////////////////////////////////
0028 void OMTFProcessor::resetConfiguration() {
0029   myResults.clear();
0030   for (auto it : theGPs)
0031     delete it.second;
0032   theGPs.clear();
0033 }
0034 ///////////////////////////////////////////////
0035 ///////////////////////////////////////////////
0036 bool OMTFProcessor::configure(const OMTFConfiguration *omtfConfig, const L1TMuonOverlapParams *omtfPatterns) {
0037   resetConfiguration();
0038 
0039   myOmtfConfig = omtfConfig;
0040 
0041   myResults.assign(myOmtfConfig->nTestRefHits(), OMTFProcessor::resultsMap());
0042 
0043   const l1t::LUT *chargeLUT = omtfPatterns->chargeLUT();
0044   const l1t::LUT *etaLUT = omtfPatterns->etaLUT();
0045   const l1t::LUT *ptLUT = omtfPatterns->ptLUT();
0046   const l1t::LUT *pdfLUT = omtfPatterns->pdfLUT();
0047   const l1t::LUT *meanDistPhiLUT = omtfPatterns->meanDistPhiLUT();
0048 
0049   unsigned int nGPs = myOmtfConfig->nGoldenPatterns();
0050   unsigned int address = 0;
0051   unsigned int iEta, iPt;
0052   int iCharge;
0053   for (unsigned int iGP = 0; iGP < nGPs; ++iGP) {
0054     address = iGP;
0055     iEta = etaLUT->data(address);
0056     iCharge = chargeLUT->data(address) == 0 ? -1 : 1;
0057     iPt = ptLUT->data(address);
0058 
0059     GoldenPattern::vector2D meanDistPhi2D(myOmtfConfig->nLayers());
0060     GoldenPattern::vector1D pdf1D(exp2(myOmtfConfig->nPdfAddrBits()));
0061     GoldenPattern::vector3D pdf3D(myOmtfConfig->nLayers());
0062     GoldenPattern::vector2D pdf2D(myOmtfConfig->nRefLayers());
0063     ///Mean dist phi data
0064     for (unsigned int iLayer = 0; iLayer < myOmtfConfig->nLayers(); ++iLayer) {
0065       GoldenPattern::vector1D meanDistPhi1D(myOmtfConfig->nRefLayers());
0066       for (unsigned int iRefLayer = 0; iRefLayer < myOmtfConfig->nRefLayers(); ++iRefLayer) {
0067         address = iRefLayer + iLayer * myOmtfConfig->nRefLayers() +
0068                   iGP * (myOmtfConfig->nRefLayers() * myOmtfConfig->nLayers());
0069         meanDistPhi1D[iRefLayer] = meanDistPhiLUT->data(address) - (1 << (meanDistPhiLUT->nrBitsData() - 1));
0070       }
0071       meanDistPhi2D[iLayer] = meanDistPhi1D;
0072       ///Pdf data
0073       for (unsigned int iRefLayer = 0; iRefLayer < myOmtfConfig->nRefLayers(); ++iRefLayer) {
0074         pdf1D.assign(1 << myOmtfConfig->nPdfAddrBits(), 0);
0075         for (unsigned int iPdf = 0; iPdf < (unsigned int)(1 << myOmtfConfig->nPdfAddrBits()); ++iPdf) {
0076           address = iPdf + iRefLayer * (1 << myOmtfConfig->nPdfAddrBits()) +
0077                     iLayer * myOmtfConfig->nRefLayers() * (1 << myOmtfConfig->nPdfAddrBits()) +
0078                     iGP * myOmtfConfig->nLayers() * myOmtfConfig->nRefLayers() * (1 << myOmtfConfig->nPdfAddrBits());
0079           pdf1D[iPdf] = pdfLUT->data(address);
0080         }
0081         pdf2D[iRefLayer] = pdf1D;
0082       }
0083       pdf3D[iLayer] = pdf2D;
0084     }
0085     Key aKey(iEta, iPt, iCharge, iGP);
0086 
0087     GoldenPattern *aGP = new GoldenPattern(aKey, myOmtfConfig);
0088     aGP->setMeanDistPhi(meanDistPhi2D);
0089     aGP->setPdf(pdf3D);
0090     addGP(aGP);
0091   }
0092   return true;
0093 }
0094 ///////////////////////////////////////////////
0095 ///////////////////////////////////////////////
0096 bool OMTFProcessor::addGP(GoldenPattern *aGP) {
0097   if (theGPs.find(aGP->key()) != theGPs.end()) {
0098     throw cms::Exception("Corrupted Golden Patterns data")
0099         << "OMTFProcessor::addGP(...) "
0100         << " Reading two Golden Patterns with the same key: " << aGP->key() << std::endl;
0101   } else
0102     theGPs[aGP->key()] = aGP;
0103 
0104   for (auto &itRegion : myResults) {
0105     OMTFResult aResult;
0106     aResult.configure(myOmtfConfig);
0107     itRegion[aGP->key()] = aResult;
0108   }
0109 
0110   return true;
0111 }
0112 ///////////////////////////////////////////////
0113 ///////////////////////////////////////////////
0114 void OMTFProcessor::averagePatterns(int charge) {
0115   Key aKey(0, 9, charge);
0116 
0117   while (theGPs.find(aKey) != theGPs.end()) {
0118     GoldenPattern *aGP1 = theGPs.find(aKey)->second;
0119     GoldenPattern *aGP2 = aGP1;
0120     GoldenPattern *aGP3 = aGP1;
0121     GoldenPattern *aGP4 = aGP1;
0122 
0123     ++aKey.thePtCode;
0124     while (theGPs.find(aKey) == theGPs.end() && aKey.thePtCode <= 401)
0125       ++aKey.thePtCode;
0126     if (aKey.thePtCode <= 401 && theGPs.find(aKey) != theGPs.end())
0127       aGP2 = theGPs.find(aKey)->second;
0128 
0129     if (aKey.thePtCode > 71) {
0130       ++aKey.thePtCode;
0131       while (theGPs.find(aKey) == theGPs.end() && aKey.thePtCode <= 401)
0132         ++aKey.thePtCode;
0133       if (aKey.thePtCode <= 401 && theGPs.find(aKey) != theGPs.end())
0134         aGP3 = theGPs.find(aKey)->second;
0135 
0136       ++aKey.thePtCode;
0137       while (theGPs.find(aKey) == theGPs.end() && aKey.thePtCode <= 401)
0138         ++aKey.thePtCode;
0139       if (aKey.thePtCode <= 401 && theGPs.find(aKey) != theGPs.end())
0140         aGP4 = theGPs.find(aKey)->second;
0141     } else {
0142       aGP3 = aGP1;
0143       aGP4 = aGP2;
0144     }
0145     //HACK. Have to clean this up.
0146     ///Previously pt codes were going by steps of 1, now this is not the case
0147     ++aKey.thePtCode;
0148     while (theGPs.find(aKey) == theGPs.end() && aKey.thePtCode <= 401)
0149       ++aKey.thePtCode;
0150     ///////////////////////////////
0151 
0152     GoldenPattern::vector2D meanDistPhi = aGP1->getMeanDistPhi();
0153 
0154     GoldenPattern::vector2D meanDistPhi1 = aGP1->getMeanDistPhi();
0155     GoldenPattern::vector2D meanDistPhi2 = aGP2->getMeanDistPhi();
0156     GoldenPattern::vector2D meanDistPhi3 = aGP3->getMeanDistPhi();
0157     GoldenPattern::vector2D meanDistPhi4 = aGP4->getMeanDistPhi();
0158 
0159     for (unsigned int iLayer = 0; iLayer < myOmtfConfig->nLayers(); ++iLayer) {
0160       for (unsigned int iRefLayer = 0; iRefLayer < myOmtfConfig->nRefLayers(); ++iRefLayer) {
0161         meanDistPhi[iLayer][iRefLayer] += meanDistPhi2[iLayer][iRefLayer];
0162         meanDistPhi[iLayer][iRefLayer] += meanDistPhi3[iLayer][iRefLayer];
0163         meanDistPhi[iLayer][iRefLayer] += meanDistPhi4[iLayer][iRefLayer];
0164         meanDistPhi[iLayer][iRefLayer] /= 4;
0165       }
0166     }
0167 
0168     aGP1->setMeanDistPhi(meanDistPhi);
0169     aGP2->setMeanDistPhi(meanDistPhi);
0170 
0171     shiftGP(aGP1, meanDistPhi, meanDistPhi1);
0172     shiftGP(aGP2, meanDistPhi, meanDistPhi2);
0173     if (aGP3 != aGP1 && aGP4 != aGP2) {
0174       aGP3->setMeanDistPhi(meanDistPhi);
0175       aGP4->setMeanDistPhi(meanDistPhi);
0176       shiftGP(aGP3, meanDistPhi, meanDistPhi3);
0177       shiftGP(aGP4, meanDistPhi, meanDistPhi4);
0178     }
0179   }
0180 }
0181 ///////////////////////////////////////////////
0182 ///////////////////////////////////////////////
0183 void OMTFProcessor::shiftGP(GoldenPattern *aGP,
0184                             const GoldenPattern::vector2D &meanDistPhiNew,
0185                             const GoldenPattern::vector2D &meanDistPhiOld) {
0186   ///Shift pdfs by differecne between original menaDistPhi, and
0187   ///the averaged value
0188   unsigned int nPdfBins = exp2(myOmtfConfig->nPdfAddrBits());
0189   GoldenPattern::vector3D pdfAllRef = aGP->getPdf();
0190 
0191   int indexShift = 0;
0192   for (unsigned int iLayer = 0; iLayer < myOmtfConfig->nLayers(); ++iLayer) {
0193     for (unsigned int iRefLayer = 0; iRefLayer < myOmtfConfig->nRefLayers(); ++iRefLayer) {
0194       indexShift = meanDistPhiOld[iLayer][iRefLayer] - meanDistPhiNew[iLayer][iRefLayer];
0195       for (unsigned int iPdfBin = 0; iPdfBin < nPdfBins; ++iPdfBin)
0196         pdfAllRef[iLayer][iRefLayer][iPdfBin] = 0;
0197       for (unsigned int iPdfBin = 0; iPdfBin < nPdfBins; ++iPdfBin) {
0198         if ((int)(iPdfBin) + indexShift >= 0 && iPdfBin + indexShift < nPdfBins)
0199           pdfAllRef[iLayer][iRefLayer][iPdfBin + indexShift] = aGP->pdfValue(iLayer, iRefLayer, iPdfBin);
0200       }
0201     }
0202   }
0203   aGP->setPdf(pdfAllRef);
0204 }
0205 ///////////////////////////////////////////////
0206 ///////////////////////////////////////////////
0207 const std::map<Key, GoldenPattern *> &OMTFProcessor::getPatterns() const { return theGPs; }
0208 ///////////////////////////////////////////////
0209 ///////////////////////////////////////////////
0210 const std::vector<OMTFProcessor::resultsMap> &OMTFProcessor::processInput(unsigned int iProcessor,
0211                                                                           const OMTFinput &aInput) {
0212   for (auto &itRegion : myResults)
0213     for (auto &itKey : itRegion)
0214       itKey.second.clear();
0215 
0216   //////////////////////////////////////
0217   //////////////////////////////////////
0218   std::bitset<128> refHitsBits = aInput.getRefHits(iProcessor);
0219   if (refHitsBits.none())
0220     return myResults;
0221 
0222   for (unsigned int iLayer = 0; iLayer < myOmtfConfig->nLayers(); ++iLayer) {
0223     const OMTFinput::vector1D &layerHits = aInput.getLayerData(iLayer);
0224     if (layerHits.empty())
0225       continue;
0226     ///Number of reference hits to be checked.
0227     unsigned int nTestedRefHits = myOmtfConfig->nTestRefHits();
0228     for (unsigned int iRefHit = 0; iRefHit < myOmtfConfig->nRefHits(); ++iRefHit) {
0229       if (!refHitsBits[iRefHit])
0230         continue;
0231       if (nTestedRefHits-- == 0)
0232         break;
0233       const RefHitDef &aRefHitDef = myOmtfConfig->getRefHitsDefs()[iProcessor][iRefHit];
0234 
0235       int phiRef = aInput.getLayerData(myOmtfConfig->getRefToLogicNumber()[aRefHitDef.iRefLayer])[aRefHitDef.iInput];
0236       int etaRef =
0237           aInput.getLayerData(myOmtfConfig->getRefToLogicNumber()[aRefHitDef.iRefLayer], true)[aRefHitDef.iInput];
0238       unsigned int iRegion = aRefHitDef.iRegion;
0239 
0240       if (myOmtfConfig->getBendingLayers().count(iLayer))
0241         phiRef = 0;
0242       const OMTFinput::vector1D restrictedLayerHits = restrictInput(iProcessor, iRegion, iLayer, layerHits);
0243       for (auto itGP : theGPs) {
0244         GoldenPattern::layerResult aLayerResult =
0245             itGP.second->process1Layer1RefLayer(aRefHitDef.iRefLayer, iLayer, phiRef, restrictedLayerHits);
0246         int phiRefSt2 = itGP.second->propagateRefPhi(phiRef, etaRef, aRefHitDef.iRefLayer);
0247         myResults[myOmtfConfig->nTestRefHits() - nTestedRefHits - 1][itGP.second->key()].setRefPhiRHits(
0248             aRefHitDef.iRefLayer, phiRef);
0249         myResults[myOmtfConfig->nTestRefHits() - nTestedRefHits - 1][itGP.second->key()].addResult(
0250             aRefHitDef.iRefLayer, iLayer, aLayerResult.first, phiRefSt2, etaRef);
0251       }
0252     }
0253   }
0254   //////////////////////////////////////
0255   //////////////////////////////////////
0256   for (auto &itRefHit : myResults)
0257     for (auto &itKey : itRefHit)
0258       itKey.second.finalise();
0259 
0260   std::ostringstream myStr;
0261   myStr << "iProcessor: " << iProcessor << std::endl;
0262   myStr << "Input: ------------" << std::endl;
0263   myStr << aInput << std::endl;
0264   edm::LogInfo("OMTF processor") << myStr.str();
0265 
0266   return myResults;
0267 }
0268 ////////////////////////////////////////////
0269 ////////////////////////////////////////////
0270 OMTFinput::vector1D OMTFProcessor::restrictInput(unsigned int iProcessor,
0271                                                  unsigned int iRegion,
0272                                                  unsigned int iLayer,
0273                                                  const OMTFinput::vector1D &layerHits) {
0274   OMTFinput::vector1D myHits = layerHits;
0275 
0276   unsigned int iStart = myOmtfConfig->getConnections()[iProcessor][iRegion][iLayer].first;
0277   unsigned int iEnd = iStart + myOmtfConfig->getConnections()[iProcessor][iRegion][iLayer].second - 1;
0278 
0279   for (unsigned int iInput = 0; iInput < 14; ++iInput) {
0280     if (iInput < iStart || iInput > iEnd)
0281       myHits[iInput] = myOmtfConfig->nPhiBins();
0282   }
0283   return myHits;
0284 }
0285 ////////////////////////////////////////////
0286 ////////////////////////////////////////////
0287 void OMTFProcessor::fillCounts(unsigned int iProcessor, const OMTFinput &aInput, const SimTrack *aSimMuon) {
0288   int theCharge = (abs(aSimMuon->type()) == 13) ? aSimMuon->type() / -13 : 0;
0289   unsigned int iPt = RPCConst::iptFromPt(aSimMuon->momentum().pt());
0290   ///Stupid conersion. Have to go through PAC pt scale, as we later
0291   ///shift resulting pt code by +1
0292   iPt += 1;
0293   if (iPt > 31)
0294     iPt = 200 * 2 + 1;
0295   else
0296     iPt = RPCConst::ptFromIpt(iPt) * 2.0 +
0297           1;  //MicroGMT has 0.5 GeV step size, with lower bin edge  (uGMT_pt_code - 1)*step_size
0298   //////
0299 
0300   //////////////////////////////////////
0301   std::bitset<128> refHitsBits = aInput.getRefHits(iProcessor);
0302   if (refHitsBits.none())
0303     return;
0304 
0305   std::ostringstream myStr;
0306   myStr << "iProcessor: " << iProcessor << std::endl;
0307   myStr << "Input: ------------" << std::endl;
0308   myStr << aInput << std::endl;
0309   edm::LogInfo("OMTF processor") << myStr.str();
0310 
0311   for (unsigned int iLayer = 0; iLayer < myOmtfConfig->nLayers(); ++iLayer) {
0312     const OMTFinput::vector1D &layerHits = aInput.getLayerData(iLayer);
0313     if (layerHits.empty())
0314       continue;
0315     ///Number of reference hits to be checked.
0316     for (unsigned int iRefHit = 0; iRefHit < myOmtfConfig->nRefHits(); ++iRefHit) {
0317       if (!refHitsBits[iRefHit])
0318         continue;
0319       const RefHitDef &aRefHitDef = myOmtfConfig->getRefHitsDefs()[iProcessor][iRefHit];
0320       int phiRef = aInput.getLayerData(myOmtfConfig->getRefToLogicNumber()[aRefHitDef.iRefLayer])[aRefHitDef.iInput];
0321       unsigned int iRegion = aRefHitDef.iRegion;
0322       if (myOmtfConfig->getBendingLayers().count(iLayer))
0323         phiRef = 0;
0324       const OMTFinput::vector1D restrictedLayerHits = restrictInput(iProcessor, iRegion, iLayer, layerHits);
0325       for (auto itGP : theGPs) {
0326         if (itGP.first.theCharge != theCharge)
0327           continue;
0328         if (itGP.first.thePtCode != iPt)
0329           continue;
0330         itGP.second->addCount(aRefHitDef.iRefLayer, iLayer, phiRef, restrictedLayerHits);
0331       }
0332     }
0333   }
0334 }
0335 ////////////////////////////////////////////
0336 ////////////////////////////////////////////