Hi, I am new to analyzing sequencing data so if the answer to my question is too obvious please let me know to delete that afterwards. I have some bam files and I know they come from paired end sequencing . first I sorted them with samtools and then deduplicated with picard MarkDuplicates and I have a reference peak set ( which is similar to bed format : first column is chrm second is start and third is end)

I want to get a count matrix for bam files i.e. I want to report how many reads fall in each region of peak set. What tools or command should I use for this purpose?

Thank you