Skip to content
Merged
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
260 changes: 171 additions & 89 deletions PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.

Check failure on line 1 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/workflow-file]

Name of a workflow file must match the name of the main struct in it (without the PWG prefix). (Class implementation files should be in "Core" directories.)
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand All @@ -18,6 +18,10 @@
#include "Common/DataModel/EventSelection.h"
#include "Common/DataModel/TrackSelectionTables.h"

#include "CCDB/BasicCCDBManager.h"
#include "DataFormatsParameters/GRPMagField.h"
#include "DataFormatsParameters/GRPObject.h"
#include "DetectorsBase/Propagator.h"
#include "Framework/ASoA.h"
#include "Framework/ASoAHelpers.h"
#include "Framework/AnalysisDataModel.h"
Expand All @@ -30,13 +34,14 @@
#include "ReconstructionDataFormats/PID.h"
#include "ReconstructionDataFormats/Track.h"

#include <TDatabasePDG.h>

Check failure on line 37 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/database]

Do not use TDatabasePDG directly. Use o2::constants::physics::Mass... or Service<o2::framework::O2DatabasePDG> instead.
#include <TMath.h>
#include <TPDGCode.h>
#include <TRandom.h>
#include <TVector2.h>
#include <TVector3.h>

#include <string>
#include <vector>

using namespace std;
Expand All @@ -46,9 +51,20 @@
using namespace o2::constants::physics;
using std::array;

using MCTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::McTrackLabels>;
namespace
{
constexpr int DeuteronPDG = o2::constants::physics::kDeuteron;
constexpr int HePDG = o2::constants::physics::kHelium3;
constexpr int HypertritonPDG = o2::constants::physics::kHyperTriton;
constexpr int HyperHelium4PDG = o2::constants::physics::kHyperHelium4;
} // namespace

struct nucleiFromHypertritonMap {

using MCCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
using MCTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::McTrackLabels>;
Preslice<aod::McParticles> mMcParticlesPerCollision = o2::aod::mcparticle::mcCollisionId;

HistogramRegistry registryMC{
"registryMC",
{},
Expand All @@ -57,12 +73,12 @@
true};

// Track Parameters
Configurable<int> min_ITS_nClusters{"min_ITS_nClusters", 7, "minimum number of found ITS clusters"};

Check failure on line 76 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<int> min_ITS_InnerBarrel_nClusters{"min_ITS_InnerBarrel_nClusters", 1, "minimum number of found ITS Inner Barrel clusters"};

Check failure on line 77 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<int> max_ITS_InnerBarrel_nClusters{"max_ITS_InnerBarrel_nClusters", 3, "maximum number of found ITS Inner Barrel clusters"};

Check failure on line 78 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<int> min_TPC_nClusters{"min_TPC_nClusters", 100, "minimum number of found TPC clusters"};

Check failure on line 79 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<int> min_TPC_nCrossedRows{"min_TPC_nCrossedRows", 70, "minimum number of TPC crossed pad rows"};

Check failure on line 80 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<float> max_chi2_TPC{"max_chi2_TPC", 4.0f, "maximum TPC chi^2/Ncls"};

Check failure on line 81 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<float> min_chi2_TPC{"min_chi2_ITS", 0.5f, "minimum TPC chi^2/Ncls"};
Configurable<float> min_eta{"min_eta", -0.8f, "minimum_eta"};
Configurable<float> max_eta{"max_eta", +0.8f, "maximum_eta"};
Expand All @@ -70,119 +86,185 @@
Configurable<float> max_dcaz{"max_dcaz", 0.05f, "Maximum DCAz"};
Configurable<float> min_nsigmaTPC{"min_nsigmaTPC", -2.0f, "Minimum nsigma TPC"};
Configurable<float> max_nsigmaTPC{"max_nsigmaTPC", +2.0f, "Maximum nsigma TPC"};
Configurable<float> min_pt{"min_pt", 0.0f, "minimum pt of the tracks"};
Configurable<float> max_pt{"max_pt", 10.0f, "maximum pt of the tracks"};
Configurable<int> nbin_pt{"nbin_pt", 50, "number of pt bins"};
Configurable<int> nbin_dca = {"nbin_dca", 50, "number of DCA bins"};

Configurable<float> settingCutVertex{"settingCutVertex", 10.0f, "Accepted z-vertex range"};
ConfigurableAxis pt_axis{"pt_axis", {50, -10.f, 10.f}, ";signed #it{p}_{T} (GeV/#it{c});Counts"};

Configurable<bool> saveHelium{"saveHelium", false, "Save helium candidates"};

int AntideuteronPDG = -1000010020;
int AntihePDG = -1000020030;
int AntiHypertritonPDG = -1010010030;
int AntiHyperHelium4PDG = -1010020040;
Configurable<int> cfgMaterialCorrection{"cfgMaterialCorrection", static_cast<int>(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Type of material correction"};
Configurable<std::string> cfgCCDBurl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
Service<o2::ccdb::BasicCCDBManager> mCcdb;
int mRunNumber = 0;
float mBz = 0.f;

int mSelectedPDG = 0;
o2::dataformats::DCA mDcaInfoCov;
o2::track::TrackParametrizationWithError<float> mTrackParCov;
o2::base::MatLayerCylSet* mLut = nullptr;

void init(InitContext const&)
{
registryMC.add("hypertritonPtGen", "hypertritonPtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
mCcdb->setURL(cfgCCDBurl);
mCcdb->setCaching(true);
mCcdb->setLocalObjectValidityChecking();
mCcdb->setFatalWhenNull(false);
mLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));

registryMC.add("hypertritonPtGen", "hypertritonPtGen", HistType::kTH1F, {pt_axis});

if (saveHelium) {
registryMC.add("he3SecPtRec_from_hypertriton", "he3SecPtRec_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("he3SecPtGen_from_hypertriton", "he3SecPtGen_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("hyperHe4PtGen", "hyperHe4PtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("he3SecPtRec_from_hyperHe4", "he3SecPtRec_from_hyperHe4", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("he3SecPtGen_from_hyperHe4", "he3SecPtGen_from_hyperHe4", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("he3PtRec", "he3PtRec", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("he3PtGen", "he3PtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("he3SecPtRec_from_hypertriton", "he3SecPtRec_from_hypertriton", HistType::kTH1F, {pt_axis});
registryMC.add("he3SecPtGen_from_hypertriton", "he3SecPtGen_from_hypertriton", HistType::kTH1F, {pt_axis});
registryMC.add("hyperHe4PtGen", "hyperHe4PtGen", HistType::kTH1F, {pt_axis});
registryMC.add("he3SecPtRec_from_hyperHe4", "he3SecPtRec_from_hyperHe4", HistType::kTH1F, {pt_axis});
registryMC.add("he3SecPtGen_from_hyperHe4", "he3SecPtGen_from_hyperHe4", HistType::kTH1F, {pt_axis});
registryMC.add("he3PtRec", "he3PtRec", HistType::kTH1F, {pt_axis});
registryMC.add("he3PtGen", "he3PtGen", HistType::kTH1F, {pt_axis});
} else {
registryMC.add("deutSecPtRec_from_hypertriton", "deutSecPtRec_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("deutSecPtGen_from_hypertriton", "deutSecPtGen_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("deutPtRec", "deutPtRec", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("deutPtGen", "deutPtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
registryMC.add("deutSecPtRec_from_hypertriton", "deutSecPtRec_from_hypertriton", HistType::kTH1F, {pt_axis});
registryMC.add("deutSecPtGen_from_hypertriton", "deutSecPtGen_from_hypertriton", HistType::kTH1F, {pt_axis});
registryMC.add("deutPtRec", "deutPtRec", HistType::kTH1F, {pt_axis});
registryMC.add("deutPtGen", "deutPtGen", HistType::kTH1F, {pt_axis});
}
}

void processMC(const aod::McParticles& mcParticles, const MCTracks& tracks)
{
int selectedPDG = 0;
if (saveHelium) {
selectedPDG = AntihePDG;
mSelectedPDG = HePDG;
} else {
selectedPDG = AntideuteronPDG;
mSelectedPDG = DeuteronPDG;
}
}

void initCCDB(const aod::BCsWithTimestamps::iterator& bc)
{
if (mRunNumber == bc.runNumber())
return;

auto timestamp = bc.timestamp();
mRunNumber = bc.runNumber();

o2::parameters::GRPMagField* grpmag = mCcdb->getForTimeStamp<o2::parameters::GRPMagField>("GLO/Config/GRPMagField", timestamp);
o2::base::Propagator::initFieldFromGRP(grpmag);
o2::base::Propagator::Instance()->setMatLUT(mLut);
mBz = static_cast<float>(grpmag->getNominalL3Field());
LOGF(info, "Retrieved GRP for timestamp %ull (%i) with magnetic field of %1.2f kZG", timestamp, mRunNumber, mBz);
}

void processMC(const MCCollisions& collisions, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles, const MCTracks& tracks)
{
for (const auto& collision : collisions) {

if (!collision.sel8() || std::abs(collision.posZ()) > settingCutVertex) {
return;
}

auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
initCCDB(bc);

const int mcCollisionId = collision.mcCollisionId();
auto mcParticlesThisCollision = mcParticles.sliceBy(mMcParticlesPerCollision, mcCollisionId);
mcParticlesThisCollision.bindExternalIndices(&mcParticles);

for (const auto& mcparticle : mcParticlesThisCollision) {

if (((std::abs(mcparticle.pdgCode()) != HypertritonPDG && std::abs(mcparticle.pdgCode()) != HyperHelium4PDG) || !mcparticle.has_daughters()) && std::abs(mcparticle.pdgCode()) != mSelectedPDG) {
continue;
}

const float particleSign = mcparticle.pdgCode() < 0 ? -1. : 1.;

if (std::abs(mcparticle.pdgCode()) == HypertritonPDG) {

for (const auto& mcparticle : mcParticles) {
if (((mcparticle.pdgCode() == AntiHypertritonPDG || mcparticle.pdgCode() == AntiHyperHelium4PDG) && mcparticle.has_daughters()) || mcparticle.pdgCode() == selectedPDG) {
if (mcparticle.pdgCode() == AntiHypertritonPDG) {
for (auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
if (daughter.pdgCode() == selectedPDG) {
registryMC.fill(HIST("hypertritonPtGen"), mcparticle.pt());
if (saveHelium) {
registryMC.fill(HIST("he3SecPtGen_from_hypertriton"), daughter.pt());
} else {
registryMC.fill(HIST("deutSecPtGen_from_hypertriton"), daughter.pt());
}
for (const auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
if (std::abs(daughter.pdgCode()) != mSelectedPDG) {
continue;
}
const float daughterSign = daughter.pdgCode() < 0 ? -1. : 1.;

registryMC.fill(HIST("hypertritonPtGen"), particleSign * mcparticle.pt());
if (saveHelium) {
registryMC.fill(HIST("he3SecPtGen_from_hypertriton"), daughterSign * daughter.pt());
} else {
registryMC.fill(HIST("deutSecPtGen_from_hypertriton"), daughterSign * daughter.pt());
}
}
}
if (mcparticle.pdgCode() == AntiHyperHelium4PDG) {
for (auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
if (daughter.pdgCode() == selectedPDG) {
registryMC.fill(HIST("hyperHe4PtGen"), mcparticle.pt());
if (saveHelium) {
registryMC.fill(HIST("he3SecPtGen_from_hyperHe4"), daughter.pt());
}
} else if (std::abs(mcparticle.pdgCode()) == HyperHelium4PDG) {
for (const auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
if (std::abs(daughter.pdgCode()) != mSelectedPDG) {
continue;
}
const float daughterSign = daughter.pdgCode() < 0 ? -1. : 1.;

registryMC.fill(HIST("hyperHe4PtGen"), particleSign * mcparticle.pt());
if (saveHelium) {
registryMC.fill(HIST("he3SecPtGen_from_hyperHe4"), daughterSign * daughter.pt());
}
}
} else if (std::abs(mcparticle.pdgCode()) == HePDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("he3PtGen"), particleSign * mcparticle.pt());
} else if (std::abs(mcparticle.pdgCode()) == DeuteronPDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("deutPtGen"), particleSign * mcparticle.pt());
}
if (mcparticle.pdgCode() == AntihePDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("he3PtGen"), mcparticle.pt());
}

const o2::math_utils::Point3D<float> collisionVertex{collision.posX(), collision.posY(), collision.posZ()};

for (const auto& track : tracks) {

if (!track.has_mcParticle()) {
continue;
}
if (mcparticle.pdgCode() == AntideuteronPDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("deutPtGen"), mcparticle.pt());

const auto& mcparticle = track.mcParticle();
if (mcparticle.pdgCode() != mSelectedPDG) {
continue;
}
}
}

for (const auto& track : tracks) {
if (!track.has_mcParticle()) {
continue;
}
auto mcparticle = track.mcParticle();
if (mcparticle.pdgCode() != selectedPDG) {
continue;
}
mDcaInfoCov.set(999, 999, 999, 999, 999);
setTrackParCov(track, mTrackParCov);
mTrackParCov.setPID(track.pidForTracking());
std::array<float, 2> dcaInfo;
o2::base::Propagator::Instance()->propagateToDCA(collisionVertex, mTrackParCov, mBz, 2.f, static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value), &dcaInfo);

if (track.itsNCls() < min_ITS_nClusters ||
track.itsNClsInnerBarrel() < min_ITS_InnerBarrel_nClusters ||
track.itsNClsInnerBarrel() > max_ITS_InnerBarrel_nClusters ||
track.tpcNClsFound() < min_TPC_nClusters ||
track.tpcNClsCrossedRows() < min_TPC_nCrossedRows ||
track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() ||
track.tpcChi2NCl() > 4.f ||
track.tpcChi2NCl() < min_chi2_TPC ||
track.eta() < min_eta || track.eta() > max_eta ||
track.dcaXY() > max_dcaxy || track.dcaXY() < -max_dcaxy ||
track.dcaZ() > max_dcaz || track.dcaZ() < -max_dcaz ||
track.itsChi2NCl() > 36.f) {
continue;
}
if (mcparticle.pdgCode() == AntideuteronPDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("deutPtRec"), track.pt());
}
if (mcparticle.pdgCode() == AntihePDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("he3PtRec"), 2 * track.pt());
}
const float dcaXY = dcaInfo[0];
const float dcaZ = dcaInfo[1];

if (track.itsNCls() < min_ITS_nClusters ||
track.itsNClsInnerBarrel() < min_ITS_InnerBarrel_nClusters ||
track.itsNClsInnerBarrel() > max_ITS_InnerBarrel_nClusters ||
track.tpcNClsFound() < min_TPC_nClusters ||
track.tpcNClsCrossedRows() < min_TPC_nCrossedRows ||
track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() ||

Check failure on line 237 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
track.tpcChi2NCl() > max_chi2_TPC ||
track.tpcChi2NCl() < min_chi2_TPC ||
track.eta() < min_eta || track.eta() > max_eta ||
dcaXY > max_dcaxy || dcaXY < -max_dcaxy ||
dcaZ > max_dcaz || dcaZ < -max_dcaz ||
track.itsChi2NCl() > 36.f) {

Check failure on line 243 in PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
continue;
}

const float particleSign = mcparticle.pdgCode() < 0 ? -1. : 1.;

for (auto& motherparticle : mcparticle.mothers_as<aod::McParticles>()) {
if (motherparticle.pdgCode() == AntiHypertritonPDG || motherparticle.pdgCode() == AntiHyperHelium4PDG) {
if (motherparticle.pdgCode() == AntiHypertritonPDG) {
if (mcparticle.pdgCode() == AntihePDG) {
registryMC.fill(HIST("he3SecPtRec_from_hypertriton"), 2 * track.pt());
if (std::abs(mcparticle.pdgCode()) == DeuteronPDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("deutPtRec"), particleSign * track.pt());
}
if (std::abs(mcparticle.pdgCode()) == HePDG && mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("he3PtRec"), particleSign * 2 * track.pt());
}

for (const auto& motherparticle : mcparticle.mothers_as<aod::McParticles>()) {

if (std::abs(motherparticle.pdgCode()) == HypertritonPDG) {
if (std::abs(mcparticle.pdgCode()) == HePDG) {
registryMC.fill(HIST("he3SecPtRec_from_hypertriton"), 2 * particleSign * track.pt());
} else {
registryMC.fill(HIST("deutSecPtRec_from_hypertriton"), track.pt());
registryMC.fill(HIST("deutSecPtRec_from_hypertriton"), particleSign * track.pt());
}
} else if (std::abs(motherparticle.pdgCode()) == HyperHelium4PDG) {
if (std::abs(mcparticle.pdgCode()) == HePDG) {
registryMC.fill(HIST("he3SecPtRec_from_hyperHe4"), 2 * particleSign * track.pt());
}
}
if (motherparticle.pdgCode() == AntiHyperHelium4PDG) {
registryMC.fill(HIST("he3SecPtRec_from_hyperHe4"), 2 * track.pt());
}
}
}
Expand Down
Loading