[PWGLF,Trigger] Trigger task to select events with 2 phi meson (#10607)

skundu692 · web-flow · commit e48b8c132728 · 2025-03-24T09:39:39.000+01:00
diff --git a/EventFiltering/CMakeLists.txt b/EventFiltering/CMakeLists.txt
@@ -112,11 +112,16 @@ o2physics_add_dpl_workflow(lf-f1proton-filter
     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase
     COMPONENT_NAME Analysis)
 
+ o2physics_add_dpl_workflow(lf-doublephi-filter
+    SOURCES PWGLF/filterdoublephi.cxx
+    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase
+    COMPONENT_NAME Analysis)
+
 o2physics_add_library(EventFilteringUtils
     SOURCES Zorro.cxx ZorroSummary.cxx
     INSTALL_HEADERS ZorroHelper.h ZorroSummary.h
     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore Arrow::arrow_shared)
 
 o2physics_target_root_dictionary(EventFilteringUtils
    HEADERS ZorroHelper.h ZorroSummary.h
-   LINKDEF EventFilteringUtilsLinkDef.h)
+   LINKDEF EventFilteringUtilsLinkDef.h)
diff --git a/EventFiltering/PWGLF/filterdoublephi.cxx b/EventFiltering/PWGLF/filterdoublephi.cxx
@@ -0,0 +1,279 @@
+// 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 filterdoublephi.cxx
+/// \brief Selection of events with triplets and pairs for femtoscopic studies
+///
+/// \author Sourav Kundu, sourav.kundu@cern.ch
+
+#include <Framework/Configurable.h>
+#include <Math/GenVector/Boost.h>
+#include <Math/Vector4D.h>
+#include <TMath.h>
+#include <fairlogger/Logger.h>
+#include <TDatabasePDG.h> // FIXME
+#include <TPDGCode.h>     // FIXME
+
+#include <iostream>
+#include <iterator>
+#include <string>
+#include <vector>
+
+#include "../filterTables.h"
+#include "PWGLF/DataModel/ReducedDoublePhiTables.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/runDataProcessing.h"
+#include "CommonConstants/MathConstants.h"
+#include "Common/Core/TrackSelection.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+#include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/Multiplicity.h"
+#include "Common/DataModel/PIDResponse.h"
+#include "DataFormatsTPC/BetheBlochAleph.h"
+#include "CCDB/BasicCCDBManager.h"
+#include "CCDB/CcdbApi.h"
+#include "Common/DataModel/PIDResponseITS.h"
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+struct filterdoublephi {
+
+  // Produce derived tables
+  Produces<aod::DoublePhiFilters> tags;
+
+  // events
+  Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"};
+  // Configurable<float> cfgCutCentralityMax{"cfgCutCentralityMax", 0.0f, "Accepted maximum Centrality"};
+  // Configurable<float> cfgCutCentralityMin{"cfgCutCentralityMin", 100.0f, "Accepted minimum Centrality"};
+  // track
+  Configurable<bool> isPtdepPID1{"isPtdepPID1", false, "use pt dep PID kplus"};
+  Configurable<bool> isPtdepPID2{"isPtdepPID2", false, "use pt dep PID kminus"};
+  Configurable<bool> useGlobalTrack{"useGlobalTrack", true, "use Global track"};
+  Configurable<float> cfgCutTOFBeta{"cfgCutTOFBeta", 0.0, "cut TOF beta"};
+  Configurable<float> cfgCutCharge{"cfgCutCharge", 0.0, "cut on Charge"};
+  Configurable<float> cfgCutPT{"cfgCutPT", 0.2, "PT cut on daughter track"};
+  Configurable<float> cfgCutEta{"cfgCutEta", 0.8, "Eta cut on daughter track"};
+  Configurable<float> cfgCutDCAxy{"cfgCutDCAxy", 2.0f, "DCAxy range for tracks"};
+  Configurable<float> cfgCutDCAz{"cfgCutDCAz", 2.0f, "DCAz range for tracks"};
+  Configurable<float> nsigmaCutTPC{"nsigmacutTPC", -2.0, "Value of the TPC Nsigma cut"};
+  Configurable<float> nsigmaCutTPCPreSel{"nsigmacutTPCPreSel", 3.0, "Value of the TPC Nsigma cut Pre selection"};
+  Configurable<float> nsigmaCutTOF{"nsigmaCutTOF", 3.0, "Value of the TOF Nsigma cut"};
+  Configurable<int> cfgITScluster{"cfgITScluster", 0, "Number of ITS cluster"};
+  Configurable<int> cfgTPCcluster{"cfgTPCcluster", 70, "Number of TPC cluster"};
+  Configurable<bool> isDeepAngle{"isDeepAngle", true, "Deep Angle cut"};
+  Configurable<double> cfgDeepAngle{"cfgDeepAngle", 0.04, "Deep Angle cut value"};
+  ConfigurableAxis configThnAxisInvMass{"configThnAxisInvMass", {120, 0.98, 1.1}, "#it{M} (GeV/#it{c}^{2})"};
+  ConfigurableAxis configThnAxisPt{"configThnAxisPt", {100, 0.0, 10.}, "#it{p}_{T} (GeV/#it{c})"};
+  Configurable<float> minPhiMass{"minPhiMass", 1.01, "Minimum phi mass"};
+  Configurable<float> maxPhiMass{"maxPhiMass", 1.03, "Maximum phi mass"};
+  Configurable<float> MinPhiPairPt{"MinPhiPairPt", 2.0, "Minimum phi pair Pt"};
+  Configurable<float> MinPhiPairMass{"MinPhiPairMass", 2.5, "Minimum phi pair mass"};
+  Configurable<float> MaxPhiPairMass{"MaxPhiPairMass", 3.0, "Max phi pair mass"};
+
+  // Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
+  // Filter centralityFilter = (nabs(aod::cent::centFT0C) < cfgCutCentralityMax && nabs(aod::cent::centFT0C) > cfgCutCentralityMin);
+  Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPT);
+  Filter DCAcutFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz);
+  Filter PIDcutFilter = nabs(aod::pidtpc::tpcNSigmaKa) < nsigmaCutTPCPreSel;
+
+  // using EventCandidates = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults>>;
+  using EventCandidates = soa::Join<aod::Collisions, aod::EvSels>;
+  using TrackCandidates = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTOFbeta, aod::pidTPCFullKa, aod::pidTOFFullKa>>;
+
+  SliceCache cache;
+  Partition<TrackCandidates> posTracks = aod::track::signed1Pt > cfgCutCharge;
+  Partition<TrackCandidates> negTracks = aod::track::signed1Pt < cfgCutCharge;
+
+  // Histogram
+  OutputObj<TH1D> hProcessedEvents{TH1D("hProcessedEvents", ";; Number of events", 3, 0.0f, 3.0f)};
+  HistogramRegistry qaRegistry{"QAHistos", {
+                                             {"hInvMassPhi", "hInvMassPhi", {HistType::kTH2F, {{40, 1.0f, 1.04f}, {100, 0.0f, 10.0f}}}},
+                                             {"hInvMassDoublePhi", "hInvMassDoublePhi", {HistType::kTH2F, {{1000, 2.0f, 3.0f}, {100, 0.0f, 10.0f}}}},
+                                             {"hNsigmaPtkaonTPC", "hNsigmaPtkaonTPC", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
+                                             {"hNsigmaPtkaonTOF", "hNsigmaPtkaonTOF", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
+                                           },
+                               OutputObjHandlingPolicy::AnalysisObject};
+
+  double massKa = o2::constants::physics::MassKPlus;
+
+  void init(o2::framework::InitContext&)
+  {
+    hProcessedEvents->GetXaxis()->SetBinLabel(1, "All events");
+    hProcessedEvents->GetXaxis()->SetBinLabel(2, "Events with double Phi without sel.");
+    hProcessedEvents->GetXaxis()->SetBinLabel(3, aod::filtering::TriggerEventDoublePhi::columnLabel());
+  }
+
+  template <typename T>
+  bool selectionTrack(const T& candidate)
+  {
+    if (useGlobalTrack && !(candidate.isGlobalTrack() && candidate.isPVContributor() && candidate.itsNCls() > cfgITScluster && candidate.tpcNClsCrossedRows() > cfgTPCcluster)) {
+      return false;
+    }
+    return true;
+  }
+  template <typename T>
+  bool selectionPID(const T& candidate)
+  {
+    if (candidate.pt() < 0.5 && candidate.tpcNSigmaKa() > nsigmaCutTPC && candidate.tpcNSigmaKa() < 3.0) {
+      return true;
+    }
+    if (candidate.pt() >= 0.5) {
+      if (!candidate.hasTOF() && candidate.tpcNSigmaKa() > nsigmaCutTPC && candidate.tpcNSigmaKa() < 3.0) {
+        return true;
+      }
+      if (candidate.hasTOF() && candidate.beta() > cfgCutTOFBeta && candidate.tpcNSigmaKa() > nsigmaCutTPC && candidate.tpcNSigmaKa() < 3.0 && TMath::Abs(candidate.tofNSigmaKa()) < nsigmaCutTOF) {
+        return true;
+      }
+    }
+    return false;
+  }
+  template <typename T>
+  bool selectionPID2(const T& candidate)
+  {
+    if (candidate.pt() < 0.5 && candidate.tpcNSigmaKa() > nsigmaCutTPC && candidate.tpcNSigmaKa() < 3.0) {
+      return true;
+    }
+    if (candidate.pt() >= 0.5) {
+      if (candidate.hasTOF() && candidate.beta() > cfgCutTOFBeta && candidate.tpcNSigmaKa() > nsigmaCutTPC && candidate.tpcNSigmaKa() < 3.0 && TMath::Abs(candidate.tofNSigmaKa()) < nsigmaCutTOF) {
+        return true;
+      }
+    }
+    return false;
+  }
+  // deep angle cut on pair to remove photon conversion
+  template <typename T1, typename T2>
+  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 = TMath::ACos((pt1 * pt2 + pz1 * pz2) / (p1 * p2));
+    if (isDeepAngle && angle < cfgDeepAngle) {
+      return false;
+    }
+    return true;
+  }
+
+  ROOT::Math::PxPyPzMVector KaonPlus, KaonMinus, PhiMesonMother, PhiVectorDummy, PhiVectorDummy2, PhiPair;
+  void processPhiReducedTable(EventCandidates::iterator const& collision, TrackCandidates const&, aod::BCsWithTimestamps const&)
+  {
+    bool keepEventDoublePhi = false;
+    int numberPhi = 0;
+    o2::aod::ITSResponse itsResponse;
+    std::vector<int64_t> Phid1Index = {};
+    std::vector<int64_t> Phid2Index = {};
+    std::vector<ROOT::Math::PtEtaPhiMVector> phiresonance, phiresonanced1, phiresonanced2;
+    int Npostrack = 0;
+    int Nnegtrack = 0;
+    hProcessedEvents->Fill(0.5);
+    if (collision.sel8()) {
+      auto posThisColl = posTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
+      auto negThisColl = negTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
+      for (auto track1 : posThisColl) {
+        // track selection
+        if (!selectionTrack(track1)) {
+          continue;
+        }
+        // PID check
+        if (isPtdepPID1 && !selectionPID2(track1)) {
+          continue;
+        }
+        if (!isPtdepPID1 && !selectionPID(track1)) {
+          continue;
+        }
+        if (track1.pt() > 0.4 && track1.pt() < 1.0 && !(itsResponse.nSigmaITS<o2::track::PID::Kaon>(track1) > -2.0 && itsResponse.nSigmaITS<o2::track::PID::Kaon>(track1) < 3.0)) {
+          continue;
+        }
+        Npostrack = Npostrack + 1;
+        qaRegistry.fill(HIST("hNsigmaPtkaonTPC"), track1.tpcNSigmaKa(), track1.pt());
+        if (track1.hasTOF()) {
+          qaRegistry.fill(HIST("hNsigmaPtkaonTOF"), track1.tofNSigmaKa(), track1.pt());
+        }
+        auto track1ID = track1.globalIndex();
+        for (auto track2 : negThisColl) {
+          // track selection
+          if (!selectionTrack(track2)) {
+            continue;
+          }
+          // PID check
+          if (isPtdepPID2 && !selectionPID2(track2)) {
+            continue;
+          }
+          if (!isPtdepPID2 && !selectionPID(track2)) {
+            continue;
+          }
+          if (track2.pt() > 0.4 && track2.pt() < 1.0 && !(itsResponse.nSigmaITS<o2::track::PID::Kaon>(track2) > -2.0 && itsResponse.nSigmaITS<o2::track::PID::Kaon>(track2) < 3.0)) {
+            continue;
+          }
+          if (Npostrack == 1) {
+            Nnegtrack = Nnegtrack + 1;
+          }
+          auto track2ID = track2.globalIndex();
+          if (track2ID == track1ID) {
+            continue;
+          }
+          if (!selectionPair(track1, track2)) {
+            continue;
+          }
+          KaonPlus = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKa);
+          KaonMinus = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massKa);
+          PhiMesonMother = KaonPlus + KaonMinus;
+          if (PhiMesonMother.M() > minPhiMass && PhiMesonMother.M() < maxPhiMass) {
+            numberPhi = numberPhi + 1;
+            ROOT::Math::PtEtaPhiMVector temp1(track1.pt(), track1.eta(), track1.phi(), massKa);
+            ROOT::Math::PtEtaPhiMVector temp2(track2.pt(), track2.eta(), track2.phi(), massKa);
+            ROOT::Math::PtEtaPhiMVector temp3(PhiMesonMother.pt(), PhiMesonMother.eta(), PhiMesonMother.phi(), PhiMesonMother.M());
+            phiresonanced1.push_back(temp1);
+            phiresonanced2.push_back(temp2);
+            phiresonance.push_back(temp3);
+            Phid1Index.push_back(track1.globalIndex());
+            Phid2Index.push_back(track2.globalIndex());
+            qaRegistry.fill(HIST("hInvMassPhi"), PhiMesonMother.M(), PhiMesonMother.Pt());
+          }
+        }
+      }
+    } // select collision
+    if (numberPhi > 1 && Npostrack > 1 && Nnegtrack > 1 && (phiresonance.size() == phiresonanced1.size()) && (phiresonance.size() == phiresonanced2.size())) {
+      hProcessedEvents->Fill(1.5);
+      for (auto if1 = phiresonance.begin(); if1 != phiresonance.end(); ++if1) {
+        auto i5 = std::distance(phiresonance.begin(), if1);
+        PhiVectorDummy = phiresonance.at(i5);
+        for (auto if2 = if1 + 1; if2 != phiresonance.end(); ++if2) {
+          auto i6 = std::distance(phiresonance.begin(), if2);
+          PhiVectorDummy2 = phiresonance.at(i6);
+          PhiPair = PhiVectorDummy + PhiVectorDummy2;
+          if ((Phid1Index.at(i5) != Phid1Index.at(i6)) && (Phid2Index.at(i5) != Phid2Index.at(i6)) && PhiPair.M() > MinPhiPairMass && PhiPair.M() < MaxPhiPairMass && PhiPair.Pt() > MinPhiPairPt) {
+            qaRegistry.fill(HIST("hInvMassDoublePhi"), PhiPair.M(), PhiPair.Pt());
+            keepEventDoublePhi = true;
+          }
+        }
+      }
+    }
+    if (keepEventDoublePhi) {
+      hProcessedEvents->Fill(2.5);
+    }
+    tags(keepEventDoublePhi);
+  } // process
+  PROCESS_SWITCH(filterdoublephi, processPhiReducedTable, "Process table creation for double phi", true);
+};
+WorkflowSpec defineDataProcessing(ConfigContext const& cfg)
+{
+  return WorkflowSpec{adaptAnalysisTask<filterdoublephi>(cfg, TaskName{"lf-doublephi-filter"})};
+}
diff --git a/EventFiltering/filterTables.h b/EventFiltering/filterTables.h
@@ -156,6 +156,9 @@ DECLARE_SOA_COLUMN(OmegaHighMult, hasOmegaHighMult, bool);       //! at least 1
 // F1-proton
 DECLARE_SOA_COLUMN(TriggerEventF1Proton, triggereventf1proton, bool); //! F1 - proton femto trigger event
 
+// Double Phi
+DECLARE_SOA_COLUMN(TriggerEventDoublePhi, triggereventdoublephi, bool); //! Double Phi trigger event
+
 // multiplicity
 DECLARE_SOA_COLUMN(HighTrackMult, hasHighTrackMult, bool);     //! high trk muliplicity
 DECLARE_SOA_COLUMN(HighMultFv0, hasHighMultFv0, bool);         //! high FV0 muliplicity
@@ -294,6 +297,11 @@ DECLARE_SOA_TABLE(F1ProtonFilters, "AOD", "F1ProtonFilters", //!
                   filtering::TriggerEventF1Proton);
 using F1ProtonFilter = F1ProtonFilters::iterator;
 
+// Double Phi
+DECLARE_SOA_TABLE(DoublePhiFilters, "AOD", "DoublePhiFilters", //!
+                  filtering::TriggerEventDoublePhi);
+using DoublePhiFilter = DoublePhiFilters::iterator;
+
 // multiplicity
 DECLARE_SOA_TABLE(MultFilters, "AOD", "MultFilters", //!
                   filtering::HighTrackMult, filtering::HighMultFv0, filtering::HighFt0Mult, filtering::HighFt0Flat, filtering::HighFt0cFv0Mult, filtering::HighFt0cFv0Flat, filtering::LeadingPtTrack);
@@ -324,11 +332,11 @@ DECLARE_SOA_TABLE(BCRanges, "AOD", "BCRanges", //!
 using BCRange = BCRanges::iterator;
 
 /// List of the available filters, the description of their tables and the name of the tasks
-constexpr int NumberOfFilters{13};
-constexpr std::array<char[32], NumberOfFilters> AvailableFilters{"NucleiFilters", "DiffractionFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "StrangenessFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "HeavyNeutralMesonFilters"};
-constexpr std::array<char[16], NumberOfFilters> FilterDescriptions{"NucleiFilters", "DiffFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "LFStrgFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "HNMesonFilters"};
-constexpr std::array<char[128], NumberOfFilters> FilteringTaskNames{"o2-analysis-nuclei-filter", "o2-analysis-diffraction-filter", "o2-analysis-dq-filter-pp-with-association", "o2-analysis-hf-filter", "o2-analysis-cf-filter", "o2-analysis-je-filter", "o2-analysis-je-hf-filter", "o2-analysis-fje-filter", "o2-analysis-lf-strangeness-filter", "o2-analysis-mult-filter", "o2-analysis-em-photon-filter", "o2-analysis-lf-f1proton-filter", "o2-analysis-heavy-meson-filter"};
-constexpr o2::framework::pack<NucleiFilters, DiffractionFilters, DqFilters, HfFilters, CFFilters, JetFilters, JetHFFilters, FullJetFilters, StrangenessFilters, MultFilters, PhotonFilters, F1ProtonFilters, HeavyNeutralMesonFilters> FiltersPack;
+constexpr int NumberOfFilters{14};
+constexpr std::array<char[32], NumberOfFilters> AvailableFilters{"NucleiFilters", "DiffractionFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "StrangenessFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "DoublePhiFilters", "HeavyNeutralMesonFilters"};
+constexpr std::array<char[16], NumberOfFilters> FilterDescriptions{"NucleiFilters", "DiffFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "LFStrgFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "2PhiFilters", "HNMesonFilters"};
+constexpr std::array<char[128], NumberOfFilters> FilteringTaskNames{"o2-analysis-nuclei-filter", "o2-analysis-diffraction-filter", "o2-analysis-dq-filter-pp-with-association", "o2-analysis-hf-filter", "o2-analysis-cf-filter", "o2-analysis-je-filter", "o2-analysis-je-hf-filter", "o2-analysis-fje-filter", "o2-analysis-lf-strangeness-filter", "o2-analysis-mult-filter", "o2-analysis-em-photon-filter", "o2-analysis-lf-f1proton-filter", "o2-analysis-lf-doublephi-filter", "o2-analysis-heavy-meson-filter"};
+constexpr o2::framework::pack<NucleiFilters, DiffractionFilters, DqFilters, HfFilters, CFFilters, JetFilters, JetHFFilters, FullJetFilters, StrangenessFilters, MultFilters, PhotonFilters, F1ProtonFilters, DoublePhiFilters, HeavyNeutralMesonFilters> FiltersPack;
 static_assert(o2::framework::pack_size(FiltersPack) == NumberOfFilters);
 
 template <typename T, typename C>
