From 20e8ac3ac6c2adccfaa7f74763653dfb52707c1e Mon Sep 17 00:00:00 2001 From: ahmadtabikh Date: Mon, 26 Jan 2026 17:20:17 +0100 Subject: [PATCH 1/4] fixe some settings and add systematics computing functions --- Jets/JetSpectrum_DrawingFunctions.C | 89 +++++ Jets/JetSpectrum_RunMacro_template.C | 66 ---- Jets/JetSpectrum_inputs_template.h | 54 --- Jets/JetSpectrum_settings_template.h | 470 --------------------------- Jets/JetSpectrum_systematics.C | 230 ++++++++++++- 5 files changed, 305 insertions(+), 604 deletions(-) delete mode 100644 Jets/JetSpectrum_RunMacro_template.C delete mode 100644 Jets/JetSpectrum_inputs_template.h delete mode 100644 Jets/JetSpectrum_settings_template.h diff --git a/Jets/JetSpectrum_DrawingFunctions.C b/Jets/JetSpectrum_DrawingFunctions.C index f7cb231..f130429 100644 --- a/Jets/JetSpectrum_DrawingFunctions.C +++ b/Jets/JetSpectrum_DrawingFunctions.C @@ -23,6 +23,10 @@ #include "RooUnfoldSvd.h" #include "TSVDUnfold.h" +#include +#include +#include + //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,79 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete // } } +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 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 +1665,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 \ No newline at end of file diff --git a/Jets/JetSpectrum_RunMacro_template.C b/Jets/JetSpectrum_RunMacro_template.C deleted file mode 100644 index 4a784e0..0000000 --- a/Jets/JetSpectrum_RunMacro_template.C +++ /dev/null @@ -1,66 +0,0 @@ -// This is a template. To use the JetSpectrum drawing functions, rename this file to JetSpectrum_RunMacro.C and edit it how you want and run "root pathToFolderAnalysisToolsO2/Jets/JetSpectrum_RunMacro.C+" - -#include "./JetSpectrum_DrawingFunctions.h" -#include "./JetSpectrum_DrawingFunctions.C" - -using namespace std; - -// Misc utilities -void SetStyle(Bool_t graypalette=kFALSE); - -///////////////////////////////////////////////////// -///////////////////// Main Macro //////////////////// -///////////////////////////////////////////////////// - -void JetSpectrum_RunMacro() { - // Set the default style - SetStyle(); - - // gathers the analysis options in a single char[] - snprintf(optionsAnalysis, sizeof(optionsAnalysis), "%s,%s", unfoldingPrior, unfoldingMethod); - cout << "Analysis options are: " << optionsAnalysis << endl; - - mcCollHistIsObsolete = inputMcCollHistIsObsolete; - - int iDataset = 0; - int iRadius = 0; - - Draw_ResponseMatrices_Fluctuations(iDataset, iRadius); - Draw_ResponseMatrices_detectorResponse(iDataset, iRadius); - Draw_ResponseMatrices_DetectorAndFluctuationsCombined(iDataset, iRadius, optionsAnalysis); - - // Draw_ResponseMatrices_Fluctuations(1, iRadius); - // Draw_ResponseMatrices_detectorResponse(1, iRadius); - // Draw_ResponseMatrices_DetectorAndFluctuationsCombined(1, iRadius, optionsAnalysis); - - // // Draw_Pt_spectrum_unfolded_FluctResponseOnly(iDataset, iRadius, optionsAnalysis); // NOT FIXED YET - result meaningless - // Draw_Pt_spectrum_raw(iDataset, iRadius, optionsAnalysis); - // Draw_Pt_spectrum_raw(iDataset, iRadius, optionsAnalysis+(std::string)"noEventNormNorBinWidthScaling"); - // Draw_Pt_spectrum_mcp(iDataset, iRadius, optionsAnalysis); - // Draw_Pt_spectrum_mcp(iDataset, iRadius, optionsAnalysis+(std::string)"noEventNormNorBinWidthScaling"); - // Draw_Pt_spectrum_mcdMatched(iDataset, iRadius, optionsAnalysis); - // Draw_Pt_spectrum_mcdMatched(iDataset, iRadius, optionsAnalysis+(std::string)"noEventNormNorBinWidthScaling"); - - // Draw_Pt_efficiency_jets(iRadius, optionsAnalysis); - // Draw_kinematicEfficiency(iRadius, optionsAnalysis); - // Draw_FakeRatio(iRadius, optionsAnalysis); - - // int unfoldParameterInput = 5; - // Draw_Pt_spectrum_unfolded_singleDataset(iDataset, iRadius, unfoldParameterInput, optionsAnalysis); - // int unfoldParameterInput1 = 6; - // Draw_Pt_spectrum_unfolded_singleDataset(iDataset, iRadius, unfoldParameterInput1, optionsAnalysis); - // int unfoldParameterInput2 = 8; - // Draw_Pt_spectrum_unfolded_singleDataset(iDataset, iRadius, unfoldParameterInput2, optionsAnalysis); - // Draw_Pt_spectrum_unfolded_datasetComparison(iRadius, unfoldParameterInput2, optionsAnalysis); - int unfoldParameterInput3 = 10; - Draw_Pt_spectrum_unfolded_singleDataset(iDataset, iRadius, unfoldParameterInput3, optionsAnalysis); - // int unfoldParameterInput4 = 12; - // Draw_Pt_spectrum_unfolded_singleDataset(iDataset, iRadius, unfoldParameterInput4, optionsAnalysis); - - // int unfoldParameterInputMin = 7; - // int unfoldParameterInputMax = 12; - // int unfoldParameterInputStep = 2; - // Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset(iDataset, iRadius, unfoldParameterInputMin, unfoldParameterInputMax, unfoldParameterInputStep, optionsAnalysis); - -} - diff --git a/Jets/JetSpectrum_inputs_template.h b/Jets/JetSpectrum_inputs_template.h deleted file mode 100644 index 1cf0f61..0000000 --- a/Jets/JetSpectrum_inputs_template.h +++ /dev/null @@ -1,54 +0,0 @@ - -// This is a template. To use the JetSpectrum_DrawingMacro.C, rename this file to JetSpectrum_inputs.h and edit it how you want. - - -#ifndef JETSPECTRUM_INPUTS_H -#define JETSPECTRUM_INPUTS_H - -////////////////////////////////////////////////////////////////////////////////////////////////// -////////////////////////////// file access choice //////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////////////////////// -TFile* file_O2Analysis_run2ComparisonFileHannaBossiLaura = new TFile("Datasets/Run2_Unfolding_AreaBased_HannahMethod_R020_Nominal_ExtendedPtRange/Unfolding_AreaBased_HannahMethod_R020_Nominal_ExtendedPtRange.root"); -TFile* file_O2Analysis_run2ComparisonFileMLPaper = new TFile("Datasets/Run2_Unfolding_MachineLearningMethod_R020/Ch-jetSuppression_PbPb502TeV.root"); - - -//////// -------- LHC23zzh pass 4 with - pp sim anchored to PbPb 10% - lead05 /////// -TString* texEnergyPbPb = new TString("#sqrt{#it{s}_{NN}} = 5.36 TeV"); -TString* texEnergy = new TString("pp, #sqrt{#it{s}} = 5.36 TeV"); -TString* texCollisionDataType = new TString("0#font[122]{-}10% Pb#font[122]{-}Pb"); -TString* texCollisionDataInfo = new TString((TString)*texCollisionDataType+", "+(TString)*texEnergyPbPb); -TString* texCollisionMCType = new TString("PYTHIA + GEANT4"); -TString* texCollisionMCInfo = new TString((TString)*texCollisionMCType+", "+(TString)*texEnergy); -const TString* texDatasetsComparisonType = new TString("0#font[122]{-}10% cent."); -const TString* texDatasetsComparisonCommonDenominator = new TString("ALICE performance"); - -const int nDatasets = 1; -const TString Datasets[nDatasets] = {"LHC25b6_pp_sim_PbPbAnchor_train420439" - }; -const TString DatasetsNames[nDatasets] = {"ppAnchorPbPb"}; -TFile* file_O2Analysis_list[nDatasets] = {new TFile("Datasets/"+Datasets[0]+"/AnalysisResults.root") - }; -TFile* file_O2Analysis_MCfileForMatrix[nDatasets] = {new TFile("Datasets/LHC25b6_pp_sim_PbPbAnchor_train420439/AnalysisResults.root") - }; -TFile* file_O2Analysis_MCfile_UnfoldingControl_input[nDatasets] = {new TFile("Datasets/LHC25b6_pp_sim_PbPbAnchor_0010centEffCorrection_train571118_MCUnfoldingInput/AnalysisResults.root") - }; // use this MC file as input to unfolding (with h_jet_pt_rhoareasubtracted distrib on file) and as comparison to gen (with h_jet_pt_part distrib on file) - -TFile* file_O2Analysis_MCfile_UnfoldingControl_response[nDatasets] = {new TFile("Datasets/LHC25b6_pp_sim_PbPbAnchor_0010centEffCorrection_train571118_MCResponse/AnalysisResults.root") - }; // use this MC file as input to unfolding (with h_jet_pt_rhoareasubtracted distrib on file) and as comparison to gen (with h_jet_pt_part distrib on file) - - -const TString trainIdData = ""; -const TString analysisWorkflowData = "jet-spectra-charged_lead_05_100"+trainIdData; - -const TString trainIdBkg = ""; -const TString analysisWorkflowBkg = "jet-background-analysis"+trainIdBkg; -const TString trainIdUnfoldingControl = ""; -const TString analysisWorkflow_unfoldingControl = "jet-spectra-charged_lead_05_100"+trainIdUnfoldingControl; - -const TString trainIdMC = ""; -const TString analysisWorkflowMC = "jet-spectra-charged_lead_05_100"+trainIdMC; -const bool etaCutOnMatchedJetsIsObsoleteVersion = false; -bool inputMcCollHistIsObsolete = true; - - -#endif diff --git a/Jets/JetSpectrum_settings_template.h b/Jets/JetSpectrum_settings_template.h deleted file mode 100644 index 10cd54c..0000000 --- a/Jets/JetSpectrum_settings_template.h +++ /dev/null @@ -1,470 +0,0 @@ -// This is a template. To use the JetSpectrum_DrawingMacro.C, rename this file to JetSpectrum_settings.h and edit it how you want. - -#ifndef JETSPECTRUM_SETTINGS_H -#define JETSPECTRUM_SETTINGS_H - -// Analysis settings -const int nJetType = 3; -const TString jetType[nJetType] = {"charged", "neutral", "full"}; -const int nJetLevel = 3; -const TString jetLevel[nJetLevel] = {"data", "mcd", "mcp"}; -const int nRadius = 3; -const TString RadiusLegend[nRadius] = {"R = 0.2", "R = 0.4", "R = 0.6"}; -const float arrayRadius[nRadius] = {0.2, 0.4, 0.6}; -const TString arrayRadiusPdfName[nRadius] = {"02", "04", "06"}; -const int nRandomConeTypes = 5; -const TString randomConeTypeList[nRandomConeTypes] = {"", "withoutleadingjet", "randomtrackdirection", "randomtrackdirectionwithoutoneleadingjets", "randomtrackdirectionwithouttwoleadingjets"}; - -const double trackEtaRange[2] = {-0.9, 0.9}; -float deltaJetEta[3] = {1.4, 1, 0.6}; - -// Choice of jet type (charged, neutral, full) and level (data, detector level, particle level) -const int iJetType = 0; -const int iJetLevel = 0; -const int randomConeType = 1; - - -const float centralityRange[2] = {0, 10}; -// const float centralityRange[2] = {50, 70}; - -//////////////////////////////////////////////// -////////// Unfolding settings - start ////////// -//////////////////////////////////////////////// - - - -char unfoldingPrior[] = "mcpPriorUnfolding"; // prior options: mcpPriorUnfolding, mcdPriorUnfolding, measuredPriorUnfolding, noPriorUnfolding, testAliPhysics /////// if using mcp as prior, should have the leading track cut like data -const bool doYSliceNormToOneDetResp = true; // should be true in Pb-Pb analysis (done by marta) (because fluct. response matrix is probability distrib; in pp where we only have detector response, it's best to keep matrix of number of events as svd paper seems to advise: less weight on matrix entries with single count; if this is set to FALSE then noPriorUnfolding should be used, as det response already comes with weighting so it'd be done twice) -const bool doYSliceNormToOneCombinedResp = false; // should be false (not done by marta); breaks unfolding with svd if true -const bool doUnfoldingPriorDivision = false; // keep false for now ; unfolding doesn't work anymore if this is done, gives almost flat pT distribution, though refolding test is good; I am not sure why, might be because roounfold already does that; one good reason to avoid it; anyway, is that roounfold already seems to deal with errors; my error propagation doesn't take into account off-diagonal covariance elements, and so can only be worse -const bool scaleRespByXYWidth = false; // keep false, does not work well if true; although necessary (along with matrixTransformationOrder 1) to get a good looking response matrix with every entry below 1, unfolding with this struggles a bit though not horrible, it's still not great looking -const bool scaleRespByYWidth = false; // keep false -const bool matrixTransformationOrder = 0; // use 0 - //0: reweight with unfoldingPrior, then rebin with merging prior, then do YSliceNorm and scaleRespByXYWidth if set to true (works well 0); - //1: rebin, then YSliceNorm and scaleRespByXYWidth, then reweight; - //2: rebin, then reweight, then YSliceNorm and scaleRespByXYWidth - -char unfoldingMethod[] = "Svd"; // unfolding method options: Bayes, Svd -char optionsAnalysis[100] = ""; - -const bool doClosure_splitMC_mcdUnfoldedVsGen = true; -const bool doClosure_splitMC_mcpFoldedWithFluct = false; // if true, uses file_O2Analysis_MCfile_UnfoldingControl_input mcp distribution, folds it with the background fluctuation matrix, unfolds it with the merged det x bkg response matrix, and compares it to the mcp distribution in file_O2Analysis_MCfile_UnfoldingControl_input -// 13/11/2025: check this comparison; for some reason it's not a ratio of 1 when running on pp even though the bkg fluctuation matrix is identity -const bool doBkgSubtractionInData = false; -const bool doBkgSubtractionInMC = false; -const bool useFactorisedMatrix = false; // use factorised response matrix for unfolding, or not; if not, the fluctuations response it replaced by the identity matrix -const bool useFactorisedMatrixInMcdUnfoldedClosure = true; // whne doing SplitMcClosure: use factorised response matrix for unfolding, or not; if not, the fluctuations response it replaced by the identity matrix -const bool mcIsWeighted = true; // use if the MC has been weighted to have more high pt jets? -bool applyFakes = true; // only applied if useManualRespMatrixSettingMethod is true; 18/03: if false? -int applyEfficiencies = 2; // 2 is best; kinematic efficiency is already be handled by roounfold (02/04/2025; one can check simply with a pp unfolding with just det matrix and fine-ish binning like "// Joonsuk binning for pp with smaller rec window to test kinematic efficiency") -//applyEfficiencies: 0: no efficiency correction, 1: kine only, 2: jet finding efficiency only, 3: both active; only applied if useManualRespMatrixSettingMethod is true - -const bool doWidthScalingEarly = false; // to avoid pT bin width having an influence on spectrum; which one should be done? early or end? for now will be done at end -const bool doWidthScalingAtEnd = true; // - -const bool normDetRespByNEvts = false; //that's what breaks svd; https://arxiv.org/pdf/hep-ph/9509307 seems to say one should use the number of events matrix (see last paragraph of conclusion) instead of a probability matrix, to further reduce errors -const bool normGenAndMeasByNEvtsForUnfoldingInput = false; // controls normalisation of input to unfolding; if false the inputs aren't normalised; if true they are normalised IF normaliseDistribsInComparisonPlots is also true (probably should remove influence of this) -const bool normaliseUnfoldingResultsAtEnd = true; // if true: Njets per event ; if false: Njets total ; both normaliseUnfoldingResultsAtEnd and normaliseDistribsInComparisonPlots should be the same, else refolding test fails; -const bool normaliseDistribsInComparisonPlots = true; // if true: Njets per event ; if false: Njets total ; both normaliseUnfoldingResultsAtEnd and normaliseDistribsInComparisonPlots should be the same, else refolding test fails; without - -const bool useManualRespMatrixSettingMethod = true; // keep true; false uses Joonsuk's resp matrix setup with roounfold methods; if set to true, both refold methods are constistent; if set to false, the roounfold one is identical as if true, but the manual one becomes different and bad; EDIT 25/04/2025: Joonsuk's method I haven't kept updated, and now gives bad results on trivial refolding test -// joonsuk's seems to have better unfolded/mcp ratio, but worse unfolded/refolded ratio; d distribution is bad with joonsuk's -const bool normaliseRespYSliceForRefold = true; // keep true; THAT IS APPARENTLY REQUIRED TO REFOLD MANUALLY! even though the initial resp matrix used for the unfolding isn't normalised like this - -const bool useMatrixOverflows = false; // false by default, haven't tried true recently - -const int usePtOverflowForKineEff = 0; // false by default, not tested yet, might be better to be at 1; only matters if applyEfficiencies is 1 or 3, and by default this is not the case - -// MC split closure test: mcd as input from file_O2Analysis_MCfile_UnfoldingControl_input with Response from file_O2Analysis_MCfile_GeneralResponse -bool controlMC = true; // use file_O2Analysis_MCfile_UnfoldingControl_input MC file as input to unfolding (with h_jet_pt(_rhoareasubtracted) distrib on file), rather than real data, and as gen in comparison to gen (with h_jet_pt_part distrib on file); - - -// Debugging and checks: -const bool doManualErrorPropagForKineEff = false; // false is likely better, but hasn't been tested yet -const bool useFineBinningTest = false; -const bool drawIntermediateResponseMatrices = false; -bool comparePbPbWithRun2 = false; // if doClosure_splitMC_mcpFoldedWithFluct == true, then do the comparison with file_O2Analysis_run2ComparisonFileHannaBossiLauraFile (Nevents for this is hardcoded to what Laura told me: see mattermost discussion) - -bool setDetectorMatrixToIdentity = false; // false by default unless doing debugging -bool foldMcpWithFluct_alsoFoldWithDetResp = false; // false by default unless doing debugging; not well implemented yet - -bool smoothenEfficiency = false; -bool smoothenMCP = false; - -bool transposeResponseHistogramsInDrawing = false; // default is false; if set to true, then one can just rotate the result matrices 90 degrees to have the correct visualisation of the response as a matrix, rather than as a histogram as is default (when false) - -bool automaticBestSvdParameter = false; // automatic function not well setup yet, should work on it; keep false for now - -bool mcpInput_useMcpCollCountForUnfoldingResultNorm = true; //if controlMC is true, mcpInput_useMcpCollCountForUnfoldingResultNorm=true will normalise by N_mccoll rather than N_coll : just unfolding function for now, maybe try getmcpdistrib as well - -float ptWindowDisplay[2] = {5, 140}; // used for drawn histograms of unfolded distrib -std::array, 2> drawnWindowUnfoldedMeasurement = {{{ptWindowDisplay[0], ptWindowDisplay[1]}, {-999, -999}}}; // {{xmin, xmax}, {ymin, ymax}} - -// TODOLIST: -// - errors on ratios in drawing functions are often not good: correlation is not taken into account; at least if using mcd as input, or data and data unfolded refolded; if using data as input and dividing by mcp, then should be fine; for data/refolded, can I get the sigma matrix from roounfold maybe? -// - why is kinematic efficiency already taken into account by roounfold? applying it myself overcorrects. Where do I pass the info about it? I don't give any overflow to the response matrix (set them to 0) -// - in unfolding control test with joonsuk pp with smaller ptrec binning than ptgen to see acceptance effects, I am losing yield before 20 GeV and a little at high pt in mcp/unfolded; aactually also true if standard pt binnning without acceptance effect -// if set useMatricOverflows to true, then mcp/unfolded now equal to 1 everywhere exactly; but then refolding test fails; probably a kine eff issue, or caused byy a bad taking into account of overflows in refolding recipe -//////////////////////////////////////////////// -////////// Unfolding settings - end //////////// -//////////////////////////////////////////////// - -// Lessons: -// - if the rec distrib is cut at some 5GeV or something else, using the prior after normYslice won't work well if I don't set the binning to start AFTER this cut, not before! it works perfectly if I start at the cut - -//////////////////////////////////////////////// -//////////////// pt binning options //////////// -//////////////////////////////////////////////// -// ML paper binning -double ptBinsJetsRec_run2[nRadius][30] = {{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.},{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.},{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.}}; -int nBinPtJetsRec_run2[nRadius] = {19,19,19}; -double ptBinsJetsGen_run2[nRadius][30] = {{20., 30., 40., 50., 60., 70., 85., 100., 120., 140.},{10., 20., 40., 60., 70., 85., 100., 120., 140., 200.},{10., 20., 40., 60., 70., 85., 100., 120., 140., 200.}}; -int nBinPtJetsGen_run2[nRadius] = {9,9,9}; -double ptMaxFit = 120; - -// // pT binning for jets - gen = rec -// double ptBinsJetsRec[nRadius][30] = {{0.0, 5., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.},{0.0, 5., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.},{0.0, 5., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.}}; -// int nBinPtJetsRec[nRadius] = {15,15,15}; -// double ptBinsJetsGen[nRadius][30] = {{0.0, 5., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.},{0.0, 5., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.},{0.0, 5., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.}}; -// int nBinPtJetsGen[nRadius] = {15,15,15}; - - -// // PbPb -// // Hannah bossi identical fo run 2 comparison -// double ptBinsJetsRec[nRadius][30] = {{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.},{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.},{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.}}; -// int nBinPtJetsRec[nRadius] = {19,19,19}; -// double ptBinsJetsGen[nRadius][30] = {{10., 20., 40., 60., 70., 85., 100., 120., 140., 200.},{10., 20., 40., 60., 70., 85., 100., 120., 140., 200.},{10., 20., 40., 60., 70., 85., 100., 120., 140., 200.}}; -// int nBinPtJetsGen[nRadius] = {9,9,9}; - - -// // PbPb -// // ML paper binning for run2 comparisons -// double ptBinsJetsRec[nRadius][30] = {{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.},{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.},{20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140.}}; -// int nBinPtJetsRec[nRadius] = {19,19,19}; -// double ptBinsJetsGen[nRadius][30] = {{5., 10., 20., 30., 40., 50., 60., 70., 85., 100., 120., 140., 200.},{10., 20., 40., 60., 70., 85., 100., 120., 140., 200.},{10., 20., 40., 60., 70., 85., 100., 120., 140., 200.}}; -// int nBinPtJetsGen[nRadius] = {12,9,9}; - - -// PbPb Aimeric default -// double ptBinsJetsRec[nRadius][30] = {{10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 130., 140.},{5., 10, 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140., 200.},{5., 10, 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140., 200.}}; -// int nBinPtJetsRec[nRadius] = {22,22,22}; -// double ptBinsJetsGen[nRadius][30] = {{5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 120., 140., 200.},{0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 120., 140., 200.},{0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 120., 140., 200.}}; -// int nBinPtJetsGen[nRadius] = {13,13,13}; - - -// PbPb Aimeric test finer - lead05 -double ptBinsJetsRec[nRadius][50] = {{10., 12, 14, 16, 18, 20., 22, 24, 26, 28, 30., 32, 34, 36, 38, 40., 42, 44, 46, 48., 50., 52, 54, 56, 58, 60., 62, 64, 66, 68, 70, 75., 80., 85., 90., 95., 100., 110., 120., 130., 140.},{5., 10, 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140., 200.},{5., 10, 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140., 200.}}; -int nBinPtJetsRec[nRadius] = {40,22,22}; -double ptBinsJetsGen[nRadius][30] = {{5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95, 100., 120., 140., 160., 200.},{0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 120., 140., 200.},{0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 120., 140., 200.}}; -int nBinPtJetsGen[nRadius] = {23,13,13}; - -// // PbPb Aimeric test finer - lead03 -// double ptBinsJetsRec[nRadius][50] = {{20., 22, 24, 26, 28, 30., 32, 34, 36, 38, 40., 42, 44, 46, 48., 50., 52, 54, 56, 58, 60., 62, 64, 66, 68, 70, 75., 80., 85., 90., 95., 100., 110., 120., 130., 140.},{5., 10, 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140., 200.},{5., 10, 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100., 110., 120., 140., 200.}}; -// int nBinPtJetsRec[nRadius] = {35,22,22}; -// double ptBinsJetsGen[nRadius][30] = {{5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95, 100., 120., 140., 200.},{0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 120., 140., 200.},{0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 120., 140., 200.}}; -// int nBinPtJetsGen[nRadius] = {22,13,13}; - -// ///////////////Wenhui binning for Data unfolding GP PbPb MC///////////////////////////// -// double ptBinsJetsRec[nRadius][30] = {{-5, 0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 200}, -// {-5, 0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 200}, -// {-5, 0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 200}}; -// int nBinPtJetsRec[nRadius] = {17,17,17}; -// double ptBinsJetsGen[nRadius][30] = {{-5, 0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 200}, -// {-5, 0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 200} , -// {-5, 0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 200}}; -// int nBinPtJetsGen[nRadius] = {17,17,17}; - - -// // Joonsuk binning for pp -// double ptBinsJetsRec[nRadius][30] = {{5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140},{0, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140, 200},{5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140}}; -// int nBinPtJetsRec[nRadius] = {19,21,19}; -// double ptBinsJetsGen[nRadius][30] = {{5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140, 200},{0, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140, 200},{5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140, 200}}; -// int nBinPtJetsGen[nRadius] = {20,21,20}; - -// // Joonsuk binning for pp with smaller rec window to test kinematic efficiency -// double ptBinsJetsRec[nRadius][30] = {{18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140},{18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140},{18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140}}; -// int nBinPtJetsRec[nRadius] = {10,10,10}; -// double ptBinsJetsGen[nRadius][30] = {{5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140, 200},{5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140, 200},{5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 85, 100, 140, 200}}; -// int nBinPtJetsGen[nRadius] = {20,20,20}; - -// // pT binning for jets - gen = rec - start at 5 but rec has a smaller window ; good to check stuff without worrying about a badly setup normalisation by pt bin width -// double ptBinsJetsRec[nRadius][30] = {{30., 40., 50., 60., 70., 80., 100., 120.},{30., 40., 50., 60., 70., 80., 100., 120.},{30., 40., 50., 60., 70., 80., 100., 120.}}; -// int nBinPtJetsRec[nRadius] = {7,7,7}; -// double ptBinsJetsGen[nRadius][30] = {{0., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.},{0., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.},{0., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 100., 120., 140., 200.}}; -// int nBinPtJetsGen[nRadius] = {14,14,14}; - - -// // fine binning for pp joonsuk test files -// // int nBinPtJetsFine[nRadius] = {120,120,120}; -// int nBinPtJetsFine[nRadius] = {115,115,115}; -// // double ptBinsJetsFine[nRadius][201] = {{05., 06., 07., 08., 09., -// double ptBinsJetsFine[nRadius][201] = {{ 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., -// 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., -// 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., -// 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., -// 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., -// 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., -// 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., -// 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., -// 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., -// 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., -// 100., 105., -// 110., 115., -// 120., 125., -// 130., 135., -// 140., 145., -// 150., 155., -// 160., 165., -// 170., 175., -// 180., 185., -// 190., 195., -// 200.}, -// // {05., 06., 07., 08., 09., -// { 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., -// 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., -// 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., -// 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., -// 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., -// 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., -// 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., -// 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., -// 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., -// 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., -// 100., 105., -// 110., 115., -// 120., 125., -// 130., 135., -// 140., 145., -// 150., 155., -// 160., 165., -// 170., 175., -// 180., 185., -// 190., 195., -// 200.}, -// // {05., 06., 07., 08., 09., -// { 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., -// 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., -// 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., -// 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., -// 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., -// 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., -// 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., -// 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., -// 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., -// 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., -// 100., 105., -// 110., 115., -// 120., 125., -// 130., 135., -// 140., 145., -// 150., 155., -// 160., 165., -// 170., 175., -// 180., 185., -// 190., 195., -// 200.},}; // shift+option+left click hold lets one edit columns in vs code - - - -// // fine binning standard, for Pb-Pb non-factorised -// // int nBinPtJetsFine[nRadius] = {120,120,120}; -// int nBinPtJetsFine[nRadius] = {240,240,240}; -// // double ptBinsJetsFine[nRadius][201] = {{05., 06., 07., 08., 09., -// double ptBinsJetsFine[nRadius][402] = { -// { -200., -195., -// -190., -185., -// -180., -175., -// -170., -165., -// -160., -155., -// -150., -145., -// -140., -135., -// -130., -125., -// -120., -115., -// -110., -105., -// -100.,-99.,-98.,-97.,-96.,-95.,-94.,-93.,-92.,-91., -// -90.,-89.,-88.,-87.,-86.,-85.,-84.,-83.,-82.,-81., -// -80.,-79.,-78.,-77.,-76.,-75.,-74.,-73.,-72.,-71., -// -70.,-69.,-68.,-67.,-66.,-65.,-64.,-63.,-62.,-61., -// -60.,-59.,-58.,-57.,-56.,-55.,-54.,-53.,-52.,-51., -// -50.,-49.,-48.,-47.,-46.,-45.,-44.,-43.,-42.,-41., -// -40.,-39.,-38.,-37.,-36.,-35.,-34.,-33.,-32.,-31., -// -30.,-29.,-28.,-27.,-26.,-25.,-24.,-23.,-22.,-21., -// -20.,-19.,-18.,-17.,-16.,-15.,-14.,-13.,-12.,-11., -// -10.,-09.,-08.,-07.,-06.,-05.,-04.,-03.,-02.,-01., -// 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., -// 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., -// 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., -// 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., -// 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., -// 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., -// 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., -// 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., -// 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., -// 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., -// 100., 105., -// 110., 115., -// 120., 125., -// 130., 135., -// 140., 145., -// 150., 155., -// 160., 165., -// 170., 175., -// 180., 185., -// 190., 195., -// 200.}, -// { -200., -195., -// -190., -185., -// -180., -175., -// -170., -165., -// -160., -155., -// -150., -145., -// -140., -135., -// -130., -125., -// -120., -115., -// -110., -105., -// -100.,-99.,-98.,-97.,-96.,-95.,-94.,-93.,-92.,-91., -// -90.,-89.,-88.,-87.,-86.,-85.,-84.,-83.,-82.,-81., -// -80.,-79.,-78.,-77.,-76.,-75.,-74.,-73.,-72.,-71., -// -70.,-69.,-68.,-67.,-66.,-65.,-64.,-63.,-62.,-61., -// -60.,-59.,-58.,-57.,-56.,-55.,-54.,-53.,-52.,-51., -// -50.,-49.,-48.,-47.,-46.,-45.,-44.,-43.,-42.,-41., -// -40.,-39.,-38.,-37.,-36.,-35.,-34.,-33.,-32.,-31., -// -30.,-29.,-28.,-27.,-26.,-25.,-24.,-23.,-22.,-21., -// -20.,-19.,-18.,-17.,-16.,-15.,-14.,-13.,-12.,-11., -// -10.,-09.,-08.,-07.,-06.,-05.,-04.,-03.,-02.,-01., -// 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., -// 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., -// 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., -// 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., -// 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., -// 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., -// 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., -// 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., -// 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., -// 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., -// 100., 105., -// 110., 115., -// 120., 125., -// 130., 135., -// 140., 145., -// 150., 155., -// 160., 165., -// 170., 175., -// 180., 185., -// 190., 195., -// 200.}, -// { -200., -195., -// -190., -185., -// -180., -175., -// -170., -165., -// -160., -155., -// -150., -145., -// -140., -135., -// -130., -125., -// -120., -115., -// -110., -105., -// -100.,-99.,-98.,-97.,-96.,-95.,-94.,-93.,-92.,-91., -// -90.,-89.,-88.,-87.,-86.,-85.,-84.,-83.,-82.,-81., -// -80.,-79.,-78.,-77.,-76.,-75.,-74.,-73.,-72.,-71., -// -70.,-69.,-68.,-67.,-66.,-65.,-64.,-63.,-62.,-61., -// -60.,-59.,-58.,-57.,-56.,-55.,-54.,-53.,-52.,-51., -// -50.,-49.,-48.,-47.,-46.,-45.,-44.,-43.,-42.,-41., -// -40.,-39.,-38.,-37.,-36.,-35.,-34.,-33.,-32.,-31., -// -30.,-29.,-28.,-27.,-26.,-25.,-24.,-23.,-22.,-21., -// -20.,-19.,-18.,-17.,-16.,-15.,-14.,-13.,-12.,-11., -// -10.,-09.,-08.,-07.,-06.,-05.,-04.,-03.,-02.,-01., -// 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., -// 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., -// 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., -// 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., -// 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., -// 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., -// 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., -// 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., -// 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., -// 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., -// 100., 105., -// 110., 115., -// 120., 125., -// 130., 135., -// 140., 145., -// 150., 155., -// 160., 165., -// 170., 175., -// 180., 185., -// 190., 195., -// 200.}}; // shift+option+left click hold lets one edit columns in vs code - - - -// int nBinPtJetsFine[nRadius] = {200,200,200}; -int nBinPtJetsFine[nRadius] = {195,195,195}; -double ptBinsJetsFine[nRadius][201] = {{05., 06., 07., 08., 09., -// double ptBinsJetsFine[nRadius][201] = {{ 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., - 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., - 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., - 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., - 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., - 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., - 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., - 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., - 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., - 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., - 100.,101.,102.,103.,104.,105.,106.,107.,108.,109., - 110.,111.,112.,113.,114.,115.,116.,117.,118.,119., - 120.,121.,122.,123.,124.,125.,126.,127.,128.,129., - 130.,131.,132.,133.,134.,135.,136.,137.,138.,139., - 140.,141.,142.,143.,144.,145.,146.,147.,148.,149., - 150.,151.,152.,153.,154.,155.,156.,157.,158.,159., - 160.,161.,162.,163.,164.,165.,166.,167.,168.,169., - 170.,171.,172.,173.,174.,175.,176.,177.,178.,179., - 180.,181.,182.,183.,184.,185.,186.,187.,188.,189., - 190.,191.,192.,193.,194.,195.,196.,197.,198.,199., - 200}, - {05., 06., 07., 08., 09., - // { 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., - 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., - 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., - 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., - 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., - 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., - 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., - 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., - 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., - 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., - 100.,101.,102.,103.,104.,105.,106.,107.,108.,109., - 110.,111.,112.,113.,114.,115.,116.,117.,118.,119., - 120.,121.,122.,123.,124.,125.,126.,127.,128.,129., - 130.,131.,132.,133.,134.,135.,136.,137.,138.,139., - 140.,141.,142.,143.,144.,145.,146.,147.,148.,149., - 150.,151.,152.,153.,154.,155.,156.,157.,158.,159., - 160.,161.,162.,163.,164.,165.,166.,167.,168.,169., - 170.,171.,172.,173.,174.,175.,176.,177.,178.,179., - 180.,181.,182.,183.,184.,185.,186.,187.,188.,189., - 190.,191.,192.,193.,194.,195.,196.,197.,198.,199., - 200}, - {05., 06., 07., 08., 09., - // { 0., 01., 02., 03., 04., 05., 06., 07., 08., 09., - 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., - 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., - 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., - 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., - 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., - 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., - 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., - 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., - 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., - 100.,101.,102.,103.,104.,105.,106.,107.,108.,109., - 110.,111.,112.,113.,114.,115.,116.,117.,118.,119., - 120.,121.,122.,123.,124.,125.,126.,127.,128.,129., - 130.,131.,132.,133.,134.,135.,136.,137.,138.,139., - 140.,141.,142.,143.,144.,145.,146.,147.,148.,149., - 150.,151.,152.,153.,154.,155.,156.,157.,158.,159., - 160.,161.,162.,163.,164.,165.,166.,167.,168.,169., - 170.,171.,172.,173.,174.,175.,176.,177.,178.,179., - 180.,181.,182.,183.,184.,185.,186.,187.,188.,189., - 190.,191.,192.,193.,194.,195.,196.,197.,198.,199., - 200}}; // shift+option+left click hold lets one edit columns in vs code - - -#endif \ No newline at end of file diff --git a/Jets/JetSpectrum_systematics.C b/Jets/JetSpectrum_systematics.C index 85bac22..61b707c 100644 --- a/Jets/JetSpectrum_systematics.C +++ b/Jets/JetSpectrum_systematics.C @@ -55,6 +55,9 @@ void LoadLibs_Systematics(); void Get_systematics_UnfoldMethod(TH1D* &hSystematicUncertainty, TH1D* &hSystematicUncertainty_PreBarlow, int iDataset, int iRadius, char** unfoldingMethodList, int* unfoldParameterInputList, int nUnfoldingMethods, std::string options); void Draw_Systematics_UnfoldMethod(int iDataset, int iRadius, char** unfoldingMethodList, int* unfoldParameterInputList, int nUnfoldingMethods, std::string options); +void Draw_Systematics_parameterVariation(int iDataset, int iRadius, int unfoldIterationMin, int unfoldIterationMax, int step, std::string options); +void Draw_Systematics_TrackEff(int iDataset, int iRadius, char** unfoldingMethodList, int* unfoldParameterInputList, int nUnfoldingMethods, std::string options); + ///////////////////////////////////////////////////// ///////////////////// Main Macro //////////////////// @@ -72,20 +75,38 @@ void JetSpectrum_systematics() { // TString* Extra = new TString(""); // gathers the analysis options in a single char[] - + mcCollHistIsObsolete = inputMcCollHistIsObsolete; + /// Do not run all functions together, select which one to run by commenting/uncommenting int iDataset = 0; - int iRadius = 0; - - + int iRadius = 1; + + //######################################################### Unfolding method Systematics ##################################################### + // char optionsAnalysis_withoutUnfoldingMethod[100] = ""; + // snprintf(optionsAnalysis_withoutUnfoldingMethod, sizeof(optionsAnalysis_withoutUnfoldingMethod), "%s", unfoldingPrior); + // const int nUnfoldingMethods = 2; + // char* unfoldingMethodList[nUnfoldingMethods] = {"Svd", "Bayes"}; // default is the first one in this list + // int unfoldParameterInputList[2] = {8, 4}; + // Draw_Systematics_UnfoldMethod(iDataset, iRadius, unfoldingMethodList, unfoldParameterInputList, nUnfoldingMethods, optionsAnalysis_withoutUnfoldingMethod); + + //######################################################### Parameter variation Systematics ##################################################### + // char optionsAnalysis[100] = ""; + // snprintf(optionsAnalysis, sizeof(optionsAnalysis), "%s,%s,%s", unfoldingPrior, unfoldingMethod); + // int unfoldParameterInputMin = 7; + // int unfoldParameterInputMax = 9; + // int unfoldParameterInputStep = 1; + // Draw_Systematics_parameterVariation(iDataset, iRadius, unfoldParameterInputMin, unfoldParameterInputMax, unfoldParameterInputStep, optionsAnalysis); + + //######################################################### Track efficiency Systematics ##################################################### char optionsAnalysis_withoutUnfoldingMethod[100] = ""; snprintf(optionsAnalysis_withoutUnfoldingMethod, sizeof(optionsAnalysis_withoutUnfoldingMethod), "%s", unfoldingPrior); + const int nUnfoldingMethods = 4; + char* unfoldingMethodList[nUnfoldingMethods] = {"Svd", "Bayes", "Svd", "Bayes"}; // first two to be with nominal efficiency, last two with efficiency varied + int unfoldParameterInputList[4] = {7, 4, 4, 2}; // first two to be with nominal efficiency, last two with efficiency varied + Draw_Systematics_TrackEff(iDataset, iRadius, unfoldingMethodList, unfoldParameterInputList, nUnfoldingMethods, optionsAnalysis_withoutUnfoldingMethod); + - const int nUnfoldingMethods = 2; - char* unfoldingMethodList[nUnfoldingMethods] = {"Svd", "Bayes"}; // default is the first one in this list - int unfoldParameterInputList[2] = {8, 10}; - Draw_Systematics_UnfoldMethod(iDataset, iRadius, unfoldingMethodList, unfoldParameterInputList, nUnfoldingMethods, optionsAnalysis_withoutUnfoldingMethod); } ///////////////////////////////////////////////////// @@ -233,8 +254,6 @@ void Get_systematics_UnfoldMethod(TH1D* &hSystematicUncertainty, TH1D* &hSystema } - - void Draw_Systematics_UnfoldMethod(int iDataset, int iRadius, char** unfoldingMethodList, int* unfoldParameterInputList, int nUnfoldingMethods, std::string options) { TH1D* hSystematicUncertainty; @@ -258,9 +277,11 @@ void Draw_Systematics_UnfoldMethod(int iDataset, int iRadius, char** unfoldingMe } } - Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded, measuredInput, iDataset, iRadius, unfoldParameterInputList[0], optionsAnalysis_withUnfoldingMethod); + Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded, measuredInput, iDataset, iRadius, unfoldParameterInputList[0], optionsAnalysis_withUnfoldingMethod); //SVD unfolding as reference hSystematicUncertainty->Divide(H1D_jetPt_unfolded); //get it as a ratio of ref corrected yield + hSystematicUncertainty->Scale(100.0); hSystematicUncertainty_PreBarlow->Divide(H1D_jetPt_unfolded); //get it as a ratio of ref corrected yield + hSystematicUncertainty_PreBarlow->Scale(100.0); TString partialUniqueSpecifier = Datasets[iDataset]+"_R="+Form("%.1f",arrayRadius[iRadius])+"]_"+unfoldingMethodList[0]+"_kUnfold="+Form("%i", unfoldParameterInputList[0]); @@ -268,8 +289,189 @@ void Draw_Systematics_UnfoldMethod(int iDataset, int iRadius, char** unfoldingMe TString* pdfName = new TString("Systematics_UnfoldMethod_"+partialUniqueSpecifier); TString* pdfName_PreBarlow = new TString("Systematics_UnfoldMethod_"+partialUniqueSpecifier+"_PreBarlow"); - TString textContext(""); + // TString textContext(""); + TString textContext = Form( + "#splitline{sys. unfolding method}" + "{k_{svd} = %d, k_{bayes} = %d}", + unfoldParameterInputList[0], + unfoldParameterInputList[1] + ); + + TString* texRelativeErrPercent = new TString ("relative error (%)"); + std::array, 2> drawnWindow = {{{5, 200}, {0, 25}}}; + Draw_TH1_Histogram(hSystematicUncertainty, textContext, pdfName, texPtJetRec, texRelativeErrPercent, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, ""); + Draw_TH1_Histogram(hSystematicUncertainty_PreBarlow, textContext, pdfName_PreBarlow, texPtJetRec, texRelativeErrPercent, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, ""); +} + +void Draw_Systematics_parameterVariation(int iDataset, int iRadius, int unfoldIterationMin, int unfoldIterationMax, int step, std::string options) { + cout << "########### Drawing systematics from parameter variation ###############" << endl; + const int nUnfoldIteration = std::floor((unfoldIterationMax - unfoldIterationMin + 1)/step); + + TH1D* H1D_jetPt_unfolded[nUnfoldIteration]; + + TString partialUniqueSpecifier; + + partialUniqueSpecifier = Datasets[iDataset]+"_R="+Form("%.1f",arrayRadius[iRadius]); + + int unfoldParameterInput = 0; + + 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); + } + } + + if (measuredInput == nullptr) { + cout << "Error: measuredInput histogram is null!" << endl; + return; + } + else { + cout << "measuredInput histogram successfully retrieved." << endl; + } + + for(int iUnfoldIteration = 0; iUnfoldIteration < nUnfoldIteration; iUnfoldIteration++){ + cout << " entering the for loop " << endl; + unfoldParameterInput = unfoldIterationMax - iUnfoldIteration * step; + + cout << "((((((((((((()))))))))))))" << endl; + cout << "Iteration "<< iUnfoldIteration << endl; + cout << "((((((((((((()))))))))))))" << endl; + Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded[iUnfoldIteration], measuredInput, iDataset, iRadius, unfoldParameterInput, options); + } + + int iNominal = nUnfoldIteration/2; // choose nominal (e.g. central iteration) + TH1D* hNom = H1D_jetPt_unfolded[iNominal]; + + TH1D* hSys_envelope = (TH1D*)hNom->Clone("hSys_envelope"); + hSys_envelope->Reset(); // will store absolute systematic (positive) + + int nBins = hNom->GetNbinsX(); + for (int ib = 1; ib <= nBins; ++ib) { + double valNom = hNom->GetBinContent(ib); + double maxAbs = 0.0; + for (int i = 0; i < nUnfoldIteration; ++i) { + double val = H1D_jetPt_unfolded[i]->GetBinContent(ib); + double d = fabs(val - valNom); + if (d > maxAbs) maxAbs = d; + } + hSys_envelope->SetBinContent(ib, maxAbs); + } + + TH1D* hSys_rel = (TH1D*)hSys_envelope->Clone("hSys_relative"); + hSys_rel->Reset(); + + for (int ib = 1; ib <= nBins; ++ib) { + double absSys = hSys_envelope->GetBinContent(ib); + double valNom = hNom->GetBinContent(ib); + + double rel = 0.0; + if (valNom > 0) rel = absSys / valNom; + + hSys_rel->SetBinContent(ib, rel*100.0); // in percent + } + TString* pdfName_envolope = new TString("Systematics_deviation_ParameterVariation_"+partialUniqueSpecifier); + TString* pdfName_relUnc = new TString("Systematics_RelativeUncertainty_ParameterVariation_"+partialUniqueSpecifier); + TString textContext("SVD unfolding"); + TString* sigma = new TString("#sigma interation variation"); + TString* relativeErrors = new TString("realtive errors (%) "); + std::array, 2> drawnWindow = {{{5, 200}, {0, 1.6}}}; + // Draw_TH1_Histogram(hSys_envelope, textContext, pdfName_envolope, texPtJetRecX, sigma, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); + Draw_TH1_Histogram(hSys_rel, textContext, pdfName_relUnc, texPtJetRec, relativeErrors, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, ""); + +} + +void Draw_Systematics_TrackEff(int iDataset, int iRadius, char** unfoldingMethodList, int* unfoldParameterInputList, int nUnfoldingMethods, std::string options) { + cout << "########### Drawing systematics from track efficiency variation ###############" << endl; + TString partialUniqueSpecifier = Datasets[iDataset]+"_R="+Form("%.1f",arrayRadius[iRadius]); + + // return histogram that has the systematics in its contents + TH1D* H1D_jetPt_unfolded[nUnfoldingMethods]; + + 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); + } + } + + char optionsAnalysis_withUnfoldingMethod[100] = ""; + for(int iMethod = 0; iMethod < nUnfoldingMethods; iMethod++){ + snprintf(optionsAnalysis_withUnfoldingMethod, sizeof(optionsAnalysis_withUnfoldingMethod), "%s,%s", options.c_str(), (const char*)unfoldingMethodList[iMethod]); + int iDatasetMC = (iMethod == 0 || iMethod == 1) ? 0 : 1; // first two to be with nominal efficiency dataset 0, last two with efficiency varied dataset 1 + Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded[iMethod], measuredInput, iDatasetMC, iRadius, unfoldParameterInputList[iMethod], optionsAnalysis_withUnfoldingMethod); + } + + // Get histograms + TH1D* H1D_Nominal_SVD = (TH1D*) H1D_jetPt_unfolded[0]->Clone("H1D_Nominal_SVD"); + TH1D* H1D_Nominal_Bayes = (TH1D*) H1D_jetPt_unfolded[1]->Clone("H1D_Nominal_Bayes"); + TH1D* H1D_Reduced_SVD = (TH1D*) H1D_jetPt_unfolded[2]->Clone("H1D_Reduced_SVD"); + TH1D* H1D_Reduced_Bayes = (TH1D*) H1D_jetPt_unfolded[3]->Clone("H1D_Reduced_Bayes"); + + // Create histograms for absolute differences and relative uncertainties + TH1D *H1D_Delta_SVD = (TH1D*)H1D_Nominal_SVD->Clone("hDiff_NominalReduced_SVD"); + H1D_Delta_SVD->Reset(); + TH1D *H1D_RelativeUncertainty_SVD = (TH1D*)H1D_Nominal_SVD->Clone("H1D_RelativeUncertainty_SVD"); + H1D_RelativeUncertainty_SVD->Reset(); + + TH1D *H1D_Delta_Bayes = (TH1D*)H1D_Nominal_Bayes->Clone("H1D_Delta_Bayes"); + H1D_Delta_Bayes->Reset(); + TH1D *H1D_RelativeUncertainty_Bayes = (TH1D*)H1D_Nominal_Bayes->Clone("H1D_RelativeUncertainty_Bayes"); + H1D_RelativeUncertainty_Bayes->Reset(); + + // Compute absolute differences bin by bin and compute relative uncertainties + for (int i = 1; i <= H1D_Nominal_SVD->GetNbinsX(); ++i) { + H1D_Delta_SVD->SetBinContent(i, abs(H1D_Nominal_SVD->GetBinContent(i) - H1D_Reduced_SVD->GetBinContent(i))); + H1D_Delta_SVD->SetBinError(i, 0.0); + + if (H1D_Nominal_SVD->GetBinContent(i) != 0) { + double relUnc = (H1D_Delta_SVD->GetBinContent(i) / H1D_Nominal_SVD->GetBinContent(i)) * 100.0; // in percent + H1D_RelativeUncertainty_SVD->SetBinContent(i, relUnc); + H1D_RelativeUncertainty_SVD->SetBinError(i, 0.0); + } else { + H1D_RelativeUncertainty_SVD->SetBinContent(i, 0.0); + H1D_RelativeUncertainty_SVD->SetBinError(i, 0.0); + } + } + + for (int i = 1; i <= H1D_Nominal_Bayes->GetNbinsX(); ++i) { + H1D_Delta_Bayes->SetBinContent(i, abs(H1D_Nominal_Bayes->GetBinContent(i) - H1D_Reduced_Bayes->GetBinContent(i))); + H1D_Delta_Bayes->SetBinError(i, 0.0); + + if (H1D_Nominal_Bayes->GetBinContent(i) != 0) { + double relUnc = (H1D_Delta_Bayes->GetBinContent(i) / H1D_Nominal_Bayes->GetBinContent(i)) * 100.0; // in percent + H1D_RelativeUncertainty_Bayes->SetBinContent(i, relUnc); + H1D_RelativeUncertainty_Bayes->SetBinError(i, 0.0); + } else { + H1D_RelativeUncertainty_Bayes->SetBinContent(i, 0.0); + H1D_RelativeUncertainty_Bayes->SetBinError(i, 0.0); + } + } + TString textContextSVD("with 3% reduced track efficiency (SVD)"); + TString textContextBayes("with 3% reduced track efficiency (Bayes)"); + + TString* pdfName_relUncSvd_logx = new TString("Relative_uncert_Svd_Data_train380686_Nominal_train533385_param7_ReducedTrackEff3perCent_train564527_param4_R=0.4_logx"); + TString* pdfName_relUncBayes_logx = new TString("Relative_uncert_Bayes_Data_train380686_Nominal_train533385_param4_ReducedTrackEff3perCent_train564527_param2_R=0.4_logx"); + TString* pdfName_relUncSvd = new TString("Relative_uncert_Svd_Data_train380686_Nominal_train533385_param7_ReducedTrackEff3perCent_train564527_param4_R=0.4"); + TString* pdfName_relUncBayes = new TString("Relative_uncert_Bayes_Data_train380686_Nominal_train533385_param4_ReducedTrackEff3perCent_train564527_param2_R=0.4"); + std::array, 2> drawnWindow = {{{5, 200}, {0, 80}}}; + + Draw_TH1_Histogram(H1D_RelativeUncertainty_SVD, textContextSVD, pdfName_relUncSvd_logx, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, "logx"); + Draw_TH1_Histogram(H1D_RelativeUncertainty_Bayes, textContextBayes, pdfName_relUncBayes_logx, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, "logx"); + Draw_TH1_Histogram(H1D_RelativeUncertainty_SVD, textContextSVD, pdfName_relUncSvd, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, ""); + Draw_TH1_Histogram(H1D_RelativeUncertainty_Bayes, textContextBayes, pdfName_relUncBayes, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, ""); + - Draw_TH1_Histogram(hSystematicUncertainty, textContext, pdfName, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); - Draw_TH1_Histogram(hSystematicUncertainty_PreBarlow, textContext, pdfName_PreBarlow, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); } \ No newline at end of file From 0d272d1157e32a58c5075a814f6ad609e5e5874f Mon Sep 17 00:00:00 2001 From: ahmadtabikh Date: Thu, 29 Jan 2026 15:01:33 +0100 Subject: [PATCH 2/4] fit --- Jets/JetSpectrum_systematics.C | 505 ++++++++++++++++++++++++++++++++- 1 file changed, 498 insertions(+), 7 deletions(-) diff --git a/Jets/JetSpectrum_systematics.C b/Jets/JetSpectrum_systematics.C index 61b707c..9451a1a 100644 --- a/Jets/JetSpectrum_systematics.C +++ b/Jets/JetSpectrum_systematics.C @@ -45,6 +45,7 @@ #include #include #include +#include using namespace std; // Misc utilities @@ -57,6 +58,7 @@ void Get_systematics_UnfoldMethod(TH1D* &hSystematicUncertainty, TH1D* &hSystema void Draw_Systematics_UnfoldMethod(int iDataset, int iRadius, char** unfoldingMethodList, int* unfoldParameterInputList, int nUnfoldingMethods, std::string options); void Draw_Systematics_parameterVariation(int iDataset, int iRadius, int unfoldIterationMin, int unfoldIterationMax, int step, std::string options); void Draw_Systematics_TrackEff(int iDataset, int iRadius, char** unfoldingMethodList, int* unfoldParameterInputList, int nUnfoldingMethods, std::string options); +void Draw_Systematics_SecondaryContamination(int iDataset, int iRadius, int unfoldParameterInput, std::string options); ///////////////////////////////////////////////////// @@ -97,13 +99,18 @@ void JetSpectrum_systematics() { // Draw_Systematics_parameterVariation(iDataset, iRadius, unfoldParameterInputMin, unfoldParameterInputMax, unfoldParameterInputStep, optionsAnalysis); //######################################################### Track efficiency Systematics ##################################################### - char optionsAnalysis_withoutUnfoldingMethod[100] = ""; - snprintf(optionsAnalysis_withoutUnfoldingMethod, sizeof(optionsAnalysis_withoutUnfoldingMethod), "%s", unfoldingPrior); - const int nUnfoldingMethods = 4; - char* unfoldingMethodList[nUnfoldingMethods] = {"Svd", "Bayes", "Svd", "Bayes"}; // first two to be with nominal efficiency, last two with efficiency varied - int unfoldParameterInputList[4] = {7, 4, 4, 2}; // first two to be with nominal efficiency, last two with efficiency varied - Draw_Systematics_TrackEff(iDataset, iRadius, unfoldingMethodList, unfoldParameterInputList, nUnfoldingMethods, optionsAnalysis_withoutUnfoldingMethod); + // char optionsAnalysis_withoutUnfoldingMethod[100] = ""; + // snprintf(optionsAnalysis_withoutUnfoldingMethod, sizeof(optionsAnalysis_withoutUnfoldingMethod), "%s", unfoldingPrior); + // const int nUnfoldingMethods = 4; + // char* unfoldingMethodList[nUnfoldingMethods] = {"Svd", "Bayes", "Svd", "Bayes"}; // first two to be with nominal efficiency, last two with efficiency varied + // int unfoldParameterInputList[4] = {7, 4, 4, 2}; // first two to be with nominal efficiency, last two with efficiency varied + // Draw_Systematics_TrackEff(iDataset, iRadius, unfoldingMethodList, unfoldParameterInputList, nUnfoldingMethods, optionsAnalysis_withoutUnfoldingMethod); + //######################################################### Secondary tracks Systematics ##################################################### + char optionsAnalysis[100] = ""; + snprintf(optionsAnalysis, sizeof(optionsAnalysis), "%s,%s,%s", unfoldingPrior, unfoldingMethod); + int unfoldParameterInput = 7; + Draw_Systematics_SecondaryContamination(iDataset, iRadius, unfoldParameterInput, optionsAnalysis); @@ -113,6 +120,130 @@ void JetSpectrum_systematics() { /////////////////// Misc utilities ////////////////// ///////////////////////////////////////////////////// + +///////////// Fit functions ////////////////////////////////// +// Convert a fitted function + fit result into a TGraphErrors that includes the 1σ confidence band of the fit. +TGraphErrors* getFunctionTGraphErrorsFromFitResult(double* xRangeFit, TF1* fitFunctionDrawn, TFitResultPtr fitResult, int nPointsGraph = 1000){ + std::vector xAxisGraph= {}; + std::vector yAxisGraph= {}; + std::vector yAxisGraphErrors= {}; + // double* ; + + for(int iPoint = 0; iPoint < nPointsGraph; iPoint++){ + xAxisGraph.push_back(xRangeFit[0]+iPoint*1./nPointsGraph*(xRangeFit[1]-xRangeFit[0])); + yAxisGraph.push_back(fitFunctionDrawn->Eval(xAxisGraph.back())); + yAxisGraphErrors.push_back(0); + } + double oneSigmaInterval = 0.683; + fitResult->GetConfidenceIntervals(nPointsGraph, 1, 1, &xAxisGraph[0], &yAxisGraphErrors[0], oneSigmaInterval, false); + TGraphErrors* fitFunctionTGraphErrors = new TGraphErrors(nPointsGraph, &xAxisGraph[0], &yAxisGraph[0], nullptr, &yAxisGraphErrors[0]); + return fitFunctionTGraphErrors; +} + +// Fit a histogram with a double Tsallis-like function and return everything needed to propagate uncertainties +std::tuple FitDoubleTsallis(TH1D* &histogramInput, int nBinsX, double* binsX, double* xRangeFit) { + TF1 *fitFunctionInit; + TF1 *fitFunctionFinal; + TF1 *fitFunctionDrawn; // drawn over the full range + TFitResultPtr fFitResult; + + double parfitFunctionInit[4]; + double parfitFunctionFinal[4]; + // double parfitFunctionInit[8]; + // double parfitFunctionFinal[8]; + // const char* doubleTsallis = "([2]+[3]*x)*pow(1 + x/([0]*[1]), -[1]) + ([6]+[7]*x)*pow(1 + x/([4]*[5]), -[5])"; + const char* doubleTsallis = "([2]+[3]*x)*pow(1 + x/([0]*[1]), -[1])"; + + + //////////////////////////////////////////////////////////////////// + //////////////////////////// Fit start ///////////////////////////// + //////////////////////////////////////////////////////////////////// + + fitFunctionInit = new TF1("fitFunctionInit_", doubleTsallis, xRangeFit[0], xRangeFit[1]); + + // Set parameter names + fitFunctionInit->SetParName(0, "p0"); + fitFunctionInit->SetParName(1, "p1"); + fitFunctionInit->SetParName(2, "p2"); + fitFunctionInit->SetParName(3, "p3"); + // fitFunctionInit->SetParName(4, "p4"); + // fitFunctionInit->SetParName(5, "p5"); + // fitFunctionInit->SetParName(6, "p6"); + // fitFunctionInit->SetParName(7, "p7"); + + // fitFunctionInit->SetParameters(0.5, 7, 50, 0, 1.2, 10, 300, 0); + // // p0, p1, p2, p3, p4, p5, p6, p7 + + fitFunctionInit->SetParameters(0.5, 7, 0, 0); + // p0, p1, p2, p3 + + fitFunctionInit->SetParLimits(0, 0.05, 1.0); + fitFunctionInit->SetParLimits(1, 3.0, 10.0); + fitFunctionInit->SetParLimits(2, -20.0, 2.0); + fitFunctionInit->SetParLimits(3, -10.0, 70.0); + // fitFunctionInit->SetParLimits(4, 0.05, 5.0); + // fitFunctionInit->SetParLimits(5, 3.0, 30.0); + // fitFunctionInit->SetParLimits(6, -50.0, 500.0); + // fitFunctionInit->SetParLimits(7, -100.0, 100.0); + + histogramInput->Fit(fitFunctionInit, "R0Q"); // R = fit range, Q = quiet, L = likelihood + fitFunctionInit->GetParameters(&parfitFunctionInit[0]); // Save initial parameters + + fitFunctionFinal = new TF1("fitFunctionFinal_", doubleTsallis, xRangeFit[0], xRangeFit[1]); + + for(int i=0; i<8; i++) fitFunctionFinal->SetParameter(i, parfitFunctionInit[i]); + + fFitResult = histogramInput->Fit(fitFunctionFinal, "RS"); + fitFunctionFinal->GetParameters(&parfitFunctionFinal[0]); + + // Check covariance availability + TMatrixDSym covMatrixFit; // default empty + if (fFitResult && fFitResult->CovMatrixStatus() == 3) { + covMatrixFit = fFitResult->GetCovarianceMatrix(); + } else { + std::cout << "Warning: Covariance matrix not available!" << std::endl; + } + + // TMatrixDSym covMatrixFit = fFitResult->GetCovarianceMatrix(); + + // Double_t *pDataSmall = covMatrixFit.GetMatrixArray(); + // for (int i = 0; i < 2*2; i++) { + // cout << "i = " << i << ", covMatrixFit[i]" << pDataSmall[i] << endl; + // } + + fitFunctionDrawn = new TF1("fitFunctionDrawn_", doubleTsallis, xRangeFit[0], xRangeFit[1]); + for(int i=0; i<8; i++) fitFunctionDrawn->SetParameter(i, parfitFunctionFinal[i]); + + std::tuple fitFunctionAndFitParams(fitFunctionDrawn, covMatrixFit, fFitResult); + return fitFunctionAndFitParams; +} + +// Use the fitted function to rebin a histogram and propagate fit uncertainties to the new bins +std::tuple RebinWithDoubleTsallisFit(TH1D* &histogramInput, int nBinsX, double* binsX, double* xRangeFit) { + std::tuple tsallisFitFunctionResult = FitDoubleTsallis(histogramInput, nBinsX, binsX, xRangeFit); + TF1* fitFunctionDrawn = std::get<0>(tsallisFitFunctionResult); + TFitResultPtr fitResult = std::get<2>(tsallisFitFunctionResult); + TGraphErrors* fitFunctionTGraphErrors = getFunctionTGraphErrorsFromFitResult(xRangeFit, fitFunctionDrawn, fitResult); + + //////////////////////////// Rebin of input histogram ///////////////////////////// + + TH1D* histogramRebinned = new TH1D("Unfolded: fit sampling", "Unfolded: fit sampling", nBinsX, binsX); + for(int iBin = 0; iBin < nBinsX; iBin++){ + double xCenter = histogramRebinned->GetXaxis()->GetBinCenter(iBin); + histogramRebinned->SetBinContent(iBin, fitFunctionDrawn->Eval(xCenter)); + double oneSigmaInterval = 0.683; + double errorEval[1] = {0}; + double xEval[1] = {xCenter}; + fitResult->GetConfidenceIntervals(1, 1, 1, xEval, errorEval, oneSigmaInterval, false); + histogramRebinned->SetBinError(iBin, errorEval[0]); + } + + std::tuple rebinResultAndFitFunction(histogramRebinned, fitFunctionTGraphErrors, fitFunctionDrawn); + return rebinResultAndFitFunction; +} + +////////////////////////////////////////////////////////////// + void LoadLibs_Systematics() { // gSystem->Load("libCore.so"); // gSystem->Load("libGeom.so"); @@ -474,4 +605,364 @@ void Draw_Systematics_TrackEff(int iDataset, int iRadius, char** unfoldingMethod Draw_TH1_Histogram(H1D_RelativeUncertainty_Bayes, textContextBayes, pdfName_relUncBayes, texPtJetRec, texSystematicsPercent, texCollisionDataInfo, drawnWindow, legendPlacementAuto, contextPlacementAuto, ""); -} \ No newline at end of file +} + +void Draw_Systematics_SecondaryContamination(int iDataset, int iRadius, int unfoldParameterInput, std::string options){ + cout << "########### Drawing systematics from secondary contamination variation ###############" << endl; + TH1D* H1D_jetPt_unfolded; + + 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); + } + } + + Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded, measuredInput, iDataset, iRadius, unfoldParameterInput, options); + + // Define your bins and fit range + int nBinsX = H1D_jetPt_unfolded->GetNbinsX(); + double* binsX = new double[nBinsX+1]; + for(int i=0; i<=nBinsX; i++) binsX[i] = H1D_jetPt_unfolded->GetBinLowEdge(i+1); + + double xRangeFit[2] = {5.0, 120.0}; // Fit range in GeV + + // Step 1: Rebin histogram using double Tsallis fit + std::tuple result = + RebinWithDoubleTsallisFit(H1D_jetPt_unfolded, nBinsX, binsX, xRangeFit); + + // Step 2: Extract outputs + TH1D* hJetPtRebinned = std::get<0>(result); + TGraphErrors* fitGraph = std::get<1>(result); + TF1* fitFunctionDrawn = std::get<2>(result); + + // Step 3: Draw original histogram and rebinned fit + TCanvas* c1 = new TCanvas("c1", "Double Tsallis Fit", 800, 600); + H1D_jetPt_unfolded->SetMarkerStyle(20); + H1D_jetPt_unfolded->SetMarkerColor(kBlack); + H1D_jetPt_unfolded->Draw("E"); // original histogram with errors + + fitGraph->SetLineColor(kRed); + fitGraph->SetLineWidth(2); + fitGraph->Draw("L SAME"); // smooth fit with ±1σ band + + hJetPtRebinned->SetMarkerStyle(24); + hJetPtRebinned->SetMarkerColor(kBlue); + hJetPtRebinned->Draw("E SAME"); // rebinned histogram + + // c1->BuildLegend(); + // ================= Legend ================= + TLegend* leg = new TLegend(0.55, 0.65, 0.85, 0.85); + leg->SetBorderSize(0); + leg->SetFillStyle(0); // transparent + leg->SetTextSize(0.035); + + leg->AddEntry(H1D_jetPt_unfolded, "Unfolded data", "lep"); + leg->AddEntry(fitGraph, "Double Tsallis fit", "l"); + leg->AddEntry(hJetPtRebinned, "Fit sampling (rebinned)", "lep"); + + leg->Draw(); + c1->Update(); + + + // TF1* ShiftTF1(TF1* f, double shift, const char* name="shifted") { + // return new TF1(name, [f, shift](double *x, double *){ return f->Eval(x[0] * shift); }, + // f->GetXmin(), f->GetXmax(), 0); + // } + + // double shiftFactor = 0.005; + // TF1* fUp = ShiftTF1(fitFunctionDrawn, 1.0 + shiftFactor, "fUp"); // 1.005 * pT + // TF1* fDown = ShiftTF1(fitFunctionDrawn, 1.0 - shiftFactor, "fDown"); // 0.995 * pT + + // TH1D* hRebinnedUp = new TH1D("hRebinnedUp", "Rebinned Up (1.005x)", nBinsX, binsX); + // TH1D* hRebinnedDown = new TH1D("hRebinnedDown", "Rebinned Down (0.995x)", nBinsX, binsX); + + // for(int iBin = 0; iBin < nBinsX; iBin++){ + // double xCenter = hRebinnedUp->GetXaxis()->GetBinCenter(iBin); + // hRebinnedUp->SetBinContent(iBin, fUp->Eval(xCenter)); + // hRebinnedDown->SetBinContent(iBin, fDown->Eval(xCenter)); + // } + + // // --- Create histograms for absolute differences --- + // TH1D* hDiffUp = new TH1D("hDiffUp", "Absolute difference Up", nBinsX, binsX); + // TH1D* hDiffDown = new TH1D("hDiffDown", "Absolute difference Down", nBinsX, binsX); + + // // --- Fill the difference histograms --- + // for(int iBin = 0; iBin < nBinsX; iBin++){ + // double nominal = hJetPtRebinned->GetBinContent(iBin); + // if(nominal == 0) nominal = 1e-12; // avoid division by zero + + // double upVal = hRebinnedUp->GetBinContent(iBin); + // double downVal = hRebinnedDown->GetBinContent(iBin); + + // hRatioUp->SetBinContent(iBin, fabs(upVal - nominal) / nominal * 100.0); + // hRatioDown->SetBinContent(iBin, fabs(downVal - nominal) / nominal * 100.0); + // } + + // TH1D** deltaHistos = new TH1D*[2]; + // deltaHistos[0] = hRatioUp; // Up variation + // deltaHistos[1] = hRatioDown; // Down variation + + // const TString Names[2] = { + // TString::Format("Up variation (+%.2f%%)", shiftFactor*100.0), + // TString::Format("Down variation (-%.2f%%)", shiftFactor*100.0) + // }; + + // TString pdfName = TString::Format("jet_spectrum_systematics_SecondaryTrackContamination_upDownVariation+%.3f.pdf", shiftFactor); + // std::array, 2> drawnWindow = {{{-25, 200},{0.85, 1.15}}}; + // std::array, 2> legendPlacement = {{{0.5, 0.7}, {0.75, 0.90}}}; // {{{x1, y1}, {x2, y2}}} + // Draw_TH1_Histograms(deltaHistos, Names, 2, textContext, pdfNamePt, texPtJetRec , texSystematicsPercent, texCollisionDataInfo, drawnWindow, legendPlacement, contextPlacementAuto, ""); +} + +/* +void Draw_Systematics_SecondaryContamination(int iDataset, int iRadius, int unfoldParameterInput, const double* xRangeFit, std::string options); + + +///////////////////////////////////////////////////// +///////////////////// Main Macro //////////////////// +///////////////////////////////////////////////////// + +void JetSpectrum_systematics() { + + int iDataset = 0; + int iRadius = 1; + + //######################################################### Secondary tracks Systematics ##################################################### + char optionsAnalysis[100] = ""; + snprintf(optionsAnalysis, sizeof(optionsAnalysis), "%s,%s,%s", unfoldingPrior, unfoldingMethod); + int unfoldParameterInput = 7; + const double xRangeFit[2] = {5.0, 120.0}; // Fit range in GeV + Draw_Systematics_SecondaryContamination(iDataset, iRadius, unfoldParameterInput, xRangeFit, optionsAnalysis); + +} + + +///////////// Fit functions ////////////////////////////////// +// Convert a fitted function + fit result into a TGraphErrors that includes the 1σ confidence band of the fit. +TGraphErrors* getFunctionTGraphErrorsFromFitResult(const double* xRangeFit, TF1* fitFunctionDrawn, TFitResultPtr fitResult, int nPointsGraph = 1000){ + std::vector xAxisGraph= {}; + std::vector yAxisGraph= {}; + std::vector yAxisGraphErrors= {}; + // double* ; + + for(int iPoint = 0; iPoint < nPointsGraph; iPoint++){ + xAxisGraph.push_back(xRangeFit[0]+iPoint*1./nPointsGraph*(xRangeFit[1]-xRangeFit[0])); + yAxisGraph.push_back(fitFunctionDrawn->Eval(xAxisGraph.back())); + yAxisGraphErrors.push_back(0); + } + double oneSigmaInterval = 0.683; + fitResult->GetConfidenceIntervals(nPointsGraph, 1, 1, &xAxisGraph[0], &yAxisGraphErrors[0], oneSigmaInterval, false); + TGraphErrors* fitFunctionTGraphErrors = new TGraphErrors(nPointsGraph, &xAxisGraph[0], &yAxisGraph[0], nullptr, &yAxisGraphErrors[0]); + return fitFunctionTGraphErrors; +} + +// Fit a histogram with a double Tsallis-like function and return everything needed to propagate uncertainties +std::tuple FitDoubleTsallis(TH1D* &histogramInput, int nBinsX, double* binsX, const double* xRangeFit) { + // ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(10000); + // ROOT::Math::MinimizerOptions::SetDefaultTolerance(1e-4); + const char* doubleTsallis = "([2]+[3]*x)*pow(1 + x/([0]*[1]), -[1]) + ([6]+[7]*x)*pow(1 + x/([4]*[5]), -[5])"; + + TF1 *fitFunctionInit; + TF1 *fitFunctionFinal; + TF1 *fitFunctionDrawn; // drawn over the full range + TFitResultPtr fFitResult; + double parInit[8]; + double parFinal[8]; + + // Initial Fit + fitFunctionInit = new TF1("fitFunctionInit_", doubleTsallis, xRangeFit[0], xRangeFit[1]); + // Set parameter names + fitFunctionInit->SetParName(0, "p0"); + fitFunctionInit->SetParName(1, "p1"); + fitFunctionInit->SetParName(2, "p2"); + fitFunctionInit->SetParName(3, "p3"); + fitFunctionInit->SetParName(4, "p4"); + fitFunctionInit->SetParName(5, "p5"); + fitFunctionInit->SetParName(6, "p6"); + fitFunctionInit->SetParName(7, "p7"); + + fitFunctionInit->SetParameters(0.5, 7, -10, 50, 0.5, 5, 300, -70); + // p0, p1, p2, p3, p4, p5, p6, p7 + fitFunctionInit->SetParLimits(0, 0.05, 5.0); + fitFunctionInit->SetParLimits(1, 3.0, 30.0); + fitFunctionInit->SetParLimits(2, -20.0, 20.0); + fitFunctionInit->SetParLimits(3, -10.0, 70.0); + fitFunctionInit->SetParLimits(4, 0.05, 5.0); + fitFunctionInit->SetParLimits(5, 3.0, 30.0); + fitFunctionInit->SetParLimits(6, -50.0, 500.0); + fitFunctionInit->SetParLimits(7, -100.0, 100.0); + + histogramInput->Fit(fitFunctionInit, "R0QL"); // R = fit range, Q = quiet, L = likelihood + fitFunctionInit->GetParameters(&parInit[0]); // Save initial parameters + + // Final Fit + fitFunctionFinal = new TF1("fitFunctionFinal_", doubleTsallis, xRangeFit[0], xRangeFit[1]); + for(int i=0; i<8; i++) fitFunctionFinal->SetParameter(i, parInit[i]); + + fFitResult = histogramInput->Fit(fitFunctionFinal, "RSLH"); // S : return a TFitResultPtr (crucial for errors, covariance ...) + fitFunctionFinal->GetParameters(&parFinal[0]); + // Check covariance availability + TMatrixDSym covMatrixFit; // default empty + if (fFitResult && fFitResult->CovMatrixStatus() == 3) { // 0 : not calculated, 1 : approximated, 2 : forced pos. def., 3 : accurate + covMatrixFit = fFitResult->GetCovarianceMatrix(); + } else { + std::cout << "Warning: Covariance matrix not available!" << std::endl; + } + // TMatrixDSym covMatrixFit = fFitResult->GetCovarianceMatrix(); + // Double_t *pDataSmall = covMatrixFit.GetMatrixArray(); + // for (int i = 0; i < 2*2; i++) { + // cout << "i = " << i << ", covMatrixFit[i]" << pDataSmall[i] << endl; + // } + fitFunctionDrawn = new TF1("fitFunctionDrawn_", doubleTsallis, xRangeFit[0], xRangeFit[1]); + for(int i=0; i<8; i++) fitFunctionDrawn->SetParameter(i, parFinal[i]); + + std::tuple fitFunctionAndFitParams(fitFunctionDrawn, covMatrixFit, fFitResult); + return fitFunctionAndFitParams; +} + +// Use the fitted function to rebin a histogram and propagate fit uncertainties to the new bins +std::tuple RebinWithDoubleTsallisFit(TH1D* &histogramInput, int nBinsX, double* binsX, const double* xRangeFit) { + std::tuple tsallisFitFunctionResult = FitDoubleTsallis(histogramInput, nBinsX, binsX, xRangeFit); + TF1* fitFunctionDrawn = std::get<0>(tsallisFitFunctionResult); + TFitResultPtr fitResult = std::get<2>(tsallisFitFunctionResult); + TGraphErrors* fitFunctionTGraphErrors = getFunctionTGraphErrorsFromFitResult(xRangeFit, fitFunctionDrawn, fitResult); + + //////////////////////////// Rebin of input histogram ///////////////////////////// + + TH1D* histogramRebinned = new TH1D("Unfolded: fit sampling", "Unfolded: fit sampling", nBinsX, binsX); + for(int iBin = 0; iBin < nBinsX; iBin++){ + double xCenter = histogramRebinned->GetXaxis()->GetBinCenter(iBin); + histogramRebinned->SetBinContent(iBin, fitFunctionDrawn->Eval(xCenter)); + double oneSigmaInterval = 0.683; + double errorEval[1] = {0}; + double xEval[1] = {xCenter}; + fitResult->GetConfidenceIntervals(1, 1, 1, xEval, errorEval, oneSigmaInterval, false); + histogramRebinned->SetBinError(iBin, errorEval[0]); + } + + std::tuple rebinResultAndFitFunction(histogramRebinned, fitFunctionTGraphErrors, fitFunctionDrawn); + return rebinResultAndFitFunction; +} + +TF1* ShiftTF1(const TF1* f, double shift, const char* name="shifted") { + return new TF1(name, [f, shift](double *x, double *){ return f->Eval(x[0] * shift); }, + f->GetXmin(), f->GetXmax(), 0); + } +///////////////////////////////////////////////////// +void Draw_Systematics_SecondaryContamination(int iDataset, int iRadius, int unfoldParameterInput, const double* xRangeFit, std::string options){ + cout << "########### Drawing systematics from secondary contamination variation ###############" << endl; + TH1D* H1D_jetPt_unfolded; + + 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); + } + } + + Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded, measuredInput, iDataset, iRadius, unfoldParameterInput, options); + + // Define your bins and fit range + int nBinsX = H1D_jetPt_unfolded->GetNbinsX(); + double* binsX = new double[nBinsX+1]; + for(int i=0; i<=nBinsX; i++) binsX[i] = H1D_jetPt_unfolded->GetBinLowEdge(i+1); + + // double xRangeFit[2] = {5.0, 120.0}; // Fit range in GeV + cout << "########### Fit range: [" << xRangeFit[0] << " , " << xRangeFit[1] << "] GeV #############" << endl; + // Step 1: Rebin histogram using double Tsallis fit + std::tuple result = RebinWithDoubleTsallisFit(H1D_jetPt_unfolded, nBinsX, binsX, xRangeFit); + + // Step 2: Extract outputs + TH1D* hJetPtRebinned = std::get<0>(result); + TGraphErrors* fitGraph = std::get<1>(result); + TF1* fitFunctionDrawn = std::get<2>(result); + + // Step 3: Draw original histogram and rebinned fit + TCanvas* c1 = new TCanvas("c1", "Double Tsallis Fit", 800, 600); + H1D_jetPt_unfolded->SetMarkerStyle(20); + H1D_jetPt_unfolded->SetMarkerColor(kBlack); + H1D_jetPt_unfolded->Draw("E"); // original histogram with errors + + fitGraph->SetLineColor(kRed); + fitGraph->SetLineWidth(2); + fitGraph->Draw("L SAME"); // smooth fit with ±1σ band + + hJetPtRebinned->SetMarkerStyle(24); + hJetPtRebinned->SetMarkerColor(kBlue); + hJetPtRebinned->Draw("E SAME"); // rebinned histogram + + // c1->BuildLegend(); + // ================= Legend ================= + TLegend* leg = new TLegend(0.55, 0.65, 0.85, 0.85); + leg->SetBorderSize(0); + leg->SetFillStyle(0); // transparent + leg->SetTextSize(0.035); + + leg->AddEntry(H1D_jetPt_unfolded, "Unfolded data", "lep"); + leg->AddEntry(fitGraph, "Double Tsallis fit", "l"); + leg->AddEntry(hJetPtRebinned, "Fit sampling (rebinned)", "lep"); + + leg->Draw(); + c1->Update(); + + // double shiftFactor = 0.005; + // TF1* fUp = ShiftTF1(fitFunctionDrawn, 1.0 + shiftFactor, "fUp"); // 1.005 * pT + // TF1* fDown = ShiftTF1(fitFunctionDrawn, 1.0 - shiftFactor, "fDown"); // 0.995 * pT + + // TString titleUp = Form("Rebinned Up (+%.2f%%)", shiftFactor * 100.0); + // TString titleDown = Form("Rebinned Down (-%.2f%%)", shiftFactor * 100.0); + + // TH1D* hRebinnedUp = new TH1D("hRebinnedUp", titleUp, nBinsX, binsX); + // TH1D* hRebinnedDown = new TH1D("hRebinnedDown", titleDown, nBinsX, binsX); + + // for(int iBin = 1; iBin < nBinsX; iBin++){ + // double xCenter = hRebinnedUp->GetXaxis()->GetBinCenter(iBin); + // hRebinnedUp->SetBinContent(iBin, fUp->Eval(xCenter)); + // hRebinnedDown->SetBinContent(iBin, fDown->Eval(xCenter)); + // } + + // // --- Create histograms for absolute differences --- + // TH1D* hRatioUp = new TH1D("hRatioUp", "Ratio difference Up", nBinsX, binsX); + // TH1D* hRatioDown = new TH1D("hRatioDown", "Ratio difference Down", nBinsX, binsX); + + // // --- Fill the difference histograms --- + // for(int iBin = 0; iBin < nBinsX; iBin++){ + // double nominal = hJetPtRebinned->GetBinContent(iBin); + // if(nominal == 0) nominal = 1e-12; // avoid division by zero + + // double upVal = hRebinnedUp->GetBinContent(iBin); + // double downVal = hRebinnedDown->GetBinContent(iBin); + + // hRatioUp->SetBinContent(iBin, fabs(upVal - nominal) / nominal * 100.0); + // hRatioDown->SetBinContent(iBin, fabs(downVal - nominal) / nominal * 100.0); + // } + + // TH1D** deltaHistos = new TH1D*[2]; + // deltaHistos[0] = hRatioUp; // Up variation + // deltaHistos[1] = hRatioDown; // Down variation + + // const TString Names[2] = { + // TString::Format("Up variation (+%.2f%%)", shiftFactor*100.0), + // TString::Format("Down variation (-%.2f%%)", shiftFactor*100.0) + // }; + + // TString pdfName = TString::Format("jet_spectrum_systematics_SecondaryTrackContamination_upDownVariation+%.3f.pdf", shiftFactor); + // std::array, 2> drawnWindow = {{{-25, 200},{0.85, 1.15}}}; + // std::array, 2> legendPlacement = {{{0.5, 0.7}, {0.75, 0.90}}}; // {{{x1, y1}, {x2, y2}}} + // TString textContext = ""; + // TString* texCollisionDataInfo = new TString("pp #sqrt{#it{s}} = 5.36 TeV"); + // Draw_TH1_Histograms(deltaHistos, Names, 2, textContext, pdfName, texPtJetRec , texSystematicsPercent, texCollisionDataInfo, drawnWindow, legendPlacement, contextPlacementAuto, ""); +} + */ \ No newline at end of file From ee8bda28902eac4e8b2221d2f0e5abfe165b9b90 Mon Sep 17 00:00:00 2001 From: ahmadtabikh Date: Mon, 2 Feb 2026 13:42:04 +0100 Subject: [PATCH 3/4] Draw unfolded spectrum with improved statistical errors --- Jets/JetSpectrum_DrawingFunctions.C | 30 +++++++++++ Jets/JetSpectrum_DrawingFunctions.h | 1 + Jets/JetSpectrum_Unfolding.C | 79 +++++++++++++++++++++++++++++ Jets/JetSpectrum_Unfolding.h | 2 + 4 files changed, 112 insertions(+) diff --git a/Jets/JetSpectrum_DrawingFunctions.C b/Jets/JetSpectrum_DrawingFunctions.C index f130429..563111f 100644 --- a/Jets/JetSpectrum_DrawingFunctions.C +++ b/Jets/JetSpectrum_DrawingFunctions.C @@ -1583,6 +1583,36 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete // } } +void Draw_Pt_spectrum_unfolded_RelativeUncertainty(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; diff --git a/Jets/JetSpectrum_DrawingFunctions.h b/Jets/JetSpectrum_DrawingFunctions.h index b7eee8c..71628b3 100644 --- a/Jets/JetSpectrum_DrawingFunctions.h +++ b/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_RelativeUncertainty(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); diff --git a/Jets/JetSpectrum_Unfolding.C b/Jets/JetSpectrum_Unfolding.C index 5cb3f82..39955e0 100644 --- a/Jets/JetSpectrum_Unfolding.C +++ b/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 \ No newline at end of file diff --git a/Jets/JetSpectrum_Unfolding.h b/Jets/JetSpectrum_Unfolding.h index 4867a44..28df392 100644 --- a/Jets/JetSpectrum_Unfolding.h +++ b/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 \ No newline at end of file From d875a9ccf3c978ea3a2f46f62915e7c98915dada Mon Sep 17 00:00:00 2001 From: ahmadtabikh Date: Mon, 2 Feb 2026 13:50:27 +0100 Subject: [PATCH 4/4] Draw unfolded spectrum with improved statistical errors --- Jets/JetSpectrum_DrawingFunctions.C | 2 +- Jets/JetSpectrum_DrawingFunctions.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Jets/JetSpectrum_DrawingFunctions.C b/Jets/JetSpectrum_DrawingFunctions.C index 563111f..0d046bc 100644 --- a/Jets/JetSpectrum_DrawingFunctions.C +++ b/Jets/JetSpectrum_DrawingFunctions.C @@ -1583,7 +1583,7 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete // } } -void Draw_Pt_spectrum_unfolded_RelativeUncertainty(int iDataset, int iRadius, int unfoldParameterInput, std::string options) { +void Draw_Pt_spectrum_unfolded_ImprovedStatErrors(int iDataset, int iRadius, int unfoldParameterInput, std::string options) { TH1D* measuredInput; if (!normGenAndMeasByNEvtsForUnfoldingInput) { diff --git a/Jets/JetSpectrum_DrawingFunctions.h b/Jets/JetSpectrum_DrawingFunctions.h index 71628b3..2a6f526 100644 --- a/Jets/JetSpectrum_DrawingFunctions.h +++ b/Jets/JetSpectrum_DrawingFunctions.h @@ -16,7 +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_RelativeUncertainty(int iDataset, int iRadius, int unfoldParameterInput, 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);