Quantifying reads mapping to multiple loci
You almost had the correct python code already, you just need to filter out secondary alignments:

def get_reads_hist(bam): 
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        if query.is_secondary:
            continue
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count            
    return counts