macro/MakePlot/MakePlot.C

0001 #include "sPhenixStyle.C"
0002
0003 const std::vector<std::string> color_code = {
0004     "#000000",
0005     "#9e0142",
0006     "#5e4fa2",
0007     "#66c2a5",
0008     "#fdae61",
0009     "#0076A7",
0010     "#fee08b",
0011
0012     "#d53e4f",
0013     "#abdda4",
0014     "#f46d43",
0015     "#e6f598",
0016     "#00A1FF",
0017
0018     "#1b9e77",  // Teal Green
0019     "#d95f02",  // Deep Orange
0020     "#7570b3",  // Muted Purple
0021     "#e7298a",  // Vivid Pink
0022     "#66a61e",  // Olive Green
0023     "#a6761d",  // Brown
0024     "#666666",  // Dark Gray
0025     "#f781bf",  // Light Pink
0026     "#0173b2",  // Rich Blue
0027     "#ff7f00"   // Bright Orange
0028
0029 };
0030
0031 double GetNonZeroMinimalBin(TH1D * hist_in)
0032 {
0033     double min = 10000000;
0034     for (int i = 1; i <= hist_in->GetNbinsX(); i++)
0035     {
0036         if (hist_in->GetBinContent(i) > 0 && hist_in->GetBinContent(i) < min)
0037         {
0038             min = hist_in->GetBinContent(i);
0039         }
0040     }
0041
0042     return min;
0043 }
0044
0045
0046 std::pair<double,double> make_comparison(
0047     bool isData_in,
0048     std::vector<std::pair<std::string, TH1D*>> data_hist_vec_in,  // note : description, hist // note : postQA
0049     std::vector<std::pair<std::string, TH1D*>> MC_hist_vec_in,    // note : description, hist // note : noQA
0050     std::string output_directory_in,
0051     std::string output_filename_in,
0052     std::string sub_label_text_in, // note : "Data", "MC", "Comp"
0053     bool DoNormalize = false,
0054     bool set_log_y = false,
0055     bool isData_more = false,
0056     double MainPlot_y_low = 0,
0057     double MainPlot_y_top = -99,
0058     double ratio_y_low = 0,
0059     double ratio_y_top = 2
0060 )
0061 {
0062     double logy_max_factor = 1000;
0063     double lineary_max_factor = 1.6;
0064     std::string sPH_labeling = (isData_in) ? "Internal" : "Simulation";
0065
0066     std::cout<<"In make_comparison, working on "<<output_filename_in<<std::endl;
0067
0068     for (auto &pair : data_hist_vec_in)
0069     {
0070         if (pair.second == nullptr)
0071         {
0072             std::cout<<"Error: "<<pair.first<<" has nullptr"<<std::endl;
0073             exit(1);
0074         }
0075     }
0076
0077     for (auto &pair : MC_hist_vec_in)
0078     {
0079         if (pair.second == nullptr)
0080         {
0081             std::cout<<"Error: "<<pair.first<<" has nullptr"<<std::endl;
0082             exit(1);
0083         }
0084     }
0085
0086     std::vector<int> data_hist_nbinsX; data_hist_nbinsX.clear();
0087     std::vector<int> MC_hist_nbinsX; MC_hist_nbinsX.clear();
0088
0089     int color_index = 0;
0090
0091     for (auto &pair : data_hist_vec_in)
0092     {
0093
0094         std::cout<<"Data: "<<pair.first<<", "<<pair.second->GetName()<<", "<<pair.second->GetNbinsX()<<std::endl;
0095
0096         data_hist_nbinsX.push_back(pair.second->GetNbinsX());
0097
0098         if (DoNormalize) {
0099             pair.second->Scale(1.0/pair.second->Integral());
0100             pair.second->GetYaxis()->SetTitle(
0101                 Form("%s (A.U.)", (pair.second->GetYaxis()->GetTitle()))
0102             );
0103         }
0104
0105         pair.second->SetLineColor(TColor::GetColor(color_code[color_index].c_str()));
0106         pair.second->SetMarkerColor(TColor::GetColor(color_code[color_index].c_str()));
0107         color_index += 1;
0108     }
0109
0110     for (auto &pair : MC_hist_vec_in)
0111     {
0112         std::cout<<"MC: "<<pair.first<<", "<<pair.second->GetName()<<", "<<pair.second->GetNbinsX()<<std::endl;
0113
0114         MC_hist_nbinsX.push_back(pair.second->GetNbinsX());
0115
0116         if (DoNormalize) {
0117             pair.second->Scale(1.0/pair.second->Integral());
0118             pair.second->GetYaxis()->SetTitle(
0119                 Form("%s (A.U.)", (pair.second->GetYaxis()->GetTitle()))
0120             );
0121         }
0122
0123         pair.second->SetLineColor(TColor::GetColor(color_code[color_index].c_str()));
0124         pair.second->SetMarkerColor(TColor::GetColor(color_code[color_index].c_str()));
0125         color_index += 1;
0126     }
0127
0128     for (int &data_nbin : data_hist_nbinsX){
0129         for (int &MC_nbin : MC_hist_nbinsX){
0130             if (data_nbin != MC_nbin){
0131                 std::cout<<"Error: Different binning for data and MC, data: "<<data_nbin<<", MC: "<<MC_nbin<<std::endl;
0132                 exit(1);
0133             }
0134         }
0135     }
0136
0137     if (MainPlot_y_top == -99)
0138     {
0139         for (auto &pair : data_hist_vec_in)
0140         {
0141             MainPlot_y_top = (pair.second->GetBinContent(pair.second->GetMaximumBin()) > MainPlot_y_top) ? pair.second->GetBinContent(pair.second->GetMaximumBin()) : MainPlot_y_top;
0142         }
0143
0144         for (auto &pair : MC_hist_vec_in)
0145         {
0146             MainPlot_y_top = (pair.second->GetBinContent(pair.second->GetMaximumBin()) > MainPlot_y_top) ? pair.second->GetBinContent(pair.second->GetMaximumBin()) : MainPlot_y_top;
0147         }
0148
0149
0150         if (set_log_y){
0151             MainPlot_y_top = MainPlot_y_top * logy_max_factor;
0152         }
0153         else{
0154             MainPlot_y_top = MainPlot_y_top * lineary_max_factor;
0155         }
0156     }
0157
0158     if (MainPlot_y_low == 0 && set_log_y)
0159     {
0160         MainPlot_y_low = 1000000;
0161
0162         for (auto &pair : data_hist_vec_in)
0163         {
0164             double min = GetNonZeroMinimalBin(pair.second);
0165
0166             MainPlot_y_low = (min < MainPlot_y_low) ? min : MainPlot_y_low;
0167         }
0168
0169         for (auto &pair : MC_hist_vec_in)
0170         {
0171             double min = GetNonZeroMinimalBin(pair.second);
0172
0173             MainPlot_y_low = (min < MainPlot_y_low) ? min : MainPlot_y_low;
0174         }
0175     }
0176
0177     TLatex * ltx = new TLatex();
0178     ltx->SetNDC();
0179     ltx->SetTextSize(0.045);
0180     ltx->SetTextAlign(31);
0181
0182     TLegend * All_leg = new TLegend(0.4, 0.82, 0.8, 0.88);
0183     All_leg -> SetBorderSize(0);
0184     All_leg -> SetTextSize(0.025);
0185
0186     TLegend * All_leg_long = new TLegend(0.7, 0.5, 0.8, 0.88);
0187     All_leg_long -> SetBorderSize(0);
0188
0189     TLatex * sub_ltx = new TLatex();
0190     sub_ltx->SetNDC();
0191     sub_ltx->SetTextSize(0.02);
0192     // ltx->SetTextAlign(31);
0193
0194     TH1D * h1D_pad1 = (TH1D*)data_hist_vec_in.front().second->Clone("h1D_pad1");
0195     TH1D * h1D_pad2 = (TH1D*)data_hist_vec_in.front().second->Clone("h1D_pad2");
0196     h1D_pad1 -> Reset("ICESM");
0197     h1D_pad2 -> Reset("ICESM");
0198     h1D_pad1 -> SetLineColorAlpha(1,0.);
0199     h1D_pad2 -> SetLineColorAlpha(1,0.);
0200     h1D_pad1 -> SetFillColorAlpha(1,0.);
0201     h1D_pad2 -> SetFillColorAlpha(1,0.);
0202     h1D_pad1 -> SetMarkerColorAlpha(1,0.);
0203     h1D_pad2 -> SetMarkerColorAlpha(1,0.);
0204
0205
0206     SetsPhenixStyle();
0207
0208     TCanvas * c1 = new TCanvas("", "c1", 950, 800);
0209
0210     TPad * pad2 = new TPad("", "pad2", 0.0, 0.0, 1.0, 0.29);
0211     pad2 -> SetBottomMargin(0.3);
0212     pad2->Draw();
0213
0214     TPad * pad1 = new TPad("", "pad1", 0.0, 0.29, 1.0, 1.0);
0215     pad1 -> SetTopMargin(0.065);
0216     pad1 -> SetBottomMargin(0.02);
0217     pad1->Draw();
0218
0219     int total_hist = data_hist_vec_in.size() + MC_hist_vec_in.size();
0220
0221     // todo: set the size fo the legend
0222     double leg_x1 = (total_hist <= 4) ? 0.45 : 0.65;
0223     double leg_y1 = (total_hist <= 4) ? 0.76 : 0.4;
0224     double leg_x2 = (total_hist <= 4) ? 0.8 : 0.8;
0225     double leg_y2 = (total_hist <= 4) ? 0.88 : 0.88;
0226
0227     TLegend * leg = new TLegend(leg_x1, leg_y1, leg_x2, leg_y2);
0228     leg -> SetBorderSize(0);
0229     leg -> SetTextSize(0.025);
0230
0231     h1D_pad1->GetXaxis()->SetLabelOffset(999);
0232     h1D_pad1->GetXaxis()->SetLabelSize(0);
0233     h1D_pad1->GetXaxis()->SetTitleOffset(999);
0234     h1D_pad1->GetYaxis()->SetLabelSize(0.045);
0235     h1D_pad1->GetYaxis()->SetTitleOffset(0.8);
0236     h1D_pad1->GetYaxis()->SetTitleSize(0.07);
0237     h1D_pad1->GetXaxis()->SetNdivisions();
0238
0239     for (auto &pair : MC_hist_vec_in){
0240         leg -> AddEntry(pair.second, Form("%s, Mean: %.2f, StdDev: %.2f", pair.first.c_str(), pair.second->GetMean(), pair.second->GetStdDev()), "f");
0241     }
0242
0243     for (auto &pair : data_hist_vec_in){
0244         // pair.second->GetXaxis()->SetNdivisions();
0245         leg -> AddEntry(pair.second, Form("%s, Mean: %.2f, StdDev: %.2f", pair.first.c_str(), pair.second->GetMean(), pair.second->GetStdDev()), "lep");
0246     }
0247
0248     sub_ltx->DrawLatex(0.45, 0.78, Form("Mean diff (PostQA - NoQA)/NoQA: %.3f", (data_hist_vec_in[0].second->GetMean() - MC_hist_vec_in[0].second->GetMean())/MC_hist_vec_in[0].second->GetMean()));
0249     sub_ltx->DrawLatex(0.45, 0.74, Form("StdDev diff (PostQA - NoQA)/NoQA: %.3f", (data_hist_vec_in[0].second->GetStdDev() - MC_hist_vec_in[0].second->GetStdDev())/MC_hist_vec_in[0].second->GetStdDev()));
0250
0251     pad1->cd();
0252
0253     std::cout<<"For "<<output_filename_in<<", MainPlot_y_low: "<<MainPlot_y_low<<", MainPlot_y_top: "<<MainPlot_y_top<<std::endl;
0254
0255     h1D_pad1 -> SetMinimum(MainPlot_y_low);
0256     h1D_pad1 -> SetMaximum(MainPlot_y_top);
0257     h1D_pad1 -> Draw("hist");
0258
0259     for (int i = 0; i < MC_hist_vec_in.size(); i++){
0260         MC_hist_vec_in[i].second -> Draw("hist same");
0261     }
0262
0263     for (int i = 0; i < data_hist_vec_in.size(); i++){
0264         data_hist_vec_in[i].second -> Draw("ep same");
0265     }
0266
0267     leg -> Draw("same");
0268
0269     pad1->SetLogy(set_log_y);
0270
0271     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", sPH_labeling.c_str()));
0272
0273     // if (labelling_vec_map.find(sub_label_text_in) != labelling_vec_map.end())
0274     // {
0275     //     for (auto &tuple : labelling_vec_map[sub_label_text_in])
0276     //     {
0277     //         double x1 = std::get<0>(tuple);
0278     //         double y1 = std::get<1>(tuple);
0279     //         std::string text = std::get<2>(tuple);
0280
0281     //         sub_ltx->DrawLatex(x1, y1, text.c_str());
0282     //     }
0283     // }
0284
0285     TLine * line = new TLine();
0286     line -> SetLineStyle(2);
0287     line -> SetLineWidth(2);
0288     line -> SetLineColor(28);
0289
0290     pad2->cd();
0291
0292     std::vector<TH1D*> ratio_hist_vec_in; ratio_hist_vec_in.clear();
0293
0294     std::vector<std::pair<std::string, TH1D*>> hist_vec_in = (isData_more) ? data_hist_vec_in : MC_hist_vec_in;
0295
0296     for (int i = 0; i < hist_vec_in.size(); i++){
0297
0298         std::pair<std::string, TH1D*> pair = hist_vec_in[i];
0299
0300         // todo: can only do the ratio w.r.t the first data hist
0301         if (isData_more){
0302             ratio_hist_vec_in.push_back( // note : hist_vec_in: data
0303                 (TH1D*) MC_hist_vec_in[0].second->Clone(Form("ratio_hist_%s", pair.first.c_str()))
0304             );
0305
0306             ratio_hist_vec_in.back() -> Sumw2(true);
0307             ratio_hist_vec_in.back() -> Divide(pair.second, ratio_hist_vec_in.back());
0308
0309         }
0310         else { // note : hist_vec_in: MC
0311             ratio_hist_vec_in.push_back(
0312                 (TH1D*) data_hist_vec_in[0].second->Clone(Form("ratio_hist_%s", pair.first.c_str()))
0313             );
0314
0315             ratio_hist_vec_in.back() -> Sumw2(true);
0316             ratio_hist_vec_in.back() -> Divide(pair.second);
0317         }
0318
0319
0320         // if (ratio_hist_vec_in.back()->GetXaxis()->GetLabelOffset() > 10)
0321         // {
0322         //     ratio_hist_vec_in.back()->GetXaxis()->SetLabelOffset(0.01);
0323         // }
0324
0325         ratio_hist_vec_in.back()->SetMarkerStyle(20);
0326         ratio_hist_vec_in.back()->SetMarkerSize(0.8);
0327         ratio_hist_vec_in.back()->SetMarkerColor(pair.second->GetLineColor());
0328         ratio_hist_vec_in.back()->SetLineColor(pair.second->GetLineColor());
0329
0330         if (i == 0){
0331             h1D_pad2->GetXaxis()->SetLabelSize(0.11);
0332             h1D_pad2->GetXaxis()->SetTitleOffset(0.8);
0333             h1D_pad2->GetXaxis()->SetTitleSize(0.14);
0334             h1D_pad2->GetXaxis()->SetNdivisions();
0335
0336             h1D_pad2->GetYaxis()->SetLabelSize(0.11);
0337             h1D_pad2->GetYaxis()->SetTitleOffset(0.4);
0338             h1D_pad2->GetYaxis()->SetTitleSize(0.14);
0339             h1D_pad2->GetYaxis()->SetNdivisions(505);
0340
0341             h1D_pad2 -> SetMinimum(ratio_y_low);
0342             h1D_pad2 -> SetMaximum(ratio_y_top);
0343
0344             h1D_pad2->GetYaxis()->SetTitle("Data/MC");
0345
0346
0347             h1D_pad2->Draw("ep");
0348             line -> DrawLine(h1D_pad2->GetXaxis()->GetXmin(), 1, h1D_pad2->GetXaxis()->GetXmax(), 1);
0349
0350             ratio_hist_vec_in.back()->Draw("ep same");
0351         }
0352         else {
0353             ratio_hist_vec_in.back()->Draw("ep same");
0354         }
0355     }
0356
0357
0358
0359     // std::string data_hist_name = data_hist->GetName();
0360
0361
0362     // if (data_hist_name.find("_l") != std::string::npos && data_hist_name.find("_phi") != std::string::npos && data_hist_name.find("_z") != std::string::npos)
0363     // {
0364     //     std::cout<<"what? : "<<data_hist_name<<std::endl;
0365
0366     //     // c1 -> Print(Form("%s/Comp_%s.pdf", (output_directory + "/" + sensor_folder).c_str(), output_filename.c_str()));
0367     // }
0368     if (true)
0369     {
0370         c1 -> Print(Form("%s/Comp_%s.pdf", output_directory_in.c_str(), output_filename_in.c_str()));
0371     }
0372
0373     return {};
0374
0375     // if (data_hist_name.find("_l") == std::string::npos || data_hist_name.find("_phi") == std::string::npos || data_hist_name.find("_z") == std::string::npos) {return {};}
0376
0377     // return GetHistYInfo(ratio_hist);
0378
0379
0380     // pad1 -> Clear();
0381     // pad2 -> Clear();
0382     // c1 -> Clear();
0383
0384     // delete pad1;
0385     // delete pad2;
0386     // delete c1;
0387     // delete ratio_hist;
0388
0389     // return;
0390 }
0391
0392 int MakePlot()
0393 {
0394     // string input_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run6_EvtZFitWidthChange/EvtVtxZ/completed/VtxZQABias/completed";
0395     // string input_file_name = "Data_InttVtxZBias_EvtBcoFullDiffCut61_00054280_merged.root";
0396     // bool isData = true;
0397
0398     string input_directory = "/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/Run7/EvtVtxZ/completed/VtxZQABias/completed";
0399     string input_file_name = "MC_InttVtxZBias_merged.root";
0400     bool isData = false;
0401
0402     std::string output_directory = input_directory + "/CompPlot";
0403
0404     system(Form("if [ ! -d %s ]; then mkdir -p %s; fi", output_directory.c_str(), output_directory.c_str()));
0405
0406     std::vector<std::string> comp_name = {
0407         "h1D_centrality",
0408         "h1D_NClus"
0409     };
0410
0411     for (int i = 0; i < 15; i++)
0412     {
0413         comp_name.push_back(Form("h1D_NClus_Mbin%d",i));
0414     }
0415
0416     TFile * file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_file_name.c_str()), "READ");
0417     if (!file_in)
0418     {
0419         std::cout<<"Error: file_in is nullptr"<<std::endl;
0420         exit(1);
0421     }
0422
0423     for (std::string hist_name : comp_name)
0424     {
0425         std::cout<<"Working on "<<hist_name<<std::endl;
0426         make_comparison(
0427             isData,
0428             {{"PostQA", (TH1D*)file_in->Get(Form("PostQA_%s", hist_name.c_str()))}},
0429             {{"NoQA", (TH1D*)file_in->Get(Form("NoQA_%s", hist_name.c_str()))}},
0430
0431             output_directory,
0432             hist_name,
0433             "", // note : "Data", "MC", "Comp"
0434             true, // note : DoNormalize
0435             false,  // note : set_log_y
0436             false, // note : isData_more
0437             0,  // note : MainPlot_y_low
0438             -99,  // note : MainPlot_y_top
0439             0.5, // note : ratio_y_low
0440             1.5 // note : ratio_y_top
0441         );
0442     }
0443
0444
0445
0446     return 888;
0447 }