diff --git a/TopAnalysis/bin/analysisWrapper.cc b/TopAnalysis/bin/analysisWrapper.cc index 6bf2fbf102c84d9b929c9f5bbca186e5fd0f0fc8..7290f734039de0f9259639bb783e4ef7c0cb0b31 100644 --- a/TopAnalysis/bin/analysisWrapper.cc +++ b/TopAnalysis/bin/analysisWrapper.cc @@ -22,7 +22,8 @@ void printHelp() << "\t --systVar - specify single systematic variation" << endl << "\t --era - era directory to use for corrections, uncertainties" << endl << "\t --normTag - normalization tag" << endl - << "\t --method - method to run" << endl; + << "\t --method - method to run" << endl + << "\t --mvatree - store selected events in a tree for mva" << endl; } // @@ -32,7 +33,7 @@ int main(int argc, char* argv[]) TString in(""),out(""),era(""),normTag(""),method(""); std::string systVar(""); bool runSysts(false); - bool debug(false); + bool debug(false), skimtree(false); int channel(0),charge(0),flag(0); for(int i=1;i<argc;i++){ string arg(argv[i]); @@ -45,6 +46,7 @@ int main(int argc, char* argv[]) else if(arg.find("--in")!=string::npos && i+1<argc) { in=argv[i+1]; i++;} else if(arg.find("--out")!=string::npos && i+1<argc) { out=argv[i+1]; i++;} else if(arg.find("--debug")!=string::npos) { debug=true; } + else if(arg.find("--mvatree")!=string::npos) { skimtree=true; } else if(arg.find("--normTag")!=string::npos && i+1<argc) { normTag=argv[i+1]; i++;} else if(arg.find("--era")!=string::npos && i+1<argc) { era=argv[i+1]; i++;} else if(arg.find("--method")!=string::npos && i+1<argc) { method=argv[i+1]; i++;} @@ -85,7 +87,10 @@ int main(int argc, char* argv[]) //check method to run if(method=="ExclusiveTop::RunExclusiveTop") RunExclusiveTop(in,out,channel,charge,normH,puH,era,debug); - else if(method=="VBFVectorBoson::RunVBFVectorBoson") RunVBFVectorBoson(in,out,flag,normH,puH,era,debug); + else if(method=="VBFVectorBoson::RunVBFVectorBoson") { + VBFVectorBoson myVBF(in,out,flag,normH,puH,era,debug,skimtree); + myVBF.RunVBFVectorBoson(); + } else { cout << "Check method=" << method <<endl; printHelp(); diff --git a/TopAnalysis/data/era2017/pileupWgts.root b/TopAnalysis/data/era2017/pileupWgts.root index c854232f8adae7d85d77d0a5bf9eebb983f83737..d2c0f6a4f67a93c12b19b49e06e49e573a09348d 100644 Binary files a/TopAnalysis/data/era2017/pileupWgts.root and b/TopAnalysis/data/era2017/pileupWgts.root differ diff --git a/TopAnalysis/interface/SelectionTools.h b/TopAnalysis/interface/SelectionTools.h index 6ca12a59378edc3b146293b70baefb94abb03e00..00d32028b52e5cd24b0c8fe0aeee1bdc3c3ec32b 100644 --- a/TopAnalysis/interface/SelectionTools.h +++ b/TopAnalysis/interface/SelectionTools.h @@ -37,7 +37,7 @@ class SelectionTool { std::vector<Particle> &getSelLeptons() { return leptons_; } std::vector<Particle> &getVetoLeptons() { return vetoLeptons_; } std::vector<Particle> flaggedPhotons(MiniEvent_t &ev); - std::vector<Particle> selPhotons(std::vector<Particle> &flaggedPhotons,int qualBit=LOOSE, double minPt = 30., double maxEta = 2.5, std::vector<Particle> veto = {}); + std::vector<Particle> selPhotons(std::vector<Particle> &flaggedPhotons,int qualBit=LOOSE, std::vector<Particle> leptons = {}, double minPt = 30., double maxEta = 2.5, std::vector<Particle> veto = {}); std::vector<Particle> &getSelPhotons() { return photons_; } std::vector<Jet> getGoodJets(MiniEvent_t &ev, double minPt = 30., double maxEta = 4.7, std::vector<Particle> leptons = {},std::vector<Particle> photons = {}); // changed for the moment to VBF std::vector<Jet> &getJets() { return jets_; } diff --git a/TopAnalysis/interface/VBFVectorBoson.h b/TopAnalysis/interface/VBFVectorBoson.h index 75dbad926f694534dbcf157840176bd9f2f0499b..53fecd556806e5378e6ccc07ba8de88a026f0447 100644 --- a/TopAnalysis/interface/VBFVectorBoson.h +++ b/TopAnalysis/interface/VBFVectorBoson.h @@ -4,15 +4,194 @@ #include "TLorentzVector.h" #include "TopLJets2015/TopAnalysis/interface/ObjectTools.h" #include "TopLJets2015/TopAnalysis/interface/SelectionTools.h" +#include <TFile.h> +#include <TROOT.h> +#include <TH1.h> +#include <TH2.h> +#include <TSystem.h> +#include <TGraph.h> +#include <TGraphAsymmErrors.h> + +#include "TopLJets2015/TopAnalysis/interface/MiniEvent.h" +#include "TopLJets2015/TopAnalysis/interface/CommonTools.h" +#include "TopLJets2015/TopAnalysis/interface/VBFVectorBoson.h" +#include "TopLJets2015/TopAnalysis/interface/EfficiencyScaleFactorsWrapper.h" + +#include "PhysicsTools/CandUtils/interface/EventShapeVariables.h" + +#include <vector> +#include <set> +#include <iostream> +#include <algorithm> +#include "TRandom3.h" +#include "TMath.h" +using namespace std; //Vector boson will be either Z or photon at the moment -void RunVBFVectorBoson(TString filename, - TString outname, - Int_t anFlag, - TH1F *normH, - TH1F *genPU, - TString era, - Bool_t debug=false); +struct Category{ + float MM,A,VBF,HighPt,HighPtVBF,V1J; + Category(float * cat){ + MM = cat[0]; + A = cat[1]; + VBF = cat[2]; + HighPt = cat[3]; + HighPtVBF = cat[4]; + V1J = cat[5]; + }; + Category(){ + MM = 0; + A = 0; + VBF = 0; + HighPt = 0; + HighPtVBF = 0; + V1J = 0; + }; + void set(float * cat){ + MM = cat[0]; + A = cat[1]; + VBF = cat[2]; + HighPt = cat[3]; + HighPtVBF = cat[4]; + V1J = cat[5]; + }; +}; + +class VBFVectorBoson{ +public: + VBFVectorBoson(TString filename_, + TString outname_, + Int_t anFlag_, + TH1F *normH_, + TH1F *genPU_, + TString era_, + Bool_t debug_=false, Bool_t skimtree_=false): + filename(filename_),outname(outname_),anFlag(anFlag_), era(era_), debug(debug_), skimtree(skimtree_) + { + normH = (TH1F*)normH_->Clone("normH_c"); + genPU = (TH1F*)genPU_->Clone("genPu_c"); + fMVATree = NULL; + newTree = NULL; + init(); + setXsecs(); + rnd.SetSeed(123456789); + }; + ~VBFVectorBoson(){} + void init(){ + this->readTree(); + cout << "...producing " << outname << " from " << nentries << " events" << endl; + this->prepareOutput(); + this->bookHistograms(); + this->loadCorrections(); + if(skimtree){ + this->addMVAvars(); + } + selector = new SelectionTool(filename, debug, triggerList,SelectionTool::VBF); + std::cout << "init done" << std::endl; + } + + void saveHistos(); + void readTree(); + void prepareOutput(); + void bookHistograms(); + void setGammaZPtWeights(); + void loadCorrections(); + void addMVAvars(); + void fill(MiniEvent_t ev, TLorentzVector boson, std::vector<Jet> jets, std::vector<double> cplotwgts, TString c); + void RunVBFVectorBoson(); + + + +private: + TString filename, outname, baseName; + Int_t anFlag, nentries; + TH1F * normH, * genPU; + TString era; + TH1* triggerList; + Bool_t debug, skimtree, isQCDEMEnriched; + TFile * f /*inFile*/, *fMVATree, *fOut; + TTree * t /*inTree*/, *newTree /*MVA*/; + HistTool * ht; + MiniEvent_t ev; + float sihih,chiso,r9,hoe, ystar; + double mindrl; + + TRandom3 rnd; + + //LUMINOSITY+PILEUP + LumiTools * lumi; + + //LEPTON EFFICIENCIES + EfficiencyScaleFactorsWrapper * gammaEffWR; + + //JEC/JER + JECTools * jec; + //Photon/Z pt weights + std::map<TString,TGraph *> photonPtWgts; + std::map<TString,std::pair<double,double> > photonPtWgtCtr; + + //Variables to be added to the MVA Tree + float centraleta, forwardeta, jjetas, centjy, ncentjj, dphivj0, dphivj1, dphivj2, dphivj3; + float evtWeight, mjj, detajj , dphijj ,jjpt; + float isotropy, circularity,sphericity, aplanarity, C, D; + float scalarht,balance, mht, training; + + ///////////////////////////////////// + // Categorie for VBF: // + // MM:A:VBF:HighPt:HighPtVBF:V1J // + ///////////////////////////////////// + Category category; + + SelectionTool * selector; + + ///////////////////////////////////////// + // To select events for training // + ///////////////////////////////////////// + + int useForTraining(){ + double myRnd = rnd.Rndm(); + if (myRnd >= 0.5) return 1; + return 0; + } + + ///////////////////////////////////////// + // Quick and ugly cross section getter // + // To be updated such as to use the // + // json file // + ///////////////////////////////////////// + std::map<TString, float> xsecRefs; + void setXsecs(){ + xsecRefs["MC13TeV_TTJets" ] = 832; + xsecRefs["MC13TeV_ZZ" ] = 0.5644; + xsecRefs["MC13TeV_WZ" ] = 47.13; + xsecRefs["MC13TeV_WW" ] = 12.178; + xsecRefs["MC13TeV_SingleTbar_tW" ] = 35.85; + xsecRefs["MC13TeV_SingleT_tW" ] = 35.85; + xsecRefs["MC13TeV_DY50toInf" ] = 5765.4; + xsecRefs["MC13TeV_QCDEM_15to20" ] = 2302200; + xsecRefs["MC13TeV_QCDEM_20to30" ] = 5352960; + xsecRefs["MC13TeV_QCDEM_30to50" ] = 9928000; + xsecRefs["MC13TeV_QCDEM_50to80" ] = 2890800; + xsecRefs["MC13TeV_QCDEM_80to120" ] = 350000; + xsecRefs["MC13TeV_QCDEM_120to170"] = 62964; + xsecRefs["MC13TeV_QCDEM_170to300"] = 18810; + xsecRefs["MC13TeV_QCDEM_300toInf"] = 1350; + xsecRefs["MC13TeV_GJets_HT40to100" ] = 20790 ; + xsecRefs["MC13TeV_GJets_HT100to200"] = 9238; + xsecRefs["MC13TeV_GJets_HT200to400"] = 2305; + xsecRefs["MC13TeV_GJets_HT400to600"] = 274.4; + xsecRefs["MC13TeV_GJets_HT600toInf"] = 93.46; + xsecRefs["MC13TeV_EWKZJJ" ] = 4.32; + } + + float getXsec(){ + for (auto const& x : xsecRefs){ + if (filename.Contains(x.first)) + return x.second; + } + return 1; + } + +}; #endif diff --git a/TopAnalysis/python/Plot.py b/TopAnalysis/python/Plot.py index 02cb734f6d360d7d3b66185acb1bfbd9aa687288..0205a832f0ce37c0109552c983e7e3350c0161c0 100644 --- a/TopAnalysis/python/Plot.py +++ b/TopAnalysis/python/Plot.py @@ -149,7 +149,6 @@ class Plot(object): pass def show(self, outDir,lumi,noStack=False,saveTeX=False,extraText=None,noRatio=False): - if len(self.mc)<1 and self.dataH is None: print '%s has 0 or 1 MC!' % self.name return @@ -184,7 +183,6 @@ class Plot(object): p1.SetTopMargin(0.05) p1.SetBottomMargin(0.12) p1.Draw() - p1.SetGridx(False) p1.SetGridy(False) #True) self._garbageList.append(p1) @@ -213,6 +211,7 @@ class Plot(object): if self.data is None: self.finalize() leg.AddEntry( self.data, self.data.GetTitle(),'p') nlegCols += 1 + for h in self.mc: #compare @@ -224,7 +223,6 @@ class Plot(object): self.mc[h].SetTitle('#splitline{%s}{#chi^{2}=%3.1f (p-val: %3.3f)}'%(self.mc[h].GetTitle(),chi2,pval)) else: refH.SetLineWidth(2) - leg.AddEntry(self.mc[h], self.mc[h].GetTitle(), 'f') nlegCols += 1 if nlegCols ==0 : @@ -238,11 +236,11 @@ class Plot(object): totalMC = None stack = ROOT.THStack('mc','mc') for h in reversed(self.mc): - + print "Integral is %f" %self.mc[h].Integral() if noStack: self.mc[h].SetFillStyle(0) self.mc[h].SetLineColor(self.mc[h].GetFillColor()) - + if lumi == 1 and self.mc[h].Integral() != 1 and self.mc[h].Integral() !=0 : self.mc[h].Scale(1./self.mc[h].Integral()) stack.Add(self.mc[h],'hist') try: @@ -365,7 +363,6 @@ class Plot(object): leg.AddEntry(self.spimpose[m],self.spimpose[m].GetTitle(),'l') if self.data is not None : self.data.Draw('p') - if (totalMC is not None and totalMC.GetMaximumBin() > totalMC.GetNbinsX()/2.): inix = 0.15 leg.SetX1(inix) @@ -400,7 +397,6 @@ class Plot(object): except: pass - #holds the ratio c.cd() if not noRatio and self.dataH and len(self.mc)>0 : @@ -487,7 +483,6 @@ class Plot(object): gr.Draw('p') except: pass - #all done if p1: p1.RedrawAxis() c.cd() @@ -495,10 +490,12 @@ class Plot(object): c.Update() #save + if lumi == 1: frame.GetYaxis().SetRangeUser(0,0.8) for ext in self.plotformats : c.SaveAs(os.path.join(outDir, self.name+'.'+ext)) if self.savelog: p1.cd() frame.GetYaxis().SetRangeUser(1,maxY*50) + if lumi == 1: frame.GetYaxis().SetRangeUser(0.001,1) p1.SetLogy() c.cd() c.Modified() @@ -508,6 +505,8 @@ class Plot(object): if saveTeX : self.convertToTeX(outDir=outDir) + + def convertToTeX(self, outDir): if len(self.mc)==0: print '%s is empty' % self.name diff --git a/TopAnalysis/scripts/plotter.py b/TopAnalysis/scripts/plotter.py index 8bf7ea5f1eadfe1ba0e51091210a66a74c4f39f9..db04a653bfdeb51f778ff63038bc068b7510ec7f 100644 --- a/TopAnalysis/scripts/plotter.py +++ b/TopAnalysis/scripts/plotter.py @@ -6,9 +6,11 @@ from collections import OrderedDict from TopLJets2015.TopAnalysis.Plot import * + """ steer the script """ + def main(): #configuration @@ -28,7 +30,7 @@ def main(): parser.add_option( '--onlyData', dest='onlyData' , help='only plots containing data', default=False, action='store_true') parser.add_option( '--saveTeX', dest='saveTeX' , help='save as tex file as well', default=False, action='store_true') 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=12900, type=float) + parser.add_option('-l', '--lumi', dest='lumi' , help='lumi to print out, if == 1 draw normalized', default=12900, type=float) parser.add_option( '--lumiSpecs', dest='lumiSpecs', help='lumi specifications for some channels [tag:lumi,tag2:lumi2,...]', default=None, type=str) parser.add_option( '--only', dest='only', help='plot only these (csv)', default='', type='string') parser.add_option( '--skip', dest='skip', help='skip these samples (csv)', default='MC13TeV_TTJets_cflip', type='string') @@ -64,6 +66,8 @@ def main(): skipList=opt.skip.split(',') + + #lumi specifications per tag lumiSpecs={} if opt.lumiSpecs: @@ -129,7 +133,6 @@ def main(): for tkey in fIn.GetListOfKeys(): keyIsSyst=False - try: key=tkey.GetName() @@ -140,7 +143,6 @@ def main(): keep=True break if not keep: continue - histos = [] obj=fIn.Get(key) if obj.InheritsFrom('TH2'): @@ -170,7 +172,8 @@ def main(): for hist in histos: if not isData and not '(data)' in sp[1]: - + print hist.GetName() + print hist.GetTitle() #check if a special scale factor needs to be applied sfVal=1.0 for procToScale in procSF: @@ -186,7 +189,6 @@ def main(): if not tag+'_' in key: continue lumi=lumiSpecs[tag] break - hist.Scale(xsec*lumi*puNormSF*sfVal) #rebin if needed @@ -209,7 +211,8 @@ def main(): else: outDir = opt.outDir os.system('mkdir -p %s' % outDir) os.system('rm %s/%s'%(outDir,opt.outName)) - for p in plots : + for p in plots: + print plots[p].name plots[p].mcUnc=opt.mcUnc if opt.saveLog : plots[p].savelog=True skipPlot=False @@ -220,7 +223,8 @@ def main(): if not tag+'_' in p: continue lumi=lumiSpecs[tag] break - if not skipPlot : plots[p].show(outDir=outDir,lumi=lumi,noStack=opt.noStack,saveTeX=opt.saveTeX) + #continue + plots[p].show(outDir=outDir,lumi=lumi,noStack=opt.noStack,saveTeX=opt.saveTeX) plots[p].appendTo('%s/%s'%(outDir,opt.outName)) plots[p].reset() diff --git a/TopAnalysis/scripts/runLocalAnalysis.py b/TopAnalysis/scripts/runLocalAnalysis.py index 058493644d5cae910c18c2e3e27ff00950b34a53..55aa7e949e4a5e78d0331623e69e2c2553f61d80 100755 --- a/TopAnalysis/scripts/runLocalAnalysis.py +++ b/TopAnalysis/scripts/runLocalAnalysis.py @@ -15,18 +15,20 @@ Wrapper to be used when run in parallel """ def RunMethodPacked(args): - method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug=args + method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug,mvatree=args print 'Running ',method,' on ',inF print 'Output file',outF print 'Selection ch=',channel,' charge=',charge,' flag=',flag,' systs=',runSysts print 'Normalization applied from tag=',tag print 'Corrections will be retrieved for era=',era + print 'Make the mva tree? ', mvatree try: cmd='analysisWrapper --era %s --normTag %s --in %s --out %s --method %s --charge %d --channel %d --flag %d --systVar %s'\ %(era, tag, inF, outF, method, charge, channel, flag, systVar) if runSysts : cmd += ' --runSysts' if debug : cmd += ' --debug' + if mvatree : cmd += ' --mvatree' print(cmd) os.system(cmd) @@ -63,6 +65,7 @@ def main(): parser.add_option( '--exactonly', dest='exactonly', help='match only exact sample tags to process [%default]', default=False, action='store_true') parser.add_option( '--outputonly', dest='outputonly', help='filter job submission for a csv list of output files [%default]', default=None, type='string') parser.add_option( '--farmappendix', dest='farmappendix', help='Appendix to condor FARM directory [%default]', default='', type='string') + parser.add_option( '--mvatree', dest='mvatree', help='make mva tree [%default]', default=False, action='store_true') (opt, args) = parser.parse_args() #parse selection lists @@ -126,7 +129,7 @@ def main(): for systVar in varList: outF=opt.output if systVar != 'nominal' and not systVar in opt.output: outF=opt.output[:-5]+'_'+systVar+'.root' - task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,opt.tag,opt.debug) ) + task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,opt.tag,opt.debug, opt.mvatree) ) else: inputTags=getEOSlslist(directory=opt.input,prepend='') @@ -166,7 +169,7 @@ def main(): continue if (len(outputOnlyList) > 1 and not outF in outputOnlyList): continue - task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,tag,opt.debug) ) + task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,tag,opt.debug, opt.mvatree) ) if (opt.skipexisting and nexisting): print '--skipexisting: %s - skipping %d of %d tasks as files already exist'%(systVar,nexisting,len(input_list)) #run the analysis jobs @@ -195,7 +198,7 @@ def main(): condor.write('+JobFlavour = "{0}"\n'.format(opt.queue)) jobNb=0 - for method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug in task_list: + for method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug,mvatree in task_list: jobNb+=1 cfgFile='%s'%(os.path.splitext(os.path.basename(outF))[0]) @@ -216,6 +219,7 @@ def main(): %(inF, localOutF, charge, channel, era, tag, flag, method, systVar) if runSysts : runOpts += ' --runSysts' if debug : runOpts += ' --debug' + if mvatree : runOpts += ' --mvatree' cfg.write('python %s/src/TopLJets2015/TopAnalysis/scripts/runLocalAnalysis.py %s\n'%(cmsswBase,runOpts)) if '/store' in outF: cfg.write('xrdcp ${WORKDIR}/%s root://eoscms//eos/cms/%s\n'%(localOutF,outF)) diff --git a/TopAnalysis/src/MiniEvent.cc b/TopAnalysis/src/MiniEvent.cc index a8fe1f92c7d3a3dab965cae7a805d889455ab1d4..436c1deccf5acd23faae8a37429d35eec118e080 100644 --- a/TopAnalysis/src/MiniEvent.cc +++ b/TopAnalysis/src/MiniEvent.cc @@ -292,6 +292,28 @@ void attachToMiniEventTree(TTree *t,MiniEvent_t &ev,bool full) t->SetBranchAddress("j_id", ev.j_id); t->SetBranchAddress("j_csv", ev.j_csv); t->SetBranchAddress("j_btag", ev.j_btag); + t->SetBranchAddress("j_c1_00", ev.j_c1_00); + t->SetBranchAddress("j_c1_02", ev.j_c1_02); + t->SetBranchAddress("j_c1_05", ev.j_c1_05); + t->SetBranchAddress("j_c1_10", ev.j_c1_10); + t->SetBranchAddress("j_c1_20", ev.j_c1_20); + t->SetBranchAddress("j_c2_00", ev.j_c2_00); + t->SetBranchAddress("j_c2_02", ev.j_c2_02); + t->SetBranchAddress("j_c2_05", ev.j_c2_05); + t->SetBranchAddress("j_c2_10", ev.j_c2_10); + t->SetBranchAddress("j_c2_20", ev.j_c2_20); + t->SetBranchAddress("j_c3_00", ev.j_c3_00); + t->SetBranchAddress("j_c3_02", ev.j_c3_02); + t->SetBranchAddress("j_c3_05", ev.j_c3_05); + t->SetBranchAddress("j_c3_10", ev.j_c3_10); + t->SetBranchAddress("j_c3_20", ev.j_c3_20); + t->SetBranchAddress("j_zg", ev.j_zg); + t->SetBranchAddress("j_mult", ev.j_mult); + t->SetBranchAddress("j_gaptd", ev.j_gaptd); + t->SetBranchAddress("j_gawidth", ev.j_gawidth); + t->SetBranchAddress("j_gathrust", ev.j_gathrust); + t->SetBranchAddress("j_tau32", ev.j_tau32); + t->SetBranchAddress("j_tau21", ev.j_tau21); t->SetBranchAddress("j_deepcsv", ev.j_deepcsv); t->SetBranchAddress("j_vtxpx", ev.j_vtxpx); t->SetBranchAddress("j_vtxpy", ev.j_vtxpy); diff --git a/TopAnalysis/src/SelectionTools.cc b/TopAnalysis/src/SelectionTools.cc index bbceac2a963dba821233665c28374cef31533bfb..af7554cf7078fee4faf784c98ce87534c5ff7243 100644 --- a/TopAnalysis/src/SelectionTools.cc +++ b/TopAnalysis/src/SelectionTools.cc @@ -41,7 +41,7 @@ TString SelectionTool::flagFinalState(MiniEvent_t &ev, std::vector<Particle> pre //decide the channel based on the lepton multiplicity and set lepton collections std::vector<Particle> tightLeptons( selLeptons(preselLeptons,TIGHT) ); - std::vector<Particle> tightPhotons( selPhotons(preselPhotons,TIGHT) ); + std::vector<Particle> tightPhotons( selPhotons(preselPhotons,TIGHT, tightLeptons) ); TString chTag(""); if(anType_==TOP) @@ -70,12 +70,12 @@ TString SelectionTool::flagFinalState(MiniEvent_t &ev, std::vector<Particle> pre if( ch==13*13 && fabs(mll-91)<15 && (tightLeptons[0].pt()>30 || tightLeptons[1].pt()>30)) chTag="MM"; leptons_=tightLeptons; } - if(tightPhotons.size()>=1) { + else if(tightPhotons.size()>=1) { chTag="A"; photons_=tightPhotons; } } - + //select jets based on the leptons and photon candidates float maxJetEta(2.4); if(anType_==VBF) maxJetEta=4.7; @@ -147,6 +147,10 @@ TString SelectionTool::flagFinalState(MiniEvent_t &ev, std::vector<Particle> pre if(chTag=="A") { if(!hasPhotonTrigger) chTag=""; + if((hasEETrigger || hasETrigger) && chTag == "A"){ + cout<< "----------------- This is EE in fact!" <<endl; + chTag = ""; + } if(ev.isData && !isPhotonPD_) chTag=""; } @@ -293,7 +297,7 @@ std::vector<Particle> SelectionTool::flaggedPhotons(MiniEvent_t &ev) } // -std::vector<Particle> SelectionTool::selPhotons(std::vector<Particle> &photons,int qualBit,double minPt, double maxEta,std::vector<Particle> veto){ +std::vector<Particle> SelectionTool::selPhotons(std::vector<Particle> &photons,int qualBit, std::vector<Particle> leptons,double minPt, double maxEta,std::vector<Particle> veto){ std::vector<Particle> selPhotons; for(size_t i =0; i<photons.size(); i++) { @@ -312,7 +316,14 @@ std::vector<Particle> SelectionTool::selPhotons(std::vector<Particle> &photons,i } if(skipThisPhoton) continue; - //lepton is selected + // cross-cleaning with leptos + bool overlapsWithPhysicsObject(false); + for (auto& lepton : leptons) { + if(photons[i].p4().DeltaR(lepton.p4())<0.4) overlapsWithPhysicsObject=true; + } + + if(overlapsWithPhysicsObject) continue; + //photon is selected selPhotons.push_back(photons[i]); } diff --git a/TopAnalysis/src/VBFVectorBoson.cc b/TopAnalysis/src/VBFVectorBoson.cc index d4f665dc4a2c0e7db7e914db3991db16e8d15abc..41c251950729b58c91db7be4fca2e5cc7ea5fa8a 100644 --- a/TopAnalysis/src/VBFVectorBoson.cc +++ b/TopAnalysis/src/VBFVectorBoson.cc @@ -18,7 +18,6 @@ #include <set> #include <iostream> #include <algorithm> - #include "TMath.h" using namespace std; @@ -28,142 +27,37 @@ using namespace std; // // -void RunVBFVectorBoson(TString filename, - TString outname, - Int_t anFlag, - TH1F *normH, - TH1F *genPU, - TString era, - Bool_t debug) +void VBFVectorBoson::RunVBFVectorBoson() { - ///////////////////// - // INITIALIZATION // - /////////////////// - - //PREPARE OUTPUT - TString baseName=gSystem->BaseName(outname); - TString dirName=gSystem->DirName(outname); - TFile *fOut=TFile::Open(dirName+"/"+baseName,"RECREATE"); - fOut->cd(); - - //read photon/Z pT weights if required - std::map<TString,TGraph *> photonPtWgts; - std::map<TString,std::pair<double,double> > photonPtWgtCtr; - if(anFlag>0) { - TString wgtUrl("${CMSSW_BASE}/src/TopLJets2015/TopAnalysis/test/analysis/VBFVectorBoson/raw/plots/ratio_plotter.root"); - gSystem->ExpandPathName(wgtUrl); - TFile *wgtF=TFile::Open(wgtUrl); - if(wgtF) { - cout << "Reading photon/Z pT weights" << endl; - TString pfix(baseName.Contains("Data13TeV_") ? "" : "_mc_MC"); - photonPtWgts["VBFA"] = new TGraph((TH1* )wgtF->Get("VBFA_vectorbosonPt_ratio/VBFA_vectorbosonPt"+pfix)); - photonPtWgts["HighPtA"] = new TGraph((TH1* )wgtF->Get("HighPtA_vectorbosonPt_ratio/HighPtA_vectorbosonPt"+pfix)); - photonPtWgtCtr["VBFA"] = std::pair<double,double>(0.0,0.0); - photonPtWgtCtr["HighPtA"] = std::pair<double,double>(0.0,0.0); - wgtF->Close(); - } else { - cout << "Requested to reweight photon spectrum but could not find " << wgtUrl << endl - << "Proceeding without" << endl; - } - } - - - //READ TREE FROM FILE - MiniEvent_t ev; - TFile *f = TFile::Open(filename); - TH1 *triggerList=(TH1 *)f->Get("analysis/triggerList"); - TTree *t = (TTree*)f->Get("analysis/data"); - attachToMiniEventTree(t,ev,true); - Int_t nentries(t->GetEntriesFast()); - if (debug) nentries = 10000; //restrict number of entries for testing - t->GetEntry(0); - - bool isQCDEMEnriched( filename.Contains("MC13TeV_QCDEM") ); - - cout << "...producing " << outname << " from " << nentries << " events" << endl; - - //LUMINOSITY+PILEUP - LumiTools lumi(era,genPU); - - //LEPTON EFFICIENCIES - EfficiencyScaleFactorsWrapper gammaEffWR(filename.Contains("Data13TeV"),era); - - //B-TAG CALIBRATION - BTagSFUtil btvSF(era,"DeepCSV",BTagEntry::OperatingPoint::OP_MEDIUM,"",0); - - //JEC/JER - JECTools jec(era); - - //BOOK HISTOGRAMS - HistTool ht; - ht.setNsyst(0); - ht.addHist("puwgtctr", new TH1F("puwgtctr", ";Weight sums;Events",2,0,2)); - ht.addHist("qscale", new TH1F("qscale", ";Q^{2} scale;Events",100,0,2000)); - ht.addHist("nvtx", new TH1F("nvtx", ";Vertex multiplicity;Events",100,-0.5,99.5)); - ht.addHist("vpt", new TH1F("vectorbosonPt", ";Boson p_{T}[GeV];Events",25,0,500)); - ht.addHist("vy", new TH1F("vectorbosony", ";Boson rapidity;Events",25,0,3)); - ht.addHist("mindrl", new TH1F("mindrl", ";min #Delta R(boson,lepton);Events",25,0,6)); - ht.addHist("sihih", new TH1F("sihih", ";#sigma(i#eta,i#eta);Events",50,0,0.1)); - ht.addHist("hoe", new TH1F("hoe", ";h/e;Events",25,0,0.1)); - ht.addHist("r9", new TH1F("r9", ";r9;Events",25,0,1.0)); - ht.addHist("chiso", new TH1F("chiso", ";Charged isolation [GeV];Events",25,0,0.10)); - ht.addHist("vystar", new TH1F("vectorbosonystar", ";y-(1/2)(y_{j1}+y_{j2});Events",25,0,5)); - ht.addHist("njets", new TH1F("njets", ";Jet multiplicity;Events",10,-0.5,9.5)); - ht.addHist("mjj", new TH1F("mjj", ";Dijet invariant mass [GeV];Events",40,0,4000)); - ht.addHist("leadpt", new TH1F("leadpt", ";Leading jet p_{T} [GeV];Events",25,0,500)); - ht.addHist("subleadpt", new TH1F("subleadpt" , ";Sub-leading jet p_{T} [GeV];Events",25,0,500)); - ht.addHist("drj1b", new TH1F("drj1b", ";#DeltaR(j_{1},boson);Events",25,0,8)); - ht.addHist("drj2b", new TH1F("drj2b" , ";#DeltaR(j_{2},boson);Events",25,0,8)); - ht.addHist("leadpumva", new TH1F("leadpumva", ";Pileup MVA;Events",25,-1,1)); - ht.addHist("subleadpumva", new TH1F("subleadpumva" , ";Pileup MVA;Events",25,-1,1)); - ht.addHist("centraleta", new TH1F("centraleta", ";Most central jet |#eta|;Events",25,0,5)); - ht.addHist("forwardeta", new TH1F("forwardeta", ";Most forward jet |#eta|;Events",25,0,5)); - ht.addHist("dijetpt", new TH1F("dijetpt", ";Dijet p_{T} [GeV];Events",20,0,1000)); - ht.addHist("detajj", new TH1F("detajj" , ";#Delta#eta(J,J);Events",20,0,8)); - ht.addHist("dphijj", new TH1F("dphijj" , ";#Delta#phi(J,J) [rad];Events",20,-3.15,3.15)); - ht.addHist("ht", new TH1F("ht", ";H_{T} [GeV];Events",20,0,4000)); - ht.addHist("mht", new TH1F("mht", ";Missing H_{T} [GeV];Events",20,0,500)); - ht.addHist("balance", new TH1F("balance", ";System p_{T} balance [GeV];Events",20,0,300)); - ht.addHist("sphericity", new TH1F("sphericity", ";Sphericity;Events",20,0,1.0)); - ht.addHist("aplanarity", new TH1F("aplanarity", ";Aplanarity;Events",20,0,1.0)); - ht.addHist("C", new TH1F("C", ";C;Events",20,0,1.0)); - ht.addHist("D", new TH1F("D", ";D;Events",20,0,1.0)); - ht.addHist("isotropy", new TH1F("isotropy", ";Isotropy;Events",20,0,1.0)); - ht.addHist("circularity", new TH1F("circularity", "Circularity;;Events",20,0,1.0)); - - std::cout << "init done" << std::endl; - /////////////////////// // LOOP OVER EVENTS // ///////////////////// - - //EVENT SELECTION WRAPPER - SelectionTool selector(filename, debug, triggerList,SelectionTool::VBF); - + float xsec = getXsec(); for (Int_t iev=0;iev<nentries;iev++) { t->GetEntry(iev); if(iev%10000==0) printf ("\r [%3.0f%%] done", 100.*(float)iev/(float)nentries); - std::vector<double>plotwgts(1,1.0); - ht.fill("qscale",ev.g_qscale,plotwgts); + ht->fill("qscale",ev.g_qscale,plotwgts); //assign randomly a run period - TString period = lumi.assignRunPeriod(); + TString period = lumi->assignRunPeriod(); ////////////////// // CORRECTIONS // //////////////// - jec.smearJetEnergies(ev); + jec->smearJetEnergies(ev); /////////////////////////// // RECO LEVEL SELECTION // ///////////////////////// - TString chTag = selector.flagFinalState(ev); - std::vector<Particle> &photons = selector.getSelPhotons(); - std::vector<Particle> &leptons = selector.getSelLeptons(); - std::vector<Jet> &alljets = selector.getJets(); + TString chTag = selector->flagFinalState(ev); + std::vector<Particle> &photons = selector->getSelPhotons(); + std::vector<Particle> &leptons = selector->getSelLeptons(); + std::vector<Jet> &alljets = selector->getJets(); std::vector<Jet> jets; + + //Pileup jet id for(auto j : alljets) { int idx=j.getJetIndex(); int jid=ev.j_id[idx]; @@ -171,20 +65,26 @@ void RunVBFVectorBoson(TString filename, if(!passLoosePu) continue; jets.push_back(j); } + //Category selection if(chTag!="A" && chTag!="MM") continue; + float cat[6] = {0,0,0,0,0,0}; + if(chTag == "A") cat[1] = 1; + if(chTag == "MM") cat[0] = 1; //jet related variables and selection - float mjj(jets.size()>=2 ? (jets[0]+jets[1]).M() : 0.); - float detajj(jets.size()>=2 ? fabs(jets[0].Eta()-jets[1].Eta()) : -1.); - float dphijj(jets.size()>=2 ? jets[0].DeltaPhi(jets[1]) : -1.); - float jjpt=(jets.size()>=2 ? (jets[0]+jets[1]).Pt() : 0.); - float scalarht(0.); + mjj = (jets.size()>=2 ? (jets[0]+jets[1]).M() : 0.); + detajj = (jets.size()>=2 ? fabs(jets[0].Eta()-jets[1].Eta()) : -1.); + dphijj = (jets.size()>=2 ? jets[0].DeltaPhi(jets[1]) : -1.); + jjpt = (jets.size()>=2 ? (jets[0]+jets[1]).Pt() : 0.); + + scalarht = 0.; TLorentzVector mhtP4(0,0,0,0); + mht = 0; for(auto j : jets) { scalarht += j.Pt(); mhtP4 += j; } - float mht(mhtP4.Pt()); + mht = mhtP4.Pt(); bool passJets(jets.size()>=2 && mjj>400); bool passVBFJetsTrigger(passJets && detajj>3.0); @@ -192,19 +92,18 @@ void RunVBFVectorBoson(TString filename, //for the photon refine also the category according to the trigger bit TLorentzVector boson(0,0,0,0); bool isHighPt(false),isVBF(false),isHighPtAndVBF(false),isBosonPlusOneJet(false); - float sihih(0),chiso(0),r9(0),hoe(0); - if(chTag=="A") { - + sihih = 0, chiso = 0 ,r9 = 0, hoe = 0; + if(chTag=="A") { boson += photons[0]; sihih = ev.gamma_sieie[photons[0].originalReference()]; chiso = ev.gamma_chargedHadronIso[photons[0].originalReference()]; r9 = ev.gamma_r9[photons[0].originalReference()]; hoe = ev.gamma_hoe[photons[0].originalReference()]; - isVBF = (selector.hasTriggerBit("HLT_Photon75_R9Id90_HE10_IsoM_EBOnly_PFJetsMJJ300DEta3_v", ev.triggerBits) + isVBF = (selector->hasTriggerBit("HLT_Photon75_R9Id90_HE10_IsoM_EBOnly_PFJetsMJJ300DEta3_v", ev.triggerBits) && photons[0].Pt()>75 && fabs(photons[0].Eta())<1.442 && passVBFJetsTrigger); - isHighPt = ( selector.hasTriggerBit("HLT_Photon200_v", ev.triggerBits) + isHighPt = ( selector->hasTriggerBit("HLT_Photon200_v", ev.triggerBits) && photons[0].Pt()>200 ); isHighPtAndVBF = (isHighPt && isVBF); isBosonPlusOneJet=(isHighPt && alljets.size()==1); @@ -216,8 +115,7 @@ void RunVBFVectorBoson(TString filename, isHighPtAndVBF = false; } - } - else { + } else { boson += leptons[0]; boson += leptons[1]; isVBF = boson.Pt()>75 && fabs(boson.Rapidity())<1.442 && passVBFJetsTrigger; @@ -225,20 +123,22 @@ void RunVBFVectorBoson(TString filename, isHighPtAndVBF = (isHighPt && isVBF); isBosonPlusOneJet=(isHighPt && alljets.size()==1); } + if(!isVBF && !isHighPt && !isBosonPlusOneJet) continue; std::vector<TString> chTags; - if(isVBF) chTags.push_back("VBF"+chTag); - if(isHighPt) chTags.push_back("HighPt"+chTag); - if(isHighPtAndVBF) chTags.push_back("HighPtVBF"+chTag); - if(isBosonPlusOneJet) chTags.push_back("V1J"+chTag); + if(isVBF) {chTags.push_back("VBF"+chTag); cat[2]=1;} + if(isHighPt) {chTags.push_back("HighPt"+chTag); cat[3]=1;} + if(isHighPtAndVBF) {chTags.push_back("HighPtVBF"+chTag); cat[4]=1;} + if(isBosonPlusOneJet) {chTags.push_back("V1J"+chTag); cat[5]=1;} + category.set(cat); //leptons and boson - double mindrl(9999.); + mindrl = 9999.; for(auto &l: leptons) mindrl=min(l.DeltaR(boson),mindrl); //system variables and event shapes - float ystar(0),balance(0); + ystar = 0,balance = 0; if(passJets) { ystar=boson.Rapidity()-0.5*(jets[0].Rapidity()+jets[1].Rapidity()); balance=(boson+jets[0]+jets[1]).Pt(); @@ -248,6 +148,13 @@ void RunVBFVectorBoson(TString filename, for(size_t ij=0; ij<min(size_t(2),jets.size());ij++) { inputVectors.push_back( math::XYZVector(jets[ij].Px(),jets[ij].Py(),jets[ij].Pz()) ); EventShapeVariables esv(inputVectors); + isotropy = esv.isotropy(); + circularity = esv.circularity(); + sphericity = esv.sphericity(1.); + aplanarity = esv.aplanarity(1.); + C = esv.C(1.); + D = esv.D(1.); + //////////////////// // EVENT WEIGHTS // @@ -259,23 +166,23 @@ void RunVBFVectorBoson(TString filename, wgt = (normH? normH->GetBinContent(1) : 1.0); // pu weight - ht.fill("puwgtctr",0,plotwgts); - double puWgt(lumi.pileupWeight(ev.g_pu,period)[0]); + ht->fill("puwgtctr",0,plotwgts); + double puWgt(lumi->pileupWeight(ev.g_pu,period)[0]); std::vector<double>puPlotWgts(1,puWgt); - ht.fill("puwgtctr",1,puPlotWgts); + ht->fill("puwgtctr",1,puPlotWgts); // photon trigger*selection weights float trigSF(1.0), selSF(1.0); if(chTag=="A") { - trigSF *= gammaEffWR.getTriggerCorrection({},photons,{}, period).first; - selSF *= gammaEffWR.getOfflineCorrection(photons[0], period).first; + trigSF *= gammaEffWR->getTriggerCorrection({},photons,{}, period).first; + selSF *= gammaEffWR->getOfflineCorrection(photons[0], period).first; } else { - trigSF *=gammaEffWR.getTriggerCorrection(leptons,{},{}, period).first; - selSF *=gammaEffWR.getOfflineCorrection(leptons[0], period).first; - selSF *=gammaEffWR.getOfflineCorrection(leptons[1], period).first; + trigSF *=gammaEffWR->getTriggerCorrection(leptons,{},{}, period).first; + selSF *=gammaEffWR->getOfflineCorrection(leptons[0], period).first; + selSF *=gammaEffWR->getOfflineCorrection(leptons[1], period).first; } wgt *= puWgt*trigSF*selSF; @@ -289,7 +196,6 @@ void RunVBFVectorBoson(TString filename, //control histograms for( auto c : chTags) { - std::vector<double> cplotwgts(plotwgts); //photon pT weighting @@ -305,48 +211,11 @@ void RunVBFVectorBoson(TString filename, cplotwgts[0]*=photonPtWgt; } - ht.fill("nvtx", ev.nvtx, cplotwgts,c); - - //boson histos - ht.fill("vpt", boson.Pt(), cplotwgts,c); - ht.fill("vy", boson.Rapidity(), cplotwgts,c); - ht.fill("sihih", sihih, cplotwgts,c); - ht.fill("r9", r9, cplotwgts,c); - ht.fill("hoe", hoe, cplotwgts,c); - ht.fill("chiso", chiso, cplotwgts,c); - ht.fill("mindrl", mindrl, cplotwgts,c); - - //jet histos - double minEta(9999),maxEta(-9999); - for(size_t ij=0; ij<min(size_t(2),jets.size());ij++) { - TString jtype(ij==0?"lead":"sublead"); - ht.fill(jtype+"pt", jets[ij].Pt(), cplotwgts,c); - ht.fill(jtype+"pumva", jets[ij].getPUMVA(), cplotwgts,c); - ht.fill("dr"+jtype+"b", jets[ij].DeltaR(boson), cplotwgts,c); - minEta=min(minEta,fabs(jets[ij].Eta())); - maxEta=max(maxEta,fabs(jets[ij].Eta())); - } - ht.fill("njets", jets.size(), cplotwgts,c); - ht.fill("ht", scalarht, cplotwgts,c); - ht.fill("mht", mht, cplotwgts,c); - ht.fill("centraleta", minEta, cplotwgts,c); - ht.fill("forwardeta", maxEta, cplotwgts,c); - ht.fill("dijetpt", jjpt, cplotwgts,c); - ht.fill("detajj", detajj, cplotwgts,c); - ht.fill("dphijj", dphijj, cplotwgts,c); - ht.fill("mjj", mjj, cplotwgts,c); - - //visible system histos - ht.fill("vystar", ystar, cplotwgts,c); - ht.fill("balance", balance, cplotwgts,c); - ht.fill("isotropy", esv.isotropy(), cplotwgts,c); - ht.fill("circularity", esv.circularity(), cplotwgts,c); - ht.fill("sphericity", esv.sphericity(1.), cplotwgts,c); - ht.fill("aplanarity", esv.aplanarity(1.), cplotwgts,c); - ht.fill("C", esv.C(1.), cplotwgts,c); - ht.fill("D", esv.D(1.), cplotwgts,c); - - } + //What is the final weight? 0 or 1 in the array? + evtWeight = cplotwgts[0]*xsec; + training = useForTraining(); + fill( ev, boson, jets, cplotwgts, c); + } } } @@ -361,23 +230,250 @@ void RunVBFVectorBoson(TString filename, } //save histos to file - fOut->cd(); - for (auto& it : ht.getPlots()) { - for(auto &wit : photonPtWgtCtr){ - if(!it.first.Contains(wit.first)) continue; - cout << "Scaling " << it.first << " by "<< wit.second.first <<endl; - it.second->Scale(wit.second.first); - break; - } - it.second->SetDirectory(fOut); it.second->Write(); - } - for (auto& it : ht.get2dPlots()) { - for(auto &wit : photonPtWgtCtr){ - if(!it.first.Contains(wit.first)) continue; - it.second->Scale(wit.second.first); - break; - } - it.second->SetDirectory(fOut); it.second->Write(); + saveHistos(); + if(skimtree){ + fMVATree->cd(); + newTree->Write(); + fMVATree->Close(); } - fOut->Close(); +} +void VBFVectorBoson::saveHistos(){ + fOut->cd(); + for (auto& it : ht->getPlots()) { + for(auto &wit : photonPtWgtCtr){ + if(!it.first.Contains(wit.first)) continue; + cout << "Scaling " << it.first << " by "<< wit.second.first <<endl; + it.second->Scale(wit.second.first); + break; + } + it.second->SetDirectory(fOut); it.second->Write(); + } + for (auto& it : ht->get2dPlots()) { + for(auto &wit : photonPtWgtCtr){ + if(!it.first.Contains(wit.first)) continue; + it.second->Scale(wit.second.first); + break; + } + it.second->SetDirectory(fOut); it.second->Write(); + } + fOut->Close(); +} + +void VBFVectorBoson::readTree(){ + f = TFile::Open(filename); + triggerList=(TH1 *)f->Get("analysis/triggerList"); + t = (TTree*)f->Get("analysis/data"); + attachToMiniEventTree(t,ev,true); + nentries = t->GetEntriesFast(); + if (debug) nentries = 10000; //restrict number of entries for testing + t->GetEntry(0); + isQCDEMEnriched = filename.Contains("MC13TeV_QCDEM"); +} + +void VBFVectorBoson::prepareOutput(){ + baseName=gSystem->BaseName(outname); + TString dirName=gSystem->DirName(outname); + + if(skimtree){ + fMVATree=TFile::Open(dirName+"/MVATree_"+baseName,"RECREATE"); + newTree = t->CloneTree(0); + } + fOut=TFile::Open(dirName+"/"+baseName,"RECREATE"); + fOut->cd(); +} + +void VBFVectorBoson::bookHistograms(){ + ht = new HistTool(0); + ht->addHist("puwgtctr", new TH1F("puwgtctr", ";Weight sums;Events",2,0,2)); + ht->addHist("qscale", new TH1F("qscale", ";Q^{2} scale;Events",100,0,2000)); + ht->addHist("nvtx", new TH1F("nvtx", ";Vertex multiplicity;Events",100,-0.5,99.5)); + ht->addHist("vpt", new TH1F("vectorbosonPt", ";Boson p_{T}[GeV];Events",25,0,500)); + ht->addHist("vy", new TH1F("vectorbosony", ";Boson rapidity;Events",25,0,3)); + ht->addHist("mindrl", new TH1F("mindrl", ";min #Delta R(boson,lepton);Events",25,0,6)); + ht->addHist("sihih", new TH1F("sihih", ";#sigma(i#eta,i#eta);Events",50,0,0.1)); + ht->addHist("hoe", new TH1F("hoe", ";h/e;Events",25,0,0.1)); + ht->addHist("r9", new TH1F("r9", ";r9;Events",25,0,1.0)); + ht->addHist("chiso", new TH1F("chiso", ";Charged isolation [GeV];Events",25,0,0.10)); + ht->addHist("vystar", new TH1F("vectorbosonystar", ";y-(1/2)(y_{j1}+y_{j2});Events",25,0,5)); + ht->addHist("njets", new TH1F("njets", ";Jet multiplicity;Events",10,-0.5,9.5)); + ht->addHist("mjj", new TH1F("mjj", ";Dijet invariant mass [GeV];Events",40,0,4000)); + ht->addHist("leadpt", new TH1F("leadpt", ";Leading jet p_{T} [GeV];Events",25,0,500)); + ht->addHist("subleadpt", new TH1F("subleadpt" , ";Sub-leading jet p_{T} [GeV];Events",25,0,500)); + ht->addHist("drj1b", new TH1F("drj1b", ";#DeltaR(j_{1},boson);Events",25,0,8)); + ht->addHist("drj2b", new TH1F("drj2b" , ";#DeltaR(j_{2},boson);Events",25,0,8)); + ht->addHist("leadpumva", new TH1F("leadpumva", ";Pileup MVA;Events",25,-1,1)); + ht->addHist("subleadpumva", new TH1F("subleadpumva" , ";Pileup MVA;Events",25,-1,1)); + ht->addHist("centraleta", new TH1F("centraleta", ";Most central jet |#eta|;Events",25,0,5)); + ht->addHist("forwardeta", new TH1F("forwardeta", ";Most forward jet |#eta|;Events",25,0,5)); + ht->addHist("dijetpt", new TH1F("dijetpt", ";Dijet p_{T} [GeV];Events",20,0,1000)); + ht->addHist("detajj", new TH1F("detajj" , ";#Delta#eta(J,J);Events",20,0,8)); + ht->addHist("dphijj", new TH1F("dphijj" , ";#Delta#phi(J,J) [rad];Events",20,-3.15,3.15)); + ht->addHist("ht", new TH1F("ht", ";H_{T} [GeV];Events",20,0,4000)); + ht->addHist("mht", new TH1F("mht", ";Missing H_{T} [GeV];Events",20,0,500)); + ht->addHist("balance", new TH1F("balance", ";System p_{T} balance [GeV];Events",20,0,300)); + ht->addHist("sphericity", new TH1F("sphericity", ";Sphericity;Events",20,0,1.0)); + ht->addHist("aplanarity", new TH1F("aplanarity", ";Aplanarity;Events",20,0,1.0)); + ht->addHist("C", new TH1F("C", ";C;Events",20,0,1.0)); + ht->addHist("D", new TH1F("D", ";D;Events",20,0,1.0)); + ht->addHist("isotropy", new TH1F("isotropy", ";Isotropy;Events",20,0,1.0)); + ht->addHist("circularity", new TH1F("circularity", "Circularity;;Events",20,0,1.0)); + // Study of jet variables + ht->addHist("jet_c1_00", new TH1F("jet_c1_00", ";Jet shape var. c1_00;Jets",100,-1,1)); + ht->addHist("jet_c1_02", new TH1F("jet_c1_02", ";Jet shape var. c1_02;Jets",100,-1,1)); + ht->addHist("jet_c1_05", new TH1F("jet_c1_05", ";Jet shape var. c1_05;Jets",100,-1,1)); + ht->addHist("jet_c2_00", new TH1F("jet_c2_00", ";Jet shape var. c2_00;Jets",100,-1,1)); + ht->addHist("jet_c2_02", new TH1F("jet_c2_02", ";Jet shape var. c2_02;Jets",100,-1,1)); + ht->addHist("jet_c2_05", new TH1F("jet_c2_05", ";Jet shape var. c2_05;Jets",100,-1,1)); + ht->addHist("jet_c3_00", new TH1F("jet_c3_00", ";Jet shape var. c3_00;Jets",100,-1,1)); + ht->addHist("jet_c3_02", new TH1F("jet_c3_02", ";Jet shape var. c3_02;Jets",100,-1,1)); + ht->addHist("jet_c3_05", new TH1F("jet_c3_05", ";Jet shape var. c3_05;Jets",100,-1,1)); + ht->addHist("jet_zg", new TH1F("jet_zg", ";Jet shape var. zg;Jets",100,-1,1)); + ht->addHist("jet_gaptd", new TH1F("jet_gaptd", ";Jet shape var. gaptd;Jets",100,-1,1)); + ht->addHist("jet_gawidth", new TH1F("jet_gawidth", ";Jet shape var. gawidth;Jets",100,-1,1)); + //additional variables from https://link.springer.com/content/pdf/10.1140/epjc/s10052-017-5315-6.pdf + ht->addHist("jjetas", new TH1F("jjetas", ";#eta_{j1}#eta_{j2};Events",200,-25,25)); + ht->addHist("centjy", new TH1F("centjy", ";Central jet rapidity;Jets",25,0,3)); + ht->addHist("ncentj", new TH1F("ncentjj", ";Number of central jets;Events",10,-0.5,9.5)); + ht->addHist("dphivj0", new TH1F("dphivj0", ";#Delta#phi(V,j0);Jets",20,0,4)); + ht->addHist("dphivj1", new TH1F("dphivj1", ";#Delta#phi(V,j1);Jets",20,0,4)); + ht->addHist("dphivj2", new TH1F("dphivj2", ";#Delta#phi(V,j2);Jets",20,0,4)); + ht->addHist("dphivj3", new TH1F("dphivj3", ";#Delta#phi(V,j3);Jets",20,0,4)); +} +void VBFVectorBoson::setGammaZPtWeights(){ + TString wgtUrl("${CMSSW_BASE}/src/TopLJets2015/TopAnalysis/test/analysis/VBFVectorBoson/raw/plots/ratio_plotter.root"); + gSystem->ExpandPathName(wgtUrl); + TFile *wgtF=TFile::Open(wgtUrl); + if(wgtF) { + cout << "Reading photon/Z pT weights" << endl; + TString pfix(baseName.Contains("Data13TeV_") ? "" : "_mc_MC"); + photonPtWgts["VBFA"] = new TGraph((TH1* )wgtF->Get("VBFA_vectorbosonPt_ratio/VBFA_vectorbosonPt"+pfix)); + photonPtWgts["HighPtA"] = new TGraph((TH1* )wgtF->Get("HighPtA_vectorbosonPt_ratio/HighPtA_vectorbosonPt"+pfix)); + photonPtWgtCtr["VBFA"] = std::pair<double,double>(0.0,0.0); + photonPtWgtCtr["HighPtA"] = std::pair<double,double>(0.0,0.0); + wgtF->Close(); + } else { + cout << "Requested to reweight photon spectrum but could not find " << wgtUrl << endl + << "Proceeding without" << endl; + } +} +void VBFVectorBoson::loadCorrections(){ + lumi = new LumiTools(era,genPU); + gammaEffWR = new EfficiencyScaleFactorsWrapper(filename.Contains("Data13TeV"),era); + jec = new JECTools(era); + if(anFlag>0) this->setGammaZPtWeights(); +} +void VBFVectorBoson::addMVAvars(){ + newTree->Branch("centralEta", ¢raleta); + newTree->Branch("mjj", &mjj); + newTree->Branch("detajj", &detajj); + newTree->Branch("jjpt", &jjpt); + newTree->Branch("dphijj", &dphijj); + newTree->Branch("forwardeta", &forwardeta); + newTree->Branch("jjetas", &jjetas); + newTree->Branch("centjy", ¢jy); + newTree->Branch("ncentjj", &ncentjj); + newTree->Branch("dphivj0", &dphivj0); + newTree->Branch("dphivj1", &dphivj1); + newTree->Branch("dphivj2", &dphivj2); + newTree->Branch("dphivj3", &dphivj3); + newTree->Branch("evtWeight", &evtWeight); + newTree->Branch("mht", &mht); + newTree->Branch("balance", &balance); + newTree->Branch("ht", &scalarht); + newTree->Branch("isotropy", &isotropy); + newTree->Branch("circularity",&circularity); + newTree->Branch("sphericity",&sphericity); + newTree->Branch("aplanarity",&aplanarity); + newTree->Branch("C",&C); + newTree->Branch("D",&D); + newTree->Branch("training",&training); + newTree->Branch("category", &category, "MM:A:VBF:HighPt:HighPtVBF:V1J"); +} + +void VBFVectorBoson::fill(MiniEvent_t ev, TLorentzVector boson, std::vector<Jet> jets, std::vector<double> cplotwgts, TString c){ + ht->fill("nvtx", ev.nvtx, cplotwgts,c); + + //boson histos + ht->fill("vpt", boson.Pt(), cplotwgts,c); + ht->fill("vy", boson.Rapidity(), cplotwgts,c); + ht->fill("sihih", sihih, cplotwgts,c); + ht->fill("r9", r9, cplotwgts,c); + ht->fill("hoe", hoe, cplotwgts,c); + ht->fill("chiso", chiso, cplotwgts,c); + ht->fill("mindrl", mindrl, cplotwgts,c); + + //jet histos + centraleta = 9999; + forwardeta = -9999; + for(size_t ij=0; ij<min(size_t(2),jets.size());ij++) { + TString jtype(ij==0?"lead":"sublead"); + ht->fill(jtype+"pt", jets[ij].Pt(), cplotwgts,c); + ht->fill(jtype+"pumva", jets[ij].getPUMVA(), cplotwgts,c); + ht->fill("dr"+jtype+"b", jets[ij].DeltaR(boson), cplotwgts,c); + ht->fill("jet_c1_00", ev.j_c1_00[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c1_02", ev.j_c1_02[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c1_05", ev.j_c1_05[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c2_00", ev.j_c2_00[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c2_02", ev.j_c2_02[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c2_05", ev.j_c2_05[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c3_00", ev.j_c3_00[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c3_02", ev.j_c3_02[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_c3_05", ev.j_c3_05[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_zg", ev.j_zg[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_gaptd", ev.j_gaptd[jets[ij].getJetIndex()] , cplotwgts,c); + ht->fill("jet_gawidth", ev.j_gawidth[jets[ij].getJetIndex()] , cplotwgts,c); + centraleta=min(centraleta,float(fabs(jets[ij].Eta()))); + forwardeta=max(forwardeta,float(fabs(jets[ij].Eta()))); + } + jjetas = 9999; + dphivj0 = 9999; dphivj1 = 9999; + if(jets.size() >= 2){ + jjetas = jets[0].Eta()*jets[1].Eta(); + dphivj0 = fabs(jets[0].DeltaPhi(boson)); + dphivj1 = fabs(jets[1].DeltaPhi(boson)); + ht->fill("jjetas", jjetas, cplotwgts,c); + ht->fill("dphivj0", dphivj0 , cplotwgts,c); + ht->fill("dphivj1", dphivj1 , cplotwgts,c); + } + dphivj2 = 9999; dphivj3 = 9999; + centjy = 9999; ncentjj = 0; + if(jets.size() > 2){ + dphivj2 = fabs(jets[2].DeltaPhi(boson)); + ht->fill("dphivj2", dphivj2 , cplotwgts,c); + for(unsigned int iJet = 2; iJet < jets.size(); iJet++){ + float dy = fabs(jets[0].Rapidity() - jets[1].Rapidity())/2; + float sumy = (jets[0].Rapidity() + jets[1].Rapidity())/2; + if(fabs(jets[iJet].Rapidity() - sumy) < dy){ + centjy = jets[iJet].Rapidity(); + ht->fill("centjy",centjy, cplotwgts,c); + ncentjj++; + } + } + ht->fill("ncentj", ncentjj, cplotwgts, c); + } + if(jets.size() > 3){ + dphivj3 = fabs(jets[3].DeltaPhi(boson)); + ht->fill("dphivj3", dphivj3 , cplotwgts,c); + } + ht->fill("njets", jets.size(), cplotwgts,c); + ht->fill("ht", scalarht, cplotwgts,c); + ht->fill("mht", mht, cplotwgts,c); + ht->fill("centraleta", centraleta, cplotwgts,c); + ht->fill("forwardeta", forwardeta, cplotwgts,c); + ht->fill("dijetpt", jjpt, cplotwgts,c); + ht->fill("detajj", detajj, cplotwgts,c); + ht->fill("dphijj", dphijj, cplotwgts,c); + ht->fill("mjj", mjj, cplotwgts,c); + + //visible system histos + ht->fill("vystar", ystar, cplotwgts,c); + ht->fill("balance", balance, cplotwgts,c); + ht->fill("isotropy", isotropy, cplotwgts,c); + ht->fill("circularity", circularity, cplotwgts,c); + ht->fill("sphericity", sphericity, cplotwgts,c); + ht->fill("aplanarity", aplanarity, cplotwgts,c); + ht->fill("C", C, cplotwgts,c); + ht->fill("D", D, cplotwgts,c); + cout<< "hist filling is done ... "<<endl; + if(skimtree) newTree->Fill(); } diff --git a/TopAnalysis/test/analysis/steerVBFVectorBoson.sh b/TopAnalysis/test/analysis/steerVBFVectorBoson.sh index b9188671a3536ddd26d788b23681ea61f2cb8925..5a1425667360ceb6b94887fee13a9960f5d704fd 100644 --- a/TopAnalysis/test/analysis/steerVBFVectorBoson.sh +++ b/TopAnalysis/test/analysis/steerVBFVectorBoson.sh @@ -27,24 +27,25 @@ NC='\e[0m' case $WHAT in TESTSEL ) - #input=${eosdir}/MC13TeV_DY50toInf/MergedMiniEvents_0_ext0.root - #output=MC13TeV_DY.root - #tag="--tag MC13TeV_DY50toInf" - input=${eosdir}/MC13TeV_GJets_HT100to200/MergedMiniEvents_0_ext0.root - #output=MC13TeV_GJets_HT100to200.root - #tag="--tag MC13TeV_GJets_HT100to200 --flag 1" - output=MC13TeV_GJets_HT100to200_raw.root - tag="--tag MC13TeV_GJets_HT100to200" + input=${eosdir}/MC13TeV_DY50toInf/MergedMiniEvents_4_ext0.root + output=MC13TeV_DY4Jets50toInf.root + tag="--tag MC13TeV_DY50toInf" + #input=${eosdir}/Data13TeV_SingleMuon_2017D/MergedMiniEvents_0_ext0.root + #output=Data13TeV_SingleMuon_2017D.root + #input=${eosdir}/MC13TeV_GJets_HT100to200_DR04/MergedMiniEvents_0_ext0.root + #output=MC13TeV_GJets_HT100to200_DR04.root \ #input=${eosdir}/Data13TeV_SingleMuon_2017C/MergedMiniEvents_0_ext0.root #output=Data13TeV_SingleMuon_2017C.root python scripts/runLocalAnalysis.py \ -i ${input} -o ${output} ${tag} \ - --njobs 1 -q local --debug \ + --njobs 1 -q local --debug --mvatree \ --era era2017 -m VBFVectorBoson::RunVBFVectorBoson --ch 0 --runSysts; ;; SEL ) python scripts/runLocalAnalysis.py -i ${eosdir} \ +# --only DY,EWKZJJ,GJets_HT,QCDEM \ +# -o ${outdir} \ -o ${outdir}/raw \ -q ${queue} \ --era era2017 -m VBFVectorBoson::RunVBFVectorBoson --ch 0 --runSysts; @@ -62,10 +63,12 @@ case $WHAT in ;; PLOT ) +# commonOpts="-i ${outdir} --puNormSF puwgtctr -l 1 --saveLog" +# python scripts/plotter.py ${commonOpts} -j data/era2017/vbf_samples.json --noStack --skip TT,ZZ,WW,WZ,Single,QCD,GJets,DY1Jets,DY2Jets,DY3Jets,DY4Jets,DY50toInf_HT -O ${outdir}/plots_DY ; commonOpts="-i ${outdir}/${EXTRA} --puNormSF puwgtctr -l ${fulllumi} --saveLog --mcUnc ${lumiUnc} --lumiSpecs VBFA:${vbflumi}" python scripts/plotter.py ${commonOpts} -j data/era2017/vbf_samples.json; - #python test/analysis/computeVBFRatios.py -i ${outdir}/${EXTRA}/plots/plotter.root -o ${outdir}/${EXTRA}/plots/ratio_plotter.root python test/analysis/computeVBFRatios.py -t -i ${outdir}/${EXTRA}/plots/plotter.root -o ${outdir}/${EXTRA}/plots/trigger_ratio_plotter.root + ;; WWW )