Post hoc analysis survival data and maximum lifespan survival correction

I every body!

Before to present my problem I apologize that I creates the same topic on the forum stackoverflow, but I think the best place for my question was here.

I'm a PhD student and I'm working on survival data currently, but I don't feel comfortable with survival analysis. I have to correlated problem I think I try to explain it in details

I my case, I'm working with samples in 6 cities which are in 2 treatments conditions (control and chemical stress). I haven't censored data.

The biological question is how genetic diversity (illustrated by the city) can interact with environment modifications (the treatment)?

To answer I created a survival parametric model with weibull distribution (AIC method to evaluate best models).

And I tested first if in control condition only, survival was different.

    data <- subset(survival.data, survival.data\$treatment=="treated")
surv <- with(data, Surv(time,event))
weibull.model <- survreg(surv~city,dist="weibull", data=data)


The ouput model confirm significant differences (p.val<0.05) and this is my first problem. If I want to compare differences between cities, how can I be sure this is a treatment effect or just because a lineage have a longer longevity in mean? My idea was to correct it by the maximum lifespan mean for every city before to use the model, but I didn't see any information for that. I discovered the relative survival but I think it's not my case. I tried it and I transformed my data which are represented as % of total life per city [treated]/[max(untreated)]*100

If I resolved this point next analysis should be differents

If not, I will work withclassical survival data but I want to test the interaction term which represent the genotype x environnement interaction With ANOVA I can have this term and use a post-hoc test as TukeyHSD function in R, but here the output of the model use a parameter of the factors as reference and I can't have the interaction.

Call:
coxph(formula = surv~treatment*city,dist="weibull", data=data)

coef exp(coef) se(coef)     z       p
treatmentpq                   5.178   177.335    0.568  9.11 < 2e-16
cityparis                    -0.674     0.510    0.866 -0.78 0.43671
citysapporo                   2.354    10.524    0.539  4.37 1.2e-05
citysokol                     0.277     1.320    0.671  0.41 0.67931
citytokyo                     1.893     6.641    0.553  3.42 0.00062
citywatsonville               1.389     4.009    0.572  2.43 0.01522
treatmentpq:cityparis         0.644     1.904    0.904  0.71 0.47630
treatmentpq:citysapporo      -1.707     0.181    0.601 -2.84 0.00450
treatmentpq:citysokol        -1.085     0.338    0.720 -1.51 0.13191
treatmentpq:citytokyo        -2.098     0.123    0.613 -3.42 0.00062
treatmentpq:citywatsonville  -0.542     0.581    0.631 -0.86 0.39009

Likelihood ratio test=488.6  on 11 df, p=<2e-16
n= 360, number of events= 249


In forum, publications or google I didn't see this type of analysis... Is it possible with survival data ? Maybe I can modelize it with glm approach but i'm not sure.

Thank you! ps: Sorry for English native speaker if I'm unclear...