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