Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:15:00

0001 #include "CalibTracker/SiStripQuality/interface/SiStripHotStripAlgorithmFromClusterOccupancy.h"
0002 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0003 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0004 
0005 SiStripHotStripAlgorithmFromClusterOccupancy::SiStripHotStripAlgorithmFromClusterOccupancy(
0006     const edm::ParameterSet& iConfig, const TrackerTopology* theTopo)
0007     : prob_(1.E-7),
0008       ratio_(1.5),
0009       MinNumEntries_(0),
0010       MinNumEntriesPerStrip_(0),
0011       Nevents_(0),
0012       occupancy_(0),
0013       OutFileName_("Occupancy.root"),
0014       tTopo(theTopo),
0015       UseInputDB_(iConfig.getUntrackedParameter<bool>("UseInputDB", false)) {
0016   minNevents_ = Nevents_ * occupancy_;
0017 }
0018 
0019 SiStripHotStripAlgorithmFromClusterOccupancy::~SiStripHotStripAlgorithmFromClusterOccupancy() {
0020   LogTrace("SiStripHotStripAlgorithmFromClusterOccupancy")
0021       << "[SiStripHotStripAlgorithmFromClusterOccupancy::~SiStripHotStripAlgorithmFromClusterOccupancy] " << std::endl;
0022 }
0023 
0024 void SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips(SiStripQuality* OutSiStripQuality,
0025                                                                     HistoMap& DM,
0026                                                                     const SiStripQuality* InSiStripQuality) {
0027   LogTrace("SiStripHotStripAlgorithmFromClusterOccupancy")
0028       << "[SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips] " << std::endl;
0029 
0030   if (WriteOutputFile_ == true) {
0031     f = new TFile(OutFileName_.c_str(), "RECREATE");
0032     f->cd();
0033 
0034     striptree = new TTree("stripOccupancy", "tree");
0035 
0036     striptree->Branch("DetRawId", &detrawid, "DetRawId/I");
0037     striptree->Branch("SubDetId", &subdetid, "SubDetId/I");
0038     striptree->Branch("Layer_Ring", &layer_ring, "Layer_Ring/I");
0039     striptree->Branch("Disc", &disc, "Disc/I");
0040     striptree->Branch("IsBack", &isback, "IsBack/I");
0041     striptree->Branch("IsExternalString", &isexternalstring, "IsExternalString/I");
0042     striptree->Branch("IsZMinusSide", &iszminusside, "IsZMinusSide/I");
0043     striptree->Branch("RodStringPetal", &rodstringpetal, "RodStringPetal/I");
0044     striptree->Branch("IsStereo", &isstereo, "IsStereo/I");
0045     striptree->Branch("ModulePosition", &module_position, "ModulePosition/I");
0046     striptree->Branch("NumberOfStrips", &number_strips, "NumberOfStrips/I");
0047     striptree->Branch("StripNumber", &strip_number, "StripNumber/I");
0048     striptree->Branch("APVChannel", &apv_channel, "APVChannel/I");
0049     striptree->Branch("StripGlobalPositionX", &global_position_x, "StripGlobalPositionX/F");
0050     striptree->Branch("StripGlobalPositionY", &global_position_y, "StripGlobalPositionY/F");
0051     striptree->Branch("StripGlobalPositionZ", &global_position_z, "StripGlobalPositionZ/F");
0052     striptree->Branch("IsHot", &isHot, "IsHot/I");
0053     striptree->Branch("HotStripsPerAPV", &hotStripsPerAPV, "HotStripsPerAPV/I");
0054     striptree->Branch("HotStripsPerModule", &hotStripsPerModule, "HotStripsPerModule/I");
0055     striptree->Branch("StripOccupancy", &stripOccupancy, "StripOccupancy/D");
0056     striptree->Branch("StripHits", &stripHits, "StripHits/I");
0057     striptree->Branch("PoissonProb", &poissonProb, "PoissonProb/D");
0058     striptree->Branch("MedianAPVHits", &medianAPVHits, "MedianAPVHits/D");
0059     striptree->Branch("AvgAPVHits", &avgAPVHits, "AvgAPVHits/D");
0060   }
0061 
0062   HistoMap::iterator it = DM.begin();
0063   HistoMap::iterator itEnd = DM.end();
0064   std::vector<unsigned int> badStripList;
0065   uint32_t detid;
0066   for (; it != itEnd; ++it) {
0067     pHisto phisto;
0068     detid = it->first;
0069 
0070     DetId detectorId = DetId(detid);
0071     phisto._SubdetId = detectorId.subdetId();
0072 
0073     if (edm::isDebugEnabled())
0074       LogTrace("SiStripHotStrip") << "Analyzing detid " << detid << std::endl;
0075 
0076     int numberAPVs = (int)(it->second.get())->GetNbinsX() / 128;
0077 
0078     // Set the values for the tree:
0079 
0080     detrawid = detid;
0081     subdetid = detectorId.subdetId();
0082     number_strips = (int)(it->second.get())->GetNbinsX();
0083     switch (detectorId.subdetId()) {
0084       case StripSubdetector::TIB:
0085         layer_ring = tTopo->tibLayer(detrawid);
0086         disc = -1;
0087         isstereo = tTopo->tibIsStereo(detrawid);
0088         isback = -1;
0089         if (tTopo->tibIsExternalString(detrawid))
0090           isexternalstring = 1;
0091         else
0092           isexternalstring = 0;
0093         if (tTopo->tibIsZMinusSide(detrawid))
0094           iszminusside = 1;
0095         else
0096           iszminusside = 0;
0097         rodstringpetal = tTopo->tibString(detrawid);
0098         module_position = tTopo->tibModule(detrawid);
0099         break;
0100 
0101       case StripSubdetector::TID:
0102         layer_ring = tTopo->tidRing(detrawid);
0103         disc = tTopo->tidWheel(detrawid);
0104         isstereo = tTopo->tidIsStereo(detrawid);
0105         if (tTopo->tidIsBackRing(detrawid))
0106           isback = 1;
0107         else
0108           isback = 0;
0109         if (tTopo->tidIsZMinusSide(detrawid))
0110           iszminusside = 1;
0111         else
0112           iszminusside = 0;
0113         isexternalstring = -1;
0114         rodstringpetal = -1;
0115         module_position = tTopo->tidModule(detrawid);
0116         break;
0117 
0118       case StripSubdetector::TOB:
0119         layer_ring = tTopo->tobLayer(detrawid);
0120         disc = -1;
0121         isstereo = tTopo->tobIsStereo(detrawid);
0122         isback = -1;
0123         if (tTopo->tobIsZMinusSide(detrawid))
0124           iszminusside = 1;
0125         else
0126           iszminusside = 0;
0127         isexternalstring = -1;
0128         rodstringpetal = tTopo->tobRod(detrawid);
0129         module_position = tTopo->tobModule(detrawid);
0130         break;
0131 
0132       case StripSubdetector::TEC:
0133         layer_ring = tTopo->tecRing(detrawid);
0134         disc = tTopo->tecWheel(detrawid);
0135         isstereo = tTopo->tecIsStereo(detrawid);
0136         if (tTopo->tecIsBackPetal(detrawid))
0137           isback = 1;
0138         else
0139           isback = 0;
0140         if (tTopo->tecIsZMinusSide(detrawid))
0141           iszminusside = 1;
0142         else
0143           iszminusside = 0;
0144         isexternalstring = -1;
0145         rodstringpetal = tTopo->tecPetalNumber(detrawid);
0146         module_position = tTopo->tecModule(detrawid);
0147         break;
0148 
0149       default:
0150         std::cout << "### Detector does not belong to TIB, TID, TOB or TEC !? ###" << std::endl;
0151         std::cout << "### DetRawId: " << detrawid << " ###" << std::endl;
0152     }
0153 
0154     // End: Set the values for the tree.
0155 
0156     pQuality = OutSiStripQuality;
0157     badStripList.clear();
0158 
0159     for (int i = 0; i < 768; i++) {
0160       ishot[i] = 0;
0161       stripoccupancy[i] = 0;
0162       striphits[i] = 0;
0163       poissonprob[i] = 0;
0164       hotstripsperapv[i / 128] = 0;
0165     }
0166 
0167     hotstripspermodule = 0;
0168 
0169     for (int apv = 0; apv < numberAPVs; apv++) {
0170       if (UseInputDB_) {
0171         if (InSiStripQuality->IsApvBad(detid, apv)) {
0172           if (edm::isDebugEnabled())
0173             LogTrace("SiStripHotStrip") << "(Module and Apv number) " << detid << " , " << apv
0174                                         << " excluded by input ESetup." << std::endl;
0175           continue;  //if the apv is already flagged as bad, continue.
0176         } else {
0177           if (edm::isDebugEnabled())
0178             LogTrace("SiStripHotStrip") << "(Module and Apv number) " << detid << " , " << apv
0179                                         << " good by input ESetup." << std::endl;
0180         }
0181       }
0182 
0183       phisto._th1f = new TH1F("tmp", "tmp", 128, 0.5, 128.5);
0184       int NumberEntriesPerAPV = 0;
0185 
0186       for (int strip = 0; strip < 128; strip++) {
0187         phisto._th1f->SetBinContent(strip + 1, (it->second.get())->GetBinContent((apv * 128) + strip + 1));
0188         NumberEntriesPerAPV += (int)(it->second.get())->GetBinContent((apv * 128) + strip + 1);
0189       }
0190 
0191       phisto._th1f->SetEntries(NumberEntriesPerAPV);
0192       phisto._NEntries = (int)phisto._th1f->GetEntries();
0193       phisto._NEmptyBins = 0;
0194 
0195       LogTrace("SiStripHotStrip") << "Number of clusters in APV " << apv << ": " << NumberEntriesPerAPV << std::endl;
0196 
0197       iterativeSearch(phisto, badStripList, apv);
0198 
0199       delete phisto._th1f;
0200     }
0201 
0202     const StripGeomDetUnit* theStripDet = dynamic_cast<const StripGeomDetUnit*>((TkGeom->idToDet(detectorId)));
0203     const StripTopology* theStripTopol = dynamic_cast<const StripTopology*>(&(theStripDet->specificTopology()));
0204 
0205     for (int strip = 0; strip < number_strips; strip++) {
0206       strip_number = strip + 1;
0207       apv_channel = (strip % 128) + 1;
0208       isHot = ishot[strip];
0209       stripOccupancy = stripoccupancy[strip];
0210       stripHits = striphits[strip];
0211       poissonProb = poissonprob[strip];
0212       medianAPVHits = medianapvhits[strip / 128];
0213       avgAPVHits = avgapvhits[strip / 128];
0214 
0215       hotStripsPerModule = hotstripspermodule;
0216       hotStripsPerAPV = hotstripsperapv[strip / 128];
0217 
0218       LocalPoint pos_strip_local = theStripTopol->localPosition(strip);
0219       GlobalPoint pos_strip_global = (TkGeom->idToDet(detectorId))->surface().toGlobal(pos_strip_local);
0220 
0221       global_position_x = pos_strip_global.x();
0222       global_position_y = pos_strip_global.y();
0223       global_position_z = pos_strip_global.z();
0224 
0225       if (WriteOutputFile_ == true)
0226         striptree->Fill();
0227     }
0228 
0229     if (badStripList.begin() == badStripList.end())
0230       continue;
0231 
0232     OutSiStripQuality->compact(detid, badStripList);
0233 
0234     SiStripQuality::Range range(badStripList.begin(), badStripList.end());
0235     if (!OutSiStripQuality->put(detid, range))
0236       edm::LogError("SiStripHotStrip")
0237           << "[SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips] detid already exists" << std::endl;
0238   }
0239   OutSiStripQuality->fillBadComponents();
0240 
0241   if (WriteOutputFile_ == true) {
0242     f->cd();
0243     striptree->Write();
0244     f->Close();
0245   }
0246 
0247   LogTrace("SiStripHotStrip") << ss.str() << std::endl;
0248 }
0249 
0250 void SiStripHotStripAlgorithmFromClusterOccupancy::iterativeSearch(pHisto& histo,
0251                                                                    std::vector<unsigned int>& vect,
0252                                                                    int apv) {
0253   if (!histo._NEntries || histo._NEntries <= MinNumEntries_ || histo._NEntries <= minNevents_)
0254     return;
0255 
0256   size_t startingSize = vect.size();
0257   long double diff = 1. - prob_;
0258 
0259   size_t Nbins = histo._th1f->GetNbinsX();
0260   size_t ibinStart = 1;
0261   size_t ibinStop = Nbins + 1;
0262   int MaxEntry = (int)histo._th1f->GetMaximum();
0263 
0264   std::vector<long double> vPoissonProbs(MaxEntry + 1, 0);
0265   long double meanVal = 1. * histo._NEntries / (1. * Nbins - histo._NEmptyBins);
0266   evaluatePoissonian(vPoissonProbs, meanVal);
0267 
0268   // Find median occupancy, taking into account only good strips
0269   unsigned int goodstripentries[128];
0270   int nGoodStrips = 0;
0271   for (size_t i = ibinStart; i < ibinStop; ++i) {
0272     if (ishot[(apv * 128) + i - 1] == 0) {
0273       goodstripentries[nGoodStrips] = (unsigned int)histo._th1f->GetBinContent(i);
0274       nGoodStrips++;
0275     }
0276   }
0277   double median = TMath::Median(nGoodStrips, goodstripentries);
0278 
0279   for (size_t i = ibinStart; i < ibinStop; ++i) {
0280     unsigned int entries = (unsigned int)histo._th1f->GetBinContent(i);
0281 
0282     if (ishot[(apv * 128) + i - 1] == 0) {
0283       stripoccupancy[(apv * 128) + i - 1] = entries / (double)Nevents_;
0284       striphits[(apv * 128) + i - 1] = entries;
0285       poissonprob[(apv * 128) + i - 1] = 1 - vPoissonProbs[entries];
0286       medianapvhits[apv] = median;
0287       avgapvhits[apv] = meanVal;
0288     }
0289     if (entries <= MinNumEntriesPerStrip_ || entries <= minNevents_ || entries / median < ratio_)
0290       continue;
0291 
0292     if (diff < vPoissonProbs[entries]) {
0293       ishot[(apv * 128) + i - 1] = 1;
0294       hotstripspermodule++;
0295       hotstripsperapv[apv]++;
0296       histo._th1f->SetBinContent(i, 0.);
0297       histo._NEntries -= entries;
0298       histo._NEmptyBins++;
0299       if (edm::isDebugEnabled())
0300         LogTrace("SiStripHotStrip") << " rejecting strip " << (apv * 128) + i - 1 << " value " << entries << " diff  "
0301                                     << diff << " prob " << vPoissonProbs[entries] << std::endl;
0302       vect.push_back(pQuality->encode((apv * 128) + i - 1, 1, 0));
0303     }
0304   }
0305   if (edm::isDebugEnabled())
0306     LogTrace("SiStripHotStrip") << " [SiStripHotStripAlgorithmFromClusterOccupancy::iterativeSearch] Nbins=" << Nbins
0307                                 << " MaxEntry=" << MaxEntry << " meanVal=" << meanVal
0308                                 << " NEmptyBins=" << histo._NEmptyBins << " NEntries=" << histo._NEntries
0309                                 << " thEntries " << histo._th1f->GetEntries() << " startingSize " << startingSize
0310                                 << " vector.size " << vect.size() << std::endl;
0311 
0312   if (vect.size() != startingSize)
0313     iterativeSearch(histo, vect, apv);
0314 }
0315 
0316 void SiStripHotStripAlgorithmFromClusterOccupancy::evaluatePoissonian(std::vector<long double>& vPoissonProbs,
0317                                                                       long double& meanVal) {
0318   for (size_t i = 0; i < vPoissonProbs.size(); ++i) {
0319     vPoissonProbs[i] = (i == 0) ? TMath::Poisson(i, meanVal) : vPoissonProbs[i - 1] + TMath::Poisson(i, meanVal);
0320   }
0321 }
0322 
0323 void SiStripHotStripAlgorithmFromClusterOccupancy::setNumberOfEvents(double Nevents) {
0324   Nevents_ = Nevents;
0325   minNevents_ = occupancy_ * Nevents_;
0326   if (edm::isDebugEnabled())
0327     LogTrace("SiStripHotStrip") << " [SiStripHotStripAlgorithmFromClusterOccupancy::setNumberOfEvents] minNumber of "
0328                                    "Events per strip used to consider a strip bad"
0329                                 << minNevents_ << " for occupancy " << occupancy_ << std::endl;
0330 }