Skip to content

Commit

Permalink
Implement particle filtering via NeulandPointFilter
Browse files Browse the repository at this point in the history
  • Loading branch information
Simon Lonzen authored and jose-luis-rs committed Dec 1, 2024
1 parent 2f6c548 commit e34265b
Show file tree
Hide file tree
Showing 18 changed files with 600 additions and 369 deletions.
1 change: 1 addition & 0 deletions CONTRIBUTORS
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ Labiche, Marc
Lagni, Andrea
Lihtar, Ivana
Loeher, Bastian
Lonzen, Simon
Milhomens da Fonseca, Leandro
Mozumdar, Nikhil
Panin, Valerii
Expand Down
7 changes: 4 additions & 3 deletions neuland/digitizing/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,8 @@ set(SRCS
R3BDigitizingPaddleNeuland.cxx
R3BNeulandHitMon.cxx
R3BNeulandDigitizer.cxx
R3BDigitizingPaddleMock.h
R3BDigitizingChannelMock.h)
R3BNeulandHitMon.cxx
NeulandPointFilter.cxx)

set(HEADERS
R3BDigitizingChannel.h
Expand All @@ -31,7 +31,8 @@ set(HEADERS
R3BDigitizingTacQuila.h
R3BDigitizingTamex.h
R3BNeulandDigitizer.h
R3BNeulandHitMon.h)
R3BNeulandHitMon.h
NeulandPointFilter.h)

add_library_with_dictionary(
LIBNAME
Expand Down
18 changes: 18 additions & 0 deletions neuland/digitizing/NeulandPointFilter.cxx
Original file line number Diff line number Diff line change
@@ -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_));
}
94 changes: 94 additions & 0 deletions neuland/digitizing/NeulandPointFilter.h
Original file line number Diff line number Diff line change
@@ -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
};
105 changes: 53 additions & 52 deletions neuland/digitizing/R3BNeulandDigitizer.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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()
Expand All @@ -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() << " "
Expand All @@ -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)
Expand All @@ -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,
Expand All @@ -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
Loading

0 comments on commit e34265b

Please sign in to comment.