plot_xsec_1d.py

Go to the documentation of this file.

#!/usr/bin/env python

from __future__ import print_function
from ROOT import *
from bisect import bisect
import argparse
import math
import os


# global variables

# file directory
infp = './'


syst_names = [
    'lightdown','lightup','ckv','calibneg','calibpos','calibshape','MuES_onlyMuKEDw','MuES_onlyMuKEUp',
    '2kA_Dw','2kA_Up','02mmBeamSpotSize_Dw','02mmBeamSpotSize_Up','1mmBeamShiftX_Dw','1mmBeamShiftX_Up',
    '1mmBeamShiftY_Dw','1mmBeamShiftY_Up','3mmHorn1X_Dw','3mmHorn1X_Up','3mmHorn1Y_Dw','3mmHorn1Y_Up',
    '3mmHorn2X_Dw','3mmHorn2X_Up','3mmHorn2Y_Dw','3mmHorn2Y_Up','7mmTargetZ_Dw','7mmTargetZ_Up',
    'MagneticFieldinDecayPipe_Dw','MagneticFieldinDecayPipe_Up','1mmHornWater_Dw','1mmHornWater_Up'
]

# names of the fractional uncertainties
fu_names = ['Flux_HP','GENIE','Light','Calibration','Cherenkov','CalibShape','MuEScale','Normalization','Focusing']
final_fu = ['Flux (HP)','GENIE','Light','Calibration','Cherenkov','Calib. Shape','Mu Energy Scale','Normalization','Focusing']
col_fe = [2,4,2,4,8,8,kOrange-3,kOrange-3,1]
col_st = [1,1,2,2,1,2,1,2,2]

def get_genie_universes(basefn):

    hl = []
    fn = os.path.join(infp, basefn)
    fgenie.append(TFile(fn))
    # find the genie universe IDs from filename
    sIdx = int(basefn.split('_')[1])
    eIdx = int(basefn.split('_')[2])+1
    for i in range(sIdx,eIdx):
        hl.append(fgenie[len(fgenie)-1].Get('xs_nominal/hXS1Dsec_genie{}'.format(i)))

    return hl

def get_ppfx_universes(fppfx_other):

    hl = []
    # by default we have 100 ppfx universes
    for i in range(0,100):
        hl.append(fppfx_other.Get('xs_nominal/hXS1Dsec_ppfx{}'.format(i)))
    
    return hl


def do_multiverse_1bin(isys, univ_vals):

    sorted_list = sorted(univ_vals)
    cv_pos = bisect(sorted_list, cv)
    nlist = len(univ_vals)
    up_pos = nlist-1 if cv_pos >= nlist else int((nlist - cv_pos)*.68 + cv_pos)
    dw_pos = 0 if cv_pos <= 0 else int(cv_pos - cv_pos*.68)

    up_shift = sorted_list[up_pos] - cv
    dw_shift = cv - sorted_list[dw_pos]

    lhup[isys].SetBinContent(b, up_shift)
    lhdw[isys].SetBinContent(b, dw_shift)
    fu = 0 if cv <= 0 else max(up_shift, dw_shift)/cv
    lhfu[isys].SetBinContent(b, fu)


def do_focusing_1bin(isys, sstart, nsys):
    
    up_shifts = [0 for k in range(nsys)]
    dw_shifts = [0 for k in range(nsys)]
    for k in range(nsys):
        dwshift = cv - hother[sstart+k*2].GetBinContent(i)
        upshift = hother[sstart+k*2+1].GetBinContent(i) - cv
        
        if upshift > 0 and dwshift > 0:
            up_shifts[k] = upshift
            dw_shifts[k] = dwshift
        elif upshift <= 0 and dwshift <= 0: # up down flip
            up_shifts[k] = -dwshift
            dw_shifts[k] = -upshift
        elif upshift > 0 and dwshift <= 0:
            up_shifts[k] = upshift
        else:
            dw_shifts[k] = dwshift
    
    up_shift = math.sqrt(sum([x**2 for x in up_shifts]))
    dw_shift = math.sqrt(sum([x**2 for x in dw_shifts]))
    lhup[isys].SetBinContent(b, up_shift)
    lhdw[isys].SetBinContent(b, dw_shift)

    fu = 0 if cv <= 0 else max(up_shift, dw_shift)/cv
    lhfu[isys].SetBinContent(b, fu)


def do_up_down_1bin(isys, upidx, dwidx):
    '''
    It could happen that in some regions both systematics go the same direction.
    In that case we take the larger value for that direction and asign zero to the other.
    '''

    # initialize
    lhup[isys].SetBinContent(b, 0)
    lhdw[isys].SetBinContent(b, 0)

    upshift = hother[upidx].GetBinContent(b)-cv
    dwshift = cv-hother[dwidx].GetBinContent(b)
    larger = max(abs(upshift), abs(dwshift))

    if upshift > 0 and dwshift > 0:
        lhup[isys].SetBinContent(b, upshift)
        lhdw[isys].SetBinContent(b, dwshift)
    elif upshift <= 0 and dwshift <= 0: # up down flip
        lhup[isys].SetBinContent(b, -dwshift)
        lhdw[isys].SetBinContent(b, -upshift)
    elif upshift > 0 and dwshift <= 0:
        lhup[isys].SetBinContent(b, upshift)
    else:
        lhdw[isys].SetBinContent(b, dwshift)
    
    fu = 0 if cv <= 0 else larger/cv
    lhfu[isys].SetBinContent(b, fu)


def do_single_shift_1bin(isys, idx):
    '''
    Each bin is initialized to 0 to allow one-sided error bars.
    '''

    lhup[isys].SetBinContent(b, 0)
    lhdw[isys].SetBinContent(b, 0)
    shift = hother[idx].GetBinContent(b)-cv
    if shift > 0:
        lhup[isys].SetBinContent(b, shift)
    else:
        lhdw[isys].SetBinContent(b, -shift)

    fu = 0 if cv <= 0 else abs(shift)/cv
    lhfu[isys].SetBinContent(b, fu)


def setUncer1D(hIn, ihist):
    hIn.SetLineColor(col_fe[ihist])
    hIn.SetLineWidth(3)
    hIn.SetLineStyle(col_st[ihist])
    hIn.SetStats(0)


if __name__ == '__main__':

    # command line arguments
    # fake data / real data switch
    parser = argparse.ArgumentParser()
    parser.add_argument('-d', '--data', action='store_true', help='Switch for using data instead of fake data.')
    args = parser.parse_args()
    use_data = args.data

    # load files into file containers
    # mc is the Monte Carlo histogram in the final plot
    fgenie = []
    fmc = TFile(os.path.join(infp, 'nominal_xsec.root'))
    # Here I load (fake) data and all other systematics.
    # It is an unforunate design to put data together with other systematics.
    # However it is too late. Get over!
    if not use_data:
        fppfx_other = TFile(os.path.join(infp, 'ppfx_others_v2.root'))
    else:
        fppfx_other = TFile(os.path.join(infp, 'ppfx_others_data.root'))
    
    # historam containers
    hgenie = []
    hppfx = []
    hother = []

    # get nominal cross section
    hmc = fmc.Get('hXS1Dsec_nominal')
    # get (fake) data cross section
    hdt = fppfx_other.Get('xs_nominal/hXS1Dsec_cv')

    # store genie histograms
    if not use_data:
        hgenie += get_genie_universes('genie_0_99_v2.root')
        hgenie += get_genie_universes('genie_100_299_v2.root')
        hgenie += get_genie_universes('genie_300_499_v2.root')
    else:
        hgenie += get_genie_universes('genie_0_99_data.root')
        hgenie += get_genie_universes('genie_100_299_data.root')
        hgenie += get_genie_universes('genie_300_499_data.root')
    
    # store ppfx histograms
    hppfx += get_ppfx_universes(fppfx_other)

    # store all other systematics
    for syst_name in syst_names:
        hother.append(fppfx_other.Get('xs_nominal/hXS1Dsec_{}'.format(syst_name)))
    
    # histogram containers: one for each systematics
    nfu = len(fu_names)
    lhup = [hmc.Clone('hUp_{}'.format(i)) for i in range(nfu)]
    lhdw = [hmc.Clone('hDw_{}'.format(i)) for i in range(nfu)]
    lhfu = [hmc.Clone('hFU_{}'.format(i)) for i in range(nfu)]

    # number of bins to loop through
    Nbins = hdt.GetNbinsX()
    
    # loop through each bin
    for b in range(1,Nbins+1):
        
        # central value of nominal MC
        cv = hdt.GetBinContent(b)

        # ppfx
        do_multiverse_1bin(0, [h.GetBinContent(b) for h in hppfx])

        # genie
        do_multiverse_1bin(1, [h.GetBinContent(b) for h in hgenie])

        # light
        do_up_down_1bin(2,1,0)
        

        # calibration
        do_up_down_1bin(3,4,3)

        # Cherenkov
        do_single_shift_1bin(4,2)

        # calibshape
        do_single_shift_1bin(5,5)

        # muon energy scale
        do_up_down_1bin(6,7,6)

        # Normalization (TN, doc-):
        # POT: 0.28% (Section 9.1)
        # Intensity: 2% (Section 9.2)
        # Particle propagation: 1-2%, so we take 2% (Section 9.5).
        # fe_norm = sqrt(pow(0.0028,2)+pow(0.02,2)+pow(0.02,2));
        fe_norm = 0.028422526
        lhup[7].SetBinContent(b, abs(cv*fe_norm))
        lhdw[7].SetBinContent(b, abs(cv*fe_norm))
        fu = 0 if cv <= 0 else fe_norm
        lhfu[7].SetBinContent(b, fu)

        # focusing
        # first argument starting index, second number of up-down groups
        do_focusing_1bin(8,8,11)


    # make sure the output directory exists
    out_subdir = 'real_data' if use_data else 'fake_data'
    dest_path = os.path.join('output1d', out_subdir)
    if not os.path.exists(dest_path):
        os.makedirs(dest_path)


    # ia = 8
    # ib = 9
    # lhup = lhup[ia:ib]
    # lhdw = lhdw[ia:ib]
    # # plot data
    # gdt = TGraphAsymmErrors()
    # set up output graph points
    # for i in range(1,Nbins+1):
    #     # sum up and down errors in quadrature
    #     qsum_up = math.sqrt(sum([h.GetBinContent(i)**2 for h in lhup]))
    #     qsum_dw = math.sqrt(sum([h.GetBinContent(i)**2 for h in lhdw]))

    #     binw  = hdt.GetXaxis().GetBinWidth(i)
    #     gdt.SetPoint(i, hdt.GetXaxis().GetBinCenter(i), hdt.GetBinContent(i))
    #     gdt.SetPointError(i, binw/2, binw/2, qsum_dw, qsum_up)
    
    # gdt.Draw('AP1')

    # plot fractional uncertainty
    cfu = TCanvas()
    lfu = TLegend(0.2,0.5,0.45,0.88)
    lhfu[0].GetXaxis().SetNdivisions(510)
    lfu.SetBorderSize(0)
    for i in range(len(lhfu)):
        setUncer1D(lhfu[i],i)
        if i == 0:
            lhfu[i].Draw('hist')
            lhfu[i].SetMaximum(1)
            lhfu[i].GetXaxis().CenterTitle()
            lhfu[i].GetXaxis().SetTitle('Neutrino Energy (GeV)')
            lhfu[i].GetYaxis().CenterTitle()
            lhfu[i].GetYaxis().SetTitle('Fractional Uncertainty')
        else:
            lhfu[i].Draw('histsame')
        lfu.AddEntry(lhfu[i],final_fu[i],'L')
    lfu.Draw()
    # save fractional uncertainty to file
    cfu.SaveAs('{}/fe1D.pdf'.format(dest_path))

    raw_input('Press enter to continue.')