diff --git a/PWGCF/FemtoDream/CMakeLists.txt b/PWGCF/FemtoDream/CMakeLists.txt index 6a7a4f19eee..1cfa1ddee34 100644 --- a/PWGCF/FemtoDream/CMakeLists.txt +++ b/PWGCF/FemtoDream/CMakeLists.txt @@ -19,6 +19,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-v0 + SOURCES femtoDreamPairTaskTrackV0.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(femtodream-hash SOURCES femtoDreamHashTask.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore diff --git a/PWGCF/FemtoDream/femtoDreamPairTaskTrackV0.cxx b/PWGCF/FemtoDream/femtoDreamPairTaskTrackV0.cxx new file mode 100644 index 00000000000..6b992509902 --- /dev/null +++ b/PWGCF/FemtoDream/femtoDreamPairTaskTrackV0.cxx @@ -0,0 +1,289 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file femtoDreamPairTaskTrackTrack.cxx +/// \brief Tasks that reads the track tables used for the pairing and builds pairs of two tracks +/// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de + +#include "PWGCF/DataModel/FemtoDerived.h" +#include "FemtoDreamParticleHisto.h" +#include "FemtoDreamEventHisto.h" +#include "FemtoDreamPairCleaner.h" +#include "FemtoDreamContainer.h" +#include "FemtoDreamDetaDphiStar.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/ASoAHelpers.h" + +using namespace o2; +using namespace o2::analysis::femtoDream; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace +{ +static constexpr int nCuts = 4; +static const std::vector cutNames{"MaxPt", "PIDthr", "nSigmaTPC", "nSigmaTPCTOF"}; +static const float cutsArray[nCuts] = {4.05f, 0.75f, 3.f, 3.f}; + +static constexpr int nNsigma = 3; +static constexpr float kNsigma[nNsigma] = {3.5f, 3.f, 2.5f}; + +enum kPIDselection { + k3d5sigma = 0, + k3sigma = 1, + k2d5sigma = 2 +}; + +enum kDetector { + kTPC = 0, + kTPCTOF = 1, + kNdetectors = 2 +}; +} // namespace + +struct femtoDreamPairTaskTrackV0 { + + /// Particle 1 (track) + Configurable ConfPDGCodePartOne{"ConfPDGCodePartOne", 2212, "Particle 1 - PDG code"}; + Configurable ConfCutPartOne{"ConfCutPartOne", 5542986, "Particle 1 - Selection bit from cutCulator"}; + Configurable> ConfPIDPartOne{"ConfPIDPartOne", std::vector{2}, "Particle 1 - Read from cutCulator"}; // we also need the possibility to specify whether the bit is true/false ->std>>vector>int>> + Configurable> cfgCutArray{"cfgCutArray", {cutsArray, nCuts, cutNames}, "Particle selections"}; + Configurable cfgNspecies{"ccfgNspecies", 4, "Number of particle spieces with PID info"}; + + /// Partition for particle 1 + Partition partsOne = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kTrack)) && + (aod::femtodreamparticle::pt < cfgCutArray->get("MaxPt")) && + ((aod::femtodreamparticle::cut & aod::femtodreamparticle::cutContainerType(ConfCutPartOne)) == aod::femtodreamparticle::cutContainerType(ConfCutPartOne)); + + /// Histogramming for particle 1 + FemtoDreamParticleHisto trackHistoPartOne; + + /// Particle 2 (V0) + Configurable ConfPDGCodePartTwo{"ConfPDGCodePartTwo", 3122, "Particle 1 - PDG code"}; + Configurable ConfCutPartTwo{"ConfCutPartTwo", 338, "Particle 2 - Selection bit"}; + + /// Partition for particle 2 + Partition partsTwo = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kV0)) && + ((aod::femtodreamparticle::cut & aod::femtodreamparticle::cutContainerType(ConfCutPartTwo)) == aod::femtodreamparticle::cutContainerType(ConfCutPartTwo)); + + /// Histogramming for particle 2 + FemtoDreamParticleHisto trackHistoPartTwo; + + /// Histogramming for Event + FemtoDreamEventHisto eventHisto; + + /// The configurables need to be passed to an std::vector + std::vector vPIDPartOne; + + /// Correlation part + ConfigurableAxis CfgMultBins{"CfgMultBins", {VARIABLE_WIDTH, 0.0f, 20.0f, 40.0f, 60.0f, 80.0f, 100.0f, 200.0f, 99999.f}, "Mixing bins - multiplicity"}; + ConfigurableAxis CfgkstarBins{"CfgkstarBins", {1500, 0., 6.}, "binning kstar"}; + ConfigurableAxis CfgkTBins{"CfgkTBins", {150, 0., 9.}, "binning kT"}; + ConfigurableAxis CfgmTBins{"CfgmTBins", {225, 0., 7.5}, "binning mT"}; + Configurable ConfNEventsMix{"ConfNEventsMix", 5, "Number of events for mixing"}; + Configurable ConfIsCPR{"ConfIsCPR", true, "Close Pair Rejection"}; + Configurable ConfBField{"ConfBField", +0.5, "Magnetic Field"}; + + FemtoDreamContainer sameEventCont; + FemtoDreamContainer mixedEventCont; + FemtoDreamPairCleaner pairCleaner; + FemtoDreamDetaDphiStar pairCloseRejection; + /// Histogram output + HistogramRegistry qaRegistry{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(InitContext&) + { + eventHisto.init(&qaRegistry); + trackHistoPartOne.init(&qaRegistry); + trackHistoPartTwo.init(&qaRegistry); + + sameEventCont.init(&resultRegistry, CfgkstarBins, CfgMultBins, CfgkTBins, CfgmTBins); + sameEventCont.setPDGCodes(ConfPDGCodePartOne, ConfPDGCodePartTwo); + mixedEventCont.init(&resultRegistry, CfgkstarBins, CfgMultBins, CfgkTBins, CfgmTBins); + mixedEventCont.setPDGCodes(ConfPDGCodePartOne, ConfPDGCodePartTwo); + pairCleaner.init(&qaRegistry); + if (ConfIsCPR) { + pairCloseRejection.init(&resultRegistry, &qaRegistry, 0.01, 0.01, ConfBField, false); /// \todo add config for Δη and ΔΦ cut values + } + + vPIDPartOne = ConfPIDPartOne; + } + + /// internal function that returns the kPIDselection element corresponding to a specifica n-sigma value + /// \param number of sigmas for PIF + /// \return kPIDselection corresponing to n-sigma + kPIDselection getPIDselection(float nSigma) + { + for (int i = 0; i < nNsigma; i++) { + if (abs(nSigma - kNsigma[i]) < 1e-3) { + return static_cast(i); + } + } + LOG(info) << "Invalid value of nSigma: " << nSigma << ". Standard 3 sigma returned." << std::endl; + return kPIDselection::k3sigma; + } + + /// function that checks whether the PID selection specified in the vectors is fulfilled + /// \param pidcut Bit-wise container for the PID + /// \param vSpecies Vector with the different selections + /// \return Whether the PID selection specified in the vectors is fulfilled + bool isPIDSelected(aod::femtodreamparticle::cutContainerType const& pidcut, std::vector const& vSpecies, float nSigma, kDetector iDet = kDetector::kTPC) + { + bool pidSelection = true; + kPIDselection kNsigma = getPIDselection(nSigma); + for (auto iSpecies : vSpecies) { + //\todo we also need the possibility to specify whether the bit is true/false ->std>>vector> + // if (!((pidcut >> it.first) & it.second)) { + int bit_to_check = cfgNspecies * kDetector::kNdetectors * kNsigma + iSpecies * kDetector::kNdetectors + iDet; + if (!(pidcut & (1UL << bit_to_check))) { + pidSelection = false; + } + } + return pidSelection; + }; + + /// function that checks whether the PID selection specified in the vectors is fulfilled, depending on the momentum TPC or TPC+TOF PID is conducted + /// \param pidcut Bit-wise container for the PID + /// \param mom Momentum of the track + /// \param pidThresh Momentum threshold that separates between TPC and TPC+TOF PID + /// \param vecTPC Vector with the different selections for the TPC PID + /// \param vecComb Vector with the different selections for the TPC+TOF PID + /// \return Whether the PID selection is fulfilled + bool isFullPIDSelected(aod::femtodreamparticle::cutContainerType const& pidCut, float const& momentum, float const& pidThresh, std::vector const& vSpecies, float nSigmaTPC = 3.5, float nSigmaTPCTOF = 3.5) + { + bool pidSelection = true; + if (momentum < pidThresh) { + /// TPC PID only + pidSelection = isPIDSelected(pidCut, vSpecies, nSigmaTPC, kDetector::kTPC); + } else { + /// TPC + TOF PID + pidSelection = isPIDSelected(pidCut, vSpecies, nSigmaTPCTOF, kDetector::kTPCTOF); + } + return pidSelection; + }; + + /// 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, ... + void processSameEvent(o2::aod::FemtoDreamCollision& col, + o2::aod::FemtoDreamParticles& parts) + { + const int multCol = col.multV0M(); + eventHisto.fillQA(col); + /// Histogramming same event + for (auto& part : partsOne) { + if (!isFullPIDSelected(part.pidcut(), part.p(), cfgCutArray->get("PIDthr"), vPIDPartOne, cfgCutArray->get("nSigmaTPC"), cfgCutArray->get("nSigmaTPCTOF"))) { + continue; + } + trackHistoPartOne.fillQA(part); + } + for (auto& part : partsTwo) { + trackHistoPartTwo.fillQA(part); + } + /// Now build the combinations + for (auto& [p1, p2] : combinations(partsOne, partsTwo)) { + if (!isFullPIDSelected(p1.pidcut(), p1.p(), cfgCutArray->get("PIDthr"), vPIDPartOne, cfgCutArray->get("nSigmaTPC"), cfgCutArray->get("nSigmaTPCTOF"))) { + continue; + } + + /// close pair rejection + if (ConfIsCPR) { + if (pairCloseRejection.isClosePair(p1, p2, parts)) { + continue; + } + } + + // track cleaning + if (!pairCleaner.isCleanPair(p1, p2, parts)) { + continue; + } + sameEventCont.setPair(p1, p2, multCol); + } + } + + PROCESS_SWITCH(femtoDreamPairTaskTrackV0, processSameEvent, "Enable processing same event", true); + + /// 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, ... + void processMixedEvent(o2::aod::FemtoDreamCollisions& cols, + o2::aod::Hashes& hashes, + o2::aod::FemtoDreamParticles& parts) + { + cols.bindExternalIndices(&parts); + auto particlesTuple = std::make_tuple(parts); + AnalysisDataProcessorBuilder::GroupSlicer slicer(cols, particlesTuple); + + for (auto& [collision1, collision2] : soa::selfCombinations("fBin", ConfNEventsMix, -1, soa::join(hashes, cols), soa::join(hashes, cols))) { + auto it1 = slicer.begin(); + auto it2 = slicer.begin(); + for (auto& slice : slicer) { + if (slice.groupingElement().index() == collision1.index()) { + it1 = slice; + break; + } + } + for (auto& slice : slicer) { + if (slice.groupingElement().index() == collision2.index()) { + it2 = slice; + break; + } + } + + auto particles1 = std::get(it1.associatedTables()); + particles1.bindExternalIndices(&cols); + auto particles2 = std::get(it2.associatedTables()); + particles2.bindExternalIndices(&cols); + + partsOne.bindTable(particles1); + partsTwo.bindTable(particles2); + + /// \todo before mixing we should check whether both collisions contain a pair of particles! + /// could work like that, but only if PID is contained within the partitioning! + // auto particlesEvent1 = std::get(it1.associatedTables()); + // particlesEvent1.bindExternalIndices(&cols); + // auto particlesEvent2 = std::get(it2.associatedTables()); + // particlesEvent2.bindExternalIndices(&cols); + /// for the x-check + // partsOne.bindTable(particlesEvent2); + // auto nPart1Evt2 = partsOne.size(); + // partsTwo.bindTable(particlesEvent1); + // auto nPart2Evt1 = partsTwo.size(); + /// for actual event mixing + // partsOne.bindTable(particlesEvent1); + // partsTwo.bindTable(particlesEvent2); + // if (partsOne.size() == 0 || nPart2Evt1 == 0 || nPart1Evt2 == 0 || partsTwo.size() == 0 ) continue; + + for (auto& [p1, p2] : combinations(partsOne, partsTwo)) { + if (!isFullPIDSelected(p1.pidcut(), p1.p(), cfgCutArray->get("PIDthr"), vPIDPartOne, cfgCutArray->get("nSigmaTPC"), cfgCutArray->get("nSigmaTPCTOF"))) { + continue; + } + if (ConfIsCPR) { + if (pairCloseRejection.isClosePair(p1, p2, parts)) { + continue; + } + } + mixedEventCont.setPair(p1, p2, collision1.multV0M()); // < \todo dirty trick, the multiplicity will be of course within the bin width used for the hashes + } + } + } + + PROCESS_SWITCH(femtoDreamPairTaskTrackV0, processMixedEvent, "Enable processing mixed events", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc), + }; + return workflow; +}