File indexing completed on 2023-03-17 10:44:44
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
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
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;
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
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 }
