Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-12 23:30:08 

0001 #include "LSTEvent.h"
0002 #include "Circle.h"
0003 
0004 #include "write_lst_ntuple.h"
0005 
0006 using namespace ALPAKA_ACCELERATOR_NAMESPACE::lst;
0007 
0008 //________________________________________________________________________________________________________________________________
0009 void createOutputBranches() {
0010   createRequiredOutputBranches();
0011   createOptionalOutputBranches();
0012 }
0013 
0014 //________________________________________________________________________________________________________________________________
0015 void fillOutputBranches(LSTEvent* event) {
0016   setOutputBranches(event);
0017   setOptionalOutputBranches(event);
0018   if (ana.gnn_ntuple)
0019     setGnnNtupleBranches(event);
0020 
0021   // Now actually fill the ttree
0022   ana.tx->fill();
0023 
0024   // Then clear the branches to default values (e.g. -999, or clear the vectors to empty vectors)
0025   ana.tx->clear();
0026 }
0027 
0028 //________________________________________________________________________________________________________________________________
0029 void createRequiredOutputBranches() {
0030   // Setup output TTree
0031   ana.tx->createBranch<std::vector<float>>("sim_pt");
0032   ana.tx->createBranch<std::vector<float>>("sim_eta");
0033   ana.tx->createBranch<std::vector<float>>("sim_phi");
0034   ana.tx->createBranch<std::vector<float>>("sim_pca_dxy");
0035   ana.tx->createBranch<std::vector<float>>("sim_pca_dz");
0036   ana.tx->createBranch<std::vector<int>>("sim_q");
0037   ana.tx->createBranch<std::vector<int>>("sim_event");
0038   ana.tx->createBranch<std::vector<int>>("sim_pdgId");
0039   ana.tx->createBranch<std::vector<float>>("sim_vx");
0040   ana.tx->createBranch<std::vector<float>>("sim_vy");
0041   ana.tx->createBranch<std::vector<float>>("sim_vz");
0042   ana.tx->createBranch<std::vector<float>>("sim_trkNtupIdx");
0043   ana.tx->createBranch<std::vector<int>>("sim_TC_matched");
0044   ana.tx->createBranch<std::vector<int>>("sim_TC_matched_mask");
0045 
0046   // Track candidates
0047   ana.tx->createBranch<std::vector<float>>("tc_pt");
0048   ana.tx->createBranch<std::vector<float>>("tc_eta");
0049   ana.tx->createBranch<std::vector<float>>("tc_phi");
0050   ana.tx->createBranch<std::vector<int>>("tc_type");
0051   ana.tx->createBranch<std::vector<int>>("tc_isFake");
0052   ana.tx->createBranch<std::vector<int>>("tc_isDuplicate");
0053   ana.tx->createBranch<std::vector<std::vector<int>>>("tc_matched_simIdx");
0054 }
0055 
0056 //________________________________________________________________________________________________________________________________
0057 void createOptionalOutputBranches() {
0058 #ifdef CUT_VALUE_DEBUG
0059   // Event-wide branches
0060   // ana.tx->createBranch<float>("evt_dummy");
0061 
0062   // Sim Track branches
0063   // NOTE: Must sync with main tc branch in length!!
0064   ana.tx->createBranch<std::vector<float>>("sim_dummy");
0065 
0066   // Track Candidate branches
0067   // NOTE: Must sync with main tc branch in length!!
0068   ana.tx->createBranch<std::vector<float>>("tc_dummy");
0069 
0070   // pT5 branches
0071   ana.tx->createBranch<std::vector<std::vector<int>>>("pT5_matched_simIdx");
0072   ana.tx->createBranch<std::vector<std::vector<int>>>("pT5_hitIdxs");
0073   ana.tx->createBranch<std::vector<int>>("sim_pT5_matched");
0074   ana.tx->createBranch<std::vector<float>>("pT5_pt");
0075   ana.tx->createBranch<std::vector<float>>("pT5_eta");
0076   ana.tx->createBranch<std::vector<float>>("pT5_phi");
0077   ana.tx->createBranch<std::vector<int>>("pT5_isFake");
0078   ana.tx->createBranch<std::vector<float>>("t5_sim_vxy");
0079   ana.tx->createBranch<std::vector<float>>("t5_sim_vz");
0080   ana.tx->createBranch<std::vector<int>>("pT5_isDuplicate");
0081   ana.tx->createBranch<std::vector<int>>("pT5_score");
0082   ana.tx->createBranch<std::vector<int>>("pT5_layer_binary");
0083   ana.tx->createBranch<std::vector<int>>("pT5_moduleType_binary");
0084   ana.tx->createBranch<std::vector<float>>("pT5_matched_pt");
0085   ana.tx->createBranch<std::vector<float>>("pT5_rzChiSquared");
0086   ana.tx->createBranch<std::vector<float>>("pT5_rPhiChiSquared");
0087   ana.tx->createBranch<std::vector<float>>("pT5_rPhiChiSquaredInwards");
0088 
0089   // pT3 branches
0090   ana.tx->createBranch<std::vector<int>>("sim_pT3_matched");
0091   ana.tx->createBranch<std::vector<float>>("pT3_pt");
0092   ana.tx->createBranch<std::vector<int>>("pT3_isFake");
0093   ana.tx->createBranch<std::vector<int>>("pT3_isDuplicate");
0094   ana.tx->createBranch<std::vector<float>>("pT3_eta");
0095   ana.tx->createBranch<std::vector<float>>("pT3_phi");
0096   ana.tx->createBranch<std::vector<float>>("pT3_score");
0097   ana.tx->createBranch<std::vector<int>>("pT3_foundDuplicate");
0098   ana.tx->createBranch<std::vector<std::vector<int>>>("pT3_matched_simIdx");
0099   ana.tx->createBranch<std::vector<std::vector<int>>>("pT3_hitIdxs");
0100   ana.tx->createBranch<std::vector<float>>("pT3_pixelRadius");
0101   ana.tx->createBranch<std::vector<float>>("pT3_pixelRadiusError");
0102   ana.tx->createBranch<std::vector<std::vector<float>>>("pT3_matched_pt");
0103   ana.tx->createBranch<std::vector<float>>("pT3_tripletRadius");
0104   ana.tx->createBranch<std::vector<float>>("pT3_rPhiChiSquared");
0105   ana.tx->createBranch<std::vector<float>>("pT3_rPhiChiSquaredInwards");
0106   ana.tx->createBranch<std::vector<float>>("pT3_rzChiSquared");
0107   ana.tx->createBranch<std::vector<int>>("pT3_layer_binary");
0108   ana.tx->createBranch<std::vector<int>>("pT3_moduleType_binary");
0109 
0110   // pLS branches
0111   ana.tx->createBranch<std::vector<int>>("sim_pLS_matched");
0112   ana.tx->createBranch<std::vector<std::vector<int>>>("pLS_matched_simIdx");
0113   ana.tx->createBranch<std::vector<int>>("pLS_isFake");
0114   ana.tx->createBranch<std::vector<int>>("pLS_isDuplicate");
0115   ana.tx->createBranch<std::vector<float>>("pLS_ptIn");
0116   ana.tx->createBranch<std::vector<float>>("pLS_ptErr");
0117   ana.tx->createBranch<std::vector<float>>("pLS_px");
0118   ana.tx->createBranch<std::vector<float>>("pLS_py");
0119   ana.tx->createBranch<std::vector<float>>("pLS_pz");
0120   ana.tx->createBranch<std::vector<float>>("pLS_eta");
0121   ana.tx->createBranch<std::vector<bool>>("pLS_isQuad");
0122   ana.tx->createBranch<std::vector<float>>("pLS_etaErr");
0123   ana.tx->createBranch<std::vector<float>>("pLS_phi");
0124   ana.tx->createBranch<std::vector<float>>("pLS_score");
0125   ana.tx->createBranch<std::vector<float>>("pLS_circleCenterX");
0126   ana.tx->createBranch<std::vector<float>>("pLS_circleCenterY");
0127   ana.tx->createBranch<std::vector<float>>("pLS_circleRadius");
0128 
0129   // T5 branches
0130   ana.tx->createBranch<std::vector<int>>("sim_T5_matched");
0131   ana.tx->createBranch<std::vector<int>>("t5_isFake");
0132   ana.tx->createBranch<std::vector<int>>("t5_isDuplicate");
0133   ana.tx->createBranch<std::vector<int>>("t5_foundDuplicate");
0134   ana.tx->createBranch<std::vector<float>>("t5_pt");
0135   ana.tx->createBranch<std::vector<float>>("t5_pMatched");
0136   ana.tx->createBranch<std::vector<float>>("t5_eta");
0137   ana.tx->createBranch<std::vector<float>>("t5_phi");
0138   ana.tx->createBranch<std::vector<float>>("t5_score_rphisum");
0139   ana.tx->createBranch<std::vector<std::vector<int>>>("t5_hitIdxs");
0140   ana.tx->createBranch<std::vector<std::vector<int>>>("t5_matched_simIdx");
0141   ana.tx->createBranch<std::vector<int>>("t5_moduleType_binary");
0142   ana.tx->createBranch<std::vector<int>>("t5_layer_binary");
0143   ana.tx->createBranch<std::vector<float>>("t5_matched_pt");
0144   ana.tx->createBranch<std::vector<float>>("t5_innerRadius");
0145   ana.tx->createBranch<std::vector<float>>("t5_outerRadius");
0146   ana.tx->createBranch<std::vector<float>>("t5_bridgeRadius");
0147   ana.tx->createBranch<std::vector<float>>("t5_chiSquared");
0148   ana.tx->createBranch<std::vector<float>>("t5_rzChiSquared");
0149   ana.tx->createBranch<std::vector<int>>("t5_isDupAlgoFlag");
0150   ana.tx->createBranch<std::vector<float>>("t5_nonAnchorChiSquared");
0151   ana.tx->createBranch<std::vector<float>>("t5_dBeta1");
0152   ana.tx->createBranch<std::vector<float>>("t5_dBeta2");
0153 
0154   // Occupancy branches
0155   ana.tx->createBranch<std::vector<int>>("module_layers");
0156   ana.tx->createBranch<std::vector<int>>("module_subdets");
0157   ana.tx->createBranch<std::vector<int>>("module_rings");
0158   ana.tx->createBranch<std::vector<int>>("module_rods");
0159   ana.tx->createBranch<std::vector<int>>("module_modules");
0160   ana.tx->createBranch<std::vector<bool>>("module_isTilted");
0161   ana.tx->createBranch<std::vector<float>>("module_eta");
0162   ana.tx->createBranch<std::vector<float>>("module_r");
0163   ana.tx->createBranch<std::vector<int>>("md_occupancies");
0164   ana.tx->createBranch<std::vector<int>>("sg_occupancies");
0165   ana.tx->createBranch<std::vector<int>>("t3_occupancies");
0166   ana.tx->createBranch<int>("tc_occupancies");
0167   ana.tx->createBranch<std::vector<int>>("t5_occupancies");
0168   ana.tx->createBranch<int>("pT3_occupancies");
0169   ana.tx->createBranch<int>("pT5_occupancies");
0170 
0171   // T5 DNN branches
0172   createT5DNNBranches();
0173   createT3DNNBranches();
0174 
0175 #endif
0176 }
0177 
0178 //________________________________________________________________________________________________________________________________
0179 void createT5DNNBranches() {
0180   // Common branches
0181   ana.tx->createBranch<std::vector<int>>("t5_t3_idx0");
0182   ana.tx->createBranch<std::vector<int>>("t5_t3_idx1");
0183   ana.tx->createBranch<std::vector<int>>("t5_tc_idx");
0184   ana.tx->createBranch<std::vector<int>>("t5_partOfTC");
0185   ana.tx->createBranch<std::vector<float>>("t5_t3_pt");
0186   ana.tx->createBranch<std::vector<float>>("t5_t3_eta");
0187   ana.tx->createBranch<std::vector<float>>("t5_t3_phi");
0188   ana.tx->createBranch<std::vector<float>>("t5_t3_fakeScore1");
0189   ana.tx->createBranch<std::vector<float>>("t5_t3_promptScore1");
0190   ana.tx->createBranch<std::vector<float>>("t5_t3_displacedScore1");
0191   ana.tx->createBranch<std::vector<float>>("t5_t3_fakeScore2");
0192   ana.tx->createBranch<std::vector<float>>("t5_t3_promptScore2");
0193   ana.tx->createBranch<std::vector<float>>("t5_t3_displacedScore2");
0194 
0195   // Hit-specific branches
0196   std::vector<std::string> hitIndices = {"0", "1", "2", "3", "4", "5"};
0197   std::vector<std::string> hitProperties = {"r", "x", "y", "z", "eta", "phi", "detId", "layer", "moduleType"};
0198 
0199   for (const auto& idx : hitIndices) {
0200     for (const auto& prop : hitProperties) {
0201       std::string branchName = "t5_t3_" + idx + "_" + prop;
0202       if (prop == "detId" || prop == "layer" || prop == "moduleType") {
0203         ana.tx->createBranch<std::vector<int>>(branchName);
0204       } else {
0205         ana.tx->createBranch<std::vector<float>>(branchName);
0206       }
0207     }
0208   }
0209 }
0210 
0211 //________________________________________________________________________________________________________________________________
0212 void createT3DNNBranches() {
0213   // Common branches for T3 properties based on TripletsSoA fields
0214   ana.tx->createBranch<std::vector<float>>("t3_betaIn");
0215   ana.tx->createBranch<std::vector<float>>("t3_centerX");
0216   ana.tx->createBranch<std::vector<float>>("t3_centerY");
0217   ana.tx->createBranch<std::vector<float>>("t3_radius");
0218   ana.tx->createBranch<std::vector<bool>>("t3_partOfPT5");
0219   ana.tx->createBranch<std::vector<bool>>("t3_partOfT5");
0220   ana.tx->createBranch<std::vector<bool>>("t3_partOfPT3");
0221   ana.tx->createBranch<std::vector<float>>("t3_pMatched");
0222   ana.tx->createBranch<std::vector<float>>("t3_sim_vxy");
0223   ana.tx->createBranch<std::vector<float>>("t3_sim_vz");
0224 
0225   // Hit-specific branches (T3 has 4 hits from two segments)
0226   std::vector<std::string> hitIndices = {"0", "1", "2", "3", "4", "5"};
0227   std::vector<std::string> hitProperties = {"r", "x", "y", "z", "eta", "phi", "detId", "layer", "moduleType"};
0228 
0229   for (const auto& idx : hitIndices) {
0230     for (const auto& prop : hitProperties) {
0231       std::string branchName = "t3_hit_" + idx + "_" + prop;
0232       if (prop == "detId" || prop == "layer" || prop == "moduleType") {
0233         ana.tx->createBranch<std::vector<int>>(branchName);
0234       } else {
0235         ana.tx->createBranch<std::vector<float>>(branchName);
0236       }
0237     }
0238   }
0239 
0240   // Additional metadata branches
0241   ana.tx->createBranch<std::vector<int>>("t3_layer_binary");
0242   ana.tx->createBranch<std::vector<std::vector<int>>>("t3_matched_simIdx");
0243 }
0244 
0245 //________________________________________________________________________________________________________________________________
0246 void createGnnNtupleBranches() {
0247   // Mini Doublets
0248   ana.tx->createBranch<std::vector<float>>("MD_pt");
0249   ana.tx->createBranch<std::vector<float>>("MD_eta");
0250   ana.tx->createBranch<std::vector<float>>("MD_phi");
0251   ana.tx->createBranch<std::vector<float>>("MD_dphichange");
0252   ana.tx->createBranch<std::vector<int>>("MD_isFake");
0253   ana.tx->createBranch<std::vector<int>>("MD_tpType");
0254   ana.tx->createBranch<std::vector<int>>("MD_detId");
0255   ana.tx->createBranch<std::vector<int>>("MD_layer");
0256   ana.tx->createBranch<std::vector<float>>("MD_0_r");
0257   ana.tx->createBranch<std::vector<float>>("MD_0_x");
0258   ana.tx->createBranch<std::vector<float>>("MD_0_y");
0259   ana.tx->createBranch<std::vector<float>>("MD_0_z");
0260   ana.tx->createBranch<std::vector<float>>("MD_1_r");
0261   ana.tx->createBranch<std::vector<float>>("MD_1_x");
0262   ana.tx->createBranch<std::vector<float>>("MD_1_y");
0263   ana.tx->createBranch<std::vector<float>>("MD_1_z");
0264 
0265   // Line Segments
0266   ana.tx->createBranch<std::vector<float>>("LS_pt");
0267   ana.tx->createBranch<std::vector<float>>("LS_eta");
0268   ana.tx->createBranch<std::vector<float>>("LS_phi");
0269   ana.tx->createBranch<std::vector<int>>("LS_isFake");
0270   ana.tx->createBranch<std::vector<int>>("LS_MD_idx0");
0271   ana.tx->createBranch<std::vector<int>>("LS_MD_idx1");
0272   ana.tx->createBranch<std::vector<float>>("LS_sim_pt");
0273   ana.tx->createBranch<std::vector<float>>("LS_sim_eta");
0274   ana.tx->createBranch<std::vector<float>>("LS_sim_phi");
0275   ana.tx->createBranch<std::vector<float>>("LS_sim_pca_dxy");
0276   ana.tx->createBranch<std::vector<float>>("LS_sim_pca_dz");
0277   ana.tx->createBranch<std::vector<int>>("LS_sim_q");
0278   ana.tx->createBranch<std::vector<int>>("LS_sim_pdgId");
0279   ana.tx->createBranch<std::vector<int>>("LS_sim_event");
0280   ana.tx->createBranch<std::vector<int>>("LS_sim_bx");
0281   ana.tx->createBranch<std::vector<float>>("LS_sim_vx");
0282   ana.tx->createBranch<std::vector<float>>("LS_sim_vy");
0283   ana.tx->createBranch<std::vector<float>>("LS_sim_vz");
0284   ana.tx->createBranch<std::vector<int>>("LS_isInTrueTC");
0285 
0286   // TC's LS
0287   ana.tx->createBranch<std::vector<std::vector<int>>>("tc_lsIdx");
0288 }
0289 
0290 //________________________________________________________________________________________________________________________________
0291 void setOutputBranches(LSTEvent* event) {
0292   // ============ Sim tracks =============
0293   int n_accepted_simtrk = 0;
0294   auto const& trk_sim_bunchCrossing = trk.getVI("sim_bunchCrossing");
0295   auto const& trk_sim_event = trk.getVI("sim_event");
0296   auto const& trk_sim_pt = trk.getVF("sim_pt");
0297   auto const& trk_sim_eta = trk.getVF("sim_eta");
0298   auto const& trk_sim_phi = trk.getVF("sim_phi");
0299   auto const& trk_sim_pca_dxy = trk.getVF("sim_pca_dxy");
0300   auto const& trk_sim_pca_dz = trk.getVF("sim_pca_dz");
0301   auto const& trk_sim_q = trk.getVI("sim_q");
0302   auto const& trk_sim_pdgId = trk.getVI("sim_pdgId");
0303   auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx");
0304   auto const& trk_simvtx_x = trk.getVF("simvtx_x");
0305   auto const& trk_simvtx_y = trk.getVF("simvtx_y");
0306   auto const& trk_simvtx_z = trk.getVF("simvtx_z");
0307   auto const& trk_ph2_x = trk.getVF("ph2_x");
0308   auto const& trk_ph2_y = trk.getVF("ph2_y");
0309   auto const& trk_ph2_z = trk.getVF("ph2_z");
0310   auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx");
0311   auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx");
0312   auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx");
0313   for (unsigned int isimtrk = 0; isimtrk < trk_sim_pt.size(); ++isimtrk) {
0314     // Skip out-of-time pileup
0315     if (trk_sim_bunchCrossing[isimtrk] != 0)
0316       continue;
0317 
0318     // Skip non-hard-scatter
0319     if (trk_sim_event[isimtrk] != 0)
0320       continue;
0321 
0322     ana.tx->pushbackToBranch<float>("sim_pt", trk_sim_pt[isimtrk]);
0323     ana.tx->pushbackToBranch<float>("sim_eta", trk_sim_eta[isimtrk]);
0324     ana.tx->pushbackToBranch<float>("sim_phi", trk_sim_phi[isimtrk]);
0325     ana.tx->pushbackToBranch<float>("sim_pca_dxy", trk_sim_pca_dxy[isimtrk]);
0326     ana.tx->pushbackToBranch<float>("sim_pca_dz", trk_sim_pca_dz[isimtrk]);
0327     ana.tx->pushbackToBranch<int>("sim_q", trk_sim_q[isimtrk]);
0328     ana.tx->pushbackToBranch<int>("sim_event", trk_sim_event[isimtrk]);
0329     ana.tx->pushbackToBranch<int>("sim_pdgId", trk_sim_pdgId[isimtrk]);
0330 
0331     // For vertex we need to look it up from simvtx info
0332     int vtxidx = trk_sim_parentVtxIdx[isimtrk];
0333     ana.tx->pushbackToBranch<float>("sim_vx", trk_simvtx_x[vtxidx]);
0334     ana.tx->pushbackToBranch<float>("sim_vy", trk_simvtx_y[vtxidx]);
0335     ana.tx->pushbackToBranch<float>("sim_vz", trk_simvtx_z[vtxidx]);
0336 
0337     // The trkNtupIdx is the idx in the trackingNtuple
0338     ana.tx->pushbackToBranch<float>("sim_trkNtupIdx", isimtrk);
0339 
0340     // Increase the counter for accepted simtrk
0341     n_accepted_simtrk++;
0342   }
0343 
0344   // Intermediate variables to keep track of matched track candidates for a given sim track
0345   std::vector<int> sim_TC_matched(n_accepted_simtrk);
0346   std::vector<int> sim_TC_matched_mask(n_accepted_simtrk);
0347   std::vector<int> sim_TC_matched_for_duplicate(trk_sim_pt.size());
0348 
0349   // Intermediate variables to keep track of matched sim tracks for a given track candidate
0350   std::vector<std::vector<int>> tc_matched_simIdx;
0351 
0352   // ============ Track candidates =============
0353   auto const& trackCandidates = event->getTrackCandidates();
0354   unsigned int nTrackCandidates = trackCandidates.nTrackCandidates();
0355   for (unsigned int idx = 0; idx < nTrackCandidates; idx++) {
0356     // Compute reco quantities of track candidate based on final object
0357     int type, isFake;
0358     float pt, eta, phi;
0359     std::vector<int> simidx;
0360     std::tie(type, pt, eta, phi, isFake, simidx) = parseTrackCandidate(
0361         event, idx, trk_ph2_x, trk_ph2_y, trk_ph2_z, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
0362     ana.tx->pushbackToBranch<float>("tc_pt", pt);
0363     ana.tx->pushbackToBranch<float>("tc_eta", eta);
0364     ana.tx->pushbackToBranch<float>("tc_phi", phi);
0365     ana.tx->pushbackToBranch<int>("tc_type", type);
0366     ana.tx->pushbackToBranch<int>("tc_isFake", isFake);
0367     tc_matched_simIdx.push_back(simidx);
0368 
0369     // Loop over matched sim idx and increase counter of TC_matched
0370     for (auto& idx : simidx) {
0371       // NOTE Important to note that the idx of the std::vector<> is same
0372       // as the tracking-ntuple's sim track idx ONLY because event==0 and bunchCrossing==0 condition is applied!!
0373       // Also do not try to access beyond the event and bunchCrossing
0374       if (idx < n_accepted_simtrk) {
0375         sim_TC_matched.at(idx) += 1;
0376         sim_TC_matched_mask.at(idx) |= (1 << type);
0377       }
0378       sim_TC_matched_for_duplicate.at(idx) += 1;
0379     }
0380   }
0381 
0382   // Using the intermedaite variables to compute whether a given track candidate is a duplicate
0383   std::vector<int> tc_isDuplicate(tc_matched_simIdx.size());
0384   // Loop over the track candidates
0385   for (unsigned int i = 0; i < tc_matched_simIdx.size(); ++i) {
0386     bool isDuplicate = false;
0387     // Loop over the sim idx matched to this track candidate
0388     for (unsigned int isim = 0; isim < tc_matched_simIdx[i].size(); ++isim) {
0389       // Using the sim_TC_matched to see whether this track candidate is matched to a sim track that is matched to more than one
0390       int simidx = tc_matched_simIdx[i][isim];
0391       if (sim_TC_matched_for_duplicate[simidx] > 1) {
0392         isDuplicate = true;
0393       }
0394     }
0395     tc_isDuplicate[i] = isDuplicate;
0396   }
0397 
0398   // Now set the last remaining branches
0399   ana.tx->setBranch<std::vector<int>>("sim_TC_matched", sim_TC_matched);
0400   ana.tx->setBranch<std::vector<int>>("sim_TC_matched_mask", sim_TC_matched_mask);
0401   ana.tx->setBranch<std::vector<std::vector<int>>>("tc_matched_simIdx", tc_matched_simIdx);
0402   ana.tx->setBranch<std::vector<int>>("tc_isDuplicate", tc_isDuplicate);
0403 }
0404 
0405 //________________________________________________________________________________________________________________________________
0406 void setOptionalOutputBranches(LSTEvent* event) {
0407 #ifdef CUT_VALUE_DEBUG
0408 
0409   setPixelQuintupletOutputBranches(event);
0410   setQuintupletOutputBranches(event);
0411   setPixelTripletOutputBranches(event);
0412   setOccupancyBranches(event);
0413   setT3DNNBranches(event);
0414   setT5DNNBranches(event);
0415   setpT3DNNBranches(event);
0416   setpLSOutputBranches(event);
0417 
0418 #endif
0419 }
0420 
0421 //________________________________________________________________________________________________________________________________
0422 void setOccupancyBranches(LSTEvent* event) {
0423   auto modules = event->getModules<ModulesSoA>();
0424   auto miniDoublets = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
0425   auto segments = event->getSegments<SegmentsOccupancySoA>();
0426   auto triplets = event->getTriplets<TripletsOccupancySoA>();
0427   auto quintuplets = event->getQuintuplets<QuintupletsOccupancySoA>();
0428   auto pixelQuintuplets = event->getPixelQuintuplets();
0429   auto pixelTriplets = event->getPixelTriplets();
0430   auto trackCandidates = event->getTrackCandidates();
0431 
0432   std::vector<int> moduleLayer;
0433   std::vector<int> moduleSubdet;
0434   std::vector<int> moduleRing;
0435   std::vector<int> moduleRod;
0436   std::vector<int> moduleModule;
0437   std::vector<float> moduleEta;
0438   std::vector<float> moduleR;
0439   std::vector<bool> moduleIsTilted;
0440   std::vector<int> trackCandidateOccupancy;
0441   std::vector<int> tripletOccupancy;
0442   std::vector<int> segmentOccupancy;
0443   std::vector<int> mdOccupancy;
0444   std::vector<int> quintupletOccupancy;
0445 
0446   for (unsigned int lowerIdx = 0; lowerIdx <= modules.nLowerModules(); lowerIdx++) {
0447     //layer = 0, subdet = 0 => pixel module
0448     moduleLayer.push_back(modules.layers()[lowerIdx]);
0449     moduleSubdet.push_back(modules.subdets()[lowerIdx]);
0450     moduleRing.push_back(modules.rings()[lowerIdx]);
0451     moduleRod.push_back(modules.rods()[lowerIdx]);
0452     moduleEta.push_back(modules.eta()[lowerIdx]);
0453     moduleR.push_back(modules.r()[lowerIdx]);
0454     bool isTilted = (modules.subdets()[lowerIdx] == 5 and modules.sides()[lowerIdx] != 3);
0455     moduleIsTilted.push_back(isTilted);
0456     moduleModule.push_back(modules.modules()[lowerIdx]);
0457     segmentOccupancy.push_back(segments.totOccupancySegments()[lowerIdx]);
0458     mdOccupancy.push_back(miniDoublets.totOccupancyMDs()[lowerIdx]);
0459 
0460     if (lowerIdx < modules.nLowerModules()) {
0461       quintupletOccupancy.push_back(quintuplets.totOccupancyQuintuplets()[lowerIdx]);
0462       tripletOccupancy.push_back(triplets.totOccupancyTriplets()[lowerIdx]);
0463     }
0464   }
0465 
0466   ana.tx->setBranch<std::vector<int>>("module_layers", moduleLayer);
0467   ana.tx->setBranch<std::vector<int>>("module_subdets", moduleSubdet);
0468   ana.tx->setBranch<std::vector<int>>("module_rings", moduleRing);
0469   ana.tx->setBranch<std::vector<int>>("module_rods", moduleRod);
0470   ana.tx->setBranch<std::vector<int>>("module_modules", moduleModule);
0471   ana.tx->setBranch<std::vector<bool>>("module_isTilted", moduleIsTilted);
0472   ana.tx->setBranch<std::vector<float>>("module_eta", moduleEta);
0473   ana.tx->setBranch<std::vector<float>>("module_r", moduleR);
0474   ana.tx->setBranch<std::vector<int>>("md_occupancies", mdOccupancy);
0475   ana.tx->setBranch<std::vector<int>>("sg_occupancies", segmentOccupancy);
0476   ana.tx->setBranch<std::vector<int>>("t3_occupancies", tripletOccupancy);
0477   ana.tx->setBranch<int>("tc_occupancies", trackCandidates.nTrackCandidates());
0478   ana.tx->setBranch<int>("pT3_occupancies", pixelTriplets.totOccupancyPixelTriplets());
0479   ana.tx->setBranch<std::vector<int>>("t5_occupancies", quintupletOccupancy);
0480   ana.tx->setBranch<int>("pT5_occupancies", pixelQuintuplets.totOccupancyPixelQuintuplets());
0481 }
0482 
0483 //________________________________________________________________________________________________________________________________
0484 void setPixelQuintupletOutputBranches(LSTEvent* event) {
0485   // ============ pT5 =============
0486   auto const pixelQuintuplets = event->getPixelQuintuplets();
0487   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
0488   auto const pixelSeeds = event->getInput<PixelSeedsSoA>();
0489   auto modules = event->getModules<ModulesSoA>();
0490   int n_accepted_simtrk = ana.tx->getBranch<std::vector<int>>("sim_TC_matched").size();
0491 
0492   unsigned int nPixelQuintuplets = pixelQuintuplets.nPixelQuintuplets();
0493   std::vector<int> sim_pT5_matched(n_accepted_simtrk);
0494   std::vector<std::vector<int>> pT5_matched_simIdx;
0495 
0496   auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx");
0497   auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx");
0498   auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx");
0499 
0500   for (unsigned int pT5 = 0; pT5 < nPixelQuintuplets; pT5++) {
0501     unsigned int T5Index = getT5FrompT5(event, pT5);
0502     unsigned int pLSIndex = getPixelLSFrompT5(event, pT5);
0503     float pt = (__H2F(quintuplets.innerRadius()[T5Index]) * k2Rinv1GeVf * 2 + pixelSeeds.ptIn()[pLSIndex]) / 2;
0504     float eta = pixelSeeds.eta()[pLSIndex];
0505     float phi = pixelSeeds.phi()[pLSIndex];
0506 
0507     std::vector<unsigned int> hit_idx = getHitIdxsFrompT5(event, pT5);
0508     std::vector<unsigned int> module_idx = getModuleIdxsFrompT5(event, pT5);
0509     std::vector<unsigned int> hit_type = getHitTypesFrompT5(event, pT5);
0510 
0511     int layer_binary = 1;
0512     int moduleType_binary = 0;
0513     for (size_t i = 0; i < module_idx.size(); i += 2) {
0514       layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4)));
0515       moduleType_binary |= (modules.moduleType()[module_idx[i]] << i);
0516     }
0517     std::vector<int> simidx =
0518         matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
0519     ana.tx->pushbackToBranch<int>("pT5_isFake", static_cast<int>(simidx.size() == 0));
0520     ana.tx->pushbackToBranch<float>("pT5_pt", pt);
0521     ana.tx->pushbackToBranch<float>("pT5_eta", eta);
0522     ana.tx->pushbackToBranch<float>("pT5_phi", phi);
0523     ana.tx->pushbackToBranch<int>("pT5_layer_binary", layer_binary);
0524     ana.tx->pushbackToBranch<int>("pT5_moduleType_binary", moduleType_binary);
0525     ana.tx->pushbackToBranch<float>("pT5_rzChiSquared", pixelQuintuplets.rzChiSquared()[pT5]);
0526 
0527     pT5_matched_simIdx.push_back(simidx);
0528 
0529     // Loop over matched sim idx and increase counter of pT5_matched
0530     for (auto& idx : simidx) {
0531       // NOTE Important to note that the idx of the std::vector<> is same
0532       // as the tracking-ntuple's sim track idx ONLY because event==0 and bunchCrossing==0 condition is applied!!
0533       // Also do not try to access beyond the event and bunchCrossing
0534       if (idx < n_accepted_simtrk) {
0535         sim_pT5_matched.at(idx) += 1;
0536       }
0537     }
0538   }
0539 
0540   // Using the intermedaite variables to compute whether a given track candidate is a duplicate
0541   std::vector<int> pT5_isDuplicate(pT5_matched_simIdx.size());
0542   // Loop over the track candidates
0543   for (unsigned int i = 0; i < pT5_matched_simIdx.size(); ++i) {
0544     bool isDuplicate = false;
0545     // Loop over the sim idx matched to this track candidate
0546     for (unsigned int isim = 0; isim < pT5_matched_simIdx[i].size(); ++isim) {
0547       // Using the sim_pT5_matched to see whether this track candidate is matched to a sim track that is matched to more than one
0548       int simidx = pT5_matched_simIdx[i][isim];
0549       if (simidx < n_accepted_simtrk) {
0550         if (sim_pT5_matched[simidx] > 1) {
0551           isDuplicate = true;
0552         }
0553       }
0554     }
0555     pT5_isDuplicate[i] = isDuplicate;
0556   }
0557 
0558   // Now set the last remaining branches
0559   ana.tx->setBranch<std::vector<int>>("sim_pT5_matched", sim_pT5_matched);
0560   ana.tx->setBranch<std::vector<std::vector<int>>>("pT5_matched_simIdx", pT5_matched_simIdx);
0561   ana.tx->setBranch<std::vector<int>>("pT5_isDuplicate", pT5_isDuplicate);
0562 }
0563 
0564 //________________________________________________________________________________________________________________________________
0565 void setQuintupletOutputBranches(LSTEvent* event) {
0566   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
0567   auto const quintupletsOccupancy = event->getQuintuplets<QuintupletsOccupancySoA>();
0568   auto ranges = event->getRanges();
0569   auto modules = event->getModules<ModulesSoA>();
0570   int n_accepted_simtrk = ana.tx->getBranch<std::vector<int>>("sim_TC_matched").size();
0571 
0572   std::vector<int> sim_t5_matched(n_accepted_simtrk);
0573   std::vector<std::vector<int>> t5_matched_simIdx;
0574 
0575   auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx");
0576   auto const& trk_simvtx_x = trk.getVF("simvtx_x");
0577   auto const& trk_simvtx_y = trk.getVF("simvtx_y");
0578   auto const& trk_simvtx_z = trk.getVF("simvtx_z");
0579   auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx");
0580   auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx");
0581   auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx");
0582 
0583   for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < modules.nLowerModules(); ++lowerModuleIdx) {
0584     int nQuintuplets = quintupletsOccupancy.nQuintuplets()[lowerModuleIdx];
0585     for (unsigned int idx = 0; idx < nQuintuplets; idx++) {
0586       unsigned int quintupletIndex = ranges.quintupletModuleIndices()[lowerModuleIdx] + idx;
0587       float pt = __H2F(quintuplets.innerRadius()[quintupletIndex]) * k2Rinv1GeVf * 2;
0588       float eta = __H2F(quintuplets.eta()[quintupletIndex]);
0589       float phi = __H2F(quintuplets.phi()[quintupletIndex]);
0590 
0591       std::vector<unsigned int> hit_idx = getHitIdxsFromT5(event, quintupletIndex);
0592       std::vector<unsigned int> hit_type = getHitTypesFromT5(event, quintupletIndex);
0593       std::vector<unsigned int> module_idx = getModuleIdxsFromT5(event, quintupletIndex);
0594 
0595       int layer_binary = 0;
0596       int moduleType_binary = 0;
0597       for (size_t i = 0; i < module_idx.size(); i += 2) {
0598         layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4)));
0599         moduleType_binary |= (modules.moduleType()[module_idx[i]] << i);
0600       }
0601 
0602       float percent_matched;
0603       std::vector<int> simidx = matchedSimTrkIdxs(
0604           hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, &percent_matched);
0605 
0606       ana.tx->pushbackToBranch<int>("t5_isFake", static_cast<int>(simidx.size() == 0));
0607       ana.tx->pushbackToBranch<float>("t5_pt", pt);
0608       ana.tx->pushbackToBranch<float>("t5_pMatched", percent_matched);
0609       ana.tx->pushbackToBranch<float>("t5_eta", eta);
0610       ana.tx->pushbackToBranch<float>("t5_phi", phi);
0611       ana.tx->pushbackToBranch<float>("t5_innerRadius", __H2F(quintuplets.innerRadius()[quintupletIndex]));
0612       ana.tx->pushbackToBranch<float>("t5_bridgeRadius", __H2F(quintuplets.bridgeRadius()[quintupletIndex]));
0613       ana.tx->pushbackToBranch<float>("t5_outerRadius", __H2F(quintuplets.outerRadius()[quintupletIndex]));
0614       ana.tx->pushbackToBranch<int>("t5_isDupAlgoFlag", quintuplets.isDup()[quintupletIndex]);
0615       ana.tx->pushbackToBranch<float>("t5_chiSquared", quintuplets.chiSquared()[quintupletIndex]);
0616       ana.tx->pushbackToBranch<float>("t5_rzChiSquared", quintuplets.rzChiSquared()[quintupletIndex]);
0617       ana.tx->pushbackToBranch<float>("t5_nonAnchorChiSquared", quintuplets.nonAnchorChiSquared()[quintupletIndex]);
0618       ana.tx->pushbackToBranch<float>("t5_dBeta1", quintuplets.dBeta1()[quintupletIndex]);
0619       ana.tx->pushbackToBranch<float>("t5_dBeta2", quintuplets.dBeta2()[quintupletIndex]);
0620       ana.tx->pushbackToBranch<int>("t5_layer_binary", layer_binary);
0621       ana.tx->pushbackToBranch<int>("t5_moduleType_binary", moduleType_binary);
0622 
0623       t5_matched_simIdx.push_back(simidx);
0624 
0625       for (auto& simtrk : simidx) {
0626         if (simtrk < n_accepted_simtrk) {
0627           sim_t5_matched.at(simtrk) += 1;
0628         }
0629       }
0630 
0631       // Avoid fakes when calculating the vertex distance, set default to 0.0.
0632       if (simidx.size() == 0) {
0633         ana.tx->pushbackToBranch<float>("t5_sim_vxy", 0.0);
0634         ana.tx->pushbackToBranch<float>("t5_sim_vz", 0.0);
0635         continue;
0636       }
0637 
0638       int vtxidx = trk_sim_parentVtxIdx[simidx[0]];
0639       float vtx_x = trk_simvtx_x[vtxidx];
0640       float vtx_y = trk_simvtx_y[vtxidx];
0641       float vtx_z = trk_simvtx_z[vtxidx];
0642 
0643       ana.tx->pushbackToBranch<float>("t5_sim_vxy", sqrt(vtx_x * vtx_x + vtx_y * vtx_y));
0644       ana.tx->pushbackToBranch<float>("t5_sim_vz", vtx_z);
0645     }
0646   }
0647 
0648   std::vector<int> t5_isDuplicate(t5_matched_simIdx.size());
0649   for (unsigned int i = 0; i < t5_matched_simIdx.size(); i++) {
0650     bool isDuplicate = false;
0651     for (unsigned int isim = 0; isim < t5_matched_simIdx[i].size(); isim++) {
0652       int simidx = t5_matched_simIdx[i][isim];
0653       if (simidx < n_accepted_simtrk) {
0654         if (sim_t5_matched[simidx] > 1) {
0655           isDuplicate = true;
0656         }
0657       }
0658     }
0659     t5_isDuplicate[i] = isDuplicate;
0660   }
0661   ana.tx->setBranch<std::vector<int>>("sim_T5_matched", sim_t5_matched);
0662   ana.tx->setBranch<std::vector<std::vector<int>>>("t5_matched_simIdx", t5_matched_simIdx);
0663   ana.tx->setBranch<std::vector<int>>("t5_isDuplicate", t5_isDuplicate);
0664 }
0665 
0666 //________________________________________________________________________________________________________________________________
0667 void setPixelTripletOutputBranches(LSTEvent* event) {
0668   auto const pixelTriplets = event->getPixelTriplets();
0669   auto modules = event->getModules<ModulesSoA>();
0670   PixelSeedsConst pixelSeeds = event->getInput<PixelSeedsSoA>();
0671   int n_accepted_simtrk = ana.tx->getBranch<std::vector<int>>("sim_TC_matched").size();
0672 
0673   unsigned int nPixelTriplets = pixelTriplets.nPixelTriplets();
0674   std::vector<int> sim_pT3_matched(n_accepted_simtrk);
0675   std::vector<std::vector<int>> pT3_matched_simIdx;
0676 
0677   auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx");
0678   auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx");
0679   auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx");
0680 
0681   for (unsigned int pT3 = 0; pT3 < nPixelTriplets; pT3++) {
0682     unsigned int T3Index = getT3FrompT3(event, pT3);
0683     unsigned int pLSIndex = getPixelLSFrompT3(event, pT3);
0684     const float pt = pixelSeeds.ptIn()[pLSIndex];
0685 
0686     float eta = pixelSeeds.eta()[pLSIndex];
0687     float phi = pixelSeeds.phi()[pLSIndex];
0688     std::vector<unsigned int> hit_idx = getHitIdxsFrompT3(event, pT3);
0689     std::vector<unsigned int> hit_type = getHitTypesFrompT3(event, pT3);
0690 
0691     std::vector<int> simidx =
0692         matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
0693     std::vector<unsigned int> module_idx = getModuleIdxsFrompT3(event, pT3);
0694     int layer_binary = 1;
0695     int moduleType_binary = 0;
0696     for (size_t i = 0; i < module_idx.size(); i += 2) {
0697       layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4)));
0698       moduleType_binary |= (modules.moduleType()[module_idx[i]] << i);
0699     }
0700     ana.tx->pushbackToBranch<int>("pT3_isFake", static_cast<int>(simidx.size() == 0));
0701     ana.tx->pushbackToBranch<float>("pT3_pt", pt);
0702     ana.tx->pushbackToBranch<float>("pT3_eta", eta);
0703     ana.tx->pushbackToBranch<float>("pT3_phi", phi);
0704     ana.tx->pushbackToBranch<int>("pT3_layer_binary", layer_binary);
0705     ana.tx->pushbackToBranch<int>("pT3_moduleType_binary", moduleType_binary);
0706 
0707     pT3_matched_simIdx.push_back(simidx);
0708 
0709     for (auto& idx : simidx) {
0710       if (idx < n_accepted_simtrk) {
0711         sim_pT3_matched.at(idx) += 1;
0712       }
0713     }
0714   }
0715 
0716   std::vector<int> pT3_isDuplicate(pT3_matched_simIdx.size());
0717   for (unsigned int i = 0; i < pT3_matched_simIdx.size(); i++) {
0718     bool isDuplicate = true;
0719     for (unsigned int isim = 0; isim < pT3_matched_simIdx[i].size(); isim++) {
0720       int simidx = pT3_matched_simIdx[i][isim];
0721       if (simidx < n_accepted_simtrk) {
0722         if (sim_pT3_matched[simidx] > 1) {
0723           isDuplicate = true;
0724         }
0725       }
0726     }
0727     pT3_isDuplicate[i] = isDuplicate;
0728   }
0729   ana.tx->setBranch<std::vector<int>>("sim_pT3_matched", sim_pT3_matched);
0730   ana.tx->setBranch<std::vector<std::vector<int>>>("pT3_matched_simIdx", pT3_matched_simIdx);
0731   ana.tx->setBranch<std::vector<int>>("pT3_isDuplicate", pT3_isDuplicate);
0732 }
0733 
0734 //________________________________________________________________________________________________________________________________
0735 void fillpT3DNNBranches(LSTEvent* event, unsigned int iPT3) {
0736   // Retrieve the pT3 object from the PixelTriplets SoA.
0737   auto pixelTriplets = event->getPixelTriplets();
0738 
0739   float pixelRadius = pixelTriplets.pixelRadius()[iPT3];
0740   float pixelRadiusError = pixelTriplets.pixelRadiusError()[iPT3];
0741   float tripletRadius = pixelTriplets.tripletRadius()[iPT3];
0742   float phi = pixelTriplets.phi()[iPT3];          // from the T3
0743   float phi_pix = pixelTriplets.phi_pix()[iPT3];  // from the pLS
0744   float rPhiChiSquared = pixelTriplets.rPhiChiSquared()[iPT3];
0745   float rPhiChiSquaredInwards = pixelTriplets.rPhiChiSquaredInwards()[iPT3];
0746   float rzChiSquared = pixelTriplets.rzChiSquared()[iPT3];
0747   float pt = pixelTriplets.pt()[iPT3];
0748   float eta = pixelTriplets.eta()[iPT3];
0749   float eta_pix = pixelTriplets.eta_pix()[iPT3];  // eta from pLS
0750   float centerX = pixelTriplets.centerX()[iPT3];  // T3-based circle center x
0751   float centerY = pixelTriplets.centerY()[iPT3];  // T3-based circle center y
0752 
0753   ana.tx->pushbackToBranch<float>("pT3_rPhiChiSquared", rPhiChiSquared);
0754   ana.tx->pushbackToBranch<float>("pT3_rPhiChiSquaredInwards", rPhiChiSquaredInwards);
0755   ana.tx->pushbackToBranch<float>("pT3_rzChiSquared", rzChiSquared);
0756   ana.tx->pushbackToBranch<float>("pT3_pixelRadius", pixelRadius);
0757   ana.tx->pushbackToBranch<float>("pT3_pixelRadiusError", pixelRadiusError);
0758   ana.tx->pushbackToBranch<float>("pT3_tripletRadius", tripletRadius);
0759 }
0760 
0761 //________________________________________________________________________________________________________________________________
0762 void fillT3DNNBranches(LSTEvent* event, unsigned int iT3) {
0763   auto hitsBase = event->getInput<HitsBaseSoA>();
0764   auto hitsExtended = event->getHits<HitsExtendedSoA>();
0765   auto modules = event->getModules<ModulesSoA>();
0766 
0767   std::vector<unsigned int> hitIdx = getHitsFromT3(event, iT3);
0768   std::vector<lst_math::Hit> hitObjects;
0769 
0770   auto const& trk_ph2_subdet = trk.getVUS("ph2_subdet");
0771   auto const& trk_ph2_layer = trk.getVUS("ph2_layer");
0772   auto const& trk_ph2_detId = trk.getVU("ph2_detId");
0773 
0774   for (int i = 0; i < hitIdx.size(); ++i) {
0775     unsigned int hit = hitIdx[i];
0776     float x = hitsBase.xs()[hit];
0777     float y = hitsBase.ys()[hit];
0778     float z = hitsBase.zs()[hit];
0779     lst_math::Hit hitObj(x, y, z);
0780     hitObjects.push_back(hitObj);
0781 
0782     std::string idx = std::to_string(i);
0783     ana.tx->pushbackToBranch<float>("t3_hit_" + idx + "_r", sqrt(x * x + y * y));
0784     ana.tx->pushbackToBranch<float>("t3_hit_" + idx + "_x", x);
0785     ana.tx->pushbackToBranch<float>("t3_hit_" + idx + "_y", y);
0786     ana.tx->pushbackToBranch<float>("t3_hit_" + idx + "_z", z);
0787     ana.tx->pushbackToBranch<float>("t3_hit_" + idx + "_eta", hitObj.eta());
0788     ana.tx->pushbackToBranch<float>("t3_hit_" + idx + "_phi", hitObj.phi());
0789 
0790     int subdet = trk_ph2_subdet[hitsBase.idxs()[hit]];
0791     int is_endcap = subdet == 4;
0792     int layer = trk_ph2_layer[hitsBase.idxs()[hit]] + 6 * is_endcap;
0793     int detId = trk_ph2_detId[hitsBase.idxs()[hit]];
0794     unsigned int module = hitsExtended.moduleIndices()[hit];
0795 
0796     ana.tx->pushbackToBranch<int>("t3_hit_" + idx + "_detId", detId);
0797     ana.tx->pushbackToBranch<int>("t3_hit_" + idx + "_layer", layer);
0798     ana.tx->pushbackToBranch<int>("t3_hit_" + idx + "_moduleType", modules.moduleType()[module]);
0799   }
0800 }
0801 
0802 //________________________________________________________________________________________________________________________________
0803 void fillT5DNNBranches(LSTEvent* event, unsigned int iT3) {
0804   auto hitsBase = event->getInput<HitsBaseSoA>();
0805   auto hitsExtended = event->getHits<HitsExtendedSoA>();
0806   auto modules = event->getModules<ModulesSoA>();
0807 
0808   std::vector<unsigned int> hitIdx = getHitsFromT3(event, iT3);
0809   std::vector<lst_math::Hit> hitObjects(hitIdx.size());
0810 
0811   auto const& trk_ph2_subdet = trk.getVUS("ph2_subdet");
0812   auto const& trk_ph2_layer = trk.getVUS("ph2_layer");
0813   auto const& trk_ph2_detId = trk.getVU("ph2_detId");
0814 
0815   for (int i = 0; i < hitIdx.size(); ++i) {
0816     unsigned int hit = hitIdx[i];
0817     float x = hitsBase.xs()[hit];
0818     float y = hitsBase.ys()[hit];
0819     float z = hitsBase.zs()[hit];
0820     hitObjects[i] = lst_math::Hit(x, y, z);
0821 
0822     std::string idx = std::to_string(i);
0823     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_r", sqrt(x * x + y * y));
0824     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_x", x);
0825     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_y", y);
0826     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_z", z);
0827     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_eta", hitObjects[i].eta());
0828     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_phi", hitObjects[i].phi());
0829 
0830     int subdet = trk_ph2_subdet[hitsBase.idxs()[hit]];
0831     int is_endcap = subdet == 4;
0832     int layer = trk_ph2_layer[hitsBase.idxs()[hit]] + 6 * is_endcap;
0833     int detId = trk_ph2_detId[hitsBase.idxs()[hit]];
0834     unsigned int module = hitsExtended.moduleIndices()[hit];
0835 
0836     ana.tx->pushbackToBranch<int>("t5_t3_" + idx + "_detId", detId);
0837     ana.tx->pushbackToBranch<int>("t5_t3_" + idx + "_layer", layer);
0838     ana.tx->pushbackToBranch<int>("t5_t3_" + idx + "_moduleType", modules.moduleType()[module]);
0839   }
0840 
0841   float radius;
0842   auto const& devHost = cms::alpakatools::host();
0843   std::tie(radius, std::ignore, std::ignore) = computeRadiusFromThreeAnchorHits(devHost,
0844                                                                                 hitObjects[0].x(),
0845                                                                                 hitObjects[0].y(),
0846                                                                                 hitObjects[1].x(),
0847                                                                                 hitObjects[1].y(),
0848                                                                                 hitObjects[2].x(),
0849                                                                                 hitObjects[2].y());
0850   ana.tx->pushbackToBranch<float>("t5_t3_pt", k2Rinv1GeVf * 2 * radius);
0851 
0852   // Angles
0853   ana.tx->pushbackToBranch<float>("t5_t3_eta", hitObjects[2].eta());
0854   ana.tx->pushbackToBranch<float>("t5_t3_phi", hitObjects[0].phi());
0855 }
0856 
0857 //________________________________________________________________________________________________________________________________
0858 void setpT3DNNBranches(LSTEvent* event) {
0859   auto pixelTriplets = event->getPixelTriplets();
0860   unsigned int nPT3 = pixelTriplets.nPixelTriplets();
0861   for (unsigned int iPT3 = 0; iPT3 < nPT3; ++iPT3) {
0862     fillpT3DNNBranches(event, iPT3);
0863   }
0864 }
0865 
0866 //________________________________________________________________________________________________________________________________
0867 void setT3DNNBranches(LSTEvent* event) {
0868   auto const triplets = event->getTriplets<TripletsSoA>();
0869   auto const tripletsOccupancy = event->getTriplets<TripletsOccupancySoA>();
0870   auto modules = event->getModules<ModulesSoA>();
0871   auto ranges = event->getRanges();
0872 
0873   auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx");
0874   auto const& trk_simvtx_x = trk.getVF("simvtx_x");
0875   auto const& trk_simvtx_y = trk.getVF("simvtx_y");
0876   auto const& trk_simvtx_z = trk.getVF("simvtx_z");
0877   auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx");
0878   auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx");
0879   auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx");
0880 
0881   for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < modules.nLowerModules(); ++lowerModuleIdx) {
0882     int nTriplets = tripletsOccupancy.nTriplets()[lowerModuleIdx];
0883     for (unsigned int idx = 0; idx < nTriplets; idx++) {
0884       unsigned int tripletIndex = ranges.tripletModuleIndices()[lowerModuleIdx] + idx;
0885 
0886       // Get hit indices and types
0887       std::vector<unsigned int> hit_idx = getHitsFromT3(event, tripletIndex);
0888       std::vector<unsigned int> hit_type = getHitTypesFromT3(event, tripletIndex);
0889       std::vector<unsigned int> module_idx = getModuleIdxsFromT3(event, tripletIndex);
0890 
0891       // Calculate layer binary representation
0892       int layer_binary = 0;
0893       for (size_t i = 0; i < module_idx.size(); i += 2) {
0894         layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4)));
0895       }
0896 
0897       // Get matching information with percent matched
0898       float percent_matched;
0899       std::vector<int> simidx = matchedSimTrkIdxs(
0900           hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, &percent_matched);
0901 
0902       // Fill the branches with T3-specific data
0903       ana.tx->pushbackToBranch<float>("t3_betaIn", triplets.betaIn()[tripletIndex]);
0904       ana.tx->pushbackToBranch<float>("t3_centerX", triplets.centerX()[tripletIndex]);
0905       ana.tx->pushbackToBranch<float>("t3_centerY", triplets.centerY()[tripletIndex]);
0906       ana.tx->pushbackToBranch<float>("t3_radius", triplets.radius()[tripletIndex]);
0907       ana.tx->pushbackToBranch<bool>("t3_partOfPT5", triplets.partOfPT5()[tripletIndex]);
0908       ana.tx->pushbackToBranch<bool>("t3_partOfT5", triplets.partOfT5()[tripletIndex]);
0909       ana.tx->pushbackToBranch<bool>("t3_partOfPT3", triplets.partOfPT3()[tripletIndex]);
0910       ana.tx->pushbackToBranch<int>("t3_layer_binary", layer_binary);
0911       ana.tx->pushbackToBranch<std::vector<int>>("t3_matched_simIdx", simidx);
0912       ana.tx->pushbackToBranch<float>("t3_pMatched", percent_matched);
0913 
0914       // Add vertex information for matched sim tracks
0915       if (simidx.size() == 0) {
0916         // No matched sim track - set default values
0917         ana.tx->pushbackToBranch<float>("t3_sim_vxy", 0.0);
0918         ana.tx->pushbackToBranch<float>("t3_sim_vz", 0.0);
0919       } else {
0920         // Get vertex information from the first matched sim track
0921         int vtxidx = trk_sim_parentVtxIdx[simidx[0]];
0922         float vtx_x = trk_simvtx_x[vtxidx];
0923         float vtx_y = trk_simvtx_y[vtxidx];
0924         float vtx_z = trk_simvtx_z[vtxidx];
0925 
0926         // Calculate transverse distance from origin
0927         float vxy = sqrt(vtx_x * vtx_x + vtx_y * vtx_y);
0928 
0929         ana.tx->pushbackToBranch<float>("t3_sim_vxy", vxy);
0930         ana.tx->pushbackToBranch<float>("t3_sim_vz", vtx_z);
0931       }
0932 
0933       // Fill hit-specific information
0934       fillT3DNNBranches(event, tripletIndex);
0935     }
0936   }
0937 }
0938 
0939 //________________________________________________________________________________________________________________________________
0940 void setT5DNNBranches(LSTEvent* event) {
0941   auto tripletsOcc = event->getTriplets<TripletsOccupancySoA>();
0942   auto tripletsSoA = event->getTriplets<TripletsSoA>();
0943   auto modules = event->getModules<ModulesSoA>();
0944   auto ranges = event->getRanges();
0945   auto const quintuplets = event->getQuintuplets<QuintupletsOccupancySoA>();
0946   auto trackCandidates = event->getTrackCandidates();
0947 
0948   std::unordered_set<unsigned int> allT3s;
0949   std::unordered_map<unsigned int, unsigned int> t3_index_map;
0950 
0951   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
0952     for (unsigned int jdx = 0; jdx < tripletsOcc.nTriplets()[idx]; ++jdx) {
0953       unsigned int t3Idx = ranges.tripletModuleIndices()[idx] + jdx;
0954       if (allT3s.insert(t3Idx).second) {
0955         t3_index_map[t3Idx] = allT3s.size() - 1;
0956         fillT5DNNBranches(event, t3Idx);
0957       }
0958     }
0959   }
0960 
0961   std::unordered_map<unsigned int, unsigned int> t5_tc_index_map;
0962   std::unordered_set<unsigned int> t5s_used_in_tc;
0963 
0964   for (unsigned int idx = 0; idx < trackCandidates.nTrackCandidates(); idx++) {
0965     if (trackCandidates.trackCandidateType()[idx] == LSTObjType::T5) {
0966       unsigned int objIdx = trackCandidates.directObjectIndices()[idx];
0967       t5s_used_in_tc.insert(objIdx);
0968       t5_tc_index_map[objIdx] = idx;
0969     }
0970   }
0971 
0972   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
0973     for (unsigned int jdx = 0; jdx < quintuplets.nQuintuplets()[idx]; ++jdx) {
0974       unsigned int t5Idx = ranges.quintupletModuleIndices()[idx] + jdx;
0975       std::vector<unsigned int> t3sIdx = getT3sFromT5(event, t5Idx);
0976 
0977       ana.tx->pushbackToBranch<int>("t5_t3_idx0", t3_index_map[t3sIdx[0]]);
0978       ana.tx->pushbackToBranch<int>("t5_t3_idx1", t3_index_map[t3sIdx[1]]);
0979 
0980       ana.tx->pushbackToBranch<float>("t5_t3_fakeScore1", tripletsSoA.fakeScore()[t3sIdx[0]]);
0981       ana.tx->pushbackToBranch<float>("t5_t3_promptScore1", tripletsSoA.promptScore()[t3sIdx[0]]);
0982       ana.tx->pushbackToBranch<float>("t5_t3_displacedScore1", tripletsSoA.displacedScore()[t3sIdx[0]]);
0983       ana.tx->pushbackToBranch<float>("t5_t3_fakeScore2", tripletsSoA.fakeScore()[t3sIdx[1]]);
0984       ana.tx->pushbackToBranch<float>("t5_t3_promptScore2", tripletsSoA.promptScore()[t3sIdx[1]]);
0985       ana.tx->pushbackToBranch<float>("t5_t3_displacedScore2", tripletsSoA.displacedScore()[t3sIdx[1]]);
0986 
0987       if (t5s_used_in_tc.find(t5Idx) != t5s_used_in_tc.end()) {
0988         ana.tx->pushbackToBranch<int>("t5_partOfTC", 1);
0989         ana.tx->pushbackToBranch<int>("t5_tc_idx", t5_tc_index_map[t5Idx]);
0990       } else {
0991         ana.tx->pushbackToBranch<int>("t5_partOfTC", 0);
0992         ana.tx->pushbackToBranch<int>("t5_tc_idx", -999);
0993       }
0994     }
0995   }
0996 }
0997 
0998 //________________________________________________________________________________________________________________________________
0999 void setGnnNtupleBranches(LSTEvent* event) {
1000   // Get relevant information
1001   SegmentsOccupancyConst segmentsOccupancy = event->getSegments<SegmentsOccupancySoA>();
1002   MiniDoubletsOccupancyConst miniDoublets = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
1003   auto hitsBase = event->getInput<HitsBaseSoA>();
1004   auto modules = event->getModules<ModulesSoA>();
1005   auto ranges = event->getRanges();
1006   auto const& trackCandidates = event->getTrackCandidates();
1007 
1008   std::set<unsigned int> mds_used_in_sg;
1009   std::map<unsigned int, unsigned int> md_index_map;
1010   std::map<unsigned int, unsigned int> sg_index_map;
1011 
1012   // Loop over modules (lower ones where the MDs are saved)
1013   unsigned int nTotalMD = 0;
1014   unsigned int nTotalLS = 0;
1015   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
1016     nTotalMD += miniDoublets.nMDs()[idx];
1017     nTotalLS += segmentsOccupancy.nSegments()[idx];
1018   }
1019 
1020   auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx");
1021   auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx");
1022   auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx");
1023 
1024   std::set<unsigned int> lss_used_in_true_tc;
1025   unsigned int nTrackCandidates = trackCandidates.nTrackCandidates();
1026   for (unsigned int idx = 0; idx < nTrackCandidates; idx++) {
1027     // Only consider true track candidates
1028     std::vector<unsigned int> hitidxs;
1029     std::vector<unsigned int> hittypes;
1030     std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromTC(event, idx);
1031     std::vector<int> simidxs =
1032         matchedSimTrkIdxs(hitidxs, hittypes, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
1033     if (simidxs.size() == 0)
1034       continue;
1035 
1036     std::vector<unsigned int> LSs = getLSsFromTC(event, idx);
1037     for (auto& LS : LSs) {
1038       if (lss_used_in_true_tc.find(LS) == lss_used_in_true_tc.end()) {
1039         lss_used_in_true_tc.insert(LS);
1040       }
1041     }
1042   }
1043 
1044   std::cout << " lss_used_in_true_tc.size(): " << lss_used_in_true_tc.size() << std::endl;
1045 
1046   // std::cout <<  " nTotalMD: " << nTotalMD <<  std::endl;
1047   // std::cout <<  " nTotalLS: " << nTotalLS <<  std::endl;
1048 
1049   auto const& trk_sim_bunchCrossing = trk.getVI("sim_bunchCrossing");
1050   auto const& trk_sim_event = trk.getVI("sim_event");
1051   auto const& trk_sim_pt = trk.getVF("sim_pt");
1052   auto const& trk_sim_eta = trk.getVF("sim_eta");
1053   auto const& trk_sim_phi = trk.getVF("sim_phi");
1054   auto const& trk_sim_pca_dxy = trk.getVF("sim_pca_dxy");
1055   auto const& trk_sim_pca_dz = trk.getVF("sim_pca_dz");
1056   auto const& trk_sim_q = trk.getVI("sim_q");
1057   auto const& trk_sim_pdgId = trk.getVI("sim_pdgId");
1058   auto const& trk_sim_parentVtxIdx = trk.getVI("sim_parentVtxIdx");
1059   auto const& trk_simvtx_x = trk.getVF("simvtx_x");
1060   auto const& trk_simvtx_y = trk.getVF("simvtx_y");
1061   auto const& trk_simvtx_z = trk.getVF("simvtx_z");
1062 
1063   // Loop over modules (lower ones where the MDs are saved)
1064   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
1065     // // Loop over minidoublets
1066     // for (unsigned int jdx = 0; jdx < miniDoublets->nMDs[idx]; jdx++)
1067     // {
1068     //     // Get the actual index to the mini-doublet using ranges
1069     //     unsigned int mdIdx = ranges->miniDoubletModuleIndices[idx] + jdx;
1070 
1071     //     setGnnNtupleMiniDoublet(event, mdIdx);
1072     // }
1073 
1074     // Loop over segments
1075     for (unsigned int jdx = 0; jdx < segmentsOccupancy.nSegments()[idx]; jdx++) {
1076       // Get the actual index to the segments using ranges
1077       unsigned int sgIdx = ranges.segmentModuleIndices()[idx] + jdx;
1078 
1079       // Get the hit indices
1080       std::vector<unsigned int> MDs = getMDsFromLS(event, sgIdx);
1081 
1082       if (mds_used_in_sg.find(MDs[0]) == mds_used_in_sg.end()) {
1083         mds_used_in_sg.insert(MDs[0]);
1084         md_index_map[MDs[0]] = mds_used_in_sg.size() - 1;
1085         setGnnNtupleMiniDoublet(event,
1086                                 MDs[0],
1087                                 trk_sim_q,
1088                                 trk_sim_pt,
1089                                 trk_sim_eta,
1090                                 trk_sim_bunchCrossing,
1091                                 trk_sim_event,
1092                                 trk_sim_parentVtxIdx,
1093                                 trk_simvtx_x,
1094                                 trk_simvtx_y,
1095                                 trk_simvtx_z,
1096                                 trk_simhit_simTrkIdx,
1097                                 trk_ph2_simHitIdx,
1098                                 trk_pix_simHitIdx);
1099       }
1100 
1101       if (mds_used_in_sg.find(MDs[1]) == mds_used_in_sg.end()) {
1102         mds_used_in_sg.insert(MDs[1]);
1103         md_index_map[MDs[1]] = mds_used_in_sg.size() - 1;
1104         setGnnNtupleMiniDoublet(event,
1105                                 MDs[1],
1106                                 trk_sim_q,
1107                                 trk_sim_pt,
1108                                 trk_sim_eta,
1109                                 trk_sim_bunchCrossing,
1110                                 trk_sim_event,
1111                                 trk_sim_parentVtxIdx,
1112                                 trk_simvtx_x,
1113                                 trk_simvtx_y,
1114                                 trk_simvtx_z,
1115                                 trk_simhit_simTrkIdx,
1116                                 trk_ph2_simHitIdx,
1117                                 trk_pix_simHitIdx);
1118       }
1119 
1120       ana.tx->pushbackToBranch<int>("LS_MD_idx0", md_index_map[MDs[0]]);
1121       ana.tx->pushbackToBranch<int>("LS_MD_idx1", md_index_map[MDs[1]]);
1122 
1123       std::vector<unsigned int> hits = getHitsFromLS(event, sgIdx);
1124 
1125       // Computing line segment pt estimate (assuming beam spot is at zero)
1126       lst_math::Hit hitA(0, 0, 0);
1127       lst_math::Hit hitB(hitsBase.xs()[hits[0]], hitsBase.ys()[hits[0]], hitsBase.zs()[hits[0]]);
1128       lst_math::Hit hitC(hitsBase.xs()[hits[2]], hitsBase.ys()[hits[2]], hitsBase.zs()[hits[2]]);
1129       lst_math::Hit center = lst_math::getCenterFromThreePoints(hitA, hitB, hitC);
1130       float pt = lst_math::ptEstimateFromRadius(center.rt());
1131       float eta = hitC.eta();
1132       float phi = hitB.phi();
1133 
1134       ana.tx->pushbackToBranch<float>("LS_pt", pt);
1135       ana.tx->pushbackToBranch<float>("LS_eta", eta);
1136       ana.tx->pushbackToBranch<float>("LS_phi", phi);
1137       // ana.tx->pushbackToBranch<int>("LS_layer0", layer0);
1138       // ana.tx->pushbackToBranch<int>("LS_layer1", layer1);
1139 
1140       std::vector<unsigned int> hitidxs;
1141       std::vector<unsigned int> hittypes;
1142       std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromLS(event, sgIdx);
1143       std::vector<int> simidxs =
1144           matchedSimTrkIdxs(hitidxs, hittypes, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
1145 
1146       ana.tx->pushbackToBranch<int>("LS_isFake", simidxs.size() == 0);
1147       ana.tx->pushbackToBranch<float>("LS_sim_pt", simidxs.size() > 0 ? trk_sim_pt[simidxs[0]] : -999);
1148       ana.tx->pushbackToBranch<float>("LS_sim_eta", simidxs.size() > 0 ? trk_sim_eta[simidxs[0]] : -999);
1149       ana.tx->pushbackToBranch<float>("LS_sim_phi", simidxs.size() > 0 ? trk_sim_phi[simidxs[0]] : -999);
1150       ana.tx->pushbackToBranch<float>("LS_sim_pca_dxy", simidxs.size() > 0 ? trk_sim_pca_dxy[simidxs[0]] : -999);
1151       ana.tx->pushbackToBranch<float>("LS_sim_pca_dz", simidxs.size() > 0 ? trk_sim_pca_dz[simidxs[0]] : -999);
1152       ana.tx->pushbackToBranch<int>("LS_sim_q", simidxs.size() > 0 ? trk_sim_q[simidxs[0]] : -999);
1153       ana.tx->pushbackToBranch<int>("LS_sim_event", simidxs.size() > 0 ? trk_sim_event[simidxs[0]] : -999);
1154       ana.tx->pushbackToBranch<int>("LS_sim_bx", simidxs.size() > 0 ? trk_sim_bunchCrossing[simidxs[0]] : -999);
1155       ana.tx->pushbackToBranch<int>("LS_sim_pdgId", simidxs.size() > 0 ? trk_sim_pdgId[simidxs[0]] : -999);
1156       ana.tx->pushbackToBranch<float>("LS_sim_vx",
1157                                       simidxs.size() > 0 ? trk_simvtx_x[trk_sim_parentVtxIdx[simidxs[0]]] : -999);
1158       ana.tx->pushbackToBranch<float>("LS_sim_vy",
1159                                       simidxs.size() > 0 ? trk_simvtx_y[trk_sim_parentVtxIdx[simidxs[0]]] : -999);
1160       ana.tx->pushbackToBranch<float>("LS_sim_vz",
1161                                       simidxs.size() > 0 ? trk_simvtx_z[trk_sim_parentVtxIdx[simidxs[0]]] : -999);
1162       ana.tx->pushbackToBranch<int>("LS_isInTrueTC", lss_used_in_true_tc.find(sgIdx) != lss_used_in_true_tc.end());
1163 
1164       sg_index_map[sgIdx] = ana.tx->getBranch<std::vector<int>>("LS_isFake").size() - 1;
1165 
1166       // // T5 eta and phi are computed using outer and innermost hits
1167       // lst_math::Hit hitA(trk_ph2_x[anchitidx], trk_ph2_y[anchitidx], trk_ph2_z[anchitidx]);
1168       // const float phi = hitA.phi();
1169       // const float eta = hitA.eta();
1170     }
1171   }
1172 
1173   for (unsigned int idx = 0; idx < nTrackCandidates; idx++) {
1174     std::vector<unsigned int> LSs = getLSsFromTC(event, idx);
1175     std::vector<int> lsIdx;
1176     for (auto& LS : LSs) {
1177       lsIdx.push_back(sg_index_map[LS]);
1178     }
1179     ana.tx->pushbackToBranch<std::vector<int>>("tc_lsIdx", lsIdx);
1180   }
1181 
1182   std::cout << " mds_used_in_sg.size(): " << mds_used_in_sg.size() << std::endl;
1183 }
1184 
1185 //________________________________________________________________________________________________________________________________
1186 void setGnnNtupleMiniDoublet(LSTEvent* event,
1187                              unsigned int MD,
1188                              std::vector<int> const& trk_sim_q,
1189                              std::vector<float> const& trk_sim_pt,
1190                              std::vector<float> const& trk_sim_eta,
1191                              std::vector<int> const& trk_sim_bunchCrossing,
1192                              std::vector<int> const& trk_sim_event,
1193                              std::vector<int> const& trk_sim_parentVtxIdx,
1194                              std::vector<float> const& trk_simvtx_x,
1195                              std::vector<float> const& trk_simvtx_y,
1196                              std::vector<float> const& trk_simvtx_z,
1197                              std::vector<int> const& trk_simhit_simTrkIdx,
1198                              std::vector<std::vector<int>> const& trk_ph2_simHitIdx,
1199                              std::vector<std::vector<int>> const& trk_pix_simHitIdx) {
1200   // Get relevant information
1201   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1202   auto hitsBase = event->getInput<HitsBaseSoA>();
1203 
1204   // Get the hit indices
1205   unsigned int hit0 = miniDoublets.anchorHitIndices()[MD];
1206   unsigned int hit1 = miniDoublets.outerHitIndices()[MD];
1207 
1208   // Get the hit infos
1209   const float hit0_x = hitsBase.xs()[hit0];
1210   const float hit0_y = hitsBase.ys()[hit0];
1211   const float hit0_z = hitsBase.zs()[hit0];
1212   const float hit0_r = sqrt(hit0_x * hit0_x + hit0_y * hit0_y);
1213   const float hit1_x = hitsBase.xs()[hit1];
1214   const float hit1_y = hitsBase.ys()[hit1];
1215   const float hit1_z = hitsBase.zs()[hit1];
1216   const float hit1_r = sqrt(hit1_x * hit1_x + hit1_y * hit1_y);
1217 
1218   // Do sim matching
1219   std::vector<unsigned int> hit_idx = {hitsBase.idxs()[hit0], hitsBase.idxs()[hit1]};
1220   std::vector<unsigned int> hit_type = {4, 4};
1221   std::vector<int> simidxs =
1222       matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
1223 
1224   bool isFake = simidxs.size() == 0;
1225   int tp_type = getDenomSimTrkType(simidxs,
1226                                    trk_sim_q,
1227                                    trk_sim_pt,
1228                                    trk_sim_eta,
1229                                    trk_sim_bunchCrossing,
1230                                    trk_sim_event,
1231                                    trk_sim_parentVtxIdx,
1232                                    trk_simvtx_x,
1233                                    trk_simvtx_y,
1234                                    trk_simvtx_z);
1235 
1236   auto const& trk_ph2_subdet = trk.getVUS("ph2_subdet");
1237   auto const& trk_ph2_layer = trk.getVUS("ph2_layer");
1238   auto const& trk_ph2_detId = trk.getVU("ph2_detId");
1239   auto const& trk_ph2_x = trk.getVF("ph2_x");
1240   auto const& trk_ph2_y = trk.getVF("ph2_y");
1241   auto const& trk_ph2_z = trk.getVF("ph2_z");
1242 
1243   // Obtain where the actual hit is located in terms of their layer, module, rod, and ring number
1244   unsigned int anchitidx = hitsBase.idxs()[hit0];
1245   int subdet = trk_ph2_subdet[hitsBase.idxs()[anchitidx]];
1246   int is_endcap = subdet == 4;
1247   // 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)
1248   int layer = trk_ph2_layer[anchitidx] + 6 * (is_endcap);
1249   int detId = trk_ph2_detId[anchitidx];
1250 
1251   // Obtaining dPhiChange
1252   float dphichange = miniDoublets.dphichanges()[MD];
1253 
1254   // Computing pt
1255   float pt = hit0_r * k2Rinv1GeVf / sin(dphichange);
1256 
1257   // T5 eta and phi are computed using outer and innermost hits
1258   lst_math::Hit hitA(trk_ph2_x[anchitidx], trk_ph2_y[anchitidx], trk_ph2_z[anchitidx]);
1259   const float phi = hitA.phi();
1260   const float eta = hitA.eta();
1261 
1262   // Mini Doublets
1263   ana.tx->pushbackToBranch<float>("MD_pt", pt);
1264   ana.tx->pushbackToBranch<float>("MD_eta", eta);
1265   ana.tx->pushbackToBranch<float>("MD_phi", phi);
1266   ana.tx->pushbackToBranch<float>("MD_dphichange", dphichange);
1267   ana.tx->pushbackToBranch<int>("MD_isFake", isFake);
1268   ana.tx->pushbackToBranch<int>("MD_tpType", tp_type);
1269   ana.tx->pushbackToBranch<int>("MD_detId", detId);
1270   ana.tx->pushbackToBranch<int>("MD_layer", layer);
1271   ana.tx->pushbackToBranch<float>("MD_0_r", hit0_r);
1272   ana.tx->pushbackToBranch<float>("MD_0_x", hit0_x);
1273   ana.tx->pushbackToBranch<float>("MD_0_y", hit0_y);
1274   ana.tx->pushbackToBranch<float>("MD_0_z", hit0_z);
1275   ana.tx->pushbackToBranch<float>("MD_1_r", hit1_r);
1276   ana.tx->pushbackToBranch<float>("MD_1_x", hit1_x);
1277   ana.tx->pushbackToBranch<float>("MD_1_y", hit1_y);
1278   ana.tx->pushbackToBranch<float>("MD_1_z", hit1_z);
1279   // ana.tx->pushbackToBranch<int>("MD_sim_idx", simidxs.size() > 0 ? simidxs[0] : -999);
1280 }
1281 
1282 //________________________________________________________________________________________________________________________________
1283 std::tuple<int, float, float, float, int, std::vector<int>> parseTrackCandidate(
1284     LSTEvent* event,
1285     unsigned int idx,
1286     std::vector<float> const& trk_ph2_x,
1287     std::vector<float> const& trk_ph2_y,
1288     std::vector<float> const& trk_ph2_z,
1289     std::vector<int> const& trk_simhit_simTrkIdx,
1290     std::vector<std::vector<int>> const& trk_ph2_simHitIdx,
1291     std::vector<std::vector<int>> const& trk_pix_simHitIdx) {
1292   // Get the type of the track candidate
1293   auto const& trackCandidates = event->getTrackCandidates();
1294   short type = trackCandidates.trackCandidateType()[idx];
1295 
1296   // Compute pt eta phi and hit indices that will be used to figure out whether the TC matched
1297   float pt, eta, phi;
1298   std::vector<unsigned int> hit_idx, hit_type;
1299   switch (type) {
1300     case LSTObjType::pT5:
1301       std::tie(pt, eta, phi, hit_idx, hit_type) = parsepT5(event, idx);
1302       break;
1303     case LSTObjType::pT3:
1304       std::tie(pt, eta, phi, hit_idx, hit_type) = parsepT3(event, idx);
1305       break;
1306     case LSTObjType::T5:
1307       std::tie(pt, eta, phi, hit_idx, hit_type) = parseT5(event, idx, trk_ph2_x, trk_ph2_y, trk_ph2_z);
1308       break;
1309     case LSTObjType::pLS:
1310       std::tie(pt, eta, phi, hit_idx, hit_type) = parsepLS(event, idx);
1311       break;
1312   }
1313 
1314   // Perform matching
1315   std::vector<int> simidx =
1316       matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
1317   int isFake = simidx.size() == 0;
1318 
1319   return {type, pt, eta, phi, isFake, simidx};
1320 }
1321 
1322 //________________________________________________________________________________________________________________________________
1323 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parsepT5(LSTEvent* event,
1324                                                                                                unsigned int idx) {
1325   // Get relevant information
1326   auto const trackCandidates = event->getTrackCandidates();
1327   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
1328   auto const pixelSeeds = event->getInput<PixelSeedsSoA>();
1329 
1330   //
1331   // pictorial representation of a pT5
1332   //
1333   // inner tracker        outer tracker
1334   // -------------  --------------------------
1335   // pLS            01    23    45    67    89   (anchor hit of a minidoublet is always the first of the pair)
1336   // ****           oo -- oo -- oo -- oo -- oo   pT5
1337   //                oo -- oo -- oo               first T3 of the T5
1338   //                            oo -- oo -- oo   second T3 of the T5
1339   unsigned int pT5 = trackCandidates.directObjectIndices()[idx];
1340   unsigned int pLS = getPixelLSFrompT5(event, pT5);
1341   unsigned int T5Index = getT5FrompT5(event, pT5);
1342 
1343   //=================================================================================
1344   // Some history and geometry lesson...
1345   // For a given T3, we compute two angles. (NOTE: This is a bit weird!)
1346   // Historically, T3 were created out of T4, which we used to build a long time ago.
1347   // So for the sake of argument let's discuss T4 first.
1348   // For a T4, we have 4 mini-doublets.
1349   // Therefore we have 4 "anchor hits".
1350   // Therefore we have 4 xyz points.
1351   //
1352   //
1353   //       *
1354   //       |\
1355     //       | \
1356     //       |1 \
1357     //       |   \
1358     //       |  * \
1359     //       |
1360   //       |
1361   //       |
1362   //       |
1363   //       |
1364   //       |  * /
1365   //       |   /
1366   //       |2 /
1367   //       | /
1368   //       |/
1369   //       *
1370   //
1371   //
1372   // Then from these 4 points, one can approximate a some sort of "best" fitted circle trajectory,
1373   // and obtain "tangential" angles from 1st and 4th hits.
1374   // See the carton below.
1375   // The "*" are the 4 physical hit points
1376   // angle 1 and 2 are the "tangential" angle for a "circle" from 4 * points.
1377   // Please note, that a straight line from first two * and the latter two * are NOT the
1378   // angle 1 and angle 2. (they were called "beta" angles)
1379   // But rather, a slightly larger angle.
1380   // Because 4 * points would be on a circle, and a tangential line on the circles
1381   // would deviate from the points on circles.
1382   //
1383   // In the early days of LST, there was an iterative algorithm (devised by Slava) to
1384   // obtain the angle beta1 and 2 _without_ actually performing a 4 point circle fit.
1385   // Hence, the beta1 and beta2 were quickly estimated without too many math operations
1386   // and afterwards (beta1-beta2) was computed to obtain what we call a "delta-beta" values.
1387   //
1388   // For a real track, the deltabeta ~ 0, for fakes, it'd have a flat distribution.
1389   //
1390   // However, after some time we abandonded the T4s, and moved to T3s.
1391   // In T3, however, now we have the following cartoon:
1392   //
1393   //       *
1394   //       |\
1395     //       | \
1396     //       |1 \
1397     //       |   \
1398     //       |  * X   (* here are "two" MDs but really just one)
1399   //       |   /
1400   //       |2 /
1401   //       | /
1402   //       |/
1403   //       *
1404   //
1405   // With the "four" *'s (really just "three") you can still run the iterative beta calculation,
1406   // which is what we still currently do, we still get two beta1 and beta2
1407   // But! high school geometry tells us that 3 points = ONLY 1 possible CIRCLE!
1408   // There is really nothing to "fit" here.
1409   // YET we still compute these in T3, out of legacy method of how we used to treat T4s.
1410   //
1411   // Hence, in the below code, "betaIn_in" and "betaOut_in" if we performed
1412   // a circle fit they would come out by definition identical values.
1413   // But due to our approximate iterative beta calculation method, they come out different values.
1414   // So if we are "cutting on" abs(deltaBeta) = abs(betaIn_in - betaOut_in) < threshold,
1415   // what does that even mean?
1416   //
1417   // Anyhow, as of now, we compute 2 beta's for T3s, and T5 has two T3s.
1418   // And from there we estimate the pt's and we compute pt_T5.
1419 
1420   // pixel pt
1421   const float pt_pLS = pixelSeeds.ptIn()[pLS];
1422   const float eta_pLS = pixelSeeds.eta()[pLS];
1423   const float phi_pLS = pixelSeeds.phi()[pLS];
1424   float pt_T5 = __H2F(quintuplets.innerRadius()[T5Index]) * 2 * k2Rinv1GeVf;
1425   const float pt = (pt_T5 + pt_pLS) / 2;
1426 
1427   // Form the hit idx/type std::vector
1428   std::vector<unsigned int> hit_idx = getHitIdxsFrompT5(event, pT5);
1429   std::vector<unsigned int> hit_type = getHitTypesFrompT5(event, pT5);
1430 
1431   return {pt, eta_pLS, phi_pLS, hit_idx, hit_type};
1432 }
1433 
1434 //________________________________________________________________________________________________________________________________
1435 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parsepT3(LSTEvent* event,
1436                                                                                                unsigned int idx) {
1437   // Get relevant information
1438   auto const trackCandidates = event->getTrackCandidates();
1439   auto const triplets = event->getTriplets<TripletsSoA>();
1440   auto const pixelSeeds = event->getInput<PixelSeedsSoA>();
1441 
1442   //
1443   // pictorial representation of a pT3
1444   //
1445   // inner tracker        outer tracker
1446   // -------------  --------------------------
1447   // pLS            01    23    45               (anchor hit of a minidoublet is always the first of the pair)
1448   // ****           oo -- oo -- oo               pT3
1449   unsigned int pT3 = trackCandidates.directObjectIndices()[idx];
1450   unsigned int pLS = getPixelLSFrompT3(event, pT3);
1451   unsigned int T3 = getT3FrompT3(event, pT3);
1452 
1453   // pixel pt
1454   const float pt_pLS = pixelSeeds.ptIn()[pLS];
1455   const float eta_pLS = pixelSeeds.eta()[pLS];
1456   const float phi_pLS = pixelSeeds.phi()[pLS];
1457   float pt_T3 = triplets.radius()[T3] * 2 * k2Rinv1GeVf;
1458 
1459   // average pt
1460   const float pt = (pt_pLS + pt_T3) / 2;
1461 
1462   // Form the hit idx/type std::vector
1463   std::vector<unsigned int> hit_idx = getHitIdxsFrompT3(event, pT3);
1464   std::vector<unsigned int> hit_type = getHitTypesFrompT3(event, pT3);
1465 
1466   return {pt, eta_pLS, phi_pLS, hit_idx, hit_type};
1467 }
1468 
1469 //________________________________________________________________________________________________________________________________
1470 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parseT5(
1471     LSTEvent* event,
1472     unsigned int idx,
1473     std::vector<float> const& trk_ph2_x,
1474     std::vector<float> const& trk_ph2_y,
1475     std::vector<float> const& trk_ph2_z) {
1476   auto const trackCandidates = event->getTrackCandidates();
1477   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
1478   unsigned int T5 = trackCandidates.directObjectIndices()[idx];
1479   std::vector<unsigned int> hits = getHitsFromT5(event, T5);
1480 
1481   //
1482   // pictorial representation of a T5
1483   //
1484   // inner tracker        outer tracker
1485   // -------------  --------------------------
1486   //                01    23    45    67    89   (anchor hit of a minidoublet is always the first of the pair)
1487   //  (none)        oo -- oo -- oo -- oo -- oo   T5
1488   unsigned int Hit_0 = hits[0];
1489   unsigned int Hit_4 = hits[4];
1490   unsigned int Hit_8 = hits[8];
1491 
1492   // T5 radius is average of the inner and outer radius
1493   const float pt = __H2F(quintuplets.innerRadius()[T5]) * k2Rinv1GeVf * 2;
1494 
1495   // T5 eta and phi are computed using outer and innermost hits
1496   lst_math::Hit hitA(trk_ph2_x[Hit_0], trk_ph2_y[Hit_0], trk_ph2_z[Hit_0]);
1497   lst_math::Hit hitB(trk_ph2_x[Hit_8], trk_ph2_y[Hit_8], trk_ph2_z[Hit_8]);
1498   const float phi = hitA.phi();
1499   const float eta = hitB.eta();
1500 
1501   std::vector<unsigned int> hit_idx = getHitIdxsFromT5(event, T5);
1502   std::vector<unsigned int> hit_type = getHitTypesFromT5(event, T5);
1503 
1504   return {pt, eta, phi, hit_idx, hit_type};
1505 }
1506 
1507 //________________________________________________________________________________________________________________________________
1508 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parsepLS(LSTEvent* event,
1509                                                                                                unsigned int idx) {
1510   auto const& trackCandidates = event->getTrackCandidates();
1511   auto pixelSeeds = event->getInput<PixelSeedsSoA>();
1512 
1513   // Getting pLS index
1514   unsigned int pLS = trackCandidates.directObjectIndices()[idx];
1515 
1516   // Getting pt eta and phi
1517   float pt = pixelSeeds.ptIn()[pLS];
1518   float eta = pixelSeeds.eta()[pLS];
1519   float phi = pixelSeeds.phi()[pLS];
1520 
1521   // Getting hit indices and types
1522   std::vector<unsigned int> hit_idx = getPixelHitIdxsFrompLS(event, pLS);
1523   std::vector<unsigned int> hit_type = getPixelHitTypesFrompLS(event, pLS);
1524 
1525   return {pt, eta, phi, hit_idx, hit_type};
1526 }
1527 
1528 void setpLSOutputBranches(LSTEvent* event) {
1529   auto const& pixelSegments = event->getPixelSegments();
1530   auto const& pixelSeeds = event->getInput<PixelSeedsSoA>();
1531   int n_accepted_simtrk = ana.tx->getBranch<std::vector<int>>("sim_TC_matched").size();
1532 
1533   auto const& trk_simhit_simTrkIdx = trk.getVI("simhit_simTrkIdx");
1534   auto const& trk_ph2_simHitIdx = trk.getVVI("ph2_simHitIdx");
1535   auto const& trk_pix_simHitIdx = trk.getVVI("pix_simHitIdx");
1536 
1537   unsigned int n_pLS = pixelSegments.metadata().size();
1538   std::vector<int> sim_pLS_matched(n_accepted_simtrk, 0);
1539   std::vector<std::vector<int>> pLS_matched_simIdx;
1540 
1541   for (unsigned int i_pLS = 0; i_pLS < n_pLS; ++i_pLS) {
1542     // Get pLS properties
1543     float pt = pixelSeeds.ptIn()[i_pLS];
1544     float px = pixelSeeds.px()[i_pLS];
1545     float py = pixelSeeds.py()[i_pLS];
1546     float pz = pixelSeeds.pz()[i_pLS];
1547     bool isQuad = static_cast<bool>(pixelSeeds.isQuad()[i_pLS]);
1548     float ptErr = pixelSeeds.ptErr()[i_pLS];
1549     float eta = pixelSeeds.eta()[i_pLS];
1550     float etaErr = pixelSeeds.etaErr()[i_pLS];
1551     float phi = pixelSeeds.phi()[i_pLS];
1552     float score = pixelSegments.score()[i_pLS];
1553     float centerX = pixelSegments.circleCenterX()[i_pLS];
1554     float centerY = pixelSegments.circleCenterY()[i_pLS];
1555     float radius = pixelSegments.circleRadius()[i_pLS];
1556 
1557     // Get hits from pLS
1558     std::vector<unsigned int> hit_idx = getPixelHitIdxsFrompLS(event, i_pLS);
1559     std::vector<unsigned int> hit_type = getPixelHitTypesFrompLS(event, i_pLS);
1560 
1561     // Match to sim tracks
1562     std::vector<int> simidx =
1563         matchedSimTrkIdxs(hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx);
1564     bool isFake = simidx.empty();
1565 
1566     // Fill branches
1567     ana.tx->pushbackToBranch<float>("pLS_ptIn", pt);
1568     ana.tx->pushbackToBranch<float>("pLS_ptErr", ptErr);
1569     ana.tx->pushbackToBranch<float>("pLS_px", px);
1570     ana.tx->pushbackToBranch<float>("pLS_py", py);
1571     ana.tx->pushbackToBranch<float>("pLS_pz", pz);
1572     ana.tx->pushbackToBranch<float>("pLS_eta", eta);
1573     ana.tx->pushbackToBranch<float>("pLS_etaErr", etaErr);
1574     ana.tx->pushbackToBranch<float>("pLS_phi", phi);
1575     ana.tx->pushbackToBranch<float>("pLS_score", score);
1576     ana.tx->pushbackToBranch<float>("pLS_circleCenterX", centerX);
1577     ana.tx->pushbackToBranch<float>("pLS_circleCenterY", centerY);
1578     ana.tx->pushbackToBranch<float>("pLS_circleRadius", radius);
1579     ana.tx->pushbackToBranch<bool>("pLS_isQuad", isQuad);
1580     ana.tx->pushbackToBranch<int>("pLS_isFake", isFake);
1581     pLS_matched_simIdx.push_back(simidx);
1582 
1583     // Count matches
1584     for (auto& idx : simidx) {
1585       if (idx < n_accepted_simtrk) {
1586         sim_pLS_matched[idx]++;
1587       }
1588     }
1589   }
1590 
1591   std::vector<int> pLS_isDuplicate(pLS_matched_simIdx.size(), 0);
1592   for (size_t i = 0; i < pLS_matched_simIdx.size(); ++i) {
1593     for (int simidx : pLS_matched_simIdx[i]) {
1594       if (simidx < n_accepted_simtrk && sim_pLS_matched[simidx] > 1) {
1595         pLS_isDuplicate[i] = 1;
1596         break;
1597       }
1598     }
1599   }
1600 
1601   ana.tx->setBranch<std::vector<int>>("sim_pLS_matched", sim_pLS_matched);
1602   ana.tx->setBranch<std::vector<std::vector<int>>>("pLS_matched_simIdx", pLS_matched_simIdx);
1603   ana.tx->setBranch<std::vector<int>>("pLS_isDuplicate", pLS_isDuplicate);
1604 }
1605 
1606 //________________________________________________________________________________________________________________________________
1607 void printHitMultiplicities(LSTEvent* event) {
1608   auto modules = event->getModules<ModulesSoA>();
1609   auto hitRanges = event->getHits<HitsRangesSoA>();
1610 
1611   int nHits = 0;
1612   for (unsigned int idx = 0; idx <= modules.nLowerModules();
1613        idx++)  // "<=" because cheating to include pixel track candidate lower module
1614   {
1615     nHits += hitRanges.hitRanges()[2 * idx][1] - hitRanges.hitRanges()[2 * idx][0] + 1;
1616     nHits += hitRanges.hitRanges()[2 * idx + 1][1] - hitRanges.hitRanges()[2 * idx + 1][0] + 1;
1617   }
1618   std::cout << " nHits: " << nHits << std::endl;
1619 }
1620 
1621 //________________________________________________________________________________________________________________________________
1622 void printMiniDoubletMultiplicities(LSTEvent* event) {
1623   MiniDoubletsOccupancyConst miniDoublets = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
1624   auto modules = event->getModules<ModulesSoA>();
1625 
1626   int nMiniDoublets = 0;
1627   int totOccupancyMiniDoublets = 0;
1628   for (unsigned int idx = 0; idx <= modules.nModules();
1629        idx++)  // "<=" because cheating to include pixel track candidate lower module
1630   {
1631     if (modules.isLower()[idx]) {
1632       nMiniDoublets += miniDoublets.nMDs()[idx];
1633       totOccupancyMiniDoublets += miniDoublets.totOccupancyMDs()[idx];
1634     }
1635   }
1636   std::cout << " nMiniDoublets: " << nMiniDoublets << std::endl;
1637   std::cout << " totOccupancyMiniDoublets (including trucated ones): " << totOccupancyMiniDoublets << std::endl;
1638 }
1639 
1640 //________________________________________________________________________________________________________________________________
1641 void printAllObjects(LSTEvent* event) {
1642   printMDs(event);
1643   printLSs(event);
1644   printpLSs(event);
1645   printT3s(event);
1646 }
1647 
1648 //________________________________________________________________________________________________________________________________
1649 void printMDs(LSTEvent* event) {
1650   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1651   MiniDoubletsOccupancyConst miniDoubletsOccupancy = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
1652   auto hitsBase = event->getInput<HitsBaseSoA>();
1653   auto modules = event->getModules<ModulesSoA>();
1654   auto ranges = event->getRanges();
1655 
1656   // Then obtain the lower module index
1657   for (unsigned int idx = 0; idx <= modules.nLowerModules(); ++idx) {
1658     for (unsigned int iMD = 0; iMD < miniDoubletsOccupancy.nMDs()[idx]; iMD++) {
1659       unsigned int mdIdx = ranges.miniDoubletModuleIndices()[idx] + iMD;
1660       unsigned int LowerHitIndex = miniDoublets.anchorHitIndices()[mdIdx];
1661       unsigned int UpperHitIndex = miniDoublets.outerHitIndices()[mdIdx];
1662       unsigned int hit0 = hitsBase.idxs()[LowerHitIndex];
1663       unsigned int hit1 = hitsBase.idxs()[UpperHitIndex];
1664       std::cout << "VALIDATION 'MD': "
1665                 << "MD"
1666                 << " hit0: " << hit0 << " hit1: " << hit1 << std::endl;
1667     }
1668   }
1669 }
1670 
1671 //________________________________________________________________________________________________________________________________
1672 void printLSs(LSTEvent* event) {
1673   SegmentsConst segments = event->getSegments<SegmentsSoA>();
1674   SegmentsOccupancyConst segmentsOccupancy = event->getSegments<SegmentsOccupancySoA>();
1675   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1676   auto hitsBase = event->getInput<HitsBaseSoA>();
1677   auto modules = event->getModules<ModulesSoA>();
1678   auto ranges = event->getRanges();
1679 
1680   int nSegments = 0;
1681   for (unsigned int i = 0; i < modules.nLowerModules(); ++i) {
1682     unsigned int idx = i;  //modules->lowerModuleIndices[i];
1683     nSegments += segmentsOccupancy.nSegments()[idx];
1684     for (unsigned int jdx = 0; jdx < segmentsOccupancy.nSegments()[idx]; jdx++) {
1685       unsigned int sgIdx = ranges.segmentModuleIndices()[idx] + jdx;
1686       unsigned int InnerMiniDoubletIndex = segments.mdIndices()[sgIdx][0];
1687       unsigned int OuterMiniDoubletIndex = segments.mdIndices()[sgIdx][1];
1688       unsigned int InnerMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[InnerMiniDoubletIndex];
1689       unsigned int InnerMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[InnerMiniDoubletIndex];
1690       unsigned int OuterMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[OuterMiniDoubletIndex];
1691       unsigned int OuterMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[OuterMiniDoubletIndex];
1692       unsigned int hit0 = hitsBase.idxs()[InnerMiniDoubletLowerHitIndex];
1693       unsigned int hit1 = hitsBase.idxs()[InnerMiniDoubletUpperHitIndex];
1694       unsigned int hit2 = hitsBase.idxs()[OuterMiniDoubletLowerHitIndex];
1695       unsigned int hit3 = hitsBase.idxs()[OuterMiniDoubletUpperHitIndex];
1696       std::cout << "VALIDATION 'LS': "
1697                 << "LS"
1698                 << " hit0: " << hit0 << " hit1: " << hit1 << " hit2: " << hit2 << " hit3: " << hit3 << std::endl;
1699     }
1700   }
1701   std::cout << "VALIDATION nSegments: " << nSegments << std::endl;
1702 }
1703 
1704 //________________________________________________________________________________________________________________________________
1705 void printpLSs(LSTEvent* event) {
1706   SegmentsConst segments = event->getSegments<SegmentsSoA>();
1707   SegmentsOccupancyConst segmentsOccupancy = event->getSegments<SegmentsOccupancySoA>();
1708   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1709   auto hitsBase = event->getInput<HitsBaseSoA>();
1710   auto modules = event->getModules<ModulesSoA>();
1711   auto ranges = event->getRanges();
1712 
1713   unsigned int i = modules.nLowerModules();
1714   unsigned int idx = i;  //modules->lowerModuleIndices[i];
1715   int npLS = segmentsOccupancy.nSegments()[idx];
1716   for (unsigned int jdx = 0; jdx < segmentsOccupancy.nSegments()[idx]; jdx++) {
1717     unsigned int sgIdx = ranges.segmentModuleIndices()[idx] + jdx;
1718     unsigned int InnerMiniDoubletIndex = segments.mdIndices()[sgIdx][0];
1719     unsigned int OuterMiniDoubletIndex = segments.mdIndices()[sgIdx][1];
1720     unsigned int InnerMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[InnerMiniDoubletIndex];
1721     unsigned int InnerMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[InnerMiniDoubletIndex];
1722     unsigned int OuterMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[OuterMiniDoubletIndex];
1723     unsigned int OuterMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[OuterMiniDoubletIndex];
1724     unsigned int hit0 = hitsBase.idxs()[InnerMiniDoubletLowerHitIndex];
1725     unsigned int hit1 = hitsBase.idxs()[InnerMiniDoubletUpperHitIndex];
1726     unsigned int hit2 = hitsBase.idxs()[OuterMiniDoubletLowerHitIndex];
1727     unsigned int hit3 = hitsBase.idxs()[OuterMiniDoubletUpperHitIndex];
1728     std::cout << "VALIDATION 'pLS': "
1729               << "pLS"
1730               << " hit0: " << hit0 << " hit1: " << hit1 << " hit2: " << hit2 << " hit3: " << hit3 << std::endl;
1731   }
1732   std::cout << "VALIDATION npLS: " << npLS << std::endl;
1733 }
1734 
1735 //________________________________________________________________________________________________________________________________
1736 void printT3s(LSTEvent* event) {
1737   auto const triplets = event->getTriplets<TripletsSoA>();
1738   auto const tripletsOccupancy = event->getTriplets<TripletsOccupancySoA>();
1739   SegmentsConst segments = event->getSegments<SegmentsSoA>();
1740   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1741   auto hitsBase = event->getInput<HitsBaseSoA>();
1742   auto modules = event->getModules<ModulesSoA>();
1743   int nTriplets = 0;
1744   for (unsigned int i = 0; i < modules.nLowerModules(); ++i) {
1745     // unsigned int idx = modules->lowerModuleIndices[i];
1746     nTriplets += tripletsOccupancy.nTriplets()[i];
1747     unsigned int idx = i;
1748     for (unsigned int jdx = 0; jdx < tripletsOccupancy.nTriplets()[idx]; jdx++) {
1749       unsigned int tpIdx = idx * 5000 + jdx;
1750       unsigned int InnerSegmentIndex = triplets.segmentIndices()[tpIdx][0];
1751       unsigned int OuterSegmentIndex = triplets.segmentIndices()[tpIdx][1];
1752       unsigned int InnerSegmentInnerMiniDoubletIndex = segments.mdIndices()[InnerSegmentIndex][0];
1753       unsigned int InnerSegmentOuterMiniDoubletIndex = segments.mdIndices()[InnerSegmentIndex][1];
1754       unsigned int OuterSegmentOuterMiniDoubletIndex = segments.mdIndices()[OuterSegmentIndex][1];
1755 
1756       unsigned int hit_idx0 = miniDoublets.anchorHitIndices()[InnerSegmentInnerMiniDoubletIndex];
1757       unsigned int hit_idx1 = miniDoublets.outerHitIndices()[InnerSegmentInnerMiniDoubletIndex];
1758       unsigned int hit_idx2 = miniDoublets.anchorHitIndices()[InnerSegmentOuterMiniDoubletIndex];
1759       unsigned int hit_idx3 = miniDoublets.outerHitIndices()[InnerSegmentOuterMiniDoubletIndex];
1760       unsigned int hit_idx4 = miniDoublets.anchorHitIndices()[OuterSegmentOuterMiniDoubletIndex];
1761       unsigned int hit_idx5 = miniDoublets.outerHitIndices()[OuterSegmentOuterMiniDoubletIndex];
1762 
1763       unsigned int hit0 = hitsBase.idxs()[hit_idx0];
1764       unsigned int hit1 = hitsBase.idxs()[hit_idx1];
1765       unsigned int hit2 = hitsBase.idxs()[hit_idx2];
1766       unsigned int hit3 = hitsBase.idxs()[hit_idx3];
1767       unsigned int hit4 = hitsBase.idxs()[hit_idx4];
1768       unsigned int hit5 = hitsBase.idxs()[hit_idx5];
1769       std::cout << "VALIDATION 'T3': "
1770                 << "T3"
1771                 << " hit0: " << hit0 << " hit1: " << hit1 << " hit2: " << hit2 << " hit3: " << hit3 << " hit4: " << hit4
1772                 << " hit5: " << hit5 << std::endl;
1773     }
1774   }
1775   std::cout << "VALIDATION nTriplets: " << nTriplets << std::endl;
1776 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!