diff --git a/CODEOWNERS b/CODEOWNERS index 523a5a8cf64..6a84d0d82e6 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -37,8 +37,8 @@ /PWGEM/Dilepton @alibuild @mikesas @rbailhac @dsekihat @ivorobye @feisenhu /PWGEM/PhotonMeson @alibuild @mikesas @rbailhac @m-c-danisch @novitzky @mhemmer-cern @dsekihat /PWGHF @alibuild @vkucera @fcolamar @fgrosa @fcatalan92 @mfaggin @mmazzilli @deepathoms @NicoleBastid @hahassan7 @jpxrk @apalasciano -/PWGLF @alibuild @ercolessi @fmazzasc @chiarapinto @BongHwi @smaff92 @mbombara @ChiaraDeMartin95 @njacazio @skundu692 -/PWGMM @alibuild @aalkin +/PWGLF @alibuild @ercolessi @fmazzasc @chiarapinto @maciacco @BongHwi @smaff92 @ChiaraDeMartin95 @njacazio @skundu692 +/PWGMM @alibuild @aalkin @njacazio @skundu692 /PWGMM/Lumi @alibuild @aalkin /PWGMM/Mult @alibuild @aalkin @aortizve @ddobrigk /PWGMM/UE @alibuild @aalkin @aortizve diff --git a/DPG/Tasks/AOTTrack/PID/TPC/qaPIDTPCMC.cxx b/DPG/Tasks/AOTTrack/PID/TPC/qaPIDTPCMC.cxx index e99f6a9116c..993c2f7b24c 100644 --- a/DPG/Tasks/AOTTrack/PID/TPC/qaPIDTPCMC.cxx +++ b/DPG/Tasks/AOTTrack/PID/TPC/qaPIDTPCMC.cxx @@ -303,6 +303,7 @@ struct pidTpcQaMc { histos.add(hnsigmaMCmat[mcID * Np + massID].data(), Form("True Secondary %s from material", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); if constexpr (mcID == massID) { + histos.add(hsignalMC[mcID].data(), Form("%s", pT[mcID]), HistType::kTH2F, {pAxis, signalAxis}); histos.add(hsignalMCprm[mcID].data(), Form("Primary %s", pT[mcID]), HistType::kTH2F, {pAxis, signalAxis}); histos.add(hsignalMCstr[mcID].data(), Form("Secondary %s from decay", pT[mcID]), HistType::kTH2F, {pAxis, signalAxis}); histos.add(hsignalMCmat[mcID].data(), Form("Secondary %s from material", pT[mcID]), HistType::kTH2F, {pAxis, signalAxis}); diff --git a/PWGCF/Flow/Tasks/FlowGFWPbPb.cxx b/PWGCF/Flow/Tasks/FlowGFWPbPb.cxx index 8ba6dc54df9..a498dcacd67 100644 --- a/PWGCF/Flow/Tasks/FlowGFWPbPb.cxx +++ b/PWGCF/Flow/Tasks/FlowGFWPbPb.cxx @@ -44,6 +44,18 @@ using namespace o2::aod::evsel; #define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; +static constexpr TrackSelectionFlags::flagtype trackSelectionITS = + TrackSelectionFlags::kITSNCls | TrackSelectionFlags::kITSChi2NDF | + TrackSelectionFlags::kITSHits; +static constexpr TrackSelectionFlags::flagtype trackSelectionTPC = + TrackSelectionFlags::kTPCNCls | + TrackSelectionFlags::kTPCCrossedRowsOverNCls | + TrackSelectionFlags::kTPCChi2NDF; +static constexpr TrackSelectionFlags::flagtype trackSelectionDCA = + TrackSelectionFlags::kDCAz | TrackSelectionFlags::kDCAxy; +static constexpr TrackSelectionFlags::flagtype trackSelectionDCAXYonly = + TrackSelectionFlags::kDCAxy; + struct FlowGFWPbPb { O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range") @@ -60,6 +72,8 @@ struct FlowGFWPbPb { O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object") O2_DEFINE_CONFIGURABLE(cfgMagnetField, std::string, "GLO/Config/GRPMagField", "CCDB path to Magnet field object") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 500, "High cut on TPC occupancy") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy") O2_DEFINE_CONFIGURABLE(dcaZ, float, 0.2f, "Custom DCA Z cut (ignored if negative)") O2_DEFINE_CONFIGURABLE(GlobalplusITS, bool, false, "Global and ITS tracks") O2_DEFINE_CONFIGURABLE(Globalonly, bool, false, "Global only tracks") @@ -77,18 +91,8 @@ struct FlowGFWPbPb { ConfigurableAxis axisT0C{"axisT0C", {70, 0, 70000}, "N_{ch} (T0C)"}; ConfigurableAxis axisT0A{"axisT0A", {200, 0, 200000}, "N_{ch} (T0A)"}; ConfigurableAxis axisNchPV{"axisNchPV", {4000, 0, 4000}, "N_{ch} (PV)"}; - - static constexpr TrackSelectionFlags::flagtype trackSelectionITS = - TrackSelectionFlags::kITSNCls | TrackSelectionFlags::kITSChi2NDF | - TrackSelectionFlags::kITSHits; - static constexpr TrackSelectionFlags::flagtype trackSelectionTPC = - TrackSelectionFlags::kTPCNCls | - TrackSelectionFlags::kTPCCrossedRowsOverNCls | - TrackSelectionFlags::kTPCChi2NDF; - static constexpr TrackSelectionFlags::flagtype trackSelectionDCA = - TrackSelectionFlags::kDCAz | TrackSelectionFlags::kDCAxy; - static constexpr TrackSelectionFlags::flagtype trackSelectionDCAXYonly = - TrackSelectionFlags::kDCAxy; + ConfigurableAxis axisDCAz{"axisDCAz", {200, -2, 2}, "DCA_{z} (cm)"}; + ConfigurableAxis axisDCAxy{"axisDCAxy", {200, -1, 1}, "DCA_{xy} (cm)"}; // Corrections TH1D* mEfficiency = nullptr; @@ -148,7 +152,8 @@ struct FlowGFWPbPb { registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(2, "After sel8"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(3, "After additional event cut"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(4, "After correction loads"); - registry.add("hPhi", "", {HistType::kTH1D, {axisPhi}}); + registry.add("hPhi", "#phi distribution", {HistType::kTH1D, {axisPhi}}); + registry.add("hPhiWeighted", "corrected #phi distribution", {HistType::kTH1D, {axisPhi}}); registry.add("hEta", "", {HistType::kTH1D, {axisEta}}); registry.add("hVtxZ", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}}); registry.add("hMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); @@ -174,9 +179,10 @@ struct FlowGFWPbPb { registry.add("multT0C_centT0C_Aft", "after cut;Centrality T0C;mulplicity T0C", {HistType::kTH2D, {axisCentForQA, axisT0C}}); // Track types - registry.add("GlobalplusITS", "Global plus ITS;Centrality FT0C;Nch", kTH1F, {axisCentrality}); - registry.add("Globalonly", "Global only;Centrality FT0C;Nch", kTH1F, {axisCentrality}); - registry.add("ITSonly", "ITS only;Centrality FT0C;Nch", kTH1F, {axisCentrality}); + registry.add("GlobalplusITS", "Global plus ITS;Centrality FT0C;No. of Events", kTH1F, {axisCentrality}); + registry.add("Globalonly", "Global only;Centrality FT0C;No. of Events", kTH1F, {axisCentrality}); + registry.add("ITSonly", "ITS only;Centrality FT0C;No. of Events", kTH1F, {axisCentrality}); + registry.add("Events_per_Centrality_Bin", "Events_per_Centrality_Bin;Centrality FT0C;No. of Events", kTH1F, {axisCentrality}); registry.add("GlobalplusITS_Nch_vs_Cent", "Global plus ITS;Centrality (%); M (|#eta| < 0.8);", {HistType::kTH2D, {axisCentrality, axisNch}}); registry.add("Globalonly_Nch_vs_Cent", "Global only;Centrality (%); M (|#eta| < 0.8);", {HistType::kTH2D, {axisCentrality, axisNch}}); registry.add("ITSonly_Nch_vs_Cent", "ITS only;Centrality (%); M (|#eta| < 0.8);", {HistType::kTH2D, {axisCentrality, axisNch}}); @@ -189,6 +195,8 @@ struct FlowGFWPbPb { registry.add("hChi2prTPCcls", "#chi^{2}/cluster for the TPC track segment", {HistType::kTH1D, {{100, 0., 5.}}}); registry.add("hnTPCClu", "Number of found TPC clusters", {HistType::kTH1D, {{100, 40, 180}}}); registry.add("hnTPCCrossedRow", "Number of crossed TPC Rows", {HistType::kTH1D, {{100, 40, 180}}}); + registry.add("hDCAz", "DCAz after cuts", {HistType::kTH1D, {{100, -3, 3}}}); + registry.add("hDCAxy", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{50, -1, 1}, {50, 0, 10}}}); // additional Output histograms registry.add("c22", ";Centrality (%) ; C_{2}{2} ", {HistType::kTProfile, {axisCentrality}}); @@ -392,25 +400,28 @@ struct FlowGFWPbPb { } auto multNTracksPV = collision.multNTracksPV(); + auto occupancy = collision.trackOccupancyInTimeRange(); if (centrality >= 70. || centrality < 0) - return false; + return 0; if (abs(vtxz) > cfgCutVertex) - return false; + return 0; if (multNTracksPV < fMultPVCutLow->Eval(centrality)) - return false; + return 0; if (multNTracksPV > fMultPVCutHigh->Eval(centrality)) - return false; + return 0; if (multTrk < fMultCutLow->Eval(centrality)) - return false; + return 0; if (multTrk > fMultCutHigh->Eval(centrality)) - return false; + return 0; + if (occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh) + return 0; // V0A T0A 5 sigma cut if (abs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > 5 * fT0AV0ASigma->Eval(collision.multFT0A())) - return false; + return 0; - return true; + return 1; } int getMagneticField(uint64_t timestamp) @@ -449,14 +460,10 @@ struct FlowGFWPbPb { // Apply process filters Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; - Filter trackSelectionProperMixed = ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::ITS) && - ncheckbit(aod::track::trackCutFlag, trackSelectionITS) && - ifnode(ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::TPC), - ncheckbit(aod::track::trackCutFlag, trackSelectionTPC), true) && - ifnode(dcaZ.node() > 0.f, nabs(aod::track::dcaZ) <= dcaZ && ncheckbit(aod::track::trackCutFlag, trackSelectionDCAXYonly), - ncheckbit(aod::track::trackCutFlag, trackSelectionDCA)) && - (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && - (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls); + Filter trackFilter = ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::ITS) && + ncheckbit(aod::track::trackCutFlag, trackSelectionITS) && + ifnode(ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::TPC), ncheckbit(aod::track::trackCutFlag, trackSelectionTPC), true) && + ifnode(dcaZ.node() > 0.f, nabs(aod::track::dcaZ) <= dcaZ && ncheckbit(aod::track::trackCutFlag, trackSelectionDCAXYonly), ncheckbit(aod::track::trackCutFlag, trackSelectionDCA)); using Colls = soa::Filtered>; // collisions filter using aodTracks = soa::Filtered>; // tracks filter @@ -526,8 +533,13 @@ struct FlowGFWPbPb { for (auto& track : tracks) { - // if (track.tpcNClsFound() < cfgCutTPCclu) - // continue; + if (std::abs(track.eta()) >= cfgCutEta) { + continue; + } + if (cfgCutPtMin > track.pt() && track.pt() > cfgCutPtMax) { + continue; + } + if (cfgUseAdditionalTrackCut && !trackSelected(track, Magnetfield)) continue; if (cfgOutputNUAWeights) @@ -540,11 +552,14 @@ struct FlowGFWPbPb { if (WithinPtRef) { registry.fill(HIST("hPhi"), track.phi()); + registry.fill(HIST("hPhiWeighted"), track.phi(), wacc); registry.fill(HIST("hEta"), track.eta()); registry.fill(HIST("hPtRef"), track.pt()); registry.fill(HIST("hChi2prTPCcls"), track.tpcChi2NCl()); registry.fill(HIST("hnTPCClu"), track.tpcNClsFound()); registry.fill(HIST("hnTPCCrossedRow"), track.tpcNClsCrossedRows()); + registry.fill(HIST("hDCAz"), track.dcaZ()); + registry.fill(HIST("hDCAxy"), track.dcaXY(), track.pt()); } globalplusits_nch++; @@ -572,6 +587,7 @@ struct FlowGFWPbPb { } // End of track loop + registry.fill(HIST("Events_per_Centrality_Bin"), cent); registry.fill(HIST("GlobalplusITS_Nch_vs_Cent"), cent, globalplusits_nch); registry.fill(HIST("Globalonly_Nch_vs_Cent"), cent, gloabalonly_nch); registry.fill(HIST("ITSonly_Nch_vs_Cent"), cent, itsonly_nch); diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index b92617fc70e..5a52ad563b5 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -765,6 +765,18 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + if (!nameStr.compare("posTrackKaonRej")) { + cut->AddCut(GetAnalysisCut("posTrack")); + cut->AddCut(GetAnalysisCut("kaonRejNsigma")); + return cut; + } + + if (!nameStr.compare("negTrackKaonRej")) { + cut->AddCut(GetAnalysisCut("negTrack")); + cut->AddCut(GetAnalysisCut("kaonRejNsigma")); + return cut; + } + if (!nameStr.compare("pTLow04")) { cut->AddCut(GetAnalysisCut("pTLow04")); return cut; @@ -4787,6 +4799,11 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("kaonRejNsigma")) { + cut->AddCut(VarManager::kTPCnSigmaKa, -3.0, 3.0, true); + return cut; + } + if (!nameStr.compare("kaonPIDnsigma2")) { cut->AddCut(VarManager::kTPCnSigmaKa, -2.0, 2.0); return cut; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index e0a6085f44e..da04adbf3bf 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -29,6 +29,7 @@ #include #include #include +#include #include #include @@ -145,6 +146,7 @@ class VarManager : public TObject kDecayToKPi, // e.g. D0 -> K+ pi- or cc. kTripleCandidateToKPiPi, // e.g. D+ -> K- pi+ pi+ kTripleCandidateToPKPi, // e.g. Lambda_c -> p K- pi+ + kTripleCandidateToKKPi, // e.g. D_s -> K+ K- pi+ kNMaxCandidateTypes }; @@ -2728,6 +2730,21 @@ void VarManager::FillTriple(T1 const& t1, T2 const& t2, T3 const& t3, float* val values[kPhi] = v123.Phi(); values[kRap] = -v123.Rapidity(); } + + if (pairType == kTripleCandidateToKKPi) { + float m1 = o2::constants::physics::MassKaonCharged; + float m2 = o2::constants::physics::MassPionCharged; + + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m1); + ROOT::Math::PtEtaPhiMVector v3(t3.pt(), t3.eta(), t3.phi(), m2); + ROOT::Math::PtEtaPhiMVector v123 = v1 + v2 + v3; + values[kMass] = v123.M(); + values[kPt] = v123.Pt(); + values[kEta] = v123.Eta(); + values[kPhi] = v123.Phi(); + values[kRap] = -v123.Rapidity(); + } } template diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index b810271792a..ac44ffd4657 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -2204,7 +2204,7 @@ struct AnalysisAsymmetricPairing { std::set> globIdxTriplets; // Based on triplet type, make suitable combinations of the partitions - if (tripletType == VarManager::kTripleCandidateToPKPi) { + if (tripletType == VarManager::kTripleCandidateToPKPi || tripletType == VarManager::kTripleCandidateToKKPi) { for (auto& [a1, a2, a3] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs, groupedLegCAssocs))) { readTriplet(a1, a2, a3, tracks, event, tripletType, histNames); } @@ -2324,6 +2324,13 @@ struct AnalysisAsymmetricPairing { runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, VarManager::kTripleCandidateToKPiPi); } + void processKaonKaonPionSkimmed(MyEventsVtxCovZdcSelected const& events, + soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) + { + runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, VarManager::kTripleCandidateToKKPi); + } + void processProtonKaonPionSkimmed(MyEventsVtxCovZdcSelected const& events, soa::Join const& barrelAssocs, MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) @@ -2338,6 +2345,7 @@ struct AnalysisAsymmetricPairing { PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionSkimmed, "Run kaon pion pairing, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionPionSkimmed, "Run kaon pion pion triplets, with skimmed tracks", false); + PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonKaonPionSkimmed, "Run kaon kaon pion triplets, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processProtonKaonPionSkimmed, "Run proton kaon pion triplets, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processDummy, "Dummy function, enabled only if none of the others are enabled", true); }; diff --git a/PWGHF/D2H/DataModel/ReducedDataModel.h b/PWGHF/D2H/DataModel/ReducedDataModel.h index 07fe57f18bc..02cdc1be0b3 100644 --- a/PWGHF/D2H/DataModel/ReducedDataModel.h +++ b/PWGHF/D2H/DataModel/ReducedDataModel.h @@ -295,6 +295,7 @@ DECLARE_SOA_TABLE(HfMcRecRedDpPis, "AOD", "HFMCRECREDDPPI", //! Table with recon hf_cand_b0_reduced::Prong0Id, hf_cand_b0_reduced::Prong1Id, hf_cand_b0::FlagMcMatchRec, + hf_cand_b0::FlagWrongCollision, hf_cand_b0::DebugMcRec, hf_b0_mc::PtMother); @@ -312,6 +313,7 @@ DECLARE_SOA_TABLE(HfMcCheckDpPis, "AOD", "HFMCCHECKDPPI", //! Table with reconst // Table with same size as HFCANDB0 DECLARE_SOA_TABLE(HfMcRecRedB0s, "AOD", "HFMCRECREDB0", //! Reconstruction-level MC information on B0 candidates for reduced workflow hf_cand_b0::FlagMcMatchRec, + hf_cand_b0::FlagWrongCollision, hf_cand_b0::DebugMcRec, hf_b0_mc::PtMother); @@ -374,6 +376,7 @@ DECLARE_SOA_TABLE(HfMcRecRedD0Pis, "AOD", "HFMCRECREDD0PI", //! Table with recon hf_cand_bplus_reduced::Prong0Id, hf_cand_bplus_reduced::Prong1Id, hf_cand_bplus::FlagMcMatchRec, + hf_cand_bplus::FlagWrongCollision, hf_cand_bplus::DebugMcRec, hf_bplus_mc::PtMother); @@ -388,6 +391,7 @@ DECLARE_SOA_TABLE(HfMcCheckD0Pis, "AOD", "HFMCCHECKD0PI", //! Table with reconst // Table with same size as HFCANDBPLUS DECLARE_SOA_TABLE(HfMcRecRedBps, "AOD", "HFMCRECREDBP", //! Reconstruction-level MC information on B+ candidates for reduced workflow hf_cand_bplus::FlagMcMatchRec, + hf_cand_bplus::FlagWrongCollision, hf_cand_bplus::DebugMcRec, hf_bplus_mc::PtMother); diff --git a/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx index 383ca9c42b1..3e354c266ea 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx @@ -293,7 +293,7 @@ struct HfCandidateCreatorB0ReducedExpressions { if ((rowDPiMcRec.prong0Id() != candB0.prong0Id()) || (rowDPiMcRec.prong1Id() != candB0.prong1Id())) { continue; } - rowB0McRec(rowDPiMcRec.flagMcMatchRec(), rowDPiMcRec.debugMcRec(), rowDPiMcRec.ptMother()); + rowB0McRec(rowDPiMcRec.flagMcMatchRec(), rowDPiMcRec.flagWrongCollision(), rowDPiMcRec.debugMcRec(), rowDPiMcRec.ptMother()); filledMcInfo = true; if constexpr (checkDecayTypeMc) { rowB0McCheck(rowDPiMcRec.pdgCodeBeautyMother(), @@ -306,7 +306,7 @@ struct HfCandidateCreatorB0ReducedExpressions { break; } if (!filledMcInfo) { // protection to get same size tables in case something went wrong: we created a candidate that was not preselected in the D-Pi creator - rowB0McRec(0, -1, -1.f); + rowB0McRec(0, -1, -1, -1.f); if constexpr (checkDecayTypeMc) { rowB0McCheck(-1, -1, -1, -1, -1, -1); } diff --git a/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx index d47ac5f3d4d..34edf46cba4 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx @@ -291,7 +291,7 @@ struct HfCandidateCreatorBplusReducedExpressions { if ((rowD0PiMcRec.prong0Id() != candBplus.prong0Id()) || (rowD0PiMcRec.prong1Id() != candBplus.prong1Id())) { continue; } - rowBplusMcRec(rowD0PiMcRec.flagMcMatchRec(), rowD0PiMcRec.debugMcRec(), rowD0PiMcRec.ptMother()); + rowBplusMcRec(rowD0PiMcRec.flagMcMatchRec(), rowD0PiMcRec.flagWrongCollision(), rowD0PiMcRec.debugMcRec(), rowD0PiMcRec.ptMother()); filledMcInfo = true; if constexpr (checkDecayTypeMc) { rowBplusMcCheck(rowD0PiMcRec.pdgCodeBeautyMother(), @@ -302,7 +302,7 @@ struct HfCandidateCreatorBplusReducedExpressions { break; } if (!filledMcInfo) { // protection to get same size tables in case something went wrong: we created a candidate that was not preselected in the D0-Pi creator - rowBplusMcRec(0, -1, -1.f); + rowBplusMcRec(0, -1, -1, -1.f); if constexpr (checkDecayTypeMc) { rowBplusMcCheck(-1, -1, -1, -1); } diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx index 7ce334dc8b9..ae794eca3c4 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx @@ -57,6 +57,12 @@ enum DecayChannel : uint8_t { BplusToD0barPi, }; +enum WrongCollisionType : uint8_t { + None = 0, + WrongAssociation, + SplitCollision, +}; + /// Creation of CharmHad-Pi pairs for Beauty hadrons struct HfDataCreatorCharmHadPiReduced { // Produces AOD tables to store track information @@ -145,6 +151,8 @@ struct HfDataCreatorCharmHadPiReduced { using CandsD0Filtered = soa::Filtered>; using CandsD0FilteredWithMl = soa::Filtered>; + using CollisionsWMcLabels = soa::Join; + Filter filterSelectDplusCandidates = (aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); Filter filterSelectDzeroCandidates = (aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar); @@ -153,6 +161,7 @@ struct HfDataCreatorCharmHadPiReduced { Preslice candsD0PerCollision = aod::track_association::collisionId; Preslice candsD0PerCollisionWithMl = aod::track_association::collisionId; Preslice trackIndicesPerCollision = aod::track_association::collisionId; + PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; std::shared_ptr hCandidatesD0, hCandidatesDPlus; HistogramRegistry registry{"registry"}; @@ -295,21 +304,62 @@ struct HfDataCreatorCharmHadPiReduced { return true; } + /// Calculates the index of the collision with the maximum number of contributions. + ///\param collisions are the collisions to search through. + ///\return The index of the collision with the maximum number of contributions. + template + int64_t getIndexCollisionMaxNumContrib(const CColl& collisions) + { + unsigned maxNumContrib = 0; + int64_t indexCollisionMaxNumContrib = -1; + for (const auto& collision : collisions) { + if (collision.numContrib() > maxNumContrib) { + maxNumContrib = collision.numContrib(); + indexCollisionMaxNumContrib = collision.globalIndex(); + } + } + return indexCollisionMaxNumContrib; + } + + /// Checks if the B meson is associated with a different collision than the one it was generated in + /// \param particleMother is the mother particle + /// \param collision is the reconstructed collision + /// \param indexCollisionMaxNumContrib is the index of the collision associated with a given MC collision with the largest number of contributors. + /// \param flagWrongCollision is the flag indicating if whether the associated collision is incorrect. + template + void checkWrongCollision(const PParticle& particleMother, + const CColl& collision, + const int64_t& indexCollisionMaxNumContrib, + int8_t& flagWrongCollision) + { + + if (particleMother.mcCollision().globalIndex() != collision.mcCollisionId()) { + flagWrongCollision = WrongCollisionType::WrongAssociation; + } else { + if (collision.globalIndex() != indexCollisionMaxNumContrib) { + flagWrongCollision = WrongCollisionType::SplitCollision; + } + } + } + /// Function for filling MC reco information in the tables /// \param particlesMc is the table with MC particles /// \param vecDaughtersB is the vector with all daughter tracks (bachelor pion in last position) /// \param indexHfCandCharm is the index of the charm-hadron candidate /// \param selectedTracksPion is the map with the indices of selected bachelor pion tracks - template - void fillMcRecoInfo(const PParticles& particlesMc, + template + void fillMcRecoInfo(const CColl& collision, + const PParticles& particlesMc, const std::vector& vecDaughtersB, int& indexHfCandCharm, - std::map selectedTracksPion) + std::map selectedTracksPion, + const int64_t indexCollisionMaxNumContrib) { // we check the MC matching to be stored int8_t sign{0}; int8_t flag{0}; + int8_t flagWrongCollision{WrongCollisionType::None}; int8_t debug{0}; int pdgCodeBeautyMother{-1}; int pdgCodeCharmMother{-1}; @@ -337,6 +387,7 @@ struct HfDataCreatorCharmHadPiReduced { if (indexMother >= 0) { auto particleMother = particlesMc.rawIteratorAt(indexMother); motherPt = particleMother.pt(); + checkWrongCollision(particleMother, collision, indexCollisionMaxNumContrib, flagWrongCollision); } } @@ -419,7 +470,7 @@ struct HfDataCreatorCharmHadPiReduced { } rowHfDPiMcCheckReduced(pdgCodeBeautyMother, pdgCodeCharmMother, pdgCodeProng0, pdgCodeProng1, pdgCodeProng2, pdgCodeProng3); } - rowHfDPiMcRecReduced(indexHfCandCharm, selectedTracksPion[vecDaughtersB.back().globalIndex()], flag, debug, motherPt); + rowHfDPiMcRecReduced(indexHfCandCharm, selectedTracksPion[vecDaughtersB.back().globalIndex()], flag, flagWrongCollision, debug, motherPt); } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // B+ → D0(bar) π+ → (K+ π-) π+ auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1], vecDaughtersB[2]}, Pdg::kBPlus, std::array{+kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); @@ -437,9 +488,10 @@ struct HfDataCreatorCharmHadPiReduced { if (indexMother >= 0) { auto particleMother = particlesMc.rawIteratorAt(indexMother); motherPt = particleMother.pt(); + checkWrongCollision(particleMother, collision, indexCollisionMaxNumContrib, flagWrongCollision); } } - rowHfD0PiMcRecReduced(indexHfCandCharm, selectedTracksPion[vecDaughtersB.back().globalIndex()], flag, debug, motherPt); + rowHfD0PiMcRecReduced(indexHfCandCharm, selectedTracksPion[vecDaughtersB.back().globalIndex()], flag, flagWrongCollision, debug, motherPt); // additional checks for correlated backgrounds if (checkDecayTypeMc) { @@ -482,6 +534,7 @@ struct HfDataCreatorCharmHadPiReduced { aod::TrackAssoc const& trackIndices, TTracks const&, PParticles const& particlesMc, + uint64_t const& indexCollisionMaxNumContrib, aod::BCsWithTimestamps const&) { // helpers for ReducedTables filling @@ -683,7 +736,7 @@ struct HfDataCreatorCharmHadPiReduced { beautyHadDauTracks.push_back(track); } beautyHadDauTracks.push_back(trackPion); - fillMcRecoInfo(particlesMc, beautyHadDauTracks, indexHfCandCharm, selectedTracksPion); + fillMcRecoInfo(collision, particlesMc, beautyHadDauTracks, indexHfCandCharm, selectedTracksPion, indexCollisionMaxNumContrib); } fillHfCandCharm = true; } // pion loop @@ -789,7 +842,7 @@ struct HfDataCreatorCharmHadPiReduced { ptProngs[1], yProngs[1], etaProngs[1]); } else if constexpr (decayChannel == DecayChannel::BplusToD0barPi) { // B+ → D0bar π+ - if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kBPlus, std::array{static_cast(Pdg::kD0), +kPiPlus}, true)) { + if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kBPlus, std::array{-static_cast(Pdg::kD0), +kPiPlus}, true)) { // Match D0bar -> π- K+ auto candD0MC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); // Printf("Checking D0bar -> π- K+"); @@ -850,7 +903,7 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsDplusPerCollision, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, -1, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); @@ -880,7 +933,7 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsDplusPerCollisionWithMl, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, -1, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); @@ -910,7 +963,7 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsD0PerCollision, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, -1, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); @@ -940,7 +993,7 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsD0PerCollisionWithMl, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, -1, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); @@ -950,12 +1003,13 @@ struct HfDataCreatorCharmHadPiReduced { //////////////////////////////////////////////////////////////////////////////////////////////////// // PROCESS FUNCTIONS FOR MC - void processDplusPiMc(soa::Join const& collisions, + void processDplusPiMc(CollisionsWMcLabels const& collisions, CandsDplusFiltered const& candsC, aod::TrackAssoc const& trackIndices, TracksPidWithSelAndMc const& tracks, aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) + aod::BCsWithTimestamps const& bcs, + McCollisions const&) { // store configurables needed for B0 workflow if (!isHfCandBhadConfigFilled) { @@ -974,7 +1028,9 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsDplusPerCollision, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + auto collsSameMcCollision = collisions.sliceBy(colPerMcCollision, collision.mcCollisionId()); + int64_t indexCollisionMaxNumContrib = getIndexCollisionMaxNumContrib(collsSameMcCollision); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, indexCollisionMaxNumContrib, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); @@ -982,12 +1038,13 @@ struct HfDataCreatorCharmHadPiReduced { } PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processDplusPiMc, "Process DplusPi with MC info and without ML info", false); - void processDplusPiMcWithMl(soa::Join const& collisions, + void processDplusPiMcWithMl(CollisionsWMcLabels const& collisions, CandsDplusFilteredWithMl const& candsC, aod::TrackAssoc const& trackIndices, TracksPidWithSelAndMc const& tracks, aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) + aod::BCsWithTimestamps const& bcs, + McCollisions const&) { // store configurables needed for B0 workflow if (!isHfCandBhadConfigFilled) { @@ -1006,7 +1063,9 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsDplusPerCollisionWithMl, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + auto collsSameMcCollision = collisions.sliceBy(colPerMcCollision, collision.mcCollisionId()); + int64_t indexCollisionMaxNumContrib = getIndexCollisionMaxNumContrib(collsSameMcCollision); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, indexCollisionMaxNumContrib, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); @@ -1014,12 +1073,13 @@ struct HfDataCreatorCharmHadPiReduced { } PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processDplusPiMcWithMl, "Process DplusPi with MC info and with ML info", false); - void processD0PiMc(soa::Join const& collisions, + void processD0PiMc(CollisionsWMcLabels const& collisions, CandsD0Filtered const& candsC, aod::TrackAssoc const& trackIndices, TracksPidWithSelAndMc const& tracks, aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) + aod::BCsWithTimestamps const& bcs, + McCollisions const&) { // store configurables needed for B+ workflow if (!isHfCandBhadConfigFilled) { @@ -1038,7 +1098,9 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsD0PerCollision, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + auto collsSameMcCollision = collisions.sliceBy(colPerMcCollision, collision.mcCollisionId()); + int64_t indexCollisionMaxNumContrib = getIndexCollisionMaxNumContrib(collsSameMcCollision); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, indexCollisionMaxNumContrib, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); @@ -1046,12 +1108,13 @@ struct HfDataCreatorCharmHadPiReduced { } PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processD0PiMc, "Process D0Pi with MC info and without ML info", false); - void processD0PiMcWithMl(soa::Join const& collisions, + void processD0PiMcWithMl(CollisionsWMcLabels const& collisions, CandsD0FilteredWithMl const& candsC, aod::TrackAssoc const& trackIndices, TracksPidWithSelAndMc const& tracks, aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) + aod::BCsWithTimestamps const& bcs, + McCollisions const&) { // store configurables needed for B+ workflow if (!isHfCandBhadConfigFilled) { @@ -1070,7 +1133,9 @@ struct HfDataCreatorCharmHadPiReduced { auto thisCollId = collision.globalIndex(); auto candsCThisColl = candsC.sliceBy(candsD0PerCollisionWithMl, thisCollId); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + auto collsSameMcCollision = collisions.sliceBy(colPerMcCollision, collision.mcCollisionId()); + int64_t indexCollisionMaxNumContrib = getIndexCollisionMaxNumContrib(collsSameMcCollision); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, indexCollisionMaxNumContrib, bcs); } // handle normalization by the right number of collisions hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); diff --git a/PWGHF/D2H/Tasks/taskB0Reduced.cxx b/PWGHF/D2H/Tasks/taskB0Reduced.cxx index 5e38055f589..991db21db2b 100644 --- a/PWGHF/D2H/Tasks/taskB0Reduced.cxx +++ b/PWGHF/D2H/Tasks/taskB0Reduced.cxx @@ -58,6 +58,7 @@ DECLARE_SOA_COLUMN(Cpa, cpa, float); //! DECLARE_SOA_COLUMN(CpaXY, cpaXY, float); //! Cosine pointing angle of candidate in transverse plane DECLARE_SOA_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, float); //! Maximum normalized difference between measured and expected impact parameter of candidate prongs DECLARE_SOA_COLUMN(MlScoreSig, mlScoreSig, float); //! ML score for signal class +DECLARE_SOA_COLUMN(FlagWrongCollision, flagWrongCollision, int8_t); //! Flag for association with wrong collision } // namespace hf_cand_b0_lite DECLARE_SOA_TABLE(HfRedCandB0Lites, "AOD", "HFREDCANDB0LITE", //! Table with some B0 properties @@ -89,10 +90,12 @@ DECLARE_SOA_TABLE(HfRedCandB0Lites, "AOD", "HFREDCANDB0LITE", //! Table with som hf_cand_b0_lite::Y, hf_cand_3prong::FlagMcMatchRec, hf_cand_3prong::OriginMcRec, + hf_cand_b0_lite::FlagWrongCollision, hf_cand_b0_lite::PtGen); DECLARE_SOA_TABLE(HfRedB0McCheck, "AOD", "HFREDB0MCCHECK", //! Table with MC decay type check hf_cand_3prong::FlagMcMatchRec, + hf_cand_b0_lite::FlagWrongCollision, hf_cand_b0_lite::MProng0, hf_cand_b0_lite::PtProng0, hf_cand_b0_lite::M, @@ -346,9 +349,11 @@ struct HfTaskB0Reduced { auto decLenXyD = RecoDecay::distanceXY(posPv, posSvD); int8_t flagMcMatchRec = 0; + int8_t flagWrongCollision = 0; bool isSignal = false; if constexpr (doMc) { flagMcMatchRec = candidate.flagMcMatchRec(); + flagWrongCollision = candidate.flagWrongCollision(); isSignal = TESTBIT(std::abs(flagMcMatchRec), hf_cand_b0::DecayTypeMc::B0ToDplusPiToPiKPiPi); } @@ -524,6 +529,7 @@ struct HfTaskB0Reduced { hfHelper.yB0(candidate), flagMcMatchRec, isSignal, + flagWrongCollision, ptMother); } if constexpr (withDecayTypeCheck) { @@ -533,6 +539,7 @@ struct HfTaskB0Reduced { } hfRedB0McCheck( flagMcMatchRec, + flagWrongCollision, invMassD, ptD, invMassB0, diff --git a/PWGHF/D2H/Tasks/taskBplusReduced.cxx b/PWGHF/D2H/Tasks/taskBplusReduced.cxx index 5b56d523944..282194f57b5 100644 --- a/PWGHF/D2H/Tasks/taskBplusReduced.cxx +++ b/PWGHF/D2H/Tasks/taskBplusReduced.cxx @@ -57,6 +57,7 @@ DECLARE_SOA_COLUMN(Cpa, cpa, float); //! DECLARE_SOA_COLUMN(CpaXY, cpaXY, float); //! Cosine pointing angle of candidate in transverse plane DECLARE_SOA_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, float); //! Maximum normalized difference between measured and expected impact parameter of candidate prongs DECLARE_SOA_COLUMN(MlScoreSig, mlScoreSig, float); //! ML score for signal class +DECLARE_SOA_COLUMN(FlagWrongCollision, flagWrongCollision, int8_t); //! Flag for association with wrong collision } // namespace hf_cand_bplus_lite DECLARE_SOA_TABLE(HfRedCandBpLites, "AOD", "HFREDCANDBPLITE", //! Table with some B+ properties @@ -88,10 +89,12 @@ DECLARE_SOA_TABLE(HfRedCandBpLites, "AOD", "HFREDCANDBPLITE", //! Table with som hf_cand_bplus_lite::Y, hf_cand_2prong::FlagMcMatchRec, hf_cand_2prong::OriginMcRec, + hf_cand_bplus_lite::FlagWrongCollision, hf_cand_bplus_lite::PtGen); DECLARE_SOA_TABLE(HfRedBpMcCheck, "AOD", "HFREDBPMCCHECK", //! Table with MC decay type check hf_cand_3prong::FlagMcMatchRec, + hf_cand_bplus_lite::FlagWrongCollision, hf_cand_bplus_lite::MProng0, hf_cand_bplus_lite::PtProng0, hf_cand_bplus_lite::M, @@ -186,10 +189,10 @@ struct HfTaskBplusReduced { registry.add("hRapidity", bPlusCandTitle + "candidate #it{y};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); registry.add("hd0d0", bPlusCandTitle + "candidate product of DCAxy to prim. vertex (cm^{2});" + stringPt, {HistType::kTH2F, {axisImpParProd, axisPtB}}); registry.add("hInvMassD0", bPlusCandTitle + "prong0, D0 inv. mass (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMassD0, axisPtD0}}); - registry.add("hDecLengthD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate decay length (cm);entries", {HistType::kTH2F, {axisPtD0, axisDecLength}}); - registry.add("hDecLengthXyD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});decay length XY (cm);entries", {HistType::kTH2F, {axisPtD0, axisDecLength}}); - registry.add("hCpaD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate cos(#vartheta_{P});entries", {HistType::kTH2F, {axisPtD0, axisCpa}}); - registry.add("hCpaXyD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate cos(#vartheta_{P}^{XY});entries", {HistType::kTH2F, {axisPtD0, axisCpa}}); + registry.add("hDecLengthD0", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate decay length (cm);entries", {HistType::kTH2F, {axisPtD0, axisDecLength}}); + registry.add("hDecLengthXyD0", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});decay length XY (cm);entries", {HistType::kTH2F, {axisPtD0, axisDecLength}}); + registry.add("hCpaD0", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate cos(#vartheta_{P});entries", {HistType::kTH2F, {axisPtD0, axisCpa}}); + registry.add("hCpaXyD0", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate cos(#vartheta_{P}^{XY});entries", {HistType::kTH2F, {axisPtD0, axisCpa}}); // ML scores of D0 daughter if (doprocessDataWithDmesMl) { @@ -340,12 +343,13 @@ struct HfTaskBplusReduced { /// \param candidatesD is the table with D0 candidates template void fillCand(Cand const& candidate, - aod::HfRed2Prongs const& /*candidatesD*/) + aod::HfRed2Prongs const& /*candidatesD*/, + TracksPion const&) { auto ptCandBplus = candidate.pt(); auto invMassBplus = hfHelper.invMassBplusToD0Pi(candidate); auto candD0 = candidate.template prong0_as(); - auto candPi = candidate.template prong1_as(); + auto candPi = candidate.template prong1_as(); auto ptD0 = candidate.ptProng0(); auto invMassD0 = (candPi.signed1Pt() < 0) ? candD0.invMassD0() : candD0.invMassD0Bar(); std::array posPv{candidate.posX(), candidate.posY(), candidate.posZ()}; @@ -357,9 +361,11 @@ struct HfTaskBplusReduced { auto decLenXyD0 = RecoDecay::distanceXY(posPv, posSvD); int8_t flagMcMatchRec = 0; + int8_t flagWrongCollision = 0; bool isSignal = false; if constexpr (doMc) { flagMcMatchRec = candidate.flagMcMatchRec(); + flagWrongCollision = candidate.flagWrongCollision(); isSignal = TESTBIT(std::abs(flagMcMatchRec), hf_cand_bplus::DecayTypeMc::BplusToD0PiToKPiPi); } @@ -537,6 +543,7 @@ struct HfTaskBplusReduced { hfHelper.yBplus(candidate), flagMcMatchRec, isSignal, + flagWrongCollision, ptMother); } if constexpr (withDecayTypeCheck) { @@ -546,6 +553,7 @@ struct HfTaskBplusReduced { } hfRedBpMcCheck( flagMcMatchRec, + flagWrongCollision, invMassD0, ptD0, invMassBplus, @@ -600,39 +608,39 @@ struct HfTaskBplusReduced { // Process functions void processData(soa::Filtered> const& candidates, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // candidate loop } // processData PROCESS_SWITCH(HfTaskBplusReduced, processData, "Process data without ML scores for D0 daughter", true); void processDataWithDmesMl(soa::Filtered> const& candidates, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // candidate loop } // processDataWithDmesMl PROCESS_SWITCH(HfTaskBplusReduced, processDataWithDmesMl, "Process data with ML scores for D0 daughter", false); void processDataWithBplusMl(soa::Filtered> const& candidates, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // candidate loop } // processDataWithBplusMl PROCESS_SWITCH(HfTaskBplusReduced, processDataWithBplusMl, "Process data with(out) ML scores for B+ (D0 daughter)", false); @@ -640,14 +648,14 @@ struct HfTaskBplusReduced { void processMc(soa::Join const& candidates, aod::HfMcGenRedBps const& mcParticles, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { // MC rec for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // rec // MC gen. level @@ -660,14 +668,14 @@ struct HfTaskBplusReduced { void processMcWithDecayTypeCheck(soa::Filtered> const& candidates, aod::HfMcGenRedBps const& mcParticles, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { // MC rec for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // rec // MC gen. level @@ -680,14 +688,14 @@ struct HfTaskBplusReduced { void processMcWithDmesMl(soa::Join const& candidates, aod::HfMcGenRedBps const& mcParticles, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { // MC rec for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // rec // MC gen. level @@ -700,14 +708,14 @@ struct HfTaskBplusReduced { void processMcWithBplusMl(soa::Filtered> const& candidates, aod::HfMcGenRedBps const& mcParticles, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { // MC rec for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // rec // MC gen. level @@ -720,14 +728,14 @@ struct HfTaskBplusReduced { void processMcWithBplusMlAndDecayTypeCheck(soa::Filtered> const& candidates, aod::HfMcGenRedBps const& mcParticles, aod::HfRed2Prongs const& candidatesD, - TracksPion const&) + TracksPion const& pionTracks) { // MC rec for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } - fillCand(candidate, candidatesD); + fillCand(candidate, candidatesD, pionTracks); } // rec // MC gen. level diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 2a886cc277a..13bdd7b378e 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -757,6 +757,7 @@ namespace hf_cand_bplus DECLARE_SOA_INDEX_COLUMN_FULL(Prong0, prong0, int, HfCand2Prong, "_0"); // D0 index // MC matching result: DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(FlagWrongCollision, flagWrongCollision, int8_t); // reconstruction level DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); // generator level DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); // particle origin, reconstruction level DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); // particle origin, generator level @@ -1659,11 +1660,12 @@ namespace hf_cand_b0 { DECLARE_SOA_INDEX_COLUMN_FULL(Prong0, prong0, int, HfCand3Prong, "_0"); // D index // MC matching result: -DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); // reconstruction level -DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); // generator level -DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); // particle origin, reconstruction level -DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); // particle origin, generator level -DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t); // debug flag for mis-association reconstruction level +DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(FlagWrongCollision, flagWrongCollision, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); // generator level +DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); // particle origin, reconstruction level +DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); // particle origin, generator level +DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t); // debug flag for mis-association reconstruction level // mapping of decay types enum DecayType { B0ToDPi }; diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index 96981187d1f..80810f73b11 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -28,6 +28,11 @@ o2physics_add_dpl_workflow(pid-creator PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(mc-pid-tof + SOURCES mcPidTof.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::TOFWorkflowUtils + COMPONENT_NAME Analysis) + # Candidate creators o2physics_add_dpl_workflow(candidate-creator-2prong diff --git a/PWGHF/TableProducer/mcPidTof.cxx b/PWGHF/TableProducer/mcPidTof.cxx new file mode 100644 index 00000000000..2248e8f0bdc --- /dev/null +++ b/PWGHF/TableProducer/mcPidTof.cxx @@ -0,0 +1,440 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file mcPidTof.cxx +/// \author Fabrizio Grosa fabrizio.grosa@cern.ch +/// \brief Task to produce PID tables for TOF split for pi, K, p, copied from https://github.com/AliceO2Group/O2Physics/blob/master/Common/TableProducer/PID/pidTOFFull.cxx +/// In addition, it applies postcalibrations for MC. +/// + +// O2 includes +#include +#include "TOFBase/EventTimeMaker.h" +#include "Framework/AnalysisTask.h" +#include "ReconstructionDataFormats/Track.h" + +// O2Physics includes +#include "TableHelper.h" +#include "Common/TableProducer/PID/pidTOFBase.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::pid; +using namespace o2::framework::expressions; +using namespace o2::track; + +void customize(std::vector& workflowOptions) +{ + std::vector options{{"add-qa", VariantType::Int, 0, {"Legacy. No effect."}}}; + std::swap(workflowOptions, options); +} + +#include "Framework/runDataProcessing.h" + +/// Task to produce the response table +struct mcPidTof { + // Tables to produce + Produces tablePIDPi; + Produces tablePIDKa; + Produces tablePIDPr; + // Detector response parameters + o2::pid::tof::TOFResoParamsV2 mRespParamsV2; + Service ccdb; + Configurable inheritFromBaseTask{"inheritFromBaseTask", true, "Flag to iherit all common configurables from the TOF base task"}; + // CCDB configuration (inherited from TOF base task) + Configurable url{"ccdb-url", "", "url of the ccdb repository"}; + Configurable timestamp{"ccdb-timestamp", -1, "timestamp of the object"}; + // TOF Calib configuration (inherited from TOF base task) + Configurable paramFileName{"paramFileName", "", "Path to the parametrization object. If empty the parametrization is not taken from file"}; + Configurable parametrizationPath{"parametrizationPath", "", "Path of the TOF parametrization on the CCDB or in the file, if the paramFileName is not empty"}; + Configurable passName{"passName", "", "Name of the pass inside of the CCDB parameter collection. If empty, the automatically deceted from metadata (to be implemented!!!)"}; + Configurable timeShiftCCDBPath{"timeShiftCCDBPath", "", "Path of the TOF time shift vs eta. If empty none is taken"}; + Configurable loadResponseFromCCDB{"loadResponseFromCCDB", false, "Flag to load the response from the CCDB"}; + Configurable enableTimeDependentResponse{"enableTimeDependentResponse", false, "Flag to use the collision timestamp to fetch the PID Response"}; + Configurable fatalOnPassNotAvailable{"fatalOnPassNotAvailable", true, "Flag to throw a fatal if the pass is not available in the retrieved CCDB object"}; + // postcalibrations to overcome MC FT0 timing issue + std::map gMcPostCalibMean{}; + std::map gMcPostCalibSigma{}; + int currentRun{0}; + struct : ConfigurableGroup { + std::string prefix = "mcRecalib"; + Configurable enable{"enable", false, "enable MC recalibration for Pi/Ka/Pr"}; + Configurable ccdbPath{"ccdbPath", "Users/f/fgrosa/RecalibMcPidTOF/", "path for MC recalibration objects in CCDB"}; + Configurable passName{"passName", "apass6", "reco pass of MC anchoring"}; + } mcRecalib; + + // Running variables + std::vector mEnabledParticles; // Vector of enabled PID hypotheses to loop on when making tables + int mLastCollisionId = -1; // Last collision ID analysed + void init(o2::framework::InitContext& initContext) + { + if (inheritFromBaseTask.value) { // Inheriting from base task + if (!getTaskOptionValue(initContext, "tof-signal", "ccdb-url", url.value, true)) { + LOG(fatal) << "Could not get ccdb-url from tof-signal task"; + } + if (!getTaskOptionValue(initContext, "tof-signal", "ccdb-timestamp", timestamp.value, true)) { + LOG(fatal) << "Could not get ccdb-timestamp from tof-signal task"; + } + if (!getTaskOptionValue(initContext, "tof-event-time", "paramFileName", paramFileName.value, true)) { + LOG(fatal) << "Could not get paramFileName from tof-event-time task"; + } + if (!getTaskOptionValue(initContext, "tof-event-time", "parametrizationPath", parametrizationPath.value, true)) { + LOG(fatal) << "Could not get parametrizationPath from tof-event-time task"; + } + if (!getTaskOptionValue(initContext, "tof-event-time", "passName", passName.value, true)) { + LOG(fatal) << "Could not get passName from tof-event-time task"; + } + if (!getTaskOptionValue(initContext, "tof-signal", "timeShiftCCDBPath", timeShiftCCDBPath.value, true)) { + LOG(fatal) << "Could not get timeShiftCCDBPath from tof-signal task"; + } + if (!getTaskOptionValue(initContext, "tof-event-time", "loadResponseFromCCDB", loadResponseFromCCDB.value, true)) { + LOG(fatal) << "Could not get loadResponseFromCCDB from tof-event-time task"; + } + if (!getTaskOptionValue(initContext, "tof-event-time", "enableTimeDependentResponse", enableTimeDependentResponse.value, true)) { + LOG(fatal) << "Could not get enableTimeDependentResponse from tof-event-time task"; + } + if (!getTaskOptionValue(initContext, "tof-event-time", "fatalOnPassNotAvailable", fatalOnPassNotAvailable.value, true)) { + LOG(fatal) << "Could not get fatalOnPassNotAvailable from tof-event-time task"; + } + } + if (doprocessWSlice == true && doprocessWoSlice == true) { + LOGF(fatal, "Cannot enable processWoSlice and processWSlice at the same time. Please choose one."); + } + if (doprocessWSlice == false && doprocessWoSlice == false) { + LOGF(fatal, "Cannot run without any of processWoSlice and processWSlice enabled. Please choose one."); + } + + // Checking the tables are requested in the workflow and enabling them (only pi, K, p) + std::array supportedSpecies = {2, 3, 4}; + for (auto iSpecie{0u}; iSpecie < supportedSpecies.size(); ++iSpecie) { + int flag = -1; + enableFlagIfTableRequired(initContext, "mcPidTof" + particleNames[supportedSpecies[iSpecie]], flag); + } + // Printing enabled tables + LOG(info) << "++ Enabled tables:"; + for (const int& pidId : mEnabledParticles) { + LOG(info) << "++ mcPidTof" << particleNames[pidId] << " is enabled"; + } + + // Getting the parametrization parameters + ccdb->setURL(url.value); + ccdb->setTimestamp(timestamp.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + // Not later than now objects + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + // + + // TODO: implement the automatic pass name detection from metadata + if (passName.value == "") { + passName.value = "unanchored"; // temporary default + LOG(warning) << "Passed autodetect mode for pass, not implemented yet, waiting for metadata. Taking '" << passName.value << "'"; + } + LOG(info) << "Using parameter collection, starting from pass '" << passName.value << "'"; + + const std::string fname = paramFileName.value; + if (!fname.empty()) { // Loading the parametrization from file + LOG(info) << "Loading exp. sigma parametrization from file " << fname << ", using param: " << parametrizationPath.value; + if (1) { + o2::tof::ParameterCollection paramCollection; + paramCollection.loadParamFromFile(fname, parametrizationPath.value); + LOG(info) << "+++ Loaded parameter collection from file +++"; + if (!paramCollection.retrieveParameters(mRespParamsV2, passName.value)) { + if (fatalOnPassNotAvailable) { + LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } else { + LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } + } else { + mRespParamsV2.setShiftParameters(paramCollection.getPars(passName.value)); + mRespParamsV2.printShiftParameters(); + } + } else { + mRespParamsV2.loadParamFromFile(fname.data(), parametrizationPath.value); + } + } else if (loadResponseFromCCDB) { // Loading it from CCDB + LOG(info) << "Loading exp. sigma parametrization from CCDB, using path: " << parametrizationPath.value << " for timestamp " << timestamp.value; + o2::tof::ParameterCollection* paramCollection = ccdb->getForTimeStamp(parametrizationPath.value, timestamp.value); + paramCollection->print(); + if (!paramCollection->retrieveParameters(mRespParamsV2, passName.value)) { // Attempt at loading the parameters with the pass defined + if (fatalOnPassNotAvailable) { + LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } else { + LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } + } else { // Pass is available, load non standard parameters + mRespParamsV2.setShiftParameters(paramCollection->getPars(passName.value)); + mRespParamsV2.printShiftParameters(); + } + } + mRespParamsV2.print(); + if (timeShiftCCDBPath.value != "") { + if (timeShiftCCDBPath.value.find(".root") != std::string::npos) { + mRespParamsV2.setTimeShiftParameters(timeShiftCCDBPath.value, "gmean_Pos", true); + mRespParamsV2.setTimeShiftParameters(timeShiftCCDBPath.value, "gmean_Neg", false); + } else { + mRespParamsV2.setTimeShiftParameters(ccdb->getForTimeStamp(Form("%s/pos", timeShiftCCDBPath.value.c_str()), timestamp.value), true); + mRespParamsV2.setTimeShiftParameters(ccdb->getForTimeStamp(Form("%s/neg", timeShiftCCDBPath.value.c_str()), timestamp.value), false); + } + } + } + + // Reserves an empty table for the given particle ID with size of the given track table + void reserveTable(const int id, const int64_t& size) + { + switch (id) { + case 2: + tablePIDPi.reserve(size); + break; + case 3: + tablePIDKa.reserve(size); + break; + case 4: + tablePIDPr.reserve(size); + break; + default: + LOG(fatal) << "Wrong particle ID in reserveTable()"; + break; + } + } + + // Makes the table empty for the given particle ID, filling it with dummy values + void makeTableEmpty(const int id) + { + switch (id) { + case 2: + tablePIDPi(-999.f, -999.f); + break; + case 3: + tablePIDKa(-999.f, -999.f); + break; + case 4: + tablePIDPr(-999.f, -999.f); + break; + default: + LOG(fatal) << "Wrong particle ID in makeTableEmpty()"; + break; + } + } + + /// Retrieve MC postcalibration objects from CCDB + /// \param timestamp timestamp + void retrieveMcPostCalibFromCcdb(int timestamp) + { + std::map metadata; + metadata["RecoPassName"] = mcRecalib.passName; + auto calibList = ccdb->getSpecific(mcRecalib.ccdbPath, timestamp, metadata); + for (auto const& pidId : mEnabledParticles) { // Loop on enabled particle hypotheses + gMcPostCalibMean[pidId] = reinterpret_cast(calibList->FindObject(Form("Mean%s", particleNames[pidId].data()))); + gMcPostCalibSigma[pidId] = reinterpret_cast(calibList->FindObject(Form("Sigma%s", particleNames[pidId].data()))); + } + } + + /// Apply MC postcalibrations + /// \param pidId particle id + /// \param pt track pT + template + T applyMcRecalib(int pidId, T trackPt, T nSigma) + { + float shift{0.f}, scaleWidth{0.f}; + int nPoints = gMcPostCalibMean[pidId]->GetN(); + double ptMin = gMcPostCalibMean[pidId]->GetX()[0]; + double ptMax = gMcPostCalibMean[pidId]->GetX()[nPoints - 1]; + if (trackPt < ptMin) { + shift = gMcPostCalibMean[pidId]->Eval(ptMin); + scaleWidth = gMcPostCalibSigma[pidId]->Eval(ptMin); + } else if (trackPt > ptMax) { + shift = gMcPostCalibMean[pidId]->Eval(ptMax); + scaleWidth = gMcPostCalibSigma[pidId]->Eval(ptMax); + } else { + shift = gMcPostCalibMean[pidId]->Eval(trackPt); + scaleWidth = gMcPostCalibSigma[pidId]->Eval(trackPt); + } + + T nSigmaCorr = (nSigma - shift) / scaleWidth; + return nSigmaCorr; + } + + using Trks = soa::Join; + // Define slice per collision + Preslice perCollision = aod::track::collisionId; + template + using ResponseImplementation = o2::pid::tof::ExpTimes; + void processWSlice(Trks const& tracks, aod::Collisions const&, aod::BCsWithTimestamps const&) + { + constexpr auto responsePi = ResponseImplementation(); + constexpr auto responseKa = ResponseImplementation(); + constexpr auto responsePr = ResponseImplementation(); + + for (auto const& pidId : mEnabledParticles) { + reserveTable(pidId, tracks.size()); + } + + int lastCollisionId = -1; // Last collision ID analysed + float resolution = 1.f; // Last resolution assigned + for (auto const& track : tracks) { // Loop on all tracks + if (!track.has_collision()) { // Track was not assigned, cannot compute NSigma (no event time) -> filling with empty table + for (auto const& pidId : mEnabledParticles) { + makeTableEmpty(pidId); + } + continue; + } + + if (track.collisionId() == lastCollisionId) { // Tracks from last collision already processed + continue; + } + + // Fill new table for the tracks in a collision + lastCollisionId = track.collisionId(); // Cache last collision ID + timestamp.value = track.collision().bc_as().timestamp(); + if (enableTimeDependentResponse) { + LOG(debug) << "Updating parametrization from path '" << parametrizationPath.value << "' and timestamp " << timestamp.value; + if (!ccdb->getForTimeStamp(parametrizationPath.value, timestamp.value)->retrieveParameters(mRespParamsV2, passName.value)) { + if (fatalOnPassNotAvailable) { + LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } else { + LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } + } + } + + // in case of MC recalibrations, check if the objects from CCDB has to be updated + if (mcRecalib.enable) { + int runNumber = track.collision().bc_as().runNumber(); + if (runNumber != currentRun) { + // update postcalibration files + retrieveMcPostCalibFromCcdb(timestamp); + } + currentRun = runNumber; + } + + const auto& tracksInCollision = tracks.sliceBy(perCollision, lastCollisionId); + for (auto const& trkInColl : tracksInCollision) { // Loop on tracks + for (auto const& pidId : mEnabledParticles) { // Loop on enabled particle hypotheses + float nSigma{-999.f}; + switch (pidId) { + case 2: + resolution = responsePi.GetExpectedSigma(mRespParamsV2, trkInColl); + nSigma = responsePi.GetSeparation(mRespParamsV2, trkInColl, resolution); + if (mcRecalib.enable) { + nSigma = applyMcRecalib(pidId, trkInColl.pt(), nSigma); + } + tablePIDPi(resolution, nSigma); + break; + case 3: + resolution = responseKa.GetExpectedSigma(mRespParamsV2, trkInColl); + nSigma = responseKa.GetSeparation(mRespParamsV2, trkInColl, resolution); + if (mcRecalib.enable) { + nSigma = applyMcRecalib(pidId, trkInColl.pt(), nSigma); + } + tablePIDKa(resolution, nSigma); + break; + case 4: + resolution = responsePr.GetExpectedSigma(mRespParamsV2, trkInColl); + nSigma = responsePr.GetSeparation(mRespParamsV2, trkInColl, resolution); + if (mcRecalib.enable) { + nSigma = applyMcRecalib(pidId, trkInColl.pt(), nSigma); + } + tablePIDPr(resolution, nSigma); + break; + default: + LOG(fatal) << "Wrong particle ID in processWSlice()"; + break; + } + } + } + } + } + PROCESS_SWITCH(mcPidTof, processWSlice, "Process with track slices", true); + + using TrksIU = soa::Join; + template + using ResponseImplementationIU = o2::pid::tof::ExpTimes; + void processWoSlice(TrksIU const& tracks, aod::Collisions const&, aod::BCsWithTimestamps const&) + { + constexpr auto responsePi = ResponseImplementationIU(); + constexpr auto responseKa = ResponseImplementationIU(); + constexpr auto responsePr = ResponseImplementationIU(); + + for (auto const& pidId : mEnabledParticles) { + reserveTable(pidId, tracks.size()); + } + float resolution = 1.f; // Last resolution assigned + for (auto const& track : tracks) { // Loop on all tracks + if (!track.has_collision()) { // Track was not assigned, cannot compute NSigma (no event time) -> filling with empty table + for (auto const& pidId : mEnabledParticles) { + makeTableEmpty(pidId); + } + continue; + } + + if (enableTimeDependentResponse && (track.collisionId() != mLastCollisionId)) { // Time dependent calib is enabled and this is a new collision + mLastCollisionId = track.collisionId(); // Cache last collision ID + timestamp.value = track.collision().bc_as().timestamp(); + LOG(debug) << "Updating parametrization from path '" << parametrizationPath.value << "' and timestamp " << timestamp.value; + if (!ccdb->getForTimeStamp(parametrizationPath.value, timestamp.value)->retrieveParameters(mRespParamsV2, passName.value)) { + if (fatalOnPassNotAvailable) { + LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } else { + LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } + } + } + + // in case of MC recalibrations, check if the objects from CCDB has to be updated + if (mcRecalib.enable) { + int runNumber = track.collision().bc_as().runNumber(); + if (runNumber != currentRun) { + // update postcalibration files + retrieveMcPostCalibFromCcdb(timestamp); + } + currentRun = runNumber; + } + + float nSigma{-999.f}; + for (auto const& pidId : mEnabledParticles) { // Loop on enabled particle hypotheses + switch (pidId) { + case 2: + resolution = responsePi.GetExpectedSigma(mRespParamsV2, track); + nSigma = responsePi.GetSeparation(mRespParamsV2, track, resolution); + if (mcRecalib.enable) { + nSigma = applyMcRecalib(pidId, track.pt(), nSigma); + } + tablePIDPi(resolution, nSigma); + break; + case 3: + resolution = responseKa.GetExpectedSigma(mRespParamsV2, track); + nSigma = responseKa.GetSeparation(mRespParamsV2, track, resolution); + if (mcRecalib.enable) { + nSigma = applyMcRecalib(pidId, track.pt(), nSigma); + } + tablePIDKa(resolution, nSigma); + break; + case 4: + resolution = responsePr.GetExpectedSigma(mRespParamsV2, track); + nSigma = responsePr.GetSeparation(mRespParamsV2, track, resolution); + if (mcRecalib.enable) { + nSigma = applyMcRecalib(pidId, track.pt(), nSigma); + } + tablePIDPr(resolution, nSigma); + break; + default: + LOG(fatal) << "Wrong particle ID in processWoSlice()"; + break; + } + } + } + } + PROCESS_SWITCH(mcPidTof, processWoSlice, "Process without track slices and on TrackIU (faster but only Run3)", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; } diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index d79d160401f..ff25b30e055 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -387,12 +387,14 @@ bool trackAcceptanceWithDca(T const& track, float trackDcaXYMax, float trackDcaZ * retrun acceptance of prong about chi2 and error of decay length due to cut for high quality secondary vertex */ template -bool prongAcceptance(T const& prong, float prongChi2PCAMin, float prongChi2PCAMax, float prongsigmaLxyMax, float prongIPxyMin, float prongIPxyMax, bool doXYZ) +bool prongAcceptance(T const& prong, float prongChi2PCAMin, float prongChi2PCAMax, float prongsigmaLxyMax, float prongIPxyMin, float prongIPxyMax, float prongDispersionMax, bool doXYZ) { if (prong.chi2PCA() < prongChi2PCAMin) return false; if (prong.chi2PCA() > prongChi2PCAMax) return false; + if (prong.dispersion() > prongDispersionMax) + return false; if (!doXYZ) { if (prong.errorDecayLengthXY() > prongsigmaLxyMax) return false; @@ -574,13 +576,13 @@ class bjetCandSV bjetCandSV(float xpv, float ypv, float zpv, float xsv, float ysv, float zsv, float pxVal, float pyVal, float pzVal, float eVal, float mVal, float chi2Val, - float errDecayLength, float errDecayLengthXY, + float dispersion, float errDecayLength, float errDecayLengthXY, float rSecVertex, float ptVal, float pVal, std::array pVec, float etaVal, float phiVal, float yVal, float decayLen, float decayLenXY, float decayLenNorm, float decayLenXYNorm, float cpaVal, float impParXY) - : m_xPVertex(xpv), m_yPVertex(ypv), m_zPVertex(zpv), m_xSecondaryVertex(xsv), m_ySecondaryVertex(ysv), m_zSecondaryVertex(zsv), m_px(pxVal), m_py(pyVal), m_pz(pzVal), m_e(eVal), m_m(mVal), m_chi2PCA(chi2Val), m_errorDecayLength(errDecayLength), m_errorDecayLengthXY(errDecayLengthXY), m_rSecondaryVertex(rSecVertex), m_pt(ptVal), m_p(pVal), m_pVector(pVec), m_eta(etaVal), m_phi(phiVal), m_y(yVal), m_decayLength(decayLen), m_decayLengthXY(decayLenXY), m_decayLengthNormalised(decayLenNorm), m_decayLengthXYNormalised(decayLenXYNorm), m_cpa(cpaVal), m_impactParameterXY(impParXY) + : m_xPVertex(xpv), m_yPVertex(ypv), m_zPVertex(zpv), m_xSecondaryVertex(xsv), m_ySecondaryVertex(ysv), m_zSecondaryVertex(zsv), m_px(pxVal), m_py(pyVal), m_pz(pzVal), m_e(eVal), m_m(mVal), m_chi2PCA(chi2Val), m_dispersion(dispersion), m_errorDecayLength(errDecayLength), m_errorDecayLengthXY(errDecayLengthXY), m_rSecondaryVertex(rSecVertex), m_pt(ptVal), m_p(pVal), m_pVector(pVec), m_eta(etaVal), m_phi(phiVal), m_y(yVal), m_decayLength(decayLen), m_decayLengthXY(decayLenXY), m_decayLengthNormalised(decayLenNorm), m_decayLengthXYNormalised(decayLenXYNorm), m_cpa(cpaVal), m_impactParameterXY(impParXY) { } @@ -598,6 +600,7 @@ class bjetCandSV float e() const { return m_e; } float m() const { return m_m; } float chi2PCA() const { return m_chi2PCA; } + float dispersion() const { return m_dispersion; } float errorDecayLength() const { return m_errorDecayLength; } float errorDecayLengthXY() const { return m_errorDecayLengthXY; } @@ -623,7 +626,7 @@ class bjetCandSV private: float m_xPVertex, m_yPVertex, m_zPVertex; float m_xSecondaryVertex, m_ySecondaryVertex, m_zSecondaryVertex; - float m_px, m_py, m_pz, m_e, m_m, m_chi2PCA; + float m_px, m_py, m_pz, m_e, m_m, m_chi2PCA, m_dispersion; float m_errorDecayLength, m_errorDecayLengthXY; float m_rSecondaryVertex, m_pt, m_p; std::array m_pVector; @@ -647,6 +650,7 @@ bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, float const& prongChi2 float e = 0.0f; float m = 0.0f; float chi2PCA = 0.0f; + float dispersion = 0.0f; float errorDecayLength = 0.0f; float errorDecayLengthXY = 0.0f; @@ -689,6 +693,7 @@ bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, float const& prongChi2 e = prong.e(); m = prong.m(); chi2PCA = prong.chi2PCA(); + dispersion = prong.dispersion(); errorDecayLength = prong.errorDecayLength(); errorDecayLengthXY = prong.errorDecayLengthXY(); rSecondaryVertex = prong.rSecondaryVertex(); @@ -711,7 +716,7 @@ bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, float const& prongChi2 xPVertex, yPVertex, zPVertex, xSecondaryVertex, ySecondaryVertex, zSecondaryVertex, px, py, pz, e, m, chi2PCA, - errorDecayLength, errorDecayLengthXY, + dispersion, errorDecayLength, errorDecayLengthXY, rSecondaryVertex, pt, p, pVector, eta, phi, y, decayLength, decayLengthXY, @@ -720,10 +725,10 @@ bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, float const& prongChi2 } template -bool isTaggedJetSV(T const jet, U const& /*prongs*/, float const& prongChi2PCAMin, float const& prongChi2PCAMax, float const& prongsigmaLxyMax, float const& prongIPxyMin, float const& prongIPxyMax, float const& doXYZ = false, float const& tagPointForSV = 15.) +bool isTaggedJetSV(T const jet, U const& /*prongs*/, float const& prongChi2PCAMin, float const& prongChi2PCAMax, float const& prongsigmaLxyMax, float const& prongIPxyMin, float const& prongIPxyMax, float prongDispersionMax, float const& doXYZ = false, float const& tagPointForSV = 15.) { auto bjetCand = jetFromProngMaxDecayLength(jet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, doXYZ); - if (!prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, doXYZ)) + if (!prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, doXYZ)) return false; if (!doXYZ) { auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index a034c8ca5c3..20183a2e149 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -49,6 +49,7 @@ struct JetTaggerHFTask { Configurable prongIPxyMax{"prongIpxyMax", 1, "minimum impact parmeter of prongs on xy plane [cm]"}; Configurable prongChi2PCAMin{"prongChi2PCAMin", 4, "minimum Chi2 PCA of decay length of prongs"}; Configurable prongChi2PCAMax{"prongChi2PCAMax", 100, "maximum Chi2 PCA of decay length of prongs"}; + Configurable prongDispersionMax{"prongDispersionMax", 1, "maximum dispersion of sv"}; // jet flavour definition Configurable maxDeltaR{"maxDeltaR", 0.25, "maximum distance of jet axis from flavour initiating parton"}; @@ -67,7 +68,7 @@ struct JetTaggerHFTask { Configurable> paramsResoFuncLfJetMC{"paramsResoFuncLfJetMC", std::vector{1539435.343, -0.061, 0.896, 13.272, 1.034, 5.884, 0.004, 7.843, 0.090}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; Configurable minSignImpXYSig{"minsIPs", -40.0, "minimum of signed impact parameter significance"}; Configurable tagPointForIP{"tagPointForIP", 2.5, "tagging working point for IP"}; - Configurable minIPCount{"minSipCount", 2, "Select at least N signed impact parameter significance in jets"}; // default 2 + Configurable minIPCount{"minIPCount", 2, "Select at least N signed impact parameter significance in jets"}; // default 2 // configuration about SV method Configurable doSV{"doSV", false, "fill table for secondary vertex algorithm"}; Configurable useXYZForTagging{"useXYZForTagging", false, "Enable tagging decision using full XYZ DCA for secondary vertex algorithm"}; @@ -270,9 +271,9 @@ struct JetTaggerHFTask { if (jettaggingutilities::isGreaterThanTaggingPoint(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) flagtaggedjetIP = true; if (!useXYZForTagging) { - flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, useXYZForTagging, tagPointForSV); + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongDispersionMax, useXYZForTagging, tagPointForSV); } else { - flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, useXYZForTagging, tagPointForSV); + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, prongDispersionMax, useXYZForTagging, tagPointForSV); } taggingTableData(0, jetProb, flagtaggedjetIP, flagtaggedjetSV); } @@ -323,9 +324,9 @@ struct JetTaggerHFTask { if (jettaggingutilities::isGreaterThanTaggingPoint(mcdjet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) flagtaggedjetIP = true; if (!useXYZForTagging) { - flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, prongIPxyMin, prongIPxyMax, useXYZForTagging, tagPointForSV); + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, useXYZForTagging, tagPointForSV); } else { - flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, prongIPxyMin, prongIPxyMax, useXYZForTagging, tagPointForSV); + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, useXYZForTagging, tagPointForSV); } taggingTableMCD(origin, jetProb, flagtaggedjetIP, flagtaggedjetSV); } diff --git a/PWGJE/Tasks/jetfragmentation.cxx b/PWGJE/Tasks/jetfragmentation.cxx index 2c6a4304958..d2512fc3890 100644 --- a/PWGJE/Tasks/jetfragmentation.cxx +++ b/PWGJE/Tasks/jetfragmentation.cxx @@ -536,7 +536,7 @@ struct JetFragmentation { registry.add("matching/jets/missPartJetPtZTheta", "Misses", HistType::kTH3D, {partJetPtAxis, partZAxis, partThetaAxis}); } // doprocessMcMatched - if (doprocessMcMatchedV0 || doprocessMcMatchedV0Frag || doprocessMcMatchedV0JetsFrag || doprocessMcV0PerpCone) { + if (doprocessMcMatchedV0 || doprocessMcMatchedV0Frag || doprocessMcMatchedV0JetsFrag || doprocessMcV0MatchedPerpCone) { registry.add("matching/V0/nV0sEvent", "nV0sDet per event", HistType::kTH1D, {v0Count}); registry.add("matching/V0/nV0sEventWeighted", "nV0sDet per event (weighted)", HistType::kTH1D, {v0Count}); registry.get(HIST("matching/V0/nV0sEventWeighted"))->Sumw2(); @@ -838,6 +838,35 @@ struct JetFragmentation { } // doprocessDataV0PerpCone if (doprocessMcV0PerpCone) { + registry.add("mcd/PC/jetPtEtaFakeV0Pt", "JetPtEtaFakeV0Pt", HistType::kTH3D, {detJetPtAxis, etaAxis, V0PtAxis}); + registry.add("mcd/PC/fakeV0PtEtaPhi", "fakeV0PtEtaPhi", HistType::kTH3D, {V0PtAxis, V0EtaAxis, V0PhiAxis}); + registry.add("mcd/PC/fakeV0PtCtau", "fakeV0PtCtau", HistType::kTHnSparseD, {V0PtAxis, V0CtauAxis, V0CtauAxis, V0CtauAxis}); + registry.add("mcd/PC/fakeV0PtMass", "fakeV0PtMass", HistType::kTHnSparseD, {V0PtAxis, K0SMassAxis, LambdaMassAxis, LambdaMassAxis}); + registry.add("mcd/PC/fakeV0PtRadiusCosPA", "fakeV0PtRadiusCosPA", HistType::kTH3D, {V0PtAxis, V0RadiusAxis, V0CosPAAxis}); + registry.add("mcd/PC/fakeV0PtDCAposneg", "fakeV0PtDCAposneg", HistType::kTH3D, {V0PtAxis, V0DCApAxis, V0DCAnAxis}); + registry.add("mcd/PC/fakeV0PtDCAd", "fakeV0PtDCAd", HistType::kTH2D, {V0PtAxis, V0DCAdAxis}); + registry.add("mcd/PC/jetPtEtaMatchedV0Pt", "JetPtEtaMatchedV0Pt", HistType::kTH3D, {detJetPtAxis, etaAxis, V0PtAxis}); + + registry.add("mcd/PC/matchedV0PtEtaPhi", "matchedV0PtEtaPhi", HistType::kTH3D, {V0PtAxis, V0EtaAxis, V0PhiAxis}); + registry.add("mcd/PC/matchedV0PtCtau", "matchedV0PtCtau", HistType::kTHnSparseD, {V0PtAxis, V0CtauAxis, V0CtauAxis, V0CtauAxis}); + registry.add("mcd/PC/matchedV0PtMass", "matchedV0PtMass", HistType::kTHnSparseD, {V0PtAxis, K0SMassAxis, LambdaMassAxis, LambdaMassAxis}); + registry.add("mcd/PC/matchedV0PtRadiusCosPA", "matchedV0PtRadiusCosPA", HistType::kTH3D, {V0PtAxis, V0RadiusAxis, V0CosPAAxis}); + registry.add("mcd/PC/matchedV0PtDCAposneg", "matchedV0PtDCAposneg", HistType::kTH3D, {V0PtAxis, V0DCApAxis, V0DCAnAxis}); + registry.add("mcd/PC/matchedV0PtDCAd", "matchedV0PtDCAd", HistType::kTH2D, {V0PtAxis, V0DCAdAxis}); + + registry.add("mcd/PC/matchedJetPtK0SPtMass", "matchedJetPtK0SPtMass", HistType::kTH3D, {detJetPtAxis, V0PtAxis, K0SMassAxis}); + registry.add("mcd/PC/matchedJetPtLambda0PtMass", "matchedJetPtLambda0PtMass", HistType::kTH3D, {detJetPtAxis, V0PtAxis, LambdaMassAxis}); + registry.add("mcd/PC/matchedJetPtAntiLambda0PtMass", "matchedJetPtAntiLambda0PtMass", HistType::kTH3D, {detJetPtAxis, V0PtAxis, LambdaMassAxis}); + registry.add("mcd/PC/matchednV0sConePtEta", "matchednV0sConePtEta", HistType::kTH3D, {v0Count, detJetPtAxis, etaAxis}); + registry.add("mcd/PC/matchedConePtEtaPhi", "matchedConePtEtaPhi", HistType::kTH3D, {detJetPtAxis, etaAxis, phiAxis}); + registry.add("mcd/PC/matchedJetPtEtaConePt", "matchedJetPtEtaConePt", HistType::kTH3D, {detJetPtAxis, etaAxis, detJetPtAxis}); + + registry.add("mcd/PC/fakenV0sConePtEta", "fakenV0sConePtEta", HistType::kTH3D, {v0Count, detJetPtAxis, etaAxis}); + registry.add("mcd/PC/fakeConePtEtaPhi", "fakeConePtEtaPhi", HistType::kTH3D, {detJetPtAxis, etaAxis, phiAxis}); + registry.add("mcd/PC/fakeJetPtEtaConePt", "fakeJetPtEtaConePt", HistType::kTH3D, {detJetPtAxis, etaAxis, detJetPtAxis}); + } // doprocessMcV0PerpCone + + if (doprocessMcV0MatchedPerpCone) { registry.add("matching/PC/jetPtEtaFakeV0Pt", "JetPtEtaFakeV0Pt", HistType::kTH3D, {detJetPtAxis, etaAxis, V0PtAxis}); registry.add("matching/PC/jetsPtFakeV0Pt", "jetsPtFakeV0Pt", HistType::kTH3D, {partJetPtAxis, detJetPtAxis, V0PtAxis}); registry.add("matching/PC/fakeV0PtEtaPhi", "fakeV0PtEtaPhi", HistType::kTH3D, {V0PtAxis, V0EtaAxis, V0PhiAxis}); @@ -872,7 +901,7 @@ struct JetFragmentation { registry.add("matching/PC/fakeConePtEtaPhi", "fakeConePtEtaPhi", HistType::kTH3D, {detJetPtAxis, etaAxis, phiAxis}); registry.add("matching/PC/fakeJetPtEtaConePt", "fakeJetPtEtaConePt", HistType::kTH3D, {detJetPtAxis, etaAxis, detJetPtAxis}); registry.add("matching/PC/fakeJetsPtEtaConePt", "fakeJetsPtEtaConePt", HistType::kTHnSparseD, {partJetPtAxis, detJetPtAxis, etaAxis, detJetPtAxis}); - } // doprocessMcV0PerpCone + } // doprocessMcV0MatchedPerpCone } // init // TODO: This should contain a lookup table or function containing the various V0 weights @@ -1957,6 +1986,76 @@ struct JetFragmentation { } } + // Version for MCD jets + template + void fillMcPerpConeHists(T const& coll, U const& mcdjet, V const& v0s, W const& /* V0 particles */, double weight = 1.) + { + double perpConeR = mcdjet.r() * 1e-2; + double conePhi[2] = {RecoDecay::constrainAngle(mcdjet.phi() - M_PI / 2, -M_PI), + RecoDecay::constrainAngle(mcdjet.phi() + M_PI / 2, -M_PI)}; + double coneMatchedPt[2] = {0., 0.}; + double coneFakePt[2] = {0., 0.}; + int nMatchedV0sinCone[2] = {0, 0}; + int nFakeV0sinCone[2] = {0, 0}; + + for (const auto& v0 : v0s) { + double dEta = v0.eta() - mcdjet.eta(); + double dPhi[2] = {RecoDecay::constrainAngle(v0.phi() - conePhi[0], -M_PI), + RecoDecay::constrainAngle(v0.phi() - conePhi[1], -M_PI)}; + for (int i = 0; i < 2; i++) { + if (TMath::Sqrt(dEta * dEta + dPhi[i] * dPhi[i]) > perpConeR) { + continue; + } + + double ctauLambda = v0.distovertotmom(coll.posX(), coll.posY(), coll.posZ()) * o2::constants::physics::MassLambda0; + double ctauAntiLambda = v0.distovertotmom(coll.posX(), coll.posY(), coll.posZ()) * o2::constants::physics::MassLambda0Bar; + double ctauK0s = v0.distovertotmom(coll.posX(), coll.posY(), coll.posZ()) * o2::constants::physics::MassK0Short; + + if (!v0.has_mcParticle()) { // The V0 is combinatorial background + coneFakePt[i] += v0.pt(); + nFakeV0sinCone[i]++; + registry.fill(HIST("mcd/PC/jetPtEtaFakeV0Pt"), mcdjet.pt(), mcdjet.eta(), v0.pt(), weight); + + registry.fill(HIST("mcd/PC/fakeV0PtEtaPhi"), v0.pt(), v0.eta(), v0.phi(), weight); + registry.fill(HIST("mcd/PC/fakeV0PtCtau"), v0.pt(), ctauK0s, ctauLambda, ctauAntiLambda, weight); + registry.fill(HIST("mcd/PC/fakeV0PtMass"), v0.pt(), v0.mK0Short(), v0.mLambda(), v0.mAntiLambda(), weight); + registry.fill(HIST("mcd/PC/fakeV0PtRadiusCosPA"), v0.pt(), v0.v0radius(), v0.v0cosPA(), weight); + registry.fill(HIST("mcd/PC/fakeV0PtDCAposneg"), v0.pt(), v0.dcapostopv(), v0.dcanegtopv(), weight); + registry.fill(HIST("mcd/PC/fakeV0PtDCAd"), v0.pt(), v0.dcaV0daughters(), weight); + } else { + coneMatchedPt[i] += v0.pt(); + nMatchedV0sinCone[i]++; + registry.fill(HIST("mcd/PC/jetPtEtaMatchedV0Pt"), mcdjet.pt(), mcdjet.eta(), v0.pt(), weight); + + registry.fill(HIST("mcd/PC/matchedV0PtEtaPhi"), v0.pt(), v0.eta(), v0.phi(), weight); + registry.fill(HIST("mcd/PC/matchedV0PtCtau"), v0.pt(), ctauK0s, ctauLambda, ctauAntiLambda, weight); + registry.fill(HIST("mcd/PC/matchedV0PtMass"), v0.pt(), v0.mK0Short(), v0.mLambda(), v0.mAntiLambda(), weight); + registry.fill(HIST("mcd/PC/matchedV0PtRadiusCosPA"), v0.pt(), v0.v0radius(), v0.v0cosPA(), weight); + registry.fill(HIST("mcd/PC/matchedV0PtDCAposneg"), v0.pt(), v0.dcapostopv(), v0.dcanegtopv(), weight); + registry.fill(HIST("mcd/PC/matchedV0PtDCAd"), v0.pt(), v0.dcaV0daughters(), weight); + + auto particle = v0.template mcParticle_as(); + if (TMath::Abs(particle.pdgCode()) == 310) { // K0S + registry.fill(HIST("mcd/PC/matchedJetPtK0SPtMass"), mcdjet.pt(), v0.pt(), v0.mK0Short(), weight); + } else if (particle.pdgCode() == 3122) { // Lambda + registry.fill(HIST("mcd/PC/matchedJetPtLambda0PtMass"), mcdjet.pt(), v0.pt(), v0.mLambda(), weight); + } else if (particle.pdgCode() == -3122) { + registry.fill(HIST("mcd/PC/matchedJetPtAntiLambda0PtMass"), mcdjet.pt(), v0.pt(), v0.mAntiLambda(), weight); + } + } // if v0 has mcParticle + } // for cone + } // for v0s + for (int i = 0; i < 2; i++) { + registry.fill(HIST("mcd/PC/matchednV0sConePtEta"), nMatchedV0sinCone[i], coneMatchedPt[i], mcdjet.eta(), weight); + registry.fill(HIST("mcd/PC/matchedConePtEtaPhi"), coneMatchedPt[i], mcdjet.eta(), conePhi[i], weight); + registry.fill(HIST("mcd/PC/matchedJetPtEtaConePt"), mcdjet.pt(), mcdjet.eta(), coneMatchedPt[i], weight); + + registry.fill(HIST("mcd/PC/fakenV0sConePtEta"), nFakeV0sinCone[i], coneFakePt[i], mcdjet.eta(), weight); + registry.fill(HIST("mcd/PC/fakeConePtEtaPhi"), coneFakePt[i], mcdjet.eta(), conePhi[i], weight); + registry.fill(HIST("mcd/PC/fakeJetPtEtaConePt"), mcdjet.pt(), mcdjet.eta(), coneFakePt[i], weight); + } + } + // Version for matched jets template void fillMcPerpConeHists(T const& coll, U const& mcdjet, V const& mcpjet, W const& v0s, X const& /* V0 particles */, double weight = 1.) { @@ -2550,11 +2649,14 @@ struct JetFragmentation { } PROCESS_SWITCH(JetFragmentation, processDataV0JetsFragWithWeights, "Data V0 jets fragmentation with weights", false); - void processDataV0PerpCone(soa::Filtered::iterator const& jcoll, soa::Join const& v0jets, CandidatesV0Data const& v0s) + void processDataV0PerpCone(soa::Filtered::iterator const& jcoll, aod::V0ChargedJets const& v0jets, CandidatesV0Data const& v0s) { if (!jetderiveddatautilities::selectCollision(jcoll, eventSelection)) { return; } + if (v0s.size() == 0) { + return; + } registry.fill(HIST("data/V0/nV0sEvent"), v0s.size()); fillDataV0Histograms(jcoll, v0s); @@ -2689,38 +2791,50 @@ struct JetFragmentation { } PROCESS_SWITCH(JetFragmentation, processMcMatchedV0JetsFrag, "Matched V0 jets fragmentation", false); - void processMcV0PerpCone(soa::Filtered::iterator const& jcoll, JetMcCollisions const&, MatchedMCDV0Jets const& v0jets, MatchedMCPV0Jets const&, MatchedMCDJets const& chjets, MatchedMCPJets const&, soa::Join const& v0s, aod::McParticles const& particles) + void processMcV0PerpCone(soa::Filtered::iterator const& jcoll, JetMcCollisions const&, MatchedMCDV0Jets const& v0jets, soa::Join const& v0s, aod::McParticles const& particles) { if (!jetderiveddatautilities::selectCollision(jcoll, eventSelection)) { return; } + if (v0s.size() == 0) { + return; + } + double weight = jcoll.mcCollision().weight(); + registry.fill(HIST("mcd/V0/nV0sEvent"), v0s.size()); + registry.fill(HIST("mcd/V0/nV0sEventWeighted"), v0s.size(), weight); + + for (const auto& mcdjet : v0jets) { + if (!jetfindingutilities::isInEtaAcceptance(mcdjet, -99., -99., v0EtaMin, v0EtaMax)) { + continue; + } + fillMcPerpConeHists(jcoll, mcdjet, v0s, particles, weight); + } + } + PROCESS_SWITCH(JetFragmentation, processMcV0PerpCone, "Perpendicular cone V0s in MC", false); + + void processMcV0MatchedPerpCone(soa::Filtered::iterator const& jcoll, JetMcCollisions const&, MatchedMCDV0Jets const& v0jets, MatchedMCPV0Jets const&, soa::Join const& v0s, aod::McParticles const& particles) + { + if (!jetderiveddatautilities::selectCollision(jcoll, eventSelection)) { + return; + } + if (v0s.size() == 0) { + return; + } double weight = jcoll.mcCollision().weight(); registry.fill(HIST("matching/V0/nV0sEvent"), v0s.size()); registry.fill(HIST("matching/V0/nV0sEventWeighted"), v0s.size(), weight); - // For events withouth V0s, we need to fill the perp cone histograms using charged jets - if (v0s.size() > 0) { - for (const auto& mcdjet : v0jets) { - if (!jetfindingutilities::isInEtaAcceptance(mcdjet, -99., -99., v0EtaMin, v0EtaMax)) { - continue; - } - for (const auto& mcpjet : mcdjet.template matchedJetGeo_as()) { - fillMcPerpConeHists(jcoll, mcdjet, mcpjet, v0s, particles, weight); - break; // Make sure we only do this once - } + for (const auto& mcdjet : v0jets) { + if (!jetfindingutilities::isInEtaAcceptance(mcdjet, -99., -99., v0EtaMin, v0EtaMax)) { + continue; } - } else { - for (const auto& mcdjet : chjets) { - if (!jetfindingutilities::isInEtaAcceptance(mcdjet, -99., -99., v0EtaMin, v0EtaMax)) { - continue; - } - for (const auto& mcpjet : mcdjet.template matchedJetGeo_as()) { - fillMcPerpConeHists(jcoll, mcdjet, mcpjet, v0s, particles, weight); - } + for (const auto& mcpjet : mcdjet.template matchedJetGeo_as()) { + fillMcPerpConeHists(jcoll, mcdjet, mcpjet, v0s, particles, weight); + break; // Make sure we only do this once } } } - PROCESS_SWITCH(JetFragmentation, processMcV0PerpCone, "Perpendicular cone V0s in MC", false); + PROCESS_SWITCH(JetFragmentation, processMcV0MatchedPerpCone, "Perpendicular cone V0s in MC, matched jets", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index 37d925313c7..56ab3075ece 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -60,7 +60,8 @@ struct JetTaggerHFQA { Configurable prongsigmaLxyMax{"prongsigmaLxyMax", 100, "maximum sigma of decay length of prongs on xy plane"}; Configurable prongsigmaLxyzMax{"prongsigmaLxyzMax", 100, "maximum sigma of decay length of prongs on xyz plane"}; Configurable prongIPxyMin{"prongIPxyMin", 0.008, "maximum impact paramter of prongs on xy plane"}; - Configurable prongIPxyMax{"prongIpxyMax", 1, "minimum impact parmeter of prongs on xy plane"}; + Configurable prongIPxyMax{"prongIPxyMax", 1, "minimum impact parmeter of prongs on xy plane"}; + Configurable prongDispersionMax{"prongDispersionMax", 1, "maximum dispersion of sv"}; Configurable numFlavourSpecies{"numFlavourSpecies", 6, "number of jet flavour species"}; Configurable numOrder{"numOrder", 6, "number of ordering"}; Configurable pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"}; @@ -68,6 +69,7 @@ struct JetTaggerHFQA { Configurable pTHatExponent{"pTHatExponent", 6.0, "exponent of the event weight for the calculation of pTHat"}; Configurable jetAreaFractionMin{"jetAreaFractionMin", -99.0, "used to make a cut on the jet areas"}; Configurable leadingConstituentPtMin{"leadingConstituentPtMin", -99.0, "minimum pT selection on jet constituent"}; + Configurable leadingConstituentPtMax{"leadingConstituentPtMax", 9999.0, "maximum pT selection on jet constituent"}; Configurable checkMcCollisionIsMatched{"checkMcCollisionIsMatched", false, "0: count whole MCcollisions, 1: select MCcollisions which only have their correspond collisions"}; Configurable trackOccupancyInTimeRangeMax{"trackOccupancyInTimeRangeMax", 999999, "maximum occupancy of tracks in neighbouring collisions in a given time range; only applied to reconstructed collisions (data and mcd jets), not mc collisions (mcp jets)"}; Configurable trackOccupancyInTimeRangeMin{"trackOccupancyInTimeRangeMin", -999999, "minimum occupancy of tracks in neighbouring collisions in a given time range; only applied to reconstructed collisions (data and mcd jets), not mc collisions (mcp jets)"}; @@ -292,6 +294,9 @@ struct JetTaggerHFQA { registry.add("h2_jet_pt_2prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); registry.add("h2_jet_pt_2prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); registry.add("h2_jet_pt_2prong_mass_N1", "", {HistType::kTH2F, {{jetPtAxis}, {massAxis}}}); + registry.add("h2_taggedjet_pt_2prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_taggedjet_pt_2prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_taggedjet_pt_2prong_mass_N1", "", {HistType::kTH2F, {{jetPtAxis}, {massAxis}}}); } if (doprocessSV3ProngData) { registry.add("h_3prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); @@ -304,6 +309,9 @@ struct JetTaggerHFQA { registry.add("h2_jet_pt_3prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); registry.add("h2_jet_pt_3prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); registry.add("h2_jet_pt_3prong_mass_N1", "", {HistType::kTH2F, {{jetPtAxis}, {massAxis}}}); + registry.add("h2_taggedjet_pt_3prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_taggedjet_pt_3prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_taggedjet_pt_3prong_mass_N1", "", {HistType::kTH2F, {{jetPtAxis}, {massAxis}}}); } if (doprocessSV2ProngMCD || doprocessSV2ProngMCDWeighted) { registry.add("h2_2prong_nprongs_flavour", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); @@ -383,18 +391,30 @@ struct JetTaggerHFQA { return false; } } - if (leadingConstituentPtMin > -98.0) { - bool isMinleadingConstituent = false; - for (auto& constituent : jet.template tracks_as()) { - if (constituent.pt() >= leadingConstituentPtMin) { - isMinleadingConstituent = true; - break; + bool checkConstituentPt = true; + bool checkConstituentMinPt = (leadingConstituentPtMin > -98.0); + bool checkConstituentMaxPt = (leadingConstituentPtMax < 9998.0); + if (!checkConstituentMinPt && !checkConstituentMaxPt) { + checkConstituentPt = false; + } + + if (checkConstituentPt) { + bool isMinLeadingConstituent = !checkConstituentMinPt; + bool isMaxLeadingConstituent = true; + + for (const auto& constituent : jet.template tracks_as()) { + double pt = constituent.pt(); + + if (checkConstituentMinPt && pt >= leadingConstituentPtMin) { + isMinLeadingConstituent = true; + } + if (checkConstituentMaxPt && pt > leadingConstituentPtMax) { + isMaxLeadingConstituent = false; } } - if (!isMinleadingConstituent) { - return false; - } + return isMinLeadingConstituent && isMaxLeadingConstituent; } + return true; } @@ -800,14 +820,20 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_jet_pt_2prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); } auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax); - if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false)) { + if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, false)) { return; } auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + auto massSV = bjetCand.m(); auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, prongIPxyMin, prongIPxyMax, true); auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h2_jet_pt_2prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_2prong_Sxyz_N1"), jet.pt(), maxSxyz); + if (!jet.flagtaggedjetSV()) + return; + registry.fill(HIST("h2_taggedjet_pt_2prong_Sxy_N1"), jet.pt(), maxSxy); + registry.fill(HIST("h2_taggedjet_pt_2prong_Sxyz_N1"), jet.pt(), maxSxyz); + registry.fill(HIST("h2_taggedjet_pt_2prong_mass_N1"), jet.pt(), massSV); } template @@ -832,14 +858,20 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_jet_pt_3prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); } auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax); - if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false)) { + if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, false)) { return; } auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + auto massSV = bjetCand.m(); auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyzMax, prongIPxyMin, prongIPxyMax, true); auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h2_jet_pt_3prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_3prong_Sxyz_N1"), jet.pt(), maxSxyz); + if (!jet.flagtaggedjetSV()) + return; + registry.fill(HIST("h2_taggedjet_pt_3prong_Sxy_N1"), jet.pt(), maxSxy); + registry.fill(HIST("h2_taggedjet_pt_3prong_Sxyz_N1"), jet.pt(), maxSxyz); + registry.fill(HIST("h2_taggedjet_pt_3prong_mass_N1"), jet.pt(), massSV); } template @@ -866,7 +898,7 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_2prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin, eventWeight); } auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax); - if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false)) { + if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, false)) { return; } auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); @@ -911,7 +943,7 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_2prong_sigmaLxyz_flavour_run2"), mcdjet.pt(), prong.errorDecayLength(), jetflavourRun2Def, eventWeight); } auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax); - if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false)) { + if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, false)) { return; } auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); @@ -926,6 +958,7 @@ struct JetTaggerHFQA { template void fillHistogramSV3ProngMCD(T const& mcdjet, U const& /*prongs*/, float eventWeight = 1.0) { + // std::cout << "weight: " << eventWeight << std::endl; float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent)); if (mcdjet.pt() > pTHatMaxMCD * pTHat) { return; @@ -947,7 +980,7 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_3prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin, eventWeight); } auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax); - if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false)) { + if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, false)) { return; } auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); @@ -992,7 +1025,7 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_3prong_sigmaLxyz_flavour_run2"), mcdjet.pt(), prong.errorDecayLength(), jetflavourRun2Def, eventWeight); } auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax); - if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false)) { + if (!jettaggingutilities::prongAcceptance(bjetCand, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, prongDispersionMax, false)) { return; } @@ -1168,7 +1201,7 @@ struct JetTaggerHFQA { } PROCESS_SWITCH(JetTaggerHFQA, processIPsMCPMCDMatchedWeighted, "Fill impact parameter imformation for mcp mcd matched jets", false); - void processJPData(soa::Filtered::iterator const& collision, soa::Join const& jets, JetTagTracksData const&) + void processJPData(soa::Filtered::iterator const& collision, soa::Join const& jets, JetTagTracksData const&) { if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) { return; diff --git a/PWGLF/DataModel/SPCalibrationTables.h b/PWGLF/DataModel/SPCalibrationTables.h index 82d61ddacad..80de9435169 100644 --- a/PWGLF/DataModel/SPCalibrationTables.h +++ b/PWGLF/DataModel/SPCalibrationTables.h @@ -31,6 +31,9 @@ namespace spcalibrationtable DECLARE_SOA_COLUMN(TriggerEvent, triggerevent, bool); DECLARE_SOA_COLUMN(TriggerEventRunNo, triggereventrunno, int); DECLARE_SOA_COLUMN(Cent, cent, float); +DECLARE_SOA_COLUMN(Vx, vx, float); +DECLARE_SOA_COLUMN(Vy, vy, float); +DECLARE_SOA_COLUMN(Vz, vz, float); DECLARE_SOA_COLUMN(QXZDCA, qxZDCA, float); DECLARE_SOA_COLUMN(QXZDCC, qxZDCC, float); DECLARE_SOA_COLUMN(QYZDCA, qyZDCA, float); @@ -42,6 +45,9 @@ DECLARE_SOA_TABLE(SPCalibrationTables, "AOD", "SPCALCOLS", spcalibrationtable::TriggerEvent, spcalibrationtable::TriggerEventRunNo, spcalibrationtable::Cent, + spcalibrationtable::Vx, + spcalibrationtable::Vy, + spcalibrationtable::Vz, spcalibrationtable::QXZDCA, spcalibrationtable::QXZDCC, spcalibrationtable::QYZDCA, diff --git a/PWGLF/DataModel/v0qaanalysis.h b/PWGLF/DataModel/v0qaanalysis.h index 05c90791bde..9e4efd7c3c6 100644 --- a/PWGLF/DataModel/v0qaanalysis.h +++ b/PWGLF/DataModel/v0qaanalysis.h @@ -58,6 +58,8 @@ DECLARE_SOA_COLUMN(IsPhysicalPrimary, isphysprimary, bool); DECLARE_SOA_COLUMN(MultFT0M, multft0m, float); DECLARE_SOA_COLUMN(MultFV0A, multfv0a, float); DECLARE_SOA_COLUMN(EvFlag, evflag, int); +DECLARE_SOA_COLUMN(Alpha, alpha, float); +DECLARE_SOA_COLUMN(QtArm, qtarm, float); } // namespace myv0candidates @@ -73,7 +75,7 @@ DECLARE_SOA_TABLE(MyV0Candidates, "AOD", "MYV0CANDIDATES", o2::soa::Index<>, myv0candidates::PosHasTOF, myv0candidates::NegHasTOF, myv0candidates::PDGCode, myv0candidates::IsPhysicalPrimary, myv0candidates::MultFT0M, myv0candidates::MultFV0A, - myv0candidates::EvFlag); + myv0candidates::EvFlag, myv0candidates::Alpha, myv0candidates::QtArm); } // namespace o2::aod diff --git a/PWGLF/TableProducer/Common/spvector.cxx b/PWGLF/TableProducer/Common/spvector.cxx index 3eb58b66514..dea10b4ce75 100644 --- a/PWGLF/TableProducer/Common/spvector.cxx +++ b/PWGLF/TableProducer/Common/spvector.cxx @@ -127,6 +127,11 @@ struct spvector { Configurable QA{"QA", false, "QA histograms"}; Configurable ispolarization{"ispolarization", false, "Flag to check polarization"}; Configurable finecorrection{"finecorrection", false, "Flag to check fine correction"}; + Configurable rejbadevent{"rejbadevent", true, "Flag to check bad events"}; + Configurable rejbadeventcent{"rejbadeventcent", true, "Flag to check bad events for cent"}; + Configurable rejbadeventvx{"rejbadeventvx", true, "Flag to check bad events for vx"}; + Configurable rejbadeventvy{"rejbadeventvy", true, "Flag to check bad events for vy"}; + Configurable rejbadeventvz{"rejbadeventvz", true, "Flag to check bad events for vz"}; Configurable usesparse{"usesparse", false, "flag to use sparse histogram"}; Configurable usenormqn{"usenormqn", true, "flag to use normalized qs"}; Configurable refsys{"refsys", true, "flag to use own reference system"}; @@ -135,6 +140,7 @@ struct spvector { Configurable useRecenterefineSp{"useRecenterefineSp", false, "use fine Recentering with Sparse or THn"}; Configurable useRecenteresqSp{"useRecenteresqSp", false, "use Recenteringsq with Sparse or THn"}; Configurable recwitherror{"recwitherror", false, "use Recentering with error"}; + Configurable recfinewitherror{"recfinewitherror", false, "use Recentering fine with error"}; Configurable ConfGainPath{"ConfGainPath", "Users/p/prottay/My/Object/NewPbPbpass4_10092024/gaincallib", "Path to gain calibration"}; Configurable ConfRecentereSp{"ConfRecentereSp", "Users/p/prottay/My/Object/Testingwithsparse/NewPbPbpass4_17092024/recenter", "Sparse or THn Path for recentere"}; Configurable ConfRecenterecentSp{"ConfRecenterecentSp", "Users/p/prottay/My/Object/Testingwithsparse/NewPbPbpass4_17092024/recenter", "Sparse or THn Path for cent recentere"}; @@ -309,7 +315,7 @@ struct spvector { auto bc = collision.foundBC_as(); if (!bc.has_zdc()) { triggerevent = false; - spcalibrationtable(triggerevent, currentRunNumber, centrality, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + spcalibrationtable(triggerevent, currentRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); return; } @@ -319,12 +325,12 @@ struct spvector { if (znaEnergy[0] < 0.0 || znaEnergy[1] < 0.0 || znaEnergy[2] < 0.0 || znaEnergy[3] < 0.0) { triggerevent = false; - spcalibrationtable(triggerevent, currentRunNumber, centrality, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + spcalibrationtable(triggerevent, currentRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); return; } if (zncEnergy[0] < 0.0 || zncEnergy[1] < 0.0 || zncEnergy[2] < 0.0 || zncEnergy[3] < 0.0) { triggerevent = false; - spcalibrationtable(triggerevent, currentRunNumber, centrality, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + spcalibrationtable(triggerevent, currentRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); return; } @@ -350,7 +356,7 @@ struct spvector { if (znaEnergy[iChA] <= 0.0) { triggerevent = false; - spcalibrationtable(triggerevent, currentRunNumber, centrality, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + spcalibrationtable(triggerevent, currentRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); return; } else { float ampl = gainequal * znaEnergy[iChA]; @@ -363,7 +369,7 @@ struct spvector { } else { if (zncEnergy[iChA - 4] <= 0.0) { triggerevent = false; - spcalibrationtable(triggerevent, currentRunNumber, centrality, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + spcalibrationtable(triggerevent, currentRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); return; } else { float ampl = gainequal * zncEnergy[iChA - 4]; @@ -401,7 +407,7 @@ struct spvector { qyZDCA = 0.0; qyZDCC = 0.0; triggerevent = false; - spcalibrationtable(triggerevent, currentRunNumber, centrality, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + spcalibrationtable(triggerevent, currentRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); return; } @@ -462,6 +468,14 @@ struct spvector { float meanyC = hrecentereSp->GetBinContent(globalBinMeanyC); float meanyCerror = hrecentereSp->GetBinError(globalBinMeanyC); + if (rejbadevent) { + if ((TMath::Abs(meanxA) > 90000.0 || TMath::Abs(meanxC) > 90000.0 || TMath::Abs(meanyA) > 90000.0 || TMath::Abs(meanyC) > 90000.0) && (TMath::Abs(meanxAerror) > 9000.0 || TMath::Abs(meanxCerror) > 9000.0 || TMath::Abs(meanyAerror) > 9000.0 || TMath::Abs(meanyCerror) > 9000.0)) { + triggerevent = false; + spcalibrationtable(triggerevent, lastRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + return; + } + } + qxZDCA = qxZDCA - meanxA; qyZDCA = qyZDCA - meanyA; qxZDCC = qxZDCC - meanxC; @@ -522,31 +536,135 @@ struct spvector { } if (useRecenterefineSp && hrecenterecentSp) { + + float meanxAcent = hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 0.5)); + float meanyAcent = hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 1.5)); + float meanxCcent = hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 2.5)); + float meanyCcent = hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 3.5)); + + float meanxAcenterror = hrecenterecentSp->GetBinError(hrecenterecentSp->FindBin(centrality, 0.5)); + float meanyAcenterror = hrecenterecentSp->GetBinError(hrecenterecentSp->FindBin(centrality, 1.5)); + float meanxCcenterror = hrecenterecentSp->GetBinError(hrecenterecentSp->FindBin(centrality, 2.5)); + float meanyCcenterror = hrecenterecentSp->GetBinError(hrecenterecentSp->FindBin(centrality, 3.5)); + + if (rejbadeventcent) { + if ((TMath::Abs(meanxAcent) > 90000.0 || TMath::Abs(meanxCcent) > 90000.0 || TMath::Abs(meanyAcent) > 90000.0 || TMath::Abs(meanyCcent) > 90000.0) && (TMath::Abs(meanxAcenterror) > 9000.0 || TMath::Abs(meanxCcenterror) > 9000.0 || TMath::Abs(meanyAcenterror) > 9000.0 || TMath::Abs(meanyCcenterror) > 9000.0)) { + triggerevent = false; + spcalibrationtable(triggerevent, lastRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + return; + } + } + qxZDCA = qxZDCA - hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 0.5)); qyZDCA = qyZDCA - hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 1.5)); qxZDCC = qxZDCC - hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 2.5)); qyZDCC = qyZDCC - hrecenterecentSp->GetBinContent(hrecenterecentSp->FindBin(centrality, 3.5)); + + if (recfinewitherror) { + qxZDCA = qxZDCA / meanxAcenterror; + qyZDCA = qyZDCA / meanyAcenterror; + qxZDCC = qxZDCC / meanxCcenterror; + qyZDCC = qyZDCC / meanyCcenterror; + } } if (useRecenterefineSp && hrecenterevxSp) { + + float meanxAvx = hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 0.5)); + float meanyAvx = hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 1.5)); + float meanxCvx = hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 2.5)); + float meanyCvx = hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 3.5)); + + float meanxAvxerror = hrecenterevxSp->GetBinError(hrecenterevxSp->FindBin(vx, 0.5)); + float meanyAvxerror = hrecenterevxSp->GetBinError(hrecenterevxSp->FindBin(vx, 1.5)); + float meanxCvxerror = hrecenterevxSp->GetBinError(hrecenterevxSp->FindBin(vx, 2.5)); + float meanyCvxerror = hrecenterevxSp->GetBinError(hrecenterevxSp->FindBin(vx, 3.5)); + + if (rejbadeventvx) { + if ((TMath::Abs(meanxAvx) > 90000.0 || TMath::Abs(meanxCvx) > 90000.0 || TMath::Abs(meanyAvx) > 90000.0 || TMath::Abs(meanyCvx) > 90000.0) && (TMath::Abs(meanxAvxerror) > 9000.0 || TMath::Abs(meanxCvxerror) > 9000.0 || TMath::Abs(meanyAvxerror) > 9000.0 || TMath::Abs(meanyCvxerror) > 9000.0)) { + triggerevent = false; + spcalibrationtable(triggerevent, lastRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + return; + } + } + qxZDCA = qxZDCA - hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 0.5)); qyZDCA = qyZDCA - hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 1.5)); qxZDCC = qxZDCC - hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 2.5)); qyZDCC = qyZDCC - hrecenterevxSp->GetBinContent(hrecenterevxSp->FindBin(vx, 3.5)); + + if (recfinewitherror) { + qxZDCA = qxZDCA / meanxAvxerror; + qyZDCA = qyZDCA / meanyAvxerror; + qxZDCC = qxZDCC / meanxCvxerror; + qyZDCC = qyZDCC / meanyCvxerror; + } } if (useRecenterefineSp && hrecenterevySp) { + + float meanxAvy = hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 0.5)); + float meanyAvy = hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 1.5)); + float meanxCvy = hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 2.5)); + float meanyCvy = hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 3.5)); + + float meanxAvyerror = hrecenterevySp->GetBinError(hrecenterevySp->FindBin(vy, 0.5)); + float meanyAvyerror = hrecenterevySp->GetBinError(hrecenterevySp->FindBin(vy, 1.5)); + float meanxCvyerror = hrecenterevySp->GetBinError(hrecenterevySp->FindBin(vy, 2.5)); + float meanyCvyerror = hrecenterevySp->GetBinError(hrecenterevySp->FindBin(vy, 3.5)); + + if (rejbadeventvy) { + if ((TMath::Abs(meanxAvy) > 90000.0 || TMath::Abs(meanxCvy) > 90000.0 || TMath::Abs(meanyAvy) > 90000.0 || TMath::Abs(meanyCvy) > 90000.0) && (TMath::Abs(meanxAvyerror) > 9000.0 || TMath::Abs(meanxCvyerror) > 9000.0 || TMath::Abs(meanyAvyerror) > 9000.0 || TMath::Abs(meanyCvyerror) > 9000.0)) { + triggerevent = false; + spcalibrationtable(triggerevent, lastRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + return; + } + } + qxZDCA = qxZDCA - hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 0.5)); qyZDCA = qyZDCA - hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 1.5)); qxZDCC = qxZDCC - hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 2.5)); qyZDCC = qyZDCC - hrecenterevySp->GetBinContent(hrecenterevySp->FindBin(vy, 3.5)); + + if (recfinewitherror) { + qxZDCA = qxZDCA / meanxAvyerror; + qyZDCA = qyZDCA / meanyAvyerror; + qxZDCC = qxZDCC / meanxCvyerror; + qyZDCC = qyZDCC / meanyCvyerror; + } } if (useRecenterefineSp && hrecenterevzSp) { + + float meanxAvz = hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 0.5)); + float meanyAvz = hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 1.5)); + float meanxCvz = hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 2.5)); + float meanyCvz = hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 3.5)); + + float meanxAvzerror = hrecenterevzSp->GetBinError(hrecenterevzSp->FindBin(vz, 0.5)); + float meanyAvzerror = hrecenterevzSp->GetBinError(hrecenterevzSp->FindBin(vz, 1.5)); + float meanxCvzerror = hrecenterevzSp->GetBinError(hrecenterevzSp->FindBin(vz, 2.5)); + float meanyCvzerror = hrecenterevzSp->GetBinError(hrecenterevzSp->FindBin(vz, 3.5)); + + if (rejbadeventvz) { + if ((TMath::Abs(meanxAvz) > 90000.0 || TMath::Abs(meanxCvz) > 90000.0 || TMath::Abs(meanyAvz) > 90000.0 || TMath::Abs(meanyCvz) > 90000.0) && (TMath::Abs(meanxAvzerror) > 9000.0 || TMath::Abs(meanxCvzerror) > 9000.0 || TMath::Abs(meanyAvzerror) > 9000.0 || TMath::Abs(meanyCvzerror) > 9000.0)) { + triggerevent = false; + spcalibrationtable(triggerevent, lastRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + return; + } + } + qxZDCA = qxZDCA - hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 0.5)); qyZDCA = qyZDCA - hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 1.5)); qxZDCC = qxZDCC - hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 2.5)); qyZDCC = qyZDCC - hrecenterevzSp->GetBinContent(hrecenterevzSp->FindBin(vz, 3.5)); + + if (recfinewitherror) { + qxZDCA = qxZDCA / meanxAvzerror; + qyZDCA = qyZDCA / meanyAvzerror; + qxZDCC = qxZDCC / meanxCvzerror; + qyZDCC = qyZDCC / meanyCvzerror; + } } psiZDCC = 1.0 * TMath::ATan2(qyZDCC, qxZDCC); @@ -603,7 +721,7 @@ struct spvector { lastRunNumber = currentRunNumber; } - spcalibrationtable(triggerevent, lastRunNumber, centrality, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); + spcalibrationtable(triggerevent, lastRunNumber, centrality, vx, vy, vz, qxZDCA, qxZDCC, qyZDCA, qyZDCC, psiZDCC, psiZDCA); } }; diff --git a/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx b/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx index 09b2e909ec1..97967ac4fbd 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx @@ -299,7 +299,7 @@ struct nucleiSpectra { Filter trackFilter = nabs(aod::track::eta) < cfgCutEta && aod::track::tpcInnerParam > cfgCutTpcMom; using TrackCandidates = soa::Filtered>; - + using TrackCandidatesMC = soa::Filtered>; // Collisions with chentrality using CollWithCent = soa::Join::iterator; @@ -783,7 +783,7 @@ struct nucleiSpectra { PROCESS_SWITCH(nucleiSpectra, processDataFlowAlternative, "Data analysis with flow - alternative framework", false); Preslice tracksPerCollisions = aod::track::collisionId; - void processMC(soa::Join const& collisions, aod::McCollisions const& mcCollisions, soa::Join const& tracks, aod::McParticles const& particlesMC, aod::BCsWithTimestamps const&) + void processMC(soa::Join const& collisions, aod::McCollisions const& mcCollisions, TrackCandidatesMC const& tracks, aod::McParticles const& particlesMC, aod::BCsWithTimestamps const&) { nuclei::candidates.clear(); for (auto& c : mcCollisions) { diff --git a/PWGLF/TableProducer/Strangeness/v0qaanalysis.cxx b/PWGLF/TableProducer/Strangeness/v0qaanalysis.cxx index 875d0f1c671..85113fe1fa1 100644 --- a/PWGLF/TableProducer/Strangeness/v0qaanalysis.cxx +++ b/PWGLF/TableProducer/Strangeness/v0qaanalysis.cxx @@ -248,7 +248,7 @@ struct LfV0qaanalysis { v0.negTrack_as().tofNSigmaPr(), v0.posTrack_as().tofNSigmaPr(), v0.negTrack_as().tofNSigmaPi(), v0.posTrack_as().tofNSigmaPi(), v0.posTrack_as().hasTOF(), v0.negTrack_as().hasTOF(), lPDG, isPhysicalPrimary, - collision.centFT0M(), collision.centFV0A(), evFlag); + collision.centFT0M(), collision.centFV0A(), evFlag, v0.alpha(), v0.qtarm()); } } } @@ -360,7 +360,7 @@ struct LfV0qaanalysis { v0.negTrack_as().tofNSigmaPr(), v0.posTrack_as().tofNSigmaPr(), v0.negTrack_as().tofNSigmaPi(), v0.posTrack_as().tofNSigmaPi(), v0.posTrack_as().hasTOF(), v0.negTrack_as().hasTOF(), lPDG, isprimary, - mcCollision.centFT0M(), cent, evFlag); + mcCollision.centFT0M(), cent, evFlag, v0.alpha(), v0.qtarm()); } } diff --git a/PWGLF/Tasks/Strangeness/strangeness_in_jets.cxx b/PWGLF/Tasks/Strangeness/strangeness_in_jets.cxx index be54eb0aebf..07b8549e29e 100644 --- a/PWGLF/Tasks/Strangeness/strangeness_in_jets.cxx +++ b/PWGLF/Tasks/Strangeness/strangeness_in_jets.cxx @@ -168,6 +168,68 @@ struct strangeness_in_jets { registryData.add("OmegaPos_in_ue", "OmegaPos_in_ue", HistType::kTH3F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {200, 1.63, 1.71, "m_{p#piK} (GeV/#it{c}^{2})"}}); registryData.add("OmegaNeg_in_jet", "OmegaNeg_in_jet", HistType::kTH3F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {200, 1.63, 1.71, "m_{p#piK} (GeV/#it{c}^{2})"}}); registryData.add("OmegaNeg_in_ue", "OmegaNeg_in_ue", HistType::kTH3F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {200, 1.63, 1.71, "m_{p#piK} (GeV/#it{c}^{2})"}}); + + // Histograms for efficiency (generated) + registryMC.add("K0s_generated", "K0s_generated", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("Lambda_generated", "Lambda_generated", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("AntiLambda_generated", "AntiLambda_generated", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("XiPos_generated", "XiPos_generated", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("XiNeg_generated", "XiNeg_generated", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("OmegaPos_generated", "OmegaPos_generated", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("OmegaNeg_generated", "OmegaNeg_generated", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + + // Histograms for efficiency (reconstructed) + registryMC.add("K0s_reconstructed", "K0s_reconstructed", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("Lambda_reconstructed", "Lambda_reconstructed", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("AntiLambda_reconstructed", "AntiLambda_reconstructed", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("XiPos_reconstructed", "XiPos_reconstructed", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("XiNeg_reconstructed", "XiNeg_reconstructed", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("OmegaPos_reconstructed", "OmegaPos_reconstructed", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("OmegaNeg_reconstructed", "OmegaNeg_reconstructed", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + + // Histograms for secondary hadrons + registryMC.add("K0s_reconstructed_incl", "K0s_reconstructed_incl", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("Lambda_reconstructed_incl", "Lambda_reconstructed_incl", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("AntiLambda_reconstructed_incl", "AntiLambda_reconstructed_incl", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + + // Histograms for 2d reweighting (pion) + registryMC.add("Pion_eta_pt_jet", "Pion_eta_pt_jet", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Pion_eta_pt_ue", "Pion_eta_pt_ue", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Pion_eta_pt_pythia", "Pion_eta_pt_pythia", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + + // Histograms for 2d reweighting (K0s) + registryMC.add("K0s_eta_pt_jet", "K0s_eta_pt_jet", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("K0s_eta_pt_ue", "K0s_eta_pt_ue", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("K0s_eta_pt_pythia", "K0s_eta_pt_pythia", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + + // Histograms for 2d reweighting (Lambda) + registryMC.add("Lambda_eta_pt_jet", "Lambda_eta_pt_jet", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Lambda_eta_pt_ue", "Lambda_eta_pt_ue", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Lambda_eta_pt_pythia", "Lambda_eta_pt_pythia", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + + // Histograms for 2d reweighting (Xi) + registryMC.add("Xi_eta_pt_jet", "Xi_eta_pt_jet", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Xi_eta_pt_ue", "Xi_eta_pt_ue", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Xi_eta_pt_pythia", "Xi_eta_pt_pythia", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + + // Histograms for 2d reweighting (Omega) + registryMC.add("Omega_eta_pt_jet", "Omega_eta_pt_jet", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Omega_eta_pt_ue", "Omega_eta_pt_ue", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + registryMC.add("Omega_eta_pt_pythia", "Omega_eta_pt_pythia", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {18, -0.9, 0.9, "#eta"}}); + + // Histograms for efficiency (pions) + registryMC.add("pi_plus_gen", "pi_plus_gen", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("pi_plus_rec_tpc", "pi_plus_rec_tpc", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("pi_plus_rec_tof", "pi_plus_rec_tof", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("pi_minus_gen", "pi_minus_gen", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("pi_minus_rec_tpc", "pi_minus_rec_tpc", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("pi_minus_rec_tof", "pi_minus_rec_tof", HistType::kTH2F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); + + // MC Templates + registryMC.add("piplus_dcaxy_prim", "piplus_dcaxy_prim", HistType::kTH3F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + registryMC.add("piminus_dcaxy_prim", "piminus_dcaxy_prim", HistType::kTH3F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + registryMC.add("piplus_dcaxy_sec", "piplus_dcaxy_sec", HistType::kTH3F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + registryMC.add("piminus_dcaxy_sec", "piminus_dcaxy_sec", HistType::kTH3F, {multBinning, {100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); } template @@ -1076,6 +1138,500 @@ struct strangeness_in_jets { } } PROCESS_SWITCH(strangeness_in_jets, processData, "Process data", true); + + Preslice perCollisionV0 = o2::aod::v0data::collisionId; + Preslice perCollisionCasc = o2::aod::cascade::collisionId; + Preslice perMCCollision = o2::aod::mcparticle::mcCollisionId; + Preslice perCollisionTrk = o2::aod::track::collisionId; + + void processMCefficiency(SimCollisions const& collisions, MCTracks const& mcTracks, aod::V0Datas const& fullV0s, aod::CascDataExt const& Cascades, const aod::McParticles& mcParticles) + { + for (const auto& collision : collisions) { + registryMC.fill(HIST("number_of_events_mc"), 0.5); + if (!collision.sel8()) + continue; + + registryMC.fill(HIST("number_of_events_mc"), 1.5); + if (abs(collision.posZ()) > 10.0) + continue; + + registryMC.fill(HIST("number_of_events_mc"), 2.5); + float multiplicity = collision.centFT0M(); + + auto v0s_per_coll = fullV0s.sliceBy(perCollisionV0, collision.globalIndex()); + auto casc_per_coll = Cascades.sliceBy(perCollisionCasc, collision.globalIndex()); + auto mcParticles_per_coll = mcParticles.sliceBy(perMCCollision, collision.globalIndex()); + auto tracks_per_coll = mcTracks.sliceBy(perCollisionTrk, collision.globalIndex()); + + for (auto& v0 : v0s_per_coll) { + + const auto& pos = v0.posTrack_as(); + const auto& neg = v0.negTrack_as(); + if (!pos.has_mcParticle()) + continue; + if (!neg.has_mcParticle()) + continue; + + auto posParticle = pos.mcParticle_as(); + auto negParticle = neg.mcParticle_as(); + if (!posParticle.has_mothers()) + continue; + if (!negParticle.has_mothers()) + continue; + + int pdg_parent(0); + bool isPhysPrim = false; + for (auto& particleMotherOfNeg : negParticle.mothers_as()) { + for (auto& particleMotherOfPos : posParticle.mothers_as()) { + if (particleMotherOfNeg == particleMotherOfPos) { + pdg_parent = particleMotherOfNeg.pdgCode(); + isPhysPrim = particleMotherOfNeg.isPhysicalPrimary(); + } + } + } + if (pdg_parent == 0) + continue; + + if (passedK0ShortSelection(v0, pos, neg) && pdg_parent == 310) { + registryMC.fill(HIST("K0s_reconstructed_incl"), multiplicity, v0.pt()); + } + if (passedLambdaSelection(v0, pos, neg) && pdg_parent == 3122) { + registryMC.fill(HIST("Lambda_reconstructed_incl"), multiplicity, v0.pt()); + } + if (passedAntiLambdaSelection(v0, pos, neg) && pdg_parent == -3122) { + registryMC.fill(HIST("AntiLambda_reconstructed_incl"), multiplicity, v0.pt()); + } + if (!isPhysPrim) + continue; + + if (passedK0ShortSelection(v0, pos, neg) && pdg_parent == 310) { + registryMC.fill(HIST("K0s_reconstructed"), multiplicity, v0.pt()); + } + if (passedLambdaSelection(v0, pos, neg) && pdg_parent == 3122) { + registryMC.fill(HIST("Lambda_reconstructed"), multiplicity, v0.pt()); + } + if (passedAntiLambdaSelection(v0, pos, neg) && pdg_parent == -3122) { + registryMC.fill(HIST("AntiLambda_reconstructed"), multiplicity, v0.pt()); + } + } + + // Cascades + for (auto& casc : casc_per_coll) { + auto bach = casc.template bachelor_as(); + auto neg = casc.template negTrack_as(); + auto pos = casc.template posTrack_as(); + + if (!bach.has_mcParticle()) + continue; + if (!pos.has_mcParticle()) + continue; + if (!neg.has_mcParticle()) + continue; + + auto posParticle = pos.mcParticle_as(); + auto negParticle = neg.mcParticle_as(); + auto bachParticle = bach.mcParticle_as(); + if (!posParticle.has_mothers()) + continue; + if (!negParticle.has_mothers()) + continue; + if (!bachParticle.has_mothers()) + continue; + + int pdg_parent(0); + for (auto& particleMotherOfNeg : negParticle.mothers_as()) { + for (auto& particleMotherOfPos : posParticle.mothers_as()) { + for (auto& particleMotherOfBach : bachParticle.mothers_as()) { + if (particleMotherOfNeg != particleMotherOfPos) + continue; + if (abs(particleMotherOfNeg.pdgCode()) != 3122) + continue; + if (!particleMotherOfBach.isPhysicalPrimary()) + continue; + pdg_parent = particleMotherOfBach.pdgCode(); + } + } + } + if (pdg_parent == 0) + continue; + + // Xi+ + if (passedXiSelection(casc, pos, neg, bach, collision) && pdg_parent == -3312) { + registryMC.fill(HIST("XiPos_reconstructed"), multiplicity, casc.pt()); + } + // Xi- + if (passedXiSelection(casc, pos, neg, bach, collision) && pdg_parent == 3312) { + registryMC.fill(HIST("XiNeg_reconstructed"), multiplicity, casc.pt()); + } + // Omega+ + if (passedOmegaSelection(casc, pos, neg, bach, collision) && pdg_parent == -3334) { + registryMC.fill(HIST("OmegaPos_reconstructed"), multiplicity, casc.pt()); + } + // Omega- + if (passedOmegaSelection(casc, pos, neg, bach, collision) && pdg_parent == 3334) { + registryMC.fill(HIST("OmegaNeg_reconstructed"), multiplicity, casc.pt()); + } + } + + // Reconstructed Tracks + for (auto track : tracks_per_coll) { + + // Get MC Particle + if (!track.has_mcParticle()) + continue; + // Track Selection + if (!passedTrackSelectionForPions(track)) + continue; + + const auto particle = track.mcParticle(); + if (abs(particle.pdgCode()) != 211) + continue; + + if (particle.isPhysicalPrimary()) { + if (track.sign() > 0) + registryMC.fill(HIST("piplus_dcaxy_prim"), multiplicity, track.pt(), track.dcaXY()); + if (track.sign() < 0) + registryMC.fill(HIST("piminus_dcaxy_prim"), multiplicity, track.pt(), track.dcaXY()); + } + if (!particle.isPhysicalPrimary()) { + if (track.sign() > 0) + registryMC.fill(HIST("piplus_dcaxy_sec"), multiplicity, track.pt(), track.dcaXY()); + if (track.sign() < 0) + registryMC.fill(HIST("piminus_dcaxy_sec"), multiplicity, track.pt(), track.dcaXY()); + } + + if (TMath::Abs(track.dcaXY()) > dcaxyMax) + continue; + + if (track.tpcNSigmaPi() < nsigmaTPCmin || track.tpcNSigmaPi() > nsigmaTPCmax) + continue; + + if (!particle.isPhysicalPrimary()) + continue; + + if (track.sign() > 0) + registryMC.fill(HIST("pi_plus_rec_tpc"), multiplicity, track.pt()); + if (track.sign() < 0) + registryMC.fill(HIST("pi_minus_rec_tpc"), multiplicity, track.pt()); + + if (!track.hasTOF()) + continue; + if (track.tofNSigmaPi() < nsigmaTOFmin || track.tofNSigmaPi() > nsigmaTOFmax) + continue; + + if (track.sign() > 0) + registryMC.fill(HIST("pi_plus_rec_tof"), multiplicity, track.pt()); + if (track.sign() < 0) + registryMC.fill(HIST("pi_minus_rec_tof"), multiplicity, track.pt()); + } + + for (auto& mcParticle : mcParticles_per_coll) { + + if (mcParticle.eta() < etaMin || mcParticle.eta() > etaMax) + continue; + if (!mcParticle.isPhysicalPrimary()) + continue; + + // Pi+ + if (mcParticle.pdgCode() == 211) { + registryMC.fill(HIST("pi_plus_gen"), multiplicity, mcParticle.pt()); + } + // Pi- + if (mcParticle.pdgCode() == -211) { + registryMC.fill(HIST("pi_minus_gen"), multiplicity, mcParticle.pt()); + } + // K0s + if (mcParticle.pdgCode() == 310) { + registryMC.fill(HIST("K0s_generated"), multiplicity, mcParticle.pt()); + registryMC.fill(HIST("K0s_eta_pt_pythia"), mcParticle.pt(), mcParticle.eta()); + } + // Lambda + if (mcParticle.pdgCode() == 3122) { + registryMC.fill(HIST("Lambda_generated"), multiplicity, mcParticle.pt()); + registryMC.fill(HIST("Lambda_eta_pt_pythia"), mcParticle.pt(), mcParticle.eta()); + } + // AntiLambda + if (mcParticle.pdgCode() == -3122) { + registryMC.fill(HIST("AntiLambda_generated"), multiplicity, mcParticle.pt()); + registryMC.fill(HIST("Lambda_eta_pt_pythia"), mcParticle.pt(), mcParticle.eta()); + } + // Xi Pos + if (mcParticle.pdgCode() == -3312) { + registryMC.fill(HIST("XiPos_generated"), multiplicity, mcParticle.pt()); + registryMC.fill(HIST("Xi_eta_pt_pythia"), mcParticle.pt(), mcParticle.eta()); + } + // Xi Neg + if (mcParticle.pdgCode() == 3312) { + registryMC.fill(HIST("XiNeg_generated"), multiplicity, mcParticle.pt()); + registryMC.fill(HIST("Xi_eta_pt_pythia"), mcParticle.pt(), mcParticle.eta()); + } + // Omega Pos + if (mcParticle.pdgCode() == -3334) { + registryMC.fill(HIST("OmegaPos_generated"), multiplicity, mcParticle.pt()); + registryMC.fill(HIST("Omega_eta_pt_pythia"), mcParticle.pt(), mcParticle.eta()); + } + // Omega Neg + if (mcParticle.pdgCode() == 3334) { + registryMC.fill(HIST("OmegaNeg_generated"), multiplicity, mcParticle.pt()); + registryMC.fill(HIST("Omega_eta_pt_pythia"), mcParticle.pt(), mcParticle.eta()); + } + } + } + } + PROCESS_SWITCH(strangeness_in_jets, processMCefficiency, "Process MC Efficiency", false); + + void processGen(o2::aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles) + { + for (const auto& mccollision : mcCollisions) { + + registryMC.fill(HIST("number_of_events_mc"), 3.5); + + // Selection on z_{vertex} + if (abs(mccollision.posZ()) > 10) + continue; + registryMC.fill(HIST("number_of_events_mc"), 4.5); + + // MC Particles per Collision + auto mcParticles_per_coll = mcParticles.sliceBy(perMCCollision, mccollision.globalIndex()); + + // List of Tracks + std::vector trk; + + for (auto& particle : mcParticles_per_coll) { + + int pdg = abs(particle.pdgCode()); + + if (particle.isPhysicalPrimary() && pdg == 211) { + registryMC.fill(HIST("Pion_eta_pt_pythia"), particle.pt(), particle.eta()); + } + if (particle.isPhysicalPrimary() && pdg == 310) { + registryMC.fill(HIST("K0s_eta_pt_pythia"), particle.pt(), particle.eta()); + } + if (particle.isPhysicalPrimary() && pdg == 3122) { + registryMC.fill(HIST("Lambda_eta_pt_pythia"), particle.pt(), particle.eta()); + } + if (particle.isPhysicalPrimary() && pdg == 3312) { + registryMC.fill(HIST("Xi_eta_pt_pythia"), particle.pt(), particle.eta()); + } + if (particle.isPhysicalPrimary() && pdg == 3334) { + registryMC.fill(HIST("Omega_eta_pt_pythia"), particle.pt(), particle.eta()); + } + + // Select Primary Particles + double dx = particle.vx() - mccollision.posX(); + double dy = particle.vy() - mccollision.posY(); + double dz = particle.vz() - mccollision.posZ(); + double dcaxy = sqrt(dx * dx + dy * dy); + double dcaz = abs(dz); + if (dcaxy > (par0 + par1 / particle.pt())) + continue; + if (dcaz > (par0 + par1 / particle.pt())) + continue; + if (abs(particle.eta()) > 0.8) + continue; + if (particle.pt() < 0.15) + continue; + + // PDG Selection + if ((pdg != 11) && (pdg != 211) && (pdg != 321) && (pdg != 2212)) + continue; + + TVector3 momentum(particle.px(), particle.py(), particle.pz()); + trk.push_back(momentum); + } + + // Anti-kt Jet Finder + int n_particles_removed(0); + std::vector jet; + std::vector ue1; + std::vector ue2; + + do { + double dij_min(1e+06), diB_min(1e+06); + int i_min(0), j_min(0), iB_min(0); + for (int i = 0; i < static_cast(trk.size()); i++) { + if (trk[i].Mag() == 0) + continue; + double diB = 1.0 / (trk[i].Pt() * trk[i].Pt()); + if (diB < diB_min) { + diB_min = diB; + iB_min = i; + } + for (int j = (i + 1); j < static_cast(trk.size()); j++) { + if (trk[j].Mag() == 0) + continue; + double dij = calculate_dij(trk[i], trk[j], Rjet); + if (dij < dij_min) { + dij_min = dij; + i_min = i; + j_min = j; + } + } + } + if (dij_min < diB_min) { + trk[i_min] = trk[i_min] + trk[j_min]; + trk[j_min].SetXYZ(0, 0, 0); + n_particles_removed++; + } + if (dij_min > diB_min) { + jet.push_back(trk[iB_min]); + trk[iB_min].SetXYZ(0, 0, 0); + n_particles_removed++; + } + } while (n_particles_removed < static_cast(trk.size())); + + // Jet Selection + std::vector isSelected; + for (int i = 0; i < static_cast(jet.size()); i++) { + isSelected.push_back(0); + } + + int n_jets_selected(0); + for (int i = 0; i < static_cast(jet.size()); i++) { + + if ((abs(jet[i].Eta()) + Rjet) > etaMax) + continue; + + // Perpendicular cones + TVector3 ue_axis1(0, 0, 0); + TVector3 ue_axis2(0, 0, 0); + get_perpendicular_axis(jet[i], ue_axis1, +1); + get_perpendicular_axis(jet[i], ue_axis2, -1); + ue1.push_back(ue_axis1); + ue2.push_back(ue_axis2); + + double nPartJetPlusUE(0); + double ptJetPlusUE(0); + double ptJet(0); + double ptUE(0); + + for (auto& particle : mcParticles_per_coll) { + + // Select Primary Particles + double dx = particle.vx() - mccollision.posX(); + double dy = particle.vy() - mccollision.posY(); + double dz = particle.vz() - mccollision.posZ(); + double dcaxy = sqrt(dx * dx + dy * dy); + double dcaz = abs(dz); + + if (dcaxy > (par0 + par1 / particle.pt())) + continue; + if (dcaz > (par0 + par1 / particle.pt())) + continue; + if (abs(particle.eta()) > 0.8) + continue; + if (particle.pt() < 0.15) + continue; + + // PDG Selection + int pdg = abs(particle.pdgCode()); + if ((pdg != 11) && (pdg != 211) && (pdg != 321) && (pdg != 2212)) + continue; + + TVector3 sel_track(particle.px(), particle.py(), particle.pz()); + + double deltaEta_jet = sel_track.Eta() - jet[i].Eta(); + double deltaPhi_jet = GetDeltaPhi(sel_track.Phi(), jet[i].Phi()); + double deltaR_jet = sqrt(deltaEta_jet * deltaEta_jet + deltaPhi_jet * deltaPhi_jet); + double deltaEta_ue1 = sel_track.Eta() - ue_axis1.Eta(); + double deltaPhi_ue1 = GetDeltaPhi(sel_track.Phi(), ue_axis1.Phi()); + double deltaR_ue1 = sqrt(deltaEta_ue1 * deltaEta_ue1 + deltaPhi_ue1 * deltaPhi_ue1); + double deltaEta_ue2 = sel_track.Eta() - ue_axis2.Eta(); + double deltaPhi_ue2 = GetDeltaPhi(sel_track.Phi(), ue_axis2.Phi()); + double deltaR_ue2 = sqrt(deltaEta_ue2 * deltaEta_ue2 + deltaPhi_ue2 * deltaPhi_ue2); + + if (deltaR_jet < Rjet) { + nPartJetPlusUE++; + ptJetPlusUE = ptJetPlusUE + sel_track.Pt(); + } + if (deltaR_ue1 < Rjet) { + ptUE = ptUE + sel_track.Pt(); + } + if (deltaR_ue2 < Rjet) { + ptUE = ptUE + sel_track.Pt(); + } + } + ptJet = ptJetPlusUE - 0.5 * ptUE; + + if (ptJet < min_jet_pt) + continue; + if (nPartJetPlusUE < min_nPartInJet) + continue; + n_jets_selected++; + isSelected[i] = 1; + } + if (n_jets_selected == 0) + continue; + + for (int i = 0; i < static_cast(jet.size()); i++) { + + if (isSelected[i] == 0) + continue; + + // Generated Particles + for (auto& particle : mcParticles_per_coll) { + + if (!particle.isPhysicalPrimary()) + continue; + + TVector3 particle_dir(particle.px(), particle.py(), particle.pz()); + double deltaEta_jet = particle_dir.Eta() - jet[i].Eta(); + double deltaPhi_jet = GetDeltaPhi(particle_dir.Phi(), jet[i].Phi()); + double deltaR_jet = sqrt(deltaEta_jet * deltaEta_jet + deltaPhi_jet * deltaPhi_jet); + double deltaEta_ue1 = particle_dir.Eta() - ue1[i].Eta(); + double deltaPhi_ue1 = GetDeltaPhi(particle_dir.Phi(), ue1[i].Phi()); + double deltaR_ue1 = sqrt(deltaEta_ue1 * deltaEta_ue1 + deltaPhi_ue1 * deltaPhi_ue1); + double deltaEta_ue2 = particle_dir.Eta() - ue2[i].Eta(); + double deltaPhi_ue2 = GetDeltaPhi(particle_dir.Phi(), ue2[i].Phi()); + double deltaR_ue2 = sqrt(deltaEta_ue2 * deltaEta_ue2 + deltaPhi_ue2 * deltaPhi_ue2); + + int pdg = abs(particle.pdgCode()); + + if (pdg == 211) { + if (deltaR_jet < Rjet) { + registryMC.fill(HIST("Pion_eta_pt_jet"), particle.pt(), particle.eta()); + } + if (deltaR_ue1 < Rjet || deltaR_ue2 < Rjet) { + registryMC.fill(HIST("Pion_eta_pt_ue"), particle.pt(), particle.eta()); + } + } + if (pdg == 310) { + if (deltaR_jet < Rjet) { + registryMC.fill(HIST("K0s_eta_pt_jet"), particle.pt(), particle.eta()); + } + if (deltaR_ue1 < Rjet || deltaR_ue2 < Rjet) { + registryMC.fill(HIST("K0s_eta_pt_ue"), particle.pt(), particle.eta()); + } + } + if (pdg == 3122) { + if (deltaR_jet < Rjet) { + registryMC.fill(HIST("Lambda_eta_pt_jet"), particle.pt(), particle.eta()); + } + if (deltaR_ue1 < Rjet || deltaR_ue2 < Rjet) { + registryMC.fill(HIST("Lambda_eta_pt_ue"), particle.pt(), particle.eta()); + } + } + if (pdg == 3312) { + if (deltaR_jet < Rjet) { + registryMC.fill(HIST("Xi_eta_pt_jet"), particle.pt(), particle.eta()); + } + if (deltaR_ue1 < Rjet || deltaR_ue2 < Rjet) { + registryMC.fill(HIST("Xi_eta_pt_ue"), particle.pt(), particle.eta()); + } + } + if (pdg == 3334) { + if (deltaR_jet < Rjet) { + registryMC.fill(HIST("Omega_eta_pt_jet"), particle.pt(), particle.eta()); + } + if (deltaR_ue1 < Rjet || deltaR_ue2 < Rjet) { + registryMC.fill(HIST("Omega_eta_pt_ue"), particle.pt(), particle.eta()); + } + } + } + } + } + } + PROCESS_SWITCH(strangeness_in_jets, processGen, "Process generated MC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGLF/Tasks/Strangeness/v0postprocessing.cxx b/PWGLF/Tasks/Strangeness/v0postprocessing.cxx index 1844390a3da..d66f2814e1e 100644 --- a/PWGLF/Tasks/Strangeness/v0postprocessing.cxx +++ b/PWGLF/Tasks/Strangeness/v0postprocessing.cxx @@ -42,11 +42,14 @@ struct v0postprocessing { Configurable v0rejK0s{"v0rejK0s", 0.005, "V0 rej K0s"}; Configurable v0rejLambda{"v0rejLambda", 0.01, "V0 rej K0s"}; Configurable ntpcsigma{"ntpcsigma", 5, "N sigma TPC"}; + Configurable ntpcsigmaMC{"ntpcsigmaMC", 100, "N sigma TPC for MC"}; Configurable etadau{"etadau", 0.8, "Eta Daughters"}; Configurable isMC{"isMC", 1, "isMC"}; Configurable evSel{"evSel", 1, "evSel"}; Configurable hasTOF2Leg{"hasTOF2Leg", 0, "hasTOF2Leg"}; Configurable hasTOF1Leg{"hasTOF1Leg", 1, "hasTOF1Leg"}; + Configurable paramArmenterosCut{"paramArmenterosCut", 0.2, "parameter Armenteros Cut"}; + Configurable doArmenterosCut{"doArmenterosCut", 1, "do Armenteros Cut"}; HistogramRegistry registry{"registry"}; @@ -66,14 +69,16 @@ struct v0postprocessing { if (isMC) { registry.add("hMassK0Short_MC", ";M_{#pi^{+}#pi^{-}} [GeV/c^{2}]", {HistType::kTH1F, {{200, 0.4f, 0.6f}}}); - registry.add("hMassVsPtK0Short_MC", ";p_{T} [GeV/c];M_{#pi^{+}#pi^{-}} [GeV/c^{2}]", {HistType::kTH3F, {{250, 0.0f, 25.0f}, {200, 0.4f, 0.6f}, {100, 0.f, 100.f}}}); + registry.add("hMassVsPtK0Short_MC", ";p_{T} [GeV/c];M_{#pi^{+}#pi^{-}} [GeV/c^{2}]", {HistType::kTH3F, {{250, 0.0f, 25.0f}, {100, 0.f, 100.f}, {200, 0.4f, 0.6f}}}); registry.add("hMassLambda_MC", "hMassLambda", {HistType::kTH1F, {{200, 1.016f, 1.216f}}}); - registry.add("hMassVsPtLambda_MC", "hMassVsPtLambda", {HistType::kTH2F, {{100, 0.0f, 10.0f}, {200, 1.016f, 1.216f}}}); + registry.add("hMassVsPtLambda_MC", ";p_{T} [GeV/c];M_{p^{+}#pi^{-}} [GeV/c^{2}]", {HistType::kTH3F, {{250, 0.0f, 25.0f}, {100, 0.f, 100.f}, {200, 1.016f, 1.216f}}}); registry.add("hMassAntiLambda_MC", "hMassAntiLambda", {HistType::kTH1F, {{200, 1.016f, 1.216f}}}); - registry.add("hMassVsPtAntiLambda_MC", "hMassVsPtAntiLambda", {HistType::kTH2F, {{100, 0.0f, 10.0f}, {200, 1.016f, 1.216f}}}); + registry.add("hMassVsPtAntiLambda_MC", ";p_{T} [GeV/c];M_{p^{-}#pi^{+}} [GeV/c^{2}]", {HistType::kTH3F, {{250, 0.0f, 25.0f}, {100, 0.f, 100.f}, {200, 1.016f, 1.216f}}}); } // QA + registry.add("hArmenterosPodolanski", "hArmenterosPodolanski", {HistType::kTH2F, {{1000, -1.0f, 1.0f, "#alpha"}, {1000, 0.0f, 0.30f, "#it{Q}_{T}"}}}); + registry.add("hArmenterosPodolanski_Sel", "hArmenterosPodolanski_Sel", {HistType::kTH2F, {{1000, -1.0f, 1.0f, "#alpha"}, {1000, 0.0f, 0.30f, "#it{Q}_{T}"}}}); registry.add("hK0sV0Radius", "hK0sV0Radius", {HistType::kTH1D, {{200, 0.0f, 40.0f}}}); registry.add("hK0sCosPA", "hK0sCosPA", {HistType::kTH1F, {{100, 0.9f, 1.0f}}}); registry.add("hK0sV0DCANegToPV", "hK0sV0DCANegToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); @@ -85,38 +90,34 @@ struct v0postprocessing { registry.add("hK0sTPCNSigmaPosPi", "hK0sTPCNSigmaPosPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); registry.add("hK0sTPCNSigmaNegPi", "hK0sTPCNSigmaNegPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - /* registry.add("hLambdaV0Radius", "hLambdaV0Radius", {HistType::kTH1D, {{200, 0.0f, 40.0f}}}); - registry.add("hLambdaCosPA", "hLambdaCosPA", {HistType::kTH1F, {{100, 0.9f, 1.0f}}}); - registry.add("hLambdaV0DCANegToPV", "hLambdaV0DCANegToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); - registry.add("hLambdaV0DCAPosToPV", "hLambdaV0DCAPosToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); - registry.add("hLambdaV0DCAV0Daughters", "hLambdaV0DCAV0Daughters", {HistType::kTH1F, {{55, 0.0f, 2.20f}}}); - registry.add("hLambdaCtau", "hLambdaCtau", {HistType::kTH1F, {{100, 0.0f, 50.0f}}}); - registry.add("hLambdaEtaDau", "hLambdaEtaDau", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); - registry.add("hLambdaRap", "hLambdaRap", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); - registry.add("hLambdaTPCNSigmaPosPi", "hLambdaTPCNSigmaPosPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("hLambdaTPCNSigmaNegPi", "hLambdaTPCNSigmaNegPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("hLambdaTPCNSigmaNegPr", "hLambdaTPCNSigmaNegPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("hLambdaTPCNSigmaPosPr", "hLambdaTPCNSigmaPosPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - - registry.add("hAntiLambdaV0Radius", "hAntiLambdaV0Radius", {HistType::kTH1D, {{200, 0.0f, 40.0f}}}); - registry.add("hAntiLambdaCosPA", "hAntiLambdaCosPA", {HistType::kTH1F, {{100, 0.9f, 1.0f}}}); - registry.add("hAntiLambdaV0DCANegToPV", "hAntiLambdaV0DCANegToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); - registry.add("hAntiLambdaV0DCAPosToPV", "hAntiLambdaV0DCAPosToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); - registry.add("hAntiLambdaV0DCAV0Daughters", "hAntiLambdaV0DCAV0Daughters", {HistType::kTH1F, {{55, 0.0f, 2.20f}}}); - registry.add("hAntiLambdaCtau", "hAntiLambdaCtau", {HistType::kTH1F, {{100, 0.0f, 50.0f}}}); - registry.add("hAntiLambdaEtaDau", "hAntiLambdaEtaDau", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); - registry.add("hAntiLambdaRap", "hAntiLambdaRap", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); - registry.add("hAntiLambdaTPCNSigmaPosPi", "hAntiLambdaTPCNSigmaPosPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("hAntiLambdaTPCNSigmaNegPi", "hAntiLambdaTPCNSigmaNegPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("hAntiLambdaTPCNSigmaNegPr", "hAntiLambdaTPCNSigmaNegPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("hAntiLambdaTPCNSigmaPosPr", "hAntiLambdaTPCNSigmaPosPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - - registry.add("TPCNSigmaPosPr", "TPCNSigmaPosPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("TPCNSigmaNegPr", "TPCNSigmaNegPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("TOFNSigmaPosPi", "TOFNSigmaPosPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("TOFNSigmaNegPi", "TOFNSigmaNegPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("TOFNSigmaPosPr", "TOFNSigmaPosPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); - registry.add("TOFNSigmaNegPr", "TOFNSigmaNegPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); */ + registry.add("hLambdaV0Radius", "hLambdaV0Radius", {HistType::kTH1D, {{200, 0.0f, 40.0f}}}); + registry.add("hLambdaCosPA", "hLambdaCosPA", {HistType::kTH1F, {{100, 0.9f, 1.0f}}}); + registry.add("hLambdaV0DCANegToPV", "hLambdaV0DCANegToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); + registry.add("hLambdaV0DCAPosToPV", "hLambdaV0DCAPosToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); + registry.add("hLambdaV0DCAV0Daughters", "hLambdaV0DCAV0Daughters", {HistType::kTH1F, {{55, 0.0f, 2.20f}}}); + registry.add("hLambdaCtau", "hLambdaCtau", {HistType::kTH1F, {{100, 0.0f, 50.0f}}}); + registry.add("hLambdaEtaDau", "hLambdaEtaDau", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); + registry.add("hLambdaRap", "hLambdaRap", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); + registry.add("hLambdaTPCNSigmaNegPi", "hLambdaTPCNSigmaNegPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + registry.add("hLambdaTPCNSigmaPosPr", "hLambdaTPCNSigmaPosPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + + registry.add("hAntiLambdaV0Radius", "hAntiLambdaV0Radius", {HistType::kTH1D, {{200, 0.0f, 40.0f}}}); + registry.add("hAntiLambdaCosPA", "hAntiLambdaCosPA", {HistType::kTH1F, {{100, 0.9f, 1.0f}}}); + registry.add("hAntiLambdaV0DCANegToPV", "hAntiLambdaV0DCANegToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); + registry.add("hAntiLambdaV0DCAPosToPV", "hAntiLambdaV0DCAPosToPV", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}); + registry.add("hAntiLambdaV0DCAV0Daughters", "hAntiLambdaV0DCAV0Daughters", {HistType::kTH1F, {{55, 0.0f, 2.20f}}}); + registry.add("hAntiLambdaCtau", "hAntiLambdaCtau", {HistType::kTH1F, {{100, 0.0f, 50.0f}}}); + registry.add("hAntiLambdaEtaDau", "hAntiLambdaEtaDau", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); + registry.add("hAntiLambdaRap", "hAntiLambdaRap", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}); + registry.add("hAntiLambdaTPCNSigmaPosPi", "hAntiLambdaTPCNSigmaPosPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + registry.add("hAntiLambdaTPCNSigmaNegPr", "hAntiLambdaTPCNSigmaNegPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + + /*registry.add("TPCNSigmaPosPr", "TPCNSigmaPosPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + registry.add("TPCNSigmaNegPr", "TPCNSigmaNegPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + registry.add("TOFNSigmaPosPi", "TOFNSigmaPosPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + registry.add("TOFNSigmaNegPi", "TOFNSigmaNegPi", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + registry.add("TOFNSigmaPosPr", "TOFNSigmaPosPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); + registry.add("TOFNSigmaNegPr", "TOFNSigmaNegPr", {HistType::kTH1F, {{100, -10.0f, 10.0f}}}); */ } void process(aod::MyV0Candidates const& myv0s) @@ -136,10 +137,6 @@ struct v0postprocessing { continue; if (candidate.v0dcav0daughters() > dcav0dau) continue; - if (TMath::Abs(candidate.ntpcsigmanegpi()) > ntpcsigma) - continue; - if (TMath::Abs(candidate.ntpcsigmapospi()) > ntpcsigma) - continue; if (evSel && candidate.evflag() < 1) continue; if (hasTOF1Leg && !candidate.poshastof() && !candidate.neghastof()) @@ -152,57 +149,88 @@ struct v0postprocessing { TMath::Abs(candidate.rapk0short()) < rap && candidate.ctauk0short() < ctauK0s && TMath::Abs(candidate.massk0short() - o2::constants::physics::MassK0Short) < 0.075 && - TMath::Abs(candidate.masslambda() - o2::constants::physics::MassLambda0) > v0rejK0s) { - - registry.fill(HIST("hMassK0Short"), candidate.massk0short()); - registry.fill(HIST("hMassVsPtK0Short"), candidate.v0pt(), candidate.massk0short()); - registry.fill(HIST("hMassVsPtK0ShortVsCentFT0M"), candidate.v0pt(), candidate.multft0m(), candidate.massk0short()); - registry.fill(HIST("hMassVsPtK0ShortVsCentFV0A"), candidate.v0pt(), candidate.multfv0a(), candidate.massk0short()); - - // QA - if (!isMC) { - registry.fill(HIST("hK0sV0Radius"), candidate.v0radius()); - registry.fill(HIST("hK0sCosPA"), candidate.v0cospa()); - registry.fill(HIST("hK0sV0DCANegToPV"), candidate.v0dcanegtopv()); - registry.fill(HIST("hK0sV0DCAPosToPV"), candidate.v0dcapostopv()); - registry.fill(HIST("hK0sV0DCAV0Daughters"), candidate.v0dcav0daughters()); - registry.fill(HIST("hK0sCtau"), candidate.ctauk0short()); - registry.fill(HIST("hK0sEtaDau"), candidate.v0poseta()); - registry.fill(HIST("hK0sRap"), candidate.rapk0short()); - registry.fill(HIST("hK0sTPCNSigmaPosPi"), candidate.ntpcsigmapospi()); - registry.fill(HIST("hK0sTPCNSigmaNegPi"), candidate.ntpcsigmanegpi()); + TMath::Abs(candidate.masslambda() - o2::constants::physics::MassLambda0) > v0rejK0s && + TMath::Abs(candidate.ntpcsigmanegpi()) <= ntpcsigma && + TMath::Abs(candidate.ntpcsigmapospi()) <= ntpcsigma) { + + registry.fill(HIST("hArmenterosPodolanski"), candidate.alpha(), candidate.qtarm()); + + if (doArmenterosCut && candidate.qtarm() > (paramArmenterosCut * TMath::Abs(candidate.alpha()))) { + registry.fill(HIST("hArmenterosPodolanski_Sel"), candidate.alpha(), candidate.qtarm()); + registry.fill(HIST("hMassK0Short"), candidate.massk0short()); + registry.fill(HIST("hMassVsPtK0Short"), candidate.v0pt(), candidate.massk0short()); + registry.fill(HIST("hMassVsPtK0ShortVsCentFT0M"), candidate.v0pt(), candidate.multft0m(), candidate.massk0short()); + registry.fill(HIST("hMassVsPtK0ShortVsCentFV0A"), candidate.v0pt(), candidate.multfv0a(), candidate.massk0short()); + + // QA + if (!isMC) { + registry.fill(HIST("hK0sV0Radius"), candidate.v0radius()); + registry.fill(HIST("hK0sCosPA"), candidate.v0cospa()); + registry.fill(HIST("hK0sV0DCANegToPV"), candidate.v0dcanegtopv()); + registry.fill(HIST("hK0sV0DCAPosToPV"), candidate.v0dcapostopv()); + registry.fill(HIST("hK0sV0DCAV0Daughters"), candidate.v0dcav0daughters()); + registry.fill(HIST("hK0sCtau"), candidate.ctauk0short()); + registry.fill(HIST("hK0sEtaDau"), candidate.v0poseta()); + registry.fill(HIST("hK0sRap"), candidate.rapk0short()); + registry.fill(HIST("hK0sTPCNSigmaPosPi"), candidate.ntpcsigmapospi()); + registry.fill(HIST("hK0sTPCNSigmaNegPi"), candidate.ntpcsigmanegpi()); + } } } // Lambda analysis - /* if (candidate.v0cospa() > cospaLambda && - TMath::Abs(candidate.raplambda()) < rap && - candidate.ctaulambda() < ctauK0s && - TMath::Abs(candidate.masslambda() - o2::constants::physics::MassLambda0) < 0.075 && - TMath::Abs(candidate.massk0short() - o2::constants::physics::MassK0Short) > v0rejLambda) { - - registry.fill(HIST("hMassLambda"), candidate.masslambda()); - registry.fill(HIST("hMassVsPtLambda"), candidate.v0pt(), candidate.masslambda()); - registry.fill(HIST("hMassVsPtLambdaVsCentFT0M"), candidate.v0pt(), candidate.multft0m(), candidate.masslambda()); - registry.fill(HIST("hMassVsPtLambdaVsCentFV0A"), candidate.v0pt(), candidate.multfv0a(), candidate.masslambda()); - - // QA - if (!isMC) { - registry.fill(HIST("hLambdaV0Radius"), candidate.v0radius()); - registry.fill(HIST("hLambdaCosPA"), candidate.v0cospa()); - registry.fill(HIST("hLambdaV0DCANegToPV"), candidate.v0dcanegtopv()); - registry.fill(HIST("hLambdaV0DCAPosToPV"), candidate.v0dcapostopv()); - registry.fill(HIST("hLambdaV0DCAV0Daughters"), candidate.v0dcav0daughters()); - registry.fill(HIST("hLambdaCtau"), candidate.ctaulambda()); - registry.fill(HIST("hLambdaEtaDau"), candidate.v0poseta()); - registry.fill(HIST("hLambdaRap"), candidate.raplambda()); - registry.fill(HIST("hLambdaTPCNSigmaPosPi"), candidate.ntpcsigmapospi()); - registry.fill(HIST("hLambdaTPCNSigmaPosPr"), candidate.ntpcsigmapospr()); - registry.fill(HIST("hLambdaTPCNSigmaNegPi"), candidate.ntpcsigmanegpi()); - registry.fill(HIST("hLambdaTPCNSigmaNegPr"), candidate.ntpcsigmanegpr()); - } - } - */ + if (candidate.v0cospa() > cospaLambda && + TMath::Abs(candidate.raplambda()) < rap && + TMath::Abs(candidate.massk0short() - o2::constants::physics::MassK0Short) > v0rejLambda) { + + // Lambda + if (TMath::Abs(candidate.ntpcsigmanegpi()) <= ntpcsigma && TMath::Abs(candidate.ntpcsigmapospr()) <= ntpcsigma && + candidate.ctaulambda() < ctauLambda && + TMath::Abs(candidate.masslambda() - o2::constants::physics::MassLambda0) < 0.075) { + + registry.fill(HIST("hMassLambda"), candidate.masslambda()); + registry.fill(HIST("hMassVsPtLambda"), candidate.v0pt(), candidate.masslambda()); + registry.fill(HIST("hMassVsPtLambdaVsCentFT0M"), candidate.v0pt(), candidate.multft0m(), candidate.masslambda()); + + // QA + if (!isMC) { + registry.fill(HIST("hLambdaV0Radius"), candidate.v0radius()); + registry.fill(HIST("hLambdaCosPA"), candidate.v0cospa()); + registry.fill(HIST("hLambdaV0DCANegToPV"), candidate.v0dcanegtopv()); + registry.fill(HIST("hLambdaV0DCAPosToPV"), candidate.v0dcapostopv()); + registry.fill(HIST("hLambdaV0DCAV0Daughters"), candidate.v0dcav0daughters()); + registry.fill(HIST("hLambdaCtau"), candidate.ctaulambda()); + registry.fill(HIST("hLambdaEtaDau"), candidate.v0poseta()); + registry.fill(HIST("hLambdaRap"), candidate.raplambda()); + registry.fill(HIST("hLambdaTPCNSigmaPosPr"), candidate.ntpcsigmapospr()); + registry.fill(HIST("hLambdaTPCNSigmaNegPi"), candidate.ntpcsigmanegpi()); + } + } + // AntiLambda + if (TMath::Abs(candidate.ntpcsigmanegpr()) <= ntpcsigma && TMath::Abs(candidate.ntpcsigmapospi()) <= ntpcsigma && + candidate.ctauantilambda() < ctauLambda && + TMath::Abs(candidate.massantilambda() - o2::constants::physics::MassLambda0) < 0.075) { + + registry.fill(HIST("hMassAntiLambda"), candidate.massantilambda()); + registry.fill(HIST("hMassVsPtAntiLambda"), candidate.v0pt(), candidate.massantilambda()); + registry.fill(HIST("hMassVsPtAntiLambdaVsCentFT0M"), candidate.v0pt(), candidate.multft0m(), candidate.massantilambda()); + + // QA + if (!isMC) { + registry.fill(HIST("hAntiLambdaV0Radius"), candidate.v0radius()); + registry.fill(HIST("hAntiLambdaCosPA"), candidate.v0cospa()); + registry.fill(HIST("hAntiLambdaV0DCANegToPV"), candidate.v0dcanegtopv()); + registry.fill(HIST("hAntiLambdaV0DCAPosToPV"), candidate.v0dcapostopv()); + registry.fill(HIST("hAntiLambdaV0DCAV0Daughters"), candidate.v0dcav0daughters()); + registry.fill(HIST("hAntiLambdaCtau"), candidate.ctauantilambda()); + registry.fill(HIST("hAntiLambdaEtaDau"), candidate.v0poseta()); + registry.fill(HIST("hAntiLambdaRap"), candidate.raplambda()); + registry.fill(HIST("hAntiLambdaTPCNSigmaNegPr"), candidate.ntpcsigmanegpr()); + registry.fill(HIST("hAntiLambdaTPCNSigmaPosPi"), candidate.ntpcsigmapospi()); + } + } + } + if (isMC) { if (candidate.isphysprimary() == 0) @@ -214,23 +242,77 @@ struct v0postprocessing { candidate.ctauk0short() < ctauK0s && TMath::Abs(candidate.massk0short() - o2::constants::physics::MassK0Short) < 0.075 && TMath::Abs(candidate.masslambda() - o2::constants::physics::MassLambda0) > v0rejK0s && + TMath::Abs(candidate.ntpcsigmanegpi()) <= ntpcsigmaMC && + TMath::Abs(candidate.ntpcsigmapospi()) <= ntpcsigmaMC && (candidate.pdgcode() == 310)) { - registry.fill(HIST("hMassK0Short_MC"), candidate.massk0short()); - registry.fill(HIST("hMassVsPtK0Short_MC"), candidate.v0pt(), candidate.multft0m(), candidate.massk0short()); - - registry.fill(HIST("hK0sV0Radius"), candidate.v0radius()); - registry.fill(HIST("hK0sCosPA"), candidate.v0cospa()); - registry.fill(HIST("hK0sV0DCANegToPV"), candidate.v0dcanegtopv()); - registry.fill(HIST("hK0sV0DCAPosToPV"), candidate.v0dcapostopv()); - registry.fill(HIST("hK0sV0DCAV0Daughters"), candidate.v0dcav0daughters()); - registry.fill(HIST("hK0sCtau"), candidate.ctauk0short()); - registry.fill(HIST("hK0sEtaDau"), candidate.v0poseta()); - registry.fill(HIST("hK0sRap"), candidate.rapk0short()); - registry.fill(HIST("hK0sTPCNSigmaPosPi"), candidate.ntpcsigmapospi()); - registry.fill(HIST("hK0sTPCNSigmaNegPi"), candidate.ntpcsigmanegpi()); - } - } + registry.fill(HIST("hArmenterosPodolanski"), candidate.alpha(), candidate.qtarm()); + + if (doArmenterosCut && candidate.qtarm() > (paramArmenterosCut * TMath::Abs(candidate.alpha()))) { + registry.fill(HIST("hArmenterosPodolanski_Sel"), candidate.alpha(), candidate.qtarm()); + registry.fill(HIST("hMassK0Short_MC"), candidate.massk0short()); + registry.fill(HIST("hMassVsPtK0Short_MC"), candidate.v0pt(), candidate.multft0m(), candidate.massk0short()); + + registry.fill(HIST("hK0sV0Radius"), candidate.v0radius()); + registry.fill(HIST("hK0sCosPA"), candidate.v0cospa()); + registry.fill(HIST("hK0sV0DCANegToPV"), candidate.v0dcanegtopv()); + registry.fill(HIST("hK0sV0DCAPosToPV"), candidate.v0dcapostopv()); + registry.fill(HIST("hK0sV0DCAV0Daughters"), candidate.v0dcav0daughters()); + registry.fill(HIST("hK0sCtau"), candidate.ctauk0short()); + registry.fill(HIST("hK0sEtaDau"), candidate.v0poseta()); + registry.fill(HIST("hK0sRap"), candidate.rapk0short()); + registry.fill(HIST("hK0sTPCNSigmaPosPi"), candidate.ntpcsigmapospi()); + registry.fill(HIST("hK0sTPCNSigmaNegPi"), candidate.ntpcsigmanegpi()); + } + } // k0 + + // Lambda analysis + if (candidate.v0cospa() > cospaLambda && + TMath::Abs(candidate.raplambda()) < rap && + TMath::Abs(candidate.massk0short() - o2::constants::physics::MassK0Short) > v0rejLambda) { + + // Lambda + if (TMath::Abs(candidate.ntpcsigmanegpi()) <= ntpcsigmaMC && TMath::Abs(candidate.ntpcsigmapospr()) <= ntpcsigmaMC && + candidate.ctaulambda() < ctauLambda && + TMath::Abs(candidate.masslambda() - o2::constants::physics::MassLambda0) < 0.075 && + candidate.pdgcode() == 3122) { + + registry.fill(HIST("hMassLambda_MC"), candidate.masslambda()); + registry.fill(HIST("hMassVsPtLambda_MC"), candidate.v0pt(), candidate.multft0m(), candidate.masslambda()); + + registry.fill(HIST("hLambdaV0Radius"), candidate.v0radius()); + registry.fill(HIST("hLambdaCosPA"), candidate.v0cospa()); + registry.fill(HIST("hLambdaV0DCANegToPV"), candidate.v0dcanegtopv()); + registry.fill(HIST("hLambdaV0DCAPosToPV"), candidate.v0dcapostopv()); + registry.fill(HIST("hLambdaV0DCAV0Daughters"), candidate.v0dcav0daughters()); + registry.fill(HIST("hLambdaCtau"), candidate.ctaulambda()); + registry.fill(HIST("hLambdaEtaDau"), candidate.v0poseta()); + registry.fill(HIST("hLambdaRap"), candidate.raplambda()); + registry.fill(HIST("hLambdaTPCNSigmaPosPr"), candidate.ntpcsigmapospr()); + registry.fill(HIST("hLambdaTPCNSigmaNegPi"), candidate.ntpcsigmanegpi()); + } + // AntiLambda + if (TMath::Abs(candidate.ntpcsigmanegpr()) <= ntpcsigmaMC && TMath::Abs(candidate.ntpcsigmapospi()) <= ntpcsigmaMC && + candidate.ctauantilambda() < ctauLambda && + TMath::Abs(candidate.massantilambda() - o2::constants::physics::MassLambda0) < 0.075 && + candidate.pdgcode() == -3122) { + + registry.fill(HIST("hMassAntiLambda_MC"), candidate.massantilambda()); + registry.fill(HIST("hMassVsPtAntiLambda_MC"), candidate.v0pt(), candidate.multft0m(), candidate.massantilambda()); + + registry.fill(HIST("hAntiLambdaV0Radius"), candidate.v0radius()); + registry.fill(HIST("hAntiLambdaCosPA"), candidate.v0cospa()); + registry.fill(HIST("hAntiLambdaV0DCANegToPV"), candidate.v0dcanegtopv()); + registry.fill(HIST("hAntiLambdaV0DCAPosToPV"), candidate.v0dcapostopv()); + registry.fill(HIST("hAntiLambdaV0DCAV0Daughters"), candidate.v0dcav0daughters()); + registry.fill(HIST("hAntiLambdaCtau"), candidate.ctauantilambda()); + registry.fill(HIST("hAntiLambdaEtaDau"), candidate.v0poseta()); + registry.fill(HIST("hAntiLambdaRap"), candidate.raplambda()); + registry.fill(HIST("hAntiLambdaTPCNSigmaPosPi"), candidate.ntpcsigmapospi()); + registry.fill(HIST("hAntiLambdaTPCNSigmaNegPr"), candidate.ntpcsigmanegpr()); + } + } // lambda + } // is MC } } };