Skip to content

[DEE] Update AnalyzeEvent in MSubModuleChargeTransport #85

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.

Sign up for GitHub

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

Open
fhagemann wants to merge 2 commits into
base: develop/em
Choose a base branch
from
Open

[DEE] Update AnalyzeEvent in MSubModuleChargeTransport #85

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7583325
Use parameterization derived from Boggs (2023)
fhagemann Oct 30, 2025
7814282
Add cross-talk based on Bobby's parameterization
fhagemann Jan 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
191 changes: 184 additions & 7 deletions src/MSubModuleChargeTransport.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,10 @@
#include "MSubModuleChargeTransport.h"

// Standard libs:
#include <cmath>

// ROOT libs:
#include "TMath.h"

// MEGAlib libs:
#include "MSubModule.h"
Expand Down Expand Up @@ -156,9 +158,23 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event)
{
// Main data analysis routine, which updates the event to a new level

// Dummy code:
// Define physical constants
constexpr double kB = TMath::K(); // unit: J/K
constexpr double ElementaryCharge = TMath::Qe(); // unit: C
constexpr double IonizationEnergy = 0.00295; // unit: keV
constexpr double Epsilon0 = 8.85418781762039e-14; // unit: F/cm
constexpr double EpsilonR = 16.0; // in germanium, unitless

// TODO: Read bias voltage and temperature of the detector from a file
constexpr double BiasVoltage = 1050.0; // unit: V
constexpr double Temperature = 87.0; // unit: K

// TODO: Implement energy-dependent initial charge-cloud sizes
constexpr double InitialChargeCloudSize = 0.02; // 200µm in cm

// Create strip hits
list<MDEEStripHit> ChargeTransportLVHits;

list<MDEEStripHit>& LVHits = Event->GetDEEStripHitLVListReference();
for (MDEEStripHit& SH: LVHits) {
MVector Pos = SH.m_SimulatedPositionInDetector;
Expand All @@ -167,25 +183,108 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event)
unsigned int DetID = SH.m_ROE.GetDetectorID();
double XWidth = m_XWidths[DetID];
double YWidth = m_YWidths[DetID];
double Thickness = m_Thicknesses[DetID];
double XPitch = m_XPitches[DetID];
double Radius = m_Radii[DetID];
int NXStrips = m_NXStrips[DetID];

// TODO: replace magic numbers for cross-talk by reading parameters from file
double MeanElectricField = BiasVoltage / Thickness; // unit: V/cm
double N = SH.m_SimulatedEnergy / IonizationEnergy;
constexpr double A = 0.020898810806838582;
constexpr double B1 = 0.2740445312420676;
constexpr double B2 = 0.22725576252704122;
constexpr double Z1 = -6.8702771019427935;
constexpr double Z2 = -6.403253988164561;
double CrossTalk = A * std::max(1.0 - std::exp(-B1 * (10.0 * Pos.Z() - Z1)), 1.0 - std::exp(B2 * (10.0 * Pos.Z() - Z2)));

// Calculate LV strip ID by rounding down intentionally to avoid truncation towards zero
int ID = static_cast<int>(std::floor((Pos.X() + XWidth/2.0) / XPitch));

// Check for strip ID and if the position is within the allowed strip length or on the guard ring
// TODO: Confirm the correct boundary of the guard ring based on SMEX detector models
if (ID >= 0 && ID < NXStrips && std::abs(Pos.Y()) <= YWidth/2.0 && std::hypot(Pos.X(), Pos.Y()) <= Radius) {
SH.m_ROE.SetStripID(ID);
SH.m_IsGuardRing = false;

// Apply charge sharing based on relative coordinate to the gap of that strip (0 <= X < XPitch)
double XFromGap = std::fmod(Pos.X() + XWidth/2.0, XPitch);

// calculate σ and µ, assuming t = z / v = z / (µ * E)
double SigmaX = std::sqrt(2.0 * kB * (Pos.Z() + Thickness / 2.0) / (ElementaryCharge * MeanElectricField)); // in cm
double EtaX = std::cbrt(std::pow(InitialChargeCloudSize, 3) + 3.0 * N * ElementaryCharge * (Pos.Z() + Thickness / 2.0) / (4.0 * TMath::Pi() * Epsilon0 * EpsilonR * MeanElectricField)); // in cm
// cout << "X: " << Pos.X() << ", Y: " << Pos.Y() << ", Z: " << Pos.Z() << ", E: " << SH.m_SimulatedEnergy << endl;
// cout << "SigmaX: " << SigmaX << ", EtaX: " << EtaX << endl;

auto Lambda = [&](double x) -> double {
double a = (x - EtaX) / (TMath::Sqrt2() * SigmaX);
double b = (x + EtaX) / (TMath::Sqrt2() * SigmaX);
return SH.m_SimulatedEnergy * (1.0 - CrossTalk) / (8.0 * std::pow(EtaX, 3)) * (
std::erf(b) * (2.0 * std::pow(EtaX, 3) + x * (3.0 * std::pow(EtaX, 2) - 3.0 * std::pow(SigmaX, 2) - std::pow(x, 2))) +
std::erf(a) * (2.0 * std::pow(EtaX, 3) - x * (3.0 * std::pow(EtaX, 2) - 3.0 * std::pow(SigmaX, 2) - std::pow(x, 2))) +
std::exp(- std::pow(b,2)) * std::sqrt(2.0 / TMath::Pi()) * SigmaX * (EtaX * x + (2.0 * std::pow(EtaX, 2) - 2.0 * std::pow(SigmaX, 2) - std::pow(x, 2))) +
std::exp(- std::pow(a,2)) * std::sqrt(2.0 / TMath::Pi()) * SigmaX * (EtaX * x - (2.0 * std::pow(EtaX, 2) - 2.0 * std::pow(SigmaX, 2) - std::pow(x, 2)))
);
};

double MainStripEnergy = Lambda(XPitch - XFromGap) - Lambda(-XFromGap) + SH.m_SimulatedEnergy * CrossTalk;
double NNLeftStripEnergy = Lambda(-XFromGap) - Lambda(-XPitch - XFromGap) + SH.m_SimulatedEnergy * CrossTalk;
double NNRightStripEnergy = Lambda(2.0*XPitch - XFromGap) - Lambda(XPitch - XFromGap) + SH.m_SimulatedEnergy * CrossTalk;

// cout << "Energy: " << SH.m_SimulatedEnergy << ", split into " << NNLeftStripEnergy << ", " << MainStripEnergy << " and " << NNRightStripEnergy << endl;

// create entry for the main hit
MDEEStripHit MainSH = SH;
MainSH.m_ROE.SetStripID(ID);
MainSH.m_Energy = MainStripEnergy;
MainSH.m_IsGuardRing = false;
ChargeTransportLVHits.push_back(MainSH);

// create MDEEStripHit for the left NN
if (NNLeftStripEnergy > IonizationEnergy) {
MDEEStripHit NNLeftSH = SH;
NNLeftSH.m_Energy = NNLeftStripEnergy;
if (ID > 0) {
NNLeftSH.m_ROE.SetStripID(ID - 1);
NNLeftSH.m_IsGuardRing = false;
// NNLeftSH.m_IsNearestNeighbor = true;
} else {
NNLeftSH.m_ROE.SetStripID(NXStrips);
NNLeftSH.m_IsGuardRing = true;
}
ChargeTransportLVHits.push_back(NNLeftSH);
}

// create MDEEStripHit for the right NN
if (NNRightStripEnergy > IonizationEnergy) {
MDEEStripHit NNRightSH = SH;
NNRightSH.m_Energy = NNRightStripEnergy;
if (ID < NXStrips - 1) {
NNRightSH.m_ROE.SetStripID(ID + 1);
NNRightSH.m_IsGuardRing = false;
// NNRightSH.m_IsNearestNeighbor = true;
} else {
NNRightSH.m_ROE.SetStripID(NXStrips);
NNRightSH.m_IsGuardRing = true;
}
ChargeTransportLVHits.push_back(NNRightSH);
}

} else {
// TODO: implement charge sharing also for GR events
SH.m_Energy = SH.m_SimulatedEnergy;
SH.m_ROE.SetStripID(NXStrips);
SH.m_IsGuardRing = true;
ChargeTransportLVHits.push_back(SH);
}
SH.m_Energy = SH.m_SimulatedEnergy;
}

// replace old list by new list
Event->GetDEEStripHitLVListReference().clear();
for (MDEEStripHit& SH: ChargeTransportLVHits) {
Event->AddDEEStripHitLV(SH);
}

list<MDEEStripHit> ChargeTransportHVHits;

list<MDEEStripHit>& HVHits = Event->GetDEEStripHitHVListReference();
for (MDEEStripHit& SH: HVHits) {
MVector Pos = SH.m_SimulatedPositionInDetector;
Expand All @@ -194,27 +293,105 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event)
unsigned int DetID = SH.m_ROE.GetDetectorID();
double XWidth = m_XWidths[DetID];
double YWidth = m_YWidths[DetID];
double Thickness = m_Thicknesses[DetID];
double YPitch = m_YPitches[DetID];
double Radius = m_Radii[DetID];
int NYStrips = m_NYStrips[DetID];

// TODO: replace magic numbers for cross-talk by reading parameters from file
double MeanElectricField = BiasVoltage / Thickness; // unit: V/cm
double N = SH.m_SimulatedEnergy / IonizationEnergy;
constexpr double A = 0.016119326831437686;
constexpr double B1 = 0.12520563605264975;
constexpr double B2 = 0.17122757421051746;
constexpr double Z1 = 5.812039914491831;
constexpr double Z2 = 6.710904485197925;
double CrossTalk = A * std::max(1.0 - std::exp(-B1 * (10.0 * Pos.Z() - Z1)), 1.0 - std::exp(B2 * (10.0 * Pos.Z() - Z2)));

// Calculate HV strip ID by rounding down intentionally to avoid truncation towards zero
// TODO: Confirm the correct strip pitch based on SMEX detector models
int ID = static_cast<int>(std::floor((Pos.Y() + YWidth/2.0) / YPitch));

// Check for strip ID and if the position is within the allowed strip length or on the guard ring
// TODO: Confirm the correct boundary of the guard ring based on SMEX detector models
if (ID >= 0 && ID < NYStrips && std::abs(Pos.X()) <= XWidth/2.0 && std::hypot(Pos.X(), Pos.Y()) <= Radius) {
SH.m_ROE.SetStripID(ID);
SH.m_IsGuardRing = false;

// Apply charge sharing based on relative coordinate to the gap of that strip (0 <= Y < YPitch)
double YFromGap = std::fmod(Pos.Y() + YWidth/2.0, YPitch);

// calculate σ and µ, assuming t = z / v = z / (µ * E)
double SigmaY = std::sqrt(2.0 * kB * Temperature * (Thickness / 2.0 - Pos.Z()) / (ElementaryCharge * MeanElectricField)); // in cm
double EtaY = std::cbrt(std::pow(InitialChargeCloudSize, 3) + 3.0 * N * ElementaryCharge * (Thickness / 2.0 - Pos.Z()) / (4.0 * TMath::Pi() * Epsilon0 * EpsilonR * MeanElectricField)); // in cm

auto Lambda = [&](double y) -> double {
double a = (y - EtaY) / (TMath::Sqrt2() * SigmaY);
double b = (y + EtaY) / (TMath::Sqrt2() * SigmaY);
return SH.m_SimulatedEnergy * (1.0 - CrossTalk) / (8.0 * std::pow(EtaY, 3)) * (
std::erf(b) * (2.0 * std::pow(EtaY, 3) + y * (3.0 * std::pow(EtaY, 2) - 3.0 * std::pow(SigmaY, 2) - std::pow(y, 2))) +
std::erf(a) * (2.0 * std::pow(EtaY, 3) - y * (3.0 * std::pow(EtaY, 2) - 3.0 * std::pow(SigmaY, 2) - std::pow(y, 2))) +
std::exp(- std::pow(b,2)) * std::sqrt(2 / TMath::Pi()) * SigmaY * (EtaY * y + (2.0 * std::pow(EtaY, 2) - 2.0 * std::pow(SigmaY, 2) - std::pow(y, 2))) +
std::exp(- std::pow(a,2)) * std::sqrt(2 / TMath::Pi()) * SigmaY * (EtaY * y - (2.0 * std::pow(EtaY, 2) - 2.0 * std::pow(SigmaY, 2) - std::pow(y, 2)))
);
};

double MainStripEnergy = Lambda(YPitch - YFromGap) - Lambda(-YFromGap) + SH.m_SimulatedEnergy * CrossTalk;
double NNLeftStripEnergy = Lambda(-YFromGap) - Lambda(-YPitch - YFromGap) + SH.m_SimulatedEnergy * CrossTalk;
double NNRightStripEnergy = Lambda(2.0*YPitch - YFromGap) - Lambda(YPitch - YFromGap) + SH.m_SimulatedEnergy * CrossTalk;

// create entry for the main hit
MDEEStripHit MainSH = SH;
MainSH.m_ROE.SetStripID(ID);
MainSH.m_Energy = MainStripEnergy;
MainSH.m_IsGuardRing = false;
ChargeTransportHVHits.push_back(MainSH);

// create MDEEStripHit for the left NN
if (NNLeftStripEnergy > IonizationEnergy) {
MDEEStripHit NNLeftSH = SH;
NNLeftSH.m_Energy = NNLeftStripEnergy;
if (ID > 0) {
NNLeftSH.m_ROE.SetStripID(ID - 1);
NNLeftSH.m_IsGuardRing = false;
// NNLeftSH.m_IsNearestNeighbor = true;
} else {
NNLeftSH.m_ROE.SetStripID(NYStrips);
NNLeftSH.m_IsGuardRing = true;
}
ChargeTransportHVHits.push_back(NNLeftSH);
}

// create MDEEStripHit for the right NN
if (NNRightStripEnergy > IonizationEnergy) {
MDEEStripHit NNRightSH = SH;
NNRightSH.m_Energy = NNRightStripEnergy;
if (ID < NYStrips - 1) {
NNRightSH.m_ROE.SetStripID(ID + 1);
NNRightSH.m_IsGuardRing = false;
// NNRightSH.m_IsNearestNeighbor = true;
} else {
NNRightSH.m_ROE.SetStripID(NYStrips);
NNRightSH.m_IsGuardRing = true;
}
ChargeTransportHVHits.push_back(NNRightSH);
}

} else {
// TODO: implement charge sharing also for GR events
SH.m_Energy = SH.m_SimulatedEnergy;
SH.m_ROE.SetStripID(NYStrips);
SH.m_IsGuardRing = true;
ChargeTransportHVHits.push_back(SH);
}
SH.m_Energy = SH.m_SimulatedEnergy;
}

// replace old list by new list
Event->GetDEEStripHitHVListReference().clear();
for (MDEEStripHit& SH: ChargeTransportHVHits) {
Event->AddDEEStripHitHV(SH);
}

// Merge hits:
// TODO: how to deal with flags like "m_IsNearestNeighbor" etc. ?
for (auto IterLV1 = LVHits.begin(); IterLV1 != LVHits.end(); ++IterLV1) {
auto IterLV2 = std::next(IterLV1);
while (IterLV2 != LVHits.end()) {
Expand Down