dNdEta_Run2023/macros/EvtInttBCODiff.C

0001 #include <numeric>
0002
0003 float TickSize = 0.03;
0004 float AxisTitleSize = 0.05;
0005 float AxisLabelSize = 0.04;
0006 float LeftMargin = 0.15;
0007 float RightMargin = 0.05;
0008 float TopMargin = 0.08;
0009 float BottomMargin = 0.13;
0010 float textscale = 2.7;
0011
0012 void draw2Dhist(TH2 *hM, const char *xtitle, const char *ytitle, bool logz, vector<const char *> plotinfo, const char *plotname)
0013 {
0014     float ytitleoffset = 1.6;
0015
0016     TCanvas *c = new TCanvas("c", "c", 800, 700);
0017     gPad->SetRightMargin(0.15);
0018     gPad->SetTopMargin(0.07);
0019     c->cd();
0020     if (logz)
0021         c->SetLogz();
0022     hM->GetXaxis()->SetTitle(xtitle);
0023     hM->GetYaxis()->SetTitle(ytitle);
0024     hM->GetYaxis()->SetTitleOffset(ytitleoffset);
0025     // only maximal 3 digits for the axis labels
0026     TGaxis::SetMaxDigits(3);
0027     hM->Draw("colz");
0028     TLatex *t = new TLatex();
0029     t->SetTextSize(0.035);
0030     for (size_t i = 0; i < plotinfo.size(); ++i)
0031     {
0032         t->DrawLatexNDC(0.2, 0.85 - i * 0.04, plotinfo[i]);
0033     }
0034     c->SaveAs(Form("%s.pdf", plotname));
0035     c->SaveAs(Form("%s.png", plotname));
0036 }
0037
0038 void EvtInttBCODiff(const char *filename = "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Data_Run54280_20241204_ProdA2024/ntuple_00000.root", //
0039                     const char *outfile = "./test.root",                                                                              //
0040                     uint64_t bcodiffcut = 61,                                                                                         //
0041                     bool debug = true,                                                                                                //
0042                     bool snapshot = false,                                                                                            //
0043                     bool diagnosticplot = true                                                                                        //
0044 )
0045 {
0046     TStopwatch timer;
0047     timer.Start();
0048
0049     // ROOT::EnableImplicitMT();
0050     ROOT::RDataFrame df("EventTree", filename);
0051
0052     auto br_INTT_BCO = df.Take<uint64_t>("INTT_BCO");
0053     auto br_GL1Packet_BCO = df.Take<uint64_t>("GL1Packet_BCO");
0054     // difference w.r.t. the previous element
0055     std::vector<uint64_t> INTT_BCO_diff_prev(br_INTT_BCO->size());
0056     std::adjacent_difference(br_INTT_BCO->begin(), br_INTT_BCO->end(), INTT_BCO_diff_prev.begin());
0057     INTT_BCO_diff_prev[0] = 0;
0058     std::vector<uint64_t> INTT_BCO_diff_next(br_INTT_BCO->size());
0059     std::copy(INTT_BCO_diff_prev.begin() + 1, INTT_BCO_diff_prev.end(), INTT_BCO_diff_next.begin());
0060     INTT_BCO_diff_next.back() = 0;
0061
0062     std::vector<uint64_t> GL1Packet_BCO_diff_prev(br_GL1Packet_BCO->size());
0063     std::adjacent_difference(br_GL1Packet_BCO->begin(), br_GL1Packet_BCO->end(), GL1Packet_BCO_diff_prev.begin());
0064     GL1Packet_BCO_diff_prev[0] = 0;
0065     std::vector<uint64_t> GL1Packet_BCO_diff_next(br_GL1Packet_BCO->size());
0066     std::copy(GL1Packet_BCO_diff_prev.begin() + 1, GL1Packet_BCO_diff_prev.end(), GL1Packet_BCO_diff_next.begin());
0067     GL1Packet_BCO_diff_next.back() = 0;
0068
0069     if (debug)
0070     {
0071         for (size_t i = 0; i < 10; ++i)
0072         {
0073             std::cout << "INTT_BCO[" << i << "] = " << br_INTT_BCO->at(i) << ", INTT_BCO_diff_prev[" << i << "] = " << INTT_BCO_diff_prev[i] << ", INTT_BCO_diff_next[" << i << "] = " << INTT_BCO_diff_next[i] << std::endl;
0074         }
0075         for (size_t i = br_INTT_BCO->size() - 10; i < br_INTT_BCO->size(); ++i)
0076         {
0077             std::cout << "INTT_BCO[" << i << "] = " << br_INTT_BCO->at(i) << ", INTT_BCO_diff_prev[" << i << "] = " << INTT_BCO_diff_prev[i] << ", INTT_BCO_diff_next[" << i << "] = " << INTT_BCO_diff_next[i] << std::endl;
0078         }
0079     }
0080
0081     auto df_with_diff = df.Define("INTT_BCO_diff_prev", [&INTT_BCO_diff_prev](ULong64_t i) { return INTT_BCO_diff_prev[i]; }, {"rdfentry_"})
0082                             .Define("INTT_BCO_diff_next", [&INTT_BCO_diff_next](ULong64_t i) { return INTT_BCO_diff_next[i]; }, {"rdfentry_"})
0083                             .Define("GL1Packet_BCO_diff_prev", [&GL1Packet_BCO_diff_prev](ULong64_t i) { return GL1Packet_BCO_diff_prev[i]; }, {"rdfentry_"})
0084                             .Define("GL1Packet_BCO_diff_next", [&GL1Packet_BCO_diff_next](ULong64_t i) { return GL1Packet_BCO_diff_next[i]; }, {"rdfentry_"});
0085
0086     // auto df_with_diff_and_flag = df_with_diff.Define("InttBco_IsToBeRemoved", [bcodiffcut](uint64_t INTT_BCO_diff_prev, uint64_t INTT_BCO_diff_next) { return ((INTT_BCO_diff_prev < bcodiffcut && INTT_BCO_diff_prev > 0) || (INTT_BCO_diff_next < bcodiffcut && INTT_BCO_diff_next > 0)); }, {"INTT_BCO_diff_prev", "INTT_BCO_diff_next"});
0087     auto df_with_diff_and_flag = df_with_diff.Define("InttBco_IsToBeRemoved", [bcodiffcut](uint64_t INTT_BCO_diff_next) { return (INTT_BCO_diff_next <= bcodiffcut && INTT_BCO_diff_next > 0); }, {"INTT_BCO_diff_next"});
0088
0089     // get all the column names
0090     auto all_columns = df_with_diff_and_flag.GetColumnNames();
0091     if (debug)
0092     {
0093         cout << "All columns: " << endl;
0094         for (auto &col : all_columns)
0095         {
0096             cout << col << ", ";
0097         }
0098         cout << endl;
0099     }
0100
0101     // remove "firedTriggers_name" from the list of all columns
0102     all_columns.erase(std::remove(all_columns.begin(), all_columns.end(), "firedTriggers_name"), all_columns.end());
0103     if (snapshot)
0104         df_with_diff_and_flag.Snapshot("EventTree", outfile, all_columns);
0105
0106     if (diagnosticplot)
0107     {
0108         std::string plotdir = "./quickcheck_InttBcoDiffEvt";
0109         system(Form("mkdir -p %s", plotdir.c_str()));
0110
0111         float ytitleoffset = 1.6;
0112
0113         auto df_isMB_Trigbit10 =
0114             df_with_diff_and_flag.Filter([](bool is_min_bias, const std::vector<int> &firedTriggers, bool InttBco_IsToBeRemoved) { return is_min_bias && std::find(firedTriggers.begin(), firedTriggers.end(), 10) != firedTriggers.end() && !InttBco_IsToBeRemoved; }, {"is_min_bias", "firedTriggers", "InttBco_IsToBeRemoved"}).Define("NOuterClus", "static_cast<int>(NClus - NClus_Layer1)");
0115         auto hM_NInnerClus_NOuterClus_Trig10_IsMB = df_isMB_Trigbit10.Histo2D({"hM_NInnerClus_NOuterClus_Trig10_IsMB", "hM_NInnerClus_NOuterClus_Trig10_IsMB", 250, 0, 5000, 250, 0, 5000}, "NClus_Layer1", "NOuterClus");
0116         auto hM_NInnerClus_MBDChargeSum_Trig10_IsMB = df_isMB_Trigbit10.Histo2D({"hM_NInnerClus_MBDChargeSum_Trig10_IsMB", "hM_NInnerClus_MBDChargeSum_Trig10_IsMB", 250, 0, 5000, 250, 0, 5000}, "NClus_Layer1", "MBD_charge_sum");
0117         auto hM_NOuterClus_MBDChargeSum_Trig10_IsMB = df_isMB_Trigbit10.Histo2D({"hM_NOuterClus_MBDChargeSum_Trig10_IsMB", "hM_NOuterClus_MBDChargeSum_Trig10_IsMB", 250, 0, 5000, 250, 0, 5000}, "NOuterClus", "MBD_charge_sum");
0118
0119         auto df_isMB_Trigbit_woBcoDiffEvtRemoval = df_with_diff_and_flag.Filter([](bool is_min_bias, const std::vector<int> &firedTriggers) { return is_min_bias && std::find(firedTriggers.begin(), firedTriggers.end(), 10) != firedTriggers.end(); }, {"is_min_bias", "firedTriggers"}).Define("NOuterClus", "static_cast<int>(NClus - NClus_Layer1)");
0120         auto hM_NInnerClus_NOuterClus_Trig10_IsMB_woBcoDiffEvtRemoval = df_isMB_Trigbit_woBcoDiffEvtRemoval.Histo2D({"hM_NInnerClus_NOuterClus_Trig10_IsMB_woBcoDiffEvtRemoval", "hM_NInnerClus_NOuterClus_Trig10_IsMB_woBcoDiffEvtRemoval", 250, 0, 5000, 250, 0, 5000}, "NClus_Layer1", "NOuterClus");
0121         auto hM_NInnerClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval = df_isMB_Trigbit_woBcoDiffEvtRemoval.Histo2D({"hM_NInnerClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval", "hM_NInnerClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval", 250, 0, 5000, 250, 0, 5000}, "NClus_Layer1", "MBD_charge_sum");
0122         auto hM_NOuterClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval = df_isMB_Trigbit_woBcoDiffEvtRemoval.Histo2D({"hM_NOuterClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval", "hM_NOuterClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval", 250, 0, 5000, 250, 0, 5000}, "NOuterClus", "MBD_charge_sum");
0123
0124         auto df_isMB_Trigbit10_removed =
0125             df_with_diff_and_flag.Filter([](bool is_min_bias, const std::vector<int> &firedTriggers, bool InttBco_IsToBeRemoved) { return is_min_bias && std::find(firedTriggers.begin(), firedTriggers.end(), 10) != firedTriggers.end() && InttBco_IsToBeRemoved; }, {"is_min_bias", "firedTriggers", "InttBco_IsToBeRemoved"}).Define("NOuterClus", "static_cast<int>(NClus - NClus_Layer1)");
0126         auto hM_NInnerClus_NOuterClus_Trig10_IsMB_removed = df_isMB_Trigbit10_removed.Histo2D({"hM_NInnerClus_NOuterClus_Trig10_IsMB_removed", "hM_NInnerClus_NOuterClus_Trig10_IsMB_removed", 250, 0, 5000, 250, 0, 5000}, "NClus_Layer1", "NOuterClus");
0127         auto hM_NInnerClus_MBDChargeSum_Trig10_IsMB_removed = df_isMB_Trigbit10_removed.Histo2D({"hM_NInnerClus_MBDChargeSum_Trig10_IsMB_removed", "hM_NInnerClus_MBDChargeSum_Trig10_IsMB_removed", 250, 0, 5000, 250, 0, 5000}, "NClus_Layer1", "MBD_charge_sum");
0128         auto hM_NOuterClus_MBDChargeSum_Trig10_IsMB_removed = df_isMB_Trigbit10_removed.Histo2D({"hM_NOuterClus_MBDChargeSum_Trig10_IsMB_removed", "hM_NOuterClus_MBDChargeSum_Trig10_IsMB_removed", 250, 0, 5000, 250, 0, 5000}, "NOuterClus", "MBD_charge_sum");
0129
0130         draw2Dhist(hM_NInnerClus_NOuterClus_Trig10_IsMB.GetPtr(), "Number of inner clusters", "Number of outer clusters", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "w. abnormal event removal"}, Form("%s/NInnerClus_vs_NOuterClus_Trig10_IsMB", plotdir.c_str()));
0131         draw2Dhist(hM_NInnerClus_MBDChargeSum_Trig10_IsMB.GetPtr(), "Number of inner clusters", "MBD Charge Sum", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "w. abnormal event removal"}, Form("%s/NInnerClus_vs_MBDChargeSum_Trig10_IsMB", plotdir.c_str()));
0132         draw2Dhist(hM_NOuterClus_MBDChargeSum_Trig10_IsMB.GetPtr(), "Number of outer clusters", "MBD Charge Sum", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "w. abnormal event removal"}, Form("%s/NOuterClus_vs_MBDChargeSum_Trig10_IsMB", plotdir.c_str()));
0133         draw2Dhist(hM_NInnerClus_NOuterClus_Trig10_IsMB_woBcoDiffEvtRemoval.GetPtr(), "Number of inner clusters", "Number of outer clusters", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "w.o abnormal event removal"}, Form("%s/NInnerClus_vs_NOuterClus_Trig10_IsMB_woBcoDiffEvtRemoval", plotdir.c_str()));
0134         draw2Dhist(hM_NInnerClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval.GetPtr(), "Number of inner clusters", "MBD Charge Sum", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "w.o abnormal event removal"}, Form("%s/NInnerClus_vs_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval", plotdir.c_str()));
0135         draw2Dhist(hM_NOuterClus_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval.GetPtr(), "Number of outer clusters", "MBD Charge Sum", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "w.o abnormal event removal"}, Form("%s/NOuterClus_vs_MBDChargeSum_Trig10_IsMB_woBcoDiffEvtRemoval", plotdir.c_str()));
0136         draw2Dhist(hM_NInnerClus_NOuterClus_Trig10_IsMB_removed.GetPtr(), "Number of inner clusters", "Number of outer clusters", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "removed events"}, Form("%s/NInnerClus_vs_NOuterClus_Trig10_IsMB_removed", plotdir.c_str()));
0137         draw2Dhist(hM_NInnerClus_MBDChargeSum_Trig10_IsMB_removed.GetPtr(), "Number of inner clusters", "MBD Charge Sum", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "removed events"}, Form("%s/NInnerClus_vs_MBDChargeSum_Trig10_IsMB_removed", plotdir.c_str()));
0138         draw2Dhist(hM_NOuterClus_MBDChargeSum_Trig10_IsMB_removed.GetPtr(), "Number of outer clusters", "MBD Charge Sum", true, std::vector<const char *>{"Is MinBias & Trigger bit 10 (MBD N&S#geq2)", "removed events"}, Form("%s/NOuterClus_vs_MBDChargeSum_Trig10_IsMB_removed", plotdir.c_str()));
0139
0140         // check number of events remaining after the cut per centrality interval
0141         TH1F *hM_cent_IsMB_all = new TH1F("hM_cent_IsMB_all", "hM_cent_IsMB_all", 10, 0, 10);
0142         TH1F *hM_cent_IsMB_OuterInnerBranch = new TH1F("hM_cent_IsMB_OuterInnerBranch", "hM_cent_IsMB_OuterInnerBranch", 10, 0, 10);
0143         TH1F *hM_ratio = new TH1F("hM_ratio", "hM_ratio", 10, 0, 10);
0144         for (int i = 0; i < 10; i++)
0145         {
0146             float low = i * 10;
0147             float high = (i + 1) * 10;
0148             int Nevt_isMB_Trigbit10_woBcoDiffEvtRemoval = *df_isMB_Trigbit_woBcoDiffEvtRemoval.Filter([low, high](float MBD_centrality) { return MBD_centrality > low && MBD_centrality <= high; }, {"MBD_centrality"}).Count();
0149             int Nevt_isMB_Trigbit10_wBcoDiffEvtRemoval = *df_isMB_Trigbit10.Filter([low, high](float MBD_centrality) { return MBD_centrality > low && MBD_centrality <= high; }, {"MBD_centrality"}).Count();
0150             hM_cent_IsMB_all->SetBinContent(i + 1, Nevt_isMB_Trigbit10_woBcoDiffEvtRemoval);
0151             hM_cent_IsMB_OuterInnerBranch->SetBinContent(i + 1, Nevt_isMB_Trigbit10_wBcoDiffEvtRemoval);
0152             hM_ratio->SetBinContent(i + 1, static_cast<double>(Nevt_isMB_Trigbit10_wBcoDiffEvtRemoval) / static_cast<double>(Nevt_isMB_Trigbit10_woBcoDiffEvtRemoval) * 100.);
0153             hM_ratio->GetXaxis()->SetBinLabel(i + 1, Form("%d-%d%%", i * 10, (i + 1) * 10));
0154         }
0155
0156         TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);
0157         TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1);
0158         pad1->SetBottomMargin(0);
0159         pad1->SetTopMargin(TopMargin);
0160         pad1->SetRightMargin(RightMargin);
0161         pad1->Draw();
0162         TPad *pad2 = new TPad("pad2", "pad2", 0, 0, 1, 0.3);
0163         pad2->SetRightMargin(RightMargin);
0164         pad2->SetTopMargin(0);
0165         pad2->SetBottomMargin(0.31);
0166         pad2->Draw();
0167         pad1->cd();
0168         pad1->SetLogy(0);
0169         hM_cent_IsMB_all->GetXaxis()->SetTitleSize(0);
0170         hM_cent_IsMB_all->GetXaxis()->SetLabelSize(0);
0171         // hM_cent_IsMB_all->GetYaxis()->SetRangeUser(hM_cent_IsMB_OuterInnerBranch->GetMinimum(0) * 0.8, hM_cent_IsMB_all->GetMaximum() * 2);
0172         hM_cent_IsMB_all->GetYaxis()->SetRangeUser(0, hM_cent_IsMB_all->GetMaximum() * 1.5);
0173         hM_cent_IsMB_all->GetYaxis()->SetMoreLogLabels(true);
0174         hM_cent_IsMB_all->GetYaxis()->SetTitle("Counts");
0175         hM_cent_IsMB_all->GetYaxis()->SetTitleSize(AxisTitleSize);
0176         hM_cent_IsMB_all->GetYaxis()->SetLabelSize(AxisTitleSize);
0177         hM_cent_IsMB_all->SetMarkerSize(1.6);
0178         hM_cent_IsMB_all->Draw("hist");
0179         hM_cent_IsMB_OuterInnerBranch->SetLineColor(kRed);
0180         hM_cent_IsMB_OuterInnerBranch->SetMarkerSize(1.6);
0181         hM_cent_IsMB_OuterInnerBranch->Draw("histsame");
0182         TLegend *leg = new TLegend(0.35, 0.72, 0.7, 0.85);
0183         leg->AddEntry(hM_cent_IsMB_all, "Is MB && Trigger bit 10", "l");
0184         leg->AddEntry(hM_cent_IsMB_OuterInnerBranch, "Is MB && Trigger bit 10 && BCO removal", "l");
0185         leg->SetBorderSize(0);
0186         leg->SetFillStyle(0);
0187         leg->SetTextSize(0.055);
0188         leg->Draw();
0189         pad2->cd();
0190         // pad2->SetLogy(1);
0191         hM_ratio->SetLineColor(kBlack);
0192         hM_ratio->SetLineWidth(2);
0193         hM_ratio->GetYaxis()->SetTitle("Fraction (%)");
0194         hM_ratio->GetYaxis()->SetNdivisions(502);
0195         hM_ratio->GetYaxis()->SetTitleSize(AxisTitleSize * textscale);
0196         hM_ratio->GetYaxis()->SetTitleOffset(0.5);
0197         hM_ratio->GetYaxis()->SetLabelSize(AxisTitleSize * textscale);
0198         hM_ratio->GetYaxis()->SetRangeUser(0, 199);
0199         hM_ratio->GetXaxis()->SetTitle("Centrality interval");
0200         hM_ratio->GetXaxis()->SetTitleSize(AxisTitleSize * textscale);
0201         hM_ratio->GetXaxis()->SetTitleOffset(1.1);
0202         hM_ratio->GetXaxis()->SetLabelSize(AxisTitleSize * textscale);
0203         gStyle->SetPaintTextFormat(".2f%%");
0204         hM_ratio->SetMarkerSize(1.5 * textscale);
0205         hM_ratio->Draw("histtext");
0206         c1->SaveAs(Form("%s/fracEvtRemoval_centrality.pdf", plotdir.c_str()));
0207         c1->SaveAs(Form("%s/fracEvtRemoval_centrality.png", plotdir.c_str()));
0208     }
0209
0210     timer.Stop();
0211     timer.Print();
0212 }