File indexing completed on 2021-04-29 02:32:49
0001
0002 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h"
0003
0004
0005 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0006 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0007
0008
0009
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015
0016
0017 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateReco.h"
0018
0019
0020
0021
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023
0024 #include <vector>
0025 #include "boost/multi_array.hpp"
0026
0027 #include <iostream>
0028
0029 using namespace SiPixelTemplateReco;
0030
0031 using namespace std;
0032
0033 namespace {
0034 constexpr float micronsToCm = 1.0e-4;
0035 constexpr int cluster_matrix_size_x = 13;
0036 constexpr int cluster_matrix_size_y = 21;
0037 }
0038
0039
0040
0041
0042
0043 PixelCPETemplateReco::PixelCPETemplateReco(edm::ParameterSet const& conf,
0044 const MagneticField* mag,
0045 const TrackerGeometry& geom,
0046 const TrackerTopology& ttopo,
0047 const SiPixelLorentzAngle* lorentzAngle,
0048 const SiPixelTemplateDBObject* templateDBobject)
0049 : PixelCPEBase(conf, mag, geom, ttopo, lorentzAngle, nullptr, templateDBobject, nullptr, 1) {
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063 if (LoadTemplatesFromDB_) {
0064
0065
0066
0067 if (!SiPixelTemplate::pushfile(*templateDBobject_, thePixelTemp_))
0068 throw cms::Exception("PixelCPETemplateReco")
0069 << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
0070 << (*templateDBobject_).version() << "\n\n";
0071 } else {
0072
0073 barrelTemplateID_ = conf.getParameter<int>("barrelTemplateID");
0074 forwardTemplateID_ = conf.getParameter<int>("forwardTemplateID");
0075 templateDir_ = conf.getParameter<int>("directoryWithTemplates");
0076
0077 if (!SiPixelTemplate::pushfile(barrelTemplateID_, thePixelTemp_, templateDir_))
0078 throw cms::Exception("PixelCPETemplateReco")
0079 << "\nERROR: Template ID " << barrelTemplateID_
0080 << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0081
0082 if (!SiPixelTemplate::pushfile(forwardTemplateID_, thePixelTemp_, templateDir_))
0083 throw cms::Exception("PixelCPETemplateReco")
0084 << "\nERROR: Template ID " << forwardTemplateID_
0085 << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0086 }
0087
0088 speed_ = conf.getParameter<int>("speed");
0089 LogDebug("PixelCPETemplateReco::PixelCPETemplateReco:") << "Template speed = " << speed_ << "\n";
0090
0091 UseClusterSplitter_ = conf.getParameter<bool>("UseClusterSplitter");
0092 }
0093
0094
0095
0096
0097 PixelCPETemplateReco::~PixelCPETemplateReco() {}
0098
0099 std::unique_ptr<PixelCPEBase::ClusterParam> PixelCPETemplateReco::createClusterParam(const SiPixelCluster& cl) const {
0100 return std::make_unique<ClusterParamTemplate>(cl);
0101 }
0102
0103
0104
0105
0106
0107
0108
0109
0110 LocalPoint PixelCPETemplateReco::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0111 ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0112
0113 if (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0114 throw cms::Exception("PixelCPETemplateReco::localPosition :") << "A non-pixel detector type in here?";
0115
0116 const bool fpix = GeomDetEnumerators::isEndcap(theDetParam.thePart);
0117
0118 int ID = -9999;
0119 if (LoadTemplatesFromDB_) {
0120 int ID0 = templateDBobject_->getTemplateID(theDetParam.theDet->geographicalId());
0121 ID = theDetParam.detTemplateId;
0122 if (ID0 != ID)
0123 edm::LogError("PixelCPETemplateReco") << " different id" << ID << " " << ID0 << endl;
0124 } else {
0125 if (!fpix)
0126 ID = barrelTemplateID_;
0127 else
0128 ID = forwardTemplateID_;
0129 }
0130
0131
0132 SiPixelTemplate templ(thePixelTemp_);
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143 int row_offset = theClusterParam.theCluster->minPixelRow();
0144 int col_offset = theClusterParam.theCluster->minPixelCol();
0145
0146
0147
0148
0149 float tmp_x = float(row_offset) + 0.5f;
0150 float tmp_y = float(col_offset) + 0.5f;
0151
0152
0153
0154
0155
0156
0157
0158
0159 float lorentzshiftX = 0.5f * theDetParam.lorentzShiftInCmX;
0160 float lorentzshiftY = 0.5f * theDetParam.lorentzShiftInCmY;
0161
0162
0163 LocalPoint lp;
0164
0165 if (theClusterParam.with_track_angle)
0166 lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y), theClusterParam.loc_trk_pred);
0167 else {
0168 edm::LogError("PixelCPETemplateReco") << "@SUB = PixelCPETemplateReco::localPosition"
0169 << "Should never be here. PixelCPETemplateReco should always be called with "
0170 "track angles. This is a bad error !!! ";
0171
0172 lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y));
0173 }
0174
0175
0176 int mrow = 0, mcol = 0;
0177 for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0178 auto pix = theClusterParam.theCluster->pixel(i);
0179 int irow = int(pix.x);
0180 int icol = int(pix.y);
0181 mrow = std::max(mrow, irow);
0182 mcol = std::max(mcol, icol);
0183 }
0184 mrow -= row_offset;
0185 mrow += 1;
0186 mrow = std::min(mrow, cluster_matrix_size_x);
0187 mcol -= col_offset;
0188 mcol += 1;
0189 mcol = std::min(mcol, cluster_matrix_size_y);
0190 assert(mrow > 0);
0191 assert(mcol > 0);
0192
0193 float clustMatrix[mrow][mcol];
0194 memset(clustMatrix, 0, sizeof(float) * mrow * mcol);
0195
0196
0197 for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0198 auto pix = theClusterParam.theCluster->pixel(i);
0199 int irow = int(pix.x) - row_offset;
0200 int icol = int(pix.y) - col_offset;
0201
0202
0203
0204 if ((irow < mrow) & (icol < mcol))
0205 clustMatrix[irow][icol] = float(pix.adc);
0206 }
0207
0208
0209 bool xdouble[mrow], ydouble[mcol];
0210
0211 for (int irow = 0; irow < mrow; ++irow)
0212 xdouble[irow] = theDetParam.theRecTopol->isItBigPixelInX(irow + row_offset);
0213
0214
0215 for (int icol = 0; icol < mcol; ++icol)
0216 ydouble[icol] = theDetParam.theRecTopol->isItBigPixelInY(icol + col_offset);
0217
0218 SiPixelTemplateReco::ClusMatrix clusterPayload{&clustMatrix[0][0], xdouble, ydouble, mrow, mcol};
0219
0220
0221 float nonsense = -99999.9f;
0222 theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ =
0223 theClusterParam.templSigmaY_ = nonsense;
0224
0225 theClusterParam.templProbY_ = theClusterParam.templProbX_ = theClusterParam.templProbQ_ = 1.0f;
0226 theClusterParam.templQbin_ = 0;
0227
0228 theClusterParam.hasFilledProb_ = false;
0229
0230 float templYrec1_ = nonsense;
0231 float templXrec1_ = nonsense;
0232 float templYrec2_ = nonsense;
0233 float templXrec2_ = nonsense;
0234
0235
0236
0237
0238 float locBz = theDetParam.bz;
0239 float locBx = theDetParam.bx;
0240
0241 theClusterParam.ierr = PixelTempReco1D(ID,
0242 theClusterParam.cotalpha,
0243 theClusterParam.cotbeta,
0244 locBz,
0245 locBx,
0246 clusterPayload,
0247 templ,
0248 theClusterParam.templYrec_,
0249 theClusterParam.templSigmaY_,
0250 theClusterParam.templProbY_,
0251 theClusterParam.templXrec_,
0252 theClusterParam.templSigmaX_,
0253 theClusterParam.templProbX_,
0254 theClusterParam.templQbin_,
0255 speed_,
0256 theClusterParam.templProbQ_);
0257
0258
0259
0260
0261 if UNLIKELY (theClusterParam.ierr != 0) {
0262 LogDebug("PixelCPETemplateReco::localPosition")
0263 << "reconstruction failed with error " << theClusterParam.ierr << "\n";
0264
0265
0266
0267
0268
0269 if (theClusterParam.with_track_angle) {
0270 theClusterParam.templXrec_ =
0271 theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) + lorentzshiftX;
0272 theClusterParam.templYrec_ =
0273 theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred) + lorentzshiftY;
0274 } else {
0275 edm::LogError("PixelCPETemplateReco") << "@SUB = PixelCPETemplateReco::localPosition"
0276 << "Should never be here. PixelCPETemplateReco should always be called "
0277 "with track angles. This is a bad error !!! ";
0278
0279 theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) + lorentzshiftX;
0280 theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y()) + lorentzshiftY;
0281 }
0282 } else if UNLIKELY (UseClusterSplitter_ && theClusterParam.templQbin_ == 0) {
0283 edm::LogError("PixelCPETemplateReco") << " PixelCPETemplateReco: Qbin = 0 but using cluster splitter, we should "
0284 "never be here !!!!!!!!!!!!!!!!!!!!!! \n"
0285 << "(int)UseClusterSplitter_ = " << (int)UseClusterSplitter_ << endl;
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300 theClusterParam.ierr = -123;
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317 if (theClusterParam.ierr != 0) {
0318
0319
0320
0321
0322 if (theClusterParam.with_track_angle) {
0323 theClusterParam.templXrec_ =
0324 theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) + lorentzshiftX;
0325 theClusterParam.templYrec_ =
0326 theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred) + lorentzshiftY;
0327 } else {
0328 edm::LogError("PixelCPETemplateReco") << "@SUB = PixelCPETemplateReco::localPosition"
0329 << "Should never be here. PixelCPETemplateReco should always be called "
0330 "with track angles. This is a bad error !!! ";
0331 theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) + lorentzshiftX;
0332 theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y()) + lorentzshiftY;
0333 }
0334 } else {
0335
0336 templXrec1_ *= micronsToCm;
0337 templYrec1_ *= micronsToCm;
0338 templXrec2_ *= micronsToCm;
0339 templYrec2_ *= micronsToCm;
0340
0341
0342 templXrec1_ += lp.x();
0343 templYrec1_ += lp.y();
0344 templXrec2_ += lp.x();
0345 templYrec2_ += lp.y();
0346
0347
0348 float distX1 = std::abs(templXrec1_ - theClusterParam.trk_lp_x);
0349 float distX2 = std::abs(templXrec2_ - theClusterParam.trk_lp_x);
0350 float distY1 = std::abs(templYrec1_ - theClusterParam.trk_lp_y);
0351 float distY2 = std::abs(templYrec2_ - theClusterParam.trk_lp_y);
0352 theClusterParam.templXrec_ = (distX1 < distX2 ? templXrec1_ : templXrec2_);
0353 theClusterParam.templYrec_ = (distY1 < distY2 ? templYrec1_ : templYrec2_);
0354 }
0355 }
0356
0357 else
0358 {
0359
0360 theClusterParam.templXrec_ *= micronsToCm;
0361 theClusterParam.templYrec_ *= micronsToCm;
0362
0363
0364 theClusterParam.templXrec_ += lp.x();
0365 theClusterParam.templYrec_ += lp.y();
0366
0367
0368 if (doLorentzFromAlignment_) {
0369
0370 if (theDetParam.lorentzShiftInCmX != 0.0 || theDetParam.lorentzShiftInCmY != 0.0) {
0371
0372
0373
0374
0375
0376 float templateLorbiasCmX = -micronsToCm * templ.lorxbias();
0377 float templateLorbiasCmY = -micronsToCm * templ.lorybias();
0378
0379
0380
0381 theClusterParam.templXrec_ += (0.5 * (theDetParam.lorentzShiftInCmX) - templateLorbiasCmX);
0382 theClusterParam.templYrec_ += (0.5 * (theDetParam.lorentzShiftInCmY) - templateLorbiasCmY);
0383
0384
0385
0386
0387
0388
0389
0390 }
0391 }
0392 }
0393
0394
0395
0396 theClusterParam.probabilityX_ = theClusterParam.templProbX_;
0397 theClusterParam.probabilityY_ = theClusterParam.templProbY_;
0398 theClusterParam.probabilityQ_ = theClusterParam.templProbQ_;
0399 theClusterParam.qBin_ = theClusterParam.templQbin_;
0400
0401 if (theClusterParam.ierr == 0)
0402 theClusterParam.hasFilledProb_ = true;
0403
0404 return LocalPoint(theClusterParam.templXrec_, theClusterParam.templYrec_);
0405 }
0406
0407
0408
0409
0410 LocalError PixelCPETemplateReco::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0411 ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423 float xerr, yerr;
0424
0425
0426 if (theClusterParam.theCluster->getSplitClusterErrorX() > 0.0f &&
0427 theClusterParam.theCluster->getSplitClusterErrorX() < clusterSplitMaxError_ &&
0428 theClusterParam.theCluster->getSplitClusterErrorY() > 0.0f &&
0429 theClusterParam.theCluster->getSplitClusterErrorY() < clusterSplitMaxError_) {
0430 xerr = theClusterParam.theCluster->getSplitClusterErrorX() * micronsToCm;
0431 yerr = theClusterParam.theCluster->getSplitClusterErrorY() * micronsToCm;
0432
0433
0434
0435
0436 } else {
0437
0438
0439
0440
0441 int maxPixelCol = theClusterParam.theCluster->maxPixelCol();
0442 int maxPixelRow = theClusterParam.theCluster->maxPixelRow();
0443 int minPixelCol = theClusterParam.theCluster->minPixelCol();
0444 int minPixelRow = theClusterParam.theCluster->minPixelRow();
0445
0446
0447 bool edgex = (theDetParam.theRecTopol->isItEdgePixelInX(minPixelRow) ||
0448 theDetParam.theRecTopol->isItEdgePixelInX(maxPixelRow));
0449 bool edgey = (theDetParam.theRecTopol->isItEdgePixelInY(minPixelCol) ||
0450 theDetParam.theRecTopol->isItEdgePixelInY(maxPixelCol));
0451
0452 if (theClusterParam.ierr != 0) {
0453
0454
0455
0456
0457
0458 if (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0459 throw cms::Exception("PixelCPETemplateReco::localPosition :") << "A non-pixel detector type in here?";
0460
0461
0462 if (GeomDetEnumerators::isBarrel(theDetParam.thePart)) {
0463 xerr = 55.0f * micronsToCm;
0464 yerr = 36.0f * micronsToCm;
0465 } else {
0466 xerr = 42.0f * micronsToCm;
0467 yerr = 39.0f * micronsToCm;
0468 }
0469
0470
0471
0472
0473
0474 } else if (edgex || edgey) {
0475
0476 if (edgex && !edgey) {
0477 xerr = xEdgeXError_ * micronsToCm;
0478 yerr = xEdgeYError_ * micronsToCm;
0479 } else if (!edgex && edgey) {
0480 xerr = yEdgeXError_ * micronsToCm;
0481 yerr = yEdgeYError_ * micronsToCm;
0482 } else if (edgex && edgey) {
0483 xerr = bothEdgeXError_ * micronsToCm;
0484 yerr = bothEdgeYError_ * micronsToCm;
0485 } else {
0486 throw cms::Exception(" PixelCPETemplateReco::localError: Something wrong with pixel edge flag !!!");
0487 }
0488
0489
0490
0491 } else {
0492
0493
0494
0495 xerr = theClusterParam.templSigmaX_ * micronsToCm;
0496 yerr = theClusterParam.templSigmaY_ * micronsToCm;
0497
0498
0499
0500
0501
0502
0503 }
0504
0505 if (theVerboseLevel > 9) {
0506 LogDebug("PixelCPETemplateReco") << " Sizex = " << theClusterParam.theCluster->sizeX()
0507 << " Sizey = " << theClusterParam.theCluster->sizeY() << " Edgex = " << edgex
0508 << " Edgey = " << edgey << " ErrX = " << xerr << " ErrY = " << yerr;
0509 }
0510
0511 }
0512
0513 if (!(xerr > 0.0f))
0514 throw cms::Exception("PixelCPETemplateReco::localError")
0515 << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
0516
0517 if (!(yerr > 0.0f))
0518 throw cms::Exception("PixelCPETemplateReco::localError")
0519 << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
0520
0521
0522
0523
0524
0525
0526
0527 return LocalError(xerr * xerr, 0, yerr * yerr);
0528 }
0529
0530 void PixelCPETemplateReco::fillPSetDescription(edm::ParameterSetDescription& desc) {
0531 desc.add<int>("barrelTemplateID", 0);
0532 desc.add<int>("forwardTemplateID", 0);
0533 desc.add<int>("directoryWithTemplates", 0);
0534 desc.add<int>("speed", -2);
0535 desc.add<bool>("UseClusterSplitter", false);
0536 }