Cox proportional hazard, entire population as reference In R when I use coxph the first factor is automatically picked as the reference group:
df <- list(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              treatment=c('A','B','B','C','C','A','A'))
summary(coxph(Surv(time, status) ~ treatment, df)) 


And I can obtain HR and CI for treatment B and C, relative to treatment A. How can I get HR, lower95, upper95 for all 3 treatments relative to the entire cohort (average survival)?

I wouldn't recommend that. In your example, no one got an "average" treatment of 3/7 A plus 2/7 B plus 2/7 C. There's no assurance, even if anyone got such an "average" treatment, that she would have a net result in the log-hazard scale that represents the weighted mean of the individual log-hazards. Or, in a perfectly sex-balanced cohort, no individual is exactly half-male and half-female. The most useful comparisons are typically relative to a realistic scenario with covariate values that make sense.
The above holds for any model. Therneau and Grambsch in Chapter 10 discuss further issues with respect to "Expected survival." They note that a survival curve based on a hypothetical patient having cohort-averaged values of covariates is not representative of the cohort as a whole.
If you nevertheless want to get hazard ratios of the type you request, the method is the same as for any regression model. You evaluate linear functions of the coefficient estimates that represent the difference between each treatment and the average over the cohort. In your case, the average is 3/7 A plus 2/7 B plus 2/7 C. Calculations are done in the log-hazard coefficient scale, in which that for the reference level A is 0 and for the other treatments are:
cox1 <- coxph(Surv(time, status) ~ treatment, df)
coef(cox1)
# treatmentB treatmentC 
#   1.483427   1.419022

Standard errors in the log-hazard scale (and corresponding 95% CI) are based on the formula for the variance of a linear combination of random variables. That requires knowing the covariance matrix of the coefficient estimates, which you can get from vcov(cox1). Post-modeling tools like those in the emmeans package help. For example:
library(emmeans)
emm1 <- emmeans(cox1,"treatment")
coxContrasts <- contrast(emm1,list(A=c(4/7,-2/7,-2/7),B= c(-3/7,5/7,-2/7), C=c(-3/7,-2/7,5/7)))

The contrasts in that last line are related to the formulas for the difference of each treatment group's coefficient from the mean of the log-hazard coefficients. Then you can get:
confint(coxContrasts)
#  contrast estimate    SE  df asymp.LCL asymp.UCL
#  A          -0.829 0.683 Inf    -2.169      0.51
#  B           0.654 0.736 Inf    -0.789      2.10
#  C           0.590 0.916 Inf    -1.206      2.39
# 
# Results are given on the log (not the response) scale. 
# Confidence level used: 0.95 

Exponentiate the estimate, asymp.LCL and asymp.UCL to get the corresponding hazard ratios relative to the (unrealistic) mean and the confidence intervals. Or specify type="response" in the call to confint().
