Skip to content

Commit 6e75b72

Browse files
ekryshenekryshen
authored
Improvements in collision-to-TVX matching (AliceO2Group#7978)
1 parent c2b0985 commit 6e75b72

File tree

1 file changed

+166
-83
lines changed

1 file changed

+166
-83
lines changed

‎Common/TableProducer/eventSelection.cxx‎

Lines changed: 166 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -36,8 +36,9 @@ MetadataHelper metadataInfo; // Metadata helper
3636
using BCsWithRun2InfosTimestampsAndMatches = soa::Join<aod::BCs, aod::Run2BCInfos, aod::Timestamps, aod::Run2MatchedToBCSparse>;
3737
using BCsWithRun3Matchings = soa::Join<aod::BCs, aod::Timestamps, aod::Run3MatchedToBCSparse>;
3838
using BCsWithBcSelsRun2 = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels, aod::Run2BCInfos, aod::Run2MatchedToBCSparse>;
39-
using BCsWithBcSelsRun3 = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels>;
39+
using BCsWithBcSelsRun3 = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels, aod::Run3MatchedToBCSparse>;
4040
using FullTracksIU = soa::Join<aod::TracksIU, aod::TracksExtra>;
41+
const double bcNS = o2::constants::lhc::LHCBunchSpacingNS;
4142

4243
struct BcSelectionTask {
4344
Produces<aod::BcSels> bcsel;
@@ -454,6 +455,7 @@ struct EventSelectionTask {
454455
Configurable<float> confTimeRangeVetoOnCollStandard{"TimeRangeVetoOnCollStandard", 10, "Exclusion of a collision if there are other collisions nearby, +/- us"};
455456
Configurable<float> confTimeRangeVetoOnCollNarrow{"TimeRangeVetoOnCollNarrow", 4, "Exclusion of a collision if there are other collisions nearby, +/- us"};
456457
Configurable<bool> confUseWeightsForOccupancyVariable{"UseWeightsForOccupancyEstimator", 1, "Use or not the delta-time weights for the occupancy estimator"};
458+
Configurable<int> confSigmaBCforHighPtTracks{"confSigmaBCforHighPtTracks", 4, "Custom sigma (in bcs) for collisions with high-pt tracks"};
457459

458460
Partition<aod::Tracks> tracklets = (aod::track::trackType == static_cast<uint8_t>(o2::aod::track::TrackTypeEnum::Run2Tracklet));
459461

@@ -480,6 +482,40 @@ struct EventSelectionTask {
480482
return (dbc1 <= dbc2) ? index1 : index2;
481483
}
482484

485+
// helper function to find median time in the vector of TOF or TRD-track times
486+
float getMedian(std::vector<float> v)
487+
{
488+
int medianIndex = v.size() / 2;
489+
std::nth_element(v.begin(), v.begin() + medianIndex, v.end());
490+
return v[medianIndex];
491+
}
492+
493+
// helper function to find closest TVX signal in time and in zVtx
494+
int64_t findBestGlobalBC(int64_t meanBC, int64_t sigmaBC, int32_t nContrib, float zVtxCol, std::map<int64_t, float>& mapGlobalBcVtxZ)
495+
{
496+
int64_t minBC = meanBC - 3 * sigmaBC;
497+
int64_t maxBC = meanBC + 3 * sigmaBC;
498+
// TODO: use ITS ROF bounds to reduce the search range?
499+
500+
float zVtxSigma = 2.7 * pow(nContrib, -0.466) + 0.024;
501+
zVtxSigma += 1.0; // additional uncertainty due to imperfectections of FT0 time calibration
502+
503+
auto itMin = mapGlobalBcVtxZ.lower_bound(minBC);
504+
auto itMax = mapGlobalBcVtxZ.upper_bound(maxBC);
505+
506+
float bestChi2 = 1e+10;
507+
int64_t bestGlobalBC = 0;
508+
for (std::map<int64_t, float>::iterator it = itMin; it != itMax; ++it) {
509+
float chi2 = pow((it->second - zVtxCol) / zVtxSigma, 2) + pow((it->first - meanBC) / sigmaBC, 2.);
510+
if (chi2 < bestChi2) {
511+
bestChi2 = chi2;
512+
bestGlobalBC = it->first;
513+
}
514+
}
515+
516+
return bestGlobalBC;
517+
}
518+
483519
void init(InitContext&)
484520
{
485521
if (metadataInfo.isFullyDefined()) { // Check if the metadata is initialized (only if not forced from the workflow configuration)
@@ -501,7 +537,6 @@ struct EventSelectionTask {
501537
}
502538
}
503539

504-
// ccdb->setURL("http://ccdb-test.cern.ch:8080");
505540
ccdb->setURL("http://alice-ccdb.cern.ch");
506541
ccdb->setCaching(true);
507542
ccdb->setLocalObjectValidityChecking();
@@ -517,7 +552,7 @@ struct EventSelectionTask {
517552
evsel.reserve(collisions.size());
518553
}
519554

520-
void processRun2(aod::Collision const& col, BCsWithBcSelsRun2 const&, aod::Tracks const&, aod::FV0Cs const&)
555+
void processRun2(aod::Collision const& col, BCsWithBcSelsRun2 const&, FullTracksIU const&, aod::FV0Cs const&)
521556
{
522557
auto bc = col.bc_as<BCsWithBcSelsRun2>();
523558
EventSelectionParams* par = ccdb->getForTimeStamp<EventSelectionParams>("EventSelection/EventSelectionParams", bc.timestamp());
@@ -589,8 +624,8 @@ struct EventSelectionTask {
589624
}
590625
PROCESS_SWITCH(EventSelectionTask, processRun2, "Process Run2 event selection", true);
591626

592-
Preslice<FullTracksIU> perCollision = aod::track::collisionId;
593-
void processRun3(aod::Collisions const& cols, FullTracksIU const& tracks, BCsWithBcSelsRun3 const& bcs, aod::FT0s const&)
627+
Partition<FullTracksIU> pvTracks = ((aod::track::flags & (uint32_t)o2::aod::track::PVContributor) == (uint32_t)o2::aod::track::PVContributor);
628+
void processRun3(aod::Collisions const& cols, FullTracksIU const&, BCsWithBcSelsRun3 const& bcs, aod::FT0s const&)
594629
{
595630
int run = bcs.iteratorAt(0).runNumber();
596631
// extract bc pattern from CCDB for data or anchored MC only
@@ -600,7 +635,6 @@ struct EventSelectionTask {
600635
auto grplhcif = ccdb->getForTimeStamp<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", ts);
601636
bcPatternB = grplhcif->getBunchFilling().getBCPattern();
602637

603-
//
604638
EventSelectionParams* par = ccdb->getForTimeStamp<EventSelectionParams>("EventSelection/EventSelectionParams", ts);
605639
// access orbit-reset timestamp
606640
auto ctpx = ccdb->getForTimeStamp<std::vector<Long64_t>>("CTP/Calib/OrbitReset", ts);
@@ -621,26 +655,24 @@ struct EventSelectionTask {
621655
nBCsPerTF = nOrbitsPerTF * o2::constants::lhc::LHCMaxBunches;
622656
}
623657

624-
// create maps from globalBC to bc index for TVX or FT0-OR fired bcs
625-
// to be used for closest TVX (FT0-OR) searches
658+
// create maps from globalBC to bc index for TVX-fired bcs
659+
// to be used for closest TVX searches
626660
std::map<int64_t, int32_t> mapGlobalBcWithTVX;
627-
std::map<int64_t, int32_t> mapGlobalBcWithTOR;
661+
std::map<int64_t, float> mapGlobalBcVtxZ;
628662
for (auto& bc : bcs) {
629663
int64_t globalBC = bc.globalBC();
630664
// skip non-colliding bcs for data and anchored runs
631665
if (run >= 500000 && bcPatternB[globalBC % o2::constants::lhc::LHCMaxBunches] == 0) {
632666
continue;
633667
}
634-
if (bc.selection_bit(kIsBBT0A) || bc.selection_bit(kIsBBT0C)) {
635-
mapGlobalBcWithTOR[globalBC] = bc.globalIndex();
636-
}
637668
if (bc.selection_bit(kIsTriggerTVX)) {
638669
mapGlobalBcWithTVX[globalBC] = bc.globalIndex();
670+
mapGlobalBcVtxZ[globalBC] = bc.has_ft0() ? bc.ft0().posZ() : 0;
639671
}
640672
}
641673

642674
// protection against empty FT0 maps
643-
if (mapGlobalBcWithTOR.size() == 0 || mapGlobalBcWithTVX.size() == 0) {
675+
if (mapGlobalBcWithTVX.size() == 0) {
644676
LOGP(error, "FT0 table is empty or corrupted. Filling evsel table with dummy values");
645677
for (auto& col : cols) {
646678
auto bc = col.bc_as<BCsWithBcSelsRun3>();
@@ -653,88 +685,139 @@ struct EventSelectionTask {
653685
}
654686
return;
655687
}
656-
657-
std::vector<int> vFoundBCindex(cols.size(), -1); // indices of found bcs
658-
std::vector<bool> vIsVertexITSTPC(cols.size(), 0); // at least one of vertex contributors is ITS-TPC track
659-
std::vector<bool> vIsVertexTOFmatched(cols.size(), 0); // at least one of vertex contributors is matched to TOF
660-
std::vector<bool> vIsVertexTRDmatched(cols.size(), 0); // at least one of vertex contributors is matched to TRD
661-
std::vector<int> vCollisionsPerBc(bcs.size(), 0); // counter of collisions per found bc for pileup checks
662-
663-
// for the occupancy study
664-
std::vector<int64_t> vFoundGlobalBC(cols.size(), 0); // global BCs for collisions
665688
std::vector<int> vTracksITS567perColl(cols.size(), 0); // counter of tracks per found bc for occupancy studies
666689
std::vector<bool> vIsFullInfoForOccupancy(cols.size(), 0); // info for occupancy in +/- windows is available (i.e. a given coll is not too close to the TF borders)
667690
const float timeWinOccupancyCalcMinNS = confTimeIntervalForOccupancyCalculationMin * 1e3; // ns
668691
const float timeWinOccupancyCalcMaxNS = confTimeIntervalForOccupancyCalculationMax * 1e3; // ns
669-
// const double timeWinOccupancyExclusionRangeNS = confExclusionIntervalForOccupancyCalculation * 1e3; // ns
670-
const double bcNS = o2::constants::lhc::LHCBunchSpacingNS;
671-
672-
// loop to find nearest bc with FT0 entry -> foundBC index
692+
std::vector<bool> vIsVertexITSTPC(cols.size(), 0); // at least one of vertex contributors is ITS-TPC track
693+
std::vector<bool> vIsVertexTOFmatched(cols.size(), 0); // at least one of vertex contributors is matched to TOF
694+
std::vector<bool> vIsVertexTRDmatched(cols.size(), 0); // at least one of vertex contributors is matched to TRD
695+
696+
std::vector<int> vCollisionsPerBc(bcs.size(), 0); // counter of collisions per found bc for pileup checks
697+
std::vector<int> vFoundBCindex(cols.size(), -1); // indices of found bcs
698+
std::vector<int64_t> vFoundGlobalBC(cols.size(), 0); // global BCs for collisions
699+
700+
std::vector<bool> vIsVertexTOF(cols.size(), 0);
701+
std::vector<bool> vIsVertexTRD(cols.size(), 0);
702+
std::vector<bool> vIsVertexTPC(cols.size(), 0);
703+
std::vector<bool> vIsVertexHighPtTPC(cols.size(), 0);
704+
std::vector<int> vNcontributors(cols.size(), 0);
705+
std::vector<float> vWeightedTimesTPCnoTOFnoTRD(cols.size(), 0);
706+
std::vector<float> vWeightedSigmaTPCnoTOFnoTRD(cols.size(), 0);
707+
708+
// temporary vectors to find tracks with median time
709+
std::vector<float> vTrackTimesTOF;
710+
std::vector<float> vTrackTimesTRDnoTOF;
711+
712+
// first loop to match collisions to TVX
673713
for (auto& col : cols) {
714+
int32_t colIndex = col.globalIndex();
674715
auto bc = col.bc_as<BCsWithBcSelsRun3>();
675-
int64_t meanBC = bc.globalBC();
676-
const double bcNS = o2::constants::lhc::LHCBunchSpacingNS;
677-
int64_t deltaBC = std::ceil(col.collisionTimeRes() / bcNS * 4);
678-
679-
// count tracks of different types
680-
int nITS567cls = 0;
681-
int nITSTPCtracks = 0;
682-
int nTOFtracks = 0;
683-
int nTRDtracks = 0;
684-
double timeFromTOFtracks = 0;
685-
auto tracksGrouped = tracks.sliceBy(perCollision, col.globalIndex());
686-
for (auto& track : tracksGrouped) {
687-
if (!track.isPVContributor()) {
688-
continue;
689-
}
690-
nITSTPCtracks += track.hasITS() && track.hasTPC();
691-
nTOFtracks += track.hasTOF();
692-
nTRDtracks += track.hasTRD();
693-
// calculate average time using TOF tracks
694-
if (track.hasTOF()) {
695-
timeFromTOFtracks += track.trackTime();
696-
}
716+
int64_t globalBC = bc.globalBC();
717+
int64_t bcInTF = (bc.globalBC() - bcSOR) % nBCsPerTF;
718+
vIsFullInfoForOccupancy[colIndex] = ((bcInTF - 300) * bcNS > -timeWinOccupancyCalcMinNS) && ((nBCsPerTF - 4000 - bcInTF) * bcNS > timeWinOccupancyCalcMaxNS) ? true : false;
697719

720+
const auto& colPvTracks = pvTracks.sliceByCached(aod::track::collisionId, col.globalIndex(), cache);
721+
vTrackTimesTOF.clear();
722+
vTrackTimesTRDnoTOF.clear();
723+
int nPvTracksTPCnoTOFnoTRD = 0;
724+
int nPvTracksHighPtTPCnoTOFnoTRD = 0;
725+
float sumTime = 0, sumW = 0, sumHighPtTime = 0, sumHighPtW = 0;
726+
for (auto& track : colPvTracks) {
727+
float trackTime = track.trackTime();
698728
if (track.itsNCls() >= 5)
699-
nITS567cls++;
729+
vTracksITS567perColl[colIndex]++;
730+
if (track.hasTRD())
731+
vIsVertexTRDmatched[colIndex] = 1;
732+
if (track.hasTPC())
733+
vIsVertexITSTPC[colIndex] = 1;
734+
if (track.hasTOF()) {
735+
vTrackTimesTOF.push_back(trackTime);
736+
vIsVertexTOFmatched[colIndex] = 1;
737+
} else if (track.hasTRD()) {
738+
vTrackTimesTRDnoTOF.push_back(trackTime);
739+
} else if (track.hasTPC()) {
740+
float trackTimeRes = track.trackTimeRes();
741+
float trackPt = track.pt();
742+
float w = 1. / (trackTimeRes * trackTimeRes);
743+
sumTime += trackTime * w;
744+
sumW += w;
745+
nPvTracksTPCnoTOFnoTRD++;
746+
if (trackPt > 1) {
747+
sumHighPtTime += trackTime * w;
748+
sumHighPtW += w;
749+
nPvTracksHighPtTPCnoTOFnoTRD++;
750+
}
751+
}
700752
}
701-
LOGP(debug, "nContrib={} nITSTPCtracks={} nTOFtracks={} nTRDtracks={}", col.numContrib(), nITSTPCtracks, nTOFtracks, nTRDtracks);
702-
703-
if (nTOFtracks > 0) {
704-
meanBC += TMath::FloorNint(timeFromTOFtracks / nTOFtracks / bcNS); // assign collision bc using TOF-matched tracks
705-
deltaBC = 4; // use precise bc from TOF tracks with +/-4 bc margin
706-
} else if (nITSTPCtracks > 0) {
707-
deltaBC += 30; // extend deltaBC for collisions built with ITS-TPC tracks only
753+
vWeightedTimesTPCnoTOFnoTRD[colIndex] = sumW > 0 ? sumTime / sumW : 0;
754+
vWeightedSigmaTPCnoTOFnoTRD[colIndex] = sumW > 0 ? sqrt(1. / sumW) : 0;
755+
vNcontributors[colIndex] = colPvTracks.size();
756+
int nPvTracksTOF = vTrackTimesTOF.size();
757+
int nPvTracksTRDnoTOF = vTrackTimesTRDnoTOF.size();
758+
// collision type
759+
vIsVertexTOF[colIndex] = nPvTracksTOF > 0;
760+
vIsVertexTRD[colIndex] = nPvTracksTRDnoTOF > 0;
761+
vIsVertexTPC[colIndex] = nPvTracksTPCnoTOFnoTRD > 0;
762+
vIsVertexHighPtTPC[colIndex] = nPvTracksHighPtTPCnoTOFnoTRD > 0;
763+
764+
int64_t foundGlobalBC = 0;
765+
int32_t foundBCindex = -1;
766+
767+
if (nPvTracksTOF > 0) {
768+
// for collisions with TOF tracks:
769+
// take bc corresponding to TOF track with median time
770+
int64_t tofGlobalBC = globalBC + TMath::Nint(getMedian(vTrackTimesTOF) / bcNS);
771+
std::map<int64_t, int32_t>::iterator it = mapGlobalBcWithTVX.find(tofGlobalBC);
772+
if (it != mapGlobalBcWithTVX.end()) {
773+
foundGlobalBC = it->first;
774+
foundBCindex = it->second;
775+
}
776+
} else if (nPvTracksTPCnoTOFnoTRD == 0 && nPvTracksTRDnoTOF > 0) {
777+
// for collisions with TRD tracks but without TOF or ITSTPC-only tracks:
778+
// take bc corresponding to TRD track with median time
779+
int64_t trdGlobalBC = globalBC + TMath::Nint(getMedian(vTrackTimesTRDnoTOF) / bcNS);
780+
std::map<int64_t, int32_t>::iterator it = mapGlobalBcWithTVX.find(trdGlobalBC);
781+
if (it != mapGlobalBcWithTVX.end()) {
782+
foundGlobalBC = it->first;
783+
foundBCindex = it->second;
784+
}
785+
} else if (nPvTracksHighPtTPCnoTOFnoTRD > 0) {
786+
// for collisions with high-pt ITSTPC-nonTOF-nonTRD tracks
787+
// search in 3*confSigmaBCforHighPtTracks range (3*4 bcs by default)
788+
int64_t meanBC = globalBC + TMath::Nint(sumHighPtTime / sumHighPtW / bcNS);
789+
int64_t bestGlobalBC = findBestGlobalBC(meanBC, confSigmaBCforHighPtTracks, vNcontributors[colIndex], col.posZ(), mapGlobalBcVtxZ);
790+
if (bestGlobalBC > 0) {
791+
foundGlobalBC = bestGlobalBC;
792+
foundBCindex = mapGlobalBcWithTVX[bestGlobalBC];
793+
}
708794
}
709795

710-
int64_t minBC = meanBC - deltaBC;
711-
int64_t maxBC = meanBC + deltaBC;
796+
// fill foundBC indices and global BCs
797+
// keep current bc if TVX matching failed at this step
798+
vFoundBCindex[colIndex] = foundBCindex >= 0 ? foundBCindex : bc.globalIndex();
799+
vFoundGlobalBC[colIndex] = foundGlobalBC > 0 ? foundGlobalBC : globalBC;
712800

713-
int32_t indexClosestTVX = findClosest(meanBC, mapGlobalBcWithTVX);
714-
int64_t tvxBC = bcs.iteratorAt(indexClosestTVX).globalBC();
715-
if (tvxBC >= minBC && tvxBC <= maxBC) { // closest TVX within search region
716-
bc.setCursor(indexClosestTVX);
717-
} else { // no TVX within search region, searching for TOR = T0A | T0C
718-
int32_t indexClosestTOR = findClosest(meanBC, mapGlobalBcWithTOR);
719-
int64_t torBC = bcs.iteratorAt(indexClosestTOR).globalBC();
720-
if (torBC >= minBC && torBC <= maxBC) {
721-
bc.setCursor(indexClosestTOR);
722-
}
723-
}
724-
int32_t foundBC = bc.globalIndex();
801+
// erase found global BC with TVX from the pool of bcs for the next loop over low-pt TPCnoTOFnoTRD collisions
802+
if (foundBCindex >= 0)
803+
mapGlobalBcVtxZ.erase(foundGlobalBC);
804+
}
805+
806+
// second loop to match remaining low-pt TPCnoTOFnoTRD collisions
807+
for (auto& col : cols) {
725808
int32_t colIndex = col.globalIndex();
726-
LOGP(debug, "foundBC = {} globalBC = {}", foundBC, bc.globalBC());
727-
vFoundBCindex[colIndex] = foundBC;
728-
vIsVertexITSTPC[colIndex] = nITSTPCtracks > 0;
729-
vIsVertexTOFmatched[colIndex] = nTOFtracks > 0;
730-
vIsVertexTRDmatched[colIndex] = nTRDtracks > 0;
731-
vCollisionsPerBc[foundBC]++;
732-
vTracksITS567perColl[colIndex] = nITS567cls;
733-
vFoundGlobalBC[colIndex] = bc.globalBC();
734-
735-
// check that this collision has full information inside the time window (taking into account TF borders)
736-
int64_t bcInTF = (bc.globalBC() - bcSOR) % nBCsPerTF;
737-
vIsFullInfoForOccupancy[colIndex] = ((bcInTF - 300) * bcNS > -timeWinOccupancyCalcMinNS) && ((nBCsPerTF - 4000 - bcInTF) * bcNS > timeWinOccupancyCalcMaxNS) ? true : false;
809+
if (vIsVertexTPC[colIndex] > 0 && vIsVertexHighPtTPC[colIndex] == 0) {
810+
float weightedTime = vWeightedTimesTPCnoTOFnoTRD[colIndex];
811+
float weightedSigma = vWeightedSigmaTPCnoTOFnoTRD[colIndex];
812+
auto bc = col.bc_as<BCsWithBcSelsRun3>();
813+
int64_t globalBC = bc.globalBC();
814+
int64_t meanBC = globalBC + TMath::Nint(weightedTime / bcNS);
815+
int64_t bestGlobalBC = findBestGlobalBC(meanBC, weightedSigma / bcNS, vNcontributors[colIndex], col.posZ(), mapGlobalBcVtxZ);
816+
vFoundGlobalBC[colIndex] = bestGlobalBC > 0 ? bestGlobalBC : globalBC;
817+
vFoundBCindex[colIndex] = bestGlobalBC > 0 ? mapGlobalBcWithTVX[bestGlobalBC] : bc.globalIndex();
818+
}
819+
// fill pileup counter
820+
vCollisionsPerBc[vFoundBCindex[colIndex]]++;
738821
}
739822

740823
// save indices of collisions in time range for occupancy calculation
@@ -817,7 +900,7 @@ struct EventSelectionTask {
817900
int nITS567tracksInTimeBins[nTimeIntervals] = {};
818901
int nITS567tracksForVetoStandard = 0; // to veto events with nearby collisions
819902
int nITS567tracksForVetoNarrow = 0; // to veto events with nearby collisions (narrower range)
820-
for (int iCol = 0; iCol < vAssocToThisCol.size(); iCol++) {
903+
for (uint32_t iCol = 0; iCol < vAssocToThisCol.size(); iCol++) {
821904
int thisColIndex = vAssocToThisCol[iCol];
822905
if (thisColIndex == colIndex) // skip the same collision
823906
continue;

0 commit comments

Comments
 (0)