With the increasing availability of functional genomic data, incorporating genomic annotations into QTL mapping has become a standard analytical procedure. However, the existing analysis methods often lack rigor and/or computational efficiency. We present a novel algorithm to perform integrative multi-SNP QTL mapping in a probabilistic hierarchical model framework that enables accurate and efficient joint enrichment analysis and the identification of multiple causal variants. The software implementation of the method is freely available at https://github.com/xqwen/dap