diff --git a/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx b/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx index c8a26100cbd..e7c6b9c5b5d 100644 --- a/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx +++ b/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx @@ -27,7 +27,11 @@ #include "Common/DataModel/TrackSelectionTables.h" #include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/GeomConstants.h" #include "CommonConstants/MathConstants.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DetectorsBase/Propagator.h" +#include "Field/MagneticField.h" #include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" @@ -37,7 +41,11 @@ #include "Framework/runDataProcessing.h" #include "MathUtils/Utils.h" #include "ReconstructionDataFormats/GlobalTrackID.h" +#include "ReconstructionDataFormats/TrackFwd.h" +#include "Math/MatrixFunctions.h" +#include "Math/SMatrix.h" +#include "TGeoGlobalMagField.h" #include "TMCProcess.h" #include "TPDGCode.h" @@ -49,6 +57,9 @@ #include #include +using SMatrix55 = ROOT::Math::SMatrix>; +using SMatrix5 = ROOT::Math::SVector; + using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -59,9 +70,12 @@ using namespace o2::constants::math; using namespace pwgmm::mult; using namespace o2::aod::rctsel; -auto static constexpr kMinCharge = 3.f; -auto static constexpr kIntZero = 0; -auto static constexpr kZero = 0.f; +auto static constexpr CminCharge = 3.f; +auto static constexpr CintZero = 0; +auto static constexpr CintTwo = 2; +auto static constexpr Czero = 0.f; +auto static constexpr Cninety = 90.f; +auto static constexpr ConeHeighty = 180.f; enum TrkSel { trkSelAll, @@ -113,13 +127,23 @@ struct DndetaMFTPbPb { false, true}; - Configurable cfgDoIR{"cfgDoIR", false, "Flag to retrieve Interaction rate from CCDB"}; - Configurable cfgUseIRCut{"cfgUseIRCut", false, "Flag to cut on IR rate"}; - Configurable cfgIRCrashOnNull{"cfgIRCrashOnNull", false, "Flag to avoid CTP RateFetcher crash"}; - Configurable cfgIRSource{"cfgIRSource", "ZNC hadronic", "Estimator of the interaction rate (Pb-Pb: ZNC hadronic)"}; - Configurable cfgUseTrackSel{"cfgUseTrackSel", false, "Flag to apply track selection"}; - Configurable cfgUseParticleSel{"cfgUseParticleSel", false, "Flag to apply particle selection"}; - Configurable cfgRemoveReassigned{"cfgRemoveReassigned", false, "Remove reassgined tracks"}; + HistogramRegistry registryAlign{"registryAlign", {}}; + std::array mftQuadrant = {"Q0", "Q1", "Q2", "Q3"}; + std::array, 3>, 4>, 2> hAlignment; + + struct : ConfigurableGroup { + Configurable cfgDoIR{"cfgDoIR", false, "Flag to retrieve Interaction rate from CCDB"}; + Configurable cfgUseIRCut{"cfgUseIRCut", false, "Flag to cut on IR rate"}; + Configurable cfgIRCrashOnNull{"cfgIRCrashOnNull", false, "Flag to avoid CTP RateFetcher crash"}; + Configurable cfgIRSource{"cfgIRSource", "ZNC hadronic", "Estimator of the interaction rate (Pb-Pb: ZNC hadronic)"}; + Configurable cfgUseTrackSel{"cfgUseTrackSel", false, "Flag to apply track selection"}; + Configurable cfgUseParticleSel{"cfgUseParticleSel", false, "Flag to apply particle selection"}; + Configurable cfgRemoveReassigned{"cfgRemoveReassigned", false, "Remove reassgined tracks"}; + Configurable cfgUseTrackParExtra{"cfgUseTrackParExtra", false, "Use table with refitted track parameters"}; + Configurable cfgUseInelgt0{"cfgUseInelgt0", false, "Use INEL > 0 condition"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable cfgDCAtype{"cfgDCAtype", 2, "DCA coordinate type [0: DCA-X, 1: DCA-Y, 2: DCA-XY]"}; + } gConf; struct : ConfigurableGroup { ConfigurableAxis interactionRateBins{"interactionRateBins", {500, 0, 50}, "Binning for the interaction rate (kHz)"}; @@ -140,6 +164,8 @@ struct DndetaMFTPbPb { ConfigurableAxis etaBins{"etaBins", {20, -4., -2.}, "#eta binning"}; ConfigurableAxis chiSqPerNclBins{"chiSqPerNclBins", {100, 0, 100}, "#chi^{2} binning"}; ConfigurableAxis nClBins{"nClBins", {10, 0.5, 10.5}, "number of clusters binning"}; + ConfigurableAxis tanLambdaBins{"tanLambdaBins", {100, -25, 0}, "binning for tan(lambda)"}; + ConfigurableAxis invQPtBins{"invQPtBins", {200, -10, 10}, "binning for q/p_{T}"}; } binOpt; struct : ConfigurableGroup { @@ -171,13 +197,14 @@ struct DndetaMFTPbPb { Configurable minZvtx{"minZvtx", -20.0f, "minimum cut on z-vtx (cm)"}; Configurable useZDiffCut{"useZDiffCut", false, "use Zvtx reco-mc diff. cut"}; Configurable maxZvtxDiff{"maxZvtxDiff", 1.0f, "max allowed Z vtx difference for reconstruced collisions (cm)"}; + Configurable useZVtxCutMC{"useZVtxCutMC", false, "use Zvtx cut in MC"}; Configurable requireIsGoodZvtxFT0VsPV{"requireIsGoodZvtxFT0VsPV", true, "require events with PV position along z consistent (within 1 cm) between PV reconstructed using tracks and PV using FT0 A-C time difference"}; Configurable requireRejectSameBunchPileup{"requireRejectSameBunchPileup", true, "reject collisions in case of pileup with another collision in the same foundBC"}; - Configurable requireNoCollInTimeRangeStrict{"requireNoCollInTimeRangeStrict", true, " requireNoCollInTimeRangeStrict"}; - Configurable requireNoCollInRofStrict{"requireNoCollInRofStrict", false, "requireNoCollInRofStrict"}; + Configurable requireNoCollInTimeRangeStrict{"requireNoCollInTimeRangeStrict", false, " requireNoCollInTimeRangeStrict"}; + Configurable requireNoCollInRofStrict{"requireNoCollInRofStrict", true, "requireNoCollInRofStrict"}; Configurable requireNoCollInRofStandard{"requireNoCollInRofStandard", false, "requireNoCollInRofStandard"}; Configurable requireNoHighMultCollInPrevRof{"requireNoHighMultCollInPrevRof", false, "requireNoHighMultCollInPrevRof"}; - Configurable requireNoCollInTimeRangeStd{"requireNoCollInTimeRangeStd", true, "reject collisions corrupted by the cannibalism, with other collisions within +/- 10 microseconds"}; + Configurable requireNoCollInTimeRangeStd{"requireNoCollInTimeRangeStd", false, "reject collisions corrupted by the cannibalism, with other collisions within +/- 10 microseconds"}; Configurable requireNoCollInTimeRangeNarrow{"requireNoCollInTimeRangeNarrow", false, "reject collisions corrupted by the cannibalism, with other collisions within +/- 10 microseconds"}; Configurable occupancyEstimator{"occupancyEstimator", 1, "Occupancy estimator: 1 = trackOccupancyInTimeRange, 2 = ft0cOccupancyInTimeRange"}; Configurable minOccupancy{"minOccupancy", -1, "minimum occupancy from neighbouring collisions"}; @@ -199,6 +226,11 @@ struct DndetaMFTPbPb { TH2* gCurrentHadronicRate; RCTFlagsChecker rctChecker; + float bZ = 0; // Magnetic field for MFT + static constexpr double CcenterMFT[3] = {0, 0, -61.4}; // Field at center of MFT + + o2::parameters::GRPMagField* grpmag = nullptr; + std::vector ambiguousTrkIds; std::vector reassignedTrkIds; std::vector ambiguousTrkIdsMC; @@ -224,6 +256,8 @@ struct DndetaMFTPbPb { const AxisSpec etaAxis = {binOpt.etaBins, "#eta axis"}; const AxisSpec chiSqAxis = {binOpt.chiSqPerNclBins, "Chi2 axis"}; const AxisSpec nclsAxis = {binOpt.nClBins, "Number of clusters axis"}; + const AxisSpec tanLambdaAxis{binOpt.tanLambdaBins, "tan(#lambda)"}; + const AxisSpec invQPtAxis{binOpt.invQPtBins, "q/p_{T} (1/GeV)"}; rctChecker.init(rctCuts.cfgEvtRCTFlagCheckerLabel, rctCuts.cfgEvtRCTFlagCheckerZDCCheck, rctCuts.cfgEvtRCTFlagCheckerLimitAcceptAsBad); @@ -368,6 +402,11 @@ struct DndetaMFTPbPb { "; N_{ch}; occupancy", {HistType::kTH2F, {multAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/TanLambda", "; TanLambda; centrality; occupancy", {HistType::kTH2F, {tanLambdaAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/InvQPt", "; InvQPt; centrality; occupancy", {HistType::kTH2F, {invQPtAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Eta", "; #eta; centrality; occupancy", {HistType::kTH2F, {etaAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Phi", "; #varphi; centrality; occupancy", {HistType::kTH2F, {phiAxis, occupancyAxis}}}); + if (doprocessDatawBestTracksInclusive) { registry.add( {"Events/NtrkZvtxBest", @@ -416,6 +455,10 @@ struct DndetaMFTPbPb { qaregistry.add({"Tracks/TrackAmbDegree", "; N_{coll}^{comp}; occupancy", {HistType::kTH2F, {{51, -0.5, 50.5}, occupancyAxis}}}); + qaregistry.add({"Tracks/TanLambdaExtra", "; TanLambda; occupancy", {HistType::kTH2F, {tanLambdaAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/InvQPtExtra", "; InvQPt; occupancy", {HistType::kTH2F, {invQPtAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/EtaExtra", "; #eta; occupancy", {HistType::kTH2F, {etaAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/PhiExtra", "; #varphi; occupancy", {HistType::kTH2F, {phiAxis, occupancyAxis}}}); } } @@ -483,6 +526,11 @@ struct DndetaMFTPbPb { {HistType::kTHnSparseF, {nclsAxis, etaAxis, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/TanLambda", "; TanLambda; centrality; occupancy", {HistType::kTHnSparseF, {tanLambdaAxis, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/InvQPt", "; InvQPt; centrality; occupancy", {HistType::kTHnSparseF, {invQPtAxis, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/Eta", "; #eta; centrality; occupancy", {HistType::kTHnSparseF, {etaAxis, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/Phi", "; #varphi; centrality; occupancy", {HistType::kTHnSparseF, {phiAxis, centralityAxis, occupancyAxis}}}); + if (doprocessDatawBestTracksCentFT0C || doprocessDatawBestTracksCentFT0CVariant1 || doprocessDatawBestTracksCentFT0M || @@ -509,6 +557,10 @@ struct DndetaMFTPbPb { "; N_{coll}^{comp}; centrality; occupancy", {HistType::kTHnSparseF, {{51, -0.5, 50.5}, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/TanLambdaExtra", "; TanLambda; centrality; occupancy", {HistType::kTHnSparseF, {tanLambdaAxis, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/InvQPtExtra", "; InvQPt; centrality; occupancy", {HistType::kTHnSparseF, {invQPtAxis, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/EtaExtra", "; #eta; centrality; occupancy", {HistType::kTHnSparseF, {etaAxis, centralityAxis, occupancyAxis}}}); + qaregistry.add({"Tracks/Centrality/PhiExtra", "; #varphi; centrality; occupancy", {HistType::kTHnSparseF, {phiAxis, centralityAxis, occupancyAxis}}}); qaregistry.add( {"Tracks/Centrality/DCA3d", "; p_{T} (GeV/c); #eta; DCA_{XY} (cm); DCA_{Z} (cm); centrality; occupancy", @@ -870,13 +922,13 @@ struct DndetaMFTPbPb { if (doprocessReAssocMC) { // tracks not associated to any collision - hReAssoc[0] = qaregistry.add("ReAssoc/hNonAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hReAssoc[0] = qaregistry.add("ReAssoc/hNonAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); // tracks associasted to a collision - hReAssoc[1] = qaregistry.add("ReAssoc/hAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hReAssoc[1] = qaregistry.add("ReAssoc/hAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); // tracks associated to the correct collision considering only first reco collision (based on the MC collision index) - hReAssoc[2] = qaregistry.add("ReAssoc/hGoodAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hReAssoc[2] = qaregistry.add("ReAssoc/hGoodAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); // tracks associated to the correct collision considering all ambiguous reco collisions (based on the MC collision index) - hReAssoc[3] = qaregistry.add("ReAssoc/hGoodAssocTracksAmb", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hReAssoc[3] = qaregistry.add("ReAssoc/hGoodAssocTracksAmb", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); } if (doprocessEfficiencyInclusive) { @@ -972,6 +1024,12 @@ struct DndetaMFTPbPb { registry.add({"Tracks/THnDCAxyBestGenTruthPrim", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis}}}); registry.add({"Tracks/THnDCAxyBestGenPrimWrongColl", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis}}}); registry.add({"Tracks/THnDCAxyBestGenTruthPrimWrongColl", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis}}}); + + registry.add({"Tracks/BestGenPrimDeltaX", ";#Delta X (cm)", {HistType::kTH1F, {deltaZAxis}}}); + registry.add({"Tracks/BestGenSecDeltaX", ";#Delta X (cm)", {HistType::kTH1F, {deltaZAxis}}}); + registry.add({"Tracks/BestGenPrimWrongCollDeltaX", ";#Delta X (cm)", {HistType::kTH1F, {deltaZAxis}}}); + registry.add({"Tracks/BestGenSecWrongCollDeltaX", ";#Delta X (cm)", {HistType::kTH1F, {deltaZAxis}}}); + registry.add({"Tracks/THnDCAxyBestGenSec", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis}}}); registry.add({"Tracks/THnDCAxyBestGenTruthSec", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis}}}); registry.add({"Tracks/THnDCAxyBestGenSecWrongColl", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis}}}); @@ -991,6 +1049,12 @@ struct DndetaMFTPbPb { registry.add({"Tracks/Centrality/THnDCAxyBestGenTruthPrim", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm); centrality", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); registry.add({"Tracks/Centrality/THnDCAxyBestGenPrimWrongColl", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm); centrality", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); registry.add({"Tracks/Centrality/THnDCAxyBestGenTruthPrimWrongColl", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm); centrality", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + + registry.add({"Tracks/Centrality/BestGenPrimDeltaX", ";#Delta X (cm)", {HistType::kTH2F, {deltaZAxis, centralityAxis}}}); + registry.add({"Tracks/Centrality/BestGenSecDeltaX", ";#Delta X (cm)", {HistType::kTH2F, {deltaZAxis, centralityAxis}}}); + registry.add({"Tracks/Centrality/BestGenPrimWrongCollDeltaX", ";#Delta X (cm)", {HistType::kTH2F, {deltaZAxis, centralityAxis}}}); + registry.add({"Tracks/Centrality/BestGenSecWrongCollDeltaX", ";#Delta X (cm)", {HistType::kTH2F, {deltaZAxis, centralityAxis}}}); + registry.add({"Tracks/Centrality/THnDCAxyBestGenSec", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm); centrality", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); registry.add({"Tracks/Centrality/THnDCAxyBestGenTruthSec", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm); centrality", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); registry.add({"Tracks/Centrality/THnDCAxyBestGenSecWrongColl", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm); centrality", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); @@ -999,12 +1063,84 @@ struct DndetaMFTPbPb { registry.add({"Tracks/Centrality/THnDCAxyBestGenSecMat", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm); centrality", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); } } + + if (doprocessAlignmentInclusive) { + for (size_t j = 0; j < mftQuadrant.size(); j++) { + const auto& quadrant = mftQuadrant[j]; + std::string hPath = std::string("Alignment/") + quadrant + "/"; + hAlignment[0][j][0]["DCA_x_vs_z"] = registryAlign.add((hPath + "DCA_x_vs_z").c_str(), std::format("DCA(x) vs. z - {}", quadrant).c_str(), {HistType::kTH2F, {zAxis, dcaxyAxis}}); + hAlignment[0][j][0]["DCA_y_vs_z"] = registryAlign.add((hPath + "DCA_y_vs_z").c_str(), std::format("DCA(y) vs. z - {}", quadrant).c_str(), {HistType::kTH2F, {zAxis, dcaxyAxis}}); + } + } + + if (doprocessReassocEfficiencyCentFT0C) { + auto hNevt = registry.add("Events/hNReEffColls", "Number of generated and reconstructed MC collisions", HistType::kTH1F, {{3, 0.5, 3.5}}); + hNevt->GetXaxis()->SetBinLabel(1, "Reconstructed collisions"); + hNevt->GetXaxis()->SetBinLabel(2, "Generated collisions"); + + registry.add("Events/Centrality/ReEffStatus", ";status;centrality", HistType::kTH2F, {{8, 0.5, 8.5}, centralityAxis}); + auto hstat = registry.get(HIST("Events/Centrality/ReEffStatus")); + hstat->GetXaxis()->SetBinLabel(1, "All tracks"); + hstat->GetXaxis()->SetBinLabel(2, "Ambiguous tracks"); + hstat->GetXaxis()->SetBinLabel(3, "Reassigned tracks"); + hstat->GetXaxis()->SetBinLabel(4, "Extra tracks"); + hstat->GetXaxis()->SetBinLabel(5, "Originally correctly reassgined"); + hstat->GetXaxis()->SetBinLabel(6, "Correctly reassigned"); + hstat->GetXaxis()->SetBinLabel(7, "Not reassigned"); + hstat->GetXaxis()->SetBinLabel(8, "Correctly assigned true"); + + registry.add({"AmbTracks/DCAXY", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); + registry.add({"AmbTracks/DCAZ", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); + registry.add({"AmbTracks/DCAXYBest", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); + registry.add({"AmbTracks/DCAZBest", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); + registry.add({"AmbTracks/DCAXYBestPrim", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); + registry.add({"AmbTracks/DCAZBestPrim", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); + registry.add({"AmbTracks/DCAXYBestTrue", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); + registry.add({"AmbTracks/DCAZBestTrue", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); + registry.add({"AmbTracks/DCAXYBestFalse", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); + registry.add({"AmbTracks/DCAZBestFalse", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); + registry.add({"AmbTracks/DCAXYBestTrueOrigAssoc", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); + registry.add({"AmbTracks/DCAXYBestTrueOrigReAssoc", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); + registry.add({"AmbTracks/DCAZBestTrueOrigAssoc", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); + registry.add({"AmbTracks/DCAZBestTrueOrigReAssoc", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); + + registry.add({"AmbTracks/Centrality/THnDCAxyBestGen", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + registry.add({"AmbTracks/Centrality/THnDCAxyBestGenPrim", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + registry.add({"AmbTracks/Centrality/THnDCAxyBestTrue", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + registry.add({"AmbTracks/Centrality/THnDCAxyBestFalse", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + registry.add({"AmbTracks/Centrality/THnDCAxyBestTrueOrigAssoc", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + registry.add({"AmbTracks/Centrality/THnDCAxyBestTrueOrigReAssoc", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + } + + if (doprocessEventAndSignalLossCentFT0C) { + auto hNevt = registry.add("Events/Centrality/hNRecCollsSigEvtLoss", "Number of reconstructed MC collisions", HistType::kTH1F, {{1, 0.5, 1.5}}); + hNevt->GetXaxis()->SetBinLabel(1, "Reconstructed collisions"); + + registry.add("Events/Centrality/EvtSigLossStatus", ";status;centrality", {HistType::kTH2F, {{3, 0.5, 3.5}, centralityAxis}}); + auto hstat = registry.get(HIST("Events/Centrality/EvtSigLossStatus")); + hstat->GetXaxis()->SetBinLabel(1, "All MC gen events"); + hstat->GetXaxis()->SetBinLabel(2, "MC gen events with rec event with event selection"); + hstat->GetXaxis()->SetBinLabel(3, "MC gen events with no rec events"); + + registry.add("Events/Centrality/hMultGenVsCent", "event mult MC gen", {HistType::kTH2F, {centralityAxis, multFT0cAxis}}); + registry.add("Events/Centrality/hMultGenVsCentNParticlesEta05", "event mult MC gen", {HistType::kTH2F, {centralityAxis, multAxis}}); + registry.add("Events/Centrality/hMultGenVsCentNParticlesEtaMFT", "event mult MC gen", {HistType::kTH2F, {centralityAxis, multAxis}}); + registry.add("Events/Centrality/hMultGenVsCentRec", "event mult MC gen vs centrality", {HistType::kTH2F, {centralityAxis, multFT0cAxis}}); + registry.add("Events/Centrality/hMultGenVsCentRecNParticlesEta05", "event mult MC gen vs centrality", {HistType::kTH2F, {centralityAxis, multAxis}}); + registry.add("Events/Centrality/hMultGenVsCentRecNParticlesEtaMFT", "event mult MC gen vs centrality", {HistType::kTH2F, {centralityAxis, multAxis}}); + + registry.add({"Tracks/Centrality/EtaCentVsMultGen_t", "; #eta; centrality; mult gen", {HistType::kTHnSparseF, {etaAxis, centralityAxis, multFT0cAxis}}}); + registry.add({"Tracks/Centrality/EtaCentVsMultGen", "; #eta; centrality; mult gen", {HistType::kTHnSparseF, {etaAxis, centralityAxis, multFT0cAxis}}}); + + registry.add({"Tracks/Centrality/EtaGen_t", "; #eta; centrality; occupancy", {HistType::kTH2F, {etaAxis, centralityAxis}}}); + registry.add({"Tracks/Centrality/EtaGen", "; #eta; centrality; occupancy", {HistType::kTH2F, {etaAxis, centralityAxis}}}); + } } /// Filters - tracks Filter filtTrkEta = (aod::fwdtrack::eta < trackCuts.maxEta) && (aod::fwdtrack::eta > trackCuts.minEta); - Filter filtATrackID = (aod::fwdtrack::bestCollisionId >= kIntZero); + Filter filtATrackID = (aod::fwdtrack::bestCollisionId >= CintZero); Filter filtATrackDCAxy = (nabs(aod::fwdtrack::bestDCAXY) < trackCuts.maxDCAxy); Filter filtATrackDCAz = (nabs(aod::fwdtrack::bestDCAZ) < trackCuts.maxDCAz); @@ -1014,6 +1150,8 @@ struct DndetaMFTPbPb { /// Joined tables using FullBCs = soa::Join; using CollBCs = soa::Join; + using ExtBCs = soa::Join; + // Collisions using Colls = soa::Join; using Coll = Colls::iterator; @@ -1028,14 +1166,14 @@ struct DndetaMFTPbPb { using CollsGenCentFT0C = soa::Join; using CollisionsWithMCLabels = soa::Join; - using CollGenCent = CollsGenCentFT0C::iterator; using CollsCorr = soa::Join; + using CollsMCExtra = soa::Join; + // Tracks using MFTTracksLabeled = soa::Join; using MftTracksWColls = soa::Join; using MftTracksWCollsMC = soa::Join; - using BestTracksMC = soa::Join; /// Filtered tables @@ -1046,13 +1184,31 @@ struct DndetaMFTPbPb { using FiltParticles = soa::Filtered; + void initCCDB(ExtBCs::iterator const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + + grpmag = ccdb->getForTimeStamp(gConf.grpmagPath, bc.timestamp()); + LOG(info) << "Setting magnetic field to current " << grpmag->getL3Current() + << " A for run " << bc.runNumber() + << " from its GRPMagField CCDB object"; + o2::base::Propagator::initFieldFromGRP(grpmag); + mRunNumber = bc.runNumber(); + + o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + bZ = field->getBz(CcenterMFT); + LOG(info) << "The field at the center of the MFT is bZ = " << bZ; + } + template bool isBestTrackSelected(const B& besttrack) { if constexpr (fillHis) { registry.fill(HIST("Tracks/hBestTrkSel"), trkBestSelAll); } - if (besttrack.bestCollisionId() < kIntZero) { + if (besttrack.bestCollisionId() < CintZero) { return false; } if constexpr (fillHis) { @@ -1153,15 +1309,23 @@ struct DndetaMFTPbPb { if (fillHis) { float phi = track.phi(); o2::math_utils::bringTo02Pi(phi); - if (phi < kZero || TwoPI < phi) { + if (phi < Czero || TwoPI < phi) { continue; } if constexpr (has_reco_cent) { registry.fill(HIST("Tracks/Centrality/EtaZvtx"), track.eta(), z, c, occ); registry.fill(HIST("Tracks/Centrality/PhiEta"), phi, track.eta(), c, occ); + qaregistry.fill(HIST("Tracks/Centrality/TanLambda"), track.tgl(), c, occ); + qaregistry.fill(HIST("Tracks/Centrality/InvQPt"), track.signed1Pt(), c, occ); + qaregistry.fill(HIST("Tracks/Centrality/Eta"), track.eta(), c, occ); + qaregistry.fill(HIST("Tracks/Centrality/Phi"), phi, c, occ); } else { registry.fill(HIST("Tracks/EtaZvtx"), track.eta(), z, occ); registry.fill(HIST("Tracks/PhiEta"), phi, track.eta(), occ); + qaregistry.fill(HIST("Tracks/TanLambda"), track.tgl(), c, occ); + qaregistry.fill(HIST("Tracks/InvQPt"), track.signed1Pt(), c, occ); + qaregistry.fill(HIST("Tracks/Eta"), track.eta(), c, occ); + qaregistry.fill(HIST("Tracks/Phi"), phi, c, occ); } } ++nTrk; @@ -1176,6 +1340,35 @@ struct DndetaMFTPbPb { return nTrk; } + template + void countBestTracksExtra(B const& besttracksExtra, float c, float occ) + { + for (auto const& etrack : besttracksExtra) { + if (fillHis) { + float phi = etrack.phis(); + o2::math_utils::bringTo02Pi(phi); + if (phi < Czero || TwoPI < phi) { + continue; + } + if constexpr (has_reco_cent) { + if (gConf.cfgUseTrackParExtra) { + qaregistry.fill(HIST("Tracks/Centrality/TanLambdaExtra"), etrack.tgl(), c, occ); + qaregistry.fill(HIST("Tracks/Centrality/InvQPtExtra"), etrack.signed1Pt(), c, occ); + qaregistry.fill(HIST("Tracks/Centrality/EtaExtra"), etrack.etas(), c, occ); + qaregistry.fill(HIST("Tracks/Centrality/PhiExtra"), phi, c, occ); + } + } else { + if (gConf.cfgUseTrackParExtra) { + qaregistry.fill(HIST("Tracks/TanLambdaExtra"), etrack.tgl(), occ); + qaregistry.fill(HIST("Tracks/InvQPtExtra"), etrack.signed1Pt(), occ); + qaregistry.fill(HIST("Tracks/EtaExtra"), etrack.etas(), occ); + qaregistry.fill(HIST("Tracks/PhiExtra"), phi, occ); + } + } + } + } + } + template int countBestTracks(T const& tracks, B const& besttracks, float z, float c, float occ) @@ -1196,7 +1389,7 @@ struct DndetaMFTPbPb { if (fillHis) { float phi = itrack.phi(); o2::math_utils::bringTo02Pi(phi); - if (phi < kZero || TwoPI < phi) { + if (phi < Czero || TwoPI < phi) { continue; } if constexpr (has_reco_cent) { @@ -1220,7 +1413,7 @@ struct DndetaMFTPbPb { registry.fill(HIST("Tracks/hBestTrkSel"), trkBestSelNumReassoc); float phi = itrack.phi(); o2::math_utils::bringTo02Pi(phi); - if (phi < kZero || TwoPI < phi) { + if (phi < Czero || TwoPI < phi) { continue; } if constexpr (has_reco_cent) { @@ -1240,7 +1433,7 @@ struct DndetaMFTPbPb { } float phi = track.phi(); o2::math_utils::bringTo02Pi(phi); - if (phi < kZero || TwoPI < phi) { + if (phi < Czero || TwoPI < phi) { continue; } if (fillHis) { @@ -1297,7 +1490,7 @@ struct DndetaMFTPbPb { if (!isChrgParticle(particle.pdgCode())) { continue; } - if (cfgUseParticleSel && !isParticleSelected(particle)) { + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { continue; } if (particle.eta() < trackCuts.minEta || particle.eta() > trackCuts.maxEta) { @@ -1466,7 +1659,7 @@ struct DndetaMFTPbPb { if (p != nullptr) { charge = p->Charge(); } - return std::abs(charge) >= kMinCharge; + return std::abs(charge) >= CminCharge; } template @@ -1477,13 +1670,13 @@ struct DndetaMFTPbPb { if (!isChrgParticle(particle.pdgCode())) { continue; } - if (cfgUseParticleSel && !isParticleSelected(particle)) { + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { continue; } float phi = particle.phi(); o2::math_utils::bringTo02Pi(phi); - if (phi < kZero || TwoPI < phi) { + if (phi < Czero || TwoPI < phi) { continue; } if constexpr (isCent) { @@ -1499,7 +1692,7 @@ struct DndetaMFTPbPb { if (gtZeroColl) { float phi = particle.phi(); o2::math_utils::bringTo02Pi(phi); - if (phi < kZero || TwoPI < phi) { + if (phi < Czero || TwoPI < phi) { continue; } if constexpr (isCent) { @@ -1565,16 +1758,16 @@ struct DndetaMFTPbPb { return; } - if (cfgDoIR) { + if (gConf.cfgDoIR) { initHadronicRate(bc); - float ir = !cfgIRSource.value.empty() ? rateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), cfgIRSource, cfgIRCrashOnNull) * 1.e-3 : -1; + float ir = !gConf.cfgIRSource.value.empty() ? rateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), gConf.cfgIRSource, gConf.cfgIRCrashOnNull) * 1.e-3 : -1; if constexpr (has_reco_cent) { registry.fill(HIST("Events/Centrality/hInteractionRate"), c, occ, ir); } else { registry.fill(HIST("Events/hInteractionRate"), occ, ir); } float seconds = bc.timestamp() * 1.e-3 - mMinSeconds; - if (cfgUseIRCut && (ir < eventCuts.minIR || ir > eventCuts.maxIR)) { // cut on hadronic rate + if (gConf.cfgUseIRCut && (ir < eventCuts.minIR || ir > eventCuts.maxIR)) { // cut on hadronic rate return; } gCurrentHadronicRate->Fill(seconds, ir); @@ -1601,7 +1794,7 @@ struct DndetaMFTPbPb { template void processDatawBestTracks( typename C::iterator const& collision, FiltMftTracks const& tracks, - soa::SmallGroups const& besttracks, CollBCs const& /*bcs*/) + soa::SmallGroups const& besttracks, aod::BestCollisionsFwd3dExtra const& besttracksExtra, CollBCs const& /*bcs*/) { auto occ = getOccupancy(collision, eventCuts.occupancyEstimator); float c = getRecoCent(collision); @@ -1616,16 +1809,16 @@ struct DndetaMFTPbPb { return; } - if (cfgDoIR) { + if (gConf.cfgDoIR) { initHadronicRate(bc); - float ir = !cfgIRSource.value.empty() ? rateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), cfgIRSource, cfgIRCrashOnNull) * 1.e-3 : -1; + float ir = !gConf.cfgIRSource.value.empty() ? rateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), gConf.cfgIRSource, gConf.cfgIRCrashOnNull) * 1.e-3 : -1; if constexpr (has_reco_cent) { registry.fill(HIST("Events/Centrality/hInteractionRate"), c, occ, ir); } else { registry.fill(HIST("Events/hInteractionRate"), occ, ir); } float seconds = bc.timestamp() * 1.e-3 - mMinSeconds; - if (cfgUseIRCut && (ir < eventCuts.minIR || ir > eventCuts.maxIR)) { // cut on hadronic rate + if (gConf.cfgUseIRCut && (ir < eventCuts.minIR || ir > eventCuts.maxIR)) { // cut on hadronic rate return; } gCurrentHadronicRate->Fill(seconds, ir); @@ -1641,6 +1834,7 @@ struct DndetaMFTPbPb { } auto nBestTrks = countBestTracks(tracks, besttracks, z, c, occ); + countBestTracksExtra(besttracksExtra, c, occ); if constexpr (has_reco_cent) { registry.fill(HIST("Events/Centrality/NtrkZvtxBest"), nBestTrks, z, c, @@ -1707,9 +1901,9 @@ struct DndetaMFTPbPb { void processDatawBestTracksInclusive( Colls::iterator const& collision, FiltMftTracks const& tracks, - soa::SmallGroups const& besttracks, CollBCs const& bcs) + soa::SmallGroups const& besttracks, aod::BestCollisionsFwd3dExtra const& besttracksExtra, CollBCs const& bcs) { - processDatawBestTracks(collision, tracks, besttracks, bcs); + processDatawBestTracks(collision, tracks, besttracks, besttracksExtra, bcs); } PROCESS_SWITCH(DndetaMFTPbPb, processDatawBestTracksInclusive, @@ -1718,9 +1912,9 @@ struct DndetaMFTPbPb { void processDatawBestTracksCentFT0C( CollsCentFT0C::iterator const& collision, FiltMftTracks const& tracks, - soa::SmallGroups const& besttracks, CollBCs const& bcs) + soa::SmallGroups const& besttracks, aod::BestCollisionsFwd3dExtra const& besttracksExtra, CollBCs const& bcs) { - processDatawBestTracks(collision, tracks, besttracks, bcs); + processDatawBestTracks(collision, tracks, besttracks, besttracksExtra, bcs); } PROCESS_SWITCH(DndetaMFTPbPb, processDatawBestTracksCentFT0C, @@ -1730,9 +1924,9 @@ struct DndetaMFTPbPb { void processDatawBestTracksCentFT0CVariant1( CollsCentFT0CVariant1::iterator const& collision, FiltMftTracks const& tracks, - soa::SmallGroups const& besttracks, CollBCs const& bcs) + soa::SmallGroups const& besttracks, aod::BestCollisionsFwd3dExtra const& besttracksExtra, CollBCs const& bcs) { - processDatawBestTracks(collision, tracks, besttracks, bcs); + processDatawBestTracks(collision, tracks, besttracks, besttracksExtra, bcs); } PROCESS_SWITCH(DndetaMFTPbPb, processDatawBestTracksCentFT0CVariant1, @@ -1742,9 +1936,9 @@ struct DndetaMFTPbPb { void processDatawBestTracksCentFT0M( CollsCentFT0M::iterator const& collision, FiltMftTracks const& tracks, - soa::SmallGroups const& besttracks, CollBCs const& bcs) + soa::SmallGroups const& besttracks, aod::BestCollisionsFwd3dExtra const& besttracksExtra, CollBCs const& bcs) { - processDatawBestTracks(collision, tracks, besttracks, bcs); + processDatawBestTracks(collision, tracks, besttracks, besttracksExtra, bcs); } PROCESS_SWITCH(DndetaMFTPbPb, processDatawBestTracksCentFT0M, @@ -1753,9 +1947,9 @@ struct DndetaMFTPbPb { void processDatawBestTracksCentNGlobal( CollsCentNGlobal::iterator const& collision, FiltMftTracks const& tracks, - soa::SmallGroups const& besttracks, CollBCs const& bcs) + soa::SmallGroups const& besttracks, aod::BestCollisionsFwd3dExtra const& besttracksExtra, CollBCs const& bcs) { - processDatawBestTracks(collision, tracks, besttracks, bcs); + processDatawBestTracks(collision, tracks, besttracks, besttracksExtra, bcs); } PROCESS_SWITCH(DndetaMFTPbPb, processDatawBestTracksCentNGlobal, @@ -1765,9 +1959,9 @@ struct DndetaMFTPbPb { void processDatawBestTracksCentMFT( CollsCentMFT::iterator const& collision, FiltMftTracks const& tracks, - soa::SmallGroups const& besttracks, CollBCs const& bcs) + soa::SmallGroups const& besttracks, aod::BestCollisionsFwd3dExtra const& besttracksExtra, CollBCs const& bcs) { - processDatawBestTracks(collision, tracks, besttracks, bcs); + processDatawBestTracks(collision, tracks, besttracks, besttracksExtra, bcs); } PROCESS_SWITCH(DndetaMFTPbPb, processDatawBestTracksCentMFT, @@ -2203,7 +2397,7 @@ struct DndetaMFTPbPb { if (!isChrgParticle(particle.pdgCode())) { continue; } - if (cfgUseParticleSel && !isParticleSelected(particle)) { + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { continue; } @@ -2348,7 +2542,7 @@ struct DndetaMFTPbPb { if (!isChrgParticle(particle.pdgCode())) { continue; } - if (cfgUseParticleSel && !isParticleSelected(particle)) { + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { continue; } @@ -2478,7 +2672,7 @@ struct DndetaMFTPbPb { if (!isChrgParticle(particle.pdgCode())) { continue; } - if (cfgUseParticleSel && !isParticleSelected(particle)) { + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { continue; } if constexpr (has_reco_cent) { @@ -2799,7 +2993,7 @@ struct DndetaMFTPbPb { } auto itrack = track.mfttrack_as(); - if (cfgRemoveReassigned) { + if (gConf.cfgRemoveReassigned) { if (itrack.collisionId() != track.bestCollisionId()) { continue; } @@ -2810,20 +3004,24 @@ struct DndetaMFTPbPb { if (itrack.has_mcParticle()) { auto particle = itrack.mcParticle_as(); + float deltaX = -999.f; + float deltaY = -999.f; float deltaZ = -999.f; if (index) { auto collision = itrack.collision_as(); auto mcCollision = particle.mcCollision_as(); + deltaX = collision.posX() - mcCollision.posX(); + deltaY = collision.posY() - mcCollision.posY(); deltaZ = collision.posZ() - mcCollision.posZ(); if (collision.has_mcCollision() && collision.mcCollisionId() == particle.mcCollisionId()) { - hReAssoc[index + 1]->Fill(itrack.pt(), itrack.eta(), deltaZ); + hReAssoc[index + 1]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); } else { - hReAssoc[index + 2]->Fill(itrack.pt(), itrack.eta(), deltaZ); + hReAssoc[index + 2]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); } } - hReAssoc[index]->Fill(itrack.pt(), itrack.eta(), deltaZ); + hReAssoc[index]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); } else { - hReAssoc[index]->Fill(itrack.pt(), itrack.eta(), -999.f); + hReAssoc[index]->Fill(itrack.pt(), itrack.eta(), -999.f, -999.f, -999.f); } } } @@ -3026,7 +3224,7 @@ struct DndetaMFTPbPb { if (particle.eta() <= trackCuts.minEta || particle.eta() >= trackCuts.maxEta) { continue; } - if (cfgUseParticleSel && !isParticleSelected(particle)) { + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { continue; } @@ -3203,14 +3401,14 @@ struct DndetaMFTPbPb { } float phi = itrack.phi(); o2::math_utils::bringTo02Pi(phi); - if (phi < kZero || TwoPI < phi) { + if (phi < Czero || TwoPI < phi) { continue; } if (!itrack.has_collision()) { continue; } - if (cfgRemoveReassigned) { + if (gConf.cfgRemoveReassigned) { if (itrack.collisionId() != atrack.bestCollisionId()) { continue; } @@ -3230,7 +3428,7 @@ struct DndetaMFTPbPb { if (particle.eta() <= trackCuts.minEta || particle.eta() >= trackCuts.maxEta) { continue; } - if (cfgUseParticleSel && !isParticleSelected(particle)) { + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { continue; } @@ -3249,9 +3447,11 @@ struct DndetaMFTPbPb { if (collision.has_mcCollision() && collision.mcCollisionId() == particle.mcCollisionId()) { if (!particle.isPhysicalPrimary()) { // Secondaries (weak decays and material) if constexpr (has_reco_cent) { + registry.fill(HIST("Tracks/Centrality/BestGenSecDeltaX"), collision.posZ() - particle.mcCollision().posZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenSec"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenTruthSec"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth, crec); } else { + registry.fill(HIST("Tracks/BestGenSecDeltaX"), collision.posZ() - particle.mcCollision().posZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenSec"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenTruthSec"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth); } @@ -3270,9 +3470,11 @@ struct DndetaMFTPbPb { } } else { // Primaries if constexpr (has_reco_cent) { + registry.fill(HIST("Tracks/Centrality/BestGenPrimDeltaX"), collision.posZ() - particle.mcCollision().posZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenPrim"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenTruthPrim"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth, crec); } else { + registry.fill(HIST("Tracks/BestGenPrimDeltaX"), collision.posZ() - particle.mcCollision().posZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenPrim"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenTruthPrim"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth); } @@ -3280,17 +3482,21 @@ struct DndetaMFTPbPb { } else { // Wrong collision if (!particle.isPhysicalPrimary()) { // Secondaries (weak decays and material) if constexpr (has_reco_cent) { + registry.fill(HIST("Tracks/Centrality/BestGenSecWrongCollDeltaX"), collision.posZ() - particle.mcCollision().posZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenSecWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenTruthSecWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth, crec); } else { + registry.fill(HIST("Tracks/BestGenSecWrongCollDeltaX"), collision.posZ() - particle.mcCollision().posZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenSecWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenTruthSecWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth); } } else { // Primaries if constexpr (has_reco_cent) { + registry.fill(HIST("Tracks/Centrality/BestGenPrimWrongCollDeltaX"), collision.posZ() - particle.mcCollision().posZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenPrimWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ(), crec); registry.fill(HIST("Tracks/Centrality/THnDCAxyBestGenTruthPrimWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth, crec); } else { + registry.fill(HIST("Tracks/BestGenPrimWrongCollDeltaX"), collision.posZ() - particle.mcCollision().posZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenPrimWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), atrack.bestDCAXY(), atrack.bestDCAZ()); registry.fill(HIST("Tracks/THnDCAxyBestGenTruthPrimWrongColl"), particle.pt(), particle.eta(), mcCollision.posZ(), dcaXYtruth, dcaZtruth); } @@ -3339,14 +3545,14 @@ struct DndetaMFTPbPb { auto nBestTrks = 0; for (auto const& atrack : besttracks) { - if (cfgUseTrackSel && !isBestTrackSelected(atrack)) { + if (gConf.cfgUseTrackSel && !isBestTrackSelected(atrack)) { continue; } auto itrack = atrack.template mfttrack_as(); if (itrack.eta() < trackCuts.minEta || itrack.eta() > trackCuts.maxEta) { continue; } - if (cfgUseTrackSel && !isTrackSelected(itrack)) { + if (gConf.cfgUseTrackSel && !isTrackSelected(itrack)) { continue; } nBestTrks++; @@ -3364,6 +3570,351 @@ struct DndetaMFTPbPb { } PROCESS_SWITCH(DndetaMFTPbPb, processCorrelationwBestTracksInclusive, "Do correlation study based on BestCollisionsFwd3d table", false); + + int getQuadrantPhi(float phi) + { + if (phi >= Czero && phi < Cninety) { + return 0; + } + if (phi >= Cninety && phi <= ConeHeighty) { + return 1; + } + if (phi >= -ConeHeighty && phi < -Cninety) { + return 2; + } + if (phi >= -Cninety && phi < Czero) { + return 3; + } + return -1; + } + + template + int getQuadrantTrack(T const& track) + { + float phi = static_cast(track.phi()) * ConeHeighty / PI; + return getQuadrantPhi(phi); + } + + /// @brief process function to check MFT alignment (based on BestCollisionsFwd3d table) + template + void processAlignment(typename C::iterator const& collision, + FiltMftTracks const& /*tracks*/, + soa::SmallGroups const& besttracks, + CollBCs const& /*bcs*/ + ) + { + auto bc = collision.template foundBC_as(); + if (!isGoodEvent(collision)) { + return; + } + + if (gConf.cfgDoIR) { + initHadronicRate(bc); + float ir = !gConf.cfgIRSource.value.empty() ? rateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), gConf.cfgIRSource, gConf.cfgIRCrashOnNull) * 1.e-3 : -1; + float seconds = bc.timestamp() * 1.e-3 - mMinSeconds; + if (gConf.cfgUseIRCut && (ir < eventCuts.minIR || ir > eventCuts.maxIR)) { // cut on hadronic rate + return; + } + gCurrentHadronicRate->Fill(seconds, ir); + } + + auto z = collision.posZ(); + for (auto const& atrack : besttracks) { + if (!isBestTrackSelected(atrack)) { + continue; + } + auto itrack = atrack.template mfttrack_as(); + if (!isTrackSelected(itrack)) { + continue; + } + + int quadrant = getQuadrantTrack(itrack); + if (quadrant < 0) { + continue; + } + std::get>(hAlignment[0][quadrant][0]["DCA_x_vs_z"])->Fill(z, atrack.bestDCAXY()); + } + } + + void processAlignmentInclusive(Colls::iterator const& collision, + FiltMftTracks const& tracks, + soa::SmallGroups const& besttracks, + CollBCs const& bcs) + { + processAlignment(collision, tracks, besttracks, bcs); + } + + PROCESS_SWITCH(DndetaMFTPbPb, processAlignmentInclusive, "Check MFT alignment using tracks based on BestCollisionsFwd3d table (inclusive)", false); + + Preslice compCollPerCol = o2::aod::fwdtrack::collisionId; + + /// @brief process function to check reassociation efficiency + template + void processReassocEfficiency(typename soa::Join const& collisions, + MftTracksWCollsMC const& tracks, + MC const& mcCollisions, + aod::McParticles const& /*particles*/, + ExtBCs const& bcs) + { + if (bcs.size() == 0) { + return; + } + auto bc = bcs.begin(); + initCCDB(bc); + + registry.fill(HIST("Events/hNReEffColls"), 1.f, collisions.size()); + registry.fill(HIST("Events/hNReEffColls"), 2.f, mcCollisions.size()); + + for (const auto& collision : collisions) { + if (!isGoodEvent(collision)) { + continue; + } + float crec = getRecoCent(collision); + registry.fill(HIST("Events/Centrality/ReEffStatus"), 1, crec); + + if (!collision.has_mcCollision()) { + continue; + } + + std::array dcaInfOrig; + std::array dcaInfo; + double bestDCA[2]; + auto trkPerColl = tracks.sliceBy(compCollPerCol, collision.globalIndex()); + + for (auto const& track : trkPerColl) { + dcaInfOrig[0] = 999.f; // original DCAx from propagation + dcaInfOrig[1] = 999.f; // original DCAy from propagation + dcaInfOrig[2] = 999.f; // original DCAz from propagation + dcaInfo[0] = 999.f; // calcualted DCAxy + dcaInfo[1] = 999.f; // calculated DCAz - same as original + bestDCA[0] = 999.f; // minimal DCAxy + bestDCA[1] = 999.f; // minimal DCAz + + if (!isTrackSelected(track)) { + continue; + } + + auto bestCol = track.collisionId(); + uint index = uint(track.collisionId() >= 0); + auto ids = track.compatibleCollIds(); + bool isAmbiguous = (ids.size() > 1); + + if (!track.has_mcParticle()) { + continue; + } + auto particle = track.template mcParticle_as(); + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (particle.eta() <= trackCuts.minEta || particle.eta() >= trackCuts.maxEta) { + continue; + } + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { + continue; + } + + int bestMCCol = -1; + std::vector v1; // Temporary null vector for the computation of the covariance matrix + SMatrix55 tcovs(v1.begin(), v1.end()); + SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + o2::track::TrackParCovFwd trackPar{track.z(), tpars, tcovs, track.chi2()}; + + if (index) { + auto mcCollision = particle.template mcCollision_as(); + if (eventCuts.useZDiffCut) { + if (std::abs(collision.posZ() - mcCollision.posZ()) > eventCuts.maxZvtxDiff) { + continue; + } + } + + if (isAmbiguous) { + for (const auto& collIdx : track.compatibleCollIds()) { + auto ambCollision = collisions.rawIteratorAt(collIdx); + if (!ambCollision.has_mcCollision()) { + continue; + } + + trackPar.propagateToDCAhelix(bZ, {ambCollision.posX(), ambCollision.posY(), ambCollision.posZ()}, dcaInfOrig); + + if (gConf.cfgDCAtype == 0) { + dcaInfo[0] = dcaInfOrig[0]; + } else if (gConf.cfgDCAtype == 1) { + dcaInfo[0] = dcaInfOrig[1]; + } else if (gConf.cfgDCAtype == CintTwo) { + dcaInfo[0] = std::sqrt(dcaInfOrig[0] * dcaInfOrig[0] + dcaInfOrig[1] * dcaInfOrig[1]); + } + dcaInfo[1] = dcaInfOrig[2]; + + // LOGP(info, "-> track {}: {}", track.globalIndex(), dcaInfo[0]); + registry.fill(HIST("AmbTracks/DCAXY"), dcaInfo[0]); + registry.fill(HIST("AmbTracks/DCAZ"), dcaInfo[1]); + + if ((std::abs(dcaInfo[0]) < std::abs(bestDCA[0])) && (std::abs(dcaInfo[1]) < std::abs(bestDCA[1]))) { + bestCol = ambCollision.globalIndex(); + bestMCCol = ambCollision.mcCollisionId(); + bestDCA[0] = dcaInfo[0]; + bestDCA[1] = dcaInfo[1]; + } + } + + registry.fill(HIST("AmbTracks/DCAXYBest"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBest"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestGen"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + + if (particle.isPhysicalPrimary()) { + registry.fill(HIST("AmbTracks/DCAXYBestPrim"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestPrim"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestGenPrim"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + } + + auto mcCollID = particle.mcCollisionId(); + registry.fill(HIST("Events/Centrality/ReEffStatus"), 2, crec); + + if (!track.has_collision()) { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 4, crec); + if (bestMCCol == mcCollID) // correctly assigned to bestCol + { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 6, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrue"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrue"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrue"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + } else { + registry.fill(HIST("AmbTracks/DCAXYBestFalse"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestFalse"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestFalse"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + } + } else if (track.collisionId() != bestCol) { // reassgined + auto collOrig = track.template collision_as>(); + registry.fill(HIST("Events/Centrality/ReEffStatus"), 3, crec); + if (bestMCCol == mcCollID) // correctly reassigned + { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 6, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrue"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrue"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrue"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + } else { // uncorrectly reassigned + registry.fill(HIST("AmbTracks/DCAXYBestFalse"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestFalse"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestFalse"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + } + if (collOrig.mcCollisionId() == mcCollID) { // initially correctly assigned + registry.fill(HIST("Events/Centrality/ReEffStatus"), 5, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrueOrigReAssoc"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrueOrigReAssoc"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrueOrigReAssoc"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + } + } else { // the track has a collision and track.collisionId() == bestCol + if (track.collisionId() != bestCol) { + // LOGP(info, "-> track id {}: bestCollid {}", track.collisionId(), bestCol); + } + registry.fill(HIST("Events/Centrality/ReEffStatus"), 7, crec); + if (bestMCCol == mcCollID) // correctly assigned + { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 8, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrueOrigAssoc"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrueOrigAssoc"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrueOrigAssoc"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + } + } + } + } + } + } + } + + void processReassocEfficiencyCentFT0C(soa::Join const& collisions, + MftTracksWCollsMC const& tracks, + aod::McCollisions const& mccollisions, + aod::McParticles const& particles, + ExtBCs const& bcs) + { + processReassocEfficiency(collisions, tracks, mccollisions, particles, bcs); + } + + PROCESS_SWITCH(DndetaMFTPbPb, processReassocEfficiencyCentFT0C, "Process collision-reassociation efficiency based on track propagation DCA 3D (in FT0C centrality bins)", false); + + /// @brief process function to calculate signal loss based on MC + void processEventAndSignalLossCentFT0C(CollsMCExtra::iterator const& mcCollision, + soa::SmallGroups> const& collisions, + FiltParticles const& particles) + { + registry.fill(HIST("Events/Centrality/hNRecCollsSigEvtLoss"), 1.f, collisions.size()); + + if (gConf.cfgUseInelgt0 && !mcCollision.isInelGt0()) { + return; + } + if (eventCuts.useZVtxCutMC && (std::abs(mcCollision.posZ()) >= eventCuts.maxZvtx)) { + return; + } + + bool gtZeroColl = false; + auto maxNcontributors = -1; + auto centrality = -1; + for (auto const& collision : collisions) { + if (!isGoodEvent(collision)) { + continue; + } + if (std::abs(collision.posZ()) >= eventCuts.maxZvtx) { + continue; + } + if (maxNcontributors < collision.numContrib()) { + maxNcontributors = collision.numContrib(); + centrality = getRecoCent(collision); + } else { + continue; + } + gtZeroColl = true; + } + + auto perCollMCsample = mcSample->sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); + auto multMCNParticlesEtaMFT = countPart(perCollMCsample); + + registry.fill(HIST("Events/Centrality/EvtSigLossStatus"), 1., centrality); + registry.fill(HIST("Events/Centrality/hMultGenVsCent"), centrality, mcCollision.multMCFT0C()); + registry.fill(HIST("Events/Centrality/hMultGenVsCentNParticlesEta05"), centrality, mcCollision.multMCNParticlesEta05()); + registry.fill(HIST("Events/Centrality/hMultGenVsCentNParticlesEtaMFT"), centrality, multMCNParticlesEtaMFT); + + if (collisions.size() == 0) { + registry.fill(HIST("Events/Centrality/EvtSigLossStatus"), 3., centrality); + } + + if (gtZeroColl) { + registry.fill(HIST("Events/Centrality/EvtSigLossStatus"), 2., centrality); + registry.fill(HIST("Events/Centrality/hMultGenVsCentRec"), centrality, mcCollision.multMCFT0C()); + registry.fill(HIST("Events/Centrality/hMultGenVsCentRecNParticlesEta05"), centrality, mcCollision.multMCNParticlesEta05()); + registry.fill(HIST("Events/Centrality/hMultGenVsCentRecNParticlesEtaMFT"), centrality, multMCNParticlesEtaMFT); + } + + for (auto const& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { + continue; + } + + float phi = particle.phi(); + o2::math_utils::bringTo02Pi(phi); + if (phi < Czero || TwoPI < phi) { + continue; + } + + registry.fill(HIST("Tracks/Centrality/EtaCentVsMultGen_t"), particle.eta(), centrality, mcCollision.multMCFT0C()); + registry.fill(HIST("Tracks/Centrality/EtaGen_t"), particle.eta(), centrality); + + if (gtZeroColl) { + float phi = particle.phi(); + o2::math_utils::bringTo02Pi(phi); + if (phi < Czero || TwoPI < phi) { + continue; + } + registry.fill(HIST("Tracks/Centrality/EtaCentVsMultGen"), particle.eta(), centrality, mcCollision.multMCFT0C()); + registry.fill(HIST("Tracks/Centrality/EtaGen"), particle.eta(), centrality); + } + } + } + + PROCESS_SWITCH(DndetaMFTPbPb, processEventAndSignalLossCentFT0C, "Signal/event loss based on MC (in FT0C centrality bins)", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)