diff --git a/PWGCF/EbyEFluctuations/Tasks/CMakeLists.txt b/PWGCF/EbyEFluctuations/Tasks/CMakeLists.txt index c65d1705a93..287c9651abc 100644 --- a/PWGCF/EbyEFluctuations/Tasks/CMakeLists.txt +++ b/PWGCF/EbyEFluctuations/Tasks/CMakeLists.txt @@ -50,7 +50,7 @@ o2physics_add_dpl_workflow(v0pt-pp-task COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(factorial-moments - SOURCES FactorialMomentsTask.cxx + SOURCES factorialMomentsTask.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore COMPONENT_NAME Analysis) diff --git a/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx b/PWGCF/EbyEFluctuations/Tasks/factorialMomentsTask.cxx similarity index 76% rename from PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx rename to PWGCF/EbyEFluctuations/Tasks/factorialMomentsTask.cxx index b544208e91e..9e333c4623a 100644 --- a/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/factorialMomentsTask.cxx @@ -4,14 +4,15 @@ // // 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 +// granted to it by virtue of its status as an intergovernmental Organization // or submit itself to any jurisdiction. /// + +/// \file factorialMomentsTask.cxx /// \brief This task is for Normalized Factorial Moments Analysis: PRC 85, 044914 (2012), nucl-ex:1411.6083 -/// \author Salman Malik -/// \author Balwan Singh +/// \author Salman Malik, Balwan Singh + #include "TRandom.h" #include @@ -42,15 +43,15 @@ TH1D* tmpFqErr[6][5][52]; struct FactorialMomentsTask { Configurable applyCheckPtForRec{"applyCheckPtForRec", false, "Apply checkpT for reconstructed tracks"}; Configurable applyCheckPtForMC{"applyCheckPtForMC", true, "Apply checkpT for MC-generated tracks"}; - Configurable centralEta{"centralEta", 0.9, "eta limit for tracks"}; - Configurable numPt{"numPt", 5, "number of pT bins"}; - Configurable ptMin{"ptMin", 0.2f, "lower pT cut"}; - Configurable dcaXY{"dcaXY", 2.4f, "DCA xy cut"}; - Configurable dcaZ{"dcaZ", 2.0f, "DCA z cut"}; - Configurable minTPCCls{"minTPCCls", 70.0f, "minimum number of TPC clusters"}; - Configurable> centLimits{"centLimits", {0, 5}, "centrality min and max"}; - Configurable> vertexXYZ{"vertexXYZ", {0.3f, 0.4f, 10.0f}, "vertex cuts"}; - Configurable> ptCuts{"ptCuts", {0.2f, 2.0f}, "pT cuts"}; + Configurable centralEta{"centralEta", 0.9, "eta limit for tracks"}; + Configurable numPt{"numPt", 5, "number of pT bins"}; + Configurable ptMin{"ptMin", 0.2f, "lower pT cut"}; + Configurable dcaXY{"dcaXY", 2.4f, "DCA xy cut"}; + Configurable dcaZ{"dcaZ", 2.0f, "DCA z cut"}; + Configurable mintPCCls{"mintPCCls", 70.0f, "minimum number of TPC clusters"}; + Configurable> centLimits{"centLimits", {0, 5}, "centrality min and max"}; + Configurable> vertexXYZ{"vertexXYZ", {0.3f, 0.4f, 10.0f}, "vertex cuts"}; + Configurable> ptCuts{"ptCuts", {0.2f, 2.0f}, "pT cuts"}; Configurable isApplySameBunchPileup{"isApplySameBunchPileup", true, "Enable SameBunchPileup cut"}; Configurable isApplyGoodZvtxFT0vsPV{"isApplyGoodZvtxFT0vsPV", true, "Enable GoodZvtxFT0vsPV cut"}; Configurable isApplyVertexITSTPC{"isApplyVertexITSTPC", true, "Enable VertexITSTPC cut"}; @@ -112,7 +113,7 @@ struct FactorialMomentsTask { {"mVertexZ", "vertexZ", {HistType::kTH1F, {{100, -20, 20}}}}, {"mEta", "#eta", {HistType::kTH1F, {{1000, -2, 2}}}}, {"mPt", "#pt", {HistType::kTH1F, {{1000, -0.01, 50}}}}, - {"mPhi", "#phi", {HistType::kTH1F, {{100, 0, TMath::TwoPi()}}}}, + {"mPhi", "#phi", {HistType::kTH1F, {{100, 0, o2::constants::math::TwoPI}}}}, {"mEvents", "events", {HistType::kTH1D, {{5, -0.5, 4.5}}}}, {"mNFindableClsTPC", "findable TPC clusters;findable clusters", {HistType::kTH1F, {{100, 0, 200}}}}, {"mNClsTPC", "number of clusters TPC; nClusters TPC", {HistType::kTH1F, {{100, 0, 200}}}}, @@ -135,18 +136,18 @@ struct FactorialMomentsTask { }, OutputObjHandlingPolicy::AnalysisObject, true}; - static const Int_t nBins = 52; + static const int nBins = 52; double kMinCharge = 1e-6; - static const Int_t nfqOrder = 6; - Int_t countSamples = 0; - Int_t testc1 = 0, testc2 = 0, testc3 = 0; - std::array binningM; - std::array countTracks{0, 0, 0, 0, 0}; - std::array, 5>, 6> fqEvent; - std::array, 5>, 6> fqEventSampled; - std::array, 5> binConEvent; - std::array, 5>, 6> binConEventSampled; - std::array, 5>, 6> errorFq = {{{{{0, 0, 0, 0, 0}}}}}; + static const int nfqOrder = 6; + int countSamples = 0; + int testc1 = 0, testc2 = 0, testc3 = 0; + std::array binningM; + std::array countTracks{0, 0, 0, 0, 0}; + std::array, 5>, 6> fqEvent; + std::array, 5>, 6> fqEventSampled; + std::array, 5> binConEvent; + std::array, 5>, 6> binConEventSampled; + std::array, 5>, 6> errorFq = {{{{{0, 0, 0, 0, 0}}}}}; std::vector> mHistArrReset; std::vector> mHistArrQA; std::vector> mFqBinFinal; @@ -185,22 +186,22 @@ struct FactorialMomentsTask { if (useMC) { auto mMcTrackSelected = std::get>(histos.add("mMcTrackSelected", "mcTrackSelection", HistType::kTH1D, {{5, 0.5, 5.5}})); } - for (Int_t iM = 0; iM < nBins; ++iM) { + for (int iM = 0; iM < nBins; ++iM) { binningM[iM] = 2 * (iM + 2); } - for (Int_t iPt = 0; iPt < numPt; ++iPt) { + for (int iPt = 0; iPt < numPt; ++iPt) { mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mEta", iPt + 1), Form("#eta for bin %.2f-%.2f;#eta", ptCuts.value[2 * iPt], ptCuts.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", ptCuts.value[2 * iPt], ptCuts.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", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, 2 * TMath::Pi()}}))); + mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mPhi", iPt + 1), Form("#phi for bin %.2f-%.2f;#phi", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, o2::constants::math::TwoPI}}))); mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mMultiplicity", iPt + 1), Form("Multiplicity for bin %.2f-%.2f;Multiplicity", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, 15000}}))); - for (Int_t 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, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH2F, {{binningM[iM], -0.8, 0.8}, {binningM[iM], 0, 2 * TMath::Pi()}})); + for (int 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, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH2F, {{binningM[iM], -0.8, 0.8}, {binningM[iM], 0, o2::constants::math::TwoPI}})); mHistArrReset.push_back(mHistsR); - for (Int_t iq = 0; iq < nfqOrder; ++iq) { + for (int iq = 0; iq < nfqOrder; ++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 (Int_t i = 0; i < nfqOrder; ++i) { + for (int i = 0; i < nfqOrder; ++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, ptCuts.value[2 * iPt], ptCuts.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, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); @@ -218,10 +219,10 @@ struct FactorialMomentsTask { template void checkpT(const T& track) { - for (Int_t iPt = 0; iPt < numPt; ++iPt) { + for (int iPt = 0; iPt < numPt; ++iPt) { if (track.pt() > ptCuts.value[2 * iPt] && track.pt() < ptCuts.value[2 * iPt + 1]) { float iphi = track.phi(); - iphi = gRandom->Gaus(iphi, TMath::TwoPi()); + iphi = gRandom->Gaus(iphi, o2::constants::math::TwoPI); iphi = RecoDecay::constrainAngle(iphi); mHistArrQA[iPt * 4]->Fill(track.eta()); @@ -229,7 +230,7 @@ struct FactorialMomentsTask { mHistArrQA[iPt * 4 + 2]->Fill(track.phi()); countTracks[iPt]++; - for (Int_t iM = 0; iM < nBins; ++iM) { + for (int iM = 0; iM < nBins; ++iM) { mHistArrReset[iPt * nBins + iM]->Fill(track.eta(), track.phi()); } } @@ -238,26 +239,26 @@ struct FactorialMomentsTask { void calculateMoments(std::vector> hist) { - Double_t binContent = 0; + double binContent = 0; countSamples++; - Bool_t compSample = kFALSE; + bool compSample = kFALSE; if (countSamples == samplesize) { compSample = kTRUE; countSamples = 0; } // Calculate the normalized factorial moments - for (Int_t iPt = 0; iPt < numPt; ++iPt) { - for (Int_t iM = 0; iM < nBins; ++iM) { + for (int iPt = 0; iPt < numPt; ++iPt) { + for (int iM = 0; iM < nBins; ++iM) { binContent = 0; - Double_t sumfqBin[6] = {0}; + double sumfqBin[6] = {0}; - for (Int_t iEta = 1; iEta <= hist[iPt * nBins + iM]->GetNbinsX(); ++iEta) { - for (Int_t iPhi = 1; iPhi <= hist[iPt * nBins + iM]->GetNbinsY(); ++iPhi) { + for (int iEta = 1; iEta <= hist[iPt * nBins + iM]->GetNbinsX(); ++iEta) { + for (int iPhi = 1; iPhi <= hist[iPt * nBins + iM]->GetNbinsY(); ++iPhi) { double binconVal = 0; binconVal = hist[iPt * nBins + iM]->GetBinContent(iEta, iPhi); binContent += binconVal; - for (Int_t iq = 0; iq < nfqOrder; ++iq) { - Double_t fqBin = 0; + for (int iq = 0; iq < nfqOrder; ++iq) { + double fqBin = 0; if (binconVal >= iq + 2) { fqBin = TMath::Factorial(binconVal) / (TMath::Factorial(binconVal - (iq + 2))); } @@ -268,10 +269,10 @@ struct FactorialMomentsTask { } } } - binConEvent[iPt][iM] = binContent / (TMath::Power(binningM[iM], 2)); - for (Int_t iq = 0; iq < nfqOrder; ++iq) { + binConEvent[iPt][iM] = binContent / (std::pow(binningM[iM], 2)); + for (int iq = 0; iq < nfqOrder; ++iq) { if (sumfqBin[iq] > 0) { - fqEvent[iq][iPt][iM] = sumfqBin[iq] / (TMath::Power(binningM[iM], 2)); + fqEvent[iq][iPt][iM] = sumfqBin[iq] / (std::pow(binningM[iM], 2)); fqEventSampled[iq][iPt][iM] += fqEvent[iq][iPt][iM]; } binConEventSampled[iq][iPt][iM] += binConEvent[iPt][iM]; @@ -280,10 +281,10 @@ struct FactorialMomentsTask { if (compSample) { mBinConFinalSampled[iPt * 6 + iq]->Fill(iM, binConEventSampled[iq][iPt][iM] / samplesize); - double tmp = (fqEventSampled[iq][iPt][iM] / (samplesize)) / (TMath::Power(binConEventSampled[iq][iPt][iM] / (samplesize), (iq + 2))); + double tmp = (fqEventSampled[iq][iPt][iM] / (samplesize)) / (std::pow(binConEventSampled[iq][iPt][iM] / (samplesize), (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] / (samplesize), 2); + errorFq[iq][iPt][iM] += std::pow(fqEventSampled[iq][iPt][iM] / (samplesize), 2); mFqError[iPt * 6 + iq]->SetBinContent(iM + 1, 0); mFqError[iPt * 6 + iq]->Fill(iM, tmpFqErr[iq][iPt][iM]->GetStdDev()); @@ -297,6 +298,7 @@ struct FactorialMomentsTask { // write a template that takes coll and tracks from processRun3 and processMCRec using TracksFMs = soa::Filtered>; + Preslice perMcCollision = aod::mcparticle::mcCollisionId; void processRun3(soa::Filtered>::iterator const& coll, TracksFMs const& tracks) { // selection of events @@ -353,7 +355,7 @@ struct FactorialMomentsTask { checkpT(track); } } - for (Int_t iPt = 0; iPt < numPt; ++iPt) { + for (int iPt = 0; iPt < numPt; ++iPt) { if (countTracks[iPt] > 0) { mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); } @@ -361,14 +363,112 @@ struct FactorialMomentsTask { calculateMoments(mHistArrReset); } PROCESS_SWITCH(FactorialMomentsTask, processRun3, "main process function", false); - using EventSelection_Run2 = soa::Join; + using CollisionCandidateMCRec = soa::Join; + using TracksMc = soa::Filtered>; + void processMCRec(soa::Filtered::iterator const& coll, TracksMc const& colltracks, aod::McParticles const& mcParticles, 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() < centLimits.value[0] || coll.centFT0C() > centLimits.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("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) { + // if (track.hasITS()) + // if (track.hasTPC()) + // if (track.isGlobalTrack()) { + 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); + //} + } + auto mcParts = mcParticles.sliceBy(perMcCollision, coll.mcCollision().globalIndex()); + for (auto const& mc : mcParts) { + int pdgCode = mc.pdgCode(); + auto pdgInfo = pdg->GetParticle(pdgCode); + if (!pdgInfo) { + continue; + } + double charge = pdgInfo->Charge(); + double physCharge = charge / 3.0; + histos.fill(HIST("mChargeBefore"), physCharge); + if (mc.isPhysicalPrimary() && std::abs(mc.eta()) < centralEta && std::abs(physCharge) >= kMinCharge) { + histos.fill(HIST("mChargeAfter"), physCharge); + histos.fill(HIST("mEta"), mc.eta()); + histos.fill(HIST("mPt"), mc.pt()); + histos.fill(HIST("mPhi"), mc.phi()); + if (applyCheckPtForMC && !applyCheckPtForRec) { + checkpT(mc); + } + } + } + for (auto iPt = 0; iPt < numPt; ++iPt) { + // if (countTracks[iPt] > 0)countTracks = {0, 0, 0, 0, 0}; + 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(FactorialMomentsTask, processMCRec, "main process function", false); + using EventSelectionrun2 = soa::Join; using TracksRecSim = soa::Join; - using CollisionRecSim_Run2 = soa::Filtered>::iterator; + using CollisionRecSimRun2 = soa::Filtered>::iterator; using BCsWithRun2Info = soa::Join; - Preslice perMcCollision = aod::mcparticle::mcCollisionId; - void processMcRun2(CollisionRecSim_Run2 const& coll, + void processMcRun2(CollisionRecSimRun2 const& coll, + aod::BCs const&, TracksRecSim const& tracks, - aod::McParticles const& mcParticles) + aod::McParticles const& mcParticles, + aod::McCollisions const&, + BCsWithRun2Info const&) { auto bc = coll.bc_as(); if (!(bc.eventCuts() & BIT(aod::Run2EventCuts::kAliEventCutsAccepted))) { @@ -429,7 +529,7 @@ struct FactorialMomentsTask { } } - for (Int_t iPt = 0; iPt < numPt; ++iPt) { + for (int iPt = 0; iPt < numPt; ++iPt) { if (countTracks[iPt] > 0) { mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); } else { @@ -461,7 +561,7 @@ struct FactorialMomentsTask { fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; binConEvent = {{{0, 0, 0, 0, 0}}}; for (auto const& track : tracks) { - if ((track.pt() < ptMin) || (!track.isGlobalTrack()) || (track.tpcNClsFindable() < minTPCCls)) { + if ((track.pt() < ptMin) || (!track.isGlobalTrack()) || (track.tpcNClsFindable() < mintPCCls)) { continue; } histos.fill(HIST("mCollID"), track.collisionId()); @@ -486,7 +586,7 @@ struct FactorialMomentsTask { histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), track.tpcFractionSharedCls() / track.tpcNClsCrossedRows()); checkpT(track); } - for (Int_t iPt = 0; iPt < numPt; ++iPt) { + for (int iPt = 0; iPt < numPt; ++iPt) { if (countTracks[iPt] > 0) { mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); } else {