from typing import List, Optional
from copy import deepcopy
import basf2
import kinfit # type: ignore
import modularAnalysis as ma # type: ignore
import variables.collections as vc # type: ignore
from stdPhotons import stdPhotons # type: ignore
from variables import variables as vm # type: ignore
import variables.utils as vu # type: ignore
import pandas as pd # type: ignore
from vibe.core.utils.misc import fancy_validation_mode_header
from vibe.core.validation_mode import ValidationModeBaseClass
from vibe.core.helper.histogram_tools import HistVariable, Histogram, HistComponent
from vibe.core.helper.root_helper import makeROOTCompatible
__all__ = [
"mumuISR",
]
[docs]
@fancy_validation_mode_header
class mumuISR(ValidationModeBaseClass):
name = "mumuISR"
latex_str = r"$e^+ e^- \rightarrow \mu^+ \mu^- \gamma$"
plotting_strategies = ["multi_comparison_binned"]
[docs]
def create_basf2_path(self, **kwargs):
main = basf2.Path()
# Load muon candidates
ma.fillParticleList(
"mu+:good",
cut="p >= 0.3 and thetaInCDCAcceptance and dr < 0.5 and abs(dz) < 2",
path=main,
)
# load photon candidates
stdPhotons(
listtype="all",
beamBackgroundMVAWeight="MC15rd",
fakePhotonMVAWeight="MC15rd",
path=main,
)
ma.cutAndCopyList(
outputListName="gamma:reco",
inputListName="gamma:all",
cut="clusterNHits >= 1.5 and abs(clusterTiming) < 400 and thetaInCDCAcceptance",
path=main,
)
# reconstruct di-muon system
ma.reconstructDecay("gamma:dimuon -> mu+:good mu-:good", "", path=main)
# reconstruct beam
ma.reconstructDecay("Upsilon(4S) -> gamma:dimuon gamma:reco", "", path=main)
# perform kinematic fit with daughters update
kinfit.fitKinematic4C("Upsilon(4S)", daughtersUpdate=True, path=main)
# reconstruct recoil gamma from now fitted muons
ma.reconstructRecoil("gamma:recoil -> mu+:good mu-:good", "", path=main)
# get the angle to recoil
vm.addAlias("angleToPMiss", "angleToClosestInList(gamma:recoil)")
ma.matchMCTruth(list_name="Upsilon(4S)", path=main)
### Variables
base_particle_vars = self.base_particle_vars
cluster_vars = self.cluster_vars
# gamma tag
daughter_vars = self.generate_variables_for_list(
"Upsilon(4S) -> [^gamma -> mu+ mu-] gamma",
["gamma_tag"],
False,
base_particle_vars + cluster_vars + vc.recoil_kinematics,
)
# gamma reco
vm.addAlias("beamBackgroundSuppressionScore", "extraInfo(beamBackgroundSuppression)")
vm.addAlias("fakePhotonSuppressionScore", "extraInfo(fakePhotonSuppression)")
vm.addAlias("eres", "formula((clusterE-mcE)/mcE)")
daughter_vars += self.generate_variables_for_list(
"Upsilon(4S) -> [gamma -> mu+ mu-] ^gamma",
["gamma_reco"],
True,
base_particle_vars
+ cluster_vars
+ vc.recoil_kinematics
+ [
"angleToPMiss",
"beamBackgroundSuppressionScore",
"fakePhotonSuppressionScore",
"eres",
],
)
# muy--superkekb-porons
daughter_vars += self.generate_variables_for_list(
"Upsilon(4S) -> [gamma -> ^mu+ ^mu-] gamma",
["mu_p", "mu_m"],
True,
base_particle_vars + cluster_vars + vc.vertex + vc.track + vc.track_hits,
)
daughter_vars += self.generate_variables_for_list(
"^Upsilon(4S) -> [gamma -> mu+ mu-] gamma",
["beam"],
True,
base_particle_vars + vc.vertex,
)
# build ROE
roe_vars = self.add_ROE(main, "Upsilon(4S)")
# write out photons to study ROE later on
vm.addAlias("beamBackgroundSuppressionScore", "extraInfo(beamBackgroundSuppression)")
vm.addAlias("fakePhotonSuppressionScore", "extraInfo(fakePhotonSuppression)")
# basic_vars += [ f"nCleanedTracks({track_cut})", f"nCleanedECLClusters({isr_cut})"]
ma.copyList("Upsilon(4S):copy", "Upsilon(4S)", path=main)
self.variables_to_validation_ntuple(
decay_str="Upsilon(4S)",
variables=daughter_vars + roe_vars,
path=main,
)
self.variables_to_validation_ntuple(
decay_str="Upsilon(4S):copy",
variables=daughter_vars + roe_vars,
path=main,
)
return main
@property
def base_particle_vars(self):
return vc.kinematics + vc.inv_mass + vc.momentum_uncertainty + vc.pid
@property
def cluster_vars(self):
return vc.cluster + [
"clusterE",
"clusterUncorrE",
"clusterThetaID",
"clusterMdstIndex",
"clusterCellID",
"clusterSecondMoment",
"clusterTrackMatch",
]
mask_dict = {
"tight_mask": (
"",
"clusterNHits > 1.5 and 0.2967 < clusterTheta < 2.6180"
" and [ "
"[clusterReg == 1 and clusterE>0.100]"
" or [clusterReg == 2 and clusterE>0.060]"
" or [clusterReg == 3 and clusterE>0.150]"
"]"
" and minC2TDist>20",
),
"loose_mask": (
"",
"[ "
"[clusterReg == 1 and clusterE>0.080] "
"or [clusterReg == 2 and clusterE>0.030] "
"or [clusterReg == 3 and clusterE>0.060] "
"]",
),
"knunu_mask": (
"",
"[ "
"[clusterReg == 1 and clusterE>0.080] "
"or [clusterReg == 2 and clusterE>0.030] "
"or [clusterReg == 3 and clusterE>0.060] "
"]"
" and [minC2TDist>50]"
" and [thetaInCDCAcceptance==1]",
),
"kstartautau_mask": (
"[dr < 2] and [abs(dz) < 4] and [pt > 0.2] and [thetaInCDCAcceptance==1]",
"["
"[clusterReg==1 and E>0.08] "
"or [clusterReg==2 and E > 0.03] "
"or [clusterReg==3 and E > 0.06]"
"]"
" and [clusterNHits > 1.5] and [abs(clusterTiming) < 200] and [0.2967 < clusterTheta < 2.6180]",
),
}
base_roe_vars = vc.roe_multiplicities + vc.roe_kinematics + vc.extra_energy
[docs]
def add_ROE(self, path: basf2.Path, particle_list: str):
"""
A convenience function to add Rest Of Event (ROE) variables to a particle list.
Additionally adds various masked versions of the ROE
Args:
- path: The BASF2 path to add the ROE variables to.
- particle_list: The name of the particle list to calculate the ROE with
"""
ma.buildRestOfEvent(particle_list, fillWithMostLikely=True, path=path)
roe_vars = vu.create_aliases(self.base_roe_vars, "{variable}", "")
for mask_name, mask_items in self.mask_dict.items():
ma.appendROEMask(
list_name=particle_list,
mask_name=mask_name,
trackSelection=mask_items[0],
eclClusterSelection=mask_items[1],
path=path,
)
roe_vars += vu.create_aliases(
[v.replace("()", f"({mask_name})") for v in self.base_roe_vars], "{variable}", ""
)
return roe_vars
[docs]
def generate_variables_for_list(
self,
particle_list: str,
names: List[str],
add_mc_matching: bool,
variables: List[str],
):
variables_to_store = deepcopy(variables)
if add_mc_matching:
variables_to_store += vc.mc_truth + vc.mc_kinematics + vc.mc_variables
return vu.create_aliases_for_selected(variables_to_store, particle_list, names)
@property
def analysis_validation_histograms(self) -> List[Histogram]:
return [
Histogram(
name="EResolution_ISRPhoton",
title="Energy Resolution of ISR Photon",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_eres"),
label=r"$(E_{\mathrm{reco}}-E_{\mathrm{MC}})/E_{\mathrm{MC}}$",
unit="",
bins=50,
scope=(-2.0, 2.0),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterPhi_ISRPhoton",
title="Cluster Phi of ISR Photon",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterPhi"),
label=r"$\phi_{\mathrm{reco}}^{\gamma_{\mathrm{ISR}}}$",
unit="rad",
bins=50,
scope=(-3.1416, 3.1416),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterTheta",
title="Cluster Theta of ISR Photon",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterTheta"),
label="clusterTheta",
unit="rad",
bins=50,
scope=(0.0, 3.1416),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterTiming",
title="Cluster Timing",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterTiming"),
label="clusterTiming",
unit="ns",
bins=50,
scope=(-100, 100),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="MinC2TDist",
title="Minimum C2T Distance",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_minC2TDist"),
label="minC2TDist",
unit="cm",
bins=50,
scope=(0, 250),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterE1E9",
title="Cluster E1/E9",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterE1E9"),
label="clusterE1E9",
unit="",
bins=50,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="Cluster9E21",
title="Cluster E9/E21",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterE9E21"),
label="cluster9E21",
unit="",
bins=50,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterPSDMVA",
title="Cluster PSD MVA",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterPulseShapeDiscriminationMVA"),
label="clusterPulseShapeDiscriminationMVA",
unit="",
bins=50,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterLAT",
title="Cluster LAT",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterLAT"),
label="clusterLAT",
unit="",
bins=50,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterE",
title="Cluster Energy",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterE"),
label="clusterE",
unit="GeV",
bins=50,
scope=(0, 2),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterNHits",
title="Cluster Number of Hits",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterNHits"),
label="clusterNHits",
unit="",
bins=30,
scope=(0, 30),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterReg",
title="Cluster Region",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterReg"),
label="clusterReg",
unit="",
bins=13,
scope=(0, 13),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterSecondMoment",
title="Cluster Second Moment",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterSecondMoment"),
label="clusterSecondMoment",
unit="",
bins=50,
scope=(0, 40),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterTrackMatch",
title="Cluster Track Match",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterTrackMatch"),
label="clusterTrackMatch",
unit="",
bins=2,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterZernikeMVA",
title="Cluster Zernike MVA",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterZernikeMVA"),
label="clusterZernikeMVA",
unit="",
bins=50,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterAbsZernikeMoment40",
title="Cluster Abs Zernike Moment 4 0",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterAbsZernikeMoment40"),
label="clusterAbsZernikeMoment40",
unit="",
bins=50,
scope=(0, 1.7),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterAbsZernikeMoment51",
title="Cluster Abs Zernike Moment 5 1",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterAbsZernikeMoment51"),
label="clusterAbsZernikeMoment51",
unit="",
bins=50,
scope=(0, 1.2),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterEoP",
title="Cluster E over P",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterEoP"),
label="clusterEoP",
unit="",
bins=50,
scope=(0, 2),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterHasFailedTiming",
title="Cluster Has Failed Timing",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterHasFailedTiming"),
label="clusterHasFailedTiming",
unit="",
bins=2,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterHasFailedErrorTiming",
title="Cluster Has Failed Error Timing",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterHasFailedErrorTiming"),
label="clusterHasFailedErrorTiming",
unit="",
bins=2,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterHasNPhotons",
title="Cluster Has N Photons",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterHasNPhotons"),
label="clusterHasNPhotons",
unit="",
bins=2,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterHasNeutralHadron",
title="Cluster Has Neutral Hadron",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterHasNeutralHadron"),
label="clusterHasNeutralHadron",
unit="",
bins=2,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="ClusterHasPulseShapeDiscrimination",
title="Cluster Has Pulse Shape Discrimination",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_clusterHasPulseShapeDiscrimination"),
label="clusterHasPulseShapeDiscrimination",
unit="",
bins=2,
scope=(0, 1),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
Histogram(
name="NECLClusterTrackMatches",
title="Number of ECL Cluster Track Matches",
particle_list="Upsilon(4S)",
hist_variable=HistVariable(
df_label=makeROOTCompatible(variable="gamma_reco_nECLClusterTrackMatches"),
label="nECLClusterTrackMatches",
unit="",
bins=10,
scope=(0, 10),
),
hist_components=[
HistComponent(
label="Signal",
additional_cut_str="gamma_reco_isSignal == 1.0",
),
],
),
]
[docs]
def get_number_of_signal_for_efficiency(self, df: pd.DataFrame, particle_list: Optional[str] = None) -> float:
return df["gamma_reco_isSignal"].sum()
