diff --git a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx index 06c8f92b139..1074a071b60 100644 --- a/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx +++ b/DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx @@ -34,6 +34,8 @@ #include #include +#include +#include #include @@ -73,6 +75,12 @@ struct HmpidTableProducer { Configurable requireTPC{"requireTPC", true, "Require TPC track"}; Configurable requireTOF{"requireTOF", true, "Require TOF track"}; + Configurable useInAbsorberGeomMethod{"useInAbsorberGeomMethod", false, "Use geometrical method to check if daughters are born in absorber"}; + + // (reference) 473 cm - was the legacy value in run2 simulation + Configurable survivalThresholdRich2{"survivalThresholdRich2", 437.5f, "survivalThresholdRich2"}; + Configurable survivalThresholdRich4{"survivalThresholdRich4", 439.0f, "survivalThresholdRich4"}; + using CollisionCandidates = o2::soa::Join; using TrackCandidates = soa::JoinGLOBAL 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. - // ----------------------------------------------------------------------- + // (reference) HMPID Detector class in O2 static constexpr double AbsThetaDeg = 33.5; - const double mAbsCosT = std::cos(AbsThetaDeg * TMath::DegToRad()); - const double mAbsSinT = std::sin(AbsThetaDeg * TMath::DegToRad()); + double mAbsCosT = 0., mAbsSinT = 0.; // Rich2 absorber: trans2 = {435.5, 0., -155.}, thickness 40mm -> halfX = 2cm static constexpr double AbsRich2CenterX = 435.5, AbsRich2CenterZ = -155.; @@ -132,6 +119,9 @@ struct HmpidTableProducer { void init(o2::framework::InitContext&) { + mAbsCosT = std::cos(AbsThetaDeg * TMath::DegToRad()); + mAbsSinT = std::sin(AbsThetaDeg * TMath::DegToRad()); + ccdb->setURL(ccdbConfig.ccdbUrl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -159,6 +149,10 @@ struct HmpidTableProducer { 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.}}); + histos.add("hDaughterRCyl_Rich2", "hDaughterRCyl_Rich2", kTH1F, {{600, 0., 600.}}); + histos.add("hDaughterRCyl_Rich4", "hDaughterRCyl_Rich4", kTH1F, {{600, 0., 600.}}); + histos.add("hDaughterRSph_Rich2", "hDaughterRSph_Rich2", kTH1F, {{600, 0., 600.}}); + histos.add("hDaughterRSph_Rich4", "hDaughterRSph_Rich4", kTH1F, {{600, 0., 600.}}); } // ----------------------------------------------------------------------- @@ -274,8 +268,9 @@ struct HmpidTableProducer { (x[1] - planePoint[1]) * planeNormal[1] + (x[2] - planePoint[2]) * planeNormal[2]; - if (std::abs(dist) >= std::abs(distPrev)) + if (std::abs(dist) >= std::abs(distPrev)) { return false; + } distPrev = dist; s -= dist; @@ -318,14 +313,16 @@ struct HmpidTableProducer { // Intersection track - radiator plane std::array xRad{}, pAtRad{}; - if (!intersectHelixPlane(bz, charge, x, p, pRad, nRad, xRad, pAtRad)) + if (!intersectHelixPlane(bz, charge, x, p, pRad, nRad, xRad, pAtRad)) { continue; + } // Intersection track - PC plane std::array xPc{}, pAtPc{}; - if (!intersectHelixPlane(bz, charge, xRad, pAtRad, pPc, nPc, xPc, pAtPc)) + if (!intersectHelixPlane(bz, charge, xRad, pAtRad, pPc, nPc, xPc, pAtPc)) { continue; + } double theta = 0., phi = 0.; param->mars2LorsVec(ch, pAtRad.data(), theta, phi); @@ -334,8 +331,9 @@ struct HmpidTableProducer { param->mars2Lors(ch, xPc.data(), xL, yL); // Use isInside to check Chamber intersected - if (param->isInside(xL, yL, param->distCut())) + if (param->isInside(xL, yL, param->distCut())) { return ch; + } } // No chamber intersected @@ -361,17 +359,32 @@ struct HmpidTableProducer { return false; } - // subtract the box center (translation is not rotated, see geometry block above) - const double rx = vx * mAbsCosT + vy * mAbsSinT; - const double ry = -vx * mAbsSinT + vy * mAbsCosT; - const double rz = vz; + // translate to box center + const double lx = vx - centerX; + const double ly = vy; // centerY = 0 + const double lz = vz - centerZ; // rotate by -theta into the box local frame - const double lx = rx - centerX; - const double ly = ry; // centerY = 0 - const double lz = rz - centerZ; + const double rx = lx * mAbsCosT + ly * mAbsSinT; + const double ry = -lx * mAbsSinT + ly * mAbsCosT; + const double rz = lz; + + return std::abs(rx) <= halfX && std::abs(ry) <= AbsHalfY && std::abs(rz) <= AbsHalfZ; + } + + bool survivedAbsorber(double vx, double vy, int chamber) + { + float thresholdR = 0.; + if (chamber == Rich2) { + thresholdR = survivalThresholdRich2; + } else if (chamber == Rich4) { + thresholdR = survivalThresholdRich4; + } else { + return false; + } - return std::abs(lx) <= halfX && std::abs(ly) <= AbsHalfY && std::abs(lz) <= AbsHalfZ; + const float r = std::hypot(vx, vy); + return r > thresholdR; } void processEvent(CollisionCandidates::iterator const& col, @@ -397,8 +410,9 @@ struct HmpidTableProducer { const auto& globalTrack = t.template track_as(); - if (!globalTrack.has_collision()) + if (!globalTrack.has_collision()) { continue; + } const auto& col = globalTrack.template collision_as(); initCCDB(col.template bc_as()); @@ -406,17 +420,20 @@ struct HmpidTableProducer { if ((requireITS && !globalTrack.hasITS()) || (requireTPC && !globalTrack.hasTPC()) || - (requireTOF && !globalTrack.hasTOF())) + (requireTOF && !globalTrack.hasTOF())) { continue; + } - if (mCollisionsWithHmpid.insert(collId).second) + if (mCollisionsWithHmpid.insert(collId).second) { histos.fill(HIST("eventsHmpid"), 0.5); + } // clusSize diagnostics histos.fill(HIST("hClusSize"), t.hmpidClusSize()); bool isCorrupt = (t.hmpidClusSize() <= 0); - if (isCorrupt) + if (isCorrupt) { histos.fill(HIST("hClusSizeCorrupt"), t.hmpidClusSize()); + } // --- M2: clusSize encoding --- int chamberM2 = t.hmpidClusSize() / 1000000; @@ -464,20 +481,22 @@ struct HmpidTableProducer { // bin 2 = clusSize <= 0, M1 recovery (corrupt, M1 ok) // bin 3 = clusSize <= 0, M1 fails (corrupt, skipped) - if (!isCorrupt && chamberM3 >= 0) + if (!isCorrupt && chamberM3 >= 0) { histos.fill(HIST("hChamberAssignment"), 0.); - else if (!isCorrupt && chamberM3 < 0) + } else if (!isCorrupt && chamberM3 < 0) { histos.fill(HIST("hChamberAssignment"), 1.); - else if (isCorrupt && chamberM3 >= 0) + } else if (isCorrupt && chamberM3 >= 0) { histos.fill(HIST("hChamberAssignment"), 2.); - else + } else { histos.fill(HIST("hChamberAssignment"), 3.); + } histos.fill(HIST("hChamberM3"), chamberM3); histos.fill(HIST("hChamberM3vsM2"), chamberM2, chamberM3); - if (chamberM3 < 0) + if (chamberM3 < 0) { continue; + } std::vector hmpidPhotsCharge2(o2::aod::kDimPhotonsCharge, 0.f); @@ -514,23 +533,78 @@ struct HmpidTableProducer { 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); - - histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz()); - - if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) { - interactionInAbsorber = true; - break; - } - } // end loop daughters + bool foundRelevantDaughter = false; // true if at least one non-delta/photon daughter was examined + + if (useInAbsorberGeomMethod) { + for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) { + auto daughter = mcParticles.rawIteratorAt(idx); + + int absPdg = std::abs(daughter.pdgCode()); + if (absPdg == kElectron || absPdg == kGamma) { + continue; + } + + foundRelevantDaughter = true; + + // diagnostics on daughters distribution + histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz()); + + double rCyl = std::hypot(daughter.vx(), daughter.vy()); + double rSph = std::hypot(daughter.vx(), daughter.vy(), daughter.vz()); + if (chamberM3 == Rich2) { + histos.fill(HIST("hDaughterRCyl_Rich2"), rCyl); + histos.fill(HIST("hDaughterRSph_Rich2"), rSph); + } else { + histos.fill(HIST("hDaughterRCyl_Rich4"), rCyl); + histos.fill(HIST("hDaughterRSph_Rich4"), rSph); + } + + if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) { + interactionInAbsorber = true; + } + } // end loop daughters + } else { + bool survived = false; + for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) { + auto daughter = mcParticles.rawIteratorAt(idx); + + // skip delta rays (e-/e+), photons, and HMPID Cherenkov/feedback + int absPdg = std::abs(daughter.pdgCode()); + if (absPdg == kElectron || absPdg == kGamma) { + continue; + } + + foundRelevantDaughter = true; + + histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz()); + + double rCyl = std::hypot(daughter.vx(), daughter.vy()); + double rSph = std::hypot(daughter.vx(), daughter.vy(), daughter.vz()); + if (chamberM3 == Rich2) { + histos.fill(HIST("hDaughterRCyl_Rich2"), rCyl); + histos.fill(HIST("hDaughterRSph_Rich2"), rSph); + } else { + histos.fill(HIST("hDaughterRCyl_Rich4"), rCyl); + histos.fill(HIST("hDaughterRSph_Rich4"), rSph); + } + + if (survivedAbsorber(daughter.vx(), daughter.vy(), chamberM3)) { + survived = true; + } + } // end loop daughters + + // No relevant daughter found (only delta rays/photons, or no + // daughters at all): no evidence of a genuine interaction -> treat + // as primary/survived, consistent with the "no daughters" case. + interactionInAbsorber = foundRelevantDaughter ? !survived : false; + } } // 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, false); + // No MC truth associated to this track + hmpidAnalysisMC(-999, -999.f, -999.f, -999.f, false, -100, false); } } // end if constexpr (isMC)