From e0b8b91b6e434662b0b3b557c3a8017c8683800c Mon Sep 17 00:00:00 2001 From: Rocco Liotino Date: Tue, 21 Jul 2026 17:52:34 +0200 Subject: [PATCH 1/2] Update deuteron analysis tasks --- .../AOTTrack/PID/HMPID/hmpidDeuteron.cxx | 448 +++++++++++++----- .../AOTTrack/PID/HMPID/hmpidTableProducer.cxx | 153 ++++-- DPG/Tasks/AOTTrack/PID/HMPID/tableHMPID.h | 9 +- 3 files changed, 471 insertions(+), 139 deletions(-) diff --git a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidDeuteron.cxx b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidDeuteron.cxx index 27451cff7e3..38348abd7bd 100644 --- a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidDeuteron.cxx +++ b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidDeuteron.cxx @@ -12,17 +12,12 @@ #include "tableHMPID.h" #include -#include #include -#include -#include #include #include #include -#include - -#include +#include using namespace o2; using namespace o2::framework; @@ -33,7 +28,7 @@ struct HmpidDeuteron { Configurable nsigmaTPCMin{"nsigmaTPCMin", -3.0, "nsigmaTPCMin"}; Configurable nsigmaTPCMax{"nsigmaTPCMax", +3.0, "nsigmaTPCMax"}; Configurable nsigmaTOFMin{"nsigmaTOFMin", -3.0, "nsigmaTOFMin"}; - Configurable nsigmaTOFMax{"nsigmaTOFMax", +3.5, "nsigmaTOFMax"}; + Configurable nsigmaTOFMax{"nsigmaTOFMax", +3.0, "nsigmaTOFMax"}; Configurable minReqClusterITS{"minReqClusterITS", 4.0, "min number of clusters required in ITS"}; Configurable minTPCnClsFound{"minTPCnClsFound", 50.0f, "minTPCnClsFound"}; Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 70.0f, "min number of crossed rows TPC"}; @@ -42,57 +37,147 @@ struct HmpidDeuteron { Configurable maxDCAxy{"maxDCAxy", 0.5f, "maxDCAxy"}; Configurable maxDCAz{"maxDCAz", 0.5f, "maxDCAz"}; + Configurable cutMinMomGlobalTrack{"cutMinMomGlobalTrack", 1.5f, "minimum momentum of global track"}; + + // Configurables for plotting + // Momentum + Configurable momNBins{"momNBins", 490, "Momentum bins"}; + Configurable momMin{"momMin", 0.1f, "Momentum min"}; + Configurable momMax{"momMax", 5.0f, "Momentum max"}; + + // DeltaR + Configurable deltaRNBins{"deltaRNBins", 300, "DeltaR bins"}; + Configurable deltaRMin{"deltaRMin", 0.0f, "DeltaR min"}; + Configurable deltaRMax{"deltaRMax", 30.0f, "DeltaR max"}; + + // nSigma TPC + Configurable nSigmaTPCNBins{"nSigmaTPCNBins", 20, "nSigmaTPC bins"}; + Configurable nSigmaTPCMinPlot{"nSigmaTPCMinPlot", -5.0f, "nSigmaTPC min"}; + Configurable nSigmaTPCMaxPlot{"nSigmaTPCMaxPlot", 5.0f, "nSigmaTPC max"}; + + // nSigma TOF + Configurable nSigmaTOFNBins{"nSigmaTOFNBins", 20, "nSigmaTOF bins"}; + Configurable nSigmaTOFMinPlot{"nSigmaTOFMinPlot", -5.0f, "nSigmaTOF min"}; + Configurable nSigmaTOFMaxPlot{"nSigmaTOFMaxPlot", 5.0f, "nSigmaTOF max"}; + + // Charge + Configurable chargeNBins{"chargeNBins", 1000, "Charge bins"}; + Configurable chargeMin{"chargeMin", 200.0f, "Charge min"}; + Configurable chargeMax{"chargeMax", 2200.0f, "Charge max"}; + + // Eta + Configurable etaNBins{"etaNBins", 100, "Eta bins"}; + Configurable etaMin{"etaMin", -5.0f, "Eta min"}; + Configurable etaMax{"etaMax", 5.0f, "Eta max"}; + + // Phi + Configurable phiNBins{"phiNBins", 100, "Phi bins"}; + Configurable phiMin{"phiMin", -10.0f, "Phi min"}; + Configurable phiMax{"phiMax", 10.0f, "Phi max"}; + + // Mass + Configurable massNBins{"massNBins", 1000, "Mass bins"}; + Configurable massMin{"massMin", -10.0f, "Mass min"}; + Configurable massMax{"massMax", 10.0f, "Mass max"}; + + // variables for chamber_number and HVs/PCs + const int rich0 = 0, rich1 = 1, rich2 = 2, rich3 = 3, rich4 = 4, rich5 = 5, rich6 = 6; + static const int nCh = 7; + const float nmean = 1.28947; // meanIdxRad(); method from param not working, using value from HMPIDBase/Param.cxx + void init(InitContext const&) { + + // Axes for histograms + AxisSpec axisMom{momNBins, momMin, momMax, "#it{p} (GeV/#it{c})", "axisMom"}; + AxisSpec axisDeltaR{deltaRNBins, deltaRMin, deltaRMax, "#Delta R (cm)", "axisDeltaR"}; + AxisSpec axisNsigmaTOF{nSigmaTOFNBins, nSigmaTOFMinPlot, nSigmaTOFMaxPlot, "n#sigma_{TOF}", "axisNsigmaTOF"}; + AxisSpec axisNsigmaTPC{nSigmaTPCNBins, nSigmaTPCMinPlot, nSigmaTPCMaxPlot, "n#sigma_{TPC}", "axisNsigmaTPC"}; + AxisSpec axisCharge{chargeNBins, chargeMin, chargeMax, "Charge (ADC)", "axisCharge"}; + AxisSpec axisEta{etaNBins, etaMin, etaMax, "#eta", "axisEta"}; + AxisSpec axisPhi{phiNBins, phiMin, phiMax, "#phi", "axisPhi"}; + AxisSpec axisMass{massNBins, massMin, massMax, "mass^{2} (GeV/#it{c}^{2})", "axisMass"}; + AxisSpec axisChamber{10, -1.5, 8.5, "Chamber number", "axisChamber"}; + // Deuteron Pos - registryDA.add("incomingDe_Pos_8cm", "incomingDe_Pos_8cm", HistType::kTH1F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}}); - registryDA.add("incomingDe_Pos_4cm", "incomingDe_Pos_4cm", HistType::kTH1F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}}); - - registryDA.add("De_Pos_deltaR_8cm", "De_Pos_deltaR_8cm", HistType::kTH1F, {{300, 0.0, 30.0, "#Delta R (cm)"}}); - registryDA.add("De_Pos_deltaR_4cm", "De_Pos_deltaR_4cm", HistType::kTH1F, {{300, 0.0, 30.0, "#Delta R (cm)"}}); - - registryDA.add("survivingDe_Pos_8cm", "survivingDe_Pos_8cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {300, 0.0, 30.0, "#Delta R (cm)"}}); - registryDA.add("survivingDe_Pos_4cm", "survivingDe_Pos_4cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {300, 0.0, 30.0, "#Delta R (cm)"}}); - registryDA.add("De_Pos_Q_8cm", "De_Pos_Q_8cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {200, 0.0, 2000.0, "Q (ADC)"}}); - registryDA.add("De_Pos_Q_4cm", "De_Pos_Q_4cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {200, 0.0, 2000.0, "Q (ADC)"}}); - registryDA.add("De_Pos_ClsSize_8cm", "De_Pos_ClsSize_8cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, 0.0, 20., "Cls size"}}); - registryDA.add("De_Pos_ClsSize_4cm", "De_Pos_ClsSize_4cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, 0.0, 20.0, "Cls size"}}); + registryDA.add("incomingDe_Pos_8cm", "incomingDe_Pos_8cm", HistType::kTH1F, {axisMom}); + registryDA.add("incomingDe_Pos_4cm", "incomingDe_Pos_4cm", HistType::kTH1F, {axisMom}); + + registryDA.add("De_Pos_deltaR_8cm", "De_Pos_deltaR_8cm", HistType::kTH1F, {axisDeltaR}); + registryDA.add("De_Pos_deltaR_4cm", "De_Pos_deltaR_4cm", HistType::kTH1F, {axisDeltaR}); + + registryDA.add("nSigmaTOF_vs_momHMPID_DePos_8cm", "nSigmaTOF_vs_momHMPID_DePos_8cm", HistType::kTH2F, {axisMom, axisNsigmaTOF}); + registryDA.add("nSigmaTOF_vs_momHMPID_DePos_4cm", "nSigmaTOF_vs_momHMPID_DePos_4cm", HistType::kTH2F, {axisMom, axisNsigmaTOF}); + + registryDA.add("survivingDe_Pos_8cm", "survivingDe_Pos_8cm", HistType::kTH2F, {axisMom, axisDeltaR}); + registryDA.add("survivingDe_Pos_4cm", "survivingDe_Pos_4cm", HistType::kTH2F, {axisMom, axisDeltaR}); + + registryDA.add("De_Pos_Q_8cm", "De_Pos_Q_8cm", HistType::kTH2F, {axisMom, axisCharge}); + registryDA.add("De_Pos_Q_4cm", "De_Pos_Q_4cm", HistType::kTH2F, {axisMom, axisCharge}); + registryDA.add("De_Pos_momentum", "De_Pos_momentum", HistType::kTH2F, {{100, 0.0, 5.0, "#it{p}_{vtx} (GeV/#it{c})"}, {100, 0.0, 5.0, "#it{p}_{hmpid} (GeV/#it{c})"}}); - registryDA.add("nSigmaTPC_vs_momHMPID_noCut_DePos", "nSigmaTPC_vs_momHMPID_noCut_DePos", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5., 5.0, "n#sigma_TPC"}}); - registryDA.add("nSigmaTOF_vs_momHMPID_noCut_DePos", "nSigmaTOF_vs_momHMPID_noCut_DePos", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5., 5.0, "n#sigma_TOF"}}); - registryDA.add("nSigmaTPC_vs_momHMPID_Cut_DePos", "nSigmaTPC_vs_momHMPID_Cut_DePos", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5., 5.0, "n#sigma_TPC"}}); - registryDA.add("nSigmaTOF_vs_momHMPID_Cut_DePos", "nSigmaTOF_vs_momHMPID_Cut_DePos", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5., 5.0, "n#sigma_TOF"}}); + // nsigma plots + registryDA.add("nSigmaTPC_vs_momHMPID_noCut_DePos", "nSigmaTPC_vs_momHMPID_noCut_DePos", HistType::kTH2F, {axisMom, axisNsigmaTPC}); + registryDA.add("nSigmaTOF_vs_momHMPID_noCut_DePos", "nSigmaTOF_vs_momHMPID_noCut_DePos", HistType::kTH2F, {axisMom, axisNsigmaTOF}); + registryDA.add("nSigmaTPC_vs_momHMPID_Cut_DePos", "nSigmaTPC_vs_momHMPID_Cut_DePos", HistType::kTH2F, {axisMom, axisNsigmaTPC}); + registryDA.add("nSigmaTOF_vs_momHMPID_Cut_DePos", "nSigmaTOF_vs_momHMPID_Cut_DePos", HistType::kTH2F, {axisMom, axisNsigmaTOF}); // Deuteron Neg - registryDA.add("incomingDe_Neg_8cm", "incomingDe_Neg_8cm", HistType::kTH1F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}}); - registryDA.add("incomingDe_Neg_4cm", "incomingDe_Neg_4cm", HistType::kTH1F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}}); - - // plot aggiunti - registryDA.add("De_Neg_deltaR_8cm", "De_Neg_deltaR_8cm", HistType::kTH1F, {{300, 0.0, 30.0, "#Delta R (cm)"}}); - registryDA.add("De_Neg_deltaR_4cm", "De_Neg_deltaR_4cm", HistType::kTH1F, {{300, 0.0, 30.0, "#Delta R (cm)"}}); - - registryDA.add("survivingDe_Neg_8cm", "survivingDe_Neg_8cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {300, 0.0, 30.0, "#Delta R (cm)"}}); - registryDA.add("survivingDe_Neg_4cm", "survivingDe_Neg_4cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {300, 0.0, 30.0, "#Delta R (cm)"}}); - registryDA.add("De_Neg_Q_8cm", "De_Neg_Q_8cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {200, 0.0, 2000.0, "Q (ADC)"}}); - registryDA.add("De_Neg_Q_4cm", "De_Neg_Q_4cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {200, 0.0, 2000.0, "Q (ADC)"}}); - registryDA.add("De_Neg_ClsSize_8cm", "De_Neg_ClsSize_8cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, 0.0, 20.0, "Cls size"}}); - registryDA.add("De_Neg_ClsSize_4cm", "De_Neg_ClsSize_4cm", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, 0.0, 20.0, "Cls size"}}); + registryDA.add("incomingDe_Neg_8cm", "incomingDe_Neg_8cm", HistType::kTH1F, {axisMom}); + registryDA.add("incomingDe_Neg_4cm", "incomingDe_Neg_4cm", HistType::kTH1F, {axisMom}); + + registryDA.add("De_Neg_deltaR_8cm", "De_Neg_deltaR_8cm", HistType::kTH1F, {axisDeltaR}); + registryDA.add("De_Neg_deltaR_4cm", "De_Neg_deltaR_4cm", HistType::kTH1F, {axisDeltaR}); + + registryDA.add("nSigmaTOF_vs_momHMPID_DeNeg_8cm", "nSigmaTOF_vs_momHMPID_DeNeg_8cm", HistType::kTH2F, {axisMom, axisNsigmaTOF}); + registryDA.add("nSigmaTOF_vs_momHMPID_DeNeg_4cm", "nSigmaTOF_vs_momHMPID_DeNeg_4cm", HistType::kTH2F, {axisMom, axisNsigmaTOF}); + + registryDA.add("survivingDe_Neg_8cm", "survivingDe_Neg_8cm", HistType::kTH2F, {axisMom, axisDeltaR}); + registryDA.add("survivingDe_Neg_4cm", "survivingDe_Neg_4cm", HistType::kTH2F, {axisMom, axisDeltaR}); + registryDA.add("De_Neg_Q_8cm", "De_Neg_Q_8cm", HistType::kTH2F, {axisMom, axisCharge}); + registryDA.add("De_Neg_Q_4cm", "De_Neg_Q_4cm", HistType::kTH2F, {axisMom, axisCharge}); registryDA.add("De_Neg_momentum", "De_Neg_momentum", HistType::kTH2F, {{100, 0.0, 5.0, "#it{p}_{vtx} (GeV/#it{c})"}, {100, 0.0, 5.0, "#it{p}_{hmpid} (GeV/#it{c})"}}); - registryDA.add("nSigmaTPC_vs_momHMPID_noCut_DeNeg", "nSigmaTPC_vs_momHMPID_noCut_DeNeg", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5.0, 5.0, "n#sigma_TPC"}}); - registryDA.add("nSigmaTOF_vs_momHMPID_noCut_DeNeg", "nSigmaTOF_vs_momHMPID_noCut_DeNeg", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5., 5.0, "n#sigma_TOF"}}); - registryDA.add("nSigmaTPC_vs_momHMPID_Cut_DeNeg", "nSigmaTPC_vs_momHMPID_Cut_DeNeg", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5., 5.0, "n#sigma_TPC"}}); - registryDA.add("nSigmaTOF_vs_momHMPID_Cut_DeNeg", "nSigmaTOF_vs_momHMPID_Cut_DeNeg", HistType::kTH2F, {{490, 0.1, 5.0, "#it{p} (GeV/#it{c})"}, {20, -5., 5.0, "n#sigma_TOF"}}); + // first step of the analysis + registryDA.add("residualsX_Vs_momentum", "residualsX_Vs_momentum", HistType::kTH2F, {axisMom, {400, -20.0f, 20.0f, "x_{mip} - x_{track} (cm)"}}); + registryDA.add("residualsY_Vs_momentum", "residualsY_Vs_momentum", HistType::kTH2F, {axisMom, {400, -20.0f, 20.0f, "y_{mip} - y_{track} (cm)"}}); + + // nsigma plots + registryDA.add("nSigmaTPC_vs_momHMPID_noCut_DeNeg", "nSigmaTPC_vs_momHMPID_noCut_DeNeg", HistType::kTH2F, {axisMom, axisNsigmaTPC}); + registryDA.add("nSigmaTOF_vs_momHMPID_noCut_DeNeg", "nSigmaTOF_vs_momHMPID_noCut_DeNeg", HistType::kTH2F, {axisMom, axisNsigmaTOF}); + registryDA.add("nSigmaTPC_vs_momHMPID_Cut_DeNeg", "nSigmaTPC_vs_momHMPID_Cut_DeNeg", HistType::kTH2F, {axisMom, axisNsigmaTPC}); + registryDA.add("nSigmaTOF_vs_momHMPID_Cut_DeNeg", "nSigmaTOF_vs_momHMPID_Cut_DeNeg", HistType::kTH2F, {axisMom, axisNsigmaTOF}); registryDA.add("hmpidCkovvsMom", "hmpidCkovvsMom", kTH2F, {{500, 0, 10., "#it{p} (GeV/#it{c})"}, {800, 0., 0.8, "#theta_{Ch} (rad)"}}); + + // general plots + registryDA.add("hEta", "hEta", kTH1F, {axisEta}); + registryDA.add("hPhi", "hPhi", kTH1F, {axisPhi}); + registryDA.add("hMomentumHmpid", "hMomentumHmpid", kTH1F, {{1000, -10.f, 10.f, "#it{p}_{hmpid} (GeV/#it{c})"}}); + + registryDA.add("hMass", "hMass", kTH1F, {axisMass}); + registryDA.add("hMass_postDeuteron", "hMass_postDeuteron", kTH1F, {axisMass}); + + // quality check for step distributions + registryDA.add("De_Pos_deltaR_precut", "De_Pos_deltaR_precut", HistType::kTH1F, {axisDeltaR}); + registryDA.add("De_Neg_deltaR_precut", "De_Neg_deltaR_precut", HistType::kTH1F, {axisDeltaR}); + + for (int iCh = 0; iCh < nCh; iCh++) { + registryDA.add(Form("De_Pos_deltaR_%d", iCh), Form("De_Pos_deltaR_%d", iCh), HistType::kTH1F, {axisDeltaR}); + registryDA.add(Form("De_Neg_deltaR_%d", iCh), Form("De_Neg_deltaR_%d", iCh), HistType::kTH1F, {axisDeltaR}); + registryDA.add(Form("hEta_%d", iCh), Form("hEta_%d", iCh), HistType::kTH1F, {axisEta}); + registryDA.add(Form("hPhi_%d", iCh), Form("hPhi_%d", iCh), HistType::kTH1F, {axisPhi}); + } } void process(aod::HmpidAnalysis const& hmpidTable) { for (const auto& hmpid : hmpidTable) { - // filters on primary tracks + // ------------------------- + // track filters + // ------------------------- if (hmpid.itsNCluster() < minReqClusterITS) continue; if (hmpid.tpcNCluster() < minTPCnClsFound) @@ -103,89 +188,234 @@ struct HmpidDeuteron { continue; if (hmpid.itsChi2() > maxChi2ITS) continue; - if (TMath::Abs(hmpid.dcaXY()) > maxDCAxy) + if (std::abs(hmpid.dcaXY()) > maxDCAxy) continue; - if (TMath::Abs(hmpid.dcaZ()) > maxDCAz) + if (std::abs(hmpid.dcaZ()) > maxDCAz) continue; - // plots nsigma before cuts - if (hmpid.momentumHmpid() > 0) { - registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_noCut_DePos"), std::fabs(hmpid.momentumHmpid()), hmpid.tpcNSigmaDe()); - registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_noCut_DePos"), std::fabs(hmpid.momentumHmpid()), hmpid.tofNSigmaDe()); + // ------------------------- + // derived quantities + // ------------------------- + const float momHmpid = hmpid.momentumHmpid(); + const float momAbs = std::abs(momHmpid); + + const float dx = hmpid.xMip() - hmpid.xTrack(); + const float dy = hmpid.yMip() - hmpid.yTrack(); + const float dr = std::hypot(dx, dy); + + const float mass = + std::pow(nmean * momAbs * std::cos(hmpid.chAngle()), 2) - + std::pow(momAbs, 2); + + const int chamber = hmpid.chamber(); + const bool isPos = momHmpid > 0; + const bool isNeg = momHmpid < 0; + + // ------------------------- + // global distributions + // ------------------------- + registryDA.fill(HIST("hEta"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi"), hmpid.phiTrack()); + registryDA.fill(HIST("hMomentumHmpid"), momHmpid); + + if (hmpid.chAngle() > 0) + registryDA.fill(HIST("hMass"), mass); + + // ------------------------- + // RICH chamber maps (0–6) + // ------------------------- + if (chamber == rich0) { + registryDA.fill(HIST("hEta_0"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi_0"), hmpid.phiTrack()); + } + if (chamber == rich1) { + registryDA.fill(HIST("hEta_1"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi_1"), hmpid.phiTrack()); + } + if (chamber == rich2) { + registryDA.fill(HIST("hEta_2"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi_2"), hmpid.phiTrack()); + } + if (chamber == rich3) { + registryDA.fill(HIST("hEta_3"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi_3"), hmpid.phiTrack()); + } + if (chamber == rich4) { + registryDA.fill(HIST("hEta_4"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi_4"), hmpid.phiTrack()); + } + if (chamber == rich5) { + registryDA.fill(HIST("hEta_5"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi_5"), hmpid.phiTrack()); + } + if (chamber == rich6) { + registryDA.fill(HIST("hEta_6"), hmpid.etaTrack()); + registryDA.fill(HIST("hPhi_6"), hmpid.phiTrack()); + } + + // ------------------------- + // precut deltaR + // ------------------------- + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_precut"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_precut"), dr); + + // ------------------------- + // nsigma pre-cut + // ------------------------- + if (isPos) { + registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_noCut_DePos"), + momAbs, hmpid.tpcNSigmaDe()); + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_noCut_DePos"), + momAbs, hmpid.tofNSigmaDe()); } - if (hmpid.momentumHmpid() < 0) { - registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_noCut_DeNeg"), std::fabs(hmpid.momentumHmpid()), hmpid.tpcNSigmaDe()); - registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_noCut_DeNeg"), std::fabs(hmpid.momentumHmpid()), hmpid.tofNSigmaDe()); + if (isNeg) { + registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_noCut_DeNeg"), + momAbs, hmpid.tpcNSigmaDe()); + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_noCut_DeNeg"), + momAbs, hmpid.tofNSigmaDe()); } - // deuteron candidate cuts - if (hmpid.tpcNSigmaDe() < nsigmaTPCMin || hmpid.tpcNSigmaDe() > nsigmaTPCMax) + // ------------------------- + // deuteron candidate cuts - TPC + // ------------------------- + if (hmpid.tpcNSigmaDe() < nsigmaTPCMin || + hmpid.tpcNSigmaDe() > nsigmaTPCMax) continue; - if (hmpid.momentumHmpid() > 0) - registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_Cut_DePos"), std::fabs(hmpid.momentumHmpid()), hmpid.tpcNSigmaDe()); - else - registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_Cut_DeNeg"), std::fabs(hmpid.momentumHmpid()), hmpid.tpcNSigmaDe()); - - if (hmpid.tofNSigmaDe() < nsigmaTOFMin || hmpid.tofNSigmaDe() > nsigmaTOFMax) + // ------------------------- + // post TPC cut + // ------------------------- + if (isPos) + registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_Cut_DePos"), + momAbs, hmpid.tpcNSigmaDe()); + if (isNeg) + registryDA.fill(HIST("nSigmaTPC_vs_momHMPID_Cut_DeNeg"), + momAbs, hmpid.tpcNSigmaDe()); + + registryDA.fill(HIST("hMass_postDeuteron"), mass); + + // ------------------------- + // deuteron candidate cuts - TOF + // ------------------------- + if (hmpid.tofNSigmaDe() < nsigmaTOFMin || + hmpid.tofNSigmaDe() > nsigmaTOFMax) continue; - if (hmpid.momentumHmpid() > 0) - registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_Cut_DePos"), std::fabs(hmpid.momentumHmpid()), hmpid.tofNSigmaDe()); - else - registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_Cut_DeNeg"), std::fabs(hmpid.momentumHmpid()), hmpid.tofNSigmaDe()); + // ------------------------- + // TOF nsigma + // ------------------------- + if (isPos) { + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_Cut_DePos"), + momAbs, hmpid.tofNSigmaDe()); + registryDA.fill(HIST("De_Pos_momentum"), + hmpid.momentumTrack(), momAbs); + } - // plot changle vs p - registryDA.fill(HIST("hmpidCkovvsMom"), hmpid.momentumHmpid(), hmpid.chAngle()); + if (isNeg) { + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_Cut_DeNeg"), + momAbs, hmpid.tofNSigmaDe()); + } + // ------------------------- + // DEUTERON CANDIDATES AFTER ALL CUTS + // ------------------------- + + // ------------------------- + // Ckov angle vs momentum + // ------------------------- + registryDA.fill(HIST("hmpidCkovvsMom"), momHmpid, hmpid.chAngle()); + + // ------------------------- + // deltaR per chamber (Pos/Neg) + // ------------------------- + if (chamber == rich0) { + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_0"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_0"), dr); + } + if (chamber == rich1) { + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_1"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_1"), dr); + } + if (chamber == rich2) { + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_2"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_2"), dr); + } + if (chamber == rich3) { + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_3"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_3"), dr); + } + if (chamber == rich4) { + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_4"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_4"), dr); + } + if (chamber == rich5) { + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_5"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_5"), dr); + } + if (chamber == rich6) { + if (isPos) + registryDA.fill(HIST("De_Pos_deltaR_6"), dr); + if (isNeg) + registryDA.fill(HIST("De_Neg_deltaR_6"), dr); + } + + // ------------------------- // absorbers - const float dx = hmpid.xTrack() - hmpid.xMip(); - const float dy = hmpid.yTrack() - hmpid.yMip(); - const float dr = TMath::Sqrt(dx * dx + dy * dy); - - int abs4cm = 2, abs8cm = 4; - bool hmpidAbs8cm = true, hmpidAbs4cm = true; - - if (hmpid.momentumHmpid() > 0) { - registryDA.fill(HIST("De_Pos_momentum"), hmpid.momentumTrack(), std::fabs(hmpid.momentumHmpid())); - - if (hmpidAbs8cm && hmpid.chamber() == abs8cm) { - registryDA.fill(HIST("incomingDe_Pos_8cm"), std::fabs(hmpid.momentumHmpid())); - registryDA.fill(HIST("survivingDe_Pos_8cm"), std::fabs(hmpid.momentumHmpid()), dr); - registryDA.fill(HIST("De_Pos_Q_8cm"), std::fabs(hmpid.momentumHmpid()), hmpid.chargeMip()); - registryDA.fill(HIST("De_Pos_ClsSize_8cm"), std::fabs(hmpid.momentumHmpid()), hmpid.clusterSize()); - registryDA.fill(HIST("De_Pos_deltaR_8cm"), dr); - } - if (hmpidAbs4cm && hmpid.chamber() == abs4cm) { - registryDA.fill(HIST("incomingDe_Pos_4cm"), std::fabs(hmpid.momentumHmpid())); - registryDA.fill(HIST("survivingDe_Pos_4cm"), std::fabs(hmpid.momentumHmpid()), dr); - registryDA.fill(HIST("De_Pos_Q_4cm"), std::fabs(hmpid.momentumHmpid()), hmpid.chargeMip()); - registryDA.fill(HIST("De_Pos_ClsSize_4cm"), std::fabs(hmpid.momentumHmpid()), hmpid.clusterSize()); - registryDA.fill(HIST("De_Pos_deltaR_4cm"), dr); - } - } - - if (hmpid.momentumHmpid() < 0) { - registryDA.fill(HIST("De_Neg_momentum"), hmpid.momentumTrack(), std::fabs(hmpid.momentumHmpid())); - - if (hmpidAbs8cm && hmpid.chamber() == abs8cm) { - registryDA.fill(HIST("incomingDe_Neg_8cm"), std::fabs(hmpid.momentumHmpid())); - registryDA.fill(HIST("survivingDe_Neg_8cm"), std::fabs(hmpid.momentumHmpid()), dr); - registryDA.fill(HIST("De_Neg_Q_8cm"), std::fabs(hmpid.momentumHmpid()), hmpid.chargeMip()); - registryDA.fill(HIST("De_Neg_ClsSize_8cm"), std::fabs(hmpid.momentumHmpid()), hmpid.clusterSize()); - registryDA.fill(HIST("De_Neg_deltaR_8cm"), dr); - } - if (hmpidAbs4cm && hmpid.chamber() == abs4cm) { - registryDA.fill(HIST("incomingDe_Neg_4cm"), std::fabs(hmpid.momentumHmpid())); - registryDA.fill(HIST("survivingDe_Neg_4cm"), std::fabs(hmpid.momentumHmpid()), dr); - registryDA.fill(HIST("De_Neg_Q_4cm"), std::fabs(hmpid.momentumHmpid()), hmpid.chargeMip()); - registryDA.fill(HIST("De_Neg_ClsSize_4cm"), std::fabs(hmpid.momentumHmpid()), hmpid.clusterSize()); - registryDA.fill(HIST("De_Neg_deltaR_4cm"), dr); - } + // ------------------------- + const int abs4cm = 2; + const int abs8cm = 4; + + if (isPos && chamber == abs8cm) { + registryDA.fill(HIST("incomingDe_Pos_8cm"), momAbs); + registryDA.fill(HIST("survivingDe_Pos_8cm"), momAbs, dr); + registryDA.fill(HIST("De_Pos_Q_8cm"), momAbs, hmpid.chargeMip()); + registryDA.fill(HIST("De_Pos_deltaR_8cm"), dr); + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_DePos_8cm"), momAbs, hmpid.tofNSigmaDe()); } - } - } + + if (isNeg && chamber == abs8cm) { + registryDA.fill(HIST("incomingDe_Neg_8cm"), momAbs); + registryDA.fill(HIST("survivingDe_Neg_8cm"), momAbs, dr); + registryDA.fill(HIST("De_Neg_Q_8cm"), momAbs, hmpid.chargeMip()); + registryDA.fill(HIST("De_Neg_deltaR_8cm"), dr); + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_DeNeg_8cm"), momAbs, hmpid.tofNSigmaDe()); + + registryDA.fill(HIST("residualsX_Vs_momentum"), momAbs, dx); + registryDA.fill(HIST("residualsY_Vs_momentum"), momAbs, dy); + } + + if (isPos && chamber == abs4cm) { + registryDA.fill(HIST("incomingDe_Pos_4cm"), momAbs); + registryDA.fill(HIST("survivingDe_Pos_4cm"), momAbs, dr); + registryDA.fill(HIST("De_Pos_Q_4cm"), momAbs, hmpid.chargeMip()); + registryDA.fill(HIST("De_Pos_deltaR_4cm"), dr); + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_DePos_4cm"), momAbs, hmpid.tofNSigmaDe()); + } + + if (isNeg && chamber == abs4cm) { + registryDA.fill(HIST("incomingDe_Neg_4cm"), momAbs); + registryDA.fill(HIST("survivingDe_Neg_4cm"), momAbs, dr); + registryDA.fill(HIST("De_Neg_Q_4cm"), momAbs, hmpid.chargeMip()); + registryDA.fill(HIST("De_Neg_deltaR_4cm"), dr); + registryDA.fill(HIST("nSigmaTOF_vs_momHMPID_DeNeg_4cm"), momAbs, hmpid.tofNSigmaDe()); + } + } // end of loop over hmpidTable + } // end of process }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx index 0e7c1e3b464..6dd0c9f6ac7 100644 --- a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx +++ b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx @@ -29,8 +29,8 @@ #include #include #include -#include #include +#include #include #include @@ -42,6 +42,7 @@ #include #include #include +#include using namespace o2; using namespace o2::framework; @@ -55,9 +56,9 @@ struct HmpidTableProducer { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - const AxisSpec axisEvtCounter{1, 0, +1, ""}; + AxisSpec axisEvtCounter{1, 0, +1, ""}; - Service ccdb; + Service ccdb{}; struct : ConfigurableGroup { Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "URL of the CCDB repository"}; Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the material LUT"}; @@ -91,6 +92,44 @@ struct HmpidTableProducer { std::unordered_set mCollisionsWithHmpid; + static constexpr int Rich2 = 2, Rich4 = 4; + + // ----------------------------------------------------------------------- + // HMPID absorber geometry (hardcoded from HMPIDSimulation/Detector.cxx, + // Detector::ConstructGeometry / Detector::createAbsorber). + // The experiment is finalised and this geometry will not change, so the + // values are copied here instead of being re-derived from TGeoManager/CCDB + // at runtime. If the detector geometry code is ever revisited, these + // constants must be updated accordingly. + // + // Each absorber is a box (TGeoBBox) whose LOCAL->GLOBAL transform is built as: + // pMatrix->SetTranslation(T); + // pMatrix->RotateZ(theta); + // which yields, for a local point p: p_glob = Rz(theta) * p_loc + T. + // In particular the box CENTER in global coordinates is exactly T (the + // rotation does not affect T, since it is applied to p_loc only, not to + // the already-set translation). Only the box AXES are rotated by theta + // with respect to the global x,y axes. + // + // To test whether a global point lies inside the box we invert the + // transform: p_loc = Rz(-theta) * (p_glob - T), then compare component-wise + // against the box half-widths. + // ----------------------------------------------------------------------- + static constexpr double AbsThetaDeg = 33.5; + const double mAbsCosT = std::cos(AbsThetaDeg * TMath::DegToRad()); + const double mAbsSinT = std::sin(AbsThetaDeg * TMath::DegToRad()); + + // Rich2 absorber: trans2 = {435.5, 0., -155.}, thickness 40mm -> halfX = 2cm + static constexpr double AbsRich2CenterX = 435.5, AbsRich2CenterZ = -155.; + static constexpr double AbsRich2HalfX = 2.0; + + // Rich4 absorber: trans4 = {435., 0., 155.}, thickness 80mm -> halfX = 4cm + static constexpr double AbsRich4CenterX = 435.0, AbsRich4CenterZ = 155.; + static constexpr double AbsRich4HalfX = 4.0; + + // common to both absorbers: box is 1300x1300mm -> halfY = halfZ = 65cm + static constexpr double AbsHalfY = 65.0, AbsHalfZ = 65.0; + void init(o2::framework::InitContext&) { ccdb->setURL(ccdbConfig.ccdbUrl); @@ -219,7 +258,8 @@ struct HmpidTableProducer { const double kConvThr = 0.00001; const int kMaxIter = 100; - std::array x, p; + std::array x{}; + std::array p{}; for (int iter = 0; iter < kMaxIter; ++iter) { @@ -249,11 +289,17 @@ struct HmpidTableProducer { } int getHmpidChamber( - std::array x, - std::array p, + const std::array& xIn, + const std::array& pIn, double bz, int charge) { + std::array x{}; + std::array p{}; + + x = xIn; + p = pIn; + auto* param = o2::hmpid::Param::instance(); for (int ch = o2::hmpid::Param::kMinCh; ch <= o2::hmpid::Param::kMaxCh; ++ch) { @@ -268,21 +314,21 @@ struct HmpidTableProducer { param->norm(ch, nPc.data()); // Intersection track - radiator plane - std::array xRad, pAtRad; + std::array xRad{}, pAtRad{}; if (!intersectHelixPlane(bz, charge, x, p, pRad, nRad, xRad, pAtRad)) continue; // Intersection track - PC plane - std::array xPc, pAtPc; + std::array xPc{}, pAtPc{}; if (!intersectHelixPlane(bz, charge, xRad, pAtRad, pPc, nPc, xPc, pAtPc)) continue; - double theta, phi; + double theta = 0., phi = 0.; param->mars2LorsVec(ch, pAtRad.data(), theta, phi); - double xL, yL; + double xL = 0., yL = 0.; param->mars2Lors(ch, xPc.data(), xL, yL); // Use isInside to check Chamber intersected @@ -294,6 +340,38 @@ struct HmpidTableProducer { return -1; } + // Checks whether a point in global (MARS/ALICE) coordinates lies inside the + // absorber box in front of the given chamber (Rich2 or Rich4). See the + // geometry block above for the derivation of the transform used here. + bool isInAbsorber(double vx, double vy, double vz, int chamber) const + { + double centerX = 0., centerZ = 0., halfX = 0.; + + if (chamber == Rich2) { + centerX = AbsRich2CenterX; + centerZ = AbsRich2CenterZ; + halfX = AbsRich2HalfX; + } else if (chamber == Rich4) { + centerX = AbsRich4CenterX; + centerZ = AbsRich4CenterZ; + halfX = AbsRich4HalfX; + } else { + return false; + } + + // subtract the box center (translation is not rotated, see geometry block above) + const double dx = vx - centerX; + const double dy = vy; // centerY = 0 + const double dz = vz - centerZ; + + // rotate by -theta into the box local frame + const double lx = dx * mAbsCosT + dy * mAbsSinT; + const double ly = -dx * mAbsSinT + dy * mAbsCosT; + const double lz = dz; + + return std::abs(lx) <= halfX && std::abs(ly) <= AbsHalfY && std::abs(lz) <= AbsHalfZ; + } + void processEvent(CollisionCandidates::iterator const& col, aod::BCsWithTimestamps const&) { @@ -305,14 +383,15 @@ struct HmpidTableProducer { } PROCESS_SWITCH(HmpidTableProducer, processEvent, "Process event level", true); - template + template void runHmpidAnalysis( aod::HMPIDs const& hmpids, TTrackTable const&, CollisionCandidates const&, - aod::BCsWithTimestamps const&) + aod::BCsWithTimestamps const&, + TMcParticles const& mcParticles = 0) { - for (auto const& t : hmpids) { + for (auto const& t : hmpids) { // begin loop over hmpids const auto& globalTrack = t.template track_as(); @@ -356,27 +435,26 @@ struct HmpidTableProducer { static_cast(globalTrack.py()), static_cast(globalTrack.pz())}; - // int charge = (globalTrack.signed1Pt() > 0) ? +1 : -1; int16_t charge = globalTrack.sign(); auto prop = o2::base::Propagator::Instance(); - double bz = static_cast(prop->getNominalBz()); + auto bz = static_cast(prop->getNominalBz()); int chamberM1 = getHmpidChamber(x, p, bz, charge); histos.fill(HIST("hChamberM1"), chamberM1); - if (!isCorrupt) { + if (!isCorrupt) { // begin if(!isCorrupt) - fill M1vsM2 histos.fill(HIST("hChamberM1vsM2"), chamberM2, chamberM1); - } + } // end if(!isCorrupt) - fill M1vsM2 // --- M3: hybrid --- int chamberM3 = -1; - if (!isCorrupt) { + if (!isCorrupt) { // begin if/else - M3 assignment chamberM3 = chamberM2; } else { chamberM3 = chamberM1; - } + } // end if/else - M3 assignment // Legend: // bin 0 = clusSize > 0, chamber found (M2 ok) @@ -399,9 +477,11 @@ struct HmpidTableProducer { if (chamberM3 < 0) continue; - float hmpidPhotsCharge2[o2::aod::kDimPhotonsCharge]; - for (int i = 0; i < o2::aod::kDimPhotonsCharge; i++) + std::vector hmpidPhotsCharge2(o2::aod::kDimPhotonsCharge, 0.f); + + for (int i = 0; i < o2::aod::kDimPhotonsCharge; i++) { hmpidPhotsCharge2[i] = t.hmpidPhotsCharge()[i]; + } // fill hmpid table hmpidAnalysis( @@ -427,14 +507,31 @@ struct HmpidTableProducer { if constexpr (isMC) { if (globalTrack.has_mcParticle()) { const auto& mc = globalTrack.mcParticle(); - hmpidAnalysisMC(mc.pdgCode(), mc.vx(), mc.vy(), mc.vz(), mc.isPhysicalPrimary(), mc.getProcess()); + + bool interactionInAbsorber = false; + + if ((chamberM3 == Rich2 || chamberM3 == Rich4) && mc.has_daughters()) { + auto dIds = mc.daughtersIds(); + + for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) { + auto daughter = mcParticles.rawIteratorAt(idx); + + if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) { + interactionInAbsorber = true; + break; + } + } // end loop daughters + } // end if has_daughters + + hmpidAnalysisMC(mc.pdgCode(), mc.vx(), mc.vy(), mc.vz(), + mc.isPhysicalPrimary(), mc.getProcess(), interactionInAbsorber); } else { - hmpidAnalysisMC(-1, 0.f, 0.f, 0.f, false, -100); + hmpidAnalysisMC(-1, 0.f, 0.f, 0.f, false, -100, false); } - } + } // end if constexpr (isMC) - } // end HMPID loop - } + } // end for - loop over hmpids + } // end runHmpidAnalysis // process real data void processHmpid(aod::HMPIDs const& hmpids, @@ -452,9 +549,9 @@ struct HmpidTableProducer { TrackCandidatesMC const& tracks, CollisionCandidates const& cols, aod::BCsWithTimestamps const& bcs, - aod::McParticles const&) + aod::McParticles const& mcParticles) { - runHmpidAnalysis(hmpids, tracks, cols, bcs); + runHmpidAnalysis(hmpids, tracks, cols, bcs, mcParticles); } PROCESS_SWITCH(HmpidTableProducer, processHmpidMC, "Process HMPID MC entries", true); }; diff --git a/DPG/Tasks/AOTTrack/PID/HMPID/tableHMPID.h b/DPG/Tasks/AOTTrack/PID/HMPID/tableHMPID.h index 1fcf9980763..1e2445d5b05 100644 --- a/DPG/Tasks/AOTTrack/PID/HMPID/tableHMPID.h +++ b/DPG/Tasks/AOTTrack/PID/HMPID/tableHMPID.h @@ -15,6 +15,8 @@ #include #include +#include + namespace o2::aod { @@ -22,6 +24,7 @@ inline constexpr int kDimPhotonsCharge = 10; namespace variables_table { + DECLARE_SOA_COLUMN(ChAngle, chAngle, float); DECLARE_SOA_COLUMN(MomentumHmpid, momentumHmpid, float); DECLARE_SOA_COLUMN(MomentumTrack, momentumTrack, float); @@ -33,7 +36,7 @@ DECLARE_SOA_COLUMN(NPhotons, nPhotons, float); DECLARE_SOA_COLUMN(ChargeMip, chargeMip, float); DECLARE_SOA_COLUMN(ClusterSize, clusterSize, float); DECLARE_SOA_COLUMN(Chamber, chamber, float); -DECLARE_SOA_COLUMN(PhotonsCharge, photonsCharge, float[kDimPhotonsCharge]); +DECLARE_SOA_COLUMN(PhotonsCharge, photonsCharge, std::vector); DECLARE_SOA_COLUMN(EtaTrack, etaTrack, float); DECLARE_SOA_COLUMN(PhiTrack, phiTrack, float); DECLARE_SOA_COLUMN(Px, px, float); @@ -103,6 +106,7 @@ DECLARE_SOA_COLUMN(McVy, mcVy, float); DECLARE_SOA_COLUMN(McVz, mcVz, float); DECLARE_SOA_COLUMN(IsPhysPrimary, isPhysPrimary, bool); DECLARE_SOA_COLUMN(ProcessCode, processCode, int); +DECLARE_SOA_COLUMN(HasInteractedInAbsorber, hasInteractedInAbsorber, bool); } // namespace hmpid_mc @@ -112,7 +116,8 @@ DECLARE_SOA_TABLE(HmpidAnalysisMC, "AOD", "HMPIDANALYSISMC", hmpid_mc::McVy, hmpid_mc::McVz, hmpid_mc::IsPhysPrimary, - hmpid_mc::ProcessCode); + hmpid_mc::ProcessCode, + hmpid_mc::HasInteractedInAbsorber); } // namespace o2::aod From 1f9265392d8cd5bc56ace0308c7d84876e5781ab Mon Sep 17 00:00:00 2001 From: Rocco Liotino Date: Wed, 22 Jul 2026 18:29:25 +0200 Subject: [PATCH 2/2] Upgrading deuteron tasks --- .../AOTTrack/PID/HMPID/hmpidTableProducer.cxx | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx index 6dd0c9f6ac7..06c8f92b139 100644 --- a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx +++ b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx @@ -157,6 +157,8 @@ struct HmpidTableProducer { histos.add("hChamberAssignment", "Chamber assignment outcome; category; counts", kTH1F, {{4, -0.5, 3.5, ""}}); + + histos.add("hProdVertex", ";X (cm);Y (cm);Z (cm)", HistType::kTH3F, {{500, -500., 500.}, {500, -500., 500.}, {500, -500., 500.}}); } // ----------------------------------------------------------------------- @@ -360,14 +362,14 @@ struct HmpidTableProducer { } // subtract the box center (translation is not rotated, see geometry block above) - const double dx = vx - centerX; - const double dy = vy; // centerY = 0 - const double dz = vz - centerZ; + const double rx = vx * mAbsCosT + vy * mAbsSinT; + const double ry = -vx * mAbsSinT + vy * mAbsCosT; + const double rz = vz; // rotate by -theta into the box local frame - const double lx = dx * mAbsCosT + dy * mAbsSinT; - const double ly = -dx * mAbsSinT + dy * mAbsCosT; - const double lz = dz; + const double lx = rx - centerX; + const double ly = ry; // centerY = 0 + const double lz = rz - centerZ; return std::abs(lx) <= halfX && std::abs(ly) <= AbsHalfY && std::abs(lz) <= AbsHalfZ; } @@ -516,6 +518,8 @@ struct HmpidTableProducer { for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) { auto daughter = mcParticles.rawIteratorAt(idx); + histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz()); + if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) { interactionInAbsorber = true; break;