diff --git a/PWGLF/DataModel/LFResonanceTables.h b/PWGLF/DataModel/LFResonanceTables.h new file mode 100644 index 00000000000..112b99ea42e --- /dev/null +++ b/PWGLF/DataModel/LFResonanceTables.h @@ -0,0 +1,139 @@ +// 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 LFResonanceTables.h +/// \brief Definitions of tables of resonance decay candidates +/// +/// Inspired by StrangenessTables.h, FemtoDerived.h +/// +/// \author Bong-Hwi Lim + +#ifndef O2_ANALYSIS_LFRESONANCETABLES_H_ +#define O2_ANALYSIS_LFRESONANCETABLES_H_ + +#include + +#include "Common/DataModel/PIDResponse.h" +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ +/// Resonance Collisions +namespace resocollision +{ +DECLARE_SOA_COLUMN(MultV0M, multV0M, float); //! V0M multiplicity +DECLARE_SOA_COLUMN(Sphericity, sphericity, float); //! Sphericity of the event +} // namespace resocollision +DECLARE_SOA_TABLE(ResoCollisions, "AOD", "RESOCOL", + o2::soa::Index<>, + o2::aod::collision::PosX, + o2::aod::collision::PosY, + o2::aod::collision::PosZ, + resocollision::MultV0M, + resocollision::Sphericity, + timestamp::Timestamp); +using ResoCollision = ResoCollisions::iterator; + +// Resonance Daughters +// inspired from PWGCF/DataModel/FemtoDerived.h +namespace resodaughter +{ +/// Distinuishes the different daughter types +enum DaughterType { + kTrack, //! Track + kV0, //! V0 + kCascade, //! Cascade + kNDaughterTypes //! Number of Daughter types +}; +static constexpr std::string_view DaughterTypeName[kNDaughterTypes] = {"Tracks", "V0", "Cascade"}; //! Naming of the different particle types + +enum PDGtype { + kPion = BIT(0), + kKaon = BIT(1), + kProton = BIT(2), + kHasTOF = BIT(6) // Save hasTOF info for TOF selection +}; + +#define requireTPCPIDCutInFilter(mask) ((aod::resodaughter::tpcPIDselectionFlag & (uint8_t)aod::resodaughter::mask) == (uint8_t)aod::resodaughter::mask) +#define requireTOFPIDCutInFilter(mask) (((aod::resodaughter::tofPIDselectionFlag & (uint8_t)aod::resodaughter::kHasTOF) != (uint8_t)aod::resodaughter::kHasTOF) || (((aod::resodaughter::tofPIDselectionFlag & (uint8_t)aod::resodaughter::mask) == (uint8_t)aod::resodaughter::mask) && ((aod::resodaughter::tofPIDselectionFlag & (uint8_t)aod::resodaughter::kHasTOF) == (uint8_t)aod::resodaughter::kHasTOF))) +#define requireTPCPIDPionCutInFilter() requireTPCPIDCutInFilter(PDGtype::kPion) +#define requireTPCPIDKaonCutInFilter() requireTPCPIDCutInFilter(PDGtype::kKaon) +#define requireTPCPIDProtonCutInFilter() requireTPCPIDCutInFilter(PDGtype::kProton) +#define requireTOFPIDPionCutInFilter() requireTOFPIDCutInFilter(PDGtype::kPion) +#define requireTOFPIDKaonCutInFilter() requireTOFPIDCutInFilter(PDGtype::kKaon) +#define requireTOFPIDProtonCutInFilter() requireTOFPIDCutInFilter(PDGtype::kProton) + +DECLARE_SOA_INDEX_COLUMN(ResoCollision, resoCollision); +DECLARE_SOA_COLUMN(Pt, pt, float); //! p_T (GeV/c) +DECLARE_SOA_COLUMN(Px, px, float); //! p_x (GeV/c) +DECLARE_SOA_COLUMN(Py, py, float); //! p_y (GeV/c) +DECLARE_SOA_COLUMN(Pz, pz, float); //! p_z (GeV/c) +DECLARE_SOA_COLUMN(Eta, eta, float); //! Eta +DECLARE_SOA_COLUMN(Phi, phi, float); //! Phi +DECLARE_SOA_COLUMN(PartType, partType, uint8_t); //! Type of the particle, according to resodaughter::ParticleType +DECLARE_SOA_COLUMN(TempFitVar, tempFitVar, float); //! Observable for the template fitting (Track: DCA_xy, V0: CPA) +DECLARE_SOA_COLUMN(Indices, indices, int[2]); //! Field for the track indices to remove auto-correlations +DECLARE_SOA_COLUMN(Sign, sign, int8_t); //! Sign of the track charge +DECLARE_SOA_COLUMN(TPCNClsCrossedRows, tpcNClsCrossedRows, uint8_t); //! Number of TPC crossed rows +DECLARE_SOA_COLUMN(TPCPIDselectionFlag, tpcPIDselectionFlag, uint8_t); //! TPC PID selection +DECLARE_SOA_COLUMN(TOFPIDselectionFlag, tofPIDselectionFlag, uint8_t); //! TOF PID selection +DECLARE_SOA_COLUMN(DaughDCA, daughDCA, float); //! DCA between daughters +DECLARE_SOA_COLUMN(MLambda, mLambda, float); //! The invariant mass of V0 candidate, assuming lambda +DECLARE_SOA_COLUMN(MAntiLambda, mAntiLambda, float); //! The invariant mass of V0 candidate, assuming antilambda +DECLARE_SOA_COLUMN(TransRadius, transRadius, float); //! Transverse radius of the decay vertex +DECLARE_SOA_COLUMN(DecayVtxX, decayVtxX, float); //! X position of the decay vertex +DECLARE_SOA_COLUMN(DecayVtxY, decayVtxY, float); //! Y position of the decay vertex +DECLARE_SOA_COLUMN(DecayVtxZ, decayVtxZ, float); //! Z position of the decay vertex +} // namespace resodaughter +DECLARE_SOA_TABLE(ResoDaughters, "AOD", "RESODAUGHTERS", + o2::soa::Index<>, + resodaughter::ResoCollisionId, + resodaughter::Pt, + resodaughter::Px, + resodaughter::Py, + resodaughter::Pz, + resodaughter::Eta, + resodaughter::Phi, + resodaughter::PartType, + resodaughter::TempFitVar, + resodaughter::Indices, + resodaughter::Sign, + resodaughter::TPCNClsCrossedRows, + o2::aod::track::DcaXY, + o2::aod::track::DcaZ, + o2::aod::track::X, + o2::aod::track::Alpha, + resodaughter::TPCPIDselectionFlag, + resodaughter::TOFPIDselectionFlag, + o2::aod::pidtpc::TPCNSigmaPi, + o2::aod::pidtpc::TPCNSigmaKa, + o2::aod::pidtpc::TPCNSigmaPr, + o2::aod::pidtof::TOFNSigmaPi, + o2::aod::pidtof::TOFNSigmaKa, + o2::aod::pidtof::TOFNSigmaPr, + resodaughter::DaughDCA, + resodaughter::MLambda, + resodaughter::MAntiLambda, + resodaughter::TransRadius, + resodaughter::DecayVtxX, + resodaughter::DecayVtxY, + resodaughter::DecayVtxZ); +using ResoDaughter = ResoDaughters::iterator; + +using Reso2TracksExt = soa::Join; +using Reso2TracksMC = soa::Join; +using Reso2TracksPID = soa::Join; +using Reso2TracksPIDExt = soa::Join; + +} // namespace o2::aod +#endif // O2_ANALYSIS_LFRESONANCETABLES_H_ \ No newline at end of file diff --git a/PWGLF/TableProducer/CMakeLists.txt b/PWGLF/TableProducer/CMakeLists.txt index 02b921804ef..11e3443808a 100644 --- a/PWGLF/TableProducer/CMakeLists.txt +++ b/PWGLF/TableProducer/CMakeLists.txt @@ -31,5 +31,10 @@ o2physics_add_dpl_workflow(cascadefinder o2physics_add_dpl_workflow(nucleustreecreator SOURCES LFTreeCreatorNuclei.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(reso2initializer + SOURCES LFResonanceInitializer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/PWGLF/TableProducer/LFResonanceInitializer.cxx b/PWGLF/TableProducer/LFResonanceInitializer.cxx new file mode 100644 index 00000000000..7f5456d3fc5 --- /dev/null +++ b/PWGLF/TableProducer/LFResonanceInitializer.cxx @@ -0,0 +1,216 @@ +// 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 LFResonanceInitializer.cxx +/// \brief Initializes variables for the resonance candidate producers +/// +/// +/// \author Bong-Hwi Lim + +#include "Common/DataModel/PIDResponse.h" +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/StrangenessTables.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "PWGLF/DataModel/LFResonanceTables.h" +#include "PWGLF/Utils/collisionCuts.h" +#include "ReconstructionDataFormats/Track.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; + +/// Initializer for the resonance candidate producers +struct reso2initializer { + + Produces resoCollisions; + Produces reso2tracks; + + // Configurables + Configurable ConfIsRun3{"ConfIsRun3", false, "Running on Pilot beam"}; // Choose if running on converted data or pilot beam + Configurable ConfStoreV0{"ConfStoreV0", true, "True: store V0s"}; + + /// Event cuts + o2::analysis::CollisonCuts colCuts; + Configurable ConfEvtZvtx{"ConfEvtZvtx", 10.f, "Evt sel: Max. z-Vertex (cm)"}; + Configurable ConfEvtTriggerCheck{"ConfEvtTriggerCheck", true, "Evt sel: check for trigger"}; + Configurable ConfEvtTriggerSel{"ConfEvtTriggerSel", kINT7, "Evt sel: trigger"}; + Configurable ConfEvtOfflineCheck{"ConfEvtOfflineCheck", false, "Evt sel: check for offline selection"}; + + // Pre-selection cuts + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable pidnSigmaPreSelectionCut{"pidnSigmaPreSelectionCut", 5.0f, "TPC and TOF PID cut (loose, improve performance)"}; + Configurable mincrossedrows{"mincrossedrows", 70, "min crossed rows"}; + Configurable isRun2{"isRun2", 0, "if Run2: demand TPC refit"}; + + /// DCA Selections for V0 + // DCAr to PV + Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 0.05, "Track DCAr cut to PV Maximum"}; + Configurable cMinV0PosDCArToPVcut{"cMinV0PosDCArToPVcut", 0.05f, "V0 Positive Track DCAr cut to PV Minimum"}; // Pre-selection + Configurable cMinV0NegDCArToPVcut{"cMinV0NegDCArToPVcut", 0.05f, "V0 Negative Track DCAr cut to PV Minimum"}; // Pre-selection + // DCAz to PV + Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; + Configurable cMinDCAzToPVcut{"cMinDCAzToPVcut", 0.0, "Track DCAz cut to PV Minimum"}; + + Configurable cMinV0Radius{"cMinV0Radius", 5.0, "Minimum V0 radius from PV"}; + Configurable cMaxV0Radius{"cMaxV0Radius", 200.0, "Maximum V0 radius from PV"}; + Configurable cMinV0CosPA{"cMinV0CosPA", 0.995, "Minimum V0 CosPA to PV"}; + + HistogramRegistry qaRegistry{"QAHistos", { + {"hGoodTrackIndices", "hGoodTrackIndices", {HistType::kTH1F, {{4, 0.0f, 4.0f}}}}, + {"hGoodV0Indices", "hGoodV0Indices", {HistType::kTH1F, {{5, 0.0f, 5.0f}}}}, + }, + OutputObjHandlingPolicy::QAObject}; + + // Pre-filters for efficient process + // Filter tofPIDFilter = aod::track::tofExpMom < 0.f || ((aod::track::tofExpMom > 0.f) && ((nabs(aod::pidtof::tofNSigmaPi) < pidnSigmaPreSelectionCut) || (nabs(aod::pidtof::tofNSigmaKa) < pidnSigmaPreSelectionCut) || (nabs(aod::pidtof::tofNSigmaPr) < pidnSigmaPreSelectionCut))); // TOF + Filter tpcPIDFilter = nabs(aod::pidtpc::tpcNSigmaPi) < pidnSigmaPreSelectionCut || nabs(aod::pidtpc::tpcNSigmaKa) < pidnSigmaPreSelectionCut || nabs(aod::pidtpc::tpcNSigmaPr) < pidnSigmaPreSelectionCut; // TPC + Filter trackFilter = nabs(aod::track::eta) < cfgCutEta; // Eta cut + Filter trackCutFilter = requireGlobalTrackInFilter(); // Global track cuts + + void init(InitContext&) + { + colCuts.setCuts(ConfEvtZvtx, ConfEvtTriggerCheck, ConfEvtTriggerSel, ConfEvtOfflineCheck, ConfIsRun3); + colCuts.init(&qaRegistry); + } + + void process(const soa::Join::iterator& collision, + soa::Filtered const& tracks, o2::aod::V0Datas const& V0s, aod::BCsWithTimestamps const&) + { + auto bc = collision.bc_as(); /// adding timestamp to access magnetic field later + // Default event selection + if (!colCuts.isSelected(collision)) + return; + colCuts.fillQA(collision); + + if (ConfIsRun3) { + resoCollisions(collision.posX(), collision.posY(), collision.posZ(), collision.multFT0M(), colCuts.computeSphericity(collision, tracks), bc.timestamp()); + } else { + resoCollisions(collision.posX(), collision.posY(), collision.posZ(), collision.multFV0M(), colCuts.computeSphericity(collision, tracks), bc.timestamp()); + } + + int childIDs[2] = {0, 0}; // these IDs are necessary to keep track of the children + // Loop over tracks + for (auto& track : tracks) { + // Tracks are already filtered by the pre-filters + qaRegistry.fill(HIST("hGoodTrackIndices"), 0.5); + + // Add PID selection criteria here + uint8_t tpcPIDselections = 0; + uint8_t tofPIDselections = 0; + // TPC PID + if (std::abs(track.tpcNSigmaPi()) < pidnSigmaPreSelectionCut) + tpcPIDselections |= aod::resodaughter::PDGtype::kPion; + if (std::abs(track.tpcNSigmaKa()) < pidnSigmaPreSelectionCut) + tpcPIDselections |= aod::resodaughter::PDGtype::kKaon; + if (std::abs(track.tpcNSigmaPr()) < pidnSigmaPreSelectionCut) + tpcPIDselections |= aod::resodaughter::PDGtype::kProton; + // TOF PID + if (track.hasTOF()) { + tofPIDselections |= aod::resodaughter::PDGtype::kHasTOF; + if (std::abs(track.tofNSigmaPi()) < pidnSigmaPreSelectionCut) + tofPIDselections |= aod::resodaughter::PDGtype::kPion; + if (std::abs(track.tofNSigmaKa()) < pidnSigmaPreSelectionCut) + tofPIDselections |= aod::resodaughter::PDGtype::kKaon; + if (std::abs(track.tofNSigmaPr()) < pidnSigmaPreSelectionCut) + tofPIDselections |= aod::resodaughter::PDGtype::kProton; + } + + reso2tracks(resoCollisions.lastIndex(), + track.pt(), + track.px(), + track.py(), + track.pz(), + track.eta(), + track.phi(), + aod::resodaughter::DaughterType::kTrack, + track.dcaXY(), + childIDs, + track.sign(), + (uint8_t)track.tpcNClsCrossedRows(), + track.dcaXY(), + track.dcaZ(), + track.x(), + track.alpha(), + tpcPIDselections, + tofPIDselections, + track.tpcNSigmaPi(), + track.tpcNSigmaKa(), + track.tpcNSigmaPr(), + track.tofNSigmaPi(), + track.tofNSigmaKa(), + track.tofNSigmaPr(), + 0, 0, 0, + 0, 0, 0, 0); + } + /// V0s + if (ConfStoreV0) { + for (auto& v0 : V0s) { + qaRegistry.fill(HIST("hGoodV0Indices"), 0.5); + auto postrack = v0.posTrack_as(); + auto negtrack = v0.negTrack_as(); + + if (postrack.tpcNClsCrossedRows() < mincrossedrows) + continue; + if (negtrack.tpcNClsCrossedRows() < mincrossedrows) + continue; + qaRegistry.fill(HIST("hGoodV0Indices"), 1.5); + + if (fabs(postrack.dcaXY()) < cMinV0PosDCArToPVcut) + continue; + if (fabs(negtrack.dcaXY()) < cMinV0NegDCArToPVcut) + continue; + qaRegistry.fill(HIST("hGoodV0Indices"), 2.5); + + if ((v0.v0radius() > cMaxV0Radius) || (v0.v0radius() < cMinV0Radius)) + continue; + qaRegistry.fill(HIST("hGoodV0Indices"), 3.5); + if (v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < cMinV0CosPA) + continue; + qaRegistry.fill(HIST("hGoodV0Indices"), 4.5); + childIDs[0] = v0.posTrackId(); + childIDs[1] = v0.negTrackId(); + reso2tracks(resoCollisions.lastIndex(), + v0.pt(), + v0.px(), + v0.py(), + v0.pz(), + v0.eta(), + v0.phi(), + aod::resodaughter::DaughterType::kV0, + v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()), + childIDs, + 0, + 0, + v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()), + 0, + v0.x(), 0, 0, 0, 0, 0, 0, 0, 0, 0, + v0.dcaV0daughters(), v0.mLambda(), v0.mAntiLambda(), + v0.v0radius(), v0.x(), v0.y(), v0.z()); + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"lf-reso2initializer"}), + }; +} \ No newline at end of file diff --git a/PWGLF/Tasks/CMakeLists.txt b/PWGLF/Tasks/CMakeLists.txt index 19856979a23..10c85e09fdd 100644 --- a/PWGLF/Tasks/CMakeLists.txt +++ b/PWGLF/Tasks/CMakeLists.txt @@ -97,4 +97,9 @@ o2physics_add_dpl_workflow(track-checks o2physics_add_dpl_workflow(rsnanalysis SOURCES rsnanalysis.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(phianalysis + SOURCES phianalysis.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/PWGLF/Tasks/phianalysis.cxx b/PWGLF/Tasks/phianalysis.cxx new file mode 100644 index 00000000000..08bb4c08464 --- /dev/null +++ b/PWGLF/Tasks/phianalysis.cxx @@ -0,0 +1,159 @@ +// 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 phianalysis.cxx +/// \brief Reconstruction of track-track decay resonance candidates +/// +/// +/// \author Bong-Hwi Lim + +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/runDataProcessing.h" +#include "PWGLF/DataModel/LFResonanceTables.h" +#include +#include "DataFormatsParameters/GRPObject.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; + +struct phianalysis { + framework::Service ccdb; /// Accessing the CCDB + ConfigurableAxis CfgMultBins{"CfgMultBins", {VARIABLE_WIDTH, 0.0f, 20.0f, 40.0f, 60.0f, 80.0f, 100.0f, 200.0f, 99999.f}, "Mixing bins - multiplicity"}; + ConfigurableAxis CfgVtxBins{"CfgVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::QAObject}; + // Configurables + + // Pre-selection cuts + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable mincrossedrows{"mincrossedrows", 70, "min crossed rows"}; + + /// DCA Selections + // DCAr to PV + Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 0.1, "Track DCAr cut to PV Maximum"}; + // DCAz to PV + Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 5.0, "Track DCAz cut to PV Maximum"}; + Configurable cMinDCAzToPVcut{"cMinDCAzToPVcut", 0.0, "Track DCAz cut to PV Minimum"}; + + /// Partition for firstTrack + Partition parts1 = (aod::resodaughter::partType == uint8_t(aod::resodaughter::DaughterType::kTrack)) && requireTPCPIDKaonCutInFilter() && requireTOFPIDKaonCutInFilter() && (nabs(o2::aod::track::dcaZ) > static_cast(cMinDCAzToPVcut)) && (nabs(o2::aod::track::dcaZ) < static_cast(cMaxDCAzToPVcut)) && (nabs(o2::aod::track::dcaXY) < static_cast(cMaxDCArToPVcut)); // Basic DCA cuts + Partition parts2 = (aod::resodaughter::partType == uint8_t(aod::resodaughter::DaughterType::kTrack)) && requireTPCPIDKaonCutInFilter() && requireTOFPIDKaonCutInFilter() && (nabs(o2::aod::track::dcaZ) > static_cast(cMinDCAzToPVcut)) && (nabs(o2::aod::track::dcaZ) < static_cast(cMaxDCAzToPVcut)) && (nabs(o2::aod::track::dcaXY) < static_cast(cMaxDCArToPVcut)); // Basic DCA cuts + void init(o2::framework::InitContext&) + { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + long now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + ccdb->setCreatedNotAfter(now); + + AxisSpec vtxZAxis = {100, -20, 20}; + + std::vector centBinning = {0., 1., 5., 10., 20., 30., 40., 50., 70., 100.}; + AxisSpec centAxis = {centBinning, "V0M (%)"}; + std::vector ptBinning = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4., 5., 10., 20.}; + AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"}; + + // Mass QA (quick check) + histos.add("phiinvmass", "Invariant mass of Phi", kTH1F, {{500, 0.8, 1.3, "Invariant Mass (GeV/#it{c}^2)"}}); + histos.add("phiinvmassME", "Invariant mass of Phi mixed event", kTH1F, {{500, 0.8, 1.3, "Invariant Mass (GeV/#it{c}^2)"}}); + + // 3d histogram + histos.add("h3phiinvmass", "Invariant mass of Phi", kTH3F, {{100, 0.0f, 100.0f}, {100, 0.0f, 10.0f}, {500, 0.8, 1.3}}); + histos.add("h3phiinvmassME", "Invariant mass of Phi mixed event", kTH3F, {{100, 0.0f, 100.0f}, {100, 0.0f, 10.0f}, {500, 0.8, 1.3}}); + } + + double massKa = TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(); + + void process(aod::ResoCollision& collision, + aod::ResoDaughters const&, aod::Reso2TracksPIDExt const&) + { + // LOGF(info, "event id: %d", collision.bcId()); + auto group1 = parts1->sliceByCached(aod::resodaughter::resoCollisionId, collision.globalIndex()); + auto group2 = parts2->sliceByCached(aod::resodaughter::resoCollisionId, collision.globalIndex()); + + for (auto& [trk1, trk2] : combinations(CombinationsStrictlyUpperIndexPolicy(group1, group2))) { + // Un-like sign pair only + if (trk1.sign() * trk2.sign() > 0) + continue; + + auto arrMom = array{ + array{trk1.px(), trk1.py(), trk1.pz()}, + array{trk2.px(), trk2.py(), trk2.pz()}}; + auto resoPt = RecoDecay::sqrtSumOfSquares(trk1.px() + trk2.px(), trk1.py() + trk2.py()); + auto arrMass = array{massKa, massKa}; + auto mass = RecoDecay::m(arrMom, arrMass); + histos.fill(HIST("phiinvmass"), mass); + histos.fill(HIST("h3phiinvmass"), collision.multV0M(), resoPt, mass); + } + } + + // Processing Event Mixing + void processME(o2::aod::ResoCollisions& collision, + o2::aod::BCsWithTimestamps const&, aod::ResoDaughters const&, aod::Reso2TracksPIDExt const&) + { + BinningPolicy colBinning{{CfgVtxBins, CfgMultBins}, true}; + + for (auto& [collision1, collision2] : soa::selfCombinations(colBinning, 5, -1, collision, collision)) { + + auto magFieldTesla1 = 0.0; + static o2::parameters::GRPObject* grpo1 = ccdb->getForTimeStamp("GLO/GRP/GRP", collision1.timestamp()); + if (!grpo1) { + magFieldTesla1 = 0; + } else { + // LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", collision1.timestamp(), grpo1->getNominalL3Field()); + magFieldTesla1 = 0.1 * (grpo1->getNominalL3Field()); + } + + auto magFieldTesla2 = 0.0; + static o2::parameters::GRPObject* grpo2 = ccdb->getForTimeStamp("GLO/GRP/GRP", collision2.timestamp()); + if (!grpo2) { + magFieldTesla2 = 0; + } else { + // LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", collision2.timestamp(), grpo2->getNominalL3Field()); + magFieldTesla2 = 0.1 * (grpo1->getNominalL3Field()); + } + + auto group1 = parts1->sliceByCached(aod::resodaughter::resoCollisionId, collision1.globalIndex()); + auto group2 = parts2->sliceByCached(aod::resodaughter::resoCollisionId, collision2.globalIndex()); + + if (magFieldTesla1 != magFieldTesla2) { + continue; + } + for (auto& [trk1, trk2] : combinations(CombinationsStrictlyUpperIndexPolicy(group1, group2))) { + // Un-like sign pair only + if (trk1.sign() * trk2.sign() > 0) + continue; + + auto arrMom = array{ + array{trk1.px(), trk1.py(), trk1.pz()}, + array{trk2.px(), trk2.py(), trk2.pz()}}; + auto resoPt = RecoDecay::sqrtSumOfSquares(trk1.px() + trk2.px(), trk1.py() + trk2.py()); + auto arrMass = array{massKa, massKa}; + auto mass = RecoDecay::m(arrMom, arrMass); + histos.fill(HIST("phiinvmassME"), mass); + histos.fill(HIST("h3phiinvmassME"), collision1.multV0M(), resoPt, mass); + } + } + }; + PROCESS_SWITCH(phianalysis, processME, "Process EventMixing", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"lf-phianalysis"})}; +} \ No newline at end of file diff --git a/PWGLF/Utils/collisionCuts.h b/PWGLF/Utils/collisionCuts.h new file mode 100644 index 00000000000..5339e263eba --- /dev/null +++ b/PWGLF/Utils/collisionCuts.h @@ -0,0 +1,135 @@ +// 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 collisionCuts.h +/// \brief Traditional event selection cuts for O2 analysis +/// +/// Simply copied from FemtoDreamCollisionSelection.h +/// original author: Laura Serksnyte, TU München +/// +/// \author Bong-Hwi Lim + +#ifndef ANALYSIS_TASKS_PWGLF_COLLISIONCUTS_H_ +#define ANALYSIS_TASKS_PWGLF_COLLISIONCUTS_H_ + +#include "Common/CCDB/TriggerAliases.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/Logger.h" + +using namespace o2::framework; + +namespace o2::analysis +{ + +class CollisonCuts +{ + public: + virtual ~CollisonCuts() = default; + + /// Pass the selection criteria to the class + /// \param zvtxMax Maximal value of the z-vertex + /// \param checkTrigger whether or not to check for the trigger alias + /// \param trig Requested trigger alias + /// \param checkOffline whether or not to check for offline selection criteria + void setCuts(float zvtxMax, bool checkTrigger, int trig, bool checkOffline, bool checkRun3) + { + mCutsSet = true; + mZvtxMax = zvtxMax; + mCheckTrigger = checkTrigger; + mTrigger = static_cast(trig); + mCheckOffline = checkOffline; + mCheckIsRun3 = checkRun3; + } + + /// Initializes histograms for the task + /// \param registry Histogram registry to be passed + void init(HistogramRegistry* registry) + { + if (!mCutsSet) { + LOGF(error, "Event selection not set - quitting!"); + } + mHistogramRegistry = registry; + mHistogramRegistry->add("Event/posZ", "; vtx_{z} (cm); Entries", kTH1F, {{300, -12.5, 12.5}}); + mHistogramRegistry->add("Event/MultFV0M", "; vMultV0M; Entries", kTH1F, {{600, 0, 600}}); + mHistogramRegistry->add("Event/MultFT0M", "; vMultT0M; Entries", kTH1F, {{600, 0, 600}}); + } + + /// Print some debug information + void printCuts() + { + std::cout << "Debug information for CollisonCuts \n Max. z-vertex: " << mZvtxMax << "\n Check trigger: " << mCheckTrigger << "\n Trigger: " << mTrigger << "\n Check offline: " << mCheckOffline << "\n"; + } + + /// Check whether the collisions fulfills the specified selections + /// \tparam T type of the collision + /// \param col Collision + /// \return whether or not the collisions fulfills the specified selections + template + bool isSelected(T const& col) + { + if (std::abs(col.posZ()) > mZvtxMax) { + return false; + } + if (mCheckIsRun3) { + if (mCheckOffline && !col.sel8()) { + return false; + } + } else { + if (mCheckTrigger && !col.alias()[mTrigger]) { + return false; + } + if (mCheckOffline && !col.sel7()) { + return false; + } + } + return true; + } + + /// Some basic QA of the event + /// \tparam T type of the collision + /// \param col Collision + template + void fillQA(T const& col) + { + if (mHistogramRegistry) { + mHistogramRegistry->fill(HIST("Event/posZ"), col.posZ()); + mHistogramRegistry->fill(HIST("Event/MultFV0M"), col.multFV0M()); + mHistogramRegistry->fill(HIST("Event/MultFT0M"), col.multFT0M()); + } + } + + /// \todo to be implemented! + /// Compute the sphericity of an event + /// Important here is that the filter on tracks does not interfere here! + /// In Run 2 we used here global tracks within |eta| < 0.8 + /// \tparam T1 type of the collision + /// \tparam T2 type of the tracks + /// \param col Collision + /// \param tracks All tracks + /// \return value of the sphericity of the event + template + float computeSphericity(T1 const& col, T2 const& tracks) + { + return 2.f; + } + + private: + HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output + bool mCutsSet = false; ///< Protection against running without cuts + bool mCheckTrigger = false; ///< Check for trigger + bool mCheckOffline = false; ///< Check for offline criteria (might change) + bool mCheckIsRun3 = false; ///< Check if running on Pilot Beam + triggerAliases mTrigger = kINT7; ///< Trigger to check for + float mZvtxMax = 999.f; ///< Maximal deviation from nominal z-vertex (cm) +}; +} // namespace o2::analysis + +#endif /* ANALYSIS_TASKS_PWGLF_COLLISIONCUTS_H_ */ \ No newline at end of file