diff --git a/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx b/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx index b5c1a5d8fde..fefac20d019 100644 --- a/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx +++ b/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx @@ -8,33 +8,61 @@ // 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. -// -// ======================== -// -// This code loops over v0 photons for studying material budget. -// Please write to: daiki.sekihata@cern.ch -#include -#include +/// \file MaterialBudget.cxx +/// \brief Task to analyse and calculate the material budget weights +/// \author S. Mrozinski, smrozins@cern.ch -#include "TString.h" -#include "Math/Vector4D.h" -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoAHelpers.h" -#include "ReconstructionDataFormats/Track.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/PIDResponse.h" -#include "Common/Core/RecoDecay.h" -#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" -#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" -#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" -#include "PWGEM/PhotonMeson/Core/PairCut.h" +#include "PWGEM/Dilepton/Utils/MCUtilities.h" #include "PWGEM/PhotonMeson/Core/CutsLibrary.h" +#include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/HistogramsLibrary.h" +#include "PWGEM/PhotonMeson/Core/PairCut.h" +#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/MCUtilities.h" +#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TPCVDriftManager.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/McCollisionExtra.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include using namespace o2; using namespace o2::aod; @@ -43,366 +71,1233 @@ using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::aod::pwgem::photonmeson::photonpair; using namespace o2::aod::pwgem::photon; +using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; +using namespace o2::aod::pwgem::dilepton::utils::mcutil; +using namespace o2::aod::pwgem::dilepton::utils; +using o2::constants::math::TwoPI; using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; +using MyCollisionsMC = soa::Join; +using MyCollisionMC = MyCollisionsMC::iterator; + +using MyMCCollisions = soa::Join; +using MyMCCollision = MyMCCollisions::iterator; + using MyV0Photons = soa::Join; using MyV0Photon = MyV0Photons::iterator; +using MyPrimaryElectrons = soa::Filtered>; +using MyPrimaryElectron = MyPrimaryElectrons::iterator; + +using MyMCV0Legs = soa::Join; +using MyMCV0Leg = MyMCV0Legs::iterator; + +using MyTracks = soa::Join; struct MaterialBudget { Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; Configurable cfgCentMin{"cfgCentMin", 0, "min. centrality"}; - Configurable cfgCentMax{"cfgCentMax", 999, "max. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + + struct SimplePhoton { + float pt; + float eta; + float phi; + float vz; + float r; + }; + + static constexpr int MaxMixEvents = 5; + std::map, std::deque>> mixingPools; - Configurable fConfigTagCuts{"cfgTagCuts", "qc", "Comma separated list of V0 photon cuts for tag"}; - Configurable fConfigProbeCuts{"cfgProbeCuts", "qc,wwire_ib", "Comma separated list of V0 photon cuts for probe"}; - Configurable fConfigPairCuts{"cfgPairCuts", "nocut", "Comma separated list of pair cuts"}; - Configurable fDoMixing{"DoMixing", false, "do event mixing"}; + std::vector centBinEdges = {0, 10, 30, 50, 90, 200}; + std::vector zvtxBinEdges = {-10, -5, 0, 5, 10}; + + std::tuple getPoolBin(float cent, float zvtx) + { + int centbin = std::lower_bound(centBinEdges.begin(), centBinEdges.end(), cent) - centBinEdges.begin() - 1; + centbin = std::max(0, std::min(centbin, static_cast(centBinEdges.size()) - 2)); + + int zbin = std::lower_bound(zvtxBinEdges.begin(), zvtxBinEdges.end(), zvtx) - zvtxBinEdges.begin() - 1; + zbin = std::max(0, std::min(zbin, static_cast(zvtxBinEdges.size()) - 2)); + return {centbin, zbin}; + } + + HistogramRegistry registry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + static constexpr std::string_view ItsClsNames[] = { + "ITSCls0", "ITSCls1", "ITSCls2", "ITSCls3", + "ITSCls4", "ITSCls5", "ITSCls6", "ITSCls7"}; + + static constexpr std::string_view EventTypes[2] = {"before/", "after/"}; + + Configurable cfgPlotBremsstrahlung{"cfgPlotBremsstrahlung", false, "produce plots to study Bremsstrahlung"}; + Configurable cfgPlotResolution{"cfgPlotResolution", false, "produce plots to study resolution"}; + Configurable cfgPlotPurity{"cfgPlotPurity", false, "produce plots to study purity"}; + Configurable cfgPlotMBDetailed{"cfgPlotMBDetailed", false, "produce plots to study material distribution distribution"}; + Configurable cfgPlotMBGeneral{"cfgPlotMBGeneral", true, "produce plots to study material distribution"}; + Configurable cfgPlotMBCollisions{"cfgPlotMBCollisions", false, "produce plots to study material distribution if collision association is wrong"}; + + ConfigurableAxis binsPt{"binsPt", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 2.0}, ""}; + const AxisSpec axisPt{binsPt, "#it{p}_{T} [GeV/c]"}; - Configurable fConfigEMEventCut{"cfgEMEventCut", "minbias", "em event cut"}; // only 1 event cut per wagon EMPhotonEventCut fEMEventCut; - static constexpr std::string_view event_types[2] = {"before", "after"}; + struct : ConfigurableGroup { + std::string prefix = "eventcut_group"; + 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", true, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", true, "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 cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. occupancy"}; + Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. occupancy"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2, "min. FT0C occupancy"}; + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; + Configurable cfgRequireNoCollInTimeRangeStrict{"cfgRequireNoCollInTimeRangeStrict", false, "require no collision in time range strict"}; + Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "require no collision in time range standard"}; + Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "require no collision in time range strict"}; + Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "require no HM collision in previous ITS ROF"}; + } eventcuts; + + DalitzEECut fDileptonCut; + struct : ConfigurableGroup { + std::string prefix = "dileptoncut_group"; + Configurable cfgMinMass{"cfgMinMass", 0.0, "min mass"}; + Configurable cfgMaxMass{"cfgMaxMass", 0.1, "max mass"}; + Configurable cfgApplyPhiv{"cfgApplyPhiv", true, "flag to apply phiv cut"}; + Configurable cfgRequireItsibAny{"cfgRequireItsibAny", false, "flag to require ITS ib any hits"}; + Configurable cfgRequireItsib1st{"cfgRequireItsib1st", true, "flag to require ITS ib 1st hit"}; + Configurable cfgPhivSlope{"cfgPhivSlope", 0.0185, "slope for m vs. phiv"}; + Configurable cfgPhivIntercept{"cfgPhivIntercept", -0.0280, "intercept for m vs. phiv"}; - OutputObj fOutputEvent{"Event"}; - OutputObj fOutputV0{"V0"}; - OutputObj fOutputPair{"Pair"}; // 2-photon pair - THashList* fMainList = new THashList(); + Configurable cfgMinPtTrack{"cfgMinPtTrack", 0.1, "min pT for single track"}; + Configurable cfgMaxEtaTrack{"cfgMaxEtaTrack", 0.8, "max eta for single track"}; + Configurable cfgMinnclusterTPC{"cfgMinnclusterTPC", 0, "min ncluster tpc"}; + Configurable cfgMinnclusterITS{"cfgMinnclusterITS", 5, "min ncluster its"}; + Configurable cfgMinncrossedrows{"cfgMinncrossedrows", 70, "min ncrossed rows"}; + Configurable cfgMaxchi2TPC{"cfgMaxchi2TPC", 4.0, "max chi2/NclsTPC"}; + Configurable cfgMaxchi2ITS{"cfgMaxchi2ITS", 36.0, "max chi2/NclsITS"}; + Configurable cfgMaxDCAxy{"cfgMaxDCAxy", 0.05, "max dca XY for single track in cm"}; + Configurable cfgMaxDCAz{"cfgMaxDCAz", 0.05, "max dca Z for single track in cm"}; + Configurable cfgMaxDCA3dsigmaTrack{"cfgMaxDCA3dsigmaTrack", 1.5, "max DCA 3D in sigma"}; + Configurable cfgMaxFracSharedClustersTPC{"cfgMaxFracSharedClustersTPC", 999.f, "max fraction of shared clusters in TPC"}; + Configurable cfgApplyCutsFromPrefilterDerived{"cfgApplyCutsFromPrefilterDerived", false, "flag to apply prefilter to electron"}; + Configurable includeITSsa{"includeITSsa", false, "Flag to enable ITSsa tracks"}; + Configurable cfgMaxpttrackITSsa{"cfgMaxpttrackITSsa", 0.15, "max pt for ITSsa tracks"}; - std::vector fTagCuts; - std::vector fProbeCuts; - std::vector fPairCuts; + Configurable cfgPIDscheme{"cfgPIDscheme", static_cast(DalitzEECut::PIDSchemes::kTOFif), "pid scheme [kTOFif : 0, kTPConly : 1]"}; + Configurable cfgMinTPCNsigmaEl{"cfgMinTPCNsigmaEl", -2.0, "min. TPC n sigma for electron inclusion"}; + Configurable cfgMaxTPCNsigmaEl{"cfgMaxTPCNsigmaEl", +3.0, "max. TPC n sigma for electron inclusion"}; + Configurable cfgMinTPCNsigmaPi{"cfgMinTPCNsigmaPi", -0.0, "min. TPC n sigma for pion exclusion"}; + Configurable cfgMaxTPCNsigmaPi{"cfgMaxTPCNsigmaPi", +0.0, "max. TPC n sigma for pion exclusion"}; + Configurable cfgMinTOFNsigmaEl{"cfgMinTOFNsigmaEl", -3.0, "min. TOF n sigma for electron inclusion"}; + Configurable cfgMaxTOFNsigmaEl{"cfgMaxTOFNsigmaEl", +3.0, "max. TOF n sigma for electron inclusion"}; + } dileptoncuts; - std::vector fPairNames; - void init(InitContext& context) + V0PhotonCut fV0PhotonCut; + struct : ConfigurableGroup { + std::string prefix = "pcmcut_group"; + Configurable cfgRequireV0WithITSTPC{"cfgRequireV0WithITSTPC", false, "flag to select V0s with ITS-TPC matched tracks"}; + Configurable cfgRequireV0WithITSOnly{"cfgRequireV0WithITSOnly", false, "flag to select V0s with ITSonly tracks"}; + Configurable cfgRequireV0WithTPCOnly{"cfgRequireV0WithTPCOnly", false, "flag to select V0s with TPConly tracks"}; + Configurable cfgMinPtV0{"cfgMinPtV0", 0.1, "min pT for v0 photons at PV"}; + Configurable cfgMaxPtV0{"cfgMaxPtV0", 1e+10, "max pT for v0 photons at PV"}; + Configurable cfgMinEtaV0{"cfgMinEtaV0", -0.8, "min eta for v0 photons at PV"}; + Configurable cfgMaxEtaV0{"cfgMaxEtaV0", +0.8, "max eta for v0 photons at PV"}; + Configurable cfgMinV0Radius{"cfgMinV0Radius", 4.0, "min v0 radius"}; + Configurable cfgMaxV0Radius{"cfgMaxV0Radius", 90.0, "max v0 radius"}; + Configurable cfgMaxAlphaAp{"cfgMaxAlphaAp", 0.95, "max alpha for AP cut"}; + Configurable cfgMaxQtAp{"cfgMaxQtAp", 0.01, "max qT for AP cut"}; + Configurable cfgMinCospa{"cfgMinCospa", 0.999, "min V0 CosPA"}; + Configurable cfgMaxPca{"cfgMaxPca", 1.5, "max distance btween 2 legs"}; + Configurable cfgMaxChi2kf{"cfgMaxChi2kf", 1e+10, "max chi2/ndf with KF"}; + Configurable cfgRejectV0OnITSib{"cfgRejectV0OnITSib", true, "flag to reject V0s on ITSib"}; + Configurable cfgMinNclusterTPC{"cfgMinNclusterTPC", 0, "min ncluster tpc"}; + Configurable cfgMinNcrossedrows{"cfgMinNcrossedrows", 40, "min ncrossed rows"}; + Configurable cfgMaxFracSharedClustersTPC{"cfgMaxFracSharedClustersTPC", 999.f, "max fraction of shared clusters in TPC"}; + Configurable cfgMaxChi2TPC{"cfgMaxChi2TPC", 4.0, "max chi2/NclsTPC"}; + Configurable cfgMaxChi2ITS{"cfgMaxChi2ITS", 36.0, "max chi2/NclsITS"}; + Configurable cfgMinTPCNsigmaEl{"cfgMinTPCNsigmaEl", -3.0, "min. TPC n sigma for electron"}; + Configurable cfgMaxTPCNsigmaEl{"cfgMaxTPCNsigmaEl", +3.0, "max. TPC n sigma for electron"}; + Configurable cfgDisableITSonlytrack{"cfgDisableITSonlytrack", false, "flag to disable ITSonly tracks"}; + Configurable cfgDisableTPCOnlytrack{"cfgDisableTPCOnlytrack", false, "flag to disable TPConly tracks"}; + } pcmcuts; + + void init(InitContext&) { - if (context.mOptions.get("processMB")) { - fPairNames.push_back("PCMPCM"); - } - DefineTagCuts(); - DefineProbeCuts(); - DefinePairCuts(); + defineEMEventCut(); + defineDileptonCut(); + definePCMCut(); addhistograms(); - TString ev_cut_name = fConfigEMEventCut.value; - fEMEventCut = *eventcuts::GetCut(ev_cut_name.Data()); + } - fOutputEvent.setObject(reinterpret_cast(fMainList->FindObject("Event"))); - fOutputV0.setObject(reinterpret_cast(fMainList->FindObject("V0"))); - fOutputPair.setObject(reinterpret_cast(fMainList->FindObject("Pair"))); + 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.SetRequireNoCollInTimeRangeStandard(eventcuts.cfgRequireNoCollInTimeRangeStandard); + fEMEventCut.SetRequireNoCollInTimeRangeStrict(eventcuts.cfgRequireNoCollInTimeRangeStrict); + fEMEventCut.SetRequireNoCollInITSROFStandard(eventcuts.cfgRequireNoCollInITSROFStandard); + fEMEventCut.SetRequireNoCollInITSROFStrict(eventcuts.cfgRequireNoCollInITSROFStrict); + fEMEventCut.SetRequireNoHighMultCollInPrevRof(eventcuts.cfgRequireNoHighMultCollInPrevRof); } - template - void add_pair_histograms(THashList* list_pair, const std::string pairname, TCuts1 const& tagcuts, TCuts2 const& probecuts, TCuts3 const& cuts3) + void defineDileptonCut() { - for (auto& cut1 : tagcuts) { - for (auto& cut2 : probecuts) { - std::string cutname1 = cut1.GetName(); - std::string cutname2 = cut2.GetName(); - - THashList* list_pair_subsys = reinterpret_cast(list_pair->FindObject(pairname.data())); - std::string photon_cut_name = cutname1 + "_" + cutname2; - o2::aod::pwgem::photon::histogram::AddHistClass(list_pair_subsys, photon_cut_name.data()); - THashList* list_pair_subsys_photoncut = reinterpret_cast(list_pair_subsys->FindObject(photon_cut_name.data())); - - for (auto& cut3 : cuts3) { - std::string pair_cut_name = cut3.GetName(); - o2::aod::pwgem::photon::histogram::AddHistClass(list_pair_subsys_photoncut, pair_cut_name.data()); - THashList* list_pair_subsys_paircut = reinterpret_cast(list_pair_subsys_photoncut->FindObject(pair_cut_name.data())); - o2::aod::pwgem::photon::histogram::DefineHistograms(list_pair_subsys_paircut, "material_budget_study", "Pair"); - } // end of cut3 loop pair cut - } // end of cut2 loop - } // end of cut1 loop + fDileptonCut = DalitzEECut("fDileptonCut", "fDileptonCut"); + + // for pair + fDileptonCut.SetMeeRange(dileptoncuts.cfgMinMass, dileptoncuts.cfgMaxMass); + fDileptonCut.SetMaxPhivPairMeeDep([&](float mll) { return (mll - dileptoncuts.cfgPhivIntercept) / dileptoncuts.cfgPhivSlope; }); + fDileptonCut.ApplyPhiV(dileptoncuts.cfgApplyPhiv); + fDileptonCut.RequireITSibAny(dileptoncuts.cfgRequireItsibAny); + fDileptonCut.RequireITSib1st(dileptoncuts.cfgRequireItsib1st); + + // for tracks + fDileptonCut.SetTrackPtRange(dileptoncuts.cfgMinPtTrack, 1e+10f); + fDileptonCut.SetTrackEtaRange(-dileptoncuts.cfgMaxEtaTrack, +dileptoncuts.cfgMaxEtaTrack); + fDileptonCut.SetMinNClustersTPC(dileptoncuts.cfgMinnclusterTPC); + fDileptonCut.SetMinNCrossedRowsTPC(dileptoncuts.cfgMinncrossedrows); + fDileptonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fDileptonCut.SetMaxFracSharedClustersTPC(dileptoncuts.cfgMaxFracSharedClustersTPC); + fDileptonCut.SetChi2PerClusterTPC(0.0, dileptoncuts.cfgMaxchi2TPC); + fDileptonCut.SetChi2PerClusterITS(0.0, dileptoncuts.cfgMaxchi2ITS); + fDileptonCut.SetNClustersITS(dileptoncuts.cfgMinnclusterITS, 7); + fDileptonCut.SetMaxDcaXY(dileptoncuts.cfgMaxDCAxy); + fDileptonCut.SetMaxDcaZ(dileptoncuts.cfgMaxDCAz); + fDileptonCut.SetTrackDca3DRange(0.f, dileptoncuts.cfgMaxDCA3dsigmaTrack); // in sigma + fDileptonCut.IncludeITSsa(dileptoncuts.includeITSsa, dileptoncuts.cfgMaxpttrackITSsa); + + // for eID + fDileptonCut.SetPIDScheme(dileptoncuts.cfgPIDscheme); + fDileptonCut.SetTPCNsigmaElRange(dileptoncuts.cfgMinTPCNsigmaEl, dileptoncuts.cfgMaxTPCNsigmaEl); + fDileptonCut.SetTPCNsigmaPiRange(dileptoncuts.cfgMinTPCNsigmaPi, dileptoncuts.cfgMaxTPCNsigmaPi); + fDileptonCut.SetTOFNsigmaElRange(dileptoncuts.cfgMinTOFNsigmaEl, dileptoncuts.cfgMaxTOFNsigmaEl); + } + + void definePCMCut() + { + fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); + + // for v0 + fV0PhotonCut.SetV0PtRange(pcmcuts.cfgMinPtV0, pcmcuts.cfgMaxPtV0); + fV0PhotonCut.SetV0EtaRange(pcmcuts.cfgMinEtaV0, pcmcuts.cfgMaxEtaV0); + fV0PhotonCut.SetMinCosPA(pcmcuts.cfgMinCospa); + fV0PhotonCut.SetMaxPCA(pcmcuts.cfgMaxPca); + fV0PhotonCut.SetMaxChi2KF(pcmcuts.cfgMaxChi2kf); + fV0PhotonCut.SetRxyRange(pcmcuts.cfgMinV0Radius, pcmcuts.cfgMaxV0Radius); + fV0PhotonCut.SetAPRange(pcmcuts.cfgMaxAlphaAp, pcmcuts.cfgMaxQtAp); + fV0PhotonCut.RejectITSib(pcmcuts.cfgRejectV0OnITSib); + + // for track + fV0PhotonCut.SetMinNClustersTPC(pcmcuts.cfgMinNclusterTPC); + fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfgMinNcrossedrows); + fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.cfgMaxFracSharedClustersTPC); + fV0PhotonCut.SetChi2PerClusterTPC(0.0, pcmcuts.cfgMaxChi2TPC); + fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.cfgMinTPCNsigmaEl, pcmcuts.cfgMaxTPCNsigmaEl); + fV0PhotonCut.SetChi2PerClusterITS(-1e+10, pcmcuts.cfgMaxChi2ITS); + fV0PhotonCut.SetNClustersITS(0, 7); + fV0PhotonCut.SetMeanClusterSizeITSob(0.0, 16.0); + fV0PhotonCut.SetDisableITSonly(pcmcuts.cfgDisableITSonlytrack); + fV0PhotonCut.SetDisableTPConly(pcmcuts.cfgDisableTPCOnlytrack); + fV0PhotonCut.SetRequireITSTPC(pcmcuts.cfgRequireV0WithITSTPC); + fV0PhotonCut.SetRequireITSonly(pcmcuts.cfgRequireV0WithITSOnly); + fV0PhotonCut.SetRequireTPConly(pcmcuts.cfgRequireV0WithTPCOnly); } - static constexpr std::string_view pairnames[9] = {"PCMPCM", "PHOSPHOS", "EMCEMC", "PCMPHOS", "PCMEMC", "PCMDalitzEE", "PCMDalitzMuMu", "PHOSEMC", "DalitzEEDalitzEE"}; void addhistograms() { - fMainList->SetOwner(true); - fMainList->SetName("fMainList"); - // create sub lists first. - o2::aod::pwgem::photon::histogram::AddHistClass(fMainList, "Event"); - THashList* list_ev = reinterpret_cast(fMainList->FindObject("Event")); + auto hCollisionCounter = registry.add("Event/before/hCollisionCounter", "collision counter;;Number of events", kTH1F, {{10, 0.5, 10.5}}, false); + hCollisionCounter->GetXaxis()->SetBinLabel(1, "all"); + hCollisionCounter->GetXaxis()->SetBinLabel(2, "No TF border"); + hCollisionCounter->GetXaxis()->SetBinLabel(3, "No ITS ROF border"); + hCollisionCounter->GetXaxis()->SetBinLabel(4, "No Same Bunch Pileup"); + hCollisionCounter->GetXaxis()->SetBinLabel(5, "Is Vertex ITSTPC"); + hCollisionCounter->GetXaxis()->SetBinLabel(6, "Is Good Zvtx FT0vsPV"); + hCollisionCounter->GetXaxis()->SetBinLabel(7, "FT0AND"); + hCollisionCounter->GetXaxis()->SetBinLabel(8, "sel8"); + hCollisionCounter->GetXaxis()->SetBinLabel(9, "|Z_{vtx}| < 10 cm"); + hCollisionCounter->GetXaxis()->SetBinLabel(10, "accepted"); + + registry.add("Event/before/hZvtx", "vertex z; Z_{vtx} (cm)", kTH1F, {{100, -50, +50}}, false); + registry.add("Event/before/hMultNTracksPV", "hMultNTracksPV; N_{track} to PV", kTH1F, {{6001, -0.5, 6000.5}}, false); + registry.add("Event/before/hMultNTracksPVeta1", "hMultNTracksPVeta1; N_{track} to PV", kTH1F, {{6001, -0.5, 6000.5}}, false); + registry.add("Event/before/hMultFT0", "hMultFT0;mult. FT0A;mult. FT0C", kTH2F, {{300, 0, 6000}, {300, 0, 6000}}, false); + registry.add("Event/before/hCentFT0A", "hCentFT0A;centrality FT0A (%)", kTH1F, {{110, 0, 110}}, false); + registry.add("Event/before/hCentFT0C", "hCentFT0C;centrality FT0C (%)", kTH1F, {{110, 0, 110}}, false); + registry.add("Event/before/hCentFT0M", "hCentFT0M;centrality FT0M (%)", kTH1F, {{110, 0, 110}}, false); + registry.add("Event/before/hCentFT0MvsMultNTracksPV", "hCentFT0MvsMultNTracksPV;centrality FT0M (%);N_{track} to PV", kTH2F, {{110, 0, 110}, {600, 0, 6000}}, false); + registry.add("Event/before/hMultFT0MvsMultNTracksPV", "hMultFT0MvsMultNTracksPV;mult. FT0M;N_{track} to PV", kTH2F, {{600, 0, 6000}, {600, 0, 6000}}, false); + registry.addClone("Event/before/", "Event/after/"); + + if (cfgPlotMBGeneral) { + + registry.add("MBGeneral", ";z_{conv} (cm); R_{conv} (cm); #eta; #varphi (rad); p_{T}", kTHnSparseF, + {{200, -100.f, 100.f}, + {100, 0.f, 100.f}, + {80, -2, +2}, + {90, 0, o2::constants::math::TwoPI}, + {1000, 0, 10}}, // pT 5 + true); - o2::aod::pwgem::photon::histogram::AddHistClass(fMainList, "Pair"); - THashList* list_pair = reinterpret_cast(fMainList->FindObject("Pair")); + registry.add("Pi0/Same", + "Pi0; m_{#gamma#gamma}; #it{p}_{T}; R_{conv} (cm); R_{conv} (cm)", + kTHnSparseF, + {{400, 0.f, 0.3}, // x 0 + {500, 0, 10}, // pT 1 + {100, 0.f, 100.f}, // R of photon 1 + {100, 0.f, 100.f}}, + true); - o2::aod::pwgem::photon::histogram::AddHistClass(fMainList, "V0"); - THashList* list_v0 = reinterpret_cast(fMainList->FindObject("V0")); + registry.add("Dalitz/Same", + "Dalitz; m_{eeγ}; p_{T}; R_{γ} (cm); R_{e+}; R_{e-}", + kTHnSparseF, + {{400, 0.f, 0.5f}, + {500, 0.f, 10.f}, + {100, 0.f, 100.f}}, + true); - // for V0s - for (const auto& cut : fProbeCuts) { - const char* cutname = cut.GetName(); - THashList* list_v0_cut = o2::aod::pwgem::photon::histogram::AddHistClass(list_v0, cutname); - o2::aod::pwgem::photon::histogram::DefineHistograms(list_v0_cut, "material_budget_study", "V0"); + registry.add("Pi0/Mix", "Pi0 mixed-event;m_{#gamma#gamma};p_{T}", kTHnSparseF, + {{400, 0, 0.3}, {500, 0, 10}, {100, 0, 100}, {100, 0, 100}}, true); } - for (auto& pairname : fPairNames) { - LOGF(info, "Enabled pairs = %s", pairname.data()); + if (cfgPlotMBDetailed) { - THashList* list_ev_pair = reinterpret_cast(o2::aod::pwgem::photon::histogram::AddHistClass(list_ev, pairname.data())); - for (const auto& evtype : event_types) { - THashList* list_ev_type = reinterpret_cast(o2::aod::pwgem::photon::histogram::AddHistClass(list_ev_pair, evtype.data())); - o2::aod::pwgem::photon::histogram::DefineHistograms(list_ev_type, "Event", evtype.data()); - } + registry.add("MBStudiesWireLeft", + "photon characteristics for MB studies; x (cm); y (cm); z (cm); #varphi (rad.); R_{conv} (cm)", + kTHnSparseF, + { + {200, -20.f, 20.f}, // x 0 + {200, -20.f, 20.f}, // y 1 + {80, -20.f, 20.f}, // z 2 + {40, 3.15, 3.4}, // phi 3 + {80, 0.f, 20.f}, // conversion radius 4 + {1000, 0, 10} // pT 5 + }, + true); + + registry.add("MBStudiesWireRight", + "photon characteristics for MB studies; x (cm); y (cm); z (cm); #varphi (rad.); R_{conv} (cm)", + kTHnSparseF, + { + {200, -20.f, 20.f}, // x 0 + {200, -20.f, 20.f}, // y 1 + {80, -20.f, 20.f}, // z 2 + {40, 6.00, 6.15}, // phi 3 + {80, 0.f, 20.f}, // conversion radius 4 + {1000, 0, 10} // pT 5 + }, + true); + + registry.add("MBStudiesWireITS", + "photon characteristics for MB studies; x (cm); y (cm); z (cm); #varphi (rad.); R_{conv} (cm)", + kTHnSparseF, + { + {160, -40.f, 40.f}, // x 0 + {160, -40.f, 40.f}, // y 1 + {80, -40.f, 40.f}, // z 2 + {90, 0, TwoPI}, // phi 3 + {150, 0.f, 60.f}, // conversion radius 4 + {1000, 0, 10} // pT 5 + }, + true); + + registry.add("MBStudiesMFT", + "photon characteristics for MB studies; x (cm); y (cm); z (cm); #varphi (rad.); R_{conv} (cm)", + kTHnSparseF, + { + {160, -80.f, 80.f}, // x 0 + {160, -80.f, 80.f}, // y 1 + {40, -40.f, 40.f}, // z 2 + {90, 0, TwoPI}, // phi 3 + {40, 40.f, 60.f}, // conversion radius 4 + {500, 0, 10} // pT 5 + }, + true); - o2::aod::pwgem::photon::histogram::AddHistClass(list_pair, pairname.data()); + registry.add("MBStudiesTPC", + "photon characteristics for MB studies; x (cm); y (cm); z (cm); #varphi (rad.); R_{conv} (cm)", + kTHnSparseF, + { + {160, -90.f, 90.f}, // x 0 + {160, -90.f, 90.f}, // y 1 + {80, -40.f, 40.f}, // z 2 + {90, 0, TwoPI}, // phi 3 + {80, 60.f, 80.f}, // conversion radius 4 + {1000, 0, 10} // pT 5 + }, + true); - if (pairname == "PCMPCM") { - add_pair_histograms(list_pair, pairname, fTagCuts, fProbeCuts, fPairCuts); + for (int nCls = 0; nCls <= 7; nCls++) { // o2-linter: disable=magic-number (just numbers for ITS cluster) + registry.add(Form("ITSHits/Neg/ITSCls%d", nCls), + Form("Conversion radius NEG legs with %d ITS clusters;R_{conv} (cm);Entries", nCls), + HistType::kTH1F, + {AxisSpec{100, 0.f, 100.f}}); + + registry.add(Form("ITSHits/Pos/ITSCls%d", nCls), + Form("Conversion radius POS legs with %d ITS clusters;R_{conv} (cm);Entries", nCls), + HistType::kTH1F, + {AxisSpec{100, 0.f, 100.f}}); } + } + + if (cfgPlotResolution) { + + registry.add("ResolutionGen/Z_res", + "Z resolution (gen ref);z_gen (cm);R_gen (cm);phi_gen;eta_gen;pT_gen;Δz (cm)", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {120, -30, 30}}, false); + + registry.add("ResolutionGen/R_res", + "R resolution (gen ref);z_gen (cm);R_gen (cm);phi_gen;eta_gen;pT_gen;ΔR (cm)", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {120, -30, 30}}, false); + + registry.add("ResolutionGen/Phi_res", + "Phi resolution (gen ref);z_gen (cm);R_gen (cm);phi_gen;eta_gen;pT_gen;Δφ", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {100, -0.2f, 0.2f}}, false); + + registry.add("ResolutionGen/Pt_res", + "pT resolution (gen ref);z_gen (cm);R_gen (cm);phi_gen;eta_gen;pT_gen;ΔpT/pT_gen", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {200, -1.0f, 1.0f}}, false); + + registry.add("ResolutionGen/Eta_res", + "Eta resolution (gen ref);z_gen (cm);R_gen (cm);phi_gen;eta_gen;pT_gen;Δη", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {100, -0.5f, 0.5f}}, false); + + registry.add("ResolutionRec/Z_res", + "Z resolution (rec ref);z_rec (cm);R_rec (cm);phi_rec;eta_rec;pT_rec;Δz (cm)", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {120, -30, 30}}, false); + + registry.add("ResolutionRec/R_res", + "R resolution (rec ref);z_rec (cm);R_rec (cm);phi_rec;eta_rec;pT_rec;ΔR (cm)", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {120, -30, 30}}, false); + + registry.add("ResolutionRec/Phi_res", + "Phi resolution (rec ref);z_rec (cm);R_rec (cm);phi_rec;eta_rec;pT_rec;Δφ", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {100, -0.2f, 0.2f}}, false); + + registry.add("ResolutionRec/Pt_res", + "pT resolution (rec ref);z_rec (cm);R_rec (cm);phi_rec;eta_rec;pT_rec;ΔpT/pT_gen", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {200, -1.0f, 1.0f}}, false); + + registry.add("ResolutionRec/Eta_res", + "Eta resolution (rec ref);z_rec (cm);R_rec (cm);phi_rec;eta_rec;pT_rec;Δη", + kTHnSparseF, + {{200, -100, 100}, {200, 0, 100}, {90, 0, TwoPI}, {200, -1.0f, 1.0f}, {500, 0, 10}, {100, -0.5f, 0.5f}}, false); + } + + if (cfgPlotPurity) { + + registry.add("PhotonPurity", + "Photon purity classification;p_{T} (GeV/c);R_{conv} (cm);#eta;#varphi;category", + kTHnSparseF, + { + {500, 0, 10}, // pT + {200, 0, 100}, // R + {100, -1.0, 1.0}, // eta + {90, 0, o2::constants::math::TwoPI}, // phi + {12, 0.5, 12.5} // purity categories + }, + true); + + auto hPurity = registry.get(HIST("PhotonPurity")); - } // end of pair name loop + auto axis = hPurity->GetAxis(4); + + axis->SetBinLabel(1, "e^{+}e^{-}"); + axis->SetBinLabel(2, "e^{-} + #pi^{+}"); + axis->SetBinLabel(3, "#pi^{-} + e^{+}"); + axis->SetBinLabel(4, "K^{-} + e^{+}"); + axis->SetBinLabel(5, "e^{-} + K^{+}"); + axis->SetBinLabel(6, "#bar{p} + p"); + axis->SetBinLabel(7, "K^{-} + K^{+}"); + axis->SetBinLabel(8, "#pi^{-} + p"); + axis->SetBinLabel(9, "#bar{p} + #pi^{+}"); + axis->SetBinLabel(10, "#mu^{-} + #mu^{+}"); + } + + if (cfgPlotBremsstrahlung) { + + registry.add("Bremsstrahlung/EnergyLossVsR", ";R (cm);E_{loss} (GeV)", + kTH2F, {{50, 0, 100}, {100, 0, 5}}); + registry.add("Bremsstrahlung/EnergyLossXY", ";X (cm);Y (cm)", + kTH2F, {{100, -100, 100}, {100, -100, 100}}); + registry.add("Bremsstrahlung/EnergyLossXYWeigh", ";X (cm);Y (cm)", + kTH2F, {{100, -100, 100}, {100, -100, 100}}); + registry.add("Bremsstrahlung/EnergyLossVsZ", ";Z (cm);E_{loss} (GeV)", + kTH2F, {{100, -100, 100}, {200, 0, 5}}); + + registry.add("Bremsstrahlung/NBrem", ";N_{Brem} per e^{-};Counts", + kTH1I, {{10, 0, 10}}); + + registry.add("Bremsstrahlung/DeltaPoverPvsR", ";R (cm);(p_{Reco}-p_{MC})/p_{MC}", + kTH2F, {{50, 0, 100}, {200, -1, 1}}); + registry.add("Bremsstrahlung/DeltaPtvsR", ";R (cm);(p_{T,Reco}-p_{T,MC})/p_{T,MC}", + kTH2F, {{50, 0, 100}, {300, -2, 2}}); + registry.add("Bremsstrahlung/DeltaPhivsR", ";R (cm);#Delta#varphi (rad)", + kTH2F, {{50, 0, 100}, {200, -0.1, 0.1}}); + registry.add("Bremsstrahlung/DeltaEtavsR", ";R (cm);#Delta#eta", + kTH2F, {{50, 0, 100}, {200, -0.1, 0.1}}); + + registry.add("Bremsstrahlung/Sigma1PtVsR", ";R (cm);#sigma(1/p_{T}) [c/GeV]", + kTH2F, {{50, 0, 100}, {100, 0, 0.1}}); + + registry.add("Bremsstrahlung/SigmaYVsR", ";R (cm);#sigma(y) [cm]", + kTH2F, {{50, 0, 100}, {100, 0, 0.5}}); + registry.add("Bremsstrahlung/SigmaZVsR", ";R (cm);#sigma(z) [cm]", + kTH2F, {{50, 0, 100}, {100, 0, 1.0}}); + + registry.add("Bremsstrahlung/relativeResoPtWBrems", + "relative #it{p}_{T} resolution; #it{p}_{T}; #sigma(#it{p}_{T})/#it{p}_{T}", + kTProfile, + {axisPt}); + + registry.add("Bremsstrahlung/relativeResoPtWOBrems", + "relative #it{p}_{T} resolution; #it{p}_{T}; #sigma(#it{p}_{T})/#it{p}_{T}", + kTProfile, + {axisPt}); + } + + if (cfgPlotMBCollisions) { + + registry.add("RightCollisions/SparseDeltaCol", + "RightCollisions;R_{rec} (cm);#it{p}_{T} (GeV/c);#Delta #it{p}_{T}/#it{p}_{T,gen};#Delta z (cm);#Delta R (cm);#Delta col ID", + kTHnSparseF, + { + {100, 0.f, 100.f}, // R_rec + {200, 0.f, 10.f}, // pT + {200, -1.f, 1.f}, // ΔpT/pT + {200, -20.f, 20.f}, // Δz + {200, -8.f, 8.f}, // ΔR + {21, -10.5, 10.5} // ΔCollisionID + }, + true); + + registry.add("WrongCollisions/SparseDeltaCol", + "WrongCollisions;R_{rec} (cm);#it{p}_{T} (GeV/c);#Delta #it{p}_{T}/#it{p}_{T,gen};#Delta z (cm);#Delta R (cm);#Delta col ID", + kTHnSparseF, + {{100, 0.f, 100.f}, + {200, 0.f, 10.f}, + {200, -1.f, 1.f}, + {200, -20.f, 20.f}, + {200, -8.f, 8.f}, + {21, -10.5, 10.5}}, + true); + } } - void DefineTagCuts() + template + void fillITSClsNeg(HistogramRegistry& reg, float value) { - TString cutNamesStr = fConfigTagCuts.value; - if (!cutNamesStr.IsNull()) { - std::unique_ptr objArray(cutNamesStr.Tokenize(",")); - for (int icut = 0; icut < objArray->GetEntries(); ++icut) { - const char* cutname = objArray->At(icut)->GetName(); - LOGF(info, "add cut : %s", cutname); - fTagCuts.push_back(*pcmcuts::GetCut(cutname)); + if constexpr (N == 0) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls0"), value); + if constexpr (N == 1) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls1"), value); + if constexpr (N == 2) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls2"), value); + if constexpr (N == 3) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls3"), value); + if constexpr (N == 4) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls4"), value); + if constexpr (N == 5) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls5"), value); + if constexpr (N == 6) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls6"), value); + if constexpr (N == 7) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Neg/ITSCls7"), value); + } + + template + void fillITSClsPos(HistogramRegistry& reg, float value) + { + if constexpr (N == 0) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls0"), value); + if constexpr (N == 1) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls1"), value); + if constexpr (N == 2) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls2"), value); + if constexpr (N == 3) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls3"), value); + if constexpr (N == 4) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls4"), value); + if constexpr (N == 5) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls5"), value); + if constexpr (N == 6) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls6"), value); + if constexpr (N == 7) // o2-linter: disable=magic-number (just numbers for ITS cluster) + reg.fill(HIST("ITSHits/Pos/ITSCls7"), value); + } + + template + void reconstructMesons(const TV0s& v0s) + { + for (auto i = 0; i < v0s.size(); i++) { + auto g1 = v0s.iteratorAt(i); + + for (auto j = i + 1; j < v0s.size(); j++) { + auto g2 = v0s.iteratorAt(j); + + ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + registry.fill(HIST("Pi0/Same"), v12.M(), v12.Pt(), + g1.v0radius(), g2.v0radius()); } } - LOGF(info, "Number of Tag cuts = %d", fTagCuts.size()); } - void DefineProbeCuts() + void fillClusterNeg(HistogramRegistry& reg, int nCls, float value) { - TString cutNamesStr = fConfigProbeCuts.value; - if (!cutNamesStr.IsNull()) { - std::unique_ptr objArray(cutNamesStr.Tokenize(",")); - for (int icut = 0; icut < objArray->GetEntries(); ++icut) { - const char* cutname = objArray->At(icut)->GetName(); - LOGF(info, "add cut : %s", cutname); - fProbeCuts.push_back(*pcmcuts::GetCut(cutname)); - } + switch (nCls) { + case 0: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<0>(reg, value); + break; + case 1: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<1>(reg, value); + break; + case 2: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<2>(reg, value); + break; + case 3: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<3>(reg, value); + break; + case 4: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<4>(reg, value); + break; + case 5: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<5>(reg, value); + break; + case 6: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<6>(reg, value); + break; + case 7: // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillITSClsNeg<7>(reg, value); + break; + default: + break; } - LOGF(info, "Number of Probe cuts = %d", fProbeCuts.size()); } - void DefinePairCuts() + void fillClusterPos(HistogramRegistry& reg, int nCls, float value) { - TString cutNamesStr = fConfigPairCuts.value; - if (!cutNamesStr.IsNull()) { - std::unique_ptr objArray(cutNamesStr.Tokenize(",")); - for (int icut = 0; icut < objArray->GetEntries(); ++icut) { - const char* cutname = objArray->At(icut)->GetName(); - LOGF(info, "add cut : %s", cutname); - fPairCuts.push_back(*paircuts::GetCut(cutname)); - } + switch (nCls) { + case 0: + fillITSClsPos<0>(reg, value); + break; + case 1: + fillITSClsPos<1>(reg, value); + break; + case 2: + fillITSClsPos<2>(reg, value); + break; + case 3: + fillITSClsPos<3>(reg, value); + break; + case 4: + fillITSClsPos<4>(reg, value); + break; + case 5: + fillITSClsPos<5>(reg, value); + break; + case 6: + fillITSClsPos<6>(reg, value); + break; + case 7: + fillITSClsPos<7>(reg, value); + break; + default: + break; + } + } + template + void fillEventInfo(TCollision const& collision, const float /*weight*/ = 1.f) + { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 1.0); + if (collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 2.0); } - LOGF(info, "Number of Pair cuts = %d", fPairCuts.size()); + if (collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 3.0); + } + if (collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 4.0); + } + if (collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 5.0); + } + if (collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 6.0); + } + if (collision.selection_bit(o2::aod::evsel::kIsTriggerTVX)) { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 7.0); + } + if (collision.sel8()) { + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 8.0); + } + if (std::fabs(collision.posZ()) < 10.0) { // o2-linter: disable=magic-number (vertex position) + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCollisionCounter"), 9.0); + } + + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hZvtx"), collision.posZ()); + + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hMultNTracksPV"), collision.multNTracksPV()); + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hMultNTracksPVeta1"), collision.multNTracksPVeta1()); + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hMultFT0"), collision.multFT0A(), collision.multFT0C()); + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCentFT0A"), collision.centFT0A()); + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCentFT0C"), collision.centFT0C()); + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCentFT0M"), collision.centFT0M()); + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hCentFT0MvsMultNTracksPV"), collision.centFT0M(), collision.multNTracksPV()); + registry.fill(HIST("Event/") + HIST(EventTypes[evID]) + HIST("hMultFT0MvsMultNTracksPV"), collision.multFT0A() + collision.multFT0C(), collision.multNTracksPV()); } + template + void reconstructMesonsMixed(const TV0s& current, + const std::deque>& pool) + { + for (const auto& prev : pool) { + for (const auto& g1 : current) { + if (!fV0PhotonCut.IsSelected(g1)) { + continue; + } + for (const auto& g2 : prev) { + ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(g2.pt, g2.eta, g2.phi, 0.); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - Preslice perCollision_pcm = aod::v0photonkf::emeventId; + registry.fill(HIST("Pi0/Mix"), v12.M(), v12.Pt(), + g1.v0radius(), g2.r); + } + } + } + } - template - bool IsSelectedPair(TG1 const& g1, TG2 const& g2, TCut1 const& cut1, TCut2 const& cut2) + template + void reconstructDalitz(const TV0s& v0s, + const TPositrons& positronsPerCollision, + const TElectrons& electronsPerCollision, + MyCollision const& collision) { - return o2::aod::pwgem::photonmeson::photonpair::IsSelectedPair(g1, g2, cut1, cut2); + + for (const auto& g1 : v0s) { + if (!fV0PhotonCut.IsSelected(g1)) { + continue; + } + ROOT::Math::PtEtaPhiMVector vGamma(g1.pt(), g1.eta(), g1.phi(), 0.); + + for (const auto& [pos, ele] : + combinations(CombinationsFullIndexPolicy(positronsPerCollision, + electronsPerCollision))) { + + if (pos.trackId() == ele.trackId()) { + continue; + } + + if (!fDileptonCut.template IsSelectedTrack(pos, collision) || + !fDileptonCut.template IsSelectedTrack(ele, collision)) { + continue; // Track-Cuts + } + + ROOT::Math::PtEtaPhiMVector vPos(pos.pt(), pos.eta(), pos.phi(), + o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector vEle(ele.pt(), ele.eta(), ele.phi(), + o2::constants::physics::MassElectron); + + auto vee = vPos + vEle; + auto veGamma = vGamma + vee; + + registry.fill(HIST("Dalitz/Same"), + veGamma.M(), vee.Pt(), + g1.v0radius()); + } + } } - template - void fillsinglephoton(TEvents const& collisions, TPhotons const& photons, TPreslice const& perCollision, TCuts const& cuts, TLegs const& /*legs*/) + SliceCache cache; + Preslice perCollision = aod::v0photonkf::emeventId; + Preslice perCollisionElectron = aod::emprimaryelectron::emeventId; + + Partition positrons = o2::aod::emprimaryelectron::sign > int8_t(0) && dileptoncuts.cfgMinPtTrack < o2::aod::track::pt&& nabs(o2::aod::track::eta) < dileptoncuts.cfgMaxEtaTrack; + Partition electrons = o2::aod::emprimaryelectron::sign < int8_t(0) && dileptoncuts.cfgMinPtTrack < o2::aod::track::pt && nabs(o2::aod::track::eta) < dileptoncuts.cfgMaxEtaTrack; + + Filter prefilterPrimaryelectron = ifnode(dileptoncuts.cfgApplyCutsFromPrefilterDerived.node(), o2::aod::emprimaryelectron::pfbderived == static_cast(0), true); + + void processRec(MyCollisions const& collisions, + MyV0Photons const& v0photons, + aod::V0Legs const&, + MyPrimaryElectrons const& electrons) { - THashList* list_ev_pair_before = static_cast(fMainList->FindObject("Event")->FindObject(pairnames[pairtype].data())->FindObject(event_types[0].data())); - THashList* list_ev_pair_after = static_cast(fMainList->FindObject("Event")->FindObject(pairnames[pairtype].data())->FindObject(event_types[1].data())); - THashList* list_v0 = static_cast(fMainList->FindObject("V0")); - double value[4] = {0.f}; - for (auto& collision : collisions) { + for (const auto& collision : collisions) { const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { continue; } - o2::aod::pwgem::photon::histogram::FillHistClass(list_ev_pair_before, "", collision); + fillEventInfo<0>(collision); if (!fEMEventCut.IsSelected(collision)) { continue; } - o2::aod::pwgem::photon::histogram::FillHistClass(list_ev_pair_after, "", collision); - reinterpret_cast(list_ev_pair_before->FindObject("hCollisionCounter"))->Fill("accepted", 1.f); - reinterpret_cast(list_ev_pair_after->FindObject("hCollisionCounter"))->Fill("accepted", 1.f); - - auto photons_coll = photons.sliceBy(perCollision, collision.globalIndex()); - for (auto& cut : cuts) { - for (auto& photon : photons_coll) { - if (!cut.template IsSelected(photon)) { - continue; - } - - float phi_cp = atan2(photon.vy(), photon.vx()); - float eta_cp = std::atanh(photon.vz() / sqrt(pow(photon.vx(), 2) + pow(photon.vy(), 2) + pow(photon.vz(), 2))); - value[0] = photon.pt(); - value[1] = photon.v0radius(); - value[2] = phi_cp > 0 ? phi_cp : phi_cp + TMath::TwoPi(); - value[3] = eta_cp; - reinterpret_cast(list_v0->FindObject(cut.GetName())->FindObject("hs_conv_point"))->Fill(value); - - } // end of photon loop - } // end of cut loop - - } // end of collision loop + fillEventInfo<1>(collision); + registry.fill(HIST("Event/before/hCollisionCounter"), 10.0); // accepted + registry.fill(HIST("Event/after/hCollisionCounter"), 10.0); // accepted + + auto v0sThisCollision = v0photons.sliceBy(perCollision, collision.globalIndex()); + + for (const auto& v0 : v0sThisCollision) { + + if (!fV0PhotonCut.IsSelected(v0)) { + continue; + } + + auto negLeg = v0.template negTrack_as(); + auto posLeg = v0.template posTrack_as(); + + int nClsNeg = 0; + int nClsPos = 0; + + if constexpr (requires { negLeg.itsNCls(); }) { + nClsNeg = negLeg.itsNCls(); + nClsPos = posLeg.itsNCls(); + } else if constexpr (requires { negLeg.itsClusterMap(); }) { + auto countBits = [](uint8_t map) { + return std::bitset<8>(map).count(); + }; + nClsNeg = countBits(negLeg.itsClusterMap()); + nClsPos = countBits(posLeg.itsClusterMap()); + } + if (nClsNeg >= 0 && nClsNeg <= 7) { // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillClusterNeg(registry, nClsNeg, v0.v0radius()); + } + if (nClsPos >= 0 && nClsPos <= 7) { // o2-linter: disable=magic-number (just numbers for ITS cluster) + fillClusterPos(registry, nClsPos, v0.v0radius()); + } + + if (cfgPlotMBGeneral) { + + registry.fill( + HIST("MBGeneral"), + v0.vz(), // 0 + v0.v0radius(), // 1 + v0.eta(), // 2 + v0.phi(), // 3 + v0.pt()); // 4 + } + + if (cfgPlotMBDetailed) { + + registry.fill( + HIST("MBStudiesWireLeft"), + v0.vx(), // 0 + v0.vy(), // 1 + v0.vz(), // 2 + v0.phi(), // 3 + v0.v0radius(), // 4 + v0.pt()); // 5 + + registry.fill( + HIST("MBStudiesWireRight"), + v0.vx(), // 0 + v0.vy(), // 1 + v0.vz(), // 2 + v0.phi(), // 3 + v0.v0radius(), // 4 + v0.pt()); // 5 + + registry.fill( + HIST("MBStudiesWireITS"), + v0.vx(), // 0 + v0.vy(), // 1 + v0.vz(), // 2 + v0.phi(), // 3 + v0.v0radius(), // 4 + v0.pt()); // 5 + + registry.fill( + HIST("MBStudiesMFT"), + v0.vx(), // 0 + v0.vy(), // 1 + v0.vz(), // 2 + v0.phi(), // 3 + v0.v0radius(), // 4 + v0.pt()); // 5 + + registry.fill( + HIST("MBStudiesTPC"), + v0.vx(), // 0 + v0.vy(), // 1 + v0.vz(), // 2 + v0.phi(), // 3 + v0.v0radius(), // 4 + v0.pt()); // 5 + } + } + + reconstructMesons(v0sThisCollision); + + auto electronsPerCollision = electrons.sliceBy(perCollisionElectron, collision.globalIndex()); + auto positronsPerCollision = positrons.sliceBy(perCollisionElectron, collision.globalIndex()); + + reconstructDalitz(v0sThisCollision, positronsPerCollision, electronsPerCollision, collision); + + auto key = getPoolBin(centralities[cfgCentEstimator], collision.posZ()); + auto& pool = mixingPools[key]; + reconstructMesonsMixed(v0sThisCollision, pool); + + std::vector eventCopy; + eventCopy.reserve(v0sThisCollision.size()); + for (const auto& v0 : v0sThisCollision) { + if (!fV0PhotonCut.IsSelected(v0)) { + continue; + } + eventCopy.push_back({v0.pt(), v0.eta(), v0.phi(), v0.vz(), v0.v0radius()}); + } + pool.emplace_back(std::move(eventCopy)); + + if (pool.size() > MaxMixEvents) { + pool.pop_front(); + } + } + } + + int classifyPurity(int pdgNeg, int pdgPos) + { + if (pdgNeg == kElectron && pdgPos == kPositron) { + return 1; + } + if (pdgNeg == kElectron && pdgPos == kPiPlus) + return 2; + if (pdgNeg == kPiMinus && pdgPos == kPositron) + return 3; + if (pdgNeg == kKMinus && pdgPos == kPositron) + return 4; + if (pdgNeg == kElectron && pdgPos == kKPlus) + return 5; + if (pdgNeg == kProtonBar && pdgPos == kProton) + return 6; + if (pdgNeg == kKMinus && pdgPos == kKPlus) + return 7; + if (pdgNeg == kPiMinus && pdgPos == kProton) + return 8; + if (pdgNeg == kProtonBar && pdgPos == kPiPlus) + return 9; + if (pdgNeg == kMuonMinus && pdgPos == kMuonPlus) + return 10; + return 11; } - template - void SameEventPairing(TEvents const& collisions, TPhotons1 const& photons1, TPhotons2 const& photons2, TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, TCuts1 const& tagcuts, TCuts2 const& probecuts, TPairCuts const& paircuts, TLegs const& /*legs*/) + Partition groupedCollisions = cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax; // this goes to same event. + Filter collisionFilterCommon = nabs(o2::aod::collision::posZ) < 10.f && o2::aod::collision::numContrib > static_cast(0) && o2::aod::evsel::sel8 == true; // o2-linter: disable=magic-number (vertex position) + Filter collisionFilterCentrality = (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); + using MyFilteredCollisions = soa::Filtered; // this goes to same event pairing. + + void processTrue(MyCollisions const&, MyFilteredCollisions const& filteredCollisions, + MyV0Photons const& v0photons, + MyMCV0Legs const&, + aod::EMMCParticles const& mcparticles, + aod::EMMCEvents const&) { - THashList* list_pair_ss = static_cast(fMainList->FindObject("Pair")->FindObject(pairnames[pairtype].data())); + for (const auto& collision : filteredCollisions) { - for (auto& collision : collisions) { - const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; - if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + const float centralities[3] = { + collision.centFT0M(), + collision.centFT0A(), + collision.centFT0C()}; + + if (centralities[cfgCentEstimator] < cfgCentMin || + cfgCentMax < centralities[cfgCentEstimator]) { continue; } + fillEventInfo<0>(collision); if (!fEMEventCut.IsSelected(collision)) { continue; } - auto photons1_coll = photons1.sliceBy(perCollision1, collision.globalIndex()); - auto photons2_coll = photons2.sliceBy(perCollision2, collision.globalIndex()); - - double value[6] = {0.f}; - float phi_cp2 = 0.f, eta_cp2 = 0.f; - for (auto& tagcut : tagcuts) { - for (auto& probecut : probecuts) { - for (auto& g1 : photons1_coll) { - for (auto& g2 : photons2_coll) { - if (g1.globalIndex() == g2.globalIndex()) { - continue; - } - if (!IsSelectedPair(g1, g2, tagcut, probecut)) { - continue; - } - for (auto& paircut : paircuts) { - if (!paircut.IsSelected(g1, g2)) { - continue; - } - - ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); // tag - ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); // probe - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - phi_cp2 = atan2(g2.vy(), g2.vx()); - eta_cp2 = std::atanh(g2.vz() / sqrt(pow(g2.vx(), 2) + pow(g2.vy(), 2) + pow(g2.vz(), 2))); - value[0] = v12.M(); - value[1] = g1.pt(); - value[2] = g2.pt(); - value[3] = g2.v0radius(); - value[4] = phi_cp2 > 0.f ? phi_cp2 : phi_cp2 + TMath::TwoPi(); - value[5] = eta_cp2; - reinterpret_cast(list_pair_ss->FindObject(Form("%s_%s", tagcut.GetName(), probecut.GetName()))->FindObject(paircut.GetName())->FindObject("hs_conv_point_same"))->Fill(value); - } // end of pair cut loop - } // end of g2 loop - } // end of g1 loop - } // end of probecut loop - } // end of tagcut loop - } // end of collision loop + fillEventInfo<1>(collision); + + registry.fill(HIST("Event/before/hCollisionCounter"), 10.0); + registry.fill(HIST("Event/after/hCollisionCounter"), 10.0); + + auto v0PhotonsColl = v0photons.sliceBy(perCollision, collision.globalIndex()); + + for (const auto& v0 : v0PhotonsColl) { + + auto pos = v0.posTrack_as(); + auto ele = v0.negTrack_as(); + + auto posmc = pos.template emmcparticle_as(); + auto elemc = ele.template emmcparticle_as(); + + if (!fV0PhotonCut.IsSelected(v0)) { + continue; + } + + int purityCat = 12; // default: unmatched + + purityCat = classifyPurity(elemc.pdgCode(), posmc.pdgCode()); + + float rConv = v0.v0radius(); + registry.fill(HIST("PhotonPurity"), + v0.pt(), rConv, v0.eta(), v0.phi(), purityCat); + + LOGP(info, "order: {}", posmc.pdgCode()); + + int photonid = FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 22, mcparticles); + if (photonid < 0) { + continue; + } + + auto mcphoton = mcparticles.iteratorAt(photonid); + + float rRec = v0.v0radius(); + float rGen = std::sqrt(mcphoton.vx() * mcphoton.vx() + + mcphoton.vy() * mcphoton.vy()); + float deltaR = rRec - rGen; + float deltaZ = v0.vz() - mcphoton.vz(); + float deltaPhi = v0.phi() - mcphoton.phi(); + float deltaEta = v0.eta() - mcphoton.eta(); + float deltaPt = v0.pt() - mcphoton.pt(); + + registry.fill(HIST("ResolutionGen/Z_res"), + mcphoton.vz(), rGen, mcphoton.phi(), + mcphoton.eta(), mcphoton.pt(), deltaZ); + + registry.fill(HIST("ResolutionGen/R_res"), + mcphoton.vz(), rGen, mcphoton.phi(), + mcphoton.eta(), mcphoton.pt(), deltaR); + + registry.fill(HIST("ResolutionGen/Phi_res"), + mcphoton.vz(), rGen, mcphoton.phi(), + mcphoton.eta(), mcphoton.pt(), deltaPhi); + + registry.fill(HIST("ResolutionGen/Eta_res"), + mcphoton.vz(), rGen, mcphoton.phi(), + mcphoton.eta(), mcphoton.pt(), deltaEta); + + registry.fill(HIST("ResolutionGen/Pt_res"), + mcphoton.vz(), rGen, mcphoton.phi(), + mcphoton.eta(), mcphoton.pt(), + deltaPt / mcphoton.pt()); + + registry.fill(HIST("ResolutionRec/Z_res"), + v0.vz(), rRec, v0.phi(), v0.eta(), + v0.pt(), deltaZ); + + registry.fill(HIST("ResolutionRec/R_res"), + v0.vz(), rRec, v0.phi(), v0.eta(), + v0.pt(), deltaR); + + registry.fill(HIST("ResolutionRec/Phi_res"), + v0.vz(), rRec, v0.phi(), v0.eta(), + v0.pt(), deltaPhi); + + registry.fill(HIST("ResolutionRec/Eta_res"), + v0.vz(), rRec, v0.phi(), v0.eta(), + v0.pt(), deltaEta); + + registry.fill(HIST("ResolutionRec/Pt_res"), + v0.vz(), rRec, v0.phi(), v0.eta(), + v0.pt(), deltaPt / mcphoton.pt()); + } + } } - Configurable ndepth{"ndepth", 10, "depth for event mixing"}; - ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; - ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 999.f}, "Mixing bins - centrality"}; - using BinningType_M = ColumnBinningPolicy; - using BinningType_A = ColumnBinningPolicy; - using BinningType_C = ColumnBinningPolicy; - BinningType_M colBinning_M{{ConfVtxBins, ConfCentBins}, true}; - BinningType_A colBinning_A{{ConfVtxBins, ConfCentBins}, true}; - BinningType_C colBinning_C{{ConfVtxBins, ConfCentBins}, true}; - - template - void MixedEventPairing(TEvents const& collisions, TPhotons1 const& photons1, TPhotons2 const& photons2, TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, TCuts1 const& tagcuts, TCuts2 const& probecuts, TPairCuts const& paircuts, TLegs const& /*legs*/, TMixedBinning const& colBinning) + void processCollAssoc(MyCollisionsMC const& collisions, + MyMCCollisions const&, + aod::EMMCParticles const& mcparticles, + MyV0Photons const& v0photons, + MyMCV0Legs const&) { - THashList* list_pair_ss = static_cast(fMainList->FindObject("Pair")->FindObject(pairnames[pairtype].data())); - for (auto& [collision1, collision2] : soa::selfCombinations(colBinning, ndepth, -1, collisions, collisions)) { // internally, CombinationsStrictlyUpperIndexPolicy(collisions, collisions) is called. + for (auto const& col : collisions) { + + auto mccollision = col.emmcevent_as(); + + auto v0s = v0photons.sliceBy(perCollision, col.globalIndex()); + + int mcColIdDerived = mccollision.globalIndex(); // MC event (derived) global index + + for (auto const& v0 : v0s) { - const float centralities1[3] = {collision1.centFT0M(), collision1.centFT0A(), collision1.centFT0C()}; - const float centralities2[3] = {collision2.centFT0M(), collision2.centFT0A(), collision2.centFT0C()}; + if (!fV0PhotonCut.IsSelected(v0)) { + continue; + } - if (centralities1[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities1[cfgCentEstimator]) { + auto negLeg = v0.negTrack_as(); + auto posLeg = v0.posTrack_as(); + + auto negMC = negLeg.template emmcparticle_as(); + auto posMC = posLeg.template emmcparticle_as(); + + // find the photon mother in MC + int photonId = FindCommonMotherFrom2Prongs(posMC, negMC, -11, 11, 22, mcparticles); + if (photonId < 0) { + continue; + } + + auto mcPhoton = mcparticles.iteratorAt(photonId); + auto trueMcCol = mcPhoton.emmcevent_as(); // MC event where the photon was generated + int trueMcColIndex = trueMcCol.globalIndex(); + + int deltaCol = mcColIdDerived - trueMcColIndex; + + float rRec = v0.v0radius(); + float rGen = std::sqrt(mcPhoton.vx() * mcPhoton.vx() + mcPhoton.vy() * mcPhoton.vy()); + float deltaR = rRec - rGen; + float deltaZ = v0.vz() - mcPhoton.vz(); + float deltaPt = v0.pt() - mcPhoton.pt(); + float relPtRes = (mcPhoton.pt() > 0.f) ? deltaPt / mcPhoton.pt() : 0.f; + + // Fill base debug + registry.fill(HIST("DeltaRecDerived"), deltaCol); + registry.fill(HIST("AllR"), rRec); + + if (deltaCol == 0) { + registry.fill(HIST("RightCollisions/SparseDeltaCol"), + rRec, + v0.pt(), + relPtRes, + deltaZ, + deltaR, + deltaCol); + } else { + registry.fill(HIST("WrongCollisions/SparseDeltaCol"), + rRec, + v0.pt(), + relPtRes, + deltaZ, + deltaR, + deltaCol); + } + + if (deltaCol == 0) { + registry.fill(HIST("RightCollisions/DeltaColvsZ"), v0.vz()); + registry.fill(HIST("RightCollisions/DeltaColvsR"), rRec); + registry.fill(HIST("RightCollisions/DeltaColvspT"), v0.pt()); + registry.fill(HIST("RightCollisions/DeltaColvsZRes"), deltaZ); + registry.fill(HIST("RightCollisions/DeltaColvsRRes"), deltaR); + registry.fill(HIST("RightCollisions/DeltaColvspTRes"), relPtRes); + } else { + registry.fill(HIST("WrongCollisions/DeltaColvsZ"), v0.vz()); + registry.fill(HIST("WrongCollisions/DeltaColvsR"), rRec); + registry.fill(HIST("WrongCollisions/DeltaColvspT"), v0.pt()); + registry.fill(HIST("WrongCollisions/DeltaColvsZRes"), deltaZ); + registry.fill(HIST("WrongCollisions/DeltaColvsRRes"), deltaR); + registry.fill(HIST("WrongCollisions/DeltaColvspTRes"), relPtRes); + } + } + } + } + + void processBremsstrahlung(MyTracks const& tracks, aod::McParticles const& mcparticles) + { + + for (const auto& trk : tracks) { + if (!trk.has_mcParticle()) { continue; } - if (centralities2[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities2[cfgCentEstimator]) { + auto mc = trk.mcParticle(); + if (std::abs(mc.pdgCode()) != kElectron) { continue; } - if (!fEMEventCut.IsSelected(collision1) || !fEMEventCut.IsSelected(collision2)) { + + if (trk.pt() < dileptoncuts.cfgMinPtTrack) { continue; } - auto photons_coll1 = photons1.sliceBy(perCollision1, collision1.globalIndex()); - auto photons_coll2 = photons2.sliceBy(perCollision2, collision2.globalIndex()); - // LOGF(info, "collision1: posZ = %f, numContrib = %d , sel8 = %d | collision2: posZ = %f, numContrib = %d , sel8 = %d", collision1.posZ(), collision1.numContrib(), collision1.sel8(), collision2.posZ(), collision2.numContrib(), collision2.sel8()); - - double value[6] = {0.f}; - float phi_cp2 = 0.f, eta_cp2 = 0.f; - for (auto& tagcut : tagcuts) { - for (auto& probecut : probecuts) { - for (auto& g1 : photons_coll1) { - for (auto& g2 : photons_coll2) { - if (!IsSelectedPair(g1, g2, tagcut, probecut)) { - continue; - } - // LOGF(info, "Mixed event photon pair: (%d, %d) from events (%d, %d), photon event: (%d, %d)", g1.index(), g2.index(), collision1.index(), collision2.index(), g1.globalIndex(), g2.globalIndex()); - - for (auto& paircut : paircuts) { - if (!paircut.IsSelected(g1, g2)) { - continue; - } - - ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); // tag - ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); // probe - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - phi_cp2 = atan2(g2.vy(), g2.vx()); - eta_cp2 = std::atanh(g2.vz() / sqrt(pow(g2.vx(), 2) + pow(g2.vy(), 2) + pow(g2.vz(), 2))); - value[0] = v12.M(); - value[1] = g1.pt(); - value[2] = g2.pt(); - value[3] = g2.v0radius(); - value[4] = phi_cp2 > 0.f ? phi_cp2 : phi_cp2 + TMath::TwoPi(); - value[5] = eta_cp2; - reinterpret_cast(list_pair_ss->FindObject(Form("%s_%s", tagcut.GetName(), probecut.GetName()))->FindObject(paircut.GetName())->FindObject("hs_conv_point_mix"))->Fill(value); - - } // end of pair cut loop - } // end of g2 loop - } // end of g1 loop - } // end of probecut loop - } // end of tagcut loop - } // end of different collision combinations - } + double pMC = mc.p(); + double ptMC = mc.pt(); + double phiMC = mc.phi(); + double etaMC = mc.eta(); - Partition grouped_collisions = cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax; // this goes to same event. - Filter collisionFilter_common = nabs(o2::aod::collision::posZ) < 10.f && o2::aod::collision::numContrib > (uint16_t)0 && o2::aod::evsel::sel8 == true; - 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); - using MyFilteredCollisions = soa::Filtered; // this goes to mixed event. + double pReco = trk.p(); + double ptReco = trk.pt(); + double phiReco = trk.phi(); + double etaReco = trk.eta(); - void processMB(MyCollisions const&, MyFilteredCollisions const& filtered_collisions, MyV0Photons const& v0photons, aod::V0Legs const& legs) - { - fillsinglephoton(grouped_collisions, v0photons, perCollision_pcm, fProbeCuts, legs); - SameEventPairing(filtered_collisions, v0photons, v0photons, perCollision_pcm, perCollision_pcm, fTagCuts, fProbeCuts, fPairCuts, legs); - if (fDoMixing) { - if (cfgCentEstimator == 0) { - MixedEventPairing(filtered_collisions, v0photons, v0photons, perCollision_pcm, perCollision_pcm, fTagCuts, fProbeCuts, fPairCuts, legs, colBinning_M); - } else if (cfgCentEstimator == 1) { - MixedEventPairing(filtered_collisions, v0photons, v0photons, perCollision_pcm, perCollision_pcm, fTagCuts, fProbeCuts, fPairCuts, legs, colBinning_A); - } else if (cfgCentEstimator == 2) { - MixedEventPairing(filtered_collisions, v0photons, v0photons, perCollision_pcm, perCollision_pcm, fTagCuts, fProbeCuts, fPairCuts, legs, colBinning_C); + int nBrem = 0; // o2-linter: disable=magic-number + + registry.fill(HIST("Bremsstrahlung/relativeResoPtWOBrems"), trk.pt(), trk.pt() * std::sqrt(trk.c1Pt21Pt2())); + + bool isFirst = true; + for (const auto& dId : mc.daughtersIds()) { + if (dId < 0 || dId >= mcparticles.size()) { + continue; + } + auto daughter = mcparticles.iteratorAt(dId); + if (daughter.getProcess() != kPBrem) { + continue; + } + + if (std::abs(daughter.eta()) > pcmcuts.cfgMaxEtaV0) { + continue; + } + + double r = std::hypot(daughter.vx(), daughter.vy()); + double z = daughter.vz(); + + nBrem++; + + registry.fill(HIST("Bremsstrahlung/EnergyLossVsR"), r, daughter.e()); + registry.fill(HIST("Bremsstrahlung/EnergyLossVsZ"), z, daughter.e()); + registry.fill(HIST("Bremsstrahlung/EnergyLossXY"), daughter.vx(), daughter.vy()); + registry.fill(HIST("Bremsstrahlung/EnergyLossXYWeigh"), daughter.vx(), daughter.vy(), daughter.e()); + + if (isFirst) { + registry.fill(HIST("Bremsstrahlung/relativeResoPtWBrems"), trk.pt(), trk.pt() * std::sqrt(trk.c1Pt21Pt2())); + } + + registry.fill(HIST("Bremsstrahlung/Sigma1PtVsR"), r, trk.sigma1Pt()); + + registry.fill(HIST("Bremsstrahlung/SigmaYVsR"), r, trk.sigmaY()); + registry.fill(HIST("Bremsstrahlung/SigmaZVsR"), r, trk.sigmaZ()); + + if (pMC > 0) { + registry.fill(HIST("Bremsstrahlung/DeltaPoverPvsR"), r, (pReco - pMC) / pMC); + } + if (ptMC > 0) { + registry.fill(HIST("Bremsstrahlung/DeltaPtvsR"), r, (ptReco - ptMC) / ptMC); + } + registry.fill(HIST("Bremsstrahlung/DeltaPhivsR"), r, phiReco - phiMC); + registry.fill(HIST("Bremsstrahlung/DeltaEtavsR"), r, etaReco - etaMC); } + + registry.fill(HIST("Bremsstrahlung/NBrem"), nBrem); + isFirst = false; } } - void processDummy(MyCollisions::iterator const&) {} - - PROCESS_SWITCH(MaterialBudget, processMB, "process material budget", false); - PROCESS_SWITCH(MaterialBudget, processDummy, "Dummy function", true); + PROCESS_SWITCH(MaterialBudget, processRec, "process material budget", true); + PROCESS_SWITCH(MaterialBudget, processTrue, "process material budget", false); + PROCESS_SWITCH(MaterialBudget, processCollAssoc, "process material budget", false); + PROCESS_SWITCH(MaterialBudget, processBremsstrahlung, "process material budget", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"material-budget"})}; + adaptAnalysisTask(cfgc)}; }