diff --git a/ComputeSystematicsCF.py b/ComputeSystematicsCF.py new file mode 100644 index 0000000..773489f --- /dev/null +++ b/ComputeSystematicsCF.py @@ -0,0 +1,156 @@ +''' +Script to compute the raw correlation function. +The output file is RawSystCF_suffix.root + +Usage: +python3 ComputeSystematicsCF.py cfg.yml + +''' +import os +import argparse +import yaml +from itertools import chain + +from ROOT import TFile, TH2F, TH1D, TH2D + +from fempy import logger as log +from fempy.utils.io import Load, GetKeyNames + +parser = argparse.ArgumentParser() +parser.add_argument('cfg', default='') +parser.add_argument('--debug', default=False, action='store_true') +args = parser.parse_args() + +if args.debug: + log.setLevel(1) + +# Load yaml file +with open(args.cfg, "r") as stream: + try: + cfg = yaml.safe_load(stream) + except yaml.YAMLError as exc: + log.critical('Yaml configuration could not be loaded. Is it properly formatted?') + +# Define the output file +oFileBaseName = 'RawSystCF' +if cfg['suffix'] != '' and cfg['suffix'] is not None: + oFileBaseName += f'_{cfg["suffix"]}' +oFileName = os.path.join(cfg['odir'], oFileBaseName + '.root') + +# Open the output file +try: + oFile = TFile.Open(oFileName, 'recreate') +except OSError: + log.critical('The output file %s is not writable', oFileName) +combs = ['p02', 'p03', 'p12', 'p13'] +for comb in combs: + oFile.mkdir(comb) + +totalSystVars = len(list(chain.from_iterable(cfg['suffixessyst']))) +hFemtoPairs = TH2D("hFemtoPairs", "Pairs in femto region", totalSystVars, 0, totalSystVars, len(combs), 0, len(combs)) +hFemtoPairs.SetStats(0) +hFemtoPairs.GetYaxis().SetLabelSize(50) +hFemtoPairs.GetXaxis().SetLabelSize(30) +hFemtoPairs.GetXaxis().CenterLabels() + +# Compute the correlation functions for all systematic variations +hSE = {} +hME = {} +hMErew = {} +for iSystFile, systFileName in enumerate(cfg['infilessyst']): + systFile = TFile(systFileName) + for iSystVar in cfg['suffixessyst'][iSystFile]: + for iComb, comb in enumerate(combs): + iPart1 = int(comb[1]) + iPart2 = int(comb[2]) + oFile.cd(comb) + + if f'HMResults{iSystVar}' in GetKeyNames(systFile): # Make correlation functions from FemtoDream + fdcomb = f'Particle{iPart1}_Particle{iPart2}' + # The histograms are casted to TH1D with TH1::Copy to avoid NotImplementedError when computing hSE/hME + hSE[comb] = TH1D() + Load(systFile, f'HMResults{iSystVar}/HMResults{iSystVar}/{fdcomb}/SEDist_{fdcomb}').Copy(hSE[comb]) + + + hFemtoPairs.Fill(iSystVar-0.5, len(combs)-iComb-0.5, + hSE[comb].Integral(hSE[comb].FindBin(0.0001), hSE[comb].FindBin(0.2*0.9999))) + hFemtoPairs.GetYaxis().SetBinLabel(len(combs)-iComb, comb) + hFemtoPairs.GetXaxis().SetBinLabel(iSystVar, str(iSystVar)) + + # The histograms are casted to TH1D with TH1::Copy to avoid NotImplementedError when computing hSE/hME + hME[comb] = TH1D() + Load(systFile, f'HMResults{iSystVar}/HMResults{iSystVar}/{fdcomb}/MEDist_{fdcomb}').Copy(hME[comb]) + + # No need to cast the these because the projection of TH2D and TH2F is always TH1D + hSEmultk = Load(systFile, f'HMResults{iSystVar}/HMResults{iSystVar}/{fdcomb}/SEMultDist_{fdcomb}') + hMEmultk = Load(systFile, f'HMResults{iSystVar}/HMResults{iSystVar}/{fdcomb}/MEMultDist_{fdcomb}') + + nbins = hMEmultk.ProjectionX().GetNbinsX() + hMEreweightk = TH1D("MErewdistr", "MErewdistr", nbins, hMEmultk.GetXaxis().GetXmin(), hMEmultk.GetXaxis().GetXmax()) + for iBin in range(hMEmultk.ProjectionY().GetNbinsX() + 2): # Loop over underflow, all bins, and overflow + hSEbinmult = hSEmultk.ProjectionX(f'{comb}SEdistr', iBin, iBin) + hMEbinmult = hMEmultk.ProjectionX(f'{comb}MEdistr', iBin, iBin) + if(hMEbinmult.Integral() > 0): + hMEreweightk.Add(hMEbinmult, hSEbinmult.Integral()/hMEbinmult.Integral()) + hMErew[comb] = hMEreweightk + + # Sum pair and antipair + for comb in combs: + hSE['p02_13'] = hSE['p02'] + hSE['p13'] + hME['p02_13'] = hME['p02'] + hME['p13'] + hMErew['p02_13'] = hMErew['p02'] + hMErew['p13'] + + hSE['p03_12'] = hSE['p03'] + hSE['p12'] + hME['p03_12'] = hME['p03'] + hME['p12'] + hMErew['p03_12'] = hMErew['p03'] + hMErew['p12'] + + # Compute the CF and write to file + for comb in combs + ['p02_13', 'p03_12']: + rebin = round(float(cfg['binwidth']) / (hSE[comb].GetBinWidth(1) * 1000)) + hSE[comb].Rebin(rebin) + hME[comb].Rebin(rebin) + hMErew[comb].Rebin(rebin) + + if cfg['norm'] is None: # if not specified, normalize to the yields + norm = hME[comb].GetEntries() / hSE[comb].GetEntries() + normrew = hMErew[comb].GetEntries() / hSE[comb].GetEntries() + else: + firstBin = hSE[comb].FindBin(cfg['norm'][0]*1.0001) + lastBin = hSE[comb].FindBin(cfg['norm'][1]*0.9999) + norm = hME[comb].Integral(firstBin, lastBin) / hSE[comb].Integral(firstBin, lastBin) + normrew = hMErew[comb].Integral(firstBin, lastBin) / hSE[comb].Integral(firstBin, lastBin) + + hSE[comb].Sumw2() + hME[comb].Sumw2() # Just to trigger the same draw option as for hSE + + hCF = norm * hSE[comb] / hME[comb] + hCFrew = normrew * hSE[comb] / hMErew[comb] + + oFile.mkdir(f'{comb}/var{iSystVar}') + oFile.cd(f'{comb}/var{iSystVar}') + + hSE[comb].SetName(f'hSE_{iSystVar}') + hSE[comb].SetTitle(';#it{k}* (GeV/#it{c});Counts') + hSE[comb].Write() + + hME[comb].SetName(f'hME_{iSystVar}') + hME[comb].SetTitle(';#it{k}* (GeV/#it{c});Counts') + hME[comb].Write() + + hMErew[comb].SetName(f'hMErew_{iSystVar}') + hMErew[comb].SetTitle(';#it{k}* (GeV/#it{c});Counts') + + hCF.SetName(f'hCF_{iSystVar}') + hCF.SetTitle(';#it{k}* (GeV/#it{c});#it{C}(#it{k}*)') + hCF.Write() + + hCFrew.SetName(f'hCFrew_{iSystVar}') + hCFrew.SetTitle(';#it{k}* (GeV/#it{c});#it{C}(#it{k}*)') + hCFrew.Write() + + print(f'Finished systematic variation {iSystVar}') + +oFile.cd() +hFemtoPairs.Write() +oFile.Close() +print(f'output saved in {oFileName}') diff --git a/SystematicsCFAnalysis.py b/SystematicsCFAnalysis.py new file mode 100644 index 0000000..77b9a3c --- /dev/null +++ b/SystematicsCFAnalysis.py @@ -0,0 +1,201 @@ +''' +Script to compute the systematic uncertainties on the CF +bin values. +The output file is SystUncCF_suffix.root + +Usage: +python3 SystematicsCFAnalysis.py cfg.yml + +''' +import os +import argparse +import yaml +import math +import numpy as np +from itertools import chain + +from ROOT import TFile, TF1, TH1D, TH2D + +from fempy import logger as log +from fempy.utils.io import Load, GetKeyNames + +parser = argparse.ArgumentParser() +parser.add_argument('cfg', default='') +parser.add_argument('--debug', default=False, action='store_true') +args = parser.parse_args() + +if args.debug: + log.setLevel(1) + +# Load yaml file +with open(args.cfg, "r") as stream: + try: + cfg = yaml.safe_load(stream) + except yaml.YAMLError as exc: + log.critical('Yaml configuration could not be loaded. Is it properly formatted?') + +# check if the CFs for systematic variations have already been calculated, +# else compute them +oFileCFBaseName = 'RawSystCF' + f'_{cfg["suffix"]}' +oFileCFName = os.path.join(cfg['odir'], oFileCFBaseName + '.root') +if os.path.exists(os.path.expanduser(oFileCFName)): + print('Systematics CFs already computed!') + inFileCFSystVar = TFile(oFileCFName) +else: + print('Systematics CFs not yet computed!') + os.system(f"python3 ComputeSystematicsCF.py {args.cfg}") + print('Systematics CFs computed!') + inFileCFSystVar = TFile(oFileCFName) +combs = [key.GetName() for key in inFileCFSystVar.GetListOfKeys() if key.GetName()[0] == 'p'] + +# Define the output file +oFileBaseName = 'SystUncCF' +if cfg['suffix'] != '' and cfg['suffix'] is not None: + oFileBaseName += f'_{cfg["suffix"]}' +oFileName = os.path.join(cfg['odir'], oFileBaseName + '.root') + +# Open the output file +try: + oFile = TFile.Open(oFileName, 'recreate') +except OSError: + log.critical('The output file %s is not writable', oFileName) + +# Load the CF selected for the analysis +inFileCFSelected = TFile(cfg['infileCFselected']) + +# TH2D to store deviations from selected CF +systVars = list(chain.from_iterable(cfg['suffixessyst'])) +nSystVars = max(systVars) - min(systVars) + 1 +nBinsKStar = len(range(cfg['systevalmaxbin'])) +uppEdgeKStar = cfg['systevalmaxbin']*cfg['binwidth'] + +# Compute the correlation functions for all systematic variations +for comb in combs: + print(f'Picking pair {comb}') + oFile.mkdir(f'{comb}') + oFile.cd(f'{comb}') + + hSystUnc = TH1D("hSystUnc", "hSystUnc", nBinsKStar, 0, uppEdgeKStar) + hRelSystUnc = TH1D("hRelSystUnc", "hRelSystUnc", nBinsKStar, 0, uppEdgeKStar) + hRelStatUnc = TH1D("hRelStatUnc", "hRelStatUnc", nBinsKStar, 0, uppEdgeKStar) + hRatioStatSystUnc = TH1D("hRatioStatSystUnc", "hRatioStatSystUnc", nBinsKStar, 0, uppEdgeKStar) + hCFWithSystUnc = TH1D("hCFWithSystUnc", "hCFWithSystUnc", nBinsKStar, 0, uppEdgeKStar) + hCFWithStatSystUnc = TH1D("hCFWithStatSystUnc", "hCFWithStatSystUnc", nBinsKStar, 0, uppEdgeKStar) + hResiduals = TH2D("hResiduals", "Syst vars residuals", nBinsKStar, 0, uppEdgeKStar, + nSystVars, min(systVars), max(systVars)) + hResiduals.SetStats(0) + + hCFSelected = Load(inFileCFSelected, f'{comb}/sgn/hCFrew') + if cfg.get('gaussianfits'): + # list of histograms each containing the CF entries for + # all variations for a specific bin + hCFBinEntries = [] + for iBin in range(cfg['systevalmaxbin']): + hCFBinEntries.append(TH1D(f"h{round(cfg['binwidth']*iBin+(cfg['binwidth']/2))}MeV", + f"h{round(cfg['binwidth']*iBin+(cfg['binwidth']/2))}MeV", + cfg['systhistobins'], hCFSelected.GetBinContent(iBin+1)-cfg['systhistointerval'], + hCFSelected.GetBinContent(iBin+1)+cfg['systhistointerval'])) + + for iSystVar in systVars: + print(f'Picking syst variation number {iSystVar}') + hSystVar = Load(inFileCFSystVar, f'{comb}/var{iSystVar}/hCFrew_{iSystVar}') + print(f'Picked syst variation number {iSystVar}') + hResiduals.GetYaxis().SetBinLabel(iSystVar-min(systVars)+1, str(iSystVar)) + + for iBin in range(nBinsKStar): + hCFBinEntries[iBin].Fill(hSystVar.GetBinContent(iBin+1)) + hResiduals.Fill(cfg['binwidth']*(iBin+1)+cfg['binwidth']/2, iSystVar, + (hCFSelected.GetBinContent(iBin+1)-hSystVar.GetBinContent(iBin+1)) / + hCFSelected.GetBinContent(iBin+1)) + + hCFYields = TH1D("hCFYields", "hCFYields", nBinsKStar, 0, uppEdgeKStar) + hCFMeans = TH1D("hCFmeans", "hCFmeans", nBinsKStar, 0, uppEdgeKStar) + + print('Evaluating the single bins ...') + oFile.mkdir(f'{comb}/binsfits') + for iBin, ihCFbin in enumerate(hCFBinEntries): + oFile.cd(f'{comb}/binsfits') + gaus = TF1(f"gaus_{iBin}", "gaus", ihCFbin.GetBinLowEdge(1), + ihCFbin.GetBinLowEdge(ihCFbin.GetNbinsX()) + ihCFbin.GetBinWidth(1)) + gaus.SetParameter(1, hCFSelected.GetBinContent(iBin+1)) + ihCFbin.Fit(gaus, "SMRL+", "") + ihCFbin.Write(f"hCFbin{round(cfg['binwidth']*iBin+(cfg['binwidth']/2))}MeV") + hCFYields.SetBinContent(iBin+1, gaus.GetParameter(0)) + hCFMeans.SetBinContent(iBin+1, gaus.GetParameter(1)) + hSystUnc.SetBinContent(iBin+1, gaus.GetParameter(2)) + hRelSystUnc.SetBinContent(iBin+1, hSystUnc.GetBinContent(iBin+1)/hCFSelected.GetBinContent(iBin+1)) + hRelStatUnc.SetBinContent(iBin+1, hCFSelected.GetBinError(iBin+1)/hCFSelected.GetBinContent(iBin+1)) + hRatioStatSystUnc.SetBinContent(iBin+1, hCFSelected.GetBinError(iBin+1)/hSystUnc.GetBinContent(iBin+1)) + hCFWithSystUnc.SetBinContent(iBin+1, hCFSelected.GetBinContent(iBin+1)) + hCFWithSystUnc.SetBinError(iBin+1, hSystUnc.GetBinContent(iBin+1)**2) + hCFWithStatSystUnc.SetBinContent(iBin+1, hCFSelected.GetBinContent(iBin+1)) + hCFWithStatSystUnc.SetBinError(iBin+1, math.sqrt(hCFSelected.GetBinError(iBin+1)**2 + + hSystUnc.GetBinContent(iBin+1)**2 ) ) + + hCFYields.Write() + hCFMeans.Write() + + else: + binsStdDevs = [] + variedCFsEntries = [] + + for iSystVar in systVars: + iSystVarBinEntries = [] + print(f'Picking syst variation number {iSystVar}') + hSystVar = Load(inFileCFSystVar, f'{comb}/var{iSystVar}/hCFrew_{iSystVar}') + print(f'Picked syst variation number {iSystVar}') + hResiduals.GetYaxis().SetBinLabel(iSystVar-min(systVars)+1, str(iSystVar)) + + for iBin in range(nBinsKStar): + iSystVarBinEntries.append(hSystVar.GetBinContent(iBin+1)) + hResiduals.Fill(cfg['binwidth']*(iBin+1)+cfg['binwidth']/2, iSystVar, + (hCFSelected.GetBinContent(iBin+1)-hSystVar.GetBinContent(iBin+1)) / + hCFSelected.GetBinContent(iBin+1)) + + variedCFsEntries.append(iSystVarBinEntries) + + binsVarEntries = [[variedCFsEntries[j][i] for j in range(len(variedCFsEntries))] for i in range(len(variedCFsEntries[0]))] + for iBin, binVar in enumerate(binsVarEntries): + binsStdDevs.append(np.std(binVar)) + hSystUnc.SetBinContent(iBin+1, binsStdDevs[-1]) + hRelSystUnc.SetBinContent(iBin+1, hSystUnc.GetBinContent(iBin+1)/hCFSelected.GetBinContent(iBin+1)) + hRelStatUnc.SetBinContent(iBin+1, hCFSelected.GetBinError(iBin+1)/hCFSelected.GetBinContent(iBin+1)) + hRatioStatSystUnc.SetBinContent(iBin+1, hCFSelected.GetBinError(iBin+1)/hSystUnc.GetBinContent(iBin+1)) + hCFWithSystUnc.SetBinContent(iBin+1, hCFSelected.GetBinContent(iBin+1)) + hCFWithSystUnc.SetBinError(iBin+1, hSystUnc.GetBinContent(iBin+1)**2) + hCFWithStatSystUnc.SetBinContent(iBin+1, hCFSelected.GetBinContent(iBin+1)) + hCFWithStatSystUnc.SetBinError(iBin+1, math.sqrt(hCFSelected.GetBinError(iBin+1)**2 + + hSystUnc.GetBinContent(iBin+1)**2 ) ) + + oFile.mkdir(f'{comb}/unc') + oFile.cd(f'{comb}/unc') + print('Writing histos with error informations ...') + hSystUnc.Write() + hRelSystUnc.Write() + hRelStatUnc.Write() + hRatioStatSystUnc.Write() + hCFWithSystUnc.Write() + hCFWithStatSystUnc.Write() + hResiduals.Write() + +oFile.mkdir('compare_residuals') +oFile.cd('compare_residuals') +for iComb in range(0, len(combs), 2): + hResiduals = TH2D(f"hRatioResiduals_{combs[iComb]}_{combs[iComb+1]}", + "Syst vars residuals", nBinsKStar, 0, uppEdgeKStar, + nSystVars, min(systVars), max(systVars)) + hResiduals.SetStats(0) + hResidualsPairOne = Load(oFile, f"{combs[iComb]}/unc/hResiduals") + hResidualsPairTwo = Load(oFile, f"{combs[iComb+1]}/unc/hResiduals") + for iBinX in range(nBinsKStar): + for iBinY in range(nSystVars): + resCombOne = hResidualsPairOne.GetBinContent(iBinX+1, iBinY+1) + resCombTwo = hResidualsPairTwo.GetBinContent(iBinX+1, iBinY+1) + if resCombOne!=0 and resCombOne!=0: + hResiduals.SetBinContent(iBinX+1, iBinY+1, resCombOne/resCombTwo) + else: + hResiduals.SetBinContent(iBinX+1, iBinY+1, 0) + hResiduals.Write() + +oFile.Close() +print(f'output saved in {oFileName}') diff --git a/cfgSystAnalysis.yml b/cfgSystAnalysis.yml new file mode 100644 index 0000000..840bd49 --- /dev/null +++ b/cfgSystAnalysis.yml @@ -0,0 +1,22 @@ +# where to save the output +odir: path/to/output/dir +suffix: suffix + +infileCFselected: /path/to/CF/used/in/data/analysis.root +binwidth: 4 # MeV/c +norm: [0.7, 0.8] # yield normalization with null + +infilessyst: ['/path/to/AnalysisResults1.root', + '/path/to/AnalysisResults2.root', + '/path/to/AnalysisResults3.root'] + +suffixessyst: [[var1_AnRes1.root, var2_AnRes1.root, var3_AnRes1.root], # each AnalysisResults can + [var1_AnRes2.root, var2_AnRes2.root, var3_AnRes2.root], # have a custom number + [var1_AnRes3.root, var2_AnRes3.root, var3_AnRes3.root]] # of variations + +gaussianfits: true # to evaluate the systematics uncertainties with gaussian + # fits to the bin entry distributions +systevalmaxbin: 150 # maximum kstar bin for systematics evaluation +systhistobins: 50 # nbins for the histo of the kstar bin entry variations +systhistointerval: 0.05 # range of the histo of the kstar bin entry variations, + # computed from the value of the selected CF for the bin \ No newline at end of file