46 pitchtree->SetBranchAddress(
"PitchV", &
m_pitchV);
49 for (
int i = 0; i < pitchtree->GetEntries(); ++i) {
50 pitchtree->GetEntry(i);
56 typedef tuple<int, int, int> bufferkey_t;
57 typedef pair<PXDClusterShapeIndexPar, PXDClusterShapeClassifierPar> buffervalue_t;
58 map<bufferkey_t, buffervalue_t> localCalibrationMap;
62 B2INFO(
"Start calibration of clusterkind=" << clusterKind <<
" ...");
64 string gridname = str(format(
"GridKind_%1%") % clusterKind);
67 for (
auto uBin = 1; uBin <= grid->GetXaxis()->GetNbins(); uBin++) {
68 for (
auto vBin = 1; vBin <= grid->GetYaxis()->GetNbins(); vBin++) {
71 auto thetaU = grid->GetXaxis()->GetBinCenter(uBin);
72 auto thetaV = grid->GetYaxis()->GetBinCenter(vBin);
75 B2INFO(
"Skip training estimator on thetaU=" << thetaU <<
", thetaV=" << thetaV);
78 B2INFO(
"Start training estimator on thetaU=" << thetaU <<
", thetaV=" << thetaV);
81 string treename = str(format(
"tree_%1%_%2%_%3%") % clusterKind % uBin % vBin);
87 bufferkey_t key = std::make_tuple(clusterKind, uBin, vBin);
88 localCalibrationMap[key] = buffervalue_t(localShapeIndexer, localShapeClassifier);
96 int globalShapeIndex = 0;
98 shapeIndexer->
addShape(*it, globalShapeIndex);
102 B2INFO(
"Number of cluster shapes is " << globalShapeIndex);
112 string gridname = str(format(
"GridKind_%1%") % clusterKind);
115 positionEstimator->
addGrid(clusterKind, *grid);
117 for (
auto uBin = 1; uBin <= grid->GetXaxis()->GetNbins(); uBin++) {
118 for (
auto vBin = 1; vBin <= grid->GetYaxis()->GetNbins(); vBin++) {
121 auto thetaU = grid->GetXaxis()->GetBinCenter(uBin);
122 auto thetaV = grid->GetYaxis()->GetBinCenter(vBin);
130 auto iter = localCalibrationMap.find(std::make_tuple(clusterKind, uBin, vBin));
131 auto localShapeIndexer = iter->second.first;
132 auto localShapeClassifer = iter->second.second;
135 auto shapeClassifier =
localToGlobal(&localShapeClassifer, &localShapeIndexer, shapeIndexer);
138 auto mirror_vBin = grid->GetYaxis()->FindBin(-thetaV);
142 B2INFO(
"Add shape classifier for angles thetaU=" << thetaU <<
", thetaV=" << thetaV <<
", clusterkind=" << clusterKind);
144 B2INFO(
"Add mirrored shape classifier for angles thetaU=" << thetaU <<
", thetaV=" << -thetaV <<
", clusterkind=" << clusterKind);
145 positionEstimator->
setShapeClassifier(mirroredClassifier, uBin, mirror_vBin, clusterKind);
153 B2INFO(
"PXDClusterPosition Calibration Successful");
166 for (
auto indexAndValue : shapeLikelyhoodMap) {
168 auto shapeIndex = indexAndValue.first;
169 auto shapeName = shapeIndexer->
getShapeName(shapeIndex);
171 auto mirroredIndex = shapeIndexer->
getShapeIndex(mirroredName);
173 mirroredShapeClassifier.addShapeLikelyhood(mirroredIndex, indexAndValue.second);
180 for (
auto indexAndValue : offsetMap) {
182 auto shapeIndex = indexAndValue.first;
183 auto shapeName = shapeIndexer->
getShapeName(shapeIndex);
185 auto mirroredIndex = shapeIndexer->
getShapeIndex(mirroredName);
187 mirroredShapeClassifier.addShape(mirroredIndex);
190 for (
auto offset : indexAndValue.second) {
192 auto percentile = percentileMap[shapeIndex][etaBin];
193 mirroredShapeClassifier.addEtaPercentile(mirroredIndex, percentile);
195 auto likelyhood = likelyhoodMap[shapeIndex][etaBin];
196 mirroredShapeClassifier.addEtaLikelyhood(mirroredIndex, likelyhood);
199 auto mirroredOffset =
PXDClusterOffsetPar(offset.getU(), shift - offset.getV(), offset.getUSigma2(), offset.getVSigma2(),
200 -offset.getUVCovariance());
201 mirroredShapeClassifier.addEtaOffset(mirroredIndex, mirroredOffset);
206 return mirroredShapeClassifier;
217 for (
auto indexAndValue : shapeLikelyhoodMap) {
219 auto shapeIndex = indexAndValue.first;
220 auto shapeName = shapeIndexer->
getShapeName(shapeIndex);
221 auto globalIndex = globalShapeIndexer->
getShapeIndex(shapeName);
223 globalShapeClassifier.addShapeLikelyhood(globalIndex, indexAndValue.second);
230 for (
auto indexAndValue : offsetMap) {
232 auto shapeIndex = indexAndValue.first;
233 auto shapeName = shapeIndexer->
getShapeName(shapeIndex);
234 auto globalIndex = globalShapeIndexer->
getShapeIndex(shapeName);
236 globalShapeClassifier.addShape(globalIndex);
239 for (
auto offset : indexAndValue.second) {
241 auto percentile = percentileMap[shapeIndex][etaBin];
242 globalShapeClassifier.addEtaPercentile(globalIndex, percentile);
244 auto likelyhood = likelyhoodMap[shapeIndex][etaBin];
245 globalShapeClassifier.addEtaLikelyhood(globalIndex, likelyhood);
247 globalShapeClassifier.addEtaOffset(globalIndex, offset);
251 return globalShapeClassifier;
260 const auto nEntries = tree->GetEntries();
261 B2INFO(
"Number of clusters is " << nEntries);
266 tree->SetBranchAddress(
"ShapeName", &shapeNamePtr);
267 tree->SetBranchAddress(
"MirroredShapeName", &mirroredShapeNamePtr);
271 tree->SetBranchAddress(
"SizeV", &
m_sizeV);
275 vector< pair<string, float> > shapeList;
277 for (
int i = 0; i < nEntries; ++i) {
282 auto it = std::find_if(shapeList.begin(), shapeList.end(),
283 [&](
const pair<string, float>& element) { return element.first == shapeName;});
286 if (it != shapeList.end()) {
293 shapeList.push_back(pair<string, int>(shapeName, 1));
305 vector< pair<string, TH1D> > etaHistos;
311 double coverage = 0.0;
313 for (
auto iter : shapeList) {
314 auto name = iter.first;
315 auto counter = iter.second;
317 double likelyhood = counter / nEntries;
321 shapeIndexer->
addShape(name, tmpIndex);
333 coverage += likelyhood;
334 string etaname = str(format(
"eta_%1%") % name);
337 if (name ==
"SD0.0") {
338 TH1D etaHisto(etaname.c_str(), etaname.c_str(), 255, 0, 255);
339 etaHisto.SetDirectory(0);
340 etaHistos.push_back(pair<string, TH1D>(name, etaHisto));
344 TH1D etaHisto(etaname.c_str(), etaname.c_str(), 301, 0, 1);
345 etaHisto.SetDirectory(0);
346 etaHistos.push_back(pair<string, TH1D>(name, etaHisto));
351 B2INFO(
"Offset coverage is " << 100 * coverage <<
"%");
356 for (
int i = 0; i < nEntries; ++i) {
360 auto it = std::find_if(etaHistos.begin(), etaHistos.end(),
361 [&](
const pair<string, TH1D>& element) { return element.first == shapeName;});
363 if (it != etaHistos.end()) {
370 vector< pair< string, vector<TH2D> > > offsetHistosVec;
372 for (
auto iter : etaHistos) {
373 auto name = iter.first;
374 auto& histo = iter.second;
375 int nClusters = histo.GetEntries();
379 shapeClassifier->
addShape(shapeIndex);
387 vector< TH2D > offsetHistos;
389 for (
int i = 0; i < nEtaBins; i++) {
391 double xq = double(i) / nEtaBins;
394 histo.GetQuantiles(1, &yq, &xq);
398 string offsetname = str(format(
"offset_%1%_%2%") % name % i);
399 TH2D offsetHisto(offsetname.c_str(), offsetname.c_str(), 1, 0, 1, 1, 0, 1);
400 offsetHisto.SetDirectory(0);
401 offsetHisto.StatOverflows();
402 offsetHistos.push_back(offsetHisto);
405 offsetHistosVec.push_back(pair<
string, vector<TH2D> >(name, offsetHistos));
410 for (
int i = 0; i < nEntries; ++i) {
414 auto it = std::find_if(offsetHistosVec.begin(), offsetHistosVec.end(),
415 [&](
const pair<
string, vector<TH2D>>& element) { return element.first == shapeName;});
417 if (it != offsetHistosVec.end()) {
427 for (
auto iter : offsetHistosVec) {
428 auto name = iter.first;
429 auto& histovec = iter.second;
434 for (
auto& histo : histovec) {
436 double etaLikelyhood = double(histo.GetEntries()) / nEntries;
437 double offsetU = histo.GetMean(1);
438 double offsetV = histo.GetMean(2);
439 double covUV = histo.GetCovariance();
440 double covU = pow(histo.GetRMS(1), 2);
441 double covV = pow(histo.GetRMS(2), 2);
443 B2INFO(
"Name " << name <<
", posU=" << offsetU <<
", posV=" << offsetV <<
", covU=" << covU <<
", covV=" << covV <<
", covUV=" <<
446 TMatrixDSym HitCov(2);
448 HitCov(1, 0) = covUV;
449 HitCov(0, 1) = covUV;
452 TMatrixDSymEigen HitCovE(HitCov);
453 TVectorD eigenval = HitCovE.GetEigenValues();
454 if (eigenval(0) <= 0 || eigenval(1) <= 0) {
455 B2ERROR(
"Estimated covariance matrix not positive definite.");
464 B2INFO(
"Added shape classifier with coverage " << 100 * coverage <<
"% on training data sample.");