There are two methods, one in sequential mode(be careful with the assumptions made!) and another in non-sequential mode.
The Strehl ratio is the operand that I can think of to calculate the peak intensity. Though it calculates the ratio, the reference, by which I mean the peak intensity of the diffraction point spread function (PSF) in the absence of aberrations, this value does not vary too much when we just put a small defocus value. Assuming the maximum intensity is A_max, the one at FWHM is A_half, the reference is A_ref, we have the strehl ratio calculated:
at max: S1 = A_max/A_ref
at half: S2 =  A_hafl/A ref
When S2 = 0.5* S1,  A_ref = A_ref
We have A_half = 0.5* A_max.
Now let’s set the system up. 
1.    Insert a dummy surface just before the surface plane.
2.    Set the thickness of this plane using multiconfiguration operand.
3.    Read the strehl ratio at the best focus plane.
4.    Calculate the half value and set it as target, define a weight.
5.    Limit the range of the thickness to one positive and one negative. If the optimization leads to two positive or negative value, change them manually to positive and negative and reoptimize
 

There is another option if you want to get the XZ distribution. That's the detector volume in non-sequential mode. 

One thing that I want to emphasis is that the tilt about Y value is set to 90 so that the detector viewer shows XZ plane instead of XY.
 

The attachment again….