BeamBackgroundFilterAndQA/src/StreakSidebandFilter.cc

0001 /// ===========================================================================
0002 /*! \file    StreakSidebandFilter.cc
0003  *  \authors Hanpu Jiang, Derek Anderson
0004  *  \date    11.01.2024
0005  *
0006  *  Part of the BeamBackgroundFilterAndQA module, this
0007  *  filter returns true if event contains a "streaK"
0008  *  as defined by Hanpu Jiang's streak identification
0009  *  algorithm.
0010  */
0011 /// ===========================================================================
0012
0013 // module components
0014 #include "StreakSidebandFilter.h"
0015
0016 // calo base
0017 #include <calobase/TowerInfoContainer.h>
0018
0019 // phool libraries
0020 #include <phool/getClass.h>
0021
0022 // root libraries
0023 #include <TH1.h>
0024 #include <TH2.h>
0025
0026 // c++ utiilites
0027 #include <algorithm>
0028 #include <iostream>
0029 #include <map>
0030 #include <memory>
0031 #include <vector>
0032
0033 // ctor/dtor ==================================================================
0034
0035 // ----------------------------------------------------------------------------
0036 //! Default ctor
0037 // ----------------------------------------------------------------------------
0038 StreakSidebandFilter::StreakSidebandFilter(const std::string& name)
0039   : BaseBeamBackgroundFilter(name)
0040
0041 {
0042 }  // end ctor()
0043
0044 // ----------------------------------------------------------------------------
0045 //! ctor accepting config struct
0046 // ----------------------------------------------------------------------------
0047 StreakSidebandFilter::StreakSidebandFilter(const Config& config, const std::string& name)
0048   : BaseBeamBackgroundFilter(name)
0049   , m_config(config)
0050 {
0051 }  // end ctor(Config&)
0052
0053 // public methods =============================================================
0054
0055 // ----------------------------------------------------------------------------
0056 // Apply filter to check for beam background or not
0057 // ----------------------------------------------------------------------------
0058 bool StreakSidebandFilter::ApplyFilter(PHCompositeNode* topNode)
0059 {
0060   // print debug message
0061   if (m_config.debug && (m_config.verbosity > 2))
0062   {
0063     std::cout << "StreakSidebandFilter::ApplyFilter() Checking if streak found in OHCal via their sidebands" << std::endl;
0064   }
0065
0066   // grab input node
0067   GrabNodes(topNode);
0068
0069   // build tower map
0070   m_ohMap.Reset();
0071   m_ohMap.Build(m_ohContainer);
0072
0073   // reset number of streaky towers per eta bin
0074   std::fill(m_ohNumStreak.begin(), m_ohNumStreak.end(), 0);
0075
0076   // lambdas to get phi +- 1 neighbors
0077   auto getAdjacentUp = [this](const std::size_t phi)
0078   { return (phi + 1) % m_ohMap.towers.front().size(); };
0079   auto getAdjacentDown = [this](const std::size_t phi)
0080   { return (phi == 0) ? m_ohMap.towers.front().size() : (phi - 1); };
0081
0082   // loop over tower (eta, phi) map to find streaks
0083   for (std::size_t iPhi = 0; iPhi < m_ohMap.towers.front().size(); ++iPhi)
0084   {
0085     for (std::size_t iEta = 0; iEta < m_ohMap.towers.size(); ++iEta)
0086     {
0087       // check if tower is a candidate for being in a streak
0088       const bool isNotStreak = IsTowerNotStreaky(m_ohMap.towers[iEta][iPhi]);
0089       if (isNotStreak)
0090       {
0091         continue;
0092       }
0093
0094       // grab adjacent towers
0095       const std::size_t iUp = getAdjacentUp(iPhi);
0096       const std::size_t iDown = getAdjacentDown(iPhi);
0097
0098       // and check if adjacent towers consistent w/ a streak
0099       const bool isUpNotStreak = IsNeighborNotStreaky(m_ohMap.towers[iEta][iUp]);
0100       const bool isDownNotStreak = IsNeighborNotStreaky(m_ohMap.towers[iEta][iDown]);
0101       if (isUpNotStreak || isDownNotStreak)
0102       {
0103         continue;
0104       }
0105
0106       // finally, increment no. of streaky towers for this phi
0107       // and this phi + 1
0108       ++m_ohNumStreak[iPhi];
0109       ++m_ohNumStreak[iUp];
0110
0111       // fill histograms
0112       m_hists["nstreakperphi"]->Fill(iPhi);
0113       m_hists["nstreakperphi"]->Fill(iPhi);
0114       m_hists["nstreaktwretavsphi"]->Fill(iEta, iPhi);
0115       m_hists["nstreaktwretavsphi"]->Fill(iEta, iUp);
0116
0117     }  // end eta loop
0118   }    // end phi loop
0119
0120   // now find longest streak
0121   const uint32_t nMaxStreak = *std::max_element(m_ohNumStreak.begin(), m_ohNumStreak.end());
0122   m_hists["nmaxstreak"]->Fill(nMaxStreak);
0123
0124   // return if streak length above threshold
0125   return (nMaxStreak > m_config.minNumTwrsInStreak);
0126
0127 }  // end 'ApplyFilter()'
0128
0129 // ----------------------------------------------------------------------------
0130 //! Construct histograms
0131 // ----------------------------------------------------------------------------
0132 void StreakSidebandFilter::BuildHistograms(const std::string& module, const std::string& tag)
0133 {
0134   // print debug message
0135   if (m_config.debug && (m_config.verbosity > 2))
0136   {
0137     std::cout << "StreakSidebandFilter::BuildHistograms(std::string) Constructing histograms" << std::endl;
0138   }
0139
0140   // make sure module name is lower case
0141   std::string moduleAndFilterName = module + "_" + m_name;
0142
0143   // names of variables to be histogramed
0144   const std::vector<std::string> varNames = {
0145       "nmaxstreak",
0146       "nstreakperphi",
0147       "nstreaktwretavsphi"};
0148
0149   // make qa-compliant hist names
0150   std::vector<std::string> histNames = BeamBackgroundFilterAndQADefs::MakeQAHistNames(varNames, moduleAndFilterName, tag);
0151
0152   // construct histograms
0153   //   - n.b. reminder that there are
0154   //       24 ohcal towers in eta
0155   //       64 ohcal towers in phi
0156   m_hists[varNames[0]] = new TH1D(histNames[0].data(), "", 25, -0.5, 24.5);
0157   m_hists[varNames[1]] = new TH1D(histNames[1].data(), "", 65, -0.5, 64.5);
0158   m_hists[varNames[2]] = new TH2D(histNames[2].data(), "", 25, -0.5, 24.5, 65, -0.5, 64.5);
0159   return;
0160
0161 }  // end 'BuildHistograms(std::string&, std::string&)'
0162
0163 // private methods ============================================================
0164
0165 // ----------------------------------------------------------------------------
0166 //! Grab input nodes
0167 // ----------------------------------------------------------------------------
0168 void StreakSidebandFilter::GrabNodes(PHCompositeNode* topNode)
0169 {
0170   // print debug message
0171   if (m_config.debug && (m_config.verbosity > 2))
0172   {
0173     std::cout << "StreakSidebandFilter::GrabNodes(PHCompositeNode*) Grabbing input nodes" << std::endl;
0174   }
0175
0176   m_ohContainer = findNode::getClass<TowerInfoContainer>(topNode, m_config.inNodeName);
0177   return;
0178
0179 }  // end 'GrabNodes(PHCompositeNode*)'
0180
0181 // ----------------------------------------------------------------------------
0182 //! Check if tower not consistent w/ being in a streak
0183 // ----------------------------------------------------------------------------
0184 bool StreakSidebandFilter::IsTowerNotStreaky(const BeamBackgroundFilterAndQADefs::Tower& tower)
0185 {
0186   // print debug message
0187   if (m_config.debug && (m_config.verbosity > 2))
0188   {
0189     std::cout << "StreakSidebandFilter::IsTowerNotStreaky() Checking if tower not consistent w/ streak" << std::endl;
0190   }
0191
0192   const bool isBadStatus = !(tower.isGood);
0193   const bool isBelowEneCut = (tower.energy < m_config.minStreakTwrEne);
0194   return (isBadStatus || isBelowEneCut);
0195
0196 }  // end 'IsTowerNotStreaky(Tower& tower)'
0197
0198 // ----------------------------------------------------------------------------
0199 //! Check if a neighboring tower consistent w/ a streak
0200 // ----------------------------------------------------------------------------
0201 bool StreakSidebandFilter::IsNeighborNotStreaky(const BeamBackgroundFilterAndQADefs::Tower& tower)
0202 {
0203   // print debug message
0204   if (m_config.debug && (m_config.verbosity > 2))
0205   {
0206     std::cout << "StreakSidebandFilter::IsNeighborNotStreaky() Checking if neighboring tower not consistent w/ streak" << std::endl;
0207   }
0208
0209   const bool isBadStatus = !(tower.isGood);
0210   const bool isAboveEneCut = (tower.energy > m_config.maxAdjacentTwrEne);
0211   return (isBadStatus || isAboveEneCut);
0212
0213 }  // end 'IsNeighborNotStreaky(Tower& tower)'
0214
0215 // end ========================================================================