diff --git a/addons/ONNXRuntime/python/jetFlavourHelper.py b/addons/ONNXRuntime/python/jetFlavourHelper.py index a8027832352..95dfdc3c09f 100644 --- a/addons/ONNXRuntime/python/jetFlavourHelper.py +++ b/addons/ONNXRuntime/python/jetFlavourHelper.py @@ -5,7 +5,17 @@ ROOT.gROOT.SetBatch(True) class JetFlavourHelper: - def __init__(self, coll, jet, jetc, tag=""): + ''' + NOTE: (May 2025) Once the full sim tagger is retrained on the new naming convention (see https://github.com/key4hep/k4MLJetTagger?tab=readme-ov-file#open-problems--further-work), the names defined here must be altered. Then, they will also nicely match the namings in ReconstructedParticle2Track. + ''' + def __init__(self, coll, jet, jetc, tag="", sim_type="fast"): + ''' + sim_type: fast or full + ''' + # check if sim_type is valid + if sim_type not in ["fast", "full"]: + print("ERROR: sim_type must be either 'fast' or 'full'") + sys.exit() self.jet = jet self.const = jetc @@ -13,25 +23,32 @@ def __init__(self, coll, jet, jetc, tag=""): self.tag = tag if tag != "": self.tag = "_{}".format(tag) + self.sim_type = sim_type self.particle = coll["GenParticles"] self.pfcand = coll["PFParticles"] self.pftrack = coll["PFTracks"] self.pfphoton = coll["PFPhotons"] self.pfnh = coll["PFNeutralHadrons"] - self.trackstate = coll["TrackState"] + + self.trackstate = coll["TrackStates"] + self.tracks = coll["Tracks"] + self.trackerhits = coll["TrackerHits"] self.calohits = coll["CalorimeterHits"] - self.dndx = coll["dNdx"] self.l = coll["PathLength"] - self.bz = coll["Bz"] - + if sim_type == "fast": + self.dndx = coll["dNdx"] + self.bz = coll["Bz"] + elif sim_type == "full": + self.bz = "2.0" # CLD #FIXME: this should be read from the geometry + self.dndx = None + self.primvertex = coll["PV"] self.definition = dict() - # ===== VERTEX - # MC primary vertex - self.definition["pv{}".format(self.tag)] = "FCCAnalyses::MCParticle::get_EventPrimaryVertexP4()( {} )".format( - self.particle + # ===== VERTEX (reconstructed) + self.definition["pv{}".format(self.tag)] = "JetConstituentsUtils::get_primary_vertex({})".format( + self.primvertex ) # build jet constituents lists @@ -69,105 +86,114 @@ def __init__(self, coll, jet, jetc, tag=""): self.definition["pfcand_phirel{}".format(self.tag)] = "JetConstituentsUtils::get_phirel_cluster({}, {})".format( jet, self.const ) + if self.sim_type == "fast": + self.definition["Bz{}".format(self.tag)] = "{}[0]".format(self.bz) + self.definition[ + "pfcand_dndx{}".format(self.tag) + ] = "JetConstituentsUtils::get_dndx({}, {}, {}, pfcand_isChargedHad{})".format( + self.const, self.dndx, self.pftrack, self.tag + ) - self.definition[ - "pfcand_dndx{}".format(self.tag) - ] = "JetConstituentsUtils::get_dndx({}, {}, {}, pfcand_isChargedHad{})".format( - self.const, self.dndx, self.pftrack, self.tag - ) + self.definition[ + "pfcand_mtof{}".format(self.tag) + ] = "JetConstituentsUtils::get_mtof({}, {}, {}, {}, {}, {}, {}, pv{})".format( + self.const, self.l, self.pftrack, self.trackerhits, self.pfphoton, self.pfnh, self.calohits, self.tag + ) + elif self.sim_type == "full": + self.definition["Bz{}".format(self.tag)] = self.bz + # fill the dNdx and mtof variables with 0 + self.definition[ + "pfcand_dndx{}".format(self.tag) + ] = "JetConstituentsUtils::get_dndx_dummy({})".format(self.const) - self.definition[ - "pfcand_mtof{}".format(self.tag) - ] = "JetConstituentsUtils::get_mtof({}, {}, {}, {}, {}, {}, {}, pv{})".format( - self.const, self.l, self.pftrack, self.trackerhits, self.pfphoton, self.pfnh, self.calohits, self.tag - ) + self.definition[ + "pfcand_mtof{}".format(self.tag) + ] = "JetConstituentsUtils::get_mtof_dummy({})".format(self.const) - self.definition["Bz{}".format(self.tag)] = "{}[0]".format(self.bz) - - self.definition[ - "pfcand_dxy{}".format(self.tag) - ] = "JetConstituentsUtils::XPtoPar_dxy({}, {}, pv{}, Bz{})".format( - self.const, self.trackstate, self.tag, self.tag + self.definition["pfcand_dxy{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_dxy({}, {}, {}, pv{}, Bz{})".format( + self.const, self.trackstate, self.tracks, self.tag, self.tag ) - self.definition["pfcand_dz{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_dz({}, {}, pv{}, Bz{})".format( - self.const, self.trackstate, self.tag, self.tag + self.definition["pfcand_dz{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_dz({}, {}, {}, pv{}, Bz{})".format( + self.const, self.trackstate, self.tracks, self.tag, self.tag ) - self.definition[ - "pfcand_phi0{}".format(self.tag) - ] = "JetConstituentsUtils::XPtoPar_phi({}, {}, pv{}, Bz{})".format( - self.const, self.trackstate, self.tag, self.tag + self.definition["pfcand_phi0{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_phi({}, {}, {}, pv{}, Bz{})".format( + self.const, self.trackstate, self.tracks, self.tag, self.tag ) - self.definition["pfcand_C{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_C({}, {}, Bz{})".format( - self.const, self.trackstate, self.tag + self.definition["pfcand_C{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_C({}, {}, {}, Bz{})".format( + self.const, self.trackstate, self.tracks, self.tag ) - self.definition["pfcand_ct{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_ct({}, {}, Bz{})".format( - self.const, self.trackstate, self.tag + self.definition["pfcand_ct{}".format(self.tag)] = "JetConstituentsUtils::XPtoPar_ct({}, {}, {}, Bz{})".format( + self.const, self.trackstate, self.tracks, self.tag ) - self.definition["pfcand_dptdpt{}".format(self.tag)] = "JetConstituentsUtils::get_omega_cov({}, {})".format( - self.const, self.trackstate + # covariance matrix (fixed track state problem) + + self.definition["pfcand_dptdpt{}".format(self.tag)] = "JetConstituentsUtils::get_omega_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition["pfcand_dxydxy{}".format(self.tag)] = "JetConstituentsUtils::get_d0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_dxydxy{}".format(self.tag)] = "JetConstituentsUtils::get_d0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition["pfcand_dzdz{}".format(self.tag)] = "JetConstituentsUtils::get_z0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_dzdz{}".format(self.tag)] = "JetConstituentsUtils::get_z0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition["pfcand_dphidphi{}".format(self.tag)] = "JetConstituentsUtils::get_phi0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_dphidphi{}".format(self.tag)] = "JetConstituentsUtils::get_phi0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition[ - "pfcand_detadeta{}".format(self.tag) - ] = "JetConstituentsUtils::get_tanlambda_cov({}, {})".format(self.const, self.trackstate) + self.definition["pfcand_detadeta{}".format(self.tag)] = "JetConstituentsUtils::get_tanlambda_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate + ) - self.definition["pfcand_dxydz{}".format(self.tag)] = "JetConstituentsUtils::get_d0_z0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_dxydz{}".format(self.tag)] = "JetConstituentsUtils::get_d0_z0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition["pfcand_dphidxy{}".format(self.tag)] = "JetConstituentsUtils::get_phi0_d0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_dphidxy{}".format(self.tag)] = "JetConstituentsUtils::get_phi0_d0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition["pfcand_phidz{}".format(self.tag)] = "JetConstituentsUtils::get_phi0_z0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_phidz{}".format(self.tag)] = "JetConstituentsUtils::get_phi0_z0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition[ - "pfcand_phictgtheta{}".format(self.tag) - ] = "JetConstituentsUtils::get_tanlambda_phi0_cov({}, {})".format(self.const, self.trackstate) + self.definition["pfcand_phictgtheta{}".format(self.tag)] = "JetConstituentsUtils::get_tanlambda_phi0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate + ) - self.definition[ - "pfcand_dxyctgtheta{}".format(self.tag) - ] = "JetConstituentsUtils::get_tanlambda_d0_cov({}, {})".format(self.const, self.trackstate) + self.definition["pfcand_dxyctgtheta{}".format(self.tag)] = "JetConstituentsUtils::get_tanlambda_d0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate + ) - self.definition[ - "pfcand_dlambdadz{}".format(self.tag) - ] = "JetConstituentsUtils::get_tanlambda_z0_cov({}, {})".format(self.const, self.trackstate) + self.definition["pfcand_dlambdadz{}".format(self.tag)] = "JetConstituentsUtils::get_tanlambda_z0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate + ) - self.definition[ - "pfcand_cctgtheta{}".format(self.tag) - ] = "JetConstituentsUtils::get_omega_tanlambda_cov({}, {})".format(self.const, self.trackstate) + self.definition["pfcand_cctgtheta{}".format(self.tag)] = "JetConstituentsUtils::get_omega_tanlambda_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate + ) - self.definition["pfcand_phic{}".format(self.tag)] = "JetConstituentsUtils::get_omega_phi0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_phic{}".format(self.tag)] = "JetConstituentsUtils::get_omega_phi0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition["pfcand_dxyc{}".format(self.tag)] = "JetConstituentsUtils::get_omega_d0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_dxyc{}".format(self.tag)] = "JetConstituentsUtils::get_omega_d0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) - self.definition["pfcand_cdz{}".format(self.tag)] = "JetConstituentsUtils::get_omega_z0_cov({}, {})".format( - self.const, self.trackstate + self.definition["pfcand_cdz{}".format(self.tag)] = "JetConstituentsUtils::get_omega_z0_cov({}, {}, {})".format( + self.const, self.tracks, self.trackstate ) + # impact parameters + self.definition[ "pfcand_btagSip2dVal{}".format(self.tag) ] = "JetConstituentsUtils::get_Sip2dVal_clusterV({}, pfcand_dxy{}, pfcand_phi0{}, Bz{})".format( @@ -176,7 +202,7 @@ def __init__(self, coll, jet, jetc, tag=""): self.definition[ "pfcand_btagSip2dSig{}".format(self.tag) - ] = "JetConstituentsUtils::get_Sip2dSig(pfcand_btagSip2dVal{}, pfcand_dxydxy{})".format(self.tag, self.tag) + ] = 'JetConstituentsUtils::get_Sip2dSig(pfcand_btagSip2dVal{}, pfcand_dxydxy{}, "{}")'.format(self.tag, self.tag, self.sim_type) self.definition[ "pfcand_btagSip3dVal{}".format(self.tag) @@ -186,8 +212,8 @@ def __init__(self, coll, jet, jetc, tag=""): self.definition[ "pfcand_btagSip3dSig{}".format(self.tag) - ] = "JetConstituentsUtils::get_Sip3dSig(pfcand_btagSip3dVal{}, pfcand_dxydxy{}, pfcand_dzdz{})".format( - self.tag, self.tag, self.tag + ] = 'JetConstituentsUtils::get_Sip3dSig(pfcand_btagSip3dVal{}, pfcand_dxydxy{}, pfcand_dzdz{}, "{}")'.format( + self.tag, self.tag, self.tag, self.sim_type ) self.definition[ @@ -198,10 +224,12 @@ def __init__(self, coll, jet, jetc, tag=""): self.definition[ "pfcand_btagJetDistSig{}".format(self.tag) - ] = "JetConstituentsUtils::get_JetDistSig(pfcand_btagJetDistVal{}, pfcand_dxydxy{}, pfcand_dzdz{})".format( - self.tag, self.tag, self.tag + ] = 'JetConstituentsUtils::get_JetDistSig(pfcand_btagJetDistVal{}, pfcand_dxydxy{}, pfcand_dzdz{}, "{}")'.format( + self.tag, self.tag, self.tag, self.sim_type ) + # count number of particles in the jet + self.definition["jet_nmu{}".format(self.tag)] = "JetConstituentsUtils::count_type(pfcand_isMu{})".format( self.tag ) @@ -249,13 +277,16 @@ def inference(self, jsonCfg, onnxCfg, df): # convert to tuple initvars = tuple(initvars) - # then funcs + # check if all variables are defined for varname in self.variables: matches = [obs for obs in self.definition.keys() if obs == varname] if len(matches) != 1: print("ERROR: {} variables was not defined.".format(varname)) sys.exit() + # check if variables are filled with values - HOW? + + self.get_weight_str = "JetFlavourUtils::get_weights(rdfslot_, " for var in self.variables: self.get_weight_str += "{},".format(var) diff --git a/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h b/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h index 3d039f56a6d..5feb4da9c6a 100644 --- a/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h +++ b/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h @@ -3,6 +3,7 @@ #include "ROOT/RVec.hxx" #include "edm4hep/ReconstructedParticle.h" +#include "edm4hep/VertexCollection.h" #include "edm4hep/MCParticle.h" #include "edm4hep/Quantity.h" #if __has_include("edm4hep/TrackerHit3DData.h") @@ -57,6 +58,12 @@ namespace FCCAnalyses { rv::RVec get_type(const rv::RVec&); rv::RVec get_charge(const rv::RVec&); + // retrieve collections from full sim + // ROOT::VecOps::RVec get_trackstate(const rv::RVec &particles); + + // Primary vertex + TLorentzVector get_primary_vertex(ROOT::VecOps::RVec& prim_vertex); + //displacement rv::RVec get_d0(const rv::RVec&, const ROOT::VecOps::RVec&); @@ -76,75 +83,112 @@ namespace FCCAnalyses { rv::RVec XPtoPar_dxy(const rv::RVec&, const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, const TLorentzVector& V, // primary vertex const float&); rv::RVec XPtoPar_dz(const rv::RVec&, - const ROOT::VecOps::RVec&, - const TLorentzVector& V, // primary vertex - const float&); + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + const TLorentzVector& V, // primary vertex + const float&); rv::RVec XPtoPar_phi(const rv::RVec&, - const ROOT::VecOps::RVec&, - const TLorentzVector& V, // primary vertex - const float&); + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + const TLorentzVector& V, // primary vertex + const float&); rv::RVec XPtoPar_C(const rv::RVec&, - const ROOT::VecOps::RVec&, - const float&); + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + const float&); rv::RVec XPtoPar_ct(const rv::RVec&, - const ROOT::VecOps::RVec&, - const float&); + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + const float&); + //covariance matrix - //diagonal + //diagonal - 0: d0d0, 2: phiphi, 5: omegaomega, 9: z0z0, 14: tanLambdatanLambda rv::RVec get_omega_cov(const rv::RVec&, - const ROOT::VecOps::RVec&); + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 5); rv::RVec get_d0_cov(const rv::RVec&, - const ROOT::VecOps::RVec& ); - - rv::RVec get_z0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_phi0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_tanlambda_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - //off-diag - rv::RVec get_d0_z0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_phi0_d0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_phi0_z0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_tanlambda_phi0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_tanlambda_d0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_tanlambda_z0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_omega_tanlambda_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_omega_phi0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_omega_d0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); - - rv::RVec get_omega_z0_cov(const rv::RVec& jcs, - const ROOT::VecOps::RVec& tracks); + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 0); + + rv::RVec get_z0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 9); + + rv::RVec get_phi0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 2); + + rv::RVec get_tanlambda_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 14); + //off-diag - 1: phid0, 3: d0omega, 4: phiomega, 6: d0z0, 7: phiz0, 8: omegaz0, 10: d0tanLambda, 11: phitanLambda, 12: omegatanLambda, 13: tanLambdaz0 + rv::RVec get_d0_z0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 6); + + rv::RVec get_phi0_d0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 1); + + rv::RVec get_phi0_z0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 7); + + rv::RVec get_tanlambda_phi0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 11); + + rv::RVec get_tanlambda_d0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 10); + + rv::RVec get_tanlambda_z0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 13); + + rv::RVec get_omega_tanlambda_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 12); + + rv::RVec get_omega_phi0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 4); + + rv::RVec get_omega_d0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 3); + + rv::RVec get_omega_z0_cov(const rv::RVec&, + const ROOT::VecOps::RVec&, + const ROOT::VecOps::RVec&, + int cov_index = 8); rv::RVec get_dndx(const rv::RVec& jcs, const rv::RVec& dNdx, const rv::RVec& trackdata, const rv::RVec JetsConstituents_isChargedHad); + rv::RVec get_dndx_dummy(const rv::RVec &jcs); // for full sim rv::RVec get_Sip2dVal(const rv::RVec& jets, const rv::RVec& jcs, @@ -162,7 +206,8 @@ namespace FCCAnalyses { rv::RVec get_Sip2dSig(const rv::RVec& Sip2dVals, - const rv::RVec& err2_D0); + const rv::RVec& err2_D0, + const std::string& sim_type); rv::RVec get_Sip3dVal(const rv::RVec& jets, const rv::RVec& jcs, @@ -181,7 +226,8 @@ namespace FCCAnalyses { rv::RVec get_Sip3dSig(const rv::RVec& Sip3dVals, const rv::RVec& err2_D0, - const rv::RVec& err2_Z0); + const rv::RVec& err2_Z0, + const std::string& sim_type); rv::RVec get_JetDistVal(const rv::RVec& jets, const rv::RVec& jcs, @@ -200,7 +246,8 @@ namespace FCCAnalyses { rv::RVec get_JetDistSig(const rv::RVec& JetDistVal, const rv::RVec& err2_D0, - const rv::RVec& err2_Z0); + const rv::RVec& err2_Z0, + const std::string& sim_type); rv::RVec get_mtof(const rv::RVec& jcs, const rv::RVec& track_L, @@ -211,6 +258,7 @@ namespace FCCAnalyses { const rv::RVec& calohits, const TLorentzVector& V // primary vertex ); + rv::RVec get_mtof_dummy(const rv::RVec &jcs); // for full sim rv::RVec get_PIDs(const ROOT::VecOps::RVec< int > recin, diff --git a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h index d090e77ec3f..104017280f9 100644 --- a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h +++ b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h @@ -45,27 +45,32 @@ namespace ReconstructedParticle2Track{ const ROOT::VecOps::RVec& tracks); //here only computed for the first charged particle encountered ROOT::VecOps::RVec XPtoPar_dxy(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, const TLorentzVector& V, // primary vertex const float& Bz); ROOT::VecOps::RVec XPtoPar_dz(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const TLorentzVector& V, // primary vertex - const float& Bz); + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, + const TLorentzVector& V, // primary vertex + const float& Bz); ROOT::VecOps::RVec XPtoPar_phi(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const TLorentzVector& V, // primary vertex - const float& Bz); + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, + const TLorentzVector& V, // primary vertex + const float& Bz); ROOT::VecOps::RVec XPtoPar_C(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const float& Bz); + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, + const float& Bz); ROOT::VecOps::RVec XPtoPar_ct(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const float& Bz); + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, + const float& Bz); /// Return the D0 of a track to a reconstructed particle ROOT::VecOps::RVec getRP2TRK_D0 (ROOT::VecOps::RVec in, @@ -95,66 +100,16 @@ namespace ReconstructedParticle2Track{ ROOT::VecOps::RVec getRP2TRK_Z0_sig (ROOT::VecOps::RVec in, ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the the D0 of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_D0_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the the Z0 of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_Z0_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the the Phi of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_phi_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the omega of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_omega_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the tanLambda of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_tanLambda_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0, phi0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_phi0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0, omega) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_omega_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0,z0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (phi0,omega) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_phi0_omega_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (phi0,z0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_phi0_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (phi0,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_phi0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) ; - - /// Return the off-diag term (omega,z0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_omega_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) ; - - /// Return the off-diag term (omega,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_omega_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) ; - - /// Return the off-diag term (z0,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_z0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); + /* + Return the covariance matrix of a track to a reconstructed particle + @param cov_index: the index of the covariance matrix element: + - Diagonal elements are: 0: d0d0, 2: phiphi, 5: omegaomega, 9: z0z0, 14: tanLambdatanLambda + - Off-diagonal elements are: 1: phid0, 3: d0omega, 4: phiomega, 6: d0z0, 7: phiz0, 8: omegaz0, 10: d0tanLambda, 11: phitanLambda, 12: omegatanLambda, 13: tanLambdaz0 + */ + ROOT::VecOps::RVec get_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks, + ROOT::VecOps::RVec trackstates, + int cov_index); /// Return the tracks associated to reco'ed particles diff --git a/analyzers/dataframe/src/JetConstituentsUtils.cc b/analyzers/dataframe/src/JetConstituentsUtils.cc index 37537e3ff49..c60632cb9e6 100644 --- a/analyzers/dataframe/src/JetConstituentsUtils.cc +++ b/analyzers/dataframe/src/JetConstituentsUtils.cc @@ -120,6 +120,16 @@ namespace FCCAnalyses return out; }; + auto cast_constituent_5 = [](const auto &jcs, const auto &coll1, const auto &coll2, const auto &coll3, const auto &coll4, auto &&meth) + { + rv::RVec out; + for (const auto &jc : jcs) + { + out.emplace_back(meth(jc, coll1, coll2, coll3, coll4)); + } + return out; + }; + rv::RVec get_Bz(const rv::RVec &jcs, const ROOT::VecOps::RVec &tracks) { @@ -194,6 +204,44 @@ namespace FCCAnalyses return out; } + // // get trackstate collection from full sim + // ROOT::VecOps::RVec get_trackstate(const rv::RVec &particles) + // { + // ROOT::VecOps::RVec tracks; + // for (const auto &particle : particles) + // { + // if (particle.getTracks().size() > 0) + // { + // tracks.push_back(particle.getTracks()[0].getTrackStates()); + // } + // } + // return tracks; + // } + + // Primary vertex + + TLorentzVector get_primary_vertex(ROOT::VecOps::RVec& prim_vertex) + { + TLorentzVector pv_pos(0, 0, 0, 0); // Initialize primary vertex position + int i = 0; + for (const auto& pv : prim_vertex) + { + if (i > 0) // only one primary vertex is expected + { + throw std::invalid_argument("More than one primary vertex found in the event."); + } else + { + pv_pos.SetXYZT(pv.position.x, pv.position.y, pv.position.z, 0.0); // position of PV in mm + i++; + } + } + if (i == 0) + { + std::cout << "No primary vertex found in the event. Using (0,0,0, 0)" << std::endl; + } + return pv_pos; + } + // displacement (wrt (0,0,0)) rv::RVec get_d0(const rv::RVec &jcs, const ROOT::VecOps::RVec &tracks) @@ -226,138 +274,173 @@ namespace FCCAnalyses } rv::RVec XPtoPar_dxy(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + const ROOT::VecOps::RVec& tracks, const TLorentzVector &V, // primary vertex posotion and time in mm const float &Bz) { - return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_dxy); + return cast_constituent_5(jcs, trackstates, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_dxy); } rv::RVec XPtoPar_dz(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + const ROOT::VecOps::RVec& tracks, const TLorentzVector &V, // primary vertex posotion and time in mm const float &Bz) { - return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_dz); + return cast_constituent_5(jcs, trackstates, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_dz); } rv::RVec XPtoPar_phi(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + const ROOT::VecOps::RVec& tracks, const TLorentzVector &V, // primary vertex posotion and time in mm const float &Bz) { - return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_phi); + return cast_constituent_5(jcs, trackstates, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_phi); } rv::RVec XPtoPar_C(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + const ROOT::VecOps::RVec& tracks, const float &Bz) { - return cast_constituent_3(jcs, tracks, Bz, ReconstructedParticle2Track::XPtoPar_C); + return cast_constituent_4(jcs, trackstates, tracks, Bz, ReconstructedParticle2Track::XPtoPar_C); } rv::RVec XPtoPar_ct(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + const ROOT::VecOps::RVec& tracks, const float &Bz) { - return cast_constituent_3(jcs, tracks, Bz, ReconstructedParticle2Track::XPtoPar_ct); + return cast_constituent_4(jcs, trackstates, tracks, Bz, ReconstructedParticle2Track::XPtoPar_ct); } // Covariance matrix elements of tracks parameters // diagonal rv::RVec get_omega_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_omega_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_d0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_z0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Z0_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_phi0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_tanlambda_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_tanLambda_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } // off-diagonal rv::RVec get_d0_z0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_z0_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_phi0_d0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_phi0_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_phi0_z0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi0_z0_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_tanlambda_phi0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi0_tanlambda_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_tanlambda_d0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_tanlambda_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_tanlambda_z0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_z0_tanlambda_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_omega_tanlambda_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_omega_tanlambda_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_omega_phi0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi0_omega_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_omega_d0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_omega_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } rv::RVec get_omega_z0_cov(const rv::RVec &jcs, - const ROOT::VecOps::RVec &tracks) + const ROOT::VecOps::RVec &tracks, + const ROOT::VecOps::RVec &trackstates, + int cov_index) { - return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_omega_z0_cov); + return cast_constituent_4(jcs, tracks, trackstates, cov_index, ReconstructedParticle2Track::get_cov); } // neutrals are set to 0; muons and electrons are also set to 0; @@ -393,6 +476,38 @@ namespace FCCAnalyses return out; } + // dummy dndx function for full sim + rv::RVec get_dndx_dummy(const rv::RVec &jcs) { + rv::RVec out; + for (int i = 0; i < jcs.size(); ++i) + { + FCCAnalysesJetConstituents ct = jcs.at(i); + FCCAnalysesJetConstituentsData tmp; + for (int j = 0; j < ct.size(); ++j) + { + tmp.push_back(0.); + } + out.push_back(tmp); + } + return out; + } + + // dummy mtof function for full sim + rv::RVec get_mtof_dummy(const rv::RVec &jcs) { + rv::RVec out; + for (int i = 0; i < jcs.size(); ++i) + { + FCCAnalysesJetConstituents ct = jcs.at(i); + FCCAnalysesJetConstituentsData tmp; + for (int j = 0; j < ct.size(); ++j) + { + tmp.push_back(0.); + } + out.push_back(tmp); + } + return out; + } + rv::RVec get_Sip2dVal(const rv::RVec &jets, const rv::RVec &jcs, const ROOT::VecOps::RVec &tracks) @@ -482,7 +597,8 @@ namespace FCCAnalyses /// The functions get_Sip2dSig and get_Sip2dVal can be made independent; /// I passed to the former the result of the latter, avoiding the recomputation rv::RVec get_Sip2dSig(const rv::RVec &Sip2dVals, - const rv::RVec &err2_D0) + const rv::RVec &err2_D0, + const std::string& sim_type) { rv::RVec out; for (int i = 0; i < Sip2dVals.size(); ++i) @@ -496,7 +612,7 @@ namespace FCCAnalyses } else { - s.push_back(-9); + s.push_back(sim_type == "fast" ? -9 : -200); // dummy value in full sim is -200 because -9 in fast sim is still inside the distribution } } out.push_back(s); @@ -595,7 +711,8 @@ namespace FCCAnalyses rv::RVec get_Sip3dSig(const rv::RVec &Sip3dVals, const rv::RVec &err2_D0, - const rv::RVec &err2_Z0) + const rv::RVec &err2_Z0, + const std::string& sim_type) { rv::RVec out; for (int i = 0; i < Sip3dVals.size(); ++i) @@ -609,7 +726,7 @@ namespace FCCAnalyses } else { - s.push_back(-9); + s.push_back(sim_type == "fast" ? -9 : -200); // dummy value in full sim is -200 because -9 in fast sim is still inside the distribution } } out.push_back(s); @@ -719,7 +836,8 @@ namespace FCCAnalyses rv::RVec get_JetDistSig(const rv::RVec &JetDistVal, const rv::RVec &err2_D0, - const rv::RVec &err2_Z0) + const rv::RVec &err2_Z0, + const std::string& sim_type) { rv::RVec out; for (int i = 0; i < JetDistVal.size(); ++i) @@ -735,7 +853,7 @@ namespace FCCAnalyses } else { - tmp.push_back(-9.); + tmp.push_back(sim_type == "fast" ? -9 : -200); // dummy value in full sim is -200 because -9 in fast sim is still inside the distribution } } out.push_back(tmp); @@ -1089,21 +1207,22 @@ namespace FCCAnalyses rv::RVec out; for (int i = 0; i < jcs.size(); ++i) { - FCCAnalysesJetConstituentsData is_El; + FCCAnalysesJetConstituentsData is_Muon; FCCAnalysesJetConstituents ct = jcs.at(i); for (int j = 0; j < ct.size(); ++j) { - if (abs(ct.at(j).charge) > 0 and abs(ct.at(j).mass - 0.000510999) < 1.e-05) +#if edm4hep_VERSION > EDM4HEP_VERSION(0, 10, 5) + if (std::abs(ct.at(j).PDG) == 11) +#else + if (std::abs(ct.at(j).type) == 11) +#endif { - is_El.push_back(1.); + is_Muon.push_back(1.); } else - { - is_El.push_back(0.); - } + is_Muon.push_back(0.); } - - out.push_back(is_El); + out.push_back(is_Muon); } return out; } @@ -1113,21 +1232,22 @@ namespace FCCAnalyses rv::RVec out; for (int i = 0; i < jcs.size(); ++i) { - FCCAnalysesJetConstituentsData is_Mu; + FCCAnalysesJetConstituentsData is_Muon; FCCAnalysesJetConstituents ct = jcs.at(i); for (int j = 0; j < ct.size(); ++j) { - if (abs(ct.at(j).charge) > 0 and abs(ct.at(j).mass - 0.105658) < 1.e-03) +#if edm4hep_VERSION > EDM4HEP_VERSION(0, 10, 5) + if (std::abs(ct.at(j).PDG) == 13) +#else + if (std::abs(ct.at(j).type) == 13) +#endif { - is_Mu.push_back(1.); + is_Muon.push_back(1.); } else - { - is_Mu.push_back(0.); - } + is_Muon.push_back(0.); } - - out.push_back(is_Mu); + out.push_back(is_Muon); } return out; } @@ -1141,7 +1261,17 @@ namespace FCCAnalyses FCCAnalysesJetConstituents ct = jcs.at(i); for (int j = 0; j < ct.size(); ++j) { - if (abs(ct.at(j).charge) > 0 and abs(ct.at(j).mass - 0.13957) < 1.e-03) + int num_tracks = ct.at(j).tracks_end - ct.at(j).tracks_begin; + // check if num_tracks is valid (e.g. 0 or 1) + if (num_tracks < 0 || num_tracks > 1) + { + throw std::invalid_argument("Invalid number of tracks for constituent. Must be 0 or 1."); + } +#if edm4hep_VERSION > EDM4HEP_VERSION(0, 10, 5) + if (std::abs(ct.at(j).PDG) != 11 && std::abs(ct.at(j).PDG) != 13 && std::abs(ct.at(j).PDG) != 22 && num_tracks == 1) +#else + if (std::abs(ct.at(j).type) != 11 && std::abs(ct.at(j).type) != 13 && std::abs(ct.at(j).type) != 22 && num_tracks == 1) +#endif { is_ChargedHad.push_back(1.); } @@ -1165,10 +1295,16 @@ namespace FCCAnalyses FCCAnalysesJetConstituents ct = jcs.at(i); for (int j = 0; j < ct.size(); ++j) { + int num_tracks = ct.at(j).tracks_end - ct.at(j).tracks_begin; + // check if num_tracks is valid (e.g. 0 or 1) + if (num_tracks < 0 || num_tracks > 1) + { + throw std::invalid_argument("Invalid number of tracks for constituent. Must be 0 or 1."); + } #if edm4hep_VERSION > EDM4HEP_VERSION(0, 10, 5) - if (ct.at(j).PDG == 130) + if (std::abs(ct.at(j).PDG) != 11 && std::abs(ct.at(j).PDG) != 13 && std::abs(ct.at(j).PDG) != 22 && num_tracks == 0) #else - if (ct.at(j).type == 130) + if (std::abs(ct.at(j).type) != 11 && std::abs(ct.at(j).type) != 13 && std::abs(ct.at(j).type) != 22 && num_tracks == 0) #endif { is_NeutralHad.push_back(1.); diff --git a/analyzers/dataframe/src/ReconstructedParticle2Track.cc b/analyzers/dataframe/src/ReconstructedParticle2Track.cc index f8aa59e8901..82eb3331d09 100644 --- a/analyzers/dataframe/src/ReconstructedParticle2Track.cc +++ b/analyzers/dataframe/src/ReconstructedParticle2Track.cc @@ -61,8 +61,10 @@ namespace ReconstructedParticle2Track{ return Bz; } + ROOT::VecOps::RVec XPtoPar_dxy(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, const TLorentzVector& V, // primary vertex const float& Bz) { @@ -71,12 +73,12 @@ namespace ReconstructedParticle2Track{ ROOT::VecOps::RVec out; for (const auto & rp: in) { + if(rp.tracks_begin - rp.tracks_end >0) { // if any tracks + auto track = tracks.at(rp.tracks_begin); - if( rp.tracks_begin < tracks.size()) { - - float D0_wrt0 = tracks.at(rp.tracks_begin).D0; - float Z0_wrt0 = tracks.at(rp.tracks_begin).Z0; - float phi0_wrt0 = tracks.at(rp.tracks_begin).phi; + float D0_wrt0 = trackstates.at(track.trackStates_begin).D0; + float Z0_wrt0 = trackstates.at(track.trackStates_begin).Z0; + float phi0_wrt0 = trackstates.at(track.trackStates_begin).phi; TVector3 X( - D0_wrt0 * TMath::Sin(phi0_wrt0) , D0_wrt0 * TMath::Cos(phi0_wrt0) , Z0_wrt0); TVector3 x = X - V.Vect(); @@ -97,16 +99,16 @@ namespace ReconstructedParticle2Track{ out.push_back(D); } else { - out.push_back(-9.); + out.push_back(-9.); } } return out; } - ROOT::VecOps::RVec XPtoPar_dz(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, const TLorentzVector& V, // primary vertex const float& Bz) { @@ -116,11 +118,12 @@ namespace ReconstructedParticle2Track{ for (const auto & rp: in) { - if( rp.tracks_begin < tracks.size()) { + if(rp.tracks_begin - rp.tracks_end >0) { // if any tracks + auto track = tracks.at(rp.tracks_begin); - float D0_wrt0 = tracks.at(rp.tracks_begin).D0; - float Z0_wrt0 = tracks.at(rp.tracks_begin).Z0; - float phi0_wrt0 = tracks.at(rp.tracks_begin).phi; + float D0_wrt0 = trackstates.at(track.trackStates_begin).D0; + float Z0_wrt0 = trackstates.at(track.trackStates_begin).Z0; + float phi0_wrt0 = trackstates.at(track.trackStates_begin).phi; TVector3 X( - D0_wrt0 * TMath::Sin(phi0_wrt0) , D0_wrt0 * TMath::Cos(phi0_wrt0) , Z0_wrt0); TVector3 x = X - V.Vect(); @@ -154,21 +157,22 @@ namespace ReconstructedParticle2Track{ } ROOT::VecOps::RVec XPtoPar_phi(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const TLorentzVector& V, // primary vertex - const float& Bz) { + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, + const TLorentzVector& V, // primary vertex + const float& Bz) { const double cSpeed = 2.99792458e8 * 1.0e-9; //Reduced speed of light ??? ROOT::VecOps::RVec out; for (const auto & rp: in) { + if(rp.tracks_begin - rp.tracks_end >0) { // if any tracks + auto track = tracks.at(rp.tracks_begin); - if( rp.tracks_begin < tracks.size()) { - - float D0_wrt0 = tracks.at(rp.tracks_begin).D0; - float Z0_wrt0 = tracks.at(rp.tracks_begin).Z0; - float phi0_wrt0 = tracks.at(rp.tracks_begin).phi; + float D0_wrt0 = trackstates.at(track.trackStates_begin).D0; + float Z0_wrt0 = trackstates.at(track.trackStates_begin).Z0; + float phi0_wrt0 = trackstates.at(track.trackStates_begin).phi; TVector3 X( - D0_wrt0 * TMath::Sin(phi0_wrt0) , D0_wrt0 * TMath::Cos(phi0_wrt0) , Z0_wrt0); TVector3 x = X - V.Vect(); @@ -182,7 +186,7 @@ namespace ReconstructedParticle2Track{ double T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2); double phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); - out.push_back(phi0); + out.push_back(phi0); } else { out.push_back(-9.); @@ -192,23 +196,24 @@ namespace ReconstructedParticle2Track{ } ROOT::VecOps::RVec XPtoPar_C(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const float& Bz) { + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, + const float& Bz) { const double cSpeed = 2.99792458e8 * 1.0e3 * 1.0e-15; ROOT::VecOps::RVec out; for (const auto & rp: in) { - - if( rp.tracks_begin < tracks.size()) { + if(rp.tracks_begin - rp.tracks_end >0) { // if any tracks + auto track = tracks.at(rp.tracks_begin); TVector3 p(rp.momentum.x, rp.momentum.y, rp.momentum.z); double a = std::copysign(1.0, rp.charge) * Bz * cSpeed; - double pt = p.Pt(); + double pt = p.Pt(); double C = a/(2 * pt); - out.push_back(C); + out.push_back(C); } else { out.push_back(-9.); } @@ -217,22 +222,21 @@ namespace ReconstructedParticle2Track{ } ROOT::VecOps::RVec XPtoPar_ct(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const float& Bz) { + const ROOT::VecOps::RVec& trackstates, + const ROOT::VecOps::RVec& tracks, + const float& Bz) { const double cSpeed = 2.99792458e8 * 1.0e-9; ROOT::VecOps::RVec out; for (const auto & rp: in) { - - if( rp.tracks_begin < tracks.size()) { + if(rp.tracks_begin - rp.tracks_end >0) { // if any tracks + auto track = tracks.at(rp.tracks_begin); TVector3 p(rp.momentum.x, rp.momentum.y, rp.momentum.z); - double pt = p.Pt(); - + double pt = p.Pt(); double ct = p(2) / pt; - - out.push_back(ct); + out.push_back(ct); } else { out.push_back(-9.); @@ -254,14 +258,23 @@ getRP2TRK_D0(ROOT::VecOps::RVec in, return result; } + +// cov_index accesses the different elements of the covariance matrix. +// diagonal elements are: 0: d0d0, 2: phiphi, 5: omegaomega, 9: z0z0, 14: tanLambdatanLambda +// off-diagonal elements are: 1: phid0, 3: d0omega, 4: phiomega, 6: d0z0, 7: phiz0, 8: omegaz0, 10: d0tanLambda, 11: phitanLambda, 12: omegatanLambda, 13: tanLambdaz0 ROOT::VecOps::RVec -getRP2TRK_D0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { +get_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks, + ROOT::VecOps::RVec trackstates, + int cov_index) { ROOT::VecOps::RVec result; for (auto & p: in) { - if (p.tracks_begin 0) { // if any tracks + auto track = tracks.at(p.tracks_begin); + result.push_back(trackstates.at(track.trackStates_begin).covMatrix[cov_index]); + } else { + result.push_back(-9.); + } } return result; } @@ -290,18 +303,6 @@ getRP2TRK_Z0(ROOT::VecOps::RVec in, return result; } -ROOT::VecOps::RVec -getRP2TRK_Z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin getRP2TRK_Z0_sig(ROOT::VecOps::RVec in, ROOT::VecOps::RVec tracks) { @@ -326,18 +327,6 @@ getRP2TRK_phi(ROOT::VecOps::RVec in, return result; } -ROOT::VecOps::RVec -getRP2TRK_phi_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin getRP2TRK_omega(ROOT::VecOps::RVec in, @@ -351,17 +340,6 @@ getRP2TRK_omega(ROOT::VecOps::RVec in, return result; } -ROOT::VecOps::RVec -getRP2TRK_omega_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin getRP2TRK_tanLambda(ROOT::VecOps::RVec in, @@ -375,137 +353,6 @@ getRP2TRK_tanLambda(ROOT::VecOps::RVec in, return result; } -ROOT::VecOps::RVec -getRP2TRK_tanLambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_d0_phi0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_d0_omega_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_d0_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_d0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_phi0_omega_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_phi0_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_phi0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_omega_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_omega_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin -getRP2TRK_z0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) { - ROOT::VecOps::RVec result; - for (auto & p: in) { - if (p.tracks_begin getRP2TRK( ROOT::VecOps::RVec in,