Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
22 changes: 18 additions & 4 deletions DPG/Tasks/AOTTrack/D0CalibTables.h
Original file line number Diff line number Diff line change
Expand Up @@ -242,6 +242,20 @@
// column labels
static const std::vector<std::string> labelsCutVarCand = {"delta inv. mass", "max d0d0", "max pointing angle", "max pointing angle XY", "min cos pointing angle", "min cos pointing angle XY", "min norm decay length", "min norm decay length XY", "min decay length", "min decay length XY"};

// var indices for the cuts
enum CutsIds : uint8_t {
DeltaMass = 0,
MaxD0D0 = 1,
MaxPointingAngle = 2,
MaxPointingAngleXY = 3,
MinCosPointingAngle = 4,
MinCosPointingAngleXY = 5,
MinNormDecayLength = 6,
MinNormDecayLengthXY = 7,
MinDecayLength = 8,
MinDecayLengthXY = 9
};

static constexpr int NBinsPtMl = 10;
// default values for the pT bin edges (can be used to configure histogram axis)
// offset by 1 from the bin numbers in cuts array
Expand Down Expand Up @@ -324,10 +338,10 @@
DECLARE_SOA_COLUMN(TpcNumSigmaKa, tpcNumSigmaKa, int8_t); //! compressed NsigmaTPC for kaons
DECLARE_SOA_COLUMN(TofNumSigmaPi, tofNumSigmaPi, int8_t); //! compressed NsigmaTOF for pions
DECLARE_SOA_COLUMN(TofNumSigmaKa, tofNumSigmaKa, int8_t); //! compressed NsigmaTOF for kaons
DECLARE_SOA_COLUMN(ITSChi2NCl, itsChi2NCl, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column
DECLARE_SOA_COLUMN(TPCChi2NCl, tpcChi2NCl, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column
DECLARE_SOA_COLUMN(TRDChi2, trdChi2, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column
DECLARE_SOA_COLUMN(TOFChi2, tofChi2, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column
DECLARE_SOA_COLUMN(ITSChi2NCl, itsChi2NCl, uint8_t); // o2-linter: disable=name/o2-column

Check failure on line 341 in DPG/Tasks/AOTTrack/D0CalibTables.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/o2-column]

Use UpperCamelCase for names of O2 columns and matching lowerCamelCase names for their getters.
DECLARE_SOA_COLUMN(TPCChi2NCl, tpcChi2NCl, uint8_t); // o2-linter: disable=name/o2-column

Check failure on line 342 in DPG/Tasks/AOTTrack/D0CalibTables.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/o2-column]

Use UpperCamelCase for names of O2 columns and matching lowerCamelCase names for their getters.
DECLARE_SOA_COLUMN(TRDChi2, trdChi2, uint8_t); // o2-linter: disable=name/o2-column

Check failure on line 343 in DPG/Tasks/AOTTrack/D0CalibTables.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/o2-column]

Use UpperCamelCase for names of O2 columns and matching lowerCamelCase names for their getters.
DECLARE_SOA_COLUMN(TOFChi2, tofChi2, uint8_t); // o2-linter: disable=name/o2-column

Check failure on line 344 in DPG/Tasks/AOTTrack/D0CalibTables.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/o2-column]

Use UpperCamelCase for names of O2 columns and matching lowerCamelCase names for their getters.
DECLARE_SOA_COLUMN(CmoPrimUnfm80, cmoPrimUnfm80, uint8_t);
DECLARE_SOA_COLUMN(CmoFV0AUnfm80, cmoFV0AUnfm80, uint8_t);
DECLARE_SOA_COLUMN(CmoFT0AUnfm80, cmoFT0AUnfm80, uint8_t);
Expand Down
110 changes: 74 additions & 36 deletions DPG/Tasks/AOTTrack/derivedDataCreatorD0Calibration.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,8 @@
#include <array>
#include <cmath>
#include <cstdint>
#include <map>
#include <string>
#include <unordered_map>
#include <vector>

using namespace o2;
Expand Down Expand Up @@ -192,8 +192,18 @@ struct DerivedDataCreatorD0Calibration {
TTrackMeanOccs const&,
TTrackQa const&)
{
std::map<int, int> selectedCollisions; // map with indices of selected collisions (key: original AOD Collision table index, value: D0 collision index)
std::map<int, int> selectedTracks; // map with indices of selected tracks (key: original AOD Track table index, value: D0 daughter track index)
std::unordered_map<int, int> selectedCollisions; // map with indices of selected collisions (key: original AOD Collision table index, value: D0 collision index)
std::unordered_map<int, int> selectedTracks; // map with indices of selected tracks (key: original AOD Track table index, value: D0 daughter track index)

std::vector<float> bdtScoresD0;
std::vector<float> bdtScoresD0bar;
std::vector<float> featuresCandD0;
std::vector<float> featuresCandD0bar;
bdtScoresD0.resize(3);
bdtScoresD0bar.resize(3);
featuresCandD0.resize(16);
featuresCandD0bar.resize(16);
auto const& cuts = cfgCandCuts.topologicalCuts;

for (auto const& collision : collisions) {

Expand All @@ -217,6 +227,13 @@ struct DerivedDataCreatorD0Calibration {
eventIR.setFromLong(bc.globalBC());

auto groupedTrackIndices = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex());
std::vector<o2::dataformats::DCA> dcaNeg; // we cache the negative track, to avoid the same computation for each positive track
dcaNeg.resize(groupedTrackIndices.size());
std::vector<bool> isSelectedNeg(groupedTrackIndices.size(), true);
std::vector<bool> isPropagatedNeg(groupedTrackIndices.size(), false);
std::vector<int> pidTrackNegKaon(groupedTrackIndices.size(), -1);
std::vector<int> pidTrackNegPion(groupedTrackIndices.size(), -1);

for (auto const& trackIndexPos : groupedTrackIndices) {
auto trackPos = trackIndexPos.template track_as<TTracks>();
// track selections
Expand Down Expand Up @@ -255,27 +272,40 @@ struct DerivedDataCreatorD0Calibration {
pidTrackPosPion = selectorPion.statusTpcAndTof(trackPos);
}

int iTrack{-1};
for (auto const& trackIndexNeg : groupedTrackIndices) {
iTrack++;
if (!isSelectedNeg[iTrack]) {
continue;
}
auto trackNeg = trackIndexNeg.template track_as<TTracks>();
// track selections
if (trackNeg.sign() > 0) { // second negative track
isSelectedNeg[iTrack] = false;
continue;
}
if (!trackNeg.isGlobalTrackWoDCA()) {
isSelectedNeg[iTrack] = false;
continue;
}
if (trackNeg.pt() < cfgTrackCuts.ptMin) {
isSelectedNeg[iTrack] = false;
continue;
}
if (std::abs(trackNeg.eta()) > cfgTrackCuts.absEtaMax) {
isSelectedNeg[iTrack] = false;
continue;
}
auto trackParCovNeg = getTrackParCov(trackNeg);
o2::dataformats::DCA dcaNeg(trackNeg.dcaXY(), trackNeg.dcaZ(), trackNeg.cYY(), trackNeg.cZY(), trackNeg.cZZ());
if (trackNeg.collisionId() != collision.globalIndex()) {
trackParCovNeg.propagateToDCA(primaryVertex, bz, &dcaNeg);
if (!isPropagatedNeg[iTrack]) {
dcaNeg[iTrack] = o2::dataformats::DCA(trackNeg.dcaXY(), trackNeg.dcaZ(), trackNeg.cYY(), trackNeg.cZY(), trackNeg.cZZ());
if (trackNeg.collisionId() != collision.globalIndex()) {
trackParCovNeg.propagateToDCA(primaryVertex, bz, &dcaNeg[iTrack]);
}
isPropagatedNeg[iTrack] = true;
}
if (!isSelectedTrackDca(cfgTrackCuts.binsPt, cfgTrackCuts.limitsDca, trackParCovNeg.getPt(), dcaNeg.getY(), dcaNeg.getZ())) {
if (!isSelectedTrackDca(cfgTrackCuts.binsPt, cfgTrackCuts.limitsDca, trackParCovNeg.getPt(), dcaNeg[iTrack].getY(), dcaNeg[iTrack].getZ())) {
isSelectedNeg[iTrack] = false;
continue;
}

Expand Down Expand Up @@ -307,47 +337,45 @@ struct DerivedDataCreatorD0Calibration {
}
}

int pidTrackNegKaon{-1};
int pidTrackNegPion{-1};
if (cfgTrackCuts.usePidTpcOnly) {
/// kaon TPC PID negative daughter
pidTrackNegKaon = selectorKaon.statusTpc(trackNeg);
pidTrackNegKaon[iTrack] = selectorKaon.statusTpc(trackNeg);
/// pion TPC PID negative daughter
pidTrackNegPion = selectorPion.statusTpc(trackNeg);
pidTrackNegPion[iTrack] = selectorPion.statusTpc(trackNeg);
} else {
/// kaon TPC, TOF PID negative daughter
pidTrackNegKaon = selectorKaon.statusTpcAndTof(trackNeg);
pidTrackNegKaon[iTrack] = selectorKaon.statusTpcAndTof(trackNeg);
/// pion TPC, TOF PID negative daughter
pidTrackNegPion = selectorPion.statusTpcAndTof(trackNeg);
pidTrackNegPion[iTrack] = selectorPion.statusTpcAndTof(trackNeg);
}

// PID
uint8_t massHypo{D0MassHypo::D0AndD0Bar}; // both mass hypotheses a priori
if (pidTrackPosPion == TrackSelectorPID::Rejected || pidTrackNegKaon == TrackSelectorPID::Rejected) {
if (pidTrackPosPion == TrackSelectorPID::Rejected || pidTrackNegKaon[iTrack] == TrackSelectorPID::Rejected) {
massHypo -= D0MassHypo::D0; // exclude D0
}
if (pidTrackNegPion == TrackSelectorPID::Rejected || pidTrackPosKaon == TrackSelectorPID::Rejected) {
if (pidTrackNegPion[iTrack] == TrackSelectorPID::Rejected || pidTrackPosKaon == TrackSelectorPID::Rejected) {
massHypo -= D0MassHypo::D0Bar; // exclude D0Bar
}
if (massHypo == 0) {
continue;
}

// d0xd0
if (dcaPos.getY() * dcaNeg.getY() > cfgCandCuts.topologicalCuts->get(ptBinNoVtxD0, "max d0d0")) {
if (dcaPos.getY() * dcaNeg[iTrack].getY() > cuts->get(ptBinNoVtxD0, CutsIds::MaxD0D0)) {
continue;
}

// invariant mass
if (massHypo == D0MassHypo::D0 || massHypo == D0MassHypo::D0AndD0Bar) {
float invMassNoVtxD0 = RecoDecay::m(std::array{trackPos.pVector(), trackNeg.pVector()}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus});
if (std::abs(invMassNoVtxD0 - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinNoVtxD0, "delta inv. mass") + invMassTolerance) {
if (std::abs(invMassNoVtxD0 - o2::constants::physics::MassD0) > cuts->get(ptBinNoVtxD0, CutsIds::DeltaMass) + invMassTolerance) {
massHypo -= D0MassHypo::D0;
}
}
if (massHypo >= D0MassHypo::D0Bar) {
float invMassNoVtxD0bar = RecoDecay::m(std::array{trackNeg.pVector(), trackPos.pVector()}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus});
if (std::abs(invMassNoVtxD0bar - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinNoVtxD0, "delta inv. mass") + invMassTolerance) {
if (std::abs(invMassNoVtxD0bar - o2::constants::physics::MassD0) > cuts->get(ptBinNoVtxD0, CutsIds::DeltaMass) + invMassTolerance) {
massHypo -= D0MassHypo::D0Bar;
}
}
Expand Down Expand Up @@ -383,64 +411,72 @@ struct DerivedDataCreatorD0Calibration {
}

// d0xd0
if (dcaPos.getY() * dcaNeg.getY() > cfgCandCuts.topologicalCuts->get(ptBinD0, "max d0d0")) {
if (dcaPos.getY() * dcaNeg[iTrack].getY() > cuts->get(ptBinD0, CutsIds::MaxD0D0)) {
continue;
}
// cospa
float cosPaD0 = RecoDecay::cpa(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, pVecD0);
if (cosPaD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min cos pointing angle")) {
if (cosPaD0 < cuts->get(ptBinD0, CutsIds::MinCosPointingAngle)) {
continue;
}
// pointing angle
float paD0 = std::acos(cosPaD0);
if (paD0 > cfgCandCuts.topologicalCuts->get(ptBinD0, "max pointing angle")) {
if (paD0 > cuts->get(ptBinD0, CutsIds::MaxPointingAngle)) {
continue;
}
// cospa XY
float cosPaXYD0 = RecoDecay::cpaXY(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, pVecD0);
if (cosPaXYD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min cos pointing angle XY")) {
if (cosPaXYD0 < cuts->get(ptBinD0, CutsIds::MinCosPointingAngleXY)) {
continue;
}
// pointing angle XY
float paXYD0 = std::acos(cosPaXYD0);
if (paXYD0 > cfgCandCuts.topologicalCuts->get(ptBinD0, "max pointing angle XY")) {
if (paXYD0 > cuts->get(ptBinD0, CutsIds::MaxPointingAngleXY)) {
continue;
}
// decay length
float decLenD0 = RecoDecay::distance(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex);
if (decLenD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min decay length")) {
if (decLenD0 < cuts->get(ptBinD0, CutsIds::MinDecayLength)) {
continue;
}
// decay length XY
float decLenXYD0 = RecoDecay::distanceXY(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex);
if (decLenXYD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min decay length XY")) {
if (decLenXYD0 < cuts->get(ptBinD0, CutsIds::MinDecayLengthXY)) {
continue;
}
// normalised decay length
float phi{0.f}, theta{0.f};
getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta);
float errorDecayLengthD0 = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta));
if (decLenD0 / errorDecayLengthD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min norm decay length")) {
if (decLenD0 / errorDecayLengthD0 < cuts->get(ptBinD0, CutsIds::MinNormDecayLength)) {
continue;
}
// normalised decay length XY
float errorDecayLengthXYD0 = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.f) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.f));
if (decLenXYD0 / errorDecayLengthXYD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min norm decay length XY")) {
if (decLenXYD0 / errorDecayLengthXYD0 < cuts->get(ptBinD0, CutsIds::MinNormDecayLengthXY)) {
continue;
}

float invMassD0{0.f}, invMassD0bar{0.f};
float cosThetaStarD0{0.f}, cosThetaStarD0bar{0.f};
std::vector<float> bdtScoresD0{0.f, 1.f, 1.f}, bdtScoresD0bar{0.f, 1.f, 1.f}; // always selected a priori
// always selected a priori
bdtScoresD0[0] = 0.f;
bdtScoresD0[1] = 1.f;
bdtScoresD0[2] = 1.f;
bdtScoresD0bar[0] = 0.f;
bdtScoresD0bar[1] = 1.f;
bdtScoresD0bar[2] = 1.f;
if (massHypo == D0MassHypo::D0 || massHypo == D0MassHypo::D0AndD0Bar) {
invMassD0 = RecoDecay::m(std::array{pVecPos, pVecNeg}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus});
if (std::abs(invMassD0 - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinD0, "delta inv. mass")) {
if (std::abs(invMassD0 - o2::constants::physics::MassD0) > cuts->get(ptBinD0, CutsIds::DeltaMass)) {
massHypo -= D0MassHypo::D0;
bdtScoresD0 = std::vector<float>{1.f, 0.f, 0.f};
bdtScoresD0[0] = 1.f;
bdtScoresD0[1] = 0.f;
bdtScoresD0[2] = 0.f;
} else {
// apply BDT models
if (cfgMl.apply) {
std::vector<float> featuresCandD0 = {dcaPos.getY(), dcaNeg.getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg.getY(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi())};
featuresCandD0 = {dcaPos.getY(), dcaNeg[iTrack].getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg[iTrack].getY(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi())};
if (!mlResponse.isSelectedMl(featuresCandD0, ptD0, bdtScoresD0)) {
massHypo -= D0MassHypo::D0;
} else { // selected, we compute cost*
Expand All @@ -457,13 +493,15 @@ struct DerivedDataCreatorD0Calibration {
}
if (massHypo >= D0MassHypo::D0Bar) {
invMassD0bar = RecoDecay::m(std::array{pVecNeg, pVecPos}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus});
if (std::abs(invMassD0bar - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinD0, "delta inv. mass")) {
if (std::abs(invMassD0bar - o2::constants::physics::MassD0) > cuts->get(ptBinD0, CutsIds::DeltaMass)) {
massHypo -= D0MassHypo::D0Bar;
bdtScoresD0bar = std::vector<float>{1.f, 0.f, 0.f};
bdtScoresD0bar[0] = 1.f;
bdtScoresD0bar[1] = 0.f;
bdtScoresD0bar[2] = 0.f;
} else {
// apply BDT models
if (cfgMl.apply) {
std::vector<float> featuresCandD0bar = {dcaPos.getY(), dcaNeg.getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg.getY(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi())};
featuresCandD0bar = {dcaPos.getY(), dcaNeg[iTrack].getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg[iTrack].getY(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma<false>(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi())};
if (!mlResponse.isSelectedMl(featuresCandD0bar, ptD0, bdtScoresD0bar)) {
massHypo -= D0MassHypo::D0Bar;
} else { // selected, we compute cost*
Expand Down Expand Up @@ -823,8 +861,8 @@ struct DerivedDataCreatorD0Calibration {
deltaRefGloParamQ2Pt,
deltaTOFdX,
deltaTOFdZ,
o2::math_utils::detail::truncateFloatFraction(dcaNeg.getY(), precisionDcaMask),
o2::math_utils::detail::truncateFloatFraction(dcaNeg.getZ(), precisionDcaMask),
o2::math_utils::detail::truncateFloatFraction(dcaNeg[iTrack].getY(), precisionDcaMask),
o2::math_utils::detail::truncateFloatFraction(dcaNeg[iTrack].getZ(), precisionDcaMask),
getCompressedNumSigmaPid(trackNeg.tpcNSigmaPi()),
getCompressedNumSigmaPid(trackNeg.tpcNSigmaKa()),
getCompressedNumSigmaPid(trackNeg.tofNSigmaPi()),
Expand Down
Loading