analysis_INTT/plot/quickdraw_genhadron.C

0001 std::vector<float> centbin = {0, 3, 6, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70};
0002 int ncent = centbin.size() - 1;
0003
0004 int GetMbinNum(const std::string &fstr)
0005 {
0006     if (fstr.find("Centrality0to3") != std::string::npos)
0007         return 0;
0008     if (fstr.find("Centrality3to6") != std::string::npos)
0009         return 1;
0010     if (fstr.find("Centrality6to10") != std::string::npos)
0011         return 2;
0012     if (fstr.find("Centrality10to15") != std::string::npos)
0013         return 3;
0014     if (fstr.find("Centrality15to20") != std::string::npos)
0015         return 4;
0016     if (fstr.find("Centrality20to25") != std::string::npos)
0017         return 5;
0018     if (fstr.find("Centrality25to30") != std::string::npos)
0019         return 6;
0020     if (fstr.find("Centrality30to35") != std::string::npos)
0021         return 7;
0022     if (fstr.find("Centrality35to40") != std::string::npos)
0023         return 8;
0024     if (fstr.find("Centrality40to45") != std::string::npos)
0025         return 9;
0026     if (fstr.find("Centrality45to50") != std::string::npos)
0027         return 10;
0028     if (fstr.find("Centrality50to55") != std::string::npos)
0029         return 11;
0030     if (fstr.find("Centrality55to60") != std::string::npos)
0031         return 12;
0032     if (fstr.find("Centrality60to65") != std::string::npos)
0033         return 13;
0034     if (fstr.find("Centrality65to70") != std::string::npos)
0035         return 14;
0036     if (fstr.find("Centrality0to70") != std::string::npos)
0037         return 70;
0038
0039     return -1; // Default case if no match is found
0040 }
0041
0042 TH1D* GetCWHist(int Mbin)
0043 {
0044     TFile *f = new TFile(Form("/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run6_EvtZFitWidthChange/EvtVtxZ/FinalResult_10cm_Pol2BkgFit/completed/vtxZ_-10_10cm_MBin%d/Final_Mbin%d_00054280/Final_Mbin%d_00054280.root", Mbin, Mbin, Mbin), "READ");
0045     TH1D *h = (TH1D *)f->Get("h1D_dNdEta_truth");
0046     h->SetDirectory(0);
0047     f->Close();
0048     return h;
0049 }
0050
0051 void quickdraw_genhadron()
0052 {
0053     // directory to save the output plots; if the directory does not exist, create it using
0054     std::string plotdir = "quickdraw_genhadron";
0055     system(Form("mkdir -p %s", plotdir.c_str()));
0056
0057     ROOT::EnableImplicitMT();
0058     ROOT::RDataFrame df("minitree", "/sphenix/tg/tg01/hf/hjheng/ppg02/minitree/TrackletMinitree_Sim_HIJING_MDC2_ana467_20250225/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0/minitree_set2.root");
0059
0060     for (int i = 0; i < ncent; i++)
0061     {
0062         auto hM_genhadron = df.Define("GenHadron_eta_le4", "GenHadron_eta[abs(GenHadron_eta) < 4]").Filter(Form("is_min_bias&&(PV_z<=10&&PV_z>=-10)&&(MBD_centrality>%f && MBD_centrality<=%f)", centbin[i], centbin[i + 1])).Histo1D({"hM_genhadron", "hM_genhadron:Gen-hadron #eta:Entries", 27, -2.7, 2.7}, "GenHadron_eta_le4", "vtxzwei");
0063         auto h1WGevent = df.Filter(Form("is_min_bias&&(abs(PV_z)<=10.)&&(MBD_centrality>%f && MBD_centrality<=%f)", centbin[i], centbin[i + 1])).Histo1D({"h1WGevent", "h1WGevent:PV_z:Entries", 20, -10, 10}, "PV_z", "vtxzwei");
0064         float nWGevent = h1WGevent->Integral(0, h1WGevent->GetNbinsX() + 1);
0065
0066         cout << "Number of entries in hM_genhadron: " << hM_genhadron->GetEntries() << "; number of entries in h1WGevent: " << h1WGevent->Integral(0, h1WGevent->GetNbinsX() + 1) << endl;
0067
0068         hM_genhadron->Scale(1. / nWGevent, "width");
0069
0070         std::string centstr = Form("Centrality%dto%d", (int)centbin[i], (int)centbin[i + 1]);
0071         TH1D *hM_CW = GetCWHist(GetMbinNum(centstr));
0072
0073         // draw the histogram
0074         TCanvas *c = new TCanvas("c", "c", 800, 700);
0075         hM_genhadron->GetYaxis()->SetTitle("Entries");
0076         hM_genhadron->GetYaxis()->SetTitleOffset(1.4);
0077         hM_genhadron->GetYaxis()->SetRangeUser(hM_genhadron->GetMinimum(0) * 0.9, hM_genhadron->GetMaximum() * 1.2);
0078         hM_genhadron->SetMarkerSize(0);
0079         hM_genhadron->GetXaxis()->SetTitle("Gen-hadron #eta");
0080         hM_genhadron->Draw("histe");
0081         hM_CW->SetLineColor(kRed);
0082         hM_CW->SetMarkerSize(0);
0083         hM_CW->Draw("histe same");
0084         TLegend *l = new TLegend(1 - gPad->GetRightMargin() - 0.35, 1 - gPad->GetTopMargin() - 0.2, 1 - gPad->GetRightMargin() - 0.05, 1 - gPad->GetTopMargin() - 0.05);
0085         l->SetHeader(Form("Centrality %d-%d%%", (int)centbin[i], (int)centbin[i + 1]));
0086         l->SetTextAlign(kHAlignLeft + kVAlignTop);
0087         l->AddEntry(hM_genhadron.GetPtr(), "CMS approach", "l");
0088         l->AddEntry(hM_CW, "PHOBOS approach", "l");
0089         l->SetTextSize(0.035);
0090         l->SetBorderSize(0);
0091         l->SetFillStyle(0);
0092         l->Draw();
0093         c->SaveAs(Form("%s/hM_genhadron_centrality%dto%d.pdf", plotdir.c_str(), (int)centbin[i], (int)centbin[i + 1]));
0094         c->SaveAs(Form("%s/hM_genhadron_centrality%dto%d.png", plotdir.c_str(), (int)centbin[i], (int)centbin[i + 1]));
0095     }
0096
0097     system(Form("cd %s; gs -dNOPAUSE -sDEVICE=pdfwrite -sOUTPUTFILE=hM_genhadron_crosscheck_merged.pdf -dBATCH hM_genhadron_centrality*.pdf", plotdir.c_str()));
0098 }