diff --git a/PWGCF/Flow/TableProducer/zdcQVectors.cxx b/PWGCF/Flow/TableProducer/zdcQVectors.cxx index c646b410b5e..07befdd4d31 100644 --- a/PWGCF/Flow/TableProducer/zdcQVectors.cxx +++ b/PWGCF/Flow/TableProducer/zdcQVectors.cxx @@ -20,6 +20,8 @@ #include "Common/CCDB/RCTSelectionFlags.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/Multiplicity.h" #include #include @@ -155,12 +157,25 @@ struct ZdcQVectors { O2_DEFINE_CONFIGURABLE(cfgFillNothing, bool, false, "Disable ALL Histograms -> ONLY use to reduce memory"); O2_DEFINE_CONFIGURABLE(cfgNoGain, bool, false, "Do not apply gain correction to ZDC energy calibration"); + O2_DEFINE_CONFIGURABLE(cfgTrackSelsDCAxy, float, 0.2, "Cut on DCA in the transverse direction (cm)"); + O2_DEFINE_CONFIGURABLE(cfgTrackSelsDCAz, float, 0.2, "Cut on DCA in the longitudinal direction (cm)"); + O2_DEFINE_CONFIGURABLE(cfgTrackSelsPtmin, float, 0.2, "minimum pt (GeV/c)"); + O2_DEFINE_CONFIGURABLE(cfgTrackSelsPtmax, float, 10, "maximum pt (GeV/c)"); + O2_DEFINE_CONFIGURABLE(cfgTrackSelsEta, float, 0.8, "eta cut"); + O2_DEFINE_CONFIGURABLE(cfgCCDBdir_Shift, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/Shift", "CCDB directory for Shift ZDC"); + Configurable> cfgSelVec{"cfgSelVec", std::vector{1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1}, "Put 1 for every event selection from SelectionCriteria that is used in flowSP"}; + Configurable> cfgEvSelsMultPv{"cfgEvSelsMultPv", std::vector{2223.49, -75.1444, 0.963572, -0.00570399, 1.34877e-05, 3790.99, -137.064, 2.13044, -0.017122, 5.82834e-05}, "Multiplicity cuts (PV) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "}; + Configurable> cfgEvSelsMult{"cfgEvSelsMult", std::vector{1301.56, -41.4615, 0.478224, -0.00239449, 4.46966e-06, 2967.6, -102.927, 1.47488, -0.0106534, 3.28622e-05}, "Multiplicity cuts (Global) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "}; // define my..... // Filter collisionFilter = nabs(aod::collision::posZ) <; - using UsedCollisions = soa::Join; + + using UsedCollisions = soa::Join; using BCsRun3 = soa::Join; + Filter trackFilter = nabs(aod::track::eta) < cfgTrackSelsEta && aod::track::pt > cfgTrackSelsPtmin && aod::track::pt < cfgTrackSelsPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && nabs(aod::track::dcaXY) < cfgTrackSelsDCAxy && nabs(aod::track::dcaZ) < cfgTrackSelsDCAz; + using UnfilteredTracks = soa::Join; + using UsedTracks = soa::Filtered; enum SelectionCriteria { evSel_FilteredEvent, @@ -178,6 +193,7 @@ struct ZdcQVectors { evSel_kIsGoodITSLayer0123, evSel_RCTFlagsZDC, evSel_CentCuts, + evSel_MultCut, nEventSelections }; @@ -191,6 +207,15 @@ struct ZdcQVectors { // Define output HistogramRegistry registry{"Registry"}; + // Event selection cuts + std::unique_ptr fPhiCutLow = nullptr; + std::unique_ptr fPhiCutHigh = nullptr; + std::unique_ptr fMultPVCutLow = nullptr; + std::unique_ptr fMultPVCutHigh = nullptr; + std::unique_ptr fMultCutLow = nullptr; + std::unique_ptr fMultCutHigh = nullptr; + std::unique_ptr fMultMultPVCut = nullptr; + Service ccdb; // keep track of calibration histos for each given step and iteration @@ -246,6 +271,7 @@ struct ZdcQVectors { registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_CentCuts + 1, "Cenrality range"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_kIsGoodITSLayersAll + 1, "kIsGoodITSLayersAll"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_kIsGoodITSLayer0123 + 1, "kIsGoodITSLayer0123"); + registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_MultCut + 1, "Mult & MultPV"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_isSelectedZDC + 1, "isSelected"); int totalTowers = 10; @@ -380,7 +406,28 @@ struct ZdcQVectors { registry.addClone("recentering/before/", "recentering/after/"); registry.addClone("QA/before/", "QA/after/"); } + + fMultPVCutLow = std::make_unique("fMultPVCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100); + fMultPVCutHigh = std::make_unique("fMultPVCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100); + fMultCutLow = std::make_unique("fMultCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100); + fMultCutHigh = std::make_unique("fMultCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100); + + std::vector paramsMultPVCut = cfgEvSelsMultPv; + std::vector paramsMultCut = cfgEvSelsMult; + + uint64_t nParams = 10; + + if (paramsMultPVCut.size() < nParams) { + LOGF(fatal, "cfg.cEvSelsMultPv not set properly.. size = %d (should be 10) --> Check your config files!", paramsMultPVCut.size()); + } else if (paramsMultCut.size() < nParams) { + LOGF(fatal, "cfg.cEvSelsMult not set properly.. size = %d (should be 10) --> Check your config files!", paramsMultCut.size()); + } else { + fMultPVCutLow->SetParameters(paramsMultPVCut[0], paramsMultPVCut[1], paramsMultPVCut[2], paramsMultPVCut[3], paramsMultPVCut[4]); + fMultPVCutHigh->SetParameters(paramsMultPVCut[5], paramsMultPVCut[6], paramsMultPVCut[7], paramsMultPVCut[8], paramsMultPVCut[9]); + fMultCutLow->SetParameters(paramsMultCut[0], paramsMultCut[1], paramsMultCut[2], paramsMultCut[3], paramsMultCut[4]); + fMultCutHigh->SetParameters(paramsMultCut[5], paramsMultCut[6], paramsMultCut[7], paramsMultCut[8], paramsMultCut[9]); } +} double rescaleTimestamp(uint64_t timestamp, int runnumber) { @@ -443,7 +490,7 @@ struct ZdcQVectors { } template - uint16_t eventSelected(TCollision collision, TBunchCrossing bunchCrossing) + uint16_t eventSelected(TCollision collision, TBunchCrossing bunchCrossing, bool& isEventSelected, const int& multTrk) { uint16_t selectionBits = 0; bool selected; @@ -457,12 +504,16 @@ struct ZdcQVectors { if (selected) { selectionBits |= static_cast(0x1u << evSel_Zvtx); fillCutAnalysis(collision, bunchCrossing, evSel_Zvtx); + } else if(cfgSelVec.value[evSel_Zvtx]){ + isEventSelected = false; } selected = collision.sel8(); if (selected) { selectionBits |= static_cast(0x1u << evSel_sel8); fillCutAnalysis(collision, bunchCrossing, evSel_sel8); + } else if(cfgSelVec.value[evSel_sel8]){ + isEventSelected = false; } auto occupancy = collision.trackOccupancyInTimeRange(); @@ -470,59 +521,107 @@ struct ZdcQVectors { if (selected) { selectionBits |= static_cast(0x1u << evSel_occupancy); fillCutAnalysis(collision, bunchCrossing, evSel_occupancy); + } else if(cfgSelVec.value[evSel_occupancy]){ + isEventSelected = false; } selected = collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup); if (selected) { selectionBits |= static_cast(0x1u << evSel_kNoSameBunchPileup); fillCutAnalysis(collision, bunchCrossing, evSel_kNoSameBunchPileup); + } else if(cfgSelVec.value[evSel_kNoSameBunchPileup]){ + isEventSelected = false; } selected = collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV); if (selected) { selectionBits |= static_cast(0x1u << evSel_kIsGoodZvtxFT0vsPV); fillCutAnalysis(collision, bunchCrossing, evSel_kIsGoodZvtxFT0vsPV); + } else if(cfgSelVec.value[evSel_kIsGoodZvtxFT0vsPV]){ + isEventSelected = false; } selected = collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard); if (selected) { selectionBits |= static_cast(0x1u << evSel_kNoCollInTimeRangeStandard); fillCutAnalysis(collision, bunchCrossing, evSel_kNoCollInTimeRangeStandard); + } else if(cfgSelVec.value[evSel_kNoCollInTimeRangeStandard]){ + isEventSelected = false; } selected = collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow); if (selected) { selectionBits |= static_cast(0x1u << evSel_kNoCollInTimeRangeNarrow); fillCutAnalysis(collision, bunchCrossing, evSel_kNoCollInTimeRangeNarrow); + } else if(cfgSelVec.value[evSel_kNoCollInTimeRangeNarrow]){ + isEventSelected = false; } selected = collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC); if (selected) { selectionBits |= static_cast(0x1u << evSel_kIsVertexITSTPC); fillCutAnalysis(collision, bunchCrossing, evSel_kIsVertexITSTPC); + } else if(cfgSelVec.value[evSel_kIsVertexITSTPC]){ + isEventSelected = false; } selected = collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll); if (selected) { selectionBits |= static_cast(0x1u << evSel_kIsGoodITSLayersAll); fillCutAnalysis(collision, bunchCrossing, evSel_kIsGoodITSLayersAll); + } else if(cfgSelVec.value[evSel_kIsGoodITSLayersAll]){ + isEventSelected = false; } selected = collision.selection_bit(o2::aod::evsel::kIsGoodITSLayer0123); if (selected) { selectionBits |= static_cast(0x1u << evSel_kIsGoodITSLayer0123); fillCutAnalysis(collision, bunchCrossing, evSel_kIsGoodITSLayer0123); + } else if(cfgSelVec.value[evSel_kIsGoodITSLayer0123]){ + isEventSelected = false; } selected = rctChecker(collision); if (selected) { selectionBits |= static_cast(0x1u << evSel_RCTFlagsZDC); fillCutAnalysis(collision, bunchCrossing, evSel_RCTFlagsZDC); + } else if(cfgSelVec.value[evSel_RCTFlagsZDC]){ + isEventSelected = false; + } + + float vtxz = -999; + if (collision.numContrib() > 1) { + vtxz = collision.posZ(); + float zRes = std::sqrt(collision.covZZ()); + float minzRes = 0.25; + int maxNumContrib = 20; + if (zRes > minzRes && collision.numContrib() < maxNumContrib) + vtxz = -999; + } + + auto multNTracksPV = collision.multNTracksPV(); + selected = true; + + if (vtxz > cfgVtxZ || vtxz < -cfgVtxZ) + selected = false; + if (multNTracksPV < fMultPVCutLow->Eval(collision.centFT0C())) + selected = false; + if (multNTracksPV > fMultPVCutHigh->Eval(collision.centFT0C())) + selected = false; + if (multTrk < fMultCutLow->Eval(collision.centFT0C())) + selected = false; + if (multTrk > fMultCutHigh->Eval(collision.centFT0C())) + selected = false; + + if (selected) { + selectionBits |= static_cast(0x1u << evSel_MultCut); + fillCutAnalysis(collision, bunchCrossing, evSel_MultCut); + } else if(cfgSelVec.value[evSel_MultCut]){ + isEventSelected = false; } // Fill for centrality estimators! fillCutAnalysis(collision, bunchCrossing, nEventSelections); - return selectionBits; } @@ -710,7 +809,8 @@ struct ZdcQVectors { void process(UsedCollisions::iterator const& collision, BCsRun3 const& /*bcs*/, - aod::Zdcs const& /*zdcs*/) + aod::Zdcs const& /*zdcs*/, + UsedTracks const& tracks) { // for Q-vector calculation // A[0] & C[1] @@ -807,7 +907,11 @@ struct ZdcQVectors { } registry.fill(HIST("hEventCount"), evSel_isSelectedZDC); - uint16_t eventSelectionFlags = eventSelected(collision, foundBC.zdc()); + // Use this bool to check for given set of event selections of Q-vectors would be selected + // Enable plotting only if event would be selected + bool isEventSelected = true; + + uint16_t eventSelectionFlags = eventSelected(collision, foundBC.zdc(), isEventSelected, tracks.size()); // ALWAYS use these event selections if (cent < EvSel.cfgCentMin || cent > EvSel.cfgCentMax || !collision.sel8() || std::abs(collision.posZ()) > cfgVtxZ) { @@ -842,10 +946,11 @@ struct ZdcQVectors { // load the calibration histos for iteration 0 step 0 (Energy Calibration) if (!cfgNoGain) loadCalibrations(cfgEnergyCal.value); + // load the calibrations for the mean v loadCalibrations(cfgMeanv.value); - if (!cfgFillNothing) { + if (!cfgFillNothing && isEventSelected) { registry.get(HIST("vmean/hvertex_vx"))->Fill(Form("%d", runnumber), v[0]); registry.get(HIST("vmean/hvertex_vy"))->Fill(Form("%d", runnumber), v[1]); registry.get(HIST("vmean/hvertex_vz"))->Fill(Form("%d", runnumber), v[2]); @@ -891,7 +996,7 @@ struct ZdcQVectors { calibtower++; } - if (cfgFillHistRegistry && !cfgFillNothing) { + if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) { for (int i = 0; i < nTowersPerSide; i++) { float bincenter = i + .5; registry.fill(HIST("QA/ZNA_Energy"), bincenter, eZN[i]); @@ -950,7 +1055,7 @@ struct ZdcQVectors { } } - if (cfgFillHistRegistry && !cfgFillNothing) { + if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) { registry.get(HIST("QA/before/ZNA_Qx"))->Fill(Form("%d", runnumber), q[0]); registry.get(HIST("QA/before/ZNA_Qy"))->Fill(Form("%d", runnumber), q[1]); registry.get(HIST("QA/before/ZNC_Qx"))->Fill(Form("%d", runnumber), q[2]); @@ -974,7 +1079,7 @@ struct ZdcQVectors { std::vector qRec(q); if (cal.atIteration == 0) { - if (isSelected && cfgFillHistRegistry) + if (isSelected && cfgFillHistRegistry && isEventSelected) fillCommonRegistry(q[0], q[1], q[2], q[3], v, centrality, rsTimestamp); spTableZDC(runnumber, cents, v, foundBC.timestamp(), q[0], q[1], q[2], q[3], isSelected, eventSelectionFlags); @@ -982,7 +1087,7 @@ struct ZdcQVectors { lastRunNumber = runnumber; return; } else { - if (cfgFillHistRegistry) + if (cfgFillHistRegistry && isEventSelected) fillCommonRegistry(q[0], q[1], q[2], q[3], v, centrality, rsTimestamp); // vector of 4 @@ -1002,7 +1107,7 @@ struct ZdcQVectors { corrQxC.push_back(getCorrection(names[0][2].Data(), it, 1)); corrQyC.push_back(getCorrection(names[0][3].Data(), it, 1)); - if (cfgFillHistRegistry && !cfgFillNothing) { + if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) { registry.get(HIST("recentering/QXA_vs_iteration"))->Fill(pb + 1, q[0] - std::accumulate(corrQxA.begin(), corrQxA.end(), 0.0)); registry.get(HIST("recentering/QYA_vs_iteration"))->Fill(pb + 1, q[1] - std::accumulate(corrQyA.begin(), corrQyA.end(), 0.0)); registry.get(HIST("recentering/QXC_vs_iteration"))->Fill(pb + 1, q[2] - std::accumulate(corrQxC.begin(), corrQxC.end(), 0.0)); @@ -1016,7 +1121,7 @@ struct ZdcQVectors { corrQxC.push_back(getCorrection(names[step - 1][2].Data(), it, step)); corrQyC.push_back(getCorrection(names[step - 1][3].Data(), it, step)); - if (cfgFillHistRegistry && !cfgFillNothing) { + if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) { registry.get(HIST("recentering/QXA_vs_iteration"))->Fill(pb + 1, q[0] - std::accumulate(corrQxA.begin(), corrQxA.end(), 0.0)); registry.get(HIST("recentering/QYA_vs_iteration"))->Fill(pb + 1, q[1] - std::accumulate(corrQyA.begin(), corrQyA.end(), 0.0)); registry.get(HIST("recentering/QXC_vs_iteration"))->Fill(pb + 1, q[2] - std::accumulate(corrQxC.begin(), corrQxC.end(), 0.0)); @@ -1032,7 +1137,7 @@ struct ZdcQVectors { corrQxC.push_back(getCorrection(namesTS[2].Data(), it, 6)); corrQyC.push_back(getCorrection(namesTS[3].Data(), it, 6)); - if (cfgFillHistRegistry && !cfgFillNothing) { + if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) { registry.get(HIST("recentering/QXA_vs_iteration"))->Fill(pb + 1, q[0] - std::accumulate(corrQxA.begin(), corrQxA.end(), 0.0)); registry.get(HIST("recentering/QYA_vs_iteration"))->Fill(pb + 1, q[1] - std::accumulate(corrQyA.begin(), corrQyA.end(), 0.0)); registry.get(HIST("recentering/QXC_vs_iteration"))->Fill(pb + 1, q[2] - std::accumulate(corrQxC.begin(), corrQxC.end(), 0.0)); @@ -1079,7 +1184,7 @@ struct ZdcQVectors { } for (int ishift = 1; ishift <= nshift; ishift++) { - if (!cfgFillNothing) { + if (!cfgFillNothing && isEventSelected) { registry.fill(HIST("shift/ShiftZDCC"), centrality, 0.5, ishift - 0.5, std::sin(ishift * 1.0 * psiZDCC)); registry.fill(HIST("shift/ShiftZDCC"), centrality, 1.5, ishift - 0.5, std::cos(ishift * 1.0 * psiZDCC)); registry.fill(HIST("shift/ShiftZDCA"), centrality, 0.5, ishift - 0.5, std::sin(ishift * 1.0 * psiZDCA)); @@ -1120,7 +1225,7 @@ struct ZdcQVectors { psiZDCCshift = std::atan2(std::sin(psiZDCCshift), std::cos(psiZDCCshift)); psiZDCAshift = std::atan2(std::sin(psiZDCAshift), std::cos(psiZDCAshift)); - if (cfgFillHistRegistry && !cfgFillNothing) { + if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) { registry.fill(HIST("QA/shift/psiZDCA"), psiZDCA, centrality); registry.fill(HIST("QA/shift/psiZDCC"), psiZDCC, centrality); registry.fill(HIST("QA/shift/psiZDCAC"), psiZDCA, psiZDCC); @@ -1137,7 +1242,7 @@ struct ZdcQVectors { double qXcShift = std::hypot(qRec[2], qRec[3]) * std::cos(psiZDCCshift); double qYcShift = std::hypot(qRec[2], qRec[3]) * std::sin(psiZDCCshift); - if (isSelected && cfgFillHistRegistry && !cfgFillNothing) { + if (isSelected && cfgFillHistRegistry && !cfgFillNothing && isEventSelected) { fillCommonRegistry(qXaShift, qYaShift, qXcShift, qYcShift, v, centrality, rsTimestamp); registry.fill(HIST("QA/centrality_after"), centrality); registry.get(HIST("QA/after/ZNA_Qx"))->Fill(Form("%d", runnumber), qXaShift);