diff --git a/analysis/bin/BuildFile.xml b/analysis/bin/BuildFile.xml index ee7e64b3..d516f604 100644 --- a/analysis/bin/BuildFile.xml +++ b/analysis/bin/BuildFile.xml @@ -68,3 +68,16 @@ + + diff --git a/analysis/bin/doubleHiggsAnalyser.cc b/analysis/bin/doubleHiggsAnalyser.cc new file mode 100644 index 00000000..e920e2d9 --- /dev/null +++ b/analysis/bin/doubleHiggsAnalyser.cc @@ -0,0 +1,339 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "TMinuit.h" +#include "TError.h" + +#include "nano/analysis/interface/doubleHiggsAnalyser.h" + +#include "nano/external/interface/HiggsnessTopness.h" +#include "nano/oxbridgekinetics/src/Mt2/Basic_Mt2_332_Calculator.h" +#include "nano/oxbridgekinetics/src/Mt2/ChengHanBisect_Mt2_332_Calculator.h" + +using namespace std; + +void doubleHiggsAnalyser::MakeOutputBranch(TTree *tree) { + // MT2 variables + tree->Branch("basic_MT2_332_bbll",&basic_MT2_332_bbll,"basic_MT2_332_bbll/F"); + tree->Branch("basic_MT2_332_blbl",&basic_MT2_332_blbl,"basic_MT2_332_blbl/F"); + tree->Branch("basic_MT2_332_b",&basic_MT2_332_b,"basic_MT2_332_b/F"); + tree->Branch("basic_MT2_332_l",&basic_MT2_332_l,"basic_MT2_332_l/F"); + tree->Branch("mT",&mT,"mT/F"); + // lepton kinematic variables + tree->Branch("lep1", "TLorentzVector", &lepton1); + tree->Branch("lep2", "TLorentzVector", &lepton2); + tree->Branch("ll", "TLorentzVector", &leptonlepton); + tree->Branch("ll_M", &ll_M, "ll_M/F"); + tree->Branch("ll_deltaR", &ll_deltaR, "ll_deltaR/F"); + tree->Branch("ll_deltaPhi", &ll_deltaPhi, "ll_deltaPhi/F"); + // bottom kinematic variables + tree->Branch("bot1", "TLorentzVector", &bottom1); + tree->Branch("bot2", "TLorentzVector", &bottom1); + tree->Branch("bb", "TLorentzVector", &bottombottom); + tree->Branch("bb_deltaR", &bb_deltaR, "&bb_deltaR/F"); + tree->Branch("bb_deltaPhi", &bb_deltaPhi, "&bb_deltaPhi/F"); + // lepton + bottom + tree->Branch("bl11", "TLorentzVector", &bottomlepton11); + tree->Branch("bl12", "TLorentzVector", &bottomlepton12); + tree->Branch("bl21", "TLorentzVector", &bottomlepton21); + tree->Branch("bl22", "TLorentzVector", &bottomlepton22); + // lepton bottom kinematic variables + tree->Branch("bl_deltaR", "vector", &bl_deltaR); + tree->Branch("bl_min_deltaR", &bl_min_deltaR, "bl_min_deltaR/F"); + tree->Branch("bbll_deltaR", &bbll_deltaR, "bbll_delatR/F"); + tree->Branch("bbll_deltaPhi", &bbll_deltaPhi, "bbll_deltaPhi/F"); + // missing et + tree->Branch("MET","TLorentzVector",&missing); + tree->Branch("bbll","TLorentzVector",&bbll); + tree->Branch("topness",&topness,"topness/F"); + tree->Branch("higgsness",&higgsness,"higgsness/F"); + + tree->Branch("step", &step, "step/I"); + tree->Branch("tmva_bdtg_output", &tmva_bdtg_output, "tmva_bdtg_output/F"); +} + +void doubleHiggsAnalyser::SetTMVA(TString weight_file_path) { + bdtg_reader = new TMVA::Reader(); + bdtg_reader->AddVariable("ll_deltaR", &ll_deltaR); + bdtg_reader->AddVariable("ll.Pt()", &ll_Pt); + bdtg_reader->AddVariable("ll.M()", &ll_M); + bdtg_reader->AddVariable("bb_deltaR", &bb_deltaR); + bdtg_reader->AddVariable("bb.Pt()", &bb_Pt); + bdtg_reader->AddVariable("bb.M()", &bb_M); + bdtg_reader->AddVariable("bl_min_deltaR", &bl_min_deltaR); + bdtg_reader->AddVariable("bbll_deltaR", &bbll_deltaR); + bdtg_reader->AddVariable("bbll_deltaPhi", &bbll_deltaPhi); + bdtg_reader->AddVariable("mT", &mT); + bdtg_reader->AddVariable("basic_MT2_332_bbll", &basic_MT2_332_bbll); + bdtg_reader->BookMVA("BDTG", weight_file_path); + tmva_flag = true; +} + +void doubleHiggsAnalyser::SetOutput(TString output_file_name) { + out_file = TFile::Open(output_file_name,"RECREATE"); + out_tree = new TTree("events","events"); +} + +void doubleHiggsAnalyser::SetBranchAddress() { +} + +void doubleHiggsAnalyser::Initiate(TString output_file_name) { + // set output file + doubleHiggsAnalyser::SetOutput(output_file_name); + // make output branch + doubleHiggsAnalyser::MakeOutputBranch(out_tree); + // set branch address + doubleHiggsAnalyser::SetBranchAddress(); +} + +void doubleHiggsAnalyser::ResetVariables() { + //// MT2 variables + basic_MT2_332_bbll = std::numeric_limits::quiet_NaN(); + basic_MT2_332_blbl = std::numeric_limits::quiet_NaN(); + basic_MT2_332_b = std::numeric_limits::quiet_NaN(); + basic_MT2_332_l = std::numeric_limits::quiet_NaN(); + mT = std::numeric_limits::quiet_NaN(); + //// lepton variables + lepton1.Clear(); + lepton2.Clear(); + leptonlepton.Clear(); + ll_M = std::numeric_limits::quiet_NaN(); + ll_Pt = std::numeric_limits::quiet_NaN(); + ll_deltaR = std::numeric_limits::quiet_NaN(); + ll_deltaPhi = std::numeric_limits::quiet_NaN(); + leptons.clear(); + //// bottom variables + bottom1.Clear(); + bottom2.Clear(); + bottombottom.Clear(); + bb_M = std::numeric_limits::quiet_NaN(); + bb_Pt = std::numeric_limits::quiet_NaN(); + bb_deltaR = std::numeric_limits::quiet_NaN(); + bb_deltaPhi = std::numeric_limits::quiet_NaN(); + bottoms.clear(); + ////lepton and bottom variables + bottomlepton11.Clear(); + bottomlepton12.Clear(); + bottomlepton21.Clear(); + bottomlepton22.Clear(); + bbll.Clear(); + bl_deltaR.clear(); + bl_min_deltaR = std::numeric_limits::quiet_NaN(); + bbll_deltaR = std::numeric_limits::quiet_NaN(); + bbll_deltaPhi = std::numeric_limits::quiet_NaN(); + //// MET variables + missing.Clear(); + //// cut variables + step = 0; + //// higgsness and topness + higgsness = std::numeric_limits::quiet_NaN(); + topness = std::numeric_limits::quiet_NaN(); + //// tmva variables + tmva_bdtg_output = std::numeric_limits::quiet_NaN(); +} + +bool doubleHiggsAnalyser::Analysis() { + doubleHiggsAnalyser::ResetVariables(); + + // Missing ET // + if (MET_pt < 20) return false; + missing.SetPtEtaPhiM(MET_pt,0,MET_phi,0); + + // map, greater> leptons : map of : of leptons sorted by pt. + + // Muon Selection // + for (UInt_t i = 0; i < nMuon; i++) { + if (Muon_pt[i] < 10) continue; + if (Muon_eta[i] > 2.4) continue; + leptons.insert(make_pair(Muon_pt[i],make_pair(13*Muon_charge[i],i))); + } + // Electron Selection // + + if (leptons.size()<2) return false; + //// find two leptons with highest pt with opposite charge + lepton_iter = leptons.begin(); + auto l1_info = lepton_iter->second; + int pid1 = l1_info.first; + int index_l1 = l1_info.second; + lepton_iter++; + auto l2_info = lepton_iter->second; + int pid2 = l2_info.first; + int index_l2 = l2_info.second; + lepton_iter++; + while (pid1*pid2 > 0 && (lepton_iter!=leptons.end())) { + l2_info = lepton_iter->second; + pid2 = l2_info.first; + index_l2 = l2_info.second; + lepton_iter++; + } + if (pid1*pid2 > 0) return false; + lepton1.SetPtEtaPhiM(Muon_pt[index_l1], Muon_eta[index_l1], Muon_phi[index_l1], Muon_mass[index_l1]); + lepton2.SetPtEtaPhiM(Muon_pt[index_l2], Muon_eta[index_l2], Muon_phi[index_l2], Muon_mass[index_l2]); + leptonlepton = lepton1 + lepton2; + ll_deltaR = fabs(lepton1.DeltaR(lepton2)); + ll_deltaPhi = fabs(lepton1.DeltaPhi(lepton2)); + + //// Cuts on leptons + if (ll_deltaR < 0.07 || ll_deltaR > 3.3) {return true;} step++; + if (leptonlepton.M() < 5 || leptonlepton.M() > 100) {return true;} step++; + + // Bottom Selection // + for (UInt_t i = 0; i < nJet; i++) { + if (Jet_pt[i] < 30) continue; + if (Jet_eta[i] > 2.4) continue; + if (Jet_jetId[i] < 1) continue; + if (Jet_btagCSVV2[i] < 0.8484) continue; // medium cut + bottoms.insert(make_pair(Jet_pt[i],i)); + } + + if (bottoms.size()<2) return false; + //// find two leptons with highest pt with opposite charge + bottom_iter = bottoms.begin(); + int index_b1 = bottom_iter->second; + bottom_iter++; + int index_b2 = bottom_iter->second; + bottom1.SetPtEtaPhiM(Jet_pt[index_b1], Jet_eta[index_b1], Jet_phi[index_b1], Jet_mass[index_b1]); + bottom2.SetPtEtaPhiM(Jet_pt[index_b2], Jet_eta[index_b2], Jet_phi[index_b2], Jet_mass[index_b2]); + bottombottom = bottom1 + bottom2; + bb_deltaR = fabs(bottom1.DeltaR(bottom2)); + bb_deltaPhi = fabs(bottom1.DeltaPhi(bottom2)); + + //// Cuts on bottoms + if (bb_deltaR > 5) {return true;} step++; + if (bottombottom.M() < 22) {return true;} step++; + + // TLorentzVector combinations + //// 1 bottom and 1 lepton + bottomlepton11 = bottom1 + lepton1; + bottomlepton12 = bottom1 + lepton2; + bottomlepton21 = bottom2 + lepton1; + bottomlepton22 = bottom2 + lepton2; + //// 2 bottoms and 2 leptons and 1 missingET object combinations + bbll = bottombottom + leptonlepton; + bl_deltaR.push_back(lepton1.DeltaR(bottom1)); + bl_deltaR.push_back(lepton2.DeltaR(bottom1)); + bl_deltaR.push_back(lepton1.DeltaR(bottom2)); + bl_deltaR.push_back(lepton2.DeltaR(bottom2)); + bl_min_deltaR = *min_element(begin(bl_deltaR),end(bl_deltaR)); + bbll_deltaR = leptonlepton.DeltaR(bottombottom); + bbll_deltaPhi = leptonlepton.DeltaPhi(bottombottom); + + // MT2 variables + mT = sqrt(2*leptonlepton.Pt()*missing.E()*(1-cos(leptonlepton.Phi()-missing.Phi()))); + //// MT2(bb+ll) + Mt2::LorentzTransverseVector vis_A(Mt2::TwoVector(leptonlepton.Px(), leptonlepton.Py()), leptonlepton.M()); + Mt2::LorentzTransverseVector vis_B(Mt2::TwoVector(bottombottom.Px(),bottombottom.Py()), bottombottom.M()); + Mt2::TwoVector pT_Miss(missing.Px(), missing.Py()); + //// MT2(bl+bl) + Mt2::LorentzTransverseVector vis_A_blbl(Mt2::TwoVector(bottomlepton12.Px(), bottomlepton12.Py()), bottomlepton12.M()); + Mt2::LorentzTransverseVector vis_B_blbl(Mt2::TwoVector(bottomlepton21.Px(), bottomlepton21.Py()), bottomlepton21.M()); + //// MT2(b) + Mt2::LorentzTransverseVector vis_A_b(Mt2::TwoVector(bottom1.Px(), bottom1.Py()), bottom1.M()); + Mt2::LorentzTransverseVector vis_B_b(Mt2::TwoVector(bottom2.Px(), bottom2.Py()), bottom2.M()); + Mt2::TwoVector pT_Miss_b(leptonlepton.Px(), leptonlepton.Py()); + //// MT2(l) + Mt2::LorentzTransverseVector vis_A_l(Mt2::TwoVector(lepton1.Px(), lepton1.Py()), lepton1.M()); + Mt2::LorentzTransverseVector vis_B_l(Mt2::TwoVector(lepton2.Px(), lepton2.Py()), lepton2.M()); + Mt2::TwoVector pT_Miss_l(missing.Px(), missing.Py()); + + basic_MT2_332_bbll = basic_mt2_332Calculator.mt2_332(vis_A, vis_B, pT_Miss, missing.M()); + basic_MT2_332_blbl = basic_mt2_332Calculator.mt2_332(vis_A_blbl, vis_B_blbl, pT_Miss, missing.M()); + basic_MT2_332_b = basic_mt2_332Calculator.mt2_332(vis_A_b, vis_B_b, pT_Miss_b, leptonlepton.M()); + basic_MT2_332_l = basic_mt2_332Calculator.mt2_332(vis_A_l, vis_B_l, pT_Miss_l, missing.M()); + + // global TLorentzVector variables for higgsness, topness library + g_lepton1 = lepton1; + g_lepton2 = lepton2; + g_bottom1 = bottom1; + g_bottom2 = bottom2; + g_missing = missing; + // higgsness and topness + higgsness = GetHiggsness(); + topness = GetTopness(); + + // tmva variables + ll_Pt = leptonlepton.Pt(); + ll_M = leptonlepton.M(); + bb_Pt = bottombottom.Pt(); + bb_M = bottombottom.M(); + if(tmva_flag) tmva_bdtg_output = bdtg_reader->EvaluateMVA("BDTG"); + + return true; +} + +void doubleHiggsAnalyser::Loop() { + TTree *t = 0; + bool keep = false; + t = fChain; + int nevents = t->GetEntries(); + int proc = 0; int temp = 0; + int nevent = 0; + for (int iev = 0; iev < nevents; iev++) { + t->GetEntry(iev); + keep = doubleHiggsAnalyser::Analysis(); + if (keep) out_tree->Fill(); + nevent++; temp = nevent*100/nevents; + if (temp != proc) { + proc ++; + cout << "#######################" << endl; + cout << " proceeding : " << proc << " %" << endl; + cout << "#######################" << endl; + } + } + doubleHiggsAnalyser::Finalize(); // Write the tree and Close the file. +} + +void doubleHiggsAnalyser::Finalize() { + out_tree->Write(); + out_file->Close(); +} + +int main(Int_t argc, Char_t** argv) +{ + string env = getenv("CMSSW_BASE"); + string username = getenv("USER"); + + if (argc != 1) { + string dirName = "root://cms-xrdr.sdfarm.kr:1094//xrd/store/user/"+username+"/nanoAOD/"+std::string(argv[1])+"/"+std::string(argv[2]); + string temp = argv[2]; + Bool_t isMC = false; + Size_t found = temp.find("Run"); + if (found == string::npos) isMC = true; + for (Int_t i = 3; iGetObject("Events",tree); + temp = argv[i]; + found = temp.find_last_of('/'); + string outPutName = dirName+temp.substr(found); + doubleHiggsAnalyser ana(tree, isMC); + ana.Initiate(outPutName); + TString weight_file_path = "/home/sunyoung/nanoAOD/src/nano/analysis/bin/weight/DoubleHiggs_BDTG.weights.xml"; + ana.SetTMVA(weight_file_path); + ana.Loop(); // Loop through the events and do doubleHIggsAnalyser::Analysis() per event. + } + } else { + TChain *tree = new TChain("Events"); + //TString data_path = "/xrootd_UOS/store/group/nanoAOD/run2_2016v5/GluGluToHHTo2B2VTo2L2Nu_node_SM_13TeV-madgraph-v2/RunIISummer16MiniAODv2-PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/181108_163530/0000/*.root"; + TString data_path = "/xrootd_UOS/store/group/nanoAOD/run2_2016v5/TTJets_Dilept_TuneCUETP8M2T4_13TeV-amcatnloFXFX-pythia8/RunIISummer16MiniAODv2-PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/180607_120252/0000/*.root"; + tree->Add(data_path); + + tree->SetBranchStatus("*",true); + // for nanoAOD format analysis + doubleHiggsAnalyser ana(tree, true); + //ana.Initiate("HH_SM.root"); // usage : Initiate("outputfilename.root") + ana.Initiate("TT.root"); // usage : Initiate("outputfilename.root") + TString weight_file_path = "/cms/ldap_home/sunyoung/nanoAOD/src/nano/analysis/bin/weight/DoubleHiggs_BDTG.weights.xml"; + //ana.SetTMVA(weight_file_path); + //ana.Loop(); // Loop through the events and do doubleHIggsAnalyser::Analysis() per event. + } + return 0; +} diff --git a/analysis/interface/doubleHiggsAnalyser.h b/analysis/interface/doubleHiggsAnalyser.h new file mode 100644 index 00000000..a2226a14 --- /dev/null +++ b/analysis/interface/doubleHiggsAnalyser.h @@ -0,0 +1,98 @@ +#ifndef doubleHiggsAnalyser_H +#define doubleHiggsAnalyser_H +#include "nanoBase.h" +#include "nano/oxbridgekinetics/src/Mt2/Basic_Mt2_332_Calculator.h" +#include "nano/oxbridgekinetics/src/Mt2/ChengHanBisect_Mt2_332_Calculator.h" +#include "TMVA/Tools.h" +#include "TMVA/Reader.h" +#include "TMVA/DataLoader.h" +#include "TMVA/MethodCuts.h" + +class doubleHiggsAnalyser : public nanoBase { +private : + // Output Variables + TFile *out_file; + TTree *out_tree; + //// MT2 variables + Float_t basic_MT2_332_bbll = -99; + Float_t basic_MT2_332_blbl = -99; + Float_t basic_MT2_332_b = -99; + Float_t basic_MT2_332_l = -99; + Float_t mT = -99; + //// lepton variables + TLorentzVector lepton1; + TLorentzVector lepton2; + TLorentzVector leptonlepton; + Float_t ll_M = -99; + Float_t ll_Pt = -99; + Float_t ll_deltaR = -99; + Float_t ll_deltaPhi = -99; + std::map, std::greater> leptons; + std::map, std::greater>::iterator lepton_iter; + //// bottom variables + TLorentzVector bottom1; + TLorentzVector bottom2; + TLorentzVector bottombottom; + Float_t bb_M = -99; + Float_t bb_Pt = -99; + Float_t bb_deltaR = -99; + Float_t bb_deltaPhi = -99; + std::map> bottoms; + std::map>::iterator bottom_iter; + ////lepton and bottom variables + TLorentzVector bottomlepton11; + TLorentzVector bottomlepton12; + TLorentzVector bottomlepton21; + TLorentzVector bottomlepton22; + TLorentzVector bbll; + std::vector bl_deltaR; + Float_t bl_min_deltaR = -99; + Float_t bbll_deltaR = -99; + Float_t bbll_deltaPhi = -99; + //// MET variables + TLorentzVector missing; + //// cut variables + Int_t step = 0; + //// higgsness and topness + Float_t higgsness = -99; + Float_t topness = -99; + //// tmva variables + Float_t tmva_bdtg_output = -99; + bool tmva_flag = false; + + // Instances + //// TMVA reader + TMVA::Reader* bdtg_reader; + /// MT2 Calculators + Mt2::Basic_Mt2_332_Calculator basic_mt2_332Calculator; + Mt2::ChengHanBisect_Mt2_332_Calculator ch_bisect_mt2_332Calculator; + +public : + // before loop + void MakeOutputBranch(TTree *tree); + void SetOutput(TString output_file_name); + void SetBranchAddress(); + void SetTMVA(TString weight_file_path); + void Initiate(TString output_file_name); + // during loop + void ResetVariables(); + bool Analysis(); + // loop + virtual void Loop(); + // after loop + void Finalize(); + + doubleHiggsAnalyser(TTree *tree=0, Bool_t isMC = false); + ~doubleHiggsAnalyser(); +}; + +// Constructors ////// +doubleHiggsAnalyser::doubleHiggsAnalyser(TTree *tree, Bool_t isMC) : nanoBase(tree, isMC) { + std::string env = getenv("CMSSW_BASE"); + //m_rocCor = new RoccoR(env+"/src/nano/analysis/data/rcdata.2016.v3/"); +} + +// deconstructors +doubleHiggsAnalyser::~doubleHiggsAnalyser() { } + +#endif diff --git a/analysis/test/hh/script/draw_kinematic_variables.py b/analysis/test/hh/script/draw_kinematic_variables.py new file mode 100644 index 00000000..46197205 --- /dev/null +++ b/analysis/test/hh/script/draw_kinematic_variables.py @@ -0,0 +1,87 @@ +import ROOT, math +from ROOT import gStyle, gROOT + +hlist = { "bfc": {"tt":dict(), "hh":dict(), "hh_B6":dict(), "hh_B11":dict()}, "afc" : {"tt":dict(), "hh":dict(), "hh_B6":dict(), "hh_B11":dict()} } +sample_list = ["tt", "hh" ] +order_list = { 'll_pt': { "bfc":[0,1], "afc":[0,1] }, + 'bb_pt': { "bfc":[0,1], "afc":[0,1] }, + 'missing_et': { "bfc":[0,1], "afc":[0,1]}, + 'missing_et_phi': { "bfc":[1,0], "afc":[0,1] }, + 'bbll_mass': { "bfc":[0,1], "afc":[0,1] }, + "bb_mass": { "bfc":[1,0], "afc":[1,0] }, + "ll_mass": { "bfc":[1,0], "afc":[1,0] }, + "basic_MT2_332_bbll": { "bfc":[1,0], "afc":[1,0] }, + "basic_MT2_332_b": { "bfc":[1,0], "afc":[0,1] }, + "basic_MT2_332_l": { "bfc":[1,0], "afc":[1,0] }, + "basic_MT2_332_blbl": { "bfc":[1,0], "afc":[0,1] }, + "bb_deltaR": { "bfc":[1,0], "afc":[1,0] }, + "ll_deltaR": { "bfc":[1,0], "afc":[1,0] }, + "bbll_deltaR": { "bfc":[1,0], "afc":[1,0] }, + "bl_deltaR": { "bfc":[1,0], "afc":[1,0] }, + "bl_min_deltaR": { "bfc":[0,1], "afc":[1,0] }, + "tmva_bdtg_output": { "bfc":[0,1], "afc":[1,0] }, + "mT": { "bfc":[1,0], "afc":[1,0] }, + "hig_top_2d": { "bfc":[1,0], "afc":[1,0] } + } +color_list = {"tt":ROOT.kGreen+2, "hh":ROOT.kRed} + +f = ROOT.TFile("plots_1.root") +cvs = ROOT.TCanvas() +cvs.SetGrid() +cvs.SetLeftMargin(0.15) +cvs.SetBottomMargin(0.12) + +for key in ["bfc","afc"]: + for key2 in ["hh","tt"]: + for key3 in order_list.keys(): + hlist[key][key2][key3] = f.Get(key2+" "+key+" "+key3.replace("_"," ")) + +scale = {"hh":4/29999, "tt":26960510/14529280} + +for key in ["bfc","afc"]: + for key3 in order_list.keys(): + legend = ROOT.TLegend(0.70,0.55,0.90,0.90) + for i in [0,1]: + #legend = ROOT.TLegend(0.78,0.5,0.98,0.75) + h = hlist[key][sample_list[i]][key3] + #scale_factor = scale[sample_list[i]] + h.SetLineColor(color_list[sample_list[i]]) + legend.SetTextSize(0.06) + legend.AddEntry(h,sample_list[i]) + h.SetTitle("") + h.SetStats(False) + h.SetLabelSize(0.05,"X") + h.SetLabelSize(0.05,"Y") + h.SetLineWidth(3) + h.SetTickLength(0.03) + if 'deltaR' in key3: + h.GetXaxis().SetTitle(key3.replace("_"," ")) + h.GetYaxis().SetTitle("Normalized to unity") + elif 'top' in key3: + h.GetXaxis().SetTitle("log(H)") + h.GetYaxis().SetTitle("log(T)") + else: + h.GetXaxis().SetTitle(key3.replace("_"," ")+" [GeV]") + h.GetYaxis().SetTitle("Normalized to unity") + h.GetXaxis().SetTickSize(0.03) + h.GetXaxis().SetTitleSize(0.06) + h.GetXaxis().SetAxisColor(ROOT.kBlack) + cvs.RedrawAxis() + h.GetYaxis().SetAxisColor(ROOT.kBlack) + h.GetYaxis().SetTitleSize(0.05) + for ik, ikey2 in enumerate(order_list[key3][key]): + key2 = sample_list[ikey2] + h = hlist[key][key2][key3] + integ = h.Integral() + if integ != 0: + #h.Scale(scale_factor) + h.Scale(1/integ) + if ik==0: + h.Draw("hist") + else: + h.Draw("hist same") + legend.Draw("same") + cvs.SetTitle(key3.replace("_"," ")) + cvs.SaveAs(key3+"_"+key+".png") + +f.Close() diff --git a/analysis/test/hh/script/external_library_recipe.sh b/analysis/test/hh/script/external_library_recipe.sh new file mode 100644 index 00000000..cb2e28df --- /dev/null +++ b/analysis/test/hh/script/external_library_recipe.sh @@ -0,0 +1,29 @@ +cd $CMSSW_BASE/src/nano/ +wget http://www.hep.phy.cam.ac.uk/~lester/dtm662/mt2/Releases/oxbridgekinetics.tar.gz +tar -xzvf oxbridgekinetics.tar.gz +rm oxbridgekinetics.tar.gz +mv oxbridgekinetics-1.2/ oxbridgekinetics +cd oxbridgekinetics +touch BuildFile.xml +##################################################################### +#make a buildfile for integrating oxbridgekinetics library into nano +#specify flags in BuildFile to ignore unnecessary warnings and errors +cat > BuildFile.xml << EOL + + + + + + + + + +EOL +##################################################################### +mv Mt2/ src +cd src/Mt2 +sed -i s:"Mt2/:": * # modify the path of Mt2 header files. +cd ../../ +scram b +#set right path in the nano/analysis/bin/BuildFile.xml +sed -i '/erase/d' $CMSSW_BASE/src/nano/analysis/bin/BuildFile.xml diff --git a/analysis/test/hh/script/kinematic_variables.py b/analysis/test/hh/script/kinematic_variables.py new file mode 100644 index 00000000..839f7039 --- /dev/null +++ b/analysis/test/hh/script/kinematic_variables.py @@ -0,0 +1,95 @@ +import ROOT, math +hlist = { "bfc": {"tt":dict(), "hh":dict()}, "afc" : {"tt":dict(), "hh":dict()} } +limit_list = { 'll_pt': { "bfc":[60,0,400], "afc":[60,0,400] }, + 'bb_pt': { "bfc":[60,0,500], "afc":[60,0,500] }, + 'missing_et': { "bfc":[60,0,400], "afc":[60,0,400]}, + 'missing_et_phi': { "bfc":[60,-4,4], "afc":[60,-4,4] }, + 'bbll_mass': { "bfc":[100,0,700], "afc":[100,0,700] }, + "bb_mass": { "bfc":[100,0,400], "afc":[100,50,400] }, + "ll_mass": { "bfc":[100,0,400], "afc":[100,0,150] }, + "basic_MT2_332_bbll" : { "bfc":[80,0,500], "afc":[80,0,400] }, + "basic_MT2_332_blbl" : { "bfc":[80,0,500], "afc":[80,0,400] }, + "basic_MT2_332_b" : { "bfc":[80,0,500], "afc":[80,0,400] }, + "basic_MT2_332_l" : { "bfc":[80,0,500], "afc":[80,0,200] }, + "bb_deltaR" : { "bfc":[50,0,5], "afc":[50,0,5] }, + "ll_deltaR" : { "bfc":[50,0,5], "afc":[50,0,5] }, + "bbll_deltaR" : { "bfc":[50,0,5], "afc":[50,0,5] }, + "bl_deltaR" : { "bfc":[50,0,5], "afc":[50,0,5] }, + "bl_min_deltaR" : { "bfc":[50,0,5], "afc":[50,0,5] }, + "mT" : { "bfc":[75,0,300], "afc":[75,0,300] }, + "tmva_bdtg_output" : { "bfc":[50,-1.2,1.2], "afc":[50,-1.2,1.2] }, + "hig_top_2d" : { "bfc":[100,-10,10], "afc":[100,-10,10]} + } +def getHistogram(title,binning,lower_limit,upper_limit): + return ROOT.TH1D(title,title,binning,lower_limit,upper_limit) + +for key in ["bfc","afc"]: + for key2 in ["hh","tt"]: + for key3 in limit_list.keys(): + if key3 is "hig_top_2d": + hlist[key][key2][key3] = ROOT.TH2F(key2+" "+key+" "+key3.replace("_"," "),key2+" "+key+" "+key3.replace("_"," "),limit_list[key3][key][0],limit_list[key3][key][1],limit_list[key3][key][2],limit_list[key3][key][0],limit_list[key3][key][1],limit_list[key3][key][2]) + else: + hlist[key][key2][key3] = getHistogram(key2+" "+key+" "+key3.replace("_"," "),limit_list[key3][key][0],limit_list[key3][key][1],limit_list[key3][key][2]) + + + +tt_c = ROOT.TChain("events") +tt_c.Add("TT.root") +hh_c = ROOT.TChain("events") +hh_c.Add("HH_SM.root") +sample_list = {"tt":tt_c,"hh":hh_c} + +f = ROOT.TFile("plots_1.root","RECREATE") + +#for sample in ["tt", "hh", "hh_B6", "hh_B11"]: +for sample in ["tt", "hh"]: + c = sample_list[sample] + for ie, e in enumerate(c): + cut = "" + if e.step == 4 and e.tmva_bdtg_output > 0: + cut = "afc" + else: cut = "bfc" + hlist[cut][sample]['hig_top_2d'].Fill(e.higgsness, e.topness) + hlist[cut][sample]['ll_pt'].Fill(e.ll.Pt()) + hlist[cut][sample]['bb_pt'].Fill(e.bb.Pt()) + hlist[cut][sample]["missing_et_phi"].Fill(e.MET.Phi()) + hlist[cut][sample]["missing_et"].Fill(e.MET.E()) + hlist[cut][sample]["bbll_mass"].Fill(e.bbll.M()) + hlist[cut][sample]["bb_mass"].Fill(e.bb.M()) + hlist[cut][sample]["ll_mass"].Fill(e.ll.M()) + hlist[cut][sample]["basic_MT2_332_bbll"].Fill(e.basic_MT2_332_bbll) + hlist[cut][sample]["basic_MT2_332_blbl"].Fill(e.basic_MT2_332_blbl) + hlist[cut][sample]["basic_MT2_332_b"].Fill(e.basic_MT2_332_b) + hlist[cut][sample]["basic_MT2_332_l"].Fill(e.basic_MT2_332_l) + hlist[cut][sample]["bb_deltaR"].Fill(e.bb_deltaR) + hlist[cut][sample]["ll_deltaR"].Fill(e.ll_deltaR) + hlist[cut][sample]["bbll_deltaR"].Fill(e.bbll_deltaR) + hlist[cut][sample]["bl_min_deltaR"].Fill(e.bl_min_deltaR) + hlist[cut][sample]["tmva_bdtg_output"].Fill(e.tmva_bdtg_output) + hlist[cut][sample]["mT"].Fill(e.mT) + for bld in e.bl_deltaR: + hlist[cut][sample]["bl_deltaR"].Fill(bld) + +num_bg = hlist["afc"]["tt"]["hig_top_2d"].GetEntries() +num_bg = num_bg * 322627 / 14529280 +num_sg = hlist["afc"]["hh"]["hig_top_2d"].GetEntries() +num_sg = num_sg * 4. / 299999 +sensitivity = num_sg/math.sqrt(num_bg+num_sg) +print "number of background : " + str(num_bg) + "\n" +print "number of signal : " + str(num_sg) + "\n" +print "sensitivity : " + str(sensitivity) + "\n" + +for key in ["bfc","afc"]: + for key2 in ["tt", "hh"]: + #for key2 in ["tt", "hh", "hh_B6", "hh_B11"]: + for key3 in hlist[key][key2].keys(): + #integ = hlist[key][key2][key3].Integral() + #if integ == 0: + # print "%s %s %s has 0 integral" % (key, key2, key3) + #else: + # hlist[key][key2][key3].Scale(1/integ) + hlist[key][key2][key3].GetXaxis().SetTitle("log H") + hlist[key][key2][key3].GetYaxis().SetTitle("log T") + hlist[key][key2][key3].Write() + +f.Close() diff --git a/analysis/test/hh/tmva/DoubleHiggs.cc b/analysis/test/hh/tmva/DoubleHiggs.cc new file mode 100644 index 00000000..4fb80862 --- /dev/null +++ b/analysis/test/hh/tmva/DoubleHiggs.cc @@ -0,0 +1,206 @@ +#include +#include +#include +#include + +#include "TChain.h" +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TObjString.h" +#include "TSystem.h" +#include "TROOT.h" + +#include "TMVA/Factory.h" +#include "TMVA/DataLoader.h" +#include "TMVA/Tools.h" +#include "TMVA/TMVAGui.h" + +int DoubleHiggs( TString myMethodList = "" ) +{ + // This loads the library + TMVA::Tools::Instance(); + + // Default MVA methods to be trained + tested + std::map Use; + + Use["MLPBFGS"] = 1; // Recommended ANN with optional training method + Use["TMlpANN"] = 1; // ROOT's own ANN + // Boosted Decision Trees + Use["BDT"] = 1; // uses Adaptive Boost + Use["BDTG"] = 1; // uses Gradient Boost + // + // --------------------------------------------------------------- + + std::cout << std::endl; + std::cout << "==> Start DoubleHiggs" << std::endl; + + // Select methods (don't look at this code - not of interest) + if (myMethodList != "") { + for (std::map::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0; + + std::vector mlist = TMVA::gTools().SplitString( myMethodList, ',' ); + for (UInt_t i=0; i::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " "; + std::cout << std::endl; + return 1; + } + Use[regMethod] = 1; + } + } + + // -------------------------------------------------------------------------------------------------- + + // Here the preparation phase begins + + // Read training and test data + // (it is also possible to use ASCII format as input -> see TMVA Users Guide) + TFile *input(0); + TFile *input2(0); + + TString directory = "/home/sunyoung/nanoAOD/src/nano/analysis/test/hh/result/"; + TString sig_file = "HH_SM.root"; + TString bg_file = "TT.root"; + + input = TFile::Open( directory+sig_file ); + input2 = TFile::Open( directory+bg_file ); + + std::cout << "double Higgs using file : " << input->GetName() << std::endl; + + // Register the training and test trees + TTree *signalTree = (TTree*)input->Get("events"); + TTree *backgroundTree = (TTree*)input2->Get("events"); + + // Create a ROOT output file where TMVA will store ntuples, histograms, etc. + TString outfileName( "doubleHiggs_tmva" ); + TFile* outputFile = TFile::Open( outfileName, "RECREATE" ); + + TMVA::Factory *factory = new TMVA::Factory( "DoubleHiggs", outputFile, + "!V:!Silent:Color:DrawProgressBar:Transformations=I:AnalysisType=Classification" ); +// "!V:!Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D:AnalysisType=Classification" ); + + TMVA::DataLoader *dataloader=new TMVA::DataLoader("dataset_"+sig_file); + TMVA::DataLoader *dataloader2=new TMVA::DataLoader("dataset_"+bg_file); + + // (TMVA::gConfig().GetVariablePlotting()).fTimesRMS = 8.0; + // (TMVA::gConfig().GetIONames()).fWeightFileDir = "myWeightDirectory"; + (TMVA::gConfig().GetVariablePlotting()).fMaxNumOfAllowedVariablesForScatterPlots = 99999; + (TMVA::gConfig().GetVariablePlotting()).fNbinsXOfROCCurve = 1000; + + // [all types of expressions that can also be parsed by TTree::Draw( "expression" )] + + dataloader->AddVariable( "ll_deltaR", 'F' ); + dataloader->AddVariable( "ll.Pt()", 'F' ); + dataloader->AddVariable( "ll.M()", 'F' ); + dataloader->AddVariable( "bb_deltaR", 'F' ); + dataloader->AddVariable( "bb.Pt()", 'F' ); + dataloader->AddVariable( "bb.M()", 'F' ); + dataloader->AddVariable( "bl_min_deltaR", 'F' ); + dataloader->AddVariable( "bbll_deltaR", 'F' ); + dataloader->AddVariable( "bbll_deltaPhi", 'F' ); + dataloader->AddVariable( "mT", 'F'); + dataloader->AddVariable( "basic_MT2_332_bbll", 'F' ); + //dataloader->AddSpectator( "mass_had", 'F' ); + //dataloader2->AddVariable( "log_d_had := log(d_had)", 'F' ); + dataloader2->AddVariable( "ll_deltaR", 'F' ); + dataloader2->AddVariable( "ll.Pt()", 'F' ); + dataloader2->AddVariable( "ll.M()", 'F' ); + dataloader2->AddVariable( "bb_deltaR", 'F' ); + dataloader2->AddVariable( "bb.Pt()", 'F' ); + dataloader2->AddVariable( "bb.M()", 'F' ); + dataloader2->AddVariable( "bl_min_deltaR", 'F' ); + dataloader2->AddVariable( "bbll_deltaR", 'F' ); + dataloader2->AddVariable( "bbll_deltaPhi", 'F' ); + dataloader2->AddVariable( "mT", 'F'); + dataloader2->AddVariable( "basic_MT2_332_bbll", 'F' ); + //dataloader2->AddSpectator( "mass_had", 'F' ); + + // global event weights per tree (see below for setting event-wise weights) + Double_t signalWeight = 1.0; + Double_t backgroundWeight = 1.0; + + // You can add an arbitrary number of signal or background trees + dataloader->AddSignalTree ( signalTree, signalWeight ); + dataloader->AddBackgroundTree( backgroundTree, backgroundWeight ); + + //dataloader->SetBackgroundWeightExpression( "weight" ); + + //TCut mycuts = "abs(isFrom_had) == 3 && fabs(mass_had - 0.5) < 0.05"; + TCut mycuts = "step==4"; + TCut mycutb = "step==4"; + //TCut mycutb = "abs(isFrom_had) != 3 && abs(isFrom_had) != 99 && fabs(mass_had - 0.5) < 0.05"; + + //TCut mycuts2 = "abs(isFrom_had) == 3"; + // TCut mycutb2 = "abs(isFrom_had) == 5"; + + + dataloader->PrepareTrainingAndTestTree( mycuts, mycutb, +// "nTrain_Signal=150:nTrain_Background=700:SplitMode=Random:NormMode=NumEvents:!V" ); + "nTrain_Signal=5000:nTrain_Background=30000:SplitMode=Random:NormMode=NumEvents:!V" ); + + //if (Use["MLPBFGS"]) + // factory->BookMethod( dataloader, TMVA::Types::kMLP, "MLPBFGS", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:!UseRegulator" ); + //if (Use["TMlpANN"]) + // factory->BookMethod( dataloader, TMVA::Types::kTMlpANN, "TMlpANN", "!H:!V:NCycles=200:HiddenLayers=N+1,N:LearningMethod=BFGS:ValidationFraction=0.3" ); + if (Use["BDTG"]) // Gradient Boost + factory->BookMethod( dataloader, TMVA::Types::kBDT, "BDTG", + "!H:!V:NTrees=700:MinNodeSize=2.5%:BoostType=Grad:Shrinkage=0.10:UseBaggedBoost:BaggedSampleFraction=0.5:nCuts=20:MaxDepth=2" ); + if (Use["BDT"]) // Adaptive Boost + factory->BookMethod( dataloader, TMVA::Types::kBDT, "BDT", + "!H:!V:NTrees=700:MinNodeSize=2.5%:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:UseBaggedBoost:BaggedSampleFraction=0.5:SeparationType=GiniIndex:nCuts=20" ); + + + // For an example of the category classifier usage, see: DoubleHiggsCategory + // + // -------------------------------------------------------------------------------------------------- + // Now you can optimize the setting (configuration) of the MVAs using the set of training events + // STILL EXPERIMENTAL and only implemented for BDT's ! + // + // factory->OptimizeAllMethods("SigEffAt001","Scan"); + // factory->OptimizeAllMethods("ROCIntegral","FitGA"); + // + // -------------------------------------------------------------------------------------------------- + + // Now you can tell the factory to train, test, and evaluate the MVAs + // + // Train MVAs using the set of training events + factory->TrainAllMethods(); + + // Evaluate all MVAs using the set of test events + factory->TestAllMethods(); + + // Evaluate and compare performance of all configured MVAs + factory->EvaluateAllMethods(); + + // -------------------------------------------------------------- + + // Save the output + outputFile->Close(); + + std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl; + std::cout << "==> DoubleHiggs is done!" << std::endl; + + delete factory; + delete dataloader; + // Launch the GUI for the root macros + if (!gROOT->IsBatch()) TMVA::TMVAGui( outfileName ); + + return 0; +} + +int main( int argc, char** argv ) +{ + // Select methods (don't look at this code - not of interest) + TString methodList; + for (int i=1; i +#include "TMinuit.h" +#include "TError.h" + +using namespace std; + +// global variables for topness higgsness +const double mtop = 173.; +const double mh = 125; +const double mw = 80.419; +const double mhw = 30; + +const double sigmahlep = 2.; +const double sigmaWon = 5.; +const double sigmaWoff = 5.; +const double sigmat = 5.; + +const int n_unknown_par = 4; + +const double minuit_pmin = -5000; +const double minuit_pmax = 5000; + +TLorentzVector g_missing, g_lepton1, g_lepton2, g_bottom1, g_bottom2; + +TMinuit *ptMinuit; +TMinuit *ptMinuitT1; +TMinuit *ptMinuitT2; +Double_t minuit_arglist[10]; +Int_t ierflag = 0; + +double GetMass(double Px, double Py, double Pz, double E) { + double f = E*E - Px*Px - Py*Py - Pz*Pz; + double fout = 0; + if ( f < 0. ) { + fout = 0; + } else { + fout = sqrt(f); + } + return fout; +} + +void fcnT1( int& npar, double* deriv, double& f, double par[], int flag) { + double chisq = 0; + double Nu1Px = par[0]; + double Nu1Py = par[1]; + double Nu1Pz = par[2]; + double Nu1E = sqrt( Nu1Px*Nu1Px + Nu1Py*Nu1Py + Nu1Pz*Nu1Pz ); + + double Nu2Px = g_missing.Px() - Nu1Px; + double Nu2Py = g_missing.Py() - Nu1Py; + double Nu2Pz = par[3]; + double Nu2E = sqrt( Nu2Px*Nu2Px + Nu2Py*Nu2Py + Nu2Pz*Nu2Pz ); + + double Lep1E = g_lepton1.E(); + double Lep1Px = g_lepton1.Px(); + double Lep1Py = g_lepton1.Py(); + double Lep1Pz = g_lepton1.Pz(); + + double Lep2E = g_lepton2.E(); + double Lep2Px = g_lepton2.Px(); + double Lep2Py = g_lepton2.Py(); + double Lep2Pz = g_lepton2.Pz(); + + double Jet1E = g_bottom1.E(); + double Jet1Px = g_bottom1.Px(); + double Jet1Py = g_bottom1.Py(); + double Jet1Pz = g_bottom1.Pz(); + + double Jet2E = g_bottom2.E(); + double Jet2Px = g_bottom2.Px(); + double Jet2Py = g_bottom2.Py(); + double Jet2Pz = g_bottom2.Pz(); + + double jln111Px = Jet1Px + Lep1Px + Nu1Px; + double jln111Py = Jet1Py + Lep1Py + Nu1Py; + double jln111Pz = Jet1Pz + Lep1Pz + Nu1Pz; + double jln111E = Jet1E + Lep1E + Nu1E; + + double mjln111 = GetMass( jln111Px, jln111Py, jln111Pz, jln111E ); + + double deltaJLN111 = (mjln111*mjln111 - mtop*mtop) / (sigmat*sigmat); + + double jln222Px = Jet2Px + Lep2Px + Nu2Px; + double jln222Py = Jet2Py + Lep2Py + Nu2Py; + double jln222Pz = Jet2Pz + Lep2Pz + Nu2Pz; + double jln222E = Jet2E + Lep2E + Nu2E; + + double mjln222 = GetMass( jln222Px, jln222Py, jln222Pz, jln222E ); + + double deltaJLN222 = (mjln222*mjln222 - mtop*mtop) / (sigmat*sigmat); + + double ln11Px = Lep1Px + Nu1Px; + double ln11Py = Lep1Py + Nu1Py; + double ln11Pz = Lep1Pz + Nu1Pz; + double ln11E = Lep1E + Nu1E; + + double mln11 = GetMass( ln11Px, ln11Py, ln11Pz, ln11E ); + double deltaLN11on = (mln11*mln11 - mw*mw) / (sigmaWon*sigmaWon); + + double ln22Px = Lep2Px + Nu2Px; + double ln22Py = Lep2Py + Nu2Py; + double ln22Pz = Lep2Pz + Nu2Pz; + double ln22E = Lep2E + Nu2E; + + double mln22 = GetMass( ln22Px, ln22Py, ln22Pz, ln22E ); + + double deltaLN22on = (mln22*mln22 - mw*mw) / (sigmaWon*sigmaWon); + + chisq = deltaJLN111*deltaJLN111 + deltaJLN222*deltaJLN222 + deltaLN11on*deltaLN11on + deltaLN22on*deltaLN22on; + + f = chisq; +} + +void fcnT2( int& npar, double* deriv, double& f, double par[], int flag) { + double chisq = 0; + double Nu1Px = par[0]; + double Nu1Py = par[1]; + double Nu1Pz = par[2]; + double Nu1E = sqrt( Nu1Px*Nu1Px + Nu1Py*Nu1Py + Nu1Pz*Nu1Pz ); + + double Nu2Px = g_missing.Px() - Nu1Px; + double Nu2Py = g_missing.Py() - Nu1Py; + double Nu2Pz = par[3]; + double Nu2E = sqrt( Nu2Px*Nu2Px + Nu2Py*Nu2Py + Nu2Pz*Nu2Pz ); + + double Lep1E = g_lepton1.E(); + double Lep1Px = g_lepton1.Px(); + double Lep1Py = g_lepton1.Py(); + double Lep1Pz = g_lepton1.Pz(); + + double Lep2E = g_lepton2.E(); + double Lep2Px = g_lepton2.Px(); + double Lep2Py = g_lepton2.Py(); + double Lep2Pz = g_lepton2.Pz(); + + double Jet1E = g_bottom1.E(); + double Jet1Px = g_bottom1.Px(); + double Jet1Py = g_bottom1.Py(); + double Jet1Pz = g_bottom1.Pz(); + + double Jet2E = g_bottom2.E(); + double Jet2Px = g_bottom2.Px(); + double Jet2Py = g_bottom2.Py(); + double Jet2Pz = g_bottom2.Pz(); + + double jln121Px = Jet1Px + Lep2Px + Nu1Px; + double jln121Py = Jet1Py + Lep2Py + Nu1Py; + double jln121Pz = Jet1Pz + Lep2Pz + Nu1Pz; + double jln121E = Jet1E + Lep2E + Nu1E; + + double mjln121 = GetMass( jln121Px, jln121Py, jln121Pz, jln121E ); + + double deltaJLN121 = (mjln121*mjln121 - mtop*mtop) / (sigmat*sigmat); + + double jln212Px = Jet2Px + Lep1Px + Nu2Px; + double jln212Py = Jet2Py + Lep1Py + Nu2Py; + double jln212Pz = Jet2Pz + Lep1Pz + Nu2Pz; + double jln212E = Jet2E + Lep1E + Nu2E; + + double mjln212 = GetMass( jln212Px, jln212Py, jln212Pz, jln212E ); + + double deltaJLN212 = (mjln212*mjln212 - mtop*mtop) / (sigmat*sigmat); + + double ln21Px = Lep2Px + Nu1Px; + double ln21Py = Lep2Py + Nu1Py; + double ln21Pz = Lep2Pz + Nu1Pz; + double ln21E = Lep2E + Nu1E; + + double mln21 = GetMass( ln21Px, ln21Py, ln21Pz, ln21E ); + double deltaLN21on = (mln21*mln21 - mw*mw) / (sigmaWon*sigmaWon); + + double ln12Px = Lep1Px + Nu2Px; + double ln12Py = Lep1Py + Nu2Py; + double ln12Pz = Lep1Pz + Nu2Pz; + double ln12E = Lep1E + Nu2E; + + double mln12 = GetMass( ln12Px, ln12Py, ln12Pz, ln12E ); + + double deltaLN12on = (mln12*mln12 - mw*mw) / (sigmaWon*sigmaWon); + + chisq = deltaJLN121*deltaJLN121 + deltaJLN212*deltaJLN212 + deltaLN21on*deltaLN21on + deltaLN12on*deltaLN12on; + + f = chisq; +} + +void fcnH( int& npar, double* deriv, double& f, double par[], int flag) { + double chisq = 0; + double Nu1Px = par[0]; + double Nu1Py = par[1]; + double Nu1Pz = par[2]; + double Nu1E = sqrt( Nu1Px*Nu1Px + Nu1Py*Nu1Py + Nu1Pz*Nu1Pz ); + + double Nu2Px = g_missing.Px() - Nu1Px; + double Nu2Py = g_missing.Py() - Nu1Py; + double Nu2Pz = par[3]; + double Nu2E = sqrt( Nu2Px*Nu2Px + Nu2Py*Nu2Py + Nu2Pz*Nu2Pz ); + + double Lep1E = g_lepton1.E(); + double Lep1Px = g_lepton1.Px(); + double Lep1Py = g_lepton1.Py(); + double Lep1Pz = g_lepton1.Pz(); + + double Lep2E = g_lepton2.E(); + double Lep2Px = g_lepton2.Px(); + double Lep2Py = g_lepton2.Py(); + double Lep2Pz = g_lepton2.Pz(); + + double llnnPx = Lep1Px + Lep2Px + Nu1Px + Nu2Px; + double llnnPy = Lep1Py + Lep2Py + Nu1Py + Nu2Py; + double llnnPz = Lep1Pz + Lep2Pz + Nu1Pz + Nu2Pz; + double llnnE = Lep1E + Lep2E + Nu1E + Nu2E; + + double mllnn = GetMass( llnnPx, llnnPy, llnnPz, llnnE ); + + double deltaLLNN = (mllnn*mllnn - mh*mh) / (sigmahlep*sigmahlep); + + double ln11Px = Lep1Px + Nu1Px; + double ln11Py = Lep1Py + Nu1Py; + double ln11Pz = Lep1Pz + Nu1Pz; + double ln11E = Lep1E + Nu1E; + + double mln11 = GetMass( ln11Px, ln11Py, ln11Pz, ln11E ); + + double deltaLN11on = (mln11*mln11 - mw*mw) / (sigmaWon*sigmaWon); + + double deltaLN11off = (mln11*mln11 - mhw*mhw) / (sigmaWoff*sigmaWoff); + + double ln22Px = Lep2Px + Nu2Px; + double ln22Py = Lep2Py + Nu2Py; + double ln22Pz = Lep2Pz + Nu2Pz; + double ln22E = Lep2E + Nu2E; + + double mln22 = GetMass( ln22Px, ln22Py, ln22Pz, ln22E ); + + double deltaLN22on = (mln22*mln22 - mw*mw) / (sigmaWon*sigmaWon); + + double deltaLN22off = (mln22*mln22 - mhw*mhw) / (sigmaWoff*sigmaWoff); + + double ln12Px = Lep1Px + Nu2Px; + double ln12Py = Lep1Py + Nu2Py; + double ln12Pz = Lep1Pz + Nu2Pz; + double ln12E = Lep1E + Nu2E; + + double mln12 = GetMass( ln12Px, ln12Py, ln12Pz, ln12E ); + + double deltaLN12on = (mln12*mln12 - mw*mw) / (sigmaWoff*sigmaWoff); + double deltaLN12off = (mln12*mln12 - mhw*mhw) / (sigmaWoff*sigmaWoff); + + double ln21Px = Lep2Px + Nu1Px; + double ln21Py = Lep2Py + Nu1Py; + double ln21Pz = Lep2Pz + Nu2Pz; + double ln21E = Lep2E + Nu1E; + + double mln21 = GetMass( ln21Px, ln21Py, ln21Pz, ln21E ); + + double deltaLN21on = (mln21*mln21 - mw*mw) / (sigmaWon*sigmaWon); + double deltaLN21off = (mln21*mln21 - mhw*mhw) / (sigmaWoff*sigmaWoff); + + double deltaCase1 = deltaLN11on*deltaLN11on + deltaLN22off*deltaLN22off; + double deltaCase2 = deltaLN12on*deltaLN12on + deltaLN21off*deltaLN21off; + double deltaCase3 = deltaLN22on*deltaLN22on + deltaLN11off*deltaLN11off; + double deltaCase4 = deltaLN21on*deltaLN21on + deltaLN12off*deltaLN12off; + + double deltaCaseTotal[] = {deltaCase1, deltaCase2, deltaCase3, deltaCase4}; + double deltaCaseMin = *min_element(deltaCaseTotal, deltaCaseTotal+4); + + //double nnPx = Nu1Px + Nu2Px; + //double nnPy = Nu1Py + Nu2Py; + //double nnPz = Nu1Pz + Nu2Pz; + //double nnE = Nu1E + Nu2E; + + //double mnn = GetMass( nnPx, nnPy, nnPz, nnE ); + + chisq = deltaLLNN*deltaLLNN + deltaCaseMin; + + f = chisq; +} + +float GetHiggsness() { + float higgsness = 0; + ptMinuit = new TMinuit(n_unknown_par); + ptMinuit->SetPrintLevel(-1); + ptMinuit->SetFCN(fcnH); + minuit_arglist[0] = 1; + ptMinuit->mnexcm("SET ERR", minuit_arglist, 1, ierflag); + Double_t vstart[4] = {0.0, 0.0, 0.0, 0.0}; + Double_t step[4] = {0.1, 0.1, 0.1, 0.1}; + TString par_name[4] = {"Nu1Px", "Nu1Py", "Nu1Pz", "Nu2Pz"}; + for (int i=0; i<4; i++) { + ptMinuit->mnparm(i, par_name[i], vstart[i], step[i], minuit_pmin, minuit_pmax, ierflag); + } + + minuit_arglist[0] = 500; + minuit_arglist[1] = 1.; + + ptMinuit->mnexcm("MIGRAD", minuit_arglist, 2, ierflag); + + double outparH[n_unknown_par], errH[n_unknown_par]; + for (int i = 0; i<4; i++) { + ptMinuit->GetParameter(i, outparH[i], errH[i]); + } + + Double_t aminH, edmH, errdefH; + Int_t nvparH, nparxH, icstatH; + + ptMinuit->mnstat(aminH, edmH, errdefH, nvparH, nparxH, icstatH); + if ( aminH > 0 ) { + higgsness = log10(aminH); + } + delete ptMinuit; + return higgsness; +} + +float GetTopness() { + float topness = 0; + ptMinuitT1 = new TMinuit(n_unknown_par); + ptMinuitT1->SetPrintLevel(-1); + ptMinuitT1->SetFCN(fcnT1); + minuit_arglist[0] = 1; + ptMinuitT1->mnexcm("SET ERR", minuit_arglist, 1, ierflag); + Double_t vstart[4] = {0.0, 0.0, 0.0, 0.0}; + Double_t step[4] = {0.1, 0.1, 0.1, 0.1}; + TString par_name[4] = {"Nu1Px", "Nu1Py", "Nu1Pz", "Nu2Pz"}; + for (int i=0; i<4; i++) { + ptMinuitT1->mnparm(i, par_name[i], vstart[i], step[i], minuit_pmin, minuit_pmax, ierflag); + } + + minuit_arglist[0] = 500; + minuit_arglist[1] = 1.; + + ptMinuitT1->mnexcm("MIGRAD", minuit_arglist, 2, ierflag); + + double outparT1[n_unknown_par], errT1[n_unknown_par]; + for (int i = 0; i<4; i++) { + ptMinuitT1->GetParameter(i, outparT1[i], errT1[i]); + } + + Double_t aminT1, edmT1, errdefT1; + Int_t nvparT1, nparxT1, icstatT1; + + ptMinuitT1->mnstat(aminT1, edmT1, errdefT1, nvparT1, nparxT1, icstatT1); + + delete ptMinuitT1; + + ptMinuitT2 = new TMinuit(n_unknown_par); + ptMinuitT2->SetPrintLevel(-1); + ptMinuitT2->SetFCN(fcnT2); + minuit_arglist[0] = 1; + ptMinuitT2->mnexcm("SET ERR", minuit_arglist, 1, ierflag); + //vstart[4] = {0.0, 0.0, 0.0, 0.0}; + //step[4] = {0.1, 0.1, 0.1, 0.1}; + //par_name[4] = {"Nu1Px", "Nu1Py", "Nu1Pz", "Nu2Pz"}; + for (int i=0; i<4; i++) { + ptMinuitT2->mnparm(i, par_name[i], vstart[i], step[i], minuit_pmin, minuit_pmax, ierflag); + } + + minuit_arglist[0] = 500; + minuit_arglist[1] = 1.; + + ptMinuitT2->mnexcm("MIGRAD", minuit_arglist, 2, ierflag); + + double outparT2[n_unknown_par], errT2[n_unknown_par]; + for (int i = 0; i<4; i++) { + ptMinuitT2->GetParameter(i, outparT2[i], errT2[i]); + } + + Double_t aminT2, edmT2, errdefT2; + Int_t nvparT2, nparxT2, icstatT2; + + ptMinuitT2->mnstat(aminT2, edmT2, errdefT2, nvparT2, nparxT2, icstatT2); + + delete ptMinuitT2; + + Double_t aminT; + aminT = min(aminT1, aminT2); + + if( aminT > 0 ) + { + topness = log10(aminT); + } + return topness; +}