diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 3bbde1dd281..3c37df50324 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1668,20 +1668,11 @@ void VarManager::FillEvent(T const& event, float* values) } if constexpr ((fillMap & CollisionMult) > 0 || (fillMap & ReducedEventExtended) > 0) { - if constexpr ((fillMap & RapidityGapFilter) > 0) { - // UPC: Use the FIT signals from the nearest BC with FIT amplitude above threshold - values[kMultFV0A] = event.newBcMultFV0A(); - values[kMultFT0A] = event.newBcMultFT0A(); - values[kMultFT0C] = event.newBcMultFT0C(); - values[kMultFDDA] = event.newBcMultFDDA(); - values[kMultFDDC] = event.newBcMultFDDC(); - } else { - values[kMultFV0A] = event.multFV0A(); - values[kMultFT0A] = event.multFT0A(); - values[kMultFT0C] = event.multFT0C(); - values[kMultFDDA] = event.multFDDA(); - values[kMultFDDC] = event.multFDDC(); - } + values[kMultFV0A] = event.multFV0A(); + values[kMultFT0A] = event.multFT0A(); + values[kMultFT0C] = event.multFT0C(); + values[kMultFDDA] = event.multFDDA(); + values[kMultFDDC] = event.multFDDC(); values[kMultTPC] = event.multTPC(); values[kMultFV0C] = event.multFV0C(); values[kMultZNA] = event.multZNA(); diff --git a/PWGDQ/DataModel/ReducedInfoTables.h b/PWGDQ/DataModel/ReducedInfoTables.h index 8fc3931d97b..ce81e96aebb 100644 --- a/PWGDQ/DataModel/ReducedInfoTables.h +++ b/PWGDQ/DataModel/ReducedInfoTables.h @@ -36,11 +36,7 @@ namespace o2::aod namespace dqppfilter { DECLARE_SOA_COLUMN(EventFilter, eventFilter, uint64_t); //! Bit-field used for the high level event triggering -DECLARE_SOA_COLUMN(NewBcMultFT0A, newBcMultFT0A, float); //! sum of amplitudes on A side of FT0 -DECLARE_SOA_COLUMN(NewBcMultFT0C, newBcMultFT0C, float); //! sum of amplitudes on C side of FT0 -DECLARE_SOA_COLUMN(NewBcMultFDDA, newBcMultFDDA, float); //! sum of amplitudes on A side of FDD -DECLARE_SOA_COLUMN(NewBcMultFDDC, newBcMultFDDC, float); //! sum of amplitudes on C side of FDD -DECLARE_SOA_COLUMN(NewBcMultFV0A, newBcMultFV0A, float); //! sum of amplitudes on A side of FDD +DECLARE_SOA_COLUMN(NewBcIndex, newBcIndex, uint64_t); //! globalIndex of the new BC determined in filterPbPb } // namespace dqppfilter DECLARE_SOA_TABLE(DQEventFilter, "AOD", "EVENTFILTER", //! Store event-level decisions (DQ high level triggers) @@ -48,19 +44,7 @@ DECLARE_SOA_TABLE(DQEventFilter, "AOD", "EVENTFILTER", //! Store event-level dec DECLARE_SOA_TABLE(DQRapidityGapFilter, "AOD", "RAPIDITYGAPFILTER", dqppfilter::EventFilter, - dqppfilter::NewBcMultFT0A, - dqppfilter::NewBcMultFT0C, - dqppfilter::NewBcMultFDDA, - dqppfilter::NewBcMultFDDC, - dqppfilter::NewBcMultFV0A, - zdc::EnergyCommonZNA, - zdc::EnergyCommonZNC, - zdc::EnergyCommonZPA, - zdc::EnergyCommonZPC, - zdc::TimeZNA, - zdc::TimeZNC, - zdc::TimeZPA, - zdc::TimeZPC); + dqppfilter::NewBcIndex); namespace reducedevent { diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index 9bdc98f2485..a0f39c9e34b 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -598,13 +598,6 @@ struct TableMakerMC { } (reinterpret_cast(fStatsList->At(0)))->Fill(1.0, static_cast(o2::aod::evsel::kNsel)); - // apply the event filter - if constexpr ((TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { - if (!collision.eventFilter()) { - continue; - } - } - auto bc = collision.template bc_as(); // store the selection decisions uint64_t tag = static_cast(0); @@ -664,20 +657,11 @@ struct TableMakerMC { multZNC = collision.multZNC(); multTracklets = collision.multTracklets(); multTracksPV = collision.multNTracksPV(); - if constexpr ((TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { - // Use the FIT signals from the nearest BC with FIT amplitude above threshold - multFV0A = collision.newBcMultFV0A(); - multFT0A = collision.newBcMultFT0A(); - multFT0C = collision.newBcMultFT0C(); - multFDDA = collision.newBcMultFDDA(); - multFDDC = collision.newBcMultFDDC(); - } else { - multFV0A = collision.multFV0A(); - multFT0A = collision.multFT0A(); - multFT0C = collision.multFT0C(); - multFDDA = collision.multFDDA(); - multFDDC = collision.multFDDC(); - } + multFV0A = collision.multFV0A(); + multFT0A = collision.multFT0A(); + multFT0C = collision.multFT0C(); + multFDDA = collision.multFDDA(); + multFDDC = collision.multFDDC(); } if constexpr ((TEventFillMap & VarManager::ObjTypes::CollisionCent) > 0) { centFT0C = collision.centFT0C(); diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index 1ca75734902..19a11da718e 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -331,6 +331,10 @@ struct TableMaker { Partition tracksPosWithCovNoTOF = (((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)) && (aod::track::tgl > static_cast(0.05))); Partition tracksNegWithCovNoTOF = (((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)) && (aod::track::tgl < static_cast(-0.05))); + Preslice presliceWithCov = aod::track::collisionId; + Partition tracksPosWithCov = (((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)) && (aod::track::tgl > static_cast(0.05))); + Partition tracksNegWithCov = (((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)) && (aod::track::tgl < static_cast(-0.05))); + struct { std::map oMeanTimeShortA; std::map oMeanTimeShortC; @@ -788,7 +792,7 @@ struct TableMaker { } template - void skimCollisions(TEvents const& collisions, TBCs const& /*bcs*/, TZdcs const& /*zdcs*/, + void skimCollisions(TEvents const& collisions, TBCs const& bcs, TZdcs const& /*zdcs*/, TTrackAssoc const& trackAssocs, TTracks const& tracks) { // Skim collisions @@ -821,7 +825,7 @@ struct TableMaker { (reinterpret_cast(fStatsList->At(kStatsEvent)))->Fill(1.0, static_cast(o2::aod::evsel::kNsel)); // apply the event filter computed by filter-PP - if constexpr ((TEventFillMap & VarManager::ObjTypes::EventFilter) > 0 || (TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { + if constexpr ((TEventFillMap & VarManager::ObjTypes::EventFilter) > 0) { if (!collision.eventFilter()) { continue; } @@ -836,8 +840,8 @@ struct TableMaker { if (bcEvSel.globalIndex() != bc.globalIndex()) { tag |= (static_cast(1) << 0); } - // Put the 8 first bits of the event filter in the last 8 bits of the tag - if constexpr ((TEventFillMap & VarManager::ObjTypes::EventFilter) > 0 || (TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { + // Put the 8 first bits of the rapidity gap filter in the last 8 bits of the tag + if constexpr ((TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { tag |= (collision.eventFilter() << 56); } @@ -846,8 +850,12 @@ struct TableMaker { VarManager::FillEvent(collision); // extract event information and place it in the fValues array if constexpr ((TEventFillMap & VarManager::ObjTypes::Zdc) > 0) { if constexpr ((TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { - // Collision table already has ZDC info - VarManager::FillZDC(collision); + // The DQRapidityGapFilter contains the index of the bc we want to get ZDC info from + auto newbc = bcs.rawIteratorAt(collision.newBcIndex()); + if (newbc.has_zdc()) { + auto newbc_zdc = newbc.zdc(); + VarManager::FillZDC(newbc_zdc); + } } else if (bcEvSel.has_zdc()) { auto bc_zdc = bcEvSel.zdc(); VarManager::FillZDC(bc_zdc); @@ -926,20 +934,11 @@ struct TableMaker { multZNC = collision.multZNC(); multTracklets = collision.multTracklets(); multTracksPV = collision.multNTracksPV(); - if constexpr ((TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { - // Use the FIT signals from the nearest BC with FIT amplitude above threshold - multFV0A = collision.newBcMultFV0A(); - multFT0A = collision.newBcMultFT0A(); - multFT0C = collision.newBcMultFT0C(); - multFDDA = collision.newBcMultFDDA(); - multFDDC = collision.newBcMultFDDC(); - } else { - multFV0A = collision.multFV0A(); - multFT0A = collision.multFT0A(); - multFT0C = collision.multFT0C(); - multFDDA = collision.multFDDA(); - multFDDC = collision.multFDDC(); - } + multFV0A = collision.multFV0A(); + multFT0A = collision.multFT0A(); + multFT0C = collision.multFT0C(); + multFDDA = collision.multFDDA(); + multFDDC = collision.multFDDC(); } if constexpr ((TEventFillMap & VarManager::ObjTypes::CollisionCent) > 0) { centFT0C = collision.centFT0C(); @@ -952,9 +951,15 @@ struct TableMaker { eventInfo(collision.globalIndex()); if constexpr ((TEventFillMap & VarManager::ObjTypes::Zdc) > 0) { if constexpr ((TEventFillMap & VarManager::ObjTypes::RapidityGapFilter) > 0) { - // ZDC information is already in the DQRapidityGapFilter - zdc(collision.energyCommonZNA(), collision.energyCommonZNC(), collision.energyCommonZPA(), collision.energyCommonZPC(), - collision.timeZNA(), collision.timeZNC(), collision.timeZPA(), collision.timeZPC()); + // The DQRapidityGapFilter contains the index of the bc we want to get ZDC info from + auto newbc = bcs.rawIteratorAt(collision.newBcIndex()); + if (newbc.has_zdc()) { + auto newbc_zdc = newbc.zdc(); + zdc(newbc_zdc.energyCommonZNA(), newbc_zdc.energyCommonZNC(), newbc_zdc.energyCommonZPA(), newbc_zdc.energyCommonZPC(), + newbc_zdc.timeZNA(), newbc_zdc.timeZNC(), newbc_zdc.timeZPA(), newbc_zdc.timeZPC()); + } else { + zdc(-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0); + } } else if (bcEvSel.has_zdc()) { auto bc_zdc = bcEvSel.zdc(); zdc(bc_zdc.energyCommonZNA(), bc_zdc.energyCommonZNC(), bc_zdc.energyCommonZPA(), bc_zdc.energyCommonZPC(), @@ -1650,6 +1655,7 @@ struct TableMaker { MyBarrelTracksWithCov const& tracksBarrel, TrackAssoc const& trackAssocs) { + computeOccupancyEstimators(collisions, tracksPosWithCov, tracksNegWithCov, presliceWithCov, bcs); fullSkimming(collisions, bcs, zdcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr); } diff --git a/PWGDQ/Tasks/filterPbPb.cxx b/PWGDQ/Tasks/filterPbPb.cxx index 1544ea14f21..4841078f61d 100644 --- a/PWGDQ/Tasks/filterPbPb.cxx +++ b/PWGDQ/Tasks/filterPbPb.cxx @@ -9,16 +9,21 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -#include -#include -#include -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "PWGDQ/DataModel/ReducedInfoTables.h" #include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "PWGUD/Core/SGSelector.h" + #include "CommonConstants/LHCConstants.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/Vertex.h" -#include "PWGUD/Core/SGSelector.h" +#include + +#include + +#include +#include +#include using namespace std; @@ -32,7 +37,6 @@ using MyBCs = soa::Join eventRapidityGapFilter; - Produces eventFilter; OutputObj fFilterOutcome{"Filter outcome"}; Configurable fConfigNDtColl{"cfgNDtColl", 4, "Number of standard deviations to consider in BC range"}; @@ -79,189 +83,6 @@ struct DQFilterPbPbTask { static_cast((fConfigUseFDD) ? fConfigFDDCAmpLimit : -1.)}); } - // Helper function for selecting double gap and single gap events - template - uint64_t rapidityGapFilter(TEvent const& collision, TBCs const& bcs, - aod::FT0s& ft0s, aod::FV0As& fv0as, aod::FDDs& fdds, - std::vector FITAmpLimits, int nDtColl, int minNBCs, int minNPVCs, int maxNPVCs, float maxFITTime, - bool useFV0, bool useFT0, bool useFDD) - { - fFilterOutcome->Fill(1., 1.); - - // Number of primary vertex contributors - if (collision.numContrib() < minNPVCs || collision.numContrib() > maxNPVCs) { - fFilterOutcome->Fill(6, 1); - return 0; - } - - // Find BC associated with collision - if (!collision.has_foundBC()) { - fFilterOutcome->Fill(7., 1); - return 0; - } - // foundBCId is stored in EvSels - auto bc = collision.template foundBC_as(); - - // Obtain slice of compatible BCs - uint64_t mostProbableBC = bc.globalBC(); - uint64_t meanBC = mostProbableBC + std::lround(collision.collisionTime() / o2::constants::lhc::LHCBunchSpacingNS); - int deltaBC = std::ceil(collision.collisionTimeRes() / o2::constants::lhc::LHCBunchSpacingNS * nDtColl); - if (deltaBC < minNBCs) { - deltaBC = minNBCs; - } - - // Range of BCs to consider - uint64_t minBC = (uint64_t)deltaBC < meanBC ? meanBC - (uint64_t)deltaBC : 0; - uint64_t maxBC = meanBC + (uint64_t)deltaBC; - - int slicemin = 0; - int slicemax = 0; - // Check if there is overlap between acceptable and possible BCs - if (maxBC > bcs.iteratorAt(0).globalBC() && minBC < bcs.iteratorAt(bcs.size() - 1).globalBC()) { - // find slice of BCs table with BC in [minBC, maxBC] - int moveCount = 0; - int64_t minBCId = bc.globalIndex(); - int64_t maxBCId = bc.globalIndex(); - // lower limit - if (bc.globalBC() < minBC) { - while (bc != bcs.end() && bc.globalBC() < minBC) { - ++bc; - ++moveCount; - minBCId = bc.globalIndex(); - } - } else { - while (bc.globalIndex() > 0 && bc.globalBC() >= minBC) { - minBCId = bc.globalIndex(); - --bc; - --moveCount; - } - } - // upper limit - if (bc.globalBC() < maxBC) { - while (bc != bcs.end() && bc.globalBC() <= maxBC) { - maxBCId = bc.globalIndex(); - ++bc; - ++moveCount; - } - } else { - while (bc.globalIndex() > 0 && bc.globalBC() > maxBC) { - --bc; - --moveCount; - maxBCId = bc.globalIndex(); - } - } - // reset bc - bc.moveByIndex(-moveCount); - // Create BC slice - slicemin = minBCId; - slicemax = maxBCId - minBCId + 1; - } - MyBCs bcrange{{bcs.asArrowTable()->Slice(slicemin, slicemax)}, (uint16_t)slicemin}; - // Rapidity gap condition: Check FIT activity in BC range - bool isSideAClean = true; - bool isSideCClean = true; - - for (auto const& bc : bcrange) { - if (useFV0) { - if (bc.has_foundFV0()) { - auto fv0a = fv0as.iteratorAt(bc.foundFV0Id()); - float FV0Amplitude = 0; - for (auto amp : fv0a.amplitude()) { - FV0Amplitude += amp; - } - float FV0Time = std::abs(fv0a.time()); - if (FV0Amplitude > FITAmpLimits[0] || FV0Time > maxFITTime) { - isSideAClean = false; - } - } - } - if (useFT0) { - if (bc.has_foundFT0()) { - auto ft0 = ft0s.iteratorAt(bc.foundFT0Id()); - float FT0AAmplitude = 0; - float FT0CAmplitude = 0; - for (auto amp : ft0.amplitudeA()) { - FT0AAmplitude += amp; - } - for (auto amp : ft0.amplitudeC()) { - FT0CAmplitude += amp; - } - float FT0ATime = std::abs(ft0.timeA()); - float FT0CTime = std::abs(ft0.timeC()); - if (FT0AAmplitude > FITAmpLimits[1] || FT0ATime > maxFITTime) { - isSideAClean = false; - } - if (FT0CAmplitude > FITAmpLimits[2] || FT0CTime > maxFITTime) { - isSideCClean = false; - } - } - } - if (useFDD) { - if (bc.has_foundFDD()) { - auto fdd = fdds.iteratorAt(bc.foundFDDId()); - float FDDAAmplitude = 0; - float FDDCAmplitude = 0; - for (auto amp : fdd.chargeA()) { - FDDAAmplitude += amp; - } - for (auto amp : fdd.chargeC()) { - FDDCAmplitude += amp; - } - float FDDATime = std::abs(fdd.timeA()); - float FDDCTime = std::abs(fdd.timeC()); - if (FDDAAmplitude > FITAmpLimits[3] || FDDATime > maxFITTime) { - isSideAClean = false; - } - if (FDDCAmplitude > FITAmpLimits[4] || FDDCTime > maxFITTime) { - isSideCClean = false; - } - } - } - } - - // Compute FIT decision - uint64_t FITDecision = 0; - if (isSideAClean && isSideCClean) { - fFilterOutcome->Fill(2, 1); - FITDecision |= (uint64_t(1) << VarManager::kDoubleGap); - } - if (isSideAClean && !isSideCClean) { - fFilterOutcome->Fill(3, 1); - FITDecision |= (uint64_t(1) << VarManager::kSingleGapA); - } else if (!isSideAClean && isSideCClean) { - fFilterOutcome->Fill(4, 1); - FITDecision |= (uint64_t(1) << VarManager::kSingleGapC); - } else if (!isSideAClean && !isSideCClean) { - fFilterOutcome->Fill(5, 1); - } - - if (!FITDecision) { - return 0; - } - - // Return filter bitmap corresponding to FIT decision - return FITDecision; - } - - void processFilterPbPb(MyEvents::iterator const& collision, MyBCs const& bcs, - aod::FT0s& ft0s, aod::FV0As& fv0as, aod::FDDs& fdds) - { - std::vector FITAmpLimits = {fConfigFV0AmpLimit, fConfigFT0AAmpLimit, fConfigFT0CAmpLimit, fConfigFDDAAmpLimit, fConfigFDDCAmpLimit}; - - uint64_t filter = rapidityGapFilter(collision, bcs, ft0s, fv0as, fdds, - FITAmpLimits, fConfigNDtColl, fConfigMinNBCs, fConfigMinNPVCs, fConfigMaxNPVCs, fConfigMaxFITTime, - fConfigUseFV0, fConfigUseFT0, fConfigUseFDD); - - // Record whether UPC settings were used for this event - if (collision.flags() & dataformats::Vertex>::Flags::UPCMode) { - filter |= (uint64_t(1) << VarManager::kITSUPCMode); - fFilterOutcome->Fill(8, 1); - } - - eventRapidityGapFilter(filter, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0); - eventFilter(filter); - } - void processUDSGSelector(MyEvents::iterator const& collision, MyBCs const& bcs, aod::FT0s& ft0s, aod::FV0As& fv0as, aod::FDDs& fdds, aod::Zdcs& /*zdcs*/) { @@ -279,13 +100,13 @@ struct DQFilterPbPbTask { int issgevent = isSGEvent.value; // Translate SGSelector values to DQEventFilter values if (issgevent == 0) { - filter |= (uint64_t(1) << VarManager::kSingleGapA); + filter |= (static_cast(1) << VarManager::kSingleGapA); fFilterOutcome->Fill(3, 1); } else if (issgevent == 1) { - filter |= (uint64_t(1) << VarManager::kSingleGapC); + filter |= (static_cast(1) << VarManager::kSingleGapC); fFilterOutcome->Fill(4, 1); } else if (issgevent == 2) { - filter |= (uint64_t(1) << VarManager::kDoubleGap); + filter |= (static_cast(1) << VarManager::kDoubleGap); fFilterOutcome->Fill(2, 1); } else if (issgevent == 3) { fFilterOutcome->Fill(5, 1); @@ -302,20 +123,11 @@ struct DQFilterPbPbTask { // Record whether UPC settings were used for this event if (collision.flags() & dataformats::Vertex>::Flags::UPCMode) { - filter |= (uint64_t(1) << VarManager::kITSUPCMode); + filter |= (static_cast(1) << VarManager::kITSUPCMode); fFilterOutcome->Fill(8, 1); } - if (newbc.has_zdc()) { - auto zdc = newbc.zdc(); - eventRapidityGapFilter(filter, fitInfo.ampFT0A, fitInfo.ampFT0C, fitInfo.ampFDDA, fitInfo.ampFDDC, fitInfo.ampFV0A, - zdc.energyCommonZNA(), zdc.energyCommonZNC(), zdc.energyCommonZPA(), zdc.energyCommonZPC(), - zdc.timeZNA(), zdc.timeZNC(), zdc.timeZPA(), zdc.timeZPC()); - } else { - eventRapidityGapFilter(filter, fitInfo.ampFT0A, fitInfo.ampFT0C, fitInfo.ampFDDA, fitInfo.ampFDDC, fitInfo.ampFV0A, - -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0); - } - eventFilter(filter); + eventRapidityGapFilter(filter, newbc.globalIndex()); } void processDummy(MyEvents&) @@ -323,7 +135,6 @@ struct DQFilterPbPbTask { // do nothing } - PROCESS_SWITCH(DQFilterPbPbTask, processFilterPbPb, "Run filter task", true); PROCESS_SWITCH(DQFilterPbPbTask, processUDSGSelector, "Run filter task with SG selector from UD", false); PROCESS_SWITCH(DQFilterPbPbTask, processDummy, "Dummy function", false); };