par(mar=c(4, 4, 0, 0) + 0.2)

R <- function (E, A, Z) 0.02442 * A / Z^2 * (E / A)^1.75
E <- function (R, A, Z) 8.3424 * A * (R * Z^2 / A)^0.57143

Th1 <- 0.3	#mm
Th2 <- 3	#mm

E1 <- function (E0, A, Z) {
	ifelse (R(E0, A, Z) > Th1, E0 - E(R(E0, A, Z) - Th1, A, Z), E0)
}

E2 <- function (E0, A, Z) {
	Er <- ifelse (R(E0, A, Z) > Th1, E(R(E0, A, Z) - Th1, A, Z), 0)
	ifelse (R(Er, A, Z) > Th2, Er - E(R(Er, A, Z) - Th2, A, Z), Er)
}

E0 <- seq (1, 2000, by=1)

plot (c(0,350), c(0,100), type="n", xlab="Thick detector (MeV)", 
	ylab="thin detector (MeV)")
grid ()

mylwd <- 1
lines (E2(E0, 12, 6), E1(E0, 12, 6), lwd=mylwd )
lines (E2(E0, 11, 6), E1(E0, 11, 6), lwd=mylwd )
lines (E2(E0, 10, 6), E1(E0, 10, 6), lwd=mylwd )

lines (E2(E0, 11, 5), E1(E0, 11, 5), col="magenta", lwd=mylwd )
lines (E2(E0, 10, 5), E1(E0, 10, 5), col="magenta", lwd=mylwd )
lines (E2(E0, 8, 5), E1(E0, 8, 5), col="magenta", lwd=mylwd )

lines (E2(E0, 10, 4), E1(E0, 10, 4), col="orange", lwd=mylwd )
lines (E2(E0, 9, 4), E1(E0, 9, 4), col="orange", lwd=mylwd )
lines (E2(E0, 7, 4), E1(E0, 7, 4), col="orange", lwd=mylwd )

lines (E2(E0, 7, 3), E1(E0, 7, 3), col="green", lwd=mylwd )
lines (E2(E0, 6, 3), E1(E0, 6, 3), col="green", lwd=mylwd )

lines (E2(E0, 4, 2), E1(E0, 4, 2), col="blue", lwd=mylwd )
lines (E2(E0, 3, 2), E1(E0, 3, 2), col="blue", lwd=mylwd )

lines (E2(E0, 3, 1), E1(E0, 3, 1), col="red", lwd=mylwd )
lines (E2(E0, 2, 1), E1(E0, 2, 1), col="red", lwd=mylwd )
lines (E2(E0, 1, 1), E1(E0, 1, 1), col="red", lwd=mylwd )

Z.list <- c(1, 1, 1, 2, 2, 3, 3, 4, 4, 4,  5, 5,  5,  6,  6,  6)
A.list <- c(1, 2, 3, 3, 4, 6, 7, 7 ,9, 10, 8, 10, 11, 10, 11, 12)

part.ID <- sample(1:length(Z.list), 1000, replace=TRUE)
E.list <- runif (length(part.ID), 0, 2000) 

Eth<- E1(E.list, A.list[part.ID], Z.list[part.ID])
ETh<- E2(E.list, A.list[part.ID], Z.list[part.ID])


plot (c(0,350), c(0,100), type="n", xlab="Thick detector (MeV)", 
	ylab="thin detector (MeV)")
grid ()

points (ETh, Eth)