Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:25:35 

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 xy_in[BXM2][BYM2], xy_rewgt[BXM2][BYM2], xy_clust[TXSIZE][TYSIZE];
0393   int denx_clust[TXSIZE][TYSIZE], deny_clust[TXSIZE][TYSIZE];
0394   int goodWeightsUsed, nearbyWeightsUsed, noWeightsUsed;
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   // Check that the cluster container is a 13x21 matrix
0457 
0458   if (cluster.num_dimensions() != 2) {
0459     edm::LogWarning("SiPixelChargeReweightingAlgorithm") << "Cluster is not 2-dimensional. Return. \n";
0460     return 3;
0461   }
0462   nclusx = (int)cluster.shape()[0];
0463   nclusy = (int)cluster.shape()[1];
0464   if (nclusx != TXSIZE || xdouble.size() != TXSIZE) {
0465     edm::LogWarning("SiPixelChargeReweightingAlgorithm")
0466         << "Sizes in x do not match: nclusx=" << nclusx << "  xdoubleSize=" << xdouble.size() << "  TXSIZE=" << TXSIZE
0467         << ". Return. \n";
0468     return 4;
0469   }
0470   if (nclusy != TYSIZE || ydouble.size() != TYSIZE) {
0471     edm::LogWarning("SiPixelChargeReweightingAlgorithm") << "Sizes in y do not match. Return. \n";
0472     return 5;
0473   }
0474 
0475   // Sum initial charge in the cluster
0476 
0477   for (i = 0; i < TYSIZE; ++i) {
0478     for (j = 0; j < TXSIZE; ++j) {
0479       xy_clust[j][i] = 0.f;
0480       denx_clust[j][i] = 0;
0481       deny_clust[j][i] = 0;
0482     }
0483   }
0484 
0485   // Next, interpolate the physical output template needed to reweight
0486 
0487   if (!templ2D.xytemp(id_rewgt, cotalpha, cotbeta, xhit2D, yhit2D, ydouble, xdouble, xy_rewgt)) {
0488     success = 1;
0489   }
0490 
0491   if (PrintTemplates) {
0492     LogDebug("SiPixelChargeReweightingAlgorithm") << "Template irrad: \n";
0493     printCluster(xy_rewgt);
0494   }
0495 
0496   q50r = templ2D.s50();
0497   q10r = 0.2f * q50r;
0498 
0499   // Find all non-zero denominator pixels in the input template and generate "inside" weights
0500 
0501   int ntpix = 0;
0502   int ncpix = 0;
0503   std::vector<int> ytclust;
0504   std::vector<int> xtclust;
0505   std::vector<int> ycclust;
0506   std::vector<int> xcclust;
0507   for (i = 0; i < TYSIZE; ++i) {
0508     for (j = 0; j < TXSIZE; ++j) {
0509       if (xy_in[j + 1][i + 1] > q100i) {
0510         ++ntpix;
0511         ytclust.push_back(i);
0512         xtclust.push_back(j);
0513         xy_clust[j][i] = xy_rewgt[j + 1][i + 1] / xy_in[j + 1][i + 1];
0514         denx_clust[j][i] = j;
0515         deny_clust[j][i] = i;
0516       }
0517     }
0518   }
0519 
0520   // Find all non-zero numerator pixels not matched to denominator in the output template and generate "inside" weights
0521 
0522   for (i = 0; i < TYSIZE; ++i) {
0523     for (j = 0; j < TXSIZE; ++j) {
0524       if (xy_rewgt[j + 1][i + 1] > q10r && xy_clust[j][i] == 0.f && ntpix > 0) {
0525         // Search for nearest denominator pixel
0526         dmin2 = 10000.f;
0527         kclose = 0;
0528         for (k = 0; k < ntpix; ++k) {
0529           dist2 = (i - ytclust[k]) * (i - ytclust[k]) + 0.44444f * (j - xtclust[k]) * (j - xtclust[k]);
0530           if (dist2 < dmin2) {
0531             dmin2 = dist2;
0532             kclose = k;
0533           }
0534         }
0535         xy_clust[j][i] = xy_rewgt[j + 1][i + 1] / xy_in[xtclust[kclose] + 1][ytclust[kclose] + 1];
0536         denx_clust[j][i] = xtclust[kclose];
0537         deny_clust[j][i] = ytclust[kclose];
0538       }
0539     }
0540   }
0541 
0542   if (PrintTemplates) {
0543     LogDebug("SiPixelChargeReweightingAlgorithm") << "Weights: \n";
0544     printCluster(xy_clust);
0545   }
0546 
0547   // Do the reweighting
0548   goodWeightsUsed = 0;
0549   nearbyWeightsUsed = 0;
0550   noWeightsUsed = 0;
0551 
0552   for (i = 0; i < TYSIZE; ++i) {
0553     for (j = 0; j < TXSIZE; ++j) {
0554       if (xy_clust[j][i] > 0.f) {
0555         cluster[j][i] = xy_clust[j][i] * clust[denx_clust[j][i]][deny_clust[j][i]];
0556         if (cluster[j][i] > 0) {
0557           goodWeightsUsed++;
0558         }
0559       } else {
0560         if (clust[j][i] > 0.f) {
0561           ++ncpix;
0562           ycclust.push_back(i);
0563           xcclust.push_back(j);
0564         }
0565       }
0566     }
0567   }
0568 
0569   // Now reweight pixels outside of template footprint using closest weights
0570 
0571   if (ncpix > 0) {
0572     for (l = 0; l < ncpix; ++l) {
0573       i = ycclust[l];
0574       j = xcclust[l];
0575       dmin2 = 10000.f;
0576       kclose = 0;
0577       for (k = 0; k < ntpix; ++k) {
0578         dist2 = (i - ytclust[k]) * (i - ytclust[k]) + 0.44444f * (j - xtclust[k]) * (j - xtclust[k]);
0579         if (dist2 < dmin2) {
0580           dmin2 = dist2;
0581           kclose = k;
0582         }
0583       }
0584       if (dmin2 < 5.f) {
0585         nearbyWeightsUsed++;
0586         cluster[j][i] *= xy_clust[xtclust[kclose]][ytclust[kclose]];
0587       } else {
0588         noWeightsUsed++;
0589         cluster[j][i] = 0.f;
0590       }
0591     }
0592   }
0593 
0594   return success;
0595 }  // PixelTempRewgt2D
0596 
0597 inline void SiPixelChargeReweightingAlgorithm::printCluster(array_2d& cluster) {
0598   for (int col = 0; col < TYSIZE; ++col) {
0599     for (int row = 0; row < TXSIZE; ++row) {
0600       LogDebug("SiPixelChargeReweightingAlgorithm") << cluster[row][col];
0601     }
0602     LogDebug("SiPixelChargeReweightingAlgorithm") << "\n";
0603   }
0604 }
0605 
0606 inline void SiPixelChargeReweightingAlgorithm::printCluster(float arr[BXM2][BYM2]) {
0607   for (int col = 0; col < BYM2; ++col) {
0608     for (int row = 0; row < BXM2; ++row) {
0609       LogDebug("SiPixelChargeReweightingAlgorithm") << arr[row][col];
0610     }
0611     LogDebug("SiPixelChargeReweightingAlgorithm") << "\n";
0612   }
0613 }
0614 
0615 inline void SiPixelChargeReweightingAlgorithm::printCluster(float arr[TXSIZE][TYSIZE]) {
0616   for (int col = 0; col < TYSIZE; ++col) {
0617     for (int row = 0; row < TXSIZE; ++row) {
0618       LogDebug("SiPixelChargeReweightingAlgorithm") << arr[row][col];
0619     }
0620     LogDebug("SiPixelChargeReweightingAlgorithm") << "\n";
0621   }
0622 }
0623 
0624 template <class AmplitudeType, typename SignalType>
0625 bool SiPixelChargeReweightingAlgorithm::lateSignalReweight(const PixelGeomDetUnit* pixdet,
0626                                                            std::vector<PixelDigi>& digis,
0627                                                            PixelSimHitExtraInfo& loopTempSH,
0628                                                            SignalType& theNewDigiSignal,
0629                                                            const TrackerTopology* tTopo,
0630                                                            CLHEP::HepRandomEngine* engine) {
0631   uint32_t detID = pixdet->geographicalId().rawId();
0632   const PixelTopology* topol = &pixdet->specificTopology();
0633 
0634   if (UseReweighting) {
0635     edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0636     edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0637     edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0638     edm::LogError("SiPixelChargeReweightingAlgorithm") << " *****  INCONSISTENCY !!!   ***** \n";
0639     edm::LogError("SiPixelChargeReweightingAlgorithm")
0640         << " applyLateReweighting_ and UseReweighting can not be true at the same time for PU ! \n";
0641     edm::LogError("SiPixelChargeReweightingAlgorithm") << " ---> DO NOT APPLY CHARGE REWEIGHTING TWICE !!! \n";
0642     edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0643     edm::LogError("SiPixelChargeReweightingAlgorithm") << " ******************************** \n";
0644     return false;
0645   }
0646 
0647   SignalType theDigiSignal;
0648 
0649   int irow_min = topol->nrows();
0650   int irow_max = 0;
0651   int icol_min = topol->ncolumns();
0652   int icol_max = 0;
0653 
0654   float chargeBefore = 0;
0655   float chargeAfter = 0;
0656 
0657   //loop on digis
0658   std::vector<PixelDigi>::const_iterator loopDigi;
0659   for (loopDigi = digis.begin(); loopDigi != digis.end(); ++loopDigi) {
0660     unsigned int chan = loopDigi->channel();
0661     if (loopTempSH.isInTheList(chan)) {
0662       float corresponding_charge = loopDigi->adc();
0663       if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0664         edm::LogError("SiPixelChargeReweightingAlgorithm")
0665             << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
0666         return false;
0667       } else {
0668         theDigiSignal[chan] += AmplitudeType(corresponding_charge, corresponding_charge);
0669       }
0670       chargeBefore += corresponding_charge;
0671       if (loopDigi->row() < irow_min)
0672         irow_min = loopDigi->row();
0673       if (loopDigi->row() > irow_max)
0674         irow_max = loopDigi->row();
0675       if (loopDigi->column() < icol_min)
0676         icol_min = loopDigi->column();
0677       if (loopDigi->column() > icol_max)
0678         icol_max = loopDigi->column();
0679     }
0680   }
0681   // end loop on digis
0682 
0683   LocalVector direction = loopTempSH.exitPoint() - loopTempSH.entryPoint();
0684   LocalPoint hitEntryPoint = loopTempSH.entryPoint();
0685   float trajectoryScaleToPosition = std::abs(hitEntryPoint.z() / direction.z());
0686   LocalPoint hitPosition = hitEntryPoint + trajectoryScaleToPosition * direction;
0687 
0688   // addition as now the hit himself is not available
0689   // https://github.com/cms-sw/cmssw/blob/master/SimDataFormats/TrackingHit/interface/PSimHit.h#L52
0690   LocalPoint hitLocalPosition = hitEntryPoint + 0.5 * direction;
0691   MeasurementPoint hitPositionPixel = topol->measurementPosition(hitLocalPosition);
0692   std::pair<int, int> hitPixel =
0693       std::pair<int, int>(int(floor(hitPositionPixel.x())), int(floor(hitPositionPixel.y())));
0694 
0695   MeasurementPoint originPixel = MeasurementPoint(hitPixel.first - THX + 0.5, hitPixel.second - THY + 0.5);
0696   LocalPoint origin = topol->localPosition(originPixel);
0697 
0698   MeasurementPoint hitEntryPointPixel = topol->measurementPosition(loopTempSH.entryPoint());
0699   MeasurementPoint hitExitPointPixel = topol->measurementPosition(loopTempSH.exitPoint());
0700   std::pair<int, int> entryPixel =
0701       std::pair<int, int>(int(floor(hitEntryPointPixel.x())), int(floor(hitEntryPointPixel.y())));
0702   std::pair<int, int> exitPixel =
0703       std::pair<int, int>(int(floor(hitExitPointPixel.x())), int(floor(hitExitPointPixel.y())));
0704 
0705   int hitcol_min, hitcol_max, hitrow_min, hitrow_max;
0706   if (entryPixel.first > exitPixel.first) {
0707     hitrow_min = exitPixel.first;
0708     hitrow_max = entryPixel.first;
0709   } else {
0710     hitrow_min = entryPixel.first;
0711     hitrow_max = exitPixel.first;
0712   }
0713 
0714   if (entryPixel.second > exitPixel.second) {
0715     hitcol_min = exitPixel.second;
0716     hitcol_max = entryPixel.second;
0717   } else {
0718     hitcol_min = entryPixel.second;
0719     hitcol_max = exitPixel.second;
0720   }
0721 
0722   if (!(irow_min <= hitrow_max && irow_max >= hitrow_min && icol_min <= hitcol_max && icol_max >= hitcol_min)) {
0723     // The clusters do not have an overlap, hence the hit is NOT reweighted
0724     return false;
0725   }
0726 
0727   track[0] = (hitPosition.x() - origin.x()) * cmToMicrons;
0728   track[1] = (hitPosition.y() - origin.y()) * cmToMicrons;
0729   track[2] = 0.0f;  //Middle of sensor is origin for Z-axis
0730   track[3] = direction.x();
0731   track[4] = direction.y();
0732   track[5] = direction.z();
0733 
0734   array_2d pixrewgt(boost::extents[TXSIZE][TYSIZE]);
0735 
0736   /*
0737   for (int row = 0; row < TXSIZE; ++row) {
0738     for (int col = 0; col < TYSIZE; ++col) {
0739       pixrewgt[row][col] = 0;
0740     }
0741   }
0742 */
0743 
0744   for (int row = 0; row < TXSIZE; ++row) {
0745     xdouble[row] = topol->isItBigPixelInX(hitPixel.first + row - THX);
0746   }
0747 
0748   for (int col = 0; col < TYSIZE; ++col) {
0749     ydouble[col] = topol->isItBigPixelInY(hitPixel.second + col - THY);
0750   }
0751 
0752   for (int row = 0; row < TXSIZE; ++row) {
0753     for (int col = 0; col < TYSIZE; ++col) {
0754       //Fill charges into 21x13 Pixel Array with hitPixel in centre
0755       pixrewgt[row][col] =
0756           theDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)];
0757     }
0758   }
0759 
0760   if (PrintClusters) {
0761     LogDebug("SiPixelChargeReweightingAlgorithm") << " Cluster before reweighting: ";
0762     printCluster(pixrewgt);
0763   }
0764 
0765   int ierr;
0766   // for unirradiated: 2nd argument is IDden
0767   // for irradiated: 2nd argument is IDnum
0768   int ID1 = dbobject_num->getTemplateID(detID);
0769   int ID0 = dbobject_den->getTemplateID(detID);
0770 
0771   if (ID0 == ID1) {
0772     LogDebug("SiPixelChargeReweightingAlgorithm") << " same template for num and den ";
0773     return false;
0774   }
0775   ierr = PixelTempRewgt2D(ID0, ID1, pixrewgt);
0776   if (ierr != 0) {
0777     edm::LogError("SiPixelChargeReweightingAlgorithm") << "Cluster Charge Reweighting did not work properly.";
0778     return false;
0779   }
0780   if (PrintClusters) {
0781     LogDebug("SiPixelChargeReweightingAlgorithm") << " Cluster after reweighting: ";
0782     printCluster(pixrewgt);
0783   }
0784 
0785   for (int row = 0; row < TXSIZE; ++row) {
0786     for (int col = 0; col < TYSIZE; ++col) {
0787       float charge = 0;
0788       charge = pixrewgt[row][col];
0789       if ((hitPixel.first + row - THX) >= 0 && (hitPixel.first + row - THX) < topol->nrows() &&
0790           (hitPixel.second + col - THY) >= 0 && (hitPixel.second + col - THY) < topol->ncolumns() && charge > 0) {
0791         chargeAfter += charge;
0792         if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0793           edm::LogError("SiPixelChargeReweightingAlgorithm")
0794               << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
0795           return false;
0796         } else {
0797           theNewDigiSignal[PixelDigi::pixelToChannel(hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
0798               AmplitudeType(charge, charge);
0799         }
0800       }
0801     }
0802   }
0803 
0804   if (chargeBefore != 0. && chargeAfter == 0.) {
0805     return false;
0806   }
0807 
0808   if (PrintClusters) {
0809     LogDebug("SiPixelChargeReweightingAlgorithm")
0810         << "Charges (before->after): " << chargeBefore << " -> " << chargeAfter;
0811     LogDebug("SiPixelChargeReweightingAlgorithm") << "Charge loss: " << (1 - chargeAfter / chargeBefore) * 100 << " %";
0812   }
0813 
0814   // need to store the digi out of the 21x13 box.
0815   for (auto& i : theDigiSignal) {
0816     // check if in the 21x13 box
0817     int chanDigi = i.first;
0818     std::pair<int, int> ip = PixelDigi::channelToPixel(chanDigi);
0819     int row_digi = ip.first;
0820     int col_digi = ip.second;
0821     int i_transformed_row = row_digi - hitPixel.first + THX;
0822     int i_transformed_col = col_digi - hitPixel.second + THY;
0823     if (i_transformed_row < 0 || i_transformed_row > TXSIZE || i_transformed_col < 0 || i_transformed_col > TYSIZE) {
0824       // not in the box
0825       if (chanDigi >= 0 && i.second > 0) {
0826         if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
0827           edm::LogError("SiPixelChargeReweightingAlgorithm")
0828               << "Phase-2 late charge reweighting not supported (not sure we need it at all)";
0829           return false;
0830         } else {
0831           theNewDigiSignal[chanDigi] += AmplitudeType(i.second, i.second);
0832         }
0833       }
0834     }
0835   }
0836 
0837   return true;
0838 }
0839 
0840 #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!