diff --git a/PWGCF/FemtoDream/CMakeLists.txt b/PWGCF/FemtoDream/CMakeLists.txt index 6cceb1d200e..383f229439d 100644 --- a/PWGCF/FemtoDream/CMakeLists.txt +++ b/PWGCF/FemtoDream/CMakeLists.txt @@ -24,6 +24,11 @@ o2physics_add_dpl_workflow(femtodream-pair-track-track PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(femtodream-pair-track-track-track + SOURCES femtoDreamPairTaskTrackTrackTrack.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(femtodream-debug-track SOURCES femtoDreamDebugTrack.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore diff --git a/PWGCF/FemtoDream/FemtoDreamContainerThreeBody.h b/PWGCF/FemtoDream/FemtoDreamContainerThreeBody.h new file mode 100644 index 00000000000..128b8a26169 --- /dev/null +++ b/PWGCF/FemtoDream/FemtoDreamContainerThreeBody.h @@ -0,0 +1,227 @@ +// Copyright 2019-2022 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 FemtoDreamContainer.h +/// \brief Definition of the FemtoDreamContainer +/// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de +/// \author Valentina Mantovani Sarti, valentina.mantovani-sarti@tum.de +/// \author Georgios Mantzaridis, TU München, georgios.mantzaridis@tum.de +/// \author Anton Riedel, TU München, anton.riedel@tum.de +/// \author Laura Serksnyte, TU München, laura.serksnyte@tum.de + +#ifndef PWGCF_FEMTODREAM_FEMTODREAMCONTAINERTHREEBODY_H_ +#define PWGCF_FEMTODREAM_FEMTODREAMCONTAINERTHREEBODY_H_ + +#include +#include +#include + +#include "Framework/HistogramRegistry.h" +#include "FemtoDreamMath.h" + +#include "Math/Vector4D.h" +#include "TMath.h" +#include "TDatabasePDG.h" + +using namespace o2::framework; + +namespace o2::analysis::femtoDream +{ + +namespace femtoDreamContainerThreeBody +{ +/// Femtoscopic observable to be computed +enum Observable { Q3 ///< Q3 +}; + +/// Type of the event processind +enum EventType { same, ///< Pair from same event + mixed ///< Pair from mixed event +}; +}; // namespace femtoDreamContainerThreeBody + +/// \class FemtoDreamContainerThreeBody +/// \brief Container for all histogramming related to the correlation function. The two +/// particles of the pair are passed here, and the correlation function and QA histograms +/// are filled according to the specified observable +/// \tparam eventType Type of the event (same/mixed) +/// \tparam obs Observable to be computed (Q_3/...) +template +class FemtoDreamContainerThreeBody +{ + public: + /// Destructor + virtual ~FemtoDreamContainerThreeBody() = default; + + /// Initializes histograms for the task in case of three-body femtoscopy + /// Called by init both in case of reconstructed data/ Monte Carlo, and for Monte Carlo Truth + /// \tparam T type of the axis Object + /// \param folderName Name of the directory in the output file (no suffix for reconstructed data/ Monte Carlo; "_MC" for Monte Carlo Truth) + /// \param femtoObs Title of the femto observable axis + /// \param femtoObsAxis axis object for the femto observable axis + /// \param multAxis axis object for the multiplicity axis + template + void init_base(std::string folderName, std::string femtoObs, T femtoObsAxis, T multAxis) + { + + mHistogramRegistry->add((folderName + "/relTripletDist").c_str(), ("; " + femtoObs + "; Entries").c_str(), kTH1F, {femtoObsAxis}); + mHistogramRegistry->add((folderName + "/relTripletQ3Mult").c_str(), ("; " + femtoObs + "; Multiplicity").c_str(), kTH2F, {femtoObsAxis, multAxis}); + } + + /// Initializes specialized Monte Carlo truth histograms for the task in case of three-body femtoscopy + /// internal function called by init only in case of Monte Carlo truth + /// \tparam T type of the xxis Object + /// \param folderName Name of the directory in the output file (no suffix for reconstructed data/ Monte Carlo; "_MC" for Monte Carlo Truth) + /// \param femtoObsAxis axis object for the femto observable axis + template + void init_MC(std::string folderName, std::string femtoObs, T femtoObsAxis, T multAxis) + { + mHistogramRegistry->add((folderName + "/relPairDist_ReconNoFake").c_str(), ("; " + femtoObs + "; Entries").c_str(), kTH1F, {femtoObsAxis}); + mHistogramRegistry->add((folderName + "/relPairQ3Mult_ReconNoFake").c_str(), ("; " + femtoObs + "; Multiplicity").c_str(), kTH2F, {femtoObsAxis, multAxis}); + mHistogramRegistry->add((folderName + "/hNoMCtruthTripletCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}}); + mHistogramRegistry->add((folderName + "/hFakeTripletCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}}); + mHistogramRegistry->add((folderName + "/Q3_resolution").c_str(), "; #it{Q}_{3} reconstructed (GeV/#it{c}); #it{Q}_{3} truth (GeV/#it{c})", kTH2F, {femtoObsAxis, femtoObsAxis}); + } + + /// Templated function to initialize the histograms for the task in case of three-body femtoscopy + /// Always calls init_base to initialize the histograms for data/ Monte Carlo reconstructed + /// In case of Monte Carlo, calls init_base again for Monte Carlo truth and the specialized function init_MC for additional histogramms + /// \tparam T type of the configurable for the axis configuration + /// \param registry Histogram registry to be passed + /// \param Q3Bins Q3 binning for the histograms + /// \param multBins multiplicity binning for the histograms + /// \param isMC add Monte Carlo truth histograms to the output file + template + void init(HistogramRegistry* registry, T& Q3Bins, T& multBins, bool isMC) + { + mHistogramRegistry = registry; + std::string femtoObs; + if constexpr (mFemtoObs == femtoDreamContainerThreeBody::Observable::Q3) { + femtoObs = "#it{Q}_{3} (GeV/#it{c})"; + } + std::vector tmpVecMult = multBins; + framework::AxisSpec multAxis = {tmpVecMult, "Multiplicity"}; + framework::AxisSpec femtoObsAxis = {Q3Bins, femtoObs.c_str()}; + + std::string folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]); + + init_base(folderName, femtoObs, femtoObsAxis, multAxis); + if (isMC) { + folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]); + init_base(folderName, femtoObs, femtoObsAxis, multAxis); + init_MC(folderName, femtoObs, femtoObsAxis, multAxis); + } + } + + /// Set the PDG codes of the three particles involved + /// \param pdg1 PDG code of particle one + /// \param pdg2 PDG code of particle two + /// \param pdg3 PDG code of particle three + void setPDGCodes(const int pdg1, const int pdg2, const int pdg3) + { + mMassOne = TDatabasePDG::Instance()->GetParticle(pdg1)->Mass(); + mMassTwo = TDatabasePDG::Instance()->GetParticle(pdg2)->Mass(); + mMassThree = TDatabasePDG::Instance()->GetParticle(pdg3)->Mass(); + mPDGOne = pdg1; + mPDGTwo = pdg2; + mPDGThree = pdg3; + } + + /// Pass a triplet to the container and compute all the relevant observables + /// Called by setTriplet both in case of data/ and Monte Carlo reconstructed and for Monte Carlo truth + /// \tparam T type of the femtodreamparticle + /// \param part1 Particle one + /// \param part2 Particle two + /// \param part3 Particle three + /// \param mult Multiplicity of the event + template + void setTriplet_base(const float femtoObs, T const& part1, T const& part2, T const& part3, const int mult) + { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/relTripletDist"), femtoObs); + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/relTripletQ3Mult"), femtoObs, mult); + } + + /// Called by setTriplet only in case of Monte Carlo truth + /// Fills MC truth specific histogramms: + /// - Q3 distribution plots with RECONSTRUCTED information but ONLY for non-fake candidates; needed for purity calculations of tracks + /// - Q3 resolution matrix + /// Note: Standard histograms with MC truth information are filled with the setPair_base function + /// \param part1 Particle one + /// \param part2 Particle two + /// \param part3 Particle three + /// \param mult Multiplicity of the event + void setTriplet_MC(const float femtoObsMC, const float femtoObs, const int mult) + { + if (mHistogramRegistry) { + // Fill the Q3 distributions with the reconstructed information but only for particles with the right PDG code + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/relPairDist_ReconNoFake"), femtoObs); + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/relPairkstarMult_ReconNoFake"), femtoObs, mult); + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/kstar_resolution"), femtoObsMC, femtoObs); + } + } + + /// Templated function to handle data/ Monte Carlo reconstructed and Monte Carlo truth + /// Always calls setTriplet_base to compute the observables with reconstructed data + /// In case of Monte Carlo, calls setTriplet_base with MC info and specialized function setTriplet_MC for additional histogramms + /// \tparam T type of the femtodreamparticle + /// \param part1 Particle one + /// \param part2 Particle two + /// \param part3 Particle three + /// \param mult Multiplicity of the event + template + void setTriplet(T const& part1, T const& part2, T const& part3, const int mult) + { + float femtoObs, femtoObsMC; + // Calculate femto observable and the mT with reconstructed information + if constexpr (mFemtoObs == femtoDreamContainerThreeBody::Observable::Q3) { + femtoObs = FemtoDreamMath::getQ3(part1, mMassOne, part2, mMassTwo, part3, mMassThree); + } + + if (mHistogramRegistry) { + setTriplet_base(femtoObs, part1, part2, part3, mult); + + if constexpr (isMC) { + if (part1.has_fdMCParticle() && part2.has_fdMCParticle() && part3.has_fdMCParticle()) { + // calculate the femto observable with MC truth information + if constexpr (mFemtoObs == femtoDreamContainerThreeBody::Observable::Q3) { + femtoObsMC = FemtoDreamMath::getQ3(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, part3.fdMCParticle(), mMassThree); + } + + if (abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo && abs(part3.fdMCParticle().pdgMCTruth()) == mPDGThree) { // Note: all triplet-histogramms are filled with MC truth information ONLY in case of non-fake candidates + setTriplet_base(femtoObsMC, part1.fdMCParticle(), part2.fdMCParticle(), part3.fdMCParticle(), mult); + setTriplet_MC(femtoObsMC, femtoObs, mult); + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakeTripletCounter"), 0); + } + + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hNoMCtruthTripletCounter"), 0); + } + } + } + } + + protected: + HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output + static constexpr std::string_view mFolderSuffix[2] = {"SameEvent", "MixedEvent"}; ///< Folder naming for the output according to mEventType + static constexpr femtoDreamContainerThreeBody::Observable mFemtoObs = obs; ///< Femtoscopic observable to be computed (according to femtoDreamContainerThreeBody::Observable) + static constexpr int mEventType = eventType; ///< Type of the event (same/mixed, according to femtoDreamContainerThreeBody::EventType) + float mMassOne = 0.f; ///< PDG mass of particle 1 + float mMassTwo = 0.f; ///< PDG mass of particle 1 + float mMassThree = 0.f; ///< PDG mass of particle 3; if relevant + int mPDGOne = 0; ///< PDG code of particle 1 + int mPDGTwo = 0; ///< PDG code of particle 2 + int mPDGThree = 0; ///< PDG code of particle 3; if relevant +}; + +} // namespace o2::analysis::femtoDream + +#endif // PWGCF_FEMTODREAM_FEMTODREAMCONTAINERTHREEBODY_H_ diff --git a/PWGCF/FemtoDream/FemtoDreamMath.h b/PWGCF/FemtoDream/FemtoDreamMath.h index 74d4c14aaac..259c59ceb87 100644 --- a/PWGCF/FemtoDream/FemtoDreamMath.h +++ b/PWGCF/FemtoDream/FemtoDreamMath.h @@ -92,13 +92,14 @@ class FemtoDreamMath template static float getQ3(const T& part1, const float mass1, const T& part2, const float mass2, const T& part3, const float mass3) { - float E1 = sqrt(pow(part1.px(), 2) + pow(part1.py(), 2) + pow(part1.pz(), 2) + pow(mass1, 2)); - float E2 = sqrt(pow(part2.px(), 2) + pow(part2.py(), 2) + pow(part2.pz(), 2) + pow(mass2, 2)); - float E3 = sqrt(pow(part3.px(), 2) + pow(part3.py(), 2) + pow(part3.pz(), 2) + pow(mass3, 2)); - const ROOT::Math::PxPyPzEVector vecpart1(part1.px(), part1.py(), part1.pz(), E1); - const ROOT::Math::PxPyPzEVector vecpart2(part2.px(), part2.py(), part2.pz(), E2); - const ROOT::Math::PxPyPzEVector vecpart3(part3.px(), part3.py(), part3.pz(), E3); + const ROOT::Math::PtEtaPhiMVector vecpart01(part1.pt(), part1.eta(), part1.phi(), mass1); + const ROOT::Math::PtEtaPhiMVector vecpart02(part2.pt(), part2.eta(), part2.phi(), mass2); + const ROOT::Math::PtEtaPhiMVector vecpart03(part3.pt(), part3.eta(), part3.phi(), mass3); + + const ROOT::Math::PxPyPzEVector vecpart1(vecpart01); + const ROOT::Math::PxPyPzEVector vecpart2(vecpart02); + const ROOT::Math::PxPyPzEVector vecpart3(vecpart03); ROOT::Math::PxPyPzEVector q12 = getqij(vecpart1, vecpart2); ROOT::Math::PxPyPzEVector q23 = getqij(vecpart2, vecpart3); diff --git a/PWGCF/FemtoDream/femtoDreamPairTaskTrackTrackTrack.cxx b/PWGCF/FemtoDream/femtoDreamPairTaskTrackTrackTrack.cxx new file mode 100644 index 00000000000..589cafc9e38 --- /dev/null +++ b/PWGCF/FemtoDream/femtoDreamPairTaskTrackTrackTrack.cxx @@ -0,0 +1,421 @@ +// Copyright 2019-2022 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 femtoDreamPairTaskTrackTrackTrack.cxx +/// \brief Tasks that reads the track tables and creates track triplets; should be currently used only for identical tracks +/// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de +/// \author Georgios Mantzaridis, TU München, georgios.mantzaridis@tum.de +/// \author Anton Riedel, TU München, anton.riedel@tum.de +/// \author Laura Serksnyte, TU München, laura.serksnyte@tum.de + +#include +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/StepTHn.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "TDatabasePDG.h" + +#include "PWGCF/DataModel/FemtoDerived.h" +#include "FemtoDreamParticleHisto.h" +#include "FemtoDreamEventHisto.h" +#include "FemtoDreamPairCleaner.h" +#include "FemtoDreamContainerThreeBody.h" +#include "FemtoDreamDetaDphiStar.h" +#include "FemtoUtils.h" + +using namespace o2; +using namespace o2::analysis::femtoDream; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; + +namespace +{ +static constexpr int nPart = 3; +static constexpr int nCuts = 5; +static const std::vector partNames{"PartOne", "PartTwo", "PartThree"}; +static const std::vector cutNames{"MaxPt", "PIDthr", "nSigmaTPC", "nSigmaTPCTOF", "MaxP"}; +static const float cutsTable[nPart][nCuts]{ + {4.05f, 1.f, 3.f, 3.f, 100.f}, + {4.05f, 1.f, 3.f, 3.f, 100.f}, + {4.05f, 1.f, 3.f, 3.f, 100.f}}; +} // namespace + +struct femtoDreamPairTaskTrackTrackTrack { + SliceCache cache; + Preslice perCol = aod::femtodreamparticle::fdCollisionId; + + /// Particle selection part + + /// Table for both particles + Configurable> ConfCutTable{"ConfCutTable", {cutsTable[0], nPart, nCuts, partNames, cutNames}, "Particle selections"}; + Configurable ConfNspecies{"ConfNspecies", 2, "Number of particle spieces with PID info"}; + Configurable ConfIsMC{"ConfIsMC", false, "Enable additional Histogramms in the case of a MonteCarlo Run"}; + Configurable> ConfTrkPIDnSigmaMax{"ConfTrkPIDnSigmaMax", std::vector{4.f, 3.f, 2.f}, "This configurable needs to be the same as the one used in the producer task"}; + Configurable ConfUse3D{"ConfUse3D", false, "Enable three dimensional histogramms (to be used only for analysis with high statistics): k* vs mT vs multiplicity"}; + + // Which particles to analyse; currently support only for same species and cuts triplets + Configurable ConfPDGCodePart{"ConfPDGCodePart", 2212, "Particle PDG code"}; + Configurable ConfPIDPart{"ConfPIDPart", 2, "Particles - Read from cutCulator"}; + Configurable ConfCutPart{"ConfCutPart", 5542474, "Particles - Selection bit from cutCulator"}; + + // Keep this format as next code update will include also different particle species and cuts: + + /// Particle 1 + + /// Partition for particle 1 + Partition partsOne = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kTrack)) && ((aod::femtodreamparticle::cut & ConfCutPart) == ConfCutPart); + Partition> partsOneMC = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kTrack)) && ((aod::femtodreamparticle::cut & ConfCutPart) == ConfCutPart); + + /// Histogramming for particle 1 + FemtoDreamParticleHisto trackHistoPartOne; + + /// Particle 2 + + /// Partition for particle 2 + Partition partsTwo = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kTrack)) && + // (aod::femtodreamparticle::pt < cfgCutTable->get("PartTwo", "MaxPt")) && + ((aod::femtodreamparticle::cut & ConfCutPart) == ConfCutPart); + Partition> partsTwoMC = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kTrack)) && + ((aod::femtodreamparticle::cut & ConfCutPart) == ConfCutPart); + + /// Particle 3 + + /// Partition for particle 3 + Partition partsThree = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kTrack)) && + // (aod::femtodreamparticle::pt < cfgCutTable->get("PartTwo", "MaxPt")) && + ((aod::femtodreamparticle::cut & ConfCutPart) == ConfCutPart); + Partition> partsThreeMC = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kTrack)) && + ((aod::femtodreamparticle::cut & ConfCutPart) == ConfCutPart); + + /// Histogramming for Event + FemtoDreamEventHisto eventHisto; + + /// The configurables need to be passed to an std::vector + int vPIDPart; + std::vector kNsigma; + + /// particle part + ConfigurableAxis ConfTempFitVarBins{"ConfDTempFitVarBins", {300, -0.15, 0.15}, "binning of the TempFitVar in the pT vs. TempFitVar plot"}; + ConfigurableAxis ConfTempFitVarpTBins{"ConfTempFitVarpTBins", {20, 0.5, 4.05}, "pT binning of the pT vs. TempFitVar plot"}; + + /// Correlation part + // ConfigurableAxis ConfMultBins{"ConfMultBins", {VARIABLE_WIDTH, 0.0f, 4.0f, 8.0f, 12.0f, 16.0f, 20.0f, 24.0f, 28.0f, 32.0f, 36.0f, 40.0f, 44.0f, 48.0f, 52.0f, 56.0f, 60.0f, 64.0f, 68.0f, 72.0f, 76.0f, 80.0f, 84.0f, 88.0f, 92.0f, 96.0f, 100.0f, 200.0f, 99999.f}, "Mixing bins - multiplicity"}; // \todo to be obtained from the hash task + ConfigurableAxis ConfMultBins{"CfgMultBins", {VARIABLE_WIDTH, 0.0f, 20.0f, 40.0f, 60.0f, 80.0f, 100.0f, 200.0f, 99999.f}, "Mixing bins - multiplicity"}; + ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + + ColumnBinningPolicy colBinning{{ConfVtxBins, ConfMultBins}, true}; + + ConfigurableAxis ConfQ3Bins{"ConfQ3Bins", {2000, 0., 8.}, "binning Q3"}; + Configurable ConfNEventsMix{"ConfNEventsMix", 5, "Number of events for mixing"}; + Configurable ConfIsCPR{"ConfIsCPR", true, "Close Pair Rejection"}; + Configurable ConfCPRPlotPerRadii{"ConfCPRPlotPerRadii", false, "Plot CPR per radii"}; + Configurable ConfCPRdeltaPhiMax{"ConfCPRdeltaPhiMax", 0.01, "Max. Delta Phi for Close Pair Rejection"}; + Configurable ConfCPRdeltaEtaMax{"ConfCPRdeltaEtaMax", 0.01, "Max. Delta Eta for Close Pair Rejection"}; + ConfigurableAxis ConfDummy{"ConfDummy", {1, 0, 1}, "Dummy axis"}; + + FemtoDreamContainerThreeBody sameEventCont; + FemtoDreamContainerThreeBody mixedEventCont; + FemtoDreamPairCleaner pairCleaner; + FemtoDreamDetaDphiStar pairCloseRejection; + /// Histogram output + HistogramRegistry qaRegistry{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry MixQaRegistry{"MixQaRegistry", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(InitContext&) + { + + eventHisto.init(&qaRegistry); + trackHistoPartOne.init(&qaRegistry, ConfTempFitVarpTBins, ConfTempFitVarBins, ConfDummy, ConfDummy, ConfDummy, ConfDummy, ConfIsMC, ConfPDGCodePart); + + MixQaRegistry.add("MixingQA/hSECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); + MixQaRegistry.add("MixingQA/hMECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); + + sameEventCont.init(&resultRegistry, ConfQ3Bins, ConfMultBins, ConfIsMC); + mixedEventCont.init(&resultRegistry, ConfQ3Bins, ConfMultBins, ConfIsMC); + sameEventCont.setPDGCodes(ConfPDGCodePart, ConfPDGCodePart, ConfPDGCodePart); + mixedEventCont.setPDGCodes(ConfPDGCodePart, ConfPDGCodePart, ConfPDGCodePart); + pairCleaner.init(&qaRegistry); // SERKSNYTE : later check if init should be updated to have 3 separate histos + if (ConfIsCPR.value) { + pairCloseRejection.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value); + } + + // CURRENTLY do only for one species + vPIDPart = ConfPIDPart.value; + kNsigma = ConfTrkPIDnSigmaMax.value; + } + + template + void fillCollision(CollisionType col) + { + MixQaRegistry.fill(HIST("MixingQA/hSECollisionBins"), colBinning.getBin({col.posZ(), col.multNtr()})); + eventHisto.fillQA(col); + } + + /// This function processes the same event and takes care of all the histogramming + /// \todo the trivial loops over the tracks should be factored out since they will be common to all combinations of T-T, T-V0, V0-V0, ... + /// @tparam PartitionType + /// @tparam PartType + /// @tparam isMC: enables Monte Carlo truth specific histograms + /// @param groupPartsOne partition for the first particle passed by the process function + /// @param groupPartsTwo partition for the second particle passed by the process function + /// @param parts femtoDreamParticles table (in case of Monte Carlo joined with FemtoDreamMCLabels) + /// @param magFieldTesla magnetic field of the collision + /// @param multCol multiplicity of the collision + template + void doSameEvent(PartitionType groupPartsOne, PartitionType groupPartsTwo, PartitionType groupPartsThree, PartType parts, float magFieldTesla, int multCol) + { + /// Histogramming same event + int NGoodTracks = 0; + for (auto& part : groupPartsOne) { + if (part.p() > ConfCutTable->get("PartOne", "MaxP") || part.pt() > ConfCutTable->get("PartOne", "MaxPt")) { + continue; + } + if (!isFullPIDSelected(part.pidcut(), + part.p(), + ConfCutTable->get("PartOne", "PIDthr"), + vPIDPart, + ConfNspecies, + kNsigma, + ConfCutTable->get("PartOne", "nSigmaTPC"), + ConfCutTable->get("PartOne", "nSigmaTPCTOF"))) { + continue; + } + NGoodTracks++; + + trackHistoPartOne.fillQA(part, aod::femtodreamparticle::kPt); + } + + /// Now build the combinations + for (auto& [p1, p2, p3] : combinations(CombinationsStrictlyUpperIndexPolicy(groupPartsOne, groupPartsTwo, groupPartsThree))) { + if (p1.p() > ConfCutTable->get("PartOne", "MaxP") || p1.pt() > ConfCutTable->get("PartOne", "MaxPt") || p2.p() > ConfCutTable->get("PartTwo", "MaxP") || p2.pt() > ConfCutTable->get("PartTwo", "MaxPt") || p3.p() > ConfCutTable->get("PartThree", "MaxP") || p3.pt() > ConfCutTable->get("PartThree", "MaxPt")) { + continue; + } + if (!isFullPIDSelected(p1.pidcut(), + p1.p(), + ConfCutTable->get("PartOne", "PIDthr"), + vPIDPart, + ConfNspecies, + kNsigma, + ConfCutTable->get("PartOne", "nSigmaTPC"), + ConfCutTable->get("PartOne", "nSigmaTPCTOF")) || + !isFullPIDSelected(p2.pidcut(), + p2.p(), + ConfCutTable->get("PartTwo", "PIDthr"), + vPIDPart, + ConfNspecies, + kNsigma, + ConfCutTable->get("PartTwo", "nSigmaTPC"), + ConfCutTable->get("PartTwo", "nSigmaTPCTOF")) || + !isFullPIDSelected(p3.pidcut(), + p3.p(), + ConfCutTable->get("PartThree", "PIDthr"), + vPIDPart, + ConfNspecies, + kNsigma, + ConfCutTable->get("PartThree", "nSigmaTPC"), + ConfCutTable->get("PartThree", "nSigmaTPCTOF"))) { + continue; + } + + if (ConfIsCPR.value) { + if (pairCloseRejection.isClosePair(p1, p2, parts, magFieldTesla)) { + continue; + } + if (pairCloseRejection.isClosePair(p2, p3, parts, magFieldTesla)) { + continue; + } + if (pairCloseRejection.isClosePair(p1, p3, parts, magFieldTesla)) { + continue; + } + } + + // track cleaning + if (!pairCleaner.isCleanPair(p1, p2, parts)) { + continue; + } + if (!pairCleaner.isCleanPair(p2, p3, parts)) { + continue; + } + if (!pairCleaner.isCleanPair(p1, p3, parts)) { + continue; + } + sameEventCont.setTriplet(p1, p2, p3, multCol); + } + } + + /// process function to call doSameEvent with Data + /// \param col subscribe to the collision table (Data) + /// \param parts subscribe to the femtoDreamParticleTable + void processSameEvent(o2::aod::FDCollision& col, + o2::aod::FDParticles& parts) + { + fillCollision(col); + + auto thegroupPartsOne = partsOne->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + auto thegroupPartsTwo = partsTwo->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + auto thegroupPartsThree = partsThree->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + + doSameEvent(thegroupPartsOne, thegroupPartsTwo, thegroupPartsTwo, parts, col.magField(), col.multNtr()); + } + PROCESS_SWITCH(femtoDreamPairTaskTrackTrackTrack, processSameEvent, "Enable processing same event", true); + + /// process function for to call doSameEvent with Monte Carlo + /// \param col subscribe to the collision table (Monte Carlo Reconstructed reconstructed) + /// \param parts subscribe to joined table FemtoDreamParticles and FemtoDreamMCLables to access Monte Carlo truth + /// \param FemtoDreamMCParticles subscribe to the Monte Carlo truth table + void processSameEventMC(o2::aod::FDCollision& col, + soa::Join& parts, + o2::aod::FDMCParticles&) + { + fillCollision(col); + + auto thegroupPartsOne = partsOneMC->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + auto thegroupPartsTwo = partsTwoMC->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + auto thegroupPartsThree = partsThreeMC->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + + doSameEvent(thegroupPartsOne, thegroupPartsTwo, thegroupPartsThree, parts, col.magField(), col.multNtr()); + } + PROCESS_SWITCH(femtoDreamPairTaskTrackTrackTrack, processSameEventMC, "Enable processing same event for Monte Carlo", false); + + /// This function processes the mixed event + /// \todo the trivial loops over the collisions and tracks should be factored out since they will be common to all combinations of T-T, T-V0, V0-V0, ... + /// \tparam PartitionType + /// \tparam PartType + /// \tparam isMC: enables Monte Carlo truth specific histograms + /// \param groupPartsOne partition for the first particle passed by the process function + /// \param groupPartsTwo partition for the second particle passed by the process function + /// \param groupPartsThree partition for the third particle passed by the process function + /// \param parts femtoDreamParticles table (in case of Monte Carlo joined with FemtoDreamMCLabels) + /// \param magFieldTesla magnetic field of the collision + /// \param multCol multiplicity of the collision + template + void doMixedEvent(PartitionType groupPartsOne, PartitionType groupPartsTwo, PartitionType groupPartsThree, PartType parts, float magFieldTesla, int multCol) + { + for (auto& [p1, p2, p3] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo, groupPartsThree))) { + if (p1.p() > ConfCutTable->get("PartOne", "MaxP") || p1.pt() > ConfCutTable->get("PartOne", "MaxPt") || p2.p() > ConfCutTable->get("PartTwo", "MaxP") || p2.pt() > ConfCutTable->get("PartTwo", "MaxPt") || p3.p() > ConfCutTable->get("PartThree", "MaxP") || p3.pt() > ConfCutTable->get("PartThree", "MaxPt")) { + continue; + } + + if (!isFullPIDSelected(p1.pidcut(), + p1.p(), + ConfCutTable->get("PartOne", "PIDthr"), + vPIDPart, + ConfNspecies, + kNsigma, + ConfCutTable->get("PartOne", "nSigmaTPC"), + ConfCutTable->get("PartOne", "nSigmaTPCTOF")) || + !isFullPIDSelected(p2.pidcut(), + p2.p(), + ConfCutTable->get("PartTwo", "PIDthr"), + vPIDPart, + ConfNspecies, + kNsigma, + ConfCutTable->get("PartTwo", "nSigmaTPC"), + ConfCutTable->get("PartTwo", "nSigmaTPCTOF")) || + !isFullPIDSelected(p3.pidcut(), + p3.p(), + ConfCutTable->get("PartThree", "PIDthr"), + vPIDPart, + ConfNspecies, + kNsigma, + ConfCutTable->get("PartThree", "nSigmaTPC"), + ConfCutTable->get("PartThree", "nSigmaTPCTOF"))) { + continue; + } + + if (ConfIsCPR.value) { + if (pairCloseRejection.isClosePair(p1, p2, parts, magFieldTesla)) { + continue; + } + if (pairCloseRejection.isClosePair(p2, p3, parts, magFieldTesla)) { + continue; + } + + if (pairCloseRejection.isClosePair(p1, p3, parts, magFieldTesla)) { + continue; + } + } + mixedEventCont.setTriplet(p1, p2, p3, multCol); + } + } + + /// process function for to call doMixedEvent with Data + /// @param cols subscribe to the collisions table (Data) + /// @param parts subscribe to the femtoDreamParticleTable + void processMixedEvent(o2::aod::FDCollisions& cols, + o2::aod::FDParticles& parts) + { + for (auto& [collision1, collision2, collision3] : soa::selfCombinations(colBinning, 5, -1, cols, cols, cols)) { + const int multiplicityCol = collision1.multNtr(); + MixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); + + auto groupPartsOne = partsOne->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision1.globalIndex(), cache); + auto groupPartsTwo = partsTwo->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision2.globalIndex(), cache); + auto groupPartsThree = partsThree->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision3.globalIndex(), cache); + + const auto& magFieldTesla1 = collision1.magField(); + const auto& magFieldTesla2 = collision2.magField(); + const auto& magFieldTesla3 = collision3.magField(); + + if ((magFieldTesla1 != magFieldTesla2) || (magFieldTesla2 != magFieldTesla3) || (magFieldTesla1 != magFieldTesla3)) { + continue; + } + // CONSIDER testing different strategies to which events to use + + doMixedEvent(groupPartsOne, groupPartsTwo, groupPartsThree, parts, magFieldTesla1, multiplicityCol); + } + } + PROCESS_SWITCH(femtoDreamPairTaskTrackTrackTrack, processMixedEvent, "Enable processing mixed events", true); + + /// brief process function for to call doMixedEvent with Monte Carlo + /// @param cols subscribe to the collisions table (Monte Carlo Reconstructed reconstructed) + /// @param parts subscribe to joined table FemtoDreamParticles and FemtoDreamMCLables to access Monte Carlo truth + /// @param FemtoDreamMCParticles subscribe to the Monte Carlo truth table + void processMixedEventMC(o2::aod::FDCollisions& cols, + soa::Join& parts, + o2::aod::FDMCParticles&) + { + for (auto& [collision1, collision2, collision3] : soa::selfCombinations(colBinning, 5, -1, cols, cols, cols)) { + + const int multiplicityCol = collision1.multNtr(); + MixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); + + auto groupPartsOne = partsOneMC->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision1.globalIndex(), cache); + auto groupPartsTwo = partsTwoMC->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision2.globalIndex(), cache); + auto groupPartsThree = partsThreeMC->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision3.globalIndex(), cache); + + const auto& magFieldTesla1 = collision1.magField(); + const auto& magFieldTesla2 = collision2.magField(); + const auto& magFieldTesla3 = collision3.magField(); + + if ((magFieldTesla1 != magFieldTesla2) || (magFieldTesla2 != magFieldTesla3) || (magFieldTesla1 != magFieldTesla3)) { + continue; + } + // CONSIDER testing different strategies to which events to use + + doMixedEvent(groupPartsOne, groupPartsTwo, groupPartsThree, parts, magFieldTesla1, multiplicityCol); + } + } + PROCESS_SWITCH(femtoDreamPairTaskTrackTrackTrack, processMixedEventMC, "Enable processing mixed events MC", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc), + }; + return workflow; +}