Browse Source

split out header files, bdt train and eval, roo fitting

pull/1/head
Marius Pfeiffer 7 months ago
parent
commit
8733cb7027
  1. 3
      .gitignore
  2. 232
      basic_analysis.h
  3. 269
      bdt_classification.h
  4. 19
      constants.h
  5. 239
      hlt1_decision_analysis.h
  6. 524
      new_analysis.cpp

3
.gitignore

@ -3,4 +3,5 @@ dataloader/**
status_report/**
branchnames.txt
images/**
output_files/**
output_files/**
BuToHpMuMu_dataloader/**

232
basic_analysis.h

@ -0,0 +1,232 @@
#ifndef BASIC_ANALYSIS
#define BASIC_ANALYSIS
#include <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <filesystem>
#include <string_view>
#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

269
bdt_classification.h

@ -0,0 +1,269 @@
#ifndef BDT_CLASSIFICATION
#define BDT_CLASSIFICATION
#include <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <filesystem>
#include <string_view>
#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<TV *> 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<TV *> 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<TV *> 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

19
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

239
hlt1_decision_analysis.h

@ -0,0 +1,239 @@
#ifndef HLT1_DECISION_ANALYSIS
#define HLT1_DECISION_ANALYSIS
#include <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <filesystem>
#include <string_view>
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<Hlt1Decision> 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<std::string> lines;
};
const std::vector<Hlt1DecisionPlot>
Hlt1DecisionSets{
Hlt1DecisionPlot{"mva_excl", true, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}},
Hlt1DecisionPlot{"mva_incl", false, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}},
Hlt1DecisionPlot{"pt_excl", true, std::set<std::string>{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}},
Hlt1DecisionPlot{"pt_incl", false, std::set<std::string>{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}},
Hlt1DecisionPlot{"dimu_excl", true, std::set<std::string>{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}},
Hlt1DecisionPlot{"dimu_incl", false, std::set<std::string>{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}},
Hlt1DecisionPlot{"mvapt_incl", false, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}},
Hlt1DecisionPlot{"mvadimu_incl", false, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}},
Hlt1DecisionPlot{"mvadimupt_incl", false, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "LowPtMuon", "LowPtDiMuon", "DisplacedLeptons", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}},
Hlt1DecisionPlot{"ptdimu_incl", false, std::set<std::string>{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons", "DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}},
};
bool CutHlt1DecisionsAnd(const std::set<std::string> &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<std::string> &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<std::string> &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<std::string, int> &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<TH1D *> CreateHlt1DecisionHistos(const char *analysis_name)
{
std::vector<TH1D *> 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<TH1D *> 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<TH1D *> 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

524
new_analysis.cpp

@ -1,3 +1,10 @@
#include <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <filesystem>
#include <string_view>
#include "TH1D.h"
#include "TH2D.h"
#include "THStack.h"
@ -28,397 +35,162 @@
#include "RooFFTConvPdf.h"
#include "RooNovosibirsk.h"
#include <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <filesystem>
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<Hlt1Decision> 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<std::string> lines;
};
const std::vector<Hlt1DecisionPlot>
Hlt1DecisionSets{
Hlt1DecisionPlot{"mva_excl", true, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}},
Hlt1DecisionPlot{"mva_incl", false, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA"}},
Hlt1DecisionPlot{"pt_excl", true, std::set<std::string>{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}},
Hlt1DecisionPlot{"pt_incl", false, std::set<std::string>{"LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}},
Hlt1DecisionPlot{"dimu_excl", true, std::set<std::string>{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}},
Hlt1DecisionPlot{"dimu_incl", false, std::set<std::string>{"DiMuonNoIP", "DiMuonLowMass", "DiMuonHighMass"}},
Hlt1DecisionPlot{"mvapt_incl", false, std::set<std::string>{"TwoTrackMVA", "TrackMuonMVA", "TrackMVA", "LowPtMuon", "LowPtDiMuon", "DisplacedLeptons"}},
Hlt1DecisionPlot{"mvadimu_incl", false, std::set<std::string>{"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<std::string>{"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<std::string>{"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<std::string> &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<std::string> &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<std::string> &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<std::string, int> 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<TV *> 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<std::string, int> &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<TH1D *> CreateHlt1DecisionHistos(const char *analysis_name)
{
std::vector<TH1D *> 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<TH1D *> 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<TH1D *> 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<std::string, int> 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);

Loading…
Cancel
Save