Multi-Dimensional Fits with pyhf: A Minimal Example Using Flattened 3D Histograms

[JieWu-GitHub]

Hello pyhf community,

I’m Jie, and I work on LHCb data analysis in high-energy physics. I’d like to share a minimal example that demonstrates how to perform template fits with 3D histograms using two different approaches:

• CERN ROOT HistFactory (documentation) (available in both C++ and Python),
• pyhf with cabinetry.

Since pyhf doesn’t natively support 3D fits, I’ve worked around this limitation by flattening 3D histograms into 1D arrays.
This example is based on a CERN HistFactory tutorial I found online, which focuses on fitting the observable RD* in the decay. The analysis involves 4 samples in a single channel, and the example compares the two approaches, yielding consistent results!

What’s in the Example:

• Histograms: The template histograms are stored in a ROOT file (DemoHistos.root) as 3D histograms, with the data histogram named h_data.
• Approaches:
  • CERN ROOT HistFactory: This approach uses the 3D histograms directly, with no need for flattening. ROOT HistFactory naturally supports multi-dimensional histograms.
  • pyhf with cabinetry: The 3D histograms are flattened into 1D arrays, which are then read into memory. With pyhf, the flattened histograms are used for fitting, while cabinetry helps streamline the management and loading of the data.

You can find the full example, along with explanations and code, in my GitHub repository here:
https://github.com/JieWu-GitHub/ROOTHisfactory-pyhf-examples

Challenges and Suggestions from this minimal example:

• Plotting with Large Numbers of Bins (Using cabinetry): One challenge I encountered was handling plots with a large number of bins, where overlapping lines made the plots hard to interpret. To address this, I adjusted the figure width, but for cases with thousands of bins, a better solution could be to automatically split the plots into sub-plots (e.g., one plot for every 100 bins). Additionally, I made projections along each axis (under physical meanings) to simplify the comparisons.

• Normalization Factors: While using cabinetry, it would be very useful to add a normalize attribute in the Sample section. This feature would allow us to easily apply fixed normalization factors to each template. Currently, I need to load the original histograms and apply normalization manually.

Looking Ahead:

I’m currently migrating my analysis work from ROOT HistFactory to pyhf + cabinetry. I will continue refining this example and greatly appreciate any feedback or suggestions. I hope this example can be helpful to others facing similar challenges when working with multi-dimensional histograms in ROOT. The inability of pyhf to handle 3D histograms directly has been a significant hurdle for me, so I’m happy to share this minimal example and emphasise that consistent results are measured.

Looking forward to hearing your thoughts and feedback!

Best regards,
Jie

[alexander-held]

Hi @JieWu-GitHub, thanks a lot for sharing this!

While using cabinetry, it would be very useful to add a normalize attribute in the Sample section. This feature would allow us to easily apply fixed normalization factors to each template.

Do you mean adding a property to the NormFactor definition that sets the parameter to be constant by default?

[JieWu-GitHub]

Hi @alexander-held
Yes, exactly! What I mean is that if we add a normalize attribute to the Sample section, and set it to True, the code would automatically calculate the NormFactor as the sum of the datasets and fix it.
This would be particularly useful in cases where we care about the absolute yields that need to be determined from the fit on the given samples.

[alexander-held]

This is an interesting idea. Would you normalize per channel and per sample? What about systematic uncertainties which could change the normalization, would they also need to be re-normalized to preserve it?

[JieWu-GitHub]

Hi Alex, Thanks for your feedback and for raising these important points!

Regarding your questions:

1. Normalization Per Channel and Per Sample:

Yes, it should be done per channel and per sample, as the templates are created independently. The idea would be to add a fixed normalization factor per channel and per sample. The corresponding values for each sample would simply be the sum of the values in the original histogram (array).

2. Systematic Uncertainties and Re-normalization:

For systematic uncertainties, in this case, we are dealing with non-data-driven templates where we primarily care about the shape of the template, not the absolute yields. Therefore, something like normsys isn’t relevant here. The only factor that might affect the normalization would be statistical uncertainties, which could be added for uncorrelated systematics (like staterror in pyhf).
As for other types of correlated shape variations (like histosys), no modification is needed, as these variations are generated in a similar way to the original template. They are compared to the original template, so they will also be scaled by the same fixed normalization factor internally.

In conclusion, for this use case, we only need to add a fixed normalization factor for the required sample.