Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:10:50 

0001 #include "DQM/L1TMonitorClient/interface/L1TOccupancyClient.h"
0002 
0003 #include "FWCore/ServiceRegistry/interface/Service.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include <cstdio>
0010 #include <sstream>
0011 #include <cmath>
0012 #include <vector>
0013 #include <TMath.h>
0014 #include <climits>
0015 #include <TFile.h>
0016 #include <TDirectory.h>
0017 #include <TProfile.h>
0018 
0019 using namespace std;
0020 using namespace edm;
0021 
0022 //____________________________________________________________________________
0023 // Function: L1TOccupancyClient
0024 // Description: This is the constructor, basic variable initialization
0025 // Inputs:
0026 // * const edm::ParameterSet& ps = Parameter for this analyzer
0027 //____________________________________________________________________________
0028 L1TOccupancyClient::L1TOccupancyClient(const edm::ParameterSet& ps) {
0029   // Get parameters
0030   parameters_ = ps;
0031   verbose_ = ps.getParameter<bool>("verbose");
0032   tests_ = ps.getParameter<std::vector<ParameterSet> >("testParams");
0033 
0034   if (verbose_) {
0035     cout << "[L1TOccupancyClient:] Called constructor" << endl;
0036   }
0037 }
0038 
0039 //____________________________________________________________________________
0040 // Function: ~L1TOccupancyClient
0041 // Description: This is the destructor, basic variable deletion
0042 //____________________________________________________________________________
0043 L1TOccupancyClient::~L1TOccupancyClient() {
0044   if (verbose_) {
0045     cout << "[L1TOccupancyClient:] Called destructor" << endl;
0046   }
0047 }
0048 
0049 //____________________________________________________________________________
0050 // Function: beginRun
0051 // Description: This is will be run at the begining of each run
0052 // Inputs:
0053 // * const Run&        r       = Run information
0054 // * const EventSetup& context = Event Setup information
0055 //____________________________________________________________________________
0056 void L1TOccupancyClient::book(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter) {
0057   hservice_ = new L1TOccupancyClientHistogramService(parameters_, ibooker, verbose_);
0058 
0059   if (verbose_) {
0060     cout << "[L1TOccupancyClient:] Called beginRun" << endl;
0061 
0062     // In verbose mode we will produce an extra output file with several tests
0063     file_ = TFile::Open("DQM_L1TOccupancyClient_Snapshots_LS.root", "RECREATE");
0064   }
0065 
0066   ibooker.setCurrentFolder("L1T/L1TOccupancy/");
0067   //dbe_->setCurrentFolder("L1T/L1TOccupancy/Results");
0068   //dbe_->setCurrentFolder("L1T/L1TOccupancy/BadCellValues");
0069   //dbe_->setCurrentFolder("L1T/L1TOccupancy/Certification");
0070 
0071   // Loop over all tests in defined
0072   for (vector<ParameterSet>::iterator it = tests_.begin(); it != tests_.end(); it++) {
0073     // If the test algorithm is XYSymmetry we create the necessary histograms
0074     if ((*it).getUntrackedParameter<string>("algoName", "XYSymmetry") == "XYSymmetry") {
0075       // Getting Parameters for the test
0076       string testName = (*it).getParameter<string>("testName");
0077       ParameterSet algoParameters = (*it).getParameter<ParameterSet>("algoParams");
0078       string histPath = algoParameters.getParameter<string>("histPath");
0079 
0080       if (verbose_) {
0081         cout << "[L1TOccupancyClient:] Monitored histogram path: " << histPath << endl;
0082 
0083         // Creating verbose file directory structure
0084         // test_name/test_name_Results,
0085         // test_name/test_name_Histos
0086         // TDirectory *td  = file_->mkdir(testName.c_str()             ,testName.c_str());
0087         //FIXME: sub never used gcc361 warning
0088         //TDirectory *sub = td   ->mkdir((testName+"_Results").c_str(),string("_Results").c_str());
0089 
0090         //sub = td->mkdir((testName+"_Histos").c_str()      ,(testName+"_Histos").c_str());
0091         //sub = td->mkdir((testName+"_Histos_AllLS").c_str(),(testName+"_Histos_AllLS").c_str());
0092       }
0093 
0094       // Load histograms in service instance
0095       if (hservice_->loadHisto(igetter, testName, histPath)) {
0096         // Mask channels specified in python file
0097         hservice_->setMaskedBins(testName, algoParameters.getParameter<vector<ParameterSet> >("maskedAreas"));
0098 
0099         // Book MonitorElements
0100         // * Test results
0101         ibooker.setCurrentFolder("L1T/L1TOccupancy/Results");
0102         string title = testName;
0103         MonitorElement* m = ibooker.book2D(title.c_str(), hservice_->getDifferentialHistogram(testName));
0104         m->setTitle(title);
0105         m->Reset();
0106         meResults[title] = m;
0107 
0108         // * Which cells are masked as bad
0109         ibooker.setCurrentFolder("L1T/L1TOccupancy/HistogramDiff");
0110         title = testName;
0111         m = ibooker.book2D(title.c_str(), hservice_->getDifferentialHistogram(testName));
0112         m->Reset();
0113         m->setTitle(title);
0114         meDifferential[title] = m;
0115 
0116         // * Fraction of bad cells
0117         ibooker.setCurrentFolder("L1T/L1TOccupancy/Certification");
0118         title = testName;
0119         m = ibooker.book1D(title.c_str(), title.c_str(), 2500, -.5, 2500. - .5);
0120         m->setTitle(title);
0121         meCertification[title] = m;
0122 
0123         mValidTests.push_back(&(*it));
0124       }
0125     }
0126   }
0127 }
0128 
0129 //____________________________________________________________________________
0130 // Function: endRun
0131 // Description: This is will be run at the end of each run
0132 // Inputs:
0133 // * const Run&        r       = Run information
0134 // * const EventSetup& context = Event Setup information
0135 //____________________________________________________________________________
0136 void L1TOccupancyClient::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter) {
0137   book(ibooker, igetter);
0138 
0139   if (verbose_) {
0140     cout << "[L1TOccupancyClient:] Called endRun()" << endl;
0141   }
0142 
0143   // Loop over every test in python
0144   for (std::vector<ParameterSet*>::iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
0145     ParameterSet& test = (**it);
0146     string algo_name = test.getUntrackedParameter<string>("algoName", "XYSymmetry");
0147     string test_name = test.getParameter<string>("testName");
0148 
0149     if (verbose_) {
0150       cout << "[L1TOccupancyClient:] Starting calculations for: " << algo_name << " on: " << test_name << endl;
0151     }
0152 
0153     if (algo_name == "XYSymmetry") {
0154       ParameterSet ps = (**it).getParameter<ParameterSet>("algoParams");
0155       string histPath = ps.getParameter<string>("histPath");
0156 
0157       vector<pair<int, double> > deadChannels;
0158       vector<pair<int, double> > statDev;
0159       bool enoughStats = false;
0160 
0161       // Make final block
0162       hservice_->updateHistogramEndRun(test_name);
0163 
0164       // Perform the test
0165       double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
0166       stringstream str;
0167       str << test_name << "_cumu_LS_EndRun";
0168 
0169       if (verbose_) {
0170         TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0171 
0172         cumulative_save->SetTitle(str.str().c_str());
0173 
0174         TDirectory* td = file_->GetDirectory(test_name.c_str());
0175 
0176         td->cd(string(test_name + "_Histos_AllLS").c_str());
0177 
0178         cumulative_save->Write();
0179       }
0180       // If we have enough statistics, we can write test result
0181       if (enoughStats) {
0182         // Make the result histogram
0183         printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
0184 
0185         if (verbose_) {
0186           TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0187           cumulative_save->SetTitle(str.str().c_str());
0188           TDirectory* td = file_->GetDirectory(("DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
0189           td->cd(string(test_name + "_Histos").c_str());
0190           cumulative_save->Write();
0191 
0192           // save the result histo
0193           TH2F* h2f = meResults[test_name]->getTH2F();
0194           stringstream str2;
0195           str2 << test_name << "_result_LS_EndRun";
0196           TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
0197 
0198           td->cd(string(test_name + "_Results").c_str());
0199           dead_save->SetTitle(str2.str().c_str());
0200           dead_save->Write();
0201         }
0202 
0203         // Updating test results
0204         meDifferential[test_name]->Reset();
0205         meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
0206 
0207         vector<int> lsCertification = hservice_->getLSCertification(test_name);
0208 
0209         // Fill fraction of dead channels
0210         for (unsigned int i = 0; i < lsCertification.size(); i++) {
0211           int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
0212           meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
0213         }
0214 
0215         // Reset differential histo
0216         hservice_->resetHisto(test_name);
0217 
0218         if (verbose_) {
0219           cout << "Now we have enough statstics for " << test_name << endl;
0220         }
0221 
0222       } else {
0223         if (verbose_) {
0224           cout << "we don't have enough statstics for " << test_name << endl;
0225         }
0226 
0227         // Getting LS which this test monitored
0228         vector<int> lsCertification = hservice_->getLSCertification(test_name);
0229 
0230         // Fill fraction of dead channels
0231         for (unsigned int i = 0; i < lsCertification.size(); i++) {
0232           int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
0233           meCertification[test_name]->getTH1()->SetBinContent(bin, -1);
0234         }
0235       }
0236     } else {
0237       if (verbose_) {
0238         cout << "No valid algorithm" << std::endl;
0239       }
0240     }
0241   }
0242 
0243   if (verbose_) {
0244     file_->Close();
0245   }
0246 
0247   delete hservice_;
0248 }
0249 
0250 //____________________________________________________________________________
0251 // Function: beginLuminosityBlock
0252 // Description: This is will be run at the begining of each luminosity block
0253 // Inputs:
0254 // * const LuminosityBlock& lumiSeg = Luminosity Block information
0255 // * const EventSetup&      context = Event Setup information
0256 //____________________________________________________________________________
0257 
0258 //____________________________________________________________________________
0259 // Function: endLuminosityBlock
0260 // Description: This is will be run at the end of each luminosity block
0261 // Inputs:
0262 // * const LuminosityBlock& lumiSeg = Luminosity Block information
0263 // * const EventSetup&      context = Event Setup information
0264 //____________________________________________________________________________
0265 void L1TOccupancyClient::dqmEndLuminosityBlock(DQMStore::IBooker& ibooker,
0266                                                DQMStore::IGetter& igetter,
0267                                                const edm::LuminosityBlock& lumiSeg,
0268                                                const edm::EventSetup& c) {
0269   book(ibooker, igetter);
0270 
0271   int eventLS = lumiSeg.id().luminosityBlock();
0272 
0273   if (verbose_) {
0274     cout << "[L1TOccupancyClient:] Called endLuminosityBlock()" << endl;
0275     cout << "[L1TOccupancyClient:] Lumisection: " << eventLS << endl;
0276   }
0277 
0278   // Loop over every test in python
0279   for (std::vector<ParameterSet*>::const_iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
0280     ParameterSet& test = (**it);
0281     string algo_name = test.getUntrackedParameter<string>("algoName", "XYSymmetry");
0282     string test_name = test.getParameter<string>("testName");
0283 
0284     if (verbose_) {
0285       cout << "[L1TOccupancyClient:] Starting calculations for " << algo_name << " on:" << test_name << endl;
0286     }
0287 
0288     if (algo_name == "XYSymmetry") {
0289       ParameterSet ps = (**it).getParameter<ParameterSet>("algoParams");
0290       string histPath = ps.getParameter<string>("histPath");
0291 
0292       vector<pair<int, double> > deadChannels;
0293       vector<pair<int, double> > statDev;
0294       bool enoughStats = false;
0295 
0296       // Update histo's data with data of this LS
0297       hservice_->updateHistogramEndLS(igetter, test_name, histPath, eventLS);
0298 
0299       // Perform the test
0300       double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
0301       stringstream str;
0302       str << test_name << "_cumu_LS_" << eventLS;
0303 
0304       if (verbose_) {
0305         TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0306         cumulative_save->SetTitle(str.str().c_str());
0307         TDirectory* td = file_->GetDirectory(test_name.c_str());
0308         td->cd(string(test_name + "_Histos_AllLS").c_str());
0309         cumulative_save->Write();
0310       }
0311 
0312       // If we have enough statistics, we can write test result
0313       if (enoughStats) {
0314         // Make the result histogram
0315         printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
0316 
0317         if (verbose_) {
0318           TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0319           cumulative_save->SetTitle(str.str().c_str());
0320           TDirectory* td = file_->GetDirectory(("DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
0321           td->cd(string(test_name + "_Histos").c_str());
0322           cumulative_save->Write();
0323 
0324           // save the result histo
0325           TH2F* h2f = meResults[test_name]->getTH2F();
0326           stringstream str2;
0327           str2 << test_name << "_result_LS_" << eventLS;
0328           TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
0329 
0330           td->cd(string(test_name + "_Results").c_str());
0331           dead_save->SetTitle(str2.str().c_str());
0332           dead_save->Write();
0333         }
0334 
0335         // Updating test results
0336         meDifferential[test_name]->Reset();
0337         meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
0338 
0339         vector<int> lsCertification = hservice_->getLSCertification(test_name);
0340 
0341         // Fill fraction of dead channels
0342         for (unsigned int i = 0; i < lsCertification.size(); i++) {
0343           int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
0344           meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
0345         }
0346 
0347         // Reset differential histo
0348         hservice_->resetHisto(test_name);
0349 
0350         if (verbose_) {
0351           cout << "Now we have enough statstics for " << test_name << endl;
0352         }
0353 
0354       } else {
0355         if (verbose_) {
0356           cout << "we don't have enough statstics for " << test_name << endl;
0357         }
0358       }
0359     } else {
0360       if (verbose_) {
0361         cout << "No valid algorithm" << std::endl;
0362       }
0363     }
0364   }
0365 }
0366 
0367 //____________________________________________________________________________
0368 // Function: analyze
0369 // Description: This is will be run for every event
0370 // Inputs:
0371 // * const Event&      e       = Event information
0372 // * const EventSetup& context = Event Setup information
0373 //____________________________________________________________________________
0374 //void L1TOccupancyClient::analyze(const Event& e, const EventSetup& context){}
0375 
0376 //____________________________________________________________________________
0377 // Function: xySymmetry
0378 // Description: This method preforms the XY Symmetry test
0379 // Inputs:
0380 // * ParameterSet                     ps           = Parameters for the test
0381 // * std::string                      test_name    = Test name of the test to be executed
0382 // * std::vector< pair<int,double> >& deadChannels = Vector of
0383 // * std::vector< pair<int,double> >& statDev      =
0384 // * bool&                            enoughStats  =
0385 // Outputs:
0386 // * double = fraction of bins that failed test, DeadChannels in vector, in: ParameterSet of test parameters
0387 //____________________________________________________________________________
0388 double L1TOccupancyClient::xySymmetry(const ParameterSet& ps,
0389                                       string iTestName,
0390                                       vector<pair<int, double> >& deadChannels,
0391                                       vector<pair<int, double> >& statDev,
0392                                       bool& enoughStats) {
0393   // Getting differential histogram for this this thes
0394   TH2F* diffHist = hservice_->getDifferentialHistogram(iTestName);
0395 
0396   int pAxis = ps.getUntrackedParameter<int>("axis", 1);
0397   int pAverageMode = ps.getUntrackedParameter<int>("averageMode", 2);  // 1=arith. mean, 2=median
0398   int nBinsX = diffHist->GetNbinsX();                                  // actual number of bins x
0399   int nBinsY = diffHist->GetNbinsY();                                  // actual number of bins y
0400 
0401   // Axis==1 : Means symmetry axis is vertical
0402   if (pAxis == 1) {
0403     int maxBinStrip, centralBinStrip;  // x-coordinate of strips
0404 
0405     maxBinStrip = nBinsX;
0406 
0407     // If takeCenter=true  determine central bin of the pAxis
0408     // If takeCenter=false determine the bin to use based user input
0409     if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
0410       centralBinStrip = nBinsX / 2 + 1;
0411     } else {
0412       double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
0413       getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 1);
0414     }
0415 
0416     // Assuming odd number of strips --> first comparison is middle strip to itself
0417     int upBinStrip = centralBinStrip;
0418     int lowBinStrip = centralBinStrip;
0419 
0420     // If even number decrease lowBinstrip by one
0421     if (nBinsX % 2 == 0) {
0422       lowBinStrip--;
0423     }
0424 
0425     // Do we have enough statistics? Min(Max(strip_i,strip_j))>threshold
0426     std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
0427 
0428     int nActualStrips = 0;  //number of strips that are not fully masked
0429     for (int i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
0430       double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsY, j, pAverageMode);
0431       double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsY, k, pAverageMode);
0432 
0433       // Protection for when both strips are masked
0434       if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0435         maxAvgs[i] = TMath::Max(avg1, avg2);
0436         nActualStrips++;
0437       }
0438     }
0439 
0440     vector<double> defaultMu0up;
0441     defaultMu0up.push_back(13.7655);
0442     defaultMu0up.push_back(184.742);
0443     defaultMu0up.push_back(50735.3);
0444     defaultMu0up.push_back(-97.6793);
0445 
0446     TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0447     vector<double> params = ps.getUntrackedParameter<vector<double> >("params_mu0_up", defaultMu0up);
0448     for (unsigned int i = 0; i < params.size(); i++) {
0449       tf.SetParameter(i, params[i]);
0450     }
0451     int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0452 
0453     vector<double> defaultMu0low;
0454     defaultMu0low.push_back(2.19664);
0455     defaultMu0low.push_back(1.94546);
0456     defaultMu0low.push_back(-99.3263);
0457     defaultMu0low.push_back(19.388);
0458 
0459     params = ps.getUntrackedParameter<vector<double> >("params_mu0_low", defaultMu0low);
0460     for (unsigned int i = 0; i < params.size(); i++) {
0461       tf.SetParameter(i, params[i]);
0462     }
0463     int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0464 
0465     if (verbose_) {
0466       cout << "nbins: " << hservice_->getNBinsHistogram(iTestName) << endl;
0467       cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0468     }
0469 
0470     enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0471     if (verbose_) {
0472       cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0473            << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0474            << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0475     }
0476 
0477     //if enough statistics
0478     //make the test
0479     if (enoughStats) {
0480       for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0481         double avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, upBinStrip, pAverageMode);
0482         compareWithStrip(
0483             diffHist, iTestName, lowBinStrip, nBinsY, pAxis, avg, ps, deadChannels);  //compare with lower side
0484 
0485         avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, lowBinStrip, pAverageMode);
0486         compareWithStrip(
0487             diffHist, iTestName, upBinStrip, nBinsY, pAxis, avg, ps, deadChannels);  //compare with upper side
0488       }
0489     }
0490   }
0491 
0492   // pAxis==2 : Means symetry pAxis is horizontal
0493   else if (pAxis == 2) {
0494     int maxBinStrip, centralBinStrip;  //x-coordinate of strips
0495 
0496     maxBinStrip = nBinsY;
0497 
0498     // Determine center of diagram: either with set pAxis or middle of diagram
0499     if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
0500       centralBinStrip = nBinsY / 2 + 1;
0501     } else {
0502       double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
0503       getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 2);
0504     }
0505 
0506     //assuming odd number of strips --> first comparison is middle strip to itself
0507     int lowBinStrip = centralBinStrip, upBinStrip = centralBinStrip;
0508 
0509     //if even number
0510     if (nBinsX % 2 == 0) {
0511       //decrease lowBinstrip by one
0512       lowBinStrip--;
0513     }
0514 
0515     //do we have enough statistics? Min(Max(strip_i,strip_j))>threshold
0516     std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
0517     int nActualStrips = 0;
0518     for (int i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
0519       double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsX, j, pAverageMode);
0520       double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsX, k, pAverageMode);
0521       if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0522         maxAvgs[i] = TMath::Max(avg1, avg2);
0523         nActualStrips++;
0524       }
0525     }
0526 
0527     vector<double> defaultMu0up;
0528     defaultMu0up.push_back(13.7655);
0529     defaultMu0up.push_back(184.742);
0530     defaultMu0up.push_back(50735.3);
0531     defaultMu0up.push_back(-97.6793);
0532 
0533     vector<double> params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_up", defaultMu0up);
0534     TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0535     for (unsigned int i = 0; i < params.size(); i++) {
0536       tf.SetParameter(i, params[i]);
0537     }
0538     int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0539 
0540     vector<double> defaultMu0low;
0541     defaultMu0low.push_back(2.19664);
0542     defaultMu0low.push_back(1.94546);
0543     defaultMu0low.push_back(-99.3263);
0544     defaultMu0low.push_back(19.388);
0545 
0546     params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_low", defaultMu0low);
0547     for (unsigned int i = 0; i < params.size(); i++) {
0548       tf.SetParameter(i, params[i]);
0549     }
0550     int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0551     if (verbose_) {
0552       cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0553     }
0554     enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0555     if (verbose_) {
0556       cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0557            << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0558            << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0559     }
0560 
0561     //if we have enough statistics
0562     //make the test
0563     if (enoughStats) {
0564       for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0565         double avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, upBinStrip, pAverageMode);
0566         compareWithStrip(
0567             diffHist, iTestName, lowBinStrip, nBinsX, pAxis, avg, ps, deadChannels);  //compare with lower side
0568 
0569         avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, lowBinStrip, pAverageMode);
0570         compareWithStrip(
0571             diffHist, iTestName, upBinStrip, nBinsX, pAxis, avg, ps, deadChannels);  //compare with upper side
0572       }
0573     }
0574   } else {
0575     if (verbose_) {
0576       cout << "Invalid axis" << endl;
0577     }
0578   }
0579 
0580   return (deadChannels.size() - hservice_->getNBinsMasked(iTestName)) * 1.0 / hservice_->getNBinsHistogram(iTestName);
0581 }
0582 
0583 //____________________________________________________________________________
0584 // Function: getAvrg
0585 // Description: Calculate strip average with method iAvgMode, where strip is
0586 // prependicular to iAxis at bin iBinStrip of histogram iHist
0587 // Inputs:
0588 // * TH2F*  iHist     = Histogram to be tested
0589 // * string iTestName = Name of the test
0590 // * int    iAxis     = Axis prependicular to plot symmetry
0591 // * int    iNBins     = Number of bins in the strip
0592 // * int    iBinStrip = Bin corresponding to the strip in iAxis
0593 // * int    iAvgMode  = Type of average mode 1) Average 2) Median
0594 // Outputs:
0595 // * double = Average of input strip
0596 //____________________________________________________________________________
0597 double L1TOccupancyClient::getAvrg(TH2F* iHist, string iTestName, int iAxis, int iNBins, int iBinStrip, int iAvgMode) {
0598   double avg = 0.0;
0599   TH1D* proj = nullptr;
0600   TH2F* histo = new TH2F(*iHist);
0601 
0602   std::vector<double> values;
0603   int marked;
0604 
0605   if (iAxis == 1) {
0606     switch (iAvgMode) {
0607       // arithmetic average
0608       case 1:
0609         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0610         proj = histo->ProjectionX();
0611         avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0612         break;
0613 
0614       // median
0615       case 2:
0616         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0617         proj = histo->ProjectionY("_py", iBinStrip, iBinStrip);
0618         for (int i = 0; i < iNBins; i++) {
0619           values.push_back(proj->GetBinContent(i + 1));
0620         }
0621         avg = TMath::Median(iNBins, &values[0]);
0622         break;
0623       default:
0624         if (verbose_) {
0625           cout << "Invalid averaging mode!" << endl;
0626         }
0627         break;
0628     }
0629   } else if (iAxis == 2) {
0630     switch (iAvgMode) {
0631       // arithmetic average
0632       case 1:
0633         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0634         proj = histo->ProjectionY();
0635         avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0636         break;
0637       // median
0638       case 2:
0639         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0640         proj = histo->ProjectionX("_px", iBinStrip, iBinStrip);
0641         for (int i = 0; i < iNBins; i++) {
0642           values.push_back(proj->GetBinContent(i + 1));
0643         }
0644 
0645         avg = TMath::Median(iNBins, &values[0]);
0646         break;
0647       default:
0648         if (verbose_) {
0649           cout << "invalid averaging mode!" << endl;
0650         }
0651         break;
0652     }
0653   } else {
0654     if (verbose_) {
0655       cout << "invalid axis" << endl;
0656     }
0657   }
0658   delete proj;
0659   delete histo;
0660   return avg;
0661 }
0662 
0663 //____________________________________________________________________________
0664 // Function: printDeadChannels
0665 // Description:
0666 // Inputs:
0667 // * vector< pair<int,double> > iDeadChannels     = List of bin that are masked of failed tthe test
0668 // * TH2F*                      oHistDeadChannels = Histogram where test results should be printed
0669 // * vector< pair<int,double> > statDev           = ???
0670 // * string                     iTestName         = Name of the test
0671 //____________________________________________________________________________
0672 void L1TOccupancyClient::printDeadChannels(const vector<pair<int, double> >& iDeadChannels,
0673                                            TH2F* oHistDeadChannels,
0674                                            const vector<std::pair<int, double> >& statDev,
0675                                            string iTestName) {
0676   // Reset the dead channels histogram
0677   oHistDeadChannels->Reset();
0678   if (verbose_) {
0679     cout << "suspect or masked channels of " << iTestName << ": ";
0680   }
0681 
0682   int x, y, z;
0683   float chi2 = 0.0;
0684 
0685   // put all bad (value=1) and masked (value=-1) cells in histo
0686   for (std::vector<pair<int, double> >::const_iterator it = iDeadChannels.begin(); it != iDeadChannels.end(); it++) {
0687     int bin = (*it).first;
0688     oHistDeadChannels->GetBinXYZ(bin, x, y, z);
0689 
0690     if (hservice_->isMasked(iTestName, x, y)) {
0691       oHistDeadChannels->SetBinContent(bin, -1);
0692       if (verbose_) {
0693         printf("(%4i,%4i) Masked\n", x, y);
0694       }
0695     } else {
0696       oHistDeadChannels->SetBinContent(bin, 1);
0697       if (verbose_) {
0698         printf("(%4i,%4i) Failed test\n", x, y);
0699       }
0700     }
0701   }
0702 
0703   // FIXME: Is this needed?
0704   for (std::vector<pair<int, double> >::const_iterator it = statDev.begin(); it != statDev.end(); it++) {
0705     double dev = (*it).second;
0706     chi2 += dev;
0707   }
0708   //put total chi2 in float
0709 
0710   if (verbose_) {
0711     cout << "total number of suspect channels: " << (iDeadChannels.size() - (hservice_->getNBinsMasked(iTestName)))
0712          << endl;
0713   }
0714 }
0715 
0716 //____________________________________________________________________________
0717 // Function: compareWithStrip
0718 // Description: Evaluates statistical compatibility of a strip (cell by cell) against a given average
0719 // Inputs:
0720 // * TH2F*                      iHist      = Histogram to be tested
0721 // * string                     iTestName  = Which test to apply
0722 // * int                        iBinStrip  = Bin Coordinate (in bin units) of the stripo
0723 // * int                        iNBins     = Number of Bins in the strip
0724 // * int                        iAxis      = Which Axis is prependicular to the plot symmetry.
0725 // * double                     iAvg       = Average of the strip
0726 // * ParameterSet               iPS        = Parameters for the test
0727 // * vector<pair<int,double> >& oChannels  = Output of bin that are masked or failed the test
0728 // Outputs:
0729 // * int = Number of dead channels
0730 //____________________________________________________________________________
0731 int L1TOccupancyClient::compareWithStrip(TH2F* iHist,
0732                                          string iTestName,
0733                                          int iBinStrip,
0734                                          int iNBins,
0735                                          int iAxis,
0736                                          double iAvg,
0737                                          const ParameterSet& iPS,
0738                                          vector<pair<int, double> >& oChannels) {
0739   int dead = 0;
0740 
0741   //
0742   if (iAxis == 1) {
0743     // Get and set parameters for working curves
0744     TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0745     TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0746     fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0747     fmuup->SetParameter(1, iAvg);
0748     fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0749     fmulow->SetParameter(1, iAvg);
0750 
0751     TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0752 
0753     // Evaluate sigma up
0754     vector<double> defaultChi2up;
0755     defaultChi2up.push_back(5.45058e-05);
0756     defaultChi2up.push_back(0.268756);
0757     defaultChi2up.push_back(-11.7515);
0758 
0759     vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0760     for (unsigned int i = 0; i < params.size(); i++) {
0761       fchi->SetParameter(i, params[i]);
0762     }
0763     double sigma_up = fchi->Eval(iAvg);
0764 
0765     // Evaluate sigma low
0766     vector<double> defaultChi2low;
0767     defaultChi2low.push_back(4.11095e-05);
0768     defaultChi2low.push_back(0.577451);
0769     defaultChi2low.push_back(-10.378);
0770 
0771     params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0772     for (unsigned int i = 0; i < params.size(); i++) {
0773       fchi->SetParameter(i, params[i]);
0774     }
0775     double sigma_low = fchi->Eval(iAvg);
0776 
0777     if (verbose_) {
0778       cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0779     }
0780 
0781     for (int i = 1; i <= iNBins; i++) {
0782       if (verbose_) {
0783         cout << "    " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(iBinStrip, i))
0784              << " low: " << fmulow->Eval(iHist->GetBinContent(iBinStrip, i)) << endl;
0785       }
0786 
0787       // Evaluate chi2 for cells
0788       double muup = fmuup->Eval(iHist->GetBinContent(iBinStrip, i));
0789       double mulow = fmulow->Eval(iHist->GetBinContent(iBinStrip, i));
0790 
0791       // If channel is masked -> set it to value -1
0792       if (hservice_->isMasked(iTestName, iBinStrip, i)) {
0793         oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0794       }
0795       //else perform test
0796       else if (muup > sigma_up || mulow > sigma_low ||
0797                ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0798                 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0799         dead++;
0800         oChannels.push_back(
0801             pair<int, double>(iHist->GetBin(iBinStrip, i), abs(iHist->GetBinContent(iBinStrip, i) - iAvg) / iAvg));
0802       }
0803     }
0804   }
0805   //
0806   else if (iAxis == 2) {
0807     //get and set parameters for working curves
0808     TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0809     TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0810     fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0811     fmuup->SetParameter(1, iAvg);
0812     fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0813     fmulow->SetParameter(1, iAvg);
0814 
0815     TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0816 
0817     // Evaluate sigma up
0818     vector<double> defaultChi2up;
0819     defaultChi2up.push_back(5.45058e-05);
0820     defaultChi2up.push_back(0.268756);
0821     defaultChi2up.push_back(-11.7515);
0822 
0823     vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0824     for (unsigned int i = 0; i < params.size(); i++) {
0825       fchi->SetParameter(i, params[i]);
0826     }
0827     double sigma_up = fchi->Eval(iAvg);
0828 
0829     // Evaluate sigma low
0830     vector<double> defaultChi2low;
0831     defaultChi2low.push_back(4.11095e-05);
0832     defaultChi2low.push_back(0.577451);
0833     defaultChi2low.push_back(-10.378);
0834 
0835     params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0836     for (unsigned int i = 0; i < params.size(); i++) {
0837       fchi->SetParameter(i, params[i]);
0838     }
0839     double sigma_low = fchi->Eval(iAvg);
0840 
0841     if (verbose_) {
0842       cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0843     }
0844 
0845     for (int i = 1; i <= iNBins; i++) {
0846       if (verbose_) {
0847         cout << "    " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(i, iBinStrip))
0848              << " low: " << fmulow->Eval(iHist->GetBinContent(i, iBinStrip)) << endl;
0849       }
0850 
0851       //evaluate chi2 for cells
0852       double muup = fmuup->Eval(iHist->GetBinContent(i, iBinStrip));
0853       double mulow = fmulow->Eval(iHist->GetBinContent(i, iBinStrip));
0854 
0855       //if channel is masked -> set it to value -1
0856       if (hservice_->isMasked(iTestName, i, iBinStrip)) {
0857         oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0858       }
0859       //else perform test
0860       else if (muup > sigma_up || mulow > sigma_low ||
0861                ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0862                 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0863         dead++;
0864         oChannels.push_back(
0865             pair<int, double>(iHist->GetBin(i, iBinStrip), abs(iHist->GetBinContent(i, iBinStrip) - iAvg) / iAvg));
0866       }
0867     }
0868   } else {
0869     if (verbose_) {
0870       cout << "invalid axis" << endl;
0871     }
0872   }
0873 
0874   return dead;
0875 }
0876 
0877 //____________________________________________________________________________
0878 // Function: getBinCoordinateOnAxisWithValue
0879 // Description: Returns the bin global bin number with the iValue in the iAxis
0880 // Inputs:
0881 // * TH2F*  iHist          = Histogram to be tested
0882 // * double iValue         = Value to be evaluated in the histogram iHist
0883 // * int&   oBinCoordinate = (output) bin number (X or Y) for iValue
0884 // * int    iAxis          = Axis to be used
0885 //____________________________________________________________________________
0886 void L1TOccupancyClient::getBinCoordinateOnAxisWithValue(TH2F* iHist, double iValue, int& oBinCoordinate, int iAxis) {
0887   int nBinsX = iHist->GetNbinsX();  //actual number of bins x
0888   int nBinsY = iHist->GetNbinsY();  //actual number of bins y
0889 
0890   if (iAxis == 1) {
0891     int global = iHist->GetXaxis()->FindFixBin(iValue);
0892 
0893     // If parameter exceeds axis' value: set to maximum number of bins in x-axis
0894     if (global > nBinsX * nBinsY) {
0895       global = iHist->GetXaxis()->GetLast();
0896     }
0897 
0898     // Get coordinates of bin
0899     int y, z;
0900     iHist->GetBinXYZ(global, oBinCoordinate, y, z);
0901   } else if (iAxis == 2) {
0902     int global = iHist->GetYaxis()->FindFixBin(iValue);
0903 
0904     // If parameter exceeds axis' value: set to maximum number of bins in x-axis
0905     if (global > nBinsX * nBinsY) {
0906       global = iHist->GetYaxis()->GetLast();
0907     }
0908 
0909     // Get coordinates of bin
0910     int x, z;
0911     iHist->GetBinXYZ(global, x, oBinCoordinate, z);
0912   }
0913 }
0914 
0915 //define this as a plug-in
0916 DEFINE_FWK_MODULE(L1TOccupancyClient);
[ 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!