This paper deals with the application of the adjoint transport theory in order to optimize Monte Carlo based radiotherapy treatment planning. The technique is applied to Boron Neutron Capture Therapy where most often mixed beams of neutrons and gammas are involved. In normal forward Monte Carlo simulations the particles start at a source and lose energy as they travel towards the region of interest, i.e., the designated point of detection. Conversely, with adjoint Monte Carlo simulations, the so-called adjoint particles start at the region of interest and gain energy as they travel towards the source where they are detected. In this respect, the particles travel backwards and the real source and real detector become the adjoint detector and adjoint source, respectively. At the adjoint detector, an adjoint function is obtained with which numerically the same result, e.g., dose or flux in the tumor, can be derived as with forward Monte Carlo. In many cases, the adjoint method is more efficient and by that is much quicker when, for example, the response in the tumor or organ at risk for many locations and orientations of the treatment beam around the patient is required. However, a problem occurs when the treatment beam is mono-directional as the probability of detecting adjoint Monte Carlo particles traversing the beam exit ͑detector plane in adjoint mode͒ in the negative direction of the incident beam is zero. This problem is addressed here and solved first with the use of next event estimators and second with the application of a Legendre expansion technique of the angular adjoint function. In the first approach, adjoint particles are tracked deterministically through a tube to a ͑adjoint͒ point detector far away from the geometric model. The adjoint particles will traverse the disk shaped entrance of this tube ͑the beam exit in the actual geometry͒ perpendicularly. This method is slow whenever many events are involved that are not contributing to the point detector, e.g., neutrons in a scattering medium. In the second approach, adjoint particles that traverse an adjoint shaped detector plane are used to estimate the Legendre coefficients for expansion of the angular adjoint function. This provides an estimate of the adjoint function for the direction normal to the detector plane. In a realistic head model, as described in this paper, which is surrounded by 1020 mono-directional neutron/gamma beams and from which the best ones are to be selected, the example calculates the neutron and gamma fluxes in ten tumors and ten organs at risk. For small diameter beams ͑5 cm͒, and with comparable relative errors, forward Monte Carlo is seen to be 1.5 times faster than the adjoint Monte Carlo techniques. For larger diameter neutron beams ͑10 and 15 cm͒, the Legendre technique is found to be 6 and 20 times faster, respectively. In the case of gammas alone, for the 10 and 15 cm diam beams, both adjoint Monte Carlo Legendre and point detector techniques are respectively 2 and 3 times faster than forward Monte Carlo.