diff --git a/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx b/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx index c76556ec26f..7e414d4cc96 100644 --- a/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx @@ -13,29 +13,44 @@ /// \author Salman Malik /// \author Balwan Singh -#include -#include -#include #include "TRandom.h" +#include + +#include // O2 includes +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" +#include "Framework/O2DatabasePDGPlugin.h" #include "Framework/runDataProcessing.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Centrality.h" -#include "Common/Core/TrackSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" #include "ReconstructionDataFormats/GlobalTrackID.h" #include "ReconstructionDataFormats/Track.h" -#include "Framework/ASoAHelpers.h" + +#include "TDatabasePDG.h" +#include + +#include +#include +#include using std::array; using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; + +using namespace std; + +TH1D* tmpFqErr[6][5][52]; + struct FactorialMoments { + Configurable confEta{"centralEta", 0.9, "eta limit for tracks"}; - Configurable confNumPt{"numPt", 1, "number of pT bins"}; + Configurable confNumPt{"numPt", 5, "number of pT bins"}; Configurable confPtMin{"ptMin", 0.2f, "lower pT cut"}; Configurable confDCAxy{"dcaXY", 2.4f, "DCA xy cut"}; Configurable confDCAz{"dcaZ", 2.0f, "DCA z cut"}; @@ -53,10 +68,43 @@ struct FactorialMoments { Configurable includeGlobalTracks{"includeGlobalTracks", false, "Enable Global Tracks"}; Configurable includeTPCTracks{"includeTPCTracks", false, "TPC Tracks"}; Configurable includeITSTracks{"includeITSTracks", false, "ITS Tracks"}; + Configurable confSamplesize{"samplesize", 100, "Sample size"}; + Configurable confUseMC{"useMC", false, "Use MC information"}; + Configurable reduceOutput{"reduce-output", 0, "Suppress info level output (0 = all output, 1 = per collision, 2 = none)"}; Filter filterTracks = (nabs(aod::track::eta) < confEta) && (aod::track::pt >= confPtMin) && (nabs(aod::track::dcaXY) < confDCAxy) && (nabs(aod::track::dcaZ) < confDCAz); Filter filterCollisions = (nabs(aod::collision::posZ) < confVertex.value[2]) && (nabs(aod::collision::posX) < confVertex.value[0]) && (nabs(aod::collision::posY) < confVertex.value[1]); + Service pdg; // Histograms + HistogramRegistry histos1{ + "histos1", + { + {"hRecoPtBefore", "Reco pT before cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, + {"hGenPtBefore", "Gen pT before cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, + {"hRecoPtAfter", "Reco pT after cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, + {"hGenPtAfter", "Gen pT after cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, + {"hRecoEtaBefore", "Reco #eta before cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, + {"mCentMCImpactFull", "Impact Parameter (MC All Events);Impact Parameter b [fm];Counts", {HistType::kTH1F, {{200, 0.0, 50.0}}}}, + {"mCentMCImpact005", "Impact Parameter (MC 0–5%);Impact Parameter b [fm];Counts", {HistType::kTH1F, {{200, 0.0, 20.0}}}}, + {"hGenEtaBefore", "Gen #eta before cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, + {"hRecoEtaAfter", "Reco #eta after cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, + {"hGenEtaAfter", "Gen #eta after cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, + {"hRecoPhiBefore", "Reco #phi before cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, + {"hGenPhiBefore", "Gen #phi before cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, + {"hRecoPhiAfter", "Reco #phi after cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, + {"hGenCharge", "Gen particle charge;Charge;Counts", {HistType::kTH1F, {{13, -6.5, 6.5}}}}, + {"hGenPhiAfter", "Gen #phi after cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, + {"mChargeDist", "MC particle charge;Charge;Counts", {HistType::kTH1F, {{9, -4.5, 4.5}}}}, + {"mDiffPhi", "Δ#phi between selected MC particles;Δ#phi;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, + {"mPrimariesPerEvent", "Primary MC particles per event;N_{primaries};Counts", {HistType::kTH1I, {{20000, 0, 20000}}}}, + {"mDiffPt", "Δp_{T} between selected MC particles;Δp_{T} (GeV/c);Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, + {"mDiffEta", "Δ#eta between selected MC particles;Δ#eta;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, + {"mADiffPhi", "Δ#phi between selected MC particles;Δ#phi;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, + {"mADiffPt", "Δp_{T} between selected MC particles;Δp_{T} (GeV/c);Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, + {"mADiffEta", "Δ#eta between selected MC particles;Δ#eta;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, + }, + }; + HistogramRegistry histos{ "histos", { @@ -94,16 +142,26 @@ struct FactorialMoments { OutputObjHandlingPolicy::AnalysisObject, true}; static const Int_t nBins = 52; + Int_t countSamples = 0; + Int_t testc1 = 0, testc2 = 0, testc3 = 0; array binningM; array countTracks{0, 0, 0, 0, 0}; array, 5>, 6> fqEvent; + array, 5>, 6> fqEventSampled; array, 5> binConEvent; + array, 5>, 6> binConEventSampled; + array, 5>, 6> errorFq = {{{{{0, 0, 0, 0, 0}}}}}; std::vector> mHistArrReset; std::vector> mHistArrQA; std::vector> mFqBinFinal; std::vector> mBinConFinal; + std::vector> mFqBinFinalSampled; + std::vector> mBinConFinalSampled; + std::vector> mFqError; + // max number of bins restricted to 5 - static constexpr array mbinNames{"bin1/", "bin2/", "bin3/", "bin4/", "bin5/"}; + static constexpr array + mbinNames{"bin1/", "bin2/", "bin3/", "bin4/", "bin5/"}; void init(o2::framework::InitContext&) { // NOTE: check to make number of pt and the vector consistent @@ -114,52 +172,159 @@ struct FactorialMoments { } } AxisSpec axisPt[5] = {{100, -0.01, 3 * confPtBins.value[1], ""}, {100, -0.01, 3 * confPtBins.value[3], ""}, {100, -0.01, 3 * confPtBins.value[5], ""}, {100, -0.01, 3 * confPtBins.value[7], ""}, {100, -0.01, 3 * confPtBins.value[9], ""}}; // pT axis + auto mEventSelected = std::get>(histos.add("mEventSelected", "eventSelected", HistType::kTH1D, {{8, 0.5, 8.5}})); + mEventSelected->GetXaxis()->SetBinLabel(1, "all"); + mEventSelected->GetXaxis()->SetBinLabel(2, "sel8"); + mEventSelected->GetXaxis()->SetBinLabel(3, "sameBunchPileup"); + mEventSelected->GetXaxis()->SetBinLabel(4, "goodZvtxFT0vsPV"); + mEventSelected->GetXaxis()->SetBinLabel(5, "vertexITSTPC"); + mEventSelected->GetXaxis()->SetBinLabel(6, "centrality"); + mEventSelected->GetXaxis()->SetBinLabel(7, "final"); + auto mTrackSelected = std::get>(histos.add( + "mTrackSelected", "Track Selection Steps", HistType::kTH1D, {{5, 0.5, 5.5}})); + + mTrackSelected->GetXaxis()->SetBinLabel(1, "all"); + mTrackSelected->GetXaxis()->SetBinLabel(2, "charge"); + mTrackSelected->GetXaxis()->SetBinLabel(3, "PIDs"); + mTrackSelected->GetXaxis()->SetBinLabel(4, "Pids exclude pions"); + mTrackSelected->GetXaxis()->SetBinLabel(5, "Pids exclude pions + kions"); + if (confUseMC) { + auto mMcTrackSelected = std::get>(histos.add("mMcTrackSelected", "mcTrackSelection", HistType::kTH1D, {{5, 0.5, 5.5}})); + } for (auto iM = 0; iM < nBins; ++iM) { binningM[iM] = 2 * (iM + 2); } + for (auto iPt = 0; iPt < confNumPt; ++iPt) { mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mEta", iPt + 1), Form("#eta for bin %.2f-%.2f;#eta", confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{1000, -2, 2}}))); mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mPt", iPt + 1), Form("pT for bin %.2f-%.2f;pT", confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {axisPt[iPt]}))); mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mPhi", iPt + 1), Form("#phi for bin %.2f-%.2f;#phi", confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, 2 * TMath::Pi()}}))); - mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mMultiplicity", iPt + 1), Form("Multiplicity for bin %.2f-%.2f;Multiplicity", confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, 8000}}))); + mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mMultiplicity", iPt + 1), Form("Multiplicity for bin %.2f-%.2f;Multiplicity", confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, 15000}}))); for (auto iM = 0; iM < nBins; ++iM) { auto mHistsR = std::get>(histos.add(Form("bin%i/Reset/mEtaPhi%i", iPt + 1, iM), Form("#eta#phi_%i for bin %.2f-%.2f;#eta;#phi", iM, confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH2F, {{binningM[iM], -0.8, 0.8}, {binningM[iM], 0, 2 * TMath::Pi()}})); mHistArrReset.push_back(mHistsR); + + for (auto iq = 0; iq < 6; ++iq) { + tmpFqErr[iq][iPt][iM] = new TH1D(Form("tmpFqErr%i%i%i", iq, iPt, iM), Form("tmpFqErr%i%i%i", iq, iPt, iM), 100, 0, 10); + } } for (auto i = 0; i < 6; ++i) { auto mHistFq = std::get>(histos.add(Form("mFinalFq%i_bin%i", i + 2, iPt + 1), Form("Final F_%i for bin %.2f-%.2f;M", i + 2, confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); mFqBinFinal.push_back(mHistFq); auto mHistAv = std::get>(histos.add(Form("mFinalAvBin%i_bin%i", i + 2, iPt + 1), Form("Final AvBin_%i for bin %.2f-%.2f;M", i + 2, confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); mBinConFinal.push_back(mHistAv); + + auto mHistFqSampled = std::get>(histos.add(Form("mFinalFq%iSampled_bin%i", i + 2, iPt + 1), Form("Final F_%i for bin %.2f-%.2f;M", i + 2, confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); + mFqBinFinalSampled.push_back(mHistFqSampled); + auto mHistAvSampled = std::get>(histos.add(Form("mFinalAvBin%iSampled_bin%i", i + 2, iPt + 1), Form("Final AvBin_%i for bin %.2f-%.2f;M", i + 2, confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); + mBinConFinalSampled.push_back(mHistAvSampled); + + auto mHistError = std::get>(histos.add(Form("mFqError%i_bin%i", i + 2, iPt + 1), Form("Error for F_%i for bin %.2f-%.2f;M", i + 2, confPtBins.value[2 * iPt], confPtBins.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); + mFqError.push_back(mHistError); + } + } + + array charge{"pos", "neg", "all"}; + array detadphiRange{"", "close"}; + array ptOrd{"pt2>pt1_dphi<0", "pt2>pt1_dphi>0", "pt20"}; + for (Int_t i = 0; i < confNumPt; ++i) { + mHistArrQA.push_back(std::get>(histos.add(Form("mDetaDphiBef%i", i), Form("dEta dPhi for ptbin %d;#Delta#eta;#Delta#phi", i), HistType::kTH2F, {{35, -1.75, 1.75}, {50, -TMath::Pi(), TMath::Pi()}}))); + mHistArrQA.push_back(std::get>(histos.add(Form("mDetaDphiAf%i", i), Form("dEta dPhi for ptbin %d;#Delta#eta;#Delta#phi", i), HistType::kTH2F, {{50, -1.75, 1.75}, {50, -TMath::Pi(), TMath::Pi()}}))); + + for (Int_t ch = 0; ch < 3; ++ch) { + std::string histName = Form("mDetaDphiCh_%.2f-%.2f_%s%s_%d", confPtBins.value[2 * i], confPtBins.value[2 * i + 1], charge[ch].c_str(), detadphiRange[0].c_str(), i); + std::replace(histName.begin(), histName.end(), '.', '_'); + + mHistArrQA.push_back(std::get>(histos.add(histName.c_str(), Form("dEta dPhi for ptbin %d %s;#Delta#eta;#Delta#phi", i, charge[ch].data()), HistType::kTH2F, {{35, -1.75, 1.75}, {50, -TMath::Pi(), TMath::Pi()}}))); + std::string histName2 = Form("mDetaDphiCh_%.2f-%.2f_%s%s_%d", confPtBins.value[2 * i], confPtBins.value[2 * i + 1], charge[ch].data(), detadphiRange[1].data(), i); + std::replace(histName2.begin(), histName2.end(), '.', '_'); + + mHistArrQA.push_back(std::get>(histos.add(histName2.c_str(), Form("dEta dPhi for ptbin %d %s %s;#Delta#eta;#Delta#phi", i, charge[ch].data(), detadphiRange[1].data()), HistType::kTH2F, {{50, -0.02, 0.02}, {50, -0.05, 0.05}}))); + + for (Int_t j = 0; j < 4; ++j) { + std::string ptOrdSafe = ptOrd[j].data(); + std::replace(ptOrdSafe.begin(), ptOrdSafe.end(), '>', '_'); + std::replace(ptOrdSafe.begin(), ptOrdSafe.end(), '<', '_'); + std::replace(ptOrdSafe.begin(), ptOrdSafe.end(), '.', '_'); + std::string histNamePtOrd = Form("mDetaDphiPtOrd_%.2f-%.2f_%s_%s%s_%d", confPtBins.value[2 * i], confPtBins.value[2 * i + 1], ptOrdSafe.c_str(), charge[ch].data(), detadphiRange[1].data(), i); + std::replace(histNamePtOrd.begin(), histNamePtOrd.end(), '.', '_'); + histNamePtOrd += std::to_string(j); // Append unique index to avoid collisions + mHistArrQA.push_back(std::get>(histos1.add(histNamePtOrd.c_str(), Form("dEta dPhi for ptbin %d %s %s %s;#Delta#eta;#Delta#phi", i, ptOrdSafe.c_str(), charge[ch].data(), detadphiRange[1].data()), HistType::kTH2F, {{50, -0.02, 0.02}, {50, -0.05, 0.05}}))); + + std::replace(ptOrdSafe.begin(), ptOrdSafe.end(), '>', '_'); + std::replace(ptOrdSafe.begin(), ptOrdSafe.end(), '<', '_'); + std::replace(ptOrdSafe.begin(), ptOrdSafe.end(), '.', '_'); + + std::replace(histNamePtOrd.begin(), histNamePtOrd.end(), '.', '_'); + histNamePtOrd += std::to_string(j); // Append unique index to avoid collisions + mHistArrQA.push_back(std::get>(histos1.add(histNamePtOrd.c_str(), Form("dEta dPhi for ptbin %d %s %s %s;#Delta#eta;#Delta#phi", i, ptOrdSafe.c_str(), charge[ch].data(), detadphiRange[1].data()), HistType::kTH2F, {{50, -0.02, 0.02}, {50, -0.05, 0.05}}))); + } } + std::string histNameChDiff = Form("mDetaDphiChDiff_%.2f-%.2f%s", confPtBins.value[2 * i], confPtBins.value[2 * i + 1], detadphiRange[0].data()); + std::replace(histNameChDiff.begin(), histNameChDiff.end(), '.', '_'); + + histNameChDiff += "_" + std::to_string(i); // Append unique index to avoid collisions + + mHistArrQA.push_back(std::get>(histos1.add(histNameChDiff.c_str(), Form("dEta dPhi for ptbin %d different ch;#Delta#eta;#Delta#phi", i), HistType::kTH2F, {{35, -1.75, 1.75}, {50, -TMath::Pi(), TMath::Pi()}}))); + + std::string histNameChDiffClose = Form("mDetaDphiChDiffClose_%.2f-%.2f%s", confPtBins.value[2 * i], confPtBins.value[2 * i + 1], detadphiRange[1].data()); + std::replace(histNameChDiffClose.begin(), histNameChDiffClose.end(), '.', '_'); + + histNameChDiffClose += "_" + std::to_string(i); // Append unique index to avoid collisions + + mHistArrQA.push_back(std::get>(histos1.add(histNameChDiffClose.c_str(), Form("dEta dPhi for ptbin %d different ch close;#Delta#eta;#Delta#phi", i), HistType::kTH2F, {{50, -0.02, 0.02}, {50, -0.05, 0.05}}))); } } - template + int chargeFromPDG(int pdg) + { + if (pdg == 0) + return 0; + TDatabasePDG* pdgDB = TDatabasePDG::Instance(); + TParticlePDG* particle = pdgDB->GetParticle(pdg); + return particle ? int(particle->Charge() / 3.0) : 0; + } + template void checkpT(const T& track) { + // int gIndex = track.globalIndex(); + // if (usedIndices.find(gIndex) != usedIndices.end()) { + // return; // Already used this track + //} + // usedIndices.insert(gIndex); // Mark this track as used + for (auto iPt = 0; iPt < confNumPt; ++iPt) { if (track.pt() > confPtBins.value[2 * iPt] && track.pt() < confPtBins.value[2 * iPt + 1]) { float iphi = track.phi(); iphi = gRandom->Gaus(iphi, TMath::TwoPi()); + if (iphi < 0) { iphi += TMath::TwoPi(); } else if (iphi > TMath::TwoPi()) { iphi -= TMath::TwoPi(); } + mHistArrQA[iPt * 4]->Fill(track.eta()); mHistArrQA[iPt * 4 + 1]->Fill(track.pt()); - mHistArrQA[iPt * 4 + 2]->Fill(iphi); + mHistArrQA[iPt * 4 + 2]->Fill(track.phi()); countTracks[iPt]++; + for (auto iM = 0; iM < nBins; ++iM) { mHistArrReset[iPt * nBins + iM]->Fill(track.eta(), track.phi()); } } } } + void calculateMoments(std::vector> hist) { Double_t binContent = 0; + countSamples++; + Bool_t compSample = kFALSE; + if (countSamples == confSamplesize) { + compSample = kTRUE; + countSamples = 0; + } // Calculate the normalized factorial moments for (auto iPt = 0; iPt < confNumPt; ++iPt) { for (auto iM = 0; iM < nBins; ++iM) { @@ -171,30 +336,46 @@ struct FactorialMoments { double binconVal = 0; binconVal = hist[iPt * nBins + iM]->GetBinContent(iEta, iPhi); binContent += binconVal; - for (auto iOrder = 0; iOrder < 6; ++iOrder) { + for (auto iq = 0; iq < 6; ++iq) { Double_t fqBin = 0; - if (binconVal >= iOrder + 2) { - fqBin = TMath::Factorial(binconVal) / (TMath::Factorial(binconVal - (iOrder + 2))); + if (binconVal >= iq + 2) { + fqBin = TMath::Factorial(binconVal) / (TMath::Factorial(binconVal - (iq + 2))); } if (std::isnan(fqBin)) { break; } - sumfqBin[iOrder] += fqBin; + sumfqBin[iq] += fqBin; } } } binConEvent[iPt][iM] = binContent / (TMath::Power(binningM[iM], 2)); - for (auto iOrder = 0; iOrder < 6; ++iOrder) { - if (sumfqBin[iOrder] > 0) { - fqEvent[iOrder][iPt][iM] = sumfqBin[iOrder] / (TMath::Power(binningM[iM], 2)); + for (auto iq = 0; iq < 6; ++iq) { + if (sumfqBin[iq] > 0) { + fqEvent[iq][iPt][iM] = sumfqBin[iq] / (TMath::Power(binningM[iM], 2)); + fqEventSampled[iq][iPt][iM] += fqEvent[iq][iPt][iM]; } - mFqBinFinal[iPt * 6 + iOrder]->Fill(iM, fqEvent[iOrder][iPt][iM]); - mBinConFinal[iPt * 6 + iOrder]->Fill(iM, binConEvent[iPt][iM]); - } + binConEventSampled[iq][iPt][iM] += binConEvent[iPt][iM]; + mFqBinFinal[iPt * 6 + iq]->Fill(iM, fqEvent[iq][iPt][iM]); + mBinConFinal[iPt * 6 + iq]->Fill(iM, binConEvent[iPt][iM]); + if (compSample) { + mBinConFinalSampled[iPt * 6 + iq]->Fill(iM, binConEventSampled[iq][iPt][iM] / confSamplesize); + + double tmp = (fqEventSampled[iq][iPt][iM] / (confSamplesize)) / (pow(binConEventSampled[iq][iPt][iM] / (confSamplesize), (iq + 2))); + mFqBinFinalSampled[iPt * 6 + iq]->Fill(iM, tmp); + tmpFqErr[iq][iPt][iM]->Fill(tmp); + errorFq[iq][iPt][iM] += TMath::Power(fqEventSampled[iq][iPt][iM] / (confSamplesize), 2); + mFqError[iPt * 6 + iq]->SetBinContent(iM + 1, 0); + mFqError[iPt * 6 + iq]->Fill(iM, tmpFqErr[iq][iPt][iM]->GetStdDev()); + + fqEventSampled[iq][iPt][iM] = 0; + binConEventSampled[iq][iPt][iM] = 0; + } + } // end of loop over order of F } // end of loop over M bins } // end of loop over pT bins } + // write a template that takes coll and tracks from processRun3 and processMCRec using TracksFMs = soa::Filtered>; void processRun3(soa::Filtered>::iterator const& coll, TracksFMs const& tracks) { @@ -253,7 +434,6 @@ struct FactorialMoments { } } for (auto iPt = 0; iPt < confNumPt; ++iPt) { - // if (countTracks[iPt] > 0) if (countTracks[iPt] > 0) { mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); } @@ -261,7 +441,287 @@ struct FactorialMoments { // Calculate the normalized factorial moments calculateMoments(mHistArrReset); } - PROCESS_SWITCH(FactorialMoments, processRun3, "main process function", true); + PROCESS_SWITCH(FactorialMoments, processRun3, "main process function", false); + + using CollisionCandidateMCRec = soa::Join; + using TracksMc = soa::Filtered>; + void processMCRec(soa::Filtered::iterator const& coll, TracksMc const& colltracks, aod::McParticles const&, aod::McCollisions const&) + { + if (!coll.has_mcCollision()) { + return; + } + histos.fill(HIST("mEventSelected"), 0); + if (!coll.sel8()) { + return; + } + histos.fill(HIST("mEventSelected"), 1); + if (IsApplySameBunchPileup && !coll.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + return; + } + histos.fill(HIST("mEventSelected"), 2); + if (IsApplyGoodZvtxFT0vsPV && !coll.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + return; + } + histos.fill(HIST("mEventSelected"), 3); + if (IsApplyVertexITSTPC && !coll.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + return; + } + histos.fill(HIST("mEventSelected"), 4); + if (coll.centFT0C() < confCentCut.value[0] || coll.centFT0C() > confCentCut.value[1]) { + return; + } + histos.fill(HIST("mEventSelected"), 5); + histos.fill(HIST("mVertexX"), coll.posX()); + histos.fill(HIST("mVertexY"), coll.posY()); + histos.fill(HIST("mVertexZ"), coll.posZ()); + // histos.fill(HIST("mCentFT0M"), coll.centFT0M()); + // histos.fill(HIST("mCentFV0A"), coll.centFV0A()); + // histos.fill(HIST("mCentFT0A"), coll.centFT0A()); + histos.fill(HIST("mCentFT0C"), coll.centFT0C()); + for (auto const& h : mHistArrReset) { + h->Reset(); + } + countTracks = {0, 0, 0, 0, 0}; + fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; + binConEvent = {{{0, 0, 0, 0, 0}}}; + for (auto const& track : colltracks) { + histos.fill(HIST("mCollID"), track.collisionId()); + histos.fill(HIST("mEta"), track.eta()); + histos.fill(HIST("mPt"), track.pt()); + histos.fill(HIST("mPhi"), track.phi()); + histos.fill(HIST("mNFindableClsTPC"), track.tpcNClsFindable()); + histos.fill(HIST("mNClsTPC"), track.tpcNClsFound()); + histos.fill(HIST("mNClsITS"), track.itsNCls()); + histos.fill(HIST("mChi2TPC"), track.tpcChi2NCl()); + histos.fill(HIST("mChi2ITS"), track.itsChi2NCl()); + histos.fill(HIST("mChi2TRD"), track.trdChi2()); + histos.fill(HIST("mDCAxy"), track.dcaXY()); + histos.fill(HIST("mDCAx"), track.dcaZ()); + histos.fill(HIST("mDCAxyPt"), track.pt(), track.dcaXY()); + histos.fill(HIST("mDCAzPt"), track.pt(), track.dcaZ()); + histos.fill(HIST("mNSharedClsTPC"), track.tpcNClsShared()); + histos.fill(HIST("mCrossedRowsTPC"), track.tpcNClsCrossedRows()); + histos.fill(HIST("mNFinClsminusCRows"), track.tpcNClsFindableMinusCrossedRows()); + histos.fill(HIST("mNFractionShClsTPC"), track.tpcFractionSharedCls()); + histos.fill(HIST("mSharedClsvsPt"), track.pt(), track.tpcNClsShared()); + histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), track.tpcFractionSharedCls() / track.tpcNClsCrossedRows()); + checkpT(track); + } + for (auto iPt = 0; iPt < confNumPt; ++iPt) { + if (countTracks[iPt] > 0) { + mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); + } + } + histos.fill(HIST("mEventSelected"), 6); + // Calculate the normalized factorial moments + calculateMoments(mHistArrReset); + } + PROCESS_SWITCH(FactorialMoments, processMCRec, "main process function", false); + + void processMCgenr3(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles) + { + + // access MC truth information with mcCollision() and mcParticle() methods + + if (reduceOutput < 2) { + LOGF(info, "MC. vtx-z = %f", mcCollision.posZ()); + LOGF(info, "First: %d | Length: %d", mcParticles.begin().index(), mcParticles.size()); + } + int count = 0; + countTracks = {0, 0, 0, 0, 0}; + fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; + binConEvent = {{{0, 0, 0, 0, 0}}}; + if (mcCollision.impactParameter() >= 3.5) { + return; + } + for (auto& mcParticle : mcParticles) { + int pdgCode = mcParticle.pdgCode(); + auto* pd = TDatabasePDG::Instance()->GetParticle(pdgCode); + double charge = (pd != nullptr) ? pd->Charge() : 0.; + if (mcParticle.isPhysicalPrimary() && std::abs(mcParticle.eta()) < 0.8 && std::abs(charge) >= 1e-6) { + histos.fill(HIST("mEta"), mcParticle.eta()); + histos.fill(HIST("mPt"), mcParticle.pt()); + histos.fill(HIST("mPhi"), mcParticle.phi()); + count++; + // Loop over mothers and daughters + if (mcParticle.has_mothers()) { + // Check first mother + auto const& mother = mcParticle.mothers_first_as(); + if (reduceOutput == 0) { + LOGF(info, "First mother: %d has pdg code %d", mother.globalIndex(), mother.pdgCode()); + } + // Loop over all mothers (needed for some MCs with junctions etc.) + for (auto& m : mcParticle.mothers_as()) { + LOGF(debug, "M2 %d %d", mcParticle.globalIndex(), m.globalIndex()); + } + } + if (mcParticle.has_daughters()) { + for (auto& d : mcParticle.daughters_as()) { + LOGF(debug, "D2 %d %d", mcParticle.globalIndex(), d.globalIndex()); + } + } + checkpT(mcParticle); + } + histos1.fill(HIST("mPrimariesPerEvent"), count); + } + for (auto iPt = 0; iPt < confNumPt; ++iPt) { + if (countTracks[iPt] > 0) { + mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); + } + } + if (reduceOutput < 2) { + LOGF(info, "Primaries for this collision: %d", count); + } + calculateMoments(mHistArrReset); + } + PROCESS_SWITCH(FactorialMoments, processMCgenr3, "main process function", false); + using EventSelection_Run2 = soa::Join; + using TracksRecSim = soa::Join; + using CollisionRecSim_Run2 = soa::Filtered>::iterator; + using BCsWithRun2Info = soa::Join; + Preslice perMcCollision = aod::mcparticle::mcCollisionId; + + void processMcRun2(CollisionRecSim_Run2 const& coll, + aod::BCs const& bcs, + TracksRecSim const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + BCsWithRun2Info const& bcsInfo) + { + auto bc = coll.bc_as(); + if (!(bc.eventCuts() & BIT(aod::Run2EventCuts::kAliEventCutsAccepted))) { + return; // Skip event if it doesn't pass cuts + } + + if (coll.centRun2V0M() < confCentCut.value[0] || coll.centRun2V0M() > confCentCut.value[1]) { + return; + } + + /// Fill basic vertex and centrality histograms + histos.fill(HIST("mVertexX"), coll.posX()); + histos.fill(HIST("mVertexY"), coll.posY()); + histos.fill(HIST("mVertexZ"), coll.posZ()); + histos.fill(HIST("mCentFT0M"), coll.centRun2V0M()); + + /// Reset histograms and counters + for (auto const& h : mHistArrReset) { + h->Reset(); + } + + countTracks = {0, 0, 0, 0, 0}; + fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; + binConEvent = {{{0, 0, 0, 0, 0}}}; + + /// Loop over reconstructed tracks + for (auto const& track : tracks) { + double recoCharge = (track.sign() != 0) ? track.sign() : 0.; + if (std::abs(track.eta()) < 0.8 && track.isGlobalTrack() && std::abs(recoCharge) >= 1e-6) { + histos.fill(HIST("mCollID"), track.collisionId()); + + histos.fill(HIST("mNFindableClsTPC"), track.tpcNClsFindable()); + histos.fill(HIST("mNClsTPC"), track.tpcNClsFound()); + histos.fill(HIST("mNClsITS"), track.itsNCls()); + histos.fill(HIST("mChi2TPC"), track.tpcChi2NCl()); + histos.fill(HIST("mChi2ITS"), track.itsChi2NCl()); + histos.fill(HIST("mChi2TRD"), track.trdChi2()); + histos.fill(HIST("mNSharedClsTPC"), track.tpcNClsShared()); + histos.fill(HIST("mCrossedRowsTPC"), track.tpcNClsCrossedRows()); + histos.fill(HIST("mNFinClsminusCRows"), track.tpcNClsFindableMinusCrossedRows()); + histos.fill(HIST("mNFractionShClsTPC"), track.tpcFractionSharedCls()); + histos.fill(HIST("mSharedClsvsPt"), track.pt(), track.tpcNClsShared()); + histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), track.tpcFractionSharedCls() / track.tpcNClsCrossedRows()); + + checkpT(track); + } + } + + /// Access associated mcCollision from index + auto mccollision = coll.mcCollision_as(); + + /// Impact parameter can be used here if needed + float impactpar = mccollision.impactParameter(); + // Optionally: fillMCCollision(coll, mcParticles, impactpar); + + /// Now slice over selectedMCParticles using mcCollision + auto mcParts = mcParticles.sliceBy(perMcCollision, coll.mcCollision().globalIndex()); + + for (auto const& mc : mcParts) { + int pdgCode = mc.pdgCode(); + auto* pd = TDatabasePDG::Instance()->GetParticle(pdgCode); + double charge = (pd != nullptr) ? pd->Charge() : 0.; + if (mc.isPhysicalPrimary() && std::abs(mc.eta()) < 0.8 && std::abs(charge) >= 1e-6) { + histos.fill(HIST("mEta"), mc.eta()); + histos.fill(HIST("mPt"), mc.pt()); + histos.fill(HIST("mPhi"), mc.phi()); + + // Use charged MC particle here + // checkpT(mc); + } + } + + /// Apply cuts on tracks per pT bin + for (auto iPt = 0; iPt < confNumPt; ++iPt) { + if (countTracks[iPt] > 0) { + mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); + } else { + return; + } + } + + /// Final calculation of normalized factorial moments + calculateMoments(mHistArrReset); + } + + PROCESS_SWITCH(FactorialMoments, processMcRun2, "process MC Run2", true); + + void processMcGenRun2(aod::McCollision const& coll, aod::McParticles const& mcParticles) + { + auto bc = coll.bc_as(); + if (!(bc.eventCuts() & BIT(aod::Run2EventCuts::kAliEventCutsAccepted))) { + return; + } + histos1.fill(HIST("mCentMCImpactFull"), coll.impactParameter()); + // if (coll.impactParameter() >= 3.5) { + // return; + //} + histos1.fill(HIST("mCentMCImpact005"), coll.impactParameter()); + histos.fill(HIST("mVertexX"), coll.posX()); + histos.fill(HIST("mVertexY"), coll.posY()); + histos.fill(HIST("mVertexZ"), coll.posZ()); + + for (auto const& h : mHistArrReset) { + h->Reset(); + } + + countTracks = {0, 0, 0, 0, 0}; + fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; + binConEvent = {{{0, 0, 0, 0, 0}}}; + for (auto const& track : mcParticles) { + histos1.fill(HIST("hRecoPtBefore"), track.pt()); + histos1.fill(HIST("hRecoEtaBefore"), track.eta()); + histos1.fill(HIST("hRecoPhiBefore"), track.phi()); + if (std::abs(track.eta()) < 0.8 && track.isPhysicalPrimary()) { + histos1.fill(HIST("hRecoPtAfter"), track.pt()); + histos1.fill(HIST("hRecoEtaAfter"), track.eta()); + histos1.fill(HIST("hRecoPhiAfter"), track.phi()); + checkpT(track); + } + histos.fill(HIST("mTrackSelected"), 1); + } + + // counter++; + for (auto iPt = 0; iPt < confNumPt; ++iPt) { + if (countTracks[iPt] > 0) { + mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); + } else { + return; + } + } + + calculateMoments(mHistArrReset); + } + + PROCESS_SWITCH(FactorialMoments, processMcGenRun2, "process MC Run2", false); void processRun2(soa::Filtered>::iterator const& coll, TracksFMs const& tracks) { @@ -311,7 +771,7 @@ struct FactorialMoments { checkpT(track); } for (auto iPt = 0; iPt < confNumPt; ++iPt) { - if (countTracks[iPt] > 500) { + if (countTracks[iPt] > 0) { mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); } else { return;