Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-12-05 02:48:11 

0001 // to use computeRadiusFromThreeAnchorHits
0002 #include "LSTEvent.h"
0003 #include "Triplet.h"
0004 
0005 #include "write_lst_ntuple.h"
0006 
0007 using namespace ALPAKA_ACCELERATOR_NAMESPACE::lst;
0008 
0009 //________________________________________________________________________________________________________________________________
0010 void createOutputBranches() {
0011   createRequiredOutputBranches();
0012   createOptionalOutputBranches();
0013 }
0014 
0015 //________________________________________________________________________________________________________________________________
0016 void fillOutputBranches(LSTEvent* event) {
0017   setOutputBranches(event);
0018   setOptionalOutputBranches(event);
0019   if (ana.gnn_ntuple)
0020     setGnnNtupleBranches(event);
0021 
0022   // Now actually fill the ttree
0023   ana.tx->fill();
0024 
0025   // Then clear the branches to default values (e.g. -999, or clear the vectors to empty vectors)
0026   ana.tx->clear();
0027 }
0028 
0029 //________________________________________________________________________________________________________________________________
0030 void createRequiredOutputBranches() {
0031   // Setup output TTree
0032   ana.tx->createBranch<std::vector<float>>("sim_pt");
0033   ana.tx->createBranch<std::vector<float>>("sim_eta");
0034   ana.tx->createBranch<std::vector<float>>("sim_phi");
0035   ana.tx->createBranch<std::vector<float>>("sim_pca_dxy");
0036   ana.tx->createBranch<std::vector<float>>("sim_pca_dz");
0037   ana.tx->createBranch<std::vector<int>>("sim_q");
0038   ana.tx->createBranch<std::vector<int>>("sim_event");
0039   ana.tx->createBranch<std::vector<int>>("sim_pdgId");
0040   ana.tx->createBranch<std::vector<float>>("sim_vx");
0041   ana.tx->createBranch<std::vector<float>>("sim_vy");
0042   ana.tx->createBranch<std::vector<float>>("sim_vz");
0043   ana.tx->createBranch<std::vector<float>>("sim_trkNtupIdx");
0044   ana.tx->createBranch<std::vector<int>>("sim_TC_matched");
0045   ana.tx->createBranch<std::vector<int>>("sim_TC_matched_mask");
0046 
0047   // Track candidates
0048   ana.tx->createBranch<std::vector<float>>("tc_pt");
0049   ana.tx->createBranch<std::vector<float>>("tc_eta");
0050   ana.tx->createBranch<std::vector<float>>("tc_phi");
0051   ana.tx->createBranch<std::vector<int>>("tc_type");
0052   ana.tx->createBranch<std::vector<int>>("tc_isFake");
0053   ana.tx->createBranch<std::vector<int>>("tc_isDuplicate");
0054   ana.tx->createBranch<std::vector<std::vector<int>>>("tc_matched_simIdx");
0055 }
0056 
0057 //________________________________________________________________________________________________________________________________
0058 void createOptionalOutputBranches() {
0059 #ifdef CUT_VALUE_DEBUG
0060   // Event-wide branches
0061   // ana.tx->createBranch<float>("evt_dummy");
0062 
0063   // Sim Track branches
0064   // NOTE: Must sync with main tc branch in length!!
0065   ana.tx->createBranch<std::vector<float>>("sim_dummy");
0066 
0067   // Track Candidate branches
0068   // NOTE: Must sync with main tc branch in length!!
0069   ana.tx->createBranch<std::vector<float>>("tc_dummy");
0070 
0071   // pT5 branches
0072   ana.tx->createBranch<std::vector<std::vector<int>>>("pT5_matched_simIdx");
0073   ana.tx->createBranch<std::vector<std::vector<int>>>("pT5_hitIdxs");
0074   ana.tx->createBranch<std::vector<int>>("sim_pT5_matched");
0075   ana.tx->createBranch<std::vector<float>>("pT5_pt");
0076   ana.tx->createBranch<std::vector<float>>("pT5_eta");
0077   ana.tx->createBranch<std::vector<float>>("pT5_phi");
0078   ana.tx->createBranch<std::vector<int>>("pT5_isFake");
0079   ana.tx->createBranch<std::vector<float>>("t5_sim_vxy");
0080   ana.tx->createBranch<std::vector<float>>("t5_sim_vz");
0081   ana.tx->createBranch<std::vector<int>>("pT5_isDuplicate");
0082   ana.tx->createBranch<std::vector<int>>("pT5_score");
0083   ana.tx->createBranch<std::vector<int>>("pT5_layer_binary");
0084   ana.tx->createBranch<std::vector<int>>("pT5_moduleType_binary");
0085   ana.tx->createBranch<std::vector<float>>("pT5_matched_pt");
0086   ana.tx->createBranch<std::vector<float>>("pT5_rzChiSquared");
0087   ana.tx->createBranch<std::vector<float>>("pT5_rPhiChiSquared");
0088   ana.tx->createBranch<std::vector<float>>("pT5_rPhiChiSquaredInwards");
0089 
0090   // pT3 branches
0091   ana.tx->createBranch<std::vector<int>>("sim_pT3_matched");
0092   ana.tx->createBranch<std::vector<float>>("pT3_pt");
0093   ana.tx->createBranch<std::vector<int>>("pT3_isFake");
0094   ana.tx->createBranch<std::vector<int>>("pT3_isDuplicate");
0095   ana.tx->createBranch<std::vector<float>>("pT3_eta");
0096   ana.tx->createBranch<std::vector<float>>("pT3_phi");
0097   ana.tx->createBranch<std::vector<float>>("pT3_score");
0098   ana.tx->createBranch<std::vector<int>>("pT3_foundDuplicate");
0099   ana.tx->createBranch<std::vector<std::vector<int>>>("pT3_matched_simIdx");
0100   ana.tx->createBranch<std::vector<std::vector<int>>>("pT3_hitIdxs");
0101   ana.tx->createBranch<std::vector<float>>("pT3_pixelRadius");
0102   ana.tx->createBranch<std::vector<float>>("pT3_pixelRadiusError");
0103   ana.tx->createBranch<std::vector<std::vector<float>>>("pT3_matched_pt");
0104   ana.tx->createBranch<std::vector<float>>("pT3_tripletRadius");
0105   ana.tx->createBranch<std::vector<float>>("pT3_rPhiChiSquared");
0106   ana.tx->createBranch<std::vector<float>>("pT3_rPhiChiSquaredInwards");
0107   ana.tx->createBranch<std::vector<float>>("pT3_rzChiSquared");
0108   ana.tx->createBranch<std::vector<int>>("pT3_layer_binary");
0109   ana.tx->createBranch<std::vector<int>>("pT3_moduleType_binary");
0110 
0111   // pLS branches
0112   ana.tx->createBranch<std::vector<int>>("sim_pLS_matched");
0113   ana.tx->createBranch<std::vector<std::vector<int>>>("sim_pLS_types");
0114   ana.tx->createBranch<std::vector<int>>("pLS_isFake");
0115   ana.tx->createBranch<std::vector<int>>("pLS_isDuplicate");
0116   ana.tx->createBranch<std::vector<float>>("pLS_pt");
0117   ana.tx->createBranch<std::vector<float>>("pLS_eta");
0118   ana.tx->createBranch<std::vector<float>>("pLS_phi");
0119   ana.tx->createBranch<std::vector<float>>("pLS_score");
0120 
0121   // T5 branches
0122   ana.tx->createBranch<std::vector<int>>("sim_T5_matched");
0123   ana.tx->createBranch<std::vector<int>>("t5_isFake");
0124   ana.tx->createBranch<std::vector<int>>("t5_isDuplicate");
0125   ana.tx->createBranch<std::vector<int>>("t5_foundDuplicate");
0126   ana.tx->createBranch<std::vector<float>>("t5_pt");
0127   ana.tx->createBranch<std::vector<float>>("t5_pMatched");
0128   ana.tx->createBranch<std::vector<float>>("t5_eta");
0129   ana.tx->createBranch<std::vector<float>>("t5_phi");
0130   ana.tx->createBranch<std::vector<float>>("t5_score_rphisum");
0131   ana.tx->createBranch<std::vector<std::vector<int>>>("t5_hitIdxs");
0132   ana.tx->createBranch<std::vector<std::vector<int>>>("t5_matched_simIdx");
0133   ana.tx->createBranch<std::vector<int>>("t5_moduleType_binary");
0134   ana.tx->createBranch<std::vector<int>>("t5_layer_binary");
0135   ana.tx->createBranch<std::vector<float>>("t5_matched_pt");
0136   ana.tx->createBranch<std::vector<float>>("t5_innerRadius");
0137   ana.tx->createBranch<std::vector<float>>("t5_outerRadius");
0138   ana.tx->createBranch<std::vector<float>>("t5_bridgeRadius");
0139   ana.tx->createBranch<std::vector<float>>("t5_chiSquared");
0140   ana.tx->createBranch<std::vector<float>>("t5_rzChiSquared");
0141   ana.tx->createBranch<std::vector<float>>("t5_nonAnchorChiSquared");
0142   ana.tx->createBranch<std::vector<float>>("t5_dBeta1");
0143   ana.tx->createBranch<std::vector<float>>("t5_dBeta2");
0144 
0145   // Occupancy branches
0146   ana.tx->createBranch<std::vector<int>>("module_layers");
0147   ana.tx->createBranch<std::vector<int>>("module_subdets");
0148   ana.tx->createBranch<std::vector<int>>("module_rings");
0149   ana.tx->createBranch<std::vector<int>>("module_rods");
0150   ana.tx->createBranch<std::vector<int>>("module_modules");
0151   ana.tx->createBranch<std::vector<bool>>("module_isTilted");
0152   ana.tx->createBranch<std::vector<float>>("module_eta");
0153   ana.tx->createBranch<std::vector<float>>("module_r");
0154   ana.tx->createBranch<std::vector<int>>("md_occupancies");
0155   ana.tx->createBranch<std::vector<int>>("sg_occupancies");
0156   ana.tx->createBranch<std::vector<int>>("t3_occupancies");
0157   ana.tx->createBranch<int>("tc_occupancies");
0158   ana.tx->createBranch<std::vector<int>>("t5_occupancies");
0159   ana.tx->createBranch<int>("pT3_occupancies");
0160   ana.tx->createBranch<int>("pT5_occupancies");
0161 
0162   // T5 DNN branches
0163   createT5DNNBranches();
0164 
0165 #endif
0166 }
0167 
0168 //________________________________________________________________________________________________________________________________
0169 void createT5DNNBranches() {
0170   // Common branches
0171   ana.tx->createBranch<std::vector<int>>("t5_t3_idx0");
0172   ana.tx->createBranch<std::vector<int>>("t5_t3_idx1");
0173   ana.tx->createBranch<std::vector<int>>("t5_tc_idx");
0174   ana.tx->createBranch<std::vector<int>>("t5_partOfTC");
0175   ana.tx->createBranch<std::vector<float>>("t5_t3_pt");
0176   ana.tx->createBranch<std::vector<float>>("t5_t3_eta");
0177   ana.tx->createBranch<std::vector<float>>("t5_t3_phi");
0178 
0179   // Hit-specific branches
0180   std::vector<std::string> hitIndices = {"0", "1", "2", "3", "4", "5"};
0181   std::vector<std::string> hitProperties = {"r", "x", "y", "z", "eta", "phi", "detId", "layer", "moduleType"};
0182 
0183   for (const auto& idx : hitIndices) {
0184     for (const auto& prop : hitProperties) {
0185       std::string branchName = "t5_t3_" + idx + "_" + prop;
0186       if (prop == "detId" || prop == "layer" || prop == "moduleType") {
0187         ana.tx->createBranch<std::vector<int>>(branchName);
0188       } else {
0189         ana.tx->createBranch<std::vector<float>>(branchName);
0190       }
0191     }
0192   }
0193 }
0194 
0195 //________________________________________________________________________________________________________________________________
0196 void createGnnNtupleBranches() {
0197   // Mini Doublets
0198   ana.tx->createBranch<std::vector<float>>("MD_pt");
0199   ana.tx->createBranch<std::vector<float>>("MD_eta");
0200   ana.tx->createBranch<std::vector<float>>("MD_phi");
0201   ana.tx->createBranch<std::vector<float>>("MD_dphichange");
0202   ana.tx->createBranch<std::vector<int>>("MD_isFake");
0203   ana.tx->createBranch<std::vector<int>>("MD_tpType");
0204   ana.tx->createBranch<std::vector<int>>("MD_detId");
0205   ana.tx->createBranch<std::vector<int>>("MD_layer");
0206   ana.tx->createBranch<std::vector<float>>("MD_0_r");
0207   ana.tx->createBranch<std::vector<float>>("MD_0_x");
0208   ana.tx->createBranch<std::vector<float>>("MD_0_y");
0209   ana.tx->createBranch<std::vector<float>>("MD_0_z");
0210   ana.tx->createBranch<std::vector<float>>("MD_1_r");
0211   ana.tx->createBranch<std::vector<float>>("MD_1_x");
0212   ana.tx->createBranch<std::vector<float>>("MD_1_y");
0213   ana.tx->createBranch<std::vector<float>>("MD_1_z");
0214 
0215   // Line Segments
0216   ana.tx->createBranch<std::vector<float>>("LS_pt");
0217   ana.tx->createBranch<std::vector<float>>("LS_eta");
0218   ana.tx->createBranch<std::vector<float>>("LS_phi");
0219   ana.tx->createBranch<std::vector<int>>("LS_isFake");
0220   ana.tx->createBranch<std::vector<int>>("LS_MD_idx0");
0221   ana.tx->createBranch<std::vector<int>>("LS_MD_idx1");
0222   ana.tx->createBranch<std::vector<float>>("LS_sim_pt");
0223   ana.tx->createBranch<std::vector<float>>("LS_sim_eta");
0224   ana.tx->createBranch<std::vector<float>>("LS_sim_phi");
0225   ana.tx->createBranch<std::vector<float>>("LS_sim_pca_dxy");
0226   ana.tx->createBranch<std::vector<float>>("LS_sim_pca_dz");
0227   ana.tx->createBranch<std::vector<int>>("LS_sim_q");
0228   ana.tx->createBranch<std::vector<int>>("LS_sim_pdgId");
0229   ana.tx->createBranch<std::vector<int>>("LS_sim_event");
0230   ana.tx->createBranch<std::vector<int>>("LS_sim_bx");
0231   ana.tx->createBranch<std::vector<float>>("LS_sim_vx");
0232   ana.tx->createBranch<std::vector<float>>("LS_sim_vy");
0233   ana.tx->createBranch<std::vector<float>>("LS_sim_vz");
0234   ana.tx->createBranch<std::vector<int>>("LS_isInTrueTC");
0235 
0236   // TC's LS
0237   ana.tx->createBranch<std::vector<std::vector<int>>>("tc_lsIdx");
0238 }
0239 
0240 //________________________________________________________________________________________________________________________________
0241 void setOutputBranches(LSTEvent* event) {
0242   // ============ Sim tracks =============
0243   int n_accepted_simtrk = 0;
0244   for (unsigned int isimtrk = 0; isimtrk < trk.sim_pt().size(); ++isimtrk) {
0245     // Skip out-of-time pileup
0246     if (trk.sim_bunchCrossing()[isimtrk] != 0)
0247       continue;
0248 
0249     // Skip non-hard-scatter
0250     if (trk.sim_event()[isimtrk] != 0)
0251       continue;
0252 
0253     ana.tx->pushbackToBranch<float>("sim_pt", trk.sim_pt()[isimtrk]);
0254     ana.tx->pushbackToBranch<float>("sim_eta", trk.sim_eta()[isimtrk]);
0255     ana.tx->pushbackToBranch<float>("sim_phi", trk.sim_phi()[isimtrk]);
0256     ana.tx->pushbackToBranch<float>("sim_pca_dxy", trk.sim_pca_dxy()[isimtrk]);
0257     ana.tx->pushbackToBranch<float>("sim_pca_dz", trk.sim_pca_dz()[isimtrk]);
0258     ana.tx->pushbackToBranch<int>("sim_q", trk.sim_q()[isimtrk]);
0259     ana.tx->pushbackToBranch<int>("sim_event", trk.sim_event()[isimtrk]);
0260     ana.tx->pushbackToBranch<int>("sim_pdgId", trk.sim_pdgId()[isimtrk]);
0261 
0262     // For vertex we need to look it up from simvtx info
0263     int vtxidx = trk.sim_parentVtxIdx()[isimtrk];
0264     ana.tx->pushbackToBranch<float>("sim_vx", trk.simvtx_x()[vtxidx]);
0265     ana.tx->pushbackToBranch<float>("sim_vy", trk.simvtx_y()[vtxidx]);
0266     ana.tx->pushbackToBranch<float>("sim_vz", trk.simvtx_z()[vtxidx]);
0267 
0268     // The trkNtupIdx is the idx in the trackingNtuple
0269     ana.tx->pushbackToBranch<float>("sim_trkNtupIdx", isimtrk);
0270 
0271     // Increase the counter for accepted simtrk
0272     n_accepted_simtrk++;
0273   }
0274 
0275   // Intermediate variables to keep track of matched track candidates for a given sim track
0276   std::vector<int> sim_TC_matched(n_accepted_simtrk);
0277   std::vector<int> sim_TC_matched_mask(n_accepted_simtrk);
0278   std::vector<int> sim_TC_matched_for_duplicate(trk.sim_pt().size());
0279 
0280   // Intermediate variables to keep track of matched sim tracks for a given track candidate
0281   std::vector<std::vector<int>> tc_matched_simIdx;
0282 
0283   // ============ Track candidates =============
0284   auto const& trackCandidates = event->getTrackCandidates();
0285   unsigned int nTrackCandidates = trackCandidates.nTrackCandidates();
0286   for (unsigned int idx = 0; idx < nTrackCandidates; idx++) {
0287     // Compute reco quantities of track candidate based on final object
0288     int type, isFake;
0289     float pt, eta, phi;
0290     std::vector<int> simidx;
0291     std::tie(type, pt, eta, phi, isFake, simidx) = parseTrackCandidate(event, idx);
0292     ana.tx->pushbackToBranch<float>("tc_pt", pt);
0293     ana.tx->pushbackToBranch<float>("tc_eta", eta);
0294     ana.tx->pushbackToBranch<float>("tc_phi", phi);
0295     ana.tx->pushbackToBranch<int>("tc_type", type);
0296     ana.tx->pushbackToBranch<int>("tc_isFake", isFake);
0297     tc_matched_simIdx.push_back(simidx);
0298 
0299     // Loop over matched sim idx and increase counter of TC_matched
0300     for (auto& idx : simidx) {
0301       // NOTE Important to note that the idx of the std::vector<> is same
0302       // as the tracking-ntuple's sim track idx ONLY because event==0 and bunchCrossing==0 condition is applied!!
0303       // Also do not try to access beyond the event and bunchCrossing
0304       if (idx < n_accepted_simtrk) {
0305         sim_TC_matched.at(idx) += 1;
0306         sim_TC_matched_mask.at(idx) |= (1 << type);
0307       }
0308       sim_TC_matched_for_duplicate.at(idx) += 1;
0309     }
0310   }
0311 
0312   // Using the intermedaite variables to compute whether a given track candidate is a duplicate
0313   std::vector<int> tc_isDuplicate(tc_matched_simIdx.size());
0314   // Loop over the track candidates
0315   for (unsigned int i = 0; i < tc_matched_simIdx.size(); ++i) {
0316     bool isDuplicate = false;
0317     // Loop over the sim idx matched to this track candidate
0318     for (unsigned int isim = 0; isim < tc_matched_simIdx[i].size(); ++isim) {
0319       // Using the sim_TC_matched to see whether this track candidate is matched to a sim track that is matched to more than one
0320       int simidx = tc_matched_simIdx[i][isim];
0321       if (sim_TC_matched_for_duplicate[simidx] > 1) {
0322         isDuplicate = true;
0323       }
0324     }
0325     tc_isDuplicate[i] = isDuplicate;
0326   }
0327 
0328   // Now set the last remaining branches
0329   ana.tx->setBranch<std::vector<int>>("sim_TC_matched", sim_TC_matched);
0330   ana.tx->setBranch<std::vector<int>>("sim_TC_matched_mask", sim_TC_matched_mask);
0331   ana.tx->setBranch<std::vector<std::vector<int>>>("tc_matched_simIdx", tc_matched_simIdx);
0332   ana.tx->setBranch<std::vector<int>>("tc_isDuplicate", tc_isDuplicate);
0333 }
0334 
0335 //________________________________________________________________________________________________________________________________
0336 void setOptionalOutputBranches(LSTEvent* event) {
0337 #ifdef CUT_VALUE_DEBUG
0338 
0339   setPixelQuintupletOutputBranches(event);
0340   setQuintupletOutputBranches(event);
0341   setPixelTripletOutputBranches(event);
0342   setOccupancyBranches(event);
0343   setT5DNNBranches(event);
0344 
0345 #endif
0346 }
0347 
0348 //________________________________________________________________________________________________________________________________
0349 void setOccupancyBranches(LSTEvent* event) {
0350   auto modules = event->getModules<ModulesSoA>();
0351   auto miniDoublets = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
0352   auto segments = event->getSegments<SegmentsOccupancySoA>();
0353   auto triplets = event->getTriplets<TripletsOccupancySoA>();
0354   auto quintuplets = event->getQuintuplets<QuintupletsOccupancySoA>();
0355   auto pixelQuintuplets = event->getPixelQuintuplets();
0356   auto pixelTriplets = event->getPixelTriplets();
0357   auto trackCandidates = event->getTrackCandidates();
0358 
0359   std::vector<int> moduleLayer;
0360   std::vector<int> moduleSubdet;
0361   std::vector<int> moduleRing;
0362   std::vector<int> moduleRod;
0363   std::vector<int> moduleModule;
0364   std::vector<float> moduleEta;
0365   std::vector<float> moduleR;
0366   std::vector<bool> moduleIsTilted;
0367   std::vector<int> trackCandidateOccupancy;
0368   std::vector<int> tripletOccupancy;
0369   std::vector<int> segmentOccupancy;
0370   std::vector<int> mdOccupancy;
0371   std::vector<int> quintupletOccupancy;
0372 
0373   for (unsigned int lowerIdx = 0; lowerIdx <= modules.nLowerModules(); lowerIdx++) {
0374     //layer = 0, subdet = 0 => pixel module
0375     moduleLayer.push_back(modules.layers()[lowerIdx]);
0376     moduleSubdet.push_back(modules.subdets()[lowerIdx]);
0377     moduleRing.push_back(modules.rings()[lowerIdx]);
0378     moduleRod.push_back(modules.rods()[lowerIdx]);
0379     moduleEta.push_back(modules.eta()[lowerIdx]);
0380     moduleR.push_back(modules.r()[lowerIdx]);
0381     bool isTilted = (modules.subdets()[lowerIdx] == 5 and modules.sides()[lowerIdx] != 3);
0382     moduleIsTilted.push_back(isTilted);
0383     moduleModule.push_back(modules.modules()[lowerIdx]);
0384     segmentOccupancy.push_back(segments.totOccupancySegments()[lowerIdx]);
0385     mdOccupancy.push_back(miniDoublets.totOccupancyMDs()[lowerIdx]);
0386 
0387     if (lowerIdx < modules.nLowerModules()) {
0388       quintupletOccupancy.push_back(quintuplets.totOccupancyQuintuplets()[lowerIdx]);
0389       tripletOccupancy.push_back(triplets.totOccupancyTriplets()[lowerIdx]);
0390     }
0391   }
0392 
0393   ana.tx->setBranch<std::vector<int>>("module_layers", moduleLayer);
0394   ana.tx->setBranch<std::vector<int>>("module_subdets", moduleSubdet);
0395   ana.tx->setBranch<std::vector<int>>("module_rings", moduleRing);
0396   ana.tx->setBranch<std::vector<int>>("module_rods", moduleRod);
0397   ana.tx->setBranch<std::vector<int>>("module_modules", moduleModule);
0398   ana.tx->setBranch<std::vector<bool>>("module_isTilted", moduleIsTilted);
0399   ana.tx->setBranch<std::vector<float>>("module_eta", moduleEta);
0400   ana.tx->setBranch<std::vector<float>>("module_r", moduleR);
0401   ana.tx->setBranch<std::vector<int>>("md_occupancies", mdOccupancy);
0402   ana.tx->setBranch<std::vector<int>>("sg_occupancies", segmentOccupancy);
0403   ana.tx->setBranch<std::vector<int>>("t3_occupancies", tripletOccupancy);
0404   ana.tx->setBranch<int>("tc_occupancies", trackCandidates.nTrackCandidates());
0405   ana.tx->setBranch<int>("pT3_occupancies", pixelTriplets.totOccupancyPixelTriplets());
0406   ana.tx->setBranch<std::vector<int>>("t5_occupancies", quintupletOccupancy);
0407   ana.tx->setBranch<int>("pT5_occupancies", pixelQuintuplets.totOccupancyPixelQuintuplets());
0408 }
0409 
0410 //________________________________________________________________________________________________________________________________
0411 void setPixelQuintupletOutputBranches(LSTEvent* event) {
0412   // ============ pT5 =============
0413   auto const pixelQuintuplets = event->getPixelQuintuplets();
0414   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
0415   auto const segmentsPixel = event->getSegments<SegmentsPixelSoA>();
0416   auto modules = event->getModules<ModulesSoA>();
0417   int n_accepted_simtrk = ana.tx->getBranch<std::vector<int>>("sim_TC_matched").size();
0418 
0419   unsigned int nPixelQuintuplets = pixelQuintuplets.nPixelQuintuplets();
0420   std::vector<int> sim_pT5_matched(n_accepted_simtrk);
0421   std::vector<std::vector<int>> pT5_matched_simIdx;
0422 
0423   for (unsigned int pT5 = 0; pT5 < nPixelQuintuplets; pT5++) {
0424     unsigned int T5Index = getT5FrompT5(event, pT5);
0425     unsigned int pLSIndex = getPixelLSFrompT5(event, pT5);
0426     float pt = (__H2F(quintuplets.innerRadius()[T5Index]) * k2Rinv1GeVf * 2 + segmentsPixel.ptIn()[pLSIndex]) / 2;
0427     float eta = segmentsPixel.eta()[pLSIndex];
0428     float phi = segmentsPixel.phi()[pLSIndex];
0429 
0430     std::vector<unsigned int> hit_idx = getHitIdxsFrompT5(event, pT5);
0431     std::vector<unsigned int> module_idx = getModuleIdxsFrompT5(event, pT5);
0432     std::vector<unsigned int> hit_type = getHitTypesFrompT5(event, pT5);
0433 
0434     int layer_binary = 1;
0435     int moduleType_binary = 0;
0436     for (size_t i = 0; i < module_idx.size(); i += 2) {
0437       layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4)));
0438       moduleType_binary |= (modules.moduleType()[module_idx[i]] << i);
0439     }
0440     std::vector<int> simidx = matchedSimTrkIdxs(hit_idx, hit_type);
0441     ana.tx->pushbackToBranch<int>("pT5_isFake", static_cast<int>(simidx.size() == 0));
0442     ana.tx->pushbackToBranch<float>("pT5_pt", pt);
0443     ana.tx->pushbackToBranch<float>("pT5_eta", eta);
0444     ana.tx->pushbackToBranch<float>("pT5_phi", phi);
0445     ana.tx->pushbackToBranch<int>("pT5_layer_binary", layer_binary);
0446     ana.tx->pushbackToBranch<int>("pT5_moduleType_binary", moduleType_binary);
0447     ana.tx->pushbackToBranch<float>("pT5_rzChiSquared", pixelQuintuplets.rzChiSquared()[pT5]);
0448 
0449     pT5_matched_simIdx.push_back(simidx);
0450 
0451     // Loop over matched sim idx and increase counter of pT5_matched
0452     for (auto& idx : simidx) {
0453       // NOTE Important to note that the idx of the std::vector<> is same
0454       // as the tracking-ntuple's sim track idx ONLY because event==0 and bunchCrossing==0 condition is applied!!
0455       // Also do not try to access beyond the event and bunchCrossing
0456       if (idx < n_accepted_simtrk) {
0457         sim_pT5_matched.at(idx) += 1;
0458       }
0459     }
0460   }
0461 
0462   // Using the intermedaite variables to compute whether a given track candidate is a duplicate
0463   std::vector<int> pT5_isDuplicate(pT5_matched_simIdx.size());
0464   // Loop over the track candidates
0465   for (unsigned int i = 0; i < pT5_matched_simIdx.size(); ++i) {
0466     bool isDuplicate = false;
0467     // Loop over the sim idx matched to this track candidate
0468     for (unsigned int isim = 0; isim < pT5_matched_simIdx[i].size(); ++isim) {
0469       // Using the sim_pT5_matched to see whether this track candidate is matched to a sim track that is matched to more than one
0470       int simidx = pT5_matched_simIdx[i][isim];
0471       if (simidx < n_accepted_simtrk) {
0472         if (sim_pT5_matched[simidx] > 1) {
0473           isDuplicate = true;
0474         }
0475       }
0476     }
0477     pT5_isDuplicate[i] = isDuplicate;
0478   }
0479 
0480   // Now set the last remaining branches
0481   ana.tx->setBranch<std::vector<int>>("sim_pT5_matched", sim_pT5_matched);
0482   ana.tx->setBranch<std::vector<std::vector<int>>>("pT5_matched_simIdx", pT5_matched_simIdx);
0483   ana.tx->setBranch<std::vector<int>>("pT5_isDuplicate", pT5_isDuplicate);
0484 }
0485 
0486 //________________________________________________________________________________________________________________________________
0487 void setQuintupletOutputBranches(LSTEvent* event) {
0488   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
0489   auto const quintupletsOccupancy = event->getQuintuplets<QuintupletsOccupancySoA>();
0490   auto ranges = event->getRanges();
0491   auto modules = event->getModules<ModulesSoA>();
0492   int n_accepted_simtrk = ana.tx->getBranch<std::vector<int>>("sim_TC_matched").size();
0493 
0494   std::vector<int> sim_t5_matched(n_accepted_simtrk);
0495   std::vector<std::vector<int>> t5_matched_simIdx;
0496 
0497   for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < modules.nLowerModules(); ++lowerModuleIdx) {
0498     int nQuintuplets = quintupletsOccupancy.nQuintuplets()[lowerModuleIdx];
0499     for (unsigned int idx = 0; idx < nQuintuplets; idx++) {
0500       unsigned int quintupletIndex = ranges.quintupletModuleIndices()[lowerModuleIdx] + idx;
0501       float pt = __H2F(quintuplets.innerRadius()[quintupletIndex]) * k2Rinv1GeVf * 2;
0502       float eta = __H2F(quintuplets.eta()[quintupletIndex]);
0503       float phi = __H2F(quintuplets.phi()[quintupletIndex]);
0504 
0505       std::vector<unsigned int> hit_idx = getHitIdxsFromT5(event, quintupletIndex);
0506       std::vector<unsigned int> hit_type = getHitTypesFromT5(event, quintupletIndex);
0507       std::vector<unsigned int> module_idx = getModuleIdxsFromT5(event, quintupletIndex);
0508 
0509       int layer_binary = 0;
0510       int moduleType_binary = 0;
0511       for (size_t i = 0; i < module_idx.size(); i += 2) {
0512         layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4)));
0513         moduleType_binary |= (modules.moduleType()[module_idx[i]] << i);
0514       }
0515 
0516       float percent_matched;
0517       std::vector<int> simidx = matchedSimTrkIdxs(hit_idx, hit_type, false, &percent_matched);
0518 
0519       ana.tx->pushbackToBranch<int>("t5_isFake", static_cast<int>(simidx.size() == 0));
0520       ana.tx->pushbackToBranch<float>("t5_pt", pt);
0521       ana.tx->pushbackToBranch<float>("t5_pMatched", percent_matched);
0522       ana.tx->pushbackToBranch<float>("t5_eta", eta);
0523       ana.tx->pushbackToBranch<float>("t5_phi", phi);
0524       ana.tx->pushbackToBranch<float>("t5_innerRadius", __H2F(quintuplets.innerRadius()[quintupletIndex]));
0525       ana.tx->pushbackToBranch<float>("t5_bridgeRadius", __H2F(quintuplets.bridgeRadius()[quintupletIndex]));
0526       ana.tx->pushbackToBranch<float>("t5_outerRadius", __H2F(quintuplets.outerRadius()[quintupletIndex]));
0527       ana.tx->pushbackToBranch<float>("t5_chiSquared", quintuplets.chiSquared()[quintupletIndex]);
0528       ana.tx->pushbackToBranch<float>("t5_rzChiSquared", quintuplets.rzChiSquared()[quintupletIndex]);
0529       ana.tx->pushbackToBranch<float>("t5_nonAnchorChiSquared", quintuplets.nonAnchorChiSquared()[quintupletIndex]);
0530       ana.tx->pushbackToBranch<float>("t5_dBeta1", quintuplets.dBeta1()[quintupletIndex]);
0531       ana.tx->pushbackToBranch<float>("t5_dBeta2", quintuplets.dBeta2()[quintupletIndex]);
0532       ana.tx->pushbackToBranch<int>("t5_layer_binary", layer_binary);
0533       ana.tx->pushbackToBranch<int>("t5_moduleType_binary", moduleType_binary);
0534 
0535       t5_matched_simIdx.push_back(simidx);
0536 
0537       for (auto& simtrk : simidx) {
0538         if (simtrk < n_accepted_simtrk) {
0539           sim_t5_matched.at(simtrk) += 1;
0540         }
0541       }
0542 
0543       // Avoid fakes when calculating the vertex distance, set default to 0.0.
0544       if (simidx.size() == 0) {
0545         ana.tx->pushbackToBranch<float>("t5_sim_vxy", 0.0);
0546         ana.tx->pushbackToBranch<float>("t5_sim_vz", 0.0);
0547         continue;
0548       }
0549 
0550       int vtxidx = trk.sim_parentVtxIdx()[simidx[0]];
0551       float vtx_x = trk.simvtx_x()[vtxidx];
0552       float vtx_y = trk.simvtx_y()[vtxidx];
0553       float vtx_z = trk.simvtx_z()[vtxidx];
0554 
0555       ana.tx->pushbackToBranch<float>("t5_sim_vxy", sqrt(vtx_x * vtx_x + vtx_y * vtx_y));
0556       ana.tx->pushbackToBranch<float>("t5_sim_vz", vtx_z);
0557     }
0558   }
0559 
0560   std::vector<int> t5_isDuplicate(t5_matched_simIdx.size());
0561   for (unsigned int i = 0; i < t5_matched_simIdx.size(); i++) {
0562     bool isDuplicate = false;
0563     for (unsigned int isim = 0; isim < t5_matched_simIdx[i].size(); isim++) {
0564       int simidx = t5_matched_simIdx[i][isim];
0565       if (simidx < n_accepted_simtrk) {
0566         if (sim_t5_matched[simidx] > 1) {
0567           isDuplicate = true;
0568         }
0569       }
0570     }
0571     t5_isDuplicate[i] = isDuplicate;
0572   }
0573   ana.tx->setBranch<std::vector<int>>("sim_T5_matched", sim_t5_matched);
0574   ana.tx->setBranch<std::vector<std::vector<int>>>("t5_matched_simIdx", t5_matched_simIdx);
0575   ana.tx->setBranch<std::vector<int>>("t5_isDuplicate", t5_isDuplicate);
0576 }
0577 
0578 //________________________________________________________________________________________________________________________________
0579 void setPixelTripletOutputBranches(LSTEvent* event) {
0580   auto const pixelTriplets = event->getPixelTriplets();
0581   auto modules = event->getModules<ModulesSoA>();
0582   SegmentsPixelConst segmentsPixel = event->getSegments<SegmentsPixelSoA>();
0583   int n_accepted_simtrk = ana.tx->getBranch<std::vector<int>>("sim_TC_matched").size();
0584 
0585   unsigned int nPixelTriplets = pixelTriplets.nPixelTriplets();
0586   std::vector<int> sim_pT3_matched(n_accepted_simtrk);
0587   std::vector<std::vector<int>> pT3_matched_simIdx;
0588 
0589   for (unsigned int pT3 = 0; pT3 < nPixelTriplets; pT3++) {
0590     unsigned int T3Index = getT3FrompT3(event, pT3);
0591     unsigned int pLSIndex = getPixelLSFrompT3(event, pT3);
0592     const float pt = segmentsPixel.ptIn()[pLSIndex];
0593 
0594     float eta = segmentsPixel.eta()[pLSIndex];
0595     float phi = segmentsPixel.phi()[pLSIndex];
0596     std::vector<unsigned int> hit_idx = getHitIdxsFrompT3(event, pT3);
0597     std::vector<unsigned int> hit_type = getHitTypesFrompT3(event, pT3);
0598 
0599     std::vector<int> simidx = matchedSimTrkIdxs(hit_idx, hit_type);
0600     std::vector<unsigned int> module_idx = getModuleIdxsFrompT3(event, pT3);
0601     int layer_binary = 1;
0602     int moduleType_binary = 0;
0603     for (size_t i = 0; i < module_idx.size(); i += 2) {
0604       layer_binary |= (1 << (modules.layers()[module_idx[i]] + 6 * (modules.subdets()[module_idx[i]] == 4)));
0605       moduleType_binary |= (modules.moduleType()[module_idx[i]] << i);
0606     }
0607     ana.tx->pushbackToBranch<int>("pT3_isFake", static_cast<int>(simidx.size() == 0));
0608     ana.tx->pushbackToBranch<float>("pT3_pt", pt);
0609     ana.tx->pushbackToBranch<float>("pT3_eta", eta);
0610     ana.tx->pushbackToBranch<float>("pT3_phi", phi);
0611     ana.tx->pushbackToBranch<int>("pT3_layer_binary", layer_binary);
0612     ana.tx->pushbackToBranch<int>("pT3_moduleType_binary", moduleType_binary);
0613 
0614     pT3_matched_simIdx.push_back(simidx);
0615 
0616     for (auto& idx : simidx) {
0617       if (idx < n_accepted_simtrk) {
0618         sim_pT3_matched.at(idx) += 1;
0619       }
0620     }
0621   }
0622 
0623   std::vector<int> pT3_isDuplicate(pT3_matched_simIdx.size());
0624   for (unsigned int i = 0; i < pT3_matched_simIdx.size(); i++) {
0625     bool isDuplicate = true;
0626     for (unsigned int isim = 0; isim < pT3_matched_simIdx[i].size(); isim++) {
0627       int simidx = pT3_matched_simIdx[i][isim];
0628       if (simidx < n_accepted_simtrk) {
0629         if (sim_pT3_matched[simidx] > 1) {
0630           isDuplicate = true;
0631         }
0632       }
0633     }
0634     pT3_isDuplicate[i] = isDuplicate;
0635   }
0636   ana.tx->setBranch<std::vector<int>>("sim_pT3_matched", sim_pT3_matched);
0637   ana.tx->setBranch<std::vector<std::vector<int>>>("pT3_matched_simIdx", pT3_matched_simIdx);
0638   ana.tx->setBranch<std::vector<int>>("pT3_isDuplicate", pT3_isDuplicate);
0639 }
0640 
0641 //________________________________________________________________________________________________________________________________
0642 void fillT5DNNBranches(LSTEvent* event, unsigned int iT3) {
0643   auto hits = event->getHits<HitsSoA>();
0644   auto modules = event->getModules<ModulesSoA>();
0645 
0646   std::vector<unsigned int> hitIdx = getHitsFromT3(event, iT3);
0647   std::vector<lst_math::Hit> hitObjects(hitIdx.size());
0648 
0649   for (int i = 0; i < hitIdx.size(); ++i) {
0650     unsigned int hit = hitIdx[i];
0651     float x = hits.xs()[hit];
0652     float y = hits.ys()[hit];
0653     float z = hits.zs()[hit];
0654     hitObjects[i] = lst_math::Hit(x, y, z);
0655 
0656     std::string idx = std::to_string(i);
0657     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_r", sqrt(x * x + y * y));
0658     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_x", x);
0659     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_y", y);
0660     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_z", z);
0661     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_eta", hitObjects[i].eta());
0662     ana.tx->pushbackToBranch<float>("t5_t3_" + idx + "_phi", hitObjects[i].phi());
0663 
0664     int subdet = trk.ph2_subdet()[hits.idxs()[hit]];
0665     int is_endcap = subdet == 4;
0666     int layer = trk.ph2_layer()[hits.idxs()[hit]] + 6 * is_endcap;
0667     int detId = trk.ph2_detId()[hits.idxs()[hit]];
0668     unsigned int module = hits.moduleIndices()[hit];
0669 
0670     ana.tx->pushbackToBranch<int>("t5_t3_" + idx + "_detId", detId);
0671     ana.tx->pushbackToBranch<int>("t5_t3_" + idx + "_layer", layer);
0672     ana.tx->pushbackToBranch<int>("t5_t3_" + idx + "_moduleType", modules.moduleType()[module]);
0673   }
0674 
0675   float g, f;
0676   auto const& devHost = cms::alpakatools::host();
0677   float radius = computeRadiusFromThreeAnchorHits(devHost,
0678                                                   hitObjects[0].x(),
0679                                                   hitObjects[0].y(),
0680                                                   hitObjects[1].x(),
0681                                                   hitObjects[1].y(),
0682                                                   hitObjects[2].x(),
0683                                                   hitObjects[2].y(),
0684                                                   g,
0685                                                   f);
0686   ana.tx->pushbackToBranch<float>("t5_t3_pt", k2Rinv1GeVf * 2 * radius);
0687 
0688   // Angles
0689   ana.tx->pushbackToBranch<float>("t5_t3_eta", hitObjects[2].eta());
0690   ana.tx->pushbackToBranch<float>("t5_t3_phi", hitObjects[0].phi());
0691 }
0692 
0693 //________________________________________________________________________________________________________________________________
0694 void setT5DNNBranches(LSTEvent* event) {
0695   auto triplets = event->getTriplets<TripletsOccupancySoA>();
0696   auto modules = event->getModules<ModulesSoA>();
0697   auto ranges = event->getRanges();
0698   auto const quintuplets = event->getQuintuplets<QuintupletsOccupancySoA>();
0699   auto trackCandidates = event->getTrackCandidates();
0700 
0701   std::unordered_set<unsigned int> allT3s;
0702   std::unordered_map<unsigned int, unsigned int> t3_index_map;
0703 
0704   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
0705     for (unsigned int jdx = 0; jdx < triplets.nTriplets()[idx]; ++jdx) {
0706       unsigned int t3Idx = ranges.tripletModuleIndices()[idx] + jdx;
0707       if (allT3s.insert(t3Idx).second) {
0708         t3_index_map[t3Idx] = allT3s.size() - 1;
0709         fillT5DNNBranches(event, t3Idx);
0710       }
0711     }
0712   }
0713 
0714   std::unordered_map<unsigned int, unsigned int> t5_tc_index_map;
0715   std::unordered_set<unsigned int> t5s_used_in_tc;
0716 
0717   for (unsigned int idx = 0; idx < trackCandidates.nTrackCandidates(); idx++) {
0718     if (trackCandidates.trackCandidateType()[idx] == LSTObjType::T5) {
0719       unsigned int objIdx = trackCandidates.directObjectIndices()[idx];
0720       t5s_used_in_tc.insert(objIdx);
0721       t5_tc_index_map[objIdx] = idx;
0722     }
0723   }
0724 
0725   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
0726     for (unsigned int jdx = 0; jdx < quintuplets.nQuintuplets()[idx]; ++jdx) {
0727       unsigned int t5Idx = ranges.quintupletModuleIndices()[idx] + jdx;
0728       std::vector<unsigned int> t3sIdx = getT3sFromT5(event, t5Idx);
0729 
0730       ana.tx->pushbackToBranch<int>("t5_t3_idx0", t3_index_map[t3sIdx[0]]);
0731       ana.tx->pushbackToBranch<int>("t5_t3_idx1", t3_index_map[t3sIdx[1]]);
0732 
0733       if (t5s_used_in_tc.find(t5Idx) != t5s_used_in_tc.end()) {
0734         ana.tx->pushbackToBranch<int>("t5_partOfTC", 1);
0735         ana.tx->pushbackToBranch<int>("t5_tc_idx", t5_tc_index_map[t5Idx]);
0736       } else {
0737         ana.tx->pushbackToBranch<int>("t5_partOfTC", 0);
0738         ana.tx->pushbackToBranch<int>("t5_tc_idx", -999);
0739       }
0740     }
0741   }
0742 }
0743 
0744 //________________________________________________________________________________________________________________________________
0745 void setGnnNtupleBranches(LSTEvent* event) {
0746   // Get relevant information
0747   SegmentsOccupancyConst segmentsOccupancy = event->getSegments<SegmentsOccupancySoA>();
0748   MiniDoubletsOccupancyConst miniDoublets = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
0749   auto hitsEvt = event->getHits<HitsSoA>();
0750   auto modules = event->getModules<ModulesSoA>();
0751   auto ranges = event->getRanges();
0752   auto const& trackCandidates = event->getTrackCandidates();
0753 
0754   std::set<unsigned int> mds_used_in_sg;
0755   std::map<unsigned int, unsigned int> md_index_map;
0756   std::map<unsigned int, unsigned int> sg_index_map;
0757 
0758   // Loop over modules (lower ones where the MDs are saved)
0759   unsigned int nTotalMD = 0;
0760   unsigned int nTotalLS = 0;
0761   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
0762     nTotalMD += miniDoublets.nMDs()[idx];
0763     nTotalLS += segmentsOccupancy.nSegments()[idx];
0764   }
0765 
0766   std::set<unsigned int> lss_used_in_true_tc;
0767   unsigned int nTrackCandidates = trackCandidates.nTrackCandidates();
0768   for (unsigned int idx = 0; idx < nTrackCandidates; idx++) {
0769     // Only consider true track candidates
0770     std::vector<unsigned int> hitidxs;
0771     std::vector<unsigned int> hittypes;
0772     std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromTC(event, idx);
0773     std::vector<int> simidxs = matchedSimTrkIdxs(hitidxs, hittypes);
0774     if (simidxs.size() == 0)
0775       continue;
0776 
0777     std::vector<unsigned int> LSs = getLSsFromTC(event, idx);
0778     for (auto& LS : LSs) {
0779       if (lss_used_in_true_tc.find(LS) == lss_used_in_true_tc.end()) {
0780         lss_used_in_true_tc.insert(LS);
0781       }
0782     }
0783   }
0784 
0785   std::cout << " lss_used_in_true_tc.size(): " << lss_used_in_true_tc.size() << std::endl;
0786 
0787   // std::cout <<  " nTotalMD: " << nTotalMD <<  std::endl;
0788   // std::cout <<  " nTotalLS: " << nTotalLS <<  std::endl;
0789 
0790   // Loop over modules (lower ones where the MDs are saved)
0791   for (unsigned int idx = 0; idx < modules.nLowerModules(); ++idx) {
0792     // // Loop over minidoublets
0793     // for (unsigned int jdx = 0; jdx < miniDoublets->nMDs[idx]; jdx++)
0794     // {
0795     //     // Get the actual index to the mini-doublet using ranges
0796     //     unsigned int mdIdx = ranges->miniDoubletModuleIndices[idx] + jdx;
0797 
0798     //     setGnnNtupleMiniDoublet(event, mdIdx);
0799     // }
0800 
0801     // Loop over segments
0802     for (unsigned int jdx = 0; jdx < segmentsOccupancy.nSegments()[idx]; jdx++) {
0803       // Get the actual index to the segments using ranges
0804       unsigned int sgIdx = ranges.segmentModuleIndices()[idx] + jdx;
0805 
0806       // Get the hit indices
0807       std::vector<unsigned int> MDs = getMDsFromLS(event, sgIdx);
0808 
0809       if (mds_used_in_sg.find(MDs[0]) == mds_used_in_sg.end()) {
0810         mds_used_in_sg.insert(MDs[0]);
0811         md_index_map[MDs[0]] = mds_used_in_sg.size() - 1;
0812         setGnnNtupleMiniDoublet(event, MDs[0]);
0813       }
0814 
0815       if (mds_used_in_sg.find(MDs[1]) == mds_used_in_sg.end()) {
0816         mds_used_in_sg.insert(MDs[1]);
0817         md_index_map[MDs[1]] = mds_used_in_sg.size() - 1;
0818         setGnnNtupleMiniDoublet(event, MDs[1]);
0819       }
0820 
0821       ana.tx->pushbackToBranch<int>("LS_MD_idx0", md_index_map[MDs[0]]);
0822       ana.tx->pushbackToBranch<int>("LS_MD_idx1", md_index_map[MDs[1]]);
0823 
0824       std::vector<unsigned int> hits = getHitsFromLS(event, sgIdx);
0825 
0826       // Computing line segment pt estimate (assuming beam spot is at zero)
0827       lst_math::Hit hitA(0, 0, 0);
0828       lst_math::Hit hitB(hitsEvt.xs()[hits[0]], hitsEvt.ys()[hits[0]], hitsEvt.zs()[hits[0]]);
0829       lst_math::Hit hitC(hitsEvt.xs()[hits[2]], hitsEvt.ys()[hits[2]], hitsEvt.zs()[hits[2]]);
0830       lst_math::Hit center = lst_math::getCenterFromThreePoints(hitA, hitB, hitC);
0831       float pt = lst_math::ptEstimateFromRadius(center.rt());
0832       float eta = hitC.eta();
0833       float phi = hitB.phi();
0834 
0835       ana.tx->pushbackToBranch<float>("LS_pt", pt);
0836       ana.tx->pushbackToBranch<float>("LS_eta", eta);
0837       ana.tx->pushbackToBranch<float>("LS_phi", phi);
0838       // ana.tx->pushbackToBranch<int>("LS_layer0", layer0);
0839       // ana.tx->pushbackToBranch<int>("LS_layer1", layer1);
0840 
0841       std::vector<unsigned int> hitidxs;
0842       std::vector<unsigned int> hittypes;
0843       std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromLS(event, sgIdx);
0844       std::vector<int> simidxs = matchedSimTrkIdxs(hitidxs, hittypes);
0845 
0846       ana.tx->pushbackToBranch<int>("LS_isFake", simidxs.size() == 0);
0847       ana.tx->pushbackToBranch<float>("LS_sim_pt", simidxs.size() > 0 ? trk.sim_pt()[simidxs[0]] : -999);
0848       ana.tx->pushbackToBranch<float>("LS_sim_eta", simidxs.size() > 0 ? trk.sim_eta()[simidxs[0]] : -999);
0849       ana.tx->pushbackToBranch<float>("LS_sim_phi", simidxs.size() > 0 ? trk.sim_phi()[simidxs[0]] : -999);
0850       ana.tx->pushbackToBranch<float>("LS_sim_pca_dxy", simidxs.size() > 0 ? trk.sim_pca_dxy()[simidxs[0]] : -999);
0851       ana.tx->pushbackToBranch<float>("LS_sim_pca_dz", simidxs.size() > 0 ? trk.sim_pca_dz()[simidxs[0]] : -999);
0852       ana.tx->pushbackToBranch<int>("LS_sim_q", simidxs.size() > 0 ? trk.sim_q()[simidxs[0]] : -999);
0853       ana.tx->pushbackToBranch<int>("LS_sim_event", simidxs.size() > 0 ? trk.sim_event()[simidxs[0]] : -999);
0854       ana.tx->pushbackToBranch<int>("LS_sim_bx", simidxs.size() > 0 ? trk.sim_bunchCrossing()[simidxs[0]] : -999);
0855       ana.tx->pushbackToBranch<int>("LS_sim_pdgId", simidxs.size() > 0 ? trk.sim_pdgId()[simidxs[0]] : -999);
0856       ana.tx->pushbackToBranch<float>("LS_sim_vx",
0857                                       simidxs.size() > 0 ? trk.simvtx_x()[trk.sim_parentVtxIdx()[simidxs[0]]] : -999);
0858       ana.tx->pushbackToBranch<float>("LS_sim_vy",
0859                                       simidxs.size() > 0 ? trk.simvtx_y()[trk.sim_parentVtxIdx()[simidxs[0]]] : -999);
0860       ana.tx->pushbackToBranch<float>("LS_sim_vz",
0861                                       simidxs.size() > 0 ? trk.simvtx_z()[trk.sim_parentVtxIdx()[simidxs[0]]] : -999);
0862       ana.tx->pushbackToBranch<int>("LS_isInTrueTC", lss_used_in_true_tc.find(sgIdx) != lss_used_in_true_tc.end());
0863 
0864       sg_index_map[sgIdx] = ana.tx->getBranch<std::vector<int>>("LS_isFake").size() - 1;
0865 
0866       // // T5 eta and phi are computed using outer and innermost hits
0867       // lst_math::Hit hitA(trk.ph2_x()[anchitidx], trk.ph2_y()[anchitidx], trk.ph2_z()[anchitidx]);
0868       // const float phi = hitA.phi();
0869       // const float eta = hitA.eta();
0870     }
0871   }
0872 
0873   for (unsigned int idx = 0; idx < nTrackCandidates; idx++) {
0874     std::vector<unsigned int> LSs = getLSsFromTC(event, idx);
0875     std::vector<int> lsIdx;
0876     for (auto& LS : LSs) {
0877       lsIdx.push_back(sg_index_map[LS]);
0878     }
0879     ana.tx->pushbackToBranch<std::vector<int>>("tc_lsIdx", lsIdx);
0880   }
0881 
0882   std::cout << " mds_used_in_sg.size(): " << mds_used_in_sg.size() << std::endl;
0883 }
0884 
0885 //________________________________________________________________________________________________________________________________
0886 void setGnnNtupleMiniDoublet(LSTEvent* event, unsigned int MD) {
0887   // Get relevant information
0888   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
0889   auto hitsEvt = event->getHits<HitsSoA>();
0890 
0891   // Get the hit indices
0892   unsigned int hit0 = miniDoublets.anchorHitIndices()[MD];
0893   unsigned int hit1 = miniDoublets.outerHitIndices()[MD];
0894 
0895   // Get the hit infos
0896   const float hit0_x = hitsEvt.xs()[hit0];
0897   const float hit0_y = hitsEvt.ys()[hit0];
0898   const float hit0_z = hitsEvt.zs()[hit0];
0899   const float hit0_r = sqrt(hit0_x * hit0_x + hit0_y * hit0_y);
0900   const float hit1_x = hitsEvt.xs()[hit1];
0901   const float hit1_y = hitsEvt.ys()[hit1];
0902   const float hit1_z = hitsEvt.zs()[hit1];
0903   const float hit1_r = sqrt(hit1_x * hit1_x + hit1_y * hit1_y);
0904 
0905   // Do sim matching
0906   std::vector<unsigned int> hit_idx = {hitsEvt.idxs()[hit0], hitsEvt.idxs()[hit1]};
0907   std::vector<unsigned int> hit_type = {4, 4};
0908   std::vector<int> simidxs = matchedSimTrkIdxs(hit_idx, hit_type);
0909 
0910   bool isFake = simidxs.size() == 0;
0911   int tp_type = getDenomSimTrkType(simidxs);
0912 
0913   // Obtain where the actual hit is located in terms of their layer, module, rod, and ring number
0914   unsigned int anchitidx = hitsEvt.idxs()[hit0];
0915   int subdet = trk.ph2_subdet()[hitsEvt.idxs()[anchitidx]];
0916   int is_endcap = subdet == 4;
0917   int layer =
0918       trk.ph2_layer()[anchitidx] +
0919       6 * (is_endcap);  // 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)
0920   int detId = trk.ph2_detId()[anchitidx];
0921 
0922   // Obtaining dPhiChange
0923   float dphichange = miniDoublets.dphichanges()[MD];
0924 
0925   // Computing pt
0926   float pt = hit0_r * k2Rinv1GeVf / sin(dphichange);
0927 
0928   // T5 eta and phi are computed using outer and innermost hits
0929   lst_math::Hit hitA(trk.ph2_x()[anchitidx], trk.ph2_y()[anchitidx], trk.ph2_z()[anchitidx]);
0930   const float phi = hitA.phi();
0931   const float eta = hitA.eta();
0932 
0933   // Mini Doublets
0934   ana.tx->pushbackToBranch<float>("MD_pt", pt);
0935   ana.tx->pushbackToBranch<float>("MD_eta", eta);
0936   ana.tx->pushbackToBranch<float>("MD_phi", phi);
0937   ana.tx->pushbackToBranch<float>("MD_dphichange", dphichange);
0938   ana.tx->pushbackToBranch<int>("MD_isFake", isFake);
0939   ana.tx->pushbackToBranch<int>("MD_tpType", tp_type);
0940   ana.tx->pushbackToBranch<int>("MD_detId", detId);
0941   ana.tx->pushbackToBranch<int>("MD_layer", layer);
0942   ana.tx->pushbackToBranch<float>("MD_0_r", hit0_r);
0943   ana.tx->pushbackToBranch<float>("MD_0_x", hit0_x);
0944   ana.tx->pushbackToBranch<float>("MD_0_y", hit0_y);
0945   ana.tx->pushbackToBranch<float>("MD_0_z", hit0_z);
0946   ana.tx->pushbackToBranch<float>("MD_1_r", hit1_r);
0947   ana.tx->pushbackToBranch<float>("MD_1_x", hit1_x);
0948   ana.tx->pushbackToBranch<float>("MD_1_y", hit1_y);
0949   ana.tx->pushbackToBranch<float>("MD_1_z", hit1_z);
0950   // ana.tx->pushbackToBranch<int>("MD_sim_idx", simidxs.size() > 0 ? simidxs[0] : -999);
0951 }
0952 
0953 //________________________________________________________________________________________________________________________________
0954 std::tuple<int, float, float, float, int, std::vector<int>> parseTrackCandidate(LSTEvent* event, unsigned int idx) {
0955   // Get the type of the track candidate
0956   auto const& trackCandidates = event->getTrackCandidates();
0957   short type = trackCandidates.trackCandidateType()[idx];
0958 
0959   // Compute pt eta phi and hit indices that will be used to figure out whether the TC matched
0960   float pt, eta, phi;
0961   std::vector<unsigned int> hit_idx, hit_type;
0962   switch (type) {
0963     case LSTObjType::pT5:
0964       std::tie(pt, eta, phi, hit_idx, hit_type) = parsepT5(event, idx);
0965       break;
0966     case LSTObjType::pT3:
0967       std::tie(pt, eta, phi, hit_idx, hit_type) = parsepT3(event, idx);
0968       break;
0969     case LSTObjType::T5:
0970       std::tie(pt, eta, phi, hit_idx, hit_type) = parseT5(event, idx);
0971       break;
0972     case LSTObjType::pLS:
0973       std::tie(pt, eta, phi, hit_idx, hit_type) = parsepLS(event, idx);
0974       break;
0975   }
0976 
0977   // Perform matching
0978   std::vector<int> simidx = matchedSimTrkIdxs(hit_idx, hit_type);
0979   int isFake = simidx.size() == 0;
0980 
0981   return {type, pt, eta, phi, isFake, simidx};
0982 }
0983 
0984 //________________________________________________________________________________________________________________________________
0985 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parsepT5(LSTEvent* event,
0986                                                                                                unsigned int idx) {
0987   // Get relevant information
0988   auto const trackCandidates = event->getTrackCandidates();
0989   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
0990   auto const segmentsPixel = event->getSegments<SegmentsPixelSoA>();
0991 
0992   //
0993   // pictorial representation of a pT5
0994   //
0995   // inner tracker        outer tracker
0996   // -------------  --------------------------
0997   // pLS            01    23    45    67    89   (anchor hit of a minidoublet is always the first of the pair)
0998   // ****           oo -- oo -- oo -- oo -- oo   pT5
0999   //                oo -- oo -- oo               first T3 of the T5
1000   //                            oo -- oo -- oo   second T3 of the T5
1001   unsigned int pT5 = trackCandidates.directObjectIndices()[idx];
1002   unsigned int pLS = getPixelLSFrompT5(event, pT5);
1003   unsigned int T5Index = getT5FrompT5(event, pT5);
1004 
1005   //=================================================================================
1006   // Some history and geometry lesson...
1007   // For a given T3, we compute two angles. (NOTE: This is a bit weird!)
1008   // Historically, T3 were created out of T4, which we used to build a long time ago.
1009   // So for the sake of argument let's discuss T4 first.
1010   // For a T4, we have 4 mini-doublets.
1011   // Therefore we have 4 "anchor hits".
1012   // Therefore we have 4 xyz points.
1013   //
1014   //
1015   //       *
1016   //       |\
1017     //       | \
1018     //       |1 \
1019     //       |   \
1020     //       |  * \
1021     //       |
1022   //       |
1023   //       |
1024   //       |
1025   //       |
1026   //       |  * /
1027   //       |   /
1028   //       |2 /
1029   //       | /
1030   //       |/
1031   //       *
1032   //
1033   //
1034   // Then from these 4 points, one can approximate a some sort of "best" fitted circle trajectory,
1035   // and obtain "tangential" angles from 1st and 4th hits.
1036   // See the carton below.
1037   // The "*" are the 4 physical hit points
1038   // angle 1 and 2 are the "tangential" angle for a "circle" from 4 * points.
1039   // Please note, that a straight line from first two * and the latter two * are NOT the
1040   // angle 1 and angle 2. (they were called "beta" angles)
1041   // But rather, a slightly larger angle.
1042   // Because 4 * points would be on a circle, and a tangential line on the circles
1043   // would deviate from the points on circles.
1044   //
1045   // In the early days of LST, there was an iterative algorithm (devised by Slava) to
1046   // obtain the angle beta1 and 2 _without_ actually performing a 4 point circle fit.
1047   // Hence, the beta1 and beta2 were quickly estimated without too many math operations
1048   // and afterwards (beta1-beta2) was computed to obtain what we call a "delta-beta" values.
1049   //
1050   // For a real track, the deltabeta ~ 0, for fakes, it'd have a flat distribution.
1051   //
1052   // However, after some time we abandonded the T4s, and moved to T3s.
1053   // In T3, however, now we have the following cartoon:
1054   //
1055   //       *
1056   //       |\
1057     //       | \
1058     //       |1 \
1059     //       |   \
1060     //       |  * X   (* here are "two" MDs but really just one)
1061   //       |   /
1062   //       |2 /
1063   //       | /
1064   //       |/
1065   //       *
1066   //
1067   // With the "four" *'s (really just "three") you can still run the iterative beta calculation,
1068   // which is what we still currently do, we still get two beta1 and beta2
1069   // But! high school geometry tells us that 3 points = ONLY 1 possible CIRCLE!
1070   // There is really nothing to "fit" here.
1071   // YET we still compute these in T3, out of legacy method of how we used to treat T4s.
1072   //
1073   // Hence, in the below code, "betaIn_in" and "betaOut_in" if we performed
1074   // a circle fit they would come out by definition identical values.
1075   // But due to our approximate iterative beta calculation method, they come out different values.
1076   // So if we are "cutting on" abs(deltaBeta) = abs(betaIn_in - betaOut_in) < threshold,
1077   // what does that even mean?
1078   //
1079   // Anyhow, as of now, we compute 2 beta's for T3s, and T5 has two T3s.
1080   // And from there we estimate the pt's and we compute pt_T5.
1081 
1082   // pixel pt
1083   const float pt_pLS = segmentsPixel.ptIn()[pLS];
1084   const float eta_pLS = segmentsPixel.eta()[pLS];
1085   const float phi_pLS = segmentsPixel.phi()[pLS];
1086   float pt_T5 = __H2F(quintuplets.innerRadius()[T5Index]) * 2 * k2Rinv1GeVf;
1087   const float pt = (pt_T5 + pt_pLS) / 2;
1088 
1089   // Form the hit idx/type std::vector
1090   std::vector<unsigned int> hit_idx = getHitIdxsFrompT5(event, pT5);
1091   std::vector<unsigned int> hit_type = getHitTypesFrompT5(event, pT5);
1092 
1093   return {pt, eta_pLS, phi_pLS, hit_idx, hit_type};
1094 }
1095 
1096 //________________________________________________________________________________________________________________________________
1097 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parsepT3(LSTEvent* event,
1098                                                                                                unsigned int idx) {
1099   // Get relevant information
1100   auto const trackCandidates = event->getTrackCandidates();
1101   auto const triplets = event->getTriplets<TripletsSoA>();
1102   auto const segmentsPixel = event->getSegments<SegmentsPixelSoA>();
1103 
1104   //
1105   // pictorial representation of a pT3
1106   //
1107   // inner tracker        outer tracker
1108   // -------------  --------------------------
1109   // pLS            01    23    45               (anchor hit of a minidoublet is always the first of the pair)
1110   // ****           oo -- oo -- oo               pT3
1111   unsigned int pT3 = trackCandidates.directObjectIndices()[idx];
1112   unsigned int pLS = getPixelLSFrompT3(event, pT3);
1113   unsigned int T3 = getT3FrompT3(event, pT3);
1114 
1115   // pixel pt
1116   const float pt_pLS = segmentsPixel.ptIn()[pLS];
1117   const float eta_pLS = segmentsPixel.eta()[pLS];
1118   const float phi_pLS = segmentsPixel.phi()[pLS];
1119   float pt_T3 = triplets.radius()[T3] * 2 * k2Rinv1GeVf;
1120 
1121   // average pt
1122   const float pt = (pt_pLS + pt_T3) / 2;
1123 
1124   // Form the hit idx/type std::vector
1125   std::vector<unsigned int> hit_idx = getHitIdxsFrompT3(event, pT3);
1126   std::vector<unsigned int> hit_type = getHitTypesFrompT3(event, pT3);
1127 
1128   return {pt, eta_pLS, phi_pLS, hit_idx, hit_type};
1129 }
1130 
1131 //________________________________________________________________________________________________________________________________
1132 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parseT5(LSTEvent* event,
1133                                                                                               unsigned int idx) {
1134   auto const trackCandidates = event->getTrackCandidates();
1135   auto const quintuplets = event->getQuintuplets<QuintupletsSoA>();
1136   unsigned int T5 = trackCandidates.directObjectIndices()[idx];
1137   std::vector<unsigned int> hits = getHitsFromT5(event, T5);
1138 
1139   //
1140   // pictorial representation of a T5
1141   //
1142   // inner tracker        outer tracker
1143   // -------------  --------------------------
1144   //                01    23    45    67    89   (anchor hit of a minidoublet is always the first of the pair)
1145   //  (none)        oo -- oo -- oo -- oo -- oo   T5
1146   unsigned int Hit_0 = hits[0];
1147   unsigned int Hit_4 = hits[4];
1148   unsigned int Hit_8 = hits[8];
1149 
1150   // T5 radius is average of the inner and outer radius
1151   const float pt = __H2F(quintuplets.innerRadius()[T5]) * k2Rinv1GeVf * 2;
1152 
1153   // T5 eta and phi are computed using outer and innermost hits
1154   lst_math::Hit hitA(trk.ph2_x()[Hit_0], trk.ph2_y()[Hit_0], trk.ph2_z()[Hit_0]);
1155   lst_math::Hit hitB(trk.ph2_x()[Hit_8], trk.ph2_y()[Hit_8], trk.ph2_z()[Hit_8]);
1156   const float phi = hitA.phi();
1157   const float eta = hitB.eta();
1158 
1159   std::vector<unsigned int> hit_idx = getHitIdxsFromT5(event, T5);
1160   std::vector<unsigned int> hit_type = getHitTypesFromT5(event, T5);
1161 
1162   return {pt, eta, phi, hit_idx, hit_type};
1163 }
1164 
1165 //________________________________________________________________________________________________________________________________
1166 std::tuple<float, float, float, std::vector<unsigned int>, std::vector<unsigned int>> parsepLS(LSTEvent* event,
1167                                                                                                unsigned int idx) {
1168   auto const& trackCandidates = event->getTrackCandidates();
1169   SegmentsPixelConst segmentsPixel = event->getSegments<SegmentsPixelSoA>();
1170 
1171   // Getting pLS index
1172   unsigned int pLS = trackCandidates.directObjectIndices()[idx];
1173 
1174   // Getting pt eta and phi
1175   float pt = segmentsPixel.ptIn()[pLS];
1176   float eta = segmentsPixel.eta()[pLS];
1177   float phi = segmentsPixel.phi()[pLS];
1178 
1179   // Getting hit indices and types
1180   std::vector<unsigned int> hit_idx = getPixelHitIdxsFrompLS(event, pLS);
1181   std::vector<unsigned int> hit_type = getPixelHitTypesFrompLS(event, pLS);
1182 
1183   return {pt, eta, phi, hit_idx, hit_type};
1184 }
1185 
1186 //________________________________________________________________________________________________________________________________
1187 void printHitMultiplicities(LSTEvent* event) {
1188   auto modules = event->getModules<ModulesSoA>();
1189   auto hitRanges = event->getHits<HitsRangesSoA>();
1190 
1191   int nHits = 0;
1192   for (unsigned int idx = 0; idx <= modules.nLowerModules();
1193        idx++)  // "<=" because cheating to include pixel track candidate lower module
1194   {
1195     nHits += hitRanges.hitRanges()[2 * idx][1] - hitRanges.hitRanges()[2 * idx][0] + 1;
1196     nHits += hitRanges.hitRanges()[2 * idx + 1][1] - hitRanges.hitRanges()[2 * idx + 1][0] + 1;
1197   }
1198   std::cout << " nHits: " << nHits << std::endl;
1199 }
1200 
1201 //________________________________________________________________________________________________________________________________
1202 void printMiniDoubletMultiplicities(LSTEvent* event) {
1203   MiniDoubletsOccupancyConst miniDoublets = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
1204   auto modules = event->getModules<ModulesSoA>();
1205 
1206   int nMiniDoublets = 0;
1207   int totOccupancyMiniDoublets = 0;
1208   for (unsigned int idx = 0; idx <= modules.nModules();
1209        idx++)  // "<=" because cheating to include pixel track candidate lower module
1210   {
1211     if (modules.isLower()[idx]) {
1212       nMiniDoublets += miniDoublets.nMDs()[idx];
1213       totOccupancyMiniDoublets += miniDoublets.totOccupancyMDs()[idx];
1214     }
1215   }
1216   std::cout << " nMiniDoublets: " << nMiniDoublets << std::endl;
1217   std::cout << " totOccupancyMiniDoublets (including trucated ones): " << totOccupancyMiniDoublets << std::endl;
1218 }
1219 
1220 //________________________________________________________________________________________________________________________________
1221 void printAllObjects(LSTEvent* event) {
1222   printMDs(event);
1223   printLSs(event);
1224   printpLSs(event);
1225   printT3s(event);
1226 }
1227 
1228 //________________________________________________________________________________________________________________________________
1229 void printMDs(LSTEvent* event) {
1230   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1231   MiniDoubletsOccupancyConst miniDoubletsOccupancy = event->getMiniDoublets<MiniDoubletsOccupancySoA>();
1232   auto hitsEvt = event->getHits<HitsSoA>();
1233   auto modules = event->getModules<ModulesSoA>();
1234   auto ranges = event->getRanges();
1235 
1236   // Then obtain the lower module index
1237   for (unsigned int idx = 0; idx <= modules.nLowerModules(); ++idx) {
1238     for (unsigned int iMD = 0; iMD < miniDoubletsOccupancy.nMDs()[idx]; iMD++) {
1239       unsigned int mdIdx = ranges.miniDoubletModuleIndices()[idx] + iMD;
1240       unsigned int LowerHitIndex = miniDoublets.anchorHitIndices()[mdIdx];
1241       unsigned int UpperHitIndex = miniDoublets.outerHitIndices()[mdIdx];
1242       unsigned int hit0 = hitsEvt.idxs()[LowerHitIndex];
1243       unsigned int hit1 = hitsEvt.idxs()[UpperHitIndex];
1244       std::cout << "VALIDATION 'MD': "
1245                 << "MD"
1246                 << " hit0: " << hit0 << " hit1: " << hit1 << std::endl;
1247     }
1248   }
1249 }
1250 
1251 //________________________________________________________________________________________________________________________________
1252 void printLSs(LSTEvent* event) {
1253   SegmentsConst segments = event->getSegments<SegmentsSoA>();
1254   SegmentsOccupancyConst segmentsOccupancy = event->getSegments<SegmentsOccupancySoA>();
1255   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1256   auto hitsEvt = event->getHits<HitsSoA>();
1257   auto modules = event->getModules<ModulesSoA>();
1258   auto ranges = event->getRanges();
1259 
1260   int nSegments = 0;
1261   for (unsigned int i = 0; i < modules.nLowerModules(); ++i) {
1262     unsigned int idx = i;  //modules->lowerModuleIndices[i];
1263     nSegments += segmentsOccupancy.nSegments()[idx];
1264     for (unsigned int jdx = 0; jdx < segmentsOccupancy.nSegments()[idx]; jdx++) {
1265       unsigned int sgIdx = ranges.segmentModuleIndices()[idx] + jdx;
1266       unsigned int InnerMiniDoubletIndex = segments.mdIndices()[sgIdx][0];
1267       unsigned int OuterMiniDoubletIndex = segments.mdIndices()[sgIdx][1];
1268       unsigned int InnerMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[InnerMiniDoubletIndex];
1269       unsigned int InnerMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[InnerMiniDoubletIndex];
1270       unsigned int OuterMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[OuterMiniDoubletIndex];
1271       unsigned int OuterMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[OuterMiniDoubletIndex];
1272       unsigned int hit0 = hitsEvt.idxs()[InnerMiniDoubletLowerHitIndex];
1273       unsigned int hit1 = hitsEvt.idxs()[InnerMiniDoubletUpperHitIndex];
1274       unsigned int hit2 = hitsEvt.idxs()[OuterMiniDoubletLowerHitIndex];
1275       unsigned int hit3 = hitsEvt.idxs()[OuterMiniDoubletUpperHitIndex];
1276       std::cout << "VALIDATION 'LS': "
1277                 << "LS"
1278                 << " hit0: " << hit0 << " hit1: " << hit1 << " hit2: " << hit2 << " hit3: " << hit3 << std::endl;
1279     }
1280   }
1281   std::cout << "VALIDATION nSegments: " << nSegments << std::endl;
1282 }
1283 
1284 //________________________________________________________________________________________________________________________________
1285 void printpLSs(LSTEvent* event) {
1286   SegmentsConst segments = event->getSegments<SegmentsSoA>();
1287   SegmentsOccupancyConst segmentsOccupancy = event->getSegments<SegmentsOccupancySoA>();
1288   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1289   auto hitsEvt = event->getHits<HitsSoA>();
1290   auto modules = event->getModules<ModulesSoA>();
1291   auto ranges = event->getRanges();
1292 
1293   unsigned int i = modules.nLowerModules();
1294   unsigned int idx = i;  //modules->lowerModuleIndices[i];
1295   int npLS = segmentsOccupancy.nSegments()[idx];
1296   for (unsigned int jdx = 0; jdx < segmentsOccupancy.nSegments()[idx]; jdx++) {
1297     unsigned int sgIdx = ranges.segmentModuleIndices()[idx] + jdx;
1298     unsigned int InnerMiniDoubletIndex = segments.mdIndices()[sgIdx][0];
1299     unsigned int OuterMiniDoubletIndex = segments.mdIndices()[sgIdx][1];
1300     unsigned int InnerMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[InnerMiniDoubletIndex];
1301     unsigned int InnerMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[InnerMiniDoubletIndex];
1302     unsigned int OuterMiniDoubletLowerHitIndex = miniDoublets.anchorHitIndices()[OuterMiniDoubletIndex];
1303     unsigned int OuterMiniDoubletUpperHitIndex = miniDoublets.outerHitIndices()[OuterMiniDoubletIndex];
1304     unsigned int hit0 = hitsEvt.idxs()[InnerMiniDoubletLowerHitIndex];
1305     unsigned int hit1 = hitsEvt.idxs()[InnerMiniDoubletUpperHitIndex];
1306     unsigned int hit2 = hitsEvt.idxs()[OuterMiniDoubletLowerHitIndex];
1307     unsigned int hit3 = hitsEvt.idxs()[OuterMiniDoubletUpperHitIndex];
1308     std::cout << "VALIDATION 'pLS': "
1309               << "pLS"
1310               << " hit0: " << hit0 << " hit1: " << hit1 << " hit2: " << hit2 << " hit3: " << hit3 << std::endl;
1311   }
1312   std::cout << "VALIDATION npLS: " << npLS << std::endl;
1313 }
1314 
1315 //________________________________________________________________________________________________________________________________
1316 void printT3s(LSTEvent* event) {
1317   auto const triplets = event->getTriplets<TripletsSoA>();
1318   auto const tripletsOccupancy = event->getTriplets<TripletsOccupancySoA>();
1319   SegmentsConst segments = event->getSegments<SegmentsSoA>();
1320   MiniDoubletsConst miniDoublets = event->getMiniDoublets<MiniDoubletsSoA>();
1321   auto hitsEvt = event->getHits<HitsSoA>();
1322   auto modules = event->getModules<ModulesSoA>();
1323   int nTriplets = 0;
1324   for (unsigned int i = 0; i < modules.nLowerModules(); ++i) {
1325     // unsigned int idx = modules->lowerModuleIndices[i];
1326     nTriplets += tripletsOccupancy.nTriplets()[i];
1327     unsigned int idx = i;
1328     for (unsigned int jdx = 0; jdx < tripletsOccupancy.nTriplets()[idx]; jdx++) {
1329       unsigned int tpIdx = idx * 5000 + jdx;
1330       unsigned int InnerSegmentIndex = triplets.segmentIndices()[tpIdx][0];
1331       unsigned int OuterSegmentIndex = triplets.segmentIndices()[tpIdx][1];
1332       unsigned int InnerSegmentInnerMiniDoubletIndex = segments.mdIndices()[InnerSegmentIndex][0];
1333       unsigned int InnerSegmentOuterMiniDoubletIndex = segments.mdIndices()[InnerSegmentIndex][1];
1334       unsigned int OuterSegmentOuterMiniDoubletIndex = segments.mdIndices()[OuterSegmentIndex][1];
1335 
1336       unsigned int hit_idx0 = miniDoublets.anchorHitIndices()[InnerSegmentInnerMiniDoubletIndex];
1337       unsigned int hit_idx1 = miniDoublets.outerHitIndices()[InnerSegmentInnerMiniDoubletIndex];
1338       unsigned int hit_idx2 = miniDoublets.anchorHitIndices()[InnerSegmentOuterMiniDoubletIndex];
1339       unsigned int hit_idx3 = miniDoublets.outerHitIndices()[InnerSegmentOuterMiniDoubletIndex];
1340       unsigned int hit_idx4 = miniDoublets.anchorHitIndices()[OuterSegmentOuterMiniDoubletIndex];
1341       unsigned int hit_idx5 = miniDoublets.outerHitIndices()[OuterSegmentOuterMiniDoubletIndex];
1342 
1343       unsigned int hit0 = hitsEvt.idxs()[hit_idx0];
1344       unsigned int hit1 = hitsEvt.idxs()[hit_idx1];
1345       unsigned int hit2 = hitsEvt.idxs()[hit_idx2];
1346       unsigned int hit3 = hitsEvt.idxs()[hit_idx3];
1347       unsigned int hit4 = hitsEvt.idxs()[hit_idx4];
1348       unsigned int hit5 = hitsEvt.idxs()[hit_idx5];
1349       std::cout << "VALIDATION 'T3': "
1350                 << "T3"
1351                 << " hit0: " << hit0 << " hit1: " << hit1 << " hit2: " << hit2 << " hit3: " << hit3 << " hit4: " << hit4
1352                 << " hit5: " << hit5 << std::endl;
1353     }
1354   }
1355   std::cout << "VALIDATION nTriplets: " << nTriplets << std::endl;
1356 }
[ 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!