diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackCascadeExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackCascadeExtended.cxx index e76f71e7910..f07bb1cab46 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackCascadeExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackCascadeExtended.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file femtoUniversePairTaskTrackCascadeExtended.cxx -/// \brief Task for cascade QA; in the future: for cascade correlations +/// \brief Task for cascade correlations and QA /// \author Barbara Chytla, WUT Warsaw, barbara.chytla@cern.ch /// \author Shirajum Monira, WUT Warsaw, shirajum.monira@cern.ch // o2-linter: disable=name/workflow-file @@ -82,6 +82,7 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st Configurable confPhiBins{"ConfPhiBins", 29, "Number of phi bins in deta dphi"}; // o2-linter: disable=name/configurable Configurable confIsMC{"ConfIsMC", false, "Enable additional Histograms in the case of a MonteCarlo Run"}; // o2-linter: disable=name/configurable Configurable confUse3D{"ConfUse3D", false, "Enable three dimensional histogramms (to be used only for analysis with high statistics): k* vs mT vs multiplicity"}; // o2-linter: disable=name/configurable + Configurable confUseCent{"confUseCent", false, "Use centrality in place of multiplicity"}; ConfigurableAxis confVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; // o2-linter: disable=name/configurable ConfigurableAxis confTrkTempFitVarpTBins{"ConfTrkTempFitVarpTBins", {20, 0.5, 4.05}, "pT binning of the pT vs. TempFitVar plot"}; // o2-linter: disable=name/configurable ConfigurableAxis confTrkTempFitVarBins{"ConfTrkDTempFitVarBins", {300, -0.15, 0.15}, "binning of the TempFitVar in the pT vs. TempFitVar plot"}; // o2-linter: disable=name/configurable @@ -105,6 +106,9 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st FemtoUniverseParticleHisto bachHistos; FemtoUniverseParticleHisto cascQAHistos; + /// Histogramming for Event + FemtoUniverseEventHisto eventHisto; + FemtoUniverseContainer sameEventCont; FemtoUniverseContainer mixedEventCont; FemtoUniversePairCleaner pairCleaner; @@ -189,6 +193,7 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st rXiQA.add("hDcaBachtoPV", "hDcaBachtoPV", {HistType::kTH1F, {aDCAToPVAxis}}); rXiQA.add("hDcaV0toPV", "hDcaV0toPV", {HistType::kTH1F, {aDCAToPVAxis}}); + eventHisto.init(&qaRegistry); qaRegistry.add("Tracks_pos/nSigmaTPC", "; #it{p} (GeV/#it{c}); n#sigma_{TPC}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); qaRegistry.add("Tracks_pos/nSigmaTOF", "; #it{p} (GeV/#it{c}); n#sigma_{TOF}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); qaRegistry.add("Tracks_neg/nSigmaTPC", "; #it{p} (GeV/#it{c}); n#sigma_{TPC}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); @@ -273,6 +278,10 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st auto groupPartsOne = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); auto groupPartsTwo = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); + eventHisto.fillQA(col); + + const int multCol = confUseCent ? col.multV0M() : col.multNtr(); + for (const auto& part : groupPartsTwo) { if (!invMCascade(part.mLambda(), part.mAntiLambda())) continue; @@ -328,7 +337,7 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st if (!isParticleTPC(posChild, CascChildTable[confCascType1][0]) || !isParticleTPC(negChild, CascChildTable[confCascType1][1]) || !isParticleTPC(bachelor, CascChildTable[confCascType1][2])) continue; - sameEventCont.setPair(p1, p2, col.multNtr(), confUse3D, 1.0f); + sameEventCont.setPair(p1, p2, multCol, confUse3D, 1.0f); } } PROCESS_SWITCH(femtoUniversePairTaskTrackCascadeExtended, processSameEvent, "Enable processing same event for track - cascade", false); @@ -337,6 +346,10 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st { auto groupPartsTwo = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); + eventHisto.fillQA(col); + + const int multCol = confUseCent ? col.multV0M() : col.multNtr(); + for (const auto& part : groupPartsTwo) { if (!invMCascade(part.mLambda(), part.mAntiLambda())) continue; @@ -382,7 +395,7 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st if (!isParticleTPC(posChild2, CascChildTable[confCascType2][0]) || !isParticleTPC(negChild2, CascChildTable[confCascType2][1]) || !isParticleTPC(bachelor2, CascChildTable[confCascType2][2])) return; - sameEventCont.setPair(p1, p2, col.multNtr(), confUse3D, 1.0f); + sameEventCont.setPair(p1, p2, multCol, confUse3D, 1.0f); }; if (confCascType1 == confCascType2) { /// Now build the combinations for identical cascades @@ -403,6 +416,7 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st ColumnBinningPolicy colBinning{{confVtxBins, confMultBins}, true}; for (const auto& [collision1, collision2] : soa::selfCombinations(colBinning, 5, -1, cols, cols)) { + const int multCol = confUseCent ? collision1.multV0M() : collision1.multNtr(); auto groupPartsOne = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); auto groupPartsTwo = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); @@ -432,7 +446,7 @@ struct femtoUniversePairTaskTrackCascadeExtended { // o2-linter: disable=name/st continue; } - mixedEventCont.setPair(p1, p2, collision1.multNtr(), confUse3D, 1.0f); + mixedEventCont.setPair(p1, p2, multCol, confUse3D, 1.0f); } } } diff --git a/PWGLF/TableProducer/Nuspex/hypKfRecoTask.cxx b/PWGLF/TableProducer/Nuspex/hypKfRecoTask.cxx index 4dcdf0f67ec..f17b5f3c5bc 100644 --- a/PWGLF/TableProducer/Nuspex/hypKfRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/hypKfRecoTask.cxx @@ -37,7 +37,7 @@ #include "DCAFitter/DCAFitterN.h" #include "Common/DataModel/PIDResponse.h" #include "PWGLF/DataModel/LFHypernucleiKfTables.h" -#include "TRandom.h" +#include "TRandom3.h" #include "Common/DataModel/CollisionAssociationTables.h" // KFParticle @@ -171,50 +171,36 @@ constexpr double preSelectionsCascades[nCascades][nSelCas]{ //---------------------------------------------------------------------------------------------------------------- struct DaughterParticle { TString name; - int pdgCode; - double mass; - int charge; - double resolution; + int pdgCode, charge; + double mass, resolution; std::vector betheParams; - - DaughterParticle(std::string name_, int pdgCode_, double mass_, int charge_, LabeledArray bethe) + bool active; + DaughterParticle(std::string name_, int pdgCode_, double mass_, int charge_, LabeledArray bethe) : name(name_), pdgCode(pdgCode_), charge(charge_), mass(mass_), active(false) { - name = TString(name_); - pdgCode = pdgCode_; - mass = mass_; - charge = charge_; resolution = bethe.get(name, "resolution"); betheParams.clear(); for (unsigned int i = 0; i < 5; i++) betheParams.push_back(bethe.get(name, i)); } -}; // class DaughterParticle +}; // struct DaughterParticle struct HyperNucleus { TString name; int pdgCode; - double massMax; - double massMin; bool active; - std::vector daughters; - std::vector daughterTrackSigns; - - HyperNucleus(std::string name_, int pdgCode_, bool active_, std::vector daughters_, std::vector daughterTrackSigns_) + std::vector daughters, daughterTrackSigns; + HyperNucleus(std::string name_, int pdgCode_, bool active_, std::vector daughters_, std::vector daughterTrackSigns_) : pdgCode(pdgCode_), active(active_) { - name = TString(name_); - pdgCode = pdgCode_; - active = active_; - for (const int& d : daughters_) - daughters.push_back(d); - for (const int& dc : daughterTrackSigns_) - daughterTrackSigns.push_back(dc); + init(name_, daughters_, daughterTrackSigns_); } - HyperNucleus(std::string name_, int pdgCode_, bool active_, int hypDaughter, std::vector daughters_, std::vector daughterTrackSigns_) + HyperNucleus(std::string name_, int pdgCode_, bool active_, int hypDaughter, std::vector daughters_, std::vector daughterTrackSigns_) : pdgCode(pdgCode_), active(active_) { daughters.push_back(hypDaughter); + init(name_, daughters_, daughterTrackSigns_); + } + void init(std::string name_, std::vector daughters_, std::vector daughterTrackSigns_) + { name = TString(name_); - pdgCode = pdgCode_; - active = active_; for (const int& d : daughters_) daughters.push_back(d); for (const int& dc : daughterTrackSigns_) @@ -223,83 +209,78 @@ struct HyperNucleus { int getNdaughters() { return static_cast(daughters.size()); } const char* motherName() { return name.Contains("->") ? ((TString)name(0, name.First("-"))).Data() : name.Data(); } const char* daughterNames() { return name.Contains("->") ? ((TString)name(name.First("-") + 2, name.Length())).Data() : ""; } -}; // class HyperNucleus +}; // struct HyperNucleus + +struct DaughterKf { + int64_t daughterTrackId; + int hypNucId; + KFParticle daughterKfp; + float dcaToPv, dcaToPvXY, dcaToPvZ; + DaughterKf(int64_t daughterTrackId_, KFParticle daughterKfp_, std::vector vtx) : daughterTrackId(daughterTrackId_), hypNucId(-1), daughterKfp(daughterKfp_) + { + dcaToPvXY = daughterKfp.GetDistanceFromVertexXY(&vtx[0]); + dcaToPv = daughterKfp.GetDistanceFromVertex(&vtx[0]); + dcaToPvZ = std::sqrt(dcaToPv * dcaToPv - dcaToPvXY * dcaToPvXY); + } + DaughterKf(KFParticle daughterKfp_, int hypNucId_) : daughterTrackId(-999), hypNucId(hypNucId_), daughterKfp(daughterKfp_), dcaToPv(-999), dcaToPvXY(-999), dcaToPvZ(-999) {} +}; // struct DaughterKf struct HyperNucCandidate { int species; KFParticle kfp; - std::vector kfpDaughters; HyperNucCandidate* hypNucDaughter; - std::vector daughterTrackIds; + std::vector daughters; std::vector recoSV; - float devToVtx; - float dcaToVtxXY; - float dcaToVtxZ; - float chi2; - bool mcTrue; - bool isPhysPrimary; - bool isPrimaryCandidate, isSecondaryCandidate, isUsedSecondary; + float mass, px, py, pz; + float devToPvXY, dcaToPvXY, dcaToPvZ, dcaToVtxXY, dcaToVtxZ, chi2; + bool mcTrue, isPhysPrimary, isPrimaryCandidate, isSecondaryCandidate, isUsedSecondary; int64_t mcParticleId; int tableId; - - HyperNucCandidate(int species_, std::vector kfpDaughters_, std::vector daughterTrackIds_) : species(species_), hypNucDaughter(0), devToVtx(999), dcaToVtxXY(999), dcaToVtxZ(999), chi2(999), mcTrue(false), isPhysPrimary(false), isPrimaryCandidate(false), isSecondaryCandidate(false), isUsedSecondary(false), mcParticleId(-1), tableId(-1) - { - for (const auto& kfd : kfpDaughters_) - kfpDaughters.push_back(kfd); - for (const auto& dt : daughterTrackIds_) - daughterTrackIds.push_back(dt); - init(); - } - HyperNucCandidate(int species_, HyperNucCandidate* hypNucDaughter_, std::vector kfpDaughters_, std::vector daughterTrackIds_) : species(species_), hypNucDaughter(hypNucDaughter_), devToVtx(999), dcaToVtxXY(999), dcaToVtxZ(999), chi2(999), mcTrue(false), isPhysPrimary(false), isPrimaryCandidate(false), isSecondaryCandidate(false), isUsedSecondary(false), mcParticleId(-1), tableId(-1) - { - for (const auto& kfd : kfpDaughters_) - kfpDaughters.push_back(kfd); - for (const auto& dt : daughterTrackIds_) - daughterTrackIds.push_back(dt); - init(); - } - - void init() + HyperNucCandidate(int species_, HyperNucCandidate* hypNucDaughter_, std::vector daughters_) : species(species_), hypNucDaughter(hypNucDaughter_), devToPvXY(999), dcaToPvXY(999), dcaToPvZ(999), dcaToVtxXY(999), dcaToVtxZ(999), chi2(999), mcTrue(false), isPhysPrimary(false), isPrimaryCandidate(false), isSecondaryCandidate(false), isUsedSecondary(false), mcParticleId(-1), tableId(-1) { + for (const auto& d : daughters_) + daughters.push_back(d); kfp.SetConstructMethod(2); - for (size_t i = 0; i < kfpDaughters.size(); i++) - kfp.AddDaughter(kfpDaughters.at(i)); + for (size_t i = 0; i < daughters.size(); i++) + kfp.AddDaughter(daughters.at(i)->daughterKfp); kfp.TransportToDecayVertex(); chi2 = kfp.GetChi2() / kfp.GetNDF(); recoSV.clear(); recoSV.push_back(kfp.GetX()); recoSV.push_back(kfp.GetY()); recoSV.push_back(kfp.GetZ()); + mass = kfp.GetMass(); + px = kfp.GetPx(); + py = kfp.GetPy(); + pz = kfp.GetPz(); } - bool checkKfp() - { - if (kfp.GetMass() == 0) - return false; - if (std::isnan(kfp.GetMass())) - return false; - return true; - } - int getDaughterTableId() + std::vector daughterTrackIds() { - if (hypNucDaughter) - return hypNucDaughter->tableId; - return -1; + std::vector trackIds; + for (const auto& daughter : daughters) { + const auto& id = daughter->daughterTrackId; + if (id >= 0) + trackIds.push_back(id); + } + return trackIds; } + bool checkKfp() { return mass != 0 && !std::isnan(mass); } + int getDaughterTableId() { return hypNucDaughter ? hypNucDaughter->tableId : -1; } bool isCascade() { return hypNucDaughter != 0; } int getSign() { if (kfp.GetQ() == 0) - return kfpDaughters.front().GetQ() / std::abs(kfpDaughters.front().GetQ()); + return daughters.front()->daughterKfp.GetQ() / std::abs(daughters.front()->daughterKfp.GetQ()); return kfp.GetQ() / std::abs(kfp.GetQ()); } - int getNdaughters() { return static_cast(kfpDaughters.size()); } + int getNdaughters() { return static_cast(daughters.size()); } float getDcaTracks() { return getNdaughters() == 2 ? getDcaTracks2() : getMaxDcaToSv(); } - float getDcaTracks2() { return kfpDaughters.at(0).GetDistanceFromParticle(kfpDaughters.at(1)); } + float getDcaTracks2() { return daughters.at(0)->daughterKfp.GetDistanceFromParticle(daughters.at(1)->daughterKfp); } float getMaxDcaToSv() { float maxDca = std::numeric_limits::lowest(); - for (const auto& daughter : kfpDaughters) { - float dca = daughter.GetDistanceFromVertex(&recoSV[0]); + for (size_t i = 0; i < daughters.size(); i++) { + float dca = daughters.at(i)->daughterKfp.GetDistanceFromVertex(&recoSV[0]); if (dca > maxDca) maxDca = dca; } @@ -309,7 +290,7 @@ struct HyperNucCandidate { double getCpa(std::vector vtx) { kfp.TransportToDecayVertex(); - return RecoDecay::cpa(std::array{vtx[0], vtx[1], vtx[2]}, std::array{recoSV[0], recoSV[1], recoSV[2]}, std::array{kfp.GetPx(), kfp.GetPy(), kfp.GetPz()}); + return RecoDecay::cpa(std::array{vtx[0], vtx[1], vtx[2]}, std::array{recoSV[0], recoSV[1], recoSV[2]}, std::array{px, py, pz}); ; } float getCt(std::vector vtx) @@ -319,7 +300,13 @@ struct HyperNucCandidate { float tmp = recoSV.at(i) - vtx.at(i); dl += (tmp * tmp); } - return std::sqrt(dl) * kfp.GetMass() / kfp.GetP(); + return std::sqrt(dl) * mass / std::sqrt(px * px + py * py + pz * pz); + } + void getDaughterPosMom(int daughter, std::vector& posMom) + { + auto kfpDaughter = daughters.at(daughter)->daughterKfp; + kfpDaughter.TransportToPoint(&recoSV[0]); + posMom.assign({kfpDaughter.GetX(), kfpDaughter.GetY(), kfpDaughter.GetZ(), kfpDaughter.GetPx(), kfpDaughter.GetPy(), kfpDaughter.GetPz()}); } float getDcaMotherToVtxXY(std::vector vtx) { return kfp.GetDistanceFromVertexXY(&vtx[0]); } float getDcaMotherToVtxZ(std::vector vtx) @@ -327,15 +314,17 @@ struct HyperNucCandidate { kfp.TransportToPoint(&vtx[0]); return kfp.GetZ() - vtx[2]; } - void getDaughterPosMom(int daughter, std::vector& posMom) + void calcDcaToVtx(KFPVertex& vtx) { - kfpDaughters.at(daughter).TransportToPoint(&recoSV[0]); - posMom.assign({kfpDaughters.at(daughter).GetX(), kfpDaughters.at(daughter).GetY(), kfpDaughters.at(daughter).GetZ(), kfpDaughters.at(daughter).GetPx(), kfpDaughters.at(daughter).GetPy(), kfpDaughters.at(daughter).GetPz()}); + if (devToPvXY != 999) + return; + devToPvXY = kfp.GetDeviationFromVertexXY(vtx); + dcaToPvXY = kfp.GetDistanceFromVertexXY(vtx); + kfp.TransportToVertex(vtx); + dcaToPvZ = kfp.GetZ() - vtx.GetZ(); } - void calcDevToVtx(KFPVertex& vtx) { devToVtx = kfp.GetDeviationFromVertexXY(vtx); } - void calcDevToVtx(HyperNucCandidate& cand) + void calcDcaToVtx(HyperNucCandidate& cand) { - devToVtx = kfp.GetDeviationFromParticleXY(cand.kfp); dcaToVtxXY = getDcaMotherToVtxXY(cand.recoSV); dcaToVtxZ = getDcaMotherToVtxZ(cand.recoSV); } @@ -343,26 +332,12 @@ struct HyperNucCandidate { { KFParticle subDaughter; subDaughter.SetConstructMethod(2); - subDaughter.AddDaughter(kfpDaughters.at(d1)); - subDaughter.AddDaughter(kfpDaughters.at(d2)); + subDaughter.AddDaughter(daughters.at(d1)->daughterKfp); + subDaughter.AddDaughter(daughters.at(d2)->daughterKfp); subDaughter.TransportToDecayVertex(); return subDaughter.GetMass(); } - KFParticle getDaughterTrackKfp(int track) - { - return kfpDaughters.at(track + (isCascade() ? 1 : 0)); - } - float getDcaTrackToVtxXY(int track, std::vector vtx) - { - return getDaughterTrackKfp(track).GetDistanceFromVertexXY(&vtx[0]); - } - float getDcaTrackToVtxZ(int track, std::vector vtx) - { - auto dKfp = getDaughterTrackKfp(track); - dKfp.TransportToPoint(&vtx[0]); - return dKfp.GetZ() - vtx[2]; - } -}; // class HyperNucCandidate +}; // struct HyperNucCandidate struct IndexPairs { std::vector> pairs; @@ -379,7 +354,7 @@ struct IndexPairs { } return false; } -}; // class IndexPairs +}; // struct IndexPairs struct McCollInfo { bool hasRecoColl; @@ -387,7 +362,7 @@ struct McCollInfo { bool hasRecoParticle; int tableIndex; McCollInfo() : hasRecoColl(false), passedEvSel(false), hasRecoParticle(false), tableIndex(-1) {} -}; // class McCollInfo +}; // struct McCollInfo //---------------------------------------------------------------------------------------------------------------- std::vector> hDeDx; @@ -396,7 +371,7 @@ std::vector> hInvMass; //---------------------------------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------------------------------- -struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] +struct HypKfRecoTask { Produces outputMcCollisionTable; Produces outputMcParticleTable; @@ -447,22 +422,18 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] std::vector daughterParticles; std::vector> foundDaughters; - std::vector> singleHyperNucCandidates; // hypernuclei candidates - std::vector> cascadeHyperNucCandidates; // cascade candidates - std::vector singleHyperNuclei; - std::vector cascadeHyperNuclei; - std::vector primVtx; - std::vector cents; + std::vector> foundDaughterKfs, hypNucDaughterKfs; + std::vector> singleHyperNucCandidates, cascadeHyperNucCandidates; + std::vector singleHyperNuclei, cascadeHyperNuclei; + std::vector primVtx, cents; std::vector mcCollInfos; - IndexPairs trackIndices; - IndexPairs mcPartIndices; + IndexPairs trackIndices, mcPartIndices; KFPVertex kfPrimVtx; - bool collHasCandidate, collHasMcTrueCandidate; - bool collPassedEvSel; + bool collHasCandidate, collHasMcTrueCandidate, collPassedEvSel, activeCascade; int64_t mcCollTableIndex; int mRunNumber, occupancy; float dBz; - TRandom rand; + TRandom3 rand; //---------------------------------------------------------------------------------------------------------------- void init(InitContext const&) @@ -482,8 +453,11 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] for (unsigned int i = 0; i < nHyperNuclei; i++) { // create hypernuclei singleHyperNuclei.push_back(HyperNucleus(hyperNucNames.at(i), cfgHyperNucPdg->get(i, 0u), cfgHyperNucsActive->get(i, 0u), getDaughterVec(i, cfgHyperNucDaughters), getDaughterSignVec(i, cfgHyperNucSigns))); } + activeCascade = false; for (unsigned int i = 0; i < nCascades; i++) { // create cascades cascadeHyperNuclei.push_back(HyperNucleus(cascadeNames.at(i), cfgCascadesPdg->get(i, 0u), cfgCascadesActive->get(i, 0u), getHypDaughterVec(i, cfgCascadeHypDaughter), getDaughterVec(i, cfgCascadeDaughters), getDaughterSignVec(i, cfgCascadeSigns))); + if (cfgCascadesActive->get(i, 0u)) + activeCascade = true; } // define histogram axes const AxisSpec axisMagField{10, -10., 10., "magnetic field"}; @@ -543,10 +517,6 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] if (std::abs(getTPCnSigma(track, daughterParticles.at(i))) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) continue; filldedx(track, i); - if (std::abs(track.dcaXY()) < cfgTrackPIDsettings->get(i, "minDcaToPvXY")) - continue; - if (std::abs(track.dcaZ()) < cfgTrackPIDsettings->get(i, "minDcaToPvZ")) - continue; if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize")) continue; if (getMeanItsClsSize(track) > cfgTrackPIDsettings->get(i, "maxITSclsSize")) @@ -563,29 +533,83 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] //---------------------------------------------------------------------------------------------------------------- void checkMCTrueTracks(aod::McTrackLabels const& trackLabels, aod::McParticles const&) { + std::vector activePdgs; + std::vector*> hypNucVectors = {&singleHyperNuclei, &cascadeHyperNuclei}; + for (int vec = 0; vec < 2; vec++) { + for (size_t hyperNucIter = 0; hyperNucIter < hypNucVectors.at(vec)->size(); hyperNucIter++) { + HyperNucleus* hyperNuc = &(hypNucVectors.at(vec)->at(hyperNucIter)); + if (!hyperNuc->active) + continue; + activePdgs.push_back(std::abs(hyperNuc->pdgCode)); + } + } for (int i = 0; i < nDaughterParticles; i++) { auto& daughterVec = foundDaughters.at(i); - for (auto it = daughterVec.begin(); it < daughterVec.end(); it++) { - bool mcTrue = true; + if (!daughterVec.size()) + continue; + for (auto it = daughterVec.end() - 1; it >= daughterVec.begin(); it--) { const auto& mcLab = trackLabels.rawIteratorAt(*it); if (!mcLab.has_mcParticle()) { - mcTrue = false; - } else { - const auto& mcPart = mcLab.mcParticle_as(); - if (std::abs(mcPart.pdgCode()) != daughterParticles.at(i).pdgCode) { - mcTrue = false; - } - if (!mcPart.has_mothers()) { - mcTrue = false; + daughterVec.erase(it); + continue; + } + const auto& mcPart = mcLab.mcParticle_as(); + if (std::abs(mcPart.pdgCode()) != daughterParticles.at(i).pdgCode) { + daughterVec.erase(it); + continue; + } + if (!mcPart.has_mothers()) { + daughterVec.erase(it); + continue; + } + bool isDaughter = false; + for (const auto& mother : mcPart.mothers_as()) { + if (std::find(activePdgs.begin(), activePdgs.end(), std::abs(mother.pdgCode())) != activePdgs.end()) { + isDaughter = true; } } - if (!mcTrue) { + if (!isDaughter) { daughterVec.erase(it); + continue; } } } } //---------------------------------------------------------------------------------------------------------------- + void createKFDaughters(TracksFull const& tracks) + { + std::vector*> hypNucVectors = {&singleHyperNuclei, &cascadeHyperNuclei}; + for (size_t vec = 0; vec < 2; vec++) { + for (const auto& hyperNuc : *(hypNucVectors.at(vec))) { + if (!hyperNuc.active) + continue; + for (size_t i = vec; i < hyperNuc.daughters.size(); i++) { + if (foundDaughters.at(hyperNuc.daughters.at(i)).size() > 0) + daughterParticles.at(hyperNuc.daughters.at(i)).active = true; + else + break; + } + } + } + for (size_t daughterCount = 0; daughterCount < daughterParticles.size(); daughterCount++) { + const auto& daughterParticle = daughterParticles.at(daughterCount); + if (!daughterParticle.active) + continue; + const auto& daughterMass = daughterParticle.mass; + const auto& daughterCharge = daughterParticle.charge; + for (const auto& daughterId : foundDaughters.at(daughterCount)) { + const auto& daughterTrack = tracks.rawIteratorAt(daughterId); + DaughterKf daughter(daughterId, createKFParticle(daughterTrack, daughterMass, daughterCharge), primVtx); + if (std::abs(daughter.dcaToPvXY) < cfgTrackPIDsettings->get(daughterCount, "minDcaToPvXY")) + continue; + if (std::abs(daughter.dcaToPvZ) < cfgTrackPIDsettings->get(daughterCount, "minDcaToPvZ")) + continue; + foundDaughterKfs.at(daughterCount).push_back(daughter); + } + } + } + //---------------------------------------------------------------------------------------------------------------- + void createKFHypernuclei(TracksFull const& tracks) { // loop over all hypernuclei that are to be reconstructed @@ -594,80 +618,69 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] if (!hyperNuc->active) continue; int nDaughters = hyperNuc->getNdaughters(); - - std::vector::iterator> it; + std::vector::iterator> it; int nCombinations = 1; for (int i = 0; i < nDaughters; i++) { - nCombinations *= foundDaughters.at(hyperNuc->daughters.at(i)).size(); - it.push_back(foundDaughters.at(hyperNuc->daughters.at(i)).begin()); + nCombinations *= foundDaughterKfs.at(hyperNuc->daughters.at(i)).size(); + it.push_back(foundDaughterKfs.at(hyperNuc->daughters.at(i)).begin()); } if (!nCombinations) continue; - while (it[0] != foundDaughters.at(hyperNuc->daughters.at(0)).end()) { - + const float reduceFactor = cfgReduce->get(hyperNucIter, 0u); + const float minMassPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "minMass"); + const float maxMassPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "maxMass"); + const float minCtPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "minCt"); + const float maxCtPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "maxCt"); + const float minCosPaPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "minCosPa"); + const float maxDcaTracksPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "maxDcaTracks"); + const float maxDcaMotherToPvXYPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "maxDcaMotherToPvXY"); + const float maxDcaMotherToPvZPrim = cfgPreSelectionsPrimaries->get(hyperNucIter, "maxDcaMotherToPvZ"); + const float minMassSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "minMass"); + const float maxMassSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "maxMass"); + const float minCtSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "minCt"); + const float maxCtSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "maxCt"); + const float maxDcaTracksSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "maxDcaTracks"); + while (it[0] != foundDaughterKfs.at(hyperNuc->daughters.at(0)).end()) { // check for correct signs, avoid double usage of tracks bool passedChecks = true; int checkSign = 0; std::vector vec; for (int i = 0; i < nDaughters; i++) { - const auto& daughterTrack = tracks.rawIteratorAt(*(it[i])); + const auto& daughterTrack = tracks.rawIteratorAt(it[i]->daughterTrackId); if (!i) checkSign = daughterTrack.sign(); - if (daughterTrack.sign() != checkSign * hyperNuc->daughterTrackSigns.at(i) || std::find(vec.begin(), vec.end(), *it[i]) != vec.end()) { + if (daughterTrack.sign() != checkSign * hyperNuc->daughterTrackSigns.at(i) || std::find(vec.begin(), vec.end(), it[i]->daughterTrackId) != vec.end()) { passedChecks = false; break; } - vec.push_back(*it[i]); + vec.push_back(it[i]->daughterTrackId); } - if (passedChecks && rand.Rndm() <= cfgReduce->get((unsigned int)hyperNucIter, 0u)) { - // create daugther KFParticles - std::vector daughterIds; - std::vector daughterKfps; + if (passedChecks && rand.Rndm() <= reduceFactor) { + std::vector daughters; for (int i = 0; i < nDaughters; i++) { - const auto& daughterTrack = tracks.rawIteratorAt(*(it[i])); - daughterIds.push_back(*(it[i])); - auto daughterMass = daughterParticles.at(hyperNuc->daughters.at(i)).mass; - auto daughterCharge = daughterParticles.at(hyperNuc->daughters.at(i)).charge; - daughterKfps.push_back(createKFParticle(daughterTrack, daughterMass, daughterCharge)); + daughters.push_back(&(*it[i])); } - - HyperNucCandidate candidate(hyperNucIter, daughterKfps, daughterIds); - bool isPrimCandidate = true, isSecCandidate = true; + HyperNucCandidate candidate(hyperNucIter, static_cast(0), daughters); + // check preselections if (candidate.checkKfp()) { - // apply pre selections - candidate.calcDevToVtx(kfPrimVtx); - if (candidate.kfp.GetMass() < cfgPreSelectionsPrimaries->get(hyperNucIter, "minMass") || candidate.kfp.GetMass() > cfgPreSelectionsPrimaries->get(hyperNucIter, "maxMass")) - isPrimCandidate = false; - if (candidate.getDcaTracks() > cfgPreSelectionsPrimaries->get(hyperNucIter, "maxDcaTracks")) - isPrimCandidate = false; - if (candidate.getCt(primVtx) < cfgPreSelectionsPrimaries->get(hyperNucIter, "minCt") || candidate.getCt(primVtx) > cfgPreSelectionsPrimaries->get(hyperNucIter, "maxCt")) - isPrimCandidate = false; - if (candidate.getCpa(primVtx) < cfgPreSelectionsPrimaries->get(hyperNucIter, "minCosPa")) - isPrimCandidate = false; - if (std::abs(candidate.getDcaMotherToVtxXY(primVtx)) > cfgPreSelectionsPrimaries->get(hyperNucIter, "maxDcaMotherToPvXY")) - isPrimCandidate = false; - if (std::abs(candidate.getDcaMotherToVtxZ(primVtx)) > cfgPreSelectionsPrimaries->get(hyperNucIter, "maxDcaMotherToPvZ")) - isPrimCandidate = false; - if (isPrimCandidate) { - candidate.isPrimaryCandidate = true; - collHasCandidate = true; + if (candidate.mass <= maxMassPrim && candidate.mass >= minMassPrim && candidate.getDcaTracks() <= maxDcaTracksPrim && candidate.getCt(primVtx) <= maxCtPrim && candidate.getCt(primVtx) >= minCtPrim && candidate.getCpa(primVtx) >= minCosPaPrim) { + candidate.calcDcaToVtx(kfPrimVtx); + if (std::abs(candidate.dcaToPvXY) <= maxDcaMotherToPvXYPrim && std::abs(candidate.dcaToPvZ) <= maxDcaMotherToPvZPrim) { + candidate.isPrimaryCandidate = true; + collHasCandidate = true; + } } - if (candidate.kfp.GetMass() < cfgPreSelectionsSecondaries->get(hyperNucIter, "minMass") || candidate.kfp.GetMass() > cfgPreSelectionsSecondaries->get(hyperNucIter, "maxMass")) - isSecCandidate = false; - if (candidate.getDcaTracks() > cfgPreSelectionsSecondaries->get(hyperNucIter, "maxDcaTracks")) - isSecCandidate = false; - if (candidate.getCt(primVtx) < cfgPreSelectionsSecondaries->get(hyperNucIter, "minCt") || candidate.getCt(primVtx) > cfgPreSelectionsSecondaries->get(hyperNucIter, "maxCt")) - isSecCandidate = false; - if (isSecCandidate) { + if (activeCascade && candidate.mass <= maxMassSec && candidate.mass >= minMassSec && candidate.getDcaTracks() <= maxDcaTracksSec && candidate.getCt(primVtx) <= maxCtSec && candidate.getCt(primVtx) >= minCtSec) { + candidate.calcDcaToVtx(kfPrimVtx); candidate.isSecondaryCandidate = true; } - if (isPrimCandidate || isSecCandidate) + if (candidate.isPrimaryCandidate || candidate.isSecondaryCandidate) singleHyperNucCandidates.at(hyperNucIter).push_back(candidate); } } it[nDaughters - 1]++; - for (int i = nDaughters - 1; i && it[i] == foundDaughters.at(hyperNuc->daughters.at(i)).end(); i--) { - it[i] = foundDaughters.at(hyperNuc->daughters.at(i)).begin(); + for (int i = nDaughters - 1; i && it[i] == foundDaughterKfs.at(hyperNuc->daughters.at(i)).end(); i--) { + it[i] = foundDaughterKfs.at(hyperNuc->daughters.at(i)).begin(); it[i - 1]++; } } @@ -676,7 +689,6 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] //---------------------------------------------------------------------------------------------------------------- void createKFCascades(TracksFull const& tracks) { - // loop over all cascade hypernuclei that are to be reconstructed for (size_t hyperNucIter = 0; hyperNucIter < cascadeHyperNuclei.size(); hyperNucIter++) { HyperNucleus* hyperNuc = &(cascadeHyperNuclei.at(hyperNucIter)); @@ -685,86 +697,75 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] int nDaughters = hyperNuc->getNdaughters(); int nHypNucDaughters = singleHyperNucCandidates.at(hyperNuc->daughters.at(0)).size(); - std::vector vecHypNucDaughers; for (int64_t i = 0; i < static_cast(nHypNucDaughters); i++) { - vecHypNucDaughers.push_back(i); + if (singleHyperNucCandidates.at(hyperNuc->daughters.at(0)).at(i).isSecondaryCandidate) { + auto hypNucDaughter = &(singleHyperNucCandidates.at(hyperNuc->daughters.at(0)).at(i)); + hypNucDaughterKfs.at(hyperNucIter).push_back(DaughterKf(hypNucDaughter->kfp, i)); + } } - - std::vector::iterator> it; - int nCombinations = 1; - nCombinations *= nHypNucDaughters; - it.push_back(vecHypNucDaughers.begin()); + int nCombinations = hypNucDaughterKfs.at(hyperNucIter).size(); + std::vector::iterator> it; + it.push_back(hypNucDaughterKfs.at(hyperNucIter).begin()); for (int i = 1; i < nDaughters; i++) { - nCombinations *= foundDaughters.at(hyperNuc->daughters.at(i)).size(); - it.push_back(foundDaughters.at(hyperNuc->daughters.at(i)).begin()); + nCombinations *= foundDaughterKfs.at(hyperNuc->daughters.at(i)).size(); + it.push_back(foundDaughterKfs.at(hyperNuc->daughters.at(i)).begin()); } if (!nCombinations) continue; - while (it[0] != vecHypNucDaughers.end()) { - std::vector daughterIds; - std::vector daughterKfps; - + const float minMassCas = cfgPreSelectionsCascades->get(hyperNucIter, "minMass"); + const float maxMassCas = cfgPreSelectionsCascades->get(hyperNucIter, "maxMass"); + const float minCtCas = cfgPreSelectionsCascades->get(hyperNucIter, "minCt"); + const float maxCtCas = cfgPreSelectionsCascades->get(hyperNucIter, "maxCt"); + const float minCosPaCas = cfgPreSelectionsCascades->get(hyperNucIter, "minCosPa"); + const float maxDcaTracksCas = cfgPreSelectionsCascades->get(hyperNucIter, "maxDcaTracks"); + const float maxDcaMotherToPvXYCas = cfgPreSelectionsCascades->get(hyperNucIter, "maxDcaMotherToPvXY"); + const float maxDcaMotherToPvZCas = cfgPreSelectionsCascades->get(hyperNucIter, "maxDcaMotherToPvZ"); + const float minCtSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "minCt"); + const float maxCtSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "maxCt"); + const float minCosPaSvSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "minCosPaSv"); + const float maxDcaMotherToSvXYSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "maxDcaMotherToSvXY"); + const float maxDcaMotherToSvZSec = cfgPreSelectionsSecondaries->get(hyperNucIter, "maxDcaMotherToSvZ"); + + while (it[0] != hypNucDaughterKfs.at(hyperNucIter).end()) { // select hypernuclei daughter KFParticle - auto hypNucDaughter = &(singleHyperNucCandidates.at(hyperNuc->daughters.at(0)).at(*it[0])); + auto hypNucDaughter = &(singleHyperNucCandidates.at(hyperNuc->daughters.at(0)).at(it[0]->hypNucId)); // check for correct signs int checkSign = hypNucDaughter->getSign(); bool passedChecks = true; - std::vector vec = hypNucDaughter->daughterTrackIds; + std::vector vec = hypNucDaughter->daughterTrackIds(); for (int i = 1; i < nDaughters; i++) { - const auto& daughterTrack = tracks.rawIteratorAt(*(it[i])); - if (!i) - checkSign = daughterTrack.sign(); - if (daughterTrack.sign() != checkSign * hyperNuc->daughterTrackSigns.at(i) || std::find(vec.begin(), vec.end(), *it[i]) != vec.end()) { + const auto& daughterTrack = tracks.rawIteratorAt(it[i]->daughterTrackId); + if (daughterTrack.sign() != checkSign * hyperNuc->daughterTrackSigns.at(i) || std::find(vec.begin(), vec.end(), it[i]->daughterTrackId) != vec.end()) { passedChecks = false; break; } - vec.push_back(*it[i]); + vec.push_back(it[i]->daughterTrackId); } - if (passedChecks && hypNucDaughter->isSecondaryCandidate) { - daughterKfps.push_back(hypNucDaughter->kfp); + if (passedChecks) { + std::vector daughters; + daughters.push_back(&(*it[0])); for (int i = 1; i < nDaughters; i++) { - daughterIds.push_back(*(it[i])); - const auto& daughterTrack = tracks.rawIteratorAt(*(it[i])); - auto daughterMass = daughterParticles.at(hyperNuc->daughters.at(i)).mass; - auto daughterCharge = daughterParticles.at(hyperNuc->daughters.at(i)).charge; - daughterKfps.push_back(createKFParticle(daughterTrack, daughterMass, daughterCharge)); + daughters.push_back(&(*it[i])); } - - HyperNucCandidate candidate(hyperNucIter, hypNucDaughter, daughterKfps, daughterIds); + HyperNucCandidate candidate(hyperNucIter, hypNucDaughter, daughters); if (candidate.checkKfp()) { - hypNucDaughter->calcDevToVtx(candidate); - bool isCandidate = true; - // apply pre selections for hypernucleus daughter - if (hypNucDaughter->getCpa(candidate.recoSV) < cfgPreSelectionsSecondaries->get(hyperNuc->daughters.at(0), "minCosPaSv")) - isCandidate = false; - if (hypNucDaughter->getDcaMotherToVtxXY(candidate.recoSV) > cfgPreSelectionsSecondaries->get(hyperNuc->daughters.at(0), "maxDcaMotherToSvXY")) - isCandidate = false; - if (hypNucDaughter->getDcaMotherToVtxZ(candidate.recoSV) > cfgPreSelectionsSecondaries->get(hyperNuc->daughters.at(0), "maxDcaMotherToSvZ")) - isCandidate = false; - // apply pre selections for cascade - if (candidate.kfp.GetMass() < cfgPreSelectionsCascades->get(hyperNucIter, "minMass") || candidate.kfp.GetMass() > cfgPreSelectionsCascades->get(hyperNucIter, "maxMass")) - isCandidate = false; - if (candidate.getDcaTracks() > cfgPreSelectionsCascades->get(hyperNucIter, "maxDcaTracks")) - isCandidate = false; - if (candidate.getCt(primVtx) < cfgPreSelectionsCascades->get(hyperNucIter, "minCt") || candidate.getCt(primVtx) > cfgPreSelectionsCascades->get(hyperNucIter, "maxCt")) - isCandidate = false; - if (candidate.getCpa(primVtx) < cfgPreSelectionsCascades->get(hyperNucIter, "minCosPa")) - isCandidate = false; - if (std::abs(candidate.getDcaMotherToVtxXY(primVtx)) > cfgPreSelectionsCascades->get(hyperNucIter, "maxDcaMotherToPvXY")) - isCandidate = false; - if (std::abs(candidate.getDcaMotherToVtxZ(primVtx)) > cfgPreSelectionsCascades->get(hyperNucIter, "maxDcaMotherToPvZ")) - isCandidate = false; - - if (isCandidate) { - collHasCandidate = true; - hypNucDaughter->isUsedSecondary = true; - cascadeHyperNucCandidates.at(hyperNucIter).push_back(candidate); + // preselections for cascade and hypernucleus daughter + if (candidate.mass <= maxMassCas && candidate.mass >= minMassCas && candidate.getDcaTracks() <= maxDcaTracksCas && candidate.getCt(primVtx) >= minCtCas && candidate.getCt(primVtx) <= maxCtCas && hypNucDaughter->getCt(candidate.recoSV) >= minCtSec && hypNucDaughter->getCt(candidate.recoSV) <= maxCtSec && candidate.getCpa(primVtx) >= minCosPaCas && hypNucDaughter->getCpa(candidate.recoSV) >= minCosPaSvSec) { + candidate.calcDcaToVtx(kfPrimVtx); + if (std::abs(candidate.dcaToPvXY) <= maxDcaMotherToPvXYCas && std::abs(candidate.dcaToPvZ) <= maxDcaMotherToPvZCas) { + hypNucDaughter->calcDcaToVtx(candidate); + if (hypNucDaughter->dcaToVtxXY <= maxDcaMotherToSvXYSec && hypNucDaughter->dcaToVtxZ <= maxDcaMotherToSvZSec) { + collHasCandidate = true; + hypNucDaughter->isUsedSecondary = true; + cascadeHyperNucCandidates.at(hyperNucIter).push_back(candidate); + } + } } } } it[nDaughters - 1]++; - for (int i = nDaughters - 1; i && it[i] == foundDaughters.at(hyperNuc->daughters.at(i)).end(); i--) { - it[i] = foundDaughters.at(hyperNuc->daughters.at(i)).begin(); + for (int i = nDaughters - 1; i && it[i] == foundDaughterKfs.at(hyperNuc->daughters.at(i)).end(); i--) { + it[i] = foundDaughterKfs.at(hyperNuc->daughters.at(i)).begin(); it[i - 1]++; } } @@ -774,7 +775,6 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] void createMCinfo(aod::McTrackLabels const& trackLabels, aod::McCollisionLabels const&, aod::McParticles const& particlesMC, aod::McCollisions const&, bool cascadesOnly = false) { // check for mcTrue: single (primary & cascade daughter) and cascade hypernuclei - std::vector*> hypNucVectors = {&singleHyperNuclei, &cascadeHyperNuclei}; std::vector>*> candidateVectors = {&singleHyperNucCandidates, &cascadeHyperNucCandidates}; const int nVecs = candidateVectors.size(); @@ -794,16 +794,16 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] const auto& mcPart = particlesMC.rawIteratorAt(hypCand.hypNucDaughter->mcParticleId); if (!mcPart.has_mothers()) continue; - daughterCount++; for (const auto& mother : mcPart.mothers_as()) { if (mother.pdgCode() == hyperNuc->pdgCode * hypCand.getSign()) { motherIds.push_back(mother.globalIndex()); break; } } + daughterCount++; } - for (const auto& daughter : hypCand.daughterTrackIds) { - const auto& mcLab = trackLabels.rawIteratorAt(daughter); + for (const auto& daughter : hypCand.daughters) { + const auto& mcLab = trackLabels.rawIteratorAt(daughter->daughterTrackId); if (!mcLab.has_mcParticle()) continue; const auto& mcPart = mcLab.mcParticle_as(); @@ -861,10 +861,11 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] continue; if (saveOnlyMcTrue && !hypCand.mcTrue) continue; - hInvMass[vec * nHyperNuclei + hyperNucIter]->Fill(hypCand.kfp.GetMass()); + hInvMass[vec * nHyperNuclei + hyperNucIter]->Fill(hypCand.mass); std::vector vecDaugtherTracks, vecAddons, vecSubDaughters; int daughterCount = 0; - for (const auto& daughterTrackId : hypCand.daughterTrackIds) { + for (const auto& daughter : hypCand.daughters) { + const auto& daughterTrackId = daughter->daughterTrackId; int trackTableId; if (!trackIndices.getIndex(daughterTrackId, trackTableId)) { auto daught = hyperNuc->daughters.at(daughterCount); @@ -872,7 +873,7 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] outputTrackTable( hyperNuc->daughters.at(daughterCount) * track.sign(), track.pt(), track.eta(), track.phi(), - hypCand.getDcaTrackToVtxXY(daughterCount, primVtx), hypCand.getDcaTrackToVtxZ(daughterCount, primVtx), + daughter->dcaToPvXY, daughter->dcaToPvZ, track.tpcNClsFound(), track.tpcChi2NCl(), track.itsClusterSizes(), track.itsChi2NCl(), getRigidity(track), track.tpcSignal(), getTPCnSigma(track, daughterParticles.at(daught)), @@ -908,10 +909,10 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] mcPartIndices.getIndex(hypCand.mcParticleId, mcPartTableId) ? mcPartTableId : -1, outputCollisionTable.lastIndex(), vecDaugtherTracks, vecAddons, hypCand.getDaughterTableId(), vecSubDaughters, (vec * nHyperNuclei + hyperNucIter + 1) * hypCand.getSign(), - hypCand.isPrimaryCandidate, hypCand.kfp.GetMass(), - hypCand.kfp.GetPx(), hypCand.kfp.GetPy(), hypCand.kfp.GetPz(), - hypCand.getDcaMotherToVtxXY(primVtx), hypCand.getDcaMotherToVtxZ(primVtx), - hypCand.devToVtx, hypCand.dcaToVtxXY, hypCand.dcaToVtxZ, hypCand.chi2, + hypCand.isPrimaryCandidate, hypCand.mass, + hypCand.px, hypCand.py, hypCand.pz, + hypCand.dcaToPvXY, hypCand.dcaToPvZ, hypCand.devToPvXY, + hypCand.dcaToVtxXY, hypCand.dcaToVtxZ, hypCand.chi2, hypCand.recoSV.at(0), hypCand.recoSV.at(1), hypCand.recoSV.at(2)); hypCand.tableId = outputHypNucTable.lastIndex(); } @@ -989,6 +990,7 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] findDaughterParticles(tracksByColl, tracks); if (cfgSaveOnlyMcTrue) checkMCTrueTracks(trackLabelsMC, particlesMC); + createKFDaughters(tracks); createKFHypernuclei(tracks); createMCinfo(trackLabelsMC, collLabels, particlesMC, mcColls); createKFCascades(tracks); @@ -1013,7 +1015,7 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] fillTree(tracks, cfgSaveOnlyMcTrue); } } - PROCESS_SWITCH(hypKfRecoTask, processMC, "MC analysis", false); + PROCESS_SWITCH(HypKfRecoTask, processMC, "MC analysis", false); //---------------------------------------------------------------------------------------------------------------- void processData(CollisionsFull const& collisions, TracksFull const& tracks, aod::BCsWithTimestamps const&, aod::TrackAssoc const& tracksColl) { @@ -1027,6 +1029,7 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] const uint64_t collIdx = collision.globalIndex(); auto tracksByColl = tracksColl.sliceBy(perCollision, collIdx); findDaughterParticles(tracksByColl, tracks); + createKFDaughters(tracks); createKFHypernuclei(tracks); createKFCascades(tracks); if (!collHasCandidate) @@ -1035,7 +1038,7 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] fillTree(tracks); } } - PROCESS_SWITCH(hypKfRecoTask, processData, "data analysis", true); + PROCESS_SWITCH(HypKfRecoTask, processData, "data analysis", true); //---------------------------------------------------------------------------------------------------------------- void initCCDB(aod::BCsWithTimestamps::iterator const& bc) { @@ -1078,6 +1081,10 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] { foundDaughters.clear(); foundDaughters.resize(nDaughterParticles); + foundDaughterKfs.clear(); + foundDaughterKfs.resize(nDaughterParticles); + hypNucDaughterKfs.clear(); + hypNucDaughterKfs.resize(nCascades); singleHyperNucCandidates.clear(); singleHyperNucCandidates.resize(nHyperNuclei); cascadeHyperNucCandidates.clear(); @@ -1235,7 +1242,7 @@ struct hypKfRecoTask { // o2-linter: disable=[name/workflow-file][name/struct] WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc)}; } //---------------------------------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------------------------------- diff --git a/PWGLF/TableProducer/Nuspex/hypKfTreeCreator.cxx b/PWGLF/TableProducer/Nuspex/hypKfTreeCreator.cxx index 1eea2943112..bcd2dbe227c 100644 --- a/PWGLF/TableProducer/Nuspex/hypKfTreeCreator.cxx +++ b/PWGLF/TableProducer/Nuspex/hypKfTreeCreator.cxx @@ -320,7 +320,7 @@ DECLARE_SOA_TABLE(HypKfMcCascadeThreeTwoCandidates, "AOD", "HYPKFMCCAND32", HYPK using HypKfMcCascadeThreeTwoCandidate = HypKfMcCascadeThreeTwoCandidates::iterator; } // namespace o2::aod -struct hypKfTreeCreator { // o2-linter: disable=[name/workflow-file][name/struct] +struct HypKfTreeCreator { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Produces outputMcGenTable; @@ -370,7 +370,7 @@ struct hypKfTreeCreator { // o2-linter: disable=[name/workflow-file][name/struct fillTable(candidate, hypDaughter); } } - PROCESS_SWITCH(hypKfTreeCreator, processData, "single tree", false); + PROCESS_SWITCH(HypKfTreeCreator, processData, "single tree", false); //___________________________________________________________________________________________________________________________________________________________ void fillTable(HyperNucleus& cand, HyperNucleus& hypDaughter) { @@ -648,7 +648,7 @@ struct hypKfTreeCreator { // o2-linter: disable=[name/workflow-file][name/struct } hPt[2]->Divide(hPt[1].get(), hPt[0].get()); } - PROCESS_SWITCH(hypKfTreeCreator, processMC, "MC Gen tree", false); + PROCESS_SWITCH(HypKfTreeCreator, processMC, "MC Gen tree", false); //___________________________________________________________________________________________________________________________________________________________ std::vector dcaTracksAll(std::vector& posVec, TString opt = "") @@ -770,5 +770,5 @@ struct hypKfTreeCreator { // o2-linter: disable=[name/workflow-file][name/struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc)}; } diff --git a/PWGLF/Tasks/Resonances/doublephimeson.cxx b/PWGLF/Tasks/Resonances/doublephimeson.cxx index 1e4aa9e0c94..097f37dab57 100644 --- a/PWGLF/Tasks/Resonances/doublephimeson.cxx +++ b/PWGLF/Tasks/Resonances/doublephimeson.cxx @@ -17,11 +17,13 @@ #include #include #include +#include #include #include #include #include #include +#include #include "Framework/AnalysisTask.h" #include "Framework/ASoAHelpers.h" @@ -41,7 +43,6 @@ struct doublephimeson { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Configurable fillRotation{"fillRotation", 1, "Fill rotation"}; - Configurable fillPhiMass{"fillPhiMass", 1, "Fill phi mass"}; Configurable strategyPID1{"strategyPID1", 0, "PID strategy 1"}; Configurable strategyPID2{"strategyPID2", 0, "PID strategy 2"}; Configurable minPhiMass{"minPhiMass", 1.01, "Minimum phi mass"}; @@ -59,9 +60,10 @@ struct doublephimeson { ConfigurableAxis configThnAxisInvMassPhi{"configThnAxisInvMassPhi", {20, 1.01, 1.03}, "#it{M} (GeV/#it{c}^{2})"}; ConfigurableAxis configThnAxisDaugherPt{"configThnAxisDaugherPt", {25, 0.0, 50.}, "#it{p}_{T} (GeV/#it{c})"}; ConfigurableAxis configThnAxisPt{"configThnAxisPt", {40, 0.0, 20.}, "#it{p}_{T} (GeV/#it{c})"}; - // ConfigurableAxis configThnAxisKstar{"configThnAxisKstar", {50, 0.0, 0.5}, "#it{k}^{*} (GeV/#it{c})"}; + ConfigurableAxis configThnAxisKstar{"configThnAxisKstar", {200, 0.0, 2.0}, "#it{k}^{*} (GeV/#it{c})"}; ConfigurableAxis configThnAxisDeltaR{"configThnAxisDeltaR", {200, 0.0, 2.0}, "#it{k}^{*} (GeV/#it{c})"}; - ConfigurableAxis configThnAxisPhiMult{"configThnAxisPhiMult", {10, 0.5, 10.5}, "#Phi Multiplicity"}; + ConfigurableAxis configThnAxisCosTheta{"configThnAxisCosTheta", {5, 0.0, 1.0}, "cos #theta{*}"}; + // ConfigurableAxis configThnAxisPhiMult{"configThnAxisPhiMult", {10, 0.5, 10.5}, "#Phi Multiplicity"}; // Initialize the ananlysis task void init(o2::framework::InitContext&) @@ -76,20 +78,14 @@ struct doublephimeson { const AxisSpec thnAxisInvMassPhi{configThnAxisInvMassPhi, "#it{M} (GeV/#it{c}^{2})"}; const AxisSpec thnAxisDaughterPt{configThnAxisDaugherPt, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec thnAxisPt{configThnAxisPt, "#it{p}_{T} (GeV/#it{c})"}; - // const AxisSpec thnAxisKstar{configThnAxisKstar, "#it{k}^{*} (GeV/#it{c})"}; + const AxisSpec thnAxisKstar{configThnAxisKstar, "#it{k}^{*} (GeV/#it{c})"}; const AxisSpec thnAxisDeltaR{configThnAxisDeltaR, "#Delta R)"}; - const AxisSpec thnAxisPhiMult{configThnAxisPhiMult, "#Phi Multiplicity)"}; + const AxisSpec thnAxisCosTheta{configThnAxisCosTheta, "cos #theta{*}"}; + // const AxisSpec thnAxisPhiMult{configThnAxisPhiMult, "#Phi Multiplicity)"}; histos.add("SEMass", "SEMass", HistType::kTHnSparseF, {thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisDaughterPt, thnAxisDaughterPt}); - if (!fillPhiMass) { - histos.add("SEMassUnlike", "SEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDaughterPt, thnAxisDaughterPt, thnAxisDeltaR, thnAxisPhiMult}); - histos.add("SEMassRot", "SEMassRot", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDaughterPt, thnAxisDaughterPt, thnAxisDeltaR, thnAxisPhiMult}); - histos.add("MEMassUnlike", "MEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDaughterPt, thnAxisDaughterPt, thnAxisDeltaR, thnAxisPhiMult}); - } - if (fillPhiMass) { - histos.add("SEMassUnlike", "SEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisDeltaR, thnAxisPhiMult}); - histos.add("SEMassRot", "SEMassRot", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisDeltaR, thnAxisPhiMult}); - histos.add("MEMassUnlike", "MEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisDeltaR, thnAxisPhiMult}); - } + histos.add("SEMassUnlike", "SEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisKstar, thnAxisCosTheta, thnAxisDeltaR}); + histos.add("SEMassRot", "SEMassRot", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisKstar, thnAxisCosTheta, thnAxisDeltaR}); + histos.add("MEMassUnlike", "MEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisKstar, thnAxisCosTheta, thnAxisDeltaR}); } // get kstar @@ -114,6 +110,25 @@ struct doublephimeson { trackRelK = PartOneCMS - PartTwoCMS; return 0.5 * trackRelK.P(); } + + // get cosTheta + TLorentzVector daughterCMS; + ROOT::Math::XYZVector threeVecDauCM, threeVecMother; + float getCosTheta(const TLorentzVector mother, + const TLorentzVector daughter) + { + threeVecMother = mother.Vect(); + const float beta = mother.Beta(); + const float betax = beta * std::cos(mother.Phi()) * std::sin(mother.Theta()); + const float betay = beta * std::sin(mother.Phi()) * std::sin(mother.Theta()); + const float betaz = beta * std::cos(mother.Theta()); + const ROOT::Math::Boost boostPRF = ROOT::Math::Boost(-betax, -betay, -betaz); + daughterCMS = boostPRF(daughter); + threeVecDauCM = daughterCMS.Vect(); + float cosThetaStar = TMath::Abs(threeVecDauCM.Dot(threeVecMother) / std::sqrt(threeVecMother.Mag2()) / std::sqrt(threeVecDauCM.Mag2())); + return cosThetaStar; + } + bool selectionPID(float nsigmaTPC, float nsigmaTOF, int TOFHit, int PIDStrategy, float ptcand) { if (PIDStrategy == 0) { @@ -202,7 +217,6 @@ struct doublephimeson { if (additionalEvsel && (collision.numPos() < 2 || collision.numNeg() < 2)) { return; } - auto phimult = phitracks.size(); for (auto phitrackd1 : phitracks) { if (phitrackd1.phiMass() < minPhiMass || phitrackd1.phiMass() > maxPhiMass) { continue; @@ -248,14 +262,10 @@ struct doublephimeson { } Phid2.SetXYZM(phitrackd2.phiPx(), phitrackd2.phiPy(), phitrackd2.phiPz(), phitrackd2.phiMass()); exotic = Phid1 + Phid2; - // auto kstar = getkstar(Phid1, Phid2); + auto cosThetaStar = getCosTheta(exotic, Phid1); + auto kstar = getkstar(Phid1, Phid2); auto deltaR = TMath::Sqrt(TMath::Power(Phid1.Phi() - Phid2.Phi(), 2.0) + TMath::Power(Phid1.Eta() - Phid2.Eta(), 2.0)); - if (!fillPhiMass) { - histos.fill(HIST("SEMassUnlike"), exotic.M(), exotic.Pt(), Phid1.Pt(), Phid2.Pt(), deltaR, phimult); - } - if (fillPhiMass) { - histos.fill(HIST("SEMassUnlike"), exotic.M(), exotic.Pt(), Phid1.M(), Phid2.M(), deltaR, phimult); - } + histos.fill(HIST("SEMassUnlike"), exotic.M(), exotic.Pt(), kstar, cosThetaStar, deltaR); histos.fill(HIST("SEMass"), Phid1.M(), Phid2.M(), Phid1.Pt(), Phid2.Pt()); if (fillRotation) { for (int nrotbkg = 0; nrotbkg < 5; nrotbkg++) { @@ -267,14 +277,10 @@ struct doublephimeson { auto rotd1py = Phid1.Px() * std::sin(rotangle) + Phid1.Py() * std::cos(rotangle); Phid1Rot.SetXYZM(rotd1px, rotd1py, Phid1.Pz(), Phid1.M()); exoticRot = Phid1Rot + Phid2; - // auto kstar_rot = getkstar(Phid1Rot, Phid2); + auto cosThetaStar_rot = getCosTheta(exoticRot, Phid1Rot); + auto kstar_rot = getkstar(Phid1Rot, Phid2); auto deltaR_rot = TMath::Sqrt(TMath::Power(Phid1Rot.Phi() - Phid2.Phi(), 2.0) + TMath::Power(Phid1Rot.Eta() - Phid2.Eta(), 2.0)); - if (!fillPhiMass) { - histos.fill(HIST("SEMassRot"), exoticRot.M(), exoticRot.Pt(), Phid1Rot.Pt(), Phid2.Pt(), deltaR_rot, phimult); - } - if (fillPhiMass) { - histos.fill(HIST("SEMassRot"), exoticRot.M(), exoticRot.Pt(), Phid1Rot.M(), Phid2.M(), deltaR_rot, phimult); - } + histos.fill(HIST("SEMassRot"), exoticRot.M(), exoticRot.Pt(), kstar_rot, cosThetaStar_rot, deltaR_rot); } } } @@ -299,7 +305,6 @@ struct doublephimeson { if (additionalEvsel && (collision2.numPos() < 2 || collision2.numNeg() < 2)) { continue; } - auto phimult = tracks1.size(); for (auto& [phitrackd1, phitrackd2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { if (phitrackd1.phiMass() < minPhiMass || phitrackd1.phiMass() > maxPhiMass) { continue; @@ -326,14 +331,10 @@ struct doublephimeson { } Phid2.SetXYZM(phitrackd2.phiPx(), phitrackd2.phiPy(), phitrackd2.phiPz(), phitrackd2.phiMass()); exotic = Phid1 + Phid2; - // auto kstar = getkstar(Phid1, Phid2); + auto cosThetaStar = getCosTheta(exotic, Phid1); + auto kstar = getkstar(Phid1, Phid2); auto deltaR = TMath::Sqrt(TMath::Power(Phid1.Phi() - Phid2.Phi(), 2.0) + TMath::Power(Phid1.Eta() - Phid2.Eta(), 2.0)); - if (!fillPhiMass) { - histos.fill(HIST("MEMassUnlike"), exotic.M(), exotic.Pt(), Phid1.Pt(), Phid2.Pt(), deltaR, phimult); - } - if (fillPhiMass) { - histos.fill(HIST("MEMassUnlike"), exotic.M(), exotic.Pt(), Phid1.M(), Phid2.M(), deltaR, phimult); - } + histos.fill(HIST("MEMassUnlike"), exotic.M(), exotic.Pt(), kstar, cosThetaStar, deltaR); } } } diff --git a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx index 9cc2bd3462e..0c43956bc71 100644 --- a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx +++ b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx @@ -250,17 +250,21 @@ struct NonPromptCascadeTask { return; } mRunNumber = bc.runNumber(); - auto timestamp = bc.timestamp(); - if (o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(cfgGRPpath, timestamp)) { + if (o2::parameters::GRPObject* grpo = ccdb->getForRun(cfgGRPpath, mRunNumber)) { o2::base::Propagator::initFieldFromGRP(grpo); bz = grpo->getNominalL3Field(); - } else if (o2::parameters::GRPMagField* grpmag = ccdb->getForTimeStamp(cfgGRPmagPath, timestamp)) { + } else if (o2::parameters::GRPMagField* grpmag = ccdb->getForRun(cfgGRPmagPath, mRunNumber)) { o2::base::Propagator::initFieldFromGRP(grpmag); bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); LOG(debug) << "bz = " << bz; } else { - LOG(fatal) << "Got nullptr from CCDB for path " << cfgGRPmagPath << " of object GRPMagField and " << cfgGRPpath << " of object GRPObject for timestamp " << timestamp; + LOG(fatal) << "Got nullptr from CCDB for path " << cfgGRPmagPath << " of object GRPMagField and " << cfgGRPpath << " of object GRPObject for run " << mRunNumber; + } + + if (static_cast(cfgMaterialCorrection.value) == o2::base::Propagator::MatCorrType::USEMatCorrLUT) { + auto* lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForRun("GLO/Param/MatLUT", mRunNumber)); + o2::base::Propagator::Instance(true)->setMatLUT(lut); } } @@ -271,11 +275,6 @@ struct NonPromptCascadeTask { ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); - if (static_cast(cfgMaterialCorrection.value) == o2::base::Propagator::MatCorrType::USEMatCorrLUT) { - auto* lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get("GLO/Param/MatLUT")); - o2::base::Propagator::Instance(true)->setMatLUT(lut); - } - std::vector ptBinning = {0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 4.4, 4.8, 5.2, 5.6, 6.0}; AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"};