Reapply "Implement particle filtering via NeulandPointFilter"

This reverts commit 7826525.

Author: YanzhaoW

.zenodo.json

@@ -324,6 +324,11 @@
             "orcid": "0000-0002-6506-8562",
             "affiliation": "STFC Daresbury Laboratory, Warrington, United Kingdom"
         },
+        {
+            "type": "Other",
+            "name": "Lonzen, Simon",
+            "affiliation": "Universität zu Köln, 50923 Köln, Germany"
+        },
         {
             "type": "Other",
             "name": "Lorenz, Enis",

CONTRIBUTORS

@@ -14,6 +14,7 @@ Jelavic-Malenica, Deša [RBI Zagreb, HR10000 Zagreb, Croatia]
 Jenegger, Tobias [Technische Universität München, 85748 Garching, Germany]
 Kripko, Aron
 Labiche, Marc [https://orcid.org/0000-0002-6506-8562] [STFC Daresbury Laboratory, Warrington, United Kingdom]
+Lonzen, Simon
 Lorenz, Enis [Technische Universität Darmstadt, Fachbereich Physik, Institut für Kernphysik, 64289 Darmstadt, Germany]
 Paschalis, Stefanos [https://orcid.org/0000-0002-9113-3778] [School of Physics, Engineering and Technology, University of York, YO10 5DD York, United Kingdom]
 Perea, Angel [Instituto de Estructura de la Materia, CSIC, 28006 Madrid, Spain]

neuland/digitizing/CMakeLists.txt

@@ -18,8 +18,8 @@ set(SRCS
     R3BDigitizingPaddleNeuland.cxx
     R3BNeulandHitMon.cxx
     R3BNeulandDigitizer.cxx
-    R3BDigitizingPaddleMock.h
-    R3BDigitizingChannelMock.h)
+    R3BNeulandHitMon.cxx
+    NeulandPointFilter.cxx)
 
 set(HEADERS
     R3BDigitizingChannel.h
@@ -31,7 +31,8 @@ set(HEADERS
     R3BDigitizingTacQuila.h
     R3BDigitizingTamex.h
     R3BNeulandDigitizer.h
-    R3BNeulandHitMon.h)
+    R3BNeulandHitMon.h
+    NeulandPointFilter.h)
 
 add_library_with_dictionary(
     LIBNAME

neuland/digitizing/NeulandPointFilter.cxx

@@ -0,0 +1,18 @@
+#include "NeulandPointFilter.h"
+
+void NeulandPointFilter::SetFilter(R3B::Neuland::BitSetParticle filtered_particles)
+{
+    filtered_particles_ = filtered_particles;
+}
+void NeulandPointFilter::SetFilter(R3B::Neuland::BitSetParticle filtered_particles, double minimum_allowed_energy)
+{
+    filtered_particles_ = filtered_particles;
+    minimum_allowed_energy_ = minimum_allowed_energy;
+}
+
+auto NeulandPointFilter::CheckFiltered(const R3BNeulandPoint& neuland_point) -> bool
+{
+    return (
+        R3B::Neuland::CheckCriteria(R3B::Neuland::PidToBitSetParticle(neuland_point.GetPID()), filtered_particles_) or
+        (neuland_point.GetEnergyLoss() < minimum_allowed_energy_));
+}

neuland/digitizing/NeulandPointFilter.h

@@ -0,0 +1,94 @@
+#pragma once
+#include "R3BNeulandPoint.h"
+#include <bitset>
+
+namespace R3B::Neuland
+{
+    enum class BitSetParticle : uint32_t
+    {
+        none = 0x0000,
+        proton = 0x0001,
+        neutron = 0x0002,
+        electron = 0x0004,
+        positron = 0x0008,
+        alpha = 0x0010,
+        gamma = 0x0020,
+        meson = 0x0040,
+        other = 0x40000000
+    };
+
+    const std::unordered_map<int, BitSetParticle> PidToBitSetParticleHash = {
+        { 2212, BitSetParticle::proton },  { 2112, BitSetParticle::neutron },     { 11, BitSetParticle::electron },
+        { -11, BitSetParticle::positron }, { 1000040020, BitSetParticle::alpha }, { 22, BitSetParticle::gamma },
+        { 0, BitSetParticle::none }
+    };
+
+    constexpr auto ParticleBitsetSize = 32U;
+
+    using ParticleUType = std::underlying_type_t<BitSetParticle>;
+
+    constexpr auto ParticleToBitSet(BitSetParticle particle)
+    {
+        return std::bitset<ParticleBitsetSize>{ static_cast<ParticleUType>(particle) };
+    }
+
+    inline auto BitSetToParticle(std::bitset<ParticleBitsetSize> bits) -> BitSetParticle
+    {
+        return static_cast<BitSetParticle>(static_cast<uint32_t>(bits.to_ulong()));
+    }
+
+    inline auto CheckCriteria(BitSetParticle particle, BitSetParticle criteria) -> bool
+    {
+        return (ParticleToBitSet(particle) & ParticleToBitSet(criteria)) == ParticleToBitSet(particle);
+    }
+
+    inline auto operator|(BitSetParticle left, BitSetParticle right) -> BitSetParticle
+    {
+        auto left_bitset = ParticleToBitSet(left);
+        auto right_bitset = ParticleToBitSet(right);
+        return BitSetToParticle(left_bitset | right_bitset);
+    }
+
+    inline auto operator&(BitSetParticle left, BitSetParticle right) -> BitSetParticle
+    {
+        auto left_bitset = ParticleToBitSet(left);
+        auto right_bitset = ParticleToBitSet(right);
+        return BitSetToParticle(left_bitset & right_bitset);
+    }
+
+    inline auto operator~(BitSetParticle particle) -> BitSetParticle
+    {
+        return BitSetToParticle(~ParticleToBitSet(particle));
+    }
+
+    inline auto PidToBitSetParticle(int pid) -> BitSetParticle
+    {
+        if (pid > 99 and pid < 1000) // NOLINT mesons have three digit pdgs
+        {
+            return BitSetParticle::meson;
+        }
+        auto pid_to_bitset_hash_iterator = PidToBitSetParticleHash.find(pid);
+
+        if (pid_to_bitset_hash_iterator == PidToBitSetParticleHash.end())
+        {
+            return BitSetParticle::other;
+        }
+
+        return pid_to_bitset_hash_iterator->second;
+    }
+
+} // namespace R3B::Neuland
+class NeulandPointFilter
+{
+  public:
+    NeulandPointFilter() = default;
+    void SetFilter(R3B::Neuland::BitSetParticle filtered_particles);
+    void SetFilter(R3B::Neuland::BitSetParticle filtered_particles, double minimum_allowed_energy);
+    [[nodiscard]] auto GetFilter() const -> R3B::Neuland::BitSetParticle { return filtered_particles_; }
+    [[nodiscard]] auto GetMinimumAllowedEnergy() const -> double { return minimum_allowed_energy_; }
+    auto CheckFiltered(const R3BNeulandPoint& neuland_point) -> bool;
+
+  private:
+    R3B::Neuland::BitSetParticle filtered_particles_ = R3B::Neuland::BitSetParticle::none;
+    double minimum_allowed_energy_ = 0; // engergy in GeV
+};

neuland/digitizing/R3BNeulandDigitizer.cxx

@@ -16,6 +16,8 @@
 #include "FairRootManager.h"
 #include "FairRunAna.h"
 #include "FairRuntimeDb.h"
+#include "NeulandPointFilter.h"
+#include "R3BDataMonitor.h"
 #include "TGeoManager.h"
 #include "TGeoNode.h"
 #include "TH1F.h"
@@ -26,28 +28,25 @@
 #include <TFile.h>
 #include <iostream>
 #include <stdexcept>
+#include <string_view>
 #include <utility>
 
-R3BNeulandDigitizer::R3BNeulandDigitizer(TString input, TString output)
-    : R3BNeulandDigitizer(Digitizing::CreateEngine(UsePaddle<NeulandPaddle>(), UseChannel<TacquilaChannel>()),
-                          std::move(input),
-                          std::move(output))
+R3BNeulandDigitizer::R3BNeulandDigitizer()
+    : R3BNeulandDigitizer(Digitizing::CreateEngine(UsePaddle<NeulandPaddle>(), UseChannel<TacquilaChannel>())
+
+      )
 {
 }
 
-R3BNeulandDigitizer::R3BNeulandDigitizer(std::unique_ptr<Digitizing::DigitizingEngineInterface> engine,
-                                         TString input,
-                                         TString output)
+R3BNeulandDigitizer::R3BNeulandDigitizer(std::unique_ptr<Digitizing::DigitizingEngineInterface> engine)
     : FairTask("R3BNeulandDigitizer")
-    , fPoints(std::move(input))
-    , fHits(std::move(output))
-    , fDigitizingEngine(std::move(engine))
+    , digitizing_engine_(std::move(engine))
 {
 }
 
 void R3BNeulandDigitizer::SetEngine(std::unique_ptr<Digitizing::DigitizingEngineInterface> engine)
 {
-    fDigitizingEngine = std::move(engine);
+    digitizing_engine_ = std::move(engine);
 }
 
 void R3BNeulandDigitizer::SetParContainers()
@@ -64,48 +63,64 @@ void R3BNeulandDigitizer::SetParContainers()
         LOG(fatal) << "R3BNeulandDigitizer::SetParContainers: No runtime database";
     }
 
-    fNeulandGeoPar = dynamic_cast<R3BNeulandGeoPar*>(rtdb->getContainer("R3BNeulandGeoPar"));
-    if (fNeulandGeoPar == nullptr)
+    neuland_geo_par_ = dynamic_cast<R3BNeulandGeoPar*>(rtdb->getContainer("R3BNeulandGeoPar"));
+    if (neuland_geo_par_ == nullptr)
     {
         LOG(fatal) << "R3BNeulandDigitizer::SetParContainers: No R3BNeulandGeoPar";
     }
 
-    fDigitizingEngine->Init();
+    digitizing_engine_->Init();
 }
 
-InitStatus R3BNeulandDigitizer::Init()
+void R3BNeulandDigitizer::SetNeulandPointFilter(R3B::Neuland::BitSetParticle particle)
+{
+    neuland_point_filter_.SetFilter(particle);
+}
+void R3BNeulandDigitizer::SetNeulandPointFilter(R3B::Neuland::BitSetParticle particle,
+                                                double minimum_allowed_energy_gev)
 {
-    fPoints.Init();
-    fHits.Init();
+    neuland_point_filter_.SetFilter(particle, minimum_allowed_energy_gev);
+}
 
+auto R3BNeulandDigitizer::Init() -> InitStatus
+{
+    neuland_points_.init();
+    neuland_hits_.init();
     // Initialize control histograms
     auto const PaddleMulSize = 3000;
-    hMultOne = R3B::root_owned<TH1I>(
+    mult_one_ = data_monitor_.add_hist<TH1I>(
         "MultiplicityOne", "Paddle multiplicity: only one PMT per paddle", PaddleMulSize, 0, PaddleMulSize);
-    hMultTwo = R3B::root_owned<TH1I>(
+
+    mult_two_ = data_monitor_.add_hist<TH1I>(
         "MultiplicityTwo", "Paddle multiplicity: both PMTs of a paddle", PaddleMulSize, 0, PaddleMulSize);
     auto const timeBinSize = 200;
-    hRLTimeToTrig = R3B::root_owned<TH1F>("hRLTimeToTrig", "R/Ltime-triggerTime", timeBinSize, -100., 100.);
+    rl_time_to_trig_ = data_monitor_.add_hist<TH1F>("hRLTimeToTrig", "R/Ltime-triggerTime", timeBinSize, -100., 100.);
 
     return kSUCCESS;
 }
 
 void R3BNeulandDigitizer::Exec(Option_t* /*option*/)
 {
-    fHits.Reset();
+    neuland_hits_.clear();
     const auto GeVToMeVFac = 1000.;
 
     std::map<UInt_t, Double_t> paddleEnergyDeposit;
     // Look at each Land Point, if it deposited energy in the scintillator, store it with reference to the bar
-    for (const auto& point : fPoints.Retrieve())
+    for (const auto& point : neuland_points_)
     {
-        if (point->GetEnergyLoss() > 0.)
+        if (((neuland_point_filter_.GetFilter() != R3B::Neuland::BitSetParticle::none) or
+             (neuland_point_filter_.GetMinimumAllowedEnergy() != 0)) and
+            neuland_point_filter_.CheckFiltered(point))
         {
-            const Int_t paddleID = point->GetPaddle();
+            continue;
+        }
+        if (point.GetEnergyLoss() > 0.)
+        {
+            const Int_t paddleID = point.GetPaddle();
 
             // Convert position of point to paddle-coordinates, including any rotation or translation
-            const TVector3 position = point->GetPosition();
-            const TVector3 converted_position = fNeulandGeoPar->ConvertToLocalCoordinates(position, paddleID);
+            const TVector3 position = point.GetPosition();
+            const TVector3 converted_position = neuland_geo_par_->ConvertToLocalCoordinates(position, paddleID);
             LOG(debug2) << "NeulandDigitizer: Point in paddle " << paddleID
                         << " with global position XYZ: " << position.X() << " " << position.Y() << " " << position.Z();
             LOG(debug2) << "NeulandDigitizer: Converted to local position XYZ: " << converted_position.X() << " "
@@ -114,22 +129,22 @@ void R3BNeulandDigitizer::Exec(Option_t* /*option*/)
             // Within the paddle frame, the relevant distance of the light from the pmt is always given by the
             // X-Coordinate
             const Double_t dist = converted_position.X();
-            fDigitizingEngine->DepositLight(paddleID, point->GetTime(), point->GetLightYield() * GeVToMeVFac, dist);
-            paddleEnergyDeposit[paddleID] += point->GetEnergyLoss() * GeVToMeVFac;
+            digitizing_engine_->DepositLight(paddleID, point.GetTime(), point.GetLightYield() * GeVToMeVFac, dist);
+            paddleEnergyDeposit[paddleID] += point.GetEnergyLoss() * GeVToMeVFac;
         } // eloss
     }     // points
 
-    const Double_t triggerTime = fDigitizingEngine->GetTriggerTime();
-    const auto paddles = fDigitizingEngine->ExtractPaddles();
+    const Double_t triggerTime = digitizing_engine_->GetTriggerTime();
+    const auto paddles = digitizing_engine_->ExtractPaddles();
 
     // Fill control histograms
-    hMultOne->Fill(static_cast<int>(std::count_if(
+    mult_one_->Fill(static_cast<int>(std::count_if(
         paddles.begin(), paddles.end(), [](const auto& keyValue) { return keyValue.second->HasHalfFired(); })));
 
-    hMultTwo->Fill(static_cast<int>(std::count_if(
+    mult_two_->Fill(static_cast<int>(std::count_if(
         paddles.begin(), paddles.end(), [](const auto& keyValue) { return keyValue.second->HasFired(); })));
 
-    hRLTimeToTrig->Fill(triggerTime);
+    rl_time_to_trig_->Fill(triggerTime);
 
     // Create Hits
     for (const auto& [paddleID, paddle] : paddles)
@@ -144,8 +159,8 @@ void R3BNeulandDigitizer::Exec(Option_t* /*option*/)
         for (const auto signal : signals)
         {
             const TVector3 hitPositionLocal = TVector3(signal.position, 0., 0.);
-            const TVector3 hitPositionGlobal = fNeulandGeoPar->ConvertToGlobalCoordinates(hitPositionLocal, paddleID);
-            const TVector3 hitPixel = fNeulandGeoPar->ConvertGlobalToPixel(hitPositionGlobal);
+            const TVector3 hitPositionGlobal = neuland_geo_par_->ConvertToGlobalCoordinates(hitPositionLocal, paddleID);
+            const TVector3 hitPixel = neuland_geo_par_->ConvertGlobalToPixel(hitPositionGlobal);
 
             R3BNeulandHit hit(paddleID,
                               signal.leftChannel.tdc,
@@ -157,30 +172,16 @@ void R3BNeulandDigitizer::Exec(Option_t* /*option*/)
                               hitPositionGlobal,
                               hitPixel);
 
-            if (fHitFilters.IsValid(hit))
+            if (neuland_hit_filters_.IsValid(hit))
             {
-                fHits.Insert(std::move(hit));
+                neuland_hits_.get().emplace_back(std::move(hit));
                 LOG(debug) << "Adding neuland hit with id = " << paddleID << ", time = " << signal.time
                            << ", energy = " << signal.energy;
             }
         } // loop over all hits for each paddle
     }     // loop over paddles
 
-    LOG(debug) << "R3BNeulandDigitizer: produced " << fHits.Size() << " hits";
-}
-
-void R3BNeulandDigitizer::Finish()
-{
-    TDirectory* tmp = gDirectory;
-    FairRootManager::Instance()->GetOutFile()->cd();
-
-    gDirectory->mkdir("R3BNeulandDigitizer");
-    gDirectory->cd("R3BNeulandDigitizer");
-
-    hMultOne->Write();
-    hMultTwo->Write();
-
-    gDirectory = tmp;
+    LOG(debug) << "R3BNeulandDigitizer: produced " << neuland_hits_.get().size() << " hits";
 }
 
 ClassImp(R3BNeulandDigitizer); // NOLINT