File indexing completed on 2023-10-25 09:55:39
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 "TStyle.h"
0028 #include <cstdlib>
0029 #include <iostream>
0030 #include <memory>
0031 #include <sstream>
0032 #include <string>
0033 #include <vector>
0034
0035 double PatternOptimizerBase::vxMuRate(double pt_GeV) {
0036 if (pt_GeV == 0)
0037 return 0.0;
0038 const double lum = 2.0e34;
0039 const double dabseta = 1.0;
0040 const double dpt = 1.0;
0041 const double afactor = 1.0e-34 * lum * dabseta * dpt;
0042 const double a = 2 * 1.3084E6;
0043 const double mu = -0.725;
0044 const double sigma = 0.4333;
0045 const double s2 = 2 * sigma * sigma;
0046
0047 double ptlog10;
0048 ptlog10 = log10(pt_GeV);
0049 double ex = (ptlog10 - mu) * (ptlog10 - mu) / s2;
0050 double rate = (a * exp(-ex) * afactor);
0051
0052 return rate;
0053 }
0054
0055 double PatternOptimizerBase::vxIntegMuRate(double pt_GeV, double dpt, double etaFrom, double etaTo) {
0056
0057 double rate = 0.5 * (vxMuRate(pt_GeV) + vxMuRate(pt_GeV + dpt)) * dpt;
0058
0059 rate = rate * (etaTo - etaFrom);
0060
0061 return rate;
0062 }
0063
0064 PatternOptimizerBase::PatternOptimizerBase(const edm::ParameterSet& edmCfg,
0065 const OMTFConfiguration* omtfConfig,
0066 GoldenPatternVec<GoldenPatternWithStat>& gps)
0067 : EmulationObserverBase(edmCfg, omtfConfig), goldenPatterns(gps) {
0068
0069
0070 simMuPt = new TH1I("simMuPt", "simMuPt", goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5);
0071 simMuFoundByOmtfPt =
0072 new TH1I("simMuFoundByOmtfPt", "simMuFoundByOmtfPt", goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5);
0073
0074 simMuPtSpectrum = new TH1F("simMuPtSpectrum", "simMuPtSpectrum", 800, 0, 400);
0075
0076 if (edmCfg.exists("simTracksTag") == false)
0077 edm::LogError("l1tOmtfEventPrint") << "simTracksTag not found !!!" << std::endl;
0078 }
0079
0080 PatternOptimizerBase::~PatternOptimizerBase() {}
0081
0082 void PatternOptimizerBase::printPatterns() {
0083 edm::LogVerbatim("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " called!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! "
0084 << std::endl;
0085 for (int patNum = goldenPatterns.size() - 1; patNum >= 0; patNum--) {
0086 double pt = omtfConfig->getPatternPtRange(patNum).ptFrom;
0087 if (pt > 0) {
0088 edm::LogVerbatim("l1tOmtfEventPrint")
0089 << "cmsRun runThresholdCalc.py " << patNum << " " << (patNum + 1) << " _" << RPCConst::iptFromPt(pt) << "_";
0090 if (goldenPatterns[patNum]->key().theCharge == -1)
0091 edm::LogVerbatim("l1tOmtfEventPrint") << "m_";
0092 else
0093 edm::LogVerbatim("l1tOmtfEventPrint") << "p_";
0094
0095 edm::LogVerbatim("l1tOmtfEventPrint") << " > out" << patNum << ".txt" << std::endl;
0096 }
0097 }
0098 }
0099
0100 void PatternOptimizerBase::observeEventEnd(const edm::Event& iEvent,
0101 std::unique_ptr<l1t::RegionalMuonCandBxCollection>& finalCandidates) {
0102 if (simMuon == nullptr || omtfCand->getGoldenPatern() == nullptr)
0103 return;
0104
0105 double ptSim = simMuon->momentum().pt();
0106 int chargeSim = (abs(simMuon->type()) == 13) ? simMuon->type() / -13 : 0;
0107
0108 unsigned int exptPatNum = omtfConfig->getPatternNum(ptSim, chargeSim);
0109 GoldenPatternWithStat* exptCandGp = goldenPatterns.at(exptPatNum).get();
0110 simMuFoundByOmtfPt->Fill(exptCandGp->key().theNumber);
0111
0112 simMuPtSpectrum->Fill(ptSim, getEventRateWeight(ptSim));
0113 }
0114
0115 void PatternOptimizerBase::endJob() {
0116 std::string fName = edmCfg.getParameter<std::string>("optimisedPatsXmlFile");
0117 edm::LogImportant("PatternOptimizer") << " Writing optimized patterns to " << fName << std::endl;
0118 XMLConfigWriter xmlWriter(omtfConfig, false, false);
0119 xmlWriter.writeGPs(goldenPatterns, fName);
0120
0121 fName.replace(fName.find('.'), fName.length(), ".root");
0122 savePatternsInRoot(fName);
0123 }
0124
0125 void PatternOptimizerBase::savePatternsInRoot(std::string rootFileName) {
0126 gStyle->SetOptStat(111111);
0127 TFile outfile(rootFileName.c_str(), "RECREATE");
0128 edm::LogVerbatim("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " out fileName " << rootFileName
0129 << " outfile->GetName() " << outfile.GetName() << " writeLayerStat "
0130 << writeLayerStat << endl;
0131
0132 outfile.cd();
0133 simMuFoundByOmtfPt->Write();
0134
0135 simMuPtSpectrum->Write();
0136
0137 outfile.mkdir("patternsPdfs")->cd();
0138 outfile.mkdir("patternsPdfs/canvases");
0139 outfile.mkdir("layerStats");
0140 ostringstream ostrName;
0141 ostringstream ostrTtle;
0142 vector<TH1F*> classProbHists;
0143 for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns[0]->getPdf()[0].size(); ++iRefLayer) {
0144 ostrName.str("");
0145 ostrName << "Neg_RefLayer_" << iRefLayer;
0146 ostrTtle.str("");
0147 ostrTtle << "Neg_RefLayer_" << iRefLayer;
0148 classProbHists.push_back(new TH1F(
0149 ostrName.str().c_str(), ostrTtle.str().c_str(), goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5));
0150
0151 ostrName.str("");
0152 ostrName << "Pos_RefLayer_" << iRefLayer;
0153 ostrTtle.str("");
0154 ostrTtle << "Pos_RefLayer_" << iRefLayer;
0155 classProbHists.push_back(new TH1F(
0156 ostrName.str().c_str(), ostrTtle.str().c_str(), goldenPatterns.size(), -0.5, goldenPatterns.size() - 0.5));
0157 }
0158
0159 for (auto& gp : goldenPatterns) {
0160 OMTFConfiguration::PatternPt patternPt = omtfConfig->getPatternPtRange(gp->key().theNumber);
0161 if (gp->key().thePt == 0)
0162 continue;
0163
0164 ostrName.str("");
0165 ostrName << "PatNum_" << gp->key().theNumber;
0166 ostrTtle.str("");
0167 ostrTtle << "PatNum_" << gp->key().theNumber << "_ptCode_" << gp->key().thePt << "_Pt_" << patternPt.ptFrom << "_"
0168 << patternPt.ptTo << "_GeV";
0169 TCanvas* canvas = new TCanvas(ostrName.str().c_str(), ostrTtle.str().c_str(), 1200, 1000);
0170 canvas->Divide(gp->getPdf().size(), gp->getPdf()[0].size(), 0, 0);
0171
0172 for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
0173 for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
0174 canvas->cd(1 + iLayer + iRefLayer * gp->getPdf().size());
0175 ostrName.str("");
0176 ostrName << "PatNum_" << gp->key().theNumber << "_refLayer_" << iRefLayer << "_Layer_" << iLayer;
0177 ostrTtle.str("");
0178 ostrTtle << "PatNum " << gp->key().theNumber << " ptCode " << gp->key().thePt << " refLayer " << iRefLayer
0179 << " Layer " << iLayer << " meanDistPhi " << gp->meanDistPhi[iLayer][iRefLayer][0]
0180 << " distPhiBitShift "
0181 << gp->getDistPhiBitShift(iLayer, iRefLayer);
0182
0183 TH1F* hist = new TH1F(
0184 ostrName.str().c_str(), ostrTtle.str().c_str(), omtfConfig->nPdfBins(), -0.5, omtfConfig->nPdfBins() - 0.5);
0185 for (unsigned int iPdf = 0; iPdf < gp->getPdf()[iLayer][iRefLayer].size(); iPdf++) {
0186 hist->Fill(iPdf, gp->pdfAllRef[iLayer][iRefLayer][iPdf]);
0187 }
0188 if ((int)iLayer == (omtfConfig->getRefToLogicNumber()[iRefLayer]))
0189 hist->SetLineColor(kGreen);
0190
0191 hist->GetYaxis()->SetRangeUser(0, omtfConfig->pdfMaxValue() + 1);
0192 hist->Draw("hist");
0193
0194 outfile.cd("patternsPdfs");
0195 hist->Write();
0196
0197
0198 if (writeLayerStat) {
0199 string histName = "histLayerStat_" + ostrName.str();
0200 unsigned int binCnt = gp->getStatistics()[iLayer][iRefLayer].size();
0201 TH1I* histLayerStat = new TH1I(histName.c_str(), histName.c_str(), binCnt, -0.5, binCnt - 0.5);
0202 for (unsigned int iBin = 0; iBin < binCnt; iBin++) {
0203 histLayerStat->Fill(iBin, gp->getStatistics()[iLayer][iRefLayer][iBin][0]);
0204 }
0205
0206 outfile.cd("layerStats");
0207 histLayerStat->Write();
0208 }
0209 }
0210 }
0211 outfile.cd("patternsPdfs/canvases");
0212 canvas->Write();
0213 delete canvas;
0214
0215 unsigned int iPdf = omtfConfig->nPdfBins() / 2;
0216 for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[0].size(); ++iRefLayer) {
0217 unsigned int refLayerLogicNumber = omtfConfig->getRefToLogicNumber()[iRefLayer];
0218 if (gp->key().theCharge == -1) {
0219 classProbHists[2 * iRefLayer]->Fill(gp->key().theNumber, gp->pdfAllRef[refLayerLogicNumber][iRefLayer][iPdf]);
0220 } else
0221 classProbHists[2 * iRefLayer + 1]->Fill(gp->key().theNumber,
0222 gp->pdfAllRef[refLayerLogicNumber][iRefLayer][iPdf]);
0223 }
0224 }
0225
0226 outfile.cd();
0227 for (auto& classProbHist : classProbHists) {
0228 classProbHist->Write();
0229 }
0230
0231 saveHists(outfile);
0232
0233 outfile.Close();
0234 }