analysis_INTT/plot/draw_dNdPt.C

0001 float TickSize = 0.03;
0002 float AxisTitleSize = 0.05;
0003 float AxisLabelSize = 0.045;
0004 float LeftMargin = 0.15;
0005 float RightMargin = 0.05;
0006 float TopMargin = 0.08;
0007 float BottomMargin = 0.13;
0008 float textscale = 2.6;
0009
0010 std::vector<TGraphAsymmErrors *> v_STAR;
0011
0012 void getHistogram(std::string HepdataFile = "./dNdPt/HEPData-ins1676541-v1-Figure_2___0.40__Mee__0.76_GeV_c2_in_Au+Au_collisions.root")
0013 {
0014     TFile *f = new TFile(HepdataFile.c_str(), "READ");
0015     TH1F *hm = (TH1F *)f->Get("Hist1D_y1");
0016     hm->SetDirectory(0);
0017     return hm;
0018 }
0019
0020 void draw_dNdPt()
0021 {
0022     // directory to save the output plots; if the directory does not exist, create it using
0023     std::string plotdir = "./dNdPt";
0024     system(Form("mkdir -p %s", plotdir.c_str()));
0025
0026     std::string filepath = "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_HIJING_MDC2_ana472_20250307/";
0027      ROOT::EnableImplicitMT();
0028     std::vector<std::string> filelist;
0029     for (int i = 0; i < 100; i++)
0030     {
0031         std::string fname = filepath + "ntuple_" + std::string(TString::Format("%05d", i).Data()) + ".root";
0032
0033         // Check if file exists and is valid
0034         TFile *f = TFile::Open(fname.c_str());
0035         if (!f || f->IsZombie() || f->GetSize() <= 0)
0036         {
0037             std::cout << "Skipping invalid or empty file: " << fname << std::endl;
0038             if (f)
0039                 f->Close();
0040             continue;
0041         }
0042
0043         // Check if file contains the required tree
0044         if (!f->GetListOfKeys()->Contains("EventTree"))
0045         {
0046             std::cout << "Skipping file without EventTree: " << fname << std::endl;
0047             f->Close();
0048             continue;
0049         }
0050
0051         f->Close();
0052         std::cout << "Adding file: " << fname << std::endl;
0053         filelist.push_back(fname);
0054     }
0055     ROOT::RDataFrame df("EventTree", filelist);
0056
0057     // first, filter the events with is_min_bias==1 and MBD_centrality > 60 and MBD_centrality < 80
0058     auto df_minbias = df.Filter("is_min_bias==1 && MBD_centrality >= 60 && MBD_centrality < 80");
0059     int Nevt = df_minbias.Count().GetValue();
0060     // Add a new column, PrimaryG4P_Pt_eta1p1, with PrimaryG4P_isChargeHadron==1 within |eta|<=1.1 for each event
0061     auto df_Nptcl = df_minbias.Define("PrimaryG4P_Pt_eta1p1",
0062                                       [](const std::vector<float> &PrimaryG4P_Pt, const std::vector<float> &PrimaryG4P_Eta, const std::vector<bool> &PrimaryG4P_isChargeHadron)
0063                                       {
0064                                           std::vector<float> Pt_eta1p1;
0065                                           for (size_t i = 0; i < PrimaryG4P_Pt.size(); i++)
0066                                           {
0067                                               if (PrimaryG4P_isChargeHadron[i] == 1 && fabs(PrimaryG4P_Eta[i]) <= 1.1)
0068                                               {
0069                                                   Pt_eta1p1.push_back(PrimaryG4P_Pt[i]);
0070                                               }
0071                                           }
0072                                           return Pt_eta1p1;
0073                                       },
0074                                       {"PrimaryG4P_Pt", "PrimaryG4P_Eta", "PrimaryG4P_isChargeHadron"});
0075
0076     // histogra of PrimaryG4P_Pt with PrimaryG4P_isChargeHadron==1 within |eta|<=1.1
0077     auto hM_dNchdPt = df_Nptcl.Histo1D({"hM_dNchdPt", "dNch/dPt;Pt (GeV/c);dNch/dPt", 50, 0, 1}, "PrimaryG4P_Pt_eta1p1");
0078
0079     // normalize the histogram
0080     hM_dNchdPt->Scale(1./Nevt, "width");
0081
0082     // get the histogram from Hepdata
0083     TH1F *h_STAR_Mee0p4to0p76 = getHistogram(Form("%s/HEPData-ins1676541-v1-Figure_2___0.40__Mee__0.76_GeV_c2_in_Au+Au_collisions.root",plotdir.c_str()));
0084
0085
0086     // draw the histogram
0087     TCanvas *c = new TCanvas("c", "c", 800, 600);
0088     c->SetLogy(1);
0089     c->cd();
0090     hM_dNchdPt->SetLineColor(kBlack);
0091     hM_dNchdPt->SetLineWidth(2);
0092     hM_dNchdPt->Draw("HIST");
0093     c->SaveAs(Form("%s/dNdPt.png", plotdir.c_str()));
0094     c->SaveAs(Form("%s/dNdPt.pdf", plotdir.c_str()));
0095 }