Skip to content
Merged
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
73 changes: 44 additions & 29 deletions EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -10,38 +10,52 @@
// or submit itself to any jurisdiction.
///
/// \file HeavyNeutralMesonFilter.cxx
///
/// \brief This code loops over collisions to filter events contaning heavy neutral mesons (omega or eta') using EMCal clusters and V0s (PCM)
///
/// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt; Maximilian Korwieser (maximilian.korwieser@cern.ch) - Technical University Munich
///

#include <vector>
#include <iterator>
#include <string>

#include "Math/GenVector/Boost.h"
#include "Math/Vector4D.h"
#include "TMath.h"
#include "TRandom3.h"

#include "EventFiltering/filterTables.h"
//
#include "PWGEM/PhotonMeson/Utils/HNMUtilities.h"
#include "PWGEM/PhotonMeson/Utils/PairUtilities.h"
#include "PWGJE/DataModel/EMCALMatchedCollisions.h"

#include "Common/DataModel/PIDResponse.h"
#include "Common/DataModel/PIDResponseITS.h"
#include "fairlogger/Logger.h"
#include "Framework/Configurable.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/CCDB/TriggerAliases.h"
#include "Common/Core/RecoDecay.h"
#include "Common/DataModel/EventSelection.h"
#include "Common/DataModel/Multiplicity.h"
#include "Common/DataModel/PIDResponseITS.h"
#include "Common/DataModel/PIDResponseTOF.h"
#include "Common/DataModel/PIDResponseTPC.h"
#include "Common/DataModel/TrackSelectionTables.h"

#include <CommonConstants/MathConstants.h>
#include <CommonConstants/PhysicsConstants.h>
#include <Framework/AnalysisDataModel.h>
#include <Framework/AnalysisHelpers.h>
#include <Framework/AnalysisTask.h>
#include <Framework/Array2D.h>
#include <Framework/Configurable.h>
#include <Framework/HistogramRegistry.h>
#include <Framework/HistogramSpec.h>
#include <Framework/InitContext.h>
#include <Framework/Logger.h>
#include <Framework/OutputObjHeader.h>
#include <Framework/runDataProcessing.h>

#include <Math/GenVector/Boost.h>
#include <Math/Vector4D.h> // IWYU pragma: keep
#include <Math/Vector4Dfwd.h>
#include <TH1.h>
#include <TString.h>

#include <sys/types.h>

#include <array>
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>

using namespace o2;
using namespace o2::framework;
Expand Down Expand Up @@ -90,21 +104,22 @@ const std::vector<std::string> pidCutsName{"TPC min", "TPC max", "TPCTOF max", "
const std::vector<std::string> femtoFilterNames{"PPOmega", "PPEtaPrime", "Omegad", "EtaPrimed", "OmegaP", "EtaPrimeP"};

// configs for tracks
const float pidcutsTable[kNFemtoPartners][kNPIDLimits]{
// these are need [[maybe_unused]] to silence a warning from clangd, since the compiler will inline them directly to the configs down below and then say: Variable 'X' is not needed and will not be emitted
[[maybe_unused]] const float pidcutsTable[kNFemtoPartners][kNPIDLimits]{
{-4.f, 4.f, 4.f, -99.f, 99.f},
{-4.f, 4.f, 4.f, -6.f, 6.f},
{-4.f, 4.f, 4.f, -99.f, 99.f}};

const float ptcutsTable[kNFemtoPartners][3]{
[[maybe_unused]] const float ptcutsTable[kNFemtoPartners][3]{
{0.35f, 6.f, 0.75f},
{0.55f, 2.f, 1.2f},
{0.35f, 6.f, 0.75f}};

const float nClusterMinTPC[1][kNFemtoPartners]{{80.0f, 80.0f, 80.0f}};
const float nClusterMinITS[1][kNFemtoPartners]{{4, 4, 4}};
[[maybe_unused]] const float nClusterMinTPC[1][kNFemtoPartners]{{80.0f, 80.0f, 80.0f}};
[[maybe_unused]] const float nClusterMinITS[1][kNFemtoPartners]{{4, 4, 4}};

static const float triggerSwitches[1][kNFemtoTriggers]{{1, 1, 1, 1, 1, 1}};
const float triggerLimits[1][kNFemtoTriggers]{{1.f, 1.f, 1.f, 1.f, 1.f, 1.f}};
[[maybe_unused]] static const float triggerSwitches[1][kNFemtoTriggers]{{1, 1, 1, 1, 1, 1}};
[[maybe_unused]] const float triggerLimits[1][kNFemtoTriggers]{{1.f, 1.f, 1.f, 1.f, 1.f, 1.f}};
} // namespace hnmtrigger

struct HeavyNeutralMesonFilter {
Expand Down
427 changes: 204 additions & 223 deletions PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h

Large diffs are not rendered by default.

439 changes: 210 additions & 229 deletions PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

Large diffs are not rendered by default.

6 changes: 6 additions & 0 deletions PWGEM/PhotonMeson/DataModel/gammaTables.h
Original file line number Diff line number Diff line change
Expand Up @@ -584,6 +584,12 @@ DECLARE_SOA_TABLE(SkimEMCClusters, "AOD", "SKIMEMCCLUSTER", //! table of skimmed
emccluster::DeltaEta, emccluster::TrackP, emccluster::TrackPt, emccluster::Pt<skimmedcluster::E, skimmedcluster::Eta>);
using SkimEMCCluster = SkimEMCClusters::iterator;

DECLARE_SOA_TABLE_VERSIONED(SkimEMCClusters_001, "AOD", "SKIMEMCCLUSTER", 1, //! table of skimmed EMCal clusters
o2::soa::Index<>, skimmedcluster::CollisionId, emccluster::Definition, skimmedcluster::E, skimmedcluster::Eta, skimmedcluster::Phi,
skimmedcluster::M02, skimmedcluster::NCells, skimmedcluster::Time, emccluster::IsExotic, emccluster::DeltaPhi,
emccluster::DeltaEta, emccluster::TrackP, emccluster::TrackPt, emccluster::Pt<skimmedcluster::E, skimmedcluster::Eta>);
using SkimEMCCluster = SkimEMCClusters::iterator;

DECLARE_SOA_TABLE(EMCEMEventIds, "AOD", "EMCEMEVENTID", emccluster::EMEventId); // To be joined with SkimEMCClusters table at analysis level.
// iterators
using EMCEMEventId = EMCEMEventIds::iterator;
Expand Down
90 changes: 53 additions & 37 deletions PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,10 @@
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
//
// ========================
// \file photonconversionbuilder.cxx
// \brief this task produces photon data table with KFParticle.
//
// \author Daiki Sekihata <daiki.sekihata@cern.ch>, Tokyo

/// \file photonconversionbuilder.cxx
/// \brief this task produces photon data table with KFParticle.
/// \author Daiki Sekihata <daiki.sekihata@cern.ch>, Tokyo

#include "PWGEM/Dilepton/Utils/PairUtilities.h"
#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
Expand All @@ -22,31 +21,42 @@

#include "Common/Core/RecoDecay.h"
#include "Common/Core/TPCVDriftManager.h"
#include "Common/Core/TableHelper.h"
#include "Common/Core/TrackSelection.h"
#include "Common/Core/trackUtilities.h"
#include "Common/DataModel/TrackSelectionTables.h"
#include "Common/DataModel/EventSelection.h"
#include "Common/DataModel/PIDResponseTPC.h"
#include "Tools/KFparticle/KFUtilities.h"

#include "CCDB/BasicCCDBManager.h"
#include "DataFormatsParameters/GRPMagField.h"
#include "DataFormatsParameters/GRPObject.h"
#include "DetectorsBase/GeometryManager.h"
#include "DetectorsBase/Propagator.h"
#include "Framework/ASoAHelpers.h"
#include "Framework/AnalysisDataModel.h"
#include "Framework/AnalysisTask.h"
#include "Framework/RunningWorkflowInfo.h"
#include "Framework/runDataProcessing.h"
#include "ReconstructionDataFormats/Track.h"

#include "Math/Vector4D.h"
#include <CCDB/BasicCCDBManager.h>
#include <CommonConstants/MathConstants.h>
#include <CommonConstants/PhysicsConstants.h>
#include <DataFormatsParameters/GRPMagField.h>
#include <DataFormatsParameters/GRPObject.h>
#include <DetectorsBase/GeometryManager.h>
#include <DetectorsBase/MatLayerCylSet.h>
#include <DetectorsBase/Propagator.h>
#include <Framework/AnalysisDataModel.h>
#include <Framework/AnalysisHelpers.h>
#include <Framework/AnalysisTask.h>
#include <Framework/Configurable.h>
#include <Framework/HistogramRegistry.h>
#include <Framework/HistogramSpec.h>
#include <Framework/InitContext.h>
#include <Framework/runDataProcessing.h>
#include <ReconstructionDataFormats/PID.h>

#include <Math/Vector4D.h> // IWYU pragma: keep
#include <Math/Vector4Dfwd.h>
#include <TPDGCode.h>

#include <KFPTrack.h>
#include <KFPVertex.h>
#include <KFParticle.h>

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iterator>
#include <map>
#include <set>
#include <string>
Expand All @@ -69,6 +79,12 @@ using MyCollisionsMC = soa::Join<MyCollisions, aod::McCollisionLabels>;
using MyTracksIU = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::pidTPCFullEl, aod::pidTPCFullPi>;
using MyTracksIUMC = soa::Join<MyTracksIU, aod::McTrackLabels, aod::mcTPCTuneOnData>;

enum MatCorrType {
None = 0,
TGeo = 1,
LUT = 2
};

struct PhotonConversionBuilder {
Produces<aod::V0PhotonsKF> v0photonskf;
Produces<aod::V0Legs> v0legs;
Expand Down Expand Up @@ -152,7 +168,7 @@ struct PhotonConversionBuilder {
{"V0/hConversionPointXY", "conversion point in XY;X (cm);Y (cm)", {HistType::kTH2F, {{400, -100.0f, 100.0f}, {400, -100.f, 100.f}}}},
{"V0/hConversionPointRZ", "conversion point in RZ;Z (cm);R_{xy} (cm)", {HistType::kTH2F, {{200, -100.0f, 100.0f}, {200, 0.f, 100.f}}}},
{"V0/hPt", "pT of V0 at PV;p_{T,#gamma} (GeV/c)", {HistType::kTH1F, {{1000, 0.0f, 10.0f}}}},
{"V0/hEtaPhi", "#eta vs. #varphi of V0 at PV;#varphi (rad.);#eta", {HistType::kTH2F, {{72, 0.0f, 2 * M_PI}, {200, -1, +1}}}},
{"V0/hEtaPhi", "#eta vs. #varphi of V0 at PV;#varphi (rad.);#eta", {HistType::kTH2F, {{72, 0.0f, o2::constants::math::TwoPI}, {200, -1, +1}}}},
{"V0/hCosPA", "cosine of pointing angle;cosine of pointing angle", {HistType::kTH1F, {{100, 0.99f, 1.f}}}},
{"V0/hCosPA_Rxy", "cosine of pointing angle;r_{xy} (cm);cosine of pointing angle", {HistType::kTH2F, {{200, 0, 100}, {100, 0.99f, 1.f}}}},
{"V0/hCosPAXY_Rxy", "cosine of pointing angle;r_{xy} (cm);cosine of pointing angle", {HistType::kTH2F, {{200, 0, 100}, {100, 0.99f, 1.f}}}},
Expand All @@ -170,7 +186,7 @@ struct PhotonConversionBuilder {
{"V0/hPCA_diffX", "PCA vs. trackiu X - R_{xy};distance btween 2 legs (cm);min trackiu X - R_{xy} (cm)", {HistType::kTH2F, {{500, 0.0f, 5.f}, {100, -50.0, 50.0f}}}},
{"V0/hPhiV", "#phi_{V}; #phi_{V} (rad.)", {HistType::kTH1F, {{500, 0.0f, 2 * M_PI}}}},
{"V0Leg/hPt", "pT of leg at SV;p_{T,e} (GeV/c)", {HistType::kTH1F, {{1000, 0.0f, 10.0f}}}},
{"V0Leg/hEtaPhi", "#eta vs. #varphi of leg at SV;#varphi (rad.);#eta", {HistType::kTH2F, {{72, 0.0f, 2 * M_PI}, {200, -1, +1}}}},
{"V0Leg/hEtaPhi", "#eta vs. #varphi of leg at SV;#varphi (rad.);#eta", {HistType::kTH2F, {{72, 0.0f, o2::constants::math::TwoPI}, {200, -1, +1}}}},
{"V0Leg/hRelDeltaPt", "pT resolution;p_{T} (GeV/c);#Deltap_{T}/p_{T}", {HistType::kTH2F, {{1000, 0.f, 10.f}, {100, 0, 1}}}},
{"V0Leg/hDCAxyz", "DCA xy vs. z to PV;DCA_{xy} (cm);DCA_{z} (cm)", {HistType::kTH2F, {{200, -50.f, 50.f}, {200, -50.f, +50.f}}}},
{"V0Leg/hdEdx_Pin", "TPC dE/dx vs. p_{in};p_{in} (GeV/c);TPC dE/dx", {HistType::kTH2F, {{1000, 0.f, 10.f}, {200, 0.f, 200.f}}}},
Expand All @@ -190,21 +206,21 @@ struct PhotonConversionBuilder {
ccdb->setLocalObjectValidityChecking();
ccdb->setFatalWhenNull(false);

if (useMatCorrType == 1) {
if (useMatCorrType == MatCorrType::TGeo) {
LOGF(info, "TGeo correction requested, loading geometry");
if (!o2::base::GeometryManager::isGeometryLoaded()) {
ccdb->get<TGeoManager>(geoPath);
}
}
if (useMatCorrType == 2) {
if (useMatCorrType == MatCorrType::LUT) {
LOGF(info, "LUT correction requested, loading LUT");
lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get<o2::base::MatLayerCylSet>(lutPath));
}

if (useMatCorrType == 1) {
if (useMatCorrType == MatCorrType::TGeo) {
matCorr = o2::base::Propagator::MatCorrType::USEMatCorrTGeo;
}
if (useMatCorrType == 2) {
if (useMatCorrType == MatCorrType::LUT) {
matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
}
}
Expand Down Expand Up @@ -467,8 +483,8 @@ struct PhotonConversionBuilder {

KFPTrack kfp_track_pos = createKFPTrackFromTrackParCov(pTrack, pos.sign(), pos.tpcNClsFound(), pos.tpcChi2NCl());
KFPTrack kfp_track_ele = createKFPTrackFromTrackParCov(nTrack, ele.sign(), ele.tpcNClsFound(), ele.tpcChi2NCl());
KFParticle kfp_pos(kfp_track_pos, -11);
KFParticle kfp_ele(kfp_track_ele, 11);
KFParticle kfp_pos(kfp_track_pos, kPositron);
KFParticle kfp_ele(kfp_track_ele, kElectron);
const KFParticle* GammaDaughters[2] = {&kfp_pos, &kfp_ele};

KFParticle gammaKF;
Expand Down Expand Up @@ -553,7 +569,7 @@ struct PhotonConversionBuilder {
gammaKF_PV.SetProductionVertex(KFPV);
float v0pt = RecoDecay::sqrtSumOfSquares(gammaKF_PV.GetPx(), gammaKF_PV.GetPy());
float v0eta = RecoDecay::eta(std::array{gammaKF_PV.GetPx(), gammaKF_PV.GetPy(), gammaKF_PV.GetPz()});
float v0phi = RecoDecay::phi(gammaKF_PV.GetPx(), gammaKF_PV.GetPy()) > 0.f ? RecoDecay::phi(gammaKF_PV.GetPx(), gammaKF_PV.GetPy()) : RecoDecay::phi(gammaKF_PV.GetPx(), gammaKF_PV.GetPy()) + TMath::TwoPi();
float v0phi = RecoDecay::constrainAngle(RecoDecay::phi(gammaKF_PV.GetPx(), gammaKF_PV.GetPy()));

// KFParticle gammaKF_DecayVtx2 = gammaKF;
// gammaKF_DecayVtx2.SetProductionVertex(KFPV);
Expand Down Expand Up @@ -657,18 +673,18 @@ struct PhotonConversionBuilder {
registry.fill(HIST("V0/hCosPAXY_Rxy"), rxy, cospaXY_kf);
registry.fill(HIST("V0/hCosPARZ_Rxy"), rxy, cospaRZ_kf);

for (auto& leg : {kfp_pos_DecayVtx, kfp_ele_DecayVtx}) {
for (const auto& leg : {kfp_pos_DecayVtx, kfp_ele_DecayVtx}) {
float legpt = RecoDecay::sqrtSumOfSquares(leg.GetPx(), leg.GetPy());
float legeta = RecoDecay::eta(std::array{leg.GetPx(), leg.GetPy(), leg.GetPz()});
float legphi = RecoDecay::phi(leg.GetPx(), leg.GetPy()) > 0.f ? RecoDecay::phi(leg.GetPx(), leg.GetPy()) : RecoDecay::phi(leg.GetPx(), leg.GetPy()) + TMath::TwoPi();
float legphi = RecoDecay::constrainAngle(RecoDecay::phi(leg.GetPx(), leg.GetPy()));
registry.fill(HIST("V0Leg/hPt"), legpt);
registry.fill(HIST("V0Leg/hEtaPhi"), legphi, legeta);
} // end of leg loop
for (auto& leg : {pos, ele}) {
for (const auto& leg : {pos, ele}) {
registry.fill(HIST("V0Leg/hdEdx_Pin"), leg.tpcInnerParam(), leg.tpcSignal());
registry.fill(HIST("V0Leg/hTPCNsigmaEl"), leg.tpcInnerParam(), leg.tpcNSigmaEl());
} // end of leg loop
for (auto& leg : {pTrack, nTrack}) {
for (const auto& leg : {pTrack, nTrack}) {
registry.fill(HIST("V0Leg/hXZ"), leg.getZ(), leg.getX());
registry.fill(HIST("V0Leg/hRelDeltaPt"), leg.getPt(), leg.getPt() * std::sqrt(leg.getSigma1Pt2()));
} // end of leg loop
Expand Down Expand Up @@ -793,7 +809,7 @@ struct PhotonConversionBuilder {
} // end of pca_map loop
// LOGF(info, "pca_map.size() = %d", pca_map.size());

for (auto& fullv0Id : stored_fullv0Ids) {
for (const auto& fullv0Id : stored_fullv0Ids) {
auto v0Id = std::get<0>(fullv0Id);
// auto collisionId = std::get<1>(fullv0Id);
// auto posId = std::get<2>(fullv0Id);
Expand All @@ -812,7 +828,7 @@ struct PhotonConversionBuilder {
fillV0Table<isMC, TBCs, TCollisions, TTracks>(v0, true);
} // end of fullv0Id loop

for (auto& collision : collisions) {
for (const auto& collision : collisions) {
if constexpr (isMC) {
if (!collision.has_mcCollision()) {
continue;
Expand Down
Loading
Loading