|
| 1 | +############################################################################# |
| 2 | +## © Copyright CERN 2018. All rights not expressly granted are reserved. ## |
| 3 | + |
| 4 | +## This program is free software: you can redistribute it and/or modify it ## |
| 5 | +## under the terms of the GNU General Public License as published by the ## |
| 6 | +## Free Software Foundation, either version 3 of the License, or (at your ## |
| 7 | +## option) any later version. This program is distributed in the hope that ## |
| 8 | +## it will be useful, but WITHOUT ANY WARRANTY; without even the implied ## |
| 9 | +## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. ## |
| 10 | +## See the GNU General Public License for more details. ## |
| 11 | +## You should have received a copy of the GNU General Public License ## |
| 12 | +## along with this program. if not, see <https://www.gnu.org/licenses/>. ## |
| 13 | +############################################################################# |
| 14 | + |
| 15 | +""" |
| 16 | +main script for doing final stage analysis |
| 17 | +
|
| 18 | +NB: Duplicate of macro in AN Note repository. Will not work here! |
| 19 | +Just as an example how functions can be used |
| 20 | +""" |
| 21 | +import os |
| 22 | +# pylint: disable=import-error, no-name-in-module, unused-import |
| 23 | +import yaml |
| 24 | +from ROOT import TFile, gStyle, gROOT, TH1F, TGraphAsymmErrors, TH1 |
| 25 | +from ROOT import kBlue, kAzure, kOrange, kGreen, kBlack, kRed, kWhite |
| 26 | +from ROOT import Double, TCanvas, gPad |
| 27 | +from machine_learning_hep.utilities_plot import Errors |
| 28 | +from machine_learning_hep.utilities_plot import average_pkpi_pk0s |
| 29 | + |
| 30 | +def extract_histo(histo_name, path): |
| 31 | + |
| 32 | + in_file = TFile.Open(path, "READ") |
| 33 | + histo = in_file.Get(histo_name) |
| 34 | + if isinstance(histo, TH1): |
| 35 | + histo.SetDirectory(0) |
| 36 | + return histo |
| 37 | + |
| 38 | +def make_standard_save_path(prefix, filepath): |
| 39 | + |
| 40 | + folder_plots = f"{filepath}" |
| 41 | + if not os.path.exists(folder_plots): |
| 42 | + print("creating folder ", folder_plots) |
| 43 | + os.makedirs(folder_plots) |
| 44 | + return f"{folder_plots}/{prefix}.eps" |
| 45 | + |
| 46 | + |
| 47 | +############################################################################# |
| 48 | +############################# Input arguments ############################## |
| 49 | +############################################################################# |
| 50 | + |
| 51 | +SIGMAV0 = 57.8e9 |
| 52 | +BRPKPI = 0.0623 |
| 53 | +BRPK0S = 0.0109 |
| 54 | +PTBINS = 6 |
| 55 | + |
| 56 | +FILE_PKPI = "inputfiles_HP20/finalcrossLcpKpippMBvspt_ntrklmult4.root" |
| 57 | +FILE_PK0S = "inputfiles_HP20/finalcrossLcpK0sppMBvspt_ntrklmult4.root" |
| 58 | +FILE_PK0S2 = "inputfiles_HP20/FeeddownSyst_NbNbx2_LcpK0sppMBvspt_ntrkl.root" |
| 59 | +FILE_OUT = "inputfiles_HP20/finalcrossLcAverageMBvspt_ntrklmult4.root" |
| 60 | + |
| 61 | +ERROR_PKPI = "syst_HP20/syst_forLcaverage/LcpKpipp/errors_histoSigmaCorr_4.yaml" |
| 62 | +ERROR_PK0S = "syst_HP20/syst_forLcaverage/LcpK0spp/errors_histoSigmaCorr_4.yaml" |
| 63 | +ERROR_OUT = "syst_HP20/syst_forLcaverage/errors_avg_histoSigmaCorr_4.yaml" |
| 64 | + |
| 65 | +gROOT.SetBatch(True) |
| 66 | + |
| 67 | + |
| 68 | +############################################################################# |
| 69 | +################################# Average ################################### |
| 70 | +############################################################################# |
| 71 | + |
| 72 | +HISTO_PKPI = extract_histo("histoSigmaCorr_rebin", FILE_PKPI) |
| 73 | +HISTO_PKPI.SetName(f"histoSigmaCorr_pKpi") |
| 74 | +HISTO_PKPI.Scale(1./BRPKPI) |
| 75 | +HISTO_PK0S = extract_histo("histoSigmaCorr", FILE_PK0S) |
| 76 | +HISTO_PK0S.SetName(f"histoSigmaCorr_pK0s") |
| 77 | +HISTO_PK0S.Scale(1./BRPK0S) |
| 78 | +HISTO_PK0S.Scale(1./SIGMAV0) |
| 79 | + |
| 80 | +GRFD_PKPI = extract_histo("gFcCorrConservative", FILE_PKPI) |
| 81 | +GRFD_PK0S = extract_histo("gNbCorrConservative", FILE_PK0S2) |
| 82 | + |
| 83 | +DICTEXTRA_PKPI = {} |
| 84 | +HISTOEFF_PKPI = extract_histo("hDirectEffpt", FILE_PKPI) |
| 85 | +ERROREFF_PKPI = [] |
| 86 | +for i in range(HISTOEFF_PKPI.GetNbinsX()): |
| 87 | + print(HISTOEFF_PKPI.GetBinCenter(i+1)) |
| 88 | + RELERREFF = HISTOEFF_PKPI.GetBinError(i+1) / HISTOEFF_PKPI.GetBinContent(i+1) |
| 89 | + if i == 0: |
| 90 | + ERROREFF_PKPI.append([0, 0, 99, 99]) |
| 91 | + else: |
| 92 | + ERROREFF_PKPI.append([0, 0, RELERREFF, RELERREFF]) |
| 93 | +DICTEXTRA_PKPI["statunceff"] = ERROREFF_PKPI |
| 94 | +ERRSPKPI = Errors(PTBINS) |
| 95 | +ERRSPKPI.read(ERROR_PKPI, DICTEXTRA_PKPI) |
| 96 | + |
| 97 | +DICTEXTRA_PK0S = {} |
| 98 | +HISTOEFF_PK0S = extract_histo("hDirectEffpt", FILE_PK0S) |
| 99 | +ERROREFF_PK0S = [] |
| 100 | +for i in range(HISTOEFF_PK0S.GetNbinsX()): |
| 101 | + RELERREFF = HISTOEFF_PK0S.GetBinError(i+1) / HISTOEFF_PK0S.GetBinContent(i+1) |
| 102 | + ERROREFF_PK0S.append([0, 0, RELERREFF, RELERREFF]) |
| 103 | +DICTEXTRA_PK0S["statunceff"] = ERROREFF_PK0S |
| 104 | +ERRSPK0S = Errors(PTBINS) |
| 105 | +ERRSPK0S.read(ERROR_PK0S, DICTEXTRA_PK0S) |
| 106 | + |
| 107 | +MATCHPKPI = [-99, 1, 2, 3, 4, 5] |
| 108 | +MATCHPK0S = [1, 2, 3, 4, 5, 6] |
| 109 | +MATCHPKPIGR = [-99, 2, 3, 4, 5, 6] #Empty bin 0-1 still in fprompt tgraph |
| 110 | +MATCHPK0SGR = [1, 2, 3, 4, 5, 6] |
| 111 | + |
| 112 | +AVGCORRYIELD, AVGSTATUNC, AVGFPROMPT, \ |
| 113 | + AVGFPROMPTLOW, AVGFPROMPTHIGH, AVGERROR = average_pkpi_pk0s(HISTO_PKPI, HISTO_PK0S, |
| 114 | + GRFD_PKPI, GRFD_PK0S, |
| 115 | + ERRSPKPI, ERRSPK0S, |
| 116 | + MATCHPKPI, MATCHPK0S, |
| 117 | + MATCHPKPIGR, MATCHPK0SGR) |
| 118 | + |
| 119 | +HISTAVG = HISTO_PK0S.Clone("histoSigmaCorr_average") |
| 120 | +GRFDAVG = TGraphAsymmErrors(PTBINS) |
| 121 | +for ipt in range(PTBINS): |
| 122 | + HISTAVG.SetBinContent(ipt+1, AVGCORRYIELD[ipt]) |
| 123 | + HISTAVG.SetBinError(ipt+1, AVGSTATUNC[ipt]) |
| 124 | + GRFDAVG.SetPoint(ipt+1, GRFD_PK0S.GetX()[ipt+1], AVGFPROMPT[ipt]) |
| 125 | + GRFDAVG.SetPointError(ipt+1, GRFD_PK0S.GetEXlow()[ipt+1], GRFD_PK0S.GetEXhigh()[ipt+1], \ |
| 126 | + AVGFPROMPTLOW[ipt], AVGFPROMPTHIGH[ipt]) |
| 127 | + |
| 128 | +AVGERROR.print() |
| 129 | +print("\n\n Store above in", ERROR_OUT) |
| 130 | + |
| 131 | +print("\n\n\nStoring ROOT objects in", FILE_OUT) |
| 132 | +FOUT = TFile(FILE_OUT, "RECREATE") |
| 133 | +FOUT.cd() |
| 134 | +HISTAVG.Write("histoSigmaCorr_average") |
| 135 | +HISTO_PKPI.Write("histoSigmaCorr_pKpi") |
| 136 | +HISTO_PK0S.Write("histoSigmaCorr_pK0s") |
| 137 | + |
| 138 | +GRFDAVG.Write("gFcCorrConservative_average") |
| 139 | +GRFD_PKPI.Write("gFcCorrConservative_pKpi") |
| 140 | +GRFD_PK0S.Write("gNbCorrConservative_pK0s") |
| 141 | +FOUT.Close() |
0 commit comments