By the way, since I "fixed" the calculation of offsets (see code above) the gyro drift appears to be gone.
I did not need to apply any of the filtering (gAlpha calculations) to achieve that.
Assuming I am correct in this, there is a difference between your offset calculation (which averages many hundreds of samples, read as quickly as possible) and mine (which reads many hundreds of sample, over a fixed time interval - 5 seconds).