diff --git a/SL_l_nu_D0toKpi_Run3_MC_data_v2/DV.py b/SL_l_nu_D0toKpi_Run3_MC_data_v2/DV.py new file mode 100644 index 0000000000000000000000000000000000000000..3caf5c0e4a3fc098d6cc0c4f5b88d7f305dac0c0 --- /dev/null +++ b/SL_l_nu_D0toKpi_Run3_MC_data_v2/DV.py @@ -0,0 +1,506 @@ +import sys,os,math +from PyConf.reading import get_pvs, get_rec_summary, get_particles +from PyConf.Algorithms import ThOrParticleSelection +import Functors as F +from FunTuple import FunctorCollection +from FunTuple import FunTuple_Particles as Funtuple +from DaVinci.algorithms import create_lines_filter +from DaVinci import Options,make_config +import FunTuple.functorcollections as FC +from DaVinciMCTools import MCTruthAndBkgCat +from Hlt2Conf.algorithms_thor import ParticleCombiner, ParticleFilter +from Hlt2Conf.algorithms_thor import ParticleContainersMerger +from Functors.grammar import Functor +from Hlt2Conf.standard_particles import make_has_rich_long_pions,make_has_rich_up_pions,make_has_rich_down_pions +import numpy as np +from .isolationMVA import mva_functor_inclusive,mva_transform_output +from .MC_Matcher import trail_seeker +from DaVinciTools import SubstitutePID +from Gaudi.Configuration import INFO + +_BPVCORRM = Functor('_BPVCORRM', 'Composite::CorrectedMass','Compute the corrected mass') +_allpv_FDVEC = (F.ENDVERTEX_POS @ F.FORWARDARG1 - F.TOLINALG @ F.POSITION @ F.FORWARDARG0) +_allpv_NORMEDDOT = F.NORMEDDOT.bind(F.THREEMOMENTUM @ F.FORWARDARG1, _allpv_FDVEC) +ALLPV_CHI2 = lambda Vertices: F.MAP(F.CHI2) @ F.TES(Vertices) +ALLPV_CHI2DOF = lambda Vertices: F.MAP(F.CHI2DOF) @ F.TES(Vertices) +ALLPV_IPCHI2 = lambda Vertices: F.MAP(F.IPCHI2).bind(F.TES(Vertices), F.FORWARDARGS) +ALLPV_FDCHI2 = lambda Vertices: F.MAP(F.VTX_FDCHI2).bind(F.TES(Vertices), F.FORWARDARGS) +ALLPV_DIRA = lambda Vertices: F.MAP(_allpv_NORMEDDOT).bind(F.TES(Vertices), F.FORWARDARGS) +ALLPV_CORRM = lambda Vertices: F.MAP(_BPVCORRM).bind(F.TES(Vertices), F.FORWARDARGS) +ALLPV_LTIME = lambda Vertices: F.MAP(F.VTX_LTIME).bind(F.TES(Vertices), F.FORWARDARGS) +ALLPV_DLS = lambda Vertices: F.MAP(F.VTX_DLS).bind(F.TES(Vertices), F.FORWARDARGS) +ALLPV_FD_COORDINATE = lambda coordinate, Vertices: F.MAP(coordinate @ _allpv_FDVEC).bind(F.TES(Vertices), F.FORWARDARGS) + +TRUE_ID_IS = lambda id, mctruth: (F.VALUE_OR(0) @ mctruth(F.ABS @ F.PARTICLE_ID) == id) + +def make_Bst(B, pions,lepton, hadron_candidate_id,v2_pvs,do_truth_matching = True): + + Bstm = SubstitutePID( + name='PIDSubstitute', + input_particles=B, + substitutions=["B-{{B*-}}","B+{{B*+}}"], + output_level=INFO).Particles + + MVA_cut = mva_transform_output(0.1) + cut_composite = F.require_all(mva_functor_inclusive(v2_pvs)[0] > MVA_cut) + + #combine to make [B*0 -> B*- pi+]cc + Bst_1 = ParticleCombiner( + Inputs=[Bstm, pions], + name=f'Bst_1_{lepton}', + DecayDescriptor='[B*0 -> B*- pi+]cc', + CombinationCut=F.ALL, + CompositeCut=cut_composite, + ParticleCombiner="ParticleVertexFitter", + PrimaryVertices=v2_pvs + #OutputLevel=1 + ) + #combine to make [B*0 -> B*- pi-]cc + Bst_2 = ParticleCombiner( + Inputs=[Bstm, pions], + name=f'Bst_2_{lepton}', + DecayDescriptor='[B*0 -> B*- pi-]cc', + CombinationCut=F.ALL, + CompositeCut=cut_composite, + ParticleCombiner="ParticleVertexFitter", + PrimaryVertices=v2_pvs + #OutputLevel=1 + ) + #merge the two Bst candidates + Bst = ParticleContainersMerger([Bst_1, Bst_2], name = f'Bst_combiner_{lepton}') + return Bst + +def get_functors(MCTRUTH, v2_pvs, rec_summary): + #define common variables for all fields (composite and basic) + vars_common = FunctorCollection() + #all pvs: ip, ipchi2. The PV positions are stored in event variables + vars_common['ALLPV_IP[pv_indx]'] = F.ALLPV_IP(v2_pvs) + vars_common['ALLPV_IPCHI2[pv_indx]'] = ALLPV_IPCHI2(v2_pvs) + #best pv: x, y, z, ip, ipchi2 + vars_common['BPV_X'] = F.BPVX(v2_pvs) + vars_common['BPV_Y'] = F.BPVY(v2_pvs) + vars_common['BPV_Z'] = F.BPVZ(v2_pvs) + vars_common['BPV_IP'] = F.BPVIP(v2_pvs) + vars_common['BPV_IPCHI2'] = F.BPVIPCHI2(v2_pvs) + vars_common['BPV_CHI2'] = F.CHI2 @ F.BPV(v2_pvs) + vars_common['BPV_CHI2DOF'] = F.CHI2DOF @ F.BPV(v2_pvs) + vars_common['CHI2DOF'] = F.VALUE_OR(Value=F.NaN) @ F.CHI2DOF + #particle id, key, truekey. The trueid is stored in MCHierarchy + vars_common['ID'] = F.PARTICLE_ID + vars_common['KEY'] = F.OBJECT_KEY + + #get charge, min ip and min ipchi2 + vars_common['CHARGE'] = F.CHARGE + vars_common['MINIP'] = F.MINIP(v2_pvs) + vars_common['MINIPCHI2'] = F.MINIPCHI2(v2_pvs) + #reco kinematics + vars_common += FC.Kinematics() + vars_common['ETA'] = F.ETA + vars_common['PHI'] = F.PHI + if MCTRUTH != 'None': + #mc vars + vars_common['TRUEKEY'] = F.VALUE_OR(-1) @ MCTRUTH(F.OBJECT_KEY) + vars_common += FC.MCKinematics(mctruth_alg = MCTRUTH) + vars_common['TRUEETA'] = MCTRUTH(F.ETA) + vars_common['TRUEPHI'] = MCTRUTH(F.PHI) + #type of the origin vertex + vars_common['MC_VTX_TYPE'] = MCTRUTH(F.VALUE_OR(-1) @ F.MC_VTX_TYPE @ F.MC_ORIGINVERTEX) + #make some helper functions + MCMOTHER_ID = lambda n: F.VALUE_OR(0) @ MCTRUTH(F.MC_MOTHER(n, F.PARTICLE_ID)) + MCMOTHER_KEY = lambda n: F.VALUE_OR(-1) @ MCTRUTH(F.MC_MOTHER(n, F.OBJECT_KEY )) + vars_common["TRUEID"] = F.VALUE_OR(0) @ MCTRUTH(F.PARTICLE_ID) + vars_common["MCKEY"] = F.VALUE_OR(0) @ MCTRUTH(F.OBJECT_KEY) + vars_common["MC_MOTHER_ID"] = MCMOTHER_ID(1) + vars_common["MC_MOTHER_KEY"] = MCMOTHER_KEY(1) + for i in range(2,14): + prefix = "MC_GD_MOTHER_{}".format(i) + vars_common[f"{prefix}_ID"] = MCMOTHER_ID(i) + vars_common[f"{prefix}_KEY"] = MCMOTHER_KEY(i) + + #variables for composite particles + vars_composite = FunctorCollection() + #end vertex position and end vertex chi2 + vars_composite['ENDVERTEX_POS_'] = F.ENDVERTEX_POS + vars_composite['ENDVERTEX_CHI2'] = F.CHI2 @ F.ENDVERTEX + vars_composite['ENDVERTEX_CHI2DOF'] = F.CHI2DOF @ F.ENDVERTEX + #all pvs: dira, fd, fdchi2 + vars_composite['ALLPV_DIRA[pv_indx]'] = ALLPV_DIRA(v2_pvs) + vars_composite['ALLPV_FD_X[pv_indx]'] = ALLPV_FD_COORDINATE(F.X_COORDINATE, v2_pvs) + vars_composite['ALLPV_FD_Y[pv_indx]'] = ALLPV_FD_COORDINATE(F.Y_COORDINATE, v2_pvs) + vars_composite['ALLPV_FD_Z[pv_indx]'] = ALLPV_FD_COORDINATE(F.Z_COORDINATE, v2_pvs) + vars_composite['ALLPV_FDCHI2[pv_indx]'] = ALLPV_FDCHI2(v2_pvs) + vars_composite['ALLPV_CORRM[pv_indx]'] = ALLPV_CORRM(v2_pvs) + vars_composite['ALLPV_LTIME[pv_indx]'] = ALLPV_LTIME(v2_pvs) + vars_composite['ALLPV_DLS[pv_indx]'] = ALLPV_DLS(v2_pvs) + #best pv: dira, fd, fdchi2, corrm, ltime, dls + vars_composite['BPV_DIRA'] = F.BPVDIRA(v2_pvs) + vars_composite['BPV_FD_'] = F.BPVFDVEC(v2_pvs) + vars_composite['BPV_FDCHI2'] = F.BPVFDCHI2(v2_pvs) + vars_composite['BPV_CORRM'] = F.BPVCORRM(v2_pvs) + vars_composite['BPV_LTIME'] = F.BPVLTIME(v2_pvs) + vars_composite['BPV_DLS'] = F.BPVDLS(v2_pvs) + + #mc composite vertex information + if MCTRUTH != 'None': vars_composite += FC.MCVertexInfo(mctruth_alg = MCTRUTH) + #Compute maximum DOCA and maximum DOCACHI2. Since there are + # only two daughters of B0 particles, the maximum DOCA/DOCACHI2 is + # the DOCA/DOCACHI2 see below. + vars_composite['MAX_DOCA'] = F.MAXDOCA + vars_composite['MAX_DOCACHI2'] = F.MAXDOCACHI2 + vars_composite['MAX_SDOCA'] = F.MAXSDOCA + vars_composite['MAX_SDOCACHI2'] = F.MAXSDOCACHI2 + + #variables for basics + vars_basic = FunctorCollection() + vars_basic['TRACKTYPE'] = F.VALUE_OR(-1) @ F.CAST_TO_INT @ F.TRACKTYPE @ F.TRACK + vars_basic['TRACKHISTORY'] = F.VALUE_OR(-1) @ F.CAST_TO_INT @ F.TRACKHISTORY @ F.TRACK + vars_basic['TRACKFLAG'] = F.VALUE_OR(-1) @ F.CAST_TO_INT @ F.TRACKFLAG @ F.TRACK + vars_basic['TRACKNDOF'] = F.VALUE_OR(-1) @ F.NDOF @ F.TRACK + vars_basic['TRACKCHI2'] = F.CHI2 @ F.TRACK + vars_basic['TRACKCHI2DOF'] = F.CHI2DOF @ F.TRACK + vars_basic['GHOSTPROB'] = F.GHOSTPROB + vars_basic['NHITS'] = F.VALUE_OR(-1) @F.NHITS @ F.TRACK + vars_basic['NUTHITS'] = F.VALUE_OR(-1) @F.NUTHITS @ F.TRACK + vars_basic['NVPHITS'] = F.VALUE_OR(-1) @F.NVPHITS @ F.TRACK + vars_basic['NFTHITS'] = F.VALUE_OR(-1) @ F.NFTHITS @ F.TRACK + vars_basic['TRACKHASUT'] = F.VALUE_OR(-1) @ F.TRACKHASUT @ F.TRACK + vars_basic['TRACKHASVELO'] = F.VALUE_OR(-1) @ F.TRACKHASVELO @ F.TRACK + vars_basic['PIDpi'] = F.PID_PI + vars_basic['PIDk'] = F.PID_K + vars_basic['PIDp'] = F.PID_P + vars_basic['PIDe'] = F.PID_E + vars_basic['PIDmu'] = F.PID_MU + vars_basic['PROBNNe'] = F.PROBNN_E + vars_basic['PROBNNpi'] = F.PROBNN_PI + vars_basic['PROBNNk'] = F.PROBNN_K + vars_basic['PROBNNmu'] = F.PROBNN_MU + vars_basic['PROBNNp'] = F.PROBNN_P + vars_basic['PROBNNghost'] = F.PROBNN_GHOST + if MCTRUTH != 'None': vars_basic += FC.MCVertexInfo(mctruth_alg = MCTRUTH) + + #variables for event + #Collection includes tis_tos, nPVs,nHits/clusters,BUNCHCROSSING and GPS var + evt_collections = [ + FC.EventInfo(), + FC.RecSummary(), + ] + evt_vars = FunctorCollection({}) + evt_vars['ALLPV_X[pv_indx]'] = F.ALLPVX(v2_pvs) + evt_vars['ALLPV_Y[pv_indx]'] = F.ALLPVY(v2_pvs) + evt_vars['ALLPV_Z[pv_indx]'] = F.ALLPVZ(v2_pvs) + evt_vars['ALLPV_CHI2[pv_indx]'] = ALLPV_CHI2(v2_pvs) + evt_vars['ALLPV_CHI2DOF[pv_indx]']= ALLPV_CHI2DOF(v2_pvs) + + for coll in evt_collections:evt_vars += coll + + vars_brems = FunctorCollection({}) + vars_brems.update({"HASBREM": F.HASBREM}) + #vars_brems.update({"HASBREMADDED": F.HASBREMADDED}) + vars_brems.update({"BREMENERGY": F.BREMENERGY}) + vars_brems.update({"BREMBENDCORR": F.BREMBENDCORR}) + vars_brems.update({"BREMPIDE": F.BREMPIDE}) + vars_brems.update({"ECALPIDE": F.ECALPIDE}) + vars_brems.update({"ECALPIDMU": F.ECALPIDMU}) + vars_brems.update({"HCALPIDE": F.HCALPIDE}) + vars_brems.update({"HCALPIDMU": F.HCALPIDMU}) + vars_brems.update({"ELECTRONSHOWEREOP": F.ELECTRONSHOWEREOP}) + vars_brems.update({"ELECTRONSHOWERDLL": F.ELECTRONSHOWERDLL}) + vars_brems.update({"CLUSTERMATCH": F.CLUSTERMATCH_CHI2}) + vars_brems.update({"ELECTRONMATCH": F.ELECTRONMATCH_CHI2}) + vars_brems.update({"BREMHYPOMATCH": F.BREMHYPOMATCH_CHI2}) + vars_brems.update({"ELECTRONENERGY": F.ELECTRONENERGY}) + vars_brems.update({"BREMHYPOENERGY": F.BREMHYPOENERGY}) + vars_brems.update({"BREMHYPODELTAX": F.BREMHYPODELTAX}) + #vars_brems.update({"BREMTRACKBASEDENERGY": F.BREMTRACKBASEDENERGY}) + vars_brems.update({"INBREM": F.INBREM}) + vars_brems.update({"INECAL": F.INECAL}) + vars_brems.update({"PT_WITH_BREM": F.PT_WITH_BREM}) + vars_brems.update({"P_WITH_BREM": F.P_WITH_BREM}) + + #return all functors + functors = (vars_common, vars_composite, vars_basic, evt_vars,vars_brems) + return functors +''' +def get_B_momentum(v2_pvs): + PVIS = F.THREEMOMENTUM + MVIS = F.MASS + EVIS = F.ENERGY + # Get the flight direction of the mother + Fnorm = F.BPVFDIR(v2_pvs) + # Calculate P_vis^perpendicular, as defined in Eq. 5.1 in JHEP(2017)2017:21 + PVIS_PERP = F.MAGNITUDE @ F.CROSS_PRODUCT.bind(PVIS, Fnorm) + # Calculate P_vis^parallel, as defined in Eq. 5.2 in JHEP(2017)2017:21 + PVIS_PARA = F.DOT_PRODUCT.bind(PVIS, Fnorm) + # Extract the known mass of the parent + parent_mass = F.PDG_MASS("B0") + # Calculate parameter a, as defined in Eq. 5.4 in JHEP(2017)2017:21 + a = (PVIS_PARA * (parent_mass**2 - MVIS**2 - 2 * PVIS_PERP**2)) / (2 * (PVIS_PARA**2 - EVIS**2)) + # Calculate r, as defined in Eq. 5.5 in arXiv:1611.08522 + r = EVIS**2 * (parent_mass**2 - MVIS**2 - 2 * PVIS_PERP**2)**2 / (4 * (PVIS_PARA**2 - EVIS**2)**2) + (EVIS * PVIS_PERP)**2 / (PVIS_PARA**2 - EVIS**2) + + return (PVIS_PARA - a + F.SQRT @ r) * Fnorm +''' +def get_fitting_variable(v2_pvs,D0,lep,Epi): + B_mass = F.PDG_MASS("B0") + B_dis = np.array([F.END_VX-F.BPVX(v2_pvs),F.END_VY-F.BPVY(v2_pvs),F.END_VZ-F.BPVZ(v2_pvs)]) + B_dis_mag = F.SQRT @ (B_dis[0]**2+B_dis[1]**2+B_dis[2]**2) + B_Mod_P = F.ABS @ ((B_mass/F.MASS)*F.PZ/(B_dis[2]/B_dis_mag)) + B_Mod_E = F.SQRT @ (B_Mod_P**2 + B_mass**2) + + Dst_PX = F.X_COORDINATE @ ((F.FOURMOMENTUM @ D0)+(F.FOURMOMENTUM @ Epi)) + Dst_PY = F.Y_COORDINATE @ ((F.FOURMOMENTUM @ D0)+(F.FOURMOMENTUM @ Epi)) + Dst_PZ = F.Z_COORDINATE @ ((F.FOURMOMENTUM @ D0)+(F.FOURMOMENTUM @ Epi)) + Dst_E = F.E_COORDINATE @ ((F.FOURMOMENTUM @ D0)+(F.FOURMOMENTUM @ Epi)) + + B_4P = np.array([B_Mod_P*B_dis[0]/B_dis_mag,B_Mod_P*B_dis[1]/B_dis_mag,B_Mod_P*B_dis[2]/B_dis_mag,B_Mod_E]) + D_4P = np.array([Dst_PX,Dst_PY,Dst_PZ,Dst_E]) + lep_4P = np.array([F.PX @ lep,F.PY @ lep,F.PZ @ lep,F.ENERGY @ lep]) + ''' + B_P_mike = get_B_momentum(v2_pvs) + B_Mod_P_mike = F.MAGNITUDE @ B_P_mike + B_E_mike = F.SQRT @ (F.DOT_PRODUCT.bind(B_P_mike, B_P_mike) + F.B_mass**2) + B_4P_mike = np.array([F.X_COORDINATE @ B_P_mike,F.Y_COORDINATE @ B_P_mike,F.Z_COORDINATE @ B_P_mike,B_E_mike]) + ''' + vars_fitting = FunctorCollection({}) + vars_fitting["Mod_P"] = B_Mod_P + #vars_fitting["Mod_P_mike"] = B_Mod_P_mike + vars_fitting['missing_m2'] = (B_4P[3]-D_4P[3]-lep_4P[3])**2-(B_4P[0]-D_4P[0]-lep_4P[0])**2-(B_4P[1]-D_4P[1]-lep_4P[1])**2-(B_4P[2]-D_4P[2]-lep_4P[2])**2 + #vars_fitting['missing_m2_mike'] = (B_4P_mike[3]-D_4P[3]-lep_4P[3])**2-(B_4P_mike[0]-D_4P[0]-lep_4P[0])**2-(B_4P_mike[1]-D_4P[1]-lep_4P[1])**2-(B_4P_mike[2]-D_4P[2]-lep_4P[2])**2 + vars_fitting['q2'] = (B_4P[3]-D_4P[3])**2-(B_4P[0]-D_4P[0])**2-(B_4P[1]-D_4P[1])**2-(B_4P[2]-D_4P[2])**2 + #vars_fitting['q2_mike'] = (B_4P_mike[3]-D_4P[3])**2-(B_4P_mike[0]-D_4P[0])**2-(B_4P_mike[1]-D_4P[1])**2-(B_4P_mike[2]-D_4P[2])**2 + B_Beta = B_Mod_P/B_Mod_E + B_gamma = 1/F.SQRT @ (1-B_Beta**2) + + B_Lorentz = [[B_gamma,-1*B_gamma*B_4P[0]/B_4P[3],-1*B_gamma*B_4P[1]/B_4P[3],-1*B_gamma*B_4P[2]/B_4P[3]],[-1*B_gamma*B_4P[0]/B_4P[3],1+(B_gamma-1)*((B_4P[0]/B_4P[3])/B_Beta)*((B_4P[0]/B_4P[3])/B_Beta),(B_gamma-1)*((B_4P[0]/B_4P[3])/B_Beta)*((B_4P[1]/B_4P[3])/B_Beta),(B_gamma-1)*((B_4P[0]/B_4P[3])/B_Beta)*((B_4P[2]/B_4P[3])/B_Beta)],[-1*B_gamma*B_4P[1]/B_4P[3],(B_gamma-1)*((B_4P[0]/B_4P[3])/B_Beta)*((B_4P[1]/B_4P[3])/B_Beta),(B_gamma-1)*((B_4P[1]/B_4P[3])/B_Beta)*((B_4P[1]/B_4P[3])/B_Beta),(B_gamma-1)*((B_4P[1]/B_4P[3])/B_Beta)*((B_4P[2]/B_4P[3])/B_Beta)],[-1*B_gamma*B_4P[2]/B_4P[3],(B_gamma-1)*((B_4P[0]/B_4P[3])/B_Beta)*((B_4P[2]/B_4P[3])/B_Beta),(B_gamma-1)*((B_4P[1]/B_4P[3])/B_Beta)*((B_4P[2]/B_4P[3])/B_Beta),(B_gamma-1)*((B_4P[2]/B_4P[3])/B_Beta)*((B_4P[2]/B_4P[3])/B_Beta)]] + vars_fitting['lep_E_Brest'] = (B_Lorentz[0][0]*lep_4P[3]+B_Lorentz[0][1]*lep_4P[0]+B_Lorentz[0][2]*lep_4P[1]+B_Lorentz[0][3]*lep_4P[2]) + ''' + B_Beta_mike = B_Mod_P_mike/B_E_mike + B_gamma_mike = 1/F.SQRT @ (1-B_Beta_mike**2) + + B_Lorentz = [[B_gamma_mike,-1*B_gamma_mike*B_4P_mike[0]/B_4P_mike[3],-1*B_gamma_mike*B_4P_mike[1]/B_4P_mike[3],-1*B_gamma_mike*B_4P_mike[2]/B_4P_mike[3]],[-1*B_gamma_mike*B_4P_mike[0]/B_4P_mike[3],1+(B_gamma_mike-1)*((B_4P_mike[0]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[0]/B_4P_mike[3])/B_Beta_mike),(B_gamma_mike-1)*((B_4P_mike[0]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[1]/B_4P_mike[3])/B_Beta_mike),(B_gamma_mike-1)*((B_4P_mike[0]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[2]/B_4P_mike[3])/B_Beta_mike)],[-1*B_gamma_mike*B_4P_mike[1]/B_4P_mike[3],(B_gamma_mike-1)*((B_4P_mike[0]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[1]/B_4P_mike[3])/B_Beta_mike),(B_gamma_mike-1)*((B_4P_mike[1]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[1]/B_4P_mike[3])/B_Beta_mike),(B_gamma_mike-1)*((B_4P_mike[1]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[2]/B_4P_mike[3])/B_Beta_mike)],[-1*B_gamma_mike*B_4P_mike[2]/B_4P_mike[3],(B_gamma_mike-1)*((B_4P_mike[0]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[2]/B_4P_mike[3])/B_Beta_mike),(B_gamma_mike-1)*((B_4P_mike[1]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[2]/B_4P_mike[3])/B_Beta_mike),(B_gamma_mike-1)*((B_4P_mike[2]/B_4P_mike[3])/B_Beta_mike)*((B_4P_mike[2]/B_4P_mike[3])/B_Beta_mike)]] + vars_fitting['lep_E_Brest_mike'] = (B_Lorentz[0][0]*lep_4P[3]+B_Lorentz[0][1]*lep_4P[0]+B_Lorentz[0][2]*lep_4P[1]+B_Lorentz[0][3]*lep_4P[2]) + ''' + return(vars_fitting) + +def tuple_Bst(Bst,lepton,line, MCTRUTH, v2_pvs, rec_summary): + #get functors + functors = get_functors(MCTRUTH, v2_pvs, rec_summary) + vars_common = functors[0] + vars_composite = functors[1] + vars_basic = functors[2] + evt_vars = functors[3] + vars_brems = functors[4] + + #variables for B0 field + vars_Bst = FunctorCollection() + B_child = F.CHILD(1, F.FORWARDARG0) + D0_child = F.CHILD(1, F.CHILD(1,F.FORWARDARG0)) + Lep_child = F.CHILD(1, F.CHILD(2,F.FORWARDARG0)) + Epi_child = F.CHILD(2, F.FORWARDARG0) + bpv_pi = F.BPV(v2_pvs).bind(Epi_child) + #diff in vtx chi2 with and without extra particle + vars_Bst['DELTA_ENDVERTEX_CHI2'] = (F.CHI2 @ F.ENDVERTEX) - (F.CHI2 @ F.ENDVERTEX.bind(B_child)) + #IP and IPChi2 of extra particle wrt to B0 end vertex + vars_Bst['Epi_IP_WRT_B0ENDVERTEX'] = F.IP.bind(F.TOLINALG @ F.POSITION @ F.ENDVERTEX @ F.FORWARDARG0 , Epi_child) + vars_Bst['Epi_IPCHI2_WRT_B0ENDVERTEX'] = F.IPCHI2.bind(F.ENDVERTEX @ F.FORWARDARG0 , Epi_child) + #IP and IPChi2 of extra particle wrt to Bm end vertex + vars_Bst['Epi_IP_WRT_BmENDVERTEX'] = F.IP.bind(F.TOLINALG @ F.POSITION @ F.ENDVERTEX @ B_child , Epi_child) + vars_Bst['Epi_IPCHI2_WRT_BmENDVERTEX'] = F.IPCHI2.bind(F.ENDVERTEX @ B_child , Epi_child) + #angle between extra particle and Bm + vars_Bst['Epi_COSANGLE_Bm'] = F.COSANGLE.bind(F.THREEMOMENTUM @ B_child, F.THREEMOMENTUM @ Epi_child) + #diff in fd chi2 with and without extra particle + vars_Bst['Delta_BPV_of_Epi_FDCHI2'] = F.VTX_FDCHI2.bind(bpv_pi, F.FORWARDARGS) - F.VTX_FDCHI2.bind(bpv_pi, B_child) + vars_Bst['BPV_of_Epi_FDCHI2'] = F.VTX_FDCHI2.bind(bpv_pi, F.FORWARDARGS) + #IP chi2 of extra particle wrt PV and SV of B0 + vars_Bst['Epi_IP_WRT_B0BPV'] = F.IP.bind(F.TOLINALG @ F.POSITION @ F.BPV(v2_pvs) @ F.FORWARDARG0, Epi_child) + vars_Bst['Epi_IPCHI2_WRT_B0BPV'] = F.IPCHI2.bind(F.BPV(v2_pvs) @ F.FORWARDARG0, Epi_child) + #IP chi2 of extra particle wrt PV and SV of Bm + vars_Bst['Epi_IP_WRT_BmBPV'] = F.IP.bind(F.TOLINALG @ F.POSITION @ F.BPV(v2_pvs) @ B_child, Epi_child) + vars_Bst['Epi_IPCHI2_WRT_BmBPV'] = F.IPCHI2.bind(F.BPV(v2_pvs) @ B_child, Epi_child) + #DOCA and DOCACHI2 b/w Lb and extra particle + vars_Bst['DOCA12'] = F.DOCA(1, 2) + vars_Bst['DOCA12_CHI2_12'] = F.DOCACHI2(1, 2) + #vars_Bst['SDOCA12'] = F.SDOCA(1, 2) + #vars_Bst['SDOCA_CHI2_12'] = F.SDOCACHI2(1, 2) + #DOCACHI2 of extra particle wrt to mother i.e. B0 + vars_Bst['MTDOCACHI2_1'] = F.MTDOCACHI2(1, v2_pvs) + vars_Bst['MTDOCACHI2_2'] = F.MTDOCACHI2(2, v2_pvs) + vars_Bst['Delta_M'] = F.ABS @ ((F.MASS @ ((F.FOURMOMENTUM @ D0_child)+(F.FOURMOMENTUM @ Epi_child))) - (F.MASS @ D0_child)) + vars_Bst['Dst_M'] = F.MASS @ ((F.FOURMOMENTUM @ D0_child)+(F.FOURMOMENTUM @ Epi_child)) + vars_Bst['Dst_PX'] = F.X_COORDINATE @ ((F.FOURMOMENTUM @ D0_child)+(F.FOURMOMENTUM @ Epi_child)) + vars_Bst['Dst_PY'] = F.Y_COORDINATE @ ((F.FOURMOMENTUM @ D0_child)+(F.FOURMOMENTUM @ Epi_child)) + vars_Bst['Dst_PZ'] = F.Z_COORDINATE @ ((F.FOURMOMENTUM @ D0_child)+(F.FOURMOMENTUM @ Epi_child)) + #vars_Bst['iso_mva'] = mva_functor_mu_inclusive(v2_pvs)[0] + if MCTRUTH != 'None': + vars_Bst += trail_seeker(MCTRUTH) + vars_Bst['iso_mva'] = mva_functor_inclusive(v2_pvs)[0] + + + #define fields + fit_vars = get_fitting_variable(v2_pvs,D0_child,Lep_child,Epi_child) + #TIS and TOS + Hlt1_decisions = ['Hlt1TrackMVA', + 'Hlt1TwoTrackMVA', + 'Hlt1D2KK', + 'Hlt1D2KPi', + 'Hlt1D2PiPi', + 'Hlt1KsToPiPi', + 'Hlt1TrackMuonMVA', + 'Hlt1SingleHighPtMuon', + 'Hlt1TrackElectronMVA', + 'Hlt1SingleHighPtElectron', + 'Hlt1DiElectronDisplaced', + 'Hlt1DiPhotonHighMass', + 'Hlt1Pi02GammaGamma', + 'Hlt1DiElectronHighMass_SS', + 'Hlt1DiElectronHighMass'] + variables_tistos_hlt1 = FC.HltTisTos(selection_type="Hlt1", trigger_lines=Hlt1_decisions, data=Bst) + evt_vars += FC.SelectionInfo(selection_type="Hlt1", trigger_lines=Hlt1_decisions) + + Hlt2_decisions = ['Hlt2SLB_BuToD0MuNu_D0ToKPi', + 'Hlt2SLB_BuToD0MuNu_D0ToKPi_FakeMuon', + 'Hlt2SLB_BuToD0ENu_D0ToKPi', + 'Hlt2SLB_BuToD0ENu_D0ToKPi_FakeElectron', + 'Hlt2SLB_B0ToDpTauNu_DpToKPiPi_TauToENuNu', + 'Hlt2SLB_B0ToDpTauNu_DpToKPiPi_TauToMuNuNu', + 'Hlt2SLB_BuToD0TauNu_D0ToKPi_TauToENuNu', + 'Hlt2SLB_BuToD0TauNu_D0ToKPi_TauToMuNuNu', + 'Hlt2SLB_BuToD0TauNu_D0ToKPi_FakeMuon', + 'Hlt2SLB_BuToD0TauNu_D0ToKPi_FakeElectron', + 'Hlt2SLB_myBuToD0ENu_D0ToKPi', + 'Hlt2Topo2Body', + 'Hlt2Topo3Body'] + + variables_tistos_hlt2 = FC.HltTisTos(selection_type="Hlt2", trigger_lines=Hlt2_decisions, data=Bst) + evt_vars += FC.SelectionInfo(selection_type="Hlt2", trigger_lines=Hlt2_decisions) + + + #define fields + fields = {} + fields['B0'] = f'[B*0 -> (B*- -> ([D0 -> K- pi+]CC) {lepton}-) [pi+]CC]CC' + fields['Bm'] = f'[B*0 -> ^(B*- -> ([D0 -> K- pi+]CC) {lepton}-) [pi+]CC]CC' + fields['D0'] = f'[B*0 -> (B*- -> ^([D0 -> K- pi+]CC) {lepton}-) [pi+]CC]CC' + fields['Lep'] = f'[B*0 -> (B*- -> ([D0 -> K- pi+]CC) ^{lepton}-) [pi+]CC]CC' + fields['Epi'] = f'[B*0 -> (B*- -> ([D0 -> K- pi+]CC) {lepton}-) ^[pi+]CC]CC' + fields['Km'] = f'[B*0 -> (B*- -> ([D0 -> ^K- pi+]CC) {lepton}-) [pi+]CC]CC' + fields['pip'] = f'[B*0 -> (B*- -> ([D0 -> K- ^pi+]CC) {lepton}-) [pi+]CC]CC' + + #add variables + variables = {} + variables["ALL"] = vars_common + variables_tistos_hlt1 + variables_tistos_hlt2 + variables["B0"] = vars_composite + vars_Bst + fit_vars + variables["Bm"] = vars_composite + variables["D0"] = vars_composite + variables["Lep"] = vars_basic + vars_brems + variables["Epi"] = vars_basic + variables["Km"] = vars_basic + variables["pip"] = vars_basic + + #define tuple + my_tuple = Funtuple( + name=f"Tuple", + tuple_name="DecayTree", + fields=fields, + variables=variables, + event_variables=evt_vars, + store_multiple_cand_info=True, + inputs=Bst) + + return my_tuple + +def get_extra_pions(): + """ + Note this function in DaVinci requires: + https://gitlab.cern.ch/lhcb/Moore/-/merge_requests/3203 + and it's related LHCb, DaVinci and Rec MRs + """ + long_pions = make_has_rich_long_pions() + up_pions = make_has_rich_up_pions() + down_pions = make_has_rich_down_pions() + return ParticleContainersMerger([long_pions, up_pions, down_pions], + name='Pions_combiner') + +#def main(options): +def main(options: Options,line = '', lep='',data_or_mc='mc'): + + #define filer + my_filter = create_lines_filter(name="Filter", lines=[line]) + + #get data and extra particles + Bm = get_particles(f"/Event/Spruce/{line}/Particles") + #extra_particles = get_extra_pions() + extra_particles = get_particles(f"/Event/Spruce/{line}/{line}_extra_tracks/Particles") + + #get v2_pvs and rec_summary + v2_pvs = get_pvs() + rec_summary = get_rec_summary() + hadron_candidate_id = 421 + #make B0 candidate + Bst = make_Bst(Bm, extra_particles,lep,hadron_candidate_id,v2_pvs,False) + #MVA_cut = mva_transform_output(0.1) + #b_cut = F.require_all(mva_functor_inclusive(v2_pvs)[0] > MVA_cut) #make a cut + #Bst_after_cut = ParticleFilter(Input=Bst, Cut=F.FILTER(b_cut), name='Bst_after_cut') #filter the B candidates + + if data_or_mc == 'mc': + MCTRUTH_Bst = MCTruthAndBkgCat(Bst, name='MCTRUTH_Bst') + tuple_file = tuple_Bst(Bst,lep,line, MCTRUTH_Bst, v2_pvs, rec_summary) + elif data_or_mc == 'data': + tuple_file = tuple_Bst(Bst,lep,line, 'None', v2_pvs, rec_summary) + else: + raise ValueError(f"Decay channel {decay_channel} not supported") + #define algorithms + user_algorithms = {} + user_algorithms['Alg'] = [my_filter, tuple_file] + + return make_config(options, user_algorithms) + +def test(options: Options): + return main(options=options, line='Hlt2_Test_line',lep='e',data_or_mc='data') + +def MC_SLB_BuToD0ENu_D0ToKPi(options: Options): + return main(options=options, line='SpruceSLB_BuToD0ENu_D0ToKPi',lep='e',data_or_mc='mc') + +def MC_SLB_BuToDststENu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0ENu_D0ToKPi',lep='e',data_or_mc='mc') + +def MC_SLB_BdToDststENu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0ENu_D0ToKPi',lep='e',data_or_mc='mc') + +def MC_SLB_BuToD0MuNu_D0ToKPi(options: Options): + return main(options=options, line='SpruceSLB_BuToD0MuNu_D0ToKPi',lep='mu',data_or_mc='mc') + +def MC_SLB_BuToD0TauNu_D0ToKPi_TauToENuNu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_TauToENuNu',lep='e',data_or_mc='mc') + +def MC_SLB_BuToDststTauNu_TauToENuNu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_TauToENuNu',lep='e',data_or_mc='mc') + +def MC_SLB_BdToDststTauNu_TauToENuNu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_TauToENuNu',lep='e',data_or_mc='mc') + +def MC_SLB_BuToD0TauNu_D0ToKPi_TauToMuNuNu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_TauToMuNuNu',lep='mu',data_or_mc='mc') + +def Data_SLB_BuToD0ENu_D0ToKPi(options: Options): + return main(options=options, line='SpruceSLB_BuToD0ENu_D0ToKPi',lep='e',data_or_mc='data') + +def Data_SLB_BuToD0MuNu_D0ToKPi(options: Options): + return main(options=options, line='SpruceSLB_BuToD0MuNu_D0ToKPi',lep='mu',data_or_mc='data') + +def Data_SLB_BuToD0TauNu_D0ToKPi_TauToENuNu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_TauToENuNu',lep='e',data_or_mc='data') + +def Data_SLB_BuToD0TauNu_D0ToKPi_TauToMuNuNu(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_TauToMuNuNu',lep='mu',data_or_mc='data') + +def Data_SLB_BuToD0ENu_D0ToKPi_FakeElectron(options: Options): + return main(options=options, line='SpruceSLB_BuToD0ENu_D0ToKPi_FakeElectron',lep='e',data_or_mc='data') + +def Data_SLB_BuToD0MuNu_D0ToKPi_FakeMuon(options: Options): + return main(options=options, line='SpruceSLB_BuToD0MuNu_D0ToKPi_FakeMuon',lep='mu',data_or_mc='data') + +def Data_SLB_BuToD0TauNu_D0ToKPi_FakeElectron(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_FakeElectron',lep='e',data_or_mc='data') + +def Data_SLB_BuToD0TauNu_D0ToKPi_FakeMuon(options: Options): + return main(options=options, line='SpruceSLB_BuToD0TauNu_D0ToKPi_FakeMuon',lep='mu',data_or_mc='data') diff --git a/SL_l_nu_D0toKpi_Run3_MC_data_v2/Hlt1.py b/SL_l_nu_D0toKpi_Run3_MC_data_v2/Hlt1.py new file mode 100644 index 0000000000000000000000000000000000000000..abab8e2f128402618b332df6c3da52583dc7dde6 --- /dev/null +++ b/SL_l_nu_D0toKpi_Run3_MC_data_v2/Hlt1.py @@ -0,0 +1,24 @@ +############################################################################### +# (c) Copyright 2022 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +from Moore import Options +from PyConf.application import configure_input, configure + +def main(options: Options): + """ + Configures algorithm running of HLT1 via Moore to be passed to Analysis Productions. + """ + + config = configure_input(options) + from Moore.production import hlt1 + hlt1(options, "--sequence=hlt1_pp_matching_no_ut_1000KHz", "--flagging") + config["Gaudi::IODataManager/IODataManager"].AgeLimit = 0 + + return config diff --git a/SL_l_nu_D0toKpi_Run3_MC_data_v2/Hlt2.py b/SL_l_nu_D0toKpi_Run3_MC_data_v2/Hlt2.py new file mode 100644 index 0000000000000000000000000000000000000000..031cb0529e07084328cf40987030edf458bbd170 --- /dev/null +++ b/SL_l_nu_D0toKpi_Run3_MC_data_v2/Hlt2.py @@ -0,0 +1,53 @@ +############################################################################### +# (c) Copyright 2022 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +""" +Stolen and adapted from Hlt/Hlt2Conf/options/hlt2_pp_expected_24_without_UT.py from v55r7. +""" +from Moore import Options,options, run_moore, config +from Hlt2Conf.settings.hlt2_binds import config_pp_2024 +from Moore.config import Hlt2Line +from Moore.streams import Stream, Streams, DETECTORS +from RecoConf.global_tools import stateProvider_with_simplified_geom, trackMasterExtrapolator_with_simplified_geom +from RecoConf.reconstruction_objects import reconstruction +from PyConf.Tools import TrackMasterExtrapolator, TrackMasterFitter +from Hlt2Conf.lines.semileptonic import all_lines as slb_lines + +def pass_through_line(name="Hlt2Passthrough"): + """Return a HLT2 line that performs no selection but runs and persists the reconstruction + """ + return Hlt2Line(name=name, prescale=1, algs=[], persistreco=True) + +from Hlt2Conf.lines.topological_b import * + +def _make_lines(): + ret = [] + for line_dict in [slb_lines]: + for line_name, builder in line_dict.items() : + if "D0ToKPi" in line_name and "TauToPiPiPi" not in line_name and "Bc" not in line_name: + print(line_name) + ret.append(builder()) + ret += [twobody_line(),threebody_line(),pass_through_line()] + return ret + +def make_streams(): + streams = [ + Stream(lines=_make_lines(), routing_bit=98, detectors=DETECTORS) + ] + return Streams(streams=streams) + +public_tools = [ + trackMasterExtrapolator_with_simplified_geom(), + stateProvider_with_simplified_geom(), +] + +def main(options: Options): + with reconstruction.bind( from_file = False) , config_pp_2024() : + return run_moore(options, make_streams, public_tools=public_tools) diff --git a/SL_l_nu_D0toKpi_Run3_MC_data_v2/MC_Matcher.py b/SL_l_nu_D0toKpi_Run3_MC_data_v2/MC_Matcher.py new file mode 100644 index 0000000000000000000000000000000000000000..e3b279805b24fb3c9f3bce2fe8827e98c23a5e25 --- /dev/null +++ b/SL_l_nu_D0toKpi_Run3_MC_data_v2/MC_Matcher.py @@ -0,0 +1,26 @@ +from FunTuple import FunctorCollection +import Functors as F +def trail_seeker(MCTRUTH): + mcMatch = lambda decay_descriptor: (F.HAS_VALUE @ MCTRUTH(F.FIND_MCDECAY(decay_descriptor))) + mcMatch_both = lambda decay_descriptor1,decay_descriptor2: (F.HAS_VALUE @ MCTRUTH(F.FIND_MCDECAY(decay_descriptor1))+F.HAS_VALUE @ MCTRUTH(F.FIND_MCDECAY(decay_descriptor2))) + + variables = FunctorCollection({ + 'IS_Dtaunu':mcMatch_both("[B0 -> D*(2010)- tau+ nu_tau]CC","[B~0 -> D*(2010)- tau+ nu_tau]CC"), + 'IS_Denu':mcMatch_both("[B0 -> D*(2010)- e+ nu_e]CC","[B~0 -> D*(2010)- e+ nu_e]CC"), + 'IS_Dstar_2460enu':mcMatch_both("[B0 -> D*(2640)- e+ nu_e]CC","[B~0 -> D*(2640)- e+ nu_e]CC"), + 'IS_D_2Senu':mcMatch_both("[B0 -> D(2S)- e+ nu_e]CC","[B~0 -> D(2S)- e+ nu_e]CC"), + 'IS_Dstar_2460_2enu':mcMatch_both("[B0 -> D*_2(2460)- e+ nu_e]CC","[B~0 -> D*_2(2460)- e+ nu_e]CC"), + 'IS_D_2460_1enu':mcMatch_both("[B0 -> D_1(2420)- e+ nu_e]CC","[B~0 -> D_1(2420)- e+ nu_e]CC"), + 'IS_Dstar_0taunu':mcMatch_both("[B0 -> D*_0- tau+ nu_tau]CC","[B~0 -> D*_0- tau+ nu_tau]CC"), + 'IS_D_H_1taunu':mcMatch_both("[B0 -> D_1(H)- tau+ nu_tau]CC","[B~0 -> D_1(H)- tau+ nu_tau]CC"), + 'IS_D_2420_1taunu':mcMatch_both("[B0 -> D_1(2420)- tau+ nu_tau]CC","[B~0 -> D_1(2420)- tau+ nu_tau]CC"), + 'IS_Dstar_2460_2taunu':mcMatch_both("[B0 -> D*_2(2460)- tau+ nu_tau]CC","[B~0 -> D*_2(2460)- tau+ nu_tau]CC"), + 'IS_Dstar0_2640enu':mcMatch("[B- -> D*(2640)0 e- nu_e~]CC"), + 'IS_D0_2Senu':mcMatch("[B- -> D(2S)0 e- nu_e~]CC"), + 'IS_Dstar0_2460_2enu':mcMatch("[B- -> D*_2(2460)0 e- nu_e~]CC"), + 'IS_D0_2460_1enu':mcMatch("[B- -> D_1(2420)0 e- nu_e~]CC"), + 'IS_D0_H_1taunu':mcMatch("[B- -> D_1(H)0 tau- nu_tau~]CC"), + 'IS_D0_2420_1taunu':mcMatch("[B- -> D_1(2420)0 tau- nu_tau~]CC"), + 'IS_Dstar0_2460_2taunu':mcMatch("[B- -> D*_2(2460)0 tau- nu_tau~]CC"), + }) + return(variables) diff --git a/SL_l_nu_D0toKpi_Run3_MC_data_v2/Spruce.py b/SL_l_nu_D0toKpi_Run3_MC_data_v2/Spruce.py new file mode 100644 index 0000000000000000000000000000000000000000..85806e3912553008cd740311905b1892f6bda7a6 --- /dev/null +++ b/SL_l_nu_D0toKpi_Run3_MC_data_v2/Spruce.py @@ -0,0 +1,66 @@ +############################################################################### +# (c) Copyright 2022 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +#Example follows from hlt2_with_hlt1_decision.py +from Moore import Options,options, run_moore, config +from Moore.lines import optimize_controlflow +from Moore.config import SpruceLine,register_line_builder +from Moore.streams import DETECTORS, Stream, Streams +from PyConf.application import configure_input, configure, default_raw_event +from RecoConf.global_tools import stateProvider_with_simplified_geom,trackMasterExtrapolator_with_simplified_geom +from PyConf.Algorithms import VeloRetinaClusterTrackingSIMD +from RecoConf.reconstruction_objects import reconstruction +from RecoConf.legacy_rec_hlt1_tracking import (make_VeloClusterTrackingSIMD, make_RetinaClusters) +from PyConf.reading import reconstruction as reco_spruce, upfront_reconstruction as upfront_spruce +from RecoConf.decoders import default_ft_decoding_version,default_VeloCluster_source +from Moore.persistence.hlt2_tistos import list_of_full_stream_lines +from Hlt2Conf.lines.semileptonic import all_lines as slb_lines +from Hlt2Conf.lines.topological_b import all_lines as topo_lines +from Hlt2Conf.lines.semileptonic import sprucing_lines as slb_sprucing_lines +from Moore.persistence.hlt2_tistos import list_of_full_stream_lines +def lines_for_tistos(): + """ + Nikole Comment : + This TISTOS part is only needed for exclusive Sprucing + and it saves the candidates of the FULL HLT2 lines that fired for each event. + You need to make sure its only FULL HLT2 lines in the lines_for_tistos + and they are the ones you care about ie. RD inclusive lines + """ + return [linename for line_dict in [slb_lines,topo_lines] for linename in line_dict.keys()] + +def make_full_streams(): + ''' + lines = [] + for line_dict in [slb_sprucing_lines]: + for line_name, builder in line_dict.items() : + if "D0ToKPi" in line_name and "TauToPiPiPi" not in line_name and "Bc" not in line_name: + print(line_name) + lines.append(builder()) + ''' + lines = [builder() for builder in slb_sprucing_lines.values()] + streams = [Stream(lines=lines, detectors=[])] + return Streams(streams=streams) + +def main(options: Options): + default_VeloCluster_source.global_bind(bank_type="VPRetinaCluster") + make_VeloClusterTrackingSIMD.global_bind( + algorithm=VeloRetinaClusterTrackingSIMD) + + public_tools = [ + trackMasterExtrapolator_with_simplified_geom(), + stateProvider_with_simplified_geom() + ] + with reconstruction.bind(from_file=True, spruce=True),\ + reco_spruce.bind(simulation=True),\ + upfront_spruce.bind(simulation=True),\ + list_of_full_stream_lines.bind(lines=lines_for_tistos()),\ + optimize_controlflow.bind(optimization="default"): + return run_moore(options, make_full_streams, public_tools=public_tools) + diff --git a/SL_l_nu_D0toKpi_Run3_MC_data_v2/info.yaml b/SL_l_nu_D0toKpi_Run3_MC_data_v2/info.yaml new file mode 100644 index 0000000000000000000000000000000000000000..95c0f05c4b67af9c8c4669f156c0b6a93ce050bb --- /dev/null +++ b/SL_l_nu_D0toKpi_Run3_MC_data_v2/info.yaml @@ -0,0 +1,156 @@ +defaults: + inform: + - ching-hua.li@cern.ch + wg: SL + +{%- set evttype_lines_and_leptons = [ + ('11584030','SLB_BuToD0ENu_D0ToKPi','b0_dstp_e'), + ('11574020','SLB_BuToD0MuNu_D0ToKPi','b0_dstp_muon'), + ('11584010','SLB_BuToD0TauNu_D0ToKPi_TauToENuNu','b0_dstp_taue'), + ('11574010','SLB_BuToD0TauNu_D0ToKPi_TauToMuNuNu','b0_dstp_taumuon'), + ('12685400','SLB_BuToDststENu','bp_dstst_e'), + ('12883000','SLB_BuToDststTauNu_TauToENuNu','bp_dstst_taue'), + ('11686000','SLB_BdToDststENu','b0_dststp_e'), + ('11883000','SLB_BdToDststTauNu_TauToENuNu','b0_dststp_taue'), +] %} + +{%- set polarity_variables = [ + ('MagDown','sim-20231017-vc-md100'), + ('MagUp','sim-20231017-vc-mu100'), +]%} +{%- for polarity, cond_tag in polarity_variables %} + {%- for evttype, line, dk in evttype_lines_and_leptons %} +HLT1_{{dk}}_{{ polarity }}: + application: "Moore/v55r11@x86_64_v3-el9-gcc13+detdesc-opt+g" + #application: "Moore/v55r8@x86_64_v3-el9-gcc13+detdesc-opt+g" + input: + bk_query: "/MC/Dev/Beam6800GeV-expected-2024-{{polarity}}-Nu7.6-25ns-Pythia8/Sim10c/{{evttype}}/DIGI" + dq_flags: + - OK + n_test_lfns: 1 + output: hlt1.dst + options: + entrypoint: SL_l_nu_D0toKpi_Run3_MC_data_v2.Hlt1:main + extra_options: + input_raw_format: 0.3 + conddb_tag: {{cond_tag}} + dddb_tag: dddb-20231017 + input_type: ROOT + output_type: ROOT + simulation: True + data_type: "Upgrade" + scheduler_legacy_mode: False + compression: + algorithm: ZSTD + level: 1 + max_buffer_size: 1048576 + #evt_max: 100 + +HLT2_{{dk}}_{{ polarity }}: + application: "Moore/v55r11@x86_64_v3-el9-gcc13+detdesc-opt+g" + #application: "Moore/v55r8@x86_64_v3-el9-gcc13+detdesc-opt+g" + input: + job_name: HLT1_{{dk}}_{{ polarity }} + output: hlt2.dst + options: + entrypoint: SL_l_nu_D0toKpi_Run3_MC_data_v2.Hlt2:main + extra_options: + input_raw_format: 0.5 + conddb_tag: {{cond_tag}} + dddb_tag: dddb-20231017 + input_type: ROOT + output_type: ROOT + simulation: True + data_type: "Upgrade" + scheduler_legacy_mode: False + output_manifest_file: "HLT2.tck.json" + compression: + algorithm: ZSTD + level: 1 + max_buffer_size: 1048576 + #evt_max: 1000 +SPRUCE_{{dk}}_{{ polarity }}: + application: "Moore/v55r11@x86_64_v3-el9-gcc13+detdesc-opt+g" + #application: "Moore/v55r8@x86_64_v3-el9-gcc13+detdesc-opt+g" + input: + job_name: HLT2_{{dk}}_{{ polarity }} + output: spruce.dst + options: + entrypoint: SL_l_nu_D0toKpi_Run3_MC_data_v2.Spruce:main + extra_options: + input_raw_format: 0.5 + conddb_tag: {{cond_tag}} + dddb_tag: dddb-20231017 + input_type: ROOT + output_type: ROOT + simulation: True + data_type: "Upgrade" + scheduler_legacy_mode: False + input_process: "Hlt2" + input_manifest_file: "HLT2.tck.json" + output_manifest_file: "SPRUCE.tck.json" + compression: + algorithm: ZSTD + level: 1 + max_buffer_size: 1048576 + #evt_max: 1000 +DV_{{dk}}_{{ polarity }}: + application: "DaVinci/v64r8@x86_64_v2-el9-clang16-opt" + input: + job_name: SPRUCE_{{dk}}_{{ polarity }} + output: NTuple.root + options: + entrypoint: SL_l_nu_D0toKpi_Run3_MC_data_v2.DV:MC_{{line}} + extra_options: + input_raw_format: 0.5 + conddb_tag: {{cond_tag}} + dddb_tag: dddb-20231017 + input_type: ROOT + simulation: True + data_type: "Upgrade" + geometry_version: run3/trunk + conditions_version: master + input_process: "Spruce" + input_manifest_file: "SPRUCE.tck.json" + #evt_max: 1000 + + {%- endfor %} +{%- endfor %} + + +{%- set lines_and_decays = [ + ('SLB_BuToD0ENu_D0ToKPi','b0_dstp_e'), + ('SLB_BuToD0MuNu_D0ToKPi','b0_dstp_muon'), + ('SLB_BuToD0TauNu_D0ToKPi_TauToENuNu','b0_dstp_taue'), + ('SLB_BuToD0TauNu_D0ToKPi_TauToMuNuNu','b0_dstp_taumuon'), + ('SLB_BuToD0ENu_D0ToKPi_FakeElectron','b0_dstp_e_fake'), + ('SLB_BuToD0MuNu_D0ToKPi_FakeMuon','b0_dstp_muon_fake'), + ('SLB_BuToD0TauNu_D0ToKPi_FakeElectron','b0_dstp_taue_fake'), + ('SLB_BuToD0TauNu_D0ToKPi_FakeMuon','b0_dstp_taumuon_fake'), +] %} + +{%- for spruce_line, decay in lines_and_decays %} +data_{{decay}}: + application: "DaVinci/v64r8@x86_64_v2-el9-clang16-opt" + input: + bk_query: "/LHCb/Collision24/Beam6800GeV-VeloClosed-MagDown/Real Data/Sprucing24c2/90000000/SL.DST" + dq_flags: + - UNCHECKED + - OK + keep_running: true + n_test_lfns: 1 + output: DATA.ROOT + options: + entrypoint: SL_l_nu_D0toKpi_Run3_MC_data_v2.DV:Data_{{spruce_line}} + extra_options: + input_raw_format: 0.5 + input_type: ROOT # ROOT for SprucingPass, RAW for RAW data (Hlt2 output) + simulation: False + #data_type: "Upgrade" + geometry_version: run3/2024.Q1.2-v00.00 + conditions_version: master + input_process: "Spruce" # Spruce for SprucingPass, "Hlt2" for RAW data (Hlt2 output) + input_stream: "sl" # for streamed data + #evt_max: 1000 +{%- endfor %} + diff --git a/SL_l_nu_D0toKpi_Run3_MC_data_v2/isolationMVA.py b/SL_l_nu_D0toKpi_Run3_MC_data_v2/isolationMVA.py new file mode 100755 index 0000000000000000000000000000000000000000..87b3f568e9a8dfa572133bb447fc83673eaecab1 --- /dev/null +++ b/SL_l_nu_D0toKpi_Run3_MC_data_v2/isolationMVA.py @@ -0,0 +1,358 @@ +############################################################################### +# (c) Copyright 2024 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +import Functors as F +from Hlt2Conf.algorithms_thor import ParticleCombiner +from Hlt2Conf.algorithms_thor import ParticleContainersMerger +from Hlt2Conf.algorithms_thor import ParticleFilter +from RecoConf.reconstruction_objects import make_pvs +import Functors.math as fmath +import math,re,os +from dataclasses import dataclass +from PyConf.Algorithms import ThOrParticleSelection +from SelAlgorithms.monitoring import monitor, histogram_1d +from Hlt2Conf.standard_particles import make_has_rich_long_pions +from Hlt2Conf.standard_particles import make_has_rich_up_pions +from Hlt2Conf.standard_particles import make_has_rich_down_pions +from DaVinciTools import SubstitutePID + +#Multiplying by this number will convert log_e to log_10 +#we multiply by log10(e) rather than dividing by ln(10) +#since multiplication is faster +LOG10_E = math.log10(math.e) +#toggle for printing debug messages +#set to False for proper productions +DEBUG = False + + +@dataclass +class MVAinput: + """Struct for storing MVA input variables""" + name: str # The name of the var (as per the model) + number: int # The number of the var + functor: F.grammar.FunctorBase # The functor to calculate that variable + x_range: (float, float) # axis range for input var monitoring histos + + +def mva_functor_inclusive(v2_pvs, + make_histos: (bool, str) = (False, ""), + useNumbers: bool = False): + """ + Compute the output of the MVA classifier (xGBoost) for charged track isolation. + Higher values of MVA classifier output indicate that the charged track is less isolated + and is more likely to be associated to be coming from the same decay vertex as the B0. + For more details: Checkout the presentation https://indico.cern.ch/event/1234758/#sc-1-4-ml-based-charged-isolat + + Args: + v2_pvs (list): TES location of v2 PVs + make_histos (tuple): tuple of bool (do you want the input vars to be monitored) and a ParticleContainersMerger (the B*-pion combinations) + useNumbers (bool): decision of whether to use the numbers or names of vars in the MVA functor + """ + #get the children of the two-body combination (B0-extraparticle) + Bstar_p_child = F.CHILD(1, F.FORWARDARGS) + ExtraParticle_child = F.CHILD(2, F.FORWARDARGS) + + #define the input variables for mva + mva_input_vars = [] #FunctorCollection() + #define Impact parameter chi2 of the extra particle wrt to the BPV associated to the two-body combination (B0-extraparticle) + mva_input_vars.append( + MVAinput("Epi_BPV_IPCHI2", 0, + F.BPVIPCHI2(v2_pvs).bind(ExtraParticle_child), (0, 20))) # PV + #define PT of two body combination (B0-extraparticle). Here is transformed to be less peaky + mva_input_vars.append( + MVAinput("Epi_PT", 1, + fmath.log(F.PT.bind(ExtraParticle_child)) * LOG10_E, (1, 5))) + #define opening angle between B0 and extra particle. Here is transformed to be less peaky + cos_angle = F.COSANGLE.bind(F.THREEMOMENTUM @ Bstar_p_child, + F.THREEMOMENTUM @ ExtraParticle_child) + mva_input_vars.append( + MVAinput("Sb_Epi_COSANGLE_Lb", 2, 1. - fmath.pow(1. - cos_angle, 0.2), + (0, 1))) + #define DIRA of the two body combination (B0-extraparticle) + mva_input_vars.append( + MVAinput("Sb_BPV_DIRA_TF", 8, fmath.pow(1. - F.BPVDIRA(v2_pvs), 0.2), + (0, 1.2))) + #define magnitude of PV distance i.e. distance between PV and vertex of two-body combination (B0-extraparticle) + # NB: the sign of the magnitude deterimined by the by the difference in z-coordinate of PV and vertex of two-body combination + mva_input_vars.append( + MVAinput("PVdis", 3, + (F.MAGNITUDE @ (F.ENDVERTEX_POS - F.BPV_POS(v2_pvs))) * + fmath.sign(F.END_VZ - F.BPVZ(v2_pvs)), (-50, 100))) + #* F.MAGNITUDE @ (F.ENDVERTEX_POS - F.BPV_POS(v2_pvs)) + #define magnitude of SV distance i.e. distance between two-body combination (B0-extraparticle) vertex and the vertex + # without the extra particle included. The sign of the magnitude is determined by the difference in z-coordinate of both vertices + mva_input_vars.append( + MVAinput("SVdis", 4, + (F.MAGNITUDE + @ (F.ENDVERTEX_POS.bind(Bstar_p_child) - F.ENDVERTEX_POS)) * + fmath.sign(F.END_VZ - F.END_VZ.bind(Bstar_p_child)), + (-50, 50))) + #define DeltaR i.e. the difference in radius of B0 and extra particle in the rapidity-azimuth plane + delta_eta = F.ETA.bind(Bstar_p_child) - F.ETA.bind(ExtraParticle_child) + delta_phi = F.PHI.bind(Bstar_p_child) - F.PHI.bind(ExtraParticle_child) + + delta_phi = fmath.where(delta_phi > math.pi, delta_phi - 2 * math.pi, + delta_phi) + delta_phi = fmath.where(delta_phi < -math.pi, delta_phi + 2 * math.pi, + delta_phi) + ''' + delta_phi = fmath.where(delta_phi > math.pi, delta_phi - math.pi, + delta_phi) + delta_phi = fmath.where(delta_phi < -math.pi, delta_phi + math.pi, + delta_phi) + ''' + mva_input_vars.append( + MVAinput( + "DeltaREpi", 5, + fmath.pow( + fmath.sqrt( + fmath.pow(delta_eta, 2.) + fmath.pow(delta_phi, 2.)), 0.2), + (0, 1.5))) + + mva_input_vars.append( + MVAinput("Sb_MAX_DOCACHI2", 6, + fmath.log(F.MAXSDOCACHI2) * LOG10_E, (-2.6, 7.5))) + mva_input_vars.append( + MVAinput( + "Sb_Epi_IPCHI2_WRT_LbENDVERTEX", 7, + fmath.log( + F.IPCHI2.bind(F.ENDVERTEX @ Bstar_p_child, + ExtraParticle_child)) * LOG10_E, (-2, 6))) # SV + + #define the mva classifier + mva = F.MVA( + MVAType='TMVA', + Config={ + 'XMLFile':'paramfile://data/xgboost_model_cocktail_inclusive_220724.xml', + #'XMLFile': + #f"file://{os.environ['ANALYSIS_PRODUCTIONS_BASE']}/SL_l_nu_D0toKpi_Run3_MC_data_v2/xgboost_model_cocktail_inclusive_220724.xml", + #"file:///eos/lhcb/wg/semileptonic/RDst_taue/xgboost_model_cocktail_inclusive_220724.xml", + #'XMLFile': 'paramfile://data/xgboost_model_cocktail_inclusive.xml', + 'Name':'BDT', + }, + Inputs={(f"f[{var.number}]" if useNumbers else var.name): var.functor + for var in mva_input_vars}) + input_monitoring = make_input_monitoring( + mva_input_vars, make_histos) if make_histos[0] else [] + return (mva, input_monitoring) + + +def make_input_monitoring(mva_input_vars: [MVAinput], + make_histos) -> [histogram_1d]: + input_monitoring = [] + + for var in mva_input_vars: + input_monitoring.append( + histogram_1d( + functor=var.functor, + label=var.name, + name=f"/{make_histos[1]}/iso_{var.name}", + title=f"isolation MVA input: {var.name}", + bins=100, + range=var.x_range)) + + return input_monitoring + + +def mva_transform_output(xgboost_style_input: float) -> float: + """ + We've converted our input from xgboost-style to TMVA-style but this gives a transformation on the cut variables + This function converts from an xgboost-style cut to the corresponding TMVA-style + taken from: https://github.com/jpata/mlglue/blob/master/mlglue/tree.py#L400-L409 + This takes the domain from [0,1] and transforms it from [-1,1] in a strange (nonlinear) way + """ + TMVA_style_output = -math.log(1. / xgboost_style_input - 1.) + TMVA_style_output = 2. / (1. + math.exp(-2. * TMVA_style_output)) - 1. + + return TMVA_style_output + + +def combine_iso_pions(BstarPlus, + extra_particles, + line_name, + v2_pvs, + apply_fiducial_cut: bool = False): + """ + Make two body combination of B0 and extra particles. + + Args: + B0 (list): TES location of B0 particles + extra_particles (list): TES location of extra particles in the event. + Pion id is assigned to the extra particle. + + Returns: + list: TES location of two body combination of B0 and extra_particles + """ + fiducial_cut = F.ALL + if apply_fiducial_cut: + #fiducial cut on pv distance + pv_dist_min = -3. #mm + fiducial_cut_pv_dist = (F.MAGNITUDE @ (F.ENDVERTEX_POS - F.BPV_POS( + v2_pvs))) * fmath.sign(F.END_VZ - F.BPVZ(v2_pvs)) > pv_dist_min + + #fiducial cut on the log IPchi2 w.r.t. SV + #See p11 : https://indico.cern.ch/event/1405184/contributions/5907467/attachments/2835541/4955107/preCuts_11_04_2024.pdf + log_ipchi2_wrtSV_max = 5 + Bstar_p_child = F.CHILD(1, F.FORWARDARGS) + ExtraParticle_child = F.CHILD(2, F.FORWARDARGS) + fiducial_cut_sv_ipchi2 = fmath.log( + F.IPCHI2.bind( + F.ENDVERTEX @ Bstar_p_child, + ExtraParticle_child)) * LOG10_E < log_ipchi2_wrtSV_max + + fiducial_cut = F.require_all(fiducial_cut_pv_dist, + fiducial_cut_sv_ipchi2) + + #first combiner: B0 pi+ + descriptor_1 = "[B*0 -> B*+ pi-]cc" + comb_1 = ParticleCombiner( + Inputs=[BstarPlus, extra_particles], + name=f'SbCombiner_onetrack_1_{line_name}', + DecayDescriptor=descriptor_1, + CombinationCut=F.ALL, + CompositeCut=fiducial_cut, + ParticleCombiner= + "ParticleVertexFitter", #NB: for neutrals need to use different combiner e.g. ParticleAdder + ) + #second combiner: B0 pi- + descriptor_2 = "[B*0 -> B*+ pi+]cc" + comb_2 = ParticleCombiner( + Inputs=[BstarPlus, extra_particles], + name=f'SbCombiner_onetrack_2_{line_name}', + DecayDescriptor=descriptor_2, + CombinationCut=F.ALL, + CompositeCut=fiducial_cut, + ParticleCombiner="ParticleVertexFitter", + ) + #merge two combiners + comb = ParticleContainersMerger([comb_1, comb_2], + name=f'bst_combiner_{line_name}') + return comb + + +def get_extra_pions(): + long_pions = make_has_rich_long_pions() + up_pions = make_has_rich_up_pions() + down_pions = make_has_rich_down_pions() + return ParticleContainersMerger([long_pions, up_pions, down_pions], + name='Pions_combiner') + + +def extract_parent_name(line_name: str) -> (str, str): + #First do some pruning + #e.g. SpruceSLB_BcToJpsiTauNu_JpsiToMuMu_FakeElectron -> JpsiTauNu + parent_string = line_name.split('_')[1].split('To')[0] + #remove any 'Fake' stuff + parent_string = re.sub(r'Fake.*$', '', parent_string) + + #define dicts to convert from our naming convention to + Bcand_dict = { + "Bu": "B+", + "B0": "B0", + "Bs": "B_s0", + "Bc": "B_c+", + "Lb": "Lambda_b0", + "Xib0": "Xi_b0", + "Xibminus": "Xi_b-", + "Omegab": "Omega_b-", + } + anti_Bcand_dict = { + "Bu": "B-", + "B0": "B~0", + "Bs": "B_s~0", + "Bc": "B_c-", + "Lb": "Lambda_b~0", + "Xib0": "Xi_b~0", + "Xibminus": "Xi_b~+", + "Omegab": "Omega_b~+", + } + #throw an error if the parent is not in the dict + if parent_string not in Bcand_dict: + raise ValueError( + f"ERROR: Couldn't automatically determine parent for line: {line_name}" + ) + + return Bcand_dict[parent_string], anti_Bcand_dict[parent_string] + + +def prepare_MVA( + line_name: str, + line_output: ParticleCombiner, + *, #All further options must be given as keyword arguments + mva_min: float = None, + parent: str = "", + anti_parent: str = "", + mva_used: str = "inclusive") -> (ThOrParticleSelection, [monitor]): + + new_parent_name = 'B*+' + anti_new_parent_name = 'B*-' + + if not parent and not anti_parent: + parent_name, anti_parent_name = extract_parent_name(line_name) + else: + parent_name = parent + anti_parent_name = anti_parent + + # Yes, we really do need all those brackets + substitutions = [ + f"{parent_name}{{{{{new_parent_name}}}}}", + f"{anti_parent_name}{{{{{anti_new_parent_name}}}}}" + ] + if DEBUG: print(f"{substitutions = }") + from Gaudi.Configuration import INFO, ALL + output_level = ALL if DEBUG else INFO + Bst_p = SubstitutePID( + name=f'PIDSubstitute_{line_name}', + input_particles=line_output, + substitutions=substitutions, + output_level=output_level).Particles + + v2_pvs = make_pvs() + pions = get_extra_pions() + Bst_z = combine_iso_pions( + Bst_p, pions, line_name, + v2_pvs) # vertex = B + each extra track - vertex 2 + + #Only the inclusive training is being maintained now, no other models are available + if mva_used == "inclusive": + mva_getter = mva_functor_inclusive + if mva_min is None: mva_min = 0.05 # default mu MVA cut + else: + raise ValueError( + f"ERROR: in {__file__}, mva requested must be the inclusive training, no other models are currently maintainted, this has gone wrong for line: {line_name}" + ) + MVA_functor, input_monitoring = mva_getter( + v2_pvs, make_histos=(True, line_name)) + MVA_cut = mva_transform_output(mva_min) + MVA_cut_functor = MVA_functor > MVA_cut + if DEBUG: print(f"for line {line_name : <56}\t MVA cut: {MVA_cut}") + # We'll monitor all the combinations from the MVA, not just those that pass + monitoring_histo = histogram_1d( + functor=MVA_functor, + label="MVA output", + name=f"/{line_name}/iso_MVA", + title=f"isolation MVA (cut at {MVA_cut} in sprucing)", + bins=200, + range=(-1, 1)) + + # However we'll only keep the subset that pass + Bst_z_subset = ParticleFilter( + Bst_z, F.FILTER(MVA_cut_functor), name=f"filter_isoPions_{line_name}" + ) # B + extra track combos passing the MVA cuts + + monitoring_histos = monitor( + data=Bst_z, histograms=input_monitoring + [monitoring_histo]) + + extra_track_cut = (F.IS_ABS_ID('pi+')) + code_extra_track = F.FILTER(extra_track_cut) @ F.GET_CHILDREN() + return (ThOrParticleSelection( + name=f"extra_track_selection_{line_name}", + InputParticles=Bst_z_subset, + Functor=code_extra_track).OutputSelection, [monitoring_histos])