diff --git a/PWGCF/Flow/Tasks/pidFlowPtCorr.cxx b/PWGCF/Flow/Tasks/pidFlowPtCorr.cxx index ea938dd55a8..cb0fe7b46d0 100644 --- a/PWGCF/Flow/Tasks/pidFlowPtCorr.cxx +++ b/PWGCF/Flow/Tasks/pidFlowPtCorr.cxx @@ -112,13 +112,20 @@ struct PidFlowPtCorr { } dcaCutOpts; O2_DEFINE_CONFIGURABLE(cfgNSigmapid, std::vector, (std::vector{3, 3, 3, 9, 9, 9, 9, 9, 9}), "tpc, tof and its NSigma for Pion Proton Kaon") - O2_DEFINE_CONFIGURABLE(cfgMeanPtcent, std::vector, (std::vector{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}), "mean Pt in different cent bin") struct : ConfigurableGroup { std::string prefix = "correctionPathOpts"; O2_DEFINE_CONFIGURABLE(cfgAcceptancePath, std::vector, (std::vector{"Users/f/fcui/NUA/NUAREFPartical", "Users/f/fcui/NUA/NUAK0s", "Users/f/fcui/NUA/NUALambda", "Users/f/fcui/NUA/NUAXi", "Users/f/fcui/NUA/NUAOmega"}), "CCDB path to acceptance object") + O2_DEFINE_CONFIGURABLE(cfgEfficiencyPath, std::vector, (std::vector{"PathtoRef"}), "CCDB path to efficiency object") O2_DEFINE_CONFIGURABLE(cfgEfficiency2DPath, std::vector, (std::vector{"PathtoRef"}), "CCDB path to efficiency(pt, cent) object") + O2_DEFINE_CONFIGURABLE(cfgPidEfficiencyPath, std::vector, (std::vector{"PathtoRef"}), "Pi, Ka, Pr, CCDB path to PID efficiency(pt, cent) object") + + O2_DEFINE_CONFIGURABLE(cfgEfficiencyPath4ITSOnly, std::vector, (std::vector{"PathtoRef"}), "CCDB path to efficiency object") + O2_DEFINE_CONFIGURABLE(cfgEfficiency2DPath4ITSOnly, std::vector, (std::vector{"PathtoRef"}), "CCDB path to efficiency(pt, cent) object") + O2_DEFINE_CONFIGURABLE(cfgPidEfficiencyPath4ITSOnly, std::vector, (std::vector{"PathtoRef"}), "Pi, Ka, Pr, CCDB path to PID efficiency(pt, cent) object") + + O2_DEFINE_CONFIGURABLE(cfgNUEOption, int, 1, "do NUE, 1: use 1D eff, 2: Using Eff(pt, cent), 3: use pid 1D eff, 4: use pid 2D eff, other: dont do NUE") } correctionPathOpts; O2_DEFINE_CONFIGURABLE(cfgRunNumbers, std::vector, (std::vector{544095, 544098, 544116, 544121, 544122, 544123, 544124}), "Preconfigured run numbers") @@ -126,15 +133,13 @@ struct PidFlowPtCorr { O2_DEFINE_CONFIGURABLE(cfgFlowNbootstrap, int, 30, "Number of subsamples for bootstrap") // switch - O2_DEFINE_CONFIGURABLE(cfgDoAccEffCorr, bool, false, "do acc and eff corr") O2_DEFINE_CONFIGURABLE(cfgDoLocDenCorr, bool, false, "do local density corr") O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights") - O2_DEFINE_CONFIGURABLE(cfgOutputrunbyrun, bool, false, "Fill and output NUA weights run by run") + O2_DEFINE_CONFIGURABLE(cfgOutputrunbyrun, bool, false, "OPEN IF USE FUNCTION(fillcorrectiongraph) Fill and output NUA weights run by run") O2_DEFINE_CONFIGURABLE(cfgOutPutMC, bool, false, "Fill MC graphs, note that if the processMCgen is open,this MUST be open") O2_DEFINE_CONFIGURABLE(cfgOutputLocDenWeights, bool, false, "Fill and output local density weights") - O2_DEFINE_CONFIGURABLE(cfgOutputQA, bool, false, "do QA") + O2_DEFINE_CONFIGURABLE(cfgOutputQA, bool, false, "OPEN IF USE FUNCTION(detectorPidQa) do QA") - O2_DEFINE_CONFIGURABLE(cfgUsePtCentNUECorr, bool, true, "do NUA NUE, Using Eff(pt, cent) to do NUE") O2_DEFINE_CONFIGURABLE(cfgDebugMyCode, bool, false, "output some graph for code debug") O2_DEFINE_CONFIGURABLE(cfgOutPutPtSpectra, bool, false, "output pt spectra for data, MC and RECO") @@ -184,7 +189,7 @@ struct PidFlowPtCorr { // filter and using // data Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; - Filter trackFilter = (nabs(aod::track::eta) < trkQualityOpts.cfgCutEta.value) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)); + Filter trackFilter = (nabs(aod::track::eta) < trkQualityOpts.cfgCutEta.value); using TracksPID = soa::Join; // data tracks filter @@ -226,12 +231,11 @@ struct PidFlowPtCorr { // define global variables GFW* fGFW = new GFW(); // GFW class used from main src std::vector corrconfigs; + std::vector cfgAcceptance; - std::vector cfgEfficiency; - std::vector cfgEfficiency2D; + std::vector cfgMultPVCutPara; std::vector cfgNSigma; - std::vector cfgMeanPt; std::vector runNumbers; std::map>> th1sList; std::map>> th3sList; @@ -263,9 +267,11 @@ struct PidFlowPtCorr { funcNumber }; + // graphs for NUE / NUA std::vector mAcceptance; - std::vector mEfficiency; - std::vector mEfficiency2D; + std::vector mEfficiency; + std::vector mEfficiency4ITSOnly; + bool correctionsLoaded = false; TF1* fMultPVCutLow = nullptr; @@ -316,10 +322,8 @@ struct PidFlowPtCorr { ccdb->setCreatedNotAfter(cfgnolaterthan.value); cfgAcceptance = correctionPathOpts.cfgAcceptancePath.value; - cfgEfficiency = correctionPathOpts.cfgEfficiencyPath.value; - cfgEfficiency2D = correctionPathOpts.cfgEfficiency2DPath.value; + cfgNSigma = cfgNSigmapid; - cfgMeanPt = cfgMeanPtcent; cfgMultPVCutPara = evtSeleOpts.cfgMultPVCut; // Set the pt, mult and phi Axis; @@ -353,6 +357,7 @@ struct PidFlowPtCorr { registry.add("hEtaPhiVtxzREF", "", {HistType::kTH3D, {cfgaxisPhi, cfgaxisEta, {20, -10, 10}}}); registry.add("hNTracksPVvsCentrality", "", {HistType::kTH2D, {{5000, 0, 5000}, axisMultiplicity}}); + registry.add("hNchUnCorrectedVSNchCorrected", "", {HistType::kTH2D, {cfgaxisNch, cfgaxisNch}}); runNumbers = cfgRunNumbers; // TPC vs TOF vs its, comparation graphs, check the PID performance in difference pt if (cfgOutputQA) { @@ -396,6 +401,8 @@ struct PidFlowPtCorr { // hist for NUE correction registry.add("correction/hPtCentMcRec", "", {HistType::kTH2D, {cfgaxisPt, axisMultiplicity}}); registry.add("correction/hPtCentMcGen", "", {HistType::kTH2D, {cfgaxisPt, axisMultiplicity}}); + // ITS only + registry.add("correction/hPtCentMcRec4ITSOnly", "", {HistType::kTH2D, {cfgaxisPt, axisMultiplicity}}); // hist for Pi eff registry.add("correction/hPtCentMcRecPi", "", {HistType::kTH2D, {cfgaxisPt, axisMultiplicity}}); @@ -604,6 +611,13 @@ struct PidFlowPtCorr { corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiKaN08 {3} poiKaP08 {-3}", "KaKa08gap22", kFALSE)); corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPrN08 {3} poiPrP08 {-3}", "PrPr08gap22", kFALSE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPiN {2} refP {-2}", "Pion0gap22a", kFALSE)); // 35 + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPiP {2} refN {-2}", "Pion0gap22b", kFALSE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiKaN {2} refP {-2}", "Kaon0gap22a", kFALSE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiKaP {2} refN {-2}", "Kaon0gap22b", kFALSE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPrN {2} refP {-2}", "Prot0gap22a", kFALSE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPrP {2} refN {-2}", "Prot0gap22b", kFALSE)); // 40 + fGFW->CreateRegions(); // finalize the initialization // used for event selection @@ -617,16 +631,10 @@ struct PidFlowPtCorr { fT0AV0ASigma = new TF1("fT0AV0ASigma", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 200000); fT0AV0ASigma->SetParameters(463.4144, 6.796509e-02, -9.097136e-07, 7.971088e-12, -2.600581e-17); - // fWeight output - if (cfgOutputNUAWeights) { - fWeightsREF->setPtBins(nPtBins, &(axisPt.binEdges)[0]); - fWeightsREF->init(true, false); - } - std::vector pTEffBins = {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0}; hFindPtBin = new TH1D("hFindPtBin", "hFindPtBin", pTEffBins.size() - 1, &pTEffBins[0]); funcEff.resize(pTEffBins.size() - 1); - // LHC24g3 Eff + std::vector f1p0 = cfgTrackDensityP0; std::vector f1p1 = cfgTrackDensityP1; for (uint ifunc = 0; ifunc < pTEffBins.size() - 1; ifunc++) { @@ -967,78 +975,217 @@ struct PidFlowPtCorr { } // end load NUA - // load NUE 1d - if (cfgEfficiency.size() == static_cast(nspecies)) { - for (int i = 0; i <= nspecies - 1; i++) { - mEfficiency.push_back(ccdb->getForTimeStamp(cfgEfficiency[i], timestamp)); - } - if (mEfficiency.size() == static_cast(nspecies)) - LOGF(info, "Loaded efficiency histogram"); - else - LOGF(fatal, "Could not load efficiency histogram"); - } - // end load NUE1d + /// @note dont modify load NUA!, just modify NUE + // here needs to load 2 NUE graph selected from 6 ccdb path + /// @note when using eff 2D, convert TH1* to (TH2D *) + // pid eff is not available now... + std::vector effPath1D = correctionPathOpts.cfgEfficiencyPath.value; + std::vector effPath2D = correctionPathOpts.cfgEfficiency2DPath.value; + std::vector effPathPid = correctionPathOpts.cfgPidEfficiencyPath.value; + + std::vector effPath1D4ITSOnly = correctionPathOpts.cfgEfficiencyPath4ITSOnly.value; + std::vector effPath2D4ITSOnly = correctionPathOpts.cfgEfficiency2DPath4ITSOnly.value; + std::vector effPathPid4ITSOnly = correctionPathOpts.cfgPidEfficiencyPath4ITSOnly.value; + + switch (correctionPathOpts.cfgNUEOption.value) { + case 1: // use 1D eff + // load 1d eff for global track + if (effPath1D.size() != static_cast(1)) { + LOGF(warning, "eff path 1d size != 1, skip eff 1d load"); + break; + } + mEfficiency.push_back(ccdb->getForTimeStamp(effPath1D[0], timestamp)); + if (mEfficiency.size() == static_cast(1)) { + LOGF(info, "Loaded efficiency histogram"); + } else { + LOGF(fatal, "Could not load efficiency histogram"); + } + // end load 1d eff for global track - // load NUE 2D - if (cfgEfficiency2D.size() == static_cast(nspecies)) { - for (int i = 0; i <= nspecies - 1; i++) { - mEfficiency2D.push_back(ccdb->getForTimeStamp(cfgEfficiency2D[i], timestamp)); - } - if (mEfficiency2D.size() == static_cast(nspecies)) - LOGF(info, "Loaded efficiency2D histogram"); - else - LOGF(fatal, "Could not load efficiency2D histogram"); - } - // end load NUE 2D + // load 1d eff for ITS only + if (effPath1D4ITSOnly.size() != static_cast(1)) { + LOGF(warning, "eff path for its 1d size != 1, skip its eff 1d load"); + break; + } + mEfficiency4ITSOnly.push_back(ccdb->getForTimeStamp(effPath1D4ITSOnly[0], timestamp)); + if (mEfficiency4ITSOnly.size() == static_cast(1)) { + LOGF(info, "Loaded ITS only efficiency histogram"); + } else { + LOGF(fatal, "Could not load ITS only efficiency histogram"); + } + // end load 1d eff for its only + + break; + // end use 1d eff + + case 2: // Use 2D eff + // load 2d eff for global track + if (effPath2D.size() != static_cast(1)) { + LOGF(warning, "eff path 2d size != 1, skip eff 2d load"); + break; + } + mEfficiency.push_back(ccdb->getForTimeStamp(effPath2D[0], timestamp)); + if (mEfficiency.size() == static_cast(1)) { + LOGF(info, "Loaded 2D efficiency histogram"); + } else { + LOGF(fatal, "Could not load 2D efficiency histogram"); + } + // end load 2d eff for global track + + // load 2d eff for ITS only + if (effPath2D4ITSOnly.size() != static_cast(1)) { + LOGF(warning, "eff path for its 2d size != 1, skip its eff 2d load"); + break; + } + mEfficiency4ITSOnly.push_back(ccdb->getForTimeStamp(effPath2D4ITSOnly[0], timestamp)); + if (mEfficiency4ITSOnly.size() == static_cast(1)) { + LOGF(info, "Loaded ITS only 2D efficiency histogram"); + } else { + LOGF(fatal, "Could not load ITS only 2D efficiency histogram"); + } + // end load 2d eff for its only + + break; + // end use 2d eff + + case 3: // use pid 1D eff + // load pid 1d eff for ITS + global track + if (effPathPid.size() != static_cast(3)) { + LOGF(warning, "eff path pid 1d size != 3, skip pid eff 1d load"); + break; + } + for (int i = 0; i < 3; i++) { + mEfficiency.push_back(ccdb->getForTimeStamp(effPathPid[i], timestamp)); + } + if (mEfficiency.size() == static_cast(3)) { + LOGF(info, "Loaded PID efficiency histogram"); + } else { + LOGF(fatal, "Could not load PID efficiency histogram"); + } + // end load pid 1d eff for ITS + global track + + // load pid 1d eff for ITS only + if (effPathPid4ITSOnly.size() != static_cast(3)) { + LOGF(warning, "eff path for its pid 1d size != 3, skip its pid eff 1d load"); + break; + } + for (int i = 0; i < 3; i++) { + mEfficiency4ITSOnly.push_back(ccdb->getForTimeStamp(effPathPid4ITSOnly[i], timestamp)); + } + if (mEfficiency4ITSOnly.size() == static_cast(3)) { + LOGF(info, "Loaded ITS only PID efficiency histogram"); + } else { + LOGF(fatal, "Could not load ITS only PID efficiency histogram"); + } + // end load pid 1d eff for its only + + break; + // end use pid 1d eff + + default: // not do NUE + break; + } // end switch load eff correctionsLoaded = true; } template - bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, TrackObject& track, float vtxz, int ispecies) + bool setParticleNUAWeight(float& weight_nua, TrackObject& track, float vtxz) { - float eff = 1.; - int nspecies = 1; - if (mEfficiency.size() == static_cast(nspecies)) - eff = mEfficiency[ispecies]->GetBinContent(mEfficiency[ispecies]->FindBin(track.pt())); - else - eff = 1.0; - if (eff == 0) { - return false; - } - weight_nue = 1. / eff; - if (mAcceptance.size() == static_cast(nspecies)) - weight_nua = mAcceptance[ispecies]->getNUA(track.phi(), track.eta(), vtxz); + if (mAcceptance.size() == static_cast(1)) + weight_nua = mAcceptance[0]->getNUA(track.phi(), track.eta(), vtxz); else weight_nua = 1; return true; } + /** + * @brief Set the Particle Nue Weight, for global track and ITS track, use different eff path + * + * @tparam TrackObject + * @param weight_nue + * @param track + * @param cent + * @param isGlobalTrack + * @return true weight setted succesfully, note that weight == 1 also return true; + * @return false eff got from graph == 0 + */ template - bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, TrackObject& track, float vtxz, int ispecies, double cent) + bool setParticleNUEWeight(float& weight_nue, TrackObject& track, double cent, bool isGlobalTrack) { float eff = 1.; - int nspecies = 1; - // eff 2d - if (mEfficiency2D.size() == static_cast(nspecies)) { - int ptBin = mEfficiency2D[ispecies]->GetXaxis()->FindBin(track.pt()); - int centBin = mEfficiency2D[ispecies]->GetYaxis()->FindBin(cent); - eff = mEfficiency2D[ispecies]->GetBinContent(ptBin, centBin); - } - if (eff == 0) { - return false; + uint64_t sizeOfEffVec = mEfficiency.size(); + uint64_t sizeOfEffVec4ITS = mEfficiency4ITSOnly.size(); + + TH2* eff2D = 0; + TH2* eff2D4ITS = 0; + + /// @note 1. size check 2. load eff 3. dividezero check + switch (correctionPathOpts.cfgNUEOption.value) { + case 1: // 1d + if (sizeOfEffVec != 1) + break; + if (sizeOfEffVec4ITS != 1) + break; + + if (isGlobalTrack) + eff = mEfficiency[0]->GetBinContent(mEfficiency[0]->FindBin(track.pt())); + else + eff = mEfficiency4ITSOnly[0]->GetBinContent(mEfficiency4ITSOnly[0]->FindBin(track.pt())); + + if (eff == 0.) + return false; + + break; + // end 1d + + case 2: // 2d + if (sizeOfEffVec != 1) + break; + if (sizeOfEffVec4ITS != 1) + break; + + eff2D = dynamic_cast(mEfficiency[0]); + eff2D4ITS = dynamic_cast(mEfficiency4ITSOnly[0]); + + if (!eff2D || !eff2D4ITS) { + LOGF(warning, "failed while converting TH1 to TH2! skip eff apply"); + break; + } + + if (isGlobalTrack) { + int ptBin = eff2D->GetXaxis()->FindBin(track.pt()); + int centBin = eff2D->GetYaxis()->FindBin(cent); + eff = eff2D->GetBinContent(ptBin, centBin); + } else { + int ptBin = eff2D4ITS->GetXaxis()->FindBin(track.pt()); + int centBin = eff2D4ITS->GetYaxis()->FindBin(cent); + eff = eff2D4ITS->GetBinContent(ptBin, centBin); + } + + if (eff == 0.) + return false; + + break; + // end 2d + + /// @todo add pid NUE eff + case 3: // pid + if (sizeOfEffVec != 3) + break; + if (sizeOfEffVec4ITS != 3) + break; + + break; + // end pid + + default: + break; } - weight_nue = 1. / eff; - // end eff 2d - // acc - if (mAcceptance.size() == static_cast(nspecies)) - weight_nua = mAcceptance[ispecies]->getNUA(track.phi(), track.eta(), vtxz); - else - weight_nua = 1; + weight_nue = 1. / eff; return true; - // end acc } /** @@ -1046,8 +1193,7 @@ struct PidFlowPtCorr { * @note include * 1. eta cut * 2. pt cut - * 3. primaryparticle check - * 4. stable partccle + * 3. stable partccle * * @tparam mcParticle * @param particle @@ -1065,9 +1211,6 @@ struct PidFlowPtCorr { return false; if (particle.pt() > trkQualityOpts.cfgCutPtMax.value) return false; - // primary particle - if (!particle.isPhysicalPrimary()) - return false; // stable particle if (!isStable(particle.pdgCode())) return false; @@ -1076,16 +1219,12 @@ struct PidFlowPtCorr { } /** - * @brief track selection + * @brief global track selection, cut for specified detector is not include * @note include: 1. dcaxy - * 2. its ncls - * 3. tpc cross row - * 4. tpc ncls - * 5. pt cut - * 6. eta - * 7. tpc chi2 - * 8. dcaz - * 9. its CHI2 + * 2. pt cut + * 3. eta + * 4. dcaz + * * * @tparam TTrack * @param track @@ -1093,33 +1232,45 @@ struct PidFlowPtCorr { * @return false reject track */ template - bool trackSelected(TTrack& track) + bool trackSelectedGlobal(TTrack& track) { // dca xy cut // note that z cut is in filter if (dcaCutOpts.cfgDCAxyNSigma && (std::fabs(track.dcaXY()) > dcaCutOpts.fPtDepDCAxy->Eval(track.pt()))) return false; - // its ncls cut - if (track.itsNCls() <= trkQualityOpts.cfgITSNCls.value) - return false; - // tpc crossedrows cut - if (track.tpcNClsCrossedRows() <= trkQualityOpts.cfgTPCCrossedRows.value) - return false; - // tpc ncls cut - if (track.tpcNClsFound() <= trkQualityOpts.cfgTPCNCls.value) - return false; + // pt if (!((track.pt() > trkQualityOpts.cfgCutPtMin.value) && (track.pt() < trkQualityOpts.cfgCutPtMax.value))) return false; + // eta if (std::fabs(track.eta()) > trkQualityOpts.cfgCutEta.value) return false; - // tpc chi2 - if (track.tpcChi2NCl() > cfgCutChi2prTPCcls) - return false; + // dca z if (std::fabs(track.dcaZ()) > dcaCutOpts.cfgCutDCAz.value) return false; + + return true; + } + + /** + * @brief cut for ITS, include: + * 1. its n cluster + * 2. its chi2 + * + * @tparam TTrack + * @param track + * @return true + * @return false + */ + template + bool trackSelected4ITS(TTrack& track) + { + // its ncls cut + if (track.itsNCls() <= trkQualityOpts.cfgITSNCls.value) + return false; + // its chi2 if (track.itsChi2NCl() > trkQualityOpts.cfgITSChi2NDF.value) return false; @@ -1127,6 +1278,33 @@ struct PidFlowPtCorr { return true; } + /** + * @brief cut for TPC, include + * 1. TPC cross rows + * 2. tpc n cluster + * 3. tpc chi2 + * + * @tparam TTrack + * @param track + * @return true + * @return false + */ + template + bool trackSelected4TPC(TTrack& track) + { + // tpc crossedrows cut + if (track.tpcNClsCrossedRows() <= trkQualityOpts.cfgTPCCrossedRows.value) + return false; + // tpc ncls cut + if (track.tpcNClsFound() <= trkQualityOpts.cfgTPCNCls.value) + return false; + // tpc chi2 + if (track.tpcChi2NCl() > cfgCutChi2prTPCcls) + return false; + + return true; + } + /** * @brief fill eventCount for different function * @@ -1335,26 +1513,71 @@ struct PidFlowPtCorr { v4 = funcV4->Eval(cent); } // cfgDoLocDenCorr - float weff = 1; - float wacc = 1; double ptSum = 0; double ptSumw2 = 0; double nch = 0; double nchSquare = 0; double ptSquareSum = 0; - // fill GFW ref flow + double nchCorrectedPassed = 0; + double nchPassedSelection = 0; + + if (cfgDebugMyCode) { + LOGF(info, "===================================="); + } + + /// @note calculate pt + /// use ITS only for (const auto& track : tracks) { - if (cfgDoAccEffCorr) { - if (cfgUsePtCentNUECorr) { - // true - if (!setCurrentParticleWeights(weff, wacc, track, vtxz, 0, cent)) - continue; - } else { - // false - if (!setCurrentParticleWeights(weff, wacc, track, vtxz, 0)) - continue; - } // cfgUsePtCentNUECorr - } // cfgDoAccEffCorr + float weff = 1; + + // do nue + setParticleNUEWeight(weff, track, cent, false); + // end do nue + + // track cut ITS only + if (!trackSelectedGlobal(track)) + continue; + if (!trackSelected4ITS(track)) + continue; + if (track.hasTPC()) + continue; + // end track cut its only + + // calculate ncharged(nch with weight) and pt + if (std::fabs(track.eta()) < trkQualityOpts.cfgRangeEta.value) { + nch += weff; + nchSquare += weff * weff; + ptSum += weff * track.pt(); + ptSumw2 += weff * weff * track.pt(); + ptSquareSum += weff * weff * track.pt() * track.pt(); + } + // end calculate nch and pt + + nchCorrectedPassed += weff; + nchPassedSelection += 1; + + if (weff == 1. && cfgDebugMyCode) { + LOGF(info, "weff is 1, if nue opt is open and this message appears a lot, check!"); + } + } + // end pt calculation using ITS only + + if (cfgDebugMyCode) { + LOGF(info, Form("its only track num %f", nchPassedSelection)); + } + + int totalGlobalTrack = 0; + + // fill v2 flow + // use global track + for (const auto& track : tracks) { + float weff = 1; + float wacc = 1; + + // do NUE && NUA + setParticleNUAWeight(wacc, track, vtxz); + setParticleNUEWeight(weff, track, cent, true); + // end do NUE && NUA if (cfgDoLocDenCorr) { bool withinPtRef = (trkQualityOpts.cfgCutPtMin.value < track.pt()) && (track.pt() < trkQualityOpts.cfgCutPtMax.value); @@ -1412,10 +1635,28 @@ struct PidFlowPtCorr { // end thn wacc graph } // cfgDebugMycode - // track cut - if (!trackSelected(track)) + // track cut, global + ITS + TPC + if (!trackSelectedGlobal(track)) + continue; + if (!track.hasITS()) + continue; + if (!track.hasTPC()) continue; + if (!trackSelected4ITS(track)) + continue; + if (!trackSelected4TPC(track)) + continue; + // end track cut + totalGlobalTrack++; + + if (cfgDebugMyCode && weff == 1.) { + LOGF(info, "weff for global track is 1, if NUE is open and this appears a lot, check!"); + } + + if (cfgDebugMyCode && wacc == 1.) { + LOGF(info, "wacc for global track is 1, if NUA is open and this appears alot, check!"); + } // fill QA hist registry.fill(HIST("hPhi"), track.phi()); @@ -1446,11 +1687,6 @@ struct PidFlowPtCorr { registry.fill(HIST("ptSpectra/hPtCentDataNegEtaNegCh"), track.pt(), cent); } - // track pt cut - if (!((track.pt() > trkQualityOpts.cfgCutPtMin.value) && (track.pt() < trkQualityOpts.cfgCutPtMax.value))) - continue; - // end track pt cut - // fill GFW // bit mask 1: fill CHARGED PARTICLES fGFW->Fill(track.eta(), 0, track.phi(), wacc * weff, 1); //(eta, ptbin, phi, wacc*weff, bitmask) @@ -1473,20 +1709,13 @@ struct PidFlowPtCorr { // fill PROTONS and overlap Protons } // end fill GFW + } // end track loop for v2 calculation - if (cfgOutputNUAWeights) - fWeightsREF->fill(track.phi(), track.eta(), vtxz, track.pt(), cent, 0); + if (cfgDebugMyCode) { + LOGF(info, Form("global track num %d", totalGlobalTrack)); + } - // calculate ncharged(nch with weight) and pt - if (std::fabs(track.eta()) < trkQualityOpts.cfgRangeEta.value) { - nch += weff; - nchSquare += weff * weff; - ptSum += weff * track.pt(); - ptSumw2 += weff * weff * track.pt(); - ptSquareSum += weff * weff * track.pt() * track.pt(); - } - // end calculate nch and pt - } // end track loop + registry.fill(HIST("hNchUnCorrectedVSNchCorrected"), nchPassedSelection, nchCorrectedPassed); // fill hist using fGFW if (nch > 0) { @@ -1496,6 +1725,13 @@ struct PidFlowPtCorr { fillFC(MyParticleType::kCharged, corrconfigs.at(3), cent, rndm, "c32"); fillFC(MyParticleType::kCharged, corrconfigs.at(4), cent, rndm, "c34"); + fillFC(MyParticleType::kPion, corrconfigs.at(35), cent, rndm, "c22Full"); + fillFC(MyParticleType::kPion, corrconfigs.at(36), cent, rndm, "c22Full"); + fillFC(MyParticleType::kKaon, corrconfigs.at(37), cent, rndm, "c22Full"); + fillFC(MyParticleType::kKaon, corrconfigs.at(38), cent, rndm, "c22Full"); + fillFC(MyParticleType::kProton, corrconfigs.at(39), cent, rndm, "c22Full"); + fillFC(MyParticleType::kProton, corrconfigs.at(40), cent, rndm, "c22Full"); + fillFC(MyParticleType::kPion, corrconfigs.at(5), cent, rndm, "c22"); fillFC(MyParticleType::kPion, corrconfigs.at(6), cent, rndm, "c22"); fillFC(MyParticleType::kKaon, corrconfigs.at(7), cent, rndm, "c22"); @@ -1626,7 +1862,15 @@ struct PidFlowPtCorr { // loop all the track for (const auto& track : tracks) { // track cut - if (!trackSelected(track)) + if (!trackSelectedGlobal(track)) + continue; + if (!track.hasITS()) + continue; + if (!track.hasTPC()) + continue; + if (!trackSelected4ITS(track)) + continue; + if (!trackSelected4TPC(track)) continue; // end track cut @@ -1699,7 +1943,15 @@ struct PidFlowPtCorr { // start filling graphs for (const auto& track : tracks) { // track cut - if (!trackSelected(track)) + if (!trackSelectedGlobal(track)) + continue; + if (!track.hasITS()) + continue; + if (!track.hasTPC()) + continue; + if (!trackSelected4ITS(track)) + continue; + if (!trackSelected4TPC(track)) continue; // end track cut @@ -1759,7 +2011,7 @@ struct PidFlowPtCorr { // loop tracks for (const auto& track : tracks) { // track cut - if (!trackSelected(track)) + if (!trackSelectedGlobal(track)) continue; // end track cut @@ -1768,20 +2020,30 @@ struct PidFlowPtCorr { auto mcParticle = track.mcParticle(); // fill graph if (particleSelected(mcParticle)) { - // graph for all particles - registry.fill(HIST("correction/hPtCentMcRec"), track.pt(), cent); - - // identify particle and fill graph - if (isPion(track)) { - registry.fill(HIST("correction/hPtCentMcRecPi"), track.pt(), cent); - } - if (isKaon(track)) { - registry.fill(HIST("correction/hPtCentMcRecKa"), track.pt(), cent); + /// @note global track, fill rec hist + if (track.hasITS() && track.hasTPC() && trackSelected4ITS(track) && trackSelected4TPC(track)) { + // graph for all particles + registry.fill(HIST("correction/hPtCentMcRec"), track.pt(), cent); + + // identify particle and fill graph + if (isPion(track)) { + registry.fill(HIST("correction/hPtCentMcRecPi"), track.pt(), cent); + } + if (isKaon(track)) { + registry.fill(HIST("correction/hPtCentMcRecKa"), track.pt(), cent); + } + if (isProton(track)) { + registry.fill(HIST("correction/hPtCentMcRecPr"), track.pt(), cent); + } + // end identify particle and fill graph } - if (isProton(track)) { - registry.fill(HIST("correction/hPtCentMcRecPr"), track.pt(), cent); + // end global track, fill rec hist + + /// @note ITS track. fill rec hist + if (track.hasITS() && !track.hasTPC() && trackSelected4ITS(track)) { + registry.fill(HIST("correction/hPtCentMcRec4ITSOnly"), track.pt(), cent); } - // end identify particle and fill graph + // end ITS track. fill rec hist } // end fill graph } @@ -1837,7 +2099,7 @@ struct PidFlowPtCorr { // loop mc particles for (const auto& mcParticle : mcParticles) { // fill graph - if (particleSelected(mcParticle)) { + if (particleSelected(mcParticle) && mcParticle.isPhysicalPrimary()) { // graph for all particles registry.fill(HIST("correction/hPtCentMcGen"), mcParticle.pt(), cent);