Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "DQM/SiStripCommissioningAnalysis/interface/PedsFullNoiseAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/PedsFullNoiseAnalysis.h"
0003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0004 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 
0007 #include "TProfile.h"
0008 #include "TH1.h"
0009 #include "TH2.h"
0010 #include "TF1.h"
0011 #include "TFitResult.h"
0012 #include "TMath.h"
0013 #include "Math/DistFunc.h"
0014 #include "Math/ProbFuncMathCore.h"
0015 #include "Fit/BinData.h"
0016 #include "HFitInterface.h"
0017 #include "Math/GoFTest.h"
0018 #include <iostream>
0019 #include <iomanip>
0020 #include <cmath>
0021 
0022 using namespace sistrip;
0023 using namespace std;
0024 // ----------------------------------------------------------------------------
0025 //
0026 PedsFullNoiseAlgorithm::PedsFullNoiseAlgorithm(const edm::ParameterSet& pset, PedsFullNoiseAnalysis* const anal)
0027     : CommissioningAlgorithm(anal),
0028       hPeds_(nullptr, ""),
0029       hNoise_(nullptr, ""),
0030       hNoise2D_(nullptr, ""),
0031       maxDriftResidualCut_(pset.getParameter<double>("MaxDriftResidualCut")),
0032       minStripNoiseCut_(pset.getParameter<double>("MinStripNoiseCut")),
0033       maxStripNoiseCut_(pset.getParameter<double>("MaxStripNoiseCut")),
0034       maxStripNoiseSignificanceCut_(pset.getParameter<double>("MaxStripNoiseSignificanceCut")),
0035       adProbabCut_(pset.getParameter<double>("AdProbabCut")),
0036       ksProbabCut_(pset.getParameter<double>("KsProbabCut")),
0037       generateRandomHisto_(pset.getParameter<bool>("GenerateRandomHisto")),
0038       jbProbabCut_(pset.getParameter<double>("JbProbabCut")),
0039       chi2ProbabCut_(pset.getParameter<double>("Chi2ProbabCut")),
0040       kurtosisCut_(pset.getParameter<double>("KurtosisCut")),
0041       integralTailCut_(pset.getParameter<double>("IntegralTailCut")),
0042       integralNsigma_(pset.getParameter<int>("IntegralNsigma")),
0043       ashmanDistance_(pset.getParameter<double>("AshmanDistance")),
0044       amplitudeRatio_(pset.getParameter<double>("AmplitudeRatio")) {
0045   LogDebug(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
0046                              << " Set maximum drift of the mean value to: " << maxDriftResidualCut_
0047                              << " Set minimum noise value to: " << minStripNoiseCut_
0048                              << " Set maximum noise value to: " << maxStripNoiseCut_
0049                              << " Set maximum noise significance value to: " << maxStripNoiseSignificanceCut_
0050                              << " Set minimum Anderson-Darling p-value to: " << adProbabCut_
0051                              << " Set minimum Kolmogorov-Smirnov p-value to: " << ksProbabCut_
0052                              << " Set minimum Jacque-Bera p-value to: " << jbProbabCut_
0053                              << " Set minimum Chi2 p-value to: " << chi2ProbabCut_
0054                              << " Set N-sigma for the integral to : " << integralNsigma_
0055                              << " Set maximum integral tail at N-sigma to : " << integralTailCut_
0056                              << " Set maximum Kurtosis to : " << kurtosisCut_;
0057 }
0058 
0059 // ----------------------------------------------------------------------------
0060 //
0061 
0062 void PedsFullNoiseAlgorithm::extract(const std::vector<TH1*>& histos) {
0063   if (!anal()) {
0064     edm::LogWarning(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
0065                                       << " NULL pointer to Analysis object!";
0066     return;
0067   }
0068 
0069   // Check number of histograms --> Pedestal, noise and noise2D
0070   if (histos.size() != 3) {
0071     anal()->addErrorCode(sistrip::numberOfHistos_);
0072   }
0073 
0074   // Extract FED key from histo title --> i.e. APV pairs or LLD channel
0075   if (!histos.empty()) {
0076     anal()->fedKey(extractFedKey(histos.front()));
0077   }
0078 
0079   // Extract 1D histograms
0080   std::vector<TH1*>::const_iterator ihis = histos.begin();
0081   for (; ihis != histos.end(); ihis++) {
0082     // Check for NULL pointer
0083     if (!(*ihis)) {
0084       continue;
0085     }
0086 
0087     SiStripHistoTitle title((*ihis)->GetName());
0088     if (title.runType() != sistrip::PEDS_FULL_NOISE) {
0089       anal()->addErrorCode(sistrip::unexpectedTask_);
0090       continue;
0091     }
0092 
0093     // Extract peds histos
0094     if (title.extraInfo().find(sistrip::extrainfo::roughPedestals_) != std::string::npos) {
0095       //@@ something here for rough peds?
0096     } else if (title.extraInfo().find(sistrip::extrainfo::pedestals_) != std::string::npos) {
0097       hPeds_.first = *ihis;
0098       hPeds_.second = (*ihis)->GetName();
0099     } else if (title.extraInfo().find(sistrip::extrainfo::commonMode_) != std::string::npos) {
0100       //@@ something here for CM plots?
0101     } else if (title.extraInfo().find(sistrip::extrainfo::noiseProfile_) != std::string::npos) {
0102       //@@ something here for noise profile plot?
0103       hNoise_.first = *ihis;
0104       hNoise_.second = (*ihis)->GetName();
0105     } else if (title.extraInfo().find(sistrip::extrainfo::noise2D_) != std::string::npos) {
0106       hNoise2D_.first = *ihis;
0107       hNoise2D_.second = (*ihis)->GetName();
0108     } else {
0109       anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
0110     }
0111   }
0112 }
0113 
0114 // resetting vectors
0115 void PedsFullNoiseAlgorithm::reset(PedsFullNoiseAnalysis* ana) {
0116   for (size_t iapv = 0; iapv < ana->peds_.size(); iapv++) {
0117     ana->pedsMean_[iapv] = 0.;
0118     ana->rawMean_[iapv] = 0.;
0119     ana->noiseMean_[iapv] = 0.;
0120     ana->pedsSpread_[iapv] = 0.;
0121     ana->noiseSpread_[iapv] = 0.;
0122     ana->rawSpread_[iapv] = 0.;
0123     ana->pedsMax_[iapv] = 0.;
0124     ana->pedsMin_[iapv] = 0.;
0125     ana->rawMax_[iapv] = 0.;
0126     ana->rawMin_[iapv] = 0.;
0127     ana->noiseMax_[iapv] = 0.;
0128     ana->noiseMin_[iapv] = 0.;
0129 
0130     for (size_t istrip = 0; istrip < ana->peds_[iapv].size(); istrip++) {
0131       ana->peds_[iapv][istrip] = 0.;
0132       ana->noise_[iapv][istrip] = 0.;
0133       ana->raw_[iapv][istrip] = 0.;
0134       ana->adProbab_[iapv][istrip] = 0.;
0135       ana->ksProbab_[iapv][istrip] = 0.;
0136       ana->jbProbab_[iapv][istrip] = 0.;
0137       ana->chi2Probab_[iapv][istrip] = 0.;
0138       ana->residualRMS_[iapv][istrip] = 0.;
0139       ana->residualSigmaGaus_[iapv][istrip] = 0.;
0140       ana->noiseSignificance_[iapv][istrip] = 0.;
0141       ana->residualMean_[iapv][istrip] = 0.;
0142       ana->residualSkewness_[iapv][istrip] = 0.;
0143       ana->residualKurtosis_[iapv][istrip] = 0.;
0144       ana->residualIntegralNsigma_[iapv][istrip] = 0.;
0145       ana->residualIntegral_[iapv][istrip] = 0.;
0146       ana->deadStripBit_[iapv][istrip] = 0;
0147       ana->badStripBit_[iapv][istrip] = 0;
0148     }
0149   }
0150 }
0151 
0152 // -----------------------------------------------------------------------------
0153 //
0154 void PedsFullNoiseAlgorithm::analyse() {
0155   // check base analysis object
0156   if (!anal()) {
0157     edm::LogWarning(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
0158                                       << " NULL pointer to base Analysis object!";
0159     return;
0160   }
0161 
0162   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0163   PedsFullNoiseAnalysis* ana = dynamic_cast<PedsFullNoiseAnalysis*>(tmp);
0164 
0165   // check PedsFullNoiseAnalysis object
0166   if (!ana) {
0167     edm::LogWarning(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
0168                                       << " NULL pointer to derived Analysis object!";
0169     return;
0170   }
0171 
0172   // check if the histograms exists
0173   if (!hPeds_.first) {
0174     ana->addErrorCode(sistrip::nullPtr_);
0175     return;
0176   }
0177 
0178   if (!hNoise_.first) {
0179     ana->addErrorCode(sistrip::nullPtr_);
0180     return;
0181   }
0182 
0183   if (!hNoise2D_.first) {
0184     ana->addErrorCode(sistrip::nullPtr_);
0185     return;
0186   }
0187 
0188   // take the histograms
0189   TProfile* histoPeds = dynamic_cast<TProfile*>(hPeds_.first);
0190   TProfile* histoNoiseMean = dynamic_cast<TProfile*>(hNoise_.first);
0191   TH2S* histoNoise = dynamic_cast<TH2S*>(hNoise2D_.first);
0192 
0193   // Make sanity checks about pointers
0194   if (not histoPeds) {
0195     ana->addErrorCode(sistrip::nullPtr_);
0196     return;
0197   }
0198 
0199   if (not histoNoiseMean) {
0200     ana->addErrorCode(sistrip::nullPtr_);
0201     return;
0202   }
0203 
0204   if (not histoNoise) {
0205     ana->addErrorCode(sistrip::nullPtr_);
0206     return;
0207   }
0208 
0209   // check the binning  --> each x-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
0210   if (histoPeds->GetNbinsX() != 256) {
0211     ana->addErrorCode(sistrip::numberOfBins_);
0212     return;
0213   }
0214 
0215   //check the binning  --> each x-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
0216   if (histoNoiseMean->GetNbinsX() != 256) {
0217     ana->addErrorCode(sistrip::numberOfBins_);
0218     return;
0219   }
0220 
0221   //check the binning  --> each y-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
0222   if (histoNoise->GetNbinsY() != 256) {
0223     ana->addErrorCode(sistrip::numberOfBins_);
0224     return;
0225   }
0226 
0227   //Reset values
0228   reset(ana);
0229 
0230   // loop on each strip
0231   uint32_t apvID = -1;
0232 
0233   // Save basic information at strip / APV level
0234   vector<float> ped_max;
0235   vector<float> ped_min;
0236   vector<float> raw_max;
0237   vector<float> raw_min;
0238   vector<float> noise_max;
0239   vector<float> noise_min;
0240 
0241   // loop on each strip in the lldChannel
0242   for (int iStrip = 0; iStrip < histoPeds->GetNbinsX(); iStrip++) {
0243     if (iStrip < histoPeds->GetNbinsX() / 2)
0244       apvID = 0;
0245     else
0246       apvID = 1;
0247 
0248     int stripBin = 0;
0249     if (iStrip >= 128)
0250       stripBin = iStrip - 128;
0251     else
0252       stripBin = iStrip;
0253 
0254     ana->peds_[apvID][stripBin] = histoPeds->GetBinContent(iStrip + 1);        // pedestal value
0255     ana->noise_[apvID][stripBin] = histoNoiseMean->GetBinContent(iStrip + 1);  // noise value
0256     ana->raw_[apvID][stripBin] = histoPeds->GetBinError(iStrip + 1);           // raw noise value
0257 
0258     ana->pedsMean_[apvID] += ana->peds_[apvID][stripBin];    // mean pedestal
0259     ana->rawMean_[apvID] += ana->raw_[apvID][stripBin];      // mean raw noise
0260     ana->noiseMean_[apvID] += ana->noise_[apvID][stripBin];  // mean noise
0261 
0262     // max pedestal
0263     if (ped_max.size() < apvID + 1)
0264       ped_max.push_back(ana->peds_[apvID][stripBin]);
0265     else {
0266       if (ana->peds_[apvID][stripBin] > ped_max.at(apvID))
0267         ped_max.at(apvID) = ana->peds_[apvID][stripBin];
0268     }
0269 
0270     // min pedestal
0271     if (ped_min.size() < apvID + 1)
0272       ped_min.push_back(ana->peds_[apvID][stripBin]);
0273     else {
0274       if (ana->peds_[apvID][stripBin] < ped_min.at(apvID))
0275         ped_min.at(apvID) = ana->peds_[apvID][stripBin];  // min pedestal
0276     }
0277 
0278     // max noise
0279     if (noise_max.size() < apvID + 1)
0280       noise_max.push_back(ana->noise_[apvID][stripBin]);
0281     else {
0282       if (ana->noise_[apvID][stripBin] > noise_max.at(apvID))
0283         noise_max.at(apvID) = ana->noise_[apvID][stripBin];
0284     }
0285 
0286     // min noise
0287     if (noise_min.size() < apvID + 1)
0288       noise_min.push_back(ana->noise_[apvID][stripBin]);
0289     else {
0290       if (ana->noise_[apvID][stripBin] < noise_min.at(apvID))
0291         noise_min.at(apvID) = ana->noise_[apvID][stripBin];
0292     }
0293 
0294     // max raw
0295     if (raw_max.size() < apvID + 1)
0296       raw_max.push_back(ana->raw_[apvID][stripBin]);
0297     else {
0298       if (ana->raw_[apvID][stripBin] > raw_max.at(apvID))
0299         raw_max.at(apvID) = ana->raw_[apvID][stripBin];
0300     }
0301 
0302     // min raw
0303     if (raw_min.size() < apvID + 1)
0304       raw_min.push_back(ana->raw_[apvID][stripBin]);
0305     else {
0306       if (ana->raw_[apvID][stripBin] < raw_min.at(apvID))
0307         raw_min.at(apvID) = ana->raw_[apvID][stripBin];
0308     }
0309   }
0310 
0311   // Mean values
0312   for (unsigned int iApv = 0; iApv < ana->pedsMean_.size(); iApv++) {
0313     ana->pedsMean_.at(iApv) /= (ana->peds_[iApv].size());    // calculate mean pedestal per APV
0314     ana->rawMean_.at(iApv) /= (ana->raw_[iApv].size());      // calculate mean raw noise per APV
0315     ana->noiseMean_.at(iApv) /= (ana->noise_[iApv].size());  // calculate mean noise per APV
0316   }
0317 
0318   // Min and Max
0319   for (unsigned int iApv = 0; iApv < ped_max.size(); iApv++) {
0320     if (ped_max.at(iApv) > sistrip::maximum_)
0321       ana->pedsMax_.at(iApv) = sistrip::maximum_;
0322     else if (ped_max.at(iApv) < -1. * sistrip::maximum_)
0323       ana->pedsMax_.at(iApv) = -1. * sistrip::maximum_;
0324     else
0325       ana->pedsMax_.at(iApv) = ped_max.at(iApv);
0326 
0327     if (ped_min.at(iApv) > sistrip::maximum_)
0328       ana->pedsMin_.at(iApv) = sistrip::maximum_;
0329     else if (ped_min.at(iApv) < -1. * sistrip::maximum_)
0330       ana->pedsMin_.at(iApv) = -1. * sistrip::maximum_;
0331     else
0332       ana->pedsMin_.at(iApv) = ped_min.at(iApv);
0333 
0334     if (noise_max.at(iApv) > sistrip::maximum_)
0335       ana->noiseMax_.at(iApv) = sistrip::maximum_;
0336     else if (noise_max.at(iApv) < -1. * sistrip::maximum_)
0337       ana->noiseMax_.at(iApv) = -1. * sistrip::maximum_;
0338     else
0339       ana->noiseMax_.at(iApv) = noise_max.at(iApv);
0340 
0341     if (noise_min.at(iApv) > sistrip::maximum_)
0342       ana->noiseMin_.at(iApv) = sistrip::maximum_;
0343     else if (noise_min.at(iApv) < -1. * sistrip::maximum_)
0344       ana->noiseMin_.at(iApv) = -1. * sistrip::maximum_;
0345     else
0346       ana->noiseMin_.at(iApv) = noise_min.at(iApv);
0347 
0348     if (raw_max.at(iApv) > sistrip::maximum_)
0349       ana->rawMax_.at(iApv) = sistrip::maximum_;
0350     else if (raw_max.at(iApv) < -1. * sistrip::maximum_)
0351       ana->rawMax_.at(iApv) = -1. * sistrip::maximum_;
0352     else
0353       ana->rawMax_.at(iApv) = raw_max.at(iApv);
0354 
0355     if (raw_min.at(iApv) > sistrip::maximum_)
0356       ana->rawMin_.at(iApv) = sistrip::maximum_;
0357     else if (raw_min.at(iApv) < -1. * sistrip::maximum_)
0358       ana->rawMin_.at(iApv) = -1. * sistrip::maximum_;
0359     else
0360       ana->rawMin_.at(iApv) = raw_min.at(iApv);
0361   }
0362 
0363   // Calculate the spread for noise and pedestal
0364   for (int iStrip = 0; iStrip < histoNoiseMean->GetNbinsX(); iStrip++) {
0365     if (iStrip < histoNoiseMean->GetNbinsX() / 2)
0366       apvID = 0;
0367     else
0368       apvID = 1;
0369     ana->pedsSpread_[apvID] += pow(histoPeds->GetBinContent(iStrip + 1) - ana->pedsMean_.at(apvID), 2);
0370     ana->noiseSpread_[apvID] += pow(histoNoiseMean->GetBinContent(iStrip + 1) - ana->noiseMean_.at(apvID), 2);
0371     ana->rawSpread_[apvID] += pow(histoPeds->GetBinError(iStrip + 1) - ana->rawMean_.at(apvID), 2);
0372   }
0373 
0374   for (unsigned int iApv = 0; iApv < ana->pedsSpread_.size(); iApv++) {
0375     ana->pedsSpread_[iApv] = sqrt(ana->pedsSpread_[iApv]) / sqrt(ana->peds_[iApv].size() - 1);
0376     ana->noiseSpread_[iApv] = sqrt(ana->noiseSpread_[iApv]) / sqrt(ana->noise_[iApv].size() - 1);
0377     ana->rawSpread_[iApv] = sqrt(ana->rawSpread_[iApv]) / sqrt(ana->raw_[iApv].size() - 1);
0378   }
0379 
0380   // loop on each strip in the lldChannel
0381   TH1S* histoResidualStrip = new TH1S("histoResidualStrip",
0382                                       "",
0383                                       histoNoise->GetNbinsX(),
0384                                       histoNoise->GetXaxis()->GetXmin(),
0385                                       histoNoise->GetXaxis()->GetXmax());
0386   histoResidualStrip->Sumw2();
0387   histoResidualStrip->SetDirectory(nullptr);
0388   TF1* fitFunc = new TF1("fitFunc", "gaus(0)", histoNoise->GetXaxis()->GetXmin(), histoNoise->GetXaxis()->GetXmax());
0389   TF1* fit2Gaus = nullptr;
0390   TH1F* randomHisto = nullptr;
0391   TFitResultPtr result;
0392 
0393   for (int iStrip = 0; iStrip < histoNoise->GetNbinsY(); iStrip++) {
0394     // tell which APV
0395     if (iStrip < histoNoise->GetNbinsY() / 2)
0396       apvID = 0;
0397     else
0398       apvID = 1;
0399     histoResidualStrip->Reset();
0400 
0401     int stripBin = 0;
0402     if (iStrip >= 128)
0403       stripBin = iStrip - 128;
0404     else
0405       stripBin = iStrip;
0406 
0407     for (int iBinX = 0; iBinX < histoNoise->GetNbinsX(); iBinX++) {
0408       histoResidualStrip->SetBinContent(iBinX + 1, histoNoise->GetBinContent(iBinX + 1, iStrip + 1));
0409       histoResidualStrip->SetBinError(iBinX + 1, histoNoise->GetBinError(iBinX + 1, iStrip + 1));
0410     }
0411 
0412     if (histoResidualStrip->Integral() == 0) {  // dead strip --> i.e. no data
0413 
0414       // set default values
0415       ana->adProbab_[apvID][stripBin] = 0;
0416       ana->ksProbab_[apvID][stripBin] = 0;
0417       ana->jbProbab_[apvID][stripBin] = 0;
0418       ana->chi2Probab_[apvID][stripBin] = 0;
0419       ana->noiseSignificance_[apvID][stripBin] = 0;
0420       ana->residualMean_[apvID][stripBin] = 0;
0421       ana->residualRMS_[apvID][stripBin] = 0;
0422       ana->residualSigmaGaus_[apvID][stripBin] = 0;
0423       ana->residualSkewness_[apvID][stripBin] = 0;
0424       ana->residualKurtosis_[apvID][stripBin] = 0;
0425       ana->residualIntegralNsigma_[apvID][stripBin] = 0;
0426       ana->residualIntegral_[apvID][stripBin] = 0;
0427       ana->deadStrip_[apvID].push_back(stripBin);
0428       ana->deadStripBit_[apvID][stripBin] = 1;
0429       ana->badStripBit_[apvID][stripBin] = 0;
0430 
0431       SiStripFecKey fec_key(ana->fecKey());
0432       LogTrace(mlDqmClient_) << "DeadStrip: fecCrate "
0433                              << " " << fec_key.fecCrate() << " fecSlot " << fec_key.fecSlot() << " fecRing "
0434                              << fec_key.fecRing() << " ccuAddr " << fec_key.ccuAddr() << " ccChan  "
0435                              << fec_key.ccuChan() << " lldChan " << fec_key.lldChan() << " apvID " << apvID
0436                              << " stripID " << iStrip;
0437 
0438       continue;
0439     }
0440 
0441     // set / calculated basic quantities
0442     ana->residualMean_[apvID][stripBin] = histoResidualStrip->GetMean();
0443     ana->residualRMS_[apvID][stripBin] = histoResidualStrip->GetRMS();
0444     ana->residualSkewness_[apvID][stripBin] = histoResidualStrip->GetSkewness();
0445     ana->residualKurtosis_[apvID][stripBin] = histoResidualStrip->GetKurtosis();
0446     ana->noiseSignificance_[apvID][stripBin] =
0447         (ana->noise_[apvID][stripBin] - ana->noiseMean_[apvID]) / ana->noiseSpread_[apvID];
0448     ana->residualIntegral_[apvID][stripBin] = histoResidualStrip->Integral();
0449     ana->residualIntegralNsigma_[apvID][stripBin] =
0450         (histoResidualStrip->Integral(histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin] +
0451                                                                   ana->residualRMS_[apvID][stripBin] * integralNsigma_),
0452                                       histoResidualStrip->GetNbinsX() + 1) +
0453          histoResidualStrip->Integral(
0454              0,
0455              histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin] -
0456                                          ana->residualRMS_[apvID][stripBin] * integralNsigma_))) /
0457         ana->residualIntegral_[apvID][stripBin];
0458 
0459     // performing a Gaussian fit of the residual distribution
0460     fitFunc->SetRange(histoNoise->GetXaxis()->GetXmin(), histoNoise->GetXaxis()->GetXmax());
0461     fitFunc->SetParameters(ana->residualIntegral_[apvID][stripBin],
0462                            ana->residualMean_[apvID][stripBin],
0463                            ana->residualRMS_[apvID][stripBin]);
0464     result = histoResidualStrip->Fit(fitFunc, "QSRN");
0465 
0466     // Good gaussian fit
0467     if (result.Get()) {
0468       ana->residualSigmaGaus_[apvID][stripBin] = fitFunc->GetParameter(2);
0469       ana->chi2Probab_[apvID][stripBin] = result->Prob();
0470 
0471       // jacque bera probability
0472       float jbVal =
0473           (ana->residualIntegral_[apvID][stripBin] / 6) *
0474           (pow(ana->residualSkewness_[apvID][stripBin], 2) + pow(ana->residualKurtosis_[apvID][stripBin], 2) / 4);
0475       ana->jbProbab_[apvID][stripBin] = ROOT::Math::chisquared_cdf_c(jbVal, 2);
0476 
0477       //Kolmogorov Smirnov and Anderson Darlong
0478       if (randomHisto == nullptr)
0479         randomHisto = (TH1F*)histoResidualStrip->Clone("randomHisto");
0480       randomHisto->Reset();
0481       randomHisto->SetDirectory(nullptr);
0482 
0483       if (generateRandomHisto_) {  ///
0484         randomHisto->FillRandom("fitFunc", histoResidualStrip->Integral());
0485         if (randomHisto->Integral() != 0) {
0486           ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto, "N");
0487           ana->adProbab_[apvID][stripBin] = histoResidualStrip->AndersonDarlingTest(randomHisto);
0488         } else {
0489           ana->ksProbab_[apvID][stripBin] = 0;
0490           ana->adProbab_[apvID][stripBin] = 0;
0491         }
0492 
0493       } else {
0494         randomHisto->Add(fitFunc);
0495         if (randomHisto->Integral() != 0) {
0496           ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto, "N");
0497           ROOT::Fit::BinData data1;
0498           ROOT::Fit::BinData data2;
0499           ROOT::Fit::FillData(data1, histoResidualStrip, nullptr);
0500           data2.Initialize(randomHisto->GetNbinsX() + 1, 1);
0501           for (int ibin = 0; ibin < randomHisto->GetNbinsX(); ibin++) {
0502             if (histoResidualStrip->GetBinContent(ibin + 1) != 0 or randomHisto->GetBinContent(ibin + 1) >= 1)
0503               data2.Add(randomHisto->GetBinCenter(ibin + 1),
0504                         randomHisto->GetBinContent(ibin + 1),
0505                         randomHisto->GetBinError(ibin + 1));
0506           }
0507 
0508           double probab, value;
0509           ROOT::Math::GoFTest::AndersonDarling2SamplesTest(data1, data2, probab, value);
0510           ana->adProbab_[apvID][stripBin] = probab;
0511         } else {
0512           ana->ksProbab_[apvID][stripBin] = 0;
0513           ana->adProbab_[apvID][stripBin] = 0;
0514         }
0515       }
0516     }
0517 
0518     // start applying selections storing output
0519     bool badStripFlag = false;
0520     ana->deadStripBit_[apvID][stripBin] = 0;  // is not dead if the strip has data
0521 
0522     if (fabs(ana->residualMean_[apvID][stripBin]) > maxDriftResidualCut_ and not badStripFlag) {  //mean value
0523       ana->shiftedStrip_[apvID].push_back(stripBin);
0524       badStripFlag = true;
0525     }
0526 
0527     if (ana->residualRMS_[apvID][stripBin] < minStripNoiseCut_ and not badStripFlag) {  // low noise
0528       ana->lowNoiseStrip_[apvID].push_back(stripBin);
0529       badStripFlag = true;
0530     }
0531 
0532     if (ana->residualRMS_[apvID][stripBin] > maxStripNoiseCut_ and not badStripFlag) {  // large noise
0533       ana->largeNoiseStrip_[apvID].push_back(stripBin);
0534       badStripFlag = true;
0535     }
0536 
0537     if (fabs(ana->noiseSignificance_[apvID][stripBin]) > maxStripNoiseSignificanceCut_ and
0538         not badStripFlag) {  // large noise significance
0539       ana->largeNoiseSignificance_[apvID].push_back(stripBin);
0540       badStripFlag = true;
0541     }
0542 
0543     if (result.Get() and result->Status() != 0)  // bad fit status
0544       ana->badFitStatus_[apvID].push_back(stripBin);
0545 
0546     if (ana->adProbab_[apvID][stripBin] < adProbabCut_ and not badStripFlag)  // bad AD p-value --> store the strip-id
0547       ana->badADProbab_[apvID].push_back(stripBin);
0548 
0549     if (ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and not badStripFlag)  // bad KS p-value --> store the strip-id
0550       ana->badKSProbab_[apvID].push_back(stripBin);
0551 
0552     if (ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and not badStripFlag)  // bad JB p-value  --> store the strip-id
0553       ana->badJBProbab_[apvID].push_back(stripBin);
0554 
0555     if (ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_ and
0556         not badStripFlag)  // bad CHI2 p-value  --> store the strip-id
0557       ana->badChi2Probab_[apvID].push_back(stripBin);
0558 
0559     if (ana->adProbab_[apvID][stripBin] < adProbabCut_ and ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and
0560         ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_)
0561       badStripFlag = true;  // bad strip is flagged as bad by all the methods
0562 
0563     if (ana->residualKurtosis_[apvID][stripBin] > kurtosisCut_ and
0564         ana->residualIntegralNsigma_[apvID][stripBin] > integralTailCut_ and not badStripFlag) {  // bad tails
0565       ana->badTailStrip_[apvID].push_back(stripBin);
0566       badStripFlag = true;
0567     }
0568 
0569     if (badStripFlag) {  // loop for double peaked
0570 
0571       fit2Gaus = new TF1("dgaus",
0572                          "[0]*exp(-((x-[1])*(x-[1]))/(2*[2]*[2]))+[3]*exp(-((x-[4])*(x-[4]))/(2*[5]*[5]))",
0573                          histoNoise->GetXaxis()->GetXmin(),
0574                          histoNoise->GetXaxis()->GetXmax());
0575       fit2Gaus->SetParameter(0, fitFunc->GetParameter(0) / 2);
0576       fit2Gaus->SetParameter(3, fitFunc->GetParameter(0) / 2);
0577       fit2Gaus->SetParameter(1, 1.);
0578       fit2Gaus->SetParameter(4, -1.);
0579       fit2Gaus->SetParameter(2, fitFunc->GetParameter(2));
0580       fit2Gaus->SetParameter(5, fitFunc->GetParameter(2));
0581       fit2Gaus->SetParLimits(1, 0., histoNoise->GetXaxis()->GetXmax());
0582       fit2Gaus->SetParLimits(4, histoNoise->GetXaxis()->GetXmin(), 0);
0583       result = histoResidualStrip->Fit(fit2Gaus, "QSR");
0584 
0585       // ashman distance
0586       float ashman = TMath::Power(2, 0.5) * abs(fit2Gaus->GetParameter(1) - fit2Gaus->GetParameter(4)) /
0587                      (sqrt(pow(fit2Gaus->GetParameter(2), 2) + pow(fit2Gaus->GetParameter(5), 2)));
0588       // amplitude
0589       float amplitudeRatio = std::min(fit2Gaus->GetParameter(0), fit2Gaus->GetParameter(3)) /
0590                              std::max(fit2Gaus->GetParameter(0), fit2Gaus->GetParameter(3));
0591 
0592       if (ashman > ashmanDistance_ and amplitudeRatio > amplitudeRatio_)
0593         ana->badDoublePeakStrip_[apvID].push_back(stripBin);
0594     }
0595 
0596     if (badStripFlag) {  // save the final bit
0597 
0598       ana->badStrip_[apvID].push_back(stripBin);
0599       ana->badStripBit_[apvID][stripBin] = 1;
0600 
0601       SiStripFecKey fec_key(ana->fecKey());
0602       LogTrace(mlDqmClient_) << "BadStrip: fecCrate "
0603                              << " " << fec_key.fecCrate() << " fecSlot " << fec_key.fecSlot() << " fecRing "
0604                              << fec_key.fecRing() << " ccuAddr " << fec_key.ccuAddr() << " ccChan  "
0605                              << fec_key.ccuChan() << " lldChan " << fec_key.lldChan() << " apvID " << apvID
0606                              << " stripID " << stripBin;
0607 
0608     } else
0609       ana->badStripBit_[apvID][stripBin] = 0;
0610   }
0611 
0612   ped_max.clear();
0613   ped_min.clear();
0614   raw_max.clear();
0615   raw_min.clear();
0616   noise_max.clear();
0617   noise_min.clear();
0618   if (histoResidualStrip)
0619     delete histoResidualStrip;
0620   if (fitFunc)
0621     delete fitFunc;
0622   if (randomHisto)
0623     delete randomHisto;
0624   if (fit2Gaus)
0625     delete fit2Gaus;
0626 }
[ 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!