diff --git a/NAMESPACE b/NAMESPACE
index c4a64e9..989308e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,11 +1,16 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(primary_censored_cdf,default)
+S3method(primary_censored_cdf,pcens_numeric)
+S3method(primary_censored_cdf,pcens_pgamma_dunif)
+S3method(primary_censored_cdf,pcens_plnorm_dunif)
 export(check_dprimary)
 export(check_pdist)
 export(dexpgrowth)
 export(dpcens)
 export(dprimarycensoreddist)
 export(fitdistdoublecens)
+export(new_primary_censored_dist)
 export(pcd_as_cmdstan_data)
 export(pcd_cmdstan_model)
 export(pcd_load_stan_functions)
@@ -15,6 +20,7 @@ export(pcd_stan_path)
 export(pexpgrowth)
 export(ppcens)
 export(pprimarycensoreddist)
+export(primary_censored_cdf)
 export(rexpgrowth)
 export(rpcens)
 export(rprimarycensoreddist)
diff --git a/NEWS.md b/NEWS.md
index 2b40e7e..626bbc2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,6 +6,7 @@ This is the current development version.
 
 * Add `{touchstone}` based benchmarks for benchmarking R utility functions, and fitting the `stan` and `fitdistplus` models.
 * Added a "How it works" vignette.
+* Added R infrastructure for analytical solutions via the `primary_censored_dist` S3 class.
 
 # primarycensoreddist 0.4.0
 
diff --git a/R/dprimarycensoreddist.R b/R/dprimarycensoreddist.R
index c860775..71a2bb4 100644
--- a/R/dprimarycensoreddist.R
+++ b/R/dprimarycensoreddist.R
@@ -70,7 +70,8 @@
 dprimarycensoreddist <- function(
     x, pdist, pwindow = 1, swindow = 1,
     D = Inf, dprimary = stats::dunif,
-    dprimary_args = list(), log = FALSE, ...) {
+    dprimary_args = list(), log = FALSE,
+    pdist_name = NULL, dprimary_name = NULL, ...) {
   check_pdist(pdist, D, ...)
   check_dprimary(dprimary, pwindow, dprimary_args)
 
@@ -82,6 +83,13 @@ dprimarycensoreddist <- function(
     )
   }
 
+  if (is.null(pdist_name)) {
+    pdist_name <- .extract_function_name(substitute(pdist))
+  }
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(substitute(dprimary))
+  }
+
   # Compute CDFs for all unique points
   unique_points <- sort(unique(c(x, x + swindow)))
   unique_points <- unique_points[unique_points > 0]
@@ -90,7 +98,8 @@ dprimarycensoreddist <- function(
   }
 
   cdfs <- pprimarycensoreddist(
-    unique_points, pdist, pwindow, Inf, dprimary, dprimary_args, ...
+    unique_points, pdist, pwindow, Inf, dprimary, dprimary_args,
+    pdist_name = pdist_name, dprimary_name = dprimary_name, ...
   )
 
   # Create a lookup table for CDFs
diff --git a/R/pprimarycensoreddist.R b/R/pprimarycensoreddist.R
index 85e3ad1..a7a9033 100644
--- a/R/pprimarycensoreddist.R
+++ b/R/pprimarycensoreddist.R
@@ -26,7 +26,19 @@
 #' @param dprimary_args List of additional arguments to be passed to
 #' dprimary. For example, when using `dexpgrowth`, you would
 #' pass `list(min = 0, max = pwindow, r = 0.2)` to set the minimum, maximum,
-#' and rate parameters.
+#' and rate parameters
+#'
+#' @param pdist_name A string specifying the name of the delay distribution
+#' function. If NULL, the function name is extracted using
+#' [.extract_function_name()]. Used to determine if a analytical solution
+#' exists for the primary censored distribution. Must be set if `pdist` is
+#' passed a pre-assigned variable rather than a function name.
+#'
+#' @param dprimary_name A string specifying the name of the primary event
+#' distribution function. If NULL, the function name is extracted using
+#' [.extract_function_name()]. Used to determine if a analytical solution
+#' exists for the primary censored distribution. Must be set if `dprimary` is
+#' passed a pre-assigned variable rather than a function name.
 #'
 #' @param ... Additional arguments to be passed to pdist
 #'
@@ -61,7 +73,19 @@
 #' }
 #' where \eqn{F_{\text{cens,norm}}(q)} is the normalized CDF.
 #'
+#' This function creates a `primary_censored_dist` object using
+#' [new_primary_censored_dist()] and then computes the primary event
+#' censored CDF using [primary_censored_cdf()]. This abstraction allows
+#' for automatic use of analytical solutions when available, while
+#' seamlessly falling back to numerical integration when necessary.
+#'
+#' Note: For analytical detection to work correctly, `pdist` and `dprimary`
+#' must be directly passed as distribution functions, not via assignment or
+#' `pdist_name` and `dprimary_name` must be used to override the default
+#' extraction of the function name.
+#'
 #' @family primarycensoreddist
+#' @seealso [new_primary_censored_dist()] and [primary_censored_cdf()]
 #'
 #' @examples
 #' # Example: Lognormal distribution with uniform primary events
@@ -75,25 +99,29 @@
 #' )
 pprimarycensoreddist <- function(
     q, pdist, pwindow = 1, D = Inf, dprimary = stats::dunif,
-    dprimary_args = list(), ...) {
+    dprimary_args = list(), pdist_name = NULL, dprimary_name = NULL, ...) {
   check_pdist(pdist, D, ...)
   check_dprimary(dprimary, pwindow, dprimary_args)
 
-  result <- vapply(q, function(d) {
-    if (d <= 0) {
-      return(0) # Return 0 for non-positive delays
-    } else {
-      integrand <- function(p) {
-        d_adj <- d - p
-        pdist(d_adj, ...) *
-          do.call(
-            dprimary, c(list(x = p, min = 0, max = pwindow), dprimary_args)
-          )
-      }
+  if (is.null(pdist_name)) {
+    pdist_name <- .extract_function_name(substitute(pdist))
+  }
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(substitute(dprimary))
+  }
 
-      stats::integrate(integrand, lower = 0, upper = pwindow)$value
-    }
-  }, numeric(1))
+  # Create a new primary_censored_dist object
+  pcens_obj <- new_primary_censored_dist(
+    pdist,
+    dprimary,
+    dprimary_args,
+    pdist_name = pdist_name,
+    dprimary_name = dprimary_name,
+    ...
+  )
+
+  # Compute the CDF using the S3 method
+  result <- primary_censored_cdf(pcens_obj, q, pwindow)
 
   if (!is.infinite(D)) {
     # Compute normalization factor for finite D
diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
new file mode 100644
index 0000000..303cfb9
--- /dev/null
+++ b/R/primary_censored_dist.R
@@ -0,0 +1,162 @@
+#' S3 class for primary event censored distribution computation
+#'
+#' @inheritParams pprimarycensoreddist
+#'
+#' @return An object of class primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+new_primary_censored_dist <- function(
+    pdist, dprimary, dprimary_args,
+    pdist_name = NULL,
+    dprimary_name = NULL, ...) {
+  if (is.null(pdist_name)) {
+    pdist_name <- .extract_function_name(substitute(pdist))
+  }
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(substitute(dprimary))
+  }
+
+  structure(
+    list(
+      pdist = pdist,
+      dprimary = dprimary,
+      dprimary_args = dprimary_args,
+      args = list(...)
+    ),
+    class = c(
+      paste0(
+        "pcens_",
+        pdist_name, "_",
+        dprimary_name
+      )
+    )
+  )
+}
+
+#' Compute primary event censored CDF
+#'
+#' @inheritParams pprimarycensoreddist
+#'
+#' @param object A `primary_censored_dist` object as created by
+#' [new_primary_censored_dist()].
+#'
+#' @param use_numeric Logical, if TRUE forces use of numeric integration
+#' even for distributions with analytical solutions. This is primarily
+#' useful for testing purposes or for settings where the analytical solution
+#' breaks down.
+#'
+#' @return Vector of primary event censored CDFs
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  UseMethod("primary_censored_cdf")
+}
+
+#' Default method for computing primary event censored CDF
+#'
+#' This method serves as a fallback for combinations of delay and primary
+#' event distributions that don't have specific implementations. It uses
+#' the numeric integration method.
+#'
+#' @inheritParams primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.default <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+}
+
+#' Numeric method for computing primary event censored CDF
+#'
+#' This method uses numerical integration to compute the primary event censored
+#' CDF for any combination of delay distribution and primary event distribution.
+#'
+#' @inheritParams primary_censored_cdf
+#' @inheritParams pprimarycensoreddist
+#'
+#' @details
+#' This method implements the numerical integration approach for computing
+#' the primary event censored CDF. It uses the same mathematical formulation
+#' as described in the details section of [pprimarycensoreddist()], but
+#' applies numerical integration instead of analytical solutions.
+#'
+#' @seealso [pprimarycensoreddist()] for the mathematical details of the
+#' primary event censored CDF computation.
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.pcens_numeric <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  result <- vapply(q, function(d) {
+    if (d <= 0) {
+      return(0) # Return 0 for non-positive delays
+    } else {
+      integrand <- function(p) {
+        d_adj <- d - p
+        do.call(object$pdist, c(list(q = d_adj), object$args)) *
+          do.call(
+            object$dprimary,
+            c(list(x = p, min = 0, max = pwindow), object$dprimary_args)
+          )
+      }
+
+      stats::integrate(integrand, lower = 0, upper = pwindow)$value
+    }
+  }, numeric(1))
+
+  return(result)
+}
+
+#' Method for Gamma delay with uniform primary
+#'
+#' @inheritParams primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.pcens_pgamma_dunif <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  use_numeric <- TRUE
+  if (isTRUE(use_numeric)) {
+    return(
+      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+    )
+  }
+
+  result <- vapply(q, function(n) {
+    # Implement analytical solution here
+  }, numeric(1))
+
+  return(result)
+}
+
+#' Method for Log-Normal delay with uniform primary
+#'
+#' @inheritParams primary_censored_cdf
+#'
+#' @family primary_censored_dist
+#'
+#' @export
+primary_censored_cdf.pcens_plnorm_dunif <- function(
+    object, q, pwindow, use_numeric = FALSE) {
+  use_numeric <- TRUE
+  if (isTRUE(use_numeric)) {
+    return(
+      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+    )
+  }
+
+  result <- vapply(q, function(n) {
+    # Implement analytical solution here
+  }, numeric(1))
+
+  return(result)
+}
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000..0914d5d
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,15 @@
+#' Extract Base Function Name
+#'
+#' This helper function extracts the base name of a function, removing an
+#' namespace prefixes.
+#'
+#' @param func The output of `substitute` on a function.
+#'
+#' @return A character string representing the base name of the function.
+#'
+#' @keywords internal
+.extract_function_name <- function(func) {
+  func_name <- deparse(func)
+  base_name <- sub("^.*::", "", func_name)
+  return(base_name)
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 0923943..24c2429 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -43,6 +43,10 @@ reference:
   desc: Probability density and random generation functions for primary event distributions
   contents:
   - has_concept("primaryeventdistributions")
+- title: Primary censored distribution class and methods
+  desc: S3 class and methods for computing primary event censored distributions, focusing on the internal machinery used by the package. Unlike the primary event distributions section which deals with specific distribution functions, this section covers the general framework for handling censored distributions.
+  contents:
+  - has_concept("primary_censored_dist")
 - title: Distribution checking functions
   desc: Functions to validate cumulative distribution functions (CDFs) and probability density functions (PDFs)
   contents:
diff --git a/inst/WORDLIST b/inst/WORDLIST
index de393d5..a1e4485 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -15,6 +15,7 @@ PMFs
 Poisson
 SamuelBrand
 athowes
+cdf
 cdot
 cens
 cmdstan
@@ -23,6 +24,7 @@ disccensprob
 disccensunifprim
 discretisation
 discretising
+dist
 dp
 dprimary
 dprimarycensoreddist
diff --git a/man/check_dprimary.Rd b/man/check_dprimary.Rd
index 5078f69..935b9eb 100644
--- a/man/check_dprimary.Rd
+++ b/man/check_dprimary.Rd
@@ -20,7 +20,7 @@ See \code{primary_dists.R} for examples.}
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
-and rate parameters.}
+and rate parameters}
 
 \item{tolerance}{The tolerance for the integral to be considered close to 1}
 }
diff --git a/man/dot-dpcens.Rd b/man/dot-dpcens.Rd
index f2c2437..7d4fc60 100644
--- a/man/dot-dpcens.Rd
+++ b/man/dot-dpcens.Rd
@@ -30,7 +30,7 @@ See \code{primary_dists.R} for examples.}
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
-and rate parameters.}
+and rate parameters}
 
 \item{...}{Additional arguments to be passed to the distribution function}
 }
diff --git a/man/dot-extract_function_name.Rd b/man/dot-extract_function_name.Rd
new file mode 100644
index 0000000..63a63c8
--- /dev/null
+++ b/man/dot-extract_function_name.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.extract_function_name}
+\alias{.extract_function_name}
+\title{Extract Base Function Name}
+\usage{
+.extract_function_name(func)
+}
+\arguments{
+\item{func}{The output of \code{substitute} on a function.}
+}
+\value{
+A character string representing the base name of the function.
+}
+\description{
+This helper function extracts the base name of a function, removing an
+namespace prefixes.
+}
+\keyword{internal}
diff --git a/man/dot-ppcens.Rd b/man/dot-ppcens.Rd
index 9001167..3988afc 100644
--- a/man/dot-ppcens.Rd
+++ b/man/dot-ppcens.Rd
@@ -27,7 +27,7 @@ See \code{primary_dists.R} for examples.}
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
-and rate parameters.}
+and rate parameters}
 
 \item{...}{Additional arguments to be passed to pdist}
 }
diff --git a/man/dprimarycensoreddist.Rd b/man/dprimarycensoreddist.Rd
index e908eec..0608eb2 100644
--- a/man/dprimarycensoreddist.Rd
+++ b/man/dprimarycensoreddist.Rd
@@ -14,6 +14,8 @@ dprimarycensoreddist(
   dprimary = stats::dunif,
   dprimary_args = list(),
   log = FALSE,
+  pdist_name = NULL,
+  dprimary_name = NULL,
   ...
 )
 
@@ -26,6 +28,8 @@ dpcens(
   dprimary = stats::dunif,
   dprimary_args = list(),
   log = FALSE,
+  pdist_name = NULL,
+  dprimary_name = NULL,
   ...
 )
 }
@@ -52,10 +56,22 @@ See \code{primary_dists.R} for examples.}
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
-and rate parameters.}
+and rate parameters}
 
 \item{log}{Logical; if TRUE, probabilities p are given as log(p)}
 
+\item{pdist_name}{A string specifying the name of the delay distribution
+function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{pdist} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
 \item{...}{Additional arguments to be passed to the distribution function}
 }
 \value{
diff --git a/man/fitdistdoublecens.Rd b/man/fitdistdoublecens.Rd
index 3bd4191..7e3f63a 100644
--- a/man/fitdistdoublecens.Rd
+++ b/man/fitdistdoublecens.Rd
@@ -38,7 +38,7 @@ See \code{primary_dists.R} for examples.}
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
-and rate parameters.}
+and rate parameters}
 
 \item{...}{Additional arguments to be passed to \code{\link[fitdistrplus:fitdist]{fitdistrplus::fitdist()}}.}
 }
diff --git a/man/new_primary_censored_dist.Rd b/man/new_primary_censored_dist.Rd
new file mode 100644
index 0000000..7c68b33
--- /dev/null
+++ b/man/new_primary_censored_dist.Rd
@@ -0,0 +1,60 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/primary_censored_dist.R
+\name{new_primary_censored_dist}
+\alias{new_primary_censored_dist}
+\title{S3 class for primary event censored distribution computation}
+\usage{
+new_primary_censored_dist(
+  pdist,
+  dprimary,
+  dprimary_args,
+  pdist_name = NULL,
+  dprimary_name = NULL,
+  ...
+)
+}
+\arguments{
+\item{pdist}{Distribution function (CDF)}
+
+\item{dprimary}{Function to generate the probability density function
+(PDF) of primary event times. This function should take a value \code{x} and a
+\code{pwindow} parameter, and return a probability density. It should be
+normalized to integrate to 1 over [0, pwindow]. Defaults to a uniform
+distribution over [0, pwindow]. Users can provide custom functions or use
+helper functions like \code{dexpgrowth} for an exponential growth distribution.
+See \code{primary_dists.R} for examples.}
+
+\item{dprimary_args}{List of additional arguments to be passed to
+dprimary. For example, when using \code{dexpgrowth}, you would
+pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
+and rate parameters}
+
+\item{pdist_name}{A string specifying the name of the delay distribution
+function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{pdist} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{...}{Additional arguments to be passed to pdist}
+}
+\value{
+An object of class primary_censored_cdf
+}
+\description{
+S3 class for primary event censored distribution computation
+}
+\seealso{
+Low level primary event censored distribution objects and methods
+\code{\link{primary_censored_cdf}()},
+\code{\link{primary_censored_cdf.default}()},
+\code{\link{primary_censored_cdf.pcens_numeric}()},
+\code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
+\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
+}
+\concept{primary_censored_dist}
diff --git a/man/pprimarycensoreddist.Rd b/man/pprimarycensoreddist.Rd
index ca0400e..bbd4abf 100644
--- a/man/pprimarycensoreddist.Rd
+++ b/man/pprimarycensoreddist.Rd
@@ -12,6 +12,8 @@ pprimarycensoreddist(
   D = Inf,
   dprimary = stats::dunif,
   dprimary_args = list(),
+  pdist_name = NULL,
+  dprimary_name = NULL,
   ...
 )
 
@@ -22,6 +24,8 @@ ppcens(
   D = Inf,
   dprimary = stats::dunif,
   dprimary_args = list(),
+  pdist_name = NULL,
+  dprimary_name = NULL,
   ...
 )
 }
@@ -46,7 +50,19 @@ See \code{primary_dists.R} for examples.}
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
-and rate parameters.}
+and rate parameters}
+
+\item{pdist_name}{A string specifying the name of the delay distribution
+function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{pdist} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
 
 \item{...}{Additional arguments to be passed to pdist}
 }
@@ -82,6 +98,17 @@ by \eqn{F_{\text{cens}}(D)}:
 F_{\text{cens,norm}}(q) = \frac{F_{\text{cens}}(q)}{F_{\text{cens}}(D)}
 }
 where \eqn{F_{\text{cens,norm}}(q)} is the normalized CDF.
+
+This function creates a \code{primary_censored_dist} object using
+\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}} and then computes the primary event
+censored CDF using \code{\link[=primary_censored_cdf]{primary_censored_cdf()}}. This abstraction allows
+for automatic use of analytical solutions when available, while
+seamlessly falling back to numerical integration when necessary.
+
+Note: For analytical detection to work correctly, \code{pdist} and \code{dprimary}
+must be directly passed as distribution functions, not via assignment or
+\code{pdist_name} and \code{dprimary_name} must be used to override the default
+extraction of the function name.
 }
 \examples{
 # Example: Lognormal distribution with uniform primary events
@@ -95,6 +122,8 @@ pprimarycensoreddist(
 )
 }
 \seealso{
+\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}} and \code{\link[=primary_censored_cdf]{primary_censored_cdf()}}
+
 Primary event censored distribution functions
 \code{\link{dprimarycensoreddist}()},
 \code{\link{rprimarycensoreddist}()}
diff --git a/man/primary_censored_cdf.Rd b/man/primary_censored_cdf.Rd
new file mode 100644
index 0000000..4dd2d7d
--- /dev/null
+++ b/man/primary_censored_cdf.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/primary_censored_dist.R
+\name{primary_censored_cdf}
+\alias{primary_censored_cdf}
+\title{Compute primary event censored CDF}
+\usage{
+primary_censored_cdf(object, q, pwindow, use_numeric = FALSE)
+}
+\arguments{
+\item{object}{A \code{primary_censored_dist} object as created by
+\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}}.}
+
+\item{q}{Vector of quantiles}
+
+\item{pwindow}{Primary event window}
+
+\item{use_numeric}{Logical, if TRUE forces use of numeric integration
+even for distributions with analytical solutions. This is primarily
+useful for testing purposes or for settings where the analytical solution
+breaks down.}
+}
+\value{
+Vector of primary event censored CDFs
+}
+\description{
+Compute primary event censored CDF
+}
+\seealso{
+Low level primary event censored distribution objects and methods
+\code{\link{new_primary_censored_dist}()},
+\code{\link{primary_censored_cdf.default}()},
+\code{\link{primary_censored_cdf.pcens_numeric}()},
+\code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
+\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
+}
+\concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.default.Rd b/man/primary_censored_cdf.default.Rd
new file mode 100644
index 0000000..2359b20
--- /dev/null
+++ b/man/primary_censored_cdf.default.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/primary_censored_dist.R
+\name{primary_censored_cdf.default}
+\alias{primary_censored_cdf.default}
+\title{Default method for computing primary event censored CDF}
+\usage{
+\method{primary_censored_cdf}{default}(object, q, pwindow, use_numeric = FALSE)
+}
+\arguments{
+\item{object}{A \code{primary_censored_dist} object as created by
+\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}}.}
+
+\item{q}{Vector of quantiles}
+
+\item{pwindow}{Primary event window}
+
+\item{use_numeric}{Logical, if TRUE forces use of numeric integration
+even for distributions with analytical solutions. This is primarily
+useful for testing purposes or for settings where the analytical solution
+breaks down.}
+}
+\description{
+This method serves as a fallback for combinations of delay and primary
+event distributions that don't have specific implementations. It uses
+the numeric integration method.
+}
+\seealso{
+Low level primary event censored distribution objects and methods
+\code{\link{new_primary_censored_dist}()},
+\code{\link{primary_censored_cdf}()},
+\code{\link{primary_censored_cdf.pcens_numeric}()},
+\code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
+\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
+}
+\concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.pcens_numeric.Rd b/man/primary_censored_cdf.pcens_numeric.Rd
new file mode 100644
index 0000000..b43974f
--- /dev/null
+++ b/man/primary_censored_cdf.pcens_numeric.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/primary_censored_dist.R
+\name{primary_censored_cdf.pcens_numeric}
+\alias{primary_censored_cdf.pcens_numeric}
+\title{Numeric method for computing primary event censored CDF}
+\usage{
+\method{primary_censored_cdf}{pcens_numeric}(object, q, pwindow, use_numeric = FALSE)
+}
+\arguments{
+\item{object}{A \code{primary_censored_dist} object as created by
+\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}}.}
+
+\item{q}{Vector of quantiles}
+
+\item{pwindow}{Primary event window}
+
+\item{use_numeric}{Logical, if TRUE forces use of numeric integration
+even for distributions with analytical solutions. This is primarily
+useful for testing purposes or for settings where the analytical solution
+breaks down.}
+}
+\description{
+This method uses numerical integration to compute the primary event censored
+CDF for any combination of delay distribution and primary event distribution.
+}
+\details{
+This method implements the numerical integration approach for computing
+the primary event censored CDF. It uses the same mathematical formulation
+as described in the details section of \code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}}, but
+applies numerical integration instead of analytical solutions.
+}
+\seealso{
+\code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}} for the mathematical details of the
+primary event censored CDF computation.
+
+Low level primary event censored distribution objects and methods
+\code{\link{new_primary_censored_dist}()},
+\code{\link{primary_censored_cdf}()},
+\code{\link{primary_censored_cdf.default}()},
+\code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
+\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
+}
+\concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.pcens_pgamma_dunif.Rd b/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
new file mode 100644
index 0000000..1ee894f
--- /dev/null
+++ b/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/primary_censored_dist.R
+\name{primary_censored_cdf.pcens_pgamma_dunif}
+\alias{primary_censored_cdf.pcens_pgamma_dunif}
+\title{Method for Gamma delay with uniform primary}
+\usage{
+\method{primary_censored_cdf}{pcens_pgamma_dunif}(object, q, pwindow, use_numeric = FALSE)
+}
+\arguments{
+\item{object}{A \code{primary_censored_dist} object as created by
+\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}}.}
+
+\item{q}{Vector of quantiles}
+
+\item{pwindow}{Primary event window}
+
+\item{use_numeric}{Logical, if TRUE forces use of numeric integration
+even for distributions with analytical solutions. This is primarily
+useful for testing purposes or for settings where the analytical solution
+breaks down.}
+}
+\description{
+Method for Gamma delay with uniform primary
+}
+\seealso{
+Low level primary event censored distribution objects and methods
+\code{\link{new_primary_censored_dist}()},
+\code{\link{primary_censored_cdf}()},
+\code{\link{primary_censored_cdf.default}()},
+\code{\link{primary_censored_cdf.pcens_numeric}()},
+\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
+}
+\concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.pcens_plnorm_dunif.Rd b/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
new file mode 100644
index 0000000..387070f
--- /dev/null
+++ b/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/primary_censored_dist.R
+\name{primary_censored_cdf.pcens_plnorm_dunif}
+\alias{primary_censored_cdf.pcens_plnorm_dunif}
+\title{Method for Log-Normal delay with uniform primary}
+\usage{
+\method{primary_censored_cdf}{pcens_plnorm_dunif}(object, q, pwindow, use_numeric = FALSE)
+}
+\arguments{
+\item{object}{A \code{primary_censored_dist} object as created by
+\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}}.}
+
+\item{q}{Vector of quantiles}
+
+\item{pwindow}{Primary event window}
+
+\item{use_numeric}{Logical, if TRUE forces use of numeric integration
+even for distributions with analytical solutions. This is primarily
+useful for testing purposes or for settings where the analytical solution
+breaks down.}
+}
+\description{
+Method for Log-Normal delay with uniform primary
+}
+\seealso{
+Low level primary event censored distribution objects and methods
+\code{\link{new_primary_censored_dist}()},
+\code{\link{primary_censored_cdf}()},
+\code{\link{primary_censored_cdf.default}()},
+\code{\link{primary_censored_cdf.pcens_numeric}()},
+\code{\link{primary_censored_cdf.pcens_pgamma_dunif}()}
+}
+\concept{primary_censored_dist}
diff --git a/man/roxygen/meta.R b/man/roxygen/meta.R
index 7aee250..ed60282 100644
--- a/man/roxygen/meta.R
+++ b/man/roxygen/meta.R
@@ -1,6 +1,8 @@
 list( # nolint
   rd_family_title = list( # nolint
     primarycensoreddist = "Primary event censored distribution functions",
+    primary_censored_dist =
+      "Low level primary event censored distribution objects and methods",
     primaryeventdistributions = "Primary event distributions",
     check = "Distribution checking functions",
     stantools = "Tools for working with package Stan functions",
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
new file mode 100644
index 0000000..921a14b
--- /dev/null
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -0,0 +1,144 @@
+test_that("new_primary_censored_dist creates object with correct structure", {
+  pdist_name <- "pgamma"
+  pdist <- pgamma
+  dprimary_name <- "dunif"
+  dprimary <- dunif
+  shape <- 2
+  rate <- 1
+
+  obj <- new_primary_censored_dist(
+    pdist,
+    dprimary, list(),
+    pdist_name, dprimary_name,
+    shape = shape, rate = rate
+  )
+
+  expect_s3_class(obj, "pcens_pgamma_dunif")
+  expect_identical(obj$pdist, pgamma)
+  expect_identical(obj$dprimary, dunif)
+  expect_identical(obj$args, list(shape = shape, rate = rate))
+
+  new_obj <- new_primary_censored_dist(
+    pgamma, dunif, list(),
+    shape = shape, rate = rate
+  )
+  expect_identical(obj, new_obj)
+})
+
+test_that("primary_censored_cdf.pcens_numeric computes correct values", {
+  pdist_name <- "pgamma"
+  pdist <- pgamma
+  dprimary_name <- "dunif"
+  dprimary <- dunif
+  shape <- 2
+  rate <- 1
+
+  obj <- new_primary_censored_dist(
+    pdist,
+    dprimary, list(),
+    pdist_name, dprimary_name,
+    shape = shape, rate = rate
+  )
+
+  q_values <- c(0, 5, 10, 15)
+  pwindow <- 10
+
+  result <- primary_censored_cdf(obj, q = q_values, pwindow = pwindow)
+
+  expect_true(all(result >= 0))
+  expect_true(all(result <= 1))
+  expect_type(result, "double")
+  expect_length(result, length(q_values))
+  expect_true(all(diff(result) >= 0)) # Ensure CDF is non-decreasing
+})
+
+test_that("primary_censored_cdf methods dispatch correctly", {
+  pdist_name <- "pgamma"
+  pdist <- pgamma
+  dprimary_name <- "dunif"
+  dprimary <- dunif
+
+  obj_gamma <- new_primary_censored_dist(
+    pdist, dprimary, list(),
+    pdist_name, dprimary_name,
+    shape = 2, rate = 1
+  )
+
+  pdist_name <- "plnorm"
+  pdist <- plnorm
+  dprimary_name <- "dunif"
+  dprimary <- dunif
+
+  obj_lnorm <- new_primary_censored_dist(
+    pdist, dprimary, list(),
+    pdist_name, dprimary_name,
+    meanlog = 0, sdlog = 1
+  )
+
+  expect_s3_class(obj_gamma, "pcens_pgamma_dunif")
+  expect_s3_class(obj_lnorm, "pcens_plnorm_dunif")
+
+  q_values <- c(5, 10)
+  pwindow <- 10
+
+  expect_no_error(
+    primary_censored_cdf(obj_gamma, q = q_values, pwindow = pwindow)
+  )
+  expect_no_error(
+    primary_censored_cdf(obj_lnorm, q = q_values, pwindow = pwindow)
+  )
+})
+
+test_that("primary_censored_cdf.pcens_pgamma_dunif uses numeric method", {
+  pdist_name <- "pgamma"
+  pdist <- pgamma
+  dprimary_name <- "dunif"
+  dprimary <- dunif
+  shape <- 2
+  rate <- 1
+
+  obj <- new_primary_censored_dist(
+    pdist,
+    dprimary, list(),
+    pdist_name, dprimary_name,
+    shape = shape, rate = rate
+  )
+
+  q_values <- c(5, 10)
+  pwindow <- 10
+
+  expect_identical(
+    primary_censored_cdf(obj, q = q_values, pwindow = pwindow),
+    primary_censored_cdf(
+      obj,
+      q = q_values, pwindow = pwindow, use_numeric = TRUE
+    )
+  )
+})
+
+test_that("primary_censored_cdf.pcens_plnorm_dunif uses numeric method", {
+  pdist_name <- "plnorm"
+  pdist <- plnorm
+  dprimary_name <- "dunif"
+  dprimary <- dunif
+  meanlog <- 0
+  sdlog <- 1
+
+  obj <- new_primary_censored_dist(
+    pdist,
+    dprimary, list(),
+    pdist_name, dprimary_name,
+    meanlog = meanlog, sdlog = sdlog
+  )
+
+  q_values <- c(5, 10)
+  pwindow <- 10
+
+  expect_identical(
+    primary_censored_cdf(obj, q = q_values, pwindow = pwindow),
+    primary_censored_cdf(
+      obj,
+      q = q_values, pwindow = pwindow, use_numeric = TRUE
+    )
+  )
+})