analysis_INTT/plot/generator_dNdEtaNpart2.C

0001 #include <iostream>
0002 #include <string>
0003 #include <regex>
0004
0005 void generator_dNdEtaNpart2(std::string filepath = "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_HIJING_MDC2_ana472_20250307/")
0006 {
0007     std::string plotdir = "./generator_dNdEtaNpart2";
0008     system(Form("mkdir -p %s", plotdir.c_str()));
0009
0010     std::regex pattern(R"(/Sim_([^/_]+))"); // Match "Sim_" followed by the generator name
0011     std::smatch match;
0012
0013     if (!std::regex_search(filepath, match, pattern))
0014     {
0015         std::cerr << "Generator name not found in the file path" << std::endl;
0016         return;
0017     }
0018
0019     std::string generator = match[1].str();
0020
0021     std::cout << "Generator: " << generator << std::endl;
0022
0023     // Set a timer to measure the time it takes to run the analysis
0024     TStopwatch timer;
0025     timer.Start();
0026
0027     ROOT::EnableImplicitMT();
0028     std::vector<std::string> filelist;
0029     for (int i = 0; i < 5000; 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
0058     auto df_minbias = df.Filter("is_min_bias==1");
0059     // get number of minimum bias events
0060     int Nminbias = df_minbias.Count().GetValue();
0061
0062     // Add a column to calculate the average number of PrimaryG4P_Eta with PrimaryG4P_isChargeHadron==1 within |eta|<=0.3
0063     auto df_truthdNdEta = df_minbias
0064                               .Define("truthdNdEta_etaleq0p3",
0065                                       [](const std::vector<float> &PrimaryG4P_Eta, const std::vector<bool> &PrimaryG4P_isChargeHadron)
0066                                       {
0067                                           double N = 0;
0068                                           for (size_t i = 0; i < PrimaryG4P_Eta.size(); i++)
0069                                           {
0070                                               if (PrimaryG4P_isChargeHadron[i] == true && fabs(PrimaryG4P_Eta[i]) <= 0.3)
0071                                               {
0072                                                   N++;
0073                                               }
0074                                           }
0075
0076                                           return N / 0.6;
0077                                       },
0078                                       {"PrimaryG4P_Eta", "PrimaryG4P_isChargeHadron"})
0079                               .Define("truthdNdEta_etaleq0p3overNpart2", [](double N, int Npart) { return static_cast<double>(static_cast<double>(N) / (0.5 * static_cast<double>(Npart))); }, {"truthdNdEta_etaleq0p3", "npart"});
0080
0081     // 2D histogram for truth (dN/dEta)/(0.5*npart) vs npart/2
0082     auto h_truthdNdEta = df_truthdNdEta.Histo2D<int, double>({"h_truthdNdEta", "h_truthdNdEta", 40, 0, 400, 30, 0, 6}, "npart", "truthdNdEta_etaleq0p3overNpart2");
0083
0084     TCanvas *c = new TCanvas("c", "c", 800, 700);
0085     gPad->SetRightMargin(0.22);
0086     gPad->SetTopMargin(0.07);
0087     gPad->SetLeftMargin(0.17);
0088     c->cd();
0089     h_truthdNdEta->GetXaxis()->SetTitle("N_{part}");
0090     h_truthdNdEta->GetYaxis()->SetTitle("(dN_{ch}/d#eta)/(0.5N_{part})|_{|#eta|#leq0.3}");
0091     h_truthdNdEta->GetZaxis()->SetTitle("Entries");
0092     h_truthdNdEta->GetZaxis()->SetTitleOffset(1.7);
0093     h_truthdNdEta->Draw("colz");
0094     // profile
0095     TProfile *p_truthdNdEta = h_truthdNdEta->ProfileX();
0096     p_truthdNdEta->SetMarkerSize(0);
0097     p_truthdNdEta->SetLineColor(kRed);
0098     p_truthdNdEta->DrawCopy("hist l same");
0099     // draw the text for generator
0100     TLatex *t = new TLatex();
0101     t->SetTextSize(0.04);
0102     t->SetTextAlign(kHAlignRight + kVAlignCenter);
0103     t->DrawLatexNDC(0.75, 0.2, generator.c_str());
0104     c->SaveAs(Form("%s/h_truthdNdEta_%s.png", plotdir.c_str(), generator.c_str()));
0105     c->SaveAs(Form("%s/h_truthdNdEta_%s.pdf", plotdir.c_str(), generator.c_str()));
0106
0107     // Set up the TGraph for the profile
0108     TGraph *g_truthdNdEta = new TGraph(p_truthdNdEta->GetNbinsX());
0109     for (int i = 0; i < p_truthdNdEta->GetNbinsX(); i++)
0110     {
0111         g_truthdNdEta->SetPoint(i, p_truthdNdEta->GetBinCenter(i + 1), p_truthdNdEta->GetBinContent(i + 1));
0112     }
0113     g_truthdNdEta->SetName("g_truthdNdEta");
0114     g_truthdNdEta->SetTitle("g_truthdNdEta");
0115
0116     // Save to file
0117     TFile *fout = new TFile(Form("%s/truthdNdEta_%s.root", plotdir.c_str(), generator.c_str()), "RECREATE");
0118     h_truthdNdEta->Write();
0119     p_truthdNdEta->Write();
0120     g_truthdNdEta->Write();
0121     fout->Close();
0122 }