00001 #include "CalibTracker/SiStripQuality/interface/SiStripHotStripAlgorithmFromClusterOccupancy.h"
00002 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
00003 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00004 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
00005 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00006
00007
00008
00009
00010 SiStripHotStripAlgorithmFromClusterOccupancy::SiStripHotStripAlgorithmFromClusterOccupancy(const edm::ParameterSet& iConfig, const TrackerTopology* theTopo):
00011 prob_(1.E-7),
00012 MinNumEntries_(0),
00013 MinNumEntriesPerStrip_(0),
00014 Nevents_(0),
00015 occupancy_(0),
00016 OutFileName_("Occupancy.root"),
00017 tTopo(theTopo),
00018 UseInputDB_(iConfig.getUntrackedParameter<bool>("UseInputDB",false))
00019 {
00020 minNevents_=Nevents_*occupancy_;
00021 }
00022
00023
00024 SiStripHotStripAlgorithmFromClusterOccupancy::~SiStripHotStripAlgorithmFromClusterOccupancy(){
00025 LogTrace("SiStripHotStripAlgorithmFromClusterOccupancy")<<"[SiStripHotStripAlgorithmFromClusterOccupancy::~SiStripHotStripAlgorithmFromClusterOccupancy] "<<std::endl;
00026 }
00027
00028 void SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips(SiStripQuality* OutSiStripQuality,HistoMap& DM,edm::ESHandle<SiStripQuality>& InSiStripQuality){
00029
00030 LogTrace("SiStripHotStripAlgorithmFromClusterOccupancy")<<"[SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips] "<<std::endl;
00031
00032
00033
00034 if (WriteOutputFile_==true){
00035 f = new TFile(OutFileName_.c_str(),"RECREATE");
00036 f->cd();
00037
00038 striptree = new TTree("stripOccupancy","tree");
00039
00040 striptree->Branch("DetRawId", &detrawid, "DetRawId/I");
00041 striptree->Branch("SubDetId", &subdetid, "SubDetId/I");
00042 striptree->Branch("Layer_Ring", &layer_ring, "Layer_Ring/I");
00043 striptree->Branch("Disc", &disc, "Disc/I");
00044 striptree->Branch("IsBack", &isback, "IsBack/I");
00045 striptree->Branch("IsExternalString", &isexternalstring, "IsExternalString/I");
00046 striptree->Branch("IsZMinusSide", &iszminusside, "IsZMinusSide/I");
00047 striptree->Branch("RodStringPetal", &rodstringpetal, "RodStringPetal/I");
00048 striptree->Branch("IsStereo", &isstereo, "IsStereo/I");
00049 striptree->Branch("ModulePosition", &module_position, "ModulePosition/I");
00050 striptree->Branch("NumberOfStrips", &number_strips, "NumberOfStrips/I");
00051 striptree->Branch("StripNumber", &strip_number, "StripNumber/I");
00052 striptree->Branch("APVChannel", &apv_channel, "APVChannel/I");
00053 striptree->Branch("StripGlobalPositionX", &global_position_x, "StripGlobalPositionX/F");
00054 striptree->Branch("StripGlobalPositionY", &global_position_y, "StripGlobalPositionY/F");
00055 striptree->Branch("StripGlobalPositionZ", &global_position_z, "StripGlobalPositionZ/F");
00056 striptree->Branch("IsHot", &isHot, "IsHot/I");
00057 striptree->Branch("HotStripsPerAPV", &hotStripsPerAPV, "HotStripsPerAPV/I");
00058 striptree->Branch("HotStripsPerModule", &hotStripsPerModule,"HotStripsPerModule/I");
00059 striptree->Branch("StripOccupancy", &stripOccupancy, "StripOccupancy/D");
00060 striptree->Branch("StripHits", &stripHits, "StripHits/I");
00061 striptree->Branch("PoissonProb", &poissonProb, "PoissonProb/D");
00062 }
00063
00064
00065 HistoMap::iterator it=DM.begin();
00066 HistoMap::iterator itEnd=DM.end();
00067 std::vector<unsigned int> badStripList;
00068 uint32_t detid;
00069 for (;it!=itEnd;++it){
00070 pHisto phisto;
00071 detid=it->first;
00072
00073 DetId detectorId=DetId(detid);
00074 phisto._SubdetId=detectorId.subdetId();
00075
00076 if (edm::isDebugEnabled())
00077 LogTrace("SiStripHotStrip") << "Analyzing detid " << detid<< std::endl;
00078
00079 int numberAPVs = (int)(it->second.get())->GetNbinsX()/128;
00080
00081
00082
00083 detrawid = detid;
00084 subdetid = detectorId.subdetId();
00085 number_strips = (int)(it->second.get())->GetNbinsX();
00086 if (SiStripDetId(detrawid).stereo() !=0 ) isstereo = 1;
00087 else isstereo = 0;
00088 switch (detectorId.subdetId())
00089 {
00090 case StripSubdetector::TIB :
00091 layer_ring = tTopo->tibLayer(detrawid);
00092 disc = -1;
00093 isback = -1;
00094 if (tTopo->tibIsExternalString(detrawid)) isexternalstring = 1;
00095 else isexternalstring = 0;
00096 if (tTopo->tibIsZMinusSide(detrawid)) iszminusside = 1;
00097 else iszminusside = 0;
00098 rodstringpetal = tTopo->tibString(detrawid);
00099 module_position = tTopo->tibModule(detrawid);
00100 break;
00101
00102 case StripSubdetector::TID :
00103 layer_ring = tTopo->tidRing(detrawid);
00104 disc = tTopo->tidWheel(detrawid);
00105 if (tTopo->tidIsBackRing(detrawid)) isback = 1;
00106 else isback = 0;
00107 if (tTopo->tidIsZMinusSide(detrawid)) iszminusside = 1;
00108 else iszminusside = 0;
00109 isexternalstring = -1;
00110 rodstringpetal = -1;
00111 module_position = tTopo->tidModule(detrawid);
00112 break;
00113
00114 case StripSubdetector::TOB :
00115 layer_ring = tTopo->tobLayer(detrawid);
00116 disc = -1;
00117 isback = -1;
00118 if (tTopo->tobIsZMinusSide(detrawid)) iszminusside = 1;
00119 else iszminusside = 0;
00120 isexternalstring = -1;
00121 rodstringpetal = tTopo->tobRod(detrawid);
00122 module_position = tTopo->tobModule(detrawid);
00123 break;
00124
00125 case StripSubdetector::TEC :
00126 layer_ring = tTopo->tecRing(detrawid);
00127 disc = tTopo->tecWheel(detrawid);
00128 if (tTopo->tecIsBackPetal(detrawid)) isback = 1;
00129 else isback = 0;
00130 if (tTopo->tecIsZMinusSide(detrawid)) iszminusside = 1;
00131 else iszminusside = 0;
00132 isexternalstring = -1;
00133 rodstringpetal = tTopo->tecPetalNumber(detrawid);
00134 module_position = tTopo->tecModule(detrawid);
00135 break;
00136
00137 default :
00138 std::cout << "### Detector does not belong to TIB, TID, TOB or TEC !? ###" << std::endl;
00139 std::cout << "### DetRawId: " << detrawid << " ###" << std::endl;
00140 }
00141
00142
00143
00144
00145 pQuality=OutSiStripQuality;
00146 badStripList.clear();
00147
00148 for (int i=0; i<768; i++){
00149 ishot[i] = 0;
00150 stripoccupancy[i] = 0;
00151 striphits[i] = 0;
00152 poissonprob[i] = 0;
00153 hotstripsperapv[i/128] = 0;
00154 }
00155
00156 hotstripspermodule = 0;
00157
00158 for (int apv=0; apv<numberAPVs; apv++){
00159 if(UseInputDB_){
00160 if(InSiStripQuality->IsApvBad(detid,apv) ){
00161 if(edm::isDebugEnabled())
00162 LogTrace("SiStripHotStrip")<<"(Module and Apv number) "<<detid<<" , "<<apv<<" excluded by input ESetup."<<std::endl;
00163 continue;
00164 }
00165 else{
00166 if(edm::isDebugEnabled())
00167 LogTrace("SiStripHotStrip")<<"(Module and Apv number) "<<detid<<" , "<<apv<<" good by input ESetup."<<std::endl;
00168 }
00169 }
00170
00171 phisto._th1f = new TH1F("tmp","tmp",128,0.5,128.5);
00172 int NumberEntriesPerAPV=0;
00173
00174 for (int strip=0; strip<128; strip++){
00175 phisto._th1f->SetBinContent(strip+1,(it->second.get())->GetBinContent((apv*128)+strip+1));
00176 NumberEntriesPerAPV += (int)(it->second.get())->GetBinContent((apv*128)+strip+1);
00177 }
00178
00179 phisto._th1f->SetEntries(NumberEntriesPerAPV);
00180 phisto._NEntries=(int)phisto._th1f->GetEntries();
00181
00182 LogTrace("SiStripHotStrip") << "Number of clusters in APV " << apv << ": " << NumberEntriesPerAPV << std::endl;
00183
00184 iterativeSearch(phisto,badStripList,apv);
00185
00186 delete phisto._th1f;
00187 }
00188
00189 const StripGeomDetUnit* theStripDet = dynamic_cast<const StripGeomDetUnit*>( (TkGeom->idToDet(detectorId)) );
00190 const StripTopology* theStripTopol = dynamic_cast<const StripTopology*>( &(theStripDet->specificTopology()) );
00191
00192 for (int strip=0; strip<number_strips; strip++)
00193 {
00194 strip_number = strip+1;
00195 apv_channel = (strip%128)+1;
00196 isHot = ishot[strip];
00197 stripOccupancy = stripoccupancy[strip];
00198 stripHits = striphits[strip];
00199 poissonProb = poissonprob[strip];
00200
00201 hotStripsPerModule = hotstripspermodule;
00202 hotStripsPerAPV = hotstripsperapv[strip/128];
00203
00204 LocalPoint pos_strip_local = theStripTopol->localPosition(strip);
00205 GlobalPoint pos_strip_global = (TkGeom->idToDet(detectorId))->surface().toGlobal(pos_strip_local);
00206
00207 global_position_x = pos_strip_global.x();
00208 global_position_y = pos_strip_global.y();
00209 global_position_z = pos_strip_global.z();
00210
00211 if (WriteOutputFile_==true) striptree->Fill();
00212 }
00213
00214 if (badStripList.begin()==badStripList.end())
00215 continue;
00216
00217 OutSiStripQuality->compact(detid,badStripList);
00218
00219 SiStripQuality::Range range(badStripList.begin(),badStripList.end());
00220 if ( ! OutSiStripQuality->put(detid,range) )
00221 edm::LogError("SiStripHotStrip")<<"[SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips] detid already exists"<<std::endl;
00222 }
00223 OutSiStripQuality->fillBadComponents();
00224
00225 if (WriteOutputFile_==true){
00226 f->cd();
00227 striptree->Write();
00228 f->Close();
00229 }
00230
00231 LogTrace("SiStripHotStrip") << ss.str() << std::endl;
00232 }
00233
00234
00235 void SiStripHotStripAlgorithmFromClusterOccupancy::iterativeSearch(pHisto& histo,std::vector<unsigned int>& vect, int apv){
00236 if (!histo._NEntries || histo._NEntries <=MinNumEntries_ || histo._NEntries <= minNevents_)
00237 return;
00238
00239 size_t startingSize=vect.size();
00240 long double diff=1.-prob_;
00241
00242 size_t Nbins = histo._th1f->GetNbinsX();
00243 size_t ibinStart = 1;
00244 size_t ibinStop = Nbins+1;
00245 int MaxEntry = (int)histo._th1f->GetMaximum();
00246
00247 std::vector<long double> vPoissonProbs(MaxEntry+1,0);
00248 long double meanVal=1.*histo._NEntries/(1.*Nbins-histo._NEmptyBins);
00249 evaluatePoissonian(vPoissonProbs,meanVal);
00250
00251 for (size_t i=ibinStart; i<ibinStop; ++i){
00252 unsigned int entries= (unsigned int)histo._th1f->GetBinContent(i);
00253
00254 if (ishot[(apv*128)+i-1]==0){
00255 stripoccupancy[(apv*128)+i-1] = entries/(double) Nevents_;
00256 striphits[(apv*128)+i-1] = entries;
00257 poissonprob[(apv*128)+i-1] = 1-vPoissonProbs[entries];
00258 }
00259
00260 if (entries<=MinNumEntriesPerStrip_ || entries <= minNevents_)
00261 continue;
00262
00263 if(diff<vPoissonProbs[entries]){
00264 ishot[(apv*128)+i-1] = 1;
00265 hotstripspermodule++;
00266 hotstripsperapv[apv]++;
00267 histo._th1f->SetBinContent(i,0.);
00268 histo._NEntries-=entries;
00269 histo._NEmptyBins++;
00270 if (edm::isDebugEnabled())
00271 LogTrace("SiStripHotStrip")<< " rejecting strip " << (apv*128)+i-1 << " value " << entries << " diff " << diff << " prob " << vPoissonProbs[entries]<< std::endl;
00272 vect.push_back(pQuality->encode((apv*128)+i-1,1,0));
00273 }
00274
00275 }
00276 if (edm::isDebugEnabled())
00277 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;
00278
00279 if (vect.size()!=startingSize)
00280 iterativeSearch(histo,vect,apv);
00281 }
00282
00283 void SiStripHotStripAlgorithmFromClusterOccupancy::evaluatePoissonian(std::vector<long double>& vPoissonProbs, long double& meanVal){
00284 for(size_t i=0;i<vPoissonProbs.size();++i){
00285 vPoissonProbs[i]= (i==0)?TMath::Poisson(i,meanVal):vPoissonProbs[i-1]+TMath::Poisson(i,meanVal);
00286 }
00287 }
00288
00289 void SiStripHotStripAlgorithmFromClusterOccupancy::setNumberOfEvents(double Nevents){
00290 Nevents_=Nevents;
00291 minNevents_=occupancy_*Nevents_;
00292 if (edm::isDebugEnabled())
00293 LogTrace("SiStripHotStrip")<<" [SiStripHotStripAlgorithmFromClusterOccupancy::setNumberOfEvents] minNumber of Events per strip used to consider a strip bad" << minNevents_ << " for occupancy " << occupancy_ << std::endl;
00294 }
1.7.3