diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 57d6a5dd0c1..5e5d34e7f3a 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1363,6 +1363,8 @@ class VarManager : public TObject static void FillGlobalMuonRefitCov(T1 const& muontrack, T2 const& mfttrack, const C& collision, C2 const& mftcov, float* values = nullptr); template static void FillPair(T1 const& t1, T2 const& t2, float* values = nullptr); + template + static void FillPairRotation(T1 const& t1, T2 const& t2, float* values = nullptr); template static void FillPairCollision(C const& collision, T1 const& t1, T2 const& t2, float* values = nullptr); template @@ -3711,6 +3713,68 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) } } +// change_start: rotation pair +template +void VarManager::FillPairRotation(T1 const& t1, T2 const& t2, float* values) +{ + if (!values) { + values = fgValues; + } + + float m1 = o2::constants::physics::MassElectron; + float m2 = o2::constants::physics::MassElectron; + if constexpr (pairType == kDecayToMuMu) { + m1 = o2::constants::physics::MassMuon; + m2 = o2::constants::physics::MassMuon; + } + + if constexpr (pairType == kDecayToPiPi) { + m1 = o2::constants::physics::MassPionCharged; + m2 = o2::constants::physics::MassPionCharged; + } + + if constexpr (pairType == kDecayToKPi) { + m1 = o2::constants::physics::MassKaonCharged; + m2 = o2::constants::physics::MassPionCharged; + // Make the TPC information of the kaon available for pair histograms + values[kPin_leg1] = t1.tpcInnerParam(); + values[kTPCnSigmaKa_leg1] = t1.tpcNSigmaKa(); + } + + if constexpr (pairType == kElectronMuon) { + m2 = o2::constants::physics::MassMuon; + } + + double dphi = gRandom->Uniform(0., 2. * TMath::Pi()); + double rotationphi2 = t2.phi() + dphi; + + if (rotationphi2 > 2. * TMath::Pi()) + rotationphi2 -= 2. * TMath::Pi(); + + values[kCharge] = t1.sign() + t2.sign(); + values[kCharge1] = t1.sign(); + values[kCharge2] = t2.sign(); + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), rotationphi2, m2); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + values[kMass] = v12.M(); + values[kPt] = v12.Pt(); + values[kEta] = v12.Eta(); + // values[kPhi] = v12.Phi(); + values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI; + values[kRap] = -v12.Rapidity(); + double Ptot1 = TMath::Sqrt(v1.Px() * v1.Px() + v1.Py() * v1.Py() + v1.Pz() * v1.Pz()); + double Ptot2 = TMath::Sqrt(v2.Px() * v2.Px() + v2.Py() * v2.Py() + v2.Pz() * v2.Pz()); + values[kDeltaPtotTracks] = Ptot1 - Ptot2; + + values[kPt1] = t1.pt(); + values[kEta1] = t1.eta(); + values[kPhi1] = t1.phi(); + values[kPt2] = t2.pt(); + values[kEta2] = t2.eta(); + values[kPhi2] = rotationphi2; +} + template void VarManager::FillPairCollision(const C& collision, T1 const& t1, T2 const& t2, float* values) { diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 625fbd77440..adc62b0ada0 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -1281,6 +1281,10 @@ struct AnalysisSameEventPairing { Configurable fConfigQA{"cfgQA", true, "If true, fill output histograms"}; Configurable fConfigAmbiguousMuonHistograms{"cfgAmbiguousMuonHistograms", true, "If true, fill ambiguous histograms"}; + // option for TR pair fill + Configurable fConfigTRPairs{"cfgFillTRPairs", false, "If true, fill Track rotation pairs"}; + Configurable fConfigNRotations{"cfgNRotations", 20, "Number of rotations for track rotation method"}; + struct : ConfigurableGroup { Configurable url{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable grpMagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; @@ -1469,19 +1473,22 @@ struct AnalysisSameEventPairing { if (fEnableBarrelHistos) { names = { + // change define histname Form("PairsBarrelSEPM_%s", objArray->At(icut)->GetName()), Form("PairsBarrelSEPP_%s", objArray->At(icut)->GetName()), - Form("PairsBarrelSEMM_%s", objArray->At(icut)->GetName())}; - histNames += Form("%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data()); + Form("PairsBarrelSEMM_%s", objArray->At(icut)->GetName()), + Form("PairsBarrelTRPM_%s", objArray->At(icut)->GetName())}; + histNames += Form("%s;%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data(), names[3].Data()); names.push_back(Form("PairsBarrelSEPM_ambiguousextra_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsBarrelSEPP_ambiguousextra_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsBarrelSEMM_ambiguousextra_%s", objArray->At(icut)->GetName())); - histNames += Form("%s;%s;%s;", names[3].Data(), names[4].Data(), names[5].Data()); + names.push_back(Form("PairsBarrelTRPM_ambiguousextra_%s", objArray->At(icut)->GetName())); + histNames += Form("%s;%s;%s;%s;", names[4].Data(), names[5].Data(), names[6].Data(), names[7].Data()); if (fEnableBarrelMixingHistos) { names.push_back(Form("PairsBarrelMEPM_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsBarrelMEPP_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsBarrelMEMM_%s", objArray->At(icut)->GetName())); - histNames += Form("%s;%s;%s;", names[6].Data(), names[7].Data(), names[8].Data()); + histNames += Form("%s;%s;%s;", names[8].Data(), names[9].Data(), names[10].Data()); } fTrackHistNames[icut] = names; @@ -2127,6 +2134,71 @@ struct AnalysisSameEventPairing { } // end loop (pair cuts) } } // end loop (cuts) + + // edit + // rotation 20 times + if (fConfigTRPairs) { + if constexpr (TPairType == VarManager::kDecayToEE) { + twoTrackFilter = a1.isBarrelSelected_raw() & a2.isBarrelSelected_raw() & a1.isBarrelSelectedPrefilter_raw() & a2.isBarrelSelectedPrefilter_raw() & fTrackFilterMask; + + if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue + continue; + } + + auto t1 = a1.template reducedtrack_as(); + auto t2 = a2.template reducedtrack_as(); + sign1 = t1.sign(); + sign2 = t2.sign(); + if (t1.barrelAmbiguityInBunch() > 1) { + twoTrackFilter |= (static_cast(1) << 28); + } + if (t2.barrelAmbiguityInBunch() > 1) { + twoTrackFilter |= (static_cast(1) << 29); + } + if (t1.barrelAmbiguityOutOfBunch() > 1) { + twoTrackFilter |= (static_cast(1) << 30); + } + if (t2.barrelAmbiguityOutOfBunch() > 1) { + twoTrackFilter |= (static_cast(1) << 31); + } + + for (int icut = 0; icut < ncuts; icut++) { + if (twoTrackFilter & (static_cast(1) << icut)) { + isAmbiInBunch = (twoTrackFilter & (static_cast(1) << 28)) || (twoTrackFilter & (static_cast(1) << 29)); + isAmbiOutOfBunch = (twoTrackFilter & (static_cast(1) << 30)) || (twoTrackFilter & (static_cast(1) << 31)); + isUnambiguous = !(isAmbiInBunch || isAmbiOutOfBunch); + isLeg1Ambi = (twoTrackFilter & (static_cast(1) << 28) || (twoTrackFilter & (static_cast(1) << 30))); + isLeg2Ambi = (twoTrackFilter & (static_cast(1) << 29) || (twoTrackFilter & (static_cast(1) << 31))); + if constexpr (TPairType == VarManager::kDecayToEE) { + if (isLeg1Ambi && isLeg2Ambi) { + std::pair iPair(a1.reducedtrackId(), a2.reducedtrackId()); + if (fAmbiguousPairs.find(iPair) != fAmbiguousPairs.end()) { + if (fAmbiguousPairs[iPair] & (static_cast(1) << icut)) { // if this pair is already stored with this cut + isAmbiExtra = true; + } else { + fAmbiguousPairs[iPair] |= static_cast(1) << icut; + } + } else { + fAmbiguousPairs[iPair] = static_cast(1) << icut; + } + } + } + if (sign1 * sign2 < 0) { + for (int i = 0; i < fConfigNRotations.value; i++) { + VarManager::FillPairRotation(t1, t2); + if constexpr (TPairType == VarManager::kDecayToEE) { + fHistMan->FillHistClass(Form("PairsBarrelTRPM_%s", fTrackCuts[icut].Data()), VarManager::fgValues); + if (isAmbiExtra) { + fHistMan->FillHistClass(Form("PairsBarrelTRPM_ambiguousextra_%s", fTrackCuts[icut].Data()), VarManager::fgValues); + } + } + } + } + } + } + } + } + // end } // end loop over pairs of track associations VarManager::fgValues[VarManager::kNPairsPerEvent] = fNPairPerEvent; if (fEnableBarrelHistos && fConfigQA) {