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
18
19#include "TFile.h"
20
22
23namespace o2::its3
24{
26{
27 mTotClusters++;
28 bool useDf = dX < IgnoreVal / 2; // we may need to account the frequency but to not update the centroid
29
30 // std::pair<unordered_map<unsigned long, itsmft::TopoStat>::iterator,bool> ret;
31 // auto ret = mTopologyMap.insert(std::make_pair(cluster.getHash(), std::make_pair(cluster, 1)));
32 auto& topoStat = mTopologyMap[cluster.getHash()];
33 topoStat.countsTotal++;
34 if (topoStat.countsTotal == 1) { // a new topology is inserted
35 topoStat.topology = cluster;
36 //___________________DEFINING_TOPOLOGY_CHARACTERISTICS__________________
38 topInf.mPattern.setPattern(cluster.getPattern().data());
39 topInf.mSizeX = cluster.getRowSpan();
40 topInf.mSizeZ = cluster.getColumnSpan();
41 //__________________COG_Determination_____________
42 topInf.mNpixels = cluster.getClusterPattern().getCOG(topInf.mCOGx, topInf.mCOGz);
43 if (useDf) {
44 topInf.mXmean = dX;
45 topInf.mZmean = dZ;
46 topoStat.countsWithBias = 1;
47 } else { // assign expected sigmas from the pixel X, Z sizes
48 topInf.mXsigma2 = 1.f / 12.f / (float)std::min(10, topInf.mSizeX);
49 topInf.mZsigma2 = 1.f / 12.f / (float)std::min(10, topInf.mSizeZ);
50 }
51 mMapInfo.emplace(cluster.getHash(), topInf);
52 } else {
53 if (useDf) {
54 auto num = topoStat.countsWithBias++;
55 auto ind = mMapInfo.find(cluster.getHash());
56 float tmpxMean = ind->second.mXmean;
57 float newxMean = ind->second.mXmean = ((tmpxMean)*num + dX) / (num + 1);
58 float tmpxSigma2 = ind->second.mXsigma2;
59 ind->second.mXsigma2 = (num * tmpxSigma2 + (dX - tmpxMean) * (dX - newxMean)) / (num + 1); // online variance algorithm
60 float tmpzMean = ind->second.mZmean;
61 float newzMean = ind->second.mZmean = ((tmpzMean)*num + dZ) / (num + 1);
62 float tmpzSigma2 = ind->second.mZsigma2;
63 ind->second.mZsigma2 = (num * tmpzSigma2 + (dZ - tmpzMean) * (dZ - newzMean)) / (num + 1); // online variance algorithm
64 }
65 }
66}
67
69{
70 mTopologyFrequency.clear();
71 for (auto&& p : mTopologyMap) { // p is pair<ulong,TopoStat>
72 mTopologyFrequency.emplace_back(p.second.countsTotal, p.first);
73 }
74 std::sort(mTopologyFrequency.begin(), mTopologyFrequency.end(),
75 [](const std::pair<unsigned long, unsigned long>& couple1,
76 const std::pair<unsigned long, unsigned long>& couple2) { return (couple1.first > couple2.first); });
77 mNCommonTopologies = 0;
78 mDictionary.mCommonMap.clear();
79 mDictionary.mGroupMap.clear();
80 mFrequencyThreshold = thr;
81 for (auto& q : mTopologyFrequency) {
82 if (((double)q.first) / mTotClusters > thr) {
83 mNCommonTopologies++;
84 } else {
85 break;
86 }
87 }
88 if (mNCommonTopologies >= itsmft::CompCluster::InvalidPatternID) {
89 mFrequencyThreshold = ((double)mTopologyFrequency[itsmft::CompCluster::InvalidPatternID - 1].first) / mTotClusters;
90 LOGP(warning, "Redefining prob. threshould from {} to {} to be below InvalidPatternID (was {})", thr, mFrequencyThreshold, mNCommonTopologies);
91 mNCommonTopologies = itsmft::CompCluster::InvalidPatternID - 1;
92 }
93}
94
95void BuildTopologyDictionary::setNCommon(unsigned int nCommon)
96{
98 LOGP(warning, "Redefining nCommon from {} to {} to be below InvalidPatternID", nCommon, itsmft::CompCluster::InvalidPatternID - 1);
100 }
101 mTopologyFrequency.clear();
102 for (auto&& p : mTopologyMap) { // p os pair<ulong,TopoStat>
103 mTopologyFrequency.emplace_back(p.second.countsTotal, p.first);
104 }
105 std::sort(mTopologyFrequency.begin(), mTopologyFrequency.end(),
106 [](const std::pair<unsigned long, unsigned long>& couple1,
107 const std::pair<unsigned long, unsigned long>& couple2) { return (couple1.first > couple2.first); });
108 mNCommonTopologies = nCommon;
109 mDictionary.mCommonMap.clear();
110 mDictionary.mGroupMap.clear();
111 mFrequencyThreshold = ((double)mTopologyFrequency[mNCommonTopologies - 1].first) / mTotClusters;
112}
113
115{
116 mTopologyFrequency.clear();
117 if (cumulative <= 0. || cumulative >= 1.) {
118 cumulative = 0.99;
119 }
120 double totFreq = 0.;
121 for (auto&& p : mTopologyMap) { // p os pair<ulong,TopoStat>
122 mTopologyFrequency.emplace_back(p.second.countsTotal, p.first);
123 }
124 std::sort(mTopologyFrequency.begin(), mTopologyFrequency.end(),
125 [](const std::pair<unsigned long, unsigned long>& couple1,
126 const std::pair<unsigned long, unsigned long>& couple2) { return (couple1.first > couple2.first); });
127 mNCommonTopologies = 0;
128 mDictionary.mCommonMap.clear();
129 mDictionary.mGroupMap.clear();
130 for (auto& q : mTopologyFrequency) {
131 totFreq += ((double)(q.first)) / mTotClusters;
132 if (totFreq < cumulative) {
133 mNCommonTopologies++;
134 if (mNCommonTopologies >= itsmft::CompCluster::InvalidPatternID) {
135 totFreq -= ((double)(q.first)) / mTotClusters;
136 mNCommonTopologies--;
137 LOGP(warning, "Redefining cumulative threshould from {} to {} to be below InvalidPatternID)", cumulative, totFreq);
138 }
139 } else {
140 break;
141 }
142 }
143 mFrequencyThreshold = ((double)(mTopologyFrequency[--mNCommonTopologies].first)) / mTotClusters;
144 while (std::fabs(((double)mTopologyFrequency[mNCommonTopologies].first) / mTotClusters - mFrequencyThreshold) < 1.e-15) {
145 mNCommonTopologies--;
146 }
147 mFrequencyThreshold = ((double)mTopologyFrequency[mNCommonTopologies++].first) / mTotClusters;
148}
149
151{
152 LOG(info) << "Dictionary finalisation";
153 LOG(info) << "Number of clusters: " << mTotClusters;
154
155 double totFreq = 0.;
156 for (unsigned int j = 0; j < mNCommonTopologies; j++) {
158 gr.mHash = mTopologyFrequency[j].second;
159 gr.mFrequency = ((double)(mTopologyFrequency[j].first)) / mTotClusters;
160 totFreq += gr.mFrequency;
161 // rough estimation for the error considering a8 uniform distribution
162 const auto& topo = mMapInfo.find(gr.mHash)->second;
163 gr.mErrX = std::sqrt(topo.mXsigma2);
164 gr.mErrZ = std::sqrt(topo.mZsigma2);
165 gr.mErr2X = topo.mXsigma2;
166 gr.mErr2Z = topo.mZsigma2;
167 gr.mXCOG = -1 * topo.mCOGx;
168 gr.mZCOG = topo.mCOGz;
169 gr.mNpixels = topo.mNpixels;
170 gr.mPattern = topo.mPattern;
171 gr.mIsGroup = false;
172 mDictionary.mVectorOfIDs.push_back(gr);
174 LOGP(warning, "Limiting N unique topologies to {}, threshold freq. to {}, cumulative freq. to {} to be below InvalidPatternID", j, gr.mFrequency, totFreq);
175 mNCommonTopologies = j;
176 mFrequencyThreshold = gr.mFrequency;
177 break;
178 }
179 }
180 // groupRareTopologies based on binning over number of rows and columns (TopologyDictionary::NumberOfRowClasses *
181 // NumberOfColClasses)
182
183 std::unordered_map<int, std::pair<itsmft::GroupStruct, unsigned long>> tmp_GroupMap; //<group ID, <Group struct, counts>>
184
185 int grNum = 0;
186 int rowBinEdge = 0;
187 int colBinEdge = 0;
188 for (int iRowClass = 0; iRowClass < its3::TopologyDictionary::MaxNumberOfRowClasses; iRowClass++) {
189 for (int iColClass = 0; iColClass < its3::TopologyDictionary::MaxNumberOfColClasses; iColClass++) {
190 rowBinEdge = (iRowClass + 1) * its3::TopologyDictionary::RowClassSpan;
191 colBinEdge = (iColClass + 1) * its3::TopologyDictionary::ColClassSpan;
192 grNum = its3::LookUp::groupFinder(rowBinEdge, colBinEdge);
193 // Create a structure for a group of rare topologies
195 gr.mHash = (((unsigned long)(grNum)) << 32) & 0xffffffff00000000;
196 gr.mErrX = its3::TopologyDictionary::RowClassSpan / std::sqrt(12 * std::min(10, rowBinEdge));
197 gr.mErrZ = its3::TopologyDictionary::ColClassSpan / std::sqrt(12 * std::min(10, colBinEdge));
198 gr.mErr2X = gr.mErrX * gr.mErrX;
199 gr.mErr2Z = gr.mErrZ * gr.mErrZ;
200 gr.mXCOG = 0;
201 gr.mZCOG = 0;
202 gr.mNpixels = rowBinEdge * colBinEdge;
203 gr.mIsGroup = true;
204 gr.mFrequency = 0.;
206 unsigned char dummyPattern[itsmft::ClusterPattern::kExtendedPatternBytes] = {0};
207 dummyPattern[0] = (unsigned char)rowBinEdge;
208 dummyPattern[1] = (unsigned char)colBinEdge;
209 int nBits = rowBinEdge * colBinEdge;
210 int nBytes = nBits / 8;
211 for (int iB = 2; iB < nBytes + 2; iB++) {
212 dummyPattern[iB] = (unsigned char)255;
213 }
214 int residualBits = nBits % 8;
215 if (residualBits != 0) {
216 unsigned char tempChar = 0;
217 while (residualBits > 0) {
218 residualBits--;
219 tempChar |= 1 << (7 - residualBits);
220 }
221 dummyPattern[nBytes + 2] = tempChar;
222 }
223 gr.mPattern.setPattern(dummyPattern);
224 // Filling the map for groups
225 tmp_GroupMap[grNum] = std::make_pair(gr, 0);
226 }
227 }
228 int rs{}, cs{}, index{};
229
230 // Updating the counts for the groups of rare topologies
231 for (auto j{mNCommonTopologies}; j < mTopologyFrequency.size(); j++) {
232 unsigned long hash1 = mTopologyFrequency[j].second;
233 rs = mTopologyMap.find(hash1)->second.topology.getRowSpan();
234 cs = mTopologyMap.find(hash1)->second.topology.getColumnSpan();
236 tmp_GroupMap[index].second += mTopologyFrequency[j].first;
237 }
238
239 for (auto&& p : tmp_GroupMap) {
240 itsmft::GroupStruct& group = p.second.first;
241 group.mFrequency = ((double)p.second.second) / mTotClusters;
242 mDictionary.mVectorOfIDs.push_back(group);
243 }
244
245 // Sorting the dictionary preserving all unique topologies
246 std::sort(mDictionary.mVectorOfIDs.begin(), mDictionary.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) {
247 return (!a.mIsGroup) && b.mIsGroup ? true : (a.mIsGroup && (!b.mIsGroup) ? false : (a.mFrequency > b.mFrequency));
248 });
249 if (mDictionary.mVectorOfIDs.size() >= itsmft::CompCluster::InvalidPatternID - 1) {
250 LOGP(warning, "Max allowed {} patterns is reached, stopping", itsmft::CompCluster::InvalidPatternID - 1);
251 mDictionary.mVectorOfIDs.resize(itsmft::CompCluster::InvalidPatternID - 1);
252 }
253 // Sorting the dictionary to final form
254 std::sort(mDictionary.mVectorOfIDs.begin(), mDictionary.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) { return a.mFrequency > b.mFrequency; });
255 // Creating the map for common topologies
256 for (int iKey = 0; iKey < mDictionary.getSize(); iKey++) {
257 itsmft::GroupStruct& gr = mDictionary.mVectorOfIDs[iKey];
258 if (!gr.mIsGroup) {
259 mDictionary.mCommonMap.emplace(gr.mHash, iKey);
260 if (gr.mPattern.getUsedBytes() == 1) {
261 mDictionary.mSmallTopologiesLUT[(gr.mPattern.getColumnSpan() - 1) * 255 + (int)gr.mPattern.getByte(2)] = iKey;
262 }
263 } else {
264 mDictionary.mGroupMap.emplace((int)(gr.mHash >> 32) & 0x00000000ffffffff, iKey);
265 }
266 }
267 LOG(info) << "Dictionay finalised";
268 LOG(info) << "Number of keys: " << mDictionary.getSize();
269 LOG(info) << "Number of common topologies: " << mDictionary.mCommonMap.size();
270 LOG(info) << "Number of groups of rare topologies: " << mDictionary.mGroupMap.size();
271}
272
273std::ostream& operator<<(std::ostream& os, const BuildTopologyDictionary& DB)
274{
275 for (unsigned int i = 0; i < DB.mNCommonTopologies; i++) {
276 const unsigned long& hash = DB.mTopologyFrequency[i].second;
277 os << "Hash: " << hash << '\n';
278 os << "counts: " << DB.mTopologyMap.find(hash)->second.countsTotal;
279 os << " (with bias provided: " << DB.mTopologyMap.find(hash)->second.countsWithBias << ")" << '\n';
280 os << "sigmaX: " << std::sqrt(DB.mMapInfo.find(hash)->second.mXsigma2) << '\n';
281 os << "sigmaZ: " << std::sqrt(DB.mMapInfo.find(hash)->second.mZsigma2) << '\n';
282 os << DB.mTopologyMap.find(hash)->second.topology;
283 }
284 return os;
285}
286
287void BuildTopologyDictionary::printDictionary(const std::string& fname)
288{
289 LOG(info) << "Saving the the dictionary in binary format: ";
290 std::ofstream out(fname);
291 out << mDictionary;
292 out.close();
293 LOG(info) << " `-> done!";
294}
295
297{
298 LOG(info) << "Printing the dictionary: ";
299 std::ofstream out(fname);
300 mDictionary.writeBinaryFile(fname);
301 out.close();
302 LOG(info) << " `-> done!";
303}
304
306{
307 LOG(info) << "Saving the the dictionary in a ROOT file: " << fname;
308 TFile output(fname.c_str(), "recreate");
309 output.WriteObjectAny(&mDictionary, mDictionary.Class(), "ccdb_object");
310 output.Close();
311 LOG(info) << " `-> done!";
312}
313
314} // 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 SegmentationSuperAlpide class.
Definition of the LookUp class for its3.
double num
void saveDictionaryRoot(const std::string &fname)
void printDictionary(const std::string &fname)
void accountTopology(const itsmft::ClusterTopology &cluster, float dX=IgnoreVal, float dZ=IgnoreVal)
void printDictionaryBinary(const std::string &fname)
static int groupFinder(int nRow, int nCol)
Definition LookUp.cxx:45
static constexpr int ColClassSpan
Column span of the classes of rare topologies.
int getSize() const
Returns the number of elements in the dicionary;.
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.
int getColumnSpan() const
Returns the number of columns.
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...
unsigned char getByte(int n) const
Returns a specific byte of the pattern.
int getUsedBytes() const
Returns the number of bytes used for the pattern.
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
GLuint index
Definition glcorearb.h:781
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLboolean GLuint group
Definition glcorearb.h:3991
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
std::ostream & operator<<(std::ostream &os, const o2::header::RDHAny &rdh)
Definition RDHManip.cxx:84
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"