SiStripHotStripAlgorithmFromClusterOccupancy.cc

Go to the documentation of this file.

#include "CalibTracker/SiStripQuality/interface/SiStripHotStripAlgorithmFromClusterOccupancy.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"

SiStripHotStripAlgorithmFromClusterOccupancy::SiStripHotStripAlgorithmFromClusterOccupancy(
    const edm::ParameterSet& iConfig, const TrackerTopology* theTopo)
    : prob_(1.E-7),
      ratio_(1.5),
      MinNumEntries_(0),
      MinNumEntriesPerStrip_(0),
      Nevents_(0),
      occupancy_(0),
      OutFileName_("Occupancy.root"),
      tTopo(theTopo),
      UseInputDB_(iConfig.getUntrackedParameter<bool>("UseInputDB", false)) {
  minNevents_ = Nevents_ * occupancy_;
}

SiStripHotStripAlgorithmFromClusterOccupancy::~SiStripHotStripAlgorithmFromClusterOccupancy() {
  LogTrace("SiStripHotStripAlgorithmFromClusterOccupancy")
      << "[SiStripHotStripAlgorithmFromClusterOccupancy::~SiStripHotStripAlgorithmFromClusterOccupancy] " << std::endl;
}

void SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips(SiStripQuality* OutSiStripQuality,
                                                                    HistoMap& DM,
                                                                    const SiStripQuality* InSiStripQuality) {
  LogTrace("SiStripHotStripAlgorithmFromClusterOccupancy")
      << "[SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips] " << std::endl;

  if (WriteOutputFile_ == true) {
    f = new TFile(OutFileName_.c_str(), "RECREATE");
    f->cd();

    striptree = new TTree("stripOccupancy", "tree");

    striptree->Branch("DetRawId", &detrawid, "DetRawId/I");
    striptree->Branch("SubDetId", &subdetid, "SubDetId/I");
    striptree->Branch("Layer_Ring", &layer_ring, "Layer_Ring/I");
    striptree->Branch("Disc", &disc, "Disc/I");
    striptree->Branch("IsBack", &isback, "IsBack/I");
    striptree->Branch("IsExternalString", &isexternalstring, "IsExternalString/I");
    striptree->Branch("IsZMinusSide", &iszminusside, "IsZMinusSide/I");
    striptree->Branch("RodStringPetal", &rodstringpetal, "RodStringPetal/I");
    striptree->Branch("IsStereo", &isstereo, "IsStereo/I");
    striptree->Branch("ModulePosition", &module_position, "ModulePosition/I");
    striptree->Branch("NumberOfStrips", &number_strips, "NumberOfStrips/I");
    striptree->Branch("StripNumber", &strip_number, "StripNumber/I");
    striptree->Branch("APVChannel", &apv_channel, "APVChannel/I");
    striptree->Branch("StripGlobalPositionX", &global_position_x, "StripGlobalPositionX/F");
    striptree->Branch("StripGlobalPositionY", &global_position_y, "StripGlobalPositionY/F");
    striptree->Branch("StripGlobalPositionZ", &global_position_z, "StripGlobalPositionZ/F");
    striptree->Branch("IsHot", &isHot, "IsHot/I");
    striptree->Branch("HotStripsPerAPV", &hotStripsPerAPV, "HotStripsPerAPV/I");
    striptree->Branch("HotStripsPerModule", &hotStripsPerModule, "HotStripsPerModule/I");
    striptree->Branch("StripOccupancy", &stripOccupancy, "StripOccupancy/D");
    striptree->Branch("StripHits", &stripHits, "StripHits/I");
    striptree->Branch("PoissonProb", &poissonProb, "PoissonProb/D");
    striptree->Branch("MedianAPVHits", &medianAPVHits, "MedianAPVHits/D");
    striptree->Branch("AvgAPVHits", &avgAPVHits, "AvgAPVHits/D");
  }

  HistoMap::iterator it = DM.begin();
  HistoMap::iterator itEnd = DM.end();
  std::vector<unsigned int> badStripList;
  uint32_t detid;
  for (; it != itEnd; ++it) {
    pHisto phisto;
    detid = it->first;

    DetId detectorId = DetId(detid);
    phisto._SubdetId = detectorId.subdetId();

    if (edm::isDebugEnabled())
      LogTrace("SiStripHotStrip") << "Analyzing detid " << detid << std::endl;

    int numberAPVs = (int)(it->second.get())->GetNbinsX() / 128;

    // Set the values for the tree:

    detrawid = detid;
    subdetid = detectorId.subdetId();
    number_strips = (int)(it->second.get())->GetNbinsX();
    switch (detectorId.subdetId()) {
      case StripSubdetector::TIB:
        layer_ring = tTopo->tibLayer(detrawid);
        disc = -1;
        isstereo = tTopo->tibIsStereo(detrawid);
        isback = -1;
        if (tTopo->tibIsExternalString(detrawid))
          isexternalstring = 1;
        else
          isexternalstring = 0;
        if (tTopo->tibIsZMinusSide(detrawid))
          iszminusside = 1;
        else
          iszminusside = 0;
        rodstringpetal = tTopo->tibString(detrawid);
        module_position = tTopo->tibModule(detrawid);
        break;

      case StripSubdetector::TID:
        layer_ring = tTopo->tidRing(detrawid);
        disc = tTopo->tidWheel(detrawid);
        isstereo = tTopo->tidIsStereo(detrawid);
        if (tTopo->tidIsBackRing(detrawid))
          isback = 1;
        else
          isback = 0;
        if (tTopo->tidIsZMinusSide(detrawid))
          iszminusside = 1;
        else
          iszminusside = 0;
        isexternalstring = -1;
        rodstringpetal = -1;
        module_position = tTopo->tidModule(detrawid);
        break;

      case StripSubdetector::TOB:
        layer_ring = tTopo->tobLayer(detrawid);
        disc = -1;
        isstereo = tTopo->tobIsStereo(detrawid);
        isback = -1;
        if (tTopo->tobIsZMinusSide(detrawid))
          iszminusside = 1;
        else
          iszminusside = 0;
        isexternalstring = -1;
        rodstringpetal = tTopo->tobRod(detrawid);
        module_position = tTopo->tobModule(detrawid);
        break;

      case StripSubdetector::TEC:
        layer_ring = tTopo->tecRing(detrawid);
        disc = tTopo->tecWheel(detrawid);
        isstereo = tTopo->tecIsStereo(detrawid);
        if (tTopo->tecIsBackPetal(detrawid))
          isback = 1;
        else
          isback = 0;
        if (tTopo->tecIsZMinusSide(detrawid))
          iszminusside = 1;
        else
          iszminusside = 0;
        isexternalstring = -1;
        rodstringpetal = tTopo->tecPetalNumber(detrawid);
        module_position = tTopo->tecModule(detrawid);
        break;

      default:
        std::cout << "### Detector does not belong to TIB, TID, TOB or TEC !? ###" << std::endl;
        std::cout << "### DetRawId: " << detrawid << " ###" << std::endl;
    }

    // End: Set the values for the tree.

    pQuality = OutSiStripQuality;
    badStripList.clear();

    for (int i = 0; i < 768; i++) {
      ishot[i] = 0;
      stripoccupancy[i] = 0;
      striphits[i] = 0;
      poissonprob[i] = 0;
      hotstripsperapv[i / 128] = 0;
    }

    hotstripspermodule = 0;

    for (int apv = 0; apv < numberAPVs; apv++) {
      if (UseInputDB_) {
        if (InSiStripQuality->IsApvBad(detid, apv)) {
          if (edm::isDebugEnabled())
            LogTrace("SiStripHotStrip") << "(Module and Apv number) " << detid << " , " << apv
                                        << " excluded by input ESetup." << std::endl;
          continue;  //if the apv is already flagged as bad, continue.
        } else {
          if (edm::isDebugEnabled())
            LogTrace("SiStripHotStrip") << "(Module and Apv number) " << detid << " , " << apv
                                        << " good by input ESetup." << std::endl;
        }
      }

      phisto._th1f = new TH1F("tmp", "tmp", 128, 0.5, 128.5);
      int NumberEntriesPerAPV = 0;

      for (int strip = 0; strip < 128; strip++) {
        phisto._th1f->SetBinContent(strip + 1, (it->second.get())->GetBinContent((apv * 128) + strip + 1));
        NumberEntriesPerAPV += (int)(it->second.get())->GetBinContent((apv * 128) + strip + 1);
      }

      phisto._th1f->SetEntries(NumberEntriesPerAPV);
      phisto._NEntries = (int)phisto._th1f->GetEntries();
      phisto._NEmptyBins = 0;

      LogTrace("SiStripHotStrip") << "Number of clusters in APV " << apv << ": " << NumberEntriesPerAPV << std::endl;

      iterativeSearch(phisto, badStripList, apv);

      delete phisto._th1f;
    }

    const StripGeomDetUnit* theStripDet = dynamic_cast<const StripGeomDetUnit*>((TkGeom->idToDet(detectorId)));
    const StripTopology* theStripTopol = dynamic_cast<const StripTopology*>(&(theStripDet->specificTopology()));

    for (int strip = 0; strip < number_strips; strip++) {
      strip_number = strip + 1;
      apv_channel = (strip % 128) + 1;
      isHot = ishot[strip];
      stripOccupancy = stripoccupancy[strip];
      stripHits = striphits[strip];
      poissonProb = poissonprob[strip];
      medianAPVHits = medianapvhits[strip / 128];
      avgAPVHits = avgapvhits[strip / 128];

      hotStripsPerModule = hotstripspermodule;
      hotStripsPerAPV = hotstripsperapv[strip / 128];

      LocalPoint pos_strip_local = theStripTopol->localPosition(strip);
      GlobalPoint pos_strip_global = (TkGeom->idToDet(detectorId))->surface().toGlobal(pos_strip_local);

      global_position_x = pos_strip_global.x();
      global_position_y = pos_strip_global.y();
      global_position_z = pos_strip_global.z();

      if (WriteOutputFile_ == true)
        striptree->Fill();
    }

    if (badStripList.begin() == badStripList.end())
      continue;

    OutSiStripQuality->compact(detid, badStripList);

    SiStripQuality::Range range(badStripList.begin(), badStripList.end());
    if (!OutSiStripQuality->put(detid, range))
      edm::LogError("SiStripHotStrip")
          << "[SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips] detid already exists" << std::endl;
  }
  OutSiStripQuality->fillBadComponents();

  if (WriteOutputFile_ == true) {
    f->cd();
    striptree->Write();
    f->Close();
  }

  LogTrace("SiStripHotStrip") << ss.str() << std::endl;
}

void SiStripHotStripAlgorithmFromClusterOccupancy::iterativeSearch(pHisto& histo,
                                                                   std::vector<unsigned int>& vect,
                                                                   int apv) {
  if (!histo._NEntries || histo._NEntries <= MinNumEntries_ || histo._NEntries <= minNevents_)
    return;

  size_t startingSize = vect.size();
  long double diff = 1. - prob_;

  size_t Nbins = histo._th1f->GetNbinsX();
  size_t ibinStart = 1;
  size_t ibinStop = Nbins + 1;
  int MaxEntry = (int)histo._th1f->GetMaximum();

  std::vector<long double> vPoissonProbs(MaxEntry + 1, 0);
  long double meanVal = 1. * histo._NEntries / (1. * Nbins - histo._NEmptyBins);
  evaluatePoissonian(vPoissonProbs, meanVal);

  // Find median occupancy, taking into account only good strips
  unsigned int goodstripentries[128];
  int nGoodStrips = 0;
  for (size_t i = ibinStart; i < ibinStop; ++i) {
    if (ishot[(apv * 128) + i - 1] == 0) {
      goodstripentries[nGoodStrips] = (unsigned int)histo._th1f->GetBinContent(i);
      nGoodStrips++;
    }
  }
  double median = TMath::Median(nGoodStrips, goodstripentries);

  for (size_t i = ibinStart; i < ibinStop; ++i) {
    unsigned int entries = (unsigned int)histo._th1f->GetBinContent(i);

    if (ishot[(apv * 128) + i - 1] == 0) {
      stripoccupancy[(apv * 128) + i - 1] = entries / (double)Nevents_;
      striphits[(apv * 128) + i - 1] = entries;
      poissonprob[(apv * 128) + i - 1] = 1 - vPoissonProbs[entries];
      medianapvhits[apv] = median;
      avgapvhits[apv] = meanVal;
    }
    if (entries <= MinNumEntriesPerStrip_ || entries <= minNevents_ || entries / median < ratio_)
      continue;

    if (diff < vPoissonProbs[entries]) {
      ishot[(apv * 128) + i - 1] = 1;
      hotstripspermodule++;
      hotstripsperapv[apv]++;
      histo._th1f->SetBinContent(i, 0.);
      histo._NEntries -= entries;
      histo._NEmptyBins++;
      if (edm::isDebugEnabled())
        LogTrace("SiStripHotStrip") << " rejecting strip " << (apv * 128) + i - 1 << " value " << entries << " diff  "
                                    << diff << " prob " << vPoissonProbs[entries] << std::endl;
      vect.push_back(pQuality->encode((apv * 128) + i - 1, 1, 0));
    }
  }
  if (edm::isDebugEnabled())
    LogTrace("SiStripHotStrip") << " [SiStripHotStripAlgorithmFromClusterOccupancy::iterativeSearch] Nbins=" << Nbins
                                << " MaxEntry=" << MaxEntry << " meanVal=" << meanVal
                                << " NEmptyBins=" << histo._NEmptyBins << " NEntries=" << histo._NEntries
                                << " thEntries " << histo._th1f->GetEntries() << " startingSize " << startingSize
                                << " vector.size " << vect.size() << std::endl;

  if (vect.size() != startingSize)
    iterativeSearch(histo, vect, apv);
}

void SiStripHotStripAlgorithmFromClusterOccupancy::evaluatePoissonian(std::vector<long double>& vPoissonProbs,
                                                                      long double& meanVal) {
  for (size_t i = 0; i < vPoissonProbs.size(); ++i) {
    vPoissonProbs[i] = (i == 0) ? TMath::Poisson(i, meanVal) : vPoissonProbs[i - 1] + TMath::Poisson(i, meanVal);
  }
}

void SiStripHotStripAlgorithmFromClusterOccupancy::setNumberOfEvents(double Nevents) {
  Nevents_ = Nevents;
  minNevents_ = occupancy_ * Nevents_;
  if (edm::isDebugEnabled())
    LogTrace("SiStripHotStrip") << " [SiStripHotStripAlgorithmFromClusterOccupancy::setNumberOfEvents] minNumber of "
                                   "Events per strip used to consider a strip bad"
                                << minNevents_ << " for occupancy " << occupancy_ << std::endl;
}