Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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