We construct a new factorized waveform including $(l,|m|)=(2,2),(2,1),(3,3),(4,4)$ modes based on effective-one-body (EOB) formalism, which is valid for spinning binary black holes (BBH) in general equatorial orbit. When combined with the dynamics of $texttt{SEOBNRv4}$, the $(l,|m|)=(2,2)$ mode waveform generated by this new waveform can fit the original $texttt{SEOBNRv4}$ waveform very well in the case of a quasi-circular orbit. We have calibrated our new waveform model to the Simulating eXtreme Spacetimes (SXS) catalog. The comparison is done for BBH with total mass in $(20,200)M_odot$ using Advanced LIGO designed sensitivity. For the quasi-circular cases we have compared our $(2,2)$ mode waveforms to the 281 numerical relativity (NR) simulations of BBH along quasi-circular orbits. All of the matching factors are bigger than 98%. For the elliptical cases, 24 numerical relativity simulations of BBH along an elliptic orbit are used. For each elliptical BBH system, we compare our modeled gravitational polarizations against the NR results for different combinations of the inclination angle, the initial orbit phase and the source localization in the sky. We use the the minimal matching factor respect to the inclination angle, the initial orbit phase and the source localization to quantify the performance of the higher modes waveform. We found that after introducing the high modes, the minimum of the minimal matching factor among the 24 tested elliptical BBHs increases from 90% to 98%. Following our previous $texttt{SEOBNRE}$ waveform model, we call our new waveform model $texttt{SEOBNREHM}$. Our $texttt{SEOBNREHM}$ waveform model can match all tested 305 SXS waveforms better than 98% including highly spinning ($chi=0.99$) BBH, highly eccentric ($eapprox0.15$) BBH and large mass ratio ($q=10$) BBH.