Modeling halo mass functions in chameleon f(R) gravity