diff --git a/PWGCF/Flow/Tasks/CMakeLists.txt b/PWGCF/Flow/Tasks/CMakeLists.txt index 34a3dafcf30..c77ae659fdb 100644 --- a/PWGCF/Flow/Tasks/CMakeLists.txt +++ b/PWGCF/Flow/Tasks/CMakeLists.txt @@ -58,3 +58,8 @@ o2physics_add_dpl_workflow(flow-sp SOURCES flowSP.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::GFWCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(resonances-gfw-flow + SOURCES resonancesGfwFlow.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::GFWCore + COMPONENT_NAME Analysis) diff --git a/PWGCF/Flow/Tasks/resonancesGfwFlow.cxx b/PWGCF/Flow/Tasks/resonancesGfwFlow.cxx new file mode 100644 index 00000000000..1b8428f0d9d --- /dev/null +++ b/PWGCF/Flow/Tasks/resonancesGfwFlow.cxx @@ -0,0 +1,341 @@ +// 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 resonancesGfwFlow.cxx +/// \brief PID flow for resonances using the generic framework +/// \author Preet Bhanjan Pati + +#include +#include +#include +#include +#include +#include + +#include "Math/Vector4D.h" + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/StepTHn.h" + +#include "Common/DataModel/EventSelection.h" +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/Core/trackUtilities.h" +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Multiplicity.h" +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "CommonConstants/PhysicsConstants.h" + +#include "ReconstructionDataFormats/Track.h" +#include "ReconstructionDataFormats/PID.h" + +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace std; + +#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; + +struct ResonancesGfwFlow { + Service ccdb; + Configurable noLaterThan{"noLaterThan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; + Configurable ccdbUrl{"ccdbUrl", "http://ccdb-test.cern.ch:8080", "url of the ccdb repository"}; + + O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range") + O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.2f, "Minimal pT for ref tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 3.0f, "Maximal pT for ref tracks") + O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks") + O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5, "Chi2 per TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Use Nch for flow observables") + O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples") + O2_DEFINE_CONFIGURABLE(cfgTpcNsigmaCut, float, 3.0f, "TPC N-sigma cut for pions, kaons, protons") + O2_DEFINE_CONFIGURABLE(cfgTofNsigmaCut, float, 3.0f, "TOF N-sigma cut for pions, kaons, protons") + O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 0.5f, "Minimum pt to use TOF N-sigma") + + O2_DEFINE_CONFIGURABLE(cfgUsePVContributor, bool, true, "Use PV contributor for tracks") + O2_DEFINE_CONFIGURABLE(cfgFakePartCut, float, 0.1f, "Maximum difference in measured momentum and TPC inner ring momentum of particle") + O2_DEFINE_CONFIGURABLE(cfgTpcCluster, int, 70, "Number of TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgUsePointingAngle, bool, false, "Use Pointing angle for resonances") + O2_DEFINE_CONFIGURABLE(cfgPointingAnglePhi, float, 0.04f, "Minimum Pointing angle for Phi") + O2_DEFINE_CONFIGURABLE(cfgPointingAngleKo, float, 0.97f, "Minimum Pointing angle for K0") + O2_DEFINE_CONFIGURABLE(cfgPointingAngleLambda, float, 0.995f, "Minimum Pointing angle for Lambda") + O2_DEFINE_CONFIGURABLE(cfgUseVoRadius, bool, true, "Use V0 radius for particle identification") + O2_DEFINE_CONFIGURABLE(cfgVoRadiusMin, float, 0.5f, "Minimum V0 radius in cm") + O2_DEFINE_CONFIGURABLE(cfgVoRadiusMax, float, 200.0f, "Maximum V0 radius in cm") + O2_DEFINE_CONFIGURABLE(cfgUseProperLifetime, bool, true, "Use proper lifetime for particle identification") + O2_DEFINE_CONFIGURABLE(cfgProperLtK0, float, 20.0f, "Minimum lifetime for K0 in cm") + O2_DEFINE_CONFIGURABLE(cfgProperLtLambda, float, 30.0f, "Minimum lifetime for Lambda in cm") + O2_DEFINE_CONFIGURABLE(cfgUseDCA, bool, true, "Use dca for daughter tracks") + O2_DEFINE_CONFIGURABLE(cfgDCAtoPV, float, 0.06f, "Minimum DCA of daughter tracks to primary vertex in cm") + O2_DEFINE_CONFIGURABLE(cfgDCABetDaug, int, 1, "Maximum DCA between daughter tracks") + + O2_DEFINE_CONFIGURABLE(cfgCutDCAxy, float, 2.0f, "DCAxy range for tracks") + O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "DCAz range for tracks") + O2_DEFINE_CONFIGURABLE(additionalEvsel, bool, false, "Additional event selcection") + O2_DEFINE_CONFIGURABLE(useGlobalTrack, bool, true, "use Global track") + O2_DEFINE_CONFIGURABLE(cfgITScluster, int, 0, "Number of ITS cluster") + O2_DEFINE_CONFIGURABLE(cfgCutTOFBeta, float, 0.0, "cut TOF beta") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancy, int, 3000, "Occupancy cut") + O2_DEFINE_CONFIGURABLE(ispTdepPID, bool, true, "pT dependent PID") + O2_DEFINE_CONFIGURABLE(removefaketrack, bool, true, "Remove fake track from momentum difference") + O2_DEFINE_CONFIGURABLE(confRapidity, float, 0.5, "Rapidity cut") + + // Defining configurable axis + ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for histograms"}; + ConfigurableAxis axisPhi{"axisPhi", {60, 0.0, constants::math::TwoPI}, "phi axis for histograms"}; + ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"}; + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.20, 1.40, 1.60, 1.80, 2.00, 2.20, 2.40, 2.60, 2.80, 3.00, 3.50, 4.00, 5.00, 6.00, 8.00, 10.00}, "pt axis for histograms"}; + ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "centrality axis for histograms"}; + ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {80, -5, 5}, "nsigmaTPC axis"}; + ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"}; + ConfigurableAxis axisParticles{"axisParticles", {3, 0, 3}, "axis for different hadrons"}; + ConfigurableAxis axisPhiMass{"axisPhiMass", {50000, 0, 5}, "axis for invariant mass distibution for Phi"}; + ConfigurableAxis axisKoMass{"axisKoMass", {50000, 0, 5}, "axis for invariant mass distibution for K0"}; + ConfigurableAxis axisLambdaMass{"axisLambdaMass", {50000, 0, 5}, "axis for invariant mass distibution for Lambda"}; + ConfigurableAxis axisTPCsignal{"axisTPCsignal", {10000, 0, 1000}, "axis for TPC signal"}; + ConfigurableAxis axisTOFsignal{"axisTOFsignal", {10000, 0, 1000}, "axis for TOF signal"}; + + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter trackFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz) && (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtPOIMin) && (aod::track::pt < cfgCutPtPOIMax); + + using AodCollisions = soa::Filtered>; + using AodTracksWithoutBayes = soa::Filtered>; + + SliceCache cache; + Partition posTracks = aod::track::signed1Pt > 0.0f; + Partition negTracks = aod::track::signed1Pt < 0.0f; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + TAxis* fPtAxis; + TRandom3* fRndm = new TRandom3(0); + + // Event selection cuts - Alex + TF1* fMultPVCutLow = nullptr; + TF1* fMultPVCutHigh = nullptr; + TF1* fMultCutLow = nullptr; + TF1* fMultCutHigh = nullptr; + TF1* fMultMultPVCut = nullptr; + + void init(InitContext const&) + { + ccdb->setURL(ccdbUrl.value); + ccdb->setCaching(true); + ccdb->setCreatedNotAfter(noLaterThan.value); + + histos.add("hPhi", "", {HistType::kTH1D, {axisPhi}}); + histos.add("hEta", "", {HistType::kTH1D, {axisEta}}); + histos.add("hVtxZ", "", {HistType::kTH1D, {axisVertex}}); + histos.add("hMult", "", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); + histos.add("hCent", "", {HistType::kTH1D, {{90, 0, 90}}}); + + histos.add("KaplusTPC", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); + histos.add("KaminusTPC", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); + histos.add("KaplusTOF", "", {HistType::kTH2D, {{axisPt, axisTOFsignal}}}); + histos.add("KaminusTOF", "", {HistType::kTH2D, {{axisPt, axisTOFsignal}}}); + histos.add("hPhiMass_sparse", "", {HistType::kTHnSparseD, {{axisPhiMass, axisPt, axisMultiplicity}}}); + + if (additionalEvsel) { + fMultPVCutLow = new TF1("fMultPVCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 2.5*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100); + fMultPVCutLow->SetParameters(2834.66, -87.0127, 0.915126, -0.00330136, 332.513, -12.3476, 0.251663, -0.00272819, 1.12242e-05); + fMultPVCutHigh = new TF1("fMultPVCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 2.5*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100); + fMultPVCutHigh->SetParameters(2834.66, -87.0127, 0.915126, -0.00330136, 332.513, -12.3476, 0.251663, -0.00272819, 1.12242e-05); + fMultCutLow = new TF1("fMultCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 2.5*([4]+[5]*x)", 0, 100); + fMultCutLow->SetParameters(1893.94, -53.86, 0.502913, -0.0015122, 109.625, -1.19253); + fMultCutHigh = new TF1("fMultCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 3.*([4]+[5]*x)", 0, 100); + fMultCutHigh->SetParameters(1893.94, -53.86, 0.502913, -0.0015122, 109.625, -1.19253); + fMultMultPVCut = new TF1("fMultMultPVCut", "[0]+[1]*x+[2]*x*x", 0, 5000); + fMultMultPVCut->SetParameters(-0.1, 0.785, -4.7e-05); + } + } + + // deep angle cut on pair to remove photon conversion + template + bool selectionPair(const T1& candidate1, const T2& candidate2) + { + double pt1, pt2, pz1, pz2, p1, p2, angle; + pt1 = candidate1.pt(); + pt2 = candidate2.pt(); + pz1 = candidate1.pz(); + pz2 = candidate2.pz(); + p1 = candidate1.p(); + p2 = candidate2.p(); + angle = std::acos((pt1 * pt2 + pz1 * pz2) / (p1 * p2)); + if (cfgUsePointingAngle && angle < cfgPointingAnglePhi) { + return false; + } + return true; + } + + template + bool eventSelected(TCollision collision, const float& centrality) + { + auto multNTracksPV = collision.multNTracksPV(); + if (multNTracksPV < fMultPVCutLow->Eval(centrality)) + return 0; + if (multNTracksPV > fMultPVCutHigh->Eval(centrality)) + return 0; + + return 1; + } + + template + bool selectionTrack(const T& candidate) + { + if (useGlobalTrack && !(candidate.isGlobalTrack() && candidate.isPVContributor() && candidate.itsNCls() > cfgITScluster && candidate.tpcNClsFound() > cfgTpcCluster)) { + return false; + } + if (!useGlobalTrack && !(candidate.isPVContributor() && candidate.itsNCls() > cfgITScluster)) { + return false; + } + return true; + } + + template + bool selectionPIDpTdependent(const T& candidate) + { + if (candidate.pt() < 0.5 && std::abs(candidate.tpcNSigmaKa()) < cfgTpcNsigmaCut) { + return true; + } + if (candidate.pt() >= 0.5 && candidate.hasTOF() && candidate.beta() > cfgCutTOFBeta && std::abs(candidate.tpcNSigmaKa()) < cfgTpcNsigmaCut && std::abs(candidate.tofNSigmaKa()) < cfgTofNsigmaCut) { + return true; + } + if (!useGlobalTrack && !candidate.hasTPC()) { + return true; + } + return false; + } + template + bool selectionPID(const T& candidate) + { + if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaKa()) < cfgTpcNsigmaCut) { + return true; + } + if (candidate.hasTOF() && candidate.beta() > cfgCutTOFBeta && std::abs(candidate.tpcNSigmaKa()) < cfgTpcNsigmaCut && std::abs(candidate.tofNSigmaKa()) < cfgTofNsigmaCut) { + return true; + } + return false; + } + + template + bool isFakeKaon(T const& track) + { + const auto pglobal = track.p(); + const auto ptpc = track.tpcInnerParam(); + if (std::abs(pglobal - ptpc) > cfgFakePartCut) { + return true; + } + return false; + } + + using BinningTypeVertexContributor = ColumnBinningPolicy; + ROOT::Math::PxPyPzMVector phiMom, kaonPlus, kaonminus; + double massKaPlus = o2::constants::physics::MassKPlus; + + void process(AodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, AodTracksWithoutBayes const& tracks) + { + if (!collision.sel8() || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(aod::evsel::kNoITSROFrameBorder) || !collision.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) || !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + return; + } + + const auto cent = collision.centFT0C(); + int nTot = tracks.size(); + float vtxz = collision.posZ(); + int occupancy = collision.trackOccupancyInTimeRange(); + + if (occupancy > cfgCutOccupancy) { + return; + } + + if (additionalEvsel && !eventSelected(collision, cent)) { + return; + } + + histos.fill(HIST("hVtxZ"), vtxz); + histos.fill(HIST("hMult"), nTot); + histos.fill(HIST("hCent"), cent); + + auto posSlicedTracks = posTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto negSlicedTracks = negTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + + for (auto const& track1 : posSlicedTracks) { + // track selection + if (!selectionTrack(track1)) { + continue; + } + // PID check + if (ispTdepPID && !selectionPIDpTdependent(track1)) { + continue; + } + if (!ispTdepPID && !selectionPID(track1)) { + continue; + } + histos.fill(HIST("KaplusTPC"), track1.pt(), track1.tpcNSigmaKa()); + histos.fill(HIST("KaplusTOF"), track1.pt(), track1.tofNSigmaKa()); + auto track1ID = track1.globalIndex(); + + for (auto const& track2 : negSlicedTracks) { + // track selection + if (!selectionTrack(track2)) { + continue; + } + // PID check + if (ispTdepPID && !selectionPIDpTdependent(track2)) { + continue; + } + if (!ispTdepPID && !selectionPID(track2)) { + continue; + } + auto track2ID = track2.globalIndex(); + if (track2ID == track1ID) { + continue; + } + if (!selectionPair(track1, track2)) { + continue; + } + if (removefaketrack && isFakeKaon(track1)) { + continue; + } + if (removefaketrack && isFakeKaon(track2)) { + continue; + } + histos.fill(HIST("KaminusTPC"), track2.pt(), track2.tpcNSigmaKa()); + histos.fill(HIST("KaminusTOF"), track2.pt(), track2.tofNSigmaKa()); + kaonPlus = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKaPlus); + kaonminus = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massKaPlus); + phiMom = kaonPlus + kaonminus; + if (std::abs(phiMom.Rapidity()) < confRapidity) { + histos.fill(HIST("hPhiMass_sparse"), phiMom.M(), phiMom.Pt(), cent); + histos.fill(HIST("hPhi"), phiMom.Phi()); + histos.fill(HIST("hEta"), phiMom.Eta()); + } + } // end of track 2 + } // end of track 1 + } // end of process +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}