diff --git a/dst_to_dee/DV.py b/dst_to_dee/DV.py
new file mode 100644
index 0000000000000000000000000000000000000000..27c3842dc99feae17d0b7622f4a412d9453a5341
--- /dev/null
+++ b/dst_to_dee/DV.py
@@ -0,0 +1,362 @@
+
+import Functors as F
+from Functors.math import in_range
+from GaudiKernel.SystemOfUnits import  MeV
+from Functors.math import log
+from Hlt2Conf.algorithms_thor import ParticleFilter
+from DaVinci import Options, make_config
+from DaVinci.algorithms import create_lines_filter
+from PyConf.reading import get_particles
+from FunTuple import FunctorCollection
+from PyConf.reading import get_particles, get_pvs
+import FunTuple.functorcollections as FC
+from FunTuple import FunTuple_Particles as Funtuple
+from DecayTreeFitter import DecayTreeFitter
+from Hlt2Conf.lines.charm.dst_to_dee_makers import (dst_BDT_functor)
+from .tupling import (
+    make_DeltaM_variable,
+    make_basic_variables,
+    make_composite_variables,
+    make_composite_variables_3body,
+    make_composite_variables_4body,
+    make_hlt2_event_variables,
+    make_basic_dtf_variables,
+    make_composite_dtf_variables,
+    make_composite_dtf_variables_3body,
+    make_composite_dtf_variables_4body,
+    Hlt1_lines,
+    photon_isolation_variables
+)
+extra_brem_variables= (FunctorCollection(
+            {
+                "BREMHYPODELTAX": F.BREMHYPODELTAX,
+                "BREMHYPOENERGY": F.BREMHYPOENERGY,
+                "BREMHYPOID": F.BREMHYPOID,
+                "BREMHYPOMATCH_CHI2": F.BREMHYPOMATCH_CHI2,
+                "BREMPIDE": F.BREMPIDE,
+                "INBREM": F.INBREM,
+                "MASS_WITH_BREM": F.MASS_WITH_BREM,
+                "PT_WITH_BREM": F.PT_WITH_BREM,
+                "P_WITH_BREM": F.P_WITH_BREM,
+                "ECALPIDE": F.ECALPIDE,
+                "HCALPIDE": F.HCALPIDE,
+                "ELECTRONSHOWEREOP": F.ELECTRONSHOWEREOP,
+                "CLUSTERMATCH_CHI2": F.CLUSTERMATCH_CHI2,
+                "ELECTRONMATCH_CHI2": F.ELECTRONMATCH_CHI2,
+                "ELECTRONENERGY": F.ELECTRONENERGY,
+                "ELECTRONID": F.ELECTRONID,
+                "HCALEOP": F.HCALEOP,
+                "CALO_NEUTRAL_SHOWER_SHAPE": F.CALO_NEUTRAL_SHOWER_SHAPE,
+                "RICH_DLL_E": F.VALUE_OR(F.NaN)@F.RICH_DLL_E,
+                "TRACKPX": F.VALUE_OR(-1) @ F.TRACK_PX,
+                "TRACKPY": F.VALUE_OR(-1) @ F.TRACK_PY,
+                "TRACKPZ": F.VALUE_OR(-1) @ F.TRACK_PZ,
+            })
+    )
+
+
+
+decay_descriptor = {
+    'dst_kpi_os' : {
+        "Dst0":   "[ D*(2007)0 ->  (D0 ->  K-  pi+)  (gamma ->  e+  e-)]CC",
+        "D0":     "[ D*(2007)0 -> ^(D0 ->  K-  pi+)  (gamma ->  e+  e-)]CC",
+        "Kminus": "[ D*(2007)0 ->  (D0 -> ^K-  pi+)  (gamma ->  e+  e-)]CC",
+        "piplus": "[ D*(2007)0 ->  (D0 ->  K- ^pi+)  (gamma ->  e+  e-)]CC",
+        "gamma":  "[ D*(2007)0 ->  (D0 ->  K-  pi+) ^(gamma ->  e+  e-)]CC",
+        "eplus":  "[ D*(2007)0 ->  (D0 ->  K-  pi+)  (gamma -> ^e+  e-)]CC",
+        "eminus": "[ D*(2007)0 ->  (D0 ->  K-  pi+)  (gamma ->  e+ ^e-)]CC",
+    },
+
+    'dst_kpi_ss' : {
+        "Dst0":   "[ D*(2007)0 ->  (D0 ->  K-  pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "D0":     "[ D*(2007)0 -> ^(D0 ->  K-  pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "Kminus": "[ D*(2007)0 ->  (D0 -> ^K-  pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "piplus": "[ D*(2007)0 ->  (D0 ->  K- ^pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "gamma":  "[ D*(2007)0 ->  (D0 ->  K-  pi+) ^([gamma ->  e+  e+]CC)]CC",
+        "eplus":  "[ D*(2007)0 ->  (D0 ->  K-  pi+)  ([gamma -> ^e+  e+]CC)]CC",
+        "eminus": "[ D*(2007)0 ->  (D0 ->  K-  pi+)  ([gamma ->  e+ ^e+]CC)]CC",
+    },
+
+    'dst_k3pi_os' :  {
+        "Dst0":    "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+)  (gamma ->  e+  e-)]CC",
+        "D0":      "[ D*(2007)0 -> ^(D0 -> K- pi- pi+ pi+)  (gamma ->  e+  e-)]CC",
+        "Kminus":  "[ D*(2007)0 ->  (D0 -> ^K- pi- pi+ pi+) (gamma ->  e+  e-)]CC",
+        "piminus": "[ D*(2007)0 ->  (D0 -> K- ^pi- pi+ pi+) (gamma ->  e+  e-)]CC",
+        "piplus1": "[ D*(2007)0 ->  (D0 -> K- pi- ^pi+ pi+) (gamma ->  e+  e-)]CC",
+        "piplus2": "[ D*(2007)0 ->  (D0 -> K- pi- pi+ ^pi+) (gamma ->  e+  e-)]CC",
+        "gamma":   "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+) ^(gamma ->  e+  e-)]CC",
+        "eplus":   "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+)  (gamma -> ^e+  e-)]CC",
+        "eminus":  "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+)  (gamma ->  e+ ^e-)]CC",
+    },
+
+    'dst_k3pi_ss' :  {
+        "Dst0":    "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+)   ([gamma ->  e+  e+]CC)]CC",
+        "D0":      "[ D*(2007)0 -> ^(D0 -> K- pi- pi+ pi+)   ([gamma ->  e+  e+]CC)]CC",
+        "Kminus":  "[ D*(2007)0 ->  (D0 -> ^K- pi- pi+ pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "piminus": "[ D*(2007)0 ->  (D0 -> K- ^pi- pi+ pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "piplus1": "[ D*(2007)0 ->  (D0 -> K- pi- ^pi+ pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "piplus2": "[ D*(2007)0 ->  (D0 -> K- pi- pi+ ^pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "gamma":   "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+)  ^([gamma ->  e+  e+]CC)]CC",
+        "eplus":   "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+)   ([gamma -> ^e+  e+]CC)]CC",
+        "eminus":  "[ D*(2007)0 ->  (D0 -> K- pi- pi+ pi+)   ([gamma ->  e+ ^e+]CC)]CC",
+    },
+
+    'dsstp_2kpi_os' : {
+        "Dstp":   "[ D*_s+ ->  (D_s+ -> K- K+ pi+)  (gamma ->  e+  e-)]CC",
+        "DpDs":   "[ D*_s+ -> ^(D_s+ -> K- K+ pi+)  (gamma ->  e+  e-)]CC",
+        "Kminus": "[ D*_s+ ->  (D_s+ -> ^K- K+ pi+) (gamma ->  e+  e-)]CC",
+        "Kplus":  "[ D*_s+ ->  (D_s+ -> K- ^K+ pi+) (gamma ->  e+  e-)]CC",
+        "piplus": "[ D*_s+ ->  (D_s+ -> K- K+ ^pi+) (gamma ->  e+  e-)]CC",
+        "gamma":  "[ D*_s+ ->  (D_s+ -> K- K+ pi+) ^(gamma ->  e+  e-)]CC",
+        "eplus":  "[ D*_s+ ->  (D_s+ -> K- K+ pi+)  (gamma -> ^e+  e-)]CC",
+        "eminus": "[ D*_s+ ->  (D_s+ -> K- K+ pi+)  (gamma ->  e+ ^e-)]CC",
+    },
+    'dsstp_2kpi_ss' : {
+        "Dstp":   "[ D*_s+ ->  (D_s+ -> K- K+ pi+)   ([gamma ->  e+  e+]CC)]CC",
+        "DpDs":   "[ D*_s+ -> ^(D_s+ -> K- K+ pi+)   ([gamma ->  e+  e+]CC)]CC",
+        "Kminus": "[ D*_s+ ->  (D_s+ -> ^K- K+ pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "Kplus":  "[ D*_s+ ->  (D_s+ -> K- ^K+ pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "piplus": "[ D*_s+ ->  (D_s+ -> K- K+ ^pi+)  ([gamma ->  e+  e+]CC)]CC",
+        "gamma":  "[ D*_s+ ->  (D_s+ -> K- K+ pi+)  ^([gamma ->  e+  e+]CC)]CC",
+        "eplus":  "[ D*_s+ ->  (D_s+ -> K- K+ pi+)   ([gamma -> ^e+  e+]CC)]CC",
+        "eminus": "[ D*_s+ ->  (D_s+ -> K- K+ pi+)   ([gamma ->  e+ ^e+]CC)]CC",
+    },
+
+}
+
+def MVA_dst_variables(line, pvs):
+    return (FunctorCollection( { "MVA": dst_BDT_functor(pvs, line=line)}))
+#make_dtf_variables(options,  input_data, ptype):
+def make_dtf_variables(options, input_data,  pvs, ptype="basic",  mass_constrain = ["D*(2007)0", "D0" ], nparticles=2):
+    
+    DTF = DecayTreeFitter(
+        name="DTF_{hash}",
+        input_particles=input_data,
+    )
+
+    DTF_DMass_BestPV = DecayTreeFitter(
+        name="DTF_DMass_BestPV_{hash}",
+        input_particles=input_data,
+        mass_constraints=mass_constrain,
+        input_pvs=pvs,
+        fit_all_pvs=False,
+    )
+    DTF_Mass_BestPV = DecayTreeFitter(
+        name="DTF_Mass_BestPV_{hash}",
+        input_particles=input_data,
+        mass_constraints=[mass_constrain[1]],
+        input_pvs=pvs,
+        fit_all_pvs=False,
+    )
+
+    DTF_BestPV = DecayTreeFitter(
+        name="DTF_BestPV_{hash}",
+        input_particles=input_data,
+        #mass_constraints=["D*(2007)0", "D0" ],
+        input_pvs=pvs,
+        fit_all_pvs=False,
+    )
+
+    DTF_DMass = DecayTreeFitter(
+        name="DTF_DMass_{hash}",
+        input_particles=input_data,
+        mass_constraints=mass_constrain,
+        #output_level=3,
+    )
+    DTF_Mass = DecayTreeFitter(
+        name="DTF_Mass_{hash}",
+        input_particles=input_data,
+        mass_constraints=[mass_constrain[1]],
+        #output_level=3,
+    )
+    
+    if ptype == "basic":
+        dtf_variables =  make_basic_dtf_variables
+    elif ptype == "composite":
+        if nparticles == 2:
+            dtf_variables = make_composite_dtf_variables 
+        elif nparticles == 3:
+            dtf_variables = make_composite_dtf_variables_3body
+        elif nparticles == 4:
+            dtf_variables = make_composite_dtf_variables_4body 
+
+    #dtf_vars = dtf_variables(options, pvs, input_data,
+#                                        DTF=DTF,
+#                                        pv_constraint=False,
+#                                        mass_constraint=False)
+    dtf_vars = dtf_variables(options, pvs, input_data,
+                                        DTF=DTF_BestPV,
+                                        pv_constraint=True,
+                                        mass_constraint=False)
+    #dtf_vars += dtf_variables(options, pvs, input_data,
+    #                                    DTF=DTF_Mass,
+    #                                    pv_constraint=False,
+    #                                    mass_constraint=True, particle_name="D")
+    #dtf_vars += dtf_variables(options, pvs, input_data,
+    #                                    DTF=DTF_Mass_BestPV,
+    #                                    pv_constraint=True,
+    #                                    mass_constraint=True, particle_name="D")
+    #dtf_vars += dtf_variables(options, pvs, input_data,
+    #                                    DTF=DTF_DMass,
+    #                                    pv_constraint=False,
+    #                                    mass_constraint=True, particle_name="DstD")
+    dtf_vars += dtf_variables(options, pvs, input_data,
+                                        DTF=DTF_DMass_BestPV,
+                                        pv_constraint=True,
+                                        mass_constraint=True, particle_name="DstD")
+
+    return dtf_vars
+
+
+CHILD_1 = lambda func: F.CHILD(1, func)
+CHILD_2 = lambda func: F.CHILD(2, func)
+
+def make_Dst0ToD0EmEp_D0ToKmPip_tuple(options, line, pvs , rec_summary, dd='dst_kpi_os'):
+
+    init_data=get_particles(f"/Event/HLT2/{line}/Particles")
+    extra_photons = get_particles(f"/Event/HLT2/{line}/{line[10:]}_prompt_Photons/Particles")
+    my_filter = create_lines_filter(f"LineFilter_{line}_{{hash}}",[line],)
+    input_data = ParticleFilter(init_data, F.FILTER(F.require_all(
+                    in_range(1900. * MeV, F.MASS, 2100. * MeV),
+                    in_range(1825. * MeV, CHILD_1(F.MASS), 1900. * MeV),
+                    F.MAXDOCA < 0.2,
+                    CHILD_1(F.MAXDOCA) < 0.2,
+                    CHILD_2(F.MAXDOCA) < 0.2,)))
+
+    fields = decay_descriptor[dd]
+    
+    extra_variables= photon_isolation_variables(input_data, extra_photons)
+    composite_tuple_variables = {
+        "Dst0"   :  make_composite_variables(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs,  "composite") +\
+             MVA_dst_variables("Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip", pvs) + make_DeltaM_variable(options), 
+        "D0"     :  make_composite_variables(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs,  "composite"),
+        "gamma"  :  make_composite_variables(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs,  "composite")+extra_variables
+    }
+    basic_tuple_variables = { 
+        "Kminus" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic"),
+        "piplus" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic"),
+        "eplus"  : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic") + extra_brem_variables,
+        "eminus" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic") + extra_brem_variables,
+    }
+
+    my_tuple = Funtuple(
+        name=f"Tuple_{line}_{dd}",
+        tuple_name="DecayTree",
+        fields=fields,
+        variables= {**composite_tuple_variables, **basic_tuple_variables},
+        event_variables=make_hlt2_event_variables(options, pvs, rec_summary), 
+        inputs=input_data,
+        store_multiple_cand_info=True
+    )
+    return [ my_filter, my_tuple]
+
+def make_Dst0ToD0EmEp_D0ToKmPimPipPip_tuple(options, line, pvs, rec_summary, dd='dst_k3pi_os'):
+
+    init_data=get_particles(f"/Event/HLT2/{line}/Particles")
+    extra_photons = get_particles(f"/Event/HLT2/{line}/{line[10:]}_prompt_Photons/Particles")
+    my_filter = create_lines_filter(f"LineFilter_{line}_{{hash}}",[line],)
+    input_data = ParticleFilter(init_data, F.FILTER(F.require_all(
+                    in_range(1900. * MeV, F.MASS, 2100. * MeV),
+                    in_range(1825. * MeV, CHILD_1(F.MASS), 1900. * MeV),
+                    F.MAXDOCA < 0.2,
+                    CHILD_1(F.MAXDOCA) < 0.2,
+                    CHILD_2(F.MAXDOCA) < 0.2,)))
+
+    fields = decay_descriptor[dd]
+    
+    extra_variables= photon_isolation_variables(input_data, extra_photons)
+    composite_tuple_variables = {
+        "Dst0"   : make_composite_variables(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs, "composite") +\
+             MVA_dst_variables("Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip", pvs) + make_DeltaM_variable(options), 
+        "D0"     : make_composite_variables_4body(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs, "composite", nparticles=4),
+        "gamma"  : make_composite_variables(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs,  "composite")+extra_variables 
+    }
+    basic_tuple_variables = { 
+        "Kminus"  : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic"),
+        "piminus" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic"),
+        "piplus1" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic"),
+        "piplus2" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic"), 
+        "eplus"   : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic") + extra_brem_variables,
+        "eminus"  : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic") + extra_brem_variables,
+    }
+
+    my_tuple = Funtuple(
+        name=f"Tuple_{line}_{dd}",
+        tuple_name="DecayTree",
+        fields=fields,
+        variables= {**composite_tuple_variables, **basic_tuple_variables},
+        event_variables=make_hlt2_event_variables(options, pvs, rec_summary),
+        inputs=input_data,
+        store_multiple_cand_info=True
+    )
+    return [ my_filter, my_tuple]
+
+
+def make_DstpToDpDspEmEp_DpDspToKmKpPip_tuple(options, line, pvs , rec_summary, dd='dsstp_2kpi_os'):
+
+    init_data=get_particles(f"/Event/HLT2/{line}/Particles")
+    extra_photons = get_particles(f"/Event/HLT2/{line}/{line[10:]}_prompt_Photons/Particles")
+    my_filter = create_lines_filter(f"LineFilter_{line}_{{hash}}",[line],)
+    input_data = ParticleFilter(init_data, F.FILTER(F.require_all(
+                    in_range(1900. * MeV, F.MASS, 2300. * MeV), 
+                    in_range(1900. * MeV, CHILD_1(F.MASS), 2049. * MeV),
+                    F.MAXDOCA < 0.2,
+                    CHILD_1(F.MAXDOCA) < 0.2,
+                    CHILD_2(F.MAXDOCA) < 0.2,)))
+    fields = decay_descriptor[dd]
+
+    extra_variables= photon_isolation_variables(input_data, extra_photons)
+    composite_tuple_variables = {
+        "Dstp"   : make_composite_variables(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs, "composite", mass_constrain = ["D*_s+", "D_s+"]) +\
+             MVA_dst_variables("Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip", pvs) + make_DeltaM_variable(options), 
+        "DpDs"   : make_composite_variables_3body(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs, "composite", mass_constrain = ["D*_s+", "D_s+"], nparticles=3),
+        "gamma"  : make_composite_variables(options, pvs, input_data) + make_dtf_variables(options, input_data, pvs, "composite", mass_constrain = ["D*_s+", "D_s+"])+extra_variables,
+    }
+    basic_tuple_variables = { 
+        "Kminus" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic",  mass_constrain = ["D*_s+", "D_s+"]),
+        "Kplus"  : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic",  mass_constrain = ["D*_s+", "D_s+"]),
+        "piplus" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic",  mass_constrain = ["D*_s+", "D_s+"]),
+        "eplus"  : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic",  mass_constrain = ["D*_s+", "D_s+"])+ extra_brem_variables,
+        "eminus" : make_basic_variables(options, pvs, input_data)+make_dtf_variables(options, input_data, pvs, "basic",  mass_constrain = ["D*_s+", "D_s+"])+ extra_brem_variables,
+    }
+
+    my_tuple = Funtuple(
+        name=f"Tuple_{line}_{dd}",
+        tuple_name="DecayTree",
+        fields=fields,
+        variables= {**composite_tuple_variables, **basic_tuple_variables},
+        event_variables=make_hlt2_event_variables(options, pvs, rec_summary),
+        inputs=input_data,
+        store_multiple_cand_info=True
+    )
+    return [ my_filter, my_tuple]
+
+from DaVinci import Options, make_config
+from PyConf.reading import get_pvs, get_rec_summary
+
+def main(options: Options):
+
+    # get ODIN and DecReports location
+    pvs = get_pvs()
+    rec_summary = get_rec_summary()
+
+    tuples = {
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_OS':             make_Dst0ToD0EmEp_D0ToKmPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_OS", pvs, rec_summary, 'dst_kpi_os'),
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_SS':             make_Dst0ToD0EmEp_D0ToKmPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_SS", pvs, rec_summary, 'dst_kpi_ss'),
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_OS_MVA':         make_Dst0ToD0EmEp_D0ToKmPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_OS_MVA", pvs, rec_summary, 'dst_kpi_os'),
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_SS_MVA':         make_Dst0ToD0EmEp_D0ToKmPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_SS_MVA", pvs, rec_summary, 'dst_kpi_ss'),
+
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_OS':       make_Dst0ToD0EmEp_D0ToKmPimPipPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_OS", pvs, rec_summary,'dst_k3pi_os' ),
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_SS':       make_Dst0ToD0EmEp_D0ToKmPimPipPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_SS", pvs, rec_summary, 'dst_k3pi_ss'),
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_OS_MVA':   make_Dst0ToD0EmEp_D0ToKmPimPipPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_OS_MVA", pvs, rec_summary, 'dst_k3pi_os'),
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_SS_MVA':   make_Dst0ToD0EmEp_D0ToKmPimPipPip_tuple(options, "Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_SS_MVA", pvs, rec_summary, 'dst_k3pi_ss'),
+
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_OS':     make_DstpToDpDspEmEp_DpDspToKmKpPip_tuple(options, "Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_OS", pvs, rec_summary, "dsstp_2kpi_os"),
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_SS':     make_DstpToDpDspEmEp_DpDspToKmKpPip_tuple(options, "Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_SS", pvs, rec_summary, "dsstp_2kpi_ss"),
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_OS_MVA': make_DstpToDpDspEmEp_DpDspToKmKpPip_tuple(options, "Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_OS_MVA", pvs, rec_summary, "dsstp_2kpi_os"),
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_SS_MVA': make_DstpToDpDspEmEp_DpDspToKmKpPip_tuple(options, "Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_SS_MVA", pvs, rec_summary, "dsstp_2kpi_ss"),
+    }
+
+    config = make_config(options, tuples)
+
+    return config
diff --git a/dst_to_dee/README.md b/dst_to_dee/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..ef730cb6d64738b78c8816fb19e70df9ba6a0f5a
--- /dev/null
+++ b/dst_to_dee/README.md
@@ -0,0 +1,5 @@
+# Analysis production for dst_to_dee charm module.
+Creates tuples for Charm 2024-validation MC and data
+
+## Analysis motivation
+The module dst_to_dee contains several HLT2 lines to search for dark photons produced by a Charm decay.
\ No newline at end of file
diff --git a/dst_to_dee/info.yaml b/dst_to_dee/info.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..ef6c74cf31720afbc501a8ae964d1fb0a349f438
--- /dev/null
+++ b/dst_to_dee/info.yaml
@@ -0,0 +1,42 @@
+defaults:
+  inform:
+    - carlos.eduardo.cocha.toapaxi@cern.ch
+  wg: Charm
+
+{%- set datainfo = [
+( 'MagUp',   '24c2'),
+( 'MagDown', '24c2'),
+( 'MagUp',   '24c3'),
+( 'MagDown', '24c3'),
+( 'MagUp',   '24c4'),
+( 'MagDown', '24c4'),
+]%}
+
+{%- for  polarity,  sprucen in datainfo %}
+DATA_2024_{{ polarity }}_{{sprucen}}:
+  application: "DaVinci/v64r10@x86_64_v2-el9-clang16-opt"
+  input:       
+    bk_query: "/LHCb/Collision24/Beam6800GeV-VeloClosed-{{polarity}}/Real Data/Sprucing{{sprucen}}/94000000/CHARM.DST"
+    dq_flags:
+      - UNCHECKED
+      - OK
+    keep_running: True
+  output: DVTUPLE.ROOT
+  options:
+    entrypoint: dst_to_dee.DV:main
+    extra_options:
+      input_raw_format: 0.5
+      input_type: ROOT
+      simulation: False
+      data_type: "Upgrade"
+      event_store: HiveWhiteBoard
+      geometry_version: 'run3/2024.Q1.2-v00.00'
+      conditions_version: master
+      input_process: "TurboPass"
+      input_stream: "charm"
+      compression:
+        algorithm: ZSTD
+        level: 1
+        max_buffer_size: 1048576
+
+{%- endfor %}
diff --git a/dst_to_dee/tupling.py b/dst_to_dee/tupling.py
new file mode 100644
index 0000000000000000000000000000000000000000..8a7a9d2730286ce2de9cf347d5ff84c03a198fb5
--- /dev/null
+++ b/dst_to_dee/tupling.py
@@ -0,0 +1,1023 @@
+##############################################################################
+# (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.                                       #
+###############################################################################
+"""Common configuration functions
+
+"""
+
+import Functors as F
+from Functors.math import log
+from Functors.math import in_range
+from GaudiKernel.SystemOfUnits import  MeV
+from FunTuple import FunctorCollection
+from FunTuple.functorcollections import (
+    MCHierarchy,
+    MCPromptDecay,
+    Kinematics,
+    SelectionInfo,
+    HltTisTos,
+    MCVertexInfo,
+    MCKinematics,
+    ParticleID, #wrong variables PID_PI = 0, PROBNN_D = nan
+    EventInfo,
+    LHCInfo,
+    ParticleIsolation,
+    MCPrimaries,
+    MCReconstructed,
+    MCReconstructible,
+)
+
+from DaVinciMCTools import MCTruthAndBkgCat, MCReconstructed, MCReconstructible
+from PyConf.Algorithms import ParticleToSubcombinationsAlg
+from DecayTreeFitter import DecayTreeFitter
+from PyConf.Algorithms import WeightedRelTableAlg, ThOrParticleSelection
+FILTER_TREE = lambda id: F.FILTER(F.IS_ABS_ID(id)) @ F.GET_ALL_DESCENDANTS()
+Hlt1_global_lines = [
+    "Hlt1GECPassthroughDecision",
+    "Hlt1BeamGasDecision",
+    "Hlt1PassthroughDecision",
+    "Hlt1NoBeamDecision",
+    "Hlt1BeamOneDecision",
+    "Hlt1BeamTwoDecision",
+    "Hlt1BothBeamsDecision",
+    "Hlt1ODINLumiDecision",
+    "Hlt1ODINVeloOpenDecision",
+    "Hlt1ODINNoBiasDecision",
+    "Hlt1VeloMicroBiasDecision",
+    "Hlt1RICH1AlignmentDecision",
+    "Hlt1RICH2AlignmentDecision",
+    "Hlt1BeamGasDecision",
+    "Hlt1L02PPiDecision",
+    "Hlt1LowMassNoipDielectron_massSlice1_promptDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice1_promptDecision",
+    "Hlt1LowMassNoipDielectron_massSlice2_promptDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice2_promptDecision",
+    "Hlt1LowMassNoipDielectron_massSlice3_promptDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice3_promptDecision",
+    "Hlt1LowMassNoipDielectron_massSlice4_promptDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice4_promptDecision",
+    "Hlt1LowMassNoipDielectron_massSlice1_displacedDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice1_displacedDecision",
+    "Hlt1LowMassNoipDielectron_massSlice2_displacedDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice2_displacedDecision",
+    "Hlt1LowMassNoipDielectron_massSlice3_displacedDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice3_displacedDecision",
+    "Hlt1LowMassNoipDielectron_massSlice4_displacedDecision",
+    "Hlt1LowMassNoipDielectron_SS_massSlice4_displacedDecision",
+]
+Hlt1_1track_lines = [
+    "Hlt1TrackMVADecision",
+    "Hlt1LowPtMuonDecision",
+    "Hlt1SingleHighPtMuonDecision",
+    "Hlt1SingleHighPtMuonNoMuIDDecision",
+    "Hlt1TrackMuonMVADecision",
+    "Hlt1OneMuonTrackLineDecision",
+    "Hlt1TrackElectronMVADecision",
+    "Hlt1SingleHighPtElectronDecision",
+    "Hlt1SingleHighEtDecision",
+]
+Hlt1_lines = Hlt1_1track_lines+[
+    "Hlt1TwoTrackMVACharmXSecDecision",
+    "Hlt1TwoTrackMVADecision",
+    "Hlt1TwoTrackKsDecision",
+    "Hlt1D2KPiDecision",
+    "Hlt1D2KKDecision",
+    "Hlt1D2PiPiDecision",
+    "Hlt1KsToPiPiDecision",
+    "Hlt1LowPtDiMuonDecision",#removed
+    "Hlt1DiMuonNoIPDecision",
+    "Hlt1DiMuonNoIP_ssDecision",
+    "Hlt1DiMuonHighMassDecision",
+    "Hlt1DiMuonLowMassDecision",#replaced by Hlt1DiMuonDisplacedDecision
+    "Hlt1DiMuonSoftDecision",
+    "Hlt1DiMuonDisplacedDecision",
+    "Hlt1TwoKsDecision",
+    "Hlt1D2KPiAlignmentDecision",
+    "Hlt1DiMuonHighMassAlignmentDecision",
+    "Hlt1DisplacedDiMuonAlignmentDecision",
+    "Hlt1DisplacedDielectronDecision",
+    "Hlt1DisplacedLeptonsDecision",#removed
+]
+
+
+Hlt2_lines = [
+    "Hlt2Charm_DstpToD0Pip_D0ToKmPip_XSec",
+    "Hlt2Charm_D0ToKmPip_XSec",
+    "Hlt2Charm_D0ToKmKp",
+    "Hlt2Charm_D0ToKmPip",
+    "Hlt2Charm_D0ToPimPip",
+    "Hlt2Charm_DpDspToKsKp_DD",
+    "Hlt2Charm_DpDspToKsKp_LD",
+    "Hlt2Charm_DpDspToKsKp_LL",
+    "Hlt2Charm_DpDspToKsPip_DD",
+    "Hlt2Charm_DpDspToKsPip_LD",
+    "Hlt2Charm_DpDspToKsPip_LL",
+    "Hlt2Charm_DpToKmPipPip",
+    "Hlt2Charm_DspToKmKpPip",
+    "Hlt2Charm_DpToKmPipPip_NoCuts",
+    "Hlt2Charm_DspToKmKpPip_NoCuts",
+    "Hlt2Charm_DpToKmPipPip_XSec",
+    "Hlt2Charm_DspToKmKpPip_XSec",
+    "Hlt2Charm_DstpToD0Pip_D0ToKmKp",
+    "Hlt2Charm_DstpToD0Pip_D0ToKmPip",
+    "Hlt2Charm_DstpToD0Pip_D0ToKpPim",
+    "Hlt2Charm_DstpToD0Pip_D0ToPimPip",
+    "Hlt2Charm_DstpToD0Pip_D0ToKmPip_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKpPim_LowBias",
+    "Hlt2Charm_D0ToKmPip_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKmKp_LowBias",
+    "Hlt2Charm_D0ToKmKp_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToPimPip_LowBias",
+    "Hlt2Charm_D0ToPimPip_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_LLLL",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_LLLL_Tight",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_LLDD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_LLDD_Tight",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_DDDD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_DDDD_Tight",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_ULLL",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_ULLL_Tight",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_ULDD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_ULDD_Tight",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_LLLD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_LLLD_Tight",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_DDLD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKs_DDLD_Tight",
+    "Hlt2Charm_D0ToKsPimPip_LL",
+    "Hlt2Charm_D0ToKsPimPip_DD",
+    "Hlt2Charm_D0ToKsPimPip_LL_LowBias",
+    "Hlt2Charm_D0ToKsPimPip_LL_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsPimPip_LL",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsPimPip_DD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsPimPip_LL_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsPimPip_DD_LowBias",
+    "Hlt2Charm_D0ToKsKmPip_LL",
+    "Hlt2Charm_D0ToKsKmPip_DD",
+    "Hlt2Charm_D0ToKsKmPip_LL_LowBias",
+    "Hlt2Charm_D0ToKsKmPip_LL_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmPip_LL",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmPip_DD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmPip_LL_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmPip_DD_LowBias",
+    "Hlt2Charm_D0ToKsKpPim_LL",
+    "Hlt2Charm_D0ToKsKpPim_DD",
+    "Hlt2Charm_D0ToKsKpPim_LL_LowBias",
+    "Hlt2Charm_D0ToKsKpPim_DD_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKpPim_LL",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKpPim_DD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKpPim_LL_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKpPim_DD_LowBias",
+    "Hlt2Charm_D0ToKsKmKp_LL",
+    "Hlt2Charm_D0ToKsKmKp_DD",
+    "Hlt2Charm_D0ToKsKmKp_LL_LowBias",
+    "Hlt2Charm_D0ToKsKmKp_DD_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmKp_LL",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmKp_DD",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmKp_LL_LowBias",
+    "Hlt2Charm_DstpToD0Pip_D0ToKsKmKp_DD_LowBias",
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_OS',
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_SS',
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_OS_MVA',
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPip_SS_MVA',
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_OS',
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_SS',
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_OS_MVA',
+    'Hlt2Charm_Dst0ToD0EmEp_D0ToKmPimPipPip_SS_MVA',
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_OS',
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_SS',
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_OS_MVA',
+    'Hlt2Charm_DstpToDpDspEmEp_DpDspToKmKpPip_SS_MVA',
+]
+
+def make_composite_variables(options, pvs, data, add_truth=True, add_Hlt1TisTos=True):
+    if not options.simulation:
+        add_truth = False
+    variables = (
+        FunctorCollection(
+            {
+                "MAXPT": F.MAX(F.PT),
+                "MINPT": F.MIN(F.PT),
+                "SUMPT": F.SUM(F.PT),
+                "MAXP": F.MAX(F.P),
+                "MINP": F.MIN(F.P),
+                "BPVDIRA": F.BPVDIRA(pvs),
+                "BPVDLS": F.BPVDLS(pvs),
+                "VCHI2DOF": F.CHI2DOF,
+                #"VNDOF": F.NDOF,
+                "BPVFDCHI2": F.BPVFDCHI2(pvs),
+                "BPVFD": F.BPVFD(pvs),
+                "BPVVDRHO": F.BPVVDRHO(pvs),
+                "BPVVDZ": F.BPVVDZ(pvs),
+                "BPVIPCHI2": F.BPVIPCHI2(pvs),
+                "BPVIP": F.BPVIP(pvs),
+                "LOGBPVIPCHI2": log(F.BPVIPCHI2(pvs)),
+                "BPVLTIME": F.BPVLTIME(pvs),
+                "MAXBPVIPCHI2": F.MAX(F.BPVIPCHI2(pvs)),
+                "MINBPVIPCHI2": F.MIN(F.BPVIPCHI2(pvs)),
+                "MAXBPVIP": F.MAX(F.BPVIP(pvs)),
+                "MINBPVIP": F.MIN(F.BPVIP(pvs)),
+                "MAXDOCACHI2": F.MAXDOCACHI2,
+                "MAXDOCA": F.MAXDOCA,
+                "MAXSDOCACHI2": F.MAXSDOCACHI2,
+                "MAXSDOCA": F.MAXSDOCA,
+                "COS12": F.ALV(1, 2),
+                "ETA": F.ETA,
+                "PHI": F.PHI,
+                "END_VX": F.END_VX,
+                "END_VY": F.END_VY,
+                "END_VZ": F.END_VZ,
+                "END_VX_ERR":F.SQRT @ F.CALL(0,0) @ F.POS_COV_MATRIX @ F.ENDVERTEX,
+                "END_VY_ERR":F.SQRT @ F.CALL(1,1) @ F.POS_COV_MATRIX @ F.ENDVERTEX,
+                "END_VZ_ERR":F.SQRT @ F.CALL(2,2) @ F.POS_COV_MATRIX @ F.ENDVERTEX,
+                "BPVX": F.BPVX(pvs),
+                "BPVY": F.BPVY(pvs),
+                "BPVZ": F.BPVZ(pvs),
+                "BPVX_ERR": F.SQRT @ F.CALL(0,0) @ F.POS_COV_MATRIX @ F.BPV(pvs),
+                "BPVY_ERR": F.SQRT @ F.CALL(1,1) @ F.POS_COV_MATRIX @ F.BPV(pvs),
+                "BPVZ_ERR": F.SQRT @ F.CALL(2,2) @ F.POS_COV_MATRIX @ F.BPV(pvs),
+                "ALLPVFD"     : F.ALLPV_FD(pvs),
+                "ALLPVIP"     : F.ALLPV_IP(pvs),
+                "OBJECT_KEY": F.OBJECT_KEY,
+            }
+        )
+        + Kinematics()
+        #+ ParticleID(extra_info=True) #only for daughters
+    )
+
+    if add_Hlt1TisTos:
+        variables += HltTisTos(
+            selection_type="Hlt1", trigger_lines=Hlt1_lines, data=data
+        )
+
+    if add_truth:
+        variables = add_truth_matching_functors(
+            options,
+            variables,
+            data,
+            hierarchy=True,
+            kinematics=True,
+            vertex_info=True,
+            bkgcat=True,
+        )
+
+    return variables
+
+def make_tistoscombinations(options, pvs, data):
+    algo_output = ParticleToSubcombinationsAlg( Input=data )
+    relations = algo_output.OutputRelations
+
+    tistos_variables_extra = HltTisTos(
+        selection_type="Hlt1", trigger_lines=["Hlt1TwoTrackMVADecision"], data=algo_output.OutputParticles
+    )
+    tistos_combinations = {}
+    for k,v in tistos_variables_extra.get_thor_functors().items():
+        tistos_combinations[k + "_SUBCOMB" ] = F.MAP_INPUT_ARRAY( v, relations )
+    tistos_combinations = FunctorCollection( tistos_combinations )
+
+    return tistos_combinations
+
+def make_composite_variables_3body(options, pvs, data, add_truth=True):
+
+    variables = (
+        FunctorCollection(
+            {
+                #12
+                "M12": F.SUBCOMB(Functor=F.MASS, Indices=(1, 2)),
+                "SDOCA12"     : F.SDOCA(Child1=1,Child2=2),
+                "SDOCACHI212" : F.SDOCACHI2(Child1=1,Child2=2),
+                "DOCA12"      : F.DOCA(Child1=1,Child2=2),
+                "DOCACHI212"  : F.DOCACHI2(Child1=1,Child2=2),
+                "COS12": F.ALV(1, 2),
+                #"ETA12": F.SUBCOMB(Functor=F.ETA, Indices=(1, 2)),
+                #"PT12": F.SUBCOMB(Functor=F.PT, Indices=(1, 2)),
+                #"VCHI212": F.SUBCOMB(Functor=F.CHI2DOF, Indices=(1, 2)),
+                #"END_VZ12": F.SUBCOMB(Functor=F.END_VZ, Indices=(1, 2)),
+                #"BPVZ12": F.SUBCOMB(Functor=F.BPVZ(pvs), Indices=(1, 2)),
+                #"BPVCORRM12": F.SUBCOMB(Functor=F.BPVCORRM(pvs), Indices=(1, 2)),
+                #13
+                "M13": F.SUBCOMB(Functor=F.MASS, Indices=(1, 3)),
+                "SDOCA13"     : F.SDOCA(Child1=1,Child2=3),
+                "SDOCACHI213" : F.SDOCACHI2(Child1=1,Child2=3),
+                "DOCA13"      : F.DOCA(Child1=1,Child2=3),
+                "DOCACHI213"  : F.DOCACHI2(Child1=1,Child2=3),
+                "COS13": F.ALV(1, 3),
+                #"ETA13": F.SUBCOMB(Functor=F.ETA, Indices=(1, 3)),
+                #"PT13": F.SUBCOMB(Functor=F.PT, Indices=(1, 3)),
+                #"VCHI213": F.SUBCOMB(Functor=F.CHI2DOF, Indices=(1, 3)),
+                #"END_VZ13": F.SUBCOMB(Functor=F.END_VZ, Indices=(1, 3)),
+                #"BPVZ13": F.SUBCOMB(Functor=F.BPVZ(pvs), Indices=(1, 3)),
+                #"BPVCORRM13": F.SUBCOMB(Functor=F.BPVCORRM(pvs), Indices=(1, 3)),
+                #23
+                "M23": F.SUBCOMB(Functor=F.MASS, Indices=(2, 3)),
+                "SDOCA23"     : F.SDOCA(Child1=2,Child2=3),
+                "SDOCACHI223" : F.SDOCACHI2(Child1=2,Child2=3),
+                "DOCA23"      : F.DOCA(Child1=2,Child2=3),
+                "DOCACHI223"  : F.DOCACHI2(Child1=2,Child2=3),
+                "COS23": F.ALV(2, 3),
+                #"ETA23": F.SUBCOMB(Functor=F.ETA, Indices=(2, 3)),
+                #"PT23": F.SUBCOMB(Functor=F.PT, Indices=(2, 3)),
+                #"VCHI223": F.SUBCOMB(Functor=F.CHI2DOF, Indices=(2, 3)),
+                #"END_VZ23": F.SUBCOMB(Functor=F.END_VZ, Indices=(2, 3)),
+                #"BPVZ23": F.SUBCOMB(Functor=F.BPVZ(pvs), Indices=(2, 3)),
+                #"BPVCORRM23": F.SUBCOMB(Functor=F.BPVCORRM(pvs), Indices=(2, 3)),
+            }
+        )
+    )
+    return make_composite_variables(options, pvs, data, add_truth)+make_tistoscombinations(options, pvs, data)+variables
+
+def make_composite_variables_4body(options, pvs, data, add_truth=True):
+
+    variables = (
+        FunctorCollection(
+            {
+                #14
+                "M14": F.SUBCOMB(Functor=F.MASS, Indices=(1, 4)),
+                "SDOCA14"     : F.SDOCA(Child1=1,Child2=4),
+                "SDOCACHI214" : F.SDOCACHI2(Child1=1,Child2=4),
+                "DOCA14"      : F.DOCA(Child1=1,Child2=4),
+                "DOCACHI214"  : F.DOCACHI2(Child1=1,Child2=4),
+                "COS14": F.ALV(1, 4),
+                #"ETA14": F.SUBCOMB(Functor=F.ETA, Indices=(1, 4)),
+                #"PT14": F.SUBCOMB(Functor=F.PT, Indices=(1, 4)),
+                #"VCHI214": F.SUBCOMB(Functor=F.CHI2DOF, Indices=(1, 4)),
+                #"END_VZ14": F.SUBCOMB(Functor=F.END_VZ, Indices=(1, 4)),
+                #"BPVZ14": F.SUBCOMB(Functor=F.BPVZ(pvs), Indices=(1, 4)),
+                #"BPVCORRM14": F.SUBCOMB(Functor=F.BPVCORRM(pvs), Indices=(1, 4)),
+                #24
+                "M24": F.SUBCOMB(Functor=F.MASS, Indices=(2, 4)),
+                "SDOCA24"     : F.SDOCA(Child1=2,Child2=4),
+                "SDOCACHI224" : F.SDOCACHI2(Child1=2,Child2=4),
+                "DOCA24"      : F.DOCA(Child1=2,Child2=4),
+                "DOCACHI224"  : F.DOCACHI2(Child1=2,Child2=4),
+                "COS24": F.ALV(2, 4),
+                #"ETA24": F.SUBCOMB(Functor=F.ETA, Indices=(2, 4)),
+                #"PT24": F.SUBCOMB(Functor=F.PT, Indices=(2, 4)),
+                #"VCHI224": F.SUBCOMB(Functor=F.CHI2DOF, Indices=(2, 4)),
+                #"END_VZ24": F.SUBCOMB(Functor=F.END_VZ, Indices=(2, 4)),
+                #"BPVZ24": F.SUBCOMB(Functor=F.BPVZ(pvs), Indices=(2, 4)),
+                #"BPVCORRM24": F.SUBCOMB(Functor=F.BPVCORRM(pvs), Indices=(2, 4)),
+                #34
+                "M34": F.SUBCOMB(Functor=F.MASS, Indices=(3, 4)),
+                "SDOCA34"     : F.SDOCA(Child1=3,Child2=4),
+                "SDOCACHI234" : F.SDOCACHI2(Child1=3,Child2=4),
+                "DOCA34"      : F.DOCA(Child1=3,Child2=4),
+                "DOCACHI234"  : F.DOCACHI2(Child1=3,Child2=4),
+                "COS34": F.ALV(3, 4),
+                #"ETA34": F.SUBCOMB(Functor=F.ETA, Indices=(3, 4)),
+                #"PT34": F.SUBCOMB(Functor=F.PT, Indices=(3, 4)),
+                #"VCHI234": F.SUBCOMB(Functor=F.CHI2DOF, Indices=(3, 4)),
+                #"END_VZ34": F.SUBCOMB(Functor=F.END_VZ, Indices=(3, 4)),
+                #"BPVZ34": F.SUBCOMB(Functor=F.BPVZ(pvs), Indices=(3, 4)),
+                #"BPVCORRM34": F.SUBCOMB(Functor=F.BPVCORRM(pvs), Indices=(3, 4)),
+            }
+        )
+    )
+    return make_composite_variables_3body(options, pvs, data, add_truth)+variables #make_tistoscombinations already included
+
+def make_DeltaM_variable(options):
+    return FunctorCollection({"DM": F.MASS - F.CHILD(1, F.MASS)})
+
+def make_basic_variables(options, pvs, data, add_truth=True):
+    if not options.simulation:
+        add_truth = False
+    variables = (
+        FunctorCollection(
+            {
+                "TRCHI2DOF": F.CHI2DOF @ F.TRACK,
+                "ETA": F.ETA,
+                "PHI": F.PHI,
+                "TRGHOSTPROB": F.GHOSTPROB,
+                "BPVIPCHI2": F.BPVIPCHI2(pvs),
+                "BPVIP": F.BPVIP(pvs),
+                "BPVX": F.BPVX(pvs),
+                "BPVY": F.BPVY(pvs),
+                "BPVZ": F.BPVZ(pvs),
+                "TX"          : F.TX,
+                "TY"          : F.TY,
+                "MINIPCHI2"   : F.MINIPCHI2(pvs),
+                "MINIP"       : F.MINIP(pvs),
+                "KEY"         : F.VALUE_OR(-1) @ F.OBJECT_KEY @ F.TRACK,
+                "CTB"         : F.POSITION @ F.CLOSESTTOBEAM @ F.TRACK,
+                "ISMUON"      : F.ISMUON,
+                "TRACK_P": F.VALUE_OR(-1) @ F.P @ F.TRACK,
+                "TRACKPT": F.VALUE_OR(-1) @ F.TRACK_PT,
+                #"TRACK_PT_WITH_BREM": F.VALUE_OR(-1) @ F.PT_WITH_BREM @ F.TRACK,
+                #"TRACK_P_WITH_BREM": F.VALUE_OR(-1) @ F.P_WITH_BREM @ F.TRACK,
+                "TRACKHISTORY": F.VALUE_OR(-1) @ F.TRACKHISTORY @ F.TRACK,
+                "QOVERP": F.QOVERP @ F.TRACK,
+                "NDOF": F.VALUE_OR(-1) @ F.NDOF @ F.TRACK,
+                "NFTHITS": F.VALUE_OR(-1) @ F.NFTHITS @ F.TRACK,
+                "NHITS": F.VALUE_OR(-1) @ F.NHITS @ F.TRACK,
+                "NUTHITS": F.VALUE_OR(-1) @ F.NUTHITS @ F.TRACK,
+                "NVPHITS": F.VALUE_OR(-1) @ F.NVPHITS @ F.TRACK,
+                #"NMUONHITS": F.VALUE_OR(-1) @ F.NHITSMUONS @ F.TRACK, #https://gitlab.cern.ch/lhcb/Rec/-/merge_requests/3756
+                "TRACKHASVELO": F.VALUE_OR(-1) @ F.TRACKHASVELO @ F.TRACK,
+                "TRACKHASUT": F.VALUE_OR(-1) @ F.TRACKHASUT @ F.TRACK,
+                "TRACKHISTORY": F.VALUE_OR(-1) @ F.TRACKHISTORY @ F.TRACK,
+                "TRACKISCLONE": F.VALUE_OR(-1) @ F.TRACKISCLONE @ F.TRACK,
+                "TRACKISDOWNSTREAM": F.VALUE_OR(-1) @ F.TRACKISDOWNSTREAM @ F.TRACK,
+                "TRACKISINVALID": F.VALUE_OR(-1) @ F.TRACKISINVALID @ F.TRACK,
+                "TRACKISLONG": F.VALUE_OR(-1) @ F.TRACKISLONG @ F.TRACK,
+                "TRACKISSELECTED": F.VALUE_OR(-1) @ F.TRACKISSELECTED @ F.TRACK,
+                "TRACKISTTRACK": F.VALUE_OR(-1) @ F.TRACKISTTRACK @ F.TRACK,
+                "TRACKISUPSTREAM": F.VALUE_OR(-1) @ F.TRACKISUPSTREAM @ F.TRACK,
+                "TRACKISVELO": F.VALUE_OR(-1) @ F.TRACKISVELO @ F.TRACK,
+                "TRACKISVELOBACKWARD": F.VALUE_OR(-1) @ F.TRACKISVELOBACKWARD @ F.TRACK,
+                "OBJECT_KEY": F.OBJECT_KEY,
+                "HASBREM": F.HASBREM,
+                "BREMENERGY": F.BREMENERGY,
+                "BREMBENDCORR": F.BREMBENDCORR,
+                "PID_P":F.VALUE_OR(F.NaN)@F.PID_P, 
+                "PID_PI": F.VALUE_OR(F.NaN)@F.PID_PI,
+                "PID_K":F.VALUE_OR(F.NaN)@F.PID_K, 
+                "PID_E":F.VALUE_OR(F.NaN)@F.PID_E, 
+                "RICH_DLL_E": F.VALUE_OR(F.NaN)@F.RICH_DLL_E, 
+                "PROBNN_E":F.VALUE_OR(F.NaN)@F.PROBNN_E,
+                "GHOSTPROB":F.GHOSTPROB, 
+                "CHARGE":F.CHARGE,
+                "NVPHITSA": F.VALUE_OR(-1) @ F.NVPHITSA @ F.TRACK,
+                "NVPHITSC": F.VALUE_OR(-1) @ F.NVPHITSC @ F.TRACK,
+                "NVPOVERLAP": F.VALUE_OR(-1) @ F.NVPOVERLAP @ F.TRACK,
+                "POSITION_STATEAT_FirstMeasurement" : F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_LastMeasurement" : F.STATE_AT("LastMeasurement") @ F.TRACK ,
+                "STATEAT_EndVelo"                : F.EXTRAPOLATE_TRACK(770.0) @ F.TRACK,
+                "POSITION_STATEAT_EndVelo_X"     : F.POSITION_X   @ F.EXTRAPOLATE_TRACK(770.0) @ F.TRACK,
+                "POSITION_STATEAT_EndVelo_Y"     : F.POSITION_Y   @ F.EXTRAPOLATE_TRACK(770.0) @ F.TRACK,
+                "POSITION_STATEAT_EndVelo_Z"     : F.POSITION_Z   @ F.EXTRAPOLATE_TRACK(770.0) @ F.TRACK,
+                "POSITION_STATEAT_EndVelo_TX"    : F.X_COORDINATE @ F.SLOPES @ F.EXTRAPOLATE_TRACK(770.0) @ F.TRACK,
+                "POSITION_STATEAT_EndVelo_TY"    : F.Y_COORDINATE @ F.SLOPES @ F.EXTRAPOLATE_TRACK(770.0) @ F.TRACK,
+                "POSITION_STATEAT_FirstMeasurement_X"     : F.POSITION_X   @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_Y"     : F.POSITION_Y   @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_Z"     : F.POSITION_Z   @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_TX"    : F.X_COORDINATE @ F.SLOPES @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_TY"    : F.Y_COORDINATE @ F.SLOPES @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+            }
+        )
+        + Kinematics()
+        + ParticleID(extra_info=True)
+    )
+
+    variables += HltTisTos(
+        selection_type="Hlt1", trigger_lines=Hlt1_1track_lines, data=data
+    )
+
+    if add_truth:
+        variables = add_truth_matching_functors(
+            options,
+            variables,
+            data,
+            hierarchy=True,
+            kinematics=True,
+            vertex_info=True,
+            bkgcat=False,
+        )
+
+    return variables
+
+def make_hlt2_event_variables(options, pvs, rec_summary):
+    # define event level variables
+    evt_variables = EventInfo()
+    if not options.simulation:
+        evt_variables += LHCInfo() #FillNumber
+    evt_variables += SelectionInfo(selection_type="Hlt2", trigger_lines=Hlt2_lines)
+    evt_variables += SelectionInfo(selection_type="Hlt1", trigger_lines=Hlt1_lines)
+    evt_variables += FunctorCollection(
+        {
+            "ALLPVX": F.ALLPVX(pvs),
+            "ALLPVY": F.ALLPVY(pvs),
+            "ALLPVZ": F.ALLPVZ(pvs),
+            "ALLPVNDOF": F.MAP(F.VALUE_OR(F.NaN) @ F.NDOF) @ F.TES(pvs),
+            "ALLPVCHI2": F.MAP(F.VALUE_OR(F.NaN) @ F.CHI2) @ F.TES(pvs),
+            "nPVs": F.SIZE(pvs),
+            "nTracks": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nTracks"),
+            "nLongTracks": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nLongTracks"),
+            "nMuonTracks": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nMuonTracks"),
+            "nFTClusters": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nFTClusters"),
+            "nVPClusters": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nVPClusters"),
+            "nUTClusters": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nUTClusters"),
+            "nSPDhits": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nSPDhits"),
+            "nEcalClusters": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nEcalClusters"),
+            "nEcalClusters": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nEcalClusters"),
+            "eCalTot": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "eCalTot"),
+            "hCalTot": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "hCalTot"),
+            "nRich1Hits": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nRich1Hits"),
+            "nRich2Hits": F.VALUE_OR(-1) @ F.RECSUMMARY_INFO(rec_summary, "nRich2Hits"),
+        }
+    )
+    return evt_variables
+
+def add_truth_matching_functors(
+    options,
+    variables,
+    data,
+    hierarchy=True,
+    kinematics=True,
+    vertex_info=True,
+    bkgcat=False,
+    prompt_info=True,
+):
+    """Add MC truth matching functors to an existing FunctorCollection.
+
+    Args:
+        options (DaVinci.Options): DaVinci options object.
+        variables (FunctorCollection): FunctorCollection to modify
+        data: data handle
+        hierarchy (bool): Add MCHierarchy info (default True)
+        kinematics (bool): Add MCKinematics info (default True)
+        vertex_info (bool): Add MCVertexInfo (default True)
+        bkgcat (bool): Add TRUEKEY and BKGCAT functors - intended for composite particles (default False)
+        prompt_info (bool): Add MCPromptDecay info (default True)
+
+    Returns:
+        FunctorCollection: modified FunctorCollection with truth matched variables.
+    """
+    assert (
+        options.simulation
+    ), "options.simulation is set to False - it doesn't make sense to add truth matching."
+
+    MCTRUTH = MCTruthAndBkgCat(data)
+
+    if bkgcat:
+        variables += FunctorCollection(
+            {
+                # Important note: specify an invalid value for integer functors if there exists no truth info.
+                #                 The invalid value for floating point functors is set to nan.
+                "TRUEKEY": F.VALUE_OR(-1) @ MCTRUTH(F.OBJECT_KEY),
+                "BKGCAT": MCTRUTH.BkgCat,
+            }
+        )
+    if hierarchy:
+        variables += MCHierarchy(mctruth_alg=MCTRUTH)  # likely to change again!!
+    if vertex_info:
+        variables += MCVertexInfo(mctruth_alg=MCTRUTH)  # likely to change again!!
+    if kinematics:
+        variables += MCKinematics(mctruth_alg=MCTRUTH)  # likely to change again!!
+    if prompt_info:
+        variables += MCPromptDecay(mctruth_alg=MCTRUTH)
+
+    return variables
+
+### DTF variables ###
+
+def make_basic_dtf_variables(options, pvs, data, DTF=None, pv_constraint=False, mass_constraint=False, particle_name=""):
+    variables = (
+        FunctorCollection(
+            {
+                "TRCHI2DOF": F.CHI2DOF @ F.TRACK,
+                "ETA": F.ETA,
+                "PHI": F.PHI,
+                "TRGHOSTPROB": F.GHOSTPROB,
+                "TX"          : F.TX,
+                "TY"          : F.TY,
+                "KEY"         : F.VALUE_OR(-1) @ F.OBJECT_KEY @ F.TRACK,
+                "CTB"         : F.POSITION @ F.CLOSESTTOBEAM @ F.TRACK,
+                "TRACKPT": F.TRACK_PT,
+                "TRACKHISTORY": F.VALUE_OR(-1) @ F.TRACKHISTORY @ F.TRACK,
+                "QOVERP": F.QOVERP @ F.TRACK,
+                "NDOF": F.VALUE_OR(-1) @ F.NDOF @ F.TRACK,
+                "POSITION_STATEAT_FirstMeasurement_X"     : F.POSITION_X   @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_Y"     : F.POSITION_Y   @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_Z"     : F.POSITION_Z   @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_TX"    : F.X_COORDINATE @ F.SLOPES @ F.STATE_AT("FirstMeasurement") @ F.TRACK ,
+                "POSITION_STATEAT_FirstMeasurement_TY"    : F.Y_COORDINATE @ F.SLOPES @ F.STATE_AT("FirstMeasurement") @ F.TRACK , 
+            }
+        )
+        + Kinematics()
+    )
+
+    if(mass_constraint):
+        if(pv_constraint): # MASS + PV
+            dtf_variables_mass_pv = FunctorCollection({
+                        'DTF_PV_M'+ particle_name + '_' + k: DTF(v)
+                        for k, v in variables.get_thor_functors().items()
+                    })
+            return dtf_variables_mass_pv
+        else: # MASS
+            dtf_variables_mass = FunctorCollection(
+                {'DTF_M'+ particle_name + '_' + k: DTF(v)
+                 for k, v in variables.get_thor_functors().items()})
+        return dtf_variables_mass
+
+    elif(pv_constraint): # PV
+        dtf_variables_pv = FunctorCollection({
+                'DTF_PV_' + k: DTF(v)
+                for k, v in variables.get_thor_functors().items()
+            })
+        return dtf_variables_pv
+
+    else: # NO MASS/PV
+        dtf_variables = FunctorCollection(
+            {'DTF_' + k: DTF(v)
+             for k, v in variables.get_thor_functors().items()})
+        return dtf_variables
+
+def make_composite_dtf_variables(options, pvs, data, DTF=None, pv_constraint=False, mass_constraint=False, particle_name=""):
+    variables = (
+        FunctorCollection(
+            {
+                "MAXPT": F.MAX(F.PT),
+                "MINPT": F.MIN(F.PT),
+                "SUMPT": F.SUM(F.PT),
+                "MAXP": F.MAX(F.P),
+                "MINP": F.MIN(F.P),
+                "BPVDIRA": F.BPVDIRA(pvs),
+                "VCHI2DOF": F.CHI2DOF, #CHI2VXNDOF
+                "BPVFDCHI2": F.BPVFDCHI2(pvs),
+                "BPVFD": F.BPVFD(pvs),
+                "BPVVDRHO": F.BPVVDRHO(pvs),
+                "BPVVDZ": F.BPVVDZ(pvs),
+                "BPVIPCHI2": F.BPVIPCHI2(pvs),
+                "BPVIP": F.BPVIP(pvs),
+                "LOGBPVIPCHI2": log(F.BPVIPCHI2(pvs)),
+                "BPVLTIME": F.BPVLTIME(pvs),
+                "MAXBPVIPCHI2": F.MAX(F.BPVIPCHI2(pvs)), #MAX_
+                "MINBPVIPCHI2": F.MIN(F.BPVIPCHI2(pvs)),
+                "MAXBPVIP": F.MAX(F.BPVIP(pvs)),
+                "MINBPVIP": F.MIN(F.BPVIP(pvs)),
+                "ETA": F.ETA,
+                "PHI": F.PHI,
+                "END_VX": F.END_VX,
+                "END_VY": F.END_VY,
+                "END_VZ": F.END_VZ,
+                "END_VX_ERR":F.SQRT @ F.CALL(0,0) @ F.POS_COV_MATRIX @ F.ENDVERTEX,
+                "END_VY_ERR":F.SQRT @ F.CALL(1,1) @ F.POS_COV_MATRIX @ F.ENDVERTEX,
+                "END_VZ_ERR":F.SQRT @ F.CALL(2,2) @ F.POS_COV_MATRIX @ F.ENDVERTEX,
+                "BPVX": F.BPVX(pvs),
+                "BPVY": F.BPVY(pvs),
+                "BPVZ": F.BPVZ(pvs),
+                "BPVX_ERR": F.SQRT @ F.CALL(0,0) @ F.POS_COV_MATRIX @ F.BPV(pvs),
+                "BPVY_ERR": F.SQRT @ F.CALL(1,1) @ F.POS_COV_MATRIX @ F.BPV(pvs),
+                "BPVZ_ERR": F.SQRT @ F.CALL(2,2) @ F.POS_COV_MATRIX @ F.BPV(pvs),
+                "ALLPVFD"     : F.ALLPV_FD(pvs),
+                "ALLPVIP"     : F.ALLPV_IP(pvs),
+                "DM": F.MASS - F.CHILD(1, F.MASS)
+                
+            }
+        )
+        + Kinematics()
+    )
+
+    addstring = "DTF"
+    if(pv_constraint):
+            addstring += '_PV'
+    if(mass_constraint):
+            addstring += '_M'
+    addstring += particle_name
+
+    DTF_chi2ndof = FunctorCollection(
+            {
+                addstring+"_DTFCHI2": DTF.CHI2,
+                addstring+"_DTFNDOF": DTF.NDOF,
+                addstring+"_CTAU": DTF.CTAU,
+                addstring+"_CTAUERR": DTF.CTAUERR,
+                addstring+"_MERR": DTF.MASSERR,
+            }
+    )
+
+    if(mass_constraint):
+        if(pv_constraint): # MASS + PV
+            dtf_variables_mass_pv = FunctorCollection({
+                        'DTF_PV_M'+ particle_name + '_' + k: DTF(v)
+                        for k, v in variables.get_thor_functors().items()
+                    })
+            return dtf_variables_mass_pv+DTF_chi2ndof
+        else: # MASS
+            dtf_variables_mass = FunctorCollection(
+                {'DTF_M'+ particle_name + '_' + k: DTF(v)
+                 for k, v in variables.get_thor_functors().items()})
+        return dtf_variables_mass+DTF_chi2ndof
+
+    elif(pv_constraint): # PV
+        dtf_variables_pv = FunctorCollection({
+                'DTF_PV_' + k: DTF(v)
+                for k, v in variables.get_thor_functors().items()
+            })
+        return dtf_variables_pv+DTF_chi2ndof
+
+    else: # NO MASS/PV
+        dtf_variables = FunctorCollection(
+            {'DTF_' + k: DTF(v)
+             for k, v in variables.get_thor_functors().items()})
+        return dtf_variables+DTF_chi2ndof
+
+def make_composite_dtf_variables_3body(options, pvs, data, DTF=None, pv_constraint=False, mass_constraint=False, particle_name=""):
+    variables = (
+        FunctorCollection(
+            {
+                "M13": F.SUBCOMB(Functor=F.MASS, Indices=(1, 3)),
+                "M23": F.SUBCOMB(Functor=F.MASS, Indices=(2, 3)),
+            }
+        )
+    )
+
+    if(mass_constraint):
+        if(pv_constraint): # MASS + PV
+            dtf_variables_mass_pv = FunctorCollection({
+                        'DTF_PV_M'+ particle_name + '_' + k: DTF(v)
+                        for k, v in variables.get_thor_functors().items()
+                    })
+            return dtf_variables_mass_pv+make_composite_dtf_variables(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+        else: # MASS
+            dtf_variables_mass = FunctorCollection(
+                {'DTF_M'+ particle_name + '_' + k: DTF(v)
+                 for k, v in variables.get_thor_functors().items()})
+        return dtf_variables_mass+make_composite_dtf_variables(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+
+    elif(pv_constraint): # PV
+        dtf_variables_pv = FunctorCollection({
+                'DTF_PV_' + k: DTF(v)
+                for k, v in variables.get_thor_functors().items()
+            })
+        return dtf_variables_pv+make_composite_dtf_variables(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+
+    else: # NO MASS/PV
+        dtf_variables = FunctorCollection(
+            {'DTF_' + k: DTF(v)
+             for k, v in variables.get_thor_functors().items()})
+
+    return dtf_variables+make_composite_dtf_variables(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+
+def make_composite_dtf_variables_4body(options, pvs, data, DTF=None, pv_constraint=False, mass_constraint=False, particle_name=""):
+    variables = (
+        FunctorCollection(
+            {
+                "M14": F.SUBCOMB(Functor=F.MASS, Indices=(1, 4)),
+                "M24": F.SUBCOMB(Functor=F.MASS, Indices=(2, 4)),
+                "M34": F.SUBCOMB(Functor=F.MASS, Indices=(3, 4)),
+            }
+        )
+    )
+
+    if(mass_constraint):
+        if(pv_constraint): # MASS + PV
+            dtf_variables_mass_pv = FunctorCollection({
+                        'DTF_PV_M'+ particle_name + '_' + k: DTF(v)
+                        for k, v in variables.get_thor_functors().items()
+                    })
+            return dtf_variables_mass_pv+make_composite_dtf_variables_3body(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+        else: # MASS
+            dtf_variables_mass = FunctorCollection(
+                {'DTF_M'+ particle_name + '_' + k: DTF(v)
+                 for k, v in variables.get_thor_functors().items()})
+        return dtf_variables_mass+make_composite_dtf_variables_3body(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+
+    elif(pv_constraint): # PV
+        dtf_variables_pv = FunctorCollection({
+                'DTF_PV_' + k: DTF(v)
+                for k, v in variables.get_thor_functors().items()
+            })
+        return dtf_variables_pv+make_composite_dtf_variables_3body(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+
+    else: # NO MASS/PV
+        dtf_variables = FunctorCollection(
+            {'DTF_' + k: DTF(v)
+             for k, v in variables.get_thor_functors().items()})
+    
+        return dtf_variables+make_composite_dtf_variables_3body(options,pvs,data,DTF,pv_constraint,mass_constraint, particle_name)
+
+
+def make_MC_basic_variables():
+    MC_MOTHER_ID = lambda gen: F.VALUE_OR(0) @ F.MC_MOTHER(gen, F.PARTICLE_ID)
+    MC_MOTHER_KEY = lambda gen: F.VALUE_OR(-1) @ F.MC_MOTHER(gen, F.OBJECT_KEY)
+
+    variables = (
+        FunctorCollection(
+            {
+                "OBJECT_KEY": F.OBJECT_KEY,
+                "ETA": F.ETA,
+                "PHI": F.PHI,
+                "ORIGIN_VX": F.ORIGIN_VX,
+                "ORIGIN_VY": F.ORIGIN_VY,
+                "ORIGIN_VZ": F.ORIGIN_VZ,
+                "FOURMOMENTUM": F.FOURMOMENTUM,
+                "TRUEID": F.PARTICLE_ID,
+                "MC_MOTHER_ID": MC_MOTHER_ID(1),
+                "MC_MOTHER_KEY": MC_MOTHER_KEY(1),
+                "MC_GD_MOTHER_ID": MC_MOTHER_ID(2),
+                "MC_GD_MOTHER_KEY": MC_MOTHER_KEY(2),
+                "MC_GD_GD_MOTHER_ID": MC_MOTHER_ID(3),
+                "MC_GD_GD_MOTHER_KEY": MC_MOTHER_KEY(3),
+            }
+        )
+        + Kinematics()
+    )
+
+
+    from PyConf.reading import get_mc_particles, get_pp2mcp_relations, get_mc_track_info
+    MC_data = get_mc_particles("/Event/HLT2/MC/Particles")
+    relations_charged = get_pp2mcp_relations("/Event/HLT2/Relations/ChargedPP2MCP")
+    relations_neutral = get_pp2mcp_relations("/Event/HLT2/Relations/NeutralPP2MCP")
+    mcreconstructed_all = MCReconstructed(
+        input_mcparticles=MC_data,
+        use_best_mcmatch=True,
+        relations_charged=relations_charged,
+        relations_neutral=relations_neutral,
+    )
+    mcreconstructible_all = MCReconstructible(
+        input_mctrackinfo=get_mc_track_info()
+    )
+    import FunTuple.functorcollections as FC
+    vars_reconstructed = FC.MCReconstructed(mcreconstructed_alg=mcreconstructed_all, extra_info=True)
+    vars_reconstructible = FC.MCReconstructible(mcreconstructible_alg=mcreconstructible_all, extra_info=True)
+
+
+    return variables+vars_reconstructed+vars_reconstructible
+
+def make_MC_composite_variables():
+    MC_MOTHER_ID = lambda gen: F.VALUE_OR(0) @ F.MC_MOTHER(gen, F.PARTICLE_ID)
+    MC_MOTHER_KEY = lambda gen: F.VALUE_OR(-1) @ F.MC_MOTHER(gen, F.OBJECT_KEY)
+
+    variables = (
+        FunctorCollection(
+            {
+                "END_VX": F.END_VX,
+                "END_VY": F.END_VY,
+                "END_VZ": F.END_VZ,
+                "LTIME": F.MC_LIFETIME,
+                "OBJECT_KEY": F.OBJECT_KEY,
+                "ETA": F.ETA,
+                "PHI": F.PHI,
+                "ORIGIN_VX": F.ORIGIN_VX,
+                "ORIGIN_VY": F.ORIGIN_VY,
+                "ORIGIN_VZ": F.ORIGIN_VZ,
+                "FOURMOMENTUM": F.FOURMOMENTUM,
+                "TRUEID": F.PARTICLE_ID,
+                "MC_MOTHER_ID": MC_MOTHER_ID(1),
+                "MC_MOTHER_KEY": MC_MOTHER_KEY(1),
+                "MC_GD_MOTHER_ID": MC_MOTHER_ID(2),
+                "MC_GD_MOTHER_KEY": MC_MOTHER_KEY(2),
+                "MC_GD_GD_MOTHER_ID": MC_MOTHER_ID(3),
+                "MC_GD_GD_MOTHER_KEY": MC_MOTHER_KEY(3),
+            }
+        )
+        + Kinematics()
+        + MCPromptDecay()
+        )
+
+    return variables
+
+def make_MC_event_variables(mc_header):
+    # define event level variables
+    evt_variables = EventInfo()
+    evt_variables += MCPrimaries(mc_header=mc_header)
+
+    return evt_variables
+
+def make_top_isolation_variables(hlt2_line, input_data, locations = ["LongTrackIso","NeutralIso"]):
+    from PyConf.reading import get_particles
+    from IsolationTools import VertexAndConeIsolation
+
+    possible_charm_locations = ["LongTrackIso", "TTrackIso", "DownstreamTrackIso", "UpstreamTrackIso", "NeutralIso", "PizIso"]
+    coneangles = [0.25,0.5,1.,1.5,2.]
+
+    count = 0
+    for location in locations:
+        extra_particles = get_particles(f"/Event/HLT2/{hlt2_line}/{location}/Particles")
+
+        for coneangle in coneangles:
+            top_RTAlg = VertexAndConeIsolation(
+                name=location+"_"+str(coneangle),
+                reference_particles=input_data,
+                related_particles=extra_particles,
+                cut=F.DR2<coneangle)
+
+            if count == 0:
+                top_iso_variables = ParticleIsolation(isolation_alg=top_RTAlg,array_indx_name='indx')
+            else:
+                top_iso_variables += ParticleIsolation(isolation_alg=top_RTAlg,array_indx_name='indx')
+            count += 1
+
+
+    return top_iso_variables
+
+def make_basic_isolation_variables(hlt2_line, input_data, locations = ["LongTrackIso","NeutralIso"]):
+    from PyConf.reading import get_particles
+    from IsolationTools import VertexAndConeIsolation
+    from PyConf.Algorithms import ThOrParticleSelection
+
+    basic_code = (F.FILTER(F.ALL) @ F.GET_ALL_BASICS())
+    basic_particles = ThOrParticleSelection(InputParticles=input_data, Functor=basic_code).OutputSelection
+
+    possible_charm_locations = ["LongTrackIso", "TTrackIso", "DownstreamTrackIso", "UpstreamTrackIso", "NeutralIso", "PizIso"]
+    coneangles = [0.25,0.5,1.,1.5,2.]
+
+    count = 0
+    for location in locations:
+        extra_particles = get_particles(f"/Event/HLT2/{hlt2_line}/{location}/Particles")
+
+        for coneangle in coneangles:
+            basic_RTAlg = VertexAndConeIsolation(
+                name=location+"_"+str(coneangle),
+                reference_particles=basic_particles,
+                related_particles=extra_particles,
+                cut=F.DR2<coneangle)
+
+            if count == 0:
+                basic_iso_variables = ParticleIsolation(isolation_alg=basic_RTAlg,array_indx_name='indx')
+            else:
+                basic_iso_variables += ParticleIsolation(isolation_alg=basic_RTAlg,array_indx_name='indx')
+            count += 1
+
+
+    return basic_iso_variables
+
+
+def make_intermediate_isolation_variables(hlt2_line, input_data, locations = ["LongTrackIso","NeutralIso"], composite_ID = "J/psi(1S)"):
+    from PyConf.reading import get_particles
+    from IsolationTools import VertexAndConeIsolation
+    from PyConf.Algorithms import ThOrParticleSelection
+
+    intermediate_code = F.FILTER(F.IS_ABS_ID(composite_ID)) @ F.GET_ALL_DESCENDANTS()
+    intermediate_particles = ThOrParticleSelection(InputParticles=input_data, Functor=intermediate_code).OutputSelection
+
+    possible_charm_locations = ["LongTrackIso", "TTrackIso", "DownstreamTrackIso", "UpstreamTrackIso", "NeutralIso", "PizIso"]
+    coneangles = [0.25,0.5,1.,1.5,2.]
+
+    count = 0
+    for location in locations:
+        extra_particles = get_particles(f"/Event/HLT2/{hlt2_line}/{location}/Particles")
+
+        for coneangle in coneangles:
+            intermediate_RTAlg = VertexAndConeIsolation(
+                name=location+"_"+str(coneangle),
+                reference_particles=intermediate_particles,
+                related_particles=extra_particles,
+                cut=F.DR2<coneangle)
+
+            if count == 0:
+                intermediate_iso_variables = ParticleIsolation(isolation_alg=intermediate_RTAlg,array_indx_name='indx')
+            else:
+                intermediate_iso_variables += ParticleIsolation(isolation_alg=intermediate_RTAlg,array_indx_name='indx')
+            count += 1
+
+
+    return intermediate_iso_variables
+
+def photon_isolation_variables(input_data, extra_photons):
+
+    gamma_data = ThOrParticleSelection(
+        InputParticles=input_data, Functor=FILTER_TREE("gamma")
+    ).OutputSelection
+
+    ###### test isolation
+    gamma_pt_nc_isoAlg = WeightedRelTableAlg(
+        ReferenceParticles=gamma_data,
+        InputCandidates=extra_photons,
+        Cut=in_range(0. * MeV, F.COMB_MASS(), 500. * MeV))
+    
+    gamma_pt_nc_isoAlg_Rels = gamma_pt_nc_isoAlg.OutputRelations
+    
+    prefix = "HEAD_NC_gammaIsoAll"
+    extra_variables=  FunctorCollection(
+        {
+            f"{prefix}_P":   F.MAP_INPUT_ARRAY(Functor=F.P, Relations=gamma_pt_nc_isoAlg_Rels),
+            f"{prefix}_PT":  F.MAP_INPUT_ARRAY(Functor=F.PT, Relations=gamma_pt_nc_isoAlg_Rels),
+            f"{prefix}_PX":  F.MAP_INPUT_ARRAY(Functor=F.PX, Relations=gamma_pt_nc_isoAlg_Rels),
+            f"{prefix}_PY":  F.MAP_INPUT_ARRAY(Functor=F.PY, Relations=gamma_pt_nc_isoAlg_Rels),
+            f"{prefix}_PZ":  F.MAP_INPUT_ARRAY(Functor=F.PZ, Relations=gamma_pt_nc_isoAlg_Rels),
+            f"{prefix}_IS_PHOTON":  F.MAP_INPUT_ARRAY(Functor=F.VALUE_OR(F.NaN) @ F.VALUE_OR(F.NaN)@ F.IS_PHOTON, Relations=gamma_pt_nc_isoAlg_Rels),
+            f"{prefix}_IS_NOT_H":  F.MAP_INPUT_ARRAY(Functor=F.VALUE_OR(F.NaN) @ F.IS_NOT_H, Relations=gamma_pt_nc_isoAlg_Rels),
+            f"{prefix}_CALO_NEUTRAL_SHOWER_SHAPE":  F.MAP_INPUT_ARRAY(Functor=F.VALUE_OR(F.NaN) @ F.CALO_NEUTRAL_SHOWER_SHAPE, Relations=gamma_pt_nc_isoAlg_Rels),
+        })
+
+    coneangles =      [ 0.1,  0.20,  0.30, 0.4,  0.5,   1.     ]
+    coneangles_name = ["10",  "20",  "30", "40",  '50', '100' ]
+
+    for i, coneangle in enumerate(coneangles):
+        gamma_pt_nc_angle_isoAlg = WeightedRelTableAlg(
+        ReferenceParticles=gamma_data,
+        InputCandidates=extra_photons,
+        Cut=F.require_all(F.DR2<coneangle, in_range(0. * MeV, F.COMB_MASS(), 500. * MeV)))
+        gamma_pt_nc_angle_isoAlg_Rels = gamma_pt_nc_angle_isoAlg.OutputRelations
+        prefix = f"HEAD_NC_gammaIso{coneangles_name[i]}"
+        extra_variables +=  FunctorCollection(
+            {
+                f"{prefix}_MAXP":   F.VALUE_OR(0) @ F.P @ F.TO @ F.ENTRY_WITH_MAX_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MAXPX":   F.VALUE_OR(0) @ F.PX @ F.TO @ F.ENTRY_WITH_MAX_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MAXPY":   F.VALUE_OR(0) @ F.PY @ F.TO @ F.ENTRY_WITH_MAX_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MAXPZ":   F.VALUE_OR(0) @ F.PZ @ F.TO @ F.ENTRY_WITH_MAX_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MAXPT":   F.VALUE_OR(0) @ F.PT @ F.TO @ F.ENTRY_WITH_MAX_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MINP":   F.VALUE_OR(0) @ F.P @ F.TO @ F.ENTRY_WITH_MIN_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MINPX":   F.VALUE_OR(0) @ F.PX @ F.TO @ F.ENTRY_WITH_MIN_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MINPY":   F.VALUE_OR(0) @ F.PY @ F.TO @ F.ENTRY_WITH_MIN_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MINPZ":   F.VALUE_OR(0) @ F.PZ @ F.TO @ F.ENTRY_WITH_MIN_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_MINPT":   F.VALUE_OR(0) @ F.PT @ F.TO @ F.ENTRY_WITH_MIN_REL_VALUE_OF(F.P@F.TO@F.FORWARDARG0).bind( F.RELATIONS.bind(F.TES(gamma_pt_nc_angle_isoAlg_Rels), F.FORWARDARGS), F.FORWARDARGS),
+                f"{prefix}_P":   F.VALUE_OR(0) @ F.MAP_INPUT_ARRAY(Functor=F.P, Relations=gamma_pt_nc_angle_isoAlg_Rels),
+                f"{prefix}_PT":   F.VALUE_OR(0) @ F.MAP_INPUT_ARRAY(Functor=F.PT, Relations=gamma_pt_nc_angle_isoAlg_Rels),
+                f"{prefix}_PX":  F.VALUE_OR(0) @ F.MAP_INPUT_ARRAY(Functor=F.PX, Relations=gamma_pt_nc_angle_isoAlg_Rels),
+                f"{prefix}_PY":  F.VALUE_OR(0) @ F.MAP_INPUT_ARRAY(Functor=F.PY, Relations=gamma_pt_nc_angle_isoAlg_Rels),
+                f"{prefix}_PZ":  F.VALUE_OR(0) @ F.MAP_INPUT_ARRAY(Functor=F.PZ, Relations=gamma_pt_nc_angle_isoAlg_Rels),
+                #f"{prefix}_DR2":   F.VALUE_OR(0) @ F.MAP_INPUT_ARRAY(Functor=F.DR2, Relations=gamma_pt_nc_isoAlg_Rels),
+                #f"{prefix}_IS_PHOTON":   F.MAP_INPUT_ARRAY(Functor=F.IS_PHOTON, Relations=gamma_pt_nc_isoAlg_Rels),
+                #f"{prefix}_IS_NOT_H":   F.VALUE_OR(0)@  F.MAP_INPUT_ARRAY(Functor=F.IS_NOT_H, Relations=gamma_pt_nc_isoAlg_Rels),
+                #f"{prefix}_CALO_NEUTRAL_SHOWER_SHAPE":   F.VALUE_OR(0)@  F.MAP_INPUT_ARRAY(Functor=F.CALO_NEUTRAL_SHOWER_SHAPE, Relations=gamma_pt_nc_isoAlg_Rels),
+            }
+            )
+    return extra_variables
\ No newline at end of file