Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:55:39 

0001 /*
0002  * PatternOptimizerBase.cc
0003  *
0004  *  Created on: Oct 17, 2018
0005  *      Author: kbunkow
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;  //defoult is 1.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   //edm::LogError("l1tOmtfEventPrint")<<ptCode<<" "<<rate;//<<<<<<<<<<<<<<<<<<<<<<<<
0052   return rate;
0053 }
0054 
0055 double PatternOptimizerBase::vxIntegMuRate(double pt_GeV, double dpt, double etaFrom, double etaTo) {
0056   //integration using trapeze method - not exact but good enough
0057   double rate = 0.5 * (vxMuRate(pt_GeV) + vxMuRate(pt_GeV + dpt)) * dpt;
0058 
0059   rate = rate * (etaTo - etaFrom);
0060   //edm::LogError("l1tOmtfEventPrint")<<ptCode<<" "<<rate;//<<<<<<<<<<<<<<<<<<<<<<<<
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   // TODO Auto-generated constructor stub
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)  //no sim muon or empty candidate
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();  // expected pattern
0110   simMuFoundByOmtfPt->Fill(exptCandGp->key().theNumber);                    //TODO add weight of the muons pt spectrum
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     //edm::LogVerbatim("l1tOmtfEventPrint") <<__FUNCTION__<<": "<<__LINE__<<" "<<gp->key()<<std::endl;
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);  //"_Pt_"<<patternPt.ptFrom<<"_"<<patternPt.ptTo<<"_GeV
0182         //edm::LogVerbatim("l1tOmtfEventPrint") <<__FUNCTION__<<": "<<__LINE__<<" creating hist "<<ostrTtle.str()<<std::endl;
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         /////////////////////// histLayerStat
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!