Skip to content

Commit

Permalink
Improvements in collision-to-TVX matching (#7978)
Browse files Browse the repository at this point in the history
  • Loading branch information
@ekryshen
ekryshen authored Oct 14, 2024
1 parent c2b0985 commit 6e75b72
Showing 1 changed file with 166 additions and 83 deletions.
249 changes: 166 additions & 83 deletions Common/TableProducer/eventSelection.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,9 @@ MetadataHelper metadataInfo; // Metadata helper
using BCsWithRun2InfosTimestampsAndMatches = soa::Join<aod::BCs, aod::Run2BCInfos, aod::Timestamps, aod::Run2MatchedToBCSparse>;
using BCsWithRun3Matchings = soa::Join<aod::BCs, aod::Timestamps, aod::Run3MatchedToBCSparse>;
using BCsWithBcSelsRun2 = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels, aod::Run2BCInfos, aod::Run2MatchedToBCSparse>;
using BCsWithBcSelsRun3 = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels>;
using BCsWithBcSelsRun3 = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels, aod::Run3MatchedToBCSparse>;
using FullTracksIU = soa::Join<aod::TracksIU, aod::TracksExtra>;
const double bcNS = o2::constants::lhc::LHCBunchSpacingNS;

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

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

Expand All @@ -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<float> 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<int64_t, float>& 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<int64_t, float>::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)
Expand All @@ -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();
Expand All @@ -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<BCsWithBcSelsRun2>();
EventSelectionParams* par = ccdb->getForTimeStamp<EventSelectionParams>("EventSelection/EventSelectionParams", bc.timestamp());
Expand Down Expand Up @@ -589,8 +624,8 @@ struct EventSelectionTask {
}
PROCESS_SWITCH(EventSelectionTask, processRun2, "Process Run2 event selection", true);

Preslice<FullTracksIU> perCollision = aod::track::collisionId;
void processRun3(aod::Collisions const& cols, FullTracksIU const& tracks, BCsWithBcSelsRun3 const& bcs, aod::FT0s const&)
Partition<FullTracksIU> 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
Expand All @@ -600,7 +635,6 @@ struct EventSelectionTask {
auto grplhcif = ccdb->getForTimeStamp<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", ts);
bcPatternB = grplhcif->getBunchFilling().getBCPattern();

//
EventSelectionParams* par = ccdb->getForTimeStamp<EventSelectionParams>("EventSelection/EventSelectionParams", ts);
// access orbit-reset timestamp
auto ctpx = ccdb->getForTimeStamp<std::vector<Long64_t>>("CTP/Calib/OrbitReset", ts);
Expand All @@ -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<int64_t, int32_t> mapGlobalBcWithTVX;
std::map<int64_t, int32_t> mapGlobalBcWithTOR;
std::map<int64_t, float> 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<BCsWithBcSelsRun3>();
Expand All @@ -653,88 +685,139 @@ struct EventSelectionTask {
}
return;
}

std::vector<int> vFoundBCindex(cols.size(), -1); // indices of found bcs
std::vector<bool> vIsVertexITSTPC(cols.size(), 0); // at least one of vertex contributors is ITS-TPC track
std::vector<bool> vIsVertexTOFmatched(cols.size(), 0); // at least one of vertex contributors is matched to TOF
std::vector<bool> vIsVertexTRDmatched(cols.size(), 0); // at least one of vertex contributors is matched to TRD
std::vector<int> vCollisionsPerBc(bcs.size(), 0); // counter of collisions per found bc for pileup checks

// for the occupancy study
std::vector<int64_t> vFoundGlobalBC(cols.size(), 0); // global BCs for collisions
std::vector<int> vTracksITS567perColl(cols.size(), 0); // counter of tracks per found bc for occupancy studies
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)
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<bool> vIsVertexITSTPC(cols.size(), 0); // at least one of vertex contributors is ITS-TPC track
std::vector<bool> vIsVertexTOFmatched(cols.size(), 0); // at least one of vertex contributors is matched to TOF
std::vector<bool> vIsVertexTRDmatched(cols.size(), 0); // at least one of vertex contributors is matched to TRD

std::vector<int> vCollisionsPerBc(bcs.size(), 0); // counter of collisions per found bc for pileup checks
std::vector<int> vFoundBCindex(cols.size(), -1); // indices of found bcs
std::vector<int64_t> vFoundGlobalBC(cols.size(), 0); // global BCs for collisions

std::vector<bool> vIsVertexTOF(cols.size(), 0);
std::vector<bool> vIsVertexTRD(cols.size(), 0);
std::vector<bool> vIsVertexTPC(cols.size(), 0);
std::vector<bool> vIsVertexHighPtTPC(cols.size(), 0);
std::vector<int> vNcontributors(cols.size(), 0);
std::vector<float> vWeightedTimesTPCnoTOFnoTRD(cols.size(), 0);
std::vector<float> vWeightedSigmaTPCnoTOFnoTRD(cols.size(), 0);

// temporary vectors to find tracks with median time
std::vector<float> vTrackTimesTOF;
std::vector<float> vTrackTimesTRDnoTOF;

// first loop to match collisions to TVX
for (auto& col : cols) {
int32_t colIndex = col.globalIndex();
auto bc = col.bc_as<BCsWithBcSelsRun3>();
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<int64_t, int32_t>::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<int64_t, int32_t>::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<BCsWithBcSelsRun3>();
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
Expand Down Expand Up @@ -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;
Expand Down

0 comments on commit 6e75b72

Please sign in to comment.