-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsample_fig._for_ed4.11.2018.R
7 lines (7 loc) · 1.82 KB
/
sample_fig._for_ed4.11.2018.R
1
2
3
4
5
6
7
# The following chunk from plots.R in Sachs' sandbox produce an actual figure for the minor paper, namely Fig. 3.2.3 for MS06.docx. The figure is for Fe56 (600 MeV/u) DER, Si28 DER, and corresponding IEA and SEA 50-50 MIXDERS. This script and looking at Fig. 3.2.3, which is very similar (except that simple effect additivity is a little higher) but used the old data, not the revised data put in .csv files last week, may help Edward write his code, including Monte Carlo with ribbon plots corresponding to the same mixture (one or two of the ribbon plots will also appear in the minor paper).
forty_cGy=0:41 # with all factors of 0.01 removed from the regression step and corresponding definitions, this definition happens to give the correct horizontal scale.
plot(x = forty_cGy, y =100* calculate_SEA(forty_cGy, c(1/2, 1/2), c(70, 195), n = 2),bty='l', type = "l", xlab = "Dose (cGy)", ylab = "HG Prevalence", col = "black", lwd = 2, lty = 2) # with all factors of 0.01 removed from the regression step and corresponding definitions, the extra factor of 100 above happens to give the correct vertical scale. Of course the final scripts will be written so that this extra factor never need appear, but for the next few weeks it should serve as a temporary man-hole cover. Same comment applies to the next few lines.
lines(x = forty_cGy, y =100* calib_HZE_nte_ider(dose = forty_cGy, L = 70), col = "cyan", lwd = 2)
lines(x = forty_cGy, y = 100*calib_HZE_nte_ider(dose = forty_cGy, L = 195), col = "orange", lwd = 2)
lines(x = forty_cGy, y =100* calculate_complex_id(r = c(0.5 , 0.5), L = c(70, 195), d = forty_cGy, model = "NTE", lowLET = FALSE)[, 2], col = "red", lwd = 2) # I(d)
# no legend is added because other fig. codes already show how to make legends and for this example Edward will be able to get all the relevant information from this script.