diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index 2a2d644e7c4..136af7ab61c 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -60,6 +60,7 @@ struct TreeWriterTpcV0 { using Colls = soa::Join; using MyBCTable = soa::Join; using V0sWithID = soa::Join; + using CascsWithID = soa::Join; /// Tables to be produced Produces rowTPCTree; @@ -76,13 +77,26 @@ struct TreeWriterTpcV0 { Configurable dwnSmplFactor_Pi{"dwnSmplFactor_Pi", 1., "downsampling factor for pions, default fraction to keep is 1."}; Configurable dwnSmplFactor_Pr{"dwnSmplFactor_Pr", 1., "downsampling factor for protons, default fraction to keep is 1."}; Configurable dwnSmplFactor_El{"dwnSmplFactor_El", 1., "downsampling factor for electrons, default fraction to keep is 1."}; + Configurable dwnSmplFactor_Ka{"dwnSmplFactor_Ka", 1., "downsampling factor for kaons, default fraction to keep is 1."}; Configurable sqrtSNN{"sqrt_s_NN", 0., "sqrt(s_NN), used for downsampling with the Tsallis distribution"}; Configurable downsamplingTsalisPions{"downsamplingTsalisPions", -1., "Downsampling factor to reduce the number of pions"}; Configurable downsamplingTsalisProtons{"downsamplingTsalisProtons", -1., "Downsampling factor to reduce the number of protons"}; Configurable downsamplingTsalisElectrons{"downsamplingTsalisElectrons", -1., "Downsampling factor to reduce the number of electrons"}; + Configurable downsamplingTsalisKaons{"downsamplingTsalisKaons", -1., "Downsampling factor to reduce the number of kaons"}; Configurable maxPt4dwnsmplTsalisPions{"maxPt4dwnsmplTsalisPions", 100., "Maximum Pt for applying downsampling factor of pions"}; Configurable maxPt4dwnsmplTsalisProtons{"maxPt4dwnsmplTsalisProtons", 100., "Maximum Pt for applying downsampling factor of protons"}; Configurable maxPt4dwnsmplTsalisElectrons{"maxPt4dwnsmplTsalisElectrons", 100., "Maximum Pt for applying downsampling factor of electrons"}; + Configurable maxPt4dwnsmplTsalisKaons{"maxPt4dwnsmplTsalisKaons", 100., "Maximum Pt for applying downsampling factor of kaons"}; + + enum { // Reconstructed V0 + kUndef = -1, + kGamma = 0, + kK0S = 1, + kLambda = 2, + kAntiLambda = 3, + kOmega = 4, + kAntiOmega = 5 + }; Filter trackFilter = (trackSelection.node() == 0) || ((trackSelection.node() == 1) && requireGlobalTrackInFilter()) || @@ -94,8 +108,8 @@ struct TreeWriterTpcV0 { ctpRateFetcher mRateFetcher; /// Funktion to fill skimmed tables - template - void fillSkimmedV0Table(V0 const& v0, T const& track, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate) + template + void fillSkimmedV0Table(V0Casc const& v0casc, T const& track, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate) { const double ncl = track.tpcNClsFound(); @@ -107,12 +121,12 @@ struct TreeWriterTpcV0 { auto trackocc = collision.trackOccupancyInTimeRange(); auto ft0occ = collision.ft0cOccupancyInTimeRange(); - const float alpha = v0.alpha(); - const float qt = v0.qtarm(); - const float cosPA = v0.v0cosPA(); - const float pT = v0.pt(); - const float v0radius = v0.v0radius(); - const float gammapsipair = v0.psipair(); + const float alpha = v0casc.alpha(); + const float qt = v0casc.qtarm(); + const float cosPA = GetCosPA(v0casc, collision); + const float pT = v0casc.pt(); + const float v0radius = GetRadius(v0casc); + const float gammapsipair = v0casc.psipair(); const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { @@ -151,8 +165,8 @@ struct TreeWriterTpcV0 { } }; - template - void fillSkimmedV0TableWithdEdxTrQA(V0 const& v0, T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate) + template + void fillSkimmedV0TableWithdEdxTrQA(V0Casc const& v0casc, T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate) { const double ncl = track.tpcNClsFound(); @@ -164,12 +178,12 @@ struct TreeWriterTpcV0 { auto trackocc = collision.trackOccupancyInTimeRange(); auto ft0occ = collision.ft0cOccupancyInTimeRange(); - const float alpha = v0.alpha(); - const float qt = v0.qtarm(); - const float cosPA = v0.v0cosPA(); - const float pT = v0.pt(); - const float v0radius = v0.v0radius(); - const float gammapsipair = v0.psipair(); + const float alpha = v0casc.alpha(); + const float qt = v0casc.qtarm(); + const float cosPA = GetCosPA(v0casc, collision); + const float pT = v0casc.pt(); + const float v0radius = GetRadius(v0casc); + const float gammapsipair = v0casc.psipair(); const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { @@ -210,8 +224,8 @@ struct TreeWriterTpcV0 { }; /// Function to fill skimmed tables - template - void fillSkimmedV0TableWithTrQA(V0 const& v0, T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate, int bcGlobalIndex, int bcTimeFrameId, int bcBcInTimeFrame) + template + void fillSkimmedV0TableWithTrQA(V0Casc const& v0casc, T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate, int bcGlobalIndex, int bcTimeFrameId, int bcBcInTimeFrame) { const double ncl = track.tpcNClsFound(); @@ -223,12 +237,12 @@ struct TreeWriterTpcV0 { auto trackocc = collision.trackOccupancyInTimeRange(); auto ft0occ = collision.ft0cOccupancyInTimeRange(); - const float alpha = v0.alpha(); - const float qt = v0.qtarm(); - const float cosPA = v0.v0cosPA(); - const float pT = v0.pt(); - const float v0radius = v0.v0radius(); - const float gammapsipair = v0.psipair(); + const float alpha = v0casc.alpha(); + const float qt = v0casc.qtarm(); + const float cosPA = GetCosPA(v0casc, collision); + const float pT = v0casc.pt(); + const float v0radius = GetRadius(v0casc); + const float gammapsipair = v0casc.psipair(); const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { @@ -335,8 +349,34 @@ struct TreeWriterTpcV0 { ccdb->setFatalWhenNull(false); } + /// Evaluate cosPA of the v0 + template + double GetCosPA(V0sWithID::iterator const& v0, CollisionType const&) + { + return v0.v0cosPA(); + } + + /// Evaluate cosPA of the cascade + template + double GetCosPA(CascsWithID::iterator const& casc, CollisionType const& collision) + { + return casc.casccosPA(collision.posX(), collision.posY(), collision.posZ()); + } + + /// Evaluate radius of the v0 + double GetRadius(V0sWithID::iterator const& v0) + { + return v0.v0radius(); + } + + /// Evaluate radius of the cascade + double GetRadius(CascsWithID::iterator const& casc) + { + return casc.cascradius(); + } + /// Apply a track quality selection with a filter! - void processStandard(Colls::iterator const& collision, soa::Filtered const& tracks, V0sWithID const& v0s, aod::BCsWithTimestamps const&) + void processStandard(Colls::iterator const& collision, soa::Filtered const& tracks, V0sWithID const& v0s, CascsWithID const& cascs, aod::BCsWithTimestamps const&) { /// Check event slection if (!isEventSelected(collision, tracks)) { @@ -396,10 +436,25 @@ struct TreeWriterTpcV0 { } } } + + /// Loop over cascade candidates + for (const auto& casc : cascs) { + auto bachTrack = casc.bachelor_as>(); + if (casc.cascaddid() == kUndef) { + continue; + } + // Omega and antiomega + if (static_cast(bachTrack.pidbit() & (1 << kOmega)) || static_cast(bachTrack.pidbit() & (1 << kAntiOmega))) { + if (downsampleTsalisCharged(bachTrack.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], maxPt4dwnsmplTsalisKaons)) { + fillSkimmedV0Table(casc, bachTrack, collision, bachTrack.tpcNSigmaKa(), bachTrack.tofNSigmaKa(), bachTrack.tpcExpSignalKa(bachTrack.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } + } + } + } /// process Standard PROCESS_SWITCH(TreeWriterTpcV0, processStandard, "Standard V0 Samples for PID", true); - void processStandardWithCorrecteddEdx(Colls::iterator const& collision, soa::Filtered const& tracks, V0sWithID const& v0s, aod::BCsWithTimestamps const&) + void processStandardWithCorrecteddEdx(Colls::iterator const& collision, soa::Filtered const& tracks, V0sWithID const& v0s, CascsWithID const& cascs, aod::BCsWithTimestamps const&) { /// Check event slection if (!isEventSelected(collision, tracks)) { @@ -459,12 +514,27 @@ struct TreeWriterTpcV0 { } } } - } /// process Standard + + /// Loop over cascade candidates + for (const auto& casc : cascs) { + auto bachTrack = casc.bachelor_as>(); + if (casc.cascaddid() == kUndef) { + continue; + } + // Omega and antiomega + if (static_cast(bachTrack.pidbit() & (1 << kOmega)) || static_cast(bachTrack.pidbit() & (1 << kAntiOmega))) { + if (downsampleTsalisCharged(bachTrack.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], maxPt4dwnsmplTsalisKaons)) { + fillSkimmedV0Table(casc, bachTrack, collision, bachTrack.tpcNSigmaKa(), bachTrack.tofNSigmaKa(), bachTrack.tpcExpSignalKa(bachTrack.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } + } + } + } /// process StandardWithCorrecteddEdx PROCESS_SWITCH(TreeWriterTpcV0, processStandardWithCorrecteddEdx, "Standard V0 Samples for PID with corrected dEdx", false); Preslice perCollisionTracks = aod::track::collisionId; Preslice perCollisionV0s = aod::v0data::collisionId; - void processWithdEdxTrQA(Colls const& collisions, Trks const& myTracks, V0sWithID const& myV0s, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) + Preslice perCollisionCascs = aod::cascdata::collisionId; + void processWithdEdxTrQA(Colls const& collisions, Trks const& myTracks, V0sWithID const& myV0s, CascsWithID const& myCascs, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) { std::vector labelTrack2TrackQA; labelTrack2TrackQA.clear(); @@ -477,6 +547,7 @@ struct TreeWriterTpcV0 { for (const auto& collision : collisions) { auto tracks = myTracks.sliceBy(perCollisionTracks, collision.globalIndex()); auto v0s = myV0s.sliceBy(perCollisionV0s, collision.globalIndex()); + auto cascs = myCascs.sliceBy(perCollisionCascs, collision.globalIndex()); /// Check event slection if (!isEventSelected(collision, tracks)) { continue; @@ -552,12 +623,37 @@ struct TreeWriterTpcV0 { } } } + + /// Loop over cascade candidates + for (const auto& casc : cascs) { + auto bachTrack = casc.bachelor_as(); + if (casc.cascaddid() == kUndef) { + continue; + } + + aod::TracksQA bachTrackQA; + bool existBachTrkQA; + if (labelTrack2TrackQA[bachTrack.globalIndex()] != -1) { + bachTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[bachTrack.globalIndex()]); + existBachTrkQA = true; + } else { + bachTrackQA = tracksQA.iteratorAt(0); + existBachTrkQA = false; + } + + // Omega and antiomega + if (static_cast(bachTrack.pidbit() & (1 << kOmega)) || static_cast(bachTrack.pidbit() & (1 << kAntiOmega))) { + if (downsampleTsalisCharged(bachTrack.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], maxPt4dwnsmplTsalisKaons)) { + fillSkimmedV0TableWithdEdxTrQA(casc, bachTrack, bachTrackQA, existBachTrkQA, collision, bachTrack.tpcNSigmaKa(), bachTrack.tofNSigmaKa(), bachTrack.tpcExpSignalKa(bachTrack.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } + } + } } } /// process with dEdx from TrackQA PROCESS_SWITCH(TreeWriterTpcV0, processWithdEdxTrQA, "Standard V0 Samples with dEdx from Track QA for PID", false); Preslice perCollisionTracksWithNewDEdx = aod::track::collisionId; - void processWithdEdxTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, V0sWithID const& myV0s, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) + void processWithdEdxTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, V0sWithID const& myV0s, CascsWithID const& myCascs, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) { std::vector labelTrack2TrackQA; labelTrack2TrackQA.clear(); @@ -570,6 +666,7 @@ struct TreeWriterTpcV0 { for (const auto& collision : collisions) { auto tracks = myTracks.sliceBy(perCollisionTracksWithNewDEdx, collision.globalIndex()); auto v0s = myV0s.sliceBy(perCollisionV0s, collision.globalIndex()); + auto cascs = myCascs.sliceBy(perCollisionCascs, collision.globalIndex()); /// Check event slection if (!isEventSelected(collision, tracks)) { continue; @@ -645,11 +742,36 @@ struct TreeWriterTpcV0 { } } } + + /// Loop over cascade candidates + for (const auto& casc : cascs) { + auto bachTrack = casc.bachelor_as(); + if (casc.cascaddid() == kUndef) { + continue; + } + + aod::TracksQA bachTrackQA; + bool existBachTrkQA; + if (labelTrack2TrackQA[bachTrack.globalIndex()] != -1) { + bachTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[bachTrack.globalIndex()]); + existBachTrkQA = true; + } else { + bachTrackQA = tracksQA.iteratorAt(0); + existBachTrkQA = false; + } + + // Omega and antiomega + if (static_cast(bachTrack.pidbit() & (1 << kOmega)) || static_cast(bachTrack.pidbit() & (1 << kAntiOmega))) { + if (downsampleTsalisCharged(bachTrack.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], maxPt4dwnsmplTsalisKaons)) { + fillSkimmedV0TableWithdEdxTrQA(casc, bachTrack, bachTrackQA, existBachTrkQA, collision, bachTrack.tpcNSigmaKa(), bachTrack.tofNSigmaKa(), bachTrack.tpcExpSignalKa(bachTrack.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } + } + } } } /// process with dEdx from TrackQA PROCESS_SWITCH(TreeWriterTpcV0, processWithdEdxTrQAWithCorrecteddEdx, "Standard V0 Samples with dEdx from Track QA for PID with corrected dEdx", false); - void processWithTrQA(Colls const& collisions, Trks const& myTracks, V0sWithID const& myV0s, MyBCTable const&, aod::TracksQAVersion const& tracksQA) + void processWithTrQA(Colls const& collisions, Trks const& myTracks, V0sWithID const& myV0s, CascsWithID const& myCascs, MyBCTable const&, aod::TracksQAVersion const& tracksQA) { std::vector labelTrack2TrackQA; labelTrack2TrackQA.clear(); @@ -662,6 +784,7 @@ struct TreeWriterTpcV0 { for (const auto& collision : collisions) { auto tracks = myTracks.sliceBy(perCollisionTracks, collision.globalIndex()); auto v0s = myV0s.sliceBy(perCollisionV0s, collision.globalIndex()); + auto cascs = myCascs.sliceBy(perCollisionCascs, collision.globalIndex()); /// Check event slection if (!isEventSelected(collision, tracks)) { continue; @@ -740,11 +863,36 @@ struct TreeWriterTpcV0 { } } } + + /// Loop over cascade candidates + for (const auto& casc : cascs) { + auto bachTrack = casc.bachelor_as(); + if (casc.cascaddid() == kUndef) { + continue; + } + + aod::TracksQA bachTrackQA; + bool existBachTrkQA; + if (labelTrack2TrackQA[bachTrack.globalIndex()] != -1) { + bachTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[bachTrack.globalIndex()]); + existBachTrkQA = true; + } else { + bachTrackQA = tracksQA.iteratorAt(0); + existBachTrkQA = false; + } + + // Omega and antiomega + if (static_cast(bachTrack.pidbit() & (1 << kOmega)) || static_cast(bachTrack.pidbit() & (1 << kAntiOmega))) { + if (downsampleTsalisCharged(bachTrack.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], maxPt4dwnsmplTsalisKaons)) { + fillSkimmedV0TableWithTrQA(casc, bachTrack, bachTrackQA, existBachTrkQA, collision, bachTrack.tpcNSigmaKa(), bachTrack.tofNSigmaKa(), bachTrack.tpcExpSignalKa(bachTrack.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } + } } } /// process with TrackQA PROCESS_SWITCH(TreeWriterTpcV0, processWithTrQA, "Standard V0 Samples with Track QA for PID", false); - void processWithTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, V0sWithID const& myV0s, MyBCTable const&, aod::TracksQAVersion const& tracksQA) + void processWithTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, V0sWithID const& myV0s, CascsWithID const& myCascs, MyBCTable const&, aod::TracksQAVersion const& tracksQA) { std::vector labelTrack2TrackQA; labelTrack2TrackQA.clear(); @@ -757,6 +905,7 @@ struct TreeWriterTpcV0 { for (const auto& collision : collisions) { auto tracks = myTracks.sliceBy(perCollisionTracksWithNewDEdx, collision.globalIndex()); auto v0s = myV0s.sliceBy(perCollisionV0s, collision.globalIndex()); + auto cascs = myCascs.sliceBy(perCollisionCascs, collision.globalIndex()); /// Check event slection if (!isEventSelected(collision, tracks)) { continue; @@ -835,6 +984,31 @@ struct TreeWriterTpcV0 { } } } + + /// Loop over cascade candidates + for (const auto& casc : cascs) { + auto bachTrack = casc.bachelor_as(); + if (casc.cascaddid() == kUndef) { + continue; + } + + aod::TracksQA bachTrackQA; + bool existBachTrkQA; + if (labelTrack2TrackQA[bachTrack.globalIndex()] != -1) { + bachTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[bachTrack.globalIndex()]); + existBachTrkQA = true; + } else { + bachTrackQA = tracksQA.iteratorAt(0); + existBachTrkQA = false; + } + + // Omega and antiomega + if (static_cast(bachTrack.pidbit() & (1 << kOmega)) || static_cast(bachTrack.pidbit() & (1 << kAntiOmega))) { + if (downsampleTsalisCharged(bachTrack.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], maxPt4dwnsmplTsalisKaons)) { + fillSkimmedV0TableWithTrQA(casc, bachTrack, bachTrackQA, existBachTrkQA, collision, bachTrack.tpcNSigmaKa(), bachTrack.tofNSigmaKa(), bachTrack.tpcExpSignalKa(bachTrack.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } + } } } /// process with TrackQA PROCESS_SWITCH(TreeWriterTpcV0, processWithTrQAWithCorrecteddEdx, "Standard V0 Samples with Track QA for PID with corrected dEdx", false); diff --git a/PWGDQ/DataModel/ReducedInfoTables.h b/PWGDQ/DataModel/ReducedInfoTables.h index 244a1f7417c..dd92f5a9c76 100644 --- a/PWGDQ/DataModel/ReducedInfoTables.h +++ b/PWGDQ/DataModel/ReducedInfoTables.h @@ -1124,6 +1124,14 @@ DECLARE_SOA_COLUMN(V0AddID, v0addid, int8_t); //! DECLARE_SOA_TABLE(V0MapID, "AOD", "V0MAPID", //! v0mapID::V0AddID); +namespace cascmapID +{ +DECLARE_SOA_COLUMN(CascAddID, cascaddid, int8_t); //! +} // namespace cascmapID + +DECLARE_SOA_TABLE(CascMapID, "AOD", "CASCMAPID", //! + cascmapID::CascAddID); + namespace DalBits { DECLARE_SOA_COLUMN(DALITZBits, dalitzBits, uint8_t); //! diff --git a/PWGDQ/Tasks/v0selector.cxx b/PWGDQ/Tasks/v0selector.cxx index 0c970d80c91..1af3fb2cf0e 100644 --- a/PWGDQ/Tasks/v0selector.cxx +++ b/PWGDQ/Tasks/v0selector.cxx @@ -17,35 +17,40 @@ // Comments, questions, complaints, suggestions? // Please write to: daiki.sekihata@cern.ch // -#include -#include -#include -#include -#include -#include "Math/Vector4D.h" -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoAHelpers.h" -#include "ReconstructionDataFormats/Track.h" -#include "Common/Core/trackUtilities.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "Common/DataModel/EventSelection.h" +#include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" -#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" #include "DCAFitter/DCAFitterN.h" -#include "PWGDQ/DataModel/ReducedInfoTables.h" -#include "PWGDQ/Core/VarManager.h" -#include "PWGDQ/Core/HistogramManager.h" -#include "PWGDQ/Core/HistogramsLibrary.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsBase/GeometryManager.h" -#include "DataFormatsParameters/GRPObject.h" #include "DataFormatsParameters/GRPMagField.h" -#include "CCDB/BasicCCDBManager.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/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include "Math/Vector4D.h" + +#include +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; @@ -85,7 +90,24 @@ struct v0selector { Configurable cutAPL1{"cutAPL1", 0.107, "cutAPL1"}; Configurable cutAPL2{"cutAPL2", -0.69, "cutAPL2"}; Configurable cutAPL3{"cutAPL3", 0.5, "cutAPL3"}; + // Omega & A-Omega cuts + Configurable cutAPOmegaUp1{"cutAPOmegaUp1", 0.25, "cutAPOmegaUp1"}; + Configurable cutAPOmegaUp2{"cutAPOmegaUp2", 0.358, "cutAPOmegaUp2"}; + Configurable cutAPOmegaUp3{"cutAPOmegaUp3", 0.35, "cutAPOmegaUp3"}; + Configurable cutAPOmegaDown1{"cutAPOmegaDown1", 0.15, "cutAPOmegaDown1"}; + Configurable cutAPOmegaDown2{"cutAPOmegaDown2", 0.358, "cutAPOmegaDown2"}; + Configurable cutAPOmegaDown3{"cutAPOmegaDown3", 0.16, "cutAPOmegaDown3"}; + Configurable cutAlphaOmegaHigh{"cutAlphaOmegaHigh", 0.358, "cutAlphaOmegaHigh"}; + Configurable cutAlphaOmegaLow{"cutAlphaOmegaLow", 0., "cutAlphaOmegaLow"}; + Configurable cutQTOmegaLowOuterArc{"cutQTOmegaLowOuterArc", 0.14, "cutQTOmegaLowOuterArc"}; + Configurable cutMassOmegaHigh{"cutMassOmegaHigh", 1.677, "cutMassOmegaHigh"}; + Configurable cutMassOmegaLow{"cutMassOmegaLow", 1.667, "cutMassOmegaLow"}; + Configurable cutMassCascV0Low{"cutMassCascV0Low", 1.110, "cutMassCascV0Low"}; + Configurable cutMassCascV0High{"cutMassCascV0High", 1.120, "cutMassCascV0High"}; + Configurable produceV0ID{"produceV0ID", false, "Produce additional V0ID table"}; + Configurable selectCascades{"selectCascades", false, "Select cascades in addition to v0s"}; + Configurable produceCascID{"produceCascID", false, "Produce additional CascID table"}; enum { // Reconstructed V0 kUndef = -1, @@ -93,11 +115,13 @@ struct v0selector { kK0S = 1, kLambda = 2, kAntiLambda = 3, - kOmega = 4 + kOmega = 4, + kAntiOmega = 5 }; Produces v0bits; Produces v0mapID; + Produces cascmapID; // int checkV0(const array& ppos, const array& pneg) int checkV0(const float alpha, const float qt) @@ -154,6 +178,29 @@ struct v0selector { return kUndef; } + int checkCascade(float alpha, float qt) + { + const bool isAlphaPos = alpha > 0; + alpha = std::fabs(alpha); + + const float qUp = std::abs(alpha - cutAPOmegaUp2) > std::abs(cutAPOmegaUp3) ? 0. : cutAPOmegaUp1 * std::sqrt(1.0f - ((alpha - cutAPOmegaUp2) * (alpha - cutAPOmegaUp2)) / (cutAPOmegaUp3 * cutAPOmegaUp3)); + const float qDown = std::abs(alpha - cutAPOmegaDown2) > std::abs(cutAPOmegaDown3) ? 0. : cutAPOmegaDown1 * std::sqrt(1.0f - ((alpha - cutAPOmegaDown2) * (alpha - cutAPOmegaDown2)) / (cutAPOmegaDown3 * cutAPOmegaDown3)); + + if (alpha < cutAlphaOmegaLow || alpha > cutAlphaOmegaHigh || qt < qDown || qt > qUp) { + return kUndef; + } + + if (alpha > cutAPOmegaUp2 && qt < cutQTOmegaLowOuterArc) { + return kUndef; + } + + if (isAlphaPos) { + return kOmega; + } else { + return kAntiOmega; + } + } + // Configurables Configurable v0max_mee{"v0max_mee", 0.1, "max mee for photon"}; Configurable maxpsipair{"maxpsipair", 1.6, "max psi_pair for photon"}; @@ -166,16 +213,28 @@ struct v0selector { Configurable mincrossedrows{"mincrossedrows", 70, "min crossed rows"}; Configurable maxchi2tpc{"maxchi2tpc", 4.0, "max chi2/NclsTPC"}; Configurable fillhisto{"fillhisto", false, "flag to fill histograms"}; - // cutNsigmaElTPC, cutNsigmaPiTPC, cutNsigmaPrTPC + // Cascade-related Configurables + Configurable cascDcaMax{"cascDcaMax", 0.3, "DCA cascade daughters"}; + Configurable cascV0DcaMax{"cascV0DcaMax", 0.3, "DCA V0 daughters of the cascade"}; + Configurable cascRadiusMin{"cascRadiusMin", 0.0, "Cascade Radius min"}; + Configurable cascRadiusMax{"cascRadiusMax", 90.0, "Cascade Radius max"}; + Configurable cascV0RadiusMin{"cascV0RadiusMin", 0.0, "V0 of the Cascade Radius min"}; + Configurable cascV0RadiusMax{"cascV0RadiusMax", 90.0, "V0 of the Cascade Radius max"}; + Configurable cascCosinePAMin{"cascCosinePAMin", 0.998, "Cascade CosPA min"}; + Configurable cascV0CosinePAMin{"cascV0CosinePAMin", 0.995, "V0 of the Cascade CosPA min"}; + Configurable cascV0CosinePAMax{"cascV0CosinePAMax", 1.000, "V0 of the Cascade CosPA max"}; + // cutNsigmaElTPC, cutNsigmaPiTPC, cutNsigmaPrTPC, cutNsigmaKaTPC Configurable cutNsigmaElTPC{"cutNsigmaElTPC", 5.0, "cutNsigmaElTPC"}; Configurable cutNsigmaPiTPC{"cutNsigmaPiTPC", 5.0, "cutNsigmaPiTPC"}; Configurable cutNsigmaPrTPC{"cutNsigmaPrTPC", 5.0, "cutNsigmaPrTPC"}; + Configurable cutNsigmaKaTPC{"cutNsigmaKaTPC", 5.0, "cutNsigmaKaTPC"}; HistogramRegistry registry{"registry"}; void init(o2::framework::InitContext&) { if (fillhisto) { registry.add("hV0Candidate", "hV0Candidate", HistType::kTH1F, {{2, 0.5f, 2.5f}}); + registry.add("hCascCandidate", "hCascCandidate", HistType::kTH1F, {{2, 0.5f, 2.5f}}); registry.add("hMassGamma", "hMassGamma", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 0.0f, 0.1f}}); registry.add("hGammaRxy", "hGammaRxy", HistType::kTH2F, {{1800, -90.0f, 90.0f}, {1800, -90.0f, 90.0f}}); registry.add("hMassK0S", "hMassK0S", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 0.45, 0.55}}); @@ -196,10 +255,39 @@ struct v0selector { registry.add("hV0APplot", "hV0APplot", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}}); registry.add("hV0APplotSelected", "hV0APplotSelected", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}}); registry.add("hV0Psi", "hV0Psi", HistType::kTH2F, {{100, 0, TMath::PiOver2()}, {100, 0, 0.1}}); + if (selectCascades) { + registry.add("hCascPt", "pT", HistType::kTH1F, {{100, 0.0f, 10}}); + registry.add("hCascEtaPhi", "#eta vs. #varphi", HistType::kTH2F, {{63, 0, 6.3}, {20, -1.0f, 1.0f}}); + registry.add("hCascDCAxyPosToPV", "hCascDCAxyPosToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}}); + registry.add("hCascDCAxyNegToPV", "hCascDCAxyNegToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}}); + registry.add("hCascDCAxyBachToPV", "hCascDCAxyBachToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}}); + registry.add("hCascDCAzPosToPV", "hCascDCAzPosToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}}); + registry.add("hCascDCAzNegToPV", "hCascDCAzNegToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}}); + registry.add("hCascDCAzBachToPV", "hCascDCAzBachToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}}); + registry.add("hCascAPplot", "hCascAPplot", HistType::kTH2F, {{200, -1.0f, +1.0f}, {300, 0.0f, 0.3f}}); + registry.add("hCascV0APplot", "hCascV0APplot", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}}); + registry.add("hCascAPplotSelected", "hCascAPplotSelected", HistType::kTH2F, {{200, -1.0f, +1.0f}, {300, 0.0f, 0.3f}}); + registry.add("hCascV0APplotSelected", "hCascV0APplotSelected", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}}); + registry.add("hCascRadius", "hCascRadius", HistType::kTH1F, {{1000, 0.0f, 100.0f}}); + registry.add("hCascV0Radius", "hCascV0Radius", HistType::kTH1F, {{1000, 0.0f, 100.0f}}); + registry.add("hCascCosPA", "hCascCosPA", HistType::kTH1F, {{50, 0.95f, 1.0f}}); + registry.add("hCascV0CosPA", "hCascV0CosPA", HistType::kTH1F, {{50, 0.95f, 1.0f}}); + registry.add("hMassOmega", "hMassOmega", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 1.62f, 1.72f}}); + registry.add("hMassAntiOmega", "hMassAntiOmega", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 1.62f, 1.72f}}); + registry.add("hCascDCADau", "hCascDCADau", HistType::kTH1F, {{1000, 0.0f, 10.0f}}); + registry.add("hCascV0DCADau", "hCascV0DCADau", HistType::kTH1F, {{1000, 0.0f, 10.0f}}); + } + } + + if (selectCascades == false && produceCascID == true) { + LOGP(error, "produceCascID is available only when selectCascades is enabled"); + } + if (cutAPOmegaUp1 < cutAPOmegaDown1) { + LOGP(error, "cutAPOmegaUp1 must be greater than cutAPOmegaDown1"); } } - void process(aod::V0Datas const& V0s, FullTracksExt const& tracks, aod::Collisions const&) + void process(aod::V0Datas const& V0s, aod::CascDatas const& Cascs, FullTracksExt const& tracks, aod::Collisions const&) { std::vector pidmap; pidmap.clear(); @@ -210,6 +298,11 @@ struct v0selector { if (produceV0ID.value) { v0pidmap.resize(V0s.size(), -1); } + std::vector cascpidmap; + cascpidmap.clear(); + if (produceCascID.value) { + cascpidmap.resize(Cascs.size(), kUndef); + } for (auto& V0 : V0s) { // if (!(V0.posTrack_as().trackType() & o2::aod::track::TPCrefit)) { // continue; @@ -222,55 +315,22 @@ struct v0selector { if (fillhisto) { registry.fill(HIST("hV0Candidate"), 1); } - if (std::fabs(V0.posTrack_as().eta()) > 0.9) { - continue; - } - if (std::fabs(V0.negTrack_as().eta()) > 0.9) { - continue; - } - if (V0.posTrack_as().tpcNClsCrossedRows() < mincrossedrows) { - continue; - } - if (V0.negTrack_as().tpcNClsCrossedRows() < mincrossedrows) { - continue; - } - if (V0.posTrack_as().tpcChi2NCl() > maxchi2tpc) { - continue; - } - if (V0.negTrack_as().tpcChi2NCl() > maxchi2tpc) { - continue; - } - if (std::fabs(V0.posTrack_as().dcaXY()) < dcamin) { - continue; - } - if (std::fabs(V0.negTrack_as().dcaXY()) < dcamin) { - continue; - } - if (std::fabs(V0.posTrack_as().dcaXY()) > dcamax) { - continue; - } - if (std::fabs(V0.negTrack_as().dcaXY()) > dcamax) { - continue; - } + const auto& posTrack = V0.posTrack_as(); + const auto& negTrack = V0.negTrack_as(); - if (V0.posTrack_as().sign() * V0.negTrack_as().sign() > 0) { // reject same sign pair - continue; + bool isRejectV0{false}; + for (const auto& prong : {posTrack, negTrack}) { + isRejectV0 = isRejectV0 || std::fabs(prong.eta()) > 0.9; + isRejectV0 = isRejectV0 || prong.tpcNClsCrossedRows() < mincrossedrows; + isRejectV0 = isRejectV0 || prong.tpcChi2NCl() > maxchi2tpc; + isRejectV0 = isRejectV0 || std::fabs(prong.dcaXY()) < dcamin; + isRejectV0 = isRejectV0 || std::fabs(prong.dcaXY()) > dcamax; } + isRejectV0 = isRejectV0 || (posTrack.sign() * negTrack.sign() > 0); - // if (V0.posTrack_as().collisionId() != V0.negTrack_as().collisionId()) { - // continue; - // } - - // if (!V0.posTrack_as().has_collision() || !V0.negTrack_as().has_collision()) { - // continue; - // } - - // auto const& collision = V0.collision_as(); - - // if (V0.collisionId() != collision.globalIndex()) { - // continue; - // } + if (isRejectV0) + continue; float V0dca = V0.dcaV0daughters(); float V0CosinePA = V0.v0cosPA(); @@ -376,6 +436,134 @@ struct v0selector { v0mapID(v0pidmap[V0.globalIndex()]); } } + + if (selectCascades) { + for (const auto& casc : Cascs) { + if (fillhisto) { + registry.fill(HIST("hCascCandidate"), 1); + } + + const auto& posTrack = casc.posTrack_as(); + const auto& negTrack = casc.negTrack_as(); + const auto& bachelor = casc.bachelor_as(); + + bool isRejectCascade{false}; + for (const auto& prong : {posTrack, negTrack, bachelor}) { + isRejectCascade = isRejectCascade || std::fabs(prong.eta()) > 0.9; + isRejectCascade = isRejectCascade || prong.tpcNClsCrossedRows() < mincrossedrows; + isRejectCascade = isRejectCascade || prong.tpcChi2NCl() > maxchi2tpc; + isRejectCascade = isRejectCascade || std::fabs(prong.dcaXY()) < dcamin; + isRejectCascade = isRejectCascade || std::fabs(prong.dcaXY()) > dcamax; + } + isRejectCascade = isRejectCascade || (posTrack.sign() * negTrack.sign() > 0); + + if (isRejectCascade) + continue; + + if (fillhisto) { + registry.fill(HIST("hCascCandidate"), 2); + } + + auto collision = casc.collision_as(); + const float collisionX = collision.posX(); + const float collisionY = collision.posY(); + const float collisionZ = collision.posZ(); + + const float cascDca = casc.dcacascdaughters(); // NOTE the name of getter is misleading. In case of no-KF this is sqrt(Chi2) + const float cascV0Dca = casc.dcaV0daughters(); // NOTE the name of getter is misleading. In case of kfDoDCAFitterPreMinimV0 this is sqrt(Chi2) + const float cascRadius = casc.cascradius(); + const float cascV0Radius = casc.dcaV0daughters(); + const float cascCosinePA = casc.casccosPA(collisionX, collisionY, collisionZ); + const float cascV0CosinePA = casc.v0cosPA(collisionX, collisionY, collisionZ); + + if (cascDca > cascDcaMax) { + continue; + } + if (cascV0Dca > cascV0DcaMax) { + continue; + } + if (cascRadius < cascRadiusMin || cascRadius > cascRadiusMax) { + continue; + } + if (cascV0Radius < cascV0RadiusMin || cascV0Radius > cascV0RadiusMax) { + continue; + } + if (cascCosinePA < cascCosinePAMin) { + continue; + } + if (cascV0CosinePA < cascV0CosinePAMin || cascV0CosinePA > cascV0CosinePAMax) { + continue; + } + + const float mOmega = casc.mOmega(); + const float mV0Lambda = casc.mLambda(); + const float alpha = casc.alpha(); + const float qt = casc.qtarm(); + const float v0Alpha = casc.v0Alpha(); + const float v0Qt = casc.v0Qtarm(); + + if (fillhisto) { + registry.fill(HIST("hCascPt"), casc.pt()); + registry.fill(HIST("hCascEtaPhi"), casc.phi(), casc.eta()); + registry.fill(HIST("hCascDCAxyPosToPV"), casc.posTrack_as().dcaXY()); + registry.fill(HIST("hCascDCAxyNegToPV"), casc.negTrack_as().dcaXY()); + registry.fill(HIST("hCascDCAxyBachToPV"), casc.bachelor_as().dcaXY()); + registry.fill(HIST("hCascDCAzPosToPV"), casc.posTrack_as().dcaZ()); + registry.fill(HIST("hCascDCAzNegToPV"), casc.negTrack_as().dcaZ()); + registry.fill(HIST("hCascDCAzBachToPV"), casc.bachelor_as().dcaZ()); + registry.fill(HIST("hCascAPplot"), alpha, qt); + registry.fill(HIST("hCascV0APplot"), v0Alpha, v0Qt); + registry.fill(HIST("hCascRadius"), cascRadius); + registry.fill(HIST("hCascV0Radius"), cascV0Radius); + registry.fill(HIST("hCascCosPA"), cascCosinePA); + registry.fill(HIST("hCascV0CosPA"), cascV0CosinePA); + registry.fill(HIST("hCascDCADau"), cascDca); + registry.fill(HIST("hCascV0DCADau"), cascV0Dca); + } + + const int cascid = checkCascade(alpha, qt); + const int v0id = checkV0(v0Alpha, v0Qt); + if (cascid < 0) { + continue; + } + if (v0id != kLambda && v0id != kAntiLambda) { + continue; + } + if (fillhisto) { + registry.fill(HIST("hCascAPplotSelected"), alpha, qt); + registry.fill(HIST("hCascV0APplotSelected"), v0Alpha, v0Qt); + } + + auto storeCascAddID = [&](auto gix, auto id) { + if (produceCascID.value) { + cascpidmap[gix] = id; + } + }; + + if (cascid == kOmega && v0id == kLambda) { + if (fillhisto) { + registry.fill(HIST("hMassOmega"), cascRadius, mOmega); + } + if (cutMassOmegaLow < mOmega && mOmega < cutMassOmegaHigh && cutMassCascV0Low < mV0Lambda && mV0Lambda < cutMassCascV0High && std::abs(casc.posTrack_as().tpcNSigmaPr()) < cutNsigmaPrTPC && std::abs(casc.negTrack_as().tpcNSigmaPi()) < cutNsigmaPiTPC && std::abs(casc.bachelor_as().tpcNSigmaKa()) < cutNsigmaKaTPC) { + pidmap[casc.bachelorId()] |= (uint8_t(1) << kOmega); + storeCascAddID(casc.globalIndex(), kOmega); + } + } else if (cascid == kAntiOmega && v0id == kAntiLambda) { + if (fillhisto) { + registry.fill(HIST("hMassAntiOmega"), cascRadius, mOmega); + } + if (cutMassOmegaLow < mOmega && mOmega < cutMassOmegaHigh && cutMassCascV0Low < mV0Lambda && mV0Lambda < cutMassCascV0High && std::abs(casc.posTrack_as().tpcNSigmaPi()) < cutNsigmaPiTPC && std::abs(casc.negTrack_as().tpcNSigmaPr()) < cutNsigmaPrTPC && std::abs(casc.bachelor_as().tpcNSigmaKa()) < cutNsigmaKaTPC) { + pidmap[casc.bachelorId()] |= (uint8_t(1) << kAntiOmega); + storeCascAddID(casc.globalIndex(), kAntiOmega); + } + } + } // end of Casc loop + if (produceCascID.value) { + for (auto& casc : Cascs) { + cascmapID(cascpidmap[casc.globalIndex()]); + } + } + } for (auto& track : tracks) { // printf("setting pidmap[%lld] = %d\n",track.globalIndex(),pidmap[track.globalIndex()]); v0bits(pidmap[track.globalIndex()]); @@ -574,7 +762,7 @@ struct trackPIDQA { } } // end of track loop - } // end of process + } // end of process void DefineHistograms(TString histClasses) { diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h index ca58973dff8..2bb5c2f06f2 100644 --- a/PWGLF/DataModel/LFStrangenessTables.h +++ b/PWGLF/DataModel/LFStrangenessTables.h @@ -713,21 +713,21 @@ DECLARE_SOA_DYNAMIC_COLUMN(MLambda, mLambda, //! mass under lambda hypothesis [](int v0type, float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { if ((v0type & (0x1 << o2::dataformats::V0Index::kPhotonOnly)) != 0) { return 0.0f; // provide mass only if NOT a photon with TPC-only tracks (special handling) - }; + } return RecoDecay::m(std::array{std::array{pxpos, pypos, pzpos}, std::array{pxneg, pyneg, pzneg}}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); }); DECLARE_SOA_DYNAMIC_COLUMN(MAntiLambda, mAntiLambda, //! mass under antilambda hypothesis [](int v0type, float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { if ((v0type & (0x1 << o2::dataformats::V0Index::kPhotonOnly)) != 0) { return 0.0f; // provide mass only if NOT a photon with TPC-only tracks (special handling) - }; + } return RecoDecay::m(std::array{std::array{pxpos, pypos, pzpos}, std::array{pxneg, pyneg, pzneg}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton}); }); DECLARE_SOA_DYNAMIC_COLUMN(MK0Short, mK0Short, //! mass under K0short hypothesis [](int v0type, float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { if ((v0type & (0x1 << o2::dataformats::V0Index::kPhotonOnly)) != 0) { return 0.0f; // provide mass only if NOT a photon with TPC-only tracks (special handling) - }; + } return RecoDecay::m(std::array{std::array{pxpos, pypos, pzpos}, std::array{pxneg, pyneg, pzneg}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassPionCharged}); }); DECLARE_SOA_DYNAMIC_COLUMN(MLambda_unchecked, mLambda_unchecked, //! mass under lambda hypothesis without v0 type check (will include TPC only and potentially duplicates! use with care) @@ -1417,6 +1417,68 @@ DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, //! Cascade phi in the range [0, 2pi) [](float px, float py) -> float { return RecoDecay::phi(px, py); }); DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, //! Cascade pseudorapidity [](float px, float py, float pz) -> float { return RecoDecay::eta(std::array{px, py, pz}); }); + +// Armenteros-Podolanski variables +DECLARE_SOA_DYNAMIC_COLUMN(Alpha, alpha, //! Armenteros Alpha + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg, float pxbach, float pybach, float pzbach, int sign) { + const float pxv0 = pxpos + pxneg; + const float pyv0 = pypos + pyneg; + const float pzv0 = pzpos + pzneg; + + // No need to divide by momentum of the cascade (as in the v0data namespace) since the ratio of lQl is evaluated + const float lQlBach = RecoDecay::dotProd(std::array{pxbach, pybach, pzbach}, std::array{pxv0 + pxbach, pyv0 + pybach, pzv0 + pzbach}); + const float lQlV0 = RecoDecay::dotProd(std::array{pxv0, pyv0, pzv0}, std::array{pxv0 + pxbach, pyv0 + pybach, pzv0 + pzbach}); + float alpha = (lQlBach - lQlV0) / (lQlV0 + lQlBach); // alphacascade + if (sign < 0) { + alpha = -alpha; + } + return alpha; + }); + +DECLARE_SOA_DYNAMIC_COLUMN(QtArm, qtarm, //! Armenteros Qt + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg, float pxbach, float pybach, float pzbach) { + const float pxv0 = pxpos + pxneg; + const float pyv0 = pypos + pyneg; + const float pzv0 = pzpos + pzneg; + + const float momTot2 = RecoDecay::p2(pxv0 + pxbach, pyv0 + pybach, pzv0 + pzbach); + const float dp = RecoDecay::dotProd(std::array{pxbach, pybach, pzbach}, std::array{pxv0 + pxbach, pyv0 + pybach, pzv0 + pzbach}); + return std::sqrt(RecoDecay::p2(pxbach, pybach, pzbach) - dp * dp / momTot2); // qtarm + }); + +// Psi pair angle: angle between the plane defined by the v0 and bachelor momenta and the xy plane +DECLARE_SOA_DYNAMIC_COLUMN(PsiPair, psipair, //! psi pair angle + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg, float pxbach, float pybach, float pzbach, int sign) { + auto clipToPM1 = [](float x) { return x < -1.f ? -1.f : (x > 1.f ? 1.f : x); }; + const float pxv0 = pxpos + pxneg; + const float pyv0 = pypos + pyneg; + const float pzv0 = pzpos + pzneg; + + const float ptot2 = RecoDecay::p2(pxbach, pybach, pzbach) * RecoDecay::p2(pxv0, pyv0, pzv0); + const float argcos = RecoDecay::dotProd(std::array{pxbach, pybach, pzbach}, std::array{pxv0, pyv0, pzv0}) / std::sqrt(ptot2); + const float thetaV0 = std::atan2(RecoDecay::sqrtSumOfSquares(pxv0, pyv0), pzv0); + const float thetaBach = std::atan2(RecoDecay::sqrtSumOfSquares(pxbach, pybach), pzbach); + float argsin = (thetaV0 - thetaBach) / std::acos(clipToPM1(argcos)); + if (sign < 0) { + argsin = -argsin; + } + return std::asin(clipToPM1(argsin)); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(V0Alpha, v0Alpha, //! Armenteros Alpha + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { + // No need to divide by momentum of the v0 (as in the v0data namespace) since the ratio of lQl is evaluated + const float lQlPos = RecoDecay::dotProd(std::array{pxpos, pypos, pzpos}, std::array{pxneg + pxpos, pyneg + pypos, pzneg + pzpos}); + const float lQlNeg = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxneg + pxpos, pyneg + pypos, pzneg + pzpos}); + return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(V0QtArm, v0Qtarm, //! Armenteros Qt + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { + const float momTot2 = RecoDecay::p2(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); + const float dp = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}); + return std::sqrt(RecoDecay::p2(pxneg, pyneg, pzneg) - dp * dp / momTot2); // qtarm + }); } // namespace cascdata //______________________________________________________ @@ -1493,7 +1555,14 @@ DECLARE_SOA_TABLE(StoredCascCores, "AOD", "CASCCORE", //! core information about cascdata::PositiveEta, cascdata::PositivePhi, cascdata::BachelorEta, - cascdata::BachelorPhi); + cascdata::BachelorPhi, + + // Armenteros-Podolanski and psi-pair + cascdata::Alpha, + cascdata::QtArm, + cascdata::PsiPair, + cascdata::V0Alpha, + cascdata::V0QtArm); DECLARE_SOA_TABLE(StoredKFCascCores, "AOD", "KFCASCCORE", //! cascdata::Sign, cascdata::MXi, cascdata::MOmega,