Project
Loading...
Searching...
No Matches
BuildTopologyDictionary.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
13
17
20
21#include "TFile.h"
22
24
25namespace o2::its3
26{
27void BuildTopologyDictionary::accountTopology(const itsmft::ClusterTopology& cluster, bool IB, float dX, float dZ)
28{
29 accountTopologyImpl(cluster,
30 ((IB) ? mMapInfoIB : mMapInfoOB),
31 ((IB) ? mTopologyMapIB : mTopologyMapOB),
32 ((IB) ? mTotClustersIB : mTotClustersOB),
35 dX, dZ);
36}
37
38void BuildTopologyDictionary::accountTopologyImpl(const itsmft::ClusterTopology& cluster, TopoInfo& tinfo, TopoStat& tstat, unsigned int& tot, float sigmaX, float sigmaZ, float dX, float dZ)
39{
40 ++tot;
41 bool useDf = dX < IgnoreVal / 2; // we may need to account the frequency but to not update the centroid
42
43 auto& topoStat = tstat[cluster.getHash()];
44 topoStat.countsTotal++;
45 if (topoStat.countsTotal == 1) { // a new topology is inserted
46 topoStat.topology = cluster;
47 //___________________DEFINING_TOPOLOGY_CHARACTERISTICS__________________
49 topInf.mPattern.setPattern(cluster.getPattern().data());
50 topInf.mSizeX = cluster.getRowSpan();
51 topInf.mSizeZ = cluster.getColumnSpan();
52 //__________________COG_Determination_____________
53 topInf.mNpixels = cluster.getClusterPattern().getCOG(topInf.mCOGx, topInf.mCOGz);
54 if (useDf) {
55 topInf.mXmean = dX;
56 topInf.mZmean = dZ;
57 topoStat.countsWithBias = 1;
58 } else { // assign expected sigmas from the pixel X, Z sizes
59 topInf.mXsigma2 = sigmaX * sigmaX / 12.f / (float)std::min(10, topInf.mSizeX);
60 topInf.mZsigma2 = sigmaZ * sigmaZ / 12.f / (float)std::min(10, topInf.mSizeZ);
61 }
62 tinfo.emplace(cluster.getHash(), topInf);
63 } else {
64 if (useDf) {
65 auto num = topoStat.countsWithBias++;
66 auto ind = tinfo.find(cluster.getHash());
67 float tmpxMean = ind->second.mXmean;
68 float newxMean = ind->second.mXmean = ((tmpxMean)*num + dX) / (num + 1);
69 float tmpxSigma2 = ind->second.mXsigma2;
70 ind->second.mXsigma2 = (num * tmpxSigma2 + (dX - tmpxMean) * (dX - newxMean)) / (num + 1); // online variance algorithm
71 float tmpzMean = ind->second.mZmean;
72 float newzMean = ind->second.mZmean = ((tmpzMean)*num + dZ) / (num + 1);
73 float tmpzSigma2 = ind->second.mZsigma2;
74 ind->second.mZsigma2 = (num * tmpzSigma2 + (dZ - tmpzMean) * (dZ - newzMean)) / (num + 1); // online variance algorithm
75 }
76 }
77}
78
79void BuildTopologyDictionary::setNCommon(unsigned int nCommon, bool IB)
80{
81 mDictionary.resetMaps(IB);
82
83 auto& freqTopo = ((IB) ? mTopologyFrequencyIB : mTopologyFrequencyOB);
84 auto& freqThres = ((IB) ? mFrequencyThresholdIB : mFrequencyThresholdOB);
85 auto& comTopo = ((IB) ? mNCommonTopologiesIB : mNCommonTopologiesOB);
86 auto ntot = ((IB) ? mTotClustersIB : mTotClustersOB);
87
88 setNCommonImpl(nCommon,
89 freqTopo,
90 ((IB) ? mTopologyMapIB : mTopologyMapOB),
91 comTopo,
92 ntot);
93 // Recaculate also the threshold
94 freqThres = ((double)freqTopo[comTopo - 1].first) / ntot;
95}
96
97void BuildTopologyDictionary::setNCommonImpl(unsigned int ncom, TopoFreq& tfreq, TopoStat& tstat, unsigned int& ncommon, unsigned int ntot)
98{
100 LOGP(warning, "Redefining nCommon from {} to {} to be below InvalidPatternID", ncom, itsmft::CompCluster::InvalidPatternID - 1);
102 }
103 tfreq.clear();
104 for (auto&& p : tstat) { // p os pair<ulong,TopoStat>
105 tfreq.emplace_back(p.second.countsTotal, p.first);
106 }
107 std::sort(tfreq.begin(), tfreq.end(),
108 [](const std::pair<unsigned long, unsigned long>& couple1,
109 const std::pair<unsigned long, unsigned long>& couple2) { return (couple1.first > couple2.first); });
110 ncommon = ncom;
111}
112
114{
115 mDictionary.resetMaps(IB);
116 setThresholdImpl(thr,
117 ((IB) ? mTopologyFrequencyIB : mTopologyFrequencyOB),
118 ((IB) ? mMapInfoIB : mMapInfoOB),
119 ((IB) ? mTopologyMapIB : mTopologyMapOB),
120 ((IB) ? mNCommonTopologiesIB : mNCommonTopologiesOB),
121 ((IB) ? mFrequencyThresholdIB : mFrequencyThresholdOB),
122 ((IB) ? mTotClustersIB : mTotClustersOB));
123}
124
125void BuildTopologyDictionary::setThresholdImpl(double thr, TopoFreq& tfreq, TopoInfo& tinfo, TopoStat& tstat, unsigned int& ncommon, double& freqthres, unsigned int ntot)
126{
127 setNCommonImpl(0, tfreq, tstat, ncommon, ntot);
128 freqthres = thr;
129 for (const auto& q : tfreq) {
130 if (((double)q.first) / ntot > thr) {
131 ++ncommon;
132 } else {
133 break;
134 }
135 }
137 freqthres = ((double)tfreq[itsmft::CompCluster::InvalidPatternID - 1].first) / ntot;
138 LOGP(warning, "Redefining prob. threshold from {} to {} to be below InvalidPatternID (was {})", thr, freqthres, ntot);
140 }
141}
142
144{
145 if (cumulative <= 0. || cumulative >= 1.) {
146 cumulative = 0.99;
147 }
148
149 auto& freqTopo = ((IB) ? mTopologyFrequencyIB : mTopologyFrequencyOB);
150 auto& freqThres = ((IB) ? mFrequencyThresholdIB : mFrequencyThresholdOB);
151 auto& statTopo = ((IB) ? mTopologyMapIB : mTopologyMapOB);
152 auto& comTopo = ((IB) ? mNCommonTopologiesIB : mNCommonTopologiesOB);
153 auto ntot = ((IB) ? mTotClustersIB : mTotClustersOB);
154
155 mDictionary.resetMaps(IB);
156 setNCommonImpl(0, freqTopo, statTopo, comTopo, ntot);
157 setThresholdCumulativeImpl(cumulative, freqTopo, comTopo, freqThres, ntot);
158}
159
160void BuildTopologyDictionary::setThresholdCumulativeImpl(double cumulative, TopoFreq& tfreq, unsigned int& ncommon, double& freqthres, unsigned int ntot)
161{
162 double totFreq = 0.;
163 for (auto& q : tfreq) {
164 totFreq += ((double)(q.first)) / ntot;
165 if (totFreq < cumulative) {
166 ++ncommon;
168 totFreq -= ((double)(q.first)) / ntot;
169 --ncommon;
170 LOGP(warning, "Redefining cumulative threshould from {} to {} to be below InvalidPatternID)", cumulative, totFreq);
171 }
172 } else {
173 break;
174 }
175 }
176 freqthres = ((double)(tfreq[--ncommon].first)) / ntot;
177 while (std::fabs(((double)tfreq[ncommon--].first) / ntot - freqthres) < 1.e-15) {
178 }
179 freqthres = ((double)tfreq[ncommon++].first) / ntot;
180}
181
183{
184 LOG(info) << "Dictionary finalisation";
185 LOG(info) << "Number of IB clusters: " << mTotClustersIB;
186 LOG(info) << "Number of OB clusters: " << mTotClustersOB;
187
188 groupRareTopologiesImpl(mTopologyFrequencyIB, mMapInfoIB, mTopologyMapIB, mNCommonTopologiesIB, mFrequencyThresholdIB, mDictionary.mDataIB, mNCommonTopologiesIB);
189 groupRareTopologiesImpl(mTopologyFrequencyOB, mMapInfoOB, mTopologyMapOB, mNCommonTopologiesOB, mFrequencyThresholdOB, mDictionary.mDataOB, mNCommonTopologiesOB);
190
191 LOG(info) << "Dictionay finalised";
192 LOG(info) << "IB:";
193 mDictionary.mDataIB.print();
194 LOG(info) << "OB:";
195 mDictionary.mDataOB.print();
196}
197
198void BuildTopologyDictionary::groupRareTopologiesImpl(TopoFreq& tfreq, TopoInfo& tinfo, TopoStat& tstat, unsigned int& ncommon, double& freqthres, TopologyDictionaryData& data, unsigned int ntot)
199{
200 double totFreq = 0.;
201 for (unsigned int j = 0; j < ncommon; j++) {
203 gr.mHash = tfreq[j].second;
204 gr.mFrequency = ((double)(tfreq[j].first)) / ntot;
205 totFreq += gr.mFrequency;
206 // rough estimation for the error considering a uniform distribution
207 const auto& topo = tinfo.find(gr.mHash)->second;
208 gr.mErrX = std::sqrt(topo.mXsigma2);
209 gr.mErrZ = std::sqrt(topo.mZsigma2);
210 gr.mErr2X = topo.mXsigma2;
211 gr.mErr2Z = topo.mZsigma2;
212 gr.mXCOG = -1 * topo.mCOGx;
213 gr.mZCOG = topo.mCOGz;
214 gr.mNpixels = topo.mNpixels;
215 gr.mPattern = topo.mPattern;
216 gr.mIsGroup = false;
217 data.mVectorOfIDs.push_back(gr);
219 LOGP(warning, "Limiting N unique topologies to {}, threshold freq. to {}, cumulative freq. to {} to be below InvalidPatternID", j, gr.mFrequency, totFreq);
220 ncommon = j;
221 freqthres = gr.mFrequency;
222 break;
223 }
224 }
225 // groupRareTopologies based on binning over number of rows and columns (TopologyDictionary::NumberOfRowClasses *
226 // NumberOfColClasses)
227
228 std::unordered_map<int, std::pair<itsmft::GroupStruct, unsigned long>> tmp_GroupMap; //<group ID, <Group struct, counts>>
229
230 int grNum = 0;
231 int rowBinEdge = 0;
232 int colBinEdge = 0;
233 for (int iRowClass = 0; iRowClass < its3::TopologyDictionary::MaxNumberOfRowClasses; iRowClass++) {
234 for (int iColClass = 0; iColClass < its3::TopologyDictionary::MaxNumberOfColClasses; iColClass++) {
235 rowBinEdge = (iRowClass + 1) * its3::TopologyDictionary::RowClassSpan;
236 colBinEdge = (iColClass + 1) * its3::TopologyDictionary::ColClassSpan;
237 grNum = its3::LookUp::groupFinder(rowBinEdge, colBinEdge);
238 // Create a structure for a group of rare topologies
239 itsmft::GroupStruct gr;
240 gr.mHash = (((unsigned long)(grNum)) << 32) & 0xffffffff00000000;
241 gr.mErrX = its3::TopologyDictionary::RowClassSpan / std::sqrt(12.f * (float)std::min(10, rowBinEdge));
242 gr.mErrZ = its3::TopologyDictionary::ColClassSpan / std::sqrt(12.f * (float)std::min(10, colBinEdge));
243 gr.mErr2X = gr.mErrX * gr.mErrX;
244 gr.mErr2Z = gr.mErrZ * gr.mErrZ;
245 gr.mXCOG = 0;
246 gr.mZCOG = 0;
247 gr.mNpixels = rowBinEdge * colBinEdge;
248 gr.mIsGroup = true;
249 gr.mFrequency = 0.;
251 unsigned char dummyPattern[itsmft::ClusterPattern::kExtendedPatternBytes] = {0};
252 dummyPattern[0] = (unsigned char)rowBinEdge;
253 dummyPattern[1] = (unsigned char)colBinEdge;
254 int nBits = rowBinEdge * colBinEdge;
255 int nBytes = nBits / 8;
256 for (int iB = 2; iB < nBytes + 2; iB++) {
257 dummyPattern[iB] = (unsigned char)255;
258 }
259 int residualBits = nBits % 8;
260 if (residualBits != 0) {
261 unsigned char tempChar = 0;
262 while (residualBits > 0) {
263 residualBits--;
264 tempChar |= 1 << (7 - residualBits);
265 }
266 dummyPattern[nBytes + 2] = tempChar;
267 }
268 gr.mPattern.setPattern(dummyPattern);
269 // Filling the map for groups
270 tmp_GroupMap[grNum] = std::make_pair(gr, 0);
271 }
272 }
273 int rs{}, cs{}, index{};
274
275 // Updating the counts for the groups of rare topologies
276 for (auto j{ncommon}; j < tfreq.size(); j++) {
277 unsigned long hash1 = tfreq[j].second;
278 rs = tstat.find(hash1)->second.topology.getRowSpan();
279 cs = tstat.find(hash1)->second.topology.getColumnSpan();
281 tmp_GroupMap[index].second += tfreq[j].first;
282 }
283
284 for (auto&& p : tmp_GroupMap) {
285 itsmft::GroupStruct& group = p.second.first;
286 group.mFrequency = ((double)p.second.second) / ntot;
287 data.mVectorOfIDs.push_back(group);
288 }
289
290 // Sorting the dictionary preserving all unique topologies
291 std::sort(data.mVectorOfIDs.begin(), data.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) {
292 return (!a.mIsGroup) && b.mIsGroup ? true : (a.mIsGroup && (!b.mIsGroup) ? false : (a.mFrequency > b.mFrequency));
293 });
294 if (data.mVectorOfIDs.size() >= itsmft::CompCluster::InvalidPatternID - 1) {
295 LOGP(warning, "Max allowed {} patterns is reached, stopping", itsmft::CompCluster::InvalidPatternID - 1);
296 data.mVectorOfIDs.resize(itsmft::CompCluster::InvalidPatternID - 1);
297 }
298 // Sorting the dictionary to final form
299 std::sort(data.mVectorOfIDs.begin(), data.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) { return a.mFrequency > b.mFrequency; });
300 // Creating the map for common topologies
301 for (int iKey = 0; iKey < data.mVectorOfIDs.size(); iKey++) {
302 itsmft::GroupStruct& gr = data.mVectorOfIDs[iKey];
303 if (!gr.mIsGroup) {
304 data.mCommonMap.emplace(gr.mHash, iKey);
305 if (gr.mPattern.getUsedBytes() == 1) {
306 data.mSmallTopologiesLUT[(gr.mPattern.getColumnSpan() - 1) * 255 + (int)gr.mPattern.getByte(2)] = iKey;
307 }
308 } else {
309 data.mGroupMap.emplace((int)(gr.mHash >> 32) & 0x00000000ffffffff, iKey);
310 }
311 }
312}
313
314std::ostream& operator<<(std::ostream& os, const BuildTopologyDictionary& DB)
315{
316 os << "--- InnerBarrel\n";
317 for (unsigned int i = 0; i < DB.mNCommonTopologiesIB; i++) {
318 const unsigned long& hash = DB.mTopologyFrequencyIB[i].second;
319 os << "Hash: " << hash << '\n';
320 os << "counts: " << DB.mTopologyMapIB.find(hash)->second.countsTotal;
321 os << " (with bias provided: " << DB.mTopologyMapIB.find(hash)->second.countsWithBias << ")" << '\n';
322 os << "sigmaX: " << std::sqrt(DB.mMapInfoIB.find(hash)->second.mXsigma2) << '\n';
323 os << "sigmaZ: " << std::sqrt(DB.mMapInfoIB.find(hash)->second.mZsigma2) << '\n';
324 os << DB.mTopologyMapIB.find(hash)->second.topology;
325 }
326 os << "--- OuterBarrel\n";
327 for (unsigned int i = 0; i < DB.mNCommonTopologiesOB; i++) {
328 const unsigned long& hash = DB.mTopologyFrequencyOB[i].second;
329 os << "Hash: " << hash << '\n';
330 os << "counts: " << DB.mTopologyMapOB.find(hash)->second.countsTotal;
331 os << " (with bias provided: " << DB.mTopologyMapOB.find(hash)->second.countsWithBias << ")" << '\n';
332 os << "sigmaX: " << std::sqrt(DB.mMapInfoOB.find(hash)->second.mXsigma2) << '\n';
333 os << "sigmaZ: " << std::sqrt(DB.mMapInfoOB.find(hash)->second.mZsigma2) << '\n';
334 os << DB.mTopologyMapOB.find(hash)->second.topology;
335 }
336 return os;
337}
338
339void BuildTopologyDictionary::printDictionary(const std::string& fname)
340{
341 LOG(info) << "Saving the the dictionary in binary format: ";
342 std::ofstream out(fname);
343 out << mDictionary;
344 out.close();
345 LOG(info) << " `-> done!";
346}
347
349{
350 LOG(info) << "Printing the dictionary: ";
351 std::ofstream out(fname);
352 mDictionary.writeBinaryFile(fname);
353 out.close();
354 LOG(info) << " `-> done!";
355}
356
358{
359 LOG(info) << "Saving the the dictionary in a ROOT file: " << fname;
360 TFile output(fname.c_str(), "recreate");
361 output.WriteObjectAny(&mDictionary, mDictionary.Class(), "ccdb_object");
362 output.Close();
363 LOG(info) << " `-> done!";
364}
365
366} // namespace o2::its3
Definition of the ITSMFT compact cluster.
int32_t i
ClassImp(o2::itsmft::BuildTopologyDictionary)
void output(const std::map< std::string, ChannelStat > &channels)
Definition rawdump.cxx:197
uint32_t j
Definition RawData.h:0
Definition of the SegmentationAlpide class.
Definition of the SegmentationMosaix class.
Definition of the LookUp class for its3.
double num
void saveDictionaryRoot(const std::string &fname)
void setNCommon(unsigned int nCommon, bool IB)
void setThresholdCumulative(double cumulative, bool IB)
void accountTopology(const itsmft::ClusterTopology &cluster, bool IB, float dX=IgnoreVal, float dZ=IgnoreVal)
void printDictionary(const std::string &fname)
void printDictionaryBinary(const std::string &fname)
static int groupFinder(int nRow, int nCol)
Definition LookUp.cxx:47
static constexpr float PitchCol
static constexpr float PitchRow
static constexpr int ColClassSpan
Column span of the classes of rare topologies.
static constexpr int MaxNumberOfRowClasses
Maximum number of row classes for the groups of rare topologies.
static constexpr int MaxNumberOfColClasses
Maximum number of col classes for the groups of rare topologies.
static constexpr int RowClassSpan
Row span of the classes of rare topologies.
void writeBinaryFile(const std::string &outputFile)
Prints the dictionary in a binary file.
void resetMaps(bool IB=true) noexcept
static int getCOG(int rowSpan, int colSpan, const unsigned char patt[MaxPatternBytes], float &xCOG, float &zCOG)
Static: Compute pattern's COG position. Returns the number of fired pixels.
void setPattern(int nRow, int nCol, const unsigned char patt[MaxPatternBytes])
Sets the pattern.
static constexpr int kExtendedPatternBytes
Maximum number of bytes for the cluster puttern + 2 bytes respectively for the number of rows and col...
int getRowSpan() const
Returns the number of rows.
int getColumnSpan() const
Returns the number of columns.
const std::array< unsigned char, ClusterPattern::kExtendedPatternBytes > & getPattern() const
Returns the pattern.
unsigned long getHash() const
Returns the hashcode.
const ClusterPattern & getClusterPattern() const
static constexpr unsigned short InvalidPatternID
Definition CompCluster.h:46
static constexpr float PitchCol
static constexpr float PitchRow
GLuint index
Definition glcorearb.h:781
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLboolean * data
Definition glcorearb.h:298
GLboolean GLuint group
Definition glcorearb.h:3991
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
std::ostream & operator<<(std::ostream &os, const BuildTopologyDictionary &DB)
Structure containing the most relevant pieces of information of a topology.
float mErrX
Error associated to the hit point in the x direction.
float mErr2X
Squared Error associated to the hit point in the x direction.
float mErrZ
Error associated to the hit point in the z direction.
unsigned long mHash
Hashcode.
float mXCOG
x position of the COG wrt the bottom left corner of the bounding box
int mNpixels
Number of fired pixels.
double mFrequency
Frequency of the topology.
bool mIsGroup
false: common topology; true: group of rare topologies
float mErr2Z
Squared Error associated to the hit point in the z direction.
float mZCOG
z position of the COG wrt the bottom left corner of the bounding box
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"