diff --git a/PWGHF/HFL/CMakeLists.txt b/PWGHF/HFL/CMakeLists.txt index bbfd7adac2b..ab5ce2329af 100644 --- a/PWGHF/HFL/CMakeLists.txt +++ b/PWGHF/HFL/CMakeLists.txt @@ -1,10 +1,12 @@ -# 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. +#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". +#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 +#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. + +add_subdirectory(Tasks) diff --git a/PWGHF/HFL/Tasks/.HFm.cxx b/PWGHF/HFL/Tasks/.HFm.cxx new file mode 100644 index 00000000000..c3a736c097e --- /dev/null +++ b/PWGHF/HFL/Tasks/.HFm.cxx @@ -0,0 +1,346 @@ +// 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 HFmAnalysis.cxx +/// \brief Task derived from the DQ framework and used to extract the observables on single muons needed for the HF-muon analysis. +/// \author Maolin Zhang , CCNU + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/AnalysisCut.h" +#include "PWGDQ/Core/AnalysisCompositeCut.h" +#include "PWGDQ/Core/CutsLibrary.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "Framework/HistogramRegistry.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "ReconstructionDataFormats/TrackFwd.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; + +namespace o2::aod +{ +namespace dq_analysis_flags +{ +DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); +} +DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", dq_analysis_flags::IsEventSelected); +} // namespace o2::aod + +using MyCollisions = o2::soa::Join; +using MyMcCollisions = o2::soa::Join; +using MyEventsSelected = o2::soa::Join; +using MyMcEventsSelected = o2::soa::Join; +//using MyMuons = o2::soa::Join; +using MyMuons = o2::soa::Join; +//using MyMcMuons = o2::soa::Join; +using MyMcMuons = o2::soa::Join; + +constexpr static uint32_t gMuonFillMap(VarManager::ObjTypes::Muon); +constexpr static uint32_t gEventFillMap(VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision); +constexpr static uint32_t gTrackMCFillMap(VarManager::ObjTypes::ParticleMC); + +/// Defines the histograms for classes required in analysis. +void defineHistograms(HistogramManager*, TString); + +struct HfmEventSelection { + Configurable applySoftwareTrigger{"appllySoftwareTrigger", false, "whether to apply the software trigger"}; + Configurable softwareTrigger{"softwareTrigger", VarManager::kIsMuonSingleLowPt7, "software trigger flag"}; + Configurable applyCutZVtx{"applyCutZVtx", false, "whether to apply the VtxZ cut"}; + Configurable zVtxMin{"zVtxMin", -10., "min. z of primary vertex [cm]"}; + Configurable zVtxMax{"zVtxMax", 10., "max. z of primary vertex [cm]"}; + + HistogramManager* histMan; + OutputObj outputList{"output"}; + + float* values; + AnalysisCompositeCut* eventCut; + + void init(o2::framework::InitContext&) + { + VarManager::SetDefaultVarNames(); + values = new float[VarManager::kNVars]; + eventCut = new AnalysisCompositeCut("event selection", "Event Selection", true); + + AnalysisCut cut; + if (applyCutZVtx) + cut.AddCut(VarManager::kVtxZ, zVtxMin, zVtxMax, false); + if (applySoftwareTrigger) + cut.AddCut(softwareTrigger, 0.5, 1.5, false); + eventCut->AddCut(&cut); + + histMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + histMan->SetUseDefaultVariableNames(kTRUE); + histMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + + defineHistograms(histMan, "EventBeforeCuts;EventAfterCuts;"); + VarManager::SetUseVars(histMan->GetUsedVars()); + outputList.setObject(histMan->GetMainHistogramList()); + } + + Produces eventSel; + + template + void runEventSel(TEvent const& event, aod::BCs const& bcs) + { + // select events + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event, values); + histMan->FillHistClass("EventBeforeCuts", values); + + if (eventCut->IsSelected(values)) { + histMan->FillHistClass("EventAfterCuts", values); + eventSel(1); + } else { + eventSel(0); + } + } + + void processEvent(MyCollisions::iterator const& event, aod::BCs const& bcs) + { + runEventSel(event, bcs); + } + + void processEventMC(MyMcCollisions::iterator const& event, aod::BCs const& bcs) + { + runEventSel(event, bcs); + } + + PROCESS_SWITCH(HfmEventSelection, processEvent, "run event selection with real data", true); + PROCESS_SWITCH(HfmEventSelection, processEventMC, "run event selection with MC data", false); +}; + +struct MuonSelection { + Configurable muonCuts{"muonCuts", "muonQualityCuts", "muon selection"}; + + float* values; + AnalysisCompositeCut* trackCut; + + o2::framework::HistogramRegistry registry{ + "registry", + {}, + OutputObjHandlingPolicy::AnalysisObject, + true, + true}; + + void init(o2::framework::InitContext&) + { + AxisSpec axispT{200, 0., 200., "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec axisEta{500, -5., 5., "#eta"}; + AxisSpec axisDCA{500, 0., 5., "DCA (cm)"}; + AxisSpec axisSign{5, -2.5, 2.5, "Charge"}; + AxisSpec axisP{500, 0., 500., "p (GeV/#it{c})"}; + AxisSpec axisVtxZ{400, -20., 20., "Vertex Z (cm)"}; + + AxisSpec axisTrkType{5, -0.5, 4.5, "TrackType"}; + AxisSpec axisMCmask{200, -0.5, 199.5, "mcMask"}; + // kinematics for MC + AxisSpec axispTGen{200, 0., 200., "#it{p}_{T} Truth (GeV/#it{c})"}; + AxisSpec axisEtaGen{500, -5., 5., "#eta_{Truth}"}; + AxisSpec axisPGen{500, 0., 500., "p_{Truth} (GeV/#it{c})"}; + AxisSpec axispTDif{200, -2., 2., "#it{p}_{T} diff (GeV/#it{c})"}; + AxisSpec axisEtaDif{200, -2., 2., "#eta diff"}; + AxisSpec axisPDif{200, -2., 2., "p diff (GeV/#it{c})"}; + + HistogramConfigSpec hTHnMu{HistType::kTHnSparseD, {axispT, axisEta, axisDCA, axisSign, axisP, axisVtxZ, axisTrkType, axisMCmask}, 8}; + HistogramConfigSpec hTHnPt{HistType::kTHnSparseD, {axispT, axispTGen, axispTDif, axisTrkType}, 4}; + HistogramConfigSpec hTHnEta{HistType::kTHnSparseD, {axisEta, axisEtaGen, axisEtaDif, axisTrkType}, 4}; + HistogramConfigSpec hTHnP{HistType::kTHnSparseD, {axisP, axisPGen, axisPDif, axisTrkType}, 4}; + + registry.add("hMuBcuts", "", hTHnMu); + registry.add("hMuAcuts", "", hTHnMu); + registry.add("hPTBcuts", "", hTHnPt); + registry.add("hPTAcuts", "", hTHnPt); + registry.add("hEtaBcuts", "", hTHnEta); + registry.add("hEtaAcuts", "", hTHnEta); + registry.add("hPBcuts", "", hTHnP); + registry.add("hPAcuts", "", hTHnP); + + VarManager::SetDefaultVarNames(); + values = new float[VarManager::kNVars]; + + trackCut = new AnalysisCompositeCut(true); + TString selectStr = muonCuts.value; + trackCut->AddCut(dqcuts::GetAnalysisCut(selectStr.Data())); + VarManager::SetUseVars(AnalysisCut::fgUsedVars); + } + + template + void runMuonSel(TEvent const& event, aod::BCs const& bcs, TMuons const& tracks) + { + // select muons in data + if (event.isEventSelected() == 0) + return; + + VarManager::ResetValues(0, VarManager::kNMuonTrackVariables, values); + VarManager::FillEvent(event, values); + + // loop over muon tracks + for (auto const& track : tracks) { + VarManager::FillTrack(track, values); + + // compute DCAXY + /* std::vector vecCovs{ + track.cXX(), track.cXY(), track.cYY(), + track.cPhiX(), track.cPhiY(), track.cPhiPhi(), + track.cTglX(), track.cTglY(), track.cTglPhi(), track.cTglTgl(), + track.c1PtX(), track.c1PtY(), track.c1PtPhi(), track.c1PtTgl(), track.c1Pt21Pt2()}; + + SMatrix55 trackCovs(vecCovs.begin(), vecCovs.end()); + SMatrix5 trackPars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + + o2::track::TrackParCovFwd trackParcov{track.z(), trackPars, trackCovs, track.chi2()}; + trackParcov.propagateToZlinear(event.posZ()); + + const auto dcaX(trackParcov.getX() - event.posX()); + const auto dcaY(trackParcov.getY() - event.posY());*/ + const auto dcaXY(std::sqrt(values[VarManager::kMuonDCAx] * values[VarManager::kMuonDCAx] + values[VarManager::kMuonDCAy] * values[VarManager::kMuonDCAy])); + + // Before Muon Cuts + registry.fill(HIST("hMuBcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], 0); + // After Muon Cuts + if (trackCut->IsSelected(values)) + registry.fill(HIST("hMuAcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], 0); + } // end loop over muon tracks + } + + template + void runMuonSelMC(TEvent const& event, aod::BCs const& bcs, TMuons const& tracks, TMC const& mc) + { + // select muons in MC + if (event.isEventSelected() == 0) + return; + + VarManager::ResetValues(0, VarManager::kNMuonTrackVariables, values); + VarManager::FillEvent(event, values); + + // loop over muon tracks + for (auto const& track : tracks) { + VarManager::FillTrack(track, values); + + if (!track.has_mcParticle()) { + continue; + } + auto mcParticle = track.mcParticle(); + + VarManager::FillTrack(mcParticle, values); + + // compute DCAXY + /* std::vector vecCovs{ + track.cXX(), track.cXY(), track.cYY(), + track.cPhiX(), track.cPhiY(), track.cPhiPhi(), + track.cTglX(), track.cTglY(), track.cTglPhi(), track.cTglTgl(), + track.c1PtX(), track.c1PtY(), track.c1PtPhi(), track.c1PtTgl(), track.c1Pt21Pt2()}; + + SMatrix55 trackCovs(vecCovs.begin(), vecCovs.end()); + SMatrix5 trackPars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + + o2::track::TrackParCovFwd trackParcov{track.z(), trackPars, trackCovs, track.chi2()}; + trackParcov.propagateToZlinear(event.posZ()); + + const auto dcaX(trackParcov.getX() - event.posX()); + const auto dcaY(trackParcov.getY() - event.posY()); + const auto dcaXY(std::sqrt(dcaX * dcaX + dcaY * dcaY));*/ + + const auto dcaXY(std::sqrt(values[VarManager::kMuonDCAx] * values[VarManager::kMuonDCAx] + values[VarManager::kMuonDCAy] * values[VarManager::kMuonDCAy])); + // Before Muon Cuts + registry.fill(HIST("hMuBcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], + track.mcMask()); + + registry.fill(HIST("hPTBcuts"), + values[VarManager::kPt], values[VarManager::kMCPt], values[VarManager::kMCPt] - values[VarManager::kPt], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hEtaBcuts"), + values[VarManager::kEta], values[VarManager::kMCEta], values[VarManager::kMCEta] - values[VarManager::kEta], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hPBcuts"), + values[VarManager::kP], values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]), values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]) - values[VarManager::kP], values[VarManager::kMuonTrackType]); + // After Muon Cuts + if (trackCut->IsSelected(values)) { + registry.fill(HIST("hMuAcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], + track.mcMask()); + + registry.fill(HIST("hPTAcuts"), + values[VarManager::kPt], values[VarManager::kMCPt], values[VarManager::kMCPt] - values[VarManager::kPt], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hEtaAcuts"), + values[VarManager::kEta], values[VarManager::kMCEta], values[VarManager::kMCEta] - values[VarManager::kEta], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hPAcuts"), + values[VarManager::kP], values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]), values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]) - values[VarManager::kP], values[VarManager::kMuonTrackType]); + } + } // end loop over muon tracks + } + + void processMuon(MyEventsSelected::iterator const& event, aod::BCs const& bcs, + MyMuons const& tracks) + { + runMuonSel(event, bcs, tracks); + } + + void processMuonMC(MyMcEventsSelected::iterator const& event, aod::BCs const& bcs, + MyMcMuons const& tracks, aod::McParticles const& mc) + { + runMuonSelMC(event, bcs, tracks, mc); + } + + PROCESS_SWITCH(MuonSelection, processMuon, "run muon selection with real data", true); + PROCESS_SWITCH(MuonSelection, processMuonMC, "run muon selection with MC data", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + }; +} + +void defineHistograms(HistogramManager* histMan, TString histClasses) +{ + // + // Define here the histograms for all the classes required in analysis. + // The histogram classes are provided in the histClasses string, separated by semicolon ";" + // The histogram classes and their components histograms are defined below depending on the name of the histogram class + // + std::unique_ptr objArray(histClasses.Tokenize(";")); + + for (Int_t iclass = 0; iclass < objArray->GetEntries(); ++iclass) { + TString classStr = objArray->At(iclass)->GetName(); + histMan->AddHistClass(classStr.Data()); + + if (classStr.Contains("Event")) + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", "trigger,all"); + } +} diff --git a/PWGHF/HFL/Tasks/.HFmAnalysis.cxx b/PWGHF/HFL/Tasks/.HFmAnalysis.cxx new file mode 100644 index 00000000000..8a1ea382727 --- /dev/null +++ b/PWGHF/HFL/Tasks/.HFmAnalysis.cxx @@ -0,0 +1,339 @@ +// 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 HFmAnalysis.cxx +/// \brief Task derived from the DQ framework and used to extract the observables on single muons needed for the HF-muon analysis. +/// \author Maolin Zhang , CCNU + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/AnalysisCut.h" +#include "PWGDQ/Core/AnalysisCompositeCut.h" +#include "PWGDQ/Core/CutsLibrary.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "Framework/HistogramRegistry.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "ReconstructionDataFormats/TrackFwd.h" +#include "Common/DataModel/EventSelection.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; + +namespace o2::aod +{ +namespace dq_analysis_flags +{ +DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); +} +DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", dq_analysis_flags::IsEventSelected); +} // namespace o2::aod + +using MyCollisions = o2::soa::Join; +using MyMcCollisions = o2::soa::Join; +using MyEventsSelected = o2::soa::Join; +using MyMcEventsSelected = o2::soa::Join; +using MyMuons = o2::soa::Join; + +constexpr static uint32_t gMuonFillMap(VarManager::ObjTypes::Muon | VarManager::ObjTypes::MuonCov); +constexpr static uint32_t gEventFillMap(VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision); +constexpr static uint32_t gTrackMCFillMap(VarManager::ObjTypes::ParticleMC); + +/// Defines the histograms for classes required in analysis. +void defineHistograms(HistogramManager*, TString); + +struct HfmEventSelection { + Configurable applySoftwareTrigger{"appllySoftwareTrigger", false, "whether to apply the software trigger"}; + Configurable softwareTrigger{"softwareTrigger", VarManager::kIsMuonSingleLowPt7, "software trigger flag"}; + Configurable applyCutZVtx{"applyCutZVtx", false, "whether to apply the VtxZ cut"}; + Configurable zVtxMin{"zVtxMin", -10., "min. z of primary vertex [cm]"}; + Configurable zVtxMax{"zVtxMax", 10., "max. z of primary vertex [cm]"}; + + HistogramManager* histMan; + OutputObj outputList{"output"}; + + float* values; + AnalysisCompositeCut* eventCut; + + void init(o2::framework::InitContext&) + { + VarManager::SetDefaultVarNames(); + values = new float[VarManager::kNVars]; + eventCut = new AnalysisCompositeCut("event selection", "Event Selection", true); + + AnalysisCut cut; + if (applyCutZVtx) + cut.AddCut(VarManager::kVtxZ, zVtxMin, zVtxMax, false); + if (applySoftwareTrigger) + cut.AddCut(softwareTrigger, 0.5, 1.5, false); + eventCut->AddCut(&cut); + + histMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + histMan->SetUseDefaultVariableNames(kTRUE); + histMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + + defineHistograms(histMan, "EventBeforeCuts;EventAfterCuts;"); + VarManager::SetUseVars(histMan->GetUsedVars()); + outputList.setObject(histMan->GetMainHistogramList()); + } + + Produces eventSel; + + template + void runEventSel(TEvent const& event, aod::BCs const& bcs) + { + // select events + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event, values); + histMan->FillHistClass("EventBeforeCuts", values); + + if (eventCut->IsSelected(values)) { + histMan->FillHistClass("EventAfterCuts", values); + eventSel(1); + } else { + eventSel(0); + } + } + + void processEvent(MyCollisions::iterator const& event, aod::BCs const& bcs) + { + runEventSel(event, bcs); + } + + void processEventMC(MyMcCollisions::iterator const& event, aod::BCs const& bcs) + { + runEventSel(event, bcs); + } + + PROCESS_SWITCH(HfmEventSelection, processEvent, "run event selection with real data", true); + PROCESS_SWITCH(HfmEventSelection, processEventMC, "run event selection with MC data", false); +}; + +struct MuonSelection { + Configurable muonCuts{"muonCuts", "muonQualityCuts", "muon selection"}; + + float* values; + AnalysisCompositeCut* trackCut; + + o2::framework::HistogramRegistry registry{ + "registry", + {}, + OutputObjHandlingPolicy::AnalysisObject, + true, + true}; + + void init(o2::framework::InitContext&) + { + AxisSpec axispT{200, 0., 200., "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec axisEta{500, -5., 5., "#eta"}; + AxisSpec axisDCA{500, 0., 5., "DCA (cm)"}; + AxisSpec axisSign{5, -2.5, 2.5, "Charge"}; + AxisSpec axisP{500, 0., 500., "p (GeV/#it{c})"}; + AxisSpec axisVtxZ{400, -20., 20., "Vertex Z (cm)"}; + + AxisSpec axisTrkType{5, -0.5, 4.5, "TrackType"}; + AxisSpec axisMCmask{200, -0.5, 199.5, "mcMask"}; + // kinematics for MC + AxisSpec axispTGen{200, 0., 200., "#it{p}_{T} Truth (GeV/#it{c})"}; + AxisSpec axisEtaGen{500, -5., 5., "#eta_{Truth}"}; + AxisSpec axisPGen{500, 0., 500., "p_{Truth} (GeV/#it{c})"}; + AxisSpec axispTDif{200, -2., 2., "#it{p}_{T} diff (GeV/#it{c})"}; + AxisSpec axisEtaDif{200, -2., 2., "#eta diff"}; + AxisSpec axisPDif{200, -2., 2., "p diff (GeV/#it{c})"}; + + HistogramConfigSpec hTHnMu{HistType::kTHnSparseD, {axispT, axisEta, axisDCA, axisSign, axisP, axisVtxZ, axisTrkType, axisMCmask}, 8}; + HistogramConfigSpec hTHnPt{HistType::kTHnSparseD, {axispT, axispTGen, axispTDif, axisTrkType}, 4}; + HistogramConfigSpec hTHnEta{HistType::kTHnSparseD, {axisEta, axisEtaGen, axisEtaDif, axisTrkType}, 4}; + HistogramConfigSpec hTHnP{HistType::kTHnSparseD, {axisP, axisPGen, axisPDif, axisTrkType}, 4}; + + registry.add("hMuBcuts", "", hTHnMu); + registry.add("hMuAcuts", "", hTHnMu); + registry.add("hPTBcuts", "", hTHnPt); + registry.add("hPTAcuts", "", hTHnPt); + registry.add("hEtaBcuts", "", hTHnEta); + registry.add("hEtaAcuts", "", hTHnEta); + registry.add("hPBcuts", "", hTHnP); + registry.add("hPAcuts", "", hTHnP); + + VarManager::SetDefaultVarNames(); + values = new float[VarManager::kNVars]; + + trackCut = new AnalysisCompositeCut(true); + TString selectStr = muonCuts.value; + trackCut->AddCut(dqcuts::GetAnalysisCut(selectStr.Data())); + VarManager::SetUseVars(AnalysisCut::fgUsedVars); + } + + template + void runMuonSel(TEvent const& event, aod::BCs const& bcs, TMuons const& tracks) + { + // select muons in data + if (event.isEventSelected() == 0) + return; + + VarManager::ResetValues(0, VarManager::kNMuonTrackVariables, values); + VarManager::FillEvent(event, values); + + // loop over muon tracks + for (auto const& track : tracks) { + VarManager::FillTrack(track, values); + + // compute DCAXY + std::vector vecCovs{ + track.cXX(), track.cXY(), track.cYY(), + track.cPhiX(), track.cPhiY(), track.cPhiPhi(), + track.cTglX(), track.cTglY(), track.cTglPhi(), track.cTglTgl(), + track.c1PtX(), track.c1PtY(), track.c1PtPhi(), track.c1PtTgl(), track.c1Pt21Pt2()}; + + SMatrix55 trackCovs(vecCovs.begin(), vecCovs.end()); + SMatrix5 trackPars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + + o2::track::TrackParCovFwd trackParcov{track.z(), trackPars, trackCovs, track.chi2()}; + trackParcov.propagateToZlinear(event.posZ()); + + const auto dcaX(trackParcov.getX() - event.posX()); + const auto dcaY(trackParcov.getY() - event.posY()); + const auto dcaXY(std::sqrt(dcaX * dcaX + dcaY * dcaY)); + // Before Muon Cuts + registry.fill(HIST("hMuBcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], 0); + // After Muon Cuts + if (trackCut->IsSelected(values)) + registry.fill(HIST("hMuAcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], 0); + } // end loop over muon tracks + } + + template + void runMuonSelMC(TEvent const& event, aod::BCs const& bcs, TMuons const& tracks, TMC const& mc) + { + // select muons in MC + if (event.isEventSelected() == 0) + return; + + VarManager::ResetValues(0, VarManager::kNMuonTrackVariables, values); + VarManager::FillEvent(event, values); + + // loop over muon tracks + for (auto const& track : tracks) { + VarManager::FillTrack(track, values); + + if (!track.has_mcParticle()) { + continue; + } + auto mcParticle = track.mcParticle(); + + VarManager::FillTrack(mcParticle, values); + + // compute DCAXY + std::vector vecCovs{ + track.cXX(), track.cXY(), track.cYY(), + track.cPhiX(), track.cPhiY(), track.cPhiPhi(), + track.cTglX(), track.cTglY(), track.cTglPhi(), track.cTglTgl(), + track.c1PtX(), track.c1PtY(), track.c1PtPhi(), track.c1PtTgl(), track.c1Pt21Pt2()}; + + SMatrix55 trackCovs(vecCovs.begin(), vecCovs.end()); + SMatrix5 trackPars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + + o2::track::TrackParCovFwd trackParcov{track.z(), trackPars, trackCovs, track.chi2()}; + trackParcov.propagateToZlinear(event.posZ()); + + const auto dcaX(trackParcov.getX() - event.posX()); + const auto dcaY(trackParcov.getY() - event.posY()); + const auto dcaXY(std::sqrt(dcaX * dcaX + dcaY * dcaY)); + // Before Muon Cuts + registry.fill(HIST("hMuBcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], + track.mcMask()); + + registry.fill(HIST("hPTBcuts"), + values[VarManager::kPt], values[VarManager::kMCPt], values[VarManager::kMCPt] - values[VarManager::kPt], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hEtaBcuts"), + values[VarManager::kEta], values[VarManager::kMCEta], values[VarManager::kMCEta] - values[VarManager::kEta], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hPBcuts"), + values[VarManager::kP], values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]), values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]) - values[VarManager::kP], values[VarManager::kMuonTrackType]); + // After Muon Cuts + if (trackCut->IsSelected(values)) { + registry.fill(HIST("hMuAcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], + track.mcMask()); + + registry.fill(HIST("hPTAcuts"), + values[VarManager::kPt], values[VarManager::kMCPt], values[VarManager::kMCPt] - values[VarManager::kPt], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hEtaAcuts"), + values[VarManager::kEta], values[VarManager::kMCEta], values[VarManager::kMCEta] - values[VarManager::kEta], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hPAcuts"), + values[VarManager::kP], values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]), values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]) - values[VarManager::kP], values[VarManager::kMuonTrackType]); + } + } // end loop over muon tracks + } + + void processMuon(MyEventsSelected::iterator const& event, aod::BCs const& bcs, + aod::FullFwdTracks const& tracks) + { + runMuonSel(event, bcs, tracks); + } + + void processMuonMC(MyMcEventsSelected::iterator const& event, aod::BCs const& bcs, + MyMuons const& tracks, aod::McParticles const& mc) + { + runMuonSelMC(event, bcs, tracks, mc); + } + + PROCESS_SWITCH(MuonSelection, processMuon, "run muon selection with real data", true); + PROCESS_SWITCH(MuonSelection, processMuonMC, "run muon selection with MC data", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + }; +} + +void defineHistograms(HistogramManager* histMan, TString histClasses) +{ + // + // Define here the histograms for all the classes required in analysis. + // The histogram classes are provided in the histClasses string, separated by semicolon ";" + // The histogram classes and their components histograms are defined below depending on the name of the histogram class + // + std::unique_ptr objArray(histClasses.Tokenize(";")); + + for (Int_t iclass = 0; iclass < objArray->GetEntries(); ++iclass) { + TString classStr = objArray->At(iclass)->GetName(); + histMan->AddHistClass(classStr.Data()); + + if (classStr.Contains("Event")) + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", "trigger,all"); + } +} diff --git a/PWGHF/HFL/Tasks/.nfs0000001807b50009000001d5 b/PWGHF/HFL/Tasks/.nfs0000001807b50009000001d5 new file mode 100644 index 00000000000..9f7e2ee7465 Binary files /dev/null and b/PWGHF/HFL/Tasks/.nfs0000001807b50009000001d5 differ diff --git a/PWGHF/HFL/Tasks/CMakeLists.txt b/PWGHF/HFL/Tasks/CMakeLists.txt new file mode 100644 index 00000000000..0d099db11ff --- /dev/null +++ b/PWGHF/HFL/Tasks/CMakeLists.txt @@ -0,0 +1,15 @@ +#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. + +o2physics_add_dpl_workflow(muon - analysis + SOURCES HFmAnalysis.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) diff --git a/PWGHF/HFL/Tasks/HFmAnalysis.cxx b/PWGHF/HFL/Tasks/HFmAnalysis.cxx new file mode 100644 index 00000000000..7ea5efdd197 --- /dev/null +++ b/PWGHF/HFL/Tasks/HFmAnalysis.cxx @@ -0,0 +1,317 @@ +// 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 HFmAnalysis.cxx +/// \brief Task derived from the DQ framework and used to extract the observables on single muons needed for the HF-muon analysis. +/// \author Maolin Zhang , CCNU + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/AnalysisCut.h" +#include "PWGDQ/Core/AnalysisCompositeCut.h" +#include "PWGDQ/Core/CutsLibrary.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "Framework/HistogramRegistry.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "ReconstructionDataFormats/TrackFwd.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; + +namespace o2::aod +{ +namespace dq_analysis_flags +{ +DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); +} +DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", dq_analysis_flags::IsEventSelected); +} // namespace o2::aod + +using MyCollisions = o2::soa::Join; +using MyMcCollisions = o2::soa::Join; +using MyEventsSelected = o2::soa::Join; +using MyMcEventsSelected = o2::soa::Join; +using MyMuons = o2::soa::Join; +using MyMcMuons = o2::soa::Join; + +constexpr static uint32_t gMuonFillMap(VarManager::ObjTypes::Muon); +constexpr static uint32_t gEventFillMap(VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision); +constexpr static uint32_t gTrackMCFillMap(VarManager::ObjTypes::ParticleMC); + +/// Defines the histograms for classes required in analysis. +void defineHistograms(HistogramManager*, TString); + +struct HfmEventSelection { + Configurable applySoftwareTrigger{"appllySoftwareTrigger", false, "whether to apply the software trigger"}; + Configurable softwareTrigger{"softwareTrigger", VarManager::kIsMuonSingleLowPt7, "software trigger flag"}; + Configurable applyCutZVtx{"applyCutZVtx", false, "whether to apply the VtxZ cut"}; + Configurable zVtxMin{"zVtxMin", -10., "min. z of primary vertex [cm]"}; + Configurable zVtxMax{"zVtxMax", 10., "max. z of primary vertex [cm]"}; + + HistogramManager* histMan; + OutputObj outputList{"output"}; + + float* values; + AnalysisCompositeCut* eventCut; + + void init(o2::framework::InitContext&) + { + VarManager::SetDefaultVarNames(); + values = new float[VarManager::kNVars]; + eventCut = new AnalysisCompositeCut("event selection", "Event Selection", true); + + AnalysisCut cut; + if (applyCutZVtx) + cut.AddCut(VarManager::kVtxZ, zVtxMin, zVtxMax, false); + if (applySoftwareTrigger) + cut.AddCut(softwareTrigger, 0.5, 1.5, false); + eventCut->AddCut(&cut); + + histMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + histMan->SetUseDefaultVariableNames(kTRUE); + histMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + + defineHistograms(histMan, "EventBeforeCuts;EventAfterCuts;"); + VarManager::SetUseVars(histMan->GetUsedVars()); + outputList.setObject(histMan->GetMainHistogramList()); + } + + Produces eventSel; + + template + void runEventSel(TEvent const& event, aod::BCs const& bcs) + { + // select events + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event, values); + histMan->FillHistClass("EventBeforeCuts", values); + + if (eventCut->IsSelected(values)) { + histMan->FillHistClass("EventAfterCuts", values); + eventSel(1); + } else { + eventSel(0); + } + } + + void processEvent(MyCollisions::iterator const& event, aod::BCs const& bcs) + { + runEventSel(event, bcs); + } + + void processEventMC(MyMcCollisions::iterator const& event, aod::BCs const& bcs) + { + runEventSel(event, bcs); + } + + PROCESS_SWITCH(HfmEventSelection, processEvent, "run event selection with real data", true); + PROCESS_SWITCH(HfmEventSelection, processEventMC, "run event selection with MC data", false); +}; + +struct MuonSelection { + Configurable muonCuts{"muonCuts", "muonQualityCuts", "muon selection"}; + + float* values; + AnalysisCompositeCut* trackCut; + + o2::framework::HistogramRegistry registry{ + "registry", + {}, + OutputObjHandlingPolicy::AnalysisObject, + true, + true}; + + void init(o2::framework::InitContext&) + { + AxisSpec axispT{200, 0., 200., "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec axisEta{500, -5., 5., "#eta"}; + AxisSpec axisDCA{500, 0., 5., "DCA (cm)"}; + AxisSpec axisSign{5, -2.5, 2.5, "Charge"}; + AxisSpec axisP{500, 0., 500., "p (GeV/#it{c})"}; + AxisSpec axisVtxZ{400, -20., 20., "Vertex Z (cm)"}; + + AxisSpec axisTrkType{5, -0.5, 4.5, "TrackType"}; + AxisSpec axisMCmask{200, -0.5, 199.5, "mcMask"}; + // kinematics for MC + AxisSpec axispTGen{200, 0., 200., "#it{p}_{T} Truth (GeV/#it{c})"}; + AxisSpec axisEtaGen{500, -5., 5., "#eta_{Truth}"}; + AxisSpec axisPGen{500, 0., 500., "p_{Truth} (GeV/#it{c})"}; + AxisSpec axispTDif{200, -2., 2., "#it{p}_{T} diff (GeV/#it{c})"}; + AxisSpec axisEtaDif{200, -2., 2., "#eta diff"}; + AxisSpec axisPDif{200, -2., 2., "p diff (GeV/#it{c})"}; + + HistogramConfigSpec hTHnMu{HistType::kTHnSparseD, {axispT, axisEta, axisDCA, axisSign, axisP, axisVtxZ, axisTrkType, axisMCmask}, 8}; + HistogramConfigSpec hTHnPt{HistType::kTHnSparseD, {axispT, axispTGen, axispTDif, axisTrkType}, 4}; + HistogramConfigSpec hTHnEta{HistType::kTHnSparseD, {axisEta, axisEtaGen, axisEtaDif, axisTrkType}, 4}; + HistogramConfigSpec hTHnP{HistType::kTHnSparseD, {axisP, axisPGen, axisPDif, axisTrkType}, 4}; + + registry.add("hMuBcuts", "", hTHnMu); + registry.add("hMuAcuts", "", hTHnMu); + registry.add("hPTBcuts", "", hTHnPt); + registry.add("hPTAcuts", "", hTHnPt); + registry.add("hEtaBcuts", "", hTHnEta); + registry.add("hEtaAcuts", "", hTHnEta); + registry.add("hPBcuts", "", hTHnP); + registry.add("hPAcuts", "", hTHnP); + + VarManager::SetDefaultVarNames(); + values = new float[VarManager::kNVars]; + + trackCut = new AnalysisCompositeCut(true); + TString selectStr = muonCuts.value; + trackCut->AddCut(dqcuts::GetAnalysisCut(selectStr.Data())); + VarManager::SetUseVars(AnalysisCut::fgUsedVars); + } + + template + void runMuonSel(TEvent const& event, aod::BCs const& bcs, TMuons const& tracks) + { + // select muons in data + if (event.isEventSelected() == 0) + return; + + VarManager::ResetValues(0, VarManager::kNMuonTrackVariables, values); + VarManager::FillEvent(event, values); + + // loop over muon tracks + for (auto const& track : tracks) { + VarManager::FillTrack(track, values); + + // compute DCAXY + // With DQ's muon DCA calculation applied, we no longer need calculate it ourselve. + // Same for MC + const auto dcaXY(std::sqrt(values[VarManager::kMuonDCAx] * values[VarManager::kMuonDCAx] + values[VarManager::kMuonDCAy] * values[VarManager::kMuonDCAy])); + + // Before Muon Cuts + registry.fill(HIST("hMuBcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], 0); + // After Muon Cuts + if (trackCut->IsSelected(values)) + registry.fill(HIST("hMuAcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], 0); + } // end loop over muon tracks + } + + template + void runMuonSelMC(TEvent const& event, aod::BCs const& bcs, TMuons const& tracks, TMC const& mc) + { + // select muons in MC + if (event.isEventSelected() == 0) + return; + + VarManager::ResetValues(0, VarManager::kNMuonTrackVariables, values); + VarManager::FillEvent(event, values); + + // loop over muon tracks + for (auto const& track : tracks) { + VarManager::FillTrack(track, values); + + if (!track.has_mcParticle()) { + continue; + } + auto mcParticle = track.mcParticle(); + + VarManager::FillTrack(mcParticle, values); + + // compute DCAXY + const auto dcaXY(std::sqrt(values[VarManager::kMuonDCAx] * values[VarManager::kMuonDCAx] + values[VarManager::kMuonDCAy] * values[VarManager::kMuonDCAy])); + + // Before Muon Cuts + registry.fill(HIST("hMuBcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], + track.mcMask()); + + registry.fill(HIST("hPTBcuts"), + values[VarManager::kPt], values[VarManager::kMCPt], values[VarManager::kMCPt] - values[VarManager::kPt], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hEtaBcuts"), + values[VarManager::kEta], values[VarManager::kMCEta], values[VarManager::kMCEta] - values[VarManager::kEta], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hPBcuts"), + values[VarManager::kP], values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]), values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]) - values[VarManager::kP], values[VarManager::kMuonTrackType]); + // After Muon Cuts + if (trackCut->IsSelected(values)) { + registry.fill(HIST("hMuAcuts"), + values[VarManager::kPt], + values[VarManager::kEta], dcaXY, + values[VarManager::kCharge], track.p(), + values[VarManager::kVtxZ], + values[VarManager::kMuonTrackType], + track.mcMask()); + + registry.fill(HIST("hPTAcuts"), + values[VarManager::kPt], values[VarManager::kMCPt], values[VarManager::kMCPt] - values[VarManager::kPt], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hEtaAcuts"), + values[VarManager::kEta], values[VarManager::kMCEta], values[VarManager::kMCEta] - values[VarManager::kEta], values[VarManager::kMuonTrackType]); + registry.fill(HIST("hPAcuts"), + values[VarManager::kP], values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]), values[VarManager::kMCPt] * std::cosh(values[VarManager::kMCEta]) - values[VarManager::kP], values[VarManager::kMuonTrackType]); + } + } // end loop over muon tracks + } + + void processMuon(MyEventsSelected::iterator const& event, aod::BCs const& bcs, + MyMuons const& tracks) + { + runMuonSel(event, bcs, tracks); + } + + void processMuonMC(MyMcEventsSelected::iterator const& event, aod::BCs const& bcs, + MyMcMuons const& tracks, aod::McParticles const& mc) + { + runMuonSelMC(event, bcs, tracks, mc); + } + + PROCESS_SWITCH(MuonSelection, processMuon, "run muon selection with real data", true); + PROCESS_SWITCH(MuonSelection, processMuonMC, "run muon selection with MC data", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + }; +} + +void defineHistograms(HistogramManager* histMan, TString histClasses) +{ + // + // Define here the histograms for all the classes required in analysis. + // The histogram classes are provided in the histClasses string, separated by semicolon ";" + // The histogram classes and their components histograms are defined below depending on the name of the histogram class + // + std::unique_ptr objArray(histClasses.Tokenize(";")); + + for (Int_t iclass = 0; iclass < objArray->GetEntries(); ++iclass) { + TString classStr = objArray->At(iclass)->GetName(); + histMan->AddHistClass(classStr.Data()); + + if (classStr.Contains("Event")) + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", "trigger,all"); + } +}