diff --git a/.zenodo.json b/.zenodo.json index ef1839f74..df964da27 100644 --- a/.zenodo.json +++ b/.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", diff --git a/CONTRIBUTORS b/CONTRIBUTORS index a08a03f25..bafd287b4 100644 --- a/CONTRIBUTORS +++ b/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] diff --git a/neuland/digitizing/CMakeLists.txt b/neuland/digitizing/CMakeLists.txt index 03e2207e0..f87504c5e 100644 --- a/neuland/digitizing/CMakeLists.txt +++ b/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 diff --git a/neuland/digitizing/NeulandPointFilter.cxx b/neuland/digitizing/NeulandPointFilter.cxx new file mode 100644 index 000000000..cd6a5b563 --- /dev/null +++ b/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_)); +} diff --git a/neuland/digitizing/NeulandPointFilter.h b/neuland/digitizing/NeulandPointFilter.h new file mode 100644 index 000000000..da44df046 --- /dev/null +++ b/neuland/digitizing/NeulandPointFilter.h @@ -0,0 +1,94 @@ +#pragma once +#include "R3BNeulandPoint.h" +#include + +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 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; + + constexpr auto ParticleToBitSet(BitSetParticle particle) + { + return std::bitset{ static_cast(particle) }; + } + + inline auto BitSetToParticle(std::bitset bits) -> BitSetParticle + { + return static_cast(static_cast(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 +}; diff --git a/neuland/digitizing/R3BNeulandDigitizer.cxx b/neuland/digitizing/R3BNeulandDigitizer.cxx index 0050be068..906f20c99 100644 --- a/neuland/digitizing/R3BNeulandDigitizer.cxx +++ b/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 #include #include +#include #include -R3BNeulandDigitizer::R3BNeulandDigitizer(TString input, TString output) - : R3BNeulandDigitizer(Digitizing::CreateEngine(UsePaddle(), UseChannel()), - std::move(input), - std::move(output)) +R3BNeulandDigitizer::R3BNeulandDigitizer() + : R3BNeulandDigitizer(Digitizing::CreateEngine(UsePaddle(), UseChannel()) + + ) { } -R3BNeulandDigitizer::R3BNeulandDigitizer(std::unique_ptr engine, - TString input, - TString output) +R3BNeulandDigitizer::R3BNeulandDigitizer(std::unique_ptr 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 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(rtdb->getContainer("R3BNeulandGeoPar")); - if (fNeulandGeoPar == nullptr) + neuland_geo_par_ = dynamic_cast(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( + mult_one_ = data_monitor_.add_hist( "MultiplicityOne", "Paddle multiplicity: only one PMT per paddle", PaddleMulSize, 0, PaddleMulSize); - hMultTwo = R3B::root_owned( + + mult_two_ = data_monitor_.add_hist( "MultiplicityTwo", "Paddle multiplicity: both PMTs of a paddle", PaddleMulSize, 0, PaddleMulSize); auto const timeBinSize = 200; - hRLTimeToTrig = R3B::root_owned("hRLTimeToTrig", "R/Ltime-triggerTime", timeBinSize, -100., 100.); + rl_time_to_trig_ = data_monitor_.add_hist("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 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(std::count_if( + mult_one_->Fill(static_cast(std::count_if( paddles.begin(), paddles.end(), [](const auto& keyValue) { return keyValue.second->HasHalfFired(); }))); - hMultTwo->Fill(static_cast(std::count_if( + mult_two_->Fill(static_cast(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 diff --git a/neuland/digitizing/R3BNeulandDigitizer.h b/neuland/digitizing/R3BNeulandDigitizer.h index d9224944c..274bd6895 100644 --- a/neuland/digitizing/R3BNeulandDigitizer.h +++ b/neuland/digitizing/R3BNeulandDigitizer.h @@ -15,10 +15,13 @@ #include "FairTask.h" #include "Filterable.h" +#include "NeulandPointFilter.h" +#include "R3BDataMonitor.h" #include "R3BDigitizingEngine.h" #include "R3BDigitizingPaddleNeuland.h" #include "R3BDigitizingTacQuila.h" #include "R3BDigitizingTamex.h" +#include "R3BIOConnector.h" #include "R3BNeulandGeoPar.h" #include "R3BNeulandHit.h" #include "R3BNeulandHitPar.h" @@ -60,42 +63,36 @@ class R3BNeulandDigitizer : public FairTask template using UsePaddle = Digitizing::UsePaddle; - explicit R3BNeulandDigitizer(TString input = "NeulandPoints", TString output = "NeulandHits"); - explicit R3BNeulandDigitizer(std::unique_ptr engine, - TString input = "NeulandPoints", - TString output = "NeulandHits"); - - ~R3BNeulandDigitizer() override = default; - - // No copy and no move is allowed (Rule of three/five) - R3BNeulandDigitizer(const R3BNeulandDigitizer&) = delete; // copy constructor - R3BNeulandDigitizer(R3BNeulandDigitizer&&) = delete; // move constructor - R3BNeulandDigitizer& operator=(const R3BNeulandDigitizer&) = delete; // copy assignment - R3BNeulandDigitizer& operator=(R3BNeulandDigitizer&&) = delete; // move assignment + R3BNeulandDigitizer(); + explicit R3BNeulandDigitizer(std::unique_ptr engine); protected: - InitStatus Init() override; - void Finish() override; + auto Init() -> InitStatus override; + void Finish() override { data_monitor_.save_to_sink(); } void SetParContainers() override; public: void Exec(Option_t* /*option*/) override; void SetEngine(std::unique_ptr engine); - void AddFilter(const Filterable::Filter& filter) { fHitFilters.Add(filter); } + void AddFilter(const Filterable::Filter& filter) { neuland_hit_filters_.Add(filter); } + void SetNeulandPointFilter(R3B::Neuland::BitSetParticle particle); + void SetNeulandPointFilter(R3B::Neuland::BitSetParticle particle, double minimum_allowed_energy_gev); private: - TCAInputConnector fPoints; - TCAOutputConnector fHits; + R3B::InputVectorConnector neuland_points_{ "NeulandPoints" }; + R3B::OutputVectorConnector neuland_hits_{ "NeulandHits" }; - std::unique_ptr fDigitizingEngine; // owning + std::unique_ptr digitizing_engine_; // owning - Filterable fHitFilters; + Filterable neuland_hit_filters_; - R3BNeulandGeoPar* fNeulandGeoPar = nullptr; // non-owning + R3BNeulandGeoPar* neuland_geo_par_ = nullptr; // non-owning + NeulandPointFilter neuland_point_filter_; - TH1I* hMultOne = nullptr; - TH1I* hMultTwo = nullptr; - TH1F* hRLTimeToTrig = nullptr; + R3B::DataMonitor data_monitor_; + TH1I* mult_one_ = nullptr; + TH1I* mult_two_ = nullptr; + TH1F* rl_time_to_trig_ = nullptr; public: template @@ -105,12 +102,12 @@ class R3BNeulandDigitizer : public FairTask switch (option) { case Options::neulandTamex: - fDigitizingEngine = Digitizing::CreateEngine(UsePaddle(), - UseChannel(std::forward(args)...)); + digitizing_engine_ = Digitizing::CreateEngine(UsePaddle(), + UseChannel(std::forward(args)...)); break; case Options::neulandTacquila: - fDigitizingEngine = Digitizing::CreateEngine(UsePaddle(), - UseChannel(std::forward(args)...)); + digitizing_engine_ = Digitizing::CreateEngine(UsePaddle(), + UseChannel(std::forward(args)...)); break; } } diff --git a/neuland/digitizing/R3BNeulandHitMon.cxx b/neuland/digitizing/R3BNeulandHitMon.cxx index 2fcff730b..9b3ac4223 100644 --- a/neuland/digitizing/R3BNeulandHitMon.cxx +++ b/neuland/digitizing/R3BNeulandHitMon.cxx @@ -26,10 +26,8 @@ #include #include -R3BNeulandHitMon::R3BNeulandHitMon(TString input, TString output, const Option_t* option) +R3BNeulandHitMon::R3BNeulandHitMon(const Option_t* option) : FairTask("R3B NeuLAND NeulandHit Monitor") - , fOutput(std::move(output)) - , fHits(std::move(input)) { LOG(info) << "Using R3B NeuLAND NeulandHit Monitor"; @@ -37,76 +35,119 @@ R3BNeulandHitMon::R3BNeulandHitMon(TString input, TString output, const Option_t opt.ToUpper(); if (opt.Contains("3DTRACK")) { - fIs3DTrackEnabled = true; + is_3d_track_enabled_ = true; LOG(info) << "... with 3D track visualization"; } } -InitStatus R3BNeulandHitMon::Init() +auto R3BNeulandHitMon::Init() -> InitStatus { - fHits.Init(); + neuland_hits_.init(); - if (fIs3DTrackEnabled) + if (is_3d_track_enabled_) { // XYZ -> ZXY (side view) const auto xbinN = 60; const auto ybinN = 50; const auto zbinN = 50; - fh3 = R3B::root_owned("hHits", "hHits", xbinN, 1400., 1700., ybinN, -125., 125., zbinN, -125., 125.); - fh3->SetTitle("NeuLAND Hits"); - fh3->GetXaxis()->SetTitle("Z"); - fh3->GetYaxis()->SetTitle("X"); - fh3->GetZaxis()->SetTitle("Y"); - FairRootManager::Instance()->Register("NeulandHitMon", "Hits in NeuLAND", fh3, kTRUE); + + // define additional histogram parameters + const auto x_bounds_lower = 1400.; + const auto x_bounds_upper = 1700.; + const auto yz_bounds_lower = -125.; + const auto yz_bounds_upper = 125.; + + hist_3_ = data_monitor_.add_hist("hHits", + "hHits", + xbinN, + x_bounds_lower, + x_bounds_upper, + ybinN, + yz_bounds_lower, + yz_bounds_upper, + zbinN, + yz_bounds_lower, + yz_bounds_upper); + hist_3_->SetTitle("NeuLAND Hits"); + hist_3_->GetXaxis()->SetTitle("Z"); + hist_3_->GetYaxis()->SetTitle("X"); + hist_3_->GetZaxis()->SetTitle("Y"); + FairRootManager::Instance()->Register("NeulandHitMon", "Hits in NeuLAND", hist_3_, kTRUE); } // define number of bins for histograms const auto maxHitNum = 200; const auto timeBinN = 30000; const auto zDepBinN = 60; - const auto energyBinN = 100; + const auto energyBinN = 500; const auto totenergyBinN = 1000; const auto posXYBinN = 300; const auto velocityBinN = 200; - hTime = R3B::root_owned("hTime", "Hit time", timeBinN, -1000., 1000.); - hTimeAdj = R3B::root_owned("hTimeAdj", "Hit Time adjusted for flight path", timeBinN, -1000., 1000.); - - hMult = R3B::root_owned("hMult", "Hit Multiplicity", maxHitNum, 0, maxHitNum); - - hDepth = R3B::root_owned("hDepth", "Maxial penetration depth", zDepBinN, 1400., 1700.); - hDepthVSForemostEnergy = R3B::root_owned( - "hDepthVSFrontEnergy", "Depth vs Foremost Energy", zDepBinN, 1400., 1700., energyBinN, 0., 100.); - hDepthVSSternmostEnergy = R3B::root_owned( - "hDepthVSSternmostEnergy", "Depth vs Sternmost Energy", zDepBinN, 1400., 1700., energyBinN, 0, 100.); - hDepthVSEtot = - R3B::root_owned("hDepthVSEtot", "Depth vs Total Energy", zDepBinN, 1400., 1700., totenergyBinN, 0, 1000.); - hForemostEnergy = R3B::root_owned("hForemostEnergy", "Foremost energy deposition", energyBinN, 0, 100.); - hSternmostEnergy = R3B::root_owned("hSternmostEnergy", "Sternmost energy deposition", energyBinN, 0, 100.); - hEtot = R3B::root_owned("hEtot", "Total Energy", totenergyBinN, 0, 10000.); - hdeltaEE = R3B::root_owned( - "hdeltaEE", "Energy in Foremost Plane vs Etot", energyBinN, 0, 2000., energyBinN, 0, 250.); - hPosVSEnergy = R3B::root_owned( - "hPosVSEnergy", "Position vs Energy deposition", zDepBinN, 1400., 1700., totenergyBinN, 0, 1000.); - hBeta = R3B::root_owned("hBeta", "Velocity", velocityBinN, 0., 1.); - hE = R3B::root_owned("hE", "Hit Energy", energyBinN, 0., 100.); - hX = R3B::root_owned("hX", "Hit X", posXYBinN, -150., 150.); - hY = R3B::root_owned("hY", "Hit Y", posXYBinN, -150., 150.); - hT = R3B::root_owned("hT", "Hit Delta T", timeBinN, -15., -15.); - hTNeigh = R3B::root_owned("hTNeigh", "Hit Neigh Delta T", timeBinN, -15., -15.); + // define additional histogram parameters + const auto z_dep_lower = 1400.; + const auto z_dep_upper = 1700.; + const auto total_energy_upper = 10000.; + const auto pos_xy_lower = -150.; + const auto pos_xy_upper = 150.; + const auto time_para = -15.; + const auto energy_upper_foremost_vs_e_tot_a = 2000.; + const auto energy_upper_foremost_vs_e_tot_b = 250.; + + hist_time_ = data_monitor_.add_hist("hTime", "Hit time", timeBinN, -1000., 1000.); + hist_time_adj_ = + data_monitor_.add_hist("hTimeAdj", "Hit Time adjusted for flight path", timeBinN, -1000., 1000.); + + hist_mult_ = data_monitor_.add_hist("hMult", "Hit Multiplicity", maxHitNum, 0, maxHitNum); + + hist_depth_ = + data_monitor_.add_hist("hDepth", "Maxial penetration depth", zDepBinN, z_dep_lower, z_dep_upper); + hist_depth_vs_foremost_energy_ = data_monitor_.add_hist( + "hDepthVSFrontEnergy", "Depth vs Foremost Energy", zDepBinN, z_dep_lower, z_dep_upper, energyBinN, 0., 100.); + hist_depth_vs_sternmost_energy_ = data_monitor_.add_hist("hDepthVSSternmostEnergy", + "Depth vs Sternmost Energy", + zDepBinN, + z_dep_lower, + z_dep_upper, + energyBinN, + 0, + 100.); + hist_depth_vs_energy_tot_ = data_monitor_.add_hist( + "hDepthVSEtot", "Depth vs Total Energy", zDepBinN, z_dep_lower, z_dep_upper, totenergyBinN, 0, 1000.); + hist_foremost_energy_ = + data_monitor_.add_hist("hForemostEnergy", "Foremost energy deposition", energyBinN, 0, 100.); + hist_sternmost_energy_ = + data_monitor_.add_hist("hSternmostEnergy", "Sternmost energy deposition", energyBinN, 0, 100.); + hist_energy_tot_ = data_monitor_.add_hist("hEtot", "Total Energy", totenergyBinN, 0, total_energy_upper); + hdeltaEE = data_monitor_.add_hist("hdeltaEE", + "Energy in Foremost Plane vs Etot", + energyBinN, + 0, + energy_upper_foremost_vs_e_tot_a, + energyBinN, + 0, + energy_upper_foremost_vs_e_tot_b); + hist_pos_vs_energy_ = data_monitor_.add_hist( + "hPosVSEnergy", "Position vs Energy deposition", zDepBinN, z_dep_lower, z_dep_upper, totenergyBinN, 0, 1000.); + hist_beta_ = data_monitor_.add_hist("hBeta", "Velocity", velocityBinN, 0., 1.); + hist_energy_ = data_monitor_.add_hist("hE", "Hit Energy", energyBinN, 0., 100.); + hist_x_ = data_monitor_.add_hist("hX", "Hit X", posXYBinN, pos_xy_lower, pos_xy_upper); + hist_y_ = data_monitor_.add_hist("hY", "Hit Y", posXYBinN, pos_xy_lower, pos_xy_upper); + hT = data_monitor_.add_hist("hT", "Hit Delta T", timeBinN, time_para, time_para); + hTNeigh = data_monitor_.add_hist("hTNeigh", "Hit Neigh Delta T", timeBinN, time_para, time_para); return kSUCCESS; } void R3BNeulandHitMon::Exec(Option_t* /*option*/) { - const auto hits = fHits.Retrieve(); + const auto hits = neuland_hits_.get(); // checking paddle multihits std::map paddlenum; for (const auto& hit : hits) { - auto result = paddlenum.insert(std::pair(hit->GetPaddle(), 1)); + auto result = paddlenum.insert(std::pair(hit.GetPaddle(), 1)); if (!result.second) { result.first->second++; @@ -118,101 +159,72 @@ void R3BNeulandHitMon::Exec(Option_t* /*option*/) { return (lhs.second < rhs.second); }); LOG(debug) << "max dupli: " << max->second; - if (fIs3DTrackEnabled) + if (is_3d_track_enabled_) { - fh3->Reset("ICES"); + hist_3_->Reset("ICES"); for (const auto& hit : hits) { - fh3->Fill(hit->GetPosition().Z(), hit->GetPosition().X(), hit->GetPosition().Y(), hit->GetE()); + hist_3_->Fill(hit.GetPosition().Z(), hit.GetPosition().X(), hit.GetPosition().Y(), hit.GetE()); } } - hMult->Fill(static_cast(hits.size())); + hist_mult_->Fill(static_cast(hits.size())); for (const auto& hit : hits) { - hPosVSEnergy->Fill(hit->GetPosition().Z(), hit->GetE()); - hTime->Fill(hit->GetT()); - hTimeAdj->Fill(fDistanceToTarget / hit->GetPosition().Mag() * hit->GetT()); - hBeta->Fill(hit->GetBeta()); - hE->Fill(hit->GetE()); - hX->Fill(hit->GetPosition().X()); - hY->Fill(hit->GetPosition().Y()); + hist_pos_vs_energy_->Fill(hit.GetPosition().Z(), hit.GetE()); + hist_time_->Fill(hit.GetT()); + hist_time_adj_->Fill(distance_to_target_ / hit.GetPosition().Mag() * hit.GetT()); + hist_beta_->Fill(hit.GetBeta()); + hist_energy_->Fill(hit.GetE()); + hist_x_->Fill(hit.GetPosition().X()); + hist_y_->Fill(hit.GetPosition().Y()); } for (auto it1 = hits.begin(); it1 != hits.end(); it1++) { for (auto it2 = it1 + 1; it2 != hits.end(); it2++) { - if (std::abs((*it1)->GetPosition().X() - (*it2)->GetPosition().X()) < 7.5 && - std::abs((*it1)->GetPosition().Y() - (*it2)->GetPosition().Y()) < 7.5 && - std::abs((*it1)->GetPosition().Z() - (*it2)->GetPosition().Z()) < 7.5) + const auto pos_delta_cutoff = 7.5; + if (std::abs((*it1).GetPosition().X() - (*it2).GetPosition().X()) < pos_delta_cutoff && + std::abs((*it1).GetPosition().Y() - (*it2).GetPosition().Y()) < pos_delta_cutoff && + std::abs((*it1).GetPosition().Z() - (*it2).GetPosition().Z()) < pos_delta_cutoff) { - hTNeigh->Fill((*it1)->GetT() - (*it2)->GetT()); + hTNeigh->Fill((*it1).GetT() - (*it2).GetT()); } - hT->Fill((*it1)->GetT() - (*it2)->GetT()); + hT->Fill((*it1).GetT() - (*it2).GetT()); } } auto maxDepthHit = std::max_element(hits.begin(), hits.end(), - [](R3BNeulandHit* one, R3BNeulandHit* another) - { return one->GetPosition().Z() < another->GetPosition().Z(); }); + [](const R3BNeulandHit& one, const R3BNeulandHit& another) + { return one.GetPosition().Z() < another.GetPosition().Z(); }); if (maxDepthHit != hits.end()) { - hDepth->Fill((*maxDepthHit)->GetPosition().Z()); - hSternmostEnergy->Fill((*maxDepthHit)->GetE()); - hDepthVSSternmostEnergy->Fill((*maxDepthHit)->GetPosition().Z(), (*maxDepthHit)->GetE()); + hist_depth_->Fill((*maxDepthHit).GetPosition().Z()); + hist_sternmost_energy_->Fill((*maxDepthHit).GetE()); + hist_depth_vs_sternmost_energy_->Fill((*maxDepthHit).GetPosition().Z(), (*maxDepthHit).GetE()); } auto minDepthHit = std::min_element(hits.begin(), hits.end(), - [](R3BNeulandHit* one, R3BNeulandHit* another) - { return one->GetPosition().Z() < another->GetPosition().Z(); }); + [](const R3BNeulandHit& one, const R3BNeulandHit& another) + { return one.GetPosition().Z() < another.GetPosition().Z(); }); auto Etot = - std::accumulate(hits.begin(), hits.end(), 0., [](double init, const auto* hit) { return init + hit->GetE(); }); + std::accumulate(hits.begin(), hits.end(), 0., [](double init, const auto& hit) { return init + hit.GetE(); }); if (minDepthHit != hits.end()) { - hForemostEnergy->Fill((*minDepthHit)->GetE()); - hDepthVSForemostEnergy->Fill((*maxDepthHit)->GetPosition().Z(), (*minDepthHit)->GetE()); - hdeltaEE->Fill(Etot, (*minDepthHit)->GetE()); + hist_foremost_energy_->Fill((*minDepthHit).GetE()); + hist_depth_vs_foremost_energy_->Fill((*maxDepthHit).GetPosition().Z(), (*minDepthHit).GetE()); + hdeltaEE->Fill(Etot, (*minDepthHit).GetE()); } - hEtot->Fill(Etot); + hist_energy_tot_->Fill(Etot); if (maxDepthHit != hits.end()) { - hDepthVSEtot->Fill((*maxDepthHit)->GetPosition().Z(), Etot); + hist_depth_vs_energy_tot_->Fill((*maxDepthHit).GetPosition().Z(), Etot); } } -void R3BNeulandHitMon::Finish() -{ - TDirectory* tmp = gDirectory; - FairRootManager::Instance()->GetOutFile()->cd(); - - gDirectory->mkdir(fOutput); - gDirectory->cd(fOutput); - - hDepth->Write(); - hMult->Write(); - hTime->Write(); - hTimeAdj->Write(); - hForemostEnergy->Write(); - hSternmostEnergy->Write(); - hDepthVSForemostEnergy->Write(); - hDepthVSSternmostEnergy->Write(); - hEtot->Write(); - hDepthVSEtot->Write(); - hPosVSEnergy->Write(); - hdeltaEE->Write(); - hBeta->Write(); - hE->Write(); - hX->Write(); - hY->Write(); - hT->Write(); - hTNeigh->Write(); - - gDirectory = tmp; -} - ClassImp(R3BNeulandHitMon) // NOLINT diff --git a/neuland/digitizing/R3BNeulandHitMon.h b/neuland/digitizing/R3BNeulandHitMon.h index b6f00f3bd..2b0113a28 100644 --- a/neuland/digitizing/R3BNeulandHitMon.h +++ b/neuland/digitizing/R3BNeulandHitMon.h @@ -24,8 +24,11 @@ */ #include "FairTask.h" +#include "R3BDataMonitor.h" +#include "R3BIOConnector.h" #include "R3BNeulandHit.h" #include "TCAConnector.h" +#include class TH1D; class TH2D; @@ -35,55 +38,46 @@ class TH1I; class R3BNeulandHitMon : public FairTask { public: - explicit R3BNeulandHitMon(TString input = "NeulandHits", - TString output = "NeulandHitMon", - const Option_t* option = ""); + explicit R3BNeulandHitMon(const Option_t* option = ""); - ~R3BNeulandHitMon() override = default; + void Exec(Option_t* /*option*/) override; - // No copy and no move is allowed (Rule of three/five) - R3BNeulandHitMon(const R3BNeulandHitMon&) = delete; // copy constructor - R3BNeulandHitMon(R3BNeulandHitMon&&) = delete; // move constructor - R3BNeulandHitMon& operator=(const R3BNeulandHitMon&) = delete; // copy assignment - R3BNeulandHitMon& operator=(R3BNeulandHitMon&&) = delete; // move assignment + void SetDistanceToTarget(double distance) { distance_to_target_ = distance; } protected: - InitStatus Init() override; - void Finish() override; + auto Init() -> InitStatus override; + void Finish() override { data_monitor_.save_to_sink(); } - public: - void Exec(Option_t* /*option*/) override; + private: + std::string output_{ "NeulandHitMon" }; - void SetDistanceToTarget(double distance) { fDistanceToTarget = distance; } + R3B::InputVectorConnector neuland_hits_{ "NeulandHits" }; - private: - TString fOutput; - - TCAInputConnector fHits; - - double fDistanceToTarget = 0.; - - Bool_t fIs3DTrackEnabled = false; - TH3D* fh3 = nullptr; - - TH1D* hTime = nullptr; - TH1D* hTimeAdj = nullptr; - TH1I* hMult = nullptr; - TH1D* hDepth = nullptr; - TH1D* hForemostEnergy = nullptr; - TH1D* hSternmostEnergy = nullptr; - TH2D* hDepthVSForemostEnergy = nullptr; - TH2D* hDepthVSSternmostEnergy = nullptr; - TH1D* hEtot = nullptr; - TH1D* hE = nullptr; - TH1D* hX = nullptr; - TH1D* hY = nullptr; + double distance_to_target_ = 0.; + + bool is_3d_track_enabled_ = false; + + R3B::DataMonitor data_monitor_; + + TH3D* hist_3_ = nullptr; + TH1D* hist_time_ = nullptr; + TH1D* hist_time_adj_ = nullptr; + TH1I* hist_mult_ = nullptr; + TH1D* hist_depth_ = nullptr; + TH1D* hist_foremost_energy_ = nullptr; + TH1D* hist_sternmost_energy_ = nullptr; + TH2D* hist_depth_vs_foremost_energy_ = nullptr; + TH2D* hist_depth_vs_sternmost_energy_ = nullptr; + TH1D* hist_energy_tot_ = nullptr; + TH1D* hist_energy_ = nullptr; + TH1D* hist_x_ = nullptr; + TH1D* hist_y_ = nullptr; TH1D* hT = nullptr; TH1D* hTNeigh = nullptr; - TH2D* hDepthVSEtot = nullptr; - TH2D* hPosVSEnergy = nullptr; + TH2D* hist_depth_vs_energy_tot_ = nullptr; + TH2D* hist_pos_vs_energy_ = nullptr; TH2D* hdeltaEE = nullptr; - TH1D* hBeta = nullptr; + TH1D* hist_beta_ = nullptr; ClassDefOverride(R3BNeulandHitMon, 0); // NOLINT }; diff --git a/neuland/executables/neulandAna.cxx b/neuland/executables/neulandAna.cxx index f488115f8..3f513b409 100644 --- a/neuland/executables/neulandAna.cxx +++ b/neuland/executables/neulandAna.cxx @@ -2,6 +2,7 @@ #include "FairRootFileSink.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" +#include "NeulandPointFilter.h" #include "R3BDigitizingChannelMock.h" #include "R3BDigitizingPaddleMock.h" #include "R3BDigitizingPaddleNeuland.h" @@ -135,9 +136,10 @@ auto main(int argc, const char** argv) -> int fileio2->open(paraFileName2->value().c_str()); run->GetRuntimeDb()->setSecondInput(fileio2.release()); } - auto digiNeuland = std::make_unique(); + double const minimum_filter_energy_gev = 0.000; digiNeuland->SetEngine((neulandEngines.at({ paddleName->value(), channelName->value() }))()); + digiNeuland->SetNeulandPointFilter(R3B::Neuland::BitSetParticle::none, minimum_filter_energy_gev); run->AddTask(digiNeuland.release()); auto hitmon = std::make_unique(); run->AddTask(hitmon.release()); diff --git a/neuland/executables/neulandSim.cxx b/neuland/executables/neulandSim.cxx index b8bab6821..c4cd1e7ca 100644 --- a/neuland/executables/neulandSim.cxx +++ b/neuland/executables/neulandSim.cxx @@ -22,7 +22,6 @@ constexpr int DEFAULT_RUNID = 999; int main(int argc, const char** argv) { auto timer = TStopwatch{}; - auto const PID = 2112; auto const defaultEventNum = 10; timer.Start(); diff --git a/neuland/simulation/R3BNeuland.cxx b/neuland/simulation/R3BNeuland.cxx index 018956849..33eea4cd7 100644 --- a/neuland/simulation/R3BNeuland.cxx +++ b/neuland/simulation/R3BNeuland.cxx @@ -23,27 +23,29 @@ #include "TParticle.h" #include "TVirtualMC.h" #include +#include // Initialize variables from Birk' s Law -static constexpr Double_t BirkdP = 1.032; -static constexpr Double_t BirkC1 = 0.013 / BirkdP; -static constexpr Double_t BirkC2 = 9.6e-6 / (BirkdP * BirkdP); +constexpr auto seconds_to_nanoseconds = 1e9; +constexpr auto BirkdP = 1.032; +constexpr auto BirkC1 = 0.013 / BirkdP; +constexpr auto BirkC2 = 9.6e-6 / (BirkdP * BirkdP); -inline Double_t GetLightYield(const Int_t charge, const Double_t length, const Double_t edep) +inline auto GetLightYield(const int charge, const double length, const double edep) -> double { // Apply Birk's law ( Adapted from G3BIRK/Geant3) if (charge != 0 && length > 0) { - Double_t birkC1Mod = BirkC1; + auto birkC1Mod = BirkC1; // Apply correction for higher charge states if (TMath::Abs(charge) >= 2) { - birkC1Mod *= 7.2 / 12.6; + birkC1Mod *= 7.2 / 12.6; // NOLINT } - Double_t dedxcm = 1000. * edep / length; - Double_t lightYield = edep / (1. + birkC1Mod * dedxcm + BirkC2 * dedxcm * dedxcm); + const double dedxcm = 1000. * edep / length; + const double lightYield = edep / (1. + birkC1Mod * dedxcm + BirkC2 * dedxcm * dedxcm); return lightYield; } return edep; // Rarely very small energy depositions have no length? @@ -61,96 +63,86 @@ R3BNeuland::R3BNeuland(const TString& geoFile, const TGeoTranslation& trans, con R3BNeuland::R3BNeuland(const TString& geoFile, const TGeoCombiTrans& combi) : R3BDetector("R3BNeuland", kNEULAND, geoFile, combi) - , fNeulandPoints(new TClonesArray("R3BNeulandPoint")) { } -R3BNeuland::R3BNeuland(Int_t nDP, const TGeoTranslation& trans, const TGeoRotation& rot) +R3BNeuland::R3BNeuland(int nDP, const TGeoTranslation& trans, const TGeoRotation& rot) : R3BNeuland(nDP, { trans, rot }) { } -R3BNeuland::R3BNeuland(const Int_t nDP, const TGeoCombiTrans& combi) - : R3BNeuland(TString::Format("neuland_v3_%ddp.geo.root", nDP), combi) +R3BNeuland::R3BNeuland(const int nDP, const TGeoCombiTrans& combi) + : R3BNeuland(fmt::format("neuland_v3_{}dp.geo.root", nDP), combi) { } -R3BNeuland::~R3BNeuland() -{ - if (fNeulandPoints) - { - fNeulandPoints->Delete(); - delete fNeulandPoints; - } -} - -void R3BNeuland::Initialize() -{ - LOG(info) << "R3BNeuland initialization ..."; - - FairDetector::Initialize(); - - WriteParameterFile(); - ResetValues(); -} - -Bool_t R3BNeuland::ProcessHits(FairVolume*) +auto R3BNeuland::ProcessHits(FairVolume* /*v*/) -> bool { // New hit in detector if (gMC->IsTrackEntering()) { - if (!fLastHitDone) + if (!is_last_hit_done_) { LOG(warn) << "R3BNeuland: Incomplete hit discarded"; ResetValues(); } - fLastHitDone = kFALSE; - fELoss = 0.; - fLightYield = 0.; - fTime = gMC->TrackTime() * 1.0e09; - fLength = gMC->TrackLength(); - gMC->TrackPosition(fPosIn); - gMC->TrackMomentum(fMomIn); - gMC->CurrentVolOffID(1, fPaddleID); + is_last_hit_done_ = kFALSE; + energy_loss_ = 0.; + light_yield_ = 0.; + time_ = gMC->TrackTime() * seconds_to_nanoseconds; + length_ = gMC->TrackLength(); + gMC->TrackPosition(pos_in_); + gMC->TrackMomentum(mom_in_); + gMC->CurrentVolOffID(1, paddle_id_); + + particle_id_ = gMC->TrackPid(); + trackid_pid_map_.emplace(gMC->GetStack()->GetCurrentTrackNumber(), gMC->TrackPid()); + if (auto search = trackid_pid_map_.find(gMC->GetStack()->GetCurrentParentTrackNumber()); + search != trackid_pid_map_.end()) + { + partent_particle_id_ = search->first; + } } // Sum energy loss for all steps in the active volume - fELoss += gMC->Edep(); - fLightYield += GetLightYield(gMC->TrackCharge(), gMC->TrackStep(), gMC->Edep()); + energy_loss_ += gMC->Edep(); + light_yield_ += GetLightYield(gMC->TrackCharge(), gMC->TrackStep(), gMC->Edep()); // Set additional parameters at exit of active volume. Create R3BNeulandPoint. if (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) { // Do not save a hit if no energy deposited - if (fELoss < 1e-20 || fLightYield < 1e-20) + constexpr auto minimum_energy_cutoff = 1e-20; + if (energy_loss_ < minimum_energy_cutoff || light_yield_ < minimum_energy_cutoff) { ResetValues(); return kTRUE; } - fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); - gMC->TrackPosition(fPosOut); - gMC->TrackMomentum(fMomOut); + track_id_ = gMC->GetStack()->GetCurrentTrackNumber(); + gMC->TrackPosition(pos_out_); + gMC->TrackMomentum(mom_out_); // Add Point - LOG(debug) << "R3BNeuland: Adding Point at (" << fPosIn.X() << ", " << fPosIn.Y() << ", " << fPosIn.Z() - << ") cm, paddle " << fPaddleID << ", track " << fTrackID << ", energy loss " << fELoss << " GeV " - << gMC->GetStack()->GetCurrentParentTrackNumber(); - - Int_t size = fNeulandPoints->GetEntriesFast(); - new ((*fNeulandPoints)[size]) R3BNeulandPoint(fTrackID, - fPaddleID, - fPosIn.Vect(), - fMomIn.Vect(), - fTime, - fLength, - fELoss, - gMC->CurrentEvent(), - fLightYield); + LOG(debug) << "R3BNeuland: Adding Point at (" << pos_in_.X() << ", " << pos_in_.Y() << ", " << pos_in_.Z() + << ") cm, paddle " << paddle_id_ << ", track " << track_id_ << ", energy loss " << energy_loss_ + << " GeV " << gMC->GetStack()->GetCurrentParentTrackNumber(); + + neuland_points_.get().emplace_back(track_id_, + paddle_id_, + pos_in_.Vect(), + mom_in_.Vect(), + time_, + length_, + energy_loss_, + gMC->CurrentEvent(), + light_yield_, + particle_id_, + partent_particle_id_); // Increment number of LandPoints for this track - auto stack = dynamic_cast(gMC->GetStack()); + auto* stack = dynamic_cast(gMC->GetStack()); stack->AddPoint(kNEULAND); ResetValues(); } @@ -158,62 +150,49 @@ Bool_t R3BNeuland::ProcessHits(FairVolume*) return kTRUE; } -Bool_t R3BNeuland::CheckIfSensitive(std::string name) { return name == "volBC408"; } +auto R3BNeuland::CheckIfSensitive(std::string name) -> bool { return name == "volBC408"; } void R3BNeuland::EndOfEvent() { - if (fVerboseLevel) + if (fVerboseLevel != 0) { Print(); } Reset(); } -TClonesArray* R3BNeuland::GetCollection(Int_t iColl) const +void R3BNeuland::Print(Option_t* /*unused*/) const { - if (iColl == 0) - { - return fNeulandPoints; - } - return nullptr; -} - -void R3BNeuland::Register() -{ - FairRootManager::Instance()->Register("NeulandPoints", GetName(), fNeulandPoints, kTRUE); -} - -void R3BNeuland::Print(Option_t*) const -{ - LOG(info) << "R3BNeuland: " << fNeulandPoints->GetEntries() << " Neuland Points registered in this event"; + LOG(info) << "R3BNeuland: " << neuland_points_.get_constref().size() << " Neuland Points registered in this event"; } void R3BNeuland::Reset() { - fNeulandPoints->Clear(); + neuland_points_.clear(); ResetValues(); + trackid_pid_map_.clear(); } void R3BNeuland::ResetValues() { - fLastHitDone = kTRUE; - fTrackID = 0; - fPaddleID = -1; - fPosIn.Clear(); - fPosOut.Clear(); - fMomIn.Clear(); - fMomOut.Clear(); - fTime = fLength = fELoss = fLightYield = 0; + is_last_hit_done_ = kTRUE; + track_id_ = 0; + paddle_id_ = -1; + pos_in_.Clear(); + pos_out_.Clear(); + mom_in_.Clear(); + mom_out_.Clear(); + time_ = length_ = energy_loss_ = light_yield_ = 0; } void R3BNeuland::WriteParameterFile() { FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb(); - fNeulandGeoPar = dynamic_cast(rtdb->getContainer("R3BNeulandGeoPar")); + neuland_geo_par_ = dynamic_cast(rtdb->getContainer("R3BNeulandGeoPar")); // Really bad way to find the Neuland *node* (not the volume!) TGeoNode* geoNodeNeuland = nullptr; - for (Int_t i = 0; i < gGeoManager->GetTopNode()->GetNdaughters(); i++) + for (int i = 0; i < gGeoManager->GetTopNode()->GetNdaughters(); i++) { if (TString(gGeoManager->GetTopNode()->GetDaughter(i)->GetVolume()->GetName()) == "volNeuland") { @@ -222,13 +201,23 @@ void R3BNeuland::WriteParameterFile() } } - if (!geoNodeNeuland) + if (geoNodeNeuland == nullptr) { LOG(fatal) << "volNeuland not found"; } - fNeulandGeoPar->SetNeulandGeoNode(geoNodeNeuland); - fNeulandGeoPar->setChanged(); + neuland_geo_par_->SetNeulandGeoNode(geoNodeNeuland); + neuland_geo_par_->setChanged(); +} + +void R3BNeuland::Register() +{ + LOG(info) << "R3BNeuland initialization ..."; + + neuland_points_.init(); + + WriteParameterFile(); + ResetValues(); } ClassImp(R3BNeuland); diff --git a/neuland/simulation/R3BNeuland.h b/neuland/simulation/R3BNeuland.h index 6dee05be7..7e92a7d93 100644 --- a/neuland/simulation/R3BNeuland.h +++ b/neuland/simulation/R3BNeuland.h @@ -15,6 +15,8 @@ #define R3BNEULAND_H #include "R3BDetector.h" +#include "R3BIOConnector.h" +#include "R3BNeulandPoint.h" #include "TLorentzVector.h" #include #include @@ -53,52 +55,46 @@ class R3BNeuland : public R3BDetector *@param nDP number of double planes *@param trans position *@param rot rotation */ - R3BNeuland(Int_t nDP, const TGeoTranslation& trans, const TGeoRotation& rot = TGeoRotation()); + R3BNeuland(int nDP, const TGeoTranslation& trans, const TGeoRotation& rot = TGeoRotation()); /** Standard constructor. *@param nDP number of double planes *@param combi position + rotation */ - explicit R3BNeuland(Int_t nDP, const TGeoCombiTrans& combi = TGeoCombiTrans()); + explicit R3BNeuland(int nDP, const TGeoCombiTrans& combi = TGeoCombiTrans()); - /** Default Destructor */ - ~R3BNeuland() override; - - void Initialize() override; - - Bool_t ProcessHits(FairVolume* = nullptr) override; + auto ProcessHits(FairVolume* /*v*/ = nullptr) -> bool override; void EndOfEvent() override; - void Register() override; - - TClonesArray* GetCollection(Int_t iColl) const override; - - void Print(Option_t* = "") const override; + void Print(Option_t* /*unused*/ = "") const override; void Reset() override; - Bool_t CheckIfSensitive(std::string name) override; + auto CheckIfSensitive(std::string name) -> bool override; + + [[nodiscard]] auto GetCollection(Int_t /*iColl*/) const -> TClonesArray* override { return nullptr; } - // No copy and no move is allowed (Rule of three/five) - R3BNeuland(const R3BNeuland&) = delete; // copy constructor - R3BNeuland(R3BNeuland&&) = delete; // move constructor - R3BNeuland& operator=(const R3BNeuland&) = delete; // copy assignment - R3BNeuland& operator=(R3BNeuland&&) = delete; // move assignment + void Register() override; private: - TClonesArray* fNeulandPoints; //! - R3BNeulandGeoPar* fNeulandGeoPar; //! + R3B::OutputVectorConnector neuland_points_{ "NeulandPoints" }; //! + R3BNeulandGeoPar* neuland_geo_par_ = nullptr; //! + std::map trackid_pid_map_; /** Track information to be stored until the track leaves the active volume. */ - Int_t fTrackID; - Int_t fPaddleID; - TLorentzVector fPosIn, fPosOut; - TLorentzVector fMomIn, fMomOut; - Double_t fTime; - Double_t fLength; - Double_t fELoss; - Double_t fLightYield; - Bool_t fLastHitDone; + int track_id_ = 0; + int paddle_id_ = 0; + TLorentzVector pos_in_{}; + TLorentzVector pos_out_{}; + TLorentzVector mom_in_{}; + TLorentzVector mom_out_{}; + double time_ = 0.; + double length_ = 0.; + double energy_loss_ = 0.; + double light_yield_ = 0.; + bool is_last_hit_done_ = false; + int particle_id_ = 0; + int partent_particle_id_ = 0; void ResetValues(); diff --git a/neuland/test/CMakeLists.txt b/neuland/test/CMakeLists.txt index f00ea8228..3242d0359 100644 --- a/neuland/test/CMakeLists.txt +++ b/neuland/test/CMakeLists.txt @@ -11,18 +11,20 @@ # or submit itself to any jurisdiction. # ############################################################################## -if(GTEST_FOUND) +if(GTest_FOUND) + include(GoogleTest) add_executable( NeulandUnitTests digitizing/testNeulandDigitizingPaddle.cxx digitizing/testNeulandDigitizingTamex.cxx + digitizing/testNeulandPointFilter.cxx NeulandUnitTests.cxx testClusteringEngine.cxx testNeulandMultiplicityCalorimetricPar.cxx) target_link_libraries(NeulandUnitTests PRIVATE GTest::gtest_main GTest::gmock_main R3BNeulandDigitizing R3BNeulandReconstruction) gtest_discover_tests(NeulandUnitTests DISCOVERY_TIMEOUT 600) -endif(GTEST_FOUND) +endif(GTest_FOUND) set(simuPars --simuFile diff --git a/neuland/test/digitizing/testNeulandPointFilter.cxx b/neuland/test/digitizing/testNeulandPointFilter.cxx new file mode 100644 index 000000000..43855f4c2 --- /dev/null +++ b/neuland/test/digitizing/testNeulandPointFilter.cxx @@ -0,0 +1,110 @@ +/****************************************************************************** + * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * + * Copyright (C) 2019-2023 Members of R3B Collaboration * + * * + * This software is distributed under the terms of the * + * GNU General Public Licence (GPL) version 3, * + * copied verbatim in the file "LICENSE". * + * * + * In applying this license GSI does not waive the privileges and immunities * + * granted to it by virtue of its status as an Intergovernmental Organization * + * or submit itself to any jurisdiction. * + ******************************************************************************/ + +#include "NeulandPointFilter.h" +#include "R3BNeulandPoint.h" +#include +#include + +namespace +{ + + namespace neuland = ::R3B::Neuland; + + using BitSetParticle = R3B::Neuland::BitSetParticle; + constexpr auto ParticleBitsetSize = R3B::Neuland::ParticleBitsetSize; + // Test for ParticleToBitSet function + + TEST(TestNeulandPointFilter, ParticleToBitSetTest) + { + auto bitset = neuland::ParticleToBitSet(BitSetParticle::proton); + EXPECT_EQ(bitset.to_ulong(), 0x0001); + + bitset = ParticleToBitSet(BitSetParticle::neutron); + EXPECT_EQ(bitset.to_ulong(), 0x0002); + + bitset = ParticleToBitSet(BitSetParticle::gamma); + EXPECT_EQ(bitset.to_ulong(), 0x0020); + } + + // Test for BitSetToParticle function + TEST(TestNeulandPointFilter, BitSetToParticleTest) + { + constexpr auto proton_bitsetid = 0x0001; + auto particle = neuland::BitSetToParticle(std::bitset{ proton_bitsetid }); + EXPECT_EQ(particle, BitSetParticle::proton); + + constexpr auto neutron_bitsetid = 0x0002; + particle = neuland::BitSetToParticle(std::bitset{ neutron_bitsetid }); + EXPECT_EQ(particle, BitSetParticle::neutron); + + constexpr auto gamma_bitsetid = 0x0020; + particle = neuland::BitSetToParticle(std::bitset{ gamma_bitsetid }); + EXPECT_EQ(particle, BitSetParticle::gamma); + } + + // Test for CheckCriteria function + TEST(TestNeulandPointFilter, CheckCriteriaTest) + { + EXPECT_TRUE(CheckCriteria(BitSetParticle::proton, BitSetParticle::proton)); + EXPECT_FALSE(CheckCriteria(BitSetParticle::proton, BitSetParticle::neutron)); + EXPECT_TRUE(CheckCriteria(BitSetParticle::proton, BitSetParticle::neutron | BitSetParticle::proton)); + } + + // Test for PidToBitSetParticle function + TEST(TestNeulandPointFilter, PidToBitSetParticleTest) + { + EXPECT_EQ(neuland::PidToBitSetParticle(2212), BitSetParticle::proton); + EXPECT_EQ(neuland::PidToBitSetParticle(2112), BitSetParticle::neutron); + EXPECT_EQ(neuland::PidToBitSetParticle(22), BitSetParticle::gamma); + } + + // Test for NeulandPointFilter class + TEST(TestNeulandPointFilter, SetFilterTest) + { + NeulandPointFilter filter; + filter.SetFilter(BitSetParticle::proton); + EXPECT_EQ(filter.GetFilter(), BitSetParticle::proton); + + filter.SetFilter(BitSetParticle::neutron, 0.5); + EXPECT_EQ(filter.GetFilter(), BitSetParticle::neutron); + EXPECT_EQ(filter.GetMinimumAllowedEnergy(), 0.5); + } + + // Test for NeulandPointFilter::ShouldNeulandPointBeFiltered function + TEST(NeulandPointFilterTest, ShouldNeulandPointBeFilteredTest) + { + constexpr auto proton_pid = 2212; + constexpr auto neutron_pid = 2112; + constexpr auto eLoss_proton = 0.7; + constexpr auto eLoss_neutron = 0.3; + const TVector3 sample_vector; + + NeulandPointFilter filter; + + const R3BNeulandPoint protonPoint( + 0, 0, sample_vector, sample_vector, 0., 0., eLoss_proton, 0, 0., proton_pid, 0); + const R3BNeulandPoint neutronPoint( + 0, 0, sample_vector, sample_vector, 0., 0., eLoss_neutron, 0, 0., neutron_pid, 0); + + // Test filtering criteria for protons + filter.SetFilter(BitSetParticle::proton); + EXPECT_TRUE(filter.CheckFiltered(protonPoint)); + EXPECT_FALSE(filter.CheckFiltered(neutronPoint)); + + // Test minimum energy filter + filter.SetFilter(BitSetParticle::neutron, 0.5); + EXPECT_TRUE(filter.CheckFiltered(neutronPoint)); + EXPECT_FALSE(filter.CheckFiltered(protonPoint)); + } +} // namespace diff --git a/neuland/test/testNeulandDigitizer.C b/neuland/test/testNeulandDigitizer.C index 09426b13f..70597ba40 100644 --- a/neuland/test/testNeulandDigitizer.C +++ b/neuland/test/testNeulandDigitizer.C @@ -27,9 +27,9 @@ void testNeulandDigitizer() run.GetRuntimeDb()->setFirstInput(io); run.AddTask(new R3BNeulandDigitizer(R3BNeulandDigitizer::Options::neulandTamex)); - run.AddTask(new R3BNeulandClusterFinder()); - run.AddTask(new R3BNeulandPrimaryInteractionFinder()); - run.AddTask(new R3BNeulandPrimaryClusterFinder()); + // run.AddTask(new R3BNeulandClusterFinder()); + // run.AddTask(new R3BNeulandPrimaryInteractionFinder()); + // run.AddTask(new R3BNeulandPrimaryClusterFinder()); run.Init(); run.Run(0, 0); diff --git a/r3bbase/R3BIOConnector.h b/r3bbase/R3BIOConnector.h index 9eaa2b1c6..791ed195f 100644 --- a/r3bbase/R3BIOConnector.h +++ b/r3bbase/R3BIOConnector.h @@ -151,6 +151,7 @@ namespace R3B } [[nodiscard]] inline auto get() -> RawDataType& { return data_; } + [[nodiscard]] inline auto get_constref() const -> const RawDataType& { return data_; } inline void clear() { data_.clear(); } diff --git a/r3bdata/DataLinkDef.h b/r3bdata/DataLinkDef.h index 4ca0d06b8..7d045a96e 100755 --- a/r3bdata/DataLinkDef.h +++ b/r3bdata/DataLinkDef.h @@ -109,6 +109,7 @@ #pragma link C++ class R3BPspDigi+; #pragma link C++ class R3BNeulandTacquilaMappedData+; #pragma link C++ class R3BNeulandPoint+; +#pragma link C++ class std::vector+; #pragma link C++ class R3BNeulandHit+; #pragma link C++ class R3BNeulandMultiplicity+; #pragma link C++ class R3BNeulandCluster+; diff --git a/r3bdata/neulandData/R3BNeulandPoint.h b/r3bdata/neulandData/R3BNeulandPoint.h index 2687a8423..5252f5573 100644 --- a/r3bdata/neulandData/R3BNeulandPoint.h +++ b/r3bdata/neulandData/R3BNeulandPoint.h @@ -24,11 +24,12 @@ class R3BNeulandPoint : public FairMCPoint public: R3BNeulandPoint() - : FairMCPoint() - , fLightYield(0) + : fLightYield{ 0 } + , particle_id_{ 0 } + , parent_particle_id_{ 0 } { } - + // NOLINTBEGIN R3BNeulandPoint(const Int_t trackID, const Int_t detID, const TVector3& pos, @@ -37,32 +38,44 @@ class R3BNeulandPoint : public FairMCPoint const Double_t length, const Double_t eLoss, const UInt_t EventId, - const Double_t lightYield) - : FairMCPoint(trackID, detID, pos, mom, tof, length, eLoss, EventId) - , fLightYield(lightYield) + const Double_t lightYield, + const int particle_id, + const int parent_particle_id) // NOLINTEND + : FairMCPoint{ trackID, detID, pos, mom, tof, length, eLoss, EventId } + , fLightYield{ lightYield } + , particle_id_{ particle_id } + , parent_particle_id_{ parent_particle_id } { } - R3BNeulandPoint(const FairMCPoint& point, const Double_t lightYield) - : FairMCPoint(point) - , fLightYield(lightYield) + R3BNeulandPoint(const FairMCPoint& point, + const Double_t lightYield, + const int particle_id, + const int parent_particle_id) + : FairMCPoint{ point } + , fLightYield{ lightYield } + , particle_id_{ particle_id } + , parent_particle_id_{ parent_particle_id } { } - TVector3 GetMomentum() const; - TVector3 GetPosition() const; - Int_t GetPaddle() const { return GetDetectorID(); } - Double_t GetLightYield() const { return fLightYield; } + [[nodiscard]] auto GetMomentum() const -> TVector3; + [[nodiscard]] auto GetPosition() const -> TVector3; + [[nodiscard]] auto GetPaddle() const -> Int_t { return GetDetectorID(); } + [[nodiscard]] auto GetLightYield() const -> Double_t { return fLightYield; } + [[nodiscard]] auto GetPID() const -> int { return particle_id_; } + ClassDefOverride(R3BNeulandPoint, 2); - void Print(const Option_t*) const override; + void Print(const Option_t* /*opt*/) const override; protected: Double_t fLightYield; - public: - ClassDefOverride(R3BNeulandPoint, 1) + private: + int particle_id_; + int parent_particle_id_; }; -std::ostream& operator<<(std::ostream&, const R3BNeulandPoint&); // Support easy printing +auto operator<<(std::ostream&, const R3BNeulandPoint&) -> std::ostream&; // Support easy printing #endif // R3BNEULANDPOINT_H