From 0bf157045d95a81813b3307391fb45b11a00fef0 Mon Sep 17 00:00:00 2001 From: Veronika Barbasova Date: Thu, 8 Jan 2026 08:59:16 +0100 Subject: [PATCH] [PWGLF] Rotational bg improvement and resolution histograms --- .../Tasks/Resonances/phianalysisTHnSparse.cxx | 327 ++++++++++++------ 1 file changed, 212 insertions(+), 115 deletions(-) diff --git a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx index 49a754b6a4b..27d4ed55de3 100644 --- a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx +++ b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx @@ -71,6 +71,7 @@ struct PhianalysisTHnSparse { Configurable globalTrack{"globalTrack", false, "Use global track selection."}; Configurable inelGrater0{"inelGrater0", true, "Select events with INEL>0."}; Configurable tpcNClsFound{"tpcNClsFound", 70, "Cut: Minimal value of found TPC clasters"}; + Configurable vzCut{"vzCut", 10.0f, "Cut: Maximal value of Z vertex position."}; } cut; struct : ConfigurableGroup { @@ -102,14 +103,19 @@ struct PhianalysisTHnSparse { // rotational Configurable numberofRotations{"numberofRotations", 1, "Number of rotations for rotational background estimation."}; + Configurable startingAngle{"startingAngle", 0, "Starting angle for rotational background estimation."}; - // Normalization factors - ConfigurableAxis axisNch{"axisNch", {100, 0.0f, +100.0f}, "Number of charged particles in |y| < 0.5"}; + // other axes + ConfigurableAxis axisNch{"axisNch", {1000, 0.0f, +1000.0f}, "Number of charged particles."}; + ConfigurableAxis axisResolutionPt{"axisResolutionPt", {1001, -1.0f, +1.0f}, "Resolution of Pt."}; + ConfigurableAxis axisResolutionPtPhi{"axisResolutionPtPhi", {1001, -0.0001f, +0.0001f}, "Resolution of Pt and Phi."}; + ConfigurableAxis axisResolutionMass{"axisResolutionMass", {1001, -0.05f, +0.05f}, "Resolution of Mass."}; + ConfigurableAxis axisResolutionVz{"axisResolutionVz", {1001, -3.0f, +3.0f}, "Resolution of Vz."}; // Axes specifications AxisSpec posZaxis = {400, -20., 20., "V_{z} (cm)"}; - AxisSpec dcaXYaxis = {800, -2.0, 2.0, "DCA_{xy} (cm)"}; - AxisSpec dcaZaxis = {800, -2.0, 2.0, "DCA_{z} (cm)"}; + AxisSpec dcaXYaxis = {1000, -1.0, 1.0, "DCA_{xy} (cm)"}; + AxisSpec dcaZaxis = {1000, -1.0, 1.0, "DCA_{z} (cm)"}; AxisSpec etaQAaxis = {1000, -1.0, 1.0, "#eta"}; AxisSpec tpcNClsFoundQAaxis = {110, 50., 160., "tpcNClsFound"}; @@ -131,18 +137,17 @@ struct PhianalysisTHnSparse { float ptTOFThreshold = 0.5f; int tpcNClsFound = 70; int dauSize = 2; - + float vzCut = 10.0f; rsn::MixingType mixingType = rsn::MixingType::none; - float vzCut = 10.0f; Filter triggerFilter = (o2::aod::evsel::sel8 == true); Filter vtxFilter = (nabs(o2::aod::collision::posZ) < vzCut); - using EventCandidates = soa::Filtered>; + using EventCandidates = soa::Join; using EventCandidate = EventCandidates::iterator; - using TrackCandidates = soa::Join; + using TrackCandidates = soa::Join; - using EventCandidatesMC = soa::Filtered>; + using EventCandidatesMC = soa::Join; using TrackCandidatesMC = soa::Join; using EventCandidatesMCGen = soa::Join; @@ -201,6 +206,7 @@ struct PhianalysisTHnSparse { combinedNSigma = static_cast(cut.combinedNSigma); tpcPidOnly = static_cast(cut.tpcPidOnly); ptTOFThreshold = static_cast(cut.ptTOFThreshold); + vzCut = static_cast(cut.vzCut); pointPair = new double[static_cast(o2::analysis::rsn::PairAxisType::unknown)]; pointSys = new double[static_cast(o2::analysis::rsn::SystematicsAxisType::unknown)]; @@ -213,6 +219,7 @@ struct PhianalysisTHnSparse { LOGF(info, "produceStats: %s", produceStats ? "true" : "false"); LOGF(info, "produceTrue: %s", static_cast(produce.produceTrue) ? "true" : "false"); LOGF(info, "produceLikesign: %s", static_cast(produce.produceLikesign) ? "true" : "false"); + LOGF(info, "produceRotational: %s", static_cast(produce.produceRotational) ? "true" : "false"); LOGF(info, "eventMixing: %s", static_cast(produce.eventMixing).c_str()); LOGF(info, "daughterPos: %d (PDG: %d)", static_cast(daughterPos), static_cast(daughterPosPDG)); LOGF(info, "daughterNeg: %d (PDG: %d)", static_cast(daughterNeg), static_cast(daughterNegPDG)); @@ -230,9 +237,11 @@ struct PhianalysisTHnSparse { LOGF(info, "globalTrack: %s", globalTrack ? "true" : "false"); LOGF(info, "inelGrater0: %s", inelGrater0 ? "true" : "false"); LOGF(info, "tpcNClsFound: %d", tpcNClsFound); + LOGF(info, "vzCut: %.2f", vzCut); LOGF(info, "mixingType: %d", static_cast(mixingType)); LOGF(info, "numberofMixedEvents: %d", static_cast(numberofMixedEvents)); LOGF(info, "numberofRotations: %d", static_cast(numberofRotations)); + LOGF(info, "startingAngle: %d", static_cast(startingAngle)); LOGF(info, "sparseAxes: "); for (const auto& axis : static_cast>(sparseAxes)) { LOGF(info, " %s", axis.c_str()); @@ -247,15 +256,16 @@ struct PhianalysisTHnSparse { // Event QA registry.add("QAEvent/hSelection", "Event selection statistics", kTH1D, {{4, 0.0f, 4.0f}}); auto hEvent = registry.get(HIST("QAEvent/hSelection")); - hEvent->GetXaxis()->SetBinLabel(1, "Full event statistics"); - hEvent->GetXaxis()->SetBinLabel(2, "Events with at least one primary vertex"); - hEvent->GetXaxis()->SetBinLabel(3, "Events with at least one #phi candidate"); - hEvent->GetXaxis()->SetBinLabel(4, "Number of #phi candidates"); + + hEvent->GetXaxis()->SetBinLabel(1, "all events"); + hEvent->GetXaxis()->SetBinLabel(2, "Events passing trigger sel8"); + hEvent->GetXaxis()->SetBinLabel(3, "Events passing |V_{z}| cut"); + hEvent->GetXaxis()->SetBinLabel(4, "Events passing INEL>0 cut"); hEvent->SetMinimum(0.1); registry.add("QAEvent/hVtxZ", "Vertex position along the z-axis", kTH1F, {posZaxis}); registry.add("QAEvent/hCent", "Distribution of multiplicity percentile", kTH1F, {{101, 0., 101.}}); - registry.add("QAEvent/hMult", "Multiplicity (amplitude of non-zero channels in the FV0A + FV0C) ", kTH1F, {{300, 0., 30000.}}); + registry.add("QAEvent/hMult", "Multiplicity (amplitude of non-zero channels in the FT0A + FT0C) ", kTH1F, {{300, 0., 30000.}}); // Track QA registry.add("QATrack/hSelection", "Tracks statistics", kTH1D, {{9, 0.0f, 9.0f}}); @@ -278,63 +288,51 @@ struct PhianalysisTHnSparse { registry.add("QATrack/hDCAz", "Distribution of DCA_{z} of K^{+} and K^{-}", kTH1F, {dcaZaxis}); registry.add("QATrack/hPt", "Distribution of p_{T} of K^{+} and K^{-}", kTH1F, {ptaxis}); - // General Phi distributions - registry.add("QAPhi/hRapidity", "Rapidity distribution of #phi candidates", kTH1F, {rapidityaxis}); - registry.add("QAPhi/hEta", "Pseudorapidity distribution of #phi candidates", kTH1F, {etaaxis}); - - registry.add("QAPhi/hdPhi", "Azimuthal distribution of #phi candidates", kTH1F, {{100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); + // Phi candidate QA + registry.add("QAPhi/hRapidity", "Rapidity distribution of #Phi candidates", kTH1F, {{200, -1, 1}}); + registry.add("QAPhi/hEta", "Pseudorapidity distribution of #Phi candidates", kTH1F, {{200, -1, 1}}); + registry.add("QAPhi/hdPhi", "Azimuthal distribution (#Delta#phi) of #Phi candidates", kTH1F, {{100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); auto hdPhi = registry.get(HIST("QAPhi/hdPhi")); - hdPhi->GetXaxis()->SetTitle("#phi (rad)"); + hdPhi->GetXaxis()->SetTitle("#Delta#phi (rad)"); - registry.add("QAPhi/h2dPhiPt", "Azimuthal distribution of #Delta#phi candidates vs p_{T}", kTH2F, {ptaxis, {100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); + registry.add("QAPhi/h2dPhiPt", "Azimuthal distribution (#Delta#phi) of #Phi candidates vs p_{T}", kTH2F, {ptaxis, {100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); auto h2dPhiPt = registry.get(HIST("QAPhi/h2dPhiPt")); h2dPhiPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); h2dPhiPt->GetYaxis()->SetTitle("#Delta#phi (rad)"); - registry.add("QAPhi/hTheta", "Polar distribution of #phi candidates", kTH1F, {{100, 0.0f, o2::constants::math::PI}}); + registry.add("QAPhi/hTheta", "Polar distribution of #Phi candidates", kTH1F, {{100, 0.0f, o2::constants::math::PI}}); auto hTheta = registry.get(HIST("QAPhi/hTheta")); hTheta->GetXaxis()->SetTitle("#theta (rad)"); - registry.add("QAPhi/h2dThetaPt", "Polar distribution of #phi candidates vs p_{T}", kTH2F, {{12, 0, 12}, {100, -o2::constants::math::PI, o2::constants::math::PI}}); + registry.add("QAPhi/h2dThetaPt", "Polar distribution (#Delta#theta) of #Phi candidates vs p_{T}", kTH2F, {ptaxis, {100, -o2::constants::math::PI, o2::constants::math::PI}}); + auto h2dThetaPt = registry.get(HIST("QAPhi/h2dThetaPt")); h2dThetaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - h2dThetaPt->GetYaxis()->SetTitle("#theta (rad)"); - - // TPC PID - registry.add("QATrack/hTPCnSigma", "Distribution of TPC nSigma of K^{+} and K^{-}", kTH1F, {{200, -10, 10}}); - - registry.add("QATrack/h2TPCnSigma", "", kTH2F, {{200, -10, 10}, {200, -10, 10}}); - auto h2TPCnSigma = registry.get(HIST("QATrack/h2TPCnSigma")); - h2TPCnSigma->GetXaxis()->SetTitle("n#sigma_{TPC} K^{+}"); - h2TPCnSigma->GetYaxis()->SetTitle("n#sigma_{TPC} K^{-}"); - - registry.add("QATrack/h2TPCnSigmaPt", "", kTH2F, {ptaxis, {200, -10, 10}}); - auto h2TPCnSigmaPt = registry.get(HIST("QATrack/h2TPCnSigmaPt")); - h2TPCnSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - h2TPCnSigmaPt->GetYaxis()->SetTitle("n#sigma_{TPC} K^{#pm}"); - - // TOF PID - registry.add("QATrack/hTOFnSigma", "Distribution of TOF nSigma of K^{+} and K^{-}", kTH1F, {{200, -10, 10}}); + h2dThetaPt->GetYaxis()->SetTitle("#Delta#theta (rad)"); - registry.add("QATrack/h2TOFnSigma", "", kTH2F, {{200, -10, 10}, {200, -10, 10}}); - auto h2TOFnSigma = registry.get(HIST("QATrack/h2TOFnSigma")); - h2TOFnSigma->GetXaxis()->SetTitle("n#sigma_{TOF} K^{+}"); - h2TOFnSigma->GetYaxis()->SetTitle("n#sigma_{TOF} K^{-}"); + // Rotational background QA + if (produceRotational) { + registry.add("QARotational/hRapidity", "Rapidity distribution of #Phi candidates from rotational background", kTH1F, {{200, -1, 1}}); + registry.add("QARotational/hEta", "Pseudorapidity distribution of #Phi candidates from rotational background", kTH1F, {{200, -1, 1}}); + registry.add("QARotational/hdPhi", "Rotational background: Azimuthal distribution (#Delta#phi)", kTH1F, {{100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); + auto hRPhi = registry.get(HIST("QARotational/hdPhi")); + hRPhi->GetXaxis()->SetTitle("#Delta#phi"); - registry.add("QATrack/h2TOFnSigmaPt", "", kTH2F, {ptaxis, {200, -10, 10}}); - auto h2TOFnSigmaPt = registry.get(HIST("QATrack/h2TOFnSigmaPt")); - h2TOFnSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - h2TOFnSigmaPt->GetYaxis()->SetTitle("n#sigma_{TOF} K^{#pm}"); + registry.add("QARotational/h2dPhiPt", "Rotational background: Azimuthal distribution (#Delta#phi) vs p_{T}", kTH2F, {ptaxis, {100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); + auto hR2dPhiPt = registry.get(HIST("QARotational/h2dPhiPt")); + hR2dPhiPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + hR2dPhiPt->GetYaxis()->SetTitle("#Delta#phi"); - // MC Truth - if (static_cast(produce.produceTrue)) { - registry.add("QAMC/Truth/hMCEvent", "MC Truth Event statistics", kTH1F, {{1, 0.0f, 1.0f}}); - auto hMCEventTruth = registry.get(HIST("QAMC/Truth/hMCEvent")); - hMCEventTruth->GetXaxis()->SetBinLabel(1, "Full MC Truth event statistics"); - hMCEventTruth->SetMinimum(0.1); + registry.add("QARotational/hTheta", "Rotational background: Polar distribution (#theta)", kTH1F, {{100, 0.0f, o2::constants::math::PI}}); + auto hRdTheta = registry.get(HIST("QARotational/hTheta")); + hRdTheta->GetXaxis()->SetTitle("#theta (rad)"); - registry.add("QAMC/Truth/hInvMassTrueFalse", "", kTH1F, {invAxis}); // not written events in True distribution due to repetition of mothers?? + registry.add("QARotational/h2dThetaPt", "Rotational background: Polar distribution (#Delta#theta) vs p_{T}", kTH2F, {ptaxis, {100, -o2::constants::math::PI, o2::constants::math::PI}}); + auto hR2dThetaPt = registry.get(HIST("QARotational/h2dThetaPt")); + hR2dThetaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + hR2dThetaPt->GetYaxis()->SetTitle("#Delta#theta"); } + // Mixing QA if (mixingType != rsn::MixingType::none) { registry.add("QAMixing/hSelection", "Event mixing selection statistics", kTH1D, {{1, 0.0f, 1.0f}}); @@ -358,47 +356,90 @@ struct PhianalysisTHnSparse { hEMTvz->GetYaxis()->SetTitle("2.Event V_{z}"); } - if (produceRotational) { - registry.add("QARotational/hSelection", "Rotational background selection statistics", kTH1D, {{1, 0.0f, 1.0f}}); - auto hRB = registry.get(HIST("QARotational/hSelection")); - hRB->GetXaxis()->SetBinLabel(1, "Full rotational background statistics"); - hRB->SetMinimum(0.1); + // PID QA + // TPC + registry.add("QAPID/hTPCnSigma", "Distribution of TPC nSigma of K^{+} and K^{-}", kTH1F, {{200, -10, 10}}); + auto hTPCnSigma = registry.get(HIST("QAPID/hTPCnSigma")); + hTPCnSigma->GetXaxis()->SetTitle("n#sigma_{TPC} K^{#pm}"); - // General rotational distributions - registry.add("QARotational/hRapidity", "Rapidity distribution of #phi candidates from rotational background", kTH1F, {rapidityaxis}); - registry.add("QARotational/hEta", "Pseudorapidity distribution of #phi candidates from rotational background", kTH1F, {etaaxis}); + registry.add("QAPID/h2TPCnSigma", "", kTH2F, {{200, -10, 10}, {200, -10, 10}}); + auto h2TPCnSigma = registry.get(HIST("QAPID/h2TPCnSigma")); + h2TPCnSigma->GetXaxis()->SetTitle("n#sigma_{TPC} K^{+}"); + h2TPCnSigma->GetYaxis()->SetTitle("n#sigma_{TPC} K^{-}"); - // Angular distributions for rotational background - registry.add("QARotational/hdPhi", "Rotational background #Delta#phi distribution", kTH1F, {{100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); - auto hRPhi = registry.get(HIST("QARotational/hdPhi")); - hRPhi->GetXaxis()->SetTitle("#Delta#phi"); - hRPhi->GetYaxis()->SetTitle("Counts"); + registry.add("QAPID/h2TPCnSigmaPt", "", kTH2F, {ptaxis, {200, -10, 10}}); + auto h2TPCnSigmaPt = registry.get(HIST("QAPID/h2TPCnSigmaPt")); + h2TPCnSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + h2TPCnSigmaPt->GetYaxis()->SetTitle("n#sigma_{TPC} K^{#pm}"); - registry.add("QARotational/h2dPhiPt", "Rotational background #Delta#phi vs p_{T}", kTH2F, {ptaxis, {100, -o2::constants::math::TwoPI, o2::constants::math::TwoPI}}); - auto hR2dPhiPt = registry.get(HIST("QARotational/h2dPhiPt")); - hR2dPhiPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - hR2dPhiPt->GetYaxis()->SetTitle("#Delta#phi"); + // TOF + registry.add("QAPID/hTOFnSigma", "Distribution of TOF nSigma of K^{+} and K^{-}", kTH1F, {{200, -10, 10}}); + auto hTOFnSigma = registry.get(HIST("QAPID/hTOFnSigma")); + hTOFnSigma->GetXaxis()->SetTitle("n#sigma_{TOF} K^{#pm}"); - registry.add("QARotational/hTheta", "Rotational background #Delta#theta distribution", kTH1F, {{100, 0.0f, o2::constants::math::PI}}); - auto hRdTheta = registry.get(HIST("QARotational/hTheta")); - hRdTheta->GetXaxis()->SetTitle("#Delta#theta"); - hRdTheta->GetYaxis()->SetTitle("Counts"); + registry.add("QAPID/h2TOFnSigma", "", kTH2F, {{200, -10, 10}, {200, -10, 10}}); + auto h2TOFnSigma = registry.get(HIST("QAPID/h2TOFnSigma")); + h2TOFnSigma->GetXaxis()->SetTitle("n#sigma_{TOF} K^{+}"); + h2TOFnSigma->GetYaxis()->SetTitle("n#sigma_{TOF} K^{-}"); - registry.add("QARotational/h2dThetaPt", "Rotational background #Delta#theta vs p_{T}", kTH2F, {{12, 0, 12}, {100, -o2::constants::math::PI, o2::constants::math::PI}}); - auto hR2dThetaPt = registry.get(HIST("QARotational/h2dThetaPt")); - hR2dThetaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - hR2dThetaPt->GetYaxis()->SetTitle("#Delta#theta"); + registry.add("QAPID/h2TOFnSigmaPt", "", kTH2F, {ptaxis, {200, -10, 10}}); + auto h2TOFnSigmaPt = registry.get(HIST("QAPID/h2TOFnSigmaPt")); + h2TOFnSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + h2TOFnSigmaPt->GetYaxis()->SetTitle("n#sigma_{TOF} K^{#pm}"); + + // MC + if (static_cast(produce.produceTrue)) { + // Rec + registry.add("QAMC/Rec/hSelection", "MC Rec Event statistics", kTH1F, {{5, 0.0f, 5.0f}}); + auto hMCEventTruth = registry.get(HIST("QAMC/Rec/hSelection")); + hMCEventTruth->GetXaxis()->SetBinLabel(1, "Full MC Rec event statistics"); + hMCEventTruth->GetXaxis()->SetBinLabel(2, "MC Rec events passing sel8 cut"); + hMCEventTruth->GetXaxis()->SetBinLabel(3, "MC Rec events passing V_{z} cut"); + hMCEventTruth->GetXaxis()->SetBinLabel(4, "MC Rec events with V_{z} cut and INEL>0"); + hMCEventTruth->GetXaxis()->SetBinLabel(5, "Reconstructed #Phi candidates matched to true #Phi"); + hMCEventTruth->SetMinimum(0.1); + + registry.add("QAMC/Rec/hInvMassTrueFalse", "", kTH1F, {invAxis}); // not written events in True distribution due to repetition of mothers + + // Gen + registry.add("QAMC/Gen/hSelection", "MC Gen Event statistics", kTH1F, {{4, 0.0f, 4.0f}}); + auto hMCEventGen = registry.get(HIST("QAMC/Gen/hSelection")); + hMCEventGen->GetXaxis()->SetBinLabel(1, "Full MC Gen event statistics"); + hMCEventGen->GetXaxis()->SetBinLabel(2, "MC Gen events within V_{z} cut"); + hMCEventGen->GetXaxis()->SetBinLabel(3, "MC Gen events with V_{z} cut and INEL>0"); + hMCEventGen->GetXaxis()->SetBinLabel(4, "Generated #Phi candidates"); + hMCEventGen->SetMinimum(0.1); + + // Resolution + registry.add("Factors/h2ResolutionVz", "Resolution of collision V_{z}", kTH2F, {vzaxis, axisResolutionVz}); + auto hResVz = registry.get(HIST("Factors/h2ResolutionVz")); + hResVz->GetXaxis()->SetTitle("V_{z}^{rec} (cm)"); + hResVz->GetYaxis()->SetTitle("#DeltaV_{z} = V_{z}^{rec} - V_{z}^{gen} (cm)"); + registry.add("Factors/h2ResolutionPt", "Resolution of charged particles p_{T}", kTH2F, {ptaxis, axisResolutionPt}); + auto hResPt = registry.get(HIST("Factors/h2ResolutionPt")); + ; + hResPt->GetXaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); + hResPt->GetYaxis()->SetTitle("#Deltap_{T} = p_{T}^{rec} - p_{T}^{gen} (GeV/c)"); + registry.add("Factors/h2ResolutionPtPhi", "Resolution of Phi p_{T}", kTH2F, {ptaxis, axisResolutionPtPhi}); + auto hResPtPhi = registry.get(HIST("Factors/h2ResolutionPtPhi")); + hResPtPhi->GetXaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); + hResPtPhi->GetYaxis()->SetTitle("#Deltap_{T} = p_{T}^{rec} - p_{T}^{gen} (GeV/c)"); + registry.add("Factors/h2ResolutionMass", "Resolution of mass vs p_{T}^{rec}", kTH2F, {ptaxis, axisResolutionMass}); + auto hResMass = registry.get(HIST("Factors/h2ResolutionMass")); + hResMass->GetXaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); + hResMass->GetYaxis()->SetTitle("#Deltam = m^{rec} - m^{gen} (GeV/c^{2})"); } } - // Event/Signal loss histograms + // Factors registry.add("Factors/hCentralityVsMultMC", "Event centrality vs MC multiplicity", kTH2F, {{101, 0.0f, 101.0f}, axisNch}); + registry.add("Factors/hCentralityVsMult", "Event centrality vs multiplicity", kTH2F, {{101, 0.0f, 101.0f}, axisNch}); registry.add("Factors/hEventCentrality", "Event centrality", kTH1F, {{101, 0, 101}}); - registry.add("Factors/hNrecInGen", "Number of collisions in MC", kTH1F, {{4, -0.5, 3.5}}); + registry.add("Factors/hNrecInGen", "Number of collisions in MC", kTH1F, {{3, -0.5, 2.5}}); registry.add("Factors/hGenEvents", "Generated events", HistType::kTH2F, {{axisNch}, {4, 0, 4}}); auto hGenEvents = registry.get(HIST("Factors/hGenEvents")); hGenEvents->GetYaxis()->SetBinLabel(1, "All generated events"); - hGenEvents->GetYaxis()->SetBinLabel(2, "Generated events with Mc collision V_{z} cut"); - hGenEvents->GetYaxis()->SetBinLabel(3, "Generated events with Mc collision V_{z} cut and INEL>0"); + hGenEvents->GetYaxis()->SetBinLabel(2, "Generated events passing V_{z} cut"); + hGenEvents->GetYaxis()->SetBinLabel(3, "Generated events passing INEL>0"); hGenEvents->GetYaxis()->SetBinLabel(4, "Generated events with at least one reconstructed event"); registry.add("Factors/h2dGenPhi", "Centrality vs p_{T}", kTH2D, {{101, 0.0f, 101.0f}, ptaxis}); registry.add("Factors/h3dGenPhiVsMultMCVsCentrality", "MC multiplicity vs centrality vs p_{T}", kTH3D, {axisNch, {101, 0.0f, 101.0f}, ptaxis}); @@ -516,31 +557,50 @@ struct PhianalysisTHnSparse { { auto posDaughters = positive->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); auto negDaughters = negative->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - n = 0; + + int nch = 0; if (produceQA) registry.fill(HIST("QAEvent/hSelection"), 0.5); + if (!collision.sel8()) + return; + + if (produceQA) + registry.fill(HIST("QAEvent/hSelection"), 1.5); + + if (std::abs(collision.posZ()) > vzCut) + return; + + if (produceQA) + registry.fill(HIST("QAEvent/hSelection"), 2.5); + if (inelGrater0 && !collision.isInelGt0()) return; registry.fill(HIST("Factors/hEventCentrality"), collision.centFT0M()); if (produceQA) { - registry.fill(HIST("QAEvent/hSelection"), 1.5); + registry.fill(HIST("QAEvent/hSelection"), 3.5); registry.fill(HIST("QAEvent/hVtxZ"), collision.posZ()); registry.fill(HIST("QAEvent/hMult"), getMultiplicity(collision)); registry.fill(HIST("QAEvent/hCent"), getCentrality(collision)); + if (produceStats) { dataQA = true; for (const auto& track : posDaughters) { + if (track.isPrimaryTrack() && std::abs(track.eta()) >= static_cast(cut.etatrack)) + nch++; selectedTrack(track, true); } for (const auto& track : negDaughters) { + if (track.isPrimaryTrack() && std::abs(track.eta()) >= static_cast(cut.etatrack)) + nch++; selectedTrack(track, false); } dataQA = false; } + registry.fill(HIST("Factors/hCentralityVsMult"), getCentrality(collision), nch); } if (static_cast(verbose.verboselevel) > 0 && static_cast(verbose.refresh) > 0 && collision.globalIndex() % static_cast(verbose.refresh) == static_cast(verbose.refreshIndex)) @@ -557,20 +617,20 @@ struct PhianalysisTHnSparse { continue; if (produceQA) { - registry.fill(HIST("QATrack/h2TPCnSigma"), track1.tpcNSigmaKa(), track2.tpcNSigmaKa()); - registry.fill(HIST("QATrack/h2TPCnSigmaPt"), track1.pt(), track1.tpcNSigmaKa()); - registry.fill(HIST("QATrack/h2TPCnSigmaPt"), track2.pt(), track2.tpcNSigmaKa()); + registry.fill(HIST("QAPID/h2TPCnSigma"), track1.tpcNSigmaKa(), track2.tpcNSigmaKa()); + registry.fill(HIST("QAPID/h2TPCnSigmaPt"), track1.pt(), track1.tpcNSigmaKa()); + registry.fill(HIST("QAPID/h2TPCnSigmaPt"), track2.pt(), track2.tpcNSigmaKa()); - registry.fill(HIST("QATrack/h2TOFnSigma"), track1.tofNSigmaKa(), track2.tofNSigmaKa()); - registry.fill(HIST("QATrack/h2TOFnSigmaPt"), track1.pt(), track1.tofNSigmaKa()); - registry.fill(HIST("QATrack/h2TOFnSigmaPt"), track2.pt(), track2.tofNSigmaKa()); + registry.fill(HIST("QAPID/h2TOFnSigma"), track1.tofNSigmaKa(), track2.tofNSigmaKa()); + registry.fill(HIST("QAPID/h2TOFnSigmaPt"), track1.pt(), track1.tofNSigmaKa()); + registry.fill(HIST("QAPID/h2TOFnSigmaPt"), track2.pt(), track2.tofNSigmaKa()); - registry.fill(HIST("QATrack/hTPCnSigma"), track1.tpcNSigmaKa()); - registry.fill(HIST("QATrack/hTPCnSigma"), track2.tpcNSigmaKa()); + registry.fill(HIST("QAPID/hTPCnSigma"), track1.tpcNSigmaKa()); + registry.fill(HIST("QAPID/hTPCnSigma"), track2.tpcNSigmaKa()); if (track1.hasTOF()) - registry.fill(HIST("QATrack/hTOFnSigma"), track1.tofNSigmaKa()); + registry.fill(HIST("QAPID/hTOFnSigma"), track1.tofNSigmaKa()); if (track2.hasTOF()) - registry.fill(HIST("QATrack/hTOFnSigma"), track2.tofNSigmaKa()); + registry.fill(HIST("QAPID/hTOFnSigma"), track2.tofNSigmaKa()); registry.fill(HIST("QATrack/hEta"), track1.eta()); registry.fill(HIST("QATrack/hEta"), track2.eta()); @@ -607,29 +667,20 @@ struct PhianalysisTHnSparse { 0); rsnOutput->fillUnlikepm(pointPair); - if (produceQA) { - registry.fill(HIST("QAEvent/hSelection"), 3.5); - if (n == 0) - registry.fill(HIST("QAEvent/hSelection"), 2.5); - } - n = n + 1; - if (produceRotational) { for (int i = 1; i <= static_cast(numberofRotations); i++) { - // compute rotation angle in radians using o2::constants::math::PI - float angleDeg = i * (360.0f / (static_cast(numberofRotations) + 1)); - float angleRad = angleDeg * (o2::constants::math::PI / 180.0f); - float px2new = track2.px() * std::cos(angleRad) - track2.py() * std::sin(angleRad); - float py2new = track2.px() * std::sin(angleRad) + track2.py() * std::cos(angleRad); + float starting = static_cast(startingAngle) * TMath::DegToRad(); + float angle = starting + i * ((o2::constants::math::TwoPI - 2 * starting) / (static_cast(numberofRotations) + 1)); + float px2new = track2.px() * std::cos(angle) - track2.py() * std::sin(angle); + float py2new = track2.px() * std::sin(angle) + track2.py() * std::cos(angle); d2 = ROOT::Math::PxPyPzMVector(px2new, py2new, track2.pz(), massNeg); mother = d1 + d2; if (produceQA) { - registry.fill(HIST("QARotational/hSelection"), 0.5); registry.fill(HIST("QARotational/hRapidity"), mother.Rapidity()); registry.fill(HIST("QARotational/hEta"), mother.Eta()); - registry.fill(HIST("QARotational/hdPhi"), track1.phi() - track2.phi()); - registry.fill(HIST("QARotational/h2dPhiPt"), mother.Pt(), track1.phi() - track2.phi()); + registry.fill(HIST("QARotational/hdPhi"), d1.Phi() - d2.Phi()); + registry.fill(HIST("QARotational/h2dPhiPt"), mother.Pt(), d1.Phi() - d2.Phi()); registry.fill(HIST("QARotational/hTheta"), mother.Theta()); registry.fill(HIST("QARotational/h2dThetaPt"), mother.Pt(), d1.Theta() - d2.Theta()); } @@ -712,12 +763,18 @@ struct PhianalysisTHnSparse { } PROCESS_SWITCH(PhianalysisTHnSparse, processData, "Process Event for Data", true); - void processTrue(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& /*tracks*/, aod::McParticles const& /*mcParticles*/, aod::McCollisions const& /*mcCollisions*/) + void processTrue(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const& /*mcParticles*/, aod::McCollisions const& /*mcCollisions*/) { if (!static_cast(produce.produceTrue)) return; - registry.fill(HIST("QAMC/Truth/hMCEvent"), 0.5); + registry.fill(HIST("QAMC/Rec/hSelection"), 0.5); + + if (!collision.sel8()) + return; + + if (produceQA) + registry.fill(HIST("QAMC/Rec/hSelection"), 1.5); auto posDaughtersMC = positiveMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); auto negDaughtersMC = negativeMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); @@ -728,11 +785,27 @@ struct PhianalysisTHnSparse { return; } auto mcCollision = collision.mcCollision(); + registry.fill(HIST("Factors/h2ResolutionVz"), collision.posZ(), (collision.posZ() - mcCollision.posZ())); + if (std::abs(mcCollision.posZ()) > vzCut) return; + + if (produceQA) + registry.fill(HIST("QAMC/Rec/hSelection"), 2.5); + if (inelGrater0 && !collision.isInelGt0()) return; + if (produceQA) + registry.fill(HIST("QAMC/Rec/hSelection"), 3.5); + + for (const auto& track : tracks) { + if (track.has_mcParticle()) { + auto mctrack = track.mcParticle(); + registry.fill(HIST("Factors/h2ResolutionPt"), track.pt(), (track.pt() - mctrack.pt())); + } + } + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersMC, negDaughtersMC))) { if (!track1.has_mcParticle()) { @@ -786,7 +859,7 @@ struct PhianalysisTHnSparse { if (n > 0) { if (produceQA) - registry.fill(HIST("QAMC/Truth/hInvMassTrueFalse"), mother.M()); + registry.fill(HIST("QAMC/Rec/hInvMassTrueFalse"), mother.M()); continue; } @@ -803,6 +876,16 @@ struct PhianalysisTHnSparse { 0, 0); + if (produceQA) + registry.fill(HIST("QAMC/Rec/hSelection"), 4.5); + + auto phiP = mothertrack1.p(); + auto phiE = mothertrack1.e(); + auto phiMass = std::sqrt(phiE * phiE - phiP * phiP); + + registry.fill(HIST("Factors/h2ResolutionPtPhi"), mother.pt(), (mother.pt() - mothertrack1.pt())); + registry.fill(HIST("Factors/h2ResolutionMass"), mother.Pt(), (mother.M() - phiMass)); + rsnOutput->fillUnliketrue(pointPair); n++; } @@ -813,12 +896,24 @@ struct PhianalysisTHnSparse { void processGen(McCollisionMults::iterator const& mcCollision, soa::SmallGroups const& collisions, LabeledTracks const& /*particles*/, aod::McParticles const& mcParticles) { + if (!static_cast(produce.produceTrue)) + return; + + if (produceQA) + registry.fill(HIST("QAMC/Gen/hSelection"), 0.5); + if (std::abs(mcCollision.posZ()) > vzCut) return; + if (produceQA) + registry.fill(HIST("QAMC/Gen/hSelection"), 1.5); + if (inelGrater0 && !mcCollision.isInelGt0()) return; + if (produceQA) + registry.fill(HIST("QAMC/Gen/hSelection"), 2.5); + if (collisions.size() == 0) return; @@ -868,13 +963,15 @@ struct PhianalysisTHnSparse { 0); rsnOutput->fillUnlikegen(pointPair); + if (produceQA) + registry.fill(HIST("QAMC/Gen/hSelection"), 3.5); } } } } PROCESS_SWITCH(PhianalysisTHnSparse, processGen, "Process MC Generated.", false); - void processMixed(EventCandidates const& collisions, TrackCandidates const& tracks) + void processMixed(soa::Filtered const& collisions, TrackCandidates const& tracks) { if (mixingType == rsn::MixingType::none) return;