Skip to content

Simplify and make sure that time resolution is kept into account #2404

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mpuccio merged 3 commits into from
Apr 14, 2023
Merged

Simplify and make sure that time resolution is kept into account #2404

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
27819ce
Simplify and make sure that time resolution is kept into account
mpuccio Apr 14, 2023
d616dba
Please consider the following formatting changes
alibuild Apr 14, 2023
7be2355
Use static cast
mpuccio Apr 14, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
235 changes: 62 additions & 173 deletions EventFiltering/selectBCRange.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
// O2 includes

#include "Common/DataModel/EventSelection.h"
#include "CommonConstants/LHCConstants.h"
Expand All @@ -22,229 +21,119 @@

using namespace o2;
using namespace o2::framework;
using namespace o2::framework::expressions;
using o2::dataformats::IRFrame;

// .............................................................................
// Run 3
struct BCRangeSelector {

Configurable<int> nTimeRes{"nTimeRes", 4, "Range to consider for search of compatible BCs in units of vertex-time-resolution."};
Configurable<int> nMinBCs{"nMinBCs", 7, "Minimum width of time window to consider for search of compatible BCs in units of 2*BunchSpacing."};
Configurable<double> fillFac{"fillFactor", 0.0, "Factor of MB events to add"};

using FDs = aod::CefpDecisions;
using CCs = soa::Join<aod::Collisions, aod::EvSels>;
using CC = CCs::iterator;
using BCs = soa::Join<aod::BCs, aod::BcSels, aod::Run3MatchedToBCSparse>;
using BC = BCs::iterator;

uint64_t nBCs, nCompBCs, nNotCompBCs;
uint64_t clast, cnew;

// buffer for task output
std::vector<o2::dataformats::IRFrame> res;
Produces<aod::BCRanges> tags;

void init(o2::framework::InitContext&)
template <typename T>
IRFrame getIRFrame(T& collision)
{
res.clear();
auto collBC = collision.bc().globalBC();
auto evSelBC = collision.has_foundBC() ? collision.foundBC().globalBC() : collBC;
int deltaBC = std::ceil(collision.collisionTimeRes() * nTimeRes / constants::lhc::LHCBunchSpacingNS);
deltaBC = std::max(deltaBC, nMinBCs.value);
IRFrame bcRange{InteractionRecord::long2IR(std::min(collBC, evSelBC)), InteractionRecord::long2IR(std::max(collBC, evSelBC))};
bcRange.getMax() += deltaBC;
bcRange.getMin() -= deltaBC;
return bcRange;
}

void run(ProcessingContext& pc)
{
// get and prepare the input tables
auto TabConsumer1 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::BCs>>::metadata::tableLabel());
auto t1{TabConsumer1->asArrowTable()};
auto TabConsumer2 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::BcSels>>::metadata::tableLabel());
auto t2{TabConsumer2->asArrowTable()};
auto TabConsumer3 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::Run3MatchedToBCSparse>>::metadata::tableLabel());
auto t3{TabConsumer3->asArrowTable()};
auto TabConsumer4 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::Collisions>>::metadata::tableLabel());
auto t4{TabConsumer4->asArrowTable()};
auto TabConsumer5 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::EvSels>>::metadata::tableLabel());
auto t5{TabConsumer5->asArrowTable()};
auto TabConsumer6 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::CefpDecisions>>::metadata::tableLabel());
auto t6{TabConsumer6->asArrowTable()};

// join tables
auto bcs = BCs({t1, t2, t3});
auto cols = CCs({t4, t5});
auto bcConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::BCs>>::metadata::tableLabel());
auto bcTable{bcConsumer->asArrowTable()};
auto collConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::Collisions>>::metadata::tableLabel());
auto collTable{collConsumer->asArrowTable()};
auto evSelConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::EvSels>>::metadata::tableLabel());
auto evSelTable{evSelConsumer->asArrowTable()};
auto cefpConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::CefpDecisions>>::metadata::tableLabel());
auto cefpTable{cefpConsumer->asArrowTable()};

auto bcs = aod::BCs({bcTable});
auto cols = CCs({collTable, evSelTable});
cols.bindExternalIndices(&bcs);
FDs fdecs{{t6}};
if (cols.size() != fdecs.size()) {
aod::CefpDecisions decisions{{cefpTable}};

if (cols.size() != decisions.size()) {
throw std::runtime_error("Collision table and CefpDecision do not have the same number of rows! ");
}
if (cols.size() == 0) {
LOGF(warning, "No collisions found!");
pc.outputs().snapshot({"PPF", "IFRAMES", 0, Lifetime::Timeframe}, res);
return;
}

// 1. loop over collisions
auto filt = fdecs.begin();
std::vector<std::pair<uint64_t, uint64_t>> bcRanges;
auto filt = decisions.begin();
int firstSelectedCollision{-1};
std::vector<IRFrame> bcRanges;
int nColl{0}, nSelected{0};
for (auto collision : cols) {
if (filt.cefpSelected()) {

LOGF(debug, "Collision time / resolution [ns]: %f / %f", collision.collisionTime(), collision.collisionTimeRes());

// return if collisions has no associated BC
if (!collision.has_foundBC()) {
LOGF(warning, "No compatible BCs found for collision that the framework assigned to BC %i", filt.globalBC());
filt++;
continue;
}

// get associated BC
auto bcIter = collision.foundBC_as<BCs>();

// due to the filling scheme the most probable BC may not be the one estimated from the collision time
InteractionRecord mostProbableBC;
mostProbableBC.setFromLong(bcIter.globalBC());
InteractionRecord meanBC = mostProbableBC + std::lround(collision.collisionTime() / o2::constants::lhc::LHCBunchSpacingNS);

// enforce minimum number for deltaBC
int deltaBC = std::ceil(collision.collisionTimeRes() * nTimeRes / o2::constants::lhc::LHCBunchSpacingNS);
if (deltaBC < nMinBCs) {
deltaBC = nMinBCs;
}
LOGF(debug, "BC %d, deltaBC %d", bcIter.globalIndex(), deltaBC);

auto minBC = meanBC - deltaBC;
auto maxBC = meanBC + deltaBC;

uint64_t minBCId = bcIter.globalIndex();
uint64_t maxBCId = bcIter.globalIndex();

auto localIter = bcIter;
while (localIter.globalIndex() > 0) {
--localIter;
if (localIter.globalBC() >= minBC.toLong()) {
minBCId = localIter.globalIndex();
} else {
break;
}
if (firstSelectedCollision < 0) {
firstSelectedCollision = nColl;
}
localIter = bcIter;
while (localIter.globalIndex() < bcs.size() - 1) {
++localIter;
if (localIter.globalBC() <= maxBC.toLong()) {
maxBCId = localIter.globalIndex();
} else {
break;
}
}

auto ao2dBC = filt.bcIndex();
if (ao2dBC < minBCId) {
LOGF(debug, "Extending the window to keep into account the AO2D/EvSel BC discrepancy by %d BCs", minBCId - ao2dBC);
minBCId = ao2dBC;
} else if (ao2dBC > maxBCId) {
LOGF(debug, "Extending the window to keep into account the AO2D/EvSel BC discrepancy by %d BCs", ao2dBC - maxBCId);
maxBCId = ao2dBC;
}

bcRanges.push_back(std::make_pair(minBCId, maxBCId));
bcRanges.push_back(getIRFrame(collision));
nSelected++;
}
nColl++;
filt++;
}

/// We cannot merge the ranges in the previous loop because while collisions are sorted by time, the corresponding minBCs can be unsorted as the collision time resolution is not constant
std::sort(bcRanges.begin(), bcRanges.end(), [](const std::pair<uint64_t, uint64_t>& a, const std::pair<uint64_t, uint64_t>& b) {
return a.first < b.first;
});
if (bcRanges.empty()) {
LOGF(warning, "No BCs selected! This should not happen! Adding a bogus BC range to avoid crashes.");
bcRanges.push_back(std::make_pair(0, 0));
LOGF(warning, "No BCs selected!");
return;
}

float fractionSelected{static_cast<float>(nSelected) / nColl};
int nMB{std::min(static_cast<int>(fillFac * nColl) - nSelected, nColl - 1)};
LOGF(info, "Selected %d collisions (%.2f%%) and %d MB events", nSelected, fractionSelected * 100, nMB);
int maxCollisionId = std::max(nMB, firstSelectedCollision);
int minCollisionId = (maxCollisionId == nMB) ? 0 : firstSelectedCollision - nMB;
auto minCollision = cols.begin() + minCollisionId;
IRFrame minFrame{getIRFrame(minCollision)};
bcRanges[0].getMin() = std::min(bcRanges[0].getMin(), minFrame.getMin());
if (maxCollisionId == nMB) {
auto maxCollision = cols.begin() + nMB;
IRFrame maxFrame{getIRFrame(maxCollision)};
bcRanges[0].getMax() = std::max(bcRanges[0].getMax(), maxFrame.getMax());
}
std::vector<std::pair<uint64_t, uint64_t>> bcRangesMerged(1, bcRanges[0]);

/// We cannot merge the ranges in the previous loop because while collisions are sorted by time, the corresponding minBCs can be unsorted as the collision time resolution is not constant
std::sort(bcRanges.begin(), bcRanges.end(), [](const IRFrame& a, const IRFrame& b) {
return a.getMin() < b.getMin();
});

std::vector<IRFrame> bcRangesMerged(1, bcRanges[0]);
for (uint64_t iR{1}; iR < bcRanges.size(); ++iR) {
if (bcRangesMerged.back().second >= bcRanges[iR].first) {
bcRangesMerged.back().second = std::max(bcRangesMerged.back().second, bcRanges[iR].second);
if (bcRangesMerged.back().getMax() >= bcRanges[iR].getMin()) {
bcRangesMerged.back().getMax() = std::max(bcRangesMerged.back().getMax(), bcRanges[iR].getMax());
} else {
bcRangesMerged.push_back(bcRanges[iR]);
}
}
bcRanges.swap(bcRangesMerged);

// 2. extend ranges
int nBCselected{0};
for (auto& range : bcRanges) {
nBCselected += range.second - range.first + 1;
tags(range.getMin().toLong(), range.getMax().toLong());
}
int nToBeAdded = std::ceil((bcs.size() - nBCselected) * fillFac);
int nToBeAddedPerRange = std::ceil(float(nToBeAdded) / bcRanges.size());
LOGF(debug, "Extending ranges by %d BCs (%d per selected range)", nToBeAdded, nToBeAddedPerRange);

InteractionRecord IR1, IR2;
uint64_t first{bcs.iteratorAt(bcRanges[0].first).globalBC()};
IR1.setFromLong(first);
while (nToBeAdded > 0 && bcRanges[0].first > 0) { /// TODO: decide if we want to extend the ranges to the beginning of the dataframe
first = bcs.iteratorAt(bcRanges[0].first - 1).globalBC();
IR2.setFromLong(first);
if (IR1.differenceInBC(IR2) > o2::constants::lhc::LHCMaxBunches) { // protection against change of orbit in the DataFrame
LOGF(debug, "Jump by more than one orbit detected!");
break;
} else {
IR1.setFromLong(first);
}
bcRanges[0].first--;
nToBeAdded--;
}
LOGF(debug, "Extending ranges by %d BCs (%d per selected range)", nToBeAdded, nToBeAddedPerRange);

for (uint64_t iR{0}; iR < bcRanges.size() && nToBeAdded > 0; ++iR) {
uint64_t second{bcs.iteratorAt(bcRanges[iR].second).globalBC()};
IR2.setFromLong(second);
for (int i{0}; i < nToBeAddedPerRange && nToBeAdded > 0; ++i) {
if (bcRanges[iR].second < bcs.size() - 1) {
second = bcs.iteratorAt(bcRanges[iR].second + 1).globalBC();
IR1.setFromLong(second);
if (IR1.differenceInBC(IR2) > o2::constants::lhc::LHCMaxBunches) { // protection against change of orbit in the DataFrame
LOGF(debug, "Jump by more than one orbit detected!");
break;
} else {
IR2.setFromLong(second);
}
bcRanges[iR].second++;
nToBeAdded--;
}
}
}
LOGF(debug, "End extension, remaining to be added %d BCs", nToBeAdded);

// fill res
LOGF(debug, "Merged and extended sorted ranges");
for (auto& range : bcRanges) {
LOGF(debug, " %i - %i", range.first, range.second);
uint64_t first{bcs.rawIteratorAt(range.first).globalBC()}, second{bcs.rawIteratorAt(range.second).globalBC()};
IR1.setFromLong(first);
IR2.setFromLong(second);
res.emplace_back(IR1, IR2);
tags(first, second);
}

// make res an output
pc.outputs().snapshot({"PPF", "IFRAMES", 0, Lifetime::Timeframe}, res);

// clean up
res.clear();
}

// need a trivial process method
// the parameters determine the tables available in the input
void process(CCs const& collisions, BCs const& bcs, FDs const& fdecisions)
// need a trivial process method: the parameters determine the tables available in the input
void process(CCs const& collisions, aod::BCs const& bcs, aod::CefpDecisions const& fdecisions)
{
}
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
DataProcessorSpec spec{adaptAnalysisTask<BCRangeSelector>(cfgc, TaskName{"bc-ranges-selector-task"})};

// add output
spec.outputs.emplace_back("PPF", "IFRAMES", 0, Lifetime::Timeframe);
LOGF(debug, "Output %i", spec.outputs.size());

return WorkflowSpec{spec};
}