analysis_INTT/plot/countLepton.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 void countLepton()
0011 {
0012     // directory to save the output plots; if the directory does not exist, create it using
0013     std::string plotdir = "./countLepton";
0014     system(Form("mkdir -p %s", plotdir.c_str()));
0015
0016     ROOT::EnableImplicitMT();
0017     ROOT::RDataFrame df("minitree", "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/test.root");
0018
0019     // define 2 new columns, one for the number of leptons and one for the number of total charged particles
0020     // -> the former counts the number of elements in the GenHadron_PID vector with a PID of 11, -11, 13, or -13
0021     // -> the latter counts the total number of elements in the GenHadron_PID vector
0022     auto df_with_lepton_count = df.Define("Nleptons",
0023                                           [](const std::vector<int> &GenHadron_PID, const std::vector<float> &GenHadron_eta)
0024                                           {
0025                                               int Nleptons = 0;
0026                                               for (size_t i = 0; i < GenHadron_PID.size(); i++)
0027                                               {
0028                                                   if ((abs(GenHadron_PID[i]) == 11 || abs(GenHadron_PID[i]) == 13) && fabs(GenHadron_eta[i]) < 1.5)
0029                                                   {
0030                                                       Nleptons++;
0031                                                   }
0032                                               }
0033                                               return Nleptons;
0034                                           },
0035                                           {"GenHadron_PID", "GenHadron_eta"})
0036                                     .Define("Ncharged", [](const std::vector<int> &GenHadron_PID) { return static_cast<int>(GenHadron_PID.size()); }, {"GenHadron_PID"})
0037                                     .Define("NStrange",
0038                                             [](const std::vector<int> &PrimaryG4P_PID, const std::vector<float> &PrimaryG4P_eta)
0039                                             {
0040                                                 int NStrange = 0;
0041                                                 for (size_t i = 0; i < PrimaryG4P_PID.size(); i++)
0042                                                 {
0043                                                     if ((PrimaryG4P_PID[i] == 310 || abs(PrimaryG4P_PID[i]) == 3122) && fabs(PrimaryG4P_eta[i]) < 1.5)
0044                                                     {
0045                                                         NStrange++;
0046                                                     }
0047                                                 }
0048                                                 return NStrange;
0049                                             },
0050                                             {"PrimaryG4P_PID", "PrimaryG4P_eta"});
0051
0052     // get the vector of Nleptons and Ncharged
0053     auto v_NLeptons = df_with_lepton_count.Take<int>("Nleptons");
0054     auto v_NCharged = df_with_lepton_count.Take<int>("Ncharged");
0055     auto v_NStrange = df_with_lepton_count.Take<int>("NStrange");
0056     // create 3 histograms, one for the number of leptons, one for the number of charged particles, and one for the ratio of the two. All these three are as a function of the event index
0057     auto h_NLeptons = new TH1F("h_NLeptons", "h_NLeptons;Event index;Nleptons", v_NLeptons->size(), 0, v_NLeptons->size());
0058     auto h_NStrange = new TH1F("h_NStrange", "h_NStrange;Event index;Nstrange", v_NLeptons->size(), 0, v_NLeptons->size());
0059     auto h_NCharged = new TH1F("h_NCharged", "h_NCharged;Event index;Ncharged", v_NLeptons->size(), 0, v_NLeptons->size());
0060     auto h_NLeptons_ratio = new TH1F("h_NLeptons_ratio", "h_NLeptons_ratio;Event index;Nleptons/Ncharged", v_NLeptons->size(), 0, v_NLeptons->size());
0061     auto h_NStrange_ratio = new TH1F("h_NStrange_ratio", "h_NStrange_ratio;Event index;Nstrange/Ncharged", v_NLeptons->size(), 0, v_NLeptons->size());
0062     float avg_lepton = 0, avg_strange = 0;
0063     for (size_t i = 0; i < v_NLeptons->size(); i++)
0064     {
0065         h_NLeptons->SetBinContent(i + 1, (*v_NLeptons)[i]);
0066         h_NStrange->SetBinContent(i + 1, (*v_NStrange)[i]);
0067         h_NCharged->SetBinContent(i + 1, (*v_NCharged)[i]);
0068         h_NStrange_ratio->SetBinContent(i + 1, (float)(*v_NStrange)[i] / (*v_NCharged)[i] * 100.);
0069         h_NLeptons_ratio->SetBinContent(i + 1, (float)(*v_NLeptons)[i] / (*v_NCharged)[i] * 100.);
0070         avg_lepton += (float)(*v_NLeptons)[i] / (*v_NCharged)[i] * 100.;
0071         avg_strange += (float)(*v_NStrange)[i] / (*v_NCharged)[i] * 100.;
0072     }
0073     avg_lepton /= v_NLeptons->size();
0074     avg_strange /= v_NStrange->size();
0075
0076     TCanvas *c = new TCanvas("c", "c", 800, 700);
0077     TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1);
0078     pad1->SetBottomMargin(0);
0079     pad1->SetTopMargin(0.2);
0080     pad1->SetRightMargin(RightMargin);
0081     pad1->Draw();
0082     TPad *pad2 = new TPad("pad2", "pad2", 0, 0, 1, 0.3);
0083     pad2->SetRightMargin(RightMargin);
0084     pad2->SetTopMargin(0);
0085     pad2->SetBottomMargin(0.33);
0086     pad2->Draw();
0087     pad1->cd();
0088     pad1->SetLogy(1);
0089     h_NCharged->SetLineColor(kBlack);
0090     h_NCharged->SetLineWidth(2);
0091     h_NCharged->GetYaxis()->SetTitle("Counts");
0092     h_NCharged->GetYaxis()->SetTitleOffset(1.4);
0093     h_NCharged->GetYaxis()->SetRangeUser(h_NLeptons->GetMinimum(0) * 0.3, h_NCharged->GetMaximum() * 1.2);
0094     h_NCharged->GetYaxis()->SetTitleSize(AxisTitleSize);
0095     h_NCharged->GetYaxis()->SetLabelSize(AxisLabelSize);
0096     h_NCharged->GetXaxis()->SetTitleSize(0);
0097     h_NCharged->GetXaxis()->SetLabelSize(0);
0098     h_NCharged->Draw();
0099     h_NStrange->SetLineColor(38);
0100     h_NStrange->SetLineWidth(2);
0101     h_NStrange->Draw("same");
0102     h_NLeptons->SetLineColor(46);
0103     h_NLeptons->SetLineWidth(2);
0104     h_NLeptons->Draw("same");
0105     TLegend *leg = new TLegend(pad1->GetLeftMargin(),           //
0106                                1 - pad1->GetTopMargin() + 0.01, //
0107                                pad1->GetLeftMargin() + 0.25,    //
0108                                0.98);
0109
0110     leg->AddEntry(h_NCharged, "Number of charged particles", "l");
0111     leg->AddEntry(h_NStrange, "Number of strange particle (K_{s}^{0}, #Lambda), |#eta|<1.5", "l");
0112     leg->AddEntry(h_NLeptons, "Number of leptons (e, #mu), |#eta|<1.5", "l");
0113     leg->SetBorderSize(0);
0114     leg->Draw();
0115     pad2->cd();
0116     // pad2->SetLogy(1);
0117     float ratiomax = h_NStrange_ratio->GetMaximum();
0118     h_NStrange_ratio->SetLineColor(38);
0119     h_NStrange_ratio->SetLineWidth(2);
0120     h_NStrange_ratio->GetYaxis()->SetTitle("Fraction (%)");
0121     h_NStrange_ratio->GetYaxis()->SetTitleOffset(0.5);
0122     h_NStrange_ratio->GetYaxis()->SetRangeUser(0, ratiomax * 1.05);
0123     h_NStrange_ratio->GetYaxis()->SetNdivisions(505);
0124     h_NStrange_ratio->GetYaxis()->SetTitleSize(AxisTitleSize * textscale);
0125     h_NStrange_ratio->GetYaxis()->SetLabelSize(AxisLabelSize * textscale);
0126     h_NStrange_ratio->GetXaxis()->SetTitleSize(AxisTitleSize * textscale);
0127     h_NStrange_ratio->GetXaxis()->SetLabelSize(AxisLabelSize * textscale);
0128     h_NStrange_ratio->GetXaxis()->SetTitle("Event index");
0129     h_NStrange_ratio->GetXaxis()->SetTitleOffset(1.2);
0130     h_NStrange_ratio->Draw("l");
0131     h_NLeptons_ratio->SetLineColor(46);
0132     h_NLeptons_ratio->SetLineWidth(2);
0133     h_NLeptons_ratio->Draw("l same");
0134     TLine *l_lepton = new TLine(0, avg_lepton, v_NLeptons->size(), avg_lepton);
0135     l_lepton->SetLineColor(46);
0136     l_lepton->SetLineStyle(2);
0137     l_lepton->Draw("same");
0138     TLine *l_strange = new TLine(0, avg_strange, v_NLeptons->size(), avg_strange);
0139     l_strange->SetLineColor(38);
0140     l_strange->SetLineStyle(2);
0141     l_strange->Draw("same");
0142     TLatex *l_avg = new TLatex();
0143     l_avg->SetTextAlign(kHAlignLeft + kVAlignTop);
0144     l_avg->SetTextSize(AxisLabelSize * textscale);
0145     l_avg->DrawLatex(1, ratiomax * 1.05, Form("Avg. lepton fraction=%.3g%%; Avg. strange fraction=%.3g%%", avg_lepton, avg_strange));
0146     c->SaveAs(Form("%s/countLepton.png", plotdir.c_str()));
0147     c->SaveAs(Form("%s/countLepton.pdf", plotdir.c_str()));
0148 }