From 54ff4c011ee8c50412282741b2b04d11fd619bfb Mon Sep 17 00:00:00 2001 From: yuanzhe <90246048+wang-yuanzhe@users.noreply.github.com> Date: Thu, 12 Sep 2024 17:02:07 +0200 Subject: [PATCH 01/14] PWGLF: Make code compatible with additional particles in PID.h (#7656) * Make code compatible with additional particles in PID.h * Fix a small bug * Disable static_assert check to make compilation work * Fix overflow of the array --- PWGLF/Utils/mcParticle.h | 103 +++++++++++++++++++++------------------ 1 file changed, 56 insertions(+), 47 deletions(-) diff --git a/PWGLF/Utils/mcParticle.h b/PWGLF/Utils/mcParticle.h index dc0f476cc8a..e85d98ca5f0 100644 --- a/PWGLF/Utils/mcParticle.h +++ b/PWGLF/Utils/mcParticle.h @@ -53,6 +53,8 @@ class PIDExtended static constexpr ID Hyperhydrog4 = 14; static constexpr ID XiMinus = 15; static constexpr ID OmegaMinus = 16; + static constexpr ID HyperHelium4 = 17; + static constexpr ID HyperHelium5 = 18; static_assert(Electron == o2::track::PID::Electron, "PID::Electron mismatch"); static_assert(Muon == o2::track::PID::Muon, "PID::Muon mismatch"); @@ -71,59 +73,60 @@ class PIDExtended static_assert(Hyperhydrog4 == o2::track::PID::Hyperhydrog4, "PID::Hyperhydrog4 mismatch"); static_assert(XiMinus == o2::track::PID::XiMinus, "PID::XiMinus mismatch"); static_assert(OmegaMinus == o2::track::PID::OmegaMinus, "PID::OmegaMinus mismatch"); + // static_assert(HyperHelium4 == o2::track::PID::HyperHelium4, "PID::HyperHelium4 mismatch"); + // static_assert(HyperHelium5 == o2::track::PID::HyperHelium5, "PID::HyperHelium5 mismatch"); static constexpr ID PIDCountsUntilAl = 9; // Number of indices defined in PID.h equivalent to o2::track::PID::NIDs - static_assert(PIDCountsUntilAl == o2::track::PID::NIDs, "PID::NIDs mismatch"); + // static_assert(PIDCountsUntilAl == o2::track::PID::NIDs, "PID::NIDs mismatch"); - static constexpr ID PIDCounts = 17; // Number of indices defined in PID.h - static_assert(PIDCounts == o2::track::PID::NIDsTot, "PID::NIDsTot mismatch"); + static constexpr ID PIDCounts = 19; // Number of indices defined in PID.h + // static_assert(PIDCounts == o2::track::PID::NIDsTot, "PID::NIDsTot mismatch"); // Define an array of IDs static constexpr std::array mIDsUntilAl = {Electron, Muon, Pion, Kaon, Proton, Deuteron, Triton, Helium3, Alpha}; - static constexpr std::array mIDs = {Electron, Muon, Pion, Kaon, Proton, Deuteron, Triton, Helium3, Alpha, PI0, Photon, K0, Lambda, HyperTriton, Hyperhydrog4, XiMinus, OmegaMinus}; + static constexpr std::array mIDs = {Electron, Muon, Pion, Kaon, Proton, Deuteron, Triton, Helium3, Alpha, PI0, Photon, K0, Lambda, HyperTriton, Hyperhydrog4, XiMinus, OmegaMinus, HyperHelium4, HyperHelium5}; // Define the antiparticles - static constexpr ID Positron = 17; - static constexpr ID MuonPlus = 18; - static constexpr ID PionMinus = 19; - static constexpr ID KaonMinus = 20; - static constexpr ID AntiProton = 21; - static constexpr ID AntiDeuteron = 22; - static constexpr ID AntiTriton = 23; - static constexpr ID AntiHelium3 = 24; - static constexpr ID AntiAlpha = 25; - static constexpr ID AntiLambda = 26; - static constexpr ID AntiHyperTriton = 27; - static constexpr ID AntiHyperhydrog4 = 28; - static constexpr ID XiPlus = 29; - static constexpr ID OmegaPlus = 30; + static constexpr ID Positron = PIDCounts; + static constexpr ID MuonPlus = PIDCounts + 1; + static constexpr ID PionMinus = PIDCounts + 2; + static constexpr ID KaonMinus = PIDCounts + 3; + static constexpr ID AntiProton = PIDCounts + 4; + static constexpr ID AntiDeuteron = PIDCounts + 5; + static constexpr ID AntiTriton = PIDCounts + 6; + static constexpr ID AntiHelium3 = PIDCounts + 7; + static constexpr ID AntiAlpha = PIDCounts + 8; + static constexpr ID AntiLambda = PIDCounts + 9; + static constexpr ID AntiHyperTriton = PIDCounts + 10; + static constexpr ID AntiHyperhydrog4 = PIDCounts + 11; + static constexpr ID XiPlus = PIDCounts + 12; + static constexpr ID OmegaPlus = PIDCounts + 13; + static constexpr ID AntiHyperHelium4 = PIDCounts + 14; + static constexpr ID AntiHyperHelium5 = PIDCounts + 15; - static constexpr ID Neutron = 31; - static constexpr ID AntiNeutron = 32; - static constexpr ID HyperHelium4 = 33; - static constexpr ID AntiHyperHelium4 = 34; - static constexpr ID Phi = 35; - - static constexpr ID BZero = 36; - static constexpr ID BPlus = 37; - static constexpr ID BS = 38; - static constexpr ID D0 = 39; - static constexpr ID DPlus = 40; - static constexpr ID DS = 41; - static constexpr ID DStar = 42; - static constexpr ID ChiC1 = 43; - static constexpr ID JPsi = 44; - static constexpr ID LambdaB0 = 45; - static constexpr ID LambdaCPlus = 46; - static constexpr ID OmegaC0 = 47; - static constexpr ID SigmaC0 = 48; - static constexpr ID SigmaCPlusPlus = 49; - static constexpr ID X3872 = 50; - static constexpr ID Xi0 = 51; - static constexpr ID XiB0 = 52; - static constexpr ID XiCCPlusPlus = 53; - static constexpr ID XiCPlus = 54; - static constexpr ID XiC0 = 55; - static constexpr ID NIDsTot = 56; + static constexpr ID Neutron = PIDCounts + 16; + static constexpr ID AntiNeutron = PIDCounts + 17; + static constexpr ID Phi = PIDCounts + 18; + static constexpr ID BZero = PIDCounts + 19; + static constexpr ID BPlus = PIDCounts + 20; + static constexpr ID BS = PIDCounts + 21; + static constexpr ID D0 = PIDCounts + 22; + static constexpr ID DPlus = PIDCounts + 23; + static constexpr ID DS = PIDCounts + 24; + static constexpr ID DStar = PIDCounts + 25; + static constexpr ID ChiC1 = PIDCounts + 26; + static constexpr ID JPsi = PIDCounts + 27; + static constexpr ID LambdaB0 = PIDCounts + 28; + static constexpr ID LambdaCPlus = PIDCounts + 29; + static constexpr ID OmegaC0 = PIDCounts + 30; + static constexpr ID SigmaC0 = PIDCounts + 31; + static constexpr ID SigmaCPlusPlus = PIDCounts + 32; + static constexpr ID X3872 = PIDCounts + 33; + static constexpr ID Xi0 = PIDCounts + 34; + static constexpr ID XiB0 = PIDCounts + 35; + static constexpr ID XiCCPlusPlus = PIDCounts + 36; + static constexpr ID XiCPlus = PIDCounts + 37; + static constexpr ID XiC0 = PIDCounts + 38; + static constexpr ID NIDsTot = PIDCounts + 39; static constexpr const char* sNames[NIDsTot + 1] = { o2::track::pid_constants::sNames[Electron], // Electron @@ -143,6 +146,8 @@ class PIDExtended o2::track::pid_constants::sNames[Hyperhydrog4], // Hyperhydrog4 o2::track::pid_constants::sNames[XiMinus], // XiMinus o2::track::pid_constants::sNames[OmegaMinus], // OmegaMinus + "HyperHelium4", // HyperHelium4 + "HyperHelium5", // HyperHelium5 "Positron", // Positron "MuonPlus", // MuonPlus "PionMinus", // PionMinus @@ -157,10 +162,10 @@ class PIDExtended "AntiHyperhydrog4", // AntiHyperhydrog4 "XiPlus", // XiPlus "OmegaPlus", // OmegaPlus + "AntiHyperHelium4", // AntiHyperHelium4 + "AntiHyperHelium5", // AntiHyperHelium5 "Neutron", // Neutron "AntiNeutron", // AntiNeutron - "HyperHelium4", // HyperHelium4 - "AntiHyperHelium4", // AntiHyperHelium4 "Phi", // Phi "BZero", // BZero "BPlus", // BPlus @@ -244,6 +249,10 @@ class PIDExtended return HyperHelium4; case -o2::constants::physics::Pdg::kHyperHelium4: return AntiHyperHelium4; + case o2::constants::physics::Pdg::kHyperHelium5: + return HyperHelium5; + case -o2::constants::physics::Pdg::kHyperHelium5: + return AntiHyperHelium5; case 111: return PI0; case 22: From d4362cd455424ed4fbbe6d4b77eeb4c7a95683e8 Mon Sep 17 00:00:00 2001 From: Lucia Anna Tarasovicova Date: Thu, 12 Sep 2024 22:31:15 +0200 Subject: [PATCH 02/14] fix mc gen (#7669) --- PWGLF/Tasks/Strangeness/hStrangeCorrelation.cxx | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/hStrangeCorrelation.cxx b/PWGLF/Tasks/Strangeness/hStrangeCorrelation.cxx index fca804b8a13..1c77862d895 100644 --- a/PWGLF/Tasks/Strangeness/hStrangeCorrelation.cxx +++ b/PWGLF/Tasks/Strangeness/hStrangeCorrelation.cxx @@ -1498,15 +1498,17 @@ struct correlateStrangeness { Double_t getaassoc = assocParticle.eta(); Double_t gphiassoc = assocParticle.phi(); Double_t ptassoc = assocParticle.pt(); + Double_t deltaphi = ComputeDeltaPhi(gphitrigger, gphiassoc); + Double_t deltaeta = getatrigger - getaassoc; // skip if basic ranges not met - if (gphiassoc < axisRanges[0][0] || gphiassoc > axisRanges[0][1]) + if (deltaphi < axisRanges[0][0] || deltaphi > axisRanges[0][1]) continue; - if (getaassoc < axisRanges[1][0] || getaassoc > axisRanges[1][1]) + if (deltaeta < axisRanges[1][0] || deltaeta > axisRanges[1][1]) continue; if (ptassoc < axisRanges[2][0] || ptassoc > axisRanges[2][1]) continue; - histos.fill(HIST("ClosureTest/sameEvent/") + HIST(particlenames[index]), ComputeDeltaPhi(gphitrigger, gphiassoc), getatrigger - getaassoc, ptassoc, pttrigger, bestCollisionVtxZ, bestCollisionFT0Mpercentile); + histos.fill(HIST("ClosureTest/sameEvent/") + HIST(particlenames[index]), ComputeDeltaPhi(gphitrigger, gphiassoc), deltaeta, ptassoc, pttrigger, bestCollisionVtxZ, bestCollisionFT0Mpercentile); } } }); From d1d741f3f837fefa127cd6ae6550a9a1f9abf604 Mon Sep 17 00:00:00 2001 From: Mario Ciacco Date: Thu, 12 Sep 2024 22:32:48 +0200 Subject: [PATCH 03/14] add default cuts + re-compute multiplicity for run 2 (#7668) --- PWGLF/TableProducer/Nuspex/ebyeMaker.cxx | 86 +++++++++++++++++------- 1 file changed, 63 insertions(+), 23 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/ebyeMaker.cxx b/PWGLF/TableProducer/Nuspex/ebyeMaker.cxx index 383c3262779..63c918ee521 100644 --- a/PWGLF/TableProducer/Nuspex/ebyeMaker.cxx +++ b/PWGLF/TableProducer/Nuspex/ebyeMaker.cxx @@ -55,7 +55,8 @@ using BCsWithRun2Info = soa::Join; namespace { constexpr int kNpart = 2; -constexpr double trackSels[10]{/* 60, */ 80, 100, 2, 3, /* 4, */ 0.05, 0.1, /* 0.15, */ 0.5, 1, /* 1.5, */ 2, 3 /* , 4 */}; +constexpr float trackSels[10]{/* 60, */ 80, 100, 2, 3, /* 4, */ 0.05, 0.1, /* 0.15, */ 0.5, 1, /* 1.5, */ 2, 3 /* , 4 */}; +constexpr float dcaSels[3]{10., 10., 10.}; constexpr double betheBlochDefault[kNpart][6]{{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}, {-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}}; constexpr double estimatorsCorrelationCoef[2]{-0.669108, 1.04489}; constexpr double estimatorsSigmaPars[4]{0.933321, 0.0416976, -0.000936344, 8.92179e-06}; @@ -65,6 +66,7 @@ constexpr double partPdg[kNpart]{2212, o2::constants::physics::kDeuteron}; static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; static const std::vector particleNamesPar{"p", "d"}; static const std::vector trackSelsNames{"tpcClsMid", "tpcClsTight", "chi2TpcTight", "chi2TpcMid", "dcaxyTight", "dcaxyMid", "dcazTight", "dcazMid", "tpcNsigmaTight", "tpcNsigmaMid"}; +static const std::vector dcaSelsNames{"dcaxy", "dcaz", "dca"}; static const std::vector particleName{"p"}; std::array, kNpart> tofMass; void momTotXYZ(std::array& momA, std::array const& momB, std::array const& momC) @@ -177,6 +179,8 @@ struct tagRun2V0MCalibration { TH1* mhVtxAmpCorrV0A = nullptr; TH1* mhVtxAmpCorrV0C = nullptr; TH1* mhMultSelCalib = nullptr; + float mMCScalePars[6] = {0.0}; + TFormula* mMCScale = nullptr; } Run2V0MInfo; struct ebyeMaker { @@ -220,6 +224,7 @@ struct ebyeMaker { Configurable etaMaxV0dau{"etaMaxV0dau", 0.8f, "maximum eta V0 daughters"}; Configurable fillOnlySignal{"fillOnlySignal", false, "fill histograms only for true signal candidates (MC)"}; + Configurable genName{"genname", "", "Genearator name: HIJING, PYTHIA8, ... Default: \"\""}; Configurable kINT7Intervals{"kINT7Intervals", false, "toggle kINT7 trigger selection in the 10-30% and 50-90% centrality intervals (2018 Pb-Pb)"}; Configurable kUseTPCPileUpCut{"kUseTPCPileUpCut", false, "toggle strong correlation cuts (Run 2)"}; @@ -239,7 +244,8 @@ struct ebyeMaker { Configurable trackNcrossedRows{"trackNcrossedRows", 70, "Minimum number of crossed TPC rows"}; Configurable trackNclusItsCut{"trackNclusITScut", 2, "Minimum number of ITS clusters"}; Configurable trackNclusTpcCut{"trackNclusTPCcut", 60, "Minimum number of TPC clusters"}; - Configurable trackDcaCut{"trackDcaCut", 0.1f, "DCA antid to PV"}; + Configurable trackChi2Cut{"trackChi2Cut", 4.f, "Maximum chi2/ncls in TPC"}; + Configurable> cfgDcaSels{"cfgDcaSels", {dcaSels, 1, 3, particleName, dcaSelsNames}, "DCA selections"}; Configurable v0trackNcrossedRows{"v0trackNcrossedRows", 100, "Minimum number of crossed TPC rows for V0 daughter"}; Configurable v0trackNclusItsCut{"v0trackNclusITScut", 0, "Minimum number of ITS clusters for V0 daughter"}; @@ -272,7 +278,7 @@ struct ebyeMaker { Configurable antidItsClsSizeCut{"antidItsClsSizeCut", 1.e-10f, "cluster size cut for antideuterons"}; Configurable antidPtItsClsSizeCut{"antidPtItsClsSizeCut", 10.f, "pt for cluster size cut for antideuterons"}; - Configurable> cfgTrackSels{"cfgTrackSels", {trackSels, 1, 10, particleName, trackSelsNames}, "Track selections"}; + Configurable> cfgTrackSels{"cfgTrackSels", {trackSels, 1, 10, particleName, trackSelsNames}, "Track selections"}; std::array ptMin; std::array ptTof; @@ -328,7 +334,7 @@ struct ebyeMaker { track.tpcNClsFound() < trackNclusTpcCut || track.tpcNClsCrossedRows() < trackNcrossedRows || track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() || - track.tpcChi2NCl() > 4.f || + track.tpcChi2NCl() > trackChi2Cut || track.itsChi2NCl() > 36.f) { return false; } @@ -374,9 +380,28 @@ struct ebyeMaker { TH1* h = reinterpret_cast(callst->FindObject(ccdbhname)); return h; }; + auto getformulaccdb = [callst](const char* ccdbhname) { + TFormula* f = reinterpret_cast(callst->FindObject(ccdbhname)); + return f; + }; Run2V0MInfo.mhVtxAmpCorrV0A = getccdb("hVtx_fAmplitude_V0A_Normalized"); Run2V0MInfo.mhVtxAmpCorrV0C = getccdb("hVtx_fAmplitude_V0C_Normalized"); Run2V0MInfo.mhMultSelCalib = getccdb("hMultSelCalib_V0M"); + Run2V0MInfo.mMCScale = getformulaccdb(TString::Format("%s-V0M", genName->c_str()).Data()); + if ((Run2V0MInfo.mhVtxAmpCorrV0A != nullptr) && (Run2V0MInfo.mhVtxAmpCorrV0C != nullptr) && (Run2V0MInfo.mhMultSelCalib != nullptr)) { + if (genName->length() != 0) { + if (Run2V0MInfo.mMCScale != nullptr) { + for (int ixpar = 0; ixpar < 6; ++ixpar) { + Run2V0MInfo.mMCScalePars[ixpar] = Run2V0MInfo.mMCScale->GetParameter(ixpar); + } + } else { + LOGF(fatal, "MC Scale information from V0M for run %d not available", bc.runNumber()); + } + } + Run2V0MInfo.mCalibrationStored = true; + } else { + LOGF(fatal, "Calibration information from V0M for run %d corrupted", bc.runNumber()); + } } else { auto grpmagPath{"GLO/Config/GRPMagField"}; grpmag = ccdb->getForTimeStamp("GLO/Config/GRPMagField", timestamp); @@ -423,9 +448,19 @@ struct ebyeMaker { multFV0C += amplitude; } - float v0m = multFV0A * Run2V0MInfo.mhVtxAmpCorrV0A->GetBinContent(Run2V0MInfo.mhVtxAmpCorrV0A->FindFixBin(zvtx)) + - multFV0C * Run2V0MInfo.mhVtxAmpCorrV0C->GetBinContent(Run2V0MInfo.mhVtxAmpCorrV0C->FindFixBin(zvtx)); + float v0m = -1; + auto scaleMC = [](float x, float pars[6]) { + return pow(((pars[0] + pars[1] * pow(x, pars[2])) - pars[3]) / pars[4], 1.0f / pars[5]); + }; + if (Run2V0MInfo.mMCScale != nullptr) { + float multFV0M = multFV0A + multFV0C; + v0m = scaleMC(multFV0M, Run2V0MInfo.mMCScalePars); + LOGF(debug, "Unscaled v0m: %f, scaled v0m: %f", multFV0M, v0m); + } else { + v0m = multFV0A * Run2V0MInfo.mhVtxAmpCorrV0A->GetBinContent(Run2V0MInfo.mhVtxAmpCorrV0A->FindFixBin(zvtx)) + + multFV0C * Run2V0MInfo.mhVtxAmpCorrV0C->GetBinContent(Run2V0MInfo.mhVtxAmpCorrV0C->FindFixBin(zvtx)); + } return v0m; } @@ -549,7 +584,10 @@ struct ebyeMaker { auto dca = std::hypot(dcaInfo[0], dcaInfo[1]); auto trackPt = trackParCov.getPt(); auto trackEta = trackParCov.getEta(); - if (dca > trackDcaCut) { + if (dca > cfgDcaSels->get("dca")) { // dca + continue; + } + if (std::abs(dcaInfo[0]) > cfgDcaSels->get("dcaxy") || std::abs(dcaInfo[1]) > cfgDcaSels->get("dcaz")) { // dcaxy and dcaz continue; } histos.fill(HIST("QA/tpcSignal"), track.tpcInnerParam(), track.tpcSignal()); @@ -780,7 +818,7 @@ struct ebyeMaker { auto mcTrack = mcLab.template mcParticle_as(); if (std::abs(mcTrack.pdgCode()) != partPdg[iP]) continue; - if (((mcTrack.flags() & 0x8) && doprocessMcRun2) || (mcTrack.flags() & 0x2) || (mcTrack.flags() & 0x1)) + if (((mcTrack.flags() & 0x8) && (doprocessMcRun2 || doprocessMiniMcRun2)) || (mcTrack.flags() & 0x2) || (mcTrack.flags() & 0x1)) continue; if (!mcTrack.isPhysicalPrimary()) continue; @@ -813,7 +851,7 @@ struct ebyeMaker { } if (!posMother.isPhysicalPrimary() && !posMother.has_mothers()) continue; - if (((posMother.flags() & 0x8) && doprocessMcRun2) || (posMother.flags() & 0x2) || (posMother.flags() & 0x1)) + if (((posMother.flags() & 0x8) && (doprocessMcRun2 || doprocessMiniMcRun2)) || (posMother.flags() & 0x2) || (posMother.flags() & 0x1)) continue; auto genPt = std::hypot(posMother.px(), posMother.py()); @@ -836,7 +874,7 @@ struct ebyeMaker { if (std::abs(genEta) > etaMax) { continue; } - if (((mcPart.flags() & 0x8) && doprocessMcRun2) || (mcPart.flags() & 0x2) || (mcPart.flags() & 0x1)) + if (((mcPart.flags() & 0x8) && (doprocessMcRun2 || doprocessMiniMcRun2)) || (mcPart.flags() & 0x2) || (mcPart.flags() & 0x1)) continue; auto pdgCode = mcPart.pdgCode(); if (std::abs(pdgCode) == 3122) { @@ -861,7 +899,7 @@ struct ebyeMaker { if (it != candidateV0s.end()) { continue; } else { - LOGF(info, "not found!"); + LOGF(debug, "not found!"); candidateV0s.emplace_back(candV0); } } else if (std::abs(pdgCode) == partPdg[0] || std::abs(pdgCode) == partPdg[1]) { @@ -948,7 +986,7 @@ struct ebyeMaker { } PROCESS_SWITCH(ebyeMaker, processRun3, "process (Run 3)", false); - void processRun2(soa::Join const& collisions, TracksFull const& tracks, aod::V0s const& V0s, BCsWithRun2Info const&) + void processRun2(soa::Join const& collisions, TracksFull const& tracks, aod::V0s const& V0s, aod::FV0As const& fv0as, aod::FV0Cs const& fv0cs, BCsWithRun2Info const&) { for (const auto& collision : collisions) { auto bc = collision.bc_as(); @@ -963,8 +1001,9 @@ struct ebyeMaker { if (kUseTPCPileUpCut && !(bc.eventCuts() & BIT(aod::Run2EventCuts::kTPCPileUp))) continue; - auto centrality = collision.centRun2V0M(); - if (!(collision.sel7() && collision.alias_bit(kINT7)) && (!kINT7Intervals || (kINT7Intervals && ((centrality >= 10 && centrality < 30) || centrality > 50)))) + float v0m = getV0M(bc.globalIndex(), collision.posZ(), fv0as, fv0cs); + float cV0M = Run2V0MInfo.mhMultSelCalib->GetBinContent(Run2V0MInfo.mhMultSelCalib->FindFixBin(v0m)); + if (!(collision.sel7() && collision.alias_bit(kINT7)) && (!kINT7Intervals || (kINT7Intervals && ((cV0M >= 10 && cV0M < 30) || cV0M > 50)))) continue; auto centralityCl0 = collision.centRun2CL0(); @@ -972,7 +1011,7 @@ struct ebyeMaker { const auto& x = centralityCl0; const double center = estimatorsCorrelationCoef[0] + estimatorsCorrelationCoef[1] * x; const double sigma = estimatorsSigmaPars[0] + estimatorsSigmaPars[1] * x + estimatorsSigmaPars[2] * std::pow(x, 2) + estimatorsSigmaPars[3] * std::pow(x, 3); - if (centrality < center - deltaEstimatorNsigma[0] * sigma || centrality > center + deltaEstimatorNsigma[1] * sigma) { + if (cV0M < center - deltaEstimatorNsigma[0] * sigma || cV0M > center + deltaEstimatorNsigma[1] * sigma) { continue; } } @@ -984,12 +1023,12 @@ struct ebyeMaker { V0Table_thisCollision.bindExternalIndices(&tracks); auto multTracklets = collision.multTracklets(); - fillRecoEvent(collision, tracks, V0Table_thisCollision, centrality); + fillRecoEvent(collision, tracks, V0Table_thisCollision, cV0M); - histos.fill(HIST("QA/V0MvsCL0"), centralityCl0, centrality); - histos.fill(HIST("QA/trackletsVsV0M"), centrality, multTracklets); + histos.fill(HIST("QA/V0MvsCL0"), centralityCl0, cV0M); + histos.fill(HIST("QA/trackletsVsV0M"), cV0M, multTracklets); - collisionEbyeTable(centrality, collision.posZ()); + collisionEbyeTable(cV0M, collision.posZ()); for (auto& candidateV0 : candidateV0s) { lambdaEbyeTable( @@ -1139,7 +1178,7 @@ struct ebyeMaker { } PROCESS_SWITCH(ebyeMaker, processMcRun3, "process MC (Run 3)", false); - void processMcRun2(soa::Join const& collisions, aod::McCollisions const& /*mcCollisions*/, TracksFull const& tracks, aod::V0s const& V0s, aod::McParticles const& mcParticles, aod::McTrackLabels const& mcLab, BCsWithRun2Info const&) + void processMcRun2(soa::Join const& collisions, aod::McCollisions const& /*mcCollisions*/, TracksFull const& tracks, aod::V0s const& V0s, aod::FV0As const& fv0as, aod::FV0Cs const& fv0cs, aod::McParticles const& mcParticles, aod::McTrackLabels const& mcLab, BCsWithRun2Info const&) { for (auto& collision : collisions) { auto bc = collision.bc_as(); @@ -1151,7 +1190,8 @@ struct ebyeMaker { if (!(bc.eventCuts() & BIT(aod::Run2EventCuts::kAliEventCutsAccepted))) continue; - auto centrality = collision.centRun2V0M(); + float v0m = getV0M(bc.globalIndex(), collision.posZ(), fv0as, fv0cs); + float cV0M = Run2V0MInfo.mhMultSelCalib->GetBinContent(Run2V0MInfo.mhMultSelCalib->FindFixBin(v0m)); histos.fill(HIST("QA/zVtx"), collision.posZ()); @@ -1159,10 +1199,10 @@ struct ebyeMaker { auto V0Table_thisCollision = V0s.sliceBy(perCollisionV0, collIdx); V0Table_thisCollision.bindExternalIndices(&tracks); - fillMcEvent(collision, tracks, V0Table_thisCollision, centrality, mcParticles, mcLab); + fillMcEvent(collision, tracks, V0Table_thisCollision, cV0M, mcParticles, mcLab); fillMcGen(mcParticles, mcLab, collision.mcCollisionId()); - collisionEbyeTable(centrality, collision.posZ()); + collisionEbyeTable(cV0M, collision.posZ()); for (auto& candidateV0 : candidateV0s) { mcLambdaEbyeTable( From 8a283f076d0a8272c940b7d27ea929b41180ac28 Mon Sep 17 00:00:00 2001 From: Giorgio Alberto Lucia <87222843+GiorgioAlbertoLucia@users.noreply.github.com> Date: Thu, 12 Sep 2024 22:34:01 +0200 Subject: [PATCH 04/14] [PWGLF] removed global track selections, removed cuts in invariant mass, removed redundant histograms fillings, added columns to table, added zorro (#7642) * optimised code, removed global track selections, removed cuts in invariant mass, removed redundant histograms fillings --------- Co-authored-by: Francesco Mazzaschi --- PWGLF/DataModel/LFLithium4Tables.h | 14 + PWGLF/TableProducer/Nuspex/CMakeLists.txt | 2 +- .../TableProducer/Nuspex/lithium4analysis.cxx | 894 ++++++++++-------- 3 files changed, 496 insertions(+), 414 deletions(-) diff --git a/PWGLF/DataModel/LFLithium4Tables.h b/PWGLF/DataModel/LFLithium4Tables.h index 50c5acb5cca..1e65ffcf23b 100644 --- a/PWGLF/DataModel/LFLithium4Tables.h +++ b/PWGLF/DataModel/LFLithium4Tables.h @@ -44,6 +44,8 @@ DECLARE_SOA_COLUMN(InnerParamTPCPr, innerParamTPCPr, float); DECLARE_SOA_COLUMN(NClsTPCHe3, nClsTPCHe3, uint8_t); DECLARE_SOA_COLUMN(NSigmaTPCHe3, nSigmaTPCHe3, float); DECLARE_SOA_COLUMN(NSigmaTPCPr, nSigmaTOFPr, float); +DECLARE_SOA_COLUMN(Chi2TPCHe3, chi2TPCHe3, float); +DECLARE_SOA_COLUMN(Chi2TPCPr, chi2TPCPr, float); DECLARE_SOA_COLUMN(MassTOFHe3, massTOFHe3, float); DECLARE_SOA_COLUMN(MassTOFPr, massTOFPr, float); DECLARE_SOA_COLUMN(PIDtrkHe3, pidTrkHe3, uint32_t); @@ -59,7 +61,11 @@ DECLARE_SOA_COLUMN(IsBkgLS, isBkgLS, bool); DECLARE_SOA_COLUMN(IsBkgEM, isBkgEM, bool); DECLARE_SOA_COLUMN(PtMCHe3, ptMCHe3, float); +DECLARE_SOA_COLUMN(EtaMCHe3, etaMCHe3, float); +DECLARE_SOA_COLUMN(PhiMCHe3, phiMCHe3, float); DECLARE_SOA_COLUMN(PtMCPr, ptMCPr, float); +DECLARE_SOA_COLUMN(EtaMCPr, etaMCPr, float); +DECLARE_SOA_COLUMN(PhiMCPr, phiMCPr, float); DECLARE_SOA_COLUMN(SignedPtMC, signedPtMC, float); DECLARE_SOA_COLUMN(MassMC, massMC, float); @@ -83,6 +89,8 @@ DECLARE_SOA_TABLE(Lithium4Table, "AOD", "LITHIUM4TABLE", Lithium4TablesNS::NClsTPCHe3, Lithium4TablesNS::NSigmaTPCHe3, Lithium4TablesNS::NSigmaTPCPr, + Lithium4TablesNS::Chi2TPCHe3, + Lithium4TablesNS::Chi2TPCPr, Lithium4TablesNS::MassTOFHe3, Lithium4TablesNS::MassTOFPr, Lithium4TablesNS::PIDtrkHe3, @@ -111,6 +119,8 @@ DECLARE_SOA_TABLE(Lithium4TableMC, "AOD", "LITHIUM4TABLEMC", Lithium4TablesNS::NClsTPCHe3, Lithium4TablesNS::NSigmaTPCHe3, Lithium4TablesNS::NSigmaTPCPr, + Lithium4TablesNS::Chi2TPCHe3, + Lithium4TablesNS::Chi2TPCPr, Lithium4TablesNS::MassTOFHe3, Lithium4TablesNS::MassTOFPr, Lithium4TablesNS::PIDtrkHe3, @@ -122,7 +132,11 @@ DECLARE_SOA_TABLE(Lithium4TableMC, "AOD", "LITHIUM4TABLEMC", Lithium4TablesNS::IsBkgLS, Lithium4TablesNS::IsBkgEM, Lithium4TablesNS::PtMCHe3, + Lithium4TablesNS::EtaMCHe3, + Lithium4TablesNS::PhiMCHe3, Lithium4TablesNS::PtMCPr, + Lithium4TablesNS::EtaMCPr, + Lithium4TablesNS::PhiMCPr, Lithium4TablesNS::SignedPtMC, Lithium4TablesNS::MassMC) diff --git a/PWGLF/TableProducer/Nuspex/CMakeLists.txt b/PWGLF/TableProducer/Nuspex/CMakeLists.txt index 3d2e57d3878..5f9f4380def 100644 --- a/PWGLF/TableProducer/Nuspex/CMakeLists.txt +++ b/PWGLF/TableProducer/Nuspex/CMakeLists.txt @@ -46,7 +46,7 @@ o2physics_add_dpl_workflow(nucleustreecreator o2physics_add_dpl_workflow(lithium4analysis SOURCES lithium4analysis.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(nuclei-spectra diff --git a/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx b/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx index 20e717836a1..39ba20d9cfd 100644 --- a/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx +++ b/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx @@ -32,22 +32,33 @@ #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" +// #include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" #include "Framework/HistogramRegistry.h" +#include "Common/DataModel/Multiplicity.h" #include "Framework/StepTHn.h" -#include "ReconstructionDataFormats/Track.h" #include "Common/DataModel/PIDResponse.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/Centrality.h" +#include "ReconstructionDataFormats/Track.h" +#include "Common/Core/TrackSelection.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/DataModel/EventSelection.h" #include "Common/Core/trackUtilities.h" #include "Common/Core/RecoDecay.h" -#include "Common/Core/TrackSelection.h" #include "Framework/ASoAHelpers.h" #include "DataFormatsTPC/BetheBlochAleph.h" #include "CCDB/BasicCCDBManager.h" +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" + +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsBase/GeometryManager.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" + +#include "Common/Core/PID/TPCPIDResponse.h" #include "PWGLF/DataModel/LFLithium4Tables.h" +#include "PWGLF/Utils/svPoolCreator.h" using namespace o2; using namespace o2::framework; @@ -61,13 +72,12 @@ static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", constexpr float he3Mass = o2::constants::physics::MassHelium3; constexpr float protonMass = o2::constants::physics::MassProton; -constexpr int lithium4PDG = 1000030040; -constexpr int protonPDG = 2212; -constexpr int he3PDG = 1000020030; +constexpr int li4PDG = 1000030040; +constexpr int prPDG = 2212; +constexpr int hePDG = 1000020030; enum Selections { kNoCuts = 0, - kGlobalTrack, kTrackCuts, kPID, kAll @@ -75,9 +85,11 @@ enum Selections { } // namespace -struct lithium4Candidate { +struct Lithium4Candidate { float sign = 0.f; + std::array momHe3 = {99.f, 99.f, 99.f}; + std::array momPr = {99.f, 99.f, 99.f}; float recoPtHe3() const { return sign * std::hypot(momHe3[0], momHe3[1]); } float recoPhiHe3() const { return std::atan2(momHe3[1], momHe3[0]); } @@ -86,33 +98,30 @@ struct lithium4Candidate { float recoPhiPr() const { return std::atan2(momPr[1], momPr[0]); } float recoEtaPr() const { return std::asinh(momPr[2] / recoPtPr()); } - std::array momHe3 = {99.f, 99.f, 99.f}; - std::array momPr = {99.f, 99.f, 99.f}; - - uint32_t PIDtrkHe3 = 0xFFFFF; // PID in tracking - uint32_t PIDtrkPr = 0xFFFFF; - - float nSigmaHe3 = -10; - float nSigmaPr = -10; - float massTOFHe3 = -10; - float massTOFPr = -10; + float DCAxyHe3 = -10.f; + float DCAzHe3 = -10.f; + float DCAxyPr = -10.f; + float DCAzPr = -10.f; - float DCAxyHe3 = -10; - float DCAzHe3 = -10; - float DCAxyPr = -10; - float DCAzPr = -10; uint16_t tpcSignalHe3 = 0u; float momHe3TPC = -99.f; uint16_t tpcSignalPr = 0u; float momPrTPC = -99.f; - float invMass = -10.f; - uint32_t itsClSizeHe3 = 0u; - uint32_t itsClSizePr = 0u; uint8_t nTPCClustersHe3 = 0u; + float nSigmaHe3 = -10.f; + float nSigmaPr = -10.f; + float chi2TPCHe3 = -10.f; + float chi2TPCPr = -10.f; - float momHe3MC = -99.f; - float momPrMC = -99.f; + float massTOFHe3 = -10; + float massTOFPr = -10; + + uint32_t PIDtrkHe3 = 0xFFFFF; // PID in tracking + uint32_t PIDtrkPr = 0xFFFFF; + + uint32_t itsClSizeHe3 = 0u; + uint32_t itsClSizePr = 0u; uint8_t sharedClustersHe3 = 0u; uint8_t sharedClustersPr = 0u; @@ -120,66 +129,153 @@ struct lithium4Candidate { bool isBkgUS = false; bool isBkgEM = false; + int trackIDHe3 = -1; + int trackIDPr = -1; + + float momHe3MC = -99.f; + float etaHe3MC = -99.f; + float phiHe3MC = -99.f; + float momPrMC = -99.f; + float etaPrMC = -99.f; + float phiPrMC = -99.f; + float l4PtMC = -99.f; float l4MassMC = -10.f; + + float invMass = -10.f; }; struct lithium4analysis { - Produces outputDataTable; - Produces outputMCTable; - - std::vector l4Candidates; + Produces m_outputDataTable; + Produces m_outputMCTable; + + // Selections + Configurable setting_cutVertex{"setting_cutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable setting_cutPT{"setting_cutPT", 0.2, "PT cut on daughter track"}; + Configurable setting_cutMaxPrPT{"setting_cutMaxPrPT", 1.8, "Max PT cut on proton"}; + Configurable setting_cutEta{"setting_cutEta", 0.8, "Eta cut on daughter track"}; + Configurable setting_cutDCAxy{"setting_cutDCAxy", 2.0f, "DCAxy range for tracks"}; + Configurable setting_cutDCAz{"setting_cutDCAz", 2.0f, "DCAz range for tracks"}; + Configurable setting_nsigmaCutTPC{"setting_nsigmaCutTPC", 3.0, "Value of the TPC Nsigma cut"}; + Configurable setting_nsigmaCutTOF{"setting_nsigmaCutTOF", 3.0, "Value of the TOF Nsigma cut"}; + Configurable setting_noMixedEvents{"setting_noMixedEvents", 5, "Number of mixed events per event"}; + Configurable setting_enableBkgUS{"setting_enableBkgUS", false, "Enable US background"}; + Configurable setting_isMC{"setting_isMC", false, "Run MC"}; + Configurable setting_skimmedProcessing{"setting_skimmedProcessing", false, "Skimmed dataset processing"}; + + // CCDB options + Configurable setting_ccdburl{"setting_ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable setting_grpPath{"setting_grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable setting_grpmagPath{"setting_grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable setting_lutPath{"setting_lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; + Configurable setting_geoPath{"setting_geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + Configurable setting_pidPath{"setting_pidPath", "", "Path to the PID response object"}; + + Configurable> setting_BetheBlochParams{"setting_BetheBlochParams", {betheBlochDefault[0], 1, 6, {"He3"}, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for He3"}; + Configurable setting_compensatePIDinTracking{"setting_compensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"}; + Configurable setting_materialCorrection{"setting_materialCorrection", static_cast(o2::base::Propagator::MatCorrType::USEMatCorrNONE), "Material correction type"}; + + using McIter = aod::McParticles::iterator; + using CollBracket = o2::math_utils::Bracket; + Filter collZfilter = nabs(aod::collision::posZ) < setting_cutVertex; + using CollisionsFull = soa::Filtered>; //, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms>>; + using CollisionsFullMC = soa::Filtered>; //, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms>>; + using TrackCandidates = soa::Join; + using TrackCandidatesMC = soa::Join; + + o2::pid::tof::Beta m_responseBeta; + o2::pid::tof::Beta m_responseBetaMC; + + Preslice m_perCol = aod::track::collisionId; + Preslice m_perColMC = aod::track::collisionId; + // binning for EM background + ConfigurableAxis axisVertex{"axisVertex", {30, -10, 10}, "vertex axis for bin"}; + using BinningType = ColumnBinningPolicy; + BinningType binningOnPositions{{axisVertex}, true}; SliceCache cache; - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - Configurable cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"}; - // events - Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; - // track - Configurable cfgCutPT{"cfgCutPT", 0.2, "PT cut on daughter track"}; - Configurable cfgCutMaxPrPT{"cfgCutMaxPrPT", 1.8, "Max PT cut on proton"}; - Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut on daughter track"}; - Configurable cfgCutDCAxy{"cfgCutDCAxy", 2.0f, "DCAxy range for tracks"}; - Configurable cfgCutDCAz{"cfgCutDCAz", 2.0f, "DCAz range for tracks"}; - Configurable nsigmaCutTPC{"nsigmacutTPC", 3.0, "Value of the TPC Nsigma cut"}; - Configurable nsigmaCutTOF{"nsigmaCutTOF", 3.0, "Value of the TOF Nsigma cut"}; - Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 5, "Number of mixed events per event"}; - Configurable cfgEnableBkgUS{"cfgEnableBkgUS", false, "Enable US background"}; - - // bethe bloch parameters - std::array mBBparamsHe; - Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {betheBlochDefault[0], 1, 6, {"He3"}, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for He3"}; - // MC - Configurable isMC{"isMC", false, "Run MC"}; + SameKindPair m_pair{binningOnPositions, setting_noMixedEvents, -1, &cache}; + + std::array m_BBparamsHe; + + std::vector m_recoCollisionIDs; + std::vector m_goodCollisions; + std::vector m_trackPairs; + + int m_runNumber; + Service m_ccdb; + Zorro m_zorro; + OutputObj m_zorroSummary{"zorroSummary"}; + + HistogramRegistry m_qaRegistry{ + "QA", + { + {"hVtxZ", "Vertex distribution in Z;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}}, + {"hNcontributor", "Number of primary vertex contributor", {HistType::kTH1F, {{2000, 0.0f, 2000.0f}}}}, + {"hTrackSel", "Accepted tracks", {HistType::kTH1F, {{Selections::kAll, -0.5, static_cast(Selections::kAll) - 0.5}}}}, + {"hEvents", "; Events;", {HistType::kTH1F, {{3, -0.5, 2.5}}}}, + {"hEmptyPool", "svPoolCreator did not find track pairs false/true", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"hDCAxyHe3", ";DCA_{xy} (cm)", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}}, + {"hDCAzHe3", ";DCA_{z} (cm)", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}}, + {"hLitInvMass", "; M(^{3}He + p) (GeV/#it{c}^{2})", {HistType::kTH1F, {{50, 3.74f, 3.85f}}}}, + {"hHe3Pt", "#it{p}_{T} distribution; #it{p}_{T} (GeV/#it{c})", {HistType::kTH1F, {{200, -6.0f, 6.0f}}}}, + {"hProtonPt", "Pt distribution; #it{p}_{T} (GeV/#it{c})", {HistType::kTH1F, {{200, -3.0f, 3.0f}}}}, + {"h2dEdxHe3candidates", "dEdx distribution; Signed #it{p} (GeV/#it{c}); dE/dx (a.u.)", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {100, 0.0f, 2000.0f}}}}, + {"h2NsigmaHe3TPC", "NsigmaHe3 TPC distribution; Signed #it{p}/#it{z} (GeV/#it{c}); n#sigma_{TPC}({}^{3}He)", {HistType::kTH2F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}}}, + {"h2NsigmaProtonTPC", "NsigmaProton TPC distribution; Signed #it{p}/#it{z} (GeV/#it{c}); n#sigma_{TPC}(p)", {HistType::kTH2F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}}}, + {"h2NsigmaProtonTOF", "NsigmaProton TOF distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TOF}(p)", {HistType::kTH2F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}}}, + }, + OutputObjHandlingPolicy::AnalysisObject, + false, + true}; + void init(o2::framework::InitContext&) { + LOG(info) << "Initializing lithium4 analysis"; + m_zorroSummary.setObject(m_zorro.getZorroSummary()); + m_runNumber = 0; - histos.add("hCentrality", "Centrality distribution", kTH1F, {{2001, -0.5, 2000.5}}); - histos.add("hVtxZ", "Vertex distribution in Z;Z (cm)", kTH1F, {{400, -20.0, 20.0}}); - histos.add("hNcontributor", "Number of primary vertex contributor", kTH1F, {{2000, 0.0f, 2000.0f}}); - histos.add("hDCAxyHe3", ";DCA_{xy} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hDCAzHe3", ";DCA_{z} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hLitInvMass", "; M(^{3}He + p) (GeV/#it{c}^{2})", kTH1F, {{50, 3.74f, 3.85f}}); - histos.add("hHe3Pt", "#it{p}_{T} distribution; #it{p}_{T} (GeV/#it{c})", kTH1F, {{200, 0.0f, 6.0f}}); - histos.add("hProtonPt", "Pt distribution; #it{p}_{T} (GeV/#it{c})", kTH1F, {{200, 0.0f, 3.0f}}); - histos.add("h2dEdxHe3candidates", "dEdx distribution; Signed #it{p} (GeV/#it{c}); dE/dx (a.u.)", kTH2F, {{200, -5.0f, 5.0f}, {100, 0.0f, 2000.0f}}); - histos.add("h2NsigmaHe3TPC", "NsigmaHe3 TPC distribution; Signed #it{p}/#it{z} (GeV/#it{c}); n#sigma_{TPC}({}^{3}He)", kTH2F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}); - histos.add("h2NsigmaProtonTPC", "NsigmaProton TPC distribution; Signed #it{p}/#it{z} (GeV/#it{c}); n#sigma_{TPC}(p)", kTH2F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}); - histos.add("h2NsigmaProtonTOF", "NsigmaProton TOF distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TOF}(p)", kTH2F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}); - histos.add("hTrackSel", "Accepted tracks", kTH1F, {{Selections::kAll, -0.5, static_cast(Selections::kAll) - 0.5}}); + LOG(info) << "Initializing CCDB"; + m_ccdb->setURL(setting_ccdburl); + m_ccdb->setCaching(true); + m_ccdb->setLocalObjectValidityChecking(); + m_ccdb->setFatalWhenNull(false); + LOG(info) << "Initializing PID response"; for (int i = 0; i < 5; i++) { - mBBparamsHe[i] = cfgBetheBlochParams->get("He3", Form("p%i", i)); + m_BBparamsHe[i] = setting_BetheBlochParams->get("He3", Form("p%i", i)); + } + m_BBparamsHe[5] = setting_BetheBlochParams->get("He3", "resolution"); + + std::vector selection_labels = {"All", "Track selection", "PID {}^{3}He"}; + for (int i = 0; i < Selections::kAll; i++) { + m_qaRegistry.get(HIST("hTrackSel"))->GetXaxis()->SetBinLabel(i + 1, selection_labels[i].c_str()); } - mBBparamsHe[5] = cfgBetheBlochParams->get("He3", "resolution"); - std::vector labels = {"All", "Global track", "Track selection", "PID {}^{3}He"}; + std::vector events_labels = {"All", "Selected", "Zorro He events"}; for (int i = 0; i < Selections::kAll; i++) { - histos.get(HIST("hTrackSel"))->GetXaxis()->SetBinLabel(i + 1, labels[i].c_str()); + m_qaRegistry.get(HIST("hEvents"))->GetXaxis()->SetBinLabel(i + 1, events_labels[i].c_str()); + } + } + + void initCCDB(const aod::BCsWithTimestamps::iterator& bc) + { + LOG(info) << "Initializing CCDB for run " << bc.runNumber(); + if (m_runNumber == bc.runNumber()) { + return; } + LOG(info) << "Initializing zorro for run " << bc.runNumber(); + if (setting_skimmedProcessing) { + m_zorro.initCCDB(m_ccdb.service, bc.runNumber(), bc.timestamp(), "fHe"); + m_zorro.populateHistRegistry(m_qaRegistry, bc.runNumber()); + } + LOG(info) << "Initializing CCDB for run " << bc.runNumber() << " done"; + m_runNumber = bc.runNumber(); } + // ================================================================================================================== + template bool selectionTrack(const T& candidate) { @@ -200,13 +296,13 @@ struct lithium4analysis { bool selectionPIDProton(const T& candidate) { if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaPr()) < nsigmaCutTOF && std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC) { - histos.fill(HIST("h2NsigmaProtonTPC"), candidate.tpcInnerParam(), candidate.tpcNSigmaPr()); - histos.fill(HIST("h2NsigmaProtonTOF"), candidate.p(), candidate.tofNSigmaPr()); + if (std::abs(candidate.tofNSigmaPr()) < setting_nsigmaCutTOF && std::abs(candidate.tpcNSigmaPr()) < setting_nsigmaCutTPC) { + m_qaRegistry.fill(HIST("h2NsigmaProtonTPC"), candidate.tpcInnerParam(), candidate.tpcNSigmaPr()); + m_qaRegistry.fill(HIST("h2NsigmaProtonTOF"), candidate.p(), candidate.tofNSigmaPr()); return true; } - } else if (std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC) { - histos.fill(HIST("h2NsigmaProtonTPC"), candidate.tpcInnerParam(), candidate.tpcNSigmaPr()); + } else if (std::abs(candidate.tpcNSigmaPr()) < setting_nsigmaCutTPC) { + m_qaRegistry.fill(HIST("h2NsigmaProtonTPC"), candidate.tpcInnerParam(), candidate.tpcNSigmaPr()); return true; } return false; @@ -216,11 +312,10 @@ struct lithium4analysis { float computeNSigmaHe3(const T& candidate) { bool heliumPID = candidate.pidForTracking() == o2::track::PID::Helium3 || candidate.pidForTracking() == o2::track::PID::Alpha; + float correctedTPCinnerParam = (heliumPID && setting_compensatePIDinTracking) ? candidate.tpcInnerParam() / 2.f : candidate.tpcInnerParam(); + float expTPCSignal = o2::tpc::BetheBlochAleph(static_cast(correctedTPCinnerParam * 2.f / constants::physics::MassHelium3), m_BBparamsHe[0], m_BBparamsHe[1], m_BBparamsHe[2], m_BBparamsHe[3], m_BBparamsHe[4]); - float correctedTPCinnerParam = (heliumPID && cfgCompensatePIDinTracking) ? candidate.tpcInnerParam() / 2.f : candidate.tpcInnerParam(); - float expTPCSignal = o2::tpc::BetheBlochAleph(static_cast(correctedTPCinnerParam * 2.f / constants::physics::MassHelium3), mBBparamsHe[0], mBBparamsHe[1], mBBparamsHe[2], mBBparamsHe[3], mBBparamsHe[4]); - - double resoTPC{expTPCSignal * mBBparamsHe[5]}; + double resoTPC{expTPCSignal * m_BBparamsHe[5]}; return static_cast((candidate.tpcSignal() - expTPCSignal) / resoTPC); } @@ -228,386 +323,388 @@ struct lithium4analysis { bool selectionPIDHe3(const T& candidate) { auto nSigmaHe3 = computeNSigmaHe3(candidate); - if (std::abs(nSigmaHe3) < nsigmaCutTPC) { + bool heliumPID = candidate.pidForTracking() == o2::track::PID::Helium3 || candidate.pidForTracking() == o2::track::PID::Alpha; + float correctedTPCinnerParam = (heliumPID && setting_compensatePIDinTracking) ? candidate.tpcInnerParam() / 2.f : candidate.tpcInnerParam(); + if (std::abs(nSigmaHe3) < setting_nsigmaCutTPC) { + m_qaRegistry.fill(HIST("h2NsigmaHe3TPC"), candidate.sign() * correctedTPCinnerParam, nSigmaHe3); return true; } return false; } - template - bool FillCandidateInfo(const T1& candidateHe3, const T2& candidatePr, bool mix, bool /*isMC*/ = false) - { - lithium4Candidate l4Cand; + // ================================================================================================================== - l4Cand.momHe3 = array{2 * candidateHe3.px(), 2 * candidateHe3.py(), 2 * candidateHe3.pz()}; - l4Cand.momPr = array{candidatePr.px(), candidatePr.py(), candidatePr.pz()}; - - float invMass = RecoDecay::m(array{l4Cand.momHe3, l4Cand.momPr}, array{he3Mass, protonMass}); - - if (invMass < 3.74 || invMass > 3.85 || candidatePr.pt() > cfgCutMaxPrPT) { - return false; + template + bool fillCandidateInfo(const T& trackHe3, const T& trackPr, Lithium4Candidate& li4cand, bool mix) + { + li4cand.momHe3 = std::array{trackHe3.px(), trackHe3.py(), trackHe3.pz()}; + li4cand.momPr = std::array{trackPr.px(), trackPr.py(), trackPr.pz()}; + float invMass = RecoDecay::m(array{li4cand.momHe3, li4cand.momPr}, std::array{o2::constants::physics::MassHelium3, o2::constants::physics::MassProton}); + + li4cand.sign = trackHe3.sign(); + + li4cand.DCAxyHe3 = trackHe3.dcaXY(); + li4cand.DCAzHe3 = trackHe3.dcaZ(); + li4cand.DCAxyPr = trackPr.dcaXY(); + li4cand.DCAzPr = trackPr.dcaZ(); + + li4cand.tpcSignalHe3 = trackHe3.tpcSignal(); + bool heliumPID = trackHe3.pidForTracking() == o2::track::PID::Helium3 || trackHe3.pidForTracking() == o2::track::PID::Alpha; + float correctedTPCinnerParamHe3 = (heliumPID && setting_compensatePIDinTracking) ? trackHe3.tpcInnerParam() / 2.f : trackHe3.tpcInnerParam(); + li4cand.momHe3TPC = correctedTPCinnerParamHe3; + li4cand.tpcSignalPr = trackPr.tpcSignal(); + li4cand.momPrTPC = trackPr.tpcInnerParam(); + + li4cand.nTPCClustersHe3 = trackHe3.tpcNClsFound(); + li4cand.nSigmaHe3 = computeNSigmaHe3(trackHe3); + li4cand.nSigmaPr = trackPr.tpcNSigmaPr(); + li4cand.chi2TPCHe3 = trackHe3.tpcChi2NCl(); + li4cand.chi2TPCPr = trackPr.tpcChi2NCl(); + + li4cand.PIDtrkHe3 = trackHe3.pidForTracking(); + li4cand.PIDtrkPr = trackPr.pidForTracking(); + + li4cand.itsClSizeHe3 = trackHe3.itsClusterSizes(); + li4cand.itsClSizePr = trackPr.itsClusterSizes(); + + li4cand.sharedClustersHe3 = trackHe3.tpcNClsShared(); + li4cand.sharedClustersPr = trackPr.tpcNClsShared(); + + li4cand.isBkgUS = trackHe3.sign() * trackPr.sign() < 0; + li4cand.isBkgEM = mix; + + li4cand.invMass = invMass; + + li4cand.trackIDHe3 = trackHe3.globalIndex(); + li4cand.trackIDPr = trackPr.globalIndex(); + + if constexpr (isMC) { + if (trackHe3.hasTOF()) { + float beta = m_responseBetaMC.GetBeta(trackHe3); + beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + bool heliumPID = trackHe3.pidForTracking() == o2::track::PID::Helium3 || trackHe3.pidForTracking() == o2::track::PID::Alpha; + float correctedTPCinnerParamHe3 = (heliumPID && setting_compensatePIDinTracking) ? trackHe3.tpcInnerParam() / 2.f : trackHe3.tpcInnerParam(); + li4cand.massTOFHe3 = correctedTPCinnerParamHe3 * 2.f * std::sqrt(1.f / (beta * beta) - 1.f); + } + if (trackPr.hasTOF()) { + float beta = m_responseBetaMC.GetBeta(trackPr); + beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + li4cand.massTOFPr = trackPr.tpcInnerParam() * std::sqrt(1.f / (beta * beta) - 1.f); + } + } else { + if (trackHe3.hasTOF()) { + float beta = m_responseBeta.GetBeta(trackHe3); + beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + bool heliumPID = trackHe3.pidForTracking() == o2::track::PID::Helium3 || trackHe3.pidForTracking() == o2::track::PID::Alpha; + float correctedTPCinnerParamHe3 = (heliumPID && setting_compensatePIDinTracking) ? trackHe3.tpcInnerParam() / 2.f : trackHe3.tpcInnerParam(); + li4cand.massTOFHe3 = correctedTPCinnerParamHe3 * 2.f * std::sqrt(1.f / (beta * beta) - 1.f); + } + if (trackPr.hasTOF()) { + float beta = m_responseBeta.GetBeta(trackPr); + beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + li4cand.massTOFPr = trackPr.tpcInnerParam() * std::sqrt(1.f / (beta * beta) - 1.f); + } } - l4Cand.PIDtrkHe3 = candidateHe3.pidForTracking(); - l4Cand.PIDtrkPr = candidatePr.pidForTracking(); - - l4Cand.sign = candidateHe3.sign(); - - l4Cand.isBkgUS = candidateHe3.sign() * candidatePr.sign() < 0; - l4Cand.isBkgEM = mix; - - l4Cand.DCAxyHe3 = candidateHe3.dcaXY(); - l4Cand.DCAzHe3 = candidateHe3.dcaZ(); - l4Cand.DCAxyPr = candidatePr.dcaXY(); - l4Cand.DCAzPr = candidatePr.dcaZ(); - - bool heliumPID = candidateHe3.pidForTracking() == o2::track::PID::Helium3 || candidateHe3.pidForTracking() == o2::track::PID::Alpha; - float correctedTPCinnerParamHe3 = (heliumPID && cfgCompensatePIDinTracking) ? candidateHe3.tpcInnerParam() / 2.f : candidateHe3.tpcInnerParam(); - - l4Cand.tpcSignalHe3 = candidateHe3.tpcSignal(); - l4Cand.momHe3TPC = correctedTPCinnerParamHe3; - l4Cand.tpcSignalPr = candidatePr.tpcSignal(); - l4Cand.momPrTPC = candidatePr.tpcInnerParam(); - l4Cand.invMass = invMass; - - l4Cand.itsClSizeHe3 = candidateHe3.itsClusterSizes(); - l4Cand.itsClSizePr = candidatePr.itsClusterSizes(); - - l4Cand.nTPCClustersHe3 = candidateHe3.tpcNClsFound(); - - l4Cand.nSigmaHe3 = computeNSigmaHe3(candidateHe3); - l4Cand.nSigmaPr = candidatePr.tpcNSigmaPr(); - - l4Cand.sharedClustersHe3 = candidateHe3.tpcNClsShared(); - l4Cand.sharedClustersPr = candidatePr.tpcNClsShared(); - - l4Candidates.push_back(l4Cand); return true; } - template - void fillHistograms(const T& l4cand) + template + void fillCandidateInfoMC(const Mc& mctrackHe3, const Mc& mctrackPr, const Mc& mctrackMother, Lithium4Candidate& li4cand) { - int candSign = l4cand.sign; - histos.fill(HIST("hHe3Pt"), l4cand.recoPtHe3()); - histos.fill(HIST("hProtonPt"), l4cand.recoPtPr()); - histos.fill(HIST("hLitInvMass"), l4cand.invMass); - histos.fill(HIST("hDCAxyHe3"), l4cand.DCAxyHe3); - histos.fill(HIST("hDCAzHe3"), l4cand.DCAzHe3); - histos.fill(HIST("h2NsigmaHe3TPC"), candSign * l4cand.momHe3TPC, l4cand.nSigmaHe3); - histos.fill(HIST("h2NsigmaProtonTPC"), candSign * l4cand.momPrTPC, l4cand.nSigmaPr); - histos.fill(HIST("h2NsigmaProtonTOF"), l4cand.recoPtPr(), l4cand.nSigmaPr); + li4cand.momHe3MC = mctrackHe3.pt() * (mctrackHe3.pdgCode() > 0 ? 1 : -1); + li4cand.etaHe3MC = mctrackHe3.eta(); + li4cand.phiHe3MC = mctrackHe3.phi(); + li4cand.momPrMC = mctrackPr.pt() * (mctrackPr.pdgCode() > 0 ? 1 : -1); + li4cand.etaPrMC = mctrackPr.eta(); + li4cand.phiPrMC = mctrackPr.phi(); + li4cand.l4PtMC = mctrackMother.pt() * (mctrackMother.pdgCode() > 0 ? 1 : -1); + const double eLit = mctrackHe3.e() + mctrackPr.e(); + li4cand.l4MassMC = std::sqrt(eLit * eLit - mctrackMother.p() * mctrackMother.p()); } - Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; - Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPT); - Filter DCAcutFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz); + template + void pairTracksSameEvent(const T& tracks) + { + for (auto track0 : tracks) { - using EventCandidates = soa::Filtered>; - using TrackCandidates = soa::Filtered>; - using TrackCandidatesMC = soa::Filtered>; - o2::pid::tof::Beta responseBeta; - o2::pid::tof::Beta responseBetaMC; + m_qaRegistry.fill(HIST("hTrackSel"), Selections::kNoCuts); + bool heliumPID = track0.pidForTracking() == o2::track::PID::Helium3 || track0.pidForTracking() == o2::track::PID::Alpha; - Preslice perCol = aod::track::collisionId; - Preslice perColMC = aod::track::collisionId; + float correctedTPCinnerParam = (heliumPID && setting_compensatePIDinTracking) ? track0.tpcInnerParam() / 2.f : track0.tpcInnerParam(); + m_qaRegistry.fill(HIST("h2dEdxHe3candidates"), correctedTPCinnerParam * 2.f, track0.tpcSignal()); - // binning for EM background - ConfigurableAxis axisVertex{"axisVertex", {30, -10, 10}, "vertex axis for bin"}; - using BinningType = ColumnBinningPolicy; - BinningType binningOnPositions{{axisVertex}, true}; - SameKindPair pair{binningOnPositions, cfgNoMixedEvents, -1, &cache}; - - void processSameEvent(soa::Join const& collisions, TrackCandidates const& tracks, aod::BCs const&) - { - l4Candidates.clear(); - - for (auto& collision : collisions) { - if (!collision.sel8() || std::abs(collision.posZ()) > cfgCutVertex) { + if (!selectionTrack(track0)) { continue; } - histos.fill(HIST("hNcontributor"), collision.numContrib()); - histos.fill(HIST("hVtxZ"), collision.posZ()); + m_qaRegistry.fill(HIST("hTrackSel"), Selections::kTrackCuts); - const uint64_t collIdx = collision.globalIndex(); - auto TrackTable_thisCollision = tracks.sliceBy(perCol, collIdx); - TrackTable_thisCollision.bindExternalIndices(&tracks); - - for (auto track1 : TrackTable_thisCollision) { - - histos.fill(HIST("hTrackSel"), Selections::kNoCuts); - bool heliumPID = track1.pidForTracking() == o2::track::PID::Helium3 || track1.pidForTracking() == o2::track::PID::Alpha; + if (!selectionPIDHe3(track0)) { + continue; + } + m_qaRegistry.fill(HIST("hTrackSel"), Selections::kPID); - float correctedTPCinnerParam = (heliumPID && cfgCompensatePIDinTracking) ? track1.tpcInnerParam() / 2.f : track1.tpcInnerParam(); - histos.fill(HIST("h2dEdxHe3candidates"), correctedTPCinnerParam * 2.f, track1.tpcSignal()); + for (auto track1 : tracks) { - if (!track1.isGlobalTrackWoDCA()) { + if (track0 == track1) { continue; } - histos.fill(HIST("hTrackSel"), Selections::kGlobalTrack); - if (!selectionTrack(track1)) { - continue; + if (!setting_enableBkgUS) { + if (track0.sign() * track1.sign() < 0) { + continue; + } } - histos.fill(HIST("hTrackSel"), Selections::kTrackCuts); - if (!selectionPIDHe3(track1)) { + if (!selectionTrack(track1) || !selectionPIDProton(track1)) { continue; } - histos.fill(HIST("hTrackSel"), Selections::kPID); - - for (auto track2 : TrackTable_thisCollision) { - - if (track1 == track2) { - continue; - } - if (!cfgEnableBkgUS) { - if (track1.sign() * track2.sign() < 0) { - continue; - } - } - - if (!track2.isGlobalTrackWoDCA()) { - continue; - } - - if (!selectionTrack(track2)) { - continue; - } - - if (!selectionPIDProton(track2)) { - continue; - } - - if (!FillCandidateInfo(track1, track2, false)) { - continue; - } - // fill TOF info outside to avoide responseBeta crash - auto& cand = l4Candidates.back(); - if (track1.hasTOF()) { - float beta = responseBeta.GetBeta(track1); - beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked - bool heliumPID = track1.pidForTracking() == o2::track::PID::Helium3 || track1.pidForTracking() == o2::track::PID::Alpha; - float correctedTPCinnerParamHe3 = (heliumPID && cfgCompensatePIDinTracking) ? track1.tpcInnerParam() / 2.f : track1.tpcInnerParam(); - cand.massTOFHe3 = correctedTPCinnerParamHe3 * 2.f * std::sqrt(1.f / (beta * beta) - 1.f); - } - if (track2.hasTOF()) { - float beta = responseBeta.GetBeta(track2); - beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked - cand.massTOFPr = track2.tpcInnerParam() * std::sqrt(1.f / (beta * beta) - 1.f); - } - fillHistograms(cand); - } + SVCand trackPair; + trackPair.tr0Idx = track0.globalIndex(); + trackPair.tr1Idx = track1.globalIndex(); + m_trackPairs.push_back(trackPair); } } - - for (auto& l4Cand : l4Candidates) { - outputDataTable(l4Cand.recoPtHe3(), l4Cand.recoEtaHe3(), l4Cand.recoPhiHe3(), - l4Cand.recoPtPr(), l4Cand.recoEtaPr(), l4Cand.recoPhiPr(), - l4Cand.DCAxyHe3, l4Cand.DCAzHe3, l4Cand.DCAxyPr, l4Cand.DCAzPr, - l4Cand.tpcSignalHe3, l4Cand.momHe3TPC, l4Cand.tpcSignalPr, l4Cand.momPrTPC, - l4Cand.nTPCClustersHe3, - l4Cand.nSigmaHe3, l4Cand.nSigmaPr, l4Cand.massTOFHe3, l4Cand.massTOFPr, - l4Cand.PIDtrkHe3, l4Cand.PIDtrkPr, l4Cand.itsClSizeHe3, l4Cand.itsClSizePr, - l4Cand.sharedClustersHe3, l4Cand.sharedClustersPr, - l4Cand.isBkgUS, l4Cand.isBkgEM); - } } - PROCESS_SWITCH(lithium4analysis, processSameEvent, "Process Same event", false); - void processMixedEvent(EventCandidates& /*collisions*/, TrackCandidates const& /*tracks*/) + void pairTracksEventMixing() { - l4Candidates.clear(); - for (auto& [c1, tracks1, c2, tracks2] : pair) { - if (!c1.sel8()) { - continue; - } - if (!c2.sel8()) { + for (auto& [c1, tracks1, c2, tracks2] : m_pair) { + if (!c1.sel8() || !c2.sel8()) { continue; } - histos.fill(HIST("hNcontributor"), c1.numContrib()); - histos.fill(HIST("hVtxZ"), c1.posZ()); + m_qaRegistry.fill(HIST("hNcontributor"), c1.numContrib()); + m_qaRegistry.fill(HIST("hVtxZ"), c1.posZ()); for (auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (!t1.isGlobalTrackWoDCA()) { - continue; - } - - if (!selectionTrack(t1)) { - continue; - } - - if (!t2.isGlobalTrackWoDCA()) { - continue; - } - - if (!selectionTrack(t2)) { + if (!selectionTrack(t1) || !selectionTrack(t2)) { continue; } TrackCandidates::iterator he3Cand, protonCand; - bool passPID = false; if (selectionPIDHe3(t1) && selectionPIDProton(t2)) { he3Cand = t1, protonCand = t2; - passPID = true; - } - if (selectionPIDHe3(t2) && selectionPIDProton(t1)) { + } else if (selectionPIDHe3(t2) && selectionPIDProton(t1)) { he3Cand = t2, protonCand = t1; - passPID = true; - } - if (!passPID) { + } else { continue; } bool heliumPID = he3Cand.pidForTracking() == o2::track::PID::Helium3 || he3Cand.pidForTracking() == o2::track::PID::Alpha; - float correctedTPCinnerParam = (heliumPID && cfgCompensatePIDinTracking) ? he3Cand.tpcInnerParam() / 2.f : he3Cand.tpcInnerParam(); - histos.fill(HIST("h2dEdxHe3candidates"), correctedTPCinnerParam * 2.f, he3Cand.tpcSignal()); + float correctedTPCinnerParam = (heliumPID && setting_compensatePIDinTracking) ? he3Cand.tpcInnerParam() / 2.f : he3Cand.tpcInnerParam(); + m_qaRegistry.fill(HIST("h2dEdxHe3candidates"), correctedTPCinnerParam * 2.f, he3Cand.tpcSignal()); - if (!FillCandidateInfo(he3Cand, protonCand, true)) { - continue; - } - // fill TOF info outside to avoide responseBeta crash - auto& cand = l4Candidates.back(); - if (he3Cand.hasTOF()) { - float beta = responseBeta.GetBeta(he3Cand); - beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked - - bool heliumPID = t1.pidForTracking() == o2::track::PID::Helium3 || he3Cand.pidForTracking() == o2::track::PID::Alpha; - float correctedTPCinnerParamHe3 = (heliumPID && cfgCompensatePIDinTracking) ? he3Cand.tpcInnerParam() / 2.f : he3Cand.tpcInnerParam(); - cand.massTOFHe3 = correctedTPCinnerParamHe3 * 2.f * std::sqrt(1.f / (beta * beta) - 1.f); + SVCand trackPair; + trackPair.tr0Idx = he3Cand.globalIndex(); + trackPair.tr1Idx = protonCand.globalIndex(); + m_trackPairs.push_back(trackPair); + } + } + } + + void fillTable(const Lithium4Candidate& li4cand, bool isMC = false) + { + if (!isMC) { + m_outputDataTable( + li4cand.recoPtHe3(), + li4cand.recoEtaHe3(), + li4cand.recoPhiHe3(), + li4cand.recoPtPr(), + li4cand.recoEtaPr(), + li4cand.recoPhiPr(), + li4cand.DCAxyHe3, + li4cand.DCAzHe3, + li4cand.DCAxyPr, + li4cand.DCAzPr, + li4cand.tpcSignalHe3, + li4cand.momHe3TPC, + li4cand.tpcSignalPr, + li4cand.momPrTPC, + li4cand.nTPCClustersHe3, + li4cand.nSigmaHe3, + li4cand.nSigmaPr, + li4cand.chi2TPCHe3, + li4cand.chi2TPCPr, + li4cand.massTOFHe3, + li4cand.massTOFPr, + li4cand.PIDtrkHe3, + li4cand.PIDtrkPr, + li4cand.itsClSizeHe3, + li4cand.itsClSizePr, + li4cand.sharedClustersHe3, + li4cand.sharedClustersPr, + li4cand.isBkgUS, + li4cand.isBkgEM); + } else { + m_outputMCTable( + li4cand.recoPtHe3(), + li4cand.recoEtaHe3(), + li4cand.recoPhiHe3(), + li4cand.recoPtPr(), + li4cand.recoEtaPr(), + li4cand.recoPhiPr(), + li4cand.DCAxyHe3, + li4cand.DCAzHe3, + li4cand.DCAxyPr, + li4cand.DCAzPr, + li4cand.tpcSignalHe3, + li4cand.momHe3TPC, + li4cand.tpcSignalPr, + li4cand.momPrTPC, + li4cand.nTPCClustersHe3, + li4cand.nSigmaHe3, + li4cand.nSigmaPr, + li4cand.chi2TPCHe3, + li4cand.chi2TPCPr, + li4cand.massTOFHe3, + li4cand.massTOFPr, + li4cand.PIDtrkHe3, + li4cand.PIDtrkPr, + li4cand.itsClSizeHe3, + li4cand.itsClSizePr, + li4cand.sharedClustersHe3, + li4cand.sharedClustersPr, + li4cand.isBkgUS, + li4cand.isBkgEM, + li4cand.momHe3MC, + li4cand.etaHe3MC, + li4cand.phiHe3MC, + li4cand.momPrMC, + li4cand.etaPrMC, + li4cand.phiPrMC, + li4cand.l4PtMC, + li4cand.l4MassMC); + } + } + + void fillHistograms(const Lithium4Candidate& li4cand) + { + m_qaRegistry.fill(HIST("hHe3Pt"), li4cand.recoPtHe3()); + m_qaRegistry.fill(HIST("hProtonPt"), li4cand.recoPtPr()); + m_qaRegistry.fill(HIST("hLitInvMass"), li4cand.invMass); + m_qaRegistry.fill(HIST("hDCAxyHe3"), li4cand.DCAxyHe3); + m_qaRegistry.fill(HIST("hDCAzHe3"), li4cand.DCAzHe3); + } + + // ================================================================================================================== + + void processSameEvent(const CollisionsFull& collisions, const TrackCandidates& tracks, const aod::BCsWithTimestamps&) + { + for (auto& collision : collisions) { + + m_trackPairs.clear(); + m_qaRegistry.fill(HIST("hEvents"), 0); + + if (!collision.sel8() || std::abs(collision.posZ()) > setting_cutVertex) { + continue; + } + if (setting_skimmedProcessing) { + bool zorroSelected = m_zorro.isSelected(collision.template bc_as().globalBC()); + if (zorroSelected) { + m_qaRegistry.fill(HIST("hEvents"), 2); } - if (protonCand.hasTOF()) { - float beta = responseBeta.GetBeta(protonCand); - beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked - cand.massTOFPr = protonCand.tpcInnerParam() * std::sqrt(1.f / (beta * beta) - 1.f); + } + m_qaRegistry.fill(HIST("hEvents"), 1); + m_qaRegistry.fill(HIST("hNcontributor"), collision.numContrib()); + m_qaRegistry.fill(HIST("hVtxZ"), collision.posZ()); + + const uint64_t collIdx = collision.globalIndex(); + auto TrackTable_thisCollision = tracks.sliceBy(m_perCol, collIdx); + TrackTable_thisCollision.bindExternalIndices(&tracks); + + pairTracksSameEvent(TrackTable_thisCollision); + + for (auto& trackPair : m_trackPairs) { + + auto heTrack = tracks.rawIteratorAt(trackPair.tr0Idx); + auto prTrack = tracks.rawIteratorAt(trackPair.tr1Idx); + + Lithium4Candidate li4cand; + if (!fillCandidateInfo(heTrack, prTrack, li4cand, /*mix*/ false)) { + continue; } - fillHistograms(cand); + fillHistograms(li4cand); + fillTable(li4cand, false); } } + } + PROCESS_SWITCH(lithium4analysis, processSameEvent, "Process Same event", false); + + void processMixedEvent(const CollisionsFull& /*collisions*/, const TrackCandidates& tracks) + { + m_trackPairs.clear(); + pairTracksEventMixing(); - for (auto& l4Cand : l4Candidates) { - outputDataTable(l4Cand.recoPtHe3(), l4Cand.recoEtaHe3(), l4Cand.recoPhiHe3(), - l4Cand.recoPtPr(), l4Cand.recoEtaPr(), l4Cand.recoPhiPr(), - l4Cand.DCAxyHe3, l4Cand.DCAzHe3, l4Cand.DCAxyPr, l4Cand.DCAzPr, - l4Cand.tpcSignalHe3, l4Cand.momHe3TPC, l4Cand.tpcSignalPr, l4Cand.momPrTPC, - l4Cand.nTPCClustersHe3, - l4Cand.nSigmaHe3, l4Cand.nSigmaPr, l4Cand.massTOFHe3, l4Cand.massTOFPr, - l4Cand.PIDtrkHe3, l4Cand.PIDtrkPr, l4Cand.itsClSizeHe3, l4Cand.itsClSizePr, - l4Cand.sharedClustersHe3, l4Cand.sharedClustersPr, - l4Cand.isBkgUS, l4Cand.isBkgEM); + for (auto& trackPair : m_trackPairs) { + + auto heTrack = tracks.rawIteratorAt(trackPair.tr0Idx); + auto prTrack = tracks.rawIteratorAt(trackPair.tr1Idx); + + Lithium4Candidate li4cand; + if (!fillCandidateInfo(heTrack, prTrack, li4cand, /*mix*/ true)) { + continue; + } + fillHistograms(li4cand); + fillTable(li4cand, false); } } PROCESS_SWITCH(lithium4analysis, processMixedEvent, "Process Mixed event", false); - void processMC(soa::Join const& collisions, aod::BCs const&, TrackCandidatesMC const& tracks, aod::McParticles const& mcParticles) + void processMC(const CollisionsFullMC& collisions, const aod::BCsWithTimestamps&, const TrackCandidatesMC& tracks, const aod::McParticles& mcParticles) { std::vector filledMothers; - l4Candidates.clear(); for (auto& collision : collisions) { - if (!collision.sel8() || std::abs(collision.posZ()) > cfgCutVertex) { + m_trackPairs.clear(); + m_qaRegistry.fill(HIST("hEvents"), 0); + if (/*!collision.sel8() ||*/ std::abs(collision.posZ()) > setting_cutVertex) { continue; } - histos.fill(HIST("hNcontributor"), collision.numContrib()); - histos.fill(HIST("hVtxZ"), collision.posZ()); + m_qaRegistry.fill(HIST("hEvents"), 1); + m_qaRegistry.fill(HIST("hNcontributor"), collision.numContrib()); + m_qaRegistry.fill(HIST("hVtxZ"), collision.posZ()); const uint64_t collIdx = collision.globalIndex(); - auto TrackTable_thisCollision = tracks.sliceBy(perColMC, collIdx); + auto TrackTable_thisCollision = tracks.sliceBy(m_perColMC, collIdx); TrackTable_thisCollision.bindExternalIndices(&tracks); - for (auto track1 : TrackTable_thisCollision) { + pairTracksSameEvent(TrackTable_thisCollision); - if (!track1.has_mcParticle()) { - continue; - } + for (auto& trackPair : m_trackPairs) { - histos.fill(HIST("hTrackSel"), Selections::kNoCuts); + auto heTrack = tracks.rawIteratorAt(trackPair.tr0Idx); + auto prTrack = tracks.rawIteratorAt(trackPair.tr1Idx); - if (!track1.isGlobalTrackWoDCA()) { + if (!heTrack.has_mcParticle() || !prTrack.has_mcParticle()) { continue; } - histos.fill(HIST("hTrackSel"), Selections::kGlobalTrack); - if (!selectionTrack(track1)) { - continue; - } - histos.fill(HIST("hTrackSel"), Selections::kTrackCuts); + auto mctrackHe3 = heTrack.mcParticle(); + auto mctrackPr = prTrack.mcParticle(); - if (!selectionPIDHe3(track1)) { + if (std::abs(mctrackHe3.pdgCode()) != hePDG || std::abs(mctrackPr.pdgCode()) != prPDG) { continue; } - histos.fill(HIST("hTrackSel"), Selections::kPID); - for (auto track2 : TrackTable_thisCollision) { - if (!track2.has_mcParticle()) { - continue; - } - - if (!track2.isGlobalTrackWoDCA()) { - continue; - } - - if (!selectionTrack(track2)) { - continue; - } - - if (!selectionPIDProton(track2)) { - continue; - } - - if (track1.sign() * track2.sign() < 0) { - continue; - } - - const auto mctrackHe3 = track1.mcParticle(); - const auto mctrackPr = track2.mcParticle(); + for (auto& mothertrack : mctrackHe3.mothers_as()) { + for (auto& mothertrackPr : mctrackPr.mothers_as()) { - if (std::abs(mctrackHe3.pdgCode()) != he3PDG || std::abs(mctrackPr.pdgCode()) != protonPDG) { - continue; - } + if (mothertrack != mothertrackPr || std::abs(mothertrack.pdgCode()) != li4PDG || std::abs(mothertrack.y()) > 1) { + continue; + } - for (auto& mothertrack : mctrackHe3.mothers_as()) { - for (auto& mothertrackPr : mctrackPr.mothers_as()) { - - if (mothertrack != mothertrackPr || std::abs(mothertrack.pdgCode()) != lithium4PDG) { - continue; - } - - if (std::abs(mothertrack.y()) > 1) { - continue; - } - - if (!FillCandidateInfo(track1, track2, false, true)) { - continue; - } - - // fill TOF info outside to avoide responseBeta crash - auto& cand = l4Candidates.back(); - if (track1.hasTOF()) { - float beta = responseBetaMC.GetBeta(track1); - beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked - bool heliumPID = track1.pidForTracking() == o2::track::PID::Helium3 || track1.pidForTracking() == o2::track::PID::Alpha; - float correctedTPCinnerParamHe3 = (heliumPID && cfgCompensatePIDinTracking) ? track1.tpcInnerParam() / 2.f : track1.tpcInnerParam(); - cand.massTOFHe3 = correctedTPCinnerParamHe3 * 2.f * std::sqrt(1.f / (beta * beta) - 1.f); - } - if (track2.hasTOF()) { - float beta = responseBetaMC.GetBeta(track2); - beta = std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked - cand.massTOFPr = track2.tpcInnerParam() * std::sqrt(1.f / (beta * beta) - 1.f); - } - - cand.momHe3MC = mctrackHe3.pt() * (mctrackHe3.pdgCode() > 0 ? 1 : -1); - cand.momPrMC = mctrackPr.pt() * (mctrackPr.pdgCode() > 0 ? 1 : -1); - cand.l4PtMC = mothertrack.pt() * (mothertrack.pdgCode() > 0 ? 1 : -1); - double eLit = mctrackHe3.e() + mctrackPr.e(); - cand.l4MassMC = std::sqrt(eLit * eLit - mothertrack.p() * mothertrack.p()); - filledMothers.push_back(mothertrack.globalIndex()); - fillHistograms(cand); + Lithium4Candidate li4cand; + if (!fillCandidateInfo(heTrack, prTrack, li4cand, /*mix*/ false)) { + continue; } + fillCandidateInfoMC(mctrackHe3, mctrackPr, mothertrack, li4cand); + fillHistograms(li4cand); + fillTable(li4cand, true); + filledMothers.push_back(mothertrack.globalIndex()); } } } @@ -615,11 +712,7 @@ struct lithium4analysis { for (auto& mcParticle : mcParticles) { - if (std::abs(mcParticle.pdgCode()) != lithium4PDG) { - continue; - } - - if (std::abs(mcParticle.y()) > 1 || mcParticle.isPhysicalPrimary() == false) { + if (std::abs(mcParticle.pdgCode()) != li4PDG || std::abs(mcParticle.y()) > 1 || mcParticle.isPhysicalPrimary() == false) { continue; } @@ -628,53 +721,28 @@ struct lithium4analysis { } auto kDaughters = mcParticle.daughters_as(); - auto daughtHe3 = false; - auto daughtPr = false; - double eLit = 0; - int signHe3 = 0, signPr = 0; - double ptHe3 = 0, ptPr = 0; + bool daughtHe3(false), daughtPr(false); + McIter mcHe3, mcPr; for (auto kCurrentDaughter : kDaughters) { - if (std::abs(kCurrentDaughter.pdgCode()) == he3PDG) { + if (std::abs(kCurrentDaughter.pdgCode()) == hePDG) { daughtHe3 = true; - signHe3 = kCurrentDaughter.pdgCode() > 0 ? 1 : -1; - ptHe3 = kCurrentDaughter.pt(); - eLit += kCurrentDaughter.e(); - } else if (std::abs(kCurrentDaughter.pdgCode()) == protonPDG) { + mcHe3 = kCurrentDaughter; + } else if (std::abs(kCurrentDaughter.pdgCode()) == prPDG) { daughtPr = true; - signPr = kCurrentDaughter.pdgCode() > 0 ? 1 : -1; - ptPr = kCurrentDaughter.pt(); - eLit += kCurrentDaughter.e(); + mcPr = kCurrentDaughter; } } if (daughtHe3 && daughtPr) { - lithium4Candidate l4Candidate; - int signLi = mcParticle.pdgCode() > 0 ? 1 : -1; - l4Candidate.l4PtMC = mcParticle.pt() * signLi; - l4Candidate.momHe3MC = ptHe3 * signHe3; - l4Candidate.momPrMC = ptPr * signPr; - l4Candidate.l4MassMC = std::sqrt(eLit * eLit - mcParticle.p() * mcParticle.p()); - l4Candidates.push_back(l4Candidate); + Lithium4Candidate li4cand; + fillCandidateInfoMC(mcHe3, mcPr, mcParticle, li4cand); + fillTable(li4cand, true); } } - - for (auto& l4Cand : l4Candidates) { - outputMCTable(l4Cand.recoPtHe3(), l4Cand.recoEtaHe3(), l4Cand.recoPhiHe3(), - l4Cand.recoPtPr(), l4Cand.recoEtaPr(), l4Cand.recoPhiPr(), - l4Cand.DCAxyHe3, l4Cand.DCAzHe3, l4Cand.DCAxyPr, l4Cand.DCAzPr, - l4Cand.tpcSignalHe3, l4Cand.momHe3TPC, l4Cand.tpcSignalPr, l4Cand.momPrTPC, - l4Cand.nTPCClustersHe3, - l4Cand.nSigmaHe3, l4Cand.nSigmaPr, l4Cand.massTOFHe3, l4Cand.massTOFPr, - l4Cand.PIDtrkHe3, l4Cand.PIDtrkPr, l4Cand.itsClSizeHe3, l4Cand.itsClSizePr, - l4Cand.sharedClustersHe3, l4Cand.sharedClustersPr, - l4Cand.isBkgUS, l4Cand.isBkgEM, - l4Cand.momHe3MC, l4Cand.momPrMC, - l4Cand.l4PtMC, l4Cand.l4MassMC); - } } PROCESS_SWITCH(lithium4analysis, processMC, "Process MC", false); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(const ConfigContext& cfgc) { return WorkflowSpec{ adaptAnalysisTask(cfgc, TaskName{"lithium4analysis"})}; From 694c5d162e677b6218e638e5b604f3f311294e9a Mon Sep 17 00:00:00 2001 From: prottayCMT <61418725+prottayCMT@users.noreply.github.com> Date: Fri, 13 Sep 2024 03:32:52 +0530 Subject: [PATCH 05/14] added vz axis for rec (#7672) --- PWGLF/TableProducer/Common/spvector.cxx | 49 ++++++++++++++++--------- 1 file changed, 31 insertions(+), 18 deletions(-) diff --git a/PWGLF/TableProducer/Common/spvector.cxx b/PWGLF/TableProducer/Common/spvector.cxx index 9ce30f60ad8..d242909b92a 100644 --- a/PWGLF/TableProducer/Common/spvector.cxx +++ b/PWGLF/TableProducer/Common/spvector.cxx @@ -152,7 +152,7 @@ struct spvector { AxisSpec channelZDCAxis = {8, 0.0, 8.0, "ZDC tower"}; AxisSpec qxZDCAxis = {QxyNbins, lbinQxy, hbinQxy, "Qx"}; AxisSpec qyZDCAxis = {QxyNbins, lbinQxy, hbinQxy, "Qy"}; - AxisSpec phiAxis = {50, -6.28, 6.28, "phi"}; + AxisSpec phiAxis = {100, -6.28, 6.28, "phi"}; AxisSpec vzAxis = {20, -10, 10, "vz"}; histos.add("hCentrality", "hCentrality", kTH1F, {{8, 0, 80.0}}); @@ -162,10 +162,10 @@ struct spvector { histos.add("hpQyZDCAC", "hpQyZDCAC", kTProfile, {centAxis}); histos.add("hpQxZDCAQyZDCC", "hpQxZDCAQyZDCC", kTProfile, {centAxis}); histos.add("hpQxZDCCQyZDCA", "hpQxZDCCQyZDCA", kTProfile, {centAxis}); - histos.add("QxZDCC", "QxZDCC", kTH2F, {centAxis, qxZDCAxis}); - histos.add("QyZDCC", "QyZDCC", kTH2F, {centAxis, qyZDCAxis}); - histos.add("QxZDCA", "QxZDCA", kTH2F, {centAxis, qxZDCAxis}); - histos.add("QyZDCA", "QyZDCA", kTH2F, {centAxis, qyZDCAxis}); + histos.add("QxZDCC", "QxZDCC", kTH3F, {centAxis, vzAxis, qxZDCAxis}); + histos.add("QyZDCC", "QyZDCC", kTH3F, {centAxis, vzAxis, qyZDCAxis}); + histos.add("QxZDCA", "QxZDCA", kTH3F, {centAxis, vzAxis, qxZDCAxis}); + histos.add("QyZDCA", "QyZDCA", kTH3F, {centAxis, vzAxis, qyZDCAxis}); histos.add("PsiZDCC", "PsiZDCC", kTH2F, {centAxis, phiAxis}); histos.add("PsiZDCA", "PsiZDCA", kTH2F, {centAxis, phiAxis}); histos.add("ZDCAmp", "ZDCAmp", kTProfile2D, {channelZDCAxis, vzAxis}); @@ -193,7 +193,7 @@ struct spvector { int currentRunNumber = -999; int lastRunNumber = -999; TH2D* gainprofile; - TH2D* hrecentere; + TH3D* hrecentere; // Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPT); // Filter DCAcutFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz); @@ -258,7 +258,10 @@ struct spvector { for (std::size_t iChA = 0; iChA < 8; iChA++) { auto chanelid = iChA; if (useGainCallib && gainprofile) { - gainequal = gainprofile->GetBinContent(gainprofile->FindBin(chanelid, vz)); + gainequal = gainprofile->GetBinContent(gainprofile->FindBin(vz, chanelid + 0.5)); + // if (chanelid==0) + // LOG(info) <<"###################Accesed#####################"<<" "<getForTimeStamp(ConfRecentere.value, bc.timestamp()); + hrecentere = ccdb->getForTimeStamp(ConfRecentere.value, bc.timestamp()); } - if (useRecentere) { - - qxZDCA = (qxZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, 0.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, 0.5)); - qyZDCA = (qyZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, 1.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, 1.5)); - qxZDCC = (qxZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, 2.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, 2.5)); - qyZDCC = (qyZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, 3.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, 3.5)); + if (useRecentere && hrecentere) { + /* + qxZDCA = (qxZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 0.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 0.5)); + qyZDCA = (qyZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 1.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 1.5)); + qxZDCC = (qxZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 2.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 2.5)); + qyZDCC = (qyZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 3.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 3.5)); + */ + + qxZDCA = (qxZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 0.5))); + qyZDCA = (qyZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 1.5))); + qxZDCC = (qxZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 2.5))); + qyZDCC = (qyZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 3.5))); + /* + if (vz>0.0) + LOG(info) <<"###################Accesed#####################"<<" "<GetBinContent(hrecentere->FindBin(centrality, vz, 0.5)); + */ } psiZDCC = 1.0 * TMath::ATan2(qyZDCC, qxZDCC); @@ -324,10 +337,10 @@ struct spvector { histos.fill(HIST("hpQyZDCAC"), centrality, (qyZDCA * qyZDCC)); histos.fill(HIST("hpQxZDCAQyZDCC"), centrality, (qxZDCA * qyZDCC)); histos.fill(HIST("hpQxZDCCQyZDCA"), centrality, (qxZDCC * qyZDCA)); - histos.fill(HIST("QxZDCC"), centrality, qxZDCC); - histos.fill(HIST("QyZDCC"), centrality, qyZDCC); - histos.fill(HIST("QxZDCA"), centrality, qxZDCA); - histos.fill(HIST("QyZDCA"), centrality, qyZDCA); + histos.fill(HIST("QxZDCC"), centrality, vz, qxZDCC); + histos.fill(HIST("QyZDCC"), centrality, vz, qyZDCC); + histos.fill(HIST("QxZDCA"), centrality, vz, qxZDCA); + histos.fill(HIST("QyZDCA"), centrality, vz, qyZDCA); histos.fill(HIST("PsiZDCA"), centrality, psiZDCA); histos.fill(HIST("PsiZDCC"), centrality, psiZDCC); From 10d3a6025d613e32ff1870c17efd286b2b62f762 Mon Sep 17 00:00:00 2001 From: lucamicheletti93 <38209984+lucamicheletti93@users.noreply.github.com> Date: Fri, 13 Sep 2024 06:59:02 +0200 Subject: [PATCH 06/14] Adding control histograms for zorro (#7674) --- PWGDQ/TableProducer/tableMaker.cxx | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/PWGDQ/TableProducer/tableMaker.cxx b/PWGDQ/TableProducer/tableMaker.cxx index 93871082aa5..e619d809572 100644 --- a/PWGDQ/TableProducer/tableMaker.cxx +++ b/PWGDQ/TableProducer/tableMaker.cxx @@ -56,6 +56,7 @@ #include "DetectorsBase/Propagator.h" #include "DetectorsBase/GeometryManager.h" #include "EventFiltering/Zorro.h" +#include "Framework/HistogramRegistry.h" using std::cout; using std::endl; @@ -172,13 +173,18 @@ struct TableMaker { Configurable fConfigDetailedQA{"cfgDetailedQA", false, "If true, include more QA histograms (BeforeCuts classes)"}; Configurable fIsRun2{"cfgIsRun2", false, "Whether we analyze Run-2 or Run-3 data"}; Configurable fIsAmbiguous{"cfgIsAmbiguous", false, "Whether we enable QA plots for ambiguous tracks"}; - Configurable fConfigRunZorro{"cfgRunZorro", false, "Enable event selection with zorro [WARNING: under debug, do not enable!]"}; - Configurable fConfigZorroTrigMask{"cfgZorroTriggerMask", "fDiMuon", "DQ Trigger masks: fSingleE,fLMeeIMR,fLMeeHMR,fDiElectron,fSingleMuLow,fSingleMuHigh,fDiMuon"}; + + struct : ConfigurableGroup { + Configurable fConfigRunZorro{"cfgRunZorro", false, "Enable event selection with zorro [WARNING: under debug, do not enable!]"}; + Configurable fConfigZorroTrigMask{"cfgZorroTriggerMask", "fDiMuon", "DQ Trigger masks: fSingleE,fLMeeIMR,fLMeeHMR,fDiElectron,fSingleMuLow,fSingleMuHigh,fDiMuon"}; + } useZorro; + struct : ConfigurableGroup { Configurable fConfigCcdbUrl{"useCCDBConfigurations.ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable fConfigCcdbPathTPC{"useCCDBConfigurations.ccdb-path-tpc", "Users/z/zhxiong/TPCPID/PostCalib", "base path to the ccdb object"}; Configurable fConfigCcdbPathZorro{"useCCDBConfigurations.ccdb-path-zorro", "Users/r/rlietava/EventFiltering/OTS/", "base path to the ccdb object for zorro"}; } useCCDBConfigurations; + Configurable fConfigNoLaterThan{"ccdb-no-later-than", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; Configurable fConfigComputeTPCpostCalib{"cfgTPCpostCalib", false, "If true, compute TPC post-calibrated n-sigmas(electrons, pions, protons)"}; Configurable fConfigComputeTPCpostCalibKaon{"cfgTPCpostCalibKaon", false, "If true, compute TPC post-calibrated n-sigmas for kaons"}; @@ -212,6 +218,8 @@ struct TableMaker { bool fDoDetailedQA = false; // Bool to set detailed QA true, if QA is set true int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. + HistogramRegistry registry{"registry"}; + // TODO: filter on TPC dedx used temporarily until electron PID will be improved Filter barrelSelectedTracks = ifnode(fIsRun2.node() == true, aod::track::trackType == uint8_t(aod::track::Run2Track), aod::track::trackType == uint8_t(aod::track::Track)) && o2::aod::track::pt >= fConfigBarrelTrackPtLow && nabs(o2::aod::track::eta) <= fConfigBarrelTrackMaxAbsEta && o2::aod::track::tpcSignal >= fConfigMinTpcSignal && o2::aod::track::tpcSignal <= fConfigMaxTpcSignal && o2::aod::track::tpcChi2NCl < 4.0f && o2::aod::track::itsChi2NCl < 36.0f; @@ -406,7 +414,7 @@ struct TableMaker { } // Put the 8 first bits of the event filter in the last 8 bits of the tag if constexpr ((TEventFillMap & VarManager::ObjTypes::EventFilter) > 0) { - if (!fConfigRunZorro) { + if (!useZorro.fConfigRunZorro) { tag |= (collision.eventFilter() << 56); } } @@ -432,9 +440,12 @@ struct TableMaker { } (reinterpret_cast(fStatsList->At(0)))->Fill(2.0, static_cast(kNaliases)); - if (fConfigRunZorro) { + if (useZorro.fConfigRunZorro) { zorro.setBaseCCDBPath(useCCDBConfigurations.fConfigCcdbPathZorro.value); - zorro.initCCDB(fCCDB.service, fCurrentRun, bc.timestamp(), fConfigZorroTrigMask.value); + zorro.initCCDB(fCCDB.service, fCurrentRun, bc.timestamp(), useZorro.fConfigZorroTrigMask.value); + + zorro.populateHistRegistry(registry, fCurrentRun); + if (zorro.isSelected(bc.globalBC())) { tag |= (static_cast(true) << 56); // the same bit is used for this zorro selections from ccdb } @@ -861,7 +872,7 @@ struct TableMaker { } // Put the 8 first bits of the event filter in the last 8 bits of the tag if constexpr ((TEventFillMap & VarManager::ObjTypes::EventFilter) > 0) { - if (!fConfigRunZorro) { + if (!useZorro.fConfigRunZorro) { tag |= (collision.eventFilter() << 56); } } @@ -886,9 +897,12 @@ struct TableMaker { } (reinterpret_cast(fStatsList->At(0)))->Fill(2.0, static_cast(kNaliases)); - if (fConfigRunZorro) { + if (useZorro.fConfigRunZorro) { zorro.setBaseCCDBPath(useCCDBConfigurations.fConfigCcdbPathZorro.value); - zorro.initCCDB(fCCDB.service, fCurrentRun, bc.timestamp(), fConfigZorroTrigMask.value); + zorro.initCCDB(fCCDB.service, fCurrentRun, bc.timestamp(), useZorro.fConfigZorroTrigMask.value); + + zorro.populateHistRegistry(registry, fCurrentRun); + if (zorro.isSelected(bc.globalBC())) { tag |= (static_cast(true) << 56); // the same bit is used for this zorro selections from ccdb } From aa69c6b82040ba63670f171e3e488899c479290b Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi <43742195+fmazzasc@users.noreply.github.com> Date: Fri, 13 Sep 2024 06:59:52 +0200 Subject: [PATCH 07/14] Restore rigidity cut for He3 (#7667) --- PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx b/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx index a66622bb5cb..b3eec29372b 100644 --- a/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx @@ -387,6 +387,8 @@ struct hyperRecoTask { hypCand.clusterSizeITSPi = piTrack.itsClusterSizes(); bool heliumPID = heTrack.pidForTracking() == o2::track::PID::Helium3 || heTrack.pidForTracking() == o2::track::PID::Alpha; hypCand.momHe3TPC = (heliumPID && cfgCompensatePIDinTracking) ? heTrack.tpcInnerParam() / 2 : heTrack.tpcInnerParam(); + if (hypCand.momHe3TPC < TPCRigidityMinHe) + return; hypCand.momPiTPC = piTrack.tpcInnerParam(); hDeDxTot->Fill(hypCand.momHe3TPC * heTrack.sign(), heTrack.tpcSignal()); hDeDxTot->Fill(hypCand.momPiTPC * piTrack.sign(), piTrack.tpcSignal()); From d8d2f6d57595106905ba9001fa08c1ea19b66af2 Mon Sep 17 00:00:00 2001 From: Giorgio Alberto Lucia <87222843+GiorgioAlbertoLucia@users.noreply.github.com> Date: Fri, 13 Sep 2024 15:16:55 +0200 Subject: [PATCH 08/14] [PWGLF] Fixed zorro, removed unnecessary prints (#7677) --- PWGLF/TableProducer/Nuspex/lithium4analysis.cxx | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx b/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx index 39ba20d9cfd..6ad6afd931a 100644 --- a/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx +++ b/PWGLF/TableProducer/Nuspex/lithium4analysis.cxx @@ -113,7 +113,6 @@ struct Lithium4Candidate { float nSigmaPr = -10.f; float chi2TPCHe3 = -10.f; float chi2TPCPr = -10.f; - float massTOFHe3 = -10; float massTOFPr = -10; @@ -122,7 +121,6 @@ struct Lithium4Candidate { uint32_t itsClSizeHe3 = 0u; uint32_t itsClSizePr = 0u; - uint8_t sharedClustersHe3 = 0u; uint8_t sharedClustersPr = 0u; @@ -232,17 +230,14 @@ struct lithium4analysis { void init(o2::framework::InitContext&) { - LOG(info) << "Initializing lithium4 analysis"; m_zorroSummary.setObject(m_zorro.getZorroSummary()); m_runNumber = 0; - LOG(info) << "Initializing CCDB"; m_ccdb->setURL(setting_ccdburl); m_ccdb->setCaching(true); m_ccdb->setLocalObjectValidityChecking(); m_ccdb->setFatalWhenNull(false); - LOG(info) << "Initializing PID response"; for (int i = 0; i < 5; i++) { m_BBparamsHe[i] = setting_BetheBlochParams->get("He3", Form("p%i", i)); } @@ -261,16 +256,16 @@ struct lithium4analysis { void initCCDB(const aod::BCsWithTimestamps::iterator& bc) { - LOG(info) << "Initializing CCDB for run " << bc.runNumber(); + if (m_runNumber == bc.runNumber()) { return; } - LOG(info) << "Initializing zorro for run " << bc.runNumber(); + if (setting_skimmedProcessing) { m_zorro.initCCDB(m_ccdb.service, bc.runNumber(), bc.timestamp(), "fHe"); m_zorro.populateHistRegistry(m_qaRegistry, bc.runNumber()); } - LOG(info) << "Initializing CCDB for run " << bc.runNumber() << " done"; + m_runNumber = bc.runNumber(); } @@ -596,6 +591,8 @@ struct lithium4analysis { m_trackPairs.clear(); m_qaRegistry.fill(HIST("hEvents"), 0); + auto bc = collision.template bc_as(); + initCCDB(bc); if (!collision.sel8() || std::abs(collision.posZ()) > setting_cutVertex) { continue; @@ -638,7 +635,6 @@ struct lithium4analysis { pairTracksEventMixing(); for (auto& trackPair : m_trackPairs) { - auto heTrack = tracks.rawIteratorAt(trackPair.tr0Idx); auto prTrack = tracks.rawIteratorAt(trackPair.tr1Idx); From 32741546b54b6baa30d77a6a56c981c9b67d233e Mon Sep 17 00:00:00 2001 From: sarjeetagami <162087855+sarjeetagami@users.noreply.github.com> Date: Sat, 14 Sep 2024 00:12:35 +0530 Subject: [PATCH 09/14] fake kaon rejection (#7683) Co-authored-by: sarjeeta gami --- PWGLF/Tasks/Resonances/kstarpbpb.cxx | 97 ++++++++++++++++++++++------ 1 file changed, 77 insertions(+), 20 deletions(-) diff --git a/PWGLF/Tasks/Resonances/kstarpbpb.cxx b/PWGLF/Tasks/Resonances/kstarpbpb.cxx index 4894d8149d4..9eb18cec5dc 100644 --- a/PWGLF/Tasks/Resonances/kstarpbpb.cxx +++ b/PWGLF/Tasks/Resonances/kstarpbpb.cxx @@ -73,8 +73,10 @@ struct kstarpbpb { // events Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable cfgCutCentrality{"cfgCutCentrality", 80.0f, "Accepted maximum Centrality"}; // track Configurable cfgCutCharge{"cfgCutCharge", 0.0, "cut on Charge"}; + Configurable additionalEvSel2{"additionalEvSel2", true, "Additional evsel2"}; Configurable cfgCutPT{"cfgCutPT", 0.2, "PT cut on daughter track"}; Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut on daughter track"}; Configurable cfgCutDCAxy{"cfgCutDCAxy", 2.0f, "DCAxy range for tracks"}; @@ -89,6 +91,8 @@ struct kstarpbpb { ConfigurableAxis configThnAxisInvMass{"configThnAxisInvMass", {180, 0.6, 1.5}, "#it{M} (GeV/#it{c}^{2})"}; ConfigurableAxis configThnAxisPt{"configThnAxisPt", {100, 0.0, 10.}, "#it{p}_{T} (GeV/#it{c})"}; ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {8, 0., 80}, "Centrality"}; + Configurable removefaketrak{"removefaketrack", true, "Remove fake track from momentum difference"}; + Configurable ConfFakeKaonCut{"ConfFakeKaonCut", 0.1, "Cut based on track from momentum difference"}; ConfigurableAxis configThnAxisPhiminusPsi{"configThnAxisPhiminusPsi", {6, 0.0, TMath::Pi()}, "#phi - #psi"}; ConfigurableAxis configThnAxisV2{"configThnAxisV2", {200, -1, 1}, "V2"}; Configurable additionalEvsel{"additionalEvsel", false, "Additional event selcection"}; @@ -99,8 +103,10 @@ struct kstarpbpb { Configurable confMaxRot{"confMaxRot", 7.0 * TMath::Pi() / 6.0, "Maximum of rotation"}; Configurable nBkgRotations{"nBkgRotations", 9, "Number of rotated copies (background) per each original candidate"}; Configurable fillRotation{"fillRotation", true, "fill rotation"}; + Configurable cfgOccupancyCut{"cfgOccupancyCut", 500, "Occupancy cut"}; Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter centralityFilter = nabs(aod::cent::centFT0C) < cfgCutCentrality; Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPT); Filter DCAcutFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz); @@ -204,7 +210,7 @@ struct kstarpbpb { double massPi = o2::constants::physics::MassPiMinus; template - bool eventSelected(TCollision collision, const int& /*multTrk*/, const float& centrality) + bool eventSelected(TCollision collision, const float& centrality) { if (collision.alias_bit(kTVXinTRD)) { // TRD triggered @@ -224,7 +230,6 @@ struct kstarpbpb { return 1; } - template bool selectionTrack(const T& candidate) { @@ -241,17 +246,17 @@ struct kstarpbpb { bool selectionPIDpTdependent(const T& candidate, int PID) { if (PID == 0) { - if (candidate.p() < 0.5 && TMath::Abs(candidate.tpcNSigmaKa()) < nsigmaCutTPC) { + if (candidate.pt() < 0.5 && TMath::Abs(candidate.tpcNSigmaKa()) < nsigmaCutTPC) { return true; } - if (candidate.p() >= 0.5 && candidate.hasTOF() && ((candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) + (candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa())) < (nsigmaCutCombined * nsigmaCutCombined)) { + if (candidate.pt() >= 0.5 && candidate.hasTOF() && ((candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) + (candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa())) < (nsigmaCutCombined * nsigmaCutCombined)) { return true; } } else if (PID == 1) { - if (candidate.p() < 0.5 && TMath::Abs(candidate.tpcNSigmaPi()) < nsigmaCutTPC) { + if (candidate.pt() < 0.5 && TMath::Abs(candidate.tpcNSigmaPi()) < nsigmaCutTPC) { return true; } - if (candidate.p() >= 0.5 && candidate.hasTOF() && ((candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) + (candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi())) < (nsigmaCutCombined * nsigmaCutCombined)) { + if (candidate.pt() >= 0.5 && candidate.hasTOF() && ((candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) + (candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi())) < (nsigmaCutCombined * nsigmaCutCombined)) { return true; } } @@ -262,10 +267,16 @@ struct kstarpbpb { bool selectionPID(const T& candidate, int PID) { if (PID == 0) { + if (!candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaKa()) < nsigmaCutTPC) { + return true; + } if (candidate.hasTOF() && ((candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) + (candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa())) < (nsigmaCutCombined * nsigmaCutCombined)) { return true; } } else if (PID == 1) { + if (!candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPi()) < nsigmaCutTPC) { + return true; + } if (candidate.hasTOF() && ((candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) + (candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi())) < (nsigmaCutCombined * nsigmaCutCombined)) { return true; } @@ -277,17 +288,17 @@ struct kstarpbpb { bool selectionPIDNew(const T& candidate, int PID) { if (PID == 0) { - if (candidate.p() < 0.6 && TMath::Abs(candidate.tpcNSigmaKa()) < 2.0) { + if (candidate.pt() < 0.6 && TMath::Abs(candidate.tpcNSigmaKa()) < 2.0) { return true; } - if (candidate.p() >= 0.6 && candidate.p() < 3.0 && candidate.hasTOF() && candidate.tpcNSigmaKa() > -2.0 && candidate.tpcNSigmaKa() < 3.0 && TMath::Abs(candidate.tofNSigmaKa()) < 2.0) { + if (candidate.pt() >= 0.6 && candidate.pt() < 3.0 && candidate.hasTOF() && candidate.tpcNSigmaKa() > -2.0 && candidate.tpcNSigmaKa() < 3.0 && TMath::Abs(candidate.tofNSigmaKa()) < 2.0) { return true; } } else if (PID == 1) { - if (candidate.p() < 1.0 && TMath::Abs(candidate.tpcNSigmaPi()) < 2.0) { + if (candidate.pt() < 1.0 && TMath::Abs(candidate.tpcNSigmaPi()) < 2.0) { return true; } - if (candidate.p() >= 1.0 && candidate.p() < 3.0 && candidate.hasTOF() && candidate.tpcNSigmaPi() > -2.0 && candidate.tpcNSigmaPi() < 3.0 && TMath::Abs(candidate.tofNSigmaPi()) < 2.0) { + if (candidate.pt() >= 1.0 && candidate.pt() < 3.0 && candidate.hasTOF() && candidate.tpcNSigmaPi() > -2.0 && candidate.tpcNSigmaPi() < 3.0 && TMath::Abs(candidate.tofNSigmaPi()) < 2.0) { return true; } } @@ -305,12 +316,22 @@ struct kstarpbpb { } return result; } - + template + bool isFakeKaon(T const& track) + { + const auto pglobal = track.p(); + const auto ptpc = track.tpcInnerParam(); + if (TMath::Abs(pglobal - ptpc) > ConfFakeKaonCut) { + return true; + } + return false; + } ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for bin"}; ConfigurableAxis axisMultiplicityClass{"axisMultiplicityClass", {20, 0, 100}, "multiplicity percentile for bin"}; ConfigurableAxis axisEPAngle{"axisEPAngle", {6, -TMath::Pi() / 2, TMath::Pi() / 2}, "event plane angle"}; + ConfigurableAxis axisOccup{"axisOccup", {20, 0.0, 40000.0}, "occupancy axis"}; - using BinningTypeVertexContributor = ColumnBinningPolicy; + using BinningTypeVertexContributor = ColumnBinningPolicy; ROOT::Math::PxPyPzMVector KstarMother, daughter1, daughter2, kaonrot, kstarrot; void processSameEvent(EventCandidates::iterator const& collision, TrackCandidates const& tracks, aod::BCs const&) @@ -319,6 +340,9 @@ struct kstarpbpb { if (!collision.sel8()) { return; } + if (additionalEvSel2 && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) { + return; + } histos.fill(HIST("hMC"), 1.5); if (!collision.triggereventep()) { return; @@ -336,12 +360,12 @@ struct kstarpbpb { auto psiFT0A = collision.psiFT0A(); auto psiTPC = collision.psiTPC(); int occupancy = collision.trackOccupancyInTimeRange(); - if (occupancy >= 1500) // occupancy info is available for this collision (*) + if (occupancy >= cfgOccupancyCut) // occupancy info is available for this collision (*) { return; } histos.fill(HIST("hFTOCvsTPC"), centrality, multTPC); - if (additionalEvsel && !eventSelected(collision, tracks.size(), centrality)) { + if (additionalEvsel && !eventSelected(collision, centrality)) { return; } histos.fill(HIST("hFTOCvsTPCSelected"), centrality, multTPC); @@ -402,20 +426,32 @@ struct kstarpbpb { if (selectionPIDNew(track1, 1) && selectionPIDNew(track2, 0)) { track1pion = true; track2kaon = true; + if (removefaketrak && isFakeKaon(track2)) { + continue; + } } if (selectionPIDNew(track2, 1) && selectionPIDNew(track1, 0)) { track2pion = true; track1kaon = true; + if (removefaketrak && isFakeKaon(track1)) { + continue; + } } } if (!ispTdepPID) { if (selectionPID(track1, 1) && selectionPID(track2, 0)) { track1pion = true; track2kaon = true; + if (removefaketrak && isFakeKaon(track2)) { + continue; + } } if (selectionPID(track2, 1) && selectionPID(track1, 0)) { track2pion = true; track1kaon = true; + if (removefaketrak && isFakeKaon(track1)) { + continue; + } } } @@ -483,7 +519,7 @@ struct kstarpbpb { void processMixedEvent(EventCandidates const& collisions, TrackCandidates const& tracks) { auto tracksTuple = std::make_tuple(tracks); - BinningTypeVertexContributor binningOnPositions{{axisVertex, axisMultiplicityClass, axisEPAngle}, true}; + BinningTypeVertexContributor binningOnPositions{{axisVertex, axisMultiplicityClass, axisOccup}, true}; SameKindPair pair{binningOnPositions, cfgNoMixedEvents, -1, collisions, tracksTuple, &cache}; for (auto& [collision1, tracks1, collision2, tracks2] : pair) { if (!collision1.sel8() || !collision2.sel8()) { @@ -492,12 +528,19 @@ struct kstarpbpb { if (!collision1.triggereventep() || !collision2.triggereventep()) { continue; } + if (additionalEvSel2 && (!collision1.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision1.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) { + continue; + } + if (additionalEvSel2 && (!collision2.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision2.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) { + continue; + } if (timFrameEvsel && (!collision1.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision2.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision1.selection_bit(aod::evsel::kNoITSROFrameBorder) || !collision2.selection_bit(aod::evsel::kNoITSROFrameBorder))) { continue; } - int occupancy = collision1.trackOccupancyInTimeRange(); - if (occupancy >= 1500) { - return; + int occupancy1 = collision1.trackOccupancyInTimeRange(); + int occupancy2 = collision2.trackOccupancyInTimeRange(); + if (occupancy1 >= cfgOccupancyCut || occupancy2 >= cfgOccupancyCut) { + continue; } auto centrality = collision1.centFT0C(); auto centrality2 = collision2.centFT0C(); @@ -507,10 +550,12 @@ struct kstarpbpb { bool track2pion = false; bool track2kaon = false; - if (additionalEvsel && !eventSelected(collision1, tracks.size(), centrality)) { + if (additionalEvsel && !eventSelected(collision1, centrality)) { + // printf("Mix = %d\n", 4); continue; } - if (additionalEvsel && !eventSelected(collision2, tracks.size(), centrality2)) { + if (additionalEvsel && !eventSelected(collision2, centrality2)) { + // printf("Mix = %d\n", 5); continue; } @@ -538,20 +583,32 @@ struct kstarpbpb { if (selectionPIDNew(track1, 1) && selectionPIDNew(track2, 0)) { track1pion = true; track2kaon = true; + if (removefaketrak && isFakeKaon(track2)) { + continue; + } } if (selectionPIDNew(track2, 1) && selectionPIDNew(track1, 0)) { track2pion = true; track1kaon = true; + if (removefaketrak && isFakeKaon(track1)) { + continue; + } } } if (!ispTdepPID) { if (selectionPID(track1, 1) && selectionPID(track2, 0)) { track1pion = true; track2kaon = true; + if (removefaketrak && isFakeKaon(track2)) { + continue; + } } if (selectionPID(track2, 1) && selectionPID(track1, 0)) { track2pion = true; track1kaon = true; + if (removefaketrak && isFakeKaon(track1)) { + continue; + } } } if (track1kaon && track2pion) { From ef8ad848efeff0f0ccd32841f6f6d6610d292960 Mon Sep 17 00:00:00 2001 From: prottayCMT <61418725+prottayCMT@users.noreply.github.com> Date: Sat, 14 Sep 2024 00:13:08 +0530 Subject: [PATCH 10/14] updated rec with its error (#7680) Co-authored-by: Prottay Das --- PWGLF/TableProducer/Common/spvector.cxx | 31 +++++++++++++------------ 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/PWGLF/TableProducer/Common/spvector.cxx b/PWGLF/TableProducer/Common/spvector.cxx index d242909b92a..2a218f651bf 100644 --- a/PWGLF/TableProducer/Common/spvector.cxx +++ b/PWGLF/TableProducer/Common/spvector.cxx @@ -101,6 +101,7 @@ struct spvector { Configurable useGainCallib{"useGainCallib", false, "use gain calibration"}; Configurable useRecentere{"useRecentere", false, "use Recentering"}; + Configurable recwitherror{"recwitherror", false, "use Recentering with error"}; Configurable useShift{"useShift", false, "use Shift"}; Configurable ConfGainPath{"ConfGainPath", "Users/p/prottay/My/Object/NewPbPbpass4_10092024/gaincallib", "Path to gain calibration"}; Configurable ConfRecentere{"ConfRecentere", "Users/p/prottay/My/Object/NewPbPbpass4_23082024/recenter", "Path for recentere"}; @@ -313,21 +314,21 @@ struct spvector { } if (useRecentere && hrecentere) { - /* - qxZDCA = (qxZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 0.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 0.5)); - qyZDCA = (qyZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 1.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 1.5)); - qxZDCC = (qxZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 2.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 2.5)); - qyZDCC = (qyZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 3.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 3.5)); - */ - - qxZDCA = (qxZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 0.5))); - qyZDCA = (qyZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 1.5))); - qxZDCC = (qxZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 2.5))); - qyZDCC = (qyZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 3.5))); - /* - if (vz>0.0) - LOG(info) <<"###################Accesed#####################"<<" "<GetBinContent(hrecentere->FindBin(centrality, vz, 0.5)); - */ + if (recwitherror) { + qxZDCA = (qxZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 0.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 0.5)); + qyZDCA = (qyZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 1.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 1.5)); + qxZDCC = (qxZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 2.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 2.5)); + qyZDCC = (qyZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 3.5))) / hrecentere->GetBinError(hrecentere->FindBin(centrality, vz, 3.5)); + } else { + + qxZDCA = (qxZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 0.5))); + qyZDCA = (qyZDCA - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 1.5))); + qxZDCC = (qxZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 2.5))); + qyZDCC = (qyZDCC - hrecentere->GetBinContent(hrecentere->FindBin(centrality, vz, 3.5))); + } + + // if (vz>0.0) + // LOG(info) <<"###################Accesed#####################"<<" "<GetBinContent(hrecentere->FindBin(centrality, vz, 0.5)); } psiZDCC = 1.0 * TMath::ATan2(qyZDCC, qxZDCC); From e6f21acd95a2258e391935af7fd8c2340be94172 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sat, 14 Sep 2024 01:14:46 +0200 Subject: [PATCH 11/14] Common: save only non-empty BCs, add BC flags creator (#7678) * Common: save only non-empty BCs, add BC flags creator * Please consider the following formatting changes (#344) * Update multiplicityExtraTable.cxx --------- Co-authored-by: ALICE Builder --- .../TableProducer/Converters/CMakeLists.txt | 5 +++ .../Converters/bcFlagsCreator.cxx | 36 +++++++++++++++++++ .../TableProducer/multiplicityExtraTable.cxx | 17 +++++++++ 3 files changed, 58 insertions(+) create mode 100644 Common/TableProducer/Converters/bcFlagsCreator.cxx diff --git a/Common/TableProducer/Converters/CMakeLists.txt b/Common/TableProducer/Converters/CMakeLists.txt index 989b58bdab8..e7b76fd1d42 100644 --- a/Common/TableProducer/Converters/CMakeLists.txt +++ b/Common/TableProducer/Converters/CMakeLists.txt @@ -49,6 +49,11 @@ o2physics_add_dpl_workflow(bc-converter PUBLIC_LINK_LIBRARIES O2::Framework COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(bc-flags-creator + SOURCES bcFlagsCreator.cxx + PUBLIC_LINK_LIBRARIES O2::Framework + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(calo-label-converter SOURCES caloLabelConverter.cxx PUBLIC_LINK_LIBRARIES diff --git a/Common/TableProducer/Converters/bcFlagsCreator.cxx b/Common/TableProducer/Converters/bcFlagsCreator.cxx new file mode 100644 index 00000000000..8ac9d6e1516 --- /dev/null +++ b/Common/TableProducer/Converters/bcFlagsCreator.cxx @@ -0,0 +1,36 @@ +// 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. +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" + +using namespace o2; +using namespace o2::framework; + +// Creates an empty BCFlags for data that doesn't have it to be used seamlessly +// n.b. this will overwrite existing BCFlags, to be discussed if data in mixed condition +struct bcFlagsCreator { + Produces bcFlags; + + void process(aod::BCs const& bcTable) + { + for (auto& bc : bcTable) { + bcFlags(0); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + }; +} diff --git a/Common/TableProducer/multiplicityExtraTable.cxx b/Common/TableProducer/multiplicityExtraTable.cxx index 7a77faa6784..6a6bc1a76a8 100644 --- a/Common/TableProducer/multiplicityExtraTable.cxx +++ b/Common/TableProducer/multiplicityExtraTable.cxx @@ -39,6 +39,7 @@ struct MultiplicityExtraTable { // Allow for downscaling of BC table for less space use in derived data Configurable bcDownscaleFactor{"bcDownscaleFactor", 2, "Downscale factor for BC table (0: save nothing, 1: save all)"}; Configurable minFT0CforBCTable{"minFT0CforBCTable", 25.0f, "Minimum FT0C amplitude to fill BC table to reduce data"}; + Configurable saveOnlyBCsWithCollisions{"saveOnlyBCsWithCollisions", true, "save only BCs with collisions in them"}; // needed for downscale unsigned int randomSeed = 0; @@ -66,13 +67,24 @@ struct MultiplicityExtraTable { { //+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+ // determine saved BCs and corresponding new BC table index + std::vector bcHasCollision(bcs.size()); std::vector newBCindex(bcs.size()); std::vector bc2multArray(bcs.size()); int atIndex = 0; for (const auto& bc : bcs) { + bcHasCollision[bc.globalIndex()] = false; newBCindex[bc.globalIndex()] = -1; bc2multArray[bc.globalIndex()] = -1; + } + + //+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+ + // tag BCs that have a collision (from evsel foundBC) + for (const auto& collision : collisions) { + bcHasCollision[collision.foundBCId()] = true; + } + //+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+ + for (const auto& bc : bcs) { // downscale if requested to do so if (bcDownscaleFactor < 1.f && (static_cast(rand_r(&randomSeed)) / static_cast(RAND_MAX)) > bcDownscaleFactor) { continue; @@ -91,6 +103,11 @@ struct MultiplicityExtraTable { if (multFT0C < minFT0CforBCTable) { continue; // skip this event } + + if (saveOnlyBCsWithCollisions && !bcHasCollision[bc.globalIndex()]) { + continue; // skip if no collision is assigned to this BC (from evSel assignment) + } + newBCindex[bc.globalIndex()] = atIndex++; } //+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+ From e5b1083d31e88761c05ca606cb5b0112722133f6 Mon Sep 17 00:00:00 2001 From: Ionut Cristian Arsene Date: Sat, 14 Sep 2024 02:39:12 +0200 Subject: [PATCH 12/14] [PWG-DQ] added an inv mass histogram with fine pt binning in the jpsi photoproduction range (#7685) * solve conflicts * clang fixes * removed duplicate function * Added a inv mass histogram with fine binning in the photoproduction pt range * clang format --------- Co-authored-by: Ionut Cristian Arsene --- PWGDQ/Core/HistogramsLibrary.cxx | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index b0f5a155762..303dbe3a6d8 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -812,12 +812,24 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "Mass_HighRange", "", false, 375, 0.0, 15.0, VarManager::kMass); hm->AddHistogram(histClass, "Pt", "", false, 2000, 0.0, 20., VarManager::kPt); hm->AddHistogram(histClass, "Mass_Pt", "", false, 125, 0.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt); + double massBins[76]; + for (int i = 0; i < 76; i++) { + massBins[i] = 1.5 + i * 0.04; + } + double ptBins[70]; + for (int i = 0; i <= 50; i++) { + ptBins[i] = i * 0.01; + } + for (int i = 1; i <= 19; i++) { + ptBins[50 + i] = 0.5 + i * 0.5; + } + hm->AddHistogram(histClass, "Mass_PtFine", "", false, 75, massBins, VarManager::kMass, 69, ptBins, VarManager::kPt); hm->AddHistogram(histClass, "Eta_Pt", "", false, 40, -2.0, 2.0, VarManager::kEta, 40, 0.0, 20.0, VarManager::kPt); hm->AddHistogram(histClass, "Mass_VtxZ", "", true, 30, -15.0, 15.0, VarManager::kVtxZ, 500, 0.0, 5.0, VarManager::kMass); if (subGroupStr.Contains("pbpb")) { hm->AddHistogram(histClass, "Mass_CentFT0C", "", false, 125, 0.0, 5.0, VarManager::kMass, 20, 0.0, 100.0, VarManager::kCentFT0C); hm->AddHistogram(histClass, "Pt_CentFT0C", "", false, 100, 0.0, 10.0, VarManager::kPt, 20, 0.0, 100.0, VarManager::kCentFT0C); - hm->AddHistogram(histClass, "Mass_Pt_CentFT0C", "", false, 125, 0.0, 5.0, VarManager::kMass, 100, 0.0, 10.0, VarManager::kPt, 10, 0.0, 100.0, VarManager::kCentFT0C); + hm->AddHistogram(histClass, "Mass_Pt_CentFT0C", "", false, 75, 1.5, 4.5, VarManager::kMass, 20, 0.0, 10.0, VarManager::kPt, 10, 0.0, 100.0, VarManager::kCentFT0C); } if (subGroupStr.Contains("mult")) { hm->AddHistogram(histClass, "Mass_Pt_MultFV0A", "", false, 200, 0.0, 5.0, VarManager::kMass, 40, 0.0, 40.0, VarManager::kPt, 100, 0.0, 25000.0, VarManager::kMultFV0A); From 224965ca88bd2949b0085b464e70395439fd92b8 Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Sat, 14 Sep 2024 03:46:03 +0200 Subject: [PATCH 13/14] PWGEM/PhotonMeson: add N TPC shared clusters in data table (#7673) --- PWGEM/Dilepton/Core/PhotonHBT.h | 156 ++++++++++++++---- PWGEM/Dilepton/Tasks/vpPairQC.cxx | 4 +- PWGEM/Dilepton/Utils/EMTrackUtilities.h | 4 +- PWGEM/PhotonMeson/DataModel/gammaTables.h | 17 +- .../TableProducer/photonconversionbuilder.cxx | 7 +- .../TableProducer/skimmerGammaConversion.cxx | 6 +- PWGEM/PhotonMeson/Tasks/pcmQC.cxx | 19 ++- PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx | 19 ++- PWGEM/PhotonMeson/Utils/PCMUtilities.h | 76 +++++++++ 9 files changed, 253 insertions(+), 55 deletions(-) diff --git a/PWGEM/Dilepton/Core/PhotonHBT.h b/PWGEM/Dilepton/Core/PhotonHBT.h index 7d8958d2372..529e62c8ffe 100644 --- a/PWGEM/Dilepton/Core/PhotonHBT.h +++ b/PWGEM/Dilepton/Core/PhotonHBT.h @@ -417,10 +417,9 @@ struct PhotonHBT { } } - if constexpr (pairtype == ggHBTPairType::kPCMPCM) { // dr, dz of conversion points - fRegistry.add("Pair/same/hDeltaRDeltaZ", "diphoton distance in RZ;#Deltar = #sqrt{(#Deltax)^{2} + (#Deltay)^{2}} (cm);#Deltaz (cm)", kTH2D, {{100, 0, +10}, {200, -10, 10}}, true); - } else { // deta, dphi of track momentum - fRegistry.add("Pair/same/hDeltaEtaDeltaPhi", "diphoton distance in #eta-#varphi plane;#Delta#varphi (rad.);#Delta#eta", kTH2D, {{200, -0.1, +0.1}, {200, -0.1, 0.1}}, true); + fRegistry.add("Pair/same/hDeltaEtaDeltaPhi", "distance between 2 LS tracks in #eta-#varphi plane;#Delta#varphi (rad.);#Delta#eta", kTH2D, {{200, -0.1, +0.1}, {200, -0.1, 0.1}}, true); // deta, dphi of track momentum + if constexpr (pairtype == ggHBTPairType::kPCMPCM) { // dr, dz of conversion points + fRegistry.add("Pair/same/hDeltaRDeltaZ", "diphoton distance in RZ;#Deltar = #sqrt{(#Deltax)^{2} + (#Deltay)^{2}} (cm);#Deltaz (cm)", kTH2D, {{100, 0, +10}, {100, 0, 10}}, true); } fRegistry.addClone("Pair/same/", "Pair/mix/"); @@ -708,10 +707,26 @@ struct PhotonHBT { float dz = g1.vz() - g2.vz(); float dr = std::sqrt(std::pow(g1.vx() - g2.vx(), 2) + std::pow(g1.vy() - g2.vy(), 2)); - if (ggpaircuts.applydR && std::pow(dz / ggpaircuts.cfgMinDeltaZ, 2) + std::pow(dr / ggpaircuts.cfgMinDeltaR, 2) < 1.f) { + // if (ggpaircuts.applydR && std::pow(dz / ggpaircuts.cfgMinDeltaZ, 2) + std::pow(dr / ggpaircuts.cfgMinDeltaR, 2) < 1.f) { + // continue; + // } + + float deta_pos = pos1.eta() - pos2.eta(); + float dphi_pos = pos1.phi() - pos2.phi(); + o2::math_utils::bringToPMPi(dphi_pos); + float deta_ele = ele1.eta() - ele2.eta(); + float dphi_ele = ele1.phi() - ele2.phi(); + o2::math_utils::bringToPMPi(dphi_ele); + if (ggpaircuts.applydR && std::pow(deta_pos / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_pos / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + continue; + } + if (ggpaircuts.applydR && std::pow(deta_ele / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_ele / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { continue; } - fRegistry.fill(HIST("Pair/same/hDeltaRDeltaZ"), dr, dz, 1.f); + + fRegistry.fill(HIST("Pair/same/hDeltaRDeltaZ"), dr, fabs(dz), 1.f); + fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), pos1.phi() - pos2.phi(), pos1.eta() - pos2.eta(), 1.f); // distance between 2 LS tracks + fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), ele1.phi() - ele2.phi(), ele1.eta() - ele2.eta(), 1.f); // distance between 2 LS tracks fillPairHistogram<0>(collision, v1, v2, 1.f); ndiphoton++; @@ -721,12 +736,16 @@ struct PhotonHBT { if (std::find(used_photonIds.begin(), used_photonIds.end(), pair_tmp_id1) == used_photonIds.end()) { EMTrack g1tmp = EMTrack(ndf, g1.globalIndex(), collision.globalIndex(), g1.globalIndex(), g1.pt(), g1.eta(), g1.phi(), 0); g1tmp.setConversionPointXYZ(g1.vx(), g1.vy(), g1.vz()); + g1tmp.setPositiveLegPtEtaPhiM(pos1.pt(), pos1.eta(), pos1.phi(), o2::constants::physics::MassElectron); + g1tmp.setNegativeLegPtEtaPhiM(ele1.pt(), ele1.eta(), ele1.phi(), o2::constants::physics::MassElectron); emh1->AddTrackToEventPool(key_df_collision, g1tmp); used_photonIds.emplace_back(pair_tmp_id1); } if (std::find(used_photonIds.begin(), used_photonIds.end(), pair_tmp_id2) == used_photonIds.end()) { EMTrack g2tmp = EMTrack(ndf, g2.globalIndex(), collision.globalIndex(), g2.globalIndex(), g2.pt(), g2.eta(), g2.phi(), 0); g2tmp.setConversionPointXYZ(g2.vx(), g2.vy(), g2.vz()); + g2tmp.setPositiveLegPtEtaPhiM(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron); + g2tmp.setNegativeLegPtEtaPhiM(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron); emh1->AddTrackToEventPool(key_df_collision, g2tmp); used_photonIds.emplace_back(pair_tmp_id2); } @@ -801,11 +820,11 @@ struct PhotonHBT { std::pair pair_tmp = std::make_pair(std::make_pair(pos1.trackId(), ele1.trackId()), std::make_pair(pos2.trackId(), ele2.trackId())); if (std::find(used_pairs_per_collision.begin(), used_pairs_per_collision.end(), pair_tmp) == used_pairs_per_collision.end()) { - float deta_pos = v_pos1.Eta() - v_pos2.Eta(); - float dphi_pos = v_pos1.Phi() - v_pos2.Phi(); + float deta_pos = pos1.eta() - pos2.eta(); + float dphi_pos = pos1.phi() - pos2.phi(); o2::math_utils::bringToPMPi(dphi_pos); - float deta_ele = v_ele1.Eta() - v_ele2.Eta(); - float dphi_ele = v_ele1.Phi() - v_ele2.Phi(); + float deta_ele = ele1.eta() - ele2.eta(); + float dphi_ele = ele1.phi() - ele2.phi(); o2::math_utils::bringToPMPi(dphi_ele); if (ggpaircuts.applydR && std::pow(deta_pos / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_pos / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { continue; @@ -814,8 +833,8 @@ struct PhotonHBT { continue; } - fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), v_pos1.Phi() - v_pos2.Phi(), v_pos1.Eta() - v_pos2.Eta(), weight1 * weight2); // distance between 2 LS tracks - fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), v_ele1.Phi() - v_ele2.Phi(), v_ele1.Eta() - v_ele2.Eta(), weight1 * weight2); // distance between 2 LS tracks + fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), dphi_pos, deta_pos, weight1 * weight2); // distance between 2 LS tracks + fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), dphi_ele, deta_ele, weight1 * weight2); // distance between 2 LS tracks fillPairHistogram<0>(collision, v1_ee, v2_ee, weight1 * weight2); ndiphoton++; used_pairs_per_collision.emplace_back(std::make_pair(pair_tmp.first, pair_tmp.second)); @@ -833,8 +852,8 @@ struct PhotonHBT { g1pair.setGlobalPosId(pos1.globalIndex()); g1pair.setGlobalNegId(ele1.globalIndex()); g1pair.setPairDca3DinSigmaOTF(dca1_3d); - g1pair.setPositiveLegPtEtaPhiM(v_pos1.Pt(), v_pos1.Eta(), v_pos1.Phi(), o2::constants::physics::MassElectron); - g1pair.setNegativeLegPtEtaPhiM(v_ele1.Pt(), v_ele1.Eta(), v_ele1.Phi(), o2::constants::physics::MassElectron); + g1pair.setPositiveLegPtEtaPhiM(pos1.pt(), pos1.eta(), pos1.phi(), o2::constants::physics::MassElectron); + g1pair.setNegativeLegPtEtaPhiM(ele1.pt(), ele1.eta(), ele1.phi(), o2::constants::physics::MassElectron); g1pair.setAmbPosLegSelfIds(possibleIds_pos1); g1pair.setAmbNegLegSelfIds(possibleIds_ele1); emh1->AddTrackToEventPool(key_df_collision, g1pair); @@ -850,8 +869,8 @@ struct PhotonHBT { g2pair.setGlobalPosId(pos2.globalIndex()); g2pair.setGlobalNegId(ele2.globalIndex()); g2pair.setPairDca3DinSigmaOTF(dca2_3d); - g2pair.setPositiveLegPtEtaPhiM(v_pos2.Pt(), v_pos2.Eta(), v_pos2.Phi(), o2::constants::physics::MassElectron); - g2pair.setNegativeLegPtEtaPhiM(v_ele2.Pt(), v_ele2.Eta(), v_ele2.Phi(), o2::constants::physics::MassElectron); + g2pair.setPositiveLegPtEtaPhiM(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron); + g2pair.setNegativeLegPtEtaPhiM(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron); g2pair.setAmbPosLegSelfIds(possibleIds_pos2); g2pair.setAmbNegLegSelfIds(possibleIds_ele2); emh1->AddTrackToEventPool(key_df_collision, g2pair); @@ -909,25 +928,40 @@ struct PhotonHBT { ROOT::Math::PtEtaPhiMVector v_pos2(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron); ROOT::Math::PtEtaPhiMVector v_ele2(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron); ROOT::Math::PtEtaPhiMVector v2_ee = v_pos2 + v_ele2; - float deta = v1_gamma.Eta() - v2_ee.Eta(); - float dphi = v1_gamma.Phi() - v2_ee.Phi(); - o2::math_utils::bringToPMPi(dphi); - if (ggpaircuts.applydR && std::pow(deta / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + + float deta_pos = pos1.eta() - pos2.eta(); + float dphi_pos = pos1.phi() - pos2.phi(); + o2::math_utils::bringToPMPi(dphi_pos); + float deta_ele = ele1.eta() - ele2.eta(); + float dphi_ele = ele1.phi() - ele2.phi(); + o2::math_utils::bringToPMPi(dphi_ele); + if (ggpaircuts.applydR && std::pow(deta_pos / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_pos / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + continue; + } + if (ggpaircuts.applydR && std::pow(deta_ele / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_ele / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { continue; } - fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), dphi, deta, weight); + + fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), dphi_pos, deta_pos, weight); + fRegistry.fill(HIST("Pair/same/hDeltaEtaDeltaPhi"), dphi_ele, deta_ele, weight); fillPairHistogram<0>(collision, v1_gamma, v2_ee, weight); ndiphoton++; std::pair pair_tmp_id1 = std::make_pair(ndf, g1.globalIndex()); std::tuple tuple_tmp_id2 = std::make_tuple(ndf, collision.globalIndex(), pos2.globalIndex(), ele2.globalIndex()); if (std::find(used_photonIds.begin(), used_photonIds.end(), pair_tmp_id1) == used_photonIds.end()) { - emh1->AddTrackToEventPool(key_df_collision, EMTrack(ndf, g1.globalIndex(), collision.globalIndex(), g1.globalIndex(), g1.pt(), g1.eta(), g1.phi(), 0)); + EMTrack g1tmp = EMTrack(ndf, g1.globalIndex(), collision.globalIndex(), g1.globalIndex(), g1.pt(), g1.eta(), g1.phi(), 0); + g1tmp.setConversionPointXYZ(g1.vx(), g1.vy(), g1.vz()); + g1tmp.setPositiveLegPtEtaPhiM(pos1.pt(), pos1.eta(), pos1.phi(), o2::constants::physics::MassElectron); + g1tmp.setNegativeLegPtEtaPhiM(ele1.pt(), ele1.eta(), ele1.phi(), o2::constants::physics::MassElectron); + emh1->AddTrackToEventPool(key_df_collision, g1tmp); used_photonIds.emplace_back(pair_tmp_id1); } if (std::find(used_dileptonIds.begin(), used_dileptonIds.end(), tuple_tmp_id2) == used_dileptonIds.end()) { EMTrack g2pair = EMTrack(ndf, -1, collision.globalIndex(), -1, v2_ee.Pt(), v2_ee.Eta(), v2_ee.Phi(), v2_ee.M()); g2pair.setPairDca3DinSigmaOTF(dca2_3d); + g2pair.setPositiveLegPtEtaPhiM(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron); + g2pair.setNegativeLegPtEtaPhiM(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron); emh2->AddTrackToEventPool(key_df_collision, g2pair); used_dileptonIds.emplace_back(tuple_tmp_id2); } @@ -971,12 +1005,33 @@ struct PhotonHBT { ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); + auto pos1 = g1.getPositiveLeg(); + auto ele1 = g1.getNegativeLeg(); + auto pos2 = g2.getPositiveLeg(); + auto ele2 = g2.getNegativeLeg(); + float dz = g1.vz() - g2.vz(); float dr = std::sqrt(std::pow(g1.vx() - g2.vx(), 2) + std::pow(g1.vy() - g2.vy(), 2)); - if (ggpaircuts.applydR && std::pow(dz / ggpaircuts.cfgMinDeltaZ, 2) + std::pow(dr / ggpaircuts.cfgMinDeltaR, 2) < 1.f) { + // if (ggpaircuts.applydR && std::pow(dz / ggpaircuts.cfgMinDeltaZ, 2) + std::pow(dr / ggpaircuts.cfgMinDeltaR, 2) < 1.f) { + // continue; + // } + + float deta_pos = pos1.Eta() - pos2.Eta(); + float dphi_pos = pos1.Phi() - pos2.Phi(); + o2::math_utils::bringToPMPi(dphi_pos); + float deta_ele = ele1.Eta() - ele2.Eta(); + float dphi_ele = ele1.Phi() - ele2.Phi(); + o2::math_utils::bringToPMPi(dphi_ele); + if (ggpaircuts.applydR && std::pow(deta_pos / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_pos / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + continue; + } + if (ggpaircuts.applydR && std::pow(deta_ele / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_ele / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { continue; } - fRegistry.fill(HIST("Pair/mix/hDeltaRDeltaZ"), dr, dz, 1.f); + + fRegistry.fill(HIST("Pair/mix/hDeltaRDeltaZ"), dr, fabs(dz), 1.f); + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_pos, deta_pos, 1.f); // distance between 2 LS tracks + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_ele, deta_ele, 1.f); // distance between 2 LS tracks fillPairHistogram<1>(collision, v1, v2, 1.f); } } @@ -1024,6 +1079,7 @@ struct PhotonHBT { ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), g1.mass()); ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), g2.mass()); + float deta_pos = pos1.Eta() - pos2.Eta(); float dphi_pos = pos1.Phi() - pos2.Phi(); o2::math_utils::bringToPMPi(dphi_pos); @@ -1036,8 +1092,8 @@ struct PhotonHBT { if (ggpaircuts.applydR && std::pow(deta_ele / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_ele / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { continue; } - fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), pos1.Phi() - pos2.Phi(), pos1.Eta() - pos2.Eta(), 1.f); - fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), ele1.Phi() - ele2.Phi(), ele1.Eta() - ele2.Eta(), 1.f); + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_pos, deta_pos, 1.f); // distance between 2 LS tracks + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_ele, deta_ele, 1.f); // distance between 2 LS tracks fillPairHistogram<1>(collision, v1, v2, 1.f); } } @@ -1064,13 +1120,27 @@ struct PhotonHBT { for (auto& g2 : photons2_from_event_pool) { // dielectron ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.0); // keep v1 for PCM ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), g2.mass()); - float deta = v1.Eta() - v2.Eta(); - float dphi = v1.Phi() - v2.Phi(); - o2::math_utils::bringToPMPi(dphi); - if (ggpaircuts.applydR && std::pow(deta / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + + auto pos1 = g1.getPositiveLeg(); + auto ele1 = g1.getNegativeLeg(); + auto pos2 = g2.getPositiveLeg(); + auto ele2 = g2.getNegativeLeg(); + + float deta_pos = pos1.Eta() - pos2.Eta(); + float dphi_pos = pos1.Phi() - pos2.Phi(); + o2::math_utils::bringToPMPi(dphi_pos); + float deta_ele = ele1.Eta() - ele2.Eta(); + float dphi_ele = ele1.Phi() - ele2.Phi(); + o2::math_utils::bringToPMPi(dphi_ele); + if (ggpaircuts.applydR && std::pow(deta_pos / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_pos / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { continue; } - fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi, deta, 1.f); + if (ggpaircuts.applydR && std::pow(deta_ele / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_ele / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + continue; + } + + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_pos, deta_pos, 1.f); // distance between 2 LS tracks + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_ele, deta_ele, 1.f); // distance between 2 LS tracks fillPairHistogram<1>(collision, v1, v2, 1.f); } } @@ -1097,13 +1167,27 @@ struct PhotonHBT { for (auto& g2 : photons1_from_event_pool) { // PCM ROOT::Math::PtEtaPhiMVector v1(g2.pt(), g2.eta(), g2.phi(), 0.0); // keep v1 for PCM ROOT::Math::PtEtaPhiMVector v2(g1.pt(), g1.eta(), g1.phi(), g1.mass()); - float deta = v1.Eta() - v2.Eta(); - float dphi = v1.Phi() - v2.Phi(); - o2::math_utils::bringToPMPi(dphi); - if (ggpaircuts.applydR && std::pow(deta / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + + auto pos1 = g1.getPositiveLeg(); + auto ele1 = g1.getNegativeLeg(); + auto pos2 = g2.getPositiveLeg(); + auto ele2 = g2.getNegativeLeg(); + + float deta_pos = pos1.Eta() - pos2.Eta(); + float dphi_pos = pos1.Phi() - pos2.Phi(); + o2::math_utils::bringToPMPi(dphi_pos); + float deta_ele = ele1.Eta() - ele2.Eta(); + float dphi_ele = ele1.Phi() - ele2.Phi(); + o2::math_utils::bringToPMPi(dphi_ele); + if (ggpaircuts.applydR && std::pow(deta_pos / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_pos / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { continue; } - fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi, deta, 1.f); + if (ggpaircuts.applydR && std::pow(deta_ele / ggpaircuts.cfgMinDeltaEta, 2) + std::pow(dphi_ele / ggpaircuts.cfgMinDeltaPhi, 2) < 1.f) { + continue; + } + + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_pos, deta_pos, 1.f); // distance between 2 LS tracks + fRegistry.fill(HIST("Pair/mix/hDeltaEtaDeltaPhi"), dphi_ele, deta_ele, 1.f); // distance between 2 LS tracks fillPairHistogram<1>(collision, v1, v2, 1.f); } } diff --git a/PWGEM/Dilepton/Tasks/vpPairQC.cxx b/PWGEM/Dilepton/Tasks/vpPairQC.cxx index 6f6f23bb45b..dd0d91e8abb 100644 --- a/PWGEM/Dilepton/Tasks/vpPairQC.cxx +++ b/PWGEM/Dilepton/Tasks/vpPairQC.cxx @@ -215,8 +215,8 @@ struct vpPairQC { fRegistry.add("Track/positive/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); fRegistry.add("Track/positive/hPResolution", "p resolution;p (GeV/c);#Deltap/p", kTH2F, {{1000, 0.0f, 10.0f}, {100, 0.0f, 0.1f}}, false); fRegistry.add("Track/positive/hPtResolution", "p_{T} resolution;p (GeV/c);#Deltap_{T}/p_{T}", kTH2F, {{1000, 0.0f, 10.0f}, {100, 0.0f, 0.1f}}, false); - fRegistry.add("Track/positive/hThetaResolution", "#theta resolution;p (GeV/c);#Delta#theta (rad.)", kTH2F, {{1000, 0.0f, 10.0f}, {200, 0.0f, 0.02f}}, false); - fRegistry.add("Track/positive/hPhiResolution", "#varphi resolution;p (GeV/c);#Delta#varphi (rad.)", kTH2F, {{1000, 0.0f, 10.0f}, {200, 0.0f, 0.02f}}, false); + fRegistry.add("Track/positive/hThetaResolution", "#theta resolution;p (GeV/c);#Delta#theta (rad.)", kTH2F, {{1000, 0.0f, 10.0f}, {100, 0.0f, 0.01f}}, false); + fRegistry.add("Track/positive/hPhiResolution", "#varphi resolution;p (GeV/c);#Delta#varphi (rad.)", kTH2F, {{1000, 0.0f, 10.0f}, {100, 0.0f, 0.01f}}, false); fRegistry.add("Track/positive/hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); fRegistry.add("Track/positive/hDCAxyzSigma", "DCA xy vs. z;DCA_{xy} (#sigma);DCA_{z} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); fRegistry.add("Track/positive/hDCAxyRes_Pt", "DCA_{xy} resolution vs. pT;p_{T} (GeV/c);DCA_{xy} resolution (#mum)", kTH2F, {{200, 0, 10}, {200, 0., 400}}, false); diff --git a/PWGEM/Dilepton/Utils/EMTrackUtilities.h b/PWGEM/Dilepton/Utils/EMTrackUtilities.h index 1609f2f5823..01fdb8d95a0 100644 --- a/PWGEM/Dilepton/Utils/EMTrackUtilities.h +++ b/PWGEM/Dilepton/Utils/EMTrackUtilities.h @@ -77,8 +77,8 @@ float sigmaLambda(T const& track) template float sigmaP(T const& track) { - // p = pT x cosh(eta); - return std::sqrt(std::pow(std::cosh(track.eta()) * sigmaPt(track), 2) + std::pow(track.pt() * std::sinh(track.eta()) * std::cosh(track.eta()) * sigmaLambda(track), 2)); + // p = 1/1/pT x 1/cos(lambda); + return std::sqrt(std::pow(1.f / track.signed1Pt(), 4) * ((1.f + std::pow(track.tgl(), 2)) * track.c1Pt21Pt2() + 1.f / (1.f + std::pow(track.tgl(), 2)) * std::pow(track.signed1Pt() * track.tgl(), 2) * track.cTglTgl() - 2.f * track.signed1Pt() * track.tgl() * track.c1PtTgl())); } //_______________________________________________________________________ } // namespace o2::aod::pwgem::dilepton::utils::emtrackutil diff --git a/PWGEM/PhotonMeson/DataModel/gammaTables.h b/PWGEM/PhotonMeson/DataModel/gammaTables.h index fb6d9a7839b..0feef4df8f4 100644 --- a/PWGEM/PhotonMeson/DataModel/gammaTables.h +++ b/PWGEM/PhotonMeson/DataModel/gammaTables.h @@ -170,7 +170,7 @@ DECLARE_SOA_TABLE(V0Legs, "AOD", "V0LEG", //! o2::soa::Index<>, v0leg::CollisionId, v0leg::TrackId, v0leg::Sign, v0leg::Px, v0leg::Py, v0leg::Pz, track::DcaXY, track::DcaZ, - track::TPCNClsFindable, track::TPCNClsFindableMinusFound, track::TPCNClsFindableMinusCrossedRows, + track::TPCNClsFindable, track::TPCNClsFindableMinusFound, track::TPCNClsFindableMinusCrossedRows, track::TPCNClsShared, track::TPCChi2NCl, track::TPCInnerParam, track::TPCSignal, pidtpc::TPCNSigmaEl, pidtpc::TPCNSigmaPi, track::ITSClusterSizes, track::ITSChi2NCl, track::DetectorMap, @@ -185,6 +185,7 @@ DECLARE_SOA_TABLE(V0Legs, "AOD", "V0LEG", //! track::TPCNClsCrossedRows, track::TPCCrossedRowsOverFindableCls, track::TPCFoundOverFindableCls, + track::TPCFractionSharedCls, track::v001::ITSClusterMap, track::v001::ITSNCls, track::v001::ITSNClsInnerBarrel, track::HasITS, track::HasTPC, track::HasTRD, track::HasTOF, @@ -207,6 +208,7 @@ namespace v0photonkf { DECLARE_SOA_INDEX_COLUMN(EMEvent, emevent); //! DECLARE_SOA_COLUMN(CollisionId, collisionId, int); //! +DECLARE_SOA_COLUMN(V0Id, v0Id, int); //! DECLARE_SOA_INDEX_COLUMN_FULL(PosTrack, posTrack, int, V0Legs, "_Pos"); //! DECLARE_SOA_INDEX_COLUMN_FULL(NegTrack, negTrack, int, V0Legs, "_Neg"); //! DECLARE_SOA_COLUMN(Vx, vx, float); //! secondary vertex x @@ -223,6 +225,12 @@ DECLARE_SOA_COLUMN(PCA, pca, float); //! DECLARE_SOA_COLUMN(Alpha, alpha, float); //! DECLARE_SOA_COLUMN(QtArm, qtarm, float); //! DECLARE_SOA_COLUMN(ChiSquareNDF, chiSquareNDF, float); // Chi2 / NDF of the reconstructed V0 +DECLARE_SOA_COLUMN(SigmaPx2, sigmaPx2, float); //! error^2 of px in covariant matrix +DECLARE_SOA_COLUMN(SigmaPy2, sigmaPy2, float); //! error^2 of py in covariant matrix +DECLARE_SOA_COLUMN(SigmaPz2, sigmaPz2, float); //! error^2 of pz in covariant matrix +DECLARE_SOA_COLUMN(SigmaPxPy, sigmaPxPy, float); //! error of px x py in covariant matrix +DECLARE_SOA_COLUMN(SigmaPyPz, sigmaPyPz, float); //! error of py x pz in covariant matrix +DECLARE_SOA_COLUMN(SigmaPzPx, sigmaPzPx, float); //! error of pz x px in covariant matrix DECLARE_SOA_DYNAMIC_COLUMN(E, e, [](float px, float py, float pz, float m = 0) -> float { return RecoDecay::sqrtSumOfSquares(px, py, pz, m); }); //! energy of v0 photn, mass to be given as argument when getter is called! DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float px, float py) -> float { return RecoDecay::sqrtSumOfSquares(px, py); }); @@ -232,7 +240,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(P, p, [](float px, float py, float pz) -> float { ret DECLARE_SOA_DYNAMIC_COLUMN(V0Radius, v0radius, [](float vx, float vy) -> float { return RecoDecay::sqrtSumOfSquares(vx, vy); }); } // namespace v0photonkf DECLARE_SOA_TABLE(V0PhotonsKF, "AOD", "V0PHOTONKF", //! - o2::soa::Index<>, v0photonkf::CollisionId, v0photonkf::PosTrackId, v0photonkf::NegTrackId, + o2::soa::Index<>, v0photonkf::CollisionId, v0photonkf::V0Id, v0photonkf::PosTrackId, v0photonkf::NegTrackId, v0photonkf::Vx, v0photonkf::Vy, v0photonkf::Vz, v0photonkf::Px, v0photonkf::Py, v0photonkf::Pz, v0photonkf::MGamma, @@ -255,6 +263,11 @@ DECLARE_SOA_TABLE(V0KFEMEventIds, "AOD", "V0KFEMEVENTID", v0photonkf::EMEventId) // iterators using V0KFEMEventId = V0KFEMEventIds::iterator; +DECLARE_SOA_TABLE(V0PhotonsKFCov, "AOD", "V0PHOTONKFCOV", //! To be joined with V0PhotonsKF table at analysis level. + v0photonkf::SigmaPx2, v0photonkf::SigmaPy2, v0photonkf::SigmaPz2, v0photonkf::SigmaPxPy, v0photonkf::SigmaPyPz, v0photonkf::SigmaPzPx); +// iterators +using V0PhotonKFCov = V0PhotonsKFCov::iterator; + DECLARE_SOA_TABLE(EMPrimaryElectronsFromDalitz, "AOD", "EMPRIMARYELDA", //! o2::soa::Index<>, emprimaryelectron::CollisionId, emprimaryelectron::TrackId, emprimaryelectron::Sign, diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx index 226801fcc61..9786329b601 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx @@ -65,6 +65,7 @@ using MyTracksIUMC = soa::Join; struct PhotonConversionBuilder { Produces v0photonskf; + Produces v0photonskfcov; Produces v0legs; Produces events_ngpcm; @@ -370,7 +371,7 @@ struct PhotonConversionBuilder { { v0legs(track.collisionId(), track.globalIndex(), track.sign(), kfp.GetPx(), kfp.GetPy(), kfp.GetPz(), dcaXY, dcaZ, - track.tpcNClsFindable(), track.tpcNClsFindableMinusFound(), track.tpcNClsFindableMinusCrossedRows(), + track.tpcNClsFindable(), track.tpcNClsFindableMinusFound(), track.tpcNClsFindableMinusCrossedRows(), track.tpcNClsShared(), track.tpcChi2NCl(), track.tpcInnerParam(), track.tpcSignal(), track.tpcNSigmaEl(), track.tpcNSigmaPi(), track.itsClusterSizes(), track.itsChi2NCl(), track.detectorMap(), @@ -660,12 +661,14 @@ struct PhotonConversionBuilder { ROOT::Math::PxPyPzMVector v0_sv = vpos_sv + vele_sv; registry.fill(HIST("V0/hMeeSV_Rxy"), rxy, v0_sv.M()); - v0photonskf(collision.globalIndex(), v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, + v0photonskf(collision.globalIndex(), v0.globalIndex(), v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, gammaKF_DecayVtx.GetX(), gammaKF_DecayVtx.GetY(), gammaKF_DecayVtx.GetZ(), gammaKF_PV.GetPx(), gammaKF_PV.GetPy(), gammaKF_PV.GetPz(), v0_sv.M(), dca_xy_v0_to_pv, dca_z_v0_to_pv, cospa_kf, pca_kf, alpha, qt, chi2kf); + v0photonskfcov(gammaKF_PV.GetCovariance(9), gammaKF_PV.GetCovariance(14), gammaKF_PV.GetCovariance(20), gammaKF_PV.GetCovariance(13), gammaKF_PV.GetCovariance(19), gammaKF_PV.GetCovariance(18)); + fillTrackTable(pos, pTrack, kfp_pos_DecayVtx, posdcaXY, posdcaZ); // positive leg first fillTrackTable(ele, nTrack, kfp_ele_DecayVtx, eledcaXY, eledcaZ); // negative leg second } // end of fill table diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversion.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversion.cxx index 40334a036c0..bed52425fb5 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversion.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversion.cxx @@ -181,7 +181,7 @@ struct skimmerGammaConversion { v0legs(theTrack.collisionId(), theTrack.globalIndex(), theTrack.sign(), kfp.GetPx(), kfp.GetPy(), kfp.GetPz(), theTrack.dcaXY(), theTrack.dcaZ(), - theTrack.tpcNClsFindable(), theTrack.tpcNClsFindableMinusFound(), theTrack.tpcNClsFindableMinusCrossedRows(), + theTrack.tpcNClsFindable(), theTrack.tpcNClsFindableMinusFound(), theTrack.tpcNClsFindableMinusCrossedRows(), theTrack.tpcNClsShared(), theTrack.tpcChi2NCl(), theTrack.tpcInnerParam(), theTrack.tpcSignal(), theTrack.tpcNSigmaEl(), theTrack.tpcNSigmaPi(), theTrack.itsClusterSizes(), theTrack.itsChi2NCl(), theTrack.detectorMap(), @@ -304,7 +304,7 @@ struct skimmerGammaConversion { float sign_tmp = dca_y_v0_to_pv > 0 ? +1 : -1; float dca_xy_v0_to_pv = RecoDecay::sqrtSumOfSquares(dca_x_v0_to_pv, dca_y_v0_to_pv) * sign_tmp; - v0photonskf(collision.globalIndex(), v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, + v0photonskf(collision.globalIndex(), v0.globalIndex(), v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, gammaKF_DecayVtx.GetX(), gammaKF_DecayVtx.GetY(), gammaKF_DecayVtx.GetZ(), gammaKF_DecayVtx.GetPx(), gammaKF_DecayVtx.GetPy(), gammaKF_DecayVtx.GetPz(), v0_sv.M(), dca_xy_v0_to_pv, dca_z_v0_to_pv, @@ -342,7 +342,7 @@ struct skimmerGammaConversion { fillV0KF(collision, v0); } // end of v0 loop - } // end of collision loop + } // end of collision loop } PROCESS_SWITCH(skimmerGammaConversion, processRec, "process reconstructed info only", true); diff --git a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx index 49f734c1a9d..bc08e6a76af 100644 --- a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx @@ -34,7 +34,7 @@ using namespace o2::aod::pwgem::photon; using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; -using MyV0Photons = soa::Join; +using MyV0Photons = soa::Join; using MyV0Photon = MyV0Photons::iterator; struct PCMQC { @@ -66,7 +66,8 @@ struct PCMQC { Configurable cfg_require_v0_with_tpconly{"cfg_require_v0_with_tpconly", false, "flag to select V0s with TPConly tracks"}; Configurable cfg_require_v0_on_wwire_ib{"cfg_require_v0_on_wwire_ib", false, "flag to select V0s on W wires ITSib"}; Configurable cfg_min_pt_v0{"cfg_min_pt_v0", 0.1, "min pT for v0 photons at PV"}; - Configurable cfg_max_eta_v0{"cfg_max_eta_v0", 0.8, "max eta for v0 photons at PV"}; + Configurable cfg_min_eta_v0{"cfg_min_eta_v0", -0.8, "min eta for v0 photons at PV"}; + Configurable cfg_max_eta_v0{"cfg_max_eta_v0", +0.8, "max eta for v0 photons at PV"}; Configurable cfg_min_v0radius{"cfg_min_v0radius", 4.0, "min v0 radius"}; Configurable cfg_max_v0radius{"cfg_max_v0radius", 90.0, "max v0 radius"}; Configurable cfg_max_alpha_ap{"cfg_max_alpha_ap", 0.95, "max alpha for AP cut"}; @@ -138,6 +139,11 @@ struct PCMQC { fRegistry.add("V0/hKFChi2vsX", "KF chi2 vs. conversion point in X;X (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); fRegistry.add("V0/hKFChi2vsY", "KF chi2 vs. conversion point in Y;Y (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); fRegistry.add("V0/hKFChi2vsZ", "KF chi2 vs. conversion point in Z;Z (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); + fRegistry.add("V0/hPResolution", "p resolution;p_{#gamma} (GeV/c);#Deltap/p", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.1}}, false); + fRegistry.add("V0/hPtResolution", "p_{T} resolution;p_{#gamma} (GeV/c);#Deltap_{T}/p_{T}", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.1}}, false); + fRegistry.add("V0/hEtaResolution", "#eta resolution;p_{#gamma} (GeV/c);#Delta#eta", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.01}}, false); + fRegistry.add("V0/hThetaResolution", "#theta resolution;p_{#gamma} (GeV/c);#Delta#theta (rad.)", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.01}}, false); + fRegistry.add("V0/hPhiResolution", "#varphi resolution;p_{#gamma} (GeV/c);#Delta#varphi (rad.)", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.01}}, false); fRegistry.add("V0/hNgamma", "Number of #gamma candidates per collision", kTH1F, {{101, -0.5f, 100.5f}}); // v0leg info @@ -183,7 +189,7 @@ struct PCMQC { // for v0 fV0PhotonCut.SetV0PtRange(pcmcuts.cfg_min_pt_v0, 1e10f); - fV0PhotonCut.SetV0EtaRange(-pcmcuts.cfg_max_eta_v0, +pcmcuts.cfg_max_eta_v0); + fV0PhotonCut.SetV0EtaRange(pcmcuts.cfg_min_eta_v0, pcmcuts.cfg_max_eta_v0); fV0PhotonCut.SetMinCosPA(pcmcuts.cfg_min_cospa); fV0PhotonCut.SetMaxPCA(pcmcuts.cfg_max_pca); fV0PhotonCut.SetRxyRange(pcmcuts.cfg_min_v0radius, pcmcuts.cfg_max_v0radius); @@ -192,7 +198,7 @@ struct PCMQC { // for track fV0PhotonCut.SetTrackPtRange(pcmcuts.cfg_min_pt_v0 * 0.4, 1e+10f); - fV0PhotonCut.SetTrackEtaRange(-pcmcuts.cfg_max_eta_v0, +pcmcuts.cfg_max_eta_v0); + fV0PhotonCut.SetTrackEtaRange(pcmcuts.cfg_min_eta_v0, pcmcuts.cfg_max_eta_v0); fV0PhotonCut.SetMinNClustersTPC(pcmcuts.cfg_min_ncluster_tpc); fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfg_min_ncrossedrows); fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); @@ -292,6 +298,11 @@ struct PCMQC { fRegistry.fill(HIST("V0/hKFChi2vsX"), v0.vx(), v0.chiSquareNDF()); fRegistry.fill(HIST("V0/hKFChi2vsY"), v0.vy(), v0.chiSquareNDF()); fRegistry.fill(HIST("V0/hKFChi2vsZ"), v0.vz(), v0.chiSquareNDF()); + fRegistry.fill(HIST("V0/hPResolution"), v0.p(), getPResolution(v0) / v0.p()); + fRegistry.fill(HIST("V0/hPtResolution"), v0.p(), getPtResolution(v0) / v0.pt()); + fRegistry.fill(HIST("V0/hEtaResolution"), v0.p(), getEtaResolution(v0)); + fRegistry.fill(HIST("V0/hThetaResolution"), v0.p(), getThetaResolution(v0)); + fRegistry.fill(HIST("V0/hPhiResolution"), v0.p(), getPhiResolution(v0)); } template diff --git a/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx b/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx index 279c0a57c80..1fc481378fe 100644 --- a/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx @@ -41,7 +41,7 @@ using MyCollision = MyCollisions::iterator; using MyMCCollisions = soa::Join; using MyMCCollision = MyMCCollisions::iterator; -using MyV0Photons = soa::Join; +using MyV0Photons = soa::Join; using MyV0Photon = MyV0Photons::iterator; using MyMCV0Legs = soa::Join; @@ -80,7 +80,8 @@ struct PCMQCMC { Configurable cfg_require_v0_with_tpconly{"cfg_require_v0_with_tpconly", false, "flag to select V0s with TPConly tracks"}; Configurable cfg_require_v0_on_wwire_ib{"cfg_require_v0_on_wwire_ib", false, "flag to select V0s on W wires ITSib"}; Configurable cfg_min_pt_v0{"cfg_min_pt_v0", 0.1, "min pT for v0 photons at PV"}; - Configurable cfg_max_eta_v0{"cfg_max_eta_v0", 0.8, "max eta for v0 photons at PV"}; + Configurable cfg_min_eta_v0{"cfg_min_eta_v0", -0.8, "min eta for v0 photons at PV"}; + Configurable cfg_max_eta_v0{"cfg_max_eta_v0", +0.8, "max eta for v0 photons at PV"}; Configurable cfg_min_v0radius{"cfg_min_v0radius", 4.0, "min v0 radius"}; Configurable cfg_max_v0radius{"cfg_max_v0radius", 90.0, "max v0 radius"}; Configurable cfg_max_alpha_ap{"cfg_max_alpha_ap", 0.95, "max alpha for AP cut"}; @@ -180,6 +181,11 @@ struct PCMQCMC { fRegistry.add("V0/primary/hKFChi2vsX", "KF chi2 vs. conversion point in X;X (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); fRegistry.add("V0/primary/hKFChi2vsY", "KF chi2 vs. conversion point in Y;Y (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); fRegistry.add("V0/primary/hKFChi2vsZ", "KF chi2 vs. conversion point in Z;Z (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); + fRegistry.add("V0/primary/hPResolution", "p resolution;p_{#gamma} (GeV/c);#Deltap/p", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.1}}, false); + fRegistry.add("V0/primary/hPtResolution", "p_{T} resolution;p_{#gamma} (GeV/c);#Deltap_{T}/p_{T}", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.1}}, false); + fRegistry.add("V0/primary/hEtaResolution", "#eta resolution;p_{#gamma} (GeV/c);#Delta#eta", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.01}}, false); + fRegistry.add("V0/primary/hThetaResolution", "#theta resolution;p_{#gamma} (GeV/c);#Delta#theta (rad.)", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.01}}, false); + fRegistry.add("V0/primary/hPhiResolution", "#varphi resolution;p_{#gamma} (GeV/c);#Delta#varphi (rad.)", kTH2F, {{1000, 0.0f, 10}, {100, 0, 0.01}}, false); fRegistry.add("V0/primary/hNgamma", "Number of true #gamma per collision", kTH1F, {{101, -0.5f, 100.5f}}); fRegistry.add("V0/primary/hConvPoint_diffX", "conversion point diff X MC;X_{MC} (cm);X_{rec} - X_{MC} (cm)", kTH2F, {{200, -100, +100}, {100, -50.0f, 50.0f}}, true); fRegistry.add("V0/primary/hConvPoint_diffY", "conversion point diff Y MC;Y_{MC} (cm);Y_{rec} - Y_{MC} (cm)", kTH2F, {{200, -100, +100}, {100, -50.0f, 50.0f}}, true); @@ -240,7 +246,7 @@ struct PCMQCMC { // for v0 fV0PhotonCut.SetV0PtRange(pcmcuts.cfg_min_pt_v0, 1e10f); - fV0PhotonCut.SetV0EtaRange(-pcmcuts.cfg_max_eta_v0, +pcmcuts.cfg_max_eta_v0); + fV0PhotonCut.SetV0EtaRange(pcmcuts.cfg_min_eta_v0, pcmcuts.cfg_max_eta_v0); fV0PhotonCut.SetMinCosPA(pcmcuts.cfg_min_cospa); fV0PhotonCut.SetMaxPCA(pcmcuts.cfg_max_pca); fV0PhotonCut.SetRxyRange(pcmcuts.cfg_min_v0radius, pcmcuts.cfg_max_v0radius); @@ -249,7 +255,7 @@ struct PCMQCMC { // for track fV0PhotonCut.SetTrackPtRange(pcmcuts.cfg_min_pt_v0 * 0.4, 1e+10f); - fV0PhotonCut.SetTrackEtaRange(-pcmcuts.cfg_max_eta_v0, +pcmcuts.cfg_max_eta_v0); + fV0PhotonCut.SetTrackEtaRange(pcmcuts.cfg_min_eta_v0, pcmcuts.cfg_max_eta_v0); fV0PhotonCut.SetMinNClustersTPC(pcmcuts.cfg_min_ncluster_tpc); fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfg_min_ncrossedrows); fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); @@ -349,6 +355,11 @@ struct PCMQCMC { fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hKFChi2vsX"), v0.vx(), v0.chiSquareNDF()); fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hKFChi2vsY"), v0.vy(), v0.chiSquareNDF()); fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hKFChi2vsZ"), v0.vz(), v0.chiSquareNDF()); + fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hPResolution"), v0.p(), getPResolution(v0) / v0.p()); + fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hPtResolution"), v0.p(), getPtResolution(v0) / v0.pt()); + fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hEtaResolution"), v0.p(), getEtaResolution(v0)); + fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hThetaResolution"), v0.p(), getThetaResolution(v0)); + fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hPhiResolution"), v0.p(), getPhiResolution(v0)); fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hPtGen_DeltaPtOverPtGen"), mcphoton.pt(), (v0.pt() - mcphoton.pt()) / mcphoton.pt()); fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hPtGen_DeltaEta"), mcphoton.pt(), v0.eta() - mcphoton.eta()); fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hPtGen_DeltaPhi"), mcphoton.pt(), v0.phi() - mcphoton.phi()); diff --git a/PWGEM/PhotonMeson/Utils/PCMUtilities.h b/PWGEM/PhotonMeson/Utils/PCMUtilities.h index 9530fa03c4f..63101c41418 100644 --- a/PWGEM/PhotonMeson/Utils/PCMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PCMUtilities.h @@ -104,4 +104,80 @@ inline void Vtx_recalculation(o2::base::Propagator* prop, T1 lTrackPos, T2 lTrac Vtx_recalculationParCov(prop, trackPosInformation, trackNegInformation, xyz, matCorr); } +//_______________________________________________________________________ +template +float getPtResolution(TV0 const& v0) +{ + float px = v0.px(); + float py = v0.py(); + float pt = v0.pt(); + float px_err = std::sqrt(fabs(v0.sigmaPx2())); + float py_err = std::sqrt(fabs(v0.sigmaPy2())); + float pxy_err = v0.sigmaPxPy(); + return std::sqrt(std::pow(px / pt * px_err, 2) + std::pow(py / pt * py_err, 2) + 2.f * px / pt * py / pt * pxy_err); +} +//_______________________________________________________________________ +template +float getPhiResolution(TV0 const& v0) +{ + float px = v0.px(); + float py = v0.py(); + float pt = v0.pt(); + float px_err = std::sqrt(fabs(v0.sigmaPx2())); + float py_err = std::sqrt(fabs(v0.sigmaPy2())); + float pxy_err = v0.sigmaPxPy(); + return std::sqrt(std::pow(px / pt / pt * py_err, 2) + std::pow(py / pt / pt * px_err, 2) - 2.f * px / pt / pt * py / pt / pt * pxy_err); +} +//_______________________________________________________________________ +template +float getThetaResolution(TV0 const& v0) +{ + float px = v0.px(); + float py = v0.py(); + float pz = v0.pz(); + float pt = v0.pt(); + float p = v0.p(); + float px_err = std::sqrt(fabs(v0.sigmaPx2())); + float py_err = std::sqrt(fabs(v0.sigmaPy2())); + float pz_err = std::sqrt(fabs(v0.sigmaPz2())); + float pxy_err = v0.sigmaPxPy(); + float pyz_err = v0.sigmaPyPz(); + float pzx_err = v0.sigmaPzPx(); + return std::sqrt(std::pow(pz * pz / p / p, 2) * (std::pow(px / pz / pt * px_err, 2) + std::pow(py / pz / pt * py_err, 2) + std::pow(pt / pz / pz * pz_err, 2) + 2.f * (px * py / pz / pz / pt / pt * pxy_err - py / pz / pz / pz * pyz_err - px / pz / pz / pz * pzx_err))); +} +//_______________________________________________________________________ +template +float getEtaResolution(TV0 const& v0) +{ + float px = v0.px(); + float py = v0.py(); + float pz = v0.pz(); + float pt = v0.pt(); + float p = v0.p(); + float px_err = std::sqrt(fabs(v0.sigmaPx2())); + float py_err = std::sqrt(fabs(v0.sigmaPy2())); + float pz_err = std::sqrt(fabs(v0.sigmaPz2())); + float pxy_err = v0.sigmaPxPy(); + float pyz_err = v0.sigmaPyPz(); + float pzx_err = v0.sigmaPzPx(); + return std::sqrt(std::pow(1.f / p / pt / pt, 2) * (std::pow(pz * px * px_err, 2) + std::pow(pz * py * py_err, 2) + std::pow(pt * pt * pz_err, 2) + 2.f * (pz * pz * px * py * pxy_err - pt * pt * py * pz * pyz_err - pt * pt * pz * px * pzx_err))); +} +//_______________________________________________________________________ +template +float getPResolution(TV0 const& v0) +{ + float px = v0.px(); + float py = v0.py(); + float pz = v0.pz(); + float p = v0.p(); + float px_err = std::sqrt(fabs(v0.sigmaPx2())); + float py_err = std::sqrt(fabs(v0.sigmaPy2())); + float pz_err = std::sqrt(fabs(v0.sigmaPz2())); + float pxy_err = v0.sigmaPxPy(); + float pyz_err = v0.sigmaPyPz(); + float pzx_err = v0.sigmaPzPx(); + return std::sqrt(std::pow(1.f / p, 2) * (std::pow(px * px_err, 2) + std::pow(py * py_err, 2) + std::pow(pz * pz_err, 2) + 2.f * (px * py * pxy_err + py * pz * pyz_err + pz * px * pzx_err))); +} +//_______________________________________________________________________ +//_______________________________________________________________________ #endif // PWGEM_PHOTONMESON_UTILS_PCMUTILITIES_H_ From b256183d4756bc9493babaf4eca1f7e36486435d Mon Sep 17 00:00:00 2001 From: feisenhu <53603353+feisenhu@users.noreply.github.com> Date: Sat, 14 Sep 2024 12:08:39 +0200 Subject: [PATCH 14/14] [PWGDQ] add LS histograms to output of MC selected signals (#7676) * [PWGDQ] add LS histograms to output of MC selected signals * [PWGDQ] fix filling of configQA histos --- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 62 +++++++++++++++----------- 1 file changed, 35 insertions(+), 27 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 0272bfb1e0e..9624b322309 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -1205,22 +1205,22 @@ struct AnalysisSameEventPairing { if (!sigNamesStr.IsNull()) { for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { auto sig = fRecMCSignals.at(isig); + names = { + Form("PairsBarrelSEPM_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), + Form("PairsBarrelSEPP_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), + Form("PairsBarrelSEMM_%s_%s", objArray->At(icut)->GetName(), sig.GetName())}; if (fConfigQA) { - names = { - Form("PairsBarrelSEPM_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPMCorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPMIncorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPM_ambiguousInBunch_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPM_ambiguousInBunchCorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPM_ambiguousInBunchIncorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPM_ambiguousOutOfBunch_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPM_ambiguousOutOfBunchCorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName()), - Form("PairsBarrelSEPM_ambiguousOutOfBunchIncorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName())}; - histNames += Form("%s;%s;%s;%s;%s;%s;%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data(), names[3].Data(), names[4].Data(), names[5].Data(), names[6].Data(), names[7].Data(), names[8].Data()); - } else { - names = { - Form("PairsBarrelSEPM_%s_%s", objArray->At(icut)->GetName(), sig.GetName())}; - histNames += Form("%s;", names[0].Data()); + names.push_back(Form("PairsBarrelSEPMCorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + names.push_back(Form("PairsBarrelSEPMIncorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + names.push_back(Form("PairsBarrelSEPM_ambiguousInBunch_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + names.push_back(Form("PairsBarrelSEPM_ambiguousInBunchCorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + names.push_back(Form("PairsBarrelSEPM_ambiguousInBunchIncorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + names.push_back(Form("PairsBarrelSEPM_ambiguousOutOfBunch_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + names.push_back(Form("PairsBarrelSEPM_ambiguousOutOfBunchCorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + names.push_back(Form("PairsBarrelSEPM_ambiguousOutOfBunchIncorrectAssoc_%s_%s", objArray->At(icut)->GetName(), sig.GetName())); + } + for (auto& n : names) { + histNames += Form("%s;", n.Data()); } fBarrelHistNamesMCmatched.try_emplace(icut * fRecMCSignals.size() + isig, names); } // end loop over MC signals @@ -1579,37 +1579,40 @@ struct AnalysisSameEventPairing { fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues); // reconstructed, unmatched for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals if (mcDecision & (uint32_t(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][0].Data(), VarManager::fgValues); // matched signal if (fConfigQA) { - fHistMan->FillHistClass(histNamesMC[icut][isig * fRecMCSignals.size()].Data(), VarManager::fgValues); // matched signal - if (isCorrectAssoc_leg1 && isCorrectAssoc_leg2) { // correct track-collision association - fHistMan->FillHistClass(histNamesMC[icut][isig * fRecMCSignals.size() + 1].Data(), VarManager::fgValues); + if (isCorrectAssoc_leg1 && isCorrectAssoc_leg2) { // correct track-collision association + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][3].Data(), VarManager::fgValues); } else { // incorrect track-collision association - fHistMan->FillHistClass(histNamesMC[icut][isig * fRecMCSignals.size() + 2].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][4].Data(), VarManager::fgValues); } if (isAmbiInBunch) { // ambiguous in bunch - fHistMan->FillHistClass(histNames[icut][isig * fRecMCSignals.size() + 3].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut * fRecMCSignals.size() + isig][5].Data(), VarManager::fgValues); if (isCorrectAssoc_leg1 && isCorrectAssoc_leg2) { - fHistMan->FillHistClass(histNamesMC[icut][isig * fRecMCSignals.size() + 4].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][6].Data(), VarManager::fgValues); } else { - fHistMan->FillHistClass(histNamesMC[icut][isig * fRecMCSignals.size() + 5].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][7].Data(), VarManager::fgValues); } } if (isAmbiOutOfBunch) { // ambiguous out of bunch - fHistMan->FillHistClass(histNames[icut][isig * fRecMCSignals.size() + 6].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut * fRecMCSignals.size() + isig][8].Data(), VarManager::fgValues); if (isCorrectAssoc_leg1 && isCorrectAssoc_leg2) { - fHistMan->FillHistClass(histNamesMC[icut][isig * fRecMCSignals.size() + 7].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][9].Data(), VarManager::fgValues); } else { - fHistMan->FillHistClass(histNamesMC[icut][isig * fRecMCSignals.size() + 8].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][10].Data(), VarManager::fgValues); } } - } else { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][0].Data(), VarManager::fgValues); // matched signal } } } } else { if (sign1 > 0) { // ++ pairs fHistMan->FillHistClass(histNames[icut][1].Data(), VarManager::fgValues); + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + if (mcDecision & (uint32_t(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][1].Data(), VarManager::fgValues); + } + } if (fConfigQA) { if (isAmbiInBunch) { fHistMan->FillHistClass(histNames[icut][4].Data(), VarManager::fgValues); @@ -1620,6 +1623,11 @@ struct AnalysisSameEventPairing { } } else { // -- pairs fHistMan->FillHistClass(histNames[icut][2].Data(), VarManager::fgValues); + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + if (mcDecision & (uint32_t(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][2].Data(), VarManager::fgValues); + } + } if (fConfigQA) { if (isAmbiInBunch) { fHistMan->FillHistClass(histNames[icut][5].Data(), VarManager::fgValues);