The best bound solution for logged best bounds in MIQCP problem

I am using a callback to print the variables of the best bound solution. As I am interested in showing reliably that in my maximisation problem (which runs too long to wait for Gurobi to reach optimality with 0 MIPGap) the optimal solution does not exceed 1, for sanity checks I would like to also retrieve the point yielding the best bound at termination (I use MIPGap = 0.01 as a criterion as I verified that it is enough, and only takes a few minutes). The callback I wrote follows the following principles, from a comment to a similar question:

“There is no function to get the best bound solution whenever a new incumbent solution has been found. However, with a bit of bookkeeping, you can still manage to have the best bound solution point at hand. In the MIPNODE callback, whenever a better best bound has been found, you can save the corresponding best bound solution point obtained by MIPNODE_REL. You then have to replace the best bound solution point and its value only when a better bound has been reported in the MIPNODE callback. This way, you will have the best bound solution at hand whenever a new incumbent has been found.”

In my model I have a list of variables s and a list of variables b, the objective is s Q b, with Q a matrix, so I implement the above as follows:

def bbsol(model, where):
   if where == GRB.Callback.MIPNODE:
       status = model.cbGet(GRB.Callback.MIPNODE_STATUS)
       if status == GRB.OPTIMAL:
           bdsols = model.cbGetNodeRel(s)
           bdsolb = model.cbGetNodeRel(b)
           val = bilinobj(Q, bdsols, bdsolb)
           objbst = model.cbGet(GRB.Callback.MIPNODE_OBJBST)
           objbnd = model.cbGet(GRB.Callback.MIPNODE_OBJBND)
           model._data1.append([val, objbnd, objbst, bdsols, bdsolb])
           if len(model._data1) > 1:
               last = model._data1[-1]
               sec_last = model._data1[-2]
               if last[1] < sec_last[1] and last[0] >= last[2]:
                   m._data2.append(last)

m = gp.Model("bilinear")
m._data1 = [] #to save all solution at nodes 
m._data2 = []#to save only solutions obtained when the best bound is improved

…

s = [m.addVar(lb = 0.0, ub = 1.0, name = "s" + str(i)) for i in range(n)]
b = [m.addVar(lb = 0.0, ub = 1.0, name = "b" + str(i)) for i in range(n)]
…

m.Params.NonConvex = 2
m.Params.MIPGap = 0.01
m.Params.MIPFocus = 3 
m.Params.Threads = 1

m.optimize(bbsol)

What happens is that I obtain 

Optimal solution found (tolerance 1.00e-02)
Best objective 9.887624893280e-01, best bound 9.986500364876e-01, gap 1.0000%

with 

model._data2[-1] = 

[0.9889000617474336,
0.9986507617178466,
0.9887624932077008,
[0.421918568925866, 0.045379488835395124, 0.0, 0.0, 0.013398837396792306, 0.019794087498992242, 0.02375423679547895, 0.003823706819213857, 0.2235221942555108, 0.008443883944574536, 0.13409883899851766, 0.02079631120600453, 0.03498870963667133, 0.05008014866548969, 0.0, 9.870214930047734e-07],
[0.0, 0.0, 0.0, 0.0, 0.06150209368528896, 0.06736965016691181, 0.07268181934955725, 0.05244509648581034, 0.416197035059513, 0.01777218674225835, 0.1920972105255384, 0.030336263356890627, 0.0424362439785228, 0.04673825836432657, 0.0, 0.00042414228538176533]]

which is not what I want, since the objective function evaluates at the best bound solution as 0.98890, which is smaller than the best bound ultimately found: 0.99865. What I got is clearly not the best bound solution for 0.99865, as it is just a leaf node which evaluates at something better than the best bound at the moment that the best bound has improved. However the improvement is not due to the current leaf, rather to the argmax amongst the leafs, if I understand correctly the Branch and Bound description. How can I get my hands on this argmax? I do not necessarily need the very last one, which might be tricky due to the final clean up phase, the second last logged best bound solution is fine too. 

It would be necessary to label the solutions appended with them being the leafs of each stage of the branch and bound procedure, so that I could run the argmax only on the leafs of the current stage: is there a method to do that? Or a way to avoid it?