taskSingleMuonMult.cxx

// 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 <string>
#include <vector>
#include <array>

#include <TPDGCode.h>
#include <TString.h>

#include "TableHelper.h"

#include "CommonConstants/MathConstants.h"
#include "Framework/AnalysisDataModel.h"
#include "Framework/AnalysisTask.h"
#include "Framework/ASoAHelpers.h"
#include "Framework/HistogramRegistry.h"
#include "Framework/RunningWorkflowInfo.h"
#include "Framework/runDataProcessing.h"
#include "Framework/StaticFor.h"
#include "ReconstructionDataFormats/TrackFwd.h"

#include "Common/CCDB/EventSelectionParams.h"
#include "Common/Core/TrackSelection.h"
#include "Common/Core/trackUtilities.h"
#include "Common/DataModel/Centrality.h"
#include "Common/DataModel/EventSelection.h"
#include "Common/DataModel/FT0Corrected.h"
#include "Common/DataModel/Multiplicity.h"
#include "Common/DataModel/TrackSelectionTables.h"

#include "PWGLF/DataModel/LFResonanceTables.h"

using namespace o2;
using namespace o2::aod;
using namespace o2::framework;
using namespace o2::framework::expressions;
using namespace o2::aod::fwdtrack;

struct HfTaskSingleMuonMult {
  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::TracksExtra, aod::TracksIU, aod::TracksDCA, aod::TrackSelection>>;

  // Filter Global Track for Multiplicty
  Filter globalTrackFilter = ((o2::aod::track::isGlobalTrack == true) && (nabs(o2::aod::track::eta) < etaTrackMax) && ((o2::aod::track::pt) > ptTrackMin));

  HistogramRegistry registry{"registry"};

  void init(InitContext&)
  {
    AxisSpec axisCent = {101, -0.5, 100.5, "centrality"};
    AxisSpec axisEvent{3, 0.5, 3.5, "Event Selection"};
    AxisSpec axisEventSize{500, 0.5, 500.5, "Event Size"};
    AxisSpec axisVtxZ{80, -20., 20., "#it{z}_{vtx} (cm)"};
    AxisSpec axisMuTrk{5, 0.5, 5.5, "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("hEvent", " Number of Events", {HistType::kTH1F, {axisEvent}});
    registry.add("hEventSize", "Event Size", {HistType::kTH1F, {axisEventSize}});
    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("hMuTrkSel", "Selection of muon tracks at various kinematic cuts", {HistType::kTH1F, {axisMuTrk}});
    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("hEvent"));
    auto* xEv = hEvstat->GetXaxis();
    xEv->SetBinLabel(1, "All events");
    xEv->SetBinLabel(2, "sel8");
    xEv->SetBinLabel(3, "VtxZAfterSel");

    auto hMustat = registry.get<TH1>(HIST("hMuTrkSel"));
    auto* xMu = hMustat->GetXaxis();
    xMu->SetBinLabel(1, "noCut");
    xMu->SetBinLabel(2, "etaCut");
    xMu->SetBinLabel(3, "RAbsorbCut");
    xMu->SetBinLabel(4, "pDcaCut");
    xMu->SetBinLabel(5, "chi2Cut");

    // Number the types of muon tracks
    constexpr uint8_t nTrackTypes{static_cast<uint8_t>(ForwardTrackTypeEnum::MCHStandaloneTrack)};
  }

  void process(MyCollisions::iterator const& collision,
               MyTracks const& tracks,
               MyMuons const& muons)
  {
    registry.fill(HIST("hEvent"), 1);

    if (!collision.sel8()) {
      return;
    }
    registry.fill(HIST("hEvent"), 2);
    registry.fill(HIST("hVtxZBeforeSel"), collision.posZ());

    if (std::abs(collision.posZ()) > zVtxMax) {
      return;
    }
    registry.fill(HIST("hEvent"), 3);
    registry.fill(HIST("hVtxZAfterSel"), collision.posZ());

    // T0M centrality
    const auto cent = collision.centFT0M();
    registry.fill(HIST("hCentrality"), cent);

    // Charged particles
    std::size_t nCh{tracks.size()};
    if (nCh < 1) {
      return;
    }
    registry.fill(HIST("hEventSize"), 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 + 1] = {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("hMuTrkSel"), 1);
      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("hMuTrkSel"), 2);
      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("hMuTrkSel"), 3);
      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("hMuTrkSel"), 4);
      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("hMuTrkSel"), 5);

      // 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
    static_for<0u, nTrackTypes>([&](auto i) {
      constexpr int iType{i.value};
      if (nMuType[iType] > 0) {
        registry.fill(std::string(HIST("hMultNchNmuTrackType")), cent, nCh, nMuType[iType], iType);
      }
    });
    chTracks.clear();
  }
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
  return WorkflowSpec{adaptAnalysisTask<HfTaskSingleMuonMult>(cfgc)};
}