You would have to probably do this in several steps to get statistics that would show correlations like you are describing. Depending on how many comparisons you want to make you may have to redo the analysis in several combinations.
I would start by using the Spatial Join tool under the Analysis Tools Toolbox, Overlay Toolset to join of the two feature classes together. The statistical results you can get will depend on which feature class is your Target FC and which is your Join FC and the Join Type you use. I would probably start by making the story category points the target and the landmarks the join and use the JOIN_ONE_TO_MANY option to assign every individual point with the landmark type it overlays.
Once you have the overlay, then you just have to keep doing different summaries using the data to see all of the relationships quantitatively. For example, you can use the Dissolve Tool or Summary Statistics Tool with the story category and landmark type as case fields. If you do Dissolve, be sure to get a count of one of the fields so that you can sort on the count to find the pair of types that get the highest number of hits. You could include the Join_FID in the case fields if you wanted to know the polygon with the highest hits for each story type.
To do a wide array of statistics without the need to create multiple outputs, you would need to learn Python. With Python it is possible to use cursors, dictionaries, lists, numpy, etc. to compile a wide array of statistics by just doing one or two passes on the overlay data and output them to several formats. But you would need to define the kinds of statistics you really are looking for and make an attempt at doing the coding yourself to get much help on a script designed for your needs.