diff --git a/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx b/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx index c4faa289abb..5e6125d0bdb 100644 --- a/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx +++ b/PWGHF/HFC/Macros/DhCorrelationExtraction.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file DhCorrelationExtraction.cxx -/// \brief Class for D-h correlation extraction +/// \brief class for D-h correlation extraction /// \author Samuele Cattaruzzi /// \author Swapnesh Santosh Khade @@ -18,9 +18,11 @@ #include #include +#include #include #include #include +#include #include #include @@ -33,10 +35,19 @@ DhCorrelationExtraction::DhCorrelationExtraction() : // default constructor fFileMass(0x0), fFileSE(0x0), fFileME(0x0), + fFileFDTemplate(0x0), + fFileFDPromptFrac(0x0), + fFileSecPart(0x0), fDirMass(0x0), fDirSE(0x0), fDirME(0x0), + fDirSecPart(0x0), + fFilePromptMc(0x0), + fFileNonPromptMc(0x0), fCorrectedCorrHisto(0x0), + fCorrectedCorrHisto_BaselineSubtr(0x0), + fCorrectedCorrHisto_Reflected(0x0), + fCorrectedCorrHisto_Reflected_BaselineSubtr(0x0), fDmesonSpecies(kDsToKKPi), fDmesonLabel("Ds"), fNpools(9), @@ -46,21 +57,47 @@ DhCorrelationExtraction::DhCorrelationExtraction() : // default constructor fSubtractSoftPiME(kFALSE), fFileNameSE(""), fFileNameME(""), + fFileSecPartName(""), + fFilePromptMcRecName(""), + fFileNonPromptMcRecName(""), fDirNameSE(""), fDirNameME(""), + fDirSecPartName(""), fMassHistoNameSgn(""), fMassHistoNameBkg(""), fMassHistoNameSBs(""), fSECorrelSignalRegionName(""), fSECorrelSidebandsName(""), + fSECorrelSidebandLeftName(""), + fSECorrelSidebandRightName(""), fMECorrelSignalRegionName(""), fMECorrelSidebandsName(""), + fMECorrelSidebandLeftName(""), + fMECorrelSidebandRightName(""), + fFileFDTemplateName(""), + fFileFDPromptFracName(""), + fHistoFDTemplatePromptName(""), + fHistoFDTemplateNonPromptName(""), + fHistoFDPromptFracName(""), + fHistoPrimaryPartName(""), + fHistoAllPartName(""), fBkgScaleFactor(1.), fSgnYieldNorm(1.), - fRebin2Dhisto(kFALSE), + fBkgYield(1.), + fRebinAngCorr(kFALSE), + fRebinFDCorr(kFALSE), + fRebinSecPart(kFALSE), + fSidebandDivided(kFALSE), + fUseSidebLeft(kFALSE), + fUseSidebRight(kFALSE), fRebinAxisDeltaEta(1), fRebinAxisDeltaPhi(1), - fDebug(0) + fBinPtCand(0), + fBinPtHad(0), + fDebug(0), + fFDsubtraction(0), + fSecPartContamination(0), + fCorrBiasBtoD(0) { } @@ -68,10 +105,17 @@ DhCorrelationExtraction::DhCorrelationExtraction(const DhCorrelationExtraction& fFileMass(source.fFileMass), fFileSE(source.fFileSE), fFileME(source.fFileME), + fFileFDTemplate(source.fFileFDTemplate), + fFileFDPromptFrac(source.fFileFDPromptFrac), + fFileSecPart(source.fFileSecPart), + fFilePromptMc(source.fFilePromptMc), + fFileNonPromptMc(source.fFileNonPromptMc), fDirMass(source.fDirMass), - fDirSE(source.fDirSE), - fDirME(source.fDirME), + fDirSecPart(source.fDirSecPart), fCorrectedCorrHisto(source.fCorrectedCorrHisto), + fCorrectedCorrHisto_BaselineSubtr(source.fCorrectedCorrHisto_BaselineSubtr), + fCorrectedCorrHisto_Reflected(source.fCorrectedCorrHisto_Reflected), + fCorrectedCorrHisto_Reflected_BaselineSubtr(source.fCorrectedCorrHisto_Reflected_BaselineSubtr), fDmesonSpecies(source.fDmesonSpecies), fDmesonLabel(source.fDmesonLabel), fNpools(source.fNpools), @@ -81,21 +125,47 @@ DhCorrelationExtraction::DhCorrelationExtraction(const DhCorrelationExtraction& fSubtractSoftPiME(source.fSubtractSoftPiME), fFileNameSE(source.fFileNameSE), fFileNameME(source.fFileNameME), + fFileSecPartName(source.fFileSecPartName), + fFilePromptMcRecName(source.fFilePromptMcRecName), + fFileNonPromptMcRecName(source.fFileNonPromptMcRecName), fDirNameSE(source.fDirNameSE), fDirNameME(source.fDirNameME), + fDirSecPartName(source.fDirSecPartName), fMassHistoNameSgn(source.fMassHistoNameSgn), fMassHistoNameBkg(source.fMassHistoNameBkg), fMassHistoNameSBs(source.fMassHistoNameSBs), fSECorrelSignalRegionName(source.fSECorrelSignalRegionName), fSECorrelSidebandsName(source.fSECorrelSidebandsName), + fSECorrelSidebandLeftName(source.fSECorrelSidebandLeftName), + fSECorrelSidebandRightName(source.fSECorrelSidebandRightName), fMECorrelSignalRegionName(source.fMECorrelSignalRegionName), fMECorrelSidebandsName(source.fMECorrelSidebandsName), + fMECorrelSidebandLeftName(source.fMECorrelSidebandLeftName), + fMECorrelSidebandRightName(source.fMECorrelSidebandRightName), + fFileFDTemplateName(source.fFileFDTemplateName), + fFileFDPromptFracName(source.fFileFDPromptFracName), + fHistoFDTemplatePromptName(source.fHistoFDTemplatePromptName), + fHistoFDTemplateNonPromptName(source.fHistoFDTemplateNonPromptName), + fHistoFDPromptFracName(source.fHistoFDPromptFracName), + fHistoPrimaryPartName(source.fHistoPrimaryPartName), + fHistoAllPartName(source.fHistoAllPartName), fBkgScaleFactor(source.fBkgScaleFactor), fSgnYieldNorm(source.fSgnYieldNorm), - fRebin2Dhisto(source.fRebin2Dhisto), + fBkgYield(source.fBkgYield), + fRebinAngCorr(source.fRebinAngCorr), + fRebinFDCorr(source.fRebinFDCorr), + fRebinSecPart(source.fRebinSecPart), + fSidebandDivided(source.fSidebandDivided), + fUseSidebLeft(source.fUseSidebLeft), + fUseSidebRight(source.fUseSidebRight), fRebinAxisDeltaEta(source.fRebinAxisDeltaEta), fRebinAxisDeltaPhi(source.fRebinAxisDeltaPhi), - fDebug(source.fDebug) + fBinPtCand(source.fBinPtCand), + fBinPtHad(source.fBinPtHad), + fDebug(source.fDebug), + fFDsubtraction(source.fFDsubtraction), + fSecPartContamination(source.fSecPartContamination), + fCorrBiasBtoD(source.fCorrBiasBtoD) { } @@ -119,7 +189,6 @@ Bool_t DhCorrelationExtraction::SetDmesonSpecie(DmesonSpecie k) } else { fDmesonLabel = "Dstar"; } - fDmesonSpecies = k; return kTRUE; } @@ -149,20 +218,43 @@ Bool_t DhCorrelationExtraction::ExtractCorrelations(Double_t PtCandMin, Double_t TH2D* h2D_Sideb; TH2D* h2D_Subtr; + TH2D* h2D_FDTemplatePrompt; + TH2D* h2D_FDTemplateNonPrompt; + TH1D* h1D_Sign; TH1D* h1D_Sideb; TH1D* h1D_Subtr; + TH1D* h1D_SignNorm; + TH1D* h1D_SidebNorm; TH1D* h1D_SubtrNorm; + TH1D* h1D_FDTemplatePrompt; + TH1D* h1D_FDTemplateNonPrompt; + TH1D* h1D_TemplateTotal; + TH1D* h1D_SubtrFDNorm; + TH1D* h1D_PrimaryPartCorr; + TH1D* h1D_AllPartCorr; + TH1D* h1D_SecPartFrac; + TH1D* h1D_SubtrNorm_SecPart; + TH1D* h1D_BaselineSubtr; + TH1D* h1D_ReflCorr; + TH1D* h1D_ReflCorr_BaselineSubtr; + TH1D* hModul; + TH1D* hBeforeModulCorr; + + Double_t FDPromptFrac; // if (fIntegratePtBins && iBinPtHad>0) continue; for (int iPool = 0; iPool < fNpools; iPool++) { - // Retrieve 2D plots for SE and ME, signal and bkg regions, for each pTbin and pool hSE_Sign[iPool] = GetCorrelHisto(kSE, kSign, iPool, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + std::cout << "Got SE histogram signal region" << std::endl; hME_Sign[iPool] = GetCorrelHisto(kME, kSign, iPool, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + std::cout << "Got ME histogram signal region" << std::endl; hSE_Sideb[iPool] = GetCorrelHisto(kSE, kSideb, iPool, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + std::cout << "Got SE histogram sdeband region" << std::endl; hME_Sideb[iPool] = GetCorrelHisto(kME, kSideb, iPool, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + std::cout << "Got ME histogram sdeband region" << std::endl; hSE_Sign[iPool]->Sumw2(); hME_Sign[iPool]->Sumw2(); @@ -170,7 +262,7 @@ Bool_t DhCorrelationExtraction::ExtractCorrelations(Double_t PtCandMin, Double_t hME_Sideb[iPool]->Sumw2(); // rebin axes deltaEta and deltaPhi - if (fRebin2Dhisto) { + if (fRebinAngCorr) { hSE_Sign[iPool]->Rebin2D(fRebinAxisDeltaEta, fRebinAxisDeltaPhi); // Xaxis: deltaEta, Yaxis: deltaPhi hSE_Sideb[iPool]->Rebin2D(fRebinAxisDeltaEta, fRebinAxisDeltaPhi); hME_Sign[iPool]->Rebin2D(fRebinAxisDeltaEta, fRebinAxisDeltaPhi); @@ -178,30 +270,39 @@ Bool_t DhCorrelationExtraction::ExtractCorrelations(Double_t PtCandMin, Double_t if (fSubtractSoftPiME) { hME_Sideb_SoftPi[iPool]->Rebin2D(fRebinAxisDeltaEta, fRebinAxisDeltaPhi); } + std::cout << "SE and ME histograms rebinned" << std::endl; } if (fDebug >= 1) { - TCanvas* c = new TCanvas(Form("cSEME_Original_%d_%1.1fto%1.1f", iPool, PtHadMin, PtHadMax), Form("cSEME_Original_%s_pool%d_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), iPool, PtCandMin, PtCandMax, PtHadMin, PtHadMax), 100, 100, 1600, 900); - c->Divide(2, 2); + TCanvas* c = new TCanvas(Form("cSE_Original_%d_%1.1fto%1.1f", iPool, PtHadMin, PtHadMax), Form("cSE_Original_%s_pool%d_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), iPool, PtCandMin, PtCandMax, PtHadMin, PtHadMax), 100, 100, 1600, 900); + c->Divide(2, 1); c->cd(1); hSE_Sign[iPool]->SetMinimum(0); hSE_Sign[iPool]->Draw("lego2"); c->cd(2); - hME_Sign[iPool]->SetMinimum(0); - hME_Sign[iPool]->Draw("lego2"); - c->cd(3); hSE_Sideb[iPool]->SetMinimum(0); hSE_Sideb[iPool]->Draw("lego2"); - c->cd(4); + c->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrSE_Original_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrSE_Original_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } + + if (fDebug >= 1) { + TCanvas* c = new TCanvas(Form("cME_Original_%d_%1.1fto%1.1f", iPool, PtHadMin, PtHadMax), Form("cME_Original_%s_pool%d_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), iPool, PtCandMin, PtCandMax, PtHadMin, PtHadMax), 100, 100, 1600, 900); + c->Divide(2, 1); + c->cd(1); + hME_Sign[iPool]->SetMinimum(0); + hME_Sign[iPool]->Draw("lego2"); + c->cd(2); hME_Sideb[iPool]->SetMinimum(0); hME_Sideb[iPool]->Draw("lego2"); - c->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrSEandME_Original_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); - c->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrSEandME_Original_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrME_Original_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrME_Original_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); } // Scale bkg plots by ratio of signal region/sidebands hSE_Sideb[iPool]->Scale(fBkgScaleFactor); - hME_Sideb[iPool]->Scale(fBkgScaleFactor); + hME_Sideb[iPool]->Scale(fBkgScaleFactor); // when normalised this factor should cancel out + std::cout << "[INFO] fBkgScaleFactor = " << fBkgScaleFactor << std::endl; hSE_Sideb[iPool]->SetEntries(hSE_Sideb[iPool]->GetEntries() * fBkgScaleFactor); hME_Sideb[iPool]->SetEntries(hME_Sideb[iPool]->GetEntries() * fBkgScaleFactor); @@ -258,8 +359,8 @@ Bool_t DhCorrelationExtraction::ExtractCorrelations(Double_t PtCandMin, Double_t c->cd(6); hCorr_Sideb[iPool]->SetMinimum(0); hCorr_Sideb[iPool]->Draw("lego2"); - c->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrSEandME_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); - c->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrSEandME_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrSEandME_%s_Canvas_PtCand%.0fto%.0f_Pool%d_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, iPool, PtHadMin, PtHadMax)); + c->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrSEandME_%s_Canvas_PtCand%.0fto%.0f_Pool%d_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, iPool, PtHadMin, PtHadMax)); } // Pools integration @@ -288,6 +389,31 @@ Bool_t DhCorrelationExtraction::ExtractCorrelations(Double_t PtCandMin, Double_t c2D->SaveAs(Form("Output_CorrelationExtraction_%s_png/h2D_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); c2D->SaveAs(Form("Output_CorrelationExtraction_%s_Root/h2D_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + // Get FD correlations for FD subtraction + if (fFDsubtraction) { + h2D_FDTemplatePrompt = GetFDTemplateHisto(kPrompt, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + h2D_FDTemplateNonPrompt = GetFDTemplateHisto(kFD, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + // h1D_BaselineSubtr + FDPromptFrac = GetFDPromptFrac(PtCandMin, PtCandMax, PtHadMin, PtHadMax); + + if (fRebinFDCorr) { + h2D_FDTemplatePrompt->Rebin2D(fRebinAxisDeltaEta, fRebinAxisDeltaPhi); + h2D_FDTemplateNonPrompt->Rebin2D(fRebinAxisDeltaEta, fRebinAxisDeltaPhi); + } + + if (fDebug >= 1) { + TCanvas* c = new TCanvas(Form("cFDTemplate_PtCand%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax), Form("cFDTemplate_%s_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->Divide(2, 1); + c->cd(1); + h2D_FDTemplatePrompt->SetMinimum(0); + h2D_FDTemplatePrompt->Draw("lego2"); + c->cd(2); + h2D_FDTemplateNonPrompt->SetMinimum(0); + h2D_FDTemplateNonPrompt->Draw("lego2"); + c->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrFDTemplate_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } + } + // Bkg subtraction (2D plot) TCanvas* c2D_Sub = new TCanvas(Form("c2D_Subtr_IntPools_PtHAd%.0fto%.0f", PtHadMin, PtHadMax), Form("c2D_%s_Subtr_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1500, 800); h2D_Subtr = reinterpret_cast(h2D_Sign->Clone("h2D_Subtr")); @@ -315,7 +441,7 @@ Bool_t DhCorrelationExtraction::ExtractCorrelations(Double_t PtCandMin, Double_t h1D_Subtr->SetTitle("Signal region after sideb. subt. corr."); // Draw 1D plots (Signal region, Sidebands, S-SB (subtr.)) - TCanvas* c1D = new TCanvas(Form("c1D_IntPools_%.0fto%.0f", PtHadMin, PtHadMax), Form("c1D_%s_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1600, 500); + TCanvas* c1D = new TCanvas(Form("c1D_IntPools_PtCand%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax), Form("c1D_%s_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1600, 500); c1D->Divide(3, 1); c1D->cd(1); h1D_Sign->Draw(); @@ -326,20 +452,378 @@ Bool_t DhCorrelationExtraction::ExtractCorrelations(Double_t PtCandMin, Double_t c1D->SaveAs(Form("Output_CorrelationExtraction_%s_png/h1D_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); c1D->SaveAs(Form("Output_CorrelationExtraction_%s_Root/h1D_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + if (fDebug >= 1) { + h1D_SignNorm = reinterpret_cast(h1D_Sign->Clone("h1D_Sign_Norm")); + h1D_SidebNorm = reinterpret_cast(h1D_Sideb->Clone("h1D_Sideb_Norm")); + h1D_SignNorm->Scale(1. / (fSgnYieldNorm + fBkgYield)); + // h1D_SidebNorm -> Scale(1./fBkgYield); + h1D_SidebNorm->Scale(1. / fBkgScaleFactor); + h1D_SidebNorm->Scale(1. / fSBYield); + h1D_SignNorm->SetMarkerStyle(kFullCircle); + h1D_SignNorm->SetMarkerSize(1.2); + h1D_SignNorm->SetLineColor(kRed); + h1D_SignNorm->SetMarkerColor(kRed); + h1D_SignNorm->SetLineWidth(2); + h1D_SidebNorm->SetMinimum(0); + h1D_SidebNorm->SetMarkerStyle(kFullSquare); + h1D_SidebNorm->SetMarkerSize(1.2); + h1D_SidebNorm->SetLineColor(kBlue); + h1D_SidebNorm->SetMarkerColor(kBlue); + h1D_SidebNorm->SetLineWidth(2); + h1D_SidebNorm->SetTitle(Form("%.0f < p_{T} < %.0f", PtCandMin, PtCandMax)); + TCanvas* c = new TCanvas(Form("c_IntPools_PtCand%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax), ""); + c->cd(); + h1D_SidebNorm->Draw(); + h1D_SignNorm->Draw("same"); + c->SaveAs(Form("Output_CorrelationExtraction_%s_png/ComparisonSignalSidebCorr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->SaveAs(Form("Output_CorrelationExtraction_%s_Root/ComparisonSignalSidebCorr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } // Apply normalization to number of triggers h1D_SubtrNorm = reinterpret_cast(h1D_Subtr->Clone("h1D_SubtrNorm")); h1D_SubtrNorm->Sumw2(); h1D_SubtrNorm->Scale(1. / fSgnYieldNorm); h1D_SubtrNorm->SetTitle("Signal region after sideb. subt. corr. - Normalized to # of triggers"); - fCorrectedCorrHisto = reinterpret_cast(h1D_SubtrNorm->Clone(Form("hCorrectedCorr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + // Correction for bias B to D topologies + if (fCorrBiasBtoD) { + hModul = EvaluateMCClosModulations(PtCandMin, PtCandMax, PtHadMin, PtHadMax); + TCanvas* c1D_corrBbias = new TCanvas(Form("c1D_corrBbias_IntPools_PtCand%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax), Form("c1D_corrBbias_%s_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1600, 500); + c1D_corrBbias->cd(); + hBeforeModulCorr = reinterpret_cast(h1D_SubtrNorm->Clone("hBeforeModulCorr")); + hBeforeModulCorr->SetLineColor(kViolet - 3); + hBeforeModulCorr->GetYaxis()->SetRangeUser(0., 5.); + hBeforeModulCorr->Draw(); + h1D_SubtrNorm->Multiply(hModul); + h1D_SubtrNorm->Draw("same"); + c1D_corrBbias->SaveAs(Form("Output_CorrelationExtraction_%s_png/ComparisonCorrBiasBtoD_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c1D_corrBbias->SaveAs(Form("Output_CorrelationExtraction_%s_Root/ComparisonCorrBiasBtoD_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + + TFile* file = new TFile(Form("Output_CorrelationExtraction_%s_Root/SystematicCorrBiasBtoD_%s_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax), "RECREATE"); // Open file in write mode + TH1D* h1D_SubtrNorm_Clone = reinterpret_cast(h1D_SubtrNorm->Clone("h1D_SubtrNorm_Clone")); + h1D_SubtrNorm_Clone = ReflectCorrHistogram(h1D_SubtrNorm_Clone); + hBeforeModulCorr = ReflectCorrHistogram(hBeforeModulCorr); + TH1D* hSystematicCorrBiasBtoD = reinterpret_cast(h1D_SubtrNorm_Clone->Clone("hSystematicCorrBiasBtoD")); + hSystematicCorrBiasBtoD->Add(h1D_SubtrNorm_Clone, hBeforeModulCorr, 1, -1); + // Set bin contents to absolute values + for (int i = 1; i <= hSystematicCorrBiasBtoD->GetNbinsX(); ++i) { + hSystematicCorrBiasBtoD->SetBinContent(i, std::abs(hSystematicCorrBiasBtoD->GetBinContent(i)) / TMath::Sqrt(12)); + hSystematicCorrBiasBtoD->SetBinError(i, 0.); + } + hSystematicCorrBiasBtoD->Write(); + file->Close(); + } + + // Secondary particle contamination + if (fSecPartContamination) { + h1D_PrimaryPartCorr = GetCorrelHistoSecondaryPart(kPrimaryPart, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + h1D_AllPartCorr = GetCorrelHistoSecondaryPart(kAllPart, PtCandMin, PtCandMax, PtHadMin, PtHadMax); + h1D_PrimaryPartCorr->Sumw2(); + h1D_AllPartCorr->Sumw2(); + if (fRebinSecPart) { + h1D_PrimaryPartCorr->RebinX(fRebinAxisDeltaPhi); // Xaxis: deltaPhi + h1D_AllPartCorr->RebinX(fRebinAxisDeltaPhi); + std::cout << "Secondary particle histogram rebinned" << std::endl; + } + h1D_SecPartFrac = reinterpret_cast(h1D_PrimaryPartCorr->Clone(Form("hCorrRatio_PtD%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + h1D_SecPartFrac->Sumw2(); + h1D_SecPartFrac->Divide(h1D_PrimaryPartCorr, h1D_AllPartCorr, 1., 1., "B"); + + TCanvas* c1D = new TCanvas(Form("c1D_CorrPrimaryPart_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax), Form("c1D_%s_CorrPrimaryPart_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax)); + c1D->cd(); + SetTH1HistoStyle(h1D_SecPartFrac, Form("%.0f < p_{T} < %.0f GeV/c", PtCandMin, PtCandMax), "#Delta#phi [rad]", "#frac{primary part.}{part. selected}"); + h1D_SecPartFrac->Draw(); + c1D->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrPrimaryPartRatio_%s_Canvas_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c1D->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrPrimaryPartRatio_%s_Canvas_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + + h1D_SubtrNorm_SecPart = reinterpret_cast(h1D_SubtrNorm->Clone("h1D_SubtrNorm_SecPart")); + h1D_SubtrNorm_SecPart->Sumw2(); + Int_t nBinsPhi = h1D_SubtrNorm_SecPart->GetNbinsX(); + if (nBinsPhi != h1D_SecPartFrac->GetNbinsX()) { + std::cout << "[ERROR]: nBinsPhi different between h1D_SubtrNorm and h1D_SecPartFrac" << std::endl; + return kFALSE; + } + h1D_SubtrNorm_SecPart->Multiply(h1D_SecPartFrac); + } + + // FD Subtraction + if (fFDsubtraction) { + h1D_FDTemplatePrompt = reinterpret_cast(h2D_FDTemplatePrompt->ProjectionY("h1D_FDTemplatePrompt")); + h1D_FDTemplateNonPrompt = reinterpret_cast(h2D_FDTemplateNonPrompt->ProjectionY("h1D_FDTemplateNonPrompt")); + + h1D_FDTemplatePrompt->Scale(1. / h1D_FDTemplatePrompt->GetXaxis()->GetBinWidth(1)); + h1D_FDTemplateNonPrompt->Scale(1. / h1D_FDTemplateNonPrompt->GetXaxis()->GetBinWidth(1)); + + h1D_TemplateTotal = reinterpret_cast(h1D_FDTemplatePrompt->Clone("h1D_TemplateTotal")); + h1D_TemplateTotal->Sumw2(); + h1D_TemplateTotal->Scale(FDPromptFrac); + h1D_TemplateTotal->Add(h1D_FDTemplateNonPrompt, 1 - FDPromptFrac); + + if (fDebug >= 1) { + TCanvas* c = new TCanvas(Form("cFDTemplate_1D_%1.1fto%1.1f", PtHadMin, PtHadMax), Form("cFDTemplate_%s_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->cd(); + h1D_TemplateTotal->SetMinimum(0); + h1D_FDTemplateNonPrompt->SetMinimum(0); + h1D_TemplateTotal->SetMarkerColor(kGreen); + h1D_TemplateTotal->SetLineColor(kGreen); + h1D_TemplateTotal->SetLineWidth(2); + h1D_TemplateTotal->SetMarkerStyle(kFullCircle); + h1D_FDTemplatePrompt->SetMarkerColor(kRed); + h1D_FDTemplatePrompt->SetLineColor(kRed); + h1D_FDTemplatePrompt->SetLineWidth(2); + h1D_FDTemplatePrompt->SetMarkerStyle(kFullCircle); + h1D_FDTemplateNonPrompt->SetMarkerColor(kBlue); + h1D_FDTemplateNonPrompt->SetLineColor(kBlue); + h1D_FDTemplateNonPrompt->SetLineWidth(2); + h1D_FDTemplateNonPrompt->SetMarkerStyle(kFullCircle); + h1D_FDTemplateNonPrompt->Draw(); + h1D_FDTemplatePrompt->Draw("same"); + h1D_TemplateTotal->Draw("same"); + TLegend* lFD = new TLegend(); + lFD->AddEntry(h1D_TemplateTotal, "Total template"); + lFD->AddEntry(h1D_FDTemplatePrompt, "Prompt Template"); + lFD->AddEntry(h1D_FDTemplateNonPrompt, "Non prompt template"); + lFD->Draw("same"); + c->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrFDTemplate_1D_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrFDTemplate_1D_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } + + Double_t BaselineFD = CalculateBaseline(h1D_TemplateTotal, kTRUE); + Double_t BaselineData; + if (fSecPartContamination) { + BaselineData = CalculateBaseline(h1D_SubtrNorm_SecPart, kTRUE); + } else { + BaselineData = CalculateBaseline(h1D_SubtrNorm, kTRUE); + } + + std::cout << "===================== " << std::endl; + std::cout << "Baseline FD: " << BaselineFD << std::endl; + std::cout << "Baseline Data: " << BaselineData << std::endl; + std::cout << "===================== " << std::endl; + std::cout << " " << std::endl; + + Double_t Baselinediff = BaselineData - BaselineFD; + TH1D* hBaselineDiff = reinterpret_cast(h1D_FDTemplateNonPrompt->Clone("hBaselineDiff")); + for (int iBin = 0; iBin < hBaselineDiff->GetNbinsX(); iBin++) { + hBaselineDiff->SetBinContent(iBin + 1, Baselinediff); + } + h1D_FDTemplateNonPrompt->Add(hBaselineDiff); + h1D_TemplateTotal->Add(hBaselineDiff); + if (fSecPartContamination) { + h1D_SubtrFDNorm = reinterpret_cast(h1D_SubtrNorm_SecPart->Clone("h1D_SubtrFDNorm")); + } else { + h1D_SubtrFDNorm = reinterpret_cast(h1D_SubtrNorm->Clone("h1D_SubtrFDNorm")); + } + h1D_FDTemplateNonPrompt->Scale(1 - FDPromptFrac); + h1D_SubtrFDNorm->Add(h1D_FDTemplateNonPrompt, -1); + h1D_SubtrFDNorm->Scale(1. / FDPromptFrac); + + if (fDebug >= 1) { + TCanvas* c1 = new TCanvas(Form("cFDTemplateSubtr_%1.1fto%1.1f", PtHadMin, PtHadMax), Form("cFDTemplateSubtr_%s_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax), 100, 100, 1600, 900); + c1->cd(); + h1D_SubtrNorm->SetLineColor(kRed); + h1D_SubtrNorm_SecPart->SetLineColor(kOrange); + h1D_FDTemplateNonPrompt->SetLineColor(kBlue); + h1D_SubtrFDNorm->SetLineColor(kGreen); + h1D_TemplateTotal->SetLineColor(kMagenta); + h1D_SubtrNorm->SetMinimum(0); + h1D_SubtrNorm_SecPart->SetMinimum(0); + h1D_FDTemplateNonPrompt->SetMinimum(0); + h1D_SubtrFDNorm->SetMinimum(0); + // h1D_SubtrNorm -> GetYaxis() -> SetRangeUser(0., 8.); + h1D_SubtrNorm->SetMarkerStyle(kFullCircle); + h1D_SubtrNorm->SetMarkerSize(1.2); + h1D_SubtrNorm->SetMarkerColor(kRed); + h1D_SubtrNorm->SetLineWidth(2); + h1D_SubtrNorm_SecPart->SetMarkerStyle(kFullCircle); + h1D_SubtrNorm_SecPart->SetMarkerSize(1.2); + h1D_SubtrNorm_SecPart->SetMarkerColor(kOrange); + h1D_SubtrNorm_SecPart->SetLineWidth(2); + h1D_SubtrFDNorm->SetMarkerStyle(kFullCircle); + h1D_SubtrFDNorm->SetMarkerSize(1.2); + h1D_SubtrFDNorm->SetMarkerColor(kGreen); + h1D_SubtrFDNorm->SetLineWidth(2); + h1D_SubtrNorm->GetYaxis()->SetTitle("#frac{1}{N_{D}} #frac{dN^{assoc. part}}{d#Delta#phi}"); + h1D_SubtrNorm_SecPart->GetYaxis()->SetTitle("#frac{1}{N_{D}} #frac{dN^{assoc. part}}{d#Delta#phi}"); + if (fSecPartContamination) { + h1D_SubtrNorm_SecPart->Draw(); + } else { + h1D_SubtrNorm->Draw(); + } + // h1D_FDTemplateNonPrompt -> Draw("same"); + h1D_SubtrFDNorm->Draw("same"); + h1D_TemplateTotal->Draw("same"); + c1->SaveAs(Form("Output_CorrelationExtraction_%s_png/CorrFDTemplateSubtr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + c1->SaveAs(Form("Output_CorrelationExtraction_%s_Root/CorrFDTemplateSubtr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } + } + + if (fFDsubtraction) { + fCorrectedCorrHisto = reinterpret_cast(h1D_SubtrFDNorm->Clone(Form("hCorrectedCorr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + } else if (fSecPartContamination) { + fCorrectedCorrHisto = reinterpret_cast(h1D_SubtrNorm_SecPart->Clone(Form("hCorrectedCorr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + } else { + fCorrectedCorrHisto = reinterpret_cast(h1D_SubtrNorm->Clone(Form("hCorrectedCorr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + } + + std::cout << "Analysis steps completed - baseline subtraction missing" << std::endl; // Draw 1D plots (Signal region, normalized) TCanvas* cFinal = new TCanvas(Form("cFinal_%.0fto%.0f", PtHadMin, PtHadMax), Form("cFinal_%s_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1200, 700); + h1D_SubtrNorm->SetLineColor(kBlue + 1); + h1D_SubtrNorm->SetMarkerColor(kBlue + 1); + h1D_SubtrNorm->SetMarkerStyle(kFullCircle); + h1D_SubtrNorm->SetMinimum(0); h1D_SubtrNorm->Draw(); + if (fSecPartContamination) { + h1D_SubtrNorm_SecPart->SetLineColor(kRed + 1); + h1D_SubtrNorm_SecPart->SetMarkerColor(kRed + 1); + h1D_SubtrNorm_SecPart->SetMarkerStyle(kFullCircle); + h1D_SubtrNorm_SecPart->Draw("same"); + } + if (fFDsubtraction) { + h1D_SubtrFDNorm->SetLineColor(kGreen + 2); + h1D_SubtrFDNorm->SetMarkerColor(kGreen + 2); + h1D_SubtrFDNorm->SetMarkerStyle(kFullCircle); + h1D_SubtrFDNorm->Draw("same"); + } + if (fFDsubtraction) + h1D_TemplateTotal->Draw("same"); + TLegend* lFinal = new TLegend(); + lFinal->AddEntry(h1D_SubtrNorm, "Corr. after bkg subtr."); + if (fFDsubtraction) + lFinal->AddEntry(h1D_TemplateTotal, "CR Mode 2 total template"); + if (fSecPartContamination) { + lFinal->AddEntry(h1D_SubtrNorm_SecPart, "Corr. after sec. part. correction"); + } + if (fFDsubtraction) { + lFinal->AddEntry(h1D_SubtrFDNorm, "Corr. FD subtr."); + } + lFinal->Draw("same"); cFinal->SaveAs(Form("Output_CorrelationExtraction_%s_png/AzimCorrDistr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); cFinal->SaveAs(Form("Output_CorrelationExtraction_%s_Root/AzimCorrDistr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + // Baseline subtraction + Double_t BaselineData, BaselineDataErr; + TH1D* hBaseline = reinterpret_cast(h1D_SubtrNorm->Clone("hBaseline")); + hBaseline->Sumw2(); + if (fFDsubtraction) { + BaselineData = CalculateBaseline(h1D_SubtrFDNorm, kTRUE, kFALSE); // introduced kFALSE + BaselineDataErr = CalculateBaselineError(h1D_SubtrFDNorm, kTRUE, kFALSE); + for (int iBin = 0; iBin < hBaseline->GetNbinsX(); iBin++) { + hBaseline->SetBinContent(iBin + 1, BaselineData); + hBaseline->SetBinError(iBin + 1, BaselineDataErr); + } + h1D_BaselineSubtr = reinterpret_cast(h1D_SubtrFDNorm->Clone("h1D_BaselineSubtr")); + h1D_BaselineSubtr->Add(hBaseline, -1.); + } else if (fSecPartContamination) { + BaselineData = CalculateBaseline(h1D_SubtrNorm_SecPart, kTRUE, kFALSE); + BaselineDataErr = CalculateBaselineError(h1D_SubtrNorm_SecPart, kTRUE, kFALSE); + for (int iBin = 0; iBin < hBaseline->GetNbinsX(); iBin++) { + hBaseline->SetBinContent(iBin + 1, BaselineData); + hBaseline->SetBinError(iBin + 1, BaselineDataErr); + } + h1D_BaselineSubtr = reinterpret_cast(h1D_SubtrNorm_SecPart->Clone("h1D_BaselineSubtr")); + h1D_BaselineSubtr->Add(hBaseline, -1.); + } else { + BaselineData = CalculateBaseline(h1D_SubtrNorm, kTRUE, kFALSE); + BaselineDataErr = CalculateBaselineError(h1D_SubtrNorm, kTRUE, kFALSE); + for (int iBin = 0; iBin < hBaseline->GetNbinsX(); iBin++) { + hBaseline->SetBinContent(iBin + 1, BaselineData); + hBaseline->SetBinError(iBin + 1, BaselineDataErr); + } + h1D_BaselineSubtr = reinterpret_cast(h1D_SubtrNorm->Clone("h1D_BaselineSubtr")); + h1D_BaselineSubtr->Add(hBaseline, -1.); + } + + fCorrectedCorrHisto_BaselineSubtr = reinterpret_cast(h1D_BaselineSubtr->Clone(Form("hCorrectedCorrBaselineSubtr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + + TCanvas* cFinal_BaselineSubtr = new TCanvas(Form("cFinal_BaselineSubtr_%.0fto%.0f", PtHadMin, PtHadMax), Form("cFinal_BaselineSubtr_%s_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1200, 700); + h1D_BaselineSubtr->SetMarkerColor(kOrange + 8); + h1D_BaselineSubtr->SetLineColor(kOrange + 8); + h1D_BaselineSubtr->GetYaxis()->SetRangeUser(-0.2, 8.); + h1D_BaselineSubtr->Draw(); + if (fFDsubtraction) { + h1D_SubtrFDNorm->Draw("same"); + } else if (fSecPartContamination) { + h1D_SubtrNorm_SecPart->Draw("same"); + } else { + h1D_SubtrNorm->Draw("same"); + } + hBaseline->SetMarkerColor(kPink - 6); + hBaseline->SetMarkerStyle(kFullSquare); + hBaseline->SetLineColor(kPink - 6); + hBaseline->Draw("same"); + + cFinal_BaselineSubtr->SaveAs(Form("Output_CorrelationExtraction_%s_png/AzimCorrDistr_BaselineSubtr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + cFinal_BaselineSubtr->SaveAs(Form("Output_CorrelationExtraction_%s_Root/AzimCorrDistr_BaselineSubtr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + + // Reflected histograms + if (fFDsubtraction) { + h1D_ReflCorr = ReflectCorrHistogram(h1D_SubtrFDNorm); + } else if (fSecPartContamination) { + h1D_ReflCorr = ReflectCorrHistogram(h1D_SubtrNorm_SecPart); + } else { + h1D_ReflCorr = ReflectCorrHistogram(h1D_SubtrNorm); + } + + /* used as control using Run2 reflection function + if (fFDsubtraction) { + h1D_ReflCorr = ReflectHistoRun2(h1D_SubtrFDNorm, 0.5); + } else if (fSecPartContamination) { + h1D_ReflCorr = ReflectHistoRun2(h1D_SubtrNorm_SecPart, 0.5); + } else { + h1D_ReflCorr = ReflectHistoRun2(h1D_SubtrNorm, 0.5); + }*/ + + TCanvas* cFinal_Reflected = new TCanvas(Form("cFinal_Reflected_%.0fto%.0f", PtHadMin, PtHadMax), Form("cFinal_Reflected_%s_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1200, 700); + cFinal_Reflected->cd(); + SetTH1HistoStyle(h1D_ReflCorr, Form("%.0f < p_{T} < %.0f GeV/c", PtCandMin, PtCandMax), "#Delta#phi [rad]", "#frac{1}{N_{D}}#frac{dN^{assoc}}{d#Delta#phi} [rad^{-1}]", kFullCircle, kOrange + 8, 1.6, kOrange + 8, 3); + h1D_ReflCorr->SetMinimum(0); + h1D_ReflCorr->Draw(); + cFinal_Reflected->SaveAs(Form("Output_CorrelationExtraction_%s_png/AzimCorrDistr_Reflected_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + cFinal_Reflected->SaveAs(Form("Output_CorrelationExtraction_%s_Root/AzimCorrDistr_Reflected_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + + // Reflected histograms baseline subtracted + TH1D* hBaseline_Refl = reinterpret_cast(h1D_ReflCorr->Clone("hBaseline_Refl")); + hBaseline_Refl->Sumw2(); + BaselineData = CalculateBaseline(h1D_ReflCorr, kFALSE, kTRUE); + BaselineDataErr = CalculateBaselineError(h1D_ReflCorr, kFALSE, kTRUE); + + for (int iBin = 0; iBin < hBaseline_Refl->GetNbinsX(); iBin++) { + hBaseline_Refl->SetBinContent(iBin + 1, BaselineData); + hBaseline_Refl->SetBinError(iBin + 1, BaselineDataErr); + } + h1D_ReflCorr_BaselineSubtr = reinterpret_cast(h1D_ReflCorr->Clone("h1D_ReflCorr_BaselineSubtr")); + h1D_ReflCorr_BaselineSubtr->Sumw2(); + h1D_ReflCorr_BaselineSubtr->Add(hBaseline_Refl, -1.); + + TF1* fConstZero = new TF1("fConstZero", "[0]", 0., TMath::Pi()); + fConstZero->SetParameter(0, 0.); + fConstZero->SetLineColor(kMagenta); + fConstZero->SetLineStyle(9); + fConstZero->SetLineWidth(4); + fConstZero->SetTitle(""); + + TCanvas* cFinal_Reflected_BaselineSubtr = new TCanvas(Form("cFinal_Reflected_BaselineSubtr_%.0fto%.0f", PtHadMin, PtHadMax), Form("cFinal_Reflected_BaselineSubtr_%s_IntPools_PtAssoc%.0fto%.0f", fDmesonLabel.Data(), PtHadMin, PtHadMax), 100, 100, 1200, 700); + SetTH1HistoStyle(h1D_ReflCorr_BaselineSubtr, Form("%.0f < p_{T} < %.0f GeV/c", PtCandMin, PtCandMax), "#Delta#phi [rad]", "#frac{1}{N_{D}}#frac{dN^{assoc}}{d#Delta#phi} [rad^{-1}]", kFullCircle, kRed + 1, 1.6, kRed + 1, 3); + hBaseline_Refl->SetMarkerColor(kOrange); + hBaseline_Refl->SetMarkerStyle(kFullSquare); + hBaseline_Refl->SetLineColor(kOrange); + cFinal_Reflected_BaselineSubtr->cd(); + h1D_ReflCorr->SetMinimum(-0.8); + h1D_ReflCorr->SetStats(0); + hBaseline_Refl->SetStats(0); + h1D_ReflCorr->Draw(); + hBaseline_Refl->Draw("same"); + h1D_ReflCorr_BaselineSubtr->SetStats(0); + h1D_ReflCorr_BaselineSubtr->Draw("same"); // then keep just this + fConstZero->Draw("same"); + cFinal_Reflected_BaselineSubtr->SaveAs(Form("Output_CorrelationExtraction_%s_png/AzimCorrDistr_Reflected_BaselineSubtr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + cFinal_Reflected_BaselineSubtr->SaveAs(Form("Output_CorrelationExtraction_%s_Root/AzimCorrDistr_Reflected_BaselineSubtr_%s_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.root", codeName.Data(), fDmesonLabel.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + + fCorrectedCorrHisto_Reflected = reinterpret_cast(h1D_ReflCorr->Clone(Form("hCorrectedCorrReflected_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + fCorrectedCorrHisto_Reflected_BaselineSubtr = reinterpret_cast(h1D_ReflCorr_BaselineSubtr->Clone(Form("hCorrectedCorrReflected_BaselineSubtr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + return kTRUE; } @@ -394,26 +878,194 @@ Bool_t DhCorrelationExtraction::ReadInputInvMass() return kTRUE; } +Bool_t DhCorrelationExtraction::ReadInputFDSubtr() +{ + + fFileFDTemplate = TFile::Open(fFileFDTemplateName.Data()); + fFileFDPromptFrac = TFile::Open(fFileFDPromptFracName.Data()); + if (!fFileFDTemplate) { + std::cout << "[ERROR] File " << fFileFDTemplateName << " cannot be opened! check your file path!" << std::endl; + return kFALSE; + } + if (!fFileFDPromptFrac) { + std::cout << "[ERROR] File " << fFileFDPromptFracName << " cannot be opened! check your file path!" << std::endl; + return kFALSE; + } + + std::cout << "===================== " << std::endl; + std::cout << "Read inputs FD template" << std::endl; + std::cout << "TFile FD template = " << fFileFDTemplateName << std::endl; + std::cout << "TFile FD Prompt Frac = " << fFileFDPromptFracName << std::endl; + std::cout << "Histo FD template Prompt = " << fHistoFDTemplatePromptName << std::endl; + std::cout << "Histo FD template Non Prompt = " << fHistoFDTemplateNonPromptName << std::endl; + std::cout << "Histo FD Prompt Frac = " << fHistoFDPromptFracName << std::endl; + std::cout << "===================== " << std::endl; + std::cout << " " << std::endl; + + return kTRUE; +} + +Bool_t DhCorrelationExtraction::ReadInputSecondaryPartContamination() +{ + + fFileSecPart = TFile::Open(fFileSecPartName.Data()); + if (!fFileSecPart) { + std::cout << "[ERROR] File " << fFileSecPartName << " cannot be opened! check your file path!" << std::endl; + return kFALSE; + } + + fDirSecPart = reinterpret_cast(fFileSecPart->Get(fDirSecPartName.Data())); + + if (!fDirSecPart) { + std::cout << "[ERROR] Directory " << fDirSecPart << " cannot be opened! check your file path!" << std::endl; + return kFALSE; + } + + std::cout << "===================== " << std::endl; + std::cout << "Read inputs SE and ME" << std::endl; + std::cout << "TFile Sec. part. = " << fFileSecPartName << std::endl; + std::cout << "TDir Sec. part. = " << fDirSecPartName << std::endl; + std::cout << "===================== " << std::endl; + std::cout << " " << std::endl; + + return kTRUE; +} + +TH1D* DhCorrelationExtraction::EvaluateMCClosModulations(Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax) +{ + + TH1D* hModul = new TH1D(); + + fFilePromptMc = TFile::Open(fFilePromptMcRecName.Data()); + fFileNonPromptMc = TFile::Open(fFileNonPromptMcRecName.Data()); + + if (!fFilePromptMc) { + std::cout << "[ERROR] File prompt MC rec cannot be opened! check your file path!" << std::endl; + } + if (!fFileNonPromptMc) { + std::cout << "[ERROR] File non-prompt MC rec cannot be opened! check your file path!" << std::endl; + } + + // TODO: generalise this part + TH1D* hRecPrompt = reinterpret_cast(fFilePromptMc->Get(Form("h1D_Rec_iPtD%d_iPtAssoc%d", fBinPtCand, fBinPtHad))); + TH1D* hRecNonPrompt = reinterpret_cast(fFileNonPromptMc->Get(Form("h1D_Rec_iPtD%d_iPtAssoc%d", fBinPtCand, fBinPtHad))); + TH1D* hGenPrompt = reinterpret_cast(fFilePromptMc->Get(Form("h1D_Gen_iPtD%d_iPtAssoc%d", fBinPtCand, fBinPtHad))); + TH1D* hGenNonPrompt = reinterpret_cast(fFileNonPromptMc->Get(Form("h1D_Gen_iPtD%d_iPtAssoc%d", fBinPtCand, fBinPtHad))); + + printf("[INFO] Bin cand %d - Bin had %d \n", fBinPtCand, fBinPtHad); + + // hRecPrompt = ReflectCorrHistogram(hRecPrompt); + // hRecNonPrompt = ReflectCorrHistogram(hRecNonPrompt); + // hGenPrompt = ReflectCorrHistogram(hGenPrompt); + // hGenNonPrompt = ReflectCorrHistogram(hGenNonPrompt); + + hRecNonPrompt->Sumw2(); + hRecNonPrompt->Sumw2(); + hGenPrompt->Sumw2(); + hGenNonPrompt->Sumw2(); + + TH1D* hRatioNonPrompt = reinterpret_cast(hRecNonPrompt->Clone("hRatioNonPrompt")); + hRatioNonPrompt->Sumw2(); + hRatioNonPrompt->Divide(hRecNonPrompt, hGenNonPrompt, 1., 1., "B"); + hModul = reinterpret_cast(hRatioNonPrompt->Clone("hModul")); + + TF1* funFit = new TF1("funFit", "[0]", TMath::Pi() * 3. / 8., TMath::Pi() * 3 / 2); + hRatioNonPrompt->Fit(funFit, "R"); + Double_t fitVal = funFit->GetParameter(0); + + TCanvas* cRatio_MCClosure = new TCanvas(Form("cRatio_MCClosure_PtCand%.0fto%.0f_Pthad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax), Form("cRatio_MCClosure_PtCand%.0fto%.0f_Pthad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax), 100, 100, 1200, 700); + cRatio_MCClosure->cd(); + hRatioNonPrompt->GetYaxis()->SetRangeUser(0.2, 1.8); + hRatioNonPrompt->Draw(); + + Double_t FPrompt = GetFDPromptFrac(PtCandMin, PtCandMax, PtHadMin, PtHadMax); + Double_t relAmplC[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + Double_t relAmplB[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + Double_t recoKineVal[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + Double_t modul[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + for (int iBin = 0; iBin < hRatioNonPrompt->GetNbinsX(); iBin++) { + if (iBin > 1 && iBin < 13) { + recoKineVal[iBin - 2] = hRatioNonPrompt->GetBinContent(iBin + 1) - (fitVal - 1); + relAmplC[iBin - 2] = hRecPrompt->GetBinContent(iBin + 1) / (hRecPrompt->GetBinContent(iBin + 1) * FPrompt + hRecNonPrompt->GetBinContent(iBin + 1) * (1 - FPrompt)); + relAmplB[iBin - 2] = hRecNonPrompt->GetBinContent(iBin + 1) / (hRecPrompt->GetBinContent(iBin + 1) * FPrompt + hRecNonPrompt->GetBinContent(iBin + 1) * (1 - FPrompt)); + modul[iBin - 2] = relAmplC[iBin - 2] * FPrompt + relAmplB[iBin - 2] * (1 - FPrompt) / recoKineVal[iBin - 2]; + hModul->SetBinContent(iBin + 1, modul[iBin - 2]); + hModul->SetBinError(iBin + 1, 0.); + + printf("[INFO] Bin%d MODUL = %1.5f\t (Reco/Kine-fitVal = %1.4f, FPrompt = %1.3f, Ampl_ratio C,B = %1.4f, %1.4f)\n", iBin + 1, modul[iBin - 2], recoKineVal[iBin - 2], FPrompt, relAmplC[iBin - 2], relAmplB[iBin - 2]); + } else { + hModul->SetBinContent(iBin + 1, 1.); + hModul->SetBinError(iBin + 1, 0.); + } + } + + hModul->SetLineColor(kMagenta); + hModul->Draw("same"); + + cRatio_MCClosure->SaveAs(Form("Output_CorrelationExtraction_Thin2023_FullAnalysis_CentralPoints_png/Ratio_MCClosure_Canvas_PtCand%.0fto%.0f_PoolInt_PtAssoc%.0fto%.0f.png", PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + + return hModul; +} + TH2D* DhCorrelationExtraction::GetCorrelHisto(Int_t SEorME, Int_t SorSB, Int_t pool, Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax) { // TODO: Subtraction of softpion TH2D* h2D = new TH2D(); // pointer to be returned THnSparseD* hSparse = 0x0; - if (SEorME == kSE) { // Same Event if (SorSB == kSign) { hSparse = reinterpret_cast(fDirSE->Get(fSECorrelSignalRegionName.Data())); - } else { + } else if (!fSidebandDivided) { hSparse = reinterpret_cast(fDirSE->Get(fSECorrelSidebandsName.Data())); + } else if (fSidebandDivided) { + if (fUseSidebLeft && !fUseSidebRight) { + hSparse = reinterpret_cast(fDirSE->Get(fSECorrelSidebandLeftName.Data())); + } else if (!fUseSidebLeft && fUseSidebRight) { + hSparse = reinterpret_cast(fDirSE->Get(fSECorrelSidebandRightName.Data())); + } else if (fUseSidebLeft && fUseSidebRight) { + hSparse = reinterpret_cast(fDirSE->Get(fSECorrelSidebandLeftName.Data())); + hSparse->SetName("hSparse"); + THnSparseD* hSparseRightSideb = reinterpret_cast(fDirSE->Get(fSECorrelSidebandRightName.Data())); + hSparse->Add(hSparseRightSideb, 1.); + } } } else { // Mixed Event if (SorSB == kSign) { hSparse = reinterpret_cast(fDirME->Get(fMECorrelSignalRegionName.Data())); - } else { + } else if (!fSidebandDivided) { hSparse = reinterpret_cast(fDirME->Get(fMECorrelSidebandsName.Data())); + } else if (fSidebandDivided) { + if (fUseSidebLeft && !fUseSidebRight) { + hSparse = reinterpret_cast(fDirME->Get(fMECorrelSidebandLeftName.Data())); + } else if (!fUseSidebLeft && fUseSidebRight) { + hSparse = reinterpret_cast(fDirME->Get(fMECorrelSidebandRightName.Data())); + } else if (fUseSidebLeft && fUseSidebRight) { + hSparse = reinterpret_cast(fDirME->Get(fMECorrelSidebandLeftName.Data())); + hSparse->SetName("hSparse"); + THnSparseD* hSparseRightSideb = reinterpret_cast(fDirME->Get(fMECorrelSidebandRightName.Data())); + hSparse->Add(hSparseRightSideb, 1.); + } } } + /*else if (fSidebandDivided) { // Mixed Event + if (SorSB == kSign) { hSparse = reinterpret_cast fDirME -> Get(fMECorrelSignalRegionName.Data()); + } else if (!fSidebandDivided) { hSparse = reinterpret_cast fDirME -> Get(fMECorrelSidebandsName.Data()); + } else { + if (fUseSidebLeft && !fUseSidebRight) { + hSparse = reinterpret_cast fDirME -> Get(fMECorrelSidebandLeftName.Data()); + } else if (!fUseSidebLeft && fUseSidebRight) { + hSparse = reinterpret_cast fDirME -> Get(fMECorrelSidebandRightName.Data()); + } else if (fUseSidebLeft && fUseSidebRight) { + hSparse = reinterpret_cast fDirME -> Get(fMECorrelSidebandLeftName.Data()); + hSparse -> SetName("hSparse"); + THnSparseD *hSparseRightSideb = reinterpret_cast fDirME -> Get(fMECorrelSidebandRightName.Data()); + hSparse -> Add(hSparseRightSideb, 1.); + } + } + }*/ + Int_t binExtPtCandMin = (Int_t)hSparse->GetAxis(2)->FindBin(PtCandMin + 0.01); // axis2: ptCand, the 0.01 to avoid bin edges! Int_t binExtPtCandMax = (Int_t)hSparse->GetAxis(2)->FindBin(PtCandMax - 0.01); Int_t binExtPtHadMin = (Int_t)hSparse->GetAxis(3)->FindBin(PtHadMin + 0.01); // axis3: ptHad @@ -421,12 +1073,11 @@ TH2D* DhCorrelationExtraction::GetCorrelHisto(Int_t SEorME, Int_t SorSB, Int_t p Int_t binExtPoolMin; Int_t binExtPoolMax; if (fCorrectPoolsSeparately) { - binExtPoolMin = (Int_t)hSparse->GetAxis(4)->FindBin(pool + 1.01); // axis4: pool bin - binExtPoolMax = (Int_t)hSparse->GetAxis(4)->FindBin(pool + 1.99); + binExtPoolMin = (Int_t)hSparse->GetAxis(4)->FindBin(pool + 0.01); // axis4: pool bin + binExtPoolMax = (Int_t)hSparse->GetAxis(4)->FindBin(pool + 0.99); } else { // merge all pools in one binExtPoolMin = 1; binExtPoolMax = (Int_t)hSparse->GetAxis(4)->GetNbins(); - // cout << "binExtPoolMax:" << binExtPoolMax <GetAxis(1)->FindBin(fDeltaEtaMin + 0.0001); @@ -435,24 +1086,38 @@ TH2D* DhCorrelationExtraction::GetCorrelHisto(Int_t SEorME, Int_t SorSB, Int_t p binExtEtaMax = hSparse->GetAxis(1)->GetNbins(); if (binExtEtaMin < 1) binExtEtaMin = 1; - hSparse->GetAxis(1)->SetRange(binExtEtaMin, binExtEtaMax); // axis1: deltaEta hSparse->GetAxis(2)->SetRange(binExtPtCandMin, binExtPtCandMax); // axis2: ptCand hSparse->GetAxis(3)->SetRange(binExtPtHadMin, binExtPtHadMax); // axis3: ptHad - // hSparse -> GetAxis(4) -> SetRange(binExtPoolMin, binExtPoolMax); // axis4: pool bin - - h2D = reinterpret_cast(hSparse->Projection(0, 1)); // axis0: deltaPhi, axis1: deltaEta - if (SEorME == kSE) { // Same Event + hSparse->GetAxis(4)->SetRange(binExtPoolMin, binExtPoolMax); // axis4: pool bin + h2D = reinterpret_cast(hSparse->Projection(0, 1)); // axis0: deltaPhi, axis1: deltaEta + if (SEorME == kSE) { // Same Event if (SorSB == kSign) { h2D->SetName(Form("hCorr_SE_Sig_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); - } else { + } else if (!fSidebandDivided) { h2D->SetName(Form("hCorr_SE_Sideb_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } else { + if (fUseSidebLeft && !fUseSidebRight) { + h2D->SetName(Form("hCorr_SE_Sideb_Left_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } else if (!fUseSidebLeft && fUseSidebRight) { + h2D->SetName(Form("hCorr_SE_Sideb_Right_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } else if (fUseSidebLeft && fUseSidebRight) { + h2D->SetName(Form("hCorr_SE_Sideb_LeftAndRight_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } } } else { // Mixed Event if (SorSB == kSign) { h2D->SetName(Form("hCorr_ME_Sig_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } else if (!fSidebandDivided) { + h2D->SetName(Form("hCorr_SE_Sideb_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); } else { - h2D->SetName(Form("hCorr_ME_Sideb_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + if (fUseSidebLeft && !fUseSidebRight) { + h2D->SetName(Form("hCorr_ME_Sideb_Left_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } else if (!fUseSidebLeft && fUseSidebRight) { + h2D->SetName(Form("hCorr_ME_Sideb_Right_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } else if (fUseSidebLeft && fUseSidebRight) { + h2D->SetName(Form("hCorr_ME_Sideb_LeftAndRight_2D_PtCandBin%d_PtHadBin%d_iPool%d", binExtPtCandMin, binExtPtHadMin, pool)); + } } } @@ -466,6 +1131,8 @@ void DhCorrelationExtraction::GetSignalAndBackgroundForNorm(Double_t PtCandMin, TH1F* hMassFitSgnYield = reinterpret_cast(fFileMass->Get(fMassHistoNameSgn.Data())); TH1F* hMassFitBkgYield = reinterpret_cast(fFileMass->Get(fMassHistoNameBkg.Data())); TH1F* hMassFitSBsYield = reinterpret_cast(fFileMass->Get(fMassHistoNameSBs.Data())); + TH1F* hMassFitSBLYield = reinterpret_cast(fFileMass->Get("hBackgroundSidebandLeft")); + TH1F* hMassFitSBRYield = reinterpret_cast(fFileMass->Get("hBackgroundSidebandRight")); Int_t PtCandBin = hMassFitSgnYield->FindBin(PtCandMin + 0.01); if (PtCandBin != hMassFitSgnYield->FindBin(PtCandMax - 0.01)) @@ -474,6 +1141,8 @@ void DhCorrelationExtraction::GetSignalAndBackgroundForNorm(Double_t PtCandMin, Float_t SgnYield = hMassFitSgnYield->GetBinContent(PtCandBin); Float_t BkgYield = hMassFitBkgYield->GetBinContent(PtCandBin); Float_t SBsYield = hMassFitSBsYield->GetBinContent(PtCandBin); + Float_t SBLYield = hMassFitSBLYield->GetBinContent(PtCandBin); + Float_t SBRYield = hMassFitSBRYield->GetBinContent(PtCandBin); std::cout << "================================= " << std::endl; std::cout << "Getting invariant mass parameters " << std::endl; @@ -481,20 +1150,193 @@ void DhCorrelationExtraction::GetSignalAndBackgroundForNorm(Double_t PtCandMin, std::cout << "Signal yield = " << SgnYield << std::endl; std::cout << "Bkg yield = " << BkgYield << std::endl; std::cout << "Sideband yield = " << SBsYield << std::endl; + std::cout << "Sideband left yield = " << SBLYield << std::endl; + std::cout << "Sideband right yield = " << SBRYield << std::endl; std::cout << "================================= " << std::endl; std::cout << " " << std::endl; SetSignalYieldforNorm(SgnYield); - SetBkgScaleFactor(BkgYield / SBsYield); + SetBkgYield(BkgYield); + if (fUseSidebLeft && fUseSidebRight) { + SetBkgScaleFactor(BkgYield / SBsYield); + SetSBYield(SBsYield); + } else if (fUseSidebLeft && !fUseSidebRight) { + SetBkgScaleFactor(BkgYield / SBLYield); + SetSBYield(SBLYield); + } else if (!fUseSidebLeft && fUseSidebRight) { + SetBkgScaleFactor(BkgYield / SBRYield); + SetSBYield(SBRYield); + } return; } +TH2D* DhCorrelationExtraction::GetFDTemplateHisto(Int_t PromptOrFD, Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax) +{ + + TH2D* h2D = new TH2D(); // pointer to be returned + + if (PromptOrFD == kPrompt) { + h2D = reinterpret_cast(fFileFDTemplate->Get(Form("%s%.0f_%.0f_ptassoc%.0f_%.0f", fHistoFDTemplatePromptName.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + } else { + h2D = reinterpret_cast(fFileFDTemplate->Get(Form("%s%.0f_%.0f_ptassoc%.0f_%.0f", fHistoFDTemplateNonPromptName.Data(), PtCandMin, PtCandMax, PtHadMin, PtHadMax))); + } + + Int_t binExtEtaMin = (Int_t)h2D->GetXaxis()->FindBin(fDeltaEtaMin + 0.000001); + Int_t binExtEtaMax = (Int_t)h2D->GetXaxis()->FindBin(fDeltaEtaMax - 0.000001); + if (binExtEtaMax > h2D->GetXaxis()->GetNbins()) + binExtEtaMax = h2D->GetXaxis()->GetNbins(); + if (binExtEtaMin < 1) + binExtEtaMin = 1; + + h2D->GetXaxis()->SetRange(binExtEtaMin, binExtEtaMax); + if (PromptOrFD == kPrompt) { + h2D->SetName(Form("hFDTemplatePrompt_2D_PtCand%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } else { + h2D->SetName(Form("hFDTemplateNonPrompt_2D_PtCand%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } + h2D->GetYaxis()->SetTitle("#Delta#phi (rad)"); + h2D->GetXaxis()->SetTitle("#Delta#eta"); + + return h2D; +} + +TH1D* DhCorrelationExtraction::GetCorrelHistoSecondaryPart(Int_t PartType, Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax) +{ + + TH1D* h1D = new TH1D(); // pointer to be returned + + THnSparseD* hSparse = 0x0; + + if (PartType == kPrimaryPart) { // primary particles + hSparse = reinterpret_cast(fDirSecPart->Get(fHistoPrimaryPartName.Data())); + } else { // all selected particles + hSparse = reinterpret_cast(fDirSecPart->Get(fHistoAllPartName.Data())); + } + Int_t binExtPtCandMin = (Int_t)hSparse->GetAxis(2)->FindBin(PtCandMin + 0.01); // axis2: ptCand, the 0.01 to avoid bin edges! + Int_t binExtPtCandMax = (Int_t)hSparse->GetAxis(2)->FindBin(PtCandMax - 0.01); + Int_t binExtPtHadMin = (Int_t)hSparse->GetAxis(3)->FindBin(PtHadMin + 0.01); // axis3: ptHad + Int_t binExtPtHadMax = (Int_t)hSparse->GetAxis(3)->FindBin(PtHadMax - 0.01); + Int_t binExtPoolMin; + Int_t binExtPoolMax; + if (PartType == kAllPart) { + binExtPoolMin = 1; + binExtPoolMax = (Int_t)hSparse->GetAxis(4)->GetNbins(); + } + // possibility to select a certain eta region + Int_t binExtEtaMin = (Int_t)hSparse->GetAxis(1)->FindBin(fDeltaEtaMin + 0.0001); + Int_t binExtEtaMax = (Int_t)hSparse->GetAxis(1)->FindBin(fDeltaEtaMax - 0.0001); + if (binExtEtaMax > hSparse->GetAxis(1)->GetNbins()) + binExtEtaMax = hSparse->GetAxis(1)->GetNbins(); + if (binExtEtaMin < 1) + binExtEtaMin = 1; + + hSparse->GetAxis(1)->SetRange(binExtEtaMin, binExtEtaMax); // axis1: deltaEta + hSparse->GetAxis(2)->SetRange(binExtPtCandMin, binExtPtCandMax); // axis2: ptCand + hSparse->GetAxis(3)->SetRange(binExtPtHadMin, binExtPtHadMax); // axis3: ptHad + if (PartType == kAllPart) { + hSparse->GetAxis(4)->SetRange(binExtPoolMin, binExtPoolMax); // axis4: pool bin + } + + h1D = reinterpret_cast(hSparse->Projection(0)); // axis0: deltaPhi + if (PartType == kPrimaryPart) { // primary particles + h1D->SetName(Form("hPrimaryPartCorr_PtD%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } else { // all selected particles + h1D->SetName(Form("hAllPartCorr_PtD%.0fto%.0f_PtHad%.0fto%.0f", PtCandMin, PtCandMax, PtHadMin, PtHadMax)); + } + + return h1D; +} + +TH1D* DhCorrelationExtraction::ReflectCorrHistogram(TH1D*& histo) +{ + + // nBinsPhi must be a multple of 4 in order to reflect correcty the histogram + Int_t nBinsPhi = histo->GetNbinsX(); + Int_t nBinsPhiRefl = nBinsPhi / 2; + Int_t bin0Phi = nBinsPhi / 4 + 1; + Int_t binPiPhi = 3 * nBinsPhi / 4; + + TH1D* h1D = new TH1D("h1D_Reflected", "", nBinsPhiRefl, 0., TMath::Pi()); // pointer to be returned + h1D->Sumw2(); + // TH1D* h1D = reinterpret_cast histo -> Clone("h1D_Reflected"); + // h1D -> GetXaxis() -> SetRange(bin0Phi, binPiPhi); + + // reflection + Double_t reflectedContent, reflectedContentError; + for (int iBin = 0; iBin < nBinsPhiRefl / 2; iBin++) { + reflectedContent = (histo->GetBinContent(bin0Phi - iBin - 1) + histo->GetBinContent(bin0Phi + iBin)) / 2; + std::cout << "[INFO] Paired bins for reflection: " << bin0Phi - iBin - 1 << " - " << bin0Phi + iBin << std::endl; + std::cout << "[INFO] Bin filled: " << iBin + 1 << std::endl; + reflectedContentError = 0.5 * TMath::Sqrt(TMath::Power(histo->GetBinError(iBin + 1), 2) + TMath::Power(histo->GetBinError(bin0Phi + iBin), 2)); + h1D->SetBinContent(iBin + 1, reflectedContent); + h1D->SetBinError(iBin + 1, reflectedContentError); + } + for (int iBin = nBinsPhiRefl / 2; iBin < nBinsPhiRefl; iBin++) { + reflectedContent = (histo->GetBinContent(bin0Phi + iBin) + histo->GetBinContent(binPiPhi + 2 * bin0Phi - iBin - 2)) / 2; + reflectedContentError = 0.5 * TMath::Sqrt(TMath::Power(histo->GetBinError(bin0Phi + iBin), 2) + TMath::Power(histo->GetBinError(binPiPhi + 2 * bin0Phi - iBin - 2), 2)); + std::cout << "[INFO] Paired bins for reflection: " << bin0Phi + iBin << " - " << binPiPhi + 2 * bin0Phi - iBin - 2 << std::endl; + std::cout << "[INFO] Bin filled: " << iBin + 1 << std::endl; + h1D->SetBinContent(iBin + 1, reflectedContent); + h1D->SetBinError(iBin + 1, reflectedContentError); + } + + return h1D; +} + +TH1D* DhCorrelationExtraction::ReflectHistoRun2(TH1D* h, Double_t scale) +{ + + TH1D* h2 = new TH1D(Form("%sReflected", h->GetName()), Form("%sReflected", h->GetName()), h->GetNbinsX() / 2., 0., TMath::Pi()); + for (Int_t j = 1; j <= h->GetNbinsX(); j++) { + Double_t x = h->GetBinCenter(j); + Double_t y0 = h->GetBinContent(j); + Double_t ey0 = h->GetBinError(j); + Int_t j2; + if (x > 0 && x < TMath::Pi()) { + j2 = h2->FindBin(x); + } else if (x < 0) { + j2 = h2->FindBin(-1. * x); + } else if (x > TMath::Pi()) { + j2 = h2->FindBin(2. * TMath::Pi() - x); + } else { + printf("Point %d excluded \n", j); + continue; + } + Double_t y = h2->GetBinContent(j2); + Double_t ey = h2->GetBinError(j2); + h2->SetBinContent(j2, (y + y0)); + h2->SetBinError(j2, TMath::Sqrt(ey0 * ey0 + ey * ey)); + } + h2->Scale(scale); + + return h2; +} + +Double_t DhCorrelationExtraction::GetFDPromptFrac(Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax) +{ + + TH1D* h1D = new TH1D(); + h1D = reinterpret_cast(fFileFDPromptFrac->Get(fHistoFDPromptFracName.Data())); + + Int_t binPtCandMin = h1D->GetXaxis()->FindBin(PtCandMin + 0.01); + Int_t binPtCandMax = h1D->GetXaxis()->FindBin(PtCandMax - 0.01); + Double_t PromptFraction; + if (binPtCandMin == binPtCandMax) { + PromptFraction = h1D->GetBinContent(binPtCandMin); + } else { + std::cout << "[ERROR] Different bin obtained from PtCandMin and PtCandMax"; + return 0.; + } + + return PromptFraction; +} + void DhCorrelationExtraction::NormalizeMEplot(TH2D*& histoME, TH2D*& histoMEsoftPi) { - Double_t bin0phi = histoME->GetYaxis()->FindBin(0.); - Double_t bin0eta = histoME->GetXaxis()->FindBin(0.); + Int_t bin0phi = histoME->GetYaxis()->FindBin(0.); + Int_t bin0eta = histoME->GetXaxis()->FindBin(0.); // evaluate the normalization (from ALL tracks, including possible fake softpions) -> **histoME indeed includes bin1+bin2 of THnSparse, i.e. all the tracks** Double_t factorNorm = 0; @@ -506,9 +1348,13 @@ void DhCorrelationExtraction::NormalizeMEplot(TH2D*& histoME, TH2D*& histoMEsoft } factorNorm /= 4.; - if (fSubtractSoftPiME) { + std::cout << "bin 0 phi: " << bin0phi << std::endl; + std::cout << "bin 0 eta: " << bin0eta << std::endl; + std::cout << "Factor norm. ME: " << factorNorm << std::endl; + std::cout << "Bin content (0,0) ME: " << histoME->GetBinContent(bin0eta, bin0phi) << std::endl; + + if (fSubtractSoftPiME) histoME->Add(histoMEsoftPi, -1); // remove the tracks compatible with soft pion (if requested) - } // apply the normalization histoME->Scale(1. / factorNorm); @@ -516,6 +1362,144 @@ void DhCorrelationExtraction::NormalizeMEplot(TH2D*& histoME, TH2D*& histoMEsoft return; } +Double_t DhCorrelationExtraction::CalculateBaseline(TH1D*& histo, Bool_t totalRange, Bool_t reflected) +{ + + // total range = 2*Pi + // half range = Pi , for histogram reflected under symmetric assumption + + Double_t baseline, errBaseline; + Int_t nBinsPhi = histo->GetNbinsX(); + Int_t binPhiHalf = nBinsPhi / 2; + Int_t binPhiHalfMinus1 = nBinsPhi / 2 - 1; + Int_t binPhiHalfPlus1 = nBinsPhi / 2 + 1; + Int_t binPhiHalfPlus2 = nBinsPhi / 2 + 1; + + if (totalRange) { + // baseline evaluated considering: the two first points, the last two points and four points in the middle (corresponding to the outer points) + if (nBinsPhi >= 32) { + baseline = + ((histo->GetBinContent(1)) * (1. / TMath::Power(histo->GetBinError(1), 2)) + + (histo->GetBinContent(2)) * (1. / TMath::Power(histo->GetBinError(2), 2)) + + (histo->GetBinContent(binPhiHalfMinus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (histo->GetBinContent(binPhiHalfPlus2)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2)) + + (histo->GetBinContent(nBinsPhi - 1)) * (1. / TMath::Power(histo->GetBinError(nBinsPhi - 1), 2)) + + (histo->GetBinContent(nBinsPhi)) * (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))) / + ((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(2), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi - 1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } else { + baseline = + ((histo->GetBinContent(1)) * (1. / TMath::Power(histo->GetBinError(1), 2)) + + (histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (histo->GetBinContent(nBinsPhi)) * (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))) / + ((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } + } else { + if (reflected) { + baseline = + ((histo->GetBinContent(binPhiHalfMinus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (histo->GetBinContent(binPhiHalfPlus2)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))) / + ((1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))); + } else { + // baseline evaluated using the 4 middle points in the transverese region + if (nBinsPhi >= 32) { + baseline = + ((histo->GetBinContent(binPhiHalfMinus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (histo->GetBinContent(binPhiHalfPlus2)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))) / + ((1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))); + } else { + baseline = + ((histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2))) / + ((1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2))); + } + } + } + + return baseline; +} + +Double_t DhCorrelationExtraction::CalculateBaselineError(TH1D*& histo, Bool_t totalRange, Bool_t reflected) +{ + + // total range = 2*Pi + // half range = Pi , for histogram reflected under symmetric assumption + + Double_t errBaseline; + Int_t nBinsPhi = histo->GetNbinsX(); + Int_t binPhiHalf = nBinsPhi / 2; + Int_t binPhiHalfMinus1 = nBinsPhi / 2 - 1; + Int_t binPhiHalfPlus1 = nBinsPhi / 2 + 1; + Int_t binPhiHalfPlus2 = nBinsPhi / 2 + 1; + + if (totalRange) { + // baseline evaluated considering: the two first points, the last two points and four points in the middle (corresponding to the outer points) + if (nBinsPhi >= 32) { + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(2), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi - 1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } else { // fon nBinsPhi = 16 (rebin 4) + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } + } else { + // baseline evaluated using the 4 middle points in the transverese region + if (reflected) { + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))); + } else { + if (nBinsPhi >= 32) { + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))); + } else { + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2))); + } + } + } + + return errBaseline; +} + void DhCorrelationExtraction::SetTH1HistoStyle(TH1D*& histo, TString hTitle, TString hXaxisTitle, TString hYaxisTitle, Style_t markerStyle, Color_t markerColor, Double_t markerSize, Color_t lineColor, Int_t lineWidth, Float_t hTitleXaxisOffset, Float_t hTitleYaxisOffset, diff --git a/PWGHF/HFC/Macros/DhCorrelationExtraction.h b/PWGHF/HFC/Macros/DhCorrelationExtraction.h index 73b07c61e83..52508ea3bb9 100644 --- a/PWGHF/HFC/Macros/DhCorrelationExtraction.h +++ b/PWGHF/HFC/Macros/DhCorrelationExtraction.h @@ -17,8 +17,6 @@ #ifndef PWGHF_HFC_MACROS_DHCORRELATIONEXTRACTION_H_ #define PWGHF_HFC_MACROS_DHCORRELATIONEXTRACTION_H_ -#include -#include #include #include #include @@ -28,6 +26,9 @@ #include #include +#include +#include + class DhCorrelationExtraction : public TObject { @@ -40,6 +41,10 @@ class DhCorrelationExtraction : public TObject kME }; enum selectInvMassRegion { kSign, kSideb }; + enum selectDmesonOrigin { kPrompt, + kFD }; + enum selectParticleType { kPrimaryPart, + kAllPart }; DhCorrelationExtraction(); // default constructor DhCorrelationExtraction(const DhCorrelationExtraction& source); @@ -51,15 +56,36 @@ class DhCorrelationExtraction : public TObject void SetInputFilenameMass(TString filenameMass) { fFileNameMass = filenameMass; } void SetInputFilenameSE(TString filenameSE) { fFileNameSE = filenameSE; } void SetInputFilenameME(TString filenameME) { fFileNameME = filenameME; } + void SetInputFilenameSecPart(TString filenameSecPart) { fFileSecPartName = filenameSecPart; } + void SetInputFilenameBiasBtoD(TString filenamePromptMcRec, TString filenameNonPromptMcRec) + { + fFilePromptMcRecName = filenamePromptMcRec; + fFileNonPromptMcRecName = filenameNonPromptMcRec; + } void SetDirNameSE(TString dirNameSE) { fDirNameSE = dirNameSE; } void SetDirNameME(TString dirNameME) { fDirNameME = dirNameME; } + void SetDirNameSecPart(TString dirNameSecPart) { fDirSecPartName = dirNameSecPart; } void SetMassHistoNameSgn(TString massHistoNameSgn) { fMassHistoNameSgn = massHistoNameSgn; } void SetMassHistoNameBkg(TString massHistoNameBkg) { fMassHistoNameBkg = massHistoNameBkg; } void SetMassHistoNameSBs(TString massHistoNameSBs) { fMassHistoNameSBs = massHistoNameSBs; } void SetSECorrelHistoSignalName(TString correlNameSigSE) { fSECorrelSignalRegionName = correlNameSigSE; } void SetSECorrelHistoSidebandName(TString correlNameSbSE) { fSECorrelSidebandsName = correlNameSbSE; } + void SetSECorrelHistoSidebandLeftName(TString correlNameSbSE) { fSECorrelSidebandLeftName = correlNameSbSE; } + void SetSECorrelHistoSidebandRightName(TString correlNameSbSE) { fSECorrelSidebandRightName = correlNameSbSE; } void SetMECorrelHistoSignalName(TString correlNameSigME) { fMECorrelSignalRegionName = correlNameSigME; } void SetMECorrelHistoSidebandName(TString correlNameSbME) { fMECorrelSidebandsName = correlNameSbME; } + void SetMECorrelHistoSidebandLeftName(TString correlNameSbME) { fMECorrelSidebandLeftName = correlNameSbME; } + void SetMECorrelHistoSidebandRightName(TString correlNameSbME) { fMECorrelSidebandRightName = correlNameSbME; } + void SetHistoSecPartName(TString histoPrimaryPartName, TString histoAllPartName) + { + fHistoPrimaryPartName = histoPrimaryPartName; + fHistoAllPartName = histoAllPartName; + } + void SetInputFilenameFDTemplate(TString filenameFDTemplate) { fFileFDTemplateName = filenameFDTemplate; } + void SetInputFilenameFDPromptFrac(TString filenameFDPromptFrac) { fFileFDPromptFracName = filenameFDPromptFrac; } + void SetInputHistoNameFDTemplatePrompt(TString hNameFDTemplatePrompt) { fHistoFDTemplatePromptName = hNameFDTemplatePrompt; } + void SetInputHistoNameFDTemplateNonPrompt(TString hNameFDTemplateNonPrompt) { fHistoFDTemplateNonPromptName = hNameFDTemplateNonPrompt; } + void SetInputHistoNameFDPromptFrac(TString hNameFDPromptFrac) { fHistoFDPromptFracName = hNameFDPromptFrac; } // Input conditions: PtCand, PtHad, PoolBins void SetNpools(Int_t npools) { fNpools = npools; } @@ -72,66 +98,136 @@ class DhCorrelationExtraction : public TObject void SetSubtractSoftPiInMEdistr(Bool_t subtractSoftPiME) { fSubtractSoftPiME = subtractSoftPiME; } void SetBkgScaleFactor(Double_t scaleFactor) { fBkgScaleFactor = scaleFactor; } void SetSignalYieldforNorm(Double_t sgnYield) { fSgnYieldNorm = sgnYield; } + void SetBkgYield(Double_t bkgYield) { fBkgYield = bkgYield; } + void SetSBYield(Double_t SBYield) { fSBYield = SBYield; } void SetRebin2DcorrelHisto(Int_t rebinDeltaEta, Int_t rebinDeltaPhi) { - fRebin2Dhisto = kTRUE; fRebinAxisDeltaEta = rebinDeltaEta; fRebinAxisDeltaPhi = rebinDeltaPhi; } + void SetRebinOptions(Bool_t rebinAngCorr, Bool_t rebinFDCorr, Bool_t rebinSecPart) + { + fRebinAngCorr = rebinAngCorr; + fRebinFDCorr = rebinFDCorr; + fRebinSecPart = rebinSecPart; + } void GetSignalAndBackgroundForNorm(Double_t PtCandMin, Double_t PtCandMax); void NormalizeMEplot(TH2D*& histoME, TH2D*& histoMEsoftPi); void SetDebugLevel(Int_t debug) { fDebug = debug; } + void SetDividedSidebands(Bool_t dividedSideb, Bool_t useSidebLeft, Bool_t useSidebRight) + { + fSidebandDivided = dividedSideb; + fUseSidebLeft = useSidebLeft; + fUseSidebRight = useSidebRight; + } + void SetFDSubtraction(Bool_t subtractFD) { fFDsubtraction = subtractFD; } + void SetSecPartContamination(Bool_t secPartContamination) { fSecPartContamination = secPartContamination; } + void SetCorrBiasBtoD(Bool_t corrbiasBtoD) { fCorrBiasBtoD = corrbiasBtoD; } + void SetBinCandAndHad(Int_t binCand, Int_t binHad) + { + fBinPtCand = binCand; + fBinPtHad = binHad; + } /// Analysis methods TH2D* GetCorrelHisto(Int_t SEorME, Int_t SorSB, Int_t pool, Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax); + TH2D* GetFDTemplateHisto(Int_t PromptOrFD, Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax); + TH1D* GetCorrelHistoSecondaryPart(Int_t PrimaryPart, Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax); + TH1D* ReflectCorrHistogram(TH1D*& histo); + TH1D* ReflectHistoRun2(TH1D* h, Double_t scale); + TH1D* EvaluateMCClosModulations(Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax); + Double_t GetFDPromptFrac(Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax); + Double_t CalculateBaseline(TH1D*& histo, Bool_t totalRange = kTRUE, Bool_t reflected = kFALSE); + Double_t CalculateBaselineError(TH1D*& histo, Bool_t totalRange = kTRUE, Bool_t reflected = kFALSE); Bool_t ReadInputSEandME(); Bool_t ReadInputInvMass(); + Bool_t ReadInputFDSubtr(); + Bool_t ReadInputSecondaryPartContamination(); Bool_t ExtractCorrelations(Double_t PtCandMin, Double_t PtCandMax, Double_t PtHadMin, Double_t PtHadMax, TString codeName); TH1D* GetCorrectedCorrHisto() { return fCorrectedCorrHisto; } + TH1D* GetCorrectedCorrHisto_BaselineSubtr() { return fCorrectedCorrHisto_BaselineSubtr; } + TH1D* GetCorrectedCorrHisto_Reflected() { return fCorrectedCorrHisto_Reflected; } + TH1D* GetCorrectedCorrHisto_Reflected_BaselineSubtr() { return fCorrectedCorrHisto_Reflected_BaselineSubtr; } /// Histogram style void SetTH1HistoStyle(TH1D*& histo, TString hTitle, TString hXaxisTitle, TString hYaxisTitle, Style_t markerStyle = kFullCircle, Color_t markerColor = kRed + 1, Double_t markerSize = 1.4, Color_t lineColor = kRed + 1, Int_t lineWidth = 3, Float_t hTitleXaxisOffset = 1.0, Float_t hTitleYaxisOffset = 1.0, Float_t hTitleXaxisSize = 0.060, Float_t hTitleYaxisSize = 0.060, Float_t hLabelXaxisSize = 0.060, Float_t hLabelYaxisSize = 0.060, Bool_t centerXaxisTitle = false, Bool_t centerYaxisTitle = false); void SetTH2HistoStyle(TH2D*& histo, TString hTitle, TString hXaxisTitle, TString hYaxisTitle, TString hZaxisTitle, Float_t hTitleXaxisOffset = 1.8, Float_t hTitleYaxisOffset = 1.8, Float_t hTitleZaxisOffset = 1.2, Float_t hTitleXaxisSize = 0.060, Float_t hTitleYaxisSize = 0.060, Float_t hTitleZaxisSize = 0.060, Float_t hLabelXaxisSize = 0.060, Float_t hLabelYaxisSize = 0.060, Float_t hLabelZaxisSize = 0.060, Bool_t centerXaxisTitle = true, Bool_t centerYaxisTitle = true); private: - TFile* fFileMass; // file containing the mass histograms - TFile* fFileSE; // file containing the Same Event (SE) output - TFile* fFileME; // file containing the Mixed Event (ME) output - - TDirectoryFile* fDirMass; // TDirectory for mass histos - TDirectoryFile* fDirSE; // TDirectory for SE info - TDirectoryFile* fDirME; // TDirectory for ME info - - TH1D* fCorrectedCorrHisto; // Corrected correlation histogram - - DmesonSpecie fDmesonSpecies; // D meson specie - TString fDmesonLabel; // D meson label - TString fFileNameMass; // File cntaining inv. mass histograms - TString fFileNameSE; // File contaning Same Event (SE) output - TString fFileNameME; // File contaning Mixed Event (ME) output - TString fDirNameSE; // Directory in the file containing SE output - TString fDirNameME; // Directory in the file containing ME output - TString fMassHistoNameSgn; // Inv. mass histo name signal yield - TString fMassHistoNameBkg; // Inv. mass histo name background yield - TString fMassHistoNameSBs; // Inv. mass histo name sideband yield - TString fSECorrelSignalRegionName; // THnSparse name containing SE output for signal region - TString fSECorrelSidebandsName; // THnSparse name containing SE output for sideband region - TString fMECorrelSignalRegionName; // THnSparse name containing ME output for signal region - TString fMECorrelSidebandsName; // THnSparse name containing ME output for sideband region - - Int_t fNpools; // number of pools used for the ME correction - Int_t fRebinAxisDeltaEta; // rebin deltaEta axis - Int_t fRebinAxisDeltaPhi; // rebin deltaPhi axis - Int_t fDebug; // debug level - - Double_t fDeltaEtaMin; // deltaEta min value - Double_t fDeltaEtaMax; // deltaEta max value + TFile* fFileMass; // File containing the mass histograms + TFile* fFileSE; // File containing the Same Event (SE) output + TFile* fFileME; // File containing the Mixed Event (ME) output + TFile* fFileFDTemplate; // File containing FD angular correlation templates + TFile* fFileFDPromptFrac; // File containing prompt fraction (used fo FD subtraction) + TFile* fFileSecPart; // File containing secondary particle contaminaion teplates + TFile* fFilePromptMc; // File containing prompt ratio taken from MC Closure test study (to use for B to D bias correction) + TFile* fFileNonPromptMc; // File containing non-prompt ratio taken from MC Closure test study (to use for B to D bias correction) + + TDirectoryFile* fDirMass; // TDirectory for mass histos + TDirectoryFile* fDirSE; // TDirectory for SE info + TDirectoryFile* fDirME; // TDirectory for ME info + TDirectoryFile* fDirSecPart; // TDirectory for seondary particle correction + + TH1D* fCorrectedCorrHisto; // Corrected correlation histogram + TH1D* fCorrectedCorrHisto_BaselineSubtr; // Corrected correlation histogram with baseline subtracion + TH1D* fCorrectedCorrHisto_Reflected; // Corrected correlation histogram relected in azimuth + TH1D* fCorrectedCorrHisto_Reflected_BaselineSubtr; // Corrected correlation histogram reflected in azimuth with baseline subtraction + + DmesonSpecie fDmesonSpecies; // D meson specie + TString fDmesonLabel; // D meson label + TString fFileNameMass; // File name containing inv. mass histograms + TString fFileNameSE; // File name contaning Same Event (SE) output + TString fFileNameME; // File name contaning Mixed Event (ME) output + TString fFileSecPartName; // File name contaning secondary particle correction output + TString fFileFDTemplateName; // File name contaning FD angular correlation templates + TString fFileFDPromptFracName; // File name contaning prompt fraction (used for FD subtraction) + TString fFilePromptMcRecName; // File name contaning prompt angular correlation (used for B to d bias correction) + TString fFileNonPromptMcRecName; // File name contaning non-prompt angular correlation (used for B to d bias correction) + TString fDirNameSE; // Directory in the file containing SE output + TString fDirNameME; // Directory in the file containing ME output + TString fDirSecPartName; // Directory in the file containing secondary particle correction output + TString fMassHistoNameSgn; // Inv. mass histo name signal yield + TString fMassHistoNameBkg; // Inv. mass histo name background yield + TString fMassHistoNameSBs; // Inv. mass histo name sideband yield + TString fSECorrelSignalRegionName; // THnSparse name containing SE output for signal region + TString fSECorrelSidebandsName; // THnSparse name containing SE output for sideband region + TString fSECorrelSidebandLeftName; // THnSparse name containing SE output for sideband left region + TString fSECorrelSidebandRightName; // THnSparse name containing SE output for sideband right region + TString fMECorrelSignalRegionName; // THnSparse name containing ME output for signal region + TString fMECorrelSidebandsName; // THnSparse name containing ME output for sideband regions + TString fMECorrelSidebandLeftName; // THnSparse name containing ME output for sideband left region + TString fMECorrelSidebandRightName; // THnSparse name containing ME output for sideband right region + TString fHistoFDTemplatePromptName; // Prompt angular correlation histogram name + TString fHistoFDTemplateNonPromptName; // FD angular correlation histogram name + TString fHistoFDPromptFracName; // Prompt fraction histogram name + TString fHistoPrimaryPartName; // Primary particle histogram (to be used for secondary particle contamination correction) + TString fHistoAllPartName; // All particle histogram (to be used for secondary particle contamination correction) + + Int_t fNpools; // Number of pools used for the ME correction + Int_t fRebinAxisDeltaEta; // Rebin deltaEta axis value + Int_t fRebinAxisDeltaPhi; // Rebin deltaPhi axis value + Int_t fDebug; // Debug level + Int_t fBinPtCand; // Pt bin of the candidate + Int_t fBinPtHad; // Pt bin of the hadron + + Double_t fDeltaEtaMin; // DeltaEta min value + Double_t fDeltaEtaMax; // DeltaEta max value Double_t fBkgScaleFactor; // Bkg/SB factor to scale correlation plots obtained in the sideband region Double_t fSgnYieldNorm; // Signal yield (used for normalize correlation plots after bkg subtraction) + Double_t fBkgYield; // Bkg yield under signal peak region + Double_t fSBYield; // Sideband yield Bool_t fCorrectPoolsSeparately; // Possibility to do the ME correction pool-by-pool (kTRUE) or merging all pools (kFALSE) Bool_t fSubtractSoftPiME; // Soft pion subtraction (for D0 case) - Bool_t fRebin2Dhisto; // Flag to rebin the 2D correlation plots + Bool_t fRebinAngCorr; // Rebin angular correlaion distributons (SE and ME) + Bool_t fRebinFDCorr; // Rebin angular correlaion distributon templates used for FD correction (theory driven) + Bool_t fRebinSecPart; // Rebin angular correlaion distributon templates used for secodary particle contamination correction + Bool_t fSidebandDivided; // To be set to TRUE if two sideband corrlaion histograms are passed inteh config file + Bool_t fUseSidebLeft; // To be set to TRUE if only sideband left is used for the bkg correction + Bool_t fUseSidebRight; // To be set to TRUE if only sideband right is used for the bkg correction + Bool_t fFDsubtraction; // Enable feed-down (FD) correction + Bool_t fSecPartContamination; // Enable seconday particle contamination correction + Bool_t fCorrBiasBtoD; // Enable bias B to D correction (perfrmed with angular correlaion templates taken from MC simulaions) }; #endif // PWGHF_HFC_MACROS_DHCORRELATIONEXTRACTION_H_ diff --git a/PWGHF/HFC/Macros/DhCorrelationFitter.cxx b/PWGHF/HFC/Macros/DhCorrelationFitter.cxx index 8882c567535..683ae0cc04a 100644 --- a/PWGHF/HFC/Macros/DhCorrelationFitter.cxx +++ b/PWGHF/HFC/Macros/DhCorrelationFitter.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file DhCorrelationFitter.cxx -/// \brief Class to perform the azimuthal correlation fit +/// \brief class for for performing the fit of azimuthal correlations /// \author Samuele Cattaruzzi /// \author Swapnesh Santosh Khade @@ -37,8 +37,6 @@ #include #include -using namespace std; - DhCorrelationFitter::DhCorrelationFitter() : // default constructor fIsReflected(kFALSE), fTypeOfFitFunc(kConstwoGaus), @@ -53,6 +51,9 @@ DhCorrelationFitter::DhCorrelationFitter() : // default constructor fExtParsLowBounds(0x0), fExtParsUppBounds(0x0), fUseExternalPars(kFALSE), + fShiftBaselineUp(kFALSE), + fShiftBaselineDown(kFALSE), + fIsTotal(kTRUE), fNbasleinePoints(0), fBinsBaseline(0x0), fHist(0x0), @@ -64,6 +65,8 @@ DhCorrelationFitter::DhCorrelationFitter() : // default constructor fGausNS(0x0), fGausAS(0x0), fPed(0x0), + fv2AssocPart(0.), + fv2Dmeson(0.), fNSyieldBinCount(0.), fErrNSyieldBinCount(0.), fASyieldBinCount(0.), @@ -85,6 +88,9 @@ DhCorrelationFitter::DhCorrelationFitter(TH1F* histoToFit, Double_t min, Double_ fExtParsLowBounds(0x0), fExtParsUppBounds(0x0), fUseExternalPars(kFALSE), + fShiftBaselineUp(kFALSE), + fShiftBaselineDown(kFALSE), + fIsTotal(kTRUE), fNbasleinePoints(0), fBinsBaseline(0x0), fHist(0x0), @@ -96,6 +102,9 @@ DhCorrelationFitter::DhCorrelationFitter(TH1F* histoToFit, Double_t min, Double_ fGausNS(0x0), fGausAS(0x0), fPed(0x0), + fBaseTransvReg(0x0), + fv2AssocPart(0.), + fv2Dmeson(0.), fNSyieldBinCount(0.), fErrNSyieldBinCount(0.), fASyieldBinCount(0.), @@ -120,6 +129,9 @@ DhCorrelationFitter::DhCorrelationFitter(const DhCorrelationFitter& source) : // fExtParsLowBounds(source.fExtParsLowBounds), fExtParsUppBounds(source.fExtParsUppBounds), fUseExternalPars(source.fUseExternalPars), + fShiftBaselineUp(source.fShiftBaselineUp), + fShiftBaselineDown(source.fShiftBaselineDown), + fIsTotal(source.fIsTotal), fNbasleinePoints(source.fNbasleinePoints), fBinsBaseline(source.fBinsBaseline), fHist(source.fHist), @@ -131,6 +143,9 @@ DhCorrelationFitter::DhCorrelationFitter(const DhCorrelationFitter& source) : // fGausNS(source.fGausNS), fGausAS(source.fGausAS), fPed(source.fPed), + fBaseTransvReg(source.fBaseTransvReg), + fv2AssocPart(source.fv2AssocPart), + fv2Dmeson(source.fv2Dmeson), fNSyieldBinCount(source.fNSyieldBinCount), fErrNSyieldBinCount(source.fErrNSyieldBinCount), fASyieldBinCount(source.fASyieldBinCount), @@ -225,6 +240,7 @@ void DhCorrelationFitter::Fitting(Bool_t drawSplitTerm, Bool_t useExternalPars) // < 0 : fix the baseline to the weighted average of the abs(fFixBaseline) lower points // = 2 : zyam at pi/2. Fix the baseline averaging the 2 points around +-pi/2 value // = 3 : fix the baseline to the weighted average of the points passed through the function SetPointsForBaseline() + // = 4 : fix the baseline to the weighted average of the points in the transverse region in default configuration (i.e. for 32 bins the first 2, the last 2 and the 4 middle ones) // // -> fFixMean = 0 : NS & AS mean free // = 1 : NS mean fixed to 0, AS mean free @@ -238,7 +254,15 @@ void DhCorrelationFitter::Fitting(Bool_t drawSplitTerm, Bool_t useExternalPars) Printf("[INFO] DhCorrelationFitter::Fitting, Finding baseline"); FindBaseline(); } + if (fFixBase == 0) { + // set initial value of the fBaseline + fBaseline = CalculateBaseline(fHist, fIsTotal); + } Printf("[INFO] DhCorrelationFitter::Fitting, Setting Function"); + if (fTypeOfFitFunc == 7) { // case for v2 modulation + FitBaselineWv2(); // to contrain the B parameter in the fit function for the pedestal + Printf("[INFO] B parameter for v2 fit: %.3f", fBaseline); + } SetFitFunction(); if (fFixBase != 0) { @@ -248,9 +272,12 @@ void DhCorrelationFitter::Fitting(Bool_t drawSplitTerm, Bool_t useExternalPars) fFit->FixParameter(2, 0.); } if (fFixMean == 2 || fFixMean == 3) { - if (fTypeOfFitFunc != 0) + if (fTypeOfFitFunc != 0 && fTypeOfFitFunc != 3) fFit->FixParameter(5, TMath::Pi()); + if (fTypeOfFitFunc == 3 || fTypeOfFitFunc == 6) + fFit->FixParameter(2, TMath::Pi()); } + Printf("[INFO] DhCorrelationFitter::Fitting, Fitting"); TVirtualFitter::SetMaxIterations(20000); TFitResultPtr fitptr = fHist->Fit(fFit, "RIMES", "", fMinCorr, fMaxCorr); @@ -270,19 +297,39 @@ void DhCorrelationFitter::Fitting(Bool_t drawSplitTerm, Bool_t useExternalPars) fHist->SetTitle(";#Delta#varphi (rad); #frac{1}{N_{D}}#frac{dN^{assoc}}{d#Delta#varphi} (rad^{-1})"); Printf("[INFO] DhCorrelationFitter::Fitting, Now drawing, if requested"); SetSingleTermsForDrawing(drawSplitTerm); + + // NS yield from bin counting + double fNSyield = 0.; + double fNSyieldErr = 0.; + double baselinBinCount = 0; + for (int iBin = 1; iBin <= 6; iBin++) { // first six bins + fNSyield += 2 * (fHist->GetBinContent(iBin) - fBaseline) * fHist->GetBinWidth(iBin); // x2 due to the fatct the histogram is reflected + fNSyieldErr += 4 * (TMath::Power(fHist->GetBinError(iBin), 2) + TMath::Power(fErrBaseline, 2)) * (fHist->GetBinWidth(iBin) * fHist->GetBinWidth(iBin)); + baselinBinCount += fBaseline * fHist->GetBinWidth(iBin); + } + fNSyieldErr = TMath::Sqrt(fNSyieldErr); + + // AS yield from bin counting + double fASyield = 0.; + double fASyieldErr = 0.; + for (int iBin = 11; iBin <= 16; iBin++) { // last six bins + fASyield += 2 * (fHist->GetBinContent(iBin) - fBaseline) * fHist->GetBinWidth(iBin); + fASyieldErr += 4 * (TMath::Power(fHist->GetBinError(iBin), 2) + TMath::Power(fErrBaseline, 2)) * (fHist->GetBinWidth(iBin) * fHist->GetBinWidth(iBin)); + } + fASyieldErr = TMath::Sqrt(fASyieldErr); + + printf("[RESULT MINE] Bin counting results: NS Yield = %.3f +- %.3f \n[RESULT MINE] Bin counting results: AS Yield: %.3f +- %.3f \n[RESULT MINE] baseline = %.3f \n", fNSyield, fNSyieldErr, fASyield, fASyieldErr, baselinBinCount); } void DhCorrelationFitter::SetFitFunction() { - // -> fitFunc = 1: const+ G NS + G AS (w/o periodicity) - // = 2: const+ G NS + G AS (w/ periodicity) - // = 3: const+ yieldNS*[fact*(G NS)+(1- fact)*(G2 NS)] + yieldAS*(G AS) (w/ periodicity) - // = 4: const +yieldNS*(G NS) + yieldAS*[fact*(G AS)+(1- fact)*(G2 AS)] (w/ periodicity) - // = 5: v2 modulation (no gaussian terms) - // = 6: v2 modulation + G NS + G AS (w/ periodicity) - // = 7: const+ GenG NS + G AS (w/ periodicity) - // = 8: const+ GenG fixBeta NS + G AS (w/ periodicity) - // = 9: const+ GenG constrBeta NS + G AS (w/ periodicity) + // -> fitFunc = 1: const + G NS + G AS (w/o periodicity) + // = 2: const + G NS + G AS (w/ periodicity) + // = 3: const + G AS + // = 4: const + GenG NS + G AS (w/ periodicity) + // = 5: const + VonMises NS + VonMises AS (w/periodicity) + // = 6: const + VonMises AS + // = 7: baseline w v2 modulation + G NS + G AS (w/ periodicity) if (fFit) { delete fFit; @@ -372,6 +419,180 @@ void DhCorrelationFitter::SetFitFunction() fFit->SetParName(6, "AS #sigma"); break; + + case 3: + fFit = new TF1("OneGausPeriodicity", "[0]+[1]/TMath::Sqrt(2.*TMath::Pi())/[3]*TMath::Exp(-(x-[2])*(x-[2])/2./([3]*[3]))+[1]/TMath::Sqrt(2.*TMath::Pi())/[3]*TMath::Exp(-(x+2.*TMath::Pi()-[2])*(x+2.*TMath::Pi()-[2])/2./([3]*[3]))+[1]/TMath::Sqrt(2.*TMath::Pi())/[3]*TMath::Exp(-(x-2.*TMath::Pi()-[2])*(x-2.*TMath::Pi()-[2])/2./([3]*[3]))", fMinCorr, fMaxCorr); + fGausAS = new TF1("fGausASper", "[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-2.*TMath::Pi()-[1])*(x-2.*TMath::Pi()-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x+2.*TMath::Pi()-[1])*(x+2.*TMath::Pi()-[1])/2./([2]*[2]))", fMinCorr, fMaxCorr); + fPed = new TF1("fPed", "[0]", fMinCorr, fMaxCorr); + + // TODO: add possibility to use external parameters + // if(!fUseExternalPars) { + fFit->SetParLimits(0, 0, 9999.); + fFit->SetParLimits(1, 0, 999.); + fFit->SetParLimits(2, 2., 4.); + fFit->SetParLimits(3, 0, 3.14 / 2.); + + fFit->SetParameter(0, 3); + fFit->SetParameter(1, 2); + fFit->SetParameter(2, 3.14); + fFit->SetParameter(3, 0.3); + + /*} else { + for(int i=0; i SetParameter(i, fExtParsVals[i]); + fFit -> SetParLimits(i, fExtParsLowBounds[i], fExtParsUppBounds[i]); + } + } */ + + fFit->SetParName(0, "ped"); + fFit->SetParName(1, "AS Y"); + fFit->SetParName(2, "AS mean"); + fFit->SetParName(3, "AS #sigma"); + + break; + + case 4: + fFit = new TF1("kModifNSGausPeriodicity", "[0]+[1]*([7]*TMath::Sqrt(TMath::Gamma(3./[7]))/(2.*[3]*TMath::Power(TMath::Gamma(1./[7]),3./2.))*TMath::Exp(-TMath::Power(TMath::Abs(x-[2])*TMath::Sqrt(TMath::Gamma(3./[7]))/([3]*TMath::Sqrt(TMath::Gamma(1./[7]))),[7])))+[4]/TMath::Sqrt(2.*TMath::Pi())/[6]*TMath::Exp(-(x-[5])*(x-[5])/2./([6]*[6]))+[1]*([7]*TMath::Sqrt(TMath::Gamma(3./[7]))/(2.*[3]*TMath::Power(TMath::Gamma(1./[7]),3./2.))*TMath::Exp(-TMath::Power(TMath::Abs(x-2*TMath::Pi()-[2])*TMath::Sqrt(TMath::Gamma(3./[7]))/([3]*TMath::Sqrt(TMath::Gamma(1./[7]))),[7])))+[1]*([7]*TMath::Sqrt(TMath::Gamma(3./[7]))/(2.*[3]*TMath::Power(TMath::Gamma(1./[7]),3./2.))*TMath::Exp(-TMath::Power(TMath::Abs(x+2*TMath::Pi()-[2])*TMath::Sqrt(TMath::Gamma(3./[7]))/([3]*TMath::Sqrt(TMath::Gamma(1./[7]))),[7])))+[4]/TMath::Sqrt(2.*TMath::Pi())/[6]*TMath::Exp(-(x-2.*TMath::Pi()-[5])*(x-2.*TMath::Pi()-[5])/2./([6]*[6]))+[4]/TMath::Sqrt(2.*TMath::Pi())/[6]*TMath::Exp(-(x+2.*TMath::Pi()-[5])*(x+2.*TMath::Pi()-[5])/2./([6]*[6]))", fMinCorr, fMaxCorr); + fGausNS = new TF1("fModGausNSper", "[0]*([3]*TMath::Sqrt(TMath::Gamma(3./[3]))/(2.*[2]*TMath::Power(TMath::Gamma(1./[3]),3./2.))*TMath::Exp(-TMath::Power(TMath::Abs(x-[1])*TMath::Sqrt(TMath::Gamma(3./[3]))/([2]*TMath::Sqrt(TMath::Gamma(1./[3]))),[3])))+[0]*([3]*TMath::Sqrt(TMath::Gamma(3./[3]))/(2.*[2]*TMath::Power(TMath::Gamma(1./[3]),3./2.))*TMath::Exp(-TMath::Power(TMath::Abs(x-2*TMath::Pi()-[1])*TMath::Sqrt(TMath::Gamma(3./[3]))/([2]*TMath::Sqrt(TMath::Gamma(1./[3]))),[3])))+[0]*([3]*TMath::Sqrt(TMath::Gamma(3./[3]))/(2.*[2]*TMath::Power(TMath::Gamma(1./[3]),3./2.))*TMath::Exp(-TMath::Power(TMath::Abs(x+2*TMath::Pi()-[1])*TMath::Sqrt(TMath::Gamma(3./[3]))/([2]*TMath::Sqrt(TMath::Gamma(1./[3]))),[3])))", fMinCorr, fMaxCorr); + fGausAS = new TF1("fGausASper", "[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-2.*TMath::Pi()-[1])*(x-2.*TMath::Pi()-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x+2.*TMath::Pi()-[1])*(x+2.*TMath::Pi()-[1])/2./([2]*[2]))", fMinCorr, fMaxCorr); + fPed = new TF1("fPed", "[0]", fMinCorr, fMaxCorr); + + fFit->SetParLimits(0, 0., 999.); + fFit->SetParLimits(1, 0.005, 25.); + fFit->SetParLimits(2, -0.55, 0.55); + fFit->SetParLimits(3, 0, 0.8); + fFit->SetParLimits(4, 0.005, 25.); + fFit->SetParLimits(5, 2.85, 3.55); + fFit->SetParLimits(6, 0.05, 3.14 / 2.); + fFit->SetParLimits(7, 0.5, 3.5); + + // default starting pars + fFit->SetParameter(0, 1.); + fFit->SetParameter(1, 1.); + fFit->SetParameter(2, 0.); + fFit->SetParameter(3, 0.3); + fFit->SetParameter(4, 0.25); + fFit->SetParameter(5, TMath::Pi()); + fFit->SetParameter(6, 0.3); + fFit->SetParameter(7, 2); + + fFit->SetParName(0, "ped"); + fFit->SetParName(1, "NS Y"); + fFit->SetParName(2, "NS mean"); + fFit->SetParName(3, "NS #sigma"); + fFit->SetParName(4, "AS Y"); + fFit->SetParName(5, "AS mean"); + fFit->SetParName(6, "AS #sigma"); + fFit->SetParName(7, "NS shape par"); // beta of the gen. gaussian + break; + + case 5: + fFit = new TF1("kVonMises", "[0] +[1]/(2*TMath::Pi()*TMath::BesselI0([3]))*TMath::Exp([3]*TMath::Cos(x- 2*TMath::Pi() - [2])) + [4]/(2*TMath::Pi()*TMath::BesselI0([6]))*TMath::Exp([6]*TMath::Cos(x- 2*TMath::Pi()-[5]))", fMinCorr, fMaxCorr); + fGausNS = new TF1("fVonMisesNS", "[0]/(2*TMath::Pi()*TMath::BesselI0([2]))*TMath::Exp([2]*TMath::Cos(x- 2*TMath::Pi() - [1]))", fMinCorr, fMaxCorr); + fGausAS = new TF1("fVonMisesAS", "[0]/(2*TMath::Pi()*TMath::BesselI0([2]))*TMath::Exp([2]*TMath::Cos(x- 2*TMath::Pi()-[1]))", fMinCorr, fMaxCorr); + fPed = new TF1("fPed", "[0]", fMinCorr, fMaxCorr); + + fFit->SetParLimits(0, 0., 999.); + fFit->SetParLimits(1, 0.005, 25.); + fFit->SetParLimits(2, -0.55, 0.55); + fFit->SetParLimits(3, 0, 15.); + fFit->SetParLimits(4, 0.005, 25.); + fFit->SetParLimits(5, 2.85, 3.55); + fFit->SetParLimits(6, 0., 15.); + + // default starting pars + fFit->SetParameter(0, 3.5); + fFit->SetParameter(1, 0.8); + fFit->SetParameter(2, 0.); + fFit->SetParameter(3, 1.); + fFit->SetParameter(4, 1.5); + fFit->SetParameter(5, TMath::Pi()); + fFit->SetParameter(6, 1.); + + fFit->SetParName(0, "ped"); + fFit->SetParName(1, "NS Y"); + fFit->SetParName(2, "NS mean"); + fFit->SetParName(3, "NS #sigma"); + fFit->SetParName(4, "AS Y"); + fFit->SetParName(5, "AS mean"); + fFit->SetParName(6, "AS #sigma"); + + break; + + case 6: + fFit = new TF1("kSingleVonMises", "[0] +[1]/(2*TMath::Pi()*TMath::BesselI0([3]))*TMath::Exp([3]*TMath::Cos(x- 2*TMath::Pi() - [2]))", fMinCorr, fMaxCorr); + fGausAS = new TF1("fVonMisesAS", "[0]/(2*TMath::Pi()*TMath::BesselI0([2]))*TMath::Exp([2]*TMath::Cos(x- 2*TMath::Pi()-[1]))", fMinCorr, fMaxCorr); + fPed = new TF1("fPed", "[0]", fMinCorr, fMaxCorr); + + fFit->SetParLimits(0, 0., 999.); + fFit->SetParLimits(1, 0.005, 25.); + fFit->SetParLimits(2, 2.85, 3.55); + fFit->SetParLimits(3, 0., 15.); + + // default starting pars + fFit->SetParameter(0, 3.5); + fFit->SetParameter(1, 1.5); + fFit->SetParameter(2, TMath::Pi()); + fFit->SetParameter(3, 1.); + + fFit->SetParName(0, "ped"); + fFit->SetParName(1, "AS Y"); + fFit->SetParName(2, "AS mean"); + fFit->SetParName(3, "AS #sigma"); + + break; + + case 7: // case 2 Gaus w periodicity + v2 modulation + + fFit = new TF1("kTwoGausPeriodicityPlusV2modulation", "[1]/TMath::Sqrt(2.*TMath::Pi())/[3]*TMath::Exp(-(x-[2])*(x-[2])/2./([3]*[3]))+[4]/TMath::Sqrt(2.*TMath::Pi())/[6]*TMath::Exp(-(x-[5])*(x-[5])/2./([6]*[6]))+[1]/TMath::Sqrt(2.*TMath::Pi())/[3]*TMath::Exp(-(x-2.*TMath::Pi()-[2])*(x-2.*TMath::Pi()-[2])/2./([3]*[3]))+[1]/TMath::Sqrt(2.*TMath::Pi())/[3]*TMath::Exp(-(x+2.*TMath::Pi()-[2])*(x+2.*TMath::Pi()-[2])/2./([3]*[3]))+[4]/TMath::Sqrt(2.*TMath::Pi())/[6]*TMath::Exp(-(x+2.*TMath::Pi()-[5])*(x+2.*TMath::Pi()-[5])/2./([6]*[6]))+[4]/TMath::Sqrt(2.*TMath::Pi())/[6]*TMath::Exp(-(x-2.*TMath::Pi()-[5])*(x-2.*TMath::Pi()-[5])/2./([6]*[6]))+[0]*(1+2*[7]*[8]*TMath::Cos(2*x))", fMinCorr, fMaxCorr); + fGausNS = new TF1("fGausNSper", "[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-2.*TMath::Pi()-[1])*(x-2.*TMath::Pi()-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x+2.*TMath::Pi()-[1])*(x+2.*TMath::Pi()-[1])/2./([2]*[2]))", fMinCorr, fMaxCorr); + fGausAS = new TF1("fGausASper", "[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-2.*TMath::Pi()-[1])*(x-2.*TMath::Pi()-[1])/2./([2]*[2]))+[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x+2.*TMath::Pi()-[1])*(x+2.*TMath::Pi()-[1])/2./([2]*[2]))", fMinCorr, fMaxCorr); + fPed = new TF1("fPedv2Mod", "[0]*(1+2*[1]*[2]*TMath::Cos(2*x))", fMinCorr, fMaxCorr); + + fFit->SetParLimits(0, 0., 999.); + fFit->SetParLimits(1, 0, 999.); + fFit->SetParLimits(2, -0.55, 0.55); + fFit->SetParLimits(3, 0, 3.14 / 3.); + fFit->SetParLimits(4, 0, 999.); + fFit->SetParLimits(5, 2.85, 3.55); + fFit->SetParLimits(6, 0, 3.14 / 2.); + fFit->SetParLimits(7, -1, 1); + fFit->SetParLimits(8, -1, 1); + + fFit->FixParameter(0, fBaseline); + fFit->SetParameter(1, 3); + fFit->SetParameter(2, 0.); + fFit->SetParameter(3, 0.3); + fFit->SetParameter(4, 2); + fFit->SetParameter(5, TMath::Pi()); + fFit->SetParameter(6, 0.3); + fFit->SetParameter(7, 0); + fFit->SetParameter(8, 0); + + if (fUseExternalPars) { // overwrites previous configuration + for (int i = 0; i < fNpars; i++) { + fFit->SetParameter(i, fExtParsVals[i]); + fFit->SetParLimits(i, fExtParsLowBounds[i], fExtParsUppBounds[i]); + } + } + + fFit->FixParameter(7, fv2AssocPart); + fFit->FixParameter(8, fv2Dmeson); + + fPed->FixParameter(0, fBaseline); + fPed->FixParameter(1, fv2AssocPart); + fPed->FixParameter(2, fv2Dmeson); + + fFit->SetParName(0, "ped"); + fFit->SetParName(1, "NS Y"); + fFit->SetParName(2, "NS mean"); + fFit->SetParName(3, "NS #sigma"); + fFit->SetParName(4, "AS Y"); + fFit->SetParName(5, "AS mean"); + fFit->SetParName(6, "AS #sigma"); + fFit->SetParName(7, "v_{2} hadron"); + fFit->SetParName(8, "v_{2} D meson"); + break; } } @@ -412,6 +633,16 @@ Double_t DhCorrelationFitter::FindBaseline() fBaseline = min; fErrBaseline = fHist->GetBinError(iBin); + if (fShiftBaselineUp) { + fBaseline += fErrBaseline; + printf("[INFO] Shift baseline up of its statistical uncertainty"); + } + + if (fShiftBaselineDown) { + fBaseline -= fErrBaseline; + printf("[INFO] Shift baseline down of its statistical uncertainty"); + } + return fBaseline; } @@ -441,6 +672,17 @@ Double_t DhCorrelationFitter::FindBaseline() printf("[RESULT] Average fBaseline: %.3f +- %.3f", Av, errAv); fBaseline = Av; fErrBaseline = errAv; + + if (fShiftBaselineUp) { + fBaseline += fErrBaseline; + printf("[INFO] Shift baseline up of its statistical uncertainty"); + } + + if (fShiftBaselineDown) { + fBaseline -= fErrBaseline; + printf("[INFO] Shift baseline down of its statistical uncertainty"); + } + return fBaseline; } @@ -465,6 +707,16 @@ Double_t DhCorrelationFitter::FindBaseline() fBaseline = Av; fErrBaseline = errAv; + if (fShiftBaselineUp) { + fBaseline += fErrBaseline; + printf("[INFO] Shift baseline up of its statistical uncertainty"); + } + + if (fShiftBaselineDown) { + fBaseline -= fErrBaseline; + printf("[INFO] Shift baseline down of its statistical uncertainty"); + } + return fBaseline; } @@ -485,6 +737,33 @@ Double_t DhCorrelationFitter::FindBaseline() fBaseline = Av; fErrBaseline = errAv; + if (fShiftBaselineUp) { + fBaseline += fErrBaseline; + printf("[INFO] Shift baseline up of its statistical uncertainty \n"); + } + + if (fShiftBaselineDown) { + fBaseline -= fErrBaseline; + printf("[INFO] Shift baseline down of its statistical uncertainty \n"); + } + + return fBaseline; + } + + if (fFixBase == 4) { + fBaseline = CalculateBaseline(fHist, fIsTotal); // TODO: add the option for total range/ reflected range to pass in input + fErrBaseline = CalculateBaselineError(fHist, fIsTotal); + + if (fShiftBaselineUp) { + fBaseline += fErrBaseline; + printf("[INFO] Shift baseline up of its statistical uncertainty \n"); + } + + if (fShiftBaselineDown) { + fBaseline -= fErrBaseline; + printf("[INFO] Shift baseline down of its statistical uncertainty \n"); + } + return fBaseline; } @@ -492,6 +771,26 @@ Double_t DhCorrelationFitter::FindBaseline() return -1.; } +void DhCorrelationFitter::FitBaselineWv2() +{ + + fBaseTransvReg = new TF1("fBaseTransvReg", [](double* x, double* p) { + double xx = x[0]; // x value + if ((xx >= -TMath::Pi()/2 && xx <= -3*TMath::Pi()/8) || (xx >= 3*TMath::Pi()/8 && xx <= 5*TMath::Pi()/8) || (xx >= 11*TMath::Pi()/8 && xx <= 3*TMath::Pi()/2)) { + // Gaussian example: p[0] = amplitude, p[1] = mean, p[2] = sigma + return p[0]*(1+2*p[1]*p[2]*TMath::Cos(2*xx)); + } + return 0.; }, -TMath::Pi() / 2, 3 * TMath::Pi() / 2, 3); // Function valid for [0,10], with 3 parameters + + fBaseTransvReg->FixParameter(1, fv2AssocPart); + fBaseTransvReg->FixParameter(2, fv2Dmeson); + + TFitResultPtr rFit = fHist->Fit(fBaseTransvReg, "RIMES", "", fMinCorr, fMaxCorr); + fBaseline = fBaseTransvReg->GetParameter(0); + + return; +} + void DhCorrelationFitter::CalculateYieldsAboveBaseline() { @@ -502,23 +801,28 @@ void DhCorrelationFitter::CalculateYieldsAboveBaseline() cout << "[RESULT] Baseline: " << fBaseline << " +- " << fErrBaseline << endl; Int_t binMinNS = fHist->FindBin(-1.5); // slightly more than -pi/2 if (binMinNS < 1) - binMinNS = 1; // with this, it is ok even in the case of a reflected fHist (range 0 - pi) - Int_t binMaxNS = fHist->FindBin(1.5); // slightly less than +pi/2 - Int_t binMinAS = fHist->FindBin(1.6); // slightly more than +pi/2 - Int_t binMaxAS = fHist->FindBin(3.14 + 1.5); // slightly less than +3pi/2 + binMinNS = 1; // with this, it is ok even in the case of a reflected fHist (range 0 - pi) + Int_t binMaxNS = 6; // fHist -> FindBin(1.5); // slightly less than +pi/2 + Int_t binMinAS = 11; // fHist -> FindBin(1.6); // slightly more than +pi/2 + Int_t binMaxAS = 16; // fHist -> FindBin(3.14+1.5); // slightly less than +3pi/2 if (binMaxAS > fHist->GetNbinsX()) - binMaxNS = fHist->GetNbinsX(); // with this, it is ok even in the case of a reflected fHist (range 0 - pi) TODO: maybe binMaxAS + binMaxAS = fHist->GetNbinsX(); // with this, it is ok even in the case of a reflected fHist (range 0 - pi) + cout << "N bins : " << fHist->GetNbinsX() << endl; + cout << "binMinNS : " << binMinNS << endl; + cout << "binMaxNS : " << binMaxNS << endl; + cout << "binMinAS : " << binMinAS << endl; + cout << "binMaxAS : " << binMaxAS << endl; // Near Side Yield from bin counting for (Int_t bmNS = binMinNS; bmNS <= binMaxNS; bmNS++) { - fNSyieldBinCount += (fHist->GetBinContent(bmNS) - fBaseline) * fHist->GetBinWidth(bmNS); - fErrNSyieldBinCount += (fHist->GetBinError(bmNS) * fHist->GetBinError(bmNS)) * fHist->GetBinWidth(bmNS) * fHist->GetBinWidth(bmNS); + fNSyieldBinCount += 2 * (fHist->GetBinContent(bmNS) - fBaseline) * fHist->GetBinWidth(bmNS); + fErrNSyieldBinCount += 4 * (fHist->GetBinError(bmNS) * fHist->GetBinError(bmNS)) * fHist->GetBinWidth(bmNS) * fHist->GetBinWidth(bmNS); } fErrNSyieldBinCount = TMath::Sqrt(fErrNSyieldBinCount); // Away Side Yield from bin counting for (Int_t bmAS = binMinAS; bmAS <= binMaxAS; bmAS++) { - fASyieldBinCount += (fHist->GetBinContent(bmAS) - fBaseline) * fHist->GetBinWidth(bmAS); - fErrASyieldBinCount += (fHist->GetBinError(bmAS) * fHist->GetBinError(bmAS)) * fHist->GetBinWidth(bmAS) * fHist->GetBinWidth(bmAS); + fASyieldBinCount += 2 * (fHist->GetBinContent(bmAS) - fBaseline) * fHist->GetBinWidth(bmAS); + fErrASyieldBinCount += 4 * (fHist->GetBinError(bmAS) * fHist->GetBinError(bmAS)) * fHist->GetBinWidth(bmAS) * fHist->GetBinWidth(bmAS); } fErrASyieldBinCount = TMath::Sqrt(fErrASyieldBinCount); @@ -527,55 +831,322 @@ void DhCorrelationFitter::CalculateYieldsAboveBaseline() return; } +Double_t DhCorrelationFitter::CalculateBaseline(TH1F*& histo, Bool_t totalRange) +{ + + // total range = 2*Pi + // half range = Pi , for histogram reflected under symmetric assumption + + Double_t baseline, errBaseline; + Int_t nBinsPhi = histo->GetNbinsX(); + Int_t binPhiHalf = nBinsPhi / 2; + Int_t binPhiHalfMinus1 = nBinsPhi / 2 - 1; + Int_t binPhiHalfPlus1 = nBinsPhi / 2 + 1; + Int_t binPhiHalfPlus2 = nBinsPhi / 2 + 1; + + if (totalRange) { + printf("[INFO] Using total deltaPhi range \n"); + // baseline evaluated considering: the two first points, the last two points and four points in the middle (corresponding to the outer points) + if (nBinsPhi >= 32) { + printf("[INFO] nBinsPhi >= 32 \n"); + baseline = + ((histo->GetBinContent(1)) * (1. / TMath::Power(histo->GetBinError(1), 2)) + + (histo->GetBinContent(2)) * (1. / TMath::Power(histo->GetBinError(2), 2)) + + (histo->GetBinContent(binPhiHalfMinus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (histo->GetBinContent(binPhiHalfPlus2)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2)) + + (histo->GetBinContent(nBinsPhi - 1)) * (1. / TMath::Power(histo->GetBinError(nBinsPhi - 1), 2)) + + (histo->GetBinContent(nBinsPhi)) * (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))) / + ((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(2), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi - 1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } else { + printf("[INFO] nBinsPhi < 32 \n"); + baseline = + ((histo->GetBinContent(1)) * (1. / TMath::Power(histo->GetBinError(1), 2)) + + (histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (histo->GetBinContent(nBinsPhi)) * (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))) / + ((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } + } else { + printf("[INFO] Using reflected deltaPhi range \n"); + // baseline evaluated using the 4 middle points in the transverese region + if (nBinsPhi >= 16) { + printf("[INFO] 4 central points in the transverse region for baseline \n"); + baseline = + ((histo->GetBinContent(binPhiHalfMinus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (histo->GetBinContent(binPhiHalfPlus2)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))) / + ((1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))); + } else { + baseline = + ((histo->GetBinContent(binPhiHalf)) * (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (histo->GetBinContent(binPhiHalfPlus1)) * (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2))) / + ((1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2))); + } + } + + return baseline; +} + +Double_t DhCorrelationFitter::CalculateBaselineError(TH1F*& histo, Bool_t totalRange) +{ + + // total range = 2*Pi + // half range = Pi , for histogram reflected under symmetric assumption + + Double_t errBaseline; + Int_t nBinsPhi = histo->GetNbinsX(); + Int_t binPhiHalf = nBinsPhi / 2; + Int_t binPhiHalfMinus1 = nBinsPhi / 2 - 1; + Int_t binPhiHalfPlus1 = nBinsPhi / 2 + 1; + Int_t binPhiHalfPlus2 = nBinsPhi / 2 + 1; + + if (totalRange) { + // baseline evaluated considering: the two first points, the last two points and four points in the middle (corresponding to the outer points) + if (nBinsPhi >= 32) { + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(2), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi - 1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } else { // fon nBinsPhi = 16 (rebin 4) + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(nBinsPhi), 2))); + } + } else { + // baseline evaluated using the 4 middle points in the transverese region + if (nBinsPhi >= 32) { + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(binPhiHalfMinus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus2), 2))); + } else { + errBaseline = 1. / + TMath::Sqrt((1. / TMath::Power(histo->GetBinError(binPhiHalf), 2)) + + (1. / TMath::Power(histo->GetBinError(binPhiHalfPlus1), 2))); + } + } + + return errBaseline; +} + void DhCorrelationFitter::SetSingleTermsForDrawing(Bool_t draw) { Double_t* par = 0; - if (fTypeOfFitFunc == 1 || fTypeOfFitFunc == 2) { + if (fTypeOfFitFunc == 1 || fTypeOfFitFunc == 2 || fTypeOfFitFunc == 5) { par = new Double_t[7]; + } else if (fTypeOfFitFunc == 3 || fTypeOfFitFunc == 6) { + par = new Double_t[4]; + } else if (fTypeOfFitFunc == 4) { + par = new Double_t[8]; + } else if (fTypeOfFitFunc == 7) { + par = new Double_t[9]; } else { Printf("[ERROR] DhCorrelationFitter::SetSingleTermsForDrawing, wrong type of function"); return; } - fFit->GetParameters(par); - fFit->SetLineWidth(4); - - fPed->SetParameter(0, par[0]); - fPed->SetLineColor(6); // pink - fPed->SetLineStyle(9); - fPed->SetLineWidth(4); - - fGausNS->SetParameter(0, par[1]); - fGausNS->SetParameter(1, par[2]); - fGausNS->SetParameter(2, par[3]); - fGausAS->SetParameter(0, par[4]); - fGausAS->SetParameter(1, par[5]); - fGausAS->SetParameter(2, par[6]); - - fGausNS->SetLineStyle(9); - fGausNS->SetLineColor(kBlue); - fGausAS->SetLineStyle(9); - fGausAS->SetLineColor(kGreen); - fGausNS->SetLineWidth(4); - fGausAS->SetLineWidth(4); - - TPaveText* pvStatTests1 = new TPaveText(0.51, 0.58, 0.85, 0.90, "NDC"); - pvStatTests1->SetFillStyle(0); - pvStatTests1->SetTextSize(0.035); - pvStatTests1->SetBorderSize(0); - TText *t0, *t1, *t2, *t3, *t4, *t5, *t6; - t0 = pvStatTests1->AddText(0., 1.00, Form("#chi^{2}/ndf = %.1f/%d ", fFit->GetChisquare(), fFit->GetNDF())); - t1 = pvStatTests1->AddText(0., 0.80, Form("Ped = %.3f#pm%.3f ", fBaseline, fErrBaseline)); - t2 = pvStatTests1->AddText(0., 0.65, Form("NS Y = %.3f#pm%.3f ", fFit->GetParameter("NS Y"), fFit->GetParError(fFit->GetParNumber("NS Y")))); - t3 = pvStatTests1->AddText(0., 0.50, Form("NS #sigma = %.3f#pm%.3f ", fFit->GetParameter("NS #sigma"), fFit->GetParError(fFit->GetParNumber("NS #sigma")))); - t4 = pvStatTests1->AddText(0., 0.35, Form("AS Y = %.3f#pm%.3f ", fFit->GetParameter("AS Y"), fFit->GetParError(fFit->GetParNumber("AS Y")))); - t5 = pvStatTests1->AddText(0., 0.20, Form("AS #sigma = %.3f#pm%.3f ", fFit->GetParameter("AS #sigma"), fFit->GetParError(fFit->GetParNumber("AS #sigma")))); - - if (draw) { - fFit->Draw("same"); - fPed->Draw("same"); - fGausAS->Draw("same"); - fGausNS->Draw("same"); - pvStatTests1->Draw("same"); + if (fTypeOfFitFunc == 3 || fTypeOfFitFunc == 6) { + fFit->GetParameters(par); + fFit->SetLineWidth(4); + + fPed->SetParameter(0, par[0]); + fPed->SetLineColor(6); // pink + fPed->SetLineStyle(9); + fPed->SetLineWidth(4); + + fGausAS->SetParameter(0, par[1]); + fGausAS->SetParameter(1, par[2]); + fGausAS->SetParameter(2, par[3]); + + fGausAS->SetLineStyle(9); + fGausAS->SetLineColor(kGreen); + fGausAS->SetLineWidth(4); + + TPaveText* pvStatTests1 = new TPaveText(0.51, 0.58, 0.85, 0.90, "NDC"); + pvStatTests1->SetFillStyle(0); + pvStatTests1->SetTextSize(0.045); + pvStatTests1->SetBorderSize(0); + TText *t0, *t1, *t2, *t3; + t0 = pvStatTests1->AddText(0., 1.00, Form("#chi^{2}/ndf = %.1f/%d ", fFit->GetChisquare(), fFit->GetNDF())); + t1 = pvStatTests1->AddText(0., 0.80, Form("Ped = %.3f#pm%.3f ", fBaseline, fErrBaseline)); + t2 = pvStatTests1->AddText(0., 0.65, Form("AS Y = %.3f#pm%.3f ", fFit->GetParameter("AS Y"), fFit->GetParError(fFit->GetParNumber("AS Y")))); + t3 = pvStatTests1->AddText(0., 0.50, Form("AS #sigma = %.3f#pm%.3f ", fFit->GetParameter("AS #sigma"), fFit->GetParError(fFit->GetParNumber("AS #sigma")))); + + if (draw) { + fFit->Draw("same"); + fPed->Draw("same"); + fGausAS->Draw("same"); + pvStatTests1->Draw("same"); + } + } else if (fTypeOfFitFunc == 4) { + fFit->GetParameters(par); + fFit->SetLineWidth(4); + + fPed->SetParameter(0, par[0]); + fPed->SetLineColor(6); // pink + fPed->SetLineStyle(9); + fPed->SetLineWidth(4); + + fGausNS->SetParameter(0, par[1]); + fGausNS->SetParameter(1, par[2]); + fGausNS->SetParameter(2, par[3]); + fGausNS->SetParameter(3, par[7]); + fGausAS->SetParameter(0, par[4]); + fGausAS->SetParameter(1, par[5]); + fGausAS->SetParameter(2, par[6]); + + fGausNS->SetLineStyle(9); + fGausNS->SetLineColor(kBlue); + fGausAS->SetLineStyle(9); + fGausAS->SetLineColor(kGreen); + fGausNS->SetLineWidth(4); + fGausAS->SetLineWidth(4); + + TPaveText* pvStatTests1 = new TPaveText(0.51, 0.58, 0.85, 0.90, "NDC"); + pvStatTests1->SetFillStyle(0); + pvStatTests1->SetTextSize(0.045); + pvStatTests1->SetBorderSize(0); + TText *t0, *t1, *t2, *t3, *t4, *t5, *t6; + t0 = pvStatTests1->AddText(0., 1.00, Form("#chi^{2}/ndf = %.1f/%d ", fFit->GetChisquare(), fFit->GetNDF())); + t1 = pvStatTests1->AddText(0., 0.80, Form("Ped = %.3f#pm%.3f ", fBaseline, fErrBaseline)); + t2 = pvStatTests1->AddText(0., 0.65, Form("NS Y = %.3f#pm%.3f ", fFit->GetParameter("NS Y"), fFit->GetParError(fFit->GetParNumber("NS Y")))); + t3 = pvStatTests1->AddText(0., 0.50, Form("NS #sigma = %.3f#pm%.3f ", fFit->GetParameter("NS #sigma"), fFit->GetParError(fFit->GetParNumber("NS #sigma")))); + t4 = pvStatTests1->AddText(0., 0.35, Form("AS Y = %.3f#pm%.3f ", fFit->GetParameter("AS Y"), fFit->GetParError(fFit->GetParNumber("AS Y")))); + t5 = pvStatTests1->AddText(0., 0.20, Form("AS #sigma = %.3f#pm%.3f ", fFit->GetParameter("AS #sigma"), fFit->GetParError(fFit->GetParNumber("AS #sigma")))); + t6 = pvStatTests1->AddText(0., 0.05, Form("#beta = %.3f#pm%.3f ", fFit->GetParameter("NS shape par"), fFit->GetParError(fFit->GetParNumber("NS shape par")))); + + if (draw) { + fFit->Draw("same"); + fPed->Draw("same"); + fGausAS->Draw("same"); + fGausNS->Draw("same"); + pvStatTests1->Draw("same"); + } + } else if (fTypeOfFitFunc == 7) { + fFit->GetParameters(par); + fFit->SetLineWidth(4); + + fBaseTransvReg->SetLineColor(15); + fBaseTransvReg->SetLineStyle(9); + fBaseTransvReg->SetLineWidth(4); + + fPed->SetParameter(0, par[0]); + fPed->SetParameter(1, par[7]); + fPed->SetParameter(2, par[8]); + fPed->SetLineColor(6); // pink + fPed->SetLineStyle(9); + fPed->SetLineWidth(4); + + fGausNS->SetParameter(0, par[1]); + fGausNS->SetParameter(1, par[2]); + fGausNS->SetParameter(2, par[3]); + // fGausNS -> SetParameter(3, par[7]); + fGausAS->SetParameter(0, par[4]); + fGausAS->SetParameter(1, par[5]); + fGausAS->SetParameter(2, par[6]); + + fGausNS->SetLineStyle(9); + fGausNS->SetLineColor(kBlue); + fGausAS->SetLineStyle(9); + fGausAS->SetLineColor(kGreen); + fGausNS->SetLineWidth(4); + fGausAS->SetLineWidth(4); + + TPaveText* pvStatTests1 = new TPaveText(0.51, 0.58, 0.85, 0.90, "NDC"); + pvStatTests1->SetFillStyle(0); + pvStatTests1->SetTextSize(0.045); + pvStatTests1->SetBorderSize(0); + TText *t0, *t1, *t2, *t3, *t4, *t5, *t6, *t7, *t8; + t0 = pvStatTests1->AddText(0., 1.00, Form("#chi^{2}/ndf = %.1f/%d ", fFit->GetChisquare(), fFit->GetNDF())); + t2 = pvStatTests1->AddText(0., 0.80, Form("NS Y = %.3f#pm%.3f ", fFit->GetParameter("NS Y"), fFit->GetParError(fFit->GetParNumber("NS Y")))); + t3 = pvStatTests1->AddText(0., 0.65, Form("NS #sigma = %.3f#pm%.3f ", fFit->GetParameter("NS #sigma"), fFit->GetParError(fFit->GetParNumber("NS #sigma")))); + t4 = pvStatTests1->AddText(0., 0.50, Form("AS Y = %.3f#pm%.3f ", fFit->GetParameter("AS Y"), fFit->GetParError(fFit->GetParNumber("AS Y")))); + t5 = pvStatTests1->AddText(0., 0.35, Form("AS #sigma = %.3f#pm%.3f ", fFit->GetParameter("AS #sigma"), fFit->GetParError(fFit->GetParNumber("AS #sigma")))); + // t6 = pvStatTests1 -> AddText(0., 0.20, Form("#beta = %.3f#pm%.3f ", fFit -> GetParameter("NS shape par"), fFit -> GetParError(fFit->GetParNumber("NS shape par")))); + + TPaveText* pvStatTests2 = new TPaveText(0.51, 0.28, 0.85, 0.60, "NDC"); + pvStatTests2->SetFillStyle(0); + pvStatTests2->SetTextSize(0.045); + pvStatTests2->SetBorderSize(0); + t1 = pvStatTests2->AddText(0., 1.00, Form("Ped = %.3f#pm%.3f ", fFit->GetParameter("ped"), fErrBaseline /*fFit -> GetParError(fFit->GetParNumber("ped")*/)); + t7 = pvStatTests2->AddText(0., 0.65, Form("v_{2}^{hadron} = %.3f#pm%.3f ", fFit->GetParameter("v_{2} hadron"), fFit->GetParError(fFit->GetParNumber("v_{2} hadron")))); + t8 = pvStatTests2->AddText(0., 0.35, Form("v_{2}^{D} = %.3f#pm%.3f ", fFit->GetParameter("v_{2} D meson"), fFit->GetParError(fFit->GetParNumber("v_{2} D meson")))); + + if (draw) { + fFit->Draw("same"); + fPed->Draw("same"); + fBaseTransvReg->Draw("same"); + fGausAS->Draw("same"); + fGausNS->Draw("same"); + pvStatTests1->Draw("same"); + pvStatTests2->Draw("same"); + } + } else { + fFit->GetParameters(par); + fFit->SetLineWidth(4); + + fPed->SetParameter(0, par[0]); + fPed->SetLineColor(6); // pink + fPed->SetLineStyle(9); + fPed->SetLineWidth(4); + + fGausNS->SetParameter(0, par[1]); + fGausNS->SetParameter(1, par[2]); + fGausNS->SetParameter(2, par[3]); + fGausAS->SetParameter(0, par[4]); + fGausAS->SetParameter(1, par[5]); + fGausAS->SetParameter(2, par[6]); + + fGausNS->SetLineStyle(9); + fGausNS->SetLineColor(kBlue); + fGausAS->SetLineStyle(9); + fGausAS->SetLineColor(kGreen); + fGausNS->SetLineWidth(4); + fGausAS->SetLineWidth(4); + + TPaveText* pvStatTests1 = new TPaveText(0.51, 0.58, 0.85, 0.90, "NDC"); + pvStatTests1->SetFillStyle(0); + pvStatTests1->SetTextSize(0.045); + pvStatTests1->SetBorderSize(0); + TText *t0, *t1, *t2, *t3, *t4, *t5, *t6; + t0 = pvStatTests1->AddText(0., 1.00, Form("#chi^{2}/ndf = %.1f/%d ", fFit->GetChisquare(), fFit->GetNDF())); + t1 = pvStatTests1->AddText(0., 0.80, Form("Ped = %.3f#pm%.3f ", fBaseline, fErrBaseline)); + t2 = pvStatTests1->AddText(0., 0.65, Form("NS Y = %.3f#pm%.3f ", fFit->GetParameter("NS Y"), fFit->GetParError(fFit->GetParNumber("NS Y")))); + t3 = pvStatTests1->AddText(0., 0.50, Form("NS #sigma = %.3f#pm%.3f ", fFit->GetParameter("NS #sigma"), fFit->GetParError(fFit->GetParNumber("NS #sigma")))); + t4 = pvStatTests1->AddText(0., 0.35, Form("AS Y = %.3f#pm%.3f ", fFit->GetParameter("AS Y"), fFit->GetParError(fFit->GetParNumber("AS Y")))); + t5 = pvStatTests1->AddText(0., 0.20, Form("AS #sigma = %.3f#pm%.3f ", fFit->GetParameter("AS #sigma"), fFit->GetParError(fFit->GetParNumber("AS #sigma")))); + + if (draw) { + fFit->Draw("same"); + fPed->Draw("same"); + fGausAS->Draw("same"); + fGausNS->Draw("same"); + pvStatTests1->Draw("same"); + } } } diff --git a/PWGHF/HFC/Macros/DhCorrelationFitter.h b/PWGHF/HFC/Macros/DhCorrelationFitter.h index a4bdf3c3bb1..8a9db7fe6e7 100644 --- a/PWGHF/HFC/Macros/DhCorrelationFitter.h +++ b/PWGHF/HFC/Macros/DhCorrelationFitter.h @@ -17,18 +17,24 @@ #ifndef PWGHF_HFC_MACROS_DHCORRELATIONFITTER_H_ #define PWGHF_HFC_MACROS_DHCORRELATIONFITTER_H_ -#include - -#include #include #include +#include + +#include + class DhCorrelationFitter { public: enum FunctionType { kConstwoGaus = 1, - kTwoGausPeriodicity = 2 }; + kTwoGausPeriodicity = 2, + kSingleGaus = 3, + kGenGaus = 4, + kVonMises = 5, + kSingleVonMises = 6, + kTwoGausPeriodicityPlusV2modulation = 7 }; /// Constructors DhCorrelationFitter(); @@ -51,11 +57,25 @@ class DhCorrelationFitter } void SetExternalValsAndBounds(Int_t nPars, Double_t* vals, Double_t* lowBounds, Double_t* uppBounds); void SetPointsForBaseline(Int_t nBaselinePoints, Int_t* binsBaseline); + void SetReflectedCorrHisto(Bool_t isReflected) { fIsTotal = !isReflected; } + void SetBaselineUpOrDown(Bool_t baseUp, Bool_t baseDown) + { + fShiftBaselineUp = baseUp; + fShiftBaselineDown = baseDown; + } + void Setv2(Double_t v2AssocPart, Double_t v2Dmeson) + { + fv2AssocPart = v2AssocPart; + fv2Dmeson = v2Dmeson; + } /// Functions for fitting void Fitting(Bool_t drawSplitTerm = kTRUE, Bool_t useExternalPars = kFALSE); void SetFitFunction(); void CalculateYieldsAboveBaseline(); + void FitBaselineWv2(); + Double_t CalculateBaseline(TH1F*& histo, Bool_t totalRange = kTRUE); + Double_t CalculateBaselineError(TH1F*& histo, Bool_t totalRange = kTRUE); void SetSingleTermsForDrawing(Bool_t draw); Double_t FindBaseline(); @@ -92,13 +112,17 @@ class DhCorrelationFitter private: TH1F* fHist; // 1D azimuthal correlation histogram - TF1* fFit; // Total fit function - TF1* fGausNS; // Near-Side (NS) Gaussian - TF1* fGausAS; // Away-Side (AS) Gaussian - TF1* fPed; // Baseline function + TF1* fFit; // Total fit function + TF1* fGausNS; // Near-Side (NS) Gaussian + TF1* fGausAS; // Away-Side (AS) Gaussian + TF1* fPed; // Baseline function + TF1* fBaseTransvReg; // Baseline function with v2 - Bool_t fIsReflected; - Bool_t fUseExternalPars; // To use external fit parameters initial values and bounds + Bool_t fIsReflected; // To use if reflected azimuthal correlation are given as input + Bool_t fUseExternalPars; // To use external fit parameters initial values and bounds + Bool_t fShiftBaselineUp; // To shift the baseline up of its statistical uncertainty + Bool_t fShiftBaselineDown; // To shift baseline down of its statistical uncertainty + Bool_t fIsTotal; // Total range of 2*pi in the azimuthal correlation distribution FunctionType fTypeOfFitFunc; // Type of fit function @@ -121,6 +145,8 @@ class DhCorrelationFitter Double_t fErrNSyieldBinCount; // NS Yield error from bin counting Double_t fASyieldBinCount; // AS Yield from bin counting Double_t fErrASyieldBinCount; // AS Yield error from bin counting + Double_t fv2AssocPart; // v2 associated particles + Double_t fv2Dmeson; // v2 of D mesons Double_t* fExtParsVals; // Fit parameters initial values Double_t* fExtParsLowBounds; // Fit parameters lower bounds diff --git a/PWGHF/HFC/Macros/ExtractOutputCorrel.C b/PWGHF/HFC/Macros/ExtractOutputCorrel.C index 301434f66e2..e263a810ff1 100644 --- a/PWGHF/HFC/Macros/ExtractOutputCorrel.C +++ b/PWGHF/HFC/Macros/ExtractOutputCorrel.C @@ -9,19 +9,26 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file ExtractOutputCorrel.C -/// \brief Macro to perform the correlation extraction +/// \file DhCorrelationExtraction.cxx +/// \brief class for D-h correlation extraction /// \usage .L DhCorrelationExtraction.cxx+ /// \usage .x ExtractOutputCorrel.C("config-file-name") /// \author Samuele Cattaruzzi /// \author Swapnesh Santosh Khade +#include "DhCorrelationExtraction.h" #include "Riostream.h" + #include #include + #include #include -#include "DhCorrelationExtraction.h" + +#include +#include +#include +#include using namespace rapidjson; @@ -34,7 +41,7 @@ void readArray(const Value& jsonArray, std::vector& output) } } -void parseStringArray(const Value& jsonArray, std::vector& output) +void parseStringArray(const Value& jsonArray, std::vector& output) { size_t arrayLength = jsonArray.Size(); for (size_t i = 0; i < arrayLength; i++) { @@ -44,10 +51,13 @@ void parseStringArray(const Value& jsonArray, std::vector& output) } } -void SetInputCorrelNames(DhCorrelationExtraction* plotter, TString pathFileMass, TString pathFileSE, TString pathFileME, TString dirSE, TString dirME, TString histoNameCorrSignal, TString histoNameCorrSideba); -void SetInputHistoInvMassNames(DhCorrelationExtraction* plotter, std::vector inputMassNames); +void SetInputCorrelNames(DhCorrelationExtraction* plotter, TString pathFileSE, TString pathFileME, TString dirSE, TString dirME, TString histoNameCorrSignal, TString histoNameCorrSideba, TString histoNameCorrSidebaLeft, TString histoNameCorrSidebaRight); +void SetInputHistoInvMassNames(DhCorrelationExtraction* plotter, TString pathFileMass, std::vector inputMassNames); +void SetInputHistoFDSubtraction(DhCorrelationExtraction* plotter, TString pathFileFDTemplate, TString pathFileFDPromptFrac, TString histoNameFDTemplatePrompt, TString histoNameFDTemplateNonPrompt, TString histoNameRawFracPrompt); +void SetInputHistoSecPart(DhCorrelationExtraction* plotter, TString pathFileSecPart, TString dirSecPartName, TString histoNamePrimaryPart, TString histoNameAllPart); +void SetInputHistoBiasBtoD(DhCorrelationExtraction* plotter, TString pathfFilePromptMcRec, TString pathfFileNonPromptMcRec); -void ExtractOutputCorrel(TString cfgFileName = "config_CorrAnalysis.json") +void ExtractOutputCorrel_Ds(const TString cfgFileName = "config_CorrAnalysis.json") { // gStyle -> SetOptStat(0); gStyle->SetPadLeftMargin(0.15); @@ -72,24 +82,37 @@ void ExtractOutputCorrel(TString cfgFileName = "config_CorrAnalysis.json") string pathFileSE = config["pathFileSE"].GetString(); string pathFileME = config["pathFileME"].GetString(); string pathFileMass = config["pathFileMass"].GetString(); + string pathFileFDTemplate = config["pathFileFDTemplate"].GetString(); + string pathFileFDPromptFrac = config["pathFileFDPromptFrac"].GetString(); + string pathFileSecPart = config["pathFileSecPart"].GetString(); + string pathfFilePromptMcRec = config["pathfFilePromptMcRec"].GetString(); + string pathfFileNonPromptMcRec = config["pathfFileNonPromptMcRec"].GetString(); string dirSE = config["InputDirSE"].GetString(); string dirME = config["InputDirME"].GetString(); + string dirSecPart = config["InputDirSecPart"].GetString(); string histoNameCorrSignal = config["InputHistoCorrSignalName"].GetString(); string histoNameCorrSideba = config["InputHistoCorrSidebaName"].GetString(); + string histoNameCorrSidebaLeft = config["InputHistoCorrSidebaLeftName"].GetString(); + string histoNameCorrSidebaRight = config["InputHistoCorrSidebaRightName"].GetString(); + string histoNameFDTemplatePrompt = config["InputHistoFDTemplatePrompt"].GetString(); + string histoNameFDTemplateNonPrompt = config["InputHistoFDTemplateNonPrompt"].GetString(); + string histoNameRawFracPrompt = config["InputHistoFDPromptFrac"].GetString(); + string histoNamePrimaryPart = config["InputHistoPrimaryPart"].GetString(); + string histoNameAllPart = config["InputHistoAllPart"].GetString(); - vector InputHistoMassName; + std::vector InputHistoMassName; const Value& inputMassNames = config["InputHistoMassName"]; parseStringArray(inputMassNames, InputHistoMassName); - cout << InputHistoMassName[0].data() << endl; - cout << InputHistoMassName[1].data() << endl; - cout << InputHistoMassName[2].data() << endl; + std::cout << InputHistoMassName[0].data() << std::endl; + std::cout << InputHistoMassName[1].data() << std::endl; + std::cout << InputHistoMassName[2].data() << std::endl; - vector binsPtCandIntervals; - vector binsPtHadIntervals; - vector deltaEtaInterval; + std::vector binsPtCandIntervals; + std::vector binsPtHadIntervals; + std::vector deltaEtaInterval; const Value& PtCandValue = config["binsPtCandIntervals"]; readArray(PtCandValue, binsPtCandIntervals); @@ -103,9 +126,24 @@ void ExtractOutputCorrel(TString cfgFileName = "config_CorrAnalysis.json") double deltaEtaMax = deltaEtaInterval[1]; int specie = config["DmesonSpecie"].GetInt(); + bool rebinAngCorr = config["RebinAngCorr"].GetBool(); + bool rebinFDCorr = config["RebinFDCorr"].GetBool(); + bool rebinSecPart = config["RebinSecPart"].GetBool(); + int rebinDeltaPhi = config["nRebinDeltaPhi"].GetInt(); + int rebinDeltaEta = config["nRebinDeltaEta"].GetInt(); int npools = config["NumberOfPools"].GetInt(); bool poolByPool = config["CorrectPoolsSeparately"].GetBool(); + bool applySecPartCorr = config["ApplySecPartCorr"].GetBool(); + bool applyBiasBtoDCorr = config["ApplyBiasBtoDCorr"].GetBool(); + bool applyFDCorr = config["ApplyFDCorr"].GetBool(); + bool isDividedSideb = config["IsDividedSideb"].GetBool(); + bool useSidebLeft = config["UseSidebLeft"].GetBool(); + bool useSidebRight = config["UseSidebRight"].GetBool(); + + if (useSidebLeft && useSidebLeft) { + std::cout << "Using left and right" << std::endl; + } std::cout << "=========================== " << std::endl; std::cout << "Input variables from config" << std::endl; @@ -121,6 +159,9 @@ void ExtractOutputCorrel(TString cfgFileName = "config_CorrAnalysis.json") const int nBinsPtHad = binsPtHadIntervals.size() - 1; TH1D* hCorrectedCorrel[nBinsPtCand][nBinsPtHad]; + TH1D* hCorrectedCorrel_BaselineSubtr[nBinsPtCand][nBinsPtHad]; + TH1D* hCorrectedCorrel_Reflected[nBinsPtCand][nBinsPtHad]; + TH1D* hCorrectedCorrel_Reflected_BaselineSubtr[nBinsPtCand][nBinsPtHad]; // Create and set the correlation plotter class DhCorrelationExtraction* plotter = new DhCorrelationExtraction(); @@ -128,30 +169,55 @@ void ExtractOutputCorrel(TString cfgFileName = "config_CorrAnalysis.json") Bool_t flagSpecie = plotter->SetDmesonSpecie(static_cast(specie)); plotter->SetNpools(npools); plotter->SetCorrectPoolsSeparately(poolByPool); // kTRUE = pool.by-pool extraction and correction; kFALSE = merged ME pools + plotter->SetFDSubtraction(applyFDCorr); + plotter->SetSecPartContamination(applySecPartCorr); plotter->SetDeltaEtaRange(deltaEtaMin, deltaEtaMax); plotter->SetSubtractSoftPiInMEdistr(kFALSE); - plotter->SetRebin2DcorrelHisto(2, 2); // Xaxis: deltaEta, Yaxis: deltaPhi + plotter->SetRebinOptions(rebinAngCorr, rebinFDCorr, rebinSecPart); + plotter->SetRebin2DcorrelHisto(rebinDeltaEta, rebinDeltaPhi); // Xaxis: deltaEta, Yaxis: deltaPhi + plotter->SetCorrBiasBtoD(applyBiasBtoDCorr); plotter->SetDebugLevel(1); if (!flagSpecie) - cout << "[ERROR] Wrong D meson flag" << endl; + std::cout << "[ERROR] Wrong D meson flag" << std::endl; // Set the input file config - SetInputCorrelNames(plotter, pathFileMass, pathFileSE, pathFileME, dirSE, dirME, histoNameCorrSignal, histoNameCorrSideba); - SetInputHistoInvMassNames(plotter, InputHistoMassName); + SetInputCorrelNames(plotter, pathFileSE, pathFileME, dirSE, dirME, histoNameCorrSignal, histoNameCorrSideba, histoNameCorrSidebaLeft, histoNameCorrSidebaRight); + SetInputHistoInvMassNames(plotter, pathFileMass, InputHistoMassName); + if (applyFDCorr) + SetInputHistoFDSubtraction(plotter, pathFileFDTemplate, pathFileFDPromptFrac, histoNameFDTemplatePrompt, histoNameFDTemplateNonPrompt, histoNameRawFracPrompt); + if (applySecPartCorr) + SetInputHistoSecPart(plotter, pathFileSecPart, dirSecPart, histoNamePrimaryPart, histoNameAllPart); + if (applyBiasBtoDCorr) + SetInputHistoBiasBtoD(plotter, pathfFilePromptMcRec, pathfFileNonPromptMcRec); Bool_t readSEandME = plotter->ReadInputSEandME(); - Bool_t readInvMass = plotter->ReadInputInvMass(); if (readSEandME) - cout << "Files SE and ME read correctly" << endl; + std::cout << "Files SE and ME read correctly" << std::endl; + Bool_t readInvMass = plotter->ReadInputInvMass(); if (readInvMass) - cout << "Files inv. mass read correctly" << endl; + std::cout << "Files inv. mass read correctly" << std::endl; + if (applyFDCorr) { + Bool_t readFDSubtr = plotter->ReadInputFDSubtr(); + if (readFDSubtr) + std::cout << "Files for FD subtr. read correctly" << std::endl; + } + if (applySecPartCorr) { + Bool_t readSecPart = plotter->ReadInputSecondaryPartContamination(); + if (readSecPart) + std::cout << "Files for secondary part. contamination read correctly" << std::endl; + } // Loop over candidate pt and assoc. particle pt for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { + plotter->SetDividedSidebands(isDividedSideb, useSidebLeft, useSidebRight); plotter->GetSignalAndBackgroundForNorm(binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1]); for (int iBinPtHad = 0; iBinPtHad < nBinsPtHad; iBinPtHad++) { + plotter->SetBinCandAndHad(iBinPtCand + 1, iBinPtHad + 1); plotter->ExtractCorrelations(binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1], binsPtHadIntervals[iBinPtHad], binsPtHadIntervals[iBinPtHad + 1], CodeNameAnalysis); - hCorrectedCorrel[iBinPtCand][iBinPtHad] = reinterpret_cast(plotter->GetCorrectedCorrHisto()); + hCorrectedCorrel[iBinPtCand][iBinPtHad] = (TH1D*)plotter->GetCorrectedCorrHisto(); + hCorrectedCorrel_BaselineSubtr[iBinPtCand][iBinPtHad] = (TH1D*)plotter->GetCorrectedCorrHisto_BaselineSubtr(); + hCorrectedCorrel_Reflected[iBinPtCand][iBinPtHad] = (TH1D*)plotter->GetCorrectedCorrHisto_Reflected(); + hCorrectedCorrel_Reflected_BaselineSubtr[iBinPtCand][iBinPtHad] = (TH1D*)plotter->GetCorrectedCorrHisto_Reflected_BaselineSubtr(); } } @@ -165,14 +231,43 @@ void ExtractOutputCorrel(TString cfgFileName = "config_CorrAnalysis.json") } outFile->Close(); + // output file baseline subtr. + TFile* outFile_BaselineSubtr = new TFile(Form("Output_CorrelationExtraction_%s_Root/ExtractCorrelationsResults_BaselineSubtr.root", CodeNameAnalysis.data()), "RECREATE"); + outFile_BaselineSubtr->cd(); + for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { + for (int iBinPtHad = 0; iBinPtHad < nBinsPtHad; iBinPtHad++) { + hCorrectedCorrel_BaselineSubtr[iBinPtCand][iBinPtHad]->Write(); + } + } + outFile_BaselineSubtr->Close(); + + // output file reflected + TFile* outFile_Reflected = new TFile(Form("Output_CorrelationExtraction_%s_Root/ExtractCorrelationsResults_Reflected.root", CodeNameAnalysis.data()), "RECREATE"); + outFile_Reflected->cd(); + for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { + for (int iBinPtHad = 0; iBinPtHad < nBinsPtHad; iBinPtHad++) { + hCorrectedCorrel_Reflected[iBinPtCand][iBinPtHad]->Write(); + } + } + outFile_Reflected->Close(); + + // output file reflected baseline subtr. + TFile* outFile_Reflected_BaselineSubtr = new TFile(Form("Output_CorrelationExtraction_%s_Root/ExtractCorrelationsResults_Reflected_BaselineSubtr.root", CodeNameAnalysis.data()), "RECREATE"); + outFile_Reflected_BaselineSubtr->cd(); + for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { + for (int iBinPtHad = 0; iBinPtHad < nBinsPtHad; iBinPtHad++) { + hCorrectedCorrel_Reflected_BaselineSubtr[iBinPtCand][iBinPtHad]->Write(); + } + } + outFile_Reflected_BaselineSubtr->Close(); + return; } -void SetInputCorrelNames(DhCorrelationExtraction* plotter, TString pathFileMass, TString pathFileSE, TString pathFileME, TString dirSE, TString dirME, TString histoNameCorrSignal, TString histoNameCorrSideba) +void SetInputCorrelNames(DhCorrelationExtraction* plotter, TString pathFileSE, TString pathFileME, TString dirSE, TString dirME, TString histoNameCorrSignal, TString histoNameCorrSideba, TString histoNameCorrSidebaLeft, TString histoNameCorrSidebaRight) { - // paths - plotter->SetInputFilenameMass(pathFileMass.Data()); + // Ds paths plotter->SetInputFilenameSE(pathFileSE.Data()); plotter->SetInputFilenameME(pathFileME.Data()); plotter->SetDirNameSE(dirSE.Data()); @@ -181,16 +276,51 @@ void SetInputCorrelNames(DhCorrelationExtraction* plotter, TString pathFileMass, plotter->SetSECorrelHistoSidebandName(histoNameCorrSideba.Data()); plotter->SetMECorrelHistoSignalName(histoNameCorrSignal.Data()); plotter->SetMECorrelHistoSidebandName(histoNameCorrSideba.Data()); + plotter->SetSECorrelHistoSidebandLeftName(histoNameCorrSidebaLeft.Data()); + plotter->SetMECorrelHistoSidebandLeftName(histoNameCorrSidebaLeft.Data()); + plotter->SetSECorrelHistoSidebandRightName(histoNameCorrSidebaRight.Data()); + plotter->SetMECorrelHistoSidebandRightName(histoNameCorrSidebaRight.Data()); return; } -void SetInputHistoInvMassNames(DhCorrelationExtraction* plotter, std::vector inputMassNames) +void SetInputHistoInvMassNames(DhCorrelationExtraction* plotter, TString pathFileMass, std::vector inputMassNames) { // to use if sgn and bkg extraction is done apart + plotter->SetInputFilenameMass(pathFileMass.Data()); plotter->SetMassHistoNameSgn(inputMassNames[0].data()); plotter->SetMassHistoNameBkg(inputMassNames[1].data()); plotter->SetMassHistoNameSBs(inputMassNames[2].data()); return; } + +void SetInputHistoFDSubtraction(DhCorrelationExtraction* plotter, TString pathFileFDTemplate, TString pathFileFDPromptFrac, TString histoNameFDTemplatePrompt, TString histoNameFDTemplateNonPrompt, TString histoNameRawFracPrompt) +{ + + plotter->SetInputFilenameFDTemplate(pathFileFDTemplate.Data()); + plotter->SetInputFilenameFDPromptFrac(pathFileFDPromptFrac.Data()); + plotter->SetInputHistoNameFDTemplatePrompt(histoNameFDTemplatePrompt.Data()); + plotter->SetInputHistoNameFDTemplateNonPrompt(histoNameFDTemplateNonPrompt.Data()); + plotter->SetInputHistoNameFDPromptFrac(histoNameRawFracPrompt.Data()); + + return; +} + +void SetInputHistoSecPart(DhCorrelationExtraction* plotter, TString pathFileSecPart, TString dirSecPartName, TString histoNamePrimaryPart, TString histoNameAllPart) +{ + + plotter->SetInputFilenameSecPart(pathFileSecPart.Data()); + plotter->SetDirNameSecPart(dirSecPartName.Data()); + plotter->SetHistoSecPartName(histoNamePrimaryPart.Data(), histoNameAllPart.Data()); + + return; +} + +void SetInputHistoBiasBtoD(DhCorrelationExtraction* plotter, TString pathfFilePromptMcRec, TString pathfFileNonPromptMcRec) +{ + + plotter->SetInputFilenameBiasBtoD(pathfFilePromptMcRec.Data(), pathfFileNonPromptMcRec.Data()); + + return; +} diff --git a/PWGHF/HFC/Macros/FitCorrel.C b/PWGHF/HFC/Macros/FitCorrel.C index 37b36a5e1b3..7c1d64ea823 100644 --- a/PWGHF/HFC/Macros/FitCorrel.C +++ b/PWGHF/HFC/Macros/FitCorrel.C @@ -9,29 +9,37 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file FitCorrel.C -/// \brief Macro to perform the azimuthal correlation fit -/// \usage .L DhCorrelationFitter.cxx+ -/// \usage .x FitCorrel.C("config-file-name") +/// \file DhCorrelationExtraction.cxx +/// \brief class for D-h correlation extraction /// \author Samuele Cattaruzzi /// \author Swapnesh Santosh Khade +#include "DhCorrelationFitter.h" #include "Riostream.h" -#include + #include #include #include #include -#include +#include #include +#include + #include #include -#include "DhCorrelationFitter.h" +#include +#include +#include +#include + +using namespace std; using namespace rapidjson; +bool removeNSPeakLowPt = false; + template -void readArray(const Value& jsonArray, std::vector& output) +void readArray(const Value& jsonArray, vector& output) { for (auto it = jsonArray.Begin(); it != jsonArray.End(); it++) { auto value = it->template Get(); @@ -41,15 +49,21 @@ void readArray(const Value& jsonArray, std::vector& output) void SetTH1HistoStyle(TH1D*& histo, TString hTitle, TString hXaxisTitle, TString hYaxisTitle, Style_t markerStyle, Color_t markerColor, Double_t markerSize, - Color_t lineColor, Int_t lineWidth, Float_t hTitleXaxisOffset = 1., Float_t hTitleYaxisOffset = 1., - Float_t hTitleXaxisSize = 0.060, Float_t hTitleYaxisSize = 0.060, Float_t hLabelXaxisSize = 0.060, Float_t hLabelYaxisSize = 0.060, + Color_t lineColor, Int_t lineWidth, Float_t hTitleXaxisOffset = 1.3, Float_t hTitleYaxisOffset = 1.3, + Float_t hTitleXaxisSize = 0.045, Float_t hTitleYaxisSize = 0.045, Float_t hLabelXaxisSize = 0.045, Float_t hLabelYaxisSize = 0.045, + Bool_t centerXaxisTitle = false, Bool_t centerYaxisTitle = false); +void SetTH1HistoStyle(TH1F*& histo, TString hTitle, TString hXaxisTitle, TString hYaxisTitle, + Style_t markerStyle, Color_t markerColor, Double_t markerSize, + Color_t lineColor, Int_t lineWidth, Float_t hTitleXaxisOffset = 1.3, Float_t hTitleYaxisOffset = 1.3, + Float_t hTitleXaxisSize = 0.045, Float_t hTitleYaxisSize = 0.045, Float_t hLabelXaxisSize = 0.045, Float_t hLabelYaxisSize = 0.045, Bool_t centerXaxisTitle = false, Bool_t centerYaxisTitle = false); -void FitCorrel(TString cfgFileName = "config_CorrAnalysis.json") +void FitCorrel_Ds(const TString cfgFileName = "config_CorrAnalysis.json") { gStyle->SetOptStat(0); - gStyle->SetPadLeftMargin(0.15); - gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadLeftMargin(0.2); + gStyle->SetPadRightMargin(0.005); + gStyle->SetPadBottomMargin(0.2); gStyle->SetFrameLineWidth(2); gStyle->SetLineWidth(2); gStyle->SetCanvasDefH(1126); @@ -67,24 +81,57 @@ void FitCorrel(TString cfgFileName = "config_CorrAnalysis.json") gSystem->Exec(Form("rm -rf Output_CorrelationFitting_%s_Root/ Output_CorrelationFitting_%s_png/", CodeNameAnalysis.data(), CodeNameAnalysis.data())); gSystem->Exec(Form("mkdir Output_CorrelationFitting_%s_Root/ Output_CorrelationFitting_%s_png/", CodeNameAnalysis.data(), CodeNameAnalysis.data())); - const TString inFileName = Form("Output_CorrelationExtraction_%s_Root/ExtractCorrelationsResults.root", CodeNameAnalysis.data()); + string inputFileNameFit = config["InputFileNameFitCorr"].GetString(); + const TString inFileName = Form("Output_CorrelationExtraction_%s_Root/%s", CodeNameAnalysis.data(), inputFileNameFit.data()); - vector binsPtCandIntervals; - vector binsPtHadIntervals; + bool isReflected = config["IsRiflected"].GetBool(); + bool drawSystematicErrors = config["DrawSystematics"].GetBool(); + bool sameSystematics = config["SameSystematics"].GetBool(); + bool shiftBaseUp = config["ShiftBaseUp"].GetBool(); + bool shiftBaseDown = config["ShiftBaseDown"].GetBool(); + + std::vector binsPtCandIntervalsVec; + std::vector binsPtHadIntervals; + std::vector fitFunc; const Value& PtCandValue = config["binsPtCandIntervals"]; - readArray(PtCandValue, binsPtCandIntervals); + readArray(PtCandValue, binsPtCandIntervalsVec); const Value& PtHadValue = config["binsPtHadIntervals"]; readArray(PtHadValue, binsPtHadIntervals); - int fitFunc = config["FitFunction"].GetInt(); + const int nBinsPtCand = binsPtCandIntervalsVec.size() - 1; + const int nBinsPtHad = binsPtHadIntervals.size() - 1; + + double binsPtCandIntervals[nBinsPtCand + 1]; + for (int i = 0; i < nBinsPtCand + 1; i++) { + binsPtCandIntervals[i] = binsPtCandIntervalsVec[i]; + } + + const Value& FitFuncValue = config["FitFunction"]; + readArray(FitFuncValue, fitFunc); + int fixBase = config["FixBaseline"].GetInt(); int fixMean = config["FixMean"].GetInt(); + int nBaselinePoints = config["nBaselinePoints"].GetInt(); + vector pointsForBaselineVec; + const Value& pointsForBaselineValue = config["binsForBaseline"]; + readArray(pointsForBaselineValue, pointsForBaselineVec); + if (pointsForBaselineVec.size() != nBaselinePoints) { + cout << "ERROR: size of the vector pointsForBaseline is different from the number of nBaselinePoints" << endl; + return; + } + int pointsForBaseline[nBaselinePoints]; + for (int i = 0; i < nBaselinePoints; i++) { + pointsForBaseline[i] = pointsForBaselineVec[i]; + } + std::cout << "=========================== " << std::endl; std::cout << "Input variables from config" << std::endl; - std::cout << "FitFunction = " << fitFunc << std::endl; + for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { + std::cout << "iPt = " << iBinPtCand + 1 << " FitFunction = " << fitFunc[iBinPtCand] << std::endl; + } std::cout << "FixBaseline = " << fixBase << std::endl; std::cout << "FixMean = " << fixMean << std::endl; std::cout << "=========================== " << std::endl; @@ -93,44 +140,65 @@ void FitCorrel(TString cfgFileName = "config_CorrAnalysis.json") // TODO: reflections bool refl = false; - const int nBinsPtCand = binsPtCandIntervals.size() - 1; - const int nBinsPtHad = binsPtHadIntervals.size() - 1; - // Input file TFile* inFile = new TFile(inFileName.Data()); + TFile* inFileSystematicErrors = new TFile("OutputSystematicUncertainties/SystematicUncertaintesAngCorrMerged.root"); + TFile* inFileFitSystematicErrors = new TFile("OutputSystematicUncertainties/SystematicUncertaintesFitPhysObsMerged.root"); // Canvas TCanvas* CanvasCorrPhi[nBinsPtHad]; // Histograms TH1D* hCorrPhi[nBinsPtCand][nBinsPtHad]; + TH1F* hSystematicErrors[nBinsPtCand][nBinsPtHad]; + TH1D* hSystematicErrorsPlot[nBinsPtCand][nBinsPtHad]; - int nBinsPhi; - double baselineFromThreePoints[nBinsPtCand][nBinsPtHad], baselineFromThreePointsError[nBinsPtCand][nBinsPtHad]; + const int nBinsPtD = 5; + if (nBinsPtD != nBinsPtCand) { + std::cout << "[ERROR]: nBinsPtD != nBinsPtCand" << std::endl; + return; + } + double SystUncCorrelatedDs[nBinsPtD] = {20, 20, 20, 10}; // % (just the MC Closure uncertainty to put in the plot) // DhCorrelationFitter const double fMin{-0.5 * TMath::Pi()}, fMax{1.5 * TMath::Pi()}; // limits for the fitting function DhCorrelationFitter* corrFitter[nBinsPtHad][nBinsPtCand]; // Input parameters for fitting - const int npars{8}; // PED NSY NSM NSW ASY ASM ASW BETA - Double_t vals[npars] = {3., 2., 0., 0.5, 2., 3.14, 0.3, 2.}; - Double_t lowBounds[npars] = {0., 0., -1., 0., 0., 2., 0., 0.5}; - Double_t uppBounds[npars] = {9999., 999., 1., 3.14 / 3., 999., 4., 3.14 / 2., 3.5}; - - const int nBaselinePoints{8}; - Int_t pointsForBaseline[nBaselinePoints] = {1, 2, 13, 14, 15, 16, 31, 32}; + const int npars{10}; // PED NSY NSM NSW ASY ASM ASW BETA v2D v2h + Double_t vals[npars] = {3., 2., 0., 0.5, 2., 3.14, 0.3, 2., 0.1, 0.1}; + Double_t lowBounds[npars] = {0., 0., -1., 0., 0., 2., 0., 0.5, 0., 0.}; + Double_t uppBounds[npars] = {9999., 999., 1., 3.14 / 3., 999., 4., 3.14 / 2., 3.5, 0.5, 0.5}; + + Double_t v2AssocPart[nBinsPtD] = {0.15, 0.15, 0.15, 0.15}; + Double_t v2Dmeson[nBinsPtD] = {0.175, 0.09, 0.04, 0.04}; + + // Output histograms + TH1D* hBaselin[nBinsPtHad]; + TH1D* hNSYield[nBinsPtHad]; + TH1D* hNSSigma[nBinsPtHad]; + TH1D* hASYield[nBinsPtHad]; + TH1D* hASSigma[nBinsPtHad]; + TH1D* hBeta[nBinsPtHad]; + TH1D* hNSYieldBinCount[nBinsPtHad]; + TH1D* hASYieldBinCount[nBinsPtHad]; // extract TH1D and prepare fit for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { for (int iBinPtHad = 0; iBinPtHad < nBinsPtHad; iBinPtHad++) { - hCorrPhi[iBinPtCand][iBinPtHad] = reinterpret_cast(inFile->Get(Form("hCorrectedCorr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1], binsPtHadIntervals[iBinPtHad], binsPtHadIntervals[iBinPtHad + 1]))); + if (isReflected) { + hCorrPhi[iBinPtCand][iBinPtHad] = reinterpret_cast(inFile->Get(Form("hCorrectedCorrReflected_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1], binsPtHadIntervals[iBinPtHad], binsPtHadIntervals[iBinPtHad + 1]))); + } else { + hCorrPhi[iBinPtCand][iBinPtHad] = reinterpret_cast(inFile->Get(Form("hCorrectedCorr_PtCand%.0fto%.0f_PtAssoc%.0fto%.0f", binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1], binsPtHadIntervals[iBinPtHad], binsPtHadIntervals[iBinPtHad + 1]))); + } - corrFitter[iBinPtHad][iBinPtCand] = new DhCorrelationFitter(reinterpret_cast(hCorrPhi[iBinPtCand][iBinPtHad], fMin, fMax)); + corrFitter[iBinPtHad][iBinPtCand] = new DhCorrelationFitter(reinterpret_cast(hCorrPhi[iBinPtCand][iBinPtHad]), fMin, fMax); corrFitter[iBinPtHad][iBinPtCand]->SetHistoIsReflected(refl); - corrFitter[iBinPtHad][iBinPtCand]->SetFuncType(static_cast(fitFunc)); corrFitter[iBinPtHad][iBinPtCand]->SetFixBaseline(fixBase); + corrFitter[iBinPtHad][iBinPtCand]->SetBaselineUpOrDown(shiftBaseUp, shiftBaseDown); corrFitter[iBinPtHad][iBinPtCand]->SetPointsForBaseline(nBaselinePoints, pointsForBaseline); + corrFitter[iBinPtHad][iBinPtCand]->Setv2(v2AssocPart[iBinPtCand], v2Dmeson[iBinPtCand]); + corrFitter[iBinPtHad][iBinPtCand]->SetReflectedCorrHisto(isReflected); corrFitter[iBinPtHad][iBinPtCand]->SetFixMean(fixMean); corrFitter[iBinPtHad][iBinPtCand]->SetPtRanges(binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1], binsPtHadIntervals[iBinPtHad], binsPtHadIntervals[iBinPtHad + 1]); @@ -145,15 +213,50 @@ void FitCorrel(TString cfgFileName = "config_CorrAnalysis.json") CanvasCorrPhi[iBinPtHad]->Divide(2, 2); } if (nBinsPtCand > 4 && nBinsPtCand <= 6) { - CanvasCorrPhi[iBinPtHad]->Divide(2, 3); + CanvasCorrPhi[iBinPtHad]->Divide(3, 2); } + // histograms with fir parameters + hBaselin[iBinPtHad] = new TH1D(Form("hBaselin_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + hNSYield[iBinPtHad] = new TH1D(Form("hNSYield_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + hNSSigma[iBinPtHad] = new TH1D(Form("hNSSigma_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + hASYield[iBinPtHad] = new TH1D(Form("hASYield_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + hASSigma[iBinPtHad] = new TH1D(Form("hASSigma_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + hBeta[iBinPtHad] = new TH1D(Form("hBeta_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + hNSYieldBinCount[iBinPtHad] = new TH1D(Form("hNSYieldBinCount_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + hASYieldBinCount[iBinPtHad] = new TH1D(Form("hASYieldBinCount_PtBinAssoc%d", iBinPtHad + 1), "", nBinsPtCand, binsPtCandIntervals); + for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { SetTH1HistoStyle(hCorrPhi[iBinPtCand][iBinPtHad], "", "#Delta#phi [rad]", "#frac{1}{N_{D_{s}}}#frac{dN^{assoc}}{d#Delta#phi} [rad^{-1}]", kFullCircle, kRed + 1, 1.4, kRed + 1, 3); CanvasCorrPhi[iBinPtHad]->cd(iBinPtCand + 1); - hCorrPhi[iBinPtCand][iBinPtHad]->Draw(); + CanvasCorrPhi[iBinPtHad]->SetTickx(); + CanvasCorrPhi[iBinPtHad]->SetTicky(); + hCorrPhi[iBinPtCand][iBinPtHad]->SetStats(0); + hCorrPhi[iBinPtCand][iBinPtHad]->SetMinimum(0); + // hCorrPhi[iBinPtCand][iBinPtHad] -> Draw(); + + // draw systematic errors + int nBinsPhi = hCorrPhi[iBinPtCand][iBinPtHad]->GetNbinsX(); + if (drawSystematicErrors) { + hSystematicErrors[iBinPtCand][iBinPtHad] = reinterpret_cast(inFileSystematicErrors->Get(Form("hSystematicErrorsMerged_PtBin%d_PtBinAssoc%d", iBinPtCand + 1, iBinPtHad + 1))); + hSystematicErrorsPlot[iBinPtCand][iBinPtHad] = reinterpret_cast(hCorrPhi[iBinPtCand][iBinPtHad]->Clone(Form("hSystematicErrorsPlot_PtBin%d_PtBinAssoc%d", iBinPtCand + 1, iBinPtHad + 1))); + for (int iPhi = 0; iPhi < nBinsPhi; iPhi++) { + if (sameSystematics) { + hSystematicErrorsPlot[iBinPtCand][iBinPtHad]->SetBinError(iPhi + 1, hSystematicErrors[iBinPtCand][iBinPtHad]->GetBinContent(1) * hCorrPhi[iBinPtCand][iBinPtHad]->GetBinContent(iPhi + 1)); + } else { + hSystematicErrorsPlot[iBinPtCand][iBinPtHad]->SetBinError(iPhi + 1, hSystematicErrors[iBinPtCand][iBinPtHad]->GetBinContent(iPhi + 1) * hCorrPhi[iBinPtCand][iBinPtHad]->GetBinContent(iPhi + 1)); + } + hSystematicErrorsPlot[iBinPtCand][iBinPtHad]->SetBinContent(iPhi + 1, hCorrPhi[iBinPtCand][iBinPtHad]->GetBinContent(iPhi + 1)); + } + SetTH1HistoStyle(hSystematicErrorsPlot[iBinPtCand][iBinPtHad], "", "#Delta#phi [rad]", "#frac{1}{N_{D_{s}}}#frac{dN^{assoc}}{d#Delta#phi} [rad^{-1}]", kFullCircle, kRed + 1, 1.4, kRed + 1, 2); + hSystematicErrorsPlot[iBinPtCand][iBinPtHad]->SetLineColor(kRed - 4); + hSystematicErrorsPlot[iBinPtCand][iBinPtHad]->SetFillStyle(0); + // hSystematicErrorsPlot[iBinPtCand][iBinPtHad] -> Draw("E2same"); + } + // hCorrPhi[iBinPtCand][iBinPtHad] -> Draw("same"); // Fit + corrFitter[iBinPtHad][iBinPtCand]->SetFuncType(static_cast(fitFunc[iBinPtCand])); corrFitter[iBinPtHad][iBinPtCand]->Fitting(kTRUE, kTRUE); // the first term is for drawing the fit functions, the second argument is useExternalParams TF1* fFit = corrFitter[iBinPtHad][iBinPtCand]->GetFitFunction(); @@ -162,11 +265,193 @@ void FitCorrel(TString cfgFileName = "config_CorrAnalysis.json") TPaveText* pttext = new TPaveText(0.15, 0.9, 0.85, 0.95, "NDC"); pttext->SetFillStyle(0); pttext->SetBorderSize(0); - TText* tpT = pttext->AddText(0., 0.8, Form("%.0f < p_{T}^{D_{s}} < %.0f GeV/c, p_{T}^{assoc} > 0.3 GeV/c", binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1])); + TText* tpT = pttext->AddText(0., 0.8, Form("%.0f < p_{T}^{D_{s}} < %.0f GeV/c, p_{T}^{assoc} > %.1f GeV/c", binsPtCandIntervals[iBinPtCand], binsPtCandIntervals[iBinPtCand + 1], binsPtHadIntervals[iBinPtHad])); + // pttext -> Draw("same"); + + // Fill the histograms with the fit parameters + hBaselin[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetPedestal()); + hBaselin[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetPedestalError()); + if (iBinPtCand == 0 && removeNSPeakLowPt) { + hNSYield[iBinPtHad]->SetBinContent(iBinPtCand + 1, -1); + hNSYield[iBinPtHad]->SetBinError(iBinPtCand + 1, 0); + + hNSSigma[iBinPtHad]->SetBinContent(iBinPtCand + 1, -1); + hNSSigma[iBinPtHad]->SetBinError(iBinPtCand + 1, 0); + + hBeta[iBinPtHad]->SetBinContent(iBinPtCand + 1, -1); + hBeta[iBinPtHad]->SetBinError(iBinPtCand + 1, 0); + } else { + hNSYield[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetNSYield()); + hNSYield[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetNSYieldError()); + + if (fitFunc[iBinPtCand] != 5 && fitFunc[iBinPtCand] != 6) { + hNSSigma[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetNSSigma()); + hNSSigma[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetNSSigmaError()); + } else { + hNSSigma[iBinPtHad]->SetBinContent(iBinPtCand + 1, TMath::Sqrt(1. / corrFitter[iBinPtHad][iBinPtCand]->GetNSSigma())); + Double_t errrel = corrFitter[iBinPtHad][iBinPtCand]->GetNSSigmaError() / corrFitter[iBinPtHad][iBinPtCand]->GetNSSigma() / 2.; + hNSSigma[iBinPtHad]->SetBinError(iBinPtCand + 1, errrel * TMath::Sqrt(1. / corrFitter[iBinPtHad][iBinPtCand]->GetNSSigma())); + } + } + hNSYieldBinCount[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetBinCountingNSYield()); + hNSYieldBinCount[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetBinCountingNSYieldErr()); + + hASYield[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetASYield()); + hASYield[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetASYieldError()); + + hASYieldBinCount[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetBinCountingASYield()); + hASYieldBinCount[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetBinCountingASYieldErr()); + if (fitFunc[iBinPtCand] != 5 && fitFunc[iBinPtCand] != 6) { + hASSigma[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetASSigma()); + hASSigma[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetASSigmaError()); + } else { + hASSigma[iBinPtHad]->SetBinContent(iBinPtCand + 1, TMath::Sqrt(1. / corrFitter[iBinPtHad][iBinPtCand]->GetASSigma())); + Double_t errrel = corrFitter[iBinPtHad][iBinPtCand]->GetASSigmaError() / corrFitter[iBinPtHad][iBinPtCand]->GetASSigma() / 2.; + hASSigma[iBinPtHad]->SetBinError(iBinPtCand + 1, errrel * TMath::Sqrt(1. / corrFitter[iBinPtHad][iBinPtCand]->GetASSigma())); + } + if (fitFunc[iBinPtCand] == 4) { // param beta for gen. gauss + hBeta[iBinPtHad]->SetBinContent(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetBeta()); + hBeta[iBinPtHad]->SetBinError(iBinPtCand + 1, corrFitter[iBinPtHad][iBinPtCand]->GetBetaError()); + } + + // Draw + TPaveText* tCorrUncDs = new TPaveText(0.413, 0.311, 0.877, 0.392, "NDC"); + tCorrUncDs->SetFillStyle(0); + tCorrUncDs->SetBorderSize(0); + tCorrUncDs->SetTextSize(0.05); + tCorrUncDs->SetTextFont(42); + tCorrUncDs->SetTextAlign(13); + tCorrUncDs->SetTextColor(kRed + 1); + tCorrUncDs->AddText(0., 0., Form("#splitline{+%.0f%%}{#minus%.0f%%}", SystUncCorrelatedDs[iBinPtCand], SystUncCorrelatedDs[iBinPtCand])); + + TPaveText* tScaleUnc = new TPaveText(0.501, 0.292, 0.968, 0.372, "NDC"); + tScaleUnc->SetFillStyle(0); + tScaleUnc->SetBorderSize(0); + tScaleUnc->SetTextSize(0.05); + tScaleUnc->SetTextFont(42); + tScaleUnc->SetTextAlign(13); + tScaleUnc->SetTextColor(kBlack); + tScaleUnc->AddText(0., 0., "corr. unc."); + + if (drawSystematicErrors) + hSystematicErrorsPlot[iBinPtCand][iBinPtHad]->Draw("E2same"); + hCorrPhi[iBinPtCand][iBinPtHad]->Draw("same"); pttext->Draw("same"); + if (drawSystematicErrors) + tCorrUncDs->Draw("same"); + if (drawSystematicErrors) + tScaleUnc->Draw("same"); } - CanvasCorrPhi[iBinPtHad]->SaveAs(Form("Output_CorrelationFitting_%s_png/CorrPhi_PtBinAssoc%d.png", CodeNameAnalysis.data(), iBinPtHad + 1)); - CanvasCorrPhi[iBinPtHad]->SaveAs(Form("Output_CorrelationFitting_%s_Root/CorrPhi_PtBinAssoc%d.root", CodeNameAnalysis.data(), iBinPtHad + 1)); + CanvasCorrPhi[iBinPtHad]->SaveAs(Form("Output_CorrelationFitting_%s_png/CorrPhiDs_PtBinAssoc%d.png", CodeNameAnalysis.data(), iBinPtHad + 1)); + CanvasCorrPhi[iBinPtHad]->SaveAs(Form("Output_CorrelationFitting_%s_Root/CorrPhiDs_PtBinAssoc%d.root", CodeNameAnalysis.data(), iBinPtHad + 1)); + } + + // histogram with fit parameter and errors + TFile* outFile = new TFile(Form("Output_CorrelationFitting_%s_Root/CorrPhiDs_FinalPlots.root", CodeNameAnalysis.data()), "RECREATE"); + outFile->cd(); + for (int iBinPtHad = 0; iBinPtHad < nBinsPtHad; iBinPtHad++) { + hBaselin[iBinPtHad]->Write(); + hNSYield[iBinPtHad]->Write(); + hNSSigma[iBinPtHad]->Write(); + hASYield[iBinPtHad]->Write(); + hASSigma[iBinPtHad]->Write(); + hBeta[iBinPtHad]->Write(); + hNSYieldBinCount[iBinPtHad]->Write(); + hASYieldBinCount[iBinPtHad]->Write(); + } + outFile->Close(); + + // Draw plots + for (int iBinPtHad = 0; iBinPtHad < nBinsPtHad; iBinPtHad++) { + TCanvas* c1 = new TCanvas(Form("NS_yield_PtAssoc%d", iBinPtHad + 1), Form("NS_yield_PtAssoc%d", iBinPtHad + 1)); + TCanvas* c2 = new TCanvas(Form("AS_yield_PtAssoc%d", iBinPtHad + 1), Form("AS_yield_PtAssoc%d", iBinPtHad + 1)); + TCanvas* c3 = new TCanvas(Form("NS_sigma_PtAssoc%d", iBinPtHad + 1), Form("AS_sigma_PtAssoc%d", iBinPtHad + 1)); + TCanvas* c4 = new TCanvas(Form("AS_sigma_PtAssoc%d", iBinPtHad + 1), Form("AS_sigma_PtAssoc%d", iBinPtHad + 1)); + TCanvas* c5 = new TCanvas(Form("Baseline_PtAssoc%d", iBinPtHad + 1), Form("Baseline_PtAssoc%d", iBinPtHad + 1)); + SetTH1HistoStyle(hBaselin[iBinPtHad], Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "Baseline", kFullSquare, kBlue, 1.8, kBlue, 2); + SetTH1HistoStyle(hNSYield[iBinPtHad], Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "Y^{NS}", kFullSquare, kRed, 1.8, kRed, 2); + SetTH1HistoStyle(hASYield[iBinPtHad], Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "Y^{AS}", kFullSquare, kMagenta, 1.8, kMagenta, 2); + SetTH1HistoStyle(hNSSigma[iBinPtHad], Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "#sigma_{NS}", kFullSquare, kOrange + 8, 1.8, kOrange + 8, 2); + SetTH1HistoStyle(hASSigma[iBinPtHad], Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "#sigma_{AS}", kFullSquare, kViolet - 5, 1.8, kViolet - 5, 2); + c1->cd(); + hNSYield[iBinPtHad]->SetMinimum(0); + hNSYield[iBinPtHad]->Draw(); + c2->cd(); + hASYield[iBinPtHad]->SetMinimum(0); + hASYield[iBinPtHad]->Draw(); + c3->cd(); + hNSSigma[iBinPtHad]->SetMinimum(0); + hNSSigma[iBinPtHad]->Draw(); + c4->cd(); + hASSigma[iBinPtHad]->SetMinimum(0); + hASSigma[iBinPtHad]->Draw(); + c5->cd(); + hBaselin[iBinPtHad]->SetMinimum(0); + hBaselin[iBinPtHad]->Draw(); + + TH1F* hBaselinSyst; + TH1F* hNSYieldSyst; + TH1F* hNSSigmaSyst; + TH1F* hASYieldSyst; + TH1F* hASSigmaSyst; + + if (drawSystematicErrors) { + hBaselinSyst = reinterpret_cast(inFileFitSystematicErrors->Get(Form("hSystematicErrorsBaselinMerged_PtBinAssoc%d", iBinPtHad + 1))); + hNSYieldSyst = reinterpret_cast(inFileFitSystematicErrors->Get(Form("hSystematicErrorsNSYieldMerged_PtBinAssoc%d", iBinPtHad + 1))); + hNSSigmaSyst = reinterpret_cast(inFileFitSystematicErrors->Get(Form("hSystematicErrorsNSSigmaMerged_PtBinAssoc%d", iBinPtHad + 1))); + hASYieldSyst = reinterpret_cast(inFileFitSystematicErrors->Get(Form("hSystematicErrorsASYieldMerged_PtBinAssoc%d", iBinPtHad + 1))); + hASSigmaSyst = reinterpret_cast(inFileFitSystematicErrors->Get(Form("hSystematicErrorsASSigmaMerged_PtBinAssoc%d", iBinPtHad + 1))); + + for (int iBinPtCand = 0; iBinPtCand < nBinsPtCand; iBinPtCand++) { + hBaselinSyst->SetBinError(iBinPtCand + 1, hBaselinSyst->GetBinContent(iBinPtCand + 1) * hBaselin[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hNSYieldSyst->SetBinError(iBinPtCand + 1, hNSYieldSyst->GetBinContent(iBinPtCand + 1) * hNSYield[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hNSSigmaSyst->SetBinError(iBinPtCand + 1, hNSSigmaSyst->GetBinContent(iBinPtCand + 1) * hNSSigma[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hASYieldSyst->SetBinError(iBinPtCand + 1, hASYieldSyst->GetBinContent(iBinPtCand + 1) * hASYield[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hASSigmaSyst->SetBinError(iBinPtCand + 1, hASSigmaSyst->GetBinContent(iBinPtCand + 1) * hASSigma[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + + hBaselinSyst->SetBinContent(iBinPtCand + 1, hBaselin[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hNSYieldSyst->SetBinContent(iBinPtCand + 1, hNSYield[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hNSSigmaSyst->SetBinContent(iBinPtCand + 1, hNSSigma[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hASYieldSyst->SetBinContent(iBinPtCand + 1, hASYield[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + hASSigmaSyst->SetBinContent(iBinPtCand + 1, hASSigma[iBinPtHad]->GetBinContent(iBinPtCand + 1)); + } + + c1->cd(); + SetTH1HistoStyle(hNSYieldSyst, Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "Y^{NS}", kFullSquare, kRed, 1.8, kRed, 2); + hNSYieldSyst->SetFillStyle(0); + hNSYieldSyst->Draw("E2same"); + + c2->cd(); + SetTH1HistoStyle(hASYieldSyst, Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "Y^{AS}", kFullSquare, kMagenta, 1.8, kMagenta, 2); + hASYieldSyst->SetFillStyle(0); + hASYieldSyst->Draw("E2same"); + + c3->cd(); + SetTH1HistoStyle(hNSSigmaSyst, Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "#sigma_{NS}", kFullSquare, kOrange + 8, 1.8, kOrange + 8, 2); + hNSSigmaSyst->SetFillStyle(0); + hNSSigmaSyst->Draw("E2same"); + + c4->cd(); + SetTH1HistoStyle(hASSigmaSyst, Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "#sigma_{AS}", kFullSquare, kViolet - 5, 1.8, kViolet - 5, 2); + hASSigmaSyst->SetFillStyle(0); + hASSigmaSyst->Draw("E2same"); + + c5->cd(); + SetTH1HistoStyle(hBaselinSyst, Form("p_{T}^{assoc} > %.1f GeV/c", binsPtHadIntervals[iBinPtHad]), "p_{T} (GeV/c)", "Baseline", kFullSquare, kBlue, 1.8, kBlue, 2); + hBaselinSyst->SetFillStyle(0); + hBaselinSyst->Draw("E2same"); + } + + c1->SaveAs(Form("Output_CorrelationFitting_%s_png/NearSideYield_PtAssoc%d.png", CodeNameAnalysis.data(), iBinPtHad + 1)); + c2->SaveAs(Form("Output_CorrelationFitting_%s_png/AwaySideYield_PtAssoc%d.png", CodeNameAnalysis.data(), iBinPtHad + 1)); + c3->SaveAs(Form("Output_CorrelationFitting_%s_png/NearSideSigma_PtAssoc%d.png", CodeNameAnalysis.data(), iBinPtHad + 1)); + c4->SaveAs(Form("Output_CorrelationFitting_%s_png/AwaySideSigma_PtAssoc%d.png", CodeNameAnalysis.data(), iBinPtHad + 1)); + c5->SaveAs(Form("Output_CorrelationFitting_%s_png/Baseline_PtAssoc%d.png", CodeNameAnalysis.data(), iBinPtHad + 1)); + c1->SaveAs(Form("Output_CorrelationFitting_%s_Root/NearSideYield_PtBinAssoc%d.root", CodeNameAnalysis.data(), iBinPtHad + 1)); + c2->SaveAs(Form("Output_CorrelationFitting_%s_Root/AwaySideYield_PtBinAssoc%d.root", CodeNameAnalysis.data(), iBinPtHad + 1)); + c3->SaveAs(Form("Output_CorrelationFitting_%s_Root/NearSideSigma_PtBinAssoc%d.root", CodeNameAnalysis.data(), iBinPtHad + 1)); + c4->SaveAs(Form("Output_CorrelationFitting_%s_Root/AwaySideSigma_PtBinAssoc%d.root", CodeNameAnalysis.data(), iBinPtHad + 1)); + c5->SaveAs(Form("Output_CorrelationFitting_%s_Root/Baseline_PtBinAssoc%d.root", CodeNameAnalysis.data(), iBinPtHad + 1)); } return; @@ -174,9 +459,36 @@ void FitCorrel(TString cfgFileName = "config_CorrAnalysis.json") void SetTH1HistoStyle(TH1D*& histo, TString hTitle, TString hXaxisTitle, TString hYaxisTitle, Style_t markerStyle, Color_t markerColor, Double_t markerSize, - Color_t lineColor, Int_t lineWidth, Float_t hTitleXaxisOffset, Float_t hTitleYaxisOffset, - Float_t hTitleXaxisSize, Float_t hTitleYaxisSize, Float_t hLabelXaxisSize, Float_t hLabelYaxisSize, - Bool_t centerXaxisTitle, Bool_t centerYaxisTitle) + Color_t lineColor, Int_t lineWidth, Float_t hTitleXaxisOffset = 1.3, Float_t hTitleYaxisOffset = 1.3, + Float_t hTitleXaxisSize = 0.045, Float_t hTitleYaxisSize = 0.045, Float_t hLabelXaxisSize = 0.045, Float_t hLabelYaxisSize = 0.045, + Bool_t centerXaxisTitle = false, Bool_t centerYaxisTitle = false) +{ + + histo->SetTitle(hTitle.Data()); + histo->GetXaxis()->SetTitle(hXaxisTitle.Data()); + histo->GetYaxis()->SetTitle(hYaxisTitle.Data()); + histo->SetMarkerStyle(markerStyle); + histo->SetMarkerColor(markerColor); + histo->SetMarkerSize(markerSize); + histo->SetLineColor(lineColor); + histo->SetLineWidth(lineWidth); + histo->GetXaxis()->SetTitleOffset(hTitleXaxisOffset); + histo->GetYaxis()->SetTitleOffset(hTitleYaxisOffset); + histo->GetXaxis()->SetTitleSize(hTitleXaxisSize); + histo->GetYaxis()->SetTitleSize(hTitleYaxisSize); + histo->GetXaxis()->SetLabelSize(hLabelXaxisSize); + histo->GetYaxis()->SetLabelSize(hLabelYaxisSize); + histo->GetXaxis()->CenterTitle(centerXaxisTitle); + histo->GetYaxis()->CenterTitle(centerYaxisTitle); + + return; +} + +void SetTH1HistoStyle(TH1F*& histo, TString hTitle, TString hXaxisTitle, TString hYaxisTitle, + Style_t markerStyle, Color_t markerColor, Double_t markerSize, + Color_t lineColor, Int_t lineWidth, Float_t hTitleXaxisOffset = 1.3, Float_t hTitleYaxisOffset = 1.3, + Float_t hTitleXaxisSize = 0.045, Float_t hTitleYaxisSize = 0.045, Float_t hLabelXaxisSize = 0.045, Float_t hLabelYaxisSize = 0.045, + Bool_t centerXaxisTitle = false, Bool_t centerYaxisTitle = false) { histo->SetTitle(hTitle.Data()); diff --git a/PWGHF/HFC/Macros/config_CorrAnalysis_DsToKKPi.json b/PWGHF/HFC/Macros/config_CorrAnalysis_DsToKKPi.json index 27e71b97514..5866d619f87 100644 --- a/PWGHF/HFC/Macros/config_CorrAnalysis_DsToKKPi.json +++ b/PWGHF/HFC/Macros/config_CorrAnalysis_DsToKKPi.json @@ -1,13 +1,23 @@ { - "CodeName": "Default", - "_CodeName": "Name to identify the analysis", - "pathFileSE": "~/cernbox/DsCorrelationAnalysis/MasterDegreeAnalysis/AnalysisResultsFinal.root", - "pathFileME": "~/cernbox/DsCorrelationAnalysis/MasterDegreeAnalysis/AnalysisResultsFinalME.root", - "pathFileMass": "/data/dataalice/scattar/RawYieldsDs.root", + "CodeName": "Test_example", + "_CodeName": "Name to identify the analysis: pass7_Full_ptAssoc0p3, Thin2023_Full_ptAssoc0p3", + "pathFileSE": "~/cernbox/AnalysisResults_SE_pp_Thin2023_FullAnalysis.root", + "_pathFileSE": "Name SE file", + "pathFileME": "~/cernbox/AnalysisResults_ME_derived_Thinned2023.root", + "_pathFileME": "Name ME file", + "pathFileMass": "~/cernbox/DsCorrelationAnalysis/ClassesForAnalysis/InvMassFitter/InvMass_pp_Thin2023_Final_Default.root", + "_pathFileMass": "Name mass file", + "pathFileFDTemplate": "~/cernbox/TemplatesCRMode2_ptAssocInt.root", + "pathFileFDPromptFrac": "~/cernbox/CutVarDs_FDVarDefaultBkg.root", + "pathFileSecPart": "~/cernbox/AnalysisResults_MC_Anchored_Thin2023_SecPart.root", + "_pathFileSecPart": "Path file for secondary particle contamination correction", "InputDirSE": "hf-task-correlation-ds-hadrons", "InputDirME": "hf-task-correlation-ds-hadrons", + "InputDirSecPart": "hf-task-correlation-ds-hadrons", "InputHistoCorrSignalName": "hCorrel2DVsPtSignalRegion", "InputHistoCorrSidebaName": "hCorrel2DVsPtSidebands", + "InputHistoCorrSidebaLeftName": "hCorrel2DVsPtSidebandLeft", + "InputHistoCorrSidebaRightName": "hCorrel2DVsPtSidebandRight", "InputHistoMassName": [ "hRawYieldsSignal", "hRawYieldsBkg", @@ -18,21 +28,34 @@ "Bkg yield vs pt", "SBs yield vs pt" ], + "InputHistoFDTemplatePrompt": "EtaPhiPromptNorm_pt", + "InputHistoFDTemplateNonPrompt": "EtaPhiNonPromptNorm_pt", + "InputHistoFDPromptFrac": "hRawFracPrompt", + "InputHistoPrimaryPart": "hCorrel2DVsPtPhysicalPrimaryMcRec", + "InputHistoAllPart": "hCorrel2DVsPtSignalRegionMcRec", + "IsDividedSideb": true, + "UseSidebLeft": true, + "UseSidebRight": false, + "ApplySecPartCorr": true, + "ApplyFDCorr": true, + "ApplyBiasBtoDCorr": true, + "pathfFilePromptMcRec": "~/cernbox/DsCorrelationAnalysis/ClassesForAnalysis/MonteCarloClosureTest/MonteCarloClosureTest_Prompt_All.root", + "pathfFileNonPromptMcRec": "~/cernbox/DsCorrelationAnalysis/ClassesForAnalysis/MonteCarloClosureTest/MonteCarloClosureTest_NonPrompt_All.root", "binsPtCandIntervals": [ - 2.0, + 1.0, 3.0, - 4.0, 5.0, 8.0, - 16.0 + 16.0, + 36.0 ], "binsPtHadIntervals": [ - 0.0, - 11.0 + 0.3, + 50.0 ], "deltaEtaInterval": [ - -0.8, - 0.8 + -1.0, + 1.0 ], "DmesonSpecie": 2, "_DmesonSpecie": [ @@ -45,7 +68,36 @@ "_NumberOfPools": "Number of pools for the event mixing", "CorrectPoolsSeparately": false, "_CorrectPoolsSeparately": "true = pool-by-pool ME correction; false = merged-pools ME correction", - "FitFunction": 1, - "FixBaseline": 0, - "FixMean": 3 + "RebinAngCorr": false, + "RebinFDCorr": true, + "RebinSecPart": false, + "nRebinDeltaPhi": 2, + "nRebinDeltaEta": 2, + "InputFileNameFitCorr": "ExtractCorrelationsResults_Reflected.root", + "_InputFileNameFitCorr": "ExtractCorrelationsResults_Reflected.root - ExtractCorrelationsResults.root", + "IsRiflected": true, + "DrawSystematics": false, + "SameSystematics": false, + "_SameSystematics": "true: systematic take from the first delta phi bin, false: systematic applied bin per bin", + "FitFunction": [ + 2, + 2, + 2, + 2, + 2 + ], + "_FitFunction": "1 : const + G NS + G AS (w/o periodicity), 2 : const + G NS + G AS (w/ periodicity), 3 : const + G AS, 4 : const + GenG NS + G AS (w/ periodicity), 5 : const + VonMises NS + VonMises AS (w/periodicity), 6 : const + VonMises AS", + "FixBaseline": 4, + "_FixBaseline": "0 : baseline free, = 1 : fix the baseline to the minimum of the histogram, < 0 : fix the baseline to the weighted average of the abs(fFixBaseline) lower points, = 3 : fix the baseline to the weighted average of the points passed through the function SetPointsForBaseline(), 4 : fix the baseline to the weighted average of default transverse region points", + "FixMean": 3, + "_FixMean": "= 0 : NS & AS mean free, = 1 : NS mean fixed to 0, AS mean free, = 2 : AS mean fixed to pi, NS mean free, = 3 : NS mean fixed to 0, AS mean to pi", + "ShiftBaseUp": false, + "ShiftBaseDown": false, + "nBaselinePoints": 4, + "binsForBaseline": [ + 6, + 7, + 8, + 9 + ] }