diff --git a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx index 5d16db9af39..14994cd0c16 100644 --- a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx +++ b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx @@ -43,9 +43,11 @@ struct PhianalysisTHnSparse { struct : ConfigurableGroup { Configurable produceQA{"produceQA", false, "Produce qa histograms."}; + Configurable produceStats{"produceStats", false, "Produce statistics histograms."}; Configurable produceTrue{"produceTrue", false, "Produce True and Gen histograms."}; Configurable produceLikesign{"produceLikesign", false, "Produce Like sign histograms."}; Configurable eventMixing{"eventMixing", "none", "Produce Event Mixing histograms of type."}; + Configurable produceRotational{"produceRotational", false, "Produce Rotational histograms."}; } produce; Configurable daughterPos{"daughterPos", 3, "Particle type of the positive dauther according to ReconstructionDataFormats/PID.h (Default = Kaon)"}; @@ -62,7 +64,6 @@ struct PhianalysisTHnSparse { Configurable ptTOFThreshold{"ptTOFThreshold", 0.5f, "Threshold for applying TOF."}; Configurable rapidity{"rapidity", 0.5f, "Rapidity cut (maximum)."}; Configurable etatrack{"etatrack", 0.8f, "Eta cut for track."}; - Configurable etapair{"etapair", 0.5f, "Eta cut for pair."}; Configurable pt{"pt", 0.15f, "Cut: Minimal value of tracks pt."}; Configurable dcaXY{"dcaXY", 1.0f, "Cut: Maximal value of tracks DCA XY."}; Configurable dcaZ{"dcaZ", 1.0f, "Cut: Maximal value of tracks DCA Z."}; @@ -97,6 +98,9 @@ struct PhianalysisTHnSparse { ConfigurableAxis axisMultiplicityMixing{"axisMultiplicityMixing", {5, 0, 5000}, "TPC multiplicity for bin"}; ConfigurableAxis axisCentralityMixing{"axisCentralityMixing", {10, 0, 100}, "Multiplicity percentil binning for mixing"}; + // rotational + Configurable numberofRotations{"numberofRotations", 1, "Number of rotations for rotational background estimation."}; + // Axes specifications AxisSpec posZaxis = {400, -20., 20., "V_{z} (cm)"}; AxisSpec dcaXYaxis = {800, -2.0, 2.0, "DCA_{xy} (cm)"}; @@ -115,7 +119,7 @@ struct PhianalysisTHnSparse { double* pointPair = nullptr; double* pointSys = nullptr; ROOT::Math::PxPyPzMVector d1, d2, mother; - bool produceQA, dataQA, MCTruthQA, globalTrack, tpcPidOnly = false; + bool produceTrue, produceLikesign, produceQA, produceStats, produceRotational, dataQA, MCTruthQA, globalTrack, tpcPidOnly = false; float tpcnSigmaPos = 100.0f; float tpcnSigmaNeg = 100.0f; float combinedNSigma = 100.0f; @@ -178,7 +182,11 @@ struct PhianalysisTHnSparse { std::vector allAxesSys = {tpcNClsFoundAxis}; produceQA = static_cast(produce.produceQA); + produceStats = static_cast(produce.produceStats); + produceTrue = static_cast(produce.produceTrue); + produceLikesign = static_cast(produce.produceLikesign); mixingType = rsn::mixingTypeName(static_cast(produce.eventMixing)); + produceRotational = static_cast(produce.produceRotational); tpcnSigmaPos = static_cast(cut.tpcnSigmaPos); tpcnSigmaNeg = static_cast(cut.tpcnSigmaNeg); tpcNClsFound = static_cast(cut.tpcNClsFound); @@ -190,11 +198,12 @@ struct PhianalysisTHnSparse { pointPair = new double[static_cast(o2::analysis::rsn::PairAxisType::unknown)]; pointSys = new double[static_cast(o2::analysis::rsn::SystematicsAxisType::unknown)]; rsnOutput = new o2::analysis::rsn::OutputSparse(); - rsnOutput->init(sparseAxes, allAxes, sysAxes, allAxesSys, static_cast(produce.produceTrue), mixingType, static_cast(produce.produceLikesign), ®istry); + rsnOutput->init(sparseAxes, allAxes, sysAxes, allAxesSys, produceTrue, mixingType, produceLikesign, produceRotational, ®istry); // Print summary of configuration LOGF(info, "=== PhianalysisTHnSparse configuration summary ==="); LOGF(info, "produceQA: %s", produceQA ? "true" : "false"); + LOGF(info, "produceStats: %s", produceStats ? "true" : "false"); LOGF(info, "produceTrue: %s", static_cast(produce.produceTrue) ? "true" : "false"); LOGF(info, "produceLikesign: %s", static_cast(produce.produceLikesign) ? "true" : "false"); LOGF(info, "eventMixing: %s", static_cast(produce.eventMixing).c_str()); @@ -208,7 +217,6 @@ struct PhianalysisTHnSparse { LOGF(info, "ptTOFThreshold: %.2f", ptTOFThreshold); LOGF(info, "rapidity: %.2f", static_cast(cut.rapidity)); LOGF(info, "etatrack: %.2f", static_cast(cut.etatrack)); - LOGF(info, "etapair: %.2f", static_cast(cut.etapair)); LOGF(info, "pt (min): %.2f", static_cast(cut.pt)); LOGF(info, "dcaXY: %.2f", static_cast(cut.dcaXY)); LOGF(info, "dcaZ: %.2f", static_cast(cut.dcaZ)); @@ -216,6 +224,7 @@ struct PhianalysisTHnSparse { LOGF(info, "tpcNClsFound: %d", tpcNClsFound); LOGF(info, "mixingType: %d", static_cast(mixingType)); LOGF(info, "numberofMixedEvents: %d", static_cast(numberofMixedEvents)); + LOGF(info, "numberofRotations: %d", static_cast(numberofRotations)); LOGF(info, "sparseAxes: "); for (const auto& axis : static_cast>(sparseAxes)) { LOGF(info, " %s", axis.c_str()); @@ -240,18 +249,17 @@ struct PhianalysisTHnSparse { registry.add("QAEvent/hMult", "Multiplicity (amplitude of non-zero channels in the FV0A + FV0C) ", kTH1F, {{300, 0., 30000.}}); // Track QA - registry.add("QATrack/hSelection", "Tracks statistics", kTH1D, {{10, 0.0f, 10.0f}}); + registry.add("QATrack/hSelection", "Tracks statistics", kTH1D, {{9, 0.0f, 9.0f}}); auto hTrack = registry.get(HIST("QATrack/hSelection")); hTrack->GetXaxis()->SetBinLabel(1, "all tracks"); hTrack->GetXaxis()->SetBinLabel(2, "passed pT cut"); hTrack->GetXaxis()->SetBinLabel(3, "passed eta cut"); hTrack->GetXaxis()->SetBinLabel(4, "passed DCA cut"); hTrack->GetXaxis()->SetBinLabel(5, "passed PID cut"); - hTrack->GetXaxis()->SetBinLabel(6, "passed rapidity cut"); - hTrack->GetXaxis()->SetBinLabel(7, "passed tpcNClsFound cut"); - hTrack->GetXaxis()->SetBinLabel(8, "passed isPrimaryTrack cut"); - hTrack->GetXaxis()->SetBinLabel(9, "passed isPVContributor cut"); - hTrack->GetXaxis()->SetBinLabel(10, "passed all cuts"); + hTrack->GetXaxis()->SetBinLabel(6, "passed tpcNClsFound cut"); + hTrack->GetXaxis()->SetBinLabel(7, "passed isPrimaryTrack cut"); + hTrack->GetXaxis()->SetBinLabel(8, "passed isPVContributor cut"); + hTrack->GetXaxis()->SetBinLabel(9, "passed all cuts"); hTrack->SetMinimum(0.1); registry.add("QATrack/hRapidity", "Rapidity distribution of K^{+} and K^{-}", kTH1F, {{200, -1, 1}}); @@ -369,18 +377,11 @@ struct PhianalysisTHnSparse { if (produceQA && dataQA) registry.fill(HIST("QATrack/hSelection"), 4.5); - // Apply rapidity cut - if (std::abs(track.rapidity(isPositive ? massPos : massNeg)) > static_cast(cut.rapidity)) { - return false; - } - if (produceQA && dataQA) - registry.fill(HIST("QATrack/hSelection"), 5.5); - // Apply tpcNClsFound cut if (track.tpcNClsFound() < tpcNClsFound) return false; if (produceQA && dataQA) - registry.fill(HIST("QATrack/hSelection"), 6.5); + registry.fill(HIST("QATrack/hSelection"), 5.5); if (globalTrack) { // Apply Global track cuts @@ -392,16 +393,16 @@ struct PhianalysisTHnSparse { return false; } if (produceQA && dataQA) - registry.fill(HIST("QATrack/hSelection"), 7.5); + registry.fill(HIST("QATrack/hSelection"), 6.5); // Apply PV Contributor cuts if (!track.isPVContributor()) return false; if (produceQA && dataQA) - registry.fill(HIST("QATrack/hSelection"), 8.5); + registry.fill(HIST("QATrack/hSelection"), 7.5); if (produceQA && dataQA) - registry.fill(HIST("QATrack/hSelection"), 9.5); + registry.fill(HIST("QATrack/hSelection"), 8.5); return true; } @@ -412,7 +413,7 @@ struct PhianalysisTHnSparse { d2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massNeg); mother = d1 + d2; - if (std::abs(mother.Eta()) > static_cast(cut.etapair)) + if (std::abs(mother.Rapidity()) > static_cast(cut.rapidity)) return false; return true; @@ -458,15 +459,16 @@ struct PhianalysisTHnSparse { registry.fill(HIST("QAEvent/hVtxZ"), collision.posZ()); registry.fill(HIST("QAEvent/hMult"), getMultiplicity(collision)); registry.fill(HIST("QAEvent/hCent"), getCentrality(collision)); - - dataQA = true; - for (const auto& track : posDaughters) { - selectedTrack(track, true); - } - for (const auto& track : negDaughters) { - selectedTrack(track, false); + if (produceStats) { + dataQA = true; + for (const auto& track : posDaughters) { + selectedTrack(track, true); + } + for (const auto& track : negDaughters) { + selectedTrack(track, false); + } + dataQA = false; } - dataQA = false; } if (static_cast(verbose.verboselevel) > 0 && static_cast(verbose.refresh) > 0 && collision.globalIndex() % static_cast(verbose.refresh) == static_cast(verbose.refreshIndex)) @@ -526,14 +528,15 @@ struct PhianalysisTHnSparse { 0); rsnOutput->fillUnlikepm(pointPair); - if (produceQA) + if (produceQA) { registry.fill(HIST("QAEvent/hSelection"), 2.5); - if (n == 0) - registry.fill(HIST("QAEvent/hSelection"), 1.5); + if (n == 0) + registry.fill(HIST("QAEvent/hSelection"), 1.5); + } n = n + 1; } - if (static_cast(produce.produceLikesign)) { + if (produceLikesign) { for (const auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posDaughters, posDaughters))) { if (!selectedTrack(track1, true)) // both positive @@ -591,6 +594,46 @@ struct PhianalysisTHnSparse { rsnOutput->fillLikemm(pointPair); } } + + if (produceRotational) { + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughters, negDaughters))) { + + if (!selectedTrack(track1, true)) + continue; + if (!selectedTrack(track2, false)) + continue; + + d1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massPos); + d2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massNeg); + mother = d1 + d2; + + if (std::abs(mother.Rapidity()) > static_cast(cut.rapidity)) + continue; + + for (int i = 1; i <= static_cast(numberofRotations); i++) { + float angle = i * (360.0f / (static_cast(numberofRotations) + 1)); + float px2new = track2.px() * std::cos(angle * TMath::DegToRad()) - track2.py() * std::sin(angle * TMath::DegToRad()); + float py2new = track2.px() * std::sin(angle * TMath::DegToRad()) + track2.py() * std::cos(angle * TMath::DegToRad()); + d2 = ROOT::Math::PxPyPzMVector(px2new, py2new, track2.pz(), massNeg); + mother = d1 + d2; + + pointPair = fillPointPair(mother.M(), + mother.Pt(), + getMultiplicity(collision), + getCentrality(collision), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), + mother.Eta(), + mother.Rapidity(), + collision.posZ(), + 0, + 0, + 0); + + rsnOutput->fillRotationpm(pointPair); + } + } + } } PROCESS_SWITCH(PhianalysisTHnSparse, processData, "Process Event for Data", true); @@ -772,10 +815,10 @@ struct PhianalysisTHnSparse { auto tracksTuple = std::make_tuple(tracks); BinningTypeVzCe binningVzCe{{axisVertexMixing, axisCentralityMixing}, true}; - SameKindPair pairVzCe{binningVzCe, numberofMixedEvents, -1, collisions, tracksTuple, &cache}; + SameKindPair pairVzCe{binningVzCe, static_cast(numberofMixedEvents), -1, collisions, tracksTuple, &cache}; BinningTypeVzMu binningVzMu{{axisVertexMixing, axisMultiplicityMixing}, true}; - SameKindPair pairVzMu{binningVzMu, numberofMixedEvents, -1, collisions, tracksTuple, &cache}; + SameKindPair pairVzMu{binningVzMu, static_cast(numberofMixedEvents), -1, collisions, tracksTuple, &cache}; if (mixingType == rsn::MixingType::ce) { for (const auto& [c1, tracks1, c2, tracks2] : pairVzCe) { @@ -819,60 +862,6 @@ struct PhianalysisTHnSparse { rsnOutput->fillMixingpm(pointPair); } - if (static_cast(produce.produceLikesign)) { - - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc1, posDaughtersc2))) { - - if (!selectedTrack(track1, true)) // track1 is positive - continue; - if (!selectedTrack(track2, true)) // track2 is positive - continue; - - if (!selectedPair(mother, track1, track2)) - continue; - - pointPair = fillPointPair(mother.M(), - mother.Pt(), - getMultiplicity(c1), - getCentrality(c1), - track1.tpcNSigmaKa(), - track2.tpcNSigmaKa(), - mother.Eta(), - mother.Rapidity(), - c1.posZ(), - getMultiplicity(c2), - getCentrality(c2), - c2.posZ()); - - rsnOutput->fillMixingpp(pointPair); - } - - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(negDaughtersc1, negDaughtersc2))) { - - if (!selectedTrack(track1, false)) - continue; - if (!selectedTrack(track2, false)) - continue; - - if (!selectedPair(mother, track1, track2)) - continue; - pointPair = fillPointPair(mother.M(), - mother.Pt(), - getMultiplicity(c1), - getCentrality(c1), - track1.tpcNSigmaKa(), - track2.tpcNSigmaKa(), - mother.Eta(), - mother.Rapidity(), - c1.posZ(), - getMultiplicity(c2), - getCentrality(c2), - c2.posZ()); - - rsnOutput->fillMixingmm(pointPair); - } - } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc2, negDaughtersc1))) { if (!selectedTrack(track1, true)) // track1 is positive @@ -943,62 +932,6 @@ struct PhianalysisTHnSparse { rsnOutput->fillMixingpm(pointPair); } - if (static_cast(produce.produceLikesign)) { - - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc1, posDaughtersc2))) { - - if (!selectedTrack(track1, true)) // track1 is positive - - continue; - if (!selectedTrack(track2, true)) // track2 is positive - continue; - - if (!selectedPair(mother, track1, track2)) - continue; - - pointPair = fillPointPair(mother.M(), - mother.Pt(), - getMultiplicity(c1), - getCentrality(c1), - track1.tpcNSigmaKa(), - track2.tpcNSigmaKa(), - mother.Eta(), - mother.Rapidity(), - c1.posZ(), - getMultiplicity(c2), - getCentrality(c2), - c2.posZ()); - - rsnOutput->fillMixingpp(pointPair); - } - - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(negDaughtersc1, negDaughtersc2))) { - - if (!selectedTrack(track1, false)) - - continue; - if (!selectedTrack(track2, false)) - continue; - - if (!selectedPair(mother, track1, track2)) - continue; - pointPair = fillPointPair(mother.M(), - mother.Pt(), - getMultiplicity(c1), - getCentrality(c1), - track1.tpcNSigmaKa(), - track2.tpcNSigmaKa(), - mother.Eta(), - mother.Rapidity(), - c1.posZ(), - getMultiplicity(c2), - getCentrality(c2), - c2.posZ()); - - rsnOutput->fillMixingmm(pointPair); - } - } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc2, negDaughtersc1))) { if (!selectedTrack(track1, true)) diff --git a/PWGLF/Utils/rsnOutput.h b/PWGLF/Utils/rsnOutput.h index 592eaef40cc..f6540cd6d48 100644 --- a/PWGLF/Utils/rsnOutput.h +++ b/PWGLF/Utils/rsnOutput.h @@ -16,13 +16,13 @@ #ifndef PWGLF_UTILS_RSNOUTPUT_H_ #define PWGLF_UTILS_RSNOUTPUT_H_ -#include -#include -#include - #include "Framework/HistogramRegistry.h" #include "Framework/Logger.h" +#include +#include +#include + namespace o2::analysis { namespace rsn @@ -48,9 +48,8 @@ enum class PairType { unlikegen, unlikegenold, mixingpm, - mixingpp, - mixingmm, mixingmp, + rotationpm, all }; @@ -105,7 +104,7 @@ class Output public: virtual ~Output() = default; - virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, std::vector const& sysAxes, std::vector const& allAxes_sys, bool /*produceTrue*/ = false, MixingType /*eventMixing*/ = MixingType::none, bool /*produceLikesign*/ = false, o2::framework::HistogramRegistry* registry = nullptr) + virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, std::vector const& sysAxes, std::vector const& allAxes_sys, bool /*produceTrue*/ = false, MixingType /*eventMixing*/ = MixingType::none, bool /*produceLikesign*/ = false, bool /*produceRotational*/ = false, o2::framework::HistogramRegistry* registry = nullptr) { mHistogramRegistry = registry; if (mHistogramRegistry == nullptr) @@ -209,9 +208,8 @@ class Output virtual void fillUnlikegen(double* point) = 0; virtual void fillUnlikegenOld(double* point) = 0; virtual void fillMixingpm(double* point) = 0; - virtual void fillMixingpp(double* point) = 0; - virtual void fillMixingmm(double* point) = 0; virtual void fillMixingmp(double* point) = 0; + virtual void fillRotationpm(double* point) = 0; virtual void fillSystematics(double* point) = 0; PairAxisType type(std::string name) @@ -265,12 +263,11 @@ class Output class OutputSparse : public Output { public: - virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, std::vector const& sysAxes, std::vector const& allAxes_sys, bool produceTrue = false, MixingType eventMixing = MixingType::none, bool produceLikesign = false, o2::framework::HistogramRegistry* registry = nullptr) + virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, std::vector const& sysAxes, std::vector const& allAxes_sys, bool produceTrue = false, MixingType eventMixing = MixingType::none, bool produceLikesign = false, bool produceRotational = false, o2::framework::HistogramRegistry* registry = nullptr) { - Output::init(sparseAxes, allAxes, sysAxes, allAxes_sys, produceTrue, eventMixing, produceLikesign, registry); + Output::init(sparseAxes, allAxes, sysAxes, allAxes_sys, produceTrue, eventMixing, produceLikesign, produceRotational, registry); mHistogramRegistry->add("unlikepm", "Unlike pm", *mPairHisto); - // mHistogramRegistry->add("unlikemp", "Unlike mp", *mPairHisto); if (produceLikesign) { mHistogramRegistry->add("likepp", "Like PP", *mPairHisto); mHistogramRegistry->add("likemm", "Like MM", *mPairHisto); @@ -282,13 +279,11 @@ class OutputSparse : public Output } if (eventMixing != MixingType::none) { mHistogramRegistry->add("mixingpm", "Event Mixing pm", *mPairHisto); - if (produceLikesign) { - mHistogramRegistry->add("mixingpp", "Event Mixing pp", *mPairHisto); - mHistogramRegistry->add("mixingmm", "Event Mixing mm", *mPairHisto); - } mHistogramRegistry->add("mixingmp", "Event Mixing mp", *mPairHisto); } - + if (produceRotational) { + mHistogramRegistry->add("rotationpm", "Rotational pm", *mPairHisto); + } mHistogramRegistry->add("Mapping/systematics", "Systematics mapping", *mPairHistoSys); } @@ -331,15 +326,12 @@ class OutputSparse : public Output case PairType::mixingpm: fillMixingpm(point); break; - case PairType::mixingpp: - fillMixingpp(point); - break; - case PairType::mixingmm: - fillMixingmm(point); - break; case PairType::mixingmp: fillMixingmp(point); break; + case PairType::rotationpm: + fillRotationpm(point); + break; default: break; } @@ -357,17 +349,14 @@ class OutputSparse : public Output { fillSparse(HIST("likepp"), point); } - virtual void fillLikemm(double* point) { fillSparse(HIST("likemm"), point); } - virtual void fillUnliketrue(double* point) { fillSparse(HIST("unliketrue"), point); } - virtual void fillUnlikegen(double* point) { fillSparse(HIST("unlikegen"), point); @@ -380,18 +369,14 @@ class OutputSparse : public Output { fillSparse(HIST("mixingpm"), point); } - virtual void fillMixingpp(double* point) - { - fillSparse(HIST("mixingpp"), point); - } - virtual void fillMixingmm(double* point) - { - fillSparse(HIST("mixingmm"), point); - } virtual void fillMixingmp(double* point) { fillSparse(HIST("mixingmp"), point); } + virtual void fillRotationpm(double* point) + { + fillSparse(HIST("rotationpm"), point); + } virtual void fillSystematics(double* point) { fillSparse(HIST("Mapping/systematics"), point);