Belle II Software light-2406-ragdoll
graFEI_ntupliser_B.py
1#!/usr/bin/env python3
2
3
10
11
12import basf2 as b2
13import modularAnalysis as ma
14from variables import variables as vm
15from stdCharged import stdK, stdPi
16
17# Necessary to run argparse
18from ROOT import PyConfig
19PyConfig.IgnoreCommandLineOptions = True # noqa
20
21import random
22import argparse
23
24from grafei import graFEI
25
26# Random seeds
27b2.set_random_seed(42)
28random.seed(42)
29
30
31def _parse_args():
32 parser = argparse.ArgumentParser()
33 parser.add_argument(
34 "-c", "--config", type=str, default=None, help="graFEI config file"
35 )
36 parser.add_argument(
37 "-w", "--weight", type=str, default=None, help="graFEI weight file"
38 )
39 parser.add_argument(
40 "-g",
41 "--globaltag",
42 type=str,
43 default="analysis_tools_light-2403-persian",
44 help="Globaltag containing graFEI model",
45 )
46 return parser.parse_args()
47
48
49if __name__ == "__main__":
50 args = _parse_args()
51
52 store_mc_truth = True
53
54 path = b2.create_path()
55 ma.inputMdst(filename=b2.find_file('mdst14.root', 'validation', False), path=path)
56
57 b2.conditions.prepend_globaltag(ma.getAnalysisGlobaltag())
58 if args.globaltag:
59 b2.conditions.prepend_globaltag(args.globaltag)
60
61 # These priors were obtained by counting truth-matched tracks in BB mixed MC
62 # It could be modified by the user if needed
63 priors = [0.068, 0.050, 0.7326, 0.1315, 0.0183, 0.00006]
64
65 cut_charged_graFEI = [
66 f"pidIsMostLikely({','.join(str(p) for p in priors)})>0",
67 "nCDCHits>20",
68 "thetaInCDCAcceptance",
69 "abs(dz)<1.0",
70 "dr<0.5",
71 "p<5",
72 "pt>0.2",
73 ]
74 cut_photons_graFEI = [
75 "beamBackgroundSuppression>0.4",
76 "fakePhotonSuppression>0.3",
77 "abs(clusterTiming)<100",
78 "abs(formula(clusterTiming/clusterErrorTiming))<2.0",
79 "[[clusterReg==1 and E>0.09] or [clusterReg==2 and E>0.09] or [clusterReg==3 and E>0.14]]",
80 ]
81 cut_KL0_graFEI = "PDG==0" # Reject KL0 for the moment
82
83 # Additional cuts from stdPhotons tight list are added
84 cut_photons_graFEI += [
85 "inCDCAcceptance",
86 "[clusterErrorTiming < 1e6 and [clusterE1E9 > 0.4 or E > 0.075]]",
87 ]
88
89 # How to deal with multiple signal side candidates
90 ranking_variable = "abs(cosThetaBetweenParticleAndNominalB)" # "random" # abs(cosThetaBetweenParticleAndNominalB)
91 # ----------------------------------------------------------------------------------------------------------
92
93 cut_charged_graFEI = " and ".join(cut_charged_graFEI)
94 cut_photons_graFEI = " and ".join(cut_photons_graFEI)
95
96 # Load MVA's for all gamma
97 ma.fillParticleList(
98 "gamma:all",
99 "",
100 path=path,
101 )
102 ma.getBeamBackgroundProbability("gamma:all", "MC15ri", path=path)
103 ma.getFakePhotonProbability("gamma:all", "MC15ri", path=path)
104
105 # -----------------------Signal side reconstruction------------------------------------
106
107 stdK("loose", path=path)
108 stdPi("loose", path=path)
109
110 ma.reconstructDecay(
111 "D-:Kpipi -> K+:loose pi-:loose pi-:loose", "1.85 < M < 1.88", path=path
112 )
113 ma.fillParticleList(
114 "mu+:tight",
115 "abs(dz)<2 and dr<0.5 and nCDCHits>20 and thetaInCDCAcceptance and pt>0.6 and E<5.5 and muonID_noSVD>0.9",
116 path=path,
117 )
118 ma.reconstructDecay(
119 "B0:sig -> D-:Kpipi mu+:tight",
120 "abs(cosThetaBetweenParticleAndNominalB)<2",
121 path=path,
122 )
123
124 # ROE
125 ma.buildRestOfEvent("B0:sig", path=path)
126
127 # Basic cuts to evaluate truth-matching on ROE
128 basicMask = (
129 "basicMask",
130 f"pidIsMostLikely({','.join(str(p) for p in priors)})>0 and nCDCHits>0 and thetaInCDCAcceptance and abs(dz)<4 and dr<2",
131 "inCDCAcceptance and clusterErrorTiming<1e6 and [clusterE1E9>0.4 or E>0.075] and [[clusterReg == 1 and E > 0.05]"
132 " or [clusterReg == 2 and E > 0.05] or [clusterReg == 3 and E > 0.075]]",
133 )
134
135 # Continuum suppression
136 csMask = (
137 "csMask",
138 "nCDCHits > 0 and useCMSFrame(p)<=3.2",
139 "p >= 0.05 and useCMSFrame(p)<=3.2",
140 )
141
142 ma.appendROEMasks(list_name="B0:sig", mask_tuples=[csMask, basicMask], path=path)
143
144 ma.buildContinuumSuppression(list_name="B0:sig", roe_mask="csMask", path=path)
145
146 ma.applyCuts("B0:sig", "cosTBz<0.9", path=path)
147
148 # Reject events with no signal candidates
149 skimfilter = b2.register_module("SkimFilter")
150 skimfilter.param("particleLists", ["B0:sig"])
151 empty_path = b2.create_path()
152 skimfilter.if_value("=0", empty_path, b2.AfterConditionPath.END)
153 path.add_module(skimfilter)
154
155 if store_mc_truth:
156 ma.matchMCTruth("B0:sig", path=path)
157
158 # perform random candidate selection
159 ma.rankByLowest(
160 "B0:sig",
161 variable=ranking_variable,
162 outputVariable="Rank",
163 numBest=1,
164 path=path,
165 )
166
167 if store_mc_truth:
168 # Fill particle list from ROE B0:sig to evaluate truth-matching condition
169 ma.fillParticleListFromROE(
170 "B0:forTM",
171 "",
172 maskName="basicMask",
173 sourceParticleListName="B0:sig",
174 path=path,
175 )
176
177 ma.matchMCTruth("B0:forTM", path=path)
178
179 ma.variablesToEventExtraInfo("B0:forTM", {"mcPDG": "Bsig_ROE_mcPDG"}, path=path)
180
181 # -----------------------------------------Tag side reconstruction----------------------------------------
182
183 graFEI_roe_mask = (
184 "graFEIROE",
185 cut_charged_graFEI,
186 cut_photons_graFEI,
187 cut_KL0_graFEI,
188 )
189
190 ma.appendROEMasks("B0:sig", [graFEI_roe_mask], path=path)
191
192 # Fill particle list from ROE B0:sig
193 ma.fillParticleListFromROE(
194 "B0:tag",
195 "",
196 maskName="graFEIROE",
197 sourceParticleListName="B0:sig",
198 path=path,
199 )
200
201 # Run graFEI
202 graFEI_vars = graFEI(
203 "B0:tag",
204 store_mc_truth=store_mc_truth,
205 cfg_path=args.config,
206 param_file=args.weight,
207 payload_config_name="graFEIConfigFile_Breco_example", # If you default payload name just remove this argument
208 payload_model_name="graFEIModelFile_Breco_example", # If you use default payload name just remove this argument
209 path=path,
210 )
211
212 # Reconstructin Upsilon(4S)
213 ma.reconstructDecay("Upsilon(4S):graFEI -> B0:tag B0:sig", "", path=path)
214
215 # ---------------- Write information to file ---------------------------
216
217 # Variables
218 default_vars = [
219 "M",
220 "Mbc",
221 "deltaE",
222 "p",
223 "px",
224 "py",
225 "pz",
226 "pt",
227 "E",
228 "phi",
229 "theta",
230 "cosTBz",
231 "cosThetaBetweenParticleAndNominalB",
232 ]
233 if store_mc_truth:
234 default_vars.extend(
235 [
236 "isSignal",
237 "isSignalAcceptMissingNeutrino",
238 "mcErrors",
239 "genMotherPDG",
240 "mcPDG",
241 ]
242 )
243
244 # Aliases
245 vm.addAlias("Bsig_Rank", "daughter(1, extraInfo(Rank))")
246
247 for var in default_vars:
248 vm.addAlias(f"Btag_{var}", f"daughter(0, {var})")
249 vm.addAlias(f"Bsig_{var}", f"daughter(1, {var})")
250 for var in graFEI_vars:
251 vm.addAlias(f"Btag_{var}", f"eventExtraInfo({var})")
252
253 graFEI_vars = [f"Btag_{var}" for var in graFEI_vars + default_vars] + \
254 ["Btag_Mbc", "Bsig_Rank"] + [f"Bsig_{var}" for var in default_vars]
255
256 vm.addAlias("Bsig_D_M", "daughter(1, daughter(0, M))")
257 vm.addAlias("Bsig_D_E", "daughter(1, daughter(0, E))")
258 vm.addAlias("Bsig_D_pt", "daughter(1, daughter(0, pt))")
259 vm.addAlias("Bsig_D_p", "daughter(1, daughter(0, p))")
260
261 graFEI_vars.extend(["Bsig_D_M", "Bsig_D_E", "Bsig_D_pt", "Bsig_D_p"])
262
263 ma.variablesToNtuple(
264 "Upsilon(4S):graFEI",
265 sorted(graFEI_vars),
266 filename="graFEI_BReco_example.root",
267 treename="tree",
268 path=path,
269 )
270
271 # Process
272 b2.process(path)
273
274 # print out the summary
275 print(b2.statistics)