diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index 673b57b502b..23ae73b811a 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -23,6 +23,7 @@ #include #include #include +#include #include "Math/Vector4D.h" #include "Math/Vector3D.h" @@ -90,7 +91,8 @@ struct TaskPi0FlowEMC { Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable cfgDoRotation{"cfgDoRotation", true, "Flag to enable rotation background method"}; Configurable cfgDownsampling{"cfgDownsampling", 1, "Calculate rotation background only for every collision"}; - Configurable cfgRotAngle{"cfgRotAngle", M_PI / 2., "Angle used for the rotation method"}; + Configurable cfgEMCalMapLevel{"cfgEMCalMapLevel", 2, "Different levels of correction for the rotation background, the smaller number includes the level of the higher number (4: none, 3: only inside EMCal, 2: exclude bad channels, 1: remove edges)"}; + Configurable cfgRotAngle{"cfgRotAngle", std::move(const_cast(o2::constants::math::PIHalf)), "Angle used for the rotation method"}; Configurable cfgDistanceToEdge{"cfgDistanceToEdge", 1, "Distance to edge in cells required for rotated cluster to be accepted"}; // configurable axis @@ -171,6 +173,8 @@ struct TaskPi0FlowEMC { SliceCache cache; EventPlaneHelper epHelper; o2::framework::Service ccdb; + int runNow = 0; + int runBefore = -1; Filter clusterFilter = aod::skimmedcluster::time >= emccuts.cfgEMCminTime && aod::skimmedcluster::time <= emccuts.cfgEMCmaxTime && aod::skimmedcluster::m02 >= emccuts.cfgEMCminM02 && aod::skimmedcluster::m02 <= emccuts.cfgEMCmaxM02 && skimmedcluster::e >= emccuts.cfgEMCminE; Filter collisionFilter = aod::evsel::sel8 && nabs(aod::collision::posZ) <= eventcuts.cfgZvtxMax && aod::evsel::trackOccupancyInTimeRange <= eventcuts.cfgTrackOccupancyMax && aod::evsel::trackOccupancyInTimeRange >= eventcuts.cfgTrackOccupancyMin && aod::evsel::ft0cOccupancyInTimeRange <= eventcuts.cfgFT0COccupancyMax && aod::evsel::ft0cOccupancyInTimeRange >= eventcuts.cfgFT0COccupancyMin; @@ -185,9 +189,38 @@ struct TaskPi0FlowEMC { HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; o2::emcal::Geometry* emcalGeom; + o2::emcal::BadChannelMap* mBadChannels; TH1D* h1SPResolution = nullptr; + // Constants for eta and phi ranges + double etaMin = -0.75, etaMax = 0.75; + int nBinsEta = 150; // 150 bins for eta + + double phiMin = 1.35, phiMax = 5.75; + int nBinsPhi = 440; // (440 bins = 0.01 step size covering most regions) + + std::vector lookupTable1D; float epsilon = 1.e-8; + // To access the 1D array + inline int getIndex(int iEta, int iPhi) + { + return iEta * nBinsPhi + iPhi; + } + + // Function to access the lookup table + inline int8_t checkEtaPhi1D(double eta, double phi) + { + if (eta < etaMin || eta > etaMax || phi < phiMin || phi > phiMax) { + return 3; // Out of bounds + } + + // Compute indices directly + int iEta = static_cast((eta - etaMin) / ((etaMax - etaMin) / nBinsEta)); + int iPhi = static_cast((phi - phiMin) / ((phiMax - phiMin) / nBinsPhi)); + + return lookupTable1D[getIndex(iEta, iPhi)]; + } + void defineEMEventCut() { fEMEventCut = EMPhotonEventCut("fEMEventCut", "fEMEventCut"); @@ -237,9 +270,6 @@ struct TaskPi0FlowEMC { fEMCCut.SetUseTM(emccuts.cfgEMCUseTM); // disables TM o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(®istry); - // Load EMCal geometry - emcalGeom = o2::emcal::Geometry::GetInstanceFromRunNumber(300000); - const AxisSpec thnAxisInvMass{thnConfigAxisInvMass, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; const AxisSpec thnAxisPt{thnConfigAxisPt, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec thnAxisCent{thnConfigAxisCent, "Centrality (%)"}; @@ -253,12 +283,13 @@ struct TaskPi0FlowEMC { const AxisSpec thAxisEnergy{1000, 0., 100., "#it{E}_{clus} (GeV)"}; const AxisSpec thAxisEnergyCalib{100, 0., 20., "#it{E}_{clus} (GeV)"}; const AxisSpec thAxisTime{1500, -600, 900, "#it{t}_{cl} (ns)"}; - const AxisSpec thAxisEta{160, -0.8, 0.8, "#eta"}; - const AxisSpec thAxisPhi{72, 0, 2 * 3.14159, "phi"}; + const AxisSpec thAxisEta{320, -0.8, 0.8, "#eta"}; + const AxisSpec thAxisPhi{500, 0, 2 * 3.14159, "phi"}; const AxisSpec thAxisNCell{17664, 0.5, +17664.5, "#it{N}_{cell}"}; const AxisSpec thAxisPsi{360 / harmonic.value, -(1. / static_cast(harmonic.value)) * std::numbers::pi_v, (1. / static_cast(harmonic.value)) * std::numbers::pi_v, Form("#Psi_{%d}", harmonic.value)}; const AxisSpec thAxisCN{8, 0.5, 8.5, "#it{c}_{n}"}; const AxisSpec thAxisSN{8, 0.5, 8.5, "#it{s}_{n}"}; + const AxisSpec thAxisCPUTime{1000, 0, 10000, "#it{t} (#mus)"}; registry.add("hSparsePi0Flow", "THn for SP", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisScalarProd}); registry.add("hSparseBkgFlow", "THn for SP", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisScalarProd}); @@ -336,6 +367,12 @@ struct TaskPi0FlowEMC { registry.add("hInvMassPt", "Histo for inv pair mass vs pt", HistType::kTH2D, {thnAxisInvMass, thnAxisPt}); registry.add("hTanThetaPhi", "Histo for identification of conversion cluster", HistType::kTH2D, {thnAxisInvMass, thAxisTanThetaPhi}); registry.add("hAlphaPt", "Histo of meson asymmetry vs pT", HistType::kTH2D, {thAxisAlpha, thnAxisPt}); + registry.add("mesonQA/hClusterEtaPhiBefore", "hClusterEtaPhiBefore", HistType::kTH2D, {thAxisPhi, thAxisEta}); + registry.add("mesonQA/hClusterEtaPhiAfter", "hClusterEtaPhiAfter", HistType::kTH2D, {thAxisPhi, thAxisEta}); + if (cfgDoRotation) { + registry.add("mesonQA/hClusterBackEtaPhiBefore", "hClusterBackEtaPhiBefore", HistType::kTH2D, {thAxisPhi, thAxisEta}); + registry.add("mesonQA/hClusterBackEtaPhiAfter", "hClusterBackEtaPhiAfter", HistType::kTH2D, {thAxisPhi, thAxisEta}); + } } if (correctionConfig.doEMCalCalib) { @@ -352,14 +389,6 @@ struct TaskPi0FlowEMC { LOG(info) << "thnConfigAxisPt.value[1] = " << thnConfigAxisPt.value[1] << " thnConfigAxisPt.value.back() = " << thnConfigAxisPt.value.back(); }; // end init - template - void initCCDB(TCollision const& collision) - { - if (correctionConfig.cfgApplySPresolution.value) { - h1SPResolution = ccdb->getForTimeStamp(correctionConfig.cfgSpresoPath.value, collision.timestamp()); - } - } - /// Change radians to degree /// \param angle in radians /// \return angle in degree @@ -374,13 +403,7 @@ struct TaskPi0FlowEMC { float getDeltaPsiInRange(float psi1, float psi2) { float deltaPsi = psi1 - psi2; - if (std::abs(deltaPsi) > constants::math::PI / harmonic) { - if (deltaPsi > 0.) - deltaPsi -= constants::math::TwoPI / harmonic; - else - deltaPsi += constants::math::TwoPI / harmonic; - } - return deltaPsi; + return RecoDecay::constrainAngle(deltaPsi, 0.f, harmonic); } /// Fill THnSparse @@ -556,6 +579,53 @@ struct TaskPi0FlowEMC { return false; } + bool isCellMasked(int cellID) + { + bool masked = false; + if (mBadChannels) { + auto maskStatus = mBadChannels->getChannelStatus(cellID); + masked = (maskStatus != o2::emcal::BadChannelMap::MaskType_t::GOOD_CELL); + } + return masked; + } + + template + void initCCDB(TCollision const& collision) + { + // Load EMCal geometry + emcalGeom = o2::emcal::Geometry::GetInstanceFromRunNumber(collision.runNumber()); + // Load Bad Channel map + mBadChannels = ccdb->getForTimeStamp("EMC/Calib/BadChannelMap", collision.timestamp()); + lookupTable1D = std::vector(nBinsEta * nBinsPhi, -1); + double binWidthEta = (etaMax - etaMin) / nBinsEta; + double binWidthPhi = (phiMax - phiMin) / nBinsPhi; + + for (int iEta = 0; iEta < nBinsEta; ++iEta) { + double etaCenter = etaMin + (iEta + 0.5) * binWidthEta; + for (int iPhi = 0; iPhi < nBinsPhi; ++iPhi) { + double phiCenter = phiMin + (iPhi + 0.5) * binWidthPhi; + try { + // Get the cell ID + int cellID = emcalGeom->GetAbsCellIdFromEtaPhi(etaCenter, phiCenter); + + // Check conditions for the cell + if (isTooCloseToEdge(cellID, 1)) { + lookupTable1D[getIndex(iEta, iPhi)] = 2; // Edge + } else if (isCellMasked(cellID)) { + lookupTable1D[getIndex(iEta, iPhi)] = 1; // Bad + } else { + lookupTable1D[getIndex(iEta, iPhi)] = 0; // Good + } + } catch (o2::emcal::InvalidPositionException& e) { + lookupTable1D[getIndex(iEta, iPhi)] = 3; // Outside geometry + } + } + } + if (correctionConfig.cfgApplySPresolution.value) { + h1SPResolution = ccdb->getForTimeStamp(correctionConfig.cfgSpresoPath.value, collision.timestamp()); + } + } + /// \brief Calculate background using rotation background method template void rotationBackground(const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const& photons_coll, unsigned int ig1, unsigned int ig2, CollsWithQvecs::iterator const& collision) @@ -575,23 +645,21 @@ struct TaskPi0FlowEMC { photon1 = rotationMatrix * photon1; photon2 = rotationMatrix * photon2; - try { - iCellIDPhoton1 = emcalGeom->GetAbsCellIdFromEtaPhi(photon1.Eta(), photon1.Phi()); - if (isTooCloseToEdge(iCellIDPhoton1, cfgDistanceToEdge.value)) { - iCellIDPhoton1 = -1; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (emccuts.cfgEnableQA) { + registry.fill(HIST("mesonQA/hClusterBackEtaPhiBefore"), RecoDecay::constrainAngle(photon1.Phi()), photon1.Eta()); // before check but after rotation + registry.fill(HIST("mesonQA/hClusterBackEtaPhiBefore"), RecoDecay::constrainAngle(photon2.Phi()), photon2.Eta()); // before check but after rotation + } + + if (checkEtaPhi1D(photon1.Eta(), RecoDecay::constrainAngle(photon1.Phi())) >= cfgEMCalMapLevel.value) { iCellIDPhoton1 = -1; + } else if (emccuts.cfgEnableQA) { + registry.fill(HIST("mesonQA/hClusterBackEtaPhiAfter"), RecoDecay::constrainAngle(photon1.Phi()), photon1.Eta()); // after check } - try { - iCellIDPhoton2 = emcalGeom->GetAbsCellIdFromEtaPhi(photon2.Eta(), photon2.Phi()); - if (isTooCloseToEdge(iCellIDPhoton2, cfgDistanceToEdge.value)) { - iCellIDPhoton2 = -1; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (checkEtaPhi1D(photon2.Eta(), RecoDecay::constrainAngle(photon2.Phi())) >= cfgEMCalMapLevel.value) { iCellIDPhoton2 = -1; + } else if (emccuts.cfgEnableQA) { + registry.fill(HIST("mesonQA/hClusterBackEtaPhiAfter"), RecoDecay::constrainAngle(photon2.Phi()), photon2.Eta()); // after check } - if (iCellIDPhoton1 == -1 && iCellIDPhoton2 == -1) { return; } @@ -603,8 +671,11 @@ struct TaskPi0FlowEMC { if (!(fEMCCut.IsSelected(photon))) { continue; } + if (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevel.value) { + continue; + } ROOT::Math::PtEtaPhiMVector photon3(photon.pt(), photon.eta(), photon.phi(), 0.); - if (iCellIDPhoton1 > 0) { + if (iCellIDPhoton1 >= 0) { ROOT::Math::PtEtaPhiMVector mother1 = photon1 + photon3; float openingAngle1 = std::acos(photon1.Vect().Dot(photon3.Vect()) / (photon1.P() * photon3.P())); float cosNPhi1 = std::cos(harmonic * mother1.Phi()); @@ -627,7 +698,7 @@ struct TaskPi0FlowEMC { } } } - if (iCellIDPhoton2 > 0) { + if (iCellIDPhoton2 >= 0) { ROOT::Math::PtEtaPhiMVector mother2 = photon2 + photon3; float openingAngle2 = std::acos(photon2.Vect().Dot(photon3.Vect()) / (photon2.P() * photon3.P())); float cosNPhi2 = std::cos(harmonic * mother2.Phi()); @@ -671,20 +742,10 @@ struct TaskPi0FlowEMC { photon1 = rotationMatrix * photon1; photon2 = rotationMatrix * photon2; - try { - iCellIDPhoton1 = emcalGeom->GetAbsCellIdFromEtaPhi(photon1.Eta(), photon1.Phi()); - if (isTooCloseToEdge(iCellIDPhoton1, cfgDistanceToEdge.value)) { - iCellIDPhoton1 = -1; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (checkEtaPhi1D(photon1.Eta(), RecoDecay::constrainAngle(photon1.Phi())) >= cfgEMCalMapLevel.value) { iCellIDPhoton1 = -1; } - try { - iCellIDPhoton2 = emcalGeom->GetAbsCellIdFromEtaPhi(photon2.Eta(), photon2.Phi()); - if (isTooCloseToEdge(iCellIDPhoton2, cfgDistanceToEdge.value)) { - iCellIDPhoton2 = -1; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (checkEtaPhi1D(photon2.Eta(), RecoDecay::constrainAngle(photon2.Phi())) >= cfgEMCalMapLevel.value) { iCellIDPhoton2 = -1; } @@ -699,9 +760,12 @@ struct TaskPi0FlowEMC { if (!(fEMCCut.IsSelected(photon))) { continue; } + if (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevel.value) { + continue; + } ROOT::Math::PtEtaPhiMVector photon3(photon.pt(), photon.eta(), photon.phi(), 0.); - if (iCellIDPhoton1 > 0) { - if ((photon1.E() - photon3.E()) / (photon1.E() + photon3.E()) > 0.9) { // only use symmetric decays + if (iCellIDPhoton1 >= 0) { + if (std::fabs((photon1.E() - photon3.E()) / (photon1.E() + photon3.E()) < 0.1)) { // only use symmetric decays ROOT::Math::PtEtaPhiMVector mother1 = photon1 + photon3; float openingAngle1 = std::acos(photon1.Vect().Dot(photon3.Vect()) / (photon1.P() * photon3.P())); @@ -718,8 +782,8 @@ struct TaskPi0FlowEMC { } } } - if (iCellIDPhoton2 > 0) { - if ((photon2.E() - photon3.E()) / (photon2.E() + photon3.E()) > 0.9) { // only use symmetric decays + if (iCellIDPhoton2 >= 0) { + if (std::fabs((photon2.E() - photon3.E()) / (photon2.E() + photon3.E()) < 0.1)) { // only use symmetric decays ROOT::Math::PtEtaPhiMVector mother2 = photon2 + photon3; float openingAngle2 = std::acos(photon2.Vect().Dot(photon3.Vect()) / (photon2.P() * photon3.P())); @@ -815,7 +879,11 @@ struct TaskPi0FlowEMC { // selection based on QVector continue; } - initCCDB(collision); + runNow = collision.runNumber(); + if (runNow != runBefore) { + initCCDB(collision); + runBefore = runNow; + } o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted @@ -823,10 +891,15 @@ struct TaskPi0FlowEMC { if (emccuts.cfgEnableQA) { for (const auto& photon : photonsPerCollision) { registry.fill(HIST("hEClusterBefore"), photon.e()); // before cuts + registry.fill(HIST("mesonQA/hClusterEtaPhiBefore"), photon.phi(), photon.eta()); // before cuts if (!(fEMCCut.IsSelected(photon))) { continue; } + if (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= 2) { + continue; + } registry.fill(HIST("hEClusterAfter"), photon.e()); // accepted after cuts + registry.fill(HIST("mesonQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // before cuts } } for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { @@ -836,22 +909,10 @@ struct TaskPi0FlowEMC { // Cut edge clusters away, similar to rotation method to ensure same acceptance is used if (cfgDistanceToEdge.value) { - int iCellIDPhoton1 = -1; - int iCellIDPhoton2 = -1; - try { - iCellIDPhoton1 = emcalGeom->GetAbsCellIdFromEtaPhi(g1.eta(), g1.phi()); - if (isTooCloseToEdge(iCellIDPhoton1, cfgDistanceToEdge.value)) { - continue; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= 2) { continue; } - try { - iCellIDPhoton2 = emcalGeom->GetAbsCellIdFromEtaPhi(g2.eta(), g2.phi()); - if (isTooCloseToEdge(iCellIDPhoton2, cfgDistanceToEdge.value)) { - continue; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (checkEtaPhi1D(g2.eta(), RecoDecay::constrainAngle(g2.phi())) >= 2) { continue; } } @@ -935,7 +996,11 @@ struct TaskPi0FlowEMC { // selection based on QVector continue; } - initCCDB(c1); + runNow = c1.runNumber(); + if (runNow != runBefore) { + initCCDB(c1); + runBefore = runNow; + } for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(clusters1, clusters2))) { if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { continue; @@ -1124,7 +1189,11 @@ struct TaskPi0FlowEMC { // selection based on QVector continue; } - initCCDB(collision); + runNow = collision.runNumber(); + if (runNow != runBefore) { + initCCDB(collision); + runBefore = runNow; + } o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted @@ -1145,22 +1214,10 @@ struct TaskPi0FlowEMC { // Cut edge clusters away, similar to rotation method to ensure same acceptance is used if (cfgDistanceToEdge.value) { - int iCellIDPhoton1 = -1; - int iCellIDPhoton2 = -1; - try { - iCellIDPhoton1 = emcalGeom->GetAbsCellIdFromEtaPhi(g1.eta(), g1.phi()); - if (isTooCloseToEdge(iCellIDPhoton1, cfgDistanceToEdge.value)) { - continue; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= 2) { continue; } - try { - iCellIDPhoton2 = emcalGeom->GetAbsCellIdFromEtaPhi(g2.eta(), g2.phi()); - if (isTooCloseToEdge(iCellIDPhoton2, cfgDistanceToEdge.value)) { - continue; - } - } catch (o2::emcal::InvalidPositionException& e) { + if (checkEtaPhi1D(g2.eta(), RecoDecay::constrainAngle(g2.phi())) >= 2) { continue; } } @@ -1198,7 +1255,7 @@ struct TaskPi0FlowEMC { registry.fill(HIST("hClusterCuts"), 5); continue; } - if ((v1.E() - v2.E()) / (v1.E() + v2.E()) > 0.9) { // only use symmetric decays + if (std::fabs((v1.E() - v2.E()) / (v1.E() + v2.E()) < 0.1)) { // only use symmetric decays registry.fill(HIST("hClusterCuts"), 6); registry.fill(HIST("hSparseCalibSE"), vMeson.M(), vMeson.E() / 2., getCentrality(collision)); }