File indexing completed on 2024-05-15 04:21:52
0001
0002
0003
0004
0005
0006
0007 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Tools/PatternOptimizerBase.h"
0008 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/AlgoMuon.h"
0009 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/GoldenPattern.h"
0010 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/GoldenPatternBase.h"
0011 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/GoldenPatternResult.h"
0012 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/OMTFConfiguration.h"
0013 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/XMLConfigWriter.h"
0014 #include "L1Trigger/RPCTrigger/interface/RPCConst.h"
0015
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "DataFormats/L1TMuon/interface/RegionalMuonCand.h"
0018 #include "DataFormats/L1TMuon/interface/RegionalMuonCandFwd.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "SimDataFormats/Track/interface/CoreSimTrack.h"
0022 #include "SimDataFormats/Track/interface/SimTrack.h"
0023
0024 #include "TCanvas.h"
0025 #include "TFile.h"
0026 #include "TH1.h"
0027 #include "TH2.h"
0028 #include "TH2F.h"
0029 #include "TStyle.h"
0030 #include <cstdlib>
0031 #include <iostream>
0032 #include <memory>
0033 #include <sstream>
0034 #include <string>
0035 #include <vector>
0036
0037 double PatternOptimizerBase::vxMuRate(double pt_GeV) {
0038 if (pt_GeV == 0)
0039 return 0.0;
0040 const double lum = 2.0e34;
0041 const double dabseta = 1.0;
0042 const double dpt = 1.0;
0043 const double afactor = 1.0e-34 * lum * dabseta * dpt;
0044 const double a = 2 * 1.3084E6;
0045 const double mu = -0.725;
0046 const double sigma = 0.4333;
0047 const double s2 = 2 * sigma * sigma;
0048
0049 double ptlog10;
0050 ptlog10 = log10(pt_GeV);
0051 double ex = (ptlog10 - mu) * (ptlog10 - mu) / s2;
0052 double rate = (a * exp(-ex) * afactor);
0053
0054 return rate;
0055 }
0056
0057 double PatternOptimizerBase::vxIntegMuRate(double pt_GeV, double dpt, double etaFrom, double etaTo) {
0058
0059 double rate = 0.5 * (vxMuRate(pt_GeV) + vxMuRate(pt_GeV + dpt)) * dpt;
0060
0061 rate = rate * (etaTo - etaFrom);
0062
0063 return rate;
0064 }
0065
0066 PatternOptimizerBase::PatternOptimizerBase(const edm::ParameterSet& edmCfg,
0067 const OMTFConfiguration* omtfConfig,
0068 GoldenPatternVec<GoldenPatternWithStat>& gps)
0069 : EmulationObserverBase(edmCfg, omtfConfig), goldenPatterns(gps) {
0070
0071
0072 simMuPt = new TH1I("simMuPt", "simMuPt", goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5);
0073 simMuFoundByOmtfPt =
0074 new TH1I("simMuFoundByOmtfPt", "simMuFoundByOmtfPt", goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5);
0075
0076 simMuEta = new TH1I("simMuEta", "simMuEta;eta;#events", 100, -2.1, 2.1);
0077 candEta = new TH1I("candEta", "candEta;eta;#events", 100, -2.1, 2.1);
0078
0079 simMuPtSpectrum = new TH1F("simMuPtSpectrum", "simMuPtSpectrum", 800, 0, 400);
0080
0081 simMuPtVsDispl = new TH2I("simMuPtVsDispl", "simMuPtVsDispl;pt [GeV];dxy [cm]", 100, 0, 400, 100, 0, 400);
0082 simMuPtVsRho = new TH2I("simMuPtVsRho", "simMuPtVsRho;pt [GeV];rho [cm]", 100, 0, 400, 100, 0, 400);
0083
0084 if (edmCfg.exists("simTracksTag") == false)
0085 edm::LogError("l1tOmtfEventPrint") << "simTracksTag not found !!!" << std::endl;
0086 }
0087
0088 PatternOptimizerBase::~PatternOptimizerBase() {}
0089
0090 void PatternOptimizerBase::printPatterns() {
0091 edm::LogVerbatim("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " called!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! "
0092 << std::endl;
0093 for (int patNum = goldenPatterns.size() - 1; patNum >= 0; patNum--) {
0094 double pt = omtfConfig->getPatternPtRange(patNum).ptFrom;
0095 if (pt > 0) {
0096 edm::LogVerbatim("l1tOmtfEventPrint")
0097 << "cmsRun runThresholdCalc.py " << patNum << " " << (patNum + 1) << " _" << RPCConst::iptFromPt(pt) << "_";
0098 if (goldenPatterns[patNum]->key().theCharge == -1)
0099 edm::LogVerbatim("l1tOmtfEventPrint") << "m_";
0100 else
0101 edm::LogVerbatim("l1tOmtfEventPrint") << "p_";
0102
0103 edm::LogVerbatim("l1tOmtfEventPrint") << " > out" << patNum << ".txt" << std::endl;
0104 }
0105 }
0106 }
0107
0108 void PatternOptimizerBase::observeEventEnd(const edm::Event& iEvent,
0109 std::unique_ptr<l1t::RegionalMuonCandBxCollection>& finalCandidates) {
0110 if (simMuon == nullptr || omtfCand->getGoldenPatern() == nullptr)
0111 return;
0112
0113 double ptSim = simMuon->momentum().pt();
0114 int chargeSim = (std::abs(simMuon->type()) == 13) ? simMuon->type() / -13 : 0;
0115
0116
0117
0118 unsigned int exptPatNum = omtfConfig->getPatternNum(ptSim, chargeSim);
0119 GoldenPatternWithStat* exptCandGp = goldenPatterns.at(exptPatNum).get();
0120 simMuFoundByOmtfPt->Fill(exptCandGp->key().theNumber);
0121
0122 simMuPtSpectrum->Fill(ptSim, getEventRateWeight(ptSim));
0123 }
0124
0125 void PatternOptimizerBase::endJob() {
0126 std::string fName = edmCfg.getParameter<std::string>("optimisedPatsXmlFile");
0127 edm::LogImportant("PatternOptimizer") << " Writing optimized patterns to " << fName << std::endl;
0128 XMLConfigWriter xmlWriter(omtfConfig, false, false);
0129 xmlWriter.writeGPs(goldenPatterns, fName);
0130
0131 fName.replace(fName.find('.'), fName.length(), ".root");
0132 savePatternsInRoot(fName);
0133 }
0134
0135 void PatternOptimizerBase::savePatternsInRoot(std::string rootFileName) {
0136 gStyle->SetOptStat(111111);
0137 TFile outfile(rootFileName.c_str(), "RECREATE");
0138 edm::LogVerbatim("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " out fileName " << rootFileName
0139 << " outfile->GetName() " << outfile.GetName() << " writeLayerStat "
0140 << writeLayerStat << endl;
0141
0142 outfile.cd();
0143 simMuFoundByOmtfPt->Write();
0144
0145 simMuEta->Write();
0146 candEta->Write();
0147
0148 simMuPtSpectrum->Write();
0149 simMuPtVsDispl->Write();
0150 simMuPtVsRho->Write();
0151
0152 outfile.mkdir("patternsPdfs")->cd();
0153 outfile.mkdir("patternsPdfs/canvases");
0154 outfile.mkdir("patternsPdfs/canvases2");
0155 outfile.mkdir("layerStats");
0156 ostringstream ostrName;
0157 ostringstream ostrTtle;
0158 vector<TH1F*> classProbHists;
0159 for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns[0]->getPdf()[0].size(); ++iRefLayer) {
0160 ostrName.str("");
0161 ostrName << "Neg_RefLayer_" << iRefLayer;
0162 ostrTtle.str("");
0163 ostrTtle << "Neg_RefLayer_" << iRefLayer;
0164 classProbHists.push_back(new TH1F(
0165 ostrName.str().c_str(), ostrTtle.str().c_str(), goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5));
0166
0167 ostrName.str("");
0168 ostrName << "Pos_RefLayer_" << iRefLayer;
0169 ostrTtle.str("");
0170 ostrTtle << "Pos_RefLayer_" << iRefLayer;
0171 classProbHists.push_back(new TH1F(
0172 ostrName.str().c_str(), ostrTtle.str().c_str(), goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5));
0173 }
0174
0175 auto gpsCnt = goldenPatterns.size();
0176 auto layerCnt = goldenPatterns[0]->getPdf().size();
0177 auto refLayerCnt = goldenPatterns[0]->getPdf().size();
0178
0179 vector<vector<TH2F*> > distPhiLayerRefLayer(layerCnt, vector<TH2F*>(refLayerCnt, nullptr));
0180 vector<vector<TH2F*> > meanDistPhiLayerRefLayer(layerCnt, vector<TH2F*>(refLayerCnt, nullptr));
0181 vector<vector<TH2F*> > pdfsLayerRefLayer(layerCnt, vector<TH2F*>(refLayerCnt, nullptr));
0182
0183 for (unsigned int iLayer = 0; iLayer < layerCnt; ++iLayer) {
0184 int rangeFactor = 1;
0185 if (iLayer == 1 || iLayer == 3 || iLayer == 5)
0186 rangeFactor = 2;
0187
0188 for (unsigned int iRefLayer = 0; iRefLayer < refLayerCnt; ++iRefLayer) {
0189 ostrName.str("");
0190 ostrName << "distPhi_refLayer_" << iRefLayer << "_layer_" << iLayer;
0191 ostrTtle.str("");
0192 ostrTtle << "distPhi refLayer " << iRefLayer << " layer " << iLayer;
0193
0194 distPhiLayerRefLayer[iLayer][iRefLayer] = new TH2F(ostrName.str().c_str(),
0195 ostrTtle.str().c_str(),
0196 gpsCnt,
0197 0,
0198 gpsCnt,
0199 omtfConfig->nPdfBins() * rangeFactor * 2,
0200 (int)(omtfConfig->nPdfBins()) * (-rangeFactor) - 0.5,
0201 omtfConfig->nPdfBins() * rangeFactor - 0.5);
0202
0203 ostrName.str("");
0204 ostrName << "meanDistPhi_refLayer_" << iRefLayer << "_Layer_" << iLayer;
0205 ostrTtle.str("");
0206 ostrTtle << "meanDistPhi refLayer " << iRefLayer << " Layer " << iLayer;
0207
0208 meanDistPhiLayerRefLayer[iLayer][iRefLayer] = new TH2F(ostrName.str().c_str(),
0209 ostrTtle.str().c_str(),
0210 gpsCnt,
0211 0,
0212 gpsCnt,
0213 omtfConfig->nPdfBins() * rangeFactor * 2,
0214 (int)(omtfConfig->nPdfBins()) * (-rangeFactor) - 0.5,
0215 omtfConfig->nPdfBins() * rangeFactor - 0.5);
0216
0217 ostrName.str("");
0218 ostrName << "pdfs_refLayer_" << iRefLayer << "_layer_" << iLayer;
0219 ostrTtle.str("");
0220 ostrTtle << "pdfs refLayer " << iRefLayer << " layer " << iLayer;
0221
0222 pdfsLayerRefLayer[iLayer][iRefLayer] = new TH2F(ostrName.str().c_str(),
0223 ostrTtle.str().c_str(),
0224 gpsCnt,
0225 0,
0226 gpsCnt,
0227 omtfConfig->nPdfBins(),
0228 -0.5,
0229 omtfConfig->nPdfBins() - 0.5);
0230 }
0231 }
0232
0233 for (auto& gp : goldenPatterns) {
0234 OMTFConfiguration::PatternPt patternPt = omtfConfig->getPatternPtRange(gp->key().theNumber);
0235 if (gp->key().thePt == 0)
0236 continue;
0237
0238 ostrName.str("");
0239 ostrName << "PatNum_" << gp->key().theNumber;
0240 ostrTtle.str("");
0241 ostrTtle << "PatNum_" << gp->key().theNumber << "_ptCode_" << gp->key().thePt << "_Pt_" << patternPt.ptFrom << "_"
0242 << patternPt.ptTo << "_GeV";
0243 TCanvas* canvas = new TCanvas(ostrName.str().c_str(), ostrTtle.str().c_str(), 1200, 1000);
0244 canvas->Divide(gp->getPdf().size(), gp->getPdf()[0].size(), 0, 0);
0245
0246 for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[0].size(); ++iRefLayer) {
0247 for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
0248 canvas->cd(1 + iLayer + iRefLayer * gp->getPdf().size());
0249 ostrName.str("");
0250 ostrName << "PatNum_" << gp->key().theNumber << "_refLayer_" << iRefLayer << "_Layer_" << iLayer;
0251 ostrTtle.str("");
0252 ostrTtle << "PatNum " << gp->key().theNumber << " ptCode " << gp->key().thePt << " refLayer " << iRefLayer
0253 << " Layer " << iLayer << " meanDistPhi " << gp->meanDistPhi[iLayer][iRefLayer][0]
0254 << " distPhiBitShift "
0255 << gp->getDistPhiBitShift(iLayer, iRefLayer);
0256
0257 TH1F* hist = new TH1F(
0258 ostrName.str().c_str(), ostrTtle.str().c_str(), omtfConfig->nPdfBins(), -0.5, omtfConfig->nPdfBins() - 0.5);
0259
0260 int pdfMiddle = gp->getPdf()[iLayer][iRefLayer].size() / 2;
0261 int shift = gp->getDistPhiBitShift(iLayer, iRefLayer);
0262
0263 for (unsigned int iPdf = 0; iPdf < gp->getPdf()[iLayer][iRefLayer].size(); iPdf++) {
0264 hist->Fill(iPdf, gp->pdfAllRef[iLayer][iRefLayer][iPdf]);
0265
0266 distPhiLayerRefLayer[iLayer][iRefLayer]->Fill(
0267 gp->key().theNumber,
0268 (((int)iPdf - pdfMiddle) << shift) + gp->meanDistPhi[iLayer][iRefLayer][0],
0269 gp->pdfAllRef[iLayer][iRefLayer][iPdf]);
0270
0271 pdfsLayerRefLayer[iLayer][iRefLayer]->Fill(
0272 gp->key().theNumber, (int)iPdf, gp->pdfAllRef[iLayer][iRefLayer][iPdf]);
0273 }
0274 if ((int)iLayer == (omtfConfig->getRefToLogicNumber()[iRefLayer]))
0275 hist->SetLineColor(kGreen);
0276
0277 meanDistPhiLayerRefLayer[iLayer][iRefLayer]->Fill(
0278 gp->key().theNumber, gp->meanDistPhi[iLayer][iRefLayer][0], 1);
0279
0280 hist->GetYaxis()->SetRangeUser(0, omtfConfig->pdfMaxValue() + 1);
0281 hist->Draw("hist");
0282
0283 outfile.cd("patternsPdfs");
0284 hist->Write();
0285
0286
0287 if (writeLayerStat) {
0288 bool saveTh2 = false;
0289 if (gp->getStatistics()[iLayer][iRefLayer][0].size() > 1)
0290 saveTh2 = true;
0291
0292 outfile.cd("layerStats");
0293
0294 string histName = "histLayerStat_" + ostrName.str();
0295 unsigned int binCnt1 = gp->getStatistics()[iLayer][iRefLayer].size();
0296 if (!saveTh2) {
0297 TH1I* histLayerStat = new TH1I(histName.c_str(), histName.c_str(), binCnt1, -0.5, binCnt1 - 0.5);
0298 for (unsigned int iBin = 0; iBin < binCnt1; iBin++) {
0299 histLayerStat->Fill(iBin, gp->getStatistics()[iLayer][iRefLayer][iBin][0]);
0300 }
0301 histLayerStat->Write();
0302 } else {
0303 if (iRefLayer == 0 || iRefLayer == 2) {
0304 unsigned int binCnt2 = gp->getStatistics()[iLayer][iRefLayer][0].size();
0305
0306 double xmin = -0.5 - binCnt2 / 2;
0307 double xmax = binCnt2 / 2 - 0.5;
0308 if (edmCfg.getParameter<string>("patternGenerator") == "deltaPhiVsPhiRef") {
0309 xmin = -0.5;
0310 xmax = binCnt2 - 0.5;
0311 }
0312 TH2I* histLayerStat = new TH2I(histName.c_str(),
0313 (histName + ";ref phiB;delta_phi").c_str(),
0314 binCnt2,
0315 xmin,
0316 xmax,
0317 binCnt1,
0318 -0.5 - binCnt1 / 2,
0319 binCnt1 / 2 - 0.5);
0320
0321 if (edmCfg.getParameter<string>("patternGenerator") == "deltaPhiVsPhiRef")
0322 histLayerStat->GetXaxis()->SetTitle("ref phi");
0323
0324 for (unsigned int iBin1 = 0; iBin1 < binCnt1; iBin1++) {
0325 for (unsigned int iBin2 = 0; iBin2 < binCnt2; iBin2++) {
0326
0327 histLayerStat->SetBinContent(
0328 iBin2 + 1, iBin1 + 1, gp->getStatistics()[iLayer][iRefLayer][iBin1][iBin2]);
0329 }
0330 }
0331 histLayerStat->Write();
0332 histLayerStat->Delete();
0333 }
0334 }
0335 }
0336 }
0337 }
0338 outfile.cd("patternsPdfs/canvases");
0339 canvas->Write();
0340 delete canvas;
0341
0342 unsigned int iPdf = omtfConfig->nPdfBins() / 2;
0343 for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[0].size(); ++iRefLayer) {
0344 unsigned int refLayerLogicNumber = omtfConfig->getRefToLogicNumber()[iRefLayer];
0345 if (gp->key().theCharge == -1) {
0346 classProbHists[2 * iRefLayer]->Fill(gp->key().theNumber, gp->pdfAllRef[refLayerLogicNumber][iRefLayer][iPdf]);
0347 } else
0348 classProbHists[2 * iRefLayer + 1]->Fill(gp->key().theNumber,
0349 gp->pdfAllRef[refLayerLogicNumber][iRefLayer][iPdf]);
0350 }
0351 }
0352
0353 outfile.cd();
0354 for (auto& classProbHist : classProbHists) {
0355 classProbHist->Write();
0356 }
0357
0358 outfile.cd("patternsPdfs/canvases2");
0359 for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns[0]->getPdf()[0].size(); ++iRefLayer) {
0360 ostrName.str("");
0361 ostrName << "distPhiForPatterns_Reflayer_" << iRefLayer;
0362 ostrTtle.str("");
0363 ostrTtle << "distPhiForPatterns Reflayer " << iRefLayer;
0364
0365 TCanvas* canvas = new TCanvas(ostrName.str().c_str(), ostrTtle.str().c_str(), 1200, 1000);
0366 canvas->Divide(6, 3, 0, 0);
0367
0368 for (unsigned int iLayer = 0; iLayer < distPhiLayerRefLayer.size(); ++iLayer) {
0369 canvas->cd(iLayer + 1);
0370 canvas->cd(iLayer + 1)->SetGridx();
0371 canvas->cd(iLayer + 1)->SetGridy();
0372
0373
0374 distPhiLayerRefLayer[iLayer][iRefLayer]->Draw("colz");
0375
0376 meanDistPhiLayerRefLayer[iLayer][iRefLayer]->SetLineColor(kRed);
0377 meanDistPhiLayerRefLayer[iLayer][iRefLayer]->Draw("boxsame");
0378
0379 distPhiLayerRefLayer[iLayer][iRefLayer]->Write();
0380 meanDistPhiLayerRefLayer[iLayer][iRefLayer]->Write();
0381 pdfsLayerRefLayer[iLayer][iRefLayer]->Write();
0382 }
0383
0384 canvas->Write();
0385
0386 ostrName.str("");
0387 ostrName << "pdfsForPatterns_Reflayer_" << iRefLayer;
0388 ostrTtle.str("");
0389 ostrTtle << "pdfsForPatterns Reflayer " << iRefLayer;
0390
0391 canvas = new TCanvas(ostrName.str().c_str(), ostrTtle.str().c_str(), 1200, 1000);
0392 canvas->Divide(6, 3, 0, 0);
0393 for (unsigned int iLayer = 0; iLayer < distPhiLayerRefLayer.size(); ++iLayer) {
0394 canvas->cd(iLayer + 1);
0395 canvas->cd(iLayer + 1)->SetGridx();
0396 canvas->cd(iLayer + 1)->SetGridy();
0397
0398 pdfsLayerRefLayer[iLayer][iRefLayer]->Draw("colz");
0399 }
0400
0401 canvas->Write();
0402 }
0403
0404 saveHists(outfile);
0405
0406 outfile.Close();
0407 }