Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-18 22:24:07

0001 #include "SimTracker/SiPhase2Digitizer/plugins/Pixel3DDigitizerAlgorithm.h"
0002 
0003 // Framework infrastructure
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 
0008 // Calibration & Conditions
0009 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
0010 
0011 // Geometry
0012 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0013 
0014 //#include <iostream>
0015 #include <cmath>
0016 #include <vector>
0017 #include <algorithm>
0018 
0019 using namespace sipixelobjects;
0020 
0021 namespace {
0022   // Analogously to CMSUnits (no um defined)
0023   constexpr double operator""_um(long double length) { return length * 1e-4; }
0024   constexpr double operator""_um_inv(long double length) { return length * 1e4; }
0025 }  // namespace
0026 
0027 Pixel3DDigitizerAlgorithm::Pixel3DDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0028     : PixelDigitizerAlgorithm(conf, iC),
0029       np_column_radius_(
0030           (conf.getParameter<edm::ParameterSet>("Pixel3DDigitizerAlgorithm").getParameter<double>("NPColumnRadius")) *
0031           1.0_um),
0032       ohm_column_radius_(
0033           (conf.getParameter<edm::ParameterSet>("Pixel3DDigitizerAlgorithm").getParameter<double>("OhmicColumnRadius")) *
0034           1.0_um),
0035       np_column_gap_(
0036           (conf.getParameter<edm::ParameterSet>("Pixel3DDigitizerAlgorithm").getParameter<double>("NPColumnGap")) *
0037           1.0_um) {
0038   // XXX - NEEDED?
0039   pixelFlag_ = true;
0040 
0041   edm::LogInfo("Pixel3DDigitizerAlgorithm")
0042       << "Algorithm constructed \n"
0043       << "Configuration parameters:\n"
0044       << "\n*** Threshold"
0045       << "\n    Endcap = " << theThresholdInE_Endcap_ << " electrons"
0046       << "\n    Barrel = " << theThresholdInE_Barrel_ << " electrons"
0047       << "\n*** Gain"
0048       << "\n    Electrons per ADC:" << theElectronPerADC_ << "\n    ADC Full Scale: " << theAdcFullScale_
0049       << "\n*** The delta cut-off is set to " << tMax_ << "\n*** Pixel-inefficiency: " << addPixelInefficiency_;
0050 }
0051 
0052 Pixel3DDigitizerAlgorithm::~Pixel3DDigitizerAlgorithm() {}
0053 
0054 //
0055 // -- Select the Hit for Digitization
0056 //
0057 bool Pixel3DDigitizerAlgorithm::select_hit(const PSimHit& hit, double tCorr, double& sigScale) const {
0058   double time = hit.tof() - tCorr;
0059   return (time >= theTofLowerCut_ && time < theTofUpperCut_);
0060 }
0061 
0062 const bool Pixel3DDigitizerAlgorithm::is_inside_n_column_(const LocalPoint& p, const float& sensor_thickness) const {
0063   // The insensitive volume of the column: sensor thickness - column gap distance
0064   return (p.perp() <= np_column_radius_ && p.z() <= (sensor_thickness - np_column_gap_));
0065 }
0066 
0067 const bool Pixel3DDigitizerAlgorithm::is_inside_ohmic_column_(const LocalPoint& p,
0068                                                               const std::pair<float, float>& half_pitch) const {
0069   // The four corners of the cell
0070   return ((p - LocalVector(half_pitch.first, half_pitch.second, 0)).perp() <= ohm_column_radius_) ||
0071          ((p - LocalVector(-half_pitch.first, half_pitch.second, 0)).perp() <= ohm_column_radius_) ||
0072          ((p - LocalVector(half_pitch.first, -half_pitch.second, 0)).perp() <= ohm_column_radius_) ||
0073          ((p - LocalVector(-half_pitch.first, -half_pitch.second, 0)).perp() <= ohm_column_radius_);
0074 }
0075 
0076 // Diffusion algorithm: Probably not needed,
0077 // Assuming the position point is given in the reference system of the proxy
0078 // cell, centered at the n-column.
0079 // The algorithm assumes only 1-axis could produce the charge migration, this assumption
0080 // could be enough given that the p-columns (5 um radius) are in the corners of the cell
0081 // (no producing charge in there)
0082 // The output is vector of newly created charge in the neighbour pixel i+1 or i-1,
0083 // defined by its position higher than abs(half_pitch) and the the sign providing
0084 // the addition or subtraction in the pixel  (i+-1)
0085 std::vector<DigitizerUtility::EnergyDepositUnit> Pixel3DDigitizerAlgorithm::diffusion(
0086     const LocalPoint& pos,
0087     const float& ncarriers,
0088     const std::function<LocalVector(float, float)>& u_drift,
0089     const std::pair<float, float> hpitches,
0090     const float& thickness) const {
0091   // FIXME -- DM : Note that with a 0.3 will be enough (if using current sigma formulae)
0092   //          With the current sigma, this value is dependent of the thickness,
0093   //          Note that this formulae is coming from planar sensors, a similar
0094   //          study with data will be needed to extract the sigma for 3D
0095   const float max_migration_radius = 0.4_um;
0096   // Need to know which axis is the relevant one
0097   int displ_ind = -1;
0098   float pitch = 0.0;
0099 
0100   // Check the group is near the edge of the pixel, so diffusion will
0101   // be relevant in order to migrate between pixel cells
0102   if (hpitches.first - std::abs(pos.x()) < max_migration_radius) {
0103     displ_ind = 0;
0104     pitch = hpitches.first;
0105   } else if (hpitches.second - std::abs(pos.y()) < max_migration_radius) {
0106     displ_ind = 1;
0107     pitch = hpitches.second;
0108   } else {
0109     // Nothing to do, too far away
0110     return std::vector<DigitizerUtility::EnergyDepositUnit>();
0111   }
0112 
0113   // The new EnergyDeposits in the neighbour pixels
0114   // (defined by +1 to the right (first axis) and +1 to the up (second axis)) <-- XXX
0115   std::vector<DigitizerUtility::EnergyDepositUnit> migrated_charge;
0116 
0117   // FIXME -- DM
0118   const float diffusion_step = 0.1_um;
0119 
0120   // The position while drifting
0121   std::vector<float> pos_moving({pos.x(), pos.y(), pos.z()});
0122   // The drifting: drift field and steps
0123   std::function<std::vector<float>(int)> do_step =
0124       [&pos_moving, &u_drift, diffusion_step](int i) -> std::vector<float> {
0125     auto dd = u_drift(pos_moving[0], pos_moving[1]);
0126     return std::vector<float>({i * diffusion_step * dd.x(), i * diffusion_step * dd.y(), i * diffusion_step * dd.z()});
0127   };
0128 
0129   LogDebug("Pixel3DDigitizerAlgorithm::diffusion")
0130       << "\nMax. radius from the pixel edge to migrate charge: " << max_migration_radius * 1.0_um_inv << " [um]"
0131       << "\nMigration axis: " << displ_ind
0132       << "\n(super-)Charge distance to the pixel edge: " << (pitch - pos_moving[displ_ind]) * 1.0_um_inv << " [um]";
0133 
0134   // How many sigmas (probably a configurable, to be decided not now)
0135   const float N_SIGMA = 3.0;
0136 
0137   // Start the drift and check every step
0138   // Some variables needed
0139   float current_carriers = ncarriers;
0140   std::vector<float> newpos({pos_moving[0], pos_moving[1], pos_moving[2]});
0141   float distance_edge = 0.0_um;
0142   // Current diffusion value
0143   const float sigma = 0.4_um;
0144   for (int i = 1;; ++i) {
0145     std::transform(pos_moving.begin(), pos_moving.end(), do_step(i).begin(), pos_moving.begin(), std::plus<float>());
0146     distance_edge = pitch - std::abs(pos_moving[displ_ind]);
0147     // Get the amount of charge on the neighbor pixel: note the
0148     // transformation to a Normal
0149     float migrated_e = current_carriers * 0.5 * (1.0 - std::erf(distance_edge / (sigma * std::sqrt(2.0))));
0150 
0151     LogDebug("(super-)charge diffusion") << "step-" << i << ", Current carriers Ne= " << current_carriers << ","
0152                                          << "r=(" << pos_moving[0] * 1.0_um_inv << ", " << pos_moving[1] * 1.0_um_inv
0153                                          << ", " << pos_moving[2] * 1.0_um_inv << ") [um], "
0154                                          << "Migrated charge: " << migrated_e;
0155 
0156     // Move the migrated charge
0157     current_carriers -= migrated_e;
0158 
0159     // Either far away from the edge or almost half of the carriers already migrated
0160     if (std::abs(distance_edge) >= max_migration_radius || current_carriers <= 0.5 * ncarriers) {
0161       break;
0162     }
0163 
0164     // Create the ionization point:
0165     // First update the newpos vector: the new charge position at the neighbouring pixel
0166     // is created in the same position as its "parent carriers"
0167     // except the direction of migration
0168     std::vector<float> newpos(pos_moving);
0169     // Let's create the new charge carriers around 3 sigmas away
0170     newpos[displ_ind] += std::copysign(N_SIGMA * sigma, newpos[displ_ind]);
0171     migrated_charge.push_back(DigitizerUtility::EnergyDepositUnit(migrated_e, newpos[0], newpos[1], newpos[2]));
0172   }
0173   return migrated_charge;
0174 }
0175 
0176 // ======================================================================
0177 //
0178 // Drift the charge segments to the column (collection surface)
0179 // Include the effect of E-field and B-field
0180 //
0181 // =====================================================================
0182 std::vector<DigitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::drift(
0183     const PSimHit& hit,
0184     const Phase2TrackerGeomDetUnit* pixdet,
0185     const GlobalVector& bfield,
0186     const std::vector<DigitizerUtility::EnergyDepositUnit>& ionization_points) const {
0187   return drift(hit, pixdet, bfield, ionization_points, true);
0188 }
0189 std::vector<DigitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::drift(
0190     const PSimHit& hit,
0191     const Phase2TrackerGeomDetUnit* pixdet,
0192     const GlobalVector& bfield,
0193     const std::vector<DigitizerUtility::EnergyDepositUnit>& ionization_points,
0194     bool diffusion_activated) const {
0195   // -- Current reference system is placed in the center of the module
0196   // -- The natural reference frame should be discribed taking advantatge of
0197   // -- the cylindrical nature of the pixel geometry -->
0198   // -- the new reference frame should be placed in the center of the n-column, and in the
0199   // -- surface of the ROC using cylindrical coordinates
0200 
0201   // Get ROC pitch, half_pitch and sensor thickness to be used to create the
0202   // proxy pixel cell reference frame
0203   const auto pitch = pixdet->specificTopology().pitch();
0204   const auto half_pitch = std::make_pair<float, float>(pitch.first * 0.5, pitch.second * 0.5);
0205   const float thickness = pixdet->specificSurface().bounds().thickness();
0206   const int nrows = pixdet->specificTopology().nrows();
0207   const int ncolumns = pixdet->specificTopology().ncolumns();
0208   const float pix_rounding = 0.99;
0209 
0210   // The maximum radial distance is going to be used to evaluate radiation damage XXX?
0211   const float max_radial_distance =
0212       std::sqrt(half_pitch.first * half_pitch.first + half_pitch.second * half_pitch.second);
0213 
0214   // All pixels are going to be translated to a proxy pixel cell (all pixels should behave
0215   // equally no matter their position w.r.t. the module) and describe the movements there
0216   // Define the center of the pixel local reference frame: the current cartesian local reference
0217   // frame is centered at half_width_x,half_width_y,half_thickness
0218   // XXX -- This info could be obtained at init/construction time?
0219   LocalPoint center_proxy_cell(half_pitch.first, half_pitch.second, -0.5 * thickness);
0220 
0221   LogDebug("Pixel3DDigitizerAlgorithm::drift")
0222       << "Pixel pitch:" << pitch.first * 1.0_um_inv << ", " << pitch.second * 1.0_um_inv << " [um]";
0223 
0224   // And the drift direction (assumed same for all the sensor)
0225   // XXX call the function which will return a functional
0226   std::function<LocalVector(float, float)> drift_direction = [](float x, float y) -> LocalVector {
0227     const float theta = std::atan2(y, x);
0228     return LocalVector(-std::cos(theta), -std::sin(theta), 0.0);
0229   };
0230   // The output
0231   std::vector<DigitizerUtility::SignalPoint> collection_points;
0232   //collection_points.resize(ionization_points.size());
0233   collection_points.reserve(ionization_points.size());
0234 
0235   // Radiation damage limit of application
0236   // (XXX: No sense for 3D, let this be until decided what algorithm to use)
0237   const float RAD_DAMAGE = 0.001;
0238 
0239   for (const auto& super_charge : ionization_points) {
0240     // Extract the pixel cell
0241     auto current_pixel = pixdet->specificTopology().pixel(LocalPoint(super_charge.x(), super_charge.y()));
0242     // `pixel` function does not check to be in the ROC bounds,
0243     // so check it here and fix potential rounding problems.
0244     // Careful, this is assuming a rounding problem (1 unit), more than 1 pixel
0245     // away is probably showing some backward problem worth it to track.
0246     // This is also correcting out of bounds migrated charge from diffusion.
0247     // The charge will be moved to the edge of the row/column.
0248     current_pixel.first = std::clamp(current_pixel.first, float(0.0), (nrows - 1) + pix_rounding);
0249     current_pixel.second = std::clamp(current_pixel.second, float(0.0), (ncolumns - 1) + pix_rounding);
0250 
0251     const auto current_pixel_int = std::make_pair(std::floor(current_pixel.first), std::floor(current_pixel.second));
0252 
0253     // Convert to the 1x1 proxy pixel cell (pc), where all calculations are going to be
0254     // performed. The pixel is scaled to the actual pitch
0255     const auto relative_position_at_pc =
0256         std::make_pair((current_pixel.first - current_pixel_int.first) * pitch.first,
0257                        (current_pixel.second - current_pixel_int.second) * pitch.second);
0258 
0259     // Changing the reference frame to the proxy pixel cell
0260     LocalPoint position_at_pc(relative_position_at_pc.first - center_proxy_cell.x(),
0261                               relative_position_at_pc.second - center_proxy_cell.y(),
0262                               super_charge.z() - center_proxy_cell.z());
0263 
0264     LogDebug("Pixel3DDigitizerAlgorithm::drift")
0265         << "(super-)Charge\nlocal position: (" << super_charge.x() * 1.0_um_inv << ", " << super_charge.y() * 1.0_um_inv
0266         << ", " << super_charge.z() * 1.0_um_inv << ") [um]"
0267         << "\nMeasurement Point (row,column) (" << current_pixel.first << ", " << current_pixel.second << ")"
0268         << "\nProxy pixel-cell frame (centered at  left-back corner): (" << relative_position_at_pc.first * 1.0_um_inv
0269         << ", " << relative_position_at_pc.second * 1.0_um_inv << ") [um]"
0270         << "\nProxy pixel-cell frame (centered at n-column): (" << position_at_pc.x() * 1.0_um_inv << ", "
0271         << position_at_pc.y() * 1.0_um_inv << ") [um] "
0272         << "\nNe=" << super_charge.energy() << " electrons";
0273 
0274     // Check if the point is inside any of the column --> no charge was actually created then
0275     if (is_inside_n_column_(position_at_pc, thickness) || is_inside_ohmic_column_(position_at_pc, half_pitch)) {
0276       LogDebug("Pixel3DDigitizerAlgorithm::drift") << "Remove charge,  inside the n-column or p-column!!";
0277       continue;
0278     }
0279 
0280     float nelectrons = super_charge.energy();
0281     // XXX -- Diffusion: using the center frame
0282     if (diffusion_activated) {
0283       auto migrated_charges = diffusion(position_at_pc, super_charge.energy(), drift_direction, half_pitch, thickness);
0284       for (auto& mc : migrated_charges) {
0285         // Remove the migrated charges
0286         nelectrons -= mc.energy();
0287         // and convert back to the pixel ref. system
0288         // Low-left origin/pitch -> relative within the pixel (a)
0289         // Adding the pixel
0290         const float pixel_x = current_pixel_int.first + (mc.x() + center_proxy_cell.x()) / pitch.first;
0291         const float pixel_y = current_pixel_int.second + (mc.y() + center_proxy_cell.y()) / pitch.second;
0292         const auto lp = pixdet->specificTopology().localPosition(MeasurementPoint(pixel_x, pixel_y));
0293         // Remember: the drift function will move the reference system to the top. We need to subtract
0294         // (center_proxy_cell.z() is a constant negative value) what we previously added in order to
0295         // avoid a double translation when calling the drift function below the drift function
0296         // initially considers the reference system centered in the module at half thickness)
0297         mc.migrate_position(LocalPoint(lp.x(), lp.y(), mc.z() + center_proxy_cell.z()));
0298       }
0299       if (!migrated_charges.empty()) {
0300         LogDebug("Pixel3DDigitizerAlgorithm::drift") << "****************"
0301                                                      << "MIGRATING (super-)charges"
0302                                                      << "****************";
0303         // Drift this charges on the other pixel
0304         auto mig_colpoints = drift(hit, pixdet, bfield, migrated_charges, false);
0305         collection_points.insert(std::end(collection_points), mig_colpoints.begin(), mig_colpoints.end());
0306         LogDebug("Pixel3DDigitizerAlgorithm::drift") << "*****************"
0307                                                      << "DOME MIGRATION"
0308                                                      << "****************";
0309       }
0310     }
0311 
0312     // Perform the drift, and check a potential lost of carriers because
0313     // they reach the pasivation region (-z < thickness/2)
0314     // XXX: not doing nothing, the carriers reach the electrode surface,
0315     const float drift_distance = position_at_pc.perp() - np_column_radius_;
0316 
0317     // Insert a charge loss due to Rad Damage here
0318     // XXX: ??
0319     float energyOnCollector = nelectrons;
0320     // FIXME: is this correct?  Not for 3D...
0321 
0322     if (pseudoRadDamage_ >= RAD_DAMAGE) {
0323       const float module_radius = pixdet->surface().position().perp();
0324       if (module_radius <= pseudoRadDamageRadius_) {
0325         const float kValue = pseudoRadDamage_ / (module_radius * module_radius);
0326         energyOnCollector = energyOnCollector * std::exp(-1.0 * kValue * drift_distance / max_radial_distance);
0327       }
0328     }
0329     LogDebug("Pixel3DDigitizerAlgorithm::drift")
0330         << "Drift distance = " << drift_distance * 1.0_um_inv << " [um], "
0331         << "Initial electrons = " << super_charge.energy()
0332         << " [electrons], Electrons after loss/diff= " << energyOnCollector << " [electrons] ";
0333     // Load the Charge distribution parameters
0334     // XXX -- probably makes no sense the SignalPoint anymore...
0335     collection_points.push_back(DigitizerUtility::SignalPoint(
0336         current_pixel_int.first, current_pixel_int.second, 0.0, 0.0, hit.tof(), energyOnCollector));
0337   }
0338 
0339   return collection_points;
0340 }
0341 
0342 // ====================================================================
0343 // Signal is already "induced" (actually electrons transported to the
0344 // n-column) at the electrode. Just collecting and adding-up all pixel
0345 // signal and linking it to the simulated energy deposit (hit)
0346 void Pixel3DDigitizerAlgorithm::induce_signal(const PSimHit& hit,
0347                                               const size_t hitIndex,
0348                                               const uint32_t tofBin,
0349                                               const Phase2TrackerGeomDetUnit* pixdet,
0350                                               const std::vector<DigitizerUtility::SignalPoint>& collection_points) {
0351   // X  - Rows, Left-Right
0352   // Y  - Columns, Down-Up
0353   const uint32_t detId = pixdet->geographicalId().rawId();
0354   // Accumulated signal at each channel of this detector
0355   signal_map_type& the_signal = _signal[detId];
0356 
0357   // Choose the proper pixel-to-channel converter
0358   std::function<int(int, int)> pixelToChannel =
0359       pixelFlag_ ? PixelDigi::pixelToChannel
0360                  : static_cast<std::function<int(int, int)> >(Phase2TrackerDigi::pixelToChannel);
0361 
0362   // Iterate over collection points on the collection plane
0363   for (const auto& pt : collection_points) {
0364     // Extract corresponding channel (position is already given in pixel indices)
0365     const int channel = pixelToChannel(pt.position().x(), pt.position().y());
0366 
0367     float corr_time = hit.tof() - pixdet->surface().toGlobal(hit.localPosition()).mag() * c_inv;
0368     if (makeDigiSimLinks_) {
0369       the_signal[channel] +=
0370           DigitizerUtility::Amplitude(pt.amplitude(), &hit, pt.amplitude(), corr_time, hitIndex, tofBin);
0371     } else {
0372       the_signal[channel] += DigitizerUtility::Amplitude(pt.amplitude(), nullptr, pt.amplitude());
0373     }
0374 
0375     LogDebug("Pixel3DDigitizerAlgorithm::induce_signal")
0376         << " Induce charge at row,col:" << pt.position() << " N_electrons:" << pt.amplitude() << " [Channel:" << channel
0377         << "]\n   [Accumulated signal in this channel:" << the_signal[channel].ampl() << "] "
0378         << " Global index linked PSimHit:" << hitIndex;
0379   }
0380 }