First, suppose $|T| = 0$. The number of times you toss every coin for each coin to turn up heads is the maximum of $N$ IID geometric random variables with success probability $1-p$. (There are two conventions for geometric random variables differing by $1$. Some people count failures before the first success, but here we need tosses up to and including the first success.) Let the average number of tosses needed be $a(N,p)$.

If we increase $N$ and $|T|$ by the same amount, then some of the tosses don't count because at least one of the specified subset doesn't come up heads. It takes an average of $p^{-|T|}$ tosses to get one which counts. The average number of tosses before you have selected $T$ is $a(n,p) p^{-|T|}$, where $n = N-|T|$. Indeed, the whole distribution can be viewed as a negative binomial distribution compounded with the maximum of IID geometric distributions: First take the maximum of IID geometric random variables $X$, and then sample from a negative binomial distribution with success rate $p^{|T|}$ and count $X$.

To optimize, we would like to know $a(n,p)$. The exact value is complicated in the way I commented, but it does have some useful estimates, and some exact formulas. See Mike Spivey's answers to "What is the expected maximum out of a sample (size N) from a geometric distribution?" A lot of useful information can be found in Bennett Eisenberg, 2008. "On the expectation of the maximum of IID geometric random variables." *Statistics & Probability Letters* 78 (2008) 135–143. This includes short proofs of many earlier results. Here is an exact formula:

$$a(n,p) = \sum_{k=1}^n (-1)^{k+1} \frac{n \choose k}{1-p^k}$$

Unfortunately, this alternating sum seems hard to analyze directly. Here is an easy estimate by viewing geometric random variables as exponential random variables with corresponding rates rounded up. The average maximum of IID exponentials is easy to calculate. Let $\lambda = -\log p, p = \exp(-\lambda)$.

$$\frac{1}{\lambda} H_n \lt a(n,p) \lt 1+\frac{1}{\lambda} H_n$$

where $H_n = \sum_{k=1}^n \frac{1}{k}$ is the expected maximum of IID standard exponential random variables. $H_n = \log n + \gamma +o(1)$ where $\gamma$ is Euler's gamma constant. In a complicated sense, $a(n,p)$ may be viewed as approximately the midpoint $\frac{1}{2} + \frac{1}{\lambda}H_n$, although the error term does not go to $0$ as $n \to \infty$. The absolute magnitude of the error is small, and has a Cesàro mean of $0$.

If we pretend that $a(n,p) = \frac{1}{2} + \frac{1}{\lambda}(\log n + \gamma)$, then we can optimize $a(n,p) p^{-|T|}$. Using Mathematica I find one maximum on $[0,1]$:

$$p = n \exp (\gamma - \sqrt{(\gamma + \log n)(\gamma +\log n + 2/|T|)})$$

For example, if $n=100, |T|=1$, this estimates the optimal value of $p$ as $0.399$. If $n=100, |T|=2$, then $p = 0.620$. It appears that the optimal value of $p$ is much more sensitive to $|T|$ than to $n$, it is roughly a value which makes the specified subset all heads with probability $1/e$, $p \sim 1-\frac{1}{|T|+1/2}$, but it is slightly higher for small values of $n$. Since the optimal $p$ is relatively insensitive to $n$, allowing yourself to change $p$ based on $n$ in future rounds should slightly decrease the optimal value of $p$ for any given $n$.

For small values of $n$ and $|T|$, we can optimize the exact value, and compare this with the approximation. For $n=|T|=10$, the optimum value of $p$ is $0.906333$. The approximation gives $p=0.906357$. These result in averages of $80.965470$ and $80.965473$ rounds of tosses, respectively, a difference of about $3\times 10^{-6}$.