diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCandidate.h b/PWGEM/PhotonMeson/Core/V0PhotonCandidate.h index 37aeb4b5425..d653194ca4d 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCandidate.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCandidate.h @@ -89,8 +89,8 @@ struct V0PhotonCandidate { dcaXYV0ToPV = RecoDecay::sqrtSumOfSquares(dcaXV0ToPV, dcaYV0ToPV) * tmpSign; dcaZV0ToPV = (v0DecayVtx.GetZ() - v0DecayVtx.GetPz() * cospa_ * length / v0mom) - collision.posZ(); - alpha = v0_alpha(posPx, posPy, posPz, elePx, elePy, elePz); - qt = v0_qt(posPx, posPy, posPz, elePx, elePy, elePz); + alpha = v0Alpha(posPx, posPy, posPz, elePx, elePy, elePz); + qt = v0Qt(posPx, posPy, posPz, elePx, elePy, elePz); cospa = cospa_; cospaRZ = cospaRZ_; diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx index 3bb208f0452..32d3d7c3407 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx @@ -75,7 +75,7 @@ struct SkimmerGammaCalo { Configurable minE{"minE", 0.5, "Minimum energy for energy cut"}; Configurable maxE{"maxE", std::numeric_limits::max(), "Maximum energy for energy cut"}; Configurable removeExotic{"removeExotic", false, "Flag to enable the removal of exotic clusters."}; - Configurable> clusterDefinitions{"clusterDefinitions", {0, 1, 2, 10, 11, 12, 13, 20, 21, 22, 30, 40, 41, 42, 43, 44, 45}, "Cluster definitions to be accepted (e.g. 13 for kV3MostSplitLowSeed)"}; + Configurable> clusterDefinitions{"clusterDefinitions", {0, 1, 2, 10, 11, 12, 13, 20, 21, 22, 30, 40, 41, 42, 43, 44, 45}, "Cluster definitions storageID to be accepted (see O2Physics/PWGJE/DataModel/EMCALClusters.h for the storageID)"}; Configurable maxdEta{"maxdEta", 0.1, "Set a maximum difference in eta for tracks and cluster to still count as matched"}; Configurable maxdPhi{"maxdPhi", 0.1, "Set a maximum difference in phi for tracks and cluster to still count as matched"}; Configurable maxEoverP{"maxEoverP", 1.5, "Set a maximum for cluster E / track p for track matching."}; @@ -93,8 +93,8 @@ struct SkimmerGammaCalo { << "cluster table and must not be enabled together — this doubles MinClusters rows " << "relative to EMCClusterMCLabels_001."; } - historeg.add("DefinitionIn", "Cluster definitions before cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{51, -0.5, 50.5}}); - historeg.add("DefinitionOut", "Cluster definitions after cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + historeg.add("DefinitionIn", "Cluster definitions before cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{53, -0.5, 52.5}}); + historeg.add("DefinitionOut", "Cluster definitions after cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{53, -0.5, 52.5}}); historeg.add("EIn", "Energy of clusters before cuts", gHistoSpecClusterE); historeg.add("EOut", "Energy of clusters after cuts", gHistoSpecClusterE); historeg.add("MTEtaPhiBeforeTM", "Eta phi of matched tracks before TM cuts", gHistoSpecClusterTMdEtadPhi); diff --git a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt index 996e8f14038..dde484f2daf 100644 --- a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt +++ b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt @@ -126,6 +126,11 @@ o2physics_add_dpl_workflow(check-mc-v0 PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(task-photon-flow + SOURCES taskPhotonFlow.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(pi0-flow-emc SOURCES taskPi0FlowEMC.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore diff --git a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx index 17c5e908189..496cf7a8067 100644 --- a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx @@ -43,6 +43,7 @@ #include "PWGEM/PhotonMeson/Utils/MCUtilities.h" #include "Common/CCDB/EventSelectionParams.h" +#include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" @@ -50,18 +51,17 @@ #include #include #include -#include #include #include #include #include +#include #include #include #include #include #include #include -#include #include #include @@ -69,6 +69,7 @@ #include #include +#include // for the concept std::ranges::contiguous_range #include #include #include @@ -84,7 +85,7 @@ using namespace o2::aod::pwgem::photon; using namespace o2::aod::pwgem::photonmeson::utils::mcutil; using namespace o2::aod::pwgem::dilepton::utils::mcutil; -using MyCollisions = soa::Join; +using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; using MyV0Photons = soa::Join; @@ -94,7 +95,7 @@ using MyV0PhotonsML = soa::Join; using MyV0PhotonML = MyV0PhotonsML::iterator; // MC Joins -using MyCollisionsMC = soa::Join; +using MyCollisionsMC = soa::Join; using MyCollisionMC = MyCollisionsMC::iterator; using MyMCCollisions = soa::Join; @@ -111,9 +112,6 @@ struct PCMQC { Configurable cfgCentMin{"cfgCentMin", 0, "min. centrality"}; Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; Configurable d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"}; EMPhotonEventCut fEMEventCut; @@ -180,6 +178,7 @@ struct PCMQC { struct : ConfigurableGroup { std::string prefix = "qaSettings_group"; Configurable cfgDoPtDependentLossQA{"cfgDoPtDependentLossQA", false, "fill the cut variables vs. pT before AND after the V0 selection - the after/before ratio shows which candidates the cuts remove"}; + Configurable cfgDoKappaAnalysis{"cfgDoKappaAnalysis", false, "fill histograms with kappa vs pT vs cent."}; } qaSettingsGroup; // PCM ML inference @@ -258,45 +257,28 @@ struct PCMQC { ccdb->setFatalWhenNull(false); } - template + template void initCCDB(TCollision const& collision) { if (mRunNumber == collision.runNumber()) { return; } + mRunNumber = collision.runNumber(); // In case override, don't proceed, please - no CCDB access required - if (d_bz_input > -990) { // o2-linter: disable=magic-number (dummy value to indicate override) + if (d_bz_input > -990) { // o2-linter: disable=magic-number (override value) d_bz = d_bz_input; o2::parameters::GRPMagField grpmag; - if (std::fabs(d_bz) > 1e-5) { // o2-linter: disable=magic-number (dummy value to indicate override) - grpmag.setL3Current(30000.f / (d_bz / 5.0f)); // o2-linter: disable=magic-number (dummy value to indicate override) + if (std::fabs(d_bz) > 1e-5) { // o2-linter: disable=magic-number (override value) + grpmag.setL3Current(30000.f / (d_bz / 5.0f)); // o2-linter: disable=magic-number (override value) } - mRunNumber = collision.runNumber(); return; } - auto run3grp_timestamp = collision.timestamp(); - o2::parameters::GRPObject* grpo = nullptr; - o2::parameters::GRPMagField* grpmag = nullptr; - if (!skipGRPOquery) { - grpo = ccdb->getForTimeStamp(grpPath, run3grp_timestamp); - } - if (grpo) { - // Fetch magnetic field from ccdb for current collision - d_bz = grpo->getNominalL3Field(); - LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; - } else { - grpmag = ccdb->getForTimeStamp(grpmagPath, run3grp_timestamp); - if (!grpmag) { - LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp; - } - // Fetch magnetic field from ccdb for current collision - d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); - LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; - } + // Fetch magnetic field from ccdb for current collision + d_bz = collision.grpMagField().getNominalL3Field(); + LOG(info) << "Retrieved GRP for timestamp " << collision.timestamp() << " with magnetic field of " << d_bz << " kZG"; fV0PhotonCut.SetD_Bz(d_bz); - mRunNumber = collision.runNumber(); } void addhistograms() @@ -573,6 +555,19 @@ struct PCMQC { fRegistry.addClone("RecoQA/LegQuality/survived/", "RecoQA/LegQuality/lost/"); } } + + if (qaSettingsGroup.cfgDoKappaAnalysis.value) { + if (doprocessQC || doprocessQCML) { + fRegistry.add("Kappa/Rec", ";#Kappa;#it{p}_{T} (GeV/#it{c});cent (%)", kTH3D, {{200, -10.f, 10.f}, {100, 0., 10.}, {20, 0., 100.}}, true); + } + if (doprocessPCMQCMC || doprocessPCMQCMCML) { + fRegistry.add("Kappa/MC/v0", "v0 ee pairs;#Kappa;#it{p}_{T} (GeV/#it{c});cent (%)", kTH3D, {{200, -10.f, 10.f}, {100, 0., 10.}, {20, 0., 100.}}, true); + fRegistry.add("Kappa/MC/ee", "ee pairs not from v0;#Kappa;#it{p}_{T} (GeV/#it{c});cent (%)", kTH3D, {{200, -10.f, 10.f}, {100, 0., 10.}, {20, 0., 100.}}, true); + fRegistry.add("Kappa/MC/pie", "pion + e pairs;#Kappa;#it{p}_{T} (GeV/#it{c});cent (%)", kTH3D, {{200, -10.f, 10.f}, {100, 0., 10.}, {20, 0., 100.}}, true); + fRegistry.add("Kappa/MC/pipi", "pion pion pairs;#Kappa;#it{p}_{T} (GeV/#it{c});cent (%)", kTH3D, {{200, -10.f, 10.f}, {100, 0., 10.}, {20, 0., 100.}}, true); + fRegistry.add("Kappa/MC/other", "other pairs;#Kappa;#it{p}_{T} (GeV/#it{c});cent (%)", kTH3D, {{200, -10.f, 10.f}, {100, 0., 10.}, {20, 0., 100.}}, true); + } + } } void DefineEMEventCut() @@ -648,7 +643,7 @@ struct PCMQC { } } - template + template void fillEventInfo(TCollision const& collision, const float /*weight*/ = 1.f) { fRegistry.fill(HIST("Event/") + HIST(event_types[ev_id]) + HIST("hCollisionCounter"), 1.0); @@ -688,7 +683,7 @@ struct PCMQC { fRegistry.fill(HIST("Event/") + HIST(event_types[ev_id]) + HIST("hMultFT0MvsMultNTracksPV"), collision.multFT0A() + collision.multFT0C(), collision.multNTracksPV()); } - template + template void fillV0Info(TV0 const& v0) { fRegistry.fill(HIST("V0/hPt"), v0.pt()); @@ -717,7 +712,7 @@ struct PCMQC { fRegistry.fill(HIST("V0/hKFChi2vsZ"), v0.vz(), v0.chiSquareNDF()); float phi_cp = std::atan2(v0.vy(), v0.vx()); - o2::math_utils::bringTo02Pi(phi_cp); + RecoDecay::constrainAngle(phi_cp); float eta_cp = std::atanh(v0.vz() / std::sqrt(std::pow(v0.vx(), 2) + std::pow(v0.vy(), 2) + std::pow(v0.vz(), 2))); fRegistry.fill(HIST("V0/hsConvPoint"), v0.v0radius(), phi_cp, eta_cp); @@ -744,7 +739,7 @@ struct PCMQC { } } - template + template void fillV0LegInfo(TLeg const& leg) { fRegistry.fill(HIST("V0Leg/hPt"), leg.pt()); @@ -769,7 +764,7 @@ struct PCMQC { fRegistry.fill(HIST("V0Leg/hTPCNsigmaElVsEta"), leg.eta(), leg.tpcNSigmaEl()); } - template + template void fillLossQAInfo(TV0 const& v0) { if (!qaSettingsGroup.cfgDoPtDependentLossQA) { @@ -782,7 +777,7 @@ struct PCMQC { fRegistry.fill(HIST("V0/LossQA/") + HIST(event_types[ev_id]) + HIST("hMeeVsPt"), v0.pt(), v0.mGamma()); } - template + template void fillMaterialBudgetInfo(TV0 const& v0) { if (materialBudgetSettingsGroup.cfgDoMaterialDistribution) { @@ -806,16 +801,57 @@ struct PCMQC { } } - Preslice perCollisionV0 = aod::v0photonkf::pmeventId; - Preslice perCollisionV0ML = aod::v0photonkf::pmeventId; + //_______________________________________________________________________ + template + void fillKappaRec(TLeg const& pos, TLeg const& ele, const float pt, const float cent) + { + if (!qaSettingsGroup.cfgDoKappaAnalysis.value) { + return; + } + const float kappa = getV0Kappa(pos, ele); + fRegistry.fill(HIST("Kappa/Rec"), kappa, pt, cent); + } + + //_______________________________________________________________________ + template + void fillKappaMC(TLeg const& pos, TLeg const& ele, TMCParticle const& posmc, TMCParticle const& elemc, const int photonid, const float pt, const float cent) + { + if (!qaSettingsGroup.cfgDoKappaAnalysis.value) { + return; + } + + const float kappa = getV0Kappa(pos, ele); + + const int pdgPos = std::abs(posmc.pdgCode()); + const int pdgEle = std::abs(elemc.pdgCode()); + + const bool posIsElectron = (pdgPos == PDG_t::kElectron); + const bool eleIsElectron = (pdgEle == PDG_t::kElectron); + const bool posIsPion = (pdgPos == PDG_t::kPiPlus); + const bool eleIsPion = (pdgEle == PDG_t::kPiPlus); + + if (photonid > 0) { + fRegistry.fill(HIST("Kappa/MC/v0"), kappa, pt, cent); + } else if (posIsElectron && eleIsElectron) { + fRegistry.fill(HIST("Kappa/MC/ee"), kappa, pt, cent); + } else if ((posIsPion && eleIsElectron) || (posIsElectron && eleIsPion)) { + fRegistry.fill(HIST("Kappa/MC/pie"), kappa, pt, cent); + } else if (posIsPion && eleIsPion) { + fRegistry.fill(HIST("Kappa/MC/pipi"), kappa, pt, cent); + } else { + fRegistry.fill(HIST("Kappa/MC/other"), kappa, pt, cent); + } + } + + Preslice perCollisionV0 = aod::v0photonkf::pmeventId; Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); Filter collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; using FilteredMyCollisions = soa::Filtered; using FilteredMyCollisionsMC = soa::Filtered; // same filters, they act column-wise - template - void process(FilteredMyCollisions const& collisions, TV0Photon const& v0photons, aod::V0Legs const&, TPerCollision const& perCollision) + template + void processRec(FilteredMyCollisions const& collisions, TV0Photon const& v0photons, aod::V0Legs const&) { for (const auto& collision : collisions) { initCCDB(collision); @@ -834,7 +870,7 @@ struct PCMQC { fV0PhotonCut.SetCentrality(centralities[cfgCentEstimator]); int nv0 = 0; - auto v0photons_coll = v0photons.sliceBy(perCollision, collision.globalIndex()); + auto v0photons_coll = v0photons.sliceBy(perCollisionV0, collision.globalIndex()); for (const auto& v0 : v0photons_coll) { auto pos = v0.template posTrack_as(); auto ele = v0.template negTrack_as(); @@ -844,6 +880,7 @@ struct PCMQC { } fillLossQAInfo<1>(v0); fillV0Info(v0); + fillKappaRec(pos, ele, v0.pt(), centralities[cfgCentEstimator]); fillMaterialBudgetInfo(v0); for (const auto& leg : {pos, ele}) { fillV0LegInfo(leg); @@ -859,15 +896,15 @@ struct PCMQC { } void processQC(FilteredMyCollisions const& collisions, MyV0Photons const& v0photons, aod::V0Legs const& v0legs) { - process(collisions, v0photons, v0legs, perCollisionV0); + processRec(collisions, v0photons, v0legs); } // end of process void processQCML(FilteredMyCollisions const& collisions, MyV0PhotonsML const& v0photonsML, aod::V0Legs const& v0legs) { - process(collisions, v0photonsML, v0legs, perCollisionV0ML); + processRec(collisions, v0photonsML, v0legs); } // end of ML process - template + template void fillV0InfoMC(TV0 const& v0, TMCV0 const& mcphoton, TMCLeg const& mcleg) { fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hPt"), v0.pt()); @@ -917,7 +954,7 @@ struct PCMQC { } float phi_cp = std::atan2(v0.vy(), v0.vx()); - o2::math_utils::bringTo02Pi(phi_cp); + RecoDecay::constrainAngle(phi_cp); float eta_cp = std::atanh(v0.vz() / std::sqrt(std::pow(v0.vx(), 2) + std::pow(v0.vy(), 2) + std::pow(v0.vz(), 2))); fRegistry.fill(HIST("V0/") + HIST(mcphoton_types[mctype]) + HIST("hsConvPoint"), v0.v0radius(), phi_cp, eta_cp); @@ -943,7 +980,7 @@ struct PCMQC { } } - template + template void fillV0LegInfoMC(TLeg const& leg) { fRegistry.fill(HIST("V0Leg/") + HIST(mcphoton_types[mctype]) + HIST("hPt"), leg.pt()); @@ -975,8 +1012,8 @@ struct PCMQC { fRegistry.fill(HIST("V0Leg/") + HIST(mcphoton_types[mctype]) + HIST("hRxyGen_DeltaPhi"), std::sqrt(std::pow(mcleg.vx(), 2) + std::pow(mcleg.vy(), 2)), leg.phi() - mcleg.phi()); } - template - void processMC(FilteredMyCollisionsMC const& collisions, TV0Photons const& v0photons, aod::EMMCParticles const& mcparticles, MyMCV0Legs const&, aod::EMMCEvents const&, TPerCollision const& percollision) + template + void processMC(FilteredMyCollisionsMC const& collisions, TV0Photons const& v0photons, aod::EMMCParticles const& mcparticles, MyMCV0Legs const&, aod::EMMCEvents const&) { for (const auto& collision : collisions) { initCCDB(collision); @@ -994,7 +1031,7 @@ struct PCMQC { fRegistry.fill(HIST("Event/after/hCollisionCounter"), 10.0); // accepted fV0PhotonCut.SetCentrality(centralities[cfgCentEstimator]); // set centrality for BDT response - auto v0photons_coll = v0photons.sliceBy(percollision, collision.globalIndex()); + auto v0photons_coll = v0photons.sliceBy(perCollisionV0, collision.globalIndex()); int ng_primary = 0, ng_wd = 0, ng_hs = 0, nee_pi0 = 0, nee_eta = 0; for (const auto& v0 : v0photons_coll) { auto pos = v0.template posTrack_as(); @@ -1020,6 +1057,9 @@ struct PCMQC { int photonid = FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 22, mcparticles); int pi0id = FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 111, mcparticles); // pi0 dalitz decay int etaid = FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 221, mcparticles); // eta dalitz decay + + fillKappaMC(pos, ele, posmc, elemc, photonid, v0.pt(), centralities[cfgCentEstimator]); + if (photonid < 0 && pi0id < 0 && etaid < 0) { continue; } @@ -1100,15 +1140,15 @@ struct PCMQC { void processPCMQCMC(FilteredMyCollisionsMC const& collisions, MyV0Photons const& v0photons, aod::EMMCParticles const& mcparticles, MyMCV0Legs const& mcv0legs, aod::EMMCEvents const& mcevents) { - processMC(collisions, v0photons, mcparticles, mcv0legs, mcevents, perCollisionV0); + processMC(collisions, v0photons, mcparticles, mcv0legs, mcevents); } // end of MC QC process void processPCMQCMCML(FilteredMyCollisionsMC const& collisions, MyV0PhotonsML const& v0photonsML, aod::EMMCParticles const& mcparticles, MyMCV0Legs const& mcv0legs, aod::EMMCEvents const& mcevents) { - processMC(collisions, v0photonsML, mcparticles, mcv0legs, mcevents, perCollisionV0ML); + processMC(collisions, v0photonsML, mcparticles, mcv0legs, mcevents); } // end of MC QC process with ML cuts - template + template void fillBinnedData(TBinnedData const& binned_data, const float weight = 1.f) { int xbin = 0, ybin = 0, zbin = 0; @@ -1166,7 +1206,7 @@ struct PCMQC { auto daughter = mcparticles.iteratorAt(mctrack.daughtersIds()[0]); // choose ele or pos. float rxy_gen_e = std::sqrt(std::pow(daughter.vx(), 2) + std::pow(daughter.vy(), 2)); float phi_cp = std::atan2(daughter.vy(), daughter.vx()); - o2::math_utils::bringTo02Pi(phi_cp); + RecoDecay::constrainAngle(phi_cp); float eta_cp = std::atanh(daughter.vz() / std::sqrt(std::pow(daughter.vx(), 2) + std::pow(daughter.vy(), 2) + std::pow(daughter.vz(), 2))); fRegistry.fill(HIST("Generated/hR_ConversionPhoton_wideR"), rxy_gen_e); @@ -1190,7 +1230,7 @@ struct PCMQC { } // end of collision loop } - template + template void fillLegQualityRecoQA(TTrack const& track) { if constexpr (survived) { diff --git a/PWGEM/PhotonMeson/Tasks/taskPhotonFlow.cxx b/PWGEM/PhotonMeson/Tasks/taskPhotonFlow.cxx new file mode 100644 index 00000000000..41f47171448 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/taskPhotonFlow.cxx @@ -0,0 +1,783 @@ +// 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 taskPhotonFlow.cxx +/// \brief Analysis task for photon flow with PCM or EMCal +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#include "PWGEM/PhotonMeson/Core/EMBitFlags.h" +#include "PWGEM/PhotonMeson/Core/EMCPhotonCut.h" +#include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" +#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" +#include "PWGEM/PhotonMeson/DataModel/EventTables.h" +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/EventHistograms.h" + +#include "Common/Core/EventPlaneHelper.h" +#include "Common/Core/RecoDecay.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::photon; + +enum QvecEstimator { + FT0M = 0, + FT0A = 1, + FT0C = 2, + TPCPos = 3, + TPCNeg = 4, + TPCTot = 5, + FV0A = 6 +}; + +enum CentralityEstimator { + None = 0, + CFT0A = 1, + CFT0C = 2, + CFT0M = 3, + NCentralityEstimators = 4 +}; + +enum Harmonics { + kNone = 0, + kDirect = 1, + kElliptic = 2, + kTriangluar = 3, + kQuadrangular = 4, + kPentagonal = 5, + kHexagonal = 6, + kHeptagonal = 7, + kOctagonal = 8 +}; + +enum class MapLevel { + kGood = 1, + kNoBad = 2, + kInEMC = 3, + kAll = 4 +}; + +struct TaskPhotonFlow { + // configurable for flow + Configurable harmonic{"harmonic", 2, "harmonic number"}; + Configurable qvecDetector{"qvecDetector", 0, "Detector for Q vector estimation (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; + Configurable qvecSubADetector{"qvecSubADetector", 3, "Sub A Detector for Q vector estimation for resolution (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; + Configurable qvecSubBDetector{"qvecSubBDetector", 4, "Sub B Detector for Q vector estimation for resolution (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; + Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3)"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable cfgEMCalMapLevelSameEvent{"cfgEMCalMapLevelSameEvent", 4, "Different levels of correction for the same event, the smaller number includes the level of the higher number (4: none, 3: only inside EMCal, 2: remove edges, 1: exclude bad channels)"}; + Configurable cfgDistanceToEdge{"cfgDistanceToEdge", 1, "Distance to edge in cells required for rotated cluster to be accepted"}; + Configurable cfgMaxQVector{"cfgMaxQVector", 20.f, "Maximum allowed absolute QVector value."}; + + // configurable axis + ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {100, 0., 20.}, "pT axis for the neutral meson"}; + ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {20, 0., 100.}, "centrality axis for the current event"}; + ConfigurableAxis thnConfigAxisCosDeltaPhi{"thnConfigAxisCosDeltaPhi", {8, -1., 1.}, "cos(delta phi) axis for the current event"}; + ConfigurableAxis thnConfigAxisM02{"thnConfigAxisM02", {200, 0., 5.}, "M02 axis for the EMCal cluster"}; + ConfigurableAxis thnConfigAxisKappa{"thnConfigAxisKappa", {200, -10., 10.}, "Kappa axis for the EMCal cluster"}; + ConfigurableAxis thnConfigAxisEnergyCalib{"thnConfigAxisEnergyCalib", {200, 0., 20.}, "energy axis for the emcal clusters for the calibration process"}; + + EMPhotonEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcuts"; + Configurable cfgZvtxMax{"cfgZvtxMax", 10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgRequireEMCReadoutInMB{"cfgRequireEMCReadoutInMB", true, "require the EMC to be read out in an MB collision (kTVXinEMC)"}; + Configurable cfgRequireEMCHardwareTriggered{"cfgRequireEMCHardwareTriggered", false, "require the EMC to be hardware triggered (kEMC7 or kDMC7)"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -1, "min. FT0C occupancy"}; + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + Configurable cfgMinCent{"cfgMinCent", 0, "min. centrality (%)"}; + Configurable cfgMaxCent{"cfgMaxCent", 90, "max. centrality (%)"}; + Configurable onlyKeepWeightedEvents{"onlyKeepWeightedEvents", false, "flag to keep only weighted events (for JJ MCs) and remove all MB events (with weight = 1)"}; + } eventcuts; + + EMCPhotonCut fEMCCut; + struct : ConfigurableGroup { + std::string prefix = "emccuts"; + Configurable clusterDefinition{"clusterDefinition", "kV3MostSplitSmallestTimeDiff", "Clusterizer to be selected, e.g. V3Default"}; + Configurable cfgEMCminTime{"cfgEMCminTime", -25., "Minimum cluster time for EMCal time cut"}; + Configurable cfgEMCmaxTime{"cfgEMCmaxTime", +30., "Maximum cluster time for EMCal time cut"}; + Configurable cfgEMCminM02{"cfgEMCminM02", 0.1, "Minimum M02 for EMCal M02 cut"}; + Configurable cfgEMCmaxM02{"cfgEMCmaxM02", 0.7, "Maximum M02 for EMCal M02 cut"}; + Configurable cfgEMCminE{"cfgEMCminE", 0.7, "Minimum cluster energy for EMCal energy cut"}; + Configurable cfgEMCminNCell{"cfgEMCminNCell", 1, "Minimum number of cells per cluster for EMCal NCell cut"}; + Configurable> cfgEMCTMEta{"cfgEMCTMEta", {0.01f, 4.07f, -2.5f}, "|eta| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable> cfgEMCTMPhi{"cfgEMCTMPhi", {0.015f, 3.65f, -2.f}, "|phi| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable> emcSecTMEta{"emcSecTMEta", {0.01f, 4.07f, -2.5f}, "|eta| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable> emcSecTMPhi{"emcSecTMPhi", {0.015f, 3.65f, -2.f}, "|phi| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable cfgEMCEoverp{"cfgEMCEoverp", 1.75, "Minimum cluster energy over track momentum for EMCal track matching"}; + Configurable cfgEMCUseExoticCut{"cfgEMCUseExoticCut", true, "FLag to use the EMCal exotic cluster cut"}; + Configurable cfgEMCUseTM{"cfgEMCUseTM", false, "flag to use EMCal track matching cut or not"}; + Configurable emcUseSecondaryTM{"emcUseSecondaryTM", false, "flag to use EMCal secondary track matching cut or not"}; + Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; + Configurable useEMCal{"useEMCal", false, "flag to use EMCal clusters"}; + Configurable useDCal{"useDCal", false, "flag to use DCal clusters"}; + } emccuts; + + V0PhotonCut fV0PhotonCut; + struct : o2::framework::ConfigurableGroup { + std::string prefix = "PCMcuts"; + o2::framework::Configurable requireV0WithITSTPC{"requireV0WithITSTPC", false, "flag to enforce V0s have ITS and TPC"}; + o2::framework::Configurable requireV0WithITSOnly{"requireV0WithITSOnly", false, "flag to select V0s with ITSonly tracks"}; + o2::framework::Configurable requireV0WithTPCOnly{"requireV0WithTPCOnly", false, "flag to select V0s with TPConly tracks"}; + o2::framework::Configurable minPtV0{"minPtV0", 0.1, "min pT for v0 photons at PV"}; + o2::framework::Configurable maxPtV0{"maxPtV0", 1e+10, "max pT for v0 photons at PV"}; + o2::framework::Configurable minEtaV0{"minEtaV0", -0.8, "min eta for v0 photons at PV"}; + o2::framework::Configurable maxEtaV0{"maxEtaV0", 0.8, "max eta for v0 photons at PV"}; + o2::framework::Configurable minRV0{"minRV0", 4.0, "min v0 radius"}; + o2::framework::Configurable maxRV0{"maxRV0", 90.0, "max v0 radius"}; + o2::framework::Configurable maxAlphaAP{"maxAlphaAP", 0.95, "max alpha for AP cut"}; + o2::framework::Configurable maxQtAP{"maxQtAP", 0.01, "max qT for AP cut"}; + o2::framework::Configurable minCosPA{"minCosPA", 0.999, "min V0 CosPA"}; + o2::framework::Configurable maxPCA{"maxPCA", 1.5, "max distance btween 2 legs"}; + o2::framework::Configurable maxChi2KF{"maxChi2KF", 1.e+10f, "max chi2/ndf with KF"}; + o2::framework::Configurable rejectV0onITSib{"rejectV0onITSib", true, "flag to reject V0s on ITSib"}; + o2::framework::Configurable applyPrefilter{"applyPrefilter", false, "flag to apply prefilter to V0"}; + + o2::framework::Configurable minNClusterTPC{"minNClusterTPC", 0, "min NCluster TPC"}; + o2::framework::Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 40, "min ncrossed rows in TPC"}; + o2::framework::Configurable minNCrossedRowsOverFindableClustersTPC{"minNCrossedRowsOverFindableClustersTPC", 0.8f, "min fraction of crossed rows over findable clusters in TPC"}; + o2::framework::Configurable maxFracSharedClustersTPC{"maxFracSharedClustersTPC", 999.f, "max fraction of shared clusters in TPC"}; + o2::framework::Configurable minNClusterITS{"minNClusterITS", 0, "min NCluster ITS"}; + o2::framework::Configurable minMeanClusterSizeITSob{"minMeanClusterSizeITSob", 0.f, "min ITSob"}; + o2::framework::Configurable maxMeanClusterSizeITSob{"maxMeanClusterSizeITSob", 16.f, "max ITSob"}; + o2::framework::Configurable macChi2TPC{"macChi2TPC", 4.f, "max chi2/NclsTPC"}; + o2::framework::Configurable macChi2ITS{"macChi2ITS", 36.f, "max chi2/NclsITS"}; + o2::framework::Configurable minTPCNSigmaEl{"minTPCNSigmaEl", -3.0f, "min. TPC n sigma for electron"}; + o2::framework::Configurable maxTPCNSigmaEl{"maxTPCNSigmaEl", +3.0f, "max. TPC n sigma for electron"}; + o2::framework::Configurable disableITSOnly{"disableITSOnly", false, "flag to disable ITSonly tracks"}; + o2::framework::Configurable disableTPCOnly{"disableTPCOnly", false, "flag to disable TPConly tracks"}; + o2::framework::Configurable doQA{"doQA", false, "flag to set QA flag."}; + } pcmcuts; + + struct : ConfigurableGroup { + std::string prefix = "correctionConfig"; + Configurable cfgSpresoPath{"cfgSpresoPath", "Users/m/mhemmer/EM/Flow/Resolution", "Path to SP resolution file"}; + Configurable cfgApplySPresolution{"cfgApplySPresolution", false, "Apply resolution correction"}; + } correctionConfig; + + SliceCache cache; + EventPlaneHelper epHelper; + o2::framework::Service ccdb{}; + int runNow = 0; + int runBefore = -1; + + static constexpr float MaxPhiEMCal = 3.9f; // exatly the middle between EMCal and DCal + static constexpr uint16_t MaxPhiEMCalUint = static_cast(39000u); // exatly the middle between EMCal and DCal but as uint16_t that is used for storing phi values in derived data + + using EMCalPhotons = soa::Join; + using PCMPhotons = soa::Join; + using Colls = soa::Join; + + Partition emcalPhotons = aod::mincluster::storedPhi < MaxPhiEMCalUint; + Partition dcalPhotons = aod::mincluster::storedPhi >= MaxPhiEMCalUint; + + static constexpr std::size_t NQVecEntries = 6; + + PresliceOptional perCollisionEMC = o2::aod::emccluster::pmeventId; + PresliceOptional perCollisionPCM = aod::v0photonkf::pmeventId; + PresliceOptional perEMCClusterMT = o2::aod::mintm::minClusterId; + PresliceOptional perEMCClusterMS = o2::aod::mintm::minClusterId; + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + o2::emcal::Geometry* emcalGeom = nullptr; + TH1D* h1SPResolution = nullptr; + // Constants for eta and phi ranges for the look up table + static constexpr double EtaMin = -0.75, etaMax = 0.75; + static constexpr int NBinsEta = 150; // 150 bins for eta + + static constexpr double PhiMin = 1.35, phiMax = 5.75; + static constexpr int NBinsPhi = 440; // (440 bins = 0.01 step size covering most regions) + + std::array lookupTable1D{}; + float epsilon = 1.e-8; + + // To access the 1D array + static 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"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(-eventcuts.cfgZvtxMax, +eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetRequireEMCReadoutInMB(eventcuts.cfgRequireEMCReadoutInMB); + fEMEventCut.SetRequireEMCHardwareTriggered(eventcuts.cfgRequireEMCHardwareTriggered); + } + + void defineEMCCut() + { + fEMCCut = EMCPhotonCut("fEMCCut", "fEMCCut"); + + fEMCCut.SetTrackMatchingEtaParams(emccuts.cfgEMCTMEta->at(0), emccuts.cfgEMCTMEta->at(1), emccuts.cfgEMCTMEta->at(2)); + fEMCCut.SetTrackMatchingPhiParams(emccuts.cfgEMCTMPhi->at(0), emccuts.cfgEMCTMPhi->at(1), emccuts.cfgEMCTMPhi->at(2)); + + fEMCCut.SetSecTrackMatchingEtaParams(emccuts.emcSecTMEta->at(0), emccuts.emcSecTMEta->at(1), emccuts.emcSecTMEta->at(2)); + fEMCCut.SetSecTrackMatchingPhiParams(emccuts.emcSecTMPhi->at(0), emccuts.emcSecTMPhi->at(1), emccuts.emcSecTMPhi->at(2)); + fEMCCut.SetMinEoverP(emccuts.cfgEMCEoverp); + + fEMCCut.SetMinE(emccuts.cfgEMCminE); + fEMCCut.SetMinNCell(emccuts.cfgEMCminNCell); + fEMCCut.SetM02Range(emccuts.cfgEMCminM02, emccuts.cfgEMCmaxM02); + fEMCCut.SetTimeRange(emccuts.cfgEMCminTime, emccuts.cfgEMCmaxTime); + fEMCCut.SetUseExoticCut(emccuts.cfgEMCUseExoticCut); + fEMCCut.SetClusterizer(emccuts.clusterDefinition); + fEMCCut.SetUseTM(emccuts.cfgEMCUseTM.value); // disables or enables TM + fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM + fEMCCut.SetDoQA(emccuts.cfgEnableQA.value); + } + + void definePCMCut() + { + fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); + + // for v0 + fV0PhotonCut.SetV0PtRange(pcmcuts.minPtV0, pcmcuts.maxPtV0); + fV0PhotonCut.SetV0EtaRange(pcmcuts.minEtaV0, pcmcuts.maxEtaV0); + fV0PhotonCut.SetMinCosPA(pcmcuts.minCosPA); + fV0PhotonCut.SetMaxPCA(pcmcuts.maxPCA); + fV0PhotonCut.SetMaxChi2KF(pcmcuts.maxChi2KF); + fV0PhotonCut.SetRxyRange(pcmcuts.minRV0, pcmcuts.maxRV0); + fV0PhotonCut.SetAPRange(pcmcuts.maxAlphaAP, pcmcuts.maxQtAP); + fV0PhotonCut.RejectITSib(pcmcuts.rejectV0onITSib); + + // for track + fV0PhotonCut.SetMinNClustersTPC(pcmcuts.minNClusterTPC); + fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.minNCrossedRowsTPC); + fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(pcmcuts.minNCrossedRowsOverFindableClustersTPC); + fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.maxFracSharedClustersTPC); + fV0PhotonCut.SetChi2PerClusterTPC(0.f, pcmcuts.macChi2TPC); + fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.minTPCNSigmaEl, pcmcuts.maxTPCNSigmaEl); + fV0PhotonCut.SetChi2PerClusterITS(0.f, pcmcuts.macChi2ITS); + fV0PhotonCut.SetNClustersITS(pcmcuts.minNClusterITS, 7); + fV0PhotonCut.SetMeanClusterSizeITSob(pcmcuts.minMeanClusterSizeITSob, pcmcuts.maxMeanClusterSizeITSob); + fV0PhotonCut.SetDisableITSonly(pcmcuts.disableITSOnly); + fV0PhotonCut.SetDisableTPConly(pcmcuts.disableTPCOnly); + fV0PhotonCut.SetRequireITSTPC(pcmcuts.requireV0WithITSTPC); + fV0PhotonCut.SetRequireITSonly(pcmcuts.requireV0WithITSOnly); + fV0PhotonCut.SetRequireTPConly(pcmcuts.requireV0WithTPCOnly); + + fV0PhotonCut.setDoQA(pcmcuts.doQA.value); + } + + void init(InitContext&) + { + if (harmonic != kElliptic && harmonic != kTriangluar) { + LOG(info) << "Harmonic was set to " << harmonic << " but can only be 2 or 3!"; + } + + defineEMEventCut(); + defineEMCCut(); + fEMCCut.addQAHistograms(®istry); + definePCMCut(); + fV0PhotonCut.addQAHistograms(®istry); + o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(®istry); + + const AxisSpec thnAxisPt{thnConfigAxisPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec thnAxisCent{thnConfigAxisCent, "Centrality (%)"}; + const AxisSpec thnAxisCosDeltaPhi{thnConfigAxisCosDeltaPhi, std::format("cos({}(#varphi - #Psi_{{sub}}))", harmonic.value)}; + const AxisSpec thnAxisM02{thnConfigAxisM02, "M_{02}"}; + const AxisSpec thnAxisKappa{thnConfigAxisKappa, "#Kappa"}; + const AxisSpec thAxisClusterEnergy{thnConfigAxisPt, "#it{E} (GeV)"}; + const AxisSpec thAxisEnergyCalib{thnConfigAxisEnergyCalib, "#it{E}_{clus} (GeV)"}; + const AxisSpec thAxisEnergy{1000, 0., 100., "#it{E}_{clus} (GeV)"}; + const AxisSpec thAxisEta{320, -0.8, 0.8, "#eta"}; + const AxisSpec thAxisPhi{500, 0, 2 * 3.14159, "phi"}; + + if (emccuts.cfgEnableQA.value) { + registry.add("clusterQA/hEClusterBefore", "Histo for cluster energy before cuts", HistType::kTH1D, {thAxisClusterEnergy}); + registry.add("clusterQA/hEClusterAfter", "Histo for cluster energy after cuts", HistType::kTH1D, {thAxisClusterEnergy}); + registry.add("clusterQA/hClusterEtaPhiBefore", "hClusterEtaPhiBefore", HistType::kTH2D, {thAxisPhi, thAxisEta}); + registry.add("clusterQA/hClusterEtaPhiAfter", "hClusterEtaPhiAfter", HistType::kTH2D, {thAxisPhi, thAxisEta}); + } + + if (doprocessEMC) { + registry.add("EMCal/hSparsePhotonFlow", " vs M_{02} vs p_T vs cent", HistType::kTProfile3D, {thnAxisM02, thnAxisPt, thnAxisCent}); + registry.add("EMCal/hSparsePhoton", "M_{02} vs p_T vs cent", HistType::kTH3D, {thnAxisM02, thnAxisPt, thnAxisCent}); + } + if (doprocessPCM) { + registry.add("PCM/hSparsePhotonFlow", " vs M_{02} vs p_T vs cent", HistType::kTProfile3D, {thnAxisM02, thnAxisPt, thnAxisCent}); + registry.add("PCM/hSparsePhoton", "M_{02} vs p_T vs cent", HistType::kTH3D, {thnAxisM02, thnAxisPt, thnAxisCent}); + } + + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + LOG(info) << "thnConfigAxisPt.value[1] = " << thnConfigAxisPt.value[1] << " thnConfigAxisPt.value.back() = " << thnConfigAxisPt.value.back(); + + }; // end init + + /// \brief Check whether a photon (by its phi) falls in the EMCal or DCal acceptance + /// \param phi azimuthal angle of the photon + static bool isEMCalRegion(float phi) + { + return phi < MaxPhiEMCal; + } + + /// Change radians to degree + /// \param angle in radians + /// \return angle in degree + static constexpr float getAngleDegree(float angle) + { + return angle * o2::constants::math::Rad2Deg; + } + + /// Compute the delta psi in the range [0, pi/harmonic] + /// \param psi1 is the first angle + /// \param psi2 is the second angle + float getDeltaPsiInRange(float psi1, float psi2) + { + float deltaPsi = psi1 - psi2; + return RecoDecay::constrainAngle(deltaPsi, 0.f, harmonic); + } + + /// Get the centrality + /// \param collision is the collision with the centrality information + template + float getCentrality(TCollision const& collision) + { + float cent = -999.; + switch (centEstimator) { + case CentralityEstimator::CFT0M: + cent = collision.centFT0M(); + break; + case CentralityEstimator::CFT0A: + cent = collision.centFT0A(); + break; + case CentralityEstimator::CFT0C: + cent = collision.centFT0C(); + break; + default: + LOG(warning) << "Centrality estimator not valid. Possible values are T0M, T0A, T0C. Fallback to T0C"; + cent = collision.centFT0C(); + break; + } + return cent; + } + + /// Get all used Q vector + /// \param collision is the collision with the Q vector information + template + std::array getAllQvec(TCollision const& collision) + { + // Retrieve the Q vectors using the helper function for each detector + auto [xQVecMain, yQVecMain] = getQvec(collision, qvecDetector); + auto [xQVecSubA, yQVecSubA] = getQvec(collision, qvecSubADetector); + auto [xQVecSubB, yQVecSubB] = getQvec(collision, qvecSubBDetector); + + return {xQVecMain, yQVecMain, xQVecSubA, yQVecSubA, xQVecSubB, yQVecSubB}; + } + + /// Get the Q vector + /// \param collision is the collision with the Q vector information + template + std::pair getQvec(TCollision const& collision, int detector) + { + float xQVec = -999.f; + float yQVec = -999.f; + + switch (detector) { + case QvecEstimator::FT0M: + if (harmonic == kElliptic) { + xQVec = collision.q2xft0m(); + yQVec = collision.q2yft0m(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0m(); + yQVec = collision.q3yft0m(); + } + break; + case QvecEstimator::FT0A: + if (harmonic == kElliptic) { + xQVec = collision.q2xft0a(); + yQVec = collision.q2yft0a(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0a(); + yQVec = collision.q3yft0a(); + } + break; + case QvecEstimator::FT0C: + if (harmonic == kElliptic) { + xQVec = collision.q2xft0c(); + yQVec = collision.q2yft0c(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0c(); + yQVec = collision.q3yft0c(); + } + break; + case QvecEstimator::TPCPos: + if (harmonic == kElliptic) { + xQVec = collision.q2xbpos(); + yQVec = collision.q2ybpos(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xbpos(); + yQVec = collision.q3ybpos(); + } + break; + case QvecEstimator::TPCNeg: + if (harmonic == kElliptic) { + xQVec = collision.q2xbneg(); + yQVec = collision.q2ybneg(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xbneg(); + yQVec = collision.q3ybneg(); + } + break; + case QvecEstimator::TPCTot: + if (harmonic == kElliptic) { + xQVec = collision.q2xbtot(); + yQVec = collision.q2ybtot(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xbtot(); + yQVec = collision.q3ybtot(); + } + break; + case QvecEstimator::FV0A: + if (harmonic == kElliptic) { + xQVec = collision.q2xfv0a(); + yQVec = collision.q2yfv0a(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xfv0a(); + yQVec = collision.q3yfv0a(); + } + break; + default: + LOG(warning) << "Q vector estimator not valid. Falling back to FT0M"; + if (harmonic == kElliptic) { + xQVec = collision.q2xft0m(); + yQVec = collision.q2yft0m(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0m(); + yQVec = collision.q3yft0m(); + } + break; + } + return {xQVec, yQVec}; + } + + /// Check if the QVector values are within reasonable range + /// \param collision is the collision with the Q vector information + bool isQvecGood(std::array const& QVecs) + { + bool isgood = true; + for (const auto& QVec : QVecs) { + if (std::fabs(QVec) > cfgMaxQVector) { + isgood = false; + break; + } + } + return isgood; + } + + bool isTooCloseToEdge(const int cellID, const int DistanceToBorder = 1) + { + if (DistanceToBorder <= 0) { + return false; + } + if (cellID < 0) { + return true; + } + + int iBadCell = -1; + + // check distance to border in case the cell is okay + auto [iSupMod, iMod, iPhi, iEta] = emcalGeom->GetCellIndex(cellID); + auto [irow, icol] = emcalGeom->GetCellPhiEtaIndexInSModule(iSupMod, iMod, iPhi, iEta); + + // Check rows/phi + int iRowLast = 24; + if (emcalGeom->GetSMType(iSupMod) == o2::emcal::EMCALSMType::EMCAL_HALF) { + iRowLast /= 2; // 2/3 sm case + } else if ((emcalGeom->GetSMType(iSupMod) == o2::emcal::EMCALSMType::EMCAL_THIRD) || + (emcalGeom->GetSMType(iSupMod) == o2::emcal::EMCALSMType::DCAL_EXT)) { + iRowLast /= 3; // 1/3 sm case + } + + if (irow < DistanceToBorder || (iRowLast - irow) <= DistanceToBorder) { + iBadCell = 1; + } + + if (iBadCell > 0) { + return true; + } + return false; + } + + template + bool isCellMasked(int cellID, TCollision const& collision) + { + bool masked = false; + auto maskStatus = collision.badChannelMap().getChannelStatus(cellID); + masked = (maskStatus != o2::emcal::BadChannelMap::MaskType_t::GOOD_CELL); + return masked; + } + + template + float getV0Kappa(TLeg const& pos, TLeg const& ele) + { + float kappa = (std::fabs(pos.tpcNSigmaEl()) + std::fabs(ele.tpcNSigmaEl())) / 2.f + (2.f * pos.tpcNSigmaEl() + ele.tpcNSigmaEl()); + return kappa; + } + + template + void initCCDB(TCollision const& collision) + { + // Load EMCal geometry + emcalGeom = o2::emcal::Geometry::GetInstanceFromRunNumber(collision.runNumber()); + lookupTable1D.fill(-1); + double binWidthEta = (etaMax - EtaMin) / NBinsEta; + double binWidthPhi = (phiMax - PhiMin) / NBinsPhi; + + if (cfgEMCalMapLevelSameEvent >= static_cast(MapLevel::kAll)) { + // in this case we do not want to check the clusters, so just say thery are all good. + lookupTable1D.fill(0); // good + } else { + 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, collision)) { + 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 check if standard event cuts + FT0 occupancy + centrality + QVec good is + /// \param collision collision that will be checked + /// \return true if collision survives all checks, otherwise false + template + bool isFullEventSelected(TCollision const& collision, bool fillHisto = false) + { + if (fillHisto) { + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); + } + if (!(fEMEventCut.IsSelected(collision))) { + // general event selection + return false; + } + if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { + // occupancy selection + return false; + } + float cent = getCentrality(collision); + if (cent < eventcuts.cfgMinCent || cent > eventcuts.cfgMaxCent) { + // event selection + return false; + } + if (!isQvecGood(getAllQvec(collision))) { + // selection based on QVector + return false; + } + if (fillHisto) { + 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 + } + return true; + } + + void processEMC(Colls const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) + { + if (clusters.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + + for (const auto& collision : collisions) { + float cent = getCentrality(collision); + if (!isFullEventSelected(collision, true)) { + continue; + } + runNow = collision.runNumber(); + if (runNow != runBefore) { + initCCDB(collision); + runBefore = runNow; + } + auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); + + for (const auto& photon : photonsPerCollision) { + if (emccuts.cfgEnableQA.value) { + registry.fill(HIST("clusterQA/hEClusterBefore"), photon.corrE()); // before cuts + registry.fill(HIST("clusterQA/hClusterEtaPhiBefore"), photon.phi(), photon.eta()); // before cuts + } + if (!(emcFlags.test(photon.globalIndex()))) { + continue; + } + if ((emccuts.useEMCal.value && !emccuts.useDCal.value && (!isEMCalRegion(photon.phi()))) || (!emccuts.useEMCal.value && emccuts.useDCal.value && (isEMCalRegion(photon.phi())))) { + continue; + } + if (cfgDistanceToEdge.value > 0 && (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelSameEvent.value)) { + continue; + } + if (emccuts.cfgEnableQA.value) { + registry.fill(HIST("clusterQA/hEClusterAfter"), photon.corrE()); // accepted after cuts + registry.fill(HIST("clusterQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // after cuts + } + + auto [xQVec, yQVec] = getQvec(collision, qvecDetector); + + float phiCand = photon.phi(); + + float cosNPhi = std::cos(harmonic * phiCand); + float sinNPhi = std::sin(harmonic * phiCand); + float scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; + + if (correctionConfig.cfgApplySPresolution.value) { + scalprodCand = scalprodCand / h1SPResolution->GetBinContent(h1SPResolution->FindBin(cent + epsilon)); + } + registry.fill(HIST("EMCal/hSparsePhotonFlow"), photon.m02(), photon.corrPt(), cent, scalprodCand); + registry.fill(HIST("EMCal/hSparsePhoton"), photon.m02(), photon.corrPt(), cent); + } // end of loop over single cluster + } // end of loop over collisions + } // processEMC + PROCESS_SWITCH(TaskPhotonFlow, processEMC, "Process single EMCal clusters as function of M02", false); + + void processPCM(Colls const& collisions, PCMPhotons const& photons, aod::V0Legs const& legs) + { + if (photons.size() <= 0 || legs.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + auto ele = legs.begin(); + auto pos = legs.begin(); + EMBitFlags v0flags(photons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, photons, ®istry); + for (const auto& collision : collisions) { + + if (!isFullEventSelected(collision, true)) { + continue; + } + + float cent = getCentrality(collision); + + runNow = collision.runNumber(); + if (runNow != runBefore) { + initCCDB(collision); + runBefore = runNow; + } + + auto photonsPerCollision = photons.sliceBy(perCollisionPCM, collision.globalIndex()); + for (const auto& photon : photonsPerCollision) { + if (!(v0flags.test(photon.globalIndex()))) { + continue; + } + ele.setCursor(photon.negTrackId()); + pos.setCursor(photon.posTrackId()); + + auto [xQVec, yQVec] = getQvec(collision, qvecDetector); + + float phiCand = photon.phi(); + + float cosNPhi = std::cos(harmonic * phiCand); + float sinNPhi = std::sin(harmonic * phiCand); + float scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; + + if (correctionConfig.cfgApplySPresolution.value) { + scalprodCand = scalprodCand / h1SPResolution->GetBinContent(h1SPResolution->FindBin(cent + epsilon)); + } + + float kappa = getV0Kappa(pos, ele); + registry.fill(HIST("PCM/hSparsePhotonFlow"), kappa, photon.corrPt(), cent, scalprodCand); + registry.fill(HIST("PCM/hSparsePhoton"), kappa, photon.corrPt(), cent); + } + } // end of loop over collisions + } + PROCESS_SWITCH(TaskPhotonFlow, processPCM, "Process PCM Pi0 candidates", true); + +}; // End struct TaskPhotonFlow + +WorkflowSpec defineDataProcessing(ConfigContext const& context) +{ + return WorkflowSpec{adaptAnalysisTask(context)}; +} diff --git a/PWGEM/PhotonMeson/Utils/PCMUtilities.h b/PWGEM/PhotonMeson/Utils/PCMUtilities.h index 8c88d20ccad..e69703a3ec5 100644 --- a/PWGEM/PhotonMeson/Utils/PCMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PCMUtilities.h @@ -43,7 +43,7 @@ inline bool checkAP(const float alpha, const float qt, const float alpha_max = 0 return (ellipse < 1.0); } //_______________________________________________________________________ -inline float v0_alpha(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) +inline float v0Alpha(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { float momTot = RecoDecay::p(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); float lQlNeg = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}) / momTot; @@ -51,7 +51,7 @@ inline float v0_alpha(float pxpos, float pypos, float pzpos, float pxneg, float return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); // longitudinal momentum asymmetry of v0 } //_______________________________________________________________________ -inline float v0_qt(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) +inline float v0Qt(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { float momTot = RecoDecay::p2(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); float dp = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}); @@ -262,4 +262,11 @@ inline V0TruthClass classifyV0Truth(TTrack const& pos, TTrack const& ele, TMCPar return kV0OtherFake; } +template +float getV0Kappa(TLeg const& pos, TLeg const& ele) +{ + float kappa = (std::fabs(pos.tpcNSigmaEl()) + std::fabs(ele.tpcNSigmaEl())) / 2.f + (2.f * (pos.tpcNSigmaEl() + ele.tpcNSigmaEl())); + return kappa; +} + #endif // PWGEM_PHOTONMESON_UTILS_PCMUTILITIES_H_