include first version of the 2TagCount code

test/ttbar/TTbarEventAnalysis.h

@@ -3,6 +3,7 @@
 
 #include "TFile.h"
 #include "TH1F.h"
+#include "TH2F.h"
 #include "TString.h"
 #include "TTree.h"
 #include "TSystem.h"
@@ -13,6 +14,7 @@
 #include <iostream>
 #include <map>
 #include <vector>
+#include <algorithm>
 
 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
 
@@ -24,6 +26,27 @@ struct LJKinematics_t
 
 bool sortLJKinematicsByDR (LJKinematics_t i,LJKinematics_t j) { return (i.dr<j.dr); }
 
+//prepare to read the tree (for jets only interested in a couple of variables)
+class MyEventInfoBranches_t 
+{
+  public:
+    Int_t Run,Evt,LumiBlock,nPV;
+    Int_t   ttbar_chan, ttbar_trigWord, ttbar_metfilterWord;
+    Int_t   ttbar_nl, ttbar_lid[10], ttbar_lgid[10], ttbar_lch[10];
+    Float_t ttbar_lpt[10], ttbar_leta[10], ttbar_lphi[10], ttbar_lm[10];
+    Float_t ttbar_metpt,ttbar_metphi;
+    Float_t ttbar_rho;
+    Int_t   ttbar_nw;
+    Float_t nPUtrue;
+    Int_t nPU;
+    Float_t ttbar_w[1095];
+    Int_t nJet;
+    Float_t Jet_pt[100],Jet_genpt[100],Jet_area[100],Jet_jes[100],Jet_eta[100],Jet_phi[100],Jet_mass[100];
+    Float_t Jet_Svx[100],Jet_CombIVF[100],Jet_Proba[100],Jet_Ip2P[100],Jet_DeepCSVBDisc[100],Jet_DeepFlavourBDisc[100];
+    Int_t Jet_nseltracks[100];
+    Int_t Jet_flavour[100];
+};
+
 class TTbarEventAnalysis
 {
  public:
@@ -33,23 +56,14 @@ class TTbarEventAnalysis
     puWgtGr_(0),puWgtDownGr_(0),puWgtUpGr_(0)
       {
 	//jet uncertainty parameterization
-	TString jecUncUrl("${CMSSW_BASE}/src/RecoBTag/PerformanceMeasurements/test/ttbar/data/Summer15_25nsV5_DATA_Uncertainty_AK4PFchs.txt");
+	TString jecUncUrl("${CMSSW_BASE}/src/RecoBTag/PerformanceMeasurements/test/ttbar/data/Autumn18_V8_MC_Uncertainty_AK4PF.txt");
 	gSystem->ExpandPathName(jecUncUrl);
 	jecUnc_ = new JetCorrectionUncertainty(jecUncUrl.Data());
 
 	//pileup weights
-	TString puWgtUrl("${CMSSW_BASE}/src/RecoBTag/PerformanceMeasurements/test/ttbar/data/pileupWgts.root");
-	gSystem->ExpandPathName(puWgtUrl);
-	TFile *fIn=TFile::Open(puWgtUrl);
-	if(fIn){
-	  puWgtGr_     = (TGraph *)fIn->Get("puwgts_nom");
-	  puWgtDownGr_ = (TGraph *)fIn->Get("puwgts_down");
-	  puWgtUpGr_   = (TGraph *)fIn->Get("puwgts_up");
-	  fIn->Close();
-	}
-	else{
-	  std::cout << "Unable to find data/pileupWgts.root, no PU reweighting will be applied" << std::endl;
-	}
+	TString stdTarget("${CMSSW_BASE}/src/RecoBTag/PerformanceMeasurements/test/ttbar/data/pileupWgts.root");
+	gSystem->ExpandPathName(stdTarget);
+	SetPUWeightTarget(stdTarget,"");
       }
   ~TTbarEventAnalysis(){}
   void setReadTTJetsGenWeights(bool readTTJetsGenWeights)     { readTTJetsGenWeights_=readTTJetsGenWeights; }
@@ -59,6 +73,28 @@ class TTbarEventAnalysis
   void prepareOutput(TString outFile);
   void processFile(TString inFile,TH1F *xsecWgt,Bool_t isData);
   void finalizeOutput();
+	void SetPUWeightTarget(TString targetFile,TString sampleName){
+      TFile *fIn=TFile::Open(targetFile);
+	    if(fIn){
+        std::string nom( "puwgts_nom");  
+        std::string up(  "puwgts_down");
+        std::string down("puwgts_up");
+        if(!sampleName.IsNull()){
+            nom.append("_");
+            nom.append(sampleName);
+            up.append("_");
+            up.append(sampleName);
+            down.append("_");
+            down.append(sampleName);
+        }
+	      puWgtGr_     = (TGraph *)fIn->Get(nom.c_str());
+	      puWgtDownGr_ = (TGraph *)fIn->Get(down.c_str());
+	      puWgtUpGr_   = (TGraph *)fIn->Get(up.c_str());
+	    }else{
+	      std::cout << "Unable to find data/pileupWgts.root, no PU reweighting will be applied" << std::endl;
+	    }
+	    fIn->Close();
+	}
 
  private:
   JetCorrectionUncertainty *jecUnc_;
@@ -67,14 +103,27 @@ class TTbarEventAnalysis
   std::vector<float> getJetEnergyScales(float pt,float eta,float rawsf,float area,float rho);
   std::vector<float> getJetResolutionScales(float pt, float eta, float genjpt);
 
+  MyEventInfoBranches_t ev;
+  std::vector<Int_t> selJets;
+  Float_t evWgt;
+
+  std::vector<std::string> wpLabel;
+  std::vector<std::string> systName;
+  std::map<std::string,double> systWeight;
+
+  void TwoTag(std::string tagName, std::string discriminator, std::pair<int, int>);  
+  std::map<std::string,std::vector<float>> btaggingWPs;
+  void GetBestJetPair(std::pair<int, int>& myIndices, std::string discriminator="deepCSV");
+  float ReturnVarAtIndex(std::string varName, unsigned int index);
+
   TGraph *puWgtGr_,*puWgtDownGr_,*puWgtUpGr_;
   bool readTTJetsGenWeights_;
   TString weightsDir_;
   std::vector<TString> tmvaVarNames_;
   TMVA::Reader *tmvaReader_;
   TFile *outF_;
   Int_t eventInfo_[3],ttbar_chan_,npv_;
-  Float_t weight_[15];
+  Float_t weight_[27];
   Int_t jetFlavour_[2],jetmult_,jetrank_;
   Float_t jetPt_[2],jetEta_[2];
   Float_t close_mlj_[5],close_deta_,close_dphi_,close_ptrel_,close_lj2ll_deta_, close_lj2ll_dphi_;
@@ -85,6 +134,8 @@ class TTbarEventAnalysis
   std::vector<std::pair<Int_t,Int_t> > triggerBits_;
   TTree *kinTree_,*ftmTree_;
   std::map<TString,TH1F *> histos_;
+  std::map<TString,TH2F *> histos2d_;
+  bool noEventsSelected;
 
 };
 

test/ttbar/data/samples_Run2018.json

@@ -0,0 +1,20 @@
+{
+"MC13TeV_TTJets_AH_training" : [380.09,   0, "/TTToHadronic_TuneCP5_13TeV-powheg-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v1/MINIAODSIM", "t#bar{t} AH",          920, false, false, 1],
+"MC13TeV_TTJets_DL_training" : [87.31,   0, "/TTTo2L2Nu_TuneCP5_13TeV-powheg-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v1/MINIAODSIM", "t#bar{t} DL",          920, false, false, 1],
+"MC13TeV_TTJets_SL_training" : [364.35,   0, "/TTToSemiLeptonic_TuneCP5_13TeV-powheg-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v1/MINIAODSIM", "t#bar{t} SL",          920, false, false, 1],
+"MC13TeV_atW"    : [34.97,     0, "/ST_tW_antitop_5f_inclusiveDecays_TuneCP5_13TeV-powheg-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15_ext1-v1/MINIAODSIM", "tW",       7, false, false, 1],
+"MC13TeV_tW"     : [34.97,     0, "/ST_tW_top_5f_inclusiveDecays_TuneCP5_13TeV-powheg-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15_ext1-v1/MINIAODSIM", "tW",        7, false, false, 1],
+"MC13TeV_DYJetsToLL_M-50-madgraphMLM" : [6225.42,   0, "/DYJetsToLL_M-50_TuneCP5_13TeV-madgraphMLM-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v1/MINIAODSIM", "DY",          8, false, false, 1],
+"MC13TeV_DYJetsToLL-M-50-amcatnloFXFX" : [6529.0,   0, "/DYJetsToLL_M-50_TuneCP5_13TeV-amcatnloFXFX-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v1/MINIAODSIM", "DY",          8, false, false, 1],
+"MC13TeV_DYJetsToLL_M-10to50-madgraphMLM" : [15820,   0, "/DYJetsToLL_M-10to50_TuneCP5_13TeV-madgraphMLM-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v2/MINIAODSIM", "DY",          8, false, false, 1],
+"MC13TeV_WW" : [75.8,   0, "/WW_TuneCP5_13TeV-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v2/MINIAODSIM", "WW",          4, false, false, 1],
+"MC13TeV_WZ" : [27.6,  0, "/WZ_TuneCP5_13TeV-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v3/MINIAODSIM", "WZ",          4, false, false, 1],
+"MC13TeV_ZZ" : [12.14,  0, "/ZZ_TuneCP5_13TeV-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v2/MINIAODSIM", "ZZ",          4, false, false, 1],
+"MC13TeV_W2JetsToLNu" : [2793.0, 0, "/W2JetsToLNu_TuneCP5_13TeV-madgraphMLM-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v2/MINIAODSIM", "other", 4],
+"MC13TeV_W3JetsToLNu" : [992.5, 0, "/W3JetsToLNu_TuneCP5_13TeV-madgraphMLM-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v2/MINIAODSIM", "other", 4],
+"MC13TeV_W4JetsToLNu" : [544.3, 0, "/W4JetsToLNu_TuneCP5_13TeV-madgraphMLM-pythia8/RunIIAutumn18MiniAOD-102X_upgrade2018_realistic_v15-v2/MINIAODSIM", "other", 4],
+"crab_Data13TeV_MuonEG_Run2018A": [1,      1, "/MuonEG/Run2018A-PromptReco/MINIAOD",       "Data", 1, false, false, 0],
+"crab_Data13TeV_MuonEG_Run2018B": [1,      1, "/MuonEG/Run2018B-PromptReco/MINIAOD",       "Data", 1, false, false, 0],
+"crab_Data13TeV_MuonEG_Run2018C": [1,      1, "/MuonEG/Run2018C-PromptReco/MINIAOD",       "Data", 1, false, false, 0],
+"crab_Data13TeV_MuonEG_Run2018D": [1,      1, "/MuonEG/Run2018D-PromptReco/MINIAOD",       "Data", 1, false, false, 0]
+}

test/ttbar/plotter.py

@@ -23,6 +23,7 @@ def __init__(self,name):
         self.plotformats = ['pdf','png']
         self.savelog = False
         self.ratiorange = (0.46,1.54)
+        #self.ratiorange = (0.75,1.25)
 
     def add(self, h, title, color, isData,isSyst):
         h.SetTitle(title)
@@ -123,7 +124,7 @@ def show(self, outDir,lumi,noScale=False,saveTeX=False):
         p1.cd()
 
         # legend
-        leg = ROOT.TLegend(0.5, 0.85-0.03*max(len(self.mc)-2,0), 0.98, 0.9)        
+        leg = ROOT.TLegend(0.5, 0.80-0.03*max(len(self.mc)-2,0), 0.98, 0.9)        
         leg.SetBorderSize(0)
         leg.SetFillStyle(0)
         leg.SetTextFont(43)
@@ -185,6 +186,12 @@ def show(self, outDir,lumi,noScale=False,saveTeX=False):
         if self.dataH:
             if maxY<self.dataH.GetMaximum():
                 maxY=self.dataH.GetMaximum()
+        if totalMC and self.dataH:
+            print self.name,"data, MC, data/MC",self.dataH.Integral(),totalMC.Integral(),
+            if totalMC.Integral()>0:
+                print self.dataH.Integral()/totalMC.Integral(),
+            print
+
         frame.GetYaxis().SetRangeUser(0.1,maxY*1.3)
         frame.SetDirectory(0)
         frame.Reset('ICE')
@@ -230,7 +237,26 @@ def show(self, outDir,lumi,noScale=False,saveTeX=False):
         p2.cd()
         ratioframe=frame.Clone('ratioframe')
         ratioframe.GetYaxis().SetTitle('Data/MC')
-        ratioframe.GetYaxis().SetRangeUser(self.ratiorange[0], self.ratiorange[1])
+        maxVal=0
+        minVal=10
+        dataOverMC=self.dataH.Clone("dataOverMC")
+        print "integrals",dataOverMC.GetName(),totalMC.GetName(),dataOverMC.Integral(),totalMC.Integral()
+        dataOverMC.Divide(totalMC)
+        for xBin in xrange(1,dataOverMC.GetXaxis().GetNbins()+1):
+            thisValue=dataOverMC.GetBinContent(xBin)
+            if thisValue==0:
+                continue
+            #print xBin,thisValue
+            if maxVal<thisValue:
+                maxVal=thisValue
+            if minVal>thisValue:
+                minVal=thisValue
+
+        print "maxVal,minVal,",maxVal,minVal
+        print "range",max(self.ratiorange[0],minVal*0.95), min(self.ratiorange[1],maxVal*1.05)
+
+        ratioframe.GetYaxis().SetRangeUser(max(self.ratiorange[0],minVal*0.95), min(self.ratiorange[1],maxVal*1.05))
+        #ratioframe.GetYaxis().SetRangeUser(self.ratiorange[0], self.ratiorange[1])
         self._garbageList.append(frame)
         ratioframe.GetYaxis().SetNdivisions(5)
         ratioframe.GetYaxis().SetLabelSize(0.18)        
@@ -371,7 +397,8 @@ def convertToPoissonErrorGr(h):
     #check https://twiki.cern.ch/twiki/bin/view/CMS/PoissonErrorBars
     alpha = 1 - 0.6827;
     grpois = ROOT.TGraphAsymmErrors(h);
-    for i in xrange(0,grpois.GetN()+1) :
+    #for i in xrange(0,grpois.GetN()+1) :
+    for i in xrange(0,grpois.GetN()) :
         N = grpois.GetY()[i]
         if N<200 :
             L = 0
@@ -402,6 +429,7 @@ def main():
     parser.add_option(      '--rebin',       dest='rebin',       help='rebin factor',                   default=1,       type=int)
     parser.add_option('-l', '--lumi',        dest='lumi' ,       help='lumi to print out',              default=41.6,    type=float)
     parser.add_option(      '--only',        dest='only',        help='plot only these (csv)',          default='',      type='string')
+    parser.add_option(      '--outLabel',    dest='outLabel',    help='appends the plots dir name',     default='',      type='string')
     (opt, args) = parser.parse_args()
 
     #read list of samples
@@ -422,11 +450,14 @@ def main():
         if slist is None : continue
         for tag,sample in slist: 
 
-            if isSyst and not 't#bar{t}' in sample[3] : continue
+            if isSyst and not 't#bar{t}' in sample[3] : 
+                print "Skipping syst sample",sample
+                continue
 
             inDir=opt.inDir
-            if isSyst : inDir += '/syst'
+            #if isSyst : inDir += '/syst'
             fIn=ROOT.TFile.Open('%s/%s.root' % ( inDir, tag) )
+            print "opening ",'%s/%s.root' % ( inDir, tag) 
             try:
                 for tkey in fIn.GetListOfKeys():
                     key=tkey.GetName()
@@ -449,6 +480,8 @@ def main():
     ROOT.gStyle.SetOptStat(0)
     ROOT.gROOT.SetBatch(True)
     outDir=opt.inDir+'/plots'
+    if opt.outLabel:
+        outDir=outDir+"_"+opt.outLabel
     os.system('mkdir -p %s' % outDir)
     for p in plots : 
         if opt.saveLog    : plots[p].savelog=True

test/ttbar/runPileupEstimation.py

@@ -3,7 +3,7 @@
 import json
 import commands
 import ROOT
-from SimGeneral.MixingModule.mix_2016_25ns_SpringMC_PUScenarioV1_PoissonOOTPU_cfi import *
+from SimGeneral.MixingModule.mix_2018_25ns_JuneProjectionFull18_PoissonOOTPU_cfi import *
 
 """
 steer the script
@@ -30,11 +30,14 @@ def main():
     #compute pileup in data assuming different xsec
     puDist=[]
     puWgts=[]
-    MINBIASXSEC={'nom':opt.mbXsec,'up':opt.mbXsec*1.1,'down':opt.mbXsec*0.9}
+    XSECERROR=1.046
+    MINBIASXSEC={'nom':opt.mbXsec,'up':opt.mbXsec*XSECERROR,'down':opt.mbXsec/XSECERROR}
     for scenario in MINBIASXSEC:
         print scenario, 'xsec=',MINBIASXSEC[scenario]
         cmd='pileupCalc.py -i %s --inputLumiJSON %s --calcMode true --minBiasXsec %f --maxPileupBin %d --numPileupBins %s Pileup.root'%(opt.inJson,opt.puJson,MINBIASXSEC[scenario],NPUBINS,NPUBINS)
         commands.getstatusoutput(cmd)
+        print cmd
+
 
         fIn=ROOT.TFile.Open('Pileup.root')
         pileupH=fIn.Get('pileup')