diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index 4e67ba3..fe9c095 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -27,8 +27,10 @@ #include "MSubModuleChargeTransport.h" // Standard libs: +#include // ROOT libs: +#include "TMath.h" // MEGAlib libs: #include "MSubModule.h" @@ -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 ChargeTransportLVHits; + list& LVHits = Event->GetDEEStripHitLVListReference(); for (MDEEStripHit& SH: LVHits) { MVector Pos = SH.m_SimulatedPositionInDetector; @@ -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(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 ChargeTransportHVHits; + list& HVHits = Event->GetDEEStripHitHVListReference(); for (MDEEStripHit& SH: HVHits) { MVector Pos = SH.m_SimulatedPositionInDetector; @@ -194,10 +293,21 @@ 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(std::floor((Pos.Y() + YWidth/2.0) / YPitch)); @@ -205,16 +315,83 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) // 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()) {