src/Tools/PatternGenerator.cc

0001 /*
0002  * PatternGenerator.cc
0003  *
0004  *  Created on: Nov 8, 2019
0005  *      Author: kbunkow
0006  */
0007
0008 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Tools/PatternGenerator.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <boost/range/adaptor/reversed.hpp>
0012
0013 #include "TFile.h"
0014 #include "TDirectory.h"
0015
0016 PatternGenerator::PatternGenerator(const edm::ParameterSet& edmCfg,
0017                                    const OMTFConfiguration* omtfConfig,
0018                                    GoldenPatternVec<GoldenPatternWithStat>& gps)
0019     : PatternOptimizerBase(edmCfg, omtfConfig, gps), eventCntPerGp(gps.size(), 0) {
0020   edm::LogImportant("l1tOmtfEventPrint") << "constructing PatternGenerator, type: "
0021                                          << edmCfg.getParameter<string>("patternGenerator") << std::endl;
0022
0023   if (edmCfg.getParameter<string>("patternGenerator") == "patternGen")
0024     initPatternGen();
0025 }
0026
0027 PatternGenerator::~PatternGenerator() {}
0028
0029 void PatternGenerator::initPatternGen() {
0030   //reseting the golden patterns
0031   unsigned int i = 0;
0032   for (auto& gp : goldenPatterns) {
0033     gp->setKeyNumber(i++);  //needed  if patterns were added
0034
0035     if (gp->key().thePt == 0)
0036       continue;
0037
0038     gp->reset();
0039
0040     int statBinsCnt = 1024;  //gp->getPdf()[0][0].size() * 8; //TODO should be big enough to comprise the pdf tails
0041     gp->iniStatisitics(statBinsCnt, 1);  //TODO
0042   }
0043
0044   edm::LogImportant("l1tOmtfEventPrint") << "PatternGenerator::initPatternGen():" << __LINE__
0045                                          << " goldenPatterns.size() " << goldenPatterns.size() << std::endl;
0046
0047   //GoldenPatternResult::setFinalizeFunction(3); TODO why it was this one????
0048   // edm::LogImportant("l1tOmtfEventPrint") << "reseting golden pattern !!!!!" << std::endl;
0049
0050   //setting all pdf to 1, this will cause that  when the OmtfProcessor process the input, the result will be based only on the number of fired layers,
0051   //and then the omtfCand will come from the processor that has the biggest number of fired layers
0052   //however, if the GoldenPatternResult::finalise3() is used - which just count the number of muonStubs (but do not check if it is valid, i.e. fired the pdf)
0053   // - the below does not matter
0054   for (auto& gp : goldenPatterns) {
0055     for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
0056       for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
0057         for (unsigned int iBin = 0; iBin < gp->getPdf()[iLayer][iRefLayer].size(); iBin++) {
0058           gp->pdfAllRef[iLayer][iRefLayer][iBin] = 1;
0059         }
0060       }
0061     }
0062   }
0063
0064   //TODO uncomment if filling the ptDeltaPhiHist is needed
0065   /*  ptDeltaPhiHists.resize(2);
0066   for(unsigned int iCharge = 0; iCharge <= 1; iCharge++) {
0067     for(unsigned int iLayer = 0; iLayer < omtfConfig->nLayers(); ++iLayer) { //for the moment filing only ref layer, remove whe
0068       if(iLayer == 0 || iLayer == 2 || iLayer == 4 || iLayer == 6 || iLayer == 7 || iLayer == 10 || iLayer == 11  || iLayer == 16 || //refLayars
0069          iLayer == 1 || iLayer == 3 || iLayer == 5  ) //bending layers
0070       {
0071         ostringstream name;
0072         name<<"ptDeltaPhiHist_ch_"<<iCharge<<"_Layer_"<<iLayer;
0073         int phiFrom = -10;
0074         int phiTo   = 300; //TODO
0075         int phiBins = phiTo - phiFrom;
0076
0077         if(iCharge == 1) {
0078           phiFrom = -300; //TODO
0079           phiTo = 10;
0080         }
0081
0082         TH2I* ptDeltaPhiHist = new TH2I(name.str().c_str(), name.str().c_str(), 400, 0, 200, phiBins, phiFrom -0.5, phiTo -0.5);
0083         //cout<<"BinLowEdge "<<ptDeltaPhiHist->GetYaxis()->GetBinLowEdge(100)<<" BinUpEdge "<<ptDeltaPhiHist->GetYaxis()->GetBinUpEdge(100);
0084         ptDeltaPhiHists[iCharge].push_back(ptDeltaPhiHist);
0085       }
0086       else
0087         ptDeltaPhiHists[iCharge].push_back(nullptr);
0088     }
0089   }*/
0090
0091   /* cannot be called  here, will cause crash
0092   edm::LogImportant("OMTFReconstruction")<<" PatternGenerator constructor - patterns after modification "<<std::endl;
0093   for(auto& gp : goldenPatterns) {
0094     edm::LogImportant("OMTFReconstruction")<<gp->key()<<" "
0095         <<omtfConfig->getPatternPtRange(gp->key().theNumber).ptFrom
0096         <<" - "<<omtfConfig->getPatternPtRange(gp->key().theNumber).ptTo<<" GeV"<<std::endl;
0097   }*/
0098 }
0099
0100 void PatternGenerator::updateStat() {
0101   //cout<<__FUNCTION__<<":"<<__LINE__<<" omtfCand "<<*omtfCand<<std::endl;;
0102   AlgoMuon* algoMuon = omtfCand.get();
0103   if (!algoMuon) {
0104     edm::LogImportant("l1tOmtfEventPrint") << ":" << __LINE__ << " algoMuon is null" << std::endl;
0105     throw runtime_error("algoMuon is null");
0106   }
0107
0108   double ptSim = simMuon->momentum().pt();
0109   int chargeSim = (abs(simMuon->type()) == 13) ? simMuon->type() / -13 : 0;
0110
0111   unsigned int exptPatNum = omtfConfig->getPatternNum(ptSim, chargeSim);
0112   GoldenPatternWithStat* exptCandGp = goldenPatterns.at(exptPatNum).get();  // expected pattern
0113
0114   eventCntPerGp[exptPatNum]++;
0115
0116   //edm::LogImportant("l1tOmtfEventPrint")<<"\n" <<__FUNCTION__<<": "<<__LINE__<<" exptCandGp "<<exptCandGp->key()<<" candProcIndx "<<candProcIndx<<" ptSim "<<ptSim<<" chargeSim "<<chargeSim<<std::endl;
0117
0118   int pdfMiddle = 1 << (omtfConfig->nPdfAddrBits() - 1);
0119
0120   //iRefHit is the index of the hit
0121   for (unsigned int iRefHit = 0; iRefHit < exptCandGp->getResults()[candProcIndx].size(); ++iRefHit) {
0122     auto& gpResult = exptCandGp->getResults()[candProcIndx][iRefHit];
0123
0124     unsigned int refLayer = gpResult.getRefLayer();
0125
0126     if (gpResult.getFiredLayerCnt() >= 3) {
0127       for (unsigned int iLayer = 0; iLayer < gpResult.getStubResults().size(); iLayer++) {
0128         //updating statistic for the gp which should have fired
0129
0130         bool fired = false;
0131         if (gpResult.getStubResults()[iLayer].getMuonStub()) {
0132           if (omtfConfig->isBendingLayer(iLayer)) {
0133             if (gpResult.getStubResults()[iLayer].getMuonStub()->qualityHw >= 4)  //TODO change quality cut if needed
0134               fired = true;
0135           } else
0136             fired = true;
0137         }
0138
0139         if (fired) {  //the result is not empty
0140           int phiDist = gpResult.getStubResults()[iLayer].getPdfBin();
0141           phiDist += exptCandGp->meanDistPhiValue(iLayer, refLayer) - pdfMiddle;
0142           //removing the shift applied in the GoldenPatternBase::process1Layer1RefLayer
0143
0144           //TODO uncomment if filling ptDeltaPhiHists is needed
0145           /*
0146           unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefHit];
0147           if(ptDeltaPhiHists[iCharge][iLayer] != nullptr &&
0148               (iLayer == refLayerLogicNum || omtfConfig->getLogicToLogic().at(iLayer) == (int)refLayerLogicNum) )
0149             ptDeltaPhiHists[iCharge][iLayer]->Fill(ttAlgoMuon->getPt(), phiDist); //TODO correct
0150            */
0151
0152           phiDist += exptCandGp->getStatistics()[iLayer][refLayer].size() / 2;
0153
0154           //edm::LogImportant("l1tOmtfEventPrint")<<__FUNCTION__<<":"<<__LINE__<<" refLayer "<<refLayer<<" iLayer "<<iLayer<<" phiDist "<<phiDist<<" getPdfBin "<<gpResult.getStubResults()[iLayer].getPdfBin()<<std::endl;
0155           if (phiDist > 0 && phiDist < (int)(exptCandGp->getStatistics()[iLayer][refLayer].size())) {
0156             //updating statistic for the gp which found the candidate
0157             //edm::LogImportant("l1tOmtfEventPrint")<<__FUNCTION__<<":"<<__LINE__<<" updating statistic "<<std::endl;
0158             exptCandGp->updateStat(iLayer, refLayer, phiDist, 0, 1);
0159           }
0160         } else {  //if there is no hit at all in a given layer, the bin = 0 is filled
0161           int phiDist = 0;
0162           exptCandGp->updateStat(iLayer, refLayer, phiDist, 0, 1);
0163         }
0164       }
0165     }
0166   }
0167 }
0168
0169 void PatternGenerator::observeEventEnd(const edm::Event& iEvent,
0170                                        std::unique_ptr<l1t::RegionalMuonCandBxCollection>& finalCandidates) {
0171   if (simMuon == nullptr || omtfCand->getGoldenPatern() == nullptr)  //no sim muon or empty candidate
0172     return;
0173
0174   if (abs(simMuon->momentum().eta()) < 0.8 || abs(simMuon->momentum().eta()) > 1.24)
0175     return;
0176
0177   PatternOptimizerBase::observeEventEnd(iEvent, finalCandidates);
0178
0179   updateStat();
0180 }
0181
0182 void PatternGenerator::endJob() {
0183   if (edmCfg.getParameter<string>("patternGenerator") == "modifyClassProb")
0184     modifyClassProb(1);
0185   else if (edmCfg.getParameter<string>("patternGenerator") == "groupPatterns")
0186     groupPatterns();
0187   else if (edmCfg.getParameter<string>("patternGenerator") == "patternGen") {
0188     upadatePdfs();
0189     writeLayerStat = true;
0190   } else if (edmCfg.getParameter<string>("patternGenerator") == "patternGenFromStat") {
0191     std::string rootFileName = edmCfg.getParameter<edm::FileInPath>("patternsROOTFile").fullPath();
0192     edm::LogImportant("l1tOmtfEventPrint") << "PatternGenerator::endJob() rootFileName " << rootFileName << std::endl;
0193     TFile inFile(rootFileName.c_str());
0194     TDirectory* curDir = (TDirectory*)inFile.Get("layerStats");
0195
0196     ostringstream ostrName;
0197     for (auto& gp : goldenPatterns) {
0198       if (gp->key().thePt == 0)
0199         continue;
0200
0201       int statBinsCnt = 1024;  //= gp->getPdf()[0][0].size() * 8; //TODO should be big enough to comprise the pdf tails
0202       gp->iniStatisitics(statBinsCnt, 1);  //TODO
0203
0204       for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
0205         for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
0206           ostrName.str("");
0207           ostrName << "histLayerStat_PatNum_" << gp->key().theNumber << "_refLayer_" << iRefLayer << "_Layer_"
0208                    << iLayer;
0209
0210           TH1I* histLayerStat = (TH1I*)curDir->Get(ostrName.str().c_str());
0211
0212           if (histLayerStat) {
0213             for (int iBin = 0; iBin < statBinsCnt; iBin++) {
0214               gp->updateStat(iLayer, iRefLayer, iBin, 0, histLayerStat->GetBinContent(iBin + 1));
0215             }
0216           } else {
0217             edm::LogImportant("l1tOmtfEventPrint")
0218                 << "PatternGenerator::endJob() - reading histLayerStat: histogram not found " << ostrName.str()
0219                 << std::endl;
0220           }
0221         }
0222       }
0223     }
0224
0225     TH1* simMuFoundByOmtfPt_fromFile = (TH1*)inFile.Get("simMuFoundByOmtfPt");
0226     for (unsigned int iGp = 0; iGp < eventCntPerGp.size(); iGp++) {
0227       eventCntPerGp[iGp] = simMuFoundByOmtfPt_fromFile->GetBinContent(simMuFoundByOmtfPt_fromFile->FindBin(iGp));
0228       edm::LogImportant("l1tOmtfEventPrint")
0229           << "PatternGenerator::endJob() - eventCntPerGp: iGp" << iGp << " - " << eventCntPerGp[iGp] << std::endl;
0230     }
0231
0232     //TODO chose the desired grouping ///////////////
0233     int group = 0;
0234     int indexInGroup = 0;
0235     for (auto& gp : goldenPatterns) {
0236       indexInGroup++;
0237       gp->key().setGroup(group);
0238       gp->key().setIndexInGroup(indexInGroup);
0239       //indexInGroup is counted from 1
0240
0241       edm::LogImportant("l1tOmtfEventPrint")
0242           << "setGroup(group): group " << group << " indexInGroup " << indexInGroup << std::endl;
0243
0244       if (gp->key().thePt <= 10 && indexInGroup == 2) {  //TODO
0245         indexInGroup = 0;
0246         group++;
0247       }
0248
0249       if (gp->key().thePt > 10 && indexInGroup == 4) {  //TODO
0250         indexInGroup = 0;
0251         group++;
0252       }
0253     }  /////////////////////////////////////////////
0254
0255     upadatePdfs();
0256
0257     modifyClassProb(1);
0258
0259     //groupPatterns(); IMPORTANT don't call grouping here, just set the groups above!!!!
0260
0261     reCalibratePt();
0262     this->writeLayerStat = true;
0263   }
0264
0265   PatternOptimizerBase::endJob();
0266 }
0267
0268 void PatternGenerator::upadatePdfs() {
0269   //TODO setting the DistPhiBitShift i.e. grouping of the pdfBins
0270   for (auto& gp : goldenPatterns) {
0271     if (gp->key().thePt == 0)
0272       continue;
0273     for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
0274       for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
0275         if (gp->getDistPhiBitShift(iLayer, iRefLayer)) {
0276           throw runtime_error(
0277               string(__FUNCTION__) + ":" + to_string(__LINE__) +
0278               "gp->getDistPhiBitShift(iLayer, iRefLayer) != 0 -  cannot change DistPhiBitShift then!!!!");
0279         }
0280
0281         if ((gp->key().thePt <= 10) && (iLayer == 1 || iLayer == 3 || iLayer == 5)) {
0282           gp->setDistPhiBitShift(1, iLayer, iRefLayer);
0283         } else
0284           gp->setDistPhiBitShift(0, iLayer, iRefLayer);
0285
0286         //watch out: the shift in a given layer must be the same for patterns in one group
0287         //todo  make the setting of the shift on the group base
0288         //TODO set the DistPhiBitShift
0289         /*if( (gp->key().thePt <= 10) && (iLayer == 3 || iLayer == 5 ) && (iRefLayer == 0 || iRefLayer == 2 || iRefLayer == 6 || iRefLayer == 7)) {
0290           gp->setDistPhiBitShift(3, iLayer, iRefLayer);
0291         }
0292         else if( (gp->key().thePt <= 10) && ( iLayer == 1 || iLayer == 3 || iLayer == 5 ) ) {
0293           gp->setDistPhiBitShift(2, iLayer, iRefLayer);
0294         }
0295         else if( ( (gp->key().thePt <= 10) && (iLayer == 7 ||iLayer == 8 || iLayer == 17 ) ) ) {
0296           gp->setDistPhiBitShift(1, iLayer, iRefLayer);
0297         }
0298         else if( (gp->key().thePt <= 10) && (iLayer == 10 || iLayer == 11 || iLayer == 12 || iLayer == 13) && (iRefLayer == 1)) {
0299           gp->setDistPhiBitShift(1, iLayer, iRefLayer);
0300         }*/
0301       }
0302     }
0303   }
0304
0305   double minHitCntThresh = 0.001;
0306   //Calculating meanDistPhi
0307   for (auto& gp : goldenPatterns) {
0308     if (gp->key().thePt == 0)
0309       continue;
0310
0311     int minHitCnt = minHitCntThresh * eventCntPerGp[gp->key().number()];  // //TODO tune threshold <<<<<<<<<<<<<<<<<<
0312     edm::LogImportant("l1tOmtfEventPrint")
0313         << "PatternGenerator::upadatePdfs() Calculating meanDistPhi " << gp->key() << " eventCnt "
0314         << eventCntPerGp[gp->key().number()] << " minHitCnt " << minHitCnt << std::endl;
0315
0316     for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
0317       for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
0318         //calculate meanDistPhi
0319         double meanDistPhi = 0;
0320         double count = 0;
0321         for (unsigned int iBin = 1; iBin < gp->getStatistics()[iLayer][iRefLayer].size(); iBin++) {
0322           //iBin = 0 is reserved for the no hit
0323           meanDistPhi += iBin * gp->getStatistics()[iLayer][iRefLayer][iBin][0];
0324           count += gp->getStatistics()[iLayer][iRefLayer][iBin][0];
0325         }
0326
0327         if (count != 0) {
0328           meanDistPhi /= count;
0329
0330           meanDistPhi -= (gp->getStatistics()[iLayer][iRefLayer].size() / 2);
0331
0332           if (count < minHitCnt)
0333             meanDistPhi = 0;
0334           else
0335             edm::LogImportant("l1tOmtfEventPrint")
0336                 << __FUNCTION__ << ": " << __LINE__ << " " << gp->key() << " iLayer " << iLayer << " iRefLayer "
0337                 << iRefLayer << " count " << count << " meanDistPhi " << meanDistPhi << endl;
0338         }
0339         gp->setMeanDistPhiValue(round(meanDistPhi), iLayer, iRefLayer);
0340       }
0341     }
0342   }
0343
0344   OMTFConfiguration::vector2D patternGroups = omtfConfig->getPatternGroups(goldenPatterns);
0345   edm::LogImportant("l1tOmtfEventPrint") << "patternGroups:" << std::endl;
0346   for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
0347     edm::LogImportant("l1tOmtfEventPrint") << "patternGroup " << std::setw(2) << iGroup << " ";
0348     for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
0349       edm::LogImportant("l1tOmtfEventPrint") << i << " patNum " << patternGroups[iGroup][i] << " ";
0350     }
0351     edm::LogImportant("l1tOmtfEventPrint") << std::endl;
0352   }
0353
0354   //averaging the meanDistPhi for the gp belonging to the same group
0355   for (unsigned int iLayer = 0; iLayer < goldenPatterns.at(0)->getPdf().size(); ++iLayer) {
0356     for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns.at(0)->getPdf()[iLayer].size(); ++iRefLayer) {
0357       //unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefLayer];
0358       //if(refLayerLogicNum == iLayer)
0359       {
0360         for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
0361           double meanDistPhi = 0;
0362           int mergedCnt = 0;
0363           for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
0364             auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
0365             meanDistPhi += gp->meanDistPhiValue(iLayer, iRefLayer);
0366             if (gp->meanDistPhiValue(iLayer, iRefLayer) != 0)
0367               mergedCnt++;
0368           }
0369
0370           if (mergedCnt) {
0371             meanDistPhi /= mergedCnt;
0372             //because for some gps the statistics can be too low, and then the meanDistPhiValue is 0, so it should not contribute
0373             for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
0374               auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
0375               gp->setMeanDistPhiValue(round(meanDistPhi), iLayer, iRefLayer);
0376               edm::LogImportant("l1tOmtfEventPrint")
0377                   << __FUNCTION__ << ": " << __LINE__ << " iGroup " << iGroup << " numInGroup " << i << " " << gp->key()
0378                   << " iLayer " << iLayer << " iRefLayer " << iRefLayer << " meanDistPhi after averaging "
0379                   << meanDistPhi << endl;
0380             }
0381           }
0382         }
0383       }
0384     }
0385   }
0386
0387   //calculating the pdfs
0388   for (auto& gp : goldenPatterns) {
0389     if (gp->key().thePt == 0)
0390       continue;
0391
0392     //TODO tune threshold <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
0393     int minHitCnt = 2 * minHitCntThresh * eventCntPerGp[gp->key().number()];
0394
0395     for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
0396       for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
0397         {
0398           double norm = 0;
0399           for (unsigned int iBin = 0; iBin < gp->getStatistics()[iLayer][iRefLayer].size();
0400                iBin++) {  //iBin = 0 i.e. no hit is included here, to have the proper norm
0401             norm += gp->getStatistics()[iLayer][iRefLayer][iBin][0];
0402           }
0403
0404           int pdfMiddle = gp->getPdf()[iLayer][iRefLayer].size() / 2;
0405           int statBinGroupSize = 1 << gp->getDistPhiBitShift(iLayer, iRefLayer);
0406           for (unsigned int iBinPdf = 0; iBinPdf < gp->getPdf()[iLayer][iRefLayer].size(); iBinPdf++) {
0407             double pdfVal = 0;
0408             if (iBinPdf > 0) {
0409               int groupedBins = 0;
0410               for (int i = 0; i < statBinGroupSize; i++) {
0411                 int iBinStat =
0412                     statBinGroupSize * ((int)(iBinPdf)-pdfMiddle) + i + gp->meanDistPhiValue(iLayer, iRefLayer);
0413
0414                 iBinStat += (gp->getStatistics()[iLayer][iRefLayer].size() / 2);
0415
0416                 if (iBinStat >= 0 && iBinStat < (int)gp->getStatistics()[iLayer][iRefLayer].size()) {
0417                   pdfVal += gp->getStatistics()[iLayer][iRefLayer][iBinStat][0];
0418                   groupedBins++;
0419                   //cout<<__FUNCTION__<<": "<<__LINE__<<" "<<gp->key()<<" iLayer "<<iLayer<<" iBinStat "<<iBinStat<<" iBinPdf "<<iBinPdf<<" statVal "<<gp->getStatistics()[iLayer][iRefLayer][iBinStat][0]<<endl;
0420                 }
0421               }
0422               if (norm > minHitCnt) {
0423                 pdfVal /= (norm * statBinGroupSize);
0424               } else
0425                 pdfVal = 0;
0426               /*edm::LogImportant("l1tOmtfEventPrint")
0427                                 << __FUNCTION__ << ": " << __LINE__ << " " << gp->key() << "calculating pdf: iLayer " << iLayer
0428                                 << " iRefLayer " << iRefLayer //<< " norm " << std::setw(5) << norm
0429                                 << " pdfVal " << pdfVal << endl;*/
0430             } else {  //iBinPdf == 0 i.e. no hit
0431               int iBinStat = 0;
0432               if (norm > 0) {
0433                 pdfVal = gp->getStatistics()[iLayer][iRefLayer][iBinStat][0] / norm;
0434               }
0435               edm::LogImportant("l1tOmtfEventPrint")
0436                   << __FUNCTION__ << ": " << __LINE__ << " " << gp->key() << "calculating pdf: iLayer " << iLayer
0437                   << " iRefLayer " << iRefLayer << " norm " << std::setw(5) << norm << " no hits cnt " << std::setw(5)
0438                   << gp->getStatistics()[iLayer][iRefLayer][iBinStat][0] << " pdfVal " << pdfVal << endl;
0439             }
0440
0441             double minPdfValFactor = 1;
0442             const double minPlog = log(omtfConfig->minPdfVal() * minPdfValFactor);
0443             const double pdfMaxVal = omtfConfig->pdfMaxValue();
0444
0445             int digitisedVal = 0;
0446             if (pdfVal >= omtfConfig->minPdfVal() * minPdfValFactor) {
0447               digitisedVal = rint(pdfMaxVal - log(pdfVal) / minPlog * pdfMaxVal);
0448             }
0449
0450             gp->setPdfValue(digitisedVal, iLayer, iRefLayer, iBinPdf);
0451             //cout<<__FUNCTION__<<": "<<__LINE__<<" "<<gp->key()<<" iLayer "<<iLayer<<" iBinPdf "<<iBinPdf<<" pdfVal "<<pdfVal<<" digitisedVal "<<digitisedVal<<endl;
0452           }
0453         }
0454       }
0455     }
0456   }
0457 }
0458
0459 void PatternGenerator::saveHists(TFile& outfile) {
0460   outfile.mkdir("ptDeltaPhiHists")->cd();
0461   //TODO uncomment if ptDeltaPhiHists are needed
0462   /*  for(unsigned int iCharge = 0; iCharge <= 1; iCharge++) {
0463     for(unsigned int iLayer = 0; iLayer < omtfConfig->nLayers(); ++iLayer) { //for the moment filing only ref layer, remove whe
0464       if(ptDeltaPhiHists[iCharge][iLayer]) {
0465         ptDeltaPhiHists[iCharge][iLayer]->Write();
0466       }
0467     }
0468   }*/
0469 }
0470
0471 void PatternGenerator::modifyClassProb(double step) {
0472   edm::LogImportant("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " Correcting P(C_k) " << std::endl;
0473   unsigned int iPdf = omtfConfig->nPdfBins() / 2;  // <<(omtfConfig->nPdfAddrBits()-1);
0474   for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns[0]->getPdf()[0].size(); ++iRefLayer) {
0475     unsigned int refLayerLogicNumber = omtfConfig->getRefToLogicNumber()[iRefLayer];
0476     if (iRefLayer == 0 || iRefLayer == 2)  //DT
0477       step = 1.5;
0478     else if (iRefLayer == 5)  //DT
0479       step = 1.5;
0480     else if (iRefLayer == 1)  //CSC
0481       step = 1.5;
0482     else if (iRefLayer == 3)  //CSC
0483       step = 1.5;
0484     else if (iRefLayer == 5)  //RE2/3
0485       step = 1.5;
0486     else if (iRefLayer == 6 || iRefLayer == 7)  //bRPC
0487       step = 1.5;
0488
0489     edm::LogImportant("l1tOmtfEventPrint")
0490         << __FUNCTION__ << ":" << __LINE__ << " RefLayer " << iRefLayer << " step " << step << std::endl;
0491     for (int sign = -1; sign <= 1; sign++) {
0492       for (auto& gp : boost::adaptors::reverse(goldenPatterns)) {
0493         if (gp->key().thePt == 0 || gp->key().theCharge != sign)
0494           continue;
0495
0496         double ptFrom = omtfConfig->getPatternPtRange(gp->key().theNumber).ptFrom;
0497         double ptTo = omtfConfig->getPatternPtRange(gp->key().theNumber).ptTo;
0498
0499         double ptRange = ptTo - ptFrom;
0500
0501         double minPdfValFactor = 0.1;
0502         double minPlog = log(omtfConfig->minPdfVal() * minPdfValFactor);
0503         double pdfMaxVal = omtfConfig->pdfMaxValue();
0504
0505         pdfMaxVal /= 3.;
0506         minPlog *= 2;
0507
0508         //last bin of the ptRange goes to 10000, so here we change it to 1000
0509         if (ptRange > 800)
0510           ptRange = 800;
0511
0512         double norm = 0.001;
0513         double classProb = vxIntegMuRate(ptFrom, ptRange, 0.82, 1.24) * norm;
0514
0515         int digitisedVal = rint(pdfMaxVal - log(classProb) / minPlog * pdfMaxVal);
0516
0517         int newPdfVal = digitisedVal;  //gp->getPdf()[refLayerLogicNumber][iRefLayer][iPdf]
0518
0519         if (ptFrom == 0)
0520           newPdfVal += 15;
0521         if (ptFrom == 3.5)
0522           newPdfVal += 15;
0523         if (ptFrom == 4)
0524           newPdfVal += 12;
0525         if (ptFrom == 4.5)
0526           newPdfVal += 9;
0527         if (ptFrom == 5)
0528           newPdfVal += 7;
0529         if (ptFrom == 6)
0530           newPdfVal += 4;
0531         if (ptFrom == 7)
0532           newPdfVal += 2;
0533
0534         if (ptFrom == 100)
0535           newPdfVal = 16;
0536         if (ptFrom == 200)
0537           newPdfVal = 22;
0538
0539         gp->setPdfValue(newPdfVal, refLayerLogicNumber, iRefLayer, iPdf);
0540
0541         edm::LogImportant("l1tOmtfEventPrint")
0542             << gp->key() << " " << omtfConfig->getPatternPtRange(gp->key().theNumber).ptFrom << " - "
0543             << omtfConfig->getPatternPtRange(gp->key().theNumber).ptTo << " GeV"
0544             << " ptRange " << ptRange << " RefLayer " << iRefLayer << " newPdfVal " << newPdfVal << std::endl;
0545       }
0546     }
0547   }
0548 }
0549
0550 void PatternGenerator::reCalibratePt() {
0551   edm::LogImportant("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " reCalibratePt" << std::endl;
0552   std::map<int, float> ptMap;
0553   //for Patterns_0x0009_oldSample_3_10Files_classProb2.xml
0554   ptMap[7] = 4.0;
0555   ptMap[8] = 4.5;
0556   ptMap[9] = 5.0;
0557   ptMap[10] = 5.5;
0558   ptMap[11] = 6.0;
0559   ptMap[13] = 7.0;
0560   ptMap[15] = 8.5;
0561   ptMap[17] = 10.0;
0562   ptMap[21] = 12.0;
0563   ptMap[25] = 14.0;
0564   ptMap[29] = 16.0;
0565   ptMap[33] = 18.5;
0566   ptMap[37] = 21.0;
0567   ptMap[41] = 23.0;
0568   ptMap[45] = 26.0;
0569   ptMap[49] = 28.0;
0570   ptMap[53] = 30.0;
0571   ptMap[57] = 32.0;
0572   ptMap[61] = 36.0;
0573   ptMap[71] = 40.0;
0574   ptMap[81] = 48.0;
0575   ptMap[91] = 54.0;
0576   ptMap[101] = 60.0;
0577   ptMap[121] = 70.0;
0578   ptMap[141] = 82.0;
0579   ptMap[161] = 96.0;
0580   ptMap[201] = 114.0;
0581   ptMap[401] = 200.0;
0582
0583   for (auto& gp : goldenPatterns) {
0584     if (gp->key().thePt == 0)
0585       continue;
0586
0587     int newPt = omtfConfig->ptGevToHw(ptMap[gp->key().thePt]);
0588     edm::LogImportant("l1tOmtfEventPrint") << gp->key().thePt << " -> " << newPt << std::endl;
0589
0590     gp->key().setPt(newPt);
0591   }
0592 }
0593
0594 void PatternGenerator::groupPatterns() {
0595   int group = 0;
0596   int indexInGroup = 0;
0597   for (auto& gp : goldenPatterns) {
0598     indexInGroup++;
0599     gp->key().setGroup(group);
0600     gp->key().setIndexInGroup(indexInGroup);
0601     //indexInGroup is counted from 1
0602
0603     edm::LogImportant("l1tOmtfEventPrint")
0604         << "setGroup(group): group " << group << " indexInGroup " << indexInGroup << std::endl;
0605
0606     if (gp->key().thePt <= 12 && indexInGroup == 2) {  //TODO
0607       indexInGroup = 0;
0608       group++;
0609     }
0610
0611     if (gp->key().thePt > 12 && indexInGroup == 4) {  //TODO
0612       indexInGroup = 0;
0613       group++;
0614     }
0615   }
0616
0617   OMTFConfiguration::vector2D patternGroups = omtfConfig->getPatternGroups(goldenPatterns);
0618   edm::LogImportant("l1tOmtfEventPrint") << "patternGroups:" << std::endl;
0619   for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
0620     edm::LogImportant("l1tOmtfEventPrint") << "patternGroup " << std::setw(2) << iGroup << " ";
0621     for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
0622       edm::LogImportant("l1tOmtfEventPrint") << i << " patNum " << patternGroups[iGroup][i] << " ";
0623     }
0624     edm::LogImportant("l1tOmtfEventPrint") << std::endl;
0625   }
0626
0627   int pdfBins = exp2(omtfConfig->nPdfAddrBits());
0628
0629   for (unsigned int iLayer = 0; iLayer < goldenPatterns.at(0)->getPdf().size(); ++iLayer) {
0630     for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns.at(0)->getPdf()[iLayer].size(); ++iRefLayer) {
0631       //unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefLayer];
0632       //if(refLayerLogicNum == iLayer)
0633       {
0634         //averaging the meanDistPhi for the gp belonging to the same group
0635         for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
0636           double meanDistPhi = 0;
0637           int mergedCnt = 0;
0638           for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
0639             auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
0640             meanDistPhi += gp->meanDistPhiValue(iLayer, iRefLayer);
0641             if (gp->meanDistPhiValue(iLayer, iRefLayer) != 0)
0642               mergedCnt++;
0643             edm::LogImportant("l1tOmtfEventPrint")
0644                 << __FUNCTION__ << ": " << __LINE__ << " iGroup " << iGroup << " numInGroup " << i << " " << gp->key()
0645                 << " iLayer " << iLayer << " iRefLayer " << iRefLayer << " old meanDistPhiValue "
0646                 << gp->meanDistPhiValue(iLayer, iRefLayer) << endl;
0647           }
0648
0649           if (mergedCnt) {
0650             meanDistPhi /= mergedCnt;
0651             meanDistPhi = (int)meanDistPhi;
0652
0653             //because for some gps the statistics can be too low, and then the meanDistPhiValue is 0, so it should not contribute
0654             for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
0655               auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
0656               unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefLayer];
0657               if (refLayerLogicNum != iLayer) {
0658                 int shift = meanDistPhi - gp->meanDistPhiValue(iLayer, iRefLayer);
0659                 edm::LogImportant("l1tOmtfEventPrint")
0660                     << __FUNCTION__ << ": " << __LINE__ << " iGroup " << iGroup << " numInGroup " << i << " "
0661                     << gp->key() << " iLayer " << iLayer << " iRefLayer " << iRefLayer
0662                     << " new meanDistPhi after averaging " << meanDistPhi << " old meanDistPhiValue "
0663                     << gp->meanDistPhiValue(iLayer, iRefLayer) << " shift " << shift << endl;
0664
0665                 if (shift < 0) {
0666                   for (int iBin = 1 - shift; iBin < pdfBins;
0667                        iBin++) {  //iBin = 0 i.e. no hit is included here, to have the proper norm
0668                     auto pdfVal = gp->pdfValue(iLayer, iRefLayer, iBin);
0669                     gp->setPdfValue(pdfVal, iLayer, iRefLayer, iBin + shift);
0670                     edm::LogImportant("l1tOmtfEventPrint")
0671                         << "     iBin " << iBin << "  iBin + shift " << iBin + shift << " pdfVal " << pdfVal << endl;
0672                   }
0673                   for (int iBin = pdfBins + shift; iBin < pdfBins; iBin++) {
0674                     gp->setPdfValue(0, iLayer, iRefLayer, iBin);
0675                   }
0676                 } else if (shift > 0) {
0677                   for (int iBin = pdfBins - 1 - shift; iBin > 0;
0678                        iBin--) {  //iBin = 0 i.e. no hit is included here, to have the proper norm
0679                     auto pdfVal = gp->pdfValue(iLayer, iRefLayer, iBin);
0680                     gp->setPdfValue(pdfVal, iLayer, iRefLayer, iBin + shift);
0681                     edm::LogImportant("l1tOmtfEventPrint")
0682                         << "     iBin " << iBin << "  iBin + shift " << iBin + shift << " pdfVal " << pdfVal << endl;
0683                   }
0684                   for (int iBin = shift; iBin > 0; iBin--) {
0685                     gp->setPdfValue(0, iLayer, iRefLayer, iBin);
0686                   }
0687                 }
0688               }
0689
0690               gp->setMeanDistPhiValue(round(meanDistPhi), iLayer, iRefLayer);
0691             }
0692           }
0693         }
0694       }
0695     }
0696   }
0697 }