One of the theoretically rigorous approaches to transform the quantum dynamics into trajectories is the Feynman's path integral formalism. Based on the stationary phase approximation of the propagator expressed in the functional integral form, a series of semiclassical dynamics methods have been developed. In the traditional framework of semiclassical dynamics, a swarm of classical trajectories with phases and monodromy matrices are generated to characterize the quantum fluctuations. Although semiclassical dynamics has a rigorous theoretical background and can provide more intuitive understanding of quantum dynamics based on classical trajectories, the superposition of trajectories with phases suffers from the "dynamical sign problem". Thereby, an enormous number of trajectories may be needed to achieve convergence of the results, which hampers the use of semiclassical dynamics in large systems and long-time simulations. To mitigate this numerical instability and derived a decoherence time formula from the quantum Liouville equation and proposed the augmented FSSH (A-FSSH) method, which is considered as one of the most robust decoherence schemes for FSSH. By comparing the equations of motion between the quantum-classical Liouville equation (QCLE) and FSSH, they have made a comprehensive derivation to interpret FSSH from a more rigorous theoretical foundation. A trajectory analysis method was also proposed to obtain coherence between states, which showed encouraging improvements in calculating diabatic populations in the spin-boson model. Recently, we have proposed the branching corrected surface hopping (BCSH) method to incorporate decoherence to the FSSH scheme without introducing additional parameters. As demonstrated in a large number of scattering models, BCSH recovers almost the same branching ratios as the exact quantum dynamics. Furthermore, BCSH gives highly accurate time-dependent populations of both adiabatic and diabatic states when incorporated with a proper trajectory analysis strategy. In this study, we propose a novel trajectory analysis method to consistently construct the full vibronic density matrix from surface hopping trajectories, which uses both the occupation of active states and the electronic wavefunctions. The validity of our new method is demonstrated in Tully's standard scattering models and a two-dimensional scattering model proposed by Subotnik. The results show that BCSH combined with the new trajectory analysis method can produce very close time-dependent vibronic density matrix compared to the exact quantum results.