Updating Hot Strip algorithms in CaliTracker/SiStripQuality putting e…

…xtra requirement on the strip Digi wrt APV median as developed by Mario Galanti + Message logger ourput made more logical in DQMOffline/CalibTracker/plugins/SiStripQualityHotStripIdentifierRoot.cc

CalibTracker/SiStripQuality/interface/SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy.h

@@ -93,6 +93,7 @@ class SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy{
   void fillStripDQMHistograms();
 
   long double prob_;
+  long double ratio_;
   unsigned short MinNumEntries_;
   unsigned short MinNumEntriesPerStrip_;
   long double lowoccupancy_;
@@ -193,6 +194,8 @@ class SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy{
   int hotStripsPerModule;
   double singleStripOccupancy;
   int stripHits;
+  double medianAPVHits;
+  double avgAPVHits;
   double poissonProb;
 
   int ishot[128];
@@ -201,6 +204,8 @@ class SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy{
   double stripoccupancy[128];
   int striphits[128];
   double poissonprob[128];
+  double medianapvhits[6];
+  double avgapvhits[6];
 
   std::stringstream ss;
 

CalibTracker/SiStripQuality/interface/SiStripHotStripAlgorithmFromClusterOccupancy.h

@@ -88,6 +88,7 @@ class SiStripHotStripAlgorithmFromClusterOccupancy{
   void evaluatePoissonian(std::vector<long double>& , long double& meanVal);
 
   long double prob_;
+  long double ratio_;
   unsigned short MinNumEntries_;
   unsigned short MinNumEntriesPerStrip_;
   double Nevents_;
@@ -126,6 +127,8 @@ class SiStripHotStripAlgorithmFromClusterOccupancy{
   int hotStripsPerModule;
   double stripOccupancy;
   int stripHits;
+  double medianAPVHits;
+  double avgAPVHits;
   double poissonProb;
 
   int ishot[768];
@@ -134,6 +137,8 @@ class SiStripHotStripAlgorithmFromClusterOccupancy{
   double stripoccupancy[768];
   int striphits[768];
   double poissonprob[768];
+  double medianapvhits[6];
+  double avgapvhits[6];
 
 
   std::stringstream ss;   

CalibTracker/SiStripQuality/src/SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy.cc

@@ -4,12 +4,14 @@
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "DataFormats/DetId/interface/DetId.h"
 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
 #include "Geometry/CommonTopologies/interface/StripTopology.h"
 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
 
 
 SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy(const edm::ParameterSet& iConfig, const TrackerTopology* theTopo):
+  ratio_(1.5),
   lowoccupancy_(0),
   highoccupancy_(100),
   absolutelow_(0),
@@ -82,6 +84,8 @@ void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::extractBadAPVSandStr
       striptree->Branch("StripOccupancy",       &singleStripOccupancy, "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();
@@ -658,17 +662,29 @@ void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::iterativeSearch(Apv&
   long double meanVal=1.*histo.NEntries[apv]/(1.*Nbins-histo.NEmptyBins[apv]); 
   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[apv]->GetBinContent(i);
+      nGoodStrips++;
+    }
+  }
+  double median = TMath::Median(nGoodStrips,goodstripentries);
+
   for (size_t i=ibinStart; i<ibinStop; ++i){
     unsigned int entries= (unsigned int)histo.th1f[apv]->GetBinContent(i);
 
     if (ishot[i-1]==0){
       stripoccupancy[i-1] = entries/(double) Nevents_;
       striphits[i-1]      = entries;
       poissonprob[i-1]    = 1-vPoissonProbs[entries];
+      medianapvhits[apv]  = median;
+      avgapvhits[apv] = meanVal;
     }
 
-    if (entries<=MinNumEntriesPerStrip_ || entries <= minNevents_)
-      continue;
+    if (entries<=MinNumEntriesPerStrip_ || entries <= minNevents_ || entries / median < ratio_) continue;
 
     if(diff<vPoissonProbs[entries]){
       ishot[i-1] = 1;

CalibTracker/SiStripQuality/src/SiStripHotStripAlgorithmFromClusterOccupancy.cc

@@ -1,11 +1,13 @@
 #include "CalibTracker/SiStripQuality/interface/SiStripHotStripAlgorithmFromClusterOccupancy.h"
 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
 
 
 
 
 SiStripHotStripAlgorithmFromClusterOccupancy::SiStripHotStripAlgorithmFromClusterOccupancy(const edm::ParameterSet& iConfig, const TrackerTopology* theTopo):
     prob_(1.E-7),
+    ratio_(1.5),
     MinNumEntries_(0),
     MinNumEntriesPerStrip_(0),
     Nevents_(0),
@@ -56,7 +58,9 @@ void SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips(SiStripQuali
   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();
@@ -197,6 +201,8 @@ void SiStripHotStripAlgorithmFromClusterOccupancy::extractBadStrips(SiStripQuali
 	stripOccupancy = stripoccupancy[strip];
 	stripHits      = striphits[strip];
 	poissonProb    = poissonprob[strip];
+	medianAPVHits  = medianapvhits[strip/128];
+	avgAPVHits     = avgapvhits[strip/128];
 
 	hotStripsPerModule = hotstripspermodule;
 	hotStripsPerAPV    = hotstripsperapv[strip/128];
@@ -248,17 +254,28 @@ void SiStripHotStripAlgorithmFromClusterOccupancy::iterativeSearch(pHisto& histo
   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_)
-      continue;
+    if (entries<=MinNumEntriesPerStrip_ || entries <= minNevents_ || entries / median < ratio_) continue;
 
     if(diff<vPoissonProbs[entries]){
       ishot[(apv*128)+i-1] = 1;

DQMOffline/CalibTracker/plugins/SiStripQualityHotStripIdentifierRoot.cc

@@ -244,6 +244,7 @@ void SiStripQualityHotStripIdentifierRoot::bookHistos(){
 						      << " occup " << parameters.getUntrackedParameter<double>("OccupancyThreshold",0)
 						      << " OccupancyHisto" << parameters.getUntrackedParameter<std::string>("OccupancyHisto") << std::endl;
 
+  // Check Number of Events
   for (; iter!=iterEnd;++iter) {
     std::string me_name = (*iter)->getName();
 
@@ -255,10 +256,15 @@ void SiStripQualityHotStripIdentifierRoot::bookHistos(){
 
       gotNentries=true;
       edm::LogInfo("SiStripQualityHotStripIdentifierRoot")<< "[SiStripQualityHotStripIdentifierRoot::bookHistos]  gotNentries flag " << gotNentries << std::endl;
-    } else {
-      edm::LogWarning("SiStripQualityHotStripIdentifierRoot") <<" [SiStripQualityHotStripIdentifierRoot::bookHistos] :: Histogram with to check # of evemnts missing" <<std::endl;
+      break;
     }
-
+  }
+  if (!gotNentries) {
+    edm::LogWarning("SiStripQualityHotStripIdentifierRoot") <<" [SiStripQualityHotStripIdentifierRoot::bookHistos] :: Histogram with to check # of evemnts missing" <<std::endl;
+  }
+  for (; iter!=iterEnd;++iter) {
+    std::string me_name = (*iter)->getName();
+
     if (strstr(me_name.c_str(),(parameters.getUntrackedParameter<std::string>("OccupancyHisto")).c_str())==NULL)
       continue;
 
@@ -280,7 +286,6 @@ void SiStripQualityHotStripIdentifierRoot::bookHistos(){
     ClusterPositionHistoMap[detid]=boost::shared_ptr<TH1F>(new TH1F(*(*iter)->getTH1F()));
 
   }
-  if(!gotNentries) edm::LogWarning("MissingNumberOfEvents") <<"Missing histogram to get the number of events";
 
 }