Belle II Software development
PXDClusterPositionCalibrationAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <pxd/calibration/PXDClusterPositionCalibrationAlgorithm.h>
10#include <pxd/dbobjects/PXDClusterShapeIndexPar.h>
11#include <pxd/dbobjects/PXDClusterPositionEstimatorPar.h>
12
13#include <string>
14#include <tuple>
15#include "TH2F.h"
16#include "TMatrixDSym.h"
17#include "TVectorD.h"
18#include "TMatrixDSymEigen.h"
19
20#include <boost/format.hpp>
21#include <cmath>
22
23using namespace std;
24using boost::format;
25using namespace Belle2;
26
28 CalibrationAlgorithm("PXDClusterPositionCollector"),
30 // Branches from TTree
33{
35 " -------------------------- PXDClusterPositionCalibrationAlgorithm ----------------------\n"
36 " \n"
37 " Algorithm for estimating cluster position offsets and shape likelyhoods. \n"
38 " ----------------------------------------------------------------------------------------\n"
39 );
40}
41
43{
44
45 // Read back the V pitch of all cluster kinds in source data
46 // This avoids relying on VXD::GeoCache.
47 auto pitchtree = getObjectPtr<TTree>("pitchtree");
48 pitchtree->SetBranchAddress("PitchV", &m_pitchV);
49 pitchtree->SetBranchAddress("ClusterKind", &m_clusterKind);
50
51 for (int i = 0; i < pitchtree->GetEntries(); ++i) {
52 pitchtree->GetEntry(i);
54 }
55
56 // Buffer temporary payloads for shape calibration for all
57 // clusterkinds and angle bins
58 typedef tuple<int, int, int> bufferkey_t;
59 typedef pair<PXDClusterShapeIndexPar, PXDClusterShapeClassifierPar> buffervalue_t;
60 map<bufferkey_t, buffervalue_t> localCalibrationMap;
61
62 for (auto clusterKind : clusterKinds) {
63
64 B2INFO("Start calibration of clusterkind=" << clusterKind << " ...");
65
66 string gridname = str(format("GridKind_%1%") % clusterKind);
67 auto grid = getObjectPtr<TH2F>(gridname);
68
69 for (auto uBin = 1; uBin <= grid->GetXaxis()->GetNbins(); uBin++) {
70 for (auto vBin = 1; vBin <= grid->GetYaxis()->GetNbins(); vBin++) {
71
72 // Bin is centered around angles
73 auto thetaU = grid->GetXaxis()->GetBinCenter(uBin);
74 auto thetaV = grid->GetYaxis()->GetBinCenter(vBin);
75
76 if (thetaV < 0) {
77 B2INFO("Skip training estimator on thetaU=" << thetaU << ", thetaV=" << thetaV);
78 continue;
79 } else {
80 B2INFO("Start training estimator on thetaU=" << thetaU << ", thetaV=" << thetaV);
81 }
82
83 string treename = str(format("tree_%1%_%2%_%3%") % clusterKind % uBin % vBin);
84
85 auto localShapeIndexer = PXDClusterShapeIndexPar();
86 auto localShapeClassifier = PXDClusterShapeClassifierPar();
87 createShapeClassifier(treename, &localShapeClassifier, &localShapeIndexer);
88
89 bufferkey_t key = std::make_tuple(clusterKind, uBin, vBin);
90 localCalibrationMap[key] = buffervalue_t(localShapeIndexer, localShapeClassifier);
91 }
92 }
93 }
94
95 // Create a ShapeIndexer payload
97
98 int globalShapeIndex = 0;
99 for (auto it = m_shapeSet.begin(); it != m_shapeSet.end(); ++it) {
100 shapeIndexer->addShape(*it, globalShapeIndex);
101 globalShapeIndex++;
102 }
103
104 B2INFO("Number of cluster shapes is " << globalShapeIndex);
105
106 // Save the cluster shape index table to database.
107 saveCalibration(shapeIndexer, "PXDClusterShapeIndexPar");
108
109 // Create position estimator
111
112 for (auto clusterKind : clusterKinds) {
113
114 string gridname = str(format("GridKind_%1%") % clusterKind);
115 auto grid = getObjectPtr<TH2F>(gridname);
116
117 positionEstimator->addGrid(clusterKind, *grid);
118
119 for (auto uBin = 1; uBin <= grid->GetXaxis()->GetNbins(); uBin++) {
120 for (auto vBin = 1; vBin <= grid->GetYaxis()->GetNbins(); vBin++) {
121
122 // Bin is centered around angles
123 auto thetaU = grid->GetXaxis()->GetBinCenter(uBin);
124 auto thetaV = grid->GetYaxis()->GetBinCenter(vBin);
125
126 if (thetaV < 0) {
127 // We skipped this part in training before
128 continue;
129 }
130
131 // Find the local calibration results
132 auto iter = localCalibrationMap.find(std::make_tuple(clusterKind, uBin, vBin));
133 auto localShapeIndexer = iter->second.first;
134 auto localShapeClassifer = iter->second.second;
135
136 // Require that all shape classifiers use a common shape indexer
137 auto shapeClassifier = localToGlobal(&localShapeClassifer, &localShapeIndexer, shapeIndexer);
138
139 // Mirror the shape classifier along v
140 auto mirror_vBin = grid->GetYaxis()->FindBin(-thetaV);
141 auto mirroredClassifier = mirrorShapeClassifier(&shapeClassifier, shapeIndexer, clusterKind);
142
143 // and fill into position estimator payload
144 B2INFO("Add shape classifier for angles thetaU=" << thetaU << ", thetaV=" << thetaV << ", clusterkind=" << clusterKind);
145 positionEstimator->setShapeClassifier(shapeClassifier, uBin, vBin, clusterKind);
146 B2INFO("Add mirrored shape classifier for angles thetaU=" << thetaU << ", thetaV=" << -thetaV << ", clusterkind=" << clusterKind);
147 positionEstimator->setShapeClassifier(mirroredClassifier, uBin, mirror_vBin, clusterKind);
148 }
149 }
150 }
151
152 // Save the cluster positions to database.
153 saveCalibration(positionEstimator, "PXDClusterPositionEstimatorPar");
154
155 B2INFO("PXDClusterPosition Calibration Successful");
156 return c_OK;
157}
158
159
161 shapeClassifier, PXDClusterShapeIndexPar* shapeIndexer, int clusterKind)
162{
163 // Create a mirrored shape classifier
164 auto mirroredShapeClassifier = PXDClusterShapeClassifierPar();
165
166 // Mirror the shape likelyhood map
167 auto shapeLikelyhoodMap = shapeClassifier->getShapeLikelyhoodMap();
168 for (auto indexAndValue : shapeLikelyhoodMap) {
169 // Compute the mirrored shape index
170 auto shapeIndex = indexAndValue.first;
171 auto shapeName = shapeIndexer->getShapeName(shapeIndex);
172 auto mirroredName = m_mirrorMap[shapeName];
173 auto mirroredIndex = shapeIndexer->getShapeIndex(mirroredName);
174 // Store the result
175 mirroredShapeClassifier.addShapeLikelyhood(mirroredIndex, indexAndValue.second);
176 }
177
178 // Mirror the offset related maps
179 auto offsetMap = shapeClassifier->getOffsetMap();
180 auto percentileMap = shapeClassifier->getPercentilesMap();
181 auto likelyhoodMap = shapeClassifier->getLikelyhoodMap();
182 for (auto indexAndValue : offsetMap) {
183 // Compute the mirrored shape index
184 auto shapeIndex = indexAndValue.first;
185 auto shapeName = shapeIndexer->getShapeName(shapeIndex);
186 auto mirroredName = m_mirrorMap[shapeName];
187 auto mirroredIndex = shapeIndexer->getShapeIndex(mirroredName);
188
189 mirroredShapeClassifier.addShape(mirroredIndex);
190
191 int etaBin = 0;
192 for (auto offset : indexAndValue.second) {
193 // Copy over percentile
194 auto percentile = percentileMap[shapeIndex][etaBin];
195 mirroredShapeClassifier.addEtaPercentile(mirroredIndex, percentile);
196 // Copy over likelyhood
197 auto likelyhood = likelyhoodMap[shapeIndex][etaBin];
198 mirroredShapeClassifier.addEtaLikelyhood(mirroredIndex, likelyhood);
199 // Mirror the offset: v offset shifts and covariance swaps sign
200 double shift = (m_sizeMap[shapeName] - 1) * m_pitchMap[clusterKind];
201 auto mirroredOffset = PXDClusterOffsetPar(offset.getU(), shift - offset.getV(), offset.getUSigma2(), offset.getVSigma2(),
202 -offset.getUVCovariance());
203 mirroredShapeClassifier.addEtaOffset(mirroredIndex, mirroredOffset);
204 etaBin++;
205 }
206 }
207
208 return mirroredShapeClassifier;
209}
210
212 PXDClusterShapeIndexPar* shapeIndexer, PXDClusterShapeIndexPar* globalShapeIndexer)
213{
214 // Create a shape classifier using global shape indices
215 auto globalShapeClassifier = PXDClusterShapeClassifierPar();
216
217 // Re-index the the shape likelyhood map
218 auto shapeLikelyhoodMap = shapeClassifier->getShapeLikelyhoodMap();
219 for (auto indexAndValue : shapeLikelyhoodMap) {
220 // Compute the global shape index
221 auto shapeIndex = indexAndValue.first;
222 auto shapeName = shapeIndexer->getShapeName(shapeIndex);
223 auto globalIndex = globalShapeIndexer->getShapeIndex(shapeName);
224 // Store the result
225 globalShapeClassifier.addShapeLikelyhood(globalIndex, indexAndValue.second);
226 }
227
228 // Re-index the offset related maps
229 auto offsetMap = shapeClassifier->getOffsetMap();
230 auto percentileMap = shapeClassifier->getPercentilesMap();
231 auto likelyhoodMap = shapeClassifier->getLikelyhoodMap();
232 for (auto indexAndValue : offsetMap) {
233 // Compute the global shape index
234 auto shapeIndex = indexAndValue.first;
235 auto shapeName = shapeIndexer->getShapeName(shapeIndex);
236 auto globalIndex = globalShapeIndexer->getShapeIndex(shapeName);
237
238 globalShapeClassifier.addShape(globalIndex);
239
240 int etaBin = 0;
241 for (auto offset : indexAndValue.second) {
242 // Copy over percentile
243 auto percentile = percentileMap[shapeIndex][etaBin];
244 globalShapeClassifier.addEtaPercentile(globalIndex, percentile);
245 // Copy over likelyhood
246 auto likelyhood = likelyhoodMap[shapeIndex][etaBin];
247 globalShapeClassifier.addEtaLikelyhood(globalIndex, likelyhood);
248 // Copy over offset
249 globalShapeClassifier.addEtaOffset(globalIndex, offset);
250 etaBin++;
251 }
252 }
253 return globalShapeClassifier;
254}
255
257 PXDClusterShapeClassifierPar* shapeClassifier, PXDClusterShapeIndexPar* shapeIndexer)
258{
259
260 auto tree = getObjectPtr<TTree>(treename);
261
262 const auto nEntries = tree->GetEntries();
263 B2INFO("Number of clusters is " << nEntries);
264
265 string* shapeNamePtr = &m_shapeName;
266 string* mirroredShapeNamePtr = &m_mirroredShapeName;
267
268 tree->SetBranchAddress("ShapeName", &shapeNamePtr);
269 tree->SetBranchAddress("MirroredShapeName", &mirroredShapeNamePtr);
270 tree->SetBranchAddress("ClusterEta", &m_clusterEta);
271 tree->SetBranchAddress("OffsetU", &m_positionOffsetU);
272 tree->SetBranchAddress("OffsetV", &m_positionOffsetV);
273 tree->SetBranchAddress("SizeV", &m_sizeV);
274
275 // Vector to enumerate all shapes by unique name and count their
276 // occurrence in training data.
277 vector< pair<string, float> > shapeList;
278
279 for (int i = 0; i < nEntries; ++i) {
280 tree->GetEntry(i);
281
282 string shapeName = m_shapeName;
283
284 auto it = std::find_if(shapeList.begin(), shapeList.end(),
285 [&](const pair<string, float>& element) { return element.first == shapeName;});
286
287 //Shape name exists in vector
288 if (it != shapeList.end()) {
289 //increment key in map
290 it->second++;
291 }
292 //Shape name does not exist
293 else {
294 //Not found, insert in vector
295 shapeList.push_back(pair<string, int>(shapeName, 1));
296 // Remember the relation between name of a shape and name of mirrored shape
297 m_mirrorMap[shapeName] = m_mirroredShapeName;
298 // Remember the relation between name of a shape and its size
299 m_sizeMap[shapeName] = m_sizeV;
300 }
301 }
302
303 // Loop over shapeList to select shapes for
304 // next calibration step
305
306 // Vector with eta histograms for selected shapes
307 vector< pair<string, TH1D> > etaHistos;
308
309 // Index for enumerating selected shapes
310 int tmpIndex = 0;
311
312 // Coverage of position offsets on training data
313 double coverage = 0.0;
314
315 for (auto iter : shapeList) {
316 auto name = iter.first;
317 auto counter = iter.second;
318
319 double likelyhood = counter / nEntries;
320
321 if (counter >= minClusterForShapeLikelyhood) {
322 //B2INFO("Adding shape " << name << " with index " << tmpIndex << " and shape likelyhood " << 100*likelyhood << "% and count " << counter);
323 shapeIndexer->addShape(name, tmpIndex);
324 shapeClassifier->addShapeLikelyhood(tmpIndex, likelyhood);
325
326 // Add name of shape to global (all clusterkinds + all angle bins) shape set
327 m_shapeSet.insert(name);
328 // Add name of mirrored shape as well
329 m_shapeSet.insert(m_mirrorMap[name]);
330 // Increment the index
331 tmpIndex++;
332 }
333
334 if (counter >= minClusterForPositionOffset) {
335 coverage += likelyhood;
336 string etaname = str(format("eta_%1%") % name);
337
338 // Single pixel case: Eta value is cluster charge
339 if (name == "SD0.0") {
340 TH1D etaHisto(etaname.c_str(), etaname.c_str(), 255, 0, 255);
341 etaHisto.SetDirectory(0);
342 etaHistos.push_back(pair<string, TH1D>(name, etaHisto));
343 }
344 // Multipixel case: Eta value is ratio head/(tail+head) of charges (to be less gain sensitive)
345 else {
346 TH1D etaHisto(etaname.c_str(), etaname.c_str(), 301, 0, 1);
347 etaHisto.SetDirectory(0);
348 etaHistos.push_back(pair<string, TH1D>(name, etaHisto));
349 }
350 }
351 }
352
353 B2INFO("Offset coverage is " << 100 * coverage << "%");
354
355 // Loop over the tree is to fill the eta histograms for
356 // selected shapes.
357
358 for (int i = 0; i < nEntries; ++i) {
359 tree->GetEntry(i);
360
361 string shapeName = m_shapeName;
362 auto it = std::find_if(etaHistos.begin(), etaHistos.end(),
363 [&](const pair<string, TH1D>& element) { return element.first == shapeName;});
364 //Item exists in map
365 if (it != etaHistos.end()) {
366 // increment key in map
367 it->second.Fill(m_clusterEta);
368 }
369 }
370
371 // Vector for offset histograms stored by offset shape name and eta bin
372 vector< pair< string, vector<TH2D> > > offsetHistosVec;
373
374 for (auto iter : etaHistos) {
375 auto name = iter.first;
376 auto& histo = iter.second;
377 int nClusters = histo.GetEntries();
378
379 // Add shape for offset correction
380 int shapeIndex = shapeIndexer->getShapeIndex(name);
381 shapeClassifier->addShape(shapeIndex);
382
383 // Try to split clusters into n bins with minClusterForPositionOffset clusters
384 int nEtaBins = std::max(int(nClusters / minClusterForPositionOffset), 1);
385 nEtaBins = std::min(nEtaBins, maxEtaBins);
386
387 //B2INFO("SHAPE NAME:" << name << " WITH BINS " << nEtaBins);
388
389 vector< TH2D > offsetHistos;
390
391 for (int i = 0; i < nEtaBins; i++) {
392 // Position where to compute the quantiles in [0,1]
393 double xq = double(i) / nEtaBins;
394 // Double to contain the quantile
395 double yq = 0;
396 histo.GetQuantiles(1, &yq, &xq);
397 //B2INFO(" Quantile at xq =" << xq << " is yq=" << yq);
398 shapeClassifier->addEtaPercentile(shapeIndex, yq);
399
400 string offsetname = str(format("offset_%1%_%2%") % name % i);
401 TH2D offsetHisto(offsetname.c_str(), offsetname.c_str(), 1, 0, 1, 1, 0, 1);
402 offsetHisto.SetDirectory(0);
403 offsetHisto.StatOverflows();
404 offsetHistos.push_back(offsetHisto);
405
406 }
407 offsetHistosVec.push_back(pair< string, vector<TH2D> >(name, offsetHistos));
408 }
409
410 // Loop over the tree is to fill offset histograms
411
412 for (int i = 0; i < nEntries; ++i) {
413 tree->GetEntry(i);
414
415 string shapeName = m_shapeName;
416 auto it = std::find_if(offsetHistosVec.begin(), offsetHistosVec.end(),
417 [&](const pair<string, vector<TH2D>>& element) { return element.first == shapeName;});
418 //Item exists in map
419 if (it != offsetHistosVec.end()) {
420 int shapeIndex = shapeIndexer->getShapeIndex(shapeName);
421 int etaBin = shapeClassifier->getEtaIndex(shapeIndex, m_clusterEta);
422 it->second.at(etaBin).Fill(m_positionOffsetU, m_positionOffsetV);
423 }
424 }
425
426 // Compute the moments of the offset histograms and finalize the shape classifier object
427
428 // Loop over shape names
429 for (auto iter : offsetHistosVec) {
430 auto name = iter.first;
431 auto& histovec = iter.second;
432
433 int shapeIndex = shapeIndexer->getShapeIndex(name);
434
435 // Loop over eta bins
436 for (auto& histo : histovec) {
437 // Compute offset moments
438 double etaLikelyhood = double(histo.GetEntries()) / nEntries;
439 double offsetU = histo.GetMean(1);
440 double offsetV = histo.GetMean(2);
441 double covUV = histo.GetCovariance();
442 double covU = pow(histo.GetRMS(1), 2);
443 double covV = pow(histo.GetRMS(2), 2);
444
445 B2INFO("Name " << name << ", posU=" << offsetU << ", posV=" << offsetV << ", covU=" << covU << ", covV=" << covV << ", covUV=" <<
446 covUV);
447
448 TMatrixDSym HitCov(2);
449 HitCov(0, 0) = covU;
450 HitCov(1, 0) = covUV;
451 HitCov(0, 1) = covUV;
452 HitCov(1, 1) = covV;
453
454 TMatrixDSymEigen HitCovE(HitCov);
455 TVectorD eigenval = HitCovE.GetEigenValues();
456 if (eigenval(0) <= 0 || eigenval(1) <= 0) {
457 B2ERROR("Estimated covariance matrix not positive definite.");
458 }
459
460 auto offset = PXDClusterOffsetPar(offsetU, offsetV, covU, covV, covUV);
461 shapeClassifier->addEtaLikelyhood(shapeIndex, etaLikelyhood);
462 shapeClassifier->addEtaOffset(shapeIndex, offset);
463 }
464 }
465
466 B2INFO("Added shape classifier with coverage " << 100 * coverage << "% on training data sample.");
467
468 return;
469}
470
471
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
The class for PXD cluster position offset payload.
void createShapeClassifier(std::string treename, PXDClusterShapeClassifierPar *shapeClassifier, PXDClusterShapeIndexPar *shapeIndexer)
Returns a new classifier and index trained on cluster tree.
std::vector< int > clusterKinds
Vector of clusterkinds to calibrate.
PXDClusterShapeClassifierPar mirrorShapeClassifier(PXDClusterShapeClassifierPar *shapeClassifier, PXDClusterShapeIndexPar *shapeIndexer, int clusterKind)
Returns a mirrored version of shape classifier.
PXDClusterShapeClassifierPar localToGlobal(PXDClusterShapeClassifierPar *localShapeClassifier, PXDClusterShapeIndexPar *localShapeIndexer, PXDClusterShapeIndexPar *globalShapeIndexer)
Returns a shape classifier using global shape indices instead of local ones.
std::map< std::string, int > m_sizeMap
Helper needed to map the name of a shape to the V size of the cluster.
int maxEtaBins
Maximum number of eta bins for estimating cluster position offsets.
int minClusterForPositionOffset
Minimum number of collected clusters for estimating cluster position offsets.
int minClusterForShapeLikelyhood
Minimum number of collected clusters for estimating shape likelyhood.
PXDClusterPositionCalibrationAlgorithm()
Constructor set the prefix to PXDClusterPositionCalibrationAlgorithm.
std::map< int, float > m_pitchMap
Helper needed to map the clusterkind to the V pitch of the sensor.
std::set< std::string > m_shapeSet
Set of unique shape names.
std::map< std::string, std::string > m_mirrorMap
Helper needed to map the name of a shape to the name of the mirrored shape.
The class for PXD cluster position lookup table payload.
void setShapeClassifier(const PXDClusterShapeClassifierPar &classifier, int uBin, int vBin, int clusterkind)
Set shape classifier.
void addGrid(int clusterkind, const TH2F &grid)
Add grid for clusterkind.
The class for PXD cluster shape classifier payload.
void addShapeLikelyhood(int shape_index, float likelyhood)
Add shape likelihood.
const std::map< int, std::vector< PXDClusterOffsetPar > > & getOffsetMap() const
Return offset map for position correction.
void addEtaLikelyhood(int shape_index, float likelyhood)
Add eta likelyhood to shape for position correction.
const std::map< int, std::vector< float > > & getPercentilesMap() const
Return percentiles map for position correction.
void addEtaPercentile(int shape_index, float percentile)
Add eta percentile to shape for position correction.
const std::map< int, std::vector< float > > & getLikelyhoodMap() const
Return likelyhood map for position correction.
unsigned int getEtaIndex(int shape_index, float eta) const
Get eta index for position correction.
void addShape(int shape_index)
Add shape for position correction.
const std::map< int, float > & getShapeLikelyhoodMap() const
Return shape likelyhood map.
void addEtaOffset(int shape_index, PXDClusterOffsetPar &offset)
Add offset to shape for position correction.
The class for PXD cluster shape index payload.
int getShapeIndex(const std::string &name) const
Returns shape index from name.
const std::string & getShapeName(int index) const
Returns shape name from index.
void addShape(const std::string &name, int index)
Add shape with name and index.
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Abstract base class for different kinds of events.
STL namespace.