Get unfolded spectrum with improved stat errors

Jets/JetQC_inputs_template.h

@@ -1965,15 +1965,15 @@ TFile* file_AliAnalysis;
 
 
 //////// -------- LHC25b6 - pp sim anchored to PbPb localTest ////////
-TString* texCollisionDataInfo = new TString("pp #sqrt{#it{s}} = 13.6 TeV");
+TString* texCollisionDataInfo = new TString("pp #sqrt{#it{s}} = 5.36 TeV");
 const TString* texDatasetsComparisonType = new TString("");
 const TString* texDatasetsComparisonCommonDenominator = new TString("");
 const int nDatasets = 2;
-const TString Datasets[nDatasets] = {"pp13TeV_sim_MB", "pp13TeV_sim_jetjet"};
+const TString Datasets[nDatasets] = {"LHC26a6_train615296_Unf", "LHC25b4ab6_train600389"};
 // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"};
-const TString DatasetsNames[nDatasets] = {"MB", "jet-jet"};
-TFile* file_O2Analysis_list[nDatasets] = {new TFile("Datasets/"+Datasets[0]+"/AnalysisResults.root"),
-                                          new TFile("Datasets/"+Datasets[1]+"/AnalysisResults.root")
+const TString DatasetsNames[nDatasets] = {"JJ LHC26a6", "MB LHC25b4ab6"};
+TFile* file_O2Analysis_list[nDatasets] = {new TFile("../Datasets/"+Datasets[0]+"/AnalysisResults.root"),
+                                          new TFile("../Datasets/"+Datasets[1]+"/AnalysisResults.root")
                                         };
 
 const TString analysisWorkflow[nDatasets] = {"jet-spectra-charged", "jet-spectra-charged",
@@ -1982,6 +1982,6 @@ const TString analysisWorkflow[nDatasets] = {"jet-spectra-charged", "jet-spectra
 const TString wagonId[nDatasets] = {"", ""
                                     };
 
-const bool isDatasetWeighted[nDatasets] = {false, true};
+const bool isDatasetWeighted[nDatasets] = {true, false};
 const std::string histDatasetComparisonStructure = "";
 const bool datasetsAreSubsetsofId0 = false;

Jets/JetQC_settings_template.h

@@ -15,8 +15,8 @@ const TString jetLevel[nJetLevel] = {"data", "mcd", "mcp"};
 // float arrayRadius[nRadius] = {0.2, 0.4, 0.6};
 // const float areaDisplayMax[nRadius] = {0.5, 1, 1.5};
 const int nRadius = 1;
-const TString RadiusLegend[nRadius] = {"R = 0.2"};
-float arrayRadius[nRadius] = {0.2};
+const TString RadiusLegend[nRadius] = {"R = 0.4"};
+float arrayRadius[nRadius] = {0.4};
 const float areaDisplayMax[nRadius] = {0.5};
 // const int nRadius = 9;
 // const TString RadiusLegend[nRadius] = {"R = 0.2", "R = 0.25", "R = 0.3", "R = 0.35", "R = 0.4", "R = 0.45", "R = 0.5", "R = 0.55", "R = 0.6"};
@@ -28,13 +28,13 @@ TString methodRandomConeHistLegend[nMethodRC] = {"Random Cones (RC)", "RC w/o le
 
 // Choice of jet type (charged, neutral, full) and level (data, detector level, particle level)
 const int iJetType = 0;
-const int iJetLevel = 0;
+const int iJetLevel = 1;
 
 // Choice of jet QA type (uncorrected jets, background corrected jet (rho area version), event wise constituent subtraction
-const int iJetFinderQaType = 1;
+const int iJetFinderQaType = 0;
 
 // Choice of Random Cone method:
-const int iMethodRandomCone = 1; 
+const int iMethodRandomCone = 0; 
 // Default window for random cone:
 std::array<std::array<float, 2>, 2> drawnWindowRCdefault = {{{-30, 60}, {5E-7, 20}}}; // {{xmin, xmax}, {ymin, ymax}}
 

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"
@@ -695,7 +699,16 @@ void Draw_Pt_spectrum_unfolded_singleDataset(int iDataset, int iRadius, int unfo
 
   unfoldParameter = Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded, measuredInput, iDataset, iRadius, unfoldParameterInput, options).first;
   // TH1D* H1D_jetPt_unfolded2 = (TH1D*)H1D_jetPt_unfolded->Clone(H1D_jetPt_unfolded->GetName()+(TString)"H1D_jetPt_unfolded2");
-
+  ////////////////////////////////////////////////////////////
+
+  if (writeOutputRootFile) {
+    cout << "######################### FILE CREATED IN PRINCIPLE##################" << endl; 
+    TFile* outFile = new TFile("output.root", "UPDATE");   // "UPDATE" Open existing or create if missing / "RECREATE" Always deletes file and creates new one
+    H1D_jetPt_unfolded->Write("MB_LHC25b4ab6_train600389_ppref_unfolded"); //  MB_LHC25b4b5_train533385_Unf_ppref_unfolded
+    outFile->Close();
+    cout << "######################### HISTO SAVED IN PRINCIPLE ##################" << endl; 
+  }
+  ////////////////////////////////////////////////////////////
   cout << "comparison with raw measured" << endl; 
   if (!useFineBinningTest) {
     Get_Pt_spectrum_bkgCorrected_genBinning(H1D_jetPt_measured_genBinning, iDataset, iRadius, options);
@@ -1112,6 +1125,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 +1283,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) {
@@ -1379,8 +1402,10 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete
     divideSuccessDatasets[iDataset] = H1D_jetPt_unfolded_ratio_datasets[iDataset]->Divide(H1D_jetPt_unfolded[0]);
     // Creating to-be-plotted histograms
     datasetNameSpecifier[iDataset] = "_"+DatasetsNames[iDataset]+Form("%.1d",iDataset);
+    // ####################
+
 
-
+    // ####################
     cout << "comparison with raw measured" << endl; 
     if (!useFineBinningTest) {
       Get_Pt_spectrum_bkgCorrected_genBinning(H1D_jetPt_measured_genBinning[iDataset], iDataset, iRadius, options);
@@ -1569,6 +1594,153 @@ 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");
+}
+
+void MakeRatio()
+{
+    // Open file in UPDATE mode (so we can write result)
+    TFile* MB = TFile::Open("../20260226_Unf600389_LHC25b4ab6/output.root", "READ");
+    TFile* JJ = TFile::Open("output.root", "READ");
+
+    if (!MB || MB->IsZombie()) {
+        std::cerr << "Error: cannot open MB file!" << std::endl;
+        return;
+    }
+    if (!JJ || JJ->IsZombie()) {
+        std::cerr << "Error: cannot open JJ file!" << std::endl;
+        return;
+    }
+    // JJ_LHC26a6_train615296_Unf_ppref_unfolded //  MB_LHC25b4b5_train533385_Unf_ppref_unfolded
+    // Retrieve histograms
+    TH1D* h1 = (TH1D*)MB->Get("MB_LHC25b4ab6_train600389_ppref_unfolded");
+    TH1D* h2 = (TH1D*)JJ->Get("JJ_LHC26a6_train615296_Unf_ppref_unfolded");
+
+    if (!h1 || !h2) {
+        std::cerr << "Error: one of the histograms not found!" << std::endl;
+        MB->Close();
+        JJ->Close();
+        return;
+    }
+
+    // Clone first histogram to store ratio
+    TH1D* hRatio = (TH1D*)h1->Clone("hRatio");
+
+    // Perform division
+    hRatio->Divide(h2);
+    TString* pdfName = new TString("Ratio of MB to JJ");
+    TString textContext(contextCustomOneField(*texDatasetsComparisonCommonDenominator, "Ratio_of_DataUnfold_w_MB_LHC25b4ab6_to_JJ_LHC26a6"));
+    TString* YLabel = new TString("Unf MB / Unf JJ");
+
+    // Write ratio to file (overwrite if exists)
+    // hRatio->Write("", TObject::kOverwrite);
+    Draw_TH1_Histogram(hRatio, textContext, pdfName, texPtJetRec, YLabel, texCollisionDataInfo, drawnWindowUnfoldedMeasurement, legendPlacementAuto, contextPlacementAuto, "ratioLine");
+    // file->Close();
+
+    std::cout << "Ratio successfully created and saved." << std::endl;
+}
+
+
 // rename refoldedUnfolded as closure test?
 // and try and spend 15 min to clean hist names for the spectrum analysis
 
@@ -1578,4 +1750,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,9 +16,11 @@ 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 MakeRatio();
 void Draw_Pt_efficiency_jets(int iRadius, std::string options);
 void Draw_kinematicEfficiency(int iRadius, std::string options);
 void Draw_FakeRatio(int iRadius, std::string options);
 
+
 #endif