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])