Systematics cf

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}')

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}')

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