From faae56d362c4faa9125d83571bcebf3b4fbb1506 Mon Sep 17 00:00:00 2001 From: creetz16 Date: Fri, 25 Oct 2024 10:59:18 +0200 Subject: [PATCH 1/2] Fix cov matrix entry and gen decay vertex array --- .../Nuspex/decay3bodybuilder.cxx | 2 +- .../TableProducer/Nuspex/threebodyKFTask.cxx | 24 ++++++------------- 2 files changed, 8 insertions(+), 18 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx index bfe1ad23f92..2da7dec2ea5 100644 --- a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx +++ b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx @@ -768,7 +768,7 @@ struct decay3bodyBuilder { auto trackParCovPVNeg = trackParCovNeg; auto trackParCovPVBach = trackParCovBach; mPV.setPos({collision.posX(), collision.posY(), collision.posZ()}); - mPV.setCov(collision.covXX(), collision.covXX(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); + mPV.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); o2::base::Propagator::Instance()->propagateToDCABxByBz(mPV, trackParCovPVPos, 2.f, matCorr, &mDcaInfoCovPos); o2::base::Propagator::Instance()->propagateToDCABxByBz(mPV, trackParCovPVNeg, 2.f, matCorr, &mDcaInfoCovNeg); o2::base::Propagator::Instance()->propagateToDCABxByBz(mPV, trackParCovPVBach, 2.f, matCorr, &mDcaInfoCovBach); diff --git a/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx b/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx index 25f2ee57181..0d560ecde46 100644 --- a/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx +++ b/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx @@ -198,9 +198,7 @@ struct threebodyKFTask { if ((MCvtx3body.pdgCode() == motherPdgCode && lMCTrack0.pdgCode() == 2212 && lMCTrack1.pdgCode() == -211 && lMCTrack2.pdgCode() == bachelorPdgCode) || (MCvtx3body.pdgCode() == -motherPdgCode && lMCTrack0.pdgCode() == 211 && lMCTrack1.pdgCode() == -2212 && lMCTrack2.pdgCode() == -bachelorPdgCode)) { vtx3bodyPDGcode = MCvtx3body.pdgCode(); - genDecVtx[0] = lMCTrack0.vx(); - genDecVtx[0] = lMCTrack0.vy(); - genDecVtx[0] = lMCTrack0.vz(); + genDecVtx = {lMCTrack0.vx(), lMCTrack0.vy(), lMCTrack0.vz()}; MClifetime = RecoDecay::sqrtSumOfSquares(lMCTrack2.vx() - MCvtx3body.vx(), lMCTrack2.vy() - MCvtx3body.vy(), lMCTrack2.vz() - MCvtx3body.vz()) * o2::constants::physics::MassHyperTriton / MCvtx3body.p(); genPhi = MCvtx3body.phi(); genEta = MCvtx3body.eta(); @@ -318,16 +316,12 @@ struct threebodyKFTask { std::array piMinusMom{0.f}; for (auto& mcparticleDaughter : mcparticle.template daughters_as()) { if (mcparticleDaughter.pdgCode() == 2212) { - protonMom[0] = mcparticleDaughter.px(); - protonMom[1] = mcparticleDaughter.py(); - protonMom[2] = mcparticleDaughter.pz(); + protonMom = {mcparticleDaughter.px(), mcparticleDaughter.py(), mcparticleDaughter.pz()}; } else if (mcparticleDaughter.pdgCode() == -211) { - piMinusMom[0] = mcparticleDaughter.px(); - piMinusMom[1] = mcparticleDaughter.py(); - piMinusMom[2] = mcparticleDaughter.pz(); + piMinusMom = {mcparticleDaughter.px(), mcparticleDaughter.py(), mcparticleDaughter.pz()}; } } - genMCmassPrPi = RecoDecay::m(array{array{protonMom[0], protonMom[1], protonMom[2]}, array{piMinusMom[0], piMinusMom[1], piMinusMom[2]}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); + genMCmassPrPi = RecoDecay::m(array{protonMom, piMinusMom}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); registry.fill(HIST("hTrueHypertritonMCMassPrPi_nonReco"), genMCmassPrPi); registry.fill(HIST("hTrueHypertritonMCPtProton_nonReco"), RecoDecay::sqrtSumOfSquares(protonMom[0], protonMom[1])); registry.fill(HIST("hTrueHypertritonMCPtPion_nonReco"), RecoDecay::sqrtSumOfSquares(piMinusMom[0], piMinusMom[1])); @@ -338,16 +332,12 @@ struct threebodyKFTask { std::array piPlusMom{0.f}; for (auto& mcparticleDaughter : mcparticle.template daughters_as()) { if (mcparticleDaughter.pdgCode() == -2212) { - antiProtonMom[0] = mcparticleDaughter.px(); - antiProtonMom[1] = mcparticleDaughter.py(); - antiProtonMom[2] = mcparticleDaughter.pz(); + antiProtonMom = {mcparticleDaughter.px(), mcparticleDaughter.py(), mcparticleDaughter.pz()}; } else if (mcparticleDaughter.pdgCode() == 211) { - piPlusMom[0] = mcparticleDaughter.px(); - piPlusMom[1] = mcparticleDaughter.py(); - piPlusMom[2] = mcparticleDaughter.pz(); + piPlusMom = {mcparticleDaughter.px(), mcparticleDaughter.py(), mcparticleDaughter.pz()}; } } - genMCmassPrPi = RecoDecay::m(array{array{antiProtonMom[0], antiProtonMom[1], antiProtonMom[2]}, array{piPlusMom[0], piPlusMom[1], piPlusMom[2]}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); + genMCmassPrPi = RecoDecay::m(array{antiProtonMom, piPlusMom}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); registry.fill(HIST("hTrueAntiHypertritonMCMassPrPi_nonReco"), genMCmassPrPi); registry.fill(HIST("hTrueAntiHypertritonMCPtProton_nonReco"), RecoDecay::sqrtSumOfSquares(antiProtonMom[0], antiProtonMom[1])); registry.fill(HIST("hTrueAntiHypertritonMCPtPion_nonReco"), RecoDecay::sqrtSumOfSquares(piPlusMom[0], piPlusMom[1])); From fe3c1224e7fbc0e8bd8c40a548633e85c0e7650b Mon Sep 17 00:00:00 2001 From: creetz16 Date: Tue, 5 Nov 2024 16:36:04 +0100 Subject: [PATCH 2/2] Add TPC signal to candidate table, add deuteron TOF PID, fix event properties filling --- PWGLF/DataModel/Vtx3BodyTables.h | 10 +++ .../Nuspex/decay3bodybuilder.cxx | 87 ++++++++++++------- .../TableProducer/Nuspex/threebodyKFTask.cxx | 4 + 3 files changed, 70 insertions(+), 31 deletions(-) diff --git a/PWGLF/DataModel/Vtx3BodyTables.h b/PWGLF/DataModel/Vtx3BodyTables.h index 26c8dff1424..b9a734e8f8a 100644 --- a/PWGLF/DataModel/Vtx3BodyTables.h +++ b/PWGLF/DataModel/Vtx3BodyTables.h @@ -386,6 +386,10 @@ DECLARE_SOA_COLUMN(Track2Sign, track2sign, float); //! sig DECLARE_SOA_COLUMN(TPCNSigmaProton, tpcnsigmaproton, float); //! nsigma of TPC PID of the proton daughter DECLARE_SOA_COLUMN(TPCNSigmaPion, tpcnsigmapion, float); //! nsigma of TPC PID of the pion daughter DECLARE_SOA_COLUMN(TPCNSigmaDeuteron, tpcnsigmadeuteron, float); //! nsigma of TPC PID of the bachelor daughter +DECLARE_SOA_COLUMN(TPCdEdxProton, tpcdedxproton, float); //! TPC dEdx of the proton daughter +DECLARE_SOA_COLUMN(TPCdEdxPion, tpcdedxpion, float); //! TPC dEdx of the pion daughter +DECLARE_SOA_COLUMN(TPCdEdxDeuteron, tpcdedxdeuteron, float); //! TPC dEdx of the bachelor daughter +DECLARE_SOA_COLUMN(TOFNSigmaDeuteron, tofnsigmadeuteron, float); //! nsigma of TOF PID of the bachelor daughter // Monte Carlo DECLARE_SOA_COLUMN(GenP, genp, float); //! generated momentum @@ -441,6 +445,8 @@ DECLARE_SOA_TABLE(KFVtx3BodyDatas, "AOD", "KFVTX3BODYDATA", kfvtx3body::DCATrackPosToPV, kfvtx3body::DCATrackNegToPV, kfvtx3body::DCATrackBachToPV, kfvtx3body::Track0Sign, kfvtx3body::Track1Sign, kfvtx3body::Track2Sign, // track sing: proton, pion, deuteron kfvtx3body::TPCNSigmaProton, kfvtx3body::TPCNSigmaPion, kfvtx3body::TPCNSigmaDeuteron, + kfvtx3body::TPCdEdxProton, kfvtx3body::TPCdEdxPion, kfvtx3body::TPCdEdxDeuteron, + kfvtx3body::TOFNSigmaDeuteron, // dynamic columns vtx3body::VtxRadius, @@ -505,6 +511,8 @@ DECLARE_SOA_TABLE(KFVtx3BodyDatasLite, "AOD", "KF3BODYLITE", kfvtx3body::DCATrackPosToPV, kfvtx3body::DCATrackNegToPV, kfvtx3body::DCATrackBachToPV, kfvtx3body::Track0Sign, kfvtx3body::Track1Sign, kfvtx3body::Track2Sign, // track sing: proton, pion, deuteron kfvtx3body::TPCNSigmaProton, kfvtx3body::TPCNSigmaPion, kfvtx3body::TPCNSigmaDeuteron, + kfvtx3body::TPCdEdxProton, kfvtx3body::TPCdEdxPion, kfvtx3body::TPCdEdxDeuteron, + kfvtx3body::TOFNSigmaDeuteron, // dynamic columns vtx3body::VtxRadius, @@ -560,6 +568,8 @@ DECLARE_SOA_TABLE(McKFVtx3BodyDatas, "AOD", "MCKF3BODYDATAS", kfvtx3body::DCATrackPosToPV, kfvtx3body::DCATrackNegToPV, kfvtx3body::DCATrackBachToPV, kfvtx3body::Track0Sign, kfvtx3body::Track1Sign, kfvtx3body::Track2Sign, // track sing: proton, pion, deuteron kfvtx3body::TPCNSigmaProton, kfvtx3body::TPCNSigmaPion, kfvtx3body::TPCNSigmaDeuteron, + kfvtx3body::TPCdEdxProton, kfvtx3body::TPCdEdxPion, kfvtx3body::TPCdEdxDeuteron, + kfvtx3body::TOFNSigmaDeuteron, // MC information kfvtx3body::GenP, diff --git a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx index 2da7dec2ea5..eabb56d36c3 100644 --- a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx +++ b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx @@ -629,10 +629,38 @@ struct decay3bodyBuilder { //------------------------------------------------------------------ // 3body candidate builder with KFParticle - template + template void buildVtx3BodyDataTableKFParticle(TCollision const& collision, aod::Decay3Bodys const& decay3bodys, int bachelorcharge = 1) { LOG(debug) << "buildVtx3BodyDataTableKFParticle called."; + + // initialise KF primary vertex + KFPVertex kfpVertex = createKFPVertexFromCollision(collision); + KFParticle kfpv(kfpVertex); + LOG(debug) << "Created KF PV."; + + // fill event QA histograms + if (kfparticleConfigurations.doVertexQA) { + registry.fill(HIST("QA/Event/hVtxXKF"), kfpv.GetX()); + registry.fill(HIST("QA/Event/hVtxYKF"), kfpv.GetY()); + registry.fill(HIST("QA/Event/hVtxZKF"), kfpv.GetZ()); + registry.fill(HIST("QA/Event/hVtxCovXXKF"), kfpv.GetCovariance(0)); + registry.fill(HIST("QA/Event/hVtxCovYYKF"), kfpv.GetCovariance(2)); + registry.fill(HIST("QA/Event/hVtxCovZZKF"), kfpv.GetCovariance(5)); + registry.fill(HIST("QA/Event/hVtxCovXYKF"), kfpv.GetCovariance(1)); + registry.fill(HIST("QA/Event/hVtxCovXZKF"), kfpv.GetCovariance(3)); + registry.fill(HIST("QA/Event/hVtxCovYZKF"), kfpv.GetCovariance(4)); + registry.fill(HIST("QA/Event/hVtxX"), collision.posX()); + registry.fill(HIST("QA/Event/hVtxY"), collision.posY()); + registry.fill(HIST("QA/Event/hVtxZ"), collision.posZ()); + registry.fill(HIST("QA/Event/hVtxCovXX"), collision.covXX()); + registry.fill(HIST("QA/Event/hVtxCovYY"), collision.covYY()); + registry.fill(HIST("QA/Event/hVtxCovZZ"), collision.covZZ()); + registry.fill(HIST("QA/Event/hVtxCovXY"), collision.covXY()); + registry.fill(HIST("QA/Event/hVtxCovXZ"), collision.covXZ()); + registry.fill(HIST("QA/Event/hVtxCovYZ"), collision.covYZ()); + } + for (auto& vtx3body : decay3bodys) { LOG(debug) << "Entered decay3bodys loop."; @@ -646,13 +674,9 @@ struct decay3bodyBuilder { auto trackParCovBach = getTrackParCov(trackBach); LOG(debug) << "Got all daughter tracks."; - KFPVertex kfpVertex = createKFPVertexFromCollision(collision); - KFParticle kfpv(kfpVertex); - LOG(debug) << "Created KF PV."; - bool isMatter = trackBach.sign() > 0 ? true : false; - // ---------- fill trackQA and vertexQA histograms + // ---------- fill track QA histograms ---------- if (kfparticleConfigurations.doTrackQA) { registry.fill(HIST("QA/Tracks/hTrackPosTPCNcls"), trackPos.tpcNClsFound()); registry.fill(HIST("QA/Tracks/hTrackNegTPCNcls"), trackNeg.tpcNClsFound()); @@ -676,27 +700,6 @@ struct decay3bodyBuilder { registry.fill(HIST("QA/Tracks/hTrackBachPt"), trackBach.pt()); } - if (kfparticleConfigurations.doVertexQA) { - registry.fill(HIST("QA/Event/hVtxXKF"), kfpv.GetX()); - registry.fill(HIST("QA/Event/hVtxYKF"), kfpv.GetY()); - registry.fill(HIST("QA/Event/hVtxZKF"), kfpv.GetZ()); - registry.fill(HIST("QA/Event/hVtxCovXXKF"), kfpv.GetCovariance(0)); - registry.fill(HIST("QA/Event/hVtxCovYYKF"), kfpv.GetCovariance(2)); - registry.fill(HIST("QA/Event/hVtxCovZZKF"), kfpv.GetCovariance(5)); - registry.fill(HIST("QA/Event/hVtxCovXYKF"), kfpv.GetCovariance(1)); - registry.fill(HIST("QA/Event/hVtxCovXZKF"), kfpv.GetCovariance(3)); - registry.fill(HIST("QA/Event/hVtxCovYZKF"), kfpv.GetCovariance(4)); - registry.fill(HIST("QA/Event/hVtxX"), collision.posX()); - registry.fill(HIST("QA/Event/hVtxY"), collision.posY()); - registry.fill(HIST("QA/Event/hVtxZ"), collision.posZ()); - registry.fill(HIST("QA/Event/hVtxCovXX"), collision.covXX()); - registry.fill(HIST("QA/Event/hVtxCovYY"), collision.covYY()); - registry.fill(HIST("QA/Event/hVtxCovZZ"), collision.covZZ()); - registry.fill(HIST("QA/Event/hVtxCovXY"), collision.covXY()); - registry.fill(HIST("QA/Event/hVtxCovXZ"), collision.covXZ()); - registry.fill(HIST("QA/Event/hVtxCovYZ"), collision.covYZ()); - } - // -------- STEP 1: track selection -------- // collision ID --> not correct? tracks can have different collisions, but belong to one 3prong vertex! // if (trackPos.collisionId() != trackNeg.collisionId() || trackPos.collisionId() != trackBach.collisionId() || trackNeg.collisionId() != trackBach.collisionId()) { @@ -742,16 +745,23 @@ struct decay3bodyBuilder { // TPC PID float tpcNsigmaProton; float tpcNsigmaPion; + float dEdxProton; + float dEdxPion; float tpcNsigmaDeuteron = trackBach.tpcNSigmaDe(); + float dEdxDeuteron = trackBach.tpcSignal(); if (isMatter) { // hypertriton (proton, pi-, deuteron) tpcNsigmaProton = trackPos.tpcNSigmaPr(); tpcNsigmaPion = trackNeg.tpcNSigmaPi(); + dEdxProton = trackPos.tpcSignal(); + dEdxPion = trackNeg.tpcSignal(); if (!selectTPCPID(trackPos, trackNeg, trackBach)) { continue; } } else if (!isMatter) { // anti-hypertriton (anti-proton, pi+, deuteron) tpcNsigmaProton = trackNeg.tpcNSigmaPr(); tpcNsigmaPion = trackPos.tpcNSigmaPi(); + dEdxProton = trackNeg.tpcSignal(); + dEdxPion = trackPos.tpcSignal(); if (!selectTPCPID(trackNeg, trackPos, trackBach)) { continue; } @@ -759,6 +769,13 @@ struct decay3bodyBuilder { registry.fill(HIST("hVtx3BodyCounterKFParticle"), kKfVtxTPCPID); LOG(debug) << "Basic track selections done."; + // TOF PID of deuteron (set motherhyp correctly) + double tofNSigmaDeuteron = -999; + if (trackBach.has_collision() && trackBach.hasTOF()) { + auto originalcol = trackBach.template collision_as(); + tofNSigmaDeuteron = bachelorTOFPID.GetTOFNSigma(trackBach, originalcol, collision); + } + // track DCAxy and DCAz to PV associated with decay3body o2::dataformats::VertexBase mPV; o2::dataformats::DCA mDcaInfoCovPos; @@ -1011,7 +1028,11 @@ struct decay3bodyBuilder { // daughter PID tpcNsigmaProton, tpcNsigmaPion, - tpcNsigmaDeuteron); + tpcNsigmaDeuteron, + dEdxProton, + dEdxPion, + dEdxDeuteron, + tofNSigmaDeuteron); if (kfparticleConfigurations.fillCandidateLiteTable) { kfvtx3bodydatalite( @@ -1061,7 +1082,11 @@ struct decay3bodyBuilder { // daughter PID tpcNsigmaProton, tpcNsigmaPion, - tpcNsigmaDeuteron); + tpcNsigmaDeuteron, + dEdxProton, + dEdxPion, + dEdxDeuteron, + tofNSigmaDeuteron); } LOG(debug) << "Table filled."; @@ -1085,7 +1110,7 @@ struct decay3bodyBuilder { } PROCESS_SWITCH(decay3bodyBuilder, processRun3, "Produce DCA fitter decay3body tables", true); - void processRun3withKFParticle(MyCollisions const& collisions, FullTracksExtPIDIU const&, aod::Decay3Bodys const& decay3bodys, aod::BCsWithTimestamps const&) + void processRun3withKFParticle(ColwithEvTimes const& collisions, TrackExtPIDIUwithEvTimes const&, aod::Decay3Bodys const& decay3bodys, aod::BCsWithTimestamps const&) { for (const auto& collision : collisions) { // event selection @@ -1108,7 +1133,7 @@ struct decay3bodyBuilder { // LOG(debug) << "Collision index: " << collIdx; auto Decay3BodyTable_thisCollision = decay3bodys.sliceBy(perCollision, collIdx); // LOG(debug) << "Decay3Body tables sliced per collision. Calling buildVtx3BodyDataTableKFParticle function..."; - buildVtx3BodyDataTableKFParticle(collision, Decay3BodyTable_thisCollision, bachelorcharge); + buildVtx3BodyDataTableKFParticle(collision, Decay3BodyTable_thisCollision, bachelorcharge); LOG(debug) << "End of processKFParticle."; } } diff --git a/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx b/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx index 0d560ecde46..558b34dae2a 100644 --- a/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx +++ b/PWGLF/TableProducer/Nuspex/threebodyKFTask.cxx @@ -256,6 +256,8 @@ struct threebodyKFTask { vtx3bodydata.dcatrackpostopv(), vtx3bodydata.dcatracknegtopv(), vtx3bodydata.dcatrackbachtopv(), vtx3bodydata.track0sign(), vtx3bodydata.track1sign(), vtx3bodydata.track2sign(), // proton, pion, deuteron vtx3bodydata.tpcnsigmaproton(), vtx3bodydata.tpcnsigmapion(), vtx3bodydata.tpcnsigmadeuteron(), + vtx3bodydata.tpcdedxproton(), vtx3bodydata.tpcdedxpion(), vtx3bodydata.tpcdedxdeuteron(), + vtx3bodydata.tofnsigmadeuteron(), // MC info (-1 if not matched to MC particle) genP, genPt, @@ -379,6 +381,8 @@ struct threebodyKFTask { -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, + -1, // gen information mcparticle.p(), mcparticle.pt(),