File indexing completed on 2025-02-05 23:51:49
0001 #ifndef SimTracker_Common_SiPixelChargeReweightingAlgorithm_h
0002 #define SimTracker_Common_SiPixelChargeReweightingAlgorithm_h
0003
0004
0005
0006
0007
0008 #include <map>
0009 #include <memory>
0010 #include <vector>
0011 #include <iostream>
0012 #include <iomanip>
0013 #include <type_traits>
0014 #include <gsl/gsl_sf_erf.h>
0015
0016 #include "boost/multi_array.hpp"
0017
0018 #include "CondFormats/DataRecord/interface/SiPixel2DTemplateDBObjectRcd.h"
0019 #include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
0020 #include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
0021 #include "CondFormats/SiPixelObjects/interface/SiPixel2DTemplateDBObject.h"
0022 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate2D.h"
0023 #include "DataFormats/DetId/interface/DetId.h"
0024 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0025 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0026 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0027 #include "FWCore/Framework/interface/ConsumesCollector.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "FWCore/Utilities/interface/Exception.h"
0032 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0033 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0034 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0035 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0036 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0037 #include "SimDataFormats/Track/interface/SimTrack.h"
0038 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0039 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0040 #include "SimTracker/Common/interface/DigitizerUtility.h"
0041 #include "SimTracker/SiPixelDigitizer/plugins/SiPixelDigitizerAlgorithm.h"
0042
0043
0044 class DetId;
0045 class GaussianTailNoiseGenerator;
0046 class PixelDigi;
0047 class PixelDigiSimLink;
0048 class PixelGeomDetUnit;
0049 class SiG4UniversalFluctuation;
0050
0051 class SiPixelChargeReweightingAlgorithm {
0052 public:
0053 SiPixelChargeReweightingAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC);
0054 ~SiPixelChargeReweightingAlgorithm();
0055
0056
0057 void init(const edm::EventSetup& es);
0058
0059
0060 template <class AmplitudeType, typename SignalType>
0061 bool hitSignalReweight(const PSimHit& hit,
0062 std::map<int, float, std::less<int> >& hit_signal,
0063 const size_t hitIndex,
0064 const size_t hitIndex4CR,
0065 const unsigned int tofBin,
0066 const PixelTopology* topol,
0067 uint32_t detID,
0068 SignalType& theSignal,
0069 unsigned short int processType,
0070 const bool& boolmakeDigiSimLinks);
0071
0072 template <class AmplitudeType, typename SignalType>
0073 bool lateSignalReweight(const PixelGeomDetUnit* pixdet,
0074 std::vector<PixelDigi>& digis,
0075 PixelSimHitExtraInfo& loopTempSH,
0076 SignalType& theNewDigiSignal,
0077 const TrackerTopology* tTopo,
0078 CLHEP::HepRandomEngine* engine);
0079
0080 template <class AmplitudeType, typename SignalType>
0081 bool lateSignalReweight(const PixelGeomDetUnit* pixdet,
0082 std::vector<PixelDigi>& digis,
0083 PixelSimHitExtraInfoLite& loopTempSH,
0084 SignalType& theNewDigiSignal,
0085 const TrackerTopology* tTopo,
0086 CLHEP::HepRandomEngine* engine);
0087 private:
0088
0089 typedef GloballyPositioned<double> Frame;
0090 typedef std::vector<edm::ParameterSet> Parameters;
0091 typedef boost::multi_array<float, 2> array_2d;
0092
0093 std::vector<SiPixelTemplateStore2D> templateStores_;
0094
0095
0096 SiPixelTemplate2D templ2D;
0097 std::vector<bool> xdouble;
0098 std::vector<bool> ydouble;
0099 std::vector<float> track;
0100 int IDnum, IDden;
0101
0102 const bool UseReweighting;
0103 bool applyLateReweighting_;
0104 const bool PrintClusters;
0105 const bool PrintTemplates;
0106
0107 static constexpr float cmToMicrons = 10000.f;
0108
0109 edm::ESGetToken<SiPixel2DTemplateDBObject, SiPixel2DTemplateDBObjectRcd> SiPixel2DTemp_den_token_;
0110 edm::ESGetToken<SiPixel2DTemplateDBObject, SiPixel2DTemplateDBObjectRcd> SiPixel2DTemp_num_token_;
0111 const SiPixel2DTemplateDBObject* dbobject_den;
0112 const SiPixel2DTemplateDBObject* dbobject_num;
0113
0114
0115 int PixelTempRewgt2D(int id_gen, int id_rewgt, array_2d& cluster);
0116 void printCluster(array_2d& cluster);
0117 void printCluster(float arr[BXM2][BYM2]);
0118 void printCluster(float arr[TXSIZE][TYSIZE]);
0119 };
0120
0121 inline void SiPixelChargeReweightingAlgorithm::init(const edm::EventSetup& es) {
0122
0123 if (UseReweighting || applyLateReweighting_) {
0124 dbobject_den = &es.getData(SiPixel2DTemp_den_token_);
0125 dbobject_num = &es.getData(SiPixel2DTemp_num_token_);
0126
0127 int numOfTemplates = dbobject_den->numOfTempl() + dbobject_num->numOfTempl();
0128 templateStores_.reserve(numOfTemplates);
0129 SiPixelTemplate2D::pushfile(*dbobject_den, templateStores_);
0130 SiPixelTemplate2D::pushfile(*dbobject_num, templateStores_);
0131
0132 track.resize(6);
0133 LogDebug("SiPixelChargeReweightingAlgorithm") << "Init SiPixelChargeReweightingAlgorithm";
0134 }
0135 }
0136
0137
0138
0139 inline SiPixelChargeReweightingAlgorithm::SiPixelChargeReweightingAlgorithm(const edm::ParameterSet& conf,
0140 edm::ConsumesCollector iC)
0141 :
0142
0143 templ2D(templateStores_),
0144 xdouble(TXSIZE),
0145 ydouble(TYSIZE),
0146 IDnum(conf.exists("TemplateIDnumerator") ? conf.getParameter<int>("TemplateIDnumerator") : 0),
0147 IDden(conf.exists("TemplateIDdenominator") ? conf.getParameter<int>("TemplateIDdenominator") : 0),
0148
0149 UseReweighting(conf.getParameter<bool>("UseReweighting")),
0150 applyLateReweighting_(conf.exists("applyLateReweighting") ? conf.getParameter<bool>("applyLateReweighting")
0151 : false),
0152 PrintClusters(conf.exists("PrintClusters") ? conf.getParameter<bool>("PrintClusters") : false),
0153 PrintTemplates(conf.exists("PrintTemplates") ? conf.getParameter<bool>("PrintTemplates") : false) {
0154 if (UseReweighting || applyLateReweighting_) {
0155 SiPixel2DTemp_den_token_ = iC.esConsumes(edm::ESInputTag("", "denominator"));
0156 SiPixel2DTemp_num_token_ = iC.esConsumes(edm::ESInputTag("", "numerator"));
0157 }
0158 LogDebug("SiPixelChargeReweightingAlgorithm")
0159 << "SiPixelChargeReweightingAlgorithm constructed"
0160 << " with UseReweighting = " << UseReweighting << " and applyLateReweighting_ = " << applyLateReweighting_;
0161 }
0162
0163
0164 inline SiPixelChargeReweightingAlgorithm::~SiPixelChargeReweightingAlgorithm() {
0165 LogDebug("SiPixelChargeReweightingAlgorithm") << "SiPixelChargeReweightingAlgorithm deleted";
0166 }
0167
0168
0169 template <class AmplitudeType, typename SignalType>
0170 bool SiPixelChargeReweightingAlgorithm::hitSignalReweight(const PSimHit& hit,
0171 std::map<int, float, std::less<int> >& hit_signal,
0172 const size_t hitIndex,
0173 const size_t hitIndex4CR,
0174 const unsigned int tofBin,
0175 const PixelTopology* topol,
0176 uint32_t detID,
0177 SignalType& theSignal,
0178 unsigned short int processType,
0179 const bool& boolmakeDigiSimLinks) {
0180 int irow_min = topol->nrows();
0181 int irow_max = 0;
0182 int icol_min = topol->ncolumns();
0183 int icol_max = 0;
0184
0185 float chargeBefore = 0;
0186 float chargeAfter = 0;
0187 SignalType hitSignal;
0188 LocalVector direction = hit.exitPoint() - hit.entryPoint();
0189
0190 for (std::map<int, float, std::less<int> >::const_iterator im = hit_signal.begin(); im != hit_signal.end(); ++im) {
0191 int chan = (*im).first;
0192 std::pair<int, int> pixelWithCharge = PixelDigi::channelToPixel(chan);
0193 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0194 hitSignal[chan] += AmplitudeType((*im).second, &hit, (*im).second, 0., hitIndex, tofBin);
0195 } else {
0196 hitSignal[chan] +=
0197 (boolmakeDigiSimLinks ? AmplitudeType((*im).second, &hit, hitIndex, hitIndex4CR, tofBin, (*im).second)
0198 : AmplitudeType((*im).second, (*im).second));
0199 }
0200 chargeBefore += (*im).second;
0201
0202 if (pixelWithCharge.first < irow_min)
0203 irow_min = pixelWithCharge.first;
0204 if (pixelWithCharge.first > irow_max)
0205 irow_max = pixelWithCharge.first;
0206 if (pixelWithCharge.second < icol_min)
0207 icol_min = pixelWithCharge.second;
0208 if (pixelWithCharge.second > icol_max)
0209 icol_max = pixelWithCharge.second;
0210 }
0211
0212 LocalPoint hitEntryPoint = hit.entryPoint();
0213
0214 float trajectoryScaleToPosition = std::abs(hitEntryPoint.z() / direction.z());
0215
0216 LocalPoint hitPosition = hitEntryPoint + trajectoryScaleToPosition * direction;
0217
0218 MeasurementPoint hitPositionPixel = topol->measurementPosition(hit.localPosition());
0219 std::pair<int, int> hitPixel =
0220 std::pair<int, int>(int(floor(hitPositionPixel.x())), int(floor(hitPositionPixel.y())));
0221
0222 MeasurementPoint originPixel = MeasurementPoint(hitPixel.first - THX + 0.5, hitPixel.second - THY + 0.5);
0223 LocalPoint origin = topol->localPosition(originPixel);
0224
0225 MeasurementPoint hitEntryPointPixel = topol->measurementPosition(hit.entryPoint());
0226 MeasurementPoint hitExitPointPixel = topol->measurementPosition(hit.exitPoint());
0227 std::pair<int, int> entryPixel =
0228 std::pair<int, int>(int(floor(hitEntryPointPixel.x())), int(floor(hitEntryPointPixel.y())));
0229 std::pair<int, int> exitPixel =
0230 std::pair<int, int>(int(floor(hitExitPointPixel.x())), int(floor(hitExitPointPixel.y())));
0231
0232 int hitcol_min, hitcol_max, hitrow_min, hitrow_max;
0233 if (entryPixel.first > exitPixel.first) {
0234 hitrow_min = exitPixel.first;
0235 hitrow_max = entryPixel.first;
0236 } else {
0237 hitrow_min = entryPixel.first;
0238 hitrow_max = exitPixel.first;
0239 }
0240
0241 if (entryPixel.second > exitPixel.second) {
0242 hitcol_min = exitPixel.second;
0243 hitcol_max = entryPixel.second;
0244 } else {
0245 hitcol_min = entryPixel.second;
0246 hitcol_max = exitPixel.second;
0247 }
0248
0249 LocalPoint CMSSWhitPosition = hit.localPosition();
0250 LogDebug("SiPixelChargeReweightingAlgorithm")
0251 << "\n"
0252 << "Particle ID is: " << hit.particleType() << "\n"
0253 << "Process type: " << hit.processType() << "\n"
0254 << "HitPosition:"
0255 << "\n"
0256 << "Hit entry x/y/z: " << hit.entryPoint().x() << " " << hit.entryPoint().y() << " " << hit.entryPoint().z()
0257 << " "
0258 << "Hit exit x/y/z: " << hit.exitPoint().x() << " " << hit.exitPoint().y() << " " << hit.exitPoint().z() << " "
0259
0260 << "Pixel Pos - X: " << hitPositionPixel.x() << " Y: " << hitPositionPixel.y() << "\n"
0261 << "Cart.Cor. - X: " << CMSSWhitPosition.x() << " Y: " << CMSSWhitPosition.y() << "\n"
0262 << "Z=0 Pos - X: " << hitPosition.x() << " Y: " << hitPosition.y() << "\n"
0263
0264 << "Origin of the template:"
0265 << "\n"
0266 << "Pixel Pos - X: " << originPixel.x() << " Y: " << originPixel.y() << "\n"
0267 << "Cart.Cor. - X: " << origin.x() << " Y: " << origin.y() << "\n"
0268 << "\n"
0269 << "Entry/Exit:"
0270 << "\n"
0271 << "Entry - X: " << hit.entryPoint().x() << " Y: " << hit.entryPoint().y() << " Z: " << hit.entryPoint().z()
0272 << "\n"
0273 << "Exit - X: " << hit.exitPoint().x() << " Y: " << hit.exitPoint().y() << " Z: " << hit.exitPoint().z() << "\n"
0274
0275 << "Entry - X Pixel: " << hitEntryPointPixel.x() << " Y Pixel: " << hitEntryPointPixel.y() << "\n"
0276 << "Exit - X Pixel: " << hitExitPointPixel.x() << " Y Pixel: " << hitExitPointPixel.y() << "\n"
0277
0278 << "row min: " << irow_min << " col min: " << icol_min << "\n";
0279
0280 if (!(irow_min <= hitrow_max && irow_max >= hitrow_min && icol_min <= hitcol_max && icol_max >= hitcol_min)) {
0281
0282 LogDebug("SiPixelChargeReweightingAlgorithm") << "Outside the row/col boundary, exit charge reweighting";
0283 return false;
0284 }
0285
0286 track[0] = (hitPosition.x() - origin.x()) * cmToMicrons;
0287 track[1] = (hitPosition.y() - origin.y()) * cmToMicrons;
0288 track[2] = 0.0f;
0289 track[3] = direction.x();
0290 track[4] = direction.y();
0291 track[5] = direction.z();
0292
0293 LogDebug("SiPixelChargeReweightingAlgorithm") << "Init array of size x = " << TXSIZE << " and y = " << TYSIZE;
0294 array_2d pixrewgt(boost::extents[TXSIZE][TYSIZE]);
0295
0296 for (int row = 0; row < TXSIZE; ++row) {
0297 for (int col = 0; col < TYSIZE; ++col) {
0298 pixrewgt[row][col] = 0;
0299 }
0300 }
0301
0302 for (int row = 0; row < TXSIZE; ++row) {
0303 xdouble[row] = topol->isItBigPixelInX(hitPixel.first + row - THX);
0304 }
0305
0306 for (int col = 0; col < TYSIZE; ++col) {
0307 ydouble[col] = topol->isItBigPixelInY(hitPixel.second + col - THY);
0308 }
0309
0310
0311
0312 int rowmin = std::max(0, THX - hitPixel.first);
0313 int rowmax = std::min(TXSIZE, topol->nrows() + THX - hitPixel.first);
0314 int colmin = std::max(0, THY - hitPixel.second);
0315 int colmax = std::min(TYSIZE, topol->ncolumns() + THY - hitPixel.second);
0316
0317 for (int row = rowmin; row < rowmax; ++row) {
0318 for (int col = colmin; col < colmax; ++col) {
0319
0320 pixrewgt[row][col] =
0321 hitSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)];
0322 }
0323 }
0324
0325 if (PrintClusters) {
0326 LogDebug("SiPixelChargeReweightingAlgorithm ") << "Cluster before reweighting: ";
0327 printCluster(pixrewgt);
0328 }
0329
0330 int ierr;
0331
0332
0333 if (UseReweighting == true) {
0334 int ID1 = dbobject_num->getTemplateID(detID);
0335 int ID0 = dbobject_den->getTemplateID(detID);
0336
0337 if (ID0 == ID1) {
0338 return false;
0339 }
0340 ierr = PixelTempRewgt2D(ID0, ID1, pixrewgt);
0341 } else {
0342 ierr = PixelTempRewgt2D(IDden, IDden, pixrewgt);
0343 }
0344 if (ierr != 0) {
0345 LogDebug("SiPixelChargeReweightingAlgorithm ") << "Cluster Charge Reweighting did not work properly.";
0346 return false;
0347 }
0348
0349 if (PrintClusters) {
0350 LogDebug("SiPixelChargeReweightingAlgorithm ") << "Cluster after reweighting: ";
0351 printCluster(pixrewgt);
0352 }
0353
0354 for (int row = rowmin; row < rowmax; ++row) {
0355 for (int col = colmin; col < colmax; ++col) {
0356 float charge = pixrewgt[row][col];
0357 if (charge > 0) {
0358 chargeAfter += charge;
0359 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0360 theSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
0361 AmplitudeType(charge, &hit, charge, 0., hitIndex, tofBin);
0362 } else {
0363 theSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
0364 (boolmakeDigiSimLinks ? AmplitudeType(charge, &hit, hitIndex, hitIndex4CR, tofBin, charge)
0365 : AmplitudeType(charge, charge));
0366 }
0367 }
0368 }
0369 }
0370
0371 if (chargeBefore != 0. && chargeAfter == 0.) {
0372 return false;
0373 }
0374
0375 if (PrintClusters) {
0376 LogDebug("SiPixelChargeReweightingAlgorithm ")
0377 << "Charges (before->after): " << chargeBefore << " -> " << chargeAfter;
0378 LogDebug("SiPixelChargeReweightingAlgorithm ")
0379 << "Charge loss: " << (1 - chargeAfter / chargeBefore) * 100 << " % \n";
0380 }
0381
0382 return true;
0383 }
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394 inline int SiPixelChargeReweightingAlgorithm::PixelTempRewgt2D(int id_in, int id_rewgt, array_2d& cluster) {
0395
0396 int i, j, k, l, kclose;
0397 int nclusx, nclusy, success;
0398 float xsize, ysize, q50i, q100i, q50r, q10r, xhit2D, yhit2D, dist2, dmin2;
0399 float qscalei, rqscale;
0400 float xy_in[BXM2][BYM2], xy_rewgt[BXM2][BYM2], xy_clust[TXSIZE][TYSIZE];
0401 int denx_clust[TXSIZE][TYSIZE], deny_clust[TXSIZE][TYSIZE];
0402 float cotalpha, cotbeta;
0403
0404 success = 0;
0405
0406
0407 array_2d clust(cluster);
0408
0409
0410 templ2D.getid(id_in);
0411 xsize = templ2D.xsize();
0412 ysize = templ2D.ysize();
0413
0414
0415
0416 if (std::abs(track[5]) > 0.f) {
0417 cotalpha = track[3] / track[5];
0418 cotbeta = track[4] / track[5];
0419 } else {
0420 LogDebug("SiPixelChargeReweightingAlgorithm") << "Reweighting angle is not good! \n";
0421 return 9;
0422 }
0423
0424
0425 if (ydouble[0]) {
0426 yhit2D = track[1] - cotbeta * track[2] + ysize;
0427 } else {
0428 yhit2D = track[1] - cotbeta * track[2] + 0.5f * ysize;
0429 }
0430 if (xdouble[0]) {
0431 xhit2D = track[0] - cotalpha * track[2] + xsize;
0432 } else {
0433 xhit2D = track[0] - cotalpha * track[2] + 0.5f * xsize;
0434 }
0435
0436
0437 for (i = 0; i < BYM2; ++i) {
0438 for (j = 0; j < BXM2; ++j) {
0439 xy_in[j][i] = 0.f;
0440 xy_rewgt[j][i] = 0.f;
0441 }
0442 }
0443
0444
0445
0446 if (!templ2D.xytemp(id_in, cotalpha, cotbeta, xhit2D, yhit2D, ydouble, xdouble, xy_in)) {
0447 success = 1;
0448 }
0449 if (success != 0) {
0450 LogDebug("SiPixelChargeReweightingAlgorithm") << "No matching template found \n";
0451 return 2;
0452 }
0453
0454 if (PrintTemplates) {
0455 LogDebug("SiPixelChargeReweightingAlgorithm") << "Template unirrad: \n";
0456 printCluster(xy_in);
0457 }
0458
0459 q50i = templ2D.s50();
0460
0461 q100i = 2.f * q50i;
0462
0463
0464 qscalei = templ2D.qscale();
0465
0466
0467
0468 if (cluster.num_dimensions() != 2) {
0469 edm::LogWarning("SiPixelChargeReweightingAlgorithm") << "Cluster is not 2-dimensional. Return. \n";
0470 return 3;
0471 }
0472 nclusx = (int)cluster.shape()[0];
0473 nclusy = (int)cluster.shape()[1];
0474 if (nclusx != TXSIZE || xdouble.size() != TXSIZE) {
0475 edm::LogWarning("SiPixelChargeReweightingAlgorithm")
0476 << "Sizes in x do not match: nclusx=" << nclusx << " xdoubleSize=" << xdouble.size() << " TXSIZE=" << TXSIZE
0477 << ". Return. \n";
0478 return 4;
0479 }
0480 if (nclusy != TYSIZE || ydouble.size() != TYSIZE) {
0481 edm::LogWarning("SiPixelChargeReweightingAlgorithm") << "Sizes in y do not match. Return. \n";
0482 return 5;
0483 }
0484
0485
0486
0487 for (i = 0; i < TYSIZE; ++i) {
0488 for (j = 0; j < TXSIZE; ++j) {
0489 xy_clust[j][i] = 0.f;
0490 denx_clust[j][i] = 0;
0491 deny_clust[j][i] = 0;
0492 }
0493 }
0494
0495
0496
0497 if (!templ2D.xytemp(id_rewgt, cotalpha, cotbeta, xhit2D, yhit2D, ydouble, xdouble, xy_rewgt)) {
0498 success = 1;
0499 }
0500
0501 if (PrintTemplates) {
0502 LogDebug("SiPixelChargeReweightingAlgorithm") << "Template irrad: \n";
0503 printCluster(xy_rewgt);
0504 }
0505
0506 q50r = templ2D.s50();
0507 q10r = 0.2f * q50r;
0508
0509
0510 rqscale = qscalei / templ2D.qscale();
0511
0512
0513
0514 int ntpix = 0;
0515 int ncpix = 0;
0516 std::vector<int> ytclust;
0517 std::vector<int> xtclust;
0518 std::vector<int> ycclust;
0519 std::vector<int> xcclust;
0520 for (i = 0; i < TYSIZE; ++i) {
0521 for (j = 0; j < TXSIZE; ++j) {
0522 if (xy_in[j + 1][i + 1] > q100i) {
0523 ++ntpix;
0524 ytclust.push_back(i);
0525 xtclust.push_back(j);
0526 xy_clust[j][i] = xy_rewgt[j + 1][i + 1] / xy_in[j + 1][i + 1];
0527 denx_clust[j][i] = j;
0528 deny_clust[j][i] = i;
0529 }
0530 }
0531 }
0532
0533
0534
0535 for (i = 0; i < TYSIZE; ++i) {
0536 for (j = 0; j < TXSIZE; ++j) {
0537 if (xy_rewgt[j + 1][i + 1] > q10r && xy_clust[j][i] == 0.f && ntpix > 0) {
0538
0539 dmin2 = 10000.f;
0540 kclose = 0;
0541 for (k = 0; k < ntpix; ++k) {
0542 dist2 = (i - ytclust[k]) * (i - ytclust[k]) + 0.44444f * (j - xtclust[k]) * (j - xtclust[k]);
0543 if (dist2 < dmin2) {
0544 dmin2 = dist2;
0545 kclose = k;
0546 }
0547 }
0548 xy_clust[j][i] = xy_rewgt[j + 1][i + 1] / xy_in[xtclust[kclose] + 1][ytclust[kclose] + 1];
0549 denx_clust[j][i] = xtclust[kclose];
0550 deny_clust[j][i] = ytclust[kclose];
0551 }
0552 }
0553 }
0554
0555 if (PrintTemplates) {
0556 LogDebug("SiPixelChargeReweightingAlgorithm") << "Weights: \n";
0557 printCluster(xy_clust);
0558 }
0559
0560 for (i = 0; i < TYSIZE; ++i) {
0561 for (j = 0; j < TXSIZE; ++j) {
0562 if (xy_clust[j][i] > 0.f) {
0563 cluster[j][i] = xy_clust[j][i] * clust[denx_clust[j][i]][deny_clust[j][i]];
0564 } else {
0565 if (clust[j][i] > 0.f) {
0566 ++ncpix;
0567 ycclust.push_back(i);
0568 xcclust.push_back(j);
0569 }
0570 }
0571 }
0572 }
0573
0574
0575
0576 if (ncpix > 0) {
0577 for (l = 0; l < ncpix; ++l) {
0578 i = ycclust[l];
0579 j = xcclust[l];
0580 dmin2 = 10000.f;
0581 kclose = 0;
0582 for (k = 0; k < ntpix; ++k) {
0583 dist2 = (i - ytclust[k]) * (i - ytclust[k]) + 0.44444f * (j - xtclust[k]) * (j - xtclust[k]);
0584 if (dist2 < dmin2) {
0585 dmin2 = dist2;
0586 kclose = k;
0587 }
0588 }
0589 if (dmin2 < 5.f) {
0590 cluster[j][i] *= xy_clust[xtclust[kclose]][ytclust[kclose]];
0591 } else {
0592 cluster[j][i] = 0.f;
0593 }
0594 }
0595 }
0596
0597
0598
0599 for (i = 0; i < TYSIZE; ++i) {
0600 for (j = 0; j < TXSIZE; ++j) {
0601 cluster[j][i] *= rqscale;
0602 }
0603 }
0604
0605 return success;
0606 }
0607
0608 inline void SiPixelChargeReweightingAlgorithm::printCluster(array_2d& cluster) {
0609 for (int col = 0; col < TYSIZE; ++col) {
0610 for (int row = 0; row < TXSIZE; ++row) {
0611 LogDebug("SiPixelChargeReweightingAlgorithm") << cluster[row][col];
0612 }
0613 LogDebug("SiPixelChargeReweightingAlgorithm") << "\n";
0614 }
0615 }
0616
0617 inline void SiPixelChargeReweightingAlgorithm::printCluster(float arr[BXM2][BYM2]) {
0618 for (int col = 0; col < BYM2; ++col) {
0619 for (int row = 0; row < BXM2; ++row) {
0620 LogDebug("SiPixelChargeReweightingAlgorithm") << arr[row][col];
0621 }
0622 LogDebug("SiPixelChargeReweightingAlgorithm") << "\n";
0623 }
0624 }
0625
0626 inline void SiPixelChargeReweightingAlgorithm::printCluster(float arr[TXSIZE][TYSIZE]) {
0627 for (int col = 0; col < TYSIZE; ++col) {
0628 for (int row = 0; row < TXSIZE; ++row) {
0629 LogDebug("SiPixelChargeReweightingAlgorithm") << arr[row][col];
0630 }
0631 LogDebug("SiPixelChargeReweightingAlgorithm") << "\n";
0632 }
0633 }
0634
0635 template <class AmplitudeType, typename SignalType>
0636 bool SiPixelChargeReweightingAlgorithm::lateSignalReweight(const PixelGeomDetUnit* pixdet,
0637 std::vector<PixelDigi>& digis,
0638 PixelSimHitExtraInfo& loopTempSH,
0639 SignalType& theNewDigiSignal,
0640 const TrackerTopology* tTopo,
0641 CLHEP::HepRandomEngine* engine) {
0642 uint32_t detID = pixdet->geographicalId().rawId();
0643 const PixelTopology* topol = &pixdet->specificTopology();
0644
0645 if (UseReweighting) {
0646 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0647 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0648 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0649 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ***** INCONSISTENCY !!! ***** \n";
0650 edm::LogError("SiPixelChargeReweightingAlgorithm")
0651 << " applyLateReweighting_ and UseReweighting can not be true at the same time for PU ! \n";
0652 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ---> DO NOT APPLY CHARGE REWEIGHTING TWICE !!! \n";
0653 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0654 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0655 return false;
0656 }
0657
0658 SignalType theDigiSignal;
0659
0660 int irow_min = topol->nrows();
0661 int irow_max = 0;
0662 int icol_min = topol->ncolumns();
0663 int icol_max = 0;
0664
0665 float chargeBefore = 0;
0666 float chargeAfter = 0;
0667
0668
0669 std::vector<PixelDigi>::const_iterator loopDigi;
0670 for (loopDigi = digis.begin(); loopDigi != digis.end(); ++loopDigi) {
0671 unsigned int chan = loopDigi->channel();
0672 if (loopTempSH.isInTheList(chan)) {
0673 float corresponding_charge = loopDigi->adc();
0674 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0675 edm::LogError("SiPixelChargeReweightingAlgorithm")
0676 << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
0677 return false;
0678 } else {
0679 theDigiSignal[chan] += AmplitudeType(corresponding_charge, corresponding_charge);
0680 }
0681 chargeBefore += corresponding_charge;
0682 if (loopDigi->row() < irow_min)
0683 irow_min = loopDigi->row();
0684 if (loopDigi->row() > irow_max)
0685 irow_max = loopDigi->row();
0686 if (loopDigi->column() < icol_min)
0687 icol_min = loopDigi->column();
0688 if (loopDigi->column() > icol_max)
0689 icol_max = loopDigi->column();
0690 }
0691 }
0692
0693
0694 LocalVector direction = loopTempSH.exitPoint() - loopTempSH.entryPoint();
0695 LocalPoint hitEntryPoint = loopTempSH.entryPoint();
0696 float trajectoryScaleToPosition = std::abs(hitEntryPoint.z() / direction.z());
0697 LocalPoint hitPosition = hitEntryPoint + trajectoryScaleToPosition * direction;
0698
0699
0700
0701 LocalPoint hitLocalPosition = hitEntryPoint + 0.5 * direction;
0702 MeasurementPoint hitPositionPixel = topol->measurementPosition(hitLocalPosition);
0703 std::pair<int, int> hitPixel =
0704 std::pair<int, int>(int(floor(hitPositionPixel.x())), int(floor(hitPositionPixel.y())));
0705
0706 MeasurementPoint originPixel = MeasurementPoint(hitPixel.first - THX + 0.5, hitPixel.second - THY + 0.5);
0707 LocalPoint origin = topol->localPosition(originPixel);
0708
0709 MeasurementPoint hitEntryPointPixel = topol->measurementPosition(loopTempSH.entryPoint());
0710 MeasurementPoint hitExitPointPixel = topol->measurementPosition(loopTempSH.exitPoint());
0711 std::pair<int, int> entryPixel =
0712 std::pair<int, int>(int(floor(hitEntryPointPixel.x())), int(floor(hitEntryPointPixel.y())));
0713 std::pair<int, int> exitPixel =
0714 std::pair<int, int>(int(floor(hitExitPointPixel.x())), int(floor(hitExitPointPixel.y())));
0715
0716 int hitcol_min, hitcol_max, hitrow_min, hitrow_max;
0717 if (entryPixel.first > exitPixel.first) {
0718 hitrow_min = exitPixel.first;
0719 hitrow_max = entryPixel.first;
0720 } else {
0721 hitrow_min = entryPixel.first;
0722 hitrow_max = exitPixel.first;
0723 }
0724
0725 if (entryPixel.second > exitPixel.second) {
0726 hitcol_min = exitPixel.second;
0727 hitcol_max = entryPixel.second;
0728 } else {
0729 hitcol_min = entryPixel.second;
0730 hitcol_max = exitPixel.second;
0731 }
0732
0733 if (!(irow_min <= hitrow_max && irow_max >= hitrow_min && icol_min <= hitcol_max && icol_max >= hitcol_min)) {
0734
0735 return false;
0736 }
0737
0738 track[0] = (hitPosition.x() - origin.x()) * cmToMicrons;
0739 track[1] = (hitPosition.y() - origin.y()) * cmToMicrons;
0740 track[2] = 0.0f;
0741 track[3] = direction.x();
0742 track[4] = direction.y();
0743 track[5] = direction.z();
0744
0745 array_2d pixrewgt(boost::extents[TXSIZE][TYSIZE]);
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755 for (int row = 0; row < TXSIZE; ++row) {
0756 xdouble[row] = topol->isItBigPixelInX(hitPixel.first + row - THX);
0757 }
0758
0759 for (int col = 0; col < TYSIZE; ++col) {
0760 ydouble[col] = topol->isItBigPixelInY(hitPixel.second + col - THY);
0761 }
0762
0763 for (int row = 0; row < TXSIZE; ++row) {
0764 for (int col = 0; col < TYSIZE; ++col) {
0765
0766 pixrewgt[row][col] =
0767 theDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)];
0768 }
0769 }
0770
0771 if (PrintClusters) {
0772 LogDebug("SiPixelChargeReweightingAlgorithm") << " Cluster before reweighting: ";
0773 printCluster(pixrewgt);
0774 }
0775
0776 int ierr;
0777
0778
0779 int ID1 = dbobject_num->getTemplateID(detID);
0780 int ID0 = dbobject_den->getTemplateID(detID);
0781
0782 if (ID0 == ID1) {
0783 LogDebug("SiPixelChargeReweightingAlgorithm") << " same template for num and den ";
0784 return false;
0785 }
0786 ierr = PixelTempRewgt2D(ID0, ID1, pixrewgt);
0787 if (ierr != 0) {
0788 LogDebug("SiPixelChargeReweightingAlgorithm ") << "Cluster Charge Reweighting did not work properly.";
0789 return false;
0790 }
0791 if (PrintClusters) {
0792 LogDebug("SiPixelChargeReweightingAlgorithm") << " Cluster after reweighting: ";
0793 printCluster(pixrewgt);
0794 }
0795
0796 for (int row = 0; row < TXSIZE; ++row) {
0797 for (int col = 0; col < TYSIZE; ++col) {
0798 float charge = 0;
0799 charge = pixrewgt[row][col];
0800 if ((hitPixel.first + row - THX) >= 0 && (hitPixel.first + row - THX) < topol->nrows() &&
0801 (hitPixel.second + col - THY) >= 0 && (hitPixel.second + col - THY) < topol->ncolumns() && charge > 0) {
0802 chargeAfter += charge;
0803 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0804 edm::LogError("SiPixelChargeReweightingAlgorithm")
0805 << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
0806 return false;
0807 } else {
0808 theNewDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
0809 AmplitudeType(charge, charge);
0810 }
0811 }
0812 }
0813 }
0814
0815 if (chargeBefore != 0. && chargeAfter == 0.) {
0816 return false;
0817 }
0818
0819 if (PrintClusters) {
0820 LogDebug("SiPixelChargeReweightingAlgorithm")
0821 << "Charges (before->after): " << chargeBefore << " -> " << chargeAfter;
0822 LogDebug("SiPixelChargeReweightingAlgorithm") << "Charge loss: " << (1 - chargeAfter / chargeBefore) * 100 << " %";
0823 }
0824
0825
0826 for (const auto& i : theDigiSignal) {
0827
0828 int chanDigi = i.first;
0829 std::pair<int, int> ip = PixelDigi::channelToPixel(chanDigi);
0830 int row_digi = ip.first;
0831 int col_digi = ip.second;
0832 int i_transformed_row = row_digi - hitPixel.first + THX;
0833 int i_transformed_col = col_digi - hitPixel.second + THY;
0834 if (i_transformed_row < 0 || i_transformed_row > TXSIZE || i_transformed_col < 0 || i_transformed_col > TYSIZE) {
0835
0836 if (chanDigi >= 0 && i.second > 0) {
0837 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0838 edm::LogError("SiPixelChargeReweightingAlgorithm")
0839 << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
0840 return false;
0841 } else {
0842 theNewDigiSignal[chanDigi] += AmplitudeType(i.second, i.second);
0843 }
0844 }
0845 }
0846 }
0847
0848 return true;
0849 }
0850
0851 template <class AmplitudeType, typename SignalType>
0852 bool SiPixelChargeReweightingAlgorithm::lateSignalReweight(const PixelGeomDetUnit* pixdet,
0853 std::vector<PixelDigi>& digis,
0854 PixelSimHitExtraInfoLite& loopTempSH,
0855 SignalType& theNewDigiSignal,
0856 const TrackerTopology* tTopo,
0857 CLHEP::HepRandomEngine* engine) {
0858 uint32_t detID = pixdet->geographicalId().rawId();
0859 const PixelTopology* topol = &pixdet->specificTopology();
0860
0861 if (UseReweighting) {
0862 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0863 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0864 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0865 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ***** INCONSISTENCY !!! ***** \n";
0866 edm::LogError("SiPixelChargeReweightingAlgorithm")
0867 << " applyLateReweighting_ and UseReweighting can not be true at the same time for PU ! \n";
0868 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ---> DO NOT APPLY CHARGE REWEIGHTING TWICE !!! \n";
0869 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0870 edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0871 return false;
0872 }
0873
0874 SignalType theDigiSignal;
0875
0876 int irow_min = topol->nrows();
0877 int irow_max = 0;
0878 int icol_min = topol->ncolumns();
0879 int icol_max = 0;
0880
0881 float chargeBefore = 0;
0882 float chargeAfter = 0;
0883
0884
0885 std::vector<PixelDigi>::const_iterator loopDigi;
0886 for (loopDigi = digis.begin(); loopDigi != digis.end(); ++loopDigi) {
0887 unsigned int chan = loopDigi->channel();
0888 if (loopTempSH.isInTheList(chan)) {
0889 float corresponding_charge = loopDigi->adc();
0890 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0891 edm::LogError("SiPixelChargeReweightingAlgorithm")
0892 << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
0893 return false;
0894 } else {
0895 theDigiSignal[chan] += AmplitudeType(corresponding_charge, corresponding_charge);
0896 }
0897 chargeBefore += corresponding_charge;
0898 if (loopDigi->row() < irow_min)
0899 irow_min = loopDigi->row();
0900 if (loopDigi->row() > irow_max)
0901 irow_max = loopDigi->row();
0902 if (loopDigi->column() < icol_min)
0903 icol_min = loopDigi->column();
0904 if (loopDigi->column() > icol_max)
0905 icol_max = loopDigi->column();
0906 }
0907 }
0908
0909
0910 LocalVector direction = loopTempSH.exitPoint() - loopTempSH.entryPoint();
0911 LocalPoint hitEntryPoint = loopTempSH.entryPoint();
0912 float trajectoryScaleToPosition = std::abs(hitEntryPoint.z() / direction.z());
0913 LocalPoint hitPosition = hitEntryPoint + trajectoryScaleToPosition * direction;
0914
0915
0916
0917 LocalPoint hitLocalPosition = hitEntryPoint + 0.5 * direction;
0918 MeasurementPoint hitPositionPixel = topol->measurementPosition(hitLocalPosition);
0919 std::pair<int, int> hitPixel =
0920 std::pair<int, int>(int(floor(hitPositionPixel.x())), int(floor(hitPositionPixel.y())));
0921
0922 MeasurementPoint originPixel = MeasurementPoint(hitPixel.first - THX + 0.5, hitPixel.second - THY + 0.5);
0923 LocalPoint origin = topol->localPosition(originPixel);
0924
0925 MeasurementPoint hitEntryPointPixel = topol->measurementPosition(loopTempSH.entryPoint());
0926 MeasurementPoint hitExitPointPixel = topol->measurementPosition(loopTempSH.exitPoint());
0927 std::pair<int, int> entryPixel =
0928 std::pair<int, int>(int(floor(hitEntryPointPixel.x())), int(floor(hitEntryPointPixel.y())));
0929 std::pair<int, int> exitPixel =
0930 std::pair<int, int>(int(floor(hitExitPointPixel.x())), int(floor(hitExitPointPixel.y())));
0931
0932 int hitcol_min, hitcol_max, hitrow_min, hitrow_max;
0933 if (entryPixel.first > exitPixel.first) {
0934 hitrow_min = exitPixel.first;
0935 hitrow_max = entryPixel.first;
0936 } else {
0937 hitrow_min = entryPixel.first;
0938 hitrow_max = exitPixel.first;
0939 }
0940
0941 if (entryPixel.second > exitPixel.second) {
0942 hitcol_min = exitPixel.second;
0943 hitcol_max = entryPixel.second;
0944 } else {
0945 hitcol_min = entryPixel.second;
0946 hitcol_max = exitPixel.second;
0947 }
0948
0949 if (!(irow_min <= hitrow_max && irow_max >= hitrow_min && icol_min <= hitcol_max && icol_max >= hitcol_min)) {
0950
0951 return false;
0952 }
0953
0954 track[0] = (hitPosition.x() - origin.x()) * cmToMicrons;
0955 track[1] = (hitPosition.y() - origin.y()) * cmToMicrons;
0956 track[2] = 0.0f;
0957 track[3] = direction.x();
0958 track[4] = direction.y();
0959 track[5] = direction.z();
0960
0961 array_2d pixrewgt(boost::extents[TXSIZE][TYSIZE]);
0962
0963
0964
0965
0966
0967
0968
0969
0970
0971 for (int row = 0; row < TXSIZE; ++row) {
0972 xdouble[row] = topol->isItBigPixelInX(hitPixel.first + row - THX);
0973 }
0974
0975 for (int col = 0; col < TYSIZE; ++col) {
0976 ydouble[col] = topol->isItBigPixelInY(hitPixel.second + col - THY);
0977 }
0978
0979 for (int row = 0; row < TXSIZE; ++row) {
0980 for (int col = 0; col < TYSIZE; ++col) {
0981
0982 pixrewgt[row][col] =
0983 theDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)];
0984 }
0985 }
0986
0987 if (PrintClusters) {
0988 LogDebug("SiPixelChargeReweightingAlgorithm") << " Cluster before reweighting: ";
0989 printCluster(pixrewgt);
0990 }
0991
0992 int ierr;
0993
0994
0995 int ID1 = dbobject_num->getTemplateID(detID);
0996 int ID0 = dbobject_den->getTemplateID(detID);
0997
0998 if (ID0 == ID1) {
0999 LogDebug("SiPixelChargeReweightingAlgorithm") << " same template for num and den ";
1000 return false;
1001 }
1002 ierr = PixelTempRewgt2D(ID0, ID1, pixrewgt);
1003 if (ierr != 0) {
1004 LogDebug("SiPixelChargeReweightingAlgorithm") << "Cluster Charge Reweighting did not work properly.";
1005 return false;
1006 }
1007 if (PrintClusters) {
1008 LogDebug("SiPixelChargeReweightingAlgorithm") << " Cluster after reweighting: ";
1009 printCluster(pixrewgt);
1010 }
1011
1012 for (int row = 0; row < TXSIZE; ++row) {
1013 for (int col = 0; col < TYSIZE; ++col) {
1014 float charge = 0;
1015 charge = pixrewgt[row][col];
1016 if ((hitPixel.first + row - THX) >= 0 && (hitPixel.first + row - THX) < topol->nrows() &&
1017 (hitPixel.second + col - THY) >= 0 && (hitPixel.second + col - THY) < topol->ncolumns() && charge > 0) {
1018 chargeAfter += charge;
1019 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
1020 edm::LogError("SiPixelChargeReweightingAlgorithm")
1021 << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
1022 return false;
1023 } else {
1024 theNewDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
1025 AmplitudeType(charge, charge);
1026 }
1027 }
1028 }
1029 }
1030
1031 if (chargeBefore != 0. && chargeAfter == 0.) {
1032 return false;
1033 }
1034
1035 if (PrintClusters) {
1036 LogDebug("SiPixelChargeReweightingAlgorithm")
1037 << "Charges (before->after): " << chargeBefore << " -> " << chargeAfter;
1038 LogDebug("SiPixelChargeReweightingAlgorithm") << "Charge loss: " << (1 - chargeAfter / chargeBefore) * 100 << " %";
1039 }
1040
1041
1042 for (const auto& i : theDigiSignal) {
1043
1044 int chanDigi = i.first;
1045 std::pair<int, int> ip = PixelDigi::channelToPixel(chanDigi);
1046 int row_digi = ip.first;
1047 int col_digi = ip.second;
1048 int i_transformed_row = row_digi - hitPixel.first + THX;
1049 int i_transformed_col = col_digi - hitPixel.second + THY;
1050 if (i_transformed_row < 0 || i_transformed_row > TXSIZE || i_transformed_col < 0 || i_transformed_col > TYSIZE) {
1051
1052 if (chanDigi >= 0 && i.second > 0) {
1053 if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
1054 edm::LogError("SiPixelChargeReweightingAlgorithm")
1055 << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
1056 return false;
1057 } else {
1058 theNewDigiSignal[chanDigi] += AmplitudeType(i.second, i.second);
1059 }
1060 }
1061 }
1062 }
1063
1064 return true;
1065 }
1066 #endif