Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:04:59 

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