Belle II Software development
caf_ecms.py
1
8
9"""
10Airflow script to perform eCMS calibration (combination of the had-B and mumu method).
11"""
12
13from prompt import CalibrationSettings, INPUT_DATA_FILTERS
14from prompt.calibrations.caf_boostvector import settings as boostvector
15
16from basf2 import create_path, register_module, get_file_metadata, B2INFO, B2WARNING
17from reconstruction import prepare_cdst_analysis
18import modularAnalysis as ma
19import vertex
20import stdCharged
21import stdPi0s
22import os
23
24
25
26settings = CalibrationSettings(
27 name="Ecms Calibrations",
28 expert_username="zlebcik",
29 subsystem="beam",
30 description=__doc__,
31 input_data_formats=["cdst"],
32 input_data_names=["hadron4S", "mumu4S", "mumuOff"],
33 input_data_filters={
34 "hadron4S": [
35 INPUT_DATA_FILTERS["Data Tag"]["btocharm_calib"],
36 INPUT_DATA_FILTERS["Run Type"]["physics"],
37 INPUT_DATA_FILTERS["Beam Energy"]["4S"],
38 INPUT_DATA_FILTERS["Data Quality Tag"]["Good Or Recoverable"],
39 INPUT_DATA_FILTERS["Magnet"]["On"]],
40 "mumu4S": [
41 INPUT_DATA_FILTERS["Data Tag"]["mumu_tight_or_highm_calib"],
42 INPUT_DATA_FILTERS["Run Type"]["physics"],
43 INPUT_DATA_FILTERS["Beam Energy"]["4S"],
44 INPUT_DATA_FILTERS["Data Quality Tag"]["Good Or Recoverable"],
45 INPUT_DATA_FILTERS["Magnet"]["On"]],
46 "mumuOff": [
47 INPUT_DATA_FILTERS["Data Tag"]["mumu_tight_or_highm_calib"],
48 INPUT_DATA_FILTERS["Run Type"]["physics"],
49 INPUT_DATA_FILTERS["Beam Energy"]["Continuum"],
50 INPUT_DATA_FILTERS['Beam Energy']['Scan'],
51 INPUT_DATA_FILTERS["Data Quality Tag"]["Good Or Recoverable"],
52 INPUT_DATA_FILTERS["Magnet"]["On"]]
53 },
54 expert_config={
55 "outerLoss": "pow(0.000010e0*rawTime, 2) + 1./nEv",
56 "innerLoss": "pow(0.000120e0*rawTime, 2) + 1./nEv",
57 "runHadB": True,
58 "eCMSmumuSpread": 5.2e-3,
59 "eCMSmumuShift": 10e-3,
60 "minPXDhits": 0},
61 depends_on=[boostvector],
62 produced_payloads=["CollisionInvariantMass"])
63
64
65
66
67def get_hadB_path(isCDST):
68 """ Selects the hadronic B decays, function returns corresponding path """
69
70 # module to be run prior the collector
71 rec_path_1 = create_path()
72 if isCDST:
73 prepare_cdst_analysis(path=rec_path_1, components=['SVD', 'CDC', 'ECL', 'KLM'])
74
75 stdCharged.stdPi(listtype='loose', path=rec_path_1)
76 stdCharged.stdK(listtype='good', path=rec_path_1)
77 stdPi0s.stdPi0s(listtype='eff40_May2020', path=rec_path_1)
78
79 ma.cutAndCopyList("pi+:my", "pi+:loose", "[abs(dz)<2.0] and [abs(dr)<0.5]", path=rec_path_1)
80 ma.cutAndCopyList("K+:my", "K+:good", "[abs(dz)<2.0] and [abs(dr)<0.5]", path=rec_path_1)
81
82 ma.cutAndCopyList("pi0:my", "pi0:eff40_May2020", "", path=rec_path_1)
83
84
87
88 DcutLoose = '1.7 < M < 2.1'
89 Dcut = '1.830 < M < 1.894'
90 # Reconstructs D0s and sets decay mode identifiers
91 ma.reconstructDecay(decayString='D0:Kpi -> K-:my pi+:my', cut=DcutLoose, dmID=1, path=rec_path_1)
92 ma.reconstructDecay(decayString='D0:Kpipi0 -> K-:my pi+:my pi0:my',
93 cut=DcutLoose, dmID=2, path=rec_path_1)
94 ma.reconstructDecay(decayString='D0:Kpipipi -> K-:my pi+:my pi-:my pi+:my',
95 cut=DcutLoose, dmID=3, path=rec_path_1)
96
97 # Performs mass constrained fit for all D0 candidates
98 vertex.kFit(list_name='D0:Kpi', conf_level=0.0, fit_type='mass', path=rec_path_1)
99 # vertex.kFit(list_name='D0:Kpipi0', conf_level=0.0, fit_type='mass', path=rec_path_1)
100 vertex.kFit(list_name='D0:Kpipipi', conf_level=0.0, fit_type='mass', path=rec_path_1)
101
102 ma.applyCuts("D0:Kpi", Dcut, path=rec_path_1)
103 ma.applyCuts("D0:Kpipi0", Dcut, path=rec_path_1)
104 ma.applyCuts("D0:Kpipipi", Dcut, path=rec_path_1)
105
106 DStarcutLoose = 'massDifference(0) < 0.16'
107
108 # Reconstructs D*-s and sets decay mode identifiers
109 ma.reconstructDecay(decayString='D*+:D0pi_Kpi -> D0:Kpi pi+:my', cut=DStarcutLoose, dmID=1, path=rec_path_1)
110 ma.reconstructDecay(decayString='D*+:D0pi_Kpipi0 -> D0:Kpipi0 pi+:my',
111 cut=DStarcutLoose, dmID=2, path=rec_path_1)
112 ma.reconstructDecay(decayString='D*+:D0pi_Kpipipi -> D0:Kpipipi pi+:my',
113 cut=DStarcutLoose, dmID=3, path=rec_path_1)
114
115 BcutLoose = '[ useCMSFrame(p) < 1.6 ] and [abs(dM) < 0.25]'
116 Bcut = '[ useCMSFrame(p) < 1.2 ] and [abs(dM) < 0.05]'
117
118 # Reconstructs the signal B0 candidates from Dstar
119 ma.reconstructDecay(decayString='B0:Dstpi_D0pi_Kpi -> D*-:D0pi_Kpi pi+:my',
120 cut=BcutLoose,
121 dmID=1, path=rec_path_1)
122 ma.reconstructDecay(decayString='B0:Dstpi_D0pi_Kpipi0 -> D*-:D0pi_Kpipi0 pi+:my',
123 cut=BcutLoose,
124 dmID=2, path=rec_path_1)
125 ma.reconstructDecay(decayString='B0:Dstpi_D0pi_Kpipipi -> D*-:D0pi_Kpipipi pi+:my',
126 cut=BcutLoose,
127 dmID=3, path=rec_path_1)
128
129 vertex.treeFit('B0:Dstpi_D0pi_Kpi', updateAllDaughters=True, ipConstraint=True, path=rec_path_1)
130 vertex.treeFit('B0:Dstpi_D0pi_Kpipi0', updateAllDaughters=True, ipConstraint=True, path=rec_path_1)
131 vertex.treeFit('B0:Dstpi_D0pi_Kpipipi', updateAllDaughters=True, ipConstraint=True, path=rec_path_1)
132
133
136
137 # Reconstructs charged D mesons and sets decay mode identifiers
138 ma.reconstructDecay(decayString='D-:Kpipi -> K+:my pi-:my pi-:my',
139 cut=DcutLoose, dmID=4, path=rec_path_1)
140
141 vertex.kFit(list_name='D-:Kpipi', conf_level=0.0, fit_type='mass', path=rec_path_1)
142 ma.applyCuts("D-:Kpipi", '1.844 < M < 1.894', path=rec_path_1)
143
144 # Reconstructs the signal B candidates
145 ma.reconstructDecay(decayString='B0:Dpi_Kpipi -> D-:Kpipi pi+:my',
146 cut=BcutLoose, dmID=4, path=rec_path_1)
147
148
151
152 # Reconstructs the signal B- candidates
153 ma.reconstructDecay(decayString='B-:D0pi_Kpi -> D0:Kpi pi-:my',
154 cut=BcutLoose,
155 dmID=5, path=rec_path_1)
156 ma.reconstructDecay(decayString='B-:D0pi_Kpipi0 -> D0:Kpipi0 pi-:my',
157 cut=BcutLoose,
158 dmID=6, path=rec_path_1)
159 ma.reconstructDecay(decayString='B-:D0pi_Kpipipi -> D0:Kpipipi pi-:my',
160 cut=BcutLoose,
161 dmID=7, path=rec_path_1)
162
163 vertex.treeFit('B-:D0pi_Kpi', updateAllDaughters=True, ipConstraint=True, path=rec_path_1)
164 vertex.treeFit('B-:D0pi_Kpipi0', updateAllDaughters=True, ipConstraint=True, path=rec_path_1)
165 vertex.treeFit('B-:D0pi_Kpipipi', updateAllDaughters=True, ipConstraint=True, path=rec_path_1)
166
167 ma.copyLists(
168 outputListName='B0:merged',
169 inputListNames=[
170 'B0:Dstpi_D0pi_Kpi',
171 'B0:Dstpi_D0pi_Kpipi0',
172 'B0:Dstpi_D0pi_Kpipipi',
173 'B0:Dpi_Kpipi'
174 ],
175 path=rec_path_1)
176
177 ma.copyLists(
178 outputListName='B-:merged',
179 inputListNames=[
180 'B-:D0pi_Kpi',
181 'B-:D0pi_Kpipi0',
182 'B-:D0pi_Kpipipi',
183 ],
184 path=rec_path_1)
185
186 # Builds the rest of event object, which contains all particles not used in the reconstruction of B0 candidates.
187 ma.buildRestOfEvent(target_list_name='B0:merged', path=rec_path_1)
188
189 # Calculates the continuum suppression variables
190 cleanMask = ('cleanMask', 'nCDCHits > 0 and useCMSFrame(p)<=3.2', 'p >= 0.05 and useCMSFrame(p)<=3.2')
191 ma.appendROEMasks(list_name='B0:merged', mask_tuples=[cleanMask], path=rec_path_1)
192 ma.buildContinuumSuppression(list_name='B0:merged', roe_mask='cleanMask', path=rec_path_1)
193
194 # Builds the rest of event object, which contains all particles not used in the reconstruction of B- candidates.
195 ma.buildRestOfEvent(target_list_name='B-:merged', path=rec_path_1)
196
197 # Calculates the continuum suppression variables
198 cleanMask = ('cleanMask', 'nCDCHits > 0 and useCMSFrame(p)<=3.2', 'p >= 0.05 and useCMSFrame(p)<=3.2')
199 ma.appendROEMasks(list_name='B-:merged', mask_tuples=[cleanMask], path=rec_path_1)
200 ma.buildContinuumSuppression(list_name='B-:merged', roe_mask='cleanMask', path=rec_path_1)
201
202 ma.applyCuts("B0:merged", "[R2 < 0.3] and " + Bcut, path=rec_path_1)
203 ma.applyCuts("B-:merged", "[R2 < 0.3] and " + Bcut, path=rec_path_1)
204
205 return rec_path_1
206
207
208def get_mumu_path(isCDST, kwargs):
209 """ Selects the ee -> mumu events, function returns corresponding path """
210
211 # module to be run prior the collector
212 rec_path_1 = create_path()
213 if isCDST:
214 prepare_cdst_analysis(path=rec_path_1, components=['SVD', 'CDC', 'ECL', 'KLM'])
215
216 minPXDhits = kwargs['expert_config']['minPXDhits']
217 muSelection = '[p>1.0]'
218 muSelection += ' and abs(dz)<2.0 and abs(dr)<0.5'
219 muSelection += f' and nPXDHits >= {minPXDhits} and nSVDHits >= 8 and nCDCHits >= 20'
220
221 ma.fillParticleList('mu+:BV', muSelection, path=rec_path_1)
222 ma.reconstructDecay('Upsilon(4S):BV -> mu+:BV mu-:BV', '9.5<M<11.5', path=rec_path_1)
223 vertex.treeFit('Upsilon(4S):BV', updateAllDaughters=True, ipConstraint=True, path=rec_path_1)
224
225 return rec_path_1
226
227
228def get_data_info(inData, kwargs):
229 """ Filter the input data and returns the IOVs """
230
231 # In this script we want to use one sources of input data.
232 # Get the input files from the input_data variable
233 file_to_iov_physics = inData
234
235 # We might have requested an enormous amount of data across a run range.
236 # There's a LOT more files than runs!
237 # Lets set some limits because this calibration doesn't need that much to run.
238 max_files_per_run = 1000000
239
240 # We filter out any more than 100 files per run. The input data files are sorted alphabetically by b2caf-prompt-run
241 # already. This procedure respects that ordering
242 from prompt.utils import filter_by_max_files_per_run
243
244 reduced_file_to_iov_physics = filter_by_max_files_per_run(file_to_iov_physics, max_files_per_run)
245 input_files_physics = list(reduced_file_to_iov_physics.keys())
246 B2INFO(f"Total number of files actually used as input = {len(input_files_physics)}")
247
248 # Get the overall IoV we our process should cover. Includes the end values that we may want to ignore since our output
249 # IoV should be open ended. We could also use this as part of the input data selection in some way.
250 requested_iov = kwargs.get("requested_iov", None)
251
252 from caf.utils import IoV
253 # The actual value our output IoV payload should have. Notice that we've set it open ended.
254 output_iov = IoV(requested_iov.exp_low, requested_iov.run_low, -1, -1)
255
256 return input_files_physics, output_iov
257
258
259def is_cDST_file(fName):
260 """ Check if the file is cDST based on the metadata """
261
262 metaData = get_file_metadata(fName)
263 description = metaData.getDataDescription()
264
265 # if dataLevel is missing, determine from file name
266 if 'dataLevel' not in description:
267 B2WARNING('The cdst/mdst info is not stored in file metadata')
268 return ('cdst' in os.path.basename(fName))
269
270 return (description['dataLevel'] == 'cdst')
271
272
273def get_calibrations(input_data, **kwargs):
274 """
275 Required function used by b2caf-prompt-run tool.
276 This function return a list of Calibration objects we assign to the CAF process.
277
278 Parameters:
279 input_data (dict): Should contain every name from the 'input_data_names' variable as a key.
280 Each value is a dictionary with {"/path/to/file_e1_r5.root": IoV(1,5,1,5), ...}. Useful for
281 assigning to calibration.files_to_iov
282
283 **kwargs: Configuration options to be sent in. Since this may change we use kwargs as a way to help prevent
284 backwards compatibility problems. But you could use the correct arguments in b2caf-prompt-run for this
285 release explicitly if you want to.
286
287 Currently only kwargs["output_iov"] is used. This is the output IoV range that your payloads should
288 correspond to. Generally your highest ExpRun payload should be open ended e.g. IoV(3,4,-1,-1)
289
290 Returns:
291 list(caf.framework.Calibration): All of the calibration objects we want to assign to the CAF process
292 """
293
294 from caf.framework import Calibration
295 from caf.strategies import SingleIOV
296
297 from ROOT import Belle2 # noqa: make the Belle2 namespace available
298 from ROOT.Belle2 import InvariantMassAlgorithm
299 from caf.framework import Collection
300
301 input_files_Had, output_iov_Had = get_data_info(input_data["hadron4S"], kwargs)
302 input_files_MuMu4S, output_iov_MuMu4S = get_data_info(input_data["mumu4S"], kwargs)
303 input_files_MuMuOff, output_iov_MuMuOff = get_data_info(input_data["mumuOff"], kwargs)
304
305 # Determine if the input files are cDSTs
306 isCDST_had = is_cDST_file(input_files_Had[0]) if len(input_files_Had) > 0 else True
307 isCDST_mumu = is_cDST_file((input_files_MuMu4S + input_files_MuMuOff)[0])
308
309 rec_path_HadB = get_hadB_path(isCDST_had)
310 rec_path_MuMu = get_mumu_path(isCDST_mumu, kwargs)
311
312 collector_HadB = register_module('EcmsCollector')
313 collector_MuMu = register_module('BoostVectorCollector', Y4SPListName='Upsilon(4S):BV')
314
315 algorithm_ecms = InvariantMassAlgorithm()
316 algorithm_ecms.setOuterLoss(kwargs['expert_config']['outerLoss'])
317 algorithm_ecms.setInnerLoss(kwargs['expert_config']['innerLoss'])
318
319 algorithm_ecms.includeHadBcalib(kwargs['expert_config']['runHadB'])
320 algorithm_ecms.setMuMuEcmsSpread(kwargs['expert_config']['eCMSmumuSpread'])
321 algorithm_ecms.setMuMuEcmsOffset(kwargs['expert_config']['eCMSmumuShift'])
322
323 calibration_ecms = Calibration('eCMS', algorithms=algorithm_ecms)
324
325 collection_HadB = Collection(collector=collector_HadB,
326 input_files=input_files_Had,
327 pre_collector_path=rec_path_HadB)
328 collection_MuMu4S = Collection(collector=collector_MuMu,
329 input_files=input_files_MuMu4S,
330 pre_collector_path=rec_path_MuMu)
331 collection_MuMuOff = Collection(collector=collector_MuMu,
332 input_files=input_files_MuMuOff,
333 pre_collector_path=rec_path_MuMu)
334
335 calibration_ecms.add_collection(name='dimuon_4S', collection=collection_MuMu4S)
336 calibration_ecms.add_collection(name='dimuon_Off', collection=collection_MuMuOff)
337 calibration_ecms.add_collection(name='hadB_4S', collection=collection_HadB)
338
339 calibration_ecms.strategies = SingleIOV
340
341 # Do this for the default AlgorithmStrategy to force the output payload IoV
342 # It may be different if you are using another strategy like SequentialRunByRun
343 for algorithm in calibration_ecms.algorithms:
344 algorithm.params = {"iov_coverage": output_iov_Had}
345
346 # Most other options like database chain and backend args will be overwritten by b2caf-prompt-run.
347 # So we don't bother setting them.
348
349 # You must return all calibrations you want to run in the prompt process, even if it's only one
350 return [calibration_ecms]
351
352
stdPi(listtype=_defaultlist, path=None, writeOut=True)
stdK(listtype=_defaultlist, path=None, writeOut=True)
stdPi0s(listtype="eff60_May2020", path=None, beamBackgroundMVAWeight="", fakePhotonMVAWeight="", biasCorrectionTable="")
Definition stdPi0s.py:22
treeFit(list_name, conf_level=0.001, massConstraint=[], ipConstraint=False, updateAllDaughters=False, massConstraintDecayString='', massConstraintMassValues=[], customOriginConstraint=False, customOriginVertex=[0.001, 0, 0.0116], customOriginCovariance=[0.0048, 0, 0, 0, 0.003567, 0, 0, 0, 0.0400], originDimension=3, treatAsInvisible='', ignoreFromVertexFit='', path=None)
Definition vertex.py:237
kFit(list_name, conf_level, fit_type='vertex', constraint='', daughtersUpdate=False, decay_string='', massConstraint=[], recoilMass=0, smearing=0, path=None)
Definition vertex.py:126