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