adding NoisyStripFinder package

Tracker/TrackerRecAlgs/NoisyStripFinder/CMakeLists.txt

+################################################################################
+# Package: NoisyStripFinder
+################################################################################
+
+# Declare the package name:
+atlas_subdir( NoisyStripFinder )
+
+# Component(s) in the package:
+atlas_add_component( NoisyStripFinder
+                     src/*.cxx src/*.h
+                     src/components/*.cxx
+                     LINK_LIBRARIES AthenaBaseComps StoreGateLib SGtests Identifier GaudiKernel TrackerRawData TrackerPrepRawData FaserDetDescr TrackerIdentifier 
+                                    TrackerReadoutGeometry AthViews FaserSCT_ConditionsData xAODFaserTrigger )
+
+# Install files from the package:
+#atlas_install_headers( TrackerPrepRawDataFormation )
+
+atlas_install_python_modules( python/*.py )
+
+#atlas_install_scripts( test/*.py )
+

Tracker/TrackerRecAlgs/NoisyStripFinder/python/NoisyStripFinderConfig.py

+"""
+Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+"""
+
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+
+
+def NoisyStripFinderBasicCfg(flags, **kwargs):
+    acc = FaserSCT_GeometryCfg(flags)
+    kwargs.setdefault("DataObjectName", "SCT_RDOs")
+    kwargs.setdefault("OutputHistRootName", "NoisyStripFinderHist.root")
+    kwargs.setdefault("TriggerMask", 0x10)
+    Tracker__NoisyStripFinder = CompFactory.Tracker.NoisyStripFinder
+    acc.addEventAlgo(Tracker__NoisyStripFinder(**kwargs))
+    return acc
+
+
+def NoisyStripFinderCfg(flags, **kwargs):
+    acc = NoisyStripFinderBasicCfg(flags, **kwargs)
+    return acc
+

Tracker/TrackerRecAlgs/NoisyStripFinder/share/NoisyStripFinderHist_Analysis.py

+#!/usr/bin/env python
+
+# This script combines the strip occupancy histograms from NoisyStripFinderHist.root files then it generates a pdf with sumarry plots and an XML file with a list of noisy strips
+
+# Set up (Py)ROOT.
+import ROOT
+import xml.etree.ElementTree as ET
+
+ROOT.gStyle.SetOptStat(0) #take away option box in histograms
+ROOT.gStyle.SetOptTitle(1)
+#ROOT.gROOT.SetStyle("ATLAS")
+
+def GetKeyNames( self ): # gets a list of object names in a root file
+    return [key.GetName() for key in f.GetListOfKeys()]
+ROOT.TFile.GetKeyNames = GetKeyNames
+
+numEvents = 0
+nfiles = 0
+HistDict = {}
+
+ROOT.TH1.AddDirectory(0) # This is necessary in order to have the histogram data after closing the file
+
+for inputfile in ["./NoisyStripFinderHist.root"]:
+    print("opening root file ",inputfile)
+    f = ROOT.TFile.Open(inputfile, "r")
+    numEvents += f.Get("numEvents").GetVal()
+    if nfiles == 0:
+        trigger = f.Get("trigger").GetVal()
+    for rootkey in f.GetKeyNames():
+        if rootkey == 'numEvents' or rootkey == 'trigger':
+            continue # skip over the root objects TParameters that store the trigger and number of events data
+        if rootkey in HistDict: # if sensor histogram has already been stored, then add to it 
+            HistDict[rootkey].Add(f.Get(rootkey),1.0)
+        else: # if sensor histogram has not already been stored, then store this histogram
+            HistDict[rootkey] = f.Get(rootkey).Clone()
+
+    nfiles += 1
+
+    # Good to cleanup
+    f.Close()
+
+print("Total # of root files analyzed = ", nfiles)
+print("Trigger mask = ", trigger)
+print("Total number of events = ", numEvents)
+
+# Now make some plots
+filename = "NoisyStripFinderHist_Analysis.pdf"
+
+c = ROOT.TCanvas()
+c.Print(filename+'[')
+
+for dictkey in HistDict:
+    c = ROOT.TCanvas()
+    c.Divide(1,2)
+    c.cd(1)
+    tempHist = HistDict[dictkey].Clone()
+    tempHist.Draw("")
+    c.cd(2)
+    HistDict[dictkey].Scale(1.0/float(numEvents)) # Normalize histrograms by total event number in order to get occupancies
+    HistDict[dictkey].GetYaxis().SetTitle("Strip Occupancy")
+    HistDict[dictkey].Draw("")
+#    ROOT.gPad.SetLogy()
+    c.Print(filename)
+
+c.Print(filename+']') # Must close file at the end
+
+# Now make an XML file to store noisy strip data
+root = ET.Element('NoisyStripData')
+for dictkey in HistDict:
+    sensorhash = ET.SubElement(root, "Sensor_Hash")
+    sensorhash.text = dictkey
+    bini = 1
+    while bini <= 768:
+        if HistDict[dictkey].GetBinContent(bini) >= 0.1 :
+            strip = ET.SubElement(sensorhash, "Strip")
+            strip.text = str(bini - 1) # strip number is offset by histogram bin number by 1 because of underflow bin
+            occupancy = ET.SubElement(strip, "Occupancy")
+            occupancy.text = str(HistDict[dictkey].GetBinContent(bini))
+        bini += 1           
+ 
+tree = ET.ElementTree(root)
+tree.write("NoisyPixels.xml")
+

Tracker/TrackerRecAlgs/NoisyStripFinder/src/NoisyStripFinder.cxx

+/*
+   Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+   */
+
+#include "NoisyStripFinder.h"
+
+#include "FaserDetDescr/FaserDetectorID.h"    
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerRawData/FaserSCT_RDO_Container.h"
+#include "TrackerRawData/FaserSCT_RDORawData.h"
+#include "StoreGate/WriteHandle.h"
+
+#include <sstream>
+#include <string.h>
+#include <list>
+#include <TH1.h>
+#include <TFile.h>
+#include <TParameter.h>
+
+namespace Tracker
+{
+// Constructor with parameters:
+NoisyStripFinder::NoisyStripFinder(const std::string& name, ISvcLocator* pSvcLocator) :
+  AthReentrantAlgorithm(name, pSvcLocator),
+  m_idHelper{nullptr}
+{
+}
+
+// Initialize method:
+StatusCode NoisyStripFinder::initialize() {
+  ATH_MSG_INFO("NoisyStripFinder::initialize()!");
+
+  ATH_CHECK(m_rdoContainerKey.initialize());
+  ATH_CHECK(m_FaserTriggerData.initialize());
+
+  // Get the SCT ID helper
+  ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+
+  return StatusCode::SUCCESS;
+}
+
+// Execute method:
+StatusCode NoisyStripFinder::execute(const EventContext& ctx) const {
+
+  SG::ReadHandle<xAOD::FaserTriggerData> xaod(m_FaserTriggerData, ctx);
+
+  int trig_int = xaod->tap();
+  int trigmask_int = m_triggerMask.value();
+
+  ATH_MSG_DEBUG("Found FaserTriggerData");
+  ATH_MSG_DEBUG("trigger = " << trig_int);
+  ATH_MSG_DEBUG("mask = " << trigmask_int);
+
+  if (!(xaod->tap() & m_triggerMask.value())) return StatusCode::SUCCESS; // only process events that pass the trigger mask
+
+  ATH_MSG_INFO("trigger passed mask");
+  ++m_numberOfEvents;                                                                
+
+  // First, we have to retrieve and access the container, not because we want to 
+  // use it, but in order to generate the proxies for the collections, if they 
+  // are being provided by a container converter.
+  SG::ReadHandle<FaserSCT_RDO_Container> rdoContainer{m_rdoContainerKey, ctx};
+  ATH_CHECK(rdoContainer.isValid());
+
+  // Anything to dereference the DataHandle will trigger the converter
+  FaserSCT_RDO_Container::const_iterator rdoCollections{rdoContainer->begin()};
+  FaserSCT_RDO_Container::const_iterator rdoCollectionsEnd{rdoContainer->end()};
+  for (; rdoCollections != rdoCollectionsEnd; ++rdoCollections) {
+    ++m_numberOfRDOCollection;
+    const TrackerRawDataCollection<FaserSCT_RDORawData>* rd{*rdoCollections};
+    ATH_MSG_DEBUG("RDO collection size=" << rd->size() << ", Hash=" << rd->identifyHash());
+    for (const FaserSCT_RDORawData* pRawData: *rd) {
+      ++m_numberOfRDO;
+      const Identifier      firstStripId(pRawData->identify());
+      const int             thisStrip(m_idHelper->strip(firstStripId));
+      ATH_MSG_DEBUG( "---------- new rdo ----------------------------- " );
+      ATH_MSG_DEBUG( "rd->identifyHash() = " << rd->identifyHash() );
+      ATH_MSG_DEBUG( "strip = " << thisStrip );
+
+      if ( NoisyStrip_histmap.count(rd->identifyHash()) == 0 ) { // map.count(key) returns 1 if the key exists and 0 if it does not
+        std::string histname_str = std::to_string(rd->identifyHash());
+        char const *histname = histname_str.c_str();
+        std::string station_num = std::to_string(m_idHelper->station(firstStripId));
+        std::string layer_num = std::to_string(m_idHelper->layer(firstStripId));
+        std::string phi_num = std::to_string(m_idHelper->phi_module(firstStripId));
+        std::string eta_num = std::to_string(m_idHelper->eta_module(firstStripId));
+        std::string side_num = std::to_string(m_idHelper->side(firstStripId));
+        std::string hist_title_str = "Station " + station_num + " | Layer " + layer_num + " | phi-module " + phi_num + " | eta-module " + eta_num + " | Side " + side_num + "; strip number; # of events";
+        char const *hist_title = hist_title_str.c_str();
+        NoisyStrip_histmap[rd->identifyHash()] = new TH1D(histname, hist_title, 768, -0.5, 767.5);
+      } 
+      NoisyStrip_histmap[rd->identifyHash()]->Fill(thisStrip); // increase the value by one so as to count the number of times the strip was active    
+    }
+  }
+  return StatusCode::SUCCESS;
+}
+
+// Finalize method:
+StatusCode NoisyStripFinder::finalize() 
+{
+  ATH_MSG_INFO("NoisyStripFinder::finalize()");
+  ATH_MSG_INFO( m_numberOfEvents << " events processed" );
+  ATH_MSG_INFO( m_numberOfRDOCollection << " RDO collections processed" );
+  ATH_MSG_INFO( m_numberOfRDO<< " RawData" );
+  ATH_MSG_INFO( "Number of sensors found = " << NoisyStrip_histmap.size() << " out of 144" );
+
+  for (int ihash = 0; ihash < 144; ++ihash){ // print out the sensors that are missing 
+    if ( NoisyStrip_histmap.count(ihash) == 0 ){
+      ATH_MSG_INFO("missing sensor # " << ihash);
+    }
+  }
+
+  const char *outputname = m_OutputRootName.value().c_str();
+
+  TFile* outputfile = new TFile(outputname,"RECREATE");
+
+  int trigmask_int = m_triggerMask.value();
+
+  TParameter("numEvents", m_numberOfEvents).Write();
+  TParameter("trigger", trigmask_int).Write();
+
+  std::map<int,TH1D*>::iterator it = NoisyStrip_histmap.begin();
+  // Iterate over the map using Iterator till end.
+  while (it != NoisyStrip_histmap.end()){
+    ATH_MSG_INFO( "" );
+    ATH_MSG_INFO( "---------- hot strip occupancy >= 0.1 for Tracker Sensor hash = "<< it->first <<" ----------" );
+    int i = 1;
+    while (i <= 768){
+      if ( it->second->GetBinContent(i)/(double)m_numberOfEvents >= 0.1 ){
+        ATH_MSG_INFO( "hot strip # = " << i-1 << ", hit occupancy = " << it->second->GetBinContent(i)/(double)m_numberOfEvents ); // print out hot strips
+      }
+      i++;
+    }
+    it->second->Write();
+    it++;
+  }
+
+  outputfile->Close();
+
+  return StatusCode::SUCCESS;
+}
+}
+

Tracker/TrackerRecAlgs/NoisyStripFinder/src/NoisyStripFinder.h

+/*
+   Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+   */
+
+#ifndef NoisyStripFinder_H
+#define NoisyStripFinder_H
+
+// Base class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+////Next contains a typedef so cannot be fwd declared
+#include "TrackerRawData/FaserSCT_RDO_Container.h"
+#include "xAODFaserTrigger/FaserTriggerData.h"
+
+//Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "StoreGate/ReadHandleKey.h"
+
+
+//STL
+#include <map>
+#include <TH1.h>
+#include <TFile.h>
+#include <TParameter.h>
+
+#include <string>
+
+class FaserSCT_ID;
+class ISvcLocator;
+class StatusCode;
+
+namespace Tracker
+{
+
+/**
+ *    @class NoisyStripFinder
+ *    @brief Creates histograms with strip occupancy data from SCT Raw Data Objects
+ *    Creates histograms with strip occupancy data from SCT Raw Data Objects. Root files containing strip occupancy histograms can then be combined and analyzed (example pyROOT script in share folder) in order to make an XML database of noisy strips.
+ */
+class NoisyStripFinder : public AthReentrantAlgorithm {
+  public:
+    /// Constructor with parameters:
+    NoisyStripFinder(const std::string& name, ISvcLocator* pSvcLocator);
+
+    /**    @name Usual algorithm methods */
+    //@{
+    ///Retrieve the tools used and initialize variables
+    virtual StatusCode initialize() override;
+    ///Form clusters and record them in StoreGate (detector store)
+    virtual StatusCode execute(const EventContext& ctx) const override;
+    ///Clean up and release the collection containers
+    virtual StatusCode finalize() override;
+    //Make this algorithm clonable.
+    virtual bool isClonable() const override { return true; };
+    //@}
+
+  private:
+    /**    @name Disallow default instantiation, copy, assignment */
+    //@{
+    //NoisyStripFinder() = delete;
+    //NoisyStripFinder(const NoisyStripFinder&) = delete;
+    //NoisyStripFinder &operator=(const NoisyStripFinder&) = delete;
+    //@}
+
+    StringProperty m_OutputRootName{this, "OutputHistRootName", "NoisyStripFinderHist.root", "Name of output histogram root file for NoisyStripFinder"};
+
+    UnsignedIntegerProperty m_triggerMask{this, "TriggerMask", 0x10, "Trigger mask (0x10 = random trig)"};
+
+    const FaserSCT_ID* m_idHelper;
+
+    SG::ReadHandleKey<FaserSCT_RDO_Container> m_rdoContainerKey{this, "DataObjectName", "FaserSCT_RDOs", "FaserSCT RDOs"};
+    SG::ReadHandleKey<xAOD::FaserTriggerData> m_FaserTriggerData{ this, "FaserTriggerDataKey", "FaserTriggerData", "ReadHandleKey for xAOD::FaserTriggerData"};
+
+    mutable int m_numberOfEvents{0};
+    mutable std::atomic<int> m_numberOfRDOCollection{0};
+    mutable std::atomic<int> m_numberOfRDO{0};
+
+    mutable std::map<int,TH1D*> NoisyStrip_histmap;
+
+};
+}
+#endif // NoisyStripFinder_H
+

Tracker/TrackerRecAlgs/NoisyStripFinder/src/components/NoisyStripFinder_entries.cxx

+#include "../NoisyStripFinder.h"
+
+DECLARE_COMPONENT( Tracker::NoisyStripFinder )
+