So, basically, what you are saying is that I have to generate signal+background and background only toys, and compare it with the data, in order to extract the maximum likelihood values, for a specific test_mu, and then proceed and calculate the corresponding p-values, according to the observed value of the test statistic?
Almost!
You have to generate signal+background where the signal is scaled by a given test \mu and background only toys, to obtain test statistic distributions for the background hypothesis and signal hypothesis for a given \mu. You then compute the test statistic for the data and then proceed and calculate the corresponding p-values, according to the observed value of the test statistic. From the p-values, you get the CLs.
But here is the kicker: to get the 95 % CL upper limit, you have to iteratively repeat this procedure for many values of the test \mu to find the value for which the CLs is 0.95, just as it says at the end of section 2.1 in the paper! That means generating tons of toys in total, which can take a very long time. That’s why people usually use the asymptotic approximations, like explained by @wernerd in this thread:
https://root-forum.cern.ch/t/cls-upper-limits-without-roofit-roostats-packages