diff --git a/Common/TableProducer/eventSelection.cxx b/Common/TableProducer/eventSelection.cxx index 9aa2a664544..be359cc2f61 100644 --- a/Common/TableProducer/eventSelection.cxx +++ b/Common/TableProducer/eventSelection.cxx @@ -36,8 +36,9 @@ MetadataHelper metadataInfo; // Metadata helper using BCsWithRun2InfosTimestampsAndMatches = soa::Join; using BCsWithRun3Matchings = soa::Join; using BCsWithBcSelsRun2 = soa::Join; -using BCsWithBcSelsRun3 = soa::Join; +using BCsWithBcSelsRun3 = soa::Join; using FullTracksIU = soa::Join; +const double bcNS = o2::constants::lhc::LHCBunchSpacingNS; struct BcSelectionTask { Produces bcsel; @@ -454,6 +455,7 @@ struct EventSelectionTask { Configurable confTimeRangeVetoOnCollStandard{"TimeRangeVetoOnCollStandard", 10, "Exclusion of a collision if there are other collisions nearby, +/- us"}; Configurable confTimeRangeVetoOnCollNarrow{"TimeRangeVetoOnCollNarrow", 4, "Exclusion of a collision if there are other collisions nearby, +/- us"}; Configurable confUseWeightsForOccupancyVariable{"UseWeightsForOccupancyEstimator", 1, "Use or not the delta-time weights for the occupancy estimator"}; + Configurable confSigmaBCforHighPtTracks{"confSigmaBCforHighPtTracks", 4, "Custom sigma (in bcs) for collisions with high-pt tracks"}; Partition tracklets = (aod::track::trackType == static_cast(o2::aod::track::TrackTypeEnum::Run2Tracklet)); @@ -480,6 +482,40 @@ struct EventSelectionTask { return (dbc1 <= dbc2) ? index1 : index2; } + // helper function to find median time in the vector of TOF or TRD-track times + float getMedian(std::vector v) + { + int medianIndex = v.size() / 2; + std::nth_element(v.begin(), v.begin() + medianIndex, v.end()); + return v[medianIndex]; + } + + // helper function to find closest TVX signal in time and in zVtx + int64_t findBestGlobalBC(int64_t meanBC, int64_t sigmaBC, int32_t nContrib, float zVtxCol, std::map& mapGlobalBcVtxZ) + { + int64_t minBC = meanBC - 3 * sigmaBC; + int64_t maxBC = meanBC + 3 * sigmaBC; + // TODO: use ITS ROF bounds to reduce the search range? + + float zVtxSigma = 2.7 * pow(nContrib, -0.466) + 0.024; + zVtxSigma += 1.0; // additional uncertainty due to imperfectections of FT0 time calibration + + auto itMin = mapGlobalBcVtxZ.lower_bound(minBC); + auto itMax = mapGlobalBcVtxZ.upper_bound(maxBC); + + float bestChi2 = 1e+10; + int64_t bestGlobalBC = 0; + for (std::map::iterator it = itMin; it != itMax; ++it) { + float chi2 = pow((it->second - zVtxCol) / zVtxSigma, 2) + pow((it->first - meanBC) / sigmaBC, 2.); + if (chi2 < bestChi2) { + bestChi2 = chi2; + bestGlobalBC = it->first; + } + } + + return bestGlobalBC; + } + void init(InitContext&) { if (metadataInfo.isFullyDefined()) { // Check if the metadata is initialized (only if not forced from the workflow configuration) @@ -501,7 +537,6 @@ struct EventSelectionTask { } } - // ccdb->setURL("http://ccdb-test.cern.ch:8080"); ccdb->setURL("http://alice-ccdb.cern.ch"); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -517,7 +552,7 @@ struct EventSelectionTask { evsel.reserve(collisions.size()); } - void processRun2(aod::Collision const& col, BCsWithBcSelsRun2 const&, aod::Tracks const&, aod::FV0Cs const&) + void processRun2(aod::Collision const& col, BCsWithBcSelsRun2 const&, FullTracksIU const&, aod::FV0Cs const&) { auto bc = col.bc_as(); EventSelectionParams* par = ccdb->getForTimeStamp("EventSelection/EventSelectionParams", bc.timestamp()); @@ -589,8 +624,8 @@ struct EventSelectionTask { } PROCESS_SWITCH(EventSelectionTask, processRun2, "Process Run2 event selection", true); - Preslice perCollision = aod::track::collisionId; - void processRun3(aod::Collisions const& cols, FullTracksIU const& tracks, BCsWithBcSelsRun3 const& bcs, aod::FT0s const&) + Partition pvTracks = ((aod::track::flags & (uint32_t)o2::aod::track::PVContributor) == (uint32_t)o2::aod::track::PVContributor); + void processRun3(aod::Collisions const& cols, FullTracksIU const&, BCsWithBcSelsRun3 const& bcs, aod::FT0s const&) { int run = bcs.iteratorAt(0).runNumber(); // extract bc pattern from CCDB for data or anchored MC only @@ -600,7 +635,6 @@ struct EventSelectionTask { auto grplhcif = ccdb->getForTimeStamp("GLO/Config/GRPLHCIF", ts); bcPatternB = grplhcif->getBunchFilling().getBCPattern(); - // EventSelectionParams* par = ccdb->getForTimeStamp("EventSelection/EventSelectionParams", ts); // access orbit-reset timestamp auto ctpx = ccdb->getForTimeStamp>("CTP/Calib/OrbitReset", ts); @@ -621,26 +655,24 @@ struct EventSelectionTask { nBCsPerTF = nOrbitsPerTF * o2::constants::lhc::LHCMaxBunches; } - // create maps from globalBC to bc index for TVX or FT0-OR fired bcs - // to be used for closest TVX (FT0-OR) searches + // create maps from globalBC to bc index for TVX-fired bcs + // to be used for closest TVX searches std::map mapGlobalBcWithTVX; - std::map mapGlobalBcWithTOR; + std::map mapGlobalBcVtxZ; for (auto& bc : bcs) { int64_t globalBC = bc.globalBC(); // skip non-colliding bcs for data and anchored runs if (run >= 500000 && bcPatternB[globalBC % o2::constants::lhc::LHCMaxBunches] == 0) { continue; } - if (bc.selection_bit(kIsBBT0A) || bc.selection_bit(kIsBBT0C)) { - mapGlobalBcWithTOR[globalBC] = bc.globalIndex(); - } if (bc.selection_bit(kIsTriggerTVX)) { mapGlobalBcWithTVX[globalBC] = bc.globalIndex(); + mapGlobalBcVtxZ[globalBC] = bc.has_ft0() ? bc.ft0().posZ() : 0; } } // protection against empty FT0 maps - if (mapGlobalBcWithTOR.size() == 0 || mapGlobalBcWithTVX.size() == 0) { + if (mapGlobalBcWithTVX.size() == 0) { LOGP(error, "FT0 table is empty or corrupted. Filling evsel table with dummy values"); for (auto& col : cols) { auto bc = col.bc_as(); @@ -653,88 +685,139 @@ struct EventSelectionTask { } return; } - - std::vector vFoundBCindex(cols.size(), -1); // indices of found bcs - std::vector vIsVertexITSTPC(cols.size(), 0); // at least one of vertex contributors is ITS-TPC track - std::vector vIsVertexTOFmatched(cols.size(), 0); // at least one of vertex contributors is matched to TOF - std::vector vIsVertexTRDmatched(cols.size(), 0); // at least one of vertex contributors is matched to TRD - std::vector vCollisionsPerBc(bcs.size(), 0); // counter of collisions per found bc for pileup checks - - // for the occupancy study - std::vector vFoundGlobalBC(cols.size(), 0); // global BCs for collisions std::vector vTracksITS567perColl(cols.size(), 0); // counter of tracks per found bc for occupancy studies std::vector vIsFullInfoForOccupancy(cols.size(), 0); // info for occupancy in +/- windows is available (i.e. a given coll is not too close to the TF borders) const float timeWinOccupancyCalcMinNS = confTimeIntervalForOccupancyCalculationMin * 1e3; // ns const float timeWinOccupancyCalcMaxNS = confTimeIntervalForOccupancyCalculationMax * 1e3; // ns - // const double timeWinOccupancyExclusionRangeNS = confExclusionIntervalForOccupancyCalculation * 1e3; // ns - const double bcNS = o2::constants::lhc::LHCBunchSpacingNS; - - // loop to find nearest bc with FT0 entry -> foundBC index + std::vector vIsVertexITSTPC(cols.size(), 0); // at least one of vertex contributors is ITS-TPC track + std::vector vIsVertexTOFmatched(cols.size(), 0); // at least one of vertex contributors is matched to TOF + std::vector vIsVertexTRDmatched(cols.size(), 0); // at least one of vertex contributors is matched to TRD + + std::vector vCollisionsPerBc(bcs.size(), 0); // counter of collisions per found bc for pileup checks + std::vector vFoundBCindex(cols.size(), -1); // indices of found bcs + std::vector vFoundGlobalBC(cols.size(), 0); // global BCs for collisions + + std::vector vIsVertexTOF(cols.size(), 0); + std::vector vIsVertexTRD(cols.size(), 0); + std::vector vIsVertexTPC(cols.size(), 0); + std::vector vIsVertexHighPtTPC(cols.size(), 0); + std::vector vNcontributors(cols.size(), 0); + std::vector vWeightedTimesTPCnoTOFnoTRD(cols.size(), 0); + std::vector vWeightedSigmaTPCnoTOFnoTRD(cols.size(), 0); + + // temporary vectors to find tracks with median time + std::vector vTrackTimesTOF; + std::vector vTrackTimesTRDnoTOF; + + // first loop to match collisions to TVX for (auto& col : cols) { + int32_t colIndex = col.globalIndex(); auto bc = col.bc_as(); - int64_t meanBC = bc.globalBC(); - const double bcNS = o2::constants::lhc::LHCBunchSpacingNS; - int64_t deltaBC = std::ceil(col.collisionTimeRes() / bcNS * 4); - - // count tracks of different types - int nITS567cls = 0; - int nITSTPCtracks = 0; - int nTOFtracks = 0; - int nTRDtracks = 0; - double timeFromTOFtracks = 0; - auto tracksGrouped = tracks.sliceBy(perCollision, col.globalIndex()); - for (auto& track : tracksGrouped) { - if (!track.isPVContributor()) { - continue; - } - nITSTPCtracks += track.hasITS() && track.hasTPC(); - nTOFtracks += track.hasTOF(); - nTRDtracks += track.hasTRD(); - // calculate average time using TOF tracks - if (track.hasTOF()) { - timeFromTOFtracks += track.trackTime(); - } + int64_t globalBC = bc.globalBC(); + int64_t bcInTF = (bc.globalBC() - bcSOR) % nBCsPerTF; + vIsFullInfoForOccupancy[colIndex] = ((bcInTF - 300) * bcNS > -timeWinOccupancyCalcMinNS) && ((nBCsPerTF - 4000 - bcInTF) * bcNS > timeWinOccupancyCalcMaxNS) ? true : false; + const auto& colPvTracks = pvTracks.sliceByCached(aod::track::collisionId, col.globalIndex(), cache); + vTrackTimesTOF.clear(); + vTrackTimesTRDnoTOF.clear(); + int nPvTracksTPCnoTOFnoTRD = 0; + int nPvTracksHighPtTPCnoTOFnoTRD = 0; + float sumTime = 0, sumW = 0, sumHighPtTime = 0, sumHighPtW = 0; + for (auto& track : colPvTracks) { + float trackTime = track.trackTime(); if (track.itsNCls() >= 5) - nITS567cls++; + vTracksITS567perColl[colIndex]++; + if (track.hasTRD()) + vIsVertexTRDmatched[colIndex] = 1; + if (track.hasTPC()) + vIsVertexITSTPC[colIndex] = 1; + if (track.hasTOF()) { + vTrackTimesTOF.push_back(trackTime); + vIsVertexTOFmatched[colIndex] = 1; + } else if (track.hasTRD()) { + vTrackTimesTRDnoTOF.push_back(trackTime); + } else if (track.hasTPC()) { + float trackTimeRes = track.trackTimeRes(); + float trackPt = track.pt(); + float w = 1. / (trackTimeRes * trackTimeRes); + sumTime += trackTime * w; + sumW += w; + nPvTracksTPCnoTOFnoTRD++; + if (trackPt > 1) { + sumHighPtTime += trackTime * w; + sumHighPtW += w; + nPvTracksHighPtTPCnoTOFnoTRD++; + } + } } - LOGP(debug, "nContrib={} nITSTPCtracks={} nTOFtracks={} nTRDtracks={}", col.numContrib(), nITSTPCtracks, nTOFtracks, nTRDtracks); - - if (nTOFtracks > 0) { - meanBC += TMath::FloorNint(timeFromTOFtracks / nTOFtracks / bcNS); // assign collision bc using TOF-matched tracks - deltaBC = 4; // use precise bc from TOF tracks with +/-4 bc margin - } else if (nITSTPCtracks > 0) { - deltaBC += 30; // extend deltaBC for collisions built with ITS-TPC tracks only + vWeightedTimesTPCnoTOFnoTRD[colIndex] = sumW > 0 ? sumTime / sumW : 0; + vWeightedSigmaTPCnoTOFnoTRD[colIndex] = sumW > 0 ? sqrt(1. / sumW) : 0; + vNcontributors[colIndex] = colPvTracks.size(); + int nPvTracksTOF = vTrackTimesTOF.size(); + int nPvTracksTRDnoTOF = vTrackTimesTRDnoTOF.size(); + // collision type + vIsVertexTOF[colIndex] = nPvTracksTOF > 0; + vIsVertexTRD[colIndex] = nPvTracksTRDnoTOF > 0; + vIsVertexTPC[colIndex] = nPvTracksTPCnoTOFnoTRD > 0; + vIsVertexHighPtTPC[colIndex] = nPvTracksHighPtTPCnoTOFnoTRD > 0; + + int64_t foundGlobalBC = 0; + int32_t foundBCindex = -1; + + if (nPvTracksTOF > 0) { + // for collisions with TOF tracks: + // take bc corresponding to TOF track with median time + int64_t tofGlobalBC = globalBC + TMath::Nint(getMedian(vTrackTimesTOF) / bcNS); + std::map::iterator it = mapGlobalBcWithTVX.find(tofGlobalBC); + if (it != mapGlobalBcWithTVX.end()) { + foundGlobalBC = it->first; + foundBCindex = it->second; + } + } else if (nPvTracksTPCnoTOFnoTRD == 0 && nPvTracksTRDnoTOF > 0) { + // for collisions with TRD tracks but without TOF or ITSTPC-only tracks: + // take bc corresponding to TRD track with median time + int64_t trdGlobalBC = globalBC + TMath::Nint(getMedian(vTrackTimesTRDnoTOF) / bcNS); + std::map::iterator it = mapGlobalBcWithTVX.find(trdGlobalBC); + if (it != mapGlobalBcWithTVX.end()) { + foundGlobalBC = it->first; + foundBCindex = it->second; + } + } else if (nPvTracksHighPtTPCnoTOFnoTRD > 0) { + // for collisions with high-pt ITSTPC-nonTOF-nonTRD tracks + // search in 3*confSigmaBCforHighPtTracks range (3*4 bcs by default) + int64_t meanBC = globalBC + TMath::Nint(sumHighPtTime / sumHighPtW / bcNS); + int64_t bestGlobalBC = findBestGlobalBC(meanBC, confSigmaBCforHighPtTracks, vNcontributors[colIndex], col.posZ(), mapGlobalBcVtxZ); + if (bestGlobalBC > 0) { + foundGlobalBC = bestGlobalBC; + foundBCindex = mapGlobalBcWithTVX[bestGlobalBC]; + } } - int64_t minBC = meanBC - deltaBC; - int64_t maxBC = meanBC + deltaBC; + // fill foundBC indices and global BCs + // keep current bc if TVX matching failed at this step + vFoundBCindex[colIndex] = foundBCindex >= 0 ? foundBCindex : bc.globalIndex(); + vFoundGlobalBC[colIndex] = foundGlobalBC > 0 ? foundGlobalBC : globalBC; - int32_t indexClosestTVX = findClosest(meanBC, mapGlobalBcWithTVX); - int64_t tvxBC = bcs.iteratorAt(indexClosestTVX).globalBC(); - if (tvxBC >= minBC && tvxBC <= maxBC) { // closest TVX within search region - bc.setCursor(indexClosestTVX); - } else { // no TVX within search region, searching for TOR = T0A | T0C - int32_t indexClosestTOR = findClosest(meanBC, mapGlobalBcWithTOR); - int64_t torBC = bcs.iteratorAt(indexClosestTOR).globalBC(); - if (torBC >= minBC && torBC <= maxBC) { - bc.setCursor(indexClosestTOR); - } - } - int32_t foundBC = bc.globalIndex(); + // erase found global BC with TVX from the pool of bcs for the next loop over low-pt TPCnoTOFnoTRD collisions + if (foundBCindex >= 0) + mapGlobalBcVtxZ.erase(foundGlobalBC); + } + + // second loop to match remaining low-pt TPCnoTOFnoTRD collisions + for (auto& col : cols) { int32_t colIndex = col.globalIndex(); - LOGP(debug, "foundBC = {} globalBC = {}", foundBC, bc.globalBC()); - vFoundBCindex[colIndex] = foundBC; - vIsVertexITSTPC[colIndex] = nITSTPCtracks > 0; - vIsVertexTOFmatched[colIndex] = nTOFtracks > 0; - vIsVertexTRDmatched[colIndex] = nTRDtracks > 0; - vCollisionsPerBc[foundBC]++; - vTracksITS567perColl[colIndex] = nITS567cls; - vFoundGlobalBC[colIndex] = bc.globalBC(); - - // check that this collision has full information inside the time window (taking into account TF borders) - int64_t bcInTF = (bc.globalBC() - bcSOR) % nBCsPerTF; - vIsFullInfoForOccupancy[colIndex] = ((bcInTF - 300) * bcNS > -timeWinOccupancyCalcMinNS) && ((nBCsPerTF - 4000 - bcInTF) * bcNS > timeWinOccupancyCalcMaxNS) ? true : false; + if (vIsVertexTPC[colIndex] > 0 && vIsVertexHighPtTPC[colIndex] == 0) { + float weightedTime = vWeightedTimesTPCnoTOFnoTRD[colIndex]; + float weightedSigma = vWeightedSigmaTPCnoTOFnoTRD[colIndex]; + auto bc = col.bc_as(); + int64_t globalBC = bc.globalBC(); + int64_t meanBC = globalBC + TMath::Nint(weightedTime / bcNS); + int64_t bestGlobalBC = findBestGlobalBC(meanBC, weightedSigma / bcNS, vNcontributors[colIndex], col.posZ(), mapGlobalBcVtxZ); + vFoundGlobalBC[colIndex] = bestGlobalBC > 0 ? bestGlobalBC : globalBC; + vFoundBCindex[colIndex] = bestGlobalBC > 0 ? mapGlobalBcWithTVX[bestGlobalBC] : bc.globalIndex(); + } + // fill pileup counter + vCollisionsPerBc[vFoundBCindex[colIndex]]++; } // save indices of collisions in time range for occupancy calculation @@ -817,7 +900,7 @@ struct EventSelectionTask { int nITS567tracksInTimeBins[nTimeIntervals] = {}; int nITS567tracksForVetoStandard = 0; // to veto events with nearby collisions int nITS567tracksForVetoNarrow = 0; // to veto events with nearby collisions (narrower range) - for (int iCol = 0; iCol < vAssocToThisCol.size(); iCol++) { + for (uint32_t iCol = 0; iCol < vAssocToThisCol.size(); iCol++) { int thisColIndex = vAssocToThisCol[iCol]; if (thisColIndex == colIndex) // skip the same collision continue;