Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:56:22 

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   apvID = -1;
0365 
0366   for (int iStrip = 0; iStrip < histoNoiseMean->GetNbinsX(); iStrip++) {
0367     if (iStrip < histoNoiseMean->GetNbinsX() / 2)
0368       apvID = 0;
0369     else
0370       apvID = 1;
0371     ana->pedsSpread_[apvID] += pow(histoPeds->GetBinContent(iStrip + 1) - ana->pedsMean_.at(apvID), 2);
0372     ana->noiseSpread_[apvID] += pow(histoNoiseMean->GetBinContent(iStrip + 1) - ana->noiseMean_.at(apvID), 2);
0373     ana->rawSpread_[apvID] += pow(histoPeds->GetBinError(iStrip + 1) - ana->rawMean_.at(apvID), 2);
0374   }
0375 
0376   for (unsigned int iApv = 0; iApv < ana->pedsSpread_.size(); iApv++) {
0377     ana->pedsSpread_[iApv] = sqrt(ana->pedsSpread_[iApv]) / sqrt(ana->peds_[iApv].size() - 1);
0378     ana->noiseSpread_[iApv] = sqrt(ana->noiseSpread_[iApv]) / sqrt(ana->noise_[iApv].size() - 1);
0379     ana->rawSpread_[iApv] = sqrt(ana->rawSpread_[iApv]) / sqrt(ana->raw_[iApv].size() - 1);
0380   }
0381 
0382   // loop on each strip in the lldChannel
0383   apvID = 0;
0384   TH1S* histoResidualStrip = new TH1S("histoResidualStrip",
0385                                       "",
0386                                       histoNoise->GetNbinsX(),
0387                                       histoNoise->GetXaxis()->GetXmin(),
0388                                       histoNoise->GetXaxis()->GetXmax());
0389   histoResidualStrip->Sumw2();
0390   histoResidualStrip->SetDirectory(nullptr);
0391   TF1* fitFunc = new TF1("fitFunc", "gaus(0)", histoNoise->GetXaxis()->GetXmin(), histoNoise->GetXaxis()->GetXmax());
0392   TF1* fit2Gaus = nullptr;
0393   TH1F* randomHisto = nullptr;
0394   TFitResultPtr result;
0395 
0396   for (int iStrip = 0; iStrip < histoNoise->GetNbinsY(); iStrip++) {
0397     // tell which APV
0398     if (iStrip < histoNoise->GetNbinsY() / 2)
0399       apvID = 0;
0400     else
0401       apvID = 1;
0402     histoResidualStrip->Reset();
0403 
0404     int stripBin = 0;
0405     if (iStrip >= 128)
0406       stripBin = iStrip - 128;
0407     else
0408       stripBin = iStrip;
0409 
0410     for (int iBinX = 0; iBinX < histoNoise->GetNbinsX(); iBinX++) {
0411       histoResidualStrip->SetBinContent(iBinX + 1, histoNoise->GetBinContent(iBinX + 1, iStrip + 1));
0412       histoResidualStrip->SetBinError(iBinX + 1, histoNoise->GetBinError(iBinX + 1, iStrip + 1));
0413     }
0414 
0415     if (histoResidualStrip->Integral() == 0) {  // dead strip --> i.e. no data
0416 
0417       // set default values
0418       ana->adProbab_[apvID][stripBin] = 0;
0419       ana->ksProbab_[apvID][stripBin] = 0;
0420       ana->jbProbab_[apvID][stripBin] = 0;
0421       ana->chi2Probab_[apvID][stripBin] = 0;
0422       ana->noiseSignificance_[apvID][stripBin] = 0;
0423       ana->residualMean_[apvID][stripBin] = 0;
0424       ana->residualRMS_[apvID][stripBin] = 0;
0425       ana->residualSigmaGaus_[apvID][stripBin] = 0;
0426       ana->residualSkewness_[apvID][stripBin] = 0;
0427       ana->residualKurtosis_[apvID][stripBin] = 0;
0428       ana->residualIntegralNsigma_[apvID][stripBin] = 0;
0429       ana->residualIntegral_[apvID][stripBin] = 0;
0430       ana->deadStrip_[apvID].push_back(stripBin);
0431       ana->deadStripBit_[apvID][stripBin] = 1;
0432       ana->badStripBit_[apvID][stripBin] = 0;
0433 
0434       SiStripFecKey fec_key(ana->fecKey());
0435       LogTrace(mlDqmClient_) << "DeadStrip: fecCrate "
0436                              << " " << fec_key.fecCrate() << " fecSlot " << fec_key.fecSlot() << " fecRing "
0437                              << fec_key.fecRing() << " ccuAddr " << fec_key.ccuAddr() << " ccChan  "
0438                              << fec_key.ccuChan() << " lldChan " << fec_key.lldChan() << " apvID " << apvID
0439                              << " stripID " << iStrip;
0440 
0441       continue;
0442     }
0443 
0444     // set / calculated basic quantities
0445     ana->residualMean_[apvID][stripBin] = histoResidualStrip->GetMean();
0446     ana->residualRMS_[apvID][stripBin] = histoResidualStrip->GetRMS();
0447     ana->residualSkewness_[apvID][stripBin] = histoResidualStrip->GetSkewness();
0448     ana->residualKurtosis_[apvID][stripBin] = histoResidualStrip->GetKurtosis();
0449     ana->noiseSignificance_[apvID][stripBin] =
0450         (ana->noise_[apvID][stripBin] - ana->noiseMean_[apvID]) / ana->noiseSpread_[apvID];
0451     ana->residualIntegral_[apvID][stripBin] = histoResidualStrip->Integral();
0452     ana->residualIntegralNsigma_[apvID][stripBin] =
0453         (histoResidualStrip->Integral(histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin] +
0454                                                                   ana->residualRMS_[apvID][stripBin] * integralNsigma_),
0455                                       histoResidualStrip->GetNbinsX() + 1) +
0456          histoResidualStrip->Integral(
0457              0,
0458              histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin] -
0459                                          ana->residualRMS_[apvID][stripBin] * integralNsigma_))) /
0460         ana->residualIntegral_[apvID][stripBin];
0461 
0462     // performing a Gaussian fit of the residual distribution
0463     fitFunc->SetRange(histoNoise->GetXaxis()->GetXmin(), histoNoise->GetXaxis()->GetXmax());
0464     fitFunc->SetParameters(ana->residualIntegral_[apvID][stripBin],
0465                            ana->residualMean_[apvID][stripBin],
0466                            ana->residualRMS_[apvID][stripBin]);
0467     result = histoResidualStrip->Fit(fitFunc, "QSRN");
0468 
0469     // Good gaussian fit
0470     if (result.Get()) {
0471       ana->residualSigmaGaus_[apvID][stripBin] = fitFunc->GetParameter(2);
0472       ana->chi2Probab_[apvID][stripBin] = result->Prob();
0473 
0474       // jacque bera probability
0475       float jbVal =
0476           (ana->residualIntegral_[apvID][stripBin] / 6) *
0477           (pow(ana->residualSkewness_[apvID][stripBin], 2) + pow(ana->residualKurtosis_[apvID][stripBin], 2) / 4);
0478       ana->jbProbab_[apvID][stripBin] = ROOT::Math::chisquared_cdf_c(jbVal, 2);
0479 
0480       //Kolmogorov Smirnov and Anderson Darlong
0481       if (randomHisto == nullptr)
0482         randomHisto = (TH1F*)histoResidualStrip->Clone("randomHisto");
0483       randomHisto->Reset();
0484       randomHisto->SetDirectory(nullptr);
0485 
0486       if (generateRandomHisto_) {  ///
0487         randomHisto->FillRandom("fitFunc", histoResidualStrip->Integral());
0488         if (randomHisto->Integral() != 0) {
0489           ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto, "N");
0490           ana->adProbab_[apvID][stripBin] = histoResidualStrip->AndersonDarlingTest(randomHisto);
0491         } else {
0492           ana->ksProbab_[apvID][stripBin] = 0;
0493           ana->adProbab_[apvID][stripBin] = 0;
0494         }
0495 
0496       } else {
0497         randomHisto->Add(fitFunc);
0498         if (randomHisto->Integral() != 0) {
0499           ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto, "N");
0500           ROOT::Fit::BinData data1;
0501           ROOT::Fit::BinData data2;
0502           ROOT::Fit::FillData(data1, histoResidualStrip, nullptr);
0503           data2.Initialize(randomHisto->GetNbinsX() + 1, 1);
0504           for (int ibin = 0; ibin < randomHisto->GetNbinsX(); ibin++) {
0505             if (histoResidualStrip->GetBinContent(ibin + 1) != 0 or randomHisto->GetBinContent(ibin + 1) >= 1)
0506               data2.Add(randomHisto->GetBinCenter(ibin + 1),
0507                         randomHisto->GetBinContent(ibin + 1),
0508                         randomHisto->GetBinError(ibin + 1));
0509           }
0510 
0511           double probab, value;
0512           ROOT::Math::GoFTest::AndersonDarling2SamplesTest(data1, data2, probab, value);
0513           ana->adProbab_[apvID][stripBin] = probab;
0514         } else {
0515           ana->ksProbab_[apvID][stripBin] = 0;
0516           ana->adProbab_[apvID][stripBin] = 0;
0517         }
0518       }
0519     }
0520 
0521     // start applying selections storing output
0522     bool badStripFlag = false;
0523     ana->deadStripBit_[apvID][stripBin] = 0;  // is not dead if the strip has data
0524 
0525     if (fabs(ana->residualMean_[apvID][stripBin]) > maxDriftResidualCut_ and not badStripFlag) {  //mean value
0526       ana->shiftedStrip_[apvID].push_back(stripBin);
0527       badStripFlag = true;
0528     }
0529 
0530     if (ana->residualRMS_[apvID][stripBin] < minStripNoiseCut_ and not badStripFlag) {  // low noise
0531       ana->lowNoiseStrip_[apvID].push_back(stripBin);
0532       badStripFlag = true;
0533     }
0534 
0535     if (ana->residualRMS_[apvID][stripBin] > maxStripNoiseCut_ and not badStripFlag) {  // large noise
0536       ana->largeNoiseStrip_[apvID].push_back(stripBin);
0537       badStripFlag = true;
0538     }
0539 
0540     if (fabs(ana->noiseSignificance_[apvID][stripBin]) > maxStripNoiseSignificanceCut_ and
0541         not badStripFlag) {  // large noise significance
0542       ana->largeNoiseSignificance_[apvID].push_back(stripBin);
0543       badStripFlag = true;
0544     }
0545 
0546     if (result.Get() and result->Status() != 0)  // bad fit status
0547       ana->badFitStatus_[apvID].push_back(stripBin);
0548 
0549     if (ana->adProbab_[apvID][stripBin] < adProbabCut_ and not badStripFlag)  // bad AD p-value --> store the strip-id
0550       ana->badADProbab_[apvID].push_back(stripBin);
0551 
0552     if (ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and not badStripFlag)  // bad KS p-value --> store the strip-id
0553       ana->badKSProbab_[apvID].push_back(stripBin);
0554 
0555     if (ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and not badStripFlag)  // bad JB p-value  --> store the strip-id
0556       ana->badJBProbab_[apvID].push_back(stripBin);
0557 
0558     if (ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_ and
0559         not badStripFlag)  // bad CHI2 p-value  --> store the strip-id
0560       ana->badChi2Probab_[apvID].push_back(stripBin);
0561 
0562     if (ana->adProbab_[apvID][stripBin] < adProbabCut_ and ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and
0563         ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_)
0564       badStripFlag = true;  // bad strip is flagged as bad by all the methods
0565 
0566     if (ana->residualKurtosis_[apvID][stripBin] > kurtosisCut_ and
0567         ana->residualIntegralNsigma_[apvID][stripBin] > integralTailCut_ and not badStripFlag) {  // bad tails
0568       ana->badTailStrip_[apvID].push_back(stripBin);
0569       badStripFlag = true;
0570     }
0571 
0572     if (badStripFlag) {  // loop for double peaked
0573 
0574       fit2Gaus = new TF1("dgaus",
0575                          "[0]*exp(-((x-[1])*(x-[1]))/(2*[2]*[2]))+[3]*exp(-((x-[4])*(x-[4]))/(2*[5]*[5]))",
0576                          histoNoise->GetXaxis()->GetXmin(),
0577                          histoNoise->GetXaxis()->GetXmax());
0578       fit2Gaus->SetParameter(0, fitFunc->GetParameter(0) / 2);
0579       fit2Gaus->SetParameter(3, fitFunc->GetParameter(0) / 2);
0580       fit2Gaus->SetParameter(1, 1.);
0581       fit2Gaus->SetParameter(4, -1.);
0582       fit2Gaus->SetParameter(2, fitFunc->GetParameter(2));
0583       fit2Gaus->SetParameter(5, fitFunc->GetParameter(2));
0584       fit2Gaus->SetParLimits(1, 0., histoNoise->GetXaxis()->GetXmax());
0585       fit2Gaus->SetParLimits(4, histoNoise->GetXaxis()->GetXmin(), 0);
0586       result = histoResidualStrip->Fit(fit2Gaus, "QSR");
0587 
0588       // ashman distance
0589       float ashman = TMath::Power(2, 0.5) * abs(fit2Gaus->GetParameter(1) - fit2Gaus->GetParameter(4)) /
0590                      (sqrt(pow(fit2Gaus->GetParameter(2), 2) + pow(fit2Gaus->GetParameter(5), 2)));
0591       // amplitude
0592       float amplitudeRatio = std::min(fit2Gaus->GetParameter(0), fit2Gaus->GetParameter(3)) /
0593                              std::max(fit2Gaus->GetParameter(0), fit2Gaus->GetParameter(3));
0594 
0595       if (ashman > ashmanDistance_ and amplitudeRatio > amplitudeRatio_)
0596         ana->badDoublePeakStrip_[apvID].push_back(stripBin);
0597     }
0598 
0599     if (badStripFlag) {  // save the final bit
0600 
0601       ana->badStrip_[apvID].push_back(stripBin);
0602       ana->badStripBit_[apvID][stripBin] = 1;
0603 
0604       SiStripFecKey fec_key(ana->fecKey());
0605       LogTrace(mlDqmClient_) << "BadStrip: fecCrate "
0606                              << " " << fec_key.fecCrate() << " fecSlot " << fec_key.fecSlot() << " fecRing "
0607                              << fec_key.fecRing() << " ccuAddr " << fec_key.ccuAddr() << " ccChan  "
0608                              << fec_key.ccuChan() << " lldChan " << fec_key.lldChan() << " apvID " << apvID
0609                              << " stripID " << stripBin;
0610 
0611     } else
0612       ana->badStripBit_[apvID][stripBin] = 0;
0613   }
0614 
0615   ped_max.clear();
0616   ped_min.clear();
0617   raw_max.clear();
0618   raw_min.clear();
0619   noise_max.clear();
0620   noise_min.clear();
0621   if (histoResidualStrip)
0622     delete histoResidualStrip;
0623   if (fitFunc)
0624     delete fitFunc;
0625   if (randomHisto)
0626     delete randomHisto;
0627   if (fit2Gaus)
0628     delete fit2Gaus;
0629 }
[ 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!