File indexing completed on 2024-06-13 03:23:11
0001 #include <cmath>
0002 #include <fstream>
0003 #include <iomanip>
0004 #include <iostream>
0005 #include <string>
0006
0007 #include "toolbox.h"
0008 #include "Options.h"
0009
0010 #include "TAxis.h"
0011 #include "TBranch.h"
0012 #include "TCanvas.h"
0013 #include "TChain.h"
0014 #include "TCut.h"
0015 #include "TF1.h"
0016 #include "TFile.h"
0017 #include "TFrame.h"
0018 #include "TH1.h"
0019 #include "TH1D.h"
0020 #include "TH1F.h"
0021 #include "TH2.h"
0022 #include "TH2F.h"
0023 #include "THStack.h"
0024 #include "TLatex.h"
0025 #include "TLegend.h"
0026 #include "TMath.h"
0027 #include "TPad.h"
0028 #include "TPaveText.h"
0029 #include "TRandom.h"
0030 #include "TString.h"
0031 #include "TStyle.h"
0032 #include "TSystem.h"
0033 #include "TTree.h"
0034
0035
0036 #include <boost/filesystem.hpp>
0037 #include <boost/property_tree/ptree.hpp>
0038 #include <boost/property_tree/json_parser.hpp>
0039 #include <boost/range/adaptor/indexed.hpp>
0040
0041 #include "Alignment/OfflineValidation/interface/FitWithRooFit.h"
0042 #include "Alignment/OfflineValidation/macros/CMS_lumi.h"
0043
0044 using namespace RooFit;
0045 using namespace std;
0046 using namespace AllInOneConfig;
0047 namespace pt = boost::property_tree;
0048
0049 void makeNicePlotStyle(RooPlot* plot) {
0050 plot->GetXaxis()->CenterTitle(true);
0051 plot->GetYaxis()->CenterTitle(true);
0052 plot->GetXaxis()->SetTitleFont(42);
0053 plot->GetYaxis()->SetTitleFont(42);
0054 plot->GetXaxis()->SetTitleSize(0.05);
0055 plot->GetYaxis()->SetTitleSize(0.05);
0056 plot->GetXaxis()->SetTitleOffset(0.9);
0057 plot->GetYaxis()->SetTitleOffset(1.3);
0058 plot->GetXaxis()->SetLabelFont(42);
0059 plot->GetYaxis()->SetLabelFont(42);
0060 plot->GetYaxis()->SetLabelSize(.05);
0061 plot->GetXaxis()->SetLabelSize(.05);
0062 }
0063
0064
0065 RooRealVar MuMu_mass("MuMu_mass", "MuMu_mass", 70, 110);
0066 static TString GT = "";
0067 TLatex* tlxg = new TLatex();
0068 class FitOut {
0069 public:
0070 double mean;
0071 double mean_err;
0072 double sigma;
0073 double sigma_err;
0074 double chi2;
0075 FitOut(double a, double b, double c, double d) : mean(a), mean_err(b), sigma(c), sigma_err(d) {}
0076 };
0077
0078 FitOut ZMassBinFit_OldTool(TH1D* th1d_input, TString s_name = "zmumu_fitting", TString output_path = "./") {
0079
0080 RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
0081
0082 double xMean = 91.1876;
0083 double xMin = th1d_input->GetXaxis()->GetXmin();
0084 double xMax = th1d_input->GetXaxis()->GetXmax();
0085
0086 double sigma(2.);
0087 double sigmaMin(0.1);
0088 double sigmaMax(10.);
0089
0090 double sigma2(0.1);
0091 double sigma2Min(0.);
0092 double sigma2Max(10.);
0093
0094 std::unique_ptr<FitWithRooFit> fitter = std::make_unique<FitWithRooFit>();
0095
0096 bool useChi2(false);
0097
0098 fitter->useChi2_ = useChi2;
0099 fitter->initMean(xMean, xMin, xMax);
0100 fitter->initSigma(sigma, sigmaMin, sigmaMax);
0101 fitter->initSigma2(sigma2, sigma2Min, sigma2Max);
0102 fitter->initAlpha(1.5, 0.05, 10.);
0103 fitter->initN(1, 0.01, 100.);
0104 fitter->initFGCB(0.4, 0., 1.);
0105 fitter->initGamma(2.4952, 0., 10.);
0106 fitter->gamma()->setConstant(kTRUE);
0107 fitter->initMean2(0., -20., 20.);
0108 fitter->mean2()->setConstant(kTRUE);
0109 fitter->initSigma(1.2, 0., 5.);
0110 fitter->initAlpha(1.5, 0.05, 10.);
0111 fitter->initN(1, 0.01, 100.);
0112 fitter->initExpCoeffA0(-1., -10., 10.);
0113 fitter->initExpCoeffA1(0., -10., 10.);
0114 fitter->initExpCoeffA2(0., -2., 2.);
0115 fitter->initFsig(0.9, 0., 1.);
0116 fitter->initA0(0., -10., 10.);
0117 fitter->initA1(0., -10., 10.);
0118 fitter->initA2(0., -10., 10.);
0119 fitter->initA3(0., -10., 10.);
0120 fitter->initA4(0., -10., 10.);
0121 fitter->initA5(0., -10., 10.);
0122 fitter->initA6(0., -10., 10.);
0123
0124 TCanvas* c1 = new TCanvas();
0125 c1->Clear();
0126 c1->SetLeftMargin(0.15);
0127 c1->SetRightMargin(0.10);
0128
0129 fitter->fit(th1d_input, "breitWignerTimesCB", "exponential", xMin, xMax, false);
0130
0131 c1->Print(Form("%s/fitResultPlot/%s_oldtool.pdf", output_path.Data(), s_name.Data()));
0132 c1->Print(Form("%s/fitResultPlot/%s_oldtool.root", output_path.Data(), s_name.Data()));
0133
0134 FitOut fitRes(
0135 fitter->mean()->getVal(), fitter->mean()->getError(), fitter->sigma()->getVal(), fitter->sigma()->getError());
0136
0137 return fitRes;
0138 }
0139 void Draw_th1d(TH1D* th1d_input, TString variable_name, TString output_path) {
0140 TCanvas* c = new TCanvas();
0141 c->cd();
0142 gStyle->SetOptStat(0);
0143 th1d_input->SetMarkerStyle(kFullCircle);
0144 th1d_input->SetMarkerColor(kRed);
0145 th1d_input->SetLineColor(kRed);
0146 th1d_input->SetMaximum(91.4);
0147 th1d_input->SetMinimum(90.85);
0148 th1d_input->GetXaxis()->SetTitle(variable_name.Data());
0149 th1d_input->GetXaxis()->SetTitleOffset(1.2);
0150 th1d_input->GetYaxis()->SetTitle("Mass mean (GeV)");
0151 th1d_input->Draw();
0152 tlxg->DrawLatexNDC(0.2, 0.8, Form("%s", GT.Data()));
0153 c->Print(Form("%s/fitResultPlot/mass_VS_%s.pdf", output_path.Data(), variable_name.Data()));
0154 }
0155
0156 const static int variables_number = 8;
0157 const TString tstring_variables_name[variables_number] = {
0158 "CosThetaCS", "DeltaEta", "EtaMinus", "EtaPlus", "PhiCS", "PhiMinus", "PhiPlus", "Pt"};
0159 const TString tstring_variables_name_label[variables_number] = {"cos #theta_{CS}",
0160 "#Delta #eta",
0161 "#eta_{#mu^{-}}",
0162 "#eta_{#mu^{+}}",
0163 "#phi_{CS}",
0164 "#phi_{#mu^{-}}",
0165 "#phi_{#mu^{+}}",
0166 "p_{T}"};
0167
0168 void Fitting_GetMassmeanVSvariables(TString inputfile_name, TString output_path) {
0169 TH2D* th2d_mass_variables[variables_number];
0170 TFile* inputfile = TFile::Open(inputfile_name.Data());
0171 TDirectoryFile* tdirectory = (TDirectoryFile*)inputfile->Get("DiMuonMassValidation");
0172 for (int i = 0; i < variables_number; i++) {
0173 TString th2d_name = Form("th2d_mass_%s", tstring_variables_name[i].Data());
0174 th2d_mass_variables[i] = (TH2D*)tdirectory->Get(th2d_name);
0175 }
0176
0177 gSystem->Exec(Form("mkdir -p %s", output_path.Data()));
0178 gSystem->Exec(Form("mkdir -p %s/fitResultPlot", output_path.Data()));
0179 TFile* outputfile = TFile::Open(Form("%s/fitting_output.root", output_path.Data()), "RECREATE");
0180 TH1D* th1d_variables_meanmass[variables_number];
0181 TH1D* th1d_variables_entries[variables_number];
0182 const int variables_rebin[variables_number] = {1, 1, 1, 1, 1, 1, 1, 1};
0183 const double xaxis_range[variables_number][2] = {
0184 {-1, 1}, {-4.8, 4.8}, {-2.4, 2.4}, {-2.4, 2.4}, {-1, 1}, {-M_PI, M_PI}, {-M_PI, M_PI}, {0, 100}};
0185 for (int i = 0; i < variables_number; i++) {
0186 TString th1d_name = Form("th1d_meanmass_%s", tstring_variables_name[i].Data());
0187
0188 th2d_mass_variables[i]->RebinY(variables_rebin[i]);
0189 th1d_variables_meanmass[i] = th2d_mass_variables[i]->ProjectionY(th1d_name, 1, 1, "d");
0190 for (int j = 0; j < th1d_variables_meanmass[i]->GetNbinsX(); j++) {
0191 if (i == 7 and j > 25) {
0192 continue;
0193 }
0194 std::cout << __PRETTY_FUNCTION__
0195 << " th1d_variables_meanmass[i]->GetNbinsX()=" << th1d_variables_meanmass[i]->GetNbinsX() << endl;
0196 std::cout << __PRETTY_FUNCTION__ << " th2d_mass_variables[i]->GetNbinsY()=" << th2d_mass_variables[i]->GetNbinsY()
0197 << endl;
0198 th1d_variables_meanmass[i]->SetBinContent(j, 0);
0199 th1d_variables_meanmass[i]->SetBinError(j, 0);
0200
0201 TString th1d_mass_temp_name = Form("th1d_mass_%s_%d", tstring_variables_name[i].Data(), j);
0202 TH1D* th1d_i = th2d_mass_variables[i]->ProjectionX(th1d_mass_temp_name, j, j, "d");
0203 th1d_i->Write(th1d_mass_temp_name);
0204 TString s_cut = Form("nocut");
0205 TString s_name = Form("%s_%d", tstring_variables_name[i].Data(), j);
0206
0207 FitOut fitR = ZMassBinFit_OldTool(th1d_i, s_name, output_path);
0208
0209 th1d_variables_meanmass[i]->SetBinContent(j, fitR.mean);
0210 th1d_variables_meanmass[i]->SetBinError(j, fitR.mean_err);
0211 }
0212
0213 th1d_variables_meanmass[i]->GetXaxis()->SetRangeUser(xaxis_range[i][0], xaxis_range[i][1]);
0214 Draw_th1d(th1d_variables_meanmass[i], tstring_variables_name[i], output_path);
0215 outputfile->cd();
0216 th1d_variables_meanmass[i]->Write(th1d_name);
0217
0218 TString th1d_name_entries = Form("th1d_entries_%s", tstring_variables_name[i].Data());
0219 th1d_variables_entries[i] = th2d_mass_variables[i]->ProjectionY(th1d_name_entries, 0, -1, "d");
0220 th1d_variables_entries[i]->GetXaxis()->SetTitle(tstring_variables_name[i].Data());
0221 th1d_variables_entries[i]->GetYaxis()->SetTitle("Entry");
0222 outputfile->cd();
0223 th1d_variables_entries[i]->Write(th1d_name_entries);
0224 }
0225
0226 if (outputfile->IsOpen()) {
0227
0228 const char* path = outputfile->GetPath();
0229
0230 if (path) {
0231 std::cout << "File is going to be written in the directory: " << path << " for input file: " << inputfile_name
0232 << std::endl;
0233 } else {
0234 std::cerr << "Error: Unable to determine the path." << std::endl;
0235 }
0236 outputfile->Close();
0237 delete outputfile;
0238 }
0239 }
0240
0241 const static int max_file_number = 10;
0242 void Draw_TH1D_forMultiRootFiles(const vector<TString>& file_names,
0243 const vector<TString>& label_names,
0244 const vector<int>& colors,
0245 const vector<int>& styles,
0246 const TString& Rlabel,
0247 const TString& th1d_name,
0248 const TString& xlabel,
0249 const TString& ylabel,
0250 const TString& output_name) {
0251 if (file_names.empty() || label_names.empty()) {
0252 std::cout << "Provided an empty list of file and label names" << std::endl;
0253 return;
0254 }
0255
0256
0257 assert(file_names.size() == label_names.size());
0258
0259 TH1D* th1d_input[max_file_number];
0260 TFile* file_input[max_file_number];
0261 for (auto const& filename : file_names | boost::adaptors::indexed(0)) {
0262 file_input[filename.index()] = TFile::Open(filename.value());
0263 th1d_input[filename.index()] = (TH1D*)file_input[filename.index()]->Get(th1d_name);
0264 th1d_input[filename.index()]->SetTitle("");
0265 }
0266
0267 int W = 800;
0268 int H = 800;
0269
0270 float T = 0.08 * H;
0271 float B = 0.12 * H;
0272 float L = 0.12 * W;
0273 float R = 0.04 * W;
0274
0275
0276 TString canvasName;
0277 canvasName.Form("canv_%s", th1d_name.Data());
0278
0279
0280 TCanvas* canv = new TCanvas(canvasName, canvasName, W, H);
0281 canv->SetFillColor(0);
0282 canv->SetBorderMode(0);
0283 canv->SetFrameFillStyle(0);
0284 canv->SetFrameBorderMode(0);
0285 canv->SetLeftMargin(L / W + 0.05);
0286 canv->SetRightMargin(R / W);
0287 canv->SetTopMargin(T / H);
0288 canv->SetBottomMargin(B / H);
0289 canv->SetTickx(0);
0290 canv->SetTicky(0);
0291 canv->SetGrid();
0292 canv->cd();
0293
0294 gStyle->SetOptStat(0);
0295
0296 TLegend* lg = new TLegend(0.3, 0.7, 0.7, 0.9);
0297 lg->SetFillStyle(0);
0298 lg->SetLineColor(0);
0299 lg->SetEntrySeparation(0.05);
0300
0301 double ymin = 0.;
0302 double ymax = 0.;
0303
0304 for (auto const& labelname : label_names | boost::adaptors::indexed(0)) {
0305 double temp_ymin = th1d_input[labelname.index()]->GetMinimum();
0306 double temp_ymax = th1d_input[labelname.index()]->GetMaximum();
0307 if (labelname.index() == 0) {
0308 ymin = temp_ymin;
0309 ymax = temp_ymax;
0310 }
0311 if (temp_ymin <= ymin) {
0312 ymin = temp_ymin;
0313 }
0314 if (temp_ymax >= ymax) {
0315 ymax = temp_ymax;
0316 }
0317 }
0318
0319 for (auto const& labelname : label_names | boost::adaptors::indexed(0)) {
0320 th1d_input[labelname.index()]->SetMarkerColor(colors[labelname.index()]);
0321 th1d_input[labelname.index()]->SetLineColor(colors[labelname.index()]);
0322 th1d_input[labelname.index()]->SetMarkerStyle(styles[labelname.index()]);
0323 th1d_input[labelname.index()]->GetXaxis()->SetTitle(xlabel);
0324 th1d_input[labelname.index()]->GetYaxis()->SetTitle(ylabel);
0325 th1d_input[labelname.index()]->GetYaxis()->SetTitleOffset(2.0);
0326 lg->AddEntry(th1d_input[labelname.index()], labelname.value());
0327
0328 TString label_meanmass_plot = "Mean M_{#mu#mu} (GeV)";
0329 if (ylabel.EqualTo(label_meanmass_plot)) {
0330 double ycenter = (ymax + ymin) / 2.0;
0331 double yrange = (ymax - ymin) * 2;
0332 double Ymin = ycenter - yrange;
0333 double Ymax = ycenter + yrange * 1.10;
0334 th1d_input[labelname.index()]->SetAxisRange(Ymin, Ymax, "Y");
0335 th1d_input[labelname.index()]->Draw("PEX0same");
0336 } else {
0337 double Ymin = ymin - ymin * 0.07;
0338 double Ymax = ymax + ymax * 0.35;
0339 th1d_input[labelname.index()]->SetAxisRange(Ymin, Ymax, "Y");
0340 th1d_input[labelname.index()]->Draw("HIST E1 same");
0341 }
0342 }
0343
0344 CMS_lumi(canv, 0, 3, Rlabel);
0345
0346 lg->Draw("same");
0347
0348 canv->Update();
0349 canv->RedrawAxis();
0350 canv->GetFrame()->Draw();
0351 canv->SaveAs(output_name);
0352
0353 if (output_name.Contains(".pdf")) {
0354 TString output_name_png(output_name);
0355 output_name_png.Replace(output_name_png.Index(".pdf"), 4, ".png");
0356 canv->SaveAs(output_name_png);
0357 }
0358 }
0359
0360 int Zmumumerge(int argc, char* argv[]) {
0361 vector<TString> vec_single_file_path;
0362 vector<TString> vec_single_file_name;
0363 vector<TString> vec_global_tag;
0364 vector<TString> vec_title;
0365 vector<int> vec_color;
0366 vector<int> vec_style;
0367 vector<TString> vec_right_title;
0368
0369 Options options;
0370 options.helper(argc, argv);
0371 options.parser(argc, argv);
0372 pt::ptree main_tree;
0373 pt::read_json(options.config, main_tree);
0374 pt::ptree alignments = main_tree.get_child("alignments");
0375 pt::ptree validation = main_tree.get_child("validation");
0376
0377
0378 std::vector<std::pair<std::string, pt::ptree>> alignmentsOrdered;
0379 for (const auto& childTree : alignments) {
0380 alignmentsOrdered.push_back(childTree);
0381 }
0382 std::sort(alignmentsOrdered.begin(),
0383 alignmentsOrdered.end(),
0384 [](const std::pair<std::string, pt::ptree>& left, const std::pair<std::string, pt::ptree>& right) {
0385 return left.second.get<int>("index") < right.second.get<int>("index");
0386 });
0387
0388 for (const auto& childTree : alignmentsOrdered) {
0389
0390 if (childTree.second.find("file") == childTree.second.not_found()) {
0391 std::cerr << "Ignoring alignment: " << childTree.second.get<std::string>("title") << ".\nNo file to merged found!"
0392 << std::endl;
0393 continue;
0394 } else {
0395 std::cout << "Storing alignment: " << childTree.second.get<std::string>("title") << std::endl;
0396 }
0397 vec_single_file_path.push_back(childTree.second.get<std::string>("file"));
0398 vec_single_file_name.push_back(childTree.second.get<std::string>("file") + "/Zmumu.root");
0399 vec_color.push_back(childTree.second.get<int>("color"));
0400 vec_style.push_back(childTree.second.get<int>("style"));
0401 if (childTree.second.find("customrighttitle") == childTree.second.not_found()) {
0402 vec_right_title.push_back("");
0403 } else {
0404 vec_right_title.push_back(childTree.second.get<std::string>("customrighttitle"));
0405 }
0406 vec_global_tag.push_back(childTree.second.get<std::string>("globaltag"));
0407 vec_title.push_back(childTree.second.get<std::string>("title"));
0408
0409
0410 }
0411
0412 TString merge_output = main_tree.get<std::string>("output");
0413
0414 vector<TString> vec_single_fittingoutput;
0415 vec_single_fittingoutput.clear();
0416 for (unsigned i = 0; i < vec_single_file_path.size(); i++) {
0417 Fitting_GetMassmeanVSvariables(vec_single_file_name[i], vec_single_file_path[i]);
0418 vec_single_fittingoutput.push_back(vec_single_file_path[i] + "/fitting_output.root");
0419 }
0420
0421 int files_number = vec_single_file_path.size();
0422 cout << "files_number=" << files_number << endl;
0423 for (int idx_variable = 0; idx_variable < variables_number; idx_variable++) {
0424 TString th1d_name = Form("th1d_meanmass_%s", tstring_variables_name[idx_variable].Data());
0425
0426 Draw_TH1D_forMultiRootFiles(
0427 vec_single_fittingoutput,
0428 vec_title,
0429 vec_color,
0430 vec_style,
0431 vec_right_title[0],
0432 th1d_name,
0433 tstring_variables_name_label[idx_variable].Data(),
0434 "Mean M_{#mu#mu} (GeV)",
0435 merge_output + Form("/meanmass_%s_GTs.pdf", tstring_variables_name[idx_variable].Data()));
0436
0437 TString th1d_name_entries = Form("th1d_entries_%s", tstring_variables_name[idx_variable].Data());
0438
0439 Draw_TH1D_forMultiRootFiles(
0440 vec_single_fittingoutput,
0441 vec_title,
0442 vec_color,
0443 vec_style,
0444 vec_right_title[0],
0445 th1d_name_entries,
0446 tstring_variables_name_label[idx_variable].Data(),
0447 "Events",
0448 merge_output + Form("/entries_%s_GTs.pdf", tstring_variables_name[idx_variable].Data()));
0449 }
0450
0451
0452 return EXIT_SUCCESS;
0453 }
0454 #ifndef DOXYGEN_SHOULD_SKIP_THIS
0455 int main(int argc, char* argv[]) { return Zmumumerge(argc, argv); }
0456 #endif
