-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrunHZZFiducialXS.py
1662 lines (1417 loc) · 117 KB
/
runHZZFiducialXS.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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/python
#-----------------------------------------------
# Latest update: 2014.10.16
#-----------------------------------------------
import ast
import json
import math
import optparse
import os
import sys
from cgi import test
from decimal import *
import yaml
# INFO: Following items are imported from either python directory or Inputs
from createXSworkspace import createXSworkspace
from higgs_xsbr_13TeV import filtereff, higgs4l_br, higgs_xs, higgsZZ_br
from Input_Info import combineOutputs, datacardInputs, channels_central
from read_bins import read_bins
from ROOT import *
from sample_shortnames import background_samples, sample_shortnames
from Utils import (GetDirectory, get_linenumber, logger, logging,
mergeDictionary_average, processCmd)
jes_sources = [
"Abs",
"Abs_year",
"BBEC1",
"BBEC1_year",
"EC2",
"EC2_year",
"FlavQCD",
"HF",
"HF_year",
"RelBal",
"RelSample_year", ]
### Define function for parsing options
def parseOptions():
global opt, args, runAllSteps, unblindString
usage = ('usage: %prog [options]\n'
+ '%prog -h for help')
parser = optparse.OptionParser(usage)
# input options
parser.add_option('-d', '--dir', dest='SOURCEDIR',type='string',default='./', help='run from the SOURCEDIR as working area, skip if SOURCEDIR is an empty string')
parser.add_option('', '--asimovModelName',dest='ASIMOVMODEL',type='string',default='SM_125', help='Name of the Asimov Model')
parser.add_option('', '--asimovMass',dest='ASIMOVMASS',type='string',default='125.38', help='Asimov Mass')
parser.add_option('', '--modelNames',dest='MODELNAMES',type='string',default="SM_125",help='Names of models for unfolding, separated by , (comma) like "SM_125,SMup_125,SMdn_125". Default is "SM_125"')
# FIXME: `FIXMASS` option should be bool. As per its name. No?
parser.add_option('', '--fixMass', dest='FIXMASS', type='string',default='125.38', help='Fix mass, default is a string "125.38" or can be changed to another string, e.g."125.6" or "False"')
parser.add_option('', '--obsName', dest='OBSNAME', type='string',default='', help='Name of the observable, supported: "inclusive", "pT4l", "eta4l", "massZ2", "nJets"')
parser.add_option('', '--obsBins', dest='OBSBINS', type='string',default='', help='Bin boundaries for the diff. measurement separated by "|", e.g. as "|0|50|100|", use the defalut if empty string')
parser.add_option('', '--fixFrac', action='store_true', dest='FIXFRAC', default=False, help='fix the fractions of 4e and 4mu when extracting the results, default is False')
# action options - "redo"
parser.add_option('', '--redoEff', action='store_true', dest='redoEff', default=False, help='Redo the eff. factors, default is False')
parser.add_option('', '--redoTemplates', action='store_true', dest='redoTemplates',default=False, help='Redo the bkg shapes and fractions, default is False')
# action options - "only"
parser.add_option('', '--effOnly', action='store_true', dest='effOnly', default=False, help='Extract the eff. factors only, default is False')
parser.add_option('', '--templatesOnly', action='store_true', dest='templatesOnly', default=False, help='Prepare the bkg shapes and fractions only, default is False')
parser.add_option('', '--uncertOnly', action='store_true', dest='uncertOnly', default=False, help='Extract the uncertanties only, default is False')
parser.add_option('', '--resultsOnly', action='store_true', dest='resultsOnly', default=False, help='Run the measurement only, default is False')
parser.add_option('', '--finalplotsOnly',action='store_true', dest='finalplotsOnly',default=False, help='Make the final plots only, default is False')
# action options - do Z4l measurement
parser.add_option('', '--doZ4l', action='store_true', dest='doZ4l', default=False, help='Perform the Z->4l measurement instead of H->4l, default is False')
parser.add_option('', '--doRatio', action='store_true', dest='doRatio', default=False, help='Do H4l/Z4l ratio, default is False')
# Unblind option
parser.add_option('', '--unblind', action='store_true', dest='UNBLIND', default=False, help='Use real data')
# Calculate Systematic Uncertainties
parser.add_option('', '--calcSys', action='store_true', dest='SYS', default=False, help='Calculate Systematic Uncertainties (in addition to stat+sys)')
parser.add_option('', '--lumiscale', type='string', dest='LUMISCALE', default='1.0', help='Scale yields')
parser.add_option('-i', '--inYAMLFile', dest='inYAMLFile', type='string', default="Inputs/observables_list.yml", help='Input YAML file having observable names and bin information')
parser.add_option("-l", "--logLevel", dest="logLevel", type = 'string', default = '2', help="Change log verbosity(WARNING: 0, INFO: 1, DEBUG: 2, ERROR: 3)")
parser.add_option('-y', '--year', dest="ERA", type = 'string', default = '2018', help='Specifies the data taking period')
# store options and arguments as global variables
global opt, args
(opt, args) = parser.parse_args()
# NOTE: append the directory `datacardInputs`, as .py files inside is going to load using import.
# load XS-specific modules
GetDirectory(datacardInputs.format(year = opt.ERA))
sys.path.append('./'+datacardInputs.format(year = opt.ERA))
log_level = logging.DEBUG # default initialization
if opt.logLevel == "0":
log_level = logging.WARNING
elif opt.logLevel == "1":
log_level = logging.INFO
elif opt.logLevel == "2":
log_level = logging.DEBUG
elif opt.logLevel == "3":
log_level = logging.ERROR
logger.setLevel( log_level)
unblindString = ""
if (opt.UNBLIND): unblindString = "_unblind"
# prepare the global flag if all the step should be run
runAllSteps = not(opt.effOnly or opt.templatesOnly or opt.uncertOnly or opt.resultsOnly or opt.finalplotsOnly)
# runAllSteps = not(opt.templatesOnly or opt.uncertOnly or opt.resultsOnly or opt.finalplotsOnly)
if (opt.OBSBINS=='' and opt.OBSNAME!='inclusive'):
parser.error('Bin boundaries not specified for differential measurement. Exiting...')
sys.exit()
# FIXME: Check why we need these directories
# The directory `combineOutputs` is used to keep the generated workspaces
GetDirectory(combineOutputs.format(year = opt.ERA))
dirToExist = ['templates', datacardInputs.format(year = opt.ERA), combineOutputs.format(year = opt.ERA)]
for dir in dirToExist:
if not os.path.isdir(os.getcwd()+'/'+dir+'/'):
parser.error(os.getcwd()+'/'+dir+'/ is not a directory. Exiting...')
sys.exit()
### Add v3 model specific commands
def runv3(years, observableBins, obsName, fitName, physicalModel, JES = False, fStates=['4e', '4mu', '2e2mu']):
card_name = 'hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.txt'
cmd_combCards = 'combineCards.py '
nBins = len(observableBins) - 1
print ("obs Bins: ")
print(observableBins)
#if not doubleDiff: nBins = nBins-1 VM: Add the option for double diff measurement
for year in years:
for cat in fStates:
for i in range(nBins):
if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName and not obsName == 'njets_pt30_eta4p7':
low = str(observableBins[i][0]).replace('.','p').replace('-','m')
high = str(observableBins[i][1]).replace('.','p').replace('-','m')
low_2nd = str(observableBins[i][2]).replace('.','p').replace('-','m')
high_2nd = str(observableBins[i][3]).replace('.','p').replace('-','m')
boundaries = low+'_'+high+'_'+low_2nd+'_'+high_2nd
else:
low = str(observableBins[i]).replace('.','p').replace('-','m')
high = str(observableBins[i+1]).replace('.','p').replace('-','m')
boundaries = low+'_'+high
if float(observableBins[i+1]) > 1000:
boundaries = 'GT'+str(int(observableBins[i]))
dc_name = 'xs_125.0_%s/hzz4l_%sS_13TeV_xs_%s_bin%d_v3.txt ' %(year,cat,obsName,i)
cmd_combCards += 'hzz_%s_%s_cat%s_%s=%s' %(fitName,boundaries,cat,year,dc_name)
cmd_combCards += '> %s' %card_name
cmd_addNuis = ''
if opt.ERA == 'allYear':
if obsName == 'mass4l_zzfloating': # Remove bkg theo nuisances in case of zz floating
cmd_addNuis = 'echo "nuis group = CMS_eff_e CMS_eff_m CMS_hzz2e2mu_Zjets_2016 CMS_hzz2e2mu_Zjets_2017 CMS_hzz2e2mu_Zjets_2018 CMS_hzz4e_Zjets_2016 CMS_hzz4e_Zjets_2017 CMS_hzz4e_Zjets_2018 CMS_hzz4mu_Zjets_2016 CMS_hzz4mu_Zjets_2017 CMS_hzz4mu_Zjets_2018 lumi_13TeV_2016_uncorrelated lumi_13TeV_2017_uncorrelated lumi_13TeV_2018_uncorrelated lumi_13TeV_correlated_16_17_18 lumi_13TeV_correlated_17_18 CMS_zz4l_sigma_e_sig CMS_zz4l_sigma_m_sig CMS_zz4l_n_sig_3_2016 CMS_zz4l_n_sig_3_2017 CMS_zz4l_n_sig_3_2018 CMS_zz4l_n_sig_2_2016 CMS_zz4l_n_sig_2_2017 CMS_zz4l_n_sig_2_2018 CMS_zz4l_n_sig_1_2016 CMS_zz4l_n_sig_1_2017 CMS_zz4l_n_sig_1_2018'
else:
cmd_addNuis = 'echo "nuis group = CMS_eff_e CMS_eff_m CMS_hzz2e2mu_Zjets_2016 CMS_hzz2e2mu_Zjets_2017 CMS_hzz2e2mu_Zjets_2018 CMS_hzz4e_Zjets_2016 CMS_hzz4e_Zjets_2017 CMS_hzz4e_Zjets_2018 CMS_hzz4mu_Zjets_2016 CMS_hzz4mu_Zjets_2017 CMS_hzz4mu_Zjets_2018 QCDscale_VV QCDscale_ggVV kfactor_ggzz lumi_13TeV_2016_uncorrelated lumi_13TeV_2017_uncorrelated lumi_13TeV_2018_uncorrelated lumi_13TeV_correlated_16_17_18 lumi_13TeV_correlated_17_18 pdf_gg pdf_qqbar CMS_zz4l_sigma_e_sig CMS_zz4l_sigma_m_sig CMS_zz4l_n_sig_3_2016 CMS_zz4l_n_sig_3_2017 CMS_zz4l_n_sig_3_2018 CMS_zz4l_n_sig_2_2016 CMS_zz4l_n_sig_2_2017 CMS_zz4l_n_sig_2_2018 CMS_zz4l_n_sig_1_2016 CMS_zz4l_n_sig_1_2017 CMS_zz4l_n_sig_1_2018'
if JES:
cmd_addNuis += ' CMS_scale_j_Abs CMS_scale_j_Abs_2016 CMS_scale_j_BBEC1 CMS_scale_j_BBEC1_2016 CMS_scale_j_EC2 CMS_scale_j_EC2_2016 CMS_scale_j_FlavQCD CMS_scale_j_HF CMS_scale_j_HF_2016 CMS_scale_j_RelBal CMS_scale_j_RelSample_2016 CMS_scale_j_Abs CMS_scale_j_Abs_2017 CMS_scale_j_BBEC1 CMS_scale_j_BBEC1_2017 CMS_scale_j_EC2 CMS_scale_j_EC2_2017 CMS_scale_j_FlavQCD CMS_scale_j_HF CMS_scale_j_HF_2017 CMS_scale_j_RelBal CMS_scale_j_RelSample_2017 CMS_scale_j_Abs CMS_scale_j_Abs_2018 CMS_scale_j_BBEC1 CMS_scale_j_BBEC1_2018 CMS_scale_j_EC2 CMS_scale_j_EC2_2018 CMS_scale_j_FlavQCD CMS_scale_j_HF CMS_scale_j_HF_2018 CMS_scale_j_RelBal CMS_scale_j_RelSample_2018'
else:
if obsName == 'mass4l_zzfloating':
cmd_addNuis = 'echo "nuis group = CMS_eff_e CMS_eff_m CMS_hzz2e2mu_Zjets_'+str(opt.ERA)+' CMS_hzz4e_Zjets_'+str(opt.ERA)+' CMS_hzz4mu_Zjets_'+str(opt.ERA)+' lumi_13TeV_'+str(opt.ERA)+' CMS_zz4l_sigma_e_sig CMS_zz4l_sigma_m_sig CMS_zz4l_n_sig_3_'+str(opt.ERA)+' CMS_zz4l_mean_e_sig CMS_zz4l_mean_m_sig'
else:
cmd_addNuis = 'echo "nuis group = CMS_eff_e CMS_eff_m CMS_hzz2e2mu_Zjets_'+str(opt.ERA)+' CMS_hzz4e_Zjets_'+str(opt.ERA)+' CMS_hzz4mu_Zjets_'+str(opt.ERA)+' QCDscale_VV QCDscale_ggVV kfactor_ggzz lumi_13TeV_'+str(opt.ERA)+' pdf_gg pdf_qqbar CMS_zz4l_sigma_e_sig CMS_zz4l_sigma_m_sig CMS_zz4l_n_sig_3_'+str(opt.ERA)+' CMS_zz4l_mean_e_sig CMS_zz4l_mean_m_sig'
if JES:
cmd_addNuis += ' CMS_scale_j_Abs CMS_scale_j_Abs_'+str(opt.ERA)+' CMS_scale_j_BBEC1 CMS_scale_j_BBEC1_'+str(opt.ERA)+' CMS_scale_j_EC2 CMS_scale_j_EC2_'+str(opt.ERA)+' CMS_scale_j_FlavQCD CMS_scale_j_HF CMS_scale_j_HF_'+str(opt.ERA)+' CMS_scale_j_RelBal CMS_scale_j_RelSample_'+str(opt.ERA)
cmd_addNuis += '" >> hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.txt'
processCmd(cmd_combCards)
processCmd(cmd_addNuis)
cmd_t2w = 'text2workspace.py %s -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel --PO verbose ' %card_name
cmd_t2w += "--PO 'higgsMassRange=123,127' "
for i in range(nBins):
if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName and not obsName == 'njets_pt30_eta4p7':
low = str(observableBins[i][0]).replace('.','p').replace('-','m')
high = str(observableBins[i][1]).replace('.','p').replace('-','m')
low_2nd = str(observableBins[i][2]).replace('.','p').replace('-','m')
high_2nd = str(observableBins[i][3]).replace('.','p').replace('-','m')
boundaries = low+'_'+high+'_'+low_2nd+'_'+high_2nd
else:
low = str(observableBins[i]).replace('.','p').replace('-','m')
high = str(observableBins[i+1]).replace('.','p').replace('-','m')
boundaries = low+'_'+high
if float(observableBins[i+1]) > 1000:
boundaries = 'GT'+str(int(observableBins[i]))
#process = 'smH_%s_%s' %(fitName, boundaries) VM testing
process = 'trueHBin%s' %(i)
#process = 'trueH'
POI = 'r_smH_%s_%d' %(fitName, i)
POI_n = 'r_smH_%d' %i
cmd_t2w += "--PO 'map=.*/%s:%s[1.0,0.0,3.0]' " %(process, POI)
print(cmd_t2w)
processCmd(cmd_t2w)
cmd = 'cp hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root ./combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root'
print cmd, '\n'
processCmd(cmd,1)
os.chdir('./combine_files/')
for i in range(nBins):
if obsName == 'dphijj' and i == 4:
upScanRange = 5
nPoints = 200
else:
upScanRange = 2.5
nPoints = 200
POI = 'r_smH_%s_%d' %(fitName, i)
POI_n = 'r_smH_%d' %i
cmd_fit = 'combine -n _%s_%s -M MultiDimFit %s ' %(obsName, POI_n, 'SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root')
cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points='+str(nPoints)+' --cminDefaultMinimizerStrategy 0 '
if not opt.UNBLIND: cmd_fit += '-t -1 --saveToys --setParameters %s=1 ' %(POI)
cmd_fit_tmp = cmd_fit + '-P %s --setParameterRanges %s=0.0,%i --redefineSignalPOI %s' %(POI, POI, upScanRange, POI)
print(cmd_fit_tmp)
processCmd(cmd_fit_tmp)
if obsName == 'mass4l_zzfloating':
for i in range(nBins):
POI = 'zz_norm_%d' %i
POI_xs = 'r_smH_%s_%d' %(fitName, i)
POI_n = 'r_smH_%d' %i
cmd_fit = 'combine -n _%s_zz_norm_0 -M MultiDimFit %s ' %(obsName, 'SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root')
cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --cminDefaultMinimizerStrategy 0 '
if not opt.UNBLIND: cmd_fit += '-t -1 --saveToys --setParameters %s=1 ' %(POI_xs)
cmd_fit_tmp = cmd_fit + '-P %s --redefineSignalPOI %s' %(POI, POI)
print(cmd_fit_tmp)
processCmd(cmd_fit_tmp)
#Stat-only
for i in range(nBins):
if obsName == 'dphijj' and i == 4:
upScanRange = 5
nPoints = 200
else:
upScanRange = 2.5
nPoints = 200
POI = 'r_smH_%s_%d' %(fitName, i)
POI_n = 'r_smH_%d' %i
cmd_fit = 'combine -n _%s_%s_NoSys -M MultiDimFit %s' %(obsName, POI_n, 'higgsCombine_'+obsName+'_'+POI_n+'.MultiDimFit.mH125.38')
if not opt.UNBLIND: cmd_fit = cmd_fit + '.123456'
cmd_fit += '.root -w w --snapshotName "MultiDimFit" -m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points='+str(nPoints)+' --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 '
if not opt.UNBLIND: cmd_fit += '-t -1 --saveToys --setParameters %s=1 ' %(POI)
cmd_fit_tmp = cmd_fit + '-P %s --setParameterRanges %s=0.0,%i --redefineSignalPOI %s' %(POI, POI, upScanRange, POI)
print cmd_fit_tmp
processCmd(cmd_fit_tmp)
if obsName == 'mass4l_zzfloating':
for i in range(nBins):
POI = 'zz_norm_%d' %i
POI_xs = 'r_smH_%s_%d' %(fitName, i)
POI_n = 'zz_norm_%d' %i
cmd_fit = 'combine -n _%s_zz_norm_0_NoSys -M MultiDimFit %s' %(obsName, 'higgsCombine_'+obsName+'_'+POI_n+'.MultiDimFit.mH125.38')
if not opt.UNBLIND: cmd_fit = cmd_fit + '.123456'
cmd_fit += '.root -w w --snapshotName "MultiDimFit" -m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 '
if not opt.UNBLIND: cmd_fit += '-t -1 --saveToys --setParameters %s=1 ' %(POI_xs)
cmd_fit_tmp = cmd_fit + '-P %s --redefineSignalPOI %s' %(POI, POI)
print(cmd_fit_tmp)
processCmd(cmd_fit_tmp)
### Extract the all efficiency factors (inclusive/differential, all bins, all final states)
def extractFiducialEfficiencies(obsName, observableBins, modelName):
"""Extract efficiencies and plot the 2D signal efficiency
Args:
obsName (str): name of the observable
observableBins (array): Array containing the bin boundaries
modelName (str): Name of model. For example SM_125, SMup_125, etc.
"""
#from inputs_bkg_{obsName} import fractionsBackground and observableBins
if (not opt.redoEff):
print ('[Skipping eff. and out.factors for '+str(obsName)+']')
return
print ('[Extracting eff. and out.factors]')
#cmd = 'python efficiencyFactors.py --dir='+opt.SOURCEDIR+' --obsName='+opt.OBSNAME+' --obsBins="'+opt.OBSBINS+'" -l -q -b'
cmd = 'python efficiencyFactors.py --dir='+opt.SOURCEDIR+' --obsName='+opt.OBSNAME+' --obsBins="'+opt.OBSBINS+'" -l -q -b --doPlots --doFit'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
print (output)
# if (not opt.OBSNAME.startswith("mass4l")):
if (not (opt.OBSNAME=="mass4l")):
cmd = 'python plot2dsigeffs.py -l -q -b --obsName='+opt.OBSNAME+' --obsBins="'+opt.OBSBINS+'"'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
### Extract the templates for given obs, for all bins and final states (differential)
def extractBackgroundTemplatesAndFractions(obsName, observableBins, year, obs_ifJES, obs_ifJES2):
global opt
logger.debug("[INFO] Obs Name: {:15}\tBins: {}".format(obsName, observableBins))
fractionBkg = {}; lambdajesdnBkg={}; lambdajesupBkg={}
if os.path.isfile(datacardInputs.format(year = year)+'/inputs_bkg_'+{0:'',1:'z4l_'}[int(opt.doZ4l)]+obsName.replace(" ","_")+'.py'):
_temp = __import__('inputs_bkg_'+{0:'',1:'z4l_'}[int(opt.doZ4l)]+obsName.replace(" ","_"), globals(), locals(), ['observableBins','fractionsBackground','lambdajesupBkg','lambdajesdnBkg'], -1)
if (hasattr(_temp,'observableBins') and _temp.observableBins == observableBins and not opt.redoTemplates):
logger.info ('[Fractions already exist for the given binning. Skipping templates/shapes... ]')
return
if (hasattr(_temp,'fractionsBackground') and hasattr(_temp,'lambdajesupBkg') and hasattr(_temp,'lambdajesdnBkg')):
fractionBkg = _temp.fractionsBackground
lambdajesupBkg = _temp.lambdajesupBkg
lambdajesdnBkg = _temp.lambdajesdnBkg
logger.debug('Preparing bkg shapes and fractions, for bins with boundaries {}'.format(observableBins))
# save/create/prepare directories and compile templates script
# FIXME: directory name hardcoded
currentDir = os.getcwd(); os.chdir('./templates/')
# Here, rm command is mandatory, just to ensure that make works. If make command files
# then this can pick older executable. So, to avoid this we first delete the executable
# logger.info("==> Remove the executable and Compile the package main_fiducialXSTemplates...")
# cmd = 'rm main_fiducialXSTemplates; make';
# processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# moved the compilation part in the "setup.sh"
# FIXME: this directory name is hardcoded in fiducialXSTemplates.C
# FIXME: Try to link the two automatically.
# FIXME: Also, this name is passed as one of arguments of `main_fiducialXSTemplates()`
DirectoryToCreate = 'templatesXS_'+str(year)+'/DTreeXS_' + obsName.replace(" ","_") + '/13TeV/'
logger.debug('Create directory: {}'.format(DirectoryToCreate))
logger.debug('compile the script inside the template directory')
cmd = 'mkdir -p '+DirectoryToCreate; processCmd(cmd, get_linenumber(), os.path.basename(__file__), 1)
# extract bkg templates and bin fractions
sZZname2e2mu = 'ZZTo2e2mu_powheg'
sZZname4mu = 'ZZTo4mu_powheg'
sZZname4e = 'ZZTo4e_powheg'
if (opt.doZ4l):
"""If Z->4l analysis is set true
"""
sZZname2e2mu = 'ZZTo2e2mu_powheg_tchan'
sZZname4mu = 'ZZTo4mu_powheg_tchan'
sZZname4e = 'ZZTo4e_powheg_tchan'
# FIXME: Explain each name for documentation
bkg_sample_tags = ['ZX4l_CR', 'ZX4l_CR_4e', 'ZX4l_CR_4mu', 'ZX4l_CR_2e2mu', sZZname2e2mu, sZZname4e, sZZname4mu,'ggZZ_2e2mu_MCFM67', 'ggZZ_4e_MCFM67', 'ggZZ_4mu_MCFM67']
bkg_samples_shorttags = {sZZname2e2mu:'qqZZ', sZZname4e:'qqZZ', sZZname4mu:'qqZZ', 'ggZZ_2e2mu_MCFM67':'ggZZ', 'ggZZ_4e_MCFM67':'ggZZ', 'ggZZ_4mu_MCFM67':'ggZZ', 'ZX4l_CR':'ZJetsCR', 'ZX4l_CR_4e':'ZJetsCR', "ZX4l_CR_4mu":'ZJetsCR', 'ZX4l_CR_2e2mu':'ZJetsCR'}
bkg_samples_fStates = {sZZname2e2mu:'2e2mu', sZZname4e:'4e', sZZname4mu:'4mu', 'ggZZ_2e2mu_MCFM67':'2e2mu', 'ggZZ_4e_MCFM67':'4e', 'ggZZ_4mu_MCFM67':'4mu', 'ZX4l_CR':'AllChans', 'ZX4l_CR_4e':'4e', 'ZX4l_CR_4mu':'4mu', 'ZX4l_CR_2e2mu':'2e2mu'}
logger.info('Loop over each background sample tags: \n\t{}\n'.format(bkg_sample_tags))
for sample_tag in bkg_sample_tags:
tmpSrcDir = opt.SOURCEDIR
if (sample_tag=='ZX4l_CR'):
# tmpSrcDir = '/eos/user/v/vmilosev/Skim_2018_HZZ/WoW'
tmpSrcDir = opt.SOURCEDIR # FIXME: if the paths for ZX CR is different then we need to update this
# FIXME: Try to understand this syntax
fitTypeZ4l = [['none','doRatio'],['doZ4l','doZ4l']][opt.doZ4l][opt.doRatio]
logger.info("observableBins: "+str(observableBins))
if (" vs " not in obsName ):
#cmd = './main_fiducialXSTemplates '+bkg_samples_shorttags[sample_tag]+' "'+tmpSrcDir+'/'+background_samples[year][sample_tag]+'" '+bkg_samples_fStates[sample_tag]+' '+obsName+' "'+opt.OBSBINS+'" "'+opt.OBSBINS+'" 13TeV templatesXS_'+str(year)+' DTreeXS ' + fitTypeZ4l+ ' 0' +' "' +str(obs_ifJES).lower() +'" '
cmd = './main_fiducialXSTemplates '+bkg_samples_shorttags[sample_tag]+' "'+tmpSrcDir+'/'+background_samples[year][sample_tag]+'" '+bkg_samples_fStates[sample_tag]+' '+obsName+' "'+opt.OBSBINS+'" "'+opt.OBSBINS+'" 13TeV templatesXS'+' DTreeXS ' + fitTypeZ4l+ ' 0 ' +str(int(obs_ifJES)) + ' "'+ str(year)+'"'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# FIXME: URGENT: Here previous command copies all the cout info in variable `output`
# then from the string it is going to extract useful information about bin fraction
tmp_fracs = output.split("[Bin fraction: ")
print('[INFO] tmp_fracs: {}'.format(tmp_fracs))
logger.debug('Length of observables bins: {}'.format(observableBins))
logger.debug('Length of observables bins: {}'.format(len(observableBins)))
# for lines_ in tmp_fracs:
if (" does not exists" in tmp_fracs[0]):
"""INFO: Sometime template maker
"""
logger.error("Seems like template maker did not finish as expected. Please check!!!")
sys.exit(0)
for obsBin in range(0,len(observableBins)-1):
fractionBkg[sample_tag+'_'+bkg_samples_fStates[sample_tag]+'_'+obsName+'_recobin'+str(obsBin)] = 0
lambdajesupBkg[sample_tag+'_'+bkg_samples_fStates[sample_tag]+'_'+obsName+'_recobin'+str(obsBin)] = 0
lambdajesdnBkg[sample_tag+'_'+bkg_samples_fStates[sample_tag]+'_'+obsName+'_recobin'+str(obsBin)] = 0
tmpFrac = float(tmp_fracs[obsBin+1].split("][end fraction]")[0])
if not math.isnan(tmpFrac):
fractionBkg[sample_tag+'_'+bkg_samples_fStates[sample_tag]+'_'+obsName+'_recobin'+str(obsBin)] = tmpFrac
#if (('jet' in obsName) and tmpFrac!=0 and not math.isnan(tmpFrac)):
if ((obs_ifJES) and tmpFrac!=0 and not math.isnan(tmpFrac)):
tmpFrac_up =float(tmp_fracs[obsBin+1].split("Bin fraction (JESup): ")[1].split("]")[0])
tmpFrac_dn =float(tmp_fracs[obsBin+1].split("Bin fraction (JESdn): ")[1].split("]")[0])
if not math.isnan(tmpFrac_up):
lambdajesupBkg[sample_tag+'_'+bkg_samples_fStates[sample_tag]+'_'+obsName+'_recobin'+str(obsBin)] = tmpFrac_up/tmpFrac - 1
if not math.isnan(tmpFrac_dn):
lambdajesdnBkg[sample_tag+'_'+bkg_samples_fStates[sample_tag]+'_'+obsName+'_recobin'+str(obsBin)] = tmpFrac_dn/tmpFrac - 1
else:
for nBins_ in range(len(observableBins)):
# First observable boundary: observableBins[nBins_][0]
# 2nd observable boundary: observableBins[nBins_][1]
bin0_ = '|' + observableBins[nBins_][0][0] + '|' + observableBins[nBins_][0][1] + '|' # FIXME: Discuss with Vukasin
bin1_ = '|' + observableBins[nBins_][1][0] + '|' + observableBins[nBins_][1][1] + '|'
ListObsName = (''.join(obsName.split())).split('vs')
logger.info(ListObsName[0]+ ' : ' + str(bin0_)+"\t"+ListObsName[1] + ' : '+str(bin1_))
cmd = './main_fiducialXSTemplates '+bkg_samples_shorttags[sample_tag]+' "'+tmpSrcDir+'/'+background_samples[year][sample_tag]+'" '+bkg_samples_fStates[sample_tag]+' '+ListObsName[0]+' "'+bin0_+'" "'+bin0_ +'" 13TeV templatesXS'+' DTreeXS ' + fitTypeZ4l+ ' 0 ' + str(int(obs_ifJES)) + ' "'+ str(year)+'"' + ' ' + ListObsName[1]+' "'+bin1_+'" "'+bin1_ +'" '+str(int(obs_ifJES2))
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
tmp_fracs = output.split("[Bin fraction: ")
logger.debug('tmp_fracs: {}'.format(tmp_fracs))
logger.debug('Length of observables bins: {}'.format(len(observableBins)))
tag_ = sample_tag+'_'+bkg_samples_fStates[sample_tag]+'_'+obsName.replace(' ','_')+'_recobin'+str(nBins_)
logger.debug('tag_ : {}'.format(tag_))
fractionBkg[tag_] = 0
lambdajesupBkg[tag_] = 0
lambdajesdnBkg[tag_] = 0
logger.debug('tmp_fracs: {}'.format(tmp_fracs))
logger.debug('tmp_fracs[1]: {}'.format(tmp_fracs[1]))
logger.debug('tmp_fracs[1].split("][end fraction]"): {}'.format(tmp_fracs[1].split("][end fraction]")))
# The current for loop depends on how many times "[Bin fraction: " appears, when we run the `main_fiducialXSTemplates`
tmpFrac = float(tmp_fracs[1].split("][end fraction]")[0])
if not math.isnan(tmpFrac):
fractionBkg[tag_] = tmpFrac
else:
logger.error('total entries in the current bin is isnan. Please check...')
#FIXME: should we exit program here or not???
#if (('jet' in ListObsName[0] or 'jet' in ListObsName[1]) and tmpFrac!=0 and not math.isnan(tmpFrac)): # FIXME: this part will have issue when we have double differential obs.
if ((obs_ifJES or obs_ifJES2) and tmpFrac!=0 and not math.isnan(tmpFrac)): # FIXME: this part will have issue when we have double differential obs.
# FIXME: Check jets information, if its passing correctly or not.
# FIXME: Again we are grabbing things from the log of `main_fiducialXSTemplates`. The next result depends on the couts of `main_fiducialXSTemplates`
logger.debug('tmp_fracs[1].split("Bin fraction (JESup): "): {}'.format(tmp_fracs[1].split("Bin fraction (JESup): ")))
logger.debug('tmp_fracs[1].split("Bin fraction (JESup): ")[1]: {}'.format(tmp_fracs[1].split("Bin fraction (JESup): ")[1]))
logger.debug('tmp_fracs[1].split("Bin fraction (JESup): ")[1].split("]"): {}'.format(tmp_fracs[1].split("Bin fraction (JESup): ")[1].split("]")))
tmpFrac_up =float(tmp_fracs[1].split("Bin fraction (JESup): ")[1].split("]")[0])
tmpFrac_dn =float(tmp_fracs[1].split("Bin fraction (JESdn): ")[1].split("]")[0])
if not math.isnan(tmpFrac_up):
lambdajesupBkg[tag_] = tmpFrac_up/tmpFrac - 1
if not math.isnan(tmpFrac_dn):
lambdajesdnBkg[tag_] = tmpFrac_dn/tmpFrac - 1
os.chdir(currentDir)
logger.debug('observableBins: {}'.format(observableBins))
logger.debug('fractionBkg: {}'.format(fractionBkg))
logger.debug('lambdajesupBkg: {}'.format(lambdajesupBkg))
logger.debug('lambdajesdnBkg: {}'.format(lambdajesdnBkg))
OutputFileName = 'inputs_bkg_'+{0:'',1:'z4l_'}[int(opt.doZ4l)]+obsName.replace(' ','_')+'.py'
logger.debug('File Name: '+OutputFileName)
with open(datacardInputs.format(year = year) + '/' + OutputFileName, 'w') as f:
f.write('observableBins = ' +json.dumps(observableBins)+';\n')
f.write('fractionsBackground = '+json.dumps(fractionBkg) +';\n')
f.write('lambdajesupBkg = ' +json.dumps(lambdajesupBkg)+';\n')
f.write('lambdajesdnBkg = ' +json.dumps(lambdajesdnBkg)+';\n')
### Extract the XS-specific uncertainties for given obs and bin, for all final states (differential)
def extractUncertainties(obsName, observableBinDn, observableBinUp):
"""This function should compute the uncertanities
THis is now computed using a separate script named `getUnc_Unc.py`
I am keeping this, thinking that later we will keep only this and remove the
script `RunEverything.py`
# FIXME
"""
print ('[Extracting uncertainties - range ('+observableBinDn+', '+observableBinUp+')]')
cmd = 'some command...with some parameters...'
#processCmd(cmd, get_linenumber(), os.path.basename(__file__))
### Produce datacards for given obs and bin, for all final states
#def produceDatacards(obsName, observableBins, modelName, physicalModel, obs_ifJES, obs_ifJES2, year):
def produceDatacards(obsName, observableBins, modelName, physicalModel, obs_ifJES, obs_ifJES2, obs_ifAbs, obs_ifAbs2, year):
"""Produce workspace/datacards for the given observable and bins
Args:
obsName (str): Name of observable
observableBins (array): Array having bin boundaries
modelName (str): Name of model. For example: SM_125
physicalModel (str): version of model, for example: v2
"""
logger.info ('Producing workspace/datacards for obsName - '+obsName+', bins - '+str(observableBins)+']')
logger.debug('obsName: {}'.format(obsName))
logger.debug('observableBins: {}'.format(observableBins))
ListObsName = (''.join(obsName.split())).split('vs')
logger.debug('ListObsName: {}'.format(ListObsName))
fStates = channels_central
# INFO: in case of 2D obs nbins is n else its n-1
logger.debug("len(observableBins): = "+str(len(observableBins)))
nBins = len(observableBins) -1
if len(ListObsName) == 2: # INFO: for 2D this list size == 2
nBins = len(observableBins)
for fState in fStates:
logger.info('''VM Creating datacards for:
\tobsName: {obsName}
\tfState: {fState}
\tnBins: {nBins}
\tobservableBins: {observableBins}
\tmodelName: {modelName}
\tphysicalModel: {physicalModel}'''.format(
obsName = ListObsName , fState = fState , nBins = nBins ,
observableBins = observableBins , modelName = modelName ,
physicalModel = physicalModel
))
if (obs_ifJES):
# Call read_jes module
command_read_jes = 'python python/read_jes.py -c {channel} -b {nbins} -v {variable} -y {year}'.format(channel = fState, nbins = nBins, variable = obsName.replace(' ','_'), year = year)
processCmd(command_read_jes, get_linenumber(), os.path.basename(__file__))
if (not (obsName=="mass4l")):
logger.debug("Running the datacard_maker.py...")
processCmd("python python/datacard_maker.py -c {} -b {} -y {} -m {}".format(fState, nBins, year,physicalModel),get_linenumber(), os.path.basename(__file__))
logger.debug("Completed the datacard_maker.py...")
for obsBin in range(nBins):
logger.debug("=="*51)
logger.debug("""
\ttmpObsName = {},
\tfState = {},
\tnBins = {}, obsBin = {},
\tobservableBins = {},
\tFalse = {}, True = {},
\tmodelName = {}, physicalModel = {},
\tobs_ifJES = {}, obs_ifJES2 = {}
""".format(ListObsName, fState, nBins, obsBin, observableBins, False, True, modelName, physicalModel, obs_ifJES, obs_ifJES2))
#ndata = createXSworkspace(ListObsName, fState, nBins, obsBin, observableBins, False, True, modelName, physicalModel, year, obs_ifJES, obs_ifJES2)
ndata = createXSworkspace(ListObsName, fState, nBins, obsBin, observableBins, False, True, modelName, physicalModel, year, obs_ifJES, obs_ifJES2, obs_ifAbs, obs_ifAbs2)
CopyDataCardToOutputDir = "cp xs_125.0_"+str(nBins)+"bins/hzz4l_"+fState+"S_13TeV_xs_bin"+str(obsBin)+".txt "+combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName.replace(' ','_')+"_bin"+str(obsBin)+"_"+physicalModel+".txt"
processCmd(CopyDataCardToOutputDir, get_linenumber(), os.path.basename(__file__))
UpdateObservationValue = "sed -i 's~observation [0-9]*~observation "+str(ndata)+"~g' "+combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName.replace(' ','_')+"_bin"+str(obsBin)+"_"+physicalModel+".txt"
processCmd(UpdateObservationValue, get_linenumber(), os.path.basename(__file__))
UpdateDatabin = "sed -i 's~_xs.Databin"+str(obsBin)+"~_xs_"+modelName+"_"+obsName.replace(' ','_')+"_"+physicalModel+".Databin"+str(obsBin)+"~g' "+combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName.replace(' ','_')+"_bin"+str(obsBin)+"_"+physicalModel+".txt"
os.system(UpdateDatabin)
if (obs_ifJES):
# INFO: This "command_append" here just to ensures appending new line from next one. Below, "command_append" is just added a blank in card
command_append = "echo '" "' >> " +combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName.replace(' ','_')+"_bin"+str(obsBin)+"_"+physicalModel+".txt"
os.system(command_append)
with open("Inputs/JES/{year}/{obsName}/datacardLines_JESnuis_{obsName}_{channel}.txt".format(year = year, obsName = obsName.replace(' ','_'), channel= fState)) as f:
data = f.read()
d = ast.literal_eval(data)
for keys, values in d.items():
if (keys == obsName.replace(' ','_') + '_' + str(obsBin)):
for value in values:
command_append = "echo '"+str(value)+" ' >> " +combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName.replace(' ','_')+"_bin"+str(obsBin)+"_"+physicalModel+".txt"
os.system(command_append)
else:
logger.debug("Running the datacard_maker.py...")
dataCardMakerCommand = "python python/datacard_maker.py -c {} -b {} -y {}".format(fState, 1, year)
processCmd(dataCardMakerCommand, get_linenumber(), os.path.basename(__file__))
logger.debug("Completed the datacard_maker.py...")
#ndata = createXSworkspace(ListObsName,fState, nBins, 0, observableBins, False, True, modelName, physicalModel, year, obs_ifJES, obs_ifJES2)
print ("TEST")
print (ListObsName)
ndata = createXSworkspace(ListObsName,fState, nBins, 0, observableBins, False, True, modelName, physicalModel, year, obs_ifJES, obs_ifJES2, obs_ifAbs, obs_ifAbs2)
if obsName=='mass4l':
CopyDataCardToOutputDir = "cp xs_125.0_1bin/hzz4l_"+fState+"S_13TeV_xs_inclusive_bin0.txt "+combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName+"_bin0_"+physicalModel+".txt"
processCmd(CopyDataCardToOutputDir, get_linenumber(), os.path.basename(__file__))
if obsName=='mass4lREFIT':
os.system("cp xs_125.0_1bin/hzz4l_"+fState+"S_13TeV_xs_inclusiveREFIT_bin0.txt "+combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName+"_bin0_"+physicalModel+".txt")
UpdateObservationValue = "sed -i 's~observation [0-9]*~observation "+str(ndata)+"~g' "+combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName+"_bin0_"+physicalModel+".txt"
processCmd(UpdateObservationValue, get_linenumber(), os.path.basename(__file__))
UpdateDatabin = "sed -i 's~_xs.Databin0~_xs_"+modelName+"_"+obsName+"_"+physicalModel+".Databin0~g' "+combineOutputs.format(year = year)+"/hzz4l_"+fState+"S_13TeV_xs_"+obsName+"_bin0_"+physicalModel+".txt"
processCmd(UpdateDatabin, get_linenumber(), os.path.basename(__file__))
def createAsimov_161718(obsName, observableBins, modelName, resultsXS, physicalModel, years = ['2016', '2017', '2018'], obs_ifJES = False):
logger.info('[Combined datacard for all 3 years for obsName "'+obsName.replace(' ','_')+'" using '+modelName+']')
logger.debug('obsName: {}'.format(obsName))
logger.debug('observableBins: {}'.format(observableBins))
logger.debug("resultsXS: {}".format(resultsXS))
ListObsName = (''.join(obsName.split())).split('vs')
logger.debug('ListObsName: {}'.format(ListObsName))
# INFO: in case of 2D obs nbins is n else its n-1
nBins = len(observableBins) -1
if len(ListObsName) == 2: # INFO: for 2D this list size == 2
nBins = len(observableBins)
logger.debug("nBins: = "+str(nBins))
GetDirectory(combineOutputs.format(year = 'allYear'))
fStates = channels_central
AllChannelCombinedDataCards = 'hzz4l_all_13TeV_xs_{obsName}_bin_{physicalModel}.txt '
cmd = 'combineCards.py '
for year in years:
pathOfCard = combineOutputs.format(year = year)
logger.info("path of datacard for year - {}: {}".format(year, pathOfCard))
# cmd += pathOfCard + '/' + AllChannelCombinedDataCards.format(obsName = obsName.replace(' ','_'), physicalModel = physicalModel)
cmd += obsName+'_4e4mu2e2mu_'+year + '=' + pathOfCard + '/' + AllChannelCombinedDataCards.format(obsName = obsName.replace(' ','_'), physicalModel = physicalModel)
# FIXME: Hardcoded dir string year
cmd += ' > ' + combineOutputs.format(year = 'allYear') + '/' + AllChannelCombinedDataCards.format(obsName = obsName.replace(' ','_'), physicalModel = physicalModel)
processCmd(cmd, get_linenumber(), os.path.basename(__file__))
if (not opt.LUMISCALE=="1.0"):
os.system('echo " " >> '+ combineOutputs.format(year = 'allYear') + '/' +'hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.txt')
os.system('echo "lumiscale rateParam * * '+opt.LUMISCALE+'" >> '+combineOutputs.format(year = 'allYear') + '/' +'hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.txt')
os.system('echo "nuisance edit freeze lumiscale" >> '+combineOutputs.format(year = 'allYear') + '/' +'hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.txt')
# text-to-workspace
if (physicalModel=="v2"):
cmd = 'text2workspace.py ' + combineOutputs.format(year = 'allYear') + '/' + AllChannelCombinedDataCards.format(obsName = obsName.replace(' ','_'), physicalModel = physicalModel) + ' -P HiggsAnalysis.CombinedLimit.HZZ4L_Fiducial_v2:differentialFiducialV2 --PO higgsMassRange=115,135 --PO nBin='+str(nBins)+' -o ' + combineOutputs.format(year = 'allYear') + '/' +(AllChannelCombinedDataCards.format(obsName = obsName.replace(' ','_'), physicalModel = physicalModel)).replace('.txt', '.root')
processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# FIXME: Can we improve this?
cmd = 'mv ' + combineOutputs.format(year = 'allYear') + '/' +(AllChannelCombinedDataCards.format(obsName = obsName.replace(' ','_'), physicalModel = physicalModel)).replace('.txt', '.root') + ' ' + combineOutputs.format(year = 'allYear') + '/' +modelName+'_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.root'
processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# Run the Combine
logger.info("Going to run the combine commands...")
cmd = ''
if (physicalModel=="v2"):
cmd = 'combine -n '+obsName.replace(' ','_')+'_'+str("allYear")+unblindString+' -M MultiDimFit '+combineOutputs.format(year = "allYear")+'/'+modelName+'_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.root -m '+opt.ASIMOVMASS+' --setParameters '
logger.debug("command:\n\t{}".format(cmd))
for fState in fStates:
for obsBin in range(nBins):
fidxs = 0
for year in years:
# import acc factors
logger.debug("Import module: {}".format((datacardInputs.format(year = year)).replace('/','.') + '.' + 'inputs_sig_'+obsName.replace(' ','_')))
_temp = __import__((datacardInputs.format(year = year)).replace('/','.') + '.' + 'inputs_sig_'+obsName.replace(' ','_'), globals(), locals(), ['acc'], -1)
acc = _temp.acc
fidxs += higgs_xs['ggH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['ggH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
#fidxs += acc['ggH_HRes_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['VBF_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['VBF_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['WH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['WH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['ZH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['ZH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['ttH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['ttH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
cmd = cmd + 'r'+fState+'Bin'+str(obsBin)+'='+str(fidxs/3.0)+',' # FIXME: Check division by 3.0
cmd = cmd+ 'MH='+opt.ASIMOVMASS
logger.debug("command+:\n\t{}".format(cmd))
for fState in fStates:
for obsBin in range(nBins):
cmd = cmd + ' -P r'+fState+'Bin'+str(obsBin)
if (opt.FIXMASS=="False"):
cmd = cmd + ' -P MH '
else:
cmd = cmd + ' --floatOtherPOIs=0'
cmd = cmd +' -t -1 --saveWorkspace --saveToys --saveFitResult '
logger.debug("command+:\n\t{}".format(cmd))
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# INFO: Move combine output to combineOutputs.format(year = opt.ERA) directory
processCmd('mv higgsCombine'+obsName.replace(' ','_')+'_'+str("allYear")+unblindString+'.MultiDimFit.mH'+opt.ASIMOVMASS.rstrip('.0')+'.123456.root '+combineOutputs.format(year = "allYear")+'/'+modelName+'_all_'+obsName.replace(' ','_')+'_13TeV_Asimov_'+physicalModel+unblindString+'.root', get_linenumber(), os.path.basename(__file__), 1)
#cmd = cmd.replace(' --freezeNuisanceGroups r4muBin0,r4eBin0,r2e2muBin0','')
#cmd = cmd.replace(' --X-rtd TMCSO_PseudoAsimov=1000000','')
cmd = cmd + ' --algo=singles --cl=0.68 --robustFit=1'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
processCmd('mv higgsCombine'+obsName.replace(' ','_')+'_'+str("allYear")+unblindString+'.MultiDimFit.mH'+opt.ASIMOVMASS.rstrip('.0')+'.123456.root '+combineOutputs.format(year = "allYear")+'/', get_linenumber(), os.path.basename(__file__))
# parse the results for all the bins and the given final state
tmp_resultsXS = {}
for fState in fStates:
for obsBin in range(nBins):
binTag = str(obsBin)
tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag] = parseXSResults(output,'r'+fState+'Bin'+str(obsBin)+' :')
# merge the results for 3 final states, for the given bins
for obsBin in range(nBins):
binTag = str(obsBin)
resultsXS['AsimovData_'+obsName+'_genbin'+binTag] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
for fState in fStates:
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
tmp_central = 0.0
tmp_uncerDn = 0.0
tmp_uncerUp = 0.0
for fState in fStates:
tmp_central += tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag]['central']
tmp_uncerDn += tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag]['uncerDn']**2
tmp_uncerUp += tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag]['uncerUp']**2
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag]['central'] = float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag]['central']))
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag]['uncerDn'] = -float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag]['uncerDn']))
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag]['uncerUp'] = +float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag]['uncerUp']))
resultsXS['AsimovData_'+obsName+'_genbin'+binTag]['central'] = float("{0:.5f}".format(tmp_central))
resultsXS['AsimovData_'+obsName+'_genbin'+binTag]['uncerDn'] = -float("{0:.5f}".format(tmp_uncerDn**0.5))
resultsXS['AsimovData_'+obsName+'_genbin'+binTag]['uncerUp'] = +float("{0:.5f}".format(tmp_uncerUp**0.5))
# run combine with no systematics for stat only uncertainty
if (opt.SYS):
cmd = cmd + ' --freezeParameters allConstrainedNuisances'
# cmd = cmd + ' --freezeParameters '
# cmd = cmd + 'CMS_fakeH_p1_1'+year+',CMS_fakeH_p1_2'+year+',CMS_fakeH_p1_3'+year+',CMS_fakeH_p3_1'+year+',CMS_fakeH_p3_2'+year+',CMS_fakeH_p3_3'+year+','
# cmd += 'CMS_eff_e,CMS_eff_m,CMS_hzz2e2mu_Zjets_'+str(year)+',CMS_hzz4e_Zjets_'+str(year)+',CMS_hzz4mu_Zjets_'+str(year)+',QCDscale_VV,QCDscale_ggVV,kfactor_ggzz,lumi_13TeV_'+str(year)+'_uncorrelated,norm_fakeH,pdf_gg,pdf_qqbar,CMS_zz4l_sigma_e_sig_'+str(year)+',CMS_zz4l_sigma_m_sig_'+str(year)+',CMS_zz4l_n_sig_1_'+str(year)+',CMS_zz4l_n_sig_2_'+str(year)+',CMS_zz4l_n_sig_3_'+str(year)+',CMS_zz4l_mean_e_sig_'+str(year)+',CMS_zz4l_mean_m_sig_'+str(year)+',CMS_zz4l_sigma_e_sig_'+str(year)+','
# if (obs_ifJES):
# jes_sources_year_updated = [x.replace('year', year) for x in jes_sources]
# for jes_source in jes_sources:
# cmd += ",CMS_scale_j_"+ jes_sources_year_updated[jes_sources.index(jes_source)]
# cmd = cmd + ',lumi_13TeV_correlated_16_17_18,lumi_13TeV_correlated_17_18'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
#for year in years:
# cmd = cmd + 'CMS_fakeH_p1_1'+year+',CMS_fakeH_p1_2'+year+',CMS_fakeH_p1_3'+year+',CMS_fakeH_p3_1'+year+',CMS_fakeH_p3_2'+year+',CMS_fakeH_p3_3'+year+','
# cmd += 'CMS_eff_e,CMS_eff_m,CMS_hzz2e2mu_Zjets_'+str(year)+',CMS_hzz4e_Zjets_'+str(year)+',CMS_hzz4mu_Zjets_'+str(year)+',QCDscale_VV,QCDscale_ggVV,kfactor_ggzz,lumi_13TeV_'+str(year)+'_uncorrelated,norm_fakeH,pdf_gg,pdf_qqbar,CMS_zz4l_sigma_e_sig_'+str(year)+',CMS_zz4l_sigma_m_sig_'+str(year)+',CMS_zz4l_n_sig_1_'+str(year)+',CMS_zz4l_n_sig_2_'+str(year)+',CMS_zz4l_n_sig_3_'+str(year)+',CMS_zz4l_mean_e_sig_'+str(year)+',CMS_zz4l_mean_m_sig_'+str(year)+',CMS_zz4l_sigma_e_sig_'+str(year)+','
# if len(years)>1:
# cmd = cmd + 'lumi_13TeV_correlated_16_17_18,lumi_13TeV_correlated_17_18'
#if (obs_ifJES):
# jes_sources_year_updated = [x.replace('year', year) for x in jes_sources]
# for jes_source in jes_sources:
# cmd += ",CMS_scale_j_"+ jes_sources_year_updated[jes_sources.index(jes_source)]
#logger.error(cmd)
#sys.exit()
# processCmd(cmd, get_linenumber(), os.path.basename(__file__))
processCmd('mv higgsCombine'+obsName.replace(' ','_')+'_'+str("allYear")+unblindString+'.MultiDimFit.mH'+opt.ASIMOVMASS.rstrip('.0')+'.123456.root '+combineOutputs.format(year = "allYear")+'/', get_linenumber(), os.path.basename(__file__))
for fState in fStates:
for obsBin in range(nBins):
binTag = str(obsBin)
tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag+'_statOnly'] = parseXSResults(output,'r'+fState+'Bin'+str(obsBin)+' :')
# merge the results for 3 final states, for the given bins
for obsBin in range(nBins):
binTag = str(obsBin)
resultsXS['AsimovData_'+obsName+'_genbin'+binTag+'_statOnly'] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
for fState in fStates:
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag+'_statOnly'] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
tmp_central = 0.0
tmp_uncerDn = 0.0
tmp_uncerUp = 0.0
for fState in fStates:
tmp_central += tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag+'_statOnly']['central']
tmp_uncerDn += tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag+'_statOnly']['uncerDn']**2
tmp_uncerUp += tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag+'_statOnly']['uncerUp']**2
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag+'_statOnly']['central'] = float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag+'_statOnly']['central']))
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag+'_statOnly']['uncerDn'] = -float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag+'_statOnly']['uncerDn']))
resultsXS['AsimovData_'+obsName+'_'+fState+'_genbin'+binTag+'_statOnly']['uncerUp'] = +float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName+'_genbin'+binTag+'_statOnly']['uncerUp']))
resultsXS['AsimovData_'+obsName+'_genbin'+binTag+'_statOnly']['central'] = float("{0:.5f}".format(tmp_central))
resultsXS['AsimovData_'+obsName+'_genbin'+binTag+'_statOnly']['uncerDn'] = -float("{0:.5f}".format(tmp_uncerDn**0.5))
resultsXS['AsimovData_'+obsName+'_genbin'+binTag+'_statOnly']['uncerUp'] = +float("{0:.5f}".format(tmp_uncerUp**0.5))
return resultsXS
### Create the asimov dataset and return fit results
#def createAsimov(obsName, observableBins, modelName, resultsXS, physicalModel, year = 2018, obs_ifJES = False):
def createAsimov(obsName, observableBins, modelName, resultsXS, physicalModel, year = 2018, obs_ifJES = False, obs_ifAbs = False):
"""Create the Asimov dataset and return the fit results.
Args:
obsName (str): Name of observable
observableBins (array): Array having bin boundaries
modelName (str): Name of model
resultsXS (_type_): _description_
physicalModel (_type_): _description_
"""
logger.info('[Producing/merging workspaces and datacards for obsName "'+obsName.replace(' ','_')+'" using '+modelName+']')
logger.debug('obsName: {}'.format(obsName))
logger.debug('observableBins: {}'.format(observableBins))
ListObsName = (''.join(obsName.split())).split('vs')
logger.debug('ListObsName: {}'.format(ListObsName))
# INFO: in case of 2D obs nbins is n else its n-1
nBins = len(observableBins) -1
if len(ListObsName) == 2: # INFO: for 2D this list size == 2
nBins = len(observableBins)
logger.debug("nBins: = "+str(nBins))
print("yuji debug before Run combineCards in asimov")
# Run combineCards and text2workspace
currentDir = os.getcwd(); os.chdir('./'+combineOutputs.format(year = year)+'/') # FIXME: Hardcoded directory
fStates = channels_central
for fState in fStates:
"""Combine cards for all bins corresponding to each final state.
"""
# if (nBins>1):
cmd = 'combineCards.py '
for obsBin in range(nBins):
cmd += obsName + '_' + fState + 'S_'+str(obsBin)+'_'+year + '=' +'hzz4l_'+fState+'S_13TeV_xs_'+obsName.replace(' ','_')+'_bin'+str(obsBin)+'_'+physicalModel+'.txt '
cmd += ' > hzz4l_'+fState+'S_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt' # Output combine card name
processCmd(cmd, get_linenumber(), os.path.basename(__file__), 1)
print("yuji debug station A")
# combine 3 final state cards
cmd = 'combineCards.py '
for channels_ in fStates:
cmd += obsName + '_' + channels_ + 'S_' + year + '=hzz4l_' + channels_ + 'S_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt '
cmd += ' > hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt'
processCmd(cmd, get_linenumber(), os.path.basename(__file__), 1)
if (not opt.LUMISCALE=="1.0"):
os.system('echo " " >> hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt')
os.system('echo "lumiscale rateParam * * '+opt.LUMISCALE+'" >> hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt')
os.system('echo "nuisance edit freeze lumiscale" >> hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt')
print("yuji debug station B")
# text-to-workspace
if (physicalModel=="v2"):
cmd = 'text2workspace.py hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt -P HiggsAnalysis.CombinedLimit.HZZ4L_Fiducial_v2:differentialFiducialV2 --PO higgsMassRange=115,135 --PO nBin='+str(nBins)+' -o hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.root'
processCmd(cmd, get_linenumber(), os.path.basename(__file__))
print("debug in v2")
if (physicalModel=="v3"):
cmd = 'text2workspace.py hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.txt -P HiggsAnalysis.CombinedLimit.HZZ4L_Fiducial:differentialFiducialV3 --PO higgsMassRange=115,135 --PO nBin='+str(nBins)+' -o hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.root'
processCmd(cmd, get_linenumber(), os.path.basename(__file__))
print("debug in v3")
print("yuji debug before mv hzz4l")
# FIXME: Can we improve this?
cmd = 'mv hzz4l_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.root '+modelName+'_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.root'
processCmd(cmd, get_linenumber(), os.path.basename(__file__))
os.chdir(currentDir)
# import acc factors
_temp = __import__('inputs_sig_'+obsName.replace(' ','_'), globals(), locals(), ['acc'], -1)
acc = _temp.acc
# Run the Combine
logger.info("Going to run the combine commands...")
if (physicalModel=="v2"):
cmd = 'combine -n '+obsName.replace(' ','_')+'_'+str(year)+unblindString+' -M MultiDimFit '+combineOutputs.format(year = year)+'/'+modelName+'_all_13TeV_xs_'+obsName.replace(' ','_')+'_bin_'+physicalModel+'.root -m '+opt.ASIMOVMASS+' --setParameters '
for fState in fStates:
for obsBin in range(nBins):
fidxs = 0
fidxs += higgs_xs['ggH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['ggH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
#fidxs += acc['ggH_HRes_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['VBF_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['VBF_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['WH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['WH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['ZH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['ZH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
fidxs += higgs_xs['ttH_'+opt.ASIMOVMASS]*higgs4l_br[opt.ASIMOVMASS+'_'+fState]*acc['ttH_powheg_JHUgen_125.38_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]
cmd = cmd + 'r'+fState+'Bin'+str(obsBin)+'='+str(fidxs)+','
cmd = cmd+ 'MH='+opt.ASIMOVMASS
for fState in fStates:
for obsBin in range(nBins):
cmd = cmd + ' -P r'+fState+'Bin'+str(obsBin)
if (opt.FIXMASS=="False"):
cmd = cmd + ' -P MH '
else:
cmd = cmd + ' --floatOtherPOIs=0'
cmd = cmd +' -t -1 --saveWorkspace --saveToys'
logger.debug("command+:\n\t{}".format(cmd))
GetDirectory(combineOutputs.format(year = year))
#cmd += ' --freezeNuisanceGroups r4muBin0,r4eBin0,r2e2muBin0'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# INFO: Move combine output to combineOutputs.format(year = opt.ERA) directory
processCmd('mv higgsCombine'+obsName.replace(' ','_')+'_'+str(year)+unblindString+'.MultiDimFit.mH'+opt.ASIMOVMASS.rstrip('.0')+'.123456.root '+combineOutputs.format(year = year)+'/'+modelName+'_all_'+obsName.replace(' ','_')+'_13TeV_Asimov_'+physicalModel+unblindString+'.root', get_linenumber(), os.path.basename(__file__), 1)
#cmd = cmd.replace(' --freezeNuisanceGroups r4muBin0,r4eBin0,r2e2muBin0','')
#cmd = cmd.replace(' --X-rtd TMCSO_PseudoAsimov=1000000','')
cmd = cmd + ' --algo=singles --cl=0.68 --robustFit=1'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# lsCMD = "ls *.root"
# processCmd(lsCMD, get_linenumber(), os.path.basename(__file__))
# INFO: Move combine output to combineOutputs.format(year = opt.ERA) directory
processCmd('mv higgsCombine'+obsName.replace(' ','_')+'_'+str(year)+unblindString+'.MultiDimFit.mH'+opt.ASIMOVMASS.rstrip('.0')+'.123456.root '+combineOutputs.format(year = year)+'/', get_linenumber(), os.path.basename(__file__))
# parse the results for all the bins and the given final state
tmp_resultsXS = {}
for fState in fStates:
for obsBin in range(nBins):
binTag = str(obsBin)
tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag] = parseXSResults(output,'r'+fState+'Bin'+str(obsBin)+' :')
# merge the results for 3 final states, for the given bins
for obsBin in range(nBins):
binTag = str(obsBin)
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
for fState in fStates:
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
tmp_central = 0.0
tmp_uncerDn = 0.0
tmp_uncerUp = 0.0
for fState in fStates:
tmp_central += tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag]['central']
tmp_uncerDn += tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag]['uncerDn']**2
tmp_uncerUp += tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag]['uncerUp']**2
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag]['central'] = float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag]['central']))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag]['uncerDn'] = -float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag]['uncerDn']))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag]['uncerUp'] = +float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag]['uncerUp']))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag]['central'] = float("{0:.5f}".format(tmp_central))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag]['uncerDn'] = -float("{0:.5f}".format(tmp_uncerDn**0.5))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag]['uncerUp'] = +float("{0:.5f}".format(tmp_uncerUp**0.5))
# run combine with no systematics for stat only uncertainty
if (opt.SYS):
# Test VM: cmd = cmd + ' --freezeNuisanceGroups CMS_fakeH_p1_1_8,CMS_fakeH_p1_2_8,CMS_fakeH_p1_3_8,CMS_fakeH_p3_1_8,CMS_fakeH_p3_2_8,CMS_fakeH_p3_3_8 --freezeParameters allConstrainedNuisances'
cmd = cmd + ' --freezeParameters allConstrainedNuisances'
# cmd = cmd + ' --freezeParameters '
# cmd = cmd + 'CMS_fakeH_p1_1'+year+',CMS_fakeH_p1_2'+year+',CMS_fakeH_p1_3'+year+',CMS_fakeH_p3_1'+year+',CMS_fakeH_p3_2'+year+',CMS_fakeH_p3_3'+year+','
# cmd += 'CMS_eff_e,CMS_eff_m,CMS_hzz2e2mu_Zjets_'+str(year)+',CMS_hzz4e_Zjets_'+str(year)+',CMS_hzz4mu_Zjets_'+str(year)+',QCDscale_VV,QCDscale_ggVV,kfactor_ggzz,lumi_13TeV_'+str(year)+'_uncorrelated,norm_fakeH,pdf_gg,pdf_qqbar,CMS_zz4l_sigma_e_sig_'+str(year)+',CMS_zz4l_sigma_m_sig_'+str(year)+',CMS_zz4l_n_sig_1_'+str(year)+',CMS_zz4l_n_sig_2_'+str(year)+',CMS_zz4l_n_sig_3_'+str(year)+',CMS_zz4l_mean_e_sig_'+str(year)+',CMS_zz4l_mean_m_sig_'+str(year)+',CMS_zz4l_sigma_e_sig_'+str(year)+','
# if (obs_ifJES):
# jes_sources_year_updated = [x.replace('year', year) for x in jes_sources]
# for jes_source in jes_sources:
# cmd += ",CMS_scale_j_"+ jes_sources_year_updated[jes_sources.index(jes_source)]
# cmd = cmd + 'lumi_13TeV_correlated_16_17_18,lumi_13TeV_correlated_17_18'
output = processCmd(cmd, get_linenumber(), os.path.basename(__file__))
# FIXME: The current combine command and the previous
# combine command generates output with exactly same name (obvious).
# But, we are not renaming them so its updated with current command.
# check if this is fine?
# print("="*51)
# lsCMD = 'echo "ls *.root";ls *.root'
# processCmd(lsCMD, get_linenumber(), os.path.basename(__file__))
processCmd('mv higgsCombine'+obsName.replace(' ','_')+'_'+str(year)+unblindString+'.MultiDimFit.mH'+opt.ASIMOVMASS.rstrip('.0')+'.123456.root '+combineOutputs.format(year = year)+'/', get_linenumber(), os.path.basename(__file__))
for fState in fStates:
for obsBin in range(nBins):
binTag = str(obsBin)
tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly'] = parseXSResults(output,'r'+fState+'Bin'+str(obsBin)+' :')
# merge the results for 3 final states, for the given bins
for obsBin in range(nBins):
binTag = str(obsBin)
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly'] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
for fState in fStates:
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag+'_statOnly'] = {'central':0.0, 'uncerDn':0.0, 'uncerUp':0.0}
tmp_central = 0.0
tmp_uncerDn = 0.0
tmp_uncerUp = 0.0
for fState in fStates:
tmp_central += tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['central']
tmp_uncerDn += tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['uncerDn']**2
tmp_uncerUp += tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['uncerUp']**2
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag+'_statOnly']['central'] = float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['central']))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag+'_statOnly']['uncerDn'] = -float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['uncerDn']))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_'+fState+'_genbin'+binTag+'_statOnly']['uncerUp'] = +float("{0:.5f}".format(tmp_resultsXS[modelName+'_'+fState+'_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['uncerUp']))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['central'] = float("{0:.5f}".format(tmp_central))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['uncerDn'] = -float("{0:.5f}".format(tmp_uncerDn**0.5))
resultsXS['AsimovData_'+obsName.replace(' ','_')+'_genbin'+binTag+'_statOnly']['uncerUp'] = +float("{0:.5f}".format(tmp_uncerUp**0.5))
return resultsXS
# parse the fit results from the MultiDim fit output "resultLog", for the bin and final state designated by "rTag"
def parseXSResults(resultLog, rTag):
"""parse the fit results from the MultiDim fit output "resultLog", for the bin and final state designated by "rTag"
Args:
resultLog (str): This contains full log output of the combine log from MultiDim filt
rTag (str): This contains the bin and final state inforamation, that will help to get useful information
Returns:
dict: dict having central, up and down values.
"""
try:
fXS_c = float(resultLog.split(rTag)[1].split(' (68%)')[0].strip().split(" ")[0])
fXS_d = float('-'+resultLog.split(rTag)[1].split(' (68%)')[0].strip().split(" -")[1].split("/+")[0])
fXS_u = float('+'+resultLog.split(rTag)[1].split(' (68%)')[0].strip().split(" -")[1].split("/+")[1])
fXS = {'central':fXS_c, 'uncerDn':fXS_d, 'uncerUp':fXS_u}
return fXS
except IndexError:
logger.error("Parsing Failed!!! Inserting dummy values!!! check log!!!")
fXS = {'central':-1.0, 'uncerDn':0.0, 'uncerUp':0.0}
return fXS
### Extract the results and do plotting
def extractResults(obsName, observableBins, modelName, physicalModel, asimovModelName, asimovPhysicalModel, resultsXS, year = 2018):
# Run combineCards and text2workspace
logger.info('[Extract the results and do plotting obsName "'+obsName.replace(' ','_')+'" using '+modelName+']')
logger.debug("physicalModel: {}".format(physicalModel))
logger.debug("resultsXS: {}".format(resultsXS))
logger.debug('obsName: {}'.format(obsName))
logger.debug('observableBins: {}'.format(observableBins))
ListObsName = (''.join(obsName.split())).split('vs')
logger.debug('ListObsName: {}'.format(ListObsName))
# INFO: in case of 2D obs nbins is n else its n-1
logger.debug("len(observableBins): = "+str(len(observableBins)))
nBins = len(observableBins) -1
if len(ListObsName) == 2: # INFO: for 2D this list size == 2
nBins = len(observableBins)
logger.debug("nBins: = "+str(nBins))
currentDir = os.getcwd(); os.chdir(combineOutputs.format(year = year))
fStates = channels_central
AllChannelCombinedDataCards = 'hzz4l_all_13TeV_xs_{obsName}_bin_{physicalModel}.txt '
for fState in fStates:
"""Combine cards for all bins corresponding to each final state.
"""
if year == "allYear":
cmd = 'combineCards.py '
for year_ in ['2016', '2017', '2018']:
pathOfCard = '../'+combineOutputs.format(year = year_)
# pathOfCard = "" # since each card already had the path attached with the root file so no need to attach the path
logger.info("path of datacard for year_ - {}: {}".format(year_, pathOfCard))