From 2e7540c874cfd4be4f3bae6c0aa363a883386154 Mon Sep 17 00:00:00 2001 From: Richard Tyson <58739511+rtysonCLAS12@users.noreply.github.com> Date: Wed, 31 Jul 2024 11:04:32 -0400 Subject: [PATCH] Sector Finder: Support for charged & neutral particles, added action function and validator (#243) --- .../clas12/SectorFinder/Algorithm.cc | 168 ++++++++++++------ .../clas12/SectorFinder/Algorithm.h | 63 ++++++- .../clas12/SectorFinder/Config.yaml | 10 +- .../clas12/SectorFinder/Validator.cc | 142 +++++++++++++++ .../clas12/SectorFinder/Validator.h | 41 +++++ src/iguana/algorithms/meson.build | 2 +- 6 files changed, 365 insertions(+), 61 deletions(-) create mode 100644 src/iguana/algorithms/clas12/SectorFinder/Validator.cc create mode 100644 src/iguana/algorithms/clas12/SectorFinder/Validator.h diff --git a/src/iguana/algorithms/clas12/SectorFinder/Algorithm.cc b/src/iguana/algorithms/clas12/SectorFinder/Algorithm.cc index 8f2e5ff1..e3ea605c 100644 --- a/src/iguana/algorithms/clas12/SectorFinder/Algorithm.cc +++ b/src/iguana/algorithms/clas12/SectorFinder/Algorithm.cc @@ -9,80 +9,136 @@ namespace iguana::clas12 { // define options, their default values, and cache them ParseYAMLConfig(); - o_bankname = GetOptionScalar("bank"); + o_bankname_charged = GetOptionScalar("bank_charged"); + o_bankname_uncharged = GetOptionScalar("bank_uncharged"); + bool setDefaultBanks=false; // get expected bank indices b_particle = GetBankIndex(banks, "REC::Particle"); - if(o_bankname != "default") { - b_user = GetBankIndex(banks, o_bankname); - userSpecifiedBank = true; + if(o_bankname_charged != "default") { + b_user_charged = GetBankIndex(banks, o_bankname_charged); + userSpecifiedBank_charged = true; } else { b_calorimeter = GetBankIndex(banks, "REC::Calorimeter"); b_track = GetBankIndex(banks, "REC::Track"); b_scint = GetBankIndex(banks, "REC::Scintillator"); - userSpecifiedBank = false; + setDefaultBanks = true; + userSpecifiedBank_charged = false; + } + + if(o_bankname_uncharged != "default") { + b_user_uncharged = GetBankIndex(banks, o_bankname_uncharged); + userSpecifiedBank_uncharged = true; + } + else { + //avoid setting default banks twice + if(!setDefaultBanks){ + b_calorimeter = GetBankIndex(banks, "REC::Calorimeter"); + b_track = GetBankIndex(banks, "REC::Track"); + b_scint = GetBankIndex(banks, "REC::Scintillator"); + setDefaultBanks = true; + } + userSpecifiedBank_uncharged = false; } // create the output bank // FIXME: generalize the groupid and itemid - auto result_schema = CreateBank(banks, b_result, "REC::Particle::Sector", {"sector/I"}, 0xF000, 2); + auto result_schema = CreateBank(banks, b_result, "REC::Particle::Sector", {"sector/I","pindex/I"}, 0xF000, 4); i_sector = result_schema.getEntryOrder("sector"); + i_pindex = result_schema.getEntryOrder("pindex"); } - void SectorFinder::Run(hipo::banklist& banks) const { auto& particleBank = GetBank(banks, b_particle, "REC::Particle"); auto& resultBank = GetBank(banks, b_result, "REC::Particle::Sector"); + std::vector sectors_uncharged; + std::vector pindices_uncharged; + if(userSpecifiedBank_uncharged){ + auto const& userBank = GetBank(banks, b_user_uncharged); + GetListsSectorPindex(userBank,sectors_uncharged,pindices_uncharged); + } + + std::vector sectors_charged; + std::vector pindices_charged; + if(userSpecifiedBank_charged){ + auto const& userBank = GetBank(banks, b_user_charged); + GetListsSectorPindex(userBank,sectors_charged,pindices_charged); + } + + std::vector sectors_track; + std::vector pindices_track; + if(!userSpecifiedBank_charged || !userSpecifiedBank_uncharged){ + auto const& trackBank = GetBank(banks, b_track); + GetListsSectorPindex(trackBank,sectors_track,pindices_track); + } + + std::vector sectors_cal; + std::vector pindices_cal; + if(!userSpecifiedBank_charged || !userSpecifiedBank_uncharged){ + auto const& calBank = GetBank(banks, b_calorimeter); + GetListsSectorPindex(calBank,sectors_cal,pindices_cal); + } + + std::vector sectors_scint; + std::vector pindices_scint; + if(!userSpecifiedBank_charged || !userSpecifiedBank_uncharged){ + auto const& scintBank = GetBank(banks, b_scint); + GetListsSectorPindex(scintBank,sectors_scint,pindices_scint); + } + + // sync new bank with particle bank, and fill it with zeroes resultBank.setRows(particleBank.getRows()); resultBank.getMutableRowList().setList(particleBank.getRowList()); - for(int row = 0; row < resultBank.getRows(); row++) + for(int row = 0; row < resultBank.getRows(); row++){ resultBank.putInt(i_sector, row, 0); - - // filter the input bank for requested PDG code(s) - if(userSpecifiedBank) { // if user specified a specific bank - auto const& userBank = GetBank(banks, b_user); - for(auto const& row : particleBank.getRowList()) { - if(userBank.getRowList().size() > 0) { - resultBank.putInt(i_sector, row, GetSector(userBank, row)); - } - } + resultBank.putInt(i_pindex, row, row); } - else { // use the standard method - auto const& calBank = GetBank(banks, b_calorimeter); - auto const& scintBank = GetBank(banks, b_scint); - auto const& trackBank = GetBank(banks, b_track); - for(auto const& row : particleBank.getRowList()) { - int trackSector = 0; - if(trackBank.getRowList().size() > 0) { - trackSector = GetSector(trackBank, row); - } - int scintSector = 0; - if(scintBank.getRowList().size() > 0) { - scintSector = GetSector(scintBank, row); - } - int calSector = 0; - if(calBank.getRowList().size() > 0) { - calSector = GetSector(calBank, row); - } + for(int row = 0; row < particleBank.getRows(); row++) { + int charge=particleBank.getInt("charge",row); - if(trackSector != 0) { - resultBank.putInt(i_sector, row, trackSector); - } - else if(scintSector != 0) { - resultBank.putInt(i_sector, row, scintSector); + bool userSp = charge==0 ? userSpecifiedBank_uncharged : userSpecifiedBank_charged; + std::vector& sct_us = charge==0 ? sectors_uncharged : sectors_charged; + std::vector& pin_us = charge==0 ? pindices_uncharged : pindices_charged; + + + + if(userSp){ + int sect=GetSector(sct_us, pin_us,row); + if (sect!=-1){ + resultBank.putInt(i_sector, row, sect); + resultBank.putInt(i_pindex, row, row); } - else { - // FIXME: add even if calSector is 0 - // need an entry per pindex - // can happen that particle not in FD - // ie sector is 0 - resultBank.putInt(i_sector, row, calSector); + } else { + enum det_enum {kTrack, kScint, kCal, nDet}; // try to get sector from these detectors, in this order + for(int d = 0; d < nDet; d++) { + int sect = -1; + std::string det_name; + switch(d) { + case kTrack: + sect=GetSector(sectors_track, pindices_track,row); + det_name="track"; + break; + case kScint: + sect=GetSector(sectors_scint, pindices_scint,row); + det_name="scint"; + break; + case kCal: + sect=GetSector(sectors_cal, pindices_cal,row); + det_name="cal"; + break; + } + if(sect != -1) { + m_log->Trace("{} pindex {} sect {}", det_name, row, sect); + resultBank.putInt(i_sector, row, sect); + resultBank.putInt(i_pindex, row, row); + break; + } } } } @@ -90,15 +146,27 @@ namespace iguana::clas12 { ShowBank(resultBank, Logger::Header("CREATED BANK")); } - int SectorFinder::GetSector(hipo::bank const& bank, int const pindex) const + void SectorFinder::GetListsSectorPindex(hipo::bank const& bank, std::vector& sectors, std::vector& pindices) const { - int sector = 0; for(auto const& row : bank.getRowList()) { - if(bank.getInt("pindex", row) == pindex) { - return bank.getInt("sector", row); + //check that we're only using FD detectors + //eg have "sectors" in CND which we don't want to add here + int det=bank.getInt("detector",row); + if (listFDDets.find(det) != listFDDets.end()) { + sectors.push_back(bank.getInt("sector", row)); + pindices.push_back(bank.getInt("pindex", row)); + } + } + } + + int SectorFinder::GetSector(std::vector const& sectors, std::vector const& pindices, int const pindex) const + { + for(std::size_t i=0;i sectors; + /// std::vector pindices + /// + /// //bank is a hipo::bank object from which we want to get the sectors + /// //for example the bank object related to REC::Calorimeter + /// for(auto const& row : bank.getRowList()) { + /// + /// int det=bank.getInt("detector",row); + /// + /// //NB: you should check you read from an FD detector + /// // eg det 7 is the ECAL + /// if(det==7){ + /// sectors.push_back(bank.getInt("sector", row)); + /// pindices.push_back(bank.getInt("pindex", row)); + /// } + /// } + /// + /// //partbank is a hipo::bank object related to REC::Particle + /// //algo_sector_finder is the iguana::clas12::SectorFinder object + /// for(auto const& row : partbank.getRowList()) { + /// int sector = algo_sector_finder.GetSector(sectors, pindices, row); + /// } + /// ``` + /// + /// @param sectors list of sectors in a detector bank + /// @param pindices list of pindices in a detector bank + /// @param pindex index in `REC::Particle` bank for which to get sector + /// @returns sector for `pindex` in list, -1 if `pindex` not in inputted list + int GetSector(std::vector const& sectors, std::vector const& pindices, int const pindex) const; + + /// fill lists of sectors and pindices present in the input bank + /// @param bank bank from which to get lists of sectors and pindices + /// @param sectors list to fill with sectors in the bank + /// @param pindices list to fill with pindices in the bank + void GetListsSectorPindex(hipo::bank const& bank, std::vector& sectors, std::vector& pindices) const; private: @@ -38,15 +78,22 @@ namespace iguana::clas12 { hipo::banklist::size_type b_calorimeter; hipo::banklist::size_type b_track; hipo::banklist::size_type b_scint; - hipo::banklist::size_type b_user; + hipo::banklist::size_type b_user_charged; + hipo::banklist::size_type b_user_uncharged; hipo::banklist::size_type b_result; - bool userSpecifiedBank{false}; + bool userSpecifiedBank_charged{false}; + bool userSpecifiedBank_uncharged{true}; // `b_result` bank item indices int i_sector; + int i_pindex; /// Configuration options - std::string o_bankname; + std::string o_bankname_charged; + std::string o_bankname_uncharged; + + //only want sectors from FD detectors + std::set const listFDDets{6,7,12,15,16,18}; }; } diff --git a/src/iguana/algorithms/clas12/SectorFinder/Config.yaml b/src/iguana/algorithms/clas12/SectorFinder/Config.yaml index 8edf1b3d..fd1aeda3 100644 --- a/src/iguana/algorithms/clas12/SectorFinder/Config.yaml +++ b/src/iguana/algorithms/clas12/SectorFinder/Config.yaml @@ -1,5 +1,11 @@ clas12::SectorFinder: # use the default banks - bank: default + # for both charged/uncharged particles + bank_charged: default + bank_uncharged: default + # or force finder to use a certain bank with eg: - #bank: REC::Calorimeter + #can have default for charged but not uncharged particles + #or vice versa + #bank_charged: default + #bank_uncharged: REC::Calorimeter diff --git a/src/iguana/algorithms/clas12/SectorFinder/Validator.cc b/src/iguana/algorithms/clas12/SectorFinder/Validator.cc new file mode 100644 index 00000000..b024e4a6 --- /dev/null +++ b/src/iguana/algorithms/clas12/SectorFinder/Validator.cc @@ -0,0 +1,142 @@ +#include "Validator.h" + +#include +#include + +namespace iguana::clas12 { + + REGISTER_IGUANA_VALIDATOR(SectorFinderValidator); + + void SectorFinderValidator::Start(hipo::banklist& banks) + { + // define the algorithm sequence + m_algo_seq = std::make_unique(); + m_algo_seq->Add("clas12::EventBuilderFilter"); + m_algo_seq->Add("clas12::SectorFinder"); + m_algo_seq->SetOption>("clas12::EventBuilderFilter", "pids", u_pdg_list); + m_algo_seq->SetOption("clas12::SectorFinder", "bank_charged", "REC::Track"); + m_algo_seq->SetOption("clas12::SectorFinder", "bank_uncharged", "default"); + m_algo_seq->Start(banks); + + + b_particle = GetBankIndex(banks, "REC::Particle"); + b_cal = GetBankIndex(banks, "REC::Calorimeter"); + b_sector = GetBankIndex(banks, "REC::Particle::Sector"); + + // set an output file + auto output_dir = GetOutputDirectory(); + if(output_dir) { + m_output_file_basename = output_dir.value() + "/sector_finder"; + m_output_file = new TFile(m_output_file_basename + ".root", "RECREATE"); + } + + // define plots + gStyle->SetOptStat(0); + for(auto const& pdg : u_pdg_list) { + std::vector YvsX; + TString particle_name = particle::name.at(particle::PDG(pdg)); + TString particle_title = particle::title.at(particle::PDG(pdg)); + for(int sec = 1; sec <= 6; sec++) { + TString sector_name = Form("sec%d", sec); + TString sector_title = Form("sector %d", sec); + YvsX.push_back(new TH2D( + "YvsX_" + particle_name + "_" + sector_name, + particle_title + " Calorimeter Hit Position, " + sector_title + ";X [cm];Y [cm]", + 50, -500, 500, + 50, -500, 500)); + } + u_YvsX.insert({pdg, YvsX}); + } + u_IsInFD = new TH1D( + "IsInFD", + "e^{-} with #theta>6.5^{o} Sector; e^{-} Sector", + 7, -0.5, 6.5); + + } + + void SectorFinderValidator::Run(hipo::banklist& banks) const + { + + auto& particle_bank = GetBank(banks, b_particle, "REC::Particle"); + auto& sector_bank = GetBank(banks, b_sector, "REC::Particle::Sector"); + auto& cal_bank = GetBank(banks, b_cal, "REC::Calorimeter"); + + + m_algo_seq->Run(banks); + + // lock the mutex, so we can mutate plots + std::scoped_lock lock(m_mutex); + // fill the plots + for(auto const& row : particle_bank.getRowList()) { + + auto pdg = particle_bank.getInt("pid", row); + auto sector = sector_bank.getInt("sector", row); + + double x=0,y=0; + for(auto const& rowcal : cal_bank.getRowList()){ + if(cal_bank.getInt("pindex", rowcal)==row){ + x=cal_bank.getFloat("x", rowcal); + y=cal_bank.getFloat("y", rowcal); + } + } + + if(pdg==11){ + double Px=particle_bank.getFloat("px", row); + double Py=particle_bank.getFloat("py", row); + double Pz=particle_bank.getFloat("pz", row); + double P=sqrt(Px*Px+Py*Py+Pz*Pz); + double Theta=acos(Pz/P) * 180./M_PI ; + //electrons are in FT or FD + //sector should always be 1 if theta is larger than 5.5 degrees + if (Theta>6.5){ + u_IsInFD->Fill(sector); + if(sector==0){ + m_log->Trace("e' with theta={} and sector==0, this should not happen", Theta); + } + } + } + + // skip central particle, or unknown sector + if(sector == 0) + continue; + m_log->Trace("Filling SectorFinder Validator, pdg {} sector {} pindex {}", pdg, sector, row); + u_YvsX.at(pdg).at(sector - 1)->Fill(x, y); + } + + } + + + void SectorFinderValidator::Stop() + { + if(GetOutputDirectory()) { + for(auto const& [pdg, plots] : u_YvsX) { + int n_cols = 3; + int n_rows = 2; + TString canv_name = Form("canv%d", pdg); + auto canv = new TCanvas(canv_name, canv_name, n_cols * 800, n_rows * 600); + canv->Divide(n_cols, n_rows); + int pad_num = 0; + for(auto const& plot : plots) { + auto pad = canv->GetPad(++pad_num); + pad->cd(); + pad->SetGrid(1, 1); + pad->SetLeftMargin(0.12); + pad->SetRightMargin(0.12); + pad->SetBottomMargin(0.12); + plot->Draw("colz"); + } + canv->SaveAs(m_output_file_basename + "_" + std::to_string(pdg) + ".png"); + + } + + auto canv1D = new TCanvas("1D canvas","1D canvas",800, 600); + u_IsInFD->Draw(); + canv1D->SaveAs(m_output_file_basename+"_elIsInFD.png"); + + m_output_file->Write(); + m_log->Info("Wrote output file {}", m_output_file->GetName()); + m_output_file->Close(); + } + } + +} diff --git a/src/iguana/algorithms/clas12/SectorFinder/Validator.h b/src/iguana/algorithms/clas12/SectorFinder/Validator.h new file mode 100644 index 00000000..089c1afa --- /dev/null +++ b/src/iguana/algorithms/clas12/SectorFinder/Validator.h @@ -0,0 +1,41 @@ +#pragma once + +#include "iguana/algorithms/TypeDefs.h" +#include "iguana/algorithms/Validator.h" + +#include +#include +#include +#include + +namespace iguana::clas12 { + + /// @brief `iguana::clas12::SectorFinder` validator + class SectorFinderValidator : public Validator + { + + DEFINE_IGUANA_VALIDATOR(SectorFinderValidator, clas12::SectorFinderValidator) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + private: + + hipo::banklist::size_type b_particle; + hipo::banklist::size_type b_sector; + hipo::banklist::size_type b_cal; + + std::vector const u_pdg_list = { + particle::PDG::electron, + particle::PDG::photon,}; + + TString m_output_file_basename; + TFile* m_output_file; + mutable std::unordered_map> u_YvsX; + TH1D* u_IsInFD; + }; + +} diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 21fab07d..e2475324 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -61,7 +61,7 @@ algo_dict += [ { 'name': 'clas12::SectorFinder', 'has_c_bindings': false, - 'has_validator': false, + 'has_validator': true, 'test_args': { 'banks': [ 'REC::Particle', 'REC::Calorimeter', 'REC::Track', 'REC::Scintillator' ], },