diff --git a/SimTracker/SiPhase2Digitizer/plugins/Pixel3DDigitizerAlgorithm.cc b/SimTracker/SiPhase2Digitizer/plugins/Pixel3DDigitizerAlgorithm.cc index 3d80a683a7cff..421425ebe40f7 100644 --- a/SimTracker/SiPhase2Digitizer/plugins/Pixel3DDigitizerAlgorithm.cc +++ b/SimTracker/SiPhase2Digitizer/plugins/Pixel3DDigitizerAlgorithm.cc @@ -158,56 +158,45 @@ std::vector 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 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()); 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 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 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 Pixel3DDigitizerAlgorithm::drift( const GlobalVector& bfield, const std::vector& 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 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 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,8 +308,8 @@ std::vector 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) @@ -328,8 +317,10 @@ std::vector Pixel3DDigitizerAlgorithm::drift( 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; - ; } } diff --git a/SimTracker/SiPhase2Digitizer/python/PixelTestBeamValidation_cfi.py b/SimTracker/SiPhase2Digitizer/python/PixelTestBeamValidation_cfi.py index ff7accf0ef1f7..400bf88f511a4 100644 --- a/SimTracker/SiPhase2Digitizer/python/PixelTestBeamValidation_cfi.py +++ b/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(), diff --git a/SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.cc b/SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.cc index d7145e4c8ef7e..484edcce93521 100644 --- a/SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.cc +++ b/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 @@ -40,6 +41,7 @@ PixelTestBeamValidation::PixelTestBeamValidation(const edm::ParameterSet& iConfi : config_(iConfig), geomType_(iConfig.getParameter("GeometryType")), //phiValues(iConfig.getParameter >("PhiAngles")), + thresholdInElectrons_(iConfig.getParameter("ThresholdInElectrons")), electronsPerADC_(iConfig.getParameter("ElectronsPerADC")), tracksEntryAngleX_( iConfig.getUntrackedParameter>("TracksEntryAngleX", std::vector())), @@ -51,6 +53,9 @@ PixelTestBeamValidation::PixelTestBeamValidation(const edm::ParameterSet& iConfi simTrackToken_(consumes(iConfig.getParameter("SimTrackSource"))) { LogDebug("PixelTestBeamValidation") << ">>> Construct PixelTestBeamValidation "; + // The value to be used for ToT == 0 in electrons + electronsAtToT0_ = 0.5 * (thresholdInElectrons_ + electronsPerADC_); + const std::vector psimhit_v(config_.getParameter>("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> 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()); diff --git a/SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.h b/SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.h index 25720d0020e37..9e85c5ad7f2c3 100644 --- a/SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.h +++ b/SimTracker/SiPhase2Digitizer/test/PixelTestBeamValidation.h @@ -145,6 +145,7 @@ class PixelTestBeamValidation : public DQMEDAnalyzer { std::map vME_track_dydzAngle_; std::map vME_dx1D_; std::map vME_dy1D_; + std::map vME_dxy2D_; std::map vME_digi_charge1D_; std::map vME_digi_chargeElec1D_; std::map 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 tracksEntryAngleX_; std::vector tracksEntryAngleY_;