Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:02:19

0001 #include "Utilities/Testing/interface/CppUnit_testdriver.icpp"
0002 #include "cppunit/extensions/HelperMacros.h"
0003 
0004 #include "FWCore/Utilities/interface/FileInPath.h"
0005 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0006 #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
0007 
0008 #include <iostream>
0009 #include <fstream>
0010 #include <vector>
0011 #include <string>
0012 #include <iomanip>
0013 #include <cstdlib>
0014 
0015 #include "TBenchmark.h"
0016 
0017 using std::cout;
0018 using std::endl;
0019 using std::flush;
0020 using std::ofstream;
0021 using std::setw;
0022 using std::string;
0023 using std::vector;
0024 
0025 class testJetCorrectorParameters : public CppUnit::TestFixture {
0026   CPPUNIT_TEST_SUITE(testJetCorrectorParameters);
0027   CPPUNIT_TEST(benchmarkBinIndex1D);
0028   CPPUNIT_TEST(benchmarkBinIndex3D);
0029   CPPUNIT_TEST(compareBinIndex1D);
0030   CPPUNIT_TEST(compareBinIndex3D);
0031   CPPUNIT_TEST_SUITE_END();
0032 
0033 public:
0034   testJetCorrectorParameters() {}
0035   ~testJetCorrectorParameters() {}
0036   void setUp();
0037   void tearDown();
0038   void setupCorrector(bool is3D);
0039   void destroyCorrector();
0040   void generateFiles();
0041 
0042   void benchmarkBinIndex(bool is3D);
0043   void benchmarkBinIndex1D() { benchmarkBinIndex(false); }
0044   void benchmarkBinIndex3D() { benchmarkBinIndex(true); }
0045   void compareBinIndex1D();
0046   void compareBinIndex3D();
0047 
0048   inline void loadbar3(
0049       unsigned int x, unsigned int n, unsigned int w = 50, unsigned int freq = 100, string prefix = "") {
0050     if ((x != n) && (x % (n / freq) != 0))
0051       return;
0052     float ratio = x / (float)n;
0053     int c = ratio * w;
0054 
0055     cout << prefix << std::fixed << setw(8) << std::setprecision(0) << (ratio * 100) << "% [";
0056     for (int x = 0; x < c; x++)
0057       cout << "=";
0058     for (unsigned int x = c; x < w; x++)
0059       cout << " ";
0060     cout << "] (" << x << "/" << n << ")\r" << flush;
0061   }
0062 
0063 private:
0064   string filename;
0065   TBenchmark* m_benchmark;
0066   JetCorrectorParameters* L1JetPar;
0067   vector<JetCorrectorParameters> vPar;
0068   FactorizedJetCorrector* jetCorrector;
0069   vector<float> fX;
0070   vector<float> veta;
0071   vector<float> vrho;
0072   vector<float> vpt;
0073 };
0074 
0075 ///registration of the test so that the runner can find it
0076 CPPUNIT_TEST_SUITE_REGISTRATION(testJetCorrectorParameters);
0077 
0078 void testJetCorrectorParameters::setUp() {
0079   m_benchmark = new TBenchmark();
0080   veta = {-5.191, -4.889, -4.716, -4.538, -4.363, -4.191, -4.013, -3.839, -3.664, -3.489, -3.314, -3.139,
0081           -2.964, -2.853, -2.65,  -2.5,   -2.322, -2.172, -2.043, -1.93,  -1.83,  -1.74,  -1.653, -1.566,
0082           -1.479, -1.392, -1.305, -1.218, -1.131, -1.044, -0.957, -0.879, -0.783, -0.696, -0.609, -0.522,
0083           -0.435, -0.348, -0.261, -0.174, -0.087, 0,      0.087,  0.174,  0.261,  0.348,  0.435,  0.522,
0084           0.609,  0.696,  0.783,  0.879,  0.957,  1.044,  1.131,  1.218,  1.305,  1.392,  1.479,  1.566,
0085           1.653,  1.74,   1.83,   1.93,   2.043,  2.172,  2.322,  2.5,    2.65,   2.853,  2.964,  3.139,
0086           3.314,  3.489,  3.664,  3.839,  4.013,  4.191,  4.363,  4.538,  4.716,  4.889,  5.191};
0087   for (unsigned int irho = 1; irho <= 50; irho++) {
0088     vrho.push_back(irho);
0089   }
0090   vpt = {1,   2,   3,   4,    5,    6,    7,    8,    9,    10,   11,   12,   13,   14,   15,
0091          17,  20,  23,  27,   30,   35,   40,   45,   57,   72,   90,   120,  150,  200,  300,
0092          400, 550, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500};
0093 
0094   generateFiles();
0095 }
0096 
0097 void testJetCorrectorParameters::tearDown() {
0098   if (m_benchmark != nullptr)
0099     delete m_benchmark;
0100 }
0101 
0102 void testJetCorrectorParameters::setupCorrector(bool is3D) {
0103   string path = "CondFormats/JetMETObjects/data/";
0104   try {
0105     string to_try = path + ((is3D) ? "testJetCorrectorParameters_3D_L1FastJet_AK4PFchs.txt"
0106                                    : "testJetCorrectorParameters_1D_L1FastJet_AK4PFchs.txt");
0107     edm::FileInPath strFIP(to_try);
0108     filename = strFIP.fullPath();
0109   } catch (edm::Exception const& ex) {
0110     throw ex;
0111   }
0112 
0113   L1JetPar = new JetCorrectorParameters(filename);
0114   vPar.push_back(*L1JetPar);
0115   jetCorrector = new FactorizedJetCorrector(vPar);
0116 }
0117 
0118 void testJetCorrectorParameters::destroyCorrector() {
0119   delete L1JetPar;
0120   delete jetCorrector;
0121 }
0122 
0123 void testJetCorrectorParameters::generateFiles() {
0124   string path = "CondFormats/JetMETObjects/data/";
0125   string name1D = "testJetCorrectorParameters_1D_L1FastJet_AK4PFchs.txt";
0126   string name3D = "testJetCorrectorParameters_3D_L1FastJet_AK4PFchs.txt";
0127 
0128   ofstream txtFile1D;
0129   txtFile1D.open((path + name1D).c_str());
0130   CPPUNIT_ASSERT(txtFile1D.is_open());
0131   txtFile1D << "{1 JetEta 3 Rho JetPt JetA 1 Correction L1FastJet}" << endl;
0132   for (unsigned int ieta = 0; ieta < veta.size() - 1; ieta++) {
0133     txtFile1D << std::setw(15) << veta[ieta] << std::setw(15) << veta[ieta + 1] << std::setw(15) << "6" << std::setw(15)
0134               << "0" << std::setw(15) << "200" << std::setw(15) << "1" << std::setw(15) << "6500" << std::setw(15)
0135               << "0" << std::setw(15) << "10" << std::setw(15) << endl;
0136   }
0137   txtFile1D.close();
0138 
0139   vector<float> ptbins = {4,  8,  10,  13,  16,  19,  22,  25,  30,  35,   40,  45,
0140                           60, 75, 100, 160, 230, 330, 440, 600, 820, 1100, 1400};
0141   ofstream txtFile3D;
0142   txtFile3D.open((path + name3D).c_str());
0143   CPPUNIT_ASSERT(txtFile3D.is_open());
0144   txtFile3D << "{3 JetEta Rho JetPt 1 JetPt 1 Correction L1FastJet}" << endl;
0145   for (unsigned int ieta = 0; ieta < veta.size() - 1; ieta++) {
0146     for (unsigned int irho = 0; irho < vrho.size() - 1; irho++) {
0147       for (unsigned int ipt = 0; ipt < ptbins.size() - 1; ipt++) {
0148         txtFile3D << std::setw(15) << veta[ieta] << std::setw(15) << veta[ieta + 1] << std::setw(15) << vrho[irho]
0149                   << std::setw(15);
0150         if (vrho[irho + 1] == 50)
0151           txtFile3D << "200" << std::setw(15);
0152         else
0153           txtFile3D << vrho[irho + 1] << std::setw(15);
0154         txtFile3D << ptbins[ipt] << std::setw(15);
0155         if (ipt + 1 == ptbins.size() - 1)
0156           txtFile3D << "6500" << std::setw(15);
0157         else
0158           txtFile3D << ptbins[ipt + 1] << std::setw(15);
0159         txtFile3D << "2" << std::setw(15) << ptbins[ipt] << std::setw(15) << ptbins[ipt + 1] << std::setw(15) << endl;
0160       }
0161     }
0162   }
0163   txtFile3D.close();
0164 }
0165 
0166 void testJetCorrectorParameters::benchmarkBinIndex(bool is3D) {
0167   float oldCPU = 0, newCPU = 0, oldReal = 0, newReal = 0;
0168   unsigned int ntests = (is3D) ? 100000 : 1000000;
0169   if (is3D)
0170     fX = {5.0, 50.0, 100.0};
0171   else
0172     fX = {5.0};
0173   setupCorrector(is3D);
0174 
0175   cout << endl << "testJetCorrectorParameters::benchmarkBinIndex NTests = " << ntests << endl;
0176   cout << "testJetCorrectorParameters::benchmarkBinIndex benchmarking binIndex with file " << filename << " ... "
0177        << flush;
0178   m_benchmark->Reset();
0179   m_benchmark->Start("event");
0180   for (unsigned int i = 0; i < ntests; i++) {
0181     L1JetPar->binIndex(fX);
0182   }
0183   m_benchmark->Stop("event");
0184   cout << "DONE" << endl
0185        << "testJetCorrectorParameters::benchmarkBinIndex" << endl
0186        << "\tCPU time = " << m_benchmark->GetCpuTime("event") / double(ntests) << " s" << endl
0187        << "\tReal time = " << m_benchmark->GetRealTime("event") / double(ntests) << " s" << endl;
0188   oldCPU = m_benchmark->GetCpuTime("event") / double(ntests);
0189   oldReal = m_benchmark->GetRealTime("event") / double(ntests);
0190 
0191   cout << "testJetCorrectorParameters::benchmarkBinIndex benchmarking binIndexN with file " << filename << " ... "
0192        << flush;
0193   m_benchmark->Reset();
0194   m_benchmark->Start("event");
0195   for (unsigned int i = 0; i < ntests; i++) {
0196     L1JetPar->binIndexN(fX);
0197   }
0198   m_benchmark->Stop("event");
0199   cout << "DONE" << endl
0200        << "testJetCorrectorParameters::benchmarkBinIndex" << endl
0201        << "\tCPU time = " << m_benchmark->GetCpuTime("event") / double(ntests) << " s" << endl
0202        << "\tReal time = " << m_benchmark->GetRealTime("event") / double(ntests) << " s" << endl;
0203   newCPU = m_benchmark->GetCpuTime("event") / double(ntests);
0204   newReal = m_benchmark->GetRealTime("event") / double(ntests);
0205 
0206   cout << "testJetCorrectorParameters::benchmarkBinIndex" << endl
0207        << "\tCPU speedup = " << oldCPU / newCPU << endl
0208        << "\tReal speedup = " << oldReal / newReal << endl;
0209 
0210   CPPUNIT_ASSERT(oldCPU / newCPU >= 1.0);
0211   //CPPUNIT_ASSERT(oldReal/newReal >= 1.0); //this might fail not due to longer L1JetPar->binIndexN(fX) execution
0212   if (oldReal / newReal >= 1.0)
0213     cout << "newReal value increased oldReal, which might be due to system load" << endl;
0214 
0215   destroyCorrector();
0216 }
0217 
0218 void testJetCorrectorParameters::compareBinIndex1D() {
0219   fX = {0.0};
0220   float eta = -9999;
0221   int oldBin = -1, newBin = -1;
0222   setupCorrector(false);
0223 
0224   cout << endl << "testJetCorrectorParameters::compareBinIndex1D" << endl;
0225   for (unsigned int ieta = 0; ieta < veta.size() - 1; ieta++) {
0226     loadbar3(ieta + 1, veta.size() - 1, 50, 10, "\tProgress:");
0227     eta = (veta[ieta] + veta[ieta + 1]) / 2.0;
0228     fX = {eta};
0229     oldBin = L1JetPar->binIndex(fX);
0230     newBin = L1JetPar->binIndexN(fX);
0231     if ((oldBin < 0 && newBin >= 0) || (oldBin >= 0 && newBin < 0)) {
0232       cout << "ERROR::testJetCorrectorParameters::compareBinIndex1D Unable to find the right bin for (eta)=(" << eta
0233            << ")" << endl
0234            << "\t(oldBin,newBin)=(" << oldBin << "," << newBin << ")" << endl;
0235       CPPUNIT_ASSERT(oldBin >= 0 && newBin >= 0);
0236     } else if (oldBin != newBin) {
0237       cout << "ERROR::testJetCorrectorParameters::compareBinIndex1D oldBin!=newBin (" << oldBin << "!=" << newBin
0238            << ") for (eta)=(" << eta << ")" << endl;
0239       CPPUNIT_ASSERT(oldBin == newBin);
0240     }
0241 
0242     jetCorrector->setJetEta(eta);
0243     jetCorrector->setRho(50.0);
0244     jetCorrector->setJetA(0.5);
0245     jetCorrector->setJetPt(100.0);
0246     CPPUNIT_ASSERT(jetCorrector->getCorrection() >= 0);
0247   }
0248   destroyCorrector();
0249   cout << endl
0250        << "testJetCorrectorParameters::compareBinIndex1D All bins match between the linear and non-linear search "
0251           "algorithms for 1D files."
0252        << endl;
0253 }
0254 
0255 void testJetCorrectorParameters::compareBinIndex3D() {
0256   fX = {0.0, 0.0, 0.0};
0257   float eta = -9999, rho = -9999, pt = -9999;
0258   int oldBin = -1, newBin = -1;
0259   setupCorrector(true);
0260 
0261   cout << endl << "testJetCorrectorParameters::compareBinIndex3D" << endl;
0262   for (unsigned int ieta = 0; ieta < veta.size() - 1; ieta++) {
0263     for (unsigned int irho = 0; irho < vrho.size() - 1; irho++) {
0264       for (unsigned int ipt = 0; ipt < vpt.size() - 1; ipt++) {
0265         loadbar3(ieta * ((vrho.size() - 1) * (vpt.size() - 1)) + irho * (vpt.size() - 1) + ipt + 1,
0266                  (veta.size() - 1) * (vrho.size() - 1) * (vpt.size() - 1),
0267                  50,
0268                  100,
0269                  "\tProgress");
0270         eta = (veta[ieta] + veta[ieta + 1]) / 2.0;
0271         rho = (vrho[irho] + vrho[irho + 1]) / 2.0;
0272         pt = (vpt[ipt] + vpt[ipt + 1]) / 2.0;
0273         fX = {eta, rho, pt};
0274         oldBin = L1JetPar->binIndex(fX);
0275         newBin = L1JetPar->binIndexN(fX);
0276         if ((oldBin < 0 && newBin >= 0) || (oldBin >= 0 && newBin < 0)) {
0277           cout << "ERROR::testJetCorrectorParameters::compareBinIndex3D Unable to find the right bin for (eta,rho,pt)=("
0278                << eta << "," << rho << "," << pt << ")" << endl
0279                << "\t(oldBin,newBin)=(" << oldBin << "," << newBin << ")" << endl;
0280           CPPUNIT_ASSERT(oldBin >= 0 && newBin >= 0);
0281         } else if (oldBin != newBin) {
0282           cout << "ERROR::testJetCorrectorParameters::compareBinIndex3D oldBin!=newBin (" << oldBin << "!=" << newBin
0283                << ") for (eta,rho,pt)=(" << eta << "," << rho << "," << pt << ")" << endl;
0284           CPPUNIT_ASSERT(oldBin == newBin);
0285         }
0286 
0287         jetCorrector->setJetEta(eta);
0288         jetCorrector->setRho(rho);
0289         jetCorrector->setJetA(0.5);
0290         jetCorrector->setJetPt(pt);
0291         CPPUNIT_ASSERT(jetCorrector->getCorrection() >= 0);
0292       }
0293     }
0294   }
0295   destroyCorrector();
0296   cout << endl
0297        << "testJetCorrectorParameters::compareBinIndex3D All bins match between the linear and non-linear search "
0298           "algorithms for 3D files."
0299        << endl;
0300 }