diff --git a/PWGEM/Dilepton/Core/DileptonProducer.h b/PWGEM/Dilepton/Core/DileptonProducer.h new file mode 100644 index 00000000000..057aedfc3bc --- /dev/null +++ b/PWGEM/Dilepton/Core/DileptonProducer.h @@ -0,0 +1,777 @@ +// 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_DILEPTONPRODUCER_H_ +#define PWGEM_DILEPTON_CORE_DILEPTONPRODUCER_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/MlResponseDielectronSingleTrack.h" +#include "PWGEM/Dilepton/Utils/PairUtilities.h" + +#include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/trackUtilities.h" +#include "Tools/ML/MlResponse.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/LHCConstants.h" +#include "DataFormatsParameters/GRPECSObject.h" +#include "DataFormatsParameters/GRPLHCIFData.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "MathUtils/Utils.h" + +#include "Math/Vector4D.h" +#include "TH1D.h" +#include "TString.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::dilepton::utils; +using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; +using namespace o2::aod::pwgem::dilepton::utils::pairutil; + +using MyCollisions = soa::Join; +using MyCollision = MyCollisions::iterator; + +using MyElectrons = soa::Join; +using MyElectron = MyElectrons::iterator; +using FilteredMyElectrons = soa::Filtered; +using FilteredMyElectron = FilteredMyElectrons::iterator; + +using MyMuons = soa::Join; +using MyMuon = MyMuons::iterator; +using FilteredMyMuons = soa::Filtered; +using FilteredMyMuon = FilteredMyMuons::iterator; + +template +struct DileptonProducer { + Produces eventTable; + Produces normTable; + Produces dileptonTable; + + // Configurables + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; + Configurable d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"}; + + Configurable cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5"}; + Configurable cfgQvecEstimator{"cfgQvecEstimator", 0, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5"}; + Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; + Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; + Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + Configurable cfgApplyWeightTTCA{"cfgApplyWeightTTCA", false, "flag to apply weighting by 1/N"}; + Configurable cfgDCAType{"cfgDCAType", 0, "type of DCA for output. 0:3D, 1:XY, 2:Z, else:3D"}; + + EMEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcut_group"; + Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; + Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + Configurable cfgRequireVertexTOFmatched{"cfgRequireVertexTOFmatched", false, "require Vertex TOFmatched in event cut"}; // ITS-TPC-TOF matched track contributes PV. + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. occupancy"}; + Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. occupancy"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2, "min. FT0C occupancy"}; + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; + Configurable cfgRequireNoCollInTimeRangeStrict{"cfgRequireNoCollInTimeRangeStrict", false, "require no collision in time range strict"}; + Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "require no collision in time range standard"}; + Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "require no collision in time range strict"}; + Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "require no HM collision in previous ITS ROF"}; + Configurable cfgRequireGoodITSLayer3{"cfgRequireGoodITSLayer3", false, "number of inactive chips on ITS layer 3 are below threshold "}; + Configurable cfgRequireGoodITSLayer0123{"cfgRequireGoodITSLayer0123", false, "number of inactive chips on ITS layers 0-3 are below threshold "}; + Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "number of inactive chips on all ITS layers are below threshold "}; + // for RCT + Configurable cfgRequireGoodRCT{"cfgRequireGoodRCT", false, "require good detector flag in run condtion table"}; + Configurable cfgRCTLabel{"cfgRCTLabel", "CBT_hadronPID", "select 1 [CBT, CBT_hadronPID, CBT_muon_glo] see O2Physics/Common/CCDB/RCTSelectionFlags.h"}; + Configurable cfgCheckZDC{"cfgCheckZDC", false, "set ZDC flag for PbPb"}; + Configurable cfgTreatLimitedAcceptanceAsBad{"cfgTreatLimitedAcceptanceAsBad", false, "reject all events where the detectors relevant for the specified Runlist are flagged as LimitedAcceptance"}; + } eventcuts; + + DielectronCut fDielectronCut; + struct : ConfigurableGroup { + std::string prefix = "dielectroncut_group"; + Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; + Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pT"}; + Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pT"}; + Configurable cfg_min_pair_y{"cfg_min_pair_y", -0.8, "min pair rapidity"}; + Configurable cfg_max_pair_y{"cfg_max_pair_y", +0.8, "max pair rapidity"}; + Configurable cfg_min_pair_dca3d{"cfg_min_pair_dca3d", 0.0, "min pair dca3d in sigma"}; + Configurable cfg_max_pair_dca3d{"cfg_max_pair_dca3d", 1e+10, "max pair dca3d in sigma"}; + Configurable cfg_apply_phiv{"cfg_apply_phiv", true, "flag to apply phiv cut"}; + Configurable cfg_phiv_slope{"cfg_phiv_slope", 0.0185, "slope for m vs. phiv"}; + Configurable cfg_phiv_intercept{"cfg_phiv_intercept", -0.0280, "intercept for m vs. phiv"}; + Configurable cfg_min_phiv{"cfg_min_phiv", 0.0, "min phiv (constant)"}; + Configurable cfg_max_phiv{"cfg_max_phiv", 3.2, "max phiv (constant)"}; + Configurable cfg_apply_detadphi{"cfg_apply_detadphi", false, "flag to apply deta-dphi elliptic cut at PV"}; + Configurable cfg_min_deta{"cfg_min_deta", 0.02, "min deta between 2 electrons (elliptic cut)"}; + Configurable cfg_min_dphi{"cfg_min_dphi", 0.2, "min dphi between 2 electrons (elliptic cut)"}; + Configurable cfg_min_opang{"cfg_min_opang", 0.0, "min opening angle"}; + Configurable cfg_max_opang{"cfg_max_opang", 6.4, "max opening angle"}; + Configurable cfg_require_diff_sides{"cfg_require_diff_sides", false, "flag to require 2 tracks are from different sides."}; + + Configurable cfg_apply_cuts_from_prefilter{"cfg_apply_cuts_from_prefilter", false, "flag to apply prefilter set when producing derived data"}; + 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 + + Configurable cfg_apply_cuts_from_prefilter_derived{"cfg_apply_cuts_from_prefilter_derived", false, "flag to apply pair cut same as prefilter set in derived data"}; + 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 + + Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"}; + Configurable cfg_max_pt_track{"cfg_max_pt_track", 1e+10, "max pT for single track"}; + Configurable cfg_min_eta_track{"cfg_min_eta_track", -0.8, "min eta for single track"}; + Configurable cfg_max_eta_track{"cfg_max_eta_track", +0.8, "max eta for single track"}; + Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.f, "min phi for single track"}; + Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for single track"}; + 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"}; + Configurable cfg_reject_phi_track{"cfg_reject_phi_track", false, "reject the phi interval"}; + Configurable cfg_min_ncluster_tpc{"cfg_min_ncluster_tpc", 0, "min ncluster tpc"}; + Configurable cfg_min_ncluster_its{"cfg_min_ncluster_its", 5, "min ncluster its"}; + Configurable cfg_min_ncrossedrows{"cfg_min_ncrossedrows", 100, "min ncrossed rows"}; + Configurable cfg_max_frac_shared_clusters_tpc{"cfg_max_frac_shared_clusters_tpc", 999.f, "max fraction of shared clusters in TPC"}; + Configurable cfg_max_chi2tpc{"cfg_max_chi2tpc", 4.0, "max chi2/NclsTPC"}; + Configurable cfg_max_chi2its{"cfg_max_chi2its", 5.0, "max chi2/NclsITS"}; + Configurable cfg_max_chi2tof{"cfg_max_chi2tof", 1e+10, "max chi2 TOF"}; + Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 0.2, "max dca XY for single track in cm"}; + Configurable cfg_max_dcaz{"cfg_max_dcaz", 0.2, "max dca Z for single track in cm"}; + Configurable cfg_require_itsib_any{"cfg_require_itsib_any", false, "flag to require ITS ib any hits"}; + Configurable cfg_require_itsib_1st{"cfg_require_itsib_1st", true, "flag to require ITS ib 1st hit"}; + Configurable cfg_min_its_cluster_size{"cfg_min_its_cluster_size", 0.f, "min ITS cluster size"}; + Configurable cfg_max_its_cluster_size{"cfg_max_its_cluster_size", 16.f, "max ITS cluster size"}; + Configurable cfg_min_rel_diff_pin{"cfg_min_rel_diff_pin", -1e+10, "min rel. diff. between pin and ppv"}; + Configurable cfg_max_rel_diff_pin{"cfg_max_rel_diff_pin", +1e+10, "max rel. diff. between pin and ppv"}; + Configurable cfgRefR{"cfgRefR", 0.50, "ref. radius (m) for calculating phi position"}; // 0.50 +/- 0.06 can be syst. unc. + Configurable cfg_min_phiposition_track{"cfg_min_phiposition_track", 0.f, "min phi position for single track at certain radius"}; + Configurable cfg_max_phiposition_track{"cfg_max_phiposition_track", 6.3, "max phi position for single track at certain radius"}; + + 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]"}; + Configurable cfg_min_TPCNsigmaEl{"cfg_min_TPCNsigmaEl", -2.0, "min. TPC n sigma for electron inclusion"}; + Configurable cfg_max_TPCNsigmaEl{"cfg_max_TPCNsigmaEl", +3.0, "max. TPC n sigma for electron inclusion"}; + // Configurable cfg_min_TPCNsigmaMu{"cfg_min_TPCNsigmaMu", -0.0, "min. TPC n sigma for muon exclusion"}; + // Configurable cfg_max_TPCNsigmaMu{"cfg_max_TPCNsigmaMu", +0.0, "max. TPC n sigma for muon exclusion"}; + Configurable cfg_min_TPCNsigmaPi{"cfg_min_TPCNsigmaPi", -1e+10, "min. TPC n sigma for pion exclusion"}; + Configurable cfg_max_TPCNsigmaPi{"cfg_max_TPCNsigmaPi", +3.0, "max. TPC n sigma for pion exclusion"}; + Configurable cfg_min_TPCNsigmaKa{"cfg_min_TPCNsigmaKa", -3.0, "min. TPC n sigma for kaon exclusion"}; + Configurable cfg_max_TPCNsigmaKa{"cfg_max_TPCNsigmaKa", +3.0, "max. TPC n sigma for kaon exclusion"}; + Configurable cfg_min_TPCNsigmaPr{"cfg_min_TPCNsigmaPr", -3.0, "min. TPC n sigma for proton exclusion"}; + Configurable cfg_max_TPCNsigmaPr{"cfg_max_TPCNsigmaPr", +3.0, "max. TPC n sigma for proton exclusion"}; + Configurable cfg_min_TOFNsigmaEl{"cfg_min_TOFNsigmaEl", -3.0, "min. TOF n sigma for electron inclusion"}; + Configurable cfg_max_TOFNsigmaEl{"cfg_max_TOFNsigmaEl", +3.0, "max. TOF n sigma for electron inclusion"}; + Configurable cfg_min_pin_pirejTPC{"cfg_min_pin_pirejTPC", 0.f, "min. pin for pion rejection in TPC"}; + Configurable cfg_max_pin_pirejTPC{"cfg_max_pin_pirejTPC", 1e+10, "max. pin for pion rejection in TPC"}; + Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; + Configurable includeITSsa{"includeITSsa", false, "Flag to enable ITSsa tracks"}; + Configurable cfg_max_pt_track_ITSsa{"cfg_max_pt_track_ITSsa", 0.15, "max pt for ITSsa tracks"}; + + // configuration for PID ML + Configurable> onnxFileNames{"onnxFileNames", std::vector{"filename"}, "ONNX file names for each bin (if not from CCDB full path)"}; + Configurable> onnxPathsCCDB{"onnxPathsCCDB", std::vector{"path"}, "Paths of models on CCDB"}; + Configurable> binsMl{"binsMl", std::vector{-999999., 999999.}, "Bin limits for ML application"}; + Configurable> cutsMl{"cutsMl", std::vector{0.95}, "ML cuts per bin"}; + Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature"}, "Names of ML model input features"}; + Configurable nameBinningFeature{"nameBinningFeature", "pt", "Names of ML model binning feature"}; + 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"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + Configurable enableOptimizations{"enableOptimizations", false, "Enables the ONNX extended model-optimization: sessionOptions.SetGraphOptimizationLevel(GraphOptimizationLevel::ORT_ENABLE_EXTENDED)"}; + } dielectroncuts; + + DimuonCut fDimuonCut; + struct : ConfigurableGroup { + std::string prefix = "dimuoncut_group"; + Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; + Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pt"}; + Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pt"}; + Configurable cfg_min_pair_y{"cfg_min_pair_y", -4.0, "min pair rapidity"}; + Configurable cfg_max_pair_y{"cfg_max_pair_y", -2.5, "max pair rapidity"}; + Configurable cfg_min_pair_dcaxy{"cfg_min_pair_dcaxy", 0.0, "min pair dca3d in sigma"}; + Configurable cfg_max_pair_dcaxy{"cfg_max_pair_dcaxy", 1e+10, "max pair dca3d in sigma"}; + Configurable cfg_apply_detadphi{"cfg_apply_detadphi", false, "flag to apply deta-dphi elliptic cut"}; + Configurable cfg_min_deta{"cfg_min_deta", 0.02, "min deta between 2 muons (elliptic cut)"}; + Configurable cfg_min_dphi{"cfg_min_dphi", 0.02, "min dphi between 2 muons (elliptic cut)"}; + + Configurable cfg_track_type{"cfg_track_type", 3, "muon track type [0: MFT-MCH-MID, 3: MCH-MID]"}; + Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"}; + Configurable cfg_max_pt_track{"cfg_max_pt_track", 1e+10, "max pT for single track"}; + Configurable cfg_min_eta_track{"cfg_min_eta_track", -4.0, "min eta for single track"}; + Configurable cfg_max_eta_track{"cfg_max_eta_track", -2.5, "max eta for single track"}; + Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.f, "min phi for single track"}; + Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for single track"}; + Configurable cfg_min_ncluster_mft{"cfg_min_ncluster_mft", 5, "min ncluster MFT"}; + Configurable cfg_min_ncluster_mch{"cfg_min_ncluster_mch", 5, "min ncluster MCH"}; + Configurable cfg_max_chi2{"cfg_max_chi2", 1e+6, "max chi2/ndf"}; + Configurable cfg_max_chi2mft{"cfg_max_chi2mft", 1e+6, "max chi2/ndf"}; + Configurable cfg_max_matching_chi2_mftmch{"cfg_max_matching_chi2_mftmch", 40, "max chi2 for MFT-MCH matching"}; + Configurable cfg_max_matching_chi2_mchmid{"cfg_max_matching_chi2_mchmid", 1e+10, "max chi2 for MCH-MID matching"}; + Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 1e+10, "max dca XY for single track in cm"}; + Configurable cfg_min_rabs{"cfg_min_rabs", 17.6, "min Radius at the absorber end"}; + Configurable cfg_max_rabs{"cfg_max_rabs", 89.5, "max Radius at the absorber end"}; + Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; + Configurable cfg_max_relDPt_wrt_matchedMCHMID{"cfg_max_relDPt_wrt_matchedMCHMID", 1e+10f, "max. relative dpt between MFT-MCH-MID and MCH-MID"}; + Configurable cfg_max_DEta_wrt_matchedMCHMID{"cfg_max_DEta_wrt_matchedMCHMID", 1e+10f, "max. deta between MFT-MCH-MID and MCH-MID"}; + Configurable cfg_max_DPhi_wrt_matchedMCHMID{"cfg_max_DPhi_wrt_matchedMCHMID", 1e+10f, "max. dphi between MFT-MCH-MID and MCH-MID"}; + Configurable requireMFTHitMap{"requireMFTHitMap", false, "flag to apply MFT hit map"}; + Configurable> requiredMFTDisks{"requiredMFTDisks", std::vector{0}, "hit map on MFT disks [0,1,2,3,4]. logical-OR of each double-sided disk"}; + } dimuoncuts; + + o2::aod::rctsel::RCTFlagsChecker rctChecker; + Service ccdb; + int mRunNumber; + float d_bz; + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"}; + + std::mt19937 engine; + float leptonM1 = 0.f; + float leptonM2 = 0.f; + + void init(InitContext& /*context*/) + { + mRunNumber = 0; + d_bz = 0; + + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + rctChecker.init(eventcuts.cfgRCTLabel.value, eventcuts.cfgCheckZDC.value, eventcuts.cfgTreatLimitedAcceptanceAsBad.value); + + std::random_device seed_gen; + engine = std::mt19937(seed_gen()); + + 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; + } + } + + template + void initCCDB(TCollision const& collision) + { + if (mRunNumber == collision.runNumber()) { + return; + } + + // In case override, don't proceed, please - no CCDB access required + if (d_bz_input > -990) { + d_bz = d_bz_input; + o2::parameters::GRPMagField grpmag; + if (std::fabs(d_bz) > 1e-5) { + grpmag.setL3Current(30000.f / (d_bz / 5.0f)); + } + o2::base::Propagator::initFieldFromGRP(&grpmag); + mRunNumber = collision.runNumber(); + return; + } + + auto run3grp_timestamp = collision.timestamp(); + o2::parameters::GRPObject* grpo = 0x0; + o2::parameters::GRPMagField* grpmag = 0x0; + if (!skipGRPOquery) + grpo = ccdb->getForTimeStamp(grpPath, run3grp_timestamp); + 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"; + } + mRunNumber = collision.runNumber(); + fDielectronCut.SetTrackPhiPositionRange(dielectroncuts.cfg_min_phiposition_track, dielectroncuts.cfg_max_phiposition_track, dielectroncuts.cfgRefR, d_bz, dielectroncuts.cfg_mirror_phi_track); + } + + ~DileptonProducer() + { + } + + void addhistograms() + { + } + + 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); + } + + o2::analysis::MlResponseDielectronSingleTrack mlResponseSingleTrack; + 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, false, 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); + + // 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.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.IncludeITSsa(dielectroncuts.includeITSsa, dielectroncuts.cfg_max_pt_track_ITSsa); + + // 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.binsMl.value.size()); + for (size_t i = 0; i < dielectroncuts.binsMl.value.size(); i++) { + binsML.emplace_back(dielectroncuts.binsMl.value[i]); + } + std::vector thresholdsML{}; + thresholdsML.reserve(dielectroncuts.cutsMl.value.size()); + for (size_t i = 0; i < dielectroncuts.cutsMl.value.size(); i++) { + thresholdsML.emplace_back(dielectroncuts.cutsMl.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.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); + } + + template + bool fillPairInfo(TCollision const&, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks) + { + 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; + } + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { + 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, d_bz, 0)) { + 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 (cfgApplyWeightTTCA) { + weight = map_weight[std::make_pair(t1.globalIndex(), t2.globalIndex())]; + } + + 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 v12 = v1 + v2; + + float dca1 = 999.f, dca2 = 999.f; + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + dca1 = dca3DinSigma(t1); + dca2 = dca3DinSigma(t2); + if (cfgDCAType == 1) { + dca1 = dcaXYinSigma(t1); + dca2 = dcaXYinSigma(t2); + } else if (cfgDCAType == 2) { + dca1 = dcaZinSigma(t1); + dca2 = dcaZinSigma(t2); + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + dca1 = fwdDcaXYinSigma(t1); + dca2 = fwdDcaXYinSigma(t2); + } + + // fill table here + dileptonTable(eventTable.lastIndex() + 1, // lastIndex starts from -1. + t1.pt(), t1.eta(), t1.phi(), t1.sign(), dca1, + t2.pt(), t2.eta(), t2.phi(), t2.sign(), dca2, + weight); + + return true; + } + + Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); + // Filter collisionFilter_numContrib = cfgNumContribMin <= o2::aod::collision::numContrib && o2::aod::collision::numContrib < cfgNumContribMax; + Filter collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; + Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; + using FilteredMyCollisions = soa::Filtered; + + SliceCache cache; + Preslice perCollision_electron = aod::emprimaryelectron::emeventId; + 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; + Filter ttcaFilter_electron = ifnode(dielectroncuts.enableTTCA.node(), o2::aod::emprimaryelectron::isAssociatedToMPC == true || o2::aod::emprimaryelectron::isAssociatedToMPC == false, o2::aod::emprimaryelectron::isAssociatedToMPC == true); + 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::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)); + + 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)); + + Partition positive_electrons = o2::aod::emprimaryelectron::sign > int8_t(0); + Partition negative_electrons = o2::aod::emprimaryelectron::sign < int8_t(0); + + Preslice perCollision_muon = aod::emprimarymuon::emeventId; + 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; + Filter ttcaFilter_muon = ifnode(dimuoncuts.enableTTCA.node(), o2::aod::emprimarymuon::isAssociatedToMPC == true || o2::aod::emprimarymuon::isAssociatedToMPC == false, o2::aod::emprimarymuon::isAssociatedToMPC == true); + Partition positive_muons = o2::aod::emprimarymuon::sign > int8_t(0); + Partition negative_muons = o2::aod::emprimarymuon::sign < int8_t(0); + + 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] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + float eventplanes_2_for_mix[6] = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg()}; + float ep2 = eventplanes_2_for_mix[cfgEP2Estimator_for_Mix]; + + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + if (eventcuts.cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { + continue; + } + + auto posTracks_per_coll = posTracks.sliceByCached(perCollision, collision.globalIndex(), cache); + auto negTracks_per_coll = negTracks.sliceByCached(perCollision, collision.globalIndex(), cache); + + int nuls = 0, nlspp = 0, nlsmm = 0; + for (const auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS + bool is_pair_ok = fillPairInfo(collision, pos, neg, cut, tracks); + if (is_pair_ok) { + nuls++; + } + } + for (const auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ + bool is_pair_ok = fillPairInfo(collision, pos1, pos2, cut, tracks); + if (is_pair_ok) { + nlspp++; + } + } + for (const auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- + bool is_pair_ok = fillPairInfo(collision, neg1, neg2, cut, tracks); + if (is_pair_ok) { + nlsmm++; + } + } + + if (nuls > 0 || nlspp > 0 || nlsmm > 0) { + eventTable(collision.runNumber(), collision.globalBC(), collision.timestamp(), collision.posZ(), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange(), collision.centFT0C(), ep2); + } + } // end of collision loop + } // end of DF + + template + bool isPairOK(TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks) + { + 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; + } + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + if (!cut.IsSelectedTrack(t1) || !cut.IsSelectedTrack(t2)) { + 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, d_bz, 0.0)) { + 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] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + if (eventcuts.cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { + 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(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(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(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)); + // LOGF(info, "std::get<0>(pairId) = %d, std::get<1>(pairId) = %d, t1.globalIndex() = %d, t2.globalIndex() = %d", std::get<0>(pairId), std::get<1>(pairId), t1.globalIndex(), t2.globalIndex()); + + 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...)); + 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(); + ndf++; + } + PROCESS_SWITCH(DileptonProducer, processAnalysis, "run dilepton analysis", true); + + Filter collisionFilter_centrality_norm = cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax; + using FilteredNormInfos = soa::Filtered; + void processNorm(FilteredNormInfos const& collisions) + { + for (const auto& collision : collisions) { + if (collision.centFT0C() < cfgCentMin || cfgCentMax < collision.centFT0C()) { + continue; + } + + normTable(collision.alias_raw(), collision.selection_raw(), collision.rct_raw(), collision.posZ(), collision.centFT0C()); + + } // end of collision loop + } + PROCESS_SWITCH(DileptonProducer, processNorm, "process normalization info", true); +}; + +#endif // PWGEM_DILEPTON_CORE_DILEPTONPRODUCER_H_ diff --git a/PWGEM/Dilepton/DataModel/dileptonTables.h b/PWGEM/Dilepton/DataModel/dileptonTables.h index 7164dba9ece..d116b92485d 100644 --- a/PWGEM/Dilepton/DataModel/dileptonTables.h +++ b/PWGEM/Dilepton/DataModel/dileptonTables.h @@ -293,7 +293,7 @@ DECLARE_SOA_TABLE(EMEoIs, "AOD", "EMEOI", //! joinable to aod::Collisions in cre using EMEoI = EMEoIs::iterator; DECLARE_SOA_TABLE(EMEventNormInfos, "AOD", "EMEVENTNORMINFO", //! event information for normalization - o2::soa::Index<>, evsel::Alias, evsel::Selection, evsel::Rct, emevent::PosZint16, cent::CentFT0C, emevent::PosZ, emevent::Sel8); + o2::soa::Index<>, evsel::Alias, evsel::Selection, evsel::Rct, emevent::PosZint16, cent::CentFT0C, emevent::PosZ, emevent::Sel8, o2::soa::Marker<1>); using EMEventNormInfo = EMEventNormInfos::iterator; namespace emmcevent @@ -937,6 +937,44 @@ DECLARE_SOA_TABLE(EMPrimaryTrackEMEventIdsTMP, "AOD", "PRMTRKEVIDTMP", track::Co // iterators using EMPrimaryTrackEMEventIdTMP = EMPrimaryTrackEMEventIdsTMP::iterator; +namespace emthinevent +{ +DECLARE_SOA_COLUMN(EP2, ep2, float); //! +} // namespace emthinevent +DECLARE_SOA_TABLE_VERSIONED(EMThinEvents_000, "AOD", "EMTHINEVENT", 0, //! Thin event information table + o2::soa::Index<>, bc::RunNumber, bc::GlobalBC, timestamp::Timestamp, collision::PosZ, + evsel::NumTracksInTimeRange, evsel::SumAmpFT0CInTimeRange, cent::CentFT0C, emthinevent::EP2); +using EMThinEvents = EMThinEvents_000; +using EMThinEvent = EMThinEvents::iterator; + +DECLARE_SOA_TABLE(EMThinEventNormInfos, "AOD", "EMTHINEVENTNORM", //! event information for normalization + o2::soa::Index<>, evsel::Alias, evsel::Selection, evsel::Rct, collision::PosZ, cent::CentFT0C, emevent::Sel8, o2::soa::Marker<2>); +using EMThinEventNormInfo = EMThinEventNormInfos::iterator; + +namespace emdilepton +{ +DECLARE_SOA_INDEX_COLUMN(EMThinEvent, emthinevent); //! +DECLARE_SOA_COLUMN(Pt1, pt1, float); //! +DECLARE_SOA_COLUMN(Eta1, eta1, float); //! +DECLARE_SOA_COLUMN(Phi1, phi1, float); //! +DECLARE_SOA_COLUMN(Sign1, sign1, short); //! +DECLARE_SOA_COLUMN(DCA1, dca1, float); //! DCA in sigma. Users should decide 3D or XY or Z +DECLARE_SOA_COLUMN(Pt2, pt2, float); //! +DECLARE_SOA_COLUMN(Eta2, eta2, float); //! +DECLARE_SOA_COLUMN(Phi2, phi2, float); //! +DECLARE_SOA_COLUMN(DCA2, dca2, float); //! DCA in sigma. Users should decide 3D or XY or Z +DECLARE_SOA_COLUMN(Sign2, sign2, short); //! +DECLARE_SOA_COLUMN(Weight, weight, float); //! possible pair weight +} // namespace emdilepton + +DECLARE_SOA_TABLE_VERSIONED(EMDileptons_000, "AOD", "EMDILEPTON", 0, + o2::soa::Index<>, emdilepton::EMThinEventId, + emdilepton::Pt1, emdilepton::Eta1, emdilepton::Phi1, emdilepton::Sign1, emdilepton::DCA1, + emdilepton::Pt2, emdilepton::Eta2, emdilepton::Phi2, emdilepton::Sign2, emdilepton::DCA2, + emdilepton::Weight); +using EMDileptons = EMDileptons_000; +using EMDilepton = EMDileptons::iterator; + // Dummy data for MC namespace emdummydata { diff --git a/PWGEM/Dilepton/TableProducer/CMakeLists.txt b/PWGEM/Dilepton/TableProducer/CMakeLists.txt index 38a01c9af2e..eba934d8b31 100644 --- a/PWGEM/Dilepton/TableProducer/CMakeLists.txt +++ b/PWGEM/Dilepton/TableProducer/CMakeLists.txt @@ -64,3 +64,13 @@ o2physics_add_dpl_workflow(filter-eoi SOURCES filterEoI.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(dielectron-producer + SOURCES dielectronProducer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMDileptonCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(dimuon-producer + SOURCES dimuonProducer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMDileptonCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/Dilepton/TableProducer/dielectronProducer.cxx b/PWGEM/Dilepton/TableProducer/dielectronProducer.cxx new file mode 100644 index 00000000000..ed9b5047afc --- /dev/null +++ b/PWGEM/Dilepton/TableProducer/dielectronProducer.cxx @@ -0,0 +1,26 @@ +// 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/DileptonProducer.h" + +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask>(cfgc, TaskName{"dielectron-producer"})}; +} diff --git a/PWGEM/Dilepton/TableProducer/dimuonProducer.cxx b/PWGEM/Dilepton/TableProducer/dimuonProducer.cxx new file mode 100644 index 00000000000..32cd1bafa72 --- /dev/null +++ b/PWGEM/Dilepton/TableProducer/dimuonProducer.cxx @@ -0,0 +1,26 @@ +// 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 dimuon analyses. +// Please write to: daiki.sekihata@cern.ch + +#include "PWGEM/Dilepton/Core/DileptonProducer.h" + +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask>(cfgc, TaskName{"dimuon-producer"})}; +} diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt b/PWGEM/Dilepton/Tasks/CMakeLists.txt index cf87a9201f8..61779f3fde4 100644 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -175,3 +175,8 @@ o2physics_add_dpl_workflow(evaluate-acceptance SOURCES evaluateAcceptance.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(dilepton-polarization + SOURCES dileptonPolarization.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/Dilepton/Tasks/dileptonPolarization.cxx b/PWGEM/Dilepton/Tasks/dileptonPolarization.cxx new file mode 100644 index 00000000000..3e6da01168b --- /dev/null +++ b/PWGEM/Dilepton/Tasks/dileptonPolarization.cxx @@ -0,0 +1,842 @@ +// 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 + +#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 "CCDB/BasicCCDBManager.h" +#include "CommonConstants/LHCConstants.h" +#include "DataFormatsParameters/GRPECSObject.h" +#include "DataFormatsParameters/GRPLHCIFData.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "MathUtils/Utils.h" + +#include "Math/Vector4D.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::dilepton::utils; +using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; +using namespace o2::aod::pwgem::dilepton::utils::pairutil; + +using MyEMH_electron = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, EMTrack>; +using MyEMH_muon = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, EMFwdTrack>; +using MyEMH_pair = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, std::tuple>; + +struct DileptonPolarization { + // Configurables + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; + Configurable d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"}; + + Configurable cfgPairType{"cfgPairType", 0, "0:dielectron:0, 1:dimuon"}; + Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; + Configurable cfgDoMix{"cfgDoMix", true, "flag for event mixing"}; + Configurable ndepth{"ndepth", 100, "depth for event mixing"}; + Configurable ndiff_bc_mix{"ndiff_bc_mix", 594, "difference in global BC required in mixed events"}; + ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f, 999.f}, "Mixing bins - centrality"}; + ConfigurableAxis ConfEPBins{"ConfEPBins", {16, -M_PI / 2, +M_PI / 2}, "Mixing bins - event plane angle"}; + ConfigurableAxis ConfOccupancyBins{"ConfOccupancyBins", {VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; + Configurable cfgPolarizationFrame{"cfgPolarizationFrame", 0, "frame of polarization. 0:CS, 1:HX, else:FATAL"}; + + ConfigurableAxis ConfMllBins{"ConfMllBins", {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, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.00, 8.10, 8.20, 8.30, 8.40, 8.50, 8.60, 8.70, 8.80, 8.90, 9.00, 9.10, 9.20, 9.30, 9.40, 9.50, 9.60, 9.70, 9.80, 9.90, 10.00, 10.10, 10.20, 10.30, 10.40, 10.50, 10.60, 10.70, 10.80, 10.90, 11.00, 11.1, 11.2, 11.3, 11.4, 11.50, 11.6, 11.7, 11.8, 11.9, 12.0}, "mll bins for output histograms"}; + ConfigurableAxis ConfPtllBins{"ConfPtllBins", {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"}; + ConfigurableAxis ConfDCAllBins{"ConfDCAllBins", {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"}; + ConfigurableAxis ConfYllBins{"ConYllBins", {1, -1.f, 1.f}, "yll bins for output histograms"}; // pair rapidity + + ConfigurableAxis ConfPolarizationCosThetaBins{"ConfPolarizationCosThetaBins", {20, -1.f, 1.f}, "cos(theta) bins for polarization analysis"}; + ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, -M_PI, M_PI}, "phi bins for polarization analysis"}; + ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {15, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2> + + struct : ConfigurableGroup { + std::string prefix = "eventcut_group"; + Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; + Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; + // Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + // Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + // Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + // Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + // Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + // Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + // Configurable cfgRequireVertexTOFmatched{"cfgRequireVertexTOFmatched", false, "require Vertex TOFmatched in event cut"}; // ITS-TPC-TOF matched track contributes PV. + // Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. occupancy"}; + Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. occupancy"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2, "min. FT0C occupancy"}; + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + } eventcuts; + + struct : ConfigurableGroup { + std::string prefix = "dileptoncut_group"; + Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; + Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pT"}; + Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pT"}; + Configurable cfg_min_pair_y{"cfg_min_pair_y", -0.9, "min pair rapidity"}; + Configurable cfg_max_pair_y{"cfg_max_pair_y", +0.9, "max pair rapidity"}; + + Configurable cfg_min_track_pt{"cfg_min_track_pt", 0.2, "min pT for single track"}; + Configurable cfg_max_track_pt{"cfg_max_track_pt", 1e+10, "max pT for single track"}; + Configurable cfg_min_track_eta{"cfg_min_track_eta", -0.9, "min eta for single track"}; + Configurable cfg_max_track_eta{"cfg_max_track_eta", +0.9, "max eta for single track"}; + } dileptoncuts; + + struct : ConfigurableGroup { + std::string prefix = "accBins"; + ConfigurableAxis ConfMllAccBins{"ConfMllAccBins", {40, 0, 4}, "mll bins for acceptance for plarization"}; + ConfigurableAxis ConfPtllAccBins{"ConfPtllAccBins", {100, 0, 10}, "pTll bins for acceptance for plarization"}; + ConfigurableAxis ConfEtallAccBins{"ConEtallAccBins", {30, -1.5f, 1.5f}, "etall bins for acceptance for plarization"}; // pair pseudo-rapidity + ConfigurableAxis ConfPhillAccBins{"ConPhillAccBins", {36, 0.f, 2 * M_PI}, "phill bins for acceptance for plarization"}; // pair pseudo-rapidity + } accBins; + + Service ccdb; + int mRunNumber; + float d_bz; + + HistogramRegistry fRegistry{"output", {}, 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; + std::vector mll_bin_edges; + std::vector ptll_bin_edges; + std::vector etall_bin_edges; + std::vector phill_bin_edges; + + 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 + + void init(InitContext& /*context*/) + { + mRunNumber = 0; + d_bz = 0; + + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + if (cfgPairType.value == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron)) { + leptonM1 = o2::constants::physics::MassElectron; + leptonM2 = o2::constants::physics::MassElectron; + } else if (cfgPairType.value == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon)) { + leptonM1 = o2::constants::physics::MassMuon; + leptonM2 = o2::constants::physics::MassMuon; + } + + if (ConfVtxBins.value[0] == 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, "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] == 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, "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] == 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, "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] == 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, "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_pair_uls = new MyEMH_pair(ndepth); + emh_pair_lspp = new MyEMH_pair(ndepth); + emh_pair_lsmm = new MyEMH_pair(ndepth); + + if (accBins.ConfMllAccBins.value[0] == VARIABLE_WIDTH) { + mll_bin_edges = std::vector(accBins.ConfMllAccBins.value.begin(), accBins.ConfMllAccBins.value.end()); + mll_bin_edges.erase(mll_bin_edges.begin()); + for (const auto& edge : mll_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: mll_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfMllAccBins.value[0]); + float xmin = static_cast(accBins.ConfMllAccBins.value[1]); + float xmax = static_cast(accBins.ConfMllAccBins.value[2]); + mll_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + mll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: mll_bin_edges[%d] = %f", i, mll_bin_edges[i]); + } + } + + if (accBins.ConfPtllAccBins.value[0] == VARIABLE_WIDTH) { + ptll_bin_edges = std::vector(accBins.ConfPtllAccBins.value.begin(), accBins.ConfPtllAccBins.value.end()); + ptll_bin_edges.erase(ptll_bin_edges.begin()); + for (const auto& edge : ptll_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: ptll_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfPtllAccBins.value[0]); + float xmin = static_cast(accBins.ConfPtllAccBins.value[1]); + float xmax = static_cast(accBins.ConfPtllAccBins.value[2]); + ptll_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + ptll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: ptll_bin_edges[%d] = %f", i, ptll_bin_edges[i]); + } + } + + if (accBins.ConfEtallAccBins.value[0] == VARIABLE_WIDTH) { + etall_bin_edges = std::vector(accBins.ConfEtallAccBins.value.begin(), accBins.ConfEtallAccBins.value.end()); + etall_bin_edges.erase(etall_bin_edges.begin()); + for (const auto& edge : etall_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: etall_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfEtallAccBins.value[0]); + float xmin = static_cast(accBins.ConfEtallAccBins.value[1]); + float xmax = static_cast(accBins.ConfEtallAccBins.value[2]); + etall_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + etall_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: etall_bin_edges[%d] = %f", i, etall_bin_edges[i]); + } + } + + if (accBins.ConfPhillAccBins.value[0] == VARIABLE_WIDTH) { + phill_bin_edges = std::vector(accBins.ConfPhillAccBins.value.begin(), accBins.ConfPhillAccBins.value.end()); + phill_bin_edges.erase(phill_bin_edges.begin()); + for (const auto& edge : phill_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: phill_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfPhillAccBins.value[0]); + float xmin = static_cast(accBins.ConfPhillAccBins.value[1]); + float xmax = static_cast(accBins.ConfPhillAccBins.value[2]); + phill_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + phill_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: phill_bin_edges[%d] = %f", i, phill_bin_edges[i]); + } + } + + std::random_device seed_gen; + engine = std::mt19937(seed_gen()); + + addhistograms(); + + fRegistry.add("Pair/mix/hDiffBC", "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", kTH1D, {{10001, -0.5, 10000.5}}, true); + } + + template + void initCCDB(TCollision const& collision) + { + if (mRunNumber == collision.runNumber()) { + return; + } + + // In case override, don't proceed, please - no CCDB access required + if (d_bz_input > -990) { + d_bz = d_bz_input; + o2::parameters::GRPMagField grpmag; + if (std::fabs(d_bz) > 1e-5) { + grpmag.setL3Current(30000.f / (d_bz / 5.0f)); + } + o2::base::Propagator::initFieldFromGRP(&grpmag); + mRunNumber = collision.runNumber(); + return; + } + + auto run3grp_timestamp = collision.timestamp(); + o2::parameters::GRPObject* grpo = 0x0; + o2::parameters::GRPMagField* grpmag = 0x0; + if (!skipGRPOquery) + grpo = ccdb->getForTimeStamp(grpPath, run3grp_timestamp); + 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"; + } + mRunNumber = collision.runNumber(); + + 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); + } + + ~DileptonPolarization() + { + delete emh_pair_uls; + emh_pair_uls = 0x0; + delete emh_pair_lspp; + emh_pair_lspp = 0x0; + delete emh_pair_lsmm; + emh_pair_lsmm = 0x0; + + map_mixed_eventId_to_globalBC.clear(); + } + + void addhistograms() + { + auto hCollisionCounter = fRegistry.add("Event/before/hCollisionCounter", "collision counter;;Number of events", kTH1D, {{9, 0.5, 9 + 0.5}}, false); + hCollisionCounter->GetXaxis()->SetBinLabel(1, "all"); + hCollisionCounter->GetXaxis()->SetBinLabel(2, "FT0AND"); + hCollisionCounter->GetXaxis()->SetBinLabel(3, "No TF border"); + hCollisionCounter->GetXaxis()->SetBinLabel(4, "No ITS ROF border"); + hCollisionCounter->GetXaxis()->SetBinLabel(5, "No Same Bunch Pileup"); + hCollisionCounter->GetXaxis()->SetBinLabel(6, "Is Good Zvtx FT0vsPV"); + hCollisionCounter->GetXaxis()->SetBinLabel(7, "sel8"); + hCollisionCounter->GetXaxis()->SetBinLabel(8, "|Z_{vtx}| < 10 cm"); + hCollisionCounter->GetXaxis()->SetBinLabel(9, "accepted"); + + fRegistry.add("Event/before/hZvtx", "vertex z; Z_{vtx} (cm)", kTH1D, {{100, -50, +50}}, false); + fRegistry.add("Event/before/hCentFT0C", "hCentFT0C;centrality FT0C (%)", kTH1D, {{110, 0, 110}}, false); + fRegistry.add("Event/before/hCorrOccupancy", "occupancy correlation;FT0C occupancy;track occupancy", kTH2D, {{200, 0, 200000}, {200, 0, 20000}}, false); + fRegistry.add("Event/before/hEP2_CentFT0C_forMix", "2nd harmonics event plane for mix;centrality FT0C (%);#Psi_{2} (rad.)", kTH2F, {{110, 0, 110}, {180, -M_PI_2, +M_PI_2}}, false); + fRegistry.addClone("Event/before/", "Event/after/"); + + 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 (cfgPairType == static_cast(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_dca_axis_title = "DCA_{ee} (#sigma)"; + pair_y_axis_title = "y_{ee}"; + } else if (cfgPairType == static_cast(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} (#sigma)"; + pair_y_axis_title = "y_{#mu#mu}"; + } + + // pair info + const AxisSpec axis_mass{ConfMllBins, mass_axis_title}; + const AxisSpec axis_pt{ConfPtllBins, pair_pt_axis_title}; + const AxisSpec axis_dca{ConfDCAllBins, pair_dca_axis_title}; + const AxisSpec axis_y{ConfYllBins, pair_y_axis_title}; + + 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 AxisSpec axis_cos_theta{ConfPolarizationCosThetaBins, Form("cos(#theta^{%s})", frameName.data())}; + const AxisSpec axis_phi{ConfPolarizationPhiBins, Form("#varphi^{%s} (rad.)", frameName.data())}; + const AxisSpec axis_quadmom{ConfPolarizationQuadMomBins, Form("#frac{3 cos^{2}(#theta^{%s}) -1}{2}", frameName.data())}; + fRegistry.add("Pair/same/uls/hs", "dilepton", 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/"); + } + + template + bool fillPairInfo(TCollision const& collision, TDilepton const& dilepton) + { + float weight = 1.f; + ROOT::Math::PtEtaPhiMVector v1(dilepton.pt1(), dilepton.eta1(), dilepton.phi1(), leptonM1); + ROOT::Math::PtEtaPhiMVector v2(dilepton.pt2(), dilepton.eta2(), dilepton.phi2(), leptonM2); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + if (v12.Rapidity() < dileptoncuts.cfg_min_pair_y || dileptoncuts.cfg_max_pair_y < v12.Rapidity()) { + return false; + } + + float pair_dca = pairDCAQuadSum(dilepton.dca1(), dilepton.dca2()); + 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 = dilepton.sign1() * dilepton.sign2() < 0 ? (dilepton.sign1() > 0 ? std::array{static_cast(v1.Px()), static_cast(v1.Py()), static_cast(v1.Pz()), leptonM1} : std::array{static_cast(v2.Px()), static_cast(v2.Py()), static_cast(v2.Pz()), leptonM2}) : (random_sign > 0 ? std::array{static_cast(v1.Px()), static_cast(v1.Py()), static_cast(v1.Pz()), leptonM1} : std::array{static_cast(v2.Px()), static_cast(v2.Py()), static_cast(v2.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 (dilepton.sign1() * dilepton.sign2() < 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 (dilepton.sign1() > 0 && dilepton.sign2() > 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 (dilepton.sign1() < 0 && dilepton.sign2() < 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); + } + + if constexpr (ev_id == 0) { // same event + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), v12.M()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), v12.Pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), v12.Eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + float phi12 = v12.Phi(); + o2::math_utils::bringTo02Pi(phi12); + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), phi12) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + auto key_df_collision = std::make_pair(ndf, collision.globalIndex()); + float phi12_tmp = v12.Phi(); + o2::math_utils::bringTo02Pi(phi12_tmp); + EMPair empair = EMPair(v12.Pt(), v12.Eta(), phi12_tmp, v12.M(), 0); + empair.setPositiveLegPxPyPzM(arrD[0], arrD[1], arrD[2], leptonM1); + // empair.setNegativeLegPtEtaPhiM(t2.pt(), t2.eta(), t2.phi(), leptonM2); + empair.setPairDCA(pair_dca); + auto pair_tmp = std::make_tuple(mbin, ptbin, etabin, phibin, empair); + if (dilepton.sign1() * dilepton.sign2() < 0) { // ULS + emh_pair_uls->AddTrackToEventPool(key_df_collision, pair_tmp); + } else if (dilepton.sign1() > 0 && dilepton.sign2() > 0) { // LS++ + emh_pair_lspp->AddTrackToEventPool(key_df_collision, pair_tmp); + } else if (dilepton.sign1() < 0 && dilepton.sign2() < 0) { // LS-- + emh_pair_lsmm->AddTrackToEventPool(key_df_collision, pair_tmp); + } + } + return true; + } + + Filter collisionFilter_centrality = eventcuts.cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < eventcuts.cfgCentMax; + Filter collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; + Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; + using filteredCollisions = soa::Filtered; + + Filter dileptonFilter_track1 = (dileptoncuts.cfg_min_track_pt < o2::aod::emdilepton::pt1 && o2::aod::emdilepton::pt1 < dileptoncuts.cfg_max_track_pt) && (dileptoncuts.cfg_min_track_eta < o2::aod::emdilepton::eta1 && o2::aod::emdilepton::eta1 < dileptoncuts.cfg_max_track_eta); + Filter dileptonFilter_track2 = (dileptoncuts.cfg_min_track_pt < o2::aod::emdilepton::pt2 && o2::aod::emdilepton::pt2 < dileptoncuts.cfg_max_track_pt) && (dileptoncuts.cfg_min_track_eta < o2::aod::emdilepton::eta2 && o2::aod::emdilepton::eta2 < dileptoncuts.cfg_max_track_eta); + using filteredDileptons = soa::Filtered; + + SliceCache cache; + Preslice perCollision = aod::emdilepton::emthineventId; + Partition dileptonsULS = (o2::aod::emdilepton::sign1 > int16_t(0) && o2::aod::emdilepton::sign2 < int16_t(0)) || (o2::aod::emdilepton::sign1 < int16_t(0) && o2::aod::emdilepton::sign2 > int16_t(0)); + Partition dileptonsLSPP = o2::aod::emdilepton::sign1 > int16_t(0) && o2::aod::emdilepton::sign2 > int16_t(0); + Partition dileptonsLSMM = o2::aod::emdilepton::sign1 < int16_t(0) && o2::aod::emdilepton::sign2 < int16_t(0); + + MyEMH_pair* emh_pair_uls = nullptr; + MyEMH_pair* emh_pair_lspp = nullptr; + MyEMH_pair* emh_pair_lsmm = nullptr; + + std::map, uint64_t> map_mixed_eventId_to_globalBC; + + int ndf = 0; + + template + void runPairing(TCollisions const& collisions, TDileptons const&) + { + for (const auto& collision : collisions) { + initCCDB(collision); + float centrality = collision.centFT0C(); + if (centrality < eventcuts.cfgCentMin || eventcuts.cfgCentMax < centrality) { + continue; + } + + float ep2 = collision.ep2(); + fRegistry.fill(HIST("Event/after/hZvtx"), collision.posZ()); + fRegistry.fill(HIST("Event/after/hCollisionCounter"), 9); // is qvector calibarated + fRegistry.fill(HIST("Event/after/hCorrOccupancy"), collision.ft0cOccupancyInTimeRange(), collision.trackOccupancyInTimeRange()); + fRegistry.fill(HIST("Event/after/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); + + auto dileptons_uls_per_coll = dileptonsULS->sliceByCached(aod::emdilepton::emthineventId, collision.globalIndex(), cache); + auto dileptons_lspp_per_coll = dileptonsLSPP->sliceByCached(aod::emdilepton::emthineventId, collision.globalIndex(), cache); + auto dileptons_lsmm_per_coll = dileptonsLSMM->sliceByCached(aod::emdilepton::emthineventId, collision.globalIndex(), cache); + // LOGF(info, "collision.globalIndex() = %d, dileptons_uls_per_coll.size() = %d, dileptons_lspp_per_coll.size() = %d, dileptons_lsmm_per_coll.size() = %d", collision.globalIndex(), dileptons_uls_per_coll.size(), dileptons_lspp_per_coll.size(), dileptons_lsmm_per_coll.size()); + + int nuls = 0, nlspp = 0, nlsmm = 0; + for (const auto& dilepton : dileptons_uls_per_coll) { // ULS + bool is_pair_ok = fillPairInfo<0>(collision, dilepton); + if (is_pair_ok) { + nuls++; + } + } + for (const auto& dilepton : dileptons_lspp_per_coll) { // LS++ + bool is_pair_ok = fillPairInfo<0>(collision, dilepton); + if (is_pair_ok) { + nlspp++; + } + } + for (const auto& dilepton : dileptons_lsmm_per_coll) { // LS-- + bool is_pair_ok = fillPairInfo<0>(collision, dilepton); + if (is_pair_ok) { + nlsmm++; + } + } + + 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.ft0cOccupancyInTimeRange() = %f, zbin = %d, centbin = %d, epbin = %d, occbin = %d", collision.globalIndex(), collision.posZ(), centrality, ep2, collision.ft0cOccupancyInTimeRange(), 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. + + auto selected_pairs_uls_in_this_event = emh_pair_uls->GetTracksPerCollision(key_df_collision); + auto selected_pairs_lspp_in_this_event = emh_pair_lspp->GetTracksPerCollision(key_df_collision); + auto selected_pairs_lsmm_in_this_event = emh_pair_lsmm->GetTracksPerCollision(key_df_collision); + auto collisionIds_in_mixing_pool = emh_pair_uls->GetCollisionIdsFromEventPool(key_bin); + float weight = 1.f; + + for (const auto& mix_dfId_collisionId : collisionIds_in_mixing_pool) { + auto pairs_uls_from_event_pool = emh_pair_uls->GetTracksPerCollision(mix_dfId_collisionId); + auto pairs_lspp_from_event_pool = emh_pair_lspp->GetTracksPerCollision(mix_dfId_collisionId); + auto pairs_lsmm_from_event_pool = emh_pair_lsmm->GetTracksPerCollision(mix_dfId_collisionId); + + 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; + } + + for (const auto& pair1 : selected_pairs_uls_in_this_event) { // ULS mix + auto empair1 = std::get<4>(pair1); + auto v_pos = empair1.getPositiveLeg(); // pt, eta, phi, M + auto arrD = std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1}; + + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), empair1.mass()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), empair1.pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), empair1.eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), empair1.phi()) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + for (const auto& pair2 : std::views::filter(pairs_uls_from_event_pool, [&mbin, &ptbin, &etabin, &phibin](std::tuple t) { return std::get<0>(t) == mbin && std::get<1>(t) == ptbin && std::get<2>(t) == etabin && std::get<3>(t) == phibin; })) { + auto empair2 = std::get<4>(pair2); + auto arrM = std::array{static_cast(empair2.px()), static_cast(empair2.py()), static_cast(empair2.pz()), static_cast(empair2.mass())}; + + float cos_thetaPol = 999, phiPol = 999.f; + 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; + fRegistry.fill(HIST("Pair/mix/uls/hs"), empair1.mass(), empair1.pt(), empair1.getPairDCA(), empair1.rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } + } // end of ULS + + for (const auto& pair1 : selected_pairs_lspp_in_this_event) { // LS++ + auto empair1 = std::get<4>(pair1); + auto v_pos = empair1.getPositiveLeg(); // pt, eta, phi, M + auto arrD = std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1}; + + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), empair1.mass()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), empair1.pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), empair1.eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), empair1.phi()) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + for (const auto& pair2 : std::views::filter(pairs_lspp_from_event_pool, [&mbin, &ptbin, &etabin, &phibin](std::tuple t) { return std::get<0>(t) == mbin && std::get<1>(t) == ptbin && std::get<2>(t) == etabin && std::get<3>(t) == phibin; })) { + auto empair2 = std::get<4>(pair2); + auto arrM = std::array{static_cast(empair2.px()), static_cast(empair2.py()), static_cast(empair2.pz()), static_cast(empair2.mass())}; + + float cos_thetaPol = 999, phiPol = 999.f; + 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; + fRegistry.fill(HIST("Pair/mix/lspp/hs"), empair1.mass(), empair1.pt(), empair1.getPairDCA(), empair1.rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } + } // end of LS++ + + for (const auto& pair1 : selected_pairs_lsmm_in_this_event) { // LS-- + auto empair1 = std::get<4>(pair1); + auto v_pos = empair1.getPositiveLeg(); // pt, eta, phi, M + auto arrD = std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1}; + + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), empair1.mass()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), empair1.pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), empair1.eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), empair1.phi()) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + for (const auto& pair2 : std::views::filter(pairs_lsmm_from_event_pool, [&mbin, &ptbin, &etabin, &phibin](std::tuple t) { return std::get<0>(t) == mbin && std::get<1>(t) == ptbin && std::get<2>(t) == etabin && std::get<3>(t) == phibin; })) { + auto empair2 = std::get<4>(pair2); + auto arrM = std::array{static_cast(empair2.px()), static_cast(empair2.py()), static_cast(empair2.pz()), static_cast(empair2.mass())}; + + float cos_thetaPol = 999, phiPol = 999.f; + 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; + fRegistry.fill(HIST("Pair/mix/lsmm/hs"), empair1.mass(), empair1.pt(), empair1.getPairDCA(), empair1.rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } + } // end of LS-- + + } // 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_pair_uls->AddCollisionIdAtLast(key_bin, key_df_collision); + emh_pair_lspp->AddCollisionIdAtLast(key_bin, key_df_collision); + emh_pair_lsmm->AddCollisionIdAtLast(key_bin, key_df_collision); + } // end of if pair exist + + } // end of collision loop + + } // end of DF + + void processAnalysis(filteredCollisions const& collisions, filteredDileptons const& dileptons) + { + runPairing(collisions, dileptons); + ndf++; + } + PROCESS_SWITCH(DileptonPolarization, processAnalysis, "run dilepton analysis", true); + + void processDummy(aod::EMThinEvents const&) {} + PROCESS_SWITCH(DileptonPolarization, processDummy, "Dummy function", false); +}; +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"dilepton-polarization"})}; +}