From e5d2148a7c2135eecb3cfc0c05b308a89389262e Mon Sep 17 00:00:00 2001 From: Jan Fiete Date: Tue, 3 Nov 2026 10:52:55 +0100 Subject: [PATCH] Revert "[PWGCF] Add MFT-MFT correlation task (#18485)" This reverts commit b533857c393e39ec85a3cebb7d0c2e5bd7b71e18. --- PWGCF/Core/CorrelationContainer.cxx | 360 ++-- .../Tasks/CMakeLists.txt | 5 - .../Tasks/mftcorrelations.cxx | 1449 ----------------- 3 files changed, 178 insertions(+), 1636 deletions(-) delete mode 100644 PWGCF/TwoParticleCorrelations/Tasks/mftcorrelations.cxx diff --git a/PWGCF/Core/CorrelationContainer.cxx b/PWGCF/Core/CorrelationContainer.cxx index e60bdfb219e..6c4a2e56d0d 100644 --- a/PWGCF/Core/CorrelationContainer.cxx +++ b/PWGCF/Core/CorrelationContainer.cxx @@ -24,12 +24,11 @@ #include #include #include -#include -#include -#include #include #include #include +#include +#include #include #include #include @@ -37,8 +36,6 @@ #include #include -#include -#include #include #include #include @@ -49,7 +46,7 @@ using namespace o2::constants::math; ClassImp(CorrelationContainer); -const int CorrelationContainer::fgkCFSteps = 11; +const Int_t CorrelationContainer::fgkCFSteps = 11; CorrelationContainer::CorrelationContainer() : TNamed(), mPairHist(nullptr), @@ -67,9 +64,9 @@ CorrelationContainer::CorrelationContainer() : TNamed(), mZVtxMax(0), mTrackEtaCut(0), mWeightPerEvent(0), - mSkipScaleMixedEvent(false), + mSkipScaleMixedEvent(kFALSE), mCache(nullptr), - mGetMultCacheOn(false), + mGetMultCacheOn(kFALSE), mGetMultCache(nullptr) { // Default constructor @@ -94,9 +91,9 @@ CorrelationContainer::CorrelationContainer(const char* name, const char* objTitl mZVtxMax(0), mTrackEtaCut(0), mWeightPerEvent(0), - mSkipScaleMixedEvent(false), + mSkipScaleMixedEvent(kFALSE), mCache(nullptr), - mGetMultCacheOn(false), + mGetMultCacheOn(kFALSE), mGetMultCache(nullptr) { // correlationAxis has to provide a 6 length list of AxisSpec which contain: @@ -128,7 +125,7 @@ CorrelationContainer::CorrelationContainer(const char* name, const char* objTitl triggerAxis.insert(triggerAxis.end(), userAxis.begin(), userAxis.end()); mTriggerHist = HistFactory::createHist({"mTriggerHist", "d^{2}N_{ch}/d#varphid#eta", {HistType::kStepTHnF, triggerAxis, fgkCFSteps}}).release(); - mTrackHistEfficiency = HistFactory::createHist({"mTrackHistEfficiency", "Tracking efficiency", {HistType::kStepTHnF, {efficiencyAxis[0], efficiencyAxis[1], {5, -0.5, 4.5, "species"}, correlationAxis[3], efficiencyAxis[2], efficiencyAxis[3]}, fgkCFSteps}}).release(); + mTrackHistEfficiency = HistFactory::createHist({"mTrackHistEfficiency", "Tracking efficiency", {HistType::kStepTHnF, {efficiencyAxis[0], efficiencyAxis[1], {5, -0.5, 4.5, "species"}, correlationAxis[3], efficiencyAxis[2]}, fgkCFSteps}}).release(); mEventCount = HistFactory::createHist({"mEventCount", ";step;centrality;count", {HistType::kTH2F, {{fgkCFSteps + 2, -2.5, -0.5 + fgkCFSteps, "step"}, correlationAxis[3]}}}).release(); } @@ -149,16 +146,16 @@ CorrelationContainer::CorrelationContainer(const CorrelationContainer& c) : TNam mZVtxMax(0), mTrackEtaCut(0), mWeightPerEvent(0), - mSkipScaleMixedEvent(false), + mSkipScaleMixedEvent(kFALSE), mCache(nullptr), - mGetMultCacheOn(false), + mGetMultCacheOn(kFALSE), mGetMultCache(nullptr) { // // CorrelationContainer copy constructor // - c.Copy(*this); + ((CorrelationContainer&)c).Copy(*this); } //____________________________________________________________________ @@ -198,7 +195,7 @@ CorrelationContainer& CorrelationContainer::operator=(const CorrelationContainer // assigment operator if (this != &c) { - c.Copy(*this); + ((CorrelationContainer&)c).Copy(*this); } return *this; @@ -209,7 +206,7 @@ void CorrelationContainer::Copy(TObject& c) const { // copy function - auto& target = dynamic_cast(c); + CorrelationContainer& target = (CorrelationContainer&)c; if (mPairHist) { target.mPairHist = dynamic_cast(mPairHist->Clone()); @@ -256,14 +253,14 @@ Long64_t CorrelationContainer::Merge(TCollection* list) TObject* obj = nullptr; // collections of objects - const unsigned int kMaxLists = 4; + const UInt_t kMaxLists = 4; TList** lists = new TList*[kMaxLists]; - for (unsigned int i = 0; i < kMaxLists; i++) { + for (UInt_t i = 0; i < kMaxLists; i++) { lists[i] = new TList; } - int count = 0; + Int_t count = 0; while ((obj = iter->Next())) { CorrelationContainer* entry = dynamic_cast(obj); @@ -289,7 +286,7 @@ Long64_t CorrelationContainer::Merge(TCollection* list) mTrackHistEfficiency->Merge(lists[2]); mEventCount->Merge(lists[3]); - for (unsigned int i = 0; i < kMaxLists; i++) { + for (UInt_t i = 0; i < kMaxLists; i++) { delete lists[i]; } @@ -320,7 +317,7 @@ void CorrelationContainer::resetBinLimits(THnBase* grid, int max_dimension) max_dimension = grid->GetNdimensions(); } - for (int i = 0; i < max_dimension; i++) { + for (Int_t i = 0; i < max_dimension; i++) { if (grid->GetAxis(i)->TestBit(TAxis::kAxisRange)) { grid->GetAxis(i)->SetRange(0, 0); // reset range } @@ -328,14 +325,14 @@ void CorrelationContainer::resetBinLimits(THnBase* grid, int max_dimension) } //____________________________________________________________________ -void CorrelationContainer::countEmptyBins(CorrelationContainer::CFStep step, float ptTriggerMin, float ptTriggerMax) +void CorrelationContainer::countEmptyBins(CorrelationContainer::CFStep step, Float_t ptTriggerMin, Float_t ptTriggerMax) { // prints the number of empty bins in the track end event histograms in the given step - int binBegin[4]; - int binEnd[4]; + Int_t binBegin[4]; + Int_t binEnd[4]; - for (int i = 0; i < 4; i++) { + for (Int_t i = 0; i < 4; i++) { binBegin[i] = 1; binEnd[i] = mPairHist->getTHn(step)->GetAxis(i)->GetNbins(); } @@ -358,11 +355,11 @@ void CorrelationContainer::countEmptyBins(CorrelationContainer::CFStep step, flo // start from multiplicity 1 binBegin[3] = mPairHist->getTHn(step)->GetAxis(3)->FindBin(1); - int total = 0; - int count = 0; - int vars[4]; + Int_t total = 0; + Int_t count = 0; + Int_t vars[4]; - for (int i = 0; i < 4; i++) { + for (Int_t i = 0; i < 4; i++) { vars[i] = binBegin[i]; } @@ -378,7 +375,7 @@ void CorrelationContainer::countEmptyBins(CorrelationContainer::CFStep step, flo } vars[3]++; - for (int i = 3; i > 0; i--) { + for (Int_t i = 3; i > 0; i--) { if (vars[i] == binEnd[i] + 1) { vars[i] = binBegin[i]; vars[i - 1]++; @@ -395,7 +392,7 @@ void CorrelationContainer::countEmptyBins(CorrelationContainer::CFStep step, flo } //____________________________________________________________________ -void CorrelationContainer::getHistsZVtxMult(CorrelationContainer::CFStep step, float ptTriggerMin, float ptTriggerMax, THnBase** trackHist, TH2** eventHist) +void CorrelationContainer::getHistsZVtxMult(CorrelationContainer::CFStep step, Float_t ptTriggerMin, Float_t ptTriggerMax, THnBase** trackHist, TH2** eventHist) { // Calculates a 4d histogram with deltaphi, deltaeta, zvtx, multiplicity on track level and // a 2d histogram on event level (as fct of zvtx, multiplicity) @@ -418,15 +415,15 @@ void CorrelationContainer::getHistsZVtxMult(CorrelationContainer::CFStep step, f sparse->GetAxis(1)->SetRangeUser(mPtMin, mPtMax); } - int firstBin = sparse->GetAxis(2)->FindBin(ptTriggerMin); - int lastBin = sparse->GetAxis(2)->FindBin(ptTriggerMax); + Int_t firstBin = sparse->GetAxis(2)->FindBin(ptTriggerMin); + Int_t lastBin = sparse->GetAxis(2)->FindBin(ptTriggerMax); LOGF(info, "Using trigger pT range %d --> %d", firstBin, lastBin); sparse->GetAxis(2)->SetRange(firstBin, lastBin); mTriggerHist->getTHn(step)->GetAxis(0)->SetRange(firstBin, lastBin); - int dimensions[] = {4, 0, 5, 3}; + Int_t dimensions[] = {4, 0, 5, 3}; THnBase* tmpTrackHist = sparse->ProjectionND(4, dimensions, "E"); - *eventHist = static_cast(mTriggerHist->getTHn(step)->Projection(1, 2)); // NOTE the syntax is Projection(Y, X) --> produces a histogram where x = axis 2 (vertex) and y = axis 1 (multiplicity) + *eventHist = (TH2*)mTriggerHist->getTHn(step)->Projection(1, 2); // NOTE the syntax is Projection(Y, X) --> produces a histogram where x = axis 2 (vertex) and y = axis 1 (multiplicity) // convert to THn *trackHist = changeToThn(tmpTrackHist); delete tmpTrackHist; @@ -435,7 +432,7 @@ void CorrelationContainer::getHistsZVtxMult(CorrelationContainer::CFStep step, f resetBinLimits(mTriggerHist->getTHn(step), 3); } -TH2* CorrelationContainer::getPerTriggerYield(CorrelationContainer::CFStep step, float ptTriggerMin, float ptTriggerMax, bool normalizePerTrigger) +TH2* CorrelationContainer::getPerTriggerYield(CorrelationContainer::CFStep step, Float_t ptTriggerMin, Float_t ptTriggerMax, Bool_t normalizePerTrigger) { // Calculate per trigger yield without considering mixed event // Sums over all vertex bins, and respects the pT and centrality limits set in the class @@ -475,7 +472,7 @@ TH2* CorrelationContainer::getPerTriggerYield(CorrelationContainer::CFStep step, } TH2* yield = trackSameAll->Projection(1, 0, "E"); - float triggers = eventSameAll->Integral(vertexBinBegin, vertexBinEnd, multBinBegin, multBinEnd); + Float_t triggers = eventSameAll->Integral(vertexBinBegin, vertexBinEnd, multBinBegin, multBinEnd); if (normalizePerTrigger) { LOGF(info, "Dividing %f tracks by %f triggers", yield->Integral(), triggers); @@ -485,7 +482,7 @@ TH2* CorrelationContainer::getPerTriggerYield(CorrelationContainer::CFStep step, } // normalizate to dphi width - float normalization = yield->GetXaxis()->GetBinWidth(1); + Float_t normalization = yield->GetXaxis()->GetBinWidth(1); yield->Scale(1.0 / normalization); delete trackSameAll; @@ -494,7 +491,7 @@ TH2* CorrelationContainer::getPerTriggerYield(CorrelationContainer::CFStep step, return yield; } -TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, CorrelationContainer::CFStep step, float ptTriggerMin, float ptTriggerMax, bool normalizePerTrigger, int stepForMixed, int* trigger) +TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, CorrelationContainer::CFStep step, Float_t ptTriggerMin, Float_t ptTriggerMax, Bool_t normalizePerTrigger, Int_t stepForMixed, Int_t* trigger) { // Extract 2D per trigger yield with mixed event correction. The quantity is calculated for *each* vertex bin and multiplicity bin and then a sum of ratios is performed: // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ] @@ -510,8 +507,8 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati // normalizePerTrigger: divide through number of triggers // do not add this hists to the directory - bool oldStatus = TH1::AddDirectoryStatus(); - TH1::AddDirectory(false); + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); TH2* totalTracks = nullptr; @@ -522,11 +519,11 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati TH2* eventMixedAll = nullptr; TH2* eventMixedAllStep6 = nullptr; - int64_t totalEvents = 0; - int nCorrelationFunctions = 0; + Long64_t totalEvents = 0; + Int_t nCorrelationFunctions = 0; getHistsZVtxMult(step, ptTriggerMin, ptTriggerMax, &trackSameAll, &eventSameAll); - mixed->getHistsZVtxMult((stepForMixed == -1) ? step : static_cast(stepForMixed), ptTriggerMin, ptTriggerMax, &trackMixedAll, &eventMixedAll); + mixed->getHistsZVtxMult((stepForMixed == -1) ? step : (CFStep)stepForMixed, ptTriggerMin, ptTriggerMax, &trackMixedAll, &eventMixedAll); // If we ask for histograms from step8 (TTR cut applied) there is a hole at 0,0; so this cannot be used for the // mixed-event normalization. If step6 is available, the normalization factor is read out from that one. @@ -552,15 +549,15 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati multBinEnd = multAxis->FindBin(mCentralityMax); } - for (int multBin = std::max(1, multBinBegin); multBin <= std::min(multAxis->GetNbins(), multBinEnd); multBin++) { + for (Int_t multBin = TMath::Max(1, multBinBegin); multBin <= TMath::Min(multAxis->GetNbins(), multBinEnd); multBin++) { trackSameAll->GetAxis(3)->SetRange(multBin, multBin); trackMixedAll->GetAxis(3)->SetRange(multBin, multBin); if (trackMixedAllStep6) { trackMixedAllStep6->GetAxis(3)->SetRange(multBin, multBin); } - double mixedNorm = 1; - double mixedNormError = 0; + Double_t mixedNorm = 1; + Double_t mixedNormError = 0; if (!mSkipScaleMixedEvent) { // get mixed normalization correction factor: is independent of vertex bin if scaled with number of triggers @@ -573,40 +570,40 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati tracksMixed = trackMixedAll->Projection(1, 0, "E"); } // Printf("%f", tracksMixed->Integral()); - float binWidthEta = tracksMixed->GetYaxis()->GetBinWidth(1); + Float_t binWidthEta = tracksMixed->GetYaxis()->GetBinWidth(1); if (stepForMixed == -1 && step == kCFStepBiasStudy && !trackMixedAllStep6) { // get mixed event normalization by assuming full acceptance at deta at 0 (integrate over dphi), excluding (0, 0) - float phiExclude = 0.41; + Float_t phiExclude = 0.41; mixedNorm = tracksMixed->IntegralAndError(1, tracksMixed->GetXaxis()->FindBin(-phiExclude) - 1, tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError); - double mixedNormError2 = 0; - double mixedNorm2 = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(phiExclude) + 1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError2); + Double_t mixedNormError2 = 0; + Double_t mixedNorm2 = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(phiExclude) + 1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError2); if (mixedNormError == 0 || mixedNormError2 == 0) { LOGF(error, "ERROR: Skipping multiplicity %d because mixed event is empty %f %f %f %f", multBin, mixedNorm, mixedNormError, mixedNorm2, mixedNormError2); continue; } - int nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(-phiExclude) - 1 - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1); + Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(-phiExclude) - 1 - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1); mixedNorm /= nBinsMixedNorm; mixedNormError /= nBinsMixedNorm; - int nBinsMixedNorm2 = (tracksMixed->GetNbinsX() - tracksMixed->GetXaxis()->FindBin(phiExclude) - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1); + Int_t nBinsMixedNorm2 = (tracksMixed->GetNbinsX() - tracksMixed->GetXaxis()->FindBin(phiExclude) - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1); mixedNorm2 /= nBinsMixedNorm2; mixedNormError2 /= nBinsMixedNorm2; mixedNorm = mixedNorm / mixedNormError / mixedNormError + mixedNorm2 / mixedNormError2 / mixedNormError2; - mixedNormError = std::sqrt(1.0 / (1.0 / mixedNormError / mixedNormError + 1.0 / mixedNormError2 / mixedNormError2)); + mixedNormError = TMath::Sqrt(1.0 / (1.0 / mixedNormError / mixedNormError + 1.0 / mixedNormError2 / mixedNormError2)); mixedNorm *= mixedNormError * mixedNormError; } else { // get mixed event normalization at (0,0) // NOTE if variable bin size is used around (0,0) to reduce two-track effect to limited bins, the normalization gets a bit tricky here (finite bin correction and normalization are made only for fixed size bins). // The normalization factor has to determined in a bin as large as the normal bin size, as a proxy the bin with index (1, 1) is used - float binWidthPhi = tracksMixed->GetXaxis()->GetBinWidth(1); + Float_t binWidthPhi = tracksMixed->GetXaxis()->GetBinWidth(1); mixedNorm = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(-binWidthPhi + 1e-4), tracksMixed->GetXaxis()->FindBin(binWidthPhi - 1e-4), tracksMixed->GetYaxis()->FindBin(-binWidthEta + 1e-4), tracksMixed->GetYaxis()->FindBin(binWidthEta - 1e-4), mixedNormError); - int nBinsMixedNorm = 4; // NOTE this is fixed on purpose, even if binning is made finer around (0,0), this corresponds to the equivalent of four "large" bins around (0,0) + Int_t nBinsMixedNorm = 4; // NOTE this is fixed on purpose, even if binning is made finer around (0,0), this corresponds to the equivalent of four "large" bins around (0,0) mixedNorm /= nBinsMixedNorm; mixedNormError /= nBinsMixedNorm; @@ -618,7 +615,7 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati // finite bin correction if (mTrackEtaCut > 0) { - double finiteBinCorrection = -1.0 / (2 * mTrackEtaCut) * binWidthEta / 2 + 1; + Double_t finiteBinCorrection = -1.0 / (2 * mTrackEtaCut) * binWidthEta / 2 + 1; LOGF(info, "Finite bin correction: %f", finiteBinCorrection); mixedNorm /= finiteBinCorrection; mixedNormError /= finiteBinCorrection; @@ -626,7 +623,7 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati LOGF(error, "ERROR: mTrackEtaCut not set. Finite bin correction cannot be applied. Continuing anyway..."); } - float triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin); + Float_t triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin); // Printf("%f +- %f | %f | %f", mixedNorm, mixedNormError, triggers, mixedNorm / triggers); if (triggers <= 0) { LOGF(error, "ERROR: Skipping multiplicity %d because mixed event is empty", multBin); @@ -634,6 +631,7 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati } mixedNorm /= triggers; + mixedNormError /= triggers; delete tracksMixed; } else { @@ -650,15 +648,15 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati // normParameters->Fill(mixedNorm); TAxis* vertexAxis = trackSameAll->GetAxis(2); - int vertexBinBegin = 1; - int vertexBinEnd = vertexAxis->GetNbins(); + Int_t vertexBinBegin = 1; + Int_t vertexBinEnd = vertexAxis->GetNbins(); if (mZVtxMax > mZVtxMin) { vertexBinBegin = vertexAxis->FindBin(mZVtxMin); vertexBinEnd = vertexAxis->FindBin(mZVtxMax); } - for (int vertexBin = vertexBinBegin; vertexBin <= vertexBinEnd; vertexBin++) { + for (Int_t vertexBin = vertexBinBegin; vertexBin <= vertexBinEnd; vertexBin++) { trackSameAll->GetAxis(2)->SetRange(vertexBin, vertexBin); trackMixedAll->GetAxis(2)->SetRange(vertexBin, vertexBin); @@ -669,13 +667,13 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati // TH1* histMixedproj = mixedTwoD->ProjectionY(); // histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX()); // - // for (int x=1; x<=mixedTwoD->GetNbinsX(); x++) - // for (int y=1; y<=mixedTwoD->GetNbinsY(); y++) + // for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++) + // for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++) // mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y)); // delete histMixedproj; - float triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin); + Float_t triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin); if (triggers2 <= 0) { LOGF(error, "ERROR: Skipping multiplicity %d vertex bin %d because mixed event is empty", multBin, vertexBin); } else { @@ -690,14 +688,14 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati // new TCanvas; tracksMixed->DrawClone("SURF1"); // some code to judge the relative contribution of the different correlation functions to the overall uncertainty - double sums[] = {0, 0, 0}; - double errors[] = {0, 0, 0}; + Double_t sums[] = {0, 0, 0}; + Double_t errors[] = {0, 0, 0}; - int checkBinYBegin = 1; // tracksSame->GetXaxis()->FindBin(-0.79); - int checkBinYEnd = tracksSame->GetNbinsY(); // tracksSame->GetXaxis()->FindBin(0.79); + Int_t checkBinYBegin = 1; // tracksSame->GetXaxis()->FindBin(-0.79); + Int_t checkBinYEnd = tracksSame->GetNbinsY(); // tracksSame->GetXaxis()->FindBin(0.79); - for (int x = 1; x <= tracksSame->GetNbinsX(); x++) { - for (int y = checkBinYBegin; y <= checkBinYEnd; y++) { + for (Int_t x = 1; x <= tracksSame->GetNbinsX(); x++) { + for (Int_t y = checkBinYBegin; y <= checkBinYEnd; y++) { sums[0] += tracksSame->GetBinContent(x, y); errors[0] += tracksSame->GetBinError(x, y); sums[1] += tracksMixed->GetBinContent(x, y); @@ -707,14 +705,14 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati tracksSame->Divide(tracksMixed); - for (int x = 1; x <= tracksSame->GetNbinsX(); x++) { - for (int y = checkBinYBegin; y <= checkBinYEnd; y++) { + for (Int_t x = 1; x <= tracksSame->GetNbinsX(); x++) { + for (Int_t y = checkBinYBegin; y <= checkBinYEnd; y++) { sums[2] += tracksSame->GetBinContent(x, y); errors[2] += tracksSame->GetBinError(x, y); } } - for (int x = 0; x < 3; x++) { + for (Int_t x = 0; x < 3; x++) { if (sums[x] > 0) { errors[x] /= sums[x]; } @@ -730,7 +728,7 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati */ if (!totalTracks) { - totalTracks = dynamic_cast(tracksSame->Clone("totalTracks")); + totalTracks = (TH2*)tracksSame->Clone("totalTracks"); } else { totalTracks->Add(tracksSame); } @@ -748,11 +746,11 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati } if (totalTracks) { - double sums[] = {0, 0, 0}; - double errors[] = {0, 0, 0}; + Double_t sums[] = {0, 0, 0}; + Double_t errors[] = {0, 0, 0}; - for (int x = 1; x <= totalTracks->GetNbinsX(); x++) { - for (int y = 1; y <= totalTracks->GetNbinsY(); y++) { + for (Int_t x = 1; x <= totalTracks->GetNbinsX(); x++) { + for (Int_t y = 1; y <= totalTracks->GetNbinsY(); y++) { sums[0] += totalTracks->GetBinContent(x, y); errors[0] += totalTracks->GetBinError(x, y); } @@ -762,17 +760,17 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati } if (normalizePerTrigger) { - LOGF(info, "Dividing %f tracks by %ld events (%d correlation function(s)) (error %f)", totalTracks->Integral(), totalEvents, nCorrelationFunctions, errors[0]); + LOGF(info, "Dividing %f tracks by %lld events (%d correlation function(s)) (error %f)", totalTracks->Integral(), totalEvents, nCorrelationFunctions, errors[0]); if (totalEvents > 0) { totalTracks->Scale(1.0 / totalEvents); } } if (trigger != nullptr) { - *trigger = static_cast(totalEvents); + *trigger = (Int_t)totalEvents; } // normalizate to dphi width - float normalization = totalTracks->GetXaxis()->GetBinWidth(1); + Float_t normalization = totalTracks->GetXaxis()->GetBinWidth(1); totalTracks->Scale(1.0 / normalization); } @@ -790,14 +788,14 @@ TH2* CorrelationContainer::getSumOfRatios(CorrelationContainer* mixed, Correlati return totalTracks; } -TH1* CorrelationContainer::getTriggersAsFunctionOfMultiplicity(CorrelationContainer::CFStep step, float ptTriggerMin, float ptTriggerMax) +TH1* CorrelationContainer::getTriggersAsFunctionOfMultiplicity(CorrelationContainer::CFStep step, Float_t ptTriggerMin, Float_t ptTriggerMax) { // returns the distribution of triggers as function of centrality/multiplicity resetBinLimits(mTriggerHist->getTHn(step), 3); - int firstBin = mTriggerHist->getTHn(step)->GetAxis(0)->FindBin(ptTriggerMin); - int lastBin = mTriggerHist->getTHn(step)->GetAxis(0)->FindBin(ptTriggerMax); + Int_t firstBin = mTriggerHist->getTHn(step)->GetAxis(0)->FindBin(ptTriggerMin); + Int_t lastBin = mTriggerHist->getTHn(step)->GetAxis(0)->FindBin(ptTriggerMax); LOGF(info, "Using pT range %d --> %d", firstBin, lastBin); mTriggerHist->getTHn(step)->GetAxis(0)->SetRange(firstBin, lastBin); @@ -829,7 +827,7 @@ THnBase* CorrelationContainer::getTrackEfficiencyND(EfficiencyStep step1, Effici sourceContainer->getTHn(step2)->GetAxis(0)->SetRangeUser(mEtaMin, mEtaMax); } - int dimensions[] = {0, 1, 3, 4}; + Int_t dimensions[] = {0, 1, 3, 4}; THnBase* generated = sourceContainer->getTHn(step1)->ProjectionND(4, dimensions); THnBase* measured = sourceContainer->getTHn(step2)->ProjectionND(4, dimensions); @@ -838,7 +836,7 @@ THnBase* CorrelationContainer::getTrackEfficiencyND(EfficiencyStep step1, Effici resetBinLimits(sourceContainer->getTHn(step1)); resetBinLimits(sourceContainer->getTHn(step2)); - auto* clone = dynamic_cast(measured->Clone()); + THnBase* clone = (THnBase*)measured->Clone(); clone->Divide(measured, generated, 1, 1, "B"); @@ -849,7 +847,7 @@ THnBase* CorrelationContainer::getTrackEfficiencyND(EfficiencyStep step1, Effici } //____________________________________________________________________ -TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencyStep step2, int axis1, int axis2, int source, int axis3) +TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencyStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3) { // creates a track-level efficiency by dividing step2 by step1 // projected to axis1 and axis2 (optional if >= 0) @@ -910,8 +908,8 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt } // check for bins with less than 50 entries, print warning - int binBegin[3]; - int binEnd[3]; + Int_t binBegin[3]; + Int_t binEnd[3]; binBegin[0] = 1; binBegin[1] = 1; @@ -938,7 +936,7 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt if (mPtMax > mPtMin) { // TODO this is just checking up to 15 for now - float ptMax = std::min(static_cast(15.), mPtMax); + Float_t ptMax = TMath::Min((Float_t)15., mPtMax); if (axis1 == 1) { binBegin[0] = generated->GetXaxis()->FindBin(mPtMin); binEnd[0] = generated->GetXaxis()->FindBin(ptMax); @@ -953,15 +951,15 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt } } - int total = 0; - int count = 0; - int vars[3]; + Int_t total = 0; + Int_t count = 0; + Int_t vars[3]; - for (int i = 0; i < 3; i++) { + for (Int_t i = 0; i < 3; i++) { vars[i] = binBegin[i]; } - const int limit = 50; + const Int_t limit = 50; while (1) { if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit) { LOGF(info, "Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0])); @@ -1009,25 +1007,23 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt // 2d rebin with variable axis does not exist in root TH1* tmp = measured; - auto* measuredRebinned = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray()); - for (int x = 1; x <= tmp->GetNbinsX(); x++) { - for (int y = 1; y <= tmp->GetNbinsY(); y++) { - measuredRebinned->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y)); - measuredRebinned->SetBinError(x, y, 0); // cannot trust bin error, set to 0 + measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray()); + for (Int_t x = 1; x <= tmp->GetNbinsX(); x++) { + for (Int_t y = 1; y <= tmp->GetNbinsY(); y++) { + ((TH2*)measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y)); + measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0 } } - measured = measuredRebinned; delete tmp; tmp = generated; - auto* generatedRebinned = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray()); - for (int x = 1; x <= tmp->GetNbinsX(); x++) { - for (int y = 1; y <= tmp->GetNbinsY(); y++) { - generatedRebinned->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y)); - generatedRebinned->SetBinError(x, y, 0); // cannot trust bin error, set to 0 + generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray()); + for (Int_t x = 1; x <= tmp->GetNbinsX(); x++) { + for (Int_t y = 1; y <= tmp->GetNbinsY(); y++) { + ((TH2*)generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y)); + generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0 } } - generated = generatedRebinned; delete tmp; } else { TH1* tmp = measured; @@ -1051,7 +1047,7 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt LOGF(fatal, "Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()); } - double newBins[] = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0}; + Double_t newBins[] = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0}; // reduce binning below 5 GeV/c TH1* tmp = measured; @@ -1061,7 +1057,7 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt // increase binning above 5 GeV/c tmp = measured; measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins())); - for (int bin = 1; bin <= measured->GetNbinsX(); bin++) { + for (Int_t bin = 1; bin <= measured->GetNbinsX(); bin++) { measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin)))); measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin)))); } @@ -1075,7 +1071,7 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt // increase binning above 5 GeV/c tmp = generated; generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins())); - for (int bin = 1; bin <= generated->GetNbinsX(); bin++) { + for (Int_t bin = 1; bin <= generated->GetNbinsX(); bin++) { generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin)))); generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin)))); } @@ -1094,7 +1090,7 @@ TH1* CorrelationContainer::getTrackEfficiency(EfficiencyStep step1, EfficiencySt } //____________________________________________________________________ -TH1* CorrelationContainer::getEventEfficiency(EfficiencyStep step1, EfficiencyStep step2, int axis1, int axis2, float ptTriggerMin, float ptTriggerMax) +TH1* CorrelationContainer::getEventEfficiency(EfficiencyStep step1, EfficiencyStep step2, Int_t axis1, Int_t axis2, Float_t ptTriggerMin, Float_t ptTriggerMax) { // creates a event-level efficiency by dividing step2 by step1 // projected to axis1 and axis2 (optional if >= 0) @@ -1137,9 +1133,9 @@ void CorrelationContainer::weightHistogram(TH3* hist1, TH1* hist2) LOGF(fatal, "Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()); } - for (int x = 1; x <= hist1->GetNbinsX(); x++) { - for (int y = 1; y <= hist1->GetNbinsY(); y++) { - for (int z = 1; z <= hist1->GetNbinsZ(); z++) { + for (Int_t x = 1; x <= hist1->GetNbinsX(); x++) { + for (Int_t y = 1; y <= hist1->GetNbinsY(); y++) { + for (Int_t z = 1; z <= hist1->GetNbinsZ(); z++) { if (hist2->GetBinContent(z) > 0) { hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z)); hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z)); @@ -1153,7 +1149,7 @@ void CorrelationContainer::weightHistogram(TH3* hist1, TH1* hist2) } //____________________________________________________________________ -TH1* CorrelationContainer::getBias(EfficiencyStep step1, EfficiencyStep step2, const char* axis, float leadPtMin, float leadPtMax, int weighting) +TH1* CorrelationContainer::getBias(EfficiencyStep step1, EfficiencyStep step2, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting) { // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1) // done by weighting the track-level distribution with the number of events as function of leading pT @@ -1174,14 +1170,14 @@ TH1* CorrelationContainer::getBias(EfficiencyStep step1, EfficiencyStep step2, c resetBinLimits(mTriggerHist->getTHn(step2), 3); setBinLimits(tmp->getTHn(step2)); - TH1D* events1 = mTriggerHist->getTHn(step1)->Projection(0); - TH3D* hist1 = tmp->getTHn(step1)->Projection(0, tmp->getNVar() - 1, 2); + TH1D* events1 = (TH1D*)mTriggerHist->getTHn(step1)->Projection(0); + TH3D* hist1 = (TH3D*)tmp->getTHn(step1)->Projection(0, tmp->getNVar() - 1, 2); if (weighting == 0) { weightHistogram(hist1, events1); } - TH1D* events2 = mTriggerHist->getTHn(step2)->Projection(0); - TH3D* hist2 = tmp->getTHn(step2)->Projection(0, tmp->getNVar() - 1, 2); + TH1D* events2 = (TH1D*)mTriggerHist->getTHn(step2)->Projection(0); + TH3D* hist2 = (TH3D*)tmp->getTHn(step2)->Projection(0, tmp->getNVar() - 1, 2); if (weighting == 0) { weightHistogram(hist2, events2); } @@ -1257,7 +1253,7 @@ TH2D* CorrelationContainer::getTrackingEfficiencyCentrality() } //____________________________________________________________________ -TH1D* CorrelationContainer::getTrackingEfficiency(int axis) +TH1D* CorrelationContainer::getTrackingEfficiency(Int_t axis) { // extracts the tracking efficiency for primary particles // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics @@ -1266,7 +1262,7 @@ TH1D* CorrelationContainer::getTrackingEfficiency(int axis) } //____________________________________________________________________ -TH1D* CorrelationContainer::getFakeRate(int axis) +TH1D* CorrelationContainer::getFakeRate(Int_t axis) { return dynamic_cast(getTrackEfficiency(EfficiencyStep::RecoAll, EfficiencyStep::Fake, axis)); } @@ -1282,7 +1278,7 @@ TH2D* CorrelationContainer::getTrackingCorrection() } //____________________________________________________________________ -TH1D* CorrelationContainer::getTrackingCorrection(int axis) +TH1D* CorrelationContainer::getTrackingCorrection(Int_t axis) { // extracts the correction for tracking efficiency and secondaries // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics @@ -1313,7 +1309,7 @@ TH2D* CorrelationContainer::getTrackingEfficiencyCorrectionCentrality() } //____________________________________________________________________ -TH1D* CorrelationContainer::getTrackingEfficiencyCorrection(int axis) +TH1D* CorrelationContainer::getTrackingEfficiencyCorrection(Int_t axis) { // extracts the correction for tracking efficiency // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics @@ -1344,7 +1340,7 @@ TH2D* CorrelationContainer::getTrackingContaminationCentrality() } //____________________________________________________________________ -TH1D* CorrelationContainer::getTrackingContamination(int axis) +TH1D* CorrelationContainer::getTrackingContamination(Int_t axis) { // extracts the tracking contamination by secondaries // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics @@ -1390,7 +1386,7 @@ void CorrelationContainer::deepCopy(CorrelationContainer* from) // copies the entries of this object's members from the object to this object // fills using the fill function and thus allows that the objects have different binning - for (int step = 0; step < mPairHist->getNSteps(); step++) { + for (Int_t step = 0; step < mPairHist->getNSteps(); step++) { LOGF(info, "Copying step %d", step); THnBase* target = mPairHist->getTHn(step); THnBase* source = from->mPairHist->getTHn(step); @@ -1399,7 +1395,7 @@ void CorrelationContainer::deepCopy(CorrelationContainer* from) target->RebinnedAdd(source); } - for (int step = 0; step < mTriggerHist->getNSteps(); step++) { + for (Int_t step = 0; step < mTriggerHist->getNSteps(); step++) { LOGF(info, "Ev: Copying step %d", step); THnBase* target = mTriggerHist->getTHn(step); THnBase* source = from->mTriggerHist->getTHn(step); @@ -1408,7 +1404,7 @@ void CorrelationContainer::deepCopy(CorrelationContainer* from) target->RebinnedAdd(source); } - for (int step = 0; step < std::min(mTrackHistEfficiency->getNSteps(), from->mTrackHistEfficiency->getNSteps()); step++) { + for (Int_t step = 0; step < TMath::Min(mTrackHistEfficiency->getNSteps(), from->mTrackHistEfficiency->getNSteps()); step++) { if (!from->mTrackHistEfficiency->getTHn(step)) { continue; } @@ -1426,7 +1422,7 @@ void CorrelationContainer::symmetrizepTBins() { // copy pt,a < pt,t bins to pt,a > pt,t (inverting deltaphi and delta eta as it should be) including symmetric bins - for (int step = 0; step < mPairHist->getNSteps(); step++) { + for (Int_t step = 0; step < mPairHist->getNSteps(); step++) { LOGF(info, "Copying step %d", step); THnBase* target = mPairHist->getTHn(step); if (target->GetEntries() == 0) { @@ -1434,47 +1430,47 @@ void CorrelationContainer::symmetrizepTBins() } // for symmetric bins - auto* source = dynamic_cast(target->Clone()); + THnBase* source = (THnBase*)target->Clone(); // axes: 0 delta eta; 1 pT,a; 2 pT,t; 3 centrality; 4 delta phi; 5 vtx-z - for (int i3 = 1; i3 <= target->GetAxis(3)->GetNbins(); i3++) { - for (int i5 = 1; i5 <= target->GetAxis(5)->GetNbins(); i5++) { - for (int i1 = 1; i1 <= target->GetAxis(1)->GetNbins(); i1++) { - for (int i2 = 1; i2 <= target->GetAxis(2)->GetNbins(); i2++) { + for (Int_t i3 = 1; i3 <= target->GetAxis(3)->GetNbins(); i3++) { + for (Int_t i5 = 1; i5 <= target->GetAxis(5)->GetNbins(); i5++) { + for (Int_t i1 = 1; i1 <= target->GetAxis(1)->GetNbins(); i1++) { + for (Int_t i2 = 1; i2 <= target->GetAxis(2)->GetNbins(); i2++) { // find source bin - int binA = target->GetAxis(1)->FindBin(target->GetAxis(2)->GetBinCenter(i2)); - int binT = target->GetAxis(2)->FindBin(target->GetAxis(1)->GetBinCenter(i1)); + Int_t binA = target->GetAxis(1)->FindBin(target->GetAxis(2)->GetBinCenter(i2)); + Int_t binT = target->GetAxis(2)->FindBin(target->GetAxis(1)->GetBinCenter(i1)); LOGF(info, "(%d %d) Copying from %d %d to %d %d", i3, i5, binA, binT, i1, i2); - for (int i0 = 1; i0 <= target->GetAxis(0)->GetNbins(); i0++) { - for (int i4 = 1; i4 <= target->GetAxis(4)->GetNbins(); i4++) { - int binEta = target->GetAxis(0)->FindBin(-target->GetAxis(0)->GetBinCenter(i0)); - double phi = -target->GetAxis(4)->GetBinCenter(i4); + for (Int_t i0 = 1; i0 <= target->GetAxis(0)->GetNbins(); i0++) { + for (Int_t i4 = 1; i4 <= target->GetAxis(4)->GetNbins(); i4++) { + Int_t binEta = target->GetAxis(0)->FindBin(-target->GetAxis(0)->GetBinCenter(i0)); + Double_t phi = -target->GetAxis(4)->GetBinCenter(i4); if (phi < -PIHalf) { phi += TwoPI; } - int binPhi = target->GetAxis(4)->FindBin(phi); + Int_t binPhi = target->GetAxis(4)->FindBin(phi); - int binSource[] = {binEta, binA, binT, i3, binPhi, i5}; - int binTarget[] = {i0, i1, i2, i3, i4, i5}; + Int_t binSource[] = {binEta, binA, binT, i3, binPhi, i5}; + Int_t binTarget[] = {i0, i1, i2, i3, i4, i5}; - double value = source->GetBinContent(binSource); - double error = source->GetBinError(binSource); + Double_t value = source->GetBinContent(binSource); + Double_t error = source->GetBinError(binSource); if (error == 0) { continue; } - double value2 = target->GetBinContent(binTarget); - double error2 = target->GetBinError(binTarget); + Double_t value2 = target->GetBinContent(binTarget); + Double_t error2 = target->GetBinError(binTarget); - double sum = value; - double err = error; + Double_t sum = value; + Double_t err = error; if (error2 > 0) { sum = value + value2; - err = std::sqrt(error * error + error2 * error2); + err = TMath::Sqrt(error * error + error2 * error2); } // Printf(" Values: %f +- %f; %f +- %f --> %f +- %f", value, error, value2, error2, sum, err); @@ -1493,13 +1489,13 @@ void CorrelationContainer::symmetrizepTBins() } //____________________________________________________________________ -void CorrelationContainer::extendTrackingEfficiency(bool verbose) +void CorrelationContainer::extendTrackingEfficiency(Bool_t verbose) { // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency - float fitRangeBegin = 5.01; - float fitRangeEnd = 14.99; - float extendRangeBegin = 10.01; + Float_t fitRangeBegin = 5.01; + Float_t fitRangeEnd = 14.99; + Float_t extendRangeBegin = 10.01; if (mTrackHistEfficiency->getNVar() == 3) { TH1* obj = getTrackingEfficiency(1); @@ -1511,7 +1507,7 @@ void CorrelationContainer::extendTrackingEfficiency(bool verbose) obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd); - float trackingEff = obj->GetFunction("pol0")->GetParameter(0); + Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0); obj = getTrackingContamination(1); @@ -1522,17 +1518,17 @@ void CorrelationContainer::extendTrackingEfficiency(bool verbose) obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd); - float trackingCont = obj->GetFunction("pol0")->GetParameter(0); + Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0); LOGF(info, "CorrelationContainer::extendTrackingEfficiency: Fitted efficiency between %f and %f and got %f tracking efficiency and %f tracking contamination correction. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, trackingEff, trackingCont, extendRangeBegin, mEtaMin, mEtaMax); // extend for full pT range - for (int x = mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMin); x <= mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMax); x++) { - for (int y = mTrackHistEfficiency->getTHn(0)->GetAxis(1)->FindBin(extendRangeBegin); y <= mTrackHistEfficiency->getTHn(0)->GetAxis(1)->GetNbins(); y++) { - for (int z = 1; z <= mTrackHistEfficiency->getTHn(0)->GetAxis(2)->GetNbins(); z++) // particle type axis + for (Int_t x = mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMin); x <= mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMax); x++) { + for (Int_t y = mTrackHistEfficiency->getTHn(0)->GetAxis(1)->FindBin(extendRangeBegin); y <= mTrackHistEfficiency->getTHn(0)->GetAxis(1)->GetNbins(); y++) { + for (Int_t z = 1; z <= mTrackHistEfficiency->getTHn(0)->GetAxis(2)->GetNbins(); z++) // particle type axis { - int bins[3]; + Int_t bins[3]; bins[0] = x; bins[1] = y; bins[2] = z; @@ -1545,28 +1541,28 @@ void CorrelationContainer::extendTrackingEfficiency(bool verbose) } } else if (mTrackHistEfficiency->getNVar() == 4) { // fit in centrality intervals of 20% for efficiency, one bin for contamination - float* trackingEff = nullptr; - float* trackingCont = nullptr; - float centralityBins[] = {0, 10, 20, 40, 60, 100}; - int nCentralityBins = 5; + Float_t* trackingEff = nullptr; + Float_t* trackingCont = nullptr; + Float_t centralityBins[] = {0, 10, 20, 40, 60, 100}; + Int_t nCentralityBins = 5; LOGF(info, "CorrelationContainer::extendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, mEtaMin, mEtaMax); // 0 = eff; 1 = cont - for (int caseNo = 0; caseNo < 2; caseNo++) { - float* target = nullptr; - int centralityBinsLocal = nCentralityBins; + for (Int_t caseNo = 0; caseNo < 2; caseNo++) { + Float_t* target = nullptr; + Int_t centralityBinsLocal = nCentralityBins; if (caseNo == 0) { - trackingEff = new float[centralityBinsLocal]; + trackingEff = new Float_t[centralityBinsLocal]; target = trackingEff; } else { centralityBinsLocal = 1; - trackingCont = new float[centralityBinsLocal]; + trackingCont = new Float_t[centralityBinsLocal]; target = trackingCont; } - for (int i = 0; i < centralityBinsLocal; i++) { + for (Int_t i = 0; i < centralityBinsLocal; i++) { if (centralityBinsLocal == 1) { setCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1); } else { @@ -1577,7 +1573,7 @@ void CorrelationContainer::extendTrackingEfficiency(bool verbose) new TCanvas; proj->DrawCopy(); } - if (static_cast(proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd)) == 0) { + if ((Int_t)proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0) { target[i] = proj->GetFunction("pol0")->GetParameter(0); } else { target[i] = 0; @@ -1587,20 +1583,20 @@ void CorrelationContainer::extendTrackingEfficiency(bool verbose) } // extend for full pT range - for (int x = mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMin); x <= mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMax); x++) { - for (int y = mTrackHistEfficiency->getTHn(0)->GetAxis(1)->FindBin(extendRangeBegin); y <= mTrackHistEfficiency->getTHn(0)->GetAxis(1)->GetNbins(); y++) { - for (int z = 1; z <= mTrackHistEfficiency->getTHn(0)->GetAxis(2)->GetNbins(); z++) // particle type axis + for (Int_t x = mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMin); x <= mTrackHistEfficiency->getTHn(0)->GetAxis(0)->FindBin(mEtaMax); x++) { + for (Int_t y = mTrackHistEfficiency->getTHn(0)->GetAxis(1)->FindBin(extendRangeBegin); y <= mTrackHistEfficiency->getTHn(0)->GetAxis(1)->GetNbins(); y++) { + for (Int_t z = 1; z <= mTrackHistEfficiency->getTHn(0)->GetAxis(2)->GetNbins(); z++) // particle type axis { - for (int z2 = 1; z2 <= mTrackHistEfficiency->getTHn(0)->GetAxis(3)->GetNbins(); z2++) // centrality axis + for (Int_t z2 = 1; z2 <= mTrackHistEfficiency->getTHn(0)->GetAxis(3)->GetNbins(); z2++) // centrality axis { - int bins[4]; + Int_t bins[4]; bins[0] = x; bins[1] = y; bins[2] = z; bins[3] = z2; - int z2Bin = 0; + Int_t z2Bin = 0; while (centralityBins[z2Bin + 1] < mTrackHistEfficiency->getTHn(0)->GetAxis(3)->GetBinCenter(z2)) { z2Bin++; } @@ -1627,11 +1623,11 @@ void CorrelationContainer::extendTrackingEfficiency(bool verbose) } /* -void CorrelationContainer::Scale(double factor) +void CorrelationContainer::Scale(Double_t factor) { // scales all contained histograms by the given factor - for (int i=0; i<4; i++) + for (Int_t i=0; i<4; i++) if (mTrackHist[i]) mTrackHist[i]->Scale(factor); @@ -1643,15 +1639,15 @@ void CorrelationContainer::Reset() { // resets all contained histograms - for (int step = 0; step < mPairHist->getNSteps(); step++) { + for (Int_t step = 0; step < mPairHist->getNSteps(); step++) { mPairHist->getTHn(step)->Reset(); } - for (int step = 0; step < mTriggerHist->getNSteps(); step++) { + for (Int_t step = 0; step < mTriggerHist->getNSteps(); step++) { mTriggerHist->getTHn(step)->Reset(); } - for (int step = 0; step < mTrackHistEfficiency->getNSteps(); step++) { + for (Int_t step = 0; step < mTrackHistEfficiency->getNSteps(); step++) { mTrackHistEfficiency->getTHn(step)->Reset(); } } @@ -1663,7 +1659,7 @@ THnBase* CorrelationContainer::changeToThn(THnBase* sparse) return THn::CreateHn(Form("%s_thn", sparse->GetName()), sparse->GetTitle(), sparse); } -void CorrelationContainer::fillEvent(float centrality, CFStep step) +void CorrelationContainer::fillEvent(Float_t centrality, CFStep step) { // Fill per-event information mEventCount->Fill(step, centrality); diff --git a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt index 83f52b206e0..8c283cecc40 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt +++ b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt @@ -78,11 +78,6 @@ o2physics_add_dpl_workflow(long-range-dihadron-cor PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore O2Physics::AnalysisCCDB O2Physics::GFWCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(mft-correlations - SOURCES mftcorrelations.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore - COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(two-particle-correlation-pp SOURCES twoParticleCorrelationPp.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGCFCore diff --git a/PWGCF/TwoParticleCorrelations/Tasks/mftcorrelations.cxx b/PWGCF/TwoParticleCorrelations/Tasks/mftcorrelations.cxx deleted file mode 100644 index e34bd64f3a9..00000000000 --- a/PWGCF/TwoParticleCorrelations/Tasks/mftcorrelations.cxx +++ /dev/null @@ -1,1449 +0,0 @@ -// 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 mftcorrelations.cxx -/// \brief task for two-particle (hadron-hadron) angular correlations built from MFT tracks - -/// \author Mintu Haldar - -#include "PWGCF/Core/CorrelationContainer.h" - -#include "Common/Core/RecoDecay.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include // Pair<>, SameKindPair<> -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace constants::math; - -// Bitmask describing which multiplicity/centrality estimators feed the "multCorrelations" QA -// histogram and, when set, the outlier-rejection formula. -namespace o2::aod -{ -namespace mftmultset -{ -DECLARE_SOA_COLUMN(Multiplicities, multiplicities, std::vector); //! List of auxiliary multiplicities -enum MultiplicityEstimators : uint16_t { - CentFT0C = 1 << 0, - MultFV0A = 1 << 1, - MultNTracksPV = 1 << 2, - MultNTracksGlobal = 1 << 3, - CentFT0M = 1 << 4 -}; -inline constexpr uint32_t NMultiplicityEstimators = __builtin_ctz(CentFT0M) + 1; -} // namespace mftmultset - -DECLARE_SOA_TABLE(MFTMultSets, "AOD", "MFTMULTSET", mftmultset::Multiplicities); //! Auxilary multiplicity set table -using MFTMultSet = MFTMultSets::iterator; -} // namespace o2::aod - -#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; - -struct Mftcorrelations { - SliceCache cache; - - // ------------------------------------------------------------------ Configurables - O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range") - - // Global acceptance - O2_DEFINE_CONFIGURABLE(cfgGlobalTrackEtaCut, float, 0.8f, "Eta cut for global (central-barrel) tracks used as MFT-Global trigger") - O2_DEFINE_CONFIGURABLE(cfgGlobalTrackPtMin, float, 0.15f, "Minimal pT for global tracks") - O2_DEFINE_CONFIGURABLE(cfgGlobalTrackPtMax, float, 10.0f, "Maximal pT for global tracks") - - // MFT acceptance is roughly -3.6 < eta < -2.5; keep as configurables since it is tag/setup dependent. - O2_DEFINE_CONFIGURABLE(cfgMFTCutEtaMin, float, -3.6f, "Minimal eta for MFT tracks") - O2_DEFINE_CONFIGURABLE(cfgMFTCutEtaMax, float, -2.5f, "Maximal eta for MFT tracks") - O2_DEFINE_CONFIGURABLE(cfgMFTCutPtMin, float, 0.5f, "Minimal pT for MFT tracks") - O2_DEFINE_CONFIGURABLE(cfgMFTCutPtMax, float, 20.0f, "Maximal pT for MFT tracks") - O2_DEFINE_CONFIGURABLE(cfgMFTMinNClusters, int, 5, "Minimal number of MFT clusters on track") - - O2_DEFINE_CONFIGURABLE(cfgPtOrder, int, 1, "Only consider pairs for which pT,1 < pT,2 (0 = OFF, 1 = ON)") - O2_DEFINE_CONFIGURABLE(cfgTriggerCharge, int, 0, "Select on charge of trigger particle: 0 = all; 1 = positive; -1 = negative") - O2_DEFINE_CONFIGURABLE(cfgAssociatedCharge, int, 0, "Select on charge of associated particle: 0 = all; 1 = positive; -1 = negative") - O2_DEFINE_CONFIGURABLE(cfgPairCharge, int, 0, "Select on charge of the pair: 0 = all; 1 = like sign; -1 = unlike sign") - - O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, int, 0, "0 = OFF and 1 = ON: read efficiency THn from a local ROOT file instead of CCDB") - O2_DEFINE_CONFIGURABLE(cfgDropStepRECO, bool, false, "drop step RECO if efficiency correction is used") - O2_DEFINE_CONFIGURABLE(cfgEfficiencyTrigger, std::string, "", "CCDB path (or local file if cfgLocalEfficiency=1) to efficiency object for trigger particles") - O2_DEFINE_CONFIGURABLE(cfgEfficiencyAssociated, std::string, "", "CCDB path (or local file if cfgLocalEfficiency=1) to efficiency object for associated particles") - // FIX #6d: efficiency for central-barrel (global) tracks, used in MFT-Global (trigger) and Global-Global channels. - O2_DEFINE_CONFIGURABLE(cfgEfficiencyGlobal, std::string, "", "CCDB path (or local file if cfgLocalEfficiency=1) to efficiency object for central-barrel (global) tracks") - - O2_DEFINE_CONFIGURABLE(cfgMultCorrelationsMask, uint16_t, 0, "Bitmask of multiplicity/centrality estimators to fill in multCorrelations QA / use for outlier rejection (see mftmultset enum)") - O2_DEFINE_CONFIGURABLE(cfgMultCutFormula, std::string, "", "Multiplicity outlier-rejection formula. A result > 0 accepts the event. Parameters: [cFT0C], [mFV0A], [mPV], [mGlob], [cFT0M]") - - O2_DEFINE_CONFIGURABLE(cfgNoMixedEvents, int, 5, "Number of mixed events per event") - O2_DEFINE_CONFIGURABLE(cfgVerbosity, int, 1, "Verbosity level (0 = major, 1 = per collision)") - - // ------------------------------------------------------------------ Phi weight - // NOTE: keep OFF (empty file) when the efficiency object already has a phi axis, otherwise phi is corrected twice. - O2_DEFINE_CONFIGURABLE(cfgPhiWeightFile, std::string, "", "Local ROOT file with the phi-weight TH1 (empty = disabled)") - O2_DEFINE_CONFIGURABLE(cfgPhiWeightName, std::string, "", "Name of the phi-weight TH1 object inside cfgPhiWeightFile") - - // ------------------------------------------------------------------ Axes - ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"}; - ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"}; - ConfigurableAxis axisDeltaEta{"axisDeltaEta", {40, -2.2, 2.2}, "delta eta axis for histograms (MFT-MFT so somewhat narrower than TPC-TPC)"}; - ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 3.0}, "pt trigger axis for histograms"}; - ConfigurableAxis axisPtAssoc{"axisPtAssoc", {VARIABLE_WIDTH, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0}, "pt associated axis for histograms"}; - ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 100.1}, "multiplicity / centrality axis for histograms"}; - - ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"}; - ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -3.6, -2.5}, "eta axis for efficiency histograms (MFT acceptance)"}; - ConfigurableAxis axisPtEfficiency{"axisPtEfficiency", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.5, 2.0, 3.0}, "pt axis for efficiency histograms"}; - ConfigurableAxis axisPhiEfficiency{"axisPhiEfficiency", {36, 0, TwoPI}, "phi axis for efficiency histograms"}; - - ConfigurableAxis axisEtaEfficiencyGlobal{"axisEtaEfficiencyGlobal", {16, -0.8, 0.8}, "eta axis for efficiency histograms (global tracks)"}; - ConfigurableAxis axisPtEfficiencyGlobal{"axisPtEfficiencyGlobal", {VARIABLE_WIDTH, 0.0, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 20.0}, "pt axis for efficiency histograms (global tracks)"}; - - ConfigurableAxis axisMultCorrCent{"axisMultCorrCent", {100, 0, 100}, "multiplicity correlation axis for centralities"}; - ConfigurableAxis axisMultCorrV0{"axisMultCorrV0", {1000, 0, 100000}, "multiplicity correlation axis for V0 multiplicities"}; - ConfigurableAxis axisMultCorrMult{"axisMultCorrMult", {1000, 0, 1000}, "multiplicity correlation axis for track multiplicities"}; - - // ------------------------------------------------------------------ QA Axes - ConfigurableAxis axisMultQA{"axisMultQA", {100, 0, 100}, "multiplicity/centrality axis for QA histograms"}; - ConfigurableAxis axisPtQA{"axisPtQA", {VARIABLE_WIDTH, 0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0}, "pT axis for QA histograms (yields)"}; - ConfigurableAxis axisEtaQA{"axisEtaQA", {50, -4, -2}, "eta axis for QA histograms"}; - ConfigurableAxis axisPhiQA{"axisPhiQA", {200, 0, TwoPI}, "phi axis for QA histograms"}; - - // Global MFT - ConfigurableAxis axisEtaGlobalQA{"axisEtaGlobalQA", {40, -1, 1}, "eta axis for global-track QA histograms"}; - ConfigurableAxis axisDeltaEtaMftGlobal{"axisDeltaEtaMftGlobal", {60, 1.0, 5.0}, "delta eta axis for MFT-Global correlations (deltaEta = etaGlobal - etaMFT is POSITIVE, roughly 1.7-4.4); override via JSON if needed"}; - // FIX #4: delta-eta axis for the global-global (central-barrel hadron-hadron) channel. - ConfigurableAxis axisDeltaEtaGlobal{"axisDeltaEtaGlobal", {32, -1.6, 1.6}, "delta eta axis for global-global (central-barrel hadron-hadron) correlations"}; - - // ------------------------------------------------------------------ Event selection - Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; - - // ------------------------------------------------------------------ MFT track selection - Filter mftTrackFilter = (aod::fwdtrack::eta > cfgMFTCutEtaMin) && (aod::fwdtrack::eta < cfgMFTCutEtaMax) && - (aod::fwdtrack::pt > cfgMFTCutPtMin) && (aod::fwdtrack::pt < cfgMFTCutPtMax); - - // ------------------------------------------------------------------ Global track selection - Filter globalTrackFilter = (nabs(aod::track::eta) < cfgGlobalTrackEtaCut) && - (aod::track::pt > cfgGlobalTrackPtMin) && (aod::track::pt < cfgGlobalTrackPtMax) && - requireGlobalTrackInFilter(); - - // ------------------------------------------------------------------ Output - OutputObj same{"sameEvent"}; - OutputObj mixed{"mixedEvent"}; - - OutputObj sameMftGlobal{"sameEventMftGlobal"}; - OutputObj mixedMftGlobal{"mixedEventMftGlobal"}; - - // FIX #4: global-global (central-barrel hadron-hadron) output objects. - OutputObj sameGlobal{"sameEventGlobalGlobal"}; - OutputObj mixedGlobal{"mixedEventGlobalGlobal"}; - - HistogramRegistry registry{"registry"}; - - Service ccdb; - Service pdg; - - struct Config { - THn* mEfficiencyTrigger = nullptr; - THn* mEfficiencyAssociated = nullptr; - THn* mEfficiencyGlobal = nullptr; // FIX #6d - TH1* mPhiWeight = nullptr; - bool efficiencyLoaded = false; - } cfg; - - std::vector efficiencyAssociatedCache; - std::vector phiWeightCache; - std::vector phiWeightCache1; - std::unique_ptr multCutFormula; - std::array multCutFormulaParamIndex{}; - - // Collisions are expected to already carry the joined centrality/multiplicity estimator tables. - using ExtendedCollisions = soa::Join; - using FilteredCollisions = soa::Filtered; - using FilteredMFTTracks = soa::Filtered; - using FilteredGlobalTracks = soa::Filtered>; - - // For MFT MC Efficiency - // NOTE SmallGroups includes soa::Filtered always - using McEffCollisions = soa::SmallGroups>; - using McEffMFTTracks = soa::Filtered>; - - // For Global MC Efficiency - using McEffGlobalTracks = soa::Filtered>; - - void init(o2::framework::InitContext&) - { - registry.add("yields", "multiplicity/centrality vs pT vs eta", {HistType::kTH3F, {{axisMultQA, "multiplicity/centrality"}, {axisPtQA, "p_{T}"}, {axisEtaQA, "#eta"}}}); - registry.add("etaphi", "multiplicity/centrality vs eta vs phi", {HistType::kTH3F, {{axisMultQA, "multiplicity/centrality"}, {axisEtaQA, "#eta"}, {axisPhiQA, "#varphi"}}}); - // FIX #6b: efficiency-corrected QA histograms (weight = efficiency correction factor) - registry.add("yields_corrected", "efficiency-corrected multiplicity/centrality vs pT vs eta", {HistType::kTH3F, {{axisMultQA, "multiplicity/centrality"}, {axisPtQA, "p_{T}"}, {axisEtaQA, "#eta"}}}); - registry.add("etaphi_corrected", "efficiency-corrected multiplicity/centrality vs eta vs phi", {HistType::kTH3F, {{axisMultQA, "multiplicity/centrality"}, {axisEtaQA, "#eta"}, {axisPhiQA, "#varphi"}}}); - registry.add("multiplicity", "event multiplicity", {HistType::kTH1F, {{axisMultQA, "multiplicity/centrality"}}}); - registry.add("nClusters", "MFT track number of clusters", {HistType::kTH1F, {{12, -0.5, 11.5, "n clusters"}}}); - registry.add("multiplicity_raw", "event multiplicity (no cuts)", {HistType::kTH1F, {{axisMultQA, "multiplicity/centrality"}}}); - - // MC Efficiency - registry.add("MCEff/hNRecoColl", "N reco collisions per MC collision", {HistType::kTH1F, {{6, -0.5, 5.5}}}); - registry.add("MCEff/hCutFlow", "1=all MC coll, 2=size==1, 3=sel8", {HistType::kTH1F, {{3, 0.5, 3.5}}}); - - registry.add("MCEffGlobal/hNRecoColl", "N reco collisions per MC collision (global)", {HistType::kTH1F, {{6, -0.5, 5.5}}}); - registry.add("MCEffGlobal/hCutFlow", "1=all MC coll, 2=size==1, 3=sel8", {HistType::kTH1F, {{3, 0.5, 3.5}}}); - - if (cfgMultCorrelationsMask != 0) { - std::vector multAxes; - if (cfgMultCorrelationsMask & aod::mftmultset::CentFT0C) - multAxes.emplace_back(axisMultCorrCent, "FT0C centrality"); - if (cfgMultCorrelationsMask & aod::mftmultset::MultFV0A) - multAxes.emplace_back(axisMultCorrV0, "V0A multiplicity"); - if (cfgMultCorrelationsMask & aod::mftmultset::MultNTracksPV) - multAxes.emplace_back(axisMultCorrMult, "Nch PV"); - if (cfgMultCorrelationsMask & aod::mftmultset::MultNTracksGlobal) - multAxes.emplace_back(axisMultCorrMult, "Nch Global"); - if (cfgMultCorrelationsMask & aod::mftmultset::CentFT0M) - multAxes.emplace_back(axisMultCorrCent, "FT0M centrality"); - registry.add("multCorrelations", "Multiplicity correlations", {HistType::kTHnSparseF, multAxes}); - - // FIX #3 - if (!(cfgMultCorrelationsMask & aod::mftmultset::CentFT0C)) { - LOGF(warning, - "cfgMultCorrelationsMask=%d does not include CentFT0C. Mixed-event binning follows " - "getMultiplicity(), so make sure the efficiency objects were produced with the same estimator.", - (int)cfgMultCorrelationsMask); - } - } - - const int maxMixBin = AxisSpec(axisMultiplicity).getNbins() * AxisSpec(axisVertex).getNbins(); - registry.add("eventcount_same", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}}); - registry.add("eventcount_mixed", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}}); - registry.add("trackcount_same", "bin", {HistType::kTH2F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {50, -0.5, 49.5}}}); - registry.add("trackcount_mixed", "bin", {HistType::kTH3F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {50, -0.5, 49.5}, {50, -0.5, 49.5}}}); - - // FIX #1 - registry.add("eventcount_mftglobal_same", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}}); - registry.add("eventcount_mftglobal_mixed", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}}); - registry.add("trackcount_mftglobal_same", "bin", {HistType::kTH3F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {50, -0.5, 49.5}, {50, -0.5, 49.5}}}); - registry.add("trackcount_mftglobal_mixed", "bin", {HistType::kTH3F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {50, -0.5, 49.5}, {50, -0.5, 49.5}}}); - - if (!cfgMultCutFormula.value.empty()) { - multCutFormula = std::make_unique("multCutFormula", cfgMultCutFormula.value.c_str()); - std::fill_n(multCutFormulaParamIndex.begin(), std::size(multCutFormulaParamIndex), ~0u); - std::array pars = {"cFT0C", "mFV0A", "mPV", "mGlob", "cFT0M"}; - for (uint32_t i = 0, n = multCutFormula->GetNpar(); i < n; ++i) { - auto m = std::find(pars.begin(), pars.end(), multCutFormula->GetParName(i)); - if (m == pars.end()) { - LOGF(warning, "Unknown parameter in cfgMultCutFormula: %s", multCutFormula->GetParName(i)); - continue; - } - multCutFormulaParamIndex[std::distance(pars.begin(), m)] = i; - } - } - - registry.add("mftmft/SE/deltaEta_deltaPhi", "same event #Delta#eta-#Delta#varphi (MFT-MFT)", {HistType::kTH2F, {{axisDeltaPhi, "#Delta#varphi"}, {axisDeltaEta, "#Delta#eta"}}}); - registry.add("mftmft/ME/deltaEta_deltaPhi", "mixed event #Delta#eta-#Delta#varphi (MFT-MFT)", {HistType::kTH2F, {{axisDeltaPhi, "#Delta#varphi"}, {axisDeltaEta, "#Delta#eta"}}}); - - registry.add("mftGlobal/SE/Trig_etavsphi", "trigger (global) eta vs phi", {HistType::kTH2F, {{axisPhiQA, "#varphi"}, {axisEtaGlobalQA, "#eta"}}}); - registry.add("mftGlobal/SE/Trig_pt", "trigger (global) pT", {HistType::kTH1F, {{axisPtTrigger, "p_{T}"}}}); - registry.add("mftGlobal/SE/Assoc_etavsphi", "associated (MFT) eta vs phi", {HistType::kTH2F, {{axisPhiQA, "#varphi"}, {axisEtaQA, "#eta"}}}); - registry.add("mftGlobal/SE/Assoc_pt", "associated (MFT) pT", {HistType::kTH1F, {{axisPtAssoc, "p_{T}"}}}); - registry.add("mftGlobal/SE/deltaEta_deltaPhi", "same event #Delta#eta-#Delta#varphi", {HistType::kTH2F, {{axisDeltaPhi, "#Delta#varphi"}, {axisDeltaEtaMftGlobal, "#Delta#eta"}}}); - registry.add("mftGlobal/ME/deltaEta_deltaPhi", "mixed event #Delta#eta-#Delta#varphi", {HistType::kTH2F, {{axisDeltaPhi, "#Delta#varphi"}, {axisDeltaEtaMftGlobal, "#Delta#eta"}}}); - - // FIX #4: QA histograms for the global-global channel. - registry.add("eventcount_global_same", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}}); - registry.add("eventcount_global_mixed", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}}); - registry.add("trackcount_global_same", "bin", {HistType::kTH2F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {50, -0.5, 49.5}}}); - registry.add("trackcount_global_mixed", "bin", {HistType::kTH3F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {50, -0.5, 49.5}, {50, -0.5, 49.5}}}); - - registry.add("global/SE/etavsphi", "global track eta vs phi", {HistType::kTH2F, {{axisPhiQA, "#varphi"}, {axisEtaGlobalQA, "#eta"}}}); - registry.add("global/SE/pt", "global track pT", {HistType::kTH1F, {{axisPtQA, "p_{T}"}}}); - registry.add("global/SE/deltaEta_deltaPhi", "same event #Delta#eta-#Delta#varphi (global-global)", {HistType::kTH2F, {{axisDeltaPhi, "#Delta#varphi"}, {axisDeltaEtaGlobal, "#Delta#eta"}}}); - registry.add("global/ME/deltaEta_deltaPhi", "mixed event #Delta#eta-#Delta#varphi (global-global)", {HistType::kTH2F, {{axisDeltaPhi, "#Delta#varphi"}, {axisDeltaEtaGlobal, "#Delta#eta"}}}); - - std::vector corrAxis = {{axisDeltaEta, "#Delta#eta"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisMultiplicity, "multiplicity / centrality"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisVertex, "z-vtx (cm)"}}; - std::vector effAxis = {{axisEtaEfficiency, "#eta"}, - {axisPtEfficiency, "p_{T} (GeV/c)"}, - {axisVertexEfficiency, "z-vtx (cm)"}, - {axisPhiEfficiency, "#varphi (rad)"}}; - - std::vector effAxisGlobal = {{axisEtaEfficiencyGlobal, "#eta"}, - {axisPtEfficiencyGlobal, "p_{T} (GeV/c)"}, - {axisVertexEfficiency, "z-vtx (cm)"}, - {axisPhiEfficiency, "#varphi (rad)"}}; - - std::vector userAxis; // no extra axis needed for hadron-hadron - - std::vector corrAxisMftGlobal = {{axisDeltaEtaMftGlobal, "#Delta#eta"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisMultiplicity, "multiplicity / centrality"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisVertex, "z-vtx (cm)"}}; - - // FIX #4 - std::vector corrAxisGlobal = {{axisDeltaEtaGlobal, "#Delta#eta"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisMultiplicity, "multiplicity / centrality"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisVertex, "z-vtx (cm)"}}; - - sameMftGlobal.setObject(new CorrelationContainer("sameEventMftGlobal", "sameEventMftGlobal", corrAxisMftGlobal, effAxis, userAxis)); - mixedMftGlobal.setObject(new CorrelationContainer("mixedEventMftGlobal", "mixedEventMftGlobal", corrAxisMftGlobal, effAxis, userAxis)); - - same.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxis, effAxis, userAxis)); - - same->setTrackEtaCut(cfgMFTCutEtaMax); // upper bound of |eta| cut used internally by CorrelationContainer - - sameGlobal.setObject(new CorrelationContainer("sameEventGlobalGlobal", "sameEventGlobalGlobal", corrAxisGlobal, effAxisGlobal, userAxis)); - mixedGlobal.setObject(new CorrelationContainer("mixedEventGlobalGlobal", "mixedEventGlobalGlobal", corrAxisGlobal, effAxisGlobal, userAxis)); - sameGlobal->setTrackEtaCut(cfgGlobalTrackEtaCut); - - if (!cfgEfficiencyAssociated.value.empty()) - efficiencyAssociatedCache.reserve(512); - - // FIX #5: load the phi weight once at init(). - loadPhiWeight(); - - ccdb->setURL("http://alice-ccdb.cern.ch"); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - auto now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); - ccdb->setCreatedNotAfter(now); - } - - // ------------------------------------------------------------------ helpers - - int getMagneticField(uint64_t timestamp) - { - static o2::parameters::GRPObject* grpo = nullptr; - if (grpo == nullptr) { - grpo = ccdb->getForTimeStamp("GLO/GRP/GRP", timestamp); - if (grpo == nullptr) { - LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); - return 0; - } - } - return grpo->getNominalL3Field(); - } - - void loadEfficiency(uint64_t timestamp) - { - if (cfg.efficiencyLoaded) { - return; - } - if (!cfgEfficiencyTrigger.value.empty()) { - if (cfgLocalEfficiency > 0) { - TFile* f = TFile::Open(cfgEfficiencyTrigger.value.c_str(), "READ"); - cfg.mEfficiencyTrigger = dynamic_cast(f->Get("ccdb_object")); - } else { - cfg.mEfficiencyTrigger = ccdb->getForTimeStamp>(cfgEfficiencyTrigger, timestamp); - } - if (cfg.mEfficiencyTrigger == nullptr) { - LOGF(fatal, "Could not load trigger efficiency from %s", cfgEfficiencyTrigger.value.c_str()); - } - LOGF(info, "Loaded MFT trigger efficiency from %s (%p)", cfgEfficiencyTrigger.value.c_str(), static_cast(cfg.mEfficiencyTrigger)); - } - if (!cfgEfficiencyAssociated.value.empty()) { - if (cfgLocalEfficiency > 0) { - TFile* f = TFile::Open(cfgEfficiencyAssociated.value.c_str(), "READ"); - cfg.mEfficiencyAssociated = dynamic_cast(f->Get("ccdb_object")); - } else { - cfg.mEfficiencyAssociated = ccdb->getForTimeStamp>(cfgEfficiencyAssociated, timestamp); - } - if (cfg.mEfficiencyAssociated == nullptr) { - LOGF(fatal, "Could not load associated efficiency from %s", cfgEfficiencyAssociated.value.c_str()); - } - LOGF(info, "Loaded MFT associated efficiency from %s (%p)", cfgEfficiencyAssociated.value.c_str(), static_cast(cfg.mEfficiencyAssociated)); - } - // FIX #6d - if (!cfgEfficiencyGlobal.value.empty()) { - if (cfgLocalEfficiency > 0) { - TFile* f = TFile::Open(cfgEfficiencyGlobal.value.c_str(), "READ"); - cfg.mEfficiencyGlobal = dynamic_cast(f->Get("ccdb_object")); - } else { - cfg.mEfficiencyGlobal = ccdb->getForTimeStamp>(cfgEfficiencyGlobal, timestamp); - } - if (cfg.mEfficiencyGlobal == nullptr) { - LOGF(fatal, "Could not load global-track efficiency from %s", cfgEfficiencyGlobal.value.c_str()); - } - LOGF(info, "Loaded global-track efficiency from %s (%p)", cfgEfficiencyGlobal.value.c_str(), static_cast(cfg.mEfficiencyGlobal)); - } - cfg.efficiencyLoaded = true; - } - - // MFT phi is in [-pi,pi], McParticle/global-track phi in [0,2pi). Use [0,2pi) everywhere. - static float wrapPhi(float phi) - { - return RecoDecay::constrainAngle(phi, 0.f); - } - - // FIX #6a: expected axis order of the efficiency object: (eta, pt, multiplicity, posZ [, phi]). - // - phi is wrapped to [0, 2pi) (MFT track.phi() is in [-pi, pi]) - // - 4-dim and 5-dim objects are supported - // - bins outside the axis range (under/overflow) return 0 = track outside efficiency acceptance - double getEfficiencyCorrection(THn* eff, float eta, float pt, float multiplicity, float posZ, float phi) - { - if (eff == nullptr) { - return 1.0; - } - constexpr int NDimWithoutPhi = 4; // eta, pt, mult, z - constexpr int NDimWithPhi = 5; // eta, pt, mult, z, phi - const int nDim = eff->GetNdimensions(); - if (nDim != NDimWithoutPhi && nDim != NDimWithPhi) { - LOGF(fatal, "Efficiency object has %d dimensions, expected %d (eta,pt,mult,z) or %d (+phi)", nDim, NDimWithoutPhi, NDimWithPhi); - } - const double vals[NDimWithPhi] = {eta, pt, multiplicity, posZ, static_cast(RecoDecay::constrainAngle(phi, 0.f))}; - int bins[NDimWithPhi]; - for (int i = 0; i < nDim; ++i) { - TAxis* ax = eff->GetAxis(i); - bins[i] = ax->FindBin(vals[i]); - if (bins[i] < 1 || bins[i] > ax->GetNbins()) { - return 0.0; - } - } - return eff->GetBinContent(bins); - } - - template - using HasSign = decltype(std::declval().sign()); - template - using HasPdgCode = decltype(std::declval().pdgCode()); - template - using HasIsPhysicalPrimary = decltype(std::declval().isPhysicalPrimary()); - - // Retrieve the multiplicity/centrality estimator selected by cfgMultCorrelationsMask. - // Falls back to FT0C centrality if no mask bit is set. - template - float getMultiplicity(TCollision const& collision) - { - if (cfgMultCorrelationsMask == 0) { - return collision.centFT0C(); - } - if (cfgMultCorrelationsMask & aod::mftmultset::CentFT0C) - return collision.centFT0C(); - if (cfgMultCorrelationsMask & aod::mftmultset::CentFT0M) - return collision.centFT0M(); - return collision.centFT0C(); - } - - template - bool passOutlier(TCollision const& collision) - { - if (cfgMultCutFormula.value.empty() || multCutFormula == nullptr) - return true; - if (multCutFormulaParamIndex[0] != ~0u) - multCutFormula->SetParameter(multCutFormulaParamIndex[0], collision.centFT0C()); - if (multCutFormulaParamIndex[1] != ~0u) - multCutFormula->SetParameter(multCutFormulaParamIndex[1], collision.multFV0A()); - if (multCutFormulaParamIndex[2] != ~0u) - multCutFormula->SetParameter(multCutFormulaParamIndex[2], collision.multNTracksPV()); - if (multCutFormulaParamIndex[3] != ~0u) - multCutFormula->SetParameter(multCutFormulaParamIndex[3], collision.multNTracksPV()); // substitute a "global" estimator if/when available - if (multCutFormulaParamIndex[4] != ~0u) - multCutFormula->SetParameter(multCutFormulaParamIndex[4], collision.centFT0M()); - return multCutFormula->Eval() > 0.0f; - } - - template - void fillMultCorrelations(TCollision const& collision) - { - if (cfgMultCorrelationsMask == 0) - return; - std::vector v; - if (cfgMultCorrelationsMask & aod::mftmultset::CentFT0C) - v.push_back(collision.centFT0C()); - if (cfgMultCorrelationsMask & aod::mftmultset::MultFV0A) - v.push_back(collision.multFV0A()); - if (cfgMultCorrelationsMask & aod::mftmultset::MultNTracksPV) - v.push_back(collision.multNTracksPV()); - if (cfgMultCorrelationsMask & aod::mftmultset::MultNTracksGlobal) - v.push_back(collision.multNTracksPV()); // substitute a "global" estimator if/when available - if (cfgMultCorrelationsMask & aod::mftmultset::CentFT0M) - v.push_back(collision.centFT0M()); - registry.get(HIST("multCorrelations")).get()->Fill(v.data()); - } - - void loadPhiWeight() - { - if (cfg.mPhiWeight != nullptr || cfgPhiWeightFile.value.empty()) { - return; // already loaded, or feature disabled - } - TFile* f = TFile::Open(cfgPhiWeightFile.value.c_str(), "READ"); - if (!f || f->IsZombie()) { - LOGF(fatal, "Could not open phi weight file %s", cfgPhiWeightFile.value.c_str()); - return; - } - cfg.mPhiWeight = dynamic_cast(f->Get(cfgPhiWeightName.value.c_str())); - if (cfg.mPhiWeight) { - cfg.mPhiWeight->SetDirectory(nullptr); // detach so it survives f->Close() - } else { - LOGF(fatal, "Could not find object %s in %s", cfgPhiWeightName.value.c_str(), cfgPhiWeightFile.value.c_str()); - } - f->Close(); - LOGF(info, "Loaded phi weight '%s' from %s (%p)", cfgPhiWeightName.value.c_str(), cfgPhiWeightFile.value.c_str(), static_cast(cfg.mPhiWeight)); - } - - // Looks up a bin-content-based weight from a 1-D phi map. - // Returns 0 if the histogram is missing so the >0 check upstream skips it cleanly. - float getPhiWeight(TH1* hPhiWeight, float phi) - { - if (!hPhiWeight) { - return 0.f; - } - int binX = hPhiWeight->GetXaxis()->FindBin(wrapPhi(phi)); - return hPhiWeight->GetBinContent(binX); - } - - // ------------------------------------------------------------------ QA - // FIX #6b: "yields"/"etaphi" are RAW. "yields_corrected"/"etaphi_corrected" are weighted with the - // efficiency correction factor (associated-particle efficiency object). - template - void fillQA(TCollision const& collision, float multiplicity, TTracks const& tracks) - { - registry.fill(HIST("multiplicity"), multiplicity); - fillMultCorrelations(collision); - const float posZ = collision.posZ(); - - for (const auto& track : tracks) { - if (track.nClusters() < cfgMFTMinNClusters) { - continue; - } - registry.fill(HIST("yields"), multiplicity, track.pt(), track.eta()); - - float phiWeight = 1.f; - if (cfg.mPhiWeight) { - float phiWgt = getPhiWeight(cfg.mPhiWeight, track.phi()); - phiWeight = (phiWgt > 0) ? 1.f / phiWgt : 1.f; - } - registry.fill(HIST("etaphi"), multiplicity, track.eta(), wrapPhi(track.phi()), phiWeight); - registry.fill(HIST("nClusters"), track.nClusters()); - - float effWeight = 1.f; - if (cfg.mEfficiencyAssociated) { - effWeight = getEfficiencyCorrection(cfg.mEfficiencyAssociated, track.eta(), track.pt(), multiplicity, posZ, track.phi()); - } - registry.fill(HIST("yields_corrected"), multiplicity, track.pt(), track.eta(), effWeight); - registry.fill(HIST("etaphi_corrected"), multiplicity, track.eta(), wrapPhi(track.phi()), effWeight); - } - } - - // ------------------------------------------------------------------ event-level acceptance - template - bool acceptCollision(TCollision const& collision) - { - if (!collision.sel8()) { - return false; - } - if (!passOutlier(collision)) { - return false; - } - return true; - } - - // ------------------------------------------------------------------ correlation filling - // Hadron-hadron (MFT-MFT). tracks1 = trigger collection, tracks2 = associated collection. - // Same-event: called as fillCorrelations(..., tracks, tracks, ...) -> self-correlation with index check. - // Mixed-event: called as fillCorrelations(..., tracks1, tracks2, ...) from two different collisions. - - template - void fillCorrelations(TTarget target, TTracks1 const& tracks1, TTracks2 const& tracks2, float multiplicity, bool mixing, float posZ, float eventWeight) - { - // ---------------------------------------------------------------- per-call caches - // Indexed by plain position in the loop (NOT filteredIndex()), and built for BOTH - // collections. This is what makes the mixed event correct: tracks1 and tracks2 have the same - // C++ type in SameKindPair, so a type-based "sameCollection" is true in ME as well. - if constexpr (step == CorrelationContainer::kCFStepCorrected) { - efficiencyAssociatedCache.clear(); - if (cfg.mEfficiencyAssociated) { - efficiencyAssociatedCache.reserve(tracks2.size()); - for (const auto& track : tracks2) { - efficiencyAssociatedCache.push_back(getEfficiencyCorrection(cfg.mEfficiencyAssociated, track.eta(), track.pt(), multiplicity, posZ, track.phi())); - } - } - } - - phiWeightCache.clear(); // associated side (tracks2) - phiWeightCache1.clear(); // trigger side (tracks1) - if (cfg.mPhiWeight) { - phiWeightCache.reserve(tracks2.size()); - for (const auto& track : tracks2) { - phiWeightCache.push_back(getPhiWeight(cfg.mPhiWeight, track.phi())); - } - phiWeightCache1.reserve(tracks1.size()); - for (const auto& track : tracks1) { - phiWeightCache1.push_back(getPhiWeight(cfg.mPhiWeight, track.phi())); - } - } - - // ---------------------------------------------------------------- trigger loop - int i1 = -1; - for (const auto& track1 : tracks1) { - ++i1; // must stay before any 'continue' - if (track1.nClusters() < cfgMFTMinNClusters) { - continue; - } - - if constexpr (std::experimental::is_detected::value) { - if constexpr (std::experimental::is_detected::value) { - if (!track1.isPhysicalPrimary()) - continue; - } - } else if constexpr (std::experimental::is_detected::value) { - if (cfgTriggerCharge != 0 && cfgTriggerCharge * track1.sign() < 0) - continue; - } - - float triggerWeight = eventWeight; - if constexpr (step == CorrelationContainer::kCFStepCorrected) { - if (cfg.mEfficiencyTrigger) { - triggerWeight *= getEfficiencyCorrection(cfg.mEfficiencyTrigger, track1.eta(), track1.pt(), multiplicity, posZ, track1.phi()); - } - } - // phi weight of the trigger: applied BEFORE the trigger histogram is filled, so the - // per-trigger normalisation and the pair weights stay consistent. - if (cfg.mPhiWeight) { - const float w1 = phiWeightCache1[i1]; - if (w1 > 0.f) { - triggerWeight *= 1.f / w1; - } - } - target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, triggerWeight); - - // ---------------------------------------------------------------- associated loop - int i2 = -1; - for (const auto& track2 : tracks2) { - ++i2; // must stay before any 'continue' - - if (!mixing && track1.globalIndex() == track2.globalIndex()) { - continue; // no self-correlation (same event only) - } - if (track2.nClusters() < cfgMFTMinNClusters) { - continue; - } - - if constexpr (std::experimental::is_detected::value) { - if (cfgAssociatedCharge != 0 && cfgAssociatedCharge * track2.sign() < 0) - continue; - } - if constexpr (std::experimental::is_detected::value && - std::experimental::is_detected::value) { - if (cfgPairCharge != 0 && cfgPairCharge * track1.sign() * track2.sign() < 0) - continue; - } - - if (cfgPtOrder != 0 && track2.pt() >= track1.pt()) { - continue; - } - - float associatedWeight = triggerWeight; // already contains event, trigger eff, trigger phi weight - if constexpr (step == CorrelationContainer::kCFStepCorrected) { - if (cfg.mEfficiencyAssociated) { - associatedWeight *= efficiencyAssociatedCache[i2]; - } - } - if (cfg.mPhiWeight) { - const float w2 = phiWeightCache[i2]; - if (w2 > 0.f) { - associatedWeight *= 1.f / w2; - } - } - - const float deltaEta = track1.eta() - track2.eta(); - const float deltaPhi = RecoDecay::constrainAngle(track1.phi() - track2.phi(), -o2::constants::math::PIHalf); - - if (mixing) { - registry.fill(HIST("mftmft/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta); - } else { - registry.fill(HIST("mftmft/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta); - } - - target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, associatedWeight); - } - } - } - // FIX #2 + FIX #6d: MFT-Global. Trigger = central-barrel track (cfg.mEfficiencyGlobal), - // associated = MFT track (cfg.mEfficiencyAssociated). - template - void fillCorrMftGlobal(TTarget target, TTriggers const& triggers, TMFTs const& mftTracks, bool mixing, float posZ, float multiplicity, float eventWeight) - { - constexpr bool TriggerHasSign = std::experimental::is_detected::value; - constexpr bool AssocHasSign = std::experimental::is_detected::value; - - if (!mixing) { - for (const auto& assocTrack : mftTracks) { - if (assocTrack.nClusters() < cfgMFTMinNClusters) { - continue; - } - registry.fill(HIST("mftGlobal/SE/Assoc_etavsphi"), assocTrack.phi(), assocTrack.eta()); - registry.fill(HIST("mftGlobal/SE/Assoc_pt"), assocTrack.pt()); - } - } - - for (const auto& triggerTrack : triggers) { - if constexpr (TriggerHasSign) { - if (cfgTriggerCharge != 0 && cfgTriggerCharge * triggerTrack.sign() < 0) - continue; - } - - float triggerWeight = eventWeight; - if constexpr (step == CorrelationContainer::kCFStepCorrected) { - if (cfg.mEfficiencyGlobal) { - triggerWeight *= getEfficiencyCorrection(cfg.mEfficiencyGlobal, triggerTrack.eta(), triggerTrack.pt(), multiplicity, posZ, triggerTrack.phi()); - } - } - if (!mixing) { - registry.fill(HIST("mftGlobal/SE/Trig_etavsphi"), triggerTrack.phi(), triggerTrack.eta()); - registry.fill(HIST("mftGlobal/SE/Trig_pt"), triggerTrack.pt()); - } - target->getTriggerHist()->Fill(step, triggerTrack.pt(), multiplicity, posZ, triggerWeight); - - for (const auto& assocTrack : mftTracks) { - if (assocTrack.nClusters() < cfgMFTMinNClusters) { - continue; - } - if constexpr (AssocHasSign) { - if (cfgAssociatedCharge != 0 && cfgAssociatedCharge * assocTrack.sign() < 0) - continue; - } - if constexpr (TriggerHasSign && AssocHasSign) { - if (cfgPairCharge != 0 && cfgPairCharge * triggerTrack.sign() * assocTrack.sign() < 0) - continue; - } - float associatedWeight = triggerWeight; - if constexpr (step == CorrelationContainer::kCFStepCorrected) { - if (cfg.mEfficiencyAssociated) { - associatedWeight *= getEfficiencyCorrection(cfg.mEfficiencyAssociated, assocTrack.eta(), assocTrack.pt(), multiplicity, posZ, assocTrack.phi()); - } - } - - float deltaPhi = RecoDecay::constrainAngle(triggerTrack.phi() - assocTrack.phi(), -o2::constants::math::PIHalf); - float deltaEta = triggerTrack.eta() - assocTrack.eta(); - if (mixing) { - registry.fill(HIST("mftGlobal/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta); - } else { - registry.fill(HIST("mftGlobal/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta); - } - target->getPairHist()->Fill(step, deltaEta, assocTrack.pt(), triggerTrack.pt(), multiplicity, deltaPhi, posZ, associatedWeight); - } - } - } - - // FIX #4 + FIX #6d: Global-Global (central-barrel hadron-hadron). Both tracks use cfg.mEfficiencyGlobal. - // No nClusters() cut: global tracks are quality-selected via FilteredGlobalTracks. - template - void fillCorrGlobalGlobal(TTarget target, TTracks1 const& tracks1, TTracks2 const& tracks2, bool mixing, float posZ, float multiplicity, float eventWeight) - { - constexpr bool Track1HasSign = std::experimental::is_detected::value; - constexpr bool Track2HasSign = std::experimental::is_detected::value; - - for (const auto& track1 : tracks1) { - if constexpr (Track1HasSign) { - if (cfgTriggerCharge != 0 && cfgTriggerCharge * track1.sign() < 0) - continue; - } - - float triggerWeight = eventWeight; - if constexpr (step == CorrelationContainer::kCFStepCorrected) { - if (cfg.mEfficiencyGlobal) { - triggerWeight *= getEfficiencyCorrection(cfg.mEfficiencyGlobal, track1.eta(), track1.pt(), multiplicity, posZ, track1.phi()); - } - } - if (!mixing) { - registry.fill(HIST("global/SE/etavsphi"), track1.phi(), track1.eta()); - registry.fill(HIST("global/SE/pt"), track1.pt()); - } - target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, triggerWeight); - - for (const auto& track2 : tracks2) { - if constexpr (std::is_same::value) { - if (track1.globalIndex() == track2.globalIndex()) { - continue; // no self-correlation - } - } - - if constexpr (Track2HasSign) { - if (cfgAssociatedCharge != 0 && cfgAssociatedCharge * track2.sign() < 0) - continue; - } - if constexpr (Track1HasSign && Track2HasSign) { - if (cfgPairCharge != 0 && cfgPairCharge * track1.sign() * track2.sign() < 0) - continue; - } - if (cfgPtOrder != 0 && track2.pt() >= track1.pt()) { - continue; - } - - float associatedWeight = triggerWeight; - if constexpr (step == CorrelationContainer::kCFStepCorrected) { - if (cfg.mEfficiencyGlobal) { - associatedWeight *= getEfficiencyCorrection(cfg.mEfficiencyGlobal, track2.eta(), track2.pt(), multiplicity, posZ, track2.phi()); - } - } - - float deltaPhi = RecoDecay::constrainAngle(track1.phi() - track2.phi(), -o2::constants::math::PIHalf); - float deltaEta = track1.eta() - track2.eta(); - if (mixing) { - registry.fill(HIST("global/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta); - } else { - registry.fill(HIST("global/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta); - } - target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, associatedWeight); - } - } - } - - // ------------------------------------------------------------------ SAME EVENT (data) - void processSameMFT(FilteredCollisions::iterator const& collision, aod::BCsWithTimestamps const&, FilteredMFTTracks const& tracks) - { - registry.fill(HIST("multiplicity_raw"), collision.centFT0C()); // unconditional, no cuts - if (!acceptCollision(collision)) { - return; - } - - loadEfficiency(collision.template bc_as().timestamp()); - const float multiplicity = getMultiplicity(collision); - - using BinningType = ColumnBinningPolicy; - BinningType binning{{axisVertex, axisMultiplicity}, true}; - int bin = binning.getBin({collision.posZ(), multiplicity}); - - if (cfgVerbosity > 0) { - LOGF(info, "processSameMFT: %d MFT tracks | Vertex: %.2f | Mult: %.2f", tracks.size(), collision.posZ(), multiplicity); - } - - registry.fill(HIST("eventcount_same"), bin); - registry.fill(HIST("trackcount_same"), bin, tracks.size()); - fillQA(collision, multiplicity, tracks); - - const bool hasEfficiency = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); - const bool fillReco = !(cfgDropStepRECO && hasEfficiency); - - same->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - if (fillReco) { - same->fillEvent(multiplicity, CorrelationContainer::kCFStepReconstructed); - fillCorrelations(same, tracks, tracks, multiplicity, false, collision.posZ(), 1.0f); - } - if (hasEfficiency) { - same->fillEvent(multiplicity, CorrelationContainer::kCFStepCorrected); - fillCorrelations(same, tracks, tracks, multiplicity, false, collision.posZ(), 1.0f); - } - } - PROCESS_SWITCH(Mftcorrelations, processSameMFT, "Process same event MFT-MFT correlations", true); - - void processMftGlobalSE(FilteredCollisions::iterator const& collision, FilteredGlobalTracks const& globalTracks, FilteredMFTTracks const& mftTracks, aod::BCsWithTimestamps const&) - { - if (!acceptCollision(collision)) { - return; - } - loadEfficiency(collision.template bc_as().timestamp()); - const float multiplicity = getMultiplicity(collision); - - using MixingBinningType = ColumnBinningPolicy; - MixingBinningType binning{{axisVertex, axisMultiplicity}, true}; - int bin = binning.getBin({collision.posZ(), multiplicity}); - registry.fill(HIST("eventcount_mftglobal_same"), bin); - registry.fill(HIST("trackcount_mftglobal_same"), bin, globalTracks.size(), mftTracks.size()); - - // FIX #6d: trigger = global track, associated = MFT track - const bool hasEfficiency = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyGlobal != nullptr); - const bool fillReco = !(cfgDropStepRECO && hasEfficiency); - - sameMftGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - if (fillReco) { - sameMftGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepReconstructed); - fillCorrMftGlobal(sameMftGlobal, globalTracks, mftTracks, false, collision.posZ(), multiplicity, 1.0f); - } - if (hasEfficiency) { - sameMftGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepCorrected); - fillCorrMftGlobal(sameMftGlobal, globalTracks, mftTracks, false, collision.posZ(), multiplicity, 1.0f); - } - } - PROCESS_SWITCH(Mftcorrelations, processMftGlobalSE, "Process same event MFT-vs-global-track correlations", false); - - // FIX #4: Same-event global-global (central-barrel hadron-hadron) correlations. - void processGlobalGlobalSE(FilteredCollisions::iterator const& collision, FilteredGlobalTracks const& tracks, aod::BCsWithTimestamps const&) - { - if (!acceptCollision(collision)) { - return; - } - loadEfficiency(collision.template bc_as().timestamp()); - const float multiplicity = getMultiplicity(collision); - - using MixingBinningType = ColumnBinningPolicy; - MixingBinningType binning{{axisVertex, axisMultiplicity}, true}; - int bin = binning.getBin({collision.posZ(), multiplicity}); - - if (cfgVerbosity > 0) { - LOGF(info, "processGlobalGlobalSE: %d global tracks | Vertex: %.2f | Mult: %.2f", tracks.size(), collision.posZ(), multiplicity); - } - - registry.fill(HIST("eventcount_global_same"), bin); - registry.fill(HIST("trackcount_global_same"), bin, tracks.size()); - - // FIX #6d - const bool hasEfficiency = (cfg.mEfficiencyGlobal != nullptr); - const bool fillReco = !(cfgDropStepRECO && hasEfficiency); - - sameGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - if (fillReco) { - sameGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepReconstructed); - fillCorrGlobalGlobal(sameGlobal, tracks, tracks, false, collision.posZ(), multiplicity, 1.0f); - } - if (hasEfficiency) { - sameGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepCorrected); - fillCorrGlobalGlobal(sameGlobal, tracks, tracks, false, collision.posZ(), multiplicity, 1.0f); - } - } - PROCESS_SWITCH(Mftcorrelations, processGlobalGlobalSE, "Process same event global-global (central-barrel) correlations", false); - - void processSameMFTMultSet(FilteredCollisions::iterator const& collision, aod::BCsWithTimestamps const& bcs, FilteredMFTTracks const& tracks) - { - if (cfgMultCorrelationsMask == 0) { - LOGF(fatal, "cfgMultCorrelationsMask can not be 0 when the MultSet process function is in use."); - } - processSameMFT(collision, bcs, tracks); - } - PROCESS_SWITCH(Mftcorrelations, processSameMFTMultSet, "Process same event MFT-MFT correlations with multi-estimator QA", false); - - // ------------------------------------------------------------------ MIXED EVENT (data) - void processMixedMFT(FilteredCollisions const& collisions, FilteredMFTTracks const& tracks, aod::BCsWithTimestamps const&) - { - // Bin mixed events on the SAME multiplicity/centrality estimator that getMultiplicity() returns - // (driven by cfgMultCorrelationsMask), instead of a hardcoded CentFT0C column. This keeps the - // mixing pool synchronized with whatever estimator is actually filled into the CorrelationContainer, - // and lets you switch estimators (e.g. FT0C -> FT0M) purely from the JSON config. - auto getMult = [this](auto const& c) { return getMultiplicity(c); }; - using BinningType = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getMult)>; - BinningType binning{{getMult}, {axisVertex, axisMultiplicity}, true}; - - auto tracksTuple = std::make_tuple(tracks); - SameKindPair pairs{binning, cfgNoMixedEvents, -1, collisions, tracksTuple, &cache}; - - for (auto it = pairs.begin(); it != pairs.end(); it++) { - auto& [collision1, tracks1, collision2, tracks2] = *it; - - if (!acceptCollision(collision1) || !collision2.sel8()) { - continue; - } - - const float multiplicity = getMultiplicity(collision1); - int bin = binning.getBin(std::tuple(collision1.posZ(), multiplicity)); - float eventWeight = 1.0f / it.currentWindowNeighbours(); - - if (cfgVerbosity > 0) { - LOGF(info, "processMixedMFT: bin %d | pair (%d,%d) | Vertex %.2f/%.2f | Mult %.2f", bin, collision1.globalIndex(), collision2.globalIndex(), collision1.posZ(), collision2.posZ(), multiplicity); - } - - if (it.isNewWindow()) { - loadEfficiency(collision1.template bc_as().timestamp()); - mixed->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - } - - registry.fill(HIST("eventcount_mixed"), bin); - registry.fill(HIST("trackcount_mixed"), bin, tracks1.size(), tracks2.size()); - - const bool hasEfficiency = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); - const bool fillReco = !(cfgDropStepRECO && hasEfficiency); - - if (it.isNewWindow() && fillReco) { - mixed->fillEvent(multiplicity, CorrelationContainer::kCFStepReconstructed); - } - if (it.isNewWindow() && hasEfficiency) { - mixed->fillEvent(multiplicity, CorrelationContainer::kCFStepCorrected); - } - - if (fillReco) { - fillCorrelations(mixed, tracks1, tracks2, multiplicity, true, collision1.posZ(), eventWeight); - } - if (hasEfficiency) { - fillCorrelations(mixed, tracks1, tracks2, multiplicity, true, collision1.posZ(), eventWeight); - } - } - } - PROCESS_SWITCH(Mftcorrelations, processMixedMFT, "Process mixed event MFT-MFT correlations", false); - - void processMftGlobalME(FilteredCollisions const& collisions, FilteredGlobalTracks const& globalTracks, FilteredMFTTracks const& mftTracks, aod::BCsWithTimestamps const&) - { - // Same fix as processMixedMFT: bin on getMultiplicity() (driven by cfgMultCorrelationsMask) - // instead of a hardcoded CentFT0C column, so the mixing pool matches whatever estimator is - // actually filled into the CorrelationContainer. - auto getMult = [this](auto const& c) { return getMultiplicity(c); }; - using MixingBinningType = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getMult)>; - MixingBinningType binning{{getMult}, {axisVertex, axisMultiplicity}, true}; - - auto tracksTuple = std::make_tuple(globalTracks, mftTracks); - Pair pairs{binning, cfgNoMixedEvents, -1, collisions, tracksTuple, &cache}; - - int pairCounter = 0; - int acceptedPairCounter = 0; - int newWindowCounter = 0; - - for (auto it = pairs.begin(); it != pairs.end(); it++) { - auto& [collision1, globalTracks1, collision2, mftTracks2] = *it; - pairCounter++; - - if (!acceptCollision(collision1) || !collision2.sel8()) { - continue; - } - acceptedPairCounter++; - - const float multiplicity = getMultiplicity(collision1); - float eventWeight = 1.0f / it.currentWindowNeighbours(); - - const bool hasEfficiency = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyGlobal != nullptr); - const bool fillReco = !(cfgDropStepRECO && hasEfficiency); - - int bin = binning.getBin(std::tuple(collision1.posZ(), multiplicity)); - - if (it.isNewWindow()) { - newWindowCounter++; - loadEfficiency(collision1.template bc_as().timestamp()); - mixedMftGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - if (fillReco) { - mixedMftGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepReconstructed); - } - if (hasEfficiency) { - mixedMftGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepCorrected); - } - registry.fill(HIST("eventcount_mftglobal_mixed"), bin); - } - registry.fill(HIST("trackcount_mftglobal_mixed"), bin, globalTracks1.size(), mftTracks2.size()); - - if (fillReco) { - fillCorrMftGlobal(mixedMftGlobal, globalTracks1, mftTracks2, true, collision1.posZ(), multiplicity, eventWeight); - } - if (hasEfficiency) { - fillCorrMftGlobal(mixedMftGlobal, globalTracks1, mftTracks2, true, collision1.posZ(), multiplicity, eventWeight); - } - } - - if (cfgVerbosity > 0) { - LOGF(info, "processMftGlobalME: collisions %d | pairs iterated %d | accepted %d | new windows %d", - collisions.size(), pairCounter, acceptedPairCounter, newWindowCounter); - } - } - PROCESS_SWITCH(Mftcorrelations, processMftGlobalME, "Process mixed event MFT-vs-global-track correlations", false); - - // FIX #4: Mixed-event global-global (central-barrel hadron-hadron) correlations. - void processGlobalGlobalME(FilteredCollisions const& collisions, FilteredGlobalTracks const& tracks, aod::BCsWithTimestamps const&) - { - auto getMult = [this](auto const& c) { return getMultiplicity(c); }; - using MixingBinningType = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getMult)>; - MixingBinningType binning{{getMult}, {axisVertex, axisMultiplicity}, true}; - - auto tracksTuple = std::make_tuple(tracks); - SameKindPair pairs{binning, cfgNoMixedEvents, -1, collisions, tracksTuple, &cache}; - - for (auto it = pairs.begin(); it != pairs.end(); it++) { - auto& [collision1, tracks1, collision2, tracks2] = *it; - - if (!acceptCollision(collision1) || !collision2.sel8()) { - continue; - } - - const float multiplicity = getMultiplicity(collision1); - int bin = binning.getBin(std::tuple(collision1.posZ(), multiplicity)); - float eventWeight = 1.0f / it.currentWindowNeighbours(); - - if (cfgVerbosity > 0) { - LOGF(info, "processGlobalGlobalME: bin %d | pair (%d,%d) | Vertex %.2f/%.2f | Mult %.2f", bin, collision1.globalIndex(), collision2.globalIndex(), collision1.posZ(), collision2.posZ(), multiplicity); - } - - const bool hasEfficiency = (cfg.mEfficiencyGlobal != nullptr); - const bool fillReco = !(cfgDropStepRECO && hasEfficiency); - - if (it.isNewWindow()) { - loadEfficiency(collision1.template bc_as().timestamp()); - mixedGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - if (fillReco) { - mixedGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepReconstructed); - } - if (hasEfficiency) { - mixedGlobal->fillEvent(multiplicity, CorrelationContainer::kCFStepCorrected); - } - registry.fill(HIST("eventcount_global_mixed"), bin); - } - registry.fill(HIST("trackcount_global_mixed"), bin, tracks1.size(), tracks2.size()); - - if (fillReco) { - fillCorrGlobalGlobal(mixedGlobal, tracks1, tracks2, true, collision1.posZ(), multiplicity, eventWeight); - } - if (hasEfficiency) { - fillCorrGlobalGlobal(mixedGlobal, tracks1, tracks2, true, collision1.posZ(), multiplicity, eventWeight); - } - } - } - PROCESS_SWITCH(Mftcorrelations, processGlobalGlobalME, "Process mixed event global-global (central-barrel) correlations", false); - - void processMixedMFTMultSet(FilteredCollisions const& collisions, FilteredMFTTracks const& tracks, aod::BCsWithTimestamps const& bcs) - { - if (cfgMultCorrelationsMask == 0) { - LOGF(fatal, "cfgMultCorrelationsMask can not be 0 when the MultSet process function is in use."); - } - processMixedMFT(collisions, tracks, bcs); - } - PROCESS_SWITCH(Mftcorrelations, processMixedMFTMultSet, "Process mixed event MFT-MFT correlations with multi-estimator QA", false); - - template - bool isChargedPrimary(TParticle const& p) - { - if (!p.isPhysicalPrimary()) { - return false; - } - auto* entry = pdg->GetParticle(p.pdgCode()); - return entry != nullptr && entry->Charge() != 0; - } - - // ------------------------------------------------------------------ MC efficiency extraction - Preslice perCollisionMFT = aod::fwdtrack::collisionId; - - void processMCEfficiencyMFT(aod::McCollisions::iterator const& mcCollision, - aod::McParticles const& mcParticles, - McEffCollisions const& collisions, - McEffMFTTracks const& tracks) - { - registry.fill(HIST("MCEff/hNRecoColl"), collisions.size()); - registry.fill(HIST("MCEff/hCutFlow"), 1); - if (collisions.size() != 1) { - return; - } - registry.fill(HIST("MCEff/hCutFlow"), 2); - if (cfgVerbosity > 0) { - LOGF(info, "processMCEfficiencyMFT: ENTER | mcCollision=%lld | nMcParticles=%d | nCollisions=%d | nTracks=%d", - mcCollision.globalIndex(), mcParticles.size(), collisions.size(), tracks.size()); - } - float multiplicity = -1; - bool passSel8 = false; - int64_t collGlobalIndex = -1; - for (const auto& collision : collisions) { - if (!collision.sel8()) { - continue; - } - registry.fill(HIST("MCEff/hCutFlow"), 3); - passSel8 = true; - multiplicity = getMultiplicity(collision); // driven by cfgMultCorrelationsMask - collGlobalIndex = collision.globalIndex(); - } - if (!passSel8) { - return; - } - - for (const auto& mcParticle : mcParticles) { - if (!isChargedPrimary(mcParticle) || mcParticle.eta() < cfgMFTCutEtaMin || mcParticle.eta() > cfgMFTCutEtaMax || mcParticle.pt() < cfgMFTCutPtMin || mcParticle.pt() > cfgMFTCutPtMax) { - continue; - } - same->getTrackHistEfficiency()->Fill(CorrelationContainer::MC, mcParticle.eta(), mcParticle.pt(), 0, multiplicity, mcCollision.posZ(), wrapPhi(mcParticle.phi())); - } - - auto groupedTracks = tracks.sliceBy(perCollisionMFT, collGlobalIndex); - - std::unordered_set countedPrimaries; - std::unordered_set countedAll; - - for (const auto& track : groupedTracks) { - if (track.nClusters() < cfgMFTMinNClusters) { - continue; - } - if (track.has_mcParticle()) { - const auto& mcParticle = track.mcParticle(); - const int64_t idx = mcParticle.globalIndex(); - - if (mcParticle.isPhysicalPrimary() && countedPrimaries.insert(idx).second) { - same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoPrimaries, - mcParticle.eta(), mcParticle.pt(), 0, multiplicity, mcCollision.posZ(), wrapPhi(mcParticle.phi())); - } - if (countedAll.insert(idx).second) { - same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoAll, - mcParticle.eta(), mcParticle.pt(), 0, multiplicity, mcCollision.posZ(), wrapPhi(mcParticle.phi())); - } - } else { - // FIX #6e: wrap phi to [0,2pi) like the MC-particle phi used for the other steps - same->getTrackHistEfficiency()->Fill(CorrelationContainer::Fake, - track.eta(), track.pt(), 0, multiplicity, mcCollision.posZ(), - RecoDecay::constrainAngle(track.phi(), 0.f)); - } - } - } - PROCESS_SWITCH(Mftcorrelations, processMCEfficiencyMFT, "MC: Extract MFT tracking efficiency", false); - - //-------------------------------------------------------------------MC Efficiency Extraction for Global - - Preslice perCollisionGlobal = aod::track::collisionId; - - void processMCEfficiencyGlobal(aod::McCollisions::iterator const& mcCollision, - aod::McParticles const& mcParticles, - McEffCollisions const& collisions, - McEffGlobalTracks const& tracks) - { - registry.fill(HIST("MCEffGlobal/hNRecoColl"), collisions.size()); - registry.fill(HIST("MCEffGlobal/hCutFlow"), 1); - if (collisions.size() != 1) { // exactly one reconstructed collision per MC collision (as in the MFT efficiency) - return; - } - registry.fill(HIST("MCEffGlobal/hCutFlow"), 2); - if (cfgVerbosity > 0) { - LOGF(info, "processMCEfficiencyGlobal: mcCollision=%lld | nMcParticles=%d | nCollisions=%d | nTracks=%d", - mcCollision.globalIndex(), mcParticles.size(), collisions.size(), tracks.size()); - } - - float multiplicity = -1; - bool passSel8 = false; - int64_t collGlobalIndex = -1; - for (const auto& collision : collisions) { - if (!collision.sel8()) { - continue; - } - registry.fill(HIST("MCEffGlobal/hCutFlow"), 3); - passSel8 = true; - multiplicity = getMultiplicity(collision); // same estimator as in data (cfgMultCorrelationsMask) - collGlobalIndex = collision.globalIndex(); - } - if (!passSel8) { - return; - } - - // Generated: charged physical primaries in the global-track acceptance - for (const auto& mcParticle : mcParticles) { - if (!isChargedPrimary(mcParticle) || - std::abs(mcParticle.eta()) >= cfgGlobalTrackEtaCut || - mcParticle.pt() < cfgGlobalTrackPtMin || mcParticle.pt() > cfgGlobalTrackPtMax) { - continue; - } - sameGlobal->getTrackHistEfficiency()->Fill(CorrelationContainer::MC, mcParticle.eta(), mcParticle.pt(), 0, multiplicity, mcCollision.posZ(), mcParticle.phi()); - } - - // Reconstructed: global tracks of the selected reco collision (globalTrackFilter already applied) - auto groupedTracks = tracks.sliceBy(perCollisionGlobal, collGlobalIndex); - - std::unordered_set countedPrimaries; // count each MC particle once (split tracks) - std::unordered_set countedAll; - - for (const auto& track : groupedTracks) { - if (track.has_mcParticle()) { - const auto& mcParticle = track.mcParticle(); - const int64_t idx = mcParticle.globalIndex(); - - if (mcParticle.isPhysicalPrimary() && countedPrimaries.insert(idx).second) { - sameGlobal->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoPrimaries, - mcParticle.eta(), mcParticle.pt(), 0, multiplicity, mcCollision.posZ(), mcParticle.phi()); - } - if (countedAll.insert(idx).second) { - sameGlobal->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoAll, - mcParticle.eta(), mcParticle.pt(), 0, multiplicity, mcCollision.posZ(), mcParticle.phi()); - } - } else { - // fake track - sameGlobal->getTrackHistEfficiency()->Fill(CorrelationContainer::Fake, - track.eta(), track.pt(), 0, multiplicity, mcCollision.posZ(), - RecoDecay::constrainAngle(track.phi(), 0.f)); - } - } - } - PROCESS_SWITCH(Mftcorrelations, processMCEfficiencyGlobal, "MC: Extract central-barrel (global track) tracking efficiency", false); - - // ------------------------------------------------------------------ MC generator-level correlations (truth reference / closure test) - using McRecoCollisions = soa::Join; - - // -1, 0, +1 from the PDG database charge (0 = neutral or unknown code) - template - int getMcCharge(TParticle const& p) - { - auto* entry = pdg->GetParticle(p.pdgCode()); - if (entry == nullptr) { - return 0; - } - const double q = entry->Charge(); - return (q > 0) - (q < 0); - } - - // cheap cuts first, PDG lookup last - template - bool acceptMftTruth(TParticle const& p) - { - return p.isPhysicalPrimary() && p.eta() > cfgMFTCutEtaMin && p.eta() < cfgMFTCutEtaMax && - p.pt() > cfgMFTCutPtMin && p.pt() < cfgMFTCutPtMax && getMcCharge(p) != 0; - } - - template - bool acceptGlobalTruth(TParticle const& p) - { - return p.isPhysicalPrimary() && std::abs(p.eta()) < cfgGlobalTrackEtaCut && - p.pt() > cfgGlobalTrackPtMin && p.pt() < cfgGlobalTrackPtMax && getMcCharge(p) != 0; - } - - template - std::vector collectAccepted(TParticles const& particles, TSelector accept) - { - std::vector idx; - idx.reserve(64); - for (const auto& p : particles) { - if (accept(p)) { - idx.push_back(p.globalIndex()); - } - } - return idx; - } - - template - void fillTruthPairs(TTarget& target, TParticles1 const& particles1, TParticles2 const& particles2, - std::vector const& acc1, std::vector const& acc2, - float multiplicity, float posZ, float weight, bool sameCollection) - { - for (const auto& i1 : acc1) { - const auto& p1 = particles1.iteratorAt(i1 - particles1.begin().globalIndex()); - target->getTriggerHist()->Fill(CorrelationContainer::kCFStepAll, p1.pt(), multiplicity, posZ, weight); - for (const auto& i2 : acc2) { - if (sameCollection && i1 == i2) { - continue; - } - const auto& p2 = particles2.iteratorAt(i2 - particles2.begin().globalIndex()); - if (cfgPtOrder != 0 && p2.pt() >= p1.pt()) { - continue; - } - float deltaPhi = RecoDecay::constrainAngle(p1.phi() - p2.phi(), -o2::constants::math::PIHalf); - target->getPairHist()->Fill(CorrelationContainer::kCFStepAll, p1.eta() - p2.eta(), p2.pt(), p1.pt(), multiplicity, deltaPhi, posZ, weight); - } - } - } - - // ---------------- same event - template - void fillTruthSame(TTarget& target, aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles, - soa::SmallGroups const& collisions, TSelector accept) - { - if (collisions.size() != 1) { // same reco-collision requirement as processMCEfficiencyMFT - return; - } - auto collision = collisions.begin(); - if (!collision.sel8()) { - return; - } - const float multiplicity = getMultiplicity(collision); - const auto acc = collectAccepted(mcParticles, accept); - - target->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - fillTruthPairs(target, mcParticles, mcParticles, acc, acc, multiplicity, mcCollision.posZ(), 1.0f, true); - } - - void processMCSameMFT(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles, - soa::SmallGroups const& collisions) - { - fillTruthSame(same, mcCollision, mcParticles, collisions, [this](auto const& p) { return acceptMftTruth(p); }); - } - PROCESS_SWITCH(Mftcorrelations, processMCSameMFT, "Process MC-truth same event MFT-acceptance correlations (charged primaries)", false); - - void processMCSameGlobal(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles, - soa::SmallGroups const& collisions) - { - fillTruthSame(sameGlobal, mcCollision, mcParticles, collisions, [this](auto const& p) { return acceptGlobalTruth(p); }); - } - PROCESS_SWITCH(Mftcorrelations, processMCSameGlobal, "Process MC-truth same event global-acceptance correlations (charged primaries)", false); - - // ---------------- mixed event - PresliceUnsorted perMcCollision = aod::mccollisionlabel::mcCollisionId; - - template - void fillTruthMixed(TTarget& target, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, - McRecoCollisions const& collisions, TSelector accept) - { - // multiplicity of an MC collision = that of its single, sel8-passing reco collision; otherwise -1 (skipped by the binning) - auto getMult = [&collisions, this](auto const& mcCol) { - auto grouped = collisions.sliceBy(perMcCollision, mcCol.globalIndex()); - if (grouped.size() != 1) { - return -1.0f; - } - auto c = grouped.begin(); - if (!c.sel8()) { - return -1.0f; - } - return static_cast(getMultiplicity(c)); - }; - - using BinningTypeMC = FlexibleBinningPolicy, aod::mccollision::PosZ, decltype(getMult)>; - BinningTypeMC binning{{getMult}, {axisVertex, axisMultiplicity}, true}; - auto tuple = std::make_tuple(mcParticles); - SameKindPair pairs{binning, cfgNoMixedEvents, -1, mcCollisions, tuple, &cache}; - - for (auto it = pairs.begin(); it != pairs.end(); it++) { - auto& [mcCollision1, particles1, mcCollision2, particles2] = *it; - const float eventWeight = 1.0f / it.currentWindowNeighbours(); - const float multiplicity = getMult(mcCollision1); - - if (it.isNewWindow()) { - target->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); - } - const auto acc1 = collectAccepted(particles1, accept); - const auto acc2 = collectAccepted(particles2, accept); // once per pair, not once per trigger - fillTruthPairs(target, particles1, particles2, acc1, acc2, multiplicity, mcCollision1.posZ(), eventWeight, false); - } - } - - void processMCMixedMFT(aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, McRecoCollisions const& collisions) - { - fillTruthMixed(mixed, mcCollisions, mcParticles, collisions, [this](auto const& p) { return acceptMftTruth(p); }); - } - PROCESS_SWITCH(Mftcorrelations, processMCMixedMFT, "Process MC-truth mixed event MFT-acceptance correlations (charged primaries)", false); - - void processMCMixedGlobal(aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, McRecoCollisions const& collisions) - { - fillTruthMixed(mixedGlobal, mcCollisions, mcParticles, collisions, [this](auto const& p) { return acceptGlobalTruth(p); }); - } - PROCESS_SWITCH(Mftcorrelations, processMCMixedGlobal, "Process MC-truth mixed event global-acceptance correlations (charged primaries)", false); -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; -}