diff --git a/Jets/JetQC.C b/Jets/JetQC.C index 40e20ec..b652018 100644 --- a/Jets/JetQC.C +++ b/Jets/JetQC.C @@ -20,8 +20,8 @@ // #include "RooUnfoldBinByBin.h" //My Libraries -#include "./JetQC_settings.h" -#include "./JetQC_inputs.h" +#include "./JetQC_settings_template.h" +#include "./JetQC_inputs_template.h" #include "../Settings/AxisTitles.h" #include "../Settings/GlobalSettings.h" #include "../Utilities/AnalysisUtilities.h" @@ -50,80 +50,80 @@ void LoadLibs(); //////////// QC plot functions // Radius comparison -void Draw_Pt_RadiusComparison(int iDataset, float* etaRange); -void Draw_Eta_RadiusComparison(int iDataset, float* PtRange); -void Draw_Phi_RadiusComparison(int iDataset, float* PtRange); -void Draw_jetNTracks_RadiusComparison_withPtRange(int iDataset, float* PtRange); -void Draw_JetArea_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); -void Draw_JetNTracks_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); -void Draw_Pt_ratio_etaNeg_etaPos_RadiusComparison(int iDataset, float* etaRange); -void Draw_JetPhi_vs_JetEta_RadiusComparison(int iDataset); -void Draw_JetEta_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); -void Draw_JetPhi_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); -void Draw_JetTRDratio_vs_JetEta(int iDataset); -void Draw_JetTRDcount_vs_JetEta(int iDataset); -void Draw_Pt_ratio_etaNeg_etaPos_TRDonly_vs_noTRD(int iDataset, float* etaRange); - -void Draw_Pt_RadiusComparison_mcp(int iDataset, float* etaRange); +// void Draw_Pt_RadiusComparison(int iDataset, float* etaRange); +// void Draw_Eta_RadiusComparison(int iDataset, float* PtRange); +// void Draw_Phi_RadiusComparison(int iDataset, float* PtRange); +// void Draw_jetNTracks_RadiusComparison_withPtRange(int iDataset, float* PtRange); +// void Draw_JetArea_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); +// void Draw_JetNTracks_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); +// void Draw_Pt_ratio_etaNeg_etaPos_RadiusComparison(int iDataset, float* etaRange); +// void Draw_JetPhi_vs_JetEta_RadiusComparison(int iDataset); +// void Draw_JetEta_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); +// void Draw_JetPhi_vs_JetPt_RadiusComparison(int iDataset, float* PtRange); +// void Draw_JetTRDratio_vs_JetEta(int iDataset); +// void Draw_JetTRDcount_vs_JetEta(int iDataset); +// void Draw_Pt_ratio_etaNeg_etaPos_TRDonly_vs_noTRD(int iDataset, float* etaRange); + +// void Draw_Pt_RadiusComparison_mcp(int iDataset, float* etaRange); // Dataset comparison void Draw_Pt_DatasetComparison(float* etaRange, std::string options, float jetRadiusForJetFinderWorkflow); -void Draw_Eta_DatasetComparison(float jetRadius, float* PtRange, std::string options); -void Draw_Phi_DatasetComparison(float jetRadius, float* PtRange, std::string options); -void Draw_Pt_ratio_etaNeg_etaPos_DatasetComparison(float jetRadius, float* etaRange); -void Draw_Area_PtIntegrated_DatasetComparison(float jetRadius, float* PtRange); -void Draw_Rho_vs_Centrality_DatasetComp(); -void Draw_Rho_vs_SelectedMultiplicity_DatasetComp(); -void Draw_Rho_vs_LeadingJetPt_DatasetComp(); -void Draw_BkgFluctuations_vs_Centrality_DatasetComp(std::array, 2> drawnWindow); -void Draw_SelectedMultiplicity_vs_Centrality_DatasetComp(); -void Draw_Rho_vs_SelectedMultiplicity_DatasetCompRatio(); -void Draw_Rho_vs_SelectedMultiplicity_DatasetComp_withCutDemarcation(); -void Draw_Ncoll_vs_centrality(std::string options); -void Draw_Constituent_Pt_DatasetComparison(float* jetPtRange, float jetRadiusForJetFinderWorkflow); +// void Draw_Eta_DatasetComparison(float jetRadius, float* PtRange, std::string options); +// void Draw_Phi_DatasetComparison(float jetRadius, float* PtRange, std::string options); +// void Draw_Pt_ratio_etaNeg_etaPos_DatasetComparison(float jetRadius, float* etaRange); +// void Draw_Area_PtIntegrated_DatasetComparison(float jetRadius, float* PtRange); +// void Draw_Rho_vs_Centrality_DatasetComp(); +// void Draw_Rho_vs_SelectedMultiplicity_DatasetComp(); +// void Draw_Rho_vs_LeadingJetPt_DatasetComp(); +// void Draw_BkgFluctuations_vs_Centrality_DatasetComp(std::array, 2> drawnWindow); +// void Draw_SelectedMultiplicity_vs_Centrality_DatasetComp(); +// void Draw_Rho_vs_SelectedMultiplicity_DatasetCompRatio(); +// void Draw_Rho_vs_SelectedMultiplicity_DatasetComp_withCutDemarcation(); +// void Draw_Ncoll_vs_centrality(std::string options); +// void Draw_Constituent_Pt_DatasetComparison(float* jetPtRange, float jetRadiusForJetFinderWorkflow); -// Rebin comparison -void Draw_Area_PtIntegrated_BinningComparison(int iDataset, float jetRadius, float* PtRange); +// // Rebin comparison +// void Draw_Area_PtIntegrated_BinningComparison(int iDataset, float jetRadius, float* PtRange); -// Centrality comparison -void Draw_Pt_CentralityComparison(float jetRadius, int iDataset); -void Draw_BkgFluctuations_CentralityProjection(int iDataset, std::array, 2> drawnWindowZoom, std::string options); -void Draw_Eta_CentralityComparison(float jetRadius, int iDataset); -void Draw_Phi_CentralityComparison(float jetRadius, int iDataset); -void Draw_BkgFluctuations_CentralityProjection_withFit_CentralityComp(int iDataset, std::array, 2> drawnWindowZoom); -void Draw_BkgFluctuations_CentralityProjection_withFit_DatasetComp(float* centRange, std::array, 2> drawnWindowZoom); -void Draw_BkgFluctuations_CentralityProjection_withFit_MethodComp(float* centRange, int iDataset, std::array, 2> drawnWindowZoom); -void Draw_Rho_CentralityProjection(int iDataset, std::string options); -void Draw_Rho_CentralityProjection_DatasetComp(float* centRange, std::string options); +// // Centrality comparison +// void Draw_Pt_CentralityComparison(float jetRadius, int iDataset); +// void Draw_BkgFluctuations_CentralityProjection(int iDataset, std::array, 2> drawnWindowZoom, std::string options); +// void Draw_Eta_CentralityComparison(float jetRadius, int iDataset); +// void Draw_Phi_CentralityComparison(float jetRadius, int iDataset); +// void Draw_BkgFluctuations_CentralityProjection_withFit_CentralityComp(int iDataset, std::array, 2> drawnWindowZoom); +// void Draw_BkgFluctuations_CentralityProjection_withFit_DatasetComp(float* centRange, std::array, 2> drawnWindowZoom); +// void Draw_BkgFluctuations_CentralityProjection_withFit_MethodComp(float* centRange, int iDataset, std::array, 2> drawnWindowZoom); +// void Draw_Rho_CentralityProjection(int iDataset, std::string options); +// void Draw_Rho_CentralityProjection_DatasetComp(float* centRange, std::string options); -void Draw_RhoMean_asFunctionOf_Centrality(int iDataset, std::string options); -// NTracks comp -void Draw_Rho_withFit_NTracksProjection(int iDataset); +// void Draw_RhoMean_asFunctionOf_Centrality(int iDataset, std::string options); +// // NTracks comp +// void Draw_Rho_withFit_NTracksProjection(int iDataset); -// Pt cut comparison -void Draw_Eta_PtCutComparison(float jetRadius, int iDataset, float* PtCuts, int nPtCut, std::string options); -void Draw_Phi_PtCutComparison(float jetRadius, int iDataset, float* PtCuts, int nPtCut, std::string options); +// // Pt cut comparison +// void Draw_Eta_PtCutComparison(float jetRadius, int iDataset, float* PtCuts, int nPtCut, std::string options); +// void Draw_Phi_PtCutComparison(float jetRadius, int iDataset, float* PtCuts, int nPtCut, std::string options); -// Run 2 vs Run 3 comparison -void Draw_Pt_Run2Run3Comparison_0010Cent_R040(int iDataset, std::string options); +// // Run 2 vs Run 3 comparison +// void Draw_Pt_Run2Run3Comparison_0010Cent_R040(int iDataset, std::string options); -// MC jet pt resolution -void Draw_jet_resolution_MC_PtRangeComparison(int iDataset, double jetRadius, std::string options); -void Draw_Efficiency_Pt_DatasetComparison(float jetRadius, float* etaRange); -void Draw_Purity_Pt_DatasetComparison(float jetRadius, float* etaRange); +// // MC jet pt resolution +// void Draw_jet_resolution_MC_PtRangeComparison(int iDataset, double jetRadius, std::string options); +// void Draw_Efficiency_Pt_DatasetComparison(float jetRadius, float* etaRange); +// void Draw_Purity_Pt_DatasetComparison(float jetRadius, float* etaRange); -void Count_Nevents_perCentClass(int iDataset, std::string options); +// void Count_Nevents_perCentClass(int iDataset, std::string options); -void Draw_PtPeakPosition_vs_leadTrackCut(float* etaRange, float jetRadiusForJetFinderWorkflow); -void Draw_PtLeadCutStudy_PtOfRatio1(float* etaRange, std::string options, float jetRadiusForJetFinderWorkflow); +// void Draw_PtPeakPosition_vs_leadTrackCut(float* etaRange, float jetRadiusForJetFinderWorkflow); +// void Draw_PtLeadCutStudy_PtOfRatio1(float* etaRange, std::string options, float jetRadiusForJetFinderWorkflow); -void Draw_Pt_PbPbToPPComparison_HARDCODED(float jetRadius, float* etaRange, std::string options); +// void Draw_Pt_PbPbToPPComparison_HARDCODED(float jetRadius, float* etaRange, std::string options); -void Count_Jets_DatasetComparison(float jetRadius, float* ptRange, float* etaRange, std::string options); +// void Count_Jets_DatasetComparison(float jetRadius, float* ptRange, float* etaRange, std::string options); -void Draw_Pt_DatasetComparison_withRun2RitsuyaHardcoded(float* etaRange, std::string options, float jetRadiusForJetFinderWorkflow); +// void Draw_Pt_DatasetComparison_withRun2RitsuyaHardcoded(float* etaRange, std::string options, float jetRadiusForJetFinderWorkflow); ///////////////////////////////////////////////////// ///////////////////// Main Macro //////////////////// @@ -142,11 +142,11 @@ void JetQC() { // int iDataset = 0; - float etaRangeSym[2] = {-0.7, 0.7}; + float etaRangeSym[2] = {-0.5, 0.5}; float etaRangeNeg[2] = {-0.5, 0}; float etaRangePos[2] = {0, 0.5}; - float jetRadiusForDataComp = 0.2; + float jetRadiusForDataComp = 0.4; float jetR02 = 0.2; float jetR06 = 0.6; @@ -166,7 +166,7 @@ void JetQC() { // Draw_PtPeakPosition_vs_leadTrackCut(etaRangeSym, jetRadiusForDataComp); // Draw_PtLeadCutStudy_PtOfRatio1(etaRangeSym, "", jetRadiusForDataComp); - // Draw_Pt_DatasetComparison(etaRangeSym, "normEvents", jetRadiusForDataComp); + Draw_Pt_DatasetComparison(etaRangeSym, "normEvents", jetRadiusForDataComp); // Draw_Pt_DatasetComparison_withRun2RitsuyaHardcoded(etaRangeSym, "normEvents", jetRadiusForDataComp); for(int iPtBin = 0; iPtBin < nPtBins; iPtBin++){ jetPtMinCut = jetPtMinCutArray[iPtBin]; @@ -184,7 +184,7 @@ void JetQC() { // Draw_Eta_DatasetComparison(jetRadiusForDataComp, ptRange, "normEvents"); // Draw_Eta_DatasetComparison(jetRadiusForDataComp, ptRange, "normEntries"); // Draw_Phi_DatasetComparison(jetRadiusForDataComp, ptRange, "normEvents"); - Draw_Phi_DatasetComparison(jetRadiusForDataComp, ptRange, "normEvents"); + // Draw_Phi_DatasetComparison(jetRadiusForDataComp, ptRange, "normEntries"); // Draw_Constituent_Pt_DatasetComparison(ptRange, jetRadiusForDataComp); @@ -417,6 +417,7 @@ void SetStyle(Bool_t graypalette) { ///////////////////////////////////////////////////////////////////////////// QC plot functions ///////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// /* void Draw_Pt_RadiusComparison(int iDataset, float* etaRange) { TH3D* H3D_jetRjetPtjetEta; @@ -748,7 +749,7 @@ void Draw_JetPhi_vs_JetEta_RadiusComparison(int iDataset) { // Draw_TH2_Histograms(H2D_jetPhijetEta, RadiusLegend, nRadius, textContext, pdfNamelogyz, texPtX, texJetArea, texCollisionDataInfo, drawnWindowAuto, th2ContoursNone, contourNumberNone, "logxlogz"); } - +// */ void Draw_Pt_DatasetComparison(float* etaRange, std::string options, float jetRadiusForJetFinderWorkflow = 0.2) { float jetRadius = jetRadiusForJetFinderWorkflow; // obsolete for new jet-spectra-charged as we don't do radii comparisons that often and so files will only have 1 radius @@ -858,22 +859,24 @@ void Draw_Pt_DatasetComparison(float* etaRange, std::string options, float jetRa // std::array, 2> legendPlacementAuto = {{{-999, -999}, {-999, -999}}}; // {{{x1, y1}, {x2, y2}}} Draw_TH1_Histograms(H1D_jetPt_rebinned, DatasetsNames, nDatasets, textContext, pdfName, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texJetPtYield_EventNorm, texCollisionDataInfo, drawnWindow, legendPlacement, contextPlacementAuto, "logy"+histDatasetComparisonStructure); - if (divideSuccess == true) { - if (histDatasetComparisonStructure.find("twoByTwoDatasetPairs") != std::string::npos) { - Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatio, texCollisionDataInfo, drawnWindowAuto, legendPlacementRatio, contextPlacementAuto); - Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio_zoom, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatio, texCollisionDataInfo, drawnWindowZoomTwoByTwo, legendPlacementRatio, contextPlacementAuto); - } else { - Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_ratio, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatioDatasets, texCollisionDataInfo, drawnWindow, legendPlacementRatio, contextPlacementAuto, "noMarkerFirst"+histDatasetComparisonStructure); - Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_ratio_zoom, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatioDatasets, texCollisionDataInfo, drawnWindowZoom, legendPlacementAuto, contextPlacementAuto, "noMarkerFirst"); - } - } - else { - cout << "Divide failed in Draw_Pt_DatasetComparison" << endl; - } + // if (divideSuccess == true) { + // if (histDatasetComparisonStructure.find("twoByTwoDatasetPairs") != std::string::npos) { + // Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatio, texCollisionDataInfo, drawnWindowAuto, legendPlacementRatio, contextPlacementAuto); + // Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio_zoom, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatio, texCollisionDataInfo, drawnWindowZoomTwoByTwo, legendPlacementRatio, contextPlacementAuto); + // } else { + // Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_ratio, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatioDatasets, texCollisionDataInfo, drawnWindow, legendPlacementRatio, contextPlacementAuto, "noMarkerFirst"+histDatasetComparisonStructure); + // Draw_TH1_Histograms(H1D_jetPt_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_ratio_zoom, iJetFinderQaType == 0 ? texPtJetRawX : texPtJetBkgCorrX, texRatioDatasets, texCollisionDataInfo, drawnWindowZoom, legendPlacementAuto, contextPlacementAuto, "noMarkerFirst"); + // } + // } + // else { + // cout << "Divide failed in Draw_Pt_DatasetComparison" << endl; + // } // Draw_TH1_Histograms_ratioInSameCanvas(H1D_jetPt_rebinned, DatasetsNames, H1D_jetPt_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_testRatioInSameWindow, texPtX, texJetPtYield_EventNorm, texCollisionDataInfo, drawnWindow, drawnWindowAuto, legendPlacement, legendPlacementRatio, contextPlacementAuto, "ratioZoomToOneExtraExtra"+histDatasetComparisonStructure); } + + void Draw_Eta_DatasetComparison(float jetRadius, float* PtRange, std::string options) { TH3D* H3D_jetRjetPtjetEta[nDatasets]; @@ -983,8 +986,8 @@ void Draw_Eta_DatasetComparison(float jetRadius, float* PtRange, std::string opt if (divideSuccess == true) { if (histDatasetComparisonStructure.find("twoByTwoDatasetPairs") != std::string::npos) { - Draw_TH1_Histograms(H1D_jetEta_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio, texEtaX, texRatio, texCollisionDataInfo, drawnWindowAuto, legendPlacementRatio, contextPlacementAuto); - Draw_TH1_Histograms(H1D_jetEta_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio_zoom, texEtaX, texRatio, texCollisionDataInfo, drawnWindowZoomTwoByTwo, legendPlacementRatio, contextPlacementAuto); + // Draw_TH1_Histograms(H1D_jetEta_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio, texEtaX, texRatio, texCollisionDataInfo, drawnWindowAuto, legendPlacementRatio, contextPlacementAuto); + // Draw_TH1_Histograms(H1D_jetEta_rebinned_ratios, DatasetsNamesPairRatio, nHistPairRatio, textContext, pdfName_ratio_zoom, texEtaX, texRatio, texCollisionDataInfo, drawnWindowZoomTwoByTwo, legendPlacementRatio, contextPlacementAuto); } else { @@ -1108,7 +1111,7 @@ void Draw_Phi_DatasetComparison(float jetRadius, float* PtRange, std::string opt cout << "Divide failed in Draw_Phi_DatasetComparison" << endl; } } - +/* void Draw_Pt_ratio_etaNeg_etaPos_RadiusComparison(int iDataset, float* etaRange) { TH3D* H3D_jetRjetPtjetEta; @@ -4262,3 +4265,4 @@ void Draw_Pt_DatasetComparison_withRun2RitsuyaHardcoded(float* etaRange, std::st } +*/ \ No newline at end of file diff --git a/Jets/JetQC_inputs_template.h b/Jets/JetQC_inputs_template.h index e705e57..97c81d1 100644 --- a/Jets/JetQC_inputs_template.h +++ b/Jets/JetQC_inputs_template.h @@ -1965,15 +1965,15 @@ TFile* file_AliAnalysis; //////// -------- LHC25b6 - pp sim anchored to PbPb localTest //////// -TString* texCollisionDataInfo = new TString("pp #sqrt{#it{s}} = 13.6 TeV"); +TString* texCollisionDataInfo = new TString("pp #sqrt{#it{s}} = 5.36 TeV"); const TString* texDatasetsComparisonType = new TString(""); const TString* texDatasetsComparisonCommonDenominator = new TString(""); const int nDatasets = 2; -const TString Datasets[nDatasets] = {"pp13TeV_sim_MB", "pp13TeV_sim_jetjet"}; +const TString Datasets[nDatasets] = {"LHC26a6_train615296_Unf", "LHC25b4ab6_train600389"}; // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"}; -const TString DatasetsNames[nDatasets] = {"MB", "jet-jet"}; -TFile* file_O2Analysis_list[nDatasets] = {new TFile("Datasets/"+Datasets[0]+"/AnalysisResults.root"), - new TFile("Datasets/"+Datasets[1]+"/AnalysisResults.root") +const TString DatasetsNames[nDatasets] = {"JJ LHC26a6", "MB LHC25b4ab6"}; +TFile* file_O2Analysis_list[nDatasets] = {new TFile("../Datasets/"+Datasets[0]+"/AnalysisResults.root"), + new TFile("../Datasets/"+Datasets[1]+"/AnalysisResults.root") }; const TString analysisWorkflow[nDatasets] = {"jet-spectra-charged", "jet-spectra-charged", @@ -1982,6 +1982,6 @@ const TString analysisWorkflow[nDatasets] = {"jet-spectra-charged", "jet-spectra const TString wagonId[nDatasets] = {"", "" }; -const bool isDatasetWeighted[nDatasets] = {false, true}; +const bool isDatasetWeighted[nDatasets] = {true, false}; const std::string histDatasetComparisonStructure = ""; const bool datasetsAreSubsetsofId0 = false; \ No newline at end of file diff --git a/Jets/JetQC_settings_template.h b/Jets/JetQC_settings_template.h index bec1a7d..fd9fdd5 100644 --- a/Jets/JetQC_settings_template.h +++ b/Jets/JetQC_settings_template.h @@ -15,8 +15,8 @@ const TString jetLevel[nJetLevel] = {"data", "mcd", "mcp"}; // float arrayRadius[nRadius] = {0.2, 0.4, 0.6}; // const float areaDisplayMax[nRadius] = {0.5, 1, 1.5}; const int nRadius = 1; -const TString RadiusLegend[nRadius] = {"R = 0.2"}; -float arrayRadius[nRadius] = {0.2}; +const TString RadiusLegend[nRadius] = {"R = 0.4"}; +float arrayRadius[nRadius] = {0.4}; const float areaDisplayMax[nRadius] = {0.5}; // const int nRadius = 9; // const TString RadiusLegend[nRadius] = {"R = 0.2", "R = 0.25", "R = 0.3", "R = 0.35", "R = 0.4", "R = 0.45", "R = 0.5", "R = 0.55", "R = 0.6"}; @@ -28,13 +28,13 @@ TString methodRandomConeHistLegend[nMethodRC] = {"Random Cones (RC)", "RC w/o le // Choice of jet type (charged, neutral, full) and level (data, detector level, particle level) const int iJetType = 0; -const int iJetLevel = 0; +const int iJetLevel = 1; // Choice of jet QA type (uncorrected jets, background corrected jet (rho area version), event wise constituent subtraction -const int iJetFinderQaType = 1; +const int iJetFinderQaType = 0; // Choice of Random Cone method: -const int iMethodRandomCone = 1; +const int iMethodRandomCone = 0; // Default window for random cone: std::array, 2> drawnWindowRCdefault = {{{-30, 60}, {5E-7, 20}}}; // {{xmin, xmax}, {ymin, ymax}} diff --git a/Jets/JetSpectrum_DrawingFunctions.C b/Jets/JetSpectrum_DrawingFunctions.C index f7cb231..fd8934a 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" @@ -695,7 +699,16 @@ void Draw_Pt_spectrum_unfolded_singleDataset(int iDataset, int iRadius, int unfo unfoldParameter = Get_Pt_spectrum_unfolded(H1D_jetPt_unfolded, measuredInput, iDataset, iRadius, unfoldParameterInput, options).first; // TH1D* H1D_jetPt_unfolded2 = (TH1D*)H1D_jetPt_unfolded->Clone(H1D_jetPt_unfolded->GetName()+(TString)"H1D_jetPt_unfolded2"); - + //////////////////////////////////////////////////////////// + + if (writeOutputRootFile) { + cout << "######################### FILE CREATED IN PRINCIPLE##################" << endl; + TFile* outFile = new TFile("output.root", "UPDATE"); // "UPDATE" Open existing or create if missing / "RECREATE" Always deletes file and creates new one + H1D_jetPt_unfolded->Write("MB_LHC25b4ab6_train600389_ppref_unfolded"); // MB_LHC25b4b5_train533385_Unf_ppref_unfolded + outFile->Close(); + cout << "######################### HISTO SAVED IN PRINCIPLE ##################" << endl; + } + //////////////////////////////////////////////////////////// cout << "comparison with raw measured" << endl; if (!useFineBinningTest) { Get_Pt_spectrum_bkgCorrected_genBinning(H1D_jetPt_measured_genBinning, iDataset, iRadius, options); @@ -1112,6 +1125,8 @@ void Draw_Pt_spectrum_unfolded_singleDataset(int iDataset, int iRadius, int unfo } } } +// declare this function here. It is used in the funciton below "Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset" +void DrawRatioWithOffset(TH1D* histList[], int nUnfoldIteration, const TString& yAxisTitle,const TString& canvasName, int unfoldIterationMax, int step, double yMin, double yMax); void Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset(int iDataset, int iRadius, int unfoldIterationMin, int unfoldIterationMax, int step, std::string options) { @@ -1268,6 +1283,14 @@ void Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset(int iDataset, in TString* pdfName_ratio_refoldedComp_zoom = new TString("jet_"+jetType[iJetType]+"_"+jetLevel[iJetLevel]+"_"+partialUniqueSpecifier+"_Pt_unfolded_"+unfoldingInfo+"_ratioRefoldedUnfolded_zoom"); Draw_TH1_Histograms(H1D_jetPt_ratio_measuredRefolded, unfoldingIterationLegend, nUnfoldIteration, textContext, pdfName_ratio_refoldedComp_zoom, texPtX, texRatioRefoldedMeasured, texCollisionDataInfo, drawnWindowUnfoldedMeasurement, legendPlacementAuto, contextPlacementAuto, "zoomToOneLarge,ratioLine,zoomToOneMedium2"); } + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////ratio refolded/measured with offset differetn k//////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + DrawRatioWithOffset(H1D_jetPt_ratio_measuredRefolded, nUnfoldIteration, "Refolded / Measured", "ratio_RefoldMeasure_different_k_withOffset", unfoldIterationMax, step, 0.8, 1.3); + DrawRatioWithOffset(H1D_jetPt_ratio_mcp, nUnfoldIteration, "unfolded / mcp", "ratio_Unfolded_mcp_different_k_withOffset", unfoldIterationMax, step, 0.7, 1.4); // last two numbers are y min and max + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + } void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParameterInput, std::string options) { @@ -1379,8 +1402,10 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete divideSuccessDatasets[iDataset] = H1D_jetPt_unfolded_ratio_datasets[iDataset]->Divide(H1D_jetPt_unfolded[0]); // Creating to-be-plotted histograms datasetNameSpecifier[iDataset] = "_"+DatasetsNames[iDataset]+Form("%.1d",iDataset); + // #################### + - + // #################### cout << "comparison with raw measured" << endl; if (!useFineBinningTest) { Get_Pt_spectrum_bkgCorrected_genBinning(H1D_jetPt_measured_genBinning[iDataset], iDataset, iRadius, options); @@ -1569,6 +1594,153 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete // } } +void Draw_Pt_spectrum_unfolded_ImprovedStatErrors(int iDataset, int iRadius, int unfoldParameterInput, std::string options) { + +TH1D* measuredInput; +if (!normGenAndMeasByNEvtsForUnfoldingInput) { + Get_Pt_spectrum_bkgCorrected_recBinning_preWidthScalingAtEndAndEvtNorm(measuredInput, iDataset, iRadius, options); + if (useFineBinningTest) { + Get_Pt_spectrum_bkgCorrected_fineBinning_preWidthScalingAtEndAndEvtNorm(measuredInput, iDataset, iRadius, options); + } +} else{ + Get_Pt_spectrum_bkgCorrected_recBinning_preWidthScalingAtEnd(measuredInput, iDataset, iRadius, options); + if (useFineBinningTest) { + Get_Pt_spectrum_bkgCorrected_fineBinning_preWidthScalingAtEnd(measuredInput, iDataset, iRadius, options); + } +} +TH1D* H1D_jetPt_unfolded_withImprovedErrors; +Get_Pt_spectrum_unfolded_ImprovedStatisticalErrors(H1D_jetPt_unfolded_withImprovedErrors, measuredInput, iDataset, iRadius, unfoldParameterInput, options); + +TString pdfname = "jet_Pt_spectrum_unfolded_RelativeUncertainty_Dataset"+DatasetsNames[iDataset]+"_R"+Form("%.1f", arrayRadius[iRadius])+"_k"+Form("%i", unfoldParameterInput); +// TString textContext = "Unfolded improved errors"; +TString textContext(contextCustomOneField(*texDatasetsComparisonCommonDenominator, "")); +// Test : +// TCanvas* c_relUnc = new TCanvas(pdfname, pdfname, 800, 800); +// H1D_jetPt_unfolded_withImprovedErrors->Draw(); +// c_relUnc->SetLogy(); + +// error with Draw_TH1_Histogram!!!? +// Draw_TH1_Histogram(H1D_jetPt_unfolded_withImprovedErrors, textContext, pdfname, texPtJetRec, texJet_d2Ndptdeta_EventNorm, texCollisionDataInfo, drawnWindowUnfoldedMeasurement, legendPlacementAuto, contextPlacementAuto, "logy"); + +} + +void DrawRatioWithOffset(TH1D* histList[], int nUnfoldIteration, const TString& yAxisTitle,const TString& canvasName, int unfoldIterationMax, int step, double yMin, double yMax){ + // DrawRatioWithOffset(..., -1, -1); + TString canvasNameFull = canvasName + "_" + unfoldingMethod; + TCanvas* c = new TCanvas(canvasNameFull, canvasNameFull, 800, 800); + TMultiGraph* mg = new TMultiGraph(); + + int customColors[] = { + kBlack, kRed + 1, kBlue + 1, kGreen + 2, kOrange + 7, kViolet + 1 + }; + int nColors = sizeof(customColors) / sizeof(customColors[0]); + + //int nUnfoldIteration = histList.size(); // infer number of iterations from input list + + for (int i = 0; i < nUnfoldIteration; ++i) { + TH1D* h = histList[i]; + int nBins = h->GetNbinsX(); + + std::vector x_vals, y_vals, ex_vals, ey_vals; + + for (int bin = 1; bin <= nBins; ++bin) { + // Replace ptBinsJetsRec[iRadius] with your own binning logic if needed: + double binLowEdge = h->GetXaxis()->GetBinLowEdge(bin); + double binUpEdge = h->GetXaxis()->GetBinUpEdge(bin); + double binWidth = binUpEdge - binLowEdge; + double offset = (i - nUnfoldIteration / 2.0) * 0.065 * binWidth; + + double x = h->GetBinCenter(bin) + offset; + double y = h->GetBinContent(bin); + double ex = 0; + double ey = h->GetBinError(bin); + + x_vals.push_back(x); + y_vals.push_back(y); + ex_vals.push_back(ex); + ey_vals.push_back(ey); + } + + TGraphErrors* gr = new TGraphErrors(nBins, &x_vals[0], &y_vals[0], &ex_vals[0], &ey_vals[0]); + int color = customColors[i % nColors]; + gr->SetMarkerStyle(20 + i); + gr->SetMarkerColor(color); + gr->SetLineColor(color); + gr->SetMarkerSize(0.8); + gr->SetTitle(Form("k_{unfold} = %d", unfoldIterationMax-step*i)); // Adjust if you have different logic + + mg->Add(gr, "P"); + } + + mg->Draw("A"); + mg->GetXaxis()->SetLimits(5.0, 140.0); + if (yMin < yMax) { + mg->GetYaxis()->SetRangeUser(yMin, yMax); + } + mg->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + mg->GetYaxis()->SetTitle(yAxisTitle); + + c->BuildLegend(); + c->Update(); + + // Draw horizontal reference line at y=1 + double xmin = mg->GetXaxis()->GetXmin(); + double xmax = mg->GetXaxis()->GetXmax(); + TLine* line = new TLine(xmin, 1.0, xmax, 1.0); + line->SetLineStyle(2); + line->SetLineColor(kGray + 2); + line->Draw("same"); + + c->Update(); + // Auto-save + c->SaveAs(canvasNameFull + ".pdf"); + c->SaveAs(canvasNameFull + ".png"); +} + +void MakeRatio() +{ + // Open file in UPDATE mode (so we can write result) + TFile* MB = TFile::Open("../20260226_Unf600389_LHC25b4ab6/output.root", "READ"); + TFile* JJ = TFile::Open("output.root", "READ"); + + if (!MB || MB->IsZombie()) { + std::cerr << "Error: cannot open MB file!" << std::endl; + return; + } + if (!JJ || JJ->IsZombie()) { + std::cerr << "Error: cannot open JJ file!" << std::endl; + return; + } + // JJ_LHC26a6_train615296_Unf_ppref_unfolded // MB_LHC25b4b5_train533385_Unf_ppref_unfolded + // Retrieve histograms + TH1D* h1 = (TH1D*)MB->Get("MB_LHC25b4ab6_train600389_ppref_unfolded"); + TH1D* h2 = (TH1D*)JJ->Get("JJ_LHC26a6_train615296_Unf_ppref_unfolded"); + + if (!h1 || !h2) { + std::cerr << "Error: one of the histograms not found!" << std::endl; + MB->Close(); + JJ->Close(); + return; + } + + // Clone first histogram to store ratio + TH1D* hRatio = (TH1D*)h1->Clone("hRatio"); + + // Perform division + hRatio->Divide(h2); + TString* pdfName = new TString("Ratio of MB to JJ"); + TString textContext(contextCustomOneField(*texDatasetsComparisonCommonDenominator, "Ratio_of_DataUnfold_w_MB_LHC25b4ab6_to_JJ_LHC26a6")); + TString* YLabel = new TString("Unf MB / Unf JJ"); + + // Write ratio to file (overwrite if exists) + // hRatio->Write("", TObject::kOverwrite); + Draw_TH1_Histogram(hRatio, textContext, pdfName, texPtJetRec, YLabel, texCollisionDataInfo, drawnWindowUnfoldedMeasurement, legendPlacementAuto, contextPlacementAuto, "ratioLine"); + // file->Close(); + + std::cout << "Ratio successfully created and saved." << std::endl; +} + + // rename refoldedUnfolded as closure test? // and try and spend 15 min to clean hist names for the spectrum analysis @@ -1578,4 +1750,6 @@ void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParamete // hMcEfficiency_vsPt->Divide(hMcSignalCount_vsPt,TrueV0PtSpectrum_AnalysisBins, 1., 1., "b"); // option b for binomial because efficiency: https://twiki.cern.ch/twiki/bin/view/ALICE/PWGLFPAGSTRANGENESSEfficiency + + #endif \ No newline at end of file diff --git a/Jets/JetSpectrum_DrawingFunctions.h b/Jets/JetSpectrum_DrawingFunctions.h index b7eee8c..4b65d15 100644 --- a/Jets/JetSpectrum_DrawingFunctions.h +++ b/Jets/JetSpectrum_DrawingFunctions.h @@ -16,9 +16,11 @@ void Draw_Pt_spectrum_unfolded_singleDataset(int iDataset, int iRadius, int unfo void Draw_Pt_spectrum_unfolded_datasetComparison(int iRadius, int unfoldParameterInput, std::string options); void Draw_Pt_TestSpectrum_unfolded(int iDataset, int iRadius, std::string options); void Draw_Pt_spectrum_unfolded_parameterVariation_singleDataset(int iDataset, int iRadius, int unfoldIterationMin, int unfoldIterationMax, int step, std::string options); - +void Draw_Pt_spectrum_unfolded_ImprovedStatErrors(int iDataset, int iRadius, int unfoldParameterInput, std::string options); +void MakeRatio(); void Draw_Pt_efficiency_jets(int iRadius, std::string options); void Draw_kinematicEfficiency(int iRadius, std::string options); void Draw_FakeRatio(int iRadius, std::string options); + #endif \ 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_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 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..76ef682 100644 --- a/Jets/JetSpectrum_systematics.C +++ b/Jets/JetSpectrum_systematics.C @@ -45,6 +45,7 @@ #include #include #include +#include using namespace std; // Misc utilities @@ -55,6 +56,10 @@ 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); +void Draw_Systematics_SecondaryContamination(int iDataset, int iRadius, int unfoldParameterInput, std::string options); + ///////////////////////////////////////////////////// ///////////////////// Main Macro //////////////////// @@ -72,26 +77,173 @@ 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); + + //######################################################### Secondary tracks Systematics ##################################################### + char optionsAnalysis[100] = ""; + snprintf(optionsAnalysis, sizeof(optionsAnalysis), "%s,%s,%s", unfoldingPrior, unfoldingMethod); + int unfoldParameterInput = 8; + Draw_Systematics_SecondaryContamination(iDataset, iRadius, unfoldParameterInput, optionsAnalysis); - 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, 10}; - - Draw_Systematics_UnfoldMethod(iDataset, iRadius, unfoldingMethodList, unfoldParameterInputList, nUnfoldingMethods, optionsAnalysis_withoutUnfoldingMethod); } ///////////////////////////////////////////////////// /////////////////// 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, 0, 7, -30, 0); + // p0, p1, p2, p3, p4, p5, p6, p7 + + // fitFunctionInit->SetParameters(0.5, 7, -10, 0.5); + // 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"); @@ -233,8 +385,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 +408,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 +420,563 @@ 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 = 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, ""); + + +} + +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 + 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); + + TH1D* H1D_FitUnf_ratio = (TH1D*) hJetPtRebinned->Clone("hRatio"); + H1D_FitUnf_ratio->SetTitle("Ratio: Fit / Unfolded; p_{T}^{jet}; Ratio"); + H1D_FitUnf_ratio->Divide(H1D_jetPt_unfolded); + + TString* pdfName_fitUnfRatio = new TString("ratio fit histogram to unfolded spectrum"); TString textContext(""); + TString* yLabel = new TString("Fit / Unfolded"); + Draw_TH1_Histogram(H1D_FitUnf_ratio, textContext, pdfName_fitUnfRatio, texPtJetRec, yLabel, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); + + + // 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) + // }; - 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 + // 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 diff --git a/Settings/AxisTitles.h b/Settings/AxisTitles.h index a127fbc..32add9d 100644 --- a/Settings/AxisTitles.h +++ b/Settings/AxisTitles.h @@ -102,6 +102,7 @@ TString* texPtPeak = new TString("#it{p}_{T,ch jet}^{rec} - #it{A}_{jet} #it{#rh TString* texPtAtThreshold = new TString("#it{p}_{T,ch jet}^{rec} - #it{A}_{jet} #it{#rho} (GeV/#it{c}) at 95% threshold"); TString* texSystematicsPercent = new TString ("systematic error (%)"); +TString* texRelativeUncertainty = new TString ("relative uncertainty"); ////////////////////////// diff --git a/Tracks/TrackMcQC.C b/Tracks/TrackMcQC.C index 9b0a8eb..1103d1e 100644 --- a/Tracks/TrackMcQC.C +++ b/Tracks/TrackMcQC.C @@ -19,8 +19,8 @@ // #include "RooUnfoldBinByBin.h" //My Libraries -#include "./TrackMcQC_settings.h" -#include "./TrackMcQC_inputs.h" +#include "./TrackMcQC_settings_template.h" +#include "./TrackMcQC_inputs_template.h" #include "../Settings/AxisTitles.h" #include "../Settings/GlobalSettings.h" #include "../Utilities/AnalysisUtilities.h" @@ -109,16 +109,16 @@ void TrackMcQC() { float etaRange[2] = {-0.9, 0.9}; float ptRange[2] = {0.15, 100}; bool useLargeLegendWindow = false; - // Draw_Efficiency_Pt_DatasetComparison(etaRange,useSplit, useLargeLegendWindow); + Draw_Efficiency_Pt_DatasetComparison(etaRange,useSplit, useLargeLegendWindow); // Draw_Efficiency_Eta_DatasetComparison(ptRange,useSplit); // Draw_Efficiency_Phi_DatasetComparison(ptRange, etaRange, useSplit); // Draw_Efficiency_Pt_ratio_etaNeg_etaPos_DatasetComparison(etaRange, useSplit); // // Draw_Efficiency_Phi_DatasetComparison_finerPhi(ptRange1, etaRange); // only works with very specific datasets created locally - Draw_Purity_Pt_DatasetComparison(etaRange, useSplit); - Draw_Purity_Eta_DatasetComparison(etaRange, useSplit); - Draw_Purity_Phi_DatasetComparison(etaRange, useSplit); + // Draw_Purity_Pt_DatasetComparison(etaRange, useSplit); + // Draw_Purity_Eta_DatasetComparison(etaRange, useSplit); + // Draw_Purity_Phi_DatasetComparison(etaRange, useSplit); // Draw_Purity_Pt_ratio_etaNeg_etaPos_DatasetComparison(etaRange, useSplit); // int nPtRanges = 10; @@ -168,10 +168,12 @@ void TrackMcQC() { // // Draw_Phi_tracksReco_fromEffWorkflow_DatasetComparison(ptRange, etaRange, "primaries, secondaries, nonassociatedtrack, ratio, entriesNorm"); // Draw_Phi_tracksReco_fromEffWorkflow_DatasetComparison(ptRange, etaRange, "primaries, secondaries, nonassociatedtrack, ratio, evtNorm"); + // ##### WARNING: I changed the projection interval from bin = 0 instead of 1 ##### + // Draw_Pt_gen_DatasetComparison_H2CentVersion("primaries,secondaries, ratio, evtNorm"); // entriesNorm + // Draw_Eta_gen_DatasetComparison_H2CentVersion("primaries,secondaries, ratio, evtNorm"); // entriesNorm + // Draw_Phi_gen_DatasetComparison_H2CentVersion("primaries,secondaries, ratio, evtNorm"); + // Draw_Phi_gen_DatasetComparison_H2CentVersion("primaries,secondaries, ratio, entriesNorm"); - // Draw_Pt_gen_DatasetComparison_H2CentVersion("primaries,secondaries, ratio"); - // Draw_Eta_gen_DatasetComparison_H2CentVersion("primaries,secondaries, ratio"); - // Draw_Phi_gen_DatasetComparison_H2CentVersion("primaries,secondaries, ratio"); // Draw_PtResolution_Residuals("ptRes_vs_pt"); } ///////////////////////////////////////////////////// @@ -547,6 +549,12 @@ void Draw_Efficiency_Pt_DatasetComparison(float* etaRange, bool useSplit, bool u if (divideSuccess == true && divideSuccess_ptHigh == true) { Draw_TH1_Histograms(H1D_trackPt_efficiency_concatenated, DatasetsNames, nDatasets, textContext, pdfNameEntriesNorm, texPtMC, texTrackEfficiency, texCollisionDataInfo, drawnWindowLogEff, useLargeLegendWindow ? legendPlacementEfficiencyLarge : legendPlacementEfficiencyNarrow, contextPlacementAuto, "logx,efficiency,150MevLine"+histDatasetComparisonStructure); Draw_TH1_Histograms(H1D_trackPt_efficiency_concatenated, DatasetsNames, nDatasets, textContext, pdfNameEntriesNorm_zoom, texPtMC, texTrackEfficiency, texCollisionDataInfo, drawnWindowLogEff_zoom, useLargeLegendWindow ? legendPlacementEfficiencyLarge : legendPlacementEfficiencyNarrow, contextPlacementAuto, "logx,efficiency,150MevLine"+histDatasetComparisonStructure); + + // // relative uncertainty: percentage + // TH1D* H1D_trackPt_efficiency_concatenated_relUncert[nDatasets]; + // Compute_TH1D_RelativeUncertainty(H1D_trackPt_efficiency_concatenated, H1D_trackPt_efficiency_concatenated_relUncert, nDatasets); + // TString* pdfNameEntriesNorm_relUncert = new TString("track_Pt_efficiency_relativeUncertainty_"+dummyName[0]+"_@eta["+Form("%.1f", etaRange[0])+","+Form("%.1f", etaRange[1])+"]"); + // Draw_TH1_Histograms(H1D_trackPt_efficiency_concatenated_relUncert, DatasetsNames, nDatasets, textContext, pdfNameEntriesNorm_relUncert, texPtMC, texRelativeUncertainty, texCollisionDataInfo, drawnWindowLogEff_zoom, useLargeLegendWindow ? legendPlacementEfficiencyLarge : legendPlacementEfficiencyNarrow, contextPlacementAuto, "logx,relativeUncertainty,150MevLine"+histDatasetComparisonStructure); } else { cout << "Divide failed in Draw_Efficiency_Pt_DatasetComparison eff" << endl; @@ -2504,32 +2512,54 @@ void Draw_Pt_gen_DatasetComparison_H2CentVersion(std::string options) { H2D_centrality_track[iDataset] = (TH2D*)((TH2D*)file_O2Analysis_list[iDataset]->Get(analysisWorkflow[iDataset]+"/h2_centrality_particle_pt"))->Clone("Draw_Pt_DatasetComparison"+Datasets[iDataset]+DatasetsNames[iDataset]); - H1D_trackPt[iDataset] = (TH1D*)H2D_centrality_track[iDataset]->ProjectionY("trackPt_"+Datasets[iDataset]+DatasetsNames[iDataset], 1, H2D_centrality_track[iDataset]->GetNbinsX(), "e"); + H1D_trackPt[iDataset] = (TH1D*)H2D_centrality_track[iDataset]->ProjectionY("trackPt_"+Datasets[iDataset]+DatasetsNames[iDataset], 0, H2D_centrality_track[iDataset]->GetNbinsX(), "e"); H1D_trackPt_rebinned[iDataset] = (TH1D*)H1D_trackPt[iDataset]->Rebin(2.,"trackPt_rebinned_"+Datasets[iDataset]+DatasetsNames[iDataset]); // NormaliseYieldToNEntries(H1D_trackPt_rebinned[iDataset]); - if (isDatasetWeighted[iDataset]) { - Nevents = GetNEventsSelected_TrackEffWorkflow_gen_weighted(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); - } else { - Nevents = GetNEventsSelected_TrackEffWorkflow_gen(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); - } - NormaliseYieldToNEvents(H1D_trackPt_rebinned[iDataset], Nevents); + // if (isDatasetWeighted[iDataset]) { + // Nevents = GetNEventsSelected_TrackEffWorkflow_gen_weighted(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + // } else { + // Nevents = GetNEventsSelected_TrackEffWorkflow_gen(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + // } + // NormaliseYieldToNEvents(H1D_trackPt_rebinned[iDataset], Nevents); + if (options.find("evtNorm") != std::string::npos) { + if (isDatasetWeighted[iDataset]) { + Nevents = GetNEventsSelected_TrackEffWorkflow_gen_weighted(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + } else { + Nevents = GetNEventsSelected_TrackEffWorkflow_gen(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + } + NormaliseYieldToNEvents(H1D_trackPt_rebinned[iDataset], Nevents); + } + if (options.find("entriesNorm") != std::string::npos) { + NormaliseYieldToIntegral(H1D_trackPt_rebinned[iDataset]); + } H1D_trackPt_rebinned_ratios[iDataset] = (TH1D*)H1D_trackPt_rebinned[iDataset]->Clone("trackPt_rebinned_ratios"+Datasets[iDataset]+DatasetsNames[iDataset]); H1D_trackPt_rebinned_ratios[iDataset]->Reset("M"); divideSuccess = H1D_trackPt_rebinned_ratios[iDataset]->Divide(H1D_trackPt_rebinned[iDataset], H1D_trackPt_rebinned[0], 1., 1., datasetsAreSubsetsofId0 ? "b" : ""); } - TString* pdfName = new TString("track_Pt_gen_DataComp"); - TString* pdfName_ratio = new TString("track_Pt_gen_DataComp_ratio"); + TString* textYaxis; + TString pdfNameNorm; + if (options.find("evtNorm") != std::string::npos) { + textYaxis = texTrackPtYield_EventNorm; + pdfNameNorm = (TString)"_EventNorm"; + } + if (options.find("entriesNorm") != std::string::npos) { + textYaxis = texTrackPtYield_EntriesNorm; + pdfNameNorm = (TString)"_EntriesNorm"; + } + + TString* pdfName = new TString("track_Pt_gen_DataComp"+pdfNameNorm); + TString* pdfName_ratio = new TString("track_Pt_gen_DataComp"+pdfNameNorm+"_ratio"); TString textContext(contextTrackDatasetComp("")); - Draw_TH1_Histograms(H1D_trackPt_rebinned, DatasetsNames, nDatasets, textContext, pdfName, texPtMC, texTrackPtYield_EventNorm, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, "logx,logy"); + Draw_TH1_Histograms(H1D_trackPt_rebinned, DatasetsNames, nDatasets, textContext, pdfName, texPtMC, textYaxis, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, "logx,logy"); if (divideSuccess == true && options.find("ratio") != std::string::npos) { Draw_TH1_Histograms(H1D_trackPt_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_ratio, texPtMC, texRatioDatasets, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, "autoratio,noMarkerFirst,logx"); } @@ -2555,32 +2585,49 @@ void Draw_Eta_gen_DatasetComparison_H2CentVersion(std::string options) { H2D_centrality_track[iDataset] = (TH2D*)((TH2D*)file_O2Analysis_list[iDataset]->Get(analysisWorkflow[iDataset]+"/h2_centrality_particle_eta"))->Clone("Draw_Eta_DatasetComparison"+Datasets[iDataset]+DatasetsNames[iDataset]); - H1D_trackEta[iDataset] = (TH1D*)H2D_centrality_track[iDataset]->ProjectionY("trackEta_"+Datasets[iDataset]+DatasetsNames[iDataset], 1, H2D_centrality_track[iDataset]->GetNbinsX(), "e"); + H1D_trackEta[iDataset] = (TH1D*)H2D_centrality_track[iDataset]->ProjectionY("trackEta_"+Datasets[iDataset]+DatasetsNames[iDataset], 0, H2D_centrality_track[iDataset]->GetNbinsX(), "e"); H1D_trackEta_rebinned[iDataset] = (TH1D*)H1D_trackEta[iDataset]->Rebin(5.,"trackEta_rebinned"+Datasets[iDataset]+DatasetsNames[iDataset]); - if (isDatasetWeighted[iDataset]) { - Nevents = GetNEventsSelected_TrackEffWorkflow_gen_weighted(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); - } else { - Nevents = GetNEventsSelected_TrackEffWorkflow_gen(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + if (options.find("evtNorm") != std::string::npos) { + if (isDatasetWeighted[iDataset]) { + Nevents = GetNEventsSelected_TrackEffWorkflow_gen_weighted(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + } else { + Nevents = GetNEventsSelected_TrackEffWorkflow_gen(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + } + NormaliseYieldToNEvents(H1D_trackEta_rebinned[iDataset], Nevents); + } + if (options.find("entriesNorm") != std::string::npos) { + NormaliseYieldToIntegral(H1D_trackEta_rebinned[iDataset]); } - NormaliseYieldToNEvents(H1D_trackEta_rebinned[iDataset], Nevents); H1D_trackEta_rebinned_ratios[iDataset] = (TH1D*)H1D_trackEta_rebinned[iDataset]->Clone("trackEta_rebinned_ratios"+Datasets[iDataset]+DatasetsNames[iDataset]); H1D_trackEta_rebinned_ratios[iDataset]->Reset("M"); divideSuccess = H1D_trackEta_rebinned_ratios[iDataset]->Divide(H1D_trackEta_rebinned[iDataset], H1D_trackEta_rebinned[0], 1., 1., datasetsAreSubsetsofId0 ? "b" : ""); } - TString* pdfNameEventNorm = new TString("track_Eta_gen_DataComp_EventNorm"); - TString* pdfNameEventNorm_ratio = new TString("track_Eta_gen_DataComp_EventNorm_ratio"); + TString* textYaxis; + TString pdfNameNorm; + if (options.find("evtNorm") != std::string::npos) { + textYaxis = texTrackPtYield_EventNorm; + pdfNameNorm = (TString)"_EventNorm"; + } + if (options.find("entriesNorm") != std::string::npos) { + textYaxis = texTrackPtYield_EntriesNorm; + pdfNameNorm = (TString)"_EntriesNorm"; + } + + TString* pdfName = new TString("track_Eta_gen_DataComp"+pdfNameNorm); + TString* pdfName_ratio = new TString("track_Eta_gen_DataComp"+pdfNameNorm+"_ratio"); + TString textContext(contextTrackDatasetComp("")); - Draw_TH1_Histograms(H1D_trackEta_rebinned, DatasetsNames, nDatasets, textContext, pdfNameEventNorm, texEtaMC, texTrackEtaYield_EventNorm, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); + Draw_TH1_Histograms(H1D_trackEta_rebinned, DatasetsNames, nDatasets, textContext, pdfName, texEtaMC, textYaxis, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); if (divideSuccess == true && options.find("ratio") != std::string::npos) { - Draw_TH1_Histograms(H1D_trackEta_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfNameEventNorm_ratio, texEtaMC, texRatioDatasets, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, "autoratio,noMarkerFirst"); + Draw_TH1_Histograms(H1D_trackEta_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_ratio, texEtaMC, texRatioDatasets, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, "autoratio,noMarkerFirst"); } else { cout << "Divide failed in Draw_Eta_DatasetComparison" << endl; @@ -2603,30 +2650,47 @@ void Draw_Phi_gen_DatasetComparison_H2CentVersion(std::string options) { H2D_centrality_track[iDataset] = (TH2D*)((TH2D*)file_O2Analysis_list[iDataset]->Get(analysisWorkflow[iDataset]+"/h2_centrality_particle_phi"))->Clone("Draw_Phi_DatasetComparison"+Datasets[iDataset]+DatasetsNames[iDataset]); - H1D_trackPhi[iDataset] = (TH1D*)H2D_centrality_track[iDataset]->ProjectionY("trackPhi_"+Datasets[iDataset]+DatasetsNames[iDataset], 1, H2D_centrality_track[iDataset]->GetNbinsX(), "e"); + H1D_trackPhi[iDataset] = (TH1D*)H2D_centrality_track[iDataset]->ProjectionY("trackPhi_"+Datasets[iDataset]+DatasetsNames[iDataset], 0, H2D_centrality_track[iDataset]->GetNbinsX(), "e"); H1D_trackPhi_rebinned[iDataset] = (TH1D*)H1D_trackPhi[iDataset]->Rebin(5.,"trackPhi_rebinned_"+Datasets[iDataset]+DatasetsNames[iDataset]); cout << "H1D_trackPhi_rebinned[iDataset]->GetEntries() preNorm = " << H1D_trackPhi_rebinned[iDataset]->Integral() << endl; - // NormaliseYieldToNEntries(H1D_trackPhi_rebinned[iDataset]); - if (isDatasetWeighted[iDataset]) { - Nevents = GetNEventsSelected_TrackEffWorkflow_gen_weighted(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); - } else { - Nevents = GetNEventsSelected_TrackEffWorkflow_gen(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + if (options.find("evtNorm") != std::string::npos) { + if (isDatasetWeighted[iDataset]) { + Nevents = GetNEventsSelected_TrackEffWorkflow_gen_weighted(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + } else { + Nevents = GetNEventsSelected_TrackEffWorkflow_gen(file_O2Analysis_list[iDataset], analysisWorkflow[iDataset]); + } + NormaliseYieldToNEvents(H1D_trackPhi_rebinned[iDataset], Nevents); } - NormaliseYieldToNEvents(H1D_trackPhi_rebinned[iDataset], Nevents); - cout << "H1D_trackPhi_rebinned[iDataset]->GetEntries() postNorm = " << H1D_trackPhi_rebinned[iDataset]->Integral() << ", Nevents = " << Nevents << endl; + if (options.find("entriesNorm") != std::string::npos) { + NormaliseYieldToIntegral(H1D_trackPhi_rebinned[iDataset]); + } + // cout << "H1D_trackPhi_rebinned[iDataset]->GetEntries() postNorm = " << H1D_trackPhi_rebinned[iDataset]->Integral() << ", Nevents = " << Nevents << endl; H1D_trackPhi_rebinned_ratios[iDataset] = (TH1D*)H1D_trackPhi_rebinned[iDataset]->Clone("trackPhi_rebinned_ratios"+Datasets[iDataset]+DatasetsNames[iDataset]); H1D_trackPhi_rebinned_ratios[iDataset]->Reset("M"); divideSuccess = H1D_trackPhi_rebinned_ratios[iDataset]->Divide(H1D_trackPhi_rebinned[iDataset], H1D_trackPhi_rebinned[0], 1., 1., datasetsAreSubsetsofId0 ? "b" : ""); } - TString* pdfName = new TString("track_Phi_gen_DataComp"); - TString* pdfName_ratio = new TString("track_Phi_gen_DataComp_ratio"); + TString* textYaxis; + TString pdfNameNorm; + if (options.find("evtNorm") != std::string::npos) { + textYaxis = texTrackPtYield_EventNorm; + pdfNameNorm = (TString)"_EventNorm"; + } + if (options.find("entriesNorm") != std::string::npos) { + textYaxis = texTrackPtYield_EntriesNorm; + pdfNameNorm = (TString)"_EntriesNorm"; + } + + TString* pdfName = new TString("track_Phi_gen_DataComp"+pdfNameNorm); + TString* pdfName_ratio = new TString("track_Phi_gen_DataComp"+pdfNameNorm+"_ratio"); + + TString textContext(contextTrackDatasetComp("")); - Draw_TH1_Histograms(H1D_trackPhi_rebinned, DatasetsNames, nDatasets, textContext, pdfName, texPhiMC, texTrackPhiYield_EventNorm, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); + Draw_TH1_Histograms(H1D_trackPhi_rebinned, DatasetsNames, nDatasets, textContext, pdfName, texPhiMC, textYaxis, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, ""); if (divideSuccess == true && options.find("ratio") != std::string::npos) { Draw_TH1_Histograms(H1D_trackPhi_rebinned_ratios, DatasetsNames, nDatasets, textContext, pdfName_ratio, texPhiMC, texRatioDatasets, texCollisionDataInfo, drawnWindowAuto, legendPlacementAuto, contextPlacementAuto, "autoratio,noMarkerFirst"); } diff --git a/Tracks/TrackMcQC_inputs_template.h b/Tracks/TrackMcQC_inputs_template.h index 6c17c4b..3989d4e 100644 --- a/Tracks/TrackMcQC_inputs_template.h +++ b/Tracks/TrackMcQC_inputs_template.h @@ -936,24 +936,83 @@ // const bool datasetsAreSubsetsofId0 = false; -//////// -------- test cleanUpCollHists //////// -TString* texCollisionDataInfo = new TString("PYTHIA MC #sqrt{#it{s}_{(NN)}} = 5.36 TeV"); -const TString* texDatasetsComparisonType = new TString("simType"); +// //////// -------- test cleanUpCollHists //////// +// TString* texCollisionDataInfo = new TString("PYTHIA pp MC #sqrt{#it{s}} = 5.36 TeV"); +// const TString* texDatasetsComparisonType = new TString(""); +// const TString* texDatasetsComparisonCommonDenominator = new TString(""); +// const int nDatasets = 2; +// const TString Datasets[nDatasets] = {"LHC26a6_train613380_TracksQA", "LHC25b4ab6_train610066_tracks"}; + +// // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"}; +// const TString DatasetsNames[nDatasets] = {"JJ LHC26a6", "MB LHC25b4ab6"}; +// // const TString DatasetsNames[nDatasets] = {"pp jet-jet LHC26a6"}; + +// TFile* file_O2Analysis_list[nDatasets] = {new TFile("../Datasets/"+Datasets[0]+"/AnalysisResults.root"), +// new TFile("../Datasets/"+Datasets[1]+"/AnalysisResults.root") +// }; + +// const TString analysisWorkflow[nDatasets] = {"track-efficiency", +// "track-efficiency" +// }; +// const TString wagonId[nDatasets] = {"", +// "" +// }; +// const bool isDatasetWeighted[nDatasets] = {true, false}; + +// const std::string histDatasetComparisonStructure = ""; +// const bool datasetsAreSubsetsofId0 = false; + + +//////// -------- LHC26a6_train621344_EffPurQA vs LHC25b4ab6_train610066_tracks //////// +TString* texCollisionDataInfo = new TString("PYTHIA pp MC #sqrt{#it{s}} = 5.36 TeV"); +const TString* texDatasetsComparisonType = new TString(""); const TString* texDatasetsComparisonCommonDenominator = new TString(""); const int nDatasets = 2; -const TString Datasets[nDatasets] = {"jetjet_ppAnchoredPbPb_5TeV_localNewCleanUpCollHists", "ppRefGenPurposeMC_5TeV_localNewCleanUpCollHists"}; +const TString Datasets[nDatasets] = {"LHC26a6_train621344_EffPurQA", "LHC25b4ab6_train610066_tracks"}; + // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"}; -const TString DatasetsNames[nDatasets] = {"pp jet-jet Pb-Pb anchor", "pp ref MB MC"}; -TFile* file_O2Analysis_list[nDatasets] = {new TFile("Datasets/"+Datasets[0]+"/AnalysisResults.root"), - new TFile("Datasets/"+Datasets[1]+"/AnalysisResults.root") +const TString DatasetsNames[nDatasets] = {"JJ LHC26a6", "MB LHC25b4ab6"}; +// const TString DatasetsNames[nDatasets] = {"pp jet-jet LHC26a6"}; + +TFile* file_O2Analysis_list[nDatasets] = {new TFile("../Datasets/"+Datasets[0]+"/AnalysisResults.root"), + new TFile("../Datasets/"+Datasets[1]+"/AnalysisResults.root") }; -const TString analysisWorkflow[nDatasets] = {"track-efficiency", +const TString analysisWorkflow[nDatasets] = {"track-efficiency_id47809", // for eff and purity : track-efficiency_id47809 / for Tracks QA : track-efficiency_id47557 "track-efficiency" }; const TString wagonId[nDatasets] = {"", "" }; const bool isDatasetWeighted[nDatasets] = {true, false}; + const std::string histDatasetComparisonStructure = ""; const bool datasetsAreSubsetsofId0 = false; + + + +// //////// -------- LHC26a6_train615296_Unf vs LHC25b4ab6_train610066_tracks //////// JUST FOR MC GEN PARTICLE TRACK SPECTRA QA BECAUSE THIS CONTAINS THE MC GEN COLLISIONS NUMBER +// TString* texCollisionDataInfo = new TString("PYTHIA pp MC #sqrt{#it{s}} = 5.36 TeV"); +// const TString* texDatasetsComparisonType = new TString(""); +// const TString* texDatasetsComparisonCommonDenominator = new TString(""); +// const int nDatasets = 2; +// const TString Datasets[nDatasets] = {"LHC26a6_train615296_Unf", "LHC25b4ab6_train610066_tracks"}; + +// // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"}; +// const TString DatasetsNames[nDatasets] = {"JJ LHC26a6", "MB LHC25b4ab6"}; +// // const TString DatasetsNames[nDatasets] = {"pp jet-jet LHC26a6"}; + +// TFile* file_O2Analysis_list[nDatasets] = {new TFile("../Datasets/"+Datasets[0]+"/AnalysisResults.root"), +// new TFile("../Datasets/"+Datasets[1]+"/AnalysisResults.root") +// }; + +// const TString analysisWorkflow[nDatasets] = {"track-efficiency", +// "track-efficiency" +// }; +// const TString wagonId[nDatasets] = {"", +// "" +// }; +// const bool isDatasetWeighted[nDatasets] = {true, false}; + +// const std::string histDatasetComparisonStructure = ""; +// const bool datasetsAreSubsetsofId0 = false; diff --git a/Tracks/TrackMcQC_settings_template.h b/Tracks/TrackMcQC_settings_template.h index 3b3387d..69e0d81 100644 --- a/Tracks/TrackMcQC_settings_template.h +++ b/Tracks/TrackMcQC_settings_template.h @@ -19,6 +19,6 @@ float arrayRadius[nRadius] = {0.2, 0.4, 0.6}; // float arrayRadius[nRadius] = {0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 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 iJetLevel = 2; const double deltaEtaMcVsTrackEfficiency = 0; \ No newline at end of file diff --git a/Tracks/TrackQC.C b/Tracks/TrackQC.C index 035b4f6..9db7801 100644 --- a/Tracks/TrackQC.C +++ b/Tracks/TrackQC.C @@ -20,8 +20,8 @@ // #include "RooUnfoldBinByBin.h" //My Libraries -#include "./TrackQC_settings.h" -#include "./TrackQC_inputs.h" +#include "./TrackQC_settings_template.h" +#include "./TrackQC_inputs_template.h" #include "../Settings/AxisTitles.h" #include "../Settings/GlobalSettings.h" #include "../Utilities/AnalysisUtilities.h" @@ -121,9 +121,9 @@ void TrackQC() { jetPtMaxCut = jetPtMinCutArray[iPtBin+1]; float ptRange[2] = {jetPtMinCut, jetPtMaxCut}; - Draw_Eta_DatasetComparison(ptRange, "evtNorm"); + // Draw_Eta_DatasetComparison(ptRange, "evtNorm"); // Draw_Eta_DatasetComparison(ptRange, "entriesNorm"); - Draw_Phi_DatasetComparison(ptRange, "evtNorm"); + // Draw_Phi_DatasetComparison(ptRange, "evtNorm"); // Draw_Phi_DatasetComparison(ptRange, "entriesNorm"); // Draw_Eta_DatasetComparison_trackSelComp(); diff --git a/Tracks/TrackQC_inputs_template.h b/Tracks/TrackQC_inputs_template.h index cfde9ea..07df3c1 100644 --- a/Tracks/TrackQC_inputs_template.h +++ b/Tracks/TrackQC_inputs_template.h @@ -1071,31 +1071,100 @@ TFile* file_AliAnalysis; //dummy // const bool datasetsAreSubsetsofId0 = false; // const bool trackHistsObsoleteVersion[nDatasets] = {true, true, true, true}; -//////// -------- Efficiency differences between pp only sims //////// -TString* texCollisionDataInfo = new TString("PYTHIA MC"); -const TString* texDatasetsComparisonType = new TString("simType"); -const TString* texDatasetsComparisonCommonDenominator = new TString("pp simulation"); -const int nDatasets = 4; -const TString Datasets[nDatasets] = {"jetjet_PbPbAnchorMC_5TeV_train428451", "ppRefGenPurposeMC_5TeV_train429398", "jetjet_ppAnchorMC_13TeV_train429402", "ppGenPurposeMC_13TeV_train428885"}; +// //////// -------- Efficiency differences between pp only sims //////// +// TString* texCollisionDataInfo = new TString("PYTHIA MC"); +// const TString* texDatasetsComparisonType = new TString("simType"); +// const TString* texDatasetsComparisonCommonDenominator = new TString("pp simulation"); +// const int nDatasets = 4; +// const TString Datasets[nDatasets] = {"jetjet_PbPbAnchorMC_5TeV_train428451", "ppRefGenPurposeMC_5TeV_train429398", "jetjet_ppAnchorMC_13TeV_train429402", "ppGenPurposeMC_13TeV_train428885"}; +// // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"}; +// const TString DatasetsNames[nDatasets] = {"jet-jet 5.36TeV fullstats", "ref gen.purp. 5.36TeV", "jet-jet 13.6TeV", "gen.purp. 13.6TeV"}; +// TFile* file_O2Analysis_list[nDatasets] = {new TFile("Datasets/"+Datasets[0]+"/AnalysisResults.root"), +// new TFile("Datasets/"+Datasets[1]+"/AnalysisResults.root"), +// new TFile("Datasets/"+Datasets[2]+"/AnalysisResults.root"), +// new TFile("Datasets/"+Datasets[3]+"/AnalysisResults.root") +// }; + +// const TString analysisWorkflow[nDatasets] = {"track-efficiency_id27663", +// "track-efficiency_id31073", +// "track-efficiency_id30281", +// "track-efficiency_id30823" +// }; +// const TString wagonId[nDatasets] = {"", +// "", +// "", +// "" +// }; +// const bool isDatasetWeighted[nDatasets] = {true, false, true, false}; +// const std::string histDatasetComparisonStructure = "twoByTwoDatasetPairs"; // +// const bool datasetsAreSubsetsofId0 = false; +// const bool trackHistsObsoleteVersion[nDatasets] = {true, true, true, true}; + + + +// //////// -------- Track QC LHC26a6 //////// +// TString* texCollisionDataInfo = new TString("PYTHIA MC"); +// const TString* texDatasetsComparisonType = new TString(""); +// const TString* texDatasetsComparisonCommonDenominator = new TString("pp simulation"); +// const int nDatasets = 2; +// const TString Datasets[nDatasets] = {"LHC26a6_train613380_TracksQA", "LHC25b4ab6_train610066_tracks"}; +// // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"}; +// const TString DatasetsNames[nDatasets] = {"JJ LHC26a6", "MB LHC25b4ab6"}; +// TFile* file_O2Analysis_list[nDatasets] = {new TFile("../Datasets/"+Datasets[0]+"/AnalysisResults.root"), +// new TFile("../Datasets/"+Datasets[1]+"/AnalysisResults.root") +// }; + +// const TString analysisWorkflow[nDatasets] = {"track-efficiency", "track-efficiency", +// }; + +// const TString wagonId[nDatasets] = {"", "" +// }; + +// const bool isDatasetWeighted[nDatasets] = {true, false}; +// const std::string histDatasetComparisonStructure = ""; // +// const bool datasetsAreSubsetsofId0 = false; +// const bool trackHistsObsoleteVersion[nDatasets] = {true, true}; + + + + + +//////// -------- LHC26a6_train621344_EffPurQA vs LHC25b4ab6_train610066_tracks //////// +TString* texCollisionDataInfo = new TString("PYTHIA pp MC #sqrt{#it{s}} = 5.36 TeV"); +const TString* texDatasetsComparisonType = new TString(""); +const TString* texDatasetsComparisonCommonDenominator = new TString(""); +const int nDatasets = 2; +const TString Datasets[nDatasets] = {"LHC26a6_train621344_EffPurQA", "LHC25b4ab6_train610066_tracks"}; + // const TString DatasetsNames[nDatasets] = {"0-10%", "50-90%"}; -const TString DatasetsNames[nDatasets] = {"jet-jet 5.36TeV fullstats", "ref gen.purp. 5.36TeV", "jet-jet 13.6TeV", "gen.purp. 13.6TeV"}; -TFile* file_O2Analysis_list[nDatasets] = {new TFile("Datasets/"+Datasets[0]+"/AnalysisResults.root"), - new TFile("Datasets/"+Datasets[1]+"/AnalysisResults.root"), - new TFile("Datasets/"+Datasets[2]+"/AnalysisResults.root"), - new TFile("Datasets/"+Datasets[3]+"/AnalysisResults.root") +const TString DatasetsNames[nDatasets] = {"JJ LHC26a6", "MB LHC25b4ab6"}; +// const TString DatasetsNames[nDatasets] = {"pp jet-jet LHC26a6"}; + +TFile* file_O2Analysis_list[nDatasets] = {new TFile("../Datasets/"+Datasets[0]+"/AnalysisResults.root"), + new TFile("../Datasets/"+Datasets[1]+"/AnalysisResults.root") }; -const TString analysisWorkflow[nDatasets] = {"track-efficiency_id27663", - "track-efficiency_id31073", - "track-efficiency_id30281", - "track-efficiency_id30823" +const TString analysisWorkflow[nDatasets] = {"track-efficiency_id47557", // for eff and purity : track-efficiency_id47809 / for Tracks QA : track-efficiency_id47557 + "track-efficiency" }; const TString wagonId[nDatasets] = {"", - "", - "", "" }; -const bool isDatasetWeighted[nDatasets] = {true, false, true, false}; -const std::string histDatasetComparisonStructure = "twoByTwoDatasetPairs"; // +const bool isDatasetWeighted[nDatasets] = {true, false}; + +const std::string histDatasetComparisonStructure = ""; const bool datasetsAreSubsetsofId0 = false; -const bool trackHistsObsoleteVersion[nDatasets] = {true, true, true, true}; +const bool trackHistsObsoleteVersion[nDatasets] = {true, true}; + + + + + + + + + + + + + diff --git a/Tracks/TrackQC_settings_template.h b/Tracks/TrackQC_settings_template.h index 9faa968..5bd1254 100644 --- a/Tracks/TrackQC_settings_template.h +++ b/Tracks/TrackQC_settings_template.h @@ -4,12 +4,14 @@ // float GLOBAL_epsilon = 0.00001; // Analysis settings -const int nCollSystems = 2; -const TString collSystems[nCollSystems] = {"Pb-Pb", "pp"}; +const int nCollSystems = 1; +const TString collSystems[nCollSystems] = {"pp"}; -const int iCollSystem = 1; +const int iCollSystem = 1; // if 1 starting from 0 instead of 1 // const int nCentralityBins = 7; // const float arrayCentralityBinning[nCentralityBins+1] = {0, 10, 20, 30, 50, 70, 100}; -const int nCentralityBins = 3; -const float arrayCentralityBinning[nCentralityBins+1] = {0, 10, 50, 70}; \ No newline at end of file +const int nCentralityBins = 1; +const float arrayCentralityBinning[nCentralityBins+1] = {-10, 0}; + + diff --git a/Utilities/HistogramUtilities.C b/Utilities/HistogramUtilities.C index 784f257..d91b0b8 100644 --- a/Utilities/HistogramUtilities.C +++ b/Utilities/HistogramUtilities.C @@ -558,6 +558,48 @@ TH1D GetMatrixVectorProductTH2xTH1(TH2D* histA, TH1D* histU){ } +void Compute_TH1D_RelativeUncertainty( + TH1D* hInput[], // input histograms + TH1D* hRelUnc[], // output histograms (must be allocated outside) + int nDatasets ){ + // Example usage: + // TH1D* H1D_trackPt_efficiency_concatenated[nDatasets]; + // TH1D* H1D_trackPt_efficiency_relUnc[nDatasets]; + // Compute_TH1D_RelativeUncertainty(H1D_trackPt_efficiency_concatenated,H1D_trackPt_efficiency_relUnc,nDatasets); + + for (int iDataset = 0; iDataset < nDatasets; ++iDataset){ + if (!hInput[iDataset]) continue; + // Clone structure (axes, binning) but reset content + hRelUnc[iDataset] = (TH1D*)hInput[iDataset]->Clone(Form("%s_relUnc", hInput[iDataset]->GetName())); + hRelUnc[iDataset]->Reset(); + hRelUnc[iDataset]->SetTitle(Form("Relative uncertainty of %s", hInput[iDataset]->GetTitle())); + + int nBins = hInput[iDataset]->GetNbinsX(); + + for (int bin = 1; bin <= nBins; ++bin){ + double content = hInput[iDataset]->GetBinContent(bin); + double error = hInput[iDataset]->GetBinError(bin); + + double relUnc = 0.0; + + if (content != 0.0) + relUnc = error / content; + + hRelUnc[iDataset]->SetBinContent(bin, relUnc*100.0); // store as percentage + hRelUnc[iDataset]->SetBinError(bin, 0.0); // usually no error on relative uncertainty + + // cout << "Dataset: " << iDataset + // // << " | Bin: " << bin + // // << " | BinCenter: " << hInput[iDataset]->GetBinCenter(bin) + // // << " | Content: " << content + // // << " | Error: " << error + // << " | RelUnc: " << relUnc*100.0 + // << endl; + } + } +} + + diff --git a/Utilities/HistogramUtilities.h b/Utilities/HistogramUtilities.h index c30e1f3..51f985e 100644 --- a/Utilities/HistogramUtilities.h +++ b/Utilities/HistogramUtilities.h @@ -26,4 +26,6 @@ TH2D GetTransposeHistogram(TH2D* inputHist); TH2D GetMatrixProductTH2xTH2(TH2D* histA, TH2D* histB); TH1D GetMatrixVectorProductTH2xTH1(TH2D* histA, TH1D* histU); +void Compute_TH1D_RelativeUncertainty(TH1D* hInput[], TH1D* hRelUnc[], int nDatasets); + #endif