Untitled

    for (auto& pair : bcsMatchedTrIdsGlobal) {
      auto globalBC = static_cast<int64_t>(pair.first);
      const auto& fwdTrackIDs = pair.second; // Forward tracks (Global with MFT)
      if (fwdTrackIDs.size() != 2) {         // ensure we have two MFT tracks
        continue;
      }

      // find the corresponding MCH-MID tracks
      auto midIt = std::find_if(bcsMatchedTrIdsMID.begin(), bcsMatchedTrIdsMID.end(),
                                [globalBC](const auto& midPair) {
                                  return midPair.first == static_cast<uint64_t>(globalBC);
                                });
      const auto* midTrackIDs = (midIt != bcsMatchedTrIdsMID.end()) ? &midIt->second : nullptr;

      // ensure MCH-MID tracks are available
      if (!midTrackIDs || midTrackIDs->size() != 2) {
        continue;
      }

      std::vector<int64_t> trkCandIDs;

      // retrieve global track eta and apply the logic
      bool firstInAcceptance = false, secondInAcceptance = false;

      const auto& trk1 = fwdTracks.iteratorAt(fwdTrackIDs[0]); // Global track 1
      const auto& trk2 = fwdTracks.iteratorAt(fwdTrackIDs[1]); // Global track 2

      // get MCH-MID tracks
      const auto& trkMID1 = fwdTracks.iteratorAt((*midTrackIDs)[0]); // MCH-MID track 1
      const auto& trkMID2 = fwdTracks.iteratorAt((*midTrackIDs)[1]); // MCH-MID track 2

      float mid_P1 = std::sqrt(trkMID1.px() * trkMID1.px() + trkMID1.py() * trkMID1.py() + trkMID1.pz() * trkMID1.pz());
      float mid_P2 = std::sqrt(trkMID2.px() * trkMID2.px() + trkMID2.py() * trkMID2.py() + trkMID2.pz() * trkMID2.pz());

      float global_tgl1 = trk1.pz() / std::sqrt(trk1.px() * trk1.px() + trk1.py() * trk1.py());
      float global_tgl2 = trk2.pz() / std::sqrt(trk2.px() * trk2.px() + trk2.py() * trk2.py());

      // calculate the azimuthal angle of the global track using px and py
      float trk1_phi =  std::atan2(trk1.py(), trk1.px());
      float trk2_phi =  std::atan2(trk2.py(), trk2.px());

      float refit_px1 = mid_P1 * sin(M_PI / 2 - atan(global_tgl1)) * cos(trk1_phi);
      float refit_py1 = mid_P1 * sin(M_PI / 2 - atan(global_tgl1)) * sin(trk1_phi);
      float refit_pz1 = mid_P1 * cos(M_PI / 2 - atan(global_tgl1));

      float refit_px2 = mid_P2 * sin(M_PI / 2 - atan(global_tgl2)) * cos(trk2_phi);
      float refit_py2 = mid_P2 * sin(M_PI / 2 - atan(global_tgl2)) * sin(trk2_phi);
      float refit_pz2 = mid_P2 * cos(M_PI / 2 - atan(global_tgl2));
      // Compute kinematic properties
      TLorentzVector p1, p2, p3, p4, p5, p6; // 4-momenta of the muons (1,2 = MFT-MCH-MID, 3,4 = MCH, 5,6 = refit)
      float massMuon = 0.10566; // GeV/c^2

      if (trk1.eta() > fMinEtaMFT && trk1.eta() < fMaxEtaMFT) {
        firstInAcceptance = true;
      }
      if (trk2.eta() > fMinEtaMFT && trk2.eta() < fMaxEtaMFT) {
        secondInAcceptance = true;
      }
      if (!firstInAcceptance && !secondInAcceptance) {
        // Both outside MFT acceptance
        p1.SetXYZM(trkMID1.px(), trkMID1.py(), trkMID1.pz(), massMuon);
        p2.SetXYZM(trkMID2.px(), trkMID2.py(), trkMID2.pz(), massMuon);
        p5.SetXYZM(trkMID1.px(), trkMID1.py(), trkMID1.pz(), massMuon);
        p6.SetXYZM(trkMID2.px(), trkMID2.py(), trkMID2.pz(), massMuon);
        trkCandIDs.insert(trkCandIDs.end(), midTrackIDs->begin(), midTrackIDs->end());
      } else if (firstInAcceptance && !secondInAcceptance) {
        // First inside, second outside
        p1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massMuon);
        p2.SetXYZM(trkMID2.px(), trkMID2.py(), trkMID2.pz(), massMuon);
        p5.SetXYZM(refit_px1, refit_py1, refit_pz1, massMuon);
        p6.SetXYZM(trkMID2.px(), trkMID2.py(), trkMID2.pz(), massMuon);
        trkCandIDs.push_back(fwdTrackIDs[0]);    // Keep first global
        trkCandIDs.push_back((*midTrackIDs)[1]); // Replace second with MCH-MID
      } else if (!firstInAcceptance && secondInAcceptance) {
        // First outside, second inside
        p1.SetXYZM(trkMID1.px(), trkMID1.py(), trkMID1.pz(), massMuon);
        p2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massMuon);
        p5.SetXYZM(trkMID1.px(), trkMID1.py(), trkMID1.pz(), massMuon);
        p6.SetXYZM(refit_px2, refit_py2, refit_pz2, massMuon);
        trkCandIDs.push_back((*midTrackIDs)[0]); // Replace first with MCH-MID
        trkCandIDs.push_back(fwdTrackIDs[1]);    // Keep second global
      } else {
        // Both inside MFT acceptance
        p1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massMuon);
        p2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massMuon);
        p5.SetXYZM(refit_px1, refit_py1, refit_pz1, massMuon);
        p6.SetXYZM(refit_px2, refit_py2, refit_pz2, massMuon);
        trkCandIDs.insert(trkCandIDs.end(), fwdTrackIDs.begin(), fwdTrackIDs.end());
      }
      p3.SetXYZM(trkMID1.px(), trkMID1.py(), trkMID1.pz(), massMuon);
      p4.SetXYZM(trkMID2.px(), trkMID2.py(), trkMID2.pz(), massMuon);
      TLorentzVector dimuonGlobal = p1 + p2;
      TLorentzVector dimuonMCHMID = p3 + p4;
      TLorentzVector dimuonRefit = p5 + p6;

      // Fill the new table, use placeholder for p5 and p6 for the moment
      dimuSel(dimuonGlobal.M(), dimuonGlobal.Pt(), dimuonGlobal.Rapidity(), dimuonGlobal.Phi(),
              dimuonMCHMID.M(), dimuonMCHMID.Pt(), dimuonMCHMID.Rapidity(), dimuonMCHMID.Phi(),
              dimuonRefit.M(), dimuonRefit.Pt(), dimuonRefit.Rapidity(), dimuonRefit.Phi());