Massive stars play crucial roles in determining the physical and chemical evolution of galaxies. They shape their environment from early in their protostellar phase when they blast the surrounding with powerful jets, up until their violent deaths in the form of supernova. However, they form deeply embedded in their parental clouds, making it challenging to directly observe these stars and immediate environments. Notwithstanding, their massive outflows can extend several parsec and since accretion and ejection processes are intrinsically related, they can provide crucial information about the processes governing massive star formation very close to the central engine and reveal the elusive massive protostar. In this talk, I will walk you through an observational journey in massive star-forming regions, focusing in the near infrared including data from the HST, VLT, and LBT. I will also present a newly created python package, called sedcreator, to fit spectral energy distributions specifically developed for massive protostars.