LightFlavorRatios/lfratio_pythia/compare_lfratio_pythia.C

0001 #include "TCanvas.h"
0002 #include "TFile.h"
0003 #include "TGraph.h"
0004 #include "TH1D.h"
0005 #include "TLegend.h"
0006 #include "TStyle.h"
0007 #include "TSystem.h"
0008
0009 #include <algorithm>
0010 #include <iostream>
0011
0012 void compare_lfratio_pythia()
0013 {
0014     const char *output_file_pattern = "/sphenix/tg/tg01/hf/hjheng/HF-analysis/simulation/Pythia_ppMinBias/lfratio/lfratio_hist_%d.root";
0015     const int n_files = 5000;
0016
0017     const int n_regions = 3;
0018     const char *region_labels[n_regions] = {"|y| #leq 0.5", "0.5 < |y| #leq 1.0", "|y| > 1.0"};
0019     const char *kshort_names[n_regions] = {"h_prompt_kshort_pt_y_le_0p5", "h_prompt_kshort_pt_y_0p5_1p0", "h_prompt_kshort_pt_y_gt_1p0"};
0020     const char *lambda_names[n_regions] = {"h_prompt_lambda_pt_y_le_0p5", "h_prompt_lambda_pt_y_0p5_1p0", "h_prompt_lambda_pt_y_gt_1p0"};
0021
0022     TH1D *sum_kshort[n_regions] = {nullptr, nullptr, nullptr};
0023     TH1D *sum_lambda[n_regions] = {nullptr, nullptr, nullptr};
0024
0025     int n_files_used = 0;
0026     for (int ifile = 0; ifile < n_files; ++ifile)
0027     {
0028         const TString file_name = Form(output_file_pattern, ifile);
0029         if (gSystem->AccessPathName(file_name))
0030             continue;
0031
0032         TFile *input = TFile::Open(file_name, "READ");
0033         if (!input || input->IsZombie())
0034         {
0035             delete input;
0036             continue;
0037         }
0038
0039         for (int iregion = 0; iregion < n_regions; ++iregion)
0040         {
0041             TH1D *kshort = static_cast<TH1D *>(input->Get(kshort_names[iregion]));
0042             TH1D *lambda = static_cast<TH1D *>(input->Get(lambda_names[iregion]));
0043             if (!kshort || !lambda)
0044                 continue;
0045
0046             if (!sum_kshort[iregion])
0047             {
0048                 sum_kshort[iregion] = static_cast<TH1D *>(kshort->Clone(Form("sum_%s", kshort_names[iregion])));
0049                 sum_lambda[iregion] = static_cast<TH1D *>(lambda->Clone(Form("sum_%s", lambda_names[iregion])));
0050                 sum_kshort[iregion]->SetDirectory(nullptr);
0051                 sum_lambda[iregion]->SetDirectory(nullptr);
0052                 sum_kshort[iregion]->Sumw2();
0053                 sum_lambda[iregion]->Sumw2();
0054             }
0055             else
0056             {
0057                 sum_kshort[iregion]->Add(kshort);
0058                 sum_lambda[iregion]->Add(lambda);
0059             }
0060         }
0061
0062         ++n_files_used;
0063         input->Close();
0064         delete input;
0065     }
0066
0067     TH1D *ratio[n_regions] = {nullptr, nullptr, nullptr};
0068     const int colors[n_regions] = {kBlack, kRed + 1, kBlue + 1};
0069     double ymax = 0.0;
0070
0071     for (int iregion = 0; iregion < n_regions; ++iregion)
0072     {
0073         if (!sum_kshort[iregion] || !sum_lambda[iregion])
0074         {
0075             std::cerr << "No summed histogram for " << region_labels[iregion] << std::endl;
0076             return;
0077         }
0078
0079         ratio[iregion] = static_cast<TH1D *>(sum_lambda[iregion]->Clone(Form("h_prompt_lfratio_yregion%d", iregion)));
0080         ratio[iregion]->SetDirectory(nullptr);
0081         ratio[iregion]->SetTitle(";p_{T} [GeV];(#Lambda+#bar{#Lambda})/(2K^{0}_{S})");
0082         ratio[iregion]->Divide(sum_lambda[iregion], sum_kshort[iregion], 1.0, 2.0);
0083         ratio[iregion]->SetLineColor(colors[iregion]);
0084         ratio[iregion]->SetLineWidth(2);
0085         ymax = std::max(ymax, ratio[iregion]->GetMaximum());
0086     }
0087
0088     TGraph *ratio_graph[n_regions] = {nullptr, nullptr, nullptr};
0089     for (int iregion = 0; iregion < n_regions; ++iregion)
0090     {
0091         ratio_graph[iregion] = new TGraph();
0092         ratio_graph[iregion]->SetName(Form("g_prompt_lfratio_yregion%d", iregion));
0093         ratio_graph[iregion]->SetLineColor(colors[iregion]);
0094         ratio_graph[iregion]->SetLineWidth(2);
0095
0096         for (int ibin = 1; ibin <= ratio[iregion]->GetNbinsX(); ++ibin)
0097         {
0098             if (sum_kshort[iregion]->GetBinContent(ibin) <= 0.0)
0099                 continue;
0100
0101             const int point = ratio_graph[iregion]->GetN();
0102             ratio_graph[iregion]->SetPoint(point, ratio[iregion]->GetBinCenter(ibin), ratio[iregion]->GetBinContent(ibin));
0103         }
0104     }
0105
0106     gStyle->SetOptStat(0);
0107
0108     TCanvas *canvas = new TCanvas("c_prompt_lfratio_y_comparison", "Prompt ratio by rapidity", 800, 700);
0109     canvas->SetLeftMargin(0.15);
0110     // canvas->SetBottomMargin(0.12);
0111
0112     ratio[0]->SetMinimum(0.0);
0113     // ratio[0]->SetMaximum(ymax > 0.0 ? 1.25 * ymax : 1.0);
0114     ratio[0]->SetMaximum(1.0);
0115     ratio[0]->Draw("AXIS");
0116     for (int iregion = 0; iregion < n_regions; ++iregion)
0117         ratio_graph[iregion]->Draw("L SAME");
0118
0119     TLegend *legend = new TLegend(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.25, gPad->GetLeftMargin() + 0.3, 1 - gPad->GetTopMargin() - 0.05);
0120     legend->SetHeader("Pythia 8 (Detroit tune, pp #sqrt{s} = 200 GeV)");
0121     legend->SetTextSize(0.04);
0122     legend->SetBorderSize(0);
0123     legend->SetFillStyle(0);
0124     for (int iregion = 0; iregion < n_regions; ++iregion)
0125         legend->AddEntry(ratio_graph[iregion], region_labels[iregion], "l");
0126     legend->Draw();
0127
0128     canvas->SaveAs("compare_lfratio_pythia.pdf");
0129     canvas->SaveAs("compare_lfratio_pythia.png");
0130
0131     std::cout << "Read " << n_files_used << " files" << std::endl;
0132 }