diff --git a/PWGEM/Dilepton/Core/DileptonSV.h b/PWGEM/Dilepton/Core/DileptonSV.h new file mode 100644 index 00000000000..8296506a5a9 --- /dev/null +++ b/PWGEM/Dilepton/Core/DileptonSV.h @@ -0,0 +1,2045 @@ +// 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. +// +// ======================== +// +// This code runs loop over leptons. +// Please write to: daiki.sekihata@cern.ch + +#ifndef PWGEM_DILEPTON_CORE_DILEPTONSV_H_ +#define PWGEM_DILEPTON_CORE_DILEPTONSV_H_ + +#include "PWGEM/Dilepton/Core/DielectronCut.h" +#include "PWGEM/Dilepton/Core/DimuonCut.h" +#include "PWGEM/Dilepton/Core/EMEventCut.h" +#include "PWGEM/Dilepton/DataModel/dileptonTables.h" +#include "PWGEM/Dilepton/Utils/EMFwdTrack.h" +#include "PWGEM/Dilepton/Utils/EMTrack.h" +#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" +#include "PWGEM/Dilepton/Utils/EventHistograms.h" +#include "PWGEM/Dilepton/Utils/EventMixingHandler.h" +#include "PWGEM/Dilepton/Utils/PairUtilities.h" + +#include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/Zorro.h" +#include "Common/Core/fwdtrackUtilities.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.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 // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using MyCollisions = o2::soa::Join; +using MyCollision = MyCollisions::iterator; + +using MyElectrons = o2::soa::Join; +using MyElectron = MyElectrons::iterator; +using FilteredMyElectrons = o2::soa::Filtered; +using FilteredMyElectron = FilteredMyElectrons::iterator; + +using MyElectronsSCT = o2::soa::Join; +using MyElectronSCT = MyElectronsSCT::iterator; +using FilteredMyElectronsSCT = o2::soa::Filtered; +using FilteredMyElectronSCT = FilteredMyElectronsSCT::iterator; + +using MyMuons = o2::soa::Join; +using MyMuon = MyMuons::iterator; +using FilteredMyMuons = o2::soa::Filtered; +using FilteredMyMuon = FilteredMyMuons::iterator; + +using MyEMH_electron = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, o2::aod::pwgem::dilepton::utils::EMTrack>; +using MyEMH_muon = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, o2::aod::pwgem::dilepton::utils::EMFwdTrack>; + +template +struct DileptonSV { + + // Configurables + o2::framework::Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + o2::framework::Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + o2::framework::Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + o2::framework::Configurable cfgApplySPresolution{"cfgApplySPresolution", false, "flag to apply resolution correction for flow analysis"}; + o2::framework::Configurable spresoPath{"spresoPath", "Users/d/dsekihat/PWGEM/dilepton/Qvector/resolution/LHC23zzh/pass3/test", "Path to SP resolution file"}; + o2::framework::Configurable spresoHistName{"spresoHistName", "h1_R2_FT0M_BPos_BNeg", "histogram name of SP resolution file"}; + + o2::framework::Configurable cfgAnalysisType{"cfgAnalysisType", static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kQC), "kQC:0, kUPC:1, kFlowV2SP:2, kFlowV3SP:3, kPolarization:4, kHFll:5, kBootstrapv2:6, kFlowV2EP:7"}; + o2::framework::Configurable cfgQvecEstimator{"cfgQvecEstimator", 2, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"}; + o2::framework::Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; + + // for mixing + o2::framework::Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; + o2::framework::Configurable cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"}; + o2::framework::Configurable cfgDoMix{"cfgDoMix", true, "flag for event mixing"}; + o2::framework::Configurable ndepth{"ndepth", 1000, "depth for event mixing"}; + o2::framework::Configurable ndiff_bc_mix{"ndiff_bc_mix", 594, "difference in global BC required in mixed events"}; + o2::framework::ConfigurableAxis ConfVtxBins{"ConfVtxBins", {o2::framework::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"}; + o2::framework::ConfigurableAxis ConfCentBins{"ConfCentBins", {o2::framework::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.f, 999.f}, "Mixing bins - centrality"}; + o2::framework::ConfigurableAxis ConfEPBins{"ConfEPBins", {16, -M_PI / 2, +M_PI / 2}, "Mixing bins - event plane angle"}; + o2::framework::ConfigurableAxis ConfOccupancyBins{"ConfOccupancyBins", {o2::framework::VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; + + o2::framework::Configurable cfgApplyWeightTTCA{"cfgApplyWeightTTCA", false, "flag to apply weighting by 1/N"}; + o2::framework::Configurable cfgPolarizationFrame{"cfgPolarizationFrame", 0, "frame of polarization. 0:CS, 1:HX, else:FATAL"}; + + o2::framework::ConfigurableAxis ConfMllBins{"ConfMllBins", {o2::framework::VARIABLE_WIDTH, 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.40, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.50, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.60, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.10, 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.75, 2.80, 2.85, 2.90, 2.95, 3.00, 3.05, 3.10, 3.15, 3.20, 3.25, 3.30, 3.35, 3.40, 3.45, 3.50, 3.55, 3.60, 3.65, 3.70, 3.75, 3.80, 3.85, 3.90, 3.95, 4.00}, "mll bins for output histograms"}; + o2::framework::ConfigurableAxis ConfPtllBins{"ConfPtllBins", {o2::framework::VARIABLE_WIDTH, 0.00, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.50, 3.00, 3.50, 4.00, 4.50, 5.00, 6.00, 7.00, 8.00, 9.00, 10.00}, "pTll bins for output histograms"}; + o2::framework::ConfigurableAxis ConfDCAllBins{"ConfDCAllBins", {o2::framework::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, 1.6, 1.7, 1.8, 1.9, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}, "DCAll bins for output histograms"}; + o2::framework::ConfigurableAxis ConfYllBins{"ConfYllBins", {1, -1.f, 1.f}, "yll bins for output histograms"}; // pair rapidity + o2::framework::ConfigurableAxis ConfLog10Chi2PCABins{"ConfLog10Chi2PCABins", {100, -10.f, 0.f}, "log10 of chi2PCA bins for output histograms"}; + o2::framework::ConfigurableAxis ConfDLBins{"ConfDLBins", {100, 0.f, 10.f}, "decay length bins for output histograms"}; + + o2::framework::ConfigurableAxis ConfSPBins{"ConfSPBins", {200, -5, 5}, "SP bins for flow analysis"}; + o2::framework::ConfigurableAxis ConfPolarizationCosThetaBins{"ConfPolarizationCosThetaBins", {20, -1.f, 1.f}, "cos(theta) bins for polarization analysis"}; + o2::framework::ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, -M_PI, M_PI}, "phi bins for polarization analysis"}; + o2::framework::ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {15, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2> + + o2::framework::Configurable cfgNumBootstrapSamples{"cfgNumBootstrapSamples", 1, "Number of Bootstrap Samples"}; + + EMEventCut fEMEventCut; + struct : o2::framework::ConfigurableGroup { + std::string prefix = "eventcut_group"; + o2::framework::Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; + o2::framework::Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; + o2::framework::Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + o2::framework::Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + o2::framework::Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + o2::framework::Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + o2::framework::Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + o2::framework::Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + o2::framework::Configurable cfgRequireVertexTOFmatched{"cfgRequireVertexTOFmatched", false, "require Vertex TOFmatched in event cut"}; // ITS-TPC-TOF matched track contributes PV. + o2::framework::Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + o2::framework::Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. occupancy"}; + o2::framework::Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. occupancy"}; + o2::framework::Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2, "min. FT0C occupancy"}; + o2::framework::Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + o2::framework::Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; + o2::framework::Configurable cfgRequireNoCollInTimeRangeStrict{"cfgRequireNoCollInTimeRangeStrict", false, "require no collision in time range strict"}; + o2::framework::Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "require no collision in time range standard"}; + o2::framework::Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "require no collision in time range strict"}; + o2::framework::Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "require no HM collision in previous ITS ROF"}; + o2::framework::Configurable cfgRequireGoodITSLayer3{"cfgRequireGoodITSLayer3", false, "number of inactive chips on ITS layer 3 are below threshold "}; + o2::framework::Configurable cfgRequireGoodITSLayer0123{"cfgRequireGoodITSLayer0123", false, "number of inactive chips on ITS layers 0-3 are below threshold "}; + o2::framework::Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "number of inactive chips on all ITS layers are below threshold "}; + // for RCT + o2::framework::Configurable cfgRequireGoodRCT{"cfgRequireGoodRCT", false, "require good detector flag in run condtion table"}; + o2::framework::Configurable cfgRCTLabel{"cfgRCTLabel", "CBT_hadronPID", "select 1 [CBT, CBT_hadronPID, CBT_muon_glo] see O2Physics/Common/CCDB/RCTSelectionFlags.h"}; + o2::framework::Configurable cfgCheckZDC{"cfgCheckZDC", false, "set ZDC flag for PbPb"}; + o2::framework::Configurable cfgTreatLimitedAcceptanceAsBad{"cfgTreatLimitedAcceptanceAsBad", false, "reject all events where the detectors relevant for the specified Runlist are flagged as LimitedAcceptance"}; + + o2::framework::Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; + o2::framework::Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + o2::framework::Configurable cfgNumContribMin{"cfgNumContribMin", 0, "min. numContrib"}; + o2::framework::Configurable cfgNumContribMax{"cfgNumContribMax", 65000, "max. numContrib"}; + } eventcuts; + + DielectronCut fDielectronCut; + struct : o2::framework::ConfigurableGroup { + std::string prefix = "dielectroncut_group"; + o2::framework::Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + o2::framework::Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; + o2::framework::Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pT"}; + o2::framework::Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pT"}; + o2::framework::Configurable cfg_min_pair_y{"cfg_min_pair_y", -0.8, "min pair rapidity"}; + o2::framework::Configurable cfg_max_pair_y{"cfg_max_pair_y", +0.8, "max pair rapidity"}; + o2::framework::Configurable cfg_min_pair_dca3d{"cfg_min_pair_dca3d", 0.0, "min pair dca3d in sigma"}; + o2::framework::Configurable cfg_max_pair_dca3d{"cfg_max_pair_dca3d", 1e+10, "max pair dca3d in sigma"}; + o2::framework::Configurable cfg_apply_phiv{"cfg_apply_phiv", true, "flag to apply phiv cut"}; + o2::framework::Configurable cfg_phiv_slope{"cfg_phiv_slope", 0.0185, "slope for m vs. phiv"}; + o2::framework::Configurable cfg_phiv_intercept{"cfg_phiv_intercept", -0.0280, "intercept for m vs. phiv"}; + o2::framework::Configurable cfg_min_phiv{"cfg_min_phiv", 0.0, "min phiv (constant)"}; + o2::framework::Configurable cfg_max_phiv{"cfg_max_phiv", 3.2, "max phiv (constant)"}; + o2::framework::Configurable cfg_apply_detadphi{"cfg_apply_detadphi", false, "flag to apply deta-dphi elliptic cut at PV"}; + o2::framework::Configurable cfg_apply_detadphiPosition{"cfg_apply_detadphiPosition", false, "flag to apply deta-dphi elliptic cut at ref. radius"}; + o2::framework::Configurable cfg_min_deta{"cfg_min_deta", 0.02, "min deta between 2 electrons (elliptic cut)"}; + o2::framework::Configurable cfg_min_dphi{"cfg_min_dphi", 0.2, "min dphi between 2 electrons (elliptic cut)"}; + o2::framework::Configurable cfg_min_opang{"cfg_min_opang", 0.0, "min opening angle"}; + o2::framework::Configurable cfg_max_opang{"cfg_max_opang", 6.4, "max opening angle"}; + // o2::framework::Configurable cfg_require_diff_sides{"cfg_require_diff_sides", false, "flag to require 2 tracks are from different sides."}; + + o2::framework::Configurable cfg_apply_cuts_from_prefilter{"cfg_apply_cuts_from_prefilter", false, "flag to apply prefilter set when producing derived data"}; + o2::framework::Configurable cfg_prefilter_bits{"cfg_prefilter_bits", 0, "prefilter bits [kNone : 0, kElFromPC : 1, kElFromPi0_20MeV : 2, kElFromPi0_40MeV : 4, kElFromPi0_60MeV : 8, kElFromPi0_80MeV : 16, kElFromPi0_100MeV : 32, kElFromPi0_120MeV : 64, kElFromPi0_140MeV : 128] Please consider logical-OR among them."}; // see PairUtilities.h + + o2::framework::Configurable cfg_apply_cuts_from_prefilter_derived{"cfg_apply_cuts_from_prefilter_derived", false, "flag to apply prefilter set in derived data"}; + o2::framework::Configurable cfg_prefilter_bits_derived{"cfg_prefilter_bits_derived", 0, "prefilter bits [kNone : 0, kMee : 1, kPhiV : 2, kSplitOrMergedTrackLS : 4, kSplitOrMergedTrackULS : 8] Please consider logical-OR among them."}; // see PairUtilities.h + + o2::framework::Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"}; + o2::framework::Configurable cfg_max_pt_track{"cfg_max_pt_track", 1e+10, "max pT for single track"}; + o2::framework::Configurable cfg_min_eta_track{"cfg_min_eta_track", -0.8, "min eta for single track"}; + o2::framework::Configurable cfg_max_eta_track{"cfg_max_eta_track", +0.8, "max eta for single track"}; + o2::framework::Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.f, "min phi for single track"}; + o2::framework::Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for single track"}; + o2::framework::Configurable cfg_mirror_phi_track{"cfg_mirror_phi_track", false, "mirror the phi cut around Pi, min and max Phi should be in 0-Pi"}; + o2::framework::Configurable cfg_reject_phi_track{"cfg_reject_phi_track", false, "reject the phi interval"}; + o2::framework::Configurable cfg_min_ncluster_tpc{"cfg_min_ncluster_tpc", 0, "min ncluster tpc"}; + o2::framework::Configurable cfg_min_ncluster_its{"cfg_min_ncluster_its", 5, "min ncluster its"}; + o2::framework::Configurable cfg_min_ncrossedrows{"cfg_min_ncrossedrows", 100, "min ncrossed rows"}; + o2::framework::Configurable cfg_max_frac_shared_clusters_tpc{"cfg_max_frac_shared_clusters_tpc", 999.f, "max fraction of shared clusters in TPC"}; + o2::framework::Configurable cfg_max_chi2tpc{"cfg_max_chi2tpc", 4.0, "max chi2/NclsTPC"}; + o2::framework::Configurable cfg_max_chi2its{"cfg_max_chi2its", 5.0, "max chi2/NclsITS"}; + o2::framework::Configurable cfg_max_chi2tof{"cfg_max_chi2tof", 1e+10, "max chi2 TOF"}; + o2::framework::Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 1.0, "max dca XY for single track in cm"}; + o2::framework::Configurable cfg_max_dcaz{"cfg_max_dcaz", 1.0, "max dca Z for single track in cm"}; + o2::framework::Configurable cfg_max_dca_sigma{"cfg_max_dca_sigma", 1e+10, "max dca for single track in sigma"}; + o2::framework::Configurable cfg_require_itsib_any{"cfg_require_itsib_any", false, "flag to require ITS ib any hits"}; + o2::framework::Configurable cfg_require_itsib_1st{"cfg_require_itsib_1st", true, "flag to require ITS ib 1st hit"}; + o2::framework::Configurable cfg_min_its_cluster_size{"cfg_min_its_cluster_size", 0.f, "min ITS cluster size"}; + o2::framework::Configurable cfg_max_its_cluster_size{"cfg_max_its_cluster_size", 16.f, "max ITS cluster size"}; + // o2::framework::Configurable cfg_min_rel_diff_pin{"cfg_min_rel_diff_pin", -1e+10, "min rel. diff. between pin and ppv"}; + // o2::framework::Configurable cfg_max_rel_diff_pin{"cfg_max_rel_diff_pin", +1e+10, "max rel. diff. between pin and ppv"}; + o2::framework::Configurable cfgRefR{"cfgRefR", 0.50, "ref. radius (m) for calculating phi position"}; // 0.50 +/- 0.06 can be syst. unc. + o2::framework::Configurable cfg_min_phiposition_track{"cfg_min_phiposition_track", 0.f, "min phi position for single track at certain radius"}; + o2::framework::Configurable cfg_max_phiposition_track{"cfg_max_phiposition_track", 6.3, "max phi position for single track at certain radius"}; + o2::framework::Configurable cfgDCAType{"cfgDCAType", 0, "type of DCA for output. 0:3D, 1:XY, 2:Z, else:3D"}; + + o2::framework::Configurable cfg_pid_scheme{"cfg_pid_scheme", static_cast(DielectronCut::PIDSchemes::kTPChadrejORTOFreq), "pid scheme [kTOFreq : 0, kTPChadrej : 1, kTPChadrejORTOFreq : 2, kTPConly : 3, kTOFif : 4, kPIDML : 5, kTPChadrejORTOFreq_woTOFif : 6]"}; + o2::framework::Configurable cfg_min_TPCNsigmaEl{"cfg_min_TPCNsigmaEl", -2.0, "min. TPC n sigma for electron inclusion"}; + o2::framework::Configurable cfg_max_TPCNsigmaEl{"cfg_max_TPCNsigmaEl", +3.0, "max. TPC n sigma for electron inclusion"}; + // o2::framework::Configurable cfg_min_TPCNsigmaMu{"cfg_min_TPCNsigmaMu", -0.0, "min. TPC n sigma for muon exclusion"}; + // o2::framework::Configurable cfg_max_TPCNsigmaMu{"cfg_max_TPCNsigmaMu", +0.0, "max. TPC n sigma for muon exclusion"}; + o2::framework::Configurable cfg_min_TPCNsigmaPi{"cfg_min_TPCNsigmaPi", -1e+10, "min. TPC n sigma for pion exclusion"}; + o2::framework::Configurable cfg_max_TPCNsigmaPi{"cfg_max_TPCNsigmaPi", +3.0, "max. TPC n sigma for pion exclusion"}; + o2::framework::Configurable cfg_min_TPCNsigmaKa{"cfg_min_TPCNsigmaKa", -3.0, "min. TPC n sigma for kaon exclusion"}; + o2::framework::Configurable cfg_max_TPCNsigmaKa{"cfg_max_TPCNsigmaKa", +3.0, "max. TPC n sigma for kaon exclusion"}; + o2::framework::Configurable cfg_min_TPCNsigmaPr{"cfg_min_TPCNsigmaPr", -3.0, "min. TPC n sigma for proton exclusion"}; + o2::framework::Configurable cfg_max_TPCNsigmaPr{"cfg_max_TPCNsigmaPr", +3.0, "max. TPC n sigma for proton exclusion"}; + o2::framework::Configurable cfg_min_TOFNsigmaEl{"cfg_min_TOFNsigmaEl", -3.0, "min. TOF n sigma for electron inclusion"}; + o2::framework::Configurable cfg_max_TOFNsigmaEl{"cfg_max_TOFNsigmaEl", +3.0, "max. TOF n sigma for electron inclusion"}; + o2::framework::Configurable cfg_min_pin_pirejTPC{"cfg_min_pin_pirejTPC", 0.f, "min. pin for pion rejection in TPC"}; + o2::framework::Configurable cfg_max_pin_pirejTPC{"cfg_max_pin_pirejTPC", 1e+10, "max. pin for pion rejection in TPC"}; + o2::framework::Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; + + // configuration for PID ML + // o2::framework::Configurable> onnxFileNames{"onnxFileNames", std::vector{"filename"}, "ONNX file names for each bin (if not from CCDB full path)"}; + // o2::framework::Configurable> onnxPathsCCDB{"onnxPathsCCDB", std::vector{"path"}, "Paths of models on CCDB"}; + o2::framework::Configurable> binsMLPID{"binsMLPID", std::vector{0.1, 0.15, 0.2, 0.25, 0.4, 0.8, 1.6, 2.0, 20.f}, "Bin limits for ML application"}; + o2::framework::Configurable> cutsMLPID{"cutsMLPID", std::vector{0.97, 0.97, 0.97, 0.8, 0.95, 0.95, 0.8, 0.8}, "ML cuts per bin"}; + // o2::framework::Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature"}, "Names of ML model input features"}; + // o2::framework::Configurable nameBinningFeature{"nameBinningFeature", "pt", "Names of ML model binning feature"}; + // o2::framework::Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB. Exceptions: > 0 for the specific timestamp, 0 gets the run dependent timestamp"}; + // o2::framework::Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + // o2::framework::Configurable enableOptimizations{"enableOptimizations", false, "Enables the ONNX extended model-optimization: sessionOptions.SetGraphOptimizationLevel(GraphOptimizationLevel::ORT_ENABLE_EXTENDED)"}; + + o2::framework::Configurable> binsMLSCT{"binsMLSCT", std::vector{0.1, 0.4, 0.6, 0.8, 1, 2, 4, 20.f}, "Bin limits for ML application"}; + o2::framework::Configurable> cutsMLSCTeT_prompt_hc{"cutsMLSCTeT_prompt_hc", std::vector{0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8}, "ML cuts per bin for prompt hc"}; + o2::framework::Configurable> cutsMLSCTeT_nonprompt_hc{"cutsMLSCTeT_nonprompt_hc", std::vector{0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8}, "ML cuts per bin for nonprompt hc"}; + o2::framework::Configurable> cutsMLSCTeT_hb{"cutsMLSCTeT_hb", std::vector{0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8}, "ML cuts per bin for hb"}; + } dielectroncuts; + + DimuonCut fDimuonCut; + struct : o2::framework::ConfigurableGroup { + std::string prefix = "dimuoncut_group"; + o2::framework::Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + o2::framework::Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; + o2::framework::Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pt"}; + o2::framework::Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pt"}; + o2::framework::Configurable cfg_min_pair_y{"cfg_min_pair_y", -4.0, "min pair rapidity"}; + o2::framework::Configurable cfg_max_pair_y{"cfg_max_pair_y", -2.5, "max pair rapidity"}; + o2::framework::Configurable cfg_min_pair_dcaxy{"cfg_min_pair_dcaxy", 0.0, "min pair dca3d in sigma"}; + o2::framework::Configurable cfg_max_pair_dcaxy{"cfg_max_pair_dcaxy", 1e+10, "max pair dca3d in sigma"}; + o2::framework::Configurable cfg_apply_detadphi{"cfg_apply_detadphi", false, "flag to apply deta-dphi elliptic cut"}; + o2::framework::Configurable cfg_min_deta{"cfg_min_deta", 0.02, "min deta between 2 muons (elliptic cut)"}; + o2::framework::Configurable cfg_min_dphi{"cfg_min_dphi", 0.02, "min dphi between 2 muons (elliptic cut)"}; + + o2::framework::Configurable cfg_apply_cuts_from_prefilter_derived{"cfg_apply_cuts_from_prefilter_derived", false, "flag to apply prefilter set in derived data"}; + o2::framework::Configurable cfg_prefilter_bits_derived{"cfg_prefilter_bits_derived", 0, "prefilter bits [kNone : 0, kSplitOrMergedTrackLS : 4, kSplitOrMergedTrackULS : 8] Please consider logical-OR among them."}; // see PairUtilities.h + + o2::framework::Configurable cfg_track_type{"cfg_track_type", 3, "muon track type [0: MFT-MCH-MID, 3: MCH-MID]"}; + o2::framework::Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"}; + o2::framework::Configurable cfg_max_pt_track{"cfg_max_pt_track", 1e+10, "max pT for single track"}; + o2::framework::Configurable cfg_min_eta_track{"cfg_min_eta_track", -4.0, "min eta for single track"}; + o2::framework::Configurable cfg_max_eta_track{"cfg_max_eta_track", -2.5, "max eta for single track"}; + o2::framework::Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.f, "min phi for single track"}; + o2::framework::Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for single track"}; + o2::framework::Configurable cfg_min_ncluster_mft{"cfg_min_ncluster_mft", 5, "min ncluster MFT"}; + o2::framework::Configurable cfg_min_ncluster_mch{"cfg_min_ncluster_mch", 5, "min ncluster MCH"}; + o2::framework::Configurable cfg_max_chi2{"cfg_max_chi2", 1e+6, "max chi2/ndf"}; + o2::framework::Configurable cfg_max_chi2mft{"cfg_max_chi2mft", 1e+6, "max chi2/ndf"}; + // o2::framework::Configurable cfg_max_matching_chi2_mftmch{"cfg_max_matching_chi2_mftmch", 40, "max chi2 for MFT-MCH matching"}; + o2::framework::Configurable cfg_border_pt_for_chi2mchmft{"cfg_border_pt_for_chi2mchmft", 0, "border pt for different max chi2 for MFT-MCH matching"}; + o2::framework::Configurable cfg_max_matching_chi2_mftmch_lowPt{"cfg_max_matching_chi2_mftmch_lowPt", 8, "max chi2 for MFT-MCH matching for low pT"}; + o2::framework::Configurable cfg_max_matching_chi2_mftmch_highPt{"cfg_max_matching_chi2_mftmch_highPt", 40, "max chi2 for MFT-MCH matching for high pT"}; + o2::framework::Configurable cfg_max_matching_chi2_mchmid{"cfg_max_matching_chi2_mchmid", 1e+10, "max chi2 for MCH-MID matching"}; + o2::framework::Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 1e+10, "max dca XY for single track in cm"}; + o2::framework::Configurable cfg_min_rabs{"cfg_min_rabs", 17.6, "min Radius at the absorber end"}; + o2::framework::Configurable cfg_max_rabs{"cfg_max_rabs", 89.5, "max Radius at the absorber end"}; + o2::framework::Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; + o2::framework::Configurable cfg_max_relDPt_wrt_matchedMCHMID{"cfg_max_relDPt_wrt_matchedMCHMID", 1e+10f, "max. relative dpt between MFT-MCH-MID and MCH-MID"}; + o2::framework::Configurable cfg_max_DEta_wrt_matchedMCHMID{"cfg_max_DEta_wrt_matchedMCHMID", 1e+10f, "max. deta between MFT-MCH-MID and MCH-MID"}; + o2::framework::Configurable cfg_max_DPhi_wrt_matchedMCHMID{"cfg_max_DPhi_wrt_matchedMCHMID", 1e+10f, "max. dphi between MFT-MCH-MID and MCH-MID"}; + o2::framework::Configurable requireMFTHitMap{"requireMFTHitMap", false, "flag to apply MFT hit map"}; + o2::framework::Configurable> requiredMFTDisks{"requiredMFTDisks", std::vector{0}, "hit map on MFT disks [0,1,2,3,4]. logical-OR of each double-sided disk"}; + } dimuoncuts; + + struct : o2::framework::ConfigurableGroup { + std::string prefix = "zorroGroup"; + o2::framework::Configurable cfg_swt_name{"cfg_swt_name", "fLMeeIMR", "desired software trigger name. 1 trigger per 1 task."}; // 1 trigger per 1 task + o2::framework::Configurable ccdbPathSoftwareTrigger{"ccdbPathSoftwareTrigger", "EventFiltering/Zorro/", "ccdb path for ZORRO objects"}; + o2::framework::Configurable bcMarginForSoftwareTrigger{"bcMarginForSoftwareTrigger", 100, "Number of BCs of margin for software triggers"}; + } zorroGroup; + + Zorro zorro; + int mToIidx = 0; + int mTOICounter = 0; + int mATCounter = 0; + + struct : o2::framework::ConfigurableGroup { + std::string prefix = "dfGroup"; + o2::framework::Configurable useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"}; + o2::framework::Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; + o2::framework::Configurable maxR{"maxR", 10.f, "reject PCA's above this radius"}; + o2::framework::Configurable maxDZIni{"maxDZIni", 4.f, "reject (if > 0) PCA candidate if tracks DZ exceeds threshold"}; + o2::framework::Configurable minParamChange{"minParamChange", 1e-3, "stop iterations if largest change of any X is smaller than this"}; + o2::framework::Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + o2::framework::Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; + o2::framework::Configurable maxLog10Chi2PCA{"maxLog10Chi2PCA", 0, "max. log10(chi2PCA) for dilepton pair"}; + } dfGroup; // for DCAFitterN + + struct : o2::framework::ConfigurableGroup { + std::string prefix = "fdfGroup"; + o2::framework::Configurable useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"}; + o2::framework::Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; + o2::framework::Configurable maxR{"maxR", 10.f, "reject PCA's above this radius"}; + o2::framework::Configurable maxDXIni{"maxDXIni", 4.f, "reject (if > 0) PCA candidate if tracks DZ exceeds threshold"}; + o2::framework::Configurable minParamChange{"minParamChange", 1e-3, "stop iterations if largest change of any X is smaller than this"}; + o2::framework::Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + o2::framework::Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; + o2::framework::Configurable maxLog10Chi2PCA{"maxLog10Chi2PCA", 0, "max. log10(chi2PCA) for dilepton pair"}; + } fdfGroup; // for FwdDCAFitterN + + o2::vertexing::DCAFitterN<2> mDCAFitter; + o2::vertexing::FwdDCAFitterN<2> mFwdDCAFitter; + + o2::aod::rctsel::RCTFlagsChecker rctChecker; + // o2::ccdb::CcdbApi ccdbApi; + o2::framework::Service ccdb; + int mRunNumber{0}; + float d_bz{0}; + + o2::framework::HistogramRegistry fRegistry{"output", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, false}; + // static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; + static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"}; + + std::mt19937 engine; + std::vector cent_bin_edges; + std::vector zvtx_bin_edges; + std::vector ep_bin_edges; + std::vector occ_bin_edges; + + int nmod = -1; // this is for flow analysis + float leptonM1 = 0.f; + float leptonM2 = 0.f; + + float beamM1 = o2::constants::physics::MassProton; // mass of beam + float beamM2 = o2::constants::physics::MassProton; // mass of beam + float beamE1 = 0.f; // beam energy + float beamE2 = 0.f; // beam energy + float beamP1 = 0.f; // beam momentum + float beamP2 = 0.f; // beam momentum + TH2D* h2sp_resolution = nullptr; + + void init(o2::framework::InitContext& /*context*/) + { + mRunNumber = 0; + d_bz = 0; + mToIidx = 0; + mTOICounter = 0; + mATCounter = 0; + + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + rctChecker.init(eventcuts.cfgRCTLabel.value, eventcuts.cfgCheckZDC.value, eventcuts.cfgTreatLimitedAcceptanceAsBad.value); + + if (ConfVtxBins.value[0] == o2::framework::VARIABLE_WIDTH) { + zvtx_bin_edges = std::vector(ConfVtxBins.value.begin(), ConfVtxBins.value.end()); + zvtx_bin_edges.erase(zvtx_bin_edges.begin()); + for (const auto& edge : zvtx_bin_edges) { + LOGF(info, "o2::framework::VARIABLE_WIDTH: zvtx_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(ConfVtxBins.value[0]); + float xmin = static_cast(ConfVtxBins.value[1]); + float xmax = static_cast(ConfVtxBins.value[2]); + zvtx_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + zvtx_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: zvtx_bin_edges[%d] = %f", i, zvtx_bin_edges[i]); + } + } + + if (ConfCentBins.value[0] == o2::framework::VARIABLE_WIDTH) { + cent_bin_edges = std::vector(ConfCentBins.value.begin(), ConfCentBins.value.end()); + cent_bin_edges.erase(cent_bin_edges.begin()); + for (const auto& edge : cent_bin_edges) { + LOGF(info, "o2::framework::VARIABLE_WIDTH: cent_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(ConfCentBins.value[0]); + float xmin = static_cast(ConfCentBins.value[1]); + float xmax = static_cast(ConfCentBins.value[2]); + cent_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + cent_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: cent_bin_edges[%d] = %f", i, cent_bin_edges[i]); + } + } + + if (ConfEPBins.value[0] == o2::framework::VARIABLE_WIDTH) { + ep_bin_edges = std::vector(ConfEPBins.value.begin(), ConfEPBins.value.end()); + ep_bin_edges.erase(ep_bin_edges.begin()); + for (const auto& edge : ep_bin_edges) { + LOGF(info, "o2::framework::VARIABLE_WIDTH: ep_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(ConfEPBins.value[0]); + float xmin = static_cast(ConfEPBins.value[1]); + float xmax = static_cast(ConfEPBins.value[2]); + ep_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + ep_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: ep_bin_edges[%d] = %f", i, ep_bin_edges[i]); + } + } + + LOGF(info, "cfgOccupancyEstimator = %d", cfgOccupancyEstimator.value); + if (ConfOccupancyBins.value[0] == o2::framework::VARIABLE_WIDTH) { + occ_bin_edges = std::vector(ConfOccupancyBins.value.begin(), ConfOccupancyBins.value.end()); + occ_bin_edges.erase(occ_bin_edges.begin()); + for (const auto& edge : occ_bin_edges) { + LOGF(info, "o2::framework::VARIABLE_WIDTH: occ_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(ConfOccupancyBins.value[0]); + float xmin = static_cast(ConfOccupancyBins.value[1]); + float xmax = static_cast(ConfOccupancyBins.value[2]); + occ_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + occ_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: occ_bin_edges[%d] = %f", i, occ_bin_edges[i]); + } + } + + emh_pos = new TEMH(ndepth); + emh_neg = new TEMH(ndepth); + + DefineEMEventCut(); + addhistograms(); + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + DefineDielectronCut(); + leptonM1 = o2::constants::physics::MassElectron; + leptonM2 = o2::constants::physics::MassElectron; + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + DefineDimuonCut(); + leptonM1 = o2::constants::physics::MassMuon; + leptonM2 = o2::constants::physics::MassMuon; + } + + fRegistry.add("Pair/mix/hDiffBC", "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", o2::framework::HistType::kTH1D, {{10001, -0.5, 10000.5}}, true); + + if (doprocessTriggerAnalysis) { + LOGF(info, "Trigger analysis is enabled. Desired trigger name = %s", zorroGroup.cfg_swt_name.value.data()); + fRegistry.add("Event/trigger/hInspectedTVX", "inspected TVX;run number;N_{TVX}", o2::framework::HistType::kTProfile, {{100000, 500000.5, 600000.5}}, true); // extend X range in Run 4/5 + fRegistry.add("Event/trigger/hScaler", "trigger counter before DS;run number;counter", o2::framework::HistType::kTProfile, {{100000, 500000.5, 600000.5}}, true); // extend X range in Run 4/5 + fRegistry.add("Event/trigger/hSelection", "trigger counter after DS;run number;counter", o2::framework::HistType::kTProfile, {{100000, 500000.5, 600000.5}}, true); // extend X range in Run 4/5 + fRegistry.add("Event/trigger/hAnalysedTrigger", Form("analysed trigger %s;run number;counter", zorroGroup.cfg_swt_name.value.data()), o2::framework::HistType::kTH1D, {{100000, 500000.5, 600000.5}}, true); // extend X range in Run 4/5 + fRegistry.add("Event/trigger/hAnalysedToI", Form("analysed ToI %s;run number;counter", zorroGroup.cfg_swt_name.value.data()), o2::framework::HistType::kTH1D, {{100000, 500000.5, 600000.5}}, true); // extend X range in Run 4/5 + } + + if (doprocessNorm) { + fRegistry.addClone("Event/before/hCollisionCounter", "Event/norm/hCollisionCounter"); + fRegistry.add("Event/norm/hZvtx", "hZvtx;Z_{vtx} (cm)", o2::framework::HistType::kTH1D, {{100, -50, +50}}, false); + } + if (doprocessBC) { + auto hTVXCounter = fRegistry.add("BC/hTVXCounter", "TVX counter", o2::framework::HistType::kTH1D, {{6, -0.5f, 5.5f}}); + hTVXCounter->GetXaxis()->SetBinLabel(1, "TVX"); + hTVXCounter->GetXaxis()->SetBinLabel(2, "TVX && NoTFB"); + hTVXCounter->GetXaxis()->SetBinLabel(3, "TVX && NoITSROFB"); + hTVXCounter->GetXaxis()->SetBinLabel(4, "TVX && GoodRCT"); + hTVXCounter->GetXaxis()->SetBinLabel(5, "TVX && NoTFB && NoITSROFB"); + hTVXCounter->GetXaxis()->SetBinLabel(6, "TVX && NoTFB && NoITSROFB && GoodRCT"); + } + + mDCAFitter.setPropagateToPCA(dfGroup.propagateToPCA); + mDCAFitter.setMaxR(dfGroup.maxR); + mDCAFitter.setMaxDZIni(dfGroup.maxDZIni); + mDCAFitter.setMinParamChange(dfGroup.minParamChange); + mDCAFitter.setMinRelChi2Change(dfGroup.minRelChi2Change); + mDCAFitter.setUseAbsDCA(dfGroup.useAbsDCA); + mDCAFitter.setWeightedFinalPCA(dfGroup.useWeightedFinalPCA); + + mFwdDCAFitter.setPropagateToPCA(fdfGroup.propagateToPCA); + mFwdDCAFitter.setMaxR(fdfGroup.maxR); + mFwdDCAFitter.setMaxDXIni(fdfGroup.maxDXIni); + mFwdDCAFitter.setMinParamChange(fdfGroup.minParamChange); + mFwdDCAFitter.setMinRelChi2Change(fdfGroup.minRelChi2Change); + mFwdDCAFitter.setUseAbsDCA(fdfGroup.useAbsDCA); + } + + template + void initCCDB(TCollision const& collision) + { + if (mRunNumber == collision.runNumber()) { + return; + } + + auto run3grp_timestamp = collision.timestamp(); + o2::parameters::GRPObject* grpo = 0x0; + o2::parameters::GRPMagField* grpmag = 0x0; + if (grpo) { + o2::base::Propagator::initFieldFromGRP(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 << " kG"; + } 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; + } + o2::base::Propagator::initFieldFromGRP(grpmag); + // 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 << " kG"; + } + + mDCAFitter.setBz(d_bz); + mFwdDCAFitter.setBz(d_bz); + + // zorro + if constexpr (isTriggerAnalysis) { + zorro.setCCDBpath(zorroGroup.ccdbPathSoftwareTrigger); + zorro.setBCtolerance(zorroGroup.bcMarginForSoftwareTrigger); // this does nothing. + mToIidx = zorro.initCCDB(ccdb.service, collision.runNumber(), collision.timestamp(), zorroGroup.cfg_swt_name.value)[0]; + zorro.populateHistRegistry(fRegistry, collision.runNumber()); + + uint64_t nInspectedTVX = zorro.getInspectedTVX()->GetBinContent(1); + uint64_t nScalers = zorro.getScalers()->GetBinContent(zorro.getScalers()->GetXaxis()->FindBin(zorroGroup.cfg_swt_name.value.data())); + uint64_t nSelections = zorro.getSelections()->GetBinContent(zorro.getSelections()->GetXaxis()->FindBin(zorroGroup.cfg_swt_name.value.data())); + LOGF(info, "run number %d: total inspected TVX events = %llu, scalers = %llu, selections = %llu", collision.runNumber(), nInspectedTVX, nScalers, nSelections); + + fRegistry.fill(HIST("Event/trigger/hInspectedTVX"), collision.runNumber(), nInspectedTVX); + fRegistry.fill(HIST("Event/trigger/hScaler"), collision.runNumber(), nScalers); + fRegistry.fill(HIST("Event/trigger/hSelection"), collision.runNumber(), nSelections); + } + + mRunNumber = collision.runNumber(); + fDielectronCut.SetTrackPhiPositionRange(dielectroncuts.cfg_min_phiposition_track, dielectroncuts.cfg_max_phiposition_track, dielectroncuts.cfgRefR, d_bz, dielectroncuts.cfg_mirror_phi_track); + + auto grplhcif = ccdb->getForTimeStamp("GLO/Config/GRPLHCIF", collision.timestamp()); + int beamZ1 = grplhcif->getBeamZ(o2::constants::lhc::BeamC); + int beamZ2 = grplhcif->getBeamZ(o2::constants::lhc::BeamA); + int beamA1 = grplhcif->getBeamA(o2::constants::lhc::BeamC); + int beamA2 = grplhcif->getBeamA(o2::constants::lhc::BeamA); + beamE1 = grplhcif->getBeamEnergyPerNucleonInGeV(o2::constants::lhc::BeamC); + beamE2 = grplhcif->getBeamEnergyPerNucleonInGeV(o2::constants::lhc::BeamA); + beamM1 = o2::constants::physics::MassProton * beamA1; + beamM2 = o2::constants::physics::MassProton * beamA2; + beamP1 = std::sqrt(std::pow(beamE1, 2) - std::pow(beamM1, 2)); + beamP2 = std::sqrt(std::pow(beamE2, 2) - std::pow(beamM2, 2)); + LOGF(info, "beamZ1 = %d, beamZ2 = %d, beamA1 = %d, beamA2 = %d, beamE1 = %f (GeV), beamE2 = %f (GeV), beamM1 = %f (GeV), beamM2 = %f (GeV), beamP1 = %f (GeV), beamP2 = %f (GeV)", beamZ1, beamZ2, beamA1, beamA2, beamE1, beamE2, beamM1, beamM2, beamP1, beamP2); + + if (cfgApplySPresolution) { + auto list = ccdb->getForTimeStamp(spresoPath, collision.timestamp()); + h2sp_resolution = reinterpret_cast(list->FindObject(spresoHistName.value.data())); + LOGF(info, "h2sp_resolution.GetBinContent(40, 1) = %f", h2sp_resolution->GetBinContent(40, 1)); + } + } + + ~DileptonSV() + { + delete emh_pos; + emh_pos = 0x0; + delete emh_neg; + emh_neg = 0x0; + + used_trackIds_per_col.clear(); + used_trackIds_per_col.shrink_to_fit(); + map_mixed_eventId_to_globalBC.clear(); + + delete h2sp_resolution; + } + + void addhistograms() + { + const std::string qvec_det_names[7] = {"FT0M", "FT0A", "FT0C", "BTot", "BPos", "BNeg", "FV0A"}; + + std::string mass_axis_title = "m_{ll} (GeV/c^{2})"; + std::string pair_pt_axis_title = "p_{T,ll} (GeV/c)"; + std::string pair_dca_axis_title = "DCA_{ll} (#sigma)"; + std::string pair_y_axis_title = "y_{ll}"; + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + mass_axis_title = "m_{ee} (GeV/c^{2})"; + pair_pt_axis_title = "p_{T,ee} (GeV/c)"; + pair_y_axis_title = "y_{ee}"; + pair_dca_axis_title = "DCA_{ee}^{3D} (#sigma)"; + if (dielectroncuts.cfgDCAType == 1) { + pair_dca_axis_title = "DCA_{ee}^{XY} (#sigma)"; + } else if (dielectroncuts.cfgDCAType == 2) { + pair_dca_axis_title = "DCA_{ee}^{Z} (#sigma)"; + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + mass_axis_title = "m_{#mu#mu} (GeV/c^{2})"; + pair_pt_axis_title = "p_{T,#mu#mu} (GeV/c)"; + pair_dca_axis_title = "DCA_{#mu#mu}^{XY} (#sigma)"; + pair_y_axis_title = "y_{#mu#mu}"; + } + + // pair info + const o2::framework::AxisSpec axis_mass{ConfMllBins, mass_axis_title}; + const o2::framework::AxisSpec axis_pt{ConfPtllBins, pair_pt_axis_title}; + const o2::framework::AxisSpec axis_dca{ConfDCAllBins, pair_dca_axis_title}; + const o2::framework::AxisSpec axis_y{ConfYllBins, pair_y_axis_title}; + const o2::framework::AxisSpec axis_chi2PCA{ConfLog10Chi2PCABins, "log_{10}(#chi^{2}_{PCA})"}; + const o2::framework::AxisSpec axis_dl{ConfDLBins, "decay length (#sigma)"}; + + if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kQC)) { + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_chi2PCA, axis_dl}, true); + fRegistry.add("Pair/same/uls/hDeltaEtaDeltaPhi", "#Delta#eta-#Delta#varphi between 2 tracks;#Delta#varphi (rad.);#Delta#eta;", o2::framework::HistType::kTH2D, {{180, -M_PI, M_PI}, {400, -2, +2}}, true); + fRegistry.add("Pair/same/uls/hDeltaEtaDeltaPhiPos", Form("#Delta#eta-#Delta#varphi* between 2 tracks at r_{xy} = %3.2f m;#Delta#varphi* (rad.);#Delta#eta;", dielectroncuts.cfgRefR.value), o2::framework::HistType::kTH2D, {{180, -M_PI, M_PI}, {400, -2, +2}}, true); + + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + fRegistry.add("Pair/same/uls/hMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", o2::framework::HistType::kTH2D, {{90, 0, M_PI}, {100, 0.0f, 1.0f}}, true); // phiv is only for dielectron + } + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + fRegistry.addClone("Pair/same/", "Pair/mix/"); + + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kUPC)) { + const o2::framework::AxisSpec axis_aco{10, 0, 1.f, "#alpha = 1 - #frac{|#varphi_{l^{+}} - #varphi_{l^{-}}|}{#pi}"}; + const o2::framework::AxisSpec axis_asym_pt{10, 0, 1.f, "A = #frac{|p_{T,l^{+}} - p_{T,l^{-}}|}{|p_{T,l^{+}} + p_{T,l^{-}}|}"}; + const o2::framework::AxisSpec axis_dphi_l_ll{18, 0, M_PI, "#Delta#varphi = #varphi_{l} - #varphi_{ll} (rad.)"}; + + std::string frameName = "CS"; + if (cfgPolarizationFrame == 0) { + frameName = "CS"; + } else if (cfgPolarizationFrame == 1) { + frameName = "HX"; + } else { + LOG(fatal) << "set 0 or 1 to cfgPolarizationFrame!"; + } + const o2::framework::AxisSpec axis_cos_theta{ConfPolarizationCosThetaBins, Form("cos(#theta^{%s})", frameName.data())}; + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_aco, axis_asym_pt, axis_dphi_l_ll, axis_cos_theta}, true); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + fRegistry.addClone("Pair/same/", "Pair/mix/"); + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV2SP)) { + if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV2SP)) { + nmod = 2; + } + + const o2::framework::AxisSpec axis_sp{ConfSPBins, Form("#vec{u}_{%d,ll} #upoint #vec{Q}_{%d}^{%s}", nmod, nmod, qvec_det_names[cfgQvecEstimator].data())}; + + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_sp}, true); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + + fRegistry.add("Pair/mix/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true); + fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lspp/"); + fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lsmm/"); + + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kPolarization)) { + std::string frameName = "CS"; + if (cfgPolarizationFrame == 0) { + frameName = "CS"; + } else if (cfgPolarizationFrame == 1) { + frameName = "HX"; + } else { + LOG(fatal) << "set 0 or 1 to cfgPolarizationFrame!"; + } + + const o2::framework::AxisSpec axis_cos_theta{ConfPolarizationCosThetaBins, Form("cos(#theta^{%s})", frameName.data())}; + const o2::framework::AxisSpec axis_phi{ConfPolarizationPhiBins, Form("#varphi^{%s} (rad.)", frameName.data())}; + const o2::framework::AxisSpec axis_quadmom{ConfPolarizationQuadMomBins, Form("#frac{3 cos^{2}(#theta^{%s}) -1}{2}", frameName.data())}; + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_cos_theta, axis_phi, axis_quadmom}, true); + + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + fRegistry.addClone("Pair/same/", "Pair/mix/"); + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kHFll)) { + const o2::framework::AxisSpec axis_dphi_ll{36, -M_PI / 2., 3. / 2. * M_PI, "#Delta#varphi = #varphi_{l1} - #varphi_{l2} (rad.)"}; // for kHFll + const o2::framework::AxisSpec axis_deta_ll{40, -2., 2., "#Delta#eta = #eta_{l1} - #eta_{l2}"}; + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_dphi_ll, axis_deta_ll}, true); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + fRegistry.addClone("Pair/same/", "Pair/mix/"); + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) { + nmod = 2; + const o2::framework::AxisSpec axis_sp{ConfSPBins, Form("#vec{u}_{%d,ll} #upoint #vec{Q}_{%d}^{%s}", nmod, nmod, qvec_det_names[cfgQvecEstimator].data())}; + const o2::framework::AxisSpec axis_bootstrap{cfgNumBootstrapSamples, 0.5, static_cast(cfgNumBootstrapSamples) + 0.5, "sample"}; // for bootstrap samples + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_sp, axis_bootstrap}, true); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + fRegistry.add("Pair/mix/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true); + fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lspp/"); + fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lsmm/"); + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistogramsBootstrap(&fRegistry, cfgNumBootstrapSamples); + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV2EP)) { + if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV2EP)) { + nmod = 2; + } + + // similar to Eq.1 in PRL 126, 162001 (2021) + const o2::framework::AxisSpec axis_ep{4, 0, 2 * M_PI / nmod, Form("#varphi_{ll} #minus #Psi_{%d}^{%s}", nmod, qvec_det_names[cfgQvecEstimator].data())}; + + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_ep}, true); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + + fRegistry.add("Pair/mix/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true); + fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lspp/"); + fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lsmm/"); + } else { // same as kQC to avoid seg. fault + fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/"); + fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/"); + fRegistry.addClone("Pair/same/", "Pair/mix/"); + } + + // event info + if (nmod == 2) { + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<2>(&fRegistry); + } else { + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<-1>(&fRegistry); + } + fRegistry.add("Event/before/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), o2::framework::HistType::kTH2F, {{110, 0, 110}, {180, -M_PI_2, +M_PI_2}}, false); + fRegistry.add("Event/after/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), o2::framework::HistType::kTH2F, {{110, 0, 110}, {180, -M_PI_2, +M_PI_2}}, false); + } + + void DefineEMEventCut() + { + fEMEventCut = EMEventCut("fEMEventCut", "fEMEventCut"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(eventcuts.cfgZvtxMin, eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireVertexTOFmatched(eventcuts.cfgRequireVertexTOFmatched); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetRequireNoCollInTimeRangeStandard(eventcuts.cfgRequireNoCollInTimeRangeStandard); + fEMEventCut.SetRequireNoCollInTimeRangeStrict(eventcuts.cfgRequireNoCollInTimeRangeStrict); + fEMEventCut.SetRequireNoCollInITSROFStandard(eventcuts.cfgRequireNoCollInITSROFStandard); + fEMEventCut.SetRequireNoCollInITSROFStrict(eventcuts.cfgRequireNoCollInITSROFStrict); + fEMEventCut.SetRequireNoHighMultCollInPrevRof(eventcuts.cfgRequireNoHighMultCollInPrevRof); + fEMEventCut.SetRequireGoodITSLayer3(eventcuts.cfgRequireGoodITSLayer3); + fEMEventCut.SetRequireGoodITSLayer0123(eventcuts.cfgRequireGoodITSLayer0123); + fEMEventCut.SetRequireGoodITSLayersAll(eventcuts.cfgRequireGoodITSLayersAll); + } + + void DefineDielectronCut() + { + fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); + + // for pair + fDielectronCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); + fDielectronCut.SetPairPtRange(dielectroncuts.cfg_min_pair_pt, dielectroncuts.cfg_max_pair_pt); + fDielectronCut.SetPairYRange(dielectroncuts.cfg_min_pair_y, dielectroncuts.cfg_max_pair_y); + fDielectronCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma + fDielectronCut.SetMaxMeePhiVDep([&](float phiv) { return dielectroncuts.cfg_phiv_intercept + phiv * dielectroncuts.cfg_phiv_slope; }, dielectroncuts.cfg_min_phiv, dielectroncuts.cfg_max_phiv); + fDielectronCut.ApplyPhiV(dielectroncuts.cfg_apply_phiv); + fDielectronCut.SetMindEtadPhi(dielectroncuts.cfg_apply_detadphi, dielectroncuts.cfg_apply_detadphiPosition, dielectroncuts.cfg_min_deta, dielectroncuts.cfg_min_dphi); + fDielectronCut.SetPairOpAng(dielectroncuts.cfg_min_opang, dielectroncuts.cfg_max_opang); + // fDielectronCut.SetRequireDifferentSides(dielectroncuts.cfg_require_diff_sides); + + if (dielectroncuts.cfg_apply_detadphi && dielectroncuts.cfg_apply_detadphiPosition) { + LOG(fatal) << "cfg_apply_detadphi and cfg_apply_detadphiPosition cannot be used simultaneously. Please select only one."; + } + + // for track + fDielectronCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, dielectroncuts.cfg_max_pt_track); + fDielectronCut.SetTrackEtaRange(dielectroncuts.cfg_min_eta_track, dielectroncuts.cfg_max_eta_track); + fDielectronCut.SetTrackPhiRange(dielectroncuts.cfg_min_phi_track, dielectroncuts.cfg_max_phi_track, dielectroncuts.cfg_mirror_phi_track, dielectroncuts.cfg_reject_phi_track); + fDielectronCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); + fDielectronCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); + fDielectronCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fDielectronCut.SetMaxFracSharedClustersTPC(dielectroncuts.cfg_max_frac_shared_clusters_tpc); + fDielectronCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); + fDielectronCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); + fDielectronCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); + fDielectronCut.SetMeanClusterSizeITS(dielectroncuts.cfg_min_its_cluster_size, dielectroncuts.cfg_max_its_cluster_size); + fDielectronCut.SetTrackMaxDcaSigma(dielectroncuts.cfg_max_dca_sigma, dielectroncuts.cfgDCAType); + fDielectronCut.SetTrackMaxDcaXY(dielectroncuts.cfg_max_dcaxy); + fDielectronCut.SetTrackMaxDcaZ(dielectroncuts.cfg_max_dcaz); + fDielectronCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); + fDielectronCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); + fDielectronCut.SetChi2TOF(0, dielectroncuts.cfg_max_chi2tof); + // fDielectronCut.SetRelDiffPin(dielectroncuts.cfg_min_rel_diff_pin, dielectroncuts.cfg_max_rel_diff_pin); + fDielectronCut.EnableTTCA(dielectroncuts.enableTTCA); + + // for eID + fDielectronCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); + fDielectronCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); + // fDielectronCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); + fDielectronCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); + fDielectronCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); + fDielectronCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); + fDielectronCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); + fDielectronCut.SetPinRangeForPionRejectionTPC(dielectroncuts.cfg_min_pin_pirejTPC, dielectroncuts.cfg_max_pin_pirejTPC); + + if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { // please call this at the end of DefineDileptonCut + std::vector binsML{}; + binsML.reserve(dielectroncuts.binsMLPID.value.size()); + for (size_t i = 0; i < dielectroncuts.binsMLPID.value.size(); i++) { + binsML.emplace_back(dielectroncuts.binsMLPID.value[i]); + } + std::vector thresholdsML{}; + thresholdsML.reserve(dielectroncuts.cutsMLPID.value.size()); + for (size_t i = 0; i < dielectroncuts.cutsMLPID.value.size(); i++) { + thresholdsML.emplace_back(dielectroncuts.cutsMLPID.value[i]); + } + fDielectronCut.SetMLThresholds(binsML, thresholdsML); + } // end of PID ML + } + + void DefineDimuonCut() + { + fDimuonCut = DimuonCut("fDimuonCut", "fDimuonCut"); + + // for pair + fDimuonCut.SetMassRange(dimuoncuts.cfg_min_mass, dimuoncuts.cfg_max_mass); + fDimuonCut.SetPairPtRange(dimuoncuts.cfg_min_pair_pt, dimuoncuts.cfg_max_pair_pt); + fDimuonCut.SetPairYRange(dimuoncuts.cfg_min_pair_y, dimuoncuts.cfg_max_pair_y); + fDimuonCut.SetPairDCAxyRange(dimuoncuts.cfg_min_pair_dcaxy, dimuoncuts.cfg_max_pair_dcaxy); + fDimuonCut.SetMindEtadPhi(dimuoncuts.cfg_apply_detadphi, dimuoncuts.cfg_min_deta, dimuoncuts.cfg_min_dphi); + + // for track + fDimuonCut.SetTrackType(dimuoncuts.cfg_track_type); + fDimuonCut.SetTrackPtRange(dimuoncuts.cfg_min_pt_track, dimuoncuts.cfg_max_pt_track); + fDimuonCut.SetTrackEtaRange(dimuoncuts.cfg_min_eta_track, dimuoncuts.cfg_max_eta_track); + fDimuonCut.SetTrackPhiRange(dimuoncuts.cfg_min_phi_track, dimuoncuts.cfg_max_phi_track); + fDimuonCut.SetNClustersMFT(dimuoncuts.cfg_min_ncluster_mft, 10); + fDimuonCut.SetNClustersMCHMID(dimuoncuts.cfg_min_ncluster_mch, 20); + fDimuonCut.SetChi2(0.f, dimuoncuts.cfg_max_chi2); + fDimuonCut.SetChi2MFT(0.f, dimuoncuts.cfg_max_chi2mft); + // fDimuonCut.SetMatchingChi2MCHMFT(0.f, dimuoncuts.cfg_max_matching_chi2_mftmch); + fDimuonCut.SetMaxMatchingChi2MCHMFTPtDep([&](float pt) { return (pt < dimuoncuts.cfg_border_pt_for_chi2mchmft ? dimuoncuts.cfg_max_matching_chi2_mftmch_lowPt : dimuoncuts.cfg_max_matching_chi2_mftmch_highPt); }); + fDimuonCut.SetMatchingChi2MCHMID(0.f, dimuoncuts.cfg_max_matching_chi2_mchmid); + fDimuonCut.SetDCAxy(0.f, dimuoncuts.cfg_max_dcaxy); + fDimuonCut.SetRabs(dimuoncuts.cfg_min_rabs, dimuoncuts.cfg_max_rabs); + fDimuonCut.SetMaxPDCARabsDep([&](float rabs) { return (rabs < 26.5 ? 594.f : 324.f); }); + fDimuonCut.SetMaxdPtdEtadPhiwrtMCHMID(dimuoncuts.cfg_max_relDPt_wrt_matchedMCHMID, dimuoncuts.cfg_max_DEta_wrt_matchedMCHMID, dimuoncuts.cfg_max_DPhi_wrt_matchedMCHMID); // this is relevant for global muons + fDimuonCut.SetMFTHitMap(dimuoncuts.requireMFTHitMap, dimuoncuts.requiredMFTDisks); + fDimuonCut.EnableTTCA(dimuoncuts.enableTTCA); + } + + template + bool foundHFSV(TTrack const& track) + { + int ptbin = lower_bound(dielectroncuts.binsMLSCT.value.begin(), dielectroncuts.binsMLSCT.value.end(), track.pt()) - dielectroncuts.binsMLSCT.value.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(dielectroncuts.binsMLSCT.value.size()) - 2 < ptbin) { + ptbin = static_cast(dielectroncuts.binsMLSCT.value.size()) - 2; + } + + for (int i = 0; i < static_cast(track.nSV()); i++) { + auto probaSCT = track.probaSCT(i); + // LOGF(info, "track.globalIndex() = %d, pt = %f, i = %d, probaSCT[0] = %f, probaSCT[1] = %f, probaSCT[2] = %f, probaSCT[3] = %f", track.globalIndex(), track.pt(), i, probaSCT[0], probaSCT[1], probaSCT[2], probaSCT[3]); + if (probaSCT[1] > dielectroncuts.cutsMLSCTeT_prompt_hc.value[ptbin] || probaSCT[2] > dielectroncuts.cutsMLSCTeT_nonprompt_hc.value[ptbin] || probaSCT[3] > dielectroncuts.cutsMLSCTeT_hb.value[ptbin]) { + return true; + } + } + return false; + } + + template + bool isGoodQvector(TQvectors const& qvectors) + { + bool is_good = true; + for (const auto& qn : qvectors[nmod]) { + if (std::fabs(qn[0]) > 20.f || std::fabs(qn[1]) > 20.f) { + is_good = false; + break; + } + } + return is_good; + } + + float getSPresolution(const float centrality, const int occupancy) + { + if (h2sp_resolution == nullptr) { + return 1.f; + } else { + int binId_cen = h2sp_resolution->GetXaxis()->FindBin(centrality); + int binId_occ = h2sp_resolution->GetYaxis()->FindBin(occupancy); + + if (centrality < h2sp_resolution->GetXaxis()->GetXmin()) { + binId_cen = 1; + } + if (h2sp_resolution->GetXaxis()->GetXmax() < centrality) { + binId_cen = h2sp_resolution->GetXaxis()->GetNbins(); + } + + if (occupancy < h2sp_resolution->GetYaxis()->GetXmin()) { + binId_occ = 1; + } + if (h2sp_resolution->GetYaxis()->GetXmax() < occupancy) { + binId_occ = h2sp_resolution->GetYaxis()->GetNbins(); + } + return h2sp_resolution->GetBinContent(binId_cen, binId_occ); + } + } + + struct dileptonSV { + float pt1{999}; + float eta1{999}; + float phi1{999}; + + float pt2{999}; + float eta2{999}; + float phi2{999}; + + // float mass{999}; + // float pt{999}; + // float rapidity{999}; + + float chi2PCA{999}; + float cpa{999}; + float lxy{999}; + float lz{999}; + float lxyz{999}; + float lxyErr{999}; + float lzErr{999}; + float lxyzErr{999}; + }; + + template + bool findSV(TCollision const& collision, TTrackParCov const& t1, TTrackParCov const& t2, dileptonSV& candidate) + { + try { + if (mDCAFitter.process(t1, t2) == 0) { + return false; + } + } catch (const std::runtime_error& error) { + LOG(info) << "Run time error found: " << error.what() << ". DCAFitterN cannot work, skipping the candidate."; + return false; + } + candidate.chi2PCA = mDCAFitter.getChi2AtPCACandidate(); + if (dfGroup.maxLog10Chi2PCA < std::log10(candidate.chi2PCA)) { + return false; + } + + mDCAFitter.propagateTracksToVertex(); + const auto& secondaryVertex = mDCAFitter.getPCACandidate(); + auto trackParCov0 = mDCAFitter.getTrack(0); + auto trackParCov1 = mDCAFitter.getTrack(1); + + // get track momenta + std::array pvertex = {collision.posX(), collision.posY(), collision.posZ()}; + std::array pvec0{}; + std::array pvec1{}; + trackParCov0.getPxPyPzGlo(pvec0); + trackParCov1.getPxPyPzGlo(pvec1); + std::array momDilepton = {pvec0[0] + pvec1[0], pvec0[1] + pvec1[1], pvec0[2] + pvec1[2]}; + candidate.cpa = RecoDecay::cpa(pvertex, secondaryVertex, momDilepton); + + candidate.lxy = std::sqrt(std::pow(secondaryVertex[0] - collision.posX(), 2) + std::pow(secondaryVertex[1] - collision.posY(), 2)); + candidate.lz = secondaryVertex[2] - collision.posZ(); + candidate.lxyz = std::sqrt(std::pow(candidate.lxy, 2) + std::pow(candidate.lz, 2)); + + auto primaryVertex = getPrimaryVertex(collision); + std::array covVtx = mDCAFitter.calcPCACovMatrixFlat(); + double phi{0}, theta{0}; + getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); + candidate.lxyErr = std::sqrt(getRotatedCovMatrixXX(primaryVertex.getCov(), phi, 0.) + getRotatedCovMatrixXX(covVtx, phi, 0.)); + candidate.lzErr = std::sqrt(getRotatedCovMatrixXX(primaryVertex.getCov(), 0, theta) + getRotatedCovMatrixXX(covVtx, 0, theta)); + candidate.lxyzErr = std::sqrt(getRotatedCovMatrixXX(primaryVertex.getCov(), phi, theta) + getRotatedCovMatrixXX(covVtx, phi, theta)); + + candidate.pt1 = trackParCov0.getPt(); + candidate.eta1 = trackParCov0.getEta(); + candidate.phi1 = RecoDecay::constrainAngle(trackParCov0.getPhi(), 0, 1U); // 0-2pi + + candidate.pt2 = trackParCov1.getPt(); + candidate.eta2 = trackParCov1.getEta(); + candidate.phi2 = RecoDecay::constrainAngle(trackParCov1.getPhi(), 0, 1U); // 0-2pi + + return true; + } + + template + bool findSVFwd(TCollision const& collision, TTrackParCov const& t1, TTrackParCov const& t2, dileptonSV& candidate) + { + try { + if (mFwdDCAFitter.process(t1, t2) == 0) { + return false; + } + } catch (const std::runtime_error& error) { + LOG(info) << "Run time error found: " << error.what() << ". FwdDCAFitterN cannot work, skipping the candidate."; + return false; + } + mFwdDCAFitter.FwdpropagateTracksToVertex(); + + std::array pvertex = {collision.posX(), collision.posY(), collision.posZ()}; + const auto& secondaryVertex = mFwdDCAFitter.getPCACandidate(); + candidate.chi2PCA = mFwdDCAFitter.getChi2AtPCACandidate(); + if (fdfGroup.maxLog10Chi2PCA < std::log10(candidate.chi2PCA)) { + return false; + } + auto trackParCov0 = mFwdDCAFitter.getTrack(0); + auto trackParCov1 = mFwdDCAFitter.getTrack(1); + std::array momDilepton = {static_cast(trackParCov0.getPx() + trackParCov1.getPx()), static_cast(trackParCov0.getPy() + trackParCov1.getPy()), static_cast(trackParCov0.getPz() + trackParCov1.getPz())}; + candidate.cpa = RecoDecay::cpa(pvertex, secondaryVertex, momDilepton); + + candidate.lxy = std::sqrt(std::pow(secondaryVertex[0] - collision.posX(), 2) + std::pow(secondaryVertex[1] - collision.posY(), 2)); + candidate.lz = secondaryVertex[2] - collision.posZ(); + candidate.lxyz = std::sqrt(std::pow(candidate.lxy, 2) + std::pow(candidate.lz, 2)); + + auto primaryVertex = getPrimaryVertex(collision); + std::array covVtx = mFwdDCAFitter.calcPCACovMatrixFlat(); + double phi{0}, theta{0}; + getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); + candidate.lxyErr = std::sqrt(getRotatedCovMatrixXX(primaryVertex.getCov(), phi, 0.) + getRotatedCovMatrixXX(covVtx, phi, 0.)); + candidate.lzErr = std::sqrt(getRotatedCovMatrixXX(primaryVertex.getCov(), 0, theta) + getRotatedCovMatrixXX(covVtx, 0, theta)); + candidate.lxyzErr = std::sqrt(getRotatedCovMatrixXX(primaryVertex.getCov(), phi, theta) + getRotatedCovMatrixXX(covVtx, phi, theta)); + + candidate.pt1 = trackParCov0.getPt(); + candidate.eta1 = trackParCov0.getEta(); + candidate.phi1 = RecoDecay::constrainAngle(trackParCov0.getPhi(), 0, 1U); // 0-2pi + + candidate.pt2 = trackParCov1.getPt(); + candidate.eta2 = trackParCov1.getEta(); + candidate.phi2 = RecoDecay::constrainAngle(trackParCov1.getPhi(), 0, 1U); // 0-2pi + + return true; + } + + template + bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&, const std::vector weightvector) + { + dileptonSV candidate; + if constexpr (ev_id == 0) { + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { + if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { + return false; + } + } else { // cut-based + if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { + return false; + } + } + auto trackParCov1 = getTrackParCov(t1); + auto trackParCov2 = getTrackParCov(t2); + if (!findSV(collision, trackParCov1, trackParCov2, candidate)) { + return false; + } + if constexpr (withSCT) { + if (foundHFSV(t1) || foundHFSV(t2)) { + return false; + } + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { + return false; + } + if (!map_best_match_globalmuon[t1.globalIndex()] || !map_best_match_globalmuon[t2.globalIndex()]) { + return false; + } + auto trackParCov1 = o2::aod::fwdtrackutils::getTrackParCovFwd(t1, t1); + auto trackParCov2 = o2::aod::fwdtrackutils::getTrackParCovFwd(t2, t2); + if (!findSVFwd(collision, trackParCov1, trackParCov2, candidate)) { + return false; + } + } + } else if constexpr (ev_id == 1) { // for mixed event + candidate.pt1 = t1.pt(); + candidate.eta1 = t1.eta(); + candidate.phi1 = t1.phi(); + candidate.pt2 = t2.pt(); + candidate.eta2 = t2.eta(); + candidate.phi2 = t2.phi(); + } + + // if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + // auto trackParCov1 = getTrackParCov(t1); + // auto trackParCov2 = getTrackParCov(t2); + // if(!findSV(collision, trackParCov1, trackParCov2, candidate)){ + // return false; + // } + // } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + // auto trackParCov1 = o2::aod::fwdtrackutils::getTrackParCovFwd(t1, t1); + // auto trackParCov2 = o2::aod::fwdtrackutils::getTrackParCovFwd(t2, t2); + // if(!findSVFwd(collision, trackParCov1, trackParCov2, candidate)){ + // return false; + // } + // } + + // if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + // if (!cut.IsSelectedPair(t1, t2)) { + // return false; + // } + // } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + // if (!cut.IsSelectedPair(t1, t2)) { + // return false; + // } + // } + + float weight = 1.f; + if constexpr (ev_id == 0) { + if (cfgApplyWeightTTCA) { + weight = map_weight[std::make_pair(t1.globalIndex(), t2.globalIndex())]; + } + // LOGF(info, "ev_id = %d, t1.sign() = %d, t2.sign() = %d, map_weight[std::make_pair(%d, %d)] = %f", ev_id, t1.sign(), t2.sign(), t1.globalIndex(), t2.globalIndex(), weight); + } + + // ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), leptonM1); + // ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), leptonM2); + ROOT::Math::PtEtaPhiMVector v1(candidate.pt1, candidate.eta1, candidate.phi1, leptonM1); + ROOT::Math::PtEtaPhiMVector v2(candidate.pt2, candidate.eta2, candidate.phi2, leptonM2); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + float pair_dca = 999.f; + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + pair_dca = o2::aod::pwgem::dilepton::utils::pairutil::pairDCAQuadSum(o2::aod::pwgem::dilepton::utils::emtrackutil::dca3DinSigma(t1), o2::aod::pwgem::dilepton::utils::emtrackutil::dca3DinSigma(t2)); + if (dielectroncuts.cfgDCAType == 1) { + pair_dca = o2::aod::pwgem::dilepton::utils::pairutil::pairDCAQuadSum(o2::aod::pwgem::dilepton::utils::emtrackutil::dcaXYinSigma(t1), o2::aod::pwgem::dilepton::utils::emtrackutil::dcaXYinSigma(t2)); + } else if (dielectroncuts.cfgDCAType == 2) { + pair_dca = o2::aod::pwgem::dilepton::utils::pairutil::pairDCAQuadSum(o2::aod::pwgem::dilepton::utils::emtrackutil::dcaZinSigma(t1), o2::aod::pwgem::dilepton::utils::emtrackutil::dcaZinSigma(t2)); + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + pair_dca = o2::aod::pwgem::dilepton::utils::pairutil::pairDCAQuadSum(o2::aod::pwgem::dilepton::utils::emtrackutil::fwdDcaXYinSigma(t1), o2::aod::pwgem::dilepton::utils::emtrackutil::fwdDcaXYinSigma(t2)); + } + + if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kQC)) { + float deta = t1.sign() * v1.Pt() > t2.sign() * v2.Pt() ? v1.Eta() - v2.Eta() : v2.Eta() - v1.Eta(); + float dphi = t1.sign() * v1.Pt() > t2.sign() * v2.Pt() ? v1.Phi() - v2.Phi() : v2.Phi() - v1.Phi(); + dphi = RecoDecay::constrainAngle(dphi, -M_PI, 1U); // -pi - +pi + + float phiPosition1 = RecoDecay::constrainAngle(t1.phi() + std::asin(t1.sign() * -0.30282 * (d_bz * 0.1) * dielectroncuts.cfgRefR / (2.f * t1.pt())), 0, 1U); // 0-2pi + float phiPosition2 = RecoDecay::constrainAngle(t2.phi() + std::asin(t2.sign() * -0.30282 * (d_bz * 0.1) * dielectroncuts.cfgRefR / (2.f * t2.pt())), 0, 1U); // 0-2pi + float dphiPosition = t1.sign() * v1.Pt() > t2.sign() * v2.Pt() ? phiPosition1 - phiPosition2 : phiPosition2 - phiPosition1; + dphiPosition = RecoDecay::constrainAngle(dphiPosition, -M_PI, 1U); // -pi - +pi + + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), d_bz); + // float opAng = o2::aod::pwgem::dilepton::utils::pairutil::getOpeningAngle(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz()); + + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), std::log10(candidate.chi2PCA), candidate.lxyz / candidate.lxyzErr, weight); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hDeltaEtaDeltaPhi"), dphi, deta, weight); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hDeltaEtaDeltaPhiPos"), dphiPosition, deta, weight); + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hMvsPhiV"), phiv, v12.M(), weight); + } + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), std::log10(candidate.chi2PCA), candidate.lxyz / candidate.lxyzErr, weight); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hDeltaEtaDeltaPhi"), dphi, deta, weight); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hDeltaEtaDeltaPhiPos"), dphiPosition, deta, weight); + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hMvsPhiV"), phiv, v12.M(), weight); + } + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), std::log10(candidate.chi2PCA), candidate.lxyz / candidate.lxyzErr, weight); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hDeltaEtaDeltaPhi"), dphi, deta, weight); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hDeltaEtaDeltaPhiPos"), dphiPosition, deta, weight); + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hMvsPhiV"), phiv, v12.M(), weight); + } + } + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kUPC)) { + float dphi = v1.Phi() - v2.Phi(); + o2::math_utils::bringToPMPi(dphi); + float aco = 1.f - std::fabs(dphi) / M_PI; + float asym = std::fabs(v1.Pt() - v2.Pt()) / (v1.Pt() + v2.Pt()); + float dphi_l_ll = v1.Phi() - v12.Phi(); + o2::math_utils::bringToPMPi(dphi_l_ll); + + float cos_thetaPol = 999, phiPol = 999.f; + auto arrM = std::array{static_cast(v12.Px()), static_cast(v12.Py()), static_cast(v12.Pz()), static_cast(v12.M())}; + auto arrD = t1.sign() * t1.pt() > t2.sign() * t2.pt() ? std::array{t1.px(), t1.py(), t1.pz(), leptonM1} : std::array{t2.px(), t2.py(), t2.pz(), leptonM2}; + + if (cfgPolarizationFrame == 0) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } else if (cfgPolarizationFrame == 1) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } + + o2::math_utils::bringToPMPi(phiPol); + + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), aco, asym, std::fabs(dphi_l_ll), cos_thetaPol, weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), aco, asym, std::fabs(dphi_l_ll), cos_thetaPol, weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), aco, asym, std::fabs(dphi_l_ll), cos_thetaPol, weight); + } + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV2SP)) { + std::array q2ft0m = {collision.q2xft0m(), collision.q2yft0m()}; + std::array q2ft0a = {collision.q2xft0a(), collision.q2yft0a()}; + std::array q2ft0c = {collision.q2xft0c(), collision.q2yft0c()}; + std::array q2btot = {collision.q2xbtot(), collision.q2ybtot()}; + std::array q2bpos = {collision.q2xbpos(), collision.q2ybpos()}; + std::array q2bneg = {collision.q2xbneg(), collision.q2ybneg()}; + std::array q2fv0a = {collision.q2xfv0a(), collision.q2yfv0a()}; + + std::vector>> qvectors = { + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 0th harmonics + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 1st harmonics + {q2ft0m, q2ft0a, q2ft0c, q2btot, q2bpos, q2bneg, q2fv0a}, // 2nd harmonics + }; + + if constexpr (ev_id == 0) { + // LOGF(info, "collision.centFT0C() = %f, collision.trackOccupancyInTimeRange() = %d, getSPresolution = %f", collision.centFT0C(), collision.trackOccupancyInTimeRange(), getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange())); + + float sp = RecoDecay::dotProd(std::array{static_cast(std::cos(nmod * v12.Phi())), static_cast(std::sin(nmod * v12.Phi()))}, qvectors[nmod][cfgQvecEstimator]) / getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange()); + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, weight); + } + } else if constexpr (ev_id == 1) { + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } + } + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kPolarization)) { + float cos_thetaPol = 999, phiPol = 999.f; + auto arrM = std::array{static_cast(v12.Px()), static_cast(v12.Py()), static_cast(v12.Pz()), static_cast(v12.M())}; + auto random_sign = std::pow(-1, engine() % 2); // -1^0 = +1 or -1^1 = -1; + auto arrD = t1.sign() * t2.sign() < 0 ? (t1.sign() > 0 ? std::array{t1.px(), t1.py(), t1.pz(), leptonM1} : std::array{t2.px(), t2.py(), t2.pz(), leptonM2}) : (random_sign > 0 ? std::array{t1.px(), t1.py(), t1.pz(), leptonM1} : std::array{t2.px(), t2.py(), t2.pz(), leptonM2}); + + if (cfgPolarizationFrame == 0) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } else if (cfgPolarizationFrame == 1) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } + + o2::math_utils::bringToPMPi(phiPol); + float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; + + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kHFll)) { + float dphi = v1.Phi() - v2.Phi(); + dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf); + float deta = v1.Eta() - v2.Eta(); + + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi, deta, weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi, deta, weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi, deta, weight); + } + + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) { + std::array q2ft0m = {collision.q2xft0m(), collision.q2yft0m()}; + std::array q2ft0a = {collision.q2xft0a(), collision.q2yft0a()}; + std::array q2ft0c = {collision.q2xft0c(), collision.q2yft0c()}; + std::array q2btot = {collision.q2xbtot(), collision.q2ybtot()}; + std::array q2bpos = {collision.q2xbpos(), collision.q2ybpos()}; + std::array q2bneg = {collision.q2xbneg(), collision.q2ybneg()}; + std::array q2fv0a = {collision.q2xfv0a(), collision.q2yfv0a()}; + + std::vector>> qvectors = { + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 0th harmonics + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 1st harmonics + {q2ft0m, q2ft0a, q2ft0c, q2btot, q2bpos, q2bneg, q2fv0a}, // 2nd harmonics + }; + + if constexpr (ev_id == 0) { + // LOGF(info, "collision.centFT0C() = %f, collision.trackOccupancyInTimeRange() = %d, getSPresolution = %f", collision.centFT0C(), collision.trackOccupancyInTimeRange(), getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange())); + + float sp = RecoDecay::dotProd(std::array{static_cast(std::cos(nmod * v12.Phi())), static_cast(std::sin(nmod * v12.Phi()))}, qvectors[nmod][cfgQvecEstimator]) / getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange()); + for (int i = 0; i < cfgNumBootstrapSamples; i++) { + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, i + 0.5, weightvector.at(i)); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, i + 0.5, weightvector.at(i)); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, i + 0.5, weightvector.at(i)); + } + } + } else if constexpr (ev_id == 1) { + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } + } + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV2EP)) { + if constexpr (ev_id == 0) { + std::vector> eventplanes = { + {0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f}, // 0th harmonics + {0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f}, // 1st harmonics + {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg(), collision.ep2fv0a()}, // 2nd harmonics + }; + + float ep = RecoDecay::constrainAngle(eventplanes[nmod][cfgQvecEstimator], -o2::constants::math::PI / nmod, static_cast(nmod)); + float phi_ll = RecoDecay::constrainAngle(v12.Phi(), -o2::constants::math::PI / nmod, 1U); + float dphi_ll_ep = std::fabs(RecoDecay::constrainAngle(phi_ll - ep, 0.f, static_cast(nmod))); + + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi_ll_ep, weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi_ll_ep, weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi_ll_ep, weight); + } + } else if constexpr (ev_id == 1) { + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } + } + } else { // same as kQC to avoid seg. fault + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight); + } + } + + // store tracks for event mixing without double counting + if constexpr (ev_id == 0) { + std::pair key_df_collision = std::make_pair(ndf, collision.globalIndex()); + + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + if (std::find(used_trackIds_per_col.begin(), used_trackIds_per_col.end(), t1.globalIndex()) == used_trackIds_per_col.end()) { + used_trackIds_per_col.emplace_back(t1.globalIndex()); + if (cfgDoMix) { + if (t1.sign() > 0) { + emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrack(candidate.pt1, candidate.eta1, candidate.phi1, leptonM1, t1.sign(), t1.dcaXY(), t1.dcaZ(), t1.cYY(), t1.cZY(), t1.cZZ())); + // emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrackWithCov(t1.pt(), t1.eta(), t1.phi(), leptonM1, t1.sign(), t1.dcaXY(), t1.dcaZ(), t1.cYY(), t1.cZY(), t1.cZZ(), + // t1.x(), t1.y(), t1.z(), t1.alpha(), t1.snp(), t1.tgl(), + // t1.cSnpY(), t1.cSnpZ(), t1.cSnpSnp(), + // t1.cTglY(), t1.cTglZ(), t1.cTglSnp(), t1.cTglTgl(), + // t1.c1PtY(), t1.c1PtZ(), t1.c1PtSnp(), t1.c1PtTgl(), t1.c1Pt21Pt2())); + } else { + emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrack(candidate.pt1, candidate.eta1, candidate.phi1, leptonM1, t1.sign(), t1.dcaXY(), t1.dcaZ(), t1.cYY(), t1.cZY(), t1.cZZ())); + // emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrackWithCov(t1.pt(), t1.eta(), t1.phi(), leptonM1, t1.sign(), t1.dcaXY(), t1.dcaZ(), t1.cYY(), t1.cZY(), t1.cZZ(), + // t1.x(), t1.y(), t1.z(), t1.alpha(), t1.snp(), t1.tgl(), + // t1.cSnpY(), t1.cSnpZ(), t1.cSnpSnp(), + // t1.cTglY(), t1.cTglZ(), t1.cTglSnp(), t1.cTglTgl(), + // t1.c1PtY(), t1.c1PtZ(), t1.c1PtSnp(), t1.c1PtTgl(), t1.c1Pt21Pt2())); + } + } + } + if (std::find(used_trackIds_per_col.begin(), used_trackIds_per_col.end(), t2.globalIndex()) == used_trackIds_per_col.end()) { + used_trackIds_per_col.emplace_back(t2.globalIndex()); + if (cfgDoMix) { + if (t2.sign() > 0) { + emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrack(candidate.pt2, candidate.eta2, candidate.phi2, leptonM2, t2.sign(), t2.dcaXY(), t2.dcaZ(), t2.cYY(), t2.cZY(), t2.cZZ())); + // emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrackWithCov(t2.pt(), t2.eta(), t2.phi(), leptonM2, t2.sign(), t2.dcaXY(), t2.dcaZ(), t2.cYY(), t2.cZY(), t2.cZZ(), + // t2.x(), t2.y(), t2.z(), t2.alpha(), t2.snp(), t2.tgl(), + // t2.cSnpY(), t2.cSnpZ(), t2.cSnpSnp(), + // t2.cTglY(), t2.cTglZ(), t2.cTglSnp(), t2.cTglTgl(), + // t2.c1PtY(), t2.c1PtZ(), t2.c1PtSnp(), t2.c1PtTgl(), t2.c1Pt21Pt2())); + } else { + emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrack(candidate.pt2, candidate.eta2, candidate.phi2, leptonM2, t2.sign(), t2.dcaXY(), t2.dcaZ(), t2.cYY(), t2.cZY(), t2.cZZ())); + // emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrackWithCov(t2.pt(), t2.eta(), t2.phi(), leptonM2, t2.sign(), t2.dcaXY(), t2.dcaZ(), t2.cYY(), t2.cZY(), t2.cZZ(), + // t2.x(), t2.y(), t2.z(), t2.alpha(), t2.snp(), t2.tgl(), + // t2.cSnpY(), t2.cSnpZ(), t2.cSnpSnp(), + // t2.cTglY(), t2.cTglZ(), t2.cTglSnp(), t2.cTglTgl(), + // t2.c1PtY(), t2.c1PtZ(), t2.c1PtSnp(), t2.c1PtTgl(), t2.c1Pt21Pt2())); + } + } + } + } else if (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + if (std::find(used_trackIds_per_col.begin(), used_trackIds_per_col.end(), t1.globalIndex()) == used_trackIds_per_col.end()) { + used_trackIds_per_col.emplace_back(t1.globalIndex()); + if (cfgDoMix) { + if (t1.sign() > 0) { + emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrack(candidate.pt1, candidate.eta1, candidate.phi1, leptonM1, t1.sign(), t1.fwdDcaX(), t1.fwdDcaY(), t1.cXX(), t1.cXY(), t1.cYY())); + // emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrackWithCov(t1.pt(), t1.eta(), t1.phi(), leptonM1, t1.sign(), t1.fwdDcaX(), t1.fwdDcaY(), t1.cXX(), t1.cXY(), t1.cYY(), + // t1.x(), t1.y(), t1.z(), t1.tgl(), + // t1.cPhiX(), t1.cPhiY(), t1.cPhiPhi(), + // t1.cTglX(), t1.cTglY(), t1.cTglPhi(), t1.cTglTgl(), + // t1.c1PtX(), t1.c1PtY(), t1.c1PtPhi(), t1.c1PtTgl(), t1.c1Pt21Pt2(), t1.chi2())); + + } else { + emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrack(candidate.pt1, candidate.eta1, candidate.phi1, leptonM1, t1.sign(), t1.fwdDcaX(), t1.fwdDcaY(), t1.cXX(), t1.cXY(), t1.cYY())); + // emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrackWithCov(t1.pt(), t1.eta(), t1.phi(), leptonM1, t1.sign(), t1.fwdDcaX(), t1.fwdDcaY(), t1.cXX(), t1.cXY(), t1.cYY(), + // t1.x(), t1.y(), t1.z(), t1.tgl(), + // t1.cPhiX(), t1.cPhiY(), t1.cPhiPhi(), + // t1.cTglX(), t1.cTglY(), t1.cTglPhi(), t1.cTglTgl(), + // t1.c1PtX(), t1.c1PtY(), t1.c1PtPhi(), t1.c1PtTgl(), t1.c1Pt21Pt2(), t1.chi2())); + } + } + } + if (std::find(used_trackIds_per_col.begin(), used_trackIds_per_col.end(), t2.globalIndex()) == used_trackIds_per_col.end()) { + used_trackIds_per_col.emplace_back(t2.globalIndex()); + if (cfgDoMix) { + if (t2.sign() > 0) { + emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrack(candidate.pt2, candidate.eta2, candidate.phi2, leptonM2, t2.sign(), t2.fwdDcaX(), t2.fwdDcaY(), t2.cXX(), t2.cXY(), t2.cYY())); + // emh_pos->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrackWithCov(t2.pt(), t2.eta(), t2.phi(), leptonM2, t2.sign(), t2.fwdDcaX(), t2.fwdDcaY(), t2.cXX(), t2.cXY(), t2.cYY(), + // t2.x(), t2.y(), t2.z(), t2.tgl(), + // t2.cPhiX(), t2.cPhiY(), t2.cPhiPhi(), + // t2.cTglX(), t2.cTglY(), t2.cTglPhi(), t2.cTglTgl(), + // t2.c1PtX(), t2.c1PtY(), t2.c1PtPhi(), t2.c1PtTgl(), t2.c1Pt21Pt2(), t2.chi2())); + } else { + emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrack(candidate.pt2, candidate.eta2, candidate.phi2, leptonM2, t2.sign(), t2.fwdDcaX(), t2.fwdDcaY(), t2.cXX(), t2.cXY(), t2.cYY())); + // emh_neg->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMFwdTrackWithCov(t2.pt(), t2.eta(), t2.phi(), leptonM2, t2.sign(), t2.fwdDcaX(), t2.fwdDcaY(), t2.cXX(), t2.cXY(), t2.cYY(), + // t2.x(), t2.y(), t2.z(), t2.tgl(), + // t2.cPhiX(), t2.cPhiY(), t2.cPhiPhi(), + // t2.cTglX(), t2.cTglY(), t2.cTglPhi(), t2.cTglTgl(), + // t2.c1PtX(), t2.c1PtY(), t2.c1PtPhi(), t2.c1PtTgl(), t2.c1Pt21Pt2(), t2.chi2())); + } + } + } + } + } + return true; + } + + o2::framework::expressions::Filter collisionFilter_centrality = (eventcuts.cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < eventcuts.cfgCentMax) || (eventcuts.cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < eventcuts.cfgCentMax) || (eventcuts.cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < eventcuts.cfgCentMax); + o2::framework::expressions::Filter collisionFilter_numContrib = eventcuts.cfgNumContribMin <= o2::aod::collision::numContrib && o2::aod::collision::numContrib < eventcuts.cfgNumContribMax; + o2::framework::expressions::Filter collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; + o2::framework::expressions::Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; + using FilteredMyCollisions = o2::soa::Filtered; + + o2::framework::SliceCache cache; + o2::framework::Preslice perCollision_electron = o2::aod::emprimaryelectron::emeventId; + o2::framework::expressions::Filter trackFilter_electron = dielectroncuts.cfg_min_pt_track < o2::aod::track::pt && dielectroncuts.cfg_min_eta_track < o2::aod::track::eta && o2::aod::track::eta < dielectroncuts.cfg_max_eta_track && nabs(o2::aod::track::dcaXY) < dielectroncuts.cfg_max_dcaxy && nabs(o2::aod::track::dcaZ) < dielectroncuts.cfg_max_dcaz && o2::aod::track::itsChi2NCl < dielectroncuts.cfg_max_chi2its && o2::aod::track::tpcChi2NCl < dielectroncuts.cfg_max_chi2tpc; + o2::framework::expressions::Filter pidFilter_electron = dielectroncuts.cfg_min_TPCNsigmaEl < o2::aod::pidtpc::tpcNSigmaEl && o2::aod::pidtpc::tpcNSigmaEl < dielectroncuts.cfg_max_TPCNsigmaEl; + o2::framework::expressions::Filter ttcaFilter_electron = ifnode(dielectroncuts.enableTTCA.node(), true, o2::aod::emprimaryelectron::isAssociatedToMPC == true); + o2::framework::expressions::Filter prefilter_derived_electron = ifnode(dielectroncuts.cfg_apply_cuts_from_prefilter_derived.node() && dielectroncuts.cfg_prefilter_bits_derived.node() >= static_cast(1), + ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kMee))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kMee))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kPhiV))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kPhiV))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kPhiVLS))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kPhiVLS))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULSAtRefR))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULSAtRefR))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLSAtRefR))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLSAtRefR))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) <= static_cast(0), true), + o2::aod::emprimaryelectron::pfbderived >= static_cast(0)); + + o2::framework::expressions::Filter prefilter_electron = ifnode(dielectroncuts.cfg_apply_cuts_from_prefilter.node() && dielectroncuts.cfg_prefilter_bits.node() >= static_cast(1), + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPC))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPC))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_20MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_20MeV))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_40MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_40MeV))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_60MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_60MeV))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_80MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_80MeV))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_100MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_100MeV))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_120MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_120MeV))) <= static_cast(0), true) && + ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_140MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_140MeV))) <= static_cast(0), true), + o2::aod::emprimaryelectron::pfb >= static_cast(0)); + + using TElectronType = std::tuple_element_t<0, std::tuple>; + o2::framework::Partition positive_electrons = o2::aod::emprimaryelectron::sign > int8_t(0); + o2::framework::Partition negative_electrons = o2::aod::emprimaryelectron::sign < int8_t(0); + + o2::framework::Preslice perCollision_muon = o2::aod::emprimarymuon::emeventId; + o2::framework::expressions::Filter trackFilter_muon = o2::aod::fwdtrack::trackType == dimuoncuts.cfg_track_type && dimuoncuts.cfg_min_pt_track < o2::aod::fwdtrack::pt && o2::aod::fwdtrack::pt < dimuoncuts.cfg_max_pt_track && dimuoncuts.cfg_min_eta_track < o2::aod::fwdtrack::eta && o2::aod::fwdtrack::eta < dimuoncuts.cfg_max_eta_track; + o2::framework::expressions::Filter ttcaFilter_muon = ifnode(dimuoncuts.enableTTCA.node(), true, o2::aod::emprimarymuon::isAssociatedToMPC == true); + o2::framework::expressions::Filter prefilter_derived_muon = ifnode(dimuoncuts.cfg_apply_cuts_from_prefilter_derived.node() && dimuoncuts.cfg_prefilter_bits_derived.node() >= static_cast(1), + ifnode((dimuoncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) > static_cast(0), (o2::aod::emprimarymuon::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) <= static_cast(0), true) && + ifnode((dimuoncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) > static_cast(0), (o2::aod::emprimarymuon::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) <= static_cast(0), true), + o2::aod::emprimarymuon::pfbderived >= static_cast(0)); + + o2::framework::Partition positive_muons = o2::aod::emprimarymuon::sign > int8_t(0); + o2::framework::Partition negative_muons = o2::aod::emprimarymuon::sign < int8_t(0); + + TEMH* emh_pos = nullptr; + TEMH* emh_neg = nullptr; + + std::map, uint64_t> map_mixed_eventId_to_globalBC; + std::unordered_map map_best_match_globalmuon; + + std::vector used_trackIds_per_col; + int ndf = 0; + + template + void runPairing(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut, TAllTracks const& tracks) + { + for (const auto& collision : collisions) { + initCCDB(collision); + + const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + float centrality = centralities[cfgCentEstimator]; + if (centralities[cfgCentEstimator] < eventcuts.cfgCentMin || eventcuts.cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + std::array q2ft0m = {collision.q2xft0m(), collision.q2yft0m()}; + std::array q2ft0a = {collision.q2xft0a(), collision.q2yft0a()}; + std::array q2ft0c = {collision.q2xft0c(), collision.q2yft0c()}; + std::array q2btot = {collision.q2xbtot(), collision.q2ybtot()}; + std::array q2bpos = {collision.q2xbpos(), collision.q2ybpos()}; + std::array q2bneg = {collision.q2xbneg(), collision.q2ybneg()}; + std::array q2fv0a = {collision.q2xfv0a(), collision.q2yfv0a()}; + const float eventplanes_2_for_mix[7] = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg(), collision.ep2fv0a()}; + float ep2 = eventplanes_2_for_mix[cfgEP2Estimator_for_Mix]; + + std::vector>> qvectors = { + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 0th harmonics + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 1st harmonics + {q2ft0m, q2ft0a, q2ft0c, q2btot, q2bpos, q2bneg, q2fv0a}, // 2nd harmonics + }; + + if (nmod == 2) { + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0, 2>(&fRegistry, collision); + } else { + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0, -1>(&fRegistry, collision); + } + if (nmod < 0 || isGoodQvector(qvectors)) { + fRegistry.fill(HIST("Event/before/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev - 1); // is qvector calibarated + } + fRegistry.fill(HIST("Event/before/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); + + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + if (eventcuts.cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { + continue; + } + + if (nmod > 0 && !isGoodQvector(qvectors)) { + continue; + } + + if constexpr (isTriggerAnalysis) { + if (!zorro.isSelected(collision.globalBC(), zorroGroup.bcMarginForSoftwareTrigger)) { // triggered event + continue; + } + + auto swt_bitset = zorro.getLastResult(); // this has to be called after zorro::isSelected + auto TOIcounter = zorro.getTOIcounters()[0]; // this has to be called after zorro::isSelected + auto ATcounter = zorro.getATcounters()[mToIidx]; // this has to be called after zorro::isSelected + + if (swt_bitset.test(mToIidx)) { + while (ATcounter > mATCounter) { + mATCounter++; + fRegistry.fill(HIST("Event/trigger/hAnalysedTrigger"), collision.runNumber()); + } + + while (TOIcounter > mTOICounter) { + fRegistry.fill(HIST("Event/trigger/hAnalysedToI"), collision.runNumber()); + mTOICounter++; // always incremented by 1 in zorro!! + } + // LOGF(info, "collision.globalIndex() = %d, collision.globalBC() = %llu, mTOICounter = %d, mATcounter = %d", collision.globalIndex(), collision.globalBC(), mTOICounter, mATCounter); + } + } + + std::vector bootstrapweights = {}; + if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) { // bootstrapping for accepted events + int randomSeed = static_cast(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + TRandom3 randomNumber(randomSeed); + for (int i = 0; i < cfgNumBootstrapSamples; i++) { + float poissonweight = 0.; + poissonweight = static_cast(randomNumber.PoissonD(1.0)); + bootstrapweights.push_back(poissonweight); + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfoBootstrap(&fRegistry, collision, i, poissonweight); + } + } else { + bootstrapweights.push_back(1.0); // to pass as non-empyt dummy to use + } + + if (nmod == 2) { + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, 2>(&fRegistry, collision); + } else { + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, -1>(&fRegistry, collision); + } + fRegistry.fill(HIST("Event/after/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev - 1); // is qvector calibarated + + fRegistry.fill(HIST("Event/before/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + fRegistry.fill(HIST("Event/after/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); + + auto posTracks_per_coll = posTracks.sliceByCached(perCollision, collision.globalIndex(), cache); + auto negTracks_per_coll = negTracks.sliceByCached(perCollision, collision.globalIndex(), cache); + // LOGF(info, "collision.globalIndex() = %d , collision.posZ() = %f , collision.numContrib() = %d, centrality = %f , posTracks_per_coll.size() = %d, negTracks_per_coll.size() = %d", collision.globalIndex(), collision.posZ(), collision.numContrib(), centralities[cfgCentEstimator], posTracks_per_coll.size(), negTracks_per_coll.size()); + + used_trackIds_per_col.reserve(posTracks_per_coll.size() + negTracks_per_coll.size()); + int nuls = 0, nlspp = 0, nlsmm = 0; + for (const auto& [pos, neg] : combinations(o2::soa::CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS + bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks, bootstrapweights); + if (is_pair_ok) { + nuls++; + } + } + for (const auto& [pos1, pos2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ + bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks, bootstrapweights); + if (is_pair_ok) { + nlspp++; + } + } + for (const auto& [neg1, neg2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- + bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks, bootstrapweights); + if (is_pair_ok) { + nlsmm++; + } + } + used_trackIds_per_col.clear(); + used_trackIds_per_col.shrink_to_fit(); + + if (!cfgDoMix || !(nuls > 0 || nlspp > 0 || nlsmm > 0)) { + continue; + } + + // event mixing + int zbin = lower_bound(zvtx_bin_edges.begin(), zvtx_bin_edges.end(), collision.posZ()) - zvtx_bin_edges.begin() - 1; + if (zbin < 0) { + zbin = 0; + } else if (static_cast(zvtx_bin_edges.size()) - 2 < zbin) { + zbin = static_cast(zvtx_bin_edges.size()) - 2; + } + + int centbin = lower_bound(cent_bin_edges.begin(), cent_bin_edges.end(), centrality) - cent_bin_edges.begin() - 1; + if (centbin < 0) { + centbin = 0; + } else if (static_cast(cent_bin_edges.size()) - 2 < centbin) { + centbin = static_cast(cent_bin_edges.size()) - 2; + } + + int epbin = lower_bound(ep_bin_edges.begin(), ep_bin_edges.end(), ep2) - ep_bin_edges.begin() - 1; + if (epbin < 0) { + epbin = 0; + } else if (static_cast(ep_bin_edges.size()) - 2 < epbin) { + epbin = static_cast(ep_bin_edges.size()) - 2; + } + + int occbin = -1; + if (cfgOccupancyEstimator == 0) { + occbin = lower_bound(occ_bin_edges.begin(), occ_bin_edges.end(), collision.ft0cOccupancyInTimeRange()) - occ_bin_edges.begin() - 1; + } else if (cfgOccupancyEstimator == 1) { + occbin = lower_bound(occ_bin_edges.begin(), occ_bin_edges.end(), collision.trackOccupancyInTimeRange()) - occ_bin_edges.begin() - 1; + } else { + occbin = lower_bound(occ_bin_edges.begin(), occ_bin_edges.end(), collision.ft0cOccupancyInTimeRange()) - occ_bin_edges.begin() - 1; + } + + if (occbin < 0) { + occbin = 0; + } else if (static_cast(occ_bin_edges.size()) - 2 < occbin) { + occbin = static_cast(occ_bin_edges.size()) - 2; + } + + // LOGF(info, "collision.globalIndex() = %d, collision.posZ() = %f, centrality = %f, ep2 = %f, collision.trackOccupancyInTimeRange() = %d, zbin = %d, centbin = %d, epbin = %d, occbin = %d", collision.globalIndex(), collision.posZ(), centrality, ep2, collision.trackOccupancyInTimeRange(), zbin, centbin, epbin, occbin); + + std::tuple key_bin = std::make_tuple(zbin, centbin, epbin, occbin); + std::pair key_df_collision = std::make_pair(ndf, collision.globalIndex()); // this gives the current event. + + // make a vector of selected photons in this collision. + auto selected_posTracks_in_this_event = emh_pos->GetTracksPerCollision(key_df_collision); + auto selected_negTracks_in_this_event = emh_neg->GetTracksPerCollision(key_df_collision); + // LOGF(info, "N selected tracks in current event (%d, %d), zvtx = %f, centrality = %f , npos = %d , nneg = %d, nuls = %d , nlspp = %d, nlsmm = %d", ndf, collision.globalIndex(), collision.posZ(), centralities[cfgCentEstimator], selected_posTracks_in_this_event.size(), selected_negTracks_in_this_event.size(), nuls, nlspp, nlsmm); + + auto collisionIds_in_mixing_pool = emh_pos->GetCollisionIdsFromEventPool(key_bin); // pos/neg does not matter. + // LOGF(info, "collisionIds_in_mixing_pool.size() = %d", collisionIds_in_mixing_pool.size()); + + for (const auto& mix_dfId_collisionId : collisionIds_in_mixing_pool) { + int mix_dfId = mix_dfId_collisionId.first; + int mix_collisionId = mix_dfId_collisionId.second; + if (collision.globalIndex() == mix_collisionId && ndf == mix_dfId) { // this never happens. only protection. + continue; + } + + auto globalBC_mix = map_mixed_eventId_to_globalBC[mix_dfId_collisionId]; + uint64_t diffBC = std::max(collision.globalBC(), globalBC_mix) - std::min(collision.globalBC(), globalBC_mix); + fRegistry.fill(HIST("Pair/mix/hDiffBC"), diffBC); + if (diffBC < ndiff_bc_mix) { + continue; + } + + auto posTracks_from_event_pool = emh_pos->GetTracksPerCollision(mix_dfId_collisionId); + auto negTracks_from_event_pool = emh_neg->GetTracksPerCollision(mix_dfId_collisionId); + // LOGF(info, "Do event mixing: current event (%d, %d) | event pool (%d, %d), npos = %d , nneg = %d", ndf, collision.globalIndex(), mix_dfId, mix_collisionId, posTracks_from_event_pool.size(), negTracks_from_event_pool.size()); + + for (const auto& pos : selected_posTracks_in_this_event) { // ULS mix + for (const auto& neg : negTracks_from_event_pool) { + fillPairInfo<1>(collision, pos, neg, cut, nullptr, bootstrapweights); + } + } + + for (const auto& neg : selected_negTracks_in_this_event) { // ULS mix + for (const auto& pos : posTracks_from_event_pool) { + fillPairInfo<1>(collision, neg, pos, cut, nullptr, bootstrapweights); + } + } + + for (const auto& pos1 : selected_posTracks_in_this_event) { // LS++ mix + for (const auto& pos2 : posTracks_from_event_pool) { + fillPairInfo<1>(collision, pos1, pos2, cut, nullptr, bootstrapweights); + } + } + + for (const auto& neg1 : selected_negTracks_in_this_event) { // LS-- mix + for (const auto& neg2 : negTracks_from_event_pool) { + fillPairInfo<1>(collision, neg1, neg2, cut, nullptr, bootstrapweights); + } + } + } // end of loop over mixed event pool + + if (nuls > 0 || nlspp > 0 || nlsmm > 0) { + map_mixed_eventId_to_globalBC[key_df_collision] = collision.globalBC(); + emh_pos->AddCollisionIdAtLast(key_bin, key_df_collision); + emh_neg->AddCollisionIdAtLast(key_bin, key_df_collision); + } // end of if pair exist + + } // end of collision loop + + } // end of DF + + template + bool isPairOK(TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&) + { + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { + if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { + return false; + } + } else { // cut-based + if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { + return false; + } + } + + if constexpr (withSCT) { + if (foundHFSV(t1) || foundHFSV(t2)) { + return false; + } + } + + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + if (!cut.IsSelectedTrack(t1) || !cut.IsSelectedTrack(t2)) { + return false; + } + if (!map_best_match_globalmuon[t1.globalIndex()] || !map_best_match_globalmuon[t2.globalIndex()]) { + return false; + } + + // if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + // return false; + // } + // if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + // return false; + // } + } + + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + if (!cut.IsSelectedPair(t1, t2)) { + return false; + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + if (!cut.IsSelectedPair(t1, t2)) { + return false; + } + } + return true; + } + + std::map, float> map_weight; // -> float + template + void fillPairWeightMap(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut, TAllTracks const& tracks) + { + std::vector> passed_pairIds; + passed_pairIds.reserve(posTracks.size() * negTracks.size()); + + for (const auto& collision : collisions) { + initCCDB(collision); + const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + if (centralities[cfgCentEstimator] < eventcuts.cfgCentMin || eventcuts.cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + // if constexpr (isTriggerAnalysis) { // don't call this here. this causes double counting. + // if (!zorro.isSelected(collision.globalBC(), zorroGroup.bcMarginForSoftwareTrigger)) { // triggered event + // continue; + // } + // } + + std::array q2ft0m = {collision.q2xft0m(), collision.q2yft0m()}; + std::array q2ft0a = {collision.q2xft0a(), collision.q2yft0a()}; + std::array q2ft0c = {collision.q2xft0c(), collision.q2yft0c()}; + std::array q2btot = {collision.q2xbtot(), collision.q2ybtot()}; + std::array q2bpos = {collision.q2xbpos(), collision.q2ybpos()}; + std::array q2bneg = {collision.q2xbneg(), collision.q2ybneg()}; + std::array q2fv0a = {collision.q2xfv0a(), collision.q2yfv0a()}; + + std::vector>> qvectors = { + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 0th harmonics + {{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 1st harmonics + {q2ft0m, q2ft0a, q2ft0c, q2btot, q2bpos, q2bneg, q2fv0a}, // 2nd harmonics + }; + + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + if (eventcuts.cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { + continue; + } + + if (nmod > 0 && !isGoodQvector(qvectors)) { + continue; + } + + auto posTracks_per_coll = posTracks.sliceByCached(perCollision, collision.globalIndex(), cache); + auto negTracks_per_coll = negTracks.sliceByCached(perCollision, collision.globalIndex(), cache); + + for (const auto& [pos, neg] : combinations(o2::soa::CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS + if (isPairOK(pos, neg, cut, tracks)) { + passed_pairIds.emplace_back(std::make_pair(pos.globalIndex(), neg.globalIndex())); + } + } + for (const auto& [pos1, pos2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ + if (isPairOK(pos1, pos2, cut, tracks)) { + passed_pairIds.emplace_back(std::make_pair(pos1.globalIndex(), pos2.globalIndex())); + } + } + for (const auto& [neg1, neg2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- + if (isPairOK(neg1, neg2, cut, tracks)) { + passed_pairIds.emplace_back(std::make_pair(neg1.globalIndex(), neg2.globalIndex())); + } + } + } // end of collision loop + + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + for (const auto& pairId : passed_pairIds) { + auto t1 = tracks.rawIteratorAt(std::get<0>(pairId)); + auto t2 = tracks.rawIteratorAt(std::get<1>(pairId)); + + float n = 1.f; // include myself. + for (const auto& ambId1 : t1.ambiguousElectronsIds()) { + for (const auto& ambId2 : t2.ambiguousElectronsIds()) { + if (std::find(passed_pairIds.begin(), passed_pairIds.end(), std::make_pair(ambId1, ambId2)) != passed_pairIds.end()) { + n += 1.f; + } + } + } + map_weight[pairId] = 1.f / n; + } // end of passed_pairIds loop + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + for (const auto& pairId : passed_pairIds) { + auto t1 = tracks.rawIteratorAt(std::get<0>(pairId)); + auto t2 = tracks.rawIteratorAt(std::get<1>(pairId)); + + float n = 1.f; // include myself. + for (const auto& ambId1 : t1.ambiguousMuonsIds()) { + for (const auto& ambId2 : t2.ambiguousMuonsIds()) { + if (std::find(passed_pairIds.begin(), passed_pairIds.end(), std::make_pair(ambId1, ambId2)) != passed_pairIds.end()) { + n += 1.f; + } + } + } + map_weight[pairId] = 1.f / n; + } // end of passed_pairIds loop + } + passed_pairIds.clear(); + passed_pairIds.shrink_to_fit(); + } + + void processAnalysis(FilteredMyCollisions const& collisions, Types const&... args) + { + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + auto electrons = std::get<0>(std::tie(args...)); + if (cfgApplyWeightTTCA) { + fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); + } + runPairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + auto muons = std::get<0>(std::tie(args...)); + map_best_match_globalmuon = o2::aod::pwgem::dilepton::utils::emtrackutil::findBestMatchMap(muons, fDimuonCut); + if (cfgApplyWeightTTCA) { + fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); + } + runPairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); + } + map_weight.clear(); + map_best_match_globalmuon.clear(); + ndf++; + } + PROCESS_SWITCH(DileptonSV, processAnalysis, "run dilepton analysis", true); + + void processTriggerAnalysis(FilteredMyCollisions const& collisions, Types const&... args) + { + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + auto electrons = std::get<0>(std::tie(args...)); + if (cfgApplyWeightTTCA) { + fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); + } + runPairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + auto muons = std::get<0>(std::tie(args...)); + map_best_match_globalmuon = o2::aod::pwgem::dilepton::utils::emtrackutil::findBestMatchMap(muons, fDimuonCut); + if (cfgApplyWeightTTCA) { + fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); + } + runPairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); + } + map_weight.clear(); + map_best_match_globalmuon.clear(); + ndf++; + } + PROCESS_SWITCH(DileptonSV, processTriggerAnalysis, "run dilepton analysis on triggered data", false); + + void processNorm(o2::aod::EMEventNormInfos const& collisions) + { + for (const auto& collision : collisions) { + if (collision.centFT0C() < eventcuts.cfgCentMin || eventcuts.cfgCentMax < collision.centFT0C()) { + continue; + } + fRegistry.fill(HIST("Event/norm/hZvtx"), collision.posZ()); + + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 1.0); + if (collision.selection_bit(o2::aod::emevsel::kIsTriggerTVX)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 2.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoTimeFrameBorder)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 3.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoITSROFrameBorder)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 4.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoSameBunchPileup)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 5.0); + } + if (collision.selection_bit(o2::aod::emevsel::kIsGoodZvtxFT0vsPV)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 6.0); + } + if (collision.selection_bit(o2::aod::emevsel::kIsVertexITSTPC)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 7.0); + } + if (collision.selection_bit(o2::aod::emevsel::kIsVertexTRDmatched)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 8.0); + } + if (collision.selection_bit(o2::aod::emevsel::kIsVertexTOFmatched)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 9.0); + } + if (collision.sel8()) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 10.0); + } + if (std::fabs(collision.posZ()) < 10.0) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 11.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoCollInTimeRangeStandard)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 12.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoCollInTimeRangeStrict)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 13.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoCollInRofStandard)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 14.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoCollInRofStrict)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 15.0); + } + if (collision.selection_bit(o2::aod::emevsel::kNoHighMultCollInPrevRof)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 16.0); + } + if (collision.selection_bit(o2::aod::emevsel::kIsGoodITSLayer3)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 17.0); + } + if (collision.selection_bit(o2::aod::emevsel::kIsGoodITSLayer0123)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 18.0); + } + if (collision.selection_bit(o2::aod::emevsel::kIsGoodITSLayersAll)) { + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), 19.0); + } + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + if (eventcuts.cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { + continue; + } + fRegistry.fill(HIST("Event/norm/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + } // end of collision loop + } + PROCESS_SWITCH(DileptonSV, processNorm, "process normalization info", false); + + void processBC(o2::aod::EMBCs const& bcs) + { + for (const auto& bc : bcs) { + if (bc.selection_bit(o2::aod::emevsel::kIsTriggerTVX)) { + fRegistry.fill(HIST("BC/hTVXCounter"), 0.f); + + if (bc.selection_bit(o2::aod::emevsel::kNoTimeFrameBorder)) { + fRegistry.fill(HIST("BC/hTVXCounter"), 1.f); + } + if (bc.selection_bit(o2::aod::emevsel::kNoITSROFrameBorder)) { + fRegistry.fill(HIST("BC/hTVXCounter"), 2.f); + } + if (rctChecker(bc)) { + fRegistry.fill(HIST("BC/hTVXCounter"), 3.f); + } + if (bc.selection_bit(o2::aod::emevsel::kNoTimeFrameBorder) && bc.selection_bit(o2::aod::emevsel::kNoITSROFrameBorder)) { + fRegistry.fill(HIST("BC/hTVXCounter"), 4.f); + } + if (bc.selection_bit(o2::aod::emevsel::kNoTimeFrameBorder) && bc.selection_bit(o2::aod::emevsel::kNoITSROFrameBorder) && rctChecker(bc)) { + fRegistry.fill(HIST("BC/hTVXCounter"), 5.f); + } + } + } + } + PROCESS_SWITCH(DileptonSV, processBC, "process BC counter", false); + + void processDummy(MyCollisions const&) {} + PROCESS_SWITCH(DileptonSV, processDummy, "Dummy function", false); +}; + +#endif // PWGEM_DILEPTON_CORE_DILEPTONSV_H_ diff --git a/PWGEM/Dilepton/TableProducer/CMakeLists.txt b/PWGEM/Dilepton/TableProducer/CMakeLists.txt index 0070a556748..2a2cbf861c5 100644 --- a/PWGEM/Dilepton/TableProducer/CMakeLists.txt +++ b/PWGEM/Dilepton/TableProducer/CMakeLists.txt @@ -31,9 +31,9 @@ o2physics_add_dpl_workflow(skimmer-primary-electron PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::MLCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(skimmer-primary-electron-sct - SOURCES skimmerPrimaryElectronSCT.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::MLCore +o2physics_add_dpl_workflow(skimmer-primary-electron-sv + SOURCES skimmerPrimaryElectronSV.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore O2Physics::MLCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(skimmer-primary-electron-qc diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectronSCT.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectronSV.cxx similarity index 93% rename from PWGEM/Dilepton/TableProducer/skimmerPrimaryElectronSCT.cxx rename to PWGEM/Dilepton/TableProducer/skimmerPrimaryElectronSV.cxx index 2fa33093663..3402faa2e42 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectronSCT.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectronSV.cxx @@ -23,6 +23,8 @@ #include #include +#include +#include #include #include #include @@ -46,7 +48,7 @@ #include -struct skimmerPrimaryElectronSCT { +struct skimmerPrimaryElectronSV { // using MyBCs = o2::soa::Join; using MyBCs = o2::soa::Join; @@ -188,7 +190,7 @@ struct skimmerPrimaryElectronSCT { initCCDB(bc); electronModule.processWithoutTTCA(bcs, collisions, tracks, v0s, cascades, nullptr, products, mRegistry); } - PROCESS_SWITCH(skimmerPrimaryElectronSCT, processRec_SA, "process reconstructed info only", true); // standalone + PROCESS_SWITCH(skimmerPrimaryElectronSV, processRec_SA, "process reconstructed info only", true); // standalone void processRec_TTCA(MyCollisions const& collisions, MyBCs const& bcs, MyTracks const& tracks, o2::aod::TrackAssoc const& trackIndices, filteredMyV0s const& v0s, filteredMyCascades const& cascades) { @@ -199,7 +201,7 @@ struct skimmerPrimaryElectronSCT { initCCDB(bc); electronModule.processWithTTCA(bcs, collisions, tracks, v0s, cascades, trackIndices, nullptr, products, mRegistry, cache, perCol_track, trackIndicesPerCollision, perCol_v0, perCol_casc); } - PROCESS_SWITCH(skimmerPrimaryElectronSCT, processRec_TTCA, "process reconstructed info only", false); // with TTCA + PROCESS_SWITCH(skimmerPrimaryElectronSV, processRec_TTCA, "process reconstructed info only", false); // with TTCA void processRec_SA_SWT(MyCollisionsWithSWT const& collisions, MyBCs const& bcs, MyTracks const& tracks, filteredMyV0s const& v0s, filteredMyCascades const& cascades) { @@ -210,7 +212,7 @@ struct skimmerPrimaryElectronSCT { initCCDB(bc); electronModule.processWithoutTTCA(bcs, collisions, tracks, v0s, cascades, nullptr, products, mRegistry); } - PROCESS_SWITCH(skimmerPrimaryElectronSCT, processRec_SA_SWT, "process reconstructed info only", false); // standalone with swt + PROCESS_SWITCH(skimmerPrimaryElectronSV, processRec_SA_SWT, "process reconstructed info only", false); // standalone with swt void processRec_TTCA_SWT(MyCollisionsWithSWT const& collisions, MyBCs const& bcs, MyTracks const& tracks, o2::aod::TrackAssoc const& trackIndices, filteredMyV0s const& v0s, filteredMyCascades const& cascades) { @@ -221,7 +223,7 @@ struct skimmerPrimaryElectronSCT { initCCDB(bc); electronModule.processWithTTCA(bcs, collisions, tracks, v0s, cascades, trackIndices, nullptr, products, mRegistry, cache, perCol_track, trackIndicesPerCollision, perCol_v0, perCol_casc); } - PROCESS_SWITCH(skimmerPrimaryElectronSCT, processRec_TTCA_SWT, "process reconstructed info only", false); // with TTCA with swt + PROCESS_SWITCH(skimmerPrimaryElectronSV, processRec_TTCA_SWT, "process reconstructed info only", false); // with TTCA with swt // ---------- for MC ---------- @@ -234,7 +236,7 @@ struct skimmerPrimaryElectronSCT { initCCDB(bc); electronModule.processWithoutTTCA(bcs, collisions, tracks, v0s, cascades, mcParticles, products, mRegistry); } - PROCESS_SWITCH(skimmerPrimaryElectronSCT, processMC_SA, "process reconstructed and MC info ", false); // without TTCA in MC + PROCESS_SWITCH(skimmerPrimaryElectronSV, processMC_SA, "process reconstructed and MC info ", false); // without TTCA in MC void processMC_TTCA(o2::soa::Join const& collisions, MyBCs const& bcs, MyTracksMC const& tracks, o2::aod::TrackAssoc const& trackIndices, filteredMyV0s const& v0s, filteredMyCascades const& cascades, o2::aod::McParticles const& mcParticles) { @@ -245,7 +247,7 @@ struct skimmerPrimaryElectronSCT { initCCDB(bc); electronModule.processWithTTCA(bcs, collisions, tracks, v0s, cascades, trackIndices, mcParticles, products, mRegistry, cache, perCol_track, trackIndicesPerCollision, perCol_v0, perCol_casc); } - PROCESS_SWITCH(skimmerPrimaryElectronSCT, processMC_TTCA, "process reconstructed info only", false); // with TTCA in MC + PROCESS_SWITCH(skimmerPrimaryElectronSV, processMC_TTCA, "process reconstructed info only", false); // with TTCA in MC }; struct associateAmbiguousElectron { @@ -276,6 +278,6 @@ o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext co { o2::pid::tof::TOFResponseImpl::metadataInfo.initMetadata(cfgc); return o2::framework::WorkflowSpec{ - adaptAnalysisTask(cfgc, o2::framework::TaskName{"skimmer-primary-electron-sct"}), + adaptAnalysisTask(cfgc, o2::framework::TaskName{"skimmer-primary-electron-sv"}), adaptAnalysisTask(cfgc, o2::framework::TaskName{"associate-ambiguous-electron"})}; } diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt b/PWGEM/Dilepton/Tasks/CMakeLists.txt index c872babef87..5a1547c8a65 100644 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -90,6 +90,11 @@ o2physics_add_dpl_workflow(dielectron-mc PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMDileptonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(dielectron-sv + SOURCES dielectronSV.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::DCAFitter O2Physics::AnalysisCore O2Physics::PWGEMDileptonCore O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(dielectron-sct SOURCES dielectronSCT.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMDileptonCore O2Physics::EventFilteringUtils diff --git a/PWGEM/Dilepton/Tasks/dielectronSV.cxx b/PWGEM/Dilepton/Tasks/dielectronSV.cxx new file mode 100644 index 00000000000..d61e9beea71 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/dielectronSV.cxx @@ -0,0 +1,29 @@ +// 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. +// +// ======================== +// +// This code is for dielectron analyses. +// Please write to: daiki.sekihata@cern.ch + +#include "PWGEM/Dilepton/Core/DileptonSV.h" +#include "PWGEM/Dilepton/Utils/PairUtilities.h" + +#include +#include + +using namespace o2::framework; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask>(cfgc, TaskName{"dielectron-sv"})}; +}