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());
Leave a Comment