Get unfolded spectrum with improved stat errors

Jets/JetSpectrum_DrawingFunctions.C

@@ -23,6 +23,10 @@
 #include "RooUnfoldSvd.h"
 #include "TSVDUnfold.h"
 
+#include <TMultiGraph.h>
+#include <TGraphErrors.h>
+#include <TAxis.h>
+
 //My Libraries
 #include "./JetSpectrum_settings.h"
 #include "./JetSpectrum_inputs.h"
@@ -1112,6 +1116,8 @@ void Draw_Pt_spectrum_unfolded_singleDataset(int iDataset, int iRadius, int unfo
     }
   }
 }
+// declare this function here. It is used in the funciton below "Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset"
+void DrawRatioWithOffset(TH1D* histList[], int nUnfoldIteration,   const TString& yAxisTitle,const TString& canvasName, int unfoldIterationMax, int step, double yMin, double yMax);
 
 void Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset(int iDataset, int iRadius, int unfoldIterationMin, int unfoldIterationMax, int step, std::string options) {
 
@@ -1268,6 +1274,14 @@ void Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset(int iDataset, in
     TString* pdfName_ratio_refoldedComp_zoom = new TString("jet_"+jetType[iJetType]+"_"+jetLevel[iJetLevel]+"_"+partialUniqueSpecifier+"_Pt_unfolded_"+unfoldingInfo+"_ratioRefoldedUnfolded_zoom");
     Draw_TH1_Histograms(H1D_jetPt_ratio_measuredRefolded, unfoldingIterationLegend, nUnfoldIteration, textContext, pdfName_ratio_refoldedComp_zoom, texPtX, texRatioRefoldedMeasured, texCollisionDataInfo, drawnWindowUnfoldedMeasurement, legendPlacementAuto, contextPlacementAuto, "zoomToOneLarge,ratioLine,zoomToOneMedium2");
   }
+
+    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////ratio refolded/measured with offset differetn k////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    DrawRatioWithOffset(H1D_jetPt_ratio_measuredRefolded, nUnfoldIteration, "Refolded / Measured", "ratio_RefoldMeasure_different_k_withOffset", unfoldIterationMax, step, 0.8, 1.3);
+    DrawRatioWithOffset(H1D_jetPt_ratio_mcp, nUnfoldIteration, "unfolded / mcp", "ratio_Unfolded_mcp_different_k_withOffset", unfoldIterationMax, step, 0.7, 1.4); // last two numbers are y min and max
+    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
 }
 
 void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParameterInput, std::string options) {
@@ -1569,6 +1583,109 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete
   // }
 }
 
+void Draw_Pt_spectrum_unfolded_ImprovedStatErrors(int iDataset, int iRadius, int unfoldParameterInput, std::string options) {
+
+TH1D* measuredInput;
+if (!normGenAndMeasByNEvtsForUnfoldingInput) {
+  Get_Pt_spectrum_bkgCorrected_recBinning_preWidthScalingAtEndAndEvtNorm(measuredInput, iDataset, iRadius, options); 
+  if (useFineBinningTest) {
+    Get_Pt_spectrum_bkgCorrected_fineBinning_preWidthScalingAtEndAndEvtNorm(measuredInput, iDataset, iRadius, options);
+  }
+} else{
+  Get_Pt_spectrum_bkgCorrected_recBinning_preWidthScalingAtEnd(measuredInput, iDataset, iRadius, options);
+  if (useFineBinningTest) {
+    Get_Pt_spectrum_bkgCorrected_fineBinning_preWidthScalingAtEnd(measuredInput, iDataset, iRadius, options);
+  }
+}
+TH1D* H1D_jetPt_unfolded_withImprovedErrors;
+Get_Pt_spectrum_unfolded_ImprovedStatisticalErrors(H1D_jetPt_unfolded_withImprovedErrors, measuredInput, iDataset, iRadius, unfoldParameterInput, options);
+
+TString pdfname = "jet_Pt_spectrum_unfolded_RelativeUncertainty_Dataset"+DatasetsNames[iDataset]+"_R"+Form("%.1f", arrayRadius[iRadius])+"_k"+Form("%i", unfoldParameterInput);
+// TString textContext = "Unfolded improved errors";
+TString textContext(contextCustomOneField(*texDatasetsComparisonCommonDenominator, ""));
+// Test : 
+// TCanvas* c_relUnc = new TCanvas(pdfname, pdfname, 800, 800);
+// H1D_jetPt_unfolded_withImprovedErrors->Draw();
+// c_relUnc->SetLogy();
+
+// error with Draw_TH1_Histogram!!!?
+// Draw_TH1_Histogram(H1D_jetPt_unfolded_withImprovedErrors, textContext, pdfname, texPtJetRec, texJet_d2Ndptdeta_EventNorm, texCollisionDataInfo, drawnWindowUnfoldedMeasurement, legendPlacementAuto, contextPlacementAuto, "logy");
+
+}
+
+void DrawRatioWithOffset(TH1D* histList[], int nUnfoldIteration, const TString& yAxisTitle,const TString& canvasName, int unfoldIterationMax, int step, double yMin, double yMax){
+    // DrawRatioWithOffset(..., -1, -1);
+    TString canvasNameFull = canvasName + "_" + unfoldingMethod;
+    TCanvas* c = new TCanvas(canvasNameFull, canvasNameFull, 800, 800);
+    TMultiGraph* mg = new TMultiGraph();
+
+    int customColors[] = {
+        kBlack, kRed + 1, kBlue + 1, kGreen + 2, kOrange + 7, kViolet + 1
+    };
+    int nColors = sizeof(customColors) / sizeof(customColors[0]);
+
+    //int nUnfoldIteration = histList.size();  // infer number of iterations from input list
+
+    for (int i = 0; i < nUnfoldIteration; ++i) {
+        TH1D* h = histList[i];
+        int nBins = h->GetNbinsX();
+
+        std::vector<double> x_vals, y_vals, ex_vals, ey_vals;
+
+        for (int bin = 1; bin <= nBins; ++bin) {
+            // Replace ptBinsJetsRec[iRadius] with your own binning logic if needed:
+            double binLowEdge = h->GetXaxis()->GetBinLowEdge(bin);
+            double binUpEdge  = h->GetXaxis()->GetBinUpEdge(bin);
+            double binWidth   = binUpEdge - binLowEdge;
+            double offset     = (i - nUnfoldIteration / 2.0) * 0.065 * binWidth;
+
+            double x  = h->GetBinCenter(bin) + offset;
+            double y  = h->GetBinContent(bin);
+            double ex = 0;
+            double ey = h->GetBinError(bin);
+
+            x_vals.push_back(x);
+            y_vals.push_back(y);
+            ex_vals.push_back(ex);
+            ey_vals.push_back(ey);
+        }
+
+        TGraphErrors* gr = new TGraphErrors(nBins, &x_vals[0], &y_vals[0], &ex_vals[0], &ey_vals[0]);
+        int color = customColors[i % nColors];
+        gr->SetMarkerStyle(20 + i);
+        gr->SetMarkerColor(color);
+        gr->SetLineColor(color);
+        gr->SetMarkerSize(0.8);
+        gr->SetTitle(Form("k_{unfold} = %d", unfoldIterationMax-step*i));  // Adjust if you have different logic
+
+        mg->Add(gr, "P");
+    }
+
+    mg->Draw("A");
+    mg->GetXaxis()->SetLimits(5.0, 140.0);
+    if (yMin < yMax) {
+    mg->GetYaxis()->SetRangeUser(yMin, yMax);
+    }
+    mg->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+    mg->GetYaxis()->SetTitle(yAxisTitle);
+
+    c->BuildLegend();
+    c->Update();
+
+    // Draw horizontal reference line at y=1
+    double xmin = mg->GetXaxis()->GetXmin();
+    double xmax = mg->GetXaxis()->GetXmax();
+    TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
+    line->SetLineStyle(2);
+    line->SetLineColor(kGray + 2);
+    line->Draw("same");
+
+    c->Update();
+    // Auto-save
+    c->SaveAs(canvasNameFull + ".pdf");
+    c->SaveAs(canvasNameFull + ".png");
+}
+
 // rename refoldedUnfolded as closure test?
 // and try and spend 15 min to clean hist names for the spectrum analysis
 
@@ -1578,4 +1695,6 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete
 // hMcEfficiency_vsPt->Divide(hMcSignalCount_vsPt,TrueV0PtSpectrum_AnalysisBins, 1., 1., "b"); // option b for binomial because efficiency: https://twiki.cern.ch/twiki/bin/view/ALICE/PWGLFPAGSTRANGENESSEfficiency
 
 
+
+
 #endif

Jets/JetSpectrum_DrawingFunctions.h

@@ -16,6 +16,7 @@ void Draw_Pt_spectrum_unfolded_singleDataset(int iDataset, int iRadius, int unfo
 void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParameterInput, std::string options);
 void Draw_Pt_TestSpectrum_unfolded(int iDataset, int iRadius, std::string options);
 void Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset(int iDataset, int iRadius, int unfoldIterationMin, int unfoldIterationMax, int step, std::string options);
+void Draw_Pt_spectrum_unfolded_ImprovedStatErrors(int iDataset, int iRadius, int unfoldParameterInput, std::string options);
 
 void Draw_Pt_efficiency_jets(int iRadius, std::string options);
 void Draw_kinematicEfficiency(int iRadius, std::string options);

Jets/JetSpectrum_Unfolding.C

@@ -19,6 +19,9 @@
 #include "../Utilities/HistogramUtilities.C"
 #include "../Utilities/HistogramPlotting.C" 
 
+#include "TRandom3.h"
+#include "TProfile.h"
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////// RooUnfold Custom Utilities ///////////////////////////////////////////////////////////////////////////////
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -861,4 +864,80 @@ void Get_Pt_spectrum_mcpFoldedWithFluctuations(TH1D* &H1D_jetPt_mcpFolded, int i
   if (showFunctionInAndOutLog) {cout << "--- OUT Get_Pt_spectrum_mcpFoldedWithFluctuations()" << endl;};
 }
 
+
+void Get_Pt_spectrum_unfolded_ImprovedStatisticalErrors(TH1D* &H1D_jetPt_unfolded, TH1D* &measuredInput, int iDataset, int iRadius, int unfoldParameterInput, std::string options){
+  cout << "--- IN  Get_Pt_spectrum_unfolded_ImprovedStatisticalError()" << endl;
+
+  TString partialUniqueSpecifier = Datasets[iDataset]+DatasetsNames[iDataset]+Form("%.1d",iDataset)+"_R="+Form("%.1f",arrayRadius[iRadius]);
+  TH1D* measured = (TH1D*)measuredInput->Clone("measured_Get_Pt_spectrum_unfolded_preWidthScalingAtEndAndEvtNorm"+partialUniqueSpecifier); // before smearing
+
+  TH1D* TH1D_UnfoldNominal;
+  int unfoldParameterNominal;
+  unfoldParameterNominal = Get_Pt_spectrum_unfolded(TH1D_UnfoldNominal, measured, iDataset, iRadius, unfoldParameterInput, options).first;
+
+  const int nToys = 5;
+  TH1D* smearedMeasured[nToys];  
+  TH1D* TH1D_unfoldSmeared[nToys];
+  int unfoldParameter[nToys];
+
+  TProfile* TProfile_JetPtVar = new TProfile("TProfile_JetPtVar","",nBinPtJetsGen[iRadius], ptBinsJetsGen[iRadius],"S");   
+  for (int itoy = 0; itoy < nToys; ++itoy) {
+    smearedMeasured[itoy] = (TH1D*)measured->Clone(Form("H1D_jetPt_smeared_measured_%d", itoy));
+    smearedMeasured[itoy]->Reset("ICE");
+    for (int i = 1; i <= measured->GetNbinsX(); ++i) {
+        double M_i     = measured->GetBinContent(i);
+        double sigma_i = measured->GetBinError(i);
+        double M_i_smeared = gRandom->Gaus(M_i, sigma_i);
+        smearedMeasured[itoy]->SetBinContent(i, M_i_smeared);
+        smearedMeasured[itoy]->SetBinError(i, sigma_i);
+    }
+    unfoldParameter[itoy] = Get_Pt_spectrum_unfolded(TH1D_unfoldSmeared[itoy], smearedMeasured[itoy], iDataset, iRadius, unfoldParameterInput, options).first;
+    int nBins = TH1D_unfoldSmeared[itoy]->GetNbinsX();
+    for (int binx = 1; binx <= nBins; binx++) {
+        double cent = TH1D_unfoldSmeared[itoy]->GetXaxis()->GetBinCenter(binx);
+        double cont = TH1D_unfoldSmeared[itoy]->GetBinContent(binx);
+        TProfile_JetPtVar->Fill(cent, cont); 
+    }
+    cout << "############## SVD unfolding: toy " << itoy+1 << " out of " << nToys << ", Done! ##############" << endl;
+  }
+
+  TH1D* TH1D_RooUnfold_RelativeError = (TH1D*) TH1D_UnfoldNominal->Clone("H1D_Unfolded_Relative_Uncert_original_error"+partialUniqueSpecifier);
+  TH1D_RooUnfold_RelativeError->Reset();
+
+  TH1D* TH1D_UnfSmeared_RelativeError = (TH1D*) TH1D_UnfoldNominal->Clone("H1D_jetPt_Relative_Uncert_SmearedThenUnfolded"+partialUniqueSpecifier);
+  TH1D_UnfSmeared_RelativeError->Reset();
+
+  H1D_jetPt_unfolded = (TH1D*)TH1D_UnfoldNominal->Clone("H1D_jetPt_unfolded_improvedStatisticalErrors"+partialUniqueSpecifier);
+
+  int nBins = TH1D_UnfoldNominal->GetNbinsX();
+  for (int bin = 1; bin <= nBins; bin++) {
+      double val = TH1D_UnfoldNominal->GetBinContent(bin);
+      double err = TH1D_UnfoldNominal->GetBinError(bin);
+      TH1D_RooUnfold_RelativeError->SetBinContent(bin, val > 0 ? err / val : 0);
+      TH1D_RooUnfold_RelativeError->SetBinError(bin, 0);
+
+      double mean_smeared = TProfile_JetPtVar->GetBinContent(bin);   // ⟨x⟩
+      double rms_smeared  = TProfile_JetPtVar->GetBinError(bin);     // RMS ("S" option)
+      TH1D_UnfSmeared_RelativeError->SetBinContent(bin, mean_smeared > 0 ? rms_smeared / mean_smeared : 0);
+      TH1D_UnfSmeared_RelativeError->SetBinError(bin, 0);
+
+      H1D_jetPt_unfolded->SetBinError(bin, rms_smeared);
+  }
+
+  // TString* pdfName_RooUnfoldUncertainty = new TString("RooUnfold_relative_uncertainty_jetPt"+partialUniqueSpecifier);
+  // // TString textContext(contextCustomOneField(*texDatasetsComparisonCommonDenominator, ""));
+  // Draw_TH1_Histogram(TH1D_UnfSmeared_RelativeError, textContext, pdfName_RooUnfoldUncertainty, texPtJetRec, texRelativeErrPercent, texCollisionDataInfo, drawnWindowUnfoldedMeasurement, legendPlacementAuto, contextPlacementAuto, "logy");
+  TString* pdfName_realtiveError = new TString("relative_uncertainty_jetPt_defaultRooUnfold_and_SmearedMeasured"+partialUniqueSpecifier);
+  TString textContext = contextCustomTwoFields(*texDatasetsComparisonCommonDenominator, contextJetRadius(arrayRadius[iRadius]), "");
+  TString LegendRelativeErrors[2] = {"Default RooUnfold", "smeared then unfolded"};
+  TH1D** RelativeErrors = new TH1D*[2];
+  RelativeErrors[0] = TH1D_RooUnfold_RelativeError;
+  RelativeErrors[1] = TH1D_UnfSmeared_RelativeError;
+
+  Draw_TH1_Histograms(RelativeErrors, LegendRelativeErrors, 2, textContext, pdfName_realtiveError, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, "logy");
+}
+
+
+
+
 #endif

Jets/JetSpectrum_Unfolding.h

@@ -17,6 +17,8 @@ void Get_Pt_spectrum_dataUnfoldedThenRefolded_RooUnfoldMethod(TH1D* &H1D_jetPt_u
 void Get_Pt_spectrum_mcpFoldedWithFluctuations_preWidthScalingAtEndAndEvtNorm(TH1D* &H1D_jetPt_mcpFolded, int iDataset, int iRadius, std::string options);
 void Get_Pt_spectrum_mcpFoldedWithFluctuations_preWidthScalingAtEnd(TH1D* &H1D_jetPt_mcpFolded, int iDataset, int iRadius, std::string options);
 void Get_Pt_spectrum_mcpFoldedWithFluctuations(TH1D* &H1D_jetPt_mcpFolded, int iDataset, int iRadius, std::string options);
+void Get_Pt_spectrum_unfolded_ImprovedStatisticalErrors(TH1D* &H1D_jetPt_unfolded, TH1D* &measuredInput, int iDataset, int iRadius, int unfoldParameterInput, std::string options);
+
 
 
 #endif