Skip to content

Commit c898fc6

Browse files
jaelparkjaelpark
committed
Multiplicity correlations and outlier cuts
1 parent 894434d commit c898fc6

File tree

1 file changed

+17
-18
lines changed

1 file changed

+17
-18
lines changed

‎PWGCF/Tasks/correlations.cxx‎

Lines changed: 17 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@
4343

4444
#include <cmath>
4545
#include <experimental/type_traits>
46+
#include <memory>
4647
#include <string>
4748
#include <tuple>
4849
#include <utility>
@@ -127,7 +128,7 @@ struct CorrelationTask {
127128
Filter collisionVertexTypeFilter = (aod::collision::flags & static_cast<uint16_t>(aod::collision::CollisionFlagsRun2::Run2VertexerTracks)) == static_cast<uint16_t>(aod::collision::CollisionFlagsRun2::Run2VertexerTracks);
128129

129130
// Track filters
130-
Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true));
131+
Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true));
131132
Filter cfTrackFilter = (nabs(aod::cftrack::eta) < cfgCutEta) && (aod::cftrack::pt > cfgCutPt) && ((aod::track::trackType & (uint8_t)cfgTrackBitMask) == (uint8_t)cfgTrackBitMask);
132133

133134
// MC filters
@@ -190,11 +191,11 @@ struct CorrelationTask {
190191
if (cfgMultCorrelationsMask & aod::cfmultset::CentFT0C)
191192
multAxes.emplace_back(100, 0, 100, "FT0C centrality");
192193
if (cfgMultCorrelationsMask & aod::cfmultset::MultFV0A)
193-
multAxes.emplace_back(100, 0, 100000, "V0A multiplicity");
194+
multAxes.emplace_back(10000, 0, 100000, "V0A multiplicity");
194195
if (cfgMultCorrelationsMask & aod::cfmultset::MultNTracksPV)
195-
multAxes.emplace_back(100, 0, 1000, "Nch PV");
196+
multAxes.emplace_back(1000, 0, 1000, "Nch PV");
196197
if (cfgMultCorrelationsMask & aod::cfmultset::MultNTracksGlobal)
197-
multAxes.emplace_back(100, 0, 1000, "Nch Global");
198+
multAxes.emplace_back(1000, 0, 1000, "Nch Global");
198199
registry.add("multCorrelations", "Multiplicity correlations", {HistType::kTHnSparseF, multAxes});
199200
}
200201
registry.add("multiplicity", "event multiplicity", {HistType::kTH1F, {{1000, 0, 100, "/multiplicity/centrality"}}});
@@ -227,18 +228,20 @@ struct CorrelationTask {
227228
if (!cfgMultCutFormula.value.empty()) {
228229
multCutFormula = std::make_unique<TFormula>("multCutFormula", cfgMultCutFormula.value.c_str());
229230
std::fill_n(multCutFormulaParamIndex.begin(), std::size(multCutFormulaParamIndex), ~0u);
230-
std::array<std::string, 4> pars = {"cFT0C", "mFV0A", "mGlob", "mPV"};
231-
for (uint i = 0; i < multCutFormula->GetNpar(); ++i) {
231+
std::array<std::string, 4> pars = {"cFT0C", "mFV0A", "mPV", "mGlob"}; // must correspond the order of MultiplicityEstimators
232+
for (uint i = 0, n = multCutFormula->GetNpar(); i < n; ++i) {
232233
auto m = std::find(pars.begin(), pars.end(), multCutFormula->GetParName(i));
233234
if (m == pars.end()) {
234235

235-
LOGF(warning, "Unknown parameter cfgMultCutFormula: %s", multCutFormula->GetParName(i));
236+
LOGF(warning, "Unknown parameter in cfgMultCutFormula: %s", multCutFormula->GetParName(i));
236237
continue;
237238
}
238-
if (cfgMultCorrelationsMask.value & (1u << i) == 0)
239+
if ((cfgMultCorrelationsMask.value & (1u << i)) == 0) {
239240
LOGF(warning, "The centrality/multiplicity estimator %s is not available to be used in cfgMultCutFormula. Ensure cfgMultCorrelationsMask is correct and matches the CFMultSets in derived data.");
240-
else
241+
} else {
241242
multCutFormulaParamIndex[std::distance(pars.begin(), m)] = i;
243+
LOGF(info, "Multiplicity cut parameter %s in use.", m->c_str());
244+
}
242245
}
243246
}
244247

@@ -451,16 +454,16 @@ struct CorrelationTask {
451454
using HasPartDaugh1Id = decltype(std::declval<T&>().cfParticleDaugh1Id());
452455

453456
/*
454-
OO outlier cut:
455-
(567.785+172.715*[mGlob]+0.77888*[mGlob]*[mGlob]+-0.00693466*[mGlob]*[mGlob]*[mGlob]+1.40564e-05*[mGlob]*[mGlob]*[mGlob]*[mGlob] + 3.5*(679.853+66.8068*[mGlob]+-0.444332*[mGlob]*[mGlob]+0.00115002*[mGlob]*[mGlob]*[mGlob]+-4.92064e-07*[mGlob]*[mGlob]*[mGlob]*[mGlob])) < [mFV0A] &&(567.785+172.715*[mGlob]+0.77888*[mGlob]*[mGlob]+-0.00693466*[mGlob]*[mGlob]*[mGlob]+1.40564e-05*[mGlob]*[mGlob]*[mGlob]*[mGlob] - 3.0*(679.853+66.8068*[mGlob]+-0.444332*[mGlob]*[mGlob]+0.00115002*[mGlob]*[mGlob]*[mGlob]+-4.92064e-07*[mGlob]*[mGlob]*[mGlob]*[mGlob])) > [mFV0A] &&(172.406 + -4.50219*[cFT0C] + 0.0543038*[cFT0C]*[cFT0C] + -0.000373213*[cFT0C]*[cFT0C]*[cFT0C] + 1.15322e-06*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] + 4.0*(49.7503 + -1.29008*[cFT0C] + 0.0160059*[cFT0C]*[cFT0C] + -7.86846e-05*[cFT0C]*[cFT0C]*[cFT0C])) < [mPV] &&(172.406 + -4.50219*[cFT0C] + 0.0543038*[cFT0C]*[cFT0C] + -0.000373213*[cFT0C]*[cFT0C]*[cFT0C] + 1.15322e-06*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] - 2.5*(49.7503 + -1.29008*[cFT0C] + 0.0160059*[cFT0C]*[cFT0C] + -7.86846e-05*[cFT0C]*[cFT0C]*[cFT0C])) > [mPV] &&(125.02 + -3.30255*[cFT0C] + 0.0398663*[cFT0C]*[cFT0C] + -0.000271942*[cFT0C]*[cFT0C]*[cFT0C] + 8.34098e-07*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] + 4.0*(37.0244 + -0.949883*[cFT0C] + 0.0116622*[cFT0C]*[cFT0C] + -5.71117e-05*[cFT0C]*[cFT0C]*[cFT0C])) < [mGlob] &&(125.02 + -3.30255*[cFT0C] + 0.0398663*[cFT0C]*[cFT0C] + -0.000271942*[cFT0C]*[cFT0C]*[cFT0C] + 8.34098e-07*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] - 2.5*(37.0244 + -0.949883*[cFT0C] + 0.0116622*[cFT0C]*[cFT0C] + -5.71117e-05*[cFT0C]*[cFT0C]*[cFT0C])) > [mGlob] &&(-0.223013 + 0.715849*[mPV] + 3*(0.664242 + 0.0829653*[mPV] + -0.000503733*[mPV]*[mPV] + 1.21185e-06*[mPV]*[mPV]*[mPV])) < [mGlob]
457+
OO outlier cut (requires mask 15):
458+
(567.785+172.715*[mGlob]+0.77888*[mGlob]*[mGlob]+-0.00693466*[mGlob]*[mGlob]*[mGlob]+1.40564e-05*[mGlob]*[mGlob]*[mGlob]*[mGlob] + 3.5*(679.853+66.8068*[mGlob]+-0.444332*[mGlob]*[mGlob]+0.00115002*[mGlob]*[mGlob]*[mGlob]+-4.92064e-07*[mGlob]*[mGlob]*[mGlob]*[mGlob])) > [mFV0A] && (567.785+172.715*[mGlob]+0.77888*[mGlob]*[mGlob]+-0.00693466*[mGlob]*[mGlob]*[mGlob]+1.40564e-05*[mGlob]*[mGlob]*[mGlob]*[mGlob] - 3.0*(679.853+66.8068*[mGlob]+-0.444332*[mGlob]*[mGlob]+0.00115002*[mGlob]*[mGlob]*[mGlob]+-4.92064e-07*[mGlob]*[mGlob]*[mGlob]*[mGlob])) < [mFV0A] && (172.406 + -4.50219*[cFT0C] + 0.0543038*[cFT0C]*[cFT0C] + -0.000373213*[cFT0C]*[cFT0C]*[cFT0C] + 1.15322e-06*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] + 4.0*(49.7503 + -1.29008*[cFT0C] + 0.0160059*[cFT0C]*[cFT0C] + -7.86846e-05*[cFT0C]*[cFT0C]*[cFT0C])) > [mPV] && (172.406 + -4.50219*[cFT0C] + 0.0543038*[cFT0C]*[cFT0C] + -0.000373213*[cFT0C]*[cFT0C]*[cFT0C] + 1.15322e-06*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] - 2.5*(49.7503 + -1.29008*[cFT0C] + 0.0160059*[cFT0C]*[cFT0C] + -7.86846e-05*[cFT0C]*[cFT0C]*[cFT0C])) < [mPV] && (125.02 + -3.30255*[cFT0C] + 0.0398663*[cFT0C]*[cFT0C] + -0.000271942*[cFT0C]*[cFT0C]*[cFT0C] + 8.34098e-07*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] + 4.0*(37.0244 + -0.949883*[cFT0C] + 0.0116622*[cFT0C]*[cFT0C] + -5.71117e-05*[cFT0C]*[cFT0C]*[cFT0C])) > [mGlob] && (125.02 + -3.30255*[cFT0C] + 0.0398663*[cFT0C]*[cFT0C] + -0.000271942*[cFT0C]*[cFT0C]*[cFT0C] + 8.34098e-07*[cFT0C]*[cFT0C]*[cFT0C]*[cFT0C] - 2.5*(37.0244 + -0.949883*[cFT0C] + 0.0116622*[cFT0C]*[cFT0C] + -5.71117e-05*[cFT0C]*[cFT0C]*[cFT0C])) < [mGlob] && (-0.223013 + 0.715849*[mPV] + 3*(0.664242 + 0.0829653*[mPV] + -0.000503733*[mPV]*[mPV] + 1.21185e-06*[mPV]*[mPV]*[mPV])) > [mGlob]
456459
*/
457460
template <class CollType>
458461
bool passOutlier(CollType const& collision)
459462
{
460463
if (cfgMultCutFormula.value.empty())
461464
return true;
462465
for (uint i = 0; i < 4; ++i) {
463-
if ((cfgMultCorrelationsMask.value & (1u << i)) == 0)
466+
if ((cfgMultCorrelationsMask.value & (1u << i)) == 0 || multCutFormulaParamIndex[i] == ~0u)
464467
continue;
465468
auto estIndex = std::popcount(cfgMultCorrelationsMask.value & ((1u << i) - 1));
466469
multCutFormula->SetParameter(multCutFormulaParamIndex[i], collision.multiplicities()[estIndex]);
@@ -809,11 +812,10 @@ struct CorrelationTask {
809812
}
810813
PROCESS_SWITCH(CorrelationTask, processSameAOD, "Process same event on AOD", true);
811814

812-
using BinningTypeDerived = ColumnBinningPolicy<aod::collision::PosZ, aod::cfcollision::Multiplicity>;
813-
814815
template <class CollType, class TTracks1, class TTracks2>
815816
void processSameDerivedT(CollType const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2)
816817
{
818+
using BinningTypeDerived = ColumnBinningPolicy<aod::collision::PosZ, aod::cfcollision::Multiplicity>;
817819
BinningTypeDerived configurableBinningDerived{{axisVertex, axisMultiplicity}, true}; // true is for 'ignore overflows' (true by default). Underflows and overflows will have bin -1.
818820
if (cfgVerbosity > 0) {
819821
LOGF(info, "processSameDerivedT: Tracks for collision: %d/%d | Vertex: %.1f | Multiplicity/Centrality: %.1f", tracks1.size(), tracks2.size(), collision.posZ(), collision.multiplicity());
@@ -937,8 +939,7 @@ struct CorrelationTask {
937939
};
938940

939941
using BinningTypeDerived = FlexibleBinningPolicy<std::tuple<decltype(getMultiplicity)>, aod::collision::PosZ, decltype(getMultiplicity)>;
940-
BinningTypeDerived configurableBinningDerived{{getMultiplicity}, {axisVertex, axisMultiplicity}, true};
941-
// BinningTypeDerived configurableBinningDerived{{axisVertex, axisMultiplicity}, true}; // true is for 'ignore overflows' (true by default). Underflows and overflows will have bin -1.
942+
BinningTypeDerived configurableBinningDerived{{getMultiplicity}, {axisVertex, axisMultiplicity}, true}; // true is for 'ignore overflows' (true by default). Underflows and overflows will have bin -1.
942943
// Strictly upper categorised collisions, for cfgNoMixedEvents combinations per bin, skipping those in entry -1
943944
auto tracksTuple = std::make_tuple(std::forward<TrackTypes>(tracks)...);
944945
using TA = std::tuple_element<0, decltype(tracksTuple)>::type;
@@ -947,7 +948,6 @@ struct CorrelationTask {
947948

948949
for (auto it = pairs.begin(); it != pairs.end(); it++) {
949950
auto& [collision1, tracks1, collision2, tracks2] = *it;
950-
// int bin = configurableBinningDerived.getBin({collision1.posZ(), collision1.multiplicity()});
951951
float multiplicity = getMultiplicity(collision1);
952952
int bin = configurableBinningDerived.getBin(std::tuple(collision1.posZ(), multiplicity));
953953
float eventWeight = 1.0f / it.currentWindowNeighbours();
@@ -989,7 +989,6 @@ struct CorrelationTask {
989989

990990
void processMixedDerivedMultSet(soa::Filtered<soa::Join<aod::CFCollisions, aod::CFMultSets>> const& collisions, DerivedTracks const& tracks)
991991
{
992-
// TODO: use FlexibleBinningPolicy to reject (by returning -1 from lambda) outliers
993992
processMixedDerivedT(collisions, tracks);
994993
}
995994
PROCESS_SWITCH(CorrelationTask, processMixedDerivedMultSet, "Process mixed events on derived data with multiplicity sets", false);

0 commit comments

Comments
 (0)