diff --git a/Analysis/CMakeLists.txt b/Analysis/CMakeLists.txt index 922465bb4fc3c..bbca46e0aca02 100644 --- a/Analysis/CMakeLists.txt +++ b/Analysis/CMakeLists.txt @@ -14,3 +14,4 @@ add_subdirectory(Tasks) add_subdirectory(Tutorials) add_subdirectory(PWGDQ) add_subdirectory(ALICE3) +add_subdirectory(EventFiltering) diff --git a/Analysis/EventFiltering/CMakeLists.txt b/Analysis/EventFiltering/CMakeLists.txt new file mode 100644 index 0000000000000..d23b2879fae2b --- /dev/null +++ b/Analysis/EventFiltering/CMakeLists.txt @@ -0,0 +1,24 @@ +# Copyright CERN and copyright holders of ALICE O2. This software is distributed +# under the terms of the GNU General Public License v3 (GPL Version 3), copied +# verbatim in the file "COPYING". +# +# See http://alice-o2.web.cern.ch/license for full licensing information. +# +# In applying this license CERN 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. + +o2_add_library( + EventFiltering + SOURCES centralEventFilterProcessor.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::AnalysisDataModel O2::AnalysisCore) + +o2_add_executable(central-event-filter-processor + SOURCES cefp.cxx + COMPONENT_NAME Analysis + PUBLIC_LINK_LIBRARIES O2::EventFiltering) + +o2_add_dpl_workflow(nuclei-filter + SOURCES nucleiFilter.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::AnalysisDataModel O2::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/Analysis/EventFiltering/cefp.cxx b/Analysis/EventFiltering/cefp.cxx new file mode 100644 index 0000000000000..2f9abbfc3850f --- /dev/null +++ b/Analysis/EventFiltering/cefp.cxx @@ -0,0 +1,42 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN 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 "CommonUtils/ConfigurableParam.h" +#include "centralEventFilterProcessor.h" + +using namespace o2::framework; + +// ------------------------------------------------------------------ + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + workflowOptions.push_back(ConfigParamSpec{"config", o2::framework::VariantType::String, "train_config.json", {"Configuration of the filtering"}}); +} + +// ------------------------------------------------------------------ + +#include "Framework/runDataProcessing.h" +#include "Framework/Logger.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + auto config = configcontext.options().get("config"); + + if (config.empty()) { + LOG(FATAL) << "We need a configuration for the centralEventFilterProcessor"; + throw std::runtime_error("incompatible options provided"); + } + + WorkflowSpec specs; + specs.emplace_back(o2::aod::filtering::getCentralEventFilterProcessorSpec(config)); + return std::move(specs); +} diff --git a/Analysis/EventFiltering/centralEventFilterProcessor.cxx b/Analysis/EventFiltering/centralEventFilterProcessor.cxx new file mode 100644 index 0000000000000..867bbc4a1ccfa --- /dev/null +++ b/Analysis/EventFiltering/centralEventFilterProcessor.cxx @@ -0,0 +1,189 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN 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. + +/// @file centralEventFilterProcessor.cxx + +#include "centralEventFilterProcessor.h" + +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/TableConsumer.h" + +#include +#include + +#include +#include +#include +#include +#include + +using namespace o2::framework; +using namespace rapidjson; + +namespace +{ +Document readJsonFile(std::string& config) +{ + FILE* fp = fopen(config.data(), "rb"); + + char readBuffer[65536]; + FileReadStream is(fp, readBuffer, sizeof(readBuffer)); + + Document d; + d.ParseStream(is); + fclose(fp); + return d; +} +} // namespace + +namespace o2::aod::filtering +{ + +void CentralEventFilterProcessor::init(framework::InitContext& ic) +{ + // JSON example + // { + // "subwagon_name" : "CentralEventFilterProcessor", + // "configuration" : { + // "NucleiFilters" : { + // "H2" : 0.1, + // "H3" : 0.3, + // "HE3" : 1., + // "HE4" : 1. + // } + // } + // } + std::cout << "Start init" << std::endl; + Document d = readJsonFile(mConfigFile); + int nCols{0}; + for (auto& workflow : d["workflows"].GetArray()) { + if (std::string_view(workflow["subwagon_name"].GetString()) == "CentralEventFilterProcessor") { + auto& config = workflow["configuration"]; + for (auto& filter : AvailableFilters) { + auto& filterConfig = config[filter]; + if (config.IsObject()) { + std::unordered_map tableMap; + for (auto& node : filterConfig.GetObject()) { + tableMap[node.name.GetString()] = node.value.GetDouble(); + nCols++; + } + mDownscaling[filter] = tableMap; + } + } + break; + } + } + std::cout << "Middle init" << std::endl; + mScalers = new TH1D("mScalers", ";;Number of events", nCols + 1, -0.5, 0.5 + nCols); + mScalers->GetXaxis()->SetBinLabel(1, "Total number of events"); + + mFiltered = new TH1D("mFiltered", ";;Number of filtered events", nCols + 1, -0.5, 0.5 + nCols); + mFiltered->GetXaxis()->SetBinLabel(1, "Total number of events"); + + int bin{2}; + for (auto& table : mDownscaling) { + for (auto& column : table.second) { + mScalers->GetXaxis()->SetBinLabel(bin, column.first.data()); + mFiltered->GetXaxis()->SetBinLabel(bin++, column.first.data()); + } + } + + TFile test("test.root", "recreate"); + mScalers->Clone()->Write(); + test.Close(); +} + +void CentralEventFilterProcessor::run(ProcessingContext& pc) +{ + int64_t nEvents{-1}; + for (auto& tableName : mDownscaling) { + auto tableConsumer = pc.inputs().get(tableName.first); + + auto tablePtr{tableConsumer->asArrowTable()}; + int64_t nRows{tablePtr->num_rows()}; + nEvents = nEvents < 0 ? nRows : nEvents; + if (nEvents != nRows) { + std::cerr << "Inconsistent number of rows across trigger tables, fatal" << std::endl; ///TODO: move it to real fatal + } + + auto schema{tablePtr->schema()}; + for (auto& colName : tableName.second) { + int bin{mScalers->GetXaxis()->FindBin(colName.first.data())}; + double binCenter{mScalers->GetXaxis()->GetBinCenter(bin)}; + auto column{tablePtr->GetColumnByName(colName.first)}; + double downscaling{colName.second}; + if (column) { + for (int64_t iC{0}; iC < column->num_chunks(); ++iC) { + auto chunk{column->chunk(iC)}; + auto boolArray = std::static_pointer_cast(chunk); + for (int64_t iS{0}; iS < chunk->length(); ++iS) { + if (boolArray->Value(iS)) { + mScalers->Fill(binCenter); + if (mUniformGenerator(mGeneratorEngine) < downscaling) { + mFiltered->Fill(binCenter); + } + } + } + } + } + } + } + mScalers->SetBinContent(1, mScalers->GetBinContent(1) + nEvents); + mFiltered->SetBinContent(1, mFiltered->GetBinContent(1) + nEvents); +} + +void CentralEventFilterProcessor::endOfStream(EndOfStreamContext& ec) +{ + TFile output("trigger.root", "recreate"); + mScalers->Write("Scalers"); + mFiltered->Write("FilteredScalers"); + if (mScalers->GetBinContent(1) > 1.e-24) { + mScalers->Scale(1. / mScalers->GetBinContent(1)); + } + if (mFiltered->GetBinContent(1) > 1.e-24) { + mFiltered->Scale(1. / mFiltered->GetBinContent(1)); + } + mScalers->Write("Fractions"); + mFiltered->Write("FractionsDownscaled"); + output.Close(); +} + +DataProcessorSpec getCentralEventFilterProcessorSpec(std::string& config) +{ + + std::vector inputs; + Document d = readJsonFile(config); + + for (auto& workflow : d["workflows"].GetArray()) { + for (unsigned int iFilter{0}; iFilter < AvailableFilters.size(); ++iFilter) { + if (std::string_view(workflow["subwagon_name"].GetString()) == std::string_view(AvailableFilters[iFilter])) { + inputs.emplace_back(std::string(AvailableFilters[iFilter]), "AOD", FilterDescriptions[iFilter], 0, Lifetime::Timeframe); + std::cout << "Adding input " << std::string_view(AvailableFilters[iFilter]) << std::endl; + break; + } + } + } + + std::vector outputs; + outputs.emplace_back("AOD", "Decision", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "o2-central-event-filter-processor", + inputs, + outputs, + AlgorithmSpec{adaptFromTask(config)}, + Options{ + {"filtering-config", VariantType::String, "", {"Path to the filtering json config file"}}}}; +} + +} // namespace o2::aod::filtering + +/// o2-analysis-cefp --config /Users/stefano.trogolo/alice/run3/O2/Analysis/EventFiltering/sample.json | o2-analysis-nuclei-filter --aod-file @listOfFiles.txt | o2-analysis-trackselection | o2-analysis-trackextension | o2-analysis-pid-tpc | o2-analysis-pid-tof | o2-analysis-event-selection | o2-analysis \ No newline at end of file diff --git a/Analysis/EventFiltering/centralEventFilterProcessor.h b/Analysis/EventFiltering/centralEventFilterProcessor.h new file mode 100644 index 0000000000000..ad2e671902099 --- /dev/null +++ b/Analysis/EventFiltering/centralEventFilterProcessor.h @@ -0,0 +1,50 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN 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. + +/// @file centralEventFilterProcessor.h + +#ifndef O2_CentralEventFilterProcessor +#define O2_CentralEventFilterProcessor + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" + +#include +#include "filterTables.h" + +class TH1D; + +namespace o2::aod::filtering +{ + +class CentralEventFilterProcessor : public framework::Task +{ + public: + CentralEventFilterProcessor(const std::string& config) : mConfigFile{config} {} + ~CentralEventFilterProcessor() override = default; + void init(framework::InitContext& ic) final; + void run(framework::ProcessingContext& pc) final; + void endOfStream(framework::EndOfStreamContext& ec) final; + + private: + TH1D* mScalers; + TH1D* mFiltered; + std::string mConfigFile; + std::unordered_map> mDownscaling; + std::mt19937_64 mGeneratorEngine; + std::uniform_real_distribution mUniformGenerator = std::uniform_real_distribution(0., 1.); +}; + +/// create a processor spec +framework::DataProcessorSpec getCentralEventFilterProcessorSpec(std::string& config); + +} // namespace o2::aod::filtering + +#endif /* O2_CentralEventFilterProcessor */ diff --git a/Analysis/EventFiltering/filterTables.h b/Analysis/EventFiltering/filterTables.h new file mode 100644 index 0000000000000..6868b69990e73 --- /dev/null +++ b/Analysis/EventFiltering/filterTables.h @@ -0,0 +1,35 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN 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. +#ifndef O2_ANALYSIS_TRIGGER_H_ +#define O2_ANALYSIS_TRIGGER_H_ + +#include +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ +namespace filtering +{ +DECLARE_SOA_COLUMN(H2, hasH2, bool); +DECLARE_SOA_COLUMN(H3, hasH3, bool); +DECLARE_SOA_COLUMN(He3, hasHe3, bool); +DECLARE_SOA_COLUMN(He4, hasHe4, bool); + +} // namespace filtering + +DECLARE_SOA_TABLE(NucleiFilters, "AOD", "Nuclei Filters", filtering::H2, filtering::H3, filtering::He3, filtering::He4); + +constexpr std::array AvailableFilters{"NucleiFilters"}; +constexpr std::array FilterDescriptions{"Nuclei Filters"}; + +using NucleiFilter = NucleiFilters::iterator; +} // namespace o2::aod + +#endif // O2_ANALYSIS_TRIGGER_H_ diff --git a/Analysis/EventFiltering/nucleiFilter.cxx b/Analysis/EventFiltering/nucleiFilter.cxx new file mode 100644 index 0000000000000..30620808d9879 --- /dev/null +++ b/Analysis/EventFiltering/nucleiFilter.cxx @@ -0,0 +1,139 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN 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. +// O2 includes + +#include "ReconstructionDataFormats/Track.h" +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "AnalysisDataModel/PID/PIDResponse.h" +#include "AnalysisDataModel/TrackSelectionTables.h" + +#include "AnalysisDataModel/EventSelection.h" +#include "filterTables.h" + +#include "Framework/HistogramRegistry.h" + +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace +{ +float rapidity(float pt, float eta, float m) +{ + return std::asinh(pt / std::hypot(m, pt) * std::sinh(eta)); +} + +static constexpr int nNuclei{4}; +static constexpr int nCutsPID{5}; +static constexpr std::array masses{ + constants::physics::MassDeuteron, constants::physics::MassTriton, + constants::physics::MassHelium3, constants::physics::MassAlpha}; +static constexpr std::array charges{1, 1, 2, 2}; +static const std::vector nucleiNames{"H2", "H3", "He3", "He4"}; +static const std::vector cutsNames{ + "TPCnSigmaMin", "TPCnSigmaMax", "TOFnSigmaMin", "TOFnSigmaMax", "TOFpidStartPt"}; +static constexpr float cutsPID[nNuclei][nCutsPID]{ + {-3.f, +3.f, -4.f, +4.f, 1.0f}, /*H2*/ + {-3.f, +3.f, -4.f, +4.f, 1.6f}, /*H3*/ + {-5.f, +5.f, -4.f, +4.f, 14000.f}, /*He3*/ + {-5.f, +5.f, -4.f, +4.f, 14000.f} /*He4*/ +}; +} // namespace + +struct nucleiFilter { + + Produces tags; + + Configurable yMin{"yMin", -0.8, "Maximum rapidity"}; + Configurable yMax{"yMax", 0.8, "Minimum rapidity"}; + Configurable yBeam{"yBeam", 0., "Beam rapidity"}; + + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + + Configurable> cfgCutsPID{"nucleiCutsPID", {cutsPID[0], nNuclei, nCutsPID, nucleiNames, cutsNames}, "Nuclei PID selections"}; + + HistogramRegistry spectra{"spectra", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + void init(o2::framework::InitContext&) + { + std::vector ptBinning = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4., 5.}; + std::vector centBinning = {0., 1., 5., 10., 20., 30., 40., 50., 70., 100.}; + + AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec centAxis = {centBinning, "V0M (%)"}; + + spectra.add("fCollZpos", "collision z position", HistType::kTH1F, {{600, -20., +20., "z position (cm)"}}); + spectra.add("fTPCsignal", "Specific energy loss", HistType::kTH2F, {{600, 0., 3, "#it{p} (GeV/#it{c})"}, {1400, 0, 1400, "d#it{E} / d#it{X} (a. u.)"}}); + spectra.add("fTPCcounts", "n-sigma TPC", HistType::kTH2F, {ptAxis, {200, -100., +100., "n#sigma_{He} (a. u.)"}}); + spectra.add("fProcessedEvents", "Nuclei - event filtered", HistType::kTH1F, {{4, -0.5, 3.5, "Event counter"}}); + } + + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::isGlobalTrack == static_cast(1u)); + + using TrackCandidates = soa::Filtered>; + void process(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, TrackCandidates const& tracks) + { + // collision process loop + bool keepEvent[nNuclei]{false}; + // + spectra.fill(HIST("fCollZpos"), collision.posZ()); + // + + for (auto track : tracks) { // start loop over tracks + + const float nSigmaTPC[nNuclei]{ + track.tpcNSigmaDe(), track.tpcNSigmaTr(), track.tpcNSigmaHe(), track.tpcNSigmaAl()}; + const float nSigmaTOF[nNuclei]{ + track.tofNSigmaDe(), track.tofNSigmaTr(), track.tofNSigmaHe(), track.tofNSigmaAl()}; + + for (int iN{0}; iN < nNuclei; ++iN) { + float y{rapidity(track.pt() * charges[iN], track.eta(), masses[iN])}; + if (y < yMin + yBeam || y > yMax + yBeam) { + continue; + } + if (nSigmaTPC[iN] < cfgCutsPID->get(iN, 0u) || nSigmaTPC[iN] > cfgCutsPID->get(iN, 1u)) { + continue; + } + if (track.pt() > cfgCutsPID->get(iN, 4u) && (nSigmaTOF[iN] < cfgCutsPID->get(iN, 2u) || nSigmaTOF[iN] > cfgCutsPID->get(iN, 3u))) { + continue; + } + keepEvent[iN] = true; + } + + // + // fill QA histograms + // + spectra.fill(HIST("fTPCsignal"), track.tpcInnerParam(), track.tpcSignal()); + spectra.fill(HIST("fTPCcounts"), track.tpcInnerParam(), track.tpcNSigmaHe()); + + } // end loop over tracks + // + for (int iDecision{0}; iDecision < 4; ++iDecision) { + if (keepEvent[iDecision]) { + spectra.fill(HIST("fProcessedEvents"), iDecision); + } + } + tags(keepEvent[0], keepEvent[1], keepEvent[2], keepEvent[3]); + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfg) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfg)}; +} diff --git a/Analysis/Tasks/CMakeLists.txt b/Analysis/Tasks/CMakeLists.txt index 39b3f394b20ef..5d4dab653051c 100644 --- a/Analysis/Tasks/CMakeLists.txt +++ b/Analysis/Tasks/CMakeLists.txt @@ -18,7 +18,6 @@ add_subdirectory(Utils) add_subdirectory(SkimmingTutorials) add_subdirectory(PID) - o2_add_dpl_workflow(trackextension SOURCES trackextension.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore