Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:32:36

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 //#include "RooGlobalFunc.h"
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 
0043 using namespace RooFit;
0044 using namespace std;
0045 using namespace AllInOneConfig;
0046 namespace pt = boost::property_tree;
0047 /*--------------------------------------------------------------------*/
0048 void makeNicePlotStyle(RooPlot* plot) {
0049   plot->GetXaxis()->CenterTitle(true);
0050   plot->GetYaxis()->CenterTitle(true);
0051   plot->GetXaxis()->SetTitleFont(42);
0052   plot->GetYaxis()->SetTitleFont(42);
0053   plot->GetXaxis()->SetTitleSize(0.05);
0054   plot->GetYaxis()->SetTitleSize(0.05);
0055   plot->GetXaxis()->SetTitleOffset(0.9);
0056   plot->GetYaxis()->SetTitleOffset(1.3);
0057   plot->GetXaxis()->SetLabelFont(42);
0058   plot->GetYaxis()->SetLabelFont(42);
0059   plot->GetYaxis()->SetLabelSize(.05);
0060   plot->GetXaxis()->SetLabelSize(.05);
0061 }
0062 /*--------------------------------------------------------------------*/
0063 
0064 RooRealVar MuMu_mass("MuMu_mass", "MuMu_mass", 70, 110);
0065 static TString GT = "";
0066 TLatex* tlxg = new TLatex();
0067 class FitOut {
0068 public:
0069   double mean;
0070   double mean_err;
0071   double sigma;
0072   double sigma_err;
0073   double chi2;
0074   FitOut(double a, double b, double c, double d) : mean(a), mean_err(b), sigma(c), sigma_err(d) {}
0075 };
0076 
0077 FitOut ZMassBinFit_OldTool(TH1D* th1d_input, TString s_name = "zmumu_fitting", TString output_path = "./") {
0078   double xMin(75), xMax(105), xMean(91);
0079   double sigma = 2;
0080   double sigmaMin = 0.1;
0081   double sigmaMax = 10;
0082 
0083   FitWithRooFit* fitter = new FitWithRooFit();
0084   double sigma2(0.1), sigma2Min(0.), sigma2Max(10.), useChi2(false);
0085   fitter->useChi2_ = useChi2;
0086   fitter->initMean(xMean, xMin, xMax);
0087   fitter->initSigma(sigma, sigmaMin, sigmaMax);
0088   fitter->initSigma2(sigma2, sigma2Min, sigma2Max);
0089   fitter->initAlpha(1.5, 0.05, 10.);
0090   fitter->initN(1, 0.01, 100.);
0091   fitter->initFGCB(0.4, 0., 1.);
0092 
0093   fitter->initMean(91.1876, xMin, xMax);
0094   fitter->initGamma(2.4952, 0., 10.);
0095   fitter->gamma()->setConstant(kTRUE);
0096   fitter->initMean2(0., -20., 20.);
0097   fitter->mean2()->setConstant(kTRUE);
0098   fitter->initSigma(1.2, 0., 5.);
0099   fitter->initAlpha(1.5, 0.05, 10.);
0100   fitter->initN(1, 0.01, 100.);
0101   fitter->initExpCoeffA0(-1., -10., 10.);
0102   fitter->initExpCoeffA1(0., -10., 10.);
0103   fitter->initExpCoeffA2(0., -2., 2.);
0104   fitter->initFsig(0.9, 0., 1.);
0105   fitter->initA0(0., -10., 10.);
0106   fitter->initA1(0., -10., 10.);
0107   fitter->initA2(0., -10., 10.);
0108   fitter->initA3(0., -10., 10.);
0109   fitter->initA4(0., -10., 10.);
0110   fitter->initA5(0., -10., 10.);
0111   fitter->initA6(0., -10., 10.);
0112   TCanvas* c1 = new TCanvas();
0113   c1->Clear();
0114   c1->SetLeftMargin(0.15);
0115   c1->SetRightMargin(0.10);
0116 
0117   fitter->fit(th1d_input, "breitWignerTimesCB", "exponential", xMin, xMax, false);
0118 
0119   c1->Print(Form("%s/fitResultPlot/%s_oldtool.pdf", output_path.Data(), s_name.Data()));
0120   c1->Print(Form("%s/fitResultPlot/%s_oldtool.root", output_path.Data(), s_name.Data()));
0121 
0122   FitOut fitRes(
0123       fitter->mean()->getValV(), fitter->mean()->getError(), fitter->sigma()->getValV(), fitter->sigma()->getError());
0124   return fitRes;
0125 }
0126 void Draw_th1d(TH1D* th1d_input, TString variable_name) {
0127   TCanvas* c = new TCanvas();
0128   c->cd();
0129   gStyle->SetOptStat(0);
0130   th1d_input->SetMarkerStyle(kFullCircle);
0131   th1d_input->SetMarkerColor(kRed);
0132   th1d_input->SetLineColor(kRed);
0133   th1d_input->SetMaximum(91.4);
0134   th1d_input->SetMinimum(90.85);
0135   th1d_input->GetXaxis()->SetTitle(variable_name.Data());
0136   th1d_input->GetXaxis()->SetTitleOffset(1.2);
0137   th1d_input->GetYaxis()->SetTitle("Mass mean (GeV)");
0138   th1d_input->Draw();
0139   tlxg->DrawLatexNDC(0.2, 0.8, Form("%s", GT.Data()));
0140   c->Print(Form("%s/fitResultPlot/mass_VS_%s.pdf", GT.Data(), variable_name.Data()));
0141 }
0142 
0143 const static int variables_number = 8;
0144 const TString tstring_variables_name[variables_number] = {
0145     "CosThetaCS", "DeltaEta", "EtaMinus", "EtaPlus", "PhiCS", "PhiMinus", "PhiPlus", "Pt"};
0146 void Fitting_GetMassmeanVSvariables(TString inputfile_name, TString output_path) {
0147   TH2D* th2d_mass_variables[variables_number];
0148   TFile* inputfile = TFile::Open(inputfile_name.Data());
0149   TDirectoryFile* tdirectory = (TDirectoryFile*)inputfile->Get("DiMuonMassValidation");
0150   for (int i = 0; i < variables_number; i++) {
0151     TString th2d_name = Form("th2d_mass_%s", tstring_variables_name[i].Data());
0152     th2d_mass_variables[i] = (TH2D*)tdirectory->Get(th2d_name);
0153   }
0154 
0155   gSystem->Exec(Form("mkdir -p %s", output_path.Data()));
0156   gSystem->Exec(Form("mkdir -p %s/fitResultPlot", output_path.Data()));
0157   TFile* outpufile = TFile::Open(Form("%s/fitting_output.root", output_path.Data()), "recreate");
0158   TH1D* th1d_variables_meanmass[variables_number];
0159   TH1D* th1d_variables_entries[variables_number];
0160   const int variables_rebin[variables_number] = {1, 1, 1, 1, 1, 1, 1, 1};
0161   const double xaxis_range[variables_number][2] = {
0162       {-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}};
0163   for (int i = 0; i < variables_number; i++) {
0164     TString th1d_name = Form("th1d_meanmass_%s", tstring_variables_name[i].Data());
0165 
0166     th2d_mass_variables[i]->RebinY(variables_rebin[i]);
0167     th1d_variables_meanmass[i] = th2d_mass_variables[i]->ProjectionY(th1d_name, 1, 1, "d");
0168     for (int j = 0; j < th1d_variables_meanmass[i]->GetNbinsX(); j++) {
0169       if (i == 7 and j > 25) {
0170         continue;
0171       }
0172       cout << "th1d_variables_meanmass[i]->GetNbinsX()=" << th1d_variables_meanmass[i]->GetNbinsX() << endl;
0173       cout << "th2d_mass_variables[i]->GetNbinsY()=" << th2d_mass_variables[i]->GetNbinsY() << endl;
0174       th1d_variables_meanmass[i]->SetBinContent(j, 0);
0175       th1d_variables_meanmass[i]->SetBinError(j, 0);
0176 
0177       TString th1d_mass_temp_name = Form("th1d_mass_%s_%d", tstring_variables_name[i].Data(), j);
0178       TH1D* th1d_i = th2d_mass_variables[i]->ProjectionX(th1d_mass_temp_name, j, j, "d");
0179       th1d_i->Write(th1d_mass_temp_name);
0180       TString s_cut = Form("nocut");
0181       TString s_name = Form("%s_%d", tstring_variables_name[i].Data(), j);
0182 
0183       FitOut fitR = ZMassBinFit_OldTool(th1d_i, s_name, output_path);
0184       th1d_variables_meanmass[i]->SetBinContent(j, fitR.mean);
0185       th1d_variables_meanmass[i]->SetBinError(j, fitR.mean_err);
0186     }
0187     th1d_variables_meanmass[i]->GetXaxis()->SetRangeUser(xaxis_range[i][0], xaxis_range[i][1]);
0188     Draw_th1d(th1d_variables_meanmass[i], tstring_variables_name[i]);
0189     outpufile->cd();
0190     th1d_variables_meanmass[i]->Write(th1d_name);
0191 
0192     TString th1d_name_entries = Form("th1d_entries_%s", tstring_variables_name[i].Data());
0193     th1d_variables_entries[i] = th2d_mass_variables[i]->ProjectionY(th1d_name_entries, 0, -1, "d");
0194     th1d_variables_entries[i]->GetXaxis()->SetTitle(tstring_variables_name[i].Data());
0195     th1d_variables_entries[i]->GetYaxis()->SetTitle("Entry");
0196     outpufile->cd();
0197     th1d_variables_entries[i]->Write(th1d_name_entries);
0198   }
0199 
0200   outpufile->Write();
0201   outpufile->Close();
0202   delete outpufile;
0203 }
0204 
0205 const static int max_file_number = 10;
0206 void Draw_TH1D_forMultiRootFiles(const vector<TString>& file_names,
0207                                  const vector<TString>& label_names,
0208                                  const vector<int>& colors,
0209                                  const vector<int>& styles,
0210                                  const TString& th1d_name,
0211                                  const TString& output_name) {
0212   if (file_names.empty() || label_names.empty()) {
0213     cout << "Provided an empty list of file and label names" << std::endl;
0214     return;
0215   }
0216 
0217   // do not allow the list of files and labels names to differ
0218   assert(file_names.size() == label_names.size());
0219 
0220   TH1D* th1d_input[max_file_number];
0221   TFile* file_input[max_file_number];
0222   for (auto const& filename : file_names | boost::adaptors::indexed(0)) {
0223     file_input[filename.index()] = TFile::Open(filename.value());
0224     th1d_input[filename.index()] = (TH1D*)file_input[filename.index()]->Get(th1d_name);
0225     th1d_input[filename.index()]->SetTitle("");
0226   }
0227 
0228   TCanvas* c = new TCanvas();
0229   TLegend* lg = new TLegend(0.2, 0.7, 0.5, 0.95);
0230   c->cd();
0231   gStyle->SetOptStat(0);
0232 
0233   for (auto const& labelname : label_names | boost::adaptors::indexed(0)) {
0234     th1d_input[labelname.index()]->SetMarkerColor(colors[labelname.index()]);
0235     th1d_input[labelname.index()]->SetLineColor(colors[labelname.index()]);
0236     th1d_input[labelname.index()]->SetMarkerStyle(styles[labelname.index()]);
0237     th1d_input[labelname.index()]->Draw("same");
0238     lg->AddEntry(th1d_input[labelname.index()], labelname.value());
0239   }
0240   lg->Draw("same");
0241   c->SaveAs(output_name);
0242   if (output_name.Contains(".pdf")) {
0243     TString output_name_png(output_name);  // output_name is const, copy to modify
0244     output_name_png.Replace(output_name_png.Index(".pdf"), 4, ".png");
0245     c->SaveAs(output_name_png);
0246   }
0247 }
0248 
0249 int Zmumumerge(int argc, char* argv[]) {
0250   vector<TString> vec_single_file_path;
0251   vector<TString> vec_single_file_name;
0252   vector<TString> vec_global_tag;
0253   vector<TString> vec_title;
0254   vector<int> vec_color;
0255   vector<int> vec_style;
0256 
0257   Options options;
0258   options.helper(argc, argv);
0259   options.parser(argc, argv);
0260   pt::ptree main_tree;
0261   pt::read_json(options.config, main_tree);
0262   pt::ptree alignments = main_tree.get_child("alignments");
0263   pt::ptree validation = main_tree.get_child("validation");
0264   for (const auto& childTree : alignments) {
0265     vec_single_file_path.push_back(childTree.second.get<std::string>("file"));
0266     vec_single_file_name.push_back(childTree.second.get<std::string>("file") + "/Zmumu.root");
0267     vec_color.push_back(childTree.second.get<int>("color"));
0268     vec_style.push_back(childTree.second.get<int>("style"));
0269     vec_global_tag.push_back(childTree.second.get<std::string>("globaltag"));
0270     vec_title.push_back(childTree.second.get<std::string>("title"));
0271 
0272     //Fitting_GetMassmeanVSvariables(childTree.second.get<std::string>("file") + "/Zmumu.root", childTree.second.get<std::string>("file"));
0273   }
0274   TString merge_output = main_tree.get<std::string>("output");
0275   //=============================================
0276   vector<TString> vec_single_fittingoutput;
0277   vec_single_fittingoutput.clear();
0278   for (unsigned i = 0; i < vec_single_file_path.size(); i++) {
0279     Fitting_GetMassmeanVSvariables(vec_single_file_name[i], vec_single_file_path[i]);
0280     vec_single_fittingoutput.push_back(vec_single_file_path[i] + "/fitting_output.root");
0281   }
0282 
0283   int files_number = vec_single_file_path.size();
0284   cout << "files_number=" << files_number << endl;
0285   for (int idx_variable = 0; idx_variable < variables_number; idx_variable++) {
0286     TString th1d_name = Form("th1d_meanmass_%s", tstring_variables_name[idx_variable].Data());
0287     Draw_TH1D_forMultiRootFiles(
0288         vec_single_fittingoutput,
0289         vec_title,
0290         vec_color,
0291         vec_style,
0292         th1d_name,
0293         merge_output + Form("/meanmass_%s_GTs.pdf", tstring_variables_name[idx_variable].Data()));
0294     TString th1d_name_entries = Form("th1d_entries_%s", tstring_variables_name[idx_variable].Data());
0295     Draw_TH1D_forMultiRootFiles(
0296         vec_single_fittingoutput,
0297         vec_title,
0298         vec_color,
0299         vec_style,
0300         th1d_name_entries,
0301         merge_output + Form("/entries_%s_GTs.pdf", tstring_variables_name[idx_variable].Data()));
0302   }
0303   //=============================================
0304   return EXIT_SUCCESS;
0305 }
0306 #ifndef DOXYGEN_SHOULD_SKIP_THIS
0307 int main(int argc, char* argv[]) { return Zmumumerge(argc, argv); }
0308 #endif