Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-04 01:26:22

0001 //class SiPixelChargeReweightingAlgorithm SimTracker/SiPixelDigitizer/src/SiPixelChargeReweightingAlgorithm.cc
0002 
0003 // Original Author Caroline Collard
0004 // September 2020 : Extraction of the code for cluster charge reweighting from SiPixelDigitizerAlgorithm to a new class
0005 //
0006 #include <iostream>
0007 #include <iomanip>
0008 
0009 #include "SimGeneral/NoiseGenerators/interface/GaussianTailNoiseGenerator.h"
0010 
0011 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0012 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0013 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0014 #include "SimTracker/Common/interface/SiG4UniversalFluctuation.h"
0015 #include "SimTracker/SiPixelDigitizer/plugins/SiPixelChargeReweightingAlgorithm.h"
0016 
0017 //#include "PixelIndices.h"
0018 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0019 
0020 #include "FWCore/Framework/interface/ConsumesCollector.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 
0024 // Accessing dead pixel modules from the DB:
0025 #include "DataFormats/DetId/interface/DetId.h"
0026 
0027 #include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
0028 
0029 #include "CondFormats/SiPixelObjects/interface/PixelIndices.h"
0030 #include "CondFormats/SiPixelObjects/interface/PixelROC.h"
0031 #include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
0032 #include "CondFormats/SiPixelObjects/interface/SiPixel2DTemplateDBObject.h"
0033 
0034 #include "CondFormats/SiPixelObjects/interface/PixelROC.h"
0035 
0036 // Geometry
0037 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0038 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0039 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0040 
0041 using namespace edm;
0042 using namespace sipixelobjects;
0043 
0044 void SiPixelChargeReweightingAlgorithm::init(const edm::EventSetup& es) {
0045   // Read template files for charge reweighting
0046   if (UseReweighting || applyLateReweighting_) {
0047     dbobject_den = &es.getData(SiPixel2DTemp_den_token_);
0048     dbobject_num = &es.getData(SiPixel2DTemp_num_token_);
0049 
0050     int numOfTemplates = dbobject_den->numOfTempl() + dbobject_num->numOfTempl();
0051     templateStores_.reserve(numOfTemplates);
0052     SiPixelTemplate2D::pushfile(*dbobject_den, templateStores_);
0053     SiPixelTemplate2D::pushfile(*dbobject_num, templateStores_);
0054 
0055     track.resize(6);
0056   }
0057 }
0058 
0059 //=========================================================================
0060 
0061 SiPixelChargeReweightingAlgorithm::SiPixelChargeReweightingAlgorithm(const edm::ParameterSet& conf,
0062                                                                      edm::ConsumesCollector iC)
0063     :
0064 
0065       templ2D(templateStores_),
0066       xdouble(TXSIZE),
0067       ydouble(TYSIZE),
0068       IDnum(conf.exists("TemplateIDnumerator") ? conf.getParameter<int>("TemplateIDnumerator") : 0),
0069       IDden(conf.exists("TemplateIDdenominator") ? conf.getParameter<int>("TemplateIDdenominator") : 0),
0070 
0071       UseReweighting(conf.getParameter<bool>("UseReweighting")),
0072       applyLateReweighting_(conf.exists("applyLateReweighting") ? conf.getParameter<bool>("applyLateReweighting")
0073                                                                 : false),
0074       PrintClusters(conf.getParameter<bool>("PrintClusters")),
0075       PrintTemplates(conf.getParameter<bool>("PrintTemplates")) {
0076   if (UseReweighting || applyLateReweighting_) {
0077     SiPixel2DTemp_den_token_ = iC.esConsumes(edm::ESInputTag("", "denominator"));
0078     SiPixel2DTemp_num_token_ = iC.esConsumes(edm::ESInputTag("", "numerator"));
0079   }
0080   edm::LogVerbatim("PixelDigitizer ") << "SiPixelChargeReweightingAlgorithm constructed"
0081                                       << " with UseReweighting = " << UseReweighting
0082                                       << " and applyLateReweighting_ = " << applyLateReweighting_;
0083 }
0084 
0085 //=========================================================================
0086 SiPixelChargeReweightingAlgorithm::~SiPixelChargeReweightingAlgorithm() {
0087   LogDebug("PixelDigitizer") << "SiPixelChargeReweightingAlgorithm deleted";
0088 }
0089 
0090 //============================================================================
0091 
0092 bool SiPixelChargeReweightingAlgorithm::hitSignalReweight(const PSimHit& hit,
0093                                                           std::map<int, float, std::less<int> >& hit_signal,
0094                                                           const size_t hitIndex,
0095                                                           const size_t hitIndex4CR,
0096                                                           const unsigned int tofBin,
0097                                                           const PixelTopology* topol,
0098                                                           uint32_t detID,
0099                                                           signal_map_type& theSignal,
0100                                                           unsigned short int processType,
0101                                                           const bool& boolmakeDigiSimLinks) {
0102   int irow_min = topol->nrows();
0103   int irow_max = 0;
0104   int icol_min = topol->ncolumns();
0105   int icol_max = 0;
0106 
0107   float chargeBefore = 0;
0108   float chargeAfter = 0;
0109   signal_map_type hitSignal;
0110   LocalVector direction = hit.exitPoint() - hit.entryPoint();
0111 
0112   for (std::map<int, float, std::less<int> >::const_iterator im = hit_signal.begin(); im != hit_signal.end(); ++im) {
0113     int chan = (*im).first;
0114     std::pair<int, int> pixelWithCharge = PixelDigi::channelToPixel(chan);
0115 
0116     hitSignal[chan] += (boolmakeDigiSimLinks ? SiPixelDigitizerAlgorithm::Amplitude(
0117                                                    (*im).second, &hit, hitIndex, hitIndex4CR, tofBin, (*im).second)
0118                                              : SiPixelDigitizerAlgorithm::Amplitude((*im).second, (*im).second));
0119     chargeBefore += (*im).second;
0120 
0121     if (pixelWithCharge.first < irow_min)
0122       irow_min = pixelWithCharge.first;
0123     if (pixelWithCharge.first > irow_max)
0124       irow_max = pixelWithCharge.first;
0125     if (pixelWithCharge.second < icol_min)
0126       icol_min = pixelWithCharge.second;
0127     if (pixelWithCharge.second > icol_max)
0128       icol_max = pixelWithCharge.second;
0129   }
0130 
0131   LocalPoint hitEntryPoint = hit.entryPoint();
0132 
0133   float trajectoryScaleToPosition = std::abs(hitEntryPoint.z() / direction.z());
0134 
0135   LocalPoint hitPosition = hitEntryPoint + trajectoryScaleToPosition * direction;
0136 
0137   MeasurementPoint hitPositionPixel = topol->measurementPosition(hit.localPosition());
0138   std::pair<int, int> hitPixel =
0139       std::pair<int, int>(int(floor(hitPositionPixel.x())), int(floor(hitPositionPixel.y())));
0140 
0141   MeasurementPoint originPixel = MeasurementPoint(hitPixel.first - THX + 0.5, hitPixel.second - THY + 0.5);
0142   LocalPoint origin = topol->localPosition(originPixel);
0143 
0144   MeasurementPoint hitEntryPointPixel = topol->measurementPosition(hit.entryPoint());
0145   MeasurementPoint hitExitPointPixel = topol->measurementPosition(hit.exitPoint());
0146   std::pair<int, int> entryPixel =
0147       std::pair<int, int>(int(floor(hitEntryPointPixel.x())), int(floor(hitEntryPointPixel.y())));
0148   std::pair<int, int> exitPixel =
0149       std::pair<int, int>(int(floor(hitExitPointPixel.x())), int(floor(hitExitPointPixel.y())));
0150 
0151   int hitcol_min, hitcol_max, hitrow_min, hitrow_max;
0152   if (entryPixel.first > exitPixel.first) {
0153     hitrow_min = exitPixel.first;
0154     hitrow_max = entryPixel.first;
0155   } else {
0156     hitrow_min = entryPixel.first;
0157     hitrow_max = exitPixel.first;
0158   }
0159 
0160   if (entryPixel.second > exitPixel.second) {
0161     hitcol_min = exitPixel.second;
0162     hitcol_max = entryPixel.second;
0163   } else {
0164     hitcol_min = entryPixel.second;
0165     hitcol_max = exitPixel.second;
0166   }
0167 
0168 #ifdef TP_DEBUG
0169   LocalPoint CMSSWhitPosition = hit.localPosition();
0170 
0171   LogDebug("Pixel Digitizer") << "\n"
0172                               << "Particle ID is: " << hit.particleType() << "\n"
0173                               << "Process type: " << hit.processType() << "\n"
0174                               << "HitPosition:"
0175                               << "\n"
0176                               << "Hit entry x/y/z: " << hit.entryPoint().x() << "  " << hit.entryPoint().y() << "  "
0177                               << hit.entryPoint().z() << "  "
0178                               << "Hit exit x/y/z: " << hit.exitPoint().x() << "  " << hit.exitPoint().y() << "  "
0179                               << hit.exitPoint().z() << "  "
0180 
0181                               << "Pixel Pos - X: " << hitPositionPixel.x() << " Y: " << hitPositionPixel.y() << "\n"
0182                               << "Cart.Cor. - X: " << CMSSWhitPosition.x() << " Y: " << CMSSWhitPosition.y() << "\n"
0183                               << "Z=0 Pos - X: " << hitPosition.x() << " Y: " << hitPosition.y() << "\n"
0184 
0185                               << "Origin of the template:"
0186                               << "\n"
0187                               << "Pixel Pos - X: " << originPixel.x() << " Y: " << originPixel.y() << "\n"
0188                               << "Cart.Cor. - X: " << origin.x() << " Y: " << origin.y() << "\n"
0189                               << "\n"
0190                               << "Entry/Exit:"
0191                               << "\n"
0192                               << "Entry - X: " << hit.entryPoint().x() << " Y: " << hit.entryPoint().y()
0193                               << " Z: " << hit.entryPoint().z() << "\n"
0194                               << "Exit - X: " << hit.exitPoint().x() << " Y: " << hit.exitPoint().y()
0195                               << " Z: " << hit.exitPoint().z() << "\n"
0196 
0197                               << "Entry - X Pixel: " << hitEntryPointPixel.x() << " Y Pixel: " << hitEntryPointPixel.y()
0198                               << "\n"
0199                               << "Exit - X Pixel: " << hitExitPointPixel.x() << " Y Pixel: " << hitExitPointPixel.y()
0200                               << "\n"
0201 
0202                               << "row min: " << irow_min << " col min: " << icol_min << "\n";
0203 #endif
0204 
0205   if (!(irow_min <= hitrow_max && irow_max >= hitrow_min && icol_min <= hitcol_max && icol_max >= hitcol_min)) {
0206     // The clusters do not have an overlap, hence the hit is NOT reweighted
0207     return false;
0208   }
0209 
0210   track[0] = (hitPosition.x() - origin.x()) * cmToMicrons;
0211   track[1] = (hitPosition.y() - origin.y()) * cmToMicrons;
0212   track[2] = 0.0f;  //Middle of sensor is origin for Z-axis
0213   track[3] = direction.x();
0214   track[4] = direction.y();
0215   track[5] = direction.z();
0216 
0217   array_2d pixrewgt(boost::extents[TXSIZE][TYSIZE]);
0218 
0219   for (int row = 0; row < TXSIZE; ++row) {
0220     for (int col = 0; col < TYSIZE; ++col) {
0221       pixrewgt[row][col] = 0;
0222     }
0223   }
0224 
0225   for (int row = 0; row < TXSIZE; ++row) {
0226     xdouble[row] = topol->isItBigPixelInX(hitPixel.first + row - THX);
0227   }
0228 
0229   for (int col = 0; col < TYSIZE; ++col) {
0230     ydouble[col] = topol->isItBigPixelInY(hitPixel.second + col - THY);
0231   }
0232 
0233   // define loop boundaries that will prevent the row and col loops
0234   // from going out of physical bounds of the pixel module
0235   int rowmin = std::max(0, THX - hitPixel.first);
0236   int rowmax = std::min(TXSIZE, topol->nrows() + THX - hitPixel.first);
0237   int colmin = std::max(0, THY - hitPixel.second);
0238   int colmax = std::min(TYSIZE, topol->ncolumns() + THY - hitPixel.second);
0239 
0240   for (int row = rowmin; row < rowmax; ++row) {
0241     for (int col = colmin; col < colmax; ++col) {
0242       //Fill charges into 21x13 Pixel Array with hitPixel in centre
0243       pixrewgt[row][col] =
0244           hitSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)];
0245     }
0246   }
0247 
0248   if (PrintClusters) {
0249     LogDebug("PixelDigitizer ") << "Cluster before reweighting: ";
0250     printCluster(pixrewgt);
0251   }
0252 
0253   int ierr;
0254   // for unirradiated: 2nd argument is IDden
0255   // for irradiated: 2nd argument is IDnum
0256   if (UseReweighting == true) {
0257     int ID1 = dbobject_num->getTemplateID(detID);
0258     int ID0 = dbobject_den->getTemplateID(detID);
0259 
0260     if (ID0 == ID1) {
0261       return false;
0262     }
0263     ierr = PixelTempRewgt2D(ID0, ID1, pixrewgt);
0264   } else {
0265     ierr = PixelTempRewgt2D(IDden, IDden, pixrewgt);
0266   }
0267   if (ierr != 0) {
0268 #ifdef TP_DEBUG
0269     LogDebug("PixelDigitizer ") << "Cluster Charge Reweighting did not work properly.";
0270 #endif
0271     return false;
0272   }
0273 
0274   if (PrintClusters) {
0275     LogDebug("PixelDigitizer ") << "Cluster after reweighting: ";
0276     printCluster(pixrewgt);
0277   }
0278 
0279   for (int row = rowmin; row < rowmax; ++row) {
0280     for (int col = colmin; col < colmax; ++col) {
0281       float charge = pixrewgt[row][col];
0282       if (charge > 0) {
0283         chargeAfter += charge;
0284         theSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
0285             (boolmakeDigiSimLinks
0286                  ? SiPixelDigitizerAlgorithm::Amplitude(charge, &hit, hitIndex, hitIndex4CR, tofBin, charge)
0287                  : SiPixelDigitizerAlgorithm::Amplitude(charge, charge));
0288       }
0289     }
0290   }
0291 
0292   if (chargeBefore != 0. && chargeAfter == 0.) {
0293     return false;
0294   }
0295 
0296   if (PrintClusters) {
0297     LogDebug("PixelDigitizer ") << "Charges (before->after): " << chargeBefore << " -> " << chargeAfter;
0298     LogDebug("PixelDigitizer ") << "Charge loss: " << (1 - chargeAfter / chargeBefore) * 100 << " % \n";
0299   }
0300 
0301   return true;
0302 }
0303 
0304 // *******************************************************************************************************
0305 //! Reweight CMSSW clusters to look like clusters corresponding to Pixelav Templates.
0306 //! \param       id_in - (input) identifier of the template corresponding to the input events
0307 //! \param    id_rewgt - (input) identifier of the template corresponding to the output events
0308 //! \param     cluster - (input/output) boost multi_array container of 7x21 array of pixel signals,
0309 //!                       origin of local coords (0,0) at center of pixel cluster[3][10].
0310 //! returns 0 if everything is OK, 1 if angles are outside template coverage (cluster is probably still
0311 //! usable, > 1 if something is wrong (no reweight done).
0312 // *******************************************************************************************************
0313 int SiPixelChargeReweightingAlgorithm::PixelTempRewgt2D(int id_in, int id_rewgt, array_2d& cluster) {
0314   // Local variables
0315   int i, j, k, l, kclose;
0316   int nclusx, nclusy, success;
0317   float xsize, ysize, q50i, q100i, q50r, q10r, q100r, xhit2D, yhit2D, qclust, dist2, dmin2;
0318   float xy_in[BXM2][BYM2], xy_rewgt[BXM2][BYM2], xy_clust[TXSIZE][TYSIZE];
0319   int denx_clust[TXSIZE][TYSIZE], deny_clust[TXSIZE][TYSIZE];
0320   int goodWeightsUsed, nearbyWeightsUsed, noWeightsUsed;
0321   float cotalpha, cotbeta;
0322   // success = 0 is returned if everthing is OK
0323   success = 0;
0324 
0325   // Copy the array to remember original charges
0326   array_2d clust(cluster);
0327 
0328   // Take the pixel dimensions from the 2D template
0329   templ2D.getid(id_in);
0330   xsize = templ2D.xsize();
0331   ysize = templ2D.ysize();
0332 
0333   // Calculate the track angles
0334 
0335   if (std::abs(track[5]) > 0.f) {
0336     cotalpha = track[3] / track[5];  //if track[5] (direction in z) is 0 the hit is not processed by re-weighting
0337     cotbeta = track[4] / track[5];
0338   } else {
0339     LogDebug("Pixel Digitizer") << "Reweighting angle is not good! \n";
0340     return 9;  //returned value here indicates that no reweighting was done in this case
0341   }
0342 
0343   // The 2-D templates are defined on a shifted coordinate system wrt the 1D templates
0344   if (ydouble[0]) {
0345     yhit2D = track[1] - cotbeta * track[2] + ysize;
0346   } else {
0347     yhit2D = track[1] - cotbeta * track[2] + 0.5f * ysize;
0348   }
0349   if (xdouble[0]) {
0350     xhit2D = track[0] - cotalpha * track[2] + xsize;
0351   } else {
0352     xhit2D = track[0] - cotalpha * track[2] + 0.5f * xsize;
0353   }
0354 
0355   // Zero the input and output templates
0356   for (i = 0; i < BYM2; ++i) {
0357     for (j = 0; j < BXM2; ++j) {
0358       xy_in[j][i] = 0.f;
0359       xy_rewgt[j][i] = 0.f;
0360     }
0361   }
0362 
0363   // Next, interpolate the CMSSW template needed to analyze this cluster
0364 
0365   if (!templ2D.xytemp(id_in, cotalpha, cotbeta, xhit2D, yhit2D, ydouble, xdouble, xy_in)) {
0366     success = 1;
0367   }
0368   if (success != 0) {
0369 #ifdef TP_DEBUG
0370     LogDebug("Pixel Digitizer") << "No matching template found \n";
0371 #endif
0372     return 2;
0373   }
0374 
0375   if (PrintTemplates) {
0376     LogDebug("Pixel Digitizer") << "Template unirrad: \n";
0377     printCluster(xy_in);
0378   }
0379 
0380   q50i = templ2D.s50();
0381   //q50i = 0;
0382   q100i = 2.f * q50i;
0383 
0384   // Check that the cluster container is a 13x21 matrix
0385 
0386   if (cluster.num_dimensions() != 2) {
0387     LogWarning("Pixel Digitizer") << "Cluster is not 2-dimensional. Return. \n";
0388     return 3;
0389   }
0390   nclusx = (int)cluster.shape()[0];
0391   nclusy = (int)cluster.shape()[1];
0392   if (nclusx != TXSIZE || xdouble.size() != TXSIZE) {
0393     LogWarning("Pixel Digitizer") << "Sizes in x do not match: nclusx=" << nclusx << "  xdoubleSize=" << xdouble.size()
0394                                   << "  TXSIZE=" << TXSIZE << ". Return. \n";
0395     return 4;
0396   }
0397   if (nclusy != TYSIZE || ydouble.size() != TYSIZE) {
0398     LogWarning("Pixel Digitizer") << "Sizes in y do not match. Return. \n";
0399     return 5;
0400   }
0401 
0402   // Sum initial charge in the cluster
0403 
0404   qclust = 0.f;
0405   for (i = 0; i < TYSIZE; ++i) {
0406     for (j = 0; j < TXSIZE; ++j) {
0407       xy_clust[j][i] = 0.f;
0408       denx_clust[j][i] = 0;
0409       deny_clust[j][i] = 0;
0410       if (cluster[j][i] > q100i) {
0411         qclust += cluster[j][i];
0412       }
0413     }
0414   }
0415 
0416   // Next, interpolate the physical output template needed to reweight
0417 
0418   if (!templ2D.xytemp(id_rewgt, cotalpha, cotbeta, xhit2D, yhit2D, ydouble, xdouble, xy_rewgt)) {
0419     success = 1;
0420   }
0421 
0422   if (PrintTemplates) {
0423     LogDebug("Pixel Digitizer") << "Template irrad: \n";
0424     printCluster(xy_rewgt);
0425   }
0426 
0427   q50r = templ2D.s50();
0428   q100r = 2.f * q50r;
0429   q10r = 0.2f * q50r;
0430 
0431   // Find all non-zero denominator pixels in the input template and generate "inside" weights
0432 
0433   int ntpix = 0;
0434   int ncpix = 0;
0435   std::vector<int> ytclust;
0436   std::vector<int> xtclust;
0437   std::vector<int> ycclust;
0438   std::vector<int> xcclust;
0439   qclust = 0.f;
0440   for (i = 0; i < TYSIZE; ++i) {
0441     for (j = 0; j < TXSIZE; ++j) {
0442       if (xy_in[j + 1][i + 1] > q100i) {
0443         ++ntpix;
0444         ytclust.push_back(i);
0445         xtclust.push_back(j);
0446         xy_clust[j][i] = xy_rewgt[j + 1][i + 1] / xy_in[j + 1][i + 1];
0447         denx_clust[j][i] = j;
0448         deny_clust[j][i] = i;
0449       }
0450     }
0451   }
0452 
0453   // Find all non-zero numerator pixels not matched to denominator in the output template and generate "inside" weights
0454 
0455   for (i = 0; i < TYSIZE; ++i) {
0456     for (j = 0; j < TXSIZE; ++j) {
0457       if (xy_rewgt[j + 1][i + 1] > q10r && xy_clust[j][i] == 0.f && ntpix > 0) {
0458         // Search for nearest denominator pixel
0459         dmin2 = 10000.f;
0460         kclose = 0;
0461         for (k = 0; k < ntpix; ++k) {
0462           dist2 = (i - ytclust[k]) * (i - ytclust[k]) + 0.44444f * (j - xtclust[k]) * (j - xtclust[k]);
0463           if (dist2 < dmin2) {
0464             dmin2 = dist2;
0465             kclose = k;
0466           }
0467         }
0468         xy_clust[j][i] = xy_rewgt[j + 1][i + 1] / xy_in[xtclust[kclose] + 1][ytclust[kclose] + 1];
0469         denx_clust[j][i] = xtclust[kclose];
0470         deny_clust[j][i] = ytclust[kclose];
0471       }
0472     }
0473   }
0474 
0475   if (PrintTemplates) {
0476     LogDebug("Pixel Digitizer") << "Weights: \n";
0477     printCluster(xy_clust);
0478   }
0479 
0480   // Do the reweighting
0481   goodWeightsUsed = 0;
0482   nearbyWeightsUsed = 0;
0483   noWeightsUsed = 0;
0484 
0485   for (i = 0; i < TYSIZE; ++i) {
0486     for (j = 0; j < TXSIZE; ++j) {
0487       if (xy_clust[j][i] > 0.f) {
0488         cluster[j][i] = xy_clust[j][i] * clust[denx_clust[j][i]][deny_clust[j][i]];
0489         if (cluster[j][i] > q100r) {
0490           qclust += cluster[j][i];
0491         }
0492         if (cluster[j][i] > 0) {
0493           goodWeightsUsed++;
0494         }
0495       } else {
0496         if (clust[j][i] > 0.f) {
0497           ++ncpix;
0498           ycclust.push_back(i);
0499           xcclust.push_back(j);
0500         }
0501       }
0502     }
0503   }
0504 
0505   // Now reweight pixels outside of template footprint using closest weights
0506 
0507   if (ncpix > 0) {
0508     for (l = 0; l < ncpix; ++l) {
0509       i = ycclust[l];
0510       j = xcclust[l];
0511       dmin2 = 10000.f;
0512       kclose = 0;
0513       for (k = 0; k < ntpix; ++k) {
0514         dist2 = (i - ytclust[k]) * (i - ytclust[k]) + 0.44444f * (j - xtclust[k]) * (j - xtclust[k]);
0515         if (dist2 < dmin2) {
0516           dmin2 = dist2;
0517           kclose = k;
0518         }
0519       }
0520       if (dmin2 < 5.f) {
0521         nearbyWeightsUsed++;
0522         cluster[j][i] *= xy_clust[xtclust[kclose]][ytclust[kclose]];
0523         if (cluster[j][i] > q100r) {
0524           qclust += cluster[j][i];
0525         }
0526       } else {
0527         noWeightsUsed++;
0528         cluster[j][i] = 0.f;
0529       }
0530     }
0531   }
0532 
0533   return success;
0534 }  // PixelTempRewgt2D
0535 
0536 void SiPixelChargeReweightingAlgorithm::printCluster(array_2d& cluster) {
0537   for (int col = 0; col < TYSIZE; ++col) {
0538     for (int row = 0; row < TXSIZE; ++row) {
0539       LogDebug("Pixel Digitizer") << cluster[row][col];
0540     }
0541     LogDebug("Pixel Digitizer") << "\n";
0542   }
0543 }
0544 
0545 void SiPixelChargeReweightingAlgorithm::printCluster(float arr[BXM2][BYM2]) {
0546   for (int col = 0; col < BYM2; ++col) {
0547     for (int row = 0; row < BXM2; ++row) {
0548       LogDebug("Pixel Digitizer") << arr[row][col];
0549     }
0550     LogDebug("Pixel Digitizer") << "\n";
0551   }
0552 }
0553 
0554 void SiPixelChargeReweightingAlgorithm::printCluster(float arr[TXSIZE][TYSIZE]) {
0555   for (int col = 0; col < TYSIZE; ++col) {
0556     for (int row = 0; row < TXSIZE; ++row) {
0557       LogDebug("Pixel Digitizer") << arr[row][col];
0558     }
0559     LogDebug("Pixel Digitizer") << "\n";
0560   }
0561 }
0562 
0563 bool SiPixelChargeReweightingAlgorithm::lateSignalReweight(const PixelGeomDetUnit* pixdet,
0564                                                            std::vector<PixelDigi>& digis,
0565                                                            PixelSimHitExtraInfo& loopTempSH,
0566                                                            signal_map_type& theNewDigiSignal,
0567                                                            const TrackerTopology* tTopo,
0568                                                            CLHEP::HepRandomEngine* engine) {
0569   uint32_t detID = pixdet->geographicalId().rawId();
0570   const PixelTopology* topol = &pixdet->specificTopology();
0571 
0572   if (UseReweighting) {
0573     LogError("Pixel Digitizer") << " ******************************** \n";
0574     LogError("Pixel Digitizer") << " ******************************** \n";
0575     LogError("Pixel Digitizer") << " ******************************** \n";
0576     LogError("Pixel Digitizer") << " *****  INCONSISTENCY !!!   ***** \n";
0577     LogError("Pixel Digitizer")
0578         << " applyLateReweighting_ and UseReweighting can not be true at the same time for PU ! \n";
0579     LogError("Pixel Digitizer") << " ---> DO NOT APPLY CHARGE REWEIGHTING TWICE !!! \n";
0580     LogError("Pixel Digitizer") << " ******************************** \n";
0581     LogError("Pixel Digitizer") << " ******************************** \n";
0582     return false;
0583   }
0584 
0585   signal_map_type theDigiSignal;
0586 
0587   int irow_min = topol->nrows();
0588   int irow_max = 0;
0589   int icol_min = topol->ncolumns();
0590   int icol_max = 0;
0591 
0592   float chargeBefore = 0;
0593   float chargeAfter = 0;
0594 
0595   //loop on digis
0596   std::vector<PixelDigi>::const_iterator loopDigi;
0597   for (loopDigi = digis.begin(); loopDigi != digis.end(); ++loopDigi) {
0598     unsigned int chan = loopDigi->channel();
0599     if (loopTempSH.isInTheList(chan)) {
0600       float corresponding_charge = loopDigi->adc();
0601       theDigiSignal[chan] += SiPixelDigitizerAlgorithm::Amplitude(corresponding_charge, corresponding_charge);
0602       chargeBefore += corresponding_charge;
0603       if (loopDigi->row() < irow_min)
0604         irow_min = loopDigi->row();
0605       if (loopDigi->row() > irow_max)
0606         irow_max = loopDigi->row();
0607       if (loopDigi->column() < icol_min)
0608         icol_min = loopDigi->column();
0609       if (loopDigi->column() > icol_max)
0610         icol_max = loopDigi->column();
0611     }
0612   }
0613   // end loop on digis
0614 
0615   LocalVector direction = loopTempSH.exitPoint() - loopTempSH.entryPoint();
0616   LocalPoint hitEntryPoint = loopTempSH.entryPoint();
0617   float trajectoryScaleToPosition = std::abs(hitEntryPoint.z() / direction.z());
0618   LocalPoint hitPosition = hitEntryPoint + trajectoryScaleToPosition * direction;
0619 
0620   // addition as now the hit himself is not available
0621   // https://github.com/cms-sw/cmssw/blob/master/SimDataFormats/TrackingHit/interface/PSimHit.h#L52
0622   LocalPoint hitLocalPosition = hitEntryPoint + 0.5 * direction;
0623   MeasurementPoint hitPositionPixel = topol->measurementPosition(hitLocalPosition);
0624   std::pair<int, int> hitPixel =
0625       std::pair<int, int>(int(floor(hitPositionPixel.x())), int(floor(hitPositionPixel.y())));
0626 
0627   MeasurementPoint originPixel = MeasurementPoint(hitPixel.first - THX + 0.5, hitPixel.second - THY + 0.5);
0628   LocalPoint origin = topol->localPosition(originPixel);
0629 
0630   MeasurementPoint hitEntryPointPixel = topol->measurementPosition(loopTempSH.entryPoint());
0631   MeasurementPoint hitExitPointPixel = topol->measurementPosition(loopTempSH.exitPoint());
0632   std::pair<int, int> entryPixel =
0633       std::pair<int, int>(int(floor(hitEntryPointPixel.x())), int(floor(hitEntryPointPixel.y())));
0634   std::pair<int, int> exitPixel =
0635       std::pair<int, int>(int(floor(hitExitPointPixel.x())), int(floor(hitExitPointPixel.y())));
0636 
0637   int hitcol_min, hitcol_max, hitrow_min, hitrow_max;
0638   if (entryPixel.first > exitPixel.first) {
0639     hitrow_min = exitPixel.first;
0640     hitrow_max = entryPixel.first;
0641   } else {
0642     hitrow_min = entryPixel.first;
0643     hitrow_max = exitPixel.first;
0644   }
0645 
0646   if (entryPixel.second > exitPixel.second) {
0647     hitcol_min = exitPixel.second;
0648     hitcol_max = entryPixel.second;
0649   } else {
0650     hitcol_min = entryPixel.second;
0651     hitcol_max = exitPixel.second;
0652   }
0653 
0654   if (!(irow_min <= hitrow_max && irow_max >= hitrow_min && icol_min <= hitcol_max && icol_max >= hitcol_min)) {
0655     // The clusters do not have an overlap, hence the hit is NOT reweighted
0656     return false;
0657   }
0658 
0659   track[0] = (hitPosition.x() - origin.x()) * cmToMicrons;
0660   track[1] = (hitPosition.y() - origin.y()) * cmToMicrons;
0661   track[2] = 0.0f;  //Middle of sensor is origin for Z-axis
0662   track[3] = direction.x();
0663   track[4] = direction.y();
0664   track[5] = direction.z();
0665 
0666   array_2d pixrewgt(boost::extents[TXSIZE][TYSIZE]);
0667 
0668   /*
0669   for (int row = 0; row < TXSIZE; ++row) {
0670     for (int col = 0; col < TYSIZE; ++col) {
0671       pixrewgt[row][col] = 0;
0672     }
0673   }
0674 */
0675 
0676   for (int row = 0; row < TXSIZE; ++row) {
0677     xdouble[row] = topol->isItBigPixelInX(hitPixel.first + row - THX);
0678   }
0679 
0680   for (int col = 0; col < TYSIZE; ++col) {
0681     ydouble[col] = topol->isItBigPixelInY(hitPixel.second + col - THY);
0682   }
0683 
0684   for (int row = 0; row < TXSIZE; ++row) {
0685     for (int col = 0; col < TYSIZE; ++col) {
0686       //Fill charges into 21x13 Pixel Array with hitPixel in centre
0687       pixrewgt[row][col] =
0688           theDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)];
0689     }
0690   }
0691 
0692   if (PrintClusters) {
0693     LogDebug("Pixel Digitizer") << " Cluster before reweighting: ";
0694     printCluster(pixrewgt);
0695   }
0696 
0697   int ierr;
0698   // for unirradiated: 2nd argument is IDden
0699   // for irradiated: 2nd argument is IDnum
0700   int ID1 = dbobject_num->getTemplateID(detID);
0701   int ID0 = dbobject_den->getTemplateID(detID);
0702 
0703   if (ID0 == ID1) {
0704     LogDebug("Pixel Digitizer") << " same template for num and den ";
0705     return false;
0706   }
0707   ierr = PixelTempRewgt2D(ID0, ID1, pixrewgt);
0708   if (ierr != 0) {
0709 #ifdef TP_DEBUG
0710     LogDebug("PixelDigitizer ") << "Cluster Charge Reweighting did not work properly.";
0711 #endif
0712     return false;
0713   }
0714   if (PrintClusters) {
0715     LogDebug("Pixel Digitizer") << " Cluster after reweighting: ";
0716     printCluster(pixrewgt);
0717   }
0718 
0719   for (int row = 0; row < TXSIZE; ++row) {
0720     for (int col = 0; col < TYSIZE; ++col) {
0721       float charge = 0;
0722       charge = pixrewgt[row][col];
0723       if ((hitPixel.first + row - THX) >= 0 && (hitPixel.first + row - THX) < topol->nrows() &&
0724           (hitPixel.second + col - THY) >= 0 && (hitPixel.second + col - THY) < topol->ncolumns() && charge > 0) {
0725         chargeAfter += charge;
0726         theNewDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
0727             SiPixelDigitizerAlgorithm::Amplitude(charge, charge);
0728       }
0729     }
0730   }
0731 
0732   if (chargeBefore != 0. && chargeAfter == 0.) {
0733     return false;
0734   }
0735 
0736   if (PrintClusters) {
0737     LogDebug("Pixel Digitizer") << "Charges (before->after): " << chargeBefore << " -> " << chargeAfter;
0738     LogDebug("Pixel Digitizer") << "Charge loss: " << (1 - chargeAfter / chargeBefore) * 100 << " %";
0739   }
0740 
0741   // need to store the digi out of the 21x13 box.
0742   for (signal_map_const_iterator i = theDigiSignal.begin(); i != theDigiSignal.end(); ++i) {
0743     // check if in the 21x13 box
0744     int chanDigi = (*i).first;
0745     std::pair<int, int> ip = PixelDigi::channelToPixel(chanDigi);
0746     int row_digi = ip.first;
0747     int col_digi = ip.second;
0748     int i_transformed_row = row_digi - hitPixel.first + THX;
0749     int i_transformed_col = col_digi - hitPixel.second + THY;
0750     if (i_transformed_row < 0 || i_transformed_row > TXSIZE || i_transformed_col < 0 || i_transformed_col > TYSIZE) {
0751       // not in the box
0752       if (chanDigi >= 0 && (*i).second > 0) {
0753         theNewDigiSignal[chanDigi] += SiPixelDigitizerAlgorithm::Amplitude((*i).second, (*i).second);
0754       }
0755     }
0756   }
0757 
0758   return true;
0759 }