analysis_INTT/plot/quickdraw_centralitysim.C

0001 #include "./plotutil.h"
0002 float ytitleoffset = 1.6;
0003
0004 void quickdraw_centralitysim()
0005 {
0006     std::string plotdir = "./quickcheck_centralitysim";
0007     system(Form("mkdir -p %s", plotdir.c_str()));
0008
0009     // Set a timer to measure the time it takes to run the analysis
0010     TStopwatch timer;
0011     timer.Start();
0012
0013     ROOT::EnableImplicitMT();
0014     std::vector<std::string> filelist;
0015     for (int i = 0; i < 5000; i++)
0016     {
0017         std::string fname = "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_AMPT_MDC2_ana472_20250310/ntuple_" + std::string(TString::Format("%05d", i).Data()) + ".root";
0018
0019         // Check if file exists and is valid
0020         TFile *f = TFile::Open(fname.c_str());
0021         if (!f || f->IsZombie() || f->GetSize() <= 0) {
0022             std::cout << "Skipping invalid or empty file: " << fname << std::endl;
0023             if (f) f->Close();
0024             continue;
0025         }
0026
0027         // Check if file contains the required tree
0028         if (!f->GetListOfKeys()->Contains("EventTree")) {
0029             std::cout << "Skipping file without EventTree: " << fname << std::endl;
0030             f->Close();
0031             continue;
0032         }
0033
0034         f->Close();
0035         std::cout << "Adding file: " << fname << std::endl;
0036         filelist.push_back(fname);
0037     }
0038     ROOT::RDataFrame df("EventTree", filelist);
0039
0040     auto count_all = df.Count();
0041     std::cout << "Number of entries in the dataframe: " << *count_all << std::endl;
0042
0043     // histogram for MBD centrality without cut
0044     auto hM_MBDcentrality_wocut = df.Histo1D({"hM_MBDcentrality_wocut", "hM_MBDcentrality_wocut", 102, -1.5, 100.5}, "MBD_centrality");
0045     auto hM_MBDcentrality_isMinBias = df.Filter("is_min_bias").Histo1D({"hM_MBDcentrality_isMinBias", "hM_MBDcentrality_isMinBias", 102, -1.5, 100.5}, "MBD_centrality");
0046     auto hM_Npart_isMinBias = df.Filter("is_min_bias").Histo1D({"hM_Npart_isMinBias", "hM_Npart_isMinBias", 401, -0.5, 400.5}, "npart");
0047     std::vector<float> centralitybins = {0, 3, 6, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100};
0048     // histogram for npart in all centrality intervals, for is_min_bias events
0049     std::vector<ROOT::RDF::RResultPtr<TH1D>> hM_Npart_centralities;
0050     for (size_t i = 0; i < centralitybins.size() - 1; i++)
0051     {
0052         auto hM = df.Filter(Form("is_min_bias && MBD_centrality >= %f && MBD_centrality < %f", centralitybins[i], centralitybins[i + 1])).Histo1D({Form("hM_Npart_centralities_%zu", i), Form("hM_Npart_centralities_%zu", i), 401, -0.5, 400.5}, "npart");
0053         hM_Npart_centralities.push_back(hM);
0054     }
0055
0056     std::vector<float> meanNpart_centralities;
0057     // Draw the histograms
0058     TCanvas *c = new TCanvas("canvas", "canvas", 800, 600);
0059     c->cd();
0060     gPad->SetTopMargin(0.25);
0061     hM_MBDcentrality_wocut->GetXaxis()->SetTitle("Centrality percentile");
0062     hM_MBDcentrality_wocut->GetYaxis()->SetTitle("Counts");
0063     hM_MBDcentrality_wocut->GetYaxis()->SetTitleOffset(ytitleoffset);
0064     hM_MBDcentrality_wocut->GetYaxis()->SetRangeUser(0, 1.1 * hM_MBDcentrality_wocut->GetMaximum());
0065     hM_MBDcentrality_wocut->SetLineColor(kBlack);
0066     hM_MBDcentrality_wocut->SetLineWidth(2);
0067     hM_MBDcentrality_isMinBias->SetLineColor(kTBriBlue);
0068     hM_MBDcentrality_isMinBias->SetLineWidth(2);
0069     hM_MBDcentrality_wocut->Draw("hist");
0070     hM_MBDcentrality_isMinBias->Draw("hist same");
0071     c->RedrawAxis();
0072     TLegend *l = new TLegend(gPad->GetLeftMargin(), 1 - gPad->GetTopMargin() + 0.02, 1 - gPad->GetRightMargin() - 0.45, 0.98);
0073     l->SetNColumns(1);
0074     l->AddEntry(hM_MBDcentrality_wocut.GetPtr(), "Without cut", "l");
0075     l->AddEntry(hM_MBDcentrality_isMinBias.GetPtr(), "Is min bias", "l");
0076     l->SetTextSize(0.035);
0077     l->SetBorderSize(0);
0078     l->SetFillStyle(0);
0079     l->Draw();
0080     c->SaveAs(Form("%s/MBDcentrality_selections.pdf", plotdir.c_str()));
0081     c->SaveAs(Form("%s/MBDcentrality_selections.png", plotdir.c_str()));
0082
0083     c->Clear();
0084     gPad->SetTopMargin(0.35);
0085     gPad->Modified();
0086     gPad->Update();
0087     c->SetLogy();
0088     hM_Npart_isMinBias->GetXaxis()->SetTitle("Number of participants, N_{part}");
0089     hM_Npart_isMinBias->GetYaxis()->SetTitle("Counts");
0090     hM_Npart_isMinBias->GetYaxis()->SetTitleOffset(ytitleoffset);
0091     hM_Npart_isMinBias->SetLineColor(kBlack);
0092     hM_Npart_isMinBias->SetLineWidth(2);
0093     hM_Npart_isMinBias->Draw("hist");
0094     for (size_t i = 0; i < hM_Npart_centralities.size(); i++)
0095     {
0096         hM_Npart_centralities[i]->SetLineWidth(2);
0097         hM_Npart_centralities[i]->Draw("plc same");
0098     }
0099     hM_Npart_isMinBias->Draw("hist same");
0100     c->RedrawAxis();
0101     TLegend *leg = new TLegend(gPad->GetLeftMargin(), 1 - gPad->GetTopMargin() + 0.01, 1 - gPad->GetRightMargin(), 0.98);
0102     leg->SetNColumns(2);
0103     leg->SetBorderSize(0);
0104     leg->SetFillStyle(0);
0105     leg->SetTextSize(0.02);
0106     for (size_t i = 0; i < hM_Npart_centralities.size(); i++)
0107     {
0108         // get the mean of Npart from the histogram
0109         float meanNpart = hM_Npart_centralities[i]->GetMean();
0110         meanNpart_centralities.push_back(meanNpart);
0111         leg->AddEntry(hM_Npart_centralities[i].GetPtr(), Form("Centrality [%d-%d]%%, <N_{part}>=%.1f", (int)centralitybins[i], (int)centralitybins[i + 1], meanNpart), "l");
0112     }
0113     leg->Draw();
0114     c->SaveAs(Form("%s/Npart_centralities.pdf", plotdir.c_str()));
0115     c->SaveAs(Form("%s/Npart_centralities.png", plotdir.c_str()));
0116
0117     // print out the centrality bins and the mean Npart in each bin in Latex table format
0118     std::cout << "Latex code for the table:\n";
0119     std::cout << "---------------------------------------------------------------------------\n";
0120     std::cout << "\\begin{table}[ht]\n";
0121     std::cout << "\\centering\n";
0122     std::cout << "\\begin{tabular}{c|c}\n";
0123     std::cout << "\\hline\n";
0124     std::cout << "Centrality interval [\\%] & $\\langle$\\npart$\\rangle$ \\\\\n";
0125     std::cout << "\\hline\n";
0126
0127     for (size_t i = 0; i < centralitybins.size() - 1; ++i)
0128     {
0129         if (centralitybins[i] >= 70)
0130             continue;
0131         std::cout << std::fixed << std::setprecision(0) << centralitybins[i] << "-" << centralitybins[i + 1] << " & ";
0132         print_with_significant_digits(meanNpart_centralities[i]);
0133         std::cout << " \\\\\n";
0134     }
0135
0136     std::cout << "\\hline\n";
0137     std::cout << "\\end{tabular}\n";
0138     std::cout << "\\caption{Centrality intervals and corresponding $\\langle$\\npart$\\rangle$ values.}\n";
0139     std::cout << "\\label{tab:centinterval_npart}\n";
0140     std::cout << "\\end{table}\n";
0141     std::cout << "---------------------------------------------------------------------------\n";
0142
0143     // simple print out of the centrality bins and the mean Npart in each bin
0144     std::cout << "Centrality interval [%]: " << std::endl;
0145     for (size_t i = 0; i < centralitybins.size() - 1; ++i)
0146     {
0147         if (centralitybins[i] >= 70)
0148             continue;
0149         std::cout << std::fixed << std::setprecision(0) << centralitybins[i] << "-" << centralitybins[i + 1] << ", ";
0150     }
0151     std::cout << std::endl;
0152     std::cout << "Mean Npart: " << std::endl;
0153     for (size_t i = 0; i < meanNpart_centralities.size(); ++i)
0154     {
0155         if (centralitybins[i] >= 70)
0156             continue;
0157         print_with_significant_digits(meanNpart_centralities[i], 6);
0158         std::cout << ", ";
0159     }
0160     std::cout << std::endl;
0161
0162     timer.Stop();
0163     timer.Print();
0164 }