Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-13 02:31:36 

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: endLuminosityBlock
0252 // Description: This is will be run at the end of each luminosity block
0253 // Inputs:
0254 // * const LuminosityBlock& lumiSeg = Luminosity Block information
0255 // * const EventSetup&      context = Event Setup information
0256 //____________________________________________________________________________
0257 void L1TOccupancyClient::dqmEndLuminosityBlock(DQMStore::IBooker& ibooker,
0258                                                DQMStore::IGetter& igetter,
0259                                                const edm::LuminosityBlock& lumiSeg,
0260                                                const edm::EventSetup& c) {
0261   book(ibooker, igetter);
0262 
0263   int eventLS = lumiSeg.id().luminosityBlock();
0264 
0265   if (verbose_) {
0266     cout << "[L1TOccupancyClient:] Called endLuminosityBlock()" << endl;
0267     cout << "[L1TOccupancyClient:] Lumisection: " << eventLS << endl;
0268   }
0269 
0270   // Loop over every test in python
0271   for (std::vector<ParameterSet*>::const_iterator it = mValidTests.begin(); it != mValidTests.end(); it++) {
0272     ParameterSet& test = (**it);
0273     string algo_name = test.getUntrackedParameter<string>("algoName", "XYSymmetry");
0274     string test_name = test.getParameter<string>("testName");
0275 
0276     if (verbose_) {
0277       cout << "[L1TOccupancyClient:] Starting calculations for " << algo_name << " on:" << test_name << endl;
0278     }
0279 
0280     if (algo_name == "XYSymmetry") {
0281       ParameterSet ps = (**it).getParameter<ParameterSet>("algoParams");
0282       string histPath = ps.getParameter<string>("histPath");
0283 
0284       vector<pair<int, double> > deadChannels;
0285       vector<pair<int, double> > statDev;
0286       bool enoughStats = false;
0287 
0288       // Update histo's data with data of this LS
0289       hservice_->updateHistogramEndLS(igetter, test_name, histPath, eventLS);
0290 
0291       // Perform the test
0292       double dead = xySymmetry(ps, test_name, deadChannels, statDev, enoughStats);
0293       stringstream str;
0294       str << test_name << "_cumu_LS_" << eventLS;
0295 
0296       if (verbose_) {
0297         TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0298         cumulative_save->SetTitle(str.str().c_str());
0299         TDirectory* td = file_->GetDirectory(test_name.c_str());
0300         td->cd(string(test_name + "_Histos_AllLS").c_str());
0301         cumulative_save->Write();
0302       }
0303 
0304       // If we have enough statistics, we can write test result
0305       if (enoughStats) {
0306         // Make the result histogram
0307         printDeadChannels(deadChannels, meResults[test_name]->getTH2F(), statDev, test_name);
0308 
0309         if (verbose_) {
0310           TH2F* cumulative_save = (TH2F*)hservice_->getDifferentialHistogram(test_name)->Clone(str.str().c_str());
0311           cumulative_save->SetTitle(str.str().c_str());
0312           TDirectory* td = file_->GetDirectory(("DQM_L1TOccupancyClient_Snapshots_LS.root:/" + test_name).c_str());
0313           td->cd(string(test_name + "_Histos").c_str());
0314           cumulative_save->Write();
0315 
0316           // save the result histo
0317           TH2F* h2f = meResults[test_name]->getTH2F();
0318           stringstream str2;
0319           str2 << test_name << "_result_LS_" << eventLS;
0320           TH2F* dead_save = (TH2F*)h2f->Clone(str2.str().c_str());
0321 
0322           td->cd(string(test_name + "_Results").c_str());
0323           dead_save->SetTitle(str2.str().c_str());
0324           dead_save->Write();
0325         }
0326 
0327         // Updating test results
0328         meDifferential[test_name]->Reset();
0329         meDifferential[test_name]->getTH2F()->Add(hservice_->getDifferentialHistogram(test_name));
0330 
0331         vector<int> lsCertification = hservice_->getLSCertification(test_name);
0332 
0333         // Fill fraction of dead channels
0334         for (unsigned int i = 0; i < lsCertification.size(); i++) {
0335           int bin = meCertification[test_name]->getTH1()->FindBin(lsCertification[i]);
0336           meCertification[test_name]->getTH1()->SetBinContent(bin, 1 - dead);
0337         }
0338 
0339         // Reset differential histo
0340         hservice_->resetHisto(test_name);
0341 
0342         if (verbose_) {
0343           cout << "Now we have enough statstics for " << test_name << endl;
0344         }
0345 
0346       } else {
0347         if (verbose_) {
0348           cout << "we don't have enough statstics for " << test_name << endl;
0349         }
0350       }
0351     } else {
0352       if (verbose_) {
0353         cout << "No valid algorithm" << std::endl;
0354       }
0355     }
0356   }
0357 }
0358 
0359 //____________________________________________________________________________
0360 // Function: analyze
0361 // Description: This is will be run for every event
0362 // Inputs:
0363 // * const Event&      e       = Event information
0364 // * const EventSetup& context = Event Setup information
0365 //____________________________________________________________________________
0366 //void L1TOccupancyClient::analyze(const Event& e, const EventSetup& context){}
0367 
0368 //____________________________________________________________________________
0369 // Function: xySymmetry
0370 // Description: This method preforms the XY Symmetry test
0371 // Inputs:
0372 // * ParameterSet                     ps           = Parameters for the test
0373 // * std::string                      test_name    = Test name of the test to be executed
0374 // * std::vector< pair<int,double> >& deadChannels = Vector of
0375 // * std::vector< pair<int,double> >& statDev      =
0376 // * bool&                            enoughStats  =
0377 // Outputs:
0378 // * double = fraction of bins that failed test, DeadChannels in vector, in: ParameterSet of test parameters
0379 //____________________________________________________________________________
0380 double L1TOccupancyClient::xySymmetry(const ParameterSet& ps,
0381                                       string iTestName,
0382                                       vector<pair<int, double> >& deadChannels,
0383                                       vector<pair<int, double> >& statDev,
0384                                       bool& enoughStats) {
0385   // Getting differential histogram for this this thes
0386   TH2F* diffHist = hservice_->getDifferentialHistogram(iTestName);
0387 
0388   int pAxis = ps.getUntrackedParameter<int>("axis", 1);
0389   int pAverageMode = ps.getUntrackedParameter<int>("averageMode", 2);  // 1=arith. mean, 2=median
0390   int nBinsX = diffHist->GetNbinsX();                                  // actual number of bins x
0391   int nBinsY = diffHist->GetNbinsY();                                  // actual number of bins y
0392 
0393   // Axis==1 : Means symmetry axis is vertical
0394   if (pAxis == 1) {
0395     int maxBinStrip, centralBinStrip;  // x-coordinate of strips
0396 
0397     maxBinStrip = nBinsX;
0398 
0399     // If takeCenter=true  determine central bin of the pAxis
0400     // If takeCenter=false determine the bin to use based user input
0401     if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
0402       centralBinStrip = nBinsX / 2 + 1;
0403     } else {
0404       double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
0405       getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 1);
0406     }
0407 
0408     // Assuming odd number of strips --> first comparison is middle strip to itself
0409     int upBinStrip = centralBinStrip;
0410     int lowBinStrip = centralBinStrip;
0411 
0412     // If even number decrease lowBinstrip by one
0413     if (nBinsX % 2 == 0) {
0414       lowBinStrip--;
0415     }
0416 
0417     // Do we have enough statistics? Min(Max(strip_i,strip_j))>threshold
0418     std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
0419 
0420     int nActualStrips = 0;  //number of strips that are not fully masked
0421     for (int j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; j++, k--) {
0422       // Protection for when both strips are masked
0423       if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0424         maxAvgs[nActualStrips] = TMath::Max(getAvrg(diffHist, iTestName, pAxis, nBinsY, j, pAverageMode),
0425                                             getAvrg(diffHist, iTestName, pAxis, nBinsY, k, pAverageMode));
0426         nActualStrips++;
0427       }
0428     }
0429 
0430     vector<double> defaultMu0up;
0431     defaultMu0up.push_back(13.7655);
0432     defaultMu0up.push_back(184.742);
0433     defaultMu0up.push_back(50735.3);
0434     defaultMu0up.push_back(-97.6793);
0435 
0436     TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0437     vector<double> params = ps.getUntrackedParameter<vector<double> >("params_mu0_up", defaultMu0up);
0438     for (unsigned int i = 0; i < params.size(); i++) {
0439       tf.SetParameter(i, params[i]);
0440     }
0441     int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0442 
0443     vector<double> defaultMu0low;
0444     defaultMu0low.push_back(2.19664);
0445     defaultMu0low.push_back(1.94546);
0446     defaultMu0low.push_back(-99.3263);
0447     defaultMu0low.push_back(19.388);
0448 
0449     params = ps.getUntrackedParameter<vector<double> >("params_mu0_low", defaultMu0low);
0450     for (unsigned int i = 0; i < params.size(); i++) {
0451       tf.SetParameter(i, params[i]);
0452     }
0453     int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0454 
0455     if (verbose_) {
0456       cout << "nbins: " << hservice_->getNBinsHistogram(iTestName) << endl;
0457       cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0458     }
0459 
0460     if (nActualStrips > 0) {
0461       enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0462     } else if (verbose_) {
0463       cout << "No valid strips found, insufficient statistics." << endl;
0464     }
0465     if (verbose_) {
0466       cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0467            << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0468            << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0469     }
0470 
0471     //if enough statistics
0472     //make the test
0473     if (enoughStats) {
0474       for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0475         double avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, upBinStrip, pAverageMode);
0476         compareWithStrip(
0477             diffHist, iTestName, lowBinStrip, nBinsY, pAxis, avg, ps, deadChannels);  //compare with lower side
0478 
0479         avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, lowBinStrip, pAverageMode);
0480         compareWithStrip(
0481             diffHist, iTestName, upBinStrip, nBinsY, pAxis, avg, ps, deadChannels);  //compare with upper side
0482       }
0483     }
0484   }
0485 
0486   // pAxis==2 : Means symetry pAxis is horizontal
0487   else if (pAxis == 2) {
0488     int maxBinStrip, centralBinStrip;  //x-coordinate of strips
0489 
0490     maxBinStrip = nBinsY;
0491 
0492     // Determine center of diagram: either with set pAxis or middle of diagram
0493     if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
0494       centralBinStrip = nBinsY / 2 + 1;
0495     } else {
0496       double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
0497       getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 2);
0498     }
0499 
0500     //assuming odd number of strips --> first comparison is middle strip to itself
0501     int lowBinStrip = centralBinStrip, upBinStrip = centralBinStrip;
0502 
0503     //if even number
0504     if (nBinsX % 2 == 0) {
0505       //decrease lowBinstrip by one
0506       lowBinStrip--;
0507     }
0508 
0509     //do we have enough statistics? Min(Max(strip_i,strip_j))>threshold
0510     std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
0511     int nActualStrips = 0;
0512     for (int j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; j++, k--) {
0513       if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0514         maxAvgs[nActualStrips] = TMath::Max(getAvrg(diffHist, iTestName, pAxis, nBinsX, j, pAverageMode),
0515                                             getAvrg(diffHist, iTestName, pAxis, nBinsX, k, pAverageMode));
0516         nActualStrips++;
0517       }
0518     }
0519 
0520     vector<double> defaultMu0up;
0521     defaultMu0up.push_back(13.7655);
0522     defaultMu0up.push_back(184.742);
0523     defaultMu0up.push_back(50735.3);
0524     defaultMu0up.push_back(-97.6793);
0525 
0526     vector<double> params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_up", defaultMu0up);
0527     TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0528     for (unsigned int i = 0; i < params.size(); i++) {
0529       tf.SetParameter(i, params[i]);
0530     }
0531     int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0532 
0533     vector<double> defaultMu0low;
0534     defaultMu0low.push_back(2.19664);
0535     defaultMu0low.push_back(1.94546);
0536     defaultMu0low.push_back(-99.3263);
0537     defaultMu0low.push_back(19.388);
0538 
0539     params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_low", defaultMu0low);
0540     for (unsigned int i = 0; i < params.size(); i++) {
0541       tf.SetParameter(i, params[i]);
0542     }
0543     int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0544     if (verbose_) {
0545       cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0546     }
0547     if (nActualStrips > 0) {
0548       enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0549     } else if (verbose_) {
0550       cout << "No valid strips found, insufficient statistics." << endl;
0551     }
0552     if (verbose_) {
0553       cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0554            << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0555            << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0556     }
0557 
0558     //if we have enough statistics
0559     //make the test
0560     if (enoughStats) {
0561       for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0562         double avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, upBinStrip, pAverageMode);
0563         compareWithStrip(
0564             diffHist, iTestName, lowBinStrip, nBinsX, pAxis, avg, ps, deadChannels);  //compare with lower side
0565 
0566         avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, lowBinStrip, pAverageMode);
0567         compareWithStrip(
0568             diffHist, iTestName, upBinStrip, nBinsX, pAxis, avg, ps, deadChannels);  //compare with upper side
0569       }
0570     }
0571   } else {
0572     if (verbose_) {
0573       cout << "Invalid axis" << endl;
0574     }
0575   }
0576 
0577   return (deadChannels.size() - hservice_->getNBinsMasked(iTestName)) * 1.0 / hservice_->getNBinsHistogram(iTestName);
0578 }
0579 
0580 //____________________________________________________________________________
0581 // Function: getAvrg
0582 // Description: Calculate strip average with method iAvgMode, where strip is
0583 // prependicular to iAxis at bin iBinStrip of histogram iHist
0584 // Inputs:
0585 // * TH2F*  iHist     = Histogram to be tested
0586 // * string iTestName = Name of the test
0587 // * int    iAxis     = Axis prependicular to plot symmetry
0588 // * int    iNBins     = Number of bins in the strip
0589 // * int    iBinStrip = Bin corresponding to the strip in iAxis
0590 // * int    iAvgMode  = Type of average mode 1) Average 2) Median
0591 // Outputs:
0592 // * double = Average of input strip
0593 //____________________________________________________________________________
0594 double L1TOccupancyClient::getAvrg(TH2F* iHist, string iTestName, int iAxis, int iNBins, int iBinStrip, int iAvgMode) {
0595   double avg = 0.0;
0596   TH1D* proj = nullptr;
0597   TH2F* histo = new TH2F(*iHist);
0598 
0599   std::vector<double> values;
0600   int marked;
0601 
0602   if (iAxis == 1) {
0603     switch (iAvgMode) {
0604       // arithmetic average
0605       case 1:
0606         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0607         proj = histo->ProjectionX();
0608         avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0609         break;
0610 
0611       // median
0612       case 2:
0613         hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0614         proj = histo->ProjectionY("_py", iBinStrip, iBinStrip);
0615         for (int i = 0; i < iNBins; i++) {
0616           values.push_back(proj->GetBinContent(i + 1));
0617         }
0618         avg = TMath::Median(iNBins, &values[0]);
0619         break;
0620       default:
0621         if (verbose_) {
0622           cout << "Invalid averaging mode!" << endl;
0623         }
0624         break;
0625     }
0626   } else if (iAxis == 2) {
0627     switch (iAvgMode) {
0628       // arithmetic average
0629       case 1:
0630         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0631         proj = histo->ProjectionY();
0632         avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0633         break;
0634       // median
0635       case 2:
0636         hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0637         proj = histo->ProjectionX("_px", iBinStrip, iBinStrip);
0638         for (int i = 0; i < iNBins; i++) {
0639           values.push_back(proj->GetBinContent(i + 1));
0640         }
0641 
0642         avg = TMath::Median(iNBins, &values[0]);
0643         break;
0644       default:
0645         if (verbose_) {
0646           cout << "invalid averaging mode!" << endl;
0647         }
0648         break;
0649     }
0650   } else {
0651     if (verbose_) {
0652       cout << "invalid axis" << endl;
0653     }
0654   }
0655   delete proj;
0656   delete histo;
0657   return avg;
0658 }
0659 
0660 //____________________________________________________________________________
0661 // Function: printDeadChannels
0662 // Description:
0663 // Inputs:
0664 // * vector< pair<int,double> > iDeadChannels     = List of bin that are masked of failed tthe test
0665 // * TH2F*                      oHistDeadChannels = Histogram where test results should be printed
0666 // * vector< pair<int,double> > statDev           = ???
0667 // * string                     iTestName         = Name of the test
0668 //____________________________________________________________________________
0669 void L1TOccupancyClient::printDeadChannels(const vector<pair<int, double> >& iDeadChannels,
0670                                            TH2F* oHistDeadChannels,
0671                                            const vector<std::pair<int, double> >& statDev,
0672                                            string iTestName) {
0673   // Reset the dead channels histogram
0674   oHistDeadChannels->Reset();
0675   if (verbose_) {
0676     cout << "suspect or masked channels of " << iTestName << ": ";
0677   }
0678 
0679   int x, y, z;
0680 
0681   // put all bad (value=1) and masked (value=-1) cells in histo
0682   for (std::vector<pair<int, double> >::const_iterator it = iDeadChannels.begin(); it != iDeadChannels.end(); it++) {
0683     int bin = (*it).first;
0684     oHistDeadChannels->GetBinXYZ(bin, x, y, z);
0685 
0686     if (hservice_->isMasked(iTestName, x, y)) {
0687       oHistDeadChannels->SetBinContent(bin, -1);
0688       if (verbose_) {
0689         printf("(%4i,%4i) Masked\n", x, y);
0690       }
0691     } else {
0692       oHistDeadChannels->SetBinContent(bin, 1);
0693       if (verbose_) {
0694         printf("(%4i,%4i) Failed test\n", x, y);
0695       }
0696     }
0697   }
0698 
0699   if (verbose_) {
0700     cout << "total number of suspect channels: " << (iDeadChannels.size() - (hservice_->getNBinsMasked(iTestName)))
0701          << endl;
0702   }
0703 }
0704 
0705 //____________________________________________________________________________
0706 // Function: compareWithStrip
0707 // Description: Evaluates statistical compatibility of a strip (cell by cell) against a given average
0708 // Inputs:
0709 // * TH2F*                      iHist      = Histogram to be tested
0710 // * string                     iTestName  = Which test to apply
0711 // * int                        iBinStrip  = Bin Coordinate (in bin units) of the stripo
0712 // * int                        iNBins     = Number of Bins in the strip
0713 // * int                        iAxis      = Which Axis is prependicular to the plot symmetry.
0714 // * double                     iAvg       = Average of the strip
0715 // * ParameterSet               iPS        = Parameters for the test
0716 // * vector<pair<int,double> >& oChannels  = Output of bin that are masked or failed the test
0717 // Outputs:
0718 // * int = Number of dead channels
0719 //____________________________________________________________________________
0720 int L1TOccupancyClient::compareWithStrip(TH2F* iHist,
0721                                          string iTestName,
0722                                          int iBinStrip,
0723                                          int iNBins,
0724                                          int iAxis,
0725                                          double iAvg,
0726                                          const ParameterSet& iPS,
0727                                          vector<pair<int, double> >& oChannels) {
0728   int dead = 0;
0729 
0730   //
0731   if (iAxis == 1) {
0732     // Get and set parameters for working curves
0733     TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0734     TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0735     fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0736     fmuup->SetParameter(1, iAvg);
0737     fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0738     fmulow->SetParameter(1, iAvg);
0739 
0740     TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0741 
0742     // Evaluate sigma up
0743     vector<double> defaultChi2up;
0744     defaultChi2up.push_back(5.45058e-05);
0745     defaultChi2up.push_back(0.268756);
0746     defaultChi2up.push_back(-11.7515);
0747 
0748     vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0749     for (unsigned int i = 0; i < params.size(); i++) {
0750       fchi->SetParameter(i, params[i]);
0751     }
0752     double sigma_up = fchi->Eval(iAvg);
0753 
0754     // Evaluate sigma low
0755     vector<double> defaultChi2low;
0756     defaultChi2low.push_back(4.11095e-05);
0757     defaultChi2low.push_back(0.577451);
0758     defaultChi2low.push_back(-10.378);
0759 
0760     params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0761     for (unsigned int i = 0; i < params.size(); i++) {
0762       fchi->SetParameter(i, params[i]);
0763     }
0764     double sigma_low = fchi->Eval(iAvg);
0765 
0766     if (verbose_) {
0767       cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0768     }
0769 
0770     for (int i = 1; i <= iNBins; i++) {
0771       if (verbose_) {
0772         cout << "    " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(iBinStrip, i))
0773              << " low: " << fmulow->Eval(iHist->GetBinContent(iBinStrip, i)) << endl;
0774       }
0775 
0776       // Evaluate chi2 for cells
0777       double muup = fmuup->Eval(iHist->GetBinContent(iBinStrip, i));
0778       double mulow = fmulow->Eval(iHist->GetBinContent(iBinStrip, i));
0779 
0780       // If channel is masked -> set it to value -1
0781       if (hservice_->isMasked(iTestName, iBinStrip, i)) {
0782         oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0783       }
0784       //else perform test
0785       else if (muup > sigma_up || mulow > sigma_low ||
0786                ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0787                 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0788         dead++;
0789         oChannels.push_back(
0790             pair<int, double>(iHist->GetBin(iBinStrip, i), abs(iHist->GetBinContent(iBinStrip, i) - iAvg) / iAvg));
0791       }
0792     }
0793   }
0794   //
0795   else if (iAxis == 2) {
0796     //get and set parameters for working curves
0797     TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0798     TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0799     fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0800     fmuup->SetParameter(1, iAvg);
0801     fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0802     fmulow->SetParameter(1, iAvg);
0803 
0804     TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0805 
0806     // Evaluate sigma up
0807     vector<double> defaultChi2up;
0808     defaultChi2up.push_back(5.45058e-05);
0809     defaultChi2up.push_back(0.268756);
0810     defaultChi2up.push_back(-11.7515);
0811 
0812     vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0813     for (unsigned int i = 0; i < params.size(); i++) {
0814       fchi->SetParameter(i, params[i]);
0815     }
0816     double sigma_up = fchi->Eval(iAvg);
0817 
0818     // Evaluate sigma low
0819     vector<double> defaultChi2low;
0820     defaultChi2low.push_back(4.11095e-05);
0821     defaultChi2low.push_back(0.577451);
0822     defaultChi2low.push_back(-10.378);
0823 
0824     params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0825     for (unsigned int i = 0; i < params.size(); i++) {
0826       fchi->SetParameter(i, params[i]);
0827     }
0828     double sigma_low = fchi->Eval(iAvg);
0829 
0830     if (verbose_) {
0831       cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0832     }
0833 
0834     for (int i = 1; i <= iNBins; i++) {
0835       if (verbose_) {
0836         cout << "    " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(i, iBinStrip))
0837              << " low: " << fmulow->Eval(iHist->GetBinContent(i, iBinStrip)) << endl;
0838       }
0839 
0840       //evaluate chi2 for cells
0841       double muup = fmuup->Eval(iHist->GetBinContent(i, iBinStrip));
0842       double mulow = fmulow->Eval(iHist->GetBinContent(i, iBinStrip));
0843 
0844       //if channel is masked -> set it to value -1
0845       if (hservice_->isMasked(iTestName, i, iBinStrip)) {
0846         oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0847       }
0848       //else perform test
0849       else if (muup > sigma_up || mulow > sigma_low ||
0850                ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0851                 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0852         dead++;
0853         oChannels.push_back(
0854             pair<int, double>(iHist->GetBin(i, iBinStrip), abs(iHist->GetBinContent(i, iBinStrip) - iAvg) / iAvg));
0855       }
0856     }
0857   } else {
0858     if (verbose_) {
0859       cout << "invalid axis" << endl;
0860     }
0861   }
0862 
0863   return dead;
0864 }
0865 
0866 //____________________________________________________________________________
0867 // Function: getBinCoordinateOnAxisWithValue
0868 // Description: Returns the bin global bin number with the iValue in the iAxis
0869 // Inputs:
0870 // * TH2F*  iHist          = Histogram to be tested
0871 // * double iValue         = Value to be evaluated in the histogram iHist
0872 // * int&   oBinCoordinate = (output) bin number (X or Y) for iValue
0873 // * int    iAxis          = Axis to be used
0874 //____________________________________________________________________________
0875 void L1TOccupancyClient::getBinCoordinateOnAxisWithValue(TH2F* iHist, double iValue, int& oBinCoordinate, int iAxis) {
0876   int nBinsX = iHist->GetNbinsX();  //actual number of bins x
0877   int nBinsY = iHist->GetNbinsY();  //actual number of bins y
0878 
0879   if (iAxis == 1) {
0880     int global = iHist->GetXaxis()->FindFixBin(iValue);
0881 
0882     // If parameter exceeds axis' value: set to maximum number of bins in x-axis
0883     if (global > nBinsX * nBinsY) {
0884       global = iHist->GetXaxis()->GetLast();
0885     }
0886 
0887     // Get coordinates of bin
0888     int y, z;
0889     iHist->GetBinXYZ(global, oBinCoordinate, y, z);
0890   } else if (iAxis == 2) {
0891     int global = iHist->GetYaxis()->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->GetYaxis()->GetLast();
0896     }
0897 
0898     // Get coordinates of bin
0899     int x, z;
0900     iHist->GetBinXYZ(global, x, oBinCoordinate, z);
0901   }
0902 }
0903 
0904 //define this as a plug-in
0905 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!