diff --git a/PWGJE/Tasks/jetSpectraEseTask.cxx b/PWGJE/Tasks/jetSpectraEseTask.cxx index 20e3f34589b..4a5240f27df 100644 --- a/PWGJE/Tasks/jetSpectraEseTask.cxx +++ b/PWGJE/Tasks/jetSpectraEseTask.cxx @@ -234,7 +234,8 @@ struct JetSpectraEseTask { static constexpr EventPlaneFiller PsiFillerEse = {true, false}; static constexpr EventPlaneFiller PsiFillerFalse = {false, false}; TRandom3* fRndm = new TRandom3(0); - static constexpr int NumSubSmpl = 10; + static constexpr int NumSubSmpl = 5; + static constexpr int NumSavedRhoFitEvents = 5; std::array, NumSubSmpl> hSameSub; void applySystematicPreset() @@ -359,6 +360,10 @@ struct JetSpectraEseTask { registry.get(HIST("trackQA/hRhoTrackCounter"))->GetXaxis()->SetBinLabel(3, "Tracks for rho(phi)"); registry.add("trackQA/hPhiPtsum", "jet sumpt;sum p_{T};entries", {HistType::kTH1F, {{40, 0., o2::constants::math::TwoPI}}}); + for (int i = 0; i < NumSavedRhoFitEvents; ++i) { + const auto eventName = fmt::format("rhoPhiFitEvents/event_{}", i); + registry.add(eventName, "event;#varphi;#Sigma #it{p}_{T} (GeV/#it{c})", HistType::kTH1F, {{1, 0., o2::constants::math::TwoPI}}); + } registry.add("eventQA/hfitPar0", "", {HistType::kTH2F, {{centAxis}, {fitAxis0}}}); registry.add("eventQA/hfitPar1", "", {HistType::kTH2F, {{centAxis}, {fitAxis13}}}); registry.add("eventQA/hfitPar2", "", {HistType::kTH2F, {{centAxis}, {fitAxis24}}}); @@ -659,7 +664,8 @@ struct JetSpectraEseTask { for (const auto& jet : jets) { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) continue; - if (!isAcceptedJet(jet)) { + // if (!isAcceptedJet(jet)) { + if (!isAcceptedJet>(jet)) { continue; } auto vCorr = corr(collision, jet); @@ -793,7 +799,7 @@ struct JetSpectraEseTask { for (const auto& jet : jets1) { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) continue; - if (!isAcceptedJet(jet)) { + if (!isAcceptedJet>(jet)) { continue; } auto vCorr = corr(c1, jet); @@ -1198,13 +1204,13 @@ struct JetSpectraEseTask { auto vec1{qVecNoESE(vec)}; epMap["FT0A"] = computeEP(vec1, 0, 2.0); ep3Map["FT0A"] = computeEP(vec1, 0, 3.0); + auto vec2{qVecNoESE(vec)}; + epMap["FT0C"] = computeEP(vec2, 0, 2.0); + ep3Map["FT0C"] = computeEP(vec2, 0, 3.0); + epMap["FV0A"] = computeEP(qVecNoESE(vec), 0, 2.0); + epMap["TPCpos"] = computeEP(qVecNoESE(vec), 1, 2.0); + epMap["TPCneg"] = computeEP(qVecNoESE(vec), 1, 2.0); if constexpr (P.psi) { - auto vec2{qVecNoESE(vec)}; - epMap["FT0C"] = computeEP(vec2, 0, 2.0); - ep3Map["FT0C"] = computeEP(vec2, 0, 3.0); - epMap["FV0A"] = computeEP(qVecNoESE(vec), 0, 2.0); - epMap["TPCpos"] = computeEP(qVecNoESE(vec), 1, 2.0); - epMap["TPCneg"] = computeEP(qVecNoESE(vec), 1, 2.0); if constexpr (P.hist) fillEP(/*std::make_index_sequence<5>{},*/ vec, epMap); auto cosPsi = [](float psiX, float psiY) { return (static_cast(psiX) == InvalidValue || static_cast(psiY) == InvalidValue) ? InvalidValue : std::cos(2.0 * (psiX - psiY)); }; @@ -1312,6 +1318,24 @@ struct JetSpectraEseTask { return true; } + std::shared_ptr getRhoPhiFitEvent(int eventIndex) + { + switch (eventIndex) { + case 0: + return registry.get(HIST("rhoPhiFitEvents/event_0")); + case 1: + return registry.get(HIST("rhoPhiFitEvents/event_1")); + case 2: + return registry.get(HIST("rhoPhiFitEvents/event_2")); + case 3: + return registry.get(HIST("rhoPhiFitEvents/event_3")); + case 4: + return registry.get(HIST("rhoPhiFitEvents/event_4")); + default: + return nullptr; + } + } + template std::unique_ptr fitRho(const Col& col, const EventPlane& ep, TTracks const& tracks, Jets const& jets) { @@ -1360,7 +1384,7 @@ struct JetSpectraEseTask { modulationFit->FixParameter(2, (ep.psi2 < 0) ? RecoDecay::constrainAngle(ep.psi2) : ep.psi2); modulationFit->FixParameter(4, (ep.psi3 < 0) ? RecoDecay::constrainAngle(ep.psi3) : ep.psi3); - hPhiPt->Fit(modulationFit.get(), "Q", "", 0, o2::constants::math::TwoPI); + hPhiPt->Fit(modulationFit.get(), "QN", "", 0, o2::constants::math::TwoPI); if constexpr (fillHist) { registry.fill(HIST("eventQA/hfitPar0"), col.centFT0M(), modulationFit->GetParameter(0)); @@ -1391,6 +1415,56 @@ struct JetSpectraEseTask { if constexpr (fillHist) { registry.fill(HIST("eventQA/hPValueCentCDF"), col.centFT0M(), cDF); registry.fill(HIST("eventQA/hCentChi2Ndf"), col.centFT0M(), chi2 / (static_cast(nDF))); + + static int numSavedRhoFitEvents{0}; + if (numSavedRhoFitEvents < NumSavedRhoFitEvents) { + auto savedEvent = getRhoPhiFitEvent(numSavedRhoFitEvents); + savedEvent->SetBins(hPhiPt->GetNbinsX(), 0., o2::constants::math::TwoPI); + for (int bin = 0; bin <= hPhiPt->GetNbinsX() + 1; ++bin) { + savedEvent->SetBinContent(bin, hPhiPt->GetBinContent(bin)); + savedEvent->SetBinError(bin, hPhiPt->GetBinError(bin)); + } + savedEvent->SetEntries(hPhiPt->GetEntries()); + + const double rho0 = modulationFit->GetParameter(0); + const double v2 = modulationFit->GetParameter(1); + const double psi2 = modulationFit->GetParameter(2); + const double v3 = modulationFit->GetParameter(3); + const double psi3 = modulationFit->GetParameter(4); + + auto rho0Function = std::make_unique("rho_0", "pol0", 0., o2::constants::math::TwoPI, TF1::EAddToList::kNo); + rho0Function->SetParameter(0, rho0); + rho0Function->SetParError(0, modulationFit->GetParError(0)); + rho0Function->SetParName(0, "rho_0"); + rho0Function->SetLineStyle(2); + + auto rhoPhiFunction = std::unique_ptr(static_cast(modulationFit->Clone("rho_phi"))); + rhoPhiFunction->SetParNames("rho_0", "v_2", "Psi_2", "v_3", "Psi_3"); + rhoPhiFunction->SetLineWidth(2); + rhoPhiFunction->ResetBit(TF1::kNotDraw); + + auto rhoV2Function = std::make_unique("rho_v2", "[0] * (1. + 2. * [1] * std::cos(2. * (x - [2])))", 0., o2::constants::math::TwoPI, TF1::EAddToList::kNo); + rhoV2Function->SetParameters(rho0, v2, psi2); + rhoV2Function->SetParError(0, modulationFit->GetParError(0)); + rhoV2Function->SetParError(1, modulationFit->GetParError(1)); + rhoV2Function->SetParError(2, modulationFit->GetParError(2)); + rhoV2Function->SetParNames("rho_0", "v_2", "Psi_2"); + rhoV2Function->SetLineStyle(7); + + auto rhoV3Function = std::make_unique("rho_v3", "[0] * (1. + 2. * [1] * std::cos(3. * (x - [2])))", 0., o2::constants::math::TwoPI, TF1::EAddToList::kNo); + rhoV3Function->SetParameters(rho0, v3, psi3); + rhoV3Function->SetParError(0, modulationFit->GetParError(0)); + rhoV3Function->SetParError(1, modulationFit->GetParError(3)); + rhoV3Function->SetParError(2, modulationFit->GetParError(4)); + rhoV3Function->SetParNames("rho_0", "v_3", "Psi_3"); + rhoV3Function->SetLineStyle(9); + + savedEvent->GetListOfFunctions()->Add(rho0Function.release()); + savedEvent->GetListOfFunctions()->Add(rhoPhiFunction.release()); + savedEvent->GetListOfFunctions()->Add(rhoV2Function.release()); + savedEvent->GetListOfFunctions()->Add(rhoV3Function.release()); + ++numSavedRhoFitEvents; + } } return modulationFit; @@ -1622,7 +1696,7 @@ struct JetSpectraEseTask { for (const auto& jet : jets) { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) continue; - if (!isAcceptedJet(jet)) { + if (!isAcceptedJet>(jet)) { continue; } auto jetPtBkg = jet.pt() - (collision.rho() * jet.area());