diff --git a/README b/README index 83ad82b..ca972ac 100644 --- a/README +++ b/README @@ -1,2 +1,13 @@ -#how to use the macro -root -l VHPlotter.C+\(\"histoname\",1\); +# how to run analysis code to read ntuples, apply selection and svae histogram + +choose flag_xxx = 1/0 to decide which sample you want to run +choose macro = xAna_ele/xAna_mu to decide electron/muon channel you want to run +choose outputFolder = output_ele/output_mu for the output file folder + +Then just one Go +./globalRun.sh + + +# how to use the make plot, just ONE GO +./runPlot.sh + diff --git a/VHPlotter.C b/VHPlotter.C index 18fa320..48c7257 100644 --- a/VHPlotter.C +++ b/VHPlotter.C @@ -35,19 +35,24 @@ #include "TLegend.h" #include "TLatex.h" #include "TVirtualFitter.h" +#include "TString.h" #include #include #include #include #include -#include "interface/setNCUStyle.h" -#include "interface/readHists.h" -string path = "/afs/cern.ch/user/v/vieri/work/ZH_Analysis/CMSSW_8_0_15/src/ZHAnalysis/DataMC_comparison/"; +//string path = "/afs/cern.ch/work/y/yuchang/public/show_vieri_DataMC_comparison/output_with_PUweight/"; + +string path = "/afs/cern.ch/work/y/yuchang/ZH_analysis_with_2016_data/CMSSW_8_0_8/src/DataMC_comparison/test_ele_with_PUweight/output_ele/"; + +//string path = "/afs/cern.ch/user/v/vieri/work/ZH_Analysis/CMSSW_8_0_15/src/ZHAnalysis/DataMC_comparison/"; + +//void VHPlotter(string title="", int plot=0) { +void VHPlotter(string title="", int plot=0, TCanvas* c1=0 ) { -void VHPlotter(string title="", int plot=0) { string subdir="0"; string postfix=""; @@ -173,8 +178,8 @@ hs->Add(h_mc3); hs->Add(h_mc2); hs->Add(h_mcDY); -TCanvas* c1 = 0; -c1 = new TCanvas("c","c",10,10,800,600); +//TCanvas* c1 = 0; +//c1 = new TCanvas("c","c",10,10,800,600); c1->cd(); TPad *pad1 = new TPad("pad1","pad1",0.0,0.3,1.0,1.0); @@ -187,10 +192,10 @@ hs->Draw("HIST"); hs->GetYaxis()->SetTitle("Events"); hs->GetYaxis()->SetTitleSize(0.05); hs->GetYaxis()->SetLabelSize(0.045); -hs->GetYaxis()->SetTitleOffset(0.7); +hs->GetYaxis()->SetTitleOffset(1.0);// 0.7 hs->SetMinimum(8); hs->SetMaximum(1.2*hs->GetMaximum()); - +if (title=="ZHmass") { hs->GetXaxis()->SetRangeUser(0, 3000); } h_data->Draw("EPX0SAMES"); h_data->SetMarkerColor(kBlack); @@ -199,7 +204,10 @@ h_data->SetMarkerSize (1.0); h_data->SetStats(0); TLegend *leg; -leg = new TLegend(0.65, 0.547, 0.91, 0.88); + +if (title=="FATjetTau2dvTau1") { leg = new TLegend(0.15, 0.547, 0.41, 0.88); } +else { leg = new TLegend(0.65, 0.547, 0.91, 0.88);} +//leg = new TLegend(0.65, 0.547, 0.91, 0.88); leg->SetBorderSize(0); leg->SetEntrySeparation(0.01); leg->SetFillColor(0); @@ -216,6 +224,13 @@ leg->AddEntry(h_mc8,"ZH powheg","f"); leg->Draw(); +TLatex *lar = new TLatex(); +lar->SetNDC(kTRUE); +lar->SetTextSize(0.04); +lar->SetLineWidth(5); +lar->DrawLatex(0.14, 0.94, "CMS #it{#bf{2016}}"); +lar->DrawLatex(0.60, 0.94, "L = 4.327 fb^{-1} at #sqrt{s} = 13 TeV"); + pad1->Update(); c1->Update(); @@ -223,7 +238,7 @@ c1->cd(); TH1F *h_ratio = (TH1F*)h_data->Clone("h_ratio"); -TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.3); +TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.3);// 0.3 pad2->SetTopMargin(0); pad2->SetBottomMargin(0.3); pad2->Draw(); @@ -265,7 +280,7 @@ if (title=="Zpt") { h_ratio->GetXaxis ()->SetTitle("FAT Jet #tau_{1}"); } else if (title=="FATjetTau2") { h_ratio->GetXaxis ()->SetTitle("FAT Jet #tau_{2}"); -} else if (title=="FATjetTau12dvTau1") { +} else if (title=="FATjetTau2dvTau1") { h_ratio->GetXaxis ()->SetTitle("FAT Jet #tau_{21}"); } else if (title=="FATsubjetPt") { h_ratio->GetXaxis ()->SetTitle("FAT SubJet p_{T} [GeV/c]"); @@ -275,9 +290,9 @@ if (title=="Zpt") { h_ratio->GetXaxis ()->SetTitle("FAT SubJet Soft Drop CSV [GeV/c]"); } else if (title=="ZHmass") { h_ratio->GetXaxis ()->SetTitle("ZH invariant mass [GeV/c^{2}]"); +} - - +if (title=="ZHmass") { h_ratio->GetXaxis()->SetRangeUser(0, 3000); } h_ratio->GetXaxis()->SetTitleSize(0.11); h_ratio->GetXaxis()->SetLabelFont(42); h_ratio->GetXaxis()->SetLabelSize(0.10); @@ -298,12 +313,14 @@ OLine->SetLineColor(kRed); OLine->SetLineWidth(2); OLine->Draw(); +/* TLatex *lar = new TLatex(); lar->SetNDC(kTRUE); lar->SetTextSize(0.04); lar->SetLineWidth(5); lar->DrawLatex(0.14, 0.94, "CMS #it{#bf{2016}}"); lar->DrawLatex(0.60, 0.94, "L = 4.327 fb^{-1} at #sqrt{s} = 13 TeV"); +*/ c1->cd(); c1->SaveAs((path + title + ".pdf").c_str()); diff --git a/XSecsLumi.h b/XSecsLumi.h index b857be1..645ca64 100644 --- a/XSecsLumi.h +++ b/XSecsLumi.h @@ -11,12 +11,14 @@ double Xsec_dy_amc1 = 147.40 *1.23; //////////////////////// DY MadGraph_aMC@NLO + Pythia 8 HT200to400 -double Ngen_dy_amc2 = 8348679; +//double Ngen_dy_amc2 = 8348679; + double Ngen_dy_amc2 = 9628177; double Xsec_dy_amc2 = 41*1.23; //////////////////////// DY MadGraph_aMC@NLO + Pythia 8 HT400to600 -double Ngen_dy_amc3 = 1070450; +//double Ngen_dy_amc3 = 1070450; + double Ngen_dy_amc3 = 1070454; double Xsec_dy_amc3 = 5.678*1.23; //////////////////////// DY MadGraph_aMC@NLO + Pythia 8 HT600toInf diff --git a/globalRun.sh b/globalRun.sh new file mode 100755 index 0000000..4fd6f24 --- /dev/null +++ b/globalRun.sh @@ -0,0 +1,167 @@ +#!/bin/sh + +echo "start" +date + +# -- decide run electron or muon variable here ---- + +channel=ele +#channel=mu + +echo "channel=" $channel + +# which sample you want to run? 1 for True and 0 for False + +flag_DATA_ele=0 +flag_DATA_mu=0 + +flag_MC_DY_100=1 +flag_MC_DY_200=1 +flag_MC_DY_400=1 +flag_MC_DY_600=1 + +flag_MC_TT=0 + +flag_MC_diboson=0 + +echo "which sample are used?" +echo "flag_DATA_ele: $flag_DATA_ele , flag_DATA_mu: $flag_DATA_mu, flag_MC_DY_100: $flag_MC_DY_100 , flag_MC_DY_200: $flag_MC_DY_200 , flag_MC_DY_400: $flag_MC_DY_400 , flag_MC_DY_600: $flag_MC_DY_600 , flag_MC_TT: $flag_MC_TT , flag_MC_diboson: $flag_MC_diboson " + +# ---------------------------------- + +if [ $channel == ele ];then + macro=xAna_ele + outputFolder=output_ele +fi + +if [ $channel == mu ];then + macro=xAna_mu + outputFolder=output_mu +fi + + +# ---------------------------------- +pcncu_path=/data7/yuchang/NCUGlobalTuples_80X + +# ---------------------------------- + + + +function Macro(){ + + # $1=macro, $2=sample_name, $3=save_name + + sample_path=$pcncu_path/$2 + + echo "root -q -b -l $1.C+\(\"$sample_path\"\,\"$outputFolder\"\,\"$3\"\) " + root -q -b -l $1.C+\(\"$sample_path\"\,\"$outputFolder\"\,\"$3\"\) +} + +# ---------------------------------- + +if [ $flag_DATA_ele == 1 ];then + + # SingleElectron + + sample_name=SingleElectron + save_name=SingleElectron-Run2016B-v2 + Macro $macro $sample_name $save_name + +fi + +if [ $flag_DATA_mu == 1 ];then + + # SingleMuon + + sample_name=SingleMuon + save_name=SingleMuon-Run2016B-v2 + Macro $macro $sample_name $save_name + +fi + +# ---------------------------------- + +if [ $flag_MC_DY_100 == 1 ];then + + # DY-Ht_100to200 + sample_name=DYJetsToLL_M-50_HT-100to200_TuneCUETP8M1_13TeV-madgraphMLM-pythia8 + save_name=DYJetsToLL_M-50_HT-100to200_13TeV + Macro $macro $sample_name $save_name +fi + +if [ $flag_MC_DY_200 == 1 ];then + # DY-Ht_200to400 + sample_name=DYJetsToLL_M-50_HT-200to400_TuneCUETP8M1_13TeV-madgraphMLM-pythia8 + save_name=DYJetsToLL_M-50_HT-200to400_13TeV + Macro $macro $sample_name $save_name +fi + +if [ $flag_MC_DY_400 == 1 ];then + # DY-Ht_400to600 + sample_name=DYJetsToLL_M-50_HT-400to600_TuneCUETP8M1_13TeV-madgraphMLM-pythia8 + save_name=DYJetsToLL_M-50_HT-400to600_13TeV + Macro $macro $sample_name $save_name +fi + + +if [ $flag_MC_DY_600 == 1 ];then + # DY-Ht_600toInf + sample_name=DYJetsToLL_M-50_HT-600toInf_TuneCUETP8M1_13TeV-madgraphMLM-pythia8 + save_name=DYJetsToLL_M-50_HT-600toInf_13TeV + Macro $macro $sample_name $save_name + +fi + +# ---------------------------------- + +if [ $flag_MC_TT == 1 ];then + + # TT + sample_name=TT_TuneCUETP8M1_13TeV-powheg-pythia8 + save_name=TT_TuneCUETP8M1_13TeV + Macro $macro $sample_name $save_name + +fi + +# ---------------------------------- + +if [ $flag_MC_diboson == 1 ];then + + # ZH + sample_name=ZH_HToBB_ZToLL_M125_13TeV_amcatnloFXFX_madspin_pythia8 + save_name=ZH_HToBB_ZToLL_M125_13TeV_amcatnlo + Macro $macro $sample_name $save_name + + sample_name=ZH_HToBB_ZToLL_M125_13TeV_powheg_pythia8 + save_name=ZH_HToBB_ZToLL_M125_13TeV_powheg + Macro $macro $sample_name $save_name + + # WW + sample_name=WW_TuneCUETP8M1_13TeV-pythia8 + save_name=WW_TuneCUETP8M1_13TeV + Macro $macro $sample_name $save_name + + # WZ + sample_name=WZ_TuneCUETP8M1_13TeV-pythia8 + save_name=WZ_TuneCUETP8M1_13TeV + Macro $macro $sample_name $save_name + + # ZZ + sample_name=ZZ_TuneCUETP8M1_13TeV-pythia8 + save_name=ZZ_TuneCUETP8M1_13TeV + Macro $macro $sample_name $save_name + +fi + +# ---------------------------------- + +# delete temp file and end the code +rm -f inputdir.txt +rm -f *.pcm *.d *.so + +echo "finish" +date + +exit +# -------- end ---------- + diff --git a/interface/dataFilter.h b/interface/dataFilter.h new file mode 100644 index 0000000..74ec8f7 --- /dev/null +++ b/interface/dataFilter.h @@ -0,0 +1,54 @@ +#include +#include +#include +#include "untuplizer.h" + +bool TriggerStatus(TreeReader &data, std::string TRIGGERNAME){ + + std::string* trigName = data.GetPtrString("hlt_trigName"); + vector& trigResult = *((vector*) data.GetPtr("hlt_trigResult")); + + bool triggerStat = false; + + for(int it = 0; it < data.GetPtrStringSize(); it++){ + + std::string thisTrig = trigName[it]; + bool results = trigResult[it]; + + if( thisTrig.find(TRIGGERNAME) != std::string::npos && results ){ + + triggerStat = true; + break; + + } + + } + + return triggerStat; + +} + +bool FilterStatus(TreeReader &data, std::string FILTERNAME){ + + std::string* filterName = data.GetPtrString("hlt_filterName"); + vector& filterResult = *((vector*) data.GetPtr("hlt_filterResult")); + + bool filterStat = false; + + for(int it = 0; it < data.GetPtrStringSize(); it++){ + + std::string thisFilter = filterName[it]; + bool results = filterResult[it]; + + if( thisFilter.find(FILTERNAME) != std::string::npos && results ){ + + filterStat = true; + break; + + } + + } + + return filterStat; + +} diff --git a/interface/ele_PU/standalone_LumiReWeighting.cc b/interface/ele_PU/standalone_LumiReWeighting.cc new file mode 100644 index 0000000..bc8eafb --- /dev/null +++ b/interface/ele_PU/standalone_LumiReWeighting.cc @@ -0,0 +1,387 @@ +/* information for how I get pileup weight + + take MC's pu_nTrueInt distribution from + https://github.com/cms-sw/cmssw/blob/CMSSW_8_0_X/SimGeneral/MixingModule/python/mix_2016_25ns_SpringMC_PUScenarioV1_PoissonOOTPU_cfi.py + + data's pu_nTrueInt distribution are computed with + My Json file: /afs/cern.ch/work/y/yuchang/produce_samples/CMSSW_8_0_8/src/DelPanj/CrabUtilities/crab_20160419/crab_SingleElectron-Run2016B-PromptReco-v2/results/processedLumis.json + + INPUTLUMIJSON file : /afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions16/13TeV/PileUp/pileup_latest.txt + minBiasXsec : 71300 + + see more in https://twiki.cern.ch/twiki/bin/viewauth/CMS/ExoDiBosonResonancesRun2#PU_weights + + finaly get PU weight in standalone_LumiReWeighting.cc + see more in https://github.com/syuvivida/xtozh_common/tree/76X_analysis/macro_examples/pileup_reweight + +*/ + +/** + \class standalone_LumiReWeighting standalone_LumiReWeighting.h "PhysicsTools/Utilities/interface/standalone_LumiReWeighting.h" + \brief Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data + + This class will trivially take two histograms: + 1. The generated "flat-to-N" distributions from a given processing (or any other generated input) + 2. A histogram generated from the "estimatePileup" macro here: + + https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup + + and produce weights to convert the input distribution (1) to the latter (2). + + \author Shin-Shan Eiko Yu, Salvatore Rappoccio, modified by Mike Hildreth + +*/ +#ifndef standalone_LumiReWeighting_cxx +#define standalone_LumiReWeighting_cxx +#include "TH1.h" +#include "TFile.h" +#include +#include "standalone_LumiReWeighting.h" + +//======================================================= +// For 2012 Data and MC +//======================================================= + + +Double_t Summer2012_S10[60] = { + 0.000829312873542, + 0.00124276120498, + 0.00339329181587, + 0.00408224735376, + 0.00383036590008, + 0.00659159288946, + 0.00816022734493, + 0.00943640833116, + 0.0137777376066, + 0.017059392038, + 0.0213193035468, + 0.0247343174676, + 0.0280848773878, + 0.0323308476564, + 0.0370394341409, + 0.0456917721191, + 0.0558762890594, + 0.0576956187107, + 0.0625325287017, + 0.0591603758776, + 0.0656650815128, + 0.0678329011676, + 0.0625142146389, + 0.0548068448797, + 0.0503893295063, + 0.040209818868, + 0.0374446988111, + 0.0299661572042, + 0.0272024759921, + 0.0219328403791, + 0.0179586571619, + 0.0142926728247, + 0.00839941654725, + 0.00522366397213, + 0.00224457976761, + 0.000779274977993, + 0.000197066585944, + 7.16031761328e-05, + }; + +double Data2012[60]={ +1363.19, +24112.4, +160204, +381537, +714502, +956099, +1.18106e+06, +3.536e+06, +1.48677e+07, +2.83833e+07, +4.19268e+07, +6.41916e+07, +9.44924e+07, +1.36196e+08, +1.94289e+08, +2.62877e+08, +3.2434e+08, +3.64754e+08, +3.84137e+08, +3.86116e+08, +3.72314e+08, +3.44915e+08, +3.07339e+08, +2.63746e+08, +2.18194e+08, +1.73797e+08, +1.32715e+08, +9.65661e+07, +6.65305e+07, +4.31757e+07, +2.6302e+07, +1.50233e+07, +8.05577e+06, +4.06978e+06, +1.9485e+06, +891383, +394079, +171222, +75165.2, +34932.4, +18357.1, +11523.6, +8626.26, +7313.52, +6655.46, +6289.87, +6073.54, +5944.48, +5869.7, +5826.94, +5798.89, +5771.83, + +}; + + +double Data2012Up[60]={ +820.261, +20621.1, +135995, +297808, +645034, +834386, +1.04507e+06, +1.90678e+06, +8.64011e+06, +2.17854e+07, +3.26788e+07, +4.84338e+07, +7.18593e+07, +1.02528e+08, +1.45086e+08, +2.01703e+08, +2.64158e+08, +3.16932e+08, +3.50622e+08, +3.66567e+08, +3.67644e+08, +3.5511e+08, +3.30759e+08, +2.9737e+08, +2.58341e+08, +2.1709e+08, +1.76352e+08, +1.38051e+08, +1.03621e+08, +7.41648e+07, +5.03637e+07, +3.23302e+07, +1.958e+07, +1.11873e+07, +6.04232e+06, +3.09687e+06, +1.51486e+06, +712727, +325989, +147260, +67417, +32628, +17684.7, +11261, +8432.19, +7113.22, +6440.98, +6063.63, +5836.9, +5697.82, +5614.14, +5565.17, + +}; + +double Data2012Down[60]={ +2144.18, +28780.9, +187588, +487719, +794113, +1.08132e+06, +1.53711e+06, +7.18794e+06, +2.25803e+07, +3.61842e+07, +5.59784e+07, +8.58628e+07, +1.26612e+08, +1.85178e+08, +2.59488e+08, +3.30831e+08, +3.79617e+08, +4.0336e+08, +4.06527e+08, +3.91256e+08, +3.60234e+08, +3.17715e+08, +2.68831e+08, +2.18411e+08, +1.70002e+08, +1.2605e+08, +8.83804e+07, +5.81902e+07, +3.57878e+07, +2.05019e+07, +1.09415e+07, +5.45622e+06, +2.55726e+06, +1.13627e+06, +484466, +201789, +84588.1, +37579.2, +19083.5, +11801.9, +8839.96, +7537.79, +6895.75, +6543.22, +6338.84, +6221.11, +6155.62, +6117.94, +6088.98, +6053.17, +5994.32, +5909.29, + +}; + + + +standalone_LumiReWeighting::standalone_LumiReWeighting(int mode) { + +// std::cout << "=======================================================================" << std::endl; + + std::vector MC_distr; + std::vector Lumi_distr; + + MC_distr.clear(); + Lumi_distr.clear(); + switch (mode) + { + case 0: +// std::cout << "Using central value " << std::endl; + break; + case 1: + std::cout << "Using +1 sigma 5% value " << std::endl; + break; + case -1: + std::cout << "Using -1 sigma 5% value " << std::endl; + break; + default: + std::cout << "Using central value " << std::endl; + break; + } // end of switch + + Int_t NBins = 60; + + for( int i=0; i< NBins; ++i) { + switch (mode){ + case 0: + Lumi_distr.push_back(Data2012[i]); + break; + case 1: + Lumi_distr.push_back(Data2012Up[i]); + break; + case -1: + Lumi_distr.push_back(Data2012Down[i]); + break; + default: + Lumi_distr.push_back(Data2012[i]); + break; + } // end of switch + + MC_distr.push_back(Summer2012_S10[i]); + } // end of loop over bins + + // no histograms for input: use vectors + + // now, make histograms out of them: + + // first, check they are the same size... + + if( MC_distr.size() != Lumi_distr.size() ){ + std::cout << "MC_distr.size() = " << MC_distr.size() << std::endl; + std::cout << "Lumi_distr.size() = " << Lumi_distr.size() << std::endl; + std::cerr <<"ERROR: standalone_LumiReWeighting: input vectors have different sizes. Quitting... \n"; + + } + + + weights_ = new TH1D(Form("luminumer_%d",mode), + Form("luminumer_%d",mode), + NBins,0.0, double(NBins)); + + TH1D* den = new TH1D(Form("lumidenom_%d",mode), + Form("lumidenom_%d",mode), + NBins,0.0, double(NBins)); + + + + for(int ibin = 1; ibinSetBinContent(ibin, Lumi_distr[ibin-1]); + den->SetBinContent(ibin,MC_distr[ibin-1]); + } + +/* + std::cout << "Data Input " << std::endl; + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } + std::cout << "MC Input " << std::endl; + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } +*/ + + // check integrals, make sure things are normalized + + double deltaH = weights_->Integral(); + if(fabs(1.0 - deltaH) > 0.02 ) { //*OOPS*... + weights_->Scale( 1.0/ weights_->Integral() ); + } + double deltaMC = den->Integral(); + if(fabs(1.0 - deltaMC) > 0.02 ) { + den->Scale(1.0/ den->Integral()); + } + + weights_->Divide( den ); // so now the average weight should be 1.0 + +/* + std::cout << "Reweighting: Computed Weights per In-Time Nint " << std::endl; + + + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } + + std::cout << "=======================================================================" << std::endl; +*/ + + delete den; +} + +standalone_LumiReWeighting::~standalone_LumiReWeighting() +{ +} + + + +double standalone_LumiReWeighting::weight( double npv ) { + int bin = weights_->GetXaxis()->FindBin( npv ); + double weight_value = weights_->GetBinContent( bin ); + delete weights_; + return weight_value; + +// return weights_->GetBinContent( bin ); +} + + +#endif diff --git a/interface/ele_PU/standalone_LumiReWeighting.h b/interface/ele_PU/standalone_LumiReWeighting.h new file mode 100644 index 0000000..380436b --- /dev/null +++ b/interface/ele_PU/standalone_LumiReWeighting.h @@ -0,0 +1,44 @@ +#ifndef standalone_LumiReWeighting_h +#define standalone_LumiReWeighting_h + + +/** + \class standalone_LumiReWeighting standalone_LumiReWeighting.h "PhysicsTools/Utilities/interface/standalone_LumiReWeighting.h" + \brief Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data + + This class will trivially take two histograms: + 1. The generated "flat-to-N" distributions from a given processing + 2. A histogram generated from the "estimatePileup" macro here: + + https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup + + \author Salvatore Rappoccio +*/ + +#include "TH1.h" +#include "TFile.h" +#include +#include +#include +#include +#include + +using namespace std; + +class standalone_LumiReWeighting { + public: + + standalone_LumiReWeighting(int mode=0); // 0: central, -1: down, +1: up + virtual ~standalone_LumiReWeighting(); + double weight( double npv) ; + + protected: + + TH1D* weights_; + + +}; + + +#endif + diff --git a/interface/ele_PU_69200/MyDataPileupHistogram_true_up.root b/interface/ele_PU_69200/MyDataPileupHistogram_true_up.root new file mode 100644 index 0000000..fbd8171 Binary files /dev/null and b/interface/ele_PU_69200/MyDataPileupHistogram_true_up.root differ diff --git a/interface/ele_PU_69200/standalone_LumiReWeighting.cc b/interface/ele_PU_69200/standalone_LumiReWeighting.cc new file mode 100644 index 0000000..1f9d305 --- /dev/null +++ b/interface/ele_PU_69200/standalone_LumiReWeighting.cc @@ -0,0 +1,396 @@ +/* information for how I get pileup weight + + take MC's pu_nTrueInt distribution from + https://github.com/cms-sw/cmssw/blob/CMSSW_8_0_X/SimGeneral/MixingModule/python/mix_2016_25ns_SpringMC_PUScenarioV1_PoissonOOTPU_cfi.py + + data's pu_nTrueInt distribution are computed with + My Json file: /afs/cern.ch/work/y/yuchang/produce_samples/CMSSW_8_0_8/src/DelPanj/CrabUtilities/crab_20160419/crab_SingleElectron-Run2016B-PromptReco-v2/results/processedLumis.json + + INPUTLUMIJSON file : /afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions16/13TeV/PileUp/pileup_latest.txt + minBiasXsec : 71300 + + see more in https://twiki.cern.ch/twiki/bin/viewauth/CMS/ExoDiBosonResonancesRun2#PU_weights + + finaly get PU weight in standalone_LumiReWeighting.cc + see more in https://github.com/syuvivida/xtozh_common/tree/76X_analysis/macro_examples/pileup_reweight + +*/ + +/** + \class standalone_LumiReWeighting standalone_LumiReWeighting.h "PhysicsTools/Utilities/interface/standalone_LumiReWeighting.h" + \brief Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data + + This class will trivially take two histograms: + 1. The generated "flat-to-N" distributions from a given processing (or any other generated input) + 2. A histogram generated from the "estimatePileup" macro here: + + https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup + + and produce weights to convert the input distribution (1) to the latter (2). + + \author Shin-Shan Eiko Yu, Salvatore Rappoccio, modified by Mike Hildreth + +*/ +#ifndef standalone_LumiReWeighting_cxx +#define standalone_LumiReWeighting_cxx +#include "TH1.h" +#include "TFile.h" +#include +#include "standalone_LumiReWeighting.h" + +//======================================================= +// For 2012 Data and MC +//======================================================= + + +Double_t Summer2012_S10[60] = { + 0.000829312873542, + 0.00124276120498, + 0.00339329181587, + 0.00408224735376, + 0.00383036590008, + 0.00659159288946, + 0.00816022734493, + 0.00943640833116, + 0.0137777376066, + 0.017059392038, + 0.0213193035468, + 0.0247343174676, + 0.0280848773878, + 0.0323308476564, + 0.0370394341409, + 0.0456917721191, + 0.0558762890594, + 0.0576956187107, + 0.0625325287017, + 0.0591603758776, + 0.0656650815128, + 0.0678329011676, + 0.0625142146389, + 0.0548068448797, + 0.0503893295063, + 0.040209818868, + 0.0374446988111, + 0.0299661572042, + 0.0272024759921, + 0.0219328403791, + 0.0179586571619, + 0.0142926728247, + 0.00839941654725, + 0.00522366397213, + 0.00224457976761, + 0.000779274977993, + 0.000197066585944, + 7.16031761328e-05, + }; + +// Minimum Bias cross section=69200 +double Data2012[60]={ +1796.51, +26622.1, +175629, +442024, +758763, +1.03048e+06, +1.34376e+06, +5.36593e+06, +1.93122e+07, +3.26895e+07, +4.95164e+07, +7.61144e+07, +1.11918e+08, +1.62854e+08, +2.30866e+08, +3.02722e+08, +3.57945e+08, +3.8878e+08, +3.9898e+08, +3.91232e+08, +3.67465e+08, +3.31032e+08, +2.86392e+08, +2.38277e+08, +1.90564e+08, +1.45936e+08, +1.0634e+08, +7.3219e+07, +4.73519e+07, +2.86461e+07, +1.61853e+07, +8.55038e+06, +4.23926e+06, +1.9851e+06, +885854, +381481, +161561, +69468.9, +32009.9, +17004.3, +10963, +8438, +7294.05, +6714.37, +6389.72, +6199.02, +6088, +6025.52, +5989.21, +5961.48, +5928.1, +5873.67, +0, +0, +0, +0, +0, +0, +0, +0, + +}; + +// Up and Down I haven't use for 69200 +double Data2012Up[60]={ +820.261, +20621.1, +135995, +297808, +645034, +834386, +1.04507e+06, +1.90678e+06, +8.64011e+06, +2.17854e+07, +3.26788e+07, +4.84338e+07, +7.18593e+07, +1.02528e+08, +1.45086e+08, +2.01703e+08, +2.64158e+08, +3.16932e+08, +3.50622e+08, +3.66567e+08, +3.67644e+08, +3.5511e+08, +3.30759e+08, +2.9737e+08, +2.58341e+08, +2.1709e+08, +1.76352e+08, +1.38051e+08, +1.03621e+08, +7.41648e+07, +5.03637e+07, +3.23302e+07, +1.958e+07, +1.11873e+07, +6.04232e+06, +3.09687e+06, +1.51486e+06, +712727, +325989, +147260, +67417, +32628, +17684.7, +11261, +8432.19, +7113.22, +6440.98, +6063.63, +5836.9, +5697.82, +5614.14, +5565.17, + +}; + +double Data2012Down[60]={ +2144.18, +28780.9, +187588, +487719, +794113, +1.08132e+06, +1.53711e+06, +7.18794e+06, +2.25803e+07, +3.61842e+07, +5.59784e+07, +8.58628e+07, +1.26612e+08, +1.85178e+08, +2.59488e+08, +3.30831e+08, +3.79617e+08, +4.0336e+08, +4.06527e+08, +3.91256e+08, +3.60234e+08, +3.17715e+08, +2.68831e+08, +2.18411e+08, +1.70002e+08, +1.2605e+08, +8.83804e+07, +5.81902e+07, +3.57878e+07, +2.05019e+07, +1.09415e+07, +5.45622e+06, +2.55726e+06, +1.13627e+06, +484466, +201789, +84588.1, +37579.2, +19083.5, +11801.9, +8839.96, +7537.79, +6895.75, +6543.22, +6338.84, +6221.11, +6155.62, +6117.94, +6088.98, +6053.17, +5994.32, +5909.29, + +}; + + + +standalone_LumiReWeighting::standalone_LumiReWeighting(int mode) { + +// std::cout << "=======================================================================" << std::endl; + + std::vector MC_distr; + std::vector Lumi_distr; + + MC_distr.clear(); + Lumi_distr.clear(); + switch (mode) + { + case 0: +// std::cout << "Using central value " << std::endl; + break; + case 1: + std::cout << "Using +1 sigma 5% value " << std::endl; + break; + case -1: + std::cout << "Using -1 sigma 5% value " << std::endl; + break; + default: + std::cout << "Using central value " << std::endl; + break; + } // end of switch + + Int_t NBins = 60; + + for( int i=0; i< NBins; ++i) { + switch (mode){ + case 0: + Lumi_distr.push_back(Data2012[i]); + break; + case 1: + Lumi_distr.push_back(Data2012Up[i]); + break; + case -1: + Lumi_distr.push_back(Data2012Down[i]); + break; + default: + Lumi_distr.push_back(Data2012[i]); + break; + } // end of switch + + MC_distr.push_back(Summer2012_S10[i]); + } // end of loop over bins + + // no histograms for input: use vectors + + // now, make histograms out of them: + + // first, check they are the same size... + + if( MC_distr.size() != Lumi_distr.size() ){ + std::cout << "MC_distr.size() = " << MC_distr.size() << std::endl; + std::cout << "Lumi_distr.size() = " << Lumi_distr.size() << std::endl; + std::cerr <<"ERROR: standalone_LumiReWeighting: input vectors have different sizes. Quitting... \n"; + + } + + + weights_ = new TH1D(Form("luminumer_%d",mode), + Form("luminumer_%d",mode), + NBins,0.0, double(NBins)); + + TH1D* den = new TH1D(Form("lumidenom_%d",mode), + Form("lumidenom_%d",mode), + NBins,0.0, double(NBins)); + + + + for(int ibin = 1; ibinSetBinContent(ibin, Lumi_distr[ibin-1]); + den->SetBinContent(ibin,MC_distr[ibin-1]); + } + +/* + std::cout << "Data Input " << std::endl; + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } + std::cout << "MC Input " << std::endl; + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } +*/ + + // check integrals, make sure things are normalized + + double deltaH = weights_->Integral(); + if(fabs(1.0 - deltaH) > 0.02 ) { //*OOPS*... + weights_->Scale( 1.0/ weights_->Integral() ); + } + double deltaMC = den->Integral(); + if(fabs(1.0 - deltaMC) > 0.02 ) { + den->Scale(1.0/ den->Integral()); + } + + weights_->Divide( den ); // so now the average weight should be 1.0 + +/* + std::cout << "Reweighting: Computed Weights per In-Time Nint " << std::endl; + + + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } + + std::cout << "=======================================================================" << std::endl; +*/ + + delete den; +} + +standalone_LumiReWeighting::~standalone_LumiReWeighting() +{ +} + + + +double standalone_LumiReWeighting::weight( double npv ) { + int bin = weights_->GetXaxis()->FindBin( npv ); + double weight_value = weights_->GetBinContent( bin ); + delete weights_; + return weight_value; + +// return weights_->GetBinContent( bin ); +} + + +#endif diff --git a/interface/ele_PU_69200/standalone_LumiReWeighting.h b/interface/ele_PU_69200/standalone_LumiReWeighting.h new file mode 100644 index 0000000..380436b --- /dev/null +++ b/interface/ele_PU_69200/standalone_LumiReWeighting.h @@ -0,0 +1,44 @@ +#ifndef standalone_LumiReWeighting_h +#define standalone_LumiReWeighting_h + + +/** + \class standalone_LumiReWeighting standalone_LumiReWeighting.h "PhysicsTools/Utilities/interface/standalone_LumiReWeighting.h" + \brief Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data + + This class will trivially take two histograms: + 1. The generated "flat-to-N" distributions from a given processing + 2. A histogram generated from the "estimatePileup" macro here: + + https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup + + \author Salvatore Rappoccio +*/ + +#include "TH1.h" +#include "TFile.h" +#include +#include +#include +#include +#include + +using namespace std; + +class standalone_LumiReWeighting { + public: + + standalone_LumiReWeighting(int mode=0); // 0: central, -1: down, +1: up + virtual ~standalone_LumiReWeighting(); + double weight( double npv) ; + + protected: + + TH1D* weights_; + + +}; + + +#endif + diff --git a/interface/ele_PU_71300/standalone_LumiReWeighting.cc b/interface/ele_PU_71300/standalone_LumiReWeighting.cc new file mode 100644 index 0000000..bc8eafb --- /dev/null +++ b/interface/ele_PU_71300/standalone_LumiReWeighting.cc @@ -0,0 +1,387 @@ +/* information for how I get pileup weight + + take MC's pu_nTrueInt distribution from + https://github.com/cms-sw/cmssw/blob/CMSSW_8_0_X/SimGeneral/MixingModule/python/mix_2016_25ns_SpringMC_PUScenarioV1_PoissonOOTPU_cfi.py + + data's pu_nTrueInt distribution are computed with + My Json file: /afs/cern.ch/work/y/yuchang/produce_samples/CMSSW_8_0_8/src/DelPanj/CrabUtilities/crab_20160419/crab_SingleElectron-Run2016B-PromptReco-v2/results/processedLumis.json + + INPUTLUMIJSON file : /afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions16/13TeV/PileUp/pileup_latest.txt + minBiasXsec : 71300 + + see more in https://twiki.cern.ch/twiki/bin/viewauth/CMS/ExoDiBosonResonancesRun2#PU_weights + + finaly get PU weight in standalone_LumiReWeighting.cc + see more in https://github.com/syuvivida/xtozh_common/tree/76X_analysis/macro_examples/pileup_reweight + +*/ + +/** + \class standalone_LumiReWeighting standalone_LumiReWeighting.h "PhysicsTools/Utilities/interface/standalone_LumiReWeighting.h" + \brief Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data + + This class will trivially take two histograms: + 1. The generated "flat-to-N" distributions from a given processing (or any other generated input) + 2. A histogram generated from the "estimatePileup" macro here: + + https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup + + and produce weights to convert the input distribution (1) to the latter (2). + + \author Shin-Shan Eiko Yu, Salvatore Rappoccio, modified by Mike Hildreth + +*/ +#ifndef standalone_LumiReWeighting_cxx +#define standalone_LumiReWeighting_cxx +#include "TH1.h" +#include "TFile.h" +#include +#include "standalone_LumiReWeighting.h" + +//======================================================= +// For 2012 Data and MC +//======================================================= + + +Double_t Summer2012_S10[60] = { + 0.000829312873542, + 0.00124276120498, + 0.00339329181587, + 0.00408224735376, + 0.00383036590008, + 0.00659159288946, + 0.00816022734493, + 0.00943640833116, + 0.0137777376066, + 0.017059392038, + 0.0213193035468, + 0.0247343174676, + 0.0280848773878, + 0.0323308476564, + 0.0370394341409, + 0.0456917721191, + 0.0558762890594, + 0.0576956187107, + 0.0625325287017, + 0.0591603758776, + 0.0656650815128, + 0.0678329011676, + 0.0625142146389, + 0.0548068448797, + 0.0503893295063, + 0.040209818868, + 0.0374446988111, + 0.0299661572042, + 0.0272024759921, + 0.0219328403791, + 0.0179586571619, + 0.0142926728247, + 0.00839941654725, + 0.00522366397213, + 0.00224457976761, + 0.000779274977993, + 0.000197066585944, + 7.16031761328e-05, + }; + +double Data2012[60]={ +1363.19, +24112.4, +160204, +381537, +714502, +956099, +1.18106e+06, +3.536e+06, +1.48677e+07, +2.83833e+07, +4.19268e+07, +6.41916e+07, +9.44924e+07, +1.36196e+08, +1.94289e+08, +2.62877e+08, +3.2434e+08, +3.64754e+08, +3.84137e+08, +3.86116e+08, +3.72314e+08, +3.44915e+08, +3.07339e+08, +2.63746e+08, +2.18194e+08, +1.73797e+08, +1.32715e+08, +9.65661e+07, +6.65305e+07, +4.31757e+07, +2.6302e+07, +1.50233e+07, +8.05577e+06, +4.06978e+06, +1.9485e+06, +891383, +394079, +171222, +75165.2, +34932.4, +18357.1, +11523.6, +8626.26, +7313.52, +6655.46, +6289.87, +6073.54, +5944.48, +5869.7, +5826.94, +5798.89, +5771.83, + +}; + + +double Data2012Up[60]={ +820.261, +20621.1, +135995, +297808, +645034, +834386, +1.04507e+06, +1.90678e+06, +8.64011e+06, +2.17854e+07, +3.26788e+07, +4.84338e+07, +7.18593e+07, +1.02528e+08, +1.45086e+08, +2.01703e+08, +2.64158e+08, +3.16932e+08, +3.50622e+08, +3.66567e+08, +3.67644e+08, +3.5511e+08, +3.30759e+08, +2.9737e+08, +2.58341e+08, +2.1709e+08, +1.76352e+08, +1.38051e+08, +1.03621e+08, +7.41648e+07, +5.03637e+07, +3.23302e+07, +1.958e+07, +1.11873e+07, +6.04232e+06, +3.09687e+06, +1.51486e+06, +712727, +325989, +147260, +67417, +32628, +17684.7, +11261, +8432.19, +7113.22, +6440.98, +6063.63, +5836.9, +5697.82, +5614.14, +5565.17, + +}; + +double Data2012Down[60]={ +2144.18, +28780.9, +187588, +487719, +794113, +1.08132e+06, +1.53711e+06, +7.18794e+06, +2.25803e+07, +3.61842e+07, +5.59784e+07, +8.58628e+07, +1.26612e+08, +1.85178e+08, +2.59488e+08, +3.30831e+08, +3.79617e+08, +4.0336e+08, +4.06527e+08, +3.91256e+08, +3.60234e+08, +3.17715e+08, +2.68831e+08, +2.18411e+08, +1.70002e+08, +1.2605e+08, +8.83804e+07, +5.81902e+07, +3.57878e+07, +2.05019e+07, +1.09415e+07, +5.45622e+06, +2.55726e+06, +1.13627e+06, +484466, +201789, +84588.1, +37579.2, +19083.5, +11801.9, +8839.96, +7537.79, +6895.75, +6543.22, +6338.84, +6221.11, +6155.62, +6117.94, +6088.98, +6053.17, +5994.32, +5909.29, + +}; + + + +standalone_LumiReWeighting::standalone_LumiReWeighting(int mode) { + +// std::cout << "=======================================================================" << std::endl; + + std::vector MC_distr; + std::vector Lumi_distr; + + MC_distr.clear(); + Lumi_distr.clear(); + switch (mode) + { + case 0: +// std::cout << "Using central value " << std::endl; + break; + case 1: + std::cout << "Using +1 sigma 5% value " << std::endl; + break; + case -1: + std::cout << "Using -1 sigma 5% value " << std::endl; + break; + default: + std::cout << "Using central value " << std::endl; + break; + } // end of switch + + Int_t NBins = 60; + + for( int i=0; i< NBins; ++i) { + switch (mode){ + case 0: + Lumi_distr.push_back(Data2012[i]); + break; + case 1: + Lumi_distr.push_back(Data2012Up[i]); + break; + case -1: + Lumi_distr.push_back(Data2012Down[i]); + break; + default: + Lumi_distr.push_back(Data2012[i]); + break; + } // end of switch + + MC_distr.push_back(Summer2012_S10[i]); + } // end of loop over bins + + // no histograms for input: use vectors + + // now, make histograms out of them: + + // first, check they are the same size... + + if( MC_distr.size() != Lumi_distr.size() ){ + std::cout << "MC_distr.size() = " << MC_distr.size() << std::endl; + std::cout << "Lumi_distr.size() = " << Lumi_distr.size() << std::endl; + std::cerr <<"ERROR: standalone_LumiReWeighting: input vectors have different sizes. Quitting... \n"; + + } + + + weights_ = new TH1D(Form("luminumer_%d",mode), + Form("luminumer_%d",mode), + NBins,0.0, double(NBins)); + + TH1D* den = new TH1D(Form("lumidenom_%d",mode), + Form("lumidenom_%d",mode), + NBins,0.0, double(NBins)); + + + + for(int ibin = 1; ibinSetBinContent(ibin, Lumi_distr[ibin-1]); + den->SetBinContent(ibin,MC_distr[ibin-1]); + } + +/* + std::cout << "Data Input " << std::endl; + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } + std::cout << "MC Input " << std::endl; + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } +*/ + + // check integrals, make sure things are normalized + + double deltaH = weights_->Integral(); + if(fabs(1.0 - deltaH) > 0.02 ) { //*OOPS*... + weights_->Scale( 1.0/ weights_->Integral() ); + } + double deltaMC = den->Integral(); + if(fabs(1.0 - deltaMC) > 0.02 ) { + den->Scale(1.0/ den->Integral()); + } + + weights_->Divide( den ); // so now the average weight should be 1.0 + +/* + std::cout << "Reweighting: Computed Weights per In-Time Nint " << std::endl; + + + for(int ibin = 1; ibinGetBinContent(ibin) << std::endl; + } + + std::cout << "=======================================================================" << std::endl; +*/ + + delete den; +} + +standalone_LumiReWeighting::~standalone_LumiReWeighting() +{ +} + + + +double standalone_LumiReWeighting::weight( double npv ) { + int bin = weights_->GetXaxis()->FindBin( npv ); + double weight_value = weights_->GetBinContent( bin ); + delete weights_; + return weight_value; + +// return weights_->GetBinContent( bin ); +} + + +#endif diff --git a/interface/ele_PU_71300/standalone_LumiReWeighting.h b/interface/ele_PU_71300/standalone_LumiReWeighting.h new file mode 100644 index 0000000..380436b --- /dev/null +++ b/interface/ele_PU_71300/standalone_LumiReWeighting.h @@ -0,0 +1,44 @@ +#ifndef standalone_LumiReWeighting_h +#define standalone_LumiReWeighting_h + + +/** + \class standalone_LumiReWeighting standalone_LumiReWeighting.h "PhysicsTools/Utilities/interface/standalone_LumiReWeighting.h" + \brief Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data + + This class will trivially take two histograms: + 1. The generated "flat-to-N" distributions from a given processing + 2. A histogram generated from the "estimatePileup" macro here: + + https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup + + \author Salvatore Rappoccio +*/ + +#include "TH1.h" +#include "TFile.h" +#include +#include +#include +#include +#include + +using namespace std; + +class standalone_LumiReWeighting { + public: + + standalone_LumiReWeighting(int mode=0); // 0: central, -1: down, +1: up + virtual ~standalone_LumiReWeighting(); + double weight( double npv) ; + + protected: + + TH1D* weights_; + + +}; + + +#endif + diff --git a/interface/isPassZee.h b/interface/isPassZee.h new file mode 100644 index 0000000..a87afa9 --- /dev/null +++ b/interface/isPassZee.h @@ -0,0 +1,71 @@ +#include +#include +#include +#include +#include +#include "untuplizer.h" + +bool isPassZee(TreeReader &data, vector& goodEleID){ + + goodEleID.clear(); + + const Int_t nEle = data.GetInt("nEle"); + const Int_t* eleCharge = data.GetPtrInt("eleCharge"); + const Float_t* eleScEta = data.GetPtrFloat("eleScEta"); + const TClonesArray* eleP4 = (TClonesArray*) data.GetPtrTObject("eleP4"); + const vector& eleIsPassLoose = *((vector*) data.GetPtr("eleIsPassLoose")); + + // select good electrons + + std::vector goodElectrons; + + for( int ie = nEle-1; ie >= 0; --ie ){ + + TLorentzVector* myEle = (TLorentzVector*)eleP4->At(ie); + + if( !eleIsPassLoose[ie] ) continue; + if( fabs(eleScEta[ie]) > 1.4442 && fabs(eleScEta[ie]) < 1.566 ) continue; + if( fabs(eleScEta[ie]) > 2.5 ) continue; + if( fabs(myEle->Eta()) > 2.5 ) continue; + if( myEle->Pt() < 20 ) continue; + + goodElectrons.push_back(ie); + + } // end of ele loop + + // select good Z boson + + bool findEPair = false; + TLorentzVector *thisEle = NULL, *thatEle = NULL; + + for( unsigned int i = 0; i < goodElectrons.size(); ++i ){ + + int ie = goodElectrons[i]; + thisEle = (TLorentzVector*)eleP4->At(ie); + + for( unsigned int j = 0; j < i; ++j ){ + + int je = goodElectrons[j]; + thatEle = (TLorentzVector*)eleP4->At(je); + + if( eleCharge[ie]*eleCharge[je] > 0 ) continue; + if( TMath::Max(thisEle->Pt(),thatEle->Pt()) < 135 ) continue; + if( (*thisEle+*thatEle).M() < 70 || (*thisEle+*thatEle).M() > 110 ) continue; + if( (*thisEle+*thatEle).Pt() < 200 ) continue; + + if( !findEPair ){ + + goodEleID.push_back( (thisEle->Pt() > thatEle->Pt()) ? ie : je ); + goodEleID.push_back( (thisEle->Pt() > thatEle->Pt()) ? je : ie ); + + } + + findEPair = true; + break; + + } + } + + return findEPair ? true : false; + +} diff --git a/interface/isPassZmumu.h b/interface/isPassZmumu.h new file mode 100644 index 0000000..d6ab1a4 --- /dev/null +++ b/interface/isPassZmumu.h @@ -0,0 +1,84 @@ +#include +#include +#include +#include +#include "untuplizer.h" + +bool isPassZmumu(TreeReader &data, vector& goodMuID){ + + goodMuID.clear(); + + const Int_t nMu = data.GetInt("nMu"); + const Int_t* muCharge = data.GetPtrInt("muCharge"); + const Float_t* muTrkIso = data.GetPtrFloat("muTrkIso"); + const Float_t* muInnerTrkPt = data.GetPtrFloat("muInnerTrkPt"); + const TClonesArray* muP4 = (TClonesArray*) data.GetPtrTObject("muP4"); + const vector& isHighPtMuon = *((vector*) data.GetPtr("isHighPtMuon")); + const vector& isCustomTrackerMuon = *((vector*) data.GetPtr("isCustomTrackerMuon")); + + // select good muons + + std::vector goodMuons; + bool hasTrigMuon = false; + + for( int im = nMu-1; im >= 0; --im ){ + + TLorentzVector* myMu = (TLorentzVector*)muP4->At(im); + + if( !isHighPtMuon[im] && !isCustomTrackerMuon[im] ) continue; + if( fabs(myMu->Eta()) > 2.4 ) continue; + if( myMu->Pt() < 20 ) continue; + if( myMu->Pt() > 55 && fabs(myMu->Eta()) < 2.4 ) hasTrigMuon = true; + + goodMuons.push_back(im); + + } + + if( !hasTrigMuon ) goodMuons.clear(); + + // select good Z boson + + bool findMPair = false; + TLorentzVector *thisMu = NULL, *thatMu = NULL; + + for( unsigned int i = 0; i < goodMuons.size(); ++i ){ + + int im = goodMuons[i]; + thisMu = (TLorentzVector*)muP4->At(im); + + for( unsigned int j = 0; j < i; ++j ){ + + int jm = goodMuons[j]; + thatMu = (TLorentzVector*)muP4->At(jm); + + // if the two muons are far away, use regular isolation + + if( thisMu->DeltaR(*thatMu) > 0.3 && muTrkIso[im]/thisMu->Pt() > 0.1 ) continue; + if( thisMu->DeltaR(*thatMu) > 0.3 && muTrkIso[jm]/thatMu->Pt() > 0.1 ) continue; + + // if the two muons are close, use corrected isolation + + if( thisMu->DeltaR(*thatMu) < 0.3 && (muTrkIso[im]-muInnerTrkPt[jm])/thisMu->Pt() > 0.1 ) continue; + if( thisMu->DeltaR(*thatMu) < 0.3 && (muTrkIso[jm]-muInnerTrkPt[im])/thatMu->Pt() > 0.1 ) continue; + + if( !isHighPtMuon[im] && !isHighPtMuon[jm] ) continue; + if( muCharge[im]*muCharge[jm] > 0 ) continue; + if( (*thisMu+*thatMu).M() < 70 || (*thisMu+*thatMu).M() > 110 ) continue; + if( (*thisMu+*thatMu).Pt() < 200 ) continue; + + if( !findMPair ){ + + goodMuID.push_back( (thisMu->Pt() > thatMu->Pt()) ? im : jm ); + goodMuID.push_back( (thisMu->Pt() > thatMu->Pt()) ? jm : im ); + + } + + findMPair = true; + break; + + } + } + + return findMPair ? true : false; + +} diff --git a/interface/readHists.h b/interface/readHists.h new file mode 100644 index 0000000..dfd4a86 --- /dev/null +++ b/interface/readHists.h @@ -0,0 +1,71 @@ +#include +#include +#include +#include +#include + +class readHist{ + + public: + + readHist(string); + TH1D* getHist(string); + + private: + + TFile* thisFile; + string thisFileName; + float crossSection(string); + +}; + +readHist::readHist(string rootFile){ + + string tmpStr = rootFile.erase(rootFile.find_last_not_of("/")+1); + thisFileName = tmpStr.substr(tmpStr.find_last_of("/")+1); + thisFile = TFile::Open(rootFile.data()); + +} + +float readHist::crossSection(string thisPath){ + + ifstream textFile("/afs/cern.ch/work/y/yuchang/ZH_analysis_with_2016_data/CMSSW_8_0_8/src/DataMC_comparison/interface/xSec.txt"); +// ifstream textFile("/afs/cern.ch/work/h/htong/ZpZHllbb_13TeV/xSec.txt"); + string token; + float crosssection = 0., thisNum = 0.; + + while( textFile >> token >> thisNum ){ + + if( thisPath.find(token) != string::npos ) + crosssection = thisNum; + + } + + return crosssection; + +} + +TH1D* readHist::getHist(string hname){ + + TH1D* thisHist = (TH1D*)(thisFile->Get(Form("%s", hname.c_str()))); + +// cout<<"here is readHist::getHist"<< endl; +// cout<<"thisFileName: "<< thisFileName << endl; + +// cout<<"histo open name: "<Get(totalEvents)))->Integral(): "<< ((TH1D*)(thisFile->Get("totalEvents")))->Integral()<< endl; + +// thisHist->Scale(thisFileName.find("Run2016") != string::npos ? 1 : 2512.*crossSection(thisFileName.data())/((TH1D*)(thisFile->Get("totalEvents")))->Integral()); + +//double lumi = 4327.13; +//double lumi = 1/178*4327.13;// approx for 1M events +double lumi = 432.713; + + thisHist->Scale(thisFileName.find("Run2016") != string::npos ? 1 : lumi*crossSection(thisFileName.data())/((TH1D*)(thisFile->Get("totalEvents")))->Integral()); + + return thisHist; + +} diff --git a/interface/readSample.h b/interface/readSample.h new file mode 100644 index 0000000..93c491f --- /dev/null +++ b/interface/readSample.h @@ -0,0 +1,46 @@ +#include +#include +#include +#include +#include +#include +#include + +//void readSample(std::string inputFile, std::vector& infiles){ +void readSample(TString inputFile, std::vector& infiles){ + + string inputTextFile = "inputdir.txt"; + gSystem->Exec(Form("rm -rf %s",inputTextFile.data())); +// gSystem->Exec(Form("ls -R %s | grep -a \"%s\" >> %s", inputFile.data(),"data7", inputTextFile.data())); + gSystem->Exec(Form("ls -R %s | grep -a \"%s\" >> %s", inputFile.Data(),"data7", inputTextFile.data())); + TSystemDirectory *base = new TSystemDirectory("root","root"); + int nfile = 0; + string tempdir; + ifstream fin; + fin.open(inputTextFile.data()); + fin >> tempdir; + TString realDirName = gSystem->GetFromPipe(Form("a=%s; echo ${a%%:*}",tempdir.data())); + while(!fin.eof()){ + if(realDirName.Contains("fail")){ + fin >> tempdir; + realDirName = gSystem->GetFromPipe(Form("a=%s; echo ${a%%:*}",tempdir.data())); + continue; + } + base->SetDirectory(realDirName.Data()); + TList *listOfFiles = base->GetListOfFiles(); + TIter fileIt(listOfFiles); + TFile *fileH = new TFile(); + while((fileH = (TFile*)fileIt())) { + std::string fileN = fileH->GetName(); + std::string baseString = "root"; + if( fileH->IsFolder()) continue; + if(fileN.find(baseString) == std::string::npos)continue; + nfile++; + string tempfile = Form("%s/%s",realDirName.Data(),fileN.data()); + infiles.push_back(tempfile); + } + fin >> tempdir; + realDirName = gSystem->GetFromPipe(Form("a=%s; echo ${a%%:*}",tempdir.data())); + } + +} diff --git a/interface/setNCUStyle.h b/interface/setNCUStyle.h new file mode 100644 index 0000000..c1a18ba --- /dev/null +++ b/interface/setNCUStyle.h @@ -0,0 +1,217 @@ +#include "TStyle.h" +#include "TPad.h" + +// fixOverlay: Redraws the axis +// void fixOverlay() { +// gPad->RedrawAxis(); +// } + +void setNCUStyle(bool gridOn=false) { + TStyle* ncuStyle = new TStyle("ncuStyle","Style for P-NCU"); + + // For the canvas: + ncuStyle->SetCanvasBorderMode(0); + ncuStyle->SetCanvasColor(kWhite); + ncuStyle->SetCanvasDefH(600); //Height of canvas + ncuStyle->SetCanvasDefW(600); //Width of canvas + ncuStyle->SetCanvasDefX(0); //POsition on screen + ncuStyle->SetCanvasDefY(0); + + // For the Pad: + ncuStyle->SetPadBorderMode(0); + // ncuStyle->SetPadBorderSize(Width_t size = 1); + ncuStyle->SetPadColor(kWhite); + ncuStyle->SetPadGridX(gridOn); + ncuStyle->SetPadGridY(gridOn); + ncuStyle->SetGridColor(0); + ncuStyle->SetGridStyle(3); + ncuStyle->SetGridWidth(1); + + // For the frame: + ncuStyle->SetFrameBorderMode(0); + ncuStyle->SetFrameBorderSize(1); + ncuStyle->SetFrameFillColor(0); + ncuStyle->SetFrameFillStyle(0); + ncuStyle->SetFrameLineColor(1); + ncuStyle->SetFrameLineStyle(1); + ncuStyle->SetFrameLineWidth(3); + + // For the Legend: + ncuStyle->SetLegendBorderSize(0); + ncuStyle->SetLegendFillColor(0); + ncuStyle->SetLegendFont(42); + //ncuStyle->SetLegendFont(62); + + // For the histo: + // ncuStyle->SetHistFillColor(1); + // ncuStyle->SetHistFillStyle(0); + ncuStyle->SetHistLineColor(1); + ncuStyle->SetHistLineStyle(0); + ncuStyle->SetHistLineWidth(3); + // ncuStyle->SetLegoInnerR(Float_t rad = 0.5); + // ncuStyle->SetNumberContours(Int_t number = 20); + + ncuStyle->SetEndErrorSize(2); + // ncuStyle->SetErrorMarker(20); + //ncuStyle->SetErrorX(0.); + + ncuStyle->SetMarkerStyle(20); + + //For the fit/function: + ncuStyle->SetOptFit(1); + ncuStyle->SetFitFormat("5.4g"); + ncuStyle->SetFuncColor(2); + ncuStyle->SetFuncStyle(1); + ncuStyle->SetFuncWidth(1); + + //For the date: + ncuStyle->SetOptDate(0); + // ncuStyle->SetDateX(Float_t x = 0.01); + // ncuStyle->SetDateY(Float_t y = 0.01); + + // For the statistics box: + ncuStyle->SetOptFile(0); + ncuStyle->SetOptStat(0); // To display the mean and RMS: SetOptStat("mr"); + ncuStyle->SetStatColor(kWhite); + ncuStyle->SetStatFont(42); + ncuStyle->SetStatFontSize(0.025); + ncuStyle->SetStatTextColor(1); + ncuStyle->SetStatFormat("6.4g"); + ncuStyle->SetStatBorderSize(1); + ncuStyle->SetStatH(0.1); + ncuStyle->SetStatW(0.15); + // ncuStyle->SetStatStyle(Style_t style = 1001); + // ncuStyle->SetStatX(Float_t x = 0); + // ncuStyle->SetStatY(Float_t y = 0); + + // Margins: + ncuStyle->SetPadTopMargin(0.08); + ncuStyle->SetPadBottomMargin(0.13); + ncuStyle->SetPadLeftMargin(0.13); + ncuStyle->SetPadRightMargin(0.04); + + // For the Global title: + + ncuStyle->SetOptTitle(0); + ncuStyle->SetTitleFont(62); + ncuStyle->SetTitleColor(1); + ncuStyle->SetTitleTextColor(1); + ncuStyle->SetTitleFillColor(10); + ncuStyle->SetTitleFontSize(0.05); + // ncuStyle->SetTitleH(0); // Set the height of the title box + // ncuStyle->SetTitleW(0); // Set the width of the title box + // ncuStyle->SetTitleX(0); // Set the position of the title box + // ncuStyle->SetTitleY(0.985); // Set the position of the title box + // ncuStyle->SetTitleStyle(Style_t style = 1001); + // ncuStyle->SetTitleBorderSize(2); + + // For the axis titles: + + ncuStyle->SetTitleColor(1, "XYZ"); + ncuStyle->SetTitleSize(0.06, "XYZ"); + // ncuStyle->SetTitleXSize(Float_t size = 0.02); // Another way to set the size? + // ncuStyle->SetTitleYSize(Float_t size = 0.02); + // the following commands are not doing any thing + // ncuStyle->SetTitleOffset(1.5, "X"); // Another way to set the Offset + // ncuStyle->SetTitleOffset(1.5, "Y"); // Another way to set the Offset + // ncuStyle->SetTitleFont(42, "XYZ"); + + // For the axis labels: + ncuStyle->SetLabelColor(1, "XYZ"); + // the following command is not doing any thing + // ncuStyle->SetLabelFont(42, "XYZ"); + ncuStyle->SetLabelOffset(0.007, "XYZ"); + ncuStyle->SetLabelSize(0.05, "XYZ"); + + // For the axis: + + ncuStyle->SetAxisColor(1, "XYZ"); + ncuStyle->SetStripDecimals(kFALSE); + ncuStyle->SetTickLength(0.03, "XYZ"); + ncuStyle->SetNdivisions(510, "XYZ"); + ncuStyle->SetPadTickX(1); // To get tick marks on the opposite side of the frame + ncuStyle->SetPadTickY(1); + + // Change for log plots: + ncuStyle->SetOptLogx(0); + ncuStyle->SetOptLogy(0); + ncuStyle->SetOptLogz(0); + + // Postscript options: + ncuStyle->SetPaperSize(20.,20.); + // ncuStyle->SetLineScalePS(Float_t scale = 3); + // ncuStyle->SetLineStyleString(Int_t i, const char* text); + // ncuStyle->SetHeaderPS(const char* header); + // ncuStyle->SetTitlePS(const char* pstitle); + + // ncuStyle->SetBarOffset(Float_t baroff = 0.5); + // ncuStyle->SetBarWidth(Float_t barwidth = 0.5); + // ncuStyle->SetPaintTextFormat(const char* format = "g"); + // ncuStyle->SetPalette(Int_t ncolors = 0, Int_t* colors = 0); + // ncuStyle->SetTimeOffset(Double_t toffset); + // ncuStyle->SetHistMinimumZero(kTRUE); + + ncuStyle->SetHatchesLineWidth(5); + ncuStyle->SetHatchesSpacing(0.05); + + ncuStyle->cd(); + +} + +/* + void additional_style(){ + + // For Legend + TLegend* legend = new TLegend(0.5, 0.7, 0.9,0.88,NULL,"brNDC"); + legend->SetLineColor(1); + legend->SetLineStyle(1); + legend->SetLineWidth(1); + legend->SetFillStyle(0); + legend->SetTextSize(0.05); + + // For Data where luminosity is stated + TLatex *t2 = new TLatex(0.13, 0.94, "CMS Preliminary"); + t2->SetNDC(kTRUE); + t2->SetTextSize(0.045); + t2->SetTextFont(62); + t2->SetLineWidth(5); + t2->Draw(""); + + float lumi = 40.7; + TString latexnamepre= "#sqrt{s} = 13 TeV, #scale[0.45]{#int}L = "; + TString latexnamemiddle; + latexnamemiddle.Form("%1.2f",lumi); + TString latexnamepost = " pb^{-1}"; + TString latexname = latexnamepre+latexnamemiddle+latexnamepost; + TLatex *t2a = new TLatex(0.72,0.95,latexname); + t2a->SetTextSize(0.045); + t2a->SetTextAlign(22); + t2a->SetNDC(kTRUE); + t2a->SetTextFont(42); + t2a->Draw(""); + + // for MC figures + TLatex *lar = new TLatex(0.13, 0.94, "CMS"); + lar->SetNDC(kTRUE); + lar->SetTextSize(0.050); + lar->SetTextFont(62); + lar->SetLineWidth(5); + lar->Draw(""); + + TLatex *lar2 = new TLatex(0.25, 0.94, "Simulation"); + lar2->SetNDC(kTRUE); + lar2->SetTextSize(0.050); + lar2->SetTextFont(52); + lar2->SetLineWidth(5); + lar2->Draw(""); + + TH1F* h1 = new TH1F("h1","",100,0,1000); + // when you want to change the x-axis range + float xmin=10; + float xmax=100; + h1->GetXaxis()->SetRangeUser(xmin,xmax); + // set title offset so that the title is not covered by labels + h1->GetXaxis()->SetTitleOffset(1.5); + h1->GetYaxis()->SetTitleOffset(1.5); + } +*/ diff --git a/interface/untuplizer.h b/interface/untuplizer.h new file mode 100644 index 0000000..133a3cd --- /dev/null +++ b/interface/untuplizer.h @@ -0,0 +1,730 @@ +/* Reader of ggNtuplizer's and other TTrees. + * + * Works both with old as well as with new ntuples. In particular, for tree + * branches containing vector<...> objects of elementary data types as well as + * two-dimensional vector > or vector > arrays, + * underlying arrays are returned (vector::front()) by TreeReader::Get*() call + * family. + * + * NOTE: vector cannot be handled in this way since booleans are packed + * as bits into bytes. + * + * NOTE: for branches containing variable length arrays, corresponding counter + * branches must be read off first. + * + * Usage of int, long, ... elementary data types for retrieving branch contents + * is not portable and should be avoided. This is due to the fact that + * sizeof(int) and sizeof(long) may differ on 32bit and 64bit architectures. + * Therefore, it is advised to use the following conventions for types of branch + * contents: + * 1. Int_t instead of int; + * 2. UInt_t instead of unsigned int; + * 3. Long64_t instead of long; + * 4. ULong64_t instead of unsigned long, etc. + * + * For instance, it is suggested to use + * Int_t nEle = reader.GetInt("nEle"); + * instead of + * int nEle = reader.GetInt("nEle"); + * + * It is generally safe, however, to use float and double in place of Float_t + * and Double_t. + * + * Some notes: + * 1. In-place modifications of arrays (within array boundaries) are allowed and + * encouraged. + * 2. The code loads only requested tree branches. Therefore, in general, an + * analysis will run much faster. + * 3. Branch contents are accessed directly by accessing leafs. Thus, it is + * assumed that all leafs have unique names (equal to corresponding names of + * branches) which is indeed the case when each branch contains only one + * leaf. + * 4. In the code, no protection against invalid code usage is provided. + * 5. Tree friends are not supported. + */ + +#ifndef UNTUPLIZER_H +#define UNTUPLIZER_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// prints a message and exits gracefully +#ifndef FATAL +#define FATAL(msg) do { fprintf(stderr, "FATAL: %s\n", msg); gSystem->Exit(1); } while (0) +#endif + +class TreeReader { + + public: + + // types of branch/leaf contents + // NOTE: some unsigned types may not be supported + enum ETypes { + kBool, // single 1-byte boolean (TLeafO) + kChar, // single 1-byte integer (TLeafB) + kShort, // single 2-byte integer (TLeafS) + kInt, // single 4-byte integer (TLeafI) + kFloat, // single 4-byte float (TLeafF) + kDouble, // single 8-byte float (TLeafD) + kLong64, // single 8-byte integer (TLeafL) + kArrBool, // array of 1-byte booleans (TLeafO only) + kArrChar, // array of 1-byte integers (TLeafB or vector) + kArrCharTLeaf, // array of 1-byte integers (TLeafB) + kArrCharVector, // array of 1-byte integers (vector); NOTE: char=signed char is assumed + kArrUCharVector, // array of 1-byte unsigned integers (vector) + kArrShort, // array of 2-byte integers (TLeafS or vector) + kArrShortTLeaf, // array of 2-byte integers (TLeafS) + kArrShortVector, // array of 2-byte integers (vector) + kArrUShortVector, // array of 2-byte unsigned integers (vector) + kArrInt, // array of 4-byte integers (TLeafI or vector) + kArrIntTLeaf, // array of 4-byte integers (TLeafI) + kArrIntVector, // array of 4-byte integers (vector) + kArrUIntVector, // array of 4-byte unsigned integers (vector) + kArrFloat, // array of 4-byte floats (TLeafF or vector) + kArrFloatTLeaf, // array of 4-byte floats (TLeafF) + kArrFloatVector, // array of 4-byte floats (vector) + kArrLong64, // array of 8-byte integers (TLeafL or vector) + kArrLong64TLeaf, // array of 8-byte integers (TLeafL) + kArrLong64Vector, // array of 8-byte integers (vector) + kArrULong64Vector,// array of 8-byte unsigned integers (vector) + kArrVectInt, // array of vector (vector >) + kArrVectFloat, // array of vector (vector >) + kArrString, // array of string + kArrStringVector, // array of string (vector) + kTObject, // general object inherited from TObject + kVoidPtr // all other data types + }; + + // TTree/TChain initializers + TreeReader(TTree* tree); + TreeReader(const char* path, const char* treename = "tree/treeMaker"); +// TreeReader(const char* path, const char* treename = "treeMaker"); + TreeReader(const char** paths, int npaths, const char* treename = "treeMaker"); + TreeReader(std::vector paths, const char* treename = "treeMaker"); + + virtual ~TreeReader(); + + // simple getters + TTree* GetTree() { return fTree; } + Long64_t GetEntriesFast() { return fTree->GetEntriesFast(); } + Bool_t HasMC() { return fkMC; } + + // useful to determine which type of variable to use for which branch + void Print(); + + // sets event number to retrieve next time TreeReader::Get*() called + void GetEntry(Long64_t entry); + + // returns either a pointer to first element of one- or multi-dimensional + // array or a pointer to an unprocessed object. + void* GetPtr(const char* branch_name, ETypes cktype = kVoidPtr, Int_t* nsize = NULL); + + // specializations of the call above; + // NOTE: use same methods with type casting for unsigned data types + Char_t* GetPtrChar (const char* bname) { return (Char_t*) GetPtr(bname, kArrChar); } + Short_t* GetPtrShort (const char* bname) { return (Short_t*) GetPtr(bname, kArrShort); } + Int_t* GetPtrInt (const char* bname) { return (Int_t*) GetPtr(bname, kArrInt); } + Float_t* GetPtrFloat (const char* bname) { return (Float_t*) GetPtr(bname, kArrFloat); } + Long64_t* GetPtrLong64 (const char* bname) { return (Long64_t*) GetPtr(bname, kArrLong64); } + TObject* GetPtrTObject(const char* bname) { return (TObject*) GetPtr(bname, kTObject); } + + // NOTE: this works only for old ntuples + Bool_t* GetPtrBool(const char* bname) { return (Bool_t*) GetPtr(bname, kArrBool); } + + // return branch values for elementary types + Bool_t GetBool (const char* bname) { return ((Bool_t*) GetPtr(bname, kBool)) [0]; } + Char_t GetChar (const char* bname) { return ((Char_t*) GetPtr(bname, kChar)) [0]; } + Short_t GetShort (const char* bname) { return ((Short_t*) GetPtr(bname, kShort)) [0]; } + Int_t GetInt (const char* bname) { return ((Int_t*) GetPtr(bname, kInt)) [0]; } + Float_t GetFloat (const char* bname) { return ((Float_t*) GetPtr(bname, kFloat)) [0]; } + Double_t GetDouble(const char* bname) { return ((Double_t*) GetPtr(bname, kDouble))[0]; } + Long64_t GetLong64(const char* bname) { return ((Long64_t*) GetPtr(bname, kLong64))[0]; } + + std::string* GetPtrString(const char* bname) { return (std::string*) GetPtr(bname, kArrString); } + Int_t GetPtrStringSize(){return fStringVectorSize;} + + + // vector > and vector > tree branches + std::vector* GetPtrVectorFloat(const char* bname, Int_t &nsize) { + return (std::vector*) GetPtr(bname, kArrVectFloat, &nsize); + } + std::vector* GetPtrVectorFloat(const char* bname) { + return (std::vector*) GetPtr(bname, kArrVectFloat); + } + std::vector* GetPtrVectorInt(const char* bname, Int_t &nsize) { + return (std::vector*) GetPtr(bname, kArrVectInt, &nsize); + } + std::vector* GetPtrVectorInt(const char* bname) { + return (std::vector*) GetPtr(bname, kArrVectInt); + } + + protected: + + void InitSingleTTree(const char* path, const char* treename); + void InitTChain(const char** paths, int npaths, const char* treename); + void FindLeaf(const char* bname); + + TFile* fFile; // file handle associated with fTree + TTree* fTree; // reference to TTree or TChain to read entries from + Int_t fTreeNum; // (for TChain) current tree number in list of TTrees + Bool_t fkMC; // if kTRUE, MC truth info is available + Long64_t fEntry; // fTree entry number to read with Get*() methods + Int_t fStringVectorSize; + + // caching + std::map fLeafIdx; // leaf name => index in arrays below + std::vector fLeafAddr; // cached leaf address + std::vector fLeafType; // cached type of leaf contents + std::vector fLeafValue; // cached payload address +}; + +//______________________________________________________________________________ +TreeReader::TreeReader(TTree* tree) : + fFile(0), + fTree(0), + fTreeNum(-1), + fkMC(kFALSE), + fEntry(-1), + fStringVectorSize(0) +{ + /* Associates an external TTree or TChain with this class. + * + * tree = reference to TTree or TChain. The caller is the owner of the + * object. + */ + + // verify sizes of elementary data types + if (sizeof(int) != 4 || sizeof(long) != 8 || sizeof(float) != 4) + FATAL("int/long/float of unsupported size"); + + fTree = tree; + + // find out availability of MC truth info (check existence of "nMC" branch) + fkMC = fTree->GetBranch("nMC") ? kTRUE : kFALSE; +} + +//______________________________________________________________________________ +TreeReader::TreeReader(const char* path, const char* treename) : + fFile(0), + fTree(0), + fTreeNum(-1), + fkMC(kFALSE), + fEntry(-1) +{ + /* Takes TTree from a root file. + * + * path = any root-supported path to a file with TTree. + */ + + InitSingleTTree(path, treename); +} + +//______________________________________________________________________________ +TreeReader::TreeReader(const char** paths, int npaths, const char* treename) : + fFile(0), + fTree(0), + fTreeNum(-1), + fkMC(kFALSE), + fEntry(-1) +{ + /* Makes TChain of trees from several root files. + * + * paths = array of any root-supported paths to root files with TTrees; + * npaths = number of elements in the array above. + */ + + InitTChain(paths, npaths, treename); +} + +//______________________________________________________________________________ +TreeReader::TreeReader(std::vector paths, const char* treename) : + fFile(0), + fTree(0), + fTreeNum(-1), + fkMC(kFALSE), + fEntry(-1) +{ + /* Makes TChain of trees from several root files. + * + * paths = array of any root-supported paths to root files with TTrees. + */ + + int npaths = (int) paths.size(); + const char* paths2[npaths]; + + // convert vector into array of null-terminated strings + for (int i = 0; i < npaths; i++) + paths2[i] = paths[i].c_str(); + + InitTChain(paths2, npaths, treename); +} + +//______________________________________________________________________________ +TreeReader::~TreeReader() +{ + /* Frees memory. + */ + + if (fTree) delete fTree; + if (fFile) delete fFile; +} + +//______________________________________________________________________________ +void TreeReader::Print() +{ + /* Prints to stdout branch names together with descriptions of their content + * as to be used in an analysis code. + * + * Useful to determine which type of variable to use for which branch. + * + * NOTE: for branches containing vector<...> objects of elementary data types + * or vector<...> objects of vector or vector data types, types + * of underlying arrays are shown (as to be returned by TreeReader::Get*() + * call family). + */ + + Printf("Branch names with content descriptions:"); + + // access leafs directly + TObjArray* leafs = fTree->GetListOfLeaves(); + + // leaf loop + for (int i = 0; i < leafs->GetEntriesFast(); i++) { + TLeaf* leaf = dynamic_cast(leafs->At(i)); + if (!leaf) + FATAL("leaf is NULL"); + + // vector<...> and other objects not inherited from TObject + if (leaf->IsA() == TLeafElement::Class()) { + // content descriptor + std::string descr(leaf->GetBranch()->GetClassName()); + + // known one- and two-dimensional vector<...> arrays + // NOTE: force showing data types of well-defined sizes + const char* types[] = { + "float", "Float_t", "char", "Char_t", "signed char", "unsigned char", "UChar_t", + "short", "short int", "signed short", "signed short int", "Short_t", + "unsigned short", "unsigned short int", "UShort_t", + "int", "signed", "signed int", "Int_t", "unsigned int", "unsigned", "UInt_t", + "long", "long int", "signed long", "signed long int", "Long64_t", "Long_t", + "unsigned long", "unsigned long int", "ULong64_t", "ULong_t", + "vector ", "vector ","vector "}; + const char* type_descr[] = { + "float", "float", "Char_t", "Char_t", "Char_t", "UChar_t", "UChar_t", + "Short_t", "Short_t", "Short_t", "Short_t", "Short_t", + "UShort_t", "UShort_t", "UShort_t", + "Int_t", "Int_t", "Int_t", "Int_t", "UInt_t", "UInt_t", "UInt_t", + "Long64_t", "Long64_t", "Long64_t", "Long64_t", "Long64_t", "Long64_t", + "ULong64_t", "ULong64_t", "ULong64_t", "Long64_t", + "vector", "vector", "vector"}; + + for (int c = 0; c < 34; c++) + if (descr.compare(Form("vector<%s>", types[c])) == 0) { + descr = Form("%s*", type_descr[c]); + break; + } + + Printf(" %-36s: %s", leaf->GetName(), descr.data()); + } + + // objects inherited from TObject + else if (leaf->IsA() == TLeafObject::Class()) { + // content descriptor + std::string descr(leaf->GetBranch()->GetClassName()); + + Printf(" %-36s: %s", leaf->GetName(), descr.data()); + } + + // elementary data types; fixed/variable length arrays of elementary types + else { + std::string descr; + + // NOTE: only necessary leaf types are listed + if (leaf->IsA() == TLeafF::Class()) + descr = "Float_t"; + else if (leaf->IsA() == TLeafI::Class()) + descr = "Int_t"; + else if (leaf->IsA() == TLeafL::Class()) + descr = "Long64_t"; + else if (leaf->IsA() == TLeafO::Class()) + descr = "Bool_t"; + else if (leaf->IsA() == TLeafD::Class()) + descr = "Double_t"; + else if (leaf->IsA() == TLeafB::Class()) + descr = "Char_t"; + else if (leaf->IsA() == TLeafS::Class()) + descr = "Short_t"; + else + FATAL(Form("unsupported leaf of class %s", leaf->ClassName())); + + // single elementary variable vs fixed/variable length array + if (!leaf->GetLeafCount() && leaf->GetLenStatic() == 1) + Printf(" %-36s: %s", leaf->GetName(), descr.data()); + else { + // remove leaf name from the description + std::string title(leaf->GetTitle()); + size_t ind = title.find('['); + if (ind >= title.npos) + FATAL("string::find('[') failed"); + + Printf(" %-36s: %s%s", leaf->GetName(), descr.data(), &title[ind]); + } + } + + } // leaf loop +} + +//______________________________________________________________________________ +void TreeReader::GetEntry(Long64_t entry) +{ + /* Sets event number to retrieve next time TreeReader::Get*() called. + */ + + if (fTree->IsA() != TChain::Class()) + fEntry = entry; + + // TChain requires special treatment + else { + fEntry = ((TChain*)fTree)->LoadTree(entry); + + // reset caches on switching to new TTree + if (fTreeNum != ((TChain*)fTree)->GetTreeNumber()) { + fLeafIdx.clear(); + fLeafType.clear(); + fLeafAddr.clear(); + fLeafValue.clear(); + + fTreeNum = ((TChain*)fTree)->GetTreeNumber(); + } + } + + // reset cache of addresses to leaf payloads + for (size_t i = 0; i < fLeafValue.size(); i++) + fLeafValue[i] = NULL; +} + +//______________________________________________________________________________ +void* TreeReader::GetPtr(const char* branch_name, ETypes cktype, Int_t* nsize) +{ + /* Returns either a pointer to first element of one- or multi-dimensional + * array or a pointer to an unprocessed object. + * + * In particular, for vector<...> arrays of elementary data types as well as + * for vector > and vector > arrays, + * vector::front() is returned. In this way, the code will work both with old + * ntuples as well as with the new ntuples (provided that underlying data + * types are the same in both cases). + * + * branch_name = name of branch to search for; + * cktype = if not kVoidPtr, an additional type verification is performed; + * nsize = if not NULL, filled with vector::size() for vector<...> branches. + */ + + // entry number in fLeafAddr, fLeafType and fLeafValue + int i; + + // find the entry number + std::map::const_iterator got = fLeafIdx.find(branch_name); + + if (got != fLeafIdx.end()) + i = got->second; + else { + // this code is executed once per branch and per root file + FindLeaf(branch_name); + i = fLeafType.size() - 1; + } + + // verify leaf type, if requested + if (cktype != kVoidPtr) { + if (cktype == kArrFloat) { + if (fLeafType[i] != kArrFloatTLeaf && fLeafType[i] != kArrFloatVector) + FATAL(Form("branch is not of type Float_t*: %s", branch_name)); + } else if (cktype == kArrInt) { + if (fLeafType[i] != kArrIntTLeaf && fLeafType[i] != kArrIntVector && + fLeafType[i] != kArrUIntVector) + FATAL(Form("branch is not of type (U)Int_t*: %s", branch_name)); + } else if (cktype == kArrChar) { + if (fLeafType[i] != kArrCharTLeaf && fLeafType[i] != kArrCharVector && + fLeafType[i] != kArrUCharVector) + FATAL(Form("branch is not of type (U)Char_t*: %s", branch_name)); + } else if (cktype == kArrShort) { + if (fLeafType[i] != kArrShortTLeaf && fLeafType[i] != kArrShortVector && + fLeafType[i] != kArrUShortVector) + FATAL(Form("branch is not of type (U)Short_t*: %s", branch_name)); + } else if (cktype == kArrLong64) { + if (fLeafType[i] != kArrLong64TLeaf && fLeafType[i] != kArrLong64Vector && + fLeafType[i] != kArrULong64Vector) + FATAL(Form("branch is not of type (U)Long64_t*: %s", branch_name)); + } else if (cktype == kTObject) { + if (fLeafType[i] != kTObject) + FATAL(Form("branch content is not inherited from TObject: %s", branch_name)); + } + else if (cktype == kArrString) { + if (fLeafType[i] != kArrStringVector) + FATAL(Form("branch is not of type string*: %s", branch_name)); + } + else + if (cktype != fLeafType[i]) + FATAL(Form("invalid branch type requested: %s", branch_name)); + } + + // load branch contents into memory, if necessary + if (!fLeafValue[i]) { + // NOTE: it is assumed here that the corresponding branch is not disabled + TBranch* br = fLeafAddr[i]->GetBranch(); + if (!br) + FATAL(Form("TLeaf::GetBranch() failed: %s", branch_name)); + if (br->GetEntry(fEntry) < 0) + FATAL(Form("TBranch::GetEntry() failed: %s", branch_name)); + + // pointer to actual leaf payload + void* ptr = fLeafAddr[i]->GetValuePointer(); + if (fLeafType[i] == kTObject) + ptr = *((void**)ptr); + + // cache address to payload + if (fLeafType[i] == kArrFloatVector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrIntVector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrUIntVector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrCharVector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrUCharVector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrShortVector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrUShortVector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrLong64Vector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrULong64Vector) + fLeafValue[i] = &((std::vector*)ptr)->front(); + else if (fLeafType[i] == kArrStringVector) + { + fLeafValue[i] = &((std::vector*)ptr)->front(); + fStringVectorSize = ((std::vector*)ptr)->size(); + } + else + fLeafValue[i] = ptr; + } + + // special case of vector > and vector > + if (fLeafType[i] == kArrVectFloat) { + if (nsize) + *nsize = (Int_t) ((std::vector >*)fLeafValue[i])->size(); + return &((std::vector >*)fLeafValue[i])->front(); + } + else if (fLeafType[i] == kArrVectInt) { + if (nsize) + *nsize = (Int_t) ((std::vector >*)fLeafValue[i])->size(); + return &((std::vector >*)fLeafValue[i])->front(); + } + + return fLeafValue[i]; +} + +//______________________________________________________________________________ +void TreeReader::InitSingleTTree(const char* path, const char* treename) +{ + /* Common code for the class constructors: open single TTree. + */ + + // verify sizes of elementary data types + if (sizeof(int) != 4 || sizeof(long) != 8 || sizeof(float) != 4) + FATAL("int/long/float of unsupported size"); + + // current working directory + TDirectory* wd = gDirectory; + + // open file with TTree + fFile = TFile::Open(path); //cout<<"hello"<IsZombie()) + FATAL("TFile::Open() failed"); + + // cd back into previous current working directory + if (wd) wd->cd(); + else gDirectory = 0; + + // get tree + fTree = dynamic_cast(fFile->Get(treename)); + if (!fTree) + FATAL(Form("TTree \"%s\" not found:", treename)); + + // be 100% sure: check explicitly object's class + if (((TObject*)fTree)->IsA() != TTree::Class()) + FATAL(Form("\"%s\" is not a TTree", treename)); + + // find out availability of MC truth info (check existence of "nMC" branch) + fkMC = fTree->GetBranch("nMC") ? kTRUE : kFALSE; +} + +//______________________________________________________________________________ +void TreeReader::InitTChain(const char** paths, int npaths, const char* treename) +{ + /* Common code for the class constructors: make TChain of several TTrees. + */ + + // verify sizes of elementary data types + if (sizeof(short) != 2 || sizeof(int) != 4 || sizeof(long) != 8 || sizeof(float) != 4) + FATAL("short/int/long/float of unsupported size"); + + // be smart + if (npaths == 1) { + InitSingleTTree(paths[0], treename); + return; + } + + fTree = new TChain(treename); + + // add root files with TTrees, reading the number of entries in each file + for (int i = 0; i < npaths; i++) + if (((TChain*)fTree)->AddFile(paths[i], 0) != 1) + FATAL("TChain::AddFile() failed"); + + // find out availability of MC truth info (check existence of "nMC" branch) + fkMC = fTree->GetBranch("nMC") ? kTRUE : kFALSE; +} + +//______________________________________________________________________________ +void TreeReader::FindLeaf(const char* bname) +{ + /* Finds leaf, determines type of its contents and fills the cache of leafs + * accordingly. + */ + + // find leaf + TLeaf* leaf = fTree->FindLeaf(bname); + if (!leaf) + FATAL(Form("leaf not found: %s", bname)); + + // actual data type of the leaf's payload + ETypes type; + + // vector<...> arrays and other objects not inherited from TObject + if (leaf->IsA() == TLeafElement::Class()) { + std::string descr(leaf->GetBranch()->GetClassName()); + + if (descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrFloatVector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrIntVector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrUIntVector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrCharVector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrUCharVector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrShortVector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrUShortVector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrLong64Vector; + else if (descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0 || + descr.compare("vector") == 0) + type = kArrULong64Vector; + else if (descr.compare("vector >") == 0 || + descr.compare("vector >") == 0) + type = kArrVectFloat; + else if (descr.compare("vector >") == 0 || + descr.compare("vector >") == 0 || + descr.compare("vector >") == 0 || + descr.compare("vector >") == 0) + type = kArrVectInt; + else if (descr.compare("vector") == 0) + type = kArrStringVector; + else + type = kVoidPtr; + } // TLeafElement + + // objects inherited from TObject + else if (leaf->IsA() == TLeafObject::Class()) + type = kTObject; + + // fixed/variable length arrays of elementary data types + else if (leaf->GetLeafCount() || leaf->GetLenStatic() > 1) { + if (leaf->IsA() == TLeafF::Class()) + type = kArrFloatTLeaf; + else if (leaf->IsA() == TLeafI::Class()) + type = kArrIntTLeaf; + else if (leaf->IsA() == TLeafB::Class()) + type = kArrCharTLeaf; + else if (leaf->IsA() == TLeafS::Class()) + type = kArrShortTLeaf; + else if (leaf->IsA() == TLeafL::Class()) + type = kArrLong64TLeaf; + else if (leaf->IsA() == TLeafO::Class()) + type = kArrBool; + else + FATAL(Form("branch contains an unknown data type: %s", bname)); + } + + // single variables of elementary data types + else if (!leaf->GetLeafCount() && leaf->GetLenStatic() == 1) { + if (leaf->IsA() == TLeafF::Class()) + type = kFloat; + else if (leaf->IsA() == TLeafI::Class()) + type = kInt; + else if (leaf->IsA() == TLeafB::Class()) + type = kChar; + else if (leaf->IsA() == TLeafS::Class()) + type = kShort; + else if (leaf->IsA() == TLeafD::Class()) + type = kDouble; + else if (leaf->IsA() == TLeafL::Class()) + type = kLong64; + else if (leaf->IsA() == TLeafO::Class()) + type = kBool; + else + FATAL(Form("branch contains an unknown data type: %s", bname)); + } + + // unknown TLeaf variant + else + FATAL(Form("branch contains an unknown data type: %s", bname)); + + // update the cache + fLeafIdx[bname] = (int)fLeafAddr.size(); + fLeafAddr.push_back(leaf); + fLeafType.push_back(type); + fLeafValue.push_back(NULL); +} + +#endif diff --git a/interface/xSec.txt b/interface/xSec.txt new file mode 100644 index 0000000..253de4f --- /dev/null +++ b/interface/xSec.txt @@ -0,0 +1,20 @@ +DYJetsToLL_M-50_HT-100to200_13TeV 139.4 +DYJetsToLL_M-50_HT-200to400_13TeV 42.75 +DYJetsToLL_M-50_HT-400to600_13TeV 5.497 +DYJetsToLL_M-50_HT-600toInf_13TeV 2.21 +TT_TuneCUETP8M1_13TeV 831.76 +WW_TuneCUETP8M1_13TeV 118.7 +WZ_TuneCUETP8M1_13TeV 47.13 +ZZ_TuneCUETP8M1_13TeV 16.523 +ZH_HToBB_ZToLL_M125_13TeV 0.05 +M-800_13TeV 0.0282665 +M-1000_13TeV 0.0153743 +M-1200_13TeV 0.00790857 +M-1400_13TeV 0.00421385 +M-1600_13TeV 0.00233319 +M-1800_13TeV 0.00133522 +M-2000_13TeV 0.000785119 +M-2500_13TeV 0.000227178 +M-3000_13TeV 0.000071426 +M-3500_13TeV 0.0000235715 +M-4000_13TeV 0.00000797489 diff --git a/plot_loop.C b/plot_loop.C new file mode 100644 index 0000000..96d8c19 --- /dev/null +++ b/plot_loop.C @@ -0,0 +1,54 @@ +#include +#include +#include "TString.h" +#include "TCanvas.h" +#include "VHPlotter.C" + +void plot_loop(){ + + const int N_histo = 22; + + TString Histo_Name[N_histo]; + + Histo_Name[0]="nVtx"; + Histo_Name[1]="leadElePt"; + Histo_Name[2]="leadEleEta"; + Histo_Name[3]="subleadElePt"; + Histo_Name[4]="subleadEleEta"; + Histo_Name[5]="Zmass"; + Histo_Name[6]="Zpt"; + Histo_Name[7]="Zeta"; + Histo_Name[8]="ZRapidity"; + Histo_Name[9]="FATjetPt"; + Histo_Name[10]="FATjetEta"; + Histo_Name[11]="FATjetCISVV2"; + Histo_Name[12]="FATjetSDmass"; + Histo_Name[13]="FATjetPRmass"; + Histo_Name[14]="FATjetPRmassCorr"; + Histo_Name[15]="FATjetTau1"; + Histo_Name[16]="FATjetTau2"; + Histo_Name[17]="FATjetTau2dvTau1"; + Histo_Name[18]="FATsubjetPt"; + Histo_Name[19]="FATsubjetEta"; + Histo_Name[20]="FATsubjetSDCSV"; + Histo_Name[21]="ZHmass"; + + TCanvas* c1 = new TCanvas("c","c",10,10,800,600); + TString save_path = "/afs/cern.ch/user/y/yuchang/www/ZH_2016/MC_Data_comparison/output_ele/ele_channel.pdf"; + + for(int i=0;iPrint(save_path + "("); + else if (i==N_histo-1) c1->Print(save_path + ")"); + else c1->Print(save_path ); + + c1->Clear(); + + } + +} + diff --git a/runPlot.sh b/runPlot.sh new file mode 100755 index 0000000..fafec0d --- /dev/null +++ b/runPlot.sh @@ -0,0 +1,9 @@ +#!/bin/sh + + +root -l -q -b plot_loop.C++ + +rm *.d *.so + + + diff --git a/xAna_ele.C b/xAna_ele.C new file mode 100644 index 0000000..7b3b239 --- /dev/null +++ b/xAna_ele.C @@ -0,0 +1,444 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "interface/untuplizer.h" +#include "interface/isPassZee.h" +#include "interface/readSample.h" +//#include "interface/ele_PU_71300/standalone_LumiReWeighting.cc" +#include "interface/ele_PU_69200/standalone_LumiReWeighting.cc" +#include "interface/dataFilter.h" + + +void xAna_ele(std::string inputFile, std::string outputFolder, std::string outputFile){ + + cout<<"finish compiling and start to run macro"<< endl; + // -------- open ntuples ------------ + + // get paths of every root files + std::vector infiles; + readSample(inputFile, infiles); + + cout<<"infiles.size():"<< infiles.size() << endl; + + // combine TTree + TChain* tree = new TChain("tree/treeMaker"); + + for(unsigned int i=0 ; iAdd( open_root.c_str() ); + } + + cout<<"tree->GetEntries(): "<< tree->GetEntries() << endl; + + // read the TTree + + TreeReader data( tree ); + cout<<"finish TreeReader data( tree );"<< endl; + + // ------- Declare the histogram ---------------- + + TH1D* h_totalEvents = new TH1D("h_totalEv", "totalEvents", 10, 0, 10); + + TH1D* h_pu_nTrueInt = new TH1D("h_pu_nTrueInt", "h_pu_nTrueInt", 60, 0, 60); + TH1D* h_nVtx = new TH1D("h_nVtx", "nVtx", 30, -0.5, 29.5); + + TH1D* h_leadElePt = new TH1D("h_leadElePt", "leadElePt", 32, 10, 810); + TH1D* h_leadEleEta = new TH1D("h_leadEleEta", "leadEleEta", 40, -4, 4); + TH1D* h_subleadElePt = new TH1D("h_subleadElePt", "subleadElePt", 25, 0, 500); + TH1D* h_subleadEleEta = new TH1D("h_subleadEleEta", "subleadEleEta", 40, -4, 4); + + TH1D* h_Zmass = new TH1D("h_Zmass", "Zmass", 30, 60, 120); + TH1D* h_Zpt = new TH1D("h_Zpt", "Zpt", 50, 0, 1000); + TH1D* h_Zeta = new TH1D("h_Zeta", "Zeta", 40, -4, 4); + TH1D* h_ZRapidity = new TH1D("h_ZRapidity", "ZRapidity", 40, -4, 4); + + TH1D* h_FATjetPt = new TH1D("h_FATjetPt", "FATjetPt", 50, 0, 1000); + TH1D* h_FATjetEta = new TH1D("h_FATjetEta", "FATjetEta", 20, -4, 4); + TH1D* h_FATjetCISVV2 = new TH1D("h_FATjetCISVV2", "FATjetCISVV2", 20, 0, 1.2); + TH1D* h_FATjetSDmass = new TH1D("h_FATjetSDmass", "FATjetSDmass", 20, 0, 200); + TH1D* h_FATjetPRmass = new TH1D("h_FATjetPRmass", "FATjetPRmass", 20, 0, 200); + TH1D* h_FATjetPRmassCorr = new TH1D("h_FATjetPRmassCorr", "FATjetPRmassCorr", 20, 0, 200); + TH1D* h_FATjetTau1 = new TH1D("h_FATjetTau1", "FATjetTau1", 20, 0, 1); + TH1D* h_FATjetTau2 = new TH1D("h_FATjetTau2", "FATjetTau2", 20, 0, 1); + TH1D* h_FATjetTau2dvTau1 = new TH1D("h_FATjetTau2dvTau1", "FATjetTau2dvTau1", 20, 0, 1); + + + TH1D* h_FATnSubjet = new TH1D("h_FATnSubjet", "FATnSubjet", 5, 0, 5); + TH1D* h_FATsubjetLeadingPt = new TH1D("h_FATsubjetLeadingPt", "FATsubjetLeadingPt", 20, 0, 800); + TH1D* h_FATsubjetLeadingEta = new TH1D("h_FATsubjetLeadingEta", "h_FATsubjetLeadingEta", 20, -4, 4); + TH1D* h_FATsubjetLeadingSDCSV = new TH1D("h_FATsubjetLeadingSDCSV", "h_FATsubjetLeadingSDCSV", 20, 0, 1.2); + TH1D* h_FATsubjetSubLeadingPt = new TH1D("h_FATsubjetSubLeadingPt", "FATsubjetSubLeadingPt", 20, 0, 800); + TH1D* h_FATsubjetSubLeadingEta = new TH1D("h_FATsubjetSubLeadingEta", "h_FATsubjetSubLeadingEta", 20, -4, 4); + TH1D* h_FATsubjetSubLeadingSDCSV = new TH1D("h_FATsubjetSubLeadingSDCSV","h_FATsubjetSubLeadingSDCSV",20, 0, 1.2); + + + TH1D* h_ZHmass = new TH1D("h_ZHmass", "ZHmass", 50, 0, 5000); + + h_pu_nTrueInt ->Sumw2(); + h_nVtx ->Sumw2(); + + h_leadElePt ->Sumw2(); + h_leadEleEta ->Sumw2(); + h_subleadElePt ->Sumw2(); + h_subleadEleEta ->Sumw2(); + + h_Zmass ->Sumw2(); + h_Zpt ->Sumw2(); + h_Zeta ->Sumw2(); + h_ZRapidity ->Sumw2(); + + h_FATjetPt ->Sumw2(); + h_FATjetEta ->Sumw2(); + h_FATjetCISVV2 ->Sumw2(); + h_FATjetSDmass ->Sumw2(); + h_FATjetPRmass ->Sumw2(); + h_FATjetPRmassCorr->Sumw2(); + h_FATjetTau1 ->Sumw2(); + h_FATjetTau2 ->Sumw2(); + h_FATjetTau2dvTau1->Sumw2(); + + h_FATnSubjet ->Sumw2(); + h_FATsubjetLeadingPt ->Sumw2(); + h_FATsubjetLeadingEta ->Sumw2(); + h_FATsubjetLeadingSDCSV ->Sumw2(); + h_FATsubjetSubLeadingPt ->Sumw2(); + h_FATsubjetSubLeadingEta ->Sumw2(); + h_FATsubjetSubLeadingSDCSV ->Sumw2(); + + h_ZHmass ->Sumw2(); + + int counter_0=0; + int counter_1=0; + int counter_2=0; + int counter_3=0; + + cout<<"finish define Histogram"<< endl; + + // ------ begin of event loop --------------- + +// int Number_to_print = 100000; + int Number_to_print = 5000; + +// int Max_number_to_read = 1000000; + int Max_number_to_read = 30000; +// int Max_number_to_read = 10; + + bool PU_weight_flag = true; +// PU_weight_flag = false; // turn off the PU weight + + if( PU_weight_flag ) cout<<"Will use Pile Up weight"< Max_number_to_read) break; + + data.GetEntry(ev); + + + // get Branch + + Bool_t isData = data.GetBool("isData"); + Float_t ntrue = data.GetFloat("pu_nTrueInt"); + Int_t nVtx = data.GetInt("nVtx"); + + const TClonesArray* eleP4 = (TClonesArray*) data.GetPtrTObject("eleP4"); + + Int_t FATnJet = data.GetInt("FATnJet"); + Int_t* FATnSubSDJet = data.GetPtrInt("FATnSubSDJet"); + Float_t* FATjetCISVV2 = data.GetPtrFloat("FATjetCISVV2"); + Float_t* FATjetSDmass = data.GetPtrFloat("FATjetSDmass"); + Float_t* FATjetPRmass = data.GetPtrFloat("FATjetPRmass"); + Float_t* FATjetPRmassCorr = data.GetPtrFloat("FATjetPRmassL2L3Corr"); + Float_t* FATjetTau1 = data.GetPtrFloat("FATjetTau1"); + Float_t* FATjetTau2 = data.GetPtrFloat("FATjetTau2"); + TClonesArray* FATjetP4 = (TClonesArray*) data.GetPtrTObject("FATjetP4"); + vector& FATjetPassIDLoose = *((vector*) data.GetPtr("FATjetPassIDLoose")); + vector* FATsubjetSDPx = data.GetPtrVectorFloat("FATsubjetSDPx", FATnJet); + vector* FATsubjetSDPy = data.GetPtrVectorFloat("FATsubjetSDPy", FATnJet); + vector* FATsubjetSDPz = data.GetPtrVectorFloat("FATsubjetSDPz", FATnJet); + vector* FATsubjetSDE = data.GetPtrVectorFloat("FATsubjetSDE", FATnJet); + vector* FATsubjetSDCSV = data.GetPtrVectorFloat("FATsubjetSDCSV", FATnJet); + + // get Pile-Up weight + + + double PU_weight_central =1;// weight=1 for data + + if(!isData && PU_weight_flag){// for MC + + standalone_LumiReWeighting LumiWeights_central(0); + PU_weight_central = LumiWeights_central.weight(ntrue); + + } + +// cout<<"PU_weight_central: "<< PU_weight_central << endl; + double eventWeight = PU_weight_central; + + h_totalEvents->Fill(1,eventWeight); + h_pu_nTrueInt ->Fill(ntrue , eventWeight); + + // data filter + + bool CSCT = FilterStatus(data, "Flag_CSCTightHaloFilter"); + bool eeBadSc = FilterStatus(data, "Flag_eeBadScFilter"); + bool Noise = FilterStatus(data, "Flag_HBHENoiseFilter"); + bool NoiseIso = FilterStatus(data, "Flag_HBHENoiseIsoFilter"); + + bool filter_pass = false; + if( isData && CSCT && eeBadSc && Noise && NoiseIso ) filter_pass = true;// Data + else if (!isData) filter_pass = true; // MC, doesn't apply data filter + + if( !filter_pass ) continue; + + + + // apply HLT + + bool eleTrigger1 = TriggerStatus(data, "HLT_Ele105_CaloIdVT_GsfTrkIdT_v"); + + bool HLT_pass = false; + if( isData && eleTrigger1 ) HLT_pass= true;// Data + else if (!isData) HLT_pass = true; // MC, doesn't apply HLT + + if( !HLT_pass ) continue; + + counter_0++; + + + + // nVtx>=1 + + if( nVtx < 1 ) continue; + + counter_1++; + + + + // select good electrons and Z boson + + + vector goodEleID; + if( !isPassZee(data, goodEleID) ) continue; + + counter_2++; + + TLorentzVector* thisEle = (TLorentzVector*)eleP4->At(goodEleID[0]); + TLorentzVector* thatEle = (TLorentzVector*)eleP4->At(goodEleID[1]); + + + + // select Fat jet for Higgs + + Int_t goodFATJetID = -1; + TLorentzVector* thisJet = NULL; + + for(Int_t ij = 0; ij < FATnJet; ++ij){ + + thisJet = (TLorentzVector*)FATjetP4->At(ij); + + if( thisJet->Pt() < 200 ) continue; + if( fabs(thisJet->Eta()) > 2.4 ) continue; + if( !FATjetPassIDLoose[ij] ) continue; + if( thisJet->DeltaR(*thisEle) < 0.8 || thisJet->DeltaR(*thatEle) < 0.8 ) continue; + if( !( FATjetPRmassCorr[ij] < 65 || FATjetPRmassCorr[ij] > 135) ) continue;// side band region + + goodFATJetID = ij; + break; + + } + + if( goodFATJetID < 0 ) continue; +// cout<<"goodFATJetID: "<< goodFATJetID << endl; + + counter_3++; + + + + // fill histogram in physical object level + + // nVtx + + h_nVtx ->Fill(nVtx,eventWeight); + + // Leading and SubLeading Electron + + if( thisEle->Pt() > thatEle->Pt() ){ + + h_leadElePt ->Fill(thisEle->Pt(),eventWeight); + h_leadEleEta ->Fill(thisEle->Eta(),eventWeight); + h_subleadElePt ->Fill(thatEle->Pt(),eventWeight); + h_subleadEleEta->Fill(thatEle->Eta(),eventWeight); + + }else{ + + h_leadElePt ->Fill(thatEle->Pt(),eventWeight); + h_leadEleEta ->Fill(thatEle->Eta(),eventWeight); + h_subleadElePt ->Fill(thisEle->Pt(),eventWeight); + h_subleadEleEta->Fill(thisEle->Eta(),eventWeight); + + } + + // Z boson + + TLorentzVector l4_Z = (*thisEle+*thatEle); + + h_Zmass ->Fill( l4_Z.M() ,eventWeight); + h_Zpt ->Fill( l4_Z.Pt() ,eventWeight); + h_Zeta ->Fill( l4_Z.Eta() ,eventWeight); + h_ZRapidity->Fill( l4_Z.Rapidity(),eventWeight); + + // Fat Jet for Higgs + + h_FATjetPt ->Fill(thisJet->Pt(),eventWeight); + h_FATjetEta ->Fill(thisJet->Eta(),eventWeight); + h_FATjetCISVV2 ->Fill(FATjetCISVV2[goodFATJetID],eventWeight); + h_FATjetSDmass ->Fill(FATjetSDmass[goodFATJetID],eventWeight); + h_FATjetPRmass ->Fill(FATjetPRmass[goodFATJetID],eventWeight); + h_FATjetPRmassCorr->Fill(FATjetPRmassCorr[goodFATJetID],eventWeight); + h_FATjetTau1 ->Fill(FATjetTau1[goodFATJetID],eventWeight); + h_FATjetTau2 ->Fill(FATjetTau2[goodFATJetID],eventWeight); + h_FATjetTau2dvTau1->Fill(FATjetTau2[goodFATJetID]/FATjetTau1[goodFATJetID],eventWeight); + + // SubJet + + int nSubJet = FATnSubSDJet[goodFATJetID]; + h_FATnSubjet->Fill( nSubJet , eventWeight ); + + std::vector subjetPt; subjetPt.clear(); + std::vector subjetEta; subjetEta.clear(); + + for(Int_t is = 0; is < FATnSubSDJet[goodFATJetID]; ++is){ + + + TLorentzVector l4_subjet(0,0,0,0); + + l4_subjet.SetPxPyPzE(FATsubjetSDPx[goodFATJetID][is], + FATsubjetSDPy[goodFATJetID][is], + FATsubjetSDPz[goodFATJetID][is], + FATsubjetSDE [goodFATJetID][is]); + + subjetPt.push_back( l4_subjet.Pt() ); + subjetEta.push_back( l4_subjet.Eta() ); + +// h_FATsubjetPt ->Fill(l4_subjet.Pt(),eventWeight); +// h_FATsubjetEta ->Fill(l4_subjet.Eta(),eventWeight); +// h_FATsubjetSDCSV->Fill(FATsubjetSDCSV[goodFATJetID][is],eventWeight); + + } + + if( nSubJet == 1 ){ + h_FATsubjetLeadingPt ->Fill( subjetPt[0] , eventWeight); + h_FATsubjetLeadingEta ->Fill( subjetEta[0] , eventWeight); + h_FATsubjetLeadingSDCSV ->Fill( FATsubjetSDCSV[goodFATJetID][0] , eventWeight); + } + + if( nSubJet > 1 ){ + + double leading_pt_temp =-99; double subleading_pt_temp =-99; + unsigned int leading_index =-1; unsigned int subleading_index =-1; + + for(unsigned int is=0; is leading_pt_temp ){leading_pt_temp = subjetPt[is] ; leading_index=is; } + } + + for(unsigned int is=0; is subleading_pt_temp ){subleading_pt_temp = subjetPt[is] ; subleading_index=is; } + } + + + h_FATsubjetLeadingPt ->Fill( subjetPt[leading_index] , eventWeight); + h_FATsubjetLeadingEta ->Fill( subjetEta[leading_index] , eventWeight); + h_FATsubjetLeadingSDCSV ->Fill( FATsubjetSDCSV[goodFATJetID][leading_index] , eventWeight); + + h_FATsubjetSubLeadingPt ->Fill( subjetPt[subleading_index] , eventWeight); + h_FATsubjetSubLeadingEta ->Fill( subjetEta[subleading_index] , eventWeight); + h_FATsubjetSubLeadingSDCSV ->Fill( FATsubjetSDCSV[goodFATJetID][subleading_index] , eventWeight); + + + } + + + + +// if(subjet_pt2>subjet_pt1)cout<<"subjet_pt1: "<Fill( l4_ZH.M() ,eventWeight); + + } // ------------------ end of event loop ------------------ + + cout<<"counter_0: "<< counter_0 << endl; + cout<<"counter_1: "<< counter_1 << endl; + cout<<"counter_2: "<< counter_2 << endl; + cout<<"counter_3: "<< counter_3 << endl; + + + fprintf(stderr, "Processed all events\n"); + + std::string root_name = Form("%s.root",outputFile.c_str() ) ; + TString save_path = outputFolder +"/"+ root_name; + + TFile* outFile = new TFile( save_path , "recreate"); + + h_pu_nTrueInt ->Write("pu_nTrueInt"); + h_nVtx ->Write("nVtx"); + + h_leadElePt ->Write("leadElePt"); + h_leadEleEta ->Write("leadEleEta"); + h_subleadElePt ->Write("subleadElePt"); + h_subleadEleEta ->Write("subleadEleEta"); + + h_Zmass ->Write("Zmass"); + h_Zpt ->Write("Zpt"); + h_Zeta ->Write("Zeta"); + h_ZRapidity ->Write("ZRapidity"); + + h_FATjetPt ->Write("FATjetPt"); + h_FATjetEta ->Write("FATjetEta"); + h_FATjetCISVV2 ->Write("FATjetCISVV2"); + h_FATjetSDmass ->Write("FATjetSDmass"); + h_FATjetPRmass ->Write("FATjetPRmass"); + h_FATjetPRmassCorr->Write("FATjetPRmassCorr"); + h_FATjetTau1 ->Write("FATjetTau1"); + h_FATjetTau2 ->Write("FATjetTau2"); + h_FATjetTau2dvTau1->Write("FATjetTau2dvTau1"); + + h_FATnSubjet ->Write("FATnSubjet"); + h_FATsubjetLeadingPt ->Write("FATsubjetLeadingPt"); + h_FATsubjetLeadingEta ->Write("FATsubjetLeadingEta"); + h_FATsubjetLeadingSDCSV ->Write("FATsubjetLeadingSDCSV"); + h_FATsubjetSubLeadingPt ->Write("FATsubjetSubLeadingPt"); + h_FATsubjetSubLeadingEta ->Write("FATsubjetSubLeadingEta"); + h_FATsubjetSubLeadingSDCSV ->Write("FATsubjetSubLeadingSDCSV"); + + h_ZHmass ->Write("ZHmass"); + + h_totalEvents ->Write("totalEvents"); + + outFile->Write(); + delete outFile; + + +} + + + + diff --git a/xAna_mu.C b/xAna_mu.C new file mode 100644 index 0000000..b29903e --- /dev/null +++ b/xAna_mu.C @@ -0,0 +1,368 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "interface/untuplizer.h" +#include "interface/isPassZmumu.h" +#include "interface/readSample.h" +#include "interface/standalone_LumiReWeighting.cc" +#include "interface/dataFilter.h" + + +void xAna_mu(std::string inputFile, std::string outputFolder, std::string outputFile){ + + cout<<"finish compiling and start to run macro"<< endl; + // -------- open ntuples ------------ + + // get paths of every root files + std::vector infiles; + readSample(inputFile, infiles); + + cout<<"infiles.size():"<< infiles.size() << endl; + + // combine TTree + TChain* tree = new TChain("tree/treeMaker"); + + for(unsigned int i=0 ; iAdd( open_root.c_str() ); + } + + cout<<"tree->GetEntries(): "<< tree->GetEntries() << endl; + + // read the TTree + + TreeReader data( tree ); + cout<<"finish TreeReader data( tree );"<< endl; + + // ------- Declare the histogram ---------------- + + TH1D* h_totalEvents = new TH1D("h_totalEv", "totalEvents", 10, 0, 10); + + TH1D* h_nVtx = new TH1D("h_nVtx", "nVtx", 30, -0.5, 29.5); + + TH1D* h_leadMuPt = new TH1D("h_leadMuPt", "leadMuPt", 50, 0, 1000); + TH1D* h_leadMuEta = new TH1D("h_leadMuEta", "leadMuEta", 40, -4, 4); + TH1D* h_subleadMuPt = new TH1D("h_subleadMuPt", "subleadMuPt", 25, 0, 500); + TH1D* h_subleadMuEta = new TH1D("h_subleadMuEta", "subleadMuEta", 40, -4, 4); + + TH1D* h_Zmass = new TH1D("h_Zmass", "Zmass", 30, 60, 120); + TH1D* h_Zpt = new TH1D("h_Zpt", "Zpt", 50, 0, 1000); + TH1D* h_Zeta = new TH1D("h_Zeta", "Zeta", 40, -4, 4); + TH1D* h_ZRapidity = new TH1D("h_ZRapidity", "ZRapidity", 40, -4, 4); + + TH1D* h_FATjetPt = new TH1D("h_FATjetPt", "FATjetPt", 20, 100, 1000); + TH1D* h_FATjetEta = new TH1D("h_FATjetEta", "FATjetEta", 20, -4, 4); + TH1D* h_FATjetCISVV2 = new TH1D("h_FATjetCISVV2", "FATjetCISVV2", 20, 0, 1.2); + TH1D* h_FATjetSDmass = new TH1D("h_FATjetSDmass", "FATjetSDmass", 20, 0, 200); + TH1D* h_FATjetPRmass = new TH1D("h_FATjetPRmass", "FATjetPRmass", 20, 0, 200); + TH1D* h_FATjetPRmassCorr = new TH1D("h_FATjetPRmassCorr", "FATjetPRmassCorr", 20, 0, 200); + TH1D* h_FATjetTau1 = new TH1D("h_FATjetTau1", "FATjetTau1", 20, 0, 1); + TH1D* h_FATjetTau2 = new TH1D("h_FATjetTau2", "FATjetTau2", 20, 0, 1); + TH1D* h_FATjetTau2dvTau1 = new TH1D("h_FATjetTau2dvTau1", "FATjetTau2dvTau1", 20, 0, 1); + + TH1D* h_FATsubjetPt = new TH1D("h_FATsubjetPt", "FATsubjetPt", 20, 0, 800); + TH1D* h_FATsubjetEta = new TH1D("h_FATsubjetEta", "FATsubjetEta", 20, -4, 4); + TH1D* h_FATsubjetSDCSV = new TH1D("h_FATsubjetSDCSV", "FATsubjetSDCSV", 20, 0, 1.2); + + TH1D* h_ZHmass = new TH1D("h_ZHmass", "ZHmass", 50, 0, 5000); + + h_nVtx ->Sumw2(); + + h_leadMuPt ->Sumw2(); + h_leadMuEta ->Sumw2(); + h_subleadMuPt ->Sumw2(); + h_subleadMuEta->Sumw2(); + + + h_Zmass ->Sumw2(); + h_Zpt ->Sumw2(); + h_Zeta ->Sumw2(); + h_ZRapidity ->Sumw2(); + + h_FATjetPt ->Sumw2(); + h_FATjetEta ->Sumw2(); + h_FATjetCISVV2 ->Sumw2(); + h_FATjetSDmass ->Sumw2(); + h_FATjetPRmass ->Sumw2(); + h_FATjetPRmassCorr->Sumw2(); + h_FATjetTau1 ->Sumw2(); + h_FATjetTau2 ->Sumw2(); + h_FATjetTau2dvTau1->Sumw2(); + h_FATsubjetPt ->Sumw2(); + h_FATsubjetEta ->Sumw2(); + h_FATsubjetSDCSV ->Sumw2(); + + h_ZHmass ->Sumw2(); + + int counter_0=0; + int counter_1=0; + int counter_2=0; + int counter_3=0; + + cout<<"finish define Histogram"<< endl; + + // ------ begin of event loop --------------- + +// int Number_to_print = 100000; + int Number_to_print = 5000; + +// int Max_number_to_read = 1000000; + int Max_number_to_read = 30000; +// int Max_number_to_read = 10; + + cout<<"staring loop events"<< endl; + + for( Long64_t ev = 0; ev < data.GetEntriesFast(); ev++ ){ + + if( ev % Number_to_print == 0 ) + fprintf(stderr, "Processing event %lli of %lli\n", ev + 1, data.GetEntriesFast()); + + if( ev > Max_number_to_read) break; + + data.GetEntry(ev); + + h_totalEvents->Fill(1,1); + + // get Branch + + Bool_t isData = data.GetBool("isData"); + + Int_t nVtx = data.GetInt("nVtx"); + + TClonesArray* muP4 = (TClonesArray*) data.GetPtrTObject("muP4"); + + Int_t FATnJet = data.GetInt("FATnJet"); + Int_t* FATnSubSDJet = data.GetPtrInt("FATnSubSDJet"); + Float_t* FATjetCISVV2 = data.GetPtrFloat("FATjetCISVV2"); + Float_t* FATjetSDmass = data.GetPtrFloat("FATjetSDmass"); + Float_t* FATjetPRmass = data.GetPtrFloat("FATjetPRmass"); + Float_t* FATjetPRmassCorr = data.GetPtrFloat("FATjetPRmassL2L3Corr"); + Float_t* FATjetTau1 = data.GetPtrFloat("FATjetTau1"); + Float_t* FATjetTau2 = data.GetPtrFloat("FATjetTau2"); + TClonesArray* FATjetP4 = (TClonesArray*) data.GetPtrTObject("FATjetP4"); + vector& FATjetPassIDLoose = *((vector*) data.GetPtr("FATjetPassIDLoose")); + vector* FATsubjetSDPx = data.GetPtrVectorFloat("FATsubjetSDPx", FATnJet); + vector* FATsubjetSDPy = data.GetPtrVectorFloat("FATsubjetSDPy", FATnJet); + vector* FATsubjetSDPz = data.GetPtrVectorFloat("FATsubjetSDPz", FATnJet); + vector* FATsubjetSDE = data.GetPtrVectorFloat("FATsubjetSDE", FATnJet); + vector* FATsubjetSDCSV = data.GetPtrVectorFloat("FATsubjetSDCSV", FATnJet); + + // get Pile-Up weight + + + double PU_weight_central =1;// weight=1 for data + + if(!isData ){// for MC + standalone_LumiReWeighting LumiWeights_central(0); + + Float_t ntrue= data.GetFloat("pu_nTrueInt"); + PU_weight_central = LumiWeights_central.weight(ntrue); + + } + + double eventWeight = PU_weight_central; + + // data filter + + bool CSCT = FilterStatus(data, "Flag_CSCTightHaloFilter"); + bool eeBadSc = FilterStatus(data, "Flag_eeBadScFilter"); + bool Noise = FilterStatus(data, "Flag_HBHENoiseFilter"); + bool NoiseIso = FilterStatus(data, "Flag_HBHENoiseIsoFilter"); + + bool filter_pass = false; + if( isData && CSCT && eeBadSc && Noise && NoiseIso ) filter_pass = true;// Data + else if (!isData) filter_pass = true; // MC, doesn't apply data filter + + if( !filter_pass ) continue; + + + + // apply HLT + + bool eleTrigger1 = TriggerStatus(data, "HLT_Ele105_CaloIdVT_GsfTrkIdT_v"); + + bool HLT_pass = false; + if( isData && eleTrigger1 ) HLT_pass= true;// Data + else if (!isData) HLT_pass = true; // MC, doesn't apply HLT + +// if( !HLT_pass ) continue; + + counter_0++; + + + + // nVtx>=1 + + if( nVtx < 1 ) continue; + + counter_1++; + + + + // select good muons and Z boson + + vector goodMuID; + if( !isPassZmumu(data, goodMuID) ) continue; + + counter_2++; + + + TLorentzVector* thisMu = (TLorentzVector*)muP4->At(goodMuID[0]); + TLorentzVector* thatMu = (TLorentzVector*)muP4->At(goodMuID[1]); + + + // select Fat jet for Higgs + + Int_t goodFATJetID = -1; + TLorentzVector* thisJet = NULL; + + for(Int_t ij = 0; ij < FATnJet; ++ij){ + + thisJet = (TLorentzVector*)FATjetP4->At(ij); + + if( thisJet->Pt() < 200 ) continue; + if( fabs(thisJet->Eta()) > 2.4 ) continue; + if( !FATjetPassIDLoose[ij] ) continue; + if( thisJet->DeltaR(*thisMu) < 0.8 || thisJet->DeltaR(*thatMu) < 0.8 ) continue; + if( !( FATjetPRmassCorr[ij] < 65 || FATjetPRmassCorr[ij] > 135) ) continue;// side band region + + goodFATJetID = ij; + break; + + } + + if( goodFATJetID < 0 ) continue; + + counter_3++; + + + + // fill histogram in physical object level + + // nVtx + + h_nVtx ->Fill(nVtx,eventWeight); + + // Leading and SubLeading Electron + + if( thisMu->Pt() > thatMu->Pt() ){ + + h_leadMuPt ->Fill(thisMu->Pt(),eventWeight); + h_leadMuEta ->Fill(thisMu->Eta(),eventWeight); + h_subleadMuPt ->Fill(thatMu->Pt(),eventWeight); + h_subleadMuEta->Fill(thatMu->Eta(),eventWeight); + + }else{ + + h_leadMuPt ->Fill(thatMu->Pt(),eventWeight); + h_leadMuEta ->Fill(thatMu->Eta(),eventWeight); + h_subleadMuPt ->Fill(thisMu->Pt(),eventWeight); + h_subleadMuEta->Fill(thisMu->Eta(),eventWeight); + + } + + // Z boson + + TLorentzVector l4_Z = (*thisMu+*thatMu); + + h_Zmass ->Fill( l4_Z.M() ,eventWeight); + h_Zpt ->Fill( l4_Z.Pt() ,eventWeight); + h_Zeta ->Fill( l4_Z.Eta() ,eventWeight); + h_ZRapidity->Fill( l4_Z.Rapidity(),eventWeight); + + // Fat Jet for Higgs + + h_FATjetPt ->Fill(thisJet->Pt(),eventWeight); + h_FATjetEta ->Fill(thisJet->Eta(),eventWeight); + h_FATjetCISVV2 ->Fill(FATjetCISVV2[goodFATJetID],eventWeight); + h_FATjetSDmass ->Fill(FATjetSDmass[goodFATJetID],eventWeight); + h_FATjetPRmass ->Fill(FATjetPRmass[goodFATJetID],eventWeight); + h_FATjetPRmassCorr->Fill(FATjetPRmassCorr[goodFATJetID],eventWeight); + h_FATjetTau1 ->Fill(FATjetTau1[goodFATJetID],eventWeight); + h_FATjetTau2 ->Fill(FATjetTau2[goodFATJetID],eventWeight); + h_FATjetTau2dvTau1->Fill(FATjetTau2[goodFATJetID]/FATjetTau1[goodFATJetID],eventWeight); + + // SubJet + + for(Int_t is = 0; is < FATnSubSDJet[goodFATJetID]; ++is){ + + TLorentzVector l4_subjet(0,0,0,0); + + l4_subjet.SetPxPyPzE(FATsubjetSDPx[goodFATJetID][is], + FATsubjetSDPy[goodFATJetID][is], + FATsubjetSDPz[goodFATJetID][is], + FATsubjetSDE [goodFATJetID][is]); + + h_FATsubjetPt ->Fill(l4_subjet.Pt(),eventWeight); + h_FATsubjetEta ->Fill(l4_subjet.Eta(),eventWeight); + h_FATsubjetSDCSV->Fill(FATsubjetSDCSV[goodFATJetID][is],eventWeight); + + } + + // ZH invariant mass + + TLorentzVector l4_ZH = l4_Z + ( *thisJet ); + + h_ZHmass ->Fill( l4_ZH.M() ,eventWeight); + + } // ------------------ end of event loop ------------------ + + cout<<"counter_0: "<< counter_0 << endl; + cout<<"counter_1: "<< counter_1 << endl; + cout<<"counter_2: "<< counter_2 << endl; + cout<<"counter_3: "<< counter_3 << endl; + + + fprintf(stderr, "Processed all events\n"); + + std::string root_name = Form("%s.root",outputFile.c_str() ) ; + TString save_path = outputFolder +"/"+ root_name; + + TFile* outFile = new TFile( save_path , "recreate"); + + h_nVtx ->Write("nVtx"); + + h_leadMuPt ->Write("leadMuPt"); + h_leadMuEta ->Write("leadMuEta"); + h_subleadMuPt ->Write("subleadMuPt"); + h_subleadMuEta->Write("subleadMuEta"); + + h_Zmass ->Write("Zmass"); + h_Zpt ->Write("Zpt"); + h_Zeta ->Write("Zeta"); + h_ZRapidity ->Write("ZRapidity"); + + h_FATjetPt ->Write("FATjetPt"); + h_FATjetEta ->Write("FATjetEta"); + h_FATjetCISVV2 ->Write("FATjetCISVV2"); + h_FATjetSDmass ->Write("FATjetSDmass"); + h_FATjetPRmass ->Write("FATjetPRmass"); + h_FATjetPRmassCorr->Write("FATjetPRmassCorr"); + h_FATjetTau1 ->Write("FATjetTau1"); + h_FATjetTau2 ->Write("FATjetTau2"); + h_FATjetTau2dvTau1->Write("FATjetTau2dvTau1"); + + h_FATsubjetPt ->Write("FATsubjetPt"); + h_FATsubjetEta ->Write("FATsubjetEta"); + h_FATsubjetSDCSV ->Write("FATsubjetSDCSV"); + + h_ZHmass ->Write("ZHmass"); + + h_totalEvents ->Write("totalEvents"); + + outFile->Write(); + delete outFile; + + +} + + + +