diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 15971c1..323237b 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1,2 +1 @@ style = "blue" - diff --git a/.github/workflows/BuildDeployDoc.yml b/.github/workflows/BuildDeployDoc.yml index 171d9f0..0b9b040 100644 --- a/.github/workflows/BuildDeployDoc.yml +++ b/.github/workflows/BuildDeployDoc.yml @@ -3,7 +3,7 @@ name: Build and Deploy Documentation on: push: branches: - - master + - main - dev tags: '*' pull_request: @@ -19,13 +19,18 @@ jobs: - uses: julia-actions/setup-julia@v1 with: version: '1.10' - - name: Add custom registry + - name: clone integration test tools run: | - julia --project=docs/ -e 'import Pkg; Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/QEDjl-project/registry.git"));' - julia --project=docs/ -e 'import Pkg; Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/JuliaRegistries/General"));' - - name: Install dependencies - run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + git clone --depth 1 -b dev https://github.com/QEDjl-project/QuantumElectrodynamics.jl.git /tmp/integration_test_tools/ + - name: set dev dependencies + run: | + $(julia --project=. /tmp/integration_test_tools/.ci/integTestGen/src/get_project_name_version_path.jl) + echo "CI_DEV_PKG_NAME -> $CI_DEV_PKG_NAME" + echo "CI_DEV_PKG_VERSION -> $CI_DEV_PKG_VERSION" + echo "CI_DEV_PKG_PATH -> $CI_DEV_PKG_PATH" + julia --project=docs/ /tmp/integration_test_tools/.ci/SetupDevEnv/src/SetupDevEnv.jl + julia --project=docs/ -e 'import Pkg; Pkg.instantiate()' - name: Build and deploy env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} run: julia --project=docs/ docs/make.jl diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 2aa8a58..ad3253b 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -27,12 +27,6 @@ jobs: ENV["JULIA_PKG_SERVER"] = "" Pkg.Registry.add("General") shell: julia --color=yes {0} - - name: "Add QED custom registry" - run: | - import Pkg - ENV["JULIA_PKG_SERVER"] = "" - Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/QEDjl-project/registry.git")) - shell: julia --color=yes {0} - name: "Install CompatHelper" run: | import Pkg diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..a3d77f7 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,31 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: + inputs: + lookback: + default: "3" +permissions: + actions: read + checks: read + contents: write + deployments: read + issues: read + discussions: read + packages: read + pages: read + pull-requests: read + repository-projects: read + security-events: read + statuses: read +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.TAGBOT_PRIV }} diff --git a/.github/workflows/formatter.yaml b/.github/workflows/formatter.yaml index 5cbe730..c5e31e4 100644 --- a/.github/workflows/formatter.yaml +++ b/.github/workflows/formatter.yaml @@ -9,7 +9,7 @@ jobs: - name: install Julia uses: julia-actions/setup-julia@v1 with: - version: 1.10 + version: "1.10" - name: Install Julia requirements run: julia --project=${GITHUB_WORKSPACE}/.formatting -e 'import Pkg; Pkg.instantiate()' - name: Check code style diff --git a/.gitignore b/.gitignore index a807c87..716a394 100644 --- a/.gitignore +++ b/.gitignore @@ -543,4 +543,4 @@ $RECYCLE.BIN/ # Windows shortcuts *.lnk -# End of https://www.toptal.com/developers/gitignore/api/julia,windows,macos,linux,vs,git,vim,emacs \ No newline at end of file +# End of https://www.toptal.com/developers/gitignore/api/julia,windows,macos,linux,vs,git,vim,emacs diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 540b1e7..6fd610f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -8,15 +8,17 @@ stages: stage: unit-test script: - apt update && apt install -y git - - git clone --depth 1 -b dev https://github.com/QEDjl-project/QED.jl.git /QEDjl - - julia --project=. -e 'import Pkg; Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/QEDjl-project/registry.git"));' - - julia --project=. -e 'import Pkg; Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/JuliaRegistries/General"));' + - git clone --depth 1 -b dev https://github.com/QEDjl-project/QuantumElectrodynamics.jl.git /tmp/integration_test_tools/ + - $(julia --project=. /tmp/integration_test_tools/.ci/integTestGen/src/get_project_name_version_path.jl) + - echo "CI_DEV_PKG_NAME -> $CI_DEV_PKG_NAME" + - echo "CI_DEV_PKG_VERSION -> $CI_DEV_PKG_VERSION" + - echo "CI_DEV_PKG_PATH -> $CI_DEV_PKG_PATH" - > if [[ $CI_COMMIT_BRANCH == "main" || $CI_COMMIT_REF_NAME == "main" || $CI_COMMIT_BRANCH == "dev" || $CI_COMMIT_REF_NAME == "dev" ]]; then # set name of the commit message from CI_COMMIT_MESSAGE to NO_MESSAGE, that the script does not read accidentally custom packages from the commit message of a merge commit - julia --project=. /QEDjl/.ci/SetupDevEnv/src/SetupDevEnv.jl ${CI_PROJECT_DIR}/Project.toml NO_MESSAGE + julia --project=. /tmp/integration_test_tools/.ci/SetupDevEnv/src/SetupDevEnv.jl ${CI_PROJECT_DIR}/Project.toml NO_MESSAGE else - julia --project=. /QEDjl/.ci/SetupDevEnv/src/SetupDevEnv.jl ${CI_PROJECT_DIR}/Project.toml + julia --project=. /tmp/integration_test_tools/.ci/SetupDevEnv/src/SetupDevEnv.jl ${CI_PROJECT_DIR}/Project.toml fi - julia --project=. -e 'import Pkg; Pkg.instantiate()' - julia --project=. -e 'import Pkg; Pkg.test(; coverage = true)' @@ -28,7 +30,7 @@ unit_tests_releases: extends: .untit_test_template parallel: matrix: - - JULIA_VERSION: ["1.6", "1.7", "1.8", "1.9", "1.10", "rc"] + - JULIA_VERSION: ["1.10", "1.11", "rc"] image: julia:$JULIA_VERSION unit_tests_nightly: @@ -69,15 +71,13 @@ generate_integration_tests: stage: generate_integration_test script: # extract package name - - export CI_DEPENDENCY_NAME=$(cat $CI_PROJECT_DIR/Project.toml | grep name | awk '{ print $3 }' | tr -d '"') - - echo "CI_DEPENDENCY_NAME -> $CI_DEPENDENCY_NAME" - apt update && apt install -y git - - git clone --depth 1 -b dev https://github.com/QEDjl-project/QED.jl.git /QEDjl + - git clone --depth 1 -b dev https://github.com/QEDjl-project/QuantumElectrodynamics.jl.git /QEDjl + - $(julia --project /QEDjl/.ci/integTestGen/src/get_project_name_version_path.jl) + - echo "CI_DEV_PKG_NAME -> $CI_DEV_PKG_NAME" + - echo "CI_DEV_PKG_VERSION -> $CI_DEV_PKG_VERSION" + - echo "CI_DEV_PKG_PATH -> $CI_DEV_PKG_PATH" - cd /QEDjl/.ci/integTestGen/ - # use local registry of the QED project - - julia --project=. -e 'import Pkg; Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/QEDjl-project/registry.git"));' - # needs to add General registry again, if local registry was added - - julia --project=. -e 'import Pkg; Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/JuliaRegistries/General"));' - julia --project=. -e 'import Pkg; Pkg.instantiate()' # paths of artifacts are relative to CI_PROJECT_DIR - > @@ -109,7 +109,7 @@ verify-unit-test-deps_julia1.10: stage: verify-unit-test-deps script: - apt update && apt install -y git - - git clone --depth 1 -b dev https://github.com/QEDjl-project/QED.jl.git /QEDjl + - git clone --depth 1 -b dev https://github.com/QEDjl-project/QuantumElectrodynamics.jl.git /QEDjl - > if [[ $CI_COMMIT_BRANCH == "main" || $CI_COMMIT_REF_NAME == "main" || $CI_COMMIT_BRANCH == "dev" || $CI_COMMIT_REF_NAME == "dev" ]]; then # does not check for custom package URLs on the main and dev branch diff --git a/CHANGELOG.md b/CHANGELOG.md index 369d9c8..7790b1e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,8 +1,44 @@ # Changelog +## Version 0.3.0 + +[diff since 0.2.2](https://github.com/QEDjl-project/QEDbase.jl/compare/v0.2.2...v0.3.0) + +This release contains the renaming of `QuantumElectrodynamics.jl` (previously `QED.jl`) and improvements of the interface and functionality of `AbstractPhaseSpacePoint` for spins and polarizations of particles. This includes the addition of new types `SyncedSpin` and `SyncedPol`. + +Since Julia 1.10 is now the LTS version, this release also drops versions below 1.10 from support. + +### Breaking Changes + +The `multiplicity` function was changed in [#112](https://github.com/QEDjl-project/QEDbase.jl/pull/112) to accept `AbstractProcessDefinition`s instead of raw spins or polarizations, breaking compatibility if this function was used. The function was used in [`QEDprocesses.jl`](https://github.com/QEDjl-project/QEDprocesses.jl), with the incompatibility fixed in the respective new version `QEDprocesses.jl-v0.3.0`. + +### New features + +- [#106](https://github.com/QEDjl-project/QEDbase.jl/pull/106): Add spin and polarization interface for `AbstractProcessDefinition`, namely functions `spin_pol`, `incoming_spin_pol` and `outgoing_spin_pol` +- [#112](https://github.com/QEDjl-project/QEDbase.jl/pull/112): Add `SyncedSpin` and `SyncedPol` types, representing a synchronization of the spins or polarizations across multiple particles in the same scattering process +- [#118](https://github.com/QEDjl-project/QEDbase.jl/pull/118): Add `spin_pols_iter`, returning an iterator for a given `AbstractProcessDefinition`, yielding all possible spin and polarization combinations for the process, including support for the `SyncedSpin` and `SyncedPol` types +- [#129](https://github.com/QEDjl-project/QEDbase.jl/pull/129): Add coordinate transformation interface + +### Fixes + +- [#110](https://github.com/QEDjl-project/QEDbase.jl/pull/110): Remove `mul` from exported functions and rename, fixing a naming collision and unnecessary export +- [#111](https://github.com/QEDjl-project/QEDbase.jl/pull/111): Fix GPU incompatibility of the `momenta` implementation +- [#115](https://github.com/QEDjl-project/QEDbase.jl/pull/115): Set QEDjl-project dependencies to dev versions for documentation build and deploy job +- [#116](https://github.com/QEDjl-project/QEDbase.jl/pull/116): Fix docs build and deployment CI job + +### Maintenance + +- [#108](https://github.com/QEDjl-project/QEDbase.jl/pull/108): Separate implementations from interfaces in the repository file structure +- [#120](https://github.com/QEDjl-project/QEDbase.jl/pull/120): Apply renaming of `QuantumElectrodynamics.jl` (previously `QED.jl`) +- [#123](https://github.com/QEDjl-project/QEDbase.jl/pull/123): Remove custom QEDjl-project registry +- [#124](https://github.com/QEDjl-project/QEDbase.jl/pull/124): Add Julia TagBot +- [#130](https://github.com/QEDjl-project/QEDbase.jl/pull/130): Drop Julia 1.9 and below from support and CI +- [#125](https://github.com/QEDjl-project/QEDbase.jl/pull/125): Update integration test suite +- [#127](https://github.com/QEDjl-project/QEDbase.jl/pull/127): Add lots of docs and tutorials + ## Version 0.2.2 -[diff since 0.2.0](https://github.com/QEDjl-project/QEDbase.jl/compare/release-0.2.0...release-0.2.2) +[diff since 0.2.0](https://github.com/QEDjl-project/QEDbase.jl/compare/v0.2.0...v0.2.2) This release adds some convenience overloads to existing functions, some code maintenance and small fixes. @@ -35,9 +71,9 @@ No further breaking or non-breaking changes. ## Version 0.2.0 -[diff since 0.1.6](https://github.com/QEDjl-project/QEDbase.jl/compare/release-0.1.6...release-0.2.0) +[diff since 0.1.6](https://github.com/QEDjl-project/QEDbase.jl/compare/v0.1.6...v0.2.0) -This release is part of the restructuring processes of QED.jl (see https://github.com/QEDjl-project/QED.jl/issues/35 for details). +This release is part of the restructuring processes of QED.jl (see https://github.com/QEDjl-project/QuantumElectrodynamics.jl/issues/35 for details). It is a breaking release, indicated by the bumped minor version because we don't have a major version for this yet. @@ -75,6 +111,8 @@ fixes ## Version 0.1.6 +[diff since 0.1.6](https://github.com/QEDjl-project/QEDbase.jl/compare/1f8b2d2340e5a92c700bff458bfc385b8904448f...v0.1.6) + This is a maintenance release, which resolves, among others, some issues with the git history. ### Breaking changes @@ -94,7 +132,7 @@ No new features. ## Version 0.1.5 -[diff since 0.1.4](https://github.com/QEDjl-project/QEDbase.jl/compare/0c70f66...release-0.1.5) +[diff since 0.1.4](https://github.com/QEDjl-project/QEDbase.jl/compare/v0.1.4...v0.1.5) ### Breaking changes diff --git a/Project.toml b/Project.toml index 3bb6429..487cb89 100644 --- a/Project.toml +++ b/Project.toml @@ -7,7 +7,7 @@ authors = [ "Tom Jungnickel", "Anton Reinhard", ] -version = "0.2.2" +version = "0.3.0" [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" @@ -25,9 +25,10 @@ ArgCheck = "2.3.0" ConstructionBase = "1" DocStringExtensions = "0.8.5, 0.9" PhysicalConstants = "0.2.1" +QEDcore = "0.1.1" SimpleTraits = "0.9.4" StaticArrays = "1.2.13" -julia = "1.6" +julia = "1.10" [extras] QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e" diff --git a/README.md b/README.md index e310949..ec9c36f 100644 --- a/README.md +++ b/README.md @@ -1,115 +1,74 @@ # QEDbase - -[![Doc Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://qedjl-project.github.io/QEDbase.jl/main) +[![Doc Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://qedjl-project.github.io/QEDbase.jl/stable) [![Doc Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://qedjl-project.github.io/QEDbase.jl/dev) -[![Build Status](https://gitlab.hzdr.de/qedjl/QEDbase.jl/badges/main/pipeline.svg)](https://gitlab.hzdr.de/qedjl/QEDbase.jl/pipelines) -[![Coverage](https://gitlab.hzdr.de/qedjl/QEDbase.jl/badges/main/coverage.svg)](https://gitlab.hzdr.de/qedjl/QEDbase.jl/commits/main) -[![Coverage](https://codecov.io/gh/qedjl/QEDbase.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/qedjl/QEDbase.jl) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) -This is `QEDbase.jl`, a julia package which provides the general data structures for calculations in relativistic particle physics. +`QEDbase.jl` is a foundational package within the [`QuantumElectrodynamics.jl`](https://qedjl-project.github.io/QuantumElectrodynamics.jl/dev/) +library. It provides essential interfaces and building blocks for the computation of +quantum electrodynamics (QED) processes, facilitating interoperability with other packages in the suite. + +For a detailed explanation of the integration with other `QuantumElectrodynamics.jl` packages, +please refer to the [documentation of `QuantumElectrodynamics.jl`](https://qedjl-project.github.io/QuantumElectrodynamics.jl/dev/). -This package is part of the `QED.jl` library. For the description of the interoperability with other packages of `QED.jl` see [docs](www.docs-to-qed.jl). +## Main interfaces -## Current features +- **Dirac Tensors**: Types that facilitate operations involving Dirac matrices and spinors. +- **Lorentz Vectors**: Types that facilitate operations involving Lorentz vectors. +- **Particle Representation**: Define particles with mass, charge, and other physical properties. +- **Computation Models**: Interfaces for implementing various physical models (e.g., perturbative or strong-field QED) for calculations. +- **Scattering Processes**: Generic descriptions of scattering processes for use in QED calculations. +- **Probabilities and Cross Sections**: Core components for calculating differential probabilities and cross-sections. +- **Phase Space Descriptions**: Utility functions to define and manage phase spaces and related points. -- Generic interface for Lorentz vectors -- concrete implementations of static and mutable Lorentz-vector/four-momentum types -- general Dirac bi-spinors, its adjoint counterpart as well as Dirac matrices -- particle spinors, i.e. solutions of Dirac's equation in momentum space -- Dirac's gamma matrices +To lern how to implement these interfaces, please refer to the tutorials in the +[documentation](https://qedjl-project.github.io/QEDbase.jl/dev). ## Installation -To install the current stable version of `QEDbase.jl` you may use the standard julia package manager within the julia REPL +To install the latest stable version of `QEDbase.jl`, use the Julia package manager within the REPL: ```julia julia> using Pkg - julia> Pkg.add("QEDbase") ``` -or you use the Pkg prompt by hitting `]` within the Julia REPL and then type +Alternatively, you can enter the Pkg prompt by pressing `]` in the Julia REPL and then run: ```julia -(@v1.7) pkg> add QEDbase +pkg> add QEDbase ``` -To install the locally downloaded package on Windows, change to the parent directory and type within the Pkg prompt - -```julia -(@v1.7) pkg> add ./QEDbase.jl -``` - -## Quickstart -#### Four momentum -One can define a static four momentum component wise: - -```julia -julia> using QEDbase; using QEDcore - -juila> mass = rand()*10 - -julia> px,py,pz = rand(3) - -julia> E = sqrt(px^2 + py^2 + pz^2 + mass^2) # on-shell condition - -julia> P = SFourMomentum(E,px,py,pz) -``` - -Such `SFourMomentum` behaves like a four element static array (with all the standard arithmetics), but with the `dot` product exchanged with the Minkowski product - -```julia -julia> @assert dot(P,P) == P*P == getMass2(P) == P[1]^2 - sum(P[1:].^2) -``` - -Furthermore, the Lorentz-vector interface provides a lot of properties for such a `SFourMomentum`, e.g. - -```julia - -julia> @assert isapprox(getRapidity(mom), 0.5*log((E+pz)/(E-pz))) -julia> @assert isapprox(getPlus(mom), 0.5*(E+pz)) -julia> @assert isapprox(getPerp(mom), px^2 + py^2) -``` - -and a lot more (see [here](www.docs-to-the-lorentz-interface-getter.jl) for a complete list). There is also a mutable version of a four vector in `QEDcore.jl`, where the Lorentz-vector interface provides setters to different properties as well (see [here](www.docs-to-the-lorentz-interface-setter.jl) for details). - -## Testing - -After installation it might be necessary to check if everything works properly. For that you can run the unittests by typing within the julia REPL - -```julia -julia> using Pkg -julia> using QEDbase -julia> Pkg.test("QEDbase") - -... +## Contributing -Testing Running tests... -Test Summary: | Pass Total -QEDbase.jl | 468 468 - Testing QEDbase tests passed -``` +Contributions are welcome! If you'd like to report a bug, suggest an enhancement, or contribute +code, please feel free to open an issue or submit a pull request. -If you see the last line, you can assume that `QEDbase.jl` works properly for you. +To ensure consistency across the `QuantumElectrodynamics.jl` ecosystem, we encourage all contributors +to review the [QuantumElectrodynamics.jl contribution guide](https://qedjl-project.github.io/QuantumElectrodynamics.jl/stable/dev_guide/#Development-Guide). -## Contributing +## Credits and contributors -If you want to contribute to `QEDbase.jl` feel free to do so by opening an issue or send me a pull request. In order to keep the packages within `QED.jl` coherent, consider visiting the general [contribution guide of `QED.jl`](www.contribution-of-qed.jl). +This work was partly funded by the Center for Advanced Systems Understanding (CASUS) that +is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon +Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget +approved by the Saxon State Parliament. -## Credits and contributers +The core code of the package `QEDbase.jl` is developed by a small team at the Center for +Advanced Systems Understanding ([CASUS](https://www.casus.science)), namely -This work was partly funded by the Center for Advanced Systems Understanding (CASUS) that is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament. +### Core Contributors -The core code of the package `QEDbase.jl` is developed by a small team at the Center for Advanced Systems Understanding ([CASUS](https://www.casus.science)), namely +- **Uwe Hernandez Acosta** (CASUS/HZDR, [u.hernandez@hzdr.de](mailto:u.hernandez@hzdr.de)) +- **Anton Reinhard** (CASUS/HZDR) +- **Simeon Ehrig** (CASUS/HZDR) +- **Klaus Steiniger** (CASUS/HZDR) -- Uwe Hernandez Acosta (u.hernandez@hzdr.de) +### Former Contributors -This package would not be possible without many contributions done from the community as well. For that we want send big thanks to: +- **Tom Jungnickel** -- my Mate supplier -- the guy who will show me how to include the most recent contributers on github +We extend our sincere thanks to all contributors who have supported this project. ## License diff --git a/docs/Project.toml b/docs/Project.toml index a82680f..dc6cb76 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,6 +2,9 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93" QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e" +QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad" diff --git a/docs/make.jl b/docs/make.jl index afc3b91..61960fd 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,42 +7,110 @@ project_path = Base.Filesystem.joinpath(Base.Filesystem.dirname(Base.source_path Pkg.develop(; path=project_path) using Documenter +using Literate +using DocumenterInterLinks +using DocumenterCitations + using QEDbase +using QEDcore +using QEDprocesses -#DocMeta.setdocmeta!(QEDbase, :DocTestSetup, :(using QEDbase); recursive=true) +bib = CitationBibliography(joinpath(@__DIR__, "Bibliography.bib")) -# TODO: remove before release -Pkg.add(; url="https://github.com/QEDjl-project/QEDcore.jl", rev="dev") +# some paths for links +readme_path = joinpath(project_path, "README.md") +index_path = joinpath(project_path, "docs/src/index.md") +license_path = "https://github.com/QEDjl-project/QEDbase.jl/blob/main/LICENSE" -using DocumenterCitations +# Copy README.md from the project base folder and use it as the start page +open(readme_path, "r") do readme_in + readme_string = read(readme_in, String) + + # replace relative links in the README.md + readme_string = replace(readme_string, "[MIT](LICENSE)" => "[MIT]($(license_path))") + + open(index_path, "w") do readme_out + write(readme_out, readme_string) + end +end + +# setup interlinks +links = InterLinks("QEDcore" => "https://qedjl-project.github.io/QEDcore.jl/dev/") +# setup Bibliography bib = CitationBibliography(joinpath(@__DIR__, "Bibliography.bib")) +# setup examples using Literate.jl +literate_paths = [ + Base.Filesystem.joinpath(project_path, "docs/src/tutorial/lorentz_vectors.jl"), + Base.Filesystem.joinpath(project_path, "docs/src/tutorial/particle.jl"), + Base.Filesystem.joinpath(project_path, "docs/src/tutorial/particle_stateful.jl"), + Base.Filesystem.joinpath(project_path, "docs/src/tutorial/process.jl"), + Base.Filesystem.joinpath(project_path, "docs/src/tutorial/phase_space_point.jl"), +] + +tutorial_output_dir = joinpath(project_path, "docs/src/generated/") +!ispath(tutorial_output_dir) && mkdir(tutorial_output_dir) +@info "Literate: create temp dir at $tutorial_output_dir" +tutorial_output_dir_name = splitpath(tutorial_output_dir)[end] + pages = [ "Home" => "index.md", - "Dirac Tensors" => "dirac_tensors.md", - "Lorentz Vectors" => "lorentz_vectors.md", - "Four Momenta" => "four_momenta.md", - "Library" => [ + "Tutorials" => [ + #"Dirac Tensors" => "dirac_tensors.md", + "Lorentz Vectors" => joinpath(tutorial_output_dir_name, "lorentz_vectors.md"), + "Particles" => joinpath(tutorial_output_dir_name, "particle.md"), + "Stateful Particles" => + joinpath(tutorial_output_dir_name, "particle_stateful.md"), + "Scattering Process" => joinpath(tutorial_output_dir_name, "process.md"), + "Phase Space Points" => + joinpath(tutorial_output_dir_name, "phase_space_point.md"), + ], + "API reference" => [ "Contents" => "library/outline.md", - "API" => "library/api.md", + "Lorentz vectors" => "library/lorentz_vector.md", + "Dirac tensors" => "library/dirac_objects.md", + "Particles" => "library/particles.md", + "Scattering process" => "library/process.md", + "Phase space description" => "library/phase_space.md", + "Probability and cross section" => "library/cross_section.md", "Function index" => "library/function_index.md", ], "refs.md", ] -makedocs(; - modules=[QEDbase], - checkdocs=:exports, - authors="Uwe Hernandez Acosta", - repo=Documenter.Remotes.GitHub("QEDjl-project", "QEDbase.jl"), - sitename="QEDbase.jl", - format=Documenter.HTML(; - prettyurls=get(ENV, "CI", "false") == "true", - canonical="https://qedjl-project.gitlab.io/QEDbase.jl", - assets=String[], - ), - pages=pages, - plugins=[bib], -) +try + # generate markdown files with Literate.jl + for file in literate_paths + Literate.markdown(file, tutorial_output_dir; documenter=true) + end + + # geneate docs with Documenter.jl + makedocs(; + modules=[QEDbase], + checkdocs=:exports, + authors="Uwe Hernandez Acosta", + repo=Documenter.Remotes.GitHub("QEDjl-project", "QEDbase.jl"), + sitename="QEDbase.jl", + format=Documenter.HTML(; + prettyurls=get(ENV, "CI", "false") == "true", + canonical="https://qedjl-project.gitlab.io/QEDbase.jl", + assets=String[], + mathengine=Documenter.MathJax2(), + collapselevel=1, + # TODO: workaround + # should be fixed: https://github.com/QEDjl-project/QEDbase.jl/issues/4 + size_threshold_ignore=["index.md"], + ), + pages=pages, + plugins=[bib, links], + ) +finally + # doing some garbage collection + @info "GarbageCollection: remove generated landing page" + rm(index_path) + @info "GarbageCollection: remove generated tutorial files" + rm(tutorial_output_dir; recursive=true) +end + deploydocs(; repo="github.com/QEDjl-project/QEDbase.jl.git", push_preview=false) diff --git a/docs/src/dirac_tensors.md b/docs/src/dirac_tensors.md deleted file mode 100644 index e946da7..0000000 --- a/docs/src/dirac_tensors.md +++ /dev/null @@ -1,3 +0,0 @@ -# Dirac tensors - -We model dirac_tensors as concrete subtypes from AbstractDiracMatrix and AbstractDiracVector, respectively. diff --git a/docs/src/four_momenta.md b/docs/src/four_momenta.md deleted file mode 100644 index 3b92d8e..0000000 --- a/docs/src/four_momenta.md +++ /dev/null @@ -1,7 +0,0 @@ -# Four Momenta - -Within QEDbase, we model four-momenta with a specialized LorentzVector, where the components are forced to be float: - -``` Julia - FourMomentum <:AbstractLorentzVector{::Float64} -``` diff --git a/docs/src/index.md b/docs/src/index.md deleted file mode 100644 index 13d7c86..0000000 --- a/docs/src/index.md +++ /dev/null @@ -1,11 +0,0 @@ -```@meta -CurrentModule = QEDbase -``` - -# QEDbase - -Documentation for [QEDbase](https://gitlab.hzdr.de/QEDjl/QEDbase.jl). This library provides (more or less) simple implementations of Lorentz vectors, Dirac spinors and Dirac matrices. These are usually used in order to do calculations in particle physics, e.g. they are part of the set of Feynman rules for a given quantum field theory (see [Peskin:1995ev](@cite)). - -```@autodocs -Modules = [QEDbase] -``` diff --git a/docs/src/library/api.md b/docs/src/library/api.md deleted file mode 100644 index cc3ffff..0000000 --- a/docs/src/library/api.md +++ /dev/null @@ -1,3 +0,0 @@ -# QEDbase - -:warning: Under construction diff --git a/docs/src/library/cross_section.md b/docs/src/library/cross_section.md new file mode 100644 index 0000000..40d71aa --- /dev/null +++ b/docs/src/library/cross_section.md @@ -0,0 +1,27 @@ +# Probability and Cross Section + +## Interface + +```@meta +CurrentModule = QEDbase +``` + +```@docs +_incident_flux +``` + +## Differential and total probability + +```@docs +differential_probability +unsafe_differential_probability +total_probability +``` + +## Differential and total cross section + +```@docs +differential_cross_section +unsafe_differential_cross_section +total_cross_section +``` diff --git a/docs/src/library/dirac_objects.md b/docs/src/library/dirac_objects.md new file mode 100644 index 0000000..39b737b --- /dev/null +++ b/docs/src/library/dirac_objects.md @@ -0,0 +1,6 @@ +# Dirac Tensors + +```@docs +AbstractDiracVector +AbstractDiracMatrix +``` diff --git a/docs/src/library/function_index.md b/docs/src/library/function_index.md index 2b38eed..f5ce0a4 100644 --- a/docs/src/library/function_index.md +++ b/docs/src/library/function_index.md @@ -1,5 +1,5 @@ ### [Index](@id main-index) ```@index -Pages = ["api.md", "function_index.md"] + ``` diff --git a/docs/src/library/lorentz_vector.md b/docs/src/library/lorentz_vector.md new file mode 100644 index 0000000..c3e91ed --- /dev/null +++ b/docs/src/library/lorentz_vector.md @@ -0,0 +1,94 @@ +# Lorentz Vector + +## Interface registry + +```@docs +register_LorentzVectorLike +``` + +## Built-in Lorentz vector types + +```@docs +AbstractLorentzVector +AbstractFourMomentum +``` + +## Accessor functions + +```@docs +minkowski_dot +mdot +getT +getX +getY +getZ +getMagnitude2 +getMag2 +getMagnitude +getMag + +getInvariantMass2 +getMass2 +getInvariantMass +getMass +getE +getEnergy +getPx +getPy +getPz +getBeta +getGamma + +getTransverseMomentum2 +getPt2 +getPerp2 +getTransverseMomentum +getPt +getPerp +getTransverseMass2 +getMt2 +getTransverseMass +getMt +getRapidity + +getRho2 +getRho +getTheta +getCosTheta +getPhi +getCosPhi +getSinPhi +getPlus +getMinus +``` + +## Setter functions + +```@docs +setE! +setEnergy! +setPx! +setPy! +setPz! +setTheta! +setCosTheta! +setRho! +setPhi! +setPlus! +setMinus! + +setTransverseMomentum! +setPerp! +setPt! +setTransverseMass! +setMt! +setRapidity! +``` + +## Utilities + +```@docs +isonshell +assert_onshell +AbstractCoordinateTransformation +``` diff --git a/docs/src/library/model.md b/docs/src/library/model.md new file mode 100644 index 0000000..f2ead35 --- /dev/null +++ b/docs/src/library/model.md @@ -0,0 +1,6 @@ +# Computational Model Interface + +```@docs +AbstractModelDefinition +fundamental_interaction_type +``` diff --git a/docs/src/library/outline.md b/docs/src/library/outline.md index ccd5160..cd6741d 100644 --- a/docs/src/library/outline.md +++ b/docs/src/library/outline.md @@ -1,5 +1,5 @@ ## API Outline ```@contents -Pages = ["api.md", "function_index.md"] +Pages = ["lorentz_vector.md", "function_index.md"] ``` diff --git a/docs/src/library/particles.md b/docs/src/library/particles.md new file mode 100644 index 0000000..d332c50 --- /dev/null +++ b/docs/src/library/particles.md @@ -0,0 +1,21 @@ +# Particle Interface + +## Mandatory Interface + +```@docs +AbstractParticle +AbstractParticleType +mass +charge +base_state +propagator +``` + +## Optional Interface + +```@docs +is_fermion +is_boson +is_particle +is_anti_particle +``` diff --git a/docs/src/library/phase_space.md b/docs/src/library/phase_space.md new file mode 100644 index 0000000..6f699ac --- /dev/null +++ b/docs/src/library/phase_space.md @@ -0,0 +1,42 @@ +# Phase Space Description + +## Frames and Coordinate Systems + +```@docs +AbstractCoordinateSystem +AbstractFrameOfReference +AbstractPhasespaceDefinition +``` + +## Stateful Particles + +```@docs +AbstractParticleStateful +momentum +``` + +## Phasespace Points + +### Types + +```@docs +AbstractPhaseSpacePoint +AbstractInPhaseSpacePoint +AbstractOutPhaseSpacePoint +``` + +### Interface + +```@docs +process +model +phase_space_definition +``` + +### Convenience Functions + +```@docs +particle_direction +particle_species +momenta +``` diff --git a/docs/src/library/process.md b/docs/src/library/process.md new file mode 100644 index 0000000..dedb130 --- /dev/null +++ b/docs/src/library/process.md @@ -0,0 +1,69 @@ +# Scattering Process Interface + +## Process Interface + +```@docs +AbstractProcessDefinition +incoming_particles +outgoing_particles +incoming_spin_pols +outgoing_spin_pols +``` + +## Particle Directions + +```@docs +ParticleDirection +Incoming +Outgoing +UnknownDirection +is_incoming +is_outgoing +``` + +## Spins and Polarizations + +```@docs +AbstractSpinOrPolarization +spin_pols_iter +``` + +### Spins + +```@docs +AbstractSpin +AbstractDefiniteSpin +AbstractIndefiniteSpin +SpinUp +SpinDown +AllSpin +SyncedSpin +``` + +### Polarization + +```@docs +AbstractPolarization +AbstractDefinitePolarization +AbstractIndefinitePolarization +PolarizationX +PolX +PolarizationY +PolY +AllPolarization +AllPol +SyncedPolarization +``` + +## Utility Functions + +```@docs +number_incoming_particles +number_outgoing_particles +particles +number_particles +spin_pols +multiplicity +incoming_multiplicity +outgoing_multiplicity +``` diff --git a/docs/src/library/utility.md b/docs/src/library/utility.md new file mode 100644 index 0000000..fcd64f3 --- /dev/null +++ b/docs/src/library/utility.md @@ -0,0 +1,19 @@ +```@meta +CurrentModule = QEDbase +``` + +# Utility + +## Errors + +```@docs + +InvalidInputError +RegistryError +``` + +## Helper functions + +```@docs +_as_svec +``` diff --git a/docs/src/lorentz_vectors.md b/docs/src/lorentz_vectors.md deleted file mode 100644 index 29a6474..0000000 --- a/docs/src/lorentz_vectors.md +++ /dev/null @@ -1,3 +0,0 @@ -# Lorentz vectors - -In oder to model Lorentz vectors, i.e. elements of a Minkowski space, we use a generic type called `LorentzVector{T}`. diff --git a/docs/src/phase_space_point.md b/docs/src/phase_space_point.md new file mode 100644 index 0000000..5c3cfe2 --- /dev/null +++ b/docs/src/phase_space_point.md @@ -0,0 +1,41 @@ +# What is a Phase Space Point? + +In high-energy physics, a **phase space point** represents a specific configuration of +particles involved in a particle interaction or scattering process. More specifically, +it refers to the set of momenta (and possibly other quantum states) of all the incoming +and outgoing particles at a particular instant in the process. + +For example, in a scattering process such as electron-positron annihilation: + +```math +e^+ + e^- \to \gamma + \gamma +``` + +a phase space point would include the momenta (four-momenta, to be precise) of both the +incoming particles (the positron \(e^+\) and electron \(e^-\)) and the outgoing particles +(the two photons \(\gamma\)). Each phase space point can be thought of as a snapshot of +the physical state of the process at a given point in time. + +## Why is it important? + +In particle physics, calculations often involve integrating over all possible configurations +of momenta that satisfy certain conservation laws (like energy and momentum conservation). +The space of all these configurations is called the **phase space**, and each valid configuration +is a phase space point. These points play a central role in computing things like cross sections, +decay rates, and scattering amplitudes. + +For any given process, we need a way to: + +- Represent the particles and their momenta in the phase space. +- Ensure that the representation is consistent with the underlying physics (e.g., the correct + number and type of particles, momentum conservation). +- Access the momenta of individual particles or compute derived quantities. + +In [this tutorial], we show how to create and use an `ExamplePhaseSpacePoint`, +following the interface specification given by [`AbstractPhaseSpacePoint']. We'll be using +particles and momenta from libraries like `QEDcore` and `QEDprocesses`, focusing on the +electron-positron annihilation process. + +## Reference implementation + +TBW diff --git a/docs/src/tutorial/lorentz_vectors.jl b/docs/src/tutorial/lorentz_vectors.jl new file mode 100644 index 0000000..e85b523 --- /dev/null +++ b/docs/src/tutorial/lorentz_vectors.jl @@ -0,0 +1,117 @@ +# # Tutorial: Custom Lorentz Vectors +# +# [Lorentz vectors](https://en.wikipedia.org/wiki/Four-vector), which represent elements of +# Minkowski space, are fundamental mathematical objects in high-energy physics for describing +# the kinematic properties of particles. In particular, they are used to express four-momentum and other +# kinematic quantities. +# +# `QEDbase.jl` provides a flexible interface for working with Lorentz vectors and four-momenta, +# allowing seamless integration with custom types that represent Lorentz vectors. This enables +# users to extend the core functionality of the package to their own data structures, while +# taking advantage of optimized kinematic computations. +# +# ## Defining a Custom Lorentz Vector Type +# +# The `LorentzVector` interface in `QEDbase.jl` is designed to be extendable. By implementing +# a simple API, you can make any custom type behave like a Lorentz vector and unlock access +# to various kinematic functions in the package. + +# ### Step 1: Define Your Custom Type +# +# Suppose you want to define a custom type that behaves like a Lorentz vector. Start by +# creating the type to hold the Cartesian components of the vector: +using QEDbase + +struct CustomLorentzVector + t::Float64 # time-like component + x::Float64 # x-component of spatial vector + y::Float64 # y-component of spatial vector + z::Float64 # z-component of spatial vector +end + +# ### Step 2: Implement the Lorentz Vector API +# +# To connect your custom type to the `LorentzVector` interface, you'll need to implement the +# required getters for the time and spatial components. These functions extract the respective +# components (`t`, `x`, `y`, and `z`) from your type: + +QEDbase.getT(lv::CustomLorentzVector) = lv.t +QEDbase.getX(lv::CustomLorentzVector) = lv.x +QEDbase.getY(lv::CustomLorentzVector) = lv.y +QEDbase.getZ(lv::CustomLorentzVector) = lv.z + +# ### Step 3: Register Your Custom Type as a `LorentzVectorLike` +# +# Once you've implemented the required API, you can register your custom type as a `LorentzVectorLike` by calling the following function: + +register_LorentzVectorLike(CustomLorentzVector) + +# If any required functions are missing or incorrectly implemented, this will raise a `RegistryError` that provides details on what needs to be fixed. + +# ### Step 4: Using the Custom Lorentz Vector +# +# With your custom type registered, you can now use it as a Lorentz vector in various kinematic calculations. +# For example, the `CustomLorentzVector` can be used to calculate angles or rapidity: + +L = CustomLorentzVector(2.0, 1.0, 0.0, -1.0) + +println("Polar angle (theta)): ", getTheta(L)) +println("Rapidity: ", getRapidity(L)) + +# ### Optional: Mutable Types +# +# If you need to modify the components of the Lorentz vector after its creation, you can +# implement setter functions for your custom type. This will allow your type to be mutable +# and provide access to additional mutating functions within the `QEDbase.jl` ecosystem. +# +# First: define another custom Lorentz vector type, but make it a mutable struct + +mutable struct CustomMutableLorentzVector + t::Float64 # time-like component + x::Float64 # x-component of spatial vector + y::Float64 # y-component of spatial vector + z::Float64 # z-component of spatial vector +end + +# Second: define the accessor functions for your mutable Lorentz vector + +QEDbase.getT(lv::CustomMutableLorentzVector) = lv.t +QEDbase.getX(lv::CustomMutableLorentzVector) = lv.x +QEDbase.getY(lv::CustomMutableLorentzVector) = lv.y +QEDbase.getZ(lv::CustomMutableLorentzVector) = lv.z + +# +# Third: enable mutability, implement the following setter methods +QEDbase.setT!(lv::CustomMutableLorentzVector, new_t) = (lv.t = new_t) +QEDbase.setX!(lv::CustomMutableLorentzVector, new_x) = (lv.x = new_x) +QEDbase.setY!(lv::CustomMutableLorentzVector, new_y) = (lv.y = new_y) +QEDbase.setZ!(lv::CustomMutableLorentzVector, new_z) = (lv.z = new_z) + +# finally register your mutable Lorentz vector + +register_LorentzVectorLike(CustomMutableLorentzVector) + +# Once these setter methods are defined, your type will support mutable operations. You can +# then register your type and take advantage of mutating functions provided by `QEDbase.jl`. + +# #### Example: Updating Rapidity +# +# After defining your custom type as mutable, you can use functions like `setRapidity!` to +# update the rapidity and other kinematic properties of the vector in place: + +L = CustomMutableLorentzVector(2.0, 1.0, 0.0, -1.0) + +println("rapidity before: ", getRapidity(L)) +setRapidity!(L, 0.5) # Updates the components accordingly +println("rapidity after: ", getRapidity(L)) + +# ## Summary +# +# In this tutorial, we learned how to: +# +# 1. Define a custom type to represent a Lorentz vector. +# 2. Implement the required API functions to integrate the type into the `LorentzVector` interface of `QEDbase.jl`. +# 3. Register the type as `LorentzVectorLike` to enable its use in kinematic calculations. +# 4. Optionally, we explored how to make the custom type mutable and enable modification of its components. + +# By following these steps, you can extend `QEDbase.jl` to work with your own Lorentz vector types while benefiting from the library's optimized calculations. diff --git a/docs/src/tutorial/particle.jl b/docs/src/tutorial/particle.jl new file mode 100644 index 0000000..f4038df --- /dev/null +++ b/docs/src/tutorial/particle.jl @@ -0,0 +1,79 @@ +# # [Tutorial: Implementing a New Particle Type](@id tutorial_particle) +# +# In this tutorial, we will implement a new particle type following the interface specification +# of `AbstractParticle` as used in `QuantumElectrodynamics.jl`. The `AbstractParticle` system +# has two main components: +# +# 1. **Static functions**: These functions determine the type of particle (fermion, boson, particle, +# or anti-particle). They provide compile-time information about the particle type. +# 2. **Property functions**: These functions define the physical properties of the particle, such as mass and charge. +# +# ## Define a New Particle Species +# +# To start, we define a new concrete particle type that is a subtype of `AbstractParticleType`. +# For this example, we'll implement a new particle called `Muon`. +using QEDbase + +# Define a new particle species as a subtype of AbstractParticleType +struct Muon <: AbstractParticleType end + +# Since Muon is a subtype of AbstractParticleType, it is a singleton (no stateful properties). +# Implementing static functions for particle classification +is_fermion(::Muon) = true # Muon is a fermion +is_boson(::Muon) = false # Muon is not a boson +is_particle(::Muon) = true # Muon is a particle (not an anti-particle) +is_anti_particle(::Muon) = false # Muon is not an anti-particle + +# ## Define Physical Properties +# +# Next, we need to define the required property functions for `mass` and `charge`. +# These functions return the mass and charge of the `Muon`. + +# Define the physical properties of the Muon +mass(::Muon) = 105.66 # Muon mass in MeV/c^2 +charge(::Muon) = -1.0 # Muon has a charge of -1 (same as electron) + +# ## Anti-Particle Implementation (Optional) +# +# If we want to define the anti-particle of the muon (anti-muon), we can do that by defining +# another particle type and changing the static functions accordingly. + +# Define the anti-particle of the Muon +struct AntiMuon <: AbstractParticleType end + +# Implement static functions for the AntiMuon +is_fermion(::AntiMuon) = true # AntiMuon is also a fermion +is_boson(::AntiMuon) = false # AntiMuon is not a boson +is_particle(::AntiMuon) = false # AntiMuon is not a regular particle (it's an anti-particle) +is_anti_particle(::AntiMuon) = true # AntiMuon is an anti-particle + +# Define the physical properties for the AntiMuon +mass(::AntiMuon) = 105.66 # AntiMuon has the same mass as Muon +charge(::AntiMuon) = 1.0 # AntiMuon has the opposite charge of Muon + +# ## Example Usage +# +# Now we have both the `Muon` and its anti-particle (`AntiMuon`) implemented with the +# required interface functions. This makes these particles usable in `QuantumElectrodynamics.jl` processes. + +# Create a muon and an anti-muon +mu = Muon() +anti_mu = AntiMuon() + +# Access particle properties +println("Is Muon a fermion? ", is_fermion(mu)) # true +println("Is AntiMuon an anti-particle? ", is_anti_particle(anti_mu)) # true +println("Muon mass: ", mass(mu), " MeV/c^2") # 105.66 MeV/c^2 +println("AntiMuon charge: ", charge(anti_mu)) # +1.0 + +# ## Summary +# +# In this tutorial, we have demonstrated how to implement a new particle type in `QuantumElectrodynamics.jl`. +# By following the interface specification for `AbstractParticle`, we created: +# +# 1. A new particle (`Muon`) and defined its static functions (`is_fermion`, `is_boson`, etc.), +# 2. The required property functions (`mass` and `charge`), +# 3. An anti-particle (`AntiMuon`) with opposite charge and appropriate static functions. +# +# This setup makes these particles compatible with the rest of the `QuantumElectrodynamics.jl` framework +# and ready for use in particle physics simulations. diff --git a/docs/src/tutorial/particle_stateful.jl b/docs/src/tutorial/particle_stateful.jl new file mode 100644 index 0000000..a545883 --- /dev/null +++ b/docs/src/tutorial/particle_stateful.jl @@ -0,0 +1,94 @@ +# # Tutorial: Implementing a Stateful Particle +# +# In this tutorial, we will extend particle types (like our `Muon` and `AntiMuon` from [the previous tutorial](@ref tutorial_particle)) +# to implement **stateful particles**. A **stateful particle** is a particle that has runtime +# properties such as momentum, in addition to its compile-time species (like mass and charge). +# +# This tutorial will guide you through the implementation of `AbstractParticleStateful` by +# creating stateful versions of particle types. +# +# ## Define a Stateful Particle Type +# +# To start, we need to define a concrete subtype of `AbstractParticleStateful`. A stateful particle must store information such as: +# +# - **Direction**: Incoming or outgoing (`ParticleDirection`). +# - **Species**: The type of particle (`AbstractParticleType`). +# - **Momentum**: The four-momentum of the particle (`AbstractFourMomentum`). +# +# We will define a general `ExampleParticleStateful` type that can take any particle species as a type parameter. +# This type will also store the direction (incoming or outgoing) and the particle's four-momentum. +using QEDbase + +# Define a general stateful particle type that works for any particle species +struct ExampleParticleStateful{ + DIR<:ParticleDirection,SPECIES<:AbstractParticleType,ELEMENT<:AbstractFourMomentum +} <: AbstractParticleStateful{DIR,SPECIES,ELEMENT} + direction::DIR # Incoming or outgoing + species::SPECIES # Particle species (e.g., Muon, AntiMuon) + mom::ELEMENT # Particle's four-momentum +end + +# ## Implement the Interface Functions +# +# Next, we implement the required interface functions for `ExampleParticleStateful`, which are: +# +# - [`particle_direction`](@ref) +# - [`particle_species`](@ref) +# - [`momentum`](@ref) +# +# These functions extract the respective properties from the `ExampleParticleStateful` object. + +# Define the particle_direction function +function particle_direction(part::ExampleParticleStateful) + return part.direction +end + +# Define the particle_species function +function particle_species(part::ExampleParticleStateful) + return part.species +end + +# Define the momentum function +function momentum(part::ExampleParticleStateful) + return part.mom +end + +# ## Example Usage +# +# We can now create instances of `ExampleParticleStateful` for, say, `Electron` and `Positron`. +using QEDcore # to get predefined particles and four-momentum + +# Create a four-momentum vector (dummy example) +momentum_electron = SFourMomentum(1.0, 0.0, 0.0, 0.0); # E, px, py, pz + +# Create an incoming Electron using ExampleParticleStateful +incoming_electron = ExampleParticleStateful(Incoming(), Electron(), momentum_electron); + +# Create an outgoing Positron using ExampleParticleStateful +outgoing_positron = ExampleParticleStateful(Outgoing(), Positron(), momentum_electron); + +# Access particle properties +println("Incoming electron mass: ", mass(particle_species(incoming_electron))) +println( + "Is the outgoing positron a fermion? ", is_fermion(particle_species(outgoing_positron)) +) + +# Access momentum +println("Incoming electron momentum: ", momentum(incoming_electron)) +println("Outgoing positron momentum: ", momentum(outgoing_positron)) + +# ## Summary +# +# In this tutorial, we created a general `ExampleParticleStateful` type that can represent any particle species (like `Electron`,`Positron` from `QEDcore`, but also `Muon` or `AntiMuon` implemented in [this tutorial](@ref tutorial_particle)) by using the species as a type parameter. This approach avoids the need to define separate stateful types for each particle, making the implementation more flexible and reusable. +# +# The key steps were: +# +# 1. **Defining the general `ExampleParticleStateful` type** with parameters for direction, species, and momentum. +# 2. **Implementing the interface functions** to retrieve direction, species, and momentum. +# 3. **Using the generalized type** for all subtypes of `AbstractParticleType`. +# +# This method makes it easy to add more particles in the future while keeping the same structure. +# !!! note "Reference implementation" +# +# A reference implementation of a stateful particle type is given by +# [`ParticleStateful`](@extref :jl:type:`QEDcore.ParticleStateful`) diff --git a/docs/src/tutorial/phase_space_point.jl b/docs/src/tutorial/phase_space_point.jl new file mode 100644 index 0000000..ed27e72 --- /dev/null +++ b/docs/src/tutorial/phase_space_point.jl @@ -0,0 +1,159 @@ +# # Define a Custom Phase Space Point +# +# In this tutorial, we will define a custom **phase space point** type following the interface +# specification used in `QuantumElectrodynamics.jl`. +# +# The `AbstractPhaseSpacePoint` system has three key components: +# +# 1. **Scattering process:** The definition of the process, which includes the incoming and outgoing particles. +# 2. **Computational model:** The physics model, such as Quantum Electrodynamics (QED), that governs the dynamics of the particles. +# 3. **Phase space definition:** The structure used to create phase space points from coordinates in the multidimensional phase space. +# +# By the end of this tutorial, we’ll implement and test a custom phase space point using the **electron-positron annihilation** process: +# +# ```math +# e^+ + e^- \rightarrow \gamma + \gamma +# ``` + +# ## Step 0: Import the Necessary Ingredients +# We'll begin by importing necessary functionality from `QEDbase`, `QEDcore`, and `QEDprocesses`. + +using QEDbase +using QEDcore # Predefined particles and four-momenta +using QEDprocesses # Physics models and processes + +# ## Step 1: Define an Example Process +# We define a process that describes electron-positron annihilation. This process will involve +# two incoming particles (an electron and a positron) and two outgoing particles (two photons). + +struct ExampleProcess <: AbstractProcessDefinition end + +# Next, we specify the incoming and outgoing particles for the process by overloading the required interface functions: + +# ### Define Incoming and Outgoing Particles +# The incoming particles are an electron and a positron, and the outgoing particles are two photons. + +function QEDbase.incoming_particles(::ExampleProcess) + return (Electron(), Positron()) # These particle types are predefined in QEDcore +end + +function QEDbase.outgoing_particles(::ExampleProcess) + return (Photon(), Photon()) # Photons are also predefined in QEDcore +end + +# At this point, we have defined our `ExampleProcess` with the required particles. + +# ## Step 2: Implement the Model +# Next, we define a simple computational model, which will govern the behavior of the particles during the interaction. + +struct ExampleModel <: AbstractModelDefinition end + +# This model can later be extended to incorporate complex QED interactions, but for now, we use it as a placeholder. + +# ## Step 3: Define Phase Space Definition +# The last part we need to implement, is a phase space definition, which describes the +# coordinate structure of the phase space: + +struct ExamplePhaseSpaceDefinition <: AbstractPhasespaceDefinition end + +# This phase space definition can be extended later to include coordinate systems and +# frames of reference, but for now, this is also a placeholder. + +# ## Step 4: Define an Example Phase Space Point +# We now define the `ExamplePhaseSpacePoint` type, which will represent a specific point +# in the phase space of the electron-positron annihilation process. This type holds the process, +# the model, the phase space definition, and the incoming and outgoing particles. + +struct ExamplePhaseSpacePoint{PROC,MODEL,PSDEF,IN_PARTICLES,OUT_PARTICLES} <: + AbstractPhaseSpacePoint{PROC,MODEL,PSDEF,IN_PARTICLES,OUT_PARTICLES} + proc::PROC + mdl::MODEL + psdef::PSDEF + in_particles::IN_PARTICLES + out_particles::OUT_PARTICLES +end + +# ## Step 5: Implement the Interface Functions +# Now, we implement the interface functions required by `AbstractPhaseSpacePoint`. These include +# functions for retrieving the process, model, phase space definition, and particle information. + +function QEDbase.process(psp::ExamplePhaseSpacePoint) + return psp.proc +end + +function QEDbase.model(psp::ExamplePhaseSpacePoint) + return psp.mdl +end + +function QEDbase.phase_space_definition(psp::ExamplePhaseSpacePoint) + return psp.psdef +end + +function QEDbase.particles(psp::ExamplePhaseSpacePoint, dir::ParticleDirection) + if dir == Incoming() + return psp.in_particles + elseif dir == Outgoing() + return psp.out_particles + else + throw(ArgumentError("Invalid particle direction")) + end +end + +function Base.getindex(psp::ExamplePhaseSpacePoint, dir::ParticleDirection, n::Int) + return particles(psp, dir)[n] +end + +# With these functions, we can now retrieve the process, model, phase space definition, +# and individual particles from our `ExamplePhaseSpacePoint`. + +# ## Step 6: Access momenta +# The particles involved in the process have associated four-momenta, which describe their +# energy and momentum in spacetime. Let’s implement the functions to retrieve the momenta +# of particles. + +# Return the momentum of a specific particle +function QEDbase.momentum(psp::ExamplePhaseSpacePoint, dir::ParticleDirection, n::Int) + return momentum(particles(psp, dir)[n]) +end + +# Return a tuple of all momenta in a given direction +function QEDbase.momenta(psp::ExamplePhaseSpacePoint, dir::ParticleDirection) + return ntuple(i -> momentum(psp, dir, i), length(particles(psp, dir))) +end + +# ## Step 7: Test the Implementation +# Finally, we’ll test our implementation using the **electron-positron annihilation** process. +# +# 1. Create an instance of the process (`ExampleProcess`). +# 2. Create an instance of the model (`ExampleModel`). +# 3. Define the phase space configuration (`ExamplePhaseSpaceDefinition`). +# 4. Define the incoming and outgoing particles, specifying their four-momenta. +# 5. Create the phase space point and compute the momenta. + +proc = ExampleProcess() +model = ExampleModel() +psdef = ExamplePhaseSpaceDefinition() + +in_particles = ( + ParticleStateful(Incoming(), Electron(), SFourMomentum(1.0, 0.0, 0.0, 1.0)), + ParticleStateful(Incoming(), Positron(), SFourMomentum(1.0, 0.0, 0.0, -1.0)), +) + +out_particles = ( + ParticleStateful(Outgoing(), Photon(), SFourMomentum(1.0, 1.0, 0.0, 0.0)), + ParticleStateful(Outgoing(), Photon(), SFourMomentum(1.0, -1.0, 0.0, 0.0)), +) + +psp = ExamplePhaseSpacePoint(proc, model, psdef, in_particles, out_particles) + +incoming_momenta = momenta(psp, Incoming()) +outgoing_momenta = momenta(psp, Outgoing()) + +println("Incoming momenta: ", incoming_momenta) +println("Outgoing momenta: ", outgoing_momenta) + +# ## Conclusion +# +# We have successfully defined a custom phase space point type, implemented the necessary interface functions, +# and verified the momenta for an electron-positron annihilation process. +# This phase space point system provides a flexible way to work with scattering processes. diff --git a/docs/src/tutorial/process.jl b/docs/src/tutorial/process.jl new file mode 100644 index 0000000..ed6d00f --- /dev/null +++ b/docs/src/tutorial/process.jl @@ -0,0 +1,136 @@ +# # Tutorial: Defining a Custom Scattering Process +# +# In this tutorial, we'll define a custom scattering process by following the interface for [`AbstractProcessDefinition`](@ref). +# We'll cover the necessary functions, including particle types, spin and polarization handling, and matrix element calculations. + +# ## Step 1: Define a Custom Process Type +# +# Start by creating a custom type that inherits from `AbstractProcessDefinition`. +# This type will represent your process, for example, electron-positron scattering into two photons. + +using QEDbase +using QEDcore # for particle definitions and phase space points + +# Define a specific process by creating a subtype of `AbstractProcessDefinition`: +struct MyProcess <: AbstractProcessDefinition + ## Add relevant information, such as particle or any other process-specific properties. + incoming_particles::Tuple + outgoing_particles::Tuple +end + +# ## Step 2: Implement Required Particle Accessor Functions +# +# Every process must define the `incoming_particles` and `outgoing_particles` functions to list the particles involved. +# These functions return tuples of particles. + +QEDbase.incoming_particles(proc::MyProcess) = proc.incoming_particles +QEDbase.outgoing_particles(proc::MyProcess) = proc.outgoing_particles + +# ## Step 3: Handle Spins and Polarizations +# +# To define the spin and polarization for the particles, overload `incoming_spin_pols` and `outgoing_spin_pols`. +# For our custom process, we assume average/summation over all spins and polarizations. +# However, you can define specific spins or polarizations if needed. + +QEDbase.incoming_spin_pols(::MyProcess) = (AllSpin(), AllSpin()) # Both incoming particles are fermions (electron and positron) +QEDbase.outgoing_spin_pols(::MyProcess) = (AllPolarization(), AllPolarization()) # Photons are boson + +# ## Step 4: Define the Matrix Element Calculation +# +# The matrix element is a central part of any scattering process. It needs to be implemented for each spin and polarization combination. +# To calculate matrix elements, you must define the `_matrix_element` function. + +function QEDbase._matrix_element(psp::PhaseSpacePoint{MyProcess}) + ## Calculate the matrix element for the specific process. + ## This is a placeholder for the actual computation. + return 1.0 # Placeholder value +end + +# ## Step 5: Define Incident Flux +# +# To compute the cross section, we need to define the incident flux. +# This function calculates the initial flux of incoming particles, based on their momenta and energies. + +function QEDbase._incident_flux(psp::InPhaseSpacePoint{MyProcess}) + ## Placeholder calculation for incident flux + return 1.0 # Placeholder value +end + +# ## Step 6: Averaging Over Spin and Polarization +# +# Define the `_averaging_norm` function to return the normalization factor used to average the squared matrix elements over spins and polarizations. + +function QEDbase._averaging_norm(proc::MyProcess) + ## For example, if both incoming particles are fermions, the normalization could be the product of their spin multiplicity, i.e. 2 times 2. + return 4 # Placeholder value +end + +# ## Step 7: Check for Physical Phase Space +# +# The `_is_in_phasespace` function verifies whether a particular combination of incoming and outgoing momenta is physically allowed (on-shell, momentum conservation, etc.). + +function QEDbase._is_in_phasespace(psp::PhaseSpacePoint{MyProcess}) + ## Implement energy-momentum conservation and on-shell conditions. + return true # Placeholder value +end + +# ## Step 8: Define the Phase Space Factor +# +# To calculate cross sections, define the `_phase_space_factor` function, which returns the pre-differential factor of the invariant phase space integral measure. + +function QEDbase._phase_space_factor(psp::PhaseSpacePoint{MyProcess}) + ## Return the phase space factor + return 1.0 # Placeholder value +end + +# ## Step 9: Optional - Total Probability Calculation +# +# Optionally, you can define the `_total_probability` function to compute the total probability of the process. +# This is especially useful when computing total cross sections. + +function QEDbase._total_probability(psp::PhaseSpacePoint{MyProcess}) + ## Calculate the total probability for the process + return 1.0 # Placeholder value +end + +# ## Putting It All Together +# +# After defining the required functions, you now have a complete process definition for `MyProcess`. +# This process can be used in phase space integration, cross section calculation, and other scattering computations in `QuantumElectrodynamics.jl`. +# +# For example, if your process is electron-positron scattering into two photons: + +using QEDprocesses # for a predefined model: PerturbativeQED + +particles_in = (Electron(), Positron()) +particles_out = (Photon(), Photon()) +process = MyProcess(particles_in, particles_out) + +# You can then access the particles as follows: + +println("incoming particles: ", incoming_particles(process)) +println("outgoing particles: ", outgoing_particles(process)) + +# You can define some momenta for the incoming and outgoing particles + +electron_mass = mass(Electron()) +electron_momentum = SFourMomentum(3.0, 0, 0, sqrt(3.0^2 - electron_mass^2)) +positron_momentum = SFourMomentum(3.0, 0, 0, -sqrt(3.0^2 - electron_mass^2)) +incoming_momenta = (electron_momentum, positron_momentum) +photon_momentum1 = SFourMomentum(3.0, 3.0, 0, 0) +photon_momentum2 = SFourMomentum(3.0, -3.0, 0, 0) +outgoing_momenta = (photon_momentum1, photon_momentum2) + +# And you can define phase space definitions and computational models (here just as +# placeholder): + +ps_def = PhasespaceDefinition(SphericalCoordinateSystem(), CenterOfMomentumFrame()) +model = PerturbativeQED() + +# Finally, you can build a phase space point for your process + +psp = PhaseSpacePoint(process, model, ps_def, incoming_momenta, outgoing_momenta) + +# For this phase space point, the differential cross section can be calculated by calling + +println("differential cross section: ", differential_cross_section(psp)) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 98ec137..d458b48 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -26,7 +26,6 @@ export AbstractFourMomentum export isonshell, assert_onshell export AbstractDiracVector, AbstractDiracMatrix -export mul export AbstractGammaRepresentation @@ -46,7 +45,8 @@ export AbstractDefinitePolarization, AbstractIndefinitePolarization export PolarizationX, PolX, PolarizationY, PolY, AllPolarization, AllPol export AbstractDefiniteSpin, AbstractIndefiniteSpin export SpinUp, SpinDown, AllSpin -export multiplicity +export SyncedSpin, SyncedPolarization, SyncedPol +export spin_pols_iter # probabilities export differential_probability, unsafe_differential_probability @@ -63,6 +63,8 @@ export AbstractModelDefinition, fundamental_interaction_type export AbstractProcessDefinition, incoming_particles, outgoing_particles export number_incoming_particles, number_outgoing_particles export particles, number_particles +export incoming_spin_pols, outgoing_spin_pols, spin_pols +export multiplicity, incoming_multiplicity, outgoing_multiplicity # Abstract phase space definition interface export AbstractCoordinateSystem, AbstractFrameOfReference, AbstractPhasespaceDefinition @@ -73,6 +75,9 @@ export particle_direction, particle_species, momentum export process, model, phase_space_definition, momenta export AbstractInPhaseSpacePoint, AbstractOutPhaseSpacePoint +# Abstract coordinate transformation interface +export AbstractCoordinateTransformation + # errors export InvalidInputError, RegistryError, OnshellError, SpinorConstructionError @@ -102,17 +107,24 @@ include("interfaces/model.jl") include("interfaces/particle.jl") -include("particles/direction.jl") -include("particles/spin_pol.jl") +include("interfaces/particles/direction.jl") +include("interfaces/particles/spin_pol.jl") include("interfaces/phase_space.jl") include("interfaces/particle_stateful.jl") include("interfaces/process.jl") include("interfaces/phase_space_point.jl") -include("cross_section/diff_probability.jl") -include("cross_section/diff_cross_section.jl") -include("cross_section/total_probability.jl") -include("cross_section/total_cross_section.jl") +include("interfaces/coordinate_transformation.jl") + +include("implementations/process/momenta.jl") +include("implementations/process/particles.jl") +include("implementations/process/spin_pols.jl") +include("implementations/process/spin_pol_iterator.jl") + +include("implementations/cross_section/diff_probability.jl") +include("implementations/cross_section/diff_cross_section.jl") +include("implementations/cross_section/total_probability.jl") +include("implementations/cross_section/total_cross_section.jl") end #QEDbase diff --git a/src/errors.jl b/src/errors.jl index 7df0b51..b0337d9 100644 --- a/src/errors.jl +++ b/src/errors.jl @@ -33,8 +33,8 @@ function Base.showerror(io::IO, err::RegistryError) end """ -Abstract base type for exceptions indicating invalid input. See [`InvalidInputError`](@ref) for a simple concrete implementation. -Concrete implementations should at least implement +Abstract base type for exceptions indicating invalid input. See [`InvalidInputError`](@ref) for a simple concrete implementation. +Concrete implementations should at least implement ```Julia diff --git a/src/cross_section/diff_cross_section.jl b/src/implementations/cross_section/diff_cross_section.jl similarity index 88% rename from src/cross_section/diff_cross_section.jl rename to src/implementations/cross_section/diff_cross_section.jl index 652035c..7af6233 100644 --- a/src/cross_section/diff_cross_section.jl +++ b/src/implementations/cross_section/diff_cross_section.jl @@ -1,7 +1,7 @@ ######################## # differential and total cross sections. # -# This file contains default implementations for differential +# This file contains default implementations for differential # cross sections based on the scattering process interface ######################## @@ -23,8 +23,8 @@ If the given phase spaces are physical, return differential cross section evalua """ function differential_cross_section(phase_space_point::AbstractPhaseSpacePoint) if !_is_in_phasespace(phase_space_point) - # TODO: use the correct type here, i.e. implement a function `eltype` for psp or - # make `momentum_type` an interface function. + # TODO: use the correct type here, i.e. implement a function `eltype` for psp or + # make `momentum_type` an interface function. return zero(Float64) end diff --git a/src/cross_section/diff_probability.jl b/src/implementations/cross_section/diff_probability.jl similarity index 100% rename from src/cross_section/diff_probability.jl rename to src/implementations/cross_section/diff_probability.jl diff --git a/src/cross_section/total_cross_section.jl b/src/implementations/cross_section/total_cross_section.jl similarity index 100% rename from src/cross_section/total_cross_section.jl rename to src/implementations/cross_section/total_cross_section.jl diff --git a/src/cross_section/total_probability.jl b/src/implementations/cross_section/total_probability.jl similarity index 100% rename from src/cross_section/total_probability.jl rename to src/implementations/cross_section/total_probability.jl diff --git a/src/implementations/process/momenta.jl b/src/implementations/process/momenta.jl new file mode 100644 index 0000000..08bbdce --- /dev/null +++ b/src/implementations/process/momenta.jl @@ -0,0 +1,25 @@ +""" + _generate_momenta( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + in_phase_space::NTuple{N,T}, + out_phase_space::NTuple{M,T}, + ) where {N,M,T<:Real} + +Return four-momenta for incoming and outgoing particles for given coordinate based phase-space points. +""" +function _generate_momenta( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + in_phase_space::NTuple{N,T}, + out_phase_space::NTuple{M,T}, +) where {N,M,T<:Real} + in_momenta = _generate_incoming_momenta(proc, model, phase_space_def, in_phase_space) + out_momenta = _generate_outgoing_momenta( + proc, model, phase_space_def, in_phase_space, out_phase_space + ) + + return in_momenta, out_momenta +end diff --git a/src/implementations/process/particles.jl b/src/implementations/process/particles.jl new file mode 100644 index 0000000..00c415a --- /dev/null +++ b/src/implementations/process/particles.jl @@ -0,0 +1,59 @@ +""" + number_incoming_particles(proc_def::AbstractProcessDefinition) + +Return the number of incoming particles of a given process. +""" +@inline function number_incoming_particles(proc_def::AbstractProcessDefinition) + return length(incoming_particles(proc_def)) +end + +""" + number_outgoing_particles(proc_def::AbstractProcessDefinition) + +Return the number of outgoing particles of a given process. +""" +@inline function number_outgoing_particles(proc_def::AbstractProcessDefinition) + return length(outgoing_particles(proc_def)) +end + +""" + particles(proc_def::AbstractProcessDefinition, ::ParticleDirection) + +Convenience function dispatching to [`incoming_particles`](@ref) or [`outgoing_particles`](@ref) depending on the given direction. +""" +@inline particles(proc_def::AbstractProcessDefinition, ::Incoming) = + incoming_particles(proc_def) +@inline particles(proc_def::AbstractProcessDefinition, ::Outgoing) = + outgoing_particles(proc_def) + +""" + number_particles(proc_def::AbstractProcessDefinition, dir::ParticleDirection) + +Convenience function dispatching to [`number_incoming_particles`](@ref) or [`number_outgoing_particles`](@ref) depending on the given direction, returning the number of incoming or outgoing particles, respectively. +""" +@inline number_particles(proc_def::AbstractProcessDefinition, ::Incoming) = + number_incoming_particles(proc_def) +@inline number_particles(proc_def::AbstractProcessDefinition, ::Outgoing) = + number_outgoing_particles(proc_def) + +""" + number_particles(proc_def::AbstractProcessDefinition, dir::ParticleDirection, species::AbstractParticleType) + +Return the number of particles of the given direction and species in the given process definition. +""" +@inline function number_particles( + proc_def::AbstractProcessDefinition, dir::DIR, species::PT +) where {DIR<:ParticleDirection,PT<:AbstractParticleType} + return count(x -> x isa PT, particles(proc_def, dir)) +end + +""" + number_particles(proc_def::AbstractProcessDefinition, particle::AbstractParticleStateful) + +Return the number of particles of the given particle's direction and species in the given process definition. +""" +@inline function number_particles( + proc_def::AbstractProcessDefinition, ps::AbstractParticleStateful +) + return number_particles(proc_def, particle_direction(ps), particle_species(ps)) +end diff --git a/src/implementations/process/spin_pol_iterator.jl b/src/implementations/process/spin_pol_iterator.jl new file mode 100644 index 0000000..ad1c2da --- /dev/null +++ b/src/implementations/process/spin_pol_iterator.jl @@ -0,0 +1,132 @@ +""" + SpinPolIter + +An iterator type to iterate over spin and polarization combinations. Should be used through [`spin_pols_iter`](@ref). +""" +struct SpinPolIter{I,O} + # product iterator doing the actual iterating + product_iter::Base.Iterators.ProductIterator + # lookup table for which indices go where, translating the base iterator to the actual result + indexing_lut::Tuple{NTuple{I,Int},NTuple{O,Int}} +end + +""" + all_spin_pols(process::AbstractProcessDefinition) + +This function returns an iterator, yielding every fully definite combination of spins and polarizations allowed by the +process' [`spin_pols`](@ref). Each returned element is a Tuple of the incoming and the outgoing spins and polarizations, +in the order of the process' own spins and polarizations. + +This works together with the definite spins and polarizations, [`AllSpin`](@ref), [`AllPolarization`](@ref), and the synced versions +[`SyncedPolarization`](@ref) and [`SyncedSpin`](@ref). + +```jldoctest +julia> using QEDbase; using QEDcore; using QEDprocesses; + +julia> proc = ScatteringProcess((Photon(), Photon(), Photon(), Electron()), (Photon(), Electron()), (SyncedPolarization(1), SyncedPolarization(2), SyncedPolarization(1), SpinUp()), (SyncedPolarization(2), AllSpin())) +generic QED process + incoming: photon (synced polarization 1), photon (synced polarization 2), photon (synced polarization 1), electron (spin up) + outgoing: photon (synced polarization 2), electron (all spins) + + +julia> for sp_combo in spin_pols_iter(proc) println(sp_combo) end +((x-polarized, x-polarized, x-polarized, spin up), (x-polarized, spin up)) +((y-polarized, x-polarized, y-polarized, spin up), (x-polarized, spin up)) +((x-polarized, y-polarized, x-polarized, spin up), (y-polarized, spin up)) +((y-polarized, y-polarized, y-polarized, spin up), (y-polarized, spin up)) +((x-polarized, x-polarized, x-polarized, spin up), (x-polarized, spin down)) +((y-polarized, x-polarized, y-polarized, spin up), (x-polarized, spin down)) +((x-polarized, y-polarized, x-polarized, spin up), (y-polarized, spin down)) +((y-polarized, y-polarized, y-polarized, spin up), (y-polarized, spin down)) + +julia> length(spin_pols_iter(proc)) +8 +``` +""" +function spin_pols_iter(process::AbstractProcessDefinition) + DEF_SPINS = (SpinUp(), SpinDown()) + DEF_POLS = (PolX(), PolY()) + + in_sp = incoming_spin_pols(process) + I = length(in_sp) + out_sp = outgoing_spin_pols(process) + O = length(out_sp) + + # concatenate for now for easier indices, split again later + sps = (in_sp..., out_sp...) + + # keep indices of first seen SyncedSpins or SyncedPols + synced_seen = Dict{AbstractSpinOrPolarization,Int}() + index = 0 + for sp in sps + index += 1 + if !(sp isa SyncedSpin || sp isa SyncedPolarization) + continue + end + if !haskey(synced_seen, sp) + synced_seen[sp] = index + end + end + + # keep indices of the synced spins/pols in the iterator (not necessarily the same as synced_seen) + synced_indices = Dict{AbstractSpinOrPolarization,Int}() + + iter_tuples = Vector() + lut = Vector{Int}() + index = 0 + for sp in sps + index += 1 + if sp isa AbstractDefiniteSpin || sp isa AbstractDefinitePolarization + push!(iter_tuples, (sp,)) + push!(lut, length(iter_tuples)) + elseif sp isa SyncedSpin + # check if it's the first synced + if index == synced_seen[sp] + push!(iter_tuples, DEF_SPINS) + synced_indices[sp] = length(iter_tuples) + end + push!(lut, synced_indices[sp]) + elseif sp isa SyncedPolarization + if index == synced_seen[sp] + push!(iter_tuples, DEF_POLS) + synced_indices[sp] = length(iter_tuples) + end + push!(lut, synced_indices[sp]) + elseif sp isa AllSpin + push!(iter_tuples, DEF_SPINS) + push!(lut, length(iter_tuples)) + elseif sp isa AllPol + push!(iter_tuples, DEF_POLS) + push!(lut, length(iter_tuples)) + end + end + + return SpinPolIter( + Iterators.product(iter_tuples...), + (tuple(lut[begin:I]...), tuple(lut[(I + 1):end]...)), + ) +end + +function Base.iterate(iterator::SpinPolIter, state=nothing) + local prod_iter_res + if isnothing(state) + prod_iter_res = iterate(iterator.product_iter) + else + prod_iter_res = iterate(iterator.product_iter, state) + end + + if isnothing(prod_iter_res) + return nothing + end + prod_iter_res, state = prod_iter_res + + # translate prod_iter_res into actual result + in_t = ((prod_iter_res[i] for i in iterator.indexing_lut[1])...,) + out_t = ((prod_iter_res[i] for i in iterator.indexing_lut[2])...,) + + return (in_t, out_t), state +end + +function Base.length(iterator::SpinPolIter) + return length(iterator.product_iter) +end diff --git a/src/implementations/process/spin_pols.jl b/src/implementations/process/spin_pols.jl new file mode 100644 index 0000000..a4e50cd --- /dev/null +++ b/src/implementations/process/spin_pols.jl @@ -0,0 +1,144 @@ + +# default implementation +function incoming_spin_pols(proc_def::AbstractProcessDefinition) + return ntuple( + x -> is_fermion(incoming_particles(proc_def)[x]) ? AllSpin() : AllPolarization(), + number_incoming_particles(proc_def), + ) +end + +# default implementation +function outgoing_spin_pols(proc_def::AbstractProcessDefinition) + return ntuple( + x -> is_fermion(outgoing_particles(proc_def)[x]) ? AllSpin() : AllPolarization(), + number_outgoing_particles(proc_def), + ) +end + +""" + spin_pols(proc_def::AbstractProcessDefinition, dir::ParticleDirection) + +Return the tuple of spins and polarizations for the process in the given direction. Dispatches to [`incoming_spin_pols`](@ref) or [`outgoing_spin_pols`](@ref). +""" +spin_pols(proc_def::AbstractProcessDefinition, dir::Incoming) = incoming_spin_pols(proc_def) +spin_pols(proc_def::AbstractProcessDefinition, dir::Outgoing) = outgoing_spin_pols(proc_def) + +@inline _multiplicity(::AbstractDefinitePolarization) = 1 +@inline _multiplicity(::AbstractDefiniteSpin) = 1 +@inline _multiplicity(::AllSpin) = 2 +@inline _multiplicity(::AllPol) = 2 + +@inline _multiplicity(::SyncedPolarization{N}, seen_spin_pols::Tuple{}) where {N} = 2 +@inline _multiplicity(::SyncedSpin{N}, seen_spin_pols::Tuple{}) where {N} = 2 + +@inline function _multiplicity( + sp::SyncedPolarization{N}, seen_pols::Tuple{Int,Vararg{Int}} +) where {N} + if seen_pols[1] == N + return 1 + else + return _multiplicity(sp, seen_pols[2:end]) + end +end + +@inline function _multiplicity( + ss::SyncedSpin{N}, seen_spins::Tuple{Int,Vararg{Int}} +) where {N} + if seen_spins[1] == N + return 1 + else + return _multiplicity(ss, seen_spins[2:end]) + end +end + +# multiplicity of the empty spin_pol set is 1 +@inline _multiplicity(::Tuple{}, _, _) = 1 + +# recursion for abstract spins or pols that are not synced +@inline function _multiplicity( + spin_pols::Tuple{AbstractSpin,Vararg{AbstractSpinOrPolarization}}, + seen_spins::NTuple{S,Int}, + seen_pols::NTuple{P,Int}, +) where {S,P} + return _multiplicity(spin_pols[1]) * + _multiplicity(spin_pols[2:end], seen_spins, seen_pols) +end +@inline function _multiplicity( + spin_pols::Tuple{AbstractPolarization,Vararg{AbstractSpinOrPolarization}}, + seen_spins::NTuple{S,Int}, + seen_pols::NTuple{P,Int}, +) where {S,P} + return _multiplicity(spin_pols[1]) * + _multiplicity(spin_pols[2:end], seen_spins, seen_pols) +end + +# recursion for synced spins or pols +@inline function _multiplicity( + spin_pols::Tuple{SyncedPolarization{N},Vararg{AbstractSpinOrPolarization}}, + seen_spins::NTuple{S,Int}, + seen_pols::NTuple{P,Int}, +) where {N,S,P} + return _multiplicity(spin_pols[1], seen_pols) * + _multiplicity(spin_pols[2:end], seen_spins, (seen_pols..., N)) +end +@inline function _multiplicity( + spin_pols::Tuple{SyncedSpin{N},Vararg{AbstractSpinOrPolarization}}, + seen_spins::NTuple{S,Int}, + seen_pols::NTuple{P,Int}, +) where {N,S,P} + return _multiplicity(spin_pols[1], seen_spins) * + _multiplicity(spin_pols[2:end], (seen_spins..., N), seen_pols) +end + +""" + multiplicity(proc::AbstractProcessDefinition) + +Return the number of spin and polarization combinations represented by `proc` total. This depends on the specific [`AbstractSpinOrPolarization`](@ref)s returned by [`spin_pols`](@ref) for `proc`. +For example, a default Compton process with four indefinite spins/polarizations has a multiplicity of 2^4 = 16. A Compton process with many incoming photons that have synced polarizations +will still have a multiplicity of 16. + +!!! note Performance + As long as [`incoming_spin_pols`](@ref) and [`outgoing_spin_pols`](@ref) can be evaluated at compile time, this function is completely compiled away. + +!!! note Incoming and Outgoing Spins/Polarizations + Note that the total multiplicity is not necessarily the incoming and outgoing multiplicities multiplied. This is the case when incoming and outgoing particles are synced with one another. + +See also: [`SyncedPolarization`](@ref), [`SyncedSpin`](@ref), [`incoming_multiplicity`](@ref), [`outgoing_multiplicity`](@ref) +""" +function multiplicity(proc::AbstractProcessDefinition) + return _multiplicity( + (incoming_spin_pols(proc)..., outgoing_spin_pols(proc)...), + NTuple{0,Int}(), + NTuple{0,Int}(), + ) +end + +""" + incoming_multiplicity(proc::AbstractProcessDefinition) + +Return the number of spin and polarization combinations represented by `proc`s incoming particles. This function only considers the incoming particles' spins and polarizations, returned by +[`incoming_spin_pols`](@ref) for `proc`. + +!!! note Incoming and Outgoing Spins/Polarizations + Note that the total multiplicity is not necessarily the incoming and outgoing multiplicities multiplied. For the total process multiplicity, see [`multiplicity`](@ref). + +See also: [`SyncedPolarization`](@ref), [`SyncedSpin`](@ref), [`multiplicity`](@ref), [`outgoing_multiplicity`](@ref) +""" +function incoming_multiplicity(proc::AbstractProcessDefinition) + return _multiplicity(incoming_spin_pols(proc), NTuple{0,Int}(), NTuple{0,Int}()) +end + +""" + outgoing_multiplicity(proc::AbstractProcessDefinition) + +Return the number of spin and polarization combinations represented by `proc`s outgoing particles. This function only considers the outgoing particles' spins and polarizations, returned by +[`outgoing_spin_pols`](@ref) for `proc`. + +!!! note Incoming and Outgoing Spins/Polarizations + Note that the total multiplicity is not necessarily the incoming and outgoing multiplicities multiplied. For the total process multiplicity, see [`multiplicity`](@ref). + +See also: [`SyncedPolarization`](@ref), [`SyncedSpin`](@ref), [`multiplicity`](@ref), [`incoming_multiplicity`](@ref) +""" +function outgoing_multiplicity(proc::AbstractProcessDefinition) + return _multiplicity(outgoing_spin_pols(proc), NTuple{0,Int}(), NTuple{0,Int}()) +end diff --git a/src/interfaces/coordinate_transformation.jl b/src/interfaces/coordinate_transformation.jl new file mode 100644 index 0000000..c157a6f --- /dev/null +++ b/src/interfaces/coordinate_transformation.jl @@ -0,0 +1,153 @@ +####### +# General coordinate transformations +####### + +""" + AbstractCoordinateTransformation + +Abstract base type for coordinate transformations supposed to be acting on four-momenta. +Every subtype of `trafo::AbstractCoordinateTransformation` should implement the following interface functions: + +* `QEDbase._transform(trafo,p)`: transforms `p` +* `Base.inv(trafo)`: returns the inverted transform + +## Example + +Implementing the interface by defining the interface functions: + +```jldoctest trafo_interface +julia> using QEDbase + +julia> using QEDcore + +julia> struct TestTrafo{T} <: AbstractCoordinateTransformation + a::T + end + +julia> QEDbase._transform(trafo::TestTrafo,p) = trafo.a*p + +julia> Base.inv(trafo::TestTrafo) = TestTrafo(inv(trafo.a)) + +``` + +The `TestTrafo` can then be used to transform four-momenta: + +```jldoctest trafo_interface +julia> trafo = TestTrafo(2.0) +TestTrafo{Float64}(2.0) + +julia> p = SFourMomentum(4,3,2,1) +4-element SFourMomentum with indices SOneTo(4): + 4.0 + 3.0 + 2.0 + 1.0 + +julia> trafo(p) # multiply every component with 2.0 +4-element SFourMomentum with indices SOneTo(4): + 8.0 + 6.0 + 4.0 + 2.0 + +julia> inv(trafo)(p) # divide every component by 2.0 +4-element SFourMomentum with indices SOneTo(4): + 2.0 + 1.5 + 1.0 + 0.5 +``` +""" +abstract type AbstractCoordinateTransformation end +Base.broadcastable(trafo::AbstractCoordinateTransformation) = Ref(trafo) + +""" + _transform(trafo::AbstractCoordinateTransformation,p::AbstractFourMomentum) + +Interface function for the application of the transformation to the four-momentum `p`. +Must return a four-momentum of the same type as `p`. +""" +function _transform end + +# make the transform callable +@inline function (trafo::AbstractCoordinateTransformation)(p::AbstractFourMomentum) + return _transform(trafo, p) +end + +@inline function (trafo::AbstractCoordinateTransformation)( + psf::PSF +) where {PSF<:AbstractParticleStateful} + p_prime = _transform(trafo, momentum(psf)) + return PSF(p_prime) +end + +@inline function (trafo::AbstractCoordinateTransformation)( + psp::PSP +) where {PSP<:AbstractPhaseSpacePoint} + in_moms = momenta(psp, Incoming()) + out_moms = momenta(psp, Outgoing()) + in_moms_prime = _transform.(trafo, in_moms) + out_moms_prime = _transform.(trafo, out_moms) + + proc = process(psp) + mod = model(psp) + ps_def = phase_space_definition(psp) + return constructorof(PSP)(proc, mod, ps_def, in_moms_prime, out_moms_prime) +end + +######### +# Abstract Lorentz Boosts +######### + +""" + + AbstractLorentzTransformation <: AbstractCoordinateTransformation + +An abstract base type representing Lorentz transformations, which are coordinate +transformations between inertial and reference frames in special relativity. + +`AbstractLorentzTransformation` extends `AbstractCoordinateTransformation` and provides +the foundational framework for all types of Lorentz transformations, including boosts. +These transformations preserve the Minkowski product of two four-vectors and are fundamental to +the description of relativistic physics, ensuring the laws of physics are the same in all +inertial frames. + +""" +abstract type AbstractLorentzTransformation <: AbstractCoordinateTransformation end + +""" + + AbstractLorentzBoost <: AbstractLorentzTransformation + +An abstract base type representing Lorentz boosts, a specific type of Lorentz transformation +associated with relative motion between inertial frames along one or more spatial directions. + +`AbstractLorentzBoost` extends `AbstractLorentzTransformation` and serves as the foundation +for all types of boost transformations in special relativity. Lorentz boosts describe how +four-vectors change when transitioning between two reference frames moving at constant velocities (in units of the speed of light) relative to each other. +""" +abstract type AbstractLorentzBoost <: AbstractLorentzTransformation end + +""" + + AbstractBoostParameter + +An abstract base type representing boost parameters used in Lorentz transformations, which +describe the relative motion between two inertial frames in special relativity. + +`AbstractBoostParameter` serves as the foundation for defining specific boost parameters +that control Lorentz boosts in different spatial directions. Boost parameters typically +represent the velocity of one reference frame relative to another, expressed as a fraction +of the speed of light (`\\beta`), and are essential for performing Lorentz transformations +on four-vectors. + +## Overview + +In the context of special relativity, a Lorentz boost is a transformation that changes the +time and spatial components of a four-vector based on the relative motion between two +inertial reference frames. For example, the boost parameter ``\\beta`` is dimensionless and represents +this velocity as a fraction of the speed of light. Depending on the frame's relative velocity, +different forms of boost parameters exist, such as those associated with a single axis or +a vector describing boosts in multiple spatial dimensions. +""" +abstract type AbstractBoostParameter end diff --git a/src/interfaces/four_momentum.jl b/src/interfaces/four_momentum.jl index f79de00..449c684 100644 --- a/src/interfaces/four_momentum.jl +++ b/src/interfaces/four_momentum.jl @@ -3,7 +3,7 @@ Abstract base type for four-momentas, representing one energy and three spacial components. -Also see: `QEDcore.SFourMomentum`, `QEDcore.MFourMomentum` +Also see: [`QEDcore.SFourMomentum`](@extref), [`QEDcore.MFourMomentum`](@extref) """ abstract type AbstractFourMomentum <: AbstractLorentzVector{Float64} end diff --git a/src/interfaces/lorentz_vectors/arithmetic.jl b/src/interfaces/lorentz_vectors/arithmetic.jl index 6f6e056..6282150 100644 --- a/src/interfaces/lorentz_vectors/arithmetic.jl +++ b/src/interfaces/lorentz_vectors/arithmetic.jl @@ -12,11 +12,11 @@ end """ minkowski_dot(v1,v2) -Return the Minkowski dot product of two `LorentzVectorLike`. +Return the Minkowski dot product of two `LorentzVectorLike`. !!! example - If `(t1,x1,y1,z1)` and `(t2,x2,y2,z2)` are two `LorentzVectorLike`, this is equivalent to + If `(t1,x1,y1,z1)` and `(t2,x2,y2,z2)` are two `LorentzVectorLike`, this is equivalent to ```julia t1*t2 - (x1*x2 + y1*y2 + z1*z2) ``` @@ -41,9 +41,9 @@ const mdot = minkowski_dot """ getMagnitude2(lv) -Return the square of the magnitude of a given `LorentzVectorLike`, i.e. the sum of the squared spatial components. +Return the square of the magnitude of a given `LorentzVectorLike`, i.e. the sum of the squared spatial components. -!!! example +!!! example If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `x^2+ y^2 + z^2`. @@ -65,9 +65,9 @@ const getMag2 = getMagnitude2 """ getMagnitude(lv) -Return the magnitude of a given `LorentzVectorLike`, i.e. the euklidian norm spatial components. +Return the magnitude of a given `LorentzVectorLike`, i.e. the euklidian norm spatial components. -!!! example +!!! example If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `sqrt(x^2 + y^2 + z^2)`. @@ -89,11 +89,11 @@ const getMag = getMagnitude """ getInvariantMass2(lv) -Return the squared invariant mass of a given `LorentzVectorLike`, i.e. the minkowski dot with itself. +Return the squared invariant mass of a given `LorentzVectorLike`, i.e. the minkowski dot with itself. -!!! example +!!! example - If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `t^2 - (x^2 + y^2 + z^2)`. + If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `t^2 - (x^2 + y^2 + z^2)`. """ @@ -107,16 +107,16 @@ const getMass2 = getInvariantMass2 """ getInvariantMass(lv) -Return the the invariant mass of a given `LorentzVectorLike`, i.e. the square root of the minkowski dot with itself. +Return the the invariant mass of a given `LorentzVectorLike`, i.e. the square root of the minkowski dot with itself. -!!! example +!!! example If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `sqrt(t^2 - (x^2 + y^2 + z^2))`. !!! note - - If the squared invariant mass `m2` is negative, `-sqrt(-m2)` is returned. + + If the squared invariant mass `m2` is negative, `-sqrt(-m2)` is returned. """ @traitfn function getInvariantMass(lv::T) where {T; IsLorentzVectorLike{T}} @@ -140,9 +140,9 @@ const getMass = getInvariantMass """ getE(lv) -Return the energy component of a given `LorentzVectorLike`, i.e. its 0-component. +Return the energy component of a given `LorentzVectorLike`, i.e. its 0-component. -!!! example +!!! example If `(E,px,py,pz)` is a `LorentzVectorLike`, this is equivalent to `E`. @@ -155,9 +155,9 @@ const getEnergy = getE """ getPx(lv) -Return the ``p_x`` component of a given `LorentzVectorLike`, i.e. its 1-component. +Return the ``p_x`` component of a given `LorentzVectorLike`, i.e. its 1-component. -!!! example +!!! example If `(E,px,py,pz)` is a `LorentzVectorLike`, this is equivalent to `px`. @@ -167,9 +167,9 @@ Return the ``p_x`` component of a given `LorentzVectorLike`, i.e. its 1-componen """ getPy(lv) -Return the ``p_y`` component of a given `LorentzVectorLike`, i.e. its 2-component. +Return the ``p_y`` component of a given `LorentzVectorLike`, i.e. its 2-component. -!!! example +!!! example If `(E,px,py,pz)` is a `LorentzVectorLike`, this is equivalent to `py`. @@ -179,9 +179,9 @@ Return the ``p_y`` component of a given `LorentzVectorLike`, i.e. its 2-componen """ getPz(lv) -Return the ``p_z`` component of a given `LorentzVectorLike`, i.e. its 3-component. +Return the ``p_z`` component of a given `LorentzVectorLike`, i.e. its 3-component. -!!! example +!!! example If `(E,px,py,pz)` is a `LorentzVectorLike`, this is equivalent to `pz`. @@ -240,7 +240,7 @@ Return the squared transverse momentum for a given `LorentzVectorLike`, i.e. the !!! note - The transverse components are defined w.r.t. to the 3-axis. + The transverse components are defined w.r.t. to the 3-axis. """ @inline @traitfn function getTransverseMomentum2(lv::T) where {T; IsLorentzVectorLike{T}} @@ -264,7 +264,7 @@ Return the transverse momentum for a given `LorentzVectorLike`, i.e. the magnitu !!! note - The transverse components are defined w.r.t. to the 3-axis. + The transverse components are defined w.r.t. to the 3-axis. """ @@ -288,7 +288,7 @@ Return the squared transverse mass for a given `LorentzVectorLike`, i.e. the dif !!! note - The transverse components are defined w.r.t. to the 3-axis. + The transverse components are defined w.r.t. to the 3-axis. """ @inline @traitfn function getTransverseMass2(lv::T) where {T; IsLorentzVectorLike{T}} @@ -309,7 +309,7 @@ Return the transverse momentum for a given `LorentzVectorLike`, i.e. the square !!! note - The transverse components are defined w.r.t. to the 3-axis. + The transverse components are defined w.r.t. to the 3-axis. !!! note @@ -340,7 +340,7 @@ Return the [rapidity](https://en.wikipedia.org/wiki/Rapidity) for a given `Loren !!! note - The transverse components are defined w.r.t. to the 3-axis. + The transverse components are defined w.r.t. to the 3-axis. """ @inline @traitfn function getRapidity(lv::T) where {T; IsLorentzVectorLike{T}} @@ -367,7 +367,7 @@ Return the theta angle for a given `LorentzVectorLike`, i.e. the polar angle of !!! note - The [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system) are defined w.r.t. to the 3-axis. + The [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system) are defined w.r.t. to the 3-axis. """ @inline @traitfn function getTheta(lv::T) where {T; IsLorentzVectorLike{T}} @@ -385,7 +385,7 @@ end Return the cosine of the theta angle for a given `LorentzVectorLike`. -!!! note +!!! note This is an equivalent but faster version of `cos(getTheta(lv))`; see [`getTheta`](@ref). @@ -406,7 +406,7 @@ Return the phi angle for a given `LorentzVectorLike`, i.e. the azimuthal angle o !!! note - The [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system) are defined w.r.t. to the 3-axis. + The [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system) are defined w.r.t. to the 3-axis. """ @traitfn function getPhi(lv::T) where {T; IsLorentzVectorLike{T}} @@ -419,7 +419,7 @@ end Return the cosine of the phi angle for a given `LorentzVectorLike`. -!!! note +!!! note This is an equivalent but faster version of `cos(getPhi(lv))`; see [`getPhi`](@ref). @@ -434,7 +434,7 @@ end Return the sine of the phi angle for a given `LorentzVectorLike`. -!!! note +!!! note This is an equivalent but faster version of `sin(getPhi(lv))`; see [`getPhi`](@ref). @@ -459,7 +459,7 @@ Return the plus component for a given `LorentzVectorLike` in [light-cone coordin !!! note The [light-cone coordinates](https://en.wikipedia.org/wiki/Light-cone_coordinates) are defined w.r.t. to the 3-axis. - + !!! warning The definition ``p^+ := (E + p_z)/2` differs from the usual definition of [light-cone coordinates](https://en.wikipedia.org/wiki/Light-cone_coordinates) in general relativity. @@ -481,7 +481,7 @@ Return the minus component for a given `LorentzVectorLike` in [light-cone coordi !!! note The [light-cone coordinates](https://en.wikipedia.org/wiki/Light-cone_coordinates) are defined w.r.t. to the 3-axis. - + !!! warning The definition ``p^- := (E - p_z)/2` differs from the usual definition of [light-cone coordinates](https://en.wikipedia.org/wiki/Light-cone_coordinates) in general relativity. diff --git a/src/interfaces/lorentz_vectors/dirac_interaction.jl b/src/interfaces/lorentz_vectors/dirac_interaction.jl index 990d949..9a83e18 100644 --- a/src/interfaces/lorentz_vectors/dirac_interaction.jl +++ b/src/interfaces/lorentz_vectors/dirac_interaction.jl @@ -8,7 +8,7 @@ Product of generic Lorentz vector with a Dirac tensor from the left. Basically, !!! note "Multiplication operator" This also overloads the `*` operator for this types. """ -function mul( +function _mul( DM::T, L::TL ) where {T<:Union{AbstractDiracMatrix,AbstractDiracVector},TL<:AbstractLorentzVector} return constructorof(TL)(DM * L[1], DM * L[2], DM * L[3], DM * L[4]) @@ -16,7 +16,7 @@ end @inline function *( DM::T, L::TL ) where {T<:Union{AbstractDiracMatrix,AbstractDiracVector},TL<:AbstractLorentzVector} - return mul(DM, L) + return _mul(DM, L) end """ @@ -28,7 +28,7 @@ Product of generic Lorentz vector with a Dirac tensor from the right. Basically, This also overloads the `*` operator for this types. """ -function mul( +function _mul( L::TL, DM::T ) where {TL<:AbstractLorentzVector,T<:Union{AbstractDiracMatrix,AbstractDiracVector}} return constructorof(TL)(L[1] * DM, L[2] * DM, L[3] * DM, L[4] * DM) @@ -36,5 +36,5 @@ end @inline function *( L::TL, DM::T ) where {TL<:AbstractLorentzVector,T<:Union{AbstractDiracMatrix,AbstractDiracVector}} - return mul(L, DM) + return _mul(L, DM) end diff --git a/src/interfaces/lorentz_vectors/fields.jl b/src/interfaces/lorentz_vectors/fields.jl index 1a7fc9d..14fdba4 100644 --- a/src/interfaces/lorentz_vectors/fields.jl +++ b/src/interfaces/lorentz_vectors/fields.jl @@ -3,7 +3,7 @@ Return the 0-component of a given `LorentzVectorLike`. -!!! example +!!! example If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `t`. @@ -15,7 +15,7 @@ function getT end Return the 1-component of a given `LorentzVectorLike`. -!!! example +!!! example If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `x`. @@ -27,7 +27,7 @@ function getX end Return the 2-component of a given `LorentzVectorLike`. -!!! example +!!! example If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `y`. @@ -39,7 +39,7 @@ function getY end Return the 3-component of a given `LorentzVectorLike`. -!!! example +!!! example If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `z`. diff --git a/src/interfaces/lorentz_vectors/registry.jl b/src/interfaces/lorentz_vectors/registry.jl index 643e7da..0e5ac18 100644 --- a/src/interfaces/lorentz_vectors/registry.jl +++ b/src/interfaces/lorentz_vectors/registry.jl @@ -13,10 +13,10 @@ end """ $(SIGNATURES) -Function to register a custom type as a LorentzVectorLike. +Function to register a custom type as a LorentzVectorLike. -Ensure the passed custom type has implemented at least the function `getT, getX, getY, getZ` -and enables getter functions of the lorentz vector library for the given type. +Ensure the passed custom type has implemented at least the function `getT, getX, getY, getZ` +and enables getter functions of the lorentz vector library for the given type. If additionally the functions `setT!, setX!, setY!, setZ!` are implemened for the passed custom type, also the setter functions of the Lorentz vector interface are enabled. """ diff --git a/src/interfaces/lorentz_vectors/types.jl b/src/interfaces/lorentz_vectors/types.jl index 456aa37..cab00e6 100644 --- a/src/interfaces/lorentz_vectors/types.jl +++ b/src/interfaces/lorentz_vectors/types.jl @@ -9,41 +9,3 @@ $(TYPEDEF) Abstract type to model generic Lorentz vectors, i.e. elements of a minkowski-like space, where the component space is arbitray. """ abstract type AbstractLorentzVector{T} <: FieldVector{4,T} end - -""" -## Definition of LorentzVector interface. - -It enables several functions related to Lorentz vectors for any custom type just by implementing the API for the respective type. -For instance, say we want to implement a type, which shall act like a Lorentz vector. Then, we start with our custom type: - -```julia -struct CustomType - t::Float64 - x::Float64 - y::Float64 - z::Float64 -end -``` -The first group of functions to be implemented for this `CustomType` in order to connect this type to the Lorentz vector interface are the getter for the cartesian components. - -```julia -QEDbase.getT(lv::CustomType) = lv.t -QEDbase.getX(lv::CustomType) = lv.x -QEDbase.getY(lv::CustomType) = lv.y -QEDbase.getZ(lv::CustomType) = lv.z -``` -Make sure that you dispatch on the getter from `QEDbase` by defining `QEDbase.getT`, etc. - -With this done, we can aleady register the `CustomType` as a `LorentzVectorLike` type using the `register_LorentzVectorLike` function -```julia -register_LorentzVectorLike(CustomType) -``` -If something goes wrong, this function call will raise an `RegistryError` indicating, what is missing. With this done, `CustomType` is ready to be used as a `LorentzVectorLike` -```julia -L = CustomType(2.0,1.0,0.0,-1.0) - -getTheta(L) # -getRapidity(L) # -``` - -""" diff --git a/src/interfaces/lorentz_vectors/utility.jl b/src/interfaces/lorentz_vectors/utility.jl index 46dd554..30dda99 100644 --- a/src/interfaces/lorentz_vectors/utility.jl +++ b/src/interfaces/lorentz_vectors/utility.jl @@ -13,11 +13,11 @@ end """ $(SIGNATURES) -On-shell check of a given four-momentum `mom` w.r.t. a given mass `mass`. +On-shell check of a given four-momentum `mom` w.r.t. a given mass `mass`. !!! note "Precision" For `AbstactFourMomentum`, the element type is fixed to `Float64`, limiting the precision of comparisons between elements. - The current implementation has been tested within the boundaries for energy scales `E` with `1e-9 <= E <= 1e5`. + The current implementation has been tested within the boundaries for energy scales `E` with `1e-9 <= E <= 1e5`. In those bounds, the mass error, which is correctly detected as off-shell, is `1e-4` times the mean value of the components, but has at most the value `0.01`, e.g. at the high energy end. The energy scales correspond to `0.5 meV` for the lower bound and `50 GeV` for the upper bound. @@ -41,7 +41,7 @@ Assertion if a FourMomentum `mom` is on-shell w.r.t a given mass `mass`. !!! note "See also" The precision of this functions is explained in [`isonshell`](@ref). - + """ function assert_onshell(mom::QEDbase.AbstractLorentzVector, mass::Real) isonshell(mom, mass) || throw(OnshellError(mom, mass)) diff --git a/src/interfaces/model.jl b/src/interfaces/model.jl index 7b98325..fce7606 100644 --- a/src/interfaces/model.jl +++ b/src/interfaces/model.jl @@ -6,7 +6,7 @@ ############### # root type for models """ -Abstract base type for all compute model definitions in the context of scattering processes. Every subtype of `AbstractModelDefinition` is associated with a fundamental interaction. +Abstract base type for all compute model definitions in the context of scattering processes. Every subtype of `AbstractModelDefinition` is associated with a fundamental interaction. Therefore, one needs to implement the following soft interface function ```Julia diff --git a/src/interfaces/particle.jl b/src/interfaces/particle.jl index b4f0e36..7abed2b 100644 --- a/src/interfaces/particle.jl +++ b/src/interfaces/particle.jl @@ -2,11 +2,11 @@ # The particle interface # # In this file, we define the interface for working with particles in a general -# sense. +# sense. ############### """ -Abstract base type for every type which might be considered a *particle* in the context of `QED.jl`. For every (concrete) subtype of `AbstractParticle`, there are two kinds of interface functions implemented: static functions and property functions. +Abstract base type for every type which might be considered a *particle* in the context of `QuantumElectrodynamics.jl`. For every (concrete) subtype of `AbstractParticle`, there are two kinds of interface functions implemented: static functions and property functions. The static functions provide information on what kind of particle it is (defaults are written in square brackets) ```julia @@ -14,7 +14,7 @@ The static functions provide information on what kind of particle it is (default is_boson(::AbstractParticle)::Bool [= false] is_particle(::AbstractParticle)::Bool [= true] is_anti_particle(::AbstractParticle)::Bool [= false] -``` +``` If the output of those functions differ from the defaults for a subtype of `AbstractParticle`, these functions need to be overwritten. The second type of functions define a hard interface for `AbstractParticle`: @@ -30,7 +30,7 @@ Base.broadcastable(part::AbstractParticle) = Ref(part) """ AbstractParticleType <: AbstractParticle -This is the abstract base type for every species of particles. All functionalities defined on subtypes of `AbstractParticleType` should be static, i.e. known at compile time. +This is the abstract base type for every species of particles. All functionalities defined on subtypes of `AbstractParticleType` should be static, i.e. known at compile time. For adding runtime information, e.g. four-momenta or particle states, to a particle, consider implementing a concrete subtype of [`AbstractParticle`](@ref) instead, which may have a type parameter `P<:AbstractParticleType`. Concrete built-in subtypes of `AbstractParticleType` are available in `QEDcore.jl` and should always be singletons.. @@ -55,7 +55,7 @@ is_boson(::AbstractParticle) = false """ $(TYPEDSIGNATURES) - + Interface function for particles. Return `true` if the passed subtype of [`AbstractParticle`](@ref) can be considered a *particle* as distinct from anti-particles, and `false` otherwise. The default implementation of `is_particle` for every subtype of [`AbstractParticle`](@ref) will always return `true`. """ @@ -88,28 +88,32 @@ This needs to be implemented for each concrete subtype of [`AbstractParticle`](@ function charge end """ - propagator(particle::AbstractParticleType, mom::QEDbase.AbstractFourMomentum, [mass::Real]) + propagator(particle::AbstractParticleType, mom::QEDbase.AbstractFourMomentum) -Return the propagator of a particle for a given four-momentum. If `mass` is passed, the respective propagator for massive particles is used, if not, it is assumed the particle passed in is massless. +Compute the propagator of a particle for a given four-momentum `mom`. -!!! note "Convention" - - There are two types of implementations for propagators given in `QEDProcesses`: - For a `BosonLike` particle with four-momentum ``k`` and mass ``m``, the propagator is given as +# Notes on Convention +The `QEDProcesses.jl` package includes two types of propagators: - ```math - D(k) = \\frac{1}{k^2 - m^2}. - ``` +**Boson-like particles**: For a `BosonLike` particle with four-momentum `k` and mass `m = QEDbase.mass(particle)`, the propagator is given by: - For a `FermionLike` particle with four-momentum ``p`` and mass ``m``, the propagator is given as +```math +D(k) = \\frac{1}{k^2 - m^2} +``` - ```math - S(p) = \\frac{\\gamma^\\mu p_\\mu + mass}{p^2 - m^2}. - ``` +**Fermion-like particles**: For a `FermionLike` particle with four-momentum `p` and mass `m = QEDbase.mass(particle)`, the propagator is defined as: + +```math +S(p) = \\frac{\\gamma^\\mu p_\\mu + m}{p^2 - m^2} +``` + +Here, ``\\gamma^\\mu`` are the gamma matrices, and ``p_\\mu`` represents the four-momentum components. !!! warning - - This function does not throw when the given particle is off-shell. If an off-shell particle is passed, the function `propagator` returns `Inf`. + + This function does **not** throw an error if the particle is off-shell (i.e., if it + does not satisfy the mass-shell condition). In such cases, the function will return `Inf`, + which indicates that the propagator is undefined for an off-shell particle. """ function propagator end @@ -174,7 +178,7 @@ electron_state = base_state(QEDcore.Electron(), Incoming(), mom, SpinUp()) ``` ```jldoctest -julia> using QEDbase; using QEDcore +julia> using QEDbase; using QEDcore; julia> mass = 1.0; px,py,pz = (0.1, 0.2, 0.3); E = sqrt(px^2 + py^2 + pz^2 + mass^2); mom = SFourMomentum(E, px, py, pz) 4-element SFourMomentum with indices SOneTo(4): @@ -244,7 +248,7 @@ julia> electron_states = base_state(Electron(), Incoming(), mom, AllSpin()) where ``v_\\sigma`` is the base state of the respective outgoing anti-fermion. For a photon with four-momentum ``k^\\mu = \\omega (1, \\cos\\phi \\sin\\theta, \\sin\\phi \\sin\\theta, \\cos\\theta)``, the two polarization vectors are given as - + ```math \\begin{align*} \\epsilon^\\mu_1 &= (0, \\cos\\theta \\cos\\phi, \\cos\\theta \\sin\\phi, -\\sin\\theta),\\\\ diff --git a/src/interfaces/particle_stateful.jl b/src/interfaces/particle_stateful.jl index 9337c4e..0a762ed 100644 --- a/src/interfaces/particle_stateful.jl +++ b/src/interfaces/particle_stateful.jl @@ -3,7 +3,7 @@ AbstractParticleStateful <: QEDbase.AbstractParticle Abstract base type for the representation of a particle with a state. It requires the following interface functions to be provided: - + - [`particle_direction`](@ref): Returning the particle's direction. - [`particle_species`](@ref): Returning the particle's species. - [`momentum`](@ref): Returning the particle's momentum. diff --git a/src/particles/direction.jl b/src/interfaces/particles/direction.jl similarity index 100% rename from src/particles/direction.jl rename to src/interfaces/particles/direction.jl diff --git a/src/particles/spin_pol.jl b/src/interfaces/particles/spin_pol.jl similarity index 77% rename from src/particles/spin_pol.jl rename to src/interfaces/particles/spin_pol.jl index 9695405..bce2f03 100644 --- a/src/particles/spin_pol.jl +++ b/src/interfaces/particles/spin_pol.jl @@ -179,13 +179,40 @@ const PolY = PolarizationY Base.show(io::IO, ::PolY) = print(io, "y-polarized") """ - multiplicity(spin_or_pol) + SyncedPolarization{N::Int} <: AbstractIndefinitePolarization -Return the number of spins or polarizations respresented by `spin_or_pol`, e.g. `multiplicity(SpinUp()) == 1`, but `multiplicity(AllSpin()) = 2`. +An indefinite polarization type, indicating that multiple particles have a synced polarization. +Two polarizations are considered synced when they have the same value for `N`. This means that +the resulting multiplicity will be 2 total for all particles with the same `SyncedPolarization`. +Having a single `SyncedPolarization{N}` in a process is legal. In this case, it behaves just +like an [`AllPolarization`](@ref) would. + +See also: [`multiplicity`](@ref) +""" +struct SyncedPolarization{N} <: AbstractIndefinitePolarization + function SyncedPolarization(N::Int) + return new{N}() + end +end +const SyncedPol = SyncedPolarization +Base.show(io::IO, ::SyncedPolarization{N}) where {N} = print(io, "synced polarization $N") + +""" + SyncedSpin{N::Int} <: AbstractIndefiniteSpin + +An indefinite spin type, indicating that multiple particles have a synced spin. +Two spins are considered synced when they have the same value for `N`. This means that +the resulting multiplicity will be 2 total for all particles with the same `SyncedSpin`. + +Having a single `SyncedSpin{N}` in a process is legal. In this case, it behaves just +like an [`AllSpin`](@ref) would. + +See also: [`multiplicity`](@ref) """ -function multiplicity end -multiplicity(::AbstractDefinitePolarization) = 1 -multiplicity(::AbstractDefiniteSpin) = 1 -multiplicity(::AbstractIndefinitePolarization) = 2 -multiplicity(::AbstractIndefiniteSpin) = 2 +struct SyncedSpin{N} <: AbstractIndefiniteSpin + function SyncedSpin(N::Int) + return new{N}() + end +end +Base.show(io::IO, ::SyncedSpin{N}) where {N} = print(io, "synced spin $N") diff --git a/src/interfaces/phase_space_point.jl b/src/interfaces/phase_space_point.jl index 2a8665a..3a7b01e 100644 --- a/src/interfaces/phase_space_point.jl +++ b/src/interfaces/phase_space_point.jl @@ -102,7 +102,7 @@ end Returns the momentum of the `n`th particle in the given [`AbstractPhaseSpacePoint`](@ref) which has direction `dir` and species `species`. If `n` is outside the valid range for this phase space point, a `BoundsError` is thrown. !!! note - This function accepts n as a `Val{N}` type, i.e., a compile-time constant value (for example a literal `1` or `2`). This allows this function to add zero overhead, but **only** if `N` is actually known at compile time. + This function accepts n as a `Val{N}` type, i.e., a compile-time constant value (for example a literal `1` or `2`). This allows this function to add zero overhead, but **only** if `N` is actually known at compile time. If it is not, use the overload of this function that uses `n::Int` instead. That function is faster than calling this one with `Val(n)`. """ function momentum( @@ -171,12 +171,8 @@ end Return a `Tuple` of all the particles' momenta for the given `ParticleDirection`. """ -function momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Incoming) - return momentum.(particles(psp, QEDbase.Incoming())) -end - -function momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Outgoing) - return momentum.(particles(psp, QEDbase.Outgoing())) +function momenta(psp::AbstractPhaseSpacePoint, dir::ParticleDirection) + return ntuple(i -> momentum(psp[dir, i]), number_particles(process(psp), dir)) end """ diff --git a/src/interfaces/process.jl b/src/interfaces/process.jl index ff226b0..650c584 100644 --- a/src/interfaces/process.jl +++ b/src/interfaces/process.jl @@ -6,9 +6,9 @@ ############### """ -Abstract base type for definitions of scattering processes. It is the root type for the +Abstract base type for definitions of scattering processes. It is the root type for the process interface, which assumes that every subtype of `AbstractProcessDefinition` -implements at least +implements at least ```Julia incoming_particles(proc_def::AbstractProcessDefinition) @@ -17,8 +17,32 @@ outgoing_particles(proc_def::AbstractProcessDefinition) which return a tuple of the incoming and outgoing particles, respectively. -Furthermore, to calculate scattering probabilities and differential cross sections, the following -interface functions need to be implemented for every combination of `CustomProcess<:AbstractProcessDefinition`, +An `AbstractProcessDefinition` is also expected to contain spin and polarization information of its particles. +For this, the functions + +```Julia +incoming_spin_pols(proc_def::AbstractProcessDefinition) +outgoing_spin_pols(proc_def::AbstractProcessDefinition) +``` + +can be overloaded. They must return a tuple of [`AbstractSpinOrPolarization`], where the order must match the order of the process' particles. +A default implementation is provided which assumes [`AllSpin`](@ref) for every [`is_fermion`](@ref) particle and [`AllPolarization`](@ref) for every [`is_boson`](@ref) particle. + +!!! note Performance + It is very beneficial for the performance of derived functions if these functions return compile-time-known values. + +On top of these spin and polarization functions, the following functions are automatically defined: + +```Julia +multiplicity(proc_def::AbstractProcessDefinition) +incoming_multiplicity(proc_def::AbstractProcessDefinition) +outgoing_multiplicity(proc_def::AbstractProcessDefinition) +``` + +Which return the number of spin and polarization combinations that should be considered for the process. For more detail, refer to the functions' documentations. + +Furthermore, to calculate scattering probabilities and differential cross sections, the following +interface functions need to be implemented for every combination of `CustomProcess<:AbstractProcessDefinition`, `CustomModel<:AbstractModelDefinition`, and `CustomPhasespaceDefinition<:AbstractPhasespaceDefinition`. ```Julia @@ -33,7 +57,7 @@ interface functions need to be implemented for every combination of `CustomProce _phase_space_factor(psp::PhaseSpacePoint{CustomProcess,CustomModel,CustomPhasespaceDefinition}) ``` -Optional is the implementation of +Optional is the implementation of ```Julia @@ -41,7 +65,6 @@ Optional is the implementation of ``` to enable the calculation of total probabilities and cross sections. - """ abstract type AbstractProcessDefinition end @@ -53,6 +76,11 @@ Broadcast.broadcastable(proc::AbstractProcessDefinition) = Ref(proc) Interface function for scattering processes. Return a tuple of the incoming particles for the given process definition. This function needs to be given to implement the scattering process interface. + +!!! note Performance + It is very beneficial for the performance of derived functions if this function returns compile-time-known values. + +See also: [`AbstractParticleType`](@ref) """ function incoming_particles end @@ -61,9 +89,40 @@ function incoming_particles end Interface function for scattering processes. Return the tuple of outgoing particles for the given process definition. This function needs to be given to implement the scattering process interface. + +!!! note Performance + It is very beneficial for the performance of derived functions if this function returns compile-time-known values. + +See also: [`AbstractParticleType`](@ref) """ function outgoing_particles end +""" + incoming_spin_pols(proc_def::AbstractProcessDefinition) + +Interface function for scattering processes. Return the tuple of spins or polarizations for the given process definition. The order must be the same as the particles returned from [`incoming_particles`](@ref). +A default implementation is provided, returning [`AllSpin`](@ref) for every [`is_fermion`](@ref) and [`AllPolarization`](@ref) for every [`is_boson`](@ref). + +!!! note Performance + It is very beneficial for the performance of derived functions if this function returns compile-time-known values. + +See also: [`AbstractSpinOrPolarization`](@ref) +""" +function incoming_spin_pols end + +""" + outgoing_spin_pols(proc_def::AbstractProcessDefinition) + +Interface function for scattering processes. Return the tuple of spins or polarizations for the given process definition. The order must be the same as the particles returned from [`outgoing_particles`](@ref). +A default implementation is provided, returning [`AllSpin`](@ref) for every [`is_fermion`](@ref) and [`AllPolarization`](@ref) for every [`is_boson`](@ref). + +!!! note Performance + It is very beneficial for the performance of derived functions if this function returns compile-time-known values. + +See also: [`AbstractSpinOrPolarization`](@ref) +""" +function outgoing_spin_pols end + """ _incident_flux(in_psp::InPhaseSpacePoint{PROC,MODEL}) where { PROC <: AbstractProcessDefinition, @@ -80,14 +139,14 @@ function _incident_flux end MODEL <: AbstractModelDefinition, } -Interface function which returns a tuple of scattering matrix elements for each spin and polarization combination of `proc`. +Interface function which returns a tuple of scattering matrix elements for each spin and polarization combination of `proc`. """ function _matrix_element end """ _averaging_norm(proc::AbstractProcessDefinition) -Interface function, which returns a normalization for the averaging of the squared matrix elements over spins and polarizations. +Interface function, which returns a normalization for the averaging of the squared matrix elements over spins and polarizations. """ function _averaging_norm end @@ -109,16 +168,16 @@ function _is_in_phasespace end PSDEF <: AbstractPhasespaceDefinition, } -Interface function, which returns the pre-differential factor of the invariant phase space intergral measure. +Interface function, which returns the pre-differential factor of the invariant phase space intergral measure. !!! note "Convention" - It is assumed, that this function returns the value of + It is assumed, that this function returns the value of ```math \\mathrm{d}\\Pi_n:= \\prod_{i=1}^N \\frac{\\mathrm{d}^3p_i}{(2\\pi)^3 2 p_i^0} H(P_t, p_1, \\dots, p_N), ``` -where ``H(\\dots)`` is a characteristic function (or distribution) which constrains the phase space, e.g. ``\\delta^{(4)}(P_t - \\sum_i p_i)``. +where ``H(\\dots)`` is a characteristic function (or distribution) which constrains the phase space, e.g. ``\\delta^{(4)}(P_t - \\sum_i p_i)``. """ function _phase_space_factor end @@ -146,82 +205,16 @@ function out_phase_space_dimension end MODEL <: AbstractModelDefinition, } -Interface function for the combination of a scattering process and a physical model. Return the total of a +Interface function for the combination of a scattering process and a physical model. Return the total of a given process and model for a passed `InPhaseSpacePoint`. !!! note "total cross section" - + Given an implementation of this method and [`_incident_flux`](@ref), the respective function for the total cross section [`total_cross_section`](@ref) is also available. """ function _total_probability end -####################### -# -# utility functions -# -####################### - -""" - number_incoming_particles(proc_def::AbstractProcessDefinition) - -Return the number of incoming particles of a given process. -""" -@inline function number_incoming_particles(proc_def::AbstractProcessDefinition) - return length(incoming_particles(proc_def)) -end - -""" - number_outgoing_particles(proc_def::AbstractProcessDefinition) - -Return the number of outgoing particles of a given process. -""" -@inline function number_outgoing_particles(proc_def::AbstractProcessDefinition) - return length(outgoing_particles(proc_def)) -end - -""" - particles(proc_def::AbstractProcessDefinition, ::ParticleDirection) - -Convenience function dispatching to [`incoming_particles`](@ref) or [`outgoing_particles`](@ref) depending on the given direction. -""" -@inline particles(proc_def::AbstractProcessDefinition, ::Incoming) = - incoming_particles(proc_def) -@inline particles(proc_def::AbstractProcessDefinition, ::Outgoing) = - outgoing_particles(proc_def) - -""" - number_particles(proc_def::AbstractProcessDefinition, dir::ParticleDirection) - -Convenience function dispatching to [`number_incoming_particles`](@ref) or [`number_outgoing_particles`](@ref) depending on the given direction, returning the number of incoming or outgoing particles, respectively. -""" -@inline number_particles(proc_def::AbstractProcessDefinition, ::Incoming) = - number_incoming_particles(proc_def) -@inline number_particles(proc_def::AbstractProcessDefinition, ::Outgoing) = - number_outgoing_particles(proc_def) - -""" - number_particles(proc_def::AbstractProcessDefinition, dir::ParticleDirection, species::AbstractParticleType) - -Return the number of particles of the given direction and species in the given process definition. -""" -@inline function number_particles( - proc_def::AbstractProcessDefinition, dir::DIR, species::PT -) where {DIR<:ParticleDirection,PT<:AbstractParticleType} - return count(x -> x isa PT, particles(proc_def, dir)) -end - -""" - number_particles(proc_def::AbstractProcessDefinition, particle::AbstractParticleStateful) - -Return the number of particles of the given particle's direction and species in the given process definition. -""" -@inline function number_particles( - proc_def::AbstractProcessDefinition, ps::AbstractParticleStateful -) - return number_particles(proc_def, particle_direction(ps), particle_species(ps)) -end - ##### # Generation of four-momenta from coordinates # @@ -253,29 +246,3 @@ function _generate_incoming_momenta end Interface function to generate the four-momenta of the outgoing particles from coordinates for a given phase-space definition. """ function _generate_outgoing_momenta end - -""" - _generate_momenta( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::NTuple{N,T}, - out_phase_space::NTuple{M,T}, - ) where {N,M,T<:Real} - -Return four-momenta for incoming and outgoing particles for given coordinate based phase-space points. -""" -function _generate_momenta( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::NTuple{N,T}, - out_phase_space::NTuple{M,T}, -) where {N,M,T<:Real} - in_momenta = _generate_incoming_momenta(proc, model, phase_space_def, in_phase_space) - out_momenta = _generate_outgoing_momenta( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) - - return in_momenta, out_momenta -end diff --git a/src/total_cross_section.jl b/src/total_cross_section.jl deleted file mode 100644 index 84f55db..0000000 --- a/src/total_cross_section.jl +++ /dev/null @@ -1,9 +0,0 @@ -""" - total_cross_section(in_psp::AbstractInPhaseSpacePoint) - -Return the total cross section for a given [`AbstractInPhaseSpacePoint`](@ref). -""" -function total_cross_section(in_psp::AbstractInPhaseSpacePoint) - I = 1 / (4 * _incident_flux(in_psp)) - return I * _total_probability(in_psp) -end diff --git a/test/interfaces/coordinate_transforms.jl b/test/interfaces/coordinate_transforms.jl new file mode 100644 index 0000000..c971655 --- /dev/null +++ b/test/interfaces/coordinate_transforms.jl @@ -0,0 +1,52 @@ +using Random +using QEDbase +using QEDcore + +RNG = MersenneTwister(137137) +ATOL = 0.0 +RTOL = sqrt(eps()) + +include("../test_implementation/TestImplementation.jl") +TESTMODEL = TestImplementation.TestModel() +TESTPSDEF = TestImplementation.TestPhasespaceDef() + +TESTTRAFO = TestImplementation.TestCoordTrafo() + +@testset "broadcast" begin + test_func(trafo) = trafo + @test test_func.(TESTTRAFO) == TESTTRAFO +end + +@testset "single momenta" begin + test_mom = rand(RNG, SFourMomentum) + + test_mom_prime = @inferred TESTTRAFO(test_mom) + + @test isapprox(test_mom_prime, TestImplementation._groundtruth_coord_trafo(test_mom)) +end +@testset "set of momenta" begin + test_moms = rand(RNG, SFourMomentum, 3) + test_moms_prime = TESTTRAFO.(test_moms) + + @test isapprox(test_moms_prime, TestImplementation._groundtruth_coord_trafo.(test_moms)) +end + +@testset "phase space points" begin + @testset "($N_INCOMING,$N_OUTGOING)" for (N_INCOMING, N_OUTGOING) in Iterators.product( + (1, rand(RNG, 2:8)), (1, rand(RNG, 2:8)) + ) + INCOMING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_INCOMING)) + OUTGOING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_OUTGOING)) + + TESTPROC = TestImplementation.TestProcess(INCOMING_PARTICLES, OUTGOING_PARTICLES) + + p_in_phys = TestImplementation._rand_momenta(RNG, N_INCOMING) + p_out_phys = TestImplementation._rand_momenta(RNG, N_OUTGOING) + + PS_POINT = PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys, p_out_phys) + + test_psp_prime = @inferred TESTTRAFO(PS_POINT) + + @test test_psp_prime == TestImplementation._groundtruth_coord_trafo(PS_POINT) + end +end diff --git a/test/interfaces/process.jl b/test/interfaces/process.jl index cc25fdc..a234d65 100644 --- a/test/interfaces/process.jl +++ b/test/interfaces/process.jl @@ -34,6 +34,8 @@ include("../test_implementation/TestImplementation.jl") @testset "failed process interface" begin @test_throws MethodError incoming_particles(TESTPROC_FAIL_ALL) @test_throws MethodError outgoing_particles(TESTPROC_FAIL_ALL) + @test_throws MethodError incoming_spin_pols(TESTPROC_FAIL_ALL) + @test_throws MethodError outgoing_spin_pols(TESTPROC_FAIL_ALL) end @testset "$PROC $MODEL" for (PROC, MODEL) in Iterators.product( @@ -84,6 +86,31 @@ include("../test_implementation/TestImplementation.jl") end end + @testset "incoming/outgoing spins and polarizations" begin + groundtruth_incoming_spin_pols = ntuple( + x -> is_fermion(INCOMING_PARTICLES[x]) ? AllSpin() : AllPolarization(), + N_INCOMING, + ) + groundtruth_outgoing_spin_pols = ntuple( + x -> is_fermion(OUTGOING_PARTICLES[x]) ? AllSpin() : AllPolarization(), + N_OUTGOING, + ) + + @test incoming_spin_pols(TESTPROC) == groundtruth_incoming_spin_pols + @test outgoing_spin_pols(TESTPROC) == groundtruth_outgoing_spin_pols + @test spin_pols(TESTPROC, Incoming()) == groundtruth_incoming_spin_pols + @test spin_pols(TESTPROC, Outgoing()) == groundtruth_outgoing_spin_pols + + for (pt, sp) in Iterators.flatten(( + Iterators.zip(incoming_particles(TESTPROC), incoming_spin_pols(TESTPROC)), + Iterators.zip(outgoing_particles(TESTPROC), outgoing_spin_pols(TESTPROC)), + )) + @test is_boson(pt) ? sp isa AbstractPolarization : true + @test is_fermion(pt) ? sp isa AbstractSpin : true + @test is_boson(pt) || is_fermion(pt) + end + end + @testset "incident flux" begin test_incident_flux = QEDbase._incident_flux( InPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS) @@ -177,3 +204,87 @@ include("../test_implementation/TestImplementation.jl") InPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, ps_in_coords) end end + +@testset "Process Multiplicity" begin + boson = TestImplementation.TestParticleBoson() + fermion = TestImplementation.TestParticleFermion() + + _mult(::AbstractDefinitePolarization) = 1 + _mult(::AbstractDefiniteSpin) = 1 + _mult(::AbstractIndefinitePolarization) = 2 + _mult(::AbstractIndefiniteSpin) = 2 + + @testset "no synced spins/pols" begin # test all possible combinations for fermion+boson->fermion+boson processes, without synced + spins = (SpinUp(), SpinDown(), AllSpin()) + pols = (PolX(), PolY(), AllPol()) + for (p1, s1, p2, s2) in Iterators.product(pols, spins, pols, spins) + proc = TestImplementation.TestProcessSP( + (boson, fermion), (boson, fermion), (p1, s1), (p2, s2) + ) + @test multiplicity(proc) == prod(_mult.((p1, s1, p2, s2))) + @test incoming_multiplicity(proc) == prod(_mult.((p1, s1))) + @test outgoing_multiplicity(proc) == prod(_mult.((p2, s2))) + end + end + + # some special cases for synced spins and pols testing + for i in 1:4 + # i synced bosons + proc = TestImplementation.TestProcessSP( + (ntuple(_ -> boson, i)..., fermion), + (boson, fermion), + (ntuple(_ -> SyncedPolarization(1), i)..., AllSpin()), + (AllPol(), AllSpin()), + ) + @test multiplicity(proc) == 16 + @test incoming_multiplicity(proc) == 4 + @test outgoing_multiplicity(proc) == 4 + + proc = TestImplementation.TestProcessSP( + (boson, fermion), + (ntuple(_ -> boson, i)..., fermion), + (AllPol(), SpinDown()), + (ntuple(_ -> SyncedPolarization(1), i)..., AllSpin()), + ) + @test multiplicity(proc) == 8 + @test incoming_multiplicity(proc) == 2 + @test outgoing_multiplicity(proc) == 4 + + for j in 1:4 + # ... with j synced fermions + proc = TestImplementation.TestProcessSP( + (ntuple(_ -> boson, i)..., fermion), + (boson, ntuple(_ -> fermion, j)), + (ntuple(_ -> SyncedPolarization(1), i)..., SpinDown()), + (PolX(), ntuple(_ -> SyncedSpin(1), j)...), + ) + @test multiplicity(proc) == 4 + @test incoming_multiplicity(proc) == 2 + @test outgoing_multiplicity(proc) == 2 + end + end + + @testset "multiple differing synced polarizations" begin + proc = TestImplementation.TestProcessSP( + (ntuple(_ -> boson, 2)..., fermion), + (ntuple(_ -> boson, 2)..., fermion), + (ntuple(_ -> SyncedPolarization(1), 2)..., SpinUp()), + (ntuple(_ -> SyncedPolarization(2), 2)..., AllSpin()), + ) + @test multiplicity(proc) == 8 + @test incoming_multiplicity(proc) == 2 + @test outgoing_multiplicity(proc) == 4 + end + + @testset "synced polarization across in and out particles" begin + proc = TestImplementation.TestProcessSP( + (ntuple(_ -> boson, 2)..., fermion), + (ntuple(_ -> boson, 2)..., fermion), + (ntuple(i -> SyncedPolarization(i), 2)..., SpinUp()), + (ntuple(i -> SyncedPolarization(i), 2)..., SpinDown()), + ) + @test multiplicity(proc) == 4 + @test incoming_multiplicity(proc) == 4 + @test outgoing_multiplicity(proc) == 4 + end +end diff --git a/test/particle_properties.jl b/test/particle_properties.jl index 74a84c8..9a2e156 100644 --- a/test/particle_properties.jl +++ b/test/particle_properties.jl @@ -1,7 +1,11 @@ using QEDbase +using QEDcore using StaticArrays using Random +include("test_implementation/TestImplementation.jl") +using .TestImplementation: TestParticleBoson, TestParticleFermion + # test function to test scalar broadcasting test_broadcast(x::AbstractParticle) = x test_broadcast(x::ParticleDirection) = x @@ -16,22 +20,53 @@ test_broadcast(x::AbstractSpinOrPolarization) = x @testset "spins and polarization" begin @testset "$spin_or_pol" for spin_or_pol in ( - SpinUp(), SpinDown(), AllSpin(), PolX(), PolY(), AllPol() + SpinUp(), + SpinDown(), + AllSpin(), + PolX(), + PolY(), + AllPol(), + SyncedSpin(1), + SyncedPolarization(1), ) @test test_broadcast.(spin_or_pol) == spin_or_pol end end end -@testset "multiplicity of spins or pols" begin - @testset "single" begin - @testset "$spin_or_pol" for spin_or_pol in (SpinUp(), SpinDown(), PolX(), PolY()) - @test multiplicity(spin_or_pol) == 1 - end - end - @testset "multiple" begin - @testset "$spin_or_pol" for spin_or_pol in (AllSpin(), AllPol()) - @test multiplicity(spin_or_pol) == 2 +TESTPROCS = ( + TestImplementation.TestProcessSP( + (TestParticleBoson(), TestParticleFermion()), + (TestParticleBoson(), TestParticleFermion()), + (AllPol(), AllSpin()), + (AllPol(), AllSpin()), + ), + TestImplementation.TestProcessSP( + (TestParticleBoson(), TestParticleBoson(), TestParticleFermion()), + (TestParticleBoson(), TestParticleFermion()), + (SyncedPol(1), SyncedPol(1), AllSpin()), + (AllPol(), AllSpin()), + ), + TestImplementation.TestProcessSP( + (TestParticleBoson(), TestParticleBoson(), TestParticleFermion()), + (TestParticleBoson(), TestParticleFermion()), + (SyncedPol(1), SyncedPol(2), SyncedSpin(2)), + (SyncedPol(2), SyncedSpin(2)), + ), +) + +@testset "spin_pol iterator ($proc)" for proc in TESTPROCS + @test length(spin_pols_iter(proc)) == multiplicity(proc) + + for combinations in spin_pols_iter(proc) + @test length(combinations) == 2 + in_comb, out_comb = combinations + + @test length(in_comb) == length(incoming_particles(proc)) + @test length(out_comb) == length(outgoing_particles(proc)) + + for sp in Iterators.flatten((in_comb, out_comb)) + @test sp isa AbstractDefiniteSpin || sp isa AbstractDefinitePolarization end end end diff --git a/test/runtests.jl b/test/runtests.jl index 760dae0..c9ef72b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,12 @@ -using QEDbase using Test using SafeTestsets begin # Interfaces + @time @safetestset "coordinate transforms" begin + include("interfaces/coordinate_transforms.jl") + end + @time @safetestset "model interface" begin include("interfaces/model.jl") end diff --git a/test/test_implementation/TestImplementation.jl b/test/test_implementation/TestImplementation.jl index 5197738..1fc49a2 100644 --- a/test/test_implementation/TestImplementation.jl +++ b/test/test_implementation/TestImplementation.jl @@ -1,6 +1,6 @@ """ -This module provides a full implementation of the model and process interface. Its purpose is only for testing and it does not reflect any -real-world physics. +This module provides a full implementation of the model and process interface. Its purpose is only for testing and it does not reflect any +real-world physics. The module exports: @@ -11,8 +11,8 @@ TestParticle3 TestParticle4 TestModel # dummy compute model TestModel_FAIL # failing compute model -TestProcess # dummy scattering process -TestProcess_FAIL # failing scattering process +TestProcess # dummy scattering process +TestProcess_FAIL # failing scattering process TestPhasespaceDef # dummy phase space definition TestPhasespaceDef_FAIL # failing phase space definition ``` @@ -24,6 +24,7 @@ export TestParticle1, TestParticle2, TestParticle3, TestParticle4, PARTICLE_SET export TestModel, TestModel_FAIL export TestProcess, TestProcess_FAIL export TestPhasespaceDef, TestPhasespaceDef_FAIL +export TestCoordTrafo using Random using QEDbase @@ -33,6 +34,7 @@ using StaticArrays include("groundtruths.jl") include("test_model.jl") include("test_process.jl") +include("test_coord_trafo.jl") include("random_momenta.jl") include("utils.jl") diff --git a/test/test_implementation/groundtruths.jl b/test/test_implementation/groundtruths.jl index 44dcf2c..51a85fa 100644 --- a/test/test_implementation/groundtruths.jl +++ b/test/test_implementation/groundtruths.jl @@ -8,7 +8,7 @@ Test implementation of the incident flux. Return the Minkowski square of the sum I = \\left(\\sum p_i\\right)^2, \\end{align} ``` -where \$p_i\\in\\mathrm{ps_in}\$. +where \$p_i\\in\\mathrm{ps_in}\$. """ function _groundtruth_incident_flux(in_ps) s = sum(in_ps) @@ -18,7 +18,7 @@ end """ _groundtruth_matrix_element(in_ps, out_ps) -Test implementation for a matrix elements. Returns a list of three complex numbers without any physical meaning. +Test implementation for a matrix elements. Returns a list of three complex numbers without any physical meaning. """ function _groundtruth_matrix_element(in_ps, out_ps) s_in = sum(in_ps) @@ -156,7 +156,7 @@ end """ _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) -Test implementation of the unsafe differential cross section. Uses the test implementations of `_groundtruth_incident_flux` and `_groundtruth_unsafe_probability`. +Test implementation of the unsafe differential cross section. Uses the test implementations of `_groundtruth_incident_flux` and `_groundtruth_unsafe_probability`. """ function _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) init_flux = _groundtruth_incident_flux(in_ps) @@ -194,7 +194,7 @@ end """ _groundtruth_safe_diffCS(proc, in_ps, out_ps) -Test implementation of the safe differential cross section. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_diffCS`. +Test implementation of the safe differential cross section. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_diffCS`. """ function _groundtruth_safe_diffCS(proc, in_ps, out_ps) if !_groundtruth_is_in_phasespace(in_ps, out_ps) @@ -261,3 +261,19 @@ function _groundtruth_total_cross_section( ) where {N,T<:AbstractFourMomentum} return _groundtruth_total_cross_section.(in_psps) end + +### Coordinate trafos + +_groundtruth_coord_trafo(p::AbstractFourMomentum) = 2 * p + +# coord trafo applied to every momentum in psp +function _groundtruth_coord_trafo(psp::PhaseSpacePoint) + in_moms = momenta(psp, Incoming()) + out_moms = momenta(psp, Outgoing()) + in_moms_prime = _groundtruth_coord_trafo.(in_moms) + out_moms_prime = _groundtruth_coord_trafo.(out_moms) + + return PhaseSpacePoint( + process(psp), model(psp), phase_space_definition(psp), in_moms_prime, out_moms_prime + ) +end diff --git a/test/test_implementation/test_coord_trafo.jl b/test/test_implementation/test_coord_trafo.jl new file mode 100644 index 0000000..a0f31fc --- /dev/null +++ b/test/test_implementation/test_coord_trafo.jl @@ -0,0 +1,3 @@ + +struct TestCoordTrafo <: AbstractCoordinateTransformation end +QEDbase._transform(::TestCoordTrafo, p) = _groundtruth_coord_trafo(p) diff --git a/test/test_implementation/test_process.jl b/test/test_implementation/test_process.jl index 8a31dce..d1706e4 100644 --- a/test/test_implementation/test_process.jl +++ b/test/test_implementation/test_process.jl @@ -22,6 +22,24 @@ end QEDbase.incoming_particles(proc::TestProcess) = proc.incoming_particles QEDbase.outgoing_particles(proc::TestProcess) = proc.outgoing_particles +""" + TestProcessSP + +Process for testing with settable spin and polarization. +""" +struct TestProcessSP{IP<:Tuple,OP<:Tuple,IN_SP<:Tuple,OUT_SP<:Tuple} <: + AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP + incoming_spin_pols::IN_SP + outgoing_spin_pols::OUT_SP +end + +QEDbase.incoming_particles(proc::TestProcessSP) = proc.incoming_particles +QEDbase.outgoing_particles(proc::TestProcessSP) = proc.outgoing_particles +QEDbase.incoming_spin_pols(proc::TestProcessSP) = proc.incoming_spin_pols +QEDbase.outgoing_spin_pols(proc::TestProcessSP) = proc.outgoing_spin_pols + struct TestProcess_FAIL{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition incoming_particles::IP outgoing_particles::OP @@ -55,7 +73,7 @@ function TestProcess_FAIL_ALL(rng::AbstractRNG, N_in::Int, N_out::Int) end """ -Test process with no implemented interface except the incoming and outgoing particles. +Test process with no implemented interface except the incoming and outgoing particles. Should fail every usage except construction of itself and the respective phase space point for given four-momenta. """ struct TestProcess_FAIL_DIFFCS{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition diff --git a/test/test_implementation/utils.jl b/test/test_implementation/utils.jl index 3a141d9..819aa70 100644 --- a/test/test_implementation/utils.jl +++ b/test/test_implementation/utils.jl @@ -1,5 +1,5 @@ -# Check if any failed type is in the input +# Check if any failed type is in the input _any_fail(x...) = true _any_fail(::TestProcess, ::TestModel) = false _any_fail(::TestProcess, ::TestModel, ::TestPhasespaceDef) = false