Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:05:00

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