Skip to content

Commit 4644239

Browse files
mpucciompuccio
authored
Simplify and make sure that time resolution is kept into account (#2404)
1 parent 46d545e commit 4644239

1 file changed

Lines changed: 62 additions & 173 deletions

File tree

‎EventFiltering/selectBCRange.cxx‎

Lines changed: 62 additions & 173 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,6 @@
88
// In applying this license CERN does not waive the privileges and immunities
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
11-
// O2 includes
1211

1312
#include "Common/DataModel/EventSelection.h"
1413
#include "CommonConstants/LHCConstants.h"
@@ -22,229 +21,119 @@
2221

2322
using namespace o2;
2423
using namespace o2::framework;
25-
using namespace o2::framework::expressions;
24+
using o2::dataformats::IRFrame;
2625

27-
// .............................................................................
28-
// Run 3
2926
struct BCRangeSelector {
3027

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

35-
using FDs = aod::CefpDecisions;
3632
using CCs = soa::Join<aod::Collisions, aod::EvSels>;
37-
using CC = CCs::iterator;
38-
using BCs = soa::Join<aod::BCs, aod::BcSels, aod::Run3MatchedToBCSparse>;
39-
using BC = BCs::iterator;
40-
41-
uint64_t nBCs, nCompBCs, nNotCompBCs;
42-
uint64_t clast, cnew;
4333

4434
// buffer for task output
45-
std::vector<o2::dataformats::IRFrame> res;
4635
Produces<aod::BCRanges> tags;
4736

48-
void init(o2::framework::InitContext&)
37+
template <typename T>
38+
IRFrame getIRFrame(T& collision)
4939
{
50-
res.clear();
40+
auto collBC = collision.bc().globalBC();
41+
auto evSelBC = collision.has_foundBC() ? collision.foundBC().globalBC() : collBC;
42+
int deltaBC = std::ceil(collision.collisionTimeRes() * nTimeRes / constants::lhc::LHCBunchSpacingNS);
43+
deltaBC = std::max(deltaBC, nMinBCs.value);
44+
IRFrame bcRange{InteractionRecord::long2IR(std::min(collBC, evSelBC)), InteractionRecord::long2IR(std::max(collBC, evSelBC))};
45+
bcRange.getMax() += deltaBC;
46+
bcRange.getMin() -= deltaBC;
47+
return bcRange;
5148
}
5249

5350
void run(ProcessingContext& pc)
5451
{
55-
// get and prepare the input tables
56-
auto TabConsumer1 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::BCs>>::metadata::tableLabel());
57-
auto t1{TabConsumer1->asArrowTable()};
58-
auto TabConsumer2 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::BcSels>>::metadata::tableLabel());
59-
auto t2{TabConsumer2->asArrowTable()};
60-
auto TabConsumer3 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::Run3MatchedToBCSparse>>::metadata::tableLabel());
61-
auto t3{TabConsumer3->asArrowTable()};
62-
auto TabConsumer4 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::Collisions>>::metadata::tableLabel());
63-
auto t4{TabConsumer4->asArrowTable()};
64-
auto TabConsumer5 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::EvSels>>::metadata::tableLabel());
65-
auto t5{TabConsumer5->asArrowTable()};
66-
auto TabConsumer6 = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::CefpDecisions>>::metadata::tableLabel());
67-
auto t6{TabConsumer6->asArrowTable()};
68-
69-
// join tables
70-
auto bcs = BCs({t1, t2, t3});
71-
auto cols = CCs({t4, t5});
52+
auto bcConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::BCs>>::metadata::tableLabel());
53+
auto bcTable{bcConsumer->asArrowTable()};
54+
auto collConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::Collisions>>::metadata::tableLabel());
55+
auto collTable{collConsumer->asArrowTable()};
56+
auto evSelConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::EvSels>>::metadata::tableLabel());
57+
auto evSelTable{evSelConsumer->asArrowTable()};
58+
auto cefpConsumer = pc.inputs().get<TableConsumer>(aod::MetadataTrait<std::decay_t<aod::CefpDecisions>>::metadata::tableLabel());
59+
auto cefpTable{cefpConsumer->asArrowTable()};
60+
61+
auto bcs = aod::BCs({bcTable});
62+
auto cols = CCs({collTable, evSelTable});
7263
cols.bindExternalIndices(&bcs);
73-
FDs fdecs{{t6}};
74-
if (cols.size() != fdecs.size()) {
64+
aod::CefpDecisions decisions{{cefpTable}};
65+
66+
if (cols.size() != decisions.size()) {
7567
throw std::runtime_error("Collision table and CefpDecision do not have the same number of rows! ");
7668
}
7769
if (cols.size() == 0) {
7870
LOGF(warning, "No collisions found!");
79-
pc.outputs().snapshot({"PPF", "IFRAMES", 0, Lifetime::Timeframe}, res);
8071
return;
8172
}
8273

83-
// 1. loop over collisions
84-
auto filt = fdecs.begin();
85-
std::vector<std::pair<uint64_t, uint64_t>> bcRanges;
74+
auto filt = decisions.begin();
75+
int firstSelectedCollision{-1};
76+
std::vector<IRFrame> bcRanges;
77+
int nColl{0}, nSelected{0};
8678
for (auto collision : cols) {
8779
if (filt.cefpSelected()) {
88-
89-
LOGF(debug, "Collision time / resolution [ns]: %f / %f", collision.collisionTime(), collision.collisionTimeRes());
90-
91-
// return if collisions has no associated BC
92-
if (!collision.has_foundBC()) {
93-
LOGF(warning, "No compatible BCs found for collision that the framework assigned to BC %i", filt.globalBC());
94-
filt++;
95-
continue;
96-
}
97-
98-
// get associated BC
99-
auto bcIter = collision.foundBC_as<BCs>();
100-
101-
// due to the filling scheme the most probable BC may not be the one estimated from the collision time
102-
InteractionRecord mostProbableBC;
103-
mostProbableBC.setFromLong(bcIter.globalBC());
104-
InteractionRecord meanBC = mostProbableBC + std::lround(collision.collisionTime() / o2::constants::lhc::LHCBunchSpacingNS);
105-
106-
// enforce minimum number for deltaBC
107-
int deltaBC = std::ceil(collision.collisionTimeRes() * nTimeRes / o2::constants::lhc::LHCBunchSpacingNS);
108-
if (deltaBC < nMinBCs) {
109-
deltaBC = nMinBCs;
110-
}
111-
LOGF(debug, "BC %d, deltaBC %d", bcIter.globalIndex(), deltaBC);
112-
113-
auto minBC = meanBC - deltaBC;
114-
auto maxBC = meanBC + deltaBC;
115-
116-
uint64_t minBCId = bcIter.globalIndex();
117-
uint64_t maxBCId = bcIter.globalIndex();
118-
119-
auto localIter = bcIter;
120-
while (localIter.globalIndex() > 0) {
121-
--localIter;
122-
if (localIter.globalBC() >= minBC.toLong()) {
123-
minBCId = localIter.globalIndex();
124-
} else {
125-
break;
126-
}
80+
if (firstSelectedCollision < 0) {
81+
firstSelectedCollision = nColl;
12782
}
128-
localIter = bcIter;
129-
while (localIter.globalIndex() < bcs.size() - 1) {
130-
++localIter;
131-
if (localIter.globalBC() <= maxBC.toLong()) {
132-
maxBCId = localIter.globalIndex();
133-
} else {
134-
break;
135-
}
136-
}
137-
138-
auto ao2dBC = filt.bcIndex();
139-
if (ao2dBC < minBCId) {
140-
LOGF(debug, "Extending the window to keep into account the AO2D/EvSel BC discrepancy by %d BCs", minBCId - ao2dBC);
141-
minBCId = ao2dBC;
142-
} else if (ao2dBC > maxBCId) {
143-
LOGF(debug, "Extending the window to keep into account the AO2D/EvSel BC discrepancy by %d BCs", ao2dBC - maxBCId);
144-
maxBCId = ao2dBC;
145-
}
146-
147-
bcRanges.push_back(std::make_pair(minBCId, maxBCId));
83+
bcRanges.push_back(getIRFrame(collision));
84+
nSelected++;
14885
}
86+
nColl++;
14987
filt++;
15088
}
15189

152-
/// 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
153-
std::sort(bcRanges.begin(), bcRanges.end(), [](const std::pair<uint64_t, uint64_t>& a, const std::pair<uint64_t, uint64_t>& b) {
154-
return a.first < b.first;
155-
});
15690
if (bcRanges.empty()) {
157-
LOGF(warning, "No BCs selected! This should not happen! Adding a bogus BC range to avoid crashes.");
158-
bcRanges.push_back(std::make_pair(0, 0));
91+
LOGF(warning, "No BCs selected!");
92+
return;
93+
}
94+
95+
float fractionSelected{static_cast<float>(nSelected) / nColl};
96+
int nMB{std::min(static_cast<int>(fillFac * nColl) - nSelected, nColl - 1)};
97+
LOGF(info, "Selected %d collisions (%.2f%%) and %d MB events", nSelected, fractionSelected * 100, nMB);
98+
int maxCollisionId = std::max(nMB, firstSelectedCollision);
99+
int minCollisionId = (maxCollisionId == nMB) ? 0 : firstSelectedCollision - nMB;
100+
auto minCollision = cols.begin() + minCollisionId;
101+
IRFrame minFrame{getIRFrame(minCollision)};
102+
bcRanges[0].getMin() = std::min(bcRanges[0].getMin(), minFrame.getMin());
103+
if (maxCollisionId == nMB) {
104+
auto maxCollision = cols.begin() + nMB;
105+
IRFrame maxFrame{getIRFrame(maxCollision)};
106+
bcRanges[0].getMax() = std::max(bcRanges[0].getMax(), maxFrame.getMax());
159107
}
160-
std::vector<std::pair<uint64_t, uint64_t>> bcRangesMerged(1, bcRanges[0]);
108+
109+
/// 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
110+
std::sort(bcRanges.begin(), bcRanges.end(), [](const IRFrame& a, const IRFrame& b) {
111+
return a.getMin() < b.getMin();
112+
});
113+
114+
std::vector<IRFrame> bcRangesMerged(1, bcRanges[0]);
161115
for (uint64_t iR{1}; iR < bcRanges.size(); ++iR) {
162-
if (bcRangesMerged.back().second >= bcRanges[iR].first) {
163-
bcRangesMerged.back().second = std::max(bcRangesMerged.back().second, bcRanges[iR].second);
116+
if (bcRangesMerged.back().getMax() >= bcRanges[iR].getMin()) {
117+
bcRangesMerged.back().getMax() = std::max(bcRangesMerged.back().getMax(), bcRanges[iR].getMax());
164118
} else {
165119
bcRangesMerged.push_back(bcRanges[iR]);
166120
}
167121
}
168122
bcRanges.swap(bcRangesMerged);
169123

170-
// 2. extend ranges
171-
int nBCselected{0};
172124
for (auto& range : bcRanges) {
173-
nBCselected += range.second - range.first + 1;
125+
tags(range.getMin().toLong(), range.getMax().toLong());
174126
}
175-
int nToBeAdded = std::ceil((bcs.size() - nBCselected) * fillFac);
176-
int nToBeAddedPerRange = std::ceil(float(nToBeAdded) / bcRanges.size());
177-
LOGF(debug, "Extending ranges by %d BCs (%d per selected range)", nToBeAdded, nToBeAddedPerRange);
178-
179-
InteractionRecord IR1, IR2;
180-
uint64_t first{bcs.iteratorAt(bcRanges[0].first).globalBC()};
181-
IR1.setFromLong(first);
182-
while (nToBeAdded > 0 && bcRanges[0].first > 0) { /// TODO: decide if we want to extend the ranges to the beginning of the dataframe
183-
first = bcs.iteratorAt(bcRanges[0].first - 1).globalBC();
184-
IR2.setFromLong(first);
185-
if (IR1.differenceInBC(IR2) > o2::constants::lhc::LHCMaxBunches) { // protection against change of orbit in the DataFrame
186-
LOGF(debug, "Jump by more than one orbit detected!");
187-
break;
188-
} else {
189-
IR1.setFromLong(first);
190-
}
191-
bcRanges[0].first--;
192-
nToBeAdded--;
193-
}
194-
LOGF(debug, "Extending ranges by %d BCs (%d per selected range)", nToBeAdded, nToBeAddedPerRange);
195-
196-
for (uint64_t iR{0}; iR < bcRanges.size() && nToBeAdded > 0; ++iR) {
197-
uint64_t second{bcs.iteratorAt(bcRanges[iR].second).globalBC()};
198-
IR2.setFromLong(second);
199-
for (int i{0}; i < nToBeAddedPerRange && nToBeAdded > 0; ++i) {
200-
if (bcRanges[iR].second < bcs.size() - 1) {
201-
second = bcs.iteratorAt(bcRanges[iR].second + 1).globalBC();
202-
IR1.setFromLong(second);
203-
if (IR1.differenceInBC(IR2) > o2::constants::lhc::LHCMaxBunches) { // protection against change of orbit in the DataFrame
204-
LOGF(debug, "Jump by more than one orbit detected!");
205-
break;
206-
} else {
207-
IR2.setFromLong(second);
208-
}
209-
bcRanges[iR].second++;
210-
nToBeAdded--;
211-
}
212-
}
213-
}
214-
LOGF(debug, "End extension, remaining to be added %d BCs", nToBeAdded);
215-
216-
// fill res
217-
LOGF(debug, "Merged and extended sorted ranges");
218-
for (auto& range : bcRanges) {
219-
LOGF(debug, " %i - %i", range.first, range.second);
220-
uint64_t first{bcs.rawIteratorAt(range.first).globalBC()}, second{bcs.rawIteratorAt(range.second).globalBC()};
221-
IR1.setFromLong(first);
222-
IR2.setFromLong(second);
223-
res.emplace_back(IR1, IR2);
224-
tags(first, second);
225-
}
226-
227-
// make res an output
228-
pc.outputs().snapshot({"PPF", "IFRAMES", 0, Lifetime::Timeframe}, res);
229-
230-
// clean up
231-
res.clear();
232127
}
233128

234-
// need a trivial process method
235-
// the parameters determine the tables available in the input
236-
void process(CCs const& collisions, BCs const& bcs, FDs const& fdecisions)
129+
// need a trivial process method: the parameters determine the tables available in the input
130+
void process(CCs const& collisions, aod::BCs const& bcs, aod::CefpDecisions const& fdecisions)
237131
{
238132
}
239133
};
240134

241135
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
242136
{
243137
DataProcessorSpec spec{adaptAnalysisTask<BCRangeSelector>(cfgc, TaskName{"bc-ranges-selector-task"})};
244-
245-
// add output
246-
spec.outputs.emplace_back("PPF", "IFRAMES", 0, Lifetime::Timeframe);
247-
LOGF(debug, "Output %i", spec.outputs.size());
248-
249138
return WorkflowSpec{spec};
250139
}

0 commit comments

Comments
 (0)