BeamBackgroundFilterAndQA/src/BeamBackgroundFilterAndQA.cc

0001 /// ===========================================================================
0002 /*! \file    BeamBackgroundFilterAndQA.h
0003  *  \authors Hanpu Jiang, Derek Anderson
0004  *  \date    10.21.2024
0005  *
0006  *  A F4A module to filter out events with significant
0007  *  beam background (the so-called "streaky events")
0008  *  and produce some relevant QA histograms.
0009  */
0010 /// ===========================================================================
0011
0012 #define BEAMBACKGROUNDFILTERANDQA_CC
0013
0014 // module components
0015 #include "BeamBackgroundFilterAndQA.h"
0016 #include "BeamBackgroundFilterAndQADefs.h"
0017
0018 // calo base
0019 #include <calobase/TowerInfoContainer.h>
0020
0021 // f4a libraries
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <fun4all/Fun4AllHistoManager.h>
0024 #include <ffaobjects/FlagSavev1.h>
0025
0026 // phool libraries
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h>
0029 #include <phool/PHCompositeNode.h>
0030
0031 // qa utilities
0032 #include <qautils/QAHistManagerDef.h>
0033
0034 // root libraries
0035 #include <TH1.h>
0036
0037 // c++ utiilites
0038 #include <cassert>
0039 #include <iostream>
0040
0041 // alias for convenience
0042 namespace bbfqd = BeamBackgroundFilterAndQADefs;
0043
0044
0045
0046 // ctor/dtor ==================================================================
0047
0048 // ----------------------------------------------------------------------------
0049 //! Default module constructor
0050 // ----------------------------------------------------------------------------
0051 BeamBackgroundFilterAndQA::BeamBackgroundFilterAndQA(const std::string& name, const bool debug)
0052   : SubsysReco(name)
0053   , m_manager(nullptr)
0054   , m_flags(name)
0055 {
0056
0057   // print debug message
0058   if (debug && (Verbosity() > 0))
0059   {
0060     std::cout << "BeamBackgroundFilterAndQA::BeamBackgroundFilterAndQA(const std::string &name) Calling ctor" << std::endl;
0061   }
0062
0063 }  // end ctor(std::string&, bool)'
0064
0065
0066
0067 // ----------------------------------------------------------------------------
0068 //! Module constructor accepting a configuration
0069 // ----------------------------------------------------------------------------
0070 BeamBackgroundFilterAndQA::BeamBackgroundFilterAndQA(const Config& config)
0071   : SubsysReco(config.moduleName)
0072   , m_manager(nullptr)
0073   , m_flags(config.moduleName)
0074   , m_config(config)
0075 {
0076
0077   // print debug message
0078   if (m_config.debug && (Verbosity() > 0))
0079   {
0080     std::cout << "BeamBackgroundFilterAndQA::BeamBackgroundFilterAndQA(const std::string &name) Calling ctor" << std::endl;
0081   }
0082
0083 }  // end ctor(BeamBackgroundFilterAndQAConfig&)'
0084
0085
0086
0087 // ----------------------------------------------------------------------------
0088 //! Module destructor
0089 // ----------------------------------------------------------------------------
0090 BeamBackgroundFilterAndQA::~BeamBackgroundFilterAndQA()
0091 {
0092
0093   // print debug message
0094   if (m_config.debug && (Verbosity() > 0))
0095   {
0096     std::cout << "BeamBackgroundFilterAndQA::~BeamBackgroundFilterAndQA() Calling dtor" << std::endl;
0097   }
0098
0099   /* nothing to do */
0100
0101 }  // end dtor()
0102
0103
0104
0105 // fun4all methods ============================================================
0106
0107 // ----------------------------------------------------------------------------
0108 //! Initialize module
0109 // ----------------------------------------------------------------------------
0110 int BeamBackgroundFilterAndQA::Init(PHCompositeNode* topNode)
0111 {
0112
0113   if (m_config.debug)
0114   {
0115     std::cout << "BeamBackgroundFilterAndQA::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0116   }
0117
0118   // initialize relevant filters, histograms
0119   InitFilters();
0120   InitFlags(topNode);
0121   BuildHistograms();
0122
0123   // if needed, initialize histograms + manager
0124   if (m_config.doQA)
0125   {
0126     InitHistManager();
0127     RegisterHistograms();
0128   }
0129   return Fun4AllReturnCodes::EVENT_OK;
0130
0131 }  // end 'Init(PHCompositeNode*)'
0132
0133
0134
0135 // ----------------------------------------------------------------------------
0136 //! Grab inputs, check for beam background, and fill histograms
0137 // ----------------------------------------------------------------------------
0138 int BeamBackgroundFilterAndQA::process_event(PHCompositeNode* topNode)
0139 {
0140
0141   if (m_config.debug)
0142   {
0143     std::cout << "BeamBackgroundFilterAndQA::process_event(PHCompositeNode *topNode) Processing event" << std::endl;
0144   }
0145
0146   // reset flags
0147   SetDefaultFlags();
0148
0149   // check for beam background
0150   const bool hasBeamBkgd = ApplyFilters(topNode);
0151
0152   // update flags on node tree
0153   UpdateFlags(topNode);
0154
0155   // if debugging, print out flags
0156   if (m_config.debug)
0157   {
0158     m_flags.printint();
0159   }
0160
0161   // if it does, abort event
0162   if (hasBeamBkgd && m_config.doEvtAbort)
0163   {
0164     return Fun4AllReturnCodes::ABORTEVENT;
0165   }
0166   else
0167   {
0168     return Fun4AllReturnCodes::EVENT_OK;
0169   }
0170
0171 }  // end 'process_event(PHCompositeNode*)'
0172
0173
0174
0175 // ----------------------------------------------------------------------------
0176 //! Run final calculations
0177 // ----------------------------------------------------------------------------
0178 int BeamBackgroundFilterAndQA::End(PHCompositeNode* /*topNode*/)
0179 {
0180
0181   if (m_config.debug)
0182   {
0183     std::cout << "BeamBackgroundFilterAndQA::End(PHCompositeNode *topNode) This is the end..." << std::endl;
0184   }
0185
0186   //... nothing to do ...//
0187   return Fun4AllReturnCodes::EVENT_OK;
0188
0189 }  // end 'End(PHCompositeNode*)'
0190
0191
0192
0193 // private methods ============================================================
0194
0195 // ----------------------------------------------------------------------------
0196 //! Initialize filters
0197 // ---------------------------------------------------------------------------
0198 void BeamBackgroundFilterAndQA::InitFilters()
0199 {
0200
0201   // print debug message
0202   if (m_config.debug && (Verbosity() > 1))
0203   {
0204     std::cout << "BeamBackgroundFilterAndQA::InitFilters() Initializing background filters" << std::endl;
0205   }
0206
0207   m_filters["Null"] = std::make_unique<NullFilter>( m_config.null, "Null" );
0208   m_filters["StreakSideband"] = std::make_unique<StreakSidebandFilter>( m_config.sideband, "StreakSideband" );
0209   //... other filters added here ...//
0210   return;
0211
0212 }  // end 'InitFilters()'
0213
0214
0215
0216 // ----------------------------------------------------------------------------
0217 //! Initialize flags
0218 // ----------------------------------------------------------------------------
0219 void BeamBackgroundFilterAndQA::InitFlags(PHCompositeNode* topNode)
0220 {
0221
0222   // print debug message
0223   if (m_config.debug && (Verbosity() > 1))
0224   {
0225     std::cout << "BeamBackgroundFilterAndQA::InitFlags() Initializing flags" << std::endl;
0226   }
0227
0228   // add node for flags
0229   PHNodeIterator itNode(topNode);
0230   PHCompositeNode* parNode = dynamic_cast<PHCompositeNode*>(itNode.findFirst("PHCompositeNode", "PAR"));
0231   if (!parNode)
0232   {
0233     std::cerr << PHWHERE << " WARNING: No PAR node found! Cannot add node for background flags to node tree!" << std::endl;
0234   }
0235   else
0236   {
0237     m_flags.SaveToNodeTree(parNode, m_config.flagPrefix);
0238   }
0239
0240   // initialize flags
0241   SetDefaultFlags();
0242   return;
0243
0244 }  // end 'InitFlags()'
0245
0246
0247
0248 // ----------------------------------------------------------------------------
0249 //! Initialize histogram manager
0250 // ----------------------------------------------------------------------------
0251 void BeamBackgroundFilterAndQA::InitHistManager()
0252 {
0253
0254   // print debug message
0255   if (m_config.debug && (Verbosity() > 0))
0256   {
0257     std::cout << "BeamBackgroundFilterAndQA::InitHistManager() Initializing histogram manager" << std::endl;
0258   }
0259
0260   m_manager = QAHistManagerDef::getHistoManager();
0261   if (!m_manager)
0262   {
0263     std::cerr << PHWHERE << ": PANIC! Couldn't grab histogram manager!" << std::endl;
0264     assert(m_manager);
0265   }
0266   return;
0267
0268 }  // end 'InitHistManager()'
0269
0270
0271
0272 // ----------------------------------------------------------------------------
0273 //! Build histograms
0274 // ----------------------------------------------------------------------------
0275 void BeamBackgroundFilterAndQA::BuildHistograms()
0276 {
0277
0278   // print debug message
0279   if (m_config.debug && (Verbosity() > 0))
0280   {
0281     std::cout << "BeamBackgroundFilterAndQA::BuildHistograms() Creating histograms" << std::endl;
0282   }
0283
0284   // construct module-wide variable names
0285   std::vector<std::string> varNames = {"nevts_overall"};
0286   for (const std::string& filterToApply : m_config.filtersToApply)
0287   {
0288     varNames.push_back("nevts_" + filterToApply);
0289   }
0290
0291   // get module-wide histogram names
0292   std::vector<std::string> histNames = bbfqd::MakeQAHistNames(varNames, m_config.moduleName, m_config.histTag);
0293
0294   // create module-wide histograms
0295   for (std::size_t iVar = 0; iVar < varNames.size(); ++iVar)
0296   {
0297     m_hists[varNames[iVar]] = new TH1D(histNames[iVar].data(), "", 3, -0.5, 2.5);
0298     m_hists[varNames[iVar]]->GetXaxis()->SetBinLabel(1, "All");
0299     m_hists[varNames[iVar]]->GetXaxis()->SetBinLabel(2, "No beam bkgd.");
0300     m_hists[varNames[iVar]]->GetXaxis()->SetBinLabel(3, "Beam bkgd.");
0301   }
0302
0303   // build filter-specific histograms
0304   for (const std::string& filterToApply : m_config.filtersToApply)
0305   {
0306     m_filters.at(filterToApply)->BuildHistograms(m_config.moduleName, m_config.histTag);
0307   }
0308   return;
0309
0310 }  // end 'BuildHistograms()'
0311
0312
0313
0314 // ----------------------------------------------------------------------------
0315 //! Register histograms
0316 // ----------------------------------------------------------------------------
0317 void BeamBackgroundFilterAndQA::RegisterHistograms()
0318 {
0319
0320   // print debug message
0321   if (m_config.debug && (Verbosity() > 0))
0322   {
0323     std::cout << "BeamBackgroundFilterAndQA::RegisterHistograms() Registering histograms w/ manager" << std::endl;
0324   }
0325
0326   // register module-wide histograms
0327   for (auto& hist : m_hists)
0328   {
0329     m_manager->registerHisto( hist.second );
0330   }
0331
0332   // register filter-specific histograms
0333   for (const std::string& filterToApply : m_config.filtersToApply)
0334   {
0335     m_filters.at(filterToApply)->RegisterHistograms(m_manager);
0336   }
0337   return;
0338
0339 }  // end 'RegisterHistograms()'
0340
0341
0342
0343 // ----------------------------------------------------------------------------
0344 //! Set default values of flags
0345 // ----------------------------------------------------------------------------
0346 void BeamBackgroundFilterAndQA::SetDefaultFlags()
0347 {
0348
0349   // print debug message
0350   if (m_config.debug && (Verbosity() > 1))
0351   {
0352     std::cout << "BeamBackgroundFilterAndQA::SetDefaultFlags() Setting defuault flag values" << std::endl;
0353   }
0354
0355   // set default values
0356   for (const std::string& filterToApply : m_config.filtersToApply)
0357   {
0358     m_flags.set_int_param(MakeFlagName(filterToApply), 0);
0359   }
0360   m_flags.set_int_param(MakeFlagName(), 0);
0361   return;
0362
0363 }  // end 'SetDefaultFlags()'
0364
0365
0366
0367 // ----------------------------------------------------------------------------
0368 //! Update flags on the node tree
0369 // ----------------------------------------------------------------------------
0370 void BeamBackgroundFilterAndQA::UpdateFlags(PHCompositeNode* topNode)
0371 {
0372
0373   // print debug message
0374   if (m_config.debug && (Verbosity() > 0))
0375   {
0376     std::cout << "BeamBackgroundFilterAndQA::UpdateFlags(PHCompositeNode*) Updating flags on the node tree" << std::endl;
0377   }
0378
0379   PHNodeIterator itNode(topNode);
0380   PHCompositeNode* parNode = dynamic_cast<PHCompositeNode*>(itNode.findFirst("PHCompositeNode", "PAR"));
0381   if (parNode)
0382   {
0383     m_flags.UpdateNodeTree(parNode, m_config.flagPrefix);
0384   }
0385   return;
0386
0387 }  // end 'UpdateFlags(PHCompositeNode*)'
0388
0389
0390
0391 // ----------------------------------------------------------------------------
0392 //! Apply relevant filters
0393 // ----------------------------------------------------------------------------
0394 bool BeamBackgroundFilterAndQA::ApplyFilters(PHCompositeNode* topNode)
0395 {
0396
0397   // print debug message
0398   if (m_config.debug && (Verbosity() > 0))
0399   {
0400     std::cout << "BeamBackgroundFilterAndQA::ApplyFilters(PHCompositeNode*) Creating histograms" << std::endl;
0401   }
0402
0403   // apply individual filters
0404   bool hasBkgd = false;
0405   for (const std::string& filterToApply : m_config.filtersToApply)
0406   {
0407     const bool filterFoundBkgd = m_filters.at(filterToApply)->ApplyFilter(topNode);
0408     if (filterFoundBkgd)
0409     {
0410       m_hists["nevts_" + filterToApply]->Fill(bbfqd::Status::HasBkgd);
0411       m_flags.set_int_param(MakeFlagName(filterToApply), 1);
0412     }
0413     else
0414     {
0415       m_hists["nevts_" + filterToApply]->Fill(bbfqd::Status::NoBkgd);
0416     }
0417     m_hists["nevts_" + filterToApply]->Fill(bbfqd::Status::Evt);
0418     hasBkgd += filterFoundBkgd;
0419   }
0420
0421   // fill overall histograms and return
0422   m_hists["nevts_overall"]->Fill(bbfqd::Status::Evt);
0423   if (hasBkgd)
0424   {
0425     m_hists["nevts_overall"]->Fill(bbfqd::Status::HasBkgd);
0426     m_flags.set_int_param(MakeFlagName(), 1);
0427   }
0428   else
0429   {
0430     m_hists["nevts_overall"]->Fill(bbfqd::Status::NoBkgd);
0431   }
0432   return hasBkgd;
0433
0434 }  // end 'ApplyFilters(PHCompositeNode*)'
0435
0436
0437
0438 // ----------------------------------------------------------------------------
0439 //! Create flag name
0440 // ----------------------------------------------------------------------------
0441 /*! Helper method to create a flag (PHParameters) name from a
0442  *  provided filter name.  If no filter is provided, method
0443  *  will return just the flag prefix.
0444  */
0445 std::string BeamBackgroundFilterAndQA::MakeFlagName(const std::string& filter)
0446 {
0447
0448   // print debug message
0449   if (m_config.debug && (Verbosity() > 2))
0450   {
0451     std::cout << "BeamBackgroundFilterAndQA::MakeFlagName(std::string&) Creating flag name" << std::endl;
0452   }
0453
0454   // by default, return flag prefix; otherwise,
0455   // combine prefix and filter name
0456   if (filter.empty())
0457   {
0458     return m_config.flagPrefix;
0459   }
0460   else
0461   {
0462     return m_config.flagPrefix + "_" + filter + "Filter";
0463   }
0464
0465 }  // end 'MakeFlagName(std::string&)'
0466
0467 // end ========================================================================