Belle II Software prerelease-10-00-00a
caf_svd_dedx.py
1#!/usr/bin/env python3
2
3
10'''
11Script to perform the SVD dE/dx calibration
12'''
13from prompt import CalibrationSettings, INPUT_DATA_FILTERS
14import basf2 as b2
15
16from softwaretrigger.constants import HLT_INPUT_OBJECTS
17import modularAnalysis as ma
18import vertex as vx
19import reconstruction as re
20from reconstruction import prepare_user_cdst_analysis, DIGITS_OBJECTS
21
22settings = CalibrationSettings(
23 name="caf_svd_dedx",
24 expert_username="lisovsky",
25 description=__doc__,
26 input_data_formats=["cdst"],
27 input_data_names=["hadron_calib"],
28 input_data_filters={"hadron_calib": [INPUT_DATA_FILTERS["Data Tag"]["hadron_calib"],
29 INPUT_DATA_FILTERS["Beam Energy"]["4S"],
30 INPUT_DATA_FILTERS["Beam Energy"]["Continuum"],
31 INPUT_DATA_FILTERS["Run Type"]["physics"],
32 INPUT_DATA_FILTERS["Magnet"]["On"]]},
33
34 expert_config={
35 "isMC": False,
36 "listOfMutedCalibrations": [], # dEdxCalibration, dEdxValidation
37 "rerun_reco": False, # need to rerun reconstruction for calibration?
38 # This is needed when using release-10 calibration on older CDSTs as the dEdx calculation needs to be rerun.
39 "rerun_reco_val": False, # need to rerun reconstruction for validation?
40 # This is needed when using release-10 calibration on older CDSTs as the dEdx calculation needs to be rerun.
41 "rerun_pid_val": True, # need to rerun PID for validation?
42 # The answer is almost always "yes", but this argument is ignored if you are already rerunning the reconstruction.
43 "validation_mode": "basic", # full or basic validation; full also produces global PID performance plots
44 "MaxFilesPerRun": 10, # 15,
45 "MaxFilesPerRunValidation": 6, # be careful in MC to not exclude certain event types
46 "MinEvtsPerFile": 1, # sanity check against empty CDSTs
47 "MaxEvtsPerFile": 20000, # only if rerun the reco, to prevent jobs >10h
48 "MinEvtsPerTree": 100, # sanity check against an empty collector-output TTree
49 "NBinsP": 69, # number of momentum bins in the payload. 69 is the default historically; non-uniform binning is used.
50 "NBinsdEdx": 100, # number of dEdx bins in the payload. 100 is the default historically.
51 "dedxCutoff": 5.e6, # max dEdx value in the payload.
52 # Switch on a more precise calculation of the 1D profile for 2D
53 # histograms? Helps with the universality of dEdx:beta*gamma curves
54 # between hadrons.
55 "CustomProfile": True,
56 # Fix one of the parameters in the dEdx:beta*gamma fit, which makes the fit much more stable.
57 "FixUnstableFitParameter": True,
58 # In case of large changes in dEdx:beta*gamma trend with time, might need to set to False.
59 "NEventsToGenerate": 5e5, # how many events to generate in each momentum bin, for the new payloads?
60 "UsePionBGFunctionForEverything": False, # if the dEdx:beta*gamma fit is unstable, use the pion fit for all hadrons?
61 # (This can save us if the calibration failed due to e.g. kaon fit issue.)
62 "UseProtonBGFunctionForEverything": False, # if the dEdx:beta*gamma fit is unstable, use the proton fit for all hadrons?
63 "NumROCpoints": 175, # only for full validation: number of pionts for the ROC scan
64 "MinROCMomentum": 0., # only for full validation: in which momentum range to calculate ROC?
65 "MaxROCMomentum": 2.5, # only for full validation: in which momentum range to calculate ROC?
66 "NumEffBins": 30, # validation: how many bins for efficiency plots?
67 "MaxEffMomentum": 2.5 # validation: upper momentum boundary for efficiency plots?
68 # (No separation from SVD above 2.5 GeV or so)
69 },
70 depends_on=[])
71
72
73def create_path(rerun_reco, rerun_pid, isMC, expert_config):
74 rec_path = b2.Path()
75
76 # expert_config = kwargs.get("expert_config")
77 max_events_per_file = expert_config["MaxEvtsPerFile"]
78
79 if rerun_reco:
80 if not isMC:
81 rec_path.add_module(
82 'RootInput',
83 branchNames=HLT_INPUT_OBJECTS,
84 entrySequences=[f'0:{max_events_per_file - 1}'],
85 logLevel=b2.LogLevel.ERROR)
86 re.add_unpackers(path=rec_path)
87 else:
88 rec_path.add_module('RootInput', branchNames=list(DIGITS_OBJECTS) + [
89 'MCParticles',
90 'EventLevelTriggerTimeInfo',
91 'SoftwareTriggerResult',
92 'TRGSummary'], entrySequences=[f'0:{max_events_per_file - 1}'])
93 rec_path.add_module("Gearbox")
94 rec_path.add_module("Geometry")
95
96 re.add_reconstruction(path=rec_path, pruneTracks=False)
97 rec_path.add_module('VXDDedxPID')
98 elif rerun_pid:
99 rec_path.add_module('RootInput', entrySequences=[f'0:{max_events_per_file - 1}'])
100 prepare_user_cdst_analysis(rec_path, mc=isMC)
101 else:
102 rec_path.add_module('RootInput')
103
104 # Fill particle lists
105 ma.fillParticleList("pi+:all", "", path=rec_path)
106 ma.fillParticleList("pi+:lambda", "nCDCHits > 0", path=rec_path) # pi without track quality for reconstructing lambda
107 ma.fillParticleList("pi+:cut", "abs(dr) < 0.5 and abs(dz) < 2 and nSVDHits > 1",
108 path=rec_path) # pions for reconstructing D and Dstar
109
110 ma.fillParticleList('K-:cut', cut='abs(dr) < 0.5 and abs(dz) < 2 and nSVDHits > 1', path=rec_path) # kaon
111 ma.fillParticleList('e+:cut', cut='nSVDHits > 0', path=rec_path) # electron
112 # proton. In data, we only see background at p<0.25 GeV which motivates adding this cut.
113 ma.fillParticleList('p+:lambda', cut='nCDCHits > 0 and nSVDHits > 0 and p > 0.135', path=rec_path)
114 # ----------------------------------------------------------------------------
115 # Reconstruct D*(D0->K-pi+)pi+ and cc.
116 ma.reconstructDecay(decayString='D0:kpi -> K-:cut pi+:cut', cut='1.7 < M < 2.', path=rec_path)
117 ma.reconstructDecay(
118 decayString='D*+:myDstar -> D0:kpi pi+:all',
119 cut='1.95 < M <2.05 and massDifference(0) < 0.16',
120 path=rec_path)
121
122 # Reconstruct Lambda->p+pi- and cc.
123 ma.reconstructDecay('Lambda0:myLambda -> p+:lambda pi-:lambda', '1.1 < M < 1.3', path=rec_path)
124
125 # Reconstruct gamma->e+e- (photon conversion)
126 ma.reconstructDecay('gamma:myGamma -> e+:cut e-:cut', '0.0 < M < 0.5', path=rec_path)
127
128 # ----------------------------------------------------------------------------
129 # vertex fits
130 vx.treeFit(list_name='D*+:myDstar', conf_level=0, ipConstraint=True, updateAllDaughters=True, path=rec_path)
131 vx.treeFit(list_name='Lambda0:myLambda', conf_level=0, ipConstraint=True, updateAllDaughters=True, path=rec_path)
132 vx.treeFit(list_name='gamma:myGamma', conf_level=0, path=rec_path)
133
134 # ----------------------------------------------------------------------------
135 # Final selections on Lambda:
136 # - a tighter InvM range
137 # - a good vertex quality and a displacement requirement
138 # - a kinematic requirement p(p)>p(pi) which should always be satisfied in a true Lambda decay
139 # - a veto on the misidentified converted photons using the convertedPhotonInvariantMass tool, m(ee)>0.02 GeV
140 # - a veto on the misidentified K0S->pipi decays, vetoeing the region 0.488<m(pipi)<0.513 GeV
141
142 ma.cutAndCopyList(
143 outputListName='Lambda0:cut',
144 inputListName='Lambda0:myLambda',
145 cut=(
146 "1.10 < InvM < 1.13 and chiProb > 0.001 and distance > 1.0 and "
147 "cosAngleBetweenMomentumAndVertexVector > 0 and "
148 "formula(daughter(0,p)) > formula(daughter(1,p)) and convertedPhotonInvariantMass(0,1) > 0.02 and "
149 "[[formula((((daughter(0, px)**2+daughter(0, py)**2+daughter(0, pz)**2 + 0.13957**2)**0.5+"
150 "daughter(1, E))*((daughter(0, px)**2+daughter(0, py)**2+daughter(0, pz)**2 + 0.13957**2)**0.5+"
151 "daughter(1, E))-(daughter(0, px)+daughter(1, px))*(daughter(0, px)+daughter(1, px))-(daughter(0, py)+"
152 "daughter(1, py))*(daughter(0, py)+daughter(1, py))-(daughter(0, pz)+daughter(1, pz))*(daughter(0, pz)+"
153 "daughter(1, pz)))**0.5) < 0.488]"
154 "or [formula((((daughter(0, px)**2+daughter(0, py)**2+daughter(0, pz)**2 + 0.13957**2)**0.5+"
155 "daughter(1, E))*((daughter(0, px)**2+daughter(0, py)**2+daughter(0, pz)**2 + 0.13957**2)**0.5+"
156 "daughter(1, E))-(daughter(0, px)+daughter(1, px))*(daughter(0, px)+daughter(1, px))-(daughter(0, py)+"
157 "daughter(1, py))*(daughter(0, py)+daughter(1, py))-(daughter(0, pz)+daughter(1, pz))*(daughter(0, pz)+"
158 "daughter(1, pz)))**0.5) > 0.513]]"
159 ),
160 path=rec_path)
161
162 # ----------------------------------------------------------------------------
163 # Selections on Dstar
164 # - a tighter InvM and deltaM range
165 # - a good vertex quality
166
167 ma.cutAndCopyList(
168 outputListName='D*+:cut',
169 inputListName='D*+:myDstar',
170 cut=('massDifference(0) < 0.151 and 1.85 < daughter(0, InvM) < 1.88 and '
171 '1.95 < InvM < 2.05 and chiProb > 0.001 and '
172 'formula(daughter(0,cosAngleBetweenMomentumAndVertexVector))>-0.75'),
173 path=rec_path)
174
175 # ----------------------------------------------------------------------------
176 # Selections on gamma
177 # - a good vertex quality
178 # - the dr cut on the origin vertex coordinate of the photon conversion, that excludes the beam background
179 # - a tighter range on both the invariant mass and the convertedPhotonInvariantMass
180 # - geometric cuts on proximity of e+ and e- with convertedPhotonDelR and convertedPhotonDelZ
181
182 ma.cutAndCopyList(
183 outputListName='gamma:cut',
184 inputListName='gamma:myGamma',
185 cut=('chiProb > 0.001 and 1 < dr < 12 and InvM < 0.01'
186 ' and cosAngleBetweenMomentumAndVertexVector > 0.995'
187 ' and convertedPhotonInvariantMass(0,1) < 0.005'
188 ' and -0.05 < convertedPhotonDelR(0,1) < 0.15'
189 ' and -0.05 < convertedPhotonDelZ(0,1) < 0.05'
190 ),
191 path=rec_path)
192 return rec_path
193
194
195def get_calibrations(input_data, **kwargs):
196 """
197 Parameters:
198 input_data (dict): Should contain every name from the 'input_data_names' variable as a key.
199 Each value is a dictionary with {"/path/to/file_e1_r5.root": IoV(1,5,1,5), ...}. Useful for
200 assigning to calibration.files_to_iov
201
202 **kwargs: Configuration options to be sent in. Since this may change we use kwargs as a way to help prevent
203 backwards compatibility problems. But you could use the correct arguments in b2caf-prompt-run for this
204 release explicitly if you want to.
205
206 Currently only kwargs["requested_iov"] and kwargs["expert_config"] are used.
207
208 "requested_iov" is the IoV range of the bucket and your payloads should correspond to this range.
209 However your highest payload IoV should be open ended e.g. IoV(3,4,-1,-1)
210
211 "expert_config" is the input configuration. It takes default values from your `CalibrationSettings` but these are
212 overwritten by values from the 'expert_config' key in your input `caf_config.json` file when running ``b2caf-prompt-run``.
213
214 Returns:
215 list(caf.framework.Calibration): All of the calibration objects we want to assign to the CAF process
216 """
217 import basf2
218 # Set up config options
219
220 # In this script we want to use one sources of input data.
221 # Get the input files from the input_data variable
222 file_to_iov_hadron_calib = input_data["hadron_calib"]
223
224 expert_config = kwargs.get("expert_config")
225
226 isMC = expert_config["isMC"]
227 listOfMutedCalibrations = expert_config["listOfMutedCalibrations"]
228 rerun_reco = expert_config["rerun_reco"]
229 rerun_reco_val = expert_config["rerun_reco_val"]
230 rerun_pid_val = expert_config["rerun_pid_val"]
231 max_files_per_run = expert_config["MaxFilesPerRun"]
232 max_files_per_run_validation = expert_config["MaxFilesPerRunValidation"]
233
234 # Choose between the basic (default) or the full validation (produces more plots but depends on the global PID)
235 validation_mode = 1 if expert_config["validation_mode"] == "full" else 0
236
237 # If you are using Raw data there's a chance that input files could have zero events.
238 # This causes a B2FATAL in basf2 RootInput so the collector job will fail.
239 # Currently we don't have a good way of filtering this on the automated side, so we can check here.
240 min_events_per_file = expert_config["MinEvtsPerFile"]
241
242 from prompt.utils import filter_by_max_files_per_run
243
244 reduced_file_to_iov_hadron_calib = filter_by_max_files_per_run(file_to_iov_hadron_calib, max_files_per_run, min_events_per_file)
245 input_files_hadron_calib = list(reduced_file_to_iov_hadron_calib.keys())
246 basf2.B2INFO(f"Total number of files actually used as input for calibration = {len(input_files_hadron_calib)}")
247
248 if "dEdxValidation" not in listOfMutedCalibrations:
249 reduced_file_to_iov_hadron_validation = filter_by_max_files_per_run(
250 file_to_iov_hadron_calib, max_files_per_run_validation, min_events_per_file)
251 input_files_hadron_validation = list(reduced_file_to_iov_hadron_validation.keys())
252 basf2.B2INFO(f"Total number of files actually used as input for validation = {len(input_files_hadron_validation)}")
253 # Get the overall IoV we our process should cover. Includes the end values that we may want to ignore since our output
254 # IoV should be open ended. We could also use this as part of the input data selection in some way.
255 requested_iov = kwargs.get("requested_iov", None)
256
257 from caf.utils import IoV
258 # The actual value our output IoV payload should have. Notice that we've set it open ended.
259 output_iov = IoV(requested_iov.exp_low, requested_iov.run_low, -1, -1)
260
261
263 from ROOT import Belle2 # noqa: make the Belle2 namespace available
264 from ROOT.Belle2 import SVDdEdxCalibrationAlgorithm, SVDdEdxValidationAlgorithm
265
266 algo = SVDdEdxCalibrationAlgorithm()
267 algo.setMonitoringPlots(True)
268 algo.setNumPBins(expert_config['NBinsP'])
269 algo.setNumDEdxBins(expert_config['NBinsdEdx'])
270 algo.setDEdxCutoff(expert_config['dedxCutoff'])
271 algo.setMinEvtsPerTree(expert_config['MinEvtsPerTree'])
272 algo.setNEventsToGenerate(int(expert_config['NEventsToGenerate']))
273 algo.setUsePionBGFunctionForEverything(expert_config['UsePionBGFunctionForEverything'])
274 algo.setUseProtonBGFunctionForEverything(expert_config['UseProtonBGFunctionForEverything'])
275 algo.setCustomProfile(expert_config['CustomProfile'])
276 algo.setFixUnstableFitParameter(expert_config['FixUnstableFitParameter'])
277
278 if "dEdxValidation" not in listOfMutedCalibrations:
279 algo_val = SVDdEdxValidationAlgorithm()
280 algo_val.setMonitoringPlots(True)
281 algo_val.setMinEvtsPerTree(expert_config['MinEvtsPerTree'])
282 algo_val.setNumROCpoints(expert_config['NumROCpoints'])
283 algo_val.setMinROCMomentum(expert_config['MinROCMomentum'])
284 algo_val.setMaxROCMomentum(expert_config['MaxROCMomentum'])
285 algo_val.setNumEffBins(expert_config['NumEffBins'])
286 algo_val.setMaxEffMomentum(expert_config['MaxEffMomentum'])
287 algo_val.validationMode(validation_mode)
288
289
291
292 from caf.framework import Calibration
293
294 rec_path = create_path(rerun_reco, False, isMC, expert_config)
295 rec_path_validation = create_path(rerun_reco_val, rerun_pid_val, isMC, expert_config)
296
297 dedx_calibration = Calibration("SVDdEdxCalibration",
298 collector="SVDdEdxCollector",
299 algorithms=[algo],
300 input_files=input_files_hadron_calib,
301 pre_collector_path=rec_path)
302
303 if "dEdxValidation" not in listOfMutedCalibrations:
304 dedx_validation = Calibration("SVDdEdxValidation",
305 collector="SVDdEdxValidationCollector",
306 algorithms=[algo_val],
307 input_files=input_files_hadron_validation,
308 pre_collector_path=rec_path_validation)
309 # Do this for the default AlgorithmStrategy to force the output payload IoV
310 # It may be different if you are using another strategy like SequentialRunByRun
311 for algorithm in dedx_calibration.algorithms:
312 algorithm.params = {"apply_iov": output_iov}
313
314 # Most other options like database chain and backend args will be overwritten by b2caf-prompt-run.
315 # So we don't bother setting them.
316
317 if "dEdxValidation" not in listOfMutedCalibrations:
318 dedx_validation.depends_on(dedx_calibration)
319 # You must return all calibrations you want to run in the prompt process, even if it's only one
320 list_of_calibrations = []
321 if "dEdxCalibration" not in listOfMutedCalibrations:
322 list_of_calibrations.append(dedx_calibration)
323 if "dEdxValidation" not in listOfMutedCalibrations:
324 list_of_calibrations.append(dedx_validation)
325
326 return list_of_calibrations
327
328