diff --git a/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx b/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx index db3a65c374e..465b0caa8a3 100644 --- a/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx +++ b/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx @@ -115,7 +115,6 @@ struct UccZdc { Configurable isOccupancyCut{"isOccupancyCut", true, "Occupancy cut?"}; Configurable isApplyFT0CbasedOccupancy{"isApplyFT0CbasedOccupancy", false, "T0C Occu cut"}; Configurable isTDCcut{"isTDCcut", false, "Use TDC cut"}; - Configurable isZEMcut{"isZEMcut", true, "Use ZEM cut"}; Configurable useMidRapNchSel{"useMidRapNchSel", true, "Use mid-rapidit Nch selection"}; Configurable applyEff{"applyEff", true, "Apply track-by-track efficiency correction"}; Configurable applyFD{"applyFD", false, "Apply track-by-track feed down correction"}; @@ -147,6 +146,7 @@ struct UccZdc { Configurable minNch{"minNch", 0, "Min Nch (|eta|<0.8)"}; Configurable minZN{"minZN", -0.5, "Min ZN signal"}; Configurable minTdc{"minTdc", -15.0, "minimum TDC"}; + Configurable arbScale{"arbScale", 100.0, "Scale factor for forward multiplicity"}; Configurable maxNch{"maxNch", 3000, "Max Nch (|eta|<0.8)"}; Configurable maxITSTrack{"maxITSTrack", 6000., "Min ITS tracks"}; @@ -498,8 +498,8 @@ struct UccZdc { if (occuValue < minOccCut || occuValue > maxOccCut) { return false; } + registry.fill(HIST("hEventCounter"), EvCutLabel::OccuCut); } - registry.fill(HIST("hEventCounter"), EvCutLabel::OccuCut); if (col.centFT0C() < minT0CcentCut || col.centFT0C() > maxT0CcentCut) { return false; @@ -530,8 +530,8 @@ struct UccZdc { registry.fill(HIST("hEventCounter"), EvCutLabel::Zdc); auto zdc = foundBC.zdc(); - double aT0A{-999.}; - double aT0C{-999.}; + double aT0A{0.}; + double aT0C{0.}; if (foundBC.has_ft0()) { for (const auto& amplitude : foundBC.ft0().amplitudeA()) { aT0A += amplitude; @@ -544,7 +544,7 @@ struct UccZdc { } registry.fill(HIST("hEventCounter"), EvCutLabel::TZero); - double aV0A{-999.}; + double aV0A{0.}; if (foundBC.has_fv0a()) { for (const auto& amplitude : foundBC.fv0a().amplitude()) { aV0A += amplitude; @@ -552,10 +552,10 @@ struct UccZdc { } const double nPV{collision.multNTracksPVeta1() / 1.}; - const double normT0M{(aT0A + aT0C) / 100.}; - const double normV0A{aV0A / 100.}; - const double normT0A{aT0A / 100.}; - const double normT0C{aT0C / 100.}; + const double normT0M{(aT0A + aT0C) / arbScale}; + const double normV0A{aV0A / arbScale}; + const double normT0A{aT0A / arbScale}; + const double normT0C{aT0C / arbScale}; float znA{zdc.amplitudeZNA()}; float znC{zdc.amplitudeZNC()}; float zpA{zdc.amplitudeZPA()}; @@ -582,14 +582,6 @@ struct UccZdc { registry.fill(HIST("hEventCounter"), EvCutLabel::Tdc); } - // ZEM cut - if (isZEMcut) { - if (sumZEMs < zemCut) { - return; - } - registry.fill(HIST("hEventCounter"), EvCutLabel::Zem); - } - int itsTracks = 0, glbTracks = 0; for (const auto& track : tracks) { if (track.hasITS() && ((track.eta() > minEta) && (track.eta() < maxEta))) { @@ -667,19 +659,7 @@ struct UccZdc { registry.fill(HIST("zPos"), collision.posZ()); registry.fill(HIST("T0Ccent"), collision.centFT0C()); - registry.fill(HIST("ZNAamp"), znA); - registry.fill(HIST("ZNCamp"), znC); - registry.fill(HIST("ZPAamp"), zpA); - registry.fill(HIST("ZPCamp"), zpC); - registry.fill(HIST("ZNAVsZNC"), znC, znA); - registry.fill(HIST("ZNAVsZPA"), zpA, znA); - registry.fill(HIST("ZNCVsZPC"), zpC, znC); - registry.fill(HIST("ZPAVsZPC"), zpC, zpA); - registry.fill(HIST("ZNVsZEM"), sumZEMs, sumZNs); registry.fill(HIST("Debunch"), tZDCdif, tZDCsum); - registry.fill(HIST("ZNVsFT0A"), normT0A, sumZNs); - registry.fill(HIST("ZNVsFT0C"), normT0C, sumZNs); - registry.fill(HIST("ZNVsFT0M"), normT0M, sumZNs); registry.fill(HIST("NchVsFV0A"), normV0A, glbTracks); registry.fill(HIST("NchVsFT0A"), normT0A, glbTracks); registry.fill(HIST("NchVsFT0C"), normT0C, glbTracks); @@ -688,12 +668,29 @@ struct UccZdc { registry.fill(HIST("Nch"), glbTracks); registry.fill(HIST("NchVsNPV"), collision.multNTracksPVeta1(), glbTracks); registry.fill(HIST("NchVsITStracks"), itsTracks, glbTracks); - registry.fill(HIST("ZNAVsNch"), glbTracks, znA); - registry.fill(HIST("ZNCVsNch"), glbTracks, znC); - registry.fill(HIST("ZNVsNch"), glbTracks, sumZNs); - registry.fill(HIST("ZNDifVsNch"), glbTracks, znA - znC); if (glbTracks >= minNchSel) registry.fill(HIST("NchVsOneParCorr"), glbTracks, sumpt / glbTracks); + + // ZEM cut + if (sumZEMs > zemCut) { + registry.fill(HIST("hEventCounter"), EvCutLabel::Zem); + registry.fill(HIST("ZNAamp"), znA); + registry.fill(HIST("ZNCamp"), znC); + registry.fill(HIST("ZPAamp"), zpA); + registry.fill(HIST("ZPCamp"), zpC); + registry.fill(HIST("ZNAVsZNC"), znC, znA); + registry.fill(HIST("ZNAVsZPA"), zpA, znA); + registry.fill(HIST("ZNCVsZPC"), zpC, znC); + registry.fill(HIST("ZPAVsZPC"), zpC, zpA); + registry.fill(HIST("ZNVsZEM"), sumZEMs, sumZNs); + registry.fill(HIST("ZNVsFT0A"), normT0A, sumZNs); + registry.fill(HIST("ZNVsFT0C"), normT0C, sumZNs); + registry.fill(HIST("ZNVsFT0M"), normT0M, sumZNs); + registry.fill(HIST("ZNAVsNch"), glbTracks, znA); + registry.fill(HIST("ZNCVsNch"), glbTracks, znC); + registry.fill(HIST("ZNVsNch"), glbTracks, sumZNs); + registry.fill(HIST("ZNDifVsNch"), glbTracks, znA - znC); + } } PROCESS_SWITCH(UccZdc, processQA, "Process QA", true); void processZdcCollAss(o2::aod::ColEvSels::iterator const& collision, o2::aod::BCsRun3 const& /*bcs*/, aod::Zdcs const& /*zdcs*/, aod::FV0As const& /*fv0as*/, aod::FT0s const& /*ft0s*/, TheFilteredTracks const& tracks) @@ -711,8 +708,8 @@ struct UccZdc { } registry.fill(HIST("hEventCounter"), EvCutLabel::Zdc); - double aT0A{-999.}; - double aT0C{-999.}; + double aT0A{0.}; + double aT0C{0.}; if (foundBC.has_ft0()) { for (const auto& amplitude : foundBC.ft0().amplitudeA()) { aT0A += amplitude; @@ -725,7 +722,7 @@ struct UccZdc { } registry.fill(HIST("hEventCounter"), EvCutLabel::TZero); - double aV0A{-999.}; + double aV0A{0.}; if (foundBC.has_fv0a()) { for (const auto& amplitude : foundBC.fv0a().amplitude()) { aV0A += amplitude; @@ -733,8 +730,8 @@ struct UccZdc { } const double nPV{collision.multNTracksPVeta1() / 1.}; - const double normT0M{(aT0A + aT0C) / 100.}; - const double normV0A{aV0A / 100.}; + const double normT0M{(aT0A + aT0C) / arbScale}; + const double normV0A{aV0A / arbScale}; float znA{foundBC.zdc().amplitudeZNA()}; float znC{foundBC.zdc().amplitudeZNC()}; float zpA{foundBC.zdc().amplitudeZPA()}; @@ -772,14 +769,6 @@ struct UccZdc { registry.fill(HIST("hEventCounter"), EvCutLabel::Tdc); } - // ZEM cut - if (isZEMcut) { - if (sumZEMs < zemCut) { - return; - } - registry.fill(HIST("hEventCounter"), EvCutLabel::Zem); - } - // Nch-based selection double glbTracks{0.0}; for (const auto& track : tracks) { @@ -961,17 +950,11 @@ struct UccZdc { registry.fill(HIST("NchUncorrected"), glbTracks); registry.fill(HIST("NchVsV0A"), nchMult, normV0A); registry.fill(HIST("NchVsT0M"), nchMult, normT0M); - registry.fill(HIST("NchVsZN"), nchMult, sumZNs); - registry.fill(HIST("NchVsZP"), nchMult, sumZPs); registry.fill(HIST("NchVsOneParCorr"), nchMult, oneParCorr); registry.fill(HIST("NchVsTwoParCorr"), nchMult, twoParCorr); registry.fill(HIST("NchVsThreeParCorr"), nchMult, threeParCorr); - registry.fill(HIST("NchVsOneParCorrVsZN"), nchMult, sumZNs, oneParCorr); - registry.fill(HIST("NchVsTwoParCorrVsZN"), nchMult, sumZNs, twoParCorr); - registry.fill(HIST("NchVsThreeParCorrVsZN"), nchMult, sumZNs, threeParCorr); - registry.fill(HIST("NchVsOneParCorrVsT0M"), nchMult, normT0M, oneParCorr); registry.fill(HIST("NchVsTwoParCorrVsT0M"), nchMult, normT0M, twoParCorr); registry.fill(HIST("NchVsThreeParCorrVsT0M"), nchMult, normT0M, threeParCorr); @@ -980,8 +963,17 @@ struct UccZdc { registry.fill(HIST("NchVsTwoParCorrVsV0A"), nchMult, normV0A, twoParCorr); registry.fill(HIST("NchVsThreeParCorrVsV0A"), nchMult, normV0A, threeParCorr); + if (sumZEMs > zemCut) { + registry.fill(HIST("hEventCounter"), EvCutLabel::Zem); + registry.fill(HIST("NchVsZN"), nchMult, sumZNs); + registry.fill(HIST("NchVsZP"), nchMult, sumZPs); + registry.fill(HIST("NchVsOneParCorrVsZN"), nchMult, sumZNs, oneParCorr); + registry.fill(HIST("NchVsTwoParCorrVsZN"), nchMult, sumZNs, twoParCorr); + registry.fill(HIST("NchVsThreeParCorrVsZN"), nchMult, sumZNs, threeParCorr); + } + const uint64_t timeStamp{foundBC.timestamp()}; - eventSampling(tracks, normV0A, normT0M, sumZNs, timeStamp); + eventSampling(tracks, normV0A, normT0M, sumZNs, sumZEMs, timeStamp); } PROCESS_SWITCH(UccZdc, processZdcCollAss, "Process ZDC W/Coll Ass.", true); @@ -1021,17 +1013,10 @@ struct UccZdc { aT0C += amplitude; } } else { - return; + continue; } - // double aV0A{-999.}; - // if (foundBC.has_fv0a()) { - // for (const auto& amplitude : foundBC.fv0a().amplitude()) { aV0A += amplitude; } - // } - - const double normT0M{(aT0A + aT0C) / 100.}; - // const double normV0A{aV0A/100.}; - + const double normT0M{(aT0A + aT0C) / arbScale}; double nchRaw{0.}; double nchMult{0.}; double nchMC{0.}; @@ -1083,18 +1068,18 @@ struct UccZdc { // Reject event if nchRaw less than a lower cutoff if (nchRaw < minNchSel) { - return; + continue; } - // Calculates the event weight, W_k const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(nchRaw)}; + // Calculates the event weight, W_k for (const auto& track : groupedTracks) { // Track Selection if (track.eta() < minEta || track.eta() > maxEta) { continue; } - if (track.pt() < minPt || track.pt() > maxPt) { + if (track.pt() < minPt || track.pt() > maxPtSpectra) { continue; } if (!track.isGlobalTrack()) { @@ -1118,8 +1103,6 @@ struct UccZdc { if (std::abs(charge) < kMinCharge) { continue; } - // Is it a primary particle? - // if (!particle.isPhysicalPrimary()) { continue; } const double pt{static_cast(track.pt())}; const int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; @@ -1160,7 +1143,7 @@ struct UccZdc { std::vector vecFullEff; std::vector vecFDEqualOne; - // Calculates the event weight, W_k + // calculates the true Nch for (const auto& particle : mcParticles) { if (particle.eta() < minEta || particle.eta() > maxEta) { continue; @@ -1186,18 +1169,45 @@ struct UccZdc { if (!particle.isPhysicalPrimary()) { continue; } - - float pt{particle.pt()}; - pTsMC.emplace_back(pt); - vecFullEff.emplace_back(1.); - vecFDEqualOne.emplace_back(1.); nchMC++; } if (nchMC < minNchSel) { continue; } - // printf("nchMult = %f | nchMC = %f | nchMult/nchMc = %f\n",nchMult,nchMC,nchMult/nchMC); + + // Calculates the event weight, W_k + for (const auto& particle : mcParticles) { + if (particle.eta() < minEta || particle.eta() > maxEta) { + continue; + } + if (particle.pt() < minPt || particle.pt() > maxPtSpectra) { + continue; + } + + auto charge{0.}; + // Get the MC particle + auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle != nullptr) { + charge = pdgParticle->Charge(); + } else { + continue; + } + + // Is it a charged particle? + if (std::abs(charge) < kMinCharge) { + continue; + } + // Is it a primary particle? + if (!particle.isPhysicalPrimary()) { + continue; + } + + const float pt{particle.pt()}; + pTsMC.emplace_back(pt); + vecFullEff.emplace_back(1.); + vecFDEqualOne.emplace_back(1.); + } double p1MC, p2MC, p3MC, p4MC, w1MC, w2MC, w3MC, w4MC; p1MC = p2MC = p3MC = p4MC = w1MC = w2MC = w3MC = w4MC = 0.0; @@ -1223,7 +1233,6 @@ struct UccZdc { } else { // Correction with the remaining half of the sample registry.fill(HIST("EvtsDivided"), 1); //----- MC reconstructed -----// - // const auto& groupedTracks{simTracks.sliceBy(perCollision, collision.globalIndex())}; for (const auto& track : groupedTracks) { // Track Selection if (track.eta() < minEta || track.eta() > maxEta) { @@ -1235,9 +1244,6 @@ struct UccZdc { if (!track.isGlobalTrack()) { continue; } - registry.fill(HIST("ZposVsEta"), collision.posZ(), track.eta()); - registry.fill(HIST("EtaVsPhi"), track.eta(), track.phi()); - registry.fill(HIST("dcaXYvspT"), track.dcaXY(), track.pt()); nchRaw++; } @@ -1246,7 +1252,7 @@ struct UccZdc { if (track.eta() < minEta || track.eta() > maxEta) { continue; } - if (track.pt() < minPt || track.pt() > maxPt) { + if (track.pt() < minPt || track.pt() > maxPtSpectra) { continue; } if (!track.isGlobalTrack()) { @@ -1270,8 +1276,13 @@ struct UccZdc { if (std::abs(charge) < kMinCharge) { continue; } + // All charged particles registry.fill(HIST("Pt_all_ch"), nchRaw, track.pt()); + registry.fill(HIST("ZposVsEta"), collision.posZ(), track.eta()); + registry.fill(HIST("EtaVsPhi"), track.eta(), track.phi()); + registry.fill(HIST("dcaXYvspT"), track.dcaXY(), track.pt()); + // Is it a primary particle? if (!particle.isPhysicalPrimary()) { continue; @@ -1298,7 +1309,7 @@ struct UccZdc { if (particle.eta() < minEta || particle.eta() > maxEta) { continue; } - if (particle.pt() < minPt || particle.pt() > maxPt) { + if (particle.pt() < minPt || particle.pt() > maxPtSpectra) { continue; } @@ -1380,7 +1391,7 @@ struct UccZdc { std::vector vecFD; std::vector vecEff; - // Calculates the event weight, W_k + // Calculates the true Nch for (const auto& particle : mcParticles) { if (particle.eta() < minEta || particle.eta() > maxEta) { continue; @@ -1406,18 +1417,45 @@ struct UccZdc { if (!particle.isPhysicalPrimary()) { continue; } - - float pt{particle.pt()}; - pTs.emplace_back(pt); - vecEff.emplace_back(1.); - vecFD.emplace_back(1.); nchMult++; } if (nchMult < minNchSel) { continue; } - // printf("nchMult = %f | nchMC = %f | nchMult/nchMc = %f\n",nchMult,nchMC,nchMult/nchMC); + + // Calculates the event weight, W_k + for (const auto& particle : mcParticles) { + if (particle.eta() < minEta || particle.eta() > maxEta) { + continue; + } + if (particle.pt() < minPt || particle.pt() > maxPtSpectra) { + continue; + } + + auto charge{0.}; + // Get the MC particle + auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle != nullptr) { + charge = pdgParticle->Charge(); + } else { + continue; + } + + // Is it a charged particle? + if (std::abs(charge) < kMinCharge) { + continue; + } + // Is it a primary particle? + if (!particle.isPhysicalPrimary()) { + continue; + } + + const float pt{particle.pt()}; + pTs.emplace_back(pt); + vecEff.emplace_back(1.); + vecFD.emplace_back(1.); + } double p1, p2, p3, p4, w1, w2, w3, w4; p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0; @@ -1464,7 +1502,7 @@ struct UccZdc { std::vector vecFD; std::vector vecEff; - // Calculates the uncorrected Nch multiplicity + // Calculates the uncorrected Nch for (const auto& track : tracks) { // Track Selection if (!track.isGlobalTrack()) { @@ -1498,10 +1536,10 @@ struct UccZdc { continue; } - float pt{track.pt()}; + const float pt{track.pt()}; double fdValue{1.}; - int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; - double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; + const int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; + const double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; if (applyFD) fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); @@ -1532,10 +1570,10 @@ struct UccZdc { continue; } - float pt{track.pt()}; + const float pt{track.pt()}; double fdValue{1.}; - int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; - double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; + const int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; + const double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; if (applyFD) fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); @@ -1588,7 +1626,7 @@ struct UccZdc { } template - void eventSampling(const T& tracks, const U& normV0A, const U& normT0M, const U& sumZNs, const V& timeStamp) + void eventSampling(const T& tracks, const U& normV0A, const U& normT0M, const U& sumZNs, const U& sumZEMs, const V& timeStamp) { TRandom3 rndGen(timeStamp); std::vector vPoisson; @@ -1687,8 +1725,6 @@ struct UccZdc { pTs.emplace_back(pt); vecEff.emplace_back(effValue); vecFD.emplace_back(fdValue); - // To calculate event-averaged - registry.fill(HIST("NchVsZNVsPt"), nchMult, sumZNs, pt * (fdValue / effValue)); } } } else { @@ -1704,9 +1740,6 @@ struct UccZdc { pTs.emplace_back(track.pt()); vecEff.emplace_back(1.); vecFD.emplace_back(1.); - - // To calculate event-averaged - registry.fill(HIST("NchVsZNVsPt"), nchMult, sumZNs, track.pt()); } } @@ -1727,14 +1760,9 @@ struct UccZdc { const double numThreeParCorr{std::pow(p1, 3.) - (3. * p2 * p1) + (2. * p3)}; const double threeParCorr{numThreeParCorr / denThreeParCorr}; - hNchVsZN[replica]->Fill(nchMult, sumZNs); hNchVsV0A[replica]->Fill(nchMult, normV0A); hNchVsT0M[replica]->Fill(nchMult, normT0M); - pNchVsOneParCorrVsZN[replica]->Fill(nchMult, sumZNs, oneParCorr); - pNchVsTwoParCorrVsZN[replica]->Fill(nchMult, sumZNs, twoParCorr); - pNchVsThreeParCorrVsZN[replica]->Fill(nchMult, sumZNs, threeParCorr); - pNchVsOneParCorrVsT0M[replica]->Fill(nchMult, normT0M, oneParCorr); pNchVsTwoParCorrVsT0M[replica]->Fill(nchMult, normT0M, twoParCorr); pNchVsThreeParCorrVsT0M[replica]->Fill(nchMult, normT0M, threeParCorr); @@ -1742,14 +1770,19 @@ struct UccZdc { pNchVsOneParCorrVsV0A[replica]->Fill(nchMult, normV0A, oneParCorr); pNchVsTwoParCorrVsV0A[replica]->Fill(nchMult, normV0A, twoParCorr); pNchVsThreeParCorrVsV0A[replica]->Fill(nchMult, normV0A, threeParCorr); + + if (sumZEMs > zemCut) { + hNchVsZN[replica]->Fill(nchMult, sumZNs); + pNchVsOneParCorrVsZN[replica]->Fill(nchMult, sumZNs, oneParCorr); + pNchVsTwoParCorrVsZN[replica]->Fill(nchMult, sumZNs, twoParCorr); + pNchVsThreeParCorrVsZN[replica]->Fill(nchMult, sumZNs, threeParCorr); + } } // event per replica } // replica's loop } void loadCorrections(uint64_t timeStamp) { - // if (cfg.correctionsLoaded) return; - if (paTHEff.value.empty() == false) { cfg.hEfficiency = ccdb->getForTimeStamp(paTHEff, timeStamp); if (cfg.hEfficiency == nullptr) {