Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-11 16:23:22

0001 #include <iostream>
0002 #include <cmath>
0003 
0004 #include "SimTracker/SiPhase2Digitizer/plugins/PixelBrickedDigitizerAlgorithm.h"
0005 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0006 
0007 #include "FWCore/Framework/interface/ConsumesCollector.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
0012 
0013 #include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
0014 #include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
0015 #include "CondFormats/DataRecord/interface/SiPixelFedCablingMapRcd.h"
0016 #include "CondFormats/DataRecord/interface/SiPixelLorentzAngleSimRcd.h"
0017 
0018 // Geometry
0019 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0020 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0021 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0022 
0023 using namespace edm;
0024 using namespace sipixelobjects;
0025 
0026 PixelBrickedDigitizerAlgorithm::PixelBrickedDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0027     : PixelDigitizerAlgorithm(conf, iC)
0028 /*    odd_row_interchannelCoupling_next_row_(conf.getParameter<ParameterSet>("PixelBrickedDigitizerAlgorithm")
0029                                                  .getParameter<double>("Odd_row_interchannelCoupling_next_row"))
0030          even_row_interchannelCoupling_next_row_(conf.getParameter<ParameterSet>("PixelBrickedDigitizerAlgorithm")
0031                                                   .getParameter<double>("Even_row_interchannelCoupling_next_row")),
0032       odd_column_interchannelCoupling_next_column_(
0033           conf.getParameter<ParameterSet>("PixelBrickedDigitizerAlgorithm")
0034               .getParameter<double>("Odd_column_interchannelCoupling_next_column")),
0035       even_column_interchannelCoupling_next_column_(
0036           conf.getParameter<ParameterSet>("PixelBrickedDigitizerAlgorithm")
0037       .getParameter<double>("Even_column_interchannelCoupling_next_column"))*/
0038 {
0039   even_row_interchannelCoupling_next_row_ = conf.getParameter<ParameterSet>("PixelBrickedDigitizerAlgorithm")
0040                                                 .getParameter<double>("Even_row_interchannelCoupling_next_row");
0041   pixelFlag_ = true;
0042   LogDebug("PixelBrickedDigitizerAlgorithm")
0043       << "Algorithm constructed "
0044       << "Configuration parameters:"
0045       << "Threshold/Gain = "
0046       << "threshold in electron Endcap = " << theThresholdInE_Endcap_
0047       << "threshold in electron Barrel = " << theThresholdInE_Barrel_ << " " << theElectronPerADC_ << " "
0048       << theAdcFullScale_ << " The delta cut-off is set to " << tMax_ << " pix-inefficiency " << addPixelInefficiency_;
0049 }
0050 PixelBrickedDigitizerAlgorithm::~PixelBrickedDigitizerAlgorithm() {
0051   LogDebug("PixelBrickedDigitizerAlgorithm") << "Algorithm deleted";
0052 }
0053 void PixelBrickedDigitizerAlgorithm::induce_signal(const PSimHit& hit,
0054                                                    const size_t hitIndex,
0055                                                    const uint32_t tofBin,
0056                                                    const Phase2TrackerGeomDetUnit* pixdet,
0057                                                    const std::vector<DigitizerUtility::SignalPoint>& collection_points) {
0058   // X  - Rows, Left-Right, 160, (1.6cm)   for barrel
0059   // Y  - Columns, Down-Up, 416, (6.4cm)
0060   const Phase2TrackerTopology* topol = &pixdet->specificTopology();
0061   uint32_t detID = pixdet->geographicalId().rawId();
0062   signal_map_type& theSignal = _signal[detID];
0063 
0064   // local map to store pixels hit by 1 Hit.
0065   using hit_map_type = std::map<int, float, std::less<int> >;
0066   hit_map_type hit_signal;
0067 
0068   // Assign signals to readout channels and store sorted by channel number
0069   // Iterate over collection points on the collection plane
0070   for (auto const& v : collection_points) {
0071     float CloudCenterX = v.position().x();  // Charge position in x
0072     float CloudCenterY = v.position().y();  //                 in y
0073     float SigmaX = v.sigma_x();             // Charge spread in x
0074     float SigmaY = v.sigma_y();             //               in y
0075     float Charge = v.amplitude();           // Charge amplitude
0076 
0077     // Find the maximum cloud spread in 2D plane , assume 3*sigma
0078     float CloudRight = CloudCenterX + clusterWidth_ * SigmaX;
0079     float CloudLeft = CloudCenterX - clusterWidth_ * SigmaX;
0080     float CloudUp = CloudCenterY + clusterWidth_ * SigmaY;
0081     float CloudDown = CloudCenterY - clusterWidth_ * SigmaY;
0082 
0083     // Define 2D cloud limit points
0084     LocalPoint PointRightUp = LocalPoint(CloudRight, CloudUp);
0085     LocalPoint PointLeftDown = LocalPoint(CloudLeft, CloudDown);
0086 
0087     // This points can be located outside the sensor area.
0088     // The conversion to measurement point does not check for that
0089     // so the returned pixel index might be wrong (outside range).
0090     // We rely on the limits check below to fix this.
0091     // But remember whatever we do here THE CHARGE OUTSIDE THE ACTIVE
0092     // PIXEL AREA IS LOST, it should not be collected.
0093 
0094     // Convert the 2D points to pixel indices
0095     MeasurementPoint mp = topol->measurementPosition(PointRightUp);
0096     //MeasurementPoint mp_bricked = topol->measurementPosition(PointRightUp);
0097     int IPixRightUpX = static_cast<int>(std::floor(mp.x()));  // cast reqd.
0098     //int IPixRightUpY = static_cast<int>(std::floor(mp.y()));
0099 
0100     int numColumns = topol->ncolumns();  // det module number of cols&rows
0101     int numRows = topol->nrows();
0102     IPixRightUpX = numRows > IPixRightUpX ? IPixRightUpX : numRows - 1;
0103 
0104     //Specific to bricked geometry
0105     int IPixRightUpY = static_cast<int>(mp.y() - 0.5 * (IPixRightUpX % 2));
0106 
0107     mp = topol->measurementPosition(PointLeftDown);
0108 
0109     int IPixLeftDownX = static_cast<int>(std::floor(mp.x()));
0110 
0111     IPixLeftDownX = 0 < IPixLeftDownX ? IPixLeftDownX : 0;
0112 
0113     //Specific to bricked geometry
0114     int IPixLeftDownY = static_cast<int>(mp.y() - 0.5 * (IPixLeftDownX % 2));  //changed in case negative value
0115 
0116     IPixRightUpY = numColumns > IPixRightUpY ? IPixRightUpY : numColumns - 1;
0117     IPixLeftDownY = 0 < IPixLeftDownY ? IPixLeftDownY : 0;
0118 
0119     // First integrate charge strips in x
0120     hit_map_type x;
0121     for (int ix = IPixLeftDownX; ix <= IPixRightUpX; ++ix) {  // loop over x index
0122       float xLB, LowerBound;
0123       // Why is set to 0 if ix=0, does it meen that we accept charge
0124       // outside the sensor?
0125       if (ix == 0 || SigmaX == 0.) {  // skip for surface segemnts
0126         LowerBound = 0.;
0127       } else {
0128         mp = MeasurementPoint(ix, 0.0);
0129         xLB = topol->localPosition(mp).x();
0130         LowerBound = 1 - calcQ((xLB - CloudCenterX) / SigmaX);
0131       }
0132 
0133       float xUB, UpperBound;
0134       if (ix == numRows - 1 || SigmaX == 0.) {
0135         UpperBound = 1.;
0136       } else {
0137         mp = MeasurementPoint(ix + 1, 0.0);
0138         xUB = topol->localPosition(mp).x();
0139         UpperBound = 1. - calcQ((xUB - CloudCenterX) / SigmaX);
0140       }
0141       float TotalIntegrationRange = UpperBound - LowerBound;  // get strip
0142       x.emplace(ix, TotalIntegrationRange);                   // save strip integral
0143     }
0144 
0145     // Now integrate strips in y. Two maps will be filled: y and y_bricked which will both be used for the induced signal.
0146     int IPixLeftDownY_bricked = IPixLeftDownY;
0147     int IPixRightUpY_bricked = IPixRightUpY;
0148 
0149     //Specific to bricked geometry
0150     IPixRightUpY = std::min(IPixRightUpY + int((IPixRightUpX % 2)), numColumns - 1);
0151 
0152     //This map will be twice as large as the non-bricked hit map in y to harbor both the integrated charge from the bricked and non-bricked columns.
0153     hit_map_type y;
0154     for (int iy = IPixLeftDownY; iy <= IPixRightUpY; ++iy) {  // loop over y index
0155       float yLB, LowerBound;
0156       if (iy == 0 || SigmaY == 0.) {
0157         LowerBound = 0.;
0158       } else {
0159         mp = MeasurementPoint(0.0, iy);
0160         yLB = topol->localPosition(mp).y();
0161         LowerBound = 1. - calcQ((yLB - CloudCenterY) / SigmaY);
0162       }
0163 
0164       float yUB, UpperBound;
0165       if (iy >= numColumns - 1 || SigmaY == 0.) {
0166         UpperBound = 1.;
0167       } else {
0168         mp = MeasurementPoint(0.0, iy + 1);
0169         yUB = topol->localPosition(mp).y();
0170         UpperBound = 1. - calcQ((yUB - CloudCenterY) / SigmaY);
0171       }
0172 
0173       float TotalIntegrationRange = UpperBound - LowerBound;
0174 
0175       //Even indices correspond to the non-bricked columns
0176       y.emplace(2 * iy, TotalIntegrationRange);  // save strip integral
0177     }
0178 
0179     IPixLeftDownY_bricked = std::max(IPixLeftDownY_bricked - int((!(IPixLeftDownX % 2))), 0);
0180     for (int iy = IPixLeftDownY_bricked; iy <= IPixRightUpY_bricked; ++iy) {  // loop over y index
0181       float yLB, LowerBound;
0182       if (iy == 0 || SigmaY == 0.) {
0183         LowerBound = 0.;
0184       } else {
0185         mp = MeasurementPoint(0.0, iy + 0.5);
0186         yLB = topol->localPosition(mp).y();
0187         LowerBound = 1. - calcQ((yLB - CloudCenterY) / SigmaY);
0188       }
0189 
0190       float yUB, UpperBound;
0191       if (iy >= numColumns || SigmaY == 0.) {  // This was changed for bricked pixels
0192         UpperBound = 1.;
0193       } else {
0194         mp = MeasurementPoint(0.0, iy + 1.5);
0195         yUB = topol->localPosition(mp).y();
0196         UpperBound = 1. - calcQ((yUB - CloudCenterY) / SigmaY);
0197       }
0198 
0199       float TotalIntegrationRange = UpperBound - LowerBound;
0200       //Odd indices correspond to bricked columns
0201       y.emplace(2 * iy + 1, TotalIntegrationRange);  // save strip integral
0202     }                                                //loop over y index
0203 
0204     // Get the 2D charge integrals by folding x and y strips
0205     for (int ix = IPixLeftDownX; ix <= IPixRightUpX; ++ix) {                                 // loop over x index
0206       for (int iy = std::max(0, IPixLeftDownY - int((ix % 2))); iy <= IPixRightUpY; ++iy) {  // loop over y index
0207         int iy_considered = iy * 2 + ix % 2;
0208         float ChargeFraction = Charge * x[ix] * y[iy_considered];
0209 
0210         int chanFired = -1;
0211         if (ChargeFraction > 0.) {
0212           chanFired =
0213               pixelFlag_ ? PixelDigi::pixelToChannel(ix, iy) : Phase2TrackerDigi::pixelToChannel(ix, iy);  // Get index
0214           // Load the ampl
0215           hit_signal[chanFired] += ChargeFraction;
0216         }
0217       }
0218     }  //x loop
0219   }    //collection loop
0220   // Fill the global map with all hit pixels from this event
0221   float corr_time = hit.tof() - pixdet->surface().toGlobal(hit.localPosition()).mag() * c_inv;
0222   for (auto const& hit_s : hit_signal) {
0223     int chan = hit_s.first;
0224     theSignal[chan] +=
0225         (makeDigiSimLinks_ ? DigitizerUtility::Amplitude(hit_s.second, &hit, hit_s.second, corr_time, hitIndex, tofBin)
0226                            : DigitizerUtility::Amplitude(hit_s.second, nullptr, hit_s.second));
0227   }
0228 }