Skip to content

PWGHF: Adding D0 to taskFlowCharmHadrons, including Ml. #5087

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wuctlby wants to merge 7 commits into from
Closed

PWGHF: Adding D0 to taskFlowCharmHadrons, including Ml. #5087

Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e7c6bb0
version 2: taskFlowCharmHadrons
wuctlby Mar 8, 2024
fd3c13c
No selection flag in Thn
wuctlby Mar 11, 2024
65681d5
No selection Flag
wuctlby Mar 11, 2024
3f356f8
Add D0 to taskFlowCharmHadrons, , including Ml
wuctlby Mar 11, 2024
1b66201
PWGHF: Adding D0 to taskFlowCharmHadrons, , including Ml
wuctlby Mar 11, 2024
0331733
Please consider the following formatting changes
alibuild Mar 11, 2024
b049baa
Merge pull request #2 from alibuild/alibot-cleanup-5087
wuctlby Mar 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 77 additions & 15 deletions PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,8 @@ using namespace o2::framework::expressions;

enum DecayChannel { DplusToPiKPi = 0,
DsToKKPi,
DsToPiKK };
DsToPiKK,
D0ToPiK };

enum centralityEstimator { V0A = 0,
T0M,
Expand All @@ -50,8 +51,9 @@ struct HfTaskFlowCharmHadrons {
Configurable<float> zVtxMax{"zVtxMax", 10., "Max vertex coordinate z"};
Configurable<int> harmonic{"harmonic", 2, "harmonic number"};
Configurable<int> qvecDetector{"qvecDetector", 0, "Detector for Q vector estimation (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3, TPC Pos: 4, TPC Neg: 5)"};
Configurable<int> centDetector{"centDetector", 0, "Detector for centrality estimation (V0A: 0, T0M: 1, T0A: 2, T0C: 3"};
Configurable<int> centDetector{"centDetector", 0, "Detector for centrality estimation (V0A: 0, T0M: 1, T0A: 2, T0C: 3)"};
Configurable<int> selectionFlag{"selectionFlag", 1, "Selection Flag for hadron (e.g. 1 for skimming, 3 for topo. and kine., 7 for PID)"};
Configurable<int> prongNum{"prongNum", 3, "Number of candidate's prong (For D0, set selectionFlag = 1 and prongNum = 2)"};
Copy link
Collaborator

@vkucera vkucera Mar 11, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Configurable<int> prongNum{"prongNum", 3, "Number of candidate's prong (For D0, set selectionFlag = 1 and prongNum = 2)"};
Configurable<int> nProngs{"nProngs", 3, "Number of candidate's prong (For D0, set selectionFlag = 1 and nProngs = 2)"};

Why does this have to be set manually?

Configurable<bool> storeMl{"storeMl", false, "Flag to store ML scores"};
Configurable<bool> saveEpResoHisto{"saveEpResoHisto", false, "Flag to save event plane resolution histogram"};
Configurable<std::vector<int>> classMl{"classMl", {0, 2}, "Indexes of BDT scores to be stored. Two indexes max."};
Expand All @@ -69,10 +71,13 @@ struct HfTaskFlowCharmHadrons {
using CandDsData = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDsToKKPi>>;
using CandDplusDatawMl = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDplusToPiKPi, aod::HfMlDplusToPiKPi>>;
using CandDplusData = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDplusToPiKPi>>;
using CandD0DatawMl = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfMlD0>>;
using CandD0Data = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0>>;
using CollsWithQvecs = soa::Join<aod::Collisions, aod::EvSels, aod::QvectorFT0Cs, aod::QvectorFT0As, aod::QvectorFT0Ms, aod::QvectorFV0As, aod::QvectorBPoss, aod::QvectorBNegs, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs>;

Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag;
Filter filterSelectDplusCandidates = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag;
Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag;

Partition<CandDsData> selectedDsToKKPi = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag;
Partition<CandDsData> selectedDsToPiKK = aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag;
Expand Down Expand Up @@ -135,7 +140,8 @@ struct HfTaskFlowCharmHadrons {
/// \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
template <typename T1>
/// \param DecayChannel
template <int DecayChannel, typename T1>
void getQvecDtracks(const T1& cand,
std::vector<float>& tracksQx,
std::vector<float>& tracksQy,
Expand All @@ -150,17 +156,20 @@ struct HfTaskFlowCharmHadrons {
float pYtrack1 = cand.pyProng1();
float pTtrack1 = cand.ptProng1();
float phiTrack1 = std::atan2(pYtrack1, pXtrack1);
float pXtrack2 = cand.pxProng2();
float pYtrack2 = cand.pyProng2();
float pTtrack2 = cand.ptProng2();
float phiTrack2 = std::atan2(pYtrack2, pXtrack2);

tracksQx.push_back(cos(harmonic * phiTrack0) * pTtrack0 / ampl);
tracksQy.push_back(sin(harmonic * phiTrack0) * pTtrack0 / ampl);
tracksQx.push_back(cos(harmonic * phiTrack1) * pTtrack1 / ampl);
tracksQy.push_back(sin(harmonic * phiTrack1) * pTtrack1 / ampl);
tracksQx.push_back(cos(harmonic * phiTrack2) * pTtrack2 / ampl);
tracksQy.push_back(sin(harmonic * phiTrack2) * pTtrack2 / ampl);

if constexpr (DecayChannel != DecayChannel::D0ToPiK) {
Copy link
Collaborator

@vkucera vkucera Mar 11, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please avoid this kind of confusing naming patterns. DecayChannel on the left is a variable, DecayChannel on the right is a type.

Copy link
Contributor Author

@wuctlby wuctlby Mar 11, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I have closed this PR due to a failed submission and will resubmit. I'll change or explain existing comments.

float pXtrack2 = cand.pxProng2();
float pYtrack2 = cand.pyProng2();
float pTtrack2 = cand.ptProng2();
float phiTrack2 = std::atan2(pYtrack2, pXtrack2);
tracksQx.push_back(cos(harmonic * phiTrack2) * pTtrack2 / ampl);
tracksQy.push_back(sin(harmonic * phiTrack2) * pTtrack2 / ampl);
Comment on lines +170 to +171
Copy link
Collaborator

@vkucera vkucera Mar 11, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explicitly state the namespace you use.

}
}

/// Compute the delta psi in the range [0, pi/harmonic]
Expand Down Expand Up @@ -188,6 +197,7 @@ struct HfTaskFlowCharmHadrons {
/// \param sp is the scalar product
/// \param evtPlReso is the event plane resolution
/// \param outputMl are the ML scores
/// \param selectionFlag for D0, only 0 or1
void fillThn(float& mass,
float& pt,
float& cent,
Expand Down Expand Up @@ -320,13 +330,22 @@ struct HfTaskFlowCharmHadrons {

// 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 - 3.;
float ampl = amplQVec - float(prongNum);
std::vector<float> tracksQx = {};
std::vector<float> tracksQy = {};
getQvecDtracks(candidate, tracksQx, tracksQy, ampl);
for (unsigned int itrack = 0; itrack < 3; itrack++) {
xQVec -= tracksQx[itrack];
yQVec -= tracksQy[itrack];

if constexpr (DecayChannel == DecayChannel::D0ToPiK) {
getQvecDtracks<DecayChannel>(candidate, tracksQx, tracksQy, ampl);
for (unsigned int itrack = 0; itrack < 2; itrack++) {
xQVec -= tracksQx[itrack];
yQVec -= tracksQy[itrack];
}
} else {
getQvecDtracks<DecayChannel>(candidate, tracksQx, tracksQy, ampl);
for (unsigned int itrack = 0; itrack < 3; itrack++) {
xQVec -= tracksQx[itrack];
yQVec -= tracksQy[itrack];
}
}
}

Expand All @@ -335,7 +354,34 @@ struct HfTaskFlowCharmHadrons {
float scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec;
float cosDeltaPhi = std::cos(harmonic * (phiCand - evtPl));

fillThn(massCand, ptCand, cent, cosNPhi, cosDeltaPhi, scalprodCand, outputMl);
if constexpr (std::is_same<T1, CandD0Data>::value || std::is_same<T1, CandD0DatawMl>::value) {
std::vector<float> outputMlD0 = {-999., -999.};
std::vector<float> outputMlD0bar = {-999., -999.};

if (candidate.isSelD0() >= selectionFlag) {
massCand = hfHelper.invMassD0ToPiK(candidate);
if constexpr (std::is_same<T1, CandD0DatawMl>::value) {
std::cout << "Here is CandD0DatawMl" << std::endl;
for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) {
outputMlD0[iclass] = candidate.mlProbD0()[classMl->at(iclass)];
std::cout << "Here is iclass " << iclass << "\t Here is outputMl " << outputMlD0[iclass] << std::endl;
}
}
fillThn(massCand, ptCand, cent, cosNPhi, cosDeltaPhi, scalprodCand, outputMlD0);
}
if (candidate.isSelD0bar() >= selectionFlag) {
massCand = hfHelper.invMassD0barToKPi(candidate);
if constexpr (std::is_same<T1, CandD0DatawMl>::value) {
for (unsigned int iclass = 0; iclass < classMl->size(); iclass++)
outputMlD0bar[iclass] = candidate.mlProbD0bar()[classMl->at(iclass)];
}
fillThn(massCand, ptCand, cent, cosNPhi, cosDeltaPhi, scalprodCand, outputMlD0bar);
}
} // TODO: Whether to put all the mass calculations here

if constexpr (DecayChannel != DecayChannel::D0ToPiK) {
fillThn(massCand, ptCand, cent, cosNPhi, cosDeltaPhi, scalprodCand, outputMl);
}
}
}

Expand Down Expand Up @@ -373,6 +419,22 @@ struct HfTaskFlowCharmHadrons {
}
PROCESS_SWITCH(HfTaskFlowCharmHadrons, processDplus, "Process Dplus candidates", true);

// D0 with ML
void processD0Ml(CollsWithQvecs::iterator const& collision,
CandD0DatawMl const& candidatesD0)
{
runFlowAnalysis<DecayChannel::D0ToPiK, CandD0DatawMl>(collision, candidatesD0);
}
PROCESS_SWITCH(HfTaskFlowCharmHadrons, processD0Ml, "Process D0 candidates with ML", false);

// D0 with rectangular cuts
void processD0(CollsWithQvecs::iterator const& collision,
CandD0Data const& candidatesD0)
{
runFlowAnalysis<DecayChannel::D0ToPiK, CandD0Data>(collision, candidatesD0);
}
PROCESS_SWITCH(HfTaskFlowCharmHadrons, processD0, "Process D0 candidates", false);

// Resolution
void processResolution(CollsWithQvecs::iterator const& collision)
{
Expand Down