Merge branch 'AliceO2Group:master' into master

Author: CyF1rmiana

ALICE3/Tasks/alice3-lutmaker.cxx

@@ -25,11 +25,11 @@ using namespace framework::expressions;
 void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
 {
   std::vector<ConfigParamSpec> options{
-    {"lut-el", VariantType::Int, 1, {"LUT input for the Electron PDG code"}},
-    {"lut-mu", VariantType::Int, 1, {"LUT input for the Muon PDG code"}},
+    {"lut-el", VariantType::Int, 0, {"LUT input for the Electron PDG code"}},
+    {"lut-mu", VariantType::Int, 0, {"LUT input for the Muon PDG code"}},
     {"lut-pi", VariantType::Int, 1, {"LUT input for the Pion PDG code"}},
-    {"lut-ka", VariantType::Int, 1, {"LUT input for the Kaon PDG code"}},
-    {"lut-pr", VariantType::Int, 1, {"LUT input for the Proton PDG code"}},
+    {"lut-ka", VariantType::Int, 0, {"LUT input for the Kaon PDG code"}},
+    {"lut-pr", VariantType::Int, 0, {"LUT input for the Proton PDG code"}},
     {"lut-tr", VariantType::Int, 0, {"LUT input for the Triton PDG code"}},
     {"lut-de", VariantType::Int, 0, {"LUT input for the Deuteron PDG code"}},
     {"lut-he", VariantType::Int, 0, {"LUT input for the Helium3 PDG code"}}};
@@ -153,6 +153,7 @@ struct Alice3LutMaker {
     histos.add("QA/CovMat_sigmaSnp", "sigmaSnp" + commonTitle, kTH3F, {axisPt, axisEta, axissigmaSnp});
     histos.add("QA/CovMat_sigmaTgl", "sigmaTgl" + commonTitle, kTH3F, {axisPt, axisEta, axissigmaTgl});
     histos.add("QA/CovMat_sigma1Pt", "sigma1Pt" + commonTitle, kTH3F, {axisPt, axisEta, axissigma1Pt});
+    histos.add("QA/sigma1Pt", "sigma1Pt" + commonTitle, kTH3F, {axisPt, axisEta, axissigma1Pt});
     histos.add("QA/CovMat_rhoZY", "rhoZY" + commonTitle, kTH3F, {axisPt, axisEta, axisrhoZY});
     histos.add("QA/CovMat_rhoSnpY", "rhoSnpY" + commonTitle, kTH3F, {axisPt, axisEta, axisrhoSnpY});
     histos.add("QA/CovMat_rhoSnpZ", "rhoSnpZ" + commonTitle, kTH3F, {axisPt, axisEta, axisrhoSnpZ});
@@ -262,6 +263,7 @@ struct Alice3LutMaker {
       histos.fill(HIST("QA/CovMat_sigmaSnp"), mcParticle.pt(), mcParticle.eta(), track.sigmaSnp());
       histos.fill(HIST("QA/CovMat_sigmaTgl"), mcParticle.pt(), mcParticle.eta(), track.sigmaTgl());
       histos.fill(HIST("QA/CovMat_sigma1Pt"), mcParticle.pt(), mcParticle.eta(), track.sigma1Pt());
+      histos.fill(HIST("QA/sigma1Pt"), mcParticle.pt(), mcParticle.eta(), std::abs(track.signed1Pt()) - 1. / mcParticle.pt());
       histos.fill(HIST("QA/CovMat_rhoZY"), mcParticle.pt(), mcParticle.eta(), track.rhoZY());
       histos.fill(HIST("QA/CovMat_rhoSnpY"), mcParticle.pt(), mcParticle.eta(), track.rhoSnpY());
       histos.fill(HIST("QA/CovMat_rhoSnpZ"), mcParticle.pt(), mcParticle.eta(), track.rhoSnpZ());

CODEOWNERS

@@ -38,7 +38,7 @@
 /PWGEM @alibuild @feisenhu @dsekihat @ivorobye
 /PWGEM/Dilepton @alibuild @mikesas @rbailhac @dsekihat @ivorobye @feisenhu @hscheid
 /PWGEM/PhotonMeson @alibuild @mikesas @rbailhac @m-c-danisch @novitzky @mhemmer-cern @dsekihat
-/PWGHF @alibuild @vkucera @fcolamar @fgrosa @fcatalan92 @mfaggin @mmazzilli @deepathoms @NicoleBastid @hahassan7 @jpxrk @apalasciano @zhangbiao-phy @gluparel
+/PWGHF @alibuild @vkucera @fcolamar @fgrosa @fcatalan92 @mfaggin @mmazzilli @deepathoms @NicoleBastid @hahassan7 @jpxrk @apalasciano @zhangbiao-phy @gluparel @stefanopolitano
 # PWG-LF
 /PWGLF @alibuild @sustripathy @skundu692
 /PWGLF/DataModel @alibuild @sustripathy @skundu692 @gbencedi @abmodak @fmazzasc @maciacco @dmallick2 @smaff92 @ercolessi @romainschotter

Common/Core/fwdtrackUtilities.h

@@ -42,6 +42,7 @@ enum class propagationPoint : int {
   kToRabs = 2,
 };
 using SMatrix55 = ROOT::Math::SMatrix<double, 5, 5, ROOT::Math::MatRepSym<double, 5>>;
+using SMatrix55Std = ROOT::Math::SMatrix<double, 5>;
 using SMatrix5 = ROOT::Math::SVector<double, 5>;
 
 /// propagate fwdtrack to a certain point.
@@ -98,6 +99,48 @@ o2::dataformats::GlobalFwdTrack propagateMuon(TFwdTrack const& muon, TCollision
 
   return propmuon;
 }
+
+template <typename TFwdTrack, typename TMFTTrack>
+o2::dataformats::GlobalFwdTrack refitGlobalMuonCov(TFwdTrack const& muon, TMFTTrack const& mft)
+{
+  auto muonCov = muon.getCovariances();
+  auto mftCov = mft.getCovariances();
+
+  SMatrix55Std jacob = ROOT::Math::SMatrixIdentity();
+  auto tl = muon.getTgl();
+  auto invQPt = muon.getInvQPt();
+  jacob(4, 3) = tl / (invQPt * std::sqrt(1 + tl * tl));
+  jacob(4, 4) = -std::sqrt(1 + tl * tl) / (invQPt * invQPt);
+
+  auto covQP = ROOT::Math::Similarity(jacob, muonCov);
+  mftCov(4, 0) = 0;
+  mftCov(4, 1) = 0;
+  mftCov(4, 2) = 0;
+  mftCov(4, 3) = 0;
+
+  mftCov(0, 4) = 0;
+  mftCov(1, 4) = 0;
+  mftCov(2, 4) = 0;
+  mftCov(3, 4) = 0;
+  mftCov(4, 4) = covQP(4, 4);
+
+  SMatrix55Std jacobInv = ROOT::Math::SMatrixIdentity();
+  auto qp = std::sqrt(1 + tl * tl) / invQPt;
+  auto tlMFT = mft.getTgl();
+  jacobInv(4, 3) = tlMFT / (qp * std::sqrt(1 + tlMFT * tlMFT));
+  jacobInv(4, 4) = -std::sqrt(1 + tlMFT * tlMFT) / (qp * qp);
+  auto globalCov = ROOT::Math::Similarity(jacobInv, mftCov);
+
+  auto invQPtGlob = std::sqrt(1 + tlMFT * tlMFT) / qp;
+
+  o2::dataformats::GlobalFwdTrack globalTrack;
+  globalTrack.setParameters(mft.getParameters());
+  globalTrack.setInvQPt(invQPtGlob);
+  globalTrack.setCovariances(globalCov);
+
+  return globalTrack;
+}
+
 } // namespace fwdtrackutils
 } // namespace o2::aod
 

Common/Tasks/centralityStudy.cxx

@@ -132,6 +132,7 @@ struct centralityStudy {
   ConfigurableAxis axisMultPVContributors{"axisMultPVContributors", {200, 0, 6000}, "Number of PV Contributors"};
   ConfigurableAxis axisMultGlobalTracks{"axisMultGlobalTracks", {500, 0, 5000}, "Number of global tracks"};
   ConfigurableAxis axisMultMFTTracks{"axisMultMFTTracks", {500, 0, 5000}, "Number of MFT tracks"};
+  ConfigurableAxis axisMultMCCounts{"axisMultMCCounts", {1000, 0, 5000}, "N_{ch}"};
 
   ConfigurableAxis axisTrackOccupancy{"axisTrackOccupancy", {50, 0, 5000}, "Track occupancy"};
   ConfigurableAxis axisFT0COccupancy{"axisFT0COccupancy", {50, 0, 80000}, "FT0C occupancy"};
@@ -242,12 +243,19 @@ struct centralityStudy {
     }
 
     if (doprocessCollisionsWithResolutionStudy) {
+      // histograms with detector signals
       histos.add("hImpactParameterVsFT0A", "hImpactParameterVsFT0A", kTH2F, {axisMultFT0A, axisImpactParameter});
       histos.add("hImpactParameterVsFT0C", "hImpactParameterVsFT0C", kTH2F, {axisMultFT0C, axisImpactParameter});
       histos.add("hImpactParameterVsFT0M", "hImpactParameterVsFT0M", kTH2F, {axisMultFT0M, axisImpactParameter});
       histos.add("hImpactParameterVsFV0A", "hImpactParameterVsFV0A", kTH2F, {axisMultFV0A, axisImpactParameter});
       histos.add("hImpactParameterVsNMFTTracks", "hImpactParameterVsNMFTTracks", kTH2F, {axisMultMFTTracks, axisImpactParameter});
       histos.add("hImpactParameterVsNTPV", "hImpactParameterVsNTPV", kTH2F, {axisMultPVContributors, axisImpactParameter});
+
+      // histograms with actual MC counts in each region
+      histos.add("hImpactParameterVsMCFT0A", "hImpactParameterVsMCFT0A", kTH2F, {axisMultMCCounts, axisImpactParameter});
+      histos.add("hImpactParameterVsMCFT0C", "hImpactParameterVsMCFT0C", kTH2F, {axisMultMCCounts, axisImpactParameter});
+      histos.add("hImpactParameterVsMCFT0M", "hImpactParameterVsMCFT0M", kTH2F, {axisMultMCCounts, axisImpactParameter});
+      histos.add("hImpactParameterVsMCFV0A", "hImpactParameterVsMCFV0A", kTH2F, {axisMultMCCounts, axisImpactParameter});
     }
 
     if (doOccupancyStudyVsRawValues2d) {
@@ -696,6 +704,11 @@ struct centralityStudy {
         histos.fill(HIST("hImpactParameterVsFV0A"), multFV0A, mcCollision.impactParameter());
         histos.fill(HIST("hImpactParameterVsNMFTTracks"), mftNtracks, mcCollision.impactParameter());
         histos.fill(HIST("hImpactParameterVsNTPV"), multNTracksPV, mcCollision.impactParameter());
+
+        histos.fill(HIST("hImpactParameterVsMCFT0A"), mcCollision.multMCFT0A(), mcCollision.impactParameter());
+        histos.fill(HIST("hImpactParameterVsMCFT0C"), mcCollision.multMCFT0C(), mcCollision.impactParameter());
+        histos.fill(HIST("hImpactParameterVsMCFT0M"), (mcCollision.multMCFT0A() + mcCollision.multMCFT0C()), mcCollision.impactParameter());
+        histos.fill(HIST("hImpactParameterVsMCFV0A"), mcCollision.multMCFV0A(), mcCollision.impactParameter());
       }
     }
 

PWGCF/Flow/Tasks/flowTask.cxx

@@ -36,6 +36,7 @@
 #include "Framework/RunningWorkflowInfo.h"
 #include "Framework/runDataProcessing.h"
 #include <CCDB/BasicCCDBManager.h>
+#include <DataFormatsParameters/GRPMagField.h>
 
 #include "TList.h"
 #include <TF1.h>
@@ -141,6 +142,13 @@ struct FlowTask {
     TF1* fMultMultV0ACutHigh = nullptr;
     TF1* fT0AV0AMean = nullptr;
     TF1* fT0AV0ASigma = nullptr;
+    // for TPC sector boundary
+    O2_DEFINE_CONFIGURABLE(cfgShowTPCsectorOverlap, bool, true, "Draw TPC sector overlap")
+    O2_DEFINE_CONFIGURABLE(cfgRejectionTPCsectorOverlap, bool, false, "rejection for TPC sector overlap")
+    O2_DEFINE_CONFIGURABLE(cfgMagnetField, std::string, "GLO/Config/GRPMagField", "CCDB path to Magnet field object")
+    ConfigurableAxis axisPhiMod{"axisPhiMod", {100, 0, constants::math::PI / 9}, "fmod(#varphi,#pi/9)"};
+    TF1* fPhiCutLow = nullptr;
+    TF1* fPhiCutHigh = nullptr;
   } cfgFuncParas;
 
   ConfigurableAxis axisPtHist{"axisPtHist", {100, 0., 10.}, "pt axis for histograms"};
@@ -283,6 +291,8 @@ struct FlowTask {
     registry.add("hEta", "#eta distribution", {HistType::kTH1D, {axisEta}});
     registry.add("hPt", "p_{T} distribution before cut", {HistType::kTH1D, {axisPtHist}});
     registry.add("hPtRef", "p_{T} distribution after cut", {HistType::kTH1D, {axisPtHist}});
+    registry.add("pt_phi_bef", "before cut;p_{T};#phi_{modn}", {HistType::kTH2D, {axisPt, cfgFuncParas.axisPhiMod}});
+    registry.add("pt_phi_aft", "after cut;p_{T};#phi_{modn}", {HistType::kTH2D, {axisPt, cfgFuncParas.axisPhiMod}});
     registry.add("hChi2prTPCcls", "#chi^{2}/cluster for the TPC track segment", {HistType::kTH1D, {{100, 0., 5.}}});
     registry.add("hChi2prITScls", "#chi^{2}/cluster for the ITS track", {HistType::kTH1D, {{100, 0., 50.}}});
     registry.add("hnTPCClu", "Number of found TPC clusters", {HistType::kTH1D, {{100, 40, 180}}});
@@ -501,6 +511,11 @@ struct FlowTask {
       cfgFuncParas.fT0AV0ASigma->SetParameters(463.4144, 6.796509e-02, -9.097136e-07, 7.971088e-12, -2.600581e-17);
     }
 
+    if (cfgFuncParas.cfgShowTPCsectorOverlap) {
+      cfgFuncParas.fPhiCutLow = new TF1("fPhiCutLow", "0.06/x+pi/18.0-0.06", 0, 100);
+      cfgFuncParas.fPhiCutHigh = new TF1("fPhiCutHigh", "0.1/x+pi/18.0+0.06", 0, 100);
+    }
+
     if (cfgTrackDensityCorrUse) {
       std::vector<double> pTEffBins = {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0};
       hFindPtBin = new TH1D("hFindPtBin", "hFindPtBin", pTEffBins.size() - 1, &pTEffBins[0]);
@@ -730,12 +745,49 @@ struct FlowTask {
     return 1;
   }
 
+  int getMagneticField(uint64_t timestamp)
+  {
+    static o2::parameters::GRPMagField* grpo = nullptr;
+    if (grpo == nullptr) {
+      grpo = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(cfgFuncParas.cfgMagnetField, timestamp);
+      if (grpo == nullptr) {
+        LOGF(fatal, "GRP object not found in %s for timestamp %llu", cfgFuncParas.cfgMagnetField.value.c_str(), timestamp);
+        return 0;
+      }
+      LOGF(info, "Retrieved GRP from %s for timestamp %llu with magnetic field of %d kG", cfgFuncParas.cfgMagnetField.value.c_str(), timestamp, grpo->getNominalL3Field());
+    }
+    return grpo->getNominalL3Field();
+  }
+
   template <typename TTrack>
   bool trackSelected(TTrack track)
   {
     return ((track.tpcNClsFound() >= cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgCutTPCCrossedRows) && (track.itsNCls() >= cfgCutITSclu));
   }
 
+  template <typename TTrack>
+  bool rejectionTPCoverlap(TTrack track, const int field)
+  {
+    double phimodn = track.phi();
+    if (field < 0) // for negative polarity field
+      phimodn = o2::constants::math::TwoPI - phimodn;
+    if (track.sign() < 0) // for negative charge
+      phimodn = o2::constants::math::TwoPI - phimodn;
+    if (phimodn < 0)
+      LOGF(warning, "phi < 0: %g", phimodn);
+
+    float middle = o2::constants::math::TwoPI / 18.0;
+    phimodn += middle; // to center gap in the middle
+    phimodn = fmod(phimodn, o2::constants::math::TwoPI / 9.0);
+    registry.fill(HIST("pt_phi_bef"), track.pt(), phimodn);
+    if (cfgFuncParas.cfgRejectionTPCsectorOverlap) {
+      if (phimodn < cfgFuncParas.fPhiCutHigh->Eval(track.pt()) && phimodn > cfgFuncParas.fPhiCutLow->Eval(track.pt()))
+        return false; // reject track
+    }
+    registry.fill(HIST("pt_phi_aft"), track.pt(), phimodn);
+    return true;
+  }
+
   void initHadronicRate(aod::BCsWithTimestamps::iterator const& bc)
   {
     if (mRunNumber == bc.runNumber()) {
@@ -844,9 +896,14 @@ struct FlowTask {
     // track weights
     float weff = 1, wacc = 1;
     double nTracksCorrected = 0;
+    int magnetfield = 0;
     float independent = cent;
     if (cfgUseNch)
       independent = static_cast<float>(tracks.size());
+    if (cfgFuncParas.cfgShowTPCsectorOverlap) {
+      // magnet field dependence cut
+      magnetfield = getMagneticField(bc.timestamp());
+    }
 
     double psi2Est = 0, psi3Est = 0, psi4Est = 0;
     float wEPeff = 1;
@@ -879,6 +936,8 @@ struct FlowTask {
     for (const auto& track : tracks) {
       if (!trackSelected(track))
         continue;
+      if (cfgFuncParas.cfgShowTPCsectorOverlap && !rejectionTPCoverlap(track, magnetfield))
+        continue;
       bool withinPtPOI = (cfgCutPtPOIMin < track.pt()) && (track.pt() < cfgCutPtPOIMax); // within POI pT range
       bool withinPtRef = (cfgCutPtRefMin < track.pt()) && (track.pt() < cfgCutPtRefMax); // within RF pT range
       if (cfgOutputNUAWeights) {

PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx

@@ -1241,7 +1241,7 @@ struct HfDataCreatorCharmHadPiReduced {
         trackParCov1.propagateTo(secondaryVertexCharm[0], bz);
         df2.getTrack(0).getPxPyPzGlo(pVec0);
         df2.getTrack(1).getPxPyPzGlo(pVec1);
-        pVecCharm = RecoDecay::pVec(pVec0, pVec1);
+        pVecCharm = RecoDecay::pVec(pVec0, pVec1, pVec2);
         trackParCovCharmHad = df2.createParentTrackParCov();
         trackParCovCharmHad.setAbsCharge(0); // to be sure
       }
@@ -1279,14 +1279,18 @@ struct HfDataCreatorCharmHadPiReduced {
         }
 
         // reject pi D with same sign as D
-        if constexpr (decChannel == DecayChannel::B0ToDminusPi || decChannel == DecayChannel::BsToDsminusPi || decChannel == DecayChannel::LbToLcplusPi || decChannel == DecayChannel::B0ToDstarPi) { // D∓ → π∓ K± π∓ and Ds∓ → K∓ K± π∓ and Lc∓ → p∓ K± π∓ and D*+ → D0 π+
+        if constexpr (decChannel == DecayChannel::B0ToDminusPi || decChannel == DecayChannel::BsToDsminusPi || decChannel == DecayChannel::LbToLcplusPi) { // D∓ → π∓ K± π∓ and Ds∓ → K∓ K± π∓ and Lc∓ → p∓ K± π∓
           if (trackPion.sign() * charmHadDauTracks[0].sign() > 0) {
             continue;
           }
         } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // D0(bar) → K± π∓
           if (!((candC.isSelD0() >= hfflagConfigurations.selectionFlagD0 && trackPion.sign() < 0) || (candC.isSelD0bar() >= hfflagConfigurations.selectionFlagD0bar && trackPion.sign() > 0))) {
             continue;
           }
+        } else if constexpr (decChannel == DecayChannel::B0ToDstarPi) { // D*+ → D0 π+
+          if (trackPion.sign() * charmHadDauTracks.back().sign() > 0) {
+            continue;
+          }
         }
 
         // apply selections on pion tracks

PWGHF/D2H/Tasks/taskB0Reduced.cxx

@@ -101,7 +101,7 @@ DECLARE_SOA_COLUMN(FlagWrongCollision, flagWrongCollision, int8_t);
 } // namespace hf_cand_b0_lite
 
 DECLARE_SOA_TABLE(HfRedCandB0Lites, "AOD", "HFREDCANDB0LITE", //! Table with some B0 properties
-                                                              // B meson features                  hf_cand_b0_lite::M,
+                                                              // B meson features
                   hf_cand_b0_lite::M,
                   hf_cand_b0_lite::Pt,
                   hf_cand_b0_lite::Eta,
@@ -616,7 +616,7 @@ struct HfTaskB0Reduced {
           candidate.decayLengthXY(),
           candidate.decayLengthNormalised(),
           candidate.decayLengthXYNormalised(),
-          candidate.impactParameterProngSqSum(),
+          candidate.impactParameterProduct(),
           candidate.maxNormalisedDeltaIP(),
           candidateMlScoreSig,
           candidate.isSelB0ToDPi(),
@@ -644,12 +644,12 @@ struct HfTaskB0Reduced {
           prong0MlScorePrompt,
           prong0MlScoreNonprompt,
           // pion features
-          candidate.ptProng1(),
+          candidate.ptProng2(),
           std::abs(RecoDecay::eta(prongBachPi.pVector())),
           prongBachPi.itsNCls(),
           prongBachPi.tpcNClsCrossedRows(),
           prongBachPi.tpcChi2NCl(),
-          candidate.impactParameter1(),
+          candidate.impactParameter2(),
           prongBachPi.tpcNSigmaPi(),
           prongBachPi.tofNSigmaPi(),
           prongBachPi.tpcTofNSigmaPi(),