Including leading jet pT measurement

LHScans/plot_LLScan.py

@@ -58,7 +58,7 @@ def BuildScan(scan, param, files, color, yvals, ycut):
     bestfit = None
 
     for i in xrange(graph.GetN()):
-        if graph.GetY()[i] < 0.005:
+        if graph.GetY()[i] < 0.003:
             bestfit = graph.GetX()[i]
 
     if graph.GetY()[i] == 0:
@@ -173,15 +173,16 @@ def BuildScan(scan, param, files, color, yvals, ycut):
 if(obsName == 'pT4l'): label = 'p_{T}^{H} (GeV)'
 if(obsName == 'massZ1'): label = 'm_{Z1} (GeV)'
 if(obsName == 'massZ2'): label = 'm_{Z2} (GeV)'
-
+if(obsName == 'njets_pt30_eta2p5'): label = 'nJets, pT>30 GeV, |#eta|<2.5'
+if(obsName == 'pTj1'): label = 'p_{T}^{(Lead. jet)} (GeV)'
 obs_bins = {0:(opt.OBSBINS.split("|")[1:(len(opt.OBSBINS.split("|"))-1)]),1:['0','inf']}[opt.OBSBINS=='inclusive']
 obs_bins = [float(i) for i in obs_bins] #Convert a list of str to a list of float
 nBins = len(obs_bins)
 
 for i in range(nBins-1):
     _bin = i
     _obs_bin = _poi+str(i)
-
+    
     if obsName == 'mass4l':
         if _bin == 1:
             _obs_bin = 'r2e2muBin0'
@@ -223,8 +224,8 @@ def BuildScan(scan, param, files, color, yvals, ycut):
                 elif _bin == 3:
                     if obsName == 'mass4l':
                         graphs[ifile].SetPoint(ipoint,entry.r4eBin0,2.0*entry.deltaNLL)
-                    else:
-                        graphs[ifile].SetPoint(ipoint,entry.SigmaBin3,2.0*entry.deltaNLL)
+                    else:                   
+                        graphs[ifile].SetPoint(ipoint,entry.SigmaBin3,2.0*entry.deltaNLL)                
                     ipoint = ipoint+1
                 elif _bin == 4:
                     graphs[ifile].SetPoint(ipoint,entry.SigmaBin4,2.0*entry.deltaNLL)
@@ -239,7 +240,7 @@ def BuildScan(scan, param, files, color, yvals, ycut):
                     graphs[ifile].SetPoint(ipoint,entry.SigmaBin7,2.0*entry.deltaNLL)
                     ipoint = ipoint+1
 
-    c=TCanvas()
+    c=TCanvas("c", "c", 1000, 800)
     c.SetLeftMargin(0.14)
     c.SetRightMargin(0.08)
     c.cd()
@@ -312,12 +313,12 @@ def BuildScan(scan, param, files, color, yvals, ycut):
     leg.Draw("SAME")
 
     poi = _obs_bin
-    fname = inputPath + "higgsCombine_pT4l_"+poi+".MultiDimFit.mH125.38.123456.root"
+    fname = inputPath + "higgsCombine_"+obsName+"_"+poi+".MultiDimFit.mH125.38.123456.root"
     exp_scan = BuildScan('scan', poi, [fname], 2, yvals, 7.)
     exp_nom = exp_scan['val']
     exp_2sig = exp_scan['val_2sig']
 
-    fname = inputPath + "higgsCombine_pT4l_"+poi+"_NoSys_exp.MultiDimFit.mH125.38.root"
+    fname = inputPath + "higgsCombine_"+obsName+"_"+poi+"_NoSys_exp.MultiDimFit.mH125.38.123456.root"
     exp_scan_stat = BuildScan('scan', poi, [fname], 2, yvals, 7.)
     exp_nom_stat = exp_scan_stat['val']
     exp_2sig_stat = exp_scan_stat['val_2sig']
@@ -326,12 +327,12 @@ def BuildScan(scan, param, files, color, yvals, ycut):
     exp_do_sys = np.sqrt(abs(exp_nom[2])**2 - abs(exp_nom_stat[2])**2)
 
     if (opt.UNBLIND):
-        fname = inputPath + "higgsCombine_pT4l_"+poi+".MultiDimFit.mH125.38.root"
+        fname = inputPath + "higgsCombine_"+obsName+"_"+poi+".MultiDimFit.mH125.38.root"
         obs_scan = BuildScan('scan', poi, [fname], 2, yvals, 7.)
         obs_nom = obs_scan['val']
         obs_2sig = obs_scan['val_2sig']
 
-        fname = inputPath + "higgsCombine_pT4l_"+poi+"_NoSys.MultiDimFit.mH125.38.root"
+        fname = inputPath + "higgsCombine_"+obsName+"_"+poi+"_NoSys.MultiDimFit.mH125.38.root"
         obs_scan_stat = BuildScan('scan', poi, [fname], 2, yvals, 7.)
         obs_nom_stat = obs_scan_stat['val']
         obs_2sig_stat = obs_scan_stat['val_2sig']
@@ -363,7 +364,7 @@ def BuildScan(scan, param, files, color, yvals, ycut):
         Text4.Draw()
 
     Text = TPaveText(0.58, 0.88,0.93,0.95,'brNDC')
-    #Text.SetNDC()
+    #Text.SetNDC() 
     Text.SetTextAlign(31);
     Text.SetTextSize(0.03)
     leftText = "CMS"
@@ -388,7 +389,13 @@ def BuildScan(scan, param, files, color, yvals, ycut):
     latex2.SetTextSize(0.04)
     latex2.SetTextFont(42)
     latex2.SetTextAlign(31) # align right
-    latex2.DrawLatex(0.45,0.65, str(obs_bins[_bin])+' < '+label+' < '+str(obs_bins[_bin+1]))
+    if not 'jet' in obsName:
+	if 'pTj1' in obsName:
+		latex2.DrawLatex(0.55,0.65, str(obs_bins[_bin])+' < '+label+' < '+str(obs_bins[_bin+1]))
+	else:
+        	latex2.DrawLatex(0.45,0.65, str(obs_bins[_bin])+' < '+label+' < '+str(obs_bins[_bin+1]))
+    else:
+        latex2.DrawLatex(0.45,0.65, str(_bin)+' jet(s)')
     latex2.DrawLatex(0.91,0.22, "#scale[0.7]{#color[12]{68% CL}}")
     latex2.DrawLatex(0.91,0.52, "#scale[0.7]{#color[12]{95% CL}}")
 
@@ -410,7 +417,7 @@ def BuildScan(scan, param, files, color, yvals, ycut):
     c.Update()
     c.SaveAs("plots/lhscan_compare_"+obsName+"_"+str(_bin)+".pdf")
     c.SaveAs("plots/lhscan_compare_"+obsName+"_"+str(_bin)+".png")
-
+    
 with open('resultsXS_LHScan_expected_'+obsName+'_v3.py', 'w') as f:
     f.write('resultsXS = '+str(resultsXS_asimov)+' \n')
 

coefficients/RunCoefficients.py

@@ -72,7 +72,6 @@ def weight(df, fail, xsec, gen, lumi, additional = None):
         df['weight_histo_reco_NNLOPS'] = weight_histo_reco_NNLOPS #NNLOPS (only ggH)
     return df
 
-
 # Uproot to generate pandas
 def prepareTrees(year):
     d_sig = {}
@@ -185,11 +184,12 @@ def createDataframe(d_sig,fail,gen,xsec,signal,lumi):
     b_sig = ['EventNumber','GENmass4l', 'GENpT4l', 'GENrapidity4l', 'GENeta4l',
              'GENlep_id', 'GENlep_MomId', 'GENlep_MomMomId', 'GENlep_Hindex',
              'GENZ_DaughtersId', 'GENZ_MomId', 'passedFiducialSelection_bbf',
-             'PUWeight', 'genHEPMCweight','GENnjets_pt30_eta2p5']
+             'PUWeight', 'genHEPMCweight','GENnjets_pt30_eta2p5', 
+             'GenCleanedJetPt', 'GenCleanedJetEta', 'GENpTj1']
     if signal == 'ggH125': b_sig.append('ggH_NNLOPS_weight') #Additional entry for the weight in case of ggH
     if not fail: b_sig.extend(['ZZMass', 'ZZPt', 'ZZy', 'Z1Mass', 'Z2Mass', 'ZZEta', 'Z1Flav', 'Z2Flav',
                           'lep_genindex', 'lep_Hindex', 'overallEventWeight', 'L1prefiringWeight','dataMCWeight', 'trigEffWeight', 'njets_pt30_eta2p5',
-                          'njets_pt30_eta2p5_jesup', 'njets_pt30_eta2p5_jesdn']) #Additioanl entries for passing events
+                          'njets_pt30_eta2p5_jesup', 'njets_pt30_eta2p5_jesdn', 'pTj1']) #Additioanl entries for passing events
     df = d_sig.pandas.df(b_sig, flatten = False)
     if fail: #Negative branches for failed events (it is useful when creating fiducial pandas)
         df['ZZMass'] = -1
@@ -230,6 +230,7 @@ def createDataframe(d_sig,fail,gen,xsec,signal,lumi):
     else:
         df = weight(df, fail, xsec, gen, lumi, 'ggH')
         df = df.drop(columns=['ggH_NNLOPS_weight'])
+
     return df
 
 
@@ -521,9 +522,12 @@ def doGetCoeff(obs_reco, obs_gen, obs_name, obs_bins, type):
 elif(obs_name == 'mass4l'):
     obs_reco = 'ZZMass'
     obs_gen = 'GENmass4l'
-if (obs_name == "njets_pt30_eta2p5"):
+elif(obs_name == "njets_pt30_eta2p5"):
     obs_reco = "njets_pt30_eta2p5"
     obs_gen = "GENnjets_pt30_eta2p5"
+elif(obs_name == 'pTj1'):
+    obs_reco = 'pTj1'
+    obs_gen = 'GENpTj1' 
 
 # Generate dataframes
 d_sig = {}
@@ -587,3 +591,4 @@ def doGetCoeff(obs_reco, obs_gen, obs_name, obs_bins, type):
 acceptance = {}
 err_acceptance = {}
 doGetCoeff(obs_reco, obs_gen, obs_name, obs_bins, 'fullNNLOPS')
+

config/skim_MC_tree.C

@@ -239,6 +239,7 @@ void add(TString input_dir, TString year, TString prod_mode, bool t_failed=true)
   float _ZZy,ZZPt,ZZEta;
   Short_t nCleanedJetsPt30,nCleanedJetsPt30_jesUp,nCleanedJetsPt30_jesDn;
   Short_t _njets_pt30_eta2p5,_njets_pt30_eta2p5_jesup,_njets_pt30_eta2p5_jesdn;
+  Float_t _pTj1, _GENpTj1;
   vector<float> *LepPt = 0;
   vector<float> *LepPhi = 0;
   vector<float> *LepEta = 0;
@@ -261,6 +262,9 @@ void add(TString input_dir, TString year, TString prod_mode, bool t_failed=true)
   TBranch *njets_pt30_eta2p5_jesup = T->Branch("njets_pt30_eta2p5_jesup",&_njets_pt30_eta2p5_jesup,"_njets_pt30_eta2p5_jesup/S");
   TBranch *njets_pt30_eta2p5_jesdn = T->Branch("njets_pt30_eta2p5_jesdn",&_njets_pt30_eta2p5_jesdn,"_njets_pt30_eta2p5_jesdn/S");
 
+  TBranch *pTj1 = T->Branch("pTj1",&_pTj1,"pTj1/F");
+  TBranch *GENpTj1 = T->Branch("GENpTj1",&_GENpTj1,"GENpTj1/F");
+
   if (!t_failed) {
     T->SetBranchAddress("ZZPt",&ZZPt);
     T->SetBranchAddress("ZZEta",&ZZEta);
@@ -326,6 +330,15 @@ void add(TString input_dir, TString year, TString prod_mode, bool t_failed=true)
       }
     }
 
+    // leading GENjet pT
+    _GENpTj1 = 0;
+    for (unsigned int i = 0; i < GenCleanedJetPt->size(); ++i)
+    {
+      if(GenCleanedJetPt->at(i) > 30. && abs(GenCleanedJetEta->at(i)) < 2.5 && GenCleanedJetPt->at(i) > _GENpTj1) {
+        _GENpTj1 = GenCleanedJetPt->at(i);
+      }
+    }
+    GENpTj1->Fill();
     GenLepPtSorted->Fill();
     GenLepEtaSorted->Fill();
     GenLepPhiSorted->Fill();
@@ -343,6 +356,15 @@ void add(TString input_dir, TString year, TString prod_mode, bool t_failed=true)
 
     if (t_failed) continue; // From now on reco-only variables
 
+    // leading jet pT
+    _pTj1 = 0;
+    for (unsigned int i = 0; i < JetPt->size(); ++i)
+    {
+      if(JetPt->at(i)>30 && abs(JetEta->at(i))<2.5 && JetPt->at(i) > _pTj1) {
+        _pTj1 = JetPt->at(i);
+      }
+    }
+
     // njets
     _njets_pt30_eta2p5 = 0;
     for(unsigned int i=0;i<JetPt->size();i++){
@@ -364,6 +386,7 @@ void add(TString input_dir, TString year, TString prod_mode, bool t_failed=true)
         _njets_pt30_eta2p5_jesup++;
       }
     }
+    pTj1->Fill();
     njets_pt30_eta2p5->Fill();
     njets_pt30_eta2p5_jesdn->Fill();
     njets_pt30_eta2p5_jesup->Fill();
@@ -736,3 +759,4 @@ void skim_MC_tree (TString prod_mode = "ZH125", TString year = "2017"){
 
   return 0;
 }
+

config/skim_data_tree.C

@@ -57,13 +57,14 @@ void skim_data_tree (int year = 2017){
   Short_t Z1Flav,Z2Flav;
   float ZZMass, ZZPt, ZZEta;
   float _chan, _CMS_zz4l_mass, _ZZy;
-  float _njets_pt30_eta2p5;
+  float _njets_pt30_eta2p5, _pTj1;
   vector<float> *JetPt = 0;
   vector<float> *JetEta = 0;
   TBranch *chan = T->Branch("chan",&_chan,"chan/F");
   TBranch *ZZy = T->Branch("ZZy",&_ZZy,"ZZy/F");
   TBranch *CMS_zz4l_mass = T->Branch("CMS_zz4l_mass",&_CMS_zz4l_mass,"CMS_zz4l_mass/F");
   TBranch *njets_pt30_eta2p5 = T->Branch("njets_pt30_eta2p5",&_njets_pt30_eta2p5,"njets_pt30_eta2p5/F");
+  TBranch *pTj1 = T->Branch("pTj1",&_pTj1,"pTj1/F");
   T->SetBranchAddress("Z1Flav",&Z1Flav);
   T->SetBranchAddress("Z2Flav",&Z2Flav);
   T->SetBranchAddress("ZZMass",&ZZMass);
@@ -95,7 +96,15 @@ void skim_data_tree (int year = 2017){
       }
     }
 
+    _pTj1 = 0;
+    for(unsigned int i=0;i<JetPt->size();i++){
+      if (JetPt->at(i)>30 && abs(JetEta->at(i))<2.5 && JetPt->at(i) > _pTj1){
+        _pTj1 = JetPt->at(i);
+      }
+    }
+
     njets_pt30_eta2p5->Fill();
+    pTj1->Fill();
     ZZy->Fill();
     chan->Fill();
     CMS_zz4l_mass->Fill();
@@ -106,3 +115,4 @@ void skim_data_tree (int year = 2017){
 
   return 0;
 }
+

env.sh

@@ -3,7 +3,7 @@ echo $base_path;
 cd LHScans; mkdir plots;
 cd $base_path; 
 mkdir combine_files;
-mkdir datacard; cd datacard; mkdir datacard_2016; mkdir datacad_2017; mkdir datacard_2018;
+mkdir datacard; cd datacard; mkdir datacard_2016; mkdir datacard_2017; mkdir datacard_2018;
 cd $base_path;
 mkdir plots;
 cd templates; mkdir 2016; mkdir 2017; mkdir 2018;

fit/RunFiducialXS.py

@@ -87,7 +87,7 @@ def produceDatacards(obsName, observableBins, ModelName, PhysicalModel):
     print '[Producing workspace/datacards for obsName '+obsName+', bins '+str(observableBins)+']'
     fStates = ['2e2mu','4mu','4e']
     nBins = len(observableBins)
-    if 'jet' in obsName: JES = True
+    if (('jet' in obsName) | (obsName == 'pTj1')): JES = True
     else: JES = False
     for year in years:
         os.chdir('../datacard/datacard_'+year)

fit/createXSworkspace.py

@@ -68,6 +68,7 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, usecfacto
             if obsName == "massZ2": obsName_help = "Z2Mass"
             if obsName == "massZ1": obsName_help = "Z1Mass"
             if obsName == "njets_pt30_eta2p5": obsName_help = "njets_pt30_eta2p5"
+            if obsName == "pTj1": obsName_help = "pTj1"
             observable = ROOT.RooRealVar(obsName_help,obsName_help,float(obs_bin_lowest),float(obs_bin_highest))
         observable.Print()
 

templates/RunTemplates.py

@@ -106,13 +106,13 @@ def add_njets(pt,eta):
         if pt[i]>30 and abs(eta[i])<2.5: n=n+1
     return n
 def add_leadjet(pt,eta):
-    pT = []
-    for i in range(len(pt)):
-        if pt[i]>30 and abs(eta[i])<2.5: pT.append(pt)
-    if len(pT)!=0:
-        return np.max(pT)
+    _pTj1 = 0.0
+    if len(pt) == 0:
+	return _pTj1
     else:
-        return -1
+        for i in range(len(pt)):
+            if (pt[i]>30 and abs(eta[i])<2.5 and pt[i] > _pTj1): _pTj1 = pt[i]
+        return _pTj1 
 
 # Rapidity
 def rapidity(p, eta):
@@ -349,7 +349,7 @@ def doTemplates(df_irr, df_red, binning, var, var_string):
                     w = df['weight'].to_numpy()
                     w = np.asarray(w).astype('float')
                     # ------
-		    if (obs_name == 'rapidity4l'):
+		    if(obs_name == 'rapidity4l'):
 			histo = ROOT.TH1D("m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), "m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), 20, 105, 140)
 		    else:
                     	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, 105, 140)
@@ -391,13 +391,13 @@ def doTemplates(df_irr, df_red, binning, var, var_string):
                 w = df['yield_SR'].to_numpy()
                 w = np.asarray(w).astype('float')
                 # ------
-                if (obs_name == 'rapidity4l'):
+                if(obs_name == 'rapidity4l'):
 			histo = ROOT.TH1D("m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), "m4l_"+var_string+"_"+str(bin_low)+"_"+str(bin_high), 20, 105, 140)
                 else:
                 	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, 105, 140)
                 histo.FillN(len(mass4l), mass4l, w)
                 smoothAndNormaliseTemplate(histo, 1)
-		if (obs_name == 'rapidity4l'):
+                if(obs_name == 'rapidity4l'):
 			outFile = ROOT.TFile.Open(str(year)+"/"+var_string+"/XSBackground_ZJetsCR_"+f+"_"+var_string+"_"+str(bin_low)+"_"+str(bin_high)+".root", "RECREATE")
 		else:
                 	outFile = ROOT.TFile.Open(str(year)+"/"+var_string+"/XSBackground_ZJetsCR_"+f+"_"+var_string+"_"+str(int(bin_low))+"_"+str(int(bin_high))+".root", "RECREATE")
@@ -433,6 +433,7 @@ def doTemplates(df_irr, df_red, binning, var, var_string):
 elif(obs_name == 'massZ1'): obs_reco = 'Z1Mass'
 elif(obs_name == 'massZ2'): obs_reco = 'Z2Mass'
 elif(obs_name == 'njets_pt30_eta2p5'): obs_reco = 'njets_pt30_eta2p5'
+elif(obs_name == 'pTj1'): obs_reco = 'pTj1'
 
 # Generate pandas for ggZZ and qqZZ
 d_bkg = {}

templates/plot_templates.cpp

@@ -145,13 +145,20 @@ int main (int argc, char ** argv){
   // TString binRangeLow[N_BINS]  = {"12", "24", "28", "32", "40", "48"};
   // TString binRangeHigh[N_BINS] = {"24", "28", "32", "40", "48", "65"};
   // TString binRangeLeg[N_BINS]  = {"12 <  m(Z2) < 24 GeV", "24 < m(Z2) < 28 GeV", "28 < m(Z2) < 32 GeV", "32 < m(Z2) < 40 GeV", "40 < m(Z2) < 48 GeV", "48 < m(Z2) < 65 GeV"};
-
-  TString obsTag = "njets_pt30_eta2p5";
-  const int N_BINS = 5;
-  TString binRange[N_BINS]     = {"0_1", "1_2", "2_3", "3_4", "4_20"};
-  TString binRangeLow[N_BINS]  = {"0", "1", "2", "3", "4"};
-  TString binRangeHigh[N_BINS] = {"1", "2", "3", "4", "20"};
-  TString binRangeLeg[N_BINS]  = {"nJet=0", "nJet=1", "nJet=2", "nJet=3", "nJet>=4"};
+
+  // TString obsTag = "njets_pt30_eta2p5";
+  // const int N_BINS = 5;
+  // TString binRange[N_BINS]     = {"0_1", "1_2", "2_3", "3_4", "4_20"};
+  // TString binRangeLow[N_BINS]  = {"0", "1", "2", "3", "4"};
+  // TString binRangeHigh[N_BINS] = {"1", "2", "3", "4", "20"};
+  // TString binRangeLeg[N_BINS]  = {"nJet=0", "nJet=1", "nJet=2", "nJet=3", "nJet>=4"};
+
+  TString obsTag = "pTj1";
+  const int N_BINS = 4;
+  TString binRange[N_BINS]     = {"30_55", "55_95", "95_200", "200_13000"};
+  TString binRangeLow[N_BINS]  = {"30", "55", "95", "200"};
+  TString binRangeHigh[N_BINS] = {"55", "95", "200", "13000"};
+  TString binRangeLeg[N_BINS]  = {"30 < pTj1 < 55 GeV", "55 < pTj1 < 95 GeV", "95 < pTj1 < 200 GeV", "200 < pTj1 < 1300 GeV"};
 
   const int N_BKGS = 3;
   TString bkgName[N_BKGS]   = {"qqzz", "ggzz", "ZJetsCR"};