text: #include "DQM/SiStripCommissioningAnalysis/interface/PedsFullNoiseAlgorithm.h"
  #include "CondFormats/SiStripObjects/interface/PedsFullNoiseAnalysis.h"
  #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
  #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
  #include "FWCore/MessageLogger/interface/MessageLogger.h"
  
  #include "TProfile.h"
  #include "TH1.h"
  #include "TH2.h"
  #include "TF1.h"
  #include "TFitResult.h"
  #include "TMath.h"
  #include "Math/DistFunc.h"
  #include "Math/ProbFuncMathCore.h"
  #include "Fit/BinData.h"
  #include "HFitInterface.h"
  #include "Math/GoFTest.h"
  #include <iostream>
  #include <iomanip>
  #include <cmath>
  
  using namespace sistrip;
  using namespace std;
  // ----------------------------------------------------------------------------
  //
  PedsFullNoiseAlgorithm::PedsFullNoiseAlgorithm(const edm::ParameterSet& pset, PedsFullNoiseAnalysis* const anal)
      : CommissioningAlgorithm(anal),
        hPeds_(nullptr, ""),
        hNoise_(nullptr, ""),
        hNoise2D_(nullptr, ""),
        maxDriftResidualCut_(pset.getParameter<double>("MaxDriftResidualCut")),
        minStripNoiseCut_(pset.getParameter<double>("MinStripNoiseCut")),
        maxStripNoiseCut_(pset.getParameter<double>("MaxStripNoiseCut")),
        maxStripNoiseSignificanceCut_(pset.getParameter<double>("MaxStripNoiseSignificanceCut")),
        adProbabCut_(pset.getParameter<double>("AdProbabCut")),
        ksProbabCut_(pset.getParameter<double>("KsProbabCut")),
        generateRandomHisto_(pset.getParameter<bool>("GenerateRandomHisto")),
        jbProbabCut_(pset.getParameter<double>("JbProbabCut")),
        chi2ProbabCut_(pset.getParameter<double>("Chi2ProbabCut")),
        kurtosisCut_(pset.getParameter<double>("KurtosisCut")),
        integralTailCut_(pset.getParameter<double>("IntegralTailCut")),
        integralNsigma_(pset.getParameter<int>("IntegralNsigma")),
        ashmanDistance_(pset.getParameter<double>("AshmanDistance")),
        amplitudeRatio_(pset.getParameter<double>("AmplitudeRatio")) {
    LogDebug(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
                               << " Set maximum drift of the mean value to: " << maxDriftResidualCut_
                               << " Set minimum noise value to: " << minStripNoiseCut_
                               << " Set maximum noise value to: " << maxStripNoiseCut_
                               << " Set maximum noise significance value to: " << maxStripNoiseSignificanceCut_
                               << " Set minimum Anderson-Darling p-value to: " << adProbabCut_
                               << " Set minimum Kolmogorov-Smirnov p-value to: " << ksProbabCut_
                               << " Set minimum Jacque-Bera p-value to: " << jbProbabCut_
                               << " Set minimum Chi2 p-value to: " << chi2ProbabCut_
                               << " Set N-sigma for the integral to : " << integralNsigma_
                               << " Set maximum integral tail at N-sigma to : " << integralTailCut_
                               << " Set maximum Kurtosis to : " << kurtosisCut_;
  }
  
  // ----------------------------------------------------------------------------
  //
  
  void PedsFullNoiseAlgorithm::extract(const std::vector<TH1*>& histos) {
    if (!anal()) {
      edm::LogWarning(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
                                        << " NULL pointer to Analysis object!";
      return;
    }
  
    // Check number of histograms --> Pedestal, noise and noise2D
    if (histos.size() != 3) {
      anal()->addErrorCode(sistrip::numberOfHistos_);
    }
  
    // Extract FED key from histo title --> i.e. APV pairs or LLD channel
    if (!histos.empty()) {
      anal()->fedKey(extractFedKey(histos.front()));
    }
  
    // Extract 1D histograms
    std::vector<TH1*>::const_iterator ihis = histos.begin();
    for (; ihis != histos.end(); ihis++) {
      // Check for NULL pointer
      if (!(*ihis)) {
        continue;
      }
  
      SiStripHistoTitle title((*ihis)->GetName());
      if (title.runType() != sistrip::PEDS_FULL_NOISE) {
        anal()->addErrorCode(sistrip::unexpectedTask_);
        continue;
      }
  
      // Extract peds histos
      if (title.extraInfo().find(sistrip::extrainfo::roughPedestals_) != std::string::npos) {
        //@@ something here for rough peds?
      } else if (title.extraInfo().find(sistrip::extrainfo::pedestals_) != std::string::npos) {
        hPeds_.first = *ihis;
        hPeds_.second = (*ihis)->GetName();
      } else if (title.extraInfo().find(sistrip::extrainfo::commonMode_) != std::string::npos) {
        //@@ something here for CM plots?
      } else if (title.extraInfo().find(sistrip::extrainfo::noiseProfile_) != std::string::npos) {
        //@@ something here for noise profile plot?
        hNoise_.first = *ihis;
        hNoise_.second = (*ihis)->GetName();
      } else if (title.extraInfo().find(sistrip::extrainfo::noise2D_) != std::string::npos) {
        hNoise2D_.first = *ihis;
        hNoise2D_.second = (*ihis)->GetName();
      } else {
        anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
      }
    }
  }
  
  // resetting vectors
  void PedsFullNoiseAlgorithm::reset(PedsFullNoiseAnalysis* ana) {
    for (size_t iapv = 0; iapv < ana->peds_.size(); iapv++) {
      ana->pedsMean_[iapv] = 0.;
      ana->rawMean_[iapv] = 0.;
      ana->noiseMean_[iapv] = 0.;
      ana->pedsSpread_[iapv] = 0.;
      ana->noiseSpread_[iapv] = 0.;
      ana->rawSpread_[iapv] = 0.;
      ana->pedsMax_[iapv] = 0.;
      ana->pedsMin_[iapv] = 0.;
      ana->rawMax_[iapv] = 0.;
      ana->rawMin_[iapv] = 0.;
      ana->noiseMax_[iapv] = 0.;
      ana->noiseMin_[iapv] = 0.;
  
      for (size_t istrip = 0; istrip < ana->peds_[iapv].size(); istrip++) {
        ana->peds_[iapv][istrip] = 0.;
        ana->noise_[iapv][istrip] = 0.;
        ana->raw_[iapv][istrip] = 0.;
        ana->adProbab_[iapv][istrip] = 0.;
        ana->ksProbab_[iapv][istrip] = 0.;
        ana->jbProbab_[iapv][istrip] = 0.;
        ana->chi2Probab_[iapv][istrip] = 0.;
        ana->residualRMS_[iapv][istrip] = 0.;
        ana->residualSigmaGaus_[iapv][istrip] = 0.;
        ana->noiseSignificance_[iapv][istrip] = 0.;
        ana->residualMean_[iapv][istrip] = 0.;
        ana->residualSkewness_[iapv][istrip] = 0.;
        ana->residualKurtosis_[iapv][istrip] = 0.;
        ana->residualIntegralNsigma_[iapv][istrip] = 0.;
        ana->residualIntegral_[iapv][istrip] = 0.;
        ana->deadStripBit_[iapv][istrip] = 0;
        ana->badStripBit_[iapv][istrip] = 0;
      }
    }
  }
  
  // -----------------------------------------------------------------------------
  //
  void PedsFullNoiseAlgorithm::analyse() {
    // check base analysis object
    if (!anal()) {
      edm::LogWarning(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
                                        << " NULL pointer to base Analysis object!";
      return;
    }
  
    CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
    PedsFullNoiseAnalysis* ana = dynamic_cast<PedsFullNoiseAnalysis*>(tmp);
  
    // check PedsFullNoiseAnalysis object
    if (!ana) {
      edm::LogWarning(mlCommissioning_) << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
                                        << " NULL pointer to derived Analysis object!";
      return;
    }
  
    // check if the histograms exists
    if (!hPeds_.first) {
      ana->addErrorCode(sistrip::nullPtr_);
      return;
    }
  
    if (!hNoise_.first) {
      ana->addErrorCode(sistrip::nullPtr_);
      return;
    }
  
    if (!hNoise2D_.first) {
      ana->addErrorCode(sistrip::nullPtr_);
      return;
    }
  
    // take the histograms
    TProfile* histoPeds = dynamic_cast<TProfile*>(hPeds_.first);
    TProfile* histoNoiseMean = dynamic_cast<TProfile*>(hNoise_.first);
    TH2S* histoNoise = dynamic_cast<TH2S*>(hNoise2D_.first);
  
    // Make sanity checks about pointers
    if (not histoPeds) {
      ana->addErrorCode(sistrip::nullPtr_);
      return;
    }
  
    if (not histoNoiseMean) {
      ana->addErrorCode(sistrip::nullPtr_);
      return;
    }
  
    if (not histoNoise) {
      ana->addErrorCode(sistrip::nullPtr_);
      return;
    }
  
    // check the binning  --> each x-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
    if (histoPeds->GetNbinsX() != 256) {
      ana->addErrorCode(sistrip::numberOfBins_);
      return;
    }
  
    //check the binning  --> each x-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
    if (histoNoiseMean->GetNbinsX() != 256) {
      ana->addErrorCode(sistrip::numberOfBins_);
      return;
    }
  
    //check the binning  --> each y-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
    if (histoNoise->GetNbinsY() != 256) {
      ana->addErrorCode(sistrip::numberOfBins_);
      return;
    }
  
    //Reset values
    reset(ana);
  
    // loop on each strip
    uint32_t apvID = -1;
  
    // Save basic information at strip / APV level
    vector<float> ped_max;
    vector<float> ped_min;
    vector<float> raw_max;
    vector<float> raw_min;
    vector<float> noise_max;
    vector<float> noise_min;
  
    // loop on each strip in the lldChannel
    for (int iStrip = 0; iStrip < histoPeds->GetNbinsX(); iStrip++) {
      if (iStrip < histoPeds->GetNbinsX() / 2)
        apvID = 0;
      else
        apvID = 1;
  
      int stripBin = 0;
      if (iStrip >= 128)
        stripBin = iStrip - 128;
      else
        stripBin = iStrip;
  
      ana->peds_[apvID][stripBin] = histoPeds->GetBinContent(iStrip + 1);        // pedestal value
      ana->noise_[apvID][stripBin] = histoNoiseMean->GetBinContent(iStrip + 1);  // noise value
      ana->raw_[apvID][stripBin] = histoPeds->GetBinError(iStrip + 1);           // raw noise value
  
      ana->pedsMean_[apvID] += ana->peds_[apvID][stripBin];    // mean pedestal
      ana->rawMean_[apvID] += ana->raw_[apvID][stripBin];      // mean raw noise
      ana->noiseMean_[apvID] += ana->noise_[apvID][stripBin];  // mean noise
  
      // max pedestal
      if (ped_max.size() < apvID + 1)
        ped_max.push_back(ana->peds_[apvID][stripBin]);
      else {
        if (ana->peds_[apvID][stripBin] > ped_max.at(apvID))
          ped_max.at(apvID) = ana->peds_[apvID][stripBin];
      }
  
      // min pedestal
      if (ped_min.size() < apvID + 1)
        ped_min.push_back(ana->peds_[apvID][stripBin]);
      else {
        if (ana->peds_[apvID][stripBin] < ped_min.at(apvID))
          ped_min.at(apvID) = ana->peds_[apvID][stripBin];  // min pedestal
      }
  
      // max noise
      if (noise_max.size() < apvID + 1)
        noise_max.push_back(ana->noise_[apvID][stripBin]);
      else {
        if (ana->noise_[apvID][stripBin] > noise_max.at(apvID))
          noise_max.at(apvID) = ana->noise_[apvID][stripBin];
      }
  
      // min noise
      if (noise_min.size() < apvID + 1)
        noise_min.push_back(ana->noise_[apvID][stripBin]);
      else {
        if (ana->noise_[apvID][stripBin] < noise_min.at(apvID))
          noise_min.at(apvID) = ana->noise_[apvID][stripBin];
      }
  
      // max raw
      if (raw_max.size() < apvID + 1)
        raw_max.push_back(ana->raw_[apvID][stripBin]);
      else {
        if (ana->raw_[apvID][stripBin] > raw_max.at(apvID))
          raw_max.at(apvID) = ana->raw_[apvID][stripBin];
      }
  
      // min raw
      if (raw_min.size() < apvID + 1)
        raw_min.push_back(ana->raw_[apvID][stripBin]);
      else {
        if (ana->raw_[apvID][stripBin] < raw_min.at(apvID))
          raw_min.at(apvID) = ana->raw_[apvID][stripBin];
      }
    }
  
    // Mean values
    for (unsigned int iApv = 0; iApv < ana->pedsMean_.size(); iApv++) {
      ana->pedsMean_.at(iApv) /= (ana->peds_[iApv].size());    // calculate mean pedestal per APV
      ana->rawMean_.at(iApv) /= (ana->raw_[iApv].size());      // calculate mean raw noise per APV
      ana->noiseMean_.at(iApv) /= (ana->noise_[iApv].size());  // calculate mean noise per APV
    }
  
    // Min and Max
    for (unsigned int iApv = 0; iApv < ped_max.size(); iApv++) {
      if (ped_max.at(iApv) > sistrip::maximum_)
        ana->pedsMax_.at(iApv) = sistrip::maximum_;
      else if (ped_max.at(iApv) < -1. * sistrip::maximum_)
        ana->pedsMax_.at(iApv) = -1. * sistrip::maximum_;
      else
        ana->pedsMax_.at(iApv) = ped_max.at(iApv);
  
      if (ped_min.at(iApv) > sistrip::maximum_)
        ana->pedsMin_.at(iApv) = sistrip::maximum_;
      else if (ped_min.at(iApv) < -1. * sistrip::maximum_)
        ana->pedsMin_.at(iApv) = -1. * sistrip::maximum_;
      else
        ana->pedsMin_.at(iApv) = ped_min.at(iApv);
  
      if (noise_max.at(iApv) > sistrip::maximum_)
        ana->noiseMax_.at(iApv) = sistrip::maximum_;
      else if (noise_max.at(iApv) < -1. * sistrip::maximum_)
        ana->noiseMax_.at(iApv) = -1. * sistrip::maximum_;
      else
        ana->noiseMax_.at(iApv) = noise_max.at(iApv);
  
      if (noise_min.at(iApv) > sistrip::maximum_)
        ana->noiseMin_.at(iApv) = sistrip::maximum_;
      else if (noise_min.at(iApv) < -1. * sistrip::maximum_)
        ana->noiseMin_.at(iApv) = -1. * sistrip::maximum_;
      else
        ana->noiseMin_.at(iApv) = noise_min.at(iApv);
  
      if (raw_max.at(iApv) > sistrip::maximum_)
        ana->rawMax_.at(iApv) = sistrip::maximum_;
      else if (raw_max.at(iApv) < -1. * sistrip::maximum_)
        ana->rawMax_.at(iApv) = -1. * sistrip::maximum_;
      else
        ana->rawMax_.at(iApv) = raw_max.at(iApv);
  
      if (raw_min.at(iApv) > sistrip::maximum_)
        ana->rawMin_.at(iApv) = sistrip::maximum_;
      else if (raw_min.at(iApv) < -1. * sistrip::maximum_)
        ana->rawMin_.at(iApv) = -1. * sistrip::maximum_;
      else
        ana->rawMin_.at(iApv) = raw_min.at(iApv);
    }
  
    // Calculate the spread for noise and pedestal
    apvID = -1;
  
    for (int iStrip = 0; iStrip < histoNoiseMean->GetNbinsX(); iStrip++) {
      if (iStrip < histoNoiseMean->GetNbinsX() / 2)
        apvID = 0;
      else
        apvID = 1;
      ana->pedsSpread_[apvID] += pow(histoPeds->GetBinContent(iStrip + 1) - ana->pedsMean_.at(apvID), 2);
      ana->noiseSpread_[apvID] += pow(histoNoiseMean->GetBinContent(iStrip + 1) - ana->noiseMean_.at(apvID), 2);
      ana->rawSpread_[apvID] += pow(histoPeds->GetBinError(iStrip + 1) - ana->rawMean_.at(apvID), 2);
    }
  
    for (unsigned int iApv = 0; iApv < ana->pedsSpread_.size(); iApv++) {
      ana->pedsSpread_[iApv] = sqrt(ana->pedsSpread_[iApv]) / sqrt(ana->peds_[iApv].size() - 1);
      ana->noiseSpread_[iApv] = sqrt(ana->noiseSpread_[iApv]) / sqrt(ana->noise_[iApv].size() - 1);
      ana->rawSpread_[iApv] = sqrt(ana->rawSpread_[iApv]) / sqrt(ana->raw_[iApv].size() - 1);
    }
  
    // loop on each strip in the lldChannel
    apvID = 0;
    TH1S* histoResidualStrip = new TH1S("histoResidualStrip",
                                        "",
                                        histoNoise->GetNbinsX(),
                                        histoNoise->GetXaxis()->GetXmin(),
                                        histoNoise->GetXaxis()->GetXmax());
    histoResidualStrip->Sumw2();
    histoResidualStrip->SetDirectory(nullptr);
    TF1* fitFunc = new TF1("fitFunc", "gaus(0)", histoNoise->GetXaxis()->GetXmin(), histoNoise->GetXaxis()->GetXmax());
    TF1* fit2Gaus = nullptr;
    TH1F* randomHisto = nullptr;
    TFitResultPtr result;
  
    for (int iStrip = 0; iStrip < histoNoise->GetNbinsY(); iStrip++) {
      // tell which APV
      if (iStrip < histoNoise->GetNbinsY() / 2)
        apvID = 0;
      else
        apvID = 1;
      histoResidualStrip->Reset();
  
      int stripBin = 0;
      if (iStrip >= 128)
        stripBin = iStrip - 128;
      else
        stripBin = iStrip;
  
      for (int iBinX = 0; iBinX < histoNoise->GetNbinsX(); iBinX++) {
        histoResidualStrip->SetBinContent(iBinX + 1, histoNoise->GetBinContent(iBinX + 1, iStrip + 1));
        histoResidualStrip->SetBinError(iBinX + 1, histoNoise->GetBinError(iBinX + 1, iStrip + 1));
      }
  
      if (histoResidualStrip->Integral() == 0) {  // dead strip --> i.e. no data
  
        // set default values
        ana->adProbab_[apvID][stripBin] = 0;
        ana->ksProbab_[apvID][stripBin] = 0;
        ana->jbProbab_[apvID][stripBin] = 0;
        ana->chi2Probab_[apvID][stripBin] = 0;
        ana->noiseSignificance_[apvID][stripBin] = 0;
        ana->residualMean_[apvID][stripBin] = 0;
        ana->residualRMS_[apvID][stripBin] = 0;
        ana->residualSigmaGaus_[apvID][stripBin] = 0;
        ana->residualSkewness_[apvID][stripBin] = 0;
        ana->residualKurtosis_[apvID][stripBin] = 0;
        ana->residualIntegralNsigma_[apvID][stripBin] = 0;
        ana->residualIntegral_[apvID][stripBin] = 0;
        ana->deadStrip_[apvID].push_back(stripBin);
        ana->deadStripBit_[apvID][stripBin] = 1;
        ana->badStripBit_[apvID][stripBin] = 0;
  
        SiStripFecKey fec_key(ana->fecKey());
        LogTrace(mlDqmClient_) << "DeadStrip: fecCrate "
                               << " " << fec_key.fecCrate() << " fecSlot " << fec_key.fecSlot() << " fecRing "
                               << fec_key.fecRing() << " ccuAddr " << fec_key.ccuAddr() << " ccChan  "
                               << fec_key.ccuChan() << " lldChan " << fec_key.lldChan() << " apvID " << apvID
                               << " stripID " << iStrip;
  
        continue;
      }
  
      // set / calculated basic quantities
      ana->residualMean_[apvID][stripBin] = histoResidualStrip->GetMean();
      ana->residualRMS_[apvID][stripBin] = histoResidualStrip->GetRMS();
      ana->residualSkewness_[apvID][stripBin] = histoResidualStrip->GetSkewness();
      ana->residualKurtosis_[apvID][stripBin] = histoResidualStrip->GetKurtosis();
      ana->noiseSignificance_[apvID][stripBin] =
          (ana->noise_[apvID][stripBin] - ana->noiseMean_[apvID]) / ana->noiseSpread_[apvID];
      ana->residualIntegral_[apvID][stripBin] = histoResidualStrip->Integral();
      ana->residualIntegralNsigma_[apvID][stripBin] =
          (histoResidualStrip->Integral(histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin] +
                                                                    ana->residualRMS_[apvID][stripBin] * integralNsigma_),
                                        histoResidualStrip->GetNbinsX() + 1) +
           histoResidualStrip->Integral(
               0,
               histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin] -
                                           ana->residualRMS_[apvID][stripBin] * integralNsigma_))) /
          ana->residualIntegral_[apvID][stripBin];
  
      // performing a Gaussian fit of the residual distribution
      fitFunc->SetRange(histoNoise->GetXaxis()->GetXmin(), histoNoise->GetXaxis()->GetXmax());
      fitFunc->SetParameters(ana->residualIntegral_[apvID][stripBin],
                             ana->residualMean_[apvID][stripBin],
                             ana->residualRMS_[apvID][stripBin]);
      result = histoResidualStrip->Fit(fitFunc, "QSRN");
  
      // Good gaussian fit
      if (result.Get()) {
        ana->residualSigmaGaus_[apvID][stripBin] = fitFunc->GetParameter(2);
        ana->chi2Probab_[apvID][stripBin] = result->Prob();
  
        // jacque bera probability
        float jbVal =
            (ana->residualIntegral_[apvID][stripBin] / 6) *
            (pow(ana->residualSkewness_[apvID][stripBin], 2) + pow(ana->residualKurtosis_[apvID][stripBin], 2) / 4);
        ana->jbProbab_[apvID][stripBin] = ROOT::Math::chisquared_cdf_c(jbVal, 2);
  
        //Kolmogorov Smirnov and Anderson Darlong
        if (randomHisto == nullptr)
          randomHisto = (TH1F*)histoResidualStrip->Clone("randomHisto");
        randomHisto->Reset();
        randomHisto->SetDirectory(nullptr);
  
        if (generateRandomHisto_) {  ///
          randomHisto->FillRandom("fitFunc", histoResidualStrip->Integral());
          if (randomHisto->Integral() != 0) {
            ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto, "N");
            ana->adProbab_[apvID][stripBin] = histoResidualStrip->AndersonDarlingTest(randomHisto);
          } else {
            ana->ksProbab_[apvID][stripBin] = 0;
            ana->adProbab_[apvID][stripBin] = 0;
          }
  
        } else {
          randomHisto->Add(fitFunc);
          if (randomHisto->Integral() != 0) {
            ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto, "N");
            ROOT::Fit::BinData data1;
            ROOT::Fit::BinData data2;
            ROOT::Fit::FillData(data1, histoResidualStrip, nullptr);
            data2.Initialize(randomHisto->GetNbinsX() + 1, 1);
            for (int ibin = 0; ibin < randomHisto->GetNbinsX(); ibin++) {
              if (histoResidualStrip->GetBinContent(ibin + 1) != 0 or randomHisto->GetBinContent(ibin + 1) >= 1)
                data2.Add(randomHisto->GetBinCenter(ibin + 1),
                          randomHisto->GetBinContent(ibin + 1),
                          randomHisto->GetBinError(ibin + 1));
            }
  
            double probab, value;
            ROOT::Math::GoFTest::AndersonDarling2SamplesTest(data1, data2, probab, value);
            ana->adProbab_[apvID][stripBin] = probab;
          } else {
            ana->ksProbab_[apvID][stripBin] = 0;
            ana->adProbab_[apvID][stripBin] = 0;
          }
        }
      }
  
      // start applying selections storing output
      bool badStripFlag = false;
      ana->deadStripBit_[apvID][stripBin] = 0;  // is not dead if the strip has data
  
      if (fabs(ana->residualMean_[apvID][stripBin]) > maxDriftResidualCut_ and not badStripFlag) {  //mean value
        ana->shiftedStrip_[apvID].push_back(stripBin);
        badStripFlag = true;
      }
  
      if (ana->residualRMS_[apvID][stripBin] < minStripNoiseCut_ and not badStripFlag) {  // low noise
        ana->lowNoiseStrip_[apvID].push_back(stripBin);
        badStripFlag = true;
      }
  
      if (ana->residualRMS_[apvID][stripBin] > maxStripNoiseCut_ and not badStripFlag) {  // large noise
        ana->largeNoiseStrip_[apvID].push_back(stripBin);
        badStripFlag = true;
      }
  
      if (fabs(ana->noiseSignificance_[apvID][stripBin]) > maxStripNoiseSignificanceCut_ and
          not badStripFlag) {  // large noise significance
        ana->largeNoiseSignificance_[apvID].push_back(stripBin);
        badStripFlag = true;
      }
  
      if (result.Get() and result->Status() != 0)  // bad fit status
        ana->badFitStatus_[apvID].push_back(stripBin);
  
      if (ana->adProbab_[apvID][stripBin] < adProbabCut_ and not badStripFlag)  // bad AD p-value --> store the strip-id
        ana->badADProbab_[apvID].push_back(stripBin);
  
      if (ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and not badStripFlag)  // bad KS p-value --> store the strip-id
        ana->badKSProbab_[apvID].push_back(stripBin);
  
      if (ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and not badStripFlag)  // bad JB p-value  --> store the strip-id
        ana->badJBProbab_[apvID].push_back(stripBin);
  
      if (ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_ and
          not badStripFlag)  // bad CHI2 p-value  --> store the strip-id
        ana->badChi2Probab_[apvID].push_back(stripBin);
  
      if (ana->adProbab_[apvID][stripBin] < adProbabCut_ and ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and
          ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_)
        badStripFlag = true;  // bad strip is flagged as bad by all the methods
  
      if (ana->residualKurtosis_[apvID][stripBin] > kurtosisCut_ and
          ana->residualIntegralNsigma_[apvID][stripBin] > integralTailCut_ and not badStripFlag) {  // bad tails
        ana->badTailStrip_[apvID].push_back(stripBin);
        badStripFlag = true;
      }
  
      if (badStripFlag) {  // loop for double peaked
  
        fit2Gaus = new TF1("dgaus",
                           "[0]*exp(-((x-[1])*(x-[1]))/(2*[2]*[2]))+[3]*exp(-((x-[4])*(x-[4]))/(2*[5]*[5]))",
                           histoNoise->GetXaxis()->GetXmin(),
                           histoNoise->GetXaxis()->GetXmax());
        fit2Gaus->SetParameter(0, fitFunc->GetParameter(0) / 2);
        fit2Gaus->SetParameter(3, fitFunc->GetParameter(0) / 2);
        fit2Gaus->SetParameter(1, 1.);
        fit2Gaus->SetParameter(4, -1.);
        fit2Gaus->SetParameter(2, fitFunc->GetParameter(2));
        fit2Gaus->SetParameter(5, fitFunc->GetParameter(2));
        fit2Gaus->SetParLimits(1, 0., histoNoise->GetXaxis()->GetXmax());
        fit2Gaus->SetParLimits(4, histoNoise->GetXaxis()->GetXmin(), 0);
        result = histoResidualStrip->Fit(fit2Gaus, "QSR");
  
        // ashman distance
        float ashman = TMath::Power(2, 0.5) * abs(fit2Gaus->GetParameter(1) - fit2Gaus->GetParameter(4)) /
                       (sqrt(pow(fit2Gaus->GetParameter(2), 2) + pow(fit2Gaus->GetParameter(5), 2)));
        // amplitude
        float amplitudeRatio = std::min(fit2Gaus->GetParameter(0), fit2Gaus->GetParameter(3)) /
                               std::max(fit2Gaus->GetParameter(0), fit2Gaus->GetParameter(3));
  
        if (ashman > ashmanDistance_ and amplitudeRatio > amplitudeRatio_)
          ana->badDoublePeakStrip_[apvID].push_back(stripBin);
      }
  
      if (badStripFlag) {  // save the final bit
  
        ana->badStrip_[apvID].push_back(stripBin);
        ana->badStripBit_[apvID][stripBin] = 1;
  
        SiStripFecKey fec_key(ana->fecKey());
        LogTrace(mlDqmClient_) << "BadStrip: fecCrate "
                               << " " << fec_key.fecCrate() << " fecSlot " << fec_key.fecSlot() << " fecRing "
                               << fec_key.fecRing() << " ccuAddr " << fec_key.ccuAddr() << " ccChan  "
                               << fec_key.ccuChan() << " lldChan " << fec_key.lldChan() << " apvID " << apvID
                               << " stripID " << stripBin;
  
      } else
        ana->badStripBit_[apvID][stripBin] = 0;
    }
  
    ped_max.clear();
    ped_min.clear();
    raw_max.clear();
    raw_min.clear();
    noise_max.clear();
    noise_min.clear();
    if (histoResidualStrip)
      delete histoResidualStrip;
    if (fitFunc)
      delete fitFunc;
    if (randomHisto)
      delete randomHisto;
    if (fit2Gaus)
      delete fit2Gaus;
  }
  