<jupyter_start><jupyter_text>Example: Higgs to 4 Muons In this example, we will follow the [cms open data analyses](https://github.com/cms-opendata-analyses/HiggsExample20112012). Repeat the Higgs discovery (Higgs to 4 Muons) with the Coffea language. Selections/cuts used are from this [Analyzer](https://github.com/cms-opendata-analyses/HiggsExample20112012/blob/master/HiggsDemoAnalyzer/src/HiggsDemoAnalyzerGit.cc).<jupyter_code>%matplotlib inline
import numpy as np
import awkward as ak
from coffea.nanoevents import NanoEventsFactory, BaseSchema, NanoAODSchema<jupyter_output><empty_output><jupyter_text>The original data format used in this analyses is imcompatible with Coffea. I manually extracted related informations with this [analyser](https://github.com/yihui-lai/HiggsExample20112012/blob/master/HiggsDemoAnalyzer/src/HiggsDemoAnalyzerGit.cc) and stored them as flat root ntuple (`samples/Higgs4L1file.root`).<jupyter_code>events = NanoEventsFactory.from_root(
    "../../samples/Higgs4L1file.root",
    "demo/ttree",
    entry_stop=10000,
    schemaclass=NanoAODSchema
).events()
events.fields<jupyter_output><empty_output><jupyter_text>Here we begin the selections on events<jupyter_code>#select global particle flow muons
selected_muons = events.Muon[events.Muon.isPFMuon & events.Muon.isPFIsolationValid & events.Muon.hasglobalTrack]
print("muon total number: ",ak.sum(ak.num(selected_muons)))

#good muons 
selected_muons = selected_muons[ (np.abs(np.sqrt(selected_muons.dxy**2+selected_muons.dz**2)/np.sqrt(selected_muons.dxyErr**2+selected_muons.dzErr**2))<4) & 
                                 (np.abs(selected_muons.dxy)<0.5) & 
                                 (np.abs(selected_muons.dz)<1) & 
                                 (np.abs(selected_muons.pfRelIso04_all)<0.4) &
                                 (selected_muons.pt>5) & 
                                 (np.abs(selected_muons.eta)<2.4) ]
print("good muon total number: ", ak.sum(ak.num(selected_muons)))


#ZZ to 4Muons
selected_muons = selected_muons[ak.num(selected_muons) >=4]
print("events (nmuon>=4) number: ",len(selected_muons))

#sort from high pt to low pt
selected_muons = selected_muons[ak.argsort(selected_muons.pt, axis=-1, ascending=False)]
#only need the first 4 muons
selected_muons = selected_muons[:,0:4]
# total charge need to be 0
selected_muons[ak.sum(selected_muons[:,:,"charge"],axis=-1) == 0]
print("events (charge==0) number: ", len(selected_muons))
# 6 combinations: mu1mu2, mu1mu3, mu1mu4, mu2mu3, mu2mu4, mu3mu4
# after charge selection, only 4 combinations preserved 
selected_muon_pair = ak.combinations(selected_muons, 2)
selected_muon_pair = selected_muon_pair[((selected_muon_pair["0"][:,"charge"]+selected_muon_pair["1"][:,"charge"])==0)]
# define a function used to find the combination with mass closest to real Z mass
def closest(muon_pair):
    delta = abs(91.1876 - muon_pair.mass[:])
    closest_masses = np.min(delta,axis=-1)
    is_closest = (delta == closest_masses)
    return is_closest
# get the mask for the closest z boson and another z boson
z_pair = (selected_muon_pair["0"] + selected_muon_pair["1"])
za_mask = closest(z_pair)
zb_mask = ak.from_iter(za_mask[:,[3,2,1,0]])
# best combination of muons, za is closer to real Z boson mass
zamuons = selected_muon_pair[za_mask]
zbmuons = selected_muon_pair[zb_mask]
zamuons = ak.flatten(zamuons)
zbmuons = ak.flatten(zbmuons)
print(zamuons.type)
print(zbmuons.type)
# requirement on mass of za, zb
isgoodmza = (((zamuons["0"]+zamuons["1"]).mass>40)&((zamuons["0"]+zamuons["1"]).mass<120))
isgoodmzb = (((zbmuons["0"]+zbmuons["1"]).mass>12)&((zbmuons["0"]+zbmuons["1"]).mass<120))
# another requirement on za 
isgoodpt = (zamuons["0"].pt>20)&(zamuons["1"].pt>10)
Goodza = zamuons[isgoodpt&isgoodmza&isgoodmzb]
Goodzb = zbmuons[isgoodpt&isgoodmza&isgoodmzb]
# add 4 muons 
higgs = Goodza['0']+Goodza['1'] + Goodzb['0']+Goodzb['1']
higgs = higgs[higgs.mass>70]
print("final events number:",len(higgs))
print(higgs.mass)
import matplotlib.pyplot as plt
bins = np.arange(100, 140, 1)
fig, ax = plt.subplots()
n, bins, patches = ax.hist(higgs.mass, bins=bins,)
ax.set_xlabel('Higgs mass [GeV]')
ax.set_ylabel('Number of events')
ax.set_title(r'Higgs -> zz -> 4muons')

fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Write the previous code again within processor framwork<jupyter_code>import coffea.processor as processor
from coffea import hist

def closest(muon_pair):
    delta = abs(91.1876 - muon_pair.mass[:])
    closest_masses = np.min(delta,axis=-1)
    is_closest = (delta == closest_masses)
    return is_closest

class Processor(processor.ProcessorABC):
    def __init__(self):
        # Bins and categories for the histogram are defined here. For format, see https://coffeateam.github.io/coffea/stubs/coffea.hist.hist_tools.Hist.html && https://coffeateam.github.io/coffea/stubs/coffea.hist.hist_tools.Bin.html
        dataset_axis = hist.Cat("dataset", "")
        Higgs_axis = hist.Bin("Higgs", "Higgs [GeV]", 40, 100, 140)
        
        # The accumulator keeps our data chunks together for histogramming. It also gives us cutflow, which can be used to keep track of data.
        self._accumulator = processor.dict_accumulator({
            'Higgs': hist.Hist("Counts", dataset_axis, Higgs_axis),
            'cutflow': processor.defaultdict_accumulator(int)
        })
    
    @property
    def accumulator(self):
        return self._accumulator
    
    def process(self, events):
        output = self.accumulator.identity()
        
        dataset = events.metadata["dataset"]
        selected_muons = events.Muon[events.Muon.isPFMuon & events.Muon.isPFIsolationValid & events.Muon.hasglobalTrack]
        
        output['cutflow']['all muons'] += ak.sum(ak.num(selected_muons))

        selected_muons = selected_muons[ (np.abs(np.sqrt(selected_muons.dxy**2+selected_muons.dz**2)/np.sqrt(selected_muons.dxyErr**2+selected_muons.dzErr**2))<4) & 
                                        (np.abs(selected_muons.dxy)<0.5) & 
                                        (np.abs(selected_muons.dz)<1) & 
                                        (np.abs(selected_muons.pfRelIso04_all)<0.4) &
                                        (selected_muons.pt>5) & 
                                        (np.abs(selected_muons.eta)<2.4) ]
        output['cutflow']['good muons'] += ak.sum(ak.num(selected_muons))

        selected_muons = selected_muons[ak.num(selected_muons) >=4]
        output['cutflow']['events with 4 good muons'] += len(selected_muons)

        selected_muons = selected_muons[ak.argsort(selected_muons.pt, axis=-1, ascending=False)]
        selected_muons = selected_muons[:,0:4]
        selected_muons[ak.sum(selected_muons[:,:,"charge"],axis=-1) == 0]

        selected_muon_pair = ak.combinations(selected_muons, 2)
        selected_muon_pair = selected_muon_pair[((selected_muon_pair["0"][:,"charge"]+selected_muon_pair["1"][:,"charge"])==0)]

        z_pair = (selected_muon_pair["0"] + selected_muon_pair["1"])
        za_mask = closest(z_pair)
        zb_mask = ak.from_iter(za_mask[:,[3,2,1,0]])

        zamuons = selected_muon_pair[za_mask]
        zbmuons = selected_muon_pair[zb_mask]
        zamuons = ak.flatten(zamuons)
        zbmuons = ak.flatten(zbmuons)

        isgoodmza = (((zamuons["0"]+zamuons["1"]).mass>40)&((zamuons["0"]+zamuons["1"]).mass<120))
        isgoodmzb = (((zbmuons["0"]+zbmuons["1"]).mass>12)&((zbmuons["0"]+zbmuons["1"]).mass<120))
        isgoodpt = (zamuons["0"].pt>20)&(zamuons["1"].pt>10)

        Goodza = zamuons[isgoodpt&isgoodmza&isgoodmzb]
        Goodzb = zbmuons[isgoodpt&isgoodmza&isgoodmzb]
        higgs = Goodza['0']+Goodza['1'] + Goodzb['0']+Goodzb['1']
        
        output['cutflow']['final events'] += len(higgs.mass)
        output['cutflow']['number of chunks'] += 1
        
        output['Higgs'].fill(dataset=dataset, Higgs=higgs.mass)
        return output

    def postprocess(self, accumulator):
        return accumulator
events = NanoEventsFactory.from_root(
    "../../samples/Higgs4L1file.root",
    "demo/ttree",
    entry_stop=10000,
    metadata={"dataset": "Higgs_to_4muons"},
    schemaclass=NanoAODSchema
).events()
p = Processor()
output = p.process(events)
output
hist.plot1d(output['Higgs'], overlay='dataset', fill_opts={'edgecolor': (0,0,0,0.3), 'alpha': 0.8})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14443651925820256, 'fraction_numerical': 0.026509747979077507, 'mean_word_length': 3.937206572769953, 'pattern_counts': {'":': 1, '<': 31, '<?xml version=': 0, '>': 33, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '107111', 'license': 'MIT', 'n_tokens_mistral': 3123, 'n_tokens_neox': 2961, 'n_words': 504, 'path': 'source/notebooks/higgsto4l.ipynb', 'repo_name': 'UMDCMS/CoffeaTutorial', 'size': '8412'}