More Dynamic Workflow

[jkguiang]

At the moment, it seems that the workflow for this analysis framework is essentially based around a configuration JSON. While this is nice for simplicity, it might also suffer from its rigidity in the future. For example, I know that we have to do many checks and small studies for VBS HWW wherein I have to kind of brutalize my code in order to get timely results. As such, I wonder if you would be open to moving towards a different workflow, one that we've been using for VBS HWW and that I've been messing around with on the side. I think it could lend some flexibility that would prove generally useful.

The credit for this workflow is due to Philip. He, for a long time, has been organizing his analyses into "cutflow" objects. These cutflows are constructed by stringing together "cuts" in a tree structure. Each cut is given a name, and two lambda functions: one that contains the logic for the cut (i.e. returns pass or fail selection) and another that returns the event weight for that cut (e.g. xsec weight, b-tagging scale factors). Overall, this workflow allows for a lot of niceties, and I have listed a few under "Pros" at the bottom that it also shares. However, while I think it is nice, I think it can also be simplified and improved.

I would like to propose the following overall structure. It seems that a binary search tree is more fitting, in general, than the multi-branch tree that we use in VBS HWW. That is, each cut is a node in the BST, and whether it returns true or false determines if it goes right or left (respectively) in the tree. The iteration in the event loop then terminates when a leaf is reached. For example, I took Philip's cutflow idea (and related tools) and put them into something I call RAPIDO. I would take a similar approach within Python, where a lot of this stuff would be way simpler due to Python's dynamic typing (a lot of acrobatics is needed to achieve a similar effect in C++, and it is needed for writing to TTrees). That is, I would use a similar structure to what you have now, where selections are separated into their own functions, etc. However, I would change the overall organization towards this cutflow structure.

Let me conclude with a quick summary of the workflow along with a few "pros" I would like to highlight.

Workflow:
I propose that you organize the analysis like a BST cutflow, where you have one (or a few) common cutflow(s) in this repo. Contributors then clone this repo and make necessary changes to these common cutflows to answer questions, do weird studies, etc. Common object IDs and other tools (e.g. Python equivalent of cmstas/NanoTools) are also kept here.

Pros:

• More explicitly (i.e. pep8/general gospel of "explicit over implicit") organized, less hidden behind nested objects
• Cutflows are easily printable
• Simple/diagnostic histograms can be filled while looping over events at different stages of the cutflow
• Different signal regions or control regions are more easily/explicitly definable in the BST framework
• The BDTs that you ultimately produce can also be exported to this BST format and analyzed as above

Apologies for the long issue. Please let me know if this is at all interesting to you; I would be happy to meet with folks to discuss further.

[sam-may]

Hi Jonathan, thanks for the suggestion & detailed explanation.

This seems like a nice way to organize things and even if we don't rework the HggAnalysisDev code to include this, it could be a really nice feature for the "official" Hgg nanoAOD framework (where I think it's more worth the effort to have things done in the most optimal way).

Let me see if I understand: you would propose replacing e.g. these cuts

HggAnalysisDev/Preselection/selections/diphoton_selections.py

Lines 12 to 54 in 45de57c

	### mgg cut ###
	resonant = options["resonant"]
	if resonant:
	mgg_mask = numpy.array(events.gg_mass > options["diphotons"]["mgg_lower"]) & numpy.array(events.gg_mass < options["diphotons"]["mgg_upper"])
	else:
	sideband_low = numpy.array(events.gg_mass > options["diphotons"]["mgg_lower"]) & numpy.array(events.gg_mass < options["diphotons"]["mgg_sideband_lower"])
	sideband_high = numpy.array(events.gg_mass > options["diphotons"]["mgg_sideband_upper"]) & numpy.array(events.gg_mass < options["diphotons"]["mgg_upper"])
	mgg_mask = sideband_low | sideband_high
	lead_pt_mgg_cut = selected_photons.pt[:,0] / events.gg_mass > options["photons"]["lead_pt_mgg_cut"]
	sublead_pt_mgg_cut = selected_photons.pt[:,1] / events.gg_mass > options["photons"]["sublead_pt_mgg_cut"]
	pt_mgg_cut = lead_pt_mgg_cut & sublead_pt_mgg_cut
	lead_idmva_cut = selected_photons.mvaID[:,0] > options["photons"]["idmva_cut"]
	sublead_idmva_cut = selected_photons.mvaID[:,1] > options["photons"]["idmva_cut"]
	idmva_cut = lead_idmva_cut & sublead_idmva_cut
	lead_eveto_cut = selected_photons.electronVeto[:,0] > options["photons"]["eveto_cut"]
	sublead_eveto_cut = selected_photons.electronVeto[:,1] > options["photons"]["eveto_cut"]
	eveto_cut = lead_eveto_cut & sublead_eveto_cut
	lead_eta_cut1 = abs(selected_photons.eta[:,0]) < options["photons"]["eta"]
	lead_eta_cut2 = abs(selected_photons.eta[:,0]) < options["photons"]["transition_region_eta"][0]
	lead_eta_cut3 = abs(selected_photons.eta[:,0]) > options["photons"]["transition_region_eta"][1]
	lead_eta_cut = lead_eta_cut1 & (lead_eta_cut2 | lead_eta_cut3)
	sublead_eta_cut1 = abs(selected_photons.eta[:,1]) < options["photons"]["eta"]
	sublead_eta_cut2 = abs(selected_photons.eta[:,1]) < options["photons"]["transition_region_eta"][0]
	sublead_eta_cut3 = abs(selected_photons.eta[:,1]) > options["photons"]["transition_region_eta"][1]
	sublead_eta_cut = sublead_eta_cut1 & (sublead_eta_cut2 | sublead_eta_cut3)
	eta_cut = lead_eta_cut & sublead_eta_cut
	if options["data"]:
	if options["year"] == 2016:
	trigger_cut = events.HLT_Diphoton30_18_R9Id_OR_IsoCaloId_AND_HE_R9Id_Mass90 == True
	elif options["year"] == 2017 or options["year"] == 2018:
	trigger_cut = events.HLT_Diphoton30_22_R9Id_OR_IsoCaloId_AND_HE_R9Id_Mass90 == True
	else:
	trigger_cut = events.gg_mass > 0
	all_cuts = mgg_mask & pt_mgg_cut & idmva_cut & eveto_cut & eta_cut & trigger_cut
	cut_diagnostics.add_cuts([mgg_mask, pt_mgg_cut, idmva_cut, eveto_cut, eta_cut, trigger_cut, all_cuts], ["mgg in [100, 180]" if resonant else "mgg in [100, 120] or [130, 180]", "lead (sublead) pt/mgg > 0.33 (0.25)", "pho idmva > -0.7", "eveto cut", "eta cut", "trigger", "all"])

with a class that gives each cut a name, a cut function, and a weight function. Is that true?

The binary search tree bit would be a bit tricky to add given that we are using columnar tools (i.e. we don't always do an explicit loop through the events). (Or maybe I miss something here?)

[jkguiang]

Yes, that's about right, though I'm envisioning defining a cut that encompasses, for example, one cut for the resonant preselection and another for the non-resonant preselection. The cut name can actually just be the function name, and more information can be contained in the doc string, so a very nice and verbose cutflow table is easier to make in python, thanks to the accessibility of these attributes in Python.

To respond to your second comment, I think the central idea would be as follows. The cutflow object would have a method called run that takes the events as an argument. This function would pass the full list of events to the root cut. Those events that pass the root cut are passed to the left node, whereas those that fail are passed to the right, and so on. Of course, with more thought, perhaps a more convenient or otherwise nice alternative to this could be designed.

[bsathian]

This sounds cool! @sam-may do we have plans of adding this in the HiggsDNA framework?