diff --git a/LHScans/plot_LLScan.py b/LHScans/plot_LLScan.py index 4fa3dd2..7e1a2c7 100644 --- a/LHScans/plot_LLScan.py +++ b/LHScans/plot_LLScan.py @@ -189,7 +189,7 @@ def BuildScan(scan, param, files, color, yvals, ycut): elif(obsName == 'mass4l_zzfloating'): label = 'm_{4l}' elif(obsName == 'njets_pt30_eta4p7'): label = 'N_{jet}, pT>30 GeV, |#eta|<4.7' elif(obsName == 'pT4l'): label = 'p_{T}^{H} (GeV)' -elif(obsName == 'pT4l_kL'): label = '' +elif(obsName == 'pT4l_kL'): label = 'k_{#lambda}' elif(obsName == 'rapidity4l'): label = '|y_{H}|' elif(obsName == 'costhetaZ1'): label = 'cos(#theta_{1})' elif(obsName == 'costhetaZ2'): label = 'cos(#theta_{2})' @@ -252,7 +252,7 @@ def BuildScan(scan, param, files, color, yvals, ycut): doubleDiff = True # _poi = 'SigmaBin' -_obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta2p5': 'NJ'} +_obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} if obsName not in _obsName: _obsName[obsName] = obsName # _poi = 'r_smH_'+_obsName[obsName]+'_' @@ -533,7 +533,7 @@ def BuildScan(scan, param, files, color, yvals, ycut): graphs[ig].SetLineWidth(3) if 'stat-only' in titles[ig]: graphs[ig].SetLineWidth(2) - graphs[ig].SetLineStyle(9) + graphs[ig].SetLineStyle(2) graphs[ig].SetTitle(titles[ig]) graphs[ig].Sort() @@ -643,6 +643,17 @@ def BuildScan(scan, param, files, color, yvals, ycut): exp_up_sys *= xsec['SigmaBin'+str(i)] exp_do_sys *= xsec['SigmaBin'+str(i)] + if opt.UNBLIND: + obs_nom = list(obs_nom) + obs_nom_stat = list(obs_nom_stat) + obs_nom[0] *= xsec['SigmaBin'+str(i)] + obs_nom[1] *= xsec['SigmaBin'+str(i)] + obs_nom[2] *= xsec['SigmaBin'+str(i)] + obs_nom_stat[1] *= xsec['SigmaBin'+str(i)] + obs_nom_stat[2] *= xsec['SigmaBin'+str(i)] + obs_up_sys *= xsec['SigmaBin'+str(i)] + obs_do_sys *= xsec['SigmaBin'+str(i)] + if(opt.UNBLIND): Text3 = TPaveText(0.18, 0.81,0.4,0.9,'brNDC') else: diff --git a/coefficients/JES/PrintJES.py b/coefficients/JES/PrintJES.py index d6f4f30..9efbf13 100644 --- a/coefficients/JES/PrintJES.py +++ b/coefficients/JES/PrintJES.py @@ -6,7 +6,7 @@ from binning import binning from tabulate import tabulate -print 'Welcome in RunJES!' +print 'Welcome in PrintJES!' def parseOptions(): @@ -65,6 +65,7 @@ def parseOptions(): if doubleDiff: nBins = len(obs_bins) else: nBins = len(obs_bins)-1 + # Tables with numerical values tables = {} inclusiveJES = {} @@ -78,50 +79,50 @@ def parseOptions(): # nominal_incl += evts['signal_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # up_incl += evts['signal_jesup_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # dn_incl += evts['signal_jesdn_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] - table.append(['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin), - evts_noWeight['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['signal_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['signal_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['signal_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['signal_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - JESNP['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)]]) + table.append(['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin), + evts_noWeight['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['signal_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['signal_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['signal_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['signal_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + JESNP['signal_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)]]) # nominal_incl += evts['qqzz_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # up_incl += evts['qqzz_jesup_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # dn_incl += evts['qqzz_jesdn_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] - table.append(['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin), - evts_noWeight['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['qqzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['qqzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['qqzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['qqzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - JESNP['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)]]) + table.append(['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin), + evts_noWeight['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['qqzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['qqzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['qqzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['qqzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + JESNP['qqzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)]]) # nominal_incl += evts['ggzz_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # up_incl += evts['ggzz_jesup_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # dn_incl += evts['ggzz_jesdn_'+jesName+'_'+year+'_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] - table.append(['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin), - evts_noWeight['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['ggzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['ggzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['ggzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['ggzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - JESNP['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)]]) + table.append(['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin), + evts_noWeight['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['ggzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['ggzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['ggzz_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['ggzz_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + JESNP['ggzz_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)]]) # nominal_incl += evts['ZX_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # up_incl += evts['ZX_jesup_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] # dn_incl += evts['ZX_jesdn_'+fState+'_'+obsname_out+'_recobin'+str(recobin)] - table.append(['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin), - evts_noWeight['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['ZX_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts_noWeight['ZX_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['ZX_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - evts['ZX_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)], - JESNP['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out+'_recobin'+str(recobin)]]) + table.append(['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin), + evts_noWeight['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['ZX_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts_noWeight['ZX_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['ZX_jesup_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + evts['ZX_jesdn_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)], + JESNP['ZX_'+jesName+'_'+fState+'_'+str(year)+'_'+obsname_out.replace('pT4l', 'ZZPt')+'_recobin'+str(recobin)]]) # table.append(['INCLUSIVE', nominal_incl, up_incl, dn_incl, str(round(dn_incl/nominal_incl,3))+'/'+str(round(up_incl/nominal_incl,3))]) table.append([]) diff --git a/coefficients/JES/RunJES.py b/coefficients/JES/RunJES.py index 29affbf..b60310e 100644 --- a/coefficients/JES/RunJES.py +++ b/coefficients/JES/RunJES.py @@ -20,6 +20,8 @@ print 'Welcome in RunJES!' jesNames = ['Total', 'Abs', 'Abs_year', 'BBEC1', 'BBEC1_year', 'EC2', 'EC2_year', 'FlavQCD', 'HF', 'HF_year', 'RelBal', 'RelSample_year'] +JESobservables = ['pTj1', 'pTHj', 'mHj', 'pTj2', 'mjj', 'absdetajj', 'dphijj', 'pTHjj', 'TCjmax', 'TBjmax', 'njets_pt30_eta4p7 vs pT4l', 'pTj1 vs pTj2', 'pT4l vs pTHj', 'TCjmax vs pT4l'] + def parseOptions(): diff --git a/coefficients/JES/zx.py b/coefficients/JES/zx.py index 50d4db2..4534872 100644 --- a/coefficients/JES/zx.py +++ b/coefficients/JES/zx.py @@ -136,26 +136,25 @@ def findFSZX(df): def comb(year): if year == 2016: cb_SS = np.array([ - 1.23628, # 4e - 0.95433, # 4mu - 1.0726, # 2e2mu - 1.0726, # 2mu2e + 1.175, # 4e + 0.975, # 4mu + 1.052, # 2e2mu + 1.148, # 2mu2e ]) elif year == 2017: cb_SS = np.array([ - 1.1934, # 4e - 0.99669, # 4mu - 1.0569, # 2e2mu - 1.0569, # 2mu2e + 1.094, # 4e + 0.948, # 4mu + 0.930, # 2e2mu + 1.139, # 2mu2e ]) else: cb_SS = np.array([ - 1.2087, # 4e - 0.9878, # 4mu - 1.0552, # 2e2mu - 1.0552, # 2mu2e + 1.157, # 4e + 0.974, # 4mu + 0.930, # 2e2mu + 1.143, # 2mu2e ]) - return cb_SS @@ -163,26 +162,25 @@ def comb(year): def ratio(year): if year == 2016: fs_ROS_SS = np.array([ - 1.00245, # 4e - 0.998863, # 4mu - 1.03338, # 2e2mu - 0.998852, # 2mu2e + 1.0039, # 4e + 0.999103, # 4mu + 1.0332, # 2e2mu + 1.00216, # 2mu2e ]) elif year == 2017: fs_ROS_SS = np.array([ - 1.01198, # 4e - 1.03949, # 4mu - 1.013128, # 2e2mu - 1.00257, # 2mu2e + 0.990314, # 4e + 1.02903, # 4mu + 1.0262, # 2e2mu + 1.00154, # 2mu2e ]) else: fs_ROS_SS = np.array([ - 1.00568, # 4e - 1.02926, # 4mu - 1.03226, # 2e2mu - 1.00432, # 2mu2e + 1.00322, # 4e + 1.0187, # 4mu + 1.04216, # 2e2mu + 0.996253, # 2mu2e ]) - return fs_ROS_SS def count_jets(pt,eta,phi,mass): @@ -218,7 +216,7 @@ def doZX(year, g_FR_mu_EB, g_FR_mu_EE, g_FR_e_EB, g_FR_e_EE, obs_reco,obs_reco_2 'helcosthetaZ1','helcosthetaZ2', 'helphi', 'costhetastar', 'phistarZ1', 'ZZPhi', 'pTj1', 'pTj2', 'absdetajj', 'JetEta','JetPhi','JetMass', - 'pTHj', 'pTHjj', 'mHj', 'mHjj', 'detajj', 'dphijj', 'mjj', 'njets_pt30_eta2p5', 'ZZy', + 'pTHj', 'pTHjj', 'mHj', 'mHjj', 'detajj', 'dphijj', 'mjj', 'njets_pt30_eta4p7', 'ZZy', 'D0m', 'Dcp', 'D0hp', 'Dint', 'DL1', 'DL1int', 'DL1Zg', 'DL1Zgint', 'TCjmax', 'TBjmax'] keyZX = 'CRZLL' @@ -253,29 +251,6 @@ def doZX(year, g_FR_mu_EB, g_FR_mu_EE, g_FR_e_EB, g_FR_e_EE, obs_reco,obs_reco_2 dfZX['j2_jesdn_'+i] = [add_subleadjet(row[0],row[1],row[2],row[3],row[4]) for row in dfZX[['JetPt_JESDown_'+i,'JetEta','JetPhi','JetMass','j1_jesdn_'+i]].values] dfZX['Higgs'] = [tetra_Higgs(row[0],row[1],row[2],row[3]) for row in dfZX[['ZZMass', 'ZZEta', 'ZZPhi', 'ZZPt']].values] - # dfZX['pTj1_jesup'] = [x.Pt() for x in dfZX['j1_jesup']] - # dfZX['pTj1_jesdn'] = [x.Pt() for x in dfZX['j1_jesdn']] - # dfZX['pTj2_jesup'] = [x.Pt() for x in dfZX['j2_jesup']] - # dfZX['pTj2_jesdn'] = [x.Pt() for x in dfZX['j2_jesdn']] - # dfZX['njets_pt30_eta2p5_jesup'] = [count_jets(row[0],row[1],row[2],row[3]) for row in dfZX[['JetPt_JESUp','JetEta','JetPhi','JetMass']].values] - # dfZX['njets_pt30_eta2p5_jesdn'] = [count_jets(row[0],row[1],row[2],row[3]) for row in dfZX[['JetPt_JESDown','JetEta','JetPhi','JetMass']].values] - # dfZX['mjj_jesup'] = [(j1+j2).M() if j2.Pt()>0 else -1 for j1,j2 in zip(dfZX['j1_jesup'],dfZX['j2_jesup'])] - # dfZX['mjj_jesdn'] = [(j1+j2).M() if j2.Pt()>0 else -1 for j1,j2 in zip(dfZX['j1_jesdn'],dfZX['j2_jesdn'])] - # dfZX['pTHj_jesup'] = [(H+j1).Pt() if j1.Pt()>0 else -1 for H,j1 in zip(dfZX['Higgs'],dfZX['j1_jesup'])] - # dfZX['pTHj_jesdn'] = [(H+j1).Pt() if j1.Pt()>0 else -1 for H,j1 in zip(dfZX['Higgs'],dfZX['j1_jesdn'])] - # dfZX['pTHjj_jesup'] = [(row[0]+row[1]+row[2]).Pt() if row[2].Pt()>0 else -1 for row in dfZX[['Higgs','j1_jesup','j2_jesup']].values] - # dfZX['pTHjj_jesdn'] = [(row[0]+row[1]+row[2]).Pt() if row[2].Pt()>0 else -1 for row in dfZX[['Higgs','j1_jesdn','j2_jesdn']].values] - # dfZX['mHj_jesup'] = [(H+j1).M() if j1.Pt()>0 else -1 for H,j1 in zip(dfZX['Higgs'],dfZX['j1_jesup'])] - # dfZX['mHj_jesdn'] = [(H+j1).M() if j1.Pt()>0 else -1 for H,j1 in zip(dfZX['Higgs'],dfZX['j1_jesdn'])] - # dfZX['mHjj_jesup'] = [(row[0]+row[1]+row[2]).M() if row[2].Pt()>0 else -1 for row in dfZX[['Higgs','j1_jesup','j2_jesup']].values] - # dfZX['mHjj_jesdn'] = [(row[0]+row[1]+row[2]).M() if row[2].Pt()>0 else -1 for row in dfZX[['Higgs','j1_jesdn','j2_jesdn']].values] - # dfZX['TCjmax_jesup'] = [tc(row[0],row[1],row[2],row[3],row[4]) for row in dfZX[['JetPt_JESUp','JetEta','JetPhi','JetMass','Higgs']].values] - # dfZX['TCjmax_jesdn'] = [tc(row[0],row[1],row[2],row[3],row[4]) for row in dfZX[['JetPt_JESDown','JetEta','JetPhi','JetMass','Higgs']].values] - # dfZX['TBjmax_jesup'] = [tb(row[0],row[1],row[2],row[3],row[4]) for row in dfZX[['JetPt_JESUp','JetEta','JetPhi','JetMass','Higgs']].values] - # dfZX['TBjmax_jesdn'] = [tb(row[0],row[1],row[2],row[3],row[4]) for row in dfZX[['JetPt_JESDown','JetEta','JetPhi','JetMass','Higgs']].values] - # dfZX['absdetajj_jesup'] = [abs(j1.Eta()-j2.Eta()) if j2.Pt()>0 else -1 for j1,j2 in zip(dfZX['j1_jesup'],dfZX['j2_jesup'])] - # dfZX['absdetajj_jesdn'] = [abs(j1.Eta()-j2.Eta()) if j2.Pt()>0 else -1 for j1,j2 in zip(dfZX['j1_jesdn'],dfZX['j2_jesdn'])] - for i in jesNames: if obs_reco == 'pTj1' or obs_reco_2nd == 'pTj1': dfZX['pTj1_jesup_'+i] = [x.Pt() for x in dfZX['j1_jesup_'+i]] @@ -304,6 +279,9 @@ def doZX(year, g_FR_mu_EB, g_FR_mu_EE, g_FR_e_EB, g_FR_e_EE, obs_reco,obs_reco_2 if obs_reco == 'absdetajj' or obs_reco_2nd == 'absdetajj': dfZX['absdetajj_jesup_'+i] = [abs(j1.Eta()-j2.Eta()) if j2.Pt()>0 else -1 for j1,j2 in zip(dfZX['j1_jesup_'+i],dfZX['j2_jesup_'+i])] dfZX['absdetajj_jesdn_'+i] = [abs(j1.Eta()-j2.Eta()) if j2.Pt()>0 else -1 for j1,j2 in zip(dfZX['j1_jesdn_'+i],dfZX['j2_jesdn_'+i])] + if obs_reco == 'dphijj' or obs_reco_2nd == 'dphijj': + dfZX['dphijj_jesup_'+i] = [math.atan2(math.sin(j1.Phi()-j2.Phi()), math.cos(j1.Phi()-j2.Phi())) if j2.Pt()>0 else -99 for j1,j2 in zip(dfZX['j1_jesup_'+i],dfZX['j2_jesup_'+i])] + dfZX['dphijj_jesdn_'+i] = [math.atan2(math.sin(j1.Phi()-j2.Phi()), math.cos(j1.Phi()-j2.Phi())) if j2.Pt()>0 else -99 for j1,j2 in zip(dfZX['j1_jesdn_'+i],dfZX['j2_jesdn_'+i])] if obs_reco == 'ZZPt' or obs_reco_2nd == 'ZZPt': dfZX['ZZPt_jesup_'+i] = dfZX['ZZPt'] dfZX['ZZPt_jesdn_'+i] = dfZX['ZZPt'] diff --git a/coefficients/RunCoefficients.py b/coefficients/RunCoefficients.py index 57d282b..48016f8 100644 --- a/coefficients/RunCoefficients.py +++ b/coefficients/RunCoefficients.py @@ -103,7 +103,10 @@ def prepareTrees(year): # fname += '/'+signal+'ext/'+signal+'ext_reducedTree_MC_'+str(year)+'.root' # else: if opt.AC==True or opt.AC_ONLYACC==True: - fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year[:-4])+'.root' #[:-4] is to cut 'post' from year + if year == '2016post': + fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year[:-4])+'.root' #[:-4] is to cut 'post' from year + else: + fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year)+'.root' #[:-4] is to cut 'post' from year else: fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year)+'.root' print fname @@ -208,7 +211,10 @@ def generators(year): # fname += '/'+signal+'ext/'+signal+'ext_reducedTree_MC_'+str(year)+'.root' # else: if opt.AC==True or opt.AC_ONLYACC==True: - fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year[:-4])+'.root' #[:-4] is to cut 'post' from year + if year == '2016post': + fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year[:-4])+'.root' #[:-4] is to cut 'post' from year + else: + fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year)+'.root' #[:-4] is to cut 'post' from year else: fname += '/'+signal+'/'+signal+'_reducedTree_MC_'+str(year)+'.root' print fname diff --git a/coefficients/RunInterpolation.py b/coefficients/RunInterpolation.py index 1992561..f4d4d9b 100644 --- a/coefficients/RunInterpolation.py +++ b/coefficients/RunInterpolation.py @@ -109,7 +109,7 @@ def parseOptions(): extrap_binfrac_wrongfrac = {} nBins = len(obs_bins) if not doubleDiff: nBins = len(obs_bins)-1 #In case of 1D measurement the number of bins is -1 the length of obs_bins(=bin boundaries) -for channel in ['2e2mu', '4e', '4mu']: +for channel in ['2e2mu', '4e', '4mu']:#, '4l']: for genBin in range(nBins): for recoBin in range(nBins): fig,axs = plt.subplots(2, 3, figsize=(30,10), dpi=80) @@ -217,25 +217,35 @@ def parseOptions(): if opt.NNLOPS: continue - diff = [(extrap_acc[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] - acc[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)]) / acc[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] - diff = np.mean(diff) - extrap_acc['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] = (1+diff) * acc[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] - - diff = [(extrap_eff[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - eff[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / eff[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] - diff = np.mean(diff) - extrap_eff['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * eff[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - - diff = [(extrap_outinratio[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - outinratio[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / outinratio[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] - diff = np.mean(diff) - extrap_outinratio['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * outinratio[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - - diff = [(extrap_inc_wrongfrac[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - inc_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / inc_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] - diff = np.mean(diff) - extrap_inc_wrongfrac['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * inc_wrongfrac[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - - diff = [(extrap_binfrac_wrongfrac[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - binfrac_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / binfrac_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] if binfrac_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]>0 else 0 for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] - diff = np.mean(diff) - extrap_binfrac_wrongfrac['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * binfrac_wrongfrac[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + # diff = [(extrap_acc[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] - acc[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)]) / acc[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] + # diff = np.mean(diff) + # extrap_acc['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] = (1+diff) * acc[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] + # + # diff = [(extrap_eff[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - eff[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / eff[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] + # diff = np.mean(diff) + # extrap_eff['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * eff[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + # + # diff = [(extrap_outinratio[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - outinratio[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / outinratio[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] + # diff = np.mean(diff) + # extrap_outinratio['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * outinratio[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + # + # diff = [(extrap_inc_wrongfrac[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - inc_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / inc_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] + # diff = np.mean(diff) + # extrap_inc_wrongfrac['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * inc_wrongfrac[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + # + # diff = [(extrap_binfrac_wrongfrac[pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] - binfrac_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]) / binfrac_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] if binfrac_wrongfrac[125][pMode+'125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)]>0 else 0 for pMode in ['ggH', 'VBFH', 'ZH', 'WH']] + # diff = np.mean(diff) + # extrap_binfrac_wrongfrac['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = (1+diff) * binfrac_wrongfrac[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + + extrap_acc['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] = acc[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(genBin)] + + extrap_eff['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = eff[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + + extrap_outinratio['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = outinratio[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + + extrap_inc_wrongfrac['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = inc_wrongfrac[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] + + extrap_binfrac_wrongfrac['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] = binfrac_wrongfrac[125]['ttH125_'+channel+'_'+obsName+'_genbin'+str(genBin)+'_recobin'+str(recoBin)] # if doubleDiff: obs_name_dic = obs_name+'_'+obs_name_2nd # else: obs_name_dic = obs_name @@ -247,6 +257,7 @@ def parseOptions(): with open('../inputs/inputs_sig_extrap_'+obsName+'_'+str(opt.YEAR)+'.py', 'w') as f: f.write('observableBins = '+str(observableBins)+' \n') f.write('acc = '+str(extrap_acc)+' \n') + # f.write('err_acc = '+str(err_acc[125])+' \n') f.write('eff = '+str(extrap_eff)+' \n') f.write('err_eff = '+str(err_eff[125])+' \n') f.write('outinratio = '+str(extrap_outinratio)+' \n') diff --git a/coefficients/RunPlotCoefficients.py b/coefficients/RunPlotCoefficients.py index d424004..69c421b 100644 --- a/coefficients/RunPlotCoefficients.py +++ b/coefficients/RunPlotCoefficients.py @@ -114,6 +114,7 @@ def matrix(obs_bins, obs_name, label): plt.text(j+0.5, i+0.5,'\n\n$\pm${:.3f}'.format(eps_err[i,j]), ha='center', va='center_baseline', size = 'x-small') elif doubleDiff: plt.text(j+0.5, i+0.5, '{:.2f}'.format(z), ha='center', va='center', size = 'x-small') + # plt.text(j+0.5, i+0.5,'\n\n$\pm${:.4f}'.format(eps_err[i,j]), ha='center', va='center_baseline', size = 'x-small') plt.clim(vmin=0.01, vmax=1) # Range colorbar plt.colorbar(label = 'EFFICIENCY (> 0.01)') plt.xticks(obs_bins_label_medium, tickLabel) diff --git a/coefficients/pdfUncROOT.py b/coefficients/pdfUncROOT.py index 44ff1b4..9522ef2 100644 --- a/coefficients/pdfUncROOT.py +++ b/coefficients/pdfUncROOT.py @@ -133,7 +133,7 @@ def getunc(channel, List, m4l_bins, m4l_low, m4l_high, obs_reco, obs_gen, obs_bi cutnoth4l_gen = "(!"+cuth4l_gen+")" shortname = sample_shortnames[Sample] - processBin = shortname+'_'+channel+'_'+opt.OBSNAME+'_genbin'+str(genbin) + processBin = shortname+'_'+channel+'_'+output_name+'_genbin'+str(genbin) # GEN level Histos[processBin+"fs"] = TH1D(processBin+"fs", processBin+"fs", 100, -1, 10000) @@ -243,11 +243,15 @@ def getunc(channel, List, m4l_bins, m4l_low, m4l_high, obs_reco, obs_gen, obs_bi obs_bins, doubleDiff = binning(opt.OBSNAME) if doubleDiff: obs_name = opt.OBSNAME.split(' vs ')[0] - obs_name_2nd = opt.OBSNAME.split(' vs ')[1] + obs_name_2nd_out = opt.OBSNAME.split(' vs ')[1] obs_name_2d = opt.OBSNAME else: obs_name = opt.OBSNAME +output_name = obs_name +if doubleDiff: output_name += '_'+obs_name_2nd_out +print output_name + _temp = __import__('observables', globals(), locals(), ['observables'], -1) observables = _temp.observables print(observables) @@ -260,7 +264,6 @@ def getunc(channel, List, m4l_bins, m4l_low, m4l_high, obs_reco, obs_gen, obs_bi obs_reco = observables[obs_name]['obs_reco'] obs_gen = observables[obs_name]['obs_gen'] - Tree = {} sample_shortnames = {} if opt.NNLOPS: @@ -288,25 +291,12 @@ def getunc(channel, List, m4l_bins, m4l_low, m4l_high, obs_reco, obs_gen, obs_bi else: getunc(chan, List, m4l_bins, m4l_low, m4l_high, obs_reco, obs_gen, obs_bins, genbin) -if (obs_reco.startswith("njets")): - for chan in chans: - for genbin in range(len(obs_bins)-2): # last bin is >=3 - for Sample in List: - shortname = sample_shortnames[Sample] - processBin = shortname+'_'+chan+'_'+obs_reco+'_genbin'+str(genbin) - processBinPlus1 = shortname+'_'+chan+'_'+obs_reco+'_genbin'+str(genbin+1) - acceptance[processBin] = acceptance[processBin]-acceptance[processBinPlus1] - qcdUncert[processBin]['uncerUp'] = sqrt(qcdUncert[processBin]['uncerUp']*qcdUncert[processBin]['uncerUp']+qcdUncert[processBinPlus1]['uncerUp']*qcdUncert[processBinPlus1]['uncerUp']) - qcdUncert[processBin]['uncerDn'] = sqrt(qcdUncert[processBin]['uncerDn']*qcdUncert[processBin]['uncerDn']+qcdUncert[processBinPlus1]['uncerDn']*qcdUncert[processBinPlus1]['uncerDn']) - if opt.NNLOPS: nnlops = "_NNLOPS" else: nnlops = "" -output_name = obs_name -if doubleDiff: output_name += '_'+obs_name_2d -with open('../inputs/accUnc_'+opt.OBSNAME+nnlops+'.py', 'w') as f: +with open('../inputs/accUnc_'+output_name+nnlops+'.py', 'w') as f: #f.write('acc = '+str(acceptance)+' \n') f.write('qcdUncert = '+str(qcdUncert)+' \n') f.write('pdfUncert = '+str(pdfUncert)+' \n') diff --git a/condor_submission/batchScript.sh b/condor_submission/batchScript.sh index a7a546f..71affeb 100644 --- a/condor_submission/batchScript.sh +++ b/condor_submission/batchScript.sh @@ -118,7 +118,9 @@ if [[ "VAR" != "pT4l_kL" ]]; then cp /home/llr/cms/tarabini/CMSSW_10_2_13/src/HiggsAnalysis/FiducialXSFWK/producePlots.py . cp /home/llr/cms/tarabini/CMSSW_10_2_13/src/HiggsAnalysis/FiducialXSFWK/producePlots_v4.py . python producePlots.py --obsName 'VAR' --year 'Full' UNBLIND + python producePlots.py --obsName 'VAR' --year 'Full' UNBLIND --setLog THIRD UNBLIND + THIRD UNBLIND --setLog fi ## ----- Moving all outputs ----- ## diff --git a/condor_submission/sub.sh b/condor_submission/sub.sh index 7392322..f441f5e 100644 --- a/condor_submission/sub.sh +++ b/condor_submission/sub.sh @@ -5,16 +5,16 @@ declare -a obs=( "mass4l noJES" "mass4l_zzfloating noJES" # "njets_pt30_eta4p7 JES" -"pT4l noJES" -"pT4l_kL noJES" -"rapidity4l noJES" -"costhetaZ1 noJES" -"costhetaZ2 noJES" -"phi noJES" -"phistar noJES" -"costhetastar noJES" -"massZ1 noJES" -"massZ2 noJES" +# "pT4l noJES" +# "pT4l_kL noJES" +# "rapidity4l noJES" +# "costhetaZ1 noJES" +# "costhetaZ2 noJES" +# "phi noJES" +# "phistar noJES" +# "costhetastar noJES" +# "massZ1 noJES" +# "massZ2 noJES" # "pTj1 JES" # "pTHj JES" # "mHj JES" @@ -22,7 +22,7 @@ declare -a obs=( # "mjj JES" # "absdetajj JES" # "dphijj JES" -# "pTHjj JES" +"pTHjj JES" # "TCjmax JES" # "TBjmax JES" "D0m noJES" @@ -113,8 +113,8 @@ for i in "${!obs[@]}"; do LLscanv4="python plot_LLScan.py --obsName 'DL1' --year 'Full' --v4" impactsv4="python impacts.py --obsName 'DL1' --year 'Full' --physicsModel 'v4'" elif [ $name_folder == DL1Zg ]; then - CoeffAC="python RunCoefficients.py --obsName 'DL1Zg' --year '2016' --AC_onlyAcc --AC_hypothesis '0L1'" - CoeffACbis="python RunCoefficients.py --obsName 'DL1Zg' --year '2016' --AC_onlyAcc --AC_hypothesis '0L1f05ph0'" + CoeffAC="python RunCoefficients.py --obsName 'DL1Zg' --year '2017' --AC_onlyAcc --AC_hypothesis '0L1Zg'" + CoeffACbis="python RunCoefficients.py --obsName 'DL1Zg' --year '2017' --AC_onlyAcc --AC_hypothesis '0L1Zgf05ph0'" plotsv4="python producePlots_v4.py --obsName 'DL1Zg' --year 'Full'" LLscanv4="python plot_LLScan.py --obsName 'DL1Zg' --year 'Full' --v4" impactsv4="python impacts.py --obsName 'DL1Zg' --year 'Full' --physicsModel 'v4'" @@ -180,8 +180,10 @@ for i in "${!obs[@]}"; do sed -i "s/FIFTH/$impactsv4/g" batchScript.sh - if [ $name_folder == mass4l_zzfloating ]; then - zzfloating="--m4lLower 105 --m4lUpper 160" + # if [ $name_folder == mass4l_zzfloating ]; then + # zzfloating="--m4lLower 105 --m4lUpper 160" + if [ $name_folder == pT4l_kL ]; then + zzfloating="--m4lLower 105 --m4lUpper 140" else zzfloating="--m4lLower 105 --m4lUpper 160" #All measurements between [105,160] fi diff --git a/copy_to_www.sh b/copy_to_www.sh index 20bb821..4b185be 100644 --- a/copy_to_www.sh +++ b/copy_to_www.sh @@ -52,20 +52,44 @@ cp coefficients/matrix_nonfid/2017/nonFid_2017_${1}_* $path/${1}/. cp coefficients/matrix_nonfid/2018/nonFid_2018_${1}_* $path/${1}/. # Move datacards -cp datacard/datacard_2016/hzz4l_*_13TeV_xs_${1}_bin*_v* $path/${1}/datacard_2016/. -cp datacard/datacard_2017/hzz4l_*_13TeV_xs_${1}_bin*_v* $path/${1}/datacard_2017/. -cp datacard/datacard_2018/hzz4l_*_13TeV_xs_${1}_bin*_v* $path/${1}/datacard_2018/. -cp datacard/hzz4l_*_13TeV_xs_${1}_bin*_v* $path/${1}/. +if [ ${1} = "njets_pt30_eta4p7" ]; then + obs_datacards="NJ" +elif [ ${1} = "pT4l" ]; then + obs_datacards="PTH" +elif [ ${1} = "rapidity4l" ]; then + obs_datacards="YH" +elif [ ${1} = "pTj1" ]; then + obs_datacards="PTJET" +else + obs_datacards=${1} +fi + +cp datacard/datacard_2016/hzz4l_*_13TeV_xs_${obs_datacards}_bin*_* $path/${1}/datacard_2016/. +cp datacard/datacard_2017/hzz4l_*_13TeV_xs_${obs_datacards}_bin*_* $path/${1}/datacard_2017/. +cp datacard/datacard_2018/hzz4l_*_13TeV_xs_${obs_datacards}_bin*_* $path/${1}/datacard_2018/. # Move ws -cp datacard/datacard_2016/hzz4l_*_13TeV_xs_SM_125_${1}_v* $path/${1}/datacard_2016/. -cp datacard/datacard_2017/hzz4l_*_13TeV_xs_SM_125_${1}_v* $path/${1}/datacard_2017/. -cp datacard/datacard_2018/hzz4l_*_13TeV_xs_SM_125_${1}_v* $path/${1}/datacard_2018/. +if [ ${1} = "pT4l_kL" ]; then + cp datacard/datacard_2016/hzz4l_*_13TeV_xs_SM_125_${obs_datacards}_kLambda* $path/${1}/datacard_2016/. + cp datacard/datacard_2017/hzz4l_*_13TeV_xs_SM_125_${obs_datacards}_kLambda* $path/${1}/datacard_2017/. + cp datacard/datacard_2018/hzz4l_*_13TeV_xs_SM_125_${obs_datacards}_kLambda* $path/${1}/datacard_2018/. + cp datacard/hzz4l_all_13TeV_xs_${1}_bin_kLambda* $path/${1}/. +else + cp datacard/datacard_2016/hzz4l_*_13TeV_xs_SM_125_${obs_datacards}_v* $path/${1}/datacard_2016/. + cp datacard/datacard_2017/hzz4l_*_13TeV_xs_SM_125_${obs_datacards}_v* $path/${1}/datacard_2017/. + cp datacard/datacard_2018/hzz4l_*_13TeV_xs_SM_125_${obs_datacards}_v* $path/${1}/datacard_2018/. + cp datacard/hzz4l_all_13TeV_xs_${1}_bin_v* $path/${1}/. +fi cp plots/${1}/asimov/${1}_unfoldwith* $path/${1}/. cp plots/${1}/asimov/corr_${1}_*.png $path/${1}/. cp plots/${1}/asimov/model/* $path/${1}/. -cp impacts/${1}/impacts_*_${1}_*_asimov* $path/${1}/. -cp LHScans/plots/lhscan_compare_${1}_* $path/${1}/. +if [ ${1} = "pT4l_kL" ]; then + cp impacts/${1}/impacts_v3_kappa_lambda* $path/${1}/. + cp LHScans/plots/lhscan_compare_${1}_kappa* $path/${1}/. +else + cp impacts/${1}/impacts_*_${1}_*_asimov* $path/${1}/. + cp LHScans/plots/lhscan_compare_${1}_r* $path/${1}/. +fi cp fit/commands_${1}.py $path/${1}/. -cp impacts/commands_impacts_${1}_v* $path/${1}/. +cp impacts/${1}/commands_impacts_${1}_v* $path/${1}/. diff --git a/fit/RunFiducialXS.py b/fit/RunFiducialXS.py index 4011c48..06cc4ed 100644 --- a/fit/RunFiducialXS.py +++ b/fit/RunFiducialXS.py @@ -131,7 +131,7 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', 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: + 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') @@ -144,7 +144,7 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', if int(observableBins[i+1]) > 1000: boundaries = 'GT'+str(int(observableBins[i])) - dc_name = 'datacard_%s/hzz4l_%sS_13TeV_xs_%s_bin%d_v3.txt ' %(year,cat,obsName,i) + dc_name = 'datacard_%s/hzz4l_%sS_13TeV_xs_%s_bin%d_v3.txt ' %(year,cat,fitName,i) cmd_combCards += 'hzz_%s_%s_cat%s_%s=%s' %(fitName,boundaries,cat,year,dc_name) cmd_combCards += '> %s' %card_name @@ -167,19 +167,20 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', cmd_addNuis += '" >> hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.txt' processCmd(cmd_combCards) + cmds.append(cmd_combCards) processCmd(cmd_addNuis) + cmds.append(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: + 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') @@ -208,11 +209,27 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', 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=100 --saveToys --cminDefaultMinimizerStrategy 0 -t -1 --setParameters ' - cmd_fit_tmp = cmd_fit + '%s=1 -P %s --setParameterRanges %s=0.0,2.0 --redefineSignalPOI %s' %(POI, POI, POI, POI) + 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) + cmd_fit_tmp = cmd_fit + '-P %s --setParameterRanges %s=0.0,2.5 --redefineSignalPOI %s' %(POI, POI, POI) print(cmd_fit_tmp) processCmd(cmd_fit_tmp) + cmds.append(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) + cmds.append(cmd_fit_tmp) if obsName == 'mass4l_zzfloating': for i in range(nBins): @@ -231,10 +248,26 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', 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, 'SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root') - cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=100 --saveToys --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 -t -1 --setParameters ' - cmd_fit_tmp = cmd_fit + '%s=1 -P %s --setParameterRanges %s=0.0,2.0 --redefineSignalPOI %s' %(POI, POI, POI, POI) + cmd_fit += '-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) + cmd_fit_tmp = cmd_fit + '-P %s --setParameterRanges %s=0.0,2.5 --redefineSignalPOI %s' %(POI, POI, POI) processCmd(cmd_fit_tmp) + cmds.append(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_NoSys -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 --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) + cmds.append(cmd_fit_tmp) if obsName == 'mass4l_zzfloating': for i in range(nBins): @@ -266,7 +299,7 @@ def runFiducialXS(): print 'Theory xsec and BR at MH = '+_th_MH print 'Current directory: python' - _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta2p5': 'NJ'} + _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} if obsName not in _obsName: _obsName[obsName] = obsName @@ -414,7 +447,7 @@ def runFiducialXS(): # nBins = len(observableBins) if physicalModel == 'v2': # In this case implemented for mass4l only for channel in ['4e', '4mu', '2e2mu']: - cmd = 'combine -n _'+obsName+'_r'+channel+'Bin0 -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v2.root -m 125.38 --freezeParameters MH -P r'+channel+'Bin0 --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r'+channel+'Bin0=0.0,2.5 --redefineSignalPOI r'+channel+'Bin0 --algo=grid --points=100 --cminDefaultMinimizerStrategy 0 --saveInactivePOI=1' + cmd = 'combine -n _'+obsName+'_r'+channel+'Bin0 -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v2.root -m 125.38 --freezeParameters MH -P r'+channel+'Bin0 --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r'+channel+'Bin0=0.0,2.5 --redefineSignalPOI r'+channel+'Bin0 --algo=grid --points=200 --cminDefaultMinimizerStrategy 0 --saveInactivePOI=1' fidxs = 0 fidxs = higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ggH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] @@ -431,7 +464,7 @@ def runFiducialXS(): # if(not opt.UNBLIND): cmd = cmd + '_exp' cmd = cmd + ' -M MultiDimFit higgsCombine_'+obsName+'_r'+channel+'Bin0.MultiDimFit.mH125.38' if(not opt.UNBLIND): cmd = cmd + '.123456' - cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P r'+channel+'Bin0 --floatOtherPOIs=1 --saveWorkspace --setParameterRanges SigmaBin0=0.0,2.5 --redefineSignalPOI r'+channel+'Bin0 --algo=grid --points=100 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' + cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P r'+channel+'Bin0 --floatOtherPOIs=1 --saveWorkspace --setParameterRanges SigmaBin0=0.0,2.5 --redefineSignalPOI r'+channel+'Bin0 --algo=grid --points=200 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' if (opt.YEAR == 'Full'): cmd = cmd + ' --freezeParameters MH' else: cmd = cmd + ' --freezeParameters MH' if(not opt.UNBLIND): cmd = cmd + ' -t -1 --saveToys --setParameters r'+channel+'Bin0='+str(round(fidxs,4)) @@ -442,7 +475,7 @@ def runFiducialXS(): if physicalModel == 'v4': for obsBin in range(nBins): # ----- 2e2mu ----- - cmd = 'combine -n _'+obsName+'_r2e2muBin'+str(obsBin)+' -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v4.root -m 125.38 --freezeParameters MH -P r2e2muBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r2e2muBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r2e2muBin'+str(obsBin)+' --algo=grid --points=100 --cminDefaultMinimizerStrategy 0 --saveInactivePOI=1' + cmd = 'combine -n _'+obsName+'_r2e2muBin'+str(obsBin)+' -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v4.root -m 125.38 --freezeParameters MH -P r2e2muBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r2e2muBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r2e2muBin'+str(obsBin)+' --algo=grid --points=200 --cminDefaultMinimizerStrategy 0 --saveInactivePOI=1' fidxs = 0 fidxs = higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_2e2mu']*acc['ggH125_2e2mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] @@ -459,7 +492,7 @@ def runFiducialXS(): # if(not opt.UNBLIND): cmd = cmd + '_exp' cmd = cmd + ' -M MultiDimFit higgsCombine_'+obsName+'_r2e2muBin'+str(obsBin)+'.MultiDimFit.mH125.38' if(not opt.UNBLIND): cmd = cmd + '.123456' - cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P r2e2muBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r2e2muBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r2e2muBin'+str(obsBin)+' --algo=grid --points=100 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' + cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P r2e2muBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r2e2muBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r2e2muBin'+str(obsBin)+' --algo=grid --points=200 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' if (opt.YEAR == 'Full'): cmd = cmd + ' --freezeParameters MH' else: cmd = cmd + ' --freezeParameters MH' if(not opt.UNBLIND): cmd = cmd + ' -t -1 --saveToys --setParameters r2e2muBin'+str(obsBin)+'='+str(round(fidxs,4)) @@ -468,7 +501,7 @@ def runFiducialXS(): cmds.append(cmd) # ----- 4e+4mu = 4l ----- - cmd = 'combine -n _'+obsName+'_r4lBin'+str(obsBin)+' -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v4.root -m 125.38 --freezeParameters MH -P r4lBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r4lBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r4lBin'+str(obsBin)+' --algo=grid --points=100 --cminDefaultMinimizerStrategy 0' + cmd = 'combine -n _'+obsName+'_r4lBin'+str(obsBin)+' -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v4.root -m 125.38 --freezeParameters MH -P r4lBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r4lBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r4lBin'+str(obsBin)+' --algo=grid --points=200 --cminDefaultMinimizerStrategy 0' fidxs = 0 # 4e @@ -492,7 +525,7 @@ def runFiducialXS(): # if(not opt.UNBLIND): cmd = cmd + '_exp' cmd = cmd + ' -M MultiDimFit higgsCombine_'+obsName+'_r4lBin'+str(obsBin)+'.MultiDimFit.mH125.38' if(not opt.UNBLIND): cmd = cmd + '.123456' - cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P r4lBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r4lBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r4lBin'+str(obsBin)+' --algo=grid --points=100 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' + cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P r4lBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r4lBin'+str(obsBin)+'=0.0,2.5 --redefineSignalPOI r4lBin'+str(obsBin)+' --algo=grid --points=200 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' if (opt.YEAR == 'Full'): cmd = cmd + ' --freezeParameters MH' else: cmd = cmd + ' --freezeParameters MH' if(not opt.UNBLIND): cmd = cmd + ' -t -1 --saveToys --setParameters r4lBin'+str(obsBin)+'='+str(round(fidxs,4)) @@ -562,7 +595,7 @@ def runFiducialXS(): ## The inclusive xsec for 2j phase space is about 2.49 fb, hence enlarge fit range if ('jj' in obsName) and (obsBin == 0): max_range = '5.0' if ('njet' in obsName) and ('pTj' in obsName) and (obsBin == 0): max_range = '5.0' - cmd = 'combine -n _'+obsName+'_SigmaBin'+str(obsBin)+' -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v3.root -m 125.38 --freezeParameters MH -P SigmaBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges SigmaBin'+str(obsBin)+'=0.0,'+max_range+' --redefineSignalPOI SigmaBin'+str(obsBin)+' --algo=grid --points=100 --cminDefaultMinimizerStrategy 0' + cmd = 'combine -n _'+obsName+'_SigmaBin'+str(obsBin)+' -M MultiDimFit SM_125_all_13TeV_xs_'+obsName+'_bin_v3.root -m 125.38 --freezeParameters MH -P SigmaBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges SigmaBin'+str(obsBin)+'=0.0,'+max_range+' --redefineSignalPOI SigmaBin'+str(obsBin)+' --algo=grid --points=200 --cminDefaultMinimizerStrategy 0' if(not opt.UNBLIND): cmd = cmd + ' -t -1 --saveToys --setParameters SigmaBin'+str(obsBin)+'='+str(round(_obsxsec,4)) if opt.FIXFRAC: @@ -589,7 +622,7 @@ def runFiducialXS(): # if(not opt.UNBLIND): cmd = cmd + '_exp' cmd = cmd + ' -M MultiDimFit higgsCombine_'+obsName+'_SigmaBin'+str(obsBin)+'.MultiDimFit.mH125.38' if(not opt.UNBLIND): cmd = cmd + '.123456' - cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P SigmaBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges SigmaBin0=0.0,'+max_range+' --redefineSignalPOI SigmaBin'+str(obsBin)+' --algo=grid --points=100 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' + cmd = cmd + '.root -w w --snapshotName "MultiDimFit" -m 125.38 -P SigmaBin'+str(obsBin)+' --floatOtherPOIs=1 --saveWorkspace --setParameterRanges SigmaBin0=0.0,'+max_range+' --redefineSignalPOI SigmaBin'+str(obsBin)+' --algo=grid --points=200 --cminDefaultMinimizerStrategy 0 --freezeNuisanceGroups nuis' if (opt.YEAR == 'Full'): cmd = cmd + ' --freezeParameters MH' else: cmd = cmd + ' --freezeParameters MH' if(not opt.UNBLIND): @@ -608,11 +641,13 @@ def runFiducialXS(): cmd = 'combine SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root -n _'+obsName+' -M MultiDimFit --algo=singles -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --freezeParameters MH,r --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1' output = processCmd(cmd) cmds.append(cmd) + print(cmd) #Stat+sys grid - cmd = 'combine SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root -n _'+obsName+'_grid -M MultiDimFit --algo=grid --points=300 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --freezeParameters MH,r --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1' + cmd = 'combine SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root -n _'+obsName+'_grid -M MultiDimFit --algo=grid --points=250 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --freezeParameters MH,r --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1' output = processCmd(cmd) cmds.append(cmd) + print(cmd) #Stat-only singles cmd = 'combine higgsCombine_'+obsName+'.MultiDimFit.mH125.38' @@ -622,16 +657,17 @@ def runFiducialXS(): else: cmd += ' --freezeParameters MH' output = processCmd(cmd) cmds.append(cmd) + print(cmd) #Stat-only grid cmd = 'combine higgsCombine_'+obsName+'_grid.MultiDimFit.mH125.38' if(not opt.UNBLIND): cmd += '.123456' - cmd += '.root -n _'+obsName+'_NoSys_grid -M MultiDimFit -w w --snapshotName "MultiDimFit" --algo=grid --points=300 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1 --freezeNuisanceGroups nuis' + cmd += '.root -n _'+obsName+'_NoSys_grid -M MultiDimFit -w w --snapshotName "MultiDimFit" --algo=grid --points=250 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1 --freezeNuisanceGroups nuis' if (opt.YEAR == 'Full'): cmd += ' --freezeParameters MH,r' else: cmd += ' --freezeParameters MH' output = processCmd(cmd) cmds.append(cmd) - + print(cmd) # ----------------- Main ----------------- _fit_dir = os.getcwd() diff --git a/fit/addConstrainedModel.py b/fit/addConstrainedModel.py index dbc3577..1489c09 100644 --- a/fit/addConstrainedModel.py +++ b/fit/addConstrainedModel.py @@ -101,6 +101,7 @@ def parseOptions(): for genbin in range(lenObsBins): + print 'ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin) ggHxs = higgs_xs['ggH_125.0']*higgs4l_br['125.0_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] #ggHxs = acc['ggH_HRes_125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] VBFxs = higgs_xs['VBF_125.0']*higgs4l_br['125.0_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] @@ -108,7 +109,6 @@ def parseOptions(): ZHxs = higgs_xs['ZH_125.0']*higgs4l_br['125.0_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] ttHxs = higgs_xs['ttH_125.0']*higgs4l_br['125.0_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - effsm = ggHxs/(ggHxs+VBFxs+WHxs+ZHxs+ttHxs)*max(eff['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(recobin)],0.0) #effsm = ggHxs/(ggHxs+VBFxs+WHxs+ZHxs+ttHxs)*max(eff['ggH_HRes_125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(recobin)],0.0) effsm += VBFxs/(ggHxs+VBFxs+WHxs+ZHxs+ttHxs)*max(eff['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(recobin)],0.0) diff --git a/fit/createDatacard.py b/fit/createDatacard.py index 45e9178..9b5327b 100644 --- a/fit/createDatacard.py +++ b/fit/createDatacard.py @@ -1,7 +1,8 @@ import os,sys + def fixJes(jesnp): - # Cases where jes are '-' or single value + # Cases where jes are '-' or single value if (len(jesnp)==1) | ('/' not in jesnp): return jesnp else: @@ -11,12 +12,16 @@ def fixJes(jesnp): if ((jesnp_tmp_dn=='1.0') & (jesnp_tmp_up!='1.0')): return jesnp_tmp_up elif ((jesnp_tmp_dn!='1.0') & (jesnp_tmp_up=='1.0')): - return jesnp_tmp_up + return jesnp_tmp_dn+'/'+ str(abs(round(2-float(jesnp_tmp_dn),3))) + elif ((float(jesnp_tmp_dn) > 1) & (float(jesnp_tmp_up) > 1)): + return jesnp_tmp_up + '/' + str(abs(round(2-float(jesnp_tmp_up),3))) + elif ((float(jesnp_tmp_dn) < 1) & (float(jesnp_tmp_up) < 1)): + return jesnp_tmp_dn + '/' + str(abs(round(2-float(jesnp_tmp_dn),3))) elif (float(jesnp_tmp_dn) > float(jesnp_tmp_up)): ''' - if 1.X/0.X cases - return a single NP (symmetric) corresponding to largest variation - return the variation always as 1.X + if 1.X/0.X cases + return a single NP (symmetric) corresponding to largest variation + return the variation always as 1.X ''' if((float(jesnp_tmp_dn)-1) > (1-float(jesnp_tmp_up))): return jesnp_tmp_dn @@ -24,11 +29,10 @@ def fixJes(jesnp): return str(1+(1-float(jesnp_tmp_up))) else: ''' - if dn/up: correct jes, use it + if dn/up: correct jes, use it ''' return jesnp - def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalModel, year, nData, jes, lowerBound, upperBound, yearSetting): # Name of the bin (aFINALSTATE_ recobinX) if(channel == '4mu'): channelNumber = 1 @@ -45,7 +49,7 @@ def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalMode # if int(observableBins[obsBin+1]) > 1000: # _recobin = 'GT'+str(int(observableBins[obsBin])) - if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName: #it means it is a double differential measurement + if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName and not obsName == 'njets_pt30_eta4p7': #it means it is a double differential measurement _recobin = str(observableBins[obsBin][0]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[obsBin][1]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[obsBin][2]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[obsBin][3]).replace('.', 'p').replace('-','m') else: _recobin = str(observableBins[obsBin]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[obsBin+1]).replace('.', 'p').replace('-','m') @@ -53,12 +57,14 @@ def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalMode _recobin = 'GT'+str(int(observableBins[obsBin])) if physicalModel == 'v3': - _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta2p5': 'NJ'} + _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} if obsName not in _obsName: _obsName[obsName] = obsName binName = 'hzz_' + _obsName[obsName] + '_' + _recobin + '_cat' + channel processName = 'smH_' + _obsName[obsName] else: + _obsName = {} + _obsName[obsName] = obsName binName = 'a'+str(channelNumber)+'_recobin'+str(obsBin) processName = 'trueH'+channel+'Bin' @@ -138,7 +144,7 @@ def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalMode # ------------------------------------------------------------------------------------------------- - file = open('../datacard/datacard_'+year+'/hzz4l_'+channel+'S_13TeV_xs_'+obsName+'_bin'+str(obsBin)+'_'+physicalModel+'.txt', 'w+') + file = open('../datacard/datacard_'+year+'/hzz4l_'+channel+'S_13TeV_xs_'+_obsName[obsName]+'_bin'+str(obsBin)+'_'+physicalModel+'.txt', 'w+') file.write('imax 1 \n') file.write('jmax * \n') @@ -146,7 +152,7 @@ def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalMode file.write('------------ \n') - file.write('shapes * * hzz4l_'+channel+'S_13TeV_xs_SM_125_'+obsName+'_'+physicalModel+'.Databin'+str(obsBin)+'.root w:$PROCESS\n') + file.write('shapes * * hzz4l_'+channel+'S_13TeV_xs_SM_125_'+_obsName[obsName]+'_'+physicalModel+'.Databin'+str(obsBin)+'.root w:$PROCESS\n') file.write('------------ \n') @@ -163,7 +169,7 @@ def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalMode file.write('process ') if physicalModel == 'v3': for i in range(nBins): - if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName: + if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName and not obsName == 'njets_pt30_eta4p7': file.write(processName+'_'+str(observableBins[i][0]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[i][1]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[i][2]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[i][3]).replace('.', 'p').replace('-','m')+' ') elif observableBins[i+1] > 1000: file.write(processName+'_GT'+str(int(observableBins[i]))+' ') @@ -310,28 +316,18 @@ def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalMode if jes == True: for index,jesName in enumerate(jesNames_datacard): - if obsName == 'pTj1': - if obsBin == 0: - print jesnp['signal_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)] - p = float(jesnp['signal_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)][:-4]) - # p = p-1 - print 2-p - - - sys.out() - file.write('CMS_scale_j_'+jesName+' lnN ') - for i in range(nBins): # Signals - file.write(fixJes(str(jesnp['fiducial_'+jesNames[index]+'_'+channel+'_'+obsName+'_genbin'+str(i)+'_recobin'+str(obsBin)]))+' ') - file.write(fixJes(str(jesnp['nonFiducial_'+jesNames[index]+'_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]))+' ') - file.write(fixJes(str(jesnp['nonResonant_'+jesNames[index]+'_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]))+' ') - file.write(fixJes(str(jesnp['qqzz_'+jesNames[index]+'_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]))+' ') - file.write(fixJes(str(jesnp['ggzz_'+jesNames[index]+'_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]))+' ') - file.write('-\n') # ZX - file.write('CMS_scale_j_ZX lnN ') - for i in range(nBins+4): # All except ZX - file.write('- ') - file.write(fixJes(str(jesnp['ZX_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)]))+'\n') + for i in range(nBins+2): # Signals + out + fake + file.write(str(fixJes(jesnp['signal_'+jesNames[index]+'_'+channel+'_'+year+'_'+obsName.replace('pT4l', 'ZZPt')+'_recobin'+str(obsBin)]))+' ') + file.write(str(fixJes(jesnp['qqzz_'+jesNames[index]+'_'+channel+'_'+year+'_'+obsName.replace('pT4l', 'ZZPt')+'_recobin'+str(obsBin)]))+' ') + file.write(str(fixJes(jesnp['ggzz_'+jesNames[index]+'_'+channel+'_'+year+'_'+obsName.replace('pT4l', 'ZZPt')+'_recobin'+str(obsBin)]))+' ') + file.write(str(fixJes(jesnp['ZX_'+jesNames[index]+'_'+channel+'_'+year+'_'+obsName.replace('pT4l', 'ZZPt')+'_recobin'+str(obsBin)]))+'\n') + # file.write('CMS_scale_j_ZX lnN ') + # for i in range(nBins+4): # All except ZX + # file.write('- ') + # file.write(str(jesnp['ZX_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)])+'\n') + + #file.write('JES param 0.0 1.0\n') file.close() @@ -428,7 +424,7 @@ def createDatacard_ggH(obsName, channel, nBins, obsBin, observableBins, physical # ------------------------------------------------------------------------------------------------- - file = open('../datacard/datacard_'+year+'/hzz4l_GGH_'+channel+'S_13TeV_xs_'+obsName+'_bin'+str(obsBin)+'_'+physicalModel+'.txt', 'w+') + file = open('../datacard/datacard_'+year+'/hzz4l_GGH_'+channel+'S_13TeV_xs_'+_obsName[obsName]+'_bin'+str(obsBin)+'_'+physicalModel+'.txt', 'w+') file.write('imax 1 \n') file.write('jmax * \n') @@ -436,7 +432,7 @@ def createDatacard_ggH(obsName, channel, nBins, obsBin, observableBins, physical file.write('------------ \n') - file.write('shapes * * hzz4l_'+channel+'S_13TeV_xs_SM_125_'+obsName+'_'+physicalModel+'.Databin'+str(obsBin)+'.root w:$PROCESS\n') + file.write('shapes * * hzz4l_'+channel+'S_13TeV_xs_SM_125_'+_obsName[obsName]+'_'+physicalModel+'.Databin'+str(obsBin)+'.root w:$PROCESS\n') file.write('------------ \n') @@ -451,7 +447,7 @@ def createDatacard_ggH(obsName, channel, nBins, obsBin, observableBins, physical file.write(binName+' ') file.write('\n') file.write('process ') - if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName: + if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName and not obsName == 'njets_pt30_eta4p7': for i in range(nBins): file.write(processName+'_'+str(observableBins[i][0]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[i][1]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[i][2]).replace('.', 'p').replace('-','m')+'_'+str(observableBins[i][3]).replace('.', 'p').replace('-','m')+' ') for i in range(nBins): @@ -613,14 +609,14 @@ def createDatacard_ggH(obsName, channel, nBins, obsBin, observableBins, physical file.write(fixJes(str(jesnp['signal_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)]))+' ') # Bkgs - file.write(fixJes(str(jesnp['qqzz_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)]))+' ') - file.write(fixJes(str(jesnp['ggzz_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)]))+' ') - file.write('-\n') # ZX + file.write(str(jesnp['qqzz_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)])+' ') + file.write(str(jesnp['ggzz_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)])+' ') + file.write(str(jesnp['ZX_'+jesNames[index]+'_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)])+'\n') #file.write('JES param 0.0 1.0\n') - file.write('CMS_scale_j_ZX lnN ') - for i in range(nBins+4): # All except ZX - file.write('- ') - file.write(fixJes(str(jesnp['ZX_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)]))+'\n') + # file.write('CMS_scale_j_ZX lnN ') + # for i in range(nBins+4): # All except ZX + # file.write('- ') + # file.write(str(jesnp['ZX_'+channel+'_'+str(year)+'_'+obsName+'_recobin'+str(obsBin)])+'\n') file.close() diff --git a/fit/createXSworkspace.py b/fit/createXSworkspace.py index 324802a..4a05b67 100644 --- a/fit/createXSworkspace.py +++ b/fit/createXSworkspace.py @@ -42,7 +42,7 @@ 'DL1': True, 'DL1Zg': True, 'rapidity4l_pT4l': True, -'njets_pt30_eta4p7 vs pT4l': False, +'njets_pt30_eta4p7_pT4l': False, 'pTj1_pTj2': False, 'pT4l_pTHj': False, 'massZ1_massZ2': False, @@ -139,7 +139,7 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, if int(observableBins[obsBin+1]) > 1000: _recobin = 'GT'+str(int(observableBins[obsBin])) - _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta2p5': 'NJ'} + _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} if obsName not in _obsName: _obsName[obsName] = obsName @@ -518,11 +518,16 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, fidxs = {} for fState in ['4e','4mu']: fidxs[fState] = 0 - fidxs[fState] += higgs_xs['ggH_125.0']*higgs4l_br['125.0_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['VBF_125.0']*higgs4l_br['125.0_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['WH_125.0']*higgs4l_br['125.0_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['ZH_125.0']*higgs4l_br['125.0_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['ttH_125.0']*higgs4l_br['125.0_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + # fidxs[fState] += higgs_xs['ggH_125.0']*higgs4l_br['125.0_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + # fidxs[fState] += higgs_xs['VBF_125.0']*higgs4l_br['125.0_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + # fidxs[fState] += higgs_xs['WH_125.0']*higgs4l_br['125.0_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + # fidxs[fState] += higgs_xs['ZH_125.0']*higgs4l_br['125.0_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + # fidxs[fState] += higgs_xs['ttH_125.0']*higgs4l_br['125.0_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['ggH_125.38']*higgs4l_br['125.38_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['VBF_125.38']*higgs4l_br['125.38_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['WH_125.38']*higgs4l_br['125.38_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['ZH_125.38']*higgs4l_br['125.38_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['ttH_125.38']*higgs4l_br['125.38_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] fidxs['4l'] = fidxs['4e'] + fidxs['4mu'] fracSM4eBin[str(genbin)] = ROOT.RooRealVar('fracSM4eBin'+str(genbin), 'fracSM4eBin'+str(genbin), fidxs['4e']/fidxs['4l']) fracSM4eBin[str(genbin)].setConstant(True) @@ -962,9 +967,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, getattr(wout,'import')(data_obs.reduce(ROOT.RooArgSet(m))) if (addfakeH): - fout = ROOT.TFile('hzz4l_'+channel+'S_13TeV_xs_'+modelName+'_'+obsName+'_'+physicalModel+'.Databin'+str(obsBin)+'.root','RECREATE') + fout = ROOT.TFile('hzz4l_'+channel+'S_13TeV_xs_'+modelName+'_'+_obsName[obsName]+'_'+physicalModel+'.Databin'+str(obsBin)+'.root','RECREATE') else: - fout = ROOT.TFile('hzz4l_'+channel+'S_13TeV_xs_'+modelName+'_'+obsName+'_'+physicalModel+'.Databin'+str(obsBin)+'.NoFakeH.root','RECREATE') + fout = ROOT.TFile('hzz4l_'+channel+'S_13TeV_xs_'+modelName+'_'+_obsName[obsName]+'_'+physicalModel+'.Databin'+str(obsBin)+'.NoFakeH.root','RECREATE') print "write ws to fout" fout.WriteTObject(wout) diff --git a/fit/expected_xsec.py b/fit/expected_xsec.py index 405bdca..b529d2b 100644 --- a/fit/expected_xsec.py +++ b/fit/expected_xsec.py @@ -76,6 +76,7 @@ def exp_xsec(): XH_fs += higgs_xs['WH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['WH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] XH_fs += higgs_xs['ZH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ZH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] XH_fs += higgs_xs['ttH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ttH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + print 'Bin ', obsBin, '\t SigmaBin', obsBin, channel, ' = ', XH_fs XH[obsBin]+=XH_fs # else: # XH_fs = higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+'4l']*acc['ggH125_4l_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] @@ -87,6 +88,7 @@ def exp_xsec(): _obsxsec = XH[obsBin] + print '\n' print 'Bin ', obsBin, '\t SigmaBin', obsBin, ' = ', _obsxsec xs['SigmaBin'+str(obsBin)] = _obsxsec diff --git a/fit/impacts.py b/fit/impacts.py index 3fb9405..2cb4a78 100644 --- a/fit/impacts.py +++ b/fit/impacts.py @@ -92,7 +92,7 @@ def impactPlots(): if not doubleDiff: nBins = nBins-1 #in case of 1D measurement the number of bins is -1 the length of the list of bin boundaries # if '_' in obsName and obsName!='mass4l_zzfloating': nBins = len(observableBins)+1 - _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta2p5': 'NJ'} + _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} if obsName in _obsName: obsName_poi = _obsName[obsName] else: diff --git a/helperstuff/binning.py b/helperstuff/binning.py index c1e2def..2efea1e 100644 --- a/helperstuff/binning.py +++ b/helperstuff/binning.py @@ -19,7 +19,7 @@ 'mjj': '|-2|0|120|300|1300|', 'absdetajj': '|-100|0|0.7|1.6|3.0|10.0|', 'dphijj': '|-100|-3.14159265359|-1.5707963267948966|0|1.5707963267948966|3.14159265359|', -'pTHjj': '|-2|0|20|40|60|110|1300|', +'pTHjj': '|-2|0|20|60|1300|', 'TCjmax': '|-2|15|20|30|50|80|1000|', 'TBjmax': '|-2|30|70|130|250|400|1000|', 'D0m': '|0.0|0.4|0.5|0.6|0.7|0.8|0.9|1.0|', @@ -31,9 +31,9 @@ 'rapidity4l vs pT4l': '|0|0.5|1.0|2.5| vs |0|40|80|150|1300| / |0|45|120|1300| / |0|45|120|1300|', 'njets_pt30_eta4p7 vs pT4l': '|0|1|2|20| vs |0|15|30|1300| / |0|60|80|120|1300| / |0|100|170|250|1300|', 'pTj1 vs pTj2': '|-2|30| vs |-2|30| / |30|60| vs |30|60| / |60|350| vs |30|60| / |60|350| vs |60|350|', -'pT4l vs pTHj': '|-2|0| vs |-2|0| / |0|85| vs |0|30| / |85|350| vs |0|45| / |0|85| vs |30|350| / |85|350| vs |45|350|', +'pT4l vs pTHj': '|0|1500| vs |-2|0| / |0|85| vs |0|30| / |85|350| vs |0|45| / |0|85| vs |30|350| / |85|350| vs |45|350|', 'massZ1 vs massZ2': '|40|85| vs |12|35| / |40|70| vs |35|65| / |70|120| vs |35|65| / |85|120| vs |30|35| / |85|120| vs |24|30| / |85|120| vs |12|24|', -'TCjmax vs pT4l': '|0|15|25|40|1300| vs |0|15|30|45|70|120|1300| / |0|120|1300| / |0|120|1300| / |0|200|1300|', +'TCjmax vs pT4l': '|-2|15|25|40|1300| vs |0|15|30|45|70|120|1300| / |0|120|1300| / |0|120|1300| / |0|200|1300|', } def binning(var): diff --git a/helperstuff/createdf_jes.py b/helperstuff/createdf_jes.py index a68cec9..ed305fe 100644 --- a/helperstuff/createdf_jes.py +++ b/helperstuff/createdf_jes.py @@ -23,7 +23,7 @@ def weight(df, xsec, gen, lumi, additional = None): #Coefficient to calculate weights for histograms coeff = (lumi * 1000 * xsec) / gen #Reco - weight_reco = (df.overallEventWeight * df.L1prefiringWeight) + weight_reco = (df.overallEventWeight * df.L1prefiringWeight * df.SFcorr) if additional == 'ggH': weight_reco *= df.ggH_NNLOPS_weight elif additional == 'qqzz': @@ -235,7 +235,7 @@ def tb(pt,eta,phi,mass,H): def createDataframe(dataFrame,isBkg,gen,xsec,signal,lumi,obs_reco,obs_reco_2nd='None'): b_sig = ['EventNumber', 'PUWeight', 'genHEPMCweight', 'ZZMass', 'ZZPt','ZZEta', 'ZZPhi', 'Z1Flav', 'Z2Flav', 'JetPt', 'JetMass', 'JetEta', 'JetPhi', - 'overallEventWeight', 'L1prefiringWeight','dataMCWeight', 'trigEffWeight', + 'overallEventWeight', 'L1prefiringWeight','dataMCWeight', 'trigEffWeight', 'SFcorr', 'pTj1', 'Mj1', 'ETAj1', 'PHIj1', 'pTj2', 'Mj2', 'ETAj2', 'PHIj2'] @@ -297,6 +297,9 @@ def createDataframe(dataFrame,isBkg,gen,xsec,signal,lumi,obs_reco,obs_reco_2nd=' if obs_reco == 'absdetajj' or obs_reco_2nd == 'absdetajj': df['absdetajj_jesup_'+i] = [abs(j1.Eta()-j2.Eta()) if j2.Pt()>0 else -1 for j1,j2 in zip(df['j1_jesup_'+i],df['j2_jesup_'+i])] df['absdetajj_jesdn_'+i] = [abs(j1.Eta()-j2.Eta()) if j2.Pt()>0 else -1 for j1,j2 in zip(df['j1_jesdn_'+i],df['j2_jesdn_'+i])] + if obs_reco == 'dphijj' or obs_reco_2nd == 'dphijj': + df['dphijj_jesup_'+i] = [math.atan2(math.sin(j1.Phi()-j2.Phi()), math.cos(j1.Phi()-j2.Phi())) if j2.Pt()>0 else -99 for j1,j2 in zip(df['j1_jesup_'+i],df['j2_jesup_'+i])] + df['dphijj_jesdn_'+i] = [math.atan2(math.sin(j1.Phi()-j2.Phi()), math.cos(j1.Phi()-j2.Phi())) if j2.Pt()>0 else -99 for j1,j2 in zip(df['j1_jesdn_'+i],df['j2_jesdn_'+i])] if obs_reco == 'ZZPt' or obs_reco_2nd == 'ZZPt': df['ZZPt_jesup_'+i] = df['ZZPt'] df['ZZPt_jesdn_'+i] = df['ZZPt'] diff --git a/helperstuff/observables.py b/helperstuff/observables.py index 565723b..a0b8447 100644 --- a/helperstuff/observables.py +++ b/helperstuff/observables.py @@ -115,8 +115,8 @@ 'obs_gen': 'GENmassZ1', 'obs_gen_2nd': 'GENmassZ2'}, 'TCjmax vs pT4l': - {'obs_reco_2nd': 'TCjmax', - 'obs_gen_2nd': 'GENTCjmax', - 'obs_reco': 'ZZPt', - 'obs_gen': 'GENpT4l'}, + {'obs_reco_2nd': 'ZZPt', + 'obs_gen_2nd': 'GENpT4l', + 'obs_reco': 'TCjmax', + 'obs_gen': 'GENTCjmax'}, } diff --git a/inputs/binning.py b/inputs/binning.py deleted file mode 100644 index 61530ba..0000000 --- a/inputs/binning.py +++ /dev/null @@ -1,110 +0,0 @@ -list = {'rapidity4l': '|0.0|0.15|0.3|0.6|0.9|1.2|2.5|',#'|0.0|0.15|0.3|0.45|0.6|0.75|0.9|1.2|1.6|2.0|2.5|', -'pT4l': '|0|10|20|30|45|60|80|120|200|13000|', -'pT4l_kL': '|0.0|45.0|80.0|120.0|200.0|1300.0|', -'massZ1': '|50|64|73|85|106|', -'massZ2': '|12|20|24|28|32|40|55|65|', -'mass4l': '|105|140|', -'njets_pt30_eta2p5': '|0|1|2|3|4|14|', -'pTj1': '|-2|30|55|95|200|1300|', -'njets_pt30_eta2p5 vs pTj1': '|0|1|10| vs |0|1| / |30|60|120|350|', -'njets_pt30_eta4p7 vs pTj1_eta4p7': '|0|1|10| vs |0|1| / |30|60|120|350|', -'njets_pt30_eta2p5 vs pTj2': '|0|2|10| vs |0|1| / |30|60|120|350|', -'costhetastar': '|0.0|0.2|0.4|0.6|0.8|1.0|', #|0.0|0.125|0.25|0.375|0.5|0.625|0.75|0.875|1.0| -'costhetaZ1': '|-1.0|-0.75|-0.50|-0.25|0.0|0.25|0.50|0.75|1.0|', #'|0.0|0.2|0.4|0.6|0.8|1.0|' -'costhetaZ2': '|-1.0|-0.75|-0.50|-0.25|0.0|0.25|0.50|0.75|1.0|', #'|0.0|0.2|0.4|0.6|0.8|1.0|' -'phi': '|-3.14159265359|-2.35619449019|-1.57079632679|-0.785398163397|0.0|0.785398163397|1.57079632679|2.35619449019|3.14159265359|', -'phistar': '|-3.14159265359|-2.35619449019|-1.57079632679|-0.785398163397|0.0|0.785398163397|1.57079632679|2.35619449019|3.14159265359|', -'TCjmax': '|0|15|30|50|70|90|1000|', ## Should treat as 2D? In principle yes, leading jet involved -'TBjmax': '|0|30|60|80|100|120|1000|', ## Should treat as 2D? In principle yes, leading jet involved -'pTHj': '|-2|0|30|70|1000|', -'rapidity4l vs pT4l': '|0.0|0.5| vs |0|40| / |0.0|0.5| vs |40|80| / |0.0|0.5| vs |80|150| / |0.0|0.5| vs |150|1300| / |0.5|1.0| vs |0|45| / |0.5|1.0| vs |45|120| / |0.5|1.0| vs |120|1300| / |1.0|2.5| vs |0|45| / |1.0|2.5| vs |45|120| / |1.0|2.5| vs |120|1300|', -'njets_pt30_eta2p5 vs pT4l': '|0|1|2|3|20| vs |0|15|30|120|350| / |0|15|30|120|350| / |0|120|350| / |0|120|350|', -'massZ1 vs massZ2': '|50|80| vs |10|30| / |50|80| vs |30|60| / |80|110| vs |10|25| / |80|110| vs |25|30|', -'pT4l vs pTj1': '|0|350| vs |0|30| / |0|80| vs |30|60| / |80|350| vs |30|60| / |0|120| vs |60|120| / |120|350| vs |60|120| / |0|120| vs |120|350| / |120|350| vs |120|350|', ## Tricky, it is actually a 3D measurement? -#'pT4l vs pTHj': '|0|350| vs |0|30|', ## Tricky, it is actually a 3D measurement? -# 'mHj vs pTHj': , ## Tricky, it is actually a 3D measurement? -'njets_pt30_eta2p5 vs pTHj': '|0|1|10| vs |0|800| / |0|60|120|350|', -'njets_pt30_eta2p5 vs pTHjj': '|0|2|10| vs |0|800| / |0|60|120|350|', -'njets_pt30_eta2p5 vs mjj': '|0|2|10| vs |0|800| / |0|120|450|3000|', # Add bins of mjj -'njets_pt30_eta2p5 vs mHj': '|0|1|10| vs |0|800| / |120|180|220|300|400|600|2000|', -'njets_pt30_eta2p5 vs mHjj': '|0|2|10| vs |0|800| / |180|320|450|600|1000|2500|', # Add bins of mHjj -'njets_pt30_eta2p5 vs absdetajj': '|0|2|10| vs |0|800| / |0.0|1.0|2.5|9.0|', # Add bins of absdetajj -'pTj1 vs pTj2': '|0|30| vs |0|30| / |30|60| vs |30|60| / |60|120| vs |60|120| / |120|350| vs |120|350|', -'njets_pt30_eta2p5 vs TBjmax': '|0|1|10| vs |-1|1| / |0|30|60|80|100|120|1000|', -'njets_pt30_eta2p5 vs TCjmax': '|0|1|10| vs |-1|1| / |0|15|30|50|70|90|1000|' -# 'njets_pt30_eta2p5 vs absdphijj': '|0|2|10| vs |0|800| / ||' # Decide if we want to go to 2PI or to PI -} - - -def binning(var): - obsBins_input = list[var] - print('input', obsBins_input) - if not 'vs' in obsBins_input: #It is not a double-differential analysis - obs_bins = {0:(obsBins_input.split("|")[1:(len(obsBins_input.split("|"))-1)]),1:['0','inf']}[obsBins_input=='inclusive'] - obs_bins = [float(i) for i in obs_bins] #Convert a list of str to a list of float - doubleDiff = False - print ('It is a single-differential measurement, binning', obs_bins) - else: #It is a double-differential analysis - doubleDiff = True - # The structure of obs_bins is: - # index of the dictionary is the number of the bin - # [obs_bins_low, obs_bins_high, obs_bins_low_2nd, obs_bins_high_2nd] - # The first two entries are the lower and upper bound of the first variable - # The second two entries are the lower and upper bound of the second variable - if obsBins_input.count('vs')==1 and obsBins_input.count('/')>=1: #Situation like this one '|0|1|2|3|20| vs |0|10|20|45|90|250| / |0|10|20|80|250| / |0|20|90|250| / |0|25|250|' - obs_bins_tmp = obsBins_input.split(" vs ") #['|0|1|2|3|20|', '|0|10|20|45|90|250| / |0|10|20|80|250| / |0|20|90|250| / |0|25|250|'] - obs_bins_1st = obs_bins_tmp[0].split('|')[1:len(obs_bins_tmp[0].split('|'))-1] #['0', '1', '2', '3', '20'] - obs_bins_1st = [float(i) for i in obs_bins_1st] #Convert a list of str to a list of float - obs_bins_tmp = obs_bins_tmp[1].split(' / ') #['|0|10|20|45|90|250|', '|0|10|20|80|250|', '|0|20|90|250|', '|0|25|250|'] - obs_bins_2nd = {} - for i in range(len(obs_bins_tmp)): #At the end of the loop -> obs_bins_2nd {0: ['0', '10', '20', '45', '90', '250'], 1: ['0', '10', '20', '80', '250'], 2: ['0', '20', '90', '250'], 3: ['0', '25', '250']} - obs_bins_2nd[i] = obs_bins_tmp[i].split('|')[1:len(obs_bins_tmp[i].split('|'))-1] - obs_bins_2nd[i] = [float(j) for j in obs_bins_2nd[i]] #Convert a list of str to a list of float - obs_bins = {} - k = 0 #Bin index - for i in range(len(obs_bins_1st)-1): - for j in range(len(obs_bins_2nd[i])-1): - obs_bins[k] = [] - obs_bins[k].append(obs_bins_1st[i]) - obs_bins[k].append(obs_bins_1st[i+1]) - obs_bins[k].append(obs_bins_2nd[i][j]) - obs_bins[k].append(obs_bins_2nd[i][j+1]) - k +=1 - elif obsBins_input.count('vs')>1 and obsBins_input.count('/')>1: #Situation like this one '|50|80| vs |10|30| / |50|80| vs |30|60| / |80|110| vs |10|25| / |80|110| vs |25|30|' - obs_bins_tmp = obsBins_input.split(' / ') #['|50|80| vs |10|30|', '|50|80| vs |30|60|', '|80|110| vs |10|25|', '|80|110| vs |25|30|'] - obs_bins_1st={} - obs_bins_2nd={} - obs_bins={} - for i in range(len(obs_bins_tmp)): #At the end of the loop -> obs_bins_1st {0: ['50', '80'], 1: ['50', '80'], 2: ['80', '110'], 3: ['80', '110']} and obs_bins_2nd {0: ['10', '30'], 1: ['30', '60'], 2: ['10', '25'], 3: ['25', '30']} - obs_bins_tmp_bis = obs_bins_tmp[i].split(' vs ') - obs_bins_1st[i] = obs_bins_tmp_bis[0].split('|')[1:len(obs_bins_tmp_bis[0].split('|'))-1] - obs_bins_1st[i] = [float(j) for j in obs_bins_1st[i]] #Convert a list of str to a list of float - obs_bins_2nd[i] = obs_bins_tmp_bis[1].split('|')[1:len(obs_bins_tmp_bis[1].split('|'))-1] - obs_bins_2nd[i] = [float(j) for j in obs_bins_2nd[i]] #Convert a list of str to a list of float - obs_bins[i] = [] - obs_bins[i].append(obs_bins_1st[i][0]) - obs_bins[i].append(obs_bins_1st[i][1]) - obs_bins[i].append(obs_bins_2nd[i][0]) - obs_bins[i].append(obs_bins_2nd[i][1]) - elif obsBins_input.count('vs')==1 and obsBins_input.count('/')==0: #Situation like this one '|0|1|2|3|20| vs |0|10|20|45|90|250|' - obs_bins_tmp = obsBins_input.split(" vs ") #['|0|1|2|3|20|', '|0|10|20|45|90|250|'] - obs_bins_1st = obs_bins_tmp[0].split('|')[1:len(obs_bins_tmp[0].split('|'))-1] #['0', '1', '2', '3', '20'] - obs_bins_1st = [float(i) for i in obs_bins_1st] #Convert a list of str to a list of float - obs_bins_2nd = obs_bins_tmp[1].split('|')[1:len(obs_bins_tmp[1].split('|'))-1] #['0', '10', '20', '45', '90', '250'] - obs_bins_2nd = [float(i) for i in obs_bins_2nd] #Convert a list of str to a list of float - obs_bins = {} - k = 0 #Bin index - for i in range(len(obs_bins_1st)-1): - for j in range(len(obs_bins_2nd)-1): - obs_bins[k] = [] - obs_bins[k].append(obs_bins_1st[i]) - obs_bins[k].append(obs_bins_1st[i+1]) - obs_bins[k].append(obs_bins_2nd[j]) - obs_bins[k].append(obs_bins_2nd[j+1]) - k +=1 - else: - print ('Problem in the definition of the binning') - quit() - print ('It is a double-differential measurement, binning for the 1st variable', obs_bins_1st, 'and for the 2nd variable', obs_bins_2nd) - print (obs_bins) - return obs_bins, doubleDiff diff --git a/inputs/observables.py b/inputs/observables.py deleted file mode 100644 index d31de98..0000000 --- a/inputs/observables.py +++ /dev/null @@ -1,178 +0,0 @@ -observables = {'rapidity4l': - {'obs_reco': 'ZZy', - 'obs_gen': 'GENrapidity4l'}, -'pT4l': - {'obs_reco': 'ZZPt', - 'obs_gen': 'GENpT4l'}, -'pT4l_kL': - {'obs_reco': 'ZZPt', - 'obs_gen': 'GENpT4l'}, -'massZ1': - {'obs_reco': 'Z1Mass', - 'obs_gen': 'GENmassZ1'}, -'massZ2': - {'obs_reco': 'Z2Mass', - 'obs_gen': 'GENmassZ2'}, -'mass4l': - {'obs_reco': 'ZZMass', - 'obs_gen': 'GENmass4l'}, -'njets_pt30_eta2p5': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_gen': 'GENnjets_pt30_eta2p5'}, -'njets_pt30_eta2p5 vs pTj1': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'pTj1', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENpTj1'}, -'njets_pt30_eta4p7 vs pTj1_eta4p7': - {'obs_reco': 'njets_pt30_eta4p7', - 'obs_reco_2nd': 'pTj1_eta4p7', - 'obs_gen': 'GENnjets_pt30_eta4p7', - 'obs_gen_2nd': 'GENpTj1_eta4p7'}, -'njets_pt30_eta2p5 vs pTj2': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'pTj2', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENpTj2'}, -'pTj1 vs pTj2': - {'obs_reco': 'pTj1', - 'obs_reco_2nd': 'pTj2', - 'obs_gen': 'GENpTj1', - 'obs_gen_2nd': 'GENpTj2'}, -'pTj1': - {'obs_reco': 'pTj1', - 'obs_gen': 'GENpTj1'}, -'pTj1_eta4p7': - {'obs_reco': 'pTj1_eta4p7', - 'obs_gen': 'GENpTj1_eta4p7'}, -'costhetastar': - {'obs_reco': 'costhetastar', - 'obs_gen': 'GENcosThetaStar'}, -'costhetaZ1': - {'obs_reco': 'helcosthetaZ1', - 'obs_gen': 'GENcosTheta1'}, -'costhetaZ2': - {'obs_reco': 'helcosthetaZ2', - 'obs_gen': 'GENcosTheta2'}, -'phi': - {'obs_reco': 'helphi', - 'obs_gen': 'GENPhi'}, -'phistar': - {'obs_reco': 'phistarZ1', - 'obs_gen': 'GENPhi1'}, -'TCjmax': { - 'obs_reco': 'TCjmax', - 'obs_gen': 'GENTCjmax' - }, -'TBjmax': { - 'obs_reco': 'TBjmax', - 'obs_gen': 'GENTBjmax' - }, -'njets_pt30_eta2p5 vs TCjmax': - {'obs_reco_2nd': 'TCjmax', - 'obs_gen_2nd': 'GENTCjmax', - 'obs_reco': 'njets_pt30_eta2p5', - 'obs_gen': 'GENnjets_pt30_eta2p5'}, -'njets_pt30_eta2p5 vs TBjmax': - {'obs_reco_2nd': 'TBjmax', - 'obs_gen_2nd': 'GENTBjmax', - 'obs_reco': 'njets_pt30_eta2p5', - 'obs_gen': 'GENnjets_pt30_eta2p5'}, -'rapidity4l vs pT4l': - {'obs_reco': 'ZZy', - 'obs_reco_2nd': 'ZZPt', - 'obs_gen': 'GENrapidity4l', - 'obs_gen_2nd': 'GENpT4l'}, -'njets_pt30_eta2p5 vs pT4l': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'ZZPt', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENpT4l'}, -'massZ1 vs massZ2': - {'obs_reco': 'Z1Mass', - 'obs_reco_2nd': 'Z2Mass', - 'obs_gen': 'GENmassZ1', - 'obs_gen_2nd': 'GENmassZ2'}, -'pT4l vs pTj1': - {'obs_reco': 'ZZPt', - 'obs_reco_2nd': 'pTj1', - 'obs_gen': 'GENpT4l', - 'obs_gen_2nd': 'GENpTj1'}, -'pT4l vs pTHj': - {'obs_reco': 'ZZPt', - 'obs_reco_2nd': 'pTHj', - 'obs_gen': 'GENpT4l', - 'obs_gen_2nd': 'GENpTHj'}, -'mHj vs pTHj': - {'obs_reco': 'mHj', - 'obs_reco_2nd': 'pTHj', - 'obs_gen': 'GENmHj', - 'obs_gen_2nd': 'GENpTHj'}, -'njets_pt30_eta2p5 vs pTHj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'pTHj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENpTHj'}, -'njets_pt30_eta2p5 vs pTHjj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'pTHjj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENpTHjj'}, -'mjj': {'obs_reco': 'mjj'}, -'njets_pt30_eta2p5 vs mjj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'mjj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENmjj'}, -'njets_pt30_eta2p5 vs mHj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'mHj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENmHj'}, -'njets_pt30_eta2p5 vs mHjj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'mHjj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENmHjj'}, -'njets_pt30_eta2p5 vs detajj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'detajj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENdetajj'}, -'njets_pt30_eta2p5 vs absdetajj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'absdetajj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENabsdetajj'}, -'njets_pt30_eta2p5 vs dphijj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'dphijj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENdphijj'}, -'njets_pt30_eta2p5 vs absdphijj': - {'obs_reco': 'njets_pt30_eta2p5', - 'obs_reco_2nd': 'absdphijj', - 'obs_gen': 'GENnjets_pt30_eta2p5', - 'obs_gen_2nd': 'GENabsdphijj'}, -'D0m': - {'obs_reco': 'D0m', - 'obs_gen': 'GEN_D0m'}, -'Dcp': - {'obs_reco': 'Dcp', - 'obs_gen': 'GEN_Dcp'}, -'D0hp': - {'obs_reco': 'D0hp', - 'obs_gen': 'GEN_D0hp'}, -'Dcp': - {'obs_reco': 'Dcp', - 'obs_gen': 'GEN_Dcp'}, -'Dint': - {'obs_reco': 'Dint', - 'obs_gen': 'GEN_Dint'}, -'DL1': - {'obs_reco': 'DL1', - 'obs_gen': 'GEN_DL1'}, -'pTHj': - {'obs_reco': 'pTHj', - 'obs_gen': 'GENpTHj'} -} diff --git a/plotShapes.py b/plotShapes.py index 891779a..6ad57c5 100644 --- a/plotShapes.py +++ b/plotShapes.py @@ -63,13 +63,13 @@ def generateName(_year, _fStateNumber, _recobin, _fState, _bin, _physicalModel, binName = "ch"+_year+"_ch"+_fStateNumber procName = _fState+"Bin"+str(_bin) else: - _obsName_v3 = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta2p5': 'NJ'} + _obsName_v3 = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} if _obsName not in _obsName_v3: _obsName_v3[_obsName] = _obsName - if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName: - _recobin_final = str(_observableBins[_recobin][0]).replace('.', 'p')+'_'+str(_observableBins[_recobin][1]).replace('.', 'p')+'_'+str(_observableBins[_recobin][2]).replace('.', 'p')+'_'+str(_observableBins[_recobin][3]).replace('.', 'p') - _genbin_final = str(_observableBins[_bin][0]).replace('.', 'p')+'_'+str(_observableBins[_bin][1]).replace('.', 'p')+'_'+str(_observableBins[_bin][2]).replace('.', 'p')+'_'+str(_observableBins[_bin][3]).replace('.', 'p') + if '_' in obsName and not 'floating' in obsName and not 'kL' in obsName and obsName != 'njets_pt30_eta4p7': + _recobin_final = str(_observableBins[_recobin][0]).replace('.', 'p').replace('-','m')+'_'+str(_observableBins[_recobin][1]).replace('.', 'p').replace('-','m')+'_'+str(_observableBins[_recobin][2]).replace('.', 'p').replace('-','m')+'_'+str(_observableBins[_recobin][3]).replace('.', 'p').replace('-','m') + _genbin_final = str(_observableBins[_bin][0]).replace('.', 'p').replace('-','m')+'_'+str(_observableBins[_bin][1]).replace('.', 'p').replace('-','m')+'_'+str(_observableBins[_bin][2]).replace('.', 'p').replace('-','m')+'_'+str(_observableBins[_bin][3]).replace('.', 'p').replace('-','m') else: _recobin_final = str(_observableBins[_recobin]).replace('.', 'p').replace('-','m')+'_'+str(_observableBins[_recobin+1]).replace('.', 'p').replace('-','m') if int(_observableBins[_recobin+1]) > 1000: @@ -352,20 +352,33 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re ##### ------------------------ Data ------------------------ ##### CMS_channel = w.cat("CMS_channel") - mass = w.var("CMS_zz4l_mass").frame(RooFit.Bins(20)) + mass = w.var("CMS_zz4l_mass").frame(RooFit.Bins(30)) - - datacut = '' - bin_name, process_name = generateName(year, channel[fState], recobin, fState, bin, physicalModel, observableBins, obsName) - for year in ["1", "2", "3"]: - if(obsName!='mass4l' and obsName!='mass4l_zzfloating'): - datacut += "CMS_channel==CMS_channel::"+bin_name+" || " - else: - datacut += "CMS_channel==CMS_channel::"+bin_name+" || " - datacut = datacut.rstrip(" || ") - data = data.reduce(RooFit.Cut(datacut)) - data.plotOn(mass) - sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.ProjWData(data,True)) + if (fstate=="4l"): + datacut = '' + for year in ["1", "2", "3"]: + for fState in fStates: + bin_name, process_name = generateName(year, channel[fState], recobin, fState, bin, physicalModel, observableBins, obsName) + if(obsName!='mass4l' and obsName!='mass4l_zzfloating'): + datacut += "CMS_channel==CMS_channel::"+bin_name+" || " + else: + datacut += "CMS_channel==CMS_channel::"+bin_name+" || " + datacut = datacut.rstrip(" || ") + data = data.reduce(RooFit.Cut(datacut)) + data.plotOn(mass) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.ProjWData(data,True)) + else: + datacut = '' + for year in ["1", "2", "3"]: + bin_name, process_name = generateName(year, channel[fstate], recobin, fstate, bin, physicalModel, observableBins, obsName) + if(obsName!='mass4l' and obsName!='mass4l_zzfloating'): + datacut += "CMS_channel==CMS_channel::"+bin_name+" || " + else: + datacut += "CMS_channel==CMS_channel::"+bin_name+" || " + datacut = datacut.rstrip(" || ") + data = data.reduce(RooFit.Cut(datacut)) + data.plotOn(mass) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.ProjWData(data,True)) ##### ------------------------ Shapes ------------------------ ##### if (fstate!="4l"): @@ -373,7 +386,7 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re for bin in range(nBins): if bin==recobin: continue for year in ["1", "2", "3"]: - bin_name, process_name = generateName(year, channel[fState], recobin, fState, bin, physicalModel, observableBins, obsName) + bin_name, process_name = generateName(year, channel[fstate], recobin, fstate, bin, physicalModel, observableBins, obsName) comp_otherfid += "shapeSig_"+SignalNames[physicalModel]+process_name+"_"+bin_name+"," comp_otherfid = comp_otherfid.rstrip(',') @@ -382,7 +395,7 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re comp_zz = '' comp_zx = '' for year in ["1", "2", "3"]: - bin_name, process_name = generateName(year, channel[fState], recobin, fState, 0, physicalModel, observableBins, obsName) + bin_name, process_name = generateName(year, channel[fstate], recobin, fstate, 0, physicalModel, observableBins, obsName) comp_out += "shapeBkg_"+OutNames[physicalModel]+"_"+bin_name+"," comp_fake += "shapeBkg_"+CombNames[physicalModel]+"_"+bin_name+"," comp_zz += "shapeBkg_bkg_ggzz_"+bin_name+",shapeBkg_bkg_qqzz_"+bin_name+"," @@ -442,16 +455,21 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re dummy.SetLineWidth(0) dummy.SetMarkerSize(0) dummy.SetMarkerColor(0) - dummy.GetYaxis().SetTitle("Events / (2.33 GeV)") + dummy.GetYaxis().SetTitle("Events / (1.83 GeV)") dummy.GetXaxis().SetTitle("m_{"+fstate.replace("mu","#mu")+"} [GeV]") if (opt.UNBLIND): - dummy.SetMaximum(max(1.5*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - # if (obsName=="massZ2" and recobin==0): dummy.SetMaximum(max(3.0*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),3.5)) - else: + # dummy.SetMaximum(max(1.5*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) if fstate=='4e': dummy.SetMaximum(max(1*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) elif fstate=='4l': dummy.SetMaximum(max(0.2*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) else: dummy.SetMaximum(max(0.5*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) if (obsName=="massZ2" and recobin==0): dummy.SetMaximum(max(3.0*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),3.5)) + else: + # if fstate=='4e': dummy.SetMaximum(max(1*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) + # elif fstate=='4l': dummy.SetMaximum(max(0.2*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) + # else: dummy.SetMaximum(max(0.5*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) + # if (obsName=="massZ2" and recobin==0): dummy.SetMaximum(max(3.0*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),3.5)) + dummy.SetMaximum(max(1*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) + if (obsName=="massZ2" and recobin==0): dummy.SetMaximum(max(3.0*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),3.5)) #dummy.SetMaximum(0.5*max(n_trueH_asimov[fstate],n_zz_asimov[fstate],2.5)) dummy.SetMinimum(0.0) dummy.Draw() @@ -529,8 +547,8 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re elif (obsName=="nJets" or obsName=="njets_reco_pt30_eta4p7"): label = "N(jets) |#eta|<4.7" unit = "" - elif (obsName=="njets_pt30_eta2p5"): - label = "N(jets) |#eta|<2.5" + elif (obsName=="njets_pt30_eta4p7"): + label = "N(jets) |#eta|<4.7" unit = "" elif (obsName=='pt_leadingjet_pt30_eta4p7'): label = "p_{T}(jet)" @@ -635,7 +653,6 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re if obsName.startswith("mass4l"): PhysicalModels = ['v2','v3'] - # PhysicalModels = ['v3'] elif obsName == 'D0m' or obsName == 'Dcp' or obsName == 'D0hp' or obsName == 'Dint' or obsName == 'DL1' or obsName == 'DL1Zg' or obsName == 'costhetaZ1' or obsName == 'costhetaZ2'or obsName == 'costhetastar' or obsName == 'phi' or obsName == 'phistar' or obsName == 'massZ1' or obsName == 'massZ2': PhysicalModels = ['v4','v3'] elif 'kL' in obsName: diff --git a/producePlots.py b/producePlots.py index 83905cb..06839bf 100644 --- a/producePlots.py +++ b/producePlots.py @@ -871,6 +871,7 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): print 'a_ggH_powheg_hi',a_ggH_powheg_unc_hi print 'a_ggH_powheg_lo',a_ggH_powheg_unc_lo + a_ggH_minloHJ = array('d',[ggH_minloHJ[i] for i in range(len(ggH_minloHJ))]) v_ggH_minloHJ = TVectorD(len(a_ggH_minloHJ),a_ggH_minloHJ) a_ggH_minloHJ_unc_hi = array('d',[ggH_minloHJ_NNLOunc_hi[i] for i in range(len(ggH_minloHJ_unc_hi))]) @@ -953,6 +954,7 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): v_ggH_powheg_unc_hi = TVectorD(len(a_ggH_powheg_unc_hi),a_ggH_powheg_unc_hi) v_ggH_powheg_unc_lo = TVectorD(len(a_ggH_powheg_unc_lo),a_ggH_powheg_unc_lo) + a_ggH_minloHJ = array('d',[ggH_minloHJ[i]/((float(obs_bins_boundaries[i][1])-float(obs_bins_boundaries[i][0]))*(float(obs_bins_boundaries[i][3])-float(obs_bins_boundaries[i][2]))) for i in range(len(ggH_minloHJ))]) v_ggH_minloHJ = TVectorD(len(a_ggH_minloHJ),a_ggH_minloHJ) a_ggH_minloHJ_unc_hi = array('d',[ggH_minloHJ_NNLOunc_hi[i]/((float(obs_bins_boundaries[i][1])-float(obs_bins_boundaries[i][0]))*(float(obs_bins_boundaries[i][3])-float(obs_bins_boundaries[i][2]))) for i in range(len(ggH_minloHJ_unc_hi))]) @@ -960,6 +962,7 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): v_ggH_minloHJ_unc_hi = TVectorD(len(a_ggH_minloHJ_unc_hi),a_ggH_minloHJ_unc_hi) v_ggH_minloHJ_unc_lo = TVectorD(len(a_ggH_minloHJ_unc_lo),a_ggH_minloHJ_unc_lo) + ''' a_ggH_HRes = array('d',[ggH_HRes[i]/((float(obs_bins_boundaries[i][1])-float(obs_bins_boundaries[i][0]))*(float(obs_bins_boundaries[i][3])-float(obs_bins_boundaries[i][2]))) for i in range(len(ggH_HRes))]) v_ggH_HRes = TVectorD(len(a_ggH_HRes),a_ggH_HRes) @@ -1011,9 +1014,8 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): elif (obsName=="massZ2"): offset=20.0 elif (obsName=="massZ1"): offset=20.0 elif (obsName=="rapidity4l"): offset=20.0 - elif (obsName=="pTj1jet_pt30_eta2p5"): offset=30.0 elif (obsName=="pTj1"): offset=30.0 - elif (obsName=="njets_pt30_eta2p5"): offset=999.0 + elif (obsName=="njets_pt30_eta4p7"): offset=999.0 else: offset = 20.0 a_observable = array('d',[0.5*(float(obs_bins[i])+float(obs_bins[i+1])) for i in range(len(obs_bins)-1)]) v_observable = TVectorD(len(a_observable),a_observable) @@ -1111,24 +1113,24 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): v_data_lo_allunc = TVectorD(len(data_lo_allunc), array('d',[data_lo_allunc[i] for i in range(len(data_lo_allunc))])) - g_ggH_powheg = TGraphAsymmErrors(v_observable_1,v_ggH_powheg,v_dobservable_1,v_dobservable_1,v_ggH_powheg_unc_lo,v_ggH_powheg_unc_hi) - g_ggH_powheg.SetFillStyle(3254); - g_ggH_powheg.SetFillColor(ROOT.kAzure+2) - g_ggH_powheg.SetLineColor(ROOT.kAzure+2) - g_ggH_powheg.SetLineWidth(2) - g_ggH_powheg.SetMarkerColor(ROOT.kAzure+2) - - g_ggH_powhegBorder = TGraphAsymmErrors(v_observable_1,v_ggH_powheg,v_dobservable_1,v_dobservable_1,v_ggH_powheg_unc_lo,v_ggH_powheg_unc_hi) - g_ggH_powhegBorder.SetFillStyle(0) - g_ggH_powhegBorder.SetFillColor(ROOT.kAzure+2) - g_ggH_powhegBorder.SetLineColor(ROOT.kAzure+2) - g_ggH_powhegBorder.SetMarkerColor(ROOT.kAzure+2) - - g_ggH_powhege0 = TGraphAsymmErrors(v_observable,v_ggH_powheg,v_dobservable,v_dobservable,v_zeros,v_zeros) - g_ggH_powhege0.SetFillStyle(3254); - g_ggH_powhege0.SetFillColor(ROOT.kAzure+2) - g_ggH_powhege0.SetLineColor(ROOT.kAzure+2) - g_ggH_powhege0.SetMarkerColor(ROOT.kAzure+2) + # g_ggH_powheg = TGraphAsymmErrors(v_observable_1,v_ggH_powheg,v_dobservable_1,v_dobservable_1,v_ggH_powheg_unc_lo,v_ggH_powheg_unc_hi) + # g_ggH_powheg.SetFillStyle(3254); + # g_ggH_powheg.SetFillColor(ROOT.kAzure+2) + # g_ggH_powheg.SetLineColor(ROOT.kAzure+2) + # g_ggH_powheg.SetLineWidth(2) + # g_ggH_powheg.SetMarkerColor(ROOT.kAzure+2) + # + # g_ggH_powhegBorder = TGraphAsymmErrors(v_observable_1,v_ggH_powheg,v_dobservable_1,v_dobservable_1,v_ggH_powheg_unc_lo,v_ggH_powheg_unc_hi) + # g_ggH_powhegBorder.SetFillStyle(0) + # g_ggH_powhegBorder.SetFillColor(ROOT.kAzure+2) + # g_ggH_powhegBorder.SetLineColor(ROOT.kAzure+2) + # g_ggH_powhegBorder.SetMarkerColor(ROOT.kAzure+2) + # + # g_ggH_powhege0 = TGraphAsymmErrors(v_observable,v_ggH_powheg,v_dobservable,v_dobservable,v_zeros,v_zeros) + # g_ggH_powhege0.SetFillStyle(3254); + # g_ggH_powhege0.SetFillColor(ROOT.kAzure+2) + # g_ggH_powhege0.SetLineColor(ROOT.kAzure+2) + # g_ggH_powhege0.SetMarkerColor(ROOT.kAzure+2) if (obsName.startswith("mass4l")): @@ -1156,14 +1158,16 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): g_ggH_powheg.SetLineWidth(2) g_ggH_powheg.SetMarkerColor(ROOT.kAzure+2) + g_ggH_powhegBorder = TGraphAsymmErrors(v_observable_1,v_ggH_powheg,v_dobservable_1,v_dobservable_1,v_ggH_powheg_unc_lo,v_ggH_powheg_unc_hi) g_ggH_powhegBorder.SetFillStyle(0) g_ggH_powhegBorder.SetFillColor(ROOT.kAzure+2) g_ggH_powhegBorder.SetLineColor(ROOT.kAzure+2) g_ggH_powhegBorder.SetMarkerColor(ROOT.kAzure+2) - h_ggH_powheg = TH1D("h_ggH_powheg","h_ggH_powheg",4, 0, 4) - for i in range(4): + # h_ggH_powheg = TH1D("h_ggH_powheg","h_ggH_powheg",4, 0, 4) + h_ggH_powheg = TH1D("h_ggH_powheg","h_ggH_powheg",nBins-1, 0, nBins-1) + for i in range(nBins-1): h_ggH_powheg.SetBinContent(i+1,v_ggH_powheg[i]) else: g_ggH_powheg = TGraphAsymmErrors(v_observable_1,v_ggH_powheg,v_dobservable_1,v_dobservable_1,v_ggH_powheg_unc_lo,v_ggH_powheg_unc_hi) @@ -1225,7 +1229,6 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): h_ggH_ACbis = TH1D("h_ggH_ACbis","h_ggH_ACbis",nBins-1, array('d',[float(obs_bins[i]) for i in range(len(obs_bins))]) ) for i in range(nBins-1): h_ggH_ACbis.SetBinContent(i+1,v_ggH_ACbis[i]) - h_ggH_powheg.SetLineColor(ROOT.kAzure+2) h_ggH_powheg.SetLineWidth(2) @@ -1288,6 +1291,7 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): for i in range(nBins-1): h_ggH_minloHJ.SetBinContent(i+1,v_ggH_minloHJ[i]) + h_ggH_minloHJ.SetLineColor(ROOT.kOrange+2) h_ggH_minloHJ.SetLineWidth(2) @@ -1468,8 +1472,8 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): g_ratio_powhegBorder.SetLineColor(ROOT.kAzure+2) g_ratio_powhegBorder.SetMarkerColor(ROOT.kAzure+2) - h_ratio_powheg = TH1D("h_ratio_powheg","h_ratio_powheg",4, 0, 4) - for i in range(4): + h_ratio_powheg = TH1D("h_ratio_powheg","h_ratio_powheg",nBins-1, 0, nBins-1) + for i in range(nBins-1): h_ratio_powheg.SetBinContent(i+1,v_ratio_powheg[i]) else: g_ratio_powheg = TGraphAsymmErrors(v_observable_1,v_ratio_powheg,v_dobservable_1,v_dobservable_1,v_ratio_powheg_lo,v_ratio_powheg_hi) @@ -1550,8 +1554,8 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): g_ratio_minloHJBorder.SetLineColor(ROOT.kOrange+2) g_ratio_minloHJBorder.SetMarkerColor(ROOT.kOrange+2) - h_ratio_minloHJ = TH1D("h_ratio_minloHJ","h_ratio_minloHJ",4, 0, 4) - for i in range(4): + h_ratio_minloHJ = TH1D("h_ratio_minloHJ","h_ratio_minloHJ",nBins-1, 0, nBins-1) + for i in range(nBins-1): h_ratio_minloHJ.SetBinContent(i+1,v_ratio_minloHJ[i]) elif acFlag: # We keep the same names as in the last else (just to reduce entropy with if statements) g_ratio_minloHJ = TGraphAsymmErrors(v_observable_2,v_ratio_AC,v_dobservable_2,v_dobservable_2,v_ratio_AC_lo,v_ratio_AC_hi) @@ -1645,6 +1649,7 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): g_data_allunc.SetMarkerStyle(20) g_data_allunc.SetMarkerSize(1.2) + if (obsName=="pT4l"): label="p_{T}(H)" unit="GeV" @@ -1657,7 +1662,7 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): elif (obsName=="nJets" or obsName=="njets_pt30_eta4p7"): label = "N(jets)" unit = "" - elif (obsName=="njets_pt30_eta2p5"): + elif (obsName=="njets_pt30_eta4p7"): label = "N(jets)" unit = "" elif (obsName=="pTj1"): @@ -1719,13 +1724,24 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): elif obsName=="DL1": label = "D_{#Lambda 1}^{dec}" unit = "" + elif obsName=="DL1Zg": + label = "D_{#Lambda 1Zg}^{dec}" + unit = "" + elif obsName=="njets_pt30_eta4p7_pT4l": + label = "N_{j}" + label_2nd = "p_{T}(H)" + unit = "" + unit_2nd = "GeV" else: label = obsName + label_2nd = "" unit = "" + unit_2nd = "" c = TCanvas("c",obsName, 1400, 1400) if(opt.SETLOG): c.SetLogy() - if (not obsName.startswith("mass4l")): c.SetBottomMargin(0.35) + # if (not obsName.startswith("mass4l")): + c.SetBottomMargin(0.35) c.SetRightMargin(0.04) c.SetTopMargin(0.07) c.SetLeftMargin(0.18) @@ -1773,9 +1789,10 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): elif (obsName=="Dcp"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) elif (obsName=="D0hp"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) elif (obsName=="DL1"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) + elif (obsName=="DL1Zg"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) else: dummy.SetMaximum(55.0*max(max(a_data),(max(a_ggH_powheg)))) else: - if (obsName.startswith("mass4l")) or doubleDiff: dummy.SetMaximum(1.6*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) + if doubleDiff: dummy.SetMaximum(1.6*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) if doubleDiff: dummy.SetMaximum(2.0*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) elif (obsName=="massZ2"): dummy.SetMaximum(2.0*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) #AT elif (obsName=="massZ1"): dummy.SetMaximum(2.0*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) #AT @@ -1783,7 +1800,8 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): elif (obsName=="Dcp"): dummy.SetMaximum(2.0*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) #AT elif (obsName=="D0hp"): dummy.SetMaximum(2.5*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) #AT elif (obsName=="DL1"): dummy.SetMaximum(2.8*max(max(a_data),(max(a_ggH_powheg)))) - else: dummy.SetMaximum(1.5*(max(max(a_ggH_powheg),(max(a_data)+max(a_data_hi))))) + elif (obsName=="DL1Zg"): dummy.SetMaximum(2.8*max(max(a_data),(max(a_ggH_powheg)))) + else: dummy.SetMaximum(1.7*(max(max(a_ggH_powheg),(max(a_data)+max(a_data_hi))))) if (opt.SETLOG): dummy.SetMinimum(0.0601*max(min(a_data),(min(a_ggH_powheg)))) if (obsName.startswith("pTj1")): dummy.SetMinimum(0.0801*max(min(a_data),(min(a_ggH_powheg)))) @@ -1794,12 +1812,12 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): dummy.SetMarkerSize(0) dummy.GetXaxis().SetLabelSize(0.0) dummy.GetYaxis().SetLabelSize(0.04) - if (opt.SETLOG and (obsName.startswith('njets') or obsName.startswith('pTj1'))): + if (opt.SETLOG and (obsName.startswith('njets') or obsName.startswith('pTj1')) and not doubleDiff): dummy.SetMaximum(200.0*max(max(a_data),(max(a_ggH_powheg)))) dummy.GetXaxis().SetTitle("") else: dummy.GetXaxis().SetTitle("") - if (obsName.startswith('njets')): + if (obsName.startswith('njets') and not doubleDiff): dummy.GetXaxis().SetTitle('') dummy.GetXaxis().SetLabelSize(0.0) dummy.GetXaxis().SetBinLabel(1,'') @@ -1807,12 +1825,14 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): dummy.GetXaxis().SetBinLabel(3,'') dummy.GetXaxis().SetBinLabel(4,'') if (obsName.startswith("mass4l")): - dummy.GetXaxis().SetLabelSize(0.08) + dummy.GetXaxis().SetLabelSize(0.0) dummy.GetXaxis().SetBinLabel(1,'4l') dummy.GetXaxis().SetBinLabel(2,'2e2#mu') dummy.GetXaxis().SetBinLabel(3,'4#mu') dummy.GetXaxis().SetBinLabel(4,'4e') if doubleDiff: + if opt.SETLOG: + dummy.SetMaximum(1000000.0*max(max(a_data),(max(a_ggH_powheg)))) dummy.GetXaxis().SetLabelSize(0.08) dummy.GetXaxis().SetTitle('') for i in range(1,nBins-1): @@ -1824,24 +1844,24 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): # dummy.GetXaxis().SetBinLabel(5,'') dummy.GetXaxis().SetTitleSize(0.0) + # dummy.GetXaxis().SetRangeUser(-3.14,3.14) if (label=="inclusive"): dummy.GetYaxis().SetTitle("#sigma_{fid} (fb)") - elif (unit==''): - dummy.GetYaxis().SetTitle("d#sigma_{fid }/d"+label+" (fb)") elif doubleDiff: if unit==unit_2nd: dummy.GetYaxis().SetTitle("d^{2}#sigma_{fid }/d"+label+"d"+label_2nd+" (fb/"+unit+"^{2})") else: dummy.GetYaxis().SetTitle("d^{2}#sigma_{fid }/d"+label+"d"+label_2nd+" (fb/"+unit+unit_2nd+")") dummy.GetYaxis().SetTitleSize(0.04) + elif (unit==''): + dummy.GetYaxis().SetTitle("d#sigma_{fid }/d"+label+" (fb)") else: dummy.GetYaxis().SetTitle("d#sigma_{fid }/d"+label+" (fb/"+unit+")") - if (obsName.startswith('njets')): + if (obsName.startswith('njets') and not doubleDiff): dummy.GetYaxis().SetTitle("#sigma_{fid} (fb)") - dummy.GetYaxis().SetTitleOffset(1.4) - #dummy.GetYaxis().SetTitleOffset(1.8) + dummy.GetYaxis().SetTitleOffset(1.2) dummy.Draw("hist") h_XH.SetFillColor(kGreen-8) @@ -1862,13 +1882,16 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): legend . AddEntry(g_ggH_powheg , "gg#rightarrowH (POWHEG + JHUGen) + XH", "lf") else: legend . AddEntry(g_ggH_powheg , "SM (POWHEG + JHUGen + Pythia)", "lf") - if (obsName=="D0m"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa3=1 + Pythia)", "lf") - if (obsName=="Dcp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa3=0.5 + Pythia)", "lf") - if (obsName=="Dint"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa2=0.5 + Pythia)", "lf") - if (obsName=="D0hp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa2=1 + Pythia)", "lf") + if (obsName=="D0m"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a3}=1 + Pythia)", "lf") + if (obsName=="Dcp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a3}=0.5 + Pythia)", "lf") + if (obsName=="Dint"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a2}=0.5 + Pythia)", "lf") + if (obsName=="D0hp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a2}=1 + Pythia)", "lf") if (obsName=="DL1"): - legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen faL1=1 + Pythia)", "lf") - legend . AddEntry(g_ggH_ACbis , "AC (POWHEG + JHUGen faL1=0.5 + Pythia)", "lf") + legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{#Lambda 1}=1 + Pythia)", "lf") + legend . AddEntry(g_ggH_ACbis , "AC (POWHEG + JHUGen f_{#Lambda 1}=0.5 + Pythia)", "lf") + if (obsName=="DL1Zg"): + legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{#Lambda 1}^{Z#gamma}=1 + Pythia)", "lf") + legend . AddEntry(g_ggH_ACbis , "AC (POWHEG + JHUGen f_{#Lambda 1}^{Z#gamma}=0.5 + Pythia)", "lf") #legend . AddEntry(g_XH , "XH = VBF + VH + ttH", "l") if not acFlag: legend . AddEntry(h_XH , "XH = VBF + VH + ttH (POWHEG + JHUGen)", "f") legend . AddEntry(dummy, "(LHCHWG YR4, m_{H}="+opt.THEORYMASS+" GeV)", "") @@ -1963,118 +1986,118 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): dummy.Draw("axissame") - if (not obsName.startswith("mass4l")): - pad = TPad("pad", "pad", 0.0, 0.0, 1.0, 1.0) - pad.SetTopMargin(0.65) - pad.SetRightMargin(0.04) - pad.SetLeftMargin(0.18) - pad.SetFillColor(0) - pad.SetGridy(1) - pad.SetFillStyle(0) - pad.Draw() - pad.cd(0) - - if (obsName.startswith("mass4l")): - dummy2 = TH1D("dummy2","dummy2", 4, 0, 4) - for i in range(1,5): - dummy2.SetBinContent(i,1.02) - dummy2.GetXaxis().SetLabelSize(0.08) - dummy2.GetXaxis().SetBinLabel(1,'4l') - dummy2.GetXaxis().SetBinLabel(2,'2e2#mu') - dummy2.GetXaxis().SetBinLabel(3,'4#mu') - dummy2.GetXaxis().SetBinLabel(4,'4e') - elif doubleDiff: - dummy2 = TH1D("dummy2","dummy2", nBins-1, 0, nBins-1) - for i in range(1,nBins-1): + # if (not obsName.startswith("mass4l")): + pad = TPad("pad", "pad", 0.0, 0.0, 1.0, 1.0) + pad.SetTopMargin(0.65) + pad.SetRightMargin(0.04) + pad.SetLeftMargin(0.18) + pad.SetFillColor(0) + pad.SetGridy(1) + pad.SetFillStyle(0) + pad.Draw() + pad.cd(0) + + if (obsName.startswith("mass4l")): + dummy2 = TH1D("dummy2","dummy2", 4, 0, 4) + for i in range(1,5): + dummy2.SetBinContent(i,1.02) + dummy2.GetXaxis().SetLabelSize(0.06) + dummy2.GetXaxis().SetBinLabel(1,'4l') + dummy2.GetXaxis().SetBinLabel(2,'2e2#mu') + dummy2.GetXaxis().SetBinLabel(3,'4#mu') + dummy2.GetXaxis().SetBinLabel(4,'4e') + elif doubleDiff: + dummy2 = TH1D("dummy2","dummy2", nBins-1, 0, nBins-1) + for i in range(1,nBins-1): + dummy2.SetBinContent(i,1.02) + dummy2.GetXaxis().SetTitle("") + dummy2.GetXaxis().SetLabelSize(0.035) + for i in range(1,nBins): + dummy2.GetXaxis().SetBinLabel(i,'Bin '+str(i-1)) + # dummy2.GetXaxis().SetBinLabel(1,'0') + # dummy2.GetXaxis().SetBinLabel(2,'1') + # dummy2.GetXaxis().SetBinLabel(3,'2') + # dummy2.GetXaxis().SetBinLabel(4,'3') + # dummy2.GetXaxis().SetBinLabel(5,'4') + else: + if ("jet" in obsName and (not obsName.startswith("njets"))): + dummy2 = TH1D("dummy2","dummy2", int(float(obs_bins[nBins-1])-float(obs_bins[1])), float(obs_bins[1]), float(obs_bins[nBins-1])) + for i in range(int(float(obs_bins[nBins-1])-float(obs_bins[1]))): dummy2.SetBinContent(i,1.02) - dummy2.GetXaxis().SetTitle("") - dummy2.GetXaxis().SetLabelSize(0.08) - for i in range(1,nBins): - dummy2.GetXaxis().SetBinLabel(i,'Bin '+str(i-1)) - # dummy2.GetXaxis().SetBinLabel(1,'0') - # dummy2.GetXaxis().SetBinLabel(2,'1') - # dummy2.GetXaxis().SetBinLabel(3,'2') - # dummy2.GetXaxis().SetBinLabel(4,'3') - # dummy2.GetXaxis().SetBinLabel(5,'4') else: - if ("jet" in obsName and (not obsName.startswith("njets"))): - dummy2 = TH1D("dummy2","dummy2", int(float(obs_bins[nBins-1])-float(obs_bins[1])), float(obs_bins[1]), float(obs_bins[nBins-1])) - for i in range(int(float(obs_bins[nBins-1])-float(obs_bins[1]))): - dummy2.SetBinContent(i,1.02) + if (obsName=="pT4l"): + dummy2 = TH1D("dummy2","dummy2", int(float(obs_bins[nBins-1])-float(obs_bins[0])), float(obs_bins[0]), float(obs_bins[nBins-1])) else: - if (obsName=="pT4l"): - dummy2 = TH1D("dummy2","dummy2", int(float(obs_bins[nBins-1])-float(obs_bins[0])), float(obs_bins[0]), float(obs_bins[nBins-1])) - else: - dummy2 = TH1D("dummy2","dummy2", int(float(obs_bins[nBins-1])-float(obs_bins[0])), float(obs_bins[0]), float(obs_bins[nBins-1])) - for i in range(int(float(obs_bins[nBins-1])-float(obs_bins[0]))): - dummy2.SetBinContent(i,1.02) - if doubleDiff:dummy2.GetXaxis().SetLabelSize(0.08) - else: dummy2.GetXaxis().SetLabelSize(0.04) - - dummy2.SetLineColor(0) - dummy2.SetMarkerColor(0) - dummy2.SetLineWidth(0) - dummy2.SetMarkerSize(0) - if (obsName.startswith('njets')): - dummy2.GetXaxis().SetTitle(label) - dummy2.GetXaxis().SetLabelSize(0.08) - dummy2.GetXaxis().SetBinLabel(1,'0') - dummy2.GetXaxis().SetBinLabel(2,'1') - dummy2.GetXaxis().SetBinLabel(3,'2') - dummy2.GetXaxis().SetBinLabel(4,' 3') - dummy2.GetXaxis().SetBinLabel(5,'#geq 4') - elif (label=="inclusive"): - dummy2.GetXaxis().SetTitle("") - elif (unit==""): - dummy2.GetXaxis().SetTitle(label) - elif doubleDiff: - dummy2.GetXaxis().SetTitle("") - else: - dummy2.GetXaxis().SetTitle(label+" ("+unit+")") - - dummy2.GetYaxis().SetLabelSize(0.03) - dummy2.GetYaxis().SetNdivisions(10); - if (obsName.startswith("pTj1")): dummy2.GetYaxis().SetNdivisions(8); - dummy2.GetXaxis().SetNdivisions(510) - ratiomax=1.0 - for i in range(len(data_hi)): - if ((data[i]+data_hi[i])/ggH_minloHJ[i]>ratiomax): ratiomax=(data[i]+data_hi[i])/ggH_minloHJ[i] - dummy2.SetMaximum(ratiomax*1.1) - if (obsName=="pT4l"): - dummy2.SetMaximum(ratiomax*1.04) - dummy2.GetYaxis().SetNdivisions(508) - dummy2.SetMinimum(0.0) - dummy2.Draw("hist") - dummy2.GetYaxis().CenterTitle() - dummy2.GetYaxis().SetTitleSize(0.04) - dummy2.GetYaxis().SetTitleOffset(1.5) - dummy2.GetYaxis().SetTitleSize(0.03) - dummy2.GetYaxis().SetTitleOffset(2.0) - if acFlag: dummy2.GetYaxis().SetTitle('Ratio to SM') - else: dummy2.GetYaxis().SetTitle('Ratio to NNLOPS') - - - h_ratio_powheg.Draw("histsame") - h_ratio_minloHJ.Draw("histsame") - if acFlagBis: h_ratioBis_minloHJ.Draw("histsame") - - g_ratio_minloHJ.Draw("5same") - g_ratio_minloHJBorder.Draw("5same") - if acFlagBis: - g_ratioBis_minloHJ.Draw("5same") - g_ratioBis_minloHJBorder.Draw("5same") - #g_ratio_minloHJe0.Draw("epsame") - - g_ratio_powheg.Draw("5same") - g_ratio_powhegBorder.Draw("5same") - #g_ratio_powhege0.Draw("epsame") - - g_ratio_data.Draw("psameZ0") - g_ratio_sys.Draw("psameZ0") - g_ratio_datae0.Draw("psame") - - dummy2.GetYaxis().SetRangeUser(0,3) - dummy2.Draw("axissame") + dummy2 = TH1D("dummy2","dummy2", int(float(obs_bins[nBins-1])-float(obs_bins[0])), float(obs_bins[0]), float(obs_bins[nBins-1])) + for i in range(int(float(obs_bins[nBins-1])-float(obs_bins[0]))): + dummy2.SetBinContent(i,1.02) + dummy2.GetXaxis().SetLabelSize(0.04) + + dummy2.SetLineColor(0) + dummy2.SetMarkerColor(0) + dummy2.SetLineWidth(0) + dummy2.SetMarkerSize(0) + if (obsName.startswith('njets') and not doubleDiff): + dummy2.GetXaxis().SetTitle(label) + dummy2.GetXaxis().SetLabelSize(0.08) + dummy2.GetXaxis().SetBinLabel(1,'0') + dummy2.GetXaxis().SetBinLabel(2,'1') + dummy2.GetXaxis().SetBinLabel(3,'2') + dummy2.GetXaxis().SetBinLabel(4,' 3') + dummy2.GetXaxis().SetBinLabel(5,'#geq 4') + elif (label=="inclusive"): + dummy2.GetXaxis().SetTitle("") + elif doubleDiff: + dummy2.GetXaxis().SetTitle("") + elif (unit==""): + dummy2.GetXaxis().SetTitle(label) + else: + dummy2.GetXaxis().SetTitle(label+" ("+unit+")") + + dummy2.GetYaxis().SetLabelSize(0.03) + dummy2.GetYaxis().SetNdivisions(10); + if (obsName.startswith("pTj1")): dummy2.GetYaxis().SetNdivisions(8); + dummy2.GetXaxis().SetNdivisions(510) + ratiomax=1.0 + for i in range(len(data_hi)): + if ((data[i]+data_hi[i])/ggH_minloHJ[i]>ratiomax): ratiomax=(data[i]+data_hi[i])/ggH_minloHJ[i] + dummy2.SetMaximum(ratiomax*1.1) + if (obsName=="pT4l"): + dummy2.SetMaximum(ratiomax*1.04) + dummy2.GetYaxis().SetNdivisions(508) + dummy2.SetMinimum(0.0) + dummy2.Draw("hist") + dummy2.GetYaxis().CenterTitle() + dummy2.GetYaxis().SetTitleSize(0.04) + dummy2.GetYaxis().SetTitleOffset(1.5) + dummy2.GetYaxis().SetTitleSize(0.03) + dummy2.GetYaxis().SetTitleOffset(2.0) + if acFlag: dummy2.GetYaxis().SetTitle('Ratio to SM') + else: dummy2.GetYaxis().SetTitle('Ratio to NNLOPS') + + + h_ratio_powheg.Draw("histsame") + h_ratio_minloHJ.Draw("histsame") + if acFlagBis: h_ratioBis_minloHJ.Draw("histsame") + + g_ratio_minloHJ.Draw("5same") + g_ratio_minloHJBorder.Draw("5same") + if acFlagBis: + g_ratioBis_minloHJ.Draw("5same") + g_ratioBis_minloHJBorder.Draw("5same") + #g_ratio_minloHJe0.Draw("epsame") + + g_ratio_powheg.Draw("5same") + g_ratio_powhegBorder.Draw("5same") + #g_ratio_powhege0.Draw("epsame") + + g_ratio_data.Draw("psameZ0") + g_ratio_sys.Draw("psameZ0") + g_ratio_datae0.Draw("psame") + + dummy2.GetYaxis().SetRangeUser(0,3) + # dummy2.GetXaxis().SetRangeUser(-3.14,3.14) + dummy2.Draw("axissame") if (obsName=="pT4l"): box = TPaveText(240,-0.4,260,-0.1) @@ -2170,12 +2193,21 @@ def plotXS(obsName, obs_bins, obs_bins_boundaries = False): acFlagBis = True acSample = '0L1' acSampleBis = '0L1f05ph0' +elif obs_name == 'DL1Zg': + acFlag = True + acFlagBis = True + acSample = '0L1Zg' + acSampleBis = '0L1Zgf05ph0' else: acFlag = False acFlagBis = False if not doubleDiff: - if float(obs_bins[len(obs_bins)-1])>300.0: + if obs_name == 'TBjmax' or obs_name == 'mjj': + obs_bins[len(obs_bins)-1]='500.0' + elif obs_name == 'absdetajj': + obs_bins[0]='-1.0' + elif float(obs_bins[len(obs_bins)-1])>300.0: obs_bins[len(obs_bins)-1]='250.0' if (opt.OBSNAME=="nJets" or opt.OBSNAME.startswith("njets")): obs_bins[len(obs_bins)-1]='5' diff --git a/producePlots_v4.py b/producePlots_v4.py index f16ff97..fab71e9 100644 --- a/producePlots_v4.py +++ b/producePlots_v4.py @@ -61,16 +61,16 @@ def parseOptions(): def plotXS(obsName, obs_bins, chan, obs_bins_boundaries = False): - _temp = __import__('inputs_sig_'+obsName+'_'+opt.YEAR, globals(), locals(), ['acc'], -1) + _temp = __import__('inputs_sig_extrap_'+obsName+'_'+opt.YEAR, globals(), locals(), ['acc'], -1) acc = _temp.acc print 'inputs_sig_'+obsName+'_'+opt.YEAR # eff = _temp.eff # outinratio = _temp.outinratio if(opt.YEAR=='Full'): - _temp = __import__('inputs_sig_'+obsName+'_NNLOPS_Full', globals(), locals(), ['acc'], -1) + _temp = __import__('inputs_sig_extrap_'+obsName+'_NNLOPS_Full', globals(), locals(), ['acc'], -1) acc_NNLOPS = _temp.acc else: - _temp = __import__('inputs_sig_'+obsName+'_NNLOPS_'+opt.YEAR, globals(), locals(), ['acc'], -1) + _temp = __import__('inputs_sig_extrap_'+obsName+'_NNLOPS_'+opt.YEAR, globals(), locals(), ['acc'], -1) acc_NNLOPS = _temp.acc if(acFlag): @@ -1721,6 +1721,9 @@ def plotXS(obsName, obs_bins, chan, obs_bins_boundaries = False): elif obsName=="DL1": label = "D_{#Lambda 1}^{dec}" unit = "" + elif obsName=="DL1Zg": + label = "D_{#Lambda 1Zg}^{dec}" + unit = "" else: label = obsName unit = "" @@ -1778,6 +1781,7 @@ def plotXS(obsName, obs_bins, chan, obs_bins_boundaries = False): elif (obsName=="Dcp"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) elif (obsName=="D0hp"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) elif (obsName=="DL1"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) + elif (obsName=="DL1Zg"): dummy.SetMaximum(1100*max(max(a_data),(max(a_ggH_powheg)))) else: dummy.SetMaximum(55.0*max(max(a_data),(max(a_ggH_powheg)))) else: if (obsName=="mass4l") or doubleDiff: dummy.SetMaximum(1.6*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) @@ -1788,6 +1792,7 @@ def plotXS(obsName, obs_bins, chan, obs_bins_boundaries = False): elif (obsName=="Dcp"): dummy.SetMaximum(2.0*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) #AT elif (obsName=="D0hp"): dummy.SetMaximum(2.5*(max(max(a_data),(max(a_ggH_powheg)))+max(a_data_hi))) #AT elif (obsName=="DL1"): dummy.SetMaximum(2.8*max(max(a_data),(max(a_ggH_powheg)))) + elif (obsName=="DL1Zg"): dummy.SetMaximum(2.8*max(max(a_data),(max(a_ggH_powheg)))) else: dummy.SetMaximum(1.5*(max(max(a_ggH_powheg),(max(a_data)+max(a_data_hi))))) if (opt.SETLOG): dummy.SetMinimum(0.0601*max(min(a_data),(min(a_ggH_powheg)))) @@ -1867,13 +1872,16 @@ def plotXS(obsName, obs_bins, chan, obs_bins_boundaries = False): legend . AddEntry(g_ggH_powheg , "gg#rightarrowH (POWHEG + JHUGen) + XH", "lf") else: legend . AddEntry(g_ggH_powheg , "SM (POWHEG + JHUGen + Pythia)", "lf") - if (obsName=="D0m"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa3=1 + Pythia)", "lf") - if (obsName=="Dcp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa3=0.5 + Pythia)", "lf") - if (obsName=="Dint"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa2=0.5 + Pythia)", "lf") - if (obsName=="D0hp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen fa2=1 + Pythia)", "lf") + if (obsName=="D0m"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a3}=1 + Pythia)", "lf") + if (obsName=="Dcp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a3}=0.5 + Pythia)", "lf") + if (obsName=="Dint"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a2}=0.5 + Pythia)", "lf") + if (obsName=="D0hp"): legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{a2}=1 + Pythia)", "lf") if (obsName=="DL1"): - legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen faL1=1 + Pythia)", "lf") - legend . AddEntry(g_ggH_ACbis , "AC (POWHEG + JHUGen faL1=0.5 + Pythia)", "lf") + legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{#Lambda 1}=1 + Pythia)", "lf") + legend . AddEntry(g_ggH_ACbis , "AC (POWHEG + JHUGen f_{#Lambda 1}=0.5 + Pythia)", "lf") + if (obsName=="DL1Zg"): + legend . AddEntry(g_ggH_AC , "AC (POWHEG + JHUGen f_{#Lambda 1}^{Zg}=1 + Pythia)", "lf") + legend . AddEntry(g_ggH_ACbis , "AC (POWHEG + JHUGen f_{#Lambda 1}^{Zg}=0.5 + Pythia)", "lf") #legend . AddEntry(g_XH , "XH = VBF + VH + ttH", "l") if not acFlag: legend . AddEntry(h_XH , "XH = VBF + VH + ttH (POWHEG + JHUGen)", "f") legend . AddEntry(dummy, "(LHCHWG YR4, m_{H}="+opt.THEORYMASS+" GeV)", "") @@ -2180,6 +2188,11 @@ def plotXS(obsName, obs_bins, chan, obs_bins_boundaries = False): acFlagBis = True acSample = '0L1' acSampleBis = '0L1f05ph0' +elif obs_name == 'DL1Zg': + acFlag = True + acFlagBis = True + acSample = '0L1Zg' + acSampleBis = '0L1Zgf05ph0' else: acFlag = False acFlagBis = False diff --git a/templates/RunTemplates.py b/templates/RunTemplates.py index 4ad46d5..b5a2ed2 100644 --- a/templates/RunTemplates.py +++ b/templates/RunTemplates.py @@ -56,7 +56,7 @@ def checkDir(folder_path): # ------------------------------- FUNCTIONS TO GENERATE DATAFRAME FOR ggZZ AND qqZZ ---------------------------------------------------- # Weights for histogram def weight(df, xsec, gen, lumi, additional = None): - weight = (lumi * 1000 * xsec * df.overallEventWeight * df.L1prefiringWeight) / gen #Common structure + weight = (lumi * 1000 * xsec * df.overallEventWeight * df.SFcorr * df.L1prefiringWeight) / gen #Common structure if additional == 'ggH': weight *= df.ggH_NNLOPS_weight elif additional == 'qqzz': @@ -113,7 +113,7 @@ def generators(year): def add_njets(pt,eta): n = 0 for i in range(len(pt)): - if pt[i]>30 and abs(eta[i])<2.5: n=n+1 + if pt[i]>30 and abs(eta[i])<4.7: n=n+1 return n def add_leadjet(pt,eta): _pTj1 = 0.0 @@ -121,7 +121,7 @@ def add_leadjet(pt,eta): return _pTj1 else: for i in range(len(pt)): - if (pt[i]>30 and abs(eta[i])<2.5 and pt[i] > _pTj1): _pTj1 = pt[i] + if (pt[i]>30 and abs(eta[i])<4.7 and pt[i] > _pTj1): _pTj1 = pt[i] return _pTj1 # Rapidity @@ -163,7 +163,7 @@ def dataframes(year): 'overallEventWeight', 'L1prefiringWeight', 'JetPt', 'JetEta', 'costhetastar', 'helcosthetaZ1','helcosthetaZ2','helphi','phistarZ1', 'pTHj', 'TCjmax', 'TBjmax', 'mjj', 'pTj1', 'pTj2', 'mHj', 'mHjj', 'pTHjj', - 'njets_pt30_eta4p7', 'absdetajj', + 'njets_pt30_eta4p7', 'absdetajj', 'SFcorr', 'dphijj', 'Dcp', 'D0m', 'D0hp', 'Dint', 'DL1', 'DL1int', 'DL1Zg', 'DL1Zgint'] if (bkg == 'ZZTo4l'): b_bkg.append('KFactor_EW_qqZZ'); b_bkg.append('KFactor_QCD_qqZZ_M') @@ -405,7 +405,7 @@ def doTemplates(df_irr, df_red, binning, var, var_string, var_2nd='None'): w = np.asarray(w).astype('float') # ------ - if((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('eta' in obs_name) | acFlag): + if((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('deta' in obs_name) | acFlag): histo = ROOT.TH1D("m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), "m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), 20, opt.LOWER_BOUND, opt.UPPER_BOUND) elif doubleDiff and 'rapidity' in var_string: histo = ROOT.TH1D("m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+"_"+str(bin_low_2nd)+"_"+str(bin_high_2nd), "m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+'_'+str(bin_low_2nd)+"_"+str(bin_high_2nd), 20, opt.LOWER_BOUND, opt.UPPER_BOUND) @@ -418,7 +418,7 @@ def doTemplates(df_irr, df_red, binning, var, var_string, var_2nd='None'): histo.FillN(len(mass4l), mass4l, w) smoothAndNormaliseTemplate(histo, 1) - if ((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('eta' in obs_name) | acFlag): + if ((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('deta' in obs_name) | acFlag): outFile = ROOT.TFile.Open(str(year)+"/"+var_string+"/XSBackground_"+bkg+"_"+f+"_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+".root", "RECREATE") elif doubleDiff and 'rapidity' in var_string: outFile = ROOT.TFile.Open(str(year)+"/"+var_string+"/XSBackground_"+bkg+"_"+f+"_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+"_"+str(bin_low_2nd)+"_"+str(bin_high_2nd)+".root", "RECREATE") @@ -475,7 +475,7 @@ def doTemplates(df_irr, df_red, binning, var, var_string, var_2nd='None'): w = df['yield_SR'].to_numpy() w = np.asarray(w).astype('float') # ------ - if((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('eta' in obs_name) | acFlag): + if((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('deta' in obs_name) | acFlag): histo = ROOT.TH1D("m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), "m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), 20, opt.LOWER_BOUND, opt.UPPER_BOUND) elif doubleDiff and 'rapidity' in var_string: histo = ROOT.TH1D("m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+"_"+str(bin_low_2nd)+"_"+str(bin_high_2nd), "m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+str(bin_low_2nd)+"_"+str(bin_high_2nd), 20, opt.LOWER_BOUND, opt.UPPER_BOUND) @@ -485,7 +485,7 @@ def doTemplates(df_irr, df_red, binning, var, var_string, var_2nd='None'): histo = ROOT.TH1D("m4l_"+var_string+"_"+str(int(bin_low))+"_"+str(int(bin_high)), "m4l_"+var_string+"_"+str(int(bin_low))+"_"+str(int(bin_high)), 20, opt.LOWER_BOUND, opt.UPPER_BOUND) histo.FillN(len(mass4l), mass4l, w) smoothAndNormaliseTemplate(histo, 1) - if((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('eta' in obs_name) | acFlag): + if((obs_name == 'rapidity4l') | ('cos' in obs_name) | ('phi' in obs_name) | ('deta' in obs_name) | acFlag): outFile = ROOT.TFile.Open(str(year)+"/"+var_string+"/XSBackground_ZJetsCR_"+f+"_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+".root", "RECREATE") elif doubleDiff and 'rapidity' in var_string: outFile = ROOT.TFile.Open(str(year)+"/"+var_string+"/XSBackground_ZJetsCR_"+f+"_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+"_"+str(bin_low_2nd)+"_"+str(bin_high_2nd)+".root", "RECREATE") @@ -571,7 +571,7 @@ def doTemplates(df_irr, df_red, binning, var, var_string, var_2nd='None'): branches_ZX = ['ZZMass', 'Z1Flav', 'Z2Flav', 'LepLepId', 'LepEta', 'LepPt', 'Z1Mass', 'Z2Mass', 'ZZPt', 'ZZEta', 'JetPt', 'JetEta', 'costhetastar', 'helcosthetaZ1','helcosthetaZ2','helphi','phistarZ1', 'pTHj', 'TCjmax', 'TBjmax', 'mjj', 'pTj1', 'pTj2', 'mHj', 'mHjj', 'pTHjj', 'njets_pt30_eta4p7', - 'Dcp', 'D0m', 'D0hp', 'Dint', 'DL1', 'DL1int', 'DL1Zg', 'DL1Zgint','absdetajj'] + 'Dcp', 'D0m', 'D0hp', 'Dint', 'DL1', 'DL1int', 'DL1Zg', 'DL1Zgint','absdetajj', 'dphijj'] dfZX={} for year in years: g_FR_mu_EB, g_FR_mu_EE, g_FR_e_EB, g_FR_e_EE = openFR(year) diff --git a/templates/plot_templates.py b/templates/plot_templates.py index dc4e0b7..6f86603 100644 --- a/templates/plot_templates.py +++ b/templates/plot_templates.py @@ -40,7 +40,7 @@ 'DL1': True, 'DL1Zg': True, 'rapidity4l vs pT4l': True, -'njets vs pt30_eta4p7_pT4l': False, +'njets_pt30_eta4p7 vs pT4l': False, 'pTj1 vs pTj2': False, 'pT4l vs pTHj': False, 'massZ1 vs massZ2': False, @@ -227,7 +227,7 @@ def setLegendProperties(leg, sHeader = "skip", fillStyle=0, fillColor=0): #### qqZZZ + ggZZ +ZX #### h1D_dummy.SetMaximum(2.0*h1D_4mu[kBkg_qqZZ,iBin].GetMaximum()) h1D_dummy.Draw() - cmsPreliminary(c1, binRangeLeg[iBin]+" 2e2#mu "+year[iYear]) + cmsPreliminary(c1, binRangeLeg[iBin]+" 4#mu "+year[iYear]) leg1 = ROOT.TLegend(leg_xl,leg_yb,leg_xr,leg_yt) setLegendProperties(leg1) setHistProperties(h1D_4mu[kBkg_qqZZ,iBin],lineWidth,1,ROOT.kBlack) @@ -247,7 +247,7 @@ def setLegendProperties(leg, sHeader = "skip", fillStyle=0, fillColor=0): #### qqZZZ + ggZZ +ZX #### h1D_dummy.SetMaximum(2.0*h1D_4e[kBkg_qqZZ,iBin].GetMaximum()) h1D_dummy.Draw() - cmsPreliminary(c1, binRangeLeg[iBin]+" 2e2#mu "+year[iYear]) + cmsPreliminary(c1, binRangeLeg[iBin]+" 4e "+year[iYear]) leg1 = ROOT.TLegend(leg_xl,leg_yb,leg_xr,leg_yt) setLegendProperties(leg1) setHistProperties(h1D_4e[kBkg_qqZZ,iBin],lineWidth,1,ROOT.kBlack)