From dfc712d2c5ff9d1b16c303b8a77c61fa9c53d5b3 Mon Sep 17 00:00:00 2001 From: Ankur Yadav Date: Thu, 23 Jul 2026 15:53:21 +0200 Subject: [PATCH 1/2] TPC: add 83mKr calibration generator with Geant4 ionisation-fluctuation support - TPCDetParam.UseGeant4Edep: use Geant4's own energy deposit for ionisation instead of Bethe-Bloch/NA49, with a configurable SpecialCutsGeV threshold - Relocates the Kr-83m decay generator (GeneratorKrDecay) into Detectors/TPC/simulation as a compiled Generator subclass, with a thin krGenerator.C macro for use as an o2-sim external generator - Adds plotCluster.C for viewing the resulting calibration spectrum --- .../base/include/TPCBase/ParameterDetector.h | 2 + Detectors/TPC/simulation/CMakeLists.txt | 12 + Detectors/TPC/simulation/README.md | 49 ++- .../include/TPCSimulation/Detector.h | 2 + .../include/TPCSimulation/GeneratorKrDecay.h | 105 +++++ Detectors/TPC/simulation/macro/krGenerator.C | 27 ++ Detectors/TPC/simulation/macro/plotCluster.C | 407 ++++++++++++++++++ Detectors/TPC/simulation/src/Detector.cxx | 99 +++-- .../TPC/simulation/src/GeneratorKrDecay.cxx | 360 ++++++++++++++++ .../TPC/simulation/src/TPCSimulationLinkDef.h | 1 + 10 files changed, 1025 insertions(+), 39 deletions(-) create mode 100644 Detectors/TPC/simulation/include/TPCSimulation/GeneratorKrDecay.h create mode 100644 Detectors/TPC/simulation/macro/krGenerator.C create mode 100644 Detectors/TPC/simulation/macro/plotCluster.C create mode 100644 Detectors/TPC/simulation/src/GeneratorKrDecay.cxx diff --git a/Detectors/TPC/base/include/TPCBase/ParameterDetector.h b/Detectors/TPC/base/include/TPCBase/ParameterDetector.h index e557a174ec70a..9daaff9c7c9ff 100644 --- a/Detectors/TPC/base/include/TPCBase/ParameterDetector.h +++ b/Detectors/TPC/base/include/TPCBase/ParameterDetector.h @@ -33,6 +33,8 @@ struct ParameterDetector : public o2::conf::ConfigurableParamHelper` to `--configKeyValues`. + +`TPCDetParam.UseGeant4Edep=1` is required: it switches `Detector::ProcessHits()` to use Geant4's +own energy deposit directly instead of the default Bethe-Bloch/NA49 sampling, which is what +correctly resolves the Kr-83m decay channels into their discrete energy peaks. diff --git a/Detectors/TPC/simulation/include/TPCSimulation/Detector.h b/Detectors/TPC/simulation/include/TPCSimulation/Detector.h index 507ba992e2715..10ff59360374c 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/Detector.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/Detector.h @@ -146,6 +146,8 @@ class Detector : public o2::base::DetImpl void PostTrack() override { ; } void PreTrack() override { ; } + void SetSpecialPhysicsCuts() override; + void SetGeoFileName(const TString file) { mGeoFileName = file; } const TString& GetGeoFileName() const { return mGeoFileName; } diff --git a/Detectors/TPC/simulation/include/TPCSimulation/GeneratorKrDecay.h b/Detectors/TPC/simulation/include/TPCSimulation/GeneratorKrDecay.h new file mode 100644 index 0000000000000..df0ffd9d93b80 --- /dev/null +++ b/Detectors/TPC/simulation/include/TPCSimulation/GeneratorKrDecay.h @@ -0,0 +1,105 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file GeneratorKrDecay.h +/// \brief Generator for 83mKr decays, for TPC gain-map calibration simulation +/// \author Ankur Yadav + +#ifndef ALICEO2_TPC_GeneratorKrDecay_H_ +#define ALICEO2_TPC_GeneratorKrDecay_H_ + +#include "Generators/Generator.h" +#include +#include + +namespace o2 +{ +namespace tpc +{ + +// O2 status encoding from MCGenProperties.h +// bits 0-8: hepmc(9), bits 9-18: gen(10), bits 19-28: reserved(10), bits 29-31: sentinel=5 +inline int krO2EncodedStatus(int hepmc, int gen = 0) +{ + return (5 << 29) | ((gen & 0x3FF) << 9) | (hepmc & 0x1FF); +} + +namespace KrGenConfig +{ +inline double rInner = 83.5; +inline double rOuter = 246.5; // TPC outermost pad row outer edge ~247 cm; stay inside +inline double halfZ = 249.7; +inline int nPerEvent = 1000; // override at runtime via KR_N_PER_EVENT env var +} // namespace KrGenConfig + +struct KrProduct { + int pdg; + double eKin; +}; + +/// Table of 83mKr internal-conversion/gamma decay channels and their +/// branching fractions, derived at runtime from Geant4's level/gamma data +/// (falls back to hardcoded PhotonEvaporation5.7/z36.a83 values if +/// unavailable). Used by GeneratorKrDecay to sample one decay channel per +/// generated vertex. +class KrDecayTable +{ + public: + struct Channel { + double fraction; + int nProducts; + KrProduct products[6]; // max 5 used; 6 for safety + }; + // Eight physically motivated channels (T1 mode x T2 mode): + // T1: ICC_total=2035 -> 99.951% IC (75.163% outer-shell, 24.788% K-shell), 0.049% gamma + // K-shell: 65.2% K-fluorescence (Kalpha), 34.8% K-Auger + // T2: ICC_total=17.09 -> 94.472% IC, 5.528% gamma + // Source: G4 PhotonEvaporation5.7/z36.a83, RadioactiveDecay5.6/z36.a83 + static const int kNChannels = 8; + Channel channels[kNChannels]; + double cumulative[kNChannels]; + KrDecayTable(); + const Channel& sample() const; +}; + +} // namespace tpc +} // namespace o2 + +namespace o2 +{ +namespace eventgen +{ + +/// FairGenerator producing 83mKr decay vertices uniformly distributed in the +/// TPC drift volume, for gain-map/energy-resolution calibration simulation. +/// Each vertex emits the conversion electrons/photons of one randomly +/// sampled o2::tpc::KrDecayTable::Channel. +class GeneratorKrDecay : public Generator +{ + public: + GeneratorKrDecay(); + ~GeneratorKrDecay() override; + Bool_t Init() override; + Bool_t generateEvent() override; + Bool_t importParticles() override; + + private: + o2::tpc::KrDecayTable* mTable; + std::vector> mVertices; + // No ClassDefOverride: the base Generator class's dictionary is sufficient + // for a runtime-only generator that is never streamed via ROOT I/O + // (matches o2::eventgen::GeneratorGeantinos and other Generator subclasses). +}; + +} // namespace eventgen +} // namespace o2 + +#endif // ALICEO2_TPC_GeneratorKrDecay_H_ diff --git a/Detectors/TPC/simulation/macro/krGenerator.C b/Detectors/TPC/simulation/macro/krGenerator.C new file mode 100644 index 0000000000000..eb65dd0f21338 --- /dev/null +++ b/Detectors/TPC/simulation/macro/krGenerator.C @@ -0,0 +1,27 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file krGenerator +/// \brief This macro instantiates the compiled 83mKr TPC calibration +/// generator (o2::eventgen::GeneratorKrDecay), for use with +/// o2-sim -g external --extGenFile krGenerator.C --extGenFunc krGenerator +/// \author Ankur Yadav + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "FairGenerator.h" +#include "TPCSimulation/GeneratorKrDecay.h" +#endif + +FairGenerator* krGenerator() +{ + auto gen = new o2::eventgen::GeneratorKrDecay(); + return gen; +} diff --git a/Detectors/TPC/simulation/macro/plotCluster.C b/Detectors/TPC/simulation/macro/plotCluster.C new file mode 100644 index 0000000000000..1cef203b84b3c --- /dev/null +++ b/Detectors/TPC/simulation/macro/plotCluster.C @@ -0,0 +1,407 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file plotCluster.C +/// \brief Plots the merged 83mKr charge spectrum from tpcBoxClusters.root with a staged multi-Gaussian fit +/// \author Ankur Yadav + +// Reads tpcBoxClusters.root (produced by o2-tpc-krypton-clusterer) and plots +// the merged 83mKr charge spectrum with a staged multi-Gaussian fit. +// +// Usage: +// root -b -q 'plotCluster.C("tpcBoxClusters.root")' +// root -b -q 'plotCluster.C("tpcBoxClusters.root","IROC")' +// +// rocSel: "IROC" | "OROC1" | "OROC2" | "OROC3" | "OROC" | "ALL" +// sectorSel: -1 = all sectors merged (default); 0-35 = one sector +// qualityCuts: true (default) -- sigma-based cluster-shape quality cuts +// outputPdf: "" = auto-named kr__.pdf +// +// ROC row boundaries (O2 TPC, local pad row within sector): +// IROC rows 0 - 62 (63 rows) +// OROC1 rows 63 - 96 (34 rows) +// OROC2 rows 97 - 126 (30 rows) +// OROC3 rows 127 - 151 (25 rows) + +#include "DataFormatsTPC/KrCluster.h" +#include "TCanvas.h" +#include "TF1.h" +#include "TFile.h" +#include "TFitResult.h" +#include "TH1F.h" +#include "TLatex.h" +#include "TLine.h" +#include "TMath.h" +#include "TSystem.h" +#include "TTree.h" +#include +#include +#include +#include + +static const int kIrocLast = 62; +static const int kOroc1Last = 96; +static const int kOroc2Last = 126; +static const int kOroc3Last = 151; + +static int getRocIndex(float meanRow) +{ + int r = TMath::Nint(meanRow); + if (r <= kIrocLast) { + return 0; + } + if (r <= kOroc1Last) { + return 1; + } + if (r <= kOroc2Last) { + return 2; + } + if (r <= kOroc3Last) { + return 3; + } + return -1; +} + +static bool inRocSel(int rocIdx, const std::string& roc) +{ + if (roc == "IROC") { + return rocIdx == 0; + } + if (roc == "OROC1") { + return rocIdx == 1; + } + if (roc == "OROC2") { + return rocIdx == 2; + } + if (roc == "OROC3") { + return rocIdx == 3; + } + if (roc == "OROC") { + return rocIdx >= 1 && rocIdx <= 3; + } + return true; // ALL +} + +// Staged multi-Gaussian + exponential+erfc background fit of the 83mKr +// spectrum's six characteristic peaks (9.4, 12.6, 19.6, 29.1, 32.2, 41.6 keV). +// Parameters [0..15]: +// [0],[1] expo background (ln A, slope) +// [2] erfc shoulder amplitude (centre/width tied to the 41.6 keV peak) +// [3],[4] 9.4 keV Gaussian (amp, mu) +// [5],[6] 12.6 keV Gaussian (amp, mu) +// [7],[8] 19.6 keV Gaussian (amp, mu) +// [9],[10] 29.1 keV Gaussian (amp, mu) +// [11],[12] 32.2 keV Gaussian (amp, mu) +// [13..15] 41.6 keV Gaussian (main; amp, mu=p[14], sigma=p[15]) +// Width constraint: sigma_i = p[15] * sqrt(E_i/41.6); only p[15] is free +// (Fano-limited resolution: sigma/E ~ 1/sqrt(E), so sigma ~ sqrt(E)). +static void fitKrSpectrum(TH1F* h, double mu41) +{ + const double ratio[6] = {9.4 / 41.6, 12.6 / 41.6, 19.6 / 41.6, 29.1 / 41.6, 32.2 / 41.6, 1.0}; + const char* lbl[6] = {"T2#gamma 9.4 keV", "K#alpha 12.6 keV", "19.6 keV", "29.1 keV", "32.2 keV", "41.6 keV (main)"}; + const int col[6] = {kOrange + 2, kGreen + 2, kViolet + 1, kMagenta, kCyan + 2, kRed}; + // Satellite i (0..4) lives at parameters [ampIdx(i)], [muIdx(i)]; + // the main (41.6 keV) peak is amp=p[13], mu=p[14], sigma=p[15]. + auto ampIdx = [](int i) { return 3 + 2 * i; }; + auto muIdx = [](int i) { return 4 + 2 * i; }; + const int kAmpMain = 13, kMuMain = 14, kSigMain = 15; + const int kNPar = 16; + + double mu[6], sig0[6]; + for (int i = 0; i < 6; i++) { + mu[i] = mu41 * ratio[i]; + sig0[i] = mu[i] * 0.08; + } + + double amp41 = TMath::Max(h->GetBinContent(h->FindBin(mu41)), 1.); + const double alo[6] = {amp41 / 50., amp41 / 50., amp41 / 50., amp41 / 50., amp41 / 50., amp41 * 0.30}; + const double ahi[6] = {amp41 * 1.0, amp41 * 1.0, amp41 * 1.0, amp41 * 1.0, amp41 * 1.0, amp41 * 3.00}; + + const double slo[6] = {mu[0] * 0.05, mu[1] * 0.05, mu[2] * 0.05, mu[3] * 0.04, mu[4] * 0.04, mu[5] * 0.03}; + const double shi[6] = {mu[0] * 0.30, mu[1] * 0.28, mu[2] * 0.30, mu[3] * 0.25, mu[4] * 0.25, mu[5] * 0.20}; + + const double xlo = mu41 * 0.11, xhi = mu41 * 1.20; + + TF1* fbg = new TF1("fbg_pre", "expo", mu41 * 0.07, mu41 * 0.17); + h->Fit(fbg, "RQN0"); + + const double rlo[6] = {mu41 * 0.175, mu41 * 0.265, mu41 * 0.370, mu41 * 0.625, mu41 * 0.720, mu41 * 0.865}; + const double rhi[6] = {mu41 * 0.280, mu41 * 0.365, mu41 * 0.570, mu41 * 0.760, mu41 * 0.835, mu41 * 1.180}; + TF1* fg[6]; + for (int i = 0; i < 6; i++) { + fg[i] = new TF1(Form("fg_pre%d", i), "gaus", rlo[i], rhi[i]); + fg[i]->SetParameters(TMath::Max(h->GetBinContent(h->FindBin(mu[i])), 1.), mu[i], sig0[i]); + fg[i]->SetParLimits(0, alo[i], ahi[i]); + fg[i]->SetParLimits(1, mu[i] * 0.90, mu[i] * 1.10); + fg[i]->SetParLimits(2, slo[i], shi[i]); + h->Fit(fg[i], "RQN0"); + } + + TF1* total = new TF1( + "fitTotal", + [](double* x, double* p) -> double { + static const double r[5] = {9.4 / 41.6, 12.6 / 41.6, 19.6 / 41.6, 29.1 / 41.6, 32.2 / 41.6}; + double val = TMath::Exp(p[0] + p[1] * x[0]); + val += p[2] * TMath::Erfc((x[0] - p[14]) / (TMath::Sqrt2() * p[15])); + for (int i = 0; i < 5; i++) { + val += p[3 + 2 * i] * TMath::Gaus(x[0], p[4 + 2 * i], p[15] * TMath::Sqrt(r[i]), false); + } + val += p[13] * TMath::Gaus(x[0], p[14], p[15], false); + return val; + }, + xlo, xhi, kNPar); + total->SetNpx(3000); + total->SetParameter(0, fbg->GetParameter(0)); + total->SetParameter(1, fbg->GetParameter(1)); + double erfcSeed = h->GetBinContent(h->FindBin(mu41 * 0.905)) * 0.4; + total->SetParameter(2, erfcSeed > 1. ? erfcSeed : 1.); + for (int i = 0; i < 6; i++) { + double aSeed = TMath::Min(TMath::Max(fg[i]->GetParameter(0), alo[i]), ahi[i]); + total->SetParameter(i < 5 ? ampIdx(i) : kAmpMain, aSeed); + total->SetParameter(i < 5 ? muIdx(i) : kMuMain, fg[i]->GetParameter(1)); + } + total->SetParameter(kSigMain, TMath::Abs(fg[5]->GetParameter(2))); + + total->SetParLimits(1, -0.02, 0.); + total->SetParLimits(2, 0., 1e9); + for (int i = 0; i < 5; i++) { + total->SetParLimits(ampIdx(i), alo[i], ahi[i]); + total->SetParLimits(muIdx(i), mu[i] * 0.90, mu[i] * 1.10); + } + total->SetParLimits(kAmpMain, alo[5], ahi[5]); + total->SetParLimits(kMuMain, mu[5] * 0.90, mu[5] * 1.10); + total->SetParLimits(kSigMain, mu[5] * 0.03, mu[5] * 0.20); + + total->SetLineColor(kBlack); + total->SetLineWidth(2); + TFitResultPtr r = h->Fit(total, "RS"); + + // Explicitly draw the combined total fit (sum of background+erfc+all + // Gaussians) as its own visible curve -- the individual component draws + // below are mathematically identical pieces of this same function, but + // without this the combined shape isn't directly checkable by eye. + total->SetNpx(3000); + total->Draw("SAME"); + + TF1* fbgDraw = new TF1("fbg_draw", "exp([0]+[1]*x)+[2]*erfc((x-[3])/(sqrt(2.0)*[4]))", xlo, xhi); + fbgDraw->SetParameters(total->GetParameter(0), total->GetParameter(1), total->GetParameter(2), + total->GetParameter(kMuMain), total->GetParameter(kSigMain)); + fbgDraw->SetNpx(3000); + fbgDraw->SetLineColor(kGray + 1); + fbgDraw->SetLineStyle(7); + fbgDraw->SetLineWidth(2); + fbgDraw->Draw("SAME"); + + const double sigma41 = TMath::Abs(total->GetParameter(kSigMain)); + for (int i = 0; i < 6; i++) { + int ai = (i < 5) ? ampIdx(i) : kAmpMain; + int mi = (i < 5) ? muIdx(i) : kMuMain; + double sigmaI = sigma41 * TMath::Sqrt(ratio[i]); // = sigma41 * sqrt(E_i/41.6) + TF1* gc = new TF1(Form("gc%d", i), "[0]*exp(-0.5*pow((x-[1])/[2],2))", xlo, xhi); + gc->SetParameters(total->GetParameter(ai), total->GetParameter(mi), sigmaI); + gc->SetNpx(3000); + gc->SetLineColor(col[i]); + gc->SetLineStyle(2); + gc->SetLineWidth(2); + gc->Draw("SAME"); + double px = total->GetParameter(mi), py = total->GetParameter(ai); + if (py > h->GetMaximum() * 0.015) { + TLatex* tx = new TLatex(px + mu41 * 0.008, py * 0.65, lbl[i]); + tx->SetTextSize(0.028); + tx->SetTextColor(col[i]); + tx->Draw(); + } + } + + double chi2ndf = (r->Ndf() > 0) ? r->Chi2() / r->Ndf() : -1.; + printf("\n==== Kr-83m spectrum fit [%s] ====\n", h->GetTitle()); + printf(" 41.6 keV seed : %.1f ADC\n", mu41); + printf(" sigma_41 (free): %.1f +/- %.1f ADC (%.2f%%)\n", sigma41, total->GetParError(kSigMain), + 100. * sigma41 / total->GetParameter(kMuMain)); + printf(" chi2/ndf : %.2f\n", chi2ndf); + printf(" %-22s %14s %14s %8s\n", "Peak", "mu_fit [ADC]", "sigma (sqrtE)", "reso [%]"); + printf(" %-22s %14s %14s %8s\n", "----", "------------", "-------------", "--------"); + for (int i = 0; i < 6; i++) { + int mi = (i < 5) ? muIdx(i) : kMuMain; + double muI = total->GetParameter(mi); + double sigmaI = sigma41 * TMath::Sqrt(ratio[i]); + double reso = (muI > 0.) ? 100. * sigmaI / muI : -1.; + printf(" %-22s %7.1f +/- %4.1f %12.1f %7.2f%%\n", lbl[i], muI, total->GetParError(mi), sigmaI, reso); + } + printf("==============================================\n\n"); +} + +void plotCluster(const char* inputFile = "tpcBoxClusters.root", const char* rocSel = "IROC", int sectorSel = -1, + bool qualityCuts = true, const char* outputPdf = "") +{ + gSystem->Load("libO2TPCReconstruction"); + gSystem->Load("libO2DataFormatsTPC"); + gStyle->SetOptStat(0); // stat box hides the peaks otherwise + + std::string roc(rocSel); + for (auto& c : roc) { + c = toupper(c); + } + if (roc != "IROC" && roc != "OROC1" && roc != "OROC2" && roc != "OROC3" && roc != "OROC" && roc != "ALL") { + std::cout << "Invalid rocSel '" << rocSel << "'. Choose: IROC OROC1 OROC2 OROC3 OROC ALL" << std::endl; + return; + } + if (sectorSel < -1 || sectorSel > 35) { + std::cout << "Invalid sectorSel " << sectorSel << ". Use -1 (all) or 0-35." << std::endl; + return; + } + + auto f = TFile::Open(inputFile); + if (!f || f->IsZombie()) { + std::cout << "Cannot open " << inputFile << std::endl; + return; + } + auto t = (TTree*)f->Get("Clusters"); + if (!t) { + std::cout << "No 'Clusters' tree in " << inputFile << std::endl; + return; + } + if (!t->GetBranch("TPCBoxCluster_0")) { + std::cout << "Expected DPL branches (TPCBoxCluster_N) not found in " << inputFile << std::endl; + return; + } + + std::cout << "Input : " << inputFile << std::endl; + std::cout << "TFs : " << t->GetEntries() << std::endl; + std::cout << "ROC sel : " << roc << std::endl; + std::cout << "Sector : " << (sectorSel < 0 ? "all" : Form("%d", sectorSel)) << std::endl; + std::cout << "QC cuts : " << (qualityCuts ? "ON" : "OFF") << std::endl; + + auto passQC = [&](const o2::tpc::KrCluster& c) -> bool { + if (!qualityCuts) { + return true; + } + return (c.sigmaTime > 0.1f && c.sigmaTime < 1.8f) && + (c.sigmaRow > 0.2f && c.sigmaRow < 0.6f + c.totCharge / 4000.f) && (c.sigmaPad > 0.1f && c.sigmaPad < 1.2f); + }; + + std::vector* secCls[36] = {}; + for (int s = 0; s < 36; s++) { + t->SetBranchAddress(Form("TPCBoxCluster_%d", s), &secCls[s]); + } + + const int nBins = 400; + const double xMax = 6000.; + const double binW = xMax / nBins; // 15 ADC + std::string mergeTitle = + (sectorSel >= 0) ? Form("Sector %d -- %s", sectorSel, roc.c_str()) : Form("All sectors -- %s", roc.c_str()); + TH1F* hMerged = new TH1F("hMerged", Form("%s;Total cluster charge (ADC counts);Entries / %.0f ADC", mergeTitle.c_str(), binW), + nBins, 0., xMax); + hMerged->SetDirectory(nullptr); + + long long nTotal = 0, nQCPass = 0, nMerged = 0; + for (Long64_t ev = 0; ev < t->GetEntries(); ++ev) { + t->GetEntry(ev); + for (int s = 0; s < 36; s++) { + if (!secCls[s]) { + continue; + } + for (auto& c : *secCls[s]) { + ++nTotal; + if (!passQC(c)) { + continue; + } + int rocIdx = getRocIndex(c.meanRow); + if (rocIdx < 0) { + continue; + } + ++nQCPass; + if (inRocSel(rocIdx, roc) && (sectorSel < 0 || s == sectorSel)) { + hMerged->Fill(c.totCharge); + ++nMerged; + } + } + } + } + + printf("Total clusters read : %lld\n", nTotal); + printf("QC-passed clusters : %lld (%.1f%%)\n", nQCPass, nTotal > 0 ? 100. * nQCPass / nTotal : 0.); + printf("In canvas selection : %lld\n", nMerged); + + if (nMerged == 0) { + std::cout << "No clusters in canvas selection." << std::endl; + return; + } + + double mu41 = -1.; + { + int mBin = hMerged->FindBin(1000.); + for (int b = mBin + 1; b <= hMerged->GetNbinsX(); b++) { + if (hMerged->GetBinContent(b) > hMerged->GetBinContent(mBin)) { + mBin = b; + } + } + mu41 = hMerged->GetBinCenter(mBin); + } + const bool goodSpectrum = (mu41 >= 1200.); + printf("41.6 keV seed peak : %.0f ADC%s\n", mu41, goodSpectrum ? " [OK]" : " [WRONG SPECTRUM]"); + + std::string pdfName; + if (std::string(outputPdf).empty()) { + std::string secStr = (sectorSel < 0) ? "allsec" : Form("s%02d", sectorSel); + pdfName = Form("kr_%s_%s.pdf", roc.c_str(), secStr.c_str()); + for (auto& c : pdfName) { + c = tolower(c); + } + } else { + pdfName = outputPdf; + } + + TCanvas* cfit = new TCanvas("cfit", mergeTitle.c_str(), 1100, 700); + cfit->SetLeftMargin(0.10); + cfit->SetRightMargin(0.05); + cfit->SetBottomMargin(0.12); + + hMerged->SetLineColor(kBlue + 1); + hMerged->SetLineWidth(2); + hMerged->Draw("HIST"); + + const double ratio[6] = {9.4 / 41.6, 12.6 / 41.6, 19.6 / 41.6, 29.1 / 41.6, 32.2 / 41.6, 1.0}; + const int pcol[6] = {kOrange + 2, kGreen + 2, kViolet + 1, kMagenta, kCyan + 2, kRed}; + if (goodSpectrum) { + fitKrSpectrum(hMerged, mu41); + } else { + for (int i = 0; i < 6; i++) { + double xexp = mu41 * ratio[i]; + if (xexp < 50 || xexp > 5900) { + continue; + } + TLine* l = new TLine(xexp, 0, xexp, hMerged->GetMaximum() * 0.8); + l->SetLineColor(pcol[i]); + l->SetLineStyle(3); + l->SetLineWidth(1); + l->Draw(); + } + TLatex* msg = new TLatex(0.15, 0.85, Form("WRONG SPECTRUM: max at %.0f ADC (expected > 1200)", mu41)); + msg->SetNDC(); + msg->SetTextColor(kRed); + msg->SetTextSize(0.038); + msg->Draw(); + } + + TLatex* tlab = new TLatex(0.12, 0.92, Form("#bf{%s}", mergeTitle.c_str())); + tlab->SetNDC(); + tlab->SetTextSize(0.038); + tlab->Draw(); + TLatex* nlab = new TLatex(0.88, 0.92, Form("N_{sel}=%lld", nMerged)); + nlab->SetNDC(); + nlab->SetTextSize(0.033); + nlab->SetTextAlign(31); + nlab->Draw(); + + cfit->SaveAs(pdfName.c_str()); + printf("Saved %s\n", pdfName.c_str()); +} diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 8c9c45b6c13a8..5f1a52c7a2f40 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -193,41 +193,52 @@ Bool_t Detector::ProcessHits(FairVolume* vol) Int_t numberOfElectrons = 0; // I.H. - the type expected in addHit is short - // ---| Stepsize in cm |--- - const double stepSize = fMC->TrackStep(); - - double betaGamma = momentum.P() / fMC->TrackMass(); - betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg - - // ---| number of primary ionisations per cm |--- - const double primaryElectronsPerCM = - gasParam.Nprim * BetheBlochAleph(static_cast(betaGamma), gasParam.BetheBlochParam[0], - gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2], - gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]); - - // ---| mean number of collisions and random for this event |--- - const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM; - const int nColl = static_cast(fMC->GetRandom()->Poisson(meanNcoll)); - - // Variables needed to generate random powerlaw distributed energy loss - const double alpha_p1 = 1. - gasParam.Exp; // NA49/G3 value - const double oneOverAlpha_p1 = 1. / alpha_p1; - const double eMin = gasParam.Ipot; - const double eMax = gasParam.Eend; - const double kMin = TMath::Power(eMin, alpha_p1); - const double kMax = TMath::Power(eMax, alpha_p1); - const double wIon = gasParam.Wion; - - for (Int_t n = 0; n < nColl; n++) { - // Use GEANT3 / NA49 expression: - // P(eDep) ~ k * edep^-gasParam.getExp() - // eMin(~I) < eDep < eMax(300 electrons) - // k fixed so that Int_Emin^EMax P(Edep) = 1. - const double rndm = fMC->GetRandom()->Rndm(); - const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1); - int nel_step = static_cast(((eDep - eMin) / wIon) + 1); - nel_step = TMath::Min(nel_step, gasParam.MaxElePerStep); // 300 electrons corresponds to 10 keV - numberOfElectrons += nel_step; + // use Geant4 energy deposit directly for ionisation (Kr-83m calibration simulations) + if (detParam.UseGeant4Edep) { + // We have multiple collisions and add fluctuations: smear nel using + // gamma distr with mean = meanIon and variance = meanIon*FanoFactorG4. + // These parameters were tuned for GEANT4. + const double meanIon = fMC->Edep() / (gasParam.Wion * gasParam.ScaleFactorG4); + if (meanIon > 0.) { + numberOfElectrons = static_cast(gasParam.FanoFactorG4 * Gamma(meanIon / gasParam.FanoFactorG4)); + } + } else { + // ---| Stepsize in cm |--- + const double stepSize = fMC->TrackStep(); + + double betaGamma = momentum.P() / fMC->TrackMass(); + betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg + + // ---| number of primary ionisations per cm |--- + const double primaryElectronsPerCM = + gasParam.Nprim * BetheBlochAleph(static_cast(betaGamma), gasParam.BetheBlochParam[0], + gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2], + gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]); + + // ---| mean number of collisions and random for this event |--- + const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM; + const int nColl = static_cast(fMC->GetRandom()->Poisson(meanNcoll)); + + // Variables needed to generate random powerlaw distributed energy loss + const double alpha_p1 = 1. - gasParam.Exp; // NA49/G3 value + const double oneOverAlpha_p1 = 1. / alpha_p1; + const double eMin = gasParam.Ipot; + const double eMax = gasParam.Eend; + const double kMin = TMath::Power(eMin, alpha_p1); + const double kMax = TMath::Power(eMax, alpha_p1); + const double wIon = gasParam.Wion; + + for (Int_t n = 0; n < nColl; n++) { + // Use GEANT3 / NA49 expression: + // P(eDep) ~ k * edep^-gasParam.getExp() + // eMin(~I) < eDep < eMax(300 electrons) + // k fixed so that Int_Emin^EMax P(Edep) = 1. + const double rndm = fMC->GetRandom()->Rndm(); + const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1); + int nel_step = static_cast(((eDep - eMin) / wIon) + 1); + nel_step = TMath::Min(nel_step, gasParam.MaxElePerStep); // 300 electrons corresponds to 10 keV + numberOfElectrons += nel_step; + } } // LOG(info) << "tpc::AddHit" << FairLogger::endl << "Eloss: " @@ -3240,6 +3251,24 @@ std::string Detector::getHitBranchNames(int probe) const return std::string(); } +void Detector::SetSpecialPhysicsCuts() +{ + // lower energy threshold to track low-energy electrons for Kr-83m calibration + auto const& detParam = ParameterDetector::Instance(); + LOG(info) << "TPC SetSpecialPhysicsCuts: UseGeant4Edep=" << detParam.UseGeant4Edep; + if (detParam.UseGeant4Edep) { + auto& matmgr = o2::base::MaterialManager::Instance(); + const float specialCut = detParam.SpecialCutsGeV; + for (int med : {(int)kDriftGas1, (int)kDriftGas2, (int)kCO2}) { + matmgr.SpecialCut(GetName(), med, o2::base::ECut::kCUTELE, specialCut); + matmgr.SpecialCut(GetName(), med, o2::base::ECut::kCUTGAM, specialCut); + matmgr.SpecialCut(GetName(), med, o2::base::ECut::kDCUTE, specialCut); + matmgr.SpecialCut(GetName(), med, o2::base::ECut::kBCUTE, specialCut); + } + } + o2::base::Detector::SetSpecialPhysicsCuts(); +} + ClassImp(o2::tpc::Detector); // Define Factory method for calling from the outside diff --git a/Detectors/TPC/simulation/src/GeneratorKrDecay.cxx b/Detectors/TPC/simulation/src/GeneratorKrDecay.cxx new file mode 100644 index 0000000000000..e731d2d70bb1f --- /dev/null +++ b/Detectors/TPC/simulation/src/GeneratorKrDecay.cxx @@ -0,0 +1,360 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file GeneratorKrDecay.cxx +/// \brief Generator for 83mKr decays, for TPC gain-map calibration simulation +/// \author Ankur Yadav + +#include "TPCSimulation/GeneratorKrDecay.h" +#include "TDatabasePDG.h" +#include "TParticle.h" +#include "TParticlePDG.h" +#include "TRandom.h" +#include "TMath.h" +#include +#include +#include +#include +#include + +namespace o2 +{ +namespace tpc +{ + +// ── 83mKr decay physics ────────────────────────────────────────────────── +// +// Energies and ICC values are read at runtime from $G4LEVELGAMMADATA/z36.a83 +// (set automatically by Geant4 in any O2 alienv session). The hardcoded +// fallback values below are taken from PhotonEvaporation5.7/z36.a83 and are +// used only if the file cannot be opened or parsed. +// +// Atomic constants (NIST) — stable across G4 data releases, always hardcoded: +// Kr K-binding = 14.3256 keV +// Kr L1-binding = 1.9210 keV +// Kr Kα X-ray = 12.6000 keV +// Kr K-shell fluorescence yield ω_K = 0.652 (Bambynek et al.) +// +// From the G4 file we read for two levels: +// Level 2 (41.5569 keV, T1): E_gamma, ICC_total, K_shell_fraction +// Level 1 ( 9.4053 keV, T2): E_gamma, ICC_total +// Everything else is derived from these five numbers + the atomic constants. +// ───────────────────────────────────────────────────────────────────────── + +// Parse $G4LEVELGAMMADATA/z36.a83. +// Returns true and fills five values (energies in keV) on success. +static bool parseG4PhotonEvap(const char* path, + double& E_T1, // T1 gamma energy [keV] + double& ICC_T1, // T1 ICC_total + double& Kfrac_T1, // T1 K-shell fraction of ICC + double& E_T2, // T2 gamma energy [keV] + double& ICC_T2) // T2 ICC_total +{ + std::ifstream f(path); + if (!f.is_open()) { + return false; + } + + bool gotT1 = false, gotT2 = false; + bool wantT1 = false, wantT2 = false; + std::string line; + + while (std::getline(f, line)) { + if (line.empty()) { + continue; + } + std::istringstream ss(line); + int idx; + std::string tok; + double eLevel; + + // Header line: " N - E_level halflife ..." + if ((ss >> idx >> tok >> eLevel) && tok == "-") { + wantT1 = (idx == 2); // 41.5569 keV metastable state -> T1 transition + wantT2 = (idx == 1); // 9.4053 keV metastable state -> T2 transition + continue; + } + + if (!wantT1 && !wantT2) { + continue; + } + + // Transition line: " daughter E_gamma intensity multipolarity delta ICC_total K_frac ..." + ss.clear(); + ss.str(line); + int daughter, multi; + double Eg, inten, delta, icc, kfrac; + if (!(ss >> daughter >> Eg >> inten >> multi >> delta >> icc >> kfrac)) { + continue; + } + + if (wantT1) { + E_T1 = Eg; + ICC_T1 = icc; + Kfrac_T1 = kfrac; + gotT1 = true; + wantT1 = false; + } + if (wantT2) { + E_T2 = Eg; + ICC_T2 = icc; + gotT2 = true; + wantT2 = false; + } + + if (gotT1 && gotT2) { + break; + } + } + + if (!gotT1 || !gotT2) { + return false; + } + + // Sanity check — values far outside these ranges indicate a corrupt or wrong file + if (E_T1 < 25. || E_T1 > 40.) { + return false; + } + if (E_T2 < 5. || E_T2 > 15.) { + return false; + } + if (ICC_T1 < 100.) { + return false; + } + if (ICC_T2 < 5.) { + return false; + } + if (Kfrac_T1 < 0.1 || Kfrac_T1 > 0.5) { + return false; + } + + return true; +} + +KrDecayTable::KrDecayTable() +{ + // ── Atomic constants (NIST, keV, converted to GeV for ROOT) ────────── + static constexpr double kKbind = 14.3256e-6; // Kr K-shell binding + static constexpr double kL1bind = 1.9210e-6; // Kr L1-shell binding + static constexpr double kKalpha = 12.6000e-6; // Kr Kα X-ray + static constexpr double kKfluY = 0.652; // Kr K-shell fluorescence yield + + // ── Fallback values from PhotonEvaporation5.7/z36.a83 ──────────────── + double E_T1 = 32.1516e-6; // [GeV] T1 gamma energy + double ICC_T1 = 2035.0; + double Kfrac_T1 = 0.248; // fraction of ICC_T1 going through K-shell + double E_T2 = 9.4053e-6; // [GeV] T2 gamma energy + double ICC_T2 = 17.09; + + // ── Try to load from installed G4 data (keV in file → convert to GeV) ─ + const char* g4dir = std::getenv("G4LEVELGAMMADATA"); + if (g4dir) { + std::string path = std::string(g4dir) + "/z36.a83"; + double fE1, fICC1, fKf1, fE2, fICC2; + if (parseG4PhotonEvap(path.c_str(), fE1, fICC1, fKf1, fE2, fICC2)) { + E_T1 = fE1 * 1e-6; + ICC_T1 = fICC1; + Kfrac_T1 = fKf1; + E_T2 = fE2 * 1e-6; + ICC_T2 = fICC2; + std::cout << "[KrDecayTable] Loaded from " << path << "\n" + << " T1: E=" << fE1 << " keV ICC=" << ICC_T1 + << " K_frac=" << Kfrac_T1 << "\n" + << " T2: E=" << fE2 << " keV ICC=" << ICC_T2 << "\n"; + } else { + std::cout << "[KrDecayTable] WARNING: could not parse " << path + << " — using hardcoded fallback values\n"; + } + } else { + std::cout << "[KrDecayTable] WARNING: G4LEVELGAMMADATA not set" + << " — using hardcoded fallback values\n"; + } + + // ── Derived probabilities ───────────────────────────────────────────── + const double P_T1_g = 1.0 / (1.0 + ICC_T1); // T1 gamma + const double P_T1_K_IC = Kfrac_T1 * ICC_T1 / (1.0 + ICC_T1); // T1 K-shell IC + const double P_T1_out = ICC_T1 / (1.0 + ICC_T1) - P_T1_K_IC; // T1 outer-shell IC + const double P_T2_g = 1.0 / (1.0 + ICC_T2); // T2 gamma + const double P_T2_IC = ICC_T2 / (1.0 + ICC_T2); // T2 IC + + const double P_T1_Kf = P_T1_K_IC * kKfluY; // T1 K-IC → K-fluorescence + const double P_T1_Ka = P_T1_K_IC * (1.0 - kKfluY); // T1 K-IC → K-Auger + + // ── Particle kinetic energies ───────────────────────────────────────── + const double E_L_CE_T1 = E_T1 - kL1bind; // T1 L1-shell CE + const double E_K_CE = E_T1 - kKbind; // T1 K-shell CE + const double E_KLL = kKbind - 2.0 * kL1bind; // KLL Auger + const double E_res_aug = kKbind - E_KLL; // residual Auger (K-Auger path) + const double E_Laug_Kf = kKbind - kKalpha; // L-Auger after Kα emission + const double E_L_CE_T2 = E_T2 - kL1bind; // T2 L1-shell CE + + // ── Eight channels (T1-mode × T2-mode) ─────────────────────────────── + int i = 0; + + // Ch 0: T1 outer-IC + T2 IC → 41.6 keV local + channels[i] = {P_T1_out * P_T2_IC, 4, {{11, E_L_CE_T1}, {11, kL1bind}, {11, E_L_CE_T2}, {11, kL1bind}}}; + i++; + + // Ch 1: T1 outer-IC + T2 γ → 32.2 keV local + γ(T2) separate + channels[i] = {P_T1_out * P_T2_g, 3, {{11, E_L_CE_T1}, {11, kL1bind}, {22, E_T2}}}; + i++; + + // Ch 2: T1 K-IC + K-Auger + T2 IC → 41.6 keV local + channels[i] = {P_T1_Ka * P_T2_IC, 5, {{11, E_K_CE}, {11, E_KLL}, {11, E_res_aug}, {11, E_L_CE_T2}, {11, kL1bind}}}; + i++; + + // Ch 3: T1 K-IC + K-fluor + T2 IC → 29.1 keV local + Kα separate + channels[i] = {P_T1_Kf * P_T2_IC, 5, {{11, E_K_CE}, {11, E_Laug_Kf}, {22, kKalpha}, {11, E_L_CE_T2}, {11, kL1bind}}}; + i++; + + // Ch 4: T1 K-IC + K-fluor + T2 γ → 19.6 keV local + Kα + γ(T2) separate + channels[i] = {P_T1_Kf * P_T2_g, 4, {{11, E_K_CE}, {11, E_Laug_Kf}, {22, kKalpha}, {22, E_T2}}}; + i++; + + // Ch 5: T1 K-IC + K-Auger + T2 γ → 32.2 keV local + γ(T2) separate + channels[i] = {P_T1_Ka * P_T2_g, 4, {{11, E_K_CE}, {11, E_KLL}, {11, E_res_aug}, {22, E_T2}}}; + i++; + + // Ch 6: T1 γ + T2 IC → 9.4 keV local + γ(T1) separate + channels[i] = {P_T1_g * P_T2_IC, 3, {{22, E_T1}, {11, E_L_CE_T2}, {11, kL1bind}}}; + i++; + + // Ch 7: T1 γ + T2 γ → both photons escape + channels[i] = {P_T1_g * P_T2_g, 2, {{22, E_T1}, {22, E_T2}}}; + i++; + + double sum = 0.; + for (int j = 0; j < kNChannels; j++) { + sum += channels[j].fraction; + } + double cum = 0.; + for (int j = 0; j < kNChannels; j++) { + cum += channels[j].fraction / sum; + cumulative[j] = cum; + } +} + +const KrDecayTable::Channel& KrDecayTable::sample() const +{ + double r = gRandom->Uniform(); + for (int i = 0; i < kNChannels; i++) { + if (r <= cumulative[i]) { + return channels[i]; + } + } + return channels[kNChannels - 1]; +} + +} // namespace tpc +} // namespace o2 + +// ── GeneratorKrDecay ───────────────────────────────────────────────────── + +namespace o2 +{ +namespace eventgen +{ + +GeneratorKrDecay::GeneratorKrDecay() + : Generator("KrDecay", "83mKr TPC calibration source"), mTable(nullptr) +{ +} + +GeneratorKrDecay::~GeneratorKrDecay() { delete mTable; } + +Bool_t GeneratorKrDecay::Init() +{ + if (const char* env = std::getenv("KR_N_PER_EVENT")) { + int n = std::atoi(env); + if (n > 0) { + o2::tpc::KrGenConfig::nPerEvent = n; + } + } + std::cout << "[GeneratorKrDecay] Init:" + << " rInner=" << o2::tpc::KrGenConfig::rInner + << " rOuter=" << o2::tpc::KrGenConfig::rOuter + << " halfZ=" << o2::tpc::KrGenConfig::halfZ + << " nPerEvent=" << o2::tpc::KrGenConfig::nPerEvent << "\n"; + + mTable = new o2::tpc::KrDecayTable(); + setPositionUnit(1.0); // coords in cm + return Generator::Init(); +} + +Bool_t GeneratorKrDecay::generateEvent() +{ + mVertices.clear(); + // 1 cm safety margin inside field cage boundaries — avoids placing + // electrons exactly on sector boundaries which can cause hit coordinate + // transformation crashes in the merger when ROOT fills the TTree. + const double rInner = o2::tpc::KrGenConfig::rInner + 1.0; + const double rOuter = o2::tpc::KrGenConfig::rOuter - 1.0; + const double halfZ = o2::tpc::KrGenConfig::halfZ - 1.0; + const double r2Min = rInner * rInner; + const double r2Max = rOuter * rOuter; + for (int i = 0; i < o2::tpc::KrGenConfig::nPerEvent; ++i) { + double r = std::sqrt(gRandom->Uniform(r2Min, r2Max)); + double phi = gRandom->Uniform(0., TMath::TwoPi()); + double z = gRandom->Uniform(-halfZ, halfZ); + mVertices.push_back({{r * std::cos(phi), r * std::sin(phi), z}}); + } + return kTRUE; +} + +Bool_t GeneratorKrDecay::importParticles() +{ + mParticles.clear(); + // Reserve before any push_back to prevent std::vector reallocation. + // TParticle inherits from TObject (ROOT memory pool) and is not safe + // to move-construct via std::vector reallocation on macOS arm64 — + // ROOT's TStorage bookkeeping gets corrupted, causing malloc failures + // ~50-100 events later. Reserving eliminates all reallocations. + mParticles.reserve(o2::tpc::KrGenConfig::nPerEvent * o2::tpc::KrDecayTable::kNChannels); + + const int status = o2::tpc::krO2EncodedStatus(1, 0); + + for (size_t iv = 0; iv < mVertices.size(); ++iv) { + double vx = mVertices[iv][0]; + double vy = mVertices[iv][1]; + double vz = mVertices[iv][2]; + + const o2::tpc::KrDecayTable::Channel& ch = mTable->sample(); + for (int ip = 0; ip < ch.nProducts; ++ip) { + int pdg = ch.products[ip].pdg; + double eKin = ch.products[ip].eKin; + if (eKin < 0.1e-6) { + continue; + } + + double mass = (pdg == 11) ? 0.000511 : 0.0; + double E = eKin + mass; + double pmag = std::sqrt(std::max(0., E * E - mass * mass)); + double cosT = gRandom->Uniform(-1., 1.); + double sinT = std::sqrt(1. - cosT * cosT); + double phi = gRandom->Uniform(0., TMath::TwoPi()); + + TParticle part(pdg, status, -1, -1, -1, -1, + pmag * sinT * std::cos(phi), + pmag * sinT * std::sin(phi), + pmag * cosT, + E, vx, vy, vz, 0.); + mParticles.push_back(part); + // kToBeDone=BIT(16), kPrimary=BIT(17) + // Must be set AFTER push_back — copy constructor resets fBits + mParticles.back().SetBit(BIT(16)); + mParticles.back().SetBit(BIT(17)); + } + } + return kTRUE; +} + +} // namespace eventgen +} // namespace o2 diff --git a/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h b/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h index 6362b32c217f8..5fbdbf8f18342 100644 --- a/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h +++ b/Detectors/TPC/simulation/src/TPCSimulationLinkDef.h @@ -32,6 +32,7 @@ #pragma link C++ class o2::tpc::HitGroup + ; #pragma link C++ class o2::tpc::SAMPAProcessing + ; #pragma link C++ class o2::tpc::IDCSim + ; +#pragma link C++ class o2::eventgen::GeneratorKrDecay + ; #pragma link C++ class std::vector < o2::tpc::HitGroup> + ; From eb48bdb7479fae31e9001c5e5be780aef6bedfd5 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 23 Jul 2026 13:55:50 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- .../TPC/simulation/src/GeneratorKrDecay.cxx | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/Detectors/TPC/simulation/src/GeneratorKrDecay.cxx b/Detectors/TPC/simulation/src/GeneratorKrDecay.cxx index e731d2d70bb1f..96e7b451a3e41 100644 --- a/Detectors/TPC/simulation/src/GeneratorKrDecay.cxx +++ b/Detectors/TPC/simulation/src/GeneratorKrDecay.cxx @@ -52,11 +52,11 @@ namespace tpc // Parse $G4LEVELGAMMADATA/z36.a83. // Returns true and fills five values (energies in keV) on success. static bool parseG4PhotonEvap(const char* path, - double& E_T1, // T1 gamma energy [keV] - double& ICC_T1, // T1 ICC_total - double& Kfrac_T1, // T1 K-shell fraction of ICC - double& E_T2, // T2 gamma energy [keV] - double& ICC_T2) // T2 ICC_total + double& E_T1, // T1 gamma energy [keV] + double& ICC_T1, // T1 ICC_total + double& Kfrac_T1, // T1 K-shell fraction of ICC + double& E_T2, // T2 gamma energy [keV] + double& ICC_T2) // T2 ICC_total { std::ifstream f(path); if (!f.is_open()) { @@ -179,18 +179,18 @@ KrDecayTable::KrDecayTable() } // ── Derived probabilities ───────────────────────────────────────────── - const double P_T1_g = 1.0 / (1.0 + ICC_T1); // T1 gamma - const double P_T1_K_IC = Kfrac_T1 * ICC_T1 / (1.0 + ICC_T1); // T1 K-shell IC - const double P_T1_out = ICC_T1 / (1.0 + ICC_T1) - P_T1_K_IC; // T1 outer-shell IC - const double P_T2_g = 1.0 / (1.0 + ICC_T2); // T2 gamma - const double P_T2_IC = ICC_T2 / (1.0 + ICC_T2); // T2 IC + const double P_T1_g = 1.0 / (1.0 + ICC_T1); // T1 gamma + const double P_T1_K_IC = Kfrac_T1 * ICC_T1 / (1.0 + ICC_T1); // T1 K-shell IC + const double P_T1_out = ICC_T1 / (1.0 + ICC_T1) - P_T1_K_IC; // T1 outer-shell IC + const double P_T2_g = 1.0 / (1.0 + ICC_T2); // T2 gamma + const double P_T2_IC = ICC_T2 / (1.0 + ICC_T2); // T2 IC - const double P_T1_Kf = P_T1_K_IC * kKfluY; // T1 K-IC → K-fluorescence + const double P_T1_Kf = P_T1_K_IC * kKfluY; // T1 K-IC → K-fluorescence const double P_T1_Ka = P_T1_K_IC * (1.0 - kKfluY); // T1 K-IC → K-Auger // ── Particle kinetic energies ───────────────────────────────────────── - const double E_L_CE_T1 = E_T1 - kL1bind; // T1 L1-shell CE - const double E_K_CE = E_T1 - kKbind; // T1 K-shell CE + const double E_L_CE_T1 = E_T1 - kL1bind; // T1 L1-shell CE + const double E_K_CE = E_T1 - kKbind; // T1 K-shell CE const double E_KLL = kKbind - 2.0 * kL1bind; // KLL Auger const double E_res_aug = kKbind - E_KLL; // residual Auger (K-Auger path) const double E_Laug_Kf = kKbind - kKalpha; // L-Auger after Kα emission