Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:08:00 

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 i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
0422       double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsY, j, pAverageMode);
0423       double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsY, k, pAverageMode);
0424 
0425       // Protection for when both strips are masked
0426       if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0427         maxAvgs[i] = TMath::Max(avg1, avg2);
0428         nActualStrips++;
0429       }
0430     }
0431 
0432     vector<double> defaultMu0up;
0433     defaultMu0up.push_back(13.7655);
0434     defaultMu0up.push_back(184.742);
0435     defaultMu0up.push_back(50735.3);
0436     defaultMu0up.push_back(-97.6793);
0437 
0438     TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0439     vector<double> params = ps.getUntrackedParameter<vector<double> >("params_mu0_up", defaultMu0up);
0440     for (unsigned int i = 0; i < params.size(); i++) {
0441       tf.SetParameter(i, params[i]);
0442     }
0443     int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0444 
0445     vector<double> defaultMu0low;
0446     defaultMu0low.push_back(2.19664);
0447     defaultMu0low.push_back(1.94546);
0448     defaultMu0low.push_back(-99.3263);
0449     defaultMu0low.push_back(19.388);
0450 
0451     params = ps.getUntrackedParameter<vector<double> >("params_mu0_low", defaultMu0low);
0452     for (unsigned int i = 0; i < params.size(); i++) {
0453       tf.SetParameter(i, params[i]);
0454     }
0455     int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0456 
0457     if (verbose_) {
0458       cout << "nbins: " << hservice_->getNBinsHistogram(iTestName) << endl;
0459       cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0460     }
0461 
0462     enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0463     if (verbose_) {
0464       cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0465            << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0466            << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0467     }
0468 
0469     //if enough statistics
0470     //make the test
0471     if (enoughStats) {
0472       for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0473         double avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, upBinStrip, pAverageMode);
0474         compareWithStrip(
0475             diffHist, iTestName, lowBinStrip, nBinsY, pAxis, avg, ps, deadChannels);  //compare with lower side
0476 
0477         avg = getAvrg(diffHist, iTestName, pAxis, nBinsY, lowBinStrip, pAverageMode);
0478         compareWithStrip(
0479             diffHist, iTestName, upBinStrip, nBinsY, pAxis, avg, ps, deadChannels);  //compare with upper side
0480       }
0481     }
0482   }
0483 
0484   // pAxis==2 : Means symetry pAxis is horizontal
0485   else if (pAxis == 2) {
0486     int maxBinStrip, centralBinStrip;  //x-coordinate of strips
0487 
0488     maxBinStrip = nBinsY;
0489 
0490     // Determine center of diagram: either with set pAxis or middle of diagram
0491     if (ps.getUntrackedParameter<bool>("takeCenter", true)) {
0492       centralBinStrip = nBinsY / 2 + 1;
0493     } else {
0494       double pAxisSymmetryValue = ps.getParameter<double>("axisSymmetryValue");
0495       getBinCoordinateOnAxisWithValue(diffHist, pAxisSymmetryValue, centralBinStrip, 2);
0496     }
0497 
0498     //assuming odd number of strips --> first comparison is middle strip to itself
0499     int lowBinStrip = centralBinStrip, upBinStrip = centralBinStrip;
0500 
0501     //if even number
0502     if (nBinsX % 2 == 0) {
0503       //decrease lowBinstrip by one
0504       lowBinStrip--;
0505     }
0506 
0507     //do we have enough statistics? Min(Max(strip_i,strip_j))>threshold
0508     std::unique_ptr<double[]> maxAvgs(new double[maxBinStrip - upBinStrip + 1]);
0509     int nActualStrips = 0;
0510     for (int i = 0, j = upBinStrip, k = lowBinStrip; j <= maxBinStrip; i++, j++, k--) {
0511       double avg1 = getAvrg(diffHist, iTestName, pAxis, nBinsX, j, pAverageMode);
0512       double avg2 = getAvrg(diffHist, iTestName, pAxis, nBinsX, k, pAverageMode);
0513       if (!hservice_->isStripMasked(iTestName, j, pAxis) && !hservice_->isStripMasked(iTestName, k, pAxis)) {
0514         maxAvgs[i] = TMath::Max(avg1, avg2);
0515         nActualStrips++;
0516       }
0517     }
0518 
0519     vector<double> defaultMu0up;
0520     defaultMu0up.push_back(13.7655);
0521     defaultMu0up.push_back(184.742);
0522     defaultMu0up.push_back(50735.3);
0523     defaultMu0up.push_back(-97.6793);
0524 
0525     vector<double> params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_up", defaultMu0up);
0526     TF1 tf("myFunc", "[0]*(TMath::Log(x*[1]+[2]))+[3]", 10., 11000.);
0527     for (unsigned int i = 0; i < params.size(); i++) {
0528       tf.SetParameter(i, params[i]);
0529     }
0530     int statsup = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0531 
0532     vector<double> defaultMu0low;
0533     defaultMu0low.push_back(2.19664);
0534     defaultMu0low.push_back(1.94546);
0535     defaultMu0low.push_back(-99.3263);
0536     defaultMu0low.push_back(19.388);
0537 
0538     params = ps.getUntrackedParameter<std::vector<double> >("params_mu0_low", defaultMu0low);
0539     for (unsigned int i = 0; i < params.size(); i++) {
0540       tf.SetParameter(i, params[i]);
0541     }
0542     int statslow = (int)tf.Eval(hservice_->getNBinsHistogram(iTestName));
0543     if (verbose_) {
0544       cout << "statsup= " << statsup << ", statslow= " << statslow << endl;
0545     }
0546     enoughStats = TMath::MinElement(nActualStrips, maxAvgs.get()) > TMath::Max(statsup, statslow);
0547     if (verbose_) {
0548       cout << "stats: " << TMath::MinElement(nActualStrips, maxAvgs.get())
0549            << ", statsAvg: " << diffHist->GetEntries() / hservice_->getNBinsHistogram(iTestName)
0550            << ", threshold: " << TMath::Max(statsup, statslow) << endl;
0551     }
0552 
0553     //if we have enough statistics
0554     //make the test
0555     if (enoughStats) {
0556       for (; upBinStrip <= maxBinStrip; upBinStrip++, lowBinStrip--) {
0557         double avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, upBinStrip, pAverageMode);
0558         compareWithStrip(
0559             diffHist, iTestName, lowBinStrip, nBinsX, pAxis, avg, ps, deadChannels);  //compare with lower side
0560 
0561         avg = getAvrg(diffHist, iTestName, pAxis, nBinsX, lowBinStrip, pAverageMode);
0562         compareWithStrip(
0563             diffHist, iTestName, upBinStrip, nBinsX, pAxis, avg, ps, deadChannels);  //compare with upper side
0564       }
0565     }
0566   } else {
0567     if (verbose_) {
0568       cout << "Invalid axis" << endl;
0569     }
0570   }
0571 
0572   return (deadChannels.size() - hservice_->getNBinsMasked(iTestName)) * 1.0 / hservice_->getNBinsHistogram(iTestName);
0573 }
0574 
0575 //____________________________________________________________________________
0576 // Function: getAvrg
0577 // Description: Calculate strip average with method iAvgMode, where strip is
0578 // prependicular to iAxis at bin iBinStrip of histogram iHist
0579 // Inputs:
0580 // * TH2F*  iHist     = Histogram to be tested
0581 // * string iTestName = Name of the test
0582 // * int    iAxis     = Axis prependicular to plot symmetry
0583 // * int    iNBins     = Number of bins in the strip
0584 // * int    iBinStrip = Bin corresponding to the strip in iAxis
0585 // * int    iAvgMode  = Type of average mode 1) Average 2) Median
0586 // Outputs:
0587 // * double = Average of input strip
0588 //____________________________________________________________________________
0589 double L1TOccupancyClient::getAvrg(TH2F* iHist, string iTestName, int iAxis, int iNBins, int iBinStrip, int iAvgMode) {
0590   double avg = 0.0;
0591   TH1D* proj = nullptr;
0592   TH2F* histo = new TH2F(*iHist);
0593 
0594   std::vector<double> values;
0595   int marked;
0596 
0597   if (iAxis == 1) {
0598     switch (iAvgMode) {
0599       // arithmetic average
0600       case 1:
0601         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0602         proj = histo->ProjectionX();
0603         avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0604         break;
0605 
0606       // median
0607       case 2:
0608         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0609         proj = histo->ProjectionY("_py", iBinStrip, iBinStrip);
0610         for (int i = 0; i < iNBins; i++) {
0611           values.push_back(proj->GetBinContent(i + 1));
0612         }
0613         avg = TMath::Median(iNBins, &values[0]);
0614         break;
0615       default:
0616         if (verbose_) {
0617           cout << "Invalid averaging mode!" << endl;
0618         }
0619         break;
0620     }
0621   } else if (iAxis == 2) {
0622     switch (iAvgMode) {
0623       // arithmetic average
0624       case 1:
0625         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0626         proj = histo->ProjectionY();
0627         avg = proj->GetBinContent(iBinStrip) / (iNBins - marked);
0628         break;
0629       // median
0630       case 2:
0631         marked = hservice_->maskBins(iTestName, histo, iBinStrip, iAxis);
0632         proj = histo->ProjectionX("_px", iBinStrip, iBinStrip);
0633         for (int i = 0; i < iNBins; i++) {
0634           values.push_back(proj->GetBinContent(i + 1));
0635         }
0636 
0637         avg = TMath::Median(iNBins, &values[0]);
0638         break;
0639       default:
0640         if (verbose_) {
0641           cout << "invalid averaging mode!" << endl;
0642         }
0643         break;
0644     }
0645   } else {
0646     if (verbose_) {
0647       cout << "invalid axis" << endl;
0648     }
0649   }
0650   delete proj;
0651   delete histo;
0652   return avg;
0653 }
0654 
0655 //____________________________________________________________________________
0656 // Function: printDeadChannels
0657 // Description:
0658 // Inputs:
0659 // * vector< pair<int,double> > iDeadChannels     = List of bin that are masked of failed tthe test
0660 // * TH2F*                      oHistDeadChannels = Histogram where test results should be printed
0661 // * vector< pair<int,double> > statDev           = ???
0662 // * string                     iTestName         = Name of the test
0663 //____________________________________________________________________________
0664 void L1TOccupancyClient::printDeadChannels(const vector<pair<int, double> >& iDeadChannels,
0665                                            TH2F* oHistDeadChannels,
0666                                            const vector<std::pair<int, double> >& statDev,
0667                                            string iTestName) {
0668   // Reset the dead channels histogram
0669   oHistDeadChannels->Reset();
0670   if (verbose_) {
0671     cout << "suspect or masked channels of " << iTestName << ": ";
0672   }
0673 
0674   int x, y, z;
0675 
0676   // put all bad (value=1) and masked (value=-1) cells in histo
0677   for (std::vector<pair<int, double> >::const_iterator it = iDeadChannels.begin(); it != iDeadChannels.end(); it++) {
0678     int bin = (*it).first;
0679     oHistDeadChannels->GetBinXYZ(bin, x, y, z);
0680 
0681     if (hservice_->isMasked(iTestName, x, y)) {
0682       oHistDeadChannels->SetBinContent(bin, -1);
0683       if (verbose_) {
0684         printf("(%4i,%4i) Masked\n", x, y);
0685       }
0686     } else {
0687       oHistDeadChannels->SetBinContent(bin, 1);
0688       if (verbose_) {
0689         printf("(%4i,%4i) Failed test\n", x, y);
0690       }
0691     }
0692   }
0693 
0694   if (verbose_) {
0695     cout << "total number of suspect channels: " << (iDeadChannels.size() - (hservice_->getNBinsMasked(iTestName)))
0696          << endl;
0697   }
0698 }
0699 
0700 //____________________________________________________________________________
0701 // Function: compareWithStrip
0702 // Description: Evaluates statistical compatibility of a strip (cell by cell) against a given average
0703 // Inputs:
0704 // * TH2F*                      iHist      = Histogram to be tested
0705 // * string                     iTestName  = Which test to apply
0706 // * int                        iBinStrip  = Bin Coordinate (in bin units) of the stripo
0707 // * int                        iNBins     = Number of Bins in the strip
0708 // * int                        iAxis      = Which Axis is prependicular to the plot symmetry.
0709 // * double                     iAvg       = Average of the strip
0710 // * ParameterSet               iPS        = Parameters for the test
0711 // * vector<pair<int,double> >& oChannels  = Output of bin that are masked or failed the test
0712 // Outputs:
0713 // * int = Number of dead channels
0714 //____________________________________________________________________________
0715 int L1TOccupancyClient::compareWithStrip(TH2F* iHist,
0716                                          string iTestName,
0717                                          int iBinStrip,
0718                                          int iNBins,
0719                                          int iAxis,
0720                                          double iAvg,
0721                                          const ParameterSet& iPS,
0722                                          vector<pair<int, double> >& oChannels) {
0723   int dead = 0;
0724 
0725   //
0726   if (iAxis == 1) {
0727     // Get and set parameters for working curves
0728     TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0729     TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0730     fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0731     fmuup->SetParameter(1, iAvg);
0732     fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0733     fmulow->SetParameter(1, iAvg);
0734 
0735     TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0736 
0737     // Evaluate sigma up
0738     vector<double> defaultChi2up;
0739     defaultChi2up.push_back(5.45058e-05);
0740     defaultChi2up.push_back(0.268756);
0741     defaultChi2up.push_back(-11.7515);
0742 
0743     vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0744     for (unsigned int i = 0; i < params.size(); i++) {
0745       fchi->SetParameter(i, params[i]);
0746     }
0747     double sigma_up = fchi->Eval(iAvg);
0748 
0749     // Evaluate sigma low
0750     vector<double> defaultChi2low;
0751     defaultChi2low.push_back(4.11095e-05);
0752     defaultChi2low.push_back(0.577451);
0753     defaultChi2low.push_back(-10.378);
0754 
0755     params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0756     for (unsigned int i = 0; i < params.size(); i++) {
0757       fchi->SetParameter(i, params[i]);
0758     }
0759     double sigma_low = fchi->Eval(iAvg);
0760 
0761     if (verbose_) {
0762       cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0763     }
0764 
0765     for (int i = 1; i <= iNBins; i++) {
0766       if (verbose_) {
0767         cout << "    " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(iBinStrip, i))
0768              << " low: " << fmulow->Eval(iHist->GetBinContent(iBinStrip, i)) << endl;
0769       }
0770 
0771       // Evaluate chi2 for cells
0772       double muup = fmuup->Eval(iHist->GetBinContent(iBinStrip, i));
0773       double mulow = fmulow->Eval(iHist->GetBinContent(iBinStrip, i));
0774 
0775       // If channel is masked -> set it to value -1
0776       if (hservice_->isMasked(iTestName, iBinStrip, i)) {
0777         oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0778       }
0779       //else perform test
0780       else if (muup > sigma_up || mulow > sigma_low ||
0781                ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0782                 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0783         dead++;
0784         oChannels.push_back(
0785             pair<int, double>(iHist->GetBin(iBinStrip, i), abs(iHist->GetBinContent(iBinStrip, i) - iAvg) / iAvg));
0786       }
0787     }
0788   }
0789   //
0790   else if (iAxis == 2) {
0791     //get and set parameters for working curves
0792     TF1* fmuup = new TF1("fmuup", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0793     TF1* fmulow = new TF1("fmulow", "TMath::Log(TMath::PoissonI(x,[0])/TMath::PoissonI(x,[1]))", -10000., 10000.);
0794     fmuup->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorup", 2.0));
0795     fmuup->SetParameter(1, iAvg);
0796     fmulow->SetParameter(0, iAvg * iPS.getUntrackedParameter<double>("factorlow", 0.1));
0797     fmulow->SetParameter(1, iAvg);
0798 
0799     TF1* fchi = new TF1("fchi", "[0]*x**2+[1]*x+[2]", 0., 1500.);
0800 
0801     // Evaluate sigma up
0802     vector<double> defaultChi2up;
0803     defaultChi2up.push_back(5.45058e-05);
0804     defaultChi2up.push_back(0.268756);
0805     defaultChi2up.push_back(-11.7515);
0806 
0807     vector<double> params = iPS.getUntrackedParameter<vector<double> >("params_chi2_up", defaultChi2up);
0808     for (unsigned int i = 0; i < params.size(); i++) {
0809       fchi->SetParameter(i, params[i]);
0810     }
0811     double sigma_up = fchi->Eval(iAvg);
0812 
0813     // Evaluate sigma low
0814     vector<double> defaultChi2low;
0815     defaultChi2low.push_back(4.11095e-05);
0816     defaultChi2low.push_back(0.577451);
0817     defaultChi2low.push_back(-10.378);
0818 
0819     params = iPS.getUntrackedParameter<vector<double> >("params_chi2_low", defaultChi2low);
0820     for (unsigned int i = 0; i < params.size(); i++) {
0821       fchi->SetParameter(i, params[i]);
0822     }
0823     double sigma_low = fchi->Eval(iAvg);
0824 
0825     if (verbose_) {
0826       cout << "binstrip= " << iBinStrip << ", sigmaup= " << sigma_up << ", sigmalow= " << sigma_low << endl;
0827     }
0828 
0829     for (int i = 1; i <= iNBins; i++) {
0830       if (verbose_) {
0831         cout << "    " << i << " binContent: up:" << fmuup->Eval(iHist->GetBinContent(i, iBinStrip))
0832              << " low: " << fmulow->Eval(iHist->GetBinContent(i, iBinStrip)) << endl;
0833       }
0834 
0835       //evaluate chi2 for cells
0836       double muup = fmuup->Eval(iHist->GetBinContent(i, iBinStrip));
0837       double mulow = fmulow->Eval(iHist->GetBinContent(i, iBinStrip));
0838 
0839       //if channel is masked -> set it to value -1
0840       if (hservice_->isMasked(iTestName, i, iBinStrip)) {
0841         oChannels.push_back(pair<int, double>(iHist->GetBin(iBinStrip, i), -1.0));
0842       }
0843       //else perform test
0844       else if (muup > sigma_up || mulow > sigma_low ||
0845                ((fabs(muup) == std::numeric_limits<double>::infinity()) &&
0846                 (fabs(mulow) == std::numeric_limits<double>::infinity()))) {
0847         dead++;
0848         oChannels.push_back(
0849             pair<int, double>(iHist->GetBin(i, iBinStrip), abs(iHist->GetBinContent(i, iBinStrip) - iAvg) / iAvg));
0850       }
0851     }
0852   } else {
0853     if (verbose_) {
0854       cout << "invalid axis" << endl;
0855     }
0856   }
0857 
0858   return dead;
0859 }
0860 
0861 //____________________________________________________________________________
0862 // Function: getBinCoordinateOnAxisWithValue
0863 // Description: Returns the bin global bin number with the iValue in the iAxis
0864 // Inputs:
0865 // * TH2F*  iHist          = Histogram to be tested
0866 // * double iValue         = Value to be evaluated in the histogram iHist
0867 // * int&   oBinCoordinate = (output) bin number (X or Y) for iValue
0868 // * int    iAxis          = Axis to be used
0869 //____________________________________________________________________________
0870 void L1TOccupancyClient::getBinCoordinateOnAxisWithValue(TH2F* iHist, double iValue, int& oBinCoordinate, int iAxis) {
0871   int nBinsX = iHist->GetNbinsX();  //actual number of bins x
0872   int nBinsY = iHist->GetNbinsY();  //actual number of bins y
0873 
0874   if (iAxis == 1) {
0875     int global = iHist->GetXaxis()->FindFixBin(iValue);
0876 
0877     // If parameter exceeds axis' value: set to maximum number of bins in x-axis
0878     if (global > nBinsX * nBinsY) {
0879       global = iHist->GetXaxis()->GetLast();
0880     }
0881 
0882     // Get coordinates of bin
0883     int y, z;
0884     iHist->GetBinXYZ(global, oBinCoordinate, y, z);
0885   } else if (iAxis == 2) {
0886     int global = iHist->GetYaxis()->FindFixBin(iValue);
0887 
0888     // If parameter exceeds axis' value: set to maximum number of bins in x-axis
0889     if (global > nBinsX * nBinsY) {
0890       global = iHist->GetYaxis()->GetLast();
0891     }
0892 
0893     // Get coordinates of bin
0894     int x, z;
0895     iHist->GetBinXYZ(global, x, oBinCoordinate, z);
0896   }
0897 }
0898 
0899 //define this as a plug-in
0900 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!