-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmakeCardInputWithUncs.py
136 lines (133 loc) · 6.64 KB
/
makeCardInputWithUncs.py
1
2
3
4
5
6
7
8
9
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
import numpy as np
import pandas as pd
def makeCardNoMB2(signalD2,signalB2,signalC2,signalA2,mcStats):
index = 0
labels = ["D2","B2","C2","A2"]
yields = [signalD2,signalB2,signalC2,signalA2]
mcLines =""
for im,mcStatUnc in enumerate(mcStats):
label = labels[im]
uncArr = ["-"]*8
if (mcStatUnc < 20):
uncName = "mcStat"+label+" lnN"
uncArr[index] = "%.2f"%(1+mcStatUnc/100.)
mcLines += uncName+"\t"+"\t".join(uncArr)+"\n"
else:
nEv = round((100./mcStatUnc)**2)
uncName = "mcStat"+label+" gmN "+str(nEv)
uncArr[index] = "%.4f"%(yields[im]/nEv)
mcLines += uncName+"\t"+"\t".join(uncArr)+"\n"
if index == 0:
index+=2
else:
index+=2
card="""
imax * number of channels
jmax * number of processes -1
kmax * number of nuisance parameters (sources of systematical uncertainties)
-------
bin D2 B2 C2 A2
observation 1 1 12 10
-------
bin D2 D2 B2 B2 C2 C2 A2 A2
process sig bkg sig bkg sig bkg sig bkg
process 0 1 0 1 0 1 0 1
rate {0} 1 {1} 1 {2} 1 {3} 1
-------
clusterSyst lnN 1.2 - - - - - - -
pileupSyst lnN 1.01 - - - - - - -
lumiSyst lnN 1.016 - - - - - - -
jesSyst lnN 1.06 - - - - - - -
xsecSyst lnN 0.923/1.046 - - - - - - -
pdfSyst lnN 1.032 - - - - - - -
ptSyst lnN 0.86/1.23 - - - - - - -
delta2 rateParam D2 bkg (@0*@1/@2) beta2,gamma2,alpha2
gamma2 rateParam C2 bkg 12 [0.01,25]
beta2 rateParam B2 bkg 1 [0.01,15]
alpha2 rateParam A2 bkg 10 [0.01,25]
""".format("%.3f"%signalD2,"%.3f"%signalB2,"%.3f"%signalC2,"%.3f"%signalA2)
card+=mcLines
return card
def makeCard(signalD1,signalB1,signalC1,signalA1,signalD2,signalB2,signalC2,signalA2,mcStats):
index = 0
labels = ["D1","B1","C1","A1","D2","B2","C2","A2"]
yields = [signalD1,signalB1,signalC1,signalA1,signalD2,signalB2,signalC2,signalA2]
mcLines =""
for im,mcStatUnc in enumerate(mcStats):
label = labels[im]
if mcStatUnc > 1 and mcStatUnc <= 100:
uncArr = ["-"]*17
if (mcStatUnc < 20):
uncName = "mcStat"+label+" lnN"
uncArr[index] = "%.2f"%(1+mcStatUnc/100.)
mcLines += uncName+"\t"+"\t".join(uncArr)+"\n"
else:
nEv = int(round((100./mcStatUnc)**2)+0.001)
uncName = "mcStat"+label+" gmN "+str(nEv)
uncArr[index] = "%.4f"%(yields[im]/nEv)
mcLines += uncName+"\t"+"\t".join(uncArr)+"\n"
if index == 0:
index+=3
else:
index+=2
card="""
imax * number of channels
jmax * number of processes -1
kmax * number of nuisance parameters (sources of systematical uncertainties)
-------
bin D1 B1 C1 A1 D2 B2 C2 A2
observation 5 1 20 32 1 1 12 10
-------
bin D1 D1 D1 B1 B1 C1 C1 A1 A1 D2 D2 B2 B2 C2 C2 A2 A2
process sig bkgPT bkg sig bkg sig bkg sig bkg sig bkg sig bkg sig bkg sig bkg
process 0 2 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
rate {0} 4.8 1 {1} 1 {2} 1 {3} 1 {4} 1 {5} 1 {6} 1 {7} 1
-------
clusterSyst lnN 1.2 - - - - - - - - 1.2 - - - - - - -
pileupSyst lnN 1.01 - - - - - - - - 1.01 - - - - - - -
lumiSyst lnN 1.016 - - - - - - - - 1.016 - - - - - - -
jesSyst lnN 1.06 - - - - - - - - 1.06 - - - - - - -
xsecSyst lnN 0.923/1.046 - - - - - - - - 0.923/1.046 - - - - - - -
pdfSyst lnN 1.032 - - - - - - - - 1.032 - - - - - - -
ptSyst lnN 0.86/1.23 - - - - - - - - 0.86/1.23 - - - - - - -
punchSyst lnN - 1.4 - - - - - - - - - - - - - - -
""".format("%.3f"%signalD1,"%.3f"%signalB1,"%.3f"%signalC1,"%.3f"%signalA1,"%.3f"%signalD2,"%.3f"%signalB2,"%.3f"%signalC2,"%.3f"%signalA2)
card+=mcLines
card += """
delta1 rateParam D1 bkg (@0*@1/@2) beta1,gamma1,alpha1
gamma1 rateParam C1 bkg 20 [1,40]
beta1 rateParam B1 bkg 1 [0.001,6]
alpha1 rateParam A1 bkg 32 [5,60]
delta2 rateParam D2 bkg (@0*@1/@2) beta2,gamma2,alpha2
gamma2 rateParam C2 bkg 12 [0.01,25]
beta2 rateParam B2 bkg 1 [0.01,15]
alpha2 rateParam A2 bkg 10 [0.01,25]
"""
return card
signalYields = pd.read_csv("signalYieldsWithUncsHV.csv",delimiter=",")
import os
outputDir = "limitsSignalUncsHVUpdatePTNoMB2/cards"
noMB2 = True
if not os.path.exists(outputDir):
os.makedirs(outputDir)
for _,signal in signalYields.iterrows():
if "portal" not in list(signalYields.columns):
# card = makeCardNoMB2(signal.sr1,signal.sr2)
mcStats = [signal.mc_uncSR1,signal.mc_uncB1,signal.mc_uncC1,signal.mc_uncA1,
signal.mc_uncSR2,signal.mc_uncB2,signal.mc_uncC2,signal.mc_uncA2]
if noMB2:
card = makeCardNoMB2(signal.sr2,signal.B2,signal.C2,signal.A2,mcStats[4:])
else:
card = makeCard(signal.sr1,signal.B1,signal.C1,signal.A1,signal.sr2,signal.B2,signal.C2,signal.A2,mcStats)
with open(outputDir+"/card_TH_"+str(int(signal.mass))+"_"+str(int(signal.lifetime))+".txt","w") as f:
f.write(card)
else:
#DBCA
mcStats = [signal.mc_unc_SR1,signal.mc_unc_B1,signal.mc_unc_C1,signal.mc_unc_A1,
signal.mc_unc_SR2,signal.mc_unc_B2,signal.mc_unc_C2,signal.mc_unc_A2]
if noMB2:
card = makeCardNoMB2(signal.SR2,signal.B2,signal.C2,signal.A2,mcStats[4:])
else:
card = makeCard(signal.SR1,signal.B1,signal.C1,signal.A1,signal.SR2,signal.B2,signal.C2,signal.A2,mcStats)
with open(outputDir+"/card_HV_"+signal.portal+"_"+str(int(signal.mX))+"_"+str(int(signal.ctau))+"_"+str(signal.xiO).replace(".","p")+"_"+str(signal.xiL).replace(".","p")+".txt","w") as f:
f.write(card)