diff --git a/.gitignore b/.gitignore index df8967a..38ba3e8 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ dataloader/** status_report/** branchnames.txt images/** -output_files/** \ No newline at end of file +output_files/** +BuToHpMuMu_dataloader/** \ No newline at end of file diff --git a/basic_analysis.h b/basic_analysis.h new file mode 100644 index 0000000..5edba4b --- /dev/null +++ b/basic_analysis.h @@ -0,0 +1,232 @@ + +#ifndef BASIC_ANALYSIS +#define BASIC_ANALYSIS + +#include +#include +#include +#include +#include +#include + +#include "TTree.h" + +class FourVect +{ +private: + Float_t px, py, pz, energy; + + Float_t *PtrPX() + { + return &px; + } + + Float_t *PtrPY() + { + return &py; + } + + Float_t *PtrPZ() + { + return &pz; + } + + Float_t *PtrEnergy() + { + return &energy; + } + +public: + static FourVect *Init(TTree *tree, const char *particle) + { + FourVect *v = new FourVect{}; + tree->SetBranchAddress(TString::Format("%s_PX", particle).Data(), v->PtrPX()); + tree->SetBranchAddress(TString::Format("%s_PY", particle).Data(), v->PtrPY()); + tree->SetBranchAddress(TString::Format("%s_PZ", particle).Data(), v->PtrPZ()); + tree->SetBranchAddress(TString::Format("%s_ENERGY", particle).Data(), v->PtrEnergy()); + return v; + } + + TLorentzVector LorentzVector() + { + return TLorentzVector(px, py, pz, energy); + } + + TLorentzVector LorentzVector(double sub_mass) + { + TVector3 momentum(px, py, pz); + float energy = TMath::Sqrt(TMath::Sq(sub_mass) + momentum.Mag2()); + return TLorentzVector(momentum, energy); + } + + std::string ToString() + { + return TString::Format("(%f, %f, %f, %f)", px, py, pz, energy).Data(); + } +}; + +void DrawInDefaultCanvas(TH1 *histogram, const char *folder, double margin_left = 0, Option_t *option = "") +{ + std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); + TString name = TString::Format("%s_canvas", histogram->GetName()); + TCanvas *c = new TCanvas(name, histogram->GetName(), 0, 0, 800, 600); + c->SetLeftMargin(margin_left); + histogram->SetStats(0); + histogram->Draw(option); + c->Draw(); + c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); +} + +void DrawInDefaultCanvas(RooPlot *rooPlot, const char *folder, double margin_left = 0, Option_t *option = "") +{ + std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); + TString name = TString::Format("%s_canvas", rooPlot->GetName()); + TCanvas *c = new TCanvas(name, rooPlot->GetName(), 0, 0, 800, 600); + c->SetLeftMargin(margin_left); + rooPlot->Draw(option); + c->Draw(); + c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); +} + +void PrintProgress(const char *title, unsigned int total, unsigned int every, unsigned int current) +{ + if ((current + 1) % every == 0 || current + 1 == total) + { + std::cout << "[" + << title + << "] Processed event: " << current + 1 << " (" << TString::Format("%.2f", ((double)(current + 1) / (double)total) * 100.) << "%)" << std::endl; + } +} + +RooPlot *CreateRooFitHistogram(TH1D *hist) +{ + RooRealVar roo_var_mass("var_mass", "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + roo_var_mass.setRange("fitting_range", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + + RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); + + RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(hist->GetTitle()), RooFit::Name(TString::Format("%s_rplt", hist->GetName()))); + roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(B_MASS_HIST_BINS), RooFit::Name("B Mass Distribution")); + roo_frame_mass->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle()); + + return roo_frame_mass; +} + +struct ShapeParamters +{ + double alpha_left; + double n_left; + double alpha_right; + double n_right; +}; + +struct RooFitSummary +{ + RooPlot *fit_histogram; + RooPlot *pull_histogram; + ShapeParamters shape_parameters; +}; + +RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool useExtShape, ShapeParamters extShape) +{ + auto suffix_name = [name = hist->GetName()](const char *text) + { + return TString::Format("%s_%s", text, name); + }; + + RooRealVar roo_var_mass(suffix_name("var_mass"), "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + const char *fitting_range_name = "fitting_range"; + roo_var_mass.setRange(fitting_range_name, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + + TString hist_name = suffix_name("roohist_B_M"); + RooDataHist roohist_B_M(hist_name, "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); + + RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(hist->GetTitle()), RooFit::Name(TString::Format("%s_rplt", hist->GetName()))); + roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(B_MASS_HIST_BINS), RooFit::Name(hist_name)); + roo_frame_mass->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle()); + + // Crystal Ball for Signal + RooRealVar var_mass_x0(suffix_name("var_mass_x0"), "#mu", 5278., 5170., 5500.); + RooRealVar var_mass_sigmaLR(suffix_name("var_mass_sigmaLR"), "#sigma_{LR}", 16., 5., 40.); + + RooRealVar var_mass_alphaL(suffix_name("var_mass_alphaL"), "#alpha_{L}", 2., 0., 4.); + RooRealVar var_mass_nL(suffix_name("var_mass_nL"), "n_{L}", 5., 0., 15.); + RooRealVar var_mass_alphaR(suffix_name("var_mass_alphaR"), "#alpha_{R}", 2., 0., 4.); + RooRealVar var_mass_nR(suffix_name("var_mass_nR"), "n_{R}", 5., 0., 15.); + + if (useExtShape) { + var_mass_alphaL.setConstant(true); + var_mass_alphaL.setVal(extShape.alpha_left); + + var_mass_nL.setConstant(true); + var_mass_nL.setVal(extShape.n_left); + + var_mass_alphaR.setConstant(true); + var_mass_alphaR.setVal(extShape.alpha_right); + + var_mass_nR.setConstant(true); + var_mass_nR.setVal(extShape.n_left); + } + + TString signal_name = suffix_name("sig_cb"); + RooCrystalBall sig_cb(signal_name, "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); + + TString pull_compare_name{}; + if (hasExpBkg) + { + // Exponential for Background + RooRealVar var_mass_bkg_c(suffix_name("var_mass_bkg_c"), "#lambda", -0.0014, -0.004, -0.000); + + TString background_name = suffix_name("bgk_exp"); + RooExponential bkg_exp(background_name, "Exp Background", roo_var_mass, var_mass_bkg_c); + + RooRealVar var_mass_nsig(suffix_name("nsig"), "N_{Sig}", 0., hist->GetEntries()); + RooRealVar var_mass_nbkg(suffix_name("nbkg"), "N_{Bkg}", 0., hist->GetEntries()); + + TString fitted_pdf_name = suffix_name("sigplusbkg"); + RooAddPdf fitted_pdf(fitted_pdf_name, "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + + RooFitResult *fitres = fitted_pdf.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); + + // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + fitted_pdf.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range(fitting_range_name), RooFit::Name(fitted_pdf_name)); + fitted_pdf.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range(fitting_range_name), RooFit::Name(background_name)); + fitted_pdf.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range(fitting_range_name), RooFit::Name(signal_name)); + fitted_pdf.paramOn(roo_frame_mass); + + pull_compare_name = fitted_pdf_name; + } + else + { + RooFitResult *fitres = sig_cb.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); + + // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + sig_cb.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range(fitting_range_name), RooFit::Name(signal_name)); + pull_compare_name = signal_name; + sig_cb.paramOn(roo_frame_mass); + } + + RooPlot *roo_frame_pull = roo_var_mass.frame(RooFit::Title("Pull Distribution"), RooFit::Name(TString::Format("%s_pull_rplt", hist->GetName()))); + roo_frame_pull->addPlotable(roo_frame_mass->pullHist(hist_name, pull_compare_name), "P"); + + return RooFitSummary{ + roo_frame_mass, + roo_frame_pull, + ShapeParamters{ + var_mass_alphaL.getVal(), + var_mass_nL.getVal(), + var_mass_alphaR.getVal(), + var_mass_nR.getVal()}}; +} + +bool InRange(double value, double center, double low_intvl, double up_intvl) +{ + return center - low_intvl < value && value < center + up_intvl; +} + +bool InRange(double value, double center, double intvl) +{ + return InRange(value, center, intvl, intvl); +} + +#endif \ No newline at end of file diff --git a/bdt_classification.h b/bdt_classification.h new file mode 100644 index 0000000..b83fad0 --- /dev/null +++ b/bdt_classification.h @@ -0,0 +1,269 @@ +#ifndef BDT_CLASSIFICATION +#define BDT_CLASSIFICATION + +#include +#include +#include +#include +#include +#include + +#include "RtypesCore.h" + +enum TVT +{ + Double, + Float, + Int +}; + +class TV +{ +private: + std::string data_name; + std::string mc_name; + std::string train_name; + TVT type; + Double_t mc_double_value; + Float_t mc_float_value; + Double_t data_double_value; + Float_t data_float_value; + + TV(std::string data_name, std::string mc_name, std::string train_name, TVT type) + : data_name{data_name}, mc_name{mc_name}, train_name{train_name}, type{type} + { + } + +public: + static TV *Float(std::string data_name, std::string mc_name, std::string train_name) + { + return new TV(data_name, mc_name, train_name, TVT::Float); + } + + static TV *Float(std::string data_name, std::string mc_name) + { + return new TV(data_name, mc_name, data_name, TVT::Float); + } + + static TV *Float(std::string data_name) + { + return new TV(data_name, data_name, data_name, TVT::Float); + } + + static TV *Double(std::string data_name, std::string mc_name, std::string train_name) + { + return new TV(data_name, mc_name, train_name, TVT::Double); + } + + static TV *Double(std::string data_name, std::string mc_name) + { + return new TV(data_name, mc_name, data_name, TVT::Double); + } + + static TV *Double(std::string data_name) + { + return new TV(data_name, data_name, data_name, TVT::Double); + } + + const char *GetDataName() + { + return data_name.c_str(); + } + + const char *GetMCName() + { + return mc_name.c_str(); + } + + const char *GetTrainName() + { + return train_name.c_str(); + } + + Double_t *GetMCDoubleRef() + { + return &mc_double_value; + } + + Float_t *GetMCFloatRef() + { + return &mc_float_value; + } + + Double_t *GetDataDoubleRef() + { + return &data_double_value; + } + + Float_t *GetDataFloatRef() + { + return &data_float_value; + } + + Double_t GetDataDouble() + { + return data_double_value; + } + + Float_t GetDataFloat() + { + return data_float_value; + } + + void PrintDataValue(int entry) + { + std::cout << data_name << " (" << entry << "): "; + if (IsDouble()) + { + std::cout << data_double_value; + } + else if (IsFloat()) + { + std::cout << data_float_value; + } + + std::cout << std::endl; + } + + void PrintMCValue(int entry) + { + std::cout << mc_name << " (" << entry << "): "; + if (IsDouble()) + { + std::cout << mc_double_value; + } + else if (IsFloat()) + { + std::cout << mc_float_value; + } + + std::cout << std::endl; + } + + bool IsDataFinite() + { + if (IsDouble()) + { + return std::isfinite(data_double_value); + } + else if (IsFloat()) + { + return std::isfinite(data_float_value); + } + + return false; + } + + bool IsMCFinite() + { + if (IsDouble()) + { + return std::isfinite(mc_double_value); + } + else if (IsFloat()) + { + return std::isfinite(mc_float_value); + } + + return false; + } + + bool IsDouble() + { + return type == TVT::Double; + } + + bool IsFloat() + { + return type == TVT::Float; + } +}; + +void ConnectVarsToData(std::vector vars, TChain *data_chain, TChain *mc_chain, TTree *sig_tree, TTree *bkg_tree) +{ + for (size_t i = 0; i < vars.size(); i++) + { + if (vars[i]->IsDouble()) + { + data_chain->SetBranchAddress(vars[i]->GetDataName(), vars[i]->GetDataDoubleRef()); + mc_chain->SetBranchAddress(vars[i]->GetMCName(), vars[i]->GetMCDoubleRef()); + sig_tree->Branch(vars[i]->GetTrainName(), vars[i]->GetMCDoubleRef(), TString::Format("%s/D", vars[i]->GetTrainName())); + bkg_tree->Branch(vars[i]->GetTrainName(), vars[i]->GetDataDoubleRef(), TString::Format("%s/D", vars[i]->GetTrainName())); + } + else if (vars[i]->IsFloat()) + { + data_chain->SetBranchAddress(vars[i]->GetDataName(), vars[i]->GetDataFloatRef()); + mc_chain->SetBranchAddress(vars[i]->GetMCName(), vars[i]->GetMCFloatRef()); + sig_tree->Branch(vars[i]->GetTrainName(), vars[i]->GetMCFloatRef(), TString::Format("%s/F", vars[i]->GetTrainName())); + bkg_tree->Branch(vars[i]->GetTrainName(), vars[i]->GetDataFloatRef(), TString::Format("%s/F", vars[i]->GetTrainName())); + } + } +} + +void TrainBDT(std::vector vars, const char* unique_id, TTree *sig_tree, TTree *bkg_tree) +{ + TString outfile_name = TString::Format("%s_out.root", unique_id); + TFile *output_file = TFile::Open(outfile_name, "RECREATE"); + + TString factory_options("V:Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D:AnalysisType=Auto"); + + TMVA::Factory *factory = new TMVA::Factory(TString::Format("%s_factory", unique_id), output_file, factory_options); + TMVA::DataLoader *data_loader = new TMVA::DataLoader(TString::Format("%s_dataloader", unique_id)); + + for (int i = 0; i < vars.size(); i++) + { + std::cout << "@TMVA: Adding Branch: " << vars[i]->GetTrainName() << std::endl; + if (vars[i]->IsDouble()) + { + data_loader->AddVariable(vars[i]->GetTrainName(), 'D'); + } + else if (vars[i]->IsFloat()) + { + data_loader->AddVariable(vars[i]->GetTrainName(), 'F'); + } + } + + Double_t signal_weight = 1.0, background_weight = 1.0; + + data_loader->AddSignalTree(sig_tree, signal_weight); + data_loader->AddBackgroundTree(bkg_tree, background_weight); + data_loader->PrepareTrainingAndTestTree("", "", "nTrain_Signal=0:nTrain_Background=0:SplitMode=Random:NormMode=EqualNumEvents:!V"); + + factory->BookMethod(data_loader, TMVA::Types::kBDT, "BDT", "!H:!V:NTrees=600:MinNodeSize=2.5%:CreateMVAPdfs:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:UseBaggedBoost:BaggedSampleFraction=0.5:SeparationType=GiniIndex:nCuts=20"); + + factory->TrainAllMethods(); + factory->TestAllMethods(); + factory->EvaluateAllMethods(); + + output_file->Close(); +} + +TMVA::Reader* SetupReader(std::vector vars, Float_t* train_vars, const char* unique_id) { + TMVA::Reader *reader = new TMVA::Reader("!Color:!Silent"); + + for (size_t i = 0; i < vars.size(); i++) + { + reader->AddVariable(vars[i]->GetTrainName(), &train_vars[i]); + } + + reader->BookMVA("BDT", TString::Format("./%s_dataloader/weights/%s_factory_BDT.weights.xml", unique_id, unique_id)); + + return reader; +} + +void DrawBDTProbs(TH1D *histogram, const double cut_value, const char *folder) +{ + std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); + TString name = TString::Format("%s_canvas", histogram->GetName()); + TCanvas *c = new TCanvas(name, histogram->GetName(), 0, 0, 800, 600); + histogram->SetStats(0); + histogram->Draw(); + TLine* line = new TLine(cut_value, 0, cut_value, histogram->GetMaximum()); + line->SetLineColor(kRed); + line->SetLineStyle(kDashed); + line->Draw(); + c->Draw(); + c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); +} + +#endif \ No newline at end of file diff --git a/constants.h b/constants.h new file mode 100644 index 0000000..9e481e8 --- /dev/null +++ b/constants.h @@ -0,0 +1,19 @@ +#ifndef CONSTANTS +#define CONSTANTS + +#include "RtypesCore.h" + +const Double_t B_MASS_VAR_MIN = 5100.; +const Double_t B_MASS_VAR_MAX = 6000.; +const Int_t B_MASS_HIST_BINS = 70; + +const Double_t K_MASS = 493.677; +const Double_t PSI2S_MASS = 3686.093; +const Double_t JPSI_MASS = 3096.900; +const Double_t KSTAR_MASS = 891.760; + +const int PID_KAON = 321; +const int PID_PION = 211; +const int PID_MUON = 13; + +#endif \ No newline at end of file diff --git a/hlt1_decision_analysis.h b/hlt1_decision_analysis.h new file mode 100644 index 0000000..2ae9ffe --- /dev/null +++ b/hlt1_decision_analysis.h @@ -0,0 +1,239 @@ +#ifndef HLT1_DECISION_ANALYSIS +#define HLT1_DECISION_ANALYSIS + +#include +#include +#include +#include +#include +#include + + +struct Hlt1Decision +{ + std::string name; + int index; + Bool_t value; + + std::string GetName() const + { + return TString::Format("Hlt1%sDecision", name.c_str()).Data(); + } + + Bool_t *GetValuePointer() + { + return &value; + } +}; + +std::vector Hlt1Decisions{ + Hlt1Decision{"DiMuonHighMass", 1}, + Hlt1Decision{"DiMuonLowMass", 2}, + Hlt1Decision{"DiMuonNoIP", 3}, + Hlt1Decision{"DiMuonSoft", 4}, + Hlt1Decision{"DisplacedLeptons", 5}, + Hlt1Decision{"LowPtDiMuon", 6}, + Hlt1Decision{"LowPtMuon", 7}, + Hlt1Decision{"OneMuonTrackLine", 8}, + Hlt1Decision{"SingleHighEt", 9}, + Hlt1Decision{"SingleHighPtMuon", 10}, + Hlt1Decision{"TrackMVA", 11}, + Hlt1Decision{"TrackMuonMVA", 12}, + Hlt1Decision{"TwoTrackMVA", 13}, +}; + +struct Hlt1DecisionPlot +{ + std::string name; + bool exclusive; + std::set lines; +}; + +const std::vector + Hlt1DecisionSets{ + Hlt1DecisionPlot{"mva_excl", true, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}}, + + Hlt1DecisionPlot{"mva_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}}, + + Hlt1DecisionPlot{"pt_excl", true, std::set{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}}, + + Hlt1DecisionPlot{"pt_incl", false, std::set{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}}, + + Hlt1DecisionPlot{"dimu_excl", true, std::set{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + + Hlt1DecisionPlot{"dimu_incl", false, std::set{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + + Hlt1DecisionPlot{"mvapt_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}}, + + Hlt1DecisionPlot{"mvadimu_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + + Hlt1DecisionPlot{"mvadimupt_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "LowPtMuon", "LowPtDiMuon", "DisplacedLeptons", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + + Hlt1DecisionPlot{"ptdimu_incl", false, std::set{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + + }; + +bool CutHlt1DecisionsAnd(const std::set &decision_list) +{ + bool okay = true; + for (const auto &var : Hlt1Decisions) + { + if (decision_list.find(var.name) != decision_list.end()) + { + okay = okay && var.value; + } + } + return okay; +} + +bool CutHlt1DecisionsOr(const std::set &decision_list) +{ + bool okay = false; + for (const auto &var : Hlt1Decisions) + { + if (decision_list.find(var.name) != decision_list.end()) + { + okay = okay || var.value; + } + } + return okay; +} + +bool CutHlt1DecisionsOrOnly(const std::set &decision_list) +{ + bool okay = false; + for (const auto &var : Hlt1Decisions) + { + if (decision_list.find(var.name) != decision_list.end()) + { + okay = okay || var.value; + } + else + { + if (var.value) + { + okay = false; + break; + } + } + } + return okay; +} + +void ConnectHlt1Decisions(TChain *chain, TH2D *incl_hist, TH2D *excl_hist) +{ + for (auto &var : Hlt1Decisions) + { + if (chain->FindBranch(var.GetName().c_str())) { + chain->SetBranchAddress(var.GetName().c_str(), var.GetValuePointer()); + incl_hist->Fill(0., var.name.c_str(), 0.); + excl_hist->Fill(0., var.name.c_str(), 0.); + } + } +} + +void CheckHlt1Decisioins(TH2D *incl_hist, TH2D *excl_hist, std::map &exclusive_hits, const double reco_mass) +{ + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + incl_hist->Fill(reco_mass, var.name.c_str(), 1); + } + } + + bool exclusive = true; + std::string line{}; + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + if (!line.empty()) + { + exclusive = false; + } + + line = var.name; + } + } + + if (!line.empty() && exclusive) + { + int &hits = exclusive_hits[line]; + if (hits) + { + hits += 1; + } + else + { + hits = 1; + } + excl_hist->Fill(reco_mass, line.c_str(), 1); + } +} + +std::vector CreateHlt1DecisionHistos(const char *analysis_name) +{ + std::vector histos{}; + for (int i = 0; i < Hlt1DecisionSets.size(); i++) + { + TH1D *h1_B_Mass_hlt1 = new TH1D(TString::Format("h1_B_Mass_hlt1_%s", Hlt1DecisionSets[i].name.c_str()), TString::Format("(%s) (%s)", Hlt1DecisionSets[i].name.c_str(), analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + histos.push_back(h1_B_Mass_hlt1); + } + + return histos; +} + +void FillHlt1DecisionHistos(std::vector histos, double reco_mass) +{ + for (int i = 0; i < Hlt1DecisionSets.size(); i++) + { + if (Hlt1DecisionSets[i].exclusive) + { + if (CutHlt1DecisionsOrOnly(Hlt1DecisionSets[i].lines)) + { + histos[i]->Fill(reco_mass); + } + } + else + { + if (CutHlt1DecisionsOr(Hlt1DecisionSets[i].lines)) + { + histos[i]->Fill(reco_mass); + } + } + } +} + +void DrawHlt1DecisionHistos(const char *folder, std::vector histos) +{ + std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); + TString name = TString::Format("%s_canvas", "hlt1_decisions"); + TCanvas *c = new TCanvas(name, "HLT 1 Decisions", 0, 0, 1200, 700); + c->Divide(4, 3); + for (int i = 0; i < histos.size(); i++) + { + c->cd(i + 1); + histos[i]->SetStats(0); + histos[i]->Draw(); + + TLegend *leg1 = new TLegend(0.58, 0.89 - (0.05 * (Hlt1DecisionSets[i].lines.size() + (int)(Hlt1DecisionSets[i].lines.size() / 3))), 0.88, 0.89); + leg1->SetFillStyle(0); + leg1->SetLineStyle(0); + leg1->SetMargin(0.1); + int index = 1; + std::for_each(Hlt1DecisionSets[i].lines.begin(), Hlt1DecisionSets[i].lines.end(), [leg1, i, &index](std::string s) + { + leg1->AddEntry((TObject *)0, s.c_str(), ""); + if (index != Hlt1DecisionSets[i].lines.size() && index % 3 == 0) { + leg1->AddEntry((TObject *)0, "", ""); + } + index++; }); + leg1->AddEntry((TObject *)0, TString::Format("# Entries: %d", (int)histos[i]->GetEntries()), ""); + leg1->Draw(); + } + c->Draw(); + c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); +} + +#endif \ No newline at end of file diff --git a/new_analysis.cpp b/new_analysis.cpp index ed05062..d2ba7a1 100644 --- a/new_analysis.cpp +++ b/new_analysis.cpp @@ -1,3 +1,10 @@ +#include +#include +#include +#include +#include +#include + #include "TH1D.h" #include "TH2D.h" #include "THStack.h" @@ -28,397 +35,162 @@ #include "RooFFTConvPdf.h" #include "RooNovosibirsk.h" -#include -#include -#include -#include -#include - -const Double_t B_MASS_VAR_MIN = 5100.; -const Double_t B_MASS_VAR_MAX = 6000.; -const Int_t B_MASS_HIST_BINS = 70; - -const Double_t K_MASS = 493.677; -const Double_t PSI2S_MASS = 3686.093; -const Double_t JPSI_MASS = 3096.900; -const Double_t KSTAR_MASS = 891.760; - -const int PID_KAON = 321; -const int PID_PION = 211; -const int PID_MUON = 13; - -class FourVect -{ -private: - Float_t px, py, pz, energy; - - Float_t *PtrPX() - { - return &px; - } - - Float_t *PtrPY() - { - return &py; - } - - Float_t *PtrPZ() - { - return &pz; - } - - Float_t *PtrEnergy() - { - return &energy; - } - -public: - static FourVect *Init(TTree *tree, const char *particle) - { - FourVect *v = new FourVect{}; - tree->SetBranchAddress(TString::Format("%s_PX", particle).Data(), v->PtrPX()); - tree->SetBranchAddress(TString::Format("%s_PY", particle).Data(), v->PtrPY()); - tree->SetBranchAddress(TString::Format("%s_PZ", particle).Data(), v->PtrPZ()); - tree->SetBranchAddress(TString::Format("%s_ENERGY", particle).Data(), v->PtrEnergy()); - return v; - } - - TLorentzVector LorentzVector() - { - return TLorentzVector(px, py, pz, energy); - } - - TLorentzVector LorentzVector(double sub_mass) - { - TVector3 momentum(px, py, pz); - float energy = TMath::Sqrt(TMath::Sq(sub_mass) + momentum.Mag2()); - return TLorentzVector(momentum, energy); - } +#include "constants.h" +#include "basic_analysis.h" +#include "hlt1_decision_analysis.h" +#include "bdt_classification.h" - std::string ToString() - { - return TString::Format("(%f, %f, %f, %f)", px, py, pz, energy).Data(); - } -}; - -struct Hlt1Decision -{ - std::string name; - int index; - Bool_t value; - - std::string GetName() const - { - return TString::Format("Hlt1%sDecision", name.c_str()).Data(); - } +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bu To Hp Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - Bool_t *GetValuePointer() - { - return &value; - } -}; - -std::vector Hlt1Decisions{ - Hlt1Decision{"DiMuonHighMass", 1}, - Hlt1Decision{"DiMuonLowMass", 2}, - Hlt1Decision{"DiMuonNoIP", 3}, - Hlt1Decision{"DiMuonSoft", 4}, - Hlt1Decision{"DisplacedLeptons", 5}, - Hlt1Decision{"LowPtDiMuon", 6}, - Hlt1Decision{"LowPtMuon", 7}, - Hlt1Decision{"OneMuonTrackLine", 8}, - Hlt1Decision{"SingleHighEt", 9}, - Hlt1Decision{"SingleHighPtMuon", 10}, - Hlt1Decision{"TrackMVA", 11}, - Hlt1Decision{"TrackMuonMVA", 12}, - Hlt1Decision{"TwoTrackMVA", 13}, -}; - -struct Hlt1DecisionPlot +int analyze_bu2hpmumu_data() { - std::string name; - bool exclusive; - std::set lines; -}; - -const std::vector - Hlt1DecisionSets{ - Hlt1DecisionPlot{"mva_excl", true, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}}, - - Hlt1DecisionPlot{"mva_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}}, - - Hlt1DecisionPlot{"pt_excl", true, std::set{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}}, - - Hlt1DecisionPlot{"pt_incl", false, std::set{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}}, - - Hlt1DecisionPlot{"dimu_excl", true, std::set{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, - - Hlt1DecisionPlot{"dimu_incl", false, std::set{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, - - Hlt1DecisionPlot{"mvapt_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}}, - - Hlt1DecisionPlot{"mvadimu_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + const char *analysis_name = "BuToHpMuMu"; + const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; + const bool retrain_bdt = false; - Hlt1DecisionPlot{"mvadimupt_incl", false, std::set{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "LowPtMuon", "LowPtDiMuon", "DisplacedLeptons", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + TChain *data_chain = new TChain(TString::Format("%s/DecayTree", analysis_name)); + data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); + data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root"); - Hlt1DecisionPlot{"ptdimu_incl", false, std::set{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}}, + FourVect *l14v_data = FourVect::Init(data_chain, "L1"); + FourVect *l24v_data = FourVect::Init(data_chain, "L2"); + FourVect *hp4v_data = FourVect::Init(data_chain, "Hp"); - }; + TChain *sim_chain = new TChain(TString::Format("%s_noPID/DecayTree", analysis_name)); + sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root"); -bool CutHlt1DecisionsAnd(const std::set &decision_list) -{ - bool okay = true; - for (const auto &var : Hlt1Decisions) - { - if (decision_list.find(var.name) != decision_list.end()) - { - okay = okay && var.value; - } - } - return okay; -} + FourVect *l14v_sim = FourVect::Init(sim_chain, "L1"); + FourVect *l24v_sim = FourVect::Init(sim_chain, "L2"); + FourVect *hp4v_sim = FourVect::Init(sim_chain, "Hp"); -bool CutHlt1DecisionsOr(const std::set &decision_list) -{ - bool okay = false; - for (const auto &var : Hlt1Decisions) - { - if (decision_list.find(var.name) != decision_list.end()) - { - okay = okay || var.value; - } - } - return okay; -} + Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID; -bool CutHlt1DecisionsOrOnly(const std::set &decision_list) -{ - bool okay = false; - for (const auto &var : Hlt1Decisions) - { - if (decision_list.find(var.name) != decision_list.end()) - { - okay = okay || var.value; - } - else - { - if (var.value) - { - okay = false; - break; - } - } - } - return okay; -} + sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); + sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); + sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); + sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); -void ConnectHlt1Decisions(TChain *chain) -{ - for (auto &var : Hlt1Decisions) - { - if (chain->FindBranch(var.GetName().c_str())) - chain->SetBranchAddress(var.GetName().c_str(), var.GetValuePointer()); - } -} + TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); + TH2D *h2_Hlt1_flags_excl_B_Mass = new TH2D("h2_Hlt1_flags_excl_B_Mass", "Excl Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); + TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -2, 2); -void DrawInDefaultCanvas(TH1 *histogram, const char *folder, double margin_left = 0, Option_t *option = "") -{ - std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); - TString name = TString::Format("%s_canvas", histogram->GetName()); - TCanvas *c = new TCanvas(name, histogram->GetName(), 0, 0, 800, 600); - c->SetLeftMargin(margin_left); - histogram->SetStats(0); - histogram->Draw(option); - c->Draw(); - c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); -} + h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); + h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); + h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); -void DrawInDefaultCanvas(RooPlot *rooPlot, const char *folder, double margin_left = 0, Option_t *option = "") -{ - std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); - TString name = TString::Format("%s_canvas", rooPlot->GetName()); - TCanvas *c = new TCanvas(name, rooPlot->GetName(), 0, 0, 800, 600); - c->SetLeftMargin(margin_left); - rooPlot->Draw(option); - c->Draw(); - c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); -} + ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); + auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); + std::map exclusive_hits{}; -void PrintProgress(const char *title, unsigned int total, unsigned int every, unsigned int current) -{ - if ((current + 1) % every == 0 || current + 1 == total) - { - std::cout << "[" - << title - << "] Processed event: " << current + 1 << " (" << TString::Format("%.2f", ((double)(current + 1) / (double)total) * 100.) << "%)" << std::endl; - } -} + std::vector vars{ + TV::Float("B_PT", "B_PT"), + TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"), + TV::Float("B_BPVDIRA", "B_BPVDIRA"), + TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), + // TV::Float("Jpsi_BPVDIRA", "Jpsi_BPVDIRA"), + TV::Float("Jpsi_PT", "Jpsi_PT"), + TV::Float("Hp_BPVIPCHI2", "Hp_BPVIPCHI2"), + TV::Float("Hp_PT", "Hp_PT"), + // TV::Double("Kplus_PID_K", "K_PID_K"), + TV::Double("Hp_PROBNN_K", "Hp_PROBNN_K"), + TV::Float("L1_BPVIPCHI2", "L1_BPVIPCHI2"), + TV::Float("L2_BPVIPCHI2", "L2_BPVIPCHI2"), + }; -RooPlot *CreateRooFitHistogram(TH1D *hist) -{ - RooRealVar roo_var_mass("var_mass", "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); - roo_var_mass.setRange("fitting_range", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); + TTree *bkg_tree = new TTree("TreeB", "tree containing background data"); - RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); + ConnectVarsToData(vars, data_chain, sim_chain, sig_tree, bkg_tree); - RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(hist->GetTitle()), RooFit::Name(TString::Format("%s_rplt", hist->GetName()))); - roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(B_MASS_HIST_BINS), RooFit::Name("B Mass Distribution")); - roo_frame_mass->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle()); + unsigned int data_entries = data_chain->GetEntries(); + unsigned int sim_entries = sim_chain->GetEntries(); - return roo_frame_mass; -} + unsigned int bkg_events = 0; -void CheckHlt1Decisioins(TH2D *incl_hist, TH2D *excl_hist, std::map &exclusive_hits, const double reco_mass) -{ - for (const auto &var : Hlt1Decisions) + if (retrain_bdt) { - if (var.value) + std::cout << "# Starting with BDT retrain." << std::endl; + for (unsigned int i = 0; i < data_entries; i++) { - incl_hist->Fill(reco_mass, var.name.c_str(), 1); - } - } + data_chain->GetEntry(i); + Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector(K_MASS) + l14v_data->LorentzVector() + l24v_data->LorentzVector()).M(); - bool exclusive = true; - std::string line{}; - for (const auto &var : Hlt1Decisions) - { - if (var.value) - { - if (!line.empty()) + if (std::all_of(vars.begin(), vars.end(), [](TV *v) + { return v->IsDataFinite(); })) { - exclusive = false; + if (reconstructed_B_Mass > 5500.) + { + bkg_tree->Fill(); + bkg_events++; + } } - line = var.name; + PrintProgress(TString::Format("%s BKG Collection", analysis_name), data_entries, 10000, i); } } - - if (!line.empty() && exclusive) - { - int &hits = exclusive_hits[line]; - if (hits) - { - hits += 1; - } - else - { - hits = 1; - } - excl_hist->Fill(reco_mass, line.c_str(), 1); + else { + std::cout << "# Starting without BDT retrain." << std::endl; } -} - -std::vector CreateHlt1DecisionHistos(const char *analysis_name) -{ - std::vector histos{}; - for (int i = 0; i < Hlt1DecisionSets.size(); i++) - { - TH1D *h1_B_Mass_hlt1 = new TH1D(TString::Format("h1_B_Mass_hlt1_%s", Hlt1DecisionSets[i].name.c_str()), TString::Format("(%s) (%s)", Hlt1DecisionSets[i].name.c_str(), analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - histos.push_back(h1_B_Mass_hlt1); - } - - return histos; -} -void FillHlt1DecisionHistos(std::vector histos, double reco_mass) -{ - for (int i = 0; i < Hlt1DecisionSets.size(); i++) + for (unsigned int i = 0; i < sim_entries; i++) { - if (Hlt1DecisionSets[i].exclusive) - { - if (CutHlt1DecisionsOrOnly(Hlt1DecisionSets[i].lines)) - { - histos[i]->Fill(reco_mass); - } - } - else + sim_chain->GetEntry(i); + Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON) { - if (CutHlt1DecisionsOr(Hlt1DecisionSets[i].lines)) + if (std::all_of(vars.begin(), vars.end(), [](TV *v) + { return v->IsMCFinite(); })) { - histos[i]->Fill(reco_mass); + if (retrain_bdt) + { + sig_tree->Fill(); + } + h1_B_Mass_sim->Fill(reconstructed_B_Mass); } } + + PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); } -} -void DrawHlt1DecisionHistos(const char *folder, std::vector histos) -{ - std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); - TString name = TString::Format("%s_canvas", "hlt1_decisions"); - TCanvas *c = new TCanvas(name, "HLT 1 Decisions", 0, 0, 1200, 700); - c->Divide(4, 3); - for (int i = 0; i < histos.size(); i++) + if (retrain_bdt) { - c->cd(i + 1); - histos[i]->SetStats(0); - histos[i]->Draw(); - - TLegend *leg1 = new TLegend(0.58, 0.89 - (0.05 * (Hlt1DecisionSets[i].lines.size() + (int)(Hlt1DecisionSets[i].lines.size() / 3))), 0.88, 0.89); - leg1->SetFillStyle(0); - leg1->SetLineStyle(0); - leg1->SetMargin(0.1); - int index = 1; - std::for_each(Hlt1DecisionSets[i].lines.begin(), Hlt1DecisionSets[i].lines.end(), [leg1, i, &index](std::string s) - { - leg1->AddEntry((TObject *)0, s.c_str(), ""); - if (index != Hlt1DecisionSets[i].lines.size() && index % 3 == 0) { - leg1->AddEntry((TObject *)0, "", ""); - } - index++; }); - leg1->AddEntry((TObject *)0, TString::Format("# Entries: %d", (int)histos[i]->GetEntries()), ""); - leg1->Draw(); + TrainBDT(vars, "BuToHpMuMu", sig_tree, bkg_tree); + std::cout << "# Finished BDT retrain." << std::endl; } - c->Draw(); - c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); -} -// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bu To Hp Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + std::cout << "# Starting evaluation of data." << std::endl; -int analyze_bu2hpmumu_data() -{ - const char *analysis_name = "BuToHpMuMu"; - const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; + Float_t *train_vars = new Float_t[vars.size()]; + auto reader = SetupReader(vars, train_vars, "BuToHpMuMu"); - TChain *data_chain = new TChain(TString::Format("%s/DecayTree", analysis_name)); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root"); - - FourVect *l14v = FourVect::Init(data_chain, "L1"); - FourVect *l24v = FourVect::Init(data_chain, "L2"); - FourVect *hp4v = FourVect::Init(data_chain, "Hp"); + const double mva_cut_value = 0; // 0.09; // -0.02; - TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - TH2D *h2_Hlt1_flags_excl_B_Mass = new TH2D("h2_Hlt1_flags_excl_B_Mass", "Excl Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - - h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); - h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); + for (unsigned int i = 0; i < data_entries; i++) + { + data_chain->GetEntry(i); - ConnectHlt1Decisions(data_chain); + TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector(); + Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector(K_MASS) + dimuon).M(); - for (const auto &var : Hlt1Decisions) - { - if (data_chain->FindBranch(var.GetName().c_str())) + if (std::all_of(vars.begin(), vars.end(), [](TV *v) + { return v->IsDataFinite(); })) { - h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); - h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); + for (size_t j = 0; j < vars.size(); j++) + { + if (vars[j]->IsDouble()) + { + train_vars[j] = vars[j]->GetDataDouble(); + } + else if (vars[j]->IsFloat()) + { + train_vars[j] = vars[j]->GetDataFloat(); + } + } } - } - - auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); - - std::map exclusive_hits{}; - - unsigned int entries = data_chain->GetEntries(); - for (unsigned int i = 0; i < entries; i++) - { - data_chain->GetEntry(i); - - TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); - Double_t reconstructed_B_Mass = (hp4v->LorentzVector(K_MASS) + dimuon).M(); if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { @@ -426,7 +198,10 @@ int analyze_bu2hpmumu_data() FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); } - if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "DiMuonHighMass"})) + double mva_response = reader->EvaluateMVA("BDT"); + h1_bdt_probs->Fill(mva_response); + + if (mva_response > mva_cut_value) { if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { @@ -438,7 +213,7 @@ int analyze_bu2hpmumu_data() } } - PrintProgress(analysis_name, entries, 10000, i); + PrintProgress(TString::Format("%s BDT Evaluation", analysis_name), data_entries, 10000, i); } std::cout << "# Exclusive Hits" << std::endl; @@ -452,13 +227,23 @@ int analyze_bu2hpmumu_data() DrawInDefaultCanvas(h1_B_Mass_jpsi, analysis_name, 0.1); DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); - auto roofit_hist_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); - auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); - DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); - DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); + auto roofit_hist_sim = CreateRooFitHistogramAndFitCB(h1_B_Mass_sim, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_jpsi, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + + DrawInDefaultCanvas(roofit_hist_jpsi_fitsum.fit_histogram, analysis_name, 0.1); + DrawInDefaultCanvas(roofit_hist_psi2s_fitsum.fit_histogram, analysis_name, 0.1); + DrawInDefaultCanvas(roofit_hist_sim.fit_histogram, analysis_name, 0.1); DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); + DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name); + + // auto test_bkg_tree_file = TFile::Open("test_bkg_tree.root", "RECREATE"); + // bkg_tree->Write(); + // sig_tree->Write(); + // test_bkg_tree_file->Close(); + return 0; } @@ -491,16 +276,7 @@ int analyze_bu2hpmumu_simulation() h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - ConnectHlt1Decisions(sim_chain); - - for (const auto &var : Hlt1Decisions) - { - if (sim_chain->FindBranch(var.GetName().c_str())) - { - h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); - h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); - } - } + ConnectHlt1Decisions(sim_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); @@ -593,13 +369,7 @@ int analyze_b02hphmmumu() h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - ConnectHlt1Decisions(data_chain); - - for (const auto &var : Hlt1Decisions) - { - h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); - h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); - } + ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); @@ -703,13 +473,7 @@ int analyze_Hlt2RD_BuToKpMuMu() h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - ConnectHlt1Decisions(data_chain); - - for (const auto &var : Hlt1Decisions) - { - h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); - h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); - } + ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); @@ -793,13 +557,7 @@ int analyze_Hlt2RD_B0ToKpPimMuMu() h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - ConnectHlt1Decisions(data_chain); - - for (const auto &var : Hlt1Decisions) - { - h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); - h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); - } + ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name);