diff --git a/PWGEM/Dilepton/Core/DileptonProducer.h b/PWGEM/Dilepton/Core/DileptonProducer.h deleted file mode 100644 index 25e25fb3e3e..00000000000 --- a/PWGEM/Dilepton/Core/DileptonProducer.h +++ /dev/null @@ -1,782 +0,0 @@ -// 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/EMTrackUtilities.h" -#include "PWGEM/Dilepton/Utils/PairUtilities.h" - -#include "Common/CCDB/RCTSelectionFlags.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/PIDResponseTPC.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using MyCollisions = o2::soa::Join; -using MyCollision = MyCollisions::iterator; - -using MyElectrons = o2::soa::Join; -using MyElectron = MyElectrons::iterator; -using FilteredMyElectrons = o2::soa::Filtered; -using FilteredMyElectron = FilteredMyElectrons::iterator; - -using MyMuons = o2::soa::Join; -using MyMuon = MyMuons::iterator; -using FilteredMyMuons = o2::soa::Filtered; -using FilteredMyMuon = FilteredMyMuons::iterator; - -template -struct DileptonProducer { - o2::framework::Produces eventTable; - o2::framework::Produces normTable; - o2::framework::Produces dileptonTable; - - // Configurables - o2::framework::Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - o2::framework::Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - o2::framework::Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - o2::framework::Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; - o2::framework::Configurable d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"}; - - o2::framework::Configurable cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5"}; - o2::framework::Configurable cfgQvecEstimator{"cfgQvecEstimator", 0, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5"}; - o2::framework::Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; - o2::framework::Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; - o2::framework::Configurable cfgApplyWeightTTCA{"cfgApplyWeightTTCA", false, "flag to apply weighting by 1/N"}; - o2::framework::Configurable cfgDCAType{"cfgDCAType", 0, "type of DCA for output. 0:3D, 1:XY, 2:Z, else:3D"}; - o2::framework::Configurable cfgStoreULS{"cfgStoreULS", true, "flag to store ULS pairs"}; - o2::framework::Configurable cfgStoreLS{"cfgStoreLS", true, "flag to store LS pairs"}; - - EMEventCut fEMEventCut; - struct : o2::framework::ConfigurableGroup { - std::string prefix = "eventcut_group"; - o2::framework::Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; - o2::framework::Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; - o2::framework::Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; - o2::framework::Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; - o2::framework::Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; - o2::framework::Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; - o2::framework::Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; - o2::framework::Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. - o2::framework::Configurable cfgRequireVertexTOFmatched{"cfgRequireVertexTOFmatched", false, "require Vertex TOFmatched in event cut"}; // ITS-TPC-TOF matched track contributes PV. - o2::framework::Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; - o2::framework::Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -2, "min. occupancy"}; - o2::framework::Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. occupancy"}; - o2::framework::Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -2, "min. FT0C occupancy"}; - o2::framework::Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; - o2::framework::Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; - o2::framework::Configurable cfgRequireNoCollInTimeRangeStrict{"cfgRequireNoCollInTimeRangeStrict", false, "require no collision in time range strict"}; - o2::framework::Configurable cfgRequireNoCollInITSROFStandard{"cfgRequireNoCollInITSROFStandard", false, "require no collision in time range standard"}; - o2::framework::Configurable cfgRequireNoCollInITSROFStrict{"cfgRequireNoCollInITSROFStrict", false, "require no collision in time range strict"}; - o2::framework::Configurable cfgRequireNoHighMultCollInPrevRof{"cfgRequireNoHighMultCollInPrevRof", false, "require no HM collision in previous ITS ROF"}; - o2::framework::Configurable cfgRequireGoodITSLayer3{"cfgRequireGoodITSLayer3", false, "number of inactive chips on ITS layer 3 are below threshold "}; - o2::framework::Configurable cfgRequireGoodITSLayer0123{"cfgRequireGoodITSLayer0123", false, "number of inactive chips on ITS layers 0-3 are below threshold "}; - o2::framework::Configurable cfgRequireGoodITSLayersAll{"cfgRequireGoodITSLayersAll", false, "number of inactive chips on all ITS layers are below threshold "}; - // for RCT - o2::framework::Configurable cfgRequireGoodRCT{"cfgRequireGoodRCT", false, "require good detector flag in run condtion table"}; - o2::framework::Configurable cfgRCTLabel{"cfgRCTLabel", "CBT_hadronPID", "select 1 [CBT, CBT_hadronPID, CBT_muon_glo] see O2Physics/Common/CCDB/RCTSelectionFlags.h"}; - o2::framework::Configurable cfgCheckZDC{"cfgCheckZDC", false, "set ZDC flag for PbPb"}; - o2::framework::Configurable cfgTreatLimitedAcceptanceAsBad{"cfgTreatLimitedAcceptanceAsBad", false, "reject all events where the detectors relevant for the specified Runlist are flagged as LimitedAcceptance"}; - - o2::framework::Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; - o2::framework::Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; - o2::framework::Configurable cfgNumContribMin{"cfgNumContribMin", 0, "min. numContrib"}; - o2::framework::Configurable cfgNumContribMax{"cfgNumContribMax", 65000, "max. numContrib"}; - } eventcuts; - - DielectronCut fDielectronCut; - struct : o2::framework::ConfigurableGroup { - std::string prefix = "dielectroncut_group"; - o2::framework::Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; - o2::framework::Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; - o2::framework::Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pT"}; - o2::framework::Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pT"}; - o2::framework::Configurable cfg_min_pair_y{"cfg_min_pair_y", -0.8, "min pair rapidity"}; - o2::framework::Configurable cfg_max_pair_y{"cfg_max_pair_y", +0.8, "max pair rapidity"}; - o2::framework::Configurable cfg_min_pair_dca3d{"cfg_min_pair_dca3d", 0.0, "min pair dca3d in sigma"}; - o2::framework::Configurable cfg_max_pair_dca3d{"cfg_max_pair_dca3d", 1e+10, "max pair dca3d in sigma"}; - o2::framework::Configurable cfg_apply_phiv{"cfg_apply_phiv", true, "flag to apply phiv cut"}; - o2::framework::Configurable cfg_phiv_slope{"cfg_phiv_slope", 0.0185, "slope for m vs. phiv"}; - o2::framework::Configurable cfg_phiv_intercept{"cfg_phiv_intercept", -0.0280, "intercept for m vs. phiv"}; - o2::framework::Configurable cfg_min_phiv{"cfg_min_phiv", 0.0, "min phiv (constant)"}; - o2::framework::Configurable cfg_max_phiv{"cfg_max_phiv", 3.2, "max phiv (constant)"}; - o2::framework::Configurable cfg_apply_detadphi{"cfg_apply_detadphi", false, "flag to apply deta-dphi elliptic cut at PV"}; - o2::framework::Configurable cfg_min_deta{"cfg_min_deta", 0.02, "min deta between 2 electrons (elliptic cut)"}; - o2::framework::Configurable cfg_min_dphi{"cfg_min_dphi", 0.2, "min dphi between 2 electrons (elliptic cut)"}; - o2::framework::Configurable cfg_min_opang{"cfg_min_opang", 0.0, "min opening angle"}; - o2::framework::Configurable cfg_max_opang{"cfg_max_opang", 6.4, "max opening angle"}; - // o2::framework::Configurable cfg_require_diff_sides{"cfg_require_diff_sides", false, "flag to require 2 tracks are from different sides."}; - - o2::framework::Configurable cfg_apply_cuts_from_prefilter{"cfg_apply_cuts_from_prefilter", false, "flag to apply prefilter set when producing derived data"}; - o2::framework::Configurable cfg_prefilter_bits{"cfg_prefilter_bits", 0, "prefilter bits [kNone : 0, kElFromPC : 1, kElFromPi0_20MeV : 2, kElFromPi0_40MeV : 4, kElFromPi0_60MeV : 8, kElFromPi0_80MeV : 16, kElFromPi0_100MeV : 32, kElFromPi0_120MeV : 64, kElFromPi0_140MeV : 128] Please consider logical-OR among them."}; // see PairUtilities.h - - o2::framework::Configurable cfg_apply_cuts_from_prefilter_derived{"cfg_apply_cuts_from_prefilter_derived", false, "flag to apply pair cut same as prefilter set in derived data"}; - o2::framework::Configurable cfg_prefilter_bits_derived{"cfg_prefilter_bits_derived", 0, "prefilter bits [kNone : 0, kMee : 1, kPhiV : 2, kSplitOrMergedTrackLS : 4, kSplitOrMergedTrackULS : 8] Please consider logical-OR among them."}; // see PairUtilities.h - - o2::framework::Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"}; - o2::framework::Configurable cfg_max_pt_track{"cfg_max_pt_track", 1e+10, "max pT for single track"}; - o2::framework::Configurable cfg_min_eta_track{"cfg_min_eta_track", -0.8, "min eta for single track"}; - o2::framework::Configurable cfg_max_eta_track{"cfg_max_eta_track", +0.8, "max eta for single track"}; - o2::framework::Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.f, "min phi for single track"}; - o2::framework::Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for single track"}; - o2::framework::Configurable cfg_mirror_phi_track{"cfg_mirror_phi_track", false, "mirror the phi cut around Pi, min and max Phi should be in 0-Pi"}; - o2::framework::Configurable cfg_reject_phi_track{"cfg_reject_phi_track", false, "reject the phi interval"}; - o2::framework::Configurable cfg_min_ncluster_tpc{"cfg_min_ncluster_tpc", 0, "min ncluster tpc"}; - o2::framework::Configurable cfg_min_ncluster_its{"cfg_min_ncluster_its", 5, "min ncluster its"}; - o2::framework::Configurable cfg_min_ncrossedrows{"cfg_min_ncrossedrows", 100, "min ncrossed rows"}; - o2::framework::Configurable cfg_max_frac_shared_clusters_tpc{"cfg_max_frac_shared_clusters_tpc", 999.f, "max fraction of shared clusters in TPC"}; - o2::framework::Configurable cfg_max_chi2tpc{"cfg_max_chi2tpc", 4.0, "max chi2/NclsTPC"}; - o2::framework::Configurable cfg_max_chi2its{"cfg_max_chi2its", 5.0, "max chi2/NclsITS"}; - o2::framework::Configurable cfg_max_chi2tof{"cfg_max_chi2tof", 1e+10, "max chi2 TOF"}; - o2::framework::Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 0.2, "max dca XY for single track in cm"}; - o2::framework::Configurable cfg_max_dcaz{"cfg_max_dcaz", 0.2, "max dca Z for single track in cm"}; - o2::framework::Configurable cfg_require_itsib_any{"cfg_require_itsib_any", false, "flag to require ITS ib any hits"}; - o2::framework::Configurable cfg_require_itsib_1st{"cfg_require_itsib_1st", true, "flag to require ITS ib 1st hit"}; - o2::framework::Configurable cfg_min_its_cluster_size{"cfg_min_its_cluster_size", 0.f, "min ITS cluster size"}; - o2::framework::Configurable cfg_max_its_cluster_size{"cfg_max_its_cluster_size", 16.f, "max ITS cluster size"}; - // o2::framework::Configurable cfg_min_rel_diff_pin{"cfg_min_rel_diff_pin", -1e+10, "min rel. diff. between pin and ppv"}; - // o2::framework::Configurable cfg_max_rel_diff_pin{"cfg_max_rel_diff_pin", +1e+10, "max rel. diff. between pin and ppv"}; - o2::framework::Configurable cfgRefR{"cfgRefR", 0.50, "ref. radius (m) for calculating phi position"}; // 0.50 +/- 0.06 can be syst. unc. - o2::framework::Configurable cfg_min_phiposition_track{"cfg_min_phiposition_track", 0.f, "min phi position for single track at certain radius"}; - o2::framework::Configurable cfg_max_phiposition_track{"cfg_max_phiposition_track", 6.3, "max phi position for single track at certain radius"}; - - o2::framework::Configurable cfg_pid_scheme{"cfg_pid_scheme", static_cast(DielectronCut::PIDSchemes::kTPChadrejORTOFreq), "pid scheme [kTOFreq : 0, kTPChadrej : 1, kTPChadrejORTOFreq : 2, kTPConly : 3, kTOFif : 4, kPIDML : 5, kTPChadrejORTOFreq_woTOFif : 6]"}; - o2::framework::Configurable cfg_min_TPCNsigmaEl{"cfg_min_TPCNsigmaEl", -2.0, "min. TPC n sigma for electron inclusion"}; - o2::framework::Configurable cfg_max_TPCNsigmaEl{"cfg_max_TPCNsigmaEl", +3.0, "max. TPC n sigma for electron inclusion"}; - // o2::framework::Configurable cfg_min_TPCNsigmaMu{"cfg_min_TPCNsigmaMu", -0.0, "min. TPC n sigma for muon exclusion"}; - // o2::framework::Configurable cfg_max_TPCNsigmaMu{"cfg_max_TPCNsigmaMu", +0.0, "max. TPC n sigma for muon exclusion"}; - o2::framework::Configurable cfg_min_TPCNsigmaPi{"cfg_min_TPCNsigmaPi", -1e+10, "min. TPC n sigma for pion exclusion"}; - o2::framework::Configurable cfg_max_TPCNsigmaPi{"cfg_max_TPCNsigmaPi", +3.0, "max. TPC n sigma for pion exclusion"}; - o2::framework::Configurable cfg_min_TPCNsigmaKa{"cfg_min_TPCNsigmaKa", -3.0, "min. TPC n sigma for kaon exclusion"}; - o2::framework::Configurable cfg_max_TPCNsigmaKa{"cfg_max_TPCNsigmaKa", +3.0, "max. TPC n sigma for kaon exclusion"}; - o2::framework::Configurable cfg_min_TPCNsigmaPr{"cfg_min_TPCNsigmaPr", -3.0, "min. TPC n sigma for proton exclusion"}; - o2::framework::Configurable cfg_max_TPCNsigmaPr{"cfg_max_TPCNsigmaPr", +3.0, "max. TPC n sigma for proton exclusion"}; - o2::framework::Configurable cfg_min_TOFNsigmaEl{"cfg_min_TOFNsigmaEl", -3.0, "min. TOF n sigma for electron inclusion"}; - o2::framework::Configurable cfg_max_TOFNsigmaEl{"cfg_max_TOFNsigmaEl", +3.0, "max. TOF n sigma for electron inclusion"}; - o2::framework::Configurable cfg_min_pin_pirejTPC{"cfg_min_pin_pirejTPC", 0.f, "min. pin for pion rejection in TPC"}; - o2::framework::Configurable cfg_max_pin_pirejTPC{"cfg_max_pin_pirejTPC", 1e+10, "max. pin for pion rejection in TPC"}; - o2::framework::Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; - - // configuration for PID ML - o2::framework::Configurable> onnxFileNames{"onnxFileNames", std::vector{"filename"}, "ONNX file names for each bin (if not from CCDB full path)"}; - o2::framework::Configurable> onnxPathsCCDB{"onnxPathsCCDB", std::vector{"path"}, "Paths of models on CCDB"}; - o2::framework::Configurable> binsMl{"binsMl", std::vector{0.1, 0.15, 0.2, 0.25, 0.4, 0.8, 1.6, 2.0, 20.f}, "Bin limits for ML application"}; - o2::framework::Configurable> cutsMl{"cutsMl", std::vector{0.98, 0.98, 0.9, 0.9, 0.95, 0.95, 0.8, 0.8}, "ML cuts per bin"}; - o2::framework::Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature"}, "Names of ML model input features"}; - o2::framework::Configurable nameBinningFeature{"nameBinningFeature", "pt", "Names of ML model binning feature"}; - o2::framework::Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB. Exceptions: > 0 for the specific timestamp, 0 gets the run dependent timestamp"}; - o2::framework::Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; - o2::framework::Configurable enableOptimizations{"enableOptimizations", false, "Enables the ONNX extended model-optimization: sessionOptions.SetGraphOptimizationLevel(GraphOptimizationLevel::ORT_ENABLE_EXTENDED)"}; - } dielectroncuts; - - DimuonCut fDimuonCut; - struct : o2::framework::ConfigurableGroup { - std::string prefix = "dimuoncut_group"; - o2::framework::Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; - o2::framework::Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; - o2::framework::Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pt"}; - o2::framework::Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pt"}; - o2::framework::Configurable cfg_min_pair_y{"cfg_min_pair_y", -4.0, "min pair rapidity"}; - o2::framework::Configurable cfg_max_pair_y{"cfg_max_pair_y", -2.5, "max pair rapidity"}; - o2::framework::Configurable cfg_min_pair_dcaxy{"cfg_min_pair_dcaxy", 0.0, "min pair dca3d in sigma"}; - o2::framework::Configurable cfg_max_pair_dcaxy{"cfg_max_pair_dcaxy", 1e+10, "max pair dca3d in sigma"}; - o2::framework::Configurable cfg_apply_detadphi{"cfg_apply_detadphi", false, "flag to apply deta-dphi elliptic cut"}; - o2::framework::Configurable cfg_min_deta{"cfg_min_deta", 0.02, "min deta between 2 muons (elliptic cut)"}; - o2::framework::Configurable cfg_min_dphi{"cfg_min_dphi", 0.02, "min dphi between 2 muons (elliptic cut)"}; - - o2::framework::Configurable cfg_apply_cuts_from_prefilter_derived{"cfg_apply_cuts_from_prefilter_derived", false, "flag to apply prefilter set in derived data"}; - o2::framework::Configurable cfg_prefilter_bits_derived{"cfg_prefilter_bits_derived", 0, "prefilter bits [kNone : 0, kSplitOrMergedTrackLS : 4, kSplitOrMergedTrackULS : 8] Please consider logical-OR among them."}; // see PairUtilities.h - - o2::framework::Configurable cfg_track_type{"cfg_track_type", 3, "muon track type [0: MFT-MCH-MID, 3: MCH-MID]"}; - o2::framework::Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"}; - o2::framework::Configurable cfg_max_pt_track{"cfg_max_pt_track", 1e+10, "max pT for single track"}; - o2::framework::Configurable cfg_min_eta_track{"cfg_min_eta_track", -4.0, "min eta for single track"}; - o2::framework::Configurable cfg_max_eta_track{"cfg_max_eta_track", -2.5, "max eta for single track"}; - o2::framework::Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.f, "min phi for single track"}; - o2::framework::Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for single track"}; - o2::framework::Configurable cfg_min_ncluster_mft{"cfg_min_ncluster_mft", 5, "min ncluster MFT"}; - o2::framework::Configurable cfg_min_ncluster_mch{"cfg_min_ncluster_mch", 5, "min ncluster MCH"}; - o2::framework::Configurable cfg_max_chi2{"cfg_max_chi2", 1e+6, "max chi2/ndf"}; - o2::framework::Configurable cfg_max_chi2mft{"cfg_max_chi2mft", 1e+6, "max chi2/ndf"}; - // o2::framework::Configurable cfg_max_matching_chi2_mftmch{"cfg_max_matching_chi2_mftmch", 40, "max chi2 for MFT-MCH matching"}; - o2::framework::Configurable cfg_border_pt_for_chi2mchmft{"cfg_border_pt_for_chi2mchmft", 0, "border pt for different max chi2 for MFT-MCH matching"}; - o2::framework::Configurable cfg_max_matching_chi2_mftmch_lowPt{"cfg_max_matching_chi2_mftmch_lowPt", 8, "max chi2 for MFT-MCH matching for low pT"}; - o2::framework::Configurable cfg_max_matching_chi2_mftmch_highPt{"cfg_max_matching_chi2_mftmch_highPt", 40, "max chi2 for MFT-MCH matching for high pT"}; - o2::framework::Configurable cfg_max_matching_chi2_mchmid{"cfg_max_matching_chi2_mchmid", 1e+10, "max chi2 for MCH-MID matching"}; - o2::framework::Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 1e+10, "max dca XY for single track in cm"}; - o2::framework::Configurable cfg_min_rabs{"cfg_min_rabs", 17.6, "min Radius at the absorber end"}; - o2::framework::Configurable cfg_max_rabs{"cfg_max_rabs", 89.5, "max Radius at the absorber end"}; - o2::framework::Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; - o2::framework::Configurable cfg_max_relDPt_wrt_matchedMCHMID{"cfg_max_relDPt_wrt_matchedMCHMID", 1e+10f, "max. relative dpt between MFT-MCH-MID and MCH-MID"}; - o2::framework::Configurable cfg_max_DEta_wrt_matchedMCHMID{"cfg_max_DEta_wrt_matchedMCHMID", 1e+10f, "max. deta between MFT-MCH-MID and MCH-MID"}; - o2::framework::Configurable cfg_max_DPhi_wrt_matchedMCHMID{"cfg_max_DPhi_wrt_matchedMCHMID", 1e+10f, "max. dphi between MFT-MCH-MID and MCH-MID"}; - o2::framework::Configurable requireMFTHitMap{"requireMFTHitMap", false, "flag to apply MFT hit map"}; - o2::framework::Configurable> requiredMFTDisks{"requiredMFTDisks", std::vector{0}, "hit map on MFT disks [0,1,2,3,4]. logical-OR of each double-sided disk"}; - } dimuoncuts; - - o2::aod::rctsel::RCTFlagsChecker rctChecker; - o2::framework::Service ccdb; - int mRunNumber; - float d_bz; - - o2::framework::HistogramRegistry fRegistry{"output", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, false}; - - float leptonM1 = 0.f; - float leptonM2 = 0.f; - - void init(o2::framework::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); - - DefineEMEventCut(); - 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 DefineEMEventCut() - { - fEMEventCut = EMEventCut("fEMEventCut", "fEMEventCut"); - fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); - fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); - fEMEventCut.SetZvtxRange(eventcuts.cfgZvtxMin, eventcuts.cfgZvtxMax); - fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); - fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); - fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); - fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); - fEMEventCut.SetRequireVertexTOFmatched(eventcuts.cfgRequireVertexTOFmatched); - fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); - fEMEventCut.SetRequireNoCollInTimeRangeStandard(eventcuts.cfgRequireNoCollInTimeRangeStandard); - fEMEventCut.SetRequireNoCollInTimeRangeStrict(eventcuts.cfgRequireNoCollInTimeRangeStrict); - fEMEventCut.SetRequireNoCollInITSROFStandard(eventcuts.cfgRequireNoCollInITSROFStandard); - fEMEventCut.SetRequireNoCollInITSROFStrict(eventcuts.cfgRequireNoCollInITSROFStrict); - fEMEventCut.SetRequireNoHighMultCollInPrevRof(eventcuts.cfgRequireNoHighMultCollInPrevRof); - fEMEventCut.SetRequireGoodITSLayer3(eventcuts.cfgRequireGoodITSLayer3); - fEMEventCut.SetRequireGoodITSLayer0123(eventcuts.cfgRequireGoodITSLayer0123); - fEMEventCut.SetRequireGoodITSLayersAll(eventcuts.cfgRequireGoodITSLayersAll); - } - - void DefineDielectronCut() - { - fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); - - // for pair - fDielectronCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); - fDielectronCut.SetPairPtRange(dielectroncuts.cfg_min_pair_pt, dielectroncuts.cfg_max_pair_pt); - fDielectronCut.SetPairYRange(dielectroncuts.cfg_min_pair_y, dielectroncuts.cfg_max_pair_y); - fDielectronCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma - fDielectronCut.SetMaxMeePhiVDep([&](float phiv) { return dielectroncuts.cfg_phiv_intercept + phiv * dielectroncuts.cfg_phiv_slope; }, dielectroncuts.cfg_min_phiv, dielectroncuts.cfg_max_phiv); - fDielectronCut.ApplyPhiV(dielectroncuts.cfg_apply_phiv); - fDielectronCut.SetMindEtadPhi(dielectroncuts.cfg_apply_detadphi, 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); - - // 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.SetMaxMatchingChi2MCHMFTPtDep([&](float pt) { return (pt < dimuoncuts.cfg_border_pt_for_chi2mchmft ? dimuoncuts.cfg_max_matching_chi2_mftmch_lowPt : dimuoncuts.cfg_max_matching_chi2_mftmch_highPt); }); - fDimuonCut.SetMatchingChi2MCHMID(0.f, dimuoncuts.cfg_max_matching_chi2_mchmid); - fDimuonCut.SetDCAxy(0.f, dimuoncuts.cfg_max_dcaxy); - fDimuonCut.SetRabs(dimuoncuts.cfg_min_rabs, dimuoncuts.cfg_max_rabs); - fDimuonCut.SetMaxPDCARabsDep([&](float rabs) { return (rabs < 26.5 ? 594.f : 324.f); }); - fDimuonCut.SetMaxdPtdEtadPhiwrtMCHMID(dimuoncuts.cfg_max_relDPt_wrt_matchedMCHMID, dimuoncuts.cfg_max_DEta_wrt_matchedMCHMID, dimuoncuts.cfg_max_DPhi_wrt_matchedMCHMID); // this is relevant for global muons - fDimuonCut.SetMFTHitMap(dimuoncuts.requireMFTHitMap, dimuoncuts.requiredMFTDisks); - } - - template - bool fillPairInfo(TCollision const&, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&) - { - if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { - if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { - return false; - } - } else { // cut-based - if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { - return false; - } - } - } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { - if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { - return false; - } - if (!map_best_match_globalmuon[t1.globalIndex()] || !map_best_match_globalmuon[t2.globalIndex()]) { - return false; - } - - // if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { - // return false; - // } - // if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { - // return false; - // } - } - - if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - if (!cut.IsSelectedPair(t1, t2)) { - return false; - } - } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { - if (!cut.IsSelectedPair(t1, t2)) { - return false; - } - } - - 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 = o2::aod::pwgem::dilepton::utils::emtrackutil::dca3DinSigma(t1); - dca2 = o2::aod::pwgem::dilepton::utils::emtrackutil::dca3DinSigma(t2); - if (cfgDCAType == 1) { - dca1 = o2::aod::pwgem::dilepton::utils::emtrackutil::dcaXYinSigma(t1); - dca2 = o2::aod::pwgem::dilepton::utils::emtrackutil::dcaXYinSigma(t2); - } else if (cfgDCAType == 2) { - dca1 = o2::aod::pwgem::dilepton::utils::emtrackutil::dcaZinSigma(t1); - dca2 = o2::aod::pwgem::dilepton::utils::emtrackutil::dcaZinSigma(t2); - } - } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { - dca1 = o2::aod::pwgem::dilepton::utils::emtrackutil::fwdDcaXYinSigma(t1); - dca2 = o2::aod::pwgem::dilepton::utils::emtrackutil::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; - } - - o2::framework::expressions::Filter collisionFilter_centrality = (eventcuts.cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < eventcuts.cfgCentMax) || (eventcuts.cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < eventcuts.cfgCentMax) || (eventcuts.cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < eventcuts.cfgCentMax); - o2::framework::expressions::Filter collisionFilter_numContrib = eventcuts.cfgNumContribMin <= o2::aod::collision::numContrib && o2::aod::collision::numContrib < eventcuts.cfgNumContribMax; - o2::framework::expressions::Filter collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; - o2::framework::expressions::Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; - using FilteredMyCollisions = o2::soa::Filtered; - - o2::framework::SliceCache cache; - o2::framework::Preslice perCollision_electron = o2::aod::emprimaryelectron::emeventId; - o2::framework::expressions::Filter trackFilter_electron = dielectroncuts.cfg_min_pt_track < o2::aod::track::pt && dielectroncuts.cfg_min_eta_track < o2::aod::track::eta && o2::aod::track::eta < dielectroncuts.cfg_max_eta_track && nabs(o2::aod::track::dcaXY) < dielectroncuts.cfg_max_dcaxy && nabs(o2::aod::track::dcaZ) < dielectroncuts.cfg_max_dcaz; - o2::framework::expressions::Filter pidFilter_electron = dielectroncuts.cfg_min_TPCNsigmaEl < o2::aod::pidtpc::tpcNSigmaEl && o2::aod::pidtpc::tpcNSigmaEl < dielectroncuts.cfg_max_TPCNsigmaEl; - o2::framework::expressions::Filter ttcaFilter_electron = ifnode(dielectroncuts.enableTTCA.node(), o2::aod::emprimaryelectron::isAssociatedToMPC == true || o2::aod::emprimaryelectron::isAssociatedToMPC == false, o2::aod::emprimaryelectron::isAssociatedToMPC == true); - o2::framework::expressions::Filter prefilter_derived_electron = ifnode(dielectroncuts.cfg_apply_cuts_from_prefilter_derived.node() && dielectroncuts.cfg_prefilter_bits_derived.node() >= static_cast(1), - ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kMee))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kMee))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kPhiV))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kPhiV))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) > static_cast(0), (o2::aod::emprimaryelectron::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) <= static_cast(0), true), - o2::aod::emprimaryelectron::pfbderived >= static_cast(0)); - - o2::framework::expressions::Filter prefilter_electron = ifnode(dielectroncuts.cfg_apply_cuts_from_prefilter.node() && dielectroncuts.cfg_prefilter_bits.node() >= static_cast(1), - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPC))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPC))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_20MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_20MeV))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_40MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_40MeV))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_60MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_60MeV))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_80MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_80MeV))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_100MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_100MeV))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_120MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_120MeV))) <= static_cast(0), true) && - ifnode((dielectroncuts.cfg_prefilter_bits.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_140MeV))) > static_cast(0), (o2::aod::emprimaryelectron::pfb & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBit::kElFromPi0_140MeV))) <= static_cast(0), true), - o2::aod::emprimaryelectron::pfb >= static_cast(0)); - - o2::framework::Partition positive_electrons = o2::aod::emprimaryelectron::sign > int8_t(0); - o2::framework::Partition negative_electrons = o2::aod::emprimaryelectron::sign < int8_t(0); - - o2::framework::Preslice perCollision_muon = o2::aod::emprimarymuon::emeventId; - o2::framework::expressions::Filter trackFilter_muon = o2::aod::fwdtrack::trackType == dimuoncuts.cfg_track_type && dimuoncuts.cfg_min_pt_track < o2::aod::fwdtrack::pt && o2::aod::fwdtrack::pt < dimuoncuts.cfg_max_pt_track && dimuoncuts.cfg_min_eta_track < o2::aod::fwdtrack::eta && o2::aod::fwdtrack::eta < dimuoncuts.cfg_max_eta_track; - o2::framework::expressions::Filter ttcaFilter_muon = ifnode(dimuoncuts.enableTTCA.node(), o2::aod::emprimarymuon::isAssociatedToMPC == true || o2::aod::emprimarymuon::isAssociatedToMPC == false, o2::aod::emprimarymuon::isAssociatedToMPC == true); - o2::framework::expressions::Filter prefilter_derived_muon = ifnode(dimuoncuts.cfg_apply_cuts_from_prefilter_derived.node() && dimuoncuts.cfg_prefilter_bits_derived.node() >= static_cast(1), - ifnode((dimuoncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) > static_cast(0), (o2::aod::emprimarymuon::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackLS))) <= static_cast(0), true) && - ifnode((dimuoncuts.cfg_prefilter_bits_derived.node() & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) > static_cast(0), (o2::aod::emprimarymuon::pfbderived & static_cast(1 << int(o2::aod::pwgem::dilepton::utils::pairutil::DileptonPrefilterBitDerived::kSplitOrMergedTrackULS))) <= static_cast(0), true), - o2::aod::emprimarymuon::pfbderived >= static_cast(0)); - o2::framework::Partition positive_muons = o2::aod::emprimarymuon::sign > int8_t(0); - o2::framework::Partition negative_muons = o2::aod::emprimarymuon::sign < int8_t(0); - - template - void runPairing(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut, TAllTracks const& tracks) - { - for (const auto& collision : collisions) { - initCCDB(collision); - const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; - // float centrality = centralities[cfgCentEstimator]; - if (centralities[cfgCentEstimator] < eventcuts.cfgCentMin || eventcuts.cfgCentMax < centralities[cfgCentEstimator]) { - continue; - } - - float eventplanes_2_for_mix[7] = {collision.ep2ft0m(), collision.ep2ft0a(), collision.ep2ft0c(), collision.ep2btot(), collision.ep2bpos(), collision.ep2bneg(), collision.ep2fv0a()}; - float ep2 = eventplanes_2_for_mix[cfgEP2Estimator_for_Mix]; - - 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; - - if (cfgStoreULS) { - for (const auto& [pos, neg] : combinations(o2::soa::CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS - bool is_pair_ok = fillPairInfo(collision, pos, neg, cut, tracks); - if (is_pair_ok) { - nuls++; - } - } - } - if (cfgStoreLS) { - for (const auto& [pos1, pos2] : combinations(o2::soa::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(o2::soa::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(), centralities[cfgCentEstimator], ep2); - } - } // end of collision loop - } // end of DF - - template - bool isPairOK(TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&) - { - if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { - if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { - return false; - } - } else { // cut-based - if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { - return false; - } - } - } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { - if (!cut.IsSelectedTrack(t1) || !cut.IsSelectedTrack(t2)) { - return false; - } - if (!map_best_match_globalmuon[t1.globalIndex()] || !map_best_match_globalmuon[t2.globalIndex()]) { - return false; - } - - // if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { - // return false; - // } - // if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { - // return false; - // } - } - - if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - if (!cut.IsSelectedPair(t1, t2)) { - return false; - } - } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { - if (!cut.IsSelectedPair(t1, t2)) { - return false; - } - } - return true; - } - - std::map, float> map_weight; // -> float - template - void fillPairWeightMap(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut, TAllTracks const& tracks) - { - std::vector> passed_pairIds; - passed_pairIds.reserve(posTracks.size() * negTracks.size()); - - for (const auto& collision : collisions) { - initCCDB(collision); - const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; - if (centralities[cfgCentEstimator] < eventcuts.cfgCentMin || eventcuts.cfgCentMax < centralities[cfgCentEstimator]) { - continue; - } - - if (!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(o2::soa::CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS - if (isPairOK(pos, neg, cut, tracks)) { - passed_pairIds.emplace_back(std::make_pair(pos.globalIndex(), neg.globalIndex())); - } - } - for (const auto& [pos1, pos2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ - if (isPairOK(pos1, pos2, cut, tracks)) { - passed_pairIds.emplace_back(std::make_pair(pos1.globalIndex(), pos2.globalIndex())); - } - } - for (const auto& [neg1, neg2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- - if (isPairOK(neg1, neg2, cut, tracks)) { - passed_pairIds.emplace_back(std::make_pair(neg1.globalIndex(), neg2.globalIndex())); - } - } - } // end of collision loop - - if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - for (const auto& pairId : passed_pairIds) { - auto t1 = tracks.rawIteratorAt(std::get<0>(pairId)); - auto t2 = tracks.rawIteratorAt(std::get<1>(pairId)); - // 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(); - } - - std::unordered_map map_best_match_globalmuon; - - void processAnalysis(FilteredMyCollisions const& collisions, Types const&... args) - { - if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - auto electrons = std::get<0>(std::tie(args...)); - if (cfgApplyWeightTTCA) { - fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); - } - runPairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); - } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { - auto muons = std::get<0>(std::tie(args...)); - map_best_match_globalmuon = o2::aod::pwgem::dilepton::utils::emtrackutil::findBestMatchMap(muons, fDimuonCut); - if (cfgApplyWeightTTCA) { - fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); - } - runPairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); - } - map_weight.clear(); - map_best_match_globalmuon.clear(); - } - PROCESS_SWITCH(DileptonProducer, processAnalysis, "run dilepton analysis", true); - - void processNorm(o2::aod::EMEventNormInfos const& collisions) - { - for (const auto& collision : collisions) { - if (collision.centFT0C() < eventcuts.cfgCentMin || eventcuts.cfgCentMax < collision.centFT0C()) { - continue; - } - - normTable(collision.selection_raw(), collision.rct_raw(), collision.posZint8(), collision.centFT0Muint8(), collision.centFT0Cuint8(), collision.centNTPVuint8() /*, collision.centNGlobaluint8()*/); - - } // end of collision loop - } - PROCESS_SWITCH(DileptonProducer, processNorm, "process normalization info", true); -}; - -#endif // PWGEM_DILEPTON_CORE_DILEPTONPRODUCER_H_ diff --git a/PWGEM/Dilepton/Core/DileptonSV.h b/PWGEM/Dilepton/Core/DileptonSV.h index 71d30ccf731..36576008f94 100644 --- a/PWGEM/Dilepton/Core/DileptonSV.h +++ b/PWGEM/Dilepton/Core/DileptonSV.h @@ -1108,15 +1108,15 @@ struct DileptonSV { candidate.phi2 = t2.phi(); } - // if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - // if (!cut.IsSelectedPair(t1, t2)) { - // return false; - // } - // } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { - // if (!cut.IsSelectedPair(t1, t2)) { - // return false; - // } - // } + if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { + if (!cut.IsSelectedPair(t1, t2)) { + return false; + } + } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { + if (!cut.IsSelectedPair(t1, t2)) { + return false; + } + } float weight = 1.f; if constexpr (ev_id == 0) { diff --git a/PWGEM/Dilepton/Core/DileptonSVMC.h b/PWGEM/Dilepton/Core/DileptonSVMC.h index dc766212e5d..b12f9602116 100644 --- a/PWGEM/Dilepton/Core/DileptonSVMC.h +++ b/PWGEM/Dilepton/Core/DileptonSVMC.h @@ -1573,9 +1573,9 @@ struct DileptonSVMC { return false; } - // if (!cut.IsSelectedPair(t1, t2)) { - // return false; - // } + if (!cut.IsSelectedPair(t1, 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; @@ -1617,9 +1617,9 @@ struct DileptonSVMC { return false; } - // if (!cut.IsSelectedPair(t1, t2)) { - // return false; - // } + if (!cut.IsSelectedPair(t1, t2)) { + return false; + } } float pt1 = 0.f, eta1 = 0.f, phi1 = 0.f, pt2 = 0.f, eta2 = 0.f, phi2 = 0.f; diff --git a/PWGEM/Dilepton/TableProducer/CMakeLists.txt b/PWGEM/Dilepton/TableProducer/CMakeLists.txt index cf28a3a8836..5d2889041ae 100644 --- a/PWGEM/Dilepton/TableProducer/CMakeLists.txt +++ b/PWGEM/Dilepton/TableProducer/CMakeLists.txt @@ -86,16 +86,6 @@ o2physics_add_dpl_workflow(filter-eoi 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) - o2physics_add_dpl_workflow(qvector-dummy-otf SOURCES qVectorDummyOTF.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore diff --git a/PWGEM/Dilepton/TableProducer/dielectronProducer.cxx b/PWGEM/Dilepton/TableProducer/dielectronProducer.cxx deleted file mode 100644 index 945ac2853c9..00000000000 --- a/PWGEM/Dilepton/TableProducer/dielectronProducer.cxx +++ /dev/null @@ -1,28 +0,0 @@ -// 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 "PWGEM/Dilepton/Utils/PairUtilities.h" - -#include -#include - -using namespace o2::framework; - -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 deleted file mode 100644 index ed298d39c39..00000000000 --- a/PWGEM/Dilepton/TableProducer/dimuonProducer.cxx +++ /dev/null @@ -1,28 +0,0 @@ -// 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 "PWGEM/Dilepton/Utils/PairUtilities.h" - -#include -#include - -using namespace o2::framework; - -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 bb0adb4d478..74c51456924 100644 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -160,11 +160,6 @@ o2physics_add_dpl_workflow(evaluate-acceptance 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) - o2physics_add_dpl_workflow(check-mc-template SOURCES checkMCTemplate.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::GlobalTracking diff --git a/PWGEM/Dilepton/Tasks/dileptonPolarization.cxx b/PWGEM/Dilepton/Tasks/dileptonPolarization.cxx deleted file mode 100644 index df4f9c14fda..00000000000 --- a/PWGEM/Dilepton/Tasks/dileptonPolarization.cxx +++ /dev/null @@ -1,650 +0,0 @@ -// 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/EMTrack.h" -#include "PWGEM/Dilepton/Utils/PairUtilities.h" - -#include "Common/Core/RecoDecay.h" -#include "Common/DataModel/EventSelection.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -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::pairutil; - -struct DileptonPolarization { - // Configurables - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - - Configurable cfgPairType{"cfgPairType", 0, "0:dielectron, 1:dimuon"}; - Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; - Configurable cfgDoMix{"cfgDoMix", true, "flag for event mixing"}; - // Configurable ndepth{"ndepth", 10000, "depth for event mixing"}; - Configurable ndiff_bc_mix{"ndiff_bc_mix", 594, "difference in global BC required in mixed events"}; - ConfigurableAxis ConfVtxBins{"ConfVtxBins", {10, -10.0f, 10.f}, "Mixing bins - z-vertex"}; - ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.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", {1, -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"}; - Configurable cfgUseAbs{"cfgUseAbs", false, "flag to use absolute value for cos_theta and phi"}; // this is to increase statistics per bin. - Configurable cfgDoULS{"cfgDoULS", true, "flag to perform ULS pairing"}; - Configurable cfgDoLS{"cfgDoLS", true, "flag to perform LS pairing"}; - - 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_pair_mass{"cfg_min_pair_mass", 0.0, "min pair mass"}; - Configurable cfg_max_pair_mass{"cfg_max_pair_mass", 1e+10, "max pair 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"; - Configurable cfgDM{"cfgDM", 0.01, "dm for lorentz boost"}; - Configurable cfgDPt{"cfgDPt", 0.1, "dpT for lorentz boost"}; - Configurable cfgDEta{"cfgDEta", 0.1, "deta for lorentz boost"}; - Configurable cfgDPhi{"cfgDPhi", 0.087, "dphi (rad.) for lorentz boost"}; - } accBins; - - 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/"}; - static constexpr std::string_view pair_sign_types[3] = {"uls/", "lspp/", "lsmm/"}; - - std::mt19937 engine; - std::vector zvtx_bin_edges; - std::vector cent_bin_edges; - std::vector ep_bin_edges; - std::vector occ_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; - } else { - LOG(fatal) << "Please select either dielectron or dimuon"; - } - - 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]); - } - } - - 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"; - // } - - 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); - - mRunNumber = collision.runNumber(); - } - - ~DileptonPolarization() {} - - 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/"); - fRegistry.add("Pair/same/uls/hEta", "#eta_{ll}", kTH1D, {{2000, -10, 10}}, true); - // fRegistry.add("Pair/mix/uls/hBeta", "#beta for Lorentz boost;#beta^{same};(#beta^{mix} - #beta^{same})/#beta^{same}", kTH2D, {{100, 0, 1}, {400, -0.2, 0.2}}, true); - } - - template - bool fillPairInfo(TCollision const& collision, TDilepton const& dilepton, TMixingBin const& mixingBin) - { - 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.M() < dileptoncuts.cfg_min_pair_mass || dileptoncuts.cfg_max_pair_mass < v12.M()) { - return false; - } - - if (v12.Pt() < dileptoncuts.cfg_min_pair_pt || dileptoncuts.cfg_max_pair_pt < v12.Pt()) { - return false; - } - - 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 (cfgUseAbs) { - cos_thetaPol = std::fabs(cos_thetaPol); - phiPol = std::fabs(phiPol); - } - - if (dilepton.sign1() * dilepton.sign2() < 0) { // ULS - fRegistry.fill(HIST("Pair/same/uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight); - fRegistry.fill(HIST("Pair/same/uls/hEta"), v12.Eta(), weight); - } else if (dilepton.sign1() > 0 && dilepton.sign2() > 0) { // LS++ - fRegistry.fill(HIST("Pair/same/lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight); - fRegistry.fill(HIST("Pair/same/lspp/hEta"), v12.Eta(), weight); - } else if (dilepton.sign1() < 0 && dilepton.sign2() < 0) { // LS-- - fRegistry.fill(HIST("Pair/same/lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight); - fRegistry.fill(HIST("Pair/same/lsmm/hEta"), v12.Eta(), weight); - } - - 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_pair(collision.globalIndex(), empair); - if (dilepton.sign1() * dilepton.sign2() < 0) { // ULS - mapMixingULS[mixingBin].emplace_back(pair_tmp); - } else if (dilepton.sign1() > 0 && dilepton.sign2() > 0) { // LS++ - mapMixingLSPP[mixingBin].emplace_back(pair_tmp); - } else if (dilepton.sign1() < 0 && dilepton.sign2() < 0) { // LS-- - mapMixingLSMM[mixingBin].emplace_back(pair_tmp); - } - return true; - } - - template - void fillMixedPairInfo(TPairs const& pairs) - { - float weight = 1.f; - - for (const auto& pair1 : pairs) { - auto globalBC1 = map_mixed_eventId_to_globalBC[std::get<0>(pair1)]; - auto empair1 = std::get<1>(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}; - - auto pairs_in_same_ptetaphim_bin = std::views::filter(pairs, [&](std::pair t) { return std::get<0>(t) != std::get<0>(pair1) && std::fabs(std::get<1>(t).mass() - empair1.mass()) < accBins.cfgDM && std::fabs(std::get<1>(t).pt() - empair1.pt()) < accBins.cfgDPt && std::fabs(std::get<1>(t).eta() - empair1.eta()) < accBins.cfgDEta && std::fabs(RecoDecay::constrainAngle(std::get<1>(t).phi() - empair1.phi(), -o2::constants::math::PI, 1U)) < accBins.cfgDPhi; }); - - for (const auto& pair2 : pairs_in_same_ptetaphim_bin) { - auto globalBC2 = map_mixed_eventId_to_globalBC[std::get<0>(pair2)]; - uint64_t diffBC = std::max(globalBC1, globalBC2) - std::min(globalBC1, globalBC2); - fRegistry.fill(HIST("Pair/mix/hDiffBC"), diffBC); - if (diffBC < ndiff_bc_mix) { - continue; - } - - auto empair2 = std::get<1>(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.f, 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; - if (cfgUseAbs) { - cos_thetaPol = std::fabs(cos_thetaPol); - phiPol = std::fabs(phiPol); - } - fRegistry.fill(HIST("Pair/mix/") + HIST(pair_sign_types[signType]) + HIST("hs"), empair1.mass(), empair1.pt(), empair1.getPairDCA(), empair1.rapidity(), cos_thetaPol, phiPol, quadmom, weight); - // fRegistry.fill(HIST("Pair/mix/") + HIST(pair_sign_types[signType]) + HIST("hBeta"), empair1.p() / empair1.e(), (empair2.p() / empair2.e() - empair1.p() / empair1.e()) / (empair1.p() / empair1.e())); - - } // end of pair2 loop - } // end of pair1 loop - } - - Filter collisionFilter_centrality = eventcuts.cfgCentMin < o2::aod::emthinevent::centrality && o2::aod::emthinevent::centrality < 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 > static_cast(0) && o2::aod::emdilepton::sign2 < static_cast(0)) || (o2::aod::emdilepton::sign1 < static_cast(0) && o2::aod::emdilepton::sign2 > static_cast(0)); - Partition dileptonsLSPP = o2::aod::emdilepton::sign1 > static_cast(0) && o2::aod::emdilepton::sign2 > static_cast(0); - Partition dileptonsLSMM = o2::aod::emdilepton::sign1 < static_cast(0) && o2::aod::emdilepton::sign2 < static_cast(0); - - std::map, std::vector>> mapMixingULS; // -> vector of for ULS pairs - std::map, std::vector>> mapMixingLSPP; // -> vector of for LSPP pairs - std::map, std::vector>> mapMixingLSMM; // -> vector of for LSMM pairs - - std::unordered_map map_mixed_eventId_to_globalBC; - - template - void runPairing(TCollisions const& collisions, TDileptons const&) - { - for (const auto& collision : collisions) { - initCCDB(collision); - float centrality = collision.centrality(); - 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); - fRegistry.fill(HIST("Event/after/hCorrOccupancy"), collision.ft0cOccupancyInTimeRange(), collision.trackOccupancyInTimeRange()); - fRegistry.fill(HIST("Event/after/hEP2_CentFT0C_forMix"), centrality, ep2); - - // 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; - } - - auto mixingBin = std::make_tuple(zbin, centbin, epbin, occbin); - - 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; - if (cfgDoULS) { - for (const auto& dilepton : dileptons_uls_per_coll) { // ULS - bool is_pair_ok = fillPairInfo(collision, dilepton, mixingBin); - if (is_pair_ok) { - nuls++; - } - } - } - - if (cfgDoLS) { - for (const auto& dilepton : dileptons_lspp_per_coll) { // LS++ - bool is_pair_ok = fillPairInfo(collision, dilepton, mixingBin); - if (is_pair_ok) { - nlspp++; - } - } - for (const auto& dilepton : dileptons_lsmm_per_coll) { // LS-- - bool is_pair_ok = fillPairInfo(collision, dilepton, mixingBin); - if (is_pair_ok) { - nlsmm++; - } - } - } - - if (!cfgDoMix || !(nuls > 0 || nlspp > 0 || nlsmm > 0)) { - continue; - } - - if (nuls > 0 || nlspp > 0 || nlsmm > 0) { - map_mixed_eventId_to_globalBC[collision.globalIndex()] = collision.globalBC(); - } // end of if pair exist - - } // end of collision loop - - for (int iz = 0; iz < static_cast(zvtx_bin_edges.size()) - 1; iz++) { - for (int icent = 0; icent < static_cast(cent_bin_edges.size()) - 1; icent++) { - for (int iep = 0; iep < static_cast(ep_bin_edges.size()) - 1; iep++) { - for (int iocc = 0; iocc < static_cast(occ_bin_edges.size()) - 1; iocc++) { - auto key_bin = std::make_tuple(iz, icent, iep, iocc); - - auto pairsULS = mapMixingULS[key_bin]; - auto pairsLSPP = mapMixingLSPP[key_bin]; - auto pairsLSMM = mapMixingLSMM[key_bin]; - LOGF(info, "iz = %d, icent = %d, iep = %d, iocc = %d, pairsULS.size() = %d, pairsLSPP.size() = %d, pairsLSMM.size() = %d", iz, icent, iep, iocc, pairsULS.size(), pairsLSPP.size(), pairsLSMM.size()); - - if (cfgDoULS) { - fillMixedPairInfo<0>(pairsULS); - } - if (cfgDoLS) { - fillMixedPairInfo<1>(pairsLSPP); - fillMixedPairInfo<2>(pairsLSMM); - } - - } // end of iocc loop - } // end of iep loop - } // end of icent loop - } // end of iz loop - - mapMixingULS.clear(); - mapMixingLSPP.clear(); - mapMixingLSMM.clear(); - map_mixed_eventId_to_globalBC.clear(); - } // end of DF - - void processAnalysis(filteredCollisions const& collisions, filteredDileptons const& dileptons) - { - runPairing(collisions, dileptons); - } - 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"})}; -}