Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Include charge reweighting for Phase-2 IT planar sensors #40417

Merged
merged 2 commits into from Jan 16, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
264 changes: 142 additions & 122 deletions SimTracker/Common/interface/SiPixelChargeReweightingAlgorithm.h

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions SimTracker/SiPhase2Digitizer/plugins/PSPDigitizerAlgorithm.h
@@ -1,5 +1,5 @@
#ifndef _SimTracker_SiPhase2Digitizer_PSPDigitizerAlgorithm_h
#define _SimTracker_SiPhase2Digitizer_PSPDigitizerAlgorithm_h
#ifndef SimTracker_SiPhase2Digitizer_PSPDigitizerAlgorithm_h
#define SimTracker_SiPhase2Digitizer_PSPDigitizerAlgorithm_h

#include "CondFormats/DataRecord/interface/SiPhase2OuterTrackerLorentzAngleRcd.h"
#include "FWCore/Utilities/interface/ESGetToken.h"
Expand Down
4 changes: 2 additions & 2 deletions SimTracker/SiPhase2Digitizer/plugins/PSSDigitizerAlgorithm.h
@@ -1,5 +1,5 @@
#ifndef _SimTracker_SiPhase2Digitizer_PSSDigitizerAlgorithm_h
#define _SimTracker_SiPhase2Digitizer_PSSDigitizerAlgorithm_h
#ifndef SimTracker_SiPhase2Digitizer_PSSDigitizerAlgorithm_h
#define SimTracker_SiPhase2Digitizer_PSSDigitizerAlgorithm_h

#include "CondFormats/DataRecord/interface/SiPhase2OuterTrackerLorentzAngleRcd.h"
#include "FWCore/Utilities/interface/ESGetToken.h"
Expand Down
@@ -1,35 +1,28 @@
#include <typeinfo>
#include <iostream>
#include <cmath>

#include "SimTracker/SiPhase2Digitizer/plugins/Phase2TrackerDigitizerAlgorithm.h"

#include "SimGeneral/NoiseGenerators/interface/GaussianTailNoiseGenerator.h"
#include "SimTracker/Common/interface/SiG4UniversalFluctuation.h"
#include <iostream>
#include <typeinfo>

#include "CLHEP/Random/RandGaussQ.h"

#include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
#include "CondFormats/SiPhase2TrackerObjects/interface/SiPhase2OuterTrackerLorentzAngle.h"
#include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
#include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
#include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
#include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
#include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
#include "CondFormats/SiPhase2TrackerObjects/interface/SiPhase2OuterTrackerLorentzAngle.h"

// Geometry
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/CommonTopologies/interface/PixelTopology.h"
#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
#include "Geometry/CommonTopologies/interface/PixelTopology.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "SimGeneral/NoiseGenerators/interface/GaussianTailNoiseGenerator.h"
#include "SimTracker/Common/interface/SiG4UniversalFluctuation.h"
#include "SimTracker/SiPhase2Digitizer/plugins/Phase2TrackerDigitizerAlgorithm.h"

using namespace edm;

Expand All @@ -41,6 +34,7 @@ namespace ph2tkdigialgo {
constexpr double m_muon = 105.658;
constexpr double m_proton = 938.272;
} // namespace ph2tkdigialgo

Phase2TrackerDigitizerAlgorithm::Phase2TrackerDigitizerAlgorithm(const edm::ParameterSet& conf_common,
const edm::ParameterSet& conf_specific,
edm::ConsumesCollector iC)
Expand Down Expand Up @@ -140,6 +134,10 @@ Phase2TrackerDigitizerAlgorithm::Phase2TrackerDigitizerAlgorithm(const edm::Para
pseudoRadDamageRadius_(conf_specific.exists("PseudoRadDamageRadius")
? conf_specific.getParameter<double>("PseudoRadDamageRadius")
: double(0.0)),
// Add pixel radiation damage
useChargeReweighting_(conf_specific.getParameter<bool>("UseReweighting")),
theSiPixelChargeReweightingAlgorithm_(
useChargeReweighting_ ? std::make_unique<SiPixelChargeReweightingAlgorithm>(conf_specific, iC) : nullptr),

// delta cutoff in MeV, has to be same as in OSCAR(0.030/cmsim=1.0 MeV
// tMax(0.030), // In MeV.
Expand Down Expand Up @@ -171,6 +169,7 @@ Phase2TrackerDigitizerAlgorithm::SubdetEfficiencies::SubdetEfficiencies(const ed
barrel_efficiencies = conf.getParameter<std::vector<double> >("EfficiencyFactors_Barrel");
endcap_efficiencies = conf.getParameter<std::vector<double> >("EfficiencyFactors_Endcap");
}

void Phase2TrackerDigitizerAlgorithm::accumulateSimHits(std::vector<PSimHit>::const_iterator inputBegin,
std::vector<PSimHit>::const_iterator inputEnd,
const size_t inputBeginGlobalIndex,
Expand Down Expand Up @@ -202,7 +201,8 @@ void Phase2TrackerDigitizerAlgorithm::accumulateSimHits(std::vector<PSimHit>::co

// compute induced signal on readout elements and add to _signal
// hit needed only for SimHit<-->Digi link
induce_signal(hit, simHitGlobalIndex, tofBin, pixdet, collection_points);
induce_signal(inputBegin, hit, simHitGlobalIndex, inputBeginGlobalIndex, tofBin, pixdet, collection_points);
LogDebug("Phase2DigitizerAlgorithm") << "Induced signal was computed";
}
++simHitGlobalIndex;
}
Expand Down Expand Up @@ -252,7 +252,8 @@ std::vector<digitizerUtility::EnergyDepositUnit> Phase2TrackerDigitizerAlgorithm
digitizerUtility::EnergyDepositUnit edu(energy, point); // define position,energy point
ionization_points.push_back(edu); // save
LogDebug("Phase2TrackerDigitizerAlgorithm")
<< i << " " << edu.x() << " " << edu.y() << " " << edu.z() << " " << edu.energy();
<< "For index = " << i << " EnergyDepositUnit-x = " << edu.x() << " EnergyDepositUnit-y = " << edu.y()
<< " EnergyDepositUnit-z = " << edu.z() << " EnergyDepositUnit-energy = " << edu.energy();
}
return ionization_points;
}
Expand Down Expand Up @@ -324,11 +325,11 @@ std::vector<digitizerUtility::SignalPoint> Phase2TrackerDigitizerAlgorithm::drif
const Phase2TrackerGeomDetUnit* pixdet,
const GlobalVector& bfield,
const std::vector<digitizerUtility::EnergyDepositUnit>& ionization_points) const {
LogDebug("Phase2TrackerDigitizerAlgorithm") << "enter drift ";
LogDebug("Phase2TrackerDigitizerAlgorithm") << "Enter drift ";

std::vector<digitizerUtility::SignalPoint> collection_points;
collection_points.reserve(ionization_points.size()); // set size
LocalVector driftDir = DriftDirection(pixdet, bfield, hit.detUnitId()); // get the charge drift direction
LocalVector driftDir = driftDirection(pixdet, bfield, hit.detUnitId()); // get the charge drift direction
if (driftDir.z() == 0.) {
LogWarning("Phase2TrackerDigitizerAlgorithm") << " pxlx: drift in z is zero ";
return collection_points;
Expand All @@ -348,8 +349,9 @@ std::vector<digitizerUtility::SignalPoint> Phase2TrackerDigitizerAlgorithm::drif
float stripPitch = pixdet->specificTopology().pitch().first;

LogDebug("Phase2TrackerDigitizerAlgorithm")
<< " Lorentz Tan " << TanLorenzAngleX << " " << TanLorenzAngleY << " " << CosLorenzAngleX << " "
<< CosLorenzAngleY << " " << moduleThickness * TanLorenzAngleX << " " << driftDir;
<< " Lorentz Tan-X " << TanLorenzAngleX << " Lorentz Tan-Y " << TanLorenzAngleY << " Lorentz Cos-X "
<< CosLorenzAngleX << " Lorentz Cos-Y " << CosLorenzAngleY
<< " ticknes * Lorentz Tan-X = " << moduleThickness * TanLorenzAngleX << " drift direction " << driftDir;

for (auto const& val : ionization_points) {
// position
Expand Down Expand Up @@ -418,8 +420,10 @@ std::vector<digitizerUtility::SignalPoint> Phase2TrackerDigitizerAlgorithm::drif
//
// Induce the signal on the collection plane of the active sensor area.
void Phase2TrackerDigitizerAlgorithm::induce_signal(
std::vector<PSimHit>::const_iterator inputBegin,
const PSimHit& hit,
const size_t hitIndex,
const size_t firstHitIndex,
const uint32_t tofBin,
const Phase2TrackerGeomDetUnit* pixdet,
const std::vector<digitizerUtility::SignalPoint>& collection_points) {
Expand All @@ -430,7 +434,7 @@ void Phase2TrackerDigitizerAlgorithm::induce_signal(
signal_map_type& theSignal = _signal[detID];

LogDebug("Phase2TrackerDigitizerAlgorithm")
<< " enter induce_signal, " << topol->pitch().first << " " << topol->pitch().second;
<< "Enter induce_signal, Pitch size is " << topol->pitch().first << " cm vs " << topol->pitch().second << " cm";

// local map to store pixels hit by 1 Hit.
using hit_map_type = std::map<int, float, std::less<int> >;
Expand Down Expand Up @@ -563,14 +567,51 @@ void Phase2TrackerDigitizerAlgorithm::induce_signal(
}
}
// Fill the global map with all hit pixels from this event
bool reweighted = false;
size_t referenceIndex4CR = 0;

if (useChargeReweighting_) {
if (hit.processType() == 0) {
referenceIndex4CR = hitIndex;
reweighted =
theSiPixelChargeReweightingAlgorithm_->hitSignalReweight<digitizerUtility::Ph2Amplitude>(hit,
hit_signal,
hitIndex,
referenceIndex4CR,
tofBin,
topol,
detID,
theSignal,
hit.processType(),
makeDigiSimLinks_);
} else {
// If it's not the primary particle, use the first hit in the collection as SimHit, which should be the corresponding primary.
referenceIndex4CR = firstHitIndex;
reweighted =
theSiPixelChargeReweightingAlgorithm_->hitSignalReweight<digitizerUtility::Ph2Amplitude>((*inputBegin),
hit_signal,
hitIndex,
referenceIndex4CR,
tofBin,
topol,
detID,
theSignal,
hit.processType(),
makeDigiSimLinks_);
}
}

float corr_time = hit.tof() - pixdet->surface().toGlobal(hit.localPosition()).mag() * c_inv;
for (auto const& hit_s : hit_signal) {
int chan = hit_s.first;
theSignal[chan] += (makeDigiSimLinks_ ? digitizerUtility::Ph2Amplitude(
hit_s.second, &hit, hit_s.second, corr_time, hitIndex, tofBin)
: digitizerUtility::Ph2Amplitude(hit_s.second, nullptr, hit_s.second));
if (!reweighted) {
for (auto const& hit_s : hit_signal) {
int chan = hit_s.first;
theSignal[chan] += (makeDigiSimLinks_ ? digitizerUtility::Ph2Amplitude(
hit_s.second, &hit, hit_s.second, corr_time, hitIndex, tofBin)
: digitizerUtility::Ph2Amplitude(hit_s.second, nullptr, hit_s.second));
}
}
}
} // end of induce_signal function

// ======================================================================
//
// Add electronic noise to pixel charge
Expand All @@ -587,6 +628,7 @@ void Phase2TrackerDigitizerAlgorithm::add_noise(const Phase2TrackerGeomDetUnit*
s.second += noise;
}
}

// ======================================================================
//
// Add Cross-talk contribution
Expand Down Expand Up @@ -719,6 +761,7 @@ void Phase2TrackerDigitizerAlgorithm::pixel_inefficiency(const SubdetEfficiencie
}
}
}

void Phase2TrackerDigitizerAlgorithm::initializeEvent(CLHEP::HepRandomEngine& eng) {
if (addNoise_ || addPixelInefficiency_ || fluctuateCharge_ || addThresholdSmearing_) {
gaussDistribution_ = std::make_unique<CLHEP::RandGaussQ>(eng, 0., theNoiseInElectrons_);
Expand All @@ -742,7 +785,7 @@ void Phase2TrackerDigitizerAlgorithm::initializeEvent(CLHEP::HepRandomEngine& en
// Configurations for barrel and foward pixels possess different tanLorentzAngleperTesla
// parameter value

LocalVector Phase2TrackerDigitizerAlgorithm::DriftDirection(const Phase2TrackerGeomDetUnit* pixdet,
LocalVector Phase2TrackerDigitizerAlgorithm::driftDirection(const Phase2TrackerGeomDetUnit* pixdet,
const GlobalVector& bfield,
const DetId& detId) const {
Frame detFrame(pixdet->surface().position(), pixdet->surface().rotation());
Expand All @@ -759,6 +802,11 @@ LocalVector Phase2TrackerDigitizerAlgorithm::DriftDirection(const Phase2TrackerG
// Read Lorentz angle from DB:
if (use_LorentzAngle_DB_) {
bool isPixel = (Sub_detid == PixelSubdetector::PixelBarrel || Sub_detid == PixelSubdetector::PixelEndcap);
if (isPixel) {
LogDebug("Phase2TrackerDigitizerAlgorithm") << "Read Lorentz angle from DB for Phase-2 IT";
} else {
LogDebug("Phase2TrackerDigitizerAlgorithm") << "Read Lorentz angle from DB for Phase-2 OT";
}

float lorentzAngle =
isPixel ? siPixelLorentzAngle_->getLorentzAngle(detId) : siPhase2OTLorentzAngle_->getLorentzAngle(detId);
Expand All @@ -770,6 +818,7 @@ LocalVector Phase2TrackerDigitizerAlgorithm::DriftDirection(const Phase2TrackerG
scale = (1 + alpha2 * std::pow(Bfield.z(), 2));
} else {
// Read Lorentz angle from cfg file:
LogDebug("Phase2TrackerDigitizerAlgorithm") << "Read Lorentz angle from cfg file";
float alpha2_Endcap = 0.0;
float alpha2_Barrel = 0.0;
if (alpha2Order_) {
Expand Down Expand Up @@ -967,6 +1016,7 @@ int Phase2TrackerDigitizerAlgorithm::convertSignalToAdc(uint32_t detID, float si
}
return signal_in_adc;
}

float Phase2TrackerDigitizerAlgorithm::calcQ(float x) {
constexpr float p1 = 12.5f;
constexpr float p2 = 0.2733f;
Expand Down
@@ -1,5 +1,5 @@
#ifndef __SimTracker_SiPhase2Digitizer_Phase2TrackerDigitizerAlgorithm_h
#define __SimTracker_SiPhase2Digitizer_Phase2TrackerDigitizerAlgorithm_h
#ifndef SimTracker_SiPhase2Digitizer_Phase2TrackerDigitizerAlgorithm_h
#define SimTracker_SiPhase2Digitizer_Phase2TrackerDigitizerAlgorithm_h

#include <map>
#include <memory>
Expand All @@ -14,6 +14,7 @@
#include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"

#include "SimTracker/Common/interface/DigitizerUtility.h"
#include "SimTracker/Common/interface/SiPixelChargeReweightingAlgorithm.h"
#include "SimTracker/SiPhase2Digitizer/plugins/Phase2TrackerDigitizerFwd.h"

// Units and Constants
Expand All @@ -39,6 +40,7 @@ class SiPixelQuality;
class SiPhase2OuterTrackerLorentzAngle;
class TrackerGeometry;
class TrackerTopology;
class SiPixelChargeReweightingAlgorithm;

// REMEMBER CMS conventions:
// -- Energy: GeV
Expand Down Expand Up @@ -168,6 +170,11 @@ class Phase2TrackerDigitizerAlgorithm {
const double pseudoRadDamage_; // Decrease the amount off freed charge that reaches the collector
const double pseudoRadDamageRadius_; // Only apply pseudoRadDamage to pixels with radius<=pseudoRadDamageRadius

// charge reweighting
const bool useChargeReweighting_;
// access 2D templates from DB. Only gets initialized if useChargeReweighting_ is set to true
const std::unique_ptr<SiPixelChargeReweightingAlgorithm> theSiPixelChargeReweightingAlgorithm_;

// The PDTable
// HepPDTable *particleTable;
// ParticleDataTable *particleTable;
Expand All @@ -190,8 +197,10 @@ class Phase2TrackerDigitizerAlgorithm {
const Phase2TrackerGeomDetUnit* pixdet,
const GlobalVector& bfield,
const std::vector<digitizerUtility::EnergyDepositUnit>& ionization_points) const;
virtual void induce_signal(const PSimHit& hit,
virtual void induce_signal(std::vector<PSimHit>::const_iterator inputBegin,
const PSimHit& hit,
const size_t hitIndex,
const size_t firstHitIndex,
const uint32_t tofBin,
const Phase2TrackerGeomDetUnit* pixdet,
const std::vector<digitizerUtility::SignalPoint>& collection_points);
Expand All @@ -208,7 +217,8 @@ class Phase2TrackerDigitizerAlgorithm {

// access to the gain calibration payloads in the db. Only gets initialized if check_dead_pixels_ is set to true.
const std::unique_ptr<SiPixelGainCalibrationOfflineSimService> theSiPixelGainCalibrationService_;
LocalVector DriftDirection(const Phase2TrackerGeomDetUnit* pixdet,

LocalVector driftDirection(const Phase2TrackerGeomDetUnit* pixdet,
const GlobalVector& bfield,
const DetId& detId) const;

Expand Down