diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 862b4910c72..03e475bc7ed 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -14,6 +14,7 @@ /// /// \author S. Politanò, INFN Torino, Italy /// \author Wu Chuntai, CUG, China +/// \author Ran Tu, Fudan University, China #include #include @@ -48,7 +49,8 @@ enum DecayChannel { DplusToPiKPi = 0, LcToPKPi, LcToPiKP, XicToPKPi, - XicToPiKP + XicToPiKP, + Xic0ToXiPi }; enum QvecEstimator { FV0A = 0, @@ -75,25 +77,11 @@ struct HfTaskFlowCharmHadrons { Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable> classMl{"classMl", {0, 2}, "Indexes of BDT scores to be stored. Two indexes max."}; - ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {100, 1.78, 2.05}, ""}; - ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {10, 0., 10.}, ""}; - ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {10000, 0., 100.}, ""}; - ConfigurableAxis thnConfigAxisCosNPhi{"thnConfigAxisCosNPhi", {100, -1., 1.}, ""}; - ConfigurableAxis thnConfigAxisPsi{"thnConfigAxisPsi", {6000, 2. * TMath::Pi(), 2. * TMath::Pi()}, ""}; - ConfigurableAxis thnConfigAxisCosDeltaPhi{"thnConfigAxisCosDeltaPhi", {100, -1., 1.}, ""}; - ConfigurableAxis thnConfigAxisScalarProd{"thnConfigAxisScalarProd", {100, 0., 1.}, ""}; - ConfigurableAxis thnConfigAxisMlOne{"thnConfigAxisMlOne", {1000, 0., 1.}, ""}; - ConfigurableAxis thnConfigAxisMlTwo{"thnConfigAxisMlTwo", {1000, 0., 1.}, ""}; - ConfigurableAxis thnConfigAxisOccupancyITS{"thnConfigAxisOccupancyITS", {14, 0, 14000}, ""}; - ConfigurableAxis thnConfigAxisOccupancyFT0C{"thnConfigAxisOccupancyFT0C", {14, 0, 140000}, ""}; - ConfigurableAxis thnConfigAxisNoSameBunchPileup{"thnConfigAxisNoSameBunchPileup", {2, 0, 2}, ""}; - ConfigurableAxis thnConfigAxisOccupancy{"thnConfigAxisOccupancy", {2, 0, 2}, ""}; - ConfigurableAxis thnConfigAxisNoCollInTimeRangeNarrow{"thnConfigAxisNoCollInTimeRangeNarrow", {2, 0, 2}, ""}; - ConfigurableAxis thnConfigAxisNoCollInTimeRangeStandard{"thnConfigAxisNoCollInTimeRangeStandard", {2, 0, 2}, ""}; - ConfigurableAxis thnConfigAxisNoCollInRofStandard{"thnConfigAxisNoCollInRofStandard", {2, 0, 2}, ""}; - ConfigurableAxis thnConfigAxisResoFT0cFV0a{"thnConfigAxisResoFT0cFV0a", {160, -8, 8}, ""}; - ConfigurableAxis thnConfigAxisResoFT0cTPCtot{"thnConfigAxisResoFT0cTPCtot", {160, -8, 8}, ""}; - ConfigurableAxis thnConfigAxisResoFV0aTPCtot{"thnConfigAxisResoFV0aTPCtot", {160, -8, 8}, ""}; + HfHelper hfHelper; + EventPlaneHelper epHelper; + HfEventSelection hfEvSel; // event selection and monitoring + o2::framework::Service ccdb; + SliceCache cache; using CandDsDataWMl = soa::Filtered>; using CandDsData = soa::Filtered>; @@ -103,6 +91,8 @@ struct HfTaskFlowCharmHadrons { using CandLcDataWMl = soa::Filtered>; using CandXicData = soa::Filtered>; using CandXicDataWMl = soa::Filtered>; + using CandXic0Data = soa::Filtered>; + using CandXic0DataWMl = soa::Filtered>; using CandD0DataWMl = soa::Filtered>; using CandD0Data = soa::Filtered>; using CollsWithQvecs = soa::Join; @@ -112,6 +102,7 @@ struct HfTaskFlowCharmHadrons { Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; Filter filterSelectLcCandidates = aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlag || aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlag; Filter filterSelectXicCandidates = aod::hf_sel_candidate_xic::isSelXicToPKPi >= selectionFlag || aod::hf_sel_candidate_xic::isSelXicToPiKP >= selectionFlag; + Filter filterSelectXic0Candidates = aod::hf_sel_toxipi::resultSelections >= selectionFlag; Partition selectedDsToKKPi = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag; Partition selectedDsToPiKK = aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; @@ -129,12 +120,28 @@ struct HfTaskFlowCharmHadrons { Partition selectedXicToPiKP = aod::hf_sel_candidate_xic::isSelXicToPiKP >= selectionFlag; Partition selectedXicToPKPiWMl = aod::hf_sel_candidate_xic::isSelXicToPKPi >= selectionFlag; Partition selectedXicToPiKPWMl = aod::hf_sel_candidate_xic::isSelXicToPiKP >= selectionFlag; + Partition selectedXic0 = aod::hf_sel_toxipi::resultSelections >= selectionFlag; + Partition selectedXic0WMl = aod::hf_sel_toxipi::resultSelections >= selectionFlag; - SliceCache cache; - HfHelper hfHelper; - EventPlaneHelper epHelper; - HfEventSelection hfEvSel; // event selection and monitoring - o2::framework::Service ccdb; + ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {100, 1.78, 2.05}, ""}; + ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {10, 0., 10.}, ""}; + ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {10000, 0., 100.}, ""}; + ConfigurableAxis thnConfigAxisCosNPhi{"thnConfigAxisCosNPhi", {100, -1., 1.}, ""}; + ConfigurableAxis thnConfigAxisPsi{"thnConfigAxisPsi", {6000, 0, constants::math::TwoPI}, ""}; + ConfigurableAxis thnConfigAxisCosDeltaPhi{"thnConfigAxisCosDeltaPhi", {100, -1., 1.}, ""}; + ConfigurableAxis thnConfigAxisScalarProd{"thnConfigAxisScalarProd", {100, 0., 1.}, ""}; + ConfigurableAxis thnConfigAxisMlOne{"thnConfigAxisMlOne", {1000, 0., 1.}, ""}; + ConfigurableAxis thnConfigAxisMlTwo{"thnConfigAxisMlTwo", {1000, 0., 1.}, ""}; + ConfigurableAxis thnConfigAxisOccupancyITS{"thnConfigAxisOccupancyITS", {14, 0, 14000}, ""}; + ConfigurableAxis thnConfigAxisOccupancyFT0C{"thnConfigAxisOccupancyFT0C", {14, 0, 140000}, ""}; + ConfigurableAxis thnConfigAxisNoSameBunchPileup{"thnConfigAxisNoSameBunchPileup", {2, 0, 2}, ""}; + ConfigurableAxis thnConfigAxisOccupancy{"thnConfigAxisOccupancy", {2, 0, 2}, ""}; + ConfigurableAxis thnConfigAxisNoCollInTimeRangeNarrow{"thnConfigAxisNoCollInTimeRangeNarrow", {2, 0, 2}, ""}; + ConfigurableAxis thnConfigAxisNoCollInTimeRangeStandard{"thnConfigAxisNoCollInTimeRangeStandard", {2, 0, 2}, ""}; + ConfigurableAxis thnConfigAxisNoCollInRofStandard{"thnConfigAxisNoCollInRofStandard", {2, 0, 2}, ""}; + ConfigurableAxis thnConfigAxisResoFT0cFV0a{"thnConfigAxisResoFT0cFV0a", {160, -8, 8}, ""}; + ConfigurableAxis thnConfigAxisResoFT0cTPCtot{"thnConfigAxisResoFT0cTPCtot", {160, -8, 8}, ""}; + ConfigurableAxis thnConfigAxisResoFV0aTPCtot{"thnConfigAxisResoFV0aTPCtot", {160, -8, 8}, ""}; HistogramRegistry registry{"registry", {}}; @@ -231,15 +238,15 @@ struct HfTaskFlowCharmHadrons { } if (storeResoOccu) { - std::vector axes_reso = {thnAxisCent, thnAxisResoFT0cFV0a, thnAxisResoFT0cTPCtot, thnAxisResoFV0aTPCtot}; + std::vector axesReso = {thnAxisCent, thnAxisResoFT0cFV0a, thnAxisResoFT0cTPCtot, thnAxisResoFV0aTPCtot}; if (occEstimator == 1) { - axes_reso.insert(axes_reso.end(), {thnAxisOccupancyITS, thnAxisNoSameBunchPileup, thnAxisOccupancy, - thnAxisNoCollInTimeRangeNarrow, thnAxisNoCollInTimeRangeStandard, thnAxisNoCollInRofStandard}); + axesReso.insert(axesReso.end(), {thnAxisOccupancyITS, thnAxisNoSameBunchPileup, thnAxisOccupancy, + thnAxisNoCollInTimeRangeNarrow, thnAxisNoCollInTimeRangeStandard, thnAxisNoCollInRofStandard}); } else { - axes_reso.insert(axes_reso.end(), {thnAxisOccupancyFT0C, thnAxisNoSameBunchPileup, thnAxisOccupancy, - thnAxisNoCollInTimeRangeNarrow, thnAxisNoCollInTimeRangeStandard, thnAxisNoCollInRofStandard}); + axesReso.insert(axesReso.end(), {thnAxisOccupancyFT0C, thnAxisNoSameBunchPileup, thnAxisOccupancy, + thnAxisNoCollInTimeRangeNarrow, thnAxisNoCollInTimeRangeStandard, thnAxisNoCollInRofStandard}); } - registry.add("spReso/hSparseReso", "THn for resolution with occupancy", HistType::kTHnSparseF, axes_reso); + registry.add("spReso/hSparseReso", "THn for resolution with occupancy", HistType::kTHnSparseF, axesReso); } hfEvSel.addHistograms(registry); // collision monitoring @@ -285,6 +292,44 @@ struct HfTaskFlowCharmHadrons { } } + /// Compute the Q vector for the candidate's tracks + /// \param cand is the candidate + /// \param tracksQx is the X component of the Q vector for the tracks + /// \param tracksQy is the Y component of the Q vector for the tracks + /// \param channel is the decay channel + template + void getQvecXic0Tracks(const T1& cand, + std::vector& tracksQx, + std::vector& tracksQy, + float ampl) + { + // add possibility to consider different weights for the tracks, at the moment only pT is considered; + float pXTrack0 = cand.pxPosV0Dau(); + float pYTrack0 = cand.pyPosV0Dau(); + float pTTrack0 = std::hypot(pXTrack0, pYTrack0); + float phiTrack0 = std::atan2(pXTrack0, pYTrack0); + float pXTrack1 = cand.pxNegV0Dau(); + float pYTrack1 = cand.pyNegV0Dau(); + float pTTrack1 = std::hypot(pXTrack1, pYTrack1); + float phiTrack1 = std::atan2(pXTrack1, pYTrack1); + float pYTrack2 = cand.pxBachFromCasc(); + float pXTrack2 = cand.pyBachFromCasc(); + float pTTrack2 = std::hypot(pXTrack2, pYTrack2); + float phiTrack2 = std::atan2(pXTrack2, pYTrack2); + float pXTrack3 = cand.pxBachFromCharmBaryon(); + float pYTrack3 = cand.pyBachFromCharmBaryon(); + float pTTrack3 = std::hypot(pXTrack3, pYTrack3); + float phiTrack3 = std::atan2(pXTrack3, pYTrack3); + + tracksQx.push_back(std::cos(harmonic * phiTrack0) * pTTrack0 / ampl); + tracksQy.push_back(std::sin(harmonic * phiTrack0) * pTTrack0 / ampl); + tracksQx.push_back(std::cos(harmonic * phiTrack1) * pTTrack1 / ampl); + tracksQy.push_back(std::sin(harmonic * phiTrack1) * pTTrack1 / ampl); + tracksQx.push_back(std::cos(harmonic * phiTrack2) * pTTrack2 / ampl); + tracksQy.push_back(std::sin(harmonic * phiTrack2) * pTTrack2 / ampl); + tracksQx.push_back(std::cos(harmonic * phiTrack3) * pTTrack3 / ampl); + tracksQy.push_back(std::sin(harmonic * phiTrack3) * pTTrack3 / ampl); + } /// Compute the delta psi in the range [0, pi/harmonic] /// \param psi1 is the first angle /// \param psi2 is the second angle @@ -292,12 +337,7 @@ struct HfTaskFlowCharmHadrons { float getDeltaPsiInRange(float psi1, float psi2) { float deltaPsi = psi1 - psi2; - if (std::abs(deltaPsi) > constants::math::PI / harmonic) { - if (deltaPsi > 0.) - deltaPsi -= constants::math::TwoPI / harmonic; - else - deltaPsi += constants::math::TwoPI / harmonic; - } + deltaPsi = RecoDecay::constrainAngle(deltaPsi, -o2::constants::math::PI / harmonic, harmonic); return deltaPsi; } @@ -552,18 +592,36 @@ struct HfTaskFlowCharmHadrons { default: break; } + } else if constexpr (std::is_same_v || std::is_same_v) { + massCand = candidate.invMassCharmBaryon(); + if constexpr (std::is_same_v) { + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) + outputMl[iclass] = candidate.mlProbToXiPi()[classMl->at(iclass)]; + } } - float ptCand = candidate.pt(); - float phiCand = candidate.phi(); + float ptCand = 0.; + float phiCand = 0.; + + if constexpr (std::is_same_v || std::is_same_v) { + ptCand = candidate.kfptXic(); + phiCand = std::atan2(candidate.pxCharmBaryon(), candidate.pyCharmBaryon()); + } else { + ptCand = candidate.pt(); + phiCand = candidate.phi(); + } // If TPC is used for the SP estimation, the tracks of the hadron candidate must be removed from the TPC Q vector to avoid double counting if (qvecDetector == QvecEstimator::TPCNeg || qvecDetector == QvecEstimator::TPCPos) { float ampl = amplQVec - static_cast(nProngs); std::vector tracksQx = {}; std::vector tracksQy = {}; - - getQvecDtracks(candidate, tracksQx, tracksQy, ampl); + if constexpr (std::is_same_v || std::is_same_v) { + // std::cout<(candidate, tracksQx, tracksQy, ampl); + } for (auto iTrack{0u}; iTrack < tracksQx.size(); ++iTrack) { xQVec -= tracksQx[iTrack]; yQVec -= tracksQy[iTrack]; @@ -683,6 +741,24 @@ struct HfTaskFlowCharmHadrons { } PROCESS_SWITCH(HfTaskFlowCharmHadrons, processXic, "Process Xic candidates", false); + // Xic0 with ML + void processXic0Ml(CollsWithQvecs::iterator const& collision, + CandXic0DataWMl const&) + { + auto candsXic0WMl = selectedXic0WMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + runFlowAnalysis(collision, candsXic0WMl); + } + PROCESS_SWITCH(HfTaskFlowCharmHadrons, processXic0Ml, "Process Xic0 candidates with ML", false); + + // Xic0 + void processXic0(CollsWithQvecs::iterator const& collision, + CandXic0Data const&) + { + auto candsXic0 = selectedXic0->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + runFlowAnalysis(collision, candsXic0); + } + PROCESS_SWITCH(HfTaskFlowCharmHadrons, processXic0, "Process Xic0 candidates", false); + // Resolution void processResolution(CollsWithQvecs::iterator const& collision, aod::BCsWithTimestamps const& bcs)