diff --git a/.github/eicrecon.json b/.github/eicrecon.json new file mode 100644 index 0000000000..a641e8b2c9 --- /dev/null +++ b/.github/eicrecon.json @@ -0,0 +1,15 @@ +{ + "problemMatcher": [ + { + "owner": "eicrecon", + "pattern": [ + { + "regexp": "^\\[(.*)\\]\\s\\[(error)\\]\\s*(.*)$", + "file": 1, + "severity": 2, + "message": 3 + } + ] + } + ] +} diff --git a/.github/ubsan.supp b/.github/ubsan.supp index 94a63e2c4c..003670a034 100644 --- a/.github/ubsan.supp +++ b/.github/ubsan.supp @@ -1,3 +1,7 @@ implicit-integer-sign-change:JsonMaterialDecorator -implicit-integer-sign-change:/usr/local/include/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h +implicit-integer-sign-change:/opt/local/include/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h implicit-integer-sign-change:TObject +implicit-signed-integer-truncation:/opt/local/include/boost/iostreams/filter/gzip.hpp +implicit-integer-sign-change:/opt/local/include/eigen3/Eigen/src/Core/arch/SSE/Complex.h +implicit-integer-sign-change:/opt/local/include/root/TStorage.h +implicit-unsigned-integer-truncation:/opt/local/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp diff --git a/.github/workflows/linux-eic-shell.yml b/.github/workflows/linux-eic-shell.yml index d7c6df253a..c0d6d5a945 100644 --- a/.github/workflows/linux-eic-shell.yml +++ b/.github/workflows/linux-eic-shell.yml @@ -43,7 +43,7 @@ env: JANA_OPTIONS: "-Pjana:ticker_interval=60000 -Pjana:extended_report=1 -Pjana:warmup_timeout=0 -Pjana:timeout=0" ASAN_OPTIONS: suppressions=${{ github.workspace }}/.github/asan.supp:malloc_context_size=20:detect_leaks=1:verify_asan_link_order=0:detect_stack_use_after_return=1:detect_odr_violation=1:new_delete_type_mismatch=0:intercept_tls_get_addr=0 LSAN_OPTIONS: suppressions=${{ github.workspace }}/.github/lsan.supp - UBSAN_OPTIONS: suppressions=${{ github.workspace }}/.github/ubsan.supp:print_stacktrace=1:silence_unsigned_overflow=1 + UBSAN_OPTIONS: suppressions=${{ github.workspace }}/.github/ubsan.supp:print_stacktrace=1:silence_unsigned_overflow=1,report_error_type=1 jobs: env: @@ -78,15 +78,21 @@ jobs: CXX: clang++ # all clang debug jobs get extra flags - CC: clang + CXX: clang++ CMAKE_BUILD_TYPE: Debug CXXFLAGS: -fprofile-instr-generate -fcoverage-mapping # include some configs explicitly (if not yet included) - CC: clang + CXX: clang++ CMAKE_BUILD_TYPE: Release release: nightly - CC: clang + CXX: clang++ CMAKE_BUILD_TYPE: Release release: 24.04.0-stable + - CC: clang + CMAKE_BUILD_TYPE: Release + release: 24.05.2-stable steps: - name: mmap rnd_bits workaround # https://github.com/actions/runner/issues/3207#issuecomment-2000066889 @@ -304,24 +310,37 @@ jobs: name: codecov_report path: build/codecov_report/ + detector-info: + runs-on: ubuntu-latest + outputs: + hash: ${{ steps.detector-info.outputs.hash }} + steps: + - uses: cvmfs-contrib/github-action-cvmfs@v4 + - name: Get detector info + id: detector-info + run: | + file=/cvmfs/singularity.opensciencegrid.org/eicweb/${{ env.platform }}:${{ env.release }}/etc/jug_info + test -f ${file} + grep ' epic@' ${file} + hash=$(grep ' epic@' ${file} | sha256sum) + echo "hash=${hash%% *}" | tee $GITHUB_OUTPUT + npsim-gun: runs-on: ubuntu-latest + needs: + - detector-info strategy: matrix: particle: [pi, e] detector_config: [craterlake] steps: - uses: cvmfs-contrib/github-action-cvmfs@v4 - - name: Get detector info - id: detector_info - run: | - grep epic/nightly /cvmfs/singularity.opensciencegrid.org/eicweb/jug_xl\:nightly/etc/jug_info | sed 's/.*: .*-\(.*\)/hash=\1/g' >> $GITHUB_OUTPUT - name: Retrieve simulation files id: retrieve_simulation_files uses: actions/cache@v4 with: path: sim_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4hep.root - key: sim_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4hep.root-${{ steps.detector_info.outputs.hash }} + key: sim_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4hep.root-${{ needs.detector-info.outputs.hash }} - name: Produce simulation files uses: eic/run-cvmfs-osg-eic-shell@main if: steps.retrieve_simulation_files.outputs.cache-hit != 'true' @@ -338,6 +357,8 @@ jobs: npsim-gun-EcalLumiSpec: runs-on: ubuntu-latest + needs: + - detector-info strategy: matrix: particle: [e] @@ -347,16 +368,12 @@ jobs: with: filter: "tree:0" - uses: cvmfs-contrib/github-action-cvmfs@v4 - - name: Get detector info - id: detector_info - run: | - grep epic/nightly /cvmfs/singularity.opensciencegrid.org/eicweb/jug_xl\:nightly/etc/jug_info | sed 's/.*: .*-\(.*\)/hash=\1/g' >> $GITHUB_OUTPUT - name: Retrieve simulation files id: retrieve_simulation_files uses: actions/cache@v4 with: path: sim_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4hep.root - key: sim_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4hep.root-${{ steps.detector_info.outputs.hash }} + key: sim_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4hep.root-${{ needs.detector-info.outputs.hash }} - name: Produce simulation files uses: eic/run-cvmfs-osg-eic-shell@main if: steps.retrieve_simulation_files.outputs.cache-hit != 'true' @@ -374,23 +391,25 @@ jobs: npsim-dis: runs-on: ubuntu-latest + needs: + - detector-info strategy: matrix: - beam: [10x100, 18x275] - minq2: [1000] - detector_config: [craterlake] + include: + - beam: 10x100 + minq2: 1000 + detector_config: craterlake_tracking_only + - beam: 18x275 + minq2: 1000 + detector_config: craterlake steps: - uses: cvmfs-contrib/github-action-cvmfs@v4 - - name: Get detector info - id: detector_info - run: | - grep epic/nightly /cvmfs/singularity.opensciencegrid.org/eicweb/jug_xl\:nightly/etc/jug_info | sed 's/.*: .*-\(.*\)/hash=\1/g' >> $GITHUB_OUTPUT - name: Retrieve simulation files id: retrieve_simulation_files uses: actions/cache@v4 with: path: sim_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4hep.root - key: sim_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4hep.root-${{ steps.detector_info.outputs.hash }} + key: sim_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4hep.root-${{ needs.detector-info.outputs.hash }} - name: Produce simulation files uses: eic/run-cvmfs-osg-eic-shell@main if: steps.retrieve_simulation_files.outputs.cache-hit != 'true' @@ -408,22 +427,20 @@ jobs: npsim-minbias: runs-on: ubuntu-latest + needs: + - detector-info strategy: matrix: beam: [5x41, 10x100, 18x275] detector_config: [craterlake] steps: - uses: cvmfs-contrib/github-action-cvmfs@v4 - - name: Get detector info - id: detector_info - run: | - grep epic/nightly /cvmfs/singularity.opensciencegrid.org/eicweb/jug_xl\:nightly/etc/jug_info | sed 's/.*: .*-\(.*\)/hash=\1/g' >> $GITHUB_OUTPUT - name: Retrieve simulation files id: retrieve_simulation_files uses: actions/cache@v4 with: path: sim_dis_${{matrix.beam}}_minQ2=0_${{ matrix.detector_config }}.edm4hep.root - key: sim_dis_${{matrix.beam}}_minQ2=0_${{ matrix.detector_config }}.edm4hep.root-${{ steps.detector_info.outputs.hash }} + key: sim_dis_${{matrix.beam}}_minQ2=0_${{ matrix.detector_config }}.edm4hep.root-${{ needs.detector-info.outputs.hash }} - name: Produce simulation files uses: eic/run-cvmfs-osg-eic-shell@main if: steps.retrieve_simulation_files.outputs.cache-hit != 'true' @@ -479,7 +496,7 @@ jobs: 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${JANA_PLUGIN_PATH:+:${JANA_PLUGIN_PATH}} - $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_include_collections=EventHeader,MCParticles,EcalBarrelScFiRawHits,EcalBarrelImagingRawHits -Ppodio:output_file=two_stage_raw_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root sim_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4hep.root -Pplugins=dump_flags,janadot -Pdump_flags:json=${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json + $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_collections=EventHeader,MCParticles,EcalBarrelScFiRawHits,EcalBarrelImagingRawHits -Ppodio:output_file=two_stage_raw_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root sim_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4hep.root -Pplugins=dump_flags,janadot -Pdump_flags:json=${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json - name: Upload digitization output uses: actions/upload-artifact@v4 with: @@ -495,7 +512,7 @@ jobs: 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${JANA_PLUGIN_PATH:+:${JANA_PLUGIN_PATH}} - $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_include_collections=EcalBarrelClusters,EcalBarrelClusterAssociations -Ppodio:output_file=two_stage_rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root two_stage_raw_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root -Pplugins=dump_flags,janadot -Pdump_flags:json=${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json + $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_collections=EcalBarrelClusters,EcalBarrelClusterAssociations -Ppodio:output_file=two_stage_rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root two_stage_raw_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root -Pplugins=dump_flags,janadot -Pdump_flags:json=${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json - name: Upload reconstruction output uses: actions/upload-artifact@v4 with: @@ -691,18 +708,27 @@ jobs: 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 $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_file=rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root sim_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4hep.root -Pplugins=dump_flags,janadot,janatop $(<${{ github.workspace }}/.github/janadot.groups) -Pdump_flags:json=${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json - sed '/podio::Frame/d;/JEventProcessorPODIO/d' jana.dot | dot -Tsvg > jana.svg - mv jana.dot rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.dot - mv jana.svg rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.svg - uses: actions/upload-artifact@v4 with: name: rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root path: rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root if-no-files-found: error + - name: Convert execution graphs + uses: eic/run-cvmfs-osg-eic-shell@main + with: + platform-release: "${{ env.platform }}:${{ env.release }}" + setup: "/opt/detector/epic-${{ env.detector-version }}/bin/thisepic.sh" + run: | + mv jana.dot rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.dot + sed '/rank=sink/s/"podio::Frame";//; /podio::Frame/d; /rank=source/s/"JEventProcessorPODIO";//; /JEventProcessorPODIO/d' *.dot | dot -Tsvg > jana.svg + mv jana.svg rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.svg + continue-on-error: true - uses: actions/upload-artifact@v4 with: name: ${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json @@ -717,7 +743,7 @@ jobs: if-no-files-found: error - name: Download previous artifact id: download_previous_artifact - uses: dawidd6/action-download-artifact@v3 + uses: dawidd6/action-download-artifact@v6 with: branch: ${{ github.base_ref || github.event.merge_group.base_ref || github.ref_name }} path: ref/ @@ -734,7 +760,9 @@ jobs: run: | export PYTHONPATH=$HOME/.local/lib/python3.10/site-packages:$PYTHONPATH mkdir capybara-reports - capybara bara ref/rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root + mkdir new + ln -sf ../rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root new/ + capybara bara ref/rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root new/rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }}.edm4eic.root mv capybara-reports rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }} touch .rec_${{ matrix.particle }}_1GeV_20GeV_${{ matrix.detector_config }} - uses: actions/upload-artifact@v4 @@ -784,15 +812,22 @@ jobs: 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 - $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_file=rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4eic.root sim_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4hep.root -Ppodio:output_include_collections=EcalLumiSpecRawHits,EcalLumiSpecRecHits,EcalLumiSpecClusters,EcalLumiSpecClusterAssociations -PLUMISPECCAL:EcalLumiSpecIslandProtoClusters:splitCluster=1 -Pplugins=dump_flags,janadot,janatop $(<${{ github.workspace }}/.github/janadot.groups) -Pdump_flags:json=${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json - sed '/podio::Frame/d;/JEventProcessorPODIO/d' jana.dot | dot -Tsvg > jana.svg - mv jana.dot rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.dot - mv jana.svg rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.svg + $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_file=rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4eic.root sim_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4hep.root -Ppodio:output_collections=EcalLumiSpecRawHits,EcalLumiSpecRecHits,EcalLumiSpecClusters,EcalLumiSpecClusterAssociations -PLUMISPECCAL:EcalLumiSpecIslandProtoClusters:splitCluster=1 -Pplugins=dump_flags,janadot,janatop $(<${{ github.workspace }}/.github/janadot.groups) -Pdump_flags:json=${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json - uses: actions/upload-artifact@v4 with: name: rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4eic.root path: rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.edm4eic.root if-no-files-found: error + - name: Convert execution graphs + uses: eic/run-cvmfs-osg-eic-shell@main + with: + platform-release: "${{ env.platform }}:${{ env.release }}" + setup: "/opt/detector/epic-${{ env.detector-version }}/bin/thisepic.sh" + run: | + mv jana.dot rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.dot + sed '/rank=sink/s/"podio::Frame";//; /podio::Frame/d; /rank=source/s/"JEventProcessorPODIO";//; /JEventProcessorPODIO/d' *.dot | dot -Tsvg > jana.svg + mv jana.svg rec_${{ matrix.particle }}_EcalLumiSpec_${{ matrix.detector_config }}.svg + continue-on-error: true - uses: actions/upload-artifact@v4 with: name: ${{ matrix.particle }}_${{ matrix.detector_config }}_flags.json @@ -816,13 +851,21 @@ jobs: matrix: CC: [clang] beam: [10x100, 18x275] - minq2: [0, 1000] + minq2: [0] detector_config: [craterlake] include: - CC: clang beam: 5x41 minq2: 0 detector_config: craterlake + - CC: clang + beam: 10x100 + minq2: 1000 + detector_config: craterlake_tracking_only + - CC: clang + beam: 18x275 + minq2: 1000 + detector_config: craterlake steps: - name: mmap rnd_bits workaround # https://github.com/actions/runner/issues/3207#issuecomment-2000066889 @@ -849,18 +892,26 @@ jobs: 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 $PWD/install/bin/eicrecon $JANA_OPTIONS -Ppodio:output_file=rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root sim_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4hep.root -Pacts:WriteObj=true -Pacts:WritePly=true -Pplugins=janadot,janatop $(<${{ github.workspace }}/.github/janadot.groups) - sed '/podio::Frame/d;/JEventProcessorPODIO/d' jana.dot | dot -Tsvg > jana.svg - mv jana.dot rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.dot - mv jana.svg rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.svg - uses: actions/upload-artifact@v4 with: name: rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root path: rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root if-no-files-found: error + - name: Convert execution graphs + uses: eic/run-cvmfs-osg-eic-shell@main + with: + platform-release: "${{ env.platform }}:${{ env.release }}" + setup: "/opt/detector/epic-${{ env.detector-version }}/bin/thisepic.sh" + run: | + mv jana.dot rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.dot + sed '/rank=sink/s/"podio::Frame";//; /podio::Frame/d; /rank=source/s/"JEventProcessorPODIO";//; /JEventProcessorPODIO/d' *.dot | dot -Tsvg > jana.svg + mv jana.svg rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.svg + continue-on-error: true - uses: actions/upload-artifact@v4 with: name: rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.dot @@ -883,7 +934,7 @@ jobs: if-no-files-found: error - name: Download previous artifact id: download_previous_artifact - uses: dawidd6/action-download-artifact@v3 + uses: dawidd6/action-download-artifact@v6 with: branch: ${{ github.base_ref || github.event.merge_group.base_ref || github.ref_name }} path: ref/ @@ -901,7 +952,9 @@ jobs: pip install 'pygithub>=2' 'bokeh>=3' export PYTHONPATH=$HOME/.local/lib/python3.10/site-packages:$PYTHONPATH mkdir capybara-reports - capybara bara ref/rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root + mkdir new + ln -sf ../rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root new/ + capybara bara ref/rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root new/rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }}.edm4eic.root mv capybara-reports rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }} touch .rec_dis_${{matrix.beam}}_minQ2=${{matrix.minq2}}_${{ matrix.detector_config }} - uses: actions/upload-artifact@v4 @@ -915,7 +968,7 @@ jobs: trigger-container: runs-on: ubuntu-latest - if: ${{ github.event_name != 'merge_group' && github.actor != 'dependabot[bot]' }} + if: ${{ github.event_name != 'merge_group' && github.event_name != 'schedule' && github.actor != 'dependabot[bot]' }} needs: - eicrecon-gun steps: @@ -1013,7 +1066,7 @@ jobs: max-parallel: 4 steps: - name: Download docs artifact (other PRs) - uses: dawidd6/action-download-artifact@v3 + uses: dawidd6/action-download-artifact@v6 if: github.event.pull_request.number != matrix.pr with: commit: ${{ matrix.head_sha }} @@ -1030,7 +1083,7 @@ jobs: path: publishing_docs/pr/${{ matrix.pr }}/ - name: Download capybara artifact (other PRs) id: download_capybara - uses: dawidd6/action-download-artifact@v3 + uses: dawidd6/action-download-artifact@v6 if: github.event.pull_request.number != matrix.pr with: commit: ${{ matrix.head_sha }} @@ -1062,7 +1115,7 @@ jobs: # - If we run this on a non-main branch, we download from main with action-download-artifact. # - If we run this on the main branch, we have to download from this pipeline with download-artifact. - name: Download docs artifact (in PR) - uses: dawidd6/action-download-artifact@v3 + uses: dawidd6/action-download-artifact@v6 if: github.ref_name != 'main' with: branch: main @@ -1079,7 +1132,7 @@ jobs: path: publishing_docs/ - name: Download capybara artifact (in PR) id: download_capybara - uses: dawidd6/action-download-artifact@v3 + uses: dawidd6/action-download-artifact@v6 if: github.ref_name != 'main' with: commit: ${{ matrix.head_sha }} @@ -1117,7 +1170,9 @@ jobs: - get-docs-from-open-prs steps: - name: Merge GitHub Pages staging artifact - uses: actions/upload-artifact/merge@v4 + # FIXME pinned per https://github.com/actions/upload-artifact/issues/485#issuecomment-2271527517 + # to avoid https://github.com/eic/EICrecon/actions/runs/10486195490/job/29045183375#step:2:111 + uses: actions/upload-artifact/merge@v4.3.5 with: name: github-pages-staging pattern: | diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6fb4cf89e0..cf215b1c65 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -7,7 +7,7 @@ repos: - id: end-of-file-fixer - id: trailing-whitespace - repo: https://github.com/codespell-project/codespell - rev: v2.2.6 + rev: v2.3.0 hooks: - id: codespell - repo: https://github.com/Lucas-C/pre-commit-hooks diff --git a/CMakeLists.txt b/CMakeLists.txt index e09306c020..0d67602baf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,7 +31,7 @@ cmake_minimum_required(VERSION 3.24) # 3.27: find_package() uses upper-case _ROOT variables cmake_policy(SET CMP0144 NEW) -project(EICrecon) +project(EICrecon LANGUAGES CXX) # CMake includes include(CheckCXXCompilerFlag) @@ -197,8 +197,7 @@ endif() find_package(DD4hep REQUIRED) # ACTS cmake-lint: disable=C0103 -find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo - PluginDD4hep) +find_package(Acts REQUIRED COMPONENTS Core PluginDD4hep PluginJson) set(Acts_VERSION_MIN "19.0.0") set(Acts_VERSION "${Acts_VERSION_MAJOR}.${Acts_VERSION_MINOR}.${Acts_VERSION_PATCH}") @@ -233,6 +232,7 @@ endif() option(USE_ONNX "Compile with ONNX support" ON) if(${USE_ONNX}) find_package(onnxruntime CONFIG) + add_compile_definitions(USE_ONNX) endif() # Add CMake additional functionality: diff --git a/cmake/jana_plugin.cmake b/cmake/jana_plugin.cmake index ea5c0b2ffa..b2616f9167 100644 --- a/cmake/jana_plugin.cmake +++ b/cmake/jana_plugin.cmake @@ -269,8 +269,7 @@ endmacro() macro(plugin_add_acts _name) if(NOT Acts_FOUND) - find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo - PluginJson PluginDD4hep) + find_package(Acts REQUIRED COMPONENTS Core PluginDD4hep PluginJson) set(Acts_VERSION_MIN "20.2.0") set(Acts_VERSION "${Acts_VERSION_MAJOR}.${Acts_VERSION_MINOR}.${Acts_VERSION_PATCH}") @@ -291,10 +290,8 @@ macro(plugin_add_acts _name) plugin_link_libraries( ${PLUGIN_NAME} ActsCore - ActsPluginIdentification - ActsPluginTGeo - ActsPluginJson ActsPluginDD4hep + ActsPluginJson ${ActsCore_PATH}/${CMAKE_SHARED_LIBRARY_PREFIX}ActsExamplesFramework${CMAKE_SHARED_LIBRARY_SUFFIX} ) if(${_name}_WITH_LIBRARY) diff --git a/docs/howto/visualize_callgraph.md b/docs/howto/visualize_callgraph.md index 19f51ecbf3..d58c965248 100644 --- a/docs/howto/visualize_callgraph.md +++ b/docs/howto/visualize_callgraph.md @@ -48,12 +48,12 @@ By default `eicrecon` activates the full reconstruction. This will result in a very busy callgraph that can be hard to read if you are only interested in one detector. (See [here](https://eic.github.io/EICrecon/howto/callgraphs/all.png).) To activate for a specific factory, you can specify the relevant output -collection using the `podio:output_include_collections` parameter. +collection using the `podio:output_collections` parameter. Here is an example for the *EcalEndcapN*: ~~~bash eicrecon -Pplugins=janadot \ - -Ppodio:output_include_collections=EcalEndcapNMergedClusters \ + -Ppodio:output_collections=EcalEndcapNMergedClusters \ sim_file.edm4hep.root ~~~ @@ -61,7 +61,7 @@ The above will produce something like the following: ![EcalEndcapNMergedClusters](callgraphs/EcalEndcapNMergedClusters.png) In the plot, the oval at the top is the JEventProcessor making the -initial request for objects. (This is where the *podio:output_include_collections* +initial request for objects. (This is where the *podio:output_collections* parameter is used). The rectangles in the middle represent factories (i.e. algorithms). The green trapezoid at the bottom is the object read from the input file. This is representative of the flow JANA implements. Namely, the requests for diff --git a/docs/tutorial/02-work-environment.md b/docs/tutorial/02-work-environment.md index 940013d55e..72bd516218 100644 --- a/docs/tutorial/02-work-environment.md +++ b/docs/tutorial/02-work-environment.md @@ -137,11 +137,11 @@ eicrecon 2022-09-26_ncdis10x100_minq2-1_100ev.edm4hep.root ## Generating a podio output file To write reconstructed values to an output file, you need to tell *eicrecon* what to write. -There are several options available, but the mosrt useful one is *podio:output_include_collections*. +There are several options available, but the mosrt useful one is *podio:output_collections*. This is a comma separated list of colelctions to write to the output file. For example: ```console -eicrecon -Ppodio:output_include_collections=ReconstructedParticles 2022-09-26_ncdis10x100_minq2-1_100ev.edm4hep.root +eicrecon -Ppodio:output_collections=ReconstructedParticles 2022-09-26_ncdis10x100_minq2-1_100ev.edm4hep.root ``` To see a list of possible collections, run *eicrecon -L* . diff --git a/docs/tutorial/04-factory.md b/docs/tutorial/04-factory.md index 1c123d51ec..bd3edb3246 100644 --- a/docs/tutorial/04-factory.md +++ b/docs/tutorial/04-factory.md @@ -287,7 +287,7 @@ extern "C" { Finally, we go ahead and trigger the factory (remember, factories won't do anything unless activated by a JEventProcessor). You can open the ~~~ bash -eicrecon in.root -Ppodio:output_file=out.root -Ppodio:output_include_collections=EcalEndcapNIslandClusters -Pjana:nevents=10 +eicrecon in.root -Ppodio:output_file=out.root -Ppodio:output_collections=EcalEndcapNIslandClusters -Pjana:nevents=10 ~~~ diff --git a/docs/tutorial/05-contributing.md b/docs/tutorial/05-contributing.md index ffd2a3e48b..75d81902fd 100644 --- a/docs/tutorial/05-contributing.md +++ b/docs/tutorial/05-contributing.md @@ -52,7 +52,7 @@ More on the EIC contribution guide is in [this tutorial][eic-environment-tutoria ## Coding style -One can find coding style and other contributins policies at [CONTRIBUTING.md](https://github.com/eic/EICrecon/blob/main/CONTRIBUTING.md). It is yet to be finished but one can find current decisions on coding style there +One can find coding style and other contributions policies at [CONTRIBUTING.md](https://github.com/eic/EICrecon/blob/main/CONTRIBUTING.md). It is yet to be finished but one can find current decisions on coding style there ## References diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index 60716fc61e..495fc927ca 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -161,14 +161,15 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co // Used to optionally constrain the cluster eta to those of the contributing hits float minHitEta = std::numeric_limits::max(); float maxHitEta = std::numeric_limits::min(); - auto time = pcl.getHits()[0].getTime(); - auto timeError = pcl.getHits()[0].getTimeError(); + auto time = 0; + auto timeError = 0; for (unsigned i = 0; i < pcl.getHits().size(); ++i) { const auto& hit = pcl.getHits()[i]; const auto weight = pcl.getWeights()[i]; debug("hit energy = {} hit weight: {}", hit.getEnergy(), weight); auto energy = hit.getEnergy() * weight; totalE += energy; + time += (hit.getTime() - time) * energy / totalE; cl.addToHits(hit); cl.addToHitContributions(energy); const float eta = edm4hep::utils::eta(hit.getPosition()); @@ -228,12 +229,6 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co // Additional convenience variables - // best estimate on the cluster direction is the cluster position - // for simple 2D CoG clustering - cl.setIntrinsicTheta(edm4hep::utils::anglePolar(cl.getPosition())); - cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(cl.getPosition())); - // TODO errors - //_______________________________________ // Calculate cluster profile: // radius, @@ -248,12 +243,11 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co Eigen::Vector3f sum1_3D = Eigen::Vector3f::Zero(); Eigen::Vector2cf eigenValues_2D = Eigen::Vector2cf::Zero(); Eigen::Vector3cf eigenValues_3D = Eigen::Vector3cf::Zero(); - //the axis is the direction of the eigenvalue corresponding to the largest eigenvalue. - double axis_x=0, axis_y=0, axis_z=0; - if (cl.getNhits() > 1) { + // the axis is the direction of the eigenvalue corresponding to the largest eigenvalue. + edm4hep::Vector3f axis; + if (cl.getNhits() > 1) { for (const auto& hit : pcl.getHits()) { - float w = weightFunc(hit.getEnergy(), totalE, logWeightBase, 0); // theta, phi @@ -306,19 +300,15 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co return std::real(a) < std::real(b); } ); - auto axis = eigenvectors.col(std::distance( + auto axis_eigen = eigenvectors.col(std::distance( eigenValues_3D.begin(), max_eigenvalue_it )); - axis_x=axis(0,0).real(); - axis_y=axis(1,0).real(); - axis_z=axis(2,0).real(); - double norm=sqrt(axis_x*axis_x+axis_y*axis_y+axis_z*axis_z); - if (norm!=0){ - axis_x/=norm; - axis_y/=norm; - axis_z/=norm; - } + axis = { + axis_eigen(0,0).real(), + axis_eigen(1,0).real(), + axis_eigen(2,0).real(), + }; } } @@ -329,10 +319,22 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co cl.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z cluster width 1 cl.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z cluster width 2 cl.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z cluster width 3 - //last 3 shape parameters are the components of the axis direction - cl.addToShapeParameters( axis_x ); - cl.addToShapeParameters( axis_y ); - cl.addToShapeParameters( axis_z ); + + + double dot_product = cl.getPosition() * axis; + if (dot_product < 0) { + axis = -1 * axis; + } + + if (m_cfg.longitudinalShowerInfoAvailable) { + cl.setIntrinsicTheta(edm4hep::utils::anglePolar(axis)); + cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(axis)); + // TODO intrinsicDirectionError + } else { + cl.setIntrinsicTheta(NAN); + cl.setIntrinsicPhi(NAN); + } + return std::move(cl); } diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h index f6ea18e2e2..1b0624462e 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h @@ -23,6 +23,8 @@ namespace eicrecon { std::vector logWeightBaseCoeffs{}; double logWeightBase_Eref = 50 * dd4hep::MeV; + bool longitudinalShowerInfoAvailable = false; + // Constrain the cluster position eta to be within // the eta of the contributing hits. This is useful to avoid edge effects // for endcaps. diff --git a/src/algorithms/calorimetry/CalorimeterHitDigi.cc b/src/algorithms/calorimetry/CalorimeterHitDigi.cc index 943b5aa4a1..d0e3824f30 100644 --- a/src/algorithms/calorimetry/CalorimeterHitDigi.cc +++ b/src/algorithms/calorimetry/CalorimeterHitDigi.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -181,10 +182,13 @@ void CalorimeterHitDigi::process( std::pow(m_cfg.eRes[2] / (edep), 2) ) : 0; - double corrMeanScale_value = corrMeanScale(leading_hit); - double ped = m_cfg.pedMeanADC + m_gaussian(m_generator) * m_cfg.pedSigmaADC; - unsigned long long adc = std::llround(ped + edep * corrMeanScale_value * ( 1.0 + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC); - unsigned long long tdc = std::llround((time + m_gaussian(m_generator) * tRes) * stepTDC); + double corrMeanScale_value = corrMeanScale(leading_hit); + + double ped = m_cfg.pedMeanADC + m_gaussian(m_generator) * m_cfg.pedSigmaADC; + + // Note: both adc and tdc values must be positive numbers to avoid integer wraparound + unsigned long long adc = std::max(std::llround(ped + edep * corrMeanScale_value * (1.0 + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC), 0LL); + unsigned long long tdc = std::llround((time + m_gaussian(m_generator) * tRes) * stepTDC); if (edep> 1.e-3) trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t corrMeanScale: {}", edep, adc, time, m_cfg.capTime, tdc, corrMeanScale_value); rawhits->create( diff --git a/src/algorithms/calorimetry/CalorimeterHitReco.cc b/src/algorithms/calorimetry/CalorimeterHitReco.cc index 0457592ae0..e150061197 100644 --- a/src/algorithms/calorimetry/CalorimeterHitReco.cc +++ b/src/algorithms/calorimetry/CalorimeterHitReco.cc @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -172,6 +173,11 @@ void CalorimeterHitReco::process( continue; } + if (rh.getAmplitude() > m_cfg.capADC) { + error("Encountered hit with amplitude {} outside of ADC capacity {}", rh.getAmplitude(), m_cfg.capADC); + continue; + } + // get layer and sector ID const int lid = id_dec != nullptr && !m_cfg.layerField.empty() ? static_cast(id_dec->get(cellID, layer_idx)) : -1; @@ -265,6 +271,9 @@ void CalorimeterHitReco::process( const decltype(edm4eic::CalorimeterHitData::local) local_position(pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm); +#if EDM4EIC_VERSION_MAJOR >= 7 + auto recohit = +#endif recohits->create( rh.getCellID(), energy, @@ -276,6 +285,9 @@ void CalorimeterHitReco::process( sid, lid, local_position); +#if EDM4EIC_VERSION_MAJOR >= 7 + recohit.setRawHit(rh); +#endif } } diff --git a/src/algorithms/calorimetry/CalorimeterIslandCluster.cc b/src/algorithms/calorimetry/CalorimeterIslandCluster.cc index 53845fd947..06aa4b811c 100644 --- a/src/algorithms/calorimetry/CalorimeterIslandCluster.cc +++ b/src/algorithms/calorimetry/CalorimeterIslandCluster.cc @@ -7,32 +7,29 @@ #include #include -#include -#include +#include #include #include #include #include #include #include -#include #include -#include #include #include #include +#include #include #include #include "CalorimeterIslandCluster.h" #include "algorithms/calorimetry/CalorimeterIslandClusterConfig.h" +#include "services/evaluator/EvaluatorSvc.h" using namespace edm4eic; namespace eicrecon { -unsigned int CalorimeterIslandCluster::function_id = 0; - static double Phi_mpi_pi(double phi) { return std::remainder(phi, 2 * M_PI); } @@ -112,7 +109,7 @@ void CalorimeterIslandCluster::init() { return true; }; - std::map> uprops{ + std::vector>> uprops{ {"localDistXY", m_cfg.localDistXY}, {"localDistXZ", m_cfg.localDistXZ}, {"localDistYZ", m_cfg.localDistYZ}, @@ -122,50 +119,34 @@ void CalorimeterIslandCluster::init() { {"dimScaledLocalDistXY", m_cfg.dimScaledLocalDistXY} }; - bool method_found = false; + auto& serviceSvc = algorithms::ServiceSvc::instance(); - // Adjacency matrix methods - if (!m_cfg.adjacencyMatrix.empty()) { - // sanity checks - if (m_cfg.readout.empty()) { - error("readoutClass is not provided, it is needed to know the fields in readout ids"); - } - m_idSpec = m_detector->readout(m_cfg.readout).idSpec(); - - std::string func_name = fmt::format("_CalorimeterIslandCluster_{}", function_id++); - std::ostringstream sstr; - sstr << "bool " << func_name << "(double params[]){"; - unsigned int param_ix = 0; + std::function hit_pair_to_map = [this](const edm4eic::CalorimeterHit &h1, const edm4eic::CalorimeterHit &h2) { + std::unordered_map params; for(const auto &p : m_idSpec.fields()) { const std::string &name = p.first; const dd4hep::IDDescriptor::Field* field = p.second; - sstr << "double " << name << "_1 = params[" << (param_ix++) << "];"; - sstr << "double " << name << "_2 = params[" << (param_ix++) << "];"; + params.emplace(name + "_1", field->value(h1.getCellID())); + params.emplace(name + "_2", field->value(h2.getCellID())); + trace("{}_1 = {}", name, field->value(h1.getCellID())); + trace("{}_2 = {}", name, field->value(h2.getCellID())); } - sstr << "return " << m_cfg.adjacencyMatrix << ";"; - sstr << "}"; - debug("Compiling {}", sstr.str()); - - TInterpreter *interp = TInterpreter::Instance(); - interp->ProcessLine(sstr.str().c_str()); - std::unique_ptr func_val { gInterpreter->MakeInterpreterValue() }; - interp->Evaluate(func_name.c_str(), *func_val); - typedef bool (*func_t)(double params[]); - func_t func = ((func_t)(func_val->GetAsPointer())); - - is_neighbour = [this, func, param_ix](const CaloHit &h1, const CaloHit &h2) { - std::vector params; - params.reserve(param_ix); - for(const auto &p : m_idSpec.fields()) { - const std::string &name = p.first; - const dd4hep::IDDescriptor::Field* field = p.second; - params.push_back(field->value(h1.getCellID())); - params.push_back(field->value(h2.getCellID())); - trace("{}_1 = {}", name, field->value(h1.getCellID())); - trace("{}_2 = {}", name, field->value(h2.getCellID())); - } - return func(params.data()); - }; + return params; + }; + + if (m_cfg.readout.empty()) { + if ((!m_cfg.adjacencyMatrix.empty()) || (!m_cfg.peakNeighbourhoodMatrix.empty())) { + throw std::runtime_error("'readout' is not provided, it is needed to know the fields in readout ids"); + } + } else { + m_idSpec = m_detector->readout(m_cfg.readout).idSpec(); + } + + bool method_found = false; + + // Adjacency matrix methods + if (!m_cfg.adjacencyMatrix.empty()) { + is_neighbour = serviceSvc.service("EvaluatorSvc")->compile(m_cfg.adjacencyMatrix, hit_pair_to_map); method_found = true; } @@ -188,7 +169,6 @@ void CalorimeterIslandCluster::init() { } }; - info("Using clustering method: {}", uprop.first); break; } } @@ -199,6 +179,12 @@ void CalorimeterIslandCluster::init() { } if (m_cfg.splitCluster) { + if (!m_cfg.peakNeighbourhoodMatrix.empty()) { + is_maximum_neighbourhood = serviceSvc.service("EvaluatorSvc")->compile(m_cfg.peakNeighbourhoodMatrix, hit_pair_to_map); + } else { + is_maximum_neighbourhood = is_neighbour; + } + auto transverseEnergyProfileMetric_it = std::find_if(distMethods.begin(), distMethods.end(), [&](auto &p) { return m_cfg.transverseEnergyProfileMetric == p.first; }); if (transverseEnergyProfileMetric_it == distMethods.end()) { throw std::runtime_error(fmt::format("Unsupported value \"{}\" for \"transverseEnergyProfileMetric\"", m_cfg.transverseEnergyProfileMetric)); diff --git a/src/algorithms/calorimetry/CalorimeterIslandCluster.h b/src/algorithms/calorimetry/CalorimeterIslandCluster.h index 2ef8620945..f70d469f0a 100644 --- a/src/algorithms/calorimetry/CalorimeterIslandCluster.h +++ b/src/algorithms/calorimetry/CalorimeterIslandCluster.h @@ -67,6 +67,9 @@ namespace eicrecon { // helper function to group hits std::function is_neighbour; + // helper function to define hit maximum + std::function is_maximum_neighbourhood; + // unitless counterparts of the input parameters std::array neighbourDist; @@ -75,13 +78,11 @@ namespace eicrecon { private: - static unsigned int function_id; - // grouping function with Breadth-First Search void bfs_group(const edm4eic::CalorimeterHitCollection &hits, std::set &group, std::size_t idx, std::vector &visits) const { visits[idx] = true; - // not a qualified hit to particpate clustering, stop here + // not a qualified hit to participate clustering, stop here if (hits[idx].getEnergy() < m_cfg.minClusterHitEdep) { return; } @@ -94,7 +95,7 @@ namespace eicrecon { for (std::size_t idx1 : group) { // check neighbours for (std::size_t idx2 = 0; idx2 < hits.size(); ++idx2) { - // not a qualified hit to particpate clustering, skip + // not a qualified hit to participate clustering, skip if (hits[idx2].getEnergy() < m_cfg.minClusterHitEdep) { continue; } @@ -140,7 +141,7 @@ namespace eicrecon { continue; } - if (is_neighbour(hits[idx1], hits[idx2]) && (hits[idx2].getEnergy() > hits[idx1].getEnergy())) { + if (is_maximum_neighbourhood(hits[idx1], hits[idx2]) && (hits[idx2].getEnergy() > hits[idx1].getEnergy())) { maximum = false; break; } diff --git a/src/algorithms/calorimetry/CalorimeterIslandClusterConfig.h b/src/algorithms/calorimetry/CalorimeterIslandClusterConfig.h index 45dcb567bc..eae2433769 100644 --- a/src/algorithms/calorimetry/CalorimeterIslandClusterConfig.h +++ b/src/algorithms/calorimetry/CalorimeterIslandClusterConfig.h @@ -10,6 +10,7 @@ namespace eicrecon { struct CalorimeterIslandClusterConfig { std::string adjacencyMatrix; + std::string peakNeighbourhoodMatrix; std::string readout; // neighbour checking distances diff --git a/src/algorithms/calorimetry/ImagingTopoCluster.h b/src/algorithms/calorimetry/ImagingTopoCluster.h index 0adde19641..7e086ecc1a 100644 --- a/src/algorithms/calorimetry/ImagingTopoCluster.h +++ b/src/algorithms/calorimetry/ImagingTopoCluster.h @@ -1,17 +1,27 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Chao Peng, Sylvester Joosten, Whitney Armstrong +// Copyright (C) 2022 Chao Peng, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck /* * Topological Cell Clustering Algorithm for Imaging Calorimetry * 1. group all the adjacent pixels * * Author: Chao Peng (ANL), 06/02/2021 - * References: https://arxiv.org/pdf/1603.02934.pdf + * Original reference: https://arxiv.org/pdf/1603.02934.pdf + * + * Modifications: + * + * Wouter Deconinck (Manitoba), 08/24/2024 + * - converted hit storage model from std::vector to std::set sorted on layer + * where only hits remaining to be assigned to a group are in the set + * - erase hits that are too low in energy to be part of a cluster + * - converted group storage model frmo std::set to std::list to allow adding + * hits while keeping iterators valid * */ #pragma once #include +#include #include #include @@ -40,16 +50,6 @@ namespace eicrecon { > >; - /** Topological Cell Clustering Algorithm. - * - * Topological Cell Clustering Algorithm for Imaging Calorimetry - * 1. group all the adjacent pixels - * - * Author: Chao Peng (ANL), 06/02/2021 - * References: https://arxiv.org/pdf/1603.02934.pdf - * - * \ingroup reco - */ class ImagingTopoCluster : public ImagingTopoClusterAlgorithm, public WithPodConfig { @@ -64,8 +64,13 @@ namespace eicrecon { private: // unitless counterparts of the input parameters - double localDistXY[2]{0, 0}, layerDistEtaPhi[2]{0, 0}, sectorDist{0}; - double minClusterHitEdep{0}, minClusterCenterEdep{0}, minClusterEdep{0}, minClusterNhits{0}; + std::array localDistXY{0, 0}; + std::array layerDistEtaPhi{0, 0}; + std::array layerDistXY{0, 0}; + double sectorDist{0}; + double minClusterHitEdep{0}; + double minClusterCenterEdep{0}; + double minClusterEdep{0}; public: void init() { @@ -79,10 +84,16 @@ namespace eicrecon { error( "Expected 2 values (eta_dist, phi_dist) for layerDistEtaPhi" ); return; } + if (m_cfg.minClusterCenterEdep < m_cfg.minClusterHitEdep) { + error( "minClusterCenterEdep must be greater than or equal to minClusterHitEdep" ); + return; + } // using juggler internal units (GeV, dd4hep::mm, dd4hep::ns, dd4hep::rad) localDistXY[0] = m_cfg.localDistXY[0] / dd4hep::mm; localDistXY[1] = m_cfg.localDistXY[1] / dd4hep::mm; + layerDistXY[0] = m_cfg.layerDistXY[0] / dd4hep::mm; + layerDistXY[1] = m_cfg.layerDistXY[1] / dd4hep::mm; layerDistEtaPhi[0] = m_cfg.layerDistEtaPhi[0]; layerDistEtaPhi[1] = m_cfg.layerDistEtaPhi[1] / dd4hep::rad; sectorDist = m_cfg.sectorDist / dd4hep::mm; @@ -95,10 +106,22 @@ namespace eicrecon { "Local [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", localDistXY[0], localDistXY[1] ); - info("Neighbour layers clustering (same sector and layer id within +- {:d}: " - "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].", - m_cfg.neighbourLayersRange, layerDistEtaPhi[0], layerDistEtaPhi[1] - ); + switch (m_cfg.layerMode) { + case ImagingTopoClusterConfig::ELayerMode::etaphi: + info("Neighbour layers clustering (same sector and layer id within +- {:d}: " + "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].", + m_cfg.neighbourLayersRange, layerDistEtaPhi[0], layerDistEtaPhi[1] + ); + break; + case ImagingTopoClusterConfig::ELayerMode::xy: + info("Neighbour layers clustering (same sector and layer id within +- {:d}: " + "Local [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", + m_cfg.neighbourLayersRange, layerDistXY[0], layerDistXY[1] + ); + break; + default: + error("Unknown layer mode."); + } info("Neighbour sectors clustering (different sector): " "Global distance between hits <= {:.4f} mm.", sectorDist @@ -110,21 +133,55 @@ namespace eicrecon { const auto [hits] = input; auto [proto] = output; - // group neighbouring hits - std::vector visits(hits->size(), false); - std::vector> groups; - for (size_t i = 0; i < hits->size(); ++i) { - debug("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})", i + 1, - (*hits)[i].getLocal().x, (*hits)[i].getLocal().y, (*hits)[i].getPosition().z, - (*hits)[i].getPosition().x, (*hits)[i].getPosition().y, (*hits)[i].getPosition().z + // Sort hit indices (podio collections do not support std::sort) + auto compare = [&hits](const auto& a, const auto& b) { + // if !(a < b) and !(b < a), then a and b are equivalent + // and only one of them will be allowed in a set + if ((*hits)[a].getLayer() == (*hits)[b].getLayer()) { + return (*hits)[a].getObjectID().index < (*hits)[b].getObjectID().index; + } + return (*hits)[a].getLayer() < (*hits)[b].getLayer(); + }; + // indices contains the remaining hit indices that have not + // been assigned to a group yet + std::set indices(compare); + // set does not have a size yet, so cannot fill with iota + for (std::size_t i = 0; i < hits->size(); ++i) { + indices.insert(i); + } + // ensure no hits were dropped due to equivalency in set + if (hits->size() != indices.size()) { + error("equivalent hits were dropped: #hits {:d}, #indices {:d}", hits->size(), indices.size()); + } + + // Group neighbouring hits + std::vector> groups; + // because indices changes, the loop over indices requires some care: + // - we must use iterators instead of range-for + // - erase returns an incremented iterator and therefore acts as idx++ + // - when the set becomes empty on erase, idx is invalid and idx++ will be too + // (also applies to loop in bfs_group below) + for (auto idx = indices.begin(); idx != indices.end(); + indices.empty() ? idx = indices.end() : idx) { + + debug("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {}), energy = {}", *idx, + (*hits)[*idx].getLocal().x, (*hits)[*idx].getLocal().y, (*hits)[*idx].getPosition().z, + (*hits)[*idx].getPosition().x, (*hits)[*idx].getPosition().y, (*hits)[*idx].getPosition().z, + (*hits)[*idx].getEnergy() ); - // already in a group, or not energetic enough to form a cluster - if (visits[i] || (*hits)[i].getEnergy() < minClusterCenterEdep) { - continue; + + // not energetic enough for cluster center, but could still be cluster hit + if ((*hits)[*idx].getEnergy() < minClusterCenterEdep) { + idx++; + continue; } + // create a new group, and group all the neighbouring hits - groups.emplace_back(); - bfs_group(*hits, groups.back(), i, visits); + groups.emplace_back(std::list{*idx}); + bfs_group(*hits, indices, groups.back(), *idx); + + // wait with erasing until after bfs_group to ensure iterator is not invalidated in bfs_group + idx = indices.erase(idx); // takes role of idx++ } debug("found {} potential clusters (groups of hits)", groups.size()); for (size_t i = 0; i < groups.size(); ++i) { @@ -133,7 +190,7 @@ namespace eicrecon { // form clusters for (const auto &group : groups) { - if (static_cast(group.size()) < m_cfg.minClusterNhits) { + if (group.size() < m_cfg.minClusterNhits) { continue; } double energy = 0.; @@ -169,45 +226,59 @@ namespace eicrecon { return (std::abs(h1.getLocal().x - h2.getLocal().x) <= localDistXY[0]) && (std::abs(h1.getLocal().y - h2.getLocal().y) <= localDistXY[1]); } else if (ldiff <= m_cfg.neighbourLayersRange) { + switch(m_cfg.layerMode){ + case eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi: return (std::abs(edm4hep::utils::eta(h1.getPosition()) - edm4hep::utils::eta(h2.getPosition())) <= layerDistEtaPhi[0]) && (std::abs(edm4hep::utils::angleAzimuthal(h1.getPosition()) - edm4hep::utils::angleAzimuthal(h2.getPosition())) <= layerDistEtaPhi[1]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::xy: + return (std::abs(h1.getPosition().x - h2.getPosition().x) <= layerDistXY[0]) && + (std::abs(h1.getPosition().y - h2.getPosition().y) <= layerDistXY[1]); + } } - // not in adjacent layers return false; } // grouping function with Breadth-First Search - void bfs_group(const edm4eic::CalorimeterHitCollection &hits, std::set &group, std::size_t idx, std::vector &visits) const { - visits[idx] = true; + // note: template to allow Compare only known in local scope of caller + template + void bfs_group(const edm4eic::CalorimeterHitCollection &hits, std::set& indices, std::list &group, const std::size_t idx) const { - // not a qualified hit to particpate clustering, stop here - if (hits[idx].getEnergy() < m_cfg.minClusterHitEdep) { - return; - } + // loop over group as it grows, until the end is stable and we reach it + for (auto idx1 = group.begin(); idx1 != group.end(); ++idx1) { + // check neighbours (note comments on loop over set above) + for (auto idx2 = indices.begin(); idx2 != indices.end(); + indices.empty() ? idx2 = indices.end() : idx2) { - group.insert(idx); - size_t prev_size = 0; + // skip idx1 and original idx + // (we cannot erase idx since it would invalidate iterator in calling scope) + if (*idx2 == *idx1 || *idx2 == idx) { + idx2++; + continue; + } - while (prev_size != group.size()) { - prev_size = group.size(); - for (std::size_t idx1 : group) { - // check neighbours - for (std::size_t idx2 = 0; idx2 < hits.size(); ++idx2) { - // not a qualified hit to particpate clustering, skip - if (hits[idx2].getEnergy() < m_cfg.minClusterHitEdep) { - continue; - } - if ((!visits[idx2]) - && is_neighbour(hits[idx1], hits[idx2])) { - group.insert(idx2); - visits[idx2] = true; - } + // skip rest of list of hits when we're past relevant layers + //if (hits[*idx2].getLayer() - hits[*idx1].getLayer() > m_cfg.neighbourLayersRange) { + // break; + //} + + // not energetic enough for cluster hit + if (hits[*idx2].getEnergy() < m_cfg.minClusterHitEdep) { + idx2 = indices.erase(idx2); + continue; + } + + if (is_neighbour(hits[*idx1], hits[*idx2])) { + group.push_back(*idx2); + idx2 = indices.erase(idx2); // takes role of idx2++ + } else { + idx2++; } } } } + }; } // namespace eicrecon diff --git a/src/algorithms/calorimetry/ImagingTopoClusterConfig.h b/src/algorithms/calorimetry/ImagingTopoClusterConfig.h index 26b80dbe2a..651a0f4619 100644 --- a/src/algorithms/calorimetry/ImagingTopoClusterConfig.h +++ b/src/algorithms/calorimetry/ImagingTopoClusterConfig.h @@ -4,6 +4,7 @@ #pragma once #include +#include namespace eicrecon { @@ -13,8 +14,13 @@ namespace eicrecon { int neighbourLayersRange = 1; // maximum distance of local (x, y) to be considered as neighbors at the same layer std::vector localDistXY = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; - // maximum distance of global (eta, phi) to be considered as neighbors at different layers + // maximum distance of global (eta, phi) to be considered as neighbors at different layers (if layerMode==etaphi) std::vector layerDistEtaPhi = {0.01, 0.01}; + // maximum distance of global (x, y) to be considered as neighbors at different layers (if layerMode==xy) + std::vector layerDistXY = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; + // determines how neighbors are determined for hits in different layers (using either eta and phi, or x and y) + enum ELayerMode {etaphi=0, xy=1} layerMode = etaphi; + // maximum global distance to be considered as neighbors in different sectors double sectorDist = 1.0 * dd4hep::cm; @@ -25,8 +31,35 @@ namespace eicrecon { // minimum cluster energy (to save this cluster) double minClusterEdep = 0.5 * dd4hep::MeV; // minimum number of hits (to save this cluster) - int minClusterNhits = 10; + std::size_t minClusterNhits = 10; }; + std::istream& operator>>(std::istream& in, ImagingTopoClusterConfig::ELayerMode& layerMode) { + std::string s; + in >> s; + // stringifying the enums causes them to be converted to integers before conversion to strings + if (s == "etaphi" or s=="0") { + layerMode = ImagingTopoClusterConfig::ELayerMode::etaphi; + } else if (s == "xy" or s=="1") { + layerMode = ImagingTopoClusterConfig::ELayerMode::xy; + } else { + in.setstate(std::ios::failbit); // Set the fail bit if the input is not valid + } + + return in; + } + std::ostream& operator<<(std::ostream& out, ImagingTopoClusterConfig::ELayerMode& layerMode) { + switch(layerMode) { + case ImagingTopoClusterConfig::ELayerMode::etaphi: + out << "etaphi"; + break; + case ImagingTopoClusterConfig::ELayerMode::xy: + out << "xy"; + break; + default: + out.setstate(std::ios::failbit); + } + return out; + } } // namespace eicrecon diff --git a/src/algorithms/calorimetry/TruthEnergyPositionClusterMerger.h b/src/algorithms/calorimetry/TruthEnergyPositionClusterMerger.h index 454cb00bd0..d41293d27c 100644 --- a/src/algorithms/calorimetry/TruthEnergyPositionClusterMerger.h +++ b/src/algorithms/calorimetry/TruthEnergyPositionClusterMerger.h @@ -8,9 +8,10 @@ #include // Event Model related classes -#include -#include -#include +#include +#include +#include +#include #include #include "algorithms/interfaces/WithPodConfig.h" diff --git a/src/algorithms/interfaces/CMakeLists.txt b/src/algorithms/interfaces/CMakeLists.txt index f481ebe153..74693faaa7 100644 --- a/src/algorithms/interfaces/CMakeLists.txt +++ b/src/algorithms/interfaces/CMakeLists.txt @@ -1,2 +1,10 @@ set(PLUGIN_NAME "algorithms_interfaces") -plugin_headers_only(${PLUGIN_NAME}) + +# Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets +# Setting default includes, libraries and installation paths +plugin_add(${PLUGIN_NAME} WITH_SHARED_LIBRARY WITHOUT_PLUGIN) + +# 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}) diff --git a/src/algorithms/interfaces/ParticleSvc.cc b/src/algorithms/interfaces/ParticleSvc.cc new file mode 100644 index 0000000000..c23128ded3 --- /dev/null +++ b/src/algorithms/interfaces/ParticleSvc.cc @@ -0,0 +1,252 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022, 2024 Sylvester Joosten, Wouter Deconinck + +#include + +#include "algorithms/interfaces/ParticleSvc.h" + +namespace algorithms { + +const std::shared_ptr ParticleSvc::kParticleMap = + std::make_shared(ParticleSvc::ParticleMap{ + { 0, { 0, 0, 0.0 , "unknown" }}, + { 11, { 11, -1, 0.000510998928, "e-" }}, + { -11, { -11, 1, 0.000510998928, "e+" }}, + { 13, { 13, -1, 0.105658357 , "mu-" }}, + { -13, { -13, 1, 0.105658357 , "mu+" }}, + { 22, { 22, 0, 0.0 , "gamma" }}, + { 111, { 111, 0, 0.1349766 , "p0" }}, + { 113, { 111, 0, 0.76850 , "rho(770)^0" }}, + { 115, { 115, 0, 1.31800 , "a_2(1320)^0" }}, + { 211, { 211, 1, 0.1395701 , "pi+" }}, + { -211, { -211, -1, 0.1395701 , "pi-" }}, + { 213, { 213, 1, 0.76690 , "rho+" }}, + { -213, { -213, -1, 0.76690 , "rho-" }}, + { 215, { 215, 1, 1.31800 , "a_2(1320)+" }}, + { -215, { -215, -1, 1.31800 , "a_2(1320)-" }}, + { 221, { 221, 0, 0.54745 , "eta" }}, + { 223, { 223, 0, 0.78194 , "omega" }}, + { 225, { 225, 0, 1.27500 , "f_2(1270)" }}, + { 130, { 130, 0, 0.49767 , "KL_0" }}, + { 310, { 310, 0, 0.49767 , "KS_0" }}, + { 311, { 311, 0, 0.49767 , "K^0" }}, + { -311, { -311, 0, 0.49767 , "K~^0" }}, + { 313, { 313, 0, 0.89610 , "K*(892)^0" }}, + { -313, { -313, 0, 0.89610 , "K*(892)~^0" }}, + { 315, { 315, 0, 1.43200 , "K*_2(1430)^0" }}, + { -315, { -315, 0, 1.43200 , "K*_2(1430)~^0" }}, + { 321, { 321, 1, 0.49360 , "K^+" }}, + { -321, { -321, -1, 0.49360 , "K^-" }}, + { 323, { 323, 1, 0.89160 , "K*(892)^+" }}, + { -323, { -323, -1, 0.89160 , "K*(892)^-" }}, + { 325, { 325, 1, 1.42500 , "K*_2(1430)^+" }}, + { -325, { -325, -1, 1.42500 , "K*_2(1430)^-" }}, + { 331, { 331, 0, 0.95777 , "eta'(958)" }}, + { 333, { 333, 0, 1.01940 , "phi(1020)" }}, + { 335, { 335, 0, 1.52500 , "f'_2(1525)" }}, + { 411, { 411, 1, 1.86930 , "D+" }}, + { -411, { 411, -1, 1.86930 , "D-" }}, + { 413, { 413, 1, 2.01000 , "D*(2010)^+" }}, + { -413, { -413, -1, 2.01000 , "D*(2010)^-" }}, + { 415, { 415, 1, 2.46000 , "D*_2(2460)^+" }}, + { -415, { -415, -1, 2.46000 , "D*_2(2460)^-" }}, + { 421, { 421, 0, 1.86450 , "D^0" }}, + { -421, { -421, 0, 1.86450 , "D~^0" }}, + { 423, { 423, 0, 2.00670 , "D*(2007)^0" }}, + { -423, { -423, 0, 2.00670 , "D*(2007)~^0" }}, + { 425, { 425, 0, 2.46000 , "D*_2(2460)^0" }}, + { -425, { -425, 0, 2.46000 , "D*_2(2460)~^0" }}, + { 431, { 431, 1, 1.96850 , "D_s^+" }}, + { -431, { -431, -1, 1.96850 , "D_s^-" }}, + { 433, { 433, 1, 2.11240 , "D*_s^+" }}, + { -433, { -433, -1, 2.11240 , "D*_s^-" }}, + { 435, { 435, 1, 2.57350 , "D*_s2(2573)^+" }}, + { -435, { -435, -1, 2.57350 , "D*_s2(2573)^-" }}, + { 441, { 441, 0, 2.97980 , "eta_c(1S)" }}, + { 443, { 443, 0, 3.09688 , "J/psi(1S)" }}, + { 445, { 445, 0, 3.55620 , "chi_c2(1P)" }}, + { 511, { 511, 0, 5.27920 , "B^0" }}, + { -511, { -511, 0, 5.27920 , "B~^0" }}, + { 513, { 513, 0, 5.32480 , "B*^0" }}, + { -513, { -513, 0, 5.32480 , "B*~^0" }}, + { 515, { 515, 0, 5.83000 , "B*_2^0" }}, + { -515, { -515, 0, 5.83000 , "B*_2~^0" }}, + { 521, { 521, 1, 5.27890 , "B^+" }}, + { -521, { -521, -1, 5.27890 , "B^-" }}, + { 523, { 523, 1, 5.32480 , "B*^+" }}, + { -523, { -523, -1, 5.32480 , "B*^-" }}, + { 525, { 525, 1, 5.83000 , "B*_2^+" }}, + { -525, { -525, -1, 5.83000 , "B*_2^-" }}, + { 531, { 531, 0, 5.36930 , "B_s^0" }}, + { -531, { -531, 0, 5.36930 , "B_s~^0" }}, + { 533, { 533, 0, 5.41630 , "B*_s^0" }}, + { -533, { -533, 0, 5.41630 , "B*_s~^0" }}, + { 535, { 535, 0, 6.07000 , "B*_s2^0" }}, + { -535, { -535, 0, 6.07000 , "B*_s2~^0" }}, + { 541, { 541, 1, 6.59400 , "B_c^+" }}, + { -541, { -541, -1, 6.59400 , "B_c^-" }}, + { 543, { 543, 1, 6.60200 , "B*_c^+" }}, + { -543, { -543, -1, 6.60200 , "B*_c^-" }}, + { 545, { 545, 1, 7.35000 , "B*_c2^+" }}, + { -545, { -545, -1, 7.35000 , "B*_c2^-" }}, + { 551, { 551, 0, 9.40000 , "eta_b(1S)" }}, + { 553, { 553, 0, 9.46030 , "Upsilon(1S)" }}, + { 555, { 555, 0, 9.91320 , "chi_b2(1P)" }}, + { 990, { 990, 0, 0.00000 , "pomeron" }}, + { 1114, { 1114, -1, 1.23400 , "Delta^-" }}, + { -1114, { -1114, 1, 1.23400 , "Delta~^+" }}, + { 2112, { 2112, 0, 0.93957 , "n" }}, + { -2112, { -2112, 0, 0.93957 , "n~^0" }}, + { 2114, { 2114, 0, 1.23300 , "Delta^0" }}, + { -2114, { -2114, 0, 1.23300 , "Delta~^0" }}, + { 2212, { 2212, 1, 0.93827 , "p^+" }}, + { -2212, { -2212, -1, 0.93827 , "p~^-" }}, + { 2214, { 2214, 1, 1.23200 , "Delta^+" }}, + { -2214, { -2214, -1, 1.23200 , "Delta~^-" }}, + { 2224, { 2224, 2, 1.23100 , "Delta^++" }}, + { -2224, { -2224, -2, 1.23100 , "Delta~^--" }}, + { 3112, { 3112, -1, 1.19744 , "Sigma^-" }}, + { -3112, { -3112, 1, 1.19744 , "Sigma~^+" }}, + { 3114, { 3114, -1, 1.38720 , "Sigma*^-" }}, + { -3114, { -3114, 1, 1.38720 , "Sigma*~^+" }}, + { 3122, { 3122, 0, 1.11568 , "Lambda^0" }}, + { -3122, { -3122, 0, 1.11568 , "Lambda~^0" }}, + { 3212, { 3212, 0, 1.19255 , "Sigma^0" }}, + { -3212, { -3212, 0, 1.19255 , "Sigma~^0" }}, + { 3214, { 3214, 0, 1.38370 , "Sigma*^0" }}, + { -3214, { -3214, 0, 1.38370 , "Sigma*~^0" }}, + { 3222, { 3222, 1, 1.18937 , "Sigma^+" }}, + { -3222, { -3222, -1, 1.18937 , "Sigma~^-" }}, + { 3224, { 3224, 1, 1.38280 , "Sigma*^+" }}, + { -3224, { -3224, -1, 1.38280 , "Sigma*~^-" }}, + { 3312, { 3312, -1, 1.32130 , "Xi^-" }}, + { -3312, { -3312, 1, 1.32130 , "Xi~^+" }}, + { 3314, { 3314, -1, 1.53500 , "Xi*^-" }}, + { -3314, { -3314, 1, 1.53500 , "Xi*~^+" }}, + { 3322, { 3322, 0, 1.31490 , "Xi^0" }}, + { -3322, { -3322, 0, 1.31490 , "Xi~^0" }}, + { 3324, { 3324, 0, 1.53180 , "Xi*^0" }}, + { -3324, { -3324, 0, 1.53180 , "Xi*~^0" }}, + { 3334, { 3334, -1, 1.67245 , "Omega^-" }}, + { -3334, { -3334, 1, 1.67245 , "Omega~^+" }}, + { 4112, { 4112, 0, 2.45210 , "Sigma_c^0" }}, + { -4112, { -4112, 0, 2.45210 , "Sigma_c~^0" }}, + { 4114, { 4114, 0, 2.50000 , "Sigma*_c^0" }}, + { -4114, { -4114, 0, 2.50000 , "Sigma*_c~^0" }}, + { 4122, { 4122, 1, 2.28490 , "Lambda_c^+" }}, + { -4122, { -4122, -1, 2.28490 , "Lambda_c~^-" }}, + { 4132, { 4132, 0, 2.47030 , "Xi_c^0" }}, + { -4132, { -4132, 0, 2.47030 , "Xi_c~^0" }}, + { 4212, { 4212, 1, 2.45350 , "Sigma_c^+" }}, + { -4212, { -4212, -1, 2.45350 , "Sigma_c~^-" }}, + { 4214, { 4214, 1, 2.50000 , "Sigma*_c^+" }}, + { -4214, { -4214, -1, 2.50000 , "Sigma*_c~^-" }}, + { 4222, { 4222, 2, 2.45290 , "Sigma_c^++" }}, + { -4222, { -4222, -2, 2.45290 , "Sigma_c~^--" }}, + { 4224, { 4224, 2, 2.50000 , "Sigma*_c^++" }}, + { -4224, { -4224, -2, 2.50000 , "Sigma*_c~^--" }}, + { 4232, { 4232, 1, 2.46560 , "Xi_c^+" }}, + { -4232, { -4232, -1, 2.46560 , "Xi_c~^-" }}, + { 4312, { 4312, 0, 2.55000 , "Xi'_c^0" }}, + { -4312, { -4312, 0, 2.55000 , "Xi'_c~^0" }}, + { 4314, { 4314, 0, 2.63000 , "Xi*_c^0" }}, + { -4314, { -4314, 0, 2.63000 , "Xi*_c~^0" }}, + { 4322, { 4322, 1, 2.55000 , "Xi'_c^+" }}, + { -4322, { -4322, -1, 2.55000 , "Xi'_c~^-" }}, + { 4324, { 4324, 1, 2.63000 , "Xi*_c^+" }}, + { -4324, { -4324, -1, 2.63000 , "Xi*_c~^-" }}, + { 4332, { 4332, 0, 2.70400 , "Omega_c^0" }}, + { -4332, { -4332, 0, 2.70400 , "Omega_c~^0" }}, + { 4334, { 4334, 0, 2.80000 , "Omega*_c^0" }}, + { -4334, { -4334, 0, 2.80000 , "Omega*_c~^0" }}, + { 4412, { 4412, 1, 3.59798 , "Xi_cc^+" }}, + { -4412, { -4412, -1, 3.59798 , "Xi_cc~^-" }}, + { 4414, { 4414, 1, 3.65648 , "Xi*_cc^+" }}, + { -4414, { -4414, -1, 3.65648 , "Xi*_cc~^-" }}, + { 4422, { 4422, 2, 3.59798 , "Xi_cc^++" }}, + { -4422, { -4422, -2, 3.59798 , "Xi_cc~^--" }}, + { 4424, { 4424, 2, 3.65648 , "Xi*_cc^++" }}, + { -4424, { -4424, -2, 3.65648 , "Xi*_cc~^--" }}, + { 4432, { 4432, 1, 3.78663 , "Omega_cc^+" }}, + { -4432, { -4432, -1, 3.78663 , "Omega_cc~^-" }}, + { 4434, { 4434, 1, 3.82466 , "Omega*_cc^+" }}, + { -4434, { -4434, -1, 3.82466 , "Omega*_cc~^-" }}, + { 4444, { 4444, 2, 4.91594 , "Omega*_ccc^++" }}, + { -4444, { -4444, -2, 4.91594 , "Omega*_ccc~^--" }}, + { 5112, { 5112, -1, 5.80000 , "Sigma_b^-" }}, + { -5112, { -5112, 1, 5.80000 , "Sigma_b~^+" }}, + { 5114, { 5114, -1, 5.81000 , "Sigma*_b^-" }}, + { -5114, { -5114, 1, 5.81000 , "Sigma*_b~^+" }}, + { 5122, { 5122, 0, 5.64100 , "Lambda_b^0" }}, + { -5122, { -5122, 0, 5.64100 , "Lambda_b~^0" }}, + { 5132, { 5132, -1, 5.84000 , "Xi_b^-" }}, + { -5132, { -5132, 1, 5.84000 , "Xi_b~^+" }}, + { 5142, { 5142, 0, 7.00575 , "Xi_bc^0" }}, + { -5142, { -5142, 0, 7.00575 , "Xi_bc~^0" }}, + { 5212, { 5212, 0, 5.80000 , "Sigma_b^0" }}, + { -5212, { -5212, 0, 5.80000 , "Sigma_b~^0" }}, + { 5214, { 5214, 0, 5.81000 , "Sigma*_b^0" }}, + { -5214, { -5214, 0, 5.81000 , "Sigma*_b~^0" }}, + { 5222, { 5222, 1, 5.80000 , "Sigma_b^+" }}, + { -5222, { -5222, -1, 5.80000 , "Sigma_b~^-" }}, + { 5224, { 5224, 1, 5.81000 , "Sigma*_b^+" }}, + { -5224, { -5224, -1, 5.81000 , "Sigma*_b~^-" }}, + { 5232, { 5232, 0, 5.84000 , "Xi_b^0" }}, + { -5232, { -5232, 0, 5.84000 , "Xi_b~^0" }}, + { 5242, { 5242, 1, 7.00575 , "Xi_bc^+" }}, + { -5242, { -5242, -1, 7.00575 , "Xi_bc~^-" }}, + { 5312, { 5312, -1, 5.96000 , "Xi'_b^-" }}, + { -5312, { -5312, 1, 5.96000 , "Xi'_b~^+" }}, + { 5314, { 5314, -1, 5.97000 , "Xi*_b^-" }}, + { -5314, { -5314, 1, 5.97000 , "Xi*_b~^+" }}, + { 5322, { 5322, 0, 5.96000 , "Xi'_b^0" }}, + { -5322, { -5322, 0, 5.96000 , "Xi'_b~^0" }}, + { 5324, { 5324, 0, 5.97000 , "Xi*_b^0" }}, + { -5324, { -5324, 0, 5.97000 , "Xi*_b~^0" }}, + { 5332, { 5332, -1, 6.12000 , "Omega_b^-" }}, + { -5332, { -5332, 1, 6.12000 , "Omega_b~^+" }}, + { 5334, { 5334, -1, 6.13000 , "Omega*_b^-" }}, + { -5334, { -5334, 1, 6.13000 , "Omega*_b~^+" }}, + { 5342, { 5342, 0, 7.19099 , "Omega_bc^0" }}, + { -5342, { -5342, 0, 7.19099 , "Omega_bc~^0" }}, + { 5412, { 5412, 0, 7.03724 , "Xi'_bc^0" }}, + { -5412, { -5412, 0, 7.03724 , "Xi'_bc~^0" }}, + { 5414, { 5414, 0, 7.04850 , "Xi*_bc^0" }}, + { -5414, { -5414, 0, 7.04850 , "Xi*_bc~^0" }}, + { 5422, { 5422, 1, 7.03724 , "Xi'_bc^+" }}, + { -5422, { -5422, -1, 7.03724 , "Xi'_bc~^-" }}, + { 5424, { 5424, 1, 7.04850 , "Xi*_bc^+" }}, + { -5424, { -5424, -1, 7.04850 , "Xi*_bc~^-" }}, + { 5432, { 5432, 0, 7.21101 , "Omega'_bc^0" }}, + { -5432, { -5432, 0, 7.21101 , "Omega'_bc~^0" }}, + { 5434, { 5434, 0, 7.21900 , "Omega*_bc^0" }}, + { -5434, { -5434, 0, 7.21900 , "Omega*_bc~^0" }}, + { 5442, { 5442, 1, 8.30945 , "Omega_bcc^+" }}, + { -5442, { -5442, -1, 8.30945 , "Omega_bcc~^-" }}, + { 5444, { 5444, 1, 8.31325 , "Omega*_bcc^+" }}, + { -5444, { -5444, -1, 8.31325 , "Omega*_bcc~^-" }}, + { 5512, { 5512, -1, 10.42272 , "Xi_bb^-" }}, + { -5512, { -5512, 1, 10.42272 , "Xi_bb~^+" }}, + { 5514, { 5514, -1, 10.44144 , "Xi*_bb^-" }}, + { -5514, { -5514, 1, 10.44144 , "Xi*_bb~^+" }}, + { 5522, { 5522, 0, 10.42272 , "Xi_bb^0" }}, + { -5522, { -5522, 0, 10.42272 , "Xi_bb~^0" }}, + { 5524, { 5524, 0, 10.44144 , "Xi*_bb^0" }}, + { -5524, { -5524, 0, 10.44144 , "Xi*_bb~^0" }}, + { 5532, { 5532, -1, 10.60209 , "Omega_bb^-" }}, + { -5532, { -5532, 1, 10.60209 , "Omega_bb~^+" }}, + { 5534, { 5534, -1, 10.61426 , "Omega*_bb^-" }}, + { -5534, { -5534, 1, 10.61426 , "Omega*_bb~^+" }}, + { 5542, { 5542, 0, 11.70767 , "Omega_bbc^0" }}, + { -5542, { -5542, 0, 11.70767 , "Omega_bbc~^0" }}, + { 5544, { 5544, 0, 11.71147 , "Omega*_bbc^0" }}, + { -5544, { -5544, 0, 11.71147 , "Omega*_bbc~^0" }}, + { 5554, { 5554, -1, 15.11061 , "Omega*_bbb^-" }}, + { -5554, { -5554, 1, 15.11061 , "Omega*_bbb~^+" }}, + { 1000010020, { 1000010020, 1, 1.87561 , "Deuterium" }}, + { 1000010030, { 1000010030, 1, 2.80925 , "Tritium" }}, + { 1000020030, { 1000020030, 2, 2.80923 , "He-3" }}, + { 1000020040, { 1000020040, 2, 3.72742 , "Alpha" }}, +}); + +} // namespace algorithms diff --git a/src/algorithms/interfaces/ParticleSvc.h b/src/algorithms/interfaces/ParticleSvc.h new file mode 100644 index 0000000000..5f3db7311c --- /dev/null +++ b/src/algorithms/interfaces/ParticleSvc.h @@ -0,0 +1,51 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Wouter Deconinck + +#pragma once + +#include +#include +#include +#include + +namespace algorithms { + +class ParticleSvc : public Service { +public: + struct ParticleData { + int pdgCode; + int charge; + double mass; + std::string name; + }; + using Particle = ParticleData; + using ParticleMap = std::map; + +private: + static const std::shared_ptr kParticleMap; + +public: + virtual void init(std::shared_ptr map = kParticleMap) { + if (map != nullptr) { + m_particleMap = map; + } + } + + virtual std::shared_ptr particleMap() const { + return m_particleMap; + }; + + virtual Particle& particle(int pdg) const { + if (m_particleMap->count(pdg) == 0) { + return m_particleMap->at(0); + } + return m_particleMap->at(pdg); + }; + +protected: + std::shared_ptr m_particleMap{nullptr}; + + ALGORITHMS_DEFINE_SERVICE(ParticleSvc) +}; + +} // namespace algorithms diff --git a/src/algorithms/onnx/InclusiveKinematicsML.cc b/src/algorithms/onnx/InclusiveKinematicsML.cc index 935efadb77..2483d6c629 100644 --- a/src/algorithms/onnx/InclusiveKinematicsML.cc +++ b/src/algorithms/onnx/InclusiveKinematicsML.cc @@ -32,26 +32,36 @@ namespace eicrecon { void InclusiveKinematicsML::init() { // onnxruntime setup - Ort::Env env(ORT_LOGGING_LEVEL_WARNING, "inclusive-kinematics-ml"); + m_env = Ort::Env(ORT_LOGGING_LEVEL_WARNING, "inclusive-kinematics-ml"); Ort::SessionOptions session_options; try { - m_session = Ort::Session(env, m_cfg.modelPath.c_str(), session_options); + m_session = Ort::Session(m_env, m_cfg.modelPath.c_str(), session_options); // print name/shape of inputs Ort::AllocatorWithDefaultOptions allocator; debug("Input Node Name/Shape:"); for (std::size_t i = 0; i < m_session.GetInputCount(); i++) { m_input_names.emplace_back(m_session.GetInputNameAllocated(i, allocator).get()); - m_input_shapes.emplace_back(m_session.GetInputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); - debug("\t{} : {}", m_input_names.at(i), print_shape(m_input_shapes.at(i))); + if (m_session.GetInputTypeInfo(i).GetONNXType() == ONNX_TYPE_TENSOR) { + m_input_shapes.emplace_back(m_session.GetInputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); + debug("\t{} : {}", m_input_names.at(i), print_shape(m_input_shapes.at(i))); + } else { + m_input_shapes.emplace_back(); + debug("\t{} : not a tensor", m_input_names.at(i)); + } } // print name/shape of outputs debug("Output Node Name/Shape:"); for (std::size_t i = 0; i < m_session.GetOutputCount(); i++) { m_output_names.emplace_back(m_session.GetOutputNameAllocated(i, allocator).get()); - m_output_shapes.emplace_back(m_session.GetOutputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); - debug("\t{} : {}", m_output_names.at(i), print_shape(m_output_shapes.at(i))); + if (m_session.GetOutputTypeInfo(i).GetONNXType() == ONNX_TYPE_TENSOR) { + m_output_shapes.emplace_back(m_session.GetOutputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); + debug("\t{} : {}", m_output_names.at(i), print_shape(m_output_shapes.at(i))); + } else { + m_output_shapes.emplace_back(); + debug("\t{} : not a tensor", m_output_names.at(i)); + } } // convert names to char* diff --git a/src/algorithms/onnx/InclusiveKinematicsML.h b/src/algorithms/onnx/InclusiveKinematicsML.h index 485fdb5d8b..8db3da2326 100644 --- a/src/algorithms/onnx/InclusiveKinematicsML.h +++ b/src/algorithms/onnx/InclusiveKinematicsML.h @@ -36,6 +36,7 @@ class InclusiveKinematicsML : public InclusiveKinematicsMLAlgorithm, void process(const Input&, const Output&) const final; private: + mutable Ort::Env m_env{nullptr}; mutable Ort::Session m_session{nullptr}; std::vector m_input_names; diff --git a/src/algorithms/pid/IrtCherenkovParticleID.cc b/src/algorithms/pid/IrtCherenkovParticleID.cc index 0337c569be..960a43ebea 100644 --- a/src/algorithms/pid/IrtCherenkovParticleID.cc +++ b/src/algorithms/pid/IrtCherenkovParticleID.cc @@ -111,11 +111,9 @@ void IrtCherenkovParticleID::init( } // get PDG info for the particles we want to identify in PID - // FIXME: cannot use `TDatabasePDG` since it is not thread safe; until we - // have a proper PDG database service, we hard-code the masses in Tools.h m_log->debug("List of particles for PID:"); for(auto pdg : m_cfg.pdgList) { - auto mass = Tools::GetPDGMass(pdg); + auto mass = m_particleSvc.particle(pdg).mass; m_pdg_mass.insert({ pdg, mass }); m_log->debug(" {:>8} M={} GeV", pdg, mass); } diff --git a/src/algorithms/pid/IrtCherenkovParticleID.h b/src/algorithms/pid/IrtCherenkovParticleID.h index 79828c2f60..0ddc648a9f 100644 --- a/src/algorithms/pid/IrtCherenkovParticleID.h +++ b/src/algorithms/pid/IrtCherenkovParticleID.h @@ -21,6 +21,7 @@ // EICrecon #include "IrtCherenkovParticleIDConfig.h" +#include "algorithms/interfaces/ParticleSvc.h" #include "algorithms/interfaces/WithPodConfig.h" namespace eicrecon { @@ -69,6 +70,8 @@ namespace eicrecon { CherenkovDetectorCollection* m_irt_det_coll; CherenkovDetector* m_irt_det; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + uint64_t m_cell_mask; std::string m_det_name; std::unordered_map m_pdg_mass; diff --git a/src/algorithms/pid/Tools.h b/src/algorithms/pid/Tools.h index b74d936fa9..b99450b122 100644 --- a/src/algorithms/pid/Tools.h +++ b/src/algorithms/pid/Tools.h @@ -62,34 +62,6 @@ namespace eicrecon { } - // ------------------------------------------------------------------------------------- - // PDG mass lookup - - // local PDG mass database - // FIXME: cannot use `TDatabasePDG` since it is not thread safe; until we - // have a proper PDG database service, we hard-code the masses we need; - // use Tools::GetPDGMass for access - static std::unordered_map GetPDGMasses() { - return std::unordered_map{ - { 11, 0.000510999 }, - { 211, 0.13957 }, - { 321, 0.493677 }, - { 2212, 0.938272 } - }; - } - - static double GetPDGMass(int pdg) { - double mass; - try { mass = GetPDGMasses().at(std::abs(pdg)); } - catch(const std::out_of_range& e) { - throw std::runtime_error(fmt::format("RUNTIME ERROR: unknown PDG={} in algorithms/pid/Tools::GetPDGMass",pdg)); - } - return mass; - } - - static int GetNumPDGs() { return GetPDGMasses().size(); }; - - // ------------------------------------------------------------------------------------- // Table rebinning and lookup diff --git a/src/algorithms/pid_lut/PIDLookup.cc b/src/algorithms/pid_lut/PIDLookup.cc index b0b35d50e8..fd132178d3 100644 --- a/src/algorithms/pid_lut/PIDLookup.cc +++ b/src/algorithms/pid_lut/PIDLookup.cc @@ -130,6 +130,7 @@ void PIDLookup::process(const Input& input, const Output& output) const { if (identified_pdg != 0) { recopart.setPDG(std::copysign(identified_pdg, (identified_pdg == 11) ? -charge : charge)); + recopart.setMass(m_particleSvc.particle(identified_pdg).mass); } if (identified_pdg != 0) { diff --git a/src/algorithms/pid_lut/PIDLookup.h b/src/algorithms/pid_lut/PIDLookup.h index 3f05c41573..82969d6207 100644 --- a/src/algorithms/pid_lut/PIDLookup.h +++ b/src/algorithms/pid_lut/PIDLookup.h @@ -10,6 +10,7 @@ #include #include "PIDLookupConfig.h" +#include "algorithms/interfaces/ParticleSvc.h" #include "algorithms/interfaces/WithPodConfig.h" #include "services/pid_lut/PIDLookupTable.h" @@ -38,6 +39,7 @@ class PIDLookup : public PIDLookupAlgorithm, public WithPodConfig m_dist{0, 1}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); const PIDLookupTable* m_lut; }; diff --git a/src/algorithms/reco/FarForwardNeutronReconstruction.cc b/src/algorithms/reco/FarForwardNeutronReconstruction.cc index b8c42ffdc5..14b63a504a 100644 --- a/src/algorithms/reco/FarForwardNeutronReconstruction.cc +++ b/src/algorithms/reco/FarForwardNeutronReconstruction.cc @@ -4,70 +4,82 @@ #include #include #include +#include #include #include -#include #include +#include #include "FarForwardNeutronReconstruction.h" /** - Creates a "neutron candidate" Reconstructed Particle consisting of all clusters in a - given ClusterCollection. Its energy is the sum of the energies of the constituent clusters + Creates a "neutron candidate" Reconstructed Particle consisting of all clusters from a + specified hadronic calorimeter and (optionally) from a specified Electromagnetic calorimeter. + Its energy is the sum of the energies of the constituent clusters times a correction factor, and its direction is the direction from the origin to the position of the most energetic cluster. The correction factor is given by 1/(1+c[0]+c[1]/sqrt(E)+c[2]/E), where c is the coefficients and E is the uncorrected energy in GeV. This form was chosen empirically based on the discrepancies in single-neutron MC simulations between the uncorrected - reconstructed energies and the truth energies of the neutrons. + reconstructed energies and the truth energies of the neutrons. Separate correction factors + are included for the Hcal and Ecal. */ namespace eicrecon { void FarForwardNeutronReconstruction::init() { - if (m_cfg.scale_corr_coeff.size() < 3) { - error("Invalid configuration. m_cfg.scale_corr_coeff should have at least 3 parameters"); - throw std::runtime_error("Invalid configuration. m_cfg.scale_corr_coeff should have at least 3 parameters"); + if (m_cfg.scale_corr_coeff_hcal.size() < 3) { + error("Invalid configuration. m_cfg.scale_corr_coeff_hcal should have at least 3 parameters"); + throw std::runtime_error("Invalid configuration. m_cfg.scale_corr_coeff_hcal should have at least 3 parameters"); + } + if (m_cfg.scale_corr_coeff_ecal.size() < 3) { + error("Invalid configuration. m_cfg.scale_corr_coeff_ecal should have at least 3 parameters"); + throw std::runtime_error("Invalid configuration. m_cfg.scale_corr_coeff_ecal should have at least 3 parameters"); } } - double FarForwardNeutronReconstruction::calc_corr(double Etot) const{ - auto coeffs=m_cfg.scale_corr_coeff; + /** calculates the correction for a given uncorrected total energy and a set of coefficients*/ + double FarForwardNeutronReconstruction::calc_corr(double Etot, const std::vector& coeffs) const{ return coeffs[0]+coeffs[1]/sqrt(Etot)+coeffs[2]/Etot; } void FarForwardNeutronReconstruction::process(const FarForwardNeutronReconstruction::Input& input, const FarForwardNeutronReconstruction::Output& output) const { - const auto [clusters] = input; + const auto [clustersHcal,clustersEcal] = input; auto [out_neutrons] = output; - double Etot=0; + double Etot_hcal=0, Etot_ecal=0; double Emax=0; - double x=0; - double y=0; - double z=0; - for (const auto& cluster : *clusters) { + edm4hep::Vector3f position; + for (const auto& cluster : *clustersHcal) { double E = cluster.getEnergy(); - Etot+=E; - if(E>Emax){ - Emax=E; - x=cluster.getPosition().x; - y=cluster.getPosition().y; - z=cluster.getPosition().z; + Etot_hcal += E; + if (E > Emax){ + Emax = E; + position = cluster.getPosition(); } } - if (Etot>0){ + for (const auto& cluster : *clustersEcal) { + double E = cluster.getEnergy(); + Etot_ecal+=E; + } + double Etot=Etot_hcal+Etot_ecal; + if (Etot > 0 && Emax > 0){ auto rec_part = out_neutrons->create(); - double corr=calc_corr(Etot); - Etot=Etot/(1+corr); + double corr=calc_corr(Etot,m_cfg.scale_corr_coeff_hcal); + Etot_hcal=Etot_hcal/(1+corr); + corr=calc_corr(Etot,m_cfg.scale_corr_coeff_ecal); + Etot_ecal=Etot_ecal/(1+corr); + Etot=Etot_hcal+Etot_ecal; rec_part.setEnergy(Etot); rec_part.setPDG(2112); - double p=sqrt(Etot*Etot-m_neutron*m_neutron); - double r=sqrt(x*x+y*y+z*z); - double px=p*x/r; - double py=p*y/r; - double pz=p*z/r; - rec_part.setMomentum({(float)px, (float)py, (float)pz}); + double p = sqrt(Etot*Etot-m_neutron*m_neutron); + double r = edm4hep::utils::magnitude(position); + edm4hep::Vector3f momentum = position * (p / r); + rec_part.setMomentum(momentum); rec_part.setCharge(0); rec_part.setMass(m_neutron); - for (const auto& cluster : *clusters){ + for (const auto& cluster : *clustersHcal){ + rec_part.addToClusters(cluster); + } + for (const auto& cluster : *clustersEcal){ rec_part.addToClusters(cluster); } } diff --git a/src/algorithms/reco/FarForwardNeutronReconstruction.h b/src/algorithms/reco/FarForwardNeutronReconstruction.h index 02bff54dfa..8767cbac3f 100644 --- a/src/algorithms/reco/FarForwardNeutronReconstruction.h +++ b/src/algorithms/reco/FarForwardNeutronReconstruction.h @@ -3,13 +3,16 @@ #pragma once +#include #include #include #include #include +#include #include // for basic_string #include // for string_view -#include +#include + #include "algorithms/interfaces/WithPodConfig.h" #include "algorithms/reco/FarForwardNeutronReconstructionConfig.h" @@ -17,7 +20,8 @@ namespace eicrecon { using FarForwardNeutronReconstructionAlgorithm = algorithms::Algorithm< algorithms::Input< - const edm4eic::ClusterCollection + const edm4eic::ClusterCollection, + std::optional >, algorithms::Output< edm4eic::ReconstructedParticleCollection @@ -29,13 +33,13 @@ using FarForwardNeutronReconstructionAlgorithm = algorithms::Algorithm< public: FarForwardNeutronReconstruction(std::string_view name) : FarForwardNeutronReconstructionAlgorithm{name, - {"inputClusters"}, + {"inputClustersHcal", "inputClustersEcal"}, {"outputNeutrons"}, - "Merges all HCAL clusters in a collection into a neutron candidate"} {} + "Merges all HCAL (and optionally also ECAL) clusters in a collection into a neutron candidate"} {} void init() final; void process(const Input&, const Output&) const final; - double calc_corr(double Etot) const; + double calc_corr(double Etot, const std::vector&) const; private: std::shared_ptr m_log; double m_neutron{0.93956542052}; diff --git a/src/algorithms/reco/FarForwardNeutronReconstructionConfig.h b/src/algorithms/reco/FarForwardNeutronReconstructionConfig.h index d73c36b9b9..f236f722cf 100644 --- a/src/algorithms/reco/FarForwardNeutronReconstructionConfig.h +++ b/src/algorithms/reco/FarForwardNeutronReconstructionConfig.h @@ -6,8 +6,10 @@ namespace eicrecon { struct FarForwardNeutronReconstructionConfig { - std::vector scale_corr_coeff={-0.0756, -1.91, 2.30}; - + /** Correction factors for the Hcal */ + std::vector scale_corr_coeff_hcal={-0.0756, -1.91, 2.30}; + /** Correction factors for the (optional) Ecal */ + std::vector scale_corr_coeff_ecal={-0.352, -1.34, 1.61}; }; } // eicrecon diff --git a/src/algorithms/reco/HadronicFinalState.cc b/src/algorithms/reco/HadronicFinalState.cc index 016770a85b..60ebb5f5a5 100644 --- a/src/algorithms/reco/HadronicFinalState.cc +++ b/src/algorithms/reco/HadronicFinalState.cc @@ -25,13 +25,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void HadronicFinalState::init() { - // m_pidSvc = service("ParticleSvc"); - // if (!m_pidSvc) { - // debug("Unable to locate Particle Service. " - // "Make sure you have ParticleSvc in the configuration."); - // } - } + void HadronicFinalState::init() { } void HadronicFinalState::process( const HadronicFinalState::Input& input, @@ -49,7 +43,7 @@ namespace eicrecon { const PxPyPzEVector ei( round_beam_four_momentum( ei_coll[0].getMomentum(), - m_electron, + m_particleSvc.particle(ei_coll[0].getPDG()).mass, {-5.0, -10.0, -18.0}, 0.0) ); @@ -63,7 +57,7 @@ namespace eicrecon { const PxPyPzEVector pi( round_beam_four_momentum( pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, + m_particleSvc.particle(pi_coll[0].getPDG()).mass, {41.0, 100.0, 275.0}, m_crossingAngle) ); diff --git a/src/algorithms/reco/HadronicFinalState.h b/src/algorithms/reco/HadronicFinalState.h index 6b25f48d7f..2ceaa8c82c 100644 --- a/src/algorithms/reco/HadronicFinalState.h +++ b/src/algorithms/reco/HadronicFinalState.h @@ -11,6 +11,8 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" + namespace eicrecon { using HadronicFinalStateAlgorithm = algorithms::Algorithm< @@ -31,7 +33,8 @@ class HadronicFinalState : public HadronicFinalStateAlgorithm { void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsDA.cc b/src/algorithms/reco/InclusiveKinematicsDA.cc index d925f07f01..b2260d268c 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.cc +++ b/src/algorithms/reco/InclusiveKinematicsDA.cc @@ -23,13 +23,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsDA::init() { - // m_pidSvc = service("ParticleSvc"); - // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " - // "Make sure you have ParticleSvc in the configuration."); - // } - } + void InclusiveKinematicsDA::init() { } void InclusiveKinematicsDA::process( const InclusiveKinematicsDA::Input& input, @@ -47,7 +41,7 @@ namespace eicrecon { const PxPyPzEVector ei( round_beam_four_momentum( ei_coll[0].getMomentum(), - m_electron, + m_particleSvc.particle(ei_coll[0].getPDG()).mass, {-5.0, -10.0, -18.0}, 0.0) ); @@ -61,7 +55,7 @@ namespace eicrecon { const PxPyPzEVector pi( round_beam_four_momentum( pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, + m_particleSvc.particle(pi_coll[0].getPDG()).mass, {41.0, 100.0, 275.0}, m_crossingAngle) ); @@ -87,6 +81,7 @@ namespace eicrecon { } // Calculate kinematic variables + static const auto m_proton = m_particleSvc.particle(2212).mass; const auto y_da = tan(gamma_h/2.) / ( tan(theta_e/2.) + tan(gamma_h/2.) ); const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(gamma_h/2.)) ); const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da); diff --git a/src/algorithms/reco/InclusiveKinematicsDA.h b/src/algorithms/reco/InclusiveKinematicsDA.h index ab6b71f7de..c93e421bb8 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.h +++ b/src/algorithms/reco/InclusiveKinematicsDA.h @@ -11,6 +11,8 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" + namespace eicrecon { using InclusiveKinematicsDAAlgorithm = algorithms::Algorithm< @@ -32,7 +34,8 @@ class InclusiveKinematicsDA : public InclusiveKinematicsDAAlgorithm { void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.cc b/src/algorithms/reco/InclusiveKinematicsESigma.cc similarity index 84% rename from src/algorithms/reco/InclusiveKinematicseSigma.cc rename to src/algorithms/reco/InclusiveKinematicsESigma.cc index 56fb00eec1..e1249172dd 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicsESigma.cc @@ -17,23 +17,17 @@ #include "Beam.h" #include "Boost.h" -#include "InclusiveKinematicseSigma.h" +#include "InclusiveKinematicsESigma.h" using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicseSigma::init() { - // m_pidSvc = service("ParticleSvc"); - // if (!m_pidSvc) { - // debug("Unable to locate Particle Service. " - // "Make sure you have ParticleSvc in the configuration."); - // } - } + void InclusiveKinematicsESigma::init() { } - void InclusiveKinematicseSigma::process( - const InclusiveKinematicseSigma::Input& input, - const InclusiveKinematicseSigma::Output& output) const { + void InclusiveKinematicsESigma::process( + const InclusiveKinematicsESigma::Input& input, + const InclusiveKinematicsESigma::Output& output) const { const auto [mcparts, escat, hfs] = input; auto [kinematics] = output; @@ -47,7 +41,7 @@ namespace eicrecon { const PxPyPzEVector ei( round_beam_four_momentum( ei_coll[0].getMomentum(), - m_electron, + m_particleSvc.particle(ei_coll[0].getPDG()).mass, {-5.0, -10.0, -18.0}, 0.0) ); @@ -61,7 +55,7 @@ namespace eicrecon { const PxPyPzEVector pi( round_beam_four_momentum( pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, + m_particleSvc.particle(pi_coll[0].getPDG()).mass, {41.0, 100.0, 275.0}, m_crossingAngle) ); @@ -96,6 +90,7 @@ namespace eicrecon { const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig); const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig); + static const auto m_proton = m_particleSvc.particle(2212).mass; const auto Q2_esig = Q2_e; const auto x_esig = x_sig; const auto y_esig = Q2_esig / (4.*ei.energy()*pi.energy()*x_esig); //equivalent to (2*ei.energy() / sigma_tot)*y_sig diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.h b/src/algorithms/reco/InclusiveKinematicsESigma.h similarity index 71% rename from src/algorithms/reco/InclusiveKinematicseSigma.h rename to src/algorithms/reco/InclusiveKinematicsESigma.h index 030a37453d..63a8c0ae7d 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.h +++ b/src/algorithms/reco/InclusiveKinematicsESigma.h @@ -11,18 +11,20 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" + namespace eicrecon { -using InclusiveKinematicseSigmaAlgorithm = algorithms::Algorithm< +using InclusiveKinematicsESigmaAlgorithm = algorithms::Algorithm< algorithms::Input, algorithms::Output>; -class InclusiveKinematicseSigma : public InclusiveKinematicseSigmaAlgorithm { +class InclusiveKinematicsESigma : public InclusiveKinematicsESigmaAlgorithm { public: - InclusiveKinematicseSigma(std::string_view name) - : InclusiveKinematicseSigmaAlgorithm{ + InclusiveKinematicsESigma(std::string_view name) + : InclusiveKinematicsESigmaAlgorithm{ name, {"MCParticles", "scatteredElectron", "hadronicFinalState"}, {"inclusiveKinematics"}, @@ -32,7 +34,8 @@ class InclusiveKinematicseSigma : public InclusiveKinematicseSigmaAlgorithm { void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.cc b/src/algorithms/reco/InclusiveKinematicsElectron.cc index 8ac014f2ff..79fa7f07da 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.cc +++ b/src/algorithms/reco/InclusiveKinematicsElectron.cc @@ -23,13 +23,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsElectron::init() { - // m_pidSvc = service("ParticleSvc"); - // if (!m_pidSvc) { - // debug("Unable to locate Particle Service. " - // "Make sure you have ParticleSvc in the configuration."); - // } - } + void InclusiveKinematicsElectron::init() { } void InclusiveKinematicsElectron::process( const InclusiveKinematicsElectron::Input& input, @@ -95,7 +89,7 @@ namespace eicrecon { const PxPyPzEVector ei( round_beam_four_momentum( ei_coll[0].getMomentum(), - m_electron, + m_particleSvc.particle(ei_coll[0].getPDG()).mass, {-5.0, -10.0, -18.0}, 0.0) ); @@ -109,7 +103,7 @@ namespace eicrecon { const PxPyPzEVector pi( round_beam_four_momentum( pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, + m_particleSvc.particle(pi_coll[0].getPDG()).mass, {41.0, 100.0, 275.0}, m_crossingAngle) ); @@ -128,6 +122,7 @@ namespace eicrecon { } // DIS kinematics calculations + static const auto m_proton = m_particleSvc.particle(2212).mass; const auto ef = electrons.front(); const auto q = ei - ef; const auto q_dot_pi = q.Dot(pi); diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.h b/src/algorithms/reco/InclusiveKinematicsElectron.h index 8b4864d14a..5a9b55875c 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.h +++ b/src/algorithms/reco/InclusiveKinematicsElectron.h @@ -11,6 +11,8 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" + namespace eicrecon { using InclusiveKinematicsElectronAlgorithm = algorithms::Algorithm< @@ -32,7 +34,8 @@ class InclusiveKinematicsElectron : public InclusiveKinematicsElectronAlgorithm void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsJB.cc b/src/algorithms/reco/InclusiveKinematicsJB.cc index 48874d79bd..67b80af52c 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.cc +++ b/src/algorithms/reco/InclusiveKinematicsJB.cc @@ -20,13 +20,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsJB::init() { - // m_pidSvc = service("ParticleSvc"); - // if (!m_pidSvc) { - // debug("Unable to locate Particle Service. " - // "Make sure you have ParticleSvc in the configuration."); - // } - } + void InclusiveKinematicsJB::init() { } void InclusiveKinematicsJB::process( const InclusiveKinematicsJB::Input& input, @@ -44,7 +38,7 @@ namespace eicrecon { const PxPyPzEVector ei( round_beam_four_momentum( ei_coll[0].getMomentum(), - m_electron, + m_particleSvc.particle(ei_coll[0].getPDG()).mass, {-5.0, -10.0, -18.0}, 0.0) ); @@ -58,7 +52,7 @@ namespace eicrecon { const PxPyPzEVector pi( round_beam_four_momentum( pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, + m_particleSvc.particle(pi_coll[0].getPDG()).mass, {41.0, 100.0, 275.0}, m_crossingAngle) ); @@ -75,6 +69,7 @@ namespace eicrecon { } // Calculate kinematic variables + static const auto m_proton = m_particleSvc.particle(2212).mass; const auto y_jb = sigma_h / (2.*ei.energy()); const auto Q2_jb = ptsum*ptsum / (1. - y_jb); const auto x_jb = Q2_jb / (4.*ei.energy()*pi.energy()*y_jb); diff --git a/src/algorithms/reco/InclusiveKinematicsJB.h b/src/algorithms/reco/InclusiveKinematicsJB.h index bacc5a69f7..3f985bba3e 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.h +++ b/src/algorithms/reco/InclusiveKinematicsJB.h @@ -11,6 +11,8 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" + namespace eicrecon { using InclusiveKinematicsJBAlgorithm = algorithms::Algorithm< @@ -32,7 +34,8 @@ class InclusiveKinematicsJB : public InclusiveKinematicsJBAlgorithm { void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.cc b/src/algorithms/reco/InclusiveKinematicsSigma.cc index c63b8f3981..d6789bf2b3 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicsSigma.cc @@ -22,13 +22,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsSigma::init() { - // m_pidSvc = service("ParticleSvc"); - // if (!m_pidSvc) { - // debug("Unable to locate Particle Service. " - // "Make sure you have ParticleSvc in the configuration."); - // } - } + void InclusiveKinematicsSigma::init() { } void InclusiveKinematicsSigma::process( const InclusiveKinematicsSigma::Input& input, @@ -46,7 +40,7 @@ namespace eicrecon { const PxPyPzEVector ei( round_beam_four_momentum( ei_coll[0].getMomentum(), - m_electron, + m_particleSvc.particle(ei_coll[0].getPDG()).mass, {-5.0, -10.0, -18.0}, 0.0) ); @@ -60,7 +54,7 @@ namespace eicrecon { const PxPyPzEVector pi( round_beam_four_momentum( pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, + m_particleSvc.particle(pi_coll[0].getPDG()).mass, {41.0, 100.0, 275.0}, m_crossingAngle) ); @@ -88,6 +82,7 @@ namespace eicrecon { auto sigma_tot = sigma_e + sigma_h; // Calculate kinematic variables + static const auto m_proton = m_particleSvc.particle(2212).mass; const auto y_sig = sigma_h / sigma_tot; const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig); const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig); diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.h b/src/algorithms/reco/InclusiveKinematicsSigma.h index 6cf3d13cbe..82bd651607 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.h +++ b/src/algorithms/reco/InclusiveKinematicsSigma.h @@ -11,6 +11,8 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" + namespace eicrecon { using InclusiveKinematicsSigmaAlgorithm = algorithms::Algorithm< @@ -32,7 +34,8 @@ class InclusiveKinematicsSigma : public InclusiveKinematicsSigmaAlgorithm { void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.cc b/src/algorithms/reco/InclusiveKinematicsTruth.cc index e48bbc0d77..d05493e832 100644 --- a/src/algorithms/reco/InclusiveKinematicsTruth.cc +++ b/src/algorithms/reco/InclusiveKinematicsTruth.cc @@ -19,15 +19,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsTruth::init() { - // m_pidSvc = service("ParticleSvc"); - // if (!m_pidSvc) { - // error() << "Unable to locate Particle Service. " - // << "Make sure you have ParticleSvc in the configuration." - // << endmsg; - // return StatusCode::FAILURE; - // } - } + void InclusiveKinematicsTruth::init() { } void InclusiveKinematicsTruth::process( const InclusiveKinematicsTruth::Input& input, @@ -50,7 +42,7 @@ namespace eicrecon { } const auto ei_p = ei_coll[0].getMomentum(); const auto ei_p_mag = edm4hep::utils::magnitude(ei_p); - const auto ei_mass = m_electron; + static const auto ei_mass = m_particleSvc.particle(11).mass; const PxPyPzEVector ei(ei_p.x, ei_p.y, ei_p.z, std::hypot(ei_p_mag, ei_mass)); // Get incoming hadron beam @@ -61,7 +53,7 @@ namespace eicrecon { } const auto pi_p = pi_coll[0].getMomentum(); const auto pi_p_mag = edm4hep::utils::magnitude(pi_p); - const auto pi_mass = pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron; + const auto pi_mass = m_particleSvc.particle(pi_coll[0].getPDG()).mass; const PxPyPzEVector pi(pi_p.x, pi_p.y, pi_p.z, std::hypot(pi_p_mag, pi_mass)); // Get first scattered electron @@ -76,7 +68,7 @@ namespace eicrecon { } const auto ef_p = ef_coll[0].getMomentum(); const auto ef_p_mag = edm4hep::utils::magnitude(ef_p); - const auto ef_mass = m_electron; + static const auto ef_mass = m_particleSvc.particle(11).mass; const PxPyPzEVector ef(ef_p.x, ef_p.y, ef_p.z, std::hypot(ef_p_mag, ef_mass)); // DIS kinematics calculations diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.h b/src/algorithms/reco/InclusiveKinematicsTruth.h index 83a3f44576..39cfe05c7c 100644 --- a/src/algorithms/reco/InclusiveKinematicsTruth.h +++ b/src/algorithms/reco/InclusiveKinematicsTruth.h @@ -9,6 +9,8 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" + namespace eicrecon { using InclusiveKinematicsTruthAlgorithm = @@ -29,7 +31,8 @@ class InclusiveKinematicsTruth : public InclusiveKinematicsTruthAlgorithm { void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; } // namespace eicrecon diff --git a/src/algorithms/reco/ScatteredElectronsTruth.cc b/src/algorithms/reco/ScatteredElectronsTruth.cc index 3fead2e236..733fb311ae 100644 --- a/src/algorithms/reco/ScatteredElectronsTruth.cc +++ b/src/algorithms/reco/ScatteredElectronsTruth.cc @@ -91,12 +91,14 @@ namespace eicrecon { if (p.getObjectID() == ef_rc_id) { output_electrons->push_back( p ); + static const auto m_electron = m_particleSvc.particle(11).mass; vScatteredElectron.SetCoordinates( p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, m_electron ); electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); // break; NOTE: if we are not computing E-Pz // we can safely break here and save precious CPUs } else { // Compute the sum hadronic final state + static const auto m_pion = m_particleSvc.particle(211).mass; vHadron.SetCoordinates( p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, m_pion ); vHadronicFinalState += vHadron; } diff --git a/src/algorithms/reco/ScatteredElectronsTruth.h b/src/algorithms/reco/ScatteredElectronsTruth.h index 291d6cc3b3..33410cc289 100644 --- a/src/algorithms/reco/ScatteredElectronsTruth.h +++ b/src/algorithms/reco/ScatteredElectronsTruth.h @@ -10,6 +10,7 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" namespace eicrecon { @@ -38,7 +39,8 @@ namespace eicrecon { void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_pion{0.13957}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + }; } // namespace eicrecon diff --git a/src/algorithms/reco/TransformBreitFrame.cc b/src/algorithms/reco/TransformBreitFrame.cc index 5e168845d3..6cf37b46e1 100644 --- a/src/algorithms/reco/TransformBreitFrame.cc +++ b/src/algorithms/reco/TransformBreitFrame.cc @@ -40,7 +40,7 @@ namespace eicrecon { const PxPyPzEVector e_initial( round_beam_four_momentum( ei_coll[0].getMomentum(), - m_electron, + m_particleSvc.particle(ei_coll[0].getPDG()).mass, {-5.0, -10.0, -18.0}, 0.0) ); @@ -54,7 +54,7 @@ namespace eicrecon { const PxPyPzEVector p_initial( round_beam_four_momentum( pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, + m_particleSvc.particle(pi_coll[0].getPDG()).mass, {41.0, 100.0, 275.0}, m_crossingAngle) ); diff --git a/src/algorithms/reco/TransformBreitFrame.h b/src/algorithms/reco/TransformBreitFrame.h index 74a89dd69b..26ab2d06f7 100644 --- a/src/algorithms/reco/TransformBreitFrame.h +++ b/src/algorithms/reco/TransformBreitFrame.h @@ -10,6 +10,7 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" #include "algorithms/interfaces/WithPodConfig.h" namespace eicrecon { @@ -36,7 +37,8 @@ class TransformBreitFrame : public TransformBreitFrameAlgorithm, public WithPodC void process(const Input&, const Output&) const final; private: - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + double m_crossingAngle{-0.025}; }; // end TransformBreitFrame definition diff --git a/src/algorithms/reco/UndoAfterBurner.cc b/src/algorithms/reco/UndoAfterBurner.cc new file mode 100644 index 0000000000..69430e97e2 --- /dev/null +++ b/src/algorithms/reco/UndoAfterBurner.cc @@ -0,0 +1,142 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Alex Jentsch, Jihee Kim, Brian Page +// + +#include "UndoAfterBurner.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "algorithms/reco/Beam.h" +#include "algorithms/reco/UndoAfterBurnerConfig.h" + +void eicrecon::UndoAfterBurner::init() { + +} + +void eicrecon::UndoAfterBurner::process( + const UndoAfterBurner::Input& input, + const UndoAfterBurner::Output& output) const { + + const auto [mcparts] = input; + auto [outputParticles] = output; + + bool pidAssumePionMass = m_cfg.m_pid_assume_pion_mass; + double crossingAngle = m_cfg.m_crossing_angle; + double pidPurity = m_cfg.m_pid_purity; + bool correctBeamFX = m_cfg.m_correct_beam_FX; + bool pidUseMCTruth = m_cfg.m_pid_use_MC_truth; + + bool hasBeamHadron = true; + bool hasBeamLepton = true; + + //read MCParticles information and "postburn" to remove the afterburner effects. + //The output is then the original MC input produced by the generator. + + ROOT::Math::PxPyPzEVector e_beam(0.,0.,0.,0.); + ROOT::Math::PxPyPzEVector h_beam(0.,0.,0.,0.); + + auto incoming_lepton = find_first_beam_electron(mcparts); + if (incoming_lepton.size() == 0) { + debug("No beam electron found -- particleGun input"); + hasBeamLepton = false; + } + + auto incoming_hadron = find_first_beam_hadron(mcparts); + if (incoming_hadron.size() == 0) { + debug("No beam hadron found -- particleGun input"); + hasBeamHadron = false; + } + + if((hasBeamHadron && !hasBeamLepton) || (!hasBeamHadron && hasBeamLepton)){ + debug("Only one beam defined! Not a possible configuration!"); + return; + } + + // Handling for FF particle gun input!! + if (!hasBeamHadron || !hasBeamLepton) { + for (const auto& p: *mcparts) { + if((p.getPDG() == 2212 || p.getPDG() == 2112)) { //look for "gun" proton/neutron + hasBeamHadron = true; + h_beam.SetPxPyPzE(crossingAngle*p.getEnergy(), 0.0, p.getEnergy(), p.getEnergy()); + if(p.getEnergy() > 270.0 && p.getEnergy() < 280.0){ + hasBeamLepton = true; + e_beam.SetPxPyPzE(0.0, 0.0, -18.0, 18.0); + } + } + } + + } else { + + if (correctBeamFX) { + + h_beam.SetPxPyPzE(incoming_hadron[0].getMomentum().x, incoming_hadron[0].getMomentum().y, incoming_hadron[0].getMomentum().z, incoming_hadron[0].getEnergy()); + e_beam.SetPxPyPzE(incoming_lepton[0].getMomentum().x, incoming_lepton[0].getMomentum().y, incoming_lepton[0].getMomentum().z, incoming_lepton[0].getEnergy()); + + } else { + + h_beam.SetPxPyPzE(crossingAngle*incoming_hadron[0].getEnergy(), 0.0, incoming_hadron[0].getEnergy(), incoming_hadron[0].getEnergy()); + e_beam.SetPxPyPzE(0.0, 0.0, -incoming_lepton[0].getEnergy(), incoming_lepton[0].getEnergy()); + + } + } + + // Bail out if still no beam particles, since this leads to division by zero + if (!hasBeamHadron && !hasBeamLepton) { + return; + } + + // Calculate boost vectors and rotations here + ROOT::Math::PxPyPzEVector cm_frame_boost = e_beam + h_beam; + ROOT::Math::Cartesian3D beta(-cm_frame_boost.Px() / cm_frame_boost.E(), -cm_frame_boost.Py() / cm_frame_boost.E(), -cm_frame_boost.Pz() / cm_frame_boost.E()); + ROOT::Math::Boost boostVector(beta); + + // Boost to CM frame + e_beam = boostVector(e_beam); + h_beam = boostVector(h_beam); + + double rotationAngleY = -1.0*TMath::ATan2(h_beam.Px(), h_beam.Pz()); + double rotationAngleX = 1.0*TMath::ATan2(h_beam.Py(), h_beam.Pz()); + + ROOT::Math::RotationY rotationAboutY(rotationAngleY); + ROOT::Math::RotationX rotationAboutX(rotationAngleX); + + // Boost back to proper head-on frame + ROOT::Math::PxPyPzEVector head_on_frame_boost(0., 0., cm_frame_boost.Pz(), cm_frame_boost.E()); + ROOT::Math::Boost headOnBoostVector(head_on_frame_boost.Px()/head_on_frame_boost.E(), head_on_frame_boost.Py()/head_on_frame_boost.E(), head_on_frame_boost.Pz()/head_on_frame_boost.E()); + + // Now, loop through events and apply operations to the MCparticles + for (const auto& p: *mcparts) { + + ROOT::Math::PxPyPzEVector mc(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); + + mc = boostVector(mc); + mc = rotationAboutY(mc); + mc = rotationAboutX(mc); + mc = headOnBoostVector(mc); + + edm4hep::Vector3f mcMom(mc.Px(), mc.Py(), mc.Pz()); + edm4hep::MutableMCParticle MCTrack(p.clone()); + MCTrack.setMomentum(mcMom); + + if(pidUseMCTruth){ + MCTrack.setPDG(p.getPDG()); + MCTrack.setMass(p.getMass()); + } + if(!pidUseMCTruth && pidAssumePionMass){ + MCTrack.setPDG(211); + MCTrack.setMass(0.13957); + } + + outputParticles->push_back(MCTrack); + } + +} diff --git a/src/algorithms/reco/UndoAfterBurner.h b/src/algorithms/reco/UndoAfterBurner.h new file mode 100644 index 0000000000..57866520f0 --- /dev/null +++ b/src/algorithms/reco/UndoAfterBurner.h @@ -0,0 +1,42 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Alex Jentsch, Jihee Kim, Brian Page +// + +#include +#include +#include +#include + +#include "UndoAfterBurnerConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + + using UndoAfterBurnerAlgorithm = algorithms::Algorithm< + algorithms::Input< + edm4hep::MCParticleCollection + >, + algorithms::Output< + edm4hep::MCParticleCollection + > + >; + + class UndoAfterBurner + : public UndoAfterBurnerAlgorithm, + public WithPodConfig { + + public: + UndoAfterBurner(std::string_view name) + : UndoAfterBurnerAlgorithm{name, + {"inputMCParticles"}, + {"outputMCParticles"}, + "Apply boosts and rotations to remove crossing angle and beam effects."} {} + + void init(); + void process(const Input&, const Output&) const final; + + private: + + + }; +} diff --git a/src/algorithms/reco/UndoAfterBurnerConfig.h b/src/algorithms/reco/UndoAfterBurnerConfig.h new file mode 100644 index 0000000000..1c331a5a45 --- /dev/null +++ b/src/algorithms/reco/UndoAfterBurnerConfig.h @@ -0,0 +1,19 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Alex Jentsch, Jihee Kim, Brian Page +// + +#pragma once + +namespace eicrecon { + + struct UndoAfterBurnerConfig { + + bool m_pid_assume_pion_mass = false; + double m_crossing_angle = -0.025; + double m_pid_purity = 0.51; + bool m_correct_beam_FX = true; + bool m_pid_use_MC_truth = true; + + }; + +} diff --git a/src/algorithms/tracking/ActsToTracks.cc b/src/algorithms/tracking/ActsToTracks.cc index 179d1541cf..79480d1884 100644 --- a/src/algorithms/tracking/ActsToTracks.cc +++ b/src/algorithms/tracking/ActsToTracks.cc @@ -8,14 +8,23 @@ #include #include #include +#include +#include +#include +#include +#include #include #include #include +#include +#include #include #include #include #include #include +#include +#include #include #include @@ -42,8 +51,8 @@ void ActsToTracks::init() { } void ActsToTracks::process(const Input& input, const Output& output) const { - const auto [meas2Ds, acts_trajectories] = input; - auto [trajectories, track_parameters, tracks] = output; + const auto [meas2Ds, acts_trajectories, raw_hit_assocs] = input; + auto [trajectories, track_parameters, tracks, tracks_assoc] = output; // Loop over trajectories for (const auto traj : acts_trajectories) { @@ -147,6 +156,9 @@ void ActsToTracks::process(const Input& input, const Output& output) const { ); track.setTrajectory(trajectory); // Trajectory of this track + // Determine track association with MCParticle, weighted by number of used measurements + std::map mcparticle_weight_by_hit_count; + // save measurement2d to good measurements or outliers according to srclink index // fix me: ideally, this should be integrated into multitrajectoryhelper // fix me: should say "OutlierMeasurements" instead of "OutlierHits" etc @@ -173,6 +185,23 @@ void ActsToTracks::process(const Input& input, const Output& output) const { debug("Measurement on geo id={}, index={}, loc={},{}", geoID, srclink_index, meas2D.getLoc().a, meas2D.getLoc().b); + // Determine track associations if hit associations provided + // FIXME: not able to check whether optional inputs were provided + //if (raw_hit_assocs->has_value()) { + #if EDM4EIC_VERSION_MAJOR >= 7 + for (auto& hit : meas2D.getHits()) { + auto raw_hit = hit.getRawHit(); + for (const auto raw_hit_assoc : *raw_hit_assocs) { + if (raw_hit_assoc.getRawHit() == raw_hit) { + auto sim_hit = raw_hit_assoc.getSimHit(); + auto mc_particle = sim_hit.getMCParticle(); + mcparticle_weight_by_hit_count[mc_particle]++; + } + } + } + #endif + //} + } else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) { trajectory.addToOutliers_deprecated(meas2D); @@ -185,6 +214,22 @@ void ActsToTracks::process(const Input& input, const Output& output) const { }); + // Store track associations if hit associations provided + // FIXME: not able to check whether optional inputs were provided + //if (raw_hit_assocs->has_value()) { + double total_weight = std::accumulate( + mcparticle_weight_by_hit_count.begin(), mcparticle_weight_by_hit_count.end(), + 0, [](const double sum, const auto& i) { return sum + i.second; }); + for (const auto& [mcparticle, weight] : mcparticle_weight_by_hit_count) { + auto track_assoc = tracks_assoc->create(); + track_assoc.setRec(track); + track_assoc.setSim(mcparticle); + double normalized_weight = weight / total_weight; + track_assoc.setWeight(normalized_weight); + debug("track {}: mcparticle {} weight {}", track.id().index, mcparticle.id().index, normalized_weight); + } + //} + } } } diff --git a/src/algorithms/tracking/ActsToTracks.h b/src/algorithms/tracking/ActsToTracks.h index 95bf88ec9f..884ab25198 100644 --- a/src/algorithms/tracking/ActsToTracks.h +++ b/src/algorithms/tracking/ActsToTracks.h @@ -5,10 +5,13 @@ #include #include +#include +#include #include #include #include #include +#include #include #include #include @@ -17,13 +20,36 @@ namespace eicrecon { using ActsToTracksAlgorithm = algorithms::Algorithm< - algorithms::Input>, - algorithms::Output + algorithms::Input< + edm4eic::Measurement2DCollection, + std::vector, + std::optional + >, + algorithms::Output< + edm4eic::TrajectoryCollection, + edm4eic::TrackParametersCollection, + edm4eic::TrackCollection, + std::optional + > >; class ActsToTracks : public ActsToTracksAlgorithm { public: - ActsToTracks(std::string_view name) : ActsToTracksAlgorithm{name, {"inputActsTrajectories", "inputActsConstTrackContainer"}, {"outputTrajectoryCollection", "outputTrackParametersCollection", "outputTrackCollection"}, "Converts ACTS trajectories to EDM4eic"} {}; + ActsToTracks(std::string_view name) + : ActsToTracksAlgorithm{ + name, + { + "inputMeasurements", + "inputActsTrajectories", + "inputRawTrackerHitAssociations", + }, + { + "outputTrajectories", + "outputTrackParameters", + "outputTracks", + "outputTrackAssociations", + }, + "Converts ACTS trajectories to EDM4eic"} {}; void init() final; void process(const Input&, const Output&) const final; diff --git a/src/algorithms/tracking/CKFTracking.cc b/src/algorithms/tracking/CKFTracking.cc index de29d15f98..93d33feff7 100644 --- a/src/algorithms/tracking/CKFTracking.cc +++ b/src/algorithms/tracking/CKFTracking.cc @@ -7,21 +7,38 @@ #include #include #include +#if Acts_VERSION_MAJOR < 36 #include +#endif #include #include +#if Acts_VERSION_MAJOR >= 32 +#include "Acts/EventData/ProxyAccessor.hpp" +#endif #include #include #include #include #include #include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Propagator/AbortList.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/MaterialInteractor.hpp" +#include "Acts/Propagator/Navigator.hpp" +#endif #include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Propagator/StandardAborters.hpp" +#endif #include #include #include #include #include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Utilities/TrackHelpers.hpp" +#endif #include #include #include @@ -129,7 +146,13 @@ namespace eicrecon { cov(0, 1) = meas2D.getCovariance().xy; cov(1, 0) = meas2D.getCovariance().xy; - auto measurement = Acts::makeMeasurement(Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); +#if Acts_VERSION_MAJOR >= 36 + auto measurement = ActsExamples::makeFixedSizeMeasurement( + Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); +#else + auto measurement = Acts::makeMeasurement( + Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); +#endif measurements->emplace_back(std::move(measurement)); hit_index++; @@ -175,7 +198,9 @@ namespace eicrecon { ActsExamples::PassThroughCalibrator pcalibrator; ActsExamples::MeasurementCalibratorAdapter calibrator(pcalibrator, *measurements); Acts::GainMatrixUpdater kfUpdater; +#if Acts_VERSION_MAJOR < 34 Acts::GainMatrixSmoother kfSmoother; +#endif Acts::MeasurementSelector measSel{m_sourcelinkSelectorCfg}; Acts::CombinatorialKalmanFilterExtensions @@ -185,9 +210,11 @@ namespace eicrecon { extensions.updater.connect< &Acts::GainMatrixUpdater::operator()>( &kfUpdater); +#if Acts_VERSION_MAJOR < 34 extensions.smoother.connect< &Acts::GainMatrixSmoother::operator()>( &kfSmoother); +#endif extensions.measurementSelector.connect< &Acts::MeasurementSelector::select>( &measSel); @@ -199,9 +226,27 @@ namespace eicrecon { slAccessorDelegate.connect<&ActsExamples::IndexSourceLinkAccessor::range>(&slAccessor); // Set the CombinatorialKalmanFilter options +#if Acts_VERSION_MAJOR < 34 CKFTracking::TrackFinderOptions options( m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate, extensions, pOptions, &(*pSurface)); +#else + CKFTracking::TrackFinderOptions options( + m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate, + extensions, pOptions); +#endif + +#if Acts_VERSION_MAJOR >= 34 + Acts::Propagator, Acts::Navigator> extrapolator( + Acts::EigenStepper<>(m_BField), + Acts::Navigator({m_geoSvc->trackingGeometry()}, + logger().cloneWithSuffix("Navigator")), + logger().cloneWithSuffix("Propagator")); + + Acts::PropagatorOptions, + Acts::AbortList> + extrapolationOptions(m_geoctx, m_fieldctx); +#endif // Create track container auto trackContainer = std::make_shared(); @@ -210,7 +255,11 @@ namespace eicrecon { // Add seed number column acts_tracks.addColumn("seed"); +#if Acts_VERSION_MAJOR >= 32 + Acts::ProxyAccessor seedNumber("seed"); +#else Acts::TrackAccessor seedNumber("seed"); +#endif // Loop over seeds for (std::size_t iseed = 0; iseed < acts_init_trk_params.size(); ++iseed) { @@ -225,6 +274,27 @@ namespace eicrecon { // Set seed number for all found tracks auto& tracksForSeed = result.value(); for (auto& track : tracksForSeed) { + +#if Acts_VERSION_MAJOR >=34 + auto smoothingResult = Acts::smoothTrack(m_geoctx, track, logger()); + if (!smoothingResult.ok()) { + ACTS_ERROR("Smoothing for seed " + << iseed << " and track " << track.index() + << " failed with error " << smoothingResult.error()); + continue; + } + + auto extrapolationResult = Acts::extrapolateTrackToReferenceSurface( + track, *pSurface, extrapolator, extrapolationOptions, + Acts::TrackExtrapolationStrategy::firstOrLast, logger()); + if (!extrapolationResult.ok()) { + ACTS_ERROR("Extrapolation for seed " + << iseed << " and track " << track.index() + << " failed with error " << extrapolationResult.error()); + continue; + } +#endif + seedNumber(track) = iseed; } } @@ -251,8 +321,11 @@ namespace eicrecon { auto& constTracks = *(constTracks_v.front()); // Seed number column accessor +#if Acts_VERSION_MAJOR >= 32 + const Acts::ConstProxyAccessor constSeedNumber("seed"); +#else const Acts::ConstTrackAccessor constSeedNumber("seed"); - +#endif // Prepare the output data with MultiTrajectory, per seed std::vector acts_trajectories; diff --git a/src/algorithms/tracking/DD4hepBField.h b/src/algorithms/tracking/DD4hepBField.h index 257666c6af..12f088c235 100644 --- a/src/algorithms/tracking/DD4hepBField.h +++ b/src/algorithms/tracking/DD4hepBField.h @@ -43,7 +43,11 @@ namespace eicrecon::BField { Acts::MagneticFieldProvider::Cache makeCache(const Acts::MagneticFieldContext& mctx) const override { +#if Acts_VERSION_MAJOR >= 32 + return Acts::MagneticFieldProvider::Cache(std::in_place_type, mctx); +#else return Acts::MagneticFieldProvider::Cache::make(mctx); +#endif } /** construct constant magnetic field from field vector. diff --git a/src/algorithms/tracking/IterativeVertexFinder.cc b/src/algorithms/tracking/IterativeVertexFinder.cc index c33488ae31..c933a96feb 100644 --- a/src/algorithms/tracking/IterativeVertexFinder.cc +++ b/src/algorithms/tracking/IterativeVertexFinder.cc @@ -13,7 +13,12 @@ #include #include #include +#include +#if Acts_VERSION_MAJOR >= 32 +#include +#else #include +#endif #include #include #include @@ -57,47 +62,108 @@ std::unique_ptr eicrecon::IterativeVertexFinder::prod using Propagator = Acts::Propagator>; using PropagatorOptions = Acts::PropagatorOptions<>; +#if Acts_VERSION_MAJOR >= 33 + using Linearizer = Acts::HelicalTrackLinearizer; + using VertexFitter = Acts::FullBilloirVertexFitter; + using ImpactPointEstimator = Acts::ImpactPointEstimator; + using VertexSeeder = Acts::ZScanVertexFinder; + using VertexFinder = Acts::IterativeVertexFinder; + using VertexFinderOptions = Acts::VertexingOptions; +#else using Linearizer = Acts::HelicalTrackLinearizer; using VertexFitter = Acts::FullBilloirVertexFitter; using ImpactPointEstimator = Acts::ImpactPointEstimator; using VertexSeeder = Acts::ZScanVertexFinder; using VertexFinder = Acts::IterativeVertexFinder; using VertexFinderOptions = Acts::VertexingOptions; +#endif ACTS_LOCAL_LOGGER(eicrecon::getSpdlogLogger("IVF", m_log)); Acts::EigenStepper<> stepper(m_BField); // Set up propagator with void navigator +#if Acts_VERSION_MAJOR >= 32 + auto propagator = std::make_shared( + stepper, Acts::VoidNavigator{}, logger().cloneWithSuffix("Prop")); +#else auto propagator = std::make_shared( stepper, Acts::detail::VoidNavigator{}, logger().cloneWithSuffix("Prop")); +#endif Acts::PropagatorOptions opts(m_geoctx, m_fieldctx); - // Setup the vertex fitter - VertexFitter::Config vertexFitterCfg; - VertexFitter vertexFitter(vertexFitterCfg); // Setup the track linearizer +#if Acts_VERSION_MAJOR >= 33 + Linearizer::Config linearizerCfg; + linearizerCfg.bField = m_BField; + linearizerCfg.propagator = propagator; +#else Linearizer::Config linearizerCfg(m_BField, propagator); +#endif Linearizer linearizer(linearizerCfg, logger().cloneWithSuffix("HelLin")); + + // Setup the vertex fitter + VertexFitter::Config vertexFitterCfg; +#if Acts_VERSION_MAJOR >= 33 + vertexFitterCfg.extractParameters + .connect<&Acts::InputTrack::extractParameters>(); + vertexFitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>( + &linearizer); +#endif + VertexFitter vertexFitter(vertexFitterCfg); + // Setup the seed finder ImpactPointEstimator::Config ipEstCfg(m_BField, propagator); ImpactPointEstimator ipEst(ipEstCfg); VertexSeeder::Config seederCfg(ipEst); +#if Acts_VERSION_MAJOR >= 33 + seederCfg.extractParameters + .connect<&Acts::InputTrack::extractParameters>(); + auto seeder = std::make_shared(seederCfg); +#else VertexSeeder seeder(seederCfg); +#endif + // Set up the actual vertex finder - VertexFinder::Config finderCfg(std::move(vertexFitter), std::move(linearizer), - std::move(seeder), std::move(ipEst)); + VertexFinder::Config finderCfg( + std::move(vertexFitter), +#if Acts_VERSION_MAJOR < 33 + std::move(linearizer), +#endif + std::move(seeder), + std::move(ipEst)); finderCfg.maxVertices = m_cfg.maxVertices; finderCfg.reassignTracksAfterFirstFit = m_cfg.reassignTracksAfterFirstFit; - #if Acts_VERSION_MAJOR >= 31 - VertexFinder finder(std::move(finderCfg)); +#if Acts_VERSION_MAJOR >= 31 + #if Acts_VERSION_MAJOR >= 33 + finderCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>(); + finderCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer); + #if Acts_VERSION_MAJOR >= 36 + finderCfg.field = m_BField; #else - VertexFinder finder(finderCfg); + finderCfg.field = std::dynamic_pointer_cast( + std::const_pointer_cast(m_BField)); #endif + #endif + VertexFinder finder(std::move(finderCfg)); +#else + VertexFinder finder(finderCfg); +#endif +#if Acts_VERSION_MAJOR >= 33 + Acts::IVertexFinder::State state( + std::in_place_type, + *m_BField, + m_fieldctx); +#else VertexFinder::State state(*m_BField, m_fieldctx); +#endif VertexFinderOptions finderOpts(m_geoctx, m_fieldctx); +#if Acts_VERSION_MAJOR >= 33 + std::vector inputTracks; +#else std::vector inputTrackPointers; +#endif for (const auto& trajectory : trajectories) { auto tips = trajectory->tips(); @@ -106,12 +172,21 @@ std::unique_ptr eicrecon::IterativeVertexFinder::prod } /// CKF can provide multiple track trajectories for a single input seed for (auto& tip : tips) { +#if Acts_VERSION_MAJOR >= 33 + inputTracks.emplace_back(&(trajectory->trackParameters(tip))); +#else inputTrackPointers.push_back(&(trajectory->trackParameters(tip))); +#endif } } +#if Acts_VERSION_MAJOR >= 33 + std::vector vertices; + auto result = finder.find(inputTracks, finderOpts, state); +#else std::vector> vertices; auto result = finder.find(inputTrackPointers, finderOpts, state); +#endif if (result.ok()) { vertices = std::move(result.value()); } diff --git a/src/algorithms/tracking/OrthogonalTrackSeedingConfig.h b/src/algorithms/tracking/OrthogonalTrackSeedingConfig.h index 755e5fa17d..048be3973a 100644 --- a/src/algorithms/tracking/OrthogonalTrackSeedingConfig.h +++ b/src/algorithms/tracking/OrthogonalTrackSeedingConfig.h @@ -19,7 +19,7 @@ namespace eicrecon { float zMax = 1700. * Acts::UnitConstants::mm; // max z to look for hits to compose seeds float zMin = -1500. * Acts::UnitConstants::mm; // min z to look for hits to compose seeds float deltaRMinTopSP = 10. * Acts::UnitConstants::mm; // Min distance in r between middle and top SP in one seed - float deltaRMaxTopSP = 200. * Acts::UnitConstants::mm; // Max distance in r between middle and top SP in one seed + float deltaRMaxTopSP = 450. * Acts::UnitConstants::mm; // Max distance in r between middle and top SP in one seed float deltaRMinBottomSP = 10. * Acts::UnitConstants::mm; // Min distance in r between middle and bottom SP in one seed float deltaRMaxBottomSP = 200. * Acts::UnitConstants::mm; // Max distance in r between middle and bottom SP in one seed float collisionRegionMin = -250 * Acts::UnitConstants::mm; // Min z for primary vertex @@ -35,7 +35,6 @@ namespace eicrecon { float beamPosX = 0; // x offset for beam position float beamPosY = 0; // y offset for beam position float impactMax = 3. * Acts::UnitConstants::mm; // Maximum transverse PCA allowed - float bFieldMin = 0.1 * Acts::UnitConstants::T; // T (in Acts units of GeV/[e*mm]) - Minimum Magnetic field strength float rMinMiddle = 20. * Acts::UnitConstants::mm; // Middle spacepoint must fall between these two radii float rMaxMiddle = 400. * Acts::UnitConstants::mm; diff --git a/src/algorithms/tracking/SpacePoint.h b/src/algorithms/tracking/SpacePoint.h index c5401bc60a..b459279313 100644 --- a/src/algorithms/tracking/SpacePoint.h +++ b/src/algorithms/tracking/SpacePoint.h @@ -11,7 +11,7 @@ class SpacePoint : public edm4eic::TrackerHit public: const Acts::Surface *m_surface = nullptr; - SpacePoint(const TrackerHit& hit) : TrackerHit(hit) {} + SpacePoint(const TrackerHit& hit) : TrackerHit(hit) {} void setSurface(std::shared_ptr m_geoSvc) { @@ -36,6 +36,9 @@ class SpacePoint : public edm4eic::TrackerHit } float varianceZ() const { return getPositionError().zz; } + float t() const { return getTime(); } + float varianceT() const { return getTimeError(); } + bool isOnSurface() const { if (m_surface == nullptr) { return false; diff --git a/src/algorithms/tracking/TrackParamTruthInit.cc b/src/algorithms/tracking/TrackParamTruthInit.cc index bfd479d396..2430c839ed 100644 --- a/src/algorithms/tracking/TrackParamTruthInit.cc +++ b/src/algorithms/tracking/TrackParamTruthInit.cc @@ -10,7 +10,6 @@ #include #include -#include #include #include #include @@ -28,9 +27,6 @@ void eicrecon::TrackParamTruthInit::init(std::shared_ptr geo_svc, const std::shared_ptr logger) { m_log = logger; m_geoSvc = geo_svc; - - // TODO make a service? - m_pdg_db = std::make_shared(); } std::unique_ptr @@ -79,14 +75,9 @@ eicrecon::TrackParamTruthInit::produce(const edm4hep::MCParticleCollection* mcpa // get the particle charge // note that we cannot trust the mcparticles charge, as DD4hep // sets this value to zero! let's lookup by PDGID instead - //const double charge = m_pidSvc->particle(mcparticle.getPDG()).charge; const auto pdg = mcparticle.getPDG(); - const auto* particle = m_pdg_db->GetParticle(pdg); - if (particle == nullptr) { - m_log->debug("particle with PDG {} not in TDatabasePDG", pdg); - continue; - } - double charge = std::copysign(1.0,particle->Charge()); + const auto particle = m_particleSvc.particle(pdg); + double charge = std::copysign(1.0, particle.charge); if (abs(charge) < std::numeric_limits::epsilon()) { m_log->trace("ignoring neutral particle"); continue; diff --git a/src/algorithms/tracking/TrackParamTruthInit.h b/src/algorithms/tracking/TrackParamTruthInit.h index 7f9de2c0e9..0af8b5980a 100644 --- a/src/algorithms/tracking/TrackParamTruthInit.h +++ b/src/algorithms/tracking/TrackParamTruthInit.h @@ -4,7 +4,6 @@ #pragma once -#include #include #include #include @@ -13,6 +12,7 @@ #include "ActsGeometryProvider.h" #include "TrackParamTruthInitConfig.h" +#include "algorithms/interfaces/ParticleSvc.h" #include "algorithms/interfaces/WithPodConfig.h" namespace eicrecon { @@ -27,9 +27,10 @@ namespace eicrecon { private: std::shared_ptr m_log; - std::shared_ptr m_pdg_db; std::shared_ptr m_geoSvc; + const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance(); + std::default_random_engine generator; // TODO: need something more appropriate here std::uniform_int_distribution m_uniformIntDist{-1, 1}; // defaults to min=-1, max=1 std::normal_distribution m_normDist; diff --git a/src/algorithms/tracking/TrackProjector.cc b/src/algorithms/tracking/TrackProjector.cc index 68eb0859cf..fbf9d8fd47 100644 --- a/src/algorithms/tracking/TrackProjector.cc +++ b/src/algorithms/tracking/TrackProjector.cc @@ -80,24 +80,37 @@ namespace eicrecon { m_nCalibrated++; } - // get track state parameters and their covariances - const auto ¶meter = trackstate.predicted(); - const auto &covariance = trackstate.predictedCovariance(); + // get track state bound parameters and their boundCovs + const auto &boundParams = trackstate.predicted(); + const auto &boundCov = trackstate.predictedCovariance(); // convert local to global auto global = trackstate.referenceSurface().localToGlobal( m_geo_provider->getActsGeometryContext(), - {parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1]}, + {boundParams[Acts::eBoundLoc0], boundParams[Acts::eBoundLoc1]}, Acts::makeDirectionFromPhiTheta( - parameter[Acts::eBoundPhi], - parameter[Acts::eBoundTheta] + boundParams[Acts::eBoundPhi], + boundParams[Acts::eBoundTheta] ) ); + +#if Acts_VERSION_MAJOR >= 34 + auto freeParams = Acts::transformBoundToFreeParameters( + trackstate.referenceSurface(), + m_geo_provider->getActsGeometryContext(), + boundParams); + auto jacobian = trackstate.referenceSurface().boundToFreeJacobian( + m_geo_provider->getActsGeometryContext(), + freeParams.template segment<3>(Acts::eFreePos0), + freeParams.template segment<3>(Acts::eFreeDir0) + ); +#else auto jacobian = trackstate.referenceSurface().boundToFreeJacobian( m_geo_provider->getActsGeometryContext(), - parameter + boundParams ); - auto free_covariance = jacobian * covariance * jacobian.transpose(); +#endif + auto freeCov = jacobian * boundCov * jacobian.transpose(); // global position const decltype(edm4eic::TrackPoint::position) position{ @@ -108,45 +121,45 @@ namespace eicrecon { // local position const decltype(edm4eic::TrackParametersData::loc) loc{ - static_cast(parameter[Acts::eBoundLoc0]), - static_cast(parameter[Acts::eBoundLoc1]) + static_cast(boundParams[Acts::eBoundLoc0]), + static_cast(boundParams[Acts::eBoundLoc1]) }; const edm4eic::Cov2f locError{ - static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), - static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), - static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1)) + static_cast(boundCov(Acts::eBoundLoc0, Acts::eBoundLoc0)), + static_cast(boundCov(Acts::eBoundLoc1, Acts::eBoundLoc1)), + static_cast(boundCov(Acts::eBoundLoc0, Acts::eBoundLoc1)) }; const decltype(edm4eic::TrackPoint::positionError) positionError{ - static_cast(free_covariance(Acts::eFreePos0, Acts::eFreePos0)), - static_cast(free_covariance(Acts::eFreePos1, Acts::eFreePos1)), - static_cast(free_covariance(Acts::eFreePos2, Acts::eFreePos2)), - static_cast(free_covariance(Acts::eFreePos0, Acts::eFreePos1)), - static_cast(free_covariance(Acts::eFreePos0, Acts::eFreePos2)), - static_cast(free_covariance(Acts::eFreePos1, Acts::eFreePos2)), + static_cast(freeCov(Acts::eFreePos0, Acts::eFreePos0)), + static_cast(freeCov(Acts::eFreePos1, Acts::eFreePos1)), + static_cast(freeCov(Acts::eFreePos2, Acts::eFreePos2)), + static_cast(freeCov(Acts::eFreePos0, Acts::eFreePos1)), + static_cast(freeCov(Acts::eFreePos0, Acts::eFreePos2)), + static_cast(freeCov(Acts::eFreePos1, Acts::eFreePos2)), }; // momentum const decltype(edm4eic::TrackPoint::momentum) momentum = edm4hep::utils::sphericalToVector( - static_cast(1.0 / std::abs(parameter[Acts::eBoundQOverP])), - static_cast(parameter[Acts::eBoundTheta]), - static_cast(parameter[Acts::eBoundPhi]) + static_cast(1.0 / std::abs(boundParams[Acts::eBoundQOverP])), + static_cast(boundParams[Acts::eBoundTheta]), + static_cast(boundParams[Acts::eBoundPhi]) ); const decltype(edm4eic::TrackPoint::momentumError) momentumError{ - static_cast(covariance(Acts::eBoundTheta, Acts::eBoundTheta)), - static_cast(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), - static_cast(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)), - static_cast(covariance(Acts::eBoundTheta, Acts::eBoundPhi)), - static_cast(covariance(Acts::eBoundTheta, Acts::eBoundQOverP)), - static_cast(covariance(Acts::eBoundPhi, Acts::eBoundQOverP)) + static_cast(boundCov(Acts::eBoundTheta, Acts::eBoundTheta)), + static_cast(boundCov(Acts::eBoundPhi, Acts::eBoundPhi)), + static_cast(boundCov(Acts::eBoundQOverP, Acts::eBoundQOverP)), + static_cast(boundCov(Acts::eBoundTheta, Acts::eBoundPhi)), + static_cast(boundCov(Acts::eBoundTheta, Acts::eBoundQOverP)), + static_cast(boundCov(Acts::eBoundPhi, Acts::eBoundQOverP)) }; - const float time{static_cast(parameter(Acts::eBoundTime))}; - const float timeError{sqrt(static_cast(covariance(Acts::eBoundTime, Acts::eBoundTime)))}; - const float theta(parameter[Acts::eBoundTheta]); - const float phi(parameter[Acts::eBoundPhi]); + const float time{static_cast(boundParams(Acts::eBoundTime))}; + const float timeError{sqrt(static_cast(boundCov(Acts::eBoundTime, Acts::eBoundTime)))}; + const float theta(boundParams[Acts::eBoundTheta]); + const float phi(boundParams[Acts::eBoundPhi]); const decltype(edm4eic::TrackPoint::directionError) directionError{ - static_cast(covariance(Acts::eBoundTheta, Acts::eBoundTheta)), - static_cast(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), - static_cast(covariance(Acts::eBoundTheta, Acts::eBoundPhi)) + static_cast(boundCov(Acts::eBoundTheta, Acts::eBoundTheta)), + static_cast(boundCov(Acts::eBoundPhi, Acts::eBoundPhi)), + static_cast(boundCov(Acts::eBoundTheta, Acts::eBoundPhi)) }; const float pathLength = static_cast(trackstate.pathLength()); const float pathLengthError = 0; @@ -191,12 +204,12 @@ namespace eicrecon { m_log->debug(" ******************************"); // Local position on the reference surface. - //m_log->debug("parameter[eBoundLoc0] = {}", parameter[Acts::eBoundLoc0]); - //m_log->debug("parameter[eBoundLoc1] = {}", parameter[Acts::eBoundLoc1]); - //m_log->debug("parameter[eBoundPhi] = {}", parameter[Acts::eBoundPhi]); - //m_log->debug("parameter[eBoundTheta] = {}", parameter[Acts::eBoundTheta]); - //m_log->debug("parameter[eBoundQOverP] = {}", parameter[Acts::eBoundQOverP]); - //m_log->debug("parameter[eBoundTime] = {}", parameter[Acts::eBoundTime]); + //m_log->debug("boundParams[eBoundLoc0] = {}", boundParams[Acts::eBoundLoc0]); + //m_log->debug("boundParams[eBoundLoc1] = {}", boundParams[Acts::eBoundLoc1]); + //m_log->debug("boundParams[eBoundPhi] = {}", boundParams[Acts::eBoundPhi]); + //m_log->debug("boundParams[eBoundTheta] = {}", boundParams[Acts::eBoundTheta]); + //m_log->debug("boundParams[eBoundQOverP] = {}", boundParams[Acts::eBoundQOverP]); + //m_log->debug("boundParams[eBoundTime] = {}", boundParams[Acts::eBoundTime]); //m_log->debug("predicted variables: {}", trackstate.predicted()); }); diff --git a/src/algorithms/tracking/TrackProjector.h b/src/algorithms/tracking/TrackProjector.h index d1aa9fa5d2..2af5a74bb9 100644 --- a/src/algorithms/tracking/TrackProjector.h +++ b/src/algorithms/tracking/TrackProjector.h @@ -11,7 +11,6 @@ #include #include "ActsGeometryProvider.h" -#include "algorithms/interfaces/WithPodConfig.h" namespace eicrecon { diff --git a/src/algorithms/tracking/TrackSeeding.cc b/src/algorithms/tracking/TrackSeeding.cc index c108764ccd..c1094e78c8 100644 --- a/src/algorithms/tracking/TrackSeeding.cc +++ b/src/algorithms/tracking/TrackSeeding.cc @@ -87,6 +87,7 @@ void eicrecon::TrackSeeding::configure() { // Finder parameters m_seedFinderConfig.seedFilter = std::make_unique>(Acts::SeedFilter(m_seedFilterConfig)); m_seedFinderConfig.rMax = m_cfg.rMax; + m_seedFinderConfig.rMin = m_cfg.rMin; m_seedFinderConfig.deltaRMinTopSP = m_cfg.deltaRMinTopSP; m_seedFinderConfig.deltaRMaxTopSP = m_cfg.deltaRMaxTopSP; m_seedFinderConfig.deltaRMinBottomSP = m_cfg.deltaRMinBottomSP; @@ -120,6 +121,15 @@ std::unique_ptr eicrecon::TrackSeeding::prod Acts::SeedFinderOrthogonal finder(m_seedFinderConfig); // FIXME move into class scope +#if Acts_VERSION_MAJOR >= 32 + std::function>( + const eicrecon::SpacePoint *sp)> + create_coordinates = [](const eicrecon::SpacePoint *sp) { + Acts::Vector3 position(sp->x(), sp->y(), sp->z()); + Acts::Vector2 variance(sp->varianceR(), sp->varianceZ()); + return std::make_tuple(position, variance, sp->t()); + }; +#else std::function( const eicrecon::SpacePoint *sp)> create_coordinates = [](const eicrecon::SpacePoint *sp) { @@ -127,6 +137,7 @@ std::unique_ptr eicrecon::TrackSeeding::prod Acts::Vector2 variance(sp->varianceR(), sp->varianceZ()); return std::make_pair(position, variance); }; +#endif eicrecon::SeedContainer seeds = finder.createSeeds(m_seedFinderOptions, spacePoints, create_coordinates); diff --git a/src/algorithms/tracking/TrackerHitReconstruction.cc b/src/algorithms/tracking/TrackerHitReconstruction.cc index 155ba7d35f..5fe85534e4 100644 --- a/src/algorithms/tracking/TrackerHitReconstruction.cc +++ b/src/algorithms/tracking/TrackerHitReconstruction.cc @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -66,7 +67,10 @@ std::unique_ptr TrackerHitReconstruction::process // - XYZ segmentation: xx -> sigma_x, yy-> sigma_y, zz -> sigma_z, tt -> 0 // This is properly in line with how we get the local coordinates for the hit // in the TrackerSourceLinker. - rec_hits->create( + #if EDM4EIC_VERSION_MAJOR >= 7 + auto rec_hit = + #endif + rec_hits->create( raw_hit.getCellID(), // Raw DD4hep cell ID edm4hep::Vector3f{static_cast(pos.x() / mm), static_cast(pos.y() / mm), static_cast(pos.z() / mm)}, // mm edm4eic::CovDiag3f{get_variance(dim[0] / mm), get_variance(dim[1] / mm), // variance (see note above) @@ -75,6 +79,9 @@ std::unique_ptr TrackerHitReconstruction::process m_cfg.timeResolution, // in ns static_cast(raw_hit.getCharge() / 1.0e6), // Collected energy (GeV) 0.0F); // Error on the energy +#if EDM4EIC_VERSION_MAJOR >= 7 + rec_hit.setRawHit(raw_hit); +#endif } diff --git a/src/algorithms/tracking/TracksToParticles.cc b/src/algorithms/tracking/TracksToParticles.cc index 12399d2737..499f65f5b8 100644 --- a/src/algorithms/tracking/TracksToParticles.cc +++ b/src/algorithms/tracking/TracksToParticles.cc @@ -1,26 +1,20 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 - 2024, Sylvester Joosten, Wouter Deconinck, Dmitry Romanov, Christopher Dilks, Dmitry Kalinkin -#include #include #include -#include +#include #include #include #include #include #include #include -#include #include -#include #include -#include -#include #include #include "TracksToParticles.h" -#include "TracksToParticlesConfig.h" namespace eicrecon { @@ -30,13 +24,8 @@ namespace eicrecon { void TracksToParticles::process( const TracksToParticles::Input& input, const TracksToParticles::Output& output ) const { - const auto [mc_particles, tracks] = input; - auto [parts, assocs] = output; - - const double sinPhiOver2Tolerance = sin(0.5 * m_cfg.phiTolerance); - tracePhiToleranceOnce(sinPhiOver2Tolerance, m_cfg.phiTolerance); - - std::vector mc_prt_is_consumed(mc_particles->size(), false); // MCParticle is already consumed flag + const auto [tracks, track_assocs] = input; + auto [parts, part_assocs] = output; for (const auto &track: *tracks) { auto trajectory = track.getTrajectory(); @@ -46,155 +35,44 @@ namespace eicrecon { const auto charge_rec = std::copysign(1., trk.getQOverP()); - debug("Match: [id] [mom] [theta] [phi] [charge] [PID]"); - debug(" Track : {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<4}", - trk.getObjectID().index, edm4hep::utils::magnitude(mom), edm4hep::utils::anglePolar(mom), edm4hep::utils::angleAzimuthal(mom), charge_rec); - - // utility variables for matching - int best_match = -1; - double best_delta = std::numeric_limits::max(); - for (size_t ip = 0; ip < mc_particles->size(); ++ip) { - const auto &mc_part = (*mc_particles)[ip]; - const auto &p = mc_part.getMomentum(); - - trace(" MCParticle with id={:<4} mom={:<8.3f} charge={}", mc_part.getObjectID().index, - edm4hep::utils::magnitude(p), mc_part.getCharge()); - - // Check if used - if (mc_prt_is_consumed[ip]) { - trace(" Ignoring. Particle is already used"); - continue; - } - - // Check if non-primary - if (mc_part.getGeneratorStatus() > 1) { - trace(" Ignoring. GeneratorStatus > 1 => Non-primary particle"); - continue; - } - - // Check if neutral - if (mc_part.getCharge() == 0) { - trace(" Ignoring. Neutral particle"); - continue; - } - - // Check opposite charge - if (mc_part.getCharge() * charge_rec < 0) { - trace(" Ignoring. Opposite charge particle"); - continue; - } - - const auto p_mag = edm4hep::utils::magnitude(p); - const auto p_phi = edm4hep::utils::angleAzimuthal(p); - const auto p_eta = edm4hep::utils::eta(p); - const double dp_rel = std::abs((edm4hep::utils::magnitude(mom) - p_mag) / p_mag); - // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow - // for phi rollovers - const double dsphi = std::abs(sin(0.5 * (edm4hep::utils::angleAzimuthal(mom) - p_phi))); - const double deta = std::abs((edm4hep::utils::eta(mom) - p_eta)); - - bool is_matching = dp_rel < m_cfg.momentumRelativeTolerance && - deta < m_cfg.etaTolerance && - dsphi < sinPhiOver2Tolerance; - - // Matching kinematics with the static variables doesn't work at low angles and within beam divergence - // TODO - Maybe reconsider variables used or divide into regions - // Backward going - if ((p_eta < -5) && (edm4hep::utils::eta(mom) < -5)) { - is_matching = true; - } - // Forward going - if ((p_eta > 5) && (edm4hep::utils::eta(mom) > 5)) { - is_matching = true; - } + debug("Converting track: index={:<4} momentum={:<8.3f} theta={:<8.3f} phi={:<8.2f} charge={:<4}", + trk.getObjectID().index, edm4hep::utils::magnitude(mom), edm4hep::utils::anglePolar(mom), edm4hep::utils::angleAzimuthal(mom), charge_rec); - trace(" Decision: {} dp: {:.4f} < {} && d_eta: {:.6f} < {} && d_sin_phi: {:.4e} < {:.4e} ", - is_matching? "Matching":"Ignoring", - dp_rel, m_cfg.momentumRelativeTolerance, - deta, m_cfg.etaTolerance, - dsphi, sinPhiOver2Tolerance); - - if (is_matching) { - const double delta = - std::hypot(dp_rel / m_cfg.momentumRelativeTolerance, deta / m_cfg.etaTolerance, - dsphi / sinPhiOver2Tolerance); - if (delta < best_delta) { - best_match = ip; - best_delta = delta; - trace(" Is the best match now"); - } - } - } auto rec_part = parts->create(); rec_part.addToTracks(track); - auto referencePoint = rec_part.getReferencePoint(); - // float time = 0; - float mass = 0; - if (best_match >= 0) { - trace("Best match is found and is: {}", best_match); - mc_prt_is_consumed[best_match] = true; - const auto &best_mc_part = (*mc_particles)[best_match]; - referencePoint = { - static_cast(best_mc_part.getVertex().x), - static_cast(best_mc_part.getVertex().y), - static_cast(best_mc_part.getVertex().z)}; // @TODO: not sure if vertex/reference point makes sense here - // time = mcpart.getTime(); - mass = best_mc_part.getMass(); - } - - rec_part.setType(static_cast(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes - rec_part.setEnergy((float) std::hypot(edm4hep::utils::magnitude(mom), mass)); + rec_part.setType(0); + rec_part.setEnergy(edm4hep::utils::magnitude(mom)); rec_part.setMomentum(mom); - rec_part.setReferencePoint(referencePoint); rec_part.setCharge(charge_rec); - rec_part.setMass(mass); + rec_part.setMass(0.); rec_part.setGoodnessOfPID(0); // assume no PID until proven otherwise // rec_part.covMatrix() // @TODO: covariance matrix on 4-momentum - // Also write MC <--> truth particle association if match was found - if (best_match >= 0) { - auto rec_assoc = assocs->create(); - rec_assoc.setRecID(rec_part.getObjectID().index); - rec_assoc.setSimID((*mc_particles)[best_match].getObjectID().index); - rec_assoc.setWeight(1); - rec_assoc.setRec(rec_part); - auto sim = (*mc_particles)[best_match]; - rec_assoc.setSim(sim); - - if (level() <= algorithms::LogLevel::kDebug) { - - const auto &mcpart = (*mc_particles)[best_match]; - const auto &p = mcpart.getMomentum(); - const auto p_mag = edm4hep::utils::magnitude(p); - const auto p_phi = edm4hep::utils::angleAzimuthal(p); - const auto p_theta = edm4hep::utils::anglePolar(p); - debug(" MCPart: {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<6}", - mcpart.getObjectID().index, p_mag, p_theta, p_phi, mcpart.getCharge(), - mcpart.getPDG()); - - debug(" Assoc: id={} SimId={} RecId={}", - rec_assoc.getObjectID().index, rec_assoc.getSim().getObjectID().index, rec_assoc.getSim().getObjectID().index); - - trace(" Assoc PDGs: sim.PDG | rec.PDG | rec.particleIDUsed.PDG = {:^6} | {:^6} | {:^6}", - rec_assoc.getSim().getPDG(), - rec_assoc.getRec().getPDG(), - rec_assoc.getRec().getParticleIDUsed().isAvailable() ? rec_assoc.getRec().getParticleIDUsed().getPDG() : 0); - - + double max_weight = -1.; + for (auto track_assoc : *track_assocs) { + if (track_assoc.getRec() == track) { + trace("Found track association: index={} -> index={}, weight={}", + track_assoc.getRec().getObjectID().index, + track_assoc.getSim().getObjectID().index, + track_assoc.getWeight()); + auto part_assoc = part_assocs->create(); + part_assoc.setRec(rec_part); + part_assoc.setSim(track_assoc.getSim()); + part_assoc.setRecID(part_assoc.getRec().getObjectID().index); + part_assoc.setSimID(part_assoc.getSim().getObjectID().index); + part_assoc.setWeight(track_assoc.getWeight()); + + if (max_weight < track_assoc.getWeight()) { + max_weight = track_assoc.getWeight(); + edm4hep::Vector3f referencePoint = { + static_cast(track_assoc.getSim().getVertex().x), + static_cast(track_assoc.getSim().getVertex().y), + static_cast(track_assoc.getSim().getVertex().z)}; // @TODO: not sure if vertex/reference point makes sense here + rec_part.setReferencePoint(referencePoint); + } } } - else { - debug(" MCPart: Did not find a good match"); - } } } } - - void TracksToParticles::tracePhiToleranceOnce(const double sinPhiOver2Tolerance, double phiTolerance) const { - // This function is called once to print tolerances useful for tracing - static std::once_flag do_it_once; - std::call_once(do_it_once, [this, sinPhiOver2Tolerance, phiTolerance]() { - trace("m_cfg.phiTolerance: {:<8.4f} => sinPhiOver2Tolerance: {:<8.4f}", sinPhiOver2Tolerance, phiTolerance); - }); - } } diff --git a/src/algorithms/tracking/TracksToParticles.h b/src/algorithms/tracking/TracksToParticles.h index e8fd034869..cf0ffc9746 100644 --- a/src/algorithms/tracking/TracksToParticles.h +++ b/src/algorithms/tracking/TracksToParticles.h @@ -6,13 +6,13 @@ #include #include +#include #include #include -#include +#include #include #include -#include "TracksToParticlesConfig.h" #include "algorithms/interfaces/WithPodConfig.h" @@ -20,21 +20,24 @@ namespace eicrecon { using TracksToParticlesAlgorithm = algorithms::Algorithm< - algorithms::Input, - algorithms::Output + algorithms::Input< + edm4eic::TrackCollection, + std::optional + >, + algorithms::Output< + edm4eic::ReconstructedParticleCollection, + std::optional + > >; -class TracksToParticles : public TracksToParticlesAlgorithm, public WithPodConfig { +class TracksToParticles : public TracksToParticlesAlgorithm, public WithPodConfig { public: - TracksToParticles(std::string_view name) : TracksToParticlesAlgorithm{name, {"inputMCParticlesCollection", "inputTracksCollection"}, {"outputReconstructedParticlesCollection", "outputAssociationsCollection"}, "Converts track to particles with associations"} {}; + TracksToParticles(std::string_view name) : TracksToParticlesAlgorithm{name, {"inputTracksCollection", "inputTrackAssociationsCollection"}, {"outputReconstructedParticlesCollection", "outputAssociationsCollection"}, "Converts track to particles with associations"} {}; void init() final; void process(const Input&, const Output&) const final; -private: - - void tracePhiToleranceOnce(const double sinPhiOver2Tolerance, double phiTolerance) const; }; } diff --git a/src/algorithms/tracking/TracksToParticlesConfig.h b/src/algorithms/tracking/TracksToParticlesConfig.h deleted file mode 100644 index 59146e43db..0000000000 --- a/src/algorithms/tracking/TracksToParticlesConfig.h +++ /dev/null @@ -1,17 +0,0 @@ -// Created by Dmitry Romanov -// Subject to the terms in the LICENSE file found in the top-level directory. -// - -#pragma once - -namespace eicrecon { - - struct TracksToParticlesConfig { - double momentumRelativeTolerance = 100.0; /// Matching momentum effectively disabled - - double phiTolerance = 0.1; /// Matching phi tolerance [rad] - - double etaTolerance = 0.2; /// Matching eta tolerance - }; - -} // eicrecon diff --git a/src/benchmarks/reconstruction/femc_studies/CMakeLists.txt b/src/benchmarks/reconstruction/femc_studies/CMakeLists.txt index 784be6c6a1..bf5afcd6d5 100644 --- a/src/benchmarks/reconstruction/femc_studies/CMakeLists.txt +++ b/src/benchmarks/reconstruction/femc_studies/CMakeLists.txt @@ -17,12 +17,3 @@ plugin_add_acts(${PLUGIN_NAME}) # correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds # headers to the correct installation directory plugin_glob_all(${PLUGIN_NAME}) - -# Add include directories (same as target_include_directories but for both -# plugin and library) -plugin_include_directories(${PLUGIN_NAME} SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}) - -# Add libraries (same as target_include_directories but for both plugin and -# library) -plugin_link_libraries(${PLUGIN_NAME} ${ROOT_LIBRARIES} ActsCore - EDM4EIC::edm4eic) diff --git a/src/benchmarks/reconstruction/lfhcal_studies/CMakeLists.txt b/src/benchmarks/reconstruction/lfhcal_studies/CMakeLists.txt index 784be6c6a1..bf5afcd6d5 100644 --- a/src/benchmarks/reconstruction/lfhcal_studies/CMakeLists.txt +++ b/src/benchmarks/reconstruction/lfhcal_studies/CMakeLists.txt @@ -17,12 +17,3 @@ plugin_add_acts(${PLUGIN_NAME}) # correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds # headers to the correct installation directory plugin_glob_all(${PLUGIN_NAME}) - -# Add include directories (same as target_include_directories but for both -# plugin and library) -plugin_include_directories(${PLUGIN_NAME} SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}) - -# Add libraries (same as target_include_directories but for both plugin and -# library) -plugin_link_libraries(${PLUGIN_NAME} ${ROOT_LIBRARIES} ActsCore - EDM4EIC::edm4eic) diff --git a/src/benchmarks/reconstruction/tof_efficiency/README.md b/src/benchmarks/reconstruction/tof_efficiency/README.md index e7c4315265..337b69be5c 100644 --- a/src/benchmarks/reconstruction/tof_efficiency/README.md +++ b/src/benchmarks/reconstruction/tof_efficiency/README.md @@ -65,7 +65,7 @@ export DETECTOR_PATH="/path/to/dd4hep/epic/" export DETECTOR="epic_tracking_only" # This makes tracking output data and input MC particles to be written to the output --Ppodio:output_include_collections="ReconstructedParticles,TrackParameters,MCParticles" +-Ppodio:output_collections="ReconstructedParticles,TrackParameters,MCParticles" # There is a centralized file where plugins can save their histograms: -Phistsfile=/home/romanov/work/data/eicrecon_test/tracking_test_gun.ana.root @@ -82,7 +82,7 @@ defining a list of objects and the output file name: ```bash # This makes tracking output data and input MC particles to be written to the output --Ppodio:output_include_collections="ReconstructedParticles,TrackParameters,MCParticles" +-Ppodio:output_collections="ReconstructedParticles,TrackParameters,MCParticles" # This sets file path containing output tree -Ppodio:output_file=/home/romanov/work/data/eicrecon_test/tracking_test_gun.edm4eic.root diff --git a/src/benchmarks/reconstruction/tracking_efficiency/README.md b/src/benchmarks/reconstruction/tracking_efficiency/README.md index 59f3ec2646..aaac842178 100644 --- a/src/benchmarks/reconstruction/tracking_efficiency/README.md +++ b/src/benchmarks/reconstruction/tracking_efficiency/README.md @@ -67,7 +67,7 @@ export DETECTOR_PATH="/path/to/dd4hep/epic/" export DETECTOR="epic_tracking_only" # This makes tracking output data and input MC particles to be written to the output --Ppodio:output_include_collections="ReconstructedParticles,TrackParameters,MCParticles" +-Ppodio:output_collections="ReconstructedParticles,TrackParameters,MCParticles" # There is a centralized file where plugins can save their histograms: -Phistsfile=/home/romanov/work/data/eicrecon_test/tracking_test_gun.ana.root @@ -84,7 +84,7 @@ defining a list of objects and the output file name: ```bash # This makes tracking output data and input MC particles to be written to the output --Ppodio:output_include_collections="ReconstructedParticles,TrackParameters,MCParticles" +-Ppodio:output_collections="ReconstructedParticles,TrackParameters,MCParticles" # This sets file path containing output tree -Ppodio:output_file=/home/romanov/work/data/eicrecon_test/tracking_test_gun.edm4eic.root diff --git a/src/benchmarks/reconstruction/tracking_occupancy/CMakeLists.txt b/src/benchmarks/reconstruction/tracking_occupancy/CMakeLists.txt index 784be6c6a1..bf5afcd6d5 100644 --- a/src/benchmarks/reconstruction/tracking_occupancy/CMakeLists.txt +++ b/src/benchmarks/reconstruction/tracking_occupancy/CMakeLists.txt @@ -17,12 +17,3 @@ plugin_add_acts(${PLUGIN_NAME}) # correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds # headers to the correct installation directory plugin_glob_all(${PLUGIN_NAME}) - -# Add include directories (same as target_include_directories but for both -# plugin and library) -plugin_include_directories(${PLUGIN_NAME} SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}) - -# Add libraries (same as target_include_directories but for both plugin and -# library) -plugin_link_libraries(${PLUGIN_NAME} ${ROOT_LIBRARIES} ActsCore - EDM4EIC::edm4eic) diff --git a/src/detectors/B0TRK/B0TRK.cc b/src/detectors/B0TRK/B0TRK.cc index 213ecc5b91..ed97109140 100644 --- a/src/detectors/B0TRK/B0TRK.cc +++ b/src/detectors/B0TRK/B0TRK.cc @@ -26,7 +26,7 @@ void InitPlugin(JApplication *app) { }, { "B0TrackerRawHits", - "B0TrackerHitAssociations" + "B0TrackerRawHitAssociations" }, { .threshold = 10.0 * dd4hep::keV, diff --git a/src/detectors/BEMC/BEMC.cc b/src/detectors/BEMC/BEMC.cc index 10dbd0bee2..2d652a5ed6 100644 --- a/src/detectors/BEMC/BEMC.cc +++ b/src/detectors/BEMC/BEMC.cc @@ -97,6 +97,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed diff --git a/src/detectors/BHCAL/BHCAL.cc b/src/detectors/BHCAL/BHCAL.cc index db2015b467..f229eae805 100644 --- a/src/detectors/BHCAL/BHCAL.cc +++ b/src/detectors/BHCAL/BHCAL.cc @@ -1,15 +1,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 - 2024 David Lawrence, Derek Anderson, Wouter Deconinck -#include -#include -#include #include #include -#include -#include #include -#include #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h" #include "algorithms/calorimetry/CalorimeterIslandClusterConfig.h" @@ -19,36 +13,9 @@ #include "factories/calorimetry/CalorimeterHitReco_factory.h" #include "factories/calorimetry/CalorimeterIslandCluster_factory.h" #include "factories/calorimetry/CalorimeterTruthClustering_factory.h" -#include "services/geometry/dd4hep/DD4hep_service.h" extern "C" { - bool UseSectorIndexedBHCalReadout(JApplication *app) { - - using namespace eicrecon; - - // grab detector & segmentation descriptor - auto detector = app->GetService()->detector(); - dd4hep::IDDescriptor descriptor; - try { - descriptor = detector->readout("HcalBarrelHits").idSpec(); - } catch(const std::runtime_error &e) { - // detector not found in the geometry, result should not matter - - // default to non-deprecated behavior - return true; - } - - // check if sector field is present - bool useSectorIndex = false; - try { - auto sector = descriptor.field("sector"); - return true; - } catch(const std::runtime_error &e) { - return false; - } - - } // end 'UseSectorIndexedBHCalReadout(JApplication*)' - void InitPlugin(JApplication *app) { using namespace eicrecon; @@ -74,19 +41,6 @@ extern "C" { " ( (abs(eta_1 - eta_2) == 0) && (abs(phi_1 - phi_2) == (320 - 1)) )" ") == 1"; - // If using readout structure with sector segmentation, - // ensure adjacency matrix uses sector indices - if ( UseSectorIndexedBHCalReadout(app) ) { - HcalBarrel_adjacencyMatrix = - "(" - " abs(fmod(tower_1, 24) - fmod(tower_2, 24))" - " + min(" - " abs((sector_1 - sector_2) * (2 * 5) + (floor(tower_1 / 24) - floor(tower_2 / 24)) * 5 + fmod(tile_1, 5) - fmod(tile_2, 5))," - " (32 * 2 * 5) - abs((sector_1 - sector_2) * (2 * 5) + (floor(tower_1 / 24) - floor(tower_2 / 24)) * 5 + fmod(tile_1, 5) - fmod(tile_2, 5))" - " )" - ") == 1"; - } - app->Add(new JOmniFactoryGeneratorT( "HcalBarrelRawHits", {"HcalBarrelHits"}, {"HcalBarrelRawHits"}, { @@ -131,8 +85,6 @@ extern "C" { .adjacencyMatrix = HcalBarrel_adjacencyMatrix, .readout = "HcalBarrelHits", .sectorDist = 5.0 * dd4hep::cm, - .localDistXY = {15*dd4hep::mm, 15*dd4hep::mm}, - .dimScaledLocalDistXY = {50.0*dd4hep::mm, 50.0*dd4hep::mm}, .splitCluster = false, .minClusterHitEdep = 5.0 * dd4hep::MeV, .minClusterCenterEdep = 30.0 * dd4hep::MeV, diff --git a/src/detectors/BTOF/BTOF.cc b/src/detectors/BTOF/BTOF.cc index ef0ba38593..d40f5cdc4e 100644 --- a/src/detectors/BTOF/BTOF.cc +++ b/src/detectors/BTOF/BTOF.cc @@ -8,15 +8,13 @@ #include #include #include -#include #include #include #include #include -#include -#include #include "algorithms/interfaces/WithPodConfig.h" +#include "algorithms/pid_lut/PIDLookupConfig.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/digi/SiliconTrackerDigi_factory.h" #include "factories/tracking/TrackerHitReconstruction_factory.h" @@ -37,7 +35,7 @@ void InitPlugin(JApplication *app) { }, { "TOFBarrelRawHit", - "TOFBarrelHitAssociations" + "TOFBarrelRawHitAssociations" }, { .threshold = 6.0 * dd4hep::keV, @@ -64,32 +62,46 @@ void InitPlugin(JApplication *app) { } catch(const std::runtime_error&) { // Nothing } - for (auto qualifier : std::vector({"", "Seeded"})) { - app->Add(new JOmniFactoryGeneratorT( - fmt::format("CombinedTOF{}LUTPID", qualifier), - { - fmt::format("Reconstructed{}ChargedWithPFRICHPIDParticles", qualifier), - fmt::format("Reconstructed{}ChargedWithPFRICHPIDParticleAssociations", qualifier), - }, - { - fmt::format("Reconstructed{}ChargedWithPFRICHTOFPIDParticles", qualifier), - fmt::format("Reconstructed{}ChargedWithPFRICHTOFPIDParticleAssociations", qualifier), - fmt::format("CombinedTOF{}ParticleIDs", qualifier), - }, - { - .filename="calibrations/tof.lut", - .system=BarrelTOF_ID, - .pdg_values={11, 211, 321, 2212}, - .charge_values={1}, - .momentum_edges={0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3.0, 3.3, 3.6, 3.9, 4.2, 4.5, 4.8, 5.1, 5.4, 5.7, 6.0}, - .polar_edges={2.50, 10.95, 19.40, 27.85, 36.30, 44.75, 53.20, 61.65, 70.10, 78.55, 87.00, 95.45, 103.90, 112.35, 120.80, 129.25, 137.70, 146.15, 154.60}, - .azimuthal_binning={0., 360., 360.}, // lower, upper, step - .momentum_bin_centers_in_lut=true, - .polar_bin_centers_in_lut=true, - }, - app - )); - } + PIDLookupConfig pid_cfg { + .filename="calibrations/tof.lut", + .system=BarrelTOF_ID, + .pdg_values={11, 211, 321, 2212}, + .charge_values={1}, + .momentum_edges={0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3.0, 3.3, 3.6, 3.9, 4.2, 4.5, 4.8, 5.1, 5.4, 5.7, 6.0}, + .polar_edges={2.50, 10.95, 19.40, 27.85, 36.30, 44.75, 53.20, 61.65, 70.10, 78.55, 87.00, 95.45, 103.90, 112.35, 120.80, 129.25, 137.70, 146.15, 154.60}, + .azimuthal_binning={0., 360., 360.}, // lower, upper, step + .momentum_bin_centers_in_lut=true, + .polar_bin_centers_in_lut=true, + }; + + app->Add(new JOmniFactoryGeneratorT( + "CombinedTOFLUTPID", + { + "ReconstructedChargedWithPFRICHPIDParticles", + "ReconstructedChargedWithPFRICHPIDParticleAssociations", + }, + { + "ReconstructedChargedWithPFRICHTOFPIDParticles", + "ReconstructedChargedWithPFRICHTOFPIDParticleAssociations", + "CombinedTOFParticleIDs", + }, + pid_cfg, + app + )); + app->Add(new JOmniFactoryGeneratorT( + "CombinedTOFSeededLUTPID", + { + "ReconstructedSeededChargedWithPFRICHPIDParticles", + "ReconstructedSeededChargedWithPFRICHPIDParticleAssociations", + }, + { + "ReconstructedSeededChargedWithPFRICHTOFPIDParticles", + "ReconstructedSeededChargedWithPFRICHTOFPIDParticleAssociations", + "CombinedTOFSeededParticleIDs", + }, + pid_cfg, + app + )); } } // extern "C" diff --git a/src/detectors/BTRK/BTRK.cc b/src/detectors/BTRK/BTRK.cc index 5bc3463e55..85b14b505a 100644 --- a/src/detectors/BTRK/BTRK.cc +++ b/src/detectors/BTRK/BTRK.cc @@ -26,7 +26,7 @@ void InitPlugin(JApplication *app) { }, { "SiBarrelRawHits", - "SiBarrelHitAssociations" + "SiBarrelRawHitAssociations" }, { .threshold = 0.54 * dd4hep::keV, diff --git a/src/detectors/BVTX/BVTX.cc b/src/detectors/BVTX/BVTX.cc index 4129e083b9..799dda2ba4 100644 --- a/src/detectors/BVTX/BVTX.cc +++ b/src/detectors/BVTX/BVTX.cc @@ -26,7 +26,7 @@ void InitPlugin(JApplication *app) { }, { "SiBarrelVertexRawHits", - "SiBarrelVertexHitAssociations" + "SiBarrelVertexRawHitAssociations" }, { .threshold = 0.54 * dd4hep::keV, diff --git a/src/detectors/DIRC/DIRC.cc b/src/detectors/DIRC/DIRC.cc index 5d0a98567f..f33691803b 100644 --- a/src/detectors/DIRC/DIRC.cc +++ b/src/detectors/DIRC/DIRC.cc @@ -8,18 +8,17 @@ #include #include -#include #include #include #include #include #include -#include #include #include // algorithm configurations #include "algorithms/digi/PhotoMultiplierHitDigiConfig.h" +#include "algorithms/pid_lut/PIDLookupConfig.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" // factories #include "global/digi/PhotoMultiplierHitDigi_factory.h" @@ -91,29 +90,44 @@ extern "C" { } catch(const std::runtime_error&) { // Nothing } - for (auto qualifier : std::vector({"", "Seeded"})) { - app->Add(new JOmniFactoryGeneratorT( - fmt::format("DIRC{}LUTPID", qualifier), - { - fmt::format("Reconstructed{}ChargedWithPFRICHTOFPIDParticles", qualifier), - fmt::format("Reconstructed{}ChargedWithPFRICHTOFPIDParticleAssociations", qualifier), - }, - { - fmt::format("Reconstructed{}ChargedWithPFRICHTOFDIRCPIDParticles", qualifier), - fmt::format("Reconstructed{}ChargedWithPFRICHTOFDIRCPIDParticleAssociations", qualifier), - fmt::format("DIRC{}ParticleIDs", qualifier), - }, - { - .filename="calibrations/hpdirc.lut.gz", - .system=BarrelDIRC_ID, - .pdg_values={11, 211, 321, 2212}, - .charge_values={-1, 1}, - .momentum_edges={0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2, 7.4, 7.6, 7.8, 8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.8, 10.0, 10.2}, - .polar_edges={25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0, 40.0, 41.0, 42.0, 43.0, 44.0, 45.0, 46.0, 47.0, 48.0, 49.0, 50.0, 51.0, 52.0, 53.0, 54.0, 55.0, 56.0, 57.0, 58.0, 59.0, 60.0, 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 67.0, 68.0, 69.0, 70.0, 71.0, 72.0, 73.0, 74.0, 75.0, 76.0, 77.0, 78.0, 79.0, 80.0, 81.0, 82.0, 83.0, 84.0, 85.0, 86.0, 87.0, 88.0, 89.0, 90.0, 91.0, 92.0, 93.0, 94.0, 95.0, 96.0, 97.0, 98.0, 99.0, 100.0, 101.0, 102.0, 103.0, 104.0, 105.0, 106.0, 107.0, 108.0, 109.0, 110.0, 111.0, 112.0, 113.0, 114.0, 115.0, 116.0, 117.0, 118.0, 119.0, 120.0, 121.0, 122.0, 123.0, 124.0, 125.0, 126.0, 127.0, 128.0, 129.0, 130.0, 131.0, 132.0, 133.0, 134.0, 135.0, 136.0, 137.0, 138.0, 139.0, 140.0, 141.0, 142.0, 143.0, 144.0, 145.0, 146.0, 147.0, 148.0, 149.0, 150.0, 151.0, 152.0, 153.0, 154.0, 155.0, 156.0, 157.0, 158.0, 159.0, 160.0}, - .azimuthal_binning={0.0, 30.5, 0.5}, // lower, upper, step - }, - app - )); - } + PIDLookupConfig pid_cfg { + .filename="calibrations/hpdirc.lut.gz", + .system=BarrelDIRC_ID, + .pdg_values={11, 211, 321, 2212}, + .charge_values={-1, 1}, + .momentum_edges={0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2, 7.4, 7.6, 7.8, 8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.8, 10.0, 10.2}, + .polar_edges={25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0, 40.0, 41.0, 42.0, 43.0, 44.0, 45.0, 46.0, 47.0, 48.0, 49.0, 50.0, 51.0, 52.0, 53.0, 54.0, 55.0, 56.0, 57.0, 58.0, 59.0, 60.0, 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 67.0, 68.0, 69.0, 70.0, 71.0, 72.0, 73.0, 74.0, 75.0, 76.0, 77.0, 78.0, 79.0, 80.0, 81.0, 82.0, 83.0, 84.0, 85.0, 86.0, 87.0, 88.0, 89.0, 90.0, 91.0, 92.0, 93.0, 94.0, 95.0, 96.0, 97.0, 98.0, 99.0, 100.0, 101.0, 102.0, 103.0, 104.0, 105.0, 106.0, 107.0, 108.0, 109.0, 110.0, 111.0, 112.0, 113.0, 114.0, 115.0, 116.0, 117.0, 118.0, 119.0, 120.0, 121.0, 122.0, 123.0, 124.0, 125.0, 126.0, 127.0, 128.0, 129.0, 130.0, 131.0, 132.0, 133.0, 134.0, 135.0, 136.0, 137.0, 138.0, 139.0, 140.0, 141.0, 142.0, 143.0, 144.0, 145.0, 146.0, 147.0, 148.0, 149.0, 150.0, 151.0, 152.0, 153.0, 154.0, 155.0, 156.0, 157.0, 158.0, 159.0, 160.0}, + .azimuthal_binning={0.0, 30.5, 0.5}, // lower, upper, step + }; + + app->Add(new JOmniFactoryGeneratorT( + "DIRCLUTPID", + { + "ReconstructedChargedWithPFRICHTOFPIDParticles", + "ReconstructedChargedWithPFRICHTOFPIDParticleAssociations", + }, + { + "ReconstructedChargedWithPFRICHTOFDIRCPIDParticles", + "ReconstructedChargedWithPFRICHTOFDIRCPIDParticleAssociations", + "DIRCParticleIDs", + }, + pid_cfg, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "DIRCSeededLUTPID", + { + "ReconstructedSeededChargedWithPFRICHTOFPIDParticles", + "ReconstructedSeededChargedWithPFRICHTOFPIDParticleAssociations", + }, + { + "ReconstructedSeededChargedWithPFRICHTOFDIRCPIDParticles", + "ReconstructedSeededChargedWithPFRICHTOFDIRCPIDParticleAssociations", + "DIRCSeededParticleIDs", + }, + pid_cfg, + app + )); } } diff --git a/src/detectors/DRICH/DRICH.cc b/src/detectors/DRICH/DRICH.cc index c7ce56daa9..69c0a288c0 100644 --- a/src/detectors/DRICH/DRICH.cc +++ b/src/detectors/DRICH/DRICH.cc @@ -8,7 +8,6 @@ #include #include #include -#include #include #include #include @@ -23,6 +22,7 @@ #include "algorithms/digi/PhotoMultiplierHitDigiConfig.h" #include "algorithms/pid/IrtCherenkovParticleIDConfig.h" #include "algorithms/pid/MergeParticleIDConfig.h" +#include "algorithms/pid_lut/PIDLookupConfig.h" #include "algorithms/tracking/TrackPropagationConfig.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" // factories @@ -195,33 +195,48 @@ extern "C" { } catch(const std::runtime_error&) { // Nothing } - for (auto qualifier : std::vector({"", "Seeded"})) { - app->Add(new JOmniFactoryGeneratorT( - fmt::format("DRICH{}LUTPID", qualifier), - { - fmt::format("Reconstructed{}ChargedWithPFRICHTOFDIRCPIDParticles", qualifier), - fmt::format("Reconstructed{}ChargedWithPFRICHTOFDIRCPIDParticleAssociations", qualifier), - }, - { - fmt::format("Reconstructed{}ChargedParticles", qualifier), - fmt::format("Reconstructed{}ChargedParticleAssociations", qualifier), - fmt::format("DRICH{}ParticleIDs", qualifier), - }, - { - .filename="calibrations/drich.lut", - .system=ForwardRICH_ID, - .pdg_values={211, 321, 2212}, - .charge_values={1}, - .momentum_edges={0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75, 10.25, 10.75, 11.25, 11.75, 12.25, 12.75, 13.25, 13.75, 14.25, 14.75, 15.25, 15.75, 16.25, 16.75, 17.25, 17.75, 18.25, 18.75, 19.25, 19.75, 20.50, 21.50, 22.50, 23.50, 24.50, 25.50, 26.50, 27.50, 28.50, 29.50, 30.50}, - .polar_edges={0.060, 0.164, 0.269, 0.439}, - .azimuthal_binning={0., 2 * M_PI, 2 * M_PI}, // lower, upper, step - .polar_bin_centers_in_lut=true, - .use_radians=true, - .missing_electron_prob=true, - }, - app - )); - } + PIDLookupConfig pid_cfg { + .filename="calibrations/drich.lut", + .system=ForwardRICH_ID, + .pdg_values={211, 321, 2212}, + .charge_values={1}, + .momentum_edges={0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75, 10.25, 10.75, 11.25, 11.75, 12.25, 12.75, 13.25, 13.75, 14.25, 14.75, 15.25, 15.75, 16.25, 16.75, 17.25, 17.75, 18.25, 18.75, 19.25, 19.75, 20.50, 21.50, 22.50, 23.50, 24.50, 25.50, 26.50, 27.50, 28.50, 29.50, 30.50}, + .polar_edges={0.060, 0.164, 0.269, 0.439}, + .azimuthal_binning={0., 2 * M_PI, 2 * M_PI}, // lower, upper, step + .polar_bin_centers_in_lut=true, + .use_radians=true, + .missing_electron_prob=true, + }; + + app->Add(new JOmniFactoryGeneratorT( + "DRICHLUTPID", + { + "ReconstructedChargedWithPFRICHTOFDIRCPIDParticles", + "ReconstructedChargedWithPFRICHTOFDIRCPIDParticleAssociations", + }, + { + "ReconstructedChargedParticles", + "ReconstructedChargedParticleAssociations", + "DRICHParticleIDs", + }, + pid_cfg, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "DRICHSeededLUTPID", + { + "ReconstructedSeededChargedWithPFRICHTOFDIRCPIDParticles", + "ReconstructedSeededChargedWithPFRICHTOFDIRCPIDParticleAssociations", + }, + { + "ReconstructedSeededChargedParticles", + "ReconstructedSeededChargedParticleAssociations", + "DRICHSeededParticleIDs", + }, + pid_cfg, + app + )); // clang-format on } } diff --git a/src/detectors/ECTOF/ECTOF.cc b/src/detectors/ECTOF/ECTOF.cc index 06d04a9e0f..e6812fc971 100644 --- a/src/detectors/ECTOF/ECTOF.cc +++ b/src/detectors/ECTOF/ECTOF.cc @@ -25,7 +25,7 @@ void InitPlugin(JApplication *app) { "TOFEndcapHits"}, { "TOFEndcapRawHits", - "TOFEndcapHitAssociations" + "TOFEndcapRawHitAssociations" }, { .threshold = 6.0 * dd4hep::keV, diff --git a/src/detectors/ECTRK/ECTRK.cc b/src/detectors/ECTRK/ECTRK.cc index 2ed05da2b7..c8ecf3b969 100644 --- a/src/detectors/ECTRK/ECTRK.cc +++ b/src/detectors/ECTRK/ECTRK.cc @@ -26,7 +26,7 @@ void InitPlugin(JApplication *app) { }, { "SiEndcapTrackerRawHits", - "SiEndcapHitAssociations" + "SiEndcapTrackerRawHitAssociations" }, { .threshold = 0.54 * dd4hep::keV, diff --git a/src/detectors/EEMC/EEMC.cc b/src/detectors/EEMC/EEMC.cc index ec3a04b669..a9afa0430c 100644 --- a/src/detectors/EEMC/EEMC.cc +++ b/src/detectors/EEMC/EEMC.cc @@ -69,6 +69,7 @@ extern "C" { "EcalEndcapNIslandProtoClusters", {"EcalEndcapNRecHits"}, {"EcalEndcapNIslandProtoClusters"}, { .adjacencyMatrix = "(abs(row_1 - row_2) + abs(column_1 - column_2)) == 1", + .peakNeighbourhoodMatrix = "max(abs(row_1 - row_2), abs(column_1 - column_2)) == 1", .readout = "EcalEndcapNHits", .sectorDist = 5.0 * dd4hep::cm, .splitCluster = true, diff --git a/src/detectors/FEMC/FEMC.cc b/src/detectors/FEMC/FEMC.cc index 606dd54bd0..d613781b5d 100644 --- a/src/detectors/FEMC/FEMC.cc +++ b/src/detectors/FEMC/FEMC.cc @@ -70,7 +70,6 @@ extern "C" { { .sectorDist = 5.0 * dd4hep::cm, .localDistXY = {10.0 * dd4hep::cm, 10.0 * dd4hep::cm}, - .dimScaledLocalDistXY = {1.5, 1.5}, .splitCluster = true, .minClusterHitEdep = 0.0 * dd4hep::MeV, .minClusterCenterEdep = 10.0 * dd4hep::MeV, diff --git a/src/detectors/FHCAL/FHCAL.cc b/src/detectors/FHCAL/FHCAL.cc index e58e146f5d..0d78240fce 100644 --- a/src/detectors/FHCAL/FHCAL.cc +++ b/src/detectors/FHCAL/FHCAL.cc @@ -8,6 +8,7 @@ #include #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h" +#include "algorithms/calorimetry/ImagingTopoClusterConfig.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h" #include "factories/calorimetry/CalorimeterHitDigi_factory.h" @@ -58,7 +59,7 @@ extern "C" { .thresholdFactor = 0., .thresholdValue = 41.0, // 0.25 MeV --> 0.25 / 200 * 32768 = 41 - .sampFrac = "0.0098", + .sampFrac = "1.0", .readout = "HcalEndcapPInsertHits", .layerField="layer", }, @@ -82,8 +83,8 @@ extern "C" { "HcalEndcapPInsertSubcellHits", {"HcalEndcapPInsertRecHits"}, {"HcalEndcapPInsertSubcellHits"}, { .MIP = 800. * dd4hep::keV, - .Emin_in_MIPs=0.1, - .tmax=50 * dd4hep::ns, + .Emin_in_MIPs=0.5, + .tmax=162 * dd4hep::ns, //150 ns + (z at front face)/(speed of light) }, app // TODO: Remove me once fixed )); @@ -94,12 +95,13 @@ extern "C" { { .neighbourLayersRange = 1, .localDistXY = {0.76*side_length, 0.76*side_length*sin(M_PI/3)}, - .layerDistEtaPhi = {17e-3, 18.1e-3}, + .layerDistXY = {0.76*side_length, 0.76*side_length*sin(M_PI/3)}, + .layerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xy, .sectorDist = 10.0 * dd4hep::cm, - .minClusterHitEdep = 100.0 * dd4hep::keV, - .minClusterCenterEdep = 11.0 * dd4hep::MeV, + .minClusterHitEdep = 5.0 * dd4hep::keV, + .minClusterCenterEdep = 5.0 * dd4hep::MeV, .minClusterEdep = 11.0 * dd4hep::MeV, - .minClusterNhits = 10, + .minClusterNhits = 100, }, app // TODO: Remove me once fixed )); @@ -113,8 +115,9 @@ extern "C" { "HcalEndcapPInsertTruthClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation { .energyWeight = "log", - .sampFrac = 1.0, + .sampFrac = 0.0257, .logWeightBase = 3.6, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = true }, app // TODO: Remove me once fixed @@ -130,8 +133,9 @@ extern "C" { "HcalEndcapPInsertClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation { .energyWeight = "log", - .sampFrac = 1.0, + .sampFrac = 0.0257, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed @@ -186,10 +190,11 @@ extern "C" { // Magic constants: // 54 - number of modules in a row/column // 2 - number of towers in a module - std::string cellIdx_1 = "(54*2-moduleIDx_1*2+towerx_1)"; - std::string cellIdx_2 = "(54*2-moduleIDx_2*2+towerx_2)"; - std::string cellIdy_1 = "(54*2-moduleIDy_1*2+towery_1)"; - std::string cellIdy_2 = "(54*2-moduleIDy_2*2+towery_2)"; + // sign for towerx and towery are *negative* to maintain linearity with global X and Y + std::string cellIdx_1 = "(54*2-moduleIDx_1*2-towerx_1)"; + std::string cellIdx_2 = "(54*2-moduleIDx_2*2-towerx_2)"; + std::string cellIdy_1 = "(54*2-moduleIDy_1*2-towery_1)"; + std::string cellIdy_2 = "(54*2-moduleIDy_2*2-towery_2)"; std::string cellIdz_1 = "rlayerz_1"; std::string cellIdz_2 = "rlayerz_2"; std::string deltaX = Form("abs(%s-%s)", cellIdx_2.data(), cellIdx_1.data()); @@ -227,6 +232,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 4.5, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed @@ -244,6 +250,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 4.5, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed diff --git a/src/detectors/FOFFMTRK/FOFFMTRK.cc b/src/detectors/FOFFMTRK/FOFFMTRK.cc index 24a37ee54c..8fe4d34323 100644 --- a/src/detectors/FOFFMTRK/FOFFMTRK.cc +++ b/src/detectors/FOFFMTRK/FOFFMTRK.cc @@ -29,7 +29,7 @@ void InitPlugin(JApplication *app) { }, { "ForwardOffMTrackerRawHits", - "ForwardOffMTrackerHitAssociations" + "ForwardOffMTrackerRawHitAssociations" }, { .threshold = 10.0 * dd4hep::keV, diff --git a/src/detectors/LOWQ2/LOWQ2.cc b/src/detectors/LOWQ2/LOWQ2.cc index 34e7a19676..1b5feb86fc 100644 --- a/src/detectors/LOWQ2/LOWQ2.cc +++ b/src/detectors/LOWQ2/LOWQ2.cc @@ -42,7 +42,7 @@ extern "C" { }, { "TaggerTrackerRawHits", - "TaggerTrackerHitAssociations" + "TaggerTrackerRawHitAssociations" }, { .threshold = 1.5 * dd4hep::keV, diff --git a/src/detectors/LUMISPECCAL/README.md b/src/detectors/LUMISPECCAL/README.md index 59f3ec2646..aaac842178 100644 --- a/src/detectors/LUMISPECCAL/README.md +++ b/src/detectors/LUMISPECCAL/README.md @@ -67,7 +67,7 @@ export DETECTOR_PATH="/path/to/dd4hep/epic/" export DETECTOR="epic_tracking_only" # This makes tracking output data and input MC particles to be written to the output --Ppodio:output_include_collections="ReconstructedParticles,TrackParameters,MCParticles" +-Ppodio:output_collections="ReconstructedParticles,TrackParameters,MCParticles" # There is a centralized file where plugins can save their histograms: -Phistsfile=/home/romanov/work/data/eicrecon_test/tracking_test_gun.ana.root @@ -84,7 +84,7 @@ defining a list of objects and the output file name: ```bash # This makes tracking output data and input MC particles to be written to the output --Ppodio:output_include_collections="ReconstructedParticles,TrackParameters,MCParticles" +-Ppodio:output_collections="ReconstructedParticles,TrackParameters,MCParticles" # This sets file path containing output tree -Ppodio:output_file=/home/romanov/work/data/eicrecon_test/tracking_test_gun.edm4eic.root diff --git a/src/detectors/MPGD/MPGD.cc b/src/detectors/MPGD/MPGD.cc index e2f4c6fdb4..d4fd16a7e6 100644 --- a/src/detectors/MPGD/MPGD.cc +++ b/src/detectors/MPGD/MPGD.cc @@ -26,7 +26,7 @@ void InitPlugin(JApplication *app) { }, { "MPGDBarrelRawHits", - "MPGDBarrelHitAssociations" + "MPGDBarrelRawHitAssociations" }, { .threshold = 0.25 * dd4hep::keV, @@ -54,7 +54,7 @@ void InitPlugin(JApplication *app) { }, { "OuterMPGDBarrelRawHits", - "OuterMPGDBarrelHitAssociations" + "OuterMPGDBarrelRawHitAssociations" }, { .threshold = 0.25 * dd4hep::keV, @@ -82,7 +82,7 @@ void InitPlugin(JApplication *app) { }, { "BackwardMPGDEndcapRawHits", - "BackwardMPGDEndcapAssociations" + "BackwardMPGDEndcapRawHitAssociations" }, { .threshold = 0.25 * dd4hep::keV, @@ -110,7 +110,7 @@ void InitPlugin(JApplication *app) { }, { "ForwardMPGDEndcapRawHits", - "ForwardMPGDHitAssociations" + "ForwardMPGDEndcapRawHitAssociations" }, { .threshold = 0.25 * dd4hep::keV, diff --git a/src/detectors/PFRICH/PFRICH.cc b/src/detectors/PFRICH/PFRICH.cc index 3a414b24e8..7dd889d2e2 100644 --- a/src/detectors/PFRICH/PFRICH.cc +++ b/src/detectors/PFRICH/PFRICH.cc @@ -7,18 +7,17 @@ #include #include #include -#include #include #include #include #include #include -#include #include #include // algorithm configurations #include "algorithms/digi/PhotoMultiplierHitDigiConfig.h" +#include "algorithms/pid_lut/PIDLookupConfig.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" // factories #include "global/digi/PhotoMultiplierHitDigi_factory.h" @@ -85,33 +84,48 @@ extern "C" { } catch(const std::runtime_error&) { // Nothing } - for (auto qualifier : std::vector({"", "Seeded"})) { - app->Add(new JOmniFactoryGeneratorT( - fmt::format("RICHEndcapN{}LUTPID", qualifier), - { - fmt::format("Reconstructed{}ChargedWithoutPIDParticles", qualifier), - fmt::format("Reconstructed{}ChargedWithoutPIDParticleAssociations", qualifier), - }, - { - fmt::format("Reconstructed{}ChargedWithPFRICHPIDParticles", qualifier), - fmt::format("Reconstructed{}ChargedWithPFRICHPIDParticleAssociations", qualifier), - fmt::format("RICHEndcapN{}ParticleIDs", qualifier), - }, - { - .filename="calibrations/pfrich.lut", - .system=BackwardRICH_ID, - .pdg_values={11, 211, 321, 2212}, - .charge_values={1}, - .momentum_edges={0.4, 0.8, 1.2, 1.6, 2, 2.4, 2.8, 3.2, 3.6, 4, 4.4, 4.8, 5.2, 5.6, 6, 6.4, 6.8, 7.2, 7.6, 8, 8.4, 8.8, 9.2, 9.6, 10, 10.4, 10.8, 11.2, 11.6, 12, 12.4, 12.8, 13.2, 13.6, 14, 14.4, 14.8, 15.2}, - .polar_edges={2.65, 2.6725, 2.695, 2.7175, 2.74, 2.7625, 2.785, 2.8075, 2.83, 2.8525, 2.875, 2.8975, 2.92, 2.9425, 2.965, 2.9875, 3.01, 3.0325, 3.055, 3.0775}, - .azimuthal_binning={0., 2 * M_PI, 2 * M_PI / 120.}, // lower, upper, step - .azimuthal_bin_centers_in_lut=true, - .momentum_bin_centers_in_lut=true, - .polar_bin_centers_in_lut=true, - .use_radians=true, - }, - app - )); - } + PIDLookupConfig pid_cfg { + .filename="calibrations/pfrich.lut", + .system=BackwardRICH_ID, + .pdg_values={11, 211, 321, 2212}, + .charge_values={1}, + .momentum_edges={0.4, 0.8, 1.2, 1.6, 2, 2.4, 2.8, 3.2, 3.6, 4, 4.4, 4.8, 5.2, 5.6, 6, 6.4, 6.8, 7.2, 7.6, 8, 8.4, 8.8, 9.2, 9.6, 10, 10.4, 10.8, 11.2, 11.6, 12, 12.4, 12.8, 13.2, 13.6, 14, 14.4, 14.8, 15.2}, + .polar_edges={2.65, 2.6725, 2.695, 2.7175, 2.74, 2.7625, 2.785, 2.8075, 2.83, 2.8525, 2.875, 2.8975, 2.92, 2.9425, 2.965, 2.9875, 3.01, 3.0325, 3.055, 3.0775}, + .azimuthal_binning={0., 2 * M_PI, 2 * M_PI / 120.}, // lower, upper, step + .azimuthal_bin_centers_in_lut=true, + .momentum_bin_centers_in_lut=true, + .polar_bin_centers_in_lut=true, + .use_radians=true, + }; + + app->Add(new JOmniFactoryGeneratorT( + "RICHEndcapNLUTPID", + { + "ReconstructedChargedWithoutPIDParticles", + "ReconstructedChargedWithoutPIDParticleAssociations", + }, + { + "ReconstructedChargedWithPFRICHPIDParticles", + "ReconstructedChargedWithPFRICHPIDParticleAssociations", + "RICHEndcapNParticleIDs", + }, + pid_cfg, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "RICHEndcapNSeededLUTPID", + { + "ReconstructedSeededChargedWithoutPIDParticles", + "ReconstructedSeededChargedWithoutPIDParticleAssociations", + }, + { + "ReconstructedSeededChargedWithPFRICHPIDParticles", + "ReconstructedSeededChargedWithPFRICHPIDParticleAssociations", + "RICHEndcapNSeededParticleIDs", + }, + pid_cfg, + app + )); } } diff --git a/src/detectors/RPOTS/RPOTS.cc b/src/detectors/RPOTS/RPOTS.cc index cf76ce7781..c374f691dc 100644 --- a/src/detectors/RPOTS/RPOTS.cc +++ b/src/detectors/RPOTS/RPOTS.cc @@ -29,7 +29,7 @@ void InitPlugin(JApplication *app) { }, { "ForwardRomanPotRawHits", - "ForwardRomanPotHitAssociations" + "ForwardRomanPotRawHitAssociations" }, { .threshold = 10.0 * dd4hep::keV, diff --git a/src/detectors/ZDC/ZDC.cc b/src/detectors/ZDC/ZDC.cc index a54235dff8..9a49cd08d8 100644 --- a/src/detectors/ZDC/ZDC.cc +++ b/src/detectors/ZDC/ZDC.cc @@ -8,7 +8,7 @@ #include #include -#include "algorithms/interfaces/WithPodConfig.h" +#include "algorithms/calorimetry/ImagingTopoClusterConfig.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h" #include "factories/calorimetry/CalorimeterHitDigi_factory.h" @@ -17,6 +17,7 @@ #include "factories/calorimetry/CalorimeterTruthClustering_factory.h" #include "factories/calorimetry/HEXPLIT_factory.h" #include "factories/calorimetry/ImagingTopoCluster_factory.h" + extern "C" { void InitPlugin(JApplication *app) { @@ -63,7 +64,6 @@ extern "C" { { .sectorDist = 5.0 * dd4hep::cm, .localDistXY = {50 * dd4hep::cm, 50 * dd4hep::cm}, - .dimScaledLocalDistXY = {50.0*dd4hep::mm, 50.0*dd4hep::mm}, .splitCluster = true, .minClusterHitEdep = 0.1 * dd4hep::MeV, .minClusterCenterEdep = 3.0 * dd4hep::MeV, @@ -84,6 +84,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 3.6, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed @@ -101,6 +102,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed @@ -144,8 +146,8 @@ extern "C" { "HcalFarForwardZDCSubcellHits", {"HcalFarForwardZDCRecHits"}, {"HcalFarForwardZDCSubcellHits"}, { .MIP = 472. * dd4hep::keV, - .Emin_in_MIPs=0.1, - .tmax=320 * dd4hep::ns, + .Emin_in_MIPs=0.5, + .tmax=269 * dd4hep::ns, }, app // TODO: Remove me once fixed )); @@ -155,13 +157,14 @@ extern "C" { "HcalFarForwardZDCImagingProtoClusters", {"HcalFarForwardZDCSubcellHits"}, {"HcalFarForwardZDCImagingProtoClusters"}, { .neighbourLayersRange = 1, - .localDistXY = {0.76*side_length, 0.76*side_length*sin(M_PI/3)}, - .layerDistEtaPhi = {17e-3, 18.1e-3}, + .localDistXY = {0.5*side_length, 0.5*side_length*sin(M_PI/3)}, + .layerDistXY = {0.25*side_length, 0.25*side_length*sin(M_PI/3)}, + .layerMode=eicrecon::ImagingTopoClusterConfig::ELayerMode::xy, .sectorDist = 10.0 * dd4hep::cm, .minClusterHitEdep = 100.0 * dd4hep::keV, - .minClusterCenterEdep = 1.0 * dd4hep::MeV, + .minClusterCenterEdep = 3.0 * dd4hep::MeV, .minClusterEdep = 11.0 * dd4hep::MeV, - .minClusterNhits = 10, + .minClusterNhits = 100, }, app // TODO: Remove me once fixed )); @@ -191,6 +194,7 @@ extern "C" { .sampFrac = 0.0203, .logWeightBaseCoeffs={5.0,0.65,0.31}, .logWeightBase_Eref=50*dd4hep::GeV, + .longitudinalShowerInfoAvailable = true, }, app // TODO: Remove me once fixed )); @@ -206,7 +210,6 @@ extern "C" { { .sectorDist = 5.0 * dd4hep::cm, .localDistXY = {50 * dd4hep::cm, 50 * dd4hep::cm}, - .dimScaledLocalDistXY = {50.0*dd4hep::mm, 50.0*dd4hep::mm}, .splitCluster = true, .minClusterHitEdep = 0.1 * dd4hep::MeV, .minClusterCenterEdep = 3.0 * dd4hep::MeV, @@ -226,6 +229,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 3.6, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed @@ -242,6 +246,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 0.0203, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed diff --git a/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h b/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h index 7d4b275e38..5dc7a6a14d 100644 --- a/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h +++ b/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h @@ -29,6 +29,7 @@ class CalorimeterClusterRecoCoG_factory : public JOmniFactory m_logWeightBase {this, "logWeightBase", config().logWeightBase}; ParameterRef> m_logWeightBaseCoeffs {this, "logWeightBaseCoeffs", config().logWeightBaseCoeffs}; ParameterRef m_logWeightBase_Eref {this, "logWeightBase_Eref", config().logWeightBase_Eref}; + ParameterRef m_longitudinalShowerInfoAvailable {this, "longitudinalShowerInfoAvailable", config().longitudinalShowerInfoAvailable}; ParameterRef m_enableEtaBounds {this, "enableEtaBounds", config().enableEtaBounds}; Service m_algorithmsInit {this}; diff --git a/src/factories/calorimetry/ImagingTopoCluster_factory.h b/src/factories/calorimetry/ImagingTopoCluster_factory.h index 1cbbfef8c2..c46782203c 100644 --- a/src/factories/calorimetry/ImagingTopoCluster_factory.h +++ b/src/factories/calorimetry/ImagingTopoCluster_factory.h @@ -21,12 +21,14 @@ class ImagingTopoCluster_factory : public JOmniFactory> m_ldxy {this, "localDistXY", config().localDistXY}; ParameterRef> m_ldep {this, "layerDistEtaPhi", config().layerDistEtaPhi}; + ParameterRef> m_ldxy_adjacent {this, "layerDistXY", config().layerDistXY}; + ParameterRef m_laymode {this, "layerMode", config().layerMode}; ParameterRef m_nlr {this, "neighbourLayersRange", config().neighbourLayersRange}; ParameterRef m_sd {this, "sectorDist", config().sectorDist}; ParameterRef m_mched {this, "minClusterHitEdep", config().minClusterHitEdep}; ParameterRef m_mcced {this, "minClusterCenterEdep", config().minClusterCenterEdep}; ParameterRef m_mced {this, "minClusterEdep", config().minClusterEdep}; - ParameterRef m_mcnh {this, "minClusterNhits", config().minClusterNhits}; + ParameterRef m_mcnh {this, "minClusterNhits", config().minClusterNhits}; Service m_algorithmsInit {this}; diff --git a/src/factories/reco/FarForwardNeutronReconstruction_factory.h b/src/factories/reco/FarForwardNeutronReconstruction_factory.h index 127d1cd10d..680d2ffb84 100644 --- a/src/factories/reco/FarForwardNeutronReconstruction_factory.h +++ b/src/factories/reco/FarForwardNeutronReconstruction_factory.h @@ -16,9 +16,11 @@ namespace eicrecon { using AlgoT = eicrecon::FarForwardNeutronReconstruction; private: std::unique_ptr m_algo; - PodioInput m_clusters_input {this}; + PodioInput m_clusters_hcal_input {this}; + PodioInput m_clusters_ecal_input {this}; PodioOutput m_neutrons_output {this}; - ParameterRef> m_scale_corr_coeff {this, "scale_corr_coeff", config().scale_corr_coeff}; + ParameterRef> m_scale_corr_coeff_hcal {this, "scale_corr_coeff_hcal", config().scale_corr_coeff_hcal}; + ParameterRef> m_scale_corr_coeff_ecal {this, "scale_corr_coeff_ecal", config().scale_corr_coeff_ecal}; Service m_algorithmsInit {this}; public: @@ -34,7 +36,7 @@ namespace eicrecon { } void Process(int64_t run_number, uint64_t event_number) { - m_algo->process({m_clusters_input()},{m_neutrons_output().get()}); + m_algo->process({m_clusters_hcal_input(), m_clusters_ecal_input()},{m_neutrons_output().get()}); } }; diff --git a/src/factories/reco/HadronicFinalState_factory.h b/src/factories/reco/HadronicFinalState_factory.h index 00d03523e3..10f7bf6905 100644 --- a/src/factories/reco/HadronicFinalState_factory.h +++ b/src/factories/reco/HadronicFinalState_factory.h @@ -11,6 +11,7 @@ #include #include "extensions/jana/JOmniFactory.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" namespace eicrecon { @@ -29,6 +30,8 @@ class HadronicFinalState_factory : typename FactoryT::template PodioInput m_rc_particles_assoc_input {this}; typename FactoryT::template PodioOutput m_hadronic_final_state_output {this}; + typename FactoryT::template Service m_algorithmsInit {this}; + public: void Configure() { m_algo = std::make_unique(this->GetPrefix()); diff --git a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h index 7c395f39dc..dbc7e922e4 100644 --- a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h +++ b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h @@ -11,6 +11,7 @@ #include #include "extensions/jana/JOmniFactory.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" namespace eicrecon { @@ -29,6 +30,8 @@ class InclusiveKinematicsReconstructed_factory : typename FactoryT::template PodioInput m_hadronic_final_state_input {this}; typename FactoryT::template PodioOutput m_inclusive_kinematics_output {this}; + typename FactoryT::template Service m_algorithmsInit {this}; + public: void Configure() { m_algo = std::make_unique(this->GetPrefix()); diff --git a/src/factories/reco/InclusiveKinematicsTruth_factory.h b/src/factories/reco/InclusiveKinematicsTruth_factory.h index b119c0541b..52adcd7d0f 100644 --- a/src/factories/reco/InclusiveKinematicsTruth_factory.h +++ b/src/factories/reco/InclusiveKinematicsTruth_factory.h @@ -12,6 +12,7 @@ #include "algorithms/reco/InclusiveKinematicsTruth.h" #include "extensions/jana/JOmniFactory.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" namespace eicrecon { @@ -26,6 +27,8 @@ class InclusiveKinematicsTruth_factory : PodioInput m_mc_particles_input {this}; PodioOutput m_inclusive_kinematics_output {this}; + Service m_algorithmsInit {this}; + public: void Configure() { m_algo = std::make_unique(GetPrefix()); diff --git a/src/factories/reco/TransformBreitFrame_factory.h b/src/factories/reco/TransformBreitFrame_factory.h index 3a9852818b..ce1ee17d5d 100644 --- a/src/factories/reco/TransformBreitFrame_factory.h +++ b/src/factories/reco/TransformBreitFrame_factory.h @@ -14,6 +14,7 @@ #include "algorithms/reco/TransformBreitFrame.h" #include "extensions/jana/JOmniFactory.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" namespace eicrecon { @@ -33,6 +34,8 @@ namespace eicrecon { // output collection PodioOutput m_out_part {this}; + Service m_algorithmsInit {this}; + public: void Configure() { diff --git a/src/factories/reco/UndoAfterBurnerMCParticles_factory.h b/src/factories/reco/UndoAfterBurnerMCParticles_factory.h new file mode 100644 index 0000000000..544aea00eb --- /dev/null +++ b/src/factories/reco/UndoAfterBurnerMCParticles_factory.h @@ -0,0 +1,50 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Alex Jentsch, Jihee Kim, Brian Page +// + +#include "algorithms/reco/UndoAfterBurner.h" +#include "algorithms/reco/UndoAfterBurnerConfig.h" + +// Event Model related classes +#include +#include +#include + +#include "extensions/jana/JOmniFactory.h" + +namespace eicrecon { + +class UndoAfterBurnerMCParticles_factory : + public JOmniFactory { + +public: + using AlgoT = eicrecon::UndoAfterBurner; +private: + std::unique_ptr m_algo; + + PodioInput m_mcparts_input {this}; + PodioOutput m_postburn_output {this}; + + ParameterRef m_pid_assume_pion_mass {this, "m_pid_assume_pion_mass", config().m_pid_assume_pion_mass}; + ParameterRef m_crossing_angle {this, "m_crossing_angle", config().m_crossing_angle}; + ParameterRef m_pid_purity {this, "m_pid_purity", config().m_pid_purity}; + ParameterRef m_correct_beam_FX {this, "m_correct_beam_FX", config().m_correct_beam_FX}; + ParameterRef m_pid_use_MC_truth {this, "m_pid_use_MC_truth", config().m_pid_use_MC_truth}; + + +public: + void Configure() { + m_algo = std::make_unique(GetPrefix()); + m_algo->applyConfig(config()); + } + + void ChangeRun(int64_t run_number) { + } + + void Process(int64_t run_number, uint64_t event_number) { + m_algo->process({m_mcparts_input()}, {m_postburn_output().get()}); + } + +}; + +} // eicrecon diff --git a/src/global/beam/beam.cc b/src/global/beam/beam.cc index dc9595a462..62435fa0a6 100644 --- a/src/global/beam/beam.cc +++ b/src/global/beam/beam.cc @@ -5,14 +5,15 @@ #include #include -#include #include +#include #include #include #include "algorithms/interfaces/WithPodConfig.h" #include "algorithms/meta/SubDivideFunctors.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" +#include "factories/meta/CollectionCollector_factory.h" #include "factories/meta/SubDivideCollection_factory.h" extern "C" { @@ -22,8 +23,8 @@ extern "C" { using namespace eicrecon; // Divide MCParticle collection based on generator status and PDG - std::vector outCollections{"MCBeamElectrons","MCBeamProtons","MCScatteredElectrons","MCScatteredProtons"}; - std::vector> values{{4,11},{4,2212},{1,11},{1,2212}}; + std::vector outCollections{"MCBeamElectrons","MCBeamProtons","MCBeamNeutrons","MCScatteredElectrons","MCScatteredProtons","MCScatteredNeutrons"}; + std::vector> values{{4,11},{4,2212},{4,2112},{1,11},{1,2212},{1,2112}}; app->Add(new JOmniFactoryGeneratorT>( "BeamParticles", @@ -36,5 +37,14 @@ extern "C" { ) ); + // Combine beam protons and neutrons into beam hadrons + app->Add(new JOmniFactoryGeneratorT>( + "MCBeamHadrons", + {"MCBeamProtons","MCBeamNeutrons"}, + {"MCBeamHadrons"}, + app + ) + ); + } } diff --git a/src/global/pid/IrtCherenkovParticleID_factory.h b/src/global/pid/IrtCherenkovParticleID_factory.h index 0ae94fe589..15018ee7d3 100644 --- a/src/global/pid/IrtCherenkovParticleID_factory.h +++ b/src/global/pid/IrtCherenkovParticleID_factory.h @@ -15,6 +15,7 @@ #include "algorithms/pid/IrtCherenkovParticleID.h" #include "algorithms/pid/IrtCherenkovParticleIDConfig.h" #include "extensions/jana/JOmniFactory.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" #include "services/geometry/richgeo/RichGeo_service.h" namespace eicrecon { @@ -50,6 +51,7 @@ class IrtCherenkovParticleID_factory : ParameterRef m_cheatPhotonVertex {this, "cheatPhotonVertex", config().cheatPhotonVertex, ""}; ParameterRef m_cheatTrueRadiator {this, "cheatTrueRadiator", config().cheatTrueRadiator, ""}; + Service m_algorithmsInit {this}; Service m_RichGeoSvc {this}; public: diff --git a/src/global/pid_lut/PIDLookup_factory.h b/src/global/pid_lut/PIDLookup_factory.h index 964fa85ddc..6c56275941 100644 --- a/src/global/pid_lut/PIDLookup_factory.h +++ b/src/global/pid_lut/PIDLookup_factory.h @@ -10,6 +10,7 @@ #include "algorithms/pid_lut/PIDLookup.h" #include "extensions/jana/JOmniFactory.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" namespace eicrecon { @@ -30,6 +31,8 @@ class PIDLookup_factory : public JOmniFactory m_system{this, "system", config().system, "For the ParticleID record"}; + Service m_algorithmsInit {this}; + public: void Configure() { m_algo = std::make_unique(this->GetPrefix()); diff --git a/src/global/reco/CMakeLists.txt b/src/global/reco/CMakeLists.txt index 90886e8147..4a12da6e82 100644 --- a/src/global/reco/CMakeLists.txt +++ b/src/global/reco/CMakeLists.txt @@ -26,6 +26,8 @@ plugin_glob_all(${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 - algorithms_onnx_library algorithms_reco_library) +plugin_link_libraries(${PLUGIN_NAME} algorithms_digi_library + algorithms_tracking_library algorithms_reco_library) +if(USE_ONNX) + plugin_link_libraries(${PLUGIN_NAME} algorithms_onnx_library) +endif() diff --git a/src/global/reco/reco.cc b/src/global/reco/reco.cc index 1590a0091d..9813baac82 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -1,12 +1,12 @@ -// Copyright 2022, Dmitry Romanov -// Subject to the terms in the LICENSE file found in the top-level directory. -// -// +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 - 2024, Dmitry Romanov, Nathan Brei, Tooba Ali, Wouter Deconinck, Dmitry Kalinkin, John Lajoie, Simon Gardner, Tristan Protzman, Daniel Brandenburg, Derek M Anderson, Sebouh Paul, Tyler Kutz, Alex Jentsch, Jihee Kim, Brian Page +#include #include #include #include +#include #include #include #include @@ -19,16 +19,18 @@ #if EDM4EIC_VERSION_MAJOR >= 6 #include "algorithms/reco/HadronicFinalState.h" #include "algorithms/reco/InclusiveKinematicsDA.h" +#include "algorithms/reco/InclusiveKinematicsESigma.h" #include "algorithms/reco/InclusiveKinematicsElectron.h" #include "algorithms/reco/InclusiveKinematicsJB.h" #include "algorithms/reco/InclusiveKinematicsSigma.h" -#include "algorithms/reco/InclusiveKinematicseSigma.h" #endif #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/meta/CollectionCollector_factory.h" #include "factories/meta/FilterMatching_factory.h" #include "factories/reco/FarForwardNeutronReconstruction_factory.h" +#ifdef USE_ONNX #include "factories/reco/InclusiveKinematicsML_factory.h" +#endif #if EDM4EIC_VERSION_MAJOR >= 6 #include "factories/reco/InclusiveKinematicsReconstructed_factory.h" #endif @@ -38,6 +40,7 @@ #if EDM4EIC_VERSION_MAJOR >= 6 #include "factories/reco/HadronicFinalState_factory.h" #endif +#include "factories/reco/UndoAfterBurnerMCParticles_factory.h" #include "global/reco/ChargedReconstructedParticleSelector_factory.h" #include "global/reco/MC2SmearedParticle_factory.h" #include "global/reco/MatchClusters_factory.h" @@ -156,13 +159,25 @@ void InitPlugin(JApplication *app) { app )); - app->Add(new JOmniFactoryGeneratorT>( - "InclusiveKinematicseSigma", + app->Add(new JOmniFactoryGeneratorT>( + "InclusiveKinematicsESigma", { "MCParticles", "ScatteredElectronsTruth", "HadronicFinalState" }, + { + "InclusiveKinematicsESigma" + }, + app + )); + + // InclusiveKinematicseSigma is deprecated and will be removed, use InclusiveKinematicsESigma instead + app->Add(new JOmniFactoryGeneratorT>( + "InclusiveKinematicseSigma_legacy", + { + "InclusiveKinematicsESigma" + }, { "InclusiveKinematicseSigma" }, @@ -183,6 +198,7 @@ void InitPlugin(JApplication *app) { app )); +#ifdef USE_ONNX app->Add(new JOmniFactoryGeneratorT( "InclusiveKinematicsML", { @@ -194,6 +210,7 @@ void InitPlugin(JApplication *app) { }, app )); +#endif #endif app->Add(new JOmniFactoryGeneratorT( @@ -292,11 +309,11 @@ void InitPlugin(JApplication *app) { )); app->Add(new JOmniFactoryGeneratorT( "ReconstructedFarForwardZDCNeutrons", - {"HcalFarForwardZDCClusters"}, // edm4eic::ClusterCollection + {"HcalFarForwardZDCClusters","EcalFarForwardZDCClusters"}, // edm4eic::ClusterCollection {"ReconstructedFarForwardZDCNeutrons"}, // edm4eic::ReconstrutedParticleCollection, { - .scale_corr_coeff={-0.0756, -1.91, 2.30} - + .scale_corr_coeff_hcal={-0.0756, -1.91, 2.30}, + .scale_corr_coeff_ecal={-0.352, -1.34, 1.61} }, app // TODO: Remove me once fixed )); @@ -347,5 +364,26 @@ void InitPlugin(JApplication *app) { app )); + + //Full correction for MCParticles --> MCParticlesHeadOnFrame + app->Add(new JOmniFactoryGeneratorT( + "MCParticlesHeadOnFrameNoBeamFX", + { + "MCParticles" + }, + { + "MCParticlesHeadOnFrameNoBeamFX" + }, + { + .m_pid_assume_pion_mass = false, + .m_crossing_angle = -0.025 * dd4hep::rad, + .m_pid_purity = 0.51, //dummy value for MC truth information + .m_correct_beam_FX = true, + .m_pid_use_MC_truth = true, + }, + app + )); + + } } // extern "C" diff --git a/src/global/tracking/ActsToTracks_factory.h b/src/global/tracking/ActsToTracks_factory.h index 8fdc30a3ce..97bb4caf71 100644 --- a/src/global/tracking/ActsToTracks_factory.h +++ b/src/global/tracking/ActsToTracks_factory.h @@ -20,9 +20,11 @@ class ActsToTracks_factory PodioInput m_measurements_input {this}; Input m_acts_trajectories_input {this}; + PodioInput m_raw_hit_assocs_input {this}; PodioOutput m_trajectories_output {this}; PodioOutput m_parameters_output {this}; PodioOutput m_tracks_output {this}; + PodioOutput m_track_assocs_output {this}; public: void Configure() { @@ -38,8 +40,19 @@ class ActsToTracks_factory for (auto acts_traj : m_acts_trajectories_input()) { acts_trajectories_input.push_back(acts_traj); } - m_algo->process({m_measurements_input(), acts_trajectories_input}, - {m_trajectories_output().get(), m_parameters_output().get(), m_tracks_output().get()}); + m_algo->process( + { + m_measurements_input(), + acts_trajectories_input, + m_raw_hit_assocs_input(), + }, + { + m_trajectories_output().get(), + m_parameters_output().get(), + m_tracks_output().get(), + m_track_assocs_output().get(), + } + ); } }; diff --git a/src/global/tracking/TrackParamTruthInit_factory.h b/src/global/tracking/TrackParamTruthInit_factory.h index 4196b66c8b..09bd09dda3 100644 --- a/src/global/tracking/TrackParamTruthInit_factory.h +++ b/src/global/tracking/TrackParamTruthInit_factory.h @@ -15,6 +15,7 @@ #include "algorithms/tracking/TrackParamTruthInit.h" #include "algorithms/tracking/TrackParamTruthInitConfig.h" #include "extensions/jana/JOmniFactory.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" namespace eicrecon { @@ -37,6 +38,7 @@ class TrackParamTruthInit_factory : ParameterRef m_momentumSmear {this, "MomentumSmear", config().momentumSmear, "Momentum magnitude fraction to use as width of gaussian smearing"}; Service m_ACTSGeoSvc {this}; + Service m_algorithmsInit {this}; public: void Configure() { diff --git a/src/global/tracking/TracksToParticles_factory.h b/src/global/tracking/TracksToParticles_factory.h index 4f02d9a368..a0e82741cb 100644 --- a/src/global/tracking/TracksToParticles_factory.h +++ b/src/global/tracking/TracksToParticles_factory.h @@ -4,7 +4,6 @@ #pragma once #include "algorithms/tracking/TracksToParticles.h" -#include "algorithms/tracking/TracksToParticlesConfig.h" #include #include #include @@ -15,23 +14,18 @@ namespace eicrecon { class TracksToParticles_factory - : public JOmniFactory { + : public JOmniFactory { public: using AlgoT = eicrecon::TracksToParticles; private: std::unique_ptr m_algo; - PodioInput m_particles_input{this}; PodioInput m_tracks_input{this}; + PodioInput m_trackassocs_input{this}; PodioOutput m_recoparticles_output{this}; PodioOutput m_recoassocs_output{this}; - ParameterRef m_momentum_relative_tolerance{this, "momentumRelativeTolerance", - config().momentumRelativeTolerance}; - ParameterRef m_phi_tolerance{this, "phiTolerance", config().phiTolerance}; - ParameterRef m_eta_tolerance{this, "etaTolerance", config().etaTolerance}; - public: void Configure() { m_algo = std::make_unique(this->GetPrefix()); @@ -43,7 +37,7 @@ class TracksToParticles_factory void ChangeRun(int64_t run_number) {}; void Process(int64_t run_number, uint64_t event_number) { - m_algo->process({m_particles_input(), m_tracks_input()}, + m_algo->process({m_tracks_input(), m_trackassocs_input()}, {m_recoparticles_output().get(), m_recoassocs_output().get()}); } }; diff --git a/src/global/tracking/tracking.cc b/src/global/tracking/tracking.cc index 06a259441c..7a916effd8 100644 --- a/src/global/tracking/tracking.cc +++ b/src/global/tracking/tracking.cc @@ -1,15 +1,17 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 - 2024, Dmitry Romanov, Tyler Kutz, Wouter Deconinck +// Copyright (C) 2022 - 2024, Dmitry Romanov, Tyler Kutz, Wouter Deconinck, Dmitry Kalinkin #include #include -#include +#include +#include +#include #include #include #include #include #include -#include +#include #include #include "ActsToTracks.h" @@ -23,7 +25,6 @@ #include "TrackPropagation_factory.h" #include "TrackSeeding_factory.h" #include "TrackerMeasurementFromHits_factory.h" -#include "TracksToParticlesConfig.h" #include "TracksToParticles_factory.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/meta/CollectionCollector_factory.h" @@ -45,37 +46,45 @@ void InitPlugin(JApplication *app) { )); // Possible collections from arches, brycecanyon and craterlake configurations - std::vector> possible_collections = { - {"SiBarrelHits", "SiBarrelTrackerRecHits"}, - {"VertexBarrelHits", "SiBarrelVertexRecHits"}, - {"TrackerEndcapHits", "SiEndcapTrackerRecHits"}, - {"TOFBarrelHits", "TOFBarrelRecHit"}, - {"TOFEndcapHits", "TOFEndcapRecHits"}, - {"MPGDBarrelHits", "MPGDBarrelRecHits"}, - {"MPDGDIRCHits", "MPDGDIRCRecHits"}, - {"OuterMPGDBarrelHits", "OuterMPGDBarrelRecHits"}, - {"BackwardMPGDEndcapHits", "BackwardMPGDEndcapRecHits"}, - {"ForwardMPGDEndcapHits", "ForwardMPGDEndcapRecHits"}, - {"B0TrackerHits", "B0TrackerRecHits"} + std::vector> possible_collections = { + {"SiBarrelHits", "SiBarrelRawHits", "SiBarrelRawHitAssociations", "SiBarrelTrackerRecHits"}, + {"VertexBarrelHits", "SiBarrelVertexRawHits", "SiBarrelVertexRawHitAssociations", "SiBarrelVertexRecHits"}, + {"TrackerEndcapHits", "SiEndcapTrackerRawHits", "SiEndcapTrackerRawHitAssociations", "SiEndcapTrackerRecHits"}, + {"TOFBarrelHits", "TOFBarrelRawHits", "TOFBarrelRawHitAssociations", "TOFBarrelRecHit"}, + {"TOFEndcapHits", "TOFEndcapRawHits", "TOFEndcapRawHitAssociations", "TOFEndcapRecHits"}, + {"MPGDBarrelHits", "MPGDBarrelRawHits", "MPGDBarrelRawHitAssociations", "MPGDBarrelRecHits"}, + {"OuterMPGDBarrelHits", "OuterMPGDBarrelRawHits", "OuterMPGDBarrelRawHitAssociations", "OuterMPGDBarrelRecHits"}, + {"BackwardMPGDEndcapHits", "BackwardMPGDEndcapRawHits", "BackwardMPGDEndcapRawHitAssociations", "BackwardMPGDEndcapRecHits"}, + {"ForwardMPGDEndcapHits", "ForwardMPGDEndcapRawHits", "ForwardMPGDEndcapRawHitAssociations", "ForwardMPGDEndcapRecHits"}, + {"B0TrackerHits", "B0TrackerRawHits", "B0TrackerRawHitAssociations", "B0TrackerRecHits"} }; // Filter out collections that are not present in the current configuration - std::vector input_collections; + std::vector input_rec_collections; + std::vector input_raw_assoc_collections; auto readouts = app->GetService()->detector()->readouts(); - for (const auto& [hit_collection, rec_collection] : possible_collections) { + for (const auto& [hit_collection, raw_collection, raw_assoc_collection, rec_collection] : possible_collections) { if (readouts.find(hit_collection) != readouts.end()) { // Add the collection to the list of input collections - input_collections.push_back(rec_collection); + input_rec_collections.push_back(rec_collection); + input_raw_assoc_collections.push_back(raw_assoc_collection); } } // Tracker hits collector app->Add(new JOmniFactoryGeneratorT>( "CentralTrackingRecHits", - input_collections, + input_rec_collections, {"CentralTrackingRecHits"}, // Output collection name app)); + // Tracker hit associations collector + app->Add(new JOmniFactoryGeneratorT>( + "CentralTrackingRawHitAssociations", + input_raw_assoc_collections, + {"CentralTrackingRawHitAssociations"}, // Output collection name + app)); + app->Add(new JOmniFactoryGeneratorT( "CentralTrackerMeasurements", {"CentralTrackingRecHits"}, @@ -101,11 +110,13 @@ void InitPlugin(JApplication *app) { { "CentralTrackerMeasurements", "CentralCKFActsTrajectoriesUnfiltered", + "CentralTrackingRawHitAssociations", }, { "CentralCKFTrajectoriesUnfiltered", "CentralCKFTrackParametersUnfiltered", "CentralCKFTracksUnfiltered", + "CentralCKFTrackUnfilteredAssociations", }, app )); @@ -128,11 +139,13 @@ void InitPlugin(JApplication *app) { { "CentralTrackerMeasurements", "CentralCKFActsTrajectories", + "CentralTrackingRawHitAssociations", }, { "CentralCKFTrajectories", "CentralCKFTrackParameters", "CentralCKFTracks", + "CentralCKFTrackAssociations", }, app )); @@ -163,11 +176,13 @@ void InitPlugin(JApplication *app) { { "CentralTrackerMeasurements", "CentralCKFSeededActsTrajectoriesUnfiltered", + "CentralTrackingRawHitAssociations", }, { "CentralCKFSeededTrajectoriesUnfiltered", "CentralCKFSeededTrackParametersUnfiltered", "CentralCKFSeededTracksUnfiltered", + "CentralCKFSeededTrackUnfilteredAssociations", }, app )); @@ -190,11 +205,13 @@ void InitPlugin(JApplication *app) { { "CentralTrackerMeasurements", "CentralCKFSeededActsTrajectories", + "CentralTrackingRawHitAssociations", }, { "CentralCKFSeededTrajectories", "CentralCKFSeededTrackParameters", "CentralCKFSeededTracks", + "CentralCKFSeededTrackAssociations", }, app )); @@ -251,35 +268,51 @@ void InitPlugin(JApplication *app) { app )); - // linking of reconstructed particles to PID objects - TracksToParticlesConfig link_cfg { - .momentumRelativeTolerance = 100.0, /// Matching momentum effectively disabled - .phiTolerance = 0.1, /// Matching phi tolerance [rad] - .etaTolerance = 0.2, /// Matching eta tolerance - }; - app->Add(new JOmniFactoryGeneratorT( - "ChargedParticlesWithAssociations", - {"MCParticles", // edm4hep::MCParticle - "CentralCKFTracks", // edm4eic::Track - }, - {"ReconstructedChargedWithoutPIDParticles", // - "ReconstructedChargedWithoutPIDParticleAssociations" // edm4eic::MCRecoParticleAssociation - }, - link_cfg, - app - )); - app->Add(new JOmniFactoryGeneratorT( - "ChargedSeededParticlesWithAssociations", - {"MCParticles", // edm4hep::MCParticle - "CentralCKFSeededTracks", // edm4eic::Track - }, - {"ReconstructedSeededChargedWithoutPIDParticles", // - "ReconstructedSeededChargedWithoutPIDParticleAssociations" // edm4eic::MCRecoParticleAssociation - }, - link_cfg, - app - )); + std::vector input_track_collections; + //Check size of input_rec_collections to determine if CentralCKFTracks should be added to the input_track_collections + if (input_rec_collections.size() > 0) { + input_track_collections.push_back("CentralCKFTracks"); + } + //Check if the TaggerTracker readout is present in the current configuration + if (readouts.find("TaggerTrackerHits") != readouts.end()) { + input_track_collections.push_back("TaggerTrackerTracks"); + } + + // Add central and other tracks + app->Add(new JOmniFactoryGeneratorT>( + "CombinedTracks", + input_track_collections, + {"CombinedTracks"}, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "ChargedParticlesWithAssociations", + { + "CombinedTracks", + "CentralCKFTrackAssociations", + }, + {"ReconstructedChargedWithoutPIDParticles", + "ReconstructedChargedWithoutPIDParticleAssociations" + }, + {}, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "ChargedSeededParticlesWithAssociations", + { + "CentralCKFSeededTracks", + "CentralCKFSeededTrackAssociations", + }, + { + "ReconstructedSeededChargedWithoutPIDParticles", + "ReconstructedSeededChargedWithoutPIDParticleAssociations" + }, + {}, + app + )); } } // extern "C" diff --git a/src/services/algorithms_init/AlgorithmsInit_service.h b/src/services/algorithms_init/AlgorithmsInit_service.h index 567799ee3f..10b68adcff 100644 --- a/src/services/algorithms_init/AlgorithmsInit_service.h +++ b/src/services/algorithms_init/AlgorithmsInit_service.h @@ -13,6 +13,7 @@ #include #include +#include "algorithms/interfaces/ParticleSvc.h" #include "services/log/Log_service.h" #include "services/geometry/dd4hep/DD4hep_service.h" @@ -50,12 +51,13 @@ class AlgorithmsInit_service : public JService logger.init([this](const algorithms::LogLevel l, std::string_view caller, std::string_view msg){ static std::mutex m; std::lock_guard lock(m); - static std::map> loggers; - if (! loggers.contains(caller)) { + // storing the string_view is unsafe since it can become invalid + static std::map> loggers; + if (! loggers.contains(std::string(caller))) { this->m_log->debug("Initializing algorithms::LogSvc logger {}", caller); - loggers[caller] = this->m_log_service->logger(std::string(caller)); + loggers[std::string(caller)] = this->m_log_service->logger(std::string(caller)); } - loggers[caller]->log(static_cast(l), msg); + loggers[std::string(caller)]->log(static_cast(l), msg); }); logger.defaultLevel(level); }); @@ -68,6 +70,9 @@ class AlgorithmsInit_service : public JService r.init(); }); + // Register a particle service + [[maybe_unused]] auto& particleSvc = algorithms::ParticleSvc::instance(); + // Finally, initialize the ServiceSvc serviceSvc.init(); } diff --git a/src/services/algorithms_init/CMakeLists.txt b/src/services/algorithms_init/CMakeLists.txt index 4b3665f433..13fca135c2 100644 --- a/src/services/algorithms_init/CMakeLists.txt +++ b/src/services/algorithms_init/CMakeLists.txt @@ -18,3 +18,4 @@ plugin_glob_all(${PLUGIN_NAME}) plugin_add_algorithms(${PLUGIN_NAME}) plugin_add_dd4hep(${PLUGIN_NAME}) plugin_add_event_model(${PLUGIN_NAME}) +plugin_link_libraries(${PLUGIN_NAME} algorithms_interfaces_library) diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 7b36d2f66f..73e3be07b5 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -9,8 +9,12 @@ #include #include #include +#include +#if podio_VERSION >= PODIO_VERSION(0, 99, 0) #include -#include +#else +#include +#endif #include #include #include @@ -55,9 +59,11 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "MCBeamProtons", "MCScatteredElectrons", "MCScatteredProtons", + "MCParticlesHeadOnFrameNoBeamFX", // All tracking hits combined "CentralTrackingRecHits", + "CentralTrackingRawHitAssociations", "CentralTrackSeedingResults", "CentralTrackerMeasurements", @@ -74,9 +80,9 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "VertexBarrelHits", "TrackerEndcapHits", - "SiBarrelHitAssociations", - "SiBarrelVertexHitAssociations", - "SiEndcapHitAssociations", + "SiBarrelRawHitAssociations", + "SiBarrelVertexRawHitAssociations", + "SiEndcapTrackerRawHitAssociations", // TOF "TOFBarrelRecHit", @@ -88,8 +94,8 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "TOFBarrelHits", "TOFEndcapHits", - "TOFBarrelHitAssociations", - "TOFEndcapHitAssociations", + "TOFBarrelRawHitAssociations", + "TOFEndcapRawHitAssociations", "CombinedTOFParticleIDs", "CombinedTOFSeededParticleIDs", @@ -126,14 +132,14 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "BackwardMPGDEndcapHits", "ForwardMPGDEndcapHits", - "MPGDBarrelHitAssociations", - "OuterMPGDBarrelHitAssociations", - "BackwardMPGDEndcapAssociations", - "ForwardMPGDHitAssociations", + "MPGDBarrelRawHitAssociations", + "OuterMPGDBarrelRawHitAssociations", + "BackwardMPGDEndcapRawHitAssociations", + "ForwardMPGDEndcapRawHitAssociations", // LOWQ2 hits "TaggerTrackerRawHits", - "TaggerTrackerHitAssociations", + "TaggerTrackerRawHitAssociations", "TaggerTrackerM1L0ClusterPositions", "TaggerTrackerM1L1ClusterPositions", "TaggerTrackerM1L2ClusterPositions", @@ -153,7 +159,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "B0TrackerRecHits", "B0TrackerRawHits", "B0TrackerHits", - "B0TrackerHitAssociations", + "B0TrackerRawHitAssociations", "ForwardRomanPotRecHits", "ForwardOffMTrackerRecHits", @@ -161,8 +167,8 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "ForwardRomanPotRecParticles", "ForwardOffMRecParticles", - "ForwardRomanPotHitAssociations", - "ForwardOffMTrackerHitAssociations", + "ForwardRomanPotRawHitAssociations", + "ForwardOffMTrackerRawHitAssociations", // Reconstructed data "GeneratedParticles", @@ -182,22 +188,28 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "CentralTrackVertices", "CentralCKFTrajectories", "CentralCKFTracks", + "CentralCKFTrackAssociations", "CentralCKFTrackParameters", "CentralCKFSeededTrajectories", "CentralCKFSeededTracks", + "CentralCKFSeededTrackAssociations", "CentralCKFSeededTrackParameters", //tracking properties - true seeding "CentralCKFTrajectoriesUnfiltered", "CentralCKFTracksUnfiltered", + "CentralCKFTrackUnfilteredAssociations", "CentralCKFTrackParametersUnfiltered", //tracking properties - realistic seeding "CentralCKFSeededTrajectoriesUnfiltered", "CentralCKFSeededTracksUnfiltered", + "CentralCKFSeededTrackUnfilteredAssociations", "CentralCKFSeededTrackParametersUnfiltered", "InclusiveKinematicsDA", "InclusiveKinematicsJB", + "InclusiveKinematicsML", "InclusiveKinematicsSigma", - "InclusiveKinematicseSigma", + "InclusiveKinematicseSigma", // Deprecated, use ESigma + "InclusiveKinematicsESigma", "InclusiveKinematicsElectron", "InclusiveKinematicsTruth", "GeneratedJets", @@ -349,8 +361,11 @@ void JEventProcessorPODIO::Init() { auto *app = GetApplication(); m_log = app->GetService()->logger("JEventProcessorPODIO"); - m_log->set_level(spdlog::level::debug); +#if podio_VERSION >= PODIO_VERSION(0, 99, 0) + m_writer = std::make_unique(m_output_file); +#else m_writer = std::make_unique(m_output_file); +#endif // TODO: NWB: Verify that output file is writable NOW, rather than after event processing completes. // I definitely don't trust PODIO to do this for me. diff --git a/src/services/io/podio/JEventProcessorPODIO.h b/src/services/io/podio/JEventProcessorPODIO.h index 77d72679de..5b59f291dc 100644 --- a/src/services/io/podio/JEventProcessorPODIO.h +++ b/src/services/io/podio/JEventProcessorPODIO.h @@ -3,7 +3,12 @@ #include #include +#include +#if podio_VERSION >= PODIO_VERSION(0, 99, 0) +#include +#else #include +#endif #include #include #include @@ -25,7 +30,11 @@ class JEventProcessorPODIO : public JEventProcessor { void FindCollectionsToWrite(const std::shared_ptr& event); +#if podio_VERSION >= PODIO_VERSION(0, 99, 0) + std::unique_ptr m_writer; +#else std::unique_ptr m_writer; +#endif std::mutex m_mutex; bool m_is_first_event = true; bool m_user_included_collections = false; diff --git a/src/services/io/podio/README.md b/src/services/io/podio/README.md index 01dabc7dc8..dd0b4bb7bb 100644 --- a/src/services/io/podio/README.md +++ b/src/services/io/podio/README.md @@ -106,7 +106,7 @@ You may specify both an include list and an exclude list. Similar to the input, you may also specify which collections to write out using the -_podio:output_include_collections_ and _podio:output_exclude_collections_ configuration +_podio:output_collections_ and _podio:output_exclude_collections_ configuration parameters. ### Testing diff --git a/src/tests/algorithms_test/CMakeLists.txt b/src/tests/algorithms_test/CMakeLists.txt index 50ee83cdd0..57ecd29c88 100644 --- a/src/tests/algorithms_test/CMakeLists.txt +++ b/src/tests/algorithms_test/CMakeLists.txt @@ -6,6 +6,7 @@ add_executable( ${TEST_NAME} algorithmsInit.cc calorimetry_CalorimeterIslandCluster.cc + calorimetry_ImagingTopoCluster.cc tracking_SiliconSimpleCluster.cc calorimetry_CalorimeterHitDigi.cc calorimetry_CalorimeterClusterRecoCoG.cc @@ -22,6 +23,7 @@ target_link_libraries( PRIVATE Catch2::Catch2WithMain algorithms_calorimetry_library algorithms_fardetectors_library + algorithms_interfaces_library # for ParticleSvc algorithms_pid_library algorithms_pid_lut_library algorithms_reco_library diff --git a/src/tests/algorithms_test/algorithmsInit.cc b/src/tests/algorithms_test/algorithmsInit.cc index cdbb8b0779..357e928001 100644 --- a/src/tests/algorithms_test/algorithmsInit.cc +++ b/src/tests/algorithms_test/algorithmsInit.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -64,6 +65,9 @@ class algorithmsInitListener : public Catch::EventListenerBase { auto& lutSvc = eicrecon::PIDLookupTableSvc::instance(); serviceSvc.add(&lutSvc); + auto& particleSvc = algorithms::ParticleSvc::instance(); + serviceSvc.add(&particleSvc); + serviceSvc.init(); } }; diff --git a/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc index 32f11942d3..86358ee4b5 100644 --- a/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc +++ b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc @@ -4,6 +4,8 @@ #include // for MeV, mm, keV, ns #include // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE +#include +#include #include // for CalorimeterHitCollection, MutableCalorimeterHit, CalorimeterHitMutableCollectionIterator #include #include @@ -11,7 +13,6 @@ #include #include // for Vector3f #include -#include #include // for level_enum #include // for logger #include // for default_logger @@ -28,6 +29,8 @@ using eicrecon::CalorimeterClusterRecoCoGConfig; using edm4eic::CalorimeterHit; TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) { + const float EPSILON = 1e-5; + CalorimeterClusterRecoCoG algo("CalorimeterClusterRecoCoG"); std::shared_ptr logger = spdlog::default_logger()->clone("CalorimeterClusterRecoCoG"); @@ -35,9 +38,10 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) CalorimeterClusterRecoCoGConfig cfg; cfg.energyWeight = "log"; - cfg.sampFrac = 0.0203, - cfg.logWeightBaseCoeffs={5.0,0.65,0.31}, - cfg.logWeightBase_Eref=50*dd4hep::GeV, + cfg.sampFrac = 0.0203; + cfg.logWeightBaseCoeffs = {5.0,0.65,0.31}; + cfg.logWeightBase_Eref = 50*dd4hep::GeV; + cfg.longitudinalShowerInfoAvailable = true; algo.applyConfig(cfg); @@ -51,7 +55,7 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) //create a protocluster with 3 hits auto pclust = pclust_coll.create(); - edm4hep::Vector3f position({0,0,0*dd4hep::mm}); + edm4hep::Vector3f position({0, 0, 1 *dd4hep::mm}); auto hit1 = hits_coll.create(); hit1.setCellID(0); @@ -63,8 +67,9 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) hit1.setDimension({0,0,0}); hit1.setLocal(position); pclust.addToHits(hit1); + pclust.addToWeights(1); - position={0,0, 1*dd4hep::mm}; + position={-1 * dd4hep::mm, 0, 2 * dd4hep::mm}; auto hit2 = hits_coll.create(); hit2.setCellID(0); hit2.setEnergy(0.1*dd4hep::GeV); @@ -75,19 +80,7 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) hit2.setDimension({0,0,0}); hit2.setLocal(position); pclust.addToHits(hit2); - - position={0,0, 2*dd4hep::mm}; - auto hit3 = hits_coll.create();//0, 0.1*dd4hep::GeV, 0,0,0,position, {0,0,0}, 0,0, position); - hit3.setCellID(0); - hit3.setEnergy(0.1*dd4hep::GeV); - hit3.setEnergyError(0); - hit3.setTime(0); - hit3.setTimeError(0); - hit3.setPosition(position); - hit3.setDimension({0,0,0}); - hit3.setLocal(position); - pclust.addToHits(hit3); - pclust.addToWeights(1);pclust.addToWeights(1);pclust.addToWeights(1); + pclust.addToWeights(1); // Constructing input and output as per the algorithm's expected signature auto input = std::make_tuple(&pclust_coll, &simhits); @@ -97,10 +90,9 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) for (auto clust : *clust_coll){ - // require that this cluster's axis is 0,0,1 - REQUIRE(clust.getShapeParameters()[7] == 0); - REQUIRE(clust.getShapeParameters()[8] == 0); - REQUIRE(abs(clust.getShapeParameters()[9]) == 1); + REQUIRE_THAT(clust.getIntrinsicTheta(), Catch::Matchers::WithinAbs(M_PI / 4, EPSILON)); + // std::abs() checks if we land on -M_PI + REQUIRE_THAT(std::abs(clust.getIntrinsicPhi()), Catch::Matchers::WithinAbs(M_PI, EPSILON)); } diff --git a/src/tests/algorithms_test/calorimetry_CalorimeterIslandCluster.cc b/src/tests/algorithms_test/calorimetry_CalorimeterIslandCluster.cc index cd4d59eafe..5b09d41259 100644 --- a/src/tests/algorithms_test/calorimetry_CalorimeterIslandCluster.cc +++ b/src/tests/algorithms_test/calorimetry_CalorimeterIslandCluster.cc @@ -151,10 +151,16 @@ TEST_CASE( "the clustering algorithm runs", "[CalorimeterIslandCluster]" ) { bool use_adjacencyMatrix = GENERATE(false, true); if (use_adjacencyMatrix) { cfg.adjacencyMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1"; - cfg.readout = "MockCalorimeterHits"; } else { cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm}; } + bool disalow_diagonal_peaks = GENERATE(false, true); + if (disalow_diagonal_peaks) { + cfg.peakNeighbourhoodMatrix = "max(abs(x_1 - x_2), abs(y_1 - y_2)) == 1"; + } else { + cfg.peakNeighbourhoodMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1"; + } + cfg.readout = "MockCalorimeterHits"; cfg.splitCluster = GENERATE(false, true); if (cfg.splitCluster) { @@ -190,8 +196,17 @@ TEST_CASE( "the clustering algorithm runs", "[CalorimeterIslandCluster]" ) { 0, // std::int32_t layer, edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0) // edm4hep::Vector3f local ); + + bool test_diagonal_cluster = GENERATE(false, true); + // If false, test a cluster with two maxima: + // XxX + // If true, test a diagonal cluster: + // Xx + // X + // The idea is to test whether peakNeighbourhoodMatrix allows increasing + // peak resolution threshold while keeping the Island algorithm the same. hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 2}, {"y", 0}}), // std::uint64_t cellID, + id_desc.encode({{"system", 255}, {"x", test_diagonal_cluster ? 1 : 2}, {"y", test_diagonal_cluster ? 1 : 0}}), // std::uint64_t cellID, 6.0, // float energy, 0.0, // float energyError, 0.0, // float time, @@ -205,7 +220,11 @@ TEST_CASE( "the clustering algorithm runs", "[CalorimeterIslandCluster]" ) { auto protoclust_coll = std::make_unique(); algo.process({&hits_coll}, {protoclust_coll.get()}); - if (cfg.splitCluster) { + bool expect_two_peaks = cfg.splitCluster; + if (cfg.splitCluster && disalow_diagonal_peaks) { + expect_two_peaks = not test_diagonal_cluster; + } + if (expect_two_peaks) { REQUIRE( (*protoclust_coll).size() == 2 ); REQUIRE( (*protoclust_coll)[0].hits_size() == 3 ); REQUIRE( (*protoclust_coll)[0].weights_size() == 3 ); diff --git a/src/tests/algorithms_test/calorimetry_ImagingTopoCluster.cc b/src/tests/algorithms_test/calorimetry_ImagingTopoCluster.cc new file mode 100644 index 0000000000..4db2abc783 --- /dev/null +++ b/src/tests/algorithms_test/calorimetry_ImagingTopoCluster.cc @@ -0,0 +1,194 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Sebouh Paul + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "algorithms/calorimetry/ImagingTopoCluster.h" +#include "algorithms/calorimetry/ImagingTopoClusterConfig.h" + +using eicrecon::ImagingTopoCluster; +using eicrecon::ImagingTopoClusterConfig; + +TEST_CASE( "the clustering algorithm runs", "[ImagingTopoCluster]" ) { + ImagingTopoCluster algo("ImagingTopoCluster"); + + std::shared_ptr logger = spdlog::default_logger()->clone("ImagingTopoCluster"); + logger->set_level(spdlog::level::trace); + + ImagingTopoClusterConfig cfg; + cfg.layerMode =eicrecon::ImagingTopoClusterConfig::ELayerMode::xy; + cfg.minClusterHitEdep = 0. * dd4hep::GeV; + cfg.minClusterCenterEdep = 0. * dd4hep::GeV; + cfg.localDistXY = {1.0* dd4hep::mm, 1.0* dd4hep::mm}; //mm + cfg.layerDistXY = {1.0* dd4hep::mm, 1.0* dd4hep::mm}; //mm + cfg.minClusterEdep = 9 * dd4hep::MeV; + // minimum number of hits (to save this cluster) + cfg.minClusterNhits = 1; + auto detector = algorithms::GeoSvc::instance().detector(); + auto id_desc = detector->readout("MockCalorimeterHits").idSpec(); + + SECTION( "without splitting" ) { + algo.applyConfig(cfg); + algo.init(); + + SECTION( "on a single cell" ) { + edm4eic::CalorimeterHitCollection hits_coll; + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}, {"layer", 0}}), // std::uint64_t cellID, + 5.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 0, // std::int32_t layer, + edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local + ); + auto protoclust_coll = std::make_unique(); + algo.process({&hits_coll}, {protoclust_coll.get()}); + + REQUIRE( (*protoclust_coll).size() == 1 ); + REQUIRE( (*protoclust_coll)[0].hits_size() == 1 ); + REQUIRE( (*protoclust_coll)[0].weights_size() == 1 ); + } + + SECTION( "on two separated cells" ) { + edm4eic::CalorimeterHitCollection hits_coll; + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}, {"layer", 0}}), // std::uint64_t cellID, + 5.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 0, // std::int32_t layer, + edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local + ); + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 2}, {"y", 2}, {"layer", 0}}), // std::uint64_t cellID, + 6.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(1.1, 1.1, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 0, // std::int32_t layer, + edm4hep::Vector3f(1.1 /* mm */, 1.1 /* mm */, 0.0) // edm4hep::Vector3f local + ); + auto protoclust_coll = std::make_unique(); + algo.process({&hits_coll}, {protoclust_coll.get()}); + + REQUIRE( (*protoclust_coll).size() == 2 ); + REQUIRE( (*protoclust_coll)[0].hits_size() == 1 ); + REQUIRE( (*protoclust_coll)[0].weights_size() == 1 ); + REQUIRE( (*protoclust_coll)[1].hits_size() == 1 ); + REQUIRE( (*protoclust_coll)[1].weights_size() == 1 ); + } + + SECTION( "on two adjacent cells (same layer)" ) { + edm4eic::CalorimeterHitCollection hits_coll; + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 0, // std::int32_t layer, + edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local + ); + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, + 6.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(0.9, 0.9, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 0, // std::int32_t layer, + edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0) // edm4hep::Vector3f local + ); + auto protoclust_coll = std::make_unique(); + algo.process({&hits_coll}, {protoclust_coll.get()}); + + REQUIRE( (*protoclust_coll).size() == 1 ); + REQUIRE( (*protoclust_coll)[0].hits_size() == 2 ); + REQUIRE( (*protoclust_coll)[0].weights_size() == 2 ); + } + } + + SECTION( "run on three cells, two of which are on the same layer, and there is a third one on another layer acting as a bridge between them" ) { + + cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm}; + algo.applyConfig(cfg); + algo.init(); + + edm4eic::CalorimeterHitCollection hits_coll; + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}, {"layer", 0}}), // std::uint64_t cellID, + 5.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 0, // std::int32_t layer, + edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local + ); + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}, {"layer",1}}), // std::uint64_t cellID, + 1.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(0.9, 0.9, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 1, // std::int32_t layer, + edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0) // edm4hep::Vector3f local + ); + hits_coll.create( + id_desc.encode({{"system", 255}, {"x", 2}, {"y", 0},{"layer",0}}), // std::uint64_t cellID, + 6.0, // float energy, + 0.0, // float energyError, + 0.0, // float time, + 0.0, // float timeError, + edm4hep::Vector3f(1.8, 1.8, 0.0), // edm4hep::Vector3f position, + edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension, + 0, // std::int32_t sector, + 0, // std::int32_t layer, + edm4hep::Vector3f(1.8 /* mm */, 1.8 /* mm */, 0.0) // edm4hep::Vector3f local + ); + auto protoclust_coll = std::make_unique(); + algo.process({&hits_coll}, {protoclust_coll.get()}); + + + REQUIRE( (*protoclust_coll).size() == 1 ); + REQUIRE( (*protoclust_coll)[0].hits_size() == 3 ); + REQUIRE( (*protoclust_coll)[0].weights_size() == 3 ); + + } +} diff --git a/src/tests/algorithms_test/reco_FarForwardNeutronReconstruction.cc b/src/tests/algorithms_test/reco_FarForwardNeutronReconstruction.cc index 839bee27ef..24e4b3a3f4 100644 --- a/src/tests/algorithms_test/reco_FarForwardNeutronReconstruction.cc +++ b/src/tests/algorithms_test/reco_FarForwardNeutronReconstruction.cc @@ -15,10 +15,13 @@ #include // for sqrt, abs #include #include // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access +#include -#include "algorithms/reco/FarForwardNeutronReconstruction.h" // for Neutronreconstruction +#include "algorithms/reco/FarForwardNeutronReconstruction.h" +#include "algorithms/reco/FarForwardNeutronReconstructionConfig.h" using eicrecon::FarForwardNeutronReconstruction; +using eicrecon::FarForwardNeutronReconstructionConfig; TEST_CASE( "the cluster merging algorithm runs", "[FarForwardNeutronReconstruction]" ) { FarForwardNeutronReconstruction algo("FarForwardNeutronReconstruction"); @@ -26,32 +29,42 @@ TEST_CASE( "the cluster merging algorithm runs", "[FarForwardNeutronReconstructi std::shared_ptr logger = spdlog::default_logger()->clone("FarForwardNeutronReconstruction"); logger->set_level(spdlog::level::trace); + FarForwardNeutronReconstructionConfig cfg; + std::vector corr_parameters={-0.0756, -1.91, 2.30}; + cfg.scale_corr_coeff_hcal=corr_parameters; + cfg.scale_corr_coeff_ecal=corr_parameters; + algo.applyConfig(cfg); algo.init(); - edm4eic::ClusterCollection clust_coll; - + edm4eic::ClusterCollection clust_coll_hcal; std::array x={30*dd4hep::mm,90*dd4hep::mm,0}; std::array y={-30*dd4hep::mm,0*dd4hep::mm, -90*dd4hep::mm}; std::array z={30*dd4hep::m,30*dd4hep::m, 30*dd4hep::m}; std::array E={80*dd4hep::GeV,5*dd4hep::GeV,5*dd4hep::GeV}; float sumEnergies=0; for(size_t i=0; i<3; i++){ - auto cluster=clust_coll.create(); + auto cluster=clust_coll_hcal.create(); cluster.setEnergy(E[i]); cluster.setPosition({x[i], y[i], z[i]}); } + + edm4eic::ClusterCollection clust_coll_ecal; + auto ecal_cluster=clust_coll_ecal.create(); + ecal_cluster.setEnergy(2); + ecal_cluster.setPosition({0, 0, 25*dd4hep::m}); + auto neutroncand_coll = std::make_unique(); - algo.process({&clust_coll}, {neutroncand_coll.get()}); + algo.process({&clust_coll_hcal, &clust_coll_ecal}, {neutroncand_coll.get()}); REQUIRE( (*neutroncand_coll).size() == 1); - double corr=algo.calc_corr(90); + double corr=algo.calc_corr(92, corr_parameters); double tol=0.001; - double E_expected=90*dd4hep::GeV*1/(1+corr); - double Px_expected=0.08999*1/(1+corr); - double Py_expected=-0.08999*1/(1+corr); - double Pz_expected=89.99*1/(1+corr); + double E_expected=92*dd4hep::GeV*1/(1+corr); + double Px_expected=0.09199*1/(1+corr); + double Py_expected=-0.09199*1/(1+corr); + double Pz_expected=91.99*1/(1+corr); //check that the correct energy and momenta are being obtained std::cout << "E, px, py, pz = " << (*neutroncand_coll)[0].getEnergy() <<" " << (*neutroncand_coll)[0].getMomentum().x << " " << (*neutroncand_coll)[0].getMomentum().y << " " << (*neutroncand_coll)[0].getMomentum().z << std::endl; REQUIRE( abs((*neutroncand_coll)[0].getEnergy()-E_expected)/E_expectedRun(); } catch (JException& e) { diff --git a/src/utilities/eicrecon/print_info.h b/src/utilities/eicrecon/print_info.h index 4060d4b9fb..67676988d2 100644 --- a/src/utilities/eicrecon/print_info.h +++ b/src/utilities/eicrecon/print_info.h @@ -33,17 +33,3 @@ void printPluginNames(std::vector const& plugin_names) { plugin_table.Render(ss); std::cout << ss.str() << std::endl; } - -void printJANAHeaderIMG() { - std::cout << " ____ _ ___ ___ _ \n" - " `MM' dM. `MM\\ `M' dM. \n" - " MM ,MMb MMM\\ M ,MMb \n" - " MM d'YM. M\\MM\\ M d'YM. ____ \n" - " MM ,P `Mb M \\MM\\ M ,P `Mb 6MMMMb \n" - " MM d' YM. M \\MM\\ M d' YM. MM' `Mb \n" - " MM ,P `Mb M \\MM\\ M ,P `Mb ,MM \n" - " MM d' YM. M \\MM\\M d' YM. ,MM' \n" - "(8) MM ,MMMMMMMMb M \\MMM ,MMMMMMMMb ,M' \n" - "(( ,M9 d' YM. M \\MM d' YM. ,M' \n" - " YMMMM9 _dM_ _dMM_M_ \\M _dM_ _dMM_MMMMMMMM " << std::endl << std::endl; -}