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