diff --git a/.github/workflows/linux-eic-shell.yml b/.github/workflows/linux-eic-shell.yml index 55be36e31f..75e2f9ff21 100644 --- a/.github/workflows/linux-eic-shell.yml +++ b/.github/workflows/linux-eic-shell.yml @@ -1349,6 +1349,92 @@ jobs: if-no-files-found: error include-hidden-files: true + eicrecon-timeframe-splitting: + runs-on: ubuntu-24.04 + needs: + - build + - npsim-dis + strategy: + matrix: + include: + - CXX: clang++ + beam: 18x275 + minq2: 100 + detector_config: craterlake_18x275 + sanitizer: ASAN + + - CXX: clang++ + beam: 18x275 + minq2: 100 + detector_config: craterlake_18x275 + sanitizer: TSAN + nthreads: 4 + + steps: + - name: Checkout .github and scripts + uses: actions/checkout@v7 + with: + sparse-checkout: | + .github + src/scripts + + - name: Download install directory + uses: actions/download-artifact@v8 + with: + name: install-${{ matrix.CXX }}-eic-shell-Release-${{ env.platform }}-${{ env.release }}-${{ matrix.sanitizer }} + + - name: Unarchive install directory + run: tar -xaf install.tar.zst + + - name: Download simulation input + uses: actions/download-artifact@v8 + with: + name: sim_dis_${{ matrix.beam }}_minQ2=${{ matrix.minq2 }}_${{ matrix.detector_config }}.edm4hep.rnt.root + + - name: Setup cvmfs + uses: cvmfs-contrib/github-action-cvmfs@v5 + + - name: Run EICrecon with timeframe splitting + uses: eic/run-cvmfs-osg-eic-shell@main + with: + organization: "${{ env.organization }}" + platform-release: "${{ env.platform }}:${{ env.release }}" + setup: "/opt/detector/epic-${{ env.detector-version }}/bin/thisepic.sh" + run: | + echo "::add-matcher::${{github.workspace}}/.github/ubsan.json" + echo "::add-matcher::${{github.workspace}}/.github/eicrecon.json" + + export DETECTOR_CONFIG=${DETECTOR}_${{ matrix.detector_config }} + export LD_LIBRARY_PATH=$PWD/install/lib:$LD_LIBRARY_PATH + export JANA_PLUGIN_PATH=$PWD/install/lib/EICrecon/plugins:/usr/local/plugins + + prmon --json-summary timeframe_splitter_${{matrix.sanitizer}}.prmon.json -- \ + $PWD/install/bin/eicrecon \ + ${{env.JANA_OPTIONS}} \ + ${{matrix.nthreads > 1 && format('-Pnthreads={0}', matrix.nthreads) || ''}} \ + -Ppodio:output_file=timeframe_splitter_${{matrix.sanitizer}}.edm4eic.root \ + sim_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4hep.rnt.root \ + -Psplit_timeframes=true + + - name: Verify output file + run: | + test -s timeframe_splitter_${{matrix.sanitizer}}.edm4eic.root + + - name: Upload output + uses: actions/upload-artifact@v7 + with: + name: timeframe_splitter_${{matrix.sanitizer}}.edm4eic.root + path: timeframe_splitter_${{matrix.sanitizer}}.edm4eic.root + if-no-files-found: error + + - name: Upload prmon report + uses: actions/upload-artifact@v7 + with: + name: timeframe_splitter_${{matrix.sanitizer}}.prmon.json + path: timeframe_splitter_${{matrix.sanitizer}}.prmon.json + if-no-files-found: error + + compare-single-multi-threaded: runs-on: ubuntu-24.04 needs: diff --git a/src/detectors/B0ECAL/B0ECAL.cc b/src/detectors/B0ECAL/B0ECAL.cc index 2f9ba8a7f6..edb50aec77 100644 --- a/src/detectors/B0ECAL/B0ECAL.cc +++ b/src/detectors/B0ECAL/B0ECAL.cc @@ -2,7 +2,9 @@ // Copyright (C) 2022 - 2025 Whitney Armstrong, Sylvester Joosten, Chao Peng, David Lawrence, Wouter Deconinck, Kolja Kauder, Nathan Brei, Dmitry Kalinkin, Derek Anderson, Michael Pitt #include +#include #include +#include #include #include #include @@ -23,6 +25,9 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; app->Add(new JOmniFactoryGeneratorT( "B0ECalRawHits", {"EventHeader", "B0ECalHits"}, @@ -40,7 +45,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "B0ECalHits", }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "B0ECalRecHits", {"B0ECalRawHits"}, {"B0ECalRecHits"}, { @@ -55,7 +60,7 @@ void InitPlugin(JApplication* app) { .readout = "B0ECalHits", .sectorField = "sector", }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "B0ECalTruthProtoClusters", {"B0ECalRecHits", "B0ECalRawHitLinks"}, {"B0ECalTruthProtoClusters"}, app)); diff --git a/src/detectors/B0TRK/B0TRK.cc b/src/detectors/B0TRK/B0TRK.cc index 9d35952cd4..c8261f7228 100644 --- a/src/detectors/B0TRK/B0TRK.cc +++ b/src/detectors/B0TRK/B0TRK.cc @@ -4,7 +4,9 @@ // #include +#include #include +#include #include #include #include @@ -19,6 +21,9 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Digitization app->Add(new JOmniFactoryGeneratorT( @@ -28,7 +33,7 @@ void InitPlugin(JApplication* app) { .threshold = 10.0 * dd4hep::keV, .timeResolution = 30 * edm4eic::unit::ps, }, - app)); + app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( @@ -36,6 +41,6 @@ void InitPlugin(JApplication* app) { { .timeResolution = 30 * edm4eic::unit::ps, }, - app)); + app, hit_level)); } } // extern "C" diff --git a/src/detectors/BEMC/BEMC.cc b/src/detectors/BEMC/BEMC.cc index e57f7f6c7f..425ad663f6 100644 --- a/src/detectors/BEMC/BEMC.cc +++ b/src/detectors/BEMC/BEMC.cc @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -43,6 +44,9 @@ extern "C" { void InitPlugin(JApplication* app) { using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; InitJANAPlugin(app); @@ -250,7 +254,7 @@ void InitPlugin(JApplication* app) { .readout = "EcalBarrelScFiHits", .fields = {"fiber", "z"}, }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalBarrelScFiRecHits", {"EcalBarrelScFiRawHits"}, {"EcalBarrelScFiRecHits"}, @@ -272,7 +276,7 @@ void InitPlugin(JApplication* app) { .maskPos = "xy", .maskPosFields = {"fiber", "z"}, }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalBarrelScFiProtoClusters", {"EcalBarrelScFiRecHits"}, {"EcalBarrelScFiProtoClusters"}, @@ -373,7 +377,7 @@ void InitPlugin(JApplication* app) { .readout = "EcalBarrelImagingHits", .timeWindow = EcalBarrelImaging_timeWindow, }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalBarrelImagingRawHits", {"EventHeader", "EcalBarrelImagingProcessedHits"}, @@ -390,7 +394,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "EcalBarrelImagingHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalBarrelImagingRecHits", {"EcalBarrelImagingRawHits"}, {"EcalBarrelImagingRecHits"}, @@ -407,7 +411,7 @@ void InitPlugin(JApplication* app) { .layerField = "layer", .sectorField = "sector", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalBarrelImagingProtoClusters", {"EcalBarrelImagingRecHits"}, diff --git a/src/detectors/BHCAL/BHCAL.cc b/src/detectors/BHCAL/BHCAL.cc index 77c983dfb9..4e7f42fb07 100644 --- a/src/detectors/BHCAL/BHCAL.cc +++ b/src/detectors/BHCAL/BHCAL.cc @@ -2,7 +2,9 @@ // Copyright (C) 2022 - 2024 David Lawrence, Derek Anderson, Wouter Deconinck #include +#include #include +#include #include #include #include @@ -27,6 +29,9 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Make sure digi and reco use the same value decltype(CalorimeterHitDigiConfig::capADC) HcalBarrel_capADC = 65536; //65536, 16bit ADC @@ -64,7 +69,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "HcalBarrelHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( @@ -82,7 +87,7 @@ void InitPlugin(JApplication* app) { .layerField = "", .sectorField = "", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); // -------------------------------------------------------------------- diff --git a/src/detectors/BTOF/BTOF.cc b/src/detectors/BTOF/BTOF.cc index 242dc74da0..ddcab94432 100644 --- a/src/detectors/BTOF/BTOF.cc +++ b/src/detectors/BTOF/BTOF.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -35,14 +36,16 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Convert raw digitized hits into calibrated hits // time walk correction is still TBD app->Add(new JOmniFactoryGeneratorT( "TOFBarrelCalibratedHits", {"TOFBarrelADCTDC"}, // Input data collection tags {"TOFBarrelCalibratedHits"}, // Output data tag - {}, - app)); // Hit reco default config for factories + {}, app, hit_level)); // Hit reco default config for factories // cluster all hits in a sensor into one hit location // Currently it's just a simple weighted average @@ -65,7 +68,7 @@ void InitPlugin(JApplication* app) { .min_edep = 6.0 * edm4eic::unit::keV, .readout = "TOFBarrelHits", }, - app)); + app, hit_level)); // temporary steps to bypass pulse digitization and jump right from ChargeSharing to clusters // Avoid efficiency loss until we can simulate hardware accurately @@ -76,14 +79,13 @@ void InitPlugin(JApplication* app) { .threshold = 0.0, .timeResolution = 0.025, // [ns] }, - app)); + app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( "TOFBarrelSharedRecHits", {"TOFBarrelSharedRawHits"}, // Input data collection tags {"TOFBarrelSharedRecHits"}, // Output data tag - {}, - app)); // Hit reco default config for factories + {}, app, hit_level)); // Hit reco default config for factories // calculation of the extreme values for Landau distribution can be found on lin 514-520 of // https://root.cern.ch/root/html524/src/TMath.cxx.html#fsokrB Landau reaches minimum for mpv = diff --git a/src/detectors/BTRK/BTRK.cc b/src/detectors/BTRK/BTRK.cc index 2b0497daad..185f14f307 100644 --- a/src/detectors/BTRK/BTRK.cc +++ b/src/detectors/BTRK/BTRK.cc @@ -4,7 +4,9 @@ // #include +#include #include +#include #include #include #include @@ -22,6 +24,9 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Digitization app->Add(new JOmniFactoryGeneratorT( @@ -30,24 +35,24 @@ void InitPlugin(JApplication* app) { { .threshold = 0.54 * dd4hep::keV, }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "SiBarrelNoiseRawHits", {"EventHeader"}, {"SiBarrelNoiseRawHits"}, {.addNoise = true, .noise_rate_per_pixel_per_event = 2.0e-7, .readout_name = "SiBarrelHits"}, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT>( "SiBarrelRawHitsWithNoise", // Name of the combiner instance {"SiBarrelRawHits", "SiBarrelNoiseRawHits"}, // Inputs: original + noise-only {"SiBarrelRawHitsWithNoise"}, // Output: merged collection {}, // default config - app)); + app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( "SiBarrelTrackerRecHits", {"SiBarrelRawHitsWithNoise"}, {"SiBarrelTrackerRecHits"}, {}, // default config - app)); + app, hit_level)); } } // extern "C" diff --git a/src/detectors/BVTX/BVTX.cc b/src/detectors/BVTX/BVTX.cc index c15a0210cf..4db57f7aad 100644 --- a/src/detectors/BVTX/BVTX.cc +++ b/src/detectors/BVTX/BVTX.cc @@ -5,7 +5,9 @@ // #include +#include #include +#include #include #include #include @@ -23,6 +25,9 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Digitization app->Add(new JOmniFactoryGeneratorT( @@ -31,21 +36,21 @@ void InitPlugin(JApplication* app) { { .threshold = 0.54 * dd4hep::keV, }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "SiBarrelVertexNoiseRawHits", {"EventHeader"}, {"SiBarrelVertexNoiseRawHits"}, {.addNoise = true, .noise_rate_per_pixel_per_event = 2.0e-7, .readout_name = "VertexBarrelHits"}, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT>( "SiBarrelVertexRawHitsWithNoise", {"SiBarrelVertexRawHits", "SiBarrelVertexNoiseRawHits"}, - {"SiBarrelVertexRawHitsWithNoise"}, {}, app)); + {"SiBarrelVertexRawHitsWithNoise"}, {}, app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( "SiBarrelVertexRecHits", {"SiBarrelVertexRawHitsWithNoise"}, {"SiBarrelVertexRecHits"}, {}, // default config - app)); + app, hit_level)); } } // extern "C" diff --git a/src/detectors/ECTOF/ECTOF.cc b/src/detectors/ECTOF/ECTOF.cc index 9bff9a0459..e17c34f30a 100644 --- a/src/detectors/ECTOF/ECTOF.cc +++ b/src/detectors/ECTOF/ECTOF.cc @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -30,6 +31,9 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // cluster all hits in a sensor into one hit location // Currently it's just a simple weighted average @@ -53,7 +57,7 @@ void InitPlugin(JApplication* app) { .min_edep = 6 * dd4hep::keV, .readout = "TOFEndcapHits", }, - app)); + app, hit_level)); // temporary steps to bypass pulse digitization and jump right from ChargeSharing to clusters // Avoid efficiency loss until we can simulate hardware accurately @@ -64,14 +68,13 @@ void InitPlugin(JApplication* app) { .threshold = 0.0, .timeResolution = 0.025, // [ns] }, - app)); + app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( "TOFEndcapSharedRecHits", {"TOFEndcapSharedRawHits"}, // Input data collection tags {"TOFEndcapSharedRecHits"}, // Output data tag - {}, - app)); // Hit reco default config for factories + {}, app, hit_level)); // Hit reco default config for factories const double x_when_landau_min = -0.22278; const double landau_min = TMath::Landau(x_when_landau_min, 0, 1, true); diff --git a/src/detectors/ECTRK/ECTRK.cc b/src/detectors/ECTRK/ECTRK.cc index b8eda7d508..11bd070fea 100644 --- a/src/detectors/ECTRK/ECTRK.cc +++ b/src/detectors/ECTRK/ECTRK.cc @@ -4,7 +4,9 @@ // #include +#include #include +#include #include #include #include @@ -22,6 +24,9 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Digitization app->Add(new JOmniFactoryGeneratorT( @@ -30,21 +35,21 @@ void InitPlugin(JApplication* app) { { .threshold = 0.54 * dd4hep::keV, }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "SiEndcapTrackerNoiseRawHits", {"EventHeader"}, {"SiEndcapTrackerNoiseRawHits"}, {.addNoise = true, .noise_rate_per_pixel_per_event = 2.0e-7, .readout_name = "TrackerEndcapHits"}, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT>( "SiEndcapTrackerRawHitsWithNoise", {"SiEndcapTrackerRawHits", "SiEndcapTrackerNoiseRawHits"}, - {"SiEndcapTrackerRawHitsWithNoise"}, {}, app)); + {"SiEndcapTrackerRawHitsWithNoise"}, {}, app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( "SiEndcapTrackerRecHits", {"SiEndcapTrackerRawHitsWithNoise"}, {"SiEndcapTrackerRecHits"}, {}, // default config - app)); + app, hit_level)); } } // extern "C" diff --git a/src/detectors/EEMC/EEMC.cc b/src/detectors/EEMC/EEMC.cc index 930b36ce71..3a9747fc38 100644 --- a/src/detectors/EEMC/EEMC.cc +++ b/src/detectors/EEMC/EEMC.cc @@ -2,7 +2,9 @@ // Copyright (C) 2022 - 2025 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong, Thomas Britton, David Lawrence, Dhevan Gangadharan, Wouter Deconinck, Dmitry Kalinkin, Derek Anderson #include +#include #include +#include #include #include #include @@ -26,6 +28,9 @@ extern "C" { void InitPlugin(JApplication* app) { using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; InitJANAPlugin(app); @@ -59,7 +64,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "EcalEndcapNHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalEndcapNRecHits", {"EcalEndcapNRawHits"}, {"EcalEndcapNRecHits"}, @@ -74,7 +79,7 @@ void InitPlugin(JApplication* app) { .sampFrac = "0.96", .readout = "EcalEndcapNHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalEndcapNTruthProtoClusters", {"EcalEndcapNRecHits", "EcalEndcapNRawHitLinks"}, diff --git a/src/detectors/EHCAL/EHCAL.cc b/src/detectors/EHCAL/EHCAL.cc index 183d5eefcd..6a8d346200 100644 --- a/src/detectors/EHCAL/EHCAL.cc +++ b/src/detectors/EHCAL/EHCAL.cc @@ -2,7 +2,9 @@ // Copyright (C) 2022 - 2025 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong, David Lawrence, Friederike Bock, Nathan Brei, Wouter Deconinck, Dmitry Kalinkin, Derek Anderson #include +#include #include +#include #include #include #include @@ -25,6 +27,9 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Make sure digi and reco use the same value decltype(CalorimeterHitDigiConfig::capADC) HcalEndcapN_capADC = 32768; // assuming 15 bit ADC like FHCal @@ -50,7 +55,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "HcalEndcapNHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "HcalEndcapNRecHits", {"HcalEndcapNRawHits"}, {"HcalEndcapNRecHits"}, @@ -67,7 +72,7 @@ void InitPlugin(JApplication* app) { "0.0095", // from latest study - implement at level of reco hits rather than clusters .readout = "HcalEndcapNHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "HcalEndcapNMergedHits", {"HcalEndcapNRecHits"}, {"HcalEndcapNMergedHits"}, diff --git a/src/detectors/FEMC/FEMC.cc b/src/detectors/FEMC/FEMC.cc index 1f19cca3f1..e2277de632 100644 --- a/src/detectors/FEMC/FEMC.cc +++ b/src/detectors/FEMC/FEMC.cc @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -33,6 +34,9 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; auto log_service = app->GetService(); auto mLog = log_service->logger("FEMC"); @@ -89,7 +93,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "EcalEndcapPHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); } else if (EcalEndcapP_homogeneousFlag == 2) { app->Add(new JOmniFactoryGeneratorT( @@ -115,7 +119,7 @@ void InitPlugin(JApplication* app) { .readout = "EcalEndcapPHits", .fields = {"fiber_x", "fiber_y"}, }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); } @@ -134,7 +138,7 @@ void InitPlugin(JApplication* app) { .sampFrac = "1.00", // already taken care in DIGI code above .readout = "EcalEndcapPHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalEndcapPTruthProtoClusters", {"EcalEndcapPRecHits", "EcalEndcapPRawHitLinks"}, diff --git a/src/detectors/FHCAL/FHCAL.cc b/src/detectors/FHCAL/FHCAL.cc index 56f5baff0c..1dc64ccb45 100644 --- a/src/detectors/FHCAL/FHCAL.cc +++ b/src/detectors/FHCAL/FHCAL.cc @@ -2,7 +2,9 @@ // Copyright (C) 2023 - 2025 Friederike Bock, Wouter Deconinck #include +#include #include +#include #include #include #include @@ -29,6 +31,9 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // Make sure digi and reco use the same value decltype(CalorimeterHitDigiConfig::capADC) HcalEndcapPInsert_capADC = 32768; @@ -53,7 +58,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "HcalEndcapPInsertHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "HcalEndcapPInsertRecHits", {"HcalEndcapPInsertRawHits"}, {"HcalEndcapPInsertRecHits"}, @@ -70,7 +75,7 @@ void InitPlugin(JApplication* app) { .readout = "HcalEndcapPInsertHits", .layerField = "layer", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "HcalEndcapPInsertMergedHits", {"HcalEndcapPInsertRecHits"}, {"HcalEndcapPInsertMergedHits"}, @@ -193,7 +198,7 @@ void InitPlugin(JApplication* app) { .readout = "LFHCALHits", .fields = {"layerz"}, }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "LFHCALRecHits", {"LFHCALRawHits"}, {"LFHCALRecHits"}, @@ -209,7 +214,7 @@ void InitPlugin(JApplication* app) { .readout = "LFHCALHits", .layerField = "rlayerz", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "LFHCALTruthProtoClusters", {"LFHCALRecHits", "LFHCALRawHitLinks"}, diff --git a/src/detectors/FOFFMTRK/FOFFMTRK.cc b/src/detectors/FOFFMTRK/FOFFMTRK.cc index 2bbf035a86..29a28e3644 100644 --- a/src/detectors/FOFFMTRK/FOFFMTRK.cc +++ b/src/detectors/FOFFMTRK/FOFFMTRK.cc @@ -4,7 +4,9 @@ // #include +#include #include +#include #include #include #include @@ -20,6 +22,9 @@ extern "C" { void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; //Digitized hits, especially for thresholds app->Add(new JOmniFactoryGeneratorT( @@ -30,14 +35,14 @@ void InitPlugin(JApplication* app) { .threshold = 10.0 * dd4hep::keV, .timeResolution = 30 * edm4eic::unit::ps, }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "ForwardOffMTrackerRecHits", {"ForwardOffMTrackerRawHits"}, {"ForwardOffMTrackerRecHits"}, { .timeResolution = 30 * edm4eic::unit::ps, }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "ForwardOffMRecParticles", {"MCParticles", "ForwardOffMTrackerRecHits"}, diff --git a/src/detectors/LOWQ2/LOWQ2.cc b/src/detectors/LOWQ2/LOWQ2.cc index 1a82cca322..65d7d385cd 100644 --- a/src/detectors/LOWQ2/LOWQ2.cc +++ b/src/detectors/LOWQ2/LOWQ2.cc @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -45,6 +46,9 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; std::string readout = "TaggerTrackerHits"; @@ -96,7 +100,7 @@ void InitPlugin(JApplication* app) { .threshold = 1.5 * edm4eic::unit::keV, .timeResolution = 2 * edm4eic::unit::ns, }, - app)); + app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( @@ -104,7 +108,7 @@ void InitPlugin(JApplication* app) { { .timeResolution = 2, }, - app)); + app, hit_level)); // Divide collection based on geometry segmentation labels // This should really be done before digitization as summing hits in the same cell couldn't even be mixed between layers. At the moment just prep for clustering. diff --git a/src/detectors/LUMISPECCAL/LUMISPECCAL.cc b/src/detectors/LUMISPECCAL/LUMISPECCAL.cc index 03e799716d..196b041b12 100644 --- a/src/detectors/LUMISPECCAL/LUMISPECCAL.cc +++ b/src/detectors/LUMISPECCAL/LUMISPECCAL.cc @@ -2,7 +2,9 @@ // Copyright (C) 2022 - 2025 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong, David Lawrence, Dhevan Gangadharan, Nathan Brei,, Wouter Deconinck, Dmitry Kalinkin, Derek Anderson #include +#include #include +#include #include #include #include @@ -23,6 +25,9 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; app->Add(new JOmniFactoryGeneratorT( "EcalLumiSpecRawHits", {"EventHeader", "EcalLumiSpecHits"}, @@ -38,7 +43,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "EcalLumiSpecHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalLumiSpecRecHits", {"EcalLumiSpecRawHits"}, {"EcalLumiSpecRecHits"}, @@ -53,7 +58,7 @@ void InitPlugin(JApplication* app) { .sampFrac = "1.0", .readout = "EcalLumiSpecHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalLumiSpecTruthProtoClusters", {"EcalLumiSpecRecHits", "EcalLumiSpecRawHitLinks"}, diff --git a/src/detectors/MPGD/MPGD.cc b/src/detectors/MPGD/MPGD.cc index 4c2879dcca..3c4a6e2c97 100644 --- a/src/detectors/MPGD/MPGD.cc +++ b/src/detectors/MPGD/MPGD.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +41,9 @@ void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // ***** PIXEL or 2DSTRIP DIGITIZATION? // - This determines which of the MPGDTrackerDigi or SiliconTrackerDigi @@ -106,7 +110,7 @@ void InitPlugin(JApplication* app) { .threshold = 100 * dd4hep::eV, .timeResolution = 10, }, - app)); + app, hit_level)); } else { // Configuration parameters MPGDTrackerDigiConfig digi_cfg; @@ -131,7 +135,7 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "MPGDBarrelRawHits", {"EventHeader", "MPGDBarrelHits"}, {"MPGDBarrelRawHits", "MPGDBarrelRawHitLinks", "MPGDBarrelRawHitAssociations"}, digi_cfg, - app)); + app, hit_level)); } // Convert raw digitized hits into hits with geometry info (ready for tracking) @@ -142,7 +146,7 @@ void InitPlugin(JApplication* app) { { .timeResolution = 10, }, - app)); + app, hit_level)); } else { MPGDHitReconstructionConfig reco_cfg; reco_cfg.readout = "MPGDBarrelHits"; @@ -151,7 +155,7 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "MPGDBarrelRecHits", {"MPGDBarrelRawHits"}, // Input data collection tags {"MPGDBarrelRecHits"}, // Output data tag - reco_cfg, app)); + reco_cfg, app, hit_level)); } // ***** OuterMPGDBarrel @@ -165,7 +169,7 @@ void InitPlugin(JApplication* app) { .threshold = 100 * dd4hep::eV, .timeResolution = 10, }, - app)); + app, hit_level)); } else { MPGDTrackerDigiConfig digi_cfg; digi_cfg.readout = "OuterMPGDBarrelHits"; @@ -186,7 +190,7 @@ void InitPlugin(JApplication* app) { "OuterMPGDBarrelRawHits", {"EventHeader", "OuterMPGDBarrelHits"}, {"OuterMPGDBarrelRawHits", "OuterMPGDBarrelRawHitLinks", "OuterMPGDBarrelRawHitAssociations"}, - digi_cfg, app)); + digi_cfg, app, hit_level)); } // Convert raw digitized hits into hits with geometry info (ready for tracking) @@ -197,7 +201,7 @@ void InitPlugin(JApplication* app) { { .timeResolution = 10, }, - app)); + app, hit_level)); } else { MPGDHitReconstructionConfig reco_cfg; reco_cfg.readout = "OuterMPGDBarrelHits"; @@ -206,7 +210,7 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "OuterMPGDBarrelRecHits", {"OuterMPGDBarrelRawHits"}, // Input data collection tags {"OuterMPGDBarrelRecHits"}, // Output data tag - reco_cfg, app)); + reco_cfg, app, hit_level)); } // ***** "BackwardMPGDEndcap" @@ -219,7 +223,7 @@ void InitPlugin(JApplication* app) { .threshold = 100 * dd4hep::eV, .timeResolution = 10, }, - app)); + app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( @@ -228,7 +232,7 @@ void InitPlugin(JApplication* app) { { .timeResolution = 10, }, - app)); + app, hit_level)); // ""ForwardMPGDEndcap" // Digitization @@ -240,7 +244,7 @@ void InitPlugin(JApplication* app) { .threshold = 100 * dd4hep::eV, .timeResolution = 10, }, - app)); + app, hit_level)); // Convert raw digitized hits into hits with geometry info (ready for tracking) app->Add(new JOmniFactoryGeneratorT( @@ -249,6 +253,6 @@ void InitPlugin(JApplication* app) { { .timeResolution = 10, }, - app)); + app, hit_level)); } } // extern "C" diff --git a/src/detectors/RPOTS/RPOTS.cc b/src/detectors/RPOTS/RPOTS.cc index e0cbb06f1b..76e5b69a2f 100644 --- a/src/detectors/RPOTS/RPOTS.cc +++ b/src/detectors/RPOTS/RPOTS.cc @@ -4,7 +4,9 @@ // #include +#include #include +#include #include #include #include @@ -22,6 +24,9 @@ extern "C" { void InitPlugin(JApplication* app) { InitJANAPlugin(app); using namespace eicrecon; + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; MatrixTransferStaticConfig recon_cfg; PolynomialMatrixReconstructionConfig recon_poly_cfg; @@ -34,14 +39,14 @@ void InitPlugin(JApplication* app) { .threshold = 10.0 * dd4hep::keV, .timeResolution = 30 * edm4eic::unit::ps, }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "ForwardRomanPotRecHits", {"ForwardRomanPotRawHits"}, {"ForwardRomanPotRecHits"}, { .timeResolution = 30 * edm4eic::unit::ps, }, - app)); + app, hit_level)); app->Add(new JOmniFactoryGeneratorT( "ForwardRomanPotStaticRecParticles", diff --git a/src/detectors/ZDC/ZDC.cc b/src/detectors/ZDC/ZDC.cc index 11e8f1344f..b9d18f7929 100644 --- a/src/detectors/ZDC/ZDC.cc +++ b/src/detectors/ZDC/ZDC.cc @@ -2,7 +2,9 @@ // Copyright (C) 2021 - 2025, Chao Peng, Sylvester Joosten, Whitney Armstrong, David Lawrence, Friederike Bock, Wouter Deconinck, Nathan Brei, Sebouh Paul, Dmitry Kalinkin, Barak Schmookler #include +#include #include +#include #include #include #include @@ -26,6 +28,9 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto hit_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; // LYSO part of the ZDC app->Add(new JOmniFactoryGeneratorT( @@ -43,7 +48,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "EcalFarForwardZDCHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalFarForwardZDCRecHits", {"EcalFarForwardZDCRawHits"}, {"EcalFarForwardZDCRecHits"}, @@ -58,7 +63,7 @@ void InitPlugin(JApplication* app) { .sampFrac = "1.0", .readout = "EcalFarForwardZDCHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( "EcalFarForwardZDCTruthProtoClusters", @@ -152,7 +157,7 @@ void InitPlugin(JApplication* app) { .corrMeanScale = "1.0", .readout = "HcalFarForwardZDCHits", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( @@ -170,7 +175,7 @@ void InitPlugin(JApplication* app) { .layerField = "layer", .sectorField = "system", }, - app // TODO: Remove me once fixed + app, hit_level // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( diff --git a/src/extensions/jana/JOmniFactoryGeneratorT.h b/src/extensions/jana/JOmniFactoryGeneratorT.h index f5846bf4fc..1326b3b51e 100644 --- a/src/extensions/jana/JOmniFactoryGeneratorT.h +++ b/src/extensions/jana/JOmniFactoryGeneratorT.h @@ -18,6 +18,7 @@ template class JOmniFactoryGeneratorT : public JFactoryGenerato std::vector m_default_input_tags; std::vector m_default_output_tags; FactoryConfigType m_default_cfg; /// Must be properly copyable! + JEventLevel m_level = JEventLevel::PhysicsEvent; }; struct UntypedWiring { @@ -30,21 +31,25 @@ template class JOmniFactoryGeneratorT : public JFactoryGenerato public: explicit JOmniFactoryGeneratorT(std::string tag, std::vector default_input_tags, std::vector default_output_tags, - FactoryConfigType cfg, JApplication* app) { + FactoryConfigType cfg, JApplication* app, + JEventLevel level = JEventLevel::PhysicsEvent) { m_app = app; m_wirings.push_back({.m_tag = tag, .m_default_input_tags = default_input_tags, .m_default_output_tags = default_output_tags, - .m_default_cfg = cfg}); + .m_default_cfg = cfg, + .m_level = level}); }; explicit JOmniFactoryGeneratorT(std::string tag, std::vector default_input_tags, - std::vector default_output_tags, JApplication* app) { + std::vector default_output_tags, JApplication* app, + JEventLevel level = JEventLevel::PhysicsEvent) { m_app = app; m_wirings.push_back({.m_tag = tag, .m_default_input_tags = default_input_tags, .m_default_output_tags = default_output_tags, - .m_default_cfg = {}}); + .m_default_cfg = {}, + .m_level = level}); } explicit JOmniFactoryGeneratorT(JApplication* app) : m_app(app) {} @@ -55,7 +60,8 @@ template class JOmniFactoryGeneratorT : public JFactoryGenerato m_wirings.push_back({.m_tag = tag, .m_default_input_tags = default_input_tags, .m_default_output_tags = default_output_tags, - .m_default_cfg = cfg}); + .m_default_cfg = cfg, + .m_level = JEventLevel::PhysicsEvent}); } void AddWiring(std::string tag, std::vector default_input_tags, @@ -70,7 +76,8 @@ template class JOmniFactoryGeneratorT : public JFactoryGenerato m_wirings.push_back({.m_tag = tag, .m_default_input_tags = default_input_tags, .m_default_output_tags = default_output_tags, - .m_default_cfg = config}); + .m_default_cfg = config, + .m_level = JEventLevel::PhysicsEvent}); } void GenerateFactories(JFactorySet* factory_set) override { @@ -82,6 +89,7 @@ template class JOmniFactoryGeneratorT : public JFactoryGenerato factory->SetPluginName(this->GetPluginName()); factory->SetFactoryName(JTypeInfo::demangle()); factory->config() = wiring.m_default_cfg; + factory->SetLevel(wiring.m_level); // Set up all of the wiring prereqs so that Init() can do its thing // Specifically, it needs valid input/output tags, a valid logger, and diff --git a/src/extensions/jana/JOmniUnfolderGeneratorT.h b/src/extensions/jana/JOmniUnfolderGeneratorT.h new file mode 100644 index 0000000000..b35f313590 --- /dev/null +++ b/src/extensions/jana/JOmniUnfolderGeneratorT.h @@ -0,0 +1,71 @@ +// Copyright 2023, Jefferson Science Associates, LLC. +// Subject to the terms in the LICENSE file found in the top-level directory. +// Created by Nathan Brei + +#pragma once + +#include +#include +#include +#include + +namespace eicrecon { + +template class JOmniUnfolderGeneratorT { +public: + using ConfigT = ComponentT::ConfigT; + + struct Wiring { + std::string tag; + JEventLevel parent_level = JEventLevel::Timeslice; + JEventLevel child_level = JEventLevel::PhysicsEvent; + std::vector input_names; + std::vector input_levels; + std::vector> variadic_input_names; + std::vector variadic_input_levels; + std::vector output_names; + std::vector> variadic_output_names; + ConfigT configs = {}; // Must be copyable! + }; + +private: + Wiring m_wiring; + +public: + explicit JOmniUnfolderGeneratorT(Wiring&& wiring) : m_wiring(wiring) {} + + void Generate(JApplication* app) { + + auto component = new ComponentT; + + // EICrecon JOFG does _not_ include plugin name in prefix + // Unlike with factories, plugin name is set automatically during app->Add() + component->SetPrefix(m_wiring.tag); + component->SetTypeName(JTypeInfo::demangle()); + component->SetParentLevel(m_wiring.parent_level); + component->SetChildLevel(m_wiring.child_level); + component->config() = m_wiring.configs; + + // Override the input and output collection names via parameters + // (This is also EICrecon-specific, the long term plan is to use the wiring file instead) + app->SetDefaultParameter(m_wiring.tag + ":InputTags", m_wiring.input_names, + "Input collection names"); + app->SetDefaultParameter(m_wiring.tag + ":OutputTags", m_wiring.output_names, + "Output collection names"); + app->SetDefaultParameter(m_wiring.tag + ":VariadicInputTags", m_wiring.variadic_input_names, + "Input collection names"); + app->SetDefaultParameter(m_wiring.tag + ":VariadicOutputTags", m_wiring.variadic_output_names, + "Output collection names"); + + component->WireInputs(m_wiring.parent_level, m_wiring.input_levels, m_wiring.input_names, + m_wiring.variadic_input_levels, m_wiring.variadic_input_names); + component->WireOutputs(m_wiring.child_level, m_wiring.output_names, + m_wiring.variadic_output_names, + false // use_short_names + ); + + app->Add(component); + } +}; + +} // namespace eicrecon diff --git a/src/global/CMakeLists.txt b/src/global/CMakeLists.txt index bdca9fd35e..81dc003eca 100644 --- a/src/global/CMakeLists.txt +++ b/src/global/CMakeLists.txt @@ -4,3 +4,6 @@ add_subdirectory(particle_flow) add_subdirectory(pid) add_subdirectory(pid_lut) add_subdirectory(beam) +if(JANA_VERSION VERSION_GREATER_EQUAL "2026.01.00") + add_subdirectory(splitting) +endif() diff --git a/src/global/splitting/CMakeLists.txt b/src/global/splitting/CMakeLists.txt new file mode 100644 index 0000000000..98423cf00d --- /dev/null +++ b/src/global/splitting/CMakeLists.txt @@ -0,0 +1,19 @@ +get_filename_component(PLUGIN_NAME ${CMAKE_CURRENT_LIST_DIR} NAME) + +# Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets +# Setting default includes, libraries and installation paths +plugin_add(${PLUGIN_NAME} PLUGIN_USE_CC_ONLY) + +# The macro grabs sources as *.cc *.cpp *.c and headers as *.h *.hh *.hpp Then +# correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds +# headers to the correct installation directory +plugin_glob_all(${PLUGIN_NAME}) + +# Find dependencies +plugin_add_cern_root(${PLUGIN_NAME}) +plugin_add_event_model(${PLUGIN_NAME}) + +# Add libraries (same as target_include_directories but for both plugin and +# library) +plugin_link_libraries(${PLUGIN_NAME} algorithms_digi_library + algorithms_tracking_library) diff --git a/src/global/splitting/TimeframeSplitter.cc b/src/global/splitting/TimeframeSplitter.cc new file mode 100644 index 0000000000..98c0e0036d --- /dev/null +++ b/src/global/splitting/TimeframeSplitter.cc @@ -0,0 +1,921 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2026 Takuya Kumaoka + +#include "TimeframeSplitter.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +TimeframeSplitter::TimeframeSplitter() { + SetTypeName(NAME_OF_THIS); + SetParentLevel(JEventLevel::Timeslice); + SetChildLevel(JEventLevel::PhysicsEvent); +} + +std::uint64_t TimeframeSplitter::objIdKey(const podio::ObjectID& objectID) { + const auto collectionID = + static_cast(static_cast(objectID.collectionID)); + const auto index = static_cast(static_cast(objectID.index)); + + return (collectionID << 32U) | index; +} + +TimeframeSplitter::TrackerAssociationIndex TimeframeSplitter::buildTrkAssoId( + const edm4eic::MCRecoTrackerHitAssociationCollection* associations) { + + TrackerAssociationIndex index; + if (associations == nullptr) { + return index; + } + index.reserve(associations->size()); + + for (size_t assoId = 0; assoId < associations->size(); ++assoId) { + const auto association = associations->at(assoId); + const auto rawHit = association.getRawHit(); + if (!rawHit.isAvailable()) { + continue; + } + index[objIdKey(rawHit.getObjectID())].push_back(assoId); + } + + return index; +} + +TimeframeSplitter::CalorimeterAssociationIndex TimeframeSplitter::buildCalAssoId( + const edm4eic::MCRecoCalorimeterHitAssociationCollection* associations) { + + CalorimeterAssociationIndex index; + if (associations == nullptr) { + return index; + } + index.reserve(associations->size()); + + for (size_t assoId = 0; assoId < associations->size(); ++assoId) { + const auto association = associations->at(assoId); + const auto rawHit = association.getRawHit(); + if (!rawHit.isAvailable()) { + continue; + } + index[objIdKey(rawHit.getObjectID())].push_back(assoId); + } + + return index; +} + +bool TimeframeSplitter::overlapsTimeWindow(double hitTime, double resolution, double window_start, + double window_end) { + return hitTime + resolution > window_start && hitTime - resolution < window_end; +} + +bool TimeframeSplitter::judgeOverTimeWindow(double hitTime, double resolution, double window_end) { + return hitTime - resolution >= window_end; +} + +bool TimeframeSplitter::isValidEtaPhiBin(int etaBin, int phiBin) { + return 0 <= etaBin && etaBin < kEtaPhiBins && 0 <= phiBin && phiBin < kEtaPhiBins; +} + +bool TimeframeSplitter::judgeHitInTimeSlice(double hitTime, double timeResolution, + double timeslice_start, double timeslice_end) { + return !(hitTime + timeResolution < timeslice_start || hitTime - timeResolution > timeslice_end); +} + +std::pair TimeframeSplitter::etaPhiBins(double hitEta, double hitPhi, double etaMin, + double etaMax, int bShift) { + const double etaBinWidth = (etaMax - etaMin) / kEtaPhiBins; + const double phiMin = -std::numbers::pi; + const double phiMax = std::numbers::pi; + const double phiBinWidth = (phiMax - phiMin) / kEtaPhiBins; + + const double halfEtaBin = 0.5 * etaBinWidth * bShift; + const double halfPhiBin = 0.5 * phiBinWidth * bShift; + const double shiftedEtaMin = etaMin + halfEtaBin; + const double shiftedEtaMax = etaMax + halfEtaBin; + const double shiftedPhiMin = phiMin + halfPhiBin; + const double shiftedPhiMax = phiMax + halfPhiBin; + + if (hitEta < shiftedEtaMin || hitEta >= shiftedEtaMax || hitPhi < shiftedPhiMin || + hitPhi >= shiftedPhiMax) { + return {kInvalidEtaPhiBin, kInvalidEtaPhiBin}; + } + + const int etaBin = static_cast(std::floor((hitEta - shiftedEtaMin) / etaBinWidth)); + const int phiBin = static_cast(std::floor((hitPhi - shiftedPhiMin) / phiBinWidth)); + return {etaBin, phiBin}; +} + +size_t TimeframeSplitter::countGridCellsWithMultiplicity(const EtaPhiGrid& grid0, + const EtaPhiGrid& gridShifted, + const EtaPhiTimeGrid& gridTime0, + const EtaPhiTimeGrid& gridShiftedTime, + int threshold, double& averageTime) { + size_t count = 0; + double timeSum = 0.0; + for (size_t iEta = 0; iEta < kEtaPhiBins; ++iEta) { + for (size_t iPhi = 0; iPhi < kEtaPhiBins; ++iPhi) { + if (grid0[iEta][iPhi] >= threshold || gridShifted[iEta][iPhi] >= threshold) { + count++; + if (grid0[iEta][iPhi] >= threshold) { + timeSum += gridTime0[iEta][iPhi] / grid0[iEta][iPhi]; + } else if (gridShifted[iEta][iPhi] >= threshold) { + timeSum += gridShiftedTime[iEta][iPhi] / gridShifted[iEta][iPhi]; + } + } + } + } + averageTime = count > 0 ? timeSum / count : 0.0; + return count; +} + +double TimeframeSplitter::averageSelectedTriggerTime(const std::array& values, + const std::array& times, + std::initializer_list indices, + double fallbackTime) { + size_t count = 0; + double timeSum = 0.0; + for (const size_t index : indices) { + if (values[index] <= 0.0) { + continue; + } + timeSum += times[index]; + count++; + } + return count > 0 ? timeSum / count : fallbackTime; +} + +std::pair TimeframeSplitter::backEndEtaPhiBins(double hitEta, double hitPhi, int bShift) { + return etaPhiBins(hitEta, hitPhi, m_backwardEtaMin(), m_backwardEtaMax(), bShift); +} + +std::pair TimeframeSplitter::barrelEtaPhiBins(double hitEta, double hitPhi, int bShift) { + return etaPhiBins(hitEta, hitPhi, m_barrelEtaMin(), m_barrelEtaMax(), bShift); +} + +std::pair TimeframeSplitter::forwardEndEtaPhiBins(double hitEta, double hitPhi, + int bShift) { + return etaPhiBins(hitEta, hitPhi, m_forwardEtaMin(), m_forwardEtaMax(), bShift); +} + +double TimeframeSplitter::trkTimeResolution(TrkCollectionIndex detectorID) { + switch (detectorID) { + case kTrkB0: + [[fallthrough]]; + case kTrkTOFBarrel: + [[fallthrough]]; + case kTrkTOFEndcap: + return m_timeResolution_ACLGad(); + + case kTrkMPGDBarrel: + [[fallthrough]]; + case kTrkOuterMPGDBarrel: + [[fallthrough]]; + case kTrkBackwardMPGD: + [[fallthrough]]; + case kTrkForwardMPGD: + return m_timeResolution_MPGD(); + + case kTrkSiBarrelVertex: + [[fallthrough]]; + case kTrkSiBarrel: + [[fallthrough]]; + case kTrkSiEndcap: + [[fallthrough]]; + case kTrkTagger: + [[fallthrough]]; + case kTrkForwardRomanPot: + [[fallthrough]]; + case kTrkForwardOffMTracker: + return m_timeResolution_SiMaps(); + + case kTrkCollectionSize: + break; + } + + throw std::runtime_error("Unknown tracker detector ID"); +} + +double TimeframeSplitter::calTimeResolution(CalCollectionIndex detectorID) { + switch (detectorID) { + case kCalB0ECal: + [[fallthrough]]; + case kCalEcalBarrelImg: + [[fallthrough]]; + case kCalEcalBarrelScFi: + [[fallthrough]]; + case kCalEcalEndcapN: + [[fallthrough]]; + case kCalEcalEndcapP: + [[fallthrough]]; + case kCalEcalZDC: + [[fallthrough]]; + case kCalEcalLumiSpec: + return m_timeResolution_EMCal(); + + case kCalHcalBarrel: + [[fallthrough]]; + case kCalHcalEndcapN: + [[fallthrough]]; + case kCalHcalEndcapPInsert: + [[fallthrough]]; + case kCalHcalZDC: + [[fallthrough]]; + case kCalLFHCAL: + return m_timeResolution_HCal(); // TODO: use dedicated HCal resolution if needed + + case kCalCollectionSize: + break; + } + + throw std::runtime_error("Unknown calorimeter detector ID"); +} + +TimeframeSplitter::Result TimeframeSplitter::Unfold(const JEvent& parent, JEvent& child, + int child_idx) { + const float timeframeWidthValue = m_timeframeWidth(); + const float timesplitWidthValue = m_timesplitWidth(); + + bool bTrigger = false; + + // Materialize these output collections even when no timeslice triggers, + // so the output schema is available from the first event. + (void)m_eventHeader_outCol(); + (void)m_eventHeaderPhy_outCols(); + (void)m_eventHeaderBkg_outCols(); + + const auto trackerHitCollsIn = m_trackerHits_inCols(); + const auto caloRecHitCollsIn = m_calorimeterHit_inCols(); + const auto trkAssoCollsIn = m_trackerHitsAsso_inCols(); + const auto richRawHitCollsIn = m_richRawHits_inCols(); + const auto richAssoCollsIn = m_richHitsAsso_inCols(); + const auto calrecAssoCollsIn = m_mcRecoCalorimeterHitAssociation_inCols(); + + // == s == Register hits of TOF and MPGD detectors in the time slice ================== + if (child_idx == 0) { + // Association collections belong to the parent Timeslice and remain valid for all + // of its children. Build each lookup once per parent instead of rebuilding the full + // index for every triggered PhysicsEvent. + m_trkAssoIds.clear(); + m_trkAssoIds.reserve(trkAssoCollsIn.size()); + for (const auto* associations : trkAssoCollsIn) { + m_trkAssoIds.push_back(buildTrkAssoId(associations)); + } + + m_richAssoIds.clear(); + m_richAssoIds.reserve(richAssoCollsIn.size()); + for (const auto* associations : richAssoCollsIn) { + m_richAssoIds.push_back(buildTrkAssoId(associations)); + } + + m_calAssoIds.clear(); + m_calAssoIds.reserve(calrecAssoCollsIn.size()); + for (const auto* associations : calrecAssoCollsIn) { + m_calAssoIds.push_back(buildCalAssoId(associations)); + } + + // == s == For MC Trigger Efficiency Estimation ~~~~~~~~ + m_vPhysCollisionTimes.clear(); + + double prevMCTime = -std::numeric_limits::max(); // temp check mc particle times + for (const auto& mcparticle : *m_mcParticles_inCol()) { + if (mcparticle.getGeneratorStatus() != 1) { + continue; + } + if (std::abs(prevMCTime - mcparticle.getTime()) < 50.) { + continue; + } + double mcCollTime = mcparticle.getTime(); + m_vPhysCollisionTimes.push_back(mcCollTime); + prevMCTime = mcCollTime; + } + std::sort(m_vPhysCollisionTimes.begin(), m_vPhysCollisionTimes.end()); + auto last = std::unique(m_vPhysCollisionTimes.begin(), m_vPhysCollisionTimes.end()); + m_vPhysCollisionTimes.erase(last, m_vPhysCollisionTimes.end()); + // == e == For MC Trigger Efficiency Estimation ~~~~~~~~ + } + // == e == Register hits of TOF and MPGD detectors in the time slice ================== + + // == s == Time frame scan loop ========================================================== + double timesliceT0 = std::numeric_limits::quiet_NaN(); + bool bTimesliceTrigger = false; + + std::array bCombineTriggers{}; + std::array combineTrigTime{}; + + if (timesplitWidthValue <= 0.0F) { + throw std::runtime_error("TimeframeSplitter: timesplitWidth must be greater than zero"); + } + if (timeframeWidthValue <= 0.0F) { + throw std::runtime_error("TimeframeSplitter: timeframeWidth must be greater than zero"); + } + + const size_t nTimeSlices = + static_cast(std::floor(timeframeWidthValue / timesplitWidthValue)); + double tsTimeS = 0.0; + double tsTimeE = 0.0; + + // Scan the timeframe one time slice at a time. + // The scan stops early when a physics trigger fires; otherwise it terminates + // after all time slices in the timeframe have been processed. + while (m_iTimeSlice < nTimeSlices) { + tsTimeS = m_iTimeSlice * timesplitWidthValue; + tsTimeE = (m_iTimeSlice + 1) * timesplitWidthValue; + if (tsTimeE > timeframeWidthValue) { + break; + } + m_iTimeSlice++; + + // == s == Multiplicity Single Triggers ======================================= + std::array singleTrig{}; + std::array singleTrigTime{}; + + // Tracker-matching thresholds for each trigger region. + const std::array trackerMatchThresholds = { + m_backwardTrackerMatchThreshold(), + m_barrelTrackerMatchThreshold(), + m_forwardTrackerMatchThreshold(), + }; + + // --------------------------------------------------------------------------- + // Backward, barrel, and forward ECal / ECal+tracker triggers + // --------------------------------------------------------------------------- + for (size_t iRegion = 0; iRegion < kNumSingleTrigRegion; ++iRegion) { + + const auto& config = m_triggerRegionConfigs.at(iRegion); + + // Select the eta-phi binning function corresponding to the trigger region. + const auto binFunc = [this, iRegion](double eta, double phi, int shift) { + switch (iRegion) { + case kSingleTrigRegionBackward: + return backEndEtaPhiBins(eta, phi, shift); + + case kSingleTrigRegionBarrel: + return barrelEtaPhiBins(eta, phi, shift); + + case kSingleTrigRegionForward: + return forwardEndEtaPhiBins(eta, phi, shift); + + default: + throw std::runtime_error("Unknown single trigger region"); + } + }; + + // ------------------------------------------------------------------------- + // ECal trigger + // ------------------------------------------------------------------------- + EtaPhiGrid calGrid{}; + EtaPhiGrid calGridShifted{}; + EtaPhiTimeGrid calTimeGrid{}; + EtaPhiTimeGrid calTimeGridShifted{}; + fillEtaPhiGrids(caloRecHitCollsIn.at(config.calDetector), + m_iniCalHitPoint[config.calDetector], calTimeResolution(config.calDetector), + tsTimeS, tsTimeE, calGrid, calGridShifted, calTimeGrid, calTimeGridShifted, + binFunc); + + singleTrig[config.calTrigger] = countGridCellsWithMultiplicity( + calGrid, calGridShifted, calTimeGrid, calTimeGridShifted, m_ecalMultiplicityThreshold(), + singleTrigTime[config.calTrigger]); + + // ------------------------------------------------------------------------- + // ECal + tracker matching trigger + // ------------------------------------------------------------------------- + EtaPhiGrid trkGrid{}; + EtaPhiGrid trkGridShifted{}; + EtaPhiTimeGrid trkTimeGrid{}; + EtaPhiTimeGrid trkTimeGridShifted{}; + for (const auto trkDetector : config.trkDetectors) { + fillEtaPhiGridsMatched(trackerHitCollsIn.at(trkDetector), m_iniTrkHitPoint[trkDetector], + trkTimeResolution(trkDetector), tsTimeS, tsTimeE, calGrid, + calGridShifted, trkGrid, trkGridShifted, + trackerMatchThresholds.at(iRegion), trkTimeGrid, trkTimeGridShifted, + binFunc); + } + + singleTrig[config.calTrkTrigger] = countGridCellsWithMultiplicity( + trkGrid, trkGridShifted, trkTimeGrid, trkTimeGridShifted, + m_trackerMultiplicityThreshold(), singleTrigTime[config.calTrkTrigger]); + } + + // --------------------------------------------------------------------------- + // B0 tracker trigger + // --------------------------------------------------------------------------- + const auto hitsB0 = + countHitsInTimeWindow(trackerHitCollsIn.at(kTrkB0), m_iniTrkHitPoint[kTrkB0], + trkTimeResolution(kTrkB0), tsTimeS, tsTimeE); + m_iniTrkHitPoint[kTrkB0] = hitsB0.nextStartID; + singleTrig[kSingleTrigB0Trk] = hitsB0.count; + singleTrigTime[kSingleTrigB0Trk] = hitsB0.average_time(); + + double totalZDCEnergy = 0.0; + double totalZDCEnergyTime = 0.0; + const auto* recHitsZDCECal = caloRecHitCollsIn.at(kCalEcalZDC); + if (recHitsZDCECal != nullptr) { + for (size_t iHit = m_iniCalHitPoint[kCalEcalZDC]; iHit < recHitsZDCECal->size(); ++iHit) { + const auto& hit = recHitsZDCECal->at(iHit); + const double hitTime = timeOfFlightCorrectedTime(hit); + if (hitTime - calTimeResolution(kCalEcalZDC) > tsTimeE) { + break; + } + if (judgeHitInTimeSlice(hitTime, calTimeResolution(kCalEcalZDC), tsTimeS, tsTimeE)) { + totalZDCEnergy += hit.getEnergy(); + totalZDCEnergyTime += hit.getEnergy() * hitTime; + m_iniCalHitPoint[kCalEcalZDC] = iHit; + } + } + } + singleTrig[kSingleTrigZDCECal] = totalZDCEnergy; + singleTrigTime[kSingleTrigZDCECal] = + totalZDCEnergy > 0.0 ? totalZDCEnergyTime / totalZDCEnergy : 0.0; + + const double etaPhiCalTriggerSum = singleTrig[kSingleTrigBackEndcapECal] + + singleTrig[kSingleTrigCentBarrelECal] + + singleTrig[kSingleTrigForwardEndcapECal]; + const double etaPhiCalTrkTriggerSum = singleTrig[kSingleTrigBackEndcapECalTrk] + + singleTrig[kSingleTrigCentBarrelECalTrk] + + singleTrig[kSingleTrigForwardEndcapECalTrk]; + bCombineTriggers[kCombTrigECalTrkAndB0Trk] = + etaPhiCalTrkTriggerSum > 0 && singleTrig[kSingleTrigB0Trk] > 4; + bCombineTriggers[kCombTrigECalTrkAndZDCEcal] = + etaPhiCalTrkTriggerSum > 0 && singleTrig[kSingleTrigZDCECal] > 50; + bCombineTriggers[kCombTrigECalAndB0Trk] = + etaPhiCalTriggerSum > 0 && singleTrig[kSingleTrigB0Trk] > 4; + bCombineTriggers[kCombTrigECalAndZDCEcal] = + etaPhiCalTriggerSum > 0 && singleTrig[kSingleTrigZDCECal] > 0.005; + bCombineTriggers[kCombTrigECalTrk] = etaPhiCalTrkTriggerSum > 1; + bCombineTriggers[kCombTrigECal] = etaPhiCalTriggerSum > 2; + + if (std::none_of(bCombineTriggers.begin(), bCombineTriggers.end(), + [](bool fired) { return fired; })) { + continue; + } + + const double fallbackTriggerTime = 0.5 * (tsTimeS + tsTimeE); + if (bCombineTriggers[kCombTrigECalTrkAndB0Trk]) + combineTrigTime[kCombTrigECalTrkAndB0Trk] = + averageSelectedTriggerTime(singleTrig, singleTrigTime, + {kSingleTrigBackEndcapECalTrk, kSingleTrigCentBarrelECalTrk, + kSingleTrigForwardEndcapECalTrk, kSingleTrigB0Trk}, + fallbackTriggerTime); + if (bCombineTriggers[kCombTrigECalTrkAndZDCEcal]) + combineTrigTime[kCombTrigECalTrkAndZDCEcal] = + averageSelectedTriggerTime(singleTrig, singleTrigTime, + {kSingleTrigBackEndcapECalTrk, kSingleTrigCentBarrelECalTrk, + kSingleTrigForwardEndcapECalTrk, kSingleTrigZDCECal}, + fallbackTriggerTime); + if (bCombineTriggers[kCombTrigECalAndB0Trk]) + combineTrigTime[kCombTrigECalAndB0Trk] = + averageSelectedTriggerTime(singleTrig, singleTrigTime, + {kSingleTrigBackEndcapECal, kSingleTrigCentBarrelECal, + kSingleTrigForwardEndcapECal, kSingleTrigB0Trk}, + fallbackTriggerTime); + if (bCombineTriggers[kCombTrigECalAndZDCEcal]) { + combineTrigTime[kCombTrigECalAndZDCEcal] = + averageSelectedTriggerTime(singleTrig, singleTrigTime, + {kSingleTrigBackEndcapECal, kSingleTrigCentBarrelECal, + kSingleTrigForwardEndcapECal, kSingleTrigZDCECal}, + fallbackTriggerTime); + } + if (bCombineTriggers[kCombTrigECalTrk]) { + combineTrigTime[kCombTrigECalTrk] = + averageSelectedTriggerTime(singleTrig, singleTrigTime, + {kSingleTrigBackEndcapECalTrk, kSingleTrigCentBarrelECalTrk, + kSingleTrigForwardEndcapECalTrk}, + fallbackTriggerTime); + } + if (bCombineTriggers[kCombTrigECal]) { + combineTrigTime[kCombTrigECal] = averageSelectedTriggerTime( + singleTrig, singleTrigTime, + {kSingleTrigBackEndcapECal, kSingleTrigCentBarrelECal, kSingleTrigForwardEndcapECal}, + fallbackTriggerTime); + } + double combineTrigCount = 0; + double totalTrigTime = 0.0; + for (size_t iTrig = 0; iTrig < kNumOfCombineTrig; ++iTrig) { + if (bCombineTriggers[iTrig]) { + totalTrigTime += combineTrigTime[iTrig]; + ++combineTrigCount; + } + } + timesliceT0 = combineTrigCount > 0 ? totalTrigTime / static_cast(combineTrigCount) + : fallbackTriggerTime; + + // == s == Multiplisity Single Triggers ======================================= + + bTimesliceTrigger = std::any_of(bCombineTriggers.begin(), bCombineTriggers.end(), + [](bool fired) { return fired; }); + if (bTimesliceTrigger) { + break; + } + } + // == e == Time frame scan loop ========================================================== + + bTrigger = bTimesliceTrigger; + if (bTimesliceTrigger) { + m_bOnceTriggered = true; + // For now, a one-to-one relationship between timeslices and events + child.SetEventNumber(m_NewEventCount); + child.SetRunNumber(parent.GetRunNumber()); + m_NewEventCount++; + // Clone truth particles before detector relations so every child SimTrackerHit can + // point to an MCParticle owned by this PhysicsEvent rather than by the parent Timeslice. + std::vector copiedMCParticles; + copiedMCParticles.reserve(m_mcParticles_inCol()->size()); + std::unordered_map copiedMCParticleMap; + copiedMCParticleMap.reserve(m_mcParticles_inCol()->size()); + for (const auto& mcparticle : *m_mcParticles_inCol()) { + auto copiedMCParticle = mcparticle.clone(false); + m_mcParticles_outCol()->push_back(copiedMCParticle); + copiedMCParticles.push_back(copiedMCParticle); + copiedMCParticleMap.emplace(objIdKey(mcparticle.getObjectID()), + m_mcParticles_outCol()->at(m_mcParticles_outCol()->size() - 1)); + } + + // Recreate the MCParticle graph only after all child particles exist. Keeping this graph is + // required by calorimeter cluster truth association, which walks from each contribution's + // particle to its primary ancestor. + for (size_t particleIndex = 0; particleIndex < m_mcParticles_inCol()->size(); ++particleIndex) { + const auto parentParticle = m_mcParticles_inCol()->at(particleIndex); + auto& copiedParticle = copiedMCParticles.at(particleIndex); + + for (const auto& parentRelation : parentParticle.getParents()) { + const auto copiedParent = copiedMCParticleMap.find(objIdKey(parentRelation.getObjectID())); + if (copiedParent == copiedMCParticleMap.end()) { + throw std::runtime_error("MCParticle parent relation cannot be remapped to child event"); + } + copiedParticle.addToParents(copiedParent->second); + } + + for (const auto& daughterRelation : parentParticle.getDaughters()) { + const auto copiedDaughter = + copiedMCParticleMap.find(objIdKey(daughterRelation.getObjectID())); + if (copiedDaughter == copiedMCParticleMap.end()) { + throw std::runtime_error( + "MCParticle daughter relation cannot be remapped to child event"); + } + copiedParticle.addToDaughters(copiedDaughter->second); + } + } + + // == s == Register Tracker Hits ======================================================= + for (size_t trkDetID = 0; trkDetID < trackerHitCollsIn.size(); ++trkDetID) { + const auto* trkCollIn = trackerHitCollsIn.at(trkDetID); + + if (trkCollIn == nullptr) { + continue; + } + auto& trkCollOut = m_trackerHits_outCols().at(trkDetID); + const auto tempDetID = static_cast(trkDetID); + const double detTimeReso = trkTimeResolution(tempDetID); + const auto* trkAssoCollIn = trkAssoCollsIn.at(trkDetID); + auto& rawCollOut = m_rawTrackerHit_outCols().at(trkDetID); + auto& trkAssoCollOut = m_trackerHitsAsso_outCols().at(trkDetID); + + for (size_t iHit = 0; iHit < trkCollIn->size(); ++iHit) { + const auto& trkHit = trkCollIn->at(iHit); + const double hitT = timeOfFlightCorrectedTime(trkHit); + if (!overlapsTimeWindow(hitT, detTimeReso, timesliceT0 - m_trigTimeWindowBef(), + timesliceT0 + m_trigTimeWindowAft())) { + continue; + } + + m_iniTrkHitPoint[trkDetID] = iHit; + copyTrkHitWithRelations(trkHit, trkAssoCollIn, m_trkAssoIds.at(trkDetID), trkCollOut, + rawCollOut, trkAssoCollOut, m_simTrackerHits_outCols().at(trkDetID), + m_recoTrackerHitLinks_outCols().at(trkDetID), + m_mcParticles_outCol()); + } + } + // == e == Register Tracker Hits ======================================================= + + // == s == Register RICH Hits ========================================================== + // RICH hits are currently not time-filtered. Copy all RICH raw hits from the parent + // Timeslice into every triggered PhysicsEvent together with their truth relations. + for (size_t richDetID = 0; richDetID < richRawHitCollsIn.size(); ++richDetID) { + const auto* rawCollIn = richRawHitCollsIn.at(richDetID); + if (rawCollIn == nullptr) { + continue; + } + + const auto* assoCollIn = richAssoCollsIn.at(richDetID); + + auto& rawCollOut = m_richRawHits_outCols().at(richDetID); + auto& assoCollOut = m_richHitsAsso_outCols().at(richDetID); + auto& linkCollOut = m_richHitLinks_outCols().at(richDetID); + auto& simCollOut = m_richSimHits_outCols().at(richDetID); + + // Multiple raw hits may refer to the same SimTrackerHit. + // Copy each simulated hit only once within this child event. + std::unordered_map copiedSimHits; + + for (size_t iRawHit = 0; iRawHit < rawCollIn->size(); ++iRawHit) { + const auto& rawHit = rawCollIn->at(iRawHit); + + auto copiedRawHit = rawHit.clone(); + rawCollOut->push_back(copiedRawHit); + + if (assoCollIn == nullptr) { + continue; + } + + const auto assocIter = m_richAssoIds.at(richDetID).find(objIdKey(rawHit.getObjectID())); + + if (assocIter == m_richAssoIds.at(richDetID).end()) { + continue; + } + + for (const size_t associationPosition : assocIter->second) { + const auto association = assoCollIn->at(associationPosition); + + if (!association.getSimHit().isAvailable()) { + continue; + } + + const auto simHit = association.getSimHit(); + const auto simHitKey = objIdKey(simHit.getObjectID()); + + auto copiedSimHitIter = copiedSimHits.find(simHitKey); + + if (copiedSimHitIter == copiedSimHits.end()) { + auto copiedSimHit = simHit.clone(false); + + // Remap SimTrackerHit -> MCParticle to the particle owned by this child event. + if (simHit.getParticle().isAvailable()) { + const auto copiedParticle = + copiedMCParticleMap.find(objIdKey(simHit.getParticle().getObjectID())); + + if (copiedParticle == copiedMCParticleMap.end()) { + throw std::runtime_error( + "RICH SimTrackerHit particle relation cannot be remapped to child event"); + } + + copiedSimHit.setParticle(copiedParticle->second); + } + + simCollOut->push_back(copiedSimHit); + + copiedSimHitIter = + copiedSimHits.emplace(simHitKey, simCollOut->at(simCollOut->size() - 1)).first; + } + + const auto copiedSimHit = copiedSimHitIter->second; + + auto copiedAssociation = assoCollOut->create(); + copiedAssociation.setWeight(association.getWeight()); + copiedAssociation.setRawHit(copiedRawHit); + copiedAssociation.setSimHit(copiedSimHit); + + auto copiedLink = linkCollOut->create(); + copiedLink.setWeight(association.getWeight()); + copiedLink.setFrom(copiedRawHit); + copiedLink.setTo(copiedSimHit); + } + } + } + // == e == Register RICH Hits ========================================================== + + // == s == Register Calo Rec Hits ======================================================= + for (size_t calDetID = 0; calDetID < caloRecHitCollsIn.size(); ++calDetID) { + const auto* caloInColl = caloRecHitCollsIn.at(calDetID); + if (caloInColl == nullptr) { + continue; + } + auto& caloOutColl = m_calorimeterHit_outCols().at(calDetID); + + const auto* caloInCollAsso = calrecAssoCollsIn.at(calDetID); + if (caloInCollAsso == nullptr) { + continue; + } + + // Several raw hits can refer to the same simulated hit. Copy each simulated hit and its + // contribution graph once per child event and detector, then reuse the child-owned proxy + // for all raw-hit links and associations. + std::unordered_map copiedSimHits; + std::unordered_map copiedContributions; + + for (size_t iCalHit = 0; iCalHit < caloInColl->size(); ++iCalHit) { + const auto& caloHit = caloInColl->at(iCalHit); + + double detTimeReso = calTimeResolution(kCalEcalEndcapN); // ??? check ECal Time resolution + const double hitT = timeOfFlightCorrectedTime(caloHit); + + if (hitT - detTimeReso > timesliceT0 + m_trigTimeWindowAft()) { + continue; + } + if (overlapsTimeWindow(hitT, detTimeReso, timesliceT0 - m_trigTimeWindowBef(), + timesliceT0 + m_trigTimeWindowAft())) { + auto copiedCaloHit = caloHit.clone(); + copiedCaloHit.setRawHit(edm4hep::RawCalorimeterHit()); + + const auto rawHitFromRec = caloHit.getRawHit(); + if (rawHitFromRec.isAvailable()) { + auto& rawCollOut = m_rawCalorimeterHit_outCols().at(calDetID); + auto copiedRawHit = rawHitFromRec.clone(); + rawCollOut->push_back(copiedRawHit); + copiedCaloHit.setRawHit(copiedRawHit); + + auto& assocCollOut = m_mcRecoCalorimeterHitAssociation_outCols().at(calDetID); + auto& linkCollOut = m_mcRecoCalorimeterHitLink_outCols().at(calDetID); + auto& simCollOut = m_simCalorimeterHit_outCols().at(calDetID); + const auto rawHitID = rawHitFromRec.getObjectID(); + + const auto& association_index = m_calAssoIds.at(calDetID); + const auto assocIterCal = association_index.find(objIdKey(rawHitID)); + + if (assocIterCal != association_index.end()) { + for (const size_t association_position : assocIterCal->second) { + const auto assoc = caloInCollAsso->at(association_position); + + if (!assoc.getSimHit().isAvailable()) { + continue; + } + + const auto simHit = assoc.getSimHit(); + const auto simHitKey = objIdKey(simHit.getObjectID()); + + auto copiedSimHitIter = copiedSimHits.find(simHitKey); + if (copiedSimHitIter == copiedSimHits.end()) { + auto copiedSimHit = simHit.clone(false); + auto& contributionCollOut = m_caloHitContribution_outCols().at(calDetID); + + for (const auto& contribution : simHit.getContributions()) { + const auto contributionKey = objIdKey(contribution.getObjectID()); + auto copiedContributionIter = copiedContributions.find(contributionKey); + + if (copiedContributionIter == copiedContributions.end()) { + auto copiedContribution = contribution.clone(false); + copiedContribution.setParticle(edm4hep::MCParticle()); + + const auto particle = contribution.getParticle(); + if (particle.isAvailable()) { + const auto copiedParticle = + copiedMCParticleMap.find(objIdKey(particle.getObjectID())); + if (copiedParticle == copiedMCParticleMap.end()) { + throw std::runtime_error( + "CaloHitContribution particle relation cannot be remapped to child " + "event"); + } + copiedContribution.setParticle(copiedParticle->second); + } + + contributionCollOut->push_back(copiedContribution); + copiedContributionIter = + copiedContributions + .emplace(contributionKey, + contributionCollOut->at(contributionCollOut->size() - 1)) + .first; + } + + copiedSimHit.addToContributions(copiedContributionIter->second); + } + + simCollOut->push_back(copiedSimHit); + copiedSimHitIter = + copiedSimHits.emplace(simHitKey, simCollOut->at(simCollOut->size() - 1)) + .first; + } + + const auto copiedSimHit = copiedSimHitIter->second; + + auto copiedAssoc = assocCollOut->create(); + copiedAssoc.setWeight(assoc.getWeight()); + copiedAssoc.setRawHit(copiedRawHit); + copiedAssoc.setSimHit(copiedSimHit); + + auto copiedLink = linkCollOut->create(); + copiedLink.setWeight(assoc.getWeight()); + copiedLink.setFrom(copiedRawHit); + copiedLink.setTo(copiedSimHit); + } + } + } + + caloOutColl->push_back(copiedCaloHit); + m_iniCalHitPoint[calDetID] = iCalHit; + } + } + } + // == e == Register Calo Rec Hits ======================================================= + + // == s == For QA relation valuables QA<><><><><><><><><><><><><><><><><><>> + // == s == For MC Trigger Efficiency Estimation ~~~~~~~~ + unsigned int physEventWeight = 2; + for (auto it = m_vPhysCollisionTimes.begin(); it != m_vPhysCollisionTimes.end(); ++it) { + const double physCollTime = *it; + if ((physCollTime + m_collisionTimeMarginAft() > timesliceT0 - m_trigTimeWindowBef()) && + (physCollTime - m_collisionTimeMarginBef() < timesliceT0 + m_trigTimeWindowAft())) { + physEventWeight = 1; + m_vPhysCollisionTimes.erase(it); + break; + } + } + + if (physEventWeight == 1) { + edm4hep::MutableEventHeader eventHeader_phy; + eventHeader_phy.setRunNumber(child_idx); + eventHeader_phy.setEventNumber(m_eventNumber_TS); + eventHeader_phy.setTimeStamp(m_iTimeSlice); + eventHeader_phy.setWeight(2); + m_eventHeaderPhy_outCols()->push_back(eventHeader_phy); + + edm4hep::MutableEventHeader eventHeader_bkg; + eventHeader_bkg.setRunNumber(child_idx); + eventHeader_bkg.setEventNumber(m_eventNumber_TS); + eventHeader_bkg.setTimeStamp(m_iTimeSlice); + eventHeader_bkg.setWeight(1); + m_eventHeaderBkg_outCols()->push_back(eventHeader_bkg); + + for (size_t iTrig = 0; iTrig < kNumOfCombineTrig; ++iTrig) { + if (bCombineTriggers[iTrig]) { + edm4hep::MutableEventHeader eventHeader_phy; + eventHeader_phy.setRunNumber(child_idx); + eventHeader_phy.setEventNumber(m_eventNumber_TS); + eventHeader_phy.setTimeStamp(m_iTimeSlice); + eventHeader_phy.setWeight(iTrig + 3); + m_eventHeaderPhy_outCols()->push_back(eventHeader_phy); + } + } + + m_PhysCount++; + } else if (physEventWeight == 2) { + edm4hep::MutableEventHeader eventHeader_phy; + eventHeader_phy.setRunNumber(child_idx); + eventHeader_phy.setEventNumber(m_eventNumber_TS); + eventHeader_phy.setTimeStamp(m_iTimeSlice); + eventHeader_phy.setWeight(1); + m_eventHeaderPhy_outCols()->push_back(eventHeader_phy); + + edm4hep::MutableEventHeader eventHeader_bkg; + eventHeader_bkg.setRunNumber(child_idx); + eventHeader_bkg.setEventNumber(m_eventNumber_TS); + eventHeader_bkg.setTimeStamp(m_iTimeSlice); + eventHeader_bkg.setWeight(2); + m_eventHeaderBkg_outCols()->push_back(eventHeader_bkg); + for (size_t iTrig = 0; iTrig < kNumOfCombineTrig; ++iTrig) { + if (bCombineTriggers[iTrig]) { + edm4hep::MutableEventHeader eventHeader_bkg; + eventHeader_bkg.setRunNumber(child_idx); + eventHeader_bkg.setEventNumber(m_eventNumber_TS); + eventHeader_bkg.setTimeStamp(m_iTimeSlice); + eventHeader_bkg.setWeight(iTrig + 3); + m_eventHeaderBkg_outCols()->push_back(eventHeader_bkg); + } + } + } + m_eventNumber_TS++; + + // Insert an independent EventHeader object into the physics event. + // A subset header would keep a reference to the parent frame collection. + edm4hep::MutableEventHeader eventHeader; + if (m_eventHeader_inCol() != nullptr && !m_eventHeader_inCol()->empty()) { + const auto& eventHeader_in = m_eventHeader_inCol()->at(0); + eventHeader.setRunNumber(eventHeader_in.getRunNumber()); + eventHeader.setEventNumber(eventHeader_in.getEventNumber()); + eventHeader.setTimeStamp(eventHeader_in.getTimeStamp()); + eventHeader.setWeight(eventHeader_in.getWeight()); + } else { + eventHeader.setRunNumber(child.GetRunNumber()); + eventHeader.setEventNumber(child.GetEventNumber()); + eventHeader.setTimeStamp(m_iTimeSlice); + eventHeader.setWeight(physEventWeight); + } + m_eventHeader_outCol()->push_back(eventHeader); + // == e == For MC Trigger Efficiency Estimation ~~~~~~~~ + + // == s == For QA relation valuables QA<><><><><><><><><><><><><><><><><>> + } + + if (m_iTimeSlice >= nTimeSlices) { + m_bScanedAllTimeWindows = true; + } + if (m_bScanedAllTimeWindows) { + m_bInitialLoop = true; + m_bOnceTriggered = false; + + m_vPhysCollisionTimes.clear(); + + m_bScanedAllTimeWindows = false; + m_iTimeSlice = 0; + m_targetDetId = 0; + for (auto& start_point : m_iniTrkHitPoint) { + start_point = 0; + } + for (auto& start_point : m_iniCalHitPoint) { + start_point = 0; + } + + if (bTrigger) { + return Result::NextChildNextParent; + } else { + return Result::KeepChildNextParent; + } + + } else if (bTrigger) { + return Result::NextChildKeepParent; + } + return Result::KeepChildNextParent; +} diff --git a/src/global/splitting/TimeframeSplitter.h b/src/global/splitting/TimeframeSplitter.h new file mode 100644 index 0000000000..c6af11bf39 --- /dev/null +++ b/src/global/splitting/TimeframeSplitter.h @@ -0,0 +1,876 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2026 Takuya Kumaoka + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +struct TimeframeSplitter : public JEventUnfolder { +public: + struct ConfigT { + float timeframeWidth = 2000.0; + float timesplitWidth = 20.0; + float timeResolution_SiMaps = 2000.0; + float timeResolution_MPGD = 30.0; + float timeResolution_ACLGad = 20.0; + float timeResolution_EMCal = 20.0; + float timeResolution_HCal = 100.0; + double refInverseVelocity = 0.0034; + double backwardEtaMin = -3.78; + double backwardEtaMax = -1.63; + double barrelEtaMin = 1.80; + double barrelEtaMax = 1.81; + double forwardEtaMin = 1.77; + double forwardEtaMax = 4.04; + size_t ecalMultiplicityThreshold = 10; + size_t backwardTrackerMatchThreshold = 10; + size_t barrelTrackerMatchThreshold = 5; + size_t forwardTrackerMatchThreshold = 5; + size_t trackerMultiplicityThreshold = 1; + double trigTimeWindowBef = 10 * edm4eic::unit::ns; + double trigTimeWindowAft = 30 * edm4eic::unit::ns; + double collisionTimeMarginBef = 10 * edm4eic::unit::ns; + double collisionTimeMarginAft = 20 * edm4eic::unit::ns; + }; + +private: + ConfigT m_config; + + ParameterRef m_timeframeWidth{this, "timeframe_width", m_config.timeframeWidth, + "Width of each timeframe in ns"}; + ParameterRef m_timesplitWidth{this, "timesplit_width", m_config.timesplitWidth, + "Width of each timeslice in ns"}; + ParameterRef m_timeResolution_SiMaps{this, "timeResolution_Silicon", + m_config.timeResolution_SiMaps, + "time resolution of Silicon detector in ns"}; + ParameterRef m_timeResolution_MPGD{this, "timeResolution_MPGD", + m_config.timeResolution_MPGD, + "time resolution of MPGD detector in ns"}; + ParameterRef m_timeResolution_ACLGad{this, "timeResolution_TOF", + m_config.timeResolution_ACLGad, + "time resolution of TOF detector in ns"}; + ParameterRef m_timeResolution_EMCal{this, "timeResolution_EMCal", + m_config.timeResolution_EMCal, + "time resolution of EMCal detector in ns"}; + ParameterRef m_timeResolution_HCal{this, "timeResolution_HCal", + m_config.timeResolution_HCal, + "time resolution of HCal detector in ns"}; + + ParameterRef m_refInverseVelocity{ + this, "refInverseVelocity", m_config.refInverseVelocity, + "ns/mm estimated by MC average time of flight"}; //< ns/mm estimated by MC average time of flight / distance from IP to calorimeter + + // MPGD backward Endcap range -3.6 < eta < -1.72, +5%: -3.78 < eta < -1.634 + ParameterRef m_backwardEtaMin{this, "backward_eta_min", m_config.backwardEtaMin, + "Minimum eta for the backward trigger region"}; + ParameterRef m_backwardEtaMax{this, "backward_eta_max", m_config.backwardEtaMax, + "Maximum eta for the backward trigger region"}; + + // Barrel trigger region covering MPGD, TOF, and ECal acceptance with margin + ParameterRef m_barrelEtaMin{this, "barrel_eta_min", m_config.barrelEtaMin, + "Minimum eta for the barrel trigger region"}; + ParameterRef m_barrelEtaMax{this, "barrel_eta_max", m_config.barrelEtaMax, + "Maximum eta for the barrel trigger region"}; + + // Forward trigger region covering MPGD, TOF, and ECal acceptance with margin + ParameterRef m_forwardEtaMin{this, "forward_eta_min", m_config.forwardEtaMin, + "Minimum eta for the forward trigger region"}; + ParameterRef m_forwardEtaMax{this, "forward_eta_max", m_config.forwardEtaMax, + "Maximum eta for the forward trigger region"}; + + std::pair backEndEtaPhiBins(double hitEta, double hitPhi, int bShift); + std::pair barrelEtaPhiBins(double hitEta, double hitPhi, int bShift); + std::pair forwardEndEtaPhiBins(double hitEta, double hitPhi, int bShift); + + ParameterRef m_ecalMultiplicityThreshold{ + this, "ecal_multiplicity_threshold", m_config.ecalMultiplicityThreshold, + "Minimum ECal grid-cell multiplicity for single triggers"}; + + ParameterRef m_backwardTrackerMatchThreshold{ + this, "backward_tracker_match_threshold", m_config.backwardTrackerMatchThreshold, + "Tracker matching threshold for the backward trigger region"}; + + ParameterRef m_barrelTrackerMatchThreshold{ + this, "barrel_tracker_match_threshold", m_config.barrelTrackerMatchThreshold, + "Tracker matching threshold for the barrel trigger region"}; + + ParameterRef m_forwardTrackerMatchThreshold{ + this, "forward_tracker_match_threshold", m_config.forwardTrackerMatchThreshold, + "Tracker matching threshold for the forward trigger region"}; + + ParameterRef m_trackerMultiplicityThreshold{ + this, "tracker_multiplicity_threshold", m_config.trackerMultiplicityThreshold, + "Minimum matched tracker grid-cell multiplicity for single triggers"}; + + ParameterRef m_trigTimeWindowBef{this, "trigger_window_before", + m_config.trigTimeWindowBef, + "Time window before the trigger time"}; + ParameterRef m_trigTimeWindowAft{this, "trigger_window_after", m_config.trigTimeWindowAft, + "Time window after the trigger time"}; + + ParameterRef m_collisionTimeMarginBef{this, "collision_time_margin_before", + m_config.collisionTimeMarginBef, + "Time margin before the collision time"}; + ParameterRef m_collisionTimeMarginAft{this, "collision_time_margin_after", + m_config.collisionTimeMarginAft, + "Time margin after the collision time"}; + + PodioInput m_eventHeader_inCol{ + this, {.name = "EventHeader", .is_optional = true}}; + PodioOutput m_eventHeader_outCol{this, "EventHeader"}; + + PodioInput m_mcParticles_inCol{this, + {.name = "MCParticles", .is_optional = true}}; + PodioOutput m_mcParticles_outCol{this, "MCParticles"}; + + // tracker collections + VariadicPodioInput m_trackerHits_inCols{ + this, {.names = getTrkCollectionNames(kTrackerHit), .is_optional = true}}; + VariadicPodioOutput m_trackerHits_outCols{ + this, getTrkCollectionNames(kTrackerHit)}; + + VariadicPodioInput m_trackerHitsAsso_inCols{ + this, {.names = getTrkCollectionNames(kTrackerHitAssociation), .is_optional = true}}; + VariadicPodioOutput m_trackerHitsAsso_outCols{ + this, getTrkCollectionNames(kTrackerHitAssociation)}; + + VariadicPodioOutput m_recoTrackerHitLinks_outCols{ + this, getTrkCollectionNames(kTrackerHitLink)}; + + VariadicPodioOutput m_simTrackerHits_outCols{ + this, getTrkCollectionNames(kSimTrackerHit)}; + + VariadicPodioInput m_rawTrackerHit_inCols{ + this, {.names = getTrkCollectionNames(kRawTrackerHit), .is_optional = true}}; + VariadicPodioOutput m_rawTrackerHit_outCols{ + this, getTrkCollectionNames(kRawTrackerHit)}; + + // RICH collections + VariadicPodioInput m_richHitsAsso_inCols{ + this, {.names = getRichCollectionNames(kRichRawHitAssociation), .is_optional = true}}; + VariadicPodioOutput m_richHitsAsso_outCols{ + this, getRichCollectionNames(kRichRawHitAssociation)}; + + VariadicPodioOutput m_richHitLinks_outCols{ + this, getRichCollectionNames(kRichRawHitLink)}; + + VariadicPodioOutput m_richSimHits_outCols{ + this, getRichCollectionNames(kRichSimTrackerHit)}; + + VariadicPodioInput m_richRawHits_inCols{ + this, {.names = getRichCollectionNames(kRichRawTrackerHit), .is_optional = true}}; + VariadicPodioOutput m_richRawHits_outCols{ + this, getRichCollectionNames(kRichRawTrackerHit)}; + + // calorimeter collections + VariadicPodioInput m_rawCalorimeterHit_inCols{ + this, {.names = getCalCollectionNames(kRawCalorimeterHit), .is_optional = true}}; + VariadicPodioOutput m_rawCalorimeterHit_outCols{ + this, getCalCollectionNames(kRawCalorimeterHit)}; + + VariadicPodioOutput m_mcRecoCalorimeterHitLink_outCols{ + this, getCalCollectionNames(kCalorimeterHitLink)}; + VariadicPodioOutput m_simCalorimeterHit_outCols{ + this, getCalCollectionNames(kSimCalorimeterHit)}; + VariadicPodioOutput m_caloHitContribution_outCols{ + this, getCalContributionCollectionNames()}; + + VariadicPodioInput m_calorimeterHit_inCols{ + this, {.names = getCalCollectionNames(kCalorimeterHit), .is_optional = true}}; + VariadicPodioOutput m_calorimeterHit_outCols{ + this, getCalCollectionNames(kCalorimeterHit)}; + + VariadicPodioInput + m_mcRecoCalorimeterHitAssociation_inCols{ + this, {.names = getCalCollectionNames(kCalorimeterHitAssociation), .is_optional = true}}; + VariadicPodioOutput + m_mcRecoCalorimeterHitAssociation_outCols{this, + getCalCollectionNames(kCalorimeterHitAssociation)}; + + PodioOutput m_eventHeaderPhy_outCols{this, "EventHeader_PHY"}; + PodioOutput m_eventHeaderBkg_outCols{this, "EventHeader_BKG"}; + + std::shared_ptr m_log; + + size_t m_eventNumber_TS = 0; // Event number for the current timeslice + + static constexpr int kEtaPhiBins = 10; + static constexpr int kInvalidEtaPhiBin = -1; + + using EtaPhiGrid = std::array, kEtaPhiBins>; + using EtaPhiTimeGrid = std::array, kEtaPhiBins>; + using EtaPhiEnergyGrid = std::array, kEtaPhiBins>; + + enum TrkCollectionType : size_t { + kTrackerHit = 0, + kTrackerHitAssociation, + kTrackerHitLink, + kSimTrackerHit, + kRawTrackerHit, + kTrkCollectionTypeSize + }; + + enum RichCollectionType : size_t { + kRichRawHitAssociation = 0, + kRichRawHitLink, + kRichSimTrackerHit, + kRichRawTrackerHit, + kRichCollectionTypeSize + }; + + enum CalCollectionType : size_t { + kCalorimeterHit = 0, + kCalorimeterHitAssociation, + kCalorimeterHitLink, + kSimCalorimeterHit, + kRawCalorimeterHit, + kCalCollectionTypeSize + }; + + enum TrkCollectionIndex : size_t { + kTrkB0 = 0, + kTrkTOFBarrel, + kTrkTOFEndcap, + kTrkMPGDBarrel, + kTrkOuterMPGDBarrel, + kTrkBackwardMPGD, + kTrkForwardMPGD, + kTrkSiBarrelVertex, + kTrkSiBarrel, + kTrkSiEndcap, + kTrkTagger, + kTrkForwardRomanPot, + kTrkForwardOffMTracker, + kTrkCollectionSize + }; + + enum RichCollectionIndex : size_t { kRICHEndcapN = 0, kDIRCBar, kDRICH, kRichCollectionSize }; + + enum CalCollectionIndex : size_t { + kCalB0ECal = 0, + kCalEcalBarrelImg, + kCalEcalBarrelScFi, + kCalEcalEndcapN, + kCalEcalEndcapP, + kCalEcalZDC, + kCalEcalLumiSpec, + kCalHcalBarrel, + kCalHcalEndcapN, + kCalHcalEndcapPInsert, + kCalHcalZDC, + kCalLFHCAL, + kCalCollectionSize + }; + + enum SingleTriggerIndex : size_t { + kSingleTrigBackEndcapECal = 0, + kSingleTrigBackEndcapECalTrk, + kSingleTrigCentBarrelECal, + kSingleTrigCentBarrelECalTrk, + kSingleTrigForwardEndcapECal, + kSingleTrigForwardEndcapECalTrk, + kSingleTrigB0Trk, + kSingleTrigZDCECal, + kNumOfSingleTrig + }; + + enum SingleTriggerRegion : size_t { + kSingleTrigRegionBackward = 0, + kSingleTrigRegionBarrel, + kSingleTrigRegionForward, + kNumSingleTrigRegion + }; + + struct TriggerRegionConfig { + CalCollectionIndex calDetector; + std::vector trkDetectors; + SingleTriggerIndex calTrigger; + SingleTriggerIndex calTrkTrigger; + }; + + const std::array m_triggerRegionConfigs = {{ + { + kCalEcalEndcapN, + {kTrkBackwardMPGD}, + kSingleTrigBackEndcapECal, + kSingleTrigBackEndcapECalTrk, + }, + { + kCalEcalBarrelScFi, + {kTrkMPGDBarrel, kTrkOuterMPGDBarrel, kTrkTOFBarrel}, + kSingleTrigCentBarrelECal, + kSingleTrigCentBarrelECalTrk, + }, + { + kCalEcalEndcapP, + {kTrkForwardMPGD, kTrkTOFEndcap}, + kSingleTrigForwardEndcapECal, + kSingleTrigForwardEndcapECalTrk, + }, + }}; + + enum CombineTriggerIndex : size_t { + kCombTrigECalTrkAndB0Trk = 0, + kCombTrigECalTrkAndZDCEcal, + kCombTrigECalAndB0Trk, + kCombTrigECalAndZDCEcal, + kCombTrigECalTrk, + kCombTrigECal, + kNumOfCombineTrig + }; + + using TrkCollNames = std::array; + std::array m_trkCollNames = {{ + { + "B0TrackerRecHits", + "B0TrackerRawHitAssociations", + "B0TrackerRawHitLinks", + "B0TrackerHits", + "B0TrackerRawHits", + }, + { + "TOFBarrelSharedRecHits", + "TOFBarrelSharedRawHitAssociations", + "TOFBarrelSharedRawHitLinks", + "TOFBarrelHits", + "TOFBarrelSharedRawHits", + }, + { + "TOFEndcapSharedRecHits", + "TOFEndcapSharedRawHitAssociations", + "TOFEndcapSharedRawHitLinks", + "TOFEndcapHits", + "TOFEndcapSharedRawHits", + }, + { + "MPGDBarrelRecHits", + "MPGDBarrelRawHitAssociations", + "MPGDBarrelRawHitLinks", + "MPGDBarrelHits", + "MPGDBarrelRawHits", + }, + { + "OuterMPGDBarrelRecHits", + "OuterMPGDBarrelRawHitAssociations", + "OuterMPGDBarrelRawHitLinks", + "OuterMPGDBarrelHits", + "OuterMPGDBarrelRawHits", + }, + { + "BackwardMPGDEndcapRecHits", + "BackwardMPGDEndcapRawHitAssociations", + "BackwardMPGDEndcapRawHitLinks", + "BackwardMPGDEndcapHits", + "BackwardMPGDEndcapRawHits", + }, + { + "ForwardMPGDEndcapRecHits", + "ForwardMPGDEndcapRawHitAssociations", + "ForwardMPGDEndcapRawHitLinks", + "ForwardMPGDEndcapHits", + "ForwardMPGDEndcapRawHits", + }, + { + "SiBarrelVertexRecHits", + "SiBarrelVertexRawHitAssociations", + "SiBarrelVertexRawHitLinks", + "VertexBarrelHits", + "SiBarrelVertexRawHits", + }, + { + "SiBarrelTrackerRecHits", + "SiBarrelRawHitAssociations", + "SiBarrelRawHitLinks", + "SiBarrelHits", + "SiBarrelRawHits", + }, + { + "SiEndcapTrackerRecHits", + "SiEndcapTrackerRawHitAssociations", + "SiEndcapTrackerRawHitLinks", + "TrackerEndcapHits", + "SiEndcapTrackerRawHits", + }, + { + "TaggerTrackerRecHits", + "TaggerTrackerRawHitAssociations", + "TaggerTrackerRawHitLinks", + "TaggerTrackerHits", + "TaggerTrackerRawHits", + }, + { + "ForwardRomanPotRecHits", + "ForwardRomanPotRawHitAssociations", + "ForwardRomanPotRawHitLinks", + "ForwardRomanPotHits", + "ForwardRomanPotRawHits", + }, + { + "ForwardOffMTrackerRecHits", + "ForwardOffMTrackerRawHitAssociations", + "ForwardOffMTrackerRawHitLinks", + "ForwardOffMTrackerHits", + "ForwardOffMTrackerRawHits", + }, + }}; + + using RichCollNames = std::array; + std::array m_richCollNames = {{ + { + "RICHEndcapNRawHitsAssociations", + "RICHEndcapNRawHitsLinks", + "PFRICHHits", + "RICHEndcapNRawHits", + }, + { + "DIRCRawHitsAssociations", + "DIRCRawHitsLinks", + "DIRCBarHits", + "DIRCRawHits", + }, + { + "DRICHRawHitsAssociations", + "DRICHRawHitsLinks", + "DRICHHits", + "DRICHRawHits", + }, + }}; + + using CalCollNames = std::array; + std::array m_calCollNames = {{ + { + "B0ECalRecHits", + "B0ECalRawHitAssociations", + "B0ECalRawHitLinks", + "B0ECalHits", + "B0ECalRawHits", + }, + { + "EcalBarrelImagingRecHits", + "EcalBarrelImagingRawHitAssociations", + "EcalBarrelImagingRawHitLinks", + "EcalBarrelImagingHits", + "EcalBarrelImagingRawHits", + }, + { + "EcalBarrelScFiRecHits", + "EcalBarrelScFiRawHitAssociations", + "EcalBarrelScFiRawHitLinks", + "EcalBarrelScFiHits", + "EcalBarrelScFiRawHits", + }, + { + "EcalEndcapNRecHits", + "EcalEndcapNRawHitAssociations", + "EcalEndcapNRawHitLinks", + "EcalEndcapNHits", + "EcalEndcapNRawHits", + }, + { + "EcalEndcapPRecHits", + "EcalEndcapPRawHitAssociations", + "EcalEndcapPRawHitLinks", + "EcalEndcapPHits", + "EcalEndcapPRawHits", + }, + { + "EcalFarForwardZDCRecHits", + "EcalFarForwardZDCRawHitAssociations", + "EcalFarForwardZDCRawHitLinks", + "EcalFarForwardZDCHits", + "EcalFarForwardZDCRawHits", + }, + { + "EcalLumiSpecRecHits", + "EcalLumiSpecRawHitAssociations", + "EcalLumiSpecRawHitLinks", + "EcalLumiSpecHits", + "EcalLumiSpecRawHits", + }, + { + "HcalBarrelRecHits", + "HcalBarrelRawHitAssociations", + "HcalBarrelRawHitLinks", + "HcalBarrelHits", + "HcalBarrelRawHits", + }, + { + "HcalEndcapNRecHits", + "HcalEndcapNRawHitAssociations", + "HcalEndcapNRawHitLinks", + "HcalEndcapNHits", + "HcalEndcapNRawHits", + }, + { + "HcalEndcapPInsertRecHits", + "HcalEndcapPInsertRawHitAssociations", + "HcalEndcapPInsertRawHitLinks", + "HcalEndcapPInsertHits", + "HcalEndcapPInsertRawHits", + }, + { + "HcalFarForwardZDCRecHits", + "HcalFarForwardZDCRawHitAssociations", + "HcalFarForwardZDCRawHitLinks", + "HcalFarForwardZDCHits", + "HcalFarForwardZDCRawHits", + }, + { + "LFHCALRecHits", + "LFHCALRawHitAssociations", + "LFHCALRawHitLinks", + "LFHCALHits", + "LFHCALRawHits", + }, + }}; + + std::vector> + m_hitStartIndices_simTracker; + std::vector> + m_hitStartIndices_simCalorimeter; + + // == Global Variables ======================= + unsigned int m_NewEventCount = 0; + unsigned int m_PhysCount = 0; + + bool m_bInitialLoop = true; + + size_t m_iniTrkHitPoint[15] = {0}; // B0Trk, + size_t m_iniCalHitPoint[15] = {0}; // B0Trk, + + bool m_bOnceTriggered = false; + bool m_bScanedAllTimeWindows = false; + + unsigned int m_targetDetId = 0; + size_t m_iTimeSlice = 0; + std::vector m_vPhysCollisionTimes = {}; + // == Global Variables ======================= + + struct TimeWindowSummary { + size_t count = 0; + double timeSum = 0.0; + size_t nextStartID = 0; + + double average_time() const { return count == 0 ? 0.0 : timeSum / count; } + }; + + using TrackerAssociationIndex = std::unordered_map>; + using CalorimeterAssociationIndex = std::unordered_map>; + + std::vector m_trkAssoIds; + std::vector m_richAssoIds; + std::vector m_calAssoIds; + +public: + TimeframeSplitter(); + + /// Retrieve reference to embedded config object + ConfigT& config() { return m_config; } + + Result Unfold(const JEvent& parent, JEvent& child, int child_idx) override; + +protected: + std::vector getTrkCollectionNames(TrkCollectionType type) const { + std::vector names; + names.reserve(m_trkCollNames.size()); + for (const auto& collections : m_trkCollNames) { + names.push_back(collections[type]); + } + return names; + } + + std::vector getRichCollectionNames(RichCollectionType type) const { + std::vector names; + names.reserve(m_richCollNames.size()); + for (const auto& collections : m_richCollNames) { + names.push_back(collections[type]); + } + return names; + } + + std::vector getCalCollectionNames(CalCollectionType type) const { + std::vector names; + names.reserve(m_calCollNames.size()); + for (const auto& collections : m_calCollNames) { + names.push_back(collections[type]); + } + return names; + } + + std::vector getCalContributionCollectionNames() const { + std::vector names; + names.reserve(m_calCollNames.size()); + for (const auto& collections : m_calCollNames) { + names.push_back(collections[kSimCalorimeterHit] + "Contributions"); + } + return names; + } + + static std::uint64_t objIdKey(const podio::ObjectID& object_id); + + static TrackerAssociationIndex + buildTrkAssoId(const edm4eic::MCRecoTrackerHitAssociationCollection* associations); + + static CalorimeterAssociationIndex + buildCalAssoId(const edm4eic::MCRecoCalorimeterHitAssociationCollection* associations); + + static bool overlapsTimeWindow(double hitTime, double resolution, double window_start, + double window_end); + static bool judgeOverTimeWindow(double hitTime, double resolution, double window_end); + + static bool isValidEtaPhiBin(int etaBin, int phiBin); + + static bool judgeHitInTimeSlice(double hitTime, double timeResolution, double timeslice_start, + double timeslice_end); + + template double timeOfFlightCorrectedTime(const HitT& hit) { + const auto& position = hit.getPosition(); + const double radius = std::sqrt(position[0] * position[0] + position[1] * position[1] + + position[2] * position[2]); + + return hit.getTime() - radius * m_refInverseVelocity(); + } + + template inline void etaPhiCalc(const HitT& hit, double& hitEta, double& hitPhi) { + const double hitX = hit.getPosition()[0]; + const double hitY = hit.getPosition()[1]; + const double hitZ = hit.getPosition()[2]; + const double hitR = std::sqrt(hitX * hitX + hitY * hitY + hitZ * hitZ); + if (hitR <= 0.0) { + hitEta = 0.0; + hitPhi = 0.0; + return; + } + + const double cosTheta = std::clamp(hitZ / hitR, -1.0, 1.0); + const double hitTheta = std::acos(cosTheta); + + hitEta = -std::log(std::tan(hitTheta / 2.0)); + hitPhi = std::atan2(hitY, hitX); + } + + static std::pair etaPhiBins(double hitEta, double hitPhi, double etaMin, double etaMax, + int bShift); + + template + void fillEtaPhiGrids(const CollectionT* hits, size_t& iniHitID, double timeResolution, + double timeSliceStart, double timeSliceEnd, EtaPhiGrid& grid, + EtaPhiGrid& gridShifted, EtaPhiTimeGrid& gridTime, + EtaPhiTimeGrid& gridShiftedTime, BinFunc binFunc) { + if (hits == nullptr) { + return; + } + + const size_t hitCount = hits->size(); + for (size_t iHit = iniHitID; iHit < hitCount; ++iHit) { + const auto& hit = hits->at(iHit); + const double hitT = timeOfFlightCorrectedTime(hit); + if (hitT - timeResolution > timeSliceEnd) { + iniHitID = iHit; + break; + } + if (!judgeHitInTimeSlice(hitT, timeResolution, timeSliceStart, timeSliceEnd)) { + continue; + } + + double hitEta = 0.0; + double hitPhi = 0.0; + etaPhiCalc(hit, hitEta, hitPhi); + + const auto [eta0, phi0] = binFunc(hitEta, hitPhi, 0); + const auto [eta1, phi1] = binFunc(hitEta, hitPhi, 1); + if (isValidEtaPhiBin(eta0, phi0)) { + grid[eta0][phi0]++; + gridTime[eta0][phi0] += hitT; + } + if (isValidEtaPhiBin(eta1, phi1)) { + gridShifted[eta1][phi1]++; + gridShiftedTime[eta1][phi1] += hitT; + } + } + } + + template + void fillEtaPhiGridsMatched(const CollectionT* collection, size_t& iniHitID, + double timeResolution, double timeSliceStart, double timeSliceEnd, + const EtaPhiGrid& baseGrid, const EtaPhiGrid& baseGridShifted, + EtaPhiGrid& compGrid, EtaPhiGrid& compGridShifted, int baseThreshold, + EtaPhiTimeGrid& compGridTime, EtaPhiTimeGrid& compGridShiftedTime, + BinFunc binFunc) { + if (collection == nullptr) { + return; + } + + const size_t hitCount = collection->size(); + for (size_t iHit = iniHitID; iHit < hitCount; ++iHit) { + const auto& hit = collection->at(iHit); + const double hitT = timeOfFlightCorrectedTime(hit); + if (hitT - timeResolution > timeSliceEnd) { + iniHitID = iHit; + break; + } + if (!judgeHitInTimeSlice(hitT, timeResolution, timeSliceStart, timeSliceEnd)) { + continue; + } + + double hitEta = 0.0; + double hitPhi = 0.0; + etaPhiCalc(collection->at(iHit), hitEta, hitPhi); + + const auto [eta0, phi0] = binFunc(hitEta, hitPhi, 0); + const auto [eta1, phi1] = binFunc(hitEta, hitPhi, 1); + if (isValidEtaPhiBin(eta0, phi0) && baseGrid[eta0][phi0] >= baseThreshold) { + compGrid[eta0][phi0]++; + compGridTime[eta0][phi0] += hitT; + } + if (isValidEtaPhiBin(eta1, phi1) && baseGridShifted[eta1][phi1] >= baseThreshold) { + compGridShifted[eta1][phi1]++; + compGridShiftedTime[eta1][phi1] += hitT; + } + } + } + + static size_t countGridCellsWithMultiplicity(const EtaPhiGrid& grid0, + const EtaPhiGrid& gridShifted, + const EtaPhiTimeGrid& gridTime0, + const EtaPhiTimeGrid& gridShiftedTime, int threshold, + double& averageTime); + + static double averageSelectedTriggerTime(const std::array& values, + const std::array& times, + std::initializer_list indices, + double fallbackTime); + + double trkTimeResolution(TrkCollectionIndex detectorID); + double calTimeResolution(CalCollectionIndex detectorID); + + template + TimeWindowSummary countHitsInTimeWindow(const CollectionT* collection, size_t startHitID, + double resolution, double window_start, + double window_end) { + TimeWindowSummary summary; + summary.nextStartID = startHitID; + if (collection == nullptr) { + return summary; + } + + for (size_t i = startHitID; i < collection->size(); ++i) { + const auto& hit = collection->at(i); + const double hitTime = timeOfFlightCorrectedTime(hit); + if (judgeOverTimeWindow(hitTime, resolution, window_end)) { + break; + } + + // Drop hits which cannot overlap the next window (which begins at window_end) + if (hitTime + resolution < window_end) { + summary.nextStartID = i + 1; + } + + if (overlapsTimeWindow(hitTime, resolution, window_start, window_end)) { + ++summary.count; + summary.timeSum += hitTime; + } + } + return summary; + } + + template + static void + copyTrkHitWithRelations(const edm4eic::TrackerHit& trackerHit, + const edm4eic::MCRecoTrackerHitAssociationCollection* associations, + const TrackerAssociationIndex& association_index, + TrackerHitOutputT& trackerHits_out, RawHitOutputT& rawHits_out, + AssociationOutputT& associations_out, + std::unique_ptr& simHits_out, + std::unique_ptr& links_out, + std::unique_ptr& mc_particles_out) { + + auto trackerHitCopied = trackerHit.clone(); + trackerHitCopied.setRawHit(edm4eic::RawTrackerHit()); + + if (associations == nullptr || !trackerHit.getRawHit().isAvailable()) { + trackerHits_out->push_back(trackerHitCopied); + return; + } + + const auto rawHitId = trackerHit.getRawHit().getObjectID(); + + auto rawHitCopied = trackerHit.getRawHit().clone(); + rawHits_out->push_back(rawHitCopied); + trackerHitCopied.setRawHit(rawHitCopied); + + const auto assocIterCal = association_index.find(objIdKey(rawHitId)); + + if (assocIterCal != association_index.end()) { + for (const size_t index : assocIterCal->second) { + const auto association = associations->at(index); + + if (!association.getSimHit().isAvailable()) { + continue; + } + + const auto simHit = association.getSimHit(); + auto simHitCopied = simHit.clone(false); + + if (simHit.getParticle().isAvailable()) { + const auto mcPIndex = simHit.getParticle().getObjectID(); + + if (mcPIndex.index >= 0 && + static_cast(mcPIndex.index) < mc_particles_out->size()) { + simHitCopied.setParticle((*mc_particles_out)[mcPIndex.index]); + } + } + + simHits_out->push_back(simHitCopied); + + auto copiedAsso = associations_out->create(); + copiedAsso.setWeight(association.getWeight()); + copiedAsso.setRawHit(rawHitCopied); + copiedAsso.setSimHit(simHitCopied); + + auto copiedLink = links_out->create(); + copiedLink.setWeight(association.getWeight()); + copiedLink.setFrom(rawHitCopied); + copiedLink.setTo(simHitCopied); + } + } + + trackerHits_out->push_back(trackerHitCopied); + } +}; diff --git a/src/global/splitting/splitting.cc b/src/global/splitting/splitting.cc new file mode 100644 index 0000000000..7354d8b867 --- /dev/null +++ b/src/global/splitting/splitting.cc @@ -0,0 +1,34 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2026 Takuya Kumaoka + +#include +#include +#include +#include + +#include "TimeframeSplitter.h" + +extern "C" { +void InitPlugin(JApplication* app) { + + InitJANAPlugin(app); + + const bool splitTimeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + if (!splitTimeframes) { + return; + } + + // Unfolder that takes timeframes and splits them into physics events. + eicrecon::JOmniUnfolderGeneratorT splitter_generator( + {.tag = "TimeframeSplitter", + .parent_level = JEventLevel::Timeslice, + .child_level = JEventLevel::PhysicsEvent, + //.input_names = {"EventHeader", "MCParticles"}, + //.variadic_input_names = {{}}, + //.output_names = {}, + //.variadic_output_names = {{}}, + .configs = {}}); + splitter_generator.Generate(app); +} +} // "C" diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index a6e960148d..d4a387d21f 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -525,7 +525,11 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "BarrelNeutralCandidateParticlesAlpha", "EndcapPNeutralCandidateParticlesAlpha", + "EventHeader_PHY", + "EventHeader_BKG", + }; + std::vector output_exclude_collections; // need to get as vector, then convert to set japp->SetDefaultParameter( "podio:output_collections", output_collections, diff --git a/src/services/io/podio/podio.cc b/src/services/io/podio/podio.cc index 47622f0609..766cb68fcc 100644 --- a/src/services/io/podio/podio.cc +++ b/src/services/io/podio/podio.cc @@ -3,9 +3,11 @@ // // +#include #include #include #include +#include #include #include "JEventProcessorManagedPODIO.h" @@ -13,17 +15,33 @@ #include "JEventSourceManagedPODIO.h" #include "JEventSourcePODIO.h" +namespace { + +template +class LeveledEventSourceGeneratorT : public JEventSourceGeneratorT { +public: + explicit LeveledEventSourceGeneratorT(JEventLevel level) { this->SetLevel(level); } +}; + +} // namespace + // Make this a JANA plugin extern "C" { void InitPlugin(JApplication* app) { InitJANAPlugin(app); + const bool split_timeframes = + app->RegisterParameter("split_timeframes", false, "Enable timeframe splitting"); + const auto source_level = split_timeframes ? JEventLevel::Timeslice : JEventLevel::PhysicsEvent; + // Check if managed mode is requested if (app->GetJParameterManager()->Exists("podio:managed_socket_path")) { - app->Add(new JEventSourceManagedPODIO("", app)); + auto* source = new JEventSourceManagedPODIO("", app); + source->SetLevel(source_level); + app->Add(source); app->Add(new JEventProcessorManagedPODIO()); } else { - app->Add(new JEventSourceGeneratorT()); + app->Add(new LeveledEventSourceGeneratorT(source_level)); app->Add(new JEventProcessorPODIO()); } } diff --git a/src/utilities/eicrecon/eicrecon.cc b/src/utilities/eicrecon/eicrecon.cc index 34a4757c84..05320613de 100644 --- a/src/utilities/eicrecon/eicrecon.cc +++ b/src/utilities/eicrecon/eicrecon.cc @@ -50,6 +50,9 @@ std::vector EICRECON_DEFAULT_PLUGINS = { "ECTOF", "LOWQ2", "LUMISPECCAL", +#if defined(JANA_VERSION_MAJOR) && JANA_VERSION_MAJOR >= 2026 + "splitting", +#endif "podio", "janatop", // clang-format on