From 30fab024c13527b2824868e987a78cd3ea466273 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 12 Oct 2025 07:56:16 +0200 Subject: [PATCH 1/8] changes of psi2s and dilepton v0 --- PWGDQ/Core/CutsLibrary.cxx | 456 +++++++++++++- PWGDQ/Core/HistogramsLibrary.cxx | 64 +- PWGDQ/Core/VarManager.h | 45 +- PWGDQ/DataModel/ReducedInfoTables.h | 21 + PWGDQ/TableProducer/tableMaker.cxx | 2 + PWGDQ/TableProducer/tableMakerMC.cxx | 60 ++ .../TableProducer/tableMakerMC_withAssoc.cxx | 17 +- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 510 ++++++++++++++- PWGDQ/Tasks/tableReader.cxx | 36 +- PWGDQ/Tasks/tableReader_withAssoc.cxx | 581 +++++++++++++++++- 10 files changed, 1722 insertions(+), 70 deletions(-) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index ee152f11b6b..6690ffa56b1 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -646,6 +646,81 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + if (!nameStr.compare("pionPIDCut3")) { + cut->AddCut(GetAnalysisCut("pionQualityCut3")); + cut->AddCut(GetAnalysisCut("pionPIDnsigma")); + return cut; + } + + if (!nameStr.compare("pionPIDCut3_tof")) { + cut->AddCut(GetAnalysisCut("pionQualityCut3")); + cut->AddCut(GetAnalysisCut("pionPID_TPCnTOF")); + return cut; + } + + if (!nameStr.compare("pionPIDCut4")) { + cut->AddCut(GetAnalysisCut("pionQualityCut4")); + cut->AddCut(GetAnalysisCut("pionPIDnsigma")); + return cut; + } + + if (!nameStr.compare("pionPIDCut4_tof")) { + cut->AddCut(GetAnalysisCut("pionQualityCut4")); + cut->AddCut(GetAnalysisCut("pionPID_TPCnTOF")); + return cut; + } + + if (!nameStr.compare("trackingEff_kine_pion")) { + cut->AddCut(GetAnalysisCut("pionKineCut")); + return cut; + } + + if (!nameStr.compare("trackingEff_its_pion")) { + cut->AddCut(GetAnalysisCut("pionKineCut")); + cut->AddCut(GetAnalysisCut("ITSibany")); + return cut; + } + + if (!nameStr.compare("trackingEff_noDCA_pion")) { + cut->AddCut(GetAnalysisCut("pionQualityCut1")); + return cut; + } + + if (!nameStr.compare("trackingEff_pion")) { + cut->AddCut(GetAnalysisCut("pionQualityCut3")); + return cut; + } + + if (!nameStr.compare("PIDEff_pion_wPID")) { + cut->AddCut(GetAnalysisCut("electronStandardQualityTPCOnly")); + cut->AddCut(GetAnalysisCut("pionPIDnsigma")); + cut->AddCut(GetAnalysisCut("pidcalib_pion")); + return cut; + } + + if (!nameStr.compare("PIDEff_pion_woPID")) { + cut->AddCut(GetAnalysisCut("electronStandardQualityTPCOnly")); + cut->AddCut(GetAnalysisCut("pidcalib_pion")); + return cut; + } + + if (!nameStr.compare("PIDEff_pionMC_woPID")) { + cut->AddCut(GetAnalysisCut("pionQualityCut3")); + return cut; + } + + if (!nameStr.compare("pionDCArej")) { + cut->AddCut(GetAnalysisCut("pionDCArej")); + cut->AddCut(GetAnalysisCut("pionPIDnsigma")); + return cut; + } + + if (!nameStr.compare("pionPIDCut2")) { + cut->AddCut(GetAnalysisCut("pionQualityCut2")); + cut->AddCut(GetAnalysisCut("pionPIDnsigma")); + return cut; + } + if (!nameStr.compare("PIDCalibElectron")) { cut->AddCut(GetAnalysisCut("pidcalib_ele")); return cut; @@ -3491,6 +3566,133 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + // Psi2S to Jpsi Pi Pi cut by yiping + if (!nameStr.compare("pairX3872Minitree")) { + cut->AddCut(GetAnalysisCut("pairX3872_minitree")); + return cut; + } + + if (!nameStr.compare("pairX3872Cut8")) { + cut->AddCut(GetAnalysisCut("pairX3872_8")); + return cut; + } + if (!nameStr.compare("pairX3872Cut9")) { + cut->AddCut(GetAnalysisCut("pairX3872_9")); + return cut; + } + if (!nameStr.compare("pairX3872Cut10")) { + cut->AddCut(GetAnalysisCut("pairX3872_10")); + return cut; + } + if (!nameStr.compare("pairX3872Cut11")) { + cut->AddCut(GetAnalysisCut("pairX3872_11")); + return cut; + } + if (!nameStr.compare("pairX3872Cut12")) { + cut->AddCut(GetAnalysisCut("pairX3872_12")); + return cut; + } + if (!nameStr.compare("pairX3872Cut15")) { + cut->AddCut(GetAnalysisCut("pairX3872_15")); + return cut; + } + if (!nameStr.compare("pairX3872CutPt1")) { + cut->AddCut(GetAnalysisCut("pairX3872_13")); + cut->AddCut(GetAnalysisCut("dileptonPt_rejection")); + cut->AddCut(GetAnalysisCut("deltaR_rejection")); + return cut; + } + if (!nameStr.compare("pairX3872CutPt2")) { + cut->AddCut(GetAnalysisCut("pairX3872_14")); + cut->AddCut(GetAnalysisCut("dileptonPt_rejection")); + cut->AddCut(GetAnalysisCut("deltaR_rejection")); + return cut; + } + if (!nameStr.compare("pairX3872CutPt3")) { + cut->AddCut(GetAnalysisCut("pairX3872_14")); + cut->AddCut(GetAnalysisCut("dileptonPt_rejection")); + cut->AddCut(GetAnalysisCut("deltaR_2sigrejection")); + return cut; + } + if (!nameStr.compare("pairX3872CutPt4")) { + cut->AddCut(GetAnalysisCut("pairX3872_14")); + cut->AddCut(GetAnalysisCut("dileptonPt_rejection")); + cut->AddCut(GetAnalysisCut("deltaR_15sigrejection")); + return cut; + } + // following cuts are used for prompt-non-prompt Psi2S to ee + if (!nameStr.compare("electronSelection1_yiping")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine4")); + cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); + cut->AddCut(GetAnalysisCut("dcaCut1_ionut")); + cut->AddCut(GetAnalysisCut("jpsi_TPCPID_debug4")); + return cut; + } + if (!nameStr.compare("electronSelection2_yiping")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); + cut->AddCut(GetAnalysisCut("dcaCut1_ionut")); + cut->AddCut(GetAnalysisCut("jpsi_TPCPID_debug4")); + return cut; + } + if (!nameStr.compare("electronSelection3_yiping")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); + cut->AddCut(GetAnalysisCut("dcaCut1_ionut")); + cut->AddCut(GetAnalysisCut("jpsi_TPCPID_debug6")); + return cut; + } + if (!nameStr.compare("electronSelection4_yiping")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); + cut->AddCut(GetAnalysisCut("dcaCut1_ionut")); + cut->AddCut(GetAnalysisCut("electronPIDnsigmaMedium")); + return cut; + } + + if (!nameStr.compare("trackingEff_ele")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); + return cut; + } + + if (!nameStr.compare("trackingEff_ele_wDCA")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); + cut->AddCut(GetAnalysisCut("dcaCut1_ionut")); + return cut; + } + + if (!nameStr.compare("PIDEff_ele4_wPID")) { + cut->AddCut(GetAnalysisCut("electronPIDnsigmaMedium")); + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityTPCOnly")); + cut->AddCut(GetAnalysisCut("pidcalib_ele")); + return cut; + } + + if (!nameStr.compare("PIDEff_ele4_woPID")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityTPCOnly")); + cut->AddCut(GetAnalysisCut("pidcalib_ele")); + return cut; + } + + if (!nameStr.compare("pairX3872hasVtx")) { + cut->AddCut(GetAnalysisCut("pairX3872_hasVtx")); + return cut; + } + + if (!nameStr.compare("pairB0test")) { + cut->AddCut(GetAnalysisCut("pair_b0test")); + return cut; + } + + if (!nameStr.compare("pairBstest")) { + cut->AddCut(GetAnalysisCut("pair_bstest")); + return cut; + } + if (!nameStr.compare("DipionPairCut1")) { cut->AddCut(GetAnalysisCut("DipionMassCut1")); return cut; @@ -4685,6 +4887,11 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("pionKineCut")) { + cut->AddCut(VarManager::kPt, 0.15, 1000.0); + return cut; + } + if (!nameStr.compare("pionQualityCut1")) { cut->AddCut(VarManager::kPt, 0.15, 1000.0); cut->AddCut(VarManager::kIsITSibAny, 0.5, 1.5); @@ -4694,10 +4901,37 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) if (!nameStr.compare("pionQualityCut2")) { cut->AddCut(VarManager::kPt, 0.15, 1000.0); - cut->AddCut(VarManager::kEta, -0.9, 0.9); + // cut->AddCut(VarManager::kEta, -0.9, 0.9); cut->AddCut(VarManager::kIsITSibAny, 0.5, 1.5); cut->AddCut(VarManager::kTPCncls, 90, 161); cut->AddCut(VarManager::kTPCnclsCR, 70, 161); + cut->AddCut(VarManager::kTrackDCAxy, -0.5, 0.5); + cut->AddCut(VarManager::kTrackDCAz, -0.5, 0.5); + return cut; + } + + if (!nameStr.compare("pionQualityCut3")) { + cut->AddCut(VarManager::kPt, 0.15, 1000.0); + cut->AddCut(VarManager::kIsITSibAny, 0.5, 1.5); + cut->AddCut(VarManager::kTPCncls, 70, 161); + cut->AddCut(VarManager::kTrackDCAxy, -0.5, 0.5); + cut->AddCut(VarManager::kTrackDCAz, -0.5, 0.5); + return cut; + } + + if (!nameStr.compare("pionQualityCut4")) { + cut->AddCut(VarManager::kPt, 0.1, 1000.0); + cut->AddCut(VarManager::kIsITSibAny, 0.5, 1.5); + cut->AddCut(VarManager::kTPCncls, 70, 161); + cut->AddCut(VarManager::kTrackDCAxy, -0.5, 0.5); + cut->AddCut(VarManager::kTrackDCAz, -0.5, 0.5); + return cut; + } + + if (!nameStr.compare("pionDCArej")) { + cut->AddCut(VarManager::kTPCncls, 120, 161); + cut->AddCut(VarManager::kTrackDCAxy, -0.1, 0.1, true); + cut->AddCut(VarManager::kTrackDCAz, -0.1, 0.1, true); return cut; } @@ -5572,6 +5806,8 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) if (!nameStr.compare("kaonPIDnsigma")) { cut->AddCut(VarManager::kTPCnSigmaKa, -3.0, 3.0); + cut->AddCut(VarManager::kTOFnSigmaKa, -3.0, 3.0, false, VarManager::kHasTOF, 0.5, 1.5, false); + cut->AddCut(VarManager::kTOFnSigmaPi, 3.0, 3000.0, false, VarManager::kHasTOF, 0.5, 1.5, false, VarManager::kPin, 0.4, 3, false); return cut; } @@ -5624,7 +5860,7 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) if (!nameStr.compare("pionPID_TPCnTOF")) { cut->AddCut(VarManager::kTPCnSigmaPi, -3.0, 3.0); - cut->AddCut(VarManager::kTOFnSigmaPi, -3.0, 3.0); + cut->AddCut(VarManager::kTOFnSigmaPi, -3.0, 3.0, false, VarManager::kHasTOF, 0.5, 1.5, false); return cut; } @@ -6503,6 +6739,222 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } +<<<<<<< HEAD +======= + if (!nameStr.compare("pairX3872_minitree")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.5, 0.5); + cut->AddCut(VarManager::kDeltaR, 0.0, 3.5); + cut->AddCut(VarManager::kDitrackMass, 0.0, 0.8); + cut->AddCut(VarManager::kQuadPt, 0.0, 1000.0); + return cut; + } + + if (!nameStr.compare("pairX3872_8")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.5, 0.5); + cut->AddCut(VarManager::kDeltaR1, 0.0, 2.0); + cut->AddCut(VarManager::kDeltaR2, 0.0, 2.0); + cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); + cut->AddCut(VarManager::kDitrackMass, 0.0, 0.8); + cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); +<<<<<<< HEAD + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + return cut; + } + + if (!nameStr.compare("pairX3872_9")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.4, 0.4); + cut->AddCut(VarManager::kDeltaR1, 0.0, 1.0); + cut->AddCut(VarManager::kDeltaR2, 0.0, 1.0); + cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); + cut->AddCut(VarManager::kDitrackMass, 0.2, 0.6); + cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); +<<<<<<< HEAD + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + return cut; + } + + if (!nameStr.compare("pairX3872_10")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.2, 0.2); + // cut->AddCut(VarManager::kDeltaR1, 0.0, 1.0); + // cut->AddCut(VarManager::kDeltaR2, 0.0, 1.0); + cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); + cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); + // cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); +<<<<<<< HEAD + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + return cut; + } + + if (!nameStr.compare("pairX3872_11")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.2, 0.2); + // cut->AddCut(VarManager::kDeltaR1, 0.0, 1.0); + // cut->AddCut(VarManager::kDeltaR2, 0.0, 1.0); + cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); + cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); + // cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); +<<<<<<< HEAD + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + cut->AddCut(VarManager::kPtOverPairPt, 1.0, 1.5); + return cut; + } + + if (!nameStr.compare("pairX3872_12")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.2, 0.2); + // cut->AddCut(VarManager::kDeltaR1, 0.0, 1.0); + // cut->AddCut(VarManager::kDeltaR2, 0.0, 1.0); + cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); + cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); + // cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); +<<<<<<< HEAD + // cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + // cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + // cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + // cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + cut->AddCut(VarManager::kPtOverPairPt, 1.0, 1.5); + return cut; + } + + if (!nameStr.compare("pairX3872_13")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.2, 0.2); +<<<<<<< HEAD + cut->AddCut(VarManager::kDeltaR, 0.0, 5.0); + cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); + cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + return cut; + } + + if (!nameStr.compare("pairX3872_14")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.4, 0.2); +<<<<<<< HEAD + cut->AddCut(VarManager::kDeltaR, 0.0, 5.0); + cut->AddCut(VarManager::kDitrackMass, 0.4, 0.58); + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); + cut->AddCut(VarManager::kDitrackMass, 0.4, 0.58); + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + return cut; + } + + if (!nameStr.compare("pairX3872_15")) { + cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); + cut->AddCut(VarManager::kQ, -0.4, 0.2); + cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); + cut->AddCut(VarManager::kDitrackMass, 0.4, 0.58); +<<<<<<< HEAD + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); + cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + cut->AddCut(VarManager::kPtOverPairPt, 1.0, 1.5); + return cut; + } + + if (!nameStr.compare("pairX3872_hasVtx")) { +<<<<<<< HEAD + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); +======= + cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + cut->AddCut(VarManager::kVertexingQuadProcCode,0.5, 1.5); + return cut; + } + + if (!nameStr.compare("dileptonPt_rejection")) { + TF1* fPairPtLow = new TF1("f1", "pol1", 0., 20.); + fPairPtLow->SetParameters(-0.419607, 1.07925); + TF1* fPairPtHigh = new TF1("f2", "pol1", 0., 20.); + fPairPtHigh->SetParameters(0.515923, 1.34061); + cut->AddCut(VarManager::kQuadPt, fPairPtLow, fPairPtHigh, false, VarManager::kPairPt, 0., 20.); + return cut; + } + + if (!nameStr.compare("deltaR_rejection")) { + TF1* fDeltaRLow = new TF1("f1", "[0]*exp([1]*x)+[2]", 0., 5.); + fDeltaRLow->SetParameters(12.3641, -6.45546, 0.332389); + TF1* fDeltaRHigh = new TF1("f2", "[0]*exp([1]*x)+[2]", 0., 5.); + fDeltaRHigh->SetParameters(28.7518, -1.25672, 1.06815); + cut->AddCut(VarManager::kQuadPt, fDeltaRLow, fDeltaRHigh, false, VarManager::kDeltaR, 0., 5.); + return cut; + } + + if (!nameStr.compare("deltaR_2sigrejection")) { + TF1* fDeltaRLow = new TF1("f1", "[0]*exp([1]*x)+[2]", 0., 5.); + fDeltaRLow->SetParameters(17.5661, -3.39083, 0.340562); + TF1* fDeltaRHigh = new TF1("f2", "[0]*exp([1]*x)+[2]", 0., 5.); + fDeltaRHigh->SetParameters(28.8933, -1.51157, 1.10834); + cut->AddCut(VarManager::kQuadPt, fDeltaRLow, fDeltaRHigh, false, VarManager::kDeltaR, 0., 5.); + cut->AddCut(VarManager::kDeltaR, 0.0, 3.0); + return cut; + } + + if (!nameStr.compare("deltaR_15sigrejection")) { + TF1* fDeltaRLow = new TF1("f1", "[0]*exp([1]*x)+[2]", 0., 5.); + fDeltaRLow->SetParameters(21.1786, -3.30186, 0.80823); + TF1* fDeltaRHigh = new TF1("f2", "[0]*exp([1]*x)+[2]", 0., 5.); + fDeltaRHigh->SetParameters(28.8547, -1.67054, 1.1059); + cut->AddCut(VarManager::kQuadPt, fDeltaRLow, fDeltaRHigh, false, VarManager::kDeltaR, 0., 5.); + cut->AddCut(VarManager::kDeltaR, 0.0, 3.0); + return cut; + } + + if (!nameStr.compare("pair_b0test")) { +<<<<<<< HEAD + cut->AddCut(VarManager::kDitrackMass, 0.49, 0.502); +======= + cut->AddCut(VarManager::kDitrackMass, 0.49, 0.51); + return cut; + } + + if (!nameStr.compare("pair_bstest")) { + cut->AddCut(VarManager::kDitrackMass, 1.0, 1.04); +>>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) + return cut; + } + +>>>>>>> 9942c0482 (updated psi2s analysis and dilepton v0) if (!nameStr.compare("pairPtLow1")) { cut->AddCut(VarManager::kPt, 2.0, 1000.0); return cut; diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index b01443375fc..64d38935f2c 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -560,6 +560,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h } else { hm->AddHistogram(histClass, "TPCdedx_pIN", "TPC dE/dx vs pIN", false, 100, 0.0, 20.0, VarManager::kPin, 150, 0.0, 150., VarManager::kTPCsignal); hm->AddHistogram(histClass, "TPCnSigEle_pIN", "TPC n-#sigma(e) vs pIN", false, 100, 0.0, 20.0, VarManager::kPin, 100, -5.0, 5.0, VarManager::kTPCnSigmaEl); + hm->AddHistogram(histClass, "TPCnSigEle_pT", "TPC n-#sigma(e) vs p_{T}", false, 100, 0.0, 20.0, VarManager::kPt, 100, -5.0, 5.0, VarManager::kTPCnSigmaEl); hm->AddHistogram(histClass, "TPCnSigEle_occupancy", "TPC n-#sigma(e) vs occupancy", false, 200, 0., 20000., VarManager::kTrackOccupancyInTimeRange, 100, -5.0, 5.0, VarManager::kTPCnSigmaEl); hm->AddHistogram(histClass, "TPCnSigEle_timeFromSOR", "TPC n-#sigma(e) vs time from SOR", true, 10000, 0.0, 1000.0, VarManager::kTimeFromSOR, 10, -5.0, 5.0, VarManager::kTPCnSigmaEl); hm->AddHistogram(histClass, "TPCnSigEle_occupTPCcontribLongA", "TPC n-#sigma(e) vs pileup n-contrib, long time range A-side", false, 20, 0.0, 10000.0, VarManager::kNTPCcontribLongA, 200, -5.0, 5.0, VarManager::kTPCnSigmaEl); @@ -575,6 +576,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "TPCnSigEle_occupTPCmeanTimeShortC", "TPC n-#sigma(e) vs pileup mean time, short time range, C-side", false, 20, -15.0, 15.0, VarManager::kNTPCmeanTimeShortC, 200, -5.0, 5.0, VarManager::kTPCnSigmaEl); hm->AddHistogram(histClass, "TPCnSigEle_occupTPCmedianTimeShortC", "TPC n-#sigma(e) vs pileup mean time, short time range, C-side", false, 20, -15.0, 15.0, VarManager::kNTPCmedianTimeShortC, 200, -5.0, 5.0, VarManager::kTPCnSigmaEl); hm->AddHistogram(histClass, "TPCnSigPi_pIN", "TPC n-#sigma(#pi) vs pIN", false, 100, 0.0, 10.0, VarManager::kPin, 100, -5.0, 5.0, VarManager::kTPCnSigmaPi); + hm->AddHistogram(histClass, "TPCnSigPi_pT", "TPC n-#sigma(#pi) vs p_{T}", false, 100, 0.0, 20.0, VarManager::kPt, 100, -5.0, 5.0, VarManager::kTPCnSigmaPi); hm->AddHistogram(histClass, "TPCnSigPi_timeFromSOR", "TPC n-#sigma(#pi) vs time from SOR", true, 1000, 0.0, 1000.0, VarManager::kTimeFromSOR, 10, -5.0, 5.0, VarManager::kTPCnSigmaPi); hm->AddHistogram(histClass, "TPCnSigPi_eta", "TPC n-#sigma(#pi) vs #eta", false, 20, -1.0, 1.0, VarManager::kEta, 200, -5.0, 5.0, VarManager::kTPCnSigmaPi); hm->AddHistogram(histClass, "TPCnSigPi_etaPin_prof", " vs (#eta,p_{IN}), --s--", true, 20, -1.0, 1.0, VarManager::kEta, 20, 0.0, 10.0, VarManager::kPin, 10, -5.0, 5.0, VarManager::kTPCnSigmaPi); @@ -587,6 +589,29 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "TPCnSigPr_pIN", "TPC n-#sigma(p) vs pIN", false, 100, 0.0, 10.0, VarManager::kPin, 100, -5.0, 5.0, VarManager::kTPCnSigmaPr); hm->AddHistogram(histClass, "TPCnSigPr_timeFromSOR", "TPC n-#sigma(p) vs time from SOR", true, 10000, 0.0, 1000.0, VarManager::kTimeFromSOR, 10, -5.0, 5.0, VarManager::kTPCnSigmaPr); hm->AddHistogram(histClass, "TPCnSigPr_occupancy", "TPC n-#sigma(p) vs. occupancy", false, 200, 0., 20000., VarManager::kTrackOccupancyInTimeRange, 100, -5.0, 5.0, VarManager::kTPCnSigmaPr); + if (subGroupStr.Contains("eff")) { + const int kNvarsPID = 3; + const int kTPCnsigmaNbins = 70; + double tpcNsigmaBinLims[kTPCnsigmaNbins + 1]; + for (int i = 0; i <= kTPCnsigmaNbins; ++i) + tpcNsigmaBinLims[i] = -7.0 + 0.2 * i; + + const int kPinEleNbins = 20; + double pinEleBinLims[kPinEleNbins + 1] = {0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 8.0, 10.0, 12.0, 16.0, 20.0}; + + const int kEtaNbins = 9; + double etaBinLimsI[kEtaNbins + 1] = {-0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9}; + + TArrayD nSigBinLimits[kNvarsPID]; + nSigBinLimits[0] = TArrayD(kTPCnsigmaNbins + 1, tpcNsigmaBinLims); + nSigBinLimits[1] = TArrayD(kPinEleNbins + 1, pinEleBinLims); + nSigBinLimits[2] = TArrayD(kEtaNbins + 1, etaBinLimsI); + + int varsPIDnSigEle[kNvarsPID] = {VarManager::kTPCnSigmaEl, VarManager::kPin, VarManager::kEta}; + int varsPIDnSigPi[kNvarsPID] = {VarManager::kTPCnSigmaPi, VarManager::kPin, VarManager::kEta}; + hm->AddHistogram(histClass, "nSigmaTPCelectron", "TPC n_{#sigma}(e) Vs Pin Vs Eta", kNvarsPID, varsPIDnSigEle, nSigBinLimits); + hm->AddHistogram(histClass, "nSigmaTPCpion", "TPC n_{#sigma}(#pi) Vs Pin Vs Eta", kNvarsPID, varsPIDnSigPi, nSigBinLimits); + } if (subGroupStr.Contains("tpcpid_corr")) { hm->AddHistogram(histClass, "TPCnSigEl_Corr_pIN", "TPC n-#sigma(e) Corr. vs pIN", false, 100, 0.0, 10.0, VarManager::kPin, 100, -5.0, 5.0, VarManager::kTPCnSigmaEl_Corr); hm->AddHistogram(histClass, "TPCnSigPi_Corr_pIN", "TPC n-#sigma(#pi) Corr. vs pIN", false, 100, 0.0, 10.0, VarManager::kPin, 100, -5.0, 5.0, VarManager::kTPCnSigmaPi_Corr); @@ -1929,25 +1954,32 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h } if (!groupStr.CompareTo("dilepton-dihadron")) { if (subGroupStr.Contains("xtojpsipipi") || subGroupStr.Contains("psi2stojpsipipi")) { - hm->AddHistogram(histClass, "hMass_X3872", "", false, 1000, 3.0, 5.0, VarManager::kQuadMass); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_X3872", "", false, 1000, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass); - hm->AddHistogram(histClass, "hPt_X3872", "", false, 150, 0.0, 15.0, VarManager::kQuadPt); - hm->AddHistogram(histClass, "hMass_Pt_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 150, 0.0, 15.0, VarManager::kQuadPt); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_Pt_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 150, 0.0, 15.0, VarManager::kQuadPt); - hm->AddHistogram(histClass, "hCostheta_Jpsi_Dihadron", "", false, 100, -1.0, 1.0, VarManager::kCosthetaDileptonDitrack); - hm->AddHistogram(histClass, "hPtDilepton_PtDihadron", "", false, 150, 0, 15.0, VarManager::kPairPt, 100, 0, 10, VarManager::kDitrackPt); - hm->AddHistogram(histClass, "hPtDilepton_MassDihadron", "", false, 150, 0, 15.0, VarManager::kPairPt, 150, 0.0, 3.0, VarManager::kDitrackMass); + // hm->AddHistogram(histClass, "hMass_X3872", "", false, 1000, 3.0, 5.0, VarManager::kQuadMass); + hm->AddHistogram(histClass, "hMass_defaultDileptonMass_X3872", "", false, 1500, 3.0, 6.0, VarManager::kQuadDefaultDileptonMass); + hm->AddHistogram(histClass, "hPt_Y_X3872", "", false, 15, 0.0, 15.0, VarManager::kQuadPt, 100, -5.0, 5.0, VarManager::kRap); + hm->AddHistogram(histClass, "hPt_X3872", "", false, 500, 0.0, 15.0, VarManager::kQuadPt); + hm->AddHistogram(histClass, "hMass_defaultDileptonMass_Pt_X3872", "", false, 500, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 150, 0.0, 15.0, VarManager::kQuadPt); + // hm->AddHistogram(histClass, "hCostheta_Jpsi_Dihadron", "", false, 100, -1.0, 1.0, VarManager::kCosthetaDileptonDitrack); + // hm->AddHistogram(histClass, "hPtDilepton_PtDihadron", "", false, 150, 0, 15.0, VarManager::kPairPt, 100, 0, 10, VarManager::kDitrackPt); + hm->AddHistogram(histClass, "hPtDilepton_PtX3872", "", false, 150, 0, 15.0, VarManager::kPairPt, 150, 0.0, 15.0, VarManager::kQuadPt); hm->AddHistogram(histClass, "hQ_X3872", "", false, 300, -3.0, 3.0, VarManager::kQ); hm->AddHistogram(histClass, "hDeltaR1_X3872", "", false, 100, 0.0, 10.0, VarManager::kDeltaR1); hm->AddHistogram(histClass, "hDeltaR2_X3872", "", false, 100, 0.0, 10.0, VarManager::kDeltaR2); hm->AddHistogram(histClass, "hDeltaR_X3872", "", false, 100, 0.0, 10.0, VarManager::kDeltaR); - hm->AddHistogram(histClass, "hMass_Q_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 300, -3.0, 3.0, VarManager::kQ); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_Q_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 300, -3.0, 3.0, VarManager::kQ); - hm->AddHistogram(histClass, "hMass_DeltaR1_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kDeltaR1); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_DeltaR1_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 100, 0.0, 10.0, VarManager::kDeltaR1); - hm->AddHistogram(histClass, "hMass_DeltaR2_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kDeltaR2); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_DeltaR2_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 100, 0.0, 10.0, VarManager::kDeltaR2); - hm->AddHistogram(histClass, "hMass_X3872_MassDihadron", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 150, 0.0, 3.0, VarManager::kDitrackMass); + hm->AddHistogram(histClass, "DeltaR_Pt_X3872", "", false, 100, 0.0, 10.0, VarManager::kDeltaR, 150, 0.0, 15.0, VarManager::kQuadPt); + hm->AddHistogram(histClass, "hPtOverPairPt_X3872", "", false, 100, 0.0, 2.0, VarManager::kPtOverPairPt); + hm->AddHistogram(histClass, "hPtOverPairPt_Pt_X3872", "", false, 100, 0.0, 2.0, VarManager::kPtOverPairPt, 150, 0.0, 15.0, VarManager::kQuadPt); + hm->AddHistogram(histClass, "hMassDihadron_X3872", "", false, 300, 0.0, 3.0, VarManager::kDitrackMass); + hm->AddHistogram(histClass, "hMassDihadron_fine", "", false, 100, 1.0, 1.1, VarManager::kDitrackMass); + hm->AddHistogram(histClass, "hMassDihadron_Pt_X3872", "", false, 300, 0.0, 3.0, VarManager::kDitrackMass, 150, 0.0, 15.0, VarManager::kQuadPt); + hm->AddHistogram(histClass, "hMassDihadron_PtDilepton", "", false, 300, 0.0, 3.0, VarManager::kDitrackMass, 150, 0.0, 15.0, VarManager::kPairPt); + hm->AddHistogram(histClass, "hMassDihadron_ptDihadron", "", false, 300, 0.0, 3.0, VarManager::kDitrackMass, 150, 0.0, 15.0, VarManager::kDitrackPt); + hm->AddHistogram(histClass, "hMassDihadron_hQ", "", false, 300, 0.0, 3.0, VarManager::kDitrackMass, 300, -3.0, 3.0, VarManager::kQ); + // hm->AddHistogram(histClass, "hMass_Q_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 300, -3.0, 3.0, VarManager::kQ); + // hm->AddHistogram(histClass, "hMass_DeltaR1_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kDeltaR1); + // hm->AddHistogram(histClass, "hMass_DeltaR2_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kDeltaR2); + hm->AddHistogram(histClass, "hMass_DeltaR_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kDeltaR); + // hm->AddHistogram(histClass, "hMass_X3872_MassDihadron", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 150, 0.0, 3.0, VarManager::kDitrackMass); hm->AddHistogram(histClass, "hMass_defaultDileptonMass_X3872_MassDihadron", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 150, 0.0, 3.0, VarManager::kDitrackMass); hm->AddHistogram(histClass, "hRap_X3872", "", false, 1000, 0.0, 5.0, VarManager::kRap); hm->AddHistogram(histClass, "hMass_Rap_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 1000, 0.0, 5.0, VarManager::kRap); @@ -1968,6 +2000,8 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "Lz", "", false, 1000, -0.2, 0.2, VarManager::kVertexingLz); hm->AddHistogram(histClass, "Lxy", "", false, 1000, -0.2, 0.2, VarManager::kVertexingLxy); hm->AddHistogram(histClass, "Lxyz", "", false, 1000, -0.2, 0.2, VarManager::kVertexingLxyz); + hm->AddHistogram(histClass, "Lxy_Mass", "", false, 500, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 1000, -0.2, 0.2, VarManager::kVertexingLxy); + hm->AddHistogram(histClass, "Lxyz_Mass", "", false, 500, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 1000, -0.2, 0.2, VarManager::kVertexingLxyz); hm->AddHistogram(histClass, "Tauz", "", false, 4000, -0.01, 0.01, VarManager::kVertexingTauz); hm->AddHistogram(histClass, "Tauxy", "", false, 4000, -0.01, 0.01, VarManager::kVertexingTauxy); hm->AddHistogram(histClass, "LxyzErr", "", false, 100, 0.0, 0.2, VarManager::kVertexingLxyzErr); diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 1f5f7d0c59c..63261d3d095 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -108,6 +108,7 @@ class VarManager : public TObject ReducedEventMultExtra = BIT(19), CollisionQvectCentr = BIT(20), RapidityGapFilter = BIT(21), + ReducedEventCollInfo = BIT(22), Track = BIT(0), TrackCov = BIT(1), TrackExtra = BIT(2), @@ -136,7 +137,8 @@ class VarManager : public TObject ReducedMuonCollInfo = BIT(25), // TODO: remove it once new reduced data tables are produced for dielectron with ReducedTracksBarrelInfo MuonRealign = BIT(26), MuonCovRealign = BIT(27), - MFTCov = BIT(28) + MFTCov = BIT(28), + MCTPCtuneOnData = BIT(29) }; enum PairCandidateType { @@ -147,6 +149,9 @@ class VarManager : public TObject kElectronMuon, // e.g. Electron - muon correlations kBcToThreeMuons, // e.g. Bc -> mu+ mu- mu+ kBtoJpsiEEK, // e.g. B+ -> e+ e- K+ + kBtoJpsiEEK0S, // e.g. B0 -> e+ e- K0s + kB0toJpsiEEPiPi, // e.g. B0 -> e+ e- pi+ pi- + kBstoJpsiEEPhi, // e.g. Bs -> e+ e- phi kJpsiEEProton, // e.g. Jpsi-proton correlation, Jpsi to e+e- kXtoJpsiPiPi, // e.g. X(3872) -> J/psi pi+ pi- kPsi2StoJpsiPiPi, // e.g. Psi(2S) -> J/psi pi+ pi- @@ -2381,6 +2386,10 @@ void VarManager::FillTrack(T const& track, float* values) values[kTRDPattern] = track.trdPattern(); values[kTPCsignal] = track.tpcSignal(); + if constexpr ((fillMap & MCTPCtuneOnData) > 0) { + // TPC signal without the gain correction for MC TPC tune on data + values[kTPCsignal] = track.mcTunedTPCSignal(); + } values[kTRDsignal] = track.trdSignal(); values[kDetectorMap] = track.detectorMap(); @@ -3905,7 +3914,7 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, values[kVertexingLxyProjected] = values[kVertexingLxyProjected] / TMath::Sqrt((KFGeoTwoProng.GetPx() * KFGeoTwoProng.GetPx()) + (KFGeoTwoProng.GetPy() * KFGeoTwoProng.GetPy())); values[kVertexingLxyzProjected] = (dxPair2PV * KFGeoTwoProng.GetPx()) + (dyPair2PV * KFGeoTwoProng.GetPy()) + (dzPair2PV * KFGeoTwoProng.GetPz()); values[kVertexingLxyzProjected] = values[kVertexingLxyzProjected] / TMath::Sqrt((KFGeoTwoProng.GetPx() * KFGeoTwoProng.GetPx()) + (KFGeoTwoProng.GetPy() * KFGeoTwoProng.GetPy()) + (KFGeoTwoProng.GetPz() * KFGeoTwoProng.GetPz())); - values[kVertexingTauxyProjected] = values[kVertexingLxyProjected] * KFGeoTwoProng.GetMass() / (KFGeoTwoProng.GetPt()); + values[kVertexingTauxyProjected] = KFGeoTwoProng.GetPseudoProperDecayTime(KFPV, KFGeoTwoProng.GetMass()) / (o2::constants::physics::LightSpeedCm2NS); values[kVertexingTauxyProjectedPoleJPsiMass] = values[kVertexingLxyProjected] * o2::constants::physics::MassJPsi / (KFGeoTwoProng.GetPt()); values[kVertexingTauxyProjectedNs] = values[kVertexingTauxyProjected] / o2::constants::physics::LightSpeedCm2NS; values[kVertexingTauzProjected] = values[kVertexingLzProjected] * KFGeoTwoProng.GetMass() / TMath::Abs(KFGeoTwoProng.GetPz()); @@ -4288,7 +4297,7 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton procCode = VarManager::fgFitterThreeProngFwd.process(pars1, pars2, pars3); procCodeJpsi = VarManager::fgFitterTwoProngFwd.process(pars1, pars2); - } else if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi) && trackHasCov) { + } else if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi || candidateType == kB0toJpsiEEPiPi) && trackHasCov) { if constexpr ((candidateType == kBtoJpsiEEK) && trackHasCov) { mlepton1 = o2::constants::physics::MassElectron; mlepton2 = o2::constants::physics::MassElectron; @@ -4297,6 +4306,10 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton mlepton1 = o2::constants::physics::MassKaonCharged; mlepton2 = o2::constants::physics::MassPionCharged; mtrack = o2::constants::physics::MassPionCharged; + } else if constexpr ((candidateType) == kB0toJpsiEEPiPi && trackHasCov) { + mlepton1 = o2::constants::physics::MassElectron; + mlepton2 = o2::constants::physics::MassElectron; + mtrack = o2::constants::physics::MassKaonNeutral; } std::array lepton1pars = {lepton1.y(), lepton1.z(), lepton1.snp(), lepton1.tgl(), lepton1.signed1Pt()}; std::array lepton1covs = {lepton1.cYY(), lepton1.cZY(), lepton1.cZZ(), lepton1.cSnpY(), lepton1.cSnpZ(), @@ -4371,7 +4384,7 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton o2::dataformats::VertexBase primaryVertex = {std::move(vtxXYZ), std::move(vtxCov)}; auto covMatrixPV = primaryVertex.getCov(); - if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi) && trackHasCov) { + if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi || candidateType == kB0toJpsiEEPiPi) && trackHasCov) { secondaryVertex = fgFitterThreeProngBarrel.getPCACandidate(); covMatrixPCA = fgFitterThreeProngBarrel.calcPCACovMatrixFlat(); } else if constexpr (candidateType == kBcToThreeMuons && muonHasCov) { @@ -5160,6 +5173,11 @@ void VarManager::FillDileptonTrackTrack(T1 const& dilepton, T2 const& hadron1, T hadronMass1 = o2::constants::physics::MassElectron; hadronMass2 = o2::constants::physics::MassElectron; } + if (candidateType == kBstoJpsiEEPhi) { + defaultDileptonMass = 3.096; + hadronMass1 = o2::constants::physics::MassKaonCharged; + hadronMass2 = o2::constants::physics::MassKaonCharged; + } ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass()); ROOT::Math::PtEtaPhiMVector v2(hadron1.pt(), hadron1.eta(), hadron1.phi(), hadronMass1); @@ -5281,8 +5299,27 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l } else { Vec3D secondaryVertex; std::array covMatrixPCA; + if (candidateType == kB0toJpsiEEPiPi) { + if (!fgFitterTwoProngBarrel.process(pars3,pars4)) + return; + else { + o2::track::TrackParCov parsV0 = fgFitterTwoProngBarrel.createParentTrackParCov(0); + procCodeDileptonTrackTrack = fgFitterThreeProngBarrel.process(pars1, pars2, parsV0); + secondaryVertex = fgFitterThreeProngBarrel.getPCACandidate(); + covMatrixPCA = fgFitterThreeProngBarrel.calcPCACovMatrixFlat(); + values[kVertexingChi2PCA] = fgFitterThreeProngBarrel.getChi2AtPCACandidate(); + } + } else if (candidateType == kXtoJpsiPiPi || candidateType == kPsi2StoJpsiPiPi) { secondaryVertex = fgFitterFourProngBarrel.getPCACandidate(); covMatrixPCA = fgFitterFourProngBarrel.calcPCACovMatrixFlat(); + values[kVertexingChi2PCA] = fgFitterFourProngBarrel.getChi2AtPCACandidate(); + } else { + return; // unsupported candidate type + } + values[kVertexingProcCode] = procCodeDilepton; + values[kVertexingQuadProcCode] = procCodeDileptonTrackTrack; + // secondaryVertex = fgFitterFourProngBarrel.getPCACandidate(); + // covMatrixPCA = fgFitterFourProngBarrel.calcPCACovMatrixFlat(); o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; diff --git a/PWGDQ/DataModel/ReducedInfoTables.h b/PWGDQ/DataModel/ReducedInfoTables.h index dd92f5a9c76..a714330284e 100644 --- a/PWGDQ/DataModel/ReducedInfoTables.h +++ b/PWGDQ/DataModel/ReducedInfoTables.h @@ -698,6 +698,17 @@ DECLARE_SOA_COLUMN(Tauxy, tauxy, float); DECLARE_SOA_COLUMN(TauxyErr, tauxyErr, float); //! Error on transverse pseudo-proper time of lepton pair (in ns) DECLARE_SOA_COLUMN(Lz, lz, float); //! Longitudinal projection of decay length DECLARE_SOA_COLUMN(Lxy, lxy, float); //! Transverse projection of decay length +DECLARE_SOA_COLUMN(LzErr, lzerr, float); //! Lz over error +DECLARE_SOA_COLUMN(LxyErr, lxyerr, float); +DECLARE_SOA_COLUMN(TauzProj, tauzproj, float); +DECLARE_SOA_COLUMN(TauxyProj, tauxyproj, float); +DECLARE_SOA_COLUMN(TauzProjErr, tauzprojerr, float); +DECLARE_SOA_COLUMN(TauxyProjErr, tauxyprojerr, float); +DECLARE_SOA_COLUMN(LzProj, lzproj, float); +DECLARE_SOA_COLUMN(LxyProj, lxyproj, float); +DECLARE_SOA_COLUMN(Lxyz, lxyz, float); +DECLARE_SOA_COLUMN(Tauxyz, tauxyz, float); +DECLARE_SOA_COLUMN(TauxyzErr, tauxyzerr, float); DECLARE_SOA_COLUMN(Chi2pca, chi2pca, float); //! Chi2 for PCA of the dilepton DECLARE_SOA_COLUMN(CosPointingAngle, cosPointingAngle, float); //! Cosine of the pointing angle DECLARE_SOA_COLUMN(U2Q2, u2q2, float); //! Scalar product between unitary vector with event flow vector (harmonic 2) @@ -837,6 +848,16 @@ DECLARE_SOA_TABLE(DielectronsAll, "AOD", "RTDIELECTRONALL", //! reducedpair::Lz, reducedpair::Lxy); +DECLARE_SOA_TABLE(DielectronsVertexing, "AOD", "RTDIELECTRONVTX", //! + reducedpair::Mass, + reducedpair::Pt, reducedpair::Eta, reducedpair::Phi, reducedpair::Sign, + reducedpair::FilterMap, + reducedpair::McDecision, + dilepton_track_index::Pt1, dilepton_track_index::Eta1, dilepton_track_index::Phi1, dilepton_track_index::ITSClusterMap1, dilepton_track_index::ITSChi2NCl1, dilepton_track_index::TPCNClsCR1, dilepton_track_index::TPCNClsFound1, dilepton_track_index::TPCChi2NCl1, dilepton_track_index::DcaXY1, dilepton_track_index::DcaZ1, dilepton_track_index::TPCSignal1, dilepton_track_index::TPCNSigmaEl1, dilepton_track_index::TPCNSigmaPi1, dilepton_track_index::TPCNSigmaPr1, dilepton_track_index::TOFBeta1, dilepton_track_index::TOFNSigmaEl1, dilepton_track_index::TOFNSigmaPi1, dilepton_track_index::TOFNSigmaPr1, + dilepton_track_index::Pt2, dilepton_track_index::Eta2, dilepton_track_index::Phi2, dilepton_track_index::ITSClusterMap2, dilepton_track_index::ITSChi2NCl2, dilepton_track_index::TPCNClsCR2, dilepton_track_index::TPCNClsFound2, dilepton_track_index::TPCChi2NCl2, dilepton_track_index::DcaXY2, dilepton_track_index::DcaZ2, dilepton_track_index::TPCSignal2, dilepton_track_index::TPCNSigmaEl2, dilepton_track_index::TPCNSigmaPi2, dilepton_track_index::TPCNSigmaPr2, dilepton_track_index::TOFBeta2, dilepton_track_index::TOFNSigmaEl2, dilepton_track_index::TOFNSigmaPi2, dilepton_track_index::TOFNSigmaPr2, + reducedpair::Lxy, reducedpair::Tauxy, reducedpair::TauxyErr, reducedpair::Lz, reducedpair::Tauz, reducedpair::TauzErr); + // reducedpair::LxyProj, reducedpair::TauxyProj, reducedpair::TauxyProjErr, reducedpair::LzProj, reducedpair::TauzProj, reducedpair::TauzProjErr); + DECLARE_SOA_TABLE(DimuonsAll, "AOD", "RTDIMUONALL", //! collision::PosX, collision::PosY, collision::PosZ, collision::NumContrib, evsel::Selection, reducedpair::EventSelection, diff --git a/PWGDQ/TableProducer/tableMaker.cxx b/PWGDQ/TableProducer/tableMaker.cxx index 7ebb5876bf5..6d32d5f4436 100644 --- a/PWGDQ/TableProducer/tableMaker.cxx +++ b/PWGDQ/TableProducer/tableMaker.cxx @@ -237,6 +237,8 @@ struct TableMaker { // TODO: filter on TPC dedx used temporarily until electron PID will be improved Filter barrelSelectedTracks = ifnode(fIsRun2.node() == true, aod::track::trackType == uint8_t(aod::track::Run2Track), aod::track::trackType == uint8_t(aod::track::Track)) && o2::aod::track::pt >= fConfigBarrelTrackPtLow && nabs(o2::aod::track::eta) <= fConfigBarrelTrackMaxAbsEta && o2::aod::track::tpcSignal >= configTpcSignal.fConfigMinTpcSignal && o2::aod::track::tpcSignal <= configTpcSignal.fConfigMaxTpcSignal && o2::aod::track::tpcChi2NCl < 4.0f && o2::aod::track::itsChi2NCl < 36.0f; + // Filter barrelSelectedTracks = o2::aod::track::tpcSignal <= configTpcSignal.fConfigMaxTpcSignal && o2::aod::track::tpcChi2NCl < 4.0f && o2::aod::track::itsChi2NCl < 36.0f; + // Filter barrelSelectedTracks = o2::aod::track::tpcChi2NCl < 4.0f && o2::aod::track::itsChi2NCl < 36.0f; Filter muonFilter = o2::aod::fwdtrack::pt >= fConfigMuonPtLow; diff --git a/PWGDQ/TableProducer/tableMakerMC.cxx b/PWGDQ/TableProducer/tableMakerMC.cxx index 1d43d196f9d..3f2de522b8e 100644 --- a/PWGDQ/TableProducer/tableMakerMC.cxx +++ b/PWGDQ/TableProducer/tableMakerMC.cxx @@ -105,6 +105,57 @@ constexpr static uint32_t gkMuonFillMapWithCovAmbi = VarManager::ObjTypes::Muon constexpr static uint32_t gkTrackFillMapWithAmbi = VarManager::ObjTypes::Track | VarManager::ObjTypes::AmbiTrack; constexpr static uint32_t gkMFTFillMap = VarManager::ObjTypes::TrackMFT; +struct TrackQA { + const char* SignalName = "eFromPromptJpsi"; + MCSignal* sig; + + HistogramRegistry registry{"registry"}; + + Filter barrelSelectedTracks = o2::aod::track::pt >= 0.15f && nabs(o2::aod::track::eta) <= 0.9f; + void init(o2::framework::InitContext&) + { + sig = o2::aod::dqmcsignals::GetMCSignal(SignalName); + registry.add("hTrackType", "Track type enum", HistType::kTH1I, {{4, 0, 4}}); + registry.add("hTrackIUPt", "Pt of TrackIU", HistType::kTH1F, {{50, 0, 10}}); + registry.add("hTrackPt", "Pt of Track", HistType::kTH1F, {{50, 0, 10}}); + LOGP(info, "intialized track qa"); + } + + void process(MyEvents::iterator const& event, soa::Filtered const& tracks, aod::McCollisions const& /*mcEvents*/, aod::McParticles_001 const& /*mcTracks*/) { + // if (!event.sel8()) { + // return; + // } + for (auto& track : tracks) { + LOGP(info, "loop over tracks"); + if (track.pt() < 0.15) { + continue; + } + if (abs(track.eta()) > 0.9) { + continue; + } + if (!track.has_mcParticle()) { + continue; + LOGP(info, "Not found MC particle"); + } + auto mctrack = track.template mcParticle_as(); + LOGP(info, "found MC particle: {}", mctrack.pdgCode()); + if (sig->CheckSignal(true, mctrack)) { + LOGP(info, "found non-prompt Jpsi signal"); + if (track.trackType() == uint8_t(aod::track::TrackIU)) { + registry.fill(HIST("hTrackType"), 0); + registry.fill(HIST("hTrackIUPt"), track.pt()); + } else if (track.trackType() == uint8_t(aod::track::Track)) { + registry.fill(HIST("hTrackType"), 1); + registry.fill(HIST("hTrackPt"), track.pt()); + } else if (track.trackType() == uint8_t(aod::track::StrangeTrack)) { + registry.fill(HIST("hTrackType"), 2); + } + registry.fill(HIST("hTrackType"), 3); + } + } + } +}; + struct TableMakerMC { Produces event; @@ -179,6 +230,9 @@ struct TableMakerMC { void init(o2::framework::InitContext& context) { + if (context.mOptions.get("processDummy")) { + return; + } fCCDB->setURL(fConfigCcdbUrl); fCCDB->setCaching(true); fCCDB->setLocalObjectValidityChecking(); @@ -1759,6 +1813,10 @@ struct TableMakerMC { (reinterpret_cast(fStatsList->At(0)))->Fill(0.0, static_cast(kNaliases)); } + void processDummy(MyEvents&) { + // do nothing + } + PROCESS_SWITCH(TableMakerMC, processFull, "Produce both barrel and muon skims", false); PROCESS_SWITCH(TableMakerMC, processFullWithCov, "Produce both barrel and muon skims, w/ track and fwdtrack cov tables", false); PROCESS_SWITCH(TableMakerMC, processBarrelOnly, "Produce barrel skims", false); @@ -1780,6 +1838,7 @@ struct TableMakerMC { PROCESS_SWITCH(TableMakerMC, processAmbiguousMuonOnly, "Build muon-only DQ skimmed data model with QA plots for ambiguous muons", false); PROCESS_SWITCH(TableMakerMC, processAmbiguousMuonOnlyWithCov, "Build muon-only with cov DQ skimmed data model with QA plots for ambiguous muons", false); PROCESS_SWITCH(TableMakerMC, processAmbiguousBarrelOnly, "Build barrel-only DQ skimmed data model with QA plots for ambiguous tracks", false); + PROCESS_SWITCH(TableMakerMC, processDummy, "DummyFunction", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) @@ -1787,5 +1846,6 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) // TODO: For now TableMakerMC works just for PbPb (cent table is present) // Implement workflow arguments for pp/PbPb and possibly merge the task with tableMaker.cxx return WorkflowSpec{ + adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; } diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index f9be380c416..8b4e42a3a49 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -88,6 +88,12 @@ using MyBarrelTracksWithCov = soa::Join; +using MyBarrelTracksTunedWithCov = soa::Join; using MyMuons = soa::Join; using MyMuonsWithCov = soa::Join; using MyMuonsRealignWithCov = soa::Join; @@ -108,6 +114,7 @@ constexpr static uint32_t gkEventFillMapWithMultsRapidityGapFilter = VarManager: // constexpr static uint32_t gkEventMCFillMap = VarManager::ObjTypes::CollisionMC; // constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackSelection | VarManager::ObjTypes::TrackPID; constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackSelection | VarManager::ObjTypes::TrackCov | VarManager::ObjTypes::TrackPID; +constexpr static uint32_t gkTrackFillMapTunedWithCov = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackSelection | VarManager::ObjTypes::TrackCov | VarManager::ObjTypes::TrackPID | VarManager::ObjTypes::MCTPCtuneOnData; // constexpr static uint32_t gkTrackFillMapWithDalitzBits = gkTrackFillMap | VarManager::ObjTypes::DalitzBits; // constexpr static uint32_t gkMuonFillMap = VarManager::ObjTypes::Muon; constexpr static uint32_t gkMuonFillMapWithCov = VarManager::ObjTypes::Muon | VarManager::ObjTypes::MuonCov; @@ -260,7 +267,7 @@ struct TableMakerMC { { // Check whether barrel or muon are enabled bool isProcessBCenabled = context.mOptions.get("processPP"); - bool isBarrelEnabled = (context.mOptions.get("processPP") || context.mOptions.get("processPPBarrelOnly") || context.mOptions.get("processPbPbBarrelOnly") || context.mOptions.get("processPbPbWithFilterBarrelOnly")); + bool isBarrelEnabled = (context.mOptions.get("processPP") || context.mOptions.get("processPPBarrelOnly") || context.mOptions.get("processPPBarrelOnlyMcTuned") || context.mOptions.get("processPbPbBarrelOnly") || context.mOptions.get("processPbPbWithFilterBarrelOnly")); bool isMuonEnabled = (context.mOptions.get("processPP") || context.mOptions.get("processPPMuonOnlyBasic") || context.mOptions.get("processPPMuonOnly") || context.mOptions.get("processPPRealignedMuonOnly") || context.mOptions.get("processPbPbMuonOnly") || context.mOptions.get("processPbPbRealignedMuonOnly")) || context.mOptions.get("processPPMuonRefit"); // Make sure at least one process function is enabled if (!(isProcessBCenabled || isBarrelEnabled || isMuonEnabled)) { @@ -1455,6 +1462,13 @@ struct TableMakerMC { fullSkimming(collisions, bcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, mcCollisions, mcParticles, nullptr); } + void processPPBarrelOnlyMcTuned(MyEventsWithMults const& collisions, aod::BCsWithTimestamps const& bcs, + MyBarrelTracksTunedWithCov const& tracksBarrel, aod::TrackAssoc const& trackAssocs, + aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles) + { + fullSkimming(collisions, bcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, mcCollisions, mcParticles, nullptr); + } + void processPPMuonOnlyBasic(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, MyMuonsWithCov const& tracksMuon, MFTTrackLabeled const& mftTracks, aod::FwdTrackAssoc const& fwdTrackAssocs, aod::MFTTrackAssoc const& mftAssocs, @@ -1538,6 +1552,7 @@ struct TableMakerMC { PROCESS_SWITCH(TableMakerMC, processPP, "Produce both barrel and muon skims, pp settings", false); PROCESS_SWITCH(TableMakerMC, processPPBarrelOnly, "Produce only barrel skims, pp settings ", false); + PROCESS_SWITCH(TableMakerMC, processPPBarrelOnlyMcTuned, "Produce only barrel skims, pp settings, with MC-tuned tracking efficiency", false); PROCESS_SWITCH(TableMakerMC, processPPMuonOnlyBasic, "Produce only muon skims, pp settings, no multiplicity", false); PROCESS_SWITCH(TableMakerMC, processPPMuonOnly, "Produce only muon skims, pp settings", false); PROCESS_SWITCH(TableMakerMC, processPPMuonRefit, "Produce only muon skims, pp settings", false); diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 475bc6d6f8e..22d5635c8eb 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -1289,6 +1289,7 @@ struct AnalysisSameEventPairing { Produces dielectronsExtraList; Produces dielectronInfoList; Produces dielectronAllList; + Produces dielectronVtx; Produces dimuonsExtraList; Produces dimuonAllList; Produces dileptonMiniTreeGen; @@ -1754,7 +1755,8 @@ struct AnalysisSameEventPairing { dielectronInfoList.reserve(1); dileptonInfoList.reserve(1); if (fConfigOptions.flatTables.value) { - dielectronAllList.reserve(1); + // dielectronAllList.reserve(1); + dielectronVtx.reserve(1); dimuonAllList.reserve(1); } if (useMiniTree.fConfigMiniTree) { @@ -1844,16 +1846,20 @@ struct AnalysisSameEventPairing { dielectronsExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingLzProjected], VarManager::fgValues[VarManager::kVertexingLxyProjected]); if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrackCollInfo) > 0) { if (fConfigOptions.flatTables.value && t1.has_reducedMCTrack() && t2.has_reducedMCTrack()) { - dielectronAllList(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, mcDecision, + // dielectronAllList(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, mcDecision, + // t1.pt(), t1.eta(), t1.phi(), t1.itsClusterMap(), t1.itsChi2NCl(), t1.tpcNClsCrossedRows(), t1.tpcNClsFound(), t1.tpcChi2NCl(), t1.dcaXY(), t1.dcaZ(), t1.tpcSignal(), t1.tpcNSigmaEl(), t1.tpcNSigmaPi(), t1.tpcNSigmaPr(), t1.beta(), t1.tofNSigmaEl(), t1.tofNSigmaPi(), t1.tofNSigmaPr(), + // t2.pt(), t2.eta(), t2.phi(), t2.itsClusterMap(), t2.itsChi2NCl(), t2.tpcNClsCrossedRows(), t2.tpcNClsFound(), t2.tpcChi2NCl(), t2.dcaXY(), t2.dcaZ(), t2.tpcSignal(), t2.tpcNSigmaEl(), t2.tpcNSigmaPi(), t2.tpcNSigmaPr(), t2.beta(), t2.tofNSigmaEl(), t2.tofNSigmaPi(), t2.tofNSigmaPr(), + // VarManager::fgValues[VarManager::kKFTrack0DCAxyz], VarManager::fgValues[VarManager::kKFTrack1DCAxyz], VarManager::fgValues[VarManager::kKFDCAxyzBetweenProngs], VarManager::fgValues[VarManager::kKFTrack0DCAxy], VarManager::fgValues[VarManager::kKFTrack1DCAxy], VarManager::fgValues[VarManager::kKFDCAxyBetweenProngs], + // VarManager::fgValues[VarManager::kKFTrack0DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack0DeviationxyFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationxyFromPV], + // VarManager::fgValues[VarManager::kKFMass], VarManager::fgValues[VarManager::kKFChi2OverNDFGeo], VarManager::fgValues[VarManager::kVertexingLxyz], VarManager::fgValues[VarManager::kVertexingLxyzOverErr], VarManager::fgValues[VarManager::kVertexingLxy], VarManager::fgValues[VarManager::kVertexingLxyOverErr], VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxyErr], VarManager::fgValues[VarManager::kKFCosPA], VarManager::fgValues[VarManager::kKFJpsiDCAxyz], VarManager::fgValues[VarManager::kKFJpsiDCAxy], + // VarManager::fgValues[VarManager::kKFPairDeviationFromPV], VarManager::fgValues[VarManager::kKFPairDeviationxyFromPV], + // VarManager::fgValues[VarManager::kKFMassGeoTop], VarManager::fgValues[VarManager::kKFChi2OverNDFGeoTop], + // VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingTauxyProjected], + // VarManager::fgValues[VarManager::kVertexingLzProjected], VarManager::fgValues[VarManager::kVertexingLxyProjected]); + dielectronVtx(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, mcDecision, t1.pt(), t1.eta(), t1.phi(), t1.itsClusterMap(), t1.itsChi2NCl(), t1.tpcNClsCrossedRows(), t1.tpcNClsFound(), t1.tpcChi2NCl(), t1.dcaXY(), t1.dcaZ(), t1.tpcSignal(), t1.tpcNSigmaEl(), t1.tpcNSigmaPi(), t1.tpcNSigmaPr(), t1.beta(), t1.tofNSigmaEl(), t1.tofNSigmaPi(), t1.tofNSigmaPr(), t2.pt(), t2.eta(), t2.phi(), t2.itsClusterMap(), t2.itsChi2NCl(), t2.tpcNClsCrossedRows(), t2.tpcNClsFound(), t2.tpcChi2NCl(), t2.dcaXY(), t2.dcaZ(), t2.tpcSignal(), t2.tpcNSigmaEl(), t2.tpcNSigmaPi(), t2.tpcNSigmaPr(), t2.beta(), t2.tofNSigmaEl(), t2.tofNSigmaPi(), t2.tofNSigmaPr(), - VarManager::fgValues[VarManager::kKFTrack0DCAxyz], VarManager::fgValues[VarManager::kKFTrack1DCAxyz], VarManager::fgValues[VarManager::kKFDCAxyzBetweenProngs], VarManager::fgValues[VarManager::kKFTrack0DCAxy], VarManager::fgValues[VarManager::kKFTrack1DCAxy], VarManager::fgValues[VarManager::kKFDCAxyBetweenProngs], - VarManager::fgValues[VarManager::kKFTrack0DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack0DeviationxyFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationxyFromPV], - VarManager::fgValues[VarManager::kKFMass], VarManager::fgValues[VarManager::kKFChi2OverNDFGeo], VarManager::fgValues[VarManager::kVertexingLxyz], VarManager::fgValues[VarManager::kVertexingLxyzOverErr], VarManager::fgValues[VarManager::kVertexingLxy], VarManager::fgValues[VarManager::kVertexingLxyOverErr], VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxyErr], VarManager::fgValues[VarManager::kKFCosPA], VarManager::fgValues[VarManager::kKFJpsiDCAxyz], VarManager::fgValues[VarManager::kKFJpsiDCAxy], - VarManager::fgValues[VarManager::kKFPairDeviationFromPV], VarManager::fgValues[VarManager::kKFPairDeviationxyFromPV], - VarManager::fgValues[VarManager::kKFMassGeoTop], VarManager::fgValues[VarManager::kKFChi2OverNDFGeoTop], - VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingTauxyProjected], - VarManager::fgValues[VarManager::kVertexingLzProjected], VarManager::fgValues[VarManager::kVertexingLxyProjected]); + VarManager::fgValues[VarManager::kVertexingLxy], VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxy]/VarManager::fgValues[VarManager::kVertexingTauxyErr], VarManager::fgValues[VarManager::kVertexingLz], VarManager::fgValues[VarManager::kVertexingTauz],VarManager::fgValues[VarManager::kVertexingTauz]/VarManager::fgValues[VarManager::kVertexingTauzErr]); } } } @@ -2268,7 +2274,8 @@ struct AnalysisSameEventPairing { continue; } if (sig->CheckSignal(true, t1_raw, t2_raw)) { - VarManager::FillPairMC(t1, t2); // NOTE: This feature will only work for muons + // mcDecision |= (static_cast(1) << isig); + VarManager::FillPairMC(t1, t2); // NOTE: This feature will only work for electrons fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), VarManager::fgValues); } } @@ -4150,6 +4157,477 @@ struct AnalysisDileptonTrack { PROCESS_SWITCH(AnalysisDileptonTrack, processDummy, "Dummy function", true); }; +struct AnalysisDileptonTrackTrack{ + OutputObj fOutputList{"output"}; + + Configurable fConfigTrackCut1{"cfgTrackCut1", "pionPIDCut1", "track1 cut"}; // used for select the tracks from SelectedTracks + Configurable fConfigTrackCut2{"cfgTrackCut2", "pionPIDCut2", "track2 cut"}; // used for select the tracks from SelectedTracks + Configurable fConfigDileptonCut{"cfgDileptonCut", "pairJpsi2", "Dilepton cut"}; + Configurable fConfigQuadrupletCuts{"cfgQuadrupletCuts", "pairX3872Cut1", "Comma separated list of Dilepton-Track-Track cut"}; + Configurable fConfigMCRecSignals{"cfgBarrelMCRecSignals", "", "Comma separated list of MC signals (reconstructed)"}; + Configurable fConfigMCGenSignals{"cfgBarrelMCGenSignals", "", "Comma separated list of MC signals (generated)"}; + Configurable fConfigDileptonMCRecSignal{"cfgDileptonMCRecSignal", "", "Comma separated list of MC signals (reconstructed)"}; + Configurable fConfigAddDileptonHistogram{"cfgAddDileptonHistogram", "barrel", "Comma separated list of histograms"}; + Configurable fConfigAddQuadrupletHistogram{"cfgAddQuadrupletHistogram", "xtojpsipipi", "Comma separated list of histograms"}; + + Configurable fConfigUseKFVertexing{"cfgUseKFVertexing", false, "Use KF Particle for secondary vertex reconstruction (DCAFitter is used by default)"}; + Configurable fConfigSetupFourProngFitter{"cfgSetupFourProngFitter", false, "Use DCA for secondary vertex reconstruction (DCAFitter is used by default)"}; + Configurable fConfigGRPmagPath{"cfgGrpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable fConfigUseRemoteField{"cfgUseRemoteField", false, "Chose whether to fetch the magnetic field from ccdb or set it manually"}; + Configurable fConfigMagField{"cfgMagField", 5.0f, "Manually set magnetic field"}; + Configurable fConfigDileptonLxyCut{"cfgDileptonLxyCut", -999.f, "Dilepton Lxy cut for secondary vertex"}; + + Produces DileptonTrackTrackTable; + + int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. + // uint32_t fTrackCutBitMap; // track cut bit mask to be used in the selection of tracks associated with dileptons + // cut name setting + TString fTrackCutName1; + TString fTrackCutName2; + bool fIsSameTrackCut = false; + AnalysisCompositeCut fDileptonCut; + std::vector fQuadrupletCutNames; + std::vector fQuadrupletCuts; + + Service fCCDB; + + Filter eventFilter = aod::dqanalysisflags::isEventSelected > static_cast(0); + Filter dileptonFilter = aod::reducedpair::sign == 0 && aod::reducedpair::lxy > fConfigDileptonLxyCut; + Filter filterBarrelTrackSelected = aod::dqanalysisflags::isBarrelSelected > static_cast(0); + + constexpr static uint32_t fgDileptonFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::Pair; // fill map + + float* fValuesQuadruplet; + HistogramManager* fHistMan; + + std::vector fRecMCSignals; + std::vector fGenMCSignals; + + void init(o2::framework::InitContext& context) + { + // return; + bool isX3872 = context.mOptions.get("processX3872"); + bool isPsi2S = context.mOptions.get("processPsi2S"); + bool isAmbi = context.mOptions.get("processPsi2SNoAmbi"); + bool isMCGen = context.mOptions.get("processMCGen") || context.mOptions.get("processMCGenWithEventSelection"); + bool isDummy = context.mOptions.get("processDummy"); + + if (isDummy) { + if (isPsi2S || isX3872 || isMCGen) { + LOG(fatal) << "Dummy function is enabled even if there are normal process functions running! Fix your config!" << endl; + return; + } else { + LOG(info) << "Dummy function is enabled. Skipping the rest of the init function" << endl; + return; + } + + } + + fCurrentRun = 0; + fValuesQuadruplet = new float[VarManager::kNVars]; + VarManager::SetDefaultVarNames(); + fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + fHistMan->SetUseDefaultVariableNames(true); + fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + + TString sigNamesStr = fConfigMCRecSignals.value; + std::unique_ptr objRecSigArray(sigNamesStr.Tokenize(",")); + if (!sigNamesStr.IsNull()) { + for (int isig = 0; isig < objRecSigArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objRecSigArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() != 4) { + LOG(fatal) << "Signal at reconstructed level requested (" << sig->GetName() << ") " << "does not have 4 prongs! Fix it"; + } + fRecMCSignals.push_back(sig); + } else { + LOG(fatal) << "Signal at reconstructed level requested (" << objRecSigArray->At(isig)->GetName() << ") " << "could not be retrieved from the library! -> skipped"; + } + } + } + + TString sigGenNamesStr = fConfigMCGenSignals.value; + std::unique_ptr objGenSigArray(sigGenNamesStr.Tokenize(",")); + for (int isig = 0; isig < objGenSigArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objGenSigArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() ==1) { // NOTE: 1-prong signals required + fGenMCSignals.push_back(sig); + } + } + } + + // ToDo: get cut names from above struct + fTrackCutName1 = fConfigTrackCut1.value; + fTrackCutName2 = fConfigTrackCut2.value; + if (fTrackCutName1 == fTrackCutName2) { + fIsSameTrackCut = true; + } + TString configDileptonCutNamesStr = fConfigDileptonCut.value; + fDileptonCut = *dqcuts::GetCompositeCut(configDileptonCutNamesStr.Data()); + TString configQuadruletCutNamesStr = fConfigQuadrupletCuts.value; + std::unique_ptr objArray(configQuadruletCutNamesStr.Tokenize(",")); + for (Int_t icut = 0; icut < objArray->GetEntries(); ++icut) { + TString cutName = objArray->At(icut)->GetName(); + fQuadrupletCutNames.push_back(cutName); + fQuadrupletCuts.push_back(*dqcuts::GetCompositeCut(cutName.Data())); + } + + if (isPsi2S || isX3872 || isAmbi) { + DefineHistograms(fHistMan, Form("Pairs_%s", configDileptonCutNamesStr.Data()), fConfigAddDileptonHistogram.value.data()); + if (!configQuadruletCutNamesStr.IsNull()) { + for (std::size_t icut = 0; icut < fQuadrupletCutNames.size(); ++icut) { + if (fIsSameTrackCut) { + DefineHistograms(fHistMan, Form("QuadrupletSEPM_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + } else { + DefineHistograms(fHistMan, Form("QuadrupletSEPM_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + DefineHistograms(fHistMan, Form("QuadrupletSEMP_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + } + DefineHistograms(fHistMan, Form("QuadrupletSEPP_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + DefineHistograms(fHistMan, Form("QuadrupletSEMM_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + for (auto& sig : fRecMCSignals) { + DefineHistograms(fHistMan, Form("QuadrupletMCMatched_%s_%s", fQuadrupletCutNames[icut].Data(), sig->GetName()), fConfigAddQuadrupletHistogram.value.data()); + } + } + } + } + + if (isMCGen) { + for (auto& sig : fGenMCSignals) { + DefineHistograms(fHistMan, Form("MCTruthGenQuad_%s", sig->GetName()), ""); + DefineHistograms(fHistMan, Form("MCTruthGenQuadSel_%s", sig->GetName()), ""); + } + // DefineHistograms(fHistMan, "MCTruthGenQuadAccepted", ""); + } + + VarManager::SetUseVars(fHistMan->GetUsedVars()); + fOutputList.setObject(fHistMan->GetMainHistogramList()); + } + + // init parameters from CCDB + void initParamsFromCCDB(uint64_t timestamp) + { + if (fConfigUseRemoteField.value) { + o2::parameters::GRPMagField* grpmag = fCCDB->getForTimeStamp(fConfigGRPmagPath.value, timestamp); + float magField = 0.0; + if (grpmag != nullptr) { + magField = grpmag->getNominalL3Field(); + } else { + LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", timestamp); + } + if (fConfigUseKFVertexing.value) { + VarManager::SetupTwoProngKFParticle(magField); + VarManager::SetupFourProngKFParticle(magField); + } else if (fConfigSetupFourProngFitter.value) { + VarManager::SetupTwoProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + VarManager::SetupFourProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } + } else { + if (fConfigUseKFVertexing.value) { + VarManager::SetupTwoProngKFParticle(fConfigMagField.value); + VarManager::SetupFourProngKFParticle(fConfigMagField.value); + } else if (fConfigSetupFourProngFitter.value) { + LOGP(info, "Setting up DCA fitter for two and four prong candidates"); + VarManager::SetupTwoProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + VarManager::SetupFourProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } + } + } + + // Template function to run pair - hadron combinations + template + void runDileptonTrackTrack(TEvent const& event, TTrackAssocs const& assocs, TTracks const& tracks, TDileptons const& dileptons, ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& /*mcTracks*/) + { + VarManager::ResetValues(0, VarManager::kNVars, fValuesQuadruplet); + VarManager::FillEvent(event, fValuesQuadruplet); + VarManager::FillEvent(event.reducedMCevent(), fValuesQuadruplet); + + uint32_t mcDecision = static_cast(0); + size_t isig = 0; + + for (auto dilepton : dileptons) { + // get full track info of tracks based on the index + + int indexLepton1 = dilepton.index0Id(); + int indexLepton2 = dilepton.index1Id(); + auto lepton1 = tracks.rawIteratorAt(dilepton.index0Id()); + auto lepton2 = tracks.rawIteratorAt(dilepton.index1Id()); + auto lepton1MC = lepton1.reducedMCTrack(); + auto lepton2MC = lepton2.reducedMCTrack(); + // Check that the dilepton has zero charge + if (dilepton.sign() != 0) { + continue; + } + VarManager::FillTrack(dilepton, fValuesQuadruplet); + + bool isAmbiguousLepton = (dilepton.filterMap_raw() & (static_cast(1) << 28)) || + (dilepton.filterMap_raw() & (static_cast(1) << 29)) || + (dilepton.filterMap_raw() & (static_cast(1) << 30)) || + (dilepton.filterMap_raw() & (static_cast(1) << 31)); + // if (isAmbi && isAmbiguousLepton) + // continue; // skip ambiguous dileptons + if constexpr ((TTrackFillMap & VarManager::ObjTypes::AmbiTrack) > 0) { + if (isAmbiguousLepton) { + // LOGP(info, "Dilepton {} is ambiguous, skipping", dilepton.index0Id()); + continue; // skip ambiguous dileptons + } + } + + + // apply the dilepton cut + if (!fDileptonCut.IsSelected(fValuesQuadruplet)) + continue; + + fHistMan->FillHistClass(Form("Pairs_%s", fDileptonCut.GetName()), fValuesQuadruplet); + + // loop over hadrons pairs + for (auto& [a1, a2] : o2::soa::combinations(assocs, assocs)) { + uint32_t trackSelection = 0; + if (fIsSameTrackCut) { + trackSelection = ((a1.isBarrelSelected_raw() & (static_cast(1) << 1)) || (a2.isBarrelSelected_raw() & (static_cast(1) << 1))); + } else { + trackSelection = ((a1.isBarrelSelected_raw() & (static_cast(1) << 1)) && (a2.isBarrelSelected_raw() & (static_cast(1) << 2))); + } + // LOGP(info, "trackSelection: {}, a1: {}, a2: {}", trackSelection, a1.isBarrelSelected_raw(), a2.isBarrelSelected_raw()); + if (!trackSelection) { + continue; + } + + // get the track from this association + auto track1 = a1.template reducedtrack_as(); + auto track2 = a2.template reducedtrack_as(); + // avoid self combinations + if (track1.globalIndex() == indexLepton1 || track1.globalIndex() == indexLepton2 || track2.globalIndex() == indexLepton1 || track2.globalIndex() == indexLepton2) { + continue; + } + + // if (isAmbi) { + // check that the tracks are not ambiguous + if constexpr ((TTrackFillMap & VarManager::ObjTypes::AmbiTrack) > 0) { + if ((track1.barrelAmbiguityInBunch() > 1) || (track2.barrelAmbiguityInBunch() > 1) || + (track1.barrelAmbiguityOutOfBunch() > 1) || (track2.barrelAmbiguityOutOfBunch() > 1)) { + // LOGP(info, "Track {} or {} is ambiguous, skipping", track1.globalIndex(), track2.globalIndex()); + continue; + } + } + // } + + // fill variables + VarManager::FillDileptonTrackTrack(dilepton, track1, track2, fValuesQuadruplet); + if (fConfigSetupFourProngFitter || fConfigUseKFVertexing) { + // LOGP(info, "Using KF or DCA fitter for secondary vertexing"); + VarManager::FillDileptonTrackTrackVertexing(event, lepton1, lepton2, track1, track2, fValuesQuadruplet); + } + + auto track1MC = track1.reducedMCTrack(); + auto track2MC = track2.reducedMCTrack(); + mcDecision = 0; + isig = 0; + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + if ((*sig)->CheckSignal(true,lepton1MC, lepton2MC, track1MC, track2MC)) { + mcDecision |= (static_cast(1) << isig); + } + } + + int iCut = 0; + uint32_t CutDecision = 0; + for (auto cutname = fQuadrupletCutNames.begin(); cutname != fQuadrupletCutNames.end(); cutname++, iCut++) { + // apply dilepton-track-track cut + if (fQuadrupletCuts[iCut].IsSelected(fValuesQuadruplet)) { + CutDecision |= (1 << iCut); + if (fIsSameTrackCut) { + if (track1.sign() * track2.sign() < 0) { + fHistMan->FillHistClass(Form("QuadrupletSEPM_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + for (uint32_t isig = 0; isig < fRecMCSignals.size(); isig++) { + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(Form("QuadrupletMCMatched_%s_%s", fQuadrupletCutNames[iCut].Data(), fRecMCSignals[isig]->GetName()), fValuesQuadruplet); + } + } + } + } else { + if ((track1.sign() < 0) && (track2.sign() > 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEMP_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + } else if ((track1.sign() > 0) && (track2.sign() < 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEPM_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + } + } + if ((track1.sign() > 0) && (track2.sign() > 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEPP_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + } else if ((track1.sign() < 0) && (track2.sign() < 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEMM_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + } + } + } // end loop over quadruplet cuts + + // fill table + if (!CutDecision) + continue; + if (!mcDecision) + continue; + DileptonTrackTrackTable(fValuesQuadruplet[VarManager::kQuadMass], fValuesQuadruplet[VarManager::kQuadPt], fValuesQuadruplet[VarManager::kQuadEta], fValuesQuadruplet[VarManager::kQuadPhi], fValuesQuadruplet[VarManager::kRap], + fValuesQuadruplet[VarManager::kQ], fValuesQuadruplet[VarManager::kDeltaR1], fValuesQuadruplet[VarManager::kDeltaR2], fValuesQuadruplet[VarManager::kDeltaR], + dilepton.mass(), dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.sign(), + lepton1.tpcNSigmaEl(), lepton1.tpcNSigmaPi(), lepton1.tpcNSigmaPr(), lepton1.tpcNClsFound(), + lepton2.tpcNSigmaEl(), lepton2.tpcNSigmaPi(), lepton2.tpcNSigmaPr(), lepton2.tpcNClsFound(), + fValuesQuadruplet[VarManager::kDitrackMass], fValuesQuadruplet[VarManager::kDitrackPt], track1.pt(), track2.pt(), track1.eta(), track2.eta(), track1.phi(), track2.phi(), track1.sign(), track2.sign(), track1.tpcNSigmaPi(), track2.tpcNSigmaPi(), track1.tpcNSigmaKa(), track2.tpcNSigmaKa(), track1.tpcNSigmaPr(), track1.tpcNSigmaPr(), track1.tpcNClsFound(), track2.tpcNClsFound(), + fValuesQuadruplet[VarManager::kKFMass], fValuesQuadruplet[VarManager::kVertexingProcCode], fValuesQuadruplet[VarManager::kVertexingChi2PCA], fValuesQuadruplet[VarManager::kCosPointingAngle], fValuesQuadruplet[VarManager::kKFDCAxyzBetweenProngs], fValuesQuadruplet[VarManager::kKFChi2OverNDFGeo], + fValuesQuadruplet[VarManager::kVertexingLz], fValuesQuadruplet[VarManager::kVertexingLxy], fValuesQuadruplet[VarManager::kVertexingLxyz], fValuesQuadruplet[VarManager::kVertexingTauz], fValuesQuadruplet[VarManager::kVertexingTauxy], fValuesQuadruplet[VarManager::kVertexingLzErr], fValuesQuadruplet[VarManager::kVertexingLxyzErr], + fValuesQuadruplet[VarManager::kVertexingTauzErr], fValuesQuadruplet[VarManager::kVertexingLzProjected], fValuesQuadruplet[VarManager::kVertexingLxyProjected], fValuesQuadruplet[VarManager::kVertexingLxyzProjected], fValuesQuadruplet[VarManager::kVertexingTauzProjected], fValuesQuadruplet[VarManager::kVertexingTauxyProjected]); + } // end loop over associations + } // end loop over dileptons + } + + Preslice trackAssocsPerCollision = aod::reducedtrack_association::reducedeventId; + Preslice dielectronsPerCollision = aod::reducedpair::reducedeventId; + Preslice ditracksPerCollision = aod::reducedpair::reducedeventId; + + void processX3872(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracksWithCov const& tracks, soa::Filtered const& dileptons, + ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks) + { + // set up KF or DCAfitter + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp()); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack(event, groupedBarrelAssocs, tracks, groupedDielectrons, mcEvents, mcTracks); + } + } + + void processPsi2S(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracksWithCov const& tracks, soa::Filtered const& dileptons, + ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks) + { + // set up KF or DCAfitter + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp()); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack(event, groupedBarrelAssocs, tracks, groupedDielectrons, mcEvents, mcTracks); + } + } + + void processPsi2SNoAmbi(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracksWithCovWithAmbiguities const& tracks, soa::Filtered const& dileptons, + ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks) + { + // set up KF or DCAfitter + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp()); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack(event, groupedBarrelAssocs, tracks, groupedDielectrons, mcEvents, mcTracks); + } + } + + void processMCGen(ReducedMCTracks const& mcTracks) + { + // loop over mc stack and fill histograms for pure MC truth signals + // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event + // auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex()); + for (auto& mctrack : mcTracks) { + VarManager::FillTrackMC(mcTracks, mctrack); + // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. + // NOTE: However, the working model is that the decisions on MC signals are precomputed during skimming and are stored in the mcReducedFlags member. + // TODO: Use the mcReducedFlags to select signals + for (auto& sig : fGenMCSignals) { + if (sig->CheckSignal(true, mctrack)) { + int daughterIdFirst = mctrack.daughtersIds()[0]; + int daughterIdEnd = mctrack.daughtersIds()[1]; + int Ndaughters = daughterIdEnd - daughterIdFirst + 1; + if (Ndaughters ==3) { + auto dilepton = mcTracks.rawIteratorAt(daughterIdFirst); + auto track1 = mcTracks.rawIteratorAt(daughterIdFirst + 1); + auto track2 = mcTracks.rawIteratorAt(daughterIdFirst + 2); + VarManager::FillQuadMC(dilepton, track1, track2); + } + fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues); + } + } + } + } + + PresliceUnsorted perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + + void processMCGenWithEventSelection(soa::Filtered const& events, + ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& mcTracks) + { + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + if (!event.has_reducedMCevent()) { + continue; + } + + auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.reducedMCeventId()); + groupedMCTracks.bindInternalIndicesTo(&mcTracks); + for (auto& track : groupedMCTracks) { + + VarManager::FillTrackMC(mcTracks, track); + + auto track_raw = groupedMCTracks.rawIteratorAt(track.globalIndex()); + for (auto& sig : fGenMCSignals) { + if (sig->CheckSignal(true, track_raw)) { + int daughterIdFirst = track_raw.daughtersIds()[0]; + int daughterIdEnd = track_raw.daughtersIds()[1]; + int Ndaughters = daughterIdEnd - daughterIdFirst + 1; + if (Ndaughters == 3) { + auto dilepton = groupedMCTracks.rawIteratorAt(daughterIdFirst); + auto track1 = groupedMCTracks.rawIteratorAt(daughterIdFirst + 1); + auto track2 = groupedMCTracks.rawIteratorAt(daughterIdFirst + 2); + VarManager::FillQuadMC(dilepton, track1, track2); + } + fHistMan->FillHistClass(Form("MCTruthGenQuad_%s", sig->GetName()), VarManager::fgValues); + } + } + } + } // end loop over events + } + + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processPsi2S, "Run psi(2S) -> e+ e- + pi+ pi- pairing, using skimmed data", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processPsi2SNoAmbi, "Run psi(2S) -> e+ e- + pi+ pi- pairing, using skimmed data without ambiguities", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processX3872, "Run X(3872) -> J/psi + pi+ pi- pairing, using skimmed data", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processMCGen, "Loop over MC particle stack and fill generator level histograms", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processMCGenWithEventSelection, "Loop over MC particle stack and fill generator level histograms with event selection", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processDummy, "Dummy function", false); +}; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ @@ -4225,6 +4703,18 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_track"); } + if (classStr.Contains("MCTruthGenAccepted")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_pair"); + } + + if (classStr.Contains("MCTruthSel")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", "mctruth_track"); + } + + if (classStr.Contains("MCTruthGenQuad")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_quad"); + } + if (classStr.Contains("DileptonsSelected")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair", "barrel,vertexing"); } diff --git a/PWGDQ/Tasks/tableReader.cxx b/PWGDQ/Tasks/tableReader.cxx index 2b35db2611c..e5fa216fc7a 100644 --- a/PWGDQ/Tasks/tableReader.cxx +++ b/PWGDQ/Tasks/tableReader.cxx @@ -2135,7 +2135,7 @@ struct AnalysisDileptonTrackTrack { fOutputList.setObject(fHistMan->GetMainHistogramList()); } // Template function to run pair - track - track combinations - template + template void runDileptonTrackTrack(TEvent const& event, TTracks const& tracks, soa::Filtered const& dileptons) { VarManager::ResetValues(0, VarManager::kNVars, fValuesQuadruplet); @@ -2144,12 +2144,25 @@ struct AnalysisDileptonTrackTrack { // LOGF(info, "Number of dileptons: %d", dileptons.size()); // set up KF or DCAfitter - if (fConfigUseDCAVertexing) { - VarManager::SetupTwoProngDCAFitter(5.0f, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables - // VarManager::SetupThreeProngDCAFitter(5.0f, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // if (fConfigUseDCAVertexing) { + // VarManager::SetupTwoProngDCAFitter(5.0f, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // // VarManager::SetupThreeProngDCAFitter(5.0f, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // VarManager::SetupFourProngDCAFitter(5.0f, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // } else if (fConfigUseKFVertexing) { + // VarManager::SetupFourProngKFParticle(5.0f); + // } + if (NProng == 3) { + if (fConfigUseKFVertexing) { + VarManager::SetupThreeProngKFParticle(5.0f); + } else { + VarManager::SetupThreeProngDCAFitter(5.0f, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } + } else if (NProng == 4) { + if (fConfigUseKFVertexing) { + VarManager::SetupFourProngKFParticle(5.0f); + } else { VarManager::SetupFourProngDCAFitter(5.0f, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables - } else if (fConfigUseKFVertexing) { - VarManager::SetupFourProngKFParticle(5.0f); + } } int indexOffset = -999; @@ -2204,7 +2217,7 @@ struct AnalysisDileptonTrackTrack { // fill variables VarManager::FillDileptonTrackTrack(dilepton, t1, t2, fValuesQuadruplet); // reconstruct the secondary vertex - if (fConfigUseDCAVertexing || fConfigUseKFVertexing) { + if constexpr (NProng > 0) { VarManager::FillDileptonTrackTrackVertexing(event, lepton1, lepton2, t1, t2, fValuesQuadruplet); } @@ -2251,12 +2264,17 @@ struct AnalysisDileptonTrackTrack { void processChicToJpsiPiPi(soa::Filtered::iterator const& event, MyBarrelTracksSelectedWithCov const& tracks, soa::Filtered const& dileptons) { - runDileptonTrackTrack(event, tracks, dileptons); + runDileptonTrackTrack<4, VarManager::kXtoJpsiPiPi, gkEventFillMapWithCov, gkTrackFillMapWithCov>(event, tracks, dileptons); } void processPsi2SToJpsiPiPi(soa::Filtered::iterator const& event, MyBarrelTracksSelectedWithCov const& tracks, soa::Filtered const& dileptons) { - runDileptonTrackTrack(event, tracks, dileptons); + runDileptonTrackTrack<4, VarManager::kPsi2StoJpsiPiPi, gkEventFillMapWithCov, gkTrackFillMapWithCov>(event, tracks, dileptons); + } + + void processB0ToJpsiPiPi(soa::Filtered::iterator const& event, MyBarrelTracksSelected const& tracks, soa::Filtered const& dileptons) + { + runDileptonTrackTrack<0, VarManager::kB0toJpsiEEPiPi, gkEventFillMap, gkTrackFillMap>(event, tracks, dileptons); } void processDummy(MyEvents&) diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index e3cf391ff44..996461765d6 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -22,6 +22,7 @@ #include "PWGDQ/Core/MixingLibrary.h" #include "PWGDQ/Core/VarManager.h" #include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/Core/TableHelper.h" @@ -216,6 +217,10 @@ using MyEventsHashSelectedQvector = soa::Join; using MyEventsQvectorCentrSelected = soa::Join; using MyEventsHashSelectedQvectorCentr = soa::Join; +using MyEventsWithColl = soa::Join; +using MyEventsVtxCovWithColl = soa::Join; +using MyEventsWithCollSelected = soa::Join; +using MyEventsVtxCovWithCollSelected = soa::Join; using MyBarrelTracks = soa::Join; using MyBarrelTracksWithAmbiguities = soa::Join; @@ -223,12 +228,17 @@ using MyBarrelTracksWithCov = soa::Join; using MyBarrelTracksWithCovWithAmbiguitiesWithColl = soa::Join; using MyDielectronCandidates = soa::Join; +using MyDielectronCandidatesWithColl = soa::Join; using MyDitrackCandidates = soa::Join; using MyDimuonCandidates = soa::Join; using MyMuonTracks = soa::Join; using MyMuonTracksWithCov = soa::Join; using MyMuonTracksWithCovWithAmbiguities = soa::Join; using MyMuonTracksSelectedWithColl = soa::Join; +using MyV0s = aod::V0Datas; +using MyV0sWithCov = soa::Join; + +using FullTracksExt = soa::Join; // bit maps used for the Fill functions of the VarManager constexpr static uint32_t gkEventFillMap = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended; @@ -242,6 +252,9 @@ constexpr static uint32_t gkEventFillMapWithMultExtraWithQVector = VarManager::O constexpr static uint32_t gkEventFillMapWithMultExtraZdc = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventMultExtra | VarManager::ReducedZdc; constexpr static uint32_t gkEventFillMapWithCovZdcMultExtra = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventVtxCov | VarManager::ReducedZdc | VarManager::ReducedEventMultExtra; constexpr static uint32_t gkEventFillMapWithQvectorCentr = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::CollisionQvect | VarManager::ObjTypes::ReducedEventMultExtra; +constexpr static uint32_t gkEventFillMapWithCovQvectorCentr = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventVtxCov | VarManager::ObjTypes::CollisionQvect | VarManager::ObjTypes::ReducedEventMultExtra; +constexpr static uint32_t gkEventFillMapWithColl = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventCollInfo; +constexpr static uint32_t gkEventFillMapWithCovWithColl = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventVtxCov | VarManager::ObjTypes::ReducedEventCollInfo; // constexpr static uint32_t gkEventFillMapWithQvector = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventQvector; // constexpr static uint32_t gkEventFillMapWithCovQvector = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventVtxCov | VarManager::ObjTypes::ReducedEventQvector; constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; @@ -252,6 +265,8 @@ constexpr static uint32_t gkTrackFillMapWithCovWithColl = VarManager::ObjTypes:: constexpr static uint32_t gkMuonFillMapWithCov = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra | VarManager::ObjTypes::ReducedMuonCov; // constexpr static uint32_t gkMuonFillMapWithColl = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra | VarManager::ObjTypes::ReducedMuonCollInfo; +constexpr static uint32_t gkV0FillMap = VarManager::ObjTypes::TrackV0Bits; + constexpr static int pairTypeEE = VarManager::kDecayToEE; constexpr static int pairTypeMuMu = VarManager::kDecayToMuMu; // constexpr static int pairTypeEMu = VarManager::kElectronMuon; @@ -1231,6 +1246,7 @@ struct AnalysisSameEventPairing { Produces dielectronInfoList; Produces dimuonsExtraList; Produces dielectronAllList; + Produces dielectronVtx; Produces dimuonAllList; Produces dileptonFlowList; Produces dileptonInfoList; @@ -1325,7 +1341,7 @@ struct AnalysisSameEventPairing { void init(o2::framework::InitContext& context) { LOG(info) << "Starting initialization of AnalysisSameEventPairing (idstoreh)"; - fEnableBarrelHistos = context.mOptions.get("processAllSkimmed") || context.mOptions.get("processBarrelOnlySkimmed") || context.mOptions.get("processBarrelOnlyWithCollSkimmed") || context.mOptions.get("processBarrelOnlySkimmedNoCov") || context.mOptions.get("processBarrelOnlySkimmedNoCovWithMultExtra") || context.mOptions.get("processBarrelOnlyWithQvectorCentrSkimmedNoCov"); + fEnableBarrelHistos = context.mOptions.get("processAllSkimmed") || context.mOptions.get("processBarrelOnlySkimmed") || context.mOptions.get("processBarrelOnlyWithCollSkimmed") || context.mOptions.get("processBarrelOnlySkimmedNoCov") || context.mOptions.get("processBarrelOnlySkimmedNoCovWithMultExtra") || context.mOptions.get("processBarrelOnlyWithQvectorCentrSkimmedNoCov") || context.mOptions.get("processPiPiSkimmedNoCov"); fEnableBarrelMixingHistos = context.mOptions.get("processMixingAllSkimmed") || context.mOptions.get("processMixingBarrelSkimmed") || context.mOptions.get("processMixingBarrelSkimmedFlow") || context.mOptions.get("processMixingBarrelWithQvectorCentrSkimmedNoCov"); fEnableMuonHistos = context.mOptions.get("processAllSkimmed") || context.mOptions.get("processMuonOnlySkimmed") || context.mOptions.get("processMuonOnlySkimmedMultExtra") || context.mOptions.get("processMixingMuonSkimmed"); fEnableMuonMixingHistos = context.mOptions.get("processMixingAllSkimmed") || context.mOptions.get("processMixingMuonSkimmed"); @@ -1710,6 +1726,7 @@ struct AnalysisSameEventPairing { dileptonFlowList.reserve(1); if (fConfigOptions.flatTables.value) { dielectronAllList.reserve(1); + dielectronVtx.reserve(1); dimuonAllList.reserve(1); } if (fConfigOptions.polarTables.value) { @@ -1786,9 +1803,9 @@ struct AnalysisSameEventPairing { VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, 0); - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrackCollInfo) > 0) { - dielectronInfoList(t1.collisionId(), t1.trackId(), t2.trackId()); - dileptonInfoList(t1.collisionId(), event.posX(), event.posY(), event.posZ()); + if constexpr ((TEventFillMap & VarManager::ObjTypes::ReducedEventCollInfo) > 0) { + dielectronInfoList(event.collisionId(), t1.trackId(), t2.trackId()); + dileptonInfoList(event.collisionId(), event.posX(), event.posY(), event.posZ()); } if (fConfigOptions.polarTables.value) { dileptonPolarList(VarManager::fgValues[VarManager::kCosThetaHE], VarManager::fgValues[VarManager::kPhiHE], VarManager::fgValues[VarManager::kPhiTildeHE], @@ -1838,17 +1855,23 @@ struct AnalysisSameEventPairing { continue; if (fConfigOptions.flatTables.value) { - dielectronAllList(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, dileptonMcDecision, + // dielectronAllList(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, dileptonMcDecision, + // t1.pt(), t1.eta(), t1.phi(), t1.itsClusterMap(), t1.itsChi2NCl(), t1.tpcNClsCrossedRows(), t1.tpcNClsFound(), t1.tpcChi2NCl(), t1.dcaXY(), t1.dcaZ(), t1.tpcSignal(), t1.tpcNSigmaEl(), t1.tpcNSigmaPi(), t1.tpcNSigmaPr(), t1.beta(), t1.tofNSigmaEl(), t1.tofNSigmaPi(), t1.tofNSigmaPr(), + // t2.pt(), t2.eta(), t2.phi(), t2.itsClusterMap(), t2.itsChi2NCl(), t2.tpcNClsCrossedRows(), t2.tpcNClsFound(), t2.tpcChi2NCl(), t2.dcaXY(), t2.dcaZ(), t2.tpcSignal(), t2.tpcNSigmaEl(), t2.tpcNSigmaPi(), t2.tpcNSigmaPr(), t2.beta(), t2.tofNSigmaEl(), t2.tofNSigmaPi(), t2.tofNSigmaPr(), + // VarManager::fgValues[VarManager::kKFTrack0DCAxyz], VarManager::fgValues[VarManager::kKFTrack1DCAxyz], VarManager::fgValues[VarManager::kKFDCAxyzBetweenProngs], VarManager::fgValues[VarManager::kKFTrack0DCAxy], VarManager::fgValues[VarManager::kKFTrack1DCAxy], VarManager::fgValues[VarManager::kKFDCAxyBetweenProngs], + // VarManager::fgValues[VarManager::kKFTrack0DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack0DeviationxyFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationxyFromPV], + // VarManager::fgValues[VarManager::kKFMass], VarManager::fgValues[VarManager::kKFChi2OverNDFGeo], VarManager::fgValues[VarManager::kVertexingLxyz], VarManager::fgValues[VarManager::kVertexingLxyzOverErr], VarManager::fgValues[VarManager::kVertexingLxy], VarManager::fgValues[VarManager::kVertexingLxyOverErr], VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxyErr], VarManager::fgValues[VarManager::kKFCosPA], VarManager::fgValues[VarManager::kKFJpsiDCAxyz], VarManager::fgValues[VarManager::kKFJpsiDCAxy], + // VarManager::fgValues[VarManager::kKFPairDeviationFromPV], VarManager::fgValues[VarManager::kKFPairDeviationxyFromPV], + // VarManager::fgValues[VarManager::kKFMassGeoTop], + // VarManager::fgValues[VarManager::kKFChi2OverNDFGeoTop], VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingTauxyProjected], VarManager::fgValues[VarManager::kVertexingLzProjected], VarManager::fgValues[VarManager::kVertexingLxyProjected]); + dielectronVtx(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, dileptonMcDecision, t1.pt(), t1.eta(), t1.phi(), t1.itsClusterMap(), t1.itsChi2NCl(), t1.tpcNClsCrossedRows(), t1.tpcNClsFound(), t1.tpcChi2NCl(), t1.dcaXY(), t1.dcaZ(), t1.tpcSignal(), t1.tpcNSigmaEl(), t1.tpcNSigmaPi(), t1.tpcNSigmaPr(), t1.beta(), t1.tofNSigmaEl(), t1.tofNSigmaPi(), t1.tofNSigmaPr(), t2.pt(), t2.eta(), t2.phi(), t2.itsClusterMap(), t2.itsChi2NCl(), t2.tpcNClsCrossedRows(), t2.tpcNClsFound(), t2.tpcChi2NCl(), t2.dcaXY(), t2.dcaZ(), t2.tpcSignal(), t2.tpcNSigmaEl(), t2.tpcNSigmaPi(), t2.tpcNSigmaPr(), t2.beta(), t2.tofNSigmaEl(), t2.tofNSigmaPi(), t2.tofNSigmaPr(), - VarManager::fgValues[VarManager::kKFTrack0DCAxyz], VarManager::fgValues[VarManager::kKFTrack1DCAxyz], VarManager::fgValues[VarManager::kKFDCAxyzBetweenProngs], VarManager::fgValues[VarManager::kKFTrack0DCAxy], VarManager::fgValues[VarManager::kKFTrack1DCAxy], VarManager::fgValues[VarManager::kKFDCAxyBetweenProngs], - VarManager::fgValues[VarManager::kKFTrack0DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationFromPV], VarManager::fgValues[VarManager::kKFTrack0DeviationxyFromPV], VarManager::fgValues[VarManager::kKFTrack1DeviationxyFromPV], - VarManager::fgValues[VarManager::kKFMass], VarManager::fgValues[VarManager::kKFChi2OverNDFGeo], VarManager::fgValues[VarManager::kVertexingLxyz], VarManager::fgValues[VarManager::kVertexingLxyzOverErr], VarManager::fgValues[VarManager::kVertexingLxy], VarManager::fgValues[VarManager::kVertexingLxyOverErr], VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxyErr], VarManager::fgValues[VarManager::kKFCosPA], VarManager::fgValues[VarManager::kKFJpsiDCAxyz], VarManager::fgValues[VarManager::kKFJpsiDCAxy], - VarManager::fgValues[VarManager::kKFPairDeviationFromPV], VarManager::fgValues[VarManager::kKFPairDeviationxyFromPV], - VarManager::fgValues[VarManager::kKFMassGeoTop], - VarManager::fgValues[VarManager::kKFChi2OverNDFGeoTop], VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingTauxyProjected], VarManager::fgValues[VarManager::kVertexingLzProjected], VarManager::fgValues[VarManager::kVertexingLxyProjected]); + VarManager::fgValues[VarManager::kVertexingLxy], VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxy]/VarManager::fgValues[VarManager::kVertexingTauxyErr], VarManager::fgValues[VarManager::kVertexingLz], VarManager::fgValues[VarManager::kVertexingTauz],VarManager::fgValues[VarManager::kVertexingTauz]/VarManager::fgValues[VarManager::kVertexingTauzErr]); } } + } else { + dielectronsExtraList(t1.globalIndex(), t2.globalIndex(), -999.f, -999.f, -999.f); } } @@ -2283,11 +2306,18 @@ struct AnalysisSameEventPairing { runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); } - void processBarrelOnlyWithCollSkimmed(MyEventsVtxCovSelected const& events, + void processBarrelOnlyWithCollSkimmed(MyEventsVtxCovWithCollSelected const& events, + soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguitiesWithColl const& barrelTracks) + { + runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); + } + + void processBarrelOnlyWithCollAllSkimmed(MyEventsVtxCovWithCollSelected const& events, soa::Join const& barrelAssocs, MyBarrelTracksWithCovWithAmbiguitiesWithColl const& barrelTracks) { - runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); + runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); } void processBarrelOnlyWithQvectorCentrSkimmedNoCov(MyEventsQvectorCentrSelected const& events, @@ -3746,6 +3776,314 @@ struct AnalysisDileptonTrack { PROCESS_SWITCH(AnalysisDileptonTrack, processDummy, "Dummy function", true); }; +struct AnalysisDileptonV0 { + OutputObj fOutputList{"output"}; + + Configurable fConfigDileptonLowMass{"cfgDileptonLowMass", 2.8, "Low mass cut for the dileptons used in analysis"}; + Configurable fConfigDileptonHighMass{"cfgDileptonHighMass", 3.2, "High mass cut for the dileptons used in analysis"}; + Configurable fConfigDileptonpTCut{"cfgDileptonpTCut", 0.0, "pT cut for dileptons used in the triplet vertexing"}; + Configurable fConfigDileptonLxyCut{"cfgDileptonLxyCut", 0.0, "Lxy cut for dileptons used in the triplet vertexing"}; + // Configurable fConfigV0Type{"cfgV0Type", -1, "V0 type"}; + Configurable fConfigUseKFVertexing{"cfgUseKFVertexing", false, "Use KF Particle for secondary vertex reconstruction (DCAFitter is used by default)"}; + + // Configurable fConfigHistogramSubgroups{"cfgDileptonTrackHistogramsSubgroups", "invmass,vertexing", "Comma separated list of dilepton-track histogram subgroups"}; + // Configurable fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"}; + Configurable fConfigMixingDepth{"cfgMixingDepth", 5, "Event mixing pool depth"}; + + struct : ConfigurableGroup { + Configurable V0MassLow{"V0MassLow", 0.48, "Lower mass cut for V0 candidates"}; + Configurable V0MassHigh{"V0MassHigh", 0.52, "Upper mass cut for V0 candidates"}; + Configurable CutTPCnSigPi{"CutTPCnSigPi", 4.0, "TPC nSigma cut for V0 daughter pions"}; + Configurable CutTPCncls{"CutTPCncls", 70.0, "TPC nClusters cut for V0 daughter tracks"}; + Configurable CutTPCChi2Ncls{"CutTPCChi2Ncls", 4.0, "TPC chi2 per nClusters cut for V0 daughter tracks"}; + Configurable CutDCADaughtersToPV{"CutDCADaughtersToPV", 0.06, "DCA cut between V0 daughter tracks"}; + Configurable CutCosPointingAngle{"CutCosPointingAngle", 0.98, "Cosine of pointing angle cut for V0 candidates"}; + // Configurable fConfigV0Type{"cfgV0Type", -1, "V0 type"}; + // Configurable DaughterEtaCut{"DaughterEtaCut", 0.9, "Eta cut for V0 daughter tracks"}; + } fConfigV0Cuts; + + Configurable fConfigUseRemoteField{"cfgUseRemoteField", false, "Chose whether to fetch the magnetic field from ccdb or set it manually"}; + Configurable fConfigGRPmagPath{"cfgGrpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable fConfigMagField{"cfgMagField", 5.0f, "Manually set magnetic field"}; + + Configurable fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable fConfigNoLaterThan{"ccdb-no-later-than", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; + Configurable fConfigGeoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + + Filter eventFilter = aod::dqanalysisflags::isEventSelected > static_cast(0); + Filter dileptonFilter = aod::reducedpair::sign == 0 && aod::reducedpair::mass > fConfigDileptonLowMass && aod::reducedpair::mass < fConfigDileptonHighMass; + // Filter filterV0Selected = aod::v0mapID::v0addid == fConfigV0Type; + // Filter filterV0Selected = aod::v0data::mK0Short > fConfigV0Cuts.V0MassLow && aod::v0data::mK0Short < fConfigV0Cuts.V0MassHigh; + + constexpr static uint32_t fgDileptonFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::Pair; // fill map + + float* fValuesV0; + float* fValuesDilepton; + + o2::vertexing::DCAFitterN<2> df2; + // HistogramManager* fHistMan; + + // Produces DileptonV0Table; + + // tempt: use registry to check inv mass of dilepton-v0 + HistogramRegistry registry{"registry"}; + void init(o2::framework::InitContext& context) + { + if (context.mOptions.get("processDummy")) { + return; + } + + df2.setPropagateToPCA(true); + df2.setBz(5.0f); // temporary, will be set properly in initParamsFromCCDB + df2.setMaxR(200.0f); + df2.setMaxDZIni(4.0f); + df2.setMinParamChange(1.0e-3f); + df2.setMinRelChi2Change(0.9f); + df2.setUseAbsDCA(true); + + fValuesV0 = new float[VarManager::kNVars]; + fValuesDilepton = new float[VarManager::kNVars]; + VarManager::SetDefaultVarNames(); + + // TO Do: add initialization of fHistMan + // fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + // fHistMan->SetUseDefaultVariableNames(kTRUE); + // fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + // --------------------------------------------- + + registry.add("hInvMassDileptonV0", "hInvMassDileptonV0", HistType::kTH1F, {{2000, 4.0, 6.0}}); + registry.add("hInvMassDilepton", "hInvMassDilepton", HistType::kTH1F, {{2000, 2.0, 4.0}}); + registry.add("hPtV0", "hPtV0", HistType::kTH1F, {{2000, 0.0, 10.0}}); + registry.add("hPtV0Pos", "hPtV0Pos", HistType::kTH1F, {{2000, 0.0, 10.0}}); + registry.add("hPosDCAToPV", "hPosDCAToPV", HistType::kTH1F, {{1000, -5.0, 5.0}}); + registry.add("hNegDCAToPV", "hNegDCAToPV", HistType::kTH1F, {{1000, -5.0, 5.0}}); + // registry.add("") + registry.add("hV0Selected", "hV0Selected", HistType::kTH1F, {{100, 0.45, 0.55}}); + registry.add("hIsV0Vertexing", "hIsV0Vertexing", HistType::kTH1F, {{2, -0.5, 1.5}}); + registry.add("hIsDileptonVertexing", "hIsDileptonVertexing", HistType::kTH1F, {{2, -0.5, 1.5}}); + registry.add("hAssociatedCollisions", "hAssociatedCollisions", HistType::kTH1F, {{100, -0.5, 99.5}}); + } + + template + void runDileptonV0(TEvent const& event, TTracks const& tracks, TV0s const& v0s, TFullTracks const& fulltracks, TDileptons const& dileptons) + { + // VarManager::ResetValues(0, VarManager::kNVars, fValuesDilepton); + // VarManager::ResetValues(0, VarManager::kNVars, fValuesV0); + // VarManager::FillEvent(event, fValuesDilepton); + // VarManager::FillEvent(event, fValuesV0); + // LOGP(info, "event in collision: {}", event.collisionId()); + + for (auto dilepton : dileptons) { + if (event.collisionId() != dilepton.collisionId()) { + continue; + } + LOGP(info, "dilepton in collision: {}", dilepton.collisionId()); + // get full track info of tracks based on the index + int indexLepton1 = dilepton.index0Id(); + int indexLepton2 = dilepton.index1Id(); + auto lepton1 = tracks.rawIteratorAt(dilepton.index0Id()); + auto lepton2 = tracks.rawIteratorAt(dilepton.index1Id()); + // Check that the dilepton has zero charge + if (dilepton.sign() != 0) { + continue; + } + // make sure lepton1 is positive and lepton2 is negative + // if (lepton1.sign() < 0 && lepton2.sign() > 0) { + // std::swap(indexLepton1, indexLepton2); + // std::swap(lepton1, lepton2); + // } + // VarManager::FillTrack(lepton1, fValuesDilepton); + registry.fill(HIST("hInvMassDilepton"), dilepton.mass()); + + for (auto v0 : v0s) { + LOGP(info, "v0 in collision: {}", v0.collisionId()); + // apply basic V0 cuts + registry.fill(HIST("hPosDCAToPV"), v0.dcapostopv()); + registry.fill(HIST("hNegDCAToPV"), v0.dcanegtopv()); + if (std::abs(v0.dcapostopv()) < fConfigV0Cuts.CutDCADaughtersToPV || std::abs(v0.dcanegtopv()) < fConfigV0Cuts.CutDCADaughtersToPV) { + continue; + } + if (v0.v0cosPA() < fConfigV0Cuts.CutCosPointingAngle) { + continue; + } + if constexpr (TCandidateType == VarManager::kBtoJpsiEEK0S) { + auto v0_pos = v0.template posTrack_as(); + auto v0_neg = v0.template negTrack_as(); + + // apply daughter track cuts + if (v0_pos.tpcNClsFound() < fConfigV0Cuts.CutTPCncls || v0_neg.tpcNClsFound() < fConfigV0Cuts.CutTPCncls) { + continue; + } + if (std::abs(v0_pos.tpcNSigmaPi()) > fConfigV0Cuts.CutTPCnSigPi || std::abs(v0_neg.tpcNSigmaPi()) > fConfigV0Cuts.CutTPCnSigPi) { + continue; + } + if (v0_pos.tpcChi2NCl() > fConfigV0Cuts.CutTPCChi2Ncls || v0_neg.tpcChi2NCl() > fConfigV0Cuts.CutTPCChi2Ncls) { + continue; + } + + + float hadronMass = o2::constants::physics::MassKaonNeutral; + float daughterMass = o2::constants::physics::MassPionCharged; + // Check V0 mass window + ROOT::Math::PtEtaPhiMVector v0vec(v0_pos.pt(), v0_pos.eta(), v0_pos.phi(), daughterMass); + ROOT::Math::PtEtaPhiMVector v1vec(v0_neg.pt(), v0_neg.eta(), v0_neg.phi(), daughterMass); + ROOT::Math::PtEtaPhiMVector v0sum = v0vec + v1vec; + // if (v0sum.M() < fConfigV0Cuts.V0MassLow || v0sum.M() > fConfigV0Cuts.V0MassHigh) { + // continue; + // } + if (v0.mK0Short() < fConfigV0Cuts.V0MassLow || v0.mK0Short() > fConfigV0Cuts.V0MassHigh) { + continue; + } + registry.fill(HIST("hV0Selected"), v0.mK0Short()); + // VarManager::FillDileptonHadron(dilepton, v0, fValuesV0, hadronMass); + + //Reconstruct directly here + ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass()); + ROOT::Math::PtEtaPhiMVector v2(v0.pt(), v0.eta(), v0.phi(), hadronMass); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float invmass = v12.M(); + // float massJpsi = 3.0969; + float mass = invmass - dilepton.mass() + 3.096; + registry.fill(HIST("hInvMassDileptonV0"), mass); + registry.fill(HIST("hPtV0"), v0.pt()); + registry.fill(HIST("hPtV0Pos"), v0_pos.pt()); + + if constexpr (TThreeProngFitter) { + LOGP(info, "Three prong fitter under development"); + int indexTrack1 = lepton1.trackId(); + int indexTrack2 = lepton2.trackId(); + auto track1 = fulltracks.rawIteratorAt(indexTrack1); + auto track2 = fulltracks.rawIteratorAt(indexTrack2); + // recontruct V0 + try { + if (!df2.process(getTrackParCov(v0_pos), getTrackParCov(v0_neg))) { + registry.fill(HIST("hIsV0Vertexing"), 0); + continue; + } + } catch (const std::runtime_error& error) { + LOG(info) << "Run time error found: " << error.what() << ". DCAFitterN for V0-bachelor cannot work, skipping the candidate."; + continue; + } + o2::track::TrackParCov v0par = df2.createParentTrackParCov(0); + registry.fill(HIST("hIsV0Vertexing"), 1); + // reconstruct dilepton + try { + if (!df2.process(getTrackParCov(lepton1), getTrackParCov(lepton2))) { + registry.fill(HIST("hIsDileptonVertexing"), 0); + continue; + } + } catch (const std::runtime_error& error) { + LOG(info) << "Run time error found: " << error.what() << ". DCAFitterN for V0-bachelor cannot work, skipping the candidate."; + continue; + } + registry.fill(HIST("hIsDileptonVertexing"), 1); + o2::track::TrackParCov dileptonpar = df2.createParentTrackParCov(0); + // reconstruct B0 + try { + if (!df2.process(v0par, dileptonpar)) { + continue; + } + } catch (const std::runtime_error& error) { + LOG(info) << "Run time error found: " << error.what() << ". DCAFitterN for V0-bachelor cannot work, skipping the candidate."; + continue; + } + o2::track::TrackParCov b0par = df2.createParentTrackParCov(0); + Vec3D secondaryVertex = df2.getPCACandidate(0); + o2::math_utils::Point3D vtxXYZ(event.posX(), event.posY(), event.posZ()); + std::array vtxCov{event.covXX(), event.covXY(), event.covYY(), event.covXZ(), event.covYZ(), event.covZZ()}; + o2::dataformats::VertexBase primaryVertex = {std::move(vtxXYZ), std::move(vtxCov)}; + auto covMatrixPV = primaryVertex.getCov(); + // TODO: calculate decay length + } // secondary vertex reconstruction + } + } // loop over v0s + } // loop over dileptons + } + + Preslice dileptonsPerCollision = aod::reducedpair::collisionId; + Preslice v0sPerCollision = aod::v0data::collisionId; + + void processB0ToJpsiK0s(soa::Filtered const& events, + MyBarrelTracks const& tracks, + MyV0s const& v0s, FullTracksExt const& fulltracks, + soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + + for (auto event : events) { + auto groupedDileptons = dileptons.sliceBy(dileptonsPerCollision, event.collisionId()); + auto groupedV0s = v0s.sliceBy(v0sPerCollision, event.collisionId()); + runDileptonV0(event, tracks, groupedV0s, fulltracks, groupedDileptons); + } + } + + void processV0(soa::Filtered const& events, + MyV0s const& v0s) + { + if (events.size() == 0) { + return; + } + + for (auto event : events) { + auto groupedV0s = v0s.sliceBy(v0sPerCollision, event.collisionId()); + for (auto v0 : groupedV0s) { + registry.fill(HIST("hPosDCAToPV"), v0.dcapostopv()); + registry.fill(HIST("hNegDCAToPV"), v0.dcanegtopv()); + if (std::abs(v0.dcapostopv()) < fConfigV0Cuts.CutDCADaughtersToPV || std::abs(v0.dcanegtopv()) < fConfigV0Cuts.CutDCADaughtersToPV) { + continue; + } + if (v0.v0cosPA() < fConfigV0Cuts.CutCosPointingAngle) { + continue; + } + if (v0.mK0Short() < fConfigV0Cuts.V0MassLow || v0.mK0Short() > fConfigV0Cuts.V0MassHigh) { + continue; + } + registry.fill(HIST("hV0Selected"), v0.mK0Short()); + } + } + } + + Preslice dielectronsPerEvent = aod::reducedpair::reducedeventId; + void processDilepton(soa::Filtered const& events, + soa::Filtered const& dileptons, + MyBarrelTracksWithCovWithAmbiguitiesWithColl const& tracks) + { + if (events.size() == 0) { + return; + } + + for (auto event : events) { + auto groupedDileptons = dileptons.sliceBy(dielectronsPerEvent, event.globalIndex()); + for (auto dilepton : groupedDileptons) { + if (dilepton.sign() != 0) { + continue; + } + auto lepton1 = tracks.rawIteratorAt(dilepton.index0Id()); + auto lepton2 = tracks.rawIteratorAt(dilepton.index1Id()); + LOGP(info, "event in collision: {}", event.collisionId()); + LOGP(info, "dilepton in collision: {}", dilepton.collisionId()); + LOGP(info, "track1 in collision: {}", lepton1.collisionId()); + LOGP(info, "track2 in collision: {}", lepton2.collisionId()); + LOGP(info, "event global index: {}", event.globalIndex()); + LOGP(info, "event Id: {} =? {}", lepton1.reducedeventId(), lepton2.reducedeventId()); + LOGP(info, "***********************************************************************************************"); + } + } + } + + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisDileptonV0, processB0ToJpsiK0s, "Run B0 to J/psi K0s analysis", false); + PROCESS_SWITCH(AnalysisDileptonV0, processV0, "Run V0 analysis", false); + PROCESS_SWITCH(AnalysisDileptonV0, processDilepton, "Run dilepton analysis", false); + PROCESS_SWITCH(AnalysisDileptonV0, processDummy, "Dummy function", true); +}; + struct AnalysisDileptonTrackTrack { OutputObj fOutputList{"output"}; @@ -3756,6 +4094,8 @@ struct AnalysisDileptonTrackTrack { Configurable fConfigAddDileptonHistogram{"cfgAddDileptonHistogram", "barrel", "Comma separated list of histograms"}; Configurable fConfigAddQuadrupletHistogram{"cfgAddQuadrupletHistogram", "xtojpsipipi", "Comma separated list of histograms"}; + Configurable fConfigMixingDepth{"cfgMixingDepth", 5, "Event mixing pool depth"}; + Configurable fConfigSetupFourProngFitter{"cfgSetupFourProngFitter", false, "Use DCA for secondary vertex reconstruction (DCAFitter is used by default)"}; Configurable fConfigUseKFVertexing{"cfgUseKFVertexing", false, "Use KF Particle for secondary vertex reconstruction (DCAFitter is used by default)"}; Configurable fConfigUseRemoteField{"cfgUseRemoteField", false, "Chose whether to fetch the magnetic field from ccdb or set it manually"}; @@ -3786,6 +4126,8 @@ struct AnalysisDileptonTrackTrack { float* fValuesQuadruplet; HistogramManager* fHistMan; + NoBinningPolicy fHashBin; + void init(o2::framework::InitContext& context) { // bool isBarrel = context.mOptions.get("processJpsiPiPi"); @@ -3832,6 +4174,9 @@ struct AnalysisDileptonTrackTrack { DefineHistograms(fHistMan, Form("QuadrupletSEMM_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); } } + if (context.mOptions.get("processBarrelMixedEvent")) { + DefineHistograms(fHistMan, Form("QuadrupletME"),"mixedevent"); + } } VarManager::SetUseVars(fHistMan->GetUsedVars()); @@ -3839,7 +4184,7 @@ struct AnalysisDileptonTrackTrack { } // init parameters from CCDB - void initParamsFromCCDB(uint64_t timestamp) + void initParamsFromCCDB(uint64_t timestamp, int NProng) { if (fConfigUseRemoteField.value) { o2::parameters::GRPMagField* grpmag = fCCDB->getForTimeStamp(fConfigGRPmagPath.value, timestamp); @@ -3849,27 +4194,55 @@ struct AnalysisDileptonTrackTrack { } else { LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", timestamp); } - if (fConfigUseKFVertexing.value) { - VarManager::SetupTwoProngKFParticle(magField); + // if (fConfigUseKFVertexing.value) { + // VarManager::SetupTwoProngKFParticle(magField); + // VarManager::SetupFourProngKFParticle(magField); + // } else if (fConfigSetupFourProngFitter.value) { + // VarManager::SetupTwoProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // VarManager::SetupThreeProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); + // VarManager::SetupFourProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // } + if (NProng == 3) { + if (fConfigUseKFVertexing) { + VarManager::SetupThreeProngKFParticle(magField); + } else { + VarManager::SetupThreeProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } + } else if (NProng == 4) { + if (fConfigUseKFVertexing) { VarManager::SetupFourProngKFParticle(magField); - } else if (fConfigSetupFourProngFitter.value) { - VarManager::SetupTwoProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } else { VarManager::SetupFourProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } } } else { - if (fConfigUseKFVertexing.value) { - VarManager::SetupTwoProngKFParticle(fConfigMagField.value); + // if (fConfigUseKFVertexing.value) { + // VarManager::SetupTwoProngKFParticle(fConfigMagField.value); + // VarManager::SetupFourProngKFParticle(fConfigMagField.value); + // } else if (fConfigSetupFourProngFitter.value) { + // LOGP(info, "Setting up DCA fitter for two and four prong candidates"); + // VarManager::SetupTwoProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // VarManager::SetupThreeProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); + // VarManager::SetupFourProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + // } + if (NProng == 3) { + if (fConfigUseKFVertexing) { + VarManager::SetupThreeProngKFParticle(fConfigMagField.value); + } else { + VarManager::SetupThreeProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } + } else if (NProng == 4) { + if (fConfigUseKFVertexing) { VarManager::SetupFourProngKFParticle(fConfigMagField.value); - } else if (fConfigSetupFourProngFitter.value) { - LOGP(info, "Setting up DCA fitter for two and four prong candidates"); - VarManager::SetupTwoProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } else { VarManager::SetupFourProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, false); // TODO: get these parameters from Configurables + } } } } // Template function to run pair - hadron combinations - template + template void runDileptonTrackTrack(TEvent const& event, TTrackAssocs const& assocs, TTracks const& tracks, TDileptons const& dileptons) { VarManager::ResetValues(0, VarManager::kNVars, fValuesQuadruplet); @@ -3888,6 +4261,10 @@ struct AnalysisDileptonTrackTrack { if (dilepton.sign() != 0) { continue; } + // if (lepton1.sign() < 0 && lepton2.sign() > 0) { + // std::swap(indexLepton1, indexLepton2); + // std::swap(lepton1, lepton2); + // } VarManager::FillTrack(dilepton, fValuesQuadruplet); // LOGP(info, "is dilepton selected: {}", fDileptonCut.IsSelected(fValuesQuadruplet)); @@ -3920,7 +4297,7 @@ struct AnalysisDileptonTrackTrack { // fill variables VarManager::FillDileptonTrackTrack(dilepton, track1, track2, fValuesQuadruplet); - if (fConfigSetupFourProngFitter || fConfigUseKFVertexing) { + if constexpr (NProng > 0) { // LOGP(info, "Using KF or DCA fitter for secondary vertexing"); VarManager::FillDileptonTrackTrackVertexing(event, lepton1, lepton2, track1, track2, fValuesQuadruplet); } @@ -3934,6 +4311,9 @@ struct AnalysisDileptonTrackTrack { if (fIsSameTrackCut) { if (track1.sign() * track2.sign() < 0) { fHistMan->FillHistClass(Form("QuadrupletSEPM_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + // if (track1.sign() < 0 && track2.sign() > 0) { + // std::swap(track1, track2); + // } } } else { if ((track1.sign() < 0) && (track2.sign() > 0)) { @@ -3970,7 +4350,7 @@ struct AnalysisDileptonTrackTrack { Preslice dielectronsPerCollision = aod::reducedpair::reducedeventId; Preslice ditracksPerCollision = aod::reducedpair::reducedeventId; - void processJpsiPiPi(soa::Filtered const& events, + void processChicToJpsiPiPi(soa::Filtered const& events, soa::Filtered> const& assocs, MyBarrelTracksWithCov const& tracks, soa::Filtered const& dileptons) { @@ -3978,13 +4358,147 @@ struct AnalysisDileptonTrackTrack { return; } if (fCurrentRun != events.begin().runNumber()) { // start: runNumber - initParamsFromCCDB(events.begin().timestamp()); + initParamsFromCCDB(events.begin().timestamp(), 4); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack<4, VarManager::kXtoJpsiPiPi, gkEventFillMapWithCov, gkTrackFillMapWithCov>(event, groupedBarrelAssocs, tracks, groupedDielectrons); + } + } + + void processPsi2SToJpsiPiPi(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracksWithCov const& tracks, soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp(), 4); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack<4, VarManager::kPsi2StoJpsiPiPi, gkEventFillMapWithCov, gkTrackFillMapWithCov>(event, groupedBarrelAssocs, tracks, groupedDielectrons); + } + } + + void processB0ToJpsiPiPi(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracks const& tracks, soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp(), 0); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack<0, VarManager::kB0toJpsiEEPiPi, gkEventFillMap, gkTrackFillMap>(event, groupedBarrelAssocs, tracks, groupedDielectrons); + } + } + + void processBsToJpsiKK(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracks const& tracks, soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp(), 0); fCurrentRun = events.begin().runNumber(); } // end: runNumber for (auto& event : events) { auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); - runDileptonTrackTrack(event, groupedBarrelAssocs, tracks, groupedDielectrons); + runDileptonTrackTrack<0, VarManager::kBstoJpsiEEPhi, gkEventFillMap, gkTrackFillMap>(event, groupedBarrelAssocs, tracks, groupedDielectrons); + } + } + + void processJpsiPiPiNoAmbi(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracksWithCovWithAmbiguities const& tracks, soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp(), 4); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack<4, VarManager::kPsi2StoJpsiPiPi, gkEventFillMapWithCov, gkTrackFillMapWithCovAmbiguities>(event, groupedBarrelAssocs, tracks, groupedDielectrons); + } + } + + void processBarrelMixedEvent(soa::Filtered& events, + soa::Filtered> const& assocs, + MyBarrelTracks const& tracks, soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + events.bindExternalIndices(&dileptons); + events.bindExternalIndices(&assocs); + + // loop over two event comibnations + for (auto& [event1, event2] : selfCombinations(fHashBin, fConfigMixingDepth.value, -1, events, events)) { + // fill event quantities + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event1, VarManager::fgValues); + + // get the dilepton slice for event1 + auto evDileptons = dileptons.sliceBy(dielectronsPerCollision, event1.globalIndex()); + evDileptons.bindExternalIndices(&events); + + // get the track associations slice for event2 + auto evAssocs = assocs.sliceBy(trackAssocsPerCollision, event2.globalIndex()); + evAssocs.bindExternalIndices(&events); + + for (auto& [a1, a2] : o2::soa::combinations(evAssocs, evAssocs)) { + uint32_t trackSelection = 0; + if (fIsSameTrackCut) { + trackSelection = ((a1.isBarrelSelected_raw() & (static_cast(1) << 1)) || (a2.isBarrelSelected_raw() & (static_cast(1) << 1))); + } else { + trackSelection = ((a1.isBarrelSelected_raw() & (static_cast(1) << 1)) && (a2.isBarrelSelected_raw() & (static_cast(1) << 2))); + } + // LOGP(info, "trackSelection: {}, a1: {}, a2: {}", trackSelection, a1.isBarrelSelected_raw(), a2.isBarrelSelected_raw()); + if (!trackSelection) { + continue; + } + + // get the track from this association + auto track1 = a1.template reducedtrack_as(); + auto track2 = a2.template reducedtrack_as(); + // avoid self combinations + if (track1.globalIndex() == track2.globalIndex()) { + continue; + } + + // loop over dileptons in event1 + for (auto dilepton : evDileptons) { + + // Check that the dilepton has zero charge + if (dilepton.sign() != 0) { + continue; + } + + if (!fDileptonCut.IsSelected(fValuesQuadruplet)) + continue; + + VarManager::FillDileptonTrackTrack(dilepton, track1, track2, VarManager::fgValues); + fHistMan->FillHistClass("QuadrupletME", VarManager::fgValues); + } + } } } @@ -3993,7 +4507,11 @@ struct AnalysisDileptonTrackTrack { // do nothing } - PROCESS_SWITCH(AnalysisDileptonTrackTrack, processJpsiPiPi, "Run barrel pairing of J/psi with pion candidate", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processPsi2SToJpsiPiPi, "Run barrel pairing of J/psi with pion candidate", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processB0ToJpsiPiPi, "Run barrel pairing of B0 to J/psi with pion candidates", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processBsToJpsiKK, "Run barrel pairing of Bs to J/psi with kaon candidates", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processJpsiPiPiNoAmbi, "Run barrel pairing of J/psi with pion candidate without ambiguities", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processBarrelMixedEvent, "Run barrel mixing", false); PROCESS_SWITCH(AnalysisDileptonTrackTrack, processDummy, "Dummy function", true); }; @@ -4007,6 +4525,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; } @@ -4088,5 +4607,9 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char if (classStr.Contains("Quadruplet")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-dihadron", histName); } + + if (classStr.Contains("QuadrupletME")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-dihadron", "mixedevent"); + } } // end loop over histogram classes } From 661c0c7d463ad827f7cca114a9b654cbc60f95be Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 12 Oct 2025 09:19:32 +0200 Subject: [PATCH 2/8] cuts --- PWGDQ/Core/CutsLibrary.cxx | 71 +++++++------------------------------- 1 file changed, 13 insertions(+), 58 deletions(-) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index 6690ffa56b1..7691e58a20c 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -3541,6 +3541,11 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + if (!nameStr.compare("pairJpsi4")) { + cut->AddCut(GetAnalysisCut("pairJpsi4")); + return cut; + } + if (!nameStr.compare("pairPsi2S")) { cut->AddCut(GetAnalysisCut("pairPsi2S")); return cut; @@ -3709,13 +3714,13 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) } if (!nameStr.compare("pairJpsiLowPt1")) { - cut->AddCut(GetAnalysisCut("pairJpsi")); + cut->AddCut(GetAnalysisCut("pairJpsi2")); cut->AddCut(GetAnalysisCut("pairPtLow1")); return cut; } if (!nameStr.compare("pairJpsiLowPt2")) { - cut->AddCut(GetAnalysisCut("pairJpsi")); + cut->AddCut(GetAnalysisCut("pairJpsi2")); cut->AddCut(GetAnalysisCut("pairPtLow2")); return cut; } @@ -4601,7 +4606,7 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) if (!nameStr.compare("jpsi_trackCut_debug6")) { cut->AddCut(VarManager::kTPCchi2, 0.0, 4.0); cut->AddCut(VarManager::kTPCncls, 120., 159); - cut->AddCut(VarManager::kTPCnclsCR, 140., 159); + // cut->AddCut(VarManager::kTPCnclsCR, 140., 159); cut->AddCut(VarManager::kIsITSibAny, 0.5, 1.5); cut->AddCut(VarManager::kIsSPDfirst, 0.5, 1.5); return cut; @@ -6708,6 +6713,11 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("pairJpsi4")) { + cut->AddCut(VarManager::kMass, 2.8, 3.2); + return cut; + } + if (!nameStr.compare("pairPsi2S")) { cut->AddCut(VarManager::kMass, 3.4, 3.9); return cut; @@ -6739,8 +6749,6 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } -<<<<<<< HEAD -======= if (!nameStr.compare("pairX3872_minitree")) { cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); cut->AddCut(VarManager::kQ, -0.5, 0.5); @@ -6758,13 +6766,8 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); cut->AddCut(VarManager::kDitrackMass, 0.0, 0.8); cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); -<<<<<<< HEAD - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) return cut; } @@ -6776,13 +6779,8 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); cut->AddCut(VarManager::kDitrackMass, 0.2, 0.6); cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); -<<<<<<< HEAD - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) return cut; } @@ -6794,13 +6792,8 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); // cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); -<<<<<<< HEAD - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) return cut; } @@ -6812,13 +6805,8 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); // cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); -<<<<<<< HEAD - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) cut->AddCut(VarManager::kPtOverPairPt, 1.0, 1.5); return cut; } @@ -6831,13 +6819,8 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); // cut->AddCut(VarManager::kQuadPt, 2.0, 1000.0); -<<<<<<< HEAD - // cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - // cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= // cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); // cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) cut->AddCut(VarManager::kPtOverPairPt, 1.0, 1.5); return cut; } @@ -6845,34 +6828,20 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) if (!nameStr.compare("pairX3872_13")) { cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); cut->AddCut(VarManager::kQ, -0.2, 0.2); -<<<<<<< HEAD - cut->AddCut(VarManager::kDeltaR, 0.0, 5.0); - cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); cut->AddCut(VarManager::kDitrackMass, 0.45, 0.55); cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) return cut; } if (!nameStr.compare("pairX3872_14")) { cut->AddCut(VarManager::kQuadDefaultDileptonMass, 3.0, 5.0); cut->AddCut(VarManager::kQ, -0.4, 0.2); -<<<<<<< HEAD - cut->AddCut(VarManager::kDeltaR, 0.0, 5.0); - cut->AddCut(VarManager::kDitrackMass, 0.4, 0.58); - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kDeltaR, 0.0, 1.5); cut->AddCut(VarManager::kDitrackMass, 0.4, 0.58); cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) return cut; } @@ -6881,23 +6850,14 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) cut->AddCut(VarManager::kQ, -0.4, 0.2); cut->AddCut(VarManager::kDeltaR, 0.0, 1.0); cut->AddCut(VarManager::kDitrackMass, 0.4, 0.58); -<<<<<<< HEAD - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); - cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); cut->AddCut(VarManager::kVertexingQuadProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) cut->AddCut(VarManager::kPtOverPairPt, 1.0, 1.5); return cut; } if (!nameStr.compare("pairX3872_hasVtx")) { -<<<<<<< HEAD - cut->AddCut(VarManager::kVertexingProcCode, 0.5, 1.5); -======= cut->AddCut(VarManager::kVertexingProcCode, 0.5, 2.5); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) cut->AddCut(VarManager::kVertexingQuadProcCode,0.5, 1.5); return cut; } @@ -6941,20 +6901,15 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) } if (!nameStr.compare("pair_b0test")) { -<<<<<<< HEAD - cut->AddCut(VarManager::kDitrackMass, 0.49, 0.502); -======= cut->AddCut(VarManager::kDitrackMass, 0.49, 0.51); return cut; } if (!nameStr.compare("pair_bstest")) { cut->AddCut(VarManager::kDitrackMass, 1.0, 1.04); ->>>>>>> 450e56187 (updated psi2s analysis and dilepton v0) return cut; } ->>>>>>> 9942c0482 (updated psi2s analysis and dilepton v0) if (!nameStr.compare("pairPtLow1")) { cut->AddCut(VarManager::kPt, 2.0, 1000.0); return cut; From 892fb2ec08f7aae84c38de52f9b6ee233e2c80bf Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 12 Oct 2025 09:20:41 +0200 Subject: [PATCH 3/8] hist --- PWGDQ/Core/HistogramsLibrary.cxx | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index 64d38935f2c..08536c46030 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -969,6 +969,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "hDiTrackMass", "", false, 300, 0.0, 3.0, VarManager::kDitrackMass); hm->AddHistogram(histClass, "hMCPt_MCRap", "", false, 200, 0.0, 20.0, VarManager::kMCPt, 100, -2.0, 2.0, VarManager::kMCY); hm->AddHistogram(histClass, "hMCPhi", "", false, 100, -TMath::Pi(), TMath::Pi(), VarManager::kMCPhi); + hm->AddHistogram(histClass, "hMCPairPt", "", false, 200, 0.0, 20.0, VarManager::kPairPt); } if (!groupStr.CompareTo("mctruth_track")) { hm->AddHistogram(histClass, "PtMC", "MC pT", false, 200, 0.0, 20.0, VarManager::kMCPt); @@ -1097,6 +1098,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "Lz", "", false, 1000, -1.0, 1.0, VarManager::kVertexingLz); hm->AddHistogram(histClass, "Lxy", "", false, 1000, -1.0, 1.0, VarManager::kVertexingLxy); hm->AddHistogram(histClass, "Lxyz", "", false, 1000, -1.0, 1.0, VarManager::kVertexingLxyz); + hm->AddHistogram(histClass, "Lxyz_Mass", "", false, 50, 2.0, 4.0, VarManager::kMass, 1000, -1.0, 1.0, VarManager::kVertexingLxyz); hm->AddHistogram(histClass, "Tauz", "", false, 1000, -0.02, 0.02, VarManager::kVertexingTauz); hm->AddHistogram(histClass, "Tauxy", "", false, 1000, -0.03, 0.03, VarManager::kVertexingTauxy); hm->AddHistogram(histClass, "Tauxy_Mass_Pt", "", false, 50, 2.0, 4.0, VarManager::kMass, 10, 0.0, 20.0, VarManager::kPt, 1000, -0.03, 0.03, VarManager::kVertexingTauxy); @@ -1963,6 +1965,8 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h // hm->AddHistogram(histClass, "hPtDilepton_PtDihadron", "", false, 150, 0, 15.0, VarManager::kPairPt, 100, 0, 10, VarManager::kDitrackPt); hm->AddHistogram(histClass, "hPtDilepton_PtX3872", "", false, 150, 0, 15.0, VarManager::kPairPt, 150, 0.0, 15.0, VarManager::kQuadPt); hm->AddHistogram(histClass, "hQ_X3872", "", false, 300, -3.0, 3.0, VarManager::kQ); + hm->AddHistogram(histClass, "hQ_Pt_X3872", "", false, 300, -3.0, 3.0, VarManager::kQ, 150, 0.0, 15.0, VarManager::kQuadPt); + hm->AddHistogram(histClass, "hQ_hMass_defaultDileptonMass_X3872", "", false, 300, -3.0, 3.0, VarManager::kQ, 1000, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass); hm->AddHistogram(histClass, "hDeltaR1_X3872", "", false, 100, 0.0, 10.0, VarManager::kDeltaR1); hm->AddHistogram(histClass, "hDeltaR2_X3872", "", false, 100, 0.0, 10.0, VarManager::kDeltaR2); hm->AddHistogram(histClass, "hDeltaR_X3872", "", false, 100, 0.0, 10.0, VarManager::kDeltaR); @@ -1981,18 +1985,18 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "hMass_DeltaR_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kDeltaR); // hm->AddHistogram(histClass, "hMass_X3872_MassDihadron", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 150, 0.0, 3.0, VarManager::kDitrackMass); hm->AddHistogram(histClass, "hMass_defaultDileptonMass_X3872_MassDihadron", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 150, 0.0, 3.0, VarManager::kDitrackMass); - hm->AddHistogram(histClass, "hRap_X3872", "", false, 1000, 0.0, 5.0, VarManager::kRap); - hm->AddHistogram(histClass, "hMass_Rap_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 1000, 0.0, 5.0, VarManager::kRap); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_Rap_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 1000, 0.0, 5.0, VarManager::kRap); - hm->AddHistogram(histClass, "hDCAxyTrack1", "", false, 100, -0.1, 0.1, VarManager::kTrackDCAxy); - hm->AddHistogram(histClass, "hDCAzTrack1", "", false, 100, -0.1, 0.1, VarManager::kTrackDCAz); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_DCAxyTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 100, -0.1, 0.1, VarManager::kTrackDCAxy); - hm->AddHistogram(histClass, "hMass_defaultDileptonMass_DCAzTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 100, -0.1, 0.1, VarManager::kTrackDCAz); - hm->AddHistogram(histClass, "hMass_DCAxyTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, -0.1, 0.1, VarManager::kTrackDCAxy); - hm->AddHistogram(histClass, "hMass_DCAzTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, -0.1, 0.1, VarManager::kTrackDCAz); - hm->AddHistogram(histClass, "hPtTrack1", "", false, 100, 0.0, 10.0, VarManager::kPt); + // hm->AddHistogram(histClass, "hRap_X3872", "", false, 1000, 0.0, 5.0, VarManager::kRap); + // hm->AddHistogram(histClass, "hMass_Rap_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 1000, 0.0, 5.0, VarManager::kRap); + // hm->AddHistogram(histClass, "hMass_defaultDileptonMass_Rap_X3872", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 1000, 0.0, 5.0, VarManager::kRap); + // hm->AddHistogram(histClass, "hDCAxyTrack1", "", false, 100, -0.1, 0.1, VarManager::kTrackDCAxy); + // hm->AddHistogram(histClass, "hDCAzTrack1", "", false, 100, -0.1, 0.1, VarManager::kTrackDCAz); + // hm->AddHistogram(histClass, "hMass_defaultDileptonMass_DCAxyTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 100, -0.1, 0.1, VarManager::kTrackDCAxy); + // hm->AddHistogram(histClass, "hMass_defaultDileptonMass_DCAzTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 100, -0.1, 0.1, VarManager::kTrackDCAz); + // hm->AddHistogram(histClass, "hMass_DCAxyTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, -0.1, 0.1, VarManager::kTrackDCAxy); + // hm->AddHistogram(histClass, "hMass_DCAzTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, -0.1, 0.1, VarManager::kTrackDCAz); + // hm->AddHistogram(histClass, "hPtTrack1", "", false, 100, 0.0, 10.0, VarManager::kPt); hm->AddHistogram(histClass, "hMass_defaultDileptonMass_PtTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass, 100, 0.0, 10.0, VarManager::kPt); - hm->AddHistogram(histClass, "hMass_PtTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kPt); + // hm->AddHistogram(histClass, "hMass_PtTrack1", "", false, 100, 3.0, 5.0, VarManager::kQuadMass, 100, 0.0, 10.0, VarManager::kPt); } if (subGroupStr.Contains("vertexing")) { hm->AddHistogram(histClass, "UsedKF", "", false, 2, -0.5, 1.5, VarManager::kUsedKF); From 6a0a3254a65e033578a0227825ed322008d560b1 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 12 Oct 2025 09:23:04 +0200 Subject: [PATCH 4/8] added vars --- PWGDQ/Core/VarManager.cxx | 8 ++++++++ PWGDQ/Core/VarManager.h | 20 ++++++++++++++++++-- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 4a5769e0560..cf89f9d29e0 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -818,6 +818,8 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kVertexingTauxyErr] = "ns"; fgVariableNames[kVertexingProcCode] = "DCAFitterN<2> processing code"; fgVariableUnits[kVertexingProcCode] = ""; + fgVariableNames[kVertexingQuadProcCode] = "DCAFitterN<4> processing code dilepton"; + fgVariableUnits[kVertexingQuadProcCode] = ""; fgVariableNames[kVertexingChi2PCA] = "Pair #chi^{2} at PCA"; fgVariableUnits[kVertexingChi2PCA] = ""; fgVariableNames[kVertexingLxyOverErr] = "Pair Lxy/DLxy"; @@ -1254,6 +1256,10 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kDeltaR1] = ""; fgVariableNames[kDeltaR2] = "angular distance prong 2"; fgVariableUnits[kDeltaR2] = ""; + fgVariableNames[kDeltaR] = "angular distance"; + fgVariableUnits[kDeltaR] = ""; + fgVariableNames[kPtOverPairPt] = "p_{T}^{Quad}/p_{T}^{dilepton}"; + fgVariableUnits[kPtOverPairPt] = ""; fgVariableNames[kV22m] = "v_{2}(2)_{#mu^{-}}"; fgVariableUnits[kV22m] = ""; fgVariableNames[kV24m] = "v_{2}(4)_{#mu^{-}}"; @@ -1730,6 +1736,7 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kVertexingPz"] = kVertexingPz; fgVarNamesMap["kVertexingSV"] = kVertexingSV; fgVarNamesMap["kVertexingProcCode"] = kVertexingProcCode; + fgVarNamesMap["kVertexingQuadProcCode"] = kVertexingQuadProcCode; fgVarNamesMap["kVertexingChi2PCA"] = kVertexingChi2PCA; fgVarNamesMap["kCosThetaHE"] = kCosThetaHE; fgVarNamesMap["kPhiHE"] = kPhiHE; @@ -1897,6 +1904,7 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kQ"] = kQ; fgVarNamesMap["kDeltaR1"] = kDeltaR1; fgVarNamesMap["kDeltaR2"] = kDeltaR2; + fgVarNamesMap["kPtOverPairPt"] = kPtOverPairPt; fgVarNamesMap["kMassCharmHadron"] = kMassCharmHadron; fgVarNamesMap["kPtCharmHadron"] = kPtCharmHadron; fgVarNamesMap["kRapCharmHadron"] = kRapCharmHadron; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 63261d3d095..ce2fb918cd5 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -656,6 +656,8 @@ class VarManager : public TObject kVertexingTauxyzProjected, kVertexingTauz, kVertexingTauzErr, + kVertexingTauxyz, + kVertexingTauxyzErr, kVertexingPz, kVertexingSV, kVertexingProcCode, @@ -839,6 +841,8 @@ class VarManager : public TObject kDeltaR1, kDeltaR2, kDeltaR, + kPtOverPairPt, + kVertexingQuadProcCode, // DQ-HF correlation variables kMassCharmHadron, @@ -3692,6 +3696,10 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, m1 = o2::constants::physics::MassMuon; m2 = o2::constants::physics::MassMuon; } + if constexpr (pairType == kDecayToPiPi) { + m1 = o2::constants::physics::MassPionCharged; + m2 = o2::constants::physics::MassPionCharged; + } ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; @@ -3896,7 +3904,8 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, if (fabs(values[kVertexingLxyz]) < 1.e-8f) values[kVertexingLxyz] = 1.e-8f; values[kVertexingLxyzErr] = values[kVertexingLxyzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyzErr]) / values[kVertexingLxyz]; - values[kVertexingTauxy] = KFGeoTwoProng.GetPseudoProperDecayTime(KFPV, KFGeoTwoProng.GetMass()) / (o2::constants::physics::LightSpeedCm2NS); + // values[kVertexingTauxy] = KFGeoTwoProng.GetPseudoProperDecayTime(KFPV, KFGeoTwoProng.GetMass()) / (o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxy] = values[kVertexingLxy] * KFGeoTwoProng.GetMass() / (KFGeoTwoProng.GetPt() * o2::constants::physics::LightSpeedCm2NS); values[kVertexingTauz] = -1 * dzPair2PV * KFGeoTwoProng.GetMass() / (TMath::Abs(KFGeoTwoProng.GetPz()) * o2::constants::physics::LightSpeedCm2NS); values[kVertexingPz] = TMath::Abs(KFGeoTwoProng.GetPz()); values[kVertexingSV] = KFGeoTwoProng.GetZ(); @@ -5201,6 +5210,7 @@ void VarManager::FillDileptonTrackTrack(T1 const& dilepton, T2 const& hadron1, T ROOT::Math::PtEtaPhiMVector v23 = v2 + v3; values[kPairMass] = v1.M(); values[kPairPt] = v1.Pt(); + values[kPtOverPairPt] = values[kPairPt] > 0 ? values[kQuadPt] / values[kPairPt] : 0; values[kDitrackMass] = v23.M(); values[kDitrackPt] = v23.Pt(); values[kCosthetaDileptonDitrack] = (v1.Px() * v123.Px() + v1.Py() * v123.Py() + v1.Pz() * v123.Pz()) / (v1.P() * v123.P()); @@ -5230,7 +5240,7 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l float mtrack1, mtrack2; float mlepton1, mlepton2; - if constexpr (candidateType == kXtoJpsiPiPi || candidateType == kPsi2StoJpsiPiPi) { + if constexpr (candidateType == kXtoJpsiPiPi || candidateType == kPsi2StoJpsiPiPi || candidateType == kB0toJpsiEEPiPi) { mlepton1 = o2::constants::physics::MassElectron; mlepton2 = o2::constants::physics::MassElectron; mtrack1 = o2::constants::physics::MassPionCharged; @@ -5289,6 +5299,8 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l values[kVertexingTauxyErr] = -999.; values[kVertexingTauz] = -999.; values[kVertexingTauzErr] = -999.; + values[kVertexingTauxyz] = -999.; + values[kVertexingTauxyzErr] = -999.; values[kVertexingLzProjected] = -999.; values[kVertexingLxyProjected] = -999.; values[kVertexingLxyzProjected] = -999.; @@ -5349,6 +5361,9 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l values[kVertexingTauzErr] = values[kVertexingLzErr] * v1234.M() / (TMath::Abs(v1234.Pz()) * o2::constants::physics::LightSpeedCm2NS); values[kVertexingTauxyErr] = values[kVertexingLxyErr] * v1234.M() / (v1234.Pt() * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxyz] = values[kVertexingLxyz] * v1234.M() / (v1234.P() * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxyzErr] = values[kVertexingLxyzErr] * v1234.M() / (v1234.P() * o2::constants::physics::LightSpeedCm2NS); + values[kCosPointingAngle] = ((collision.posX() - secondaryVertex[0]) * v1234.Px() + (collision.posY() - secondaryVertex[1]) * v1234.Py() + (collision.posZ() - secondaryVertex[2]) * v1234.Pz()) / @@ -5524,6 +5539,7 @@ void VarManager::FillQuadMC(T1 const& dilepton, T2 const& track1, T2 const& trac values[kDeltaR] = sqrt(pow(values[kDeltaR1], 2) + pow(values[kDeltaR2], 2)); values[kDitrackMass] = v23.M(); values[kDitrackPt] = v23.Pt(); + values[kPairPt] = v1.Pt(); } //__________________________________________________________________ From f3a7a2c9efa66ee083894577199076ad87d1723b Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 12 Oct 2025 09:24:42 +0200 Subject: [PATCH 5/8] tasks --- PWGDQ/Tasks/dqEfficiency.cxx | 9 +++++- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 8 ++++- PWGDQ/Tasks/tableReader.cxx | 3 ++ PWGDQ/Tasks/tableReader_withAssoc.cxx | 43 +++++++++++++++++++++++--- 4 files changed, 56 insertions(+), 7 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency.cxx b/PWGDQ/Tasks/dqEfficiency.cxx index 3f97216f52d..de5cadfe2aa 100644 --- a/PWGDQ/Tasks/dqEfficiency.cxx +++ b/PWGDQ/Tasks/dqEfficiency.cxx @@ -1640,6 +1640,13 @@ struct AnalysisDileptonTrackTrack { groupedMCTracks.bindInternalIndicesTo(&tracksMC); runMCGen(groupedMCTracks); } + void processPsi2S(soa::Filtered::iterator const& event, MyBarrelTracksSelectedWithCov const& tracks, soa::Join const& dileptons, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + { + runDileptonTrackTrack(event, tracks, dileptons, eventsMC, tracksMC); + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&tracksMC); + runMCGen(groupedMCTracks); + } void processDummy(MyEvents&) { @@ -1723,7 +1730,7 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses) dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-hadron-mass"); } if (classStr.Contains("Quadruplet") || classStr.Contains("MCTruthRecQuad")) { - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-dihadron", "xtojpsipipi"); + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-dihadron", "xtojpsipipi,vertexing"); } } // end loop over histogram classes diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 22d5635c8eb..67dd4d04108 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -228,6 +228,7 @@ constexpr static uint32_t gkEventFillMapWithCovMultExtra = VarManager::ObjTypes: constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; constexpr static uint32_t gkTrackFillMapWithCovWithColl = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::ReducedTrackCollInfo; +constexpr static uint32_t gkTrackFillMapWithCovAmbiguities = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::AmbiTrack; // constexpr static uint32_t gkTrackFillMapWithColl = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::ReducedTrackCollInfo; constexpr static uint32_t gkMuonFillMap = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra; @@ -4637,7 +4638,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; } void DefineHistograms(HistogramManager* histMan, TString histClasses, const char* histGroups) @@ -4738,5 +4740,9 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char if (classStr.Contains("DileptonHadronCorrelation")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-hadron-correlation"); } + + if (classStr.Contains("Quadruplet")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-dihadron", histName); + } } // end loop over histogram classes } diff --git a/PWGDQ/Tasks/tableReader.cxx b/PWGDQ/Tasks/tableReader.cxx index e5fa216fc7a..72ce5e48df8 100644 --- a/PWGDQ/Tasks/tableReader.cxx +++ b/PWGDQ/Tasks/tableReader.cxx @@ -1514,6 +1514,8 @@ struct AnalysisSameEventPairing { constexpr bool trackHasCov = ((TTrackFillMap & VarManager::ObjTypes::TrackCov) > 0 || (TTrackFillMap & VarManager::ObjTypes::ReducedTrackBarrelCov) > 0); if constexpr ((TPairType == pairTypeEE) && trackHasCov && (TTwoProngFitter == true)) { dielectronExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauz], VarManager::fgValues[VarManager::kVertexingLz], VarManager::fgValues[VarManager::kVertexingLxy]); + } else { + dielectronExtraList(t1.globalIndex(), t2.globalIndex(), -999.f, -999.f, -999.f); } if constexpr ((TPairType == pairTypeMuMu) && (TTwoProngFitter == true)) { // LOGP(info, "mu1 collId = {}, mu2 collId = {}", t1.collisionId(), t2.collisionId()); @@ -2284,6 +2286,7 @@ struct AnalysisDileptonTrackTrack { PROCESS_SWITCH(AnalysisDileptonTrackTrack, processChicToJpsiPiPi, "Run dilepton-dihadron pairing to study X(3872), using skimmed data", false); PROCESS_SWITCH(AnalysisDileptonTrackTrack, processPsi2SToJpsiPiPi, "Run dilepton-dihadron pairing to study Psi(2S), using skimmed data", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processB0ToJpsiPiPi, "Run dilepton-dihadron pairing to study B0->Jpsi+Pi+Pi, using skimmed data", false); PROCESS_SWITCH(AnalysisDileptonTrackTrack, processDummy, "Dummy function", false); }; diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 996461765d6..d406ec123a5 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -259,7 +259,10 @@ constexpr static uint32_t gkEventFillMapWithCovWithColl = VarManager::ObjTypes:: // constexpr static uint32_t gkEventFillMapWithCovQvector = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended | VarManager::ObjTypes::ReducedEventVtxCov | VarManager::ObjTypes::ReducedEventQvector; constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; +constexpr static uint32_t gkTrackFillMapWithCovV0 = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::TrackV0Bits; constexpr static uint32_t gkTrackFillMapWithCovWithColl = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::ReducedTrackCollInfo; +constexpr static uint32_t gkTrackFillMapWithCovAmbiguities = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::AmbiTrack; +constexpr static uint32_t gkTrackFillMapWithCovV0Ambiguities = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::TrackV0Bits | VarManager::ObjTypes::AmbiTrack; // constexpr static uint32_t gkMuonFillMap = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra; constexpr static uint32_t gkMuonFillMapWithCov = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra | VarManager::ObjTypes::ReducedMuonCov; @@ -269,6 +272,7 @@ constexpr static uint32_t gkV0FillMap = VarManager::ObjTypes::TrackV0Bits; constexpr static int pairTypeEE = VarManager::kDecayToEE; constexpr static int pairTypeMuMu = VarManager::kDecayToMuMu; +constexpr static int pairTypePiPi = VarManager::kDecayToPiPi; // constexpr static int pairTypeEMu = VarManager::kElectronMuon; // Global function used to define needed histogram classes @@ -1758,7 +1762,7 @@ struct AnalysisSameEventPairing { bool isFirst = true; for (auto& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi) { 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 @@ -1996,7 +2000,7 @@ struct AnalysisSameEventPairing { 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 constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi) { if (isLeg1Ambi && isLeg2Ambi) { std::pair iPair(a1.reducedtrackId(), a2.reducedtrackId()); if (fAmbiguousPairs.find(iPair) != fAmbiguousPairs.end()) { @@ -2024,7 +2028,7 @@ struct AnalysisSameEventPairing { fHistMan->FillHistClass(histNames[icut][3 + histIdxOffset + 6].Data(), VarManager::fgValues); } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi) { fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s", fTrackCuts[icut].Data()), VarManager::fgValues); if (isAmbiExtra) { fHistMan->FillHistClass(Form("PairsBarrelSEPM_ambiguousextra_%s", fTrackCuts[icut].Data()), VarManager::fgValues); @@ -2044,7 +2048,7 @@ struct AnalysisSameEventPairing { fHistMan->FillHistClass(histNames[icut][4 + histIdxOffset + 6].Data(), VarManager::fgValues); } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi) { fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s", fTrackCuts[icut].Data()), VarManager::fgValues); if (isAmbiExtra) { fHistMan->FillHistClass(Form("PairsBarrelSEPP_ambiguousextra_%s", fTrackCuts[icut].Data()), VarManager::fgValues); @@ -2063,7 +2067,7 @@ struct AnalysisSameEventPairing { fHistMan->FillHistClass(histNames[icut][5 + histIdxOffset + 6].Data(), VarManager::fgValues); } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi) { fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s", fTrackCuts[icut].Data()), VarManager::fgValues); if (isAmbiExtra) { fHistMan->FillHistClass(Form("PairsBarrelSEMM_ambiguousextra_%s", fTrackCuts[icut].Data()), VarManager::fgValues); @@ -2299,6 +2303,14 @@ struct AnalysisSameEventPairing { runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); } + // void processPiPiSkimmed + void processPiPiSkimmedNoCov(MyEventsSelected const& events, + soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) + { + runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); + } + void processBarrelOnlySkimmedNoCovWithMultExtra(MyEventsMultExtraSelected const& events, soa::Join const& barrelAssocs, MyBarrelTracksWithAmbiguities const& barrelTracks) @@ -2380,6 +2392,7 @@ struct AnalysisSameEventPairing { PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlySkimmed, "Run barrel only pairing, with skimmed tracks", false); PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyWithCollSkimmed, "Run barrel only pairing, with skimmed tracks and with collision information", false); PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlySkimmedNoCov, "Run barrel only pairing (no covariances), with skimmed tracks and with collision information", false); + PROCESS_SWITCH(AnalysisSameEventPairing, processPiPiSkimmedNoCov, "Run barrel only pairing (no covariances), with skimmed tracks, with PiPi decay channel", false); PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlySkimmedNoCovWithMultExtra, "Run barrel only pairing (no covariances), with skimmed tracks, with collision information, with MultsExtra", false); PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyWithQvectorCentrSkimmedNoCov, "Run barrel only pairing (no covariances), with skimmed tracks, with Qvector from central framework", false); PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlySkimmedFlow, "Run barrel only pairing, with skimmed tracks and with flow", false); @@ -4268,6 +4281,19 @@ struct AnalysisDileptonTrackTrack { VarManager::FillTrack(dilepton, fValuesQuadruplet); // LOGP(info, "is dilepton selected: {}", fDileptonCut.IsSelected(fValuesQuadruplet)); + bool isAmbiguousLepton = (dilepton.filterMap_raw() & (static_cast(1) << 28)) || + (dilepton.filterMap_raw() & (static_cast(1) << 29)) || + (dilepton.filterMap_raw() & (static_cast(1) << 30)) || + (dilepton.filterMap_raw() & (static_cast(1) << 31)); + // if (isAmbi && isAmbiguousLepton) + // continue; // skip ambiguous dileptons + if constexpr ((TTrackFillMap & VarManager::ObjTypes::AmbiTrack)> 0) { + if (isAmbiguousLepton) { + // LOGP(info, "Dilepton {} is ambiguous, skipping", dilepton.index0Id()); + continue; // skip ambiguous dileptons + } + } + // apply the dilepton cut if (!fDileptonCut.IsSelected(fValuesQuadruplet)) continue; @@ -4295,6 +4321,13 @@ struct AnalysisDileptonTrackTrack { continue; } + if constexpr ((TTrackFillMap & VarManager::ObjTypes::AmbiTrack)> 0) { + if ((track1.barrelAmbiguityInBunch() > 1) || (track2.barrelAmbiguityInBunch() > 1) || + (track1.barrelAmbiguityOutOfBunch() > 1) || (track2.barrelAmbiguityOutOfBunch() > 1)) { + continue; + } + } + // fill variables VarManager::FillDileptonTrackTrack(dilepton, track1, track2, fValuesQuadruplet); if constexpr (NProng > 0) { From 28ce7547ec8bb199da0472c08aa042c6cf15ea03 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 12 Oct 2025 09:58:12 +0200 Subject: [PATCH 6/8] remove dilepton k0s in dilepton-track in varmanager.h --- PWGDQ/Core/VarManager.h | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index ce2fb918cd5..7125a18e102 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -4306,7 +4306,7 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton procCode = VarManager::fgFitterThreeProngFwd.process(pars1, pars2, pars3); procCodeJpsi = VarManager::fgFitterTwoProngFwd.process(pars1, pars2); - } else if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi || candidateType == kB0toJpsiEEPiPi) && trackHasCov) { + } else if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi) && trackHasCov) { if constexpr ((candidateType == kBtoJpsiEEK) && trackHasCov) { mlepton1 = o2::constants::physics::MassElectron; mlepton2 = o2::constants::physics::MassElectron; @@ -4315,11 +4315,7 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton mlepton1 = o2::constants::physics::MassKaonCharged; mlepton2 = o2::constants::physics::MassPionCharged; mtrack = o2::constants::physics::MassPionCharged; - } else if constexpr ((candidateType) == kB0toJpsiEEPiPi && trackHasCov) { - mlepton1 = o2::constants::physics::MassElectron; - mlepton2 = o2::constants::physics::MassElectron; - mtrack = o2::constants::physics::MassKaonNeutral; - } + } std::array lepton1pars = {lepton1.y(), lepton1.z(), lepton1.snp(), lepton1.tgl(), lepton1.signed1Pt()}; std::array lepton1covs = {lepton1.cYY(), lepton1.cZY(), lepton1.cZZ(), lepton1.cSnpY(), lepton1.cSnpZ(), lepton1.cSnpSnp(), lepton1.cTglY(), lepton1.cTglZ(), lepton1.cTglSnp(), lepton1.cTglTgl(), @@ -4393,7 +4389,7 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton o2::dataformats::VertexBase primaryVertex = {std::move(vtxXYZ), std::move(vtxCov)}; auto covMatrixPV = primaryVertex.getCov(); - if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi || candidateType == kB0toJpsiEEPiPi) && trackHasCov) { + if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi) && trackHasCov) { secondaryVertex = fgFitterThreeProngBarrel.getPCACandidate(); covMatrixPCA = fgFitterThreeProngBarrel.calcPCACovMatrixFlat(); } else if constexpr (candidateType == kBcToThreeMuons && muonHasCov) { From 1517c4bbd3e1a4209cc523c4a5d6d95b14951963 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 12 Oct 2025 10:00:41 +0200 Subject: [PATCH 7/8] =?UTF-8?q?=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- PWGDQ/Core/VarManager.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 7125a18e102..bc10126d179 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -4315,7 +4315,7 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton mlepton1 = o2::constants::physics::MassKaonCharged; mlepton2 = o2::constants::physics::MassPionCharged; mtrack = o2::constants::physics::MassPionCharged; - } + } std::array lepton1pars = {lepton1.y(), lepton1.z(), lepton1.snp(), lepton1.tgl(), lepton1.signed1Pt()}; std::array lepton1covs = {lepton1.cYY(), lepton1.cZY(), lepton1.cZZ(), lepton1.cSnpY(), lepton1.cSnpZ(), lepton1.cSnpSnp(), lepton1.cTglY(), lepton1.cTglZ(), lepton1.cTglSnp(), lepton1.cTglTgl(), From a4a12b2d7aa2a7c63d0662732cd9126982e6b517 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sun, 19 Oct 2025 19:11:53 +0200 Subject: [PATCH 8/8] tmp update --- PWGDQ/Core/CutsLibrary.cxx | 50 +++++++++++++ PWGDQ/Core/MCSignalLibrary.cxx | 78 ++++++++++++++++++++ PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 33 ++++++++- PWGDQ/Tasks/tableReader_withAssoc.cxx | 98 +++++++++++++++++++++----- 4 files changed, 242 insertions(+), 17 deletions(-) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index 7691e58a20c..673c64bd660 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -3655,6 +3655,24 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + for (unsigned int i = 0; i < 10; ++i) { + for (unsigned int j = 0; j < 5; ++j) { + if (!nameStr.compare(Form("electronSelection_cutvar_yiping%i%i", i, j))) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); + cut->AddCut(GetAnalysisCut("dcaCut1_ionut")); + cut->AddCut(GetAnalysisCut(Form("pidJpsiEle_yiping_%i", i))); + cut->AddCut(GetAnalysisCut(Form("pidJpsiEle_yiping_hadronRej_%i", j))); + return cut; + } + } + } + + if (!nameStr.compare("trackingEff_kine_ele")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine")); + return cut; + } + if (!nameStr.compare("trackingEff_ele")) { cut->AddCut(GetAnalysisCut("jpsiStandardKine")); cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug4")); @@ -4825,6 +4843,20 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("electronStandardQualityForO2MCdebug5")) { + cut->AddCut(VarManager::kIsITSibAny, 0.5, 1.5); + cut->AddCut(VarManager::kTPCchi2, 0.0, 4.0); + cut->AddCut(VarManager::kTPCncls, 100.0, 161.); + return cut; + } + + if (!nameStr.compare("electronStandardQualityForO2MCdebug6")) { + cut->AddCut(VarManager::kIsITSibAny, 0.5, 10); + cut->AddCut(VarManager::kTPCchi2, 0.0, 4.0); + cut->AddCut(VarManager::kTPCncls, 120.0, 161.); + return cut; + } + if (!nameStr.compare("electronStandardQualityITSOnly")) { cut->AddCut(VarManager::kIsSPDany, 0.5, 1.5); cut->AddCut(VarManager::kITSchi2, 0.0, 5.0); @@ -5281,6 +5313,24 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + // yiping's test: cut variation study for jpsi electron PID + std::vector yiping_TPCnSigmaEl_low = {-4.0, -3.0, -2.7, -2.5, -2.0, -4.0, -3.0, -2.7, -2.5, -2.0}; + std::vector yiping_TPCnSigmaEl_high = {4.0, 4.0, 4.0, 4.0, 4.0, 3.0, 3.0, 3.0, 3.0, 3.0}; + for (unsigned int i = 0; i < yiping_TPCnSigmaEl_low.size(); i++) { + if (!nameStr.compare(Form("pidJpsiEle_yiping_%i", i))) { + cut->AddCut(VarManager::kTPCnSigmaEl, yiping_TPCnSigmaEl_low.at(i), yiping_TPCnSigmaEl_high.at(i)); + return cut; + } + } + std::vector yiping_hadronRej = {2.0, 2.5, 2.7, 3.0, 3.5}; + for (unsigned int i = 0; i < yiping_hadronRej.size(); i++) { + if (!nameStr.compare(Form("pidJpsiEle_yiping_hadronRej_%i", i))) { + cut->AddCut(VarManager::kTPCnSigmaPi, yiping_hadronRej.at(i), 999.0); + cut->AddCut(VarManager::kTPCnSigmaPr, yiping_hadronRej.at(i), 999.0); + return cut; + } + } + if (!nameStr.compare("pidCut_lowP_Corr")) { cut->AddCut(VarManager::kTPCnSigmaEl_Corr, -3.0, 3.0, false, VarManager::kP, 0.0, 5.0); cut->AddCut(VarManager::kTPCnSigmaPi_Corr, 3.0, 999, false, VarManager::kP, 0.0, 5.0); diff --git a/PWGDQ/Core/MCSignalLibrary.cxx b/PWGDQ/Core/MCSignalLibrary.cxx index d32054afb5a..2f199d73d98 100644 --- a/PWGDQ/Core/MCSignalLibrary.cxx +++ b/PWGDQ/Core/MCSignalLibrary.cxx @@ -1834,18 +1834,55 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) return signal; } + if (!nameStr.compare("JpsiFromPromptPsi2S")) { + MCProng prong(2,{443, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {true}); + signal = new MCSignal(name, "Jpsi from prompt Psi2S", {prong}, {-1}); + return signal; + } + + if (!nameStr.compare("JpsiFromNonpromptPsi2S")) { + MCProng prong(2,{443, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {false}); + signal = new MCSignal(name, "Jpsi from non-prompt Psi2S", {prong}, {-1}); + return signal; + } + + if (!nameStr.compare("eFromJpsiFromPsi2S")) { MCProng prong(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); signal = new MCSignal(name, "Electron from Jpsi from Psi2S", {prong}, {1}); return signal; } + if (!nameStr.compare("eFromJpsiFromPromptPsi2S")) { + MCProng prong(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}, false, {503}, {true}); + signal = new MCSignal(name, "Electron from Jpsi from prompt Psi2S", {prong}, {1}); + return signal; + } + + if (!nameStr.compare("eFromJpsiFromNonpromptPsi2S")) { + MCProng prong(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}, false, {503}, {false}); + signal = new MCSignal(name, "Electron from Jpsi from non-prompt Psi2S", {prong}, {1}); + return signal; + } + if (!nameStr.compare("PionFromPsi2S")) { MCProng prong(1, {211}, {true}, {false}, {0}, {0}, {false}, false, {100443}, {false}); signal = new MCSignal(name, "Pion from Jpsi from Psi2S", {prong}, {-1}); return signal; } + if (!nameStr.compare("PionFromPromptPsi2S")) { + MCProng prong(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {true}); + signal = new MCSignal(name, "Pion from prompt Psi2S", {prong}, {-1}); + return signal; + } + + if (!nameStr.compare("PionFromNonpromptPsi2S")) { + MCProng prong(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {false}); + signal = new MCSignal(name, "Pion from non-prompt Psi2S", {prong}, {-1}); + return signal; + } + if (!nameStr.compare("eeFromJpsiFromX3872")) { MCProng prong(2, {11, 443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {9920443}, {false}); signal = new MCSignal(name, "Electron pair from Jpsi from X3872", {prong, prong}, {1, 1}); @@ -1872,6 +1909,18 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) return signal; } + if (!nameStr.compare("eeFromJpsiFromPromptPsi2S")) { + MCProng prong(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}, false, {503}, {true}); + signal = new MCSignal(name, "Electron pair from Jpsi from prompt Psi2S", {prong, prong}, {1, 1}); + return signal; + } + + if (!nameStr.compare("eeFromJpsiFromNonpromptPsi2S")) { + MCProng prong(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}, false, {503}, {false}); + signal = new MCSignal(name, "Electron pair from Jpsi from non-prompt Psi2S", {prong, prong}, {1, 1}); + return signal; + } + if (!nameStr.compare("JpsiPiPiFromPsi2S")) { MCProng prongJpsi(2, {443, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); MCProng prongPi(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); @@ -1879,12 +1928,41 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) return signal; } + if (!nameStr.compare("JpsiPiPiFromPromptPsi2S")) { + MCProng prongJpsi(2, {443, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {true}); + MCProng prongPi(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {true}); + signal = new MCSignal(name, "Jpsi and pion pair from prompt Psi2S", {prongJpsi, prongPi, prongPi}, {1, 1, 1}); + return signal; + } + + if (!nameStr.compare("JpsiPiPiFromNonpromptPsi2S")) { + MCProng prongJpsi(2, {443, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {false}); + MCProng prongPi(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {false}); + signal = new MCSignal(name, "Jpsi and pion pair from non-prompt Psi2S", {prongJpsi, prongPi, prongPi}, {1, 1, 1}); + return signal; + } + if (!nameStr.compare("eePiPiFromPsi2S")) { MCProng pronge(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); MCProng prongPi(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); signal = new MCSignal(name, "Electron pair and pion pair from Psi2S", {pronge, pronge, prongPi, prongPi}, {2, 2, 1, 1}); return signal; } + + if (!nameStr.compare("eePiPiFromPromptPsi2S")) { + MCProng pronge(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}, false, {503}, {true}); + MCProng prongPi(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {true}); + signal = new MCSignal(name, "Electron pair and pion pair from prompt Psi2S", {pronge, pronge, prongPi, prongPi}, {2, 2, 1, 1}); + return signal; + } + + if (!nameStr.compare("eePiPiFromNonpromptPsi2S")) { + MCProng pronge(3, {11, 443, 100443}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}, false, {503}, {false}); + MCProng prongPi(2, {211, 100443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {false}); + signal = new MCSignal(name, "Electron pair and pion pair from non-prompt Psi2S", {pronge, pronge, prongPi, prongPi}, {2, 2, 1, 1}); + return signal; + } + return nullptr; } diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 67dd4d04108..7d53745303f 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -1341,6 +1341,8 @@ struct AnalysisSameEventPairing { Configurable recSignals{"cfgBarrelMCRecSignals", "", "Comma separated list of MC signals (reconstructed)"}; Configurable recSignalsJSON{"cfgMCRecSignalsJSON", "", "Comma separated list of MC signals (reconstructed) via JSON"}; Configurable skimSignalOnly{"cfgSkimSignalOnly", false, "Configurable to select only matched candidates"}; + Configurable genPtMin{"cfgGenPtMin", 1.0f, "Min pt cut for generated particles"}; + Configurable genEtaMax{"cfgGenEtaMax", 0.9f, "Max eta cut for generated particles"}; // Configurable runMCGenPair{"cfgRunMCGenPair", false, "Do pairing of true MC particles"}; } fConfigMC; @@ -1655,6 +1657,7 @@ struct AnalysisSameEventPairing { } else if (sig->GetNProngs() == 2) { histNames += Form("MCTruthGenPair_%s;", sig->GetName()); histNames += Form("MCTruthGenPairSel_%s;", sig->GetName()); + histNames += Form("MCTruthGenPairAccepted_%s;", sig->GetName()); fHasTwoProngGenMCsignals = true; } } @@ -2155,6 +2158,11 @@ struct AnalysisSameEventPairing { // WARNING! To be checked dileptonMiniTreeGen(mcDecision, -999, t1.pt(), t1.eta(), t1.phi(), t2.pt(), t2.eta(), t2.phi()); } + if (t1.pt() < fConfigMC.genPtMin) continue; + if (TMath::Abs(t1.eta()) > fConfigMC.genEtaMax) continue; + if (t2.pt() < fConfigMC.genPtMin) continue; + if (TMath::Abs(t2.eta()) > fConfigMC.genEtaMax) continue; + fHistMan->FillHistClass(Form("MCTruthGenPairAccepted_%s", sig->GetName()), VarManager::fgValues); } isig++; } @@ -4168,6 +4176,10 @@ struct AnalysisDileptonTrackTrack{ Configurable fConfigMCRecSignals{"cfgBarrelMCRecSignals", "", "Comma separated list of MC signals (reconstructed)"}; Configurable fConfigMCGenSignals{"cfgBarrelMCGenSignals", "", "Comma separated list of MC signals (generated)"}; Configurable fConfigDileptonMCRecSignal{"cfgDileptonMCRecSignal", "", "Comma separated list of MC signals (reconstructed)"}; + Configurable fConfigTrackPtMin{"cfgTrackPtMin", 0.15f, "Minimum pt of tracks to be used in the quadruplet"}; + Configurable fConfigTrackEtaAbs{"cfgTrackEtaAbs", 0.9f, "Maximum absolute eta of tracks to be used in the quadruplet"}; + Configurable fConfigDileptonPtMin{"cfgDileptonPtMin", 1.0f, "Minimum pt of dileptons to be used in the quadruplet"}; + Configurable fConfigDileptonEtaAbs{"cfgDileptonEtaAbs", 0.9f, "Maximum absolute eta of dileptons to be used in the quadruplet"}; Configurable fConfigAddDileptonHistogram{"cfgAddDileptonHistogram", "barrel", "Comma separated list of histograms"}; Configurable fConfigAddQuadrupletHistogram{"cfgAddQuadrupletHistogram", "xtojpsipipi", "Comma separated list of histograms"}; @@ -4296,7 +4308,7 @@ struct AnalysisDileptonTrackTrack{ if (isMCGen) { for (auto& sig : fGenMCSignals) { DefineHistograms(fHistMan, Form("MCTruthGenQuad_%s", sig->GetName()), ""); - DefineHistograms(fHistMan, Form("MCTruthGenQuadSel_%s", sig->GetName()), ""); + DefineHistograms(fHistMan, Form("MCTruthGenQuadAccepted_%s", sig->GetName()), ""); } // DefineHistograms(fHistMan, "MCTruthGenQuadAccepted", ""); } @@ -4570,7 +4582,26 @@ struct AnalysisDileptonTrackTrack{ auto dilepton = mcTracks.rawIteratorAt(daughterIdFirst); auto track1 = mcTracks.rawIteratorAt(daughterIdFirst + 1); auto track2 = mcTracks.rawIteratorAt(daughterIdFirst + 2); + // LOGP(info, "PDG of dilepton: {}, track1: {}, track2: {}", dilepton.pdgCode(), track1.pdgCode(), track2.pdgCode()); VarManager::FillQuadMC(dilepton, track1, track2); + int daughterIdFirst2 = dilepton.daughtersIds()[0]; + int daughterIdEnd2 = dilepton.daughtersIds()[1]; + int Ndaughters2 = daughterIdEnd2 - daughterIdFirst2 + 1; + if (Ndaughters2 == 2) { + auto lepton1 = mcTracks.rawIteratorAt(daughterIdFirst2); + auto lepton2 = mcTracks.rawIteratorAt(daughterIdFirst2 + 1); + VarManager::FillPairMC(lepton1, lepton2); + // LOGP(info, "PDG of lepton1: {}, lepton2: {}", lepton1.pdgCode(), lepton2.pdgCode()); + if (lepton1.pt() > fConfigDileptonPtMin && lepton2.pt() > fConfigDileptonPtMin) { + if (TMath::Abs(lepton1.eta()) < fConfigDileptonEtaAbs && TMath::Abs(lepton2.eta()) < fConfigDileptonEtaAbs) { + if (track1.pt() > fConfigTrackPtMin && track2.pt() > fConfigTrackPtMin) { + if (TMath::Abs(track1.eta()) < fConfigTrackEtaAbs && TMath::Abs(track2.eta()) < fConfigTrackEtaAbs) { + fHistMan->FillHistClass(Form("MCTruthGenQuadAccepted_%s", sig->GetName()), VarManager::fgValues); + } + } + } + } + } } fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues); } diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index d406ec123a5..528a7bf367f 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -239,6 +239,7 @@ using MyV0s = aod::V0Datas; using MyV0sWithCov = soa::Join; using FullTracksExt = soa::Join; +using FullTracksExtEl = soa::Join; // bit maps used for the Fill functions of the VarManager constexpr static uint32_t gkEventFillMap = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended; @@ -3792,10 +3793,10 @@ struct AnalysisDileptonTrack { struct AnalysisDileptonV0 { OutputObj fOutputList{"output"}; - Configurable fConfigDileptonLowMass{"cfgDileptonLowMass", 2.8, "Low mass cut for the dileptons used in analysis"}; - Configurable fConfigDileptonHighMass{"cfgDileptonHighMass", 3.2, "High mass cut for the dileptons used in analysis"}; - Configurable fConfigDileptonpTCut{"cfgDileptonpTCut", 0.0, "pT cut for dileptons used in the triplet vertexing"}; - Configurable fConfigDileptonLxyCut{"cfgDileptonLxyCut", 0.0, "Lxy cut for dileptons used in the triplet vertexing"}; + // Configurable fConfigDileptonLowMass{"cfgDileptonLowMass", 2.8, "Low mass cut for the dileptons used in analysis"}; + // Configurable fConfigDileptonHighMass{"cfgDileptonHighMass", 3.2, "High mass cut for the dileptons used in analysis"}; + // Configurable fConfigDileptonpTCut{"cfgDileptonpTCut", 0.0, "pT cut for dileptons used in the triplet vertexing"}; + // Configurable fConfigDileptonLxyCut{"cfgDileptonLxyCut", 0.0, "Lxy cut for dileptons used in the triplet vertexing"}; // Configurable fConfigV0Type{"cfgV0Type", -1, "V0 type"}; Configurable fConfigUseKFVertexing{"cfgUseKFVertexing", false, "Use KF Particle for secondary vertex reconstruction (DCAFitter is used by default)"}; @@ -3815,6 +3816,13 @@ struct AnalysisDileptonV0 { // Configurable DaughterEtaCut{"DaughterEtaCut", 0.9, "Eta cut for V0 daughter tracks"}; } fConfigV0Cuts; + struct : ConfigurableGroup { + Configurable cfgDileptonMassLow{"DileptonMassLow", 2.8, "Low mass cut for the dileptons used in analysis"}; + Configurable cfgDileptonMassHigh{"DileptonMassHigh", 3.2, "High mass cut for the dileptons used in analysis"}; + Configurable cfgDileptonpTCut{"DileptonpTCut", 0.0, "pT cut for dileptons used in the triplet vertexing"}; + Configurable cfgDileptonLxyCut{"DileptonLxyCut", 0.0, "Lxy cut for dileptons used in the triplet vertexing"}; + } fConfigDileptonCuts; + Configurable fConfigUseRemoteField{"cfgUseRemoteField", false, "Chose whether to fetch the magnetic field from ccdb or set it manually"}; Configurable fConfigGRPmagPath{"cfgGrpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable fConfigMagField{"cfgMagField", 5.0f, "Manually set magnetic field"}; @@ -3824,7 +3832,7 @@ struct AnalysisDileptonV0 { Configurable fConfigGeoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; Filter eventFilter = aod::dqanalysisflags::isEventSelected > static_cast(0); - Filter dileptonFilter = aod::reducedpair::sign == 0 && aod::reducedpair::mass > fConfigDileptonLowMass && aod::reducedpair::mass < fConfigDileptonHighMass; + Filter dileptonFilter = aod::reducedpair::sign == 0 && aod::reducedpair::mass > fConfigDileptonCuts.cfgDileptonMassLow && aod::reducedpair::mass < fConfigDileptonCuts.cfgDileptonMassHigh && aod::reducedpair::pt > fConfigDileptonCuts.cfgDileptonpTCut && aod::reducedpair::lxy > fConfigDileptonCuts.cfgDileptonLxyCut; // Filter filterV0Selected = aod::v0mapID::v0addid == fConfigV0Type; // Filter filterV0Selected = aod::v0data::mK0Short > fConfigV0Cuts.V0MassLow && aod::v0data::mK0Short < fConfigV0Cuts.V0MassHigh; @@ -3864,14 +3872,17 @@ struct AnalysisDileptonV0 { // fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); // --------------------------------------------- - registry.add("hInvMassDileptonV0", "hInvMassDileptonV0", HistType::kTH1F, {{2000, 4.0, 6.0}}); + registry.add("hMassB0", "hMassB0", HistType::kTH1F, {{2000, 4.0, 6.0}}); + registry.add("h3MassB0PtB0PtDilepton", "hMassB0PtB0PtDilepton", HistType::kTH3F, {{2000, 4.0, 6.0}, {200, 0.0, 10.0}, {200, 0.0, 10.0}}); + registry.add("hMassChic", "hMassChic", HistType::kTH1F, {{2000, 2.5, 4.5}}); registry.add("hInvMassDilepton", "hInvMassDilepton", HistType::kTH1F, {{2000, 2.0, 4.0}}); registry.add("hPtV0", "hPtV0", HistType::kTH1F, {{2000, 0.0, 10.0}}); registry.add("hPtV0Pos", "hPtV0Pos", HistType::kTH1F, {{2000, 0.0, 10.0}}); registry.add("hPosDCAToPV", "hPosDCAToPV", HistType::kTH1F, {{1000, -5.0, 5.0}}); registry.add("hNegDCAToPV", "hNegDCAToPV", HistType::kTH1F, {{1000, -5.0, 5.0}}); // registry.add("") - registry.add("hV0Selected", "hV0Selected", HistType::kTH1F, {{100, 0.45, 0.55}}); + registry.add("hMassK0s", "hMassK0s", HistType::kTH1F, {{100, 0.45, 0.55}}); + registry.add("hMassGamma", "hMassGamma", HistType::kTH1F, {{100, 0.0, 0.1}}); registry.add("hIsV0Vertexing", "hIsV0Vertexing", HistType::kTH1F, {{2, -0.5, 1.5}}); registry.add("hIsDileptonVertexing", "hIsDileptonVertexing", HistType::kTH1F, {{2, -0.5, 1.5}}); registry.add("hAssociatedCollisions", "hAssociatedCollisions", HistType::kTH1F, {{100, -0.5, 99.5}}); @@ -3890,7 +3901,7 @@ struct AnalysisDileptonV0 { if (event.collisionId() != dilepton.collisionId()) { continue; } - LOGP(info, "dilepton in collision: {}", dilepton.collisionId()); + // LOGP(info, "dilepton in collision: {}", dilepton.collisionId()); // get full track info of tracks based on the index int indexLepton1 = dilepton.index0Id(); int indexLepton2 = dilepton.index1Id(); @@ -3909,7 +3920,7 @@ struct AnalysisDileptonV0 { registry.fill(HIST("hInvMassDilepton"), dilepton.mass()); for (auto v0 : v0s) { - LOGP(info, "v0 in collision: {}", v0.collisionId()); + // LOGP(info, "v0 in collision: {}", v0.collisionId()); // apply basic V0 cuts registry.fill(HIST("hPosDCAToPV"), v0.dcapostopv()); registry.fill(HIST("hNegDCAToPV"), v0.dcanegtopv()); @@ -3935,19 +3946,20 @@ struct AnalysisDileptonV0 { } - float hadronMass = o2::constants::physics::MassKaonNeutral; + // float hadronMass = o2::constants::physics::MassKaonNeutral; + float hadronMass = v0.mK0Short(); float daughterMass = o2::constants::physics::MassPionCharged; // Check V0 mass window - ROOT::Math::PtEtaPhiMVector v0vec(v0_pos.pt(), v0_pos.eta(), v0_pos.phi(), daughterMass); - ROOT::Math::PtEtaPhiMVector v1vec(v0_neg.pt(), v0_neg.eta(), v0_neg.phi(), daughterMass); - ROOT::Math::PtEtaPhiMVector v0sum = v0vec + v1vec; + // ROOT::Math::PtEtaPhiMVector v0vec(v0_pos.pt(), v0_pos.eta(), v0_pos.phi(), daughterMass); + // ROOT::Math::PtEtaPhiMVector v1vec(v0_neg.pt(), v0_neg.eta(), v0_neg.phi(), daughterMass); + // ROOT::Math::PtEtaPhiMVector v0sum = v0vec + v1vec; // if (v0sum.M() < fConfigV0Cuts.V0MassLow || v0sum.M() > fConfigV0Cuts.V0MassHigh) { // continue; // } if (v0.mK0Short() < fConfigV0Cuts.V0MassLow || v0.mK0Short() > fConfigV0Cuts.V0MassHigh) { continue; } - registry.fill(HIST("hV0Selected"), v0.mK0Short()); + registry.fill(HIST("hMassK0s"), v0.mK0Short()); // VarManager::FillDileptonHadron(dilepton, v0, fValuesV0, hadronMass); //Reconstruct directly here @@ -3957,7 +3969,9 @@ struct AnalysisDileptonV0 { float invmass = v12.M(); // float massJpsi = 3.0969; float mass = invmass - dilepton.mass() + 3.096; - registry.fill(HIST("hInvMassDileptonV0"), mass); + float pt = v12.Pt(); + registry.fill(HIST("hMassB0"), mass); + registry.fill(HIST("h3MassB0PtB0PtDilepton"), mass, pt, dilepton.pt()); registry.fill(HIST("hPtV0"), v0.pt()); registry.fill(HIST("hPtV0Pos"), v0_pos.pt()); @@ -4008,6 +4022,41 @@ struct AnalysisDileptonV0 { auto covMatrixPV = primaryVertex.getCov(); // TODO: calculate decay length } // secondary vertex reconstruction + } else if constexpr (TCandidateType == VarManager::kChictoJpsiEE) { + auto v0_pos = v0.template posTrack_as(); + auto v0_neg = v0.template negTrack_as(); + + // apply daughter track cuts + if (v0_pos.tpcNClsFound() < fConfigV0Cuts.CutTPCncls || v0_neg.tpcNClsFound() < fConfigV0Cuts.CutTPCncls) { + continue; + } + if (std::abs(v0_pos.tpcNSigmaEl()) > fConfigV0Cuts.CutTPCnSigPi || std::abs(v0_neg.tpcNSigmaEl()) > fConfigV0Cuts.CutTPCnSigPi) { + continue; + } + if (v0_pos.tpcChi2NCl() > fConfigV0Cuts.CutTPCChi2Ncls || v0_neg.tpcChi2NCl() > fConfigV0Cuts.CutTPCChi2Ncls) { + continue; + } + + + // float hadronMass = o2::constants::physics::MassKaonNeutral; + float hadronMass = v0.mGamma(); + // Check V0 mass window + if (v0.mGamma() < fConfigV0Cuts.V0MassLow || v0.mGamma() > fConfigV0Cuts.V0MassHigh) { + continue; + } + registry.fill(HIST("hMassGamma"), v0.mGamma()); + + //Reconstruct directly here + ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass()); + ROOT::Math::PtEtaPhiMVector v2(v0.pt(), v0.eta(), v0.phi(), hadronMass); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float invmass = v12.M(); + // float massJpsi = 3.0969; + float mass = invmass - dilepton.mass() + 3.096; + float pt = v12.Pt(); + registry.fill(HIST("hMassChic"), mass); + registry.fill(HIST("hPtV0"), v0.pt()); + registry.fill(HIST("hPtV0Pos"), v0_pos.pt()); } } // loop over v0s } // loop over dileptons @@ -4032,6 +4081,22 @@ struct AnalysisDileptonV0 { } } + void processChicToJpsiGamma(soa::Filtered const& events, + MyBarrelTracks const& tracks, + MyV0s const& v0s, FullTracksExtEl const& fulltracks, + soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + + for (auto event : events) { + auto groupedDileptons = dileptons.sliceBy(dileptonsPerCollision, event.collisionId()); + auto groupedV0s = v0s.sliceBy(v0sPerCollision, event.collisionId()); + runDileptonV0(event, tracks, groupedV0s, fulltracks, groupedDileptons); + } + } + void processV0(soa::Filtered const& events, MyV0s const& v0s) { @@ -4053,7 +4118,7 @@ struct AnalysisDileptonV0 { if (v0.mK0Short() < fConfigV0Cuts.V0MassLow || v0.mK0Short() > fConfigV0Cuts.V0MassHigh) { continue; } - registry.fill(HIST("hV0Selected"), v0.mK0Short()); + registry.fill(HIST("hMassK0s"), v0.mK0Short()); } } } @@ -4092,6 +4157,7 @@ struct AnalysisDileptonV0 { } PROCESS_SWITCH(AnalysisDileptonV0, processB0ToJpsiK0s, "Run B0 to J/psi K0s analysis", false); + PROCESS_SWITCH(AnalysisDileptonV0, processChicToJpsiGamma, "Run Chic to J/psi Gamma analysis", false); PROCESS_SWITCH(AnalysisDileptonV0, processV0, "Run V0 analysis", false); PROCESS_SWITCH(AnalysisDileptonV0, processDilepton, "Run dilepton analysis", false); PROCESS_SWITCH(AnalysisDileptonV0, processDummy, "Dummy function", true);