diff --git a/RecoTracker/LSTCore/interface/SegmentsSoA.h b/RecoTracker/LSTCore/interface/SegmentsSoA.h index 891603f8f0ced..e274e8923a1db 100644 --- a/RecoTracker/LSTCore/interface/SegmentsSoA.h +++ b/RecoTracker/LSTCore/interface/SegmentsSoA.h @@ -15,6 +15,15 @@ namespace lst { SOA_COLUMN(FPX, dPhiChanges), SOA_COLUMN(FPX, dPhiChangeMins), SOA_COLUMN(FPX, dPhiChangeMaxs), +#ifdef CUT_VALUE_DEBUG + SOA_COLUMN(FPX, zHis), + SOA_COLUMN(FPX, zLos), + SOA_COLUMN(FPX, rtHis), + SOA_COLUMN(FPX, rtLos), + SOA_COLUMN(FPX, dAlphaInners), + SOA_COLUMN(FPX, dAlphaOuters), + SOA_COLUMN(FPX, dAlphaInnerOuters), +#endif SOA_COLUMN(uint16_t, innerLowerModuleIndices), SOA_COLUMN(uint16_t, outerLowerModuleIndices), SOA_COLUMN(Params_LS::ArrayUxLayers, mdIndices), diff --git a/RecoTracker/LSTCore/src/alpaka/Segment.h b/RecoTracker/LSTCore/src/alpaka/Segment.h index 65625ce52ce2b..44867aed92dd4 100644 --- a/RecoTracker/LSTCore/src/alpaka/Segment.h +++ b/RecoTracker/LSTCore/src/alpaka/Segment.h @@ -192,6 +192,15 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float dPhiChange, float dPhiChangeMin, float dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + float zHi, + float zLo, + float rtHi, + float rtLo, + float dAlphaInner, + float dAlphaOuter, + float dAlphaInnerOuter, +#endif unsigned int idx) { segments.mdIndices()[idx][0] = lowerMDIndex; segments.mdIndices()[idx][1] = upperMDIndex; @@ -206,6 +215,16 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { segments.dPhiChanges()[idx] = __F2H(dPhiChange); segments.dPhiChangeMins()[idx] = __F2H(dPhiChangeMin); segments.dPhiChangeMaxs()[idx] = __F2H(dPhiChangeMax); + +#ifdef CUT_VALUE_DEBUG + segments.zHis()[idx] = __F2H(zHi); + segments.zLos()[idx] = __F2H(zLo); + segments.rtHis()[idx] = __F2H(rtHi); + segments.rtLos()[idx] = __F2H(rtLo); + segments.dAlphaInners()[idx] = __F2H(dAlphaInner); + segments.dAlphaOuters()[idx] = __F2H(dAlphaOuter); + segments.dAlphaInnerOuters()[idx] = __F2H(dAlphaInnerOuter); +#endif } template @@ -356,7 +375,18 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float& dPhiChange, float& dPhiChangeMin, float& dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + float& dAlphaInnerMDSegment, + float& dAlphaOuterMDSegment, + float& dAlphaInnerMDOuterMD, + float& zLo, + float& zHi, +#endif const float ptCut) { +#ifndef CUT_VALUE_DEBUG + float dAlphaInnerMDSegment, dAlphaOuterMDSegment, dAlphaInnerMDOuterMD; + float zLo, zHi; +#endif float xIn, yIn, zIn, rtIn, xOut, yOut, zOut, rtOut; xIn = mds.anchorX()[innerMDIndex]; @@ -374,9 +404,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { const float zGeom = modules.layers()[innerLowerModuleIndex] <= 2 ? 2.f * kPixelPSZpitch : 2.f * kStrip2SZpitch; - float zLo = zIn + (zIn - kDeltaZLum) * (rtOut / rtIn - 1.f) * (zIn > 0.f ? 1.f : dzDrtScale) - - zGeom; //slope-correction only on outer end - float zHi = zIn + (zIn + kDeltaZLum) * (rtOut / rtIn - 1.f) * (zIn < 0.f ? 1.f : dzDrtScale) + zGeom; + //slope-correction only on outer end + zLo = zIn + (zIn - kDeltaZLum) * (rtOut / rtIn - 1.f) * (zIn > 0.f ? 1.f : dzDrtScale) - zGeom; + zHi = zIn + (zIn + kDeltaZLum) * (rtOut / rtIn - 1.f) * (zIn < 0.f ? 1.f : dzDrtScale) + zGeom; if ((zOut < zLo) || (zOut > zHi)) return false; @@ -420,9 +450,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float innerMDAlpha = mds.dphichanges()[innerMDIndex]; float outerMDAlpha = mds.dphichanges()[outerMDIndex]; - float dAlphaInnerMDSegment = innerMDAlpha - dPhiChange; - float dAlphaOuterMDSegment = outerMDAlpha - dPhiChange; - float dAlphaInnerMDOuterMD = innerMDAlpha - outerMDAlpha; + dAlphaInnerMDSegment = innerMDAlpha - dPhiChange; + dAlphaOuterMDSegment = outerMDAlpha - dPhiChange; + dAlphaInnerMDOuterMD = innerMDAlpha - outerMDAlpha; float dAlphaInnerMDSegmentThreshold = dAlphaThresholdValues[0]; float dAlphaOuterMDSegmentThreshold = dAlphaThresholdValues[1]; @@ -449,7 +479,18 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float& dPhiChange, float& dPhiChangeMin, float& dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + float& dAlphaInnerMDSegment, + float& dAlphaOuterMDSegment, + float& dAlphaInnerMDOuterMD, + float& rtLo, + float& rtHi, +#endif const float ptCut) { +#ifndef CUT_VALUE_DEBUG + float dAlphaInnerMDSegment, dAlphaOuterMDSegment, dAlphaInnerMDOuterMD; + float rtLo, rtHi; +#endif float xIn, yIn, zIn, rtIn, xOut, yOut, zOut, rtOut; xIn = mds.anchorX()[innerMDIndex]; @@ -481,10 +522,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float dLum = alpaka::math::copysign(acc, kDeltaZLum, zIn); float drtDzScale = sdSlope / alpaka::math::tan(acc, sdSlope); - float rtLo = alpaka::math::max( - acc, rtIn * (1.f + dz / (zIn + dLum) * drtDzScale) - rtGeom, rtIn - 0.5f * rtGeom); //rt should increase - float rtHi = rtIn * (zOut - dLum) / (zIn - dLum) + - rtGeom; //dLum for luminous; rGeom for measurement size; no tanTheta_loc(pt) correction + //rt should increase + rtLo = alpaka::math::max(acc, rtIn * (1.f + dz / (zIn + dLum) * drtDzScale) - rtGeom, rtIn - 0.5f * rtGeom); + //dLum for luminous; rGeom for measurement size; no tanTheta_loc(pt) correction + rtHi = rtIn * (zOut - dLum) / (zIn - dLum) + rtGeom; // Completeness if ((rtOut < rtLo) || (rtOut > rtHi)) @@ -544,9 +585,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float innerMDAlpha = mds.dphichanges()[innerMDIndex]; float outerMDAlpha = mds.dphichanges()[outerMDIndex]; - float dAlphaInnerMDSegment = innerMDAlpha - dPhiChange; - float dAlphaOuterMDSegment = outerMDAlpha - dPhiChange; - float dAlphaInnerMDOuterMD = innerMDAlpha - outerMDAlpha; + dAlphaInnerMDSegment = innerMDAlpha - dPhiChange; + dAlphaOuterMDSegment = outerMDAlpha - dPhiChange; + dAlphaInnerMDOuterMD = innerMDAlpha - outerMDAlpha; float dAlphaInnerMDSegmentThreshold = dAlphaThresholdValues[0]; float dAlphaOuterMDSegmentThreshold = dAlphaThresholdValues[1]; @@ -573,8 +614,21 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float& dPhiChange, float& dPhiChangeMin, float& dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + float& dAlphaInnerMDSegment, + float& dAlphaOuterMDSegment, + float& dAlphaInnerMDOuterMD, + float& zLo, + float& zHi, + float& rtLo, + float& rtHi, +#endif const float ptCut) { if (modules.subdets()[innerLowerModuleIndex] == Barrel and modules.subdets()[outerLowerModuleIndex] == Barrel) { +#ifdef CUT_VALUE_DEBUG + rtLo = -999.f; + rtHi = -999.f; +#endif return runSegmentDefaultAlgoBarrel(acc, modules, mds, @@ -588,8 +642,19 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { dPhiChange, dPhiChangeMin, dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + dAlphaInnerMDSegment, + dAlphaOuterMDSegment, + dAlphaInnerMDOuterMD, + zLo, + zHi, +#endif ptCut); } else { +#ifdef CUT_VALUE_DEBUG + zLo = -999.f; + zHi = -999.f; +#endif return runSegmentDefaultAlgoEndcap(acc, modules, mds, @@ -603,6 +668,13 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { dPhiChange, dPhiChangeMin, dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + dAlphaInnerMDSegment, + dAlphaOuterMDSegment, + dAlphaInnerMDOuterMD, + rtLo, + rtHi, +#endif ptCut); } } @@ -643,6 +715,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float dPhi, dPhiMin, dPhiMax, dPhiChange, dPhiChangeMin, dPhiChangeMax; +#ifdef CUT_VALUE_DEBUG + float zLo, zHi, rtLo, rtHi, dAlphaInnerMDSegment, dAlphaOuterMDSegment, dAlphaInnerMDOuterMD; +#endif + unsigned int innerMiniDoubletAnchorHitIndex = mds.anchorHitIndices()[innerMDIndex]; unsigned int outerMiniDoubletAnchorHitIndex = mds.anchorHitIndices()[outerMDIndex]; dPhiMin = 0; @@ -662,6 +738,15 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { dPhiChange, dPhiChangeMin, dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + dAlphaInnerMDSegment, + dAlphaOuterMDSegment, + dAlphaInnerMDOuterMD, + zLo, + zHi, + rtLo, + rtHi, +#endif ptCut)) { unsigned int totOccupancySegments = alpaka::atomicAdd(acc, @@ -692,6 +777,15 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { dPhiChange, dPhiChangeMin, dPhiChangeMax, +#ifdef CUT_VALUE_DEBUG + zHi, + zLo, + rtHi, + rtLo, + dAlphaInnerMDSegment, + dAlphaOuterMDSegment, + dAlphaInnerMDOuterMD, +#endif segmentIdx); } } diff --git a/RecoTracker/LSTCore/standalone/analysis/DNN/embed_train.ipynb b/RecoTracker/LSTCore/standalone/analysis/DNN/embed_train.ipynb index 1aa64d3ea3c84..e3c35887f046b 100644 --- a/RecoTracker/LSTCore/standalone/analysis/DNN/embed_train.ipynb +++ b/RecoTracker/LSTCore/standalone/analysis/DNN/embed_train.ipynb @@ -67,18 +67,18 @@ " 't5_pMatched',\n", " 't5_sim_vxy',\n", " 't5_sim_vz',\n", - " 't5_matched_simIdx'\n", + " 't5_simIdxAll'\n", "]\n", "\n", "branches_list += [\n", " 'pLS_eta',\n", " 'pLS_etaErr',\n", " 'pLS_phi',\n", - " 'pLS_matched_simIdx',\n", + " 'pLS_simIdxAll',\n", " 'pLS_circleCenterX',\n", " 'pLS_circleCenterY',\n", " 'pLS_circleRadius',\n", - " 'pLS_ptIn',\n", + " 'pLS_pt',\n", " 'pLS_ptErr',\n", " 'pLS_px',\n", " 'pLS_py',\n", @@ -236,7 +236,7 @@ " disp_evt.append(branches['t5_sim_vxy'][ev][i])\n", "\n", " # first (or only) matched sim-index, -1 if none -----------------------\n", - " simIdx_list = branches['t5_matched_simIdx'][ev][i]\n", + " simIdx_list = branches['t5_simIdxAll'][ev][i]\n", " sim_evt.append(simIdx_list[0] if len(simIdx_list) else -1)\n", "\n", " # push to global containers ----------------------------------------------\n", @@ -287,7 +287,7 @@ " circleCenterX = np.abs(branches['pLS_circleCenterX'][ev][i])\n", " circleCenterY = np.abs(branches['pLS_circleCenterY'][ev][i])\n", " circleRadius = branches['pLS_circleRadius'][ev][i]\n", - " ptIn = branches['pLS_ptIn'][ev][i]\n", + " ptIn = branches['pLS_pt'][ev][i]\n", " ptErr = branches['pLS_ptErr'][ev][i]\n", " isQuad = branches['pLS_isQuad'][ev][i]\n", "\n", @@ -308,7 +308,7 @@ " feat_evt.append(f)\n", " eta_evt.append(eta)\n", "\n", - " sim_list = branches['pLS_matched_simIdx'][ev][i]\n", + " sim_list = branches['pLS_simIdxAll'][ev][i]\n", " sim_evt.append(sim_list[0] if len(sim_list) else -1)\n", "\n", " # ――― store per‑event containers -----------------------------------------\n", diff --git a/RecoTracker/LSTCore/standalone/analysis/DNN/train_pT3_DNN.ipynb b/RecoTracker/LSTCore/standalone/analysis/DNN/train_pT3_DNN.ipynb index 8488c9b0a766a..35051dbe5a047 100644 --- a/RecoTracker/LSTCore/standalone/analysis/DNN/train_pT3_DNN.ipynb +++ b/RecoTracker/LSTCore/standalone/analysis/DNN/train_pT3_DNN.ipynb @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -49,17 +49,12 @@ " 'pT3_eta',\n", " 'pT3_phi',\n", " 'pT3_score',\n", - " 'pT3_foundDuplicate',\n", - " 'pT3_matched_simIdx',\n", - " 'pT3_hitIdxs',\n", " 'pT3_pixelRadius',\n", " 'pT3_pixelRadiusError',\n", " 'pT3_tripletRadius',\n", " 'pT3_rPhiChiSquared',\n", " 'pT3_rPhiChiSquaredInwards',\n", " 'pT3_rzChiSquared',\n", - " 'pT3_layer_binary',\n", - " 'pT3_moduleType_binary'\n", "]\n", "\n", "file_path = \"pt3_500_fixed_new_2.root\"\n", diff --git a/RecoTracker/LSTCore/standalone/bin/lst.cc b/RecoTracker/LSTCore/standalone/bin/lst.cc index bf7f64a66643f..3fcecc7aa6262 100644 --- a/RecoTracker/LSTCore/standalone/bin/lst.cc +++ b/RecoTracker/LSTCore/standalone/bin/lst.cc @@ -65,13 +65,20 @@ int main(int argc, char **argv) { "w,write_ntuple", "Write Ntuple", cxxopts::value()->default_value("1"))( "s,streams", "Set number of streams", cxxopts::value()->default_value("1"))( "d,debug", "Run debug job. i.e. overrides output option to 'debug.root' and 'recreate's the file.")( - "l,lower_level", "write lower level objects ntuple results")("G,gnn_ntuple", "write gnn input variable ntuple")( + "l,lower_level", "write lower level objects ntuple results")( "j,nsplit_jobs", "Enable splitting jobs by N blocks (--job_index must be set)", cxxopts::value())( "I,job_index", "job_index of split jobs (--nsplit_jobs must be set. index starts from 0. i.e. 0, 1, 2, 3, etc...)", cxxopts::value())("3,tc_pls_triplets", "Allow triplet pLSs in TC collection")( "2,no_pls_dupclean", "Disable pLS duplicate cleaning (both steps)")("h,help", "Print help")( - "J,jet_branches", "Accounts for specific jet branches in input root file for testing"); + "md", "Write MD branches in output ntuple.")("ls", "Write LS branches in output ntuple.")( + "t3", "Write T3 branches in output ntuple.")("t5", "Write T5 branches in output ntuple.")( + "pls", "Write pLS branches in output ntuple.")("pt3", "Write pT3 branches in output ntuple.")( + "pt5", "Write pT5 branches in output ntuple.")("occ", "Write occupancy branches in output ntuple.")( + "t5dnn", "Write T5 DNN branches in output ntuple.")("t3dnn", "Write T3 DNN branches in output ntuple.")( + "pt3dnn", "Write pT3 DNN branches in output ntuple.")("allobj", "Write all object branches in output ntuple.")( + "J,jet", "Accounts for specific jet branches in input root file for testing")( + "sim", "Write extra sim branches in output ntuple"); auto result = options.parse(argc, argv); @@ -241,20 +248,6 @@ int main(int argc, char **argv) { ana.do_lower_level = false; } - //_______________________________________________________________________________ - // --gnn_ntuple - if (result.count("gnn_ntuple")) { - ana.gnn_ntuple = true; - // If one is not provided then throw error - if (not ana.do_write_ntuple) { - std::cout << options.help() << std::endl; - std::cout << "ERROR: option string --write_ntuple 1 and --gnn_ntuple must be set at the same time!" << std::endl; - exit(1); - } - } else { - ana.gnn_ntuple = false; - } - //_______________________________________________________________________________ // --tc_pls_triplets ana.tc_pls_triplets = result["tc_pls_triplets"].as(); @@ -264,8 +257,56 @@ int main(int argc, char **argv) { ana.no_pls_dupclean = result["no_pls_dupclean"].as(); //_______________________________________________________________________________ - // --jet_branches - ana.jet_branches = result["jet_branches"].as(); + // --md + ana.md_branches = result["md"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --ls + ana.ls_branches = result["ls"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --t3 + ana.t3_branches = result["t3"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --t5 + ana.t5_branches = result["t5"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --pls + ana.pls_branches = result["pls"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --pt3 + ana.pt3_branches = result["pt3"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --pt5 + ana.pt5_branches = result["pt5"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --occ + ana.occ_branches = result["occ"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --t5dnn + ana.t5dnn_branches = result["t5dnn"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --t3dnn + ana.t3dnn_branches = result["t3dnn"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --pt3dnn + ana.pt3dnn_branches = result["pt3dnn"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --jet + ana.jet_branches = result["jet"].as() || result["allobj"].as(); + + //_______________________________________________________________________________ + // --sim + ana.extra_sim_branches = result["sim"].as() || result["allobj"].as(); // Printing out the option settings overview std::cout << "=========================================================" << std::endl; @@ -330,9 +371,6 @@ void run_lst() { if (ana.do_write_ntuple) { createOutputBranches(); - if (ana.gnn_ntuple) { - createGnnNtupleBranches(); - } } std::vector out_lstInputHC; diff --git a/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h b/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h index 7fc4713a2dc61..3b712d25ab705 100644 --- a/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h +++ b/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h @@ -114,9 +114,6 @@ class AnalysisConfig { // Boolean to write lower level objects bool do_lower_level; - // Boolean to write gnn ntuple - bool gnn_ntuple; - // String to hold the MAKETARGET setting from compile std::string compilation_target; @@ -132,8 +129,44 @@ class AnalysisConfig { // Boolean to disable pLS duplicate cleaning bool no_pls_dupclean; + // Boolean to enable MD branches + bool md_branches; + + // Boolean to enable LS branches + bool ls_branches; + + // Boolean to enable T3 branches + bool t3_branches; + + // Boolean to enable T5 branches + bool t5_branches; + + // Boolean to enable pLS branches + bool pls_branches; + + // Boolean to enable pT3 branches + bool pt3_branches; + + // Boolean to enable pT5 branches + bool pt5_branches; + + // Boolean to enable occupancy branches + bool occ_branches; + + // Boolean to enable T3 DNN branches + bool t3dnn_branches; + + // Boolean to enable T5 DNN branches + bool t5dnn_branches; + + // Boolean to enable pT3 DNN branches + bool pt3dnn_branches; + // Boolean to enable jet branches bool jet_branches; + + // Boolean to enable extra sim branches + bool extra_sim_branches; }; extern AnalysisConfig ana; diff --git a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc index d05d28f431ee5..633b223ce425d 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc @@ -289,7 +289,25 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector> const& trk_ph2_simHitIdx, std::vector> const& trk_pix_simHitIdx, bool verbose, + float matchfrac, float* pmatched) { + std::vector matched_idxs; + std::vector matched_idx_fracs; + std::tie(matched_idxs, matched_idx_fracs) = matchedSimTrkIdxsAndFracs( + hitidxs, hittypes, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, verbose, matchfrac, pmatched); + return matched_idxs; +} + +//___________________________________________________________________________________________________________________________________________________________________________________________ +std::tuple, std::vector> matchedSimTrkIdxsAndFracs( + std::vector hitidxs, + std::vector hittypes, + std::vector const& trk_simhit_simTrkIdx, + std::vector> const& trk_ph2_simHitIdx, + std::vector> const& trk_pix_simHitIdx, + bool verbose, + float matchfrac, + float* pmatched) { if (hitidxs.size() != hittypes.size()) { std::cout << "Error: matched_sim_trk_idxs() hitidxs and hittypes have different lengths" << std::endl; std::cout << "hitidxs.size(): " << hitidxs.size() << std::endl; @@ -395,11 +413,8 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector intermediate, size_t n, std::vector>& va) { - // std::cout << " 'called': " << "called" << std::endl; if (va.size() > n) { for (auto x : va[n]) { - // std::cout << " n: " << n << std::endl; - // std::cout << " intermediate.size(): " << intermediate.size() << std::endl; std::vector copy_intermediate(intermediate); copy_intermediate.push_back(x); perm(result, copy_intermediate, n + 1, va); @@ -423,6 +438,7 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, } int maxHitMatchCount = 0; // ultimate maximum of the number of matched hits std::vector matched_sim_trk_idxs; + std::vector matched_sim_trk_idxs_frac; float max_percent_matched = 0.0f; for (auto& trkidx_perm : allperms) { std::vector counts; @@ -436,8 +452,13 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, if (trkidx < 0) continue; float percent_matched = static_cast(counts[rawidx]) / nhits_input; - if (percent_matched > 0.75f) + if (verbose) { + std::cout << " fr: " << percent_matched << std::endl; + } + if (percent_matched > matchfrac) { matched_sim_trk_idxs.push_back(trkidx); + matched_sim_trk_idxs_frac.push_back(percent_matched); + } maxHitMatchCount = std::max(maxHitMatchCount, *std::max_element(counts.begin(), counts.end())); max_percent_matched = std::max(max_percent_matched, percent_matched); } @@ -447,12 +468,35 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, *pmatched = max_percent_matched; } - std::set s; + std::map pairs; unsigned size = matched_sim_trk_idxs.size(); - for (unsigned i = 0; i < size; ++i) - s.insert(matched_sim_trk_idxs[i]); - matched_sim_trk_idxs.assign(s.begin(), s.end()); - return matched_sim_trk_idxs; + for (unsigned i = 0; i < size; ++i) { + int idx = matched_sim_trk_idxs[i]; + float frac = matched_sim_trk_idxs_frac[i]; + if (pairs.find(idx) != pairs.end()) { + if (pairs[idx] < frac) + pairs[idx] = frac; + } else { + pairs[idx] = frac; + } + } + std::vector result; + std::vector result_frac; + // Loop over the map using range-based for loop + for (const auto& pair : pairs) { + result.push_back(pair.first); + result_frac.push_back(pair.second); + } + // Sort indices based on 'values' + auto indices = sort_indices(result_frac); + // Reorder 'vec1' and 'vec2' based on the sorted indices + std::vector sorted_result(result.size()); + std::vector sorted_result_frac(result_frac.size()); + for (size_t i = 0; i < indices.size(); ++i) { + sorted_result[i] = result[indices[i]]; + sorted_result_frac[i] = result_frac[indices[i]]; + } + return std::make_tuple(sorted_result, sorted_result_frac); } //__________________________________________________________________________________________ diff --git a/RecoTracker/LSTCore/standalone/code/core/trkCore.h b/RecoTracker/LSTCore/standalone/code/core/trkCore.h index d4938490053be..e1d572bfcbb49 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.h +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.h @@ -39,7 +39,17 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector> const& trk_ph2_simHitIdx, std::vector> const& trk_pix_simHitIdx, bool verbose = false, + float matchfrac = 0.75, float* pmatched = nullptr); +std::tuple, std::vector> matchedSimTrkIdxsAndFracs( + std::vector hitidxs, + std::vector hittypes, + std::vector const& trk_simhit_simTrkIdx, + std::vector> const& trk_ph2_simHitIdx, + std::vector> const& trk_pix_simHitIdx, + bool verbose = false, + float matchfrac = 0.75, + float* pmatched = nullptr); int getDenomSimTrkType(int isimtrk, std::vector const& trk_sim_q, std::vector const& trk_sim_pt, @@ -113,6 +123,14 @@ void printTimingInformation(std::vector>& timing_information, TString get_absolute_path_after_check_file_exists(const std::string name); void writeMetaData(); +template +std::vector sort_indices(const std::vector& vec) { + std::vector indices(vec.size()); + std::iota(indices.begin(), indices.end(), 0); + std::sort(indices.begin(), indices.end(), [&vec](size_t i1, size_t i2) { return vec[i1] > vec[i2]; }); + return indices; +} + // --------------------- ======================== --------------------- // DEPRECATED FUNCTION diff --git a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc index dba2886716d1c..2ec878ccbf80a 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc +++ b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc @@ -7,16 +7,73 @@ using namespace ALPAKA_ACCELERATOR_NAMESPACE::lst; //________________________________________________________________________________________________________________________________ void createOutputBranches() { - createRequiredOutputBranches(); - createOptionalOutputBranches(); + createSimTrackContainerBranches(); + createTrackCandidateBranches(); + + if (ana.jet_branches) + createJetBranches(); + + if (ana.md_branches) + createMiniDoubletBranches(); + if (ana.ls_branches) + createLineSegmentBranches(); + if (ana.t3_branches) + createTripletBranches(); + if (ana.t5_branches) + createQuintupletBranches(); + if (ana.pls_branches) + createPixelLineSegmentBranches(); + if (ana.pt3_branches) + createPixelTripletBranches(); + if (ana.pt5_branches) + createPixelQuintupletBranches(); + + if (ana.occ_branches) + createOccupancyBranches(); + + // DNN branches + if (ana.t5dnn_branches) + createT5DNNBranches(); + if (ana.t3dnn_branches) + createT3DNNBranches(); + if (ana.pt3dnn_branches) + createpT3DNNBranches(); } //________________________________________________________________________________________________________________________________ void fillOutputBranches(LSTEvent* event) { - setOutputBranches(event); - setOptionalOutputBranches(event); - if (ana.gnn_ntuple) - setGnnNtupleBranches(event); + float matchfrac = 0.75; + + unsigned int n_accepted_simtrk = setSimTrackContainerBranches(event); + + if (ana.occ_branches) + setOccupancyBranches(event); + + if (ana.t3dnn_branches) + setT3DNNBranches(event, matchfrac); + if (ana.t5dnn_branches) + setT5DNNBranches(event); + if (ana.pt3dnn_branches) + setpT3DNNBranches(event); + + auto const md_idx_map = (ana.md_branches ? setMiniDoubletBranches(event, n_accepted_simtrk, matchfrac) + : std::map()); + auto const ls_idx_map = (ana.ls_branches ? setLineSegmentBranches(event, n_accepted_simtrk, matchfrac, md_idx_map) + : std::map()); + auto const t3_idx_map = (ana.t3_branches ? setTripletBranches(event, n_accepted_simtrk, matchfrac, ls_idx_map) + : std::map()); + auto const t5_idx_map = (ana.t5_branches ? setQuintupletBranches(event, n_accepted_simtrk, matchfrac, t3_idx_map) + : std::map()); + auto const pls_idx_map = (ana.pls_branches ? setPixelLineSegmentBranches(event, n_accepted_simtrk, matchfrac) + : std::map()); + auto const pt3_idx_map = + (ana.pt3_branches ? setPixelTripletBranches(event, n_accepted_simtrk, matchfrac, pls_idx_map, t3_idx_map) + : std::map()); + auto const pt5_idx_map = + (ana.pt5_branches ? setPixelQuintupletBranches(event, n_accepted_simtrk, matchfrac, pls_idx_map, t5_idx_map) + : std::map()); + + setTrackCandidateBranches(event, n_accepted_simtrk, t5_idx_map, pls_idx_map, pt3_idx_map, pt5_idx_map, matchfrac); // Now actually fill the ttree ana.tx->fill(); @@ -25,165 +82,6 @@ void fillOutputBranches(LSTEvent* event) { ana.tx->clear(); } -//________________________________________________________________________________________________________________________________ -void createRequiredOutputBranches() { - // Setup output TTree - - if (ana.jet_branches) { - ana.tx->createBranch>("sim_deltaEta"); - ana.tx->createBranch>("sim_deltaPhi"); - ana.tx->createBranch>("sim_deltaR"); - ana.tx->createBranch>("sim_jet_eta"); - ana.tx->createBranch>("sim_jet_phi"); - ana.tx->createBranch>("sim_jet_pt"); - } - ana.tx->createBranch>("sim_pt"); - ana.tx->createBranch>("sim_eta"); - ana.tx->createBranch>("sim_phi"); - ana.tx->createBranch>("sim_pca_dxy"); - ana.tx->createBranch>("sim_pca_dz"); - ana.tx->createBranch>("sim_q"); - ana.tx->createBranch>("sim_event"); - ana.tx->createBranch>("sim_pdgId"); - ana.tx->createBranch>("sim_vx"); - ana.tx->createBranch>("sim_vy"); - ana.tx->createBranch>("sim_vz"); - ana.tx->createBranch>("sim_trkNtupIdx"); - ana.tx->createBranch>("sim_TC_matched"); - ana.tx->createBranch>("sim_TC_matched_mask"); - - // Track candidates - ana.tx->createBranch>("tc_pt"); - ana.tx->createBranch>("tc_eta"); - ana.tx->createBranch>("tc_phi"); - ana.tx->createBranch>("tc_type"); - ana.tx->createBranch>("tc_isFake"); - ana.tx->createBranch>("tc_isDuplicate"); - ana.tx->createBranch>>("tc_matched_simIdx"); -} - -//________________________________________________________________________________________________________________________________ -void createOptionalOutputBranches() { -#ifdef CUT_VALUE_DEBUG - // Event-wide branches - // ana.tx->createBranch("evt_dummy"); - - // Sim Track branches - // NOTE: Must sync with main tc branch in length!! - ana.tx->createBranch>("sim_dummy"); - - // Track Candidate branches - // NOTE: Must sync with main tc branch in length!! - ana.tx->createBranch>("tc_dummy"); - - // pT5 branches - ana.tx->createBranch>>("pT5_matched_simIdx"); - ana.tx->createBranch>>("pT5_hitIdxs"); - ana.tx->createBranch>("sim_pT5_matched"); - ana.tx->createBranch>("pT5_pt"); - ana.tx->createBranch>("pT5_eta"); - ana.tx->createBranch>("pT5_phi"); - ana.tx->createBranch>("pT5_isFake"); - ana.tx->createBranch>("t5_sim_vxy"); - ana.tx->createBranch>("t5_sim_vz"); - ana.tx->createBranch>("pT5_isDuplicate"); - ana.tx->createBranch>("pT5_score"); - ana.tx->createBranch>("pT5_layer_binary"); - ana.tx->createBranch>("pT5_moduleType_binary"); - ana.tx->createBranch>("pT5_matched_pt"); - ana.tx->createBranch>("pT5_rzChiSquared"); - ana.tx->createBranch>("pT5_rPhiChiSquared"); - ana.tx->createBranch>("pT5_rPhiChiSquaredInwards"); - - // pT3 branches - ana.tx->createBranch>("sim_pT3_matched"); - ana.tx->createBranch>("pT3_pt"); - ana.tx->createBranch>("pT3_isFake"); - ana.tx->createBranch>("pT3_isDuplicate"); - ana.tx->createBranch>("pT3_eta"); - ana.tx->createBranch>("pT3_phi"); - ana.tx->createBranch>("pT3_score"); - ana.tx->createBranch>("pT3_foundDuplicate"); - ana.tx->createBranch>>("pT3_matched_simIdx"); - ana.tx->createBranch>>("pT3_hitIdxs"); - ana.tx->createBranch>("pT3_pixelRadius"); - ana.tx->createBranch>("pT3_pixelRadiusError"); - ana.tx->createBranch>>("pT3_matched_pt"); - ana.tx->createBranch>("pT3_tripletRadius"); - ana.tx->createBranch>("pT3_rPhiChiSquared"); - ana.tx->createBranch>("pT3_rPhiChiSquaredInwards"); - ana.tx->createBranch>("pT3_rzChiSquared"); - ana.tx->createBranch>("pT3_layer_binary"); - ana.tx->createBranch>("pT3_moduleType_binary"); - - // pLS branches - ana.tx->createBranch>("sim_pLS_matched"); - ana.tx->createBranch>>("pLS_matched_simIdx"); - ana.tx->createBranch>("pLS_isFake"); - ana.tx->createBranch>("pLS_isDuplicate"); - ana.tx->createBranch>("pLS_ptIn"); - ana.tx->createBranch>("pLS_ptErr"); - ana.tx->createBranch>("pLS_px"); - ana.tx->createBranch>("pLS_py"); - ana.tx->createBranch>("pLS_pz"); - ana.tx->createBranch>("pLS_eta"); - ana.tx->createBranch>("pLS_isQuad"); - ana.tx->createBranch>("pLS_etaErr"); - ana.tx->createBranch>("pLS_phi"); - ana.tx->createBranch>("pLS_score"); - ana.tx->createBranch>("pLS_circleCenterX"); - ana.tx->createBranch>("pLS_circleCenterY"); - ana.tx->createBranch>("pLS_circleRadius"); - - // T5 branches - ana.tx->createBranch>("sim_T5_matched"); - ana.tx->createBranch>("t5_isFake"); - ana.tx->createBranch>("t5_isDuplicate"); - ana.tx->createBranch>("t5_foundDuplicate"); - ana.tx->createBranch>("t5_pt"); - ana.tx->createBranch>("t5_pMatched"); - ana.tx->createBranch>("t5_eta"); - ana.tx->createBranch>("t5_phi"); - ana.tx->createBranch>("t5_score_rphisum"); - ana.tx->createBranch>>("t5_hitIdxs"); - ana.tx->createBranch>>("t5_matched_simIdx"); - ana.tx->createBranch>("t5_moduleType_binary"); - ana.tx->createBranch>("t5_layer_binary"); - ana.tx->createBranch>("t5_matched_pt"); - ana.tx->createBranch>("t5_innerRadius"); - ana.tx->createBranch>("t5_outerRadius"); - ana.tx->createBranch>("t5_bridgeRadius"); - ana.tx->createBranch>("t5_chiSquared"); - ana.tx->createBranch>("t5_rzChiSquared"); - ana.tx->createBranch>("t5_isDupAlgoFlag"); - ana.tx->createBranch>("t5_nonAnchorChiSquared"); - ana.tx->createBranch>("t5_dBeta1"); - ana.tx->createBranch>("t5_dBeta2"); - - // Occupancy branches - ana.tx->createBranch>("module_layers"); - ana.tx->createBranch>("module_subdets"); - ana.tx->createBranch>("module_rings"); - ana.tx->createBranch>("module_rods"); - ana.tx->createBranch>("module_modules"); - ana.tx->createBranch>("module_isTilted"); - ana.tx->createBranch>("module_eta"); - ana.tx->createBranch>("module_r"); - ana.tx->createBranch>("md_occupancies"); - ana.tx->createBranch>("sg_occupancies"); - ana.tx->createBranch>("t3_occupancies"); - ana.tx->createBranch("tc_occupancies"); - ana.tx->createBranch>("t5_occupancies"); - ana.tx->createBranch("pT3_occupancies"); - ana.tx->createBranch("pT5_occupancies"); - - // T5 DNN branches - createT5DNNBranches(); - createT3DNNBranches(); - -#endif -} - //________________________________________________________________________________________________________________________________ void createT5DNNBranches() { // Common branches @@ -217,6 +115,16 @@ void createT5DNNBranches() { } } +//________________________________________________________________________________________________________________________________ +void createpT3DNNBranches() { + ana.tx->createBranch>("pT3_pixelRadius"); + ana.tx->createBranch>("pT3_pixelRadiusError"); + ana.tx->createBranch>("pT3_tripletRadius"); + ana.tx->createBranch>("pT3_rPhiChiSquared"); + ana.tx->createBranch>("pT3_rPhiChiSquaredInwards"); + ana.tx->createBranch>("pT3_rzChiSquared"); +} + //________________________________________________________________________________________________________________________________ void createT3DNNBranches() { // Common branches for T3 properties based on TripletsSoA fields @@ -252,59 +160,380 @@ void createT3DNNBranches() { } //________________________________________________________________________________________________________________________________ -void createGnnNtupleBranches() { - // Mini Doublets - ana.tx->createBranch>("MD_pt"); - ana.tx->createBranch>("MD_eta"); - ana.tx->createBranch>("MD_phi"); - ana.tx->createBranch>("MD_dphichange"); - ana.tx->createBranch>("MD_isFake"); - ana.tx->createBranch>("MD_tpType"); - ana.tx->createBranch>("MD_detId"); - ana.tx->createBranch>("MD_layer"); - ana.tx->createBranch>("MD_0_r"); - ana.tx->createBranch>("MD_0_x"); - ana.tx->createBranch>("MD_0_y"); - ana.tx->createBranch>("MD_0_z"); - ana.tx->createBranch>("MD_1_r"); - ana.tx->createBranch>("MD_1_x"); - ana.tx->createBranch>("MD_1_y"); - ana.tx->createBranch>("MD_1_z"); +void createJetBranches() { + ana.tx->createBranch>("sim_deltaEta"); + ana.tx->createBranch>("sim_deltaPhi"); + ana.tx->createBranch>("sim_deltaR"); + ana.tx->createBranch>("sim_jet_eta"); + ana.tx->createBranch>("sim_jet_phi"); + ana.tx->createBranch>("sim_jet_pt"); +} - // Line Segments - ana.tx->createBranch>("LS_pt"); - ana.tx->createBranch>("LS_eta"); - ana.tx->createBranch>("LS_phi"); - ana.tx->createBranch>("LS_isFake"); - ana.tx->createBranch>("LS_MD_idx0"); - ana.tx->createBranch>("LS_MD_idx1"); - ana.tx->createBranch>("LS_sim_pt"); - ana.tx->createBranch>("LS_sim_eta"); - ana.tx->createBranch>("LS_sim_phi"); - ana.tx->createBranch>("LS_sim_pca_dxy"); - ana.tx->createBranch>("LS_sim_pca_dz"); - ana.tx->createBranch>("LS_sim_q"); - ana.tx->createBranch>("LS_sim_pdgId"); - ana.tx->createBranch>("LS_sim_event"); - ana.tx->createBranch>("LS_sim_bx"); - ana.tx->createBranch>("LS_sim_vx"); - ana.tx->createBranch>("LS_sim_vy"); - ana.tx->createBranch>("LS_sim_vz"); - ana.tx->createBranch>("LS_isInTrueTC"); - - // TC's LS - ana.tx->createBranch>>("tc_lsIdx"); +//________________________________________________________________________________________________________________________________ +void createSimTrackContainerBranches() { + // Simulated Track Container + // + // The container will hold per entry a simulated track in the event. Only the current bunch crossing, and + // primary vertex (hard-scattered) tracks will be saved to reduce the size of the output. + // + ana.tx->createBranch>("sim_pt"); // pt + ana.tx->createBranch>("sim_eta"); // eta + ana.tx->createBranch>("sim_phi"); // phi + ana.tx->createBranch>("sim_pca_dxy"); // dxy of point of closest approach + ana.tx->createBranch>("sim_pca_dz"); // dz of point of clossest approach + ana.tx->createBranch>("sim_q"); // charge +1, -1, 0 + ana.tx->createBranch>("sim_pdgId"); // pdgId + // production vertex x position (values are derived from simvtx_* and sim_parentVtxIdx branches in the tracking ntuple) + ana.tx->createBranch>("sim_vx"); + // production vertex y position (values are derived from simvtx_* and sim_parentVtxIdx branches in the tracking ntuple) + ana.tx->createBranch>("sim_vy"); + // production vertex z position (values are derived from simvtx_* and sim_parentVtxIdx branches in the tracking ntuple) + ana.tx->createBranch>("sim_vz"); + // production vertex r (sqrt(x**2 + y**2)) position (values are derived from simvtx_* and sim_parentVtxIdx branches in the tracking ntuple) + ana.tx->createBranch>("sim_vtxperp"); + // idx of sim_* in the tracking ntuple (N.B. this may be redundant) + ana.tx->createBranch>("sim_trkNtupIdx"); + // idx to the best match (highest nhit match) tc_* container + ana.tx->createBranch>("sim_tcIdxBest"); + // match fraction to the best match (highest nhit match) tc_* container + ana.tx->createBranch>("sim_tcIdxBestFrac"); + // idx to the best match (highest nhit match and > 75%) tc_* container + ana.tx->createBranch>("sim_tcIdx"); + // list of idx to any matches (> 0%) to tc_* container + ana.tx->createBranch>>("sim_tcIdxAll"); + // list of match fraction for each match (> 0%) to tc_* container + ana.tx->createBranch>>("sim_tcIdxAllFrac"); + + if (ana.extra_sim_branches) { + ana.tx->createBranch>>("sim_simHitX"); // list of simhit's X positions + ana.tx->createBranch>>("sim_simHitY"); // list of simhit's Y positions + ana.tx->createBranch>>("sim_simHitZ"); // list of simhit's Z positions + ana.tx->createBranch>>("sim_simHitDetId"); // list of simhit's detId + // list of simhit's layers (N.B. layer is numbered 1 2 3 4 5 6 for barrel, 7 8 9 10 11 for endcaps) + ana.tx->createBranch>>("sim_simHitLayer"); + // list of simhit's distance in xy-plane to the expected point based on simhit's z position and helix formed from pt,eta,phi,vx,vy,vz,q of the simulated track + ana.tx->createBranch>>("sim_simHitDistxyHelix"); + // length of 11 float numbers with min(simHitDistxyHelix) value for each layer. Useful for finding e.g. "sim tracks that traversed barrel detector entirelyand left a reasonable hit in layer 1 2 3 4 5 6 layers." + ana.tx->createBranch>>("sim_simHitLayerMinDistxyHelix"); + // length of 11 float numbers with min(simHitDistxyHelix) value for each layer. Useful for finding e.g. "sim tracks that traversed barrel detector entirelyand left a reasonable hit in layer 1 2 3 4 5 6 layers." + ana.tx->createBranch>>("sim_simHitLayerMinDistxyPrevHit"); + ana.tx->createBranch>>("sim_recoHitX"); // list of recohit's X positions + ana.tx->createBranch>>("sim_recoHitY"); // list of recohit's Y positions + ana.tx->createBranch>>("sim_recoHitZ"); // list of recohit's Z positions + ana.tx->createBranch>>("sim_recoHitDetId"); // list of recohit's detId + } + + if (ana.md_branches) { + // list of idx to matches (> 0%) to md_* container + ana.tx->createBranch>>("sim_mdIdxAll"); + // list of match fraction for each match (> 0%) to md_* container + ana.tx->createBranch>>("sim_mdIdxAllFrac"); + } + if (ana.ls_branches) { + // list of idx to matches (> 0%) to ls_* container + ana.tx->createBranch>>("sim_lsIdxAll"); + // list of match fraction for each match (> 0%) to ls_* container + ana.tx->createBranch>>("sim_lsIdxAllFrac"); + } + if (ana.t3_branches) { + // list of idx to matches (> 0%) to t3_* container + ana.tx->createBranch>>("sim_t3IdxAll"); + // list of match fraction for each match (> 0%) to t3_* container + ana.tx->createBranch>>("sim_t3IdxAllFrac"); + } + if (ana.t5_branches) { + // list of idx to matches (> 0%) to t5_* container + ana.tx->createBranch>>("sim_t5IdxAll"); + // list of match fraction for each match (> 0%) to t5_* container + ana.tx->createBranch>>("sim_t5IdxAllFrac"); + } + if (ana.pls_branches) { + // list of idx to matches (> 0%) to pls_* container + ana.tx->createBranch>>("sim_plsIdxAll"); + // list of match fraction for each match (> 0%) to pls_* container + ana.tx->createBranch>>("sim_plsIdxAllFrac"); + } + if (ana.pt3_branches) { + // list of idx to matches (> 0%) to pt3_* container + ana.tx->createBranch>>("sim_pt3IdxAll"); + // list of match fraction for each match (> 0%) to pt3_* container + ana.tx->createBranch>>("sim_pt3IdxAllFrac"); + } + if (ana.pt5_branches) { + // list of idx to matches (> 0%) to pt5_* container + ana.tx->createBranch>>("sim_pt5IdxAll"); + // list of match fraction for each match (> 0%) to pt5_* container + ana.tx->createBranch>>("sim_pt5IdxAllFrac"); + } } //________________________________________________________________________________________________________________________________ -void setOutputBranches(LSTEvent* event) { - // ============ Sim tracks ============= - int n_accepted_simtrk = 0; - auto const& trk_sim_bunchCrossing = trk.getVI("sim_bunchCrossing"); - auto const& trk_sim_event = trk.getVI("sim_event"); +void createTrackCandidateBranches() { + // Track Candidates + // + // The container will hold per entry a track candidate built by LST in the event. + // + ana.tx->createBranch>("tc_pt"); // pt + ana.tx->createBranch>("tc_eta"); // eta + ana.tx->createBranch>("tc_phi"); // phi + ana.tx->createBranch>("tc_type"); // type = 7 (pT5), 5 (pT3), 4 (T5), 8 (pLS) + ana.tx->createBranch>("tc_isFake"); // 1 if tc is fake 0 other if not + ana.tx->createBranch>("tc_isDuplicate"); // 1 if tc is duplicate 0 other if not + ana.tx->createBranch>("tc_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("tc_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("tc_simIdxAllFrac"); + if (ana.pt5_branches) + ana.tx->createBranch>( + "tc_pt5Idx"); // index to the pt5_* if it is the said type, if not set to -999 + if (ana.pt3_branches) + ana.tx->createBranch>( + "tc_pt3Idx"); // index to the pt3_* if it is the said type, if not set to -999 + if (ana.t5_branches) + ana.tx->createBranch>( + "tc_t5Idx"); // index to the t5_* if it is the said type, if not set to -999 + if (ana.pls_branches) + ana.tx->createBranch>( + "tc_plsIdx"); // index to the pls_* if it is the said type, if not set to -999 +} + +//________________________________________________________________________________________________________________________________ +void createMiniDoubletBranches() { + // Mini-Doublets (i.e. Two reco hits paired in a single pT-module of Outer Tracker of CMS, a.k.a. MD) + // + // The container will hold per entry a mini-doublet built by LST in the event. + // + ana.tx->createBranch>("md_pt"); // pt (computed based on delta phi change) + ana.tx->createBranch>("md_eta"); // eta (computed based on anchor hit's eta) + ana.tx->createBranch>("md_phi"); // phi (computed based on anchor hit's phi) +#ifdef CUT_VALUE_DEBUG + ana.tx->createBranch>("md_dphi"); + ana.tx->createBranch>("md_dphichange"); + ana.tx->createBranch>("md_dz"); +#endif + ana.tx->createBranch>("md_anchor_x"); // anchor hit x + ana.tx->createBranch>("md_anchor_y"); // anchor hit y + ana.tx->createBranch>("md_anchor_z"); // anchor hit z + ana.tx->createBranch>("md_other_x"); // other hit x + ana.tx->createBranch>("md_other_y"); // other hit y + ana.tx->createBranch>("md_other_z"); // other hit z + // type of the module where the mini-doublet sit (type = 1 (PS), 0 (2S)) + ana.tx->createBranch>("md_type"); + // layer index of the module where the mini-doublet sit (layer = 1 2 3 4 5 6 (barrel) 7 8 9 10 11 (endcap)) + ana.tx->createBranch>("md_layer"); + // detId = detector unique ID that contains a lot of information that can be parsed later if needed + ana.tx->createBranch>("md_detId"); + ana.tx->createBranch>("md_isFake"); // 1 if md is fake 0 other if not + ana.tx->createBranch>("md_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("md_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("md_simIdxAllFrac"); +} + +//________________________________________________________________________________________________________________________________ +void createLineSegmentBranches() { + // Line Segments (i.e. Two mini-doublets, a.k.a. LS) + // + // The container will hold per entry a line-segment built by LST in the event. + // + // pt (computed based on radius of the circle formed by three points: (origin), (anchor hit 1), (anchor hit 2)) + ana.tx->createBranch>("ls_pt"); + ana.tx->createBranch>("ls_eta"); // eta (computed based on last anchor hit's eta) + ana.tx->createBranch>("ls_phi"); // phi (computed based on first anchor hit's phi) + ana.tx->createBranch>("ls_mdIdx0"); // index to the first MD + ana.tx->createBranch>("ls_mdIdx1"); // index to the second MD + ana.tx->createBranch>("ls_isFake"); // 1 if md is fake 0 other if not + ana.tx->createBranch>("ls_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track +#ifdef CUT_VALUE_DEBUG + ana.tx->createBranch>("ls_zLos"); + ana.tx->createBranch>("ls_zHis"); + ana.tx->createBranch>("ls_rtLos"); + ana.tx->createBranch>("ls_rtHis"); + ana.tx->createBranch>("ls_dPhis"); + ana.tx->createBranch>("ls_dPhiMins"); + ana.tx->createBranch>("ls_dPhiMaxs"); + ana.tx->createBranch>("ls_dPhiChanges"); + ana.tx->createBranch>("ls_dPhiChangeMins"); + ana.tx->createBranch>("ls_dPhiChangeMaxs"); + ana.tx->createBranch>("ls_dAlphaInners"); + ana.tx->createBranch>("ls_dAlphaOuters"); + ana.tx->createBranch>("ls_dAlphaInnerOuters"); +#endif + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("ls_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("ls_simIdxAllFrac"); +} + +//________________________________________________________________________________________________________________________________ +void createTripletBranches() { + // Triplets (i.e. Three mini-doublets, a.k.a. T3) + // + // The container will hold per entry a triplets built by LST in the event. + // + // pt (computed based on radius of the circle formed by three points: anchor hit 1, 2, 3 + ana.tx->createBranch>("t3_pt"); + ana.tx->createBranch>("t3_eta"); // eta (computed based on last anchor hit's eta) + ana.tx->createBranch>("t3_phi"); // phi (computed based on first anchor hit's phi) + ana.tx->createBranch>("t3_lsIdx0"); // index to the first LS + ana.tx->createBranch>("t3_lsIdx1"); // index to the second LS + ana.tx->createBranch>("t3_isFake"); // 1 if t3 is fake 0 other if not + ana.tx->createBranch>("t3_isDuplicate"); // 1 if t3 is duplicate 0 other if not + ana.tx->createBranch>("t3_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("t3_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("t3_simIdxAllFrac"); +} + +//________________________________________________________________________________________________________________________________ +void createQuintupletBranches() { + // Quintuplets (i.e. Five mini-doublets, a.k.a. T5) + // + // The container will hold per entry a quintuplet built by LST in the event. + // + // pt (computed based on average of the 4 circles formed by, (1, 2, 3), (2, 3, 4), (3, 4, 5), (1, 3, 5) + ana.tx->createBranch>("t5_pt"); + ana.tx->createBranch>("t5_eta"); // eta (computed based on last anchor hit's eta) + ana.tx->createBranch>("t5_phi"); // phi (computed based on first anchor hit's phi) + ana.tx->createBranch>("t5_t3Idx0"); // index of first T3 + ana.tx->createBranch>("t5_t3Idx1"); // index of second T3 + ana.tx->createBranch>("t5_isFake"); // 1 if t5 is fake 0 other if not + ana.tx->createBranch>("t5_isDuplicate"); // 1 if t5 is duplicate 0 other if not + ana.tx->createBranch>("t5_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("t5_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("t5_simIdxAllFrac"); + ana.tx->createBranch>("t5_innerRadius"); + ana.tx->createBranch>("t5_outerRadius"); + ana.tx->createBranch>("t5_bridgeRadius"); + ana.tx->createBranch>("t5_pMatched"); + ana.tx->createBranch>("t5_sim_vxy"); + ana.tx->createBranch>("t5_sim_vz"); +} + +//________________________________________________________________________________________________________________________________ +void createPixelLineSegmentBranches() { + // Pixel Line Segments (a.k.a pLS) + // + // The container will hold per entry a pixel line segment (built by an external algo, e.g. patatrack) accepted by LST in the event. + // + // pt (taken from pt of the 3-vector from see_stateTrajGlbPx/Py/Pz) + ana.tx->createBranch>("pLS_pt"); + ana.tx->createBranch>("pLS_ptErr"); + // eta (taken from eta of the 3-vector from see_stateTrajGlbPx/Py/Pz) + ana.tx->createBranch>("pLS_eta"); + ana.tx->createBranch>("pLS_etaErr"); + // phi (taken from phi of the 3-vector from see_stateTrajGlbPx/Py/Pz) + ana.tx->createBranch>("pLS_phi"); + ana.tx->createBranch>("pLS_nhit"); // Number of actual hit: 3 if triplet, 4 if quadruplet + ana.tx->createBranch>("pLS_hit0_x"); // pLS's reco hit0 x + ana.tx->createBranch>("pLS_hit0_y"); // pLS's reco hit0 y + ana.tx->createBranch>("pLS_hit0_z"); // pLS's reco hit0 z + ana.tx->createBranch>("pLS_hit1_x"); // pLS's reco hit1 x + ana.tx->createBranch>("pLS_hit1_y"); // pLS's reco hit1 y + ana.tx->createBranch>("pLS_hit1_z"); // pLS's reco hit1 z + ana.tx->createBranch>("pLS_hit2_x"); // pLS's reco hit2 x + ana.tx->createBranch>("pLS_hit2_y"); // pLS's reco hit2 y + ana.tx->createBranch>("pLS_hit2_z"); // pLS's reco hit2 z + ana.tx->createBranch>("pLS_hit3_x"); // pLS's reco hit3 x (if triplet, this is set to -999) + ana.tx->createBranch>("pLS_hit3_y"); // pLS's reco hit3 y (if triplet, this is set to -999) + ana.tx->createBranch>("pLS_hit3_z"); // pLS's reco hit3 z (if triplet, this is set to -999) + ana.tx->createBranch>("pLS_isFake"); // 1 if pLS is fake 0 other if not + ana.tx->createBranch>("pLS_isDuplicate"); // 1 if pLS is duplicate 0 other if not + ana.tx->createBranch>("pLS_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("pLS_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("pLS_simIdxAllFrac"); + ana.tx->createBranch>("pLS_circleCenterX"); + ana.tx->createBranch>("pLS_circleCenterY"); + ana.tx->createBranch>("pLS_circleRadius"); + ana.tx->createBranch>("pLS_px"); + ana.tx->createBranch>("pLS_py"); + ana.tx->createBranch>("pLS_pz"); + ana.tx->createBranch>("pLS_isQuad"); +} + +//________________________________________________________________________________________________________________________________ +void createPixelTripletBranches() { + // pLS + T3 (i.e. an object where a pLS is linked with a T3, a.k.a. pT3) + // + // The container will hold per entry a pT3 built by LST in the event. + // + ana.tx->createBranch>("sim_pT3_matched"); + ana.tx->createBranch>("pT3_score"); + ana.tx->createBranch>("pT3_pt"); // pt (taken from the pLS) + ana.tx->createBranch>("pT3_eta"); // eta (taken from the pLS) + ana.tx->createBranch>("pT3_phi"); // phi (taken from the pLS) + ana.tx->createBranch>("pT3_plsIdx"); // idx to pLS + ana.tx->createBranch>("pT3_t3Idx"); // idx to T3 + ana.tx->createBranch>("pT3_isFake"); // 1 if pT3 is fake 0 other if not + ana.tx->createBranch>("pT3_isDuplicate"); // 1 if pT3 is duplicate 0 other if not + ana.tx->createBranch>("pT3_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("pT3_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("pT3_simIdxAllFrac"); +} + +//________________________________________________________________________________________________________________________________ +void createPixelQuintupletBranches() { + // pLS + T5 (i.e. an object where a pLS is linked with a T5, a.k.a. pT5) + // + // The container will hold per entry a pT5 built by LST in the event. + // + ana.tx->createBranch>("pT5_pt"); // pt (taken from the pLS) + ana.tx->createBranch>("pT5_eta"); // eta (taken from the pLS) + ana.tx->createBranch>("pT5_phi"); // phi (taken from the pLS) + ana.tx->createBranch>("pT5_plsIdx"); // idx to pLS + ana.tx->createBranch>("pT5_t5Idx"); // idx to T5 + ana.tx->createBranch>("pT5_isFake"); // 1 if pT5 is fake 0 other if not + ana.tx->createBranch>("pT5_isDuplicate"); // 1 if pT5 is duplicate 0 other if not + ana.tx->createBranch>("pT5_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("pT5_simIdxAll"); + // list of idx of all matched (> 0%) simulated track + ana.tx->createBranch>>("pT5_simIdxAllFrac"); +} + +//________________________________________________________________________________________________________________________________ +void createOccupancyBranches() { + ana.tx->createBranch>("module_layers"); + ana.tx->createBranch>("module_subdets"); + ana.tx->createBranch>("module_rings"); + ana.tx->createBranch>("module_rods"); + ana.tx->createBranch>("module_modules"); + ana.tx->createBranch>("module_isTilted"); + ana.tx->createBranch>("module_eta"); + ana.tx->createBranch>("module_r"); + ana.tx->createBranch>("md_occupancies"); + ana.tx->createBranch>("sg_occupancies"); + ana.tx->createBranch>("t3_occupancies"); + ana.tx->createBranch("tc_occupancies"); + ana.tx->createBranch>("t5_occupancies"); + ana.tx->createBranch("pT3_occupancies"); + ana.tx->createBranch("pT5_occupancies"); +} + +//________________________________________________________________________________________________________________________________ +unsigned int setSimTrackContainerBranches(LSTEvent* event) { + //-------------------------------------------- + // + // + // Sim Tracks + // + // + //-------------------------------------------- + auto const& trk_sim_pt = trk.getVF("sim_pt"); auto const& trk_sim_eta = trk.getVF("sim_eta"); auto const& trk_sim_phi = trk.getVF("sim_phi"); + auto const& trk_sim_bunchCrossing = trk.getVI("sim_bunchCrossing"); + auto const& trk_sim_event = trk.getVI("sim_event"); auto const& trk_sim_pca_dxy = trk.getVF("sim_pca_dxy"); auto const& trk_sim_pca_dz = trk.getVF("sim_pca_dz"); auto const& trk_sim_q = trk.getVI("sim_q"); @@ -313,29 +542,42 @@ void setOutputBranches(LSTEvent* event) { auto const& trk_simvtx_x = trk.getVF("simvtx_x"); auto const& trk_simvtx_y = trk.getVF("simvtx_y"); auto const& trk_simvtx_z = trk.getVF("simvtx_z"); + auto const& trk_sim_simHitIdx = trk.getVVI("sim_simHitIdx"); + auto const& trk_simhit_subdet = trk.getVUS("simhit_subdet"); + auto const& trk_simhit_layer = trk.getVUS("simhit_layer"); + auto const& trk_simhit_x = trk.getVF("simhit_x"); + auto const& trk_simhit_y = trk.getVF("simhit_y"); + auto const& trk_simhit_z = trk.getVF("simhit_z"); + auto const& trk_simhit_detId = trk.getVU("simhit_detId"); + auto const& trk_simhit_hitIdx = trk.getVVI("simhit_hitIdx"); auto const& trk_ph2_x = trk.getVF("ph2_x"); auto const& trk_ph2_y = trk.getVF("ph2_y"); auto const& trk_ph2_z = trk.getVF("ph2_z"); - auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); - auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); - auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + auto const& trk_ph2_detId = trk.getVU("ph2_detId"); + + // Total number of simulated tracks with the condition that the simulated track came from a particle produced in the hard scattering and from the current bunch-crossing) + // "accepted" here would mean that in the tracking ntuple (sim_bunchCrossing == 0 and sim_event == 0) + unsigned int n_accepted_simtrk = 0; - if (ana.jet_branches) { - auto const& trk_sim_deltaEta = trk.getVF("sim_deltaEta"); - auto const& trk_sim_deltaPhi = trk.getVF("sim_deltaPhi"); - auto const& trk_sim_deltaR = trk.getVF("sim_deltaR"); - auto const& trk_sim_jet_eta = trk.getVF("sim_jet_eta"); - auto const& trk_sim_jet_phi = trk.getVF("sim_jet_phi"); - auto const& trk_sim_jet_pt = trk.getVF("sim_jet_pt"); + // Looping over the simulated tracks in the tracking ntuple + for (unsigned int isimtrk = 0; isimtrk < trk_sim_pt.size(); ++isimtrk) { + // Skip out-of-time pileup + if (trk_sim_bunchCrossing[isimtrk] != 0) + continue; + + // Skip non-hard-scatter + if (trk_sim_event[isimtrk] != 0) + continue; - for (unsigned int isimtrk = 0; isimtrk < trk_sim_pt.size(); ++isimtrk) { - // Skip out-of-time pileup - if (trk_sim_bunchCrossing[isimtrk] != 0) - continue; + // Now we have a list of "accepted" tracks (no condition on vtx_z/perp, nor pt, eta etc are applied yet) - // Skip non-hard-scatter - if (trk_sim_event[isimtrk] != 0) - continue; + if (ana.jet_branches) { + auto const& trk_sim_deltaEta = trk.getVF("sim_deltaEta"); + auto const& trk_sim_deltaPhi = trk.getVF("sim_deltaPhi"); + auto const& trk_sim_deltaR = trk.getVF("sim_deltaR"); + auto const& trk_sim_jet_eta = trk.getVF("sim_jet_eta"); + auto const& trk_sim_jet_phi = trk.getVF("sim_jet_phi"); + auto const& trk_sim_jet_pt = trk.getVF("sim_jet_pt"); ana.tx->pushbackToBranch("sim_deltaEta", trk_sim_deltaEta[isimtrk]); ana.tx->pushbackToBranch("sim_deltaPhi", trk_sim_deltaPhi[isimtrk]); @@ -343,266 +585,1215 @@ void setOutputBranches(LSTEvent* event) { ana.tx->pushbackToBranch("sim_jet_eta", trk_sim_jet_eta[isimtrk]); ana.tx->pushbackToBranch("sim_jet_phi", trk_sim_jet_phi[isimtrk]); ana.tx->pushbackToBranch("sim_jet_pt", trk_sim_jet_pt[isimtrk]); - - ana.tx->pushbackToBranch("sim_pt", trk_sim_pt[isimtrk]); - ana.tx->pushbackToBranch("sim_eta", trk_sim_eta[isimtrk]); - ana.tx->pushbackToBranch("sim_phi", trk_sim_phi[isimtrk]); - ana.tx->pushbackToBranch("sim_pca_dxy", trk_sim_pca_dxy[isimtrk]); - ana.tx->pushbackToBranch("sim_pca_dz", trk_sim_pca_dz[isimtrk]); - ana.tx->pushbackToBranch("sim_q", trk_sim_q[isimtrk]); - ana.tx->pushbackToBranch("sim_event", trk_sim_event[isimtrk]); - ana.tx->pushbackToBranch("sim_pdgId", trk_sim_pdgId[isimtrk]); - - // For vertex we need to look it up from simvtx info - int vtxidx = trk_sim_parentVtxIdx[isimtrk]; - ana.tx->pushbackToBranch("sim_vx", trk_simvtx_x[vtxidx]); - ana.tx->pushbackToBranch("sim_vy", trk_simvtx_y[vtxidx]); - ana.tx->pushbackToBranch("sim_vz", trk_simvtx_z[vtxidx]); - - // The trkNtupIdx is the idx in the trackingNtuple - ana.tx->pushbackToBranch("sim_trkNtupIdx", isimtrk); - - // Increase the counter for accepted simtrk - n_accepted_simtrk++; } - } else { - for (unsigned int isimtrk = 0; isimtrk < trk_sim_pt.size(); ++isimtrk) { - // Skip out-of-time pileup - if (trk_sim_bunchCrossing[isimtrk] != 0) - continue; - - // Skip non-hard-scatter - if (trk_sim_event[isimtrk] != 0) - continue; - - ana.tx->pushbackToBranch("sim_pt", trk_sim_pt[isimtrk]); - ana.tx->pushbackToBranch("sim_eta", trk_sim_eta[isimtrk]); - ana.tx->pushbackToBranch("sim_phi", trk_sim_phi[isimtrk]); - ana.tx->pushbackToBranch("sim_pca_dxy", trk_sim_pca_dxy[isimtrk]); - ana.tx->pushbackToBranch("sim_pca_dz", trk_sim_pca_dz[isimtrk]); - ana.tx->pushbackToBranch("sim_q", trk_sim_q[isimtrk]); - ana.tx->pushbackToBranch("sim_event", trk_sim_event[isimtrk]); - ana.tx->pushbackToBranch("sim_pdgId", trk_sim_pdgId[isimtrk]); - - // For vertex we need to look it up from simvtx info - int vtxidx = trk_sim_parentVtxIdx[isimtrk]; - ana.tx->pushbackToBranch("sim_vx", trk_simvtx_x[vtxidx]); - ana.tx->pushbackToBranch("sim_vy", trk_simvtx_y[vtxidx]); - ana.tx->pushbackToBranch("sim_vz", trk_simvtx_z[vtxidx]); - - // The trkNtupIdx is the idx in the trackingNtuple - ana.tx->pushbackToBranch("sim_trkNtupIdx", isimtrk); - - // Increase the counter for accepted simtrk - n_accepted_simtrk++; - } - } - // Intermediate variables to keep track of matched track candidates for a given sim track - std::vector sim_TC_matched(n_accepted_simtrk); - std::vector sim_TC_matched_mask(n_accepted_simtrk); - std::vector sim_TC_matched_for_duplicate(trk_sim_pt.size()); + // Fill the branch with simulated tracks. + // N.B. these simulated tracks are looser than MTV denominator + ana.tx->pushbackToBranch("sim_pt", trk_sim_pt[isimtrk]); + ana.tx->pushbackToBranch("sim_eta", trk_sim_eta[isimtrk]); + ana.tx->pushbackToBranch("sim_phi", trk_sim_phi[isimtrk]); + ana.tx->pushbackToBranch("sim_pca_dxy", trk_sim_pca_dxy[isimtrk]); + ana.tx->pushbackToBranch("sim_pca_dz", trk_sim_pca_dz[isimtrk]); + ana.tx->pushbackToBranch("sim_q", trk_sim_q[isimtrk]); + ana.tx->pushbackToBranch("sim_pdgId", trk_sim_pdgId[isimtrk]); + + // For vertex we need to look it up from simvtx info for the given simtrack + // for each simulated track, there is an index that points to the production vertex + int vtxidx = trk_sim_parentVtxIdx[isimtrk]; + ana.tx->pushbackToBranch("sim_vx", trk_simvtx_x[vtxidx]); // using the index we retrieve xyz position + ana.tx->pushbackToBranch("sim_vy", trk_simvtx_y[vtxidx]); + ana.tx->pushbackToBranch("sim_vz", trk_simvtx_z[vtxidx]); + ana.tx->pushbackToBranch( + "sim_vtxperp", sqrt(trk_simvtx_x[vtxidx] * trk_simvtx_x[vtxidx] + trk_simvtx_y[vtxidx] * trk_simvtx_y[vtxidx])); + + // The trkNtupIdx is the idx in the trackingNtuple + ana.tx->pushbackToBranch("sim_trkNtupIdx", isimtrk); + + if (ana.extra_sim_branches) { + // Retrieve some track parameter information so we can build a helix + float pt = trk_sim_pt[isimtrk]; + float eta = trk_sim_eta[isimtrk]; + float phi = trk_sim_phi[isimtrk]; + float vx = trk_simvtx_x[vtxidx]; + float vy = trk_simvtx_y[vtxidx]; + float vz = trk_simvtx_z[vtxidx]; + float charge = trk_sim_q[isimtrk]; + + // Build the helix model. This model is useful to compute some specific expected hits. + lst_math::Helix helix(pt, eta, phi, vx, vy, vz, charge); + + // Information to keep track of so we can save to output + std::vector simHitLayer; + std::vector simHitDistxyHelix; + std::vector simHitX; + std::vector simHitY; + std::vector simHitZ; + std::vector simHitDetId; + std::vector recoHitX; + std::vector recoHitY; + std::vector recoHitZ; + std::vector recoHitDetId; + std::vector simHitLayerMinDistxyHelix(11, 999); + + std::vector> simHitIdxs(11); + float k2Rinv1GeVf = (2.99792458e-3 * 3.8) / 2; + + // Loop over the simhits (truth hits) + for (size_t isimhit = 0; isimhit < trk_sim_simHitIdx[isimtrk].size(); ++isimhit) { + // Retrieve the actual index to the simhit_* container of the tracking ntuple + int isimhitidx = trk_sim_simHitIdx[isimtrk][isimhit]; + + // Following computes the distance of the simhit's actual positionin xy to the "expected" xy position based on simhit's z position. + // i.e. Take simhit's z position -> plug them into helix parametric function to obtain the xy position for that given z. + // Then compare the computed xy position from the helix to the simhit's actualy xy position. + // This is a measure of "how off from the original trajectory the simhits are?" + // For example, if the particle got deflected early on due to material, then the xy position distance would be large. + float distxyconsistent = + distxySimHitConsistentWithHelix(helix, isimhitidx, trk_simhit_x, trk_simhit_y, trk_simhit_z); + + // Also retrieve some basic information about the simhit's location (layers, isbarrel?, etc.) + // subdet == 4 means endcap of the outer tracker, subdet == 5 means barrel of the outer tracker) + int subdet = trk_simhit_subdet[isimhitidx]; + int is_endcap = subdet == 4; + + // Now compute "logical layer" index + // N.B. if a hit is in the inner tracker, layer would be staying at layer = 0 + int layer = 0; + if (subdet == 4 or subdet == 5) // this is not an outer tracker hit + // this accounting makes it so that you have layer 1 2 3 4 5 6 in the barrel, and 7 8 9 10 11 in the endcap. (becuase endcap is ph2_subdet == 4) + layer = trk_simhit_layer[isimhitidx] + 6 * (is_endcap); + + // keep track of isimhits in each layers so we can compute mindistxy from previous hit in previous layer + if (subdet == 4 or subdet == 5) + simHitIdxs[layer - 1].push_back(isimhitidx); + + // For this hit, now we push back to the vector that we are keeping track of + simHitLayer.push_back(layer); + simHitDistxyHelix.push_back(distxyconsistent); + simHitX.push_back(trk_simhit_x[isimhitidx]); + simHitY.push_back(trk_simhit_y[isimhitidx]); + simHitZ.push_back(trk_simhit_z[isimhitidx]); + simHitDetId.push_back(trk_simhit_detId[isimhitidx]); + + // Also retrieve all the reco-hits matched to this simhit and also aggregate them + for (size_t irecohit = 0; irecohit < trk_simhit_hitIdx[isimhitidx].size(); ++irecohit) { + recoHitX.push_back(trk_ph2_x[trk_simhit_hitIdx[isimhitidx][irecohit]]); + recoHitY.push_back(trk_ph2_y[trk_simhit_hitIdx[isimhitidx][irecohit]]); + recoHitZ.push_back(trk_ph2_z[trk_simhit_hitIdx[isimhitidx][irecohit]]); + recoHitDetId.push_back(trk_ph2_detId[trk_simhit_hitIdx[isimhitidx][irecohit]]); + } - // Intermediate variables to keep track of matched sim tracks for a given track candidate - std::vector> tc_matched_simIdx; + // If the given simhit that we are dealing with is not in the outer tracker (i.e. layer == 0. see few lines above.) + // then, skip this simhit and go to the next hit. + if (layer == 0) + continue; - // ============ Track candidates ============= - auto const& trackCandidatesBase = event->getTrackCandidatesBase(); - unsigned int nTrackCandidates = trackCandidatesBase.nTrackCandidates(); - for (unsigned int idx = 0; idx < nTrackCandidates; idx++) { - // Compute reco quantities of track candidate based on final object - int type, isFake; - float pt, eta, phi; - std::vector simidx; - std::tie(type, pt, eta, phi, isFake, simidx) = parseTrackCandidate( - event, idx, trk_ph2_x, trk_ph2_y, trk_ph2_z, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); - ana.tx->pushbackToBranch("tc_pt", pt); - ana.tx->pushbackToBranch("tc_eta", eta); - ana.tx->pushbackToBranch("tc_phi", phi); - ana.tx->pushbackToBranch("tc_type", type); - ana.tx->pushbackToBranch("tc_isFake", isFake); - tc_matched_simIdx.push_back(simidx); - - // Loop over matched sim idx and increase counter of TC_matched - for (auto& idx : simidx) { - // NOTE Important to note that the idx of the std::vector<> is same - // as the tracking-ntuple's sim track idx ONLY because event==0 and bunchCrossing==0 condition is applied!! - // Also do not try to access beyond the event and bunchCrossing - if (idx < n_accepted_simtrk) { - sim_TC_matched.at(idx) += 1; - sim_TC_matched_mask.at(idx) |= (1 << type); + // If it is a outer tracker hit, then we keep track of out of the 11 layers, what is the minimum "DistxyHelix" (distance to the expected point in the helix in xy) + // This variable will have a fixed 11 float numbers, and using this to restrict "at least one hit that is not too far from the expected helix" can be useful to select some interesting denominator tracks. + if (distxyconsistent < simHitLayerMinDistxyHelix[layer - 1]) { + simHitLayerMinDistxyHelix[layer - 1] = distxyconsistent; + } } - sim_TC_matched_for_duplicate.at(idx) += 1; - } - } - // Using the intermedaite variables to compute whether a given track candidate is a duplicate - std::vector tc_isDuplicate(tc_matched_simIdx.size()); - // Loop over the track candidates - for (unsigned int i = 0; i < tc_matched_simIdx.size(); ++i) { - bool isDuplicate = false; - // Loop over the sim idx matched to this track candidate - for (unsigned int isim = 0; isim < tc_matched_simIdx[i].size(); ++isim) { - // Using the sim_TC_matched to see whether this track candidate is matched to a sim track that is matched to more than one - int simidx = tc_matched_simIdx[i][isim]; - if (sim_TC_matched_for_duplicate[simidx] > 1) { - isDuplicate = true; - } + std::vector simHitLayerMinDistxyHelixPrevHit(11, 999); + std::vector simHitLayeriSimHitMinDixtxyHelixPrevHit(11, -999); + // // The algorithm will be to start with the main helix from the sim information and get the isimhit with least distxy. + // // Then, from that you find the min distxy and repeat + // for (int ilogicallayer = 0; ilogicallayer < 11; ++ilogicallayer) + // { + // int ilayer = ilogicallayer - 1; + // float prev_pt, prev_eta, prev_phi, prev_vx, prev_vy, prev_vz; + // if (ilayer == 0) + // { + // prev_pt = pt; + // prev_eta = eta; + // prev_phi = phi; + // prev_vx = vx; + // prev_vy = vy; + // prev_vz = vz; + // } + // else + // { + // int isimhitidx = simHitLayeriSimHitMinDixtxyHelixPrevHit[ilayer - 1]; + // TVector3 pp(trk.simhit_px()[isimhitidx], trk.simhit_py()[isimhitidx], trk.simhit_pz()[isimhitidx]); + // prev_pt = pp.Pt(); + // prev_eta = pp.Eta(); + // prev_phi = pp.Phi(); + // prev_vx = trk.simhit_x()[isimhitidx]; + // prev_vy = trk.simhit_y()[isimhitidx]; + // prev_vz = trk.simhit_z()[isimhitidx]; + // } + // SDLMath::Helix prev_helix(prev_pt, prev_eta, prev_phi, prev_vx, prev_vy, prev_vz, charge); + // for (int isimhit = 0; isimhit < simHitIdxs[ilayer].size(); ++isimhit) + // { + // int isimhitidx = simHitIdxs[ilayer][isimhit]; + // float distxyconsistent = distxySimHitConsistentWithHelix(prev_helix, isimhitidx); + // if (simHitLayerMinDistxyHelixPrevHit[ilayer] > distxyconsistent) + // { + // simHitLayerMinDistxyHelixPrevHit[ilayer] = distxyconsistent; + // simHitLayeriSimHitMinDixtxyHelixPrevHit[ilayer] = isimhitidx; + // } + // } + // } + + // Now we fill the branch + ana.tx->pushbackToBranch>("sim_simHitLayer", simHitLayer); + ana.tx->pushbackToBranch>("sim_simHitDistxyHelix", simHitDistxyHelix); + ana.tx->pushbackToBranch>("sim_simHitLayerMinDistxyHelix", simHitLayerMinDistxyHelix); + ana.tx->pushbackToBranch>("sim_simHitLayerMinDistxyPrevHit", simHitLayerMinDistxyHelixPrevHit); + ana.tx->pushbackToBranch>("sim_simHitX", simHitX); + ana.tx->pushbackToBranch>("sim_simHitY", simHitY); + ana.tx->pushbackToBranch>("sim_simHitZ", simHitZ); + ana.tx->pushbackToBranch>("sim_simHitDetId", simHitDetId); + ana.tx->pushbackToBranch>("sim_recoHitX", recoHitX); + ana.tx->pushbackToBranch>("sim_recoHitY", recoHitY); + ana.tx->pushbackToBranch>("sim_recoHitZ", recoHitZ); + ana.tx->pushbackToBranch>("sim_recoHitDetId", recoHitDetId); } - tc_isDuplicate[i] = isDuplicate; + + // Increase the counter for accepted simtrk + n_accepted_simtrk++; } - // Now set the last remaining branches - ana.tx->setBranch>("sim_TC_matched", sim_TC_matched); - ana.tx->setBranch>("sim_TC_matched_mask", sim_TC_matched_mask); - ana.tx->setBranch>>("tc_matched_simIdx", tc_matched_simIdx); - ana.tx->setBranch>("tc_isDuplicate", tc_isDuplicate); + return n_accepted_simtrk; } //________________________________________________________________________________________________________________________________ -void setOptionalOutputBranches(LSTEvent* event) { -#ifdef CUT_VALUE_DEBUG +std::map setMiniDoubletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac) { + //-------------------------------------------- + // + // + // Mini-Doublets + // + // + //-------------------------------------------- + + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_ph2_subdet = trk.getVUS("ph2_subdet"); + auto const& trk_ph2_layer = trk.getVUS("ph2_layer"); + auto const& trk_ph2_detId = trk.getVU("ph2_detId"); + auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); + auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); + auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + + auto const& hitsBase = event->getInput(); + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& miniDoublets = event->getMiniDoublets(); + auto const& miniDoubletsOccupancy = event->getMiniDoublets(); + + // Following are some vectors to keep track of the information to write to the ntuple + // N.B. following two branches have a length for the entire sim track, but what actually will be written in sim_mdIdxAll branch is NOT that long + // Later in the code, it will restrict to only the ones to write out. + // The reason at this stage, the entire mdIdxAll is being tracked is to compute duplicate properly later on + // When computing a duplicate object it is important to consider all simulated tracks including pileup tracks + int n_total_simtrk = trk_sim_pt.size(); + std::vector> sim_mdIdxAll(n_total_simtrk); + std::vector> sim_mdIdxAllFrac(n_total_simtrk); + std::vector> md_simIdxAll; + std::vector> md_simIdxAllFrac; + + // global md index that will be used to keep track of md being outputted to the ntuple + // each time a md is written out the following will be counted up + unsigned int md_idx = 0; + + // map to keep track of (GPU mdIdx) -> (md_idx in ntuple output) + // There is a specific mdIdx used to navigate the GPU array of mini-doublets + std::map md_idx_map; + + // First loop over the modules (roughly there are ~13k pair of pt modules) + for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) { + // For each pt module pair, we loop over mini-doublets created + for (unsigned int iMD = 0; iMD < miniDoubletsOccupancy.nMDs()[idx]; iMD++) { + // Compute the specific MD index to access specific spot in the array of GPU memory + unsigned int mdIdx = ranges.miniDoubletModuleIndices()[idx] + iMD; + + // From that gpu memory index "mdIdx" -> output ntuple's md index is mapped + // This is useful later when connecting higher level objects to point to specific one in the ntuple + md_idx_map[mdIdx] = md_idx; - setPixelQuintupletOutputBranches(event); - setQuintupletOutputBranches(event); - setPixelTripletOutputBranches(event); - setOccupancyBranches(event); - setT3DNNBranches(event); - setT5DNNBranches(event); - setpT3DNNBranches(event); - setpLSOutputBranches(event); + // Access the list of hits in the mini-doublets (there are only two in this case) + std::vector hit_idx, hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFromMD(event, mdIdx); + // And then compute matching between simtrack and the mini-doublets + std::vector simidx; + std::vector simidxfrac; + std::tie(simidx, simidxfrac) = + matchedSimTrkIdxsAndFracs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); + + // Obtain the lower and upper hit information to compute some basic property of the mini-doublets + unsigned int LowerHitIndex = miniDoublets.anchorHitIndices()[mdIdx]; + unsigned int UpperHitIndex = miniDoublets.outerHitIndices()[mdIdx]; + unsigned int hit0 = hitsBase.idxs()[LowerHitIndex]; + unsigned int hit1 = hitsBase.idxs()[UpperHitIndex]; + float anchor_x = hitsBase.xs()[LowerHitIndex]; + float anchor_y = hitsBase.ys()[LowerHitIndex]; + float anchor_z = hitsBase.zs()[LowerHitIndex]; + float other_x = hitsBase.xs()[UpperHitIndex]; + float other_y = hitsBase.ys()[UpperHitIndex]; + float other_z = hitsBase.zs()[UpperHitIndex]; + + // Construct the anchor hit 3 vector + lst_math::Hit anchor_hit(anchor_x, anchor_y, anchor_z, LowerHitIndex); + + // Pt is computed via dphichange and the eta and phi are computed based on anchor hit + float dphichange = miniDoublets.dphichanges()[mdIdx]; + float dphi = miniDoublets.dphis()[mdIdx]; + float dz = miniDoublets.dzs()[mdIdx]; + float k2Rinv1GeVf = (2.99792458e-3 * 3.8) / 2; + float pt = anchor_hit.rt() * k2Rinv1GeVf / sin(dphichange); + float eta = anchor_hit.eta(); + float phi = anchor_hit.phi(); + + // Obtain where the actual hit is located in terms of their layer, module, rod, and ring number + int subdet = trk_ph2_subdet[hit0]; + int is_endcap = subdet == 4; + // this accounting makes it so that you have layer 1 2 3 4 5 6 in the barrel, and 7 8 9 10 11 in the endcap. (becuase endcap is ph2_subdet == 4) + int layer = trk_ph2_layer[hit0] + 6 * (is_endcap); + int detId = trk_ph2_detId[hit0]; + // See https://github.com/SegmentLinking/TrackLooper/blob/158804cab7fd0976264a7bc4cee236f4986328c2/SDL/Module.cc and Module.h + int ring = (detId & (15 << 12)) >> 12; + int isPS = is_endcap ? (layer <= 2 ? ring <= 10 : ring <= 7) : layer <= 3; + + // Write out the ntuple + ana.tx->pushbackToBranch("md_pt", pt); + ana.tx->pushbackToBranch("md_eta", eta); + ana.tx->pushbackToBranch("md_phi", phi); +#ifdef CUT_VALUE_DEBUG + ana.tx->pushbackToBranch("md_dphichange", dphichange); + ana.tx->pushbackToBranch("md_dphi", dphi); + ana.tx->pushbackToBranch("md_dz", dz); #endif + ana.tx->pushbackToBranch("md_anchor_x", anchor_x); + ana.tx->pushbackToBranch("md_anchor_y", anchor_y); + ana.tx->pushbackToBranch("md_anchor_z", anchor_z); + ana.tx->pushbackToBranch("md_other_x", other_x); + ana.tx->pushbackToBranch("md_other_y", other_y); + ana.tx->pushbackToBranch("md_other_z", other_z); + ana.tx->pushbackToBranch("md_type", isPS); + ana.tx->pushbackToBranch("md_layer", layer); + ana.tx->pushbackToBranch("md_detId", detId); + + // Compute whether this is a fake + bool isfake = true; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + if (simidxfrac[isim] > matchfrac) { + isfake = false; + break; + } + } + ana.tx->pushbackToBranch("md_isFake", isfake); + + // For this md, keep track of all the simidx that are matched + md_simIdxAll.push_back(simidx); + md_simIdxAllFrac.push_back(simidxfrac); + + // The book keeping of opposite mapping is done here + // For each matched sim idx, we go back and keep track of which obj it is matched to. + // Loop over all the matched sim idx + for (size_t is = 0; is < simidx.size(); ++is) { + // For this matched sim index keep track (sim -> md) mapping + int sim_idx = simidx.at(is); + float sim_idx_frac = simidxfrac.at(is); + sim_mdIdxAll.at(sim_idx).push_back(md_idx); + sim_mdIdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } + + // Also, among the simidx matches, find the best match (highest fractional match) + // N.B. the simidx is already returned sorted by highest number of "nhits" match + // So as it loops over, the condition will ensure that the highest fraction with highest nhits will be matched with the priority given to highest fraction + int md_simIdx = -999; + float md_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > md_simIdxBestFrac and thisfrac > matchfrac) { + md_simIdxBestFrac = thisfrac; + md_simIdx = thisidx; + } + } + + // the best match index will then be saved here + ana.tx->pushbackToBranch("md_simIdx", md_simIdx); + + // Count up the md_idx + md_idx++; + } + } + + // Now save the (obj -> simidx) mapping + ana.tx->setBranch>>("md_simIdxAll", md_simIdxAll); + ana.tx->setBranch>>("md_simIdxAllFrac", md_simIdxAllFrac); + + // Not all (sim->objIdx) will be saved but only for the sim that is from hard scatter and current bunch crossing + // So a restriction up to only "n_accepted_simtrk" done by chopping off the rest + // N.B. the reason we can simply take the first "n_accepted_simtrk" is because the tracking ntuple is organized such that those sim tracks show up on the first "n_accepted_simtrk" of tracks. + std::vector> sim_mdIdxAll_to_write; + std::vector> sim_mdIdxAllFrac_to_write; + std::copy(sim_mdIdxAll.begin(), sim_mdIdxAll.begin() + n_accepted_simtrk, std::back_inserter(sim_mdIdxAll_to_write)); + std::copy(sim_mdIdxAllFrac.begin(), + sim_mdIdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_mdIdxAllFrac_to_write)); + ana.tx->setBranch>>("sim_mdIdxAll", sim_mdIdxAll_to_write); + ana.tx->setBranch>>("sim_mdIdxAllFrac", sim_mdIdxAllFrac_to_write); + + return md_idx_map; } //________________________________________________________________________________________________________________________________ -void setOccupancyBranches(LSTEvent* event) { - auto modules = event->getModules(); - auto miniDoublets = event->getMiniDoublets(); - auto segments = event->getSegments(); - auto triplets = event->getTriplets(); - auto quintuplets = event->getQuintuplets(); - auto pixelQuintuplets = event->getPixelQuintuplets(); - auto pixelTriplets = event->getPixelTriplets(); - auto trackCandidatesBase = event->getTrackCandidatesBase(); - - std::vector moduleLayer; - std::vector moduleSubdet; - std::vector moduleRing; - std::vector moduleRod; - std::vector moduleModule; - std::vector moduleEta; - std::vector moduleR; - std::vector moduleIsTilted; - std::vector trackCandidateOccupancy; - std::vector tripletOccupancy; - std::vector segmentOccupancy; - std::vector mdOccupancy; - std::vector quintupletOccupancy; +std::map setLineSegmentBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& md_idx_map) { + //-------------------------------------------- + // + // + // Line Segments + // + // + //-------------------------------------------- - for (unsigned int lowerIdx = 0; lowerIdx <= modules.nLowerModules(); lowerIdx++) { - //layer = 0, subdet = 0 => pixel module - moduleLayer.push_back(modules.layers()[lowerIdx]); - moduleSubdet.push_back(modules.subdets()[lowerIdx]); - moduleRing.push_back(modules.rings()[lowerIdx]); - moduleRod.push_back(modules.rods()[lowerIdx]); - moduleEta.push_back(modules.eta()[lowerIdx]); - moduleR.push_back(modules.r()[lowerIdx]); - bool isTilted = (modules.subdets()[lowerIdx] == 5 and modules.sides()[lowerIdx] != 3); - moduleIsTilted.push_back(isTilted); - moduleModule.push_back(modules.modules()[lowerIdx]); - segmentOccupancy.push_back(segments.totOccupancySegments()[lowerIdx]); - mdOccupancy.push_back(miniDoublets.totOccupancyMDs()[lowerIdx]); + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); + auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); + auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); - if (lowerIdx < modules.nLowerModules()) { - quintupletOccupancy.push_back(quintuplets.totOccupancyQuintuplets()[lowerIdx]); - tripletOccupancy.push_back(triplets.totOccupancyTriplets()[lowerIdx]); + auto const& hitsBase = event->getInput(); + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& segments = event->getSegments(); + auto const& segmentsOccupancy = event->getSegments(); + + // Following are some vectors to keep track of the information to write to the ntuple + // N.B. following two branches have a length for the entire sim track, but what actually will be written in sim_objIdxAll branch is NOT that long + // Later in the code, it will restrict to only the ones to write out. + // The reason at this stage, the entire objIdxAll is being tracked is to compute duplicate properly later on + // When computing a duplicate object it is important to consider all simulated tracks including pileup tracks + int n_total_simtrk = trk_sim_pt.size(); + std::vector> sim_lsIdxAll(n_total_simtrk); + std::vector> sim_lsIdxAllFrac(n_total_simtrk); + std::vector> ls_simIdxAll; + std::vector> ls_simIdxAllFrac; + + // global index that will be used to keep track of obj being outputted to the ntuple + // each time a obj is written out the following will be counted up + unsigned int ls_idx = 0; + + // map to keep track of (GPU objIdx) -> (obj_idx in ntuple output) + // There is a specific objIdx used to navigate the GPU array of mini-doublets + std::map ls_idx_map; + + // First loop over the modules (roughly there are ~13k pair of pt modules) + for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) { + // For each pt module pair, we loop over objects created + for (unsigned int iLS = 0; iLS < segmentsOccupancy.nSegments()[idx]; iLS++) { + // Compute the specific obj index to access specific spot in the array of GPU memory + unsigned int lsIdx = ranges.segmentModuleIndices()[idx] + iLS; + + // From that gpu memory index "objIdx" -> output ntuple's obj index is mapped + // This is useful later when connecting higher level objects to point to specific one in the ntuple + ls_idx_map[lsIdx] = ls_idx; + + // Access the list of hits in the objects (there are only two in this case) + std::vector hit_idx, hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFromLS(event, lsIdx); + + // And then compute matching between simtrack and the objects + std::vector simidx; + std::vector simidxfrac; + std::tie(simidx, simidxfrac) = + matchedSimTrkIdxsAndFracs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); + std::vector mdIdxs = getMDsFromLS(event, lsIdx); + + // Computing line segment pt estimate (assuming beam spot is at zero) + lst_math::Hit hitA(0, 0, 0); + lst_math::Hit hitB(hitsBase.xs()[hit_idx[0]], hitsBase.ys()[hit_idx[0]], hitsBase.zs()[hit_idx[0]]); + lst_math::Hit hitC(hitsBase.xs()[hit_idx[2]], hitsBase.ys()[hit_idx[2]], hitsBase.zs()[hit_idx[2]]); + lst_math::Hit center = lst_math::getCenterFromThreePoints(hitA, hitB, hitC); + float pt = lst_math::ptEstimateFromRadius(center.rt()); + float eta = hitC.eta(); + float phi = hitB.phi(); + +#ifdef CUT_VALUE_DEBUG + float zHi = segments.zHis()[lsIdx]; + float zLo = segments.zLos()[lsIdx]; + float rtHi = segments.rtHis()[lsIdx]; + float rtLo = segments.rtLos()[lsIdx]; + float dAlphaInner = segments.dAlphaInners()[lsIdx]; + float dAlphaOuter = segments.dAlphaOuters()[lsIdx]; + float dAlphaInnerOuter = segments.dAlphaInnerOuters()[lsIdx]; + float dPhi = segments.dPhis()[lsIdx]; + float dPhiMin = segments.dPhiMins()[lsIdx]; + float dPhiMax = segments.dPhiMaxs()[lsIdx]; + float dPhiChange = segments.dPhiChanges()[lsIdx]; + float dPhiChangeMin = segments.dPhiChangeMins()[lsIdx]; + float dPhiChangeMax = segments.dPhiChangeMaxs()[lsIdx]; +#endif + + // Write out the ntuple + ana.tx->pushbackToBranch("ls_pt", pt); + ana.tx->pushbackToBranch("ls_eta", eta); + ana.tx->pushbackToBranch("ls_phi", phi); +#ifdef CUT_VALUE_DEBUG + ana.tx->pushbackToBranch("ls_zHis", zHi); + ana.tx->pushbackToBranch("ls_zLos", zLo); + ana.tx->pushbackToBranch("ls_rtHis", rtHi); + ana.tx->pushbackToBranch("ls_rtLos", rtLo); + ana.tx->pushbackToBranch("ls_dPhis", dPhi); + ana.tx->pushbackToBranch("ls_dPhiMins", dPhiMin); + ana.tx->pushbackToBranch("ls_dPhiMaxs", dPhiMax); + ana.tx->pushbackToBranch("ls_dPhiChanges", dPhiChange); + ana.tx->pushbackToBranch("ls_dPhiChangeMins", dPhiChangeMin); + ana.tx->pushbackToBranch("ls_dPhiChangeMaxs", dPhiChangeMax); + ana.tx->pushbackToBranch("ls_dAlphaInners", dAlphaInner); + ana.tx->pushbackToBranch("ls_dAlphaOuters", dAlphaOuter); + ana.tx->pushbackToBranch("ls_dAlphaInnerOuters", dAlphaInnerOuter); + +#endif + if (ana.md_branches) { + ana.tx->pushbackToBranch("ls_mdIdx0", md_idx_map.at(mdIdxs[0])); + ana.tx->pushbackToBranch("ls_mdIdx1", md_idx_map.at(mdIdxs[1])); + } + + // Compute whether this is a fake + bool isfake = true; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + if (simidxfrac[isim] > matchfrac) { + isfake = false; + break; + } + } + ana.tx->pushbackToBranch("ls_isFake", isfake); + + // For this obj, keep track of all the simidx that are matched + ls_simIdxAll.push_back(simidx); + ls_simIdxAllFrac.push_back(simidxfrac); + + // The book keeping of opposite mapping is done here + // For each matched sim idx, we go back and keep track of which obj it is matched to. + // Loop over all the matched sim idx + for (size_t is = 0; is < simidx.size(); ++is) { + int sim_idx = simidx.at(is); + float sim_idx_frac = simidxfrac.at(is); + if (sim_idx < n_total_simtrk) { + sim_lsIdxAll.at(sim_idx).push_back(ls_idx); + sim_lsIdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } + } + + // Also, among the simidx matches, find the best match (highest fractional match) + // N.B. the simidx is already returned sorted by highest number of "nhits" match + // So as it loops over, the condition will ensure that the highest fraction with highest nhits will be matched with the priority given to highest fraction + int ls_simIdx = -999; + float ls_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > ls_simIdxBestFrac and thisfrac > matchfrac) { + ls_simIdxBestFrac = thisfrac; + ls_simIdx = thisidx; + } + } + + // the best match index will then be saved here + ana.tx->pushbackToBranch("ls_simIdx", ls_simIdx); + + // Count up the index + ls_idx++; } } - ana.tx->setBranch>("module_layers", moduleLayer); - ana.tx->setBranch>("module_subdets", moduleSubdet); - ana.tx->setBranch>("module_rings", moduleRing); - ana.tx->setBranch>("module_rods", moduleRod); - ana.tx->setBranch>("module_modules", moduleModule); - ana.tx->setBranch>("module_isTilted", moduleIsTilted); - ana.tx->setBranch>("module_eta", moduleEta); - ana.tx->setBranch>("module_r", moduleR); - ana.tx->setBranch>("md_occupancies", mdOccupancy); - ana.tx->setBranch>("sg_occupancies", segmentOccupancy); - ana.tx->setBranch>("t3_occupancies", tripletOccupancy); - ana.tx->setBranch("tc_occupancies", trackCandidatesBase.nTrackCandidates()); - ana.tx->setBranch("pT3_occupancies", pixelTriplets.totOccupancyPixelTriplets()); - ana.tx->setBranch>("t5_occupancies", quintupletOccupancy); - ana.tx->setBranch("pT5_occupancies", pixelQuintuplets.totOccupancyPixelQuintuplets()); + // Now save the (obj -> simidx) mapping + ana.tx->setBranch>>("ls_simIdxAll", ls_simIdxAll); + ana.tx->setBranch>>("ls_simIdxAllFrac", ls_simIdxAllFrac); + + // Not all (sim->objIdx) will be saved but only for the sim that is from hard scatter and current bunch crossing + // So a restriction up to only "n_accepted_simtrk" done by chopping off the rest + // N.B. the reason we can simply take the first "n_accepted_simtrk" is because the tracking ntuple is organized such that those sim tracks show up on the first "n_accepted_simtrk" of tracks. + std::vector> sim_lsIdxAll_to_write; + std::vector> sim_lsIdxAllFrac_to_write; + std::copy(sim_lsIdxAll.begin(), sim_lsIdxAll.begin() + n_accepted_simtrk, std::back_inserter(sim_lsIdxAll_to_write)); + std::copy(sim_lsIdxAllFrac.begin(), + sim_lsIdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_lsIdxAllFrac_to_write)); + ana.tx->setBranch>>("sim_lsIdxAll", sim_lsIdxAll_to_write); + ana.tx->setBranch>>("sim_lsIdxAllFrac", sim_lsIdxAllFrac_to_write); + + return ls_idx_map; } //________________________________________________________________________________________________________________________________ -void setPixelQuintupletOutputBranches(LSTEvent* event) { - // ============ pT5 ============= - auto const pixelQuintuplets = event->getPixelQuintuplets(); - auto const quintuplets = event->getQuintuplets(); - auto const pixelSeeds = event->getInput(); - auto modules = event->getModules(); - int n_accepted_simtrk = ana.tx->getBranch>("sim_TC_matched").size(); +std::map setTripletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& ls_idx_map) { + //-------------------------------------------- + // + // + // Triplet + // + // + //-------------------------------------------- - unsigned int nPixelQuintuplets = pixelQuintuplets.nPixelQuintuplets(); - std::vector sim_pT5_matched(n_accepted_simtrk); - std::vector> pT5_matched_simIdx; + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); + auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); + auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + + auto const& hitsBase = event->getInput(); + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& triplets = event->getTriplets(); + auto const& tripletOccupancies = event->getTriplets(); + + int n_total_simtrk = trk_sim_pt.size(); + std::vector sim_t3_matched(n_accepted_simtrk); + std::vector> sim_t3IdxAll(n_total_simtrk); + std::vector> sim_t3IdxAllFrac(n_total_simtrk); + std::vector> t3_simIdxAll; + std::vector> t3_simIdxAllFrac; + // Then obtain the lower module index + unsigned int t3_idx = 0; // global t3 index that will be used to keep track of t3 being outputted to the ntuple + // map to keep track of (GPU t3Idx) -> (t3_idx in ntuple output) + std::map t3_idx_map; + // printT3s(event); + for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) { + unsigned int nmods = modules.nLowerModules(); + for (unsigned int iT3 = 0; iT3 < tripletOccupancies.nTriplets()[idx]; iT3++) { + unsigned int t3Idx = ranges.tripletModuleIndices()[idx] + iT3; + t3_idx_map[t3Idx] = t3_idx; + std::vector hit_idx, hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFromT3(event, t3Idx); + std::vector simidx; + std::vector simidxfrac; + std::tie(simidx, simidxfrac) = + matchedSimTrkIdxsAndFracs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); + std::vector lsIdxs = getLSsFromT3(event, t3Idx); + if (ana.ls_branches) { + ana.tx->pushbackToBranch("t3_lsIdx0", ls_idx_map.at(lsIdxs[0])); + ana.tx->pushbackToBranch("t3_lsIdx1", ls_idx_map.at(lsIdxs[1])); + } + // Computing line segment pt estimate (assuming beam spot is at zero) + lst_math::Hit hitA(hitsBase.xs()[hit_idx[0]], hitsBase.ys()[hit_idx[0]], hitsBase.zs()[hit_idx[0]]); + lst_math::Hit hitC(hitsBase.xs()[hit_idx[4]], hitsBase.ys()[hit_idx[4]], hitsBase.zs()[hit_idx[4]]); + float pt = __H2F(triplets.radius()[t3Idx]) * k2Rinv1GeVf * 2; + float eta = hitC.eta(); + float phi = hitA.phi(); + ana.tx->pushbackToBranch("t3_pt", pt); + ana.tx->pushbackToBranch("t3_eta", eta); + ana.tx->pushbackToBranch("t3_phi", phi); + bool isfake = true; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + if (simidxfrac[isim] > matchfrac) { + isfake = false; + break; + } + } + ana.tx->pushbackToBranch("t3_isFake", isfake); + t3_simIdxAll.push_back(simidx); + t3_simIdxAllFrac.push_back(simidxfrac); + for (size_t is = 0; is < simidx.size(); ++is) { + int sim_idx = simidx.at(is); + if (sim_idx < n_accepted_simtrk) { + sim_t3_matched.at(sim_idx) += 1; + } + float sim_idx_frac = simidxfrac.at(is); + if (sim_idx < n_total_simtrk) { + sim_t3IdxAll.at(sim_idx).push_back(t3_idx); + sim_t3IdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } + } + int t3_simIdx = -999; + float t3_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > t3_simIdxBestFrac and thisfrac > matchfrac) { + t3_simIdxBestFrac = thisfrac; + t3_simIdx = thisidx; + } + } + ana.tx->pushbackToBranch("t3_simIdx", t3_simIdx); + // count global + t3_idx++; + } + } + ana.tx->setBranch>>("t3_simIdxAll", t3_simIdxAll); + ana.tx->setBranch>>("t3_simIdxAllFrac", t3_simIdxAllFrac); + std::vector> sim_t3IdxAll_to_write; + std::vector> sim_t3IdxAllFrac_to_write; + std::copy(sim_t3IdxAll.begin(), sim_t3IdxAll.begin() + n_accepted_simtrk, std::back_inserter(sim_t3IdxAll_to_write)); + std::copy(sim_t3IdxAllFrac.begin(), + sim_t3IdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_t3IdxAllFrac_to_write)); + ana.tx->setBranch>>("sim_t3IdxAll", sim_t3IdxAll_to_write); + ana.tx->setBranch>>("sim_t3IdxAllFrac", sim_t3IdxAllFrac_to_write); + + // Using the intermedaite variables to compute whether a given object is a duplicate + std::vector t3_isDuplicate(t3_simIdxAll.size()); + for (unsigned int i = 0; i < t3_simIdxAll.size(); i++) { + bool isDuplicate = true; + for (unsigned int isim = 0; isim < t3_simIdxAll[i].size(); isim++) { + int simidx = t3_simIdxAll[i][isim]; + if (simidx < n_accepted_simtrk) { + if (sim_t3_matched[simidx] > 1) { + isDuplicate = true; + } + } + } + t3_isDuplicate[i] = isDuplicate; + } + ana.tx->setBranch>("t3_isDuplicate", t3_isDuplicate); + + return t3_idx_map; +} + +//________________________________________________________________________________________________________________________________ +std::map setQuintupletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& t3_idx_map) { + //-------------------------------------------- + // + // + // Quintuplet + // + // + //-------------------------------------------- + + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx"); + auto const& trk_simvtx_x = trk.getVF("simvtx_x"); + auto const& trk_simvtx_y = trk.getVF("simvtx_y"); + auto const& trk_simvtx_z = trk.getVF("simvtx_z"); + auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); + auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); + auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + + auto const& hitsBase = event->getInput(); + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& quintuplets = event->getQuintuplets(); + auto const& quintupletOccupancies = event->getQuintuplets(); + + int n_total_simtrk = trk_sim_pt.size(); + std::vector sim_t5_matched(n_accepted_simtrk); + std::vector> sim_t5IdxAll(n_total_simtrk); + std::vector> sim_t5IdxAllFrac(n_total_simtrk); + std::vector> t5_simIdxAll; + std::vector> t5_simIdxAllFrac; + // Then obtain the lower module index + unsigned int t5_idx = 0; // global t5 index that will be used to keep track of t5 being outputted to the ntuple + // map to keep track of (GPU t5Idx) -> (t5_idx in ntuple output) + std::map t5_idx_map; + // printT3s(event); + for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) { + unsigned int nmods = modules.nLowerModules(); + for (unsigned int iT5 = 0; iT5 < quintupletOccupancies.nQuintuplets()[idx]; iT5++) { + unsigned int t5Idx = ranges.quintupletModuleIndices()[idx] + iT5; + t5_idx_map[t5Idx] = t5_idx; + std::vector hit_idx, hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFromT5(event, t5Idx); + std::vector simidx; + std::vector simidxfrac; + float percent_matched; + std::tie(simidx, simidxfrac) = matchedSimTrkIdxsAndFracs(hit_idx, + hit_type, + trk_simhit_simTrkIdx, + trk_ph2_simHitIdx, + trk_pix_simHitIdx, + false, + matchfrac, + &percent_matched); + std::vector t3Idxs = getT3sFromT5(event, t5Idx); + if (ana.t3_branches) { + ana.tx->pushbackToBranch("t5_t3Idx0", t3_idx_map.at(t3Idxs[0])); + ana.tx->pushbackToBranch("t5_t3Idx1", t3_idx_map.at(t3Idxs[1])); + } + float pt = __H2F(quintuplets.innerRadius()[t5Idx]) * k2Rinv1GeVf * 2; + float eta = __H2F(quintuplets.eta()[t5Idx]); + float phi = __H2F(quintuplets.phi()[t5Idx]); + ana.tx->pushbackToBranch("t5_pt", pt); + ana.tx->pushbackToBranch("t5_eta", eta); + ana.tx->pushbackToBranch("t5_phi", phi); + ana.tx->pushbackToBranch("t5_innerRadius", __H2F(quintuplets.innerRadius()[t5Idx])); + ana.tx->pushbackToBranch("t5_bridgeRadius", __H2F(quintuplets.bridgeRadius()[t5Idx])); + ana.tx->pushbackToBranch("t5_outerRadius", __H2F(quintuplets.outerRadius()[t5Idx])); + ana.tx->pushbackToBranch("t5_pMatched", percent_matched); + bool isfake = true; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + if (simidxfrac[isim] > matchfrac) { + isfake = false; + break; + } + } + ana.tx->pushbackToBranch("t5_isFake", isfake); + t5_simIdxAll.push_back(simidx); + t5_simIdxAllFrac.push_back(simidxfrac); + for (size_t is = 0; is < simidx.size(); ++is) { + int sim_idx = simidx.at(is); + if (sim_idx < n_accepted_simtrk) { + sim_t5_matched.at(sim_idx) += 1; + } + float sim_idx_frac = simidxfrac.at(is); + if (sim_idx < n_total_simtrk) { + sim_t5IdxAll.at(sim_idx).push_back(t5_idx); + sim_t5IdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } + } + int t5_simIdx = -999; + float t5_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > t5_simIdxBestFrac and thisfrac > matchfrac) { + t5_simIdxBestFrac = thisfrac; + t5_simIdx = thisidx; + } + } + ana.tx->pushbackToBranch("t5_simIdx", t5_simIdx); + // count global + t5_idx++; + + // Avoid fakes when calculating the vertex distance, set default to 0.0. + if (simidx.size() == 0) { + ana.tx->pushbackToBranch("t5_sim_vxy", 0.0); + ana.tx->pushbackToBranch("t5_sim_vz", 0.0); + } else { + int vtxidx = trk_sim_parentVtxIdx[simidx[0]]; + float vtx_x = trk_simvtx_x[vtxidx]; + float vtx_y = trk_simvtx_y[vtxidx]; + float vtx_z = trk_simvtx_z[vtxidx]; + + ana.tx->pushbackToBranch("t5_sim_vxy", sqrt(vtx_x * vtx_x + vtx_y * vtx_y)); + ana.tx->pushbackToBranch("t5_sim_vz", vtx_z); + } + } + } + ana.tx->setBranch>>("t5_simIdxAll", t5_simIdxAll); + ana.tx->setBranch>>("t5_simIdxAllFrac", t5_simIdxAllFrac); + std::vector> sim_t5IdxAll_to_write; + std::vector> sim_t5IdxAllFrac_to_write; + std::copy(sim_t5IdxAll.begin(), sim_t5IdxAll.begin() + n_accepted_simtrk, std::back_inserter(sim_t5IdxAll_to_write)); + std::copy(sim_t5IdxAllFrac.begin(), + sim_t5IdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_t5IdxAllFrac_to_write)); + ana.tx->setBranch>>("sim_t5IdxAll", sim_t5IdxAll_to_write); + ana.tx->setBranch>>("sim_t5IdxAllFrac", sim_t5IdxAllFrac_to_write); + + std::vector t5_isDuplicate(t5_simIdxAll.size()); + for (unsigned int i = 0; i < t5_simIdxAll.size(); i++) { + bool isDuplicate = false; + for (unsigned int isim = 0; isim < t5_simIdxAll[i].size(); isim++) { + int simidx = t5_simIdxAll[i][isim]; + if (simidx < n_accepted_simtrk) { + if (sim_t5_matched[simidx] > 1) { + isDuplicate = true; + } + } + } + t5_isDuplicate[i] = isDuplicate; + } + ana.tx->setBranch>("t5_isDuplicate", t5_isDuplicate); + + return t5_idx_map; +} +//________________________________________________________________________________________________________________________________ +std::map setPixelLineSegmentBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac) { + //-------------------------------------------- + // + // + // pLS + // + // + //-------------------------------------------- + + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_see_pt = trk.getVF("see_pt"); + auto const& trk_see_eta = trk.getVF("see_eta"); + auto const& trk_see_phi = trk.getVF("see_phi"); + auto const& trk_see_hitIdx = trk.getVVI("see_hitIdx"); + auto const& trk_see_hitType = trk.getVVI("see_hitType"); + auto const& trk_pix_x = trk.getVF("pix_x"); + auto const& trk_pix_y = trk.getVF("pix_y"); + auto const& trk_pix_z = trk.getVF("pix_z"); auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); - for (unsigned int pT5 = 0; pT5 < nPixelQuintuplets; pT5++) { - unsigned int T5Index = getT5FrompT5(event, pT5); - unsigned int pLSIndex = getPixelLSFrompT5(event, pT5); - float pt = (__H2F(quintuplets.innerRadius()[T5Index]) * k2Rinv1GeVf * 2 + pixelSeeds.ptIn()[pLSIndex]) / 2; - float eta = pixelSeeds.eta()[pLSIndex]; - float phi = pixelSeeds.phi()[pLSIndex]; - - std::vector hit_idx = getHitIdxsFrompT5(event, pT5); - std::vector module_idx = getModuleIdxsFrompT5(event, pT5); - std::vector hit_type = getHitTypesFrompT5(event, pT5); - - int layer_binary = 1; - int moduleType_binary = 0; - for (size_t i = 0; i < module_idx.size(); i += 2) { - layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4))); - moduleType_binary |= (modules.moduleType()[module_idx[i]] << i); + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& pixelSeeds = event->getInput(); + auto const& pixelSegments = event->getPixelSegments(); + auto const& segmentsOccupancy = event->getSegments(); + + int n_total_simtrk = trk_sim_pt.size(); + std::vector sim_pLS_matched(n_accepted_simtrk, 0); + std::vector> sim_plsIdxAll(n_total_simtrk); + std::vector> sim_plsIdxAllFrac(n_total_simtrk); + std::vector> pls_simIdxAll; + std::vector> pls_simIdxAllFrac; + // Then obtain the lower module index + unsigned int pls_idx = 0; // global pls index that will be used to keep track of pls being outputted to the ntuple + // map to keep track of (GPU plsIdx) -> (pls_idx in ntuple output) + std::map pls_idx_map; + unsigned int n_pls = segmentsOccupancy.nSegments()[modules.nLowerModules()]; + unsigned int pls_range_start = ranges.segmentModuleIndices()[modules.nLowerModules()]; + for (unsigned int ipLS = 0; ipLS < n_pls; ipLS++) { + unsigned int plsIdx = pls_range_start + ipLS; + pls_idx_map[plsIdx] = pls_idx; + std::vector hit_idx, hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFrompLS(event, ipLS); + std::vector simidx; + std::vector simidxfrac; + std::tie(simidx, simidxfrac) = + matchedSimTrkIdxsAndFracs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); + ana.tx->pushbackToBranch("pLS_pt", pixelSeeds.ptIn()[ipLS]); + ana.tx->pushbackToBranch("pLS_ptErr", pixelSeeds.ptErr()[ipLS]); + ana.tx->pushbackToBranch("pLS_eta", pixelSeeds.eta()[ipLS]); + ana.tx->pushbackToBranch("pLS_etaErr", pixelSeeds.etaErr()[ipLS]); + ana.tx->pushbackToBranch("pLS_phi", pixelSeeds.phi()[ipLS]); + ana.tx->pushbackToBranch("pLS_circleCenterX", pixelSegments.circleCenterX()[ipLS]); + ana.tx->pushbackToBranch("pLS_circleCenterY", pixelSegments.circleCenterY()[ipLS]); + ana.tx->pushbackToBranch("pLS_circleRadius", pixelSegments.circleRadius()[ipLS]); + ana.tx->pushbackToBranch("pLS_px", pixelSeeds.px()[ipLS]); + ana.tx->pushbackToBranch("pLS_py", pixelSeeds.py()[ipLS]); + ana.tx->pushbackToBranch("pLS_pz", pixelSeeds.pz()[ipLS]); + ana.tx->pushbackToBranch("pLS_isQuad", static_cast(pixelSeeds.isQuad()[ipLS])); + ana.tx->pushbackToBranch("pLS_nhit", hit_idx.size()); + for (size_t ihit = 0; ihit < trk_see_hitIdx[ipLS].size(); ++ihit) { + int hitidx = trk_see_hitIdx[ipLS][ihit]; + int hittype = trk_see_hitType[ipLS][ihit]; + int x = trk_pix_x[hitidx]; + int y = trk_pix_y[hitidx]; + int z = trk_pix_z[hitidx]; + ana.tx->pushbackToBranch(TString::Format("pLS_hit%d_x", ihit), x); + ana.tx->pushbackToBranch(TString::Format("pLS_hit%d_y", ihit), y); + ana.tx->pushbackToBranch(TString::Format("pLS_hit%d_z", ihit), z); + } + if (trk_see_hitIdx[ipLS].size() == 3) { + ana.tx->pushbackToBranch("pLS_hit3_x", -999); + ana.tx->pushbackToBranch("pLS_hit3_y", -999); + ana.tx->pushbackToBranch("pLS_hit3_z", -999); + } + bool isfake = true; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + if (simidxfrac[isim] > matchfrac) { + isfake = false; + break; + } + } + ana.tx->pushbackToBranch("pLS_isFake", isfake); + pls_simIdxAll.push_back(simidx); + pls_simIdxAllFrac.push_back(simidxfrac); + for (size_t is = 0; is < simidx.size(); ++is) { + int sim_idx = simidx.at(is); + if (sim_idx < n_accepted_simtrk) { + sim_pLS_matched[sim_idx]++; + } + float sim_idx_frac = simidxfrac.at(is); + if (sim_idx < n_total_simtrk) { + sim_plsIdxAll.at(sim_idx).push_back(pls_idx); + sim_plsIdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } + } + int pls_simIdx = -999; + float pls_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > pls_simIdxBestFrac and thisfrac > matchfrac) { + pls_simIdxBestFrac = thisfrac; + pls_simIdx = thisidx; + } + } + ana.tx->pushbackToBranch("pLS_simIdx", pls_simIdx); + // count global + pls_idx++; + } + ana.tx->setBranch>>("pLS_simIdxAll", pls_simIdxAll); + ana.tx->setBranch>>("pLS_simIdxAllFrac", pls_simIdxAllFrac); + std::vector> sim_plsIdxAll_to_write; + std::vector> sim_plsIdxAllFrac_to_write; + std::copy( + sim_plsIdxAll.begin(), sim_plsIdxAll.begin() + n_accepted_simtrk, std::back_inserter(sim_plsIdxAll_to_write)); + std::copy(sim_plsIdxAllFrac.begin(), + sim_plsIdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_plsIdxAllFrac_to_write)); + ana.tx->setBranch>>("sim_plsIdxAll", sim_plsIdxAll_to_write); + ana.tx->setBranch>>("sim_plsIdxAllFrac", sim_plsIdxAllFrac_to_write); + + std::vector pLS_isDuplicate(pls_simIdxAll.size(), 0); + for (size_t i = 0; i < pls_simIdxAll.size(); ++i) { + for (int simidx : pls_simIdxAll[i]) { + if (simidx < n_accepted_simtrk && sim_pLS_matched[simidx] > 1) { + pLS_isDuplicate[i] = 1; + break; + } + } + } + ana.tx->setBranch>("pLS_isDuplicate", pLS_isDuplicate); + + return pls_idx_map; +} + +//________________________________________________________________________________________________________________________________ +std::map setPixelTripletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& pls_idx_map, + std::map const& t3_idx_map) { + //-------------------------------------------- + // + // + // pT3 + // + // + //-------------------------------------------- + + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); + auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); + auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& pixelSeeds = event->getInput(); + auto const& pixelTriplets = event->getPixelTriplets(); + + int n_total_simtrk = trk_sim_pt.size(); + std::vector sim_pT3_matched(n_accepted_simtrk, 0); + std::vector> sim_pt3IdxAll(n_total_simtrk); + std::vector> sim_pt3IdxAllFrac(n_total_simtrk); + std::vector> pt3_simIdxAll; + std::vector> pt3_simIdxAllFrac; + // Then obtain the lower module index + unsigned int pt3_idx = 0; // global pt3 index that will be used to keep track of pt3 being outputted to the ntuple + // map to keep track of (GPU pt3Idx) -> (pt3_idx in ntuple output) + std::map pt3_idx_map; + // printT3s(event); + unsigned int nPixelTriplets = pixelTriplets.nPixelTriplets(); + for (unsigned int ipT3 = 0; ipT3 < nPixelTriplets; ipT3++) { + unsigned int pt3Idx = ipT3; + pt3_idx_map[pt3Idx] = pt3_idx; + std::vector hit_idx, hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFrompT3(event, ipT3); + std::vector simidx; + std::vector simidxfrac; + std::tie(simidx, simidxfrac) = + matchedSimTrkIdxsAndFracs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); + // // Computing line segment pt estimate (assuming beam spot is at zero) + unsigned int ipLS = getPixelLSFrompT3(event, ipT3); + float pt = pixelSeeds.ptIn()[ipLS]; + float eta = pixelSeeds.eta()[ipLS]; + float phi = pixelSeeds.phi()[ipLS]; + ana.tx->pushbackToBranch("pT3_pt", pt); + ana.tx->pushbackToBranch("pT3_eta", eta); + ana.tx->pushbackToBranch("pT3_phi", phi); + ana.tx->pushbackToBranch("pT3_score", pixelTriplets.score()[ipT3]); + if (ana.pls_branches) { + unsigned int plsIdx = ranges.segmentModuleIndices()[modules.nLowerModules()] + ipLS; + unsigned int pls_idx = pls_idx_map.at(plsIdx); + ana.tx->pushbackToBranch("pT3_plsIdx", pls_idx); + } + if (ana.t3_branches) { + unsigned int t3Idx = getT3FrompT3(event, ipT3); + unsigned int t3_idx = t3_idx_map.at(t3Idx); + ana.tx->pushbackToBranch("pT3_t3Idx", t3_idx); + } + bool isfake = true; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + if (simidxfrac[isim] > matchfrac) { + isfake = false; + break; + } + } + ana.tx->pushbackToBranch("pT3_isFake", isfake); + pt3_simIdxAll.push_back(simidx); + pt3_simIdxAllFrac.push_back(simidxfrac); + for (size_t is = 0; is < simidx.size(); ++is) { + int sim_idx = simidx.at(is); + if (sim_idx < n_accepted_simtrk) { + sim_pT3_matched.at(sim_idx) += 1; + } + float sim_idx_frac = simidxfrac.at(is); + if (sim_idx < n_total_simtrk) { + sim_pt3IdxAll.at(sim_idx).push_back(pt3_idx); + sim_pt3IdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } } - std::vector simidx = - matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); - ana.tx->pushbackToBranch("pT5_isFake", static_cast(simidx.size() == 0)); + int pt3_simIdx = -999; + float pt3_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > pt3_simIdxBestFrac and thisfrac > matchfrac) { + pt3_simIdxBestFrac = thisfrac; + pt3_simIdx = thisidx; + } + } + ana.tx->pushbackToBranch("pT3_simIdx", pt3_simIdx); + // count global + pt3_idx++; + } + ana.tx->setBranch>>("pT3_simIdxAll", pt3_simIdxAll); + ana.tx->setBranch>>("pT3_simIdxAllFrac", pt3_simIdxAllFrac); + std::vector> sim_pt3IdxAll_to_write; + std::vector> sim_pt3IdxAllFrac_to_write; + std::copy( + sim_pt3IdxAll.begin(), sim_pt3IdxAll.begin() + n_accepted_simtrk, std::back_inserter(sim_pt3IdxAll_to_write)); + std::copy(sim_pt3IdxAllFrac.begin(), + sim_pt3IdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_pt3IdxAllFrac_to_write)); + ana.tx->setBranch>>("sim_pt3IdxAll", sim_pt3IdxAll_to_write); + ana.tx->setBranch>>("sim_pt3IdxAllFrac", sim_pt3IdxAllFrac_to_write); + + std::vector pT3_isDuplicate(pt3_simIdxAll.size()); + for (unsigned int i = 0; i < pt3_simIdxAll.size(); i++) { + bool isDuplicate = true; + for (unsigned int isim = 0; isim < pt3_simIdxAll[i].size(); isim++) { + int simidx = pt3_simIdxAll[i][isim]; + if (simidx < n_accepted_simtrk) { + if (sim_pT3_matched[simidx] > 1) { + isDuplicate = true; + } + } + } + pT3_isDuplicate[i] = isDuplicate; + } + ana.tx->setBranch>("sim_pT3_matched", sim_pT3_matched); + ana.tx->setBranch>("pT3_isDuplicate", pT3_isDuplicate); + + return pt3_idx_map; +} + +//________________________________________________________________________________________________________________________________ +std::map setPixelQuintupletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& pls_idx_map, + std::map const& t5_idx_map) { + //-------------------------------------------- + // + // + // pT5 + // + // + //-------------------------------------------- + + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); + auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); + auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& pixelSeeds = event->getInput(); + auto const& quintuplets = event->getQuintuplets(); + auto const& pixelQuintuplets = event->getPixelQuintuplets(); + + int n_total_simtrk = trk_sim_pt.size(); + std::vector sim_pT5_matched(n_accepted_simtrk); + std::vector> sim_pt5IdxAll(n_total_simtrk); + std::vector> sim_pt5IdxAllFrac(n_total_simtrk); + std::vector> pt5_simIdxAll; + std::vector> pt5_simIdxAllFrac; + // Then obtain the lower module index + unsigned int pt5_idx = 0; // global pt5 index that will be used to keep track of pt5 being outputted to the ntuple + // map to keep track of (GPU pt5Idx) -> (pt5_idx in ntuple output) + std::map pt5_idx_map; + // printT5s(event); + unsigned int nPixelQuintuplets = pixelQuintuplets.nPixelQuintuplets(); + for (unsigned int ipT5 = 0; ipT5 < nPixelQuintuplets; ipT5++) { + unsigned int pt5Idx = ipT5; + pt5_idx_map[pt5Idx] = pt5_idx; + std::vector hit_idx, hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFrompT5(event, ipT5); + std::vector simidx; + std::vector simidxfrac; + std::tie(simidx, simidxfrac) = + matchedSimTrkIdxsAndFracs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); + // // Computing line segment pt estimate (assuming beam spot is at zero) + unsigned int T5Index = getT5FrompT5(event, ipT5); + unsigned int ipLS = getPixelLSFrompT5(event, ipT5); + float pt = (__H2F(quintuplets.innerRadius()[T5Index]) * k2Rinv1GeVf * 2 + pixelSeeds.ptIn()[ipLS]) / 2; + float eta = pixelSeeds.eta()[ipLS]; + float phi = pixelSeeds.phi()[ipLS]; ana.tx->pushbackToBranch("pT5_pt", pt); ana.tx->pushbackToBranch("pT5_eta", eta); ana.tx->pushbackToBranch("pT5_phi", phi); - ana.tx->pushbackToBranch("pT5_layer_binary", layer_binary); - ana.tx->pushbackToBranch("pT5_moduleType_binary", moduleType_binary); - ana.tx->pushbackToBranch("pT5_rzChiSquared", pixelQuintuplets.rzChiSquared()[pT5]); - - pT5_matched_simIdx.push_back(simidx); - - // Loop over matched sim idx and increase counter of pT5_matched - for (auto& idx : simidx) { - // NOTE Important to note that the idx of the std::vector<> is same - // as the tracking-ntuple's sim track idx ONLY because event==0 and bunchCrossing==0 condition is applied!! - // Also do not try to access beyond the event and bunchCrossing - if (idx < n_accepted_simtrk) { - sim_pT5_matched.at(idx) += 1; + if (ana.pls_branches) { + unsigned int plsIdx = ranges.segmentModuleIndices()[modules.nLowerModules()] + ipLS; + unsigned int pls_idx = pls_idx_map.at(plsIdx); + ana.tx->pushbackToBranch("pT5_plsIdx", pls_idx); + } + if (ana.t5_branches) { + unsigned int t5Idx = getT5FrompT5(event, ipT5); + unsigned int t5_idx = t5_idx_map.at(t5Idx); + ana.tx->pushbackToBranch("pT5_t5Idx", t5_idx); + } + bool isfake = true; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + if (simidxfrac[isim] > matchfrac) { + isfake = false; + break; } } + ana.tx->pushbackToBranch("pT5_isFake", isfake); + pt5_simIdxAll.push_back(simidx); + pt5_simIdxAllFrac.push_back(simidxfrac); + for (size_t is = 0; is < simidx.size(); ++is) { + int sim_idx = simidx.at(is); + if (sim_idx < n_accepted_simtrk) { + sim_pT5_matched.at(sim_idx) += 1; + } + float sim_idx_frac = simidxfrac.at(is); + if (sim_idx < n_total_simtrk) { + sim_pt5IdxAll.at(sim_idx).push_back(pt5_idx); + sim_pt5IdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } + } + int pt5_simIdx = -999; + float pt5_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > pt5_simIdxBestFrac and thisfrac > matchfrac) { + pt5_simIdxBestFrac = thisfrac; + pt5_simIdx = thisidx; + } + } + ana.tx->pushbackToBranch("pT5_simIdx", pt5_simIdx); + // count global + pt5_idx++; } + ana.tx->setBranch>>("pT5_simIdxAll", pt5_simIdxAll); + ana.tx->setBranch>>("pT5_simIdxAllFrac", pt5_simIdxAllFrac); + std::vector> sim_pt5IdxAll_to_write; + std::vector> sim_pt5IdxAllFrac_to_write; + std::copy( + sim_pt5IdxAll.begin(), sim_pt5IdxAll.begin() + n_accepted_simtrk, std::back_inserter(sim_pt5IdxAll_to_write)); + std::copy(sim_pt5IdxAllFrac.begin(), + sim_pt5IdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_pt5IdxAllFrac_to_write)); + ana.tx->setBranch>>("sim_pt5IdxAll", sim_pt5IdxAll_to_write); + ana.tx->setBranch>>("sim_pt5IdxAllFrac", sim_pt5IdxAllFrac_to_write); // Using the intermedaite variables to compute whether a given track candidate is a duplicate - std::vector pT5_isDuplicate(pT5_matched_simIdx.size()); + std::vector pT5_isDuplicate(pt5_simIdxAll.size()); // Loop over the track candidates - for (unsigned int i = 0; i < pT5_matched_simIdx.size(); ++i) { + for (unsigned int i = 0; i < pt5_simIdxAll.size(); ++i) { bool isDuplicate = false; // Loop over the sim idx matched to this track candidate - for (unsigned int isim = 0; isim < pT5_matched_simIdx[i].size(); ++isim) { + for (unsigned int isim = 0; isim < pt5_simIdxAll[i].size(); ++isim) { // Using the sim_pT5_matched to see whether this track candidate is matched to a sim track that is matched to more than one - int simidx = pT5_matched_simIdx[i][isim]; + int simidx = pt5_simIdxAll[i][isim]; if (simidx < n_accepted_simtrk) { if (sim_pT5_matched[simidx] > 1) { isDuplicate = true; @@ -611,181 +1802,288 @@ void setPixelQuintupletOutputBranches(LSTEvent* event) { } pT5_isDuplicate[i] = isDuplicate; } - - // Now set the last remaining branches - ana.tx->setBranch>("sim_pT5_matched", sim_pT5_matched); - ana.tx->setBranch>>("pT5_matched_simIdx", pT5_matched_simIdx); ana.tx->setBranch>("pT5_isDuplicate", pT5_isDuplicate); + + return pt5_idx_map; } //________________________________________________________________________________________________________________________________ -void setQuintupletOutputBranches(LSTEvent* event) { - auto const quintuplets = event->getQuintuplets(); - auto const quintupletsOccupancy = event->getQuintuplets(); - auto ranges = event->getRanges(); - auto modules = event->getModules(); - int n_accepted_simtrk = ana.tx->getBranch>("sim_TC_matched").size(); - - std::vector sim_t5_matched(n_accepted_simtrk); - std::vector> t5_matched_simIdx; +void setTrackCandidateBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + std::map t5_idx_map, + std::map pls_idx_map, + std::map pt3_idx_map, + std::map pt5_idx_map, + float matchfrac) { + //-------------------------------------------- + // + // + // Track Candidates + // + // + //-------------------------------------------- - auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx"); - auto const& trk_simvtx_x = trk.getVF("simvtx_x"); - auto const& trk_simvtx_y = trk.getVF("simvtx_y"); - auto const& trk_simvtx_z = trk.getVF("simvtx_z"); + auto const& trk_sim_pt = trk.getVF("sim_pt"); + auto const& trk_ph2_x = trk.getVF("ph2_x"); + auto const& trk_ph2_y = trk.getVF("ph2_y"); + auto const& trk_ph2_z = trk.getVF("ph2_z"); auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); - for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < modules.nLowerModules(); ++lowerModuleIdx) { - int nQuintuplets = quintupletsOccupancy.nQuintuplets()[lowerModuleIdx]; - for (unsigned int idx = 0; idx < nQuintuplets; idx++) { - unsigned int quintupletIndex = ranges.quintupletModuleIndices()[lowerModuleIdx] + idx; - float pt = __H2F(quintuplets.innerRadius()[quintupletIndex]) * k2Rinv1GeVf * 2; - float eta = __H2F(quintuplets.eta()[quintupletIndex]); - float phi = __H2F(quintuplets.phi()[quintupletIndex]); + auto const& ranges = event->getRanges(); + auto const& modules = event->getModules(); + auto const& trackCandidatesBase = event->getTrackCandidatesBase(); + auto const& trackCandidatesExtended = event->getTrackCandidatesExtended(); + + // Following are some vectors to keep track of the information to write to the ntuple + // N.B. following two branches have a length for the entire sim track, but what actually will be written in sim_tcIdxAll branch is NOT that long + // Later in the code, it will restrict to only the ones to write out. + // The reason at this stage, the entire tcIdxAll is being tracked is to compute duplicate properly later on + // When computing a duplicate object it is important to consider all simulated tracks including pileup tracks + int n_total_simtrk = trk_sim_pt.size(); + std::vector> sim_tcIdxAll(n_total_simtrk); + std::vector> sim_tcIdxAllFrac(n_total_simtrk); + std::vector> tc_simIdxAll; + std::vector> tc_simIdxAllFrac; + + // Number of total track candidates created in this event + unsigned int nTrackCandidates = trackCandidatesBase.nTrackCandidates(); + + // Looping over each track candidate + for (unsigned int tc_idx = 0; tc_idx < nTrackCandidates; tc_idx++) { + // Compute reco quantities of track candidate based on final object + int type, isFake; + float pt, eta, phi; + std::vector simidx; // list of all the matched sim idx + std::vector simidxfrac; // list of match fraction for each matched sim idx + + // The following function reads off and computes the matched sim track indices + std::tie(type, pt, eta, phi, isFake, simidx, simidxfrac) = parseTrackCandidateAllMatch(event, + tc_idx, + trk_ph2_x, + trk_ph2_y, + trk_ph2_z, + trk_simhit_simTrkIdx, + trk_ph2_simHitIdx, + trk_pix_simHitIdx, + matchfrac); + + // Fill some branches for this track candidate + ana.tx->pushbackToBranch("tc_pt", pt); + ana.tx->pushbackToBranch("tc_eta", eta); + ana.tx->pushbackToBranch("tc_phi", phi); + ana.tx->pushbackToBranch("tc_type", type); + if (type == LSTObjType::pT5) { + if (ana.pt5_branches) + ana.tx->pushbackToBranch( + "tc_pt5Idx", + (ana.pt5_branches ? pt5_idx_map[trackCandidatesExtended.directObjectIndices()[tc_idx]] : -999)); + if (ana.pt3_branches) + ana.tx->pushbackToBranch("tc_pt3Idx", -999); + if (ana.t5_branches) + ana.tx->pushbackToBranch("tc_t5Idx", -999); + if (ana.pls_branches) + ana.tx->pushbackToBranch("tc_plsIdx", -999); + } else if (type == LSTObjType::pT3) { + if (ana.pt5_branches) + ana.tx->pushbackToBranch("tc_pt5Idx", -999); + if (ana.pt3_branches) + ana.tx->pushbackToBranch( + "tc_pt3Idx", + (ana.pt3_branches ? pt3_idx_map[trackCandidatesExtended.directObjectIndices()[tc_idx]] : -999)); + if (ana.t5_branches) + ana.tx->pushbackToBranch("tc_t5Idx", -999); + if (ana.pls_branches) + ana.tx->pushbackToBranch("tc_plsIdx", -999); + } else if (type == LSTObjType::T5) { + if (ana.pt5_branches) + ana.tx->pushbackToBranch("tc_pt5Idx", -999); + if (ana.pt3_branches) + ana.tx->pushbackToBranch("tc_pt3Idx", -999); + if (ana.t5_branches) + ana.tx->pushbackToBranch( + "tc_t5Idx", (ana.t5_branches ? t5_idx_map[trackCandidatesExtended.directObjectIndices()[tc_idx]] : -999)); + if (ana.pls_branches) + ana.tx->pushbackToBranch("tc_plsIdx", -999); + } else if (type == LSTObjType::pLS) { + if (ana.pt5_branches) + ana.tx->pushbackToBranch("tc_pt5Idx", -999); + if (ana.pt3_branches) + ana.tx->pushbackToBranch("tc_pt3Idx", -999); + if (ana.t5_branches) + ana.tx->pushbackToBranch("tc_t5Idx", -999); + if (ana.pls_branches) + ana.tx->pushbackToBranch( + "tc_plsIdx", + (ana.pls_branches ? pls_idx_map[ranges.segmentModuleIndices()[modules.nLowerModules()] + + trackCandidatesExtended.directObjectIndices()[tc_idx]] + : -999)); + } + + ana.tx->pushbackToBranch("tc_isFake", isFake); - std::vector hit_idx = getHitIdxsFromT5(event, quintupletIndex); - std::vector hit_type = getHitTypesFromT5(event, quintupletIndex); - std::vector module_idx = getModuleIdxsFromT5(event, quintupletIndex); + // For this tc, keep track of all the simidx that are matched + tc_simIdxAll.push_back(simidx); + tc_simIdxAllFrac.push_back(simidxfrac); + + // The book keeping of opposite mapping is done here + // For each matched sim idx, we go back and keep track of which tc it is matched to. + // Loop over all the matched sim idx + for (size_t is = 0; is < simidx.size(); ++is) { + // For this matched sim index keep track (sim -> tc) mapping + int sim_idx = simidx.at(is); + float sim_idx_frac = simidxfrac.at(is); + sim_tcIdxAll.at(sim_idx).push_back(tc_idx); + sim_tcIdxAllFrac.at(sim_idx).push_back(sim_idx_frac); + } - int layer_binary = 0; - int moduleType_binary = 0; - for (size_t i = 0; i < module_idx.size(); i += 2) { - layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4))); - moduleType_binary |= (modules.moduleType()[module_idx[i]] << i); + // Also, among the simidx matches, find the best match (highest fractional match) + // N.B. the simidx is already returned sorted by highest number of "nhits" match + // So as it loops over, the condition will ensure that the highest fraction with highest nhits will be matched with the priority given to highest fraction + int tc_simIdx = -999; + float tc_simIdxBestFrac = 0; + for (size_t isim = 0; isim < simidx.size(); ++isim) { + int thisidx = simidx[isim]; + float thisfrac = simidxfrac[isim]; + if (thisfrac > tc_simIdxBestFrac and thisfrac > matchfrac) { + tc_simIdxBestFrac = thisfrac; + tc_simIdx = thisidx; } + } - float percent_matched; - std::vector simidx = matchedSimTrkIdxs( - hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, &percent_matched); - - ana.tx->pushbackToBranch("t5_isFake", static_cast(simidx.size() == 0)); - ana.tx->pushbackToBranch("t5_pt", pt); - ana.tx->pushbackToBranch("t5_pMatched", percent_matched); - ana.tx->pushbackToBranch("t5_eta", eta); - ana.tx->pushbackToBranch("t5_phi", phi); - ana.tx->pushbackToBranch("t5_innerRadius", __H2F(quintuplets.innerRadius()[quintupletIndex])); - ana.tx->pushbackToBranch("t5_bridgeRadius", __H2F(quintuplets.bridgeRadius()[quintupletIndex])); - ana.tx->pushbackToBranch("t5_outerRadius", __H2F(quintuplets.outerRadius()[quintupletIndex])); - ana.tx->pushbackToBranch("t5_isDupAlgoFlag", quintuplets.isDup()[quintupletIndex]); - ana.tx->pushbackToBranch("t5_chiSquared", quintuplets.chiSquared()[quintupletIndex]); - ana.tx->pushbackToBranch("t5_rzChiSquared", quintuplets.rzChiSquared()[quintupletIndex]); - ana.tx->pushbackToBranch("t5_nonAnchorChiSquared", quintuplets.nonAnchorChiSquared()[quintupletIndex]); - ana.tx->pushbackToBranch("t5_dBeta1", quintuplets.dBeta1()[quintupletIndex]); - ana.tx->pushbackToBranch("t5_dBeta2", quintuplets.dBeta2()[quintupletIndex]); - ana.tx->pushbackToBranch("t5_layer_binary", layer_binary); - ana.tx->pushbackToBranch("t5_moduleType_binary", moduleType_binary); - - t5_matched_simIdx.push_back(simidx); - - for (auto& simtrk : simidx) { - if (simtrk < n_accepted_simtrk) { - sim_t5_matched.at(simtrk) += 1; - } - } + // the best match index will then be saved here + ana.tx->pushbackToBranch("tc_simIdx", tc_simIdx); + } - // Avoid fakes when calculating the vertex distance, set default to 0.0. - if (simidx.size() == 0) { - ana.tx->pushbackToBranch("t5_sim_vxy", 0.0); - ana.tx->pushbackToBranch("t5_sim_vz", 0.0); - continue; - } + // Now save the (tc -> simidx) mapping + ana.tx->setBranch>>("tc_simIdxAll", tc_simIdxAll); + ana.tx->setBranch>>("tc_simIdxAllFrac", tc_simIdxAllFrac); + + // Not all (sim->tcIdx) will be saved but only for the sim that is from hard scatter and current bunch crossing + // So a restriction up to only "n_accepted_simtrk" done by chopping off the rest + // N.B. the reason we can simply take the first "n_accepted_simtrk" is because the tracking ntuple is organized such that those sim tracks show up on the first "n_accepted_simtrk" of tracks. + std::vector> sim_tcIdxAll_to_write; + std::vector> sim_tcIdxAllFrac_to_write; + std::copy(sim_tcIdxAll.begin(), + sim_tcIdxAll.begin() + n_accepted_simtrk, + std::back_inserter( + sim_tcIdxAll_to_write)); // this is where the vector is only copying the first "n_accepted_simtrk" + std::copy(sim_tcIdxAllFrac.begin(), + sim_tcIdxAllFrac.begin() + n_accepted_simtrk, + std::back_inserter(sim_tcIdxAllFrac_to_write)); // ditto + ana.tx->setBranch>>("sim_tcIdxAll", sim_tcIdxAll_to_write); + ana.tx->setBranch>>("sim_tcIdxAllFrac", sim_tcIdxAllFrac_to_write); - int vtxidx = trk_sim_parentVtxIdx[simidx[0]]; - float vtx_x = trk_simvtx_x[vtxidx]; - float vtx_y = trk_simvtx_y[vtxidx]; - float vtx_z = trk_simvtx_z[vtxidx]; + // Using the intermedaite variables to compute whether a given track candidate is a duplicate + std::vector tc_isDuplicate(tc_simIdxAll.size()); - ana.tx->pushbackToBranch("t5_sim_vxy", sqrt(vtx_x * vtx_x + vtx_y * vtx_y)); - ana.tx->pushbackToBranch("t5_sim_vz", vtx_z); + // Loop over the track candidates + for (unsigned int tc_idx = 0; tc_idx < tc_simIdxAll.size(); ++tc_idx) { + bool isDuplicate = false; + // Loop over the sim idx matched to this track candidate + for (unsigned int isim = 0; isim < tc_simIdxAll[tc_idx].size(); ++isim) { + int sim_idx = tc_simIdxAll[tc_idx][isim]; + int n_sim_matched = 0; + for (size_t ism = 0; ism < sim_tcIdxAll.at(sim_idx).size(); ++ism) { + if (sim_tcIdxAllFrac.at(sim_idx).at(ism) > matchfrac) { + n_sim_matched += 1; + if (n_sim_matched > 1) { + isDuplicate = true; + break; + } + } + } } + tc_isDuplicate[tc_idx] = isDuplicate; } + ana.tx->setBranch>("tc_isDuplicate", tc_isDuplicate); - std::vector t5_isDuplicate(t5_matched_simIdx.size()); - for (unsigned int i = 0; i < t5_matched_simIdx.size(); i++) { - bool isDuplicate = false; - for (unsigned int isim = 0; isim < t5_matched_simIdx[i].size(); isim++) { - int simidx = t5_matched_simIdx[i][isim]; - if (simidx < n_accepted_simtrk) { - if (sim_t5_matched[simidx] > 1) { - isDuplicate = true; - } + // Similarly, the best match for the (sim -> tc is computed) + // TODO: Is this redundant? I am not sure if it is guaranteed that sim_tcIdx will have same result with tc_simIdx. + // I think it will be, but I have not rigorously checked. I only checked about first few thousands and it was all true. as long as tc->sim was pointing to a sim that is among the n_accepted. + // For the most part I think this won't be a problem. + for (size_t i = 0; i < sim_tcIdxAll_to_write.size(); ++i) { + // bestmatch is not always the first one + int bestmatch_idx = -999; + float bestmatch_frac = -999; + for (size_t jj = 0; jj < sim_tcIdxAll_to_write.at(i).size(); ++jj) { + int idx = sim_tcIdxAll_to_write.at(i).at(jj); + float frac = sim_tcIdxAllFrac_to_write.at(i).at(jj); + if (bestmatch_frac < frac) { + bestmatch_idx = idx; + bestmatch_frac = frac; } } - t5_isDuplicate[i] = isDuplicate; + ana.tx->pushbackToBranch("sim_tcIdxBest", bestmatch_idx); + ana.tx->pushbackToBranch("sim_tcIdxBestFrac", bestmatch_frac); + if (bestmatch_frac > matchfrac) // then this is a good match according to MTV + ana.tx->pushbackToBranch("sim_tcIdx", bestmatch_idx); + else + ana.tx->pushbackToBranch("sim_tcIdx", -999); } - ana.tx->setBranch>("sim_T5_matched", sim_t5_matched); - ana.tx->setBranch>>("t5_matched_simIdx", t5_matched_simIdx); - ana.tx->setBranch>("t5_isDuplicate", t5_isDuplicate); } //________________________________________________________________________________________________________________________________ -void setPixelTripletOutputBranches(LSTEvent* event) { - auto const pixelTriplets = event->getPixelTriplets(); +void setOccupancyBranches(LSTEvent* event) { auto modules = event->getModules(); - PixelSeedsConst pixelSeeds = event->getInput(); - int n_accepted_simtrk = ana.tx->getBranch>("sim_TC_matched").size(); - - unsigned int nPixelTriplets = pixelTriplets.nPixelTriplets(); - std::vector sim_pT3_matched(n_accepted_simtrk); - std::vector> pT3_matched_simIdx; - - auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); - auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); - auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + auto miniDoublets = event->getMiniDoublets(); + auto segments = event->getSegments(); + auto triplets = event->getTriplets(); + auto quintuplets = event->getQuintuplets(); + auto pixelQuintuplets = event->getPixelQuintuplets(); + auto pixelTriplets = event->getPixelTriplets(); + auto trackCandidatesBase = event->getTrackCandidatesBase(); - for (unsigned int pT3 = 0; pT3 < nPixelTriplets; pT3++) { - unsigned int T3Index = getT3FrompT3(event, pT3); - unsigned int pLSIndex = getPixelLSFrompT3(event, pT3); - const float pt = pixelSeeds.ptIn()[pLSIndex]; - - float eta = pixelSeeds.eta()[pLSIndex]; - float phi = pixelSeeds.phi()[pLSIndex]; - std::vector hit_idx = getHitIdxsFrompT3(event, pT3); - std::vector hit_type = getHitTypesFrompT3(event, pT3); - - std::vector simidx = - matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); - std::vector module_idx = getModuleIdxsFrompT3(event, pT3); - int layer_binary = 1; - int moduleType_binary = 0; - for (size_t i = 0; i < module_idx.size(); i += 2) { - layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4))); - moduleType_binary |= (modules.moduleType()[module_idx[i]] << i); - } - ana.tx->pushbackToBranch("pT3_isFake", static_cast(simidx.size() == 0)); - ana.tx->pushbackToBranch("pT3_pt", pt); - ana.tx->pushbackToBranch("pT3_eta", eta); - ana.tx->pushbackToBranch("pT3_phi", phi); - ana.tx->pushbackToBranch("pT3_layer_binary", layer_binary); - ana.tx->pushbackToBranch("pT3_moduleType_binary", moduleType_binary); + std::vector moduleLayer; + std::vector moduleSubdet; + std::vector moduleRing; + std::vector moduleRod; + std::vector moduleModule; + std::vector moduleEta; + std::vector moduleR; + std::vector moduleIsTilted; + std::vector trackCandidateOccupancy; + std::vector tripletOccupancy; + std::vector segmentOccupancy; + std::vector mdOccupancy; + std::vector quintupletOccupancy; - pT3_matched_simIdx.push_back(simidx); + for (unsigned int lowerIdx = 0; lowerIdx <= modules.nLowerModules(); lowerIdx++) { + //layer = 0, subdet = 0 => pixel module + moduleLayer.push_back(modules.layers()[lowerIdx]); + moduleSubdet.push_back(modules.subdets()[lowerIdx]); + moduleRing.push_back(modules.rings()[lowerIdx]); + moduleRod.push_back(modules.rods()[lowerIdx]); + moduleEta.push_back(modules.eta()[lowerIdx]); + moduleR.push_back(modules.r()[lowerIdx]); + bool isTilted = (modules.subdets()[lowerIdx] == 5 and modules.sides()[lowerIdx] != 3); + moduleIsTilted.push_back(isTilted); + moduleModule.push_back(modules.modules()[lowerIdx]); + segmentOccupancy.push_back(segments.totOccupancySegments()[lowerIdx]); + mdOccupancy.push_back(miniDoublets.totOccupancyMDs()[lowerIdx]); - for (auto& idx : simidx) { - if (idx < n_accepted_simtrk) { - sim_pT3_matched.at(idx) += 1; - } + if (lowerIdx < modules.nLowerModules()) { + quintupletOccupancy.push_back(quintuplets.totOccupancyQuintuplets()[lowerIdx]); + tripletOccupancy.push_back(triplets.totOccupancyTriplets()[lowerIdx]); } } - std::vector pT3_isDuplicate(pT3_matched_simIdx.size()); - for (unsigned int i = 0; i < pT3_matched_simIdx.size(); i++) { - bool isDuplicate = true; - for (unsigned int isim = 0; isim < pT3_matched_simIdx[i].size(); isim++) { - int simidx = pT3_matched_simIdx[i][isim]; - if (simidx < n_accepted_simtrk) { - if (sim_pT3_matched[simidx] > 1) { - isDuplicate = true; - } - } - } - pT3_isDuplicate[i] = isDuplicate; - } - ana.tx->setBranch>("sim_pT3_matched", sim_pT3_matched); - ana.tx->setBranch>>("pT3_matched_simIdx", pT3_matched_simIdx); - ana.tx->setBranch>("pT3_isDuplicate", pT3_isDuplicate); + ana.tx->setBranch>("module_layers", moduleLayer); + ana.tx->setBranch>("module_subdets", moduleSubdet); + ana.tx->setBranch>("module_rings", moduleRing); + ana.tx->setBranch>("module_rods", moduleRod); + ana.tx->setBranch>("module_modules", moduleModule); + ana.tx->setBranch>("module_isTilted", moduleIsTilted); + ana.tx->setBranch>("module_eta", moduleEta); + ana.tx->setBranch>("module_r", moduleR); + ana.tx->setBranch>("md_occupancies", mdOccupancy); + ana.tx->setBranch>("sg_occupancies", segmentOccupancy); + ana.tx->setBranch>("t3_occupancies", tripletOccupancy); + ana.tx->setBranch("tc_occupancies", trackCandidatesBase.nTrackCandidates()); + ana.tx->setBranch("pT3_occupancies", pixelTriplets.totOccupancyPixelTriplets()); + ana.tx->setBranch>("t5_occupancies", quintupletOccupancy); + ana.tx->setBranch("pT5_occupancies", pixelQuintuplets.totOccupancyPixelQuintuplets()); } //________________________________________________________________________________________________________________________________ @@ -817,17 +2115,17 @@ void fillpT3DNNBranches(LSTEvent* event, unsigned int iPT3) { //________________________________________________________________________________________________________________________________ void fillT3DNNBranches(LSTEvent* event, unsigned int iT3) { - auto hitsBase = event->getInput(); - auto hitsExtended = event->getHits(); - auto modules = event->getModules(); - - std::vector hitIdx = getHitsFromT3(event, iT3); - std::vector hitObjects; - auto const& trk_ph2_subdet = trk.getVUS("ph2_subdet"); auto const& trk_ph2_layer = trk.getVUS("ph2_layer"); auto const& trk_ph2_detId = trk.getVU("ph2_detId"); + auto const& hitsBase = event->getInput(); + auto const& hitsExtended = event->getHits(); + auto const& modules = event->getModules(); + + std::vector hitIdx = getHitsFromT3(event, iT3); + std::vector hitObjects; + for (int i = 0; i < hitIdx.size(); ++i) { unsigned int hit = hitIdx[i]; float x = hitsBase.xs()[hit]; @@ -921,12 +2219,7 @@ void setpT3DNNBranches(LSTEvent* event) { } //________________________________________________________________________________________________________________________________ -void setT3DNNBranches(LSTEvent* event) { - auto const triplets = event->getTriplets(); - auto const tripletsOccupancy = event->getTriplets(); - auto modules = event->getModules(); - auto ranges = event->getRanges(); - +void setT3DNNBranches(LSTEvent* event, float matchfrac) { auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx"); auto const& trk_simvtx_x = trk.getVF("simvtx_x"); auto const& trk_simvtx_y = trk.getVF("simvtx_y"); @@ -935,6 +2228,11 @@ void setT3DNNBranches(LSTEvent* event) { auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); + auto const& triplets = event->getTriplets(); + auto const& tripletsOccupancy = event->getTriplets(); + auto const& modules = event->getModules(); + auto const& ranges = event->getRanges(); + for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < modules.nLowerModules(); ++lowerModuleIdx) { int nTriplets = tripletsOccupancy.nTriplets()[lowerModuleIdx]; for (unsigned int idx = 0; idx < nTriplets; idx++) { @@ -953,8 +2251,14 @@ void setT3DNNBranches(LSTEvent* event) { // Get matching information with percent matched float percent_matched; - std::vector simidx = matchedSimTrkIdxs( - hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, &percent_matched); + std::vector simidx = matchedSimTrkIdxs(hit_idx, + hit_type, + trk_simhit_simTrkIdx, + trk_ph2_simHitIdx, + trk_pix_simHitIdx, + false, + matchfrac, + &percent_matched); // Fill the branches with T3-specific data ana.tx->pushbackToBranch("t3_betaIn", triplets.betaIn()[tripletIndex]); @@ -1054,291 +2358,48 @@ void setT5DNNBranches(LSTEvent* event) { } //________________________________________________________________________________________________________________________________ -void setGnnNtupleBranches(LSTEvent* event) { - // Get relevant information - SegmentsOccupancyConst segmentsOccupancy = event->getSegments(); - MiniDoubletsOccupancyConst miniDoublets = event->getMiniDoublets(); - auto hitsBase = event->getInput(); - auto modules = event->getModules(); - auto ranges = event->getRanges(); +std::tuple> parseTrackCandidate( + LSTEvent* event, + unsigned int idx, + std::vector const& trk_ph2_x, + std::vector const& trk_ph2_y, + std::vector const& trk_ph2_z, + std::vector const& trk_simhit_simTrkIdx, + std::vector> const& trk_ph2_simHitIdx, + std::vector> const& trk_pix_simHitIdx, + float matchfrac) { + // Get the type of the track candidate auto const& trackCandidatesBase = event->getTrackCandidatesBase(); + short type = trackCandidatesBase.trackCandidateType()[idx]; - std::set mds_used_in_sg; - std::map md_index_map; - std::map sg_index_map; - - // Loop over modules (lower ones where the MDs are saved) - unsigned int nTotalMD = 0; - unsigned int nTotalLS = 0; - for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) { - nTotalMD += miniDoublets.nMDs()[idx]; - nTotalLS += segmentsOccupancy.nSegments()[idx]; - } - - auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); - auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); - auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); - - std::set lss_used_in_true_tc; - unsigned int nTrackCandidates = trackCandidatesBase.nTrackCandidates(); - for (unsigned int idx = 0; idx < nTrackCandidates; idx++) { - // Only consider true track candidates - std::vector hitidxs; - std::vector hittypes; - std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromTC(event, idx); - std::vector simidxs = - matchedSimTrkIdxs(hitidxs, hittypes, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); - if (simidxs.size() == 0) - continue; - - std::vector LSs = getLSsFromTC(event, idx); - for (auto& LS : LSs) { - if (lss_used_in_true_tc.find(LS) == lss_used_in_true_tc.end()) { - lss_used_in_true_tc.insert(LS); - } - } - } - - std::cout << " lss_used_in_true_tc.size(): " << lss_used_in_true_tc.size() << std::endl; - - // std::cout << " nTotalMD: " << nTotalMD << std::endl; - // std::cout << " nTotalLS: " << nTotalLS << std::endl; - - auto const& trk_sim_bunchCrossing = trk.getVI("sim_bunchCrossing"); - auto const& trk_sim_event = trk.getVI("sim_event"); - auto const& trk_sim_pt = trk.getVF("sim_pt"); - auto const& trk_sim_eta = trk.getVF("sim_eta"); - auto const& trk_sim_phi = trk.getVF("sim_phi"); - auto const& trk_sim_pca_dxy = trk.getVF("sim_pca_dxy"); - auto const& trk_sim_pca_dz = trk.getVF("sim_pca_dz"); - auto const& trk_sim_q = trk.getVI("sim_q"); - auto const& trk_sim_pdgId = trk.getVI("sim_pdgId"); - auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx"); - auto const& trk_simvtx_x = trk.getVF("simvtx_x"); - auto const& trk_simvtx_y = trk.getVF("simvtx_y"); - auto const& trk_simvtx_z = trk.getVF("simvtx_z"); - - // Loop over modules (lower ones where the MDs are saved) - for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) { - // // Loop over minidoublets - // for (unsigned int jdx = 0; jdx < miniDoublets->nMDs[idx]; jdx++) - // { - // // Get the actual index to the mini-doublet using ranges - // unsigned int mdIdx = ranges->miniDoubletModuleIndices[idx] + jdx; - - // setGnnNtupleMiniDoublet(event, mdIdx); - // } - - // Loop over segments - for (unsigned int jdx = 0; jdx < segmentsOccupancy.nSegments()[idx]; jdx++) { - // Get the actual index to the segments using ranges - unsigned int sgIdx = ranges.segmentModuleIndices()[idx] + jdx; - - // Get the hit indices - std::vector MDs = getMDsFromLS(event, sgIdx); - - if (mds_used_in_sg.find(MDs[0]) == mds_used_in_sg.end()) { - mds_used_in_sg.insert(MDs[0]); - md_index_map[MDs[0]] = mds_used_in_sg.size() - 1; - setGnnNtupleMiniDoublet(event, - MDs[0], - trk_sim_q, - trk_sim_pt, - trk_sim_eta, - trk_sim_bunchCrossing, - trk_sim_event, - trk_sim_parentVtxIdx, - trk_simvtx_x, - trk_simvtx_y, - trk_simvtx_z, - trk_simhit_simTrkIdx, - trk_ph2_simHitIdx, - trk_pix_simHitIdx); - } - - if (mds_used_in_sg.find(MDs[1]) == mds_used_in_sg.end()) { - mds_used_in_sg.insert(MDs[1]); - md_index_map[MDs[1]] = mds_used_in_sg.size() - 1; - setGnnNtupleMiniDoublet(event, - MDs[1], - trk_sim_q, - trk_sim_pt, - trk_sim_eta, - trk_sim_bunchCrossing, - trk_sim_event, - trk_sim_parentVtxIdx, - trk_simvtx_x, - trk_simvtx_y, - trk_simvtx_z, - trk_simhit_simTrkIdx, - trk_ph2_simHitIdx, - trk_pix_simHitIdx); - } - - ana.tx->pushbackToBranch("LS_MD_idx0", md_index_map[MDs[0]]); - ana.tx->pushbackToBranch("LS_MD_idx1", md_index_map[MDs[1]]); - - std::vector hits = getHitsFromLS(event, sgIdx); - - // Computing line segment pt estimate (assuming beam spot is at zero) - lst_math::Hit hitA(0, 0, 0); - lst_math::Hit hitB(hitsBase.xs()[hits[0]], hitsBase.ys()[hits[0]], hitsBase.zs()[hits[0]]); - lst_math::Hit hitC(hitsBase.xs()[hits[2]], hitsBase.ys()[hits[2]], hitsBase.zs()[hits[2]]); - lst_math::Hit center = lst_math::getCenterFromThreePoints(hitA, hitB, hitC); - float pt = lst_math::ptEstimateFromRadius(center.rt()); - float eta = hitC.eta(); - float phi = hitB.phi(); - - ana.tx->pushbackToBranch("LS_pt", pt); - ana.tx->pushbackToBranch("LS_eta", eta); - ana.tx->pushbackToBranch("LS_phi", phi); - // ana.tx->pushbackToBranch("LS_layer0", layer0); - // ana.tx->pushbackToBranch("LS_layer1", layer1); - - std::vector hitidxs; - std::vector hittypes; - std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromLS(event, sgIdx); - std::vector simidxs = - matchedSimTrkIdxs(hitidxs, hittypes, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); - - ana.tx->pushbackToBranch("LS_isFake", simidxs.size() == 0); - ana.tx->pushbackToBranch("LS_sim_pt", simidxs.size() > 0 ? trk_sim_pt[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_eta", simidxs.size() > 0 ? trk_sim_eta[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_phi", simidxs.size() > 0 ? trk_sim_phi[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_pca_dxy", simidxs.size() > 0 ? trk_sim_pca_dxy[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_pca_dz", simidxs.size() > 0 ? trk_sim_pca_dz[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_q", simidxs.size() > 0 ? trk_sim_q[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_event", simidxs.size() > 0 ? trk_sim_event[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_bx", simidxs.size() > 0 ? trk_sim_bunchCrossing[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_pdgId", simidxs.size() > 0 ? trk_sim_pdgId[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_vx", - simidxs.size() > 0 ? trk_simvtx_x[trk_sim_parentVtxIdx[simidxs[0]]] : -999); - ana.tx->pushbackToBranch("LS_sim_vy", - simidxs.size() > 0 ? trk_simvtx_y[trk_sim_parentVtxIdx[simidxs[0]]] : -999); - ana.tx->pushbackToBranch("LS_sim_vz", - simidxs.size() > 0 ? trk_simvtx_z[trk_sim_parentVtxIdx[simidxs[0]]] : -999); - ana.tx->pushbackToBranch("LS_isInTrueTC", lss_used_in_true_tc.find(sgIdx) != lss_used_in_true_tc.end()); - - sg_index_map[sgIdx] = ana.tx->getBranch>("LS_isFake").size() - 1; - - // // T5 eta and phi are computed using outer and innermost hits - // lst_math::Hit hitA(trk_ph2_x[anchitidx], trk_ph2_y[anchitidx], trk_ph2_z[anchitidx]); - // const float phi = hitA.phi(); - // const float eta = hitA.eta(); - } - } - - for (unsigned int idx = 0; idx < nTrackCandidates; idx++) { - std::vector LSs = getLSsFromTC(event, idx); - std::vector lsIdx; - for (auto& LS : LSs) { - lsIdx.push_back(sg_index_map[LS]); - } - ana.tx->pushbackToBranch>("tc_lsIdx", lsIdx); + // Compute pt eta phi and hit indices that will be used to figure out whether the TC matched + float pt, eta, phi; + std::vector hit_idx, hit_type; + switch (type) { + case LSTObjType::pT5: + std::tie(pt, eta, phi, hit_idx, hit_type) = parsepT5(event, idx); + break; + case LSTObjType::pT3: + std::tie(pt, eta, phi, hit_idx, hit_type) = parsepT3(event, idx); + break; + case LSTObjType::T5: + std::tie(pt, eta, phi, hit_idx, hit_type) = parseT5(event, idx, trk_ph2_x, trk_ph2_y, trk_ph2_z); + break; + case LSTObjType::pLS: + std::tie(pt, eta, phi, hit_idx, hit_type) = parsepLS(event, idx); + break; } - std::cout << " mds_used_in_sg.size(): " << mds_used_in_sg.size() << std::endl; -} - -//________________________________________________________________________________________________________________________________ -void setGnnNtupleMiniDoublet(LSTEvent* event, - unsigned int MD, - std::vector const& trk_sim_q, - std::vector const& trk_sim_pt, - std::vector const& trk_sim_eta, - std::vector const& trk_sim_bunchCrossing, - std::vector const& trk_sim_event, - std::vector const& trk_sim_parentVtxIdx, - std::vector const& trk_simvtx_x, - std::vector const& trk_simvtx_y, - std::vector const& trk_simvtx_z, - std::vector const& trk_simhit_simTrkIdx, - std::vector> const& trk_ph2_simHitIdx, - std::vector> const& trk_pix_simHitIdx) { - // Get relevant information - MiniDoubletsConst miniDoublets = event->getMiniDoublets(); - auto hitsBase = event->getInput(); - - // Get the hit indices - unsigned int hit0 = miniDoublets.anchorHitIndices()[MD]; - unsigned int hit1 = miniDoublets.outerHitIndices()[MD]; - - // Get the hit infos - const float hit0_x = hitsBase.xs()[hit0]; - const float hit0_y = hitsBase.ys()[hit0]; - const float hit0_z = hitsBase.zs()[hit0]; - const float hit0_r = sqrt(hit0_x * hit0_x + hit0_y * hit0_y); - const float hit1_x = hitsBase.xs()[hit1]; - const float hit1_y = hitsBase.ys()[hit1]; - const float hit1_z = hitsBase.zs()[hit1]; - const float hit1_r = sqrt(hit1_x * hit1_x + hit1_y * hit1_y); - - // Do sim matching - std::vector hit_idx = {hitsBase.idxs()[hit0], hitsBase.idxs()[hit1]}; - std::vector hit_type = {4, 4}; - std::vector simidxs = - matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); - - bool isFake = simidxs.size() == 0; - int tp_type = getDenomSimTrkType(simidxs, - trk_sim_q, - trk_sim_pt, - trk_sim_eta, - trk_sim_bunchCrossing, - trk_sim_event, - trk_sim_parentVtxIdx, - trk_simvtx_x, - trk_simvtx_y, - trk_simvtx_z); - - auto const& trk_ph2_subdet = trk.getVUS("ph2_subdet"); - auto const& trk_ph2_layer = trk.getVUS("ph2_layer"); - auto const& trk_ph2_detId = trk.getVU("ph2_detId"); - auto const& trk_ph2_x = trk.getVF("ph2_x"); - auto const& trk_ph2_y = trk.getVF("ph2_y"); - auto const& trk_ph2_z = trk.getVF("ph2_z"); - - // Obtain where the actual hit is located in terms of their layer, module, rod, and ring number - unsigned int anchitidx = hitsBase.idxs()[hit0]; - int subdet = trk_ph2_subdet[hitsBase.idxs()[anchitidx]]; - int is_endcap = subdet == 4; - // this accounting makes it so that you have layer 1 2 3 4 5 6 in the barrel, and 7 8 9 10 11 in the endcap. (becuase endcap is ph2_subdet == 4) - int layer = trk_ph2_layer[anchitidx] + 6 * (is_endcap); - int detId = trk_ph2_detId[anchitidx]; - - // Obtaining dPhiChange - float dphichange = miniDoublets.dphichanges()[MD]; - - // Computing pt - float pt = hit0_r * k2Rinv1GeVf / sin(dphichange); + // Perform matching + std::vector simidx = matchedSimTrkIdxs( + hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, matchfrac); + int isFake = simidx.size() == 0; - // T5 eta and phi are computed using outer and innermost hits - lst_math::Hit hitA(trk_ph2_x[anchitidx], trk_ph2_y[anchitidx], trk_ph2_z[anchitidx]); - const float phi = hitA.phi(); - const float eta = hitA.eta(); - - // Mini Doublets - ana.tx->pushbackToBranch("MD_pt", pt); - ana.tx->pushbackToBranch("MD_eta", eta); - ana.tx->pushbackToBranch("MD_phi", phi); - ana.tx->pushbackToBranch("MD_dphichange", dphichange); - ana.tx->pushbackToBranch("MD_isFake", isFake); - ana.tx->pushbackToBranch("MD_tpType", tp_type); - ana.tx->pushbackToBranch("MD_detId", detId); - ana.tx->pushbackToBranch("MD_layer", layer); - ana.tx->pushbackToBranch("MD_0_r", hit0_r); - ana.tx->pushbackToBranch("MD_0_x", hit0_x); - ana.tx->pushbackToBranch("MD_0_y", hit0_y); - ana.tx->pushbackToBranch("MD_0_z", hit0_z); - ana.tx->pushbackToBranch("MD_1_r", hit1_r); - ana.tx->pushbackToBranch("MD_1_x", hit1_x); - ana.tx->pushbackToBranch("MD_1_y", hit1_y); - ana.tx->pushbackToBranch("MD_1_z", hit1_z); - // ana.tx->pushbackToBranch("MD_sim_idx", simidxs.size() > 0 ? simidxs[0] : -999); + return {type, pt, eta, phi, isFake, simidx}; } //________________________________________________________________________________________________________________________________ -std::tuple> parseTrackCandidate( +std::tuple, std::vector> parseTrackCandidateAllMatch( LSTEvent* event, unsigned int idx, std::vector const& trk_ph2_x, @@ -1346,7 +2407,8 @@ std::tuple> parseTrackCandidate( std::vector const& trk_ph2_z, std::vector const& trk_simhit_simTrkIdx, std::vector> const& trk_ph2_simHitIdx, - std::vector> const& trk_pix_simHitIdx) { + std::vector> const& trk_pix_simHitIdx, + float matchfrac) { // Get the type of the track candidate auto const& trackCandidatesBase = event->getTrackCandidatesBase(); short type = trackCandidatesBase.trackCandidateType()[idx]; @@ -1370,11 +2432,13 @@ std::tuple> parseTrackCandidate( } // Perform matching - std::vector simidx = - matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); + std::vector simidx; + std::vector simidxfrac; + std::tie(simidx, simidxfrac) = matchedSimTrkIdxsAndFracs( + hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, matchfrac); int isFake = simidx.size() == 0; - return {type, pt, eta, phi, isFake, simidx}; + return {type, pt, eta, phi, isFake, simidx, simidxfrac}; } //________________________________________________________________________________________________________________________________ @@ -1583,84 +2647,6 @@ std::tuple, std::vectorgetPixelSegments(); - auto const& pixelSeeds = event->getInput(); - int n_accepted_simtrk = ana.tx->getBranch>("sim_TC_matched").size(); - - auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx"); - auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx"); - auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx"); - - unsigned int n_pLS = pixelSegments.metadata().size(); - std::vector sim_pLS_matched(n_accepted_simtrk, 0); - std::vector> pLS_matched_simIdx; - - for (unsigned int i_pLS = 0; i_pLS < n_pLS; ++i_pLS) { - // Get pLS properties - float pt = pixelSeeds.ptIn()[i_pLS]; - float px = pixelSeeds.px()[i_pLS]; - float py = pixelSeeds.py()[i_pLS]; - float pz = pixelSeeds.pz()[i_pLS]; - bool isQuad = static_cast(pixelSeeds.isQuad()[i_pLS]); - float ptErr = pixelSeeds.ptErr()[i_pLS]; - float eta = pixelSeeds.eta()[i_pLS]; - float etaErr = pixelSeeds.etaErr()[i_pLS]; - float phi = pixelSeeds.phi()[i_pLS]; - float score = pixelSegments.score()[i_pLS]; - float centerX = pixelSegments.circleCenterX()[i_pLS]; - float centerY = pixelSegments.circleCenterY()[i_pLS]; - float radius = pixelSegments.circleRadius()[i_pLS]; - - // Get hits from pLS - std::vector hit_idx = getPixelHitIdxsFrompLS(event, i_pLS); - std::vector hit_type = getPixelHitTypesFrompLS(event, i_pLS); - - // Match to sim tracks - std::vector simidx = - matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx); - bool isFake = simidx.empty(); - - // Fill branches - ana.tx->pushbackToBranch("pLS_ptIn", pt); - ana.tx->pushbackToBranch("pLS_ptErr", ptErr); - ana.tx->pushbackToBranch("pLS_px", px); - ana.tx->pushbackToBranch("pLS_py", py); - ana.tx->pushbackToBranch("pLS_pz", pz); - ana.tx->pushbackToBranch("pLS_eta", eta); - ana.tx->pushbackToBranch("pLS_etaErr", etaErr); - ana.tx->pushbackToBranch("pLS_phi", phi); - ana.tx->pushbackToBranch("pLS_score", score); - ana.tx->pushbackToBranch("pLS_circleCenterX", centerX); - ana.tx->pushbackToBranch("pLS_circleCenterY", centerY); - ana.tx->pushbackToBranch("pLS_circleRadius", radius); - ana.tx->pushbackToBranch("pLS_isQuad", isQuad); - ana.tx->pushbackToBranch("pLS_isFake", isFake); - pLS_matched_simIdx.push_back(simidx); - - // Count matches - for (auto& idx : simidx) { - if (idx < n_accepted_simtrk) { - sim_pLS_matched[idx]++; - } - } - } - - std::vector pLS_isDuplicate(pLS_matched_simIdx.size(), 0); - for (size_t i = 0; i < pLS_matched_simIdx.size(); ++i) { - for (int simidx : pLS_matched_simIdx[i]) { - if (simidx < n_accepted_simtrk && sim_pLS_matched[simidx] > 1) { - pLS_isDuplicate[i] = 1; - break; - } - } - } - - ana.tx->setBranch>("sim_pLS_matched", sim_pLS_matched); - ana.tx->setBranch>>("pLS_matched_simIdx", pLS_matched_simIdx); - ana.tx->setBranch>("pLS_isDuplicate", pLS_isDuplicate); -} - //________________________________________________________________________________________________________________________________ void printHitMultiplicities(LSTEvent* event) { auto modules = event->getModules(); diff --git a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h index 64327f970b1ee..83c591bc7c5e6 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h +++ b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h @@ -15,42 +15,66 @@ using LSTEvent = ALPAKA_ACCELERATOR_NAMESPACE::lst::LSTEvent; // Common void createOutputBranches(); -void createRequiredOutputBranches(); -void createOptionalOutputBranches(); -void createGnnNtupleBranches(); +void createJetBranches(); void createT5DNNBranches(); void createT3DNNBranches(); void createpT3DNNBranches(); +void createSimTrackContainerBranches(); +void createTrackCandidateBranches(); +void createMiniDoubletBranches(); +void createLineSegmentBranches(); +void createTripletBranches(); +void createQuintupletBranches(); +void createPixelLineSegmentBranches(); +void createPixelTripletBranches(); +void createPixelQuintupletBranches(); +void createOccupancyBranches(); void fillOutputBranches(LSTEvent* event); -void setOutputBranches(LSTEvent* event); -void setOptionalOutputBranches(LSTEvent* event); void setOccupancyBranches(LSTEvent* event); -void setPixelQuintupletOutputBranches(LSTEvent* event); -void setQuintupletOutputBranches(LSTEvent* event); -void setPixelTripletOutputBranches(LSTEvent* event); -void setGnnNtupleBranches(LSTEvent* event); -void setGnnNtupleMiniDoublet(LSTEvent* event, - unsigned int MD, - std::vector const& trk_sim_q, - std::vector const& trk_sim_pt, - std::vector const& trk_sim_eta, - std::vector const& trk_sim_bunchCrossing, - std::vector const& trk_sim_event, - std::vector const& trk_sim_parentVtxIdx, - std::vector const& trk_simvtx_x, - std::vector const& trk_simvtx_y, - std::vector const& trk_simvtx_z, - std::vector const& trk_simhit_simTrkIdx, - std::vector> const& trk_ph2_simHitIdx, - std::vector> const& trk_pix_simHitIdx); +unsigned int setSimTrackContainerBranches(LSTEvent* event); +void setTrackCandidateBranches(LSTEvent* event, + unsigned int n_accepted_tracks, + std::map t5_idx_map, + std::map pls_idx_map, + std::map pt3_idx_map, + std::map pt5_idx_map, + float matchfrac); +std::map setMiniDoubletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac); +std::map setLineSegmentBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& md_idx_map); +std::map setTripletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& ls_idx_map); +std::map setQuintupletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& t3_idx_map); +std::map setPixelLineSegmentBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac); +std::map setPixelTripletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& pls_idx_map, + std::map const& t3_idx_map); +std::map setPixelQuintupletBranches(LSTEvent* event, + unsigned int n_accepted_simtrk, + float matchfrac, + std::map const& pls_idx_map, + std::map const& t5_idx_map); + void fillT5DNNBranches(LSTEvent* event, unsigned int T3); void fillT3DNNBranches(LSTEvent* event, unsigned int iT3); void fillpT3DNNBranches(LSTEvent* event, unsigned int iPT3); void setT5DNNBranches(LSTEvent* event); -void setT3DNNBranches(LSTEvent* event); +void setT3DNNBranches(LSTEvent* event, float matchfrac = 0.75); void setpT3DNNBranches(LSTEvent* event); -void setpLSOutputBranches(LSTEvent* event); std::tuple> parseTrackCandidate( LSTEvent* event, @@ -60,7 +84,18 @@ std::tuple> parseTrackCandidate( std::vector const& trk_ph2_z, std::vector const& trk_simhit_simTrkIdx, std::vector> const& trk_ph2_simHitIdx, - std::vector> const& trk_pix_simHitIdx); + std::vector> const& trk_pix_simHitIdx, + float matchfrac = 0.75); +std::tuple, std::vector> parseTrackCandidateAllMatch( + LSTEvent* event, + unsigned int, + std::vector const& trk_ph2_x, + std::vector const& trk_ph2_y, + std::vector const& trk_ph2_z, + std::vector const& trk_simhit_simTrkIdx, + std::vector> const& trk_ph2_simHitIdx, + std::vector> const& trk_pix_simHitIdx, + float matchfrac = 0.75); std::tuple, std::vector> parsepT5(LSTEvent* event, unsigned int); std::tuple, std::vector> parsepT3(LSTEvent* event, diff --git a/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.cc b/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.cc index 5e18067733bfd..f10faa486bdf7 100644 --- a/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.cc +++ b/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.cc @@ -15,7 +15,7 @@ void TreeUtil::Init(TTree* treeIn) { void TreeUtil::GetEntry(unsigned int idx) { index = idx; for (auto& pair : data) { - pair.second.isLoaded = false; + pair.second->isLoaded = false; } } @@ -82,46 +82,52 @@ void TreeUtil::loadAllBranches() { } template -T const& TreeUtil::get(std::string name) { - auto search = data.find(name); - if (search == data.end()) { +const T& TreeUtil::get(const std::string& name) { + auto it = data.find(name); + if (it == data.end()) { tree->SetMakeClass(1); - TBranch* branch = tree->GetBranch(name.c_str()); - if (branch == nullptr) + TBranch* br = tree->GetBranch(name.c_str()); + if (!br) throw std::out_of_range("Branch " + name + " does not exist!"); - search = data.emplace(name, std::move(BranchData{branch, false, static_cast(nullptr)})).first; - branch->SetAddress(&search->second.ptr); + it = data.emplace(name, std::make_unique>(br)).first; tree->SetMakeClass(0); } - if (!search->second.isLoaded) { - search->second.branch->GetEntry(index); - search->second.isLoaded = true; + + auto* bd = static_cast*>(it->second.get()); + if (!bd->isLoaded) { + bd->branch->GetEntry(index); + bd->adoptIfNeeded(); + bd->isLoaded = true; } - return *std::get(search->second.ptr); + return *static_cast(bd->getRaw()); } -short const& TreeUtil::getS(std::string name) { return get(name); } -unsigned short const& TreeUtil::getUS(std::string name) { return get(name); } -int const& TreeUtil::getI(std::string name) { return get(name); } -unsigned int const& TreeUtil::getU(std::string name) { return get(name); } -float const& TreeUtil::getF(std::string name) { return get(name); } -std::vector const& TreeUtil::getVS(std::string name) { return get>(name); } -std::vector const& TreeUtil::getVUS(std::string name) { return get>(name); } -std::vector const& TreeUtil::getVI(std::string name) { return get>(name); } -std::vector const& TreeUtil::getVU(std::string name) { return get>(name); } -std::vector const& TreeUtil::getVF(std::string name) { return get>(name); } -std::vector> const& TreeUtil::getVVS(std::string name) { +short const& TreeUtil::getS(std::string const& name) { return get(name); } +unsigned short const& TreeUtil::getUS(std::string const& name) { return get(name); } +int const& TreeUtil::getI(std::string const& name) { return get(name); } +unsigned int const& TreeUtil::getU(std::string const& name) { return get(name); } +float const& TreeUtil::getF(std::string const& name) { return get(name); } +std::vector const& TreeUtil::getVS(std::string const& name) { return get>(name); } +std::vector const& TreeUtil::getVUS(std::string const& name) { + return get>(name); +} +std::vector const& TreeUtil::getVI(std::string const& name) { return get>(name); } +std::vector const& TreeUtil::getVU(std::string const& name) { + return get>(name); +} +std::vector const& TreeUtil::getVF(std::string const& name) { return get>(name); } +std::vector> const& TreeUtil::getVVS(std::string const& name) { return get>>(name); } -std::vector> const& TreeUtil::getVVUS(std::string name) { +std::vector> const& TreeUtil::getVVUS(std::string const& name) { return get>>(name); } -std::vector> const& TreeUtil::getVVI(std::string name) { +std::vector> const& TreeUtil::getVVI(std::string const& name) { return get>>(name); } -std::vector> const& TreeUtil::getVVU(std::string name) { +std::vector> const& TreeUtil::getVVU(std::string const& name) { return get>>(name); } -std::vector> const& TreeUtil::getVVF(std::string name) { +std::vector> const& TreeUtil::getVVF(std::string const& name) { return get>>(name); } diff --git a/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.h b/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.h index 1e1e65058dc8d..d4cdae1222c3b 100644 --- a/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.h +++ b/RecoTracker/LSTCore/standalone/code/rooutil/treeutil.h @@ -11,31 +11,51 @@ class TreeUtil { private: - using DataTypes = std::variant*, - std::vector*, - std::vector*, - std::vector*, - std::vector*, - std::vector>*, - std::vector>*, - std::vector>*, - std::vector>*, - std::vector>*>; - - struct BranchData { + struct BranchDataBase { TBranch* branch; - bool isLoaded; - DataTypes ptr; + bool isLoaded = false; + explicit BranchDataBase(TBranch* b) : branch(b) {} + virtual ~BranchDataBase() = default; + virtual const void* getRaw() const = 0; + virtual void adoptIfNeeded() = 0; + }; + + template + struct BranchDataHolder : BranchDataBase { + std::unique_ptr buffer; + mutable T* raw = nullptr; + mutable std::unique_ptr owner; + + explicit BranchDataHolder(TBranch* b) : BranchDataBase(b) { + if constexpr (std::is_fundamental_v) { + buffer = std::make_unique(); + branch->SetAddress(buffer.get()); + } else { + branch->SetAddress(&raw); + branch->SetAutoDelete(false); + } + } + + const void* getRaw() const override { + if constexpr (std::is_fundamental_v) { + return buffer.get(); + } else { + return owner ? static_cast(owner.get()) : static_cast(raw); + } + } + + void adoptIfNeeded() override { + if constexpr (!std::is_fundamental_v) { + if (!owner && raw) { + owner.reset(raw); + } + } + } }; TTree* tree; unsigned int index; - std::unordered_map data; + std::unordered_map> data; public: void Init(TTree* tree); @@ -43,22 +63,22 @@ class TreeUtil { static void progress(int nEventsTotal, int nEventsChain); void loadAllBranches(); template - T const& get(std::string name); - short const& getS(std::string name); - unsigned short const& getUS(std::string name); - int const& getI(std::string name); - unsigned int const& getU(std::string name); - float const& getF(std::string name); - std::vector const& getVS(std::string name); - std::vector const& getVUS(std::string name); - std::vector const& getVI(std::string name); - std::vector const& getVU(std::string name); - std::vector const& getVF(std::string name); - std::vector> const& getVVS(std::string name); - std::vector> const& getVVUS(std::string name); - std::vector> const& getVVI(std::string name); - std::vector> const& getVVU(std::string name); - std::vector> const& getVVF(std::string name); + T const& get(std::string const& name); + short const& getS(std::string const& name); + unsigned short const& getUS(std::string const& name); + int const& getI(std::string const& name); + unsigned int const& getU(std::string const& name); + float const& getF(std::string const& name); + std::vector const& getVS(std::string const& name); + std::vector const& getVUS(std::string const& name); + std::vector const& getVI(std::string const& name); + std::vector const& getVU(std::string const& name); + std::vector const& getVF(std::string const& name); + std::vector> const& getVVS(std::string const& name); + std::vector> const& getVVUS(std::string const& name); + std::vector> const& getVVI(std::string const& name); + std::vector> const& getVVU(std::string const& name); + std::vector> const& getVVF(std::string const& name); }; #endif diff --git a/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py b/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py index 3c9064b52d8c6..c422cd79c45c1 100755 --- a/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py +++ b/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py @@ -10,7 +10,7 @@ sel_choices = ["base", "loweta", "xtr", "vtr", "none"] metric_choices = ["eff", "fakerate", "duplrate"] variable_choices = ["pt", "ptmtv", "ptlow", "eta", "phi", "dxy", "dz", "vxy", "deltaEta", "deltaPhi", "deltaR", "jet_eta", "jet_phi", "jet_pt"] -objecttype_choices = ["TC", "pT5", "T5", "pT3", "pLS", "pT5_lower", "pT3_lower", "T5_lower"] +objecttype_choices = ["TC", "pT5", "T5", "pT3", "pLS", "pT5_lower", "pT3_lower", "T5_lower", "pLS_lower"] #lowerObjectType = ["pT5_lower", "pT3_lower", "T5_lower"] r.gROOT.SetBatch(True) @@ -696,6 +696,7 @@ def plot_standard_performance_plots(args): "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], + "pLS_lower":[False], }, "fakerate":{ "TC": [True, False], @@ -706,6 +707,7 @@ def plot_standard_performance_plots(args): "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], + "pLS_lower":[False], }, "duplrate":{ "TC": [True, False], @@ -716,6 +718,7 @@ def plot_standard_performance_plots(args): "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], + "pLS_lower":[False], }, } pdgids = { @@ -764,16 +767,16 @@ def plot_standard_performance_plots(args): if args.individual: # Only eff / TC matters here - breakdowns = {"eff":{"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False]}, - "fakerate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False]}, - "duplrate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False]}} + breakdowns = {"eff":{"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}, + "fakerate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}, + "duplrate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}} else: # Only eff / TC matters here - breakdowns = {"eff":{"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False]}, - "fakerate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False]}, - "duplrate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False]}} + breakdowns = {"eff":{"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}, + "fakerate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}, + "duplrate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}} if args.yzoom: args.yzooms = [args.yzoom] diff --git a/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc b/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc index 0bac1d51ad7cf..ed625a8807031 100644 --- a/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc +++ b/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc @@ -46,8 +46,7 @@ int main(int argc, char** argv) { TString("TC_") + selnames[isel], /* pdgid */ pdgid, /* q */ charge, - /* pass */ - [&](unsigned int isim) { return lstEff.getVI("sim_TC_matched").at(isim) > 0; }, + /* pass */ [&](unsigned int isim) { return lstEff.getVI("sim_tcIdx").at(isim) >= 0; }, /* sel */ sels[isel])); list_effSetDef.push_back(SimTrackSetDefinition( /* name */ @@ -55,7 +54,11 @@ int main(int argc, char** argv) { /* pdgid */ pdgid, /* q */ charge, /* pass */ - [&](unsigned int isim) { return lstEff.getVI("sim_TC_matched_mask").at(isim) & (1 << pT5); }, + [&](unsigned int isim) { + auto& lstEff_sim_tcIdx = lstEff.getVI("sim_tcIdx"); + auto& lstEff_tc_type = lstEff.getVI("tc_type"); + return lstEff_sim_tcIdx.at(isim) >= 0 ? lstEff_tc_type.at(lstEff_sim_tcIdx.at(isim)) == pT5 : false; + }, /* sel */ sels[isel])); list_effSetDef.push_back(SimTrackSetDefinition( /* name */ @@ -63,7 +66,11 @@ int main(int argc, char** argv) { /* pdgid */ pdgid, /* q */ charge, /* pass */ - [&](unsigned int isim) { return lstEff.getVI("sim_TC_matched_mask").at(isim) & (1 << pT3); }, + [&](unsigned int isim) { + auto& lstEff_sim_tcIdx = lstEff.getVI("sim_tcIdx"); + auto& lstEff_tc_type = lstEff.getVI("tc_type"); + return lstEff_sim_tcIdx.at(isim) >= 0 ? lstEff_tc_type.at(lstEff_sim_tcIdx.at(isim)) == pT3 : false; + }, /* sel */ sels[isel])); list_effSetDef.push_back(SimTrackSetDefinition( /* name */ @@ -71,7 +78,11 @@ int main(int argc, char** argv) { /* pdgid */ pdgid, /* q */ charge, /* pass */ - [&](unsigned int isim) { return lstEff.getVI("sim_TC_matched_mask").at(isim) & (1 << T5); }, + [&](unsigned int isim) { + auto& lstEff_sim_tcIdx = lstEff.getVI("sim_tcIdx"); + auto& lstEff_tc_type = lstEff.getVI("tc_type"); + return lstEff_sim_tcIdx.at(isim) >= 0 ? lstEff_tc_type.at(lstEff_sim_tcIdx.at(isim)) == T5 : false; + }, /* sel */ sels[isel])); list_effSetDef.push_back(SimTrackSetDefinition( /* name */ @@ -79,7 +90,11 @@ int main(int argc, char** argv) { /* pdgid */ pdgid, /* q */ charge, /* pass */ - [&](unsigned int isim) { return lstEff.getVI("sim_TC_matched_mask").at(isim) & (1 << pLS); }, + [&](unsigned int isim) { + auto& lstEff_sim_tcIdx = lstEff.getVI("sim_tcIdx"); + auto& lstEff_tc_type = lstEff.getVI("tc_type"); + return lstEff_sim_tcIdx.at(isim) >= 0 ? lstEff_tc_type.at(lstEff_sim_tcIdx.at(isim)) == pLS : false; + }, /* sel */ sels[isel])); if (ana.do_lower_level) { @@ -89,7 +104,15 @@ int main(int argc, char** argv) { TString("pT5_lower_") + selnames[isel], /* pdgid */ pdgid, /* q */ charge, - /* pass */ [&](unsigned int isim) { return lstEff.getVI("sim_pT5_matched").at(isim) > 0; }, + /* pass */ + [&](unsigned int isim) { + auto& lstEff_sim_pt5IdxAllFrac = lstEff.getVVF("sim_pt5IdxAllFrac"); + for (size_t i = 0; i < lstEff_sim_pt5IdxAllFrac.at(isim).size(); ++i) { + if (lstEff_sim_pt5IdxAllFrac.at(isim).at(i) > 0.75) + return true; + } + return false; + }, /* sel */ sels[isel])); list_effSetDef.push_back( SimTrackSetDefinition(/* name */ @@ -97,14 +120,44 @@ int main(int argc, char** argv) { /* pdgid */ pdgid, /* q */ charge, /* pass */ - [&](unsigned int isim) { return lstEff.getVI("sim_T5_matched").at(isim) > 0; }, + [&](unsigned int isim) { + auto& lstEff_sim_t5IdxAllFrac = lstEff.getVVF("sim_t5IdxAllFrac"); + for (size_t i = 0; i < lstEff_sim_t5IdxAllFrac.at(isim).size(); ++i) { + if (lstEff_sim_t5IdxAllFrac.at(isim).at(i) > 0.75) + return true; + } + return false; + }, /* sel */ sels[isel])); list_effSetDef.push_back(SimTrackSetDefinition( /* name */ TString("pT3_lower_") + selnames[isel], /* pdgid */ pdgid, /* q */ charge, - /* pass */ [&](unsigned int isim) { return lstEff.getVI("sim_pT3_matched").at(isim) > 0; }, + /* pass */ + [&](unsigned int isim) { + auto& lstEff_sim_pt3IdxAllFrac = lstEff.getVVF("sim_pt3IdxAllFrac"); + for (size_t i = 0; i < lstEff_sim_pt3IdxAllFrac.at(isim).size(); ++i) { + if (lstEff_sim_pt3IdxAllFrac.at(isim).at(i) > 0.75) + return true; + } + return false; + }, + /* sel */ sels[isel])); + list_effSetDef.push_back(SimTrackSetDefinition( + /* name */ + TString("pLS_lower_") + selnames[isel], + /* pdgid */ pdgid, + /* q */ charge, + /* pass */ + [&](unsigned int isim) { + auto& lstEff_sim_plsIdxAllFrac = lstEff.getVVF("sim_plsIdxAllFrac"); + for (size_t i = 0; i < lstEff_sim_plsIdxAllFrac.at(isim).size(); ++i) { + if (lstEff_sim_plsIdxAllFrac.at(isim).at(i) > 0.75) + return true; + } + return false; + }, /* sel */ sels[isel])); } } @@ -198,6 +251,15 @@ int main(int argc, char** argv) { /* eta */ [&]() { return lstEff.getVF("pT3_eta"); }, /* phi */ [&]() { return lstEff.getVF("pT3_phi"); }, /* type */ [&]() { return std::vector(lstEff.getVF("pT3_pt").size(), 1); })); + list_FRSetDef.push_back(RecoTrackSetDefinition( + /* name */ + "pLS_lower", + /* pass */ [&](unsigned int ipLS) { return lstEff.getVI("pLS_isFake").at(ipLS) > 0; }, + /* sel */ [&](unsigned int ipLS) { return 1; }, + /* pt */ [&]() { return lstEff.getVF("pLS_pt"); }, + /* eta */ [&]() { return lstEff.getVF("pLS_eta"); }, + /* phi */ [&]() { return lstEff.getVF("pLS_phi"); }, + /* type */ [&]() { return std::vector(lstEff.getVF("pLS_pt").size(), 1); })); } bookFakeRateSets(list_FRSetDef); @@ -287,6 +349,15 @@ int main(int argc, char** argv) { /* eta */ [&]() { return lstEff.getVF("pT3_eta"); }, /* phi */ [&]() { return lstEff.getVF("pT3_phi"); }, /* type */ [&]() { return std::vector(lstEff.getVF("pT3_pt").size(), 1); })); + list_DRSetDef.push_back(RecoTrackSetDefinition( + /* name */ + "pLS_lower", + /* pass */ [&](unsigned int ipLS) { return lstEff.getVI("pLS_isDuplicate").at(ipLS) > 0; }, + /* sel */ [&](unsigned int ipLS) { return 1; }, + /* pt */ [&]() { return lstEff.getVF("pLS_pt"); }, + /* eta */ [&]() { return lstEff.getVF("pLS_eta"); }, + /* phi */ [&]() { return lstEff.getVF("pLS_phi"); }, + /* type */ [&]() { return std::vector(lstEff.getVF("pLS_pt").size(), 1); })); } bookDuplicateRateSets(list_DRSetDef);