File indexing completed on 2023-01-21 00:19:45
0001 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEGeneric.h"
0002
0003 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0004 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006
0007
0008 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h"
0009
0010
0011 #include "CondFormats/SiPixelTransient/interface/SiPixelUtils.h"
0012
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "MagneticField/Engine/interface/MagneticField.h"
0016
0017 #include "boost/multi_array.hpp"
0018
0019 #include <iostream>
0020 using namespace std;
0021
0022 namespace {
0023 constexpr float micronsToCm = 1.0e-4;
0024 }
0025
0026
0027
0028
0029 PixelCPEGeneric::PixelCPEGeneric(edm::ParameterSet const& conf,
0030 const MagneticField* mag,
0031 const TrackerGeometry& geom,
0032 const TrackerTopology& ttopo,
0033 const SiPixelLorentzAngle* lorentzAngle,
0034 const SiPixelGenErrorDBObject* genErrorDBObject,
0035 const SiPixelLorentzAngle* lorentzAngleWidth = nullptr)
0036 : PixelCPEGenericBase(conf, mag, geom, ttopo, lorentzAngle, genErrorDBObject, lorentzAngleWidth) {
0037 if (theVerboseLevel > 0)
0038 LogDebug("PixelCPEGeneric") << " constructing a generic algorithm for ideal pixel detector.\n"
0039 << " CPEGeneric:: VerboseLevel = " << theVerboseLevel;
0040
0041
0042 the_eff_charge_cut_lowX = conf.getParameter<double>("eff_charge_cut_lowX");
0043 the_eff_charge_cut_lowY = conf.getParameter<double>("eff_charge_cut_lowY");
0044 the_eff_charge_cut_highX = conf.getParameter<double>("eff_charge_cut_highX");
0045 the_eff_charge_cut_highY = conf.getParameter<double>("eff_charge_cut_highY");
0046 the_size_cutX = conf.getParameter<double>("size_cutX");
0047 the_size_cutY = conf.getParameter<double>("size_cutY");
0048
0049
0050 inflate_errors = conf.getParameter<bool>("inflate_errors");
0051 inflate_all_errors_no_trk_angle = conf.getParameter<bool>("inflate_all_errors_no_trk_angle");
0052
0053 NoTemplateErrorsWhenNoTrkAngles_ = conf.getParameter<bool>("NoTemplateErrorsWhenNoTrkAngles");
0054 IrradiationBiasCorrection_ = conf.getParameter<bool>("IrradiationBiasCorrection");
0055 DoCosmics_ = conf.getParameter<bool>("DoCosmics");
0056
0057 isPhase2_ = conf.getParameter<bool>("isPhase2");
0058
0059
0060 if ((DoCosmics_))
0061 useErrorsFromTemplates_ = false;
0062
0063 if (!useErrorsFromTemplates_ && (truncatePixelCharge_ || IrradiationBiasCorrection_ || LoadTemplatesFromDB_)) {
0064 throw cms::Exception("PixelCPEGeneric::PixelCPEGeneric: ")
0065 << "\nERROR: useErrorsFromTemplates_ is set to False in PixelCPEGeneric_cfi.py. "
0066 << " In this case it does not make sense to set any of the following to True: "
0067 << " truncatePixelCharge_, IrradiationBiasCorrection_, DoCosmics_, LoadTemplatesFromDB_ !!!"
0068 << "\n\n";
0069 }
0070
0071
0072 if (useErrorsFromTemplates_) {
0073 if (LoadTemplatesFromDB_) {
0074 if (!SiPixelGenError::pushfile(*genErrorDBObject_, thePixelGenError_))
0075 throw cms::Exception("InvalidCalibrationLoaded")
0076 << "ERROR: GenErrors not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
0077 << (*genErrorDBObject_).version();
0078 LogDebug("PixelCPEGeneric") << "Loaded genErrorDBObject v" << (*genErrorDBObject_).version();
0079 } else {
0080 if (!SiPixelGenError::pushfile(-999, thePixelGenError_))
0081 throw cms::Exception("InvalidCalibrationLoaded")
0082 << "ERROR: GenErrors not loaded correctly from text file. Reconstruction will fail.";
0083 }
0084
0085 } else {
0086 #ifdef EDM_ML_DEBUG
0087 cout << " Use simple parametrised errors " << endl;
0088 #endif
0089 }
0090
0091 #ifdef EDM_ML_DEBUG
0092 cout << "From PixelCPEGeneric::PixelCPEGeneric(...)" << endl;
0093 cout << "(int)useErrorsFromTemplates_ = " << (int)useErrorsFromTemplates_ << endl;
0094 cout << "truncatePixelCharge_ = " << (int)truncatePixelCharge_ << endl;
0095 cout << "IrradiationBiasCorrection_ = " << (int)IrradiationBiasCorrection_ << endl;
0096 cout << "(int)DoCosmics_ = " << (int)DoCosmics_ << endl;
0097 cout << "(int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl;
0098 #endif
0099 }
0100
0101
0102
0103
0104
0105
0106 LocalPoint PixelCPEGeneric::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0107 ClusterParamGeneric& theClusterParam = static_cast<ClusterParamGeneric&>(theClusterParamBase);
0108
0109
0110
0111 float chargeWidthX = (theDetParam.lorentzShiftInCmX * theDetParam.widthLAFractionX);
0112 float chargeWidthY = (theDetParam.lorentzShiftInCmY * theDetParam.widthLAFractionY);
0113 float shiftX = 0.5f * theDetParam.lorentzShiftInCmX;
0114 float shiftY = 0.5f * theDetParam.lorentzShiftInCmY;
0115
0116
0117
0118 if (useErrorsFromTemplates_) {
0119 float qclus = theClusterParam.theCluster->charge();
0120 float locBz = theDetParam.bz;
0121 float locBx = theDetParam.bx;
0122
0123
0124 theClusterParam.pixmx = -999;
0125 theClusterParam.sigmay = -999.9;
0126 theClusterParam.deltay = -999.9;
0127 theClusterParam.sigmax = -999.9;
0128 theClusterParam.deltax = -999.9;
0129 theClusterParam.sy1 = -999.9;
0130 theClusterParam.dy1 = -999.9;
0131 theClusterParam.sy2 = -999.9;
0132 theClusterParam.dy2 = -999.9;
0133 theClusterParam.sx1 = -999.9;
0134 theClusterParam.dx1 = -999.9;
0135 theClusterParam.sx2 = -999.9;
0136 theClusterParam.dx2 = -999.9;
0137
0138 SiPixelGenError gtempl(thePixelGenError_);
0139 int gtemplID_ = theDetParam.detTemplateId;
0140
0141
0142
0143
0144 theClusterParam.qBin_ = gtempl.qbin(gtemplID_,
0145 theClusterParam.cotalpha,
0146 theClusterParam.cotbeta,
0147 locBz,
0148 locBx,
0149 qclus,
0150 IrradiationBiasCorrection_,
0151 theClusterParam.pixmx,
0152 theClusterParam.sigmay,
0153 theClusterParam.deltay,
0154 theClusterParam.sigmax,
0155 theClusterParam.deltax,
0156 theClusterParam.sy1,
0157 theClusterParam.dy1,
0158 theClusterParam.sy2,
0159 theClusterParam.dy2,
0160 theClusterParam.sx1,
0161 theClusterParam.dx1,
0162 theClusterParam.sx2,
0163 theClusterParam.dx2);
0164
0165
0166
0167 bool useLAWidthFromGenError = false;
0168 if (useLAWidthFromGenError) {
0169 chargeWidthX = (-micronsToCm * gtempl.lorxwidth());
0170 chargeWidthY = (-micronsToCm * gtempl.lorywidth());
0171 LogDebug("PixelCPE localPosition():") << "redefine la width (gen-error)" << chargeWidthX << chargeWidthY;
0172 }
0173 LogDebug("PixelCPE localPosition():") << "GenError:" << gtemplID_;
0174
0175
0176 theClusterParam.deltax = theClusterParam.deltax * micronsToCm;
0177 theClusterParam.dx1 = theClusterParam.dx1 * micronsToCm;
0178 theClusterParam.dx2 = theClusterParam.dx2 * micronsToCm;
0179
0180 theClusterParam.deltay = theClusterParam.deltay * micronsToCm;
0181 theClusterParam.dy1 = theClusterParam.dy1 * micronsToCm;
0182 theClusterParam.dy2 = theClusterParam.dy2 * micronsToCm;
0183
0184 theClusterParam.sigmax = theClusterParam.sigmax * micronsToCm;
0185 theClusterParam.sx1 = theClusterParam.sx1 * micronsToCm;
0186 theClusterParam.sx2 = theClusterParam.sx2 * micronsToCm;
0187
0188 theClusterParam.sigmay = theClusterParam.sigmay * micronsToCm;
0189 theClusterParam.sy1 = theClusterParam.sy1 * micronsToCm;
0190 theClusterParam.sy2 = theClusterParam.sy2 * micronsToCm;
0191
0192 }
0193 else {
0194 theClusterParam.qBin_ = 0;
0195 }
0196
0197 int q_f_X;
0198 int q_l_X;
0199 int q_f_Y;
0200 int q_l_Y;
0201 collect_edge_charges(theClusterParam, q_f_X, q_l_X, q_f_Y, q_l_Y, useErrorsFromTemplates_ && truncatePixelCharge_);
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211 MeasurementPoint meas_URcorn_LLpix(theClusterParam.theCluster->minPixelRow() + 1.0,
0212 theClusterParam.theCluster->minPixelCol() + 1.0);
0213
0214
0215 MeasurementPoint meas_LLcorn_URpix(theClusterParam.theCluster->maxPixelRow(),
0216 theClusterParam.theCluster->maxPixelCol());
0217
0218
0219 LocalPoint local_URcorn_LLpix;
0220 LocalPoint local_LLcorn_URpix;
0221
0222
0223
0224 if (theClusterParam.with_track_angle) {
0225 local_URcorn_LLpix = theDetParam.theTopol->localPosition(meas_URcorn_LLpix, theClusterParam.loc_trk_pred);
0226 local_LLcorn_URpix = theDetParam.theTopol->localPosition(meas_LLcorn_URpix, theClusterParam.loc_trk_pred);
0227 } else {
0228 local_URcorn_LLpix = theDetParam.theTopol->localPosition(meas_URcorn_LLpix);
0229 local_LLcorn_URpix = theDetParam.theTopol->localPosition(meas_LLcorn_URpix);
0230 }
0231
0232 #ifdef EDM_ML_DEBUG
0233 if (theVerboseLevel > 20) {
0234 cout << "\n\t >>> theClusterParam.theCluster->x = " << theClusterParam.theCluster->x()
0235 << "\n\t >>> theClusterParam.theCluster->y = " << theClusterParam.theCluster->y()
0236 << "\n\t >>> cluster: minRow = " << theClusterParam.theCluster->minPixelRow()
0237 << " minCol = " << theClusterParam.theCluster->minPixelCol()
0238 << "\n\t >>> cluster: maxRow = " << theClusterParam.theCluster->maxPixelRow()
0239 << " maxCol = " << theClusterParam.theCluster->maxPixelCol()
0240 << "\n\t >>> meas: inner lower left = " << meas_URcorn_LLpix.x() << "," << meas_URcorn_LLpix.y()
0241 << "\n\t >>> meas: inner upper right = " << meas_LLcorn_URpix.x() << "," << meas_LLcorn_URpix.y() << endl;
0242 }
0243 #endif
0244
0245
0246
0247
0248
0249
0250
0251 #ifdef EDM_ML_DEBUG
0252 if (theVerboseLevel > 20)
0253 cout << "\t >>> Generic:: processing X" << endl;
0254 #endif
0255
0256 float xPos = siPixelUtils::generic_position_formula(
0257 theClusterParam.theCluster->sizeX(),
0258 q_f_X,
0259 q_l_X,
0260 local_URcorn_LLpix.x(),
0261 local_LLcorn_URpix.x(),
0262 chargeWidthX,
0263 theDetParam.theThickness,
0264 theClusterParam.cotalpha,
0265 theDetParam.thePitchX,
0266 theDetParam.theRecTopol->isItBigPixelInX(theClusterParam.theCluster->minPixelRow()),
0267 theDetParam.theRecTopol->isItBigPixelInX(theClusterParam.theCluster->maxPixelRow()),
0268 the_eff_charge_cut_lowX,
0269 the_eff_charge_cut_highX,
0270 the_size_cutX);
0271
0272
0273 xPos = xPos + shiftX;
0274
0275 #ifdef EDM_ML_DEBUG
0276 if (theVerboseLevel > 20)
0277 cout << "\t >>> Generic:: processing Y" << endl;
0278 #endif
0279
0280 float yPos = siPixelUtils::generic_position_formula(
0281 theClusterParam.theCluster->sizeY(),
0282 q_f_Y,
0283 q_l_Y,
0284 local_URcorn_LLpix.y(),
0285 local_LLcorn_URpix.y(),
0286 chargeWidthY,
0287 theDetParam.theThickness,
0288 theClusterParam.cotbeta,
0289 theDetParam.thePitchY,
0290 theDetParam.theRecTopol->isItBigPixelInY(theClusterParam.theCluster->minPixelCol()),
0291 theDetParam.theRecTopol->isItBigPixelInY(theClusterParam.theCluster->maxPixelCol()),
0292 the_eff_charge_cut_lowY,
0293 the_eff_charge_cut_highY,
0294 the_size_cutY);
0295
0296
0297 yPos = yPos + shiftY;
0298
0299
0300 if (IrradiationBiasCorrection_) {
0301 if (theClusterParam.theCluster->sizeX() == 1) {
0302
0303
0304
0305 xPos = xPos - (0.5f * theDetParam.lorentzShiftInCmX);
0306
0307 bool bigInX = theDetParam.theRecTopol->isItBigPixelInX(theClusterParam.theCluster->maxPixelRow());
0308 if (!bigInX)
0309 xPos -= theClusterParam.dx1;
0310 else
0311 xPos -= theClusterParam.dx2;
0312
0313 } else {
0314
0315 xPos -= theClusterParam.deltax;
0316
0317 }
0318
0319 if (theClusterParam.theCluster->sizeY() == 1) {
0320
0321 yPos = yPos - (0.5f * theDetParam.lorentzShiftInCmY);
0322
0323
0324 bool bigInY = theDetParam.theRecTopol->isItBigPixelInY(theClusterParam.theCluster->maxPixelCol());
0325 if (!bigInY)
0326 yPos -= theClusterParam.dy1;
0327 else
0328 yPos -= theClusterParam.dy2;
0329
0330 } else {
0331
0332 yPos -= theClusterParam.deltay;
0333 }
0334
0335 }
0336
0337
0338
0339
0340 LocalPoint pos_in_local(xPos, yPos);
0341 return pos_in_local;
0342 }
0343
0344
0345
0346
0347
0348
0349 LocalError PixelCPEGeneric::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0350 ClusterParamGeneric& theClusterParam = static_cast<ClusterParamGeneric&>(theClusterParamBase);
0351
0352
0353 float xerr, yerr;
0354 bool edgex, edgey, bigInX, bigInY;
0355 int maxPixelCol, maxPixelRow, minPixelCol, minPixelRow;
0356 uint sizex, sizey;
0357
0358 initializeLocalErrorVariables(xerr,
0359 yerr,
0360 edgex,
0361 edgey,
0362 bigInX,
0363 bigInY,
0364 maxPixelCol,
0365 maxPixelRow,
0366 minPixelCol,
0367 minPixelRow,
0368 sizex,
0369 sizey,
0370 theDetParam,
0371 theClusterParam);
0372
0373 bool useTempErrors =
0374 useErrorsFromTemplates_ && (!NoTemplateErrorsWhenNoTrkAngles_ || theClusterParam.with_track_angle);
0375
0376 if (int(sizex) != (maxPixelRow - minPixelRow + 1))
0377 LogDebug("PixelCPEGeneric") << " wrong x";
0378 if (int(sizey) != (maxPixelCol - minPixelCol + 1))
0379 LogDebug("PixelCPEGeneric") << " wrong y";
0380
0381 LogDebug("PixelCPEGeneric") << " edge clus " << xerr << " " << yerr;
0382 if (bigInX || bigInY)
0383 LogDebug("PixelCPEGeneric") << " big " << bigInX << " " << bigInY;
0384 if (edgex || edgey)
0385 LogDebug("PixelCPEGeneric") << " edge " << edgex << " " << edgey;
0386 LogDebug("PixelCPEGeneric") << " before if " << useErrorsFromTemplates_ << " " << theClusterParam.qBin_;
0387 if (theClusterParam.qBin_ == 0)
0388 LogDebug("PixelCPEGeneric") << " qbin 0! " << edgex << " " << edgey << " " << bigInX << " " << bigInY << " "
0389 << sizex << " " << sizey;
0390
0391
0392 setXYErrors(xerr, yerr, edgex, edgey, sizex, sizey, bigInX, bigInY, useTempErrors, theDetParam, theClusterParam);
0393
0394 if (!useTempErrors) {
0395 LogDebug("PixelCPEGeneric") << "Track angles are not known.\n"
0396 << "Default angle estimation which assumes track from PV (0,0,0) does not work.";
0397 }
0398
0399 if (!useTempErrors && inflate_errors) {
0400 int n_bigx = 0;
0401 int n_bigy = 0;
0402
0403 for (int irow = 0; irow < 7; ++irow) {
0404 if (theDetParam.theRecTopol->isItBigPixelInX(irow + minPixelRow))
0405 ++n_bigx;
0406 }
0407
0408 for (int icol = 0; icol < 21; ++icol) {
0409 if (theDetParam.theRecTopol->isItBigPixelInY(icol + minPixelCol))
0410 ++n_bigy;
0411 }
0412
0413 xerr = (float)(sizex + n_bigx) * theDetParam.thePitchX / std::sqrt(12.0f);
0414 yerr = (float)(sizey + n_bigy) * theDetParam.thePitchY / std::sqrt(12.0f);
0415 }
0416
0417 #ifdef EDM_ML_DEBUG
0418 if (!(xerr > 0.0))
0419 throw cms::Exception("PixelCPEGeneric::localError") << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
0420
0421 if (!(yerr > 0.0))
0422 throw cms::Exception("PixelCPEGeneric::localError") << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
0423 #endif
0424
0425 LogDebug("PixelCPEGeneric") << " errors " << xerr << " " << yerr;
0426 if (theClusterParam.qBin_ == 0)
0427 LogDebug("PixelCPEGeneric") << " qbin 0 " << xerr << " " << yerr;
0428
0429 auto xerr_sq = xerr * xerr;
0430 auto yerr_sq = yerr * yerr;
0431
0432 return LocalError(xerr_sq, 0, yerr_sq);
0433 }
0434
0435 void PixelCPEGeneric::fillPSetDescription(edm::ParameterSetDescription& desc) {
0436 PixelCPEGenericBase::fillPSetDescription(desc);
0437 desc.add<double>("eff_charge_cut_highX", 1.0);
0438 desc.add<double>("eff_charge_cut_highY", 1.0);
0439 desc.add<double>("eff_charge_cut_lowX", 0.0);
0440 desc.add<double>("eff_charge_cut_lowY", 0.0);
0441 desc.add<double>("size_cutX", 3.0);
0442 desc.add<double>("size_cutY", 3.0);
0443 desc.add<double>("EdgeClusterErrorX", 50.0);
0444 desc.add<double>("EdgeClusterErrorY", 85.0);
0445 desc.add<bool>("inflate_errors", false);
0446 desc.add<bool>("inflate_all_errors_no_trk_angle", false);
0447 desc.add<bool>("NoTemplateErrorsWhenNoTrkAngles", false);
0448 desc.add<bool>("UseErrorsFromTemplates", true);
0449 desc.add<bool>("TruncatePixelCharge", true);
0450 desc.add<bool>("IrradiationBiasCorrection", false);
0451 desc.add<bool>("DoCosmics", false);
0452 desc.add<bool>("isPhase2", false);
0453 desc.add<bool>("SmallPitch", false);
0454 }