Skip to content

Commit

Permalink
minimal fixes to make Zmumumerge run in unit tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@mmusich
mmusich committed Feb 28, 2024
1 parent d862fb4 commit df0dced
Show file tree
Hide file tree
Showing 2 changed files with 69 additions and 23 deletions.
76 changes: 57 additions & 19 deletions Alignment/OfflineValidation/bin/Zmumumerge.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,22 +75,32 @@ class FitOut {
};

FitOut ZMassBinFit_OldTool(TH1D* th1d_input, TString s_name = "zmumu_fitting", TString output_path = "./") {
double xMin(75), xMax(105), xMean(91);
double sigma = 2;
double sigmaMin = 0.1;
double sigmaMax = 10;
// silence messages
RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);

double xMean = 91.1876;
double xMin = th1d_input->GetXaxis()->GetXmin();
double xMax = th1d_input->GetXaxis()->GetXmax();

double sigma(2.);
double sigmaMin(0.1);
double sigmaMax(10.);

double sigma2(0.1);
double sigma2Min(0.);
double sigma2Max(10.);

std::unique_ptr<FitWithRooFit> fitter = std::make_unique<FitWithRooFit>();

bool useChi2(false);

FitWithRooFit* fitter = new FitWithRooFit();
double sigma2(0.1), sigma2Min(0.), sigma2Max(10.), useChi2(false);
fitter->useChi2_ = useChi2;
fitter->initMean(xMean, xMin, xMax);
fitter->initSigma(sigma, sigmaMin, sigmaMax);
fitter->initSigma2(sigma2, sigma2Min, sigma2Max);
fitter->initAlpha(1.5, 0.05, 10.);
fitter->initN(1, 0.01, 100.);
fitter->initFGCB(0.4, 0., 1.);

fitter->initMean(91.1876, xMin, xMax);
fitter->initGamma(2.4952, 0., 10.);
fitter->gamma()->setConstant(kTRUE);
fitter->initMean2(0., -20., 20.);
Expand All @@ -109,6 +119,7 @@ FitOut ZMassBinFit_OldTool(TH1D* th1d_input, TString s_name = "zmumu_fitting", T
fitter->initA4(0., -10., 10.);
fitter->initA5(0., -10., 10.);
fitter->initA6(0., -10., 10.);

TCanvas* c1 = new TCanvas();
c1->Clear();
c1->SetLeftMargin(0.15);
Expand All @@ -121,9 +132,10 @@ FitOut ZMassBinFit_OldTool(TH1D* th1d_input, TString s_name = "zmumu_fitting", T

FitOut fitRes(
fitter->mean()->getValV(), fitter->mean()->getError(), fitter->sigma()->getValV(), fitter->sigma()->getError());

return fitRes;
}
void Draw_th1d(TH1D* th1d_input, TString variable_name) {
void Draw_th1d(TH1D* th1d_input, TString variable_name, TString output_path) {
TCanvas* c = new TCanvas();
c->cd();
gStyle->SetOptStat(0);
Expand All @@ -137,12 +149,13 @@ void Draw_th1d(TH1D* th1d_input, TString variable_name) {
th1d_input->GetYaxis()->SetTitle("Mass mean (GeV)");
th1d_input->Draw();
tlxg->DrawLatexNDC(0.2, 0.8, Form("%s", GT.Data()));
c->Print(Form("%s/fitResultPlot/mass_VS_%s.pdf", GT.Data(), variable_name.Data()));
c->Print(Form("%s/fitResultPlot/mass_VS_%s.pdf", output_path.Data(), variable_name.Data()));
}

const static int variables_number = 8;
const TString tstring_variables_name[variables_number] = {
"CosThetaCS", "DeltaEta", "EtaMinus", "EtaPlus", "PhiCS", "PhiMinus", "PhiPlus", "Pt"};

void Fitting_GetMassmeanVSvariables(TString inputfile_name, TString output_path) {
TH2D* th2d_mass_variables[variables_number];
TFile* inputfile = TFile::Open(inputfile_name.Data());
Expand All @@ -154,7 +167,7 @@ void Fitting_GetMassmeanVSvariables(TString inputfile_name, TString output_path)

gSystem->Exec(Form("mkdir -p %s", output_path.Data()));
gSystem->Exec(Form("mkdir -p %s/fitResultPlot", output_path.Data()));
TFile* outpufile = TFile::Open(Form("%s/fitting_output.root", output_path.Data()), "recreate");
TFile* outputfile = TFile::Open(Form("%s/fitting_output.root", output_path.Data()), "RECREATE");
TH1D* th1d_variables_meanmass[variables_number];
TH1D* th1d_variables_entries[variables_number];
const int variables_rebin[variables_number] = {1, 1, 1, 1, 1, 1, 1, 1};
Expand All @@ -169,8 +182,8 @@ void Fitting_GetMassmeanVSvariables(TString inputfile_name, TString output_path)
if (i == 7 and j > 25) {
continue;
}
cout << "th1d_variables_meanmass[i]->GetNbinsX()=" << th1d_variables_meanmass[i]->GetNbinsX() << endl;
cout << "th2d_mass_variables[i]->GetNbinsY()=" << th2d_mass_variables[i]->GetNbinsY() << endl;
//cout << __LINE__ << " th1d_variables_meanmass[i]->GetNbinsX()=" << th1d_variables_meanmass[i]->GetNbinsX() << endl;
//cout << __LINE__ << " th2d_mass_variables[i]->GetNbinsY()=" << th2d_mass_variables[i]->GetNbinsY() << endl;
th1d_variables_meanmass[i]->SetBinContent(j, 0);
th1d_variables_meanmass[i]->SetBinError(j, 0);

Expand All @@ -181,25 +194,39 @@ void Fitting_GetMassmeanVSvariables(TString inputfile_name, TString output_path)
TString s_name = Form("%s_%d", tstring_variables_name[i].Data(), j);

FitOut fitR = ZMassBinFit_OldTool(th1d_i, s_name, output_path);

th1d_variables_meanmass[i]->SetBinContent(j, fitR.mean);
th1d_variables_meanmass[i]->SetBinError(j, fitR.mean_err);
}

th1d_variables_meanmass[i]->GetXaxis()->SetRangeUser(xaxis_range[i][0], xaxis_range[i][1]);
Draw_th1d(th1d_variables_meanmass[i], tstring_variables_name[i]);
outpufile->cd();
Draw_th1d(th1d_variables_meanmass[i], tstring_variables_name[i], output_path);
outputfile->cd();
th1d_variables_meanmass[i]->Write(th1d_name);

TString th1d_name_entries = Form("th1d_entries_%s", tstring_variables_name[i].Data());
th1d_variables_entries[i] = th2d_mass_variables[i]->ProjectionY(th1d_name_entries, 0, -1, "d");
th1d_variables_entries[i]->GetXaxis()->SetTitle(tstring_variables_name[i].Data());
th1d_variables_entries[i]->GetYaxis()->SetTitle("Entry");
outpufile->cd();
outputfile->cd();
th1d_variables_entries[i]->Write(th1d_name_entries);
}

outpufile->Write();
outpufile->Close();
delete outpufile;
if (outputfile->IsOpen()) {
// Get the path (current working directory) in which the file is going to be written
const char* path = outputfile->GetPath();

if (path) {
std::cout << "File is going to be written in the directory: " << path << " for input file: " << inputfile_name
<< std::endl;
} else {
std::cerr << "Error: Unable to determine the path." << std::endl;
}

//outputfile->Write();
outputfile->Close();
delete outputfile;
}
}

const static int max_file_number = 10;
Expand Down Expand Up @@ -261,7 +288,16 @@ int Zmumumerge(int argc, char* argv[]) {
pt::read_json(options.config, main_tree);
pt::ptree alignments = main_tree.get_child("alignments");
pt::ptree validation = main_tree.get_child("validation");

for (const auto& childTree : alignments) {
// do not consider the nodes with a "file" to merge
if (childTree.second.find("file") == childTree.second.not_found()) {
std::cerr << "Ignoring alignment: " << childTree.second.get<std::string>("title") << ".\nNo file to merged found!"
<< std::endl;
continue;
} else {
std::cout << "Storing alignment: " << childTree.second.get<std::string>("title") << std::endl;
}
vec_single_file_path.push_back(childTree.second.get<std::string>("file"));
vec_single_file_name.push_back(childTree.second.get<std::string>("file") + "/Zmumu.root");
vec_color.push_back(childTree.second.get<int>("color"));
Expand All @@ -271,6 +307,7 @@ int Zmumumerge(int argc, char* argv[]) {

//Fitting_GetMassmeanVSvariables(childTree.second.get<std::string>("file") + "/Zmumu.root", childTree.second.get<std::string>("file"));
}

TString merge_output = main_tree.get<std::string>("output");
//=============================================
vector<TString> vec_single_fittingoutput;
Expand Down Expand Up @@ -300,6 +337,7 @@ int Zmumumerge(int argc, char* argv[]) {
th1d_name_entries,
merge_output + Form("/entries_%s_GTs.pdf", tstring_variables_name[idx_variable].Data()));
}

//=============================================
return EXIT_SUCCESS;
}
Expand Down
16 changes: 12 additions & 4 deletions Alignment/OfflineValidation/src/FitWithRooFit.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ void FitWithRooFit::fit(
reinitializeParameters();

std::unique_ptr<RooDataHist> dh = importTH1(histo, xMin, xMax);
RooRealVar x(*static_cast<RooRealVar*>(dh->get()->find("x")));
RooRealVar x(*static_cast<RooRealVar*>(dh->get()->find("InvMass")));

// Make plot of binned dataset showing Poisson error bars (RooFit default)
RooPlot* frame = x.frame(RooFit::Title("di-muon mass M(#mu^{+}#mu^{-}) [GeV]"));
Expand All @@ -36,7 +36,8 @@ void FitWithRooFit::fit(
// -----------------------
// Fit with likelihood
if (!useChi2_) {
model->fitTo(*dh, RooFit::SumW2Error(sumW2Error));
// use RooFit::PrintLevel(-1) to reduce output to screen
model->fitTo(*dh, RooFit::SumW2Error(sumW2Error), RooFit::PrintLevel(-1));
}
// Fit with chi^2
else {
Expand All @@ -46,8 +47,12 @@ void FitWithRooFit::fit(
m.hesse();
// RooFitResult* r_chi2_wgt = m.save();
}

model->plotOn(frame, RooFit::LineColor(kRed));
model->plotOn(frame, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed));

// this causes segfaults when the fit causes a Warning in <ROOT::Math::Fitter::CalculateHessErrors>: Error when calculating Hessian
//model->plotOn(frame, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed));

model->paramOn(frame,
RooFit::Label("fit result"),
RooFit::Layout(0.65, 0.90, 0.90),
Expand All @@ -64,10 +69,13 @@ void FitWithRooFit::fit(
// If histogram has custom error (i.e. its contents is does not originate from a Poisson process
// but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead
// (same error bars as shown by ROOT)

RooPlot* frame2 = x.frame(RooFit::Title("di-muon mass M(#mu^{+}#mu^{-}) [GeV]"));
dh->plotOn(frame2, RooFit::DataError(RooAbsData::SumW2));
model->plotOn(frame2, RooFit::LineColor(kRed));
model->plotOn(frame2, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed));

// this causes segfaults when the fit causes a Warning in <ROOT::Math::Fitter::CalculateHessErrors>: Error when calculating Hessian
//model->plotOn(frame2, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed));
model->paramOn(frame2,
RooFit::Label("fit result"),
RooFit::Layout(0.65, 0.90, 0.90),
Expand Down

0 comments on commit df0dced

Please sign in to comment.