[PWGHF] addition of task: single muon as a function of multiplicity

PWGHF/HFL/Tasks/CMakeLists.txt

@@ -24,6 +24,11 @@ o2physics_add_dpl_workflow(task-single-muon
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(task-single-muon-mult
+                    SOURCES taskSingleMuonMult.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(task-single-muon-reader
                     SOURCES taskSingleMuonReader.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGDQCore

PWGHF/HFL/Tasks/taskSingleMuonMult.cxx

@@ -0,0 +1,265 @@
+// 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.
+
+/// \file taskSingleMuonMult.cxx
+/// \brief Task used to study the Open heavy flavour decay muon production as a function of multiplicity.
+/// \author Md Samsul Islam <md.samsul.islam@cern.ch>, IITB
+
+#include <cmath>
+
+#include "CommonConstants/MathConstants.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/runDataProcessing.h"
+#include "ReconstructionDataFormats/TrackFwd.h"
+
+#include "Common/Core/TrackSelection.h"
+#include "Common/Core/trackUtilities.h"
+#include "Common/DataModel/Centrality.h"
+#include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/Multiplicity.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+
+using namespace o2;
+using namespace o2::aod;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::aod::fwdtrack;
+
+struct HfTaskSingleMuonMult {
+
+  // enum for event selection bins
+  enum EventSelection {
+    AllEvents = 0,
+    Sel8,
+    VtxZAfterSel,
+    NEventSelection
+  };
+
+  // enum for muon track selection bins
+  enum MuonSelection {
+    NoCut = 0,
+    EtaCut,
+    RAbsorbCut,
+    PDcaCut,
+    Chi2Cut,
+    NMuonSelection
+  };
+
+  Configurable<float> zVtxMax{"zVtxMax", 10., "maxium z of primary vertex [cm]"};
+  Configurable<float> ptTrackMin{"ptTrackMin", 0.15, "minimum pt of tracks"};
+  Configurable<float> etaTrackMax{"etaTrackMax", 0.8, "maximum pseudorapidity of tracks"};
+  Configurable<float> etaMin{"etaMin", -3.6, "minimum pseudorapidity"};
+  Configurable<float> etaMax{"etaMax", -2.5, "maximum pseudorapidity"};
+  Configurable<float> pDcaMin{"pDcaMin", 324., "p*DCA value for small RAbsorb"};
+  Configurable<float> pDcaMax{"pDcaMax", 594., "p*DCA value for large RAbsorb"};
+  Configurable<float> rAbsorbMin{"rAbsorbMin", 17.6, "R at absorber end minimum value"};
+  Configurable<float> rAbsorbMax{"rAbsorbMax", 89.5, "R at absorber end maximum value"};
+  Configurable<float> rAbsorbMid{"rAbsorbMid", 26.5, "R at absorber end split point for different p*DCA selections"};
+  Configurable<bool> reduceOrphMft{"reduceOrphMft", true, "reduce orphan MFT tracks"};
+
+  using MyCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Ms>;
+  using MyMuons = soa::Join<aod::FwdTracks, aod::FwdTracksDCA>;
+  using MyMcMuons = soa::Join<aod::FwdTracks, aod::McFwdTrackLabels, aod::FwdTracksDCA>;
+  using MyTracks = soa::Filtered<soa::Join<aod::FullTracks, aod::TracksIU, aod::TracksDCA, aod::TrackSelection>>;
+
+  // Filter Global Track for Multiplicty
+  Filter trackFilter = ((nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin));
+
+  // Number the types of muon tracks
+  static constexpr uint8_t NTrackTypes{ForwardTrackTypeEnum::MCHStandaloneTrack + 1};
+
+  HistogramRegistry registry{"registry"};
+
+  void init(InitContext&)
+  {
+    AxisSpec axisCent = {101, -0.5, 100.5, "centrality"};
+    AxisSpec axisEvent{NEventSelection, 0, NEventSelection, "Event Selection"};
+    AxisSpec axisVtxZ{80, -20., 20., "#it{z}_{vtx} (cm)"};
+    AxisSpec axisMuon{NMuonSelection, 0, NMuonSelection, "Muon Selection"};
+    AxisSpec axisNCh{500, 0.5, 500.5, "#it{N}_{ch}"};
+    AxisSpec axisNMu{20, -0.5, 19.5, "#it{N}_{#mu}"};
+    AxisSpec axisPt{1000, 0., 500., "#it{p}_{T} (GeV/#it{c})"};
+    AxisSpec axisEta{250, -5., 5., "#it{#eta}"};
+    AxisSpec axisTheta{500, 170., 180., "#it{#theta}"};
+    AxisSpec axisRAbsorb{1000, 0., 100., "#it{R}_{Absorb} (cm)"};
+    AxisSpec axisDCA{500, 0., 5., "#it{DCA}_{xy} (cm)"};
+    AxisSpec axisChi2MatchMCHMFT{1000, 0., 1000., "MCH-MFT matching #chi^{2}"};
+    AxisSpec axisSign{5, -2.5, 2.5, "Charge"};
+    AxisSpec axisPDca{100000, 0, 100000, "#it{p} #times DCA (GeV/#it{c} * cm)"};
+    AxisSpec axisDCAx{1000, -5., 5., "#it{DCA}_{x or y} (cm)"};
+    AxisSpec axisEtaDif{200, -2., 2., "#it{#eta} diff"};
+    AxisSpec axisDeltaPt{10000, -50, 50, "#Delta #it{p}_{T} (GeV/#it{c})"};
+    AxisSpec axisTrackType{8, -1.5, 6.5, "TrackType"};
+    AxisSpec axisPtDif{200, -2., 2., "#it{p}_{T} diff (GeV/#it{c})"};
+
+    registry.add("hCentrality", "Centrality Percentile", {HistType::kTH1F, {axisCent}});
+    registry.add("hEventSel", " Number of Events", {HistType::kTH1F, {axisEvent}});
+    registry.add("hNch", "Charged Particle Multiplicity", {HistType::kTH1F, {axisNCh}});
+    registry.add("hVtxZBeforeSel", "Z-vertex distribution before zVtx Cut", {HistType::kTH1F, {axisVtxZ}});
+    registry.add("hVtxZAfterSel", "Z-vertex distribution after zVtx Cut", {HistType::kTH1F, {axisVtxZ}});
+
+    registry.add("hMuonSel", "Selection of muon tracks at various kinematic cuts", {HistType::kTH1F, {axisMuon}});
+    registry.add("hMuBeforeMatchMFT", "Muon information before any Kinemeatic cuts applied", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisTrackType}, 10});
+    registry.add("hMuBeforeAccCuts", "Muon information before applying Acceptance cuts", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisTrackType}, 10});
+    registry.add("h3DCABeforeAccCuts", "DCAx,DCAy,DCAz information before Acceptance cuts", {HistType::kTH3F, {axisDCAx, axisDCAx, axisTrackType}});
+    registry.add("hMuDeltaPtBeforeAccCuts", "Muon information with DeltaPt before applying Acceptance cuts", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisDeltaPt}, 10});
+    registry.add("hMuAfterEtaCuts", "Muon information after applying Eta cuts", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisTrackType}, 10});
+    registry.add("hMuAfterRAbsorbCuts", "Muon information after applying RAbsorb cuts", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisTrackType}, 10});
+    registry.add("hMuAfterPdcaCuts", "Muon information after applying Pdca cuts", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisTrackType}, 10});
+    registry.add("hMuAfterAccCuts", "Muon information after applying all Kinematic cuts", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisTrackType}, 10});
+    registry.add("h3DCAAfterAccCuts", "DCAx,DCAy,DCAz information after Acceptance cuts", {HistType::kTH3F, {axisDCAx, axisDCAx, axisTrackType}});
+    registry.add("hMuDeltaPtAfterAccCuts", "Muon information with DeltaPt after applying Acceptance cuts", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisTheta, axisRAbsorb, axisDCA, axisPDca, axisChi2MatchMCHMFT, axisDeltaPt}, 10});
+
+    registry.add("hTHnTrk", "Muon information with multiplicity", {HistType::kTHnSparseF, {axisCent, axisNCh, axisPt, axisEta, axisSign}, 5});
+    registry.add("h3MultNchNmu", "Number of muons and multiplicity", {HistType::kTH3F, {axisCent, axisNCh, axisNMu}});
+    registry.add("hMultNchNmuTrackType", "Number of muons with different types and multiplicity", {HistType::kTHnSparseF, {axisCent, axisNCh, axisNMu, axisTrackType}, 4});
+
+    auto hEvstat = registry.get<TH1>(HIST("hEventSel"));
+    auto* xEv = hEvstat->GetXaxis();
+    xEv->SetBinLabel(AllEvents + 1, "All events");
+    xEv->SetBinLabel(Sel8 + 1, "sel8");
+    xEv->SetBinLabel(VtxZAfterSel + 1, "VtxZAfterSel");
+
+    auto hMustat = registry.get<TH1>(HIST("hMuonSel"));
+    auto* xMu = hMustat->GetXaxis();
+    xMu->SetBinLabel(NoCut + 1, "noCut");
+    xMu->SetBinLabel(EtaCut + 1, "etaCut");
+    xMu->SetBinLabel(RAbsorbCut + 1, "RAbsorbCut");
+    xMu->SetBinLabel(PDcaCut + 1, "pDcaCut");
+    xMu->SetBinLabel(Chi2Cut + 1, "chi2Cut");
+  }
+
+  void process(MyCollisions::iterator const& collision,
+               MyTracks const& tracks,
+               MyMuons const& muons)
+  {
+    registry.fill(HIST("hEventSel"), AllEvents);
+
+    if (!collision.sel8()) {
+      return;
+    }
+    registry.fill(HIST("hEventSel"), Sel8);
+    registry.fill(HIST("hVtxZBeforeSel"), collision.posZ());
+
+    if (std::abs(collision.posZ()) > zVtxMax) {
+      return;
+    }
+    registry.fill(HIST("hEventSel"), VtxZAfterSel);
+    registry.fill(HIST("hVtxZAfterSel"), collision.posZ());
+
+    // T0M centrality
+    const auto cent = collision.centFT0M();
+    registry.fill(HIST("hCentrality"), cent);
+
+    // Charged particles
+    for (const auto& track : tracks) {
+      if (!track.isGlobalTrack()) {
+        continue;
+      }
+    }
+
+    auto nCh{tracks.size()};
+    if (nCh < 1) {
+      return;
+    }
+    registry.fill(HIST("hNch"), nCh);
+
+    for (const auto& track : tracks) {
+      registry.fill(HIST("hTHnTrk"), cent, nCh, track.pt(), track.eta(), track.sign());
+    }
+
+    // muons per event
+    int nMu{0};
+    int nMuType[NTrackTypes] = {0};
+
+    for (const auto& muon : muons) {
+      const auto pt{muon.pt()}, eta{muon.eta()}, theta{90.0f - ((std::atan(muon.tgl())) * constants::math::Rad2Deg)}, pDca{muon.pDca()}, rAbsorb{muon.rAtAbsorberEnd()}, chi2{muon.chi2MatchMCHMFT()};
+      const auto dcaXY{RecoDecay::sqrtSumOfSquares(muon.fwdDcaX(), muon.fwdDcaY())};
+      const auto muTrackType{muon.trackType()};
+
+      registry.fill(HIST("hMuBeforeMatchMFT"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, muTrackType);
+
+      // histograms before the acceptance cuts
+      registry.fill(HIST("hMuonSel"), NoCut);
+      registry.fill(HIST("hMuBeforeAccCuts"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, muTrackType);
+      registry.fill(HIST("h3DCABeforeAccCuts"), muon.fwdDcaX(), muon.fwdDcaY(), muTrackType);
+
+      if (muon.has_matchMCHTrack()) {
+        auto muonType3 = muon.template matchMCHTrack_as<MyMuons>();
+        auto dpt = muonType3.pt() - pt;
+        if (muTrackType == ForwardTrackTypeEnum::GlobalMuonTrack) {
+          registry.fill(HIST("hMuDeltaPtBeforeAccCuts"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, dpt);
+        }
+      }
+
+      // Apply various standard muon acceptance cuts
+      // eta cuts
+      if ((eta >= etaMax) || (eta < etaMin)) {
+        continue;
+      }
+      registry.fill(HIST("hMuonSel"), EtaCut);
+      registry.fill(HIST("hMuAfterEtaCuts"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, muTrackType);
+
+      // Rabsorb cuts
+      if ((rAbsorb < rAbsorbMin) || (rAbsorb >= rAbsorbMax)) {
+        continue;
+      }
+      registry.fill(HIST("hMuonSel"), RAbsorbCut);
+      registry.fill(HIST("hMuAfterRAbsorbCuts"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, muTrackType);
+
+      if ((rAbsorb < rAbsorbMid) && (pDca >= pDcaMin)) {
+        continue;
+      }
+      if ((rAbsorb >= rAbsorbMid) && (pDca >= pDcaMax)) {
+        continue;
+      }
+      registry.fill(HIST("hMuonSel"), PDcaCut);
+      registry.fill(HIST("hMuAfterPdcaCuts"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, muTrackType);
+
+      //  MCH-MFT matching chi2
+      if (muon.chi2() >= 1e6) {
+        continue;
+      }
+      registry.fill(HIST("hMuonSel"), Chi2Cut);
+
+      // histograms after acceptance cuts
+      registry.fill(HIST("hMuAfterAccCuts"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, muTrackType);
+      registry.fill(HIST("h3DCAAfterAccCuts"), muon.fwdDcaX(), muon.fwdDcaY(), muTrackType);
+      nMu++;
+      nMuType[muTrackType]++;
+
+      if (muon.has_matchMCHTrack()) {
+        auto muonType3 = muon.template matchMCHTrack_as<MyMuons>();
+        auto dpt = muonType3.pt() - pt;
+
+        if (muTrackType == ForwardTrackTypeEnum::GlobalMuonTrack) {
+          registry.fill(HIST("hMuDeltaPtAfterAccCuts"), cent, nCh, pt, eta, theta, rAbsorb, dcaXY, pDca, chi2, dpt);
+        }
+      }
+    }
+
+    registry.fill(HIST("h3MultNchNmu"), cent, nCh, nMu);
+
+    // Fill number of muons of various types with multiplicity
+    for (auto indexType{0u}; indexType < NTrackTypes; ++indexType) {
+      if (nMuType[indexType] > 0) {
+        registry.fill(HIST("hMultNchNmuTrackType"), cent, nCh, nMuType[indexType], indexType);
+      }
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{adaptAnalysisTask<HfTaskSingleMuonMult>(cfgc)};
+}