diff --git a/PWGDQ/Tasks/DalitzSelection.cxx b/PWGDQ/Tasks/DalitzSelection.cxx index 70eb8415d27..317015a044f 100644 --- a/PWGDQ/Tasks/DalitzSelection.cxx +++ b/PWGDQ/Tasks/DalitzSelection.cxx @@ -24,10 +24,10 @@ #include "PWGDQ/Core/VarManager.h" #include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "Common/DataModel/CollisionAssociationTables.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" -#include "Common/DataModel/TrackSelectionTables.h" #include #include @@ -52,6 +52,7 @@ #include #include #include +#include #include using namespace o2; @@ -60,20 +61,23 @@ using namespace o2::framework::expressions; using namespace o2::aod; using namespace o2::soa; -using MyEvents = soa::Join; +// using MyEvents = soa::Join +using MyEventsWithCent = soa::Join; -using MyBarrelTracks = soa::Join; -constexpr static uint32_t EventFillMap = VarManager::ObjTypes::Collision; -constexpr static uint32_t TrackFillMap = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackSelection | VarManager::ObjTypes::TrackPID; +// constexpr static uint32_t EventFillMap = VarManager::ObjTypes::Collision; +constexpr static uint32_t EventFillMapWithCent = VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCent; +constexpr static uint32_t TrackFillMap = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackPID; struct DalitzSelection { Produces dalitzbits; Preslice perCollision = aod::track::collisionId; + Preslice trackIndicesPerCollision = aod::track_association::collisionId; // Configurables // cuts @@ -102,6 +106,7 @@ struct DalitzSelection { struct : ConfigurableGroup { Configurable fConfigEnableLikeSign{"cfgEnableLikeSign", false, "Whether or not also add like-sign pairs (for studying combinatorial background which might contain misID or non-primary electrons)"}; Configurable fQA{"cfgQA", true, "QA histograms"}; + Configurable fRemoveDoubleCounting{"cfgRemoveDoubleCounting", true, "If a track/pair is selected for several collisions, we still count it only once"}; // Configurables for TPC post-calibration maps Configurable fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable fConfigCcdbPathTPC{"ccdb-path-tpc", "Users/i/iarsene/Calib/TPCpostCalib", "base path to the ccdb object"}; @@ -126,6 +131,11 @@ struct DalitzSelection { std::map fDalitzmap; // whether it is selected as dalitz decay daughter with symmetric or tag cut std::map fDalitzmapProbe; // whether it is selected as dalitz decay daughter with probe cut + // maps to remove ambiguities + std::map fDalitzmapAmbiguity; + std::map fDalitzmapProbeAmbiguity; + std::map, uint8_t> fAmbiguousPairs; + AnalysisCompositeCut* fEventCut; std::vector fTrackCuts; std::vector fTrackCutsProbe; @@ -295,13 +305,22 @@ struct DalitzSelection { fOutputList.setObject(fHistMan->GetMainHistogramList()); } - template - void runTrackSelection(TTracks const& tracksBarrel) + template + void runTrackSelection(TTracks const& tracksBarrel, TEvent const& collision) { - for (const auto& track : tracksBarrel) { + for (const auto& track1 : tracksBarrel) { uint8_t filterMap = uint8_t(0); uint8_t filterMapProbe = uint8_t(0); - VarManager::FillTrack(track); + int fullTrackIdx = 0; + if constexpr (isReassoc) { + auto const& fullTrack = track1.template track_as(); + VarManager::FillTrack(fullTrack); + VarManager::FillTrackCollision(fullTrack, collision); + fullTrackIdx = fullTrack.globalIndex(); + } else { + VarManager::FillTrack(track1); + fullTrackIdx = track1.globalIndex(); + } int i = 0; for (auto cut = fTrackCuts.begin(); cut != fTrackCuts.end(); ++cut, ++i) { if ((*cut).IsSelected(VarManager::fgValues)) { @@ -317,35 +336,63 @@ struct DalitzSelection { } } if (filterMap) { - fTrackmap[track.globalIndex()] = filterMap; + fTrackmap[fullTrackIdx] = filterMap; } if (filterMapProbe) { - fTrackmapProbe[track.globalIndex()] = filterMapProbe; + fTrackmapProbe[fullTrackIdx] = filterMapProbe; } } // end loop over tracks } - template - void runDalitzPairing(TTracks const& tracks1, TTracks const& tracks2) + template + void runDalitzPairing(TTracks const& tracks1, TTracks const& tracks2, TEvent collision) { + if (isReassoc & fConfigOptions.fRemoveDoubleCounting) { + // keep track of selections from previous events to check if the track was already selected or not + fDalitzmapAmbiguity = fDalitzmap; + fDalitzmapProbeAmbiguity = fDalitzmapProbe; + } + for (const auto& [track1, track2] : o2::soa::combinations(CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (track1.globalIndex() == track2.globalIndex()) { + + int trackIdx1; + int trackIdx2; + bool isLikeSign; + if constexpr (isReassoc) { + auto const& fullTrack1 = track1.template track_as(); + auto const& fullTrack2 = track2.template track_as(); + trackIdx1 = fullTrack1.globalIndex(); + trackIdx2 = fullTrack2.globalIndex(); + isLikeSign = (fullTrack1.sign() * fullTrack2.sign() > 0); + } else { + trackIdx1 = track1.globalIndex(); + trackIdx2 = track2.globalIndex(); + isLikeSign = (track1.sign() * track2.sign() > 0); + } + + if (trackIdx1 == trackIdx2) { continue; } - if (!fIsTagAndProbe && track1.globalIndex() >= track2.globalIndex()) { + if (!fIsTagAndProbe && trackIdx1 >= trackIdx2) { continue; } - if (!fConfigOptions.fConfigEnableLikeSign && (track1.sign() * track2.sign() > 0)) { + if (!fConfigOptions.fConfigEnableLikeSign && isLikeSign) { continue; } - uint8_t twoTracksFilterMap = fTrackmap[track1.globalIndex()] & (fIsTagAndProbe ? fTrackmapProbe[track2.globalIndex()] : fTrackmap[track2.globalIndex()]); + uint8_t twoTracksFilterMap = fTrackmap[trackIdx1] & (fIsTagAndProbe ? fTrackmapProbe[trackIdx2] : fTrackmap[trackIdx2]); if (!twoTracksFilterMap) { continue; } // pairing - VarManager::FillPair(track1, track2); + if constexpr (isReassoc) { + auto const& fullTrack1 = track1.template track_as(); + auto const& fullTrack2 = track2.template track_as(); + VarManager::FillPair(fullTrack1, fullTrack2); + } else { + VarManager::FillPair(track1, track2); + } // Fill pair selection map and fill pair histogram int icut = 0; @@ -355,21 +402,35 @@ struct DalitzSelection { continue; } if ((*pairCut).IsSelected(VarManager::fgValues)) { - if (track1.sign() * track2.sign() < 0) { - fDalitzmap[track1.globalIndex()] |= (uint8_t(1) << icut); + bool isPairAlreadySelected = false; + if (isReassoc & fConfigOptions.fRemoveDoubleCounting) { + // if we remove double counting and the pair is already selected, we don't fill the histograms + std::pair iPair(trackIdx1, trackIdx2); + if (fAmbiguousPairs.find(iPair) != fAmbiguousPairs.end()) { + if (fAmbiguousPairs[iPair] & (static_cast(1) << icut)) { // if this pair is already stored with this cut + isPairAlreadySelected = true; + } else { + fAmbiguousPairs[iPair] |= static_cast(1) << icut; + } + } else { + fAmbiguousPairs[iPair] = static_cast(1) << icut; + } + } + if (!isLikeSign) { + fDalitzmap[trackIdx1] |= (uint8_t(1) << icut); if (fIsTagAndProbe) { - fDalitzmapProbe[track2.globalIndex()] |= (uint8_t(1) << icut); - if (fConfigOptions.fQA) { + fDalitzmapProbe[trackIdx2] |= (uint8_t(1) << icut); + if (fConfigOptions.fQA && !isPairAlreadySelected) { fHistMan->FillHistClass(Form("Pair_%s_%s_%s", (*trackCut).GetName(), fTrackCutsProbe.at(icut).GetName(), (*pairCut).GetName()), VarManager::fgValues); } } else { - fDalitzmap[track2.globalIndex()] |= (uint8_t(1) << icut); - if (fConfigOptions.fQA) { + fDalitzmap[trackIdx2] |= (uint8_t(1) << icut); + if (fConfigOptions.fQA && !isPairAlreadySelected) { fHistMan->FillHistClass(Form("Pair_%s_%s", (*trackCut).GetName(), (*pairCut).GetName()), VarManager::fgValues); } } } else { - if (fConfigOptions.fQA) { + if (fConfigOptions.fQA && !isPairAlreadySelected) { fHistMan->FillHistClass(fIsTagAndProbe ? Form("PairLS_%s_%s_%s", (*trackCut).GetName(), fTrackCutsProbe.at(icut).GetName(), (*pairCut).GetName()) : Form("PairLS_%s_%s", (*trackCut).GetName(), (*pairCut).GetName()), VarManager::fgValues); } } // end if like-sign @@ -380,12 +441,37 @@ struct DalitzSelection { // Fill Hists if (fConfigOptions.fQA) { for (const auto& track : tracks1) { - uint8_t filterMap = fDalitzmap[track.globalIndex()]; - uint8_t filterMapProbe = fDalitzmapProbe[track.globalIndex()]; - if (!filterMap && !filterMapProbe) { - continue; + uint8_t filterMap; + uint8_t filterMapProbe; + if constexpr (isReassoc) { + auto const& fullTrack = track.template track_as(); + filterMap = fDalitzmap[fullTrack.globalIndex()]; + filterMapProbe = fDalitzmapProbe[fullTrack.globalIndex()]; + if (fConfigOptions.fRemoveDoubleCounting) { + // we remove track selections which were already selected before + uint8_t previousFilterMap = fDalitzmapAmbiguity[fullTrack.globalIndex()]; + uint8_t previousFilterMapProbe = fDalitzmapProbeAmbiguity[fullTrack.globalIndex()]; + filterMap &= ~previousFilterMap; + filterMapProbe &= ~previousFilterMapProbe; + } + + if (!filterMap && !filterMapProbe) { + continue; + } + + VarManager::FillTrack(fullTrack); + // The caveat here is that we only fill the DCA to the first selected collision + VarManager::FillTrackCollision(fullTrack, collision); + } else { + filterMap = fDalitzmap[track.globalIndex()]; + filterMapProbe = fDalitzmapProbe[track.globalIndex()]; + + if (!filterMap && !filterMapProbe) { + continue; + } + + VarManager::FillTrack(track); } - VarManager::FillTrack(track); int icut = 0; auto trackCut = fTrackCuts.begin(); @@ -403,7 +489,41 @@ struct DalitzSelection { } } - void processFullTracks(MyEvents const& collisions, aod::BCsWithTimestamps const&, soa::Filtered const& filteredTracks, MyBarrelTracks const& tracks) + void initNewRun(int64_t timestamp) + { + + VarManager::ResetValues(0, VarManager::kNRunWiseVariables); + + // We setup the magnetic field, because the conversion rejection cut might depend on it + float magField = 0.; + if (fConfigOptions.fUseRemoteField.value) { + grpmag = fCCDB->getForTimeStamp(fConfigOptions.grpmagPath, timestamp); + if (grpmag != nullptr) { + magField = grpmag->getNominalL3Field(); + } else { + LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", timestamp); + } + } else { + magField = fConfigOptions.fConfigMagField.value; + } + LOGF(info, "setting mag field to %f", magField); + if (magField == 0.) { + LOGF(fatal, "magnetic field not set correctly, please check"); + } + VarManager::SetMagneticField(magField); + + if (fConfigOptions.fConfigComputeTPCpostCalib) { + auto calibList = fCCDB->getForTimeStamp(fConfigOptions.fConfigCcdbPathTPC.value, timestamp); + VarManager::SetCalibrationObject(VarManager::kTPCElectronMean, calibList->FindObject("mean_map_electron")); + VarManager::SetCalibrationObject(VarManager::kTPCElectronSigma, calibList->FindObject("sigma_map_electron")); + VarManager::SetCalibrationObject(VarManager::kTPCPionMean, calibList->FindObject("mean_map_pion")); + VarManager::SetCalibrationObject(VarManager::kTPCPionSigma, calibList->FindObject("sigma_map_pion")); + VarManager::SetCalibrationObject(VarManager::kTPCProtonMean, calibList->FindObject("mean_map_proton")); + VarManager::SetCalibrationObject(VarManager::kTPCProtonSigma, calibList->FindObject("sigma_map_proton")); + } + } + + void processFullTracks(MyEventsWithCent const& collisions, aod::BCsWithTimestamps const&, soa::Filtered const& filteredTracks, MyBarrelTracks const& tracks) { const int pairType = VarManager::kDecayToEE; fDalitzmap.clear(); @@ -413,7 +533,7 @@ struct DalitzSelection { fTrackmap.clear(); fTrackmapProbe.clear(); VarManager::ResetValues(VarManager::kNRunWiseVariables, VarManager::kNBarrelTrackVariables); - VarManager::FillEvent(collision); + VarManager::FillEvent(collision); bool isEventSelected = fEventCut->IsSelected(VarManager::fgValues); if (isEventSelected) { @@ -423,41 +543,52 @@ struct DalitzSelection { auto bc = collision.template bc_as(); if (fCurrentRun != bc.runNumber()) { - VarManager::ResetValues(0, VarManager::kNRunWiseVariables); - - // We setup the magnetic field, because the conversion rejection cut might depend on it - float magField = 0.; - if (fConfigOptions.fUseRemoteField.value) { - grpmag = fCCDB->getForTimeStamp(fConfigOptions.grpmagPath, bc.timestamp()); - if (grpmag != nullptr) { - magField = grpmag->getNominalL3Field(); - } else { - LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", bc.timestamp()); - } - } else { - magField = fConfigOptions.fConfigMagField.value; - } - LOGF(info, "setting mag field to %f", magField); - if (magField == 0.) { - LOGF(fatal, "magnetic field not set correctly, please check"); - } - VarManager::SetMagneticField(magField); - - if (fConfigOptions.fConfigComputeTPCpostCalib) { - auto calibList = fCCDB->getForTimeStamp(fConfigOptions.fConfigCcdbPathTPC.value, bc.timestamp()); - VarManager::SetCalibrationObject(VarManager::kTPCElectronMean, calibList->FindObject("mean_map_electron")); - VarManager::SetCalibrationObject(VarManager::kTPCElectronSigma, calibList->FindObject("sigma_map_electron")); - VarManager::SetCalibrationObject(VarManager::kTPCPionMean, calibList->FindObject("mean_map_pion")); - VarManager::SetCalibrationObject(VarManager::kTPCPionSigma, calibList->FindObject("sigma_map_pion")); - VarManager::SetCalibrationObject(VarManager::kTPCProtonMean, calibList->FindObject("mean_map_proton")); - VarManager::SetCalibrationObject(VarManager::kTPCProtonSigma, calibList->FindObject("sigma_map_proton")); - } + initNewRun(bc.timestamp()); fCurrentRun = bc.runNumber(); } auto groupedFilteredTracks = filteredTracks.sliceBy(perCollision, collision.globalIndex()); - runTrackSelection(groupedFilteredTracks); - runDalitzPairing(groupedFilteredTracks, groupedFilteredTracks); + runTrackSelection(groupedFilteredTracks, nullptr); + runDalitzPairing(groupedFilteredTracks, groupedFilteredTracks, nullptr); + } + } + + for (const auto& track : tracks) { // Fill dalitz bits + dalitzbits(fIsTagAndProbe ? fDalitzmapProbe[track.globalIndex()] : fDalitzmap[track.globalIndex()]); + } + } + + void processFullTracksWithAssoc(MyEventsWithCent const& collisions, aod::BCsWithTimestamps const&, MyBarrelTracks const& tracks, TrackAssoc const& trackAssocs) + { + + const int pairType = VarManager::kDecayToEE; + fDalitzmap.clear(); + fDalitzmapProbe.clear(); + fDalitzmapAmbiguity.clear(); + fDalitzmapProbeAmbiguity.clear(); + fAmbiguousPairs.clear(); + + for (const auto& collision : collisions) { + fTrackmap.clear(); + fTrackmapProbe.clear(); + VarManager::ResetValues(VarManager::kNRunWiseVariables, VarManager::kNBarrelTrackVariables); + VarManager::FillEvent(collision); + bool isEventSelected = fEventCut->IsSelected(VarManager::fgValues); + + if (isEventSelected) { + + reinterpret_cast(fStatsList->At(0))->Fill(0); + + auto bc = collision.template bc_as(); + + if (fCurrentRun != bc.runNumber()) { + initNewRun(bc.timestamp()); + fCurrentRun = bc.runNumber(); + } + + auto groupedTracksAssoc = trackAssocs.sliceBy(trackIndicesPerCollision, collision.globalIndex()); + runTrackSelection(groupedTracksAssoc, collision); + runDalitzPairing(groupedTracksAssoc, groupedTracksAssoc, collision); } } @@ -466,10 +597,11 @@ struct DalitzSelection { } } - void processDummy(MyEvents const&) + void processDummy(MyEventsWithCent const&) { } + PROCESS_SWITCH(DalitzSelection, processFullTracksWithAssoc, "Run Dalitz selection on AO2D tables with reassociation", false); PROCESS_SWITCH(DalitzSelection, processFullTracks, "Run Dalitz selection on AO2D tables", false); PROCESS_SWITCH(DalitzSelection, processDummy, "Do nothing", false); };