diff --git a/Fitter/scripts/EFTFitter.py b/Fitter/scripts/EFTFitter.py index 10b54ce..4805c0d 100644 --- a/Fitter/scripts/EFTFitter.py +++ b/Fitter/scripts/EFTFitter.py @@ -25,9 +25,10 @@ def __init__(self): # WCs lists for easy use # Full list of opeators - self.wcs = ['ctW','ctZ','ctp','cpQM','ctG','cbW','cpQ3','cptb','cpt','cQl3i','cQlMi','cQei','ctli','ctei','ctlSi','ctlTi', 'cQq13', 'cQq83', 'cQq11', 'ctq1', 'cQq81', 'ctq8', 'ctt1', 'cQQ1', 'cQt8', 'cQt1', ] #TOP-22-006 - #self.wcs = ['ctp', 'cpQM', 'cpQ3', 'cpt', 'cptb', 'ctZ', 'ctW', 'cbW'] #TOP-24-004 - #self.wcs = ['ctW','ctZ','ctp','cpQM','ctG','cbW','cpQ3','cptb','cpt','cQl3i','cQlMi','cQei','ctli','ctei','ctlSi','ctlTi'] #TOP-19-001 + # TODO add new WCs + self.wcs = ['cHt', 'ctHRe', 'cHtbRe', 'cQl1i', 'cQl3i', 'cQj18', 'cQj11', 'ctj8', 'cleQt3Rei', 'ctj1', 'ctli', 'cQj31', 'cbW', 'cHQ1', 'cHQ3', 'ctei', 'cQei', 'cleQt1Rei', 'cQj38', 'ctBRe', 'ctWRe', 'ctGRe', 'ctt', 'cQt1', 'cQt8', 'cQQ1'] + # TOP-22-006 + #self.wcs = ['ctW', 'ctZ', 'ctp', 'cpQM', 'ctG', 'cbW', 'cpQ3', 'cptb', 'cpt', 'cQl3i', 'cQlMi', 'cQei', 'ctli', 'ctei', 'ctlSi', 'ctlTi', 'cQq13', 'cQq83', 'cQq11', 'ctq1', 'cQq81', 'ctq8', 'ctt1', 'cQQ1', 'cQt8', 'cQt1'] # Default pair of wcs for 2D scans # Scan ranges of the wcs self.at23v01_2sig_prof = {'ctlTi': [-0.37460753448249934, 0.37456553009578053], 'ctq1': [-0.2184374526296913, 0.21149076882524775], 'ctq8': [-0.6849733172669766, 0.25913443066575303], 'cQq83': [-0.1729820133118581, 0.1643935607326637], 'cQQ1': [-3.0622943484810774, 3.3295794607259963], 'cQt1': [-2.7549317478099273, 2.709140633280104], 'cQt8': [-5.208519171130861, 5.846070290900961], 'ctli': [-1.7912919548903847, 2.1313275751809675], 'cQq81': [-0.6931605230972014, 0.22472203594018714], 'cQlMi': [-1.5701580939064013, 2.308562271259127], 'cbW': [-0.774498659891837, 0.7785760997417999], 'cpQ3': [-0.7999267419016539, 2.1053510234908166], 'ctei': [-1.789099599899366, 2.219641594315739], 'ctlSi': [-2.6297986071455313, 2.6317949240743284], 'ctW': [-0.5595801509610011, 0.46706568379424357], 'cpQM': [-6.143737116235924, 8.096702605664662], 'cQei': [-1.9244486948695827, 1.9614340380338822], 'ctZ': [-0.7286708753806125, 0.6450395866915755], 'cQl3i': [-2.920952087582164, 2.639050916924742], 'ctG': [-0.2780854299079727, 0.23945699061451634], 'cQq13': [-0.07660660914160784, 0.07102679687100005], 'cQq11': [-0.19470826306682623, 0.1947277372883481], 'cptb': [-3.362661720898583, 3.361818811283507], 'ctt1': [-1.5812961023069028, 1.6227282264447938], 'ctp': [-9.34032472114984, 2.2960707194931254], 'cpt': [-10.507390689909828, 7.938575701774626]} @@ -38,56 +39,92 @@ def __init__(self): 'cQei' : (-4.0,4.0), 'cQl3i': (-5.5,5.5), 'cQlMi': (-4.0,4.0), + 'cQl1i': (-4.0,4.0), 'cQq11': (-0.7,0.7), + 'cQj11': (-0.7,0.7), 'cQq13': (-0.35,0.35), + 'cQj31': (-0.35,0.35), 'cQq81': (-1.7,1.5), + 'cQj18': (-1.7,1.5), 'cQq83': (-0.6,0.6), + 'cQj38': (-0.6,0.6), 'cQt1' : (-4.0,4.0), 'cQt8' : (-8.0,8.0), 'cbW' : (-3.0,3.0), 'cpQ3' : (-4.0,4.0), + 'cHQ3' : (-4.0,4.0), 'cpQM' : (-10.0,17.0), + 'cHQ1' : (-10.0,17.0), 'cpt' : (-15.0,15.0), + 'cHt' : (-15.0,15.0), 'cptb' : (-9.0,9.0), + 'cHtbRe' : (-9.0,9.0), 'ctG' : (-0.8,0.8), + 'ctGRe' : (-0.8,0.8), 'ctW' : (-1.5,1.5), + 'ctWRe' : (-1.5,1.5), 'ctZ' : (-2.0,2.0), + 'ctBRe' : (-2.0,2.0), 'ctei' : (-4.0,4.0), 'ctlSi': (-5.0,5.0), + 'cleQt1Rei': (-5.0,5.0), 'ctlTi': (-0.9,0.9), + 'cleQt3Rei': (-0.9,0.9), 'ctli' : (-4.0,4.0), 'ctp' : (-15.0,40.0), + 'ctHRe' : (-15.0,40.0), 'ctq1' : (-0.6,0.6), + 'ctj1' : (-0.6,0.6), 'ctq8' : (-1.4,1.4), + 'ctj8' : (-1.4,1.4), 'ctt1' : (-2.6,2.6), + 'ctt' : (-2.6,2.6), } self.wc_ranges_differential = { 'cQQ1' : (-6.0*2,6.0*2), 'cQei' : (-4.0*2,4.0*2), 'cQl3i': (-5.5*2,5.5*2), 'cQlMi': (-4.0*2,4.0*2), + 'cQl1i': (-4.0*2,4.0*2), 'cQq11': (-0.7*2,0.7*2), + 'cQj11': (-0.7*2,0.7*2), 'cQq13': (-0.35*2,0.35*2), + 'cQj31': (-0.35*2,0.35*2), 'cQq81': (-1.7*2,1.5*2), + 'cQj18': (-1.7*2,1.5*2), 'cQq83': (-0.6*2,0.6*2), + 'cQj38': (-0.6*2,0.6*2), 'cQt1' : (-6.0*2,6.0*2), 'cQt8' : (-10.0*2,10.0*2), 'cbW' : (-3.0*2,3.0*2), 'cpQ3' : (-4.0*2,4.0*2), + 'cHQ3' : (-4.0*2,4.0*2), 'cpQM' : (-15.0*2,20.0*2), + 'cHQ1' : (-15.0*2,20.0*2), 'cpt' : (-15.0*2,15.0*2), + 'cHt' : (-15.0*2,15.0*2), 'cptb' : (-9.0*2,9.0*2), + 'cHtb' : (-9.0*2,9.0*2), 'ctG' : (-0.8*2,0.8*2), + 'ctGRe' : (-0.8*2,0.8*2), 'ctW' : (-1.5*2,1.5*2), + 'ctWRe' : (-1.5*2,1.5*2), 'ctZ' : (-2.0*2,2.0*2), + 'ctBRe' : (-2.0*2,2.0*2), 'ctei' : (-4.0*2,4.0*2), 'ctlSi': (-5.0*2,5.0*2), - 'ctlTi': (-0.9*2,0.9*2), + 'cleQt1Rei': (-5.0*2,5.0*2), + 'cleQ': (-0.9*2,0.9*2), 'ctli' : (-4.0*2,4.0*2), 'ctp' : (-15.0*2,40.0*2), + 'ctH' : (-15.0*2,40.0*2), + 'ctHRe' : (-15.0*2,40.0*2), 'ctq1' : (-0.6*2,0.6*2), + 'ctj1' : (-0.6*2,0.6*2), 'ctq8' : (-1.4*2,1.4*2), + 'ctj8' : (-1.4*2,1.4*2), 'ctt1' : (-2.6*2,2.6*2), + 'ctt' : (-2.6*2,2.6*2), } # Limits appropriate for asimov njets fits (for prof, but can be used for frozen too) @@ -96,28 +133,50 @@ def __init__(self): 'cQei' : (-7.0,7.0), 'cQl3i': (-10.0,10.0), 'cQlMi': (-8.0,8.0), + 'cQl1i': (-8.0,8.0), 'cQq11': (-1.5,1.5), + 'cQj11': (-1.5,1.5), 'cQq13': (-0.6,0.6), + 'cQj31': (-0.6,0.6), 'cQq81': (-4.0,3.0), + 'cQj18': (-4.0,3.0), 'cQq83': (-1.2,1.2), + 'cQj38': (-1.2,1.2), 'cQt1' : (-5.0,5.0), 'cQt8' : (-10.0,10.0), 'cbW' : (-5.0,5.0), 'cpQ3' : (-10.0,7.0), + 'cHQ3' : (-10.0,7.0), 'cpQM' : (-11.0,30.0), + 'cpQM' : (-5.0,5.0), + 'cHQ1' : (-11.0,30.0), + 'cHQ1' : (-5.0,5.0), 'cpt' : (-25.0,20.0), + 'cpt' : (-5.0,5.0), + 'cHt' : (-25.0,20.0), + 'cHt' : (-5.0,5.0), 'cptb' : (-17.0,17.0), + 'cHtbRe' : (-17.0,17.0), 'ctG' : (-1.5,1.5), + 'ctGRe' : (-1.5,1.5), 'ctW' : (-4.0,3.0), + 'ctWRe' : (-4.0,3.0), 'ctZ' : (-4.0,4.0), + 'ctBRe' : (-4.0,4.0), 'ctei' : (-8.0,8.0), 'ctlSi': (-8.0,8.0), + 'cleQt1Rei': (-8.0,8.0), 'ctlTi': (-1.4,1.4), + 'cleQt3Rei': (-1.4,1.4), 'ctli' : (-8.0,8.0), 'ctp' : (-11.0,35.0), + 'ctHRe' : (-11.0,35.0), 'ctq1' : (-1.4,1.4), + 'ctj1' : (-1.4,1.4), 'ctq8' : (-3.0,3.0), + 'ctj8' : (-3.0,3.0), 'ctt1' : (-3.0,3.0), + 'ctt' : (-3.0,3.0), } @@ -454,6 +513,7 @@ def gridScan(self, name='.test', batch='', freeze=False, scan_params=['ctW','ctZ if not freeze: wall_time /= 2 # profiled scans take longer, so submit less points per job if batch=='crab': args.extend(['--job-mode','crab3','--task-name',name.replace('.',''),'--custom-crab','custom_crab.py','--split-points',str(int(round(wall_time*point_scale)))]) if batch=='condor' and freeze==False and points>3000: args.extend(['--job-mode','condor','--task-name',name.replace('.',''),'--split-points','3000','--dry-run']) + elif batch=='condor' and freeze==True and points>3000: args.extend(['--job-mode','condor','--task-name',name.replace('.',''),'--split-points','300','--dry-run']) elif batch=='condor' and freeze==False: args.extend(['--job-mode','condor','--task-name',name.replace('.',''),'--split-points','10','--dry-run']) elif batch=='condor': args.extend(['--job-mode','condor','--task-name',name.replace('.',''),'--split-points','10','--dry-run']) logging.info(' '.join(args)) @@ -573,8 +633,9 @@ def getBestValues1DEFT(self, basename, wcs=[]): return ','.join(startValues) - def retrieveGridScan(self, name='.test', batch='crab', user='apiccine'):#getpass.getuser()): + def retrieveGridScan(self, name='.test', batch='crab', user=getpass.getuser(), remove=True): ### Retrieves finished grid jobs, extracts, and hadd's into a single file ### + user = user[:-1] if user[-1].isdigit() else user # Remove trailing number to match CRAB username taskname = name.replace('.','') logging.info("Retrieving gridScan files. Task name: "+taskname) @@ -604,7 +665,7 @@ def retrieveGridScan(self, name='.test', batch='crab', user='apiccine'):#getpass if tarfile.endswith('.tar'): print(tarfiles[0]+'/'+tarfile) sp.call(['tar', '-xf', tarfiles[0]+'/'+tarfile,'-C', taskname+'tmp']) - haddargs = ['hadd','-f','../fit_files/higgsCombine'+name+'.MultiDimFit.root']+['{}tmp/{}'.format(taskname,rootfile) for rootfile in os.listdir(taskname+'tmp') if rootfile.endswith('.root')] + haddargs = ['hadd','-f','-j','../fit_files/higgsCombine'+name+'.MultiDimFit.root']+['{}tmp/{}'.format(taskname,rootfile) for rootfile in os.listdir(taskname+'tmp') if rootfile.endswith('.root')] process = sp.Popen(haddargs, stdout=sp.PIPE, stderr=sp.PIPE) with process.stdout,process.stderr: self.log_subprocess_output(process.stdout,'info') @@ -619,19 +680,21 @@ def retrieveGridScan(self, name='.test', batch='crab', user='apiccine'):#getpass logging.info("No files to hadd. Returning.") return #haddargs = ['hadd','-f','higgsCombine'+name+'.MultiDimFit.root']+sorted(glob.glob('higgsCombine{}.POINTS*.root'.format(name))) - haddargs = ['hadd','-f','-k','../fit_files/higgsCombine'+name+'.MultiDimFit.root']+sorted(glob.glob('higgsCombine{}.POINTS*.root'.format(name))) - print((['hadd','-f','../../fit_files/higgsCombine'+name+'.MultiDimFit.root']+sorted(glob.glob('higgsCombine{}.POINTS*.root'.format(name))))) + haddargs = ['hadd','-f','-j','-k','../fit_files/higgsCombine'+name+'.MultiDimFit.root']+sorted(glob.glob('higgsCombine{}.POINTS*.root'.format(name))) + print(' '.join(haddargs)) + #print((['hadd','-f','../../fit_files/higgsCombine'+name+'.MultiDimFit.root']+sorted(glob.glob('higgsCombine{}.POINTS*.root'.format(name))))) process = sp.Popen(haddargs, stdout=sp.PIPE, stderr=sp.PIPE) with process.stdout,process.stderr: self.log_subprocess_output(process.stdout,'info') self.log_subprocess_output(process.stderr,'err') process.wait() - for rootfile in glob.glob('higgsCombine{}.POINTS*.root'.format(name)): - os.remove(rootfile) - if os.path.isfile('condor_{}.sh'.format(name.replace('.',''))): - os.rename('condor_{}.sh'.format(name.replace('.','')),'condor{0}/condor_{0}.sh'.format(name)) - if os.path.isfile('condor_{}.sub'.format(name.replace('.',''))): - os.rename('condor_{}.sub'.format(name.replace('.','')),'condor{0}/condor_{0}.sub'.format(name)) + if remove: + for rootfile in glob.glob('higgsCombine{}.POINTS*.root'.format(name)): + os.remove(rootfile) + if os.path.isfile('condor_{}.sh'.format(name.replace('.',''))): + os.rename('condor_{}.sh'.format(name.replace('.','')),'condor{0}/condor_{0}.sh'.format(name)) + if os.path.isfile('condor_{}.sub'.format(name.replace('.',''))): + os.rename('condor_{}.sub'.format(name.replace('.','')),'condor{0}/condor_{0}.sub'.format(name)) def submitEFTWilks(self, name='.test', limits='/afs/crc.nd.edu/user/b/byates2/Public/wc_top22006_a24_prof_2sigma.json', workspace='ptz-lj0pt_fullR2_anatest24v01_withAutostats_withSys.root', doBest=False, asimov=False, fixed=False, wc=None, sig=0, batch='condor'): ''' @@ -874,14 +937,24 @@ def batch1DScanEFT(self, basename='.test', batch='crab', freeze=False, scan_wcs= else: masks = ','.join(mask) mask = [] - self.gridScan('{}.{}'.format(basename,wc), batch, freeze, [wc], [wcs for wcs in self.wcs if wcs != wc], points, ['--setParameterRanges {}={},{}'.format(wc,wc_ranges[wc][0],wc_ranges[wc][1])]+zero_ignore+freeze_ignore+other+['--setParameters', params+','+masks], mask, mask_syst, workspace) + other_ranges = ':'.join(['{}=0,0'.format(wc_other) for wc_other in self.wcs if wc != wc_other]) + ranges = ['--setParameterRanges {}={},{}'.format(wc,wc_ranges[wc][0],wc_ranges[wc][1])] + if freeze: + ''' + This fixes an issue with the rotaiton scheme. + If the new WCs (e.g. ctGRe) are included for dim6top rotations only (i.e. you don't have any pure SMEFTsim bins), + then these are "internal parameters" and wind up getting set to their best fit values even when fronzen. + To fix this, we set their ranges to [0,0] to force combine to fix them to the SM only when it should be a frozen WC + ''' + ranges = ['--setParameterRanges {}={},{}'.format(wc,wc_ranges[wc][0],wc_ranges[wc][1]) + ':' + other_ranges] + self.gridScan('{}.{}'.format(basename,wc), batch, freeze, [wc], [wcs for wcs in self.wcs if wcs != wc], points, ranges+zero_ignore+freeze_ignore+other+['--setParameters', params+','+masks], mask, mask_syst, workspace) ''' example: `fitter.batch2DScanEFT('.test.ctZ', batch='crab', wcs=['ctZ'], workspace='wps_njet_runII.root')` example: `fitter.batch2DScanEFT('.test.ctZ', batch='crab', wcs='ctZ', workspace='wps_njet_runII.root')` ''' - def batch2DScanEFT(self, basename='.EFT.gridScan', batch='crab', freeze=False, points=90000, allPairs=False, other=[], mask=[], mask_syst=[], wcs=[], workspace='EFTWorkspace.root', differential=None): + def batch2DScanEFT(self, basename='.EFT.gridScan', batch='crab', freeze=False, points=90000, allPairs=False, other=[], mask=[], mask_syst=[], wcs=[], workspace='EFTWorkspace.root', differential=None, skip_internal=True): ### For pairs of wcs, runs deltaNLL Scan in two wcs using CRAB or Condor ### if differential is None: differential = 'njets' not in workspace @@ -893,12 +966,17 @@ def batch2DScanEFT(self, basename='.EFT.gridScan', batch='crab', freeze=False, p # Use EVERY combination of wcs if allPairs: scan_wcs = self.wcs + if wcs is not []: scan_wcs = wcs params = ','.join(['{}=0'.format(wc) for wc in scan_wcs]) for wcs in itertools.combinations(scan_wcs,2): wcs_tracked = [wc for wc in self.wcs if wc not in wcs] #print(pois, wcs_tracked) - self.gridScan(name='{}.{}{}'.format(basename,wcs[0],wcs[1]), batch=batch, freeze=freeze, scan_params=list(wcs), params_tracked=wcs_tracked, points=points, other=['--setParameterRanges {}={},{}:{}={},{}'.format(wcs[0],wc_ranges[wcs[0]][0],wc_ranges[wcs[0]][1],wcs[1],wc_ranges[wcs[1]][0],wc_ranges[wcs[1]][1])]+other+['--setParameters', params], mask=mask, mask_syst=mask_syst, workspace=workspace) + other_ranges = ':'.join(['{}=0,0'.format(wc_other) for wc_other in wcs_tracked]) + ranges = ['--setParameterRanges {}={},{}:{}={},{}'.format(wcs[0],wc_ranges[wcs[0]][0],wc_ranges[wcs[0]][1],wcs[1],wc_ranges[wcs[1]][0],wc_ranges[wcs[1]][1])] + if freeze and skip_internal: + ranges = ['--setParameterRanges {}={},{}:{}={},{}'.format(wcs[0],wc_ranges[wcs[0]][0],wc_ranges[wcs[0]][1],wcs[1],wc_ranges[wcs[1]][0],wc_ranges[wcs[1]][1]) + other_ranges] + self.gridScan(name='{}.{}{}'.format(basename,wcs[0],wcs[1]), batch=batch, freeze=freeze, scan_params=list(wcs), params_tracked=wcs_tracked, points=points, other=ranges+other+['--setParameters', params], mask=mask, mask_syst=mask_syst, workspace=workspace) # Use each wc only once if not allPairs: @@ -916,8 +994,12 @@ def batch2DScanEFT(self, basename='.EFT.gridScan', batch='crab', freeze=False, p params = ','.join(['{}=0'.format(w) for wc in scan_wcs for w in wc]) for wcs in scan_wcs: wcs_tracked = [wc for wc in self.wcs if wc not in wcs] + other_ranges = ':'.join(['{}=0,0'.format(wc_other) for wc_other in wcs_tracked]) + ranges = ['--setParameterRanges {}={},{}:{}={},{}'.format(wcs[0],wc_ranges[wcs[0]][0],wc_ranges[wcs[0]][1],wcs[1],wc_ranges[wcs[1]][0],wc_ranges[wcs[1]][1])] + if freeze and skip_internal: + ranges = ['--setParameterRanges {}={},{}:{}={},{}'.format(wcs[0],wc_ranges[wcs[0]][0],wc_ranges[wcs[0]][1],wcs[1],wc_ranges[wcs[1]][0],wc_ranges[wcs[1]][1]) + other_ranges] #print(scan_wcs, wcs_tracked) - self.gridScan(name='{}.{}{}'.format(basename,wcs[0],wcs[1]), batch=batch, freeze=freeze, scan_params=list(wcs), params_tracked=wcs_tracked, points=points, other=['--setParameterRanges {}={},{}:{}={},{}'.format(wcs[0],wc_ranges[wcs[0]][0],wc_ranges[wcs[0]][1],wcs[1],wc_ranges[wcs[1]][0],wc_ranges[wcs[1]][1])]+other+['--setParameters', params], mask=mask, mask_syst=mask_syst, workspace=workspace) + self.gridScan(name='{}.{}{}'.format(basename,wcs[0],wcs[1]), batch=batch, freeze=freeze, scan_params=list(wcs), params_tracked=wcs_tracked, points=points, other=ranges+other+['--setParameters', params], mask=mask, mask_syst=mask_syst, workspace=workspace) def batch3DScanEFT(self, basename='.EFT.gridScan', batch='crab', freeze=False, points=27000000, allPairs=False, other=[], wc_triplet=[], mask=[], mask_syst=[]): ### For pairs of wcs, runs deltaNLL Scan in two wcs using CRAB or Condor ### @@ -978,21 +1060,21 @@ def batchResubmit2DScansEFT(self, basename='.EFT.gridScan', allPairs=False): self.log_subprocess_output(process.stderr,'err') process.wait() - def batchRetrieve1DScansEFT(self, basename='.test', batch='crab', scan_wcs=[]): + def batchRetrieve1DScansEFT(self, basename='.test', batch='crab', scan_wcs=[], remove=True): ### For each wc, retrieves finished 1D deltaNLL grid jobs, extracts, and hadd's into a single file ### if not scan_wcs: scan_wcs = self.wcs for wc in scan_wcs: - self.retrieveGridScan('{}.{}'.format(basename,wc),batch) + self.retrieveGridScan('{}.{}'.format(basename,wc),batch, remove=remove) - def batchRetrieve2DScansEFT(self, basename='.EFT.gridScan', batch='crab', allPairs=False, wcs=[]): + def batchRetrieve2DScansEFT(self, basename='.EFT.gridScan', batch='crab', allPairs=False, wcs=[], remove=True): ### For pairs of wcs, retrieves finished grid jobs, extracts, and hadd's into a single file ### # Use EVERY combination of wcs if allPairs: scan_wcs = self.wcs for wcs in itertools.combinations(scan_wcs,2): - self.retrieveGridScan('{}.{}{}'.format(basename,wcs[0],wcs[1]),batch) + self.retrieveGridScan('{}.{}{}'.format(basename,wcs[0],wcs[1]),batch,remove=remove) # Use each wc only once if not allPairs: diff --git a/Fitter/scripts/EFTPlotter.py b/Fitter/scripts/EFTPlotter.py index b356482..37c362f 100644 --- a/Fitter/scripts/EFTPlotter.py +++ b/Fitter/scripts/EFTPlotter.py @@ -19,9 +19,12 @@ def __init__(self,wc_ranges=None): self.ContourHelper = ContourHelper() self.SMMus = ['mu_ttll','mu_ttlnu','mu_ttH','mu_tllq'] - self.wcs = ['ctW','ctZ','ctp','cpQM','ctG','cbW','cpQ3','cptb','cpt','cQl3i','cQlMi','cQei','ctli','ctei','ctlSi','ctlTi'] - self.wcs_pairs = [('ctZ','ctW'),('ctp','cpt'),('ctlSi','ctli'),('cptb','cQl3i'),('ctG','cpQM'),('ctei','ctlTi'),('cQlMi','cQei'),('cpQ3','cbW')] - self.wcs = ['ctW','ctZ','ctp','cpQM','ctG','cbW','cpQ3','cptb','cpt','cQl3i','cQlMi','cQei','ctli','ctei','ctlSi','ctlTi', 'cQq13', 'cQq83', 'cQq11', 'ctq1', 'cQq81', 'ctq8', 'cQt1', 'cQt8', 'cQQ1', 'ctt1'] #TOP-22-006 + self.wcs = ['ctW', 'ctZ', 'ctp', 'cpQM', 'ctG', 'cbW', 'cpQ3', 'cptb', 'cpt', 'cQl3i', 'cQlMi', 'cQei', 'ctli', 'ctei', 'ctlSi', 'ctlTi', 'cQq13', 'cQq83', 'cQq11', 'ctq1', 'cQq81', 'ctq8', 'ctt1', 'cQQ1', 'cQt8', 'cQt1'] + self.rotate = {'ctW': 'ctWRe', 'ctZ': 'ctBRe', 'ctp': 'ctHRe', 'cpQM': 'cHQ1', 'ctG': 'ctGRe', 'cbW': 'cbW', 'cpQ3': 'cHQ3', 'cptb': 'cHtbRe', 'cpt': 'cHt', 'cQl3i': 'cQl3i', 'cQlMi': 'cQl1i', 'cQei': 'cQei', 'ctli': 'ctli', 'ctei': 'ctei', 'ctlSi': 'cleQt1Rei', 'ctlTi': 'cleQt3Rei', 'cQq13': 'cQj31', 'cQq83': 'cQj38', 'cQq11': 'cQj11', 'ctq1': 'ctj1', 'cQq81': 'cQj18', 'ctq8': 'ctj8', 'ctt1': 'ctt', 'cQQ1': 'cQQ1', 'cQt8': 'cQt8', 'cQt1': 'cQt1'} + #self.rotate = {'ctW': 'ctBRe', 'ctZ': 'ctZRe', 'ctp': 'ctH', 'cpQM': 'cHQ1', 'ctG': 'ctGRe', 'cbW': 'cbW', 'cpQ3': 'cHQ3', 'cptb': 'cHtbRe', 'cpt': 'cHt', 'cQl3i': 'cQl3i', 'cQlMi': 'cQl1i', 'cQei': 'cQei', 'ctli': 'ctli', 'ctei': 'ctei', 'ctlSi': 'cleQt1Rei', 'ctlTi': 'cleQt3Rei'} + # TODO update others like 'ctp' -> 'ctH' when ws is ready + self.wcs_pairs = [('ctBRe','ctWRe'),('ctp','cpt'),('ctlSi','ctli'),('cptb','cQl3i'),('ctGRe','cpQM'),('ctei','ctlTi'),('cQlMi','cQei'),('cpQ3','cbWRe')] + #self.wcs = ['ctWRe','ctBRe','ctp','cpQM','ctGRe','cbW','cpQ3','cptb','cpt','cQl3i','cQlMi','cQei','ctli','ctei','ctlSi','ctlTi', 'cQq13', 'cQq83', 'cQq11', 'ctq1', 'cQq81', 'ctq8', 'cQt1', 'cQt8', 'cQQ1', 'ctt1'] #TOP-22-006 #self.wcs_pairs = [('ctW','ctG'),('ctZ','ctG'),('ctp','ctG'),('cpQM','ctG'),('cbW','ctG'),('cpQ3','ctG'),('cptb','ctG'),('cpt','ctG'),('cQl3i','ctG'),('cQlMi','ctG'),('cQei','ctG'),('ctli','ctG'),('ctei','ctG'),('ctlSi','ctG'),('ctlTi','ctG')] # Set the WC ranges (if not specified, just use some numbers that generally work for njets) self.wc_ranges = { @@ -29,28 +32,46 @@ def __init__(self,wc_ranges=None): 'cQei' : (-4.0,4.0), 'cQl3i': (-5.5,5.5), 'cQlMi': (-4.0,4.0), + 'cQl1i': (-4.0,4.0), 'cQq11': (-0.7,0.7), + 'cQj11': (-0.7,0.7), 'cQq13': (-0.35,0.35), + 'cQj31': (-0.35,0.35), 'cQq81': (-1.7,1.5), 'cQq83': (-0.6,0.6), + 'cQj38': (-0.6,0.6), 'cQt1' : (-6.0,6.0), 'cQt8' : (-10.0,10.0), 'cbW' : (-3.0,3.0), 'cpQ3' : (-4.0,4.0), - 'cpQM' : (-15.0,20.0), + 'chQ3' : (-4.0,4.0), + 'cpQM' : (-2.0,2.0), + 'cHQ1' : (-15.0,20.0), 'cpt' : (-15.0,15.0), + 'cpt' : (-5.0,5.0), + 'cHt' : (-15.0,15.0), 'cptb' : (-9.0,9.0), + 'cHtbRe' : (-9.0,9.0), 'ctG' : (-0.8,0.8), + 'ctGRe' : (-0.8,0.8), 'ctW' : (-1.5,1.5), - 'ctZ' : (-2.0,2.0), + 'ctWRe' : (-1.5,1.5), + 'ctZ' : (-3.0,3.0), + 'ctBRe' : (-3.0,3.0), 'ctei' : (-4.0,4.0), 'ctlSi': (-5.0,5.0), + 'cleQt1Rei': (-5.0,5.0), 'ctlTi': (-0.9,0.9), + 'cleQt3Rei': (-0.9,0.9), 'ctli' : (-4.0,4.0), 'ctp' : (-15.0,40.0), + 'ctHRe' : (-15.0,40.0), 'ctq1' : (-0.6,0.6), + 'ctj1' : (-0.6,0.6), 'ctq8' : (-1.4,1.4), + 'ctj8' : (-1.4,1.4), 'ctt1' : (-2.6,2.6), + 'ctt' : (-2.6,2.6), } if wc_ranges is not None: self.wc_ranges = wc_ranges @@ -60,10 +81,13 @@ def __init__(self,wc_ranges=None): self.histosFileName = 'Histos.root' self.texdic = { 'ctW': '\it{c}_{\mathrm{tW}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', + 'ctWRe': '\it{c}_{\mathrm{tW}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', 'ctZ': '\it{c}_{\mathrm{tZ}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', + 'ctBRe': '\it{c}_{\mathrm{tZ}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', 'ctp': '\it{c}_{\mathrm{t} \\varphi}/\mathrm{\Lambda^{2} [TeV^{-2}]}', 'cpQM': '\it{c}^{-}_{\\varphi \mathrm{Q}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', 'ctG': '\it{c}_{\mathrm{tG}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', + 'ctGRe': '\it{c}_{\mathrm{tG}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', 'cbW': '\it{c}_{\mathrm{bW}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', 'cpQ3': '\it{c}^{3}_{\\varphi \mathrm{Q}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', 'cptb': '\it{c}_{\\varphi \mathrm{tb}}/\mathrm{\Lambda^{2} [TeV^{-2}]}', @@ -88,10 +112,13 @@ def __init__(self,wc_ranges=None): } self.texdicfrac = { 'ctW': '\it{c}_{\mathrm{tW}}', + 'ctWRe': '\it{c}_{\mathrm{tWRe}}', 'ctZ': '\it{c}_{\mathrm{tZ}}', + 'ctBRe': '\it{c}_{\mathrm{tBRe}}', 'ctp': '\it{c}_{\mathrm{t} \\varphi}', 'cpQM': '\it{c}^{-}_{\\varphi \mathrm{Q}}', 'ctG': '\it{c}_{\mathrm{tG}}', + 'ctGRe': '\it{c}_{\mathrm{tGRe}}', 'cbW': '\it{c}_{\mathrm{bW}}', 'cpQ3': '\it{c}^{3}_{\\varphi \mathrm{Q}}', 'cptb': '\it{c}_{\\varphi \mathrm{tb}}', @@ -116,10 +143,13 @@ def __init__(self,wc_ranges=None): } self.texdicmacro = { 'ctW': '\ctW', + 'ctWRe': '\ctWRe', 'ctZ': '\ctZ', + 'ctBRe': '\ctBRe', 'ctp': '\ctp', 'cpQM': '\cpQM', 'ctG': '\ctG', + 'ctGRe': '\ctGRe', 'cbW': '\cbW', 'cpQ3': '\cpQa', 'cptb': '\cptb', @@ -187,6 +217,9 @@ def GetWCsNLLFromRoot(self,base_name_lst,wc,unique=False,**kwargs): graphwcs = [] graphnlls = [] for name in base_name_lst: + if not os.path.exists('{}/higgsCombine{}.MultiDimFit.root'.format(dir_path,name)) and wc in self.rotate and os.path.exists('{}/higgsCombine{}.MultiDimFit.root'.format(dir_path,name.replace(wc, self.rotate[wc]))): + name = name.replace(wc, self.rotate[wc]) + wc = self.rotate[wc] if not os.path.exists('{}/higgsCombine{}.MultiDimFit.root'.format(dir_path,name)): logging.error("File {}/higgsCombine{}.MultiDimFit.root does not exist!".format(dir_path,name)) return [graphwcs,graphnlls] @@ -360,14 +393,52 @@ def OverlayLLPlot1DEFT(self,**kwargs): if not wc: logging.error("No wc specified!") return - for name1 in name1_lst: + + # Hacks for rotations + # This allows the new SMEFTsim WCs to play nice with older TOP-22-006 files. + wc1 = wc + wc2 = wc + for iname,name1 in enumerate(name1_lst): + old = None + if not os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d1,name1,pf1)) and wc in self.rotate and os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d1,name1.replace(wc, self.rotate[wc]),pf1)): + name1 = name1.replace(wc, self.rotate[wc]) + old = wc + wc = self.rotate[wc] + #if not os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d1,name1,pf1)): + # if 'B' in name1: + # name1 = name1.replace('ctBRe', 'ctZ') + # wc = 'ctZ' + # wc1 = 'ctZ' + # else: + # name1 = name1.replace('Re', '') + # wc = wc.replace('Re', '') + # wc1 = wc.replace('Re', '') + # print(f'Trying {name1}') + # name1_lst[iname] = name1 if not os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d1,name1,pf1)): logging.error("File higgsCombine{}.MultiDimFit{}.root does not exist!".format(name1,pf1)) return - for name2 in name1_lst: + for iname,name2 in enumerate(name2_lst): + if not os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d2,name2,pf2)) and wc in self.rotate and os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d2,name2.replace(wc, self.rotate[wc]),pf2)): + name2 = name2.replace(wc, self.rotate[wc]) + old = wc + wc = self.rotate[wc] + #if not os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d2,name2,pf1)): + # if 'B' in name2: + # name2 = name1.replace('ctBRe', 'ctZ') + # wc = 'ctZ' + # wc2 = 'ctZ' + # else: + # name2 = name2.replace('Re', '') + # wc = wc.replace('Re', '') + # wc2 = wc.replace('Re', '') + # print(f'Trying {name2}') + # name2_lst[iname] = name2 if not os.path.exists('{}/higgsCombine{}.MultiDimFit{}.root'.format(d2,name2,pf2)): logging.error("File higgsCombine{}.MultiDimFit{}.root does not exist!".format(name2,pf2)) return + if old is not None: + wc = old ROOT.gROOT.SetBatch(True) @@ -377,8 +448,8 @@ def OverlayLLPlot1DEFT(self,**kwargs): p1.cd() # Get coordinates for TGraphs - graph1wcs,graph1nlls = self.GetWCsNLLFromRoot(name1_lst,wc,unique=True,dir_path=d1) - graph2wcs,graph2nlls = self.GetWCsNLLFromRoot(name2_lst,wc,unique=True,dir_path=d2) + graph1wcs,graph1nlls = self.GetWCsNLLFromRoot(name1_lst,wc1,unique=True,dir_path=d1) + graph2wcs,graph2nlls = self.GetWCsNLLFromRoot(name2_lst,wc2,unique=True,dir_path=d2) # Rezero the y axis and make the tgraphs #zero = graph1nlls.index(0) @@ -412,8 +483,8 @@ def OverlayLLPlot1DEFT(self,**kwargs): #multigraph.GetXaxis().SetNdivisions(7) # Squeeze X down to whatever range captures the float points - xmin = self.wc_ranges[wc][0] - xmax = self.wc_ranges[wc][1] + xmin = self.wc_ranges[wc1][0] + xmax = self.wc_ranges[wc1][1] for idx in range(graph1.GetN()): if graph1.GetY()[idx] < ceiling and graph1.GetX()[idx] < xmin: xmin = graph1.GetX()[idx] @@ -1853,7 +1924,6 @@ def getIntervalFits(self,**kwargs): ### Use 1D scans instead of regular MultiDimFit ### if not params: params = self.wcs - ROOT.gROOT.SetBatch(True) @@ -1863,9 +1933,28 @@ def getIntervalFits(self,**kwargs): # This is mostly used to compare TOP-19-001 to Run II, it will skip the 10 WCs not in TOP-19-001 only and set them to +/- 999 missing_wc = False + # Hacks for rotations + # This allows the new SMEFTsim WCs to play nice with older TOP-22-006 files. + old = None for basename in basename_lst: logging.debug("Obtaining result of scan: higgsCombine{}.{}.MultiDimFit{}.root".format(basename,param,postfix)) - fit_file = ROOT.TFile.Open('{}/higgsCombine{}.{}.MultiDimFit{}.root'.format(dir_path,basename,param,postfix)) + if not os.path.exists('{}/higgsCombine{}.{}.MultiDimFit{}.root'.format(dir_path,basename,param,postfix)) and os.path.exists('{}/higgsCombine{}.{}.MultiDimFit{}.root'.format(dir_path,basename,self.rotate[param],postfix)): + #basename = basename.replace(param, self.rotate[param]) + old = param + param = self.rotate[param] + #if not os.path.exists('{}/higgsCombine{}.{}.MultiDimFit{}.root'.format(dir_path,basename,param,postfix)): + # if 'B' in param: + # param = 'ctZ' + # else: + # param = param.replace('Re', '') + if os.path.exists('{}/higgsCombine{}.{}.MultiDimFit{}.root'.format(dir_path,basename,param,postfix)): + #fit_array.append([param,0,[-999 ],[999]]) + fit_file = ROOT.TFile.Open('{}/higgsCombine{}.{}.MultiDimFit{}.root'.format(dir_path,basename,param,postfix)) + else: + print('WARNING,', '{}/higgsCombine{}.{}.MultiDimFit{}.root'.format(dir_path,basename,param,postfix), 'does not exist! Setting to +/- 999') + missing_wc = True + if param not in map(itemgetter(0),fit_array): + fit_array.append([param,0,[-999 ],[999]]) try: _ = fit_file.Get('limit') fit_file.Close() @@ -1874,6 +1963,8 @@ def getIntervalFits(self,**kwargs): if param not in map(itemgetter(0),fit_array): fit_array.append([param,0,[-999 ],[999]]) if missing_wc: continue + if old is not None: + param = old basename_lst_with_wcs_appended = self.AppendStrToItemsInLst(basename_lst,"."+param) wc_values, nll_values = self.GetWCsNLLFromRoot(basename_lst_with_wcs_appended,param,unique=True,dir_path=dir_path) @@ -2077,7 +2168,7 @@ def BestScanPlot(self, basename_float_lst=[''], basename_freeze_lst=[''], final= print('\\end{table}') for idx,line in enumerate(fits_float): - if line[0]=='ctG': + if line[0]=='ctG' or line[0]=='ctGRe': line[0] = 'ctG#times5' line[1] = line[1]*5 line[2] = [val*5 for val in line[2]] @@ -2119,7 +2210,7 @@ def BestScanPlot(self, basename_float_lst=[''], basename_freeze_lst=[''], final= line[3] = [val/2 for val in line[3]] for idx,line in enumerate(fits_freeze): - if line[0]=='ctG': + if line[0]=='ctG' or line[0]=='ctGRe': line[0] = 'ctG#times5' line[1] = line[1]*5 line[2] = [val*5 for val in line[2]] @@ -2161,7 +2252,7 @@ def BestScanPlot(self, basename_float_lst=[''], basename_freeze_lst=[''], final= line[3] = [val/2 for val in line[3]] for idx,line in enumerate(fits_float1sigma): - if line[0]=='ctG': + if line[0]=='ctG' or line[0]=='ctGRe': line[0] = 'ctG#times5' line[1] = line[1]*5 line[2] = [val*5 for val in line[2]] @@ -2203,7 +2294,7 @@ def BestScanPlot(self, basename_float_lst=[''], basename_freeze_lst=[''], final= line[3] = [val/2 for val in line[3]] for idx,line in enumerate(fits_freeze1sigma): - if line[0]=='ctG': + if line[0]=='ctG' or line[0]=='ctGRe': line[0] = 'ctG#times5' line[1] = line[1]*5 line[2] = [val*5 for val in line[2]]