Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions EventFiltering/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -111,3 +111,8 @@ o2physics_add_dpl_workflow(lf-f1proton-filter
SOURCES PWGLF/filterdoublephi.cxx
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase
COMPONENT_NAME Analysis)

o2physics_add_dpl_workflow(h2fromlbfilter
SOURCES PWGHF/H2fromLbFilter.cxx
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
COMPONENT_NAME Analysis)
247 changes: 247 additions & 0 deletions EventFiltering/PWGHF/H2fromLbFilter.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,247 @@
// 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.
//
/// \brief A filter task for non prompt deuterons
/// \author Marta Razza marta.razza@cern.ch
/// \author Francesca Ercolessi francesca.ercolessi@cern.ch
/// \since Dec 17, 2025

// #include "../filterTables.h"

#include "PWGLF/DataModel/LFParticleIdentification.h"

#include "Common/Core/RecoDecay.h"
#include "Common/Core/TrackSelection.h"
#include "Common/Core/trackUtilities.h"
#include "Common/DataModel/Centrality.h"
#include "Common/DataModel/CollisionAssociationTables.h"
#include "Common/DataModel/EventSelection.h"
#include "Common/DataModel/Multiplicity.h"
#include "Common/DataModel/PIDResponseTOF.h"
#include "Common/DataModel/TrackSelectionTables.h"

#include "CCDB/BasicCCDBManager.h"
#include "CommonConstants/PhysicsConstants.h"
#include "DCAFitter/DCAFitterN.h"
#include "DataFormatsParameters/GRPMagField.h"
#include "DataFormatsParameters/GRPObject.h"
#include "DetectorsBase/Propagator.h"
#include "Framework/ASoAHelpers.h"
#include "Framework/AnalysisDataModel.h"
#include "Framework/AnalysisTask.h"
#include "Framework/runDataProcessing.h"
#include "ReconstructionDataFormats/Track.h"
#include "ReconstructionDataFormats/TrackParametrization.h"

#include "TVector3.h"

#include <cmath>

struct H2fromLbFilter {

// Recall the output table
// Produces<aod::H2fromLbFilters> table;

o2::framework::Configurable<int> nBinsPt{"nBinsPt", 100, "N bins in pT histo"};

o2::base::Propagator::MatCorrType noMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;

// Define a histograms and registries
o2::framework::HistogramRegistry QAHistos{"QAHistos", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, true};
o2::framework::OutputObj<TH1D> hProcessedEvents{TH1D("hProcessedEvents", "Event filtered;; Number of events", 6, 0., 6.)};

o2::framework::ConfigurableAxis pTAxis{"pTAxis", {200, -10.f, 10.f}, "p_{T} GeV/c"};
o2::framework::ConfigurableAxis nSigmaAxis{"nSigmaAxis", {200, -10.f, 10.f}, "p_{T} GeV/c"};
o2::framework::ConfigurableAxis DCAxyAxis{"DCAxyAxis", {1000, -0.2f, 0.2f}, "DCA xy (cm)"};
o2::framework::ConfigurableAxis DCAzAxis{"DCAzAxis", {1000, -0.2f, 0.2f}, "DCA z (cm)"};

o2::framework::Configurable<float> cutzvertex{"cutzvertex", 20.0f, "Accepted z-vertex range"}; // 20 cm
o2::framework::Configurable<bool> isTVX{"isTVX", true, "isTVX event selection"};
o2::framework::Configurable<bool> isNoTimeFrameBorder{"isNoTimeFrameBorder", true, "isNoTimeFrameBorder event selection"};
o2::framework::Configurable<bool> isNoITSROFrameBorder{"isNoITSROFrameBorder", true, "isNoITSROFrameBorder event selection"};
o2::framework::Configurable<float> cfgTPCNsigma{"cfgTPCNsigma", 4.0f, "TPC n sigma for deuteron PID"};
o2::framework::Configurable<float> cfgTOFNsigma{"cfgTOFNsigma", 4.0f, "TOF n sigma for deuteron PID"};
o2::framework::Configurable<float> cfgEta{"cfgEta", 0.8f, "Track eta selection"};
o2::framework::Configurable<int> cfgTPCNclsFound{"cfgTPCNclsFound", 100, "Minimum TPC clusters found"};
o2::framework::Configurable<float> cfgTPCChi2Ncl{"cfgTPCChi2Ncl", 4.0f, "Maximum TPC chi2 per N clusters"};
o2::framework::Configurable<float> cfgITSChi2Ncl{"cfgITSChi2Ncl", 36.0f, "Maximum ITS chi2 per N clusters"};
o2::framework::Configurable<int> cfgITScls{"cfgITScls", 2, "Minimum ITS clusters"};
o2::framework::Configurable<float> cfgMaxPt{"cfgMaxPt", 5.0f, "Maximum pT cut"};
o2::framework::Configurable<float> cfgMinPt{"cfgMinPt", 1.0f, "Maximum pT cut"};
o2::framework::Configurable<float> cfgDCAcut{"cfgDCAcut", 0.003f, "DCA cut for non prompt deuteron"};
using myCompleteTracks = o2::soa::Join<o2::aod::Tracks, o2::aod::TracksExtra, o2::aod::TracksDCA, o2::aod::McTrackLabels>;
o2::framework::Service<o2::ccdb::BasicCCDBManager> ccdb;

void init(o2::framework::InitContext&)
{
ccdb->setURL("http://alice-ccdb.cern.ch");
ccdb->setCaching(true);
ccdb->setLocalObjectValidityChecking();
ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());

o2::framework::AxisSpec ptAxis = {100, 0.0f, 10.0f, "#it{p}_{T} (GeV/#it{c})"};

// general QA histograms
QAHistos.add("hVtxZ", "Z-Vertex distribution after selection;Z (cm)", o2::framework::HistType::kTH1F, {{100, -50, 50}});
QAHistos.add("hDCAxyVsPt", "DCAxy #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAxyAxis}});
QAHistos.add("hDCAzVsPt", "DCAz #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAzAxis}});
QAHistos.add("hnSigmaTPCVsPt", "n#sigma TPC vs p_{T} for #bar{d} hypothesis; p_{T} (GeV/c); n#sigma TPC", {o2::framework::HistType::kTH2D, {pTAxis, nSigmaAxis}});
QAHistos.add("hnSigmaTOFVsPt", "n#sigma TOF vs p_{T} for #bar{d} hypothesis; p_{T} (GeV/c); n#sigma TOF", {o2::framework::HistType::kTH2D, {pTAxis, nSigmaAxis}});
QAHistos.add("ptAntiDeuteron", "ptAntiDeuteron", {o2::framework::HistType::kTH1F, {ptAxis}});
QAHistos.add("etaAntideuteron", "etaAntideuteron", {o2::framework::HistType::kTH1F, {{100, -1.0f, 1.0f, "eta #bar{d}"}}});
QAHistos.add("hDCAxyVsPt-pre_selection", "DCAxy #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAxyAxis}});
QAHistos.add("hDCAzVsPt-pre-selection", "DCAz #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAzAxis}});

// processed events
hProcessedEvents->GetXaxis()->SetBinLabel(1, "Events processed");
hProcessedEvents->GetXaxis()->SetBinLabel(2, "TVX");
hProcessedEvents->GetXaxis()->SetBinLabel(3, "TF border");
hProcessedEvents->GetXaxis()->SetBinLabel(4, "ITS-RO border");
hProcessedEvents->GetXaxis()->SetBinLabel(5, "z-vertex");
hProcessedEvents->GetXaxis()->SetBinLabel(6, "Trigger");
}

// Tables
using CollisionCandidates = o2::soa::Join<o2::aod::Collisions, o2::aod::EvSels>;
using TrackCandidates = o2::soa::Join<o2::aod::Tracks, o2::aod::TracksCov, o2::aod::TracksExtra, o2::aod::TracksDCA, o2::aod::TrackSelection, o2::aod::pidTPCFullDe, o2::aod::pidTOFFullDe>;

// Single-Track Selection
template <typename T1>
bool passedSingleTrackSelection(const T1& track)
{
// Single-Track Selections
if (std::abs(track.eta()) > cfgEta)
return false;
if (!track.hasITS())
return false;
if (!track.hasTPC())
return false;
if (!track.hasTOF())
return false;
if (track.tpcNClsFound() < cfgTPCNclsFound)
return false;
if (track.tpcChi2NCl() > cfgTPCChi2Ncl)
return false;
if (track.itsChi2NCl() > cfgITSChi2Ncl)
return false;
if (track.itsNCls() < cfgITScls)
return false;
if (track.pt() > cfgMaxPt)
return false;
if (track.pt() < cfgMinPt)
return false;
if (track.sign() > 0)
return false;

return true;
}

/*void fillTriggerTable(bool keepEvent[])
{
table(keepEvent[0]);
}*/

o2::framework::Preslice<o2::aod::TrackAssoc> trackIndicesPerCollision = o2::aod::track_association::collisionId;
int mCurrentRun = -1;
void process(CollisionCandidates const& collisions,
o2::aod::TrackAssoc const& trackIndices,
TrackCandidates const& tracks,
o2::aod::BCsWithTimestamps const&)
{
for (const auto& collision : collisions) // start loop over collisions
{
LOG(debug) << "Processing collision with ID: " << collision.globalIndex();
if (mCurrentRun != collision.bc_as<o2::aod::BCsWithTimestamps>().runNumber()) { // If the run is new then we need to initialize the propagator field
o2::parameters::GRPMagField* grpo = ccdb->getForTimeStamp<o2::parameters::GRPMagField>("GLO/Config/GRPMagField", collision.bc_as<o2::aod::BCsWithTimestamps>().timestamp());
o2::base::Propagator::initFieldFromGRP(grpo);
mCurrentRun = collision.bc_as<o2::aod::BCsWithTimestamps>().runNumber();
}

// Is event good? keepEvent[0] = non promp deuteron
bool keepEvent[1]{}; // explicitly zero-initialised
hProcessedEvents->Fill(0.5);
// Event selection cuts
if (isTVX && !collision.selection_bit(o2::aod::evsel::kIsTriggerTVX)) {
// fillTriggerTable(keepEvent);
continue;
}
hProcessedEvents->Fill(1.5);
if (isNoTimeFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) {
// fillTriggerTable(keepEvent);
continue;
}
hProcessedEvents->Fill(2.5);
if (isNoITSROFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
// fillTriggerTable(keepEvent);
continue;
}
hProcessedEvents->Fill(3.5);
if (std::abs(collision.posZ()) > cutzvertex) {
// fillTriggerTable(keepEvent);
continue;
}
LOG(debug) << "Accepted collision!";
hProcessedEvents->Fill(4.5);
QAHistos.fill(HIST("hVtxZ"), collision.posZ());

// Loop over tracks
LOG(debug) << "Trying the slice";
const auto& trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex());
LOG(debug) << "Sliced!";
for (const auto& trackId : trackIdsThisCollision) { // start loop over tracks
LOG(debug) << "Processing track with ID: " << trackId.trackId();
const auto& track = tracks.rawIteratorAt(trackId.trackId());

std::array<float, 2> dca{track.dcaXY(), track.dcaZ()};
std::array<float, 3> pVec = track.pVector();

if (track.collisionId() != collision.globalIndex()) {
auto trackPar = getTrackParCov(track);
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackPar, 2.f, noMatCorr, &dca);
getPxPyPz(trackPar, pVec);
}

if (!passedSingleTrackSelection(track)) {
continue;
}

const bool isTOFDe = std::abs(track.tofNSigmaDe()) < cfgTOFNsigma;
const bool isTPCDe = std::abs(track.tpcNSigmaDe()) < cfgTPCNsigma;

if (isTPCDe && isTOFDe) {

QAHistos.fill(HIST("hDCAxyVsPt-pre_selection"), track.pt(), dca[0]);
QAHistos.fill(HIST("hDCAzVsPt-pre-selection"), track.pt(), dca[1]);
if (std::abs(dca[0]) < cfgDCAcut) {
continue;
}
keepEvent[0] = true;
QAHistos.fill(HIST("ptAntiDeuteron"), track.pt());
QAHistos.fill(HIST("etaAntideuteron"), track.eta());
QAHistos.fill(HIST("hDCAxyVsPt"), track.pt(), dca[0]);
QAHistos.fill(HIST("hDCAzVsPt"), track.pt(), dca[1]);
QAHistos.fill(HIST("hnSigmaTPCVsPt"), track.pt(), track.tpcNSigmaDe());
QAHistos.fill(HIST("hnSigmaTOFVsPt"), track.pt(), track.tofNSigmaDe());
}

} // end track loop
if (keepEvent[0]) {
hProcessedEvents->Fill(5.5);
}
// fillTriggerTable(keepEvent);
}
}
};

o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc)
{
return o2::framework::WorkflowSpec{o2::framework::adaptAnalysisTask<H2fromLbFilter>(cfgc)};
}
Loading