Merge pull request #33300 from claralasa/charge_migration_fixes

[Upgrade Phase-2 Pixel 3D Digitizer] Fix misevaluation of charge when migrating due to diffusion

SimTracker/SiPhase2Digitizer/plugins/Pixel3DDigitizerAlgorithm.cc

@@ -158,56 +158,45 @@ std::vector<DigitizerUtility::EnergyDepositUnit> Pixel3DDigitizerAlgorithm::diff
       << "\nMigration axis: " << displ_ind
       << "\n(super-)Charge distance to the pixel edge: " << (pitch - pos_moving[displ_ind]) * 1.0_um_inv << " [um]";
 
-  // FIXME -- Sigma reference, DM?
-  const float distance0 = 300.0_um;
-  const float sigma0 = 3.4_um;
-  // FIXME -- Tolerance, DM?
-  const float TOL = 1e-6;
   // How many sigmas (probably a configurable, to be decided not now)
   const float N_SIGMA = 3.0;
 
   // Start the drift and check every step
-  // initial position
-  int i = 0;
   // Some variables needed
   float current_carriers = ncarriers;
   std::vector<float> newpos({pos_moving[0], pos_moving[1], pos_moving[2]});
   float distance_edge = 0.0_um;
-  do {
+  // Current diffusion value
+  const float sigma = 0.4_um;
+  for (int i = 1;; ++i) {
     std::transform(pos_moving.begin(), pos_moving.end(), do_step(i).begin(), pos_moving.begin(), std::plus<float>());
     distance_edge = pitch - std::abs(pos_moving[displ_ind]);
-    // current diffusion value
-    double sigma = std::sqrt(i * diffusion_step / distance0) * (distance0 / thickness) * sigma0;
     // Get the amount of charge on the neighbor pixel: note the
     // transformation to a Normal
-    float migrated_e = current_carriers * 0.5 * (1.0 - std::erf(distance_edge / sigma));
+    float migrated_e = current_carriers * 0.5 * (1.0 - std::erf(distance_edge / (sigma * std::sqrt(2.0))));
 
-    LogDebug("(super-)charge diffusion") << "step-" << i << ", Initial Ne= " << ncarriers << ", "
+    LogDebug("(super-)charge diffusion") << "step-" << i << ", Current carriers Ne= " << current_carriers << ","
                                          << "r=(" << pos_moving[0] * 1.0_um_inv << ", " << pos_moving[1] * 1.0_um_inv
                                          << ", " << pos_moving[2] * 1.0_um_inv << ") [um], "
                                          << "Migrated charge: " << migrated_e;
 
-    // No charge was migrated (ignore creation time)
-    if (i != 0) {
-      // At least 1 electron migrated
-      if ((migrated_e - TOL) < 1.0) {
-        break;
-      }
-      // Move the migrated charge
-      current_carriers -= migrated_e;
-      // Create the ionization point:
-      // First update the newpos vector: the new charge positions at the neighbourg pixels
-      // are created in the same position that its "parent carriers"
-      // except the direction of migration
-      std::vector<float> newpos(pos_moving);
-      // Lest create the new charges around 3 sigmas away
-      newpos[displ_ind] += std::copysign(N_SIGMA * sigma, newpos[displ_ind]);
-      migrated_charge.push_back(DigitizerUtility::EnergyDepositUnit(migrated_e, newpos[0], newpos[1], newpos[2]));
+    // Move the migrated charge
+    current_carriers -= migrated_e;
+
+    // Either far away from the edge or almost half of the carriers already migrated
+    if (std::abs(distance_edge) >= max_migration_radius || current_carriers <= 0.5 * ncarriers) {
+      break;
     }
-    // Next step
-    ++i;
-  } while (std::abs(distance_edge) < max_migration_radius && current_carriers > 0.5 * ncarriers);
 
+    // Create the ionization point:
+    // First update the newpos vector: the new charge position at the neighbouring pixel
+    // is created in the same position as its "parent carriers"
+    // except the direction of migration
+    std::vector<float> newpos(pos_moving);
+    // Let's create the new charge carriers around 3 sigmas away
+    newpos[displ_ind] += std::copysign(N_SIGMA * sigma, newpos[displ_ind]);
+    migrated_charge.push_back(DigitizerUtility::EnergyDepositUnit(migrated_e, newpos[0], newpos[1], newpos[2]));
+  }
   return migrated_charge;
 }
 
@@ -230,10 +219,10 @@ std::vector<DigitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::drift(
     const GlobalVector& bfield,
     const std::vector<DigitizerUtility::EnergyDepositUnit>& ionization_points,
     bool diffusion_activated) const {
-  // -- Current reference system is placed in the center on the module
+  // -- Current reference system is placed in the center of the module
   // -- The natural reference frame should be discribed taking advantatge of
   // -- the cylindrical nature of the pixel geometry -->
-  // -- the new reference frame should be place in the center of the columncy, and in the
+  // -- the new reference frame should be placed in the center of the n-column, and in the
   // -- surface of the ROC using cylindrical coordinates
 
   // Get ROC pitch, half_pitch and sensor thickness to be used to create the
@@ -245,7 +234,7 @@ std::vector<DigitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::drift(
   const int ncolumns = pixdet->specificTopology().ncolumns();
   const float pix_rounding = 0.99;
 
-  // the maximum radial distance is going to be use to evaluate radiation damage XXX?
+  // The maximum radial distance is going to be used to evaluate radiation damage XXX?
   const float max_radial_distance =
       std::sqrt(half_pitch.first * half_pitch.first + half_pitch.second * half_pitch.second);
 
@@ -303,7 +292,7 @@ std::vector<DigitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::drift(
         << "(super-)Charge\nlocal position: (" << super_charge.x() * 1.0_um_inv << ", " << super_charge.y() * 1.0_um_inv
         << ", " << super_charge.z() * 1.0_um_inv << ") [um]"
         << "\nMeasurement Point (row,column) (" << current_pixel.first << ", " << current_pixel.second << ")"
-        << "\nProxy pixel-cell frame (centered at  left-down corner): (" << relative_position_at_pc.first * 1.0_um_inv
+        << "\nProxy pixel-cell frame (centered at  left-back corner): (" << relative_position_at_pc.first * 1.0_um_inv
         << ", " << relative_position_at_pc.second * 1.0_um_inv << ") [um]"
         << "\nProxy pixel-cell frame (centered at n-column): (" << position_at_pc.x() * 1.0_um_inv << ", "
         << position_at_pc.y() * 1.0_um_inv << ") [um] "
@@ -319,17 +308,19 @@ std::vector<DigitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::drift(
     // XXX -- Diffusion: using the center frame
     if (diffusion_activated) {
       auto migrated_charges = diffusion(position_at_pc, super_charge.energy(), drift_direction, half_pitch, thickness);
-      // remove the migrated charges
       for (auto& mc : migrated_charges) {
+        // Remove the migrated charges
         nelectrons -= mc.energy();
         // and convert back to the pixel ref. system
         // Low-left origin/pitch -> relative within the pixel (a)
         // Adding the pixel
         const float pixel_x = current_pixel_int.first + (mc.x() + center_proxy_cell.x()) / pitch.first;
         const float pixel_y = current_pixel_int.second + (mc.y() + center_proxy_cell.y()) / pitch.second;
         const auto lp = pixdet->specificTopology().localPosition(MeasurementPoint(pixel_x, pixel_y));
-        //Remember: the drift function will move the reference system to the bottom. We need to add what we previously subtract
-        //in order to avoid a double translation when calling the drift function once again below
+        // Remember: the drift function will move the reference system to the top. We need to subtract
+        // (center_proxy_cell.z() is a constant negative value) what we previously added in order to
+        // avoid a double translation when calling the drift function below the drift function
+        // initially considers the reference system centered in the module at half thickness)
         mc.migrate_position(LocalPoint(lp.x(), lp.y(), mc.z() + center_proxy_cell.z()));
       }
       if (!migrated_charges.empty()) {
@@ -412,6 +403,5 @@ void Pixel3DDigitizerAlgorithm::induce_signal(const PSimHit& hit,
         << " Induce charge at row,col:" << pt.position() << " N_electrons:" << pt.amplitude() << " [Channel:" << channel
         << "]\n   [Accumulated signal in this channel:" << the_signal[channel].ampl() << "] "
         << " Global index linked PSimHit:" << hitIndex;
-    ;
   }
 }

SimTracker/SiPhase2Digitizer/python/PixelTestBeamValidation_cfi.py

@@ -48,6 +48,14 @@
     xmin = cms.double(-150),
     xmax = cms.double(150)
     )
+Dxy2D = cms.PSet(
+    Nxbins = cms.int32(Nx),
+    Nybins = cms.int32(Ny),
+    xmin   = cms.double(-50),
+    xmax   = cms.double(50),
+    ymin   = cms.double(-150),
+    ymax   = cms.double(150)
+    )
 SimClusterCharge = cms.PSet(
     Nxbins = cms.int32(300),
     Nybins = cms.int32(201),
@@ -109,6 +117,7 @@
 
 dqmcell = DQMEDAnalyzer('PixelTestBeamValidation',
     # WARNING: be sure it is the same value used with the Pixel3DDigitizer
+    ThresholdInElectrons = cms.double(1000.0),
     ElectronsPerADC   = cms.double(1600.0),
     TracksEntryAngleX = cms.untracked.vdouble(-np.radians(2.0),np.radians(2.0)),
     TracksEntryAngleY = cms.untracked.vdouble(-np.radians(2.0),np.radians(2.0)),
@@ -143,6 +152,7 @@
     DigiRZ = TrackRZ.clone(),
     Dx1D = Dx1D.clone(),
     Dy1D = Dy1D.clone(),
+    Dxy2D = Dxy2D.clone(),
     DigiCharge1D = DigiCharge1D.clone(),
     SimClusterCharge = SimClusterCharge.clone(),
 

SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.cc

@@ -20,6 +20,7 @@
 #include "DataFormats/Math/interface/CMSUnits.h"
 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
 #include "DataFormats/Common/interface/Handle.h"
+#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
 
 // system
 #include <algorithm>
@@ -40,6 +41,7 @@ PixelTestBeamValidation::PixelTestBeamValidation(const edm::ParameterSet& iConfi
     : config_(iConfig),
       geomType_(iConfig.getParameter<std::string>("GeometryType")),
       //phiValues(iConfig.getParameter<std::vector<double> >("PhiAngles")),
+      thresholdInElectrons_(iConfig.getParameter<double>("ThresholdInElectrons")),
       electronsPerADC_(iConfig.getParameter<double>("ElectronsPerADC")),
       tracksEntryAngleX_(
           iConfig.getUntrackedParameter<std::vector<double>>("TracksEntryAngleX", std::vector<double>())),
@@ -51,6 +53,9 @@ PixelTestBeamValidation::PixelTestBeamValidation(const edm::ParameterSet& iConfi
       simTrackToken_(consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("SimTrackSource"))) {
   LogDebug("PixelTestBeamValidation") << ">>> Construct PixelTestBeamValidation ";
 
+  // The value to be used for ToT == 0 in electrons
+  electronsAtToT0_ = 0.5 * (thresholdInElectrons_ + electronsPerADC_);
+
   const std::vector<edm::InputTag> psimhit_v(config_.getParameter<std::vector<edm::InputTag>>("PSimHitSource"));
 
   for (const auto& itag : psimhit_v) {
@@ -265,11 +270,14 @@ void PixelTestBeamValidation::analyze(const edm::Event& iEvent, const edm::Event
         vME_digi_RZMap_->Fill(digi_global_pos.z(), std::hypot(digi_global_pos.x(), digi_global_pos.y()));
         // Create the MC-cluster
         cluster_tot += current_digi.adc();
-        // Add 0.5 to allow ToT = 0 (valid value)
-        cluster_tot_elec += (current_digi.adc() + 0.5) * electronsPerADC_;
+        // Assign the middle value between the threshold and the ElectronPerADC parameter
+        // to the first bin in order to allow ToT = 0 (valid value)
+        const double pixel_charge_elec =
+            (bool(current_digi.adc()) ? (current_digi.adc() * electronsPerADC_) : electronsAtToT0_);
+        cluster_tot_elec += pixel_charge_elec;
         // Use the center of the pixel
-        cluster_position.first += current_digi.adc() * (current_digi.row() + 0.5);
-        cluster_position.second += current_digi.adc() * (current_digi.column() + 0.5);
+        cluster_position.first += pixel_charge_elec * (current_digi.row() + 0.5);
+        cluster_position.second += pixel_charge_elec * (current_digi.column() + 0.5);
         // Size
         cluster_size_xy.first.insert(current_digi.row());
         cluster_size_xy.second.insert(current_digi.column());
@@ -282,8 +290,8 @@ void PixelTestBeamValidation::analyze(const edm::Event& iEvent, const edm::Event
       vME_clsize1Dy_[me_unit]->Fill(cluster_size_xy.second.size());
 
       // mean weighted
-      cluster_position.first /= double(cluster_tot);
-      cluster_position.second /= double(cluster_tot);
+      cluster_position.first /= cluster_tot_elec;
+      cluster_position.second /= cluster_tot_elec;
 
       // -- XXX Be careful, secondaries with already used the digis
       //        are going the be lost (then lost on efficiency)
@@ -302,6 +310,7 @@ void PixelTestBeamValidation::analyze(const edm::Event& iEvent, const edm::Event
         vME_charge_elec1D_[me_unit]->Fill(cluster_tot_elec);
         vME_dx1D_[me_unit]->Fill(dx_um);
         vME_dy1D_[me_unit]->Fill(dy_um);
+        vME_dxy2D_[me_unit]->Fill(dx_um, dy_um);
         // The track energy loss corresponding to that cluster
         vME_sim_cluster_charge_[me_unit]->Fill(psh->energyLoss() * 1.0_inv_keV, cluster_tot_elec);
       }
@@ -430,7 +439,11 @@ void PixelTestBeamValidation::bookHistograms(DQMStore::IBooker& ibooker,
       vME_dx1D_[me_unit] = setupH1D_(
           ibooker, "Dx1D", "MC-truth DIGI cluster residuals X;x_{PSimHit}-x^{cluster}_{digi} [#mum];N_{digi clusters}");
       vME_dy1D_[me_unit] = setupH1D_(
-          ibooker, "Dy1D", "MC-truth DIGI cluster residual Ys;y_{PSimHit}-y^{cluster}_{digi} [#mum];N_{digi clusters}");
+          ibooker, "Dy1D", "MC-truth DIGI cluster residual Y;y_{PSimHit}-y^{cluster}_{digi} [#mum];N_{digi clusters}");
+      vME_dxy2D_[me_unit] = setupH2D_(ibooker,
+                                      "Dxy2D",
+                                      "MC-truth DIGI cluster residuals;x_{PSimHit}-x^{cluster}_{digi} "
+                                      "[#mum];y_{PSimHit}-y^{cluster}_{digi} [#mum];N_{digi clusters}");
       vME_digi_charge1D_[me_unit] = setupH1D_(ibooker, "DigiCharge1D", "Digi charge;digi charge [ToT];N_{digi}");
       vME_sim_cluster_charge_[me_unit] =
           setupH2D_(ibooker,
@@ -712,8 +725,12 @@ std::set<std::pair<int, int>> PixelTestBeamValidation::get_illuminated_pixels_(c
   const double max_x = std::max(ps.entryPoint().x(), ps.exitPoint().x());
   const double max_y = std::max(ps.entryPoint().y(), ps.exitPoint().y());
   // Get the position in readout units for each point
-  const auto min_pos = tkDetUnit->specificTopology().measurementPosition(LocalPoint(min_x, min_y));
-  const auto max_pos = tkDetUnit->specificTopology().measurementPosition(LocalPoint(max_x, max_y));
+  const auto min_pos_pre = tkDetUnit->specificTopology().measurementPosition(LocalPoint(min_x, min_y));
+  const auto max_pos_pre = tkDetUnit->specificTopology().measurementPosition(LocalPoint(max_x, max_y));
+  // Adding a unit at each side in order to include charge migration
+  const MeasurementPoint min_pos(std::max(0.0, min_pos_pre.x() - 1.0), std::max(0.0, min_pos_pre.y() - 1.0));
+  const MeasurementPoint max_pos(std::min(max_pos_pre.x() + 1.0, tkDetUnit->specificTopology().nrows() - 1.0),
+                                 std::min(max_pos_pre.y() + 1.0, tkDetUnit->specificTopology().ncolumns() - 1.0));
   // Count how many cells has passed. Use the most conservative rounding:
   // round for maximums and int (floor) for minimums
   //const int ncells_x = std::round(max_pos.x())-std::floor(min_pos.x());

SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.h

@@ -145,6 +145,7 @@ class PixelTestBeamValidation : public DQMEDAnalyzer {
   std::map<int, MonitorElement *> vME_track_dydzAngle_;
   std::map<int, MonitorElement *> vME_dx1D_;
   std::map<int, MonitorElement *> vME_dy1D_;
+  std::map<int, MonitorElement *> vME_dxy2D_;
   std::map<int, MonitorElement *> vME_digi_charge1D_;
   std::map<int, MonitorElement *> vME_digi_chargeElec1D_;
   std::map<int, MonitorElement *> vME_sim_cluster_charge_;
@@ -167,9 +168,11 @@ class PixelTestBeamValidation : public DQMEDAnalyzer {
   // Geometry to use
   std::string geomType_;
 
-  // The conversion between ToT to electrons (Be carefull, this should
+  // The threshold and the conversion between ToT to electrons (Be carefull, this should
   // be using the same value used in the digitization module)
+  double thresholdInElectrons_;
   double electronsPerADC_;
+  double electronsAtToT0_;
   // The tracks entry angle to accept (if any)
   std::vector<double> tracksEntryAngleX_;
   std::vector<double> tracksEntryAngleY_;