Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-02-05 23:51:49 

0001 #ifndef SimTracker_Common_SiPixelChargeReweightingAlgorithm_h
0002 #define SimTracker_Common_SiPixelChargeReweightingAlgorithm_h
0003 
0004 // Original Author Caroline Collard
0005 // September 2020 : Extraction of the code for cluster charge reweighting from SiPixelDigitizerAlgorithm to a new class
0006 // Jan 2023: Generalize so it can be used for Phase-2 IT as well (M. Musich, T.A. Vami)
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 // forward declarations
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   // initialization that cannot be done in the constructor
0057   void init(const edm::EventSetup& es);
0058 
0059   // template this so it can be used for Phase-2 as well
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);  //  could be templated in future ?
0087 private:
0088   // Internal typedef
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   // Variables and objects for the charge reweighting using 2D templates
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   // methods for charge reweighting in irradiated sensors
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   // Read template files for charge reweighting
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     // The clusters do not have an overlap, hence the hit is NOT reweighted
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;  //Middle of sensor is origin for Z-axis
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   // define loop boundaries that will prevent the row and col loops
0311   // from going out of physical bounds of the pixel module
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       //Fill charges into 21x13 Pixel Array with hitPixel in centre
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   // for unirradiated: 2nd argument is IDden
0332   // for irradiated: 2nd argument is IDnum
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 //! Reweight CMSSW clusters to look like clusters corresponding to Pixelav Templates.
0387 //! \param       id_in - (input) identifier of the template corresponding to the input events
0388 //! \param    id_rewgt - (input) identifier of the template corresponding to the output events
0389 //! \param     cluster - (input/output) boost multi_array container of 7x21 array of pixel signals,
0390 //!                       origin of local coords (0,0) at center of pixel cluster[3][10].
0391 //! returns 0 if everything is OK, 1 if angles are outside template coverage (cluster is probably still
0392 //! usable, > 1 if something is wrong (no reweight done).
0393 // *******************************************************************************************************
0394 inline int SiPixelChargeReweightingAlgorithm::PixelTempRewgt2D(int id_in, int id_rewgt, array_2d& cluster) {
0395   // Local variables
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   // success = 0 is returned if everthing is OK
0404   success = 0;
0405 
0406   // Copy the array to remember original charges
0407   array_2d clust(cluster);
0408 
0409   // Take the pixel dimensions from the 2D template
0410   templ2D.getid(id_in);
0411   xsize = templ2D.xsize();
0412   ysize = templ2D.ysize();
0413 
0414   // Calculate the track angles
0415 
0416   if (std::abs(track[5]) > 0.f) {
0417     cotalpha = track[3] / track[5];  //if track[5] (direction in z) is 0 the hit is not processed by re-weighting
0418     cotbeta = track[4] / track[5];
0419   } else {
0420     LogDebug("SiPixelChargeReweightingAlgorithm") << "Reweighting angle is not good! \n";
0421     return 9;  //returned value here indicates that no reweighting was done in this case
0422   }
0423 
0424   // The 2-D templates are defined on a shifted coordinate system wrt the 1D templates
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   // Zero the input and output templates
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   // Next, interpolate the CMSSW template needed to analyze this cluster
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   //q50i = 0;
0461   q100i = 2.f * q50i;
0462 
0463   // save input charge scale (14/4/2023)
0464   qscalei = templ2D.qscale();
0465 
0466   // Check that the cluster container is a 13x21 matrix
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   // Sum initial charge in the cluster
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   // Next, interpolate the physical output template needed to reweight
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   // calculate ratio of charge scaling factors (14/4/2023)
0510   rqscale = qscalei / templ2D.qscale();
0511 
0512   // Find all non-zero denominator pixels in the input template and generate "inside" weights
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   // Find all non-zero numerator pixels not matched to denominator in the output template and generate "inside" weights
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         // Search for nearest denominator pixel
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   // Now reweight pixels outside of template footprint using closest weights
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   // final rescaling by the ratio of charge scaling factors (14/4/2023)
0598   // put this here to avoid changing the threshold tests above and to be vectorizable
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 }  // PixelTempRewgt2D
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   //loop on digis
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   // end loop on digis
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   // addition as now the hit himself is not available
0700   // https://github.com/cms-sw/cmssw/blob/master/SimDataFormats/TrackingHit/interface/PSimHit.h#L52
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     // The clusters do not have an overlap, hence the hit is NOT reweighted
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;  //Middle of sensor is origin for Z-axis
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   for (int row = 0; row < TXSIZE; ++row) {
0749     for (int col = 0; col < TYSIZE; ++col) {
0750       pixrewgt[row][col] = 0;
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       //Fill charges into 21x13 Pixel Array with hitPixel in centre
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   // for unirradiated: 2nd argument is IDden
0778   // for irradiated: 2nd argument is IDnum
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   // need to store the digi out of the 21x13 box.
0826   for (const auto& i : theDigiSignal) {
0827     // check if in the 21x13 box
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       // not in the box
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   //loop on digis
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   // end loop on digis
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   // addition as now the hit himself is not available
0916   // https://github.com/cms-sw/cmssw/blob/master/SimDataFormats/TrackingHit/interface/PSimHit.h#L52
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     // The clusters do not have an overlap, hence the hit is NOT reweighted
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;  //Middle of sensor is origin for Z-axis
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   for (int row = 0; row < TXSIZE; ++row) {
0965     for (int col = 0; col < TYSIZE; ++col) {
0966       pixrewgt[row][col] = 0;
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       //Fill charges into 21x13 Pixel Array with hitPixel in centre
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   // for unirradiated: 2nd argument is IDden
0994   // for irradiated: 2nd argument is IDnum
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   // need to store the digi out of the 21x13 box.
1042   for (const auto& i : theDigiSignal) {
1043     // check if in the 21x13 box
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       // not in the box
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
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!