NewCode2024/macros/TwoPlot_ratio.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 TCanvas * make_comparison(
0047     std::vector<std::pair<std::string, TH1D*>> data_hist_vec_in,  // note : description, hist // note : postQA
0048     std::vector<std::pair<std::string, TH1D*>> MC_hist_vec_in,    // note : description, hist // note : noQA
0049     std::string output_directory_in,
0050     std::string output_filename_in,
0051     std::string sub_label_text_in, // note : "Data", "MC", "Comp"
0052     bool DoNormalize = false,
0053     bool set_log_y = false,
0054     bool isData_more = false,
0055     double MainPlot_y_low = 0,
0056     double MainPlot_y_top = -99,
0057     double ratio_y_low = 0,
0058     double ratio_y_top = 2
0059 )
0060 {
0061     double logy_max_factor = 1000;
0062     double lineary_max_factor = 1.6;
0063     std::string sPH_labeling = "Internal";
0064
0065     std::cout<<"In make_comparison, working on "<<output_filename_in<<std::endl;
0066
0067     for (auto &pair : data_hist_vec_in)
0068     {
0069         if (pair.second == nullptr)
0070         {
0071             std::cout<<"Error: "<<pair.first<<" has nullptr"<<std::endl;
0072             exit(1);
0073         }
0074     }
0075
0076     for (auto &pair : MC_hist_vec_in)
0077     {
0078         if (pair.second == nullptr)
0079         {
0080             std::cout<<"Error: "<<pair.first<<" has nullptr"<<std::endl;
0081             exit(1);
0082         }
0083     }
0084
0085     std::vector<int> data_hist_nbinsX; data_hist_nbinsX.clear();
0086     std::vector<int> MC_hist_nbinsX; MC_hist_nbinsX.clear();
0087
0088     int color_index = 0;
0089
0090     for (auto &pair : data_hist_vec_in)
0091     {
0092
0093         std::cout<<"Data: "<<pair.first<<", "<<pair.second->GetName()<<", "<<pair.second->GetNbinsX()<<std::endl;
0094
0095         pair.second -> SetStats(0);
0096
0097         data_hist_nbinsX.push_back(pair.second->GetNbinsX());
0098
0099         if (DoNormalize) {
0100             pair.second->Scale(1.0/pair.second->Integral());
0101             pair.second->GetYaxis()->SetTitle(
0102                 Form("%s (A.U.)", (pair.second->GetYaxis()->GetTitle()))
0103             );
0104         }
0105
0106         pair.second->SetLineColor(TColor::GetColor(color_code[color_index].c_str()));
0107         pair.second->SetMarkerColor(TColor::GetColor(color_code[color_index].c_str()));
0108         pair.second->SetFillColorAlpha(2,0);
0109         color_index += 1;
0110     }
0111
0112     for (auto &pair : MC_hist_vec_in)
0113     {
0114         std::cout<<"MC: "<<pair.first<<", "<<pair.second->GetName()<<", "<<pair.second->GetNbinsX()<<std::endl;
0115
0116         MC_hist_nbinsX.push_back(pair.second->GetNbinsX());
0117
0118         pair.second -> SetStats(0);
0119
0120         if (DoNormalize) {
0121             pair.second->Scale(1.0/pair.second->Integral());
0122             pair.second->GetYaxis()->SetTitle(
0123                 Form("%s (A.U.)", (pair.second->GetYaxis()->GetTitle()))
0124             );
0125         }
0126
0127         pair.second->SetLineColor(TColor::GetColor(color_code[color_index].c_str()));
0128         pair.second->SetMarkerColor(TColor::GetColor(color_code[color_index].c_str()));
0129         pair.second->SetFillColorAlpha(2,0);
0130         color_index += 1;
0131     }
0132
0133     for (int &data_nbin : data_hist_nbinsX){
0134         for (int &MC_nbin : MC_hist_nbinsX){
0135             if (data_nbin != MC_nbin){
0136                 std::cout<<"Error: Different binning for data and MC, data: "<<data_nbin<<", MC: "<<MC_nbin<<std::endl;
0137                 exit(1);
0138             }
0139         }
0140     }
0141
0142     bool bins_are_same = true;
0143     for (int i = 1; i <= data_hist_vec_in[0].second->GetNbinsX(); i++)
0144     {
0145         if (data_hist_vec_in[0].second->GetBinCenter(i) != MC_hist_vec_in[0].second->GetBinCenter(i))
0146         {
0147             bins_are_same = false;
0148             break;
0149         }
0150     }
0151
0152     if (MainPlot_y_top == -99)
0153     {
0154         for (auto &pair : data_hist_vec_in)
0155         {
0156             MainPlot_y_top = (pair.second->GetBinContent(pair.second->GetMaximumBin()) > MainPlot_y_top) ? pair.second->GetBinContent(pair.second->GetMaximumBin()) : MainPlot_y_top;
0157         }
0158
0159         for (auto &pair : MC_hist_vec_in)
0160         {
0161             MainPlot_y_top = (pair.second->GetBinContent(pair.second->GetMaximumBin()) > MainPlot_y_top) ? pair.second->GetBinContent(pair.second->GetMaximumBin()) : MainPlot_y_top;
0162         }
0163
0164
0165         if (set_log_y){
0166             MainPlot_y_top = MainPlot_y_top * logy_max_factor;
0167         }
0168         else{
0169             MainPlot_y_top = MainPlot_y_top * lineary_max_factor;
0170         }
0171     }
0172
0173     if (MainPlot_y_low == 0 && set_log_y)
0174     {
0175         MainPlot_y_low = 1000000;
0176
0177         for (auto &pair : data_hist_vec_in)
0178         {
0179             double min = GetNonZeroMinimalBin(pair.second);
0180
0181             MainPlot_y_low = (min < MainPlot_y_low) ? min : MainPlot_y_low;
0182         }
0183
0184         for (auto &pair : MC_hist_vec_in)
0185         {
0186             double min = GetNonZeroMinimalBin(pair.second);
0187
0188             MainPlot_y_low = (min < MainPlot_y_low) ? min : MainPlot_y_low;
0189         }
0190     }
0191
0192     TLatex * ltx = new TLatex();
0193     ltx->SetNDC();
0194     ltx->SetTextSize(0.045);
0195     ltx->SetTextAlign(31);
0196
0197     TLegend * All_leg = new TLegend(0.4, 0.82, 0.8, 0.88);
0198     All_leg -> SetBorderSize(0);
0199     All_leg -> SetTextSize(0.025);
0200
0201     TLegend * All_leg_long = new TLegend(0.7, 0.5, 0.8, 0.88);
0202     All_leg_long -> SetBorderSize(0);
0203
0204     TLatex * sub_ltx = new TLatex();
0205     sub_ltx->SetNDC();
0206     sub_ltx->SetTextSize(0.02);
0207     // ltx->SetTextAlign(31);
0208
0209     TH1D * h1D_pad1 = (TH1D*)data_hist_vec_in.front().second->Clone("h1D_pad1");
0210     TH1D * h1D_pad2 = (TH1D*)data_hist_vec_in.front().second->Clone("h1D_pad2");
0211     h1D_pad1 -> Reset("ICESM");
0212     h1D_pad2 -> Reset("ICESM");
0213     h1D_pad1 -> SetLineColorAlpha(1,0.);
0214     h1D_pad2 -> SetLineColorAlpha(1,0.);
0215     h1D_pad1 -> SetFillColorAlpha(1,0.);
0216     h1D_pad2 -> SetFillColorAlpha(1,0.);
0217     h1D_pad1 -> SetMarkerColorAlpha(1,0.);
0218     h1D_pad2 -> SetMarkerColorAlpha(1,0.);
0219
0220
0221     SetsPhenixStyle();
0222
0223     TCanvas * c1 = new TCanvas("", "c1", 950, 800);
0224
0225     TPad * pad2 = new TPad("", "pad2", 0.0, 0.0, 1.0, 0.29);
0226     pad2 -> SetBottomMargin(0.3);
0227     pad2->Draw();
0228
0229     TPad * pad1 = new TPad("", "pad1", 0.0, 0.29, 1.0, 1.0);
0230     pad1 -> SetTopMargin(0.065);
0231     pad1 -> SetBottomMargin(0.02);
0232     pad1->Draw();
0233
0234     int total_hist = data_hist_vec_in.size() + MC_hist_vec_in.size();
0235
0236     // todo: set the size fo the legend
0237     double leg_x1 = (total_hist <= 4) ? 0.45 : 0.65;
0238     double leg_y1 = (total_hist <= 4) ? 0.76 : 0.4;
0239     double leg_x2 = (total_hist <= 4) ? 0.8 : 0.8;
0240     double leg_y2 = (total_hist <= 4) ? 0.88 : 0.88;
0241
0242     TLegend * leg = new TLegend(leg_x1, leg_y1, leg_x2, leg_y2);
0243     leg -> SetBorderSize(0);
0244     leg -> SetTextSize(0.025);
0245
0246     h1D_pad1->GetXaxis()->SetLabelOffset(999);
0247     h1D_pad1->GetXaxis()->SetLabelSize(0);
0248     h1D_pad1->GetXaxis()->SetTitleOffset(999);
0249     h1D_pad1->GetYaxis()->SetLabelSize(0.045);
0250     h1D_pad1->GetYaxis()->SetTitleOffset(0.8);
0251     h1D_pad1->GetYaxis()->SetTitleSize(0.07);
0252     h1D_pad1->GetXaxis()->SetNdivisions();
0253
0254     for (auto &pair : MC_hist_vec_in){
0255         leg -> AddEntry(pair.second, Form("%s, Mean: %.2f, StdDev: %.2f", pair.first.c_str(), pair.second->GetMean(), pair.second->GetStdDev()), "f");
0256     }
0257
0258     for (auto &pair : data_hist_vec_in){
0259         // pair.second->GetXaxis()->SetNdivisions();
0260         leg -> AddEntry(pair.second, Form("%s, Mean: %.2f, StdDev: %.2f", pair.first.c_str(), pair.second->GetMean(), pair.second->GetStdDev()), "lep");
0261     }
0262
0263     sub_ltx->DrawLatex(0.45, 0.78, Form("Mean diff (A - B)/B: %.3f", (data_hist_vec_in[0].second->GetMean() - MC_hist_vec_in[0].second->GetMean())/MC_hist_vec_in[0].second->GetMean()));
0264     sub_ltx->DrawLatex(0.45, 0.74, Form("StdDev diff (A - B)/B: %.3f", (data_hist_vec_in[0].second->GetStdDev() - MC_hist_vec_in[0].second->GetStdDev())/MC_hist_vec_in[0].second->GetStdDev()));
0265
0266     pad1->cd();
0267
0268     std::cout<<"For "<<output_filename_in<<", MainPlot_y_low: "<<MainPlot_y_low<<", MainPlot_y_top: "<<MainPlot_y_top<<std::endl;
0269
0270     h1D_pad1 -> SetMinimum(MainPlot_y_low);
0271     h1D_pad1 -> SetMaximum(MainPlot_y_top);
0272     h1D_pad1 -> Draw("hist");
0273
0274     for (int i = 0; i < MC_hist_vec_in.size(); i++){
0275         MC_hist_vec_in[i].second -> Draw("hist same");
0276     }
0277
0278     for (int i = 0; i < data_hist_vec_in.size(); i++){
0279         data_hist_vec_in[i].second -> Draw("ep same");
0280     }
0281
0282     leg -> Draw("same");
0283
0284     pad1->SetLogy(set_log_y);
0285
0286     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", sPH_labeling.c_str()));
0287
0288     // if (labelling_vec_map.find(sub_label_text_in) != labelling_vec_map.end())
0289     // {
0290     //     for (auto &tuple : labelling_vec_map[sub_label_text_in])
0291     //     {
0292     //         double x1 = std::get<0>(tuple);
0293     //         double y1 = std::get<1>(tuple);
0294     //         std::string text = std::get<2>(tuple);
0295
0296     //         sub_ltx->DrawLatex(x1, y1, text.c_str());
0297     //     }
0298     // }
0299
0300     TLine * line = new TLine();
0301     line -> SetLineStyle(2);
0302     line -> SetLineWidth(2);
0303     line -> SetLineColor(28);
0304
0305     pad2->cd();
0306
0307     std::vector<TH1D*> ratio_hist_vec_in; ratio_hist_vec_in.clear();
0308
0309     std::vector<std::pair<std::string, TH1D*>> hist_vec_in = (isData_more) ? data_hist_vec_in : MC_hist_vec_in;
0310
0311     for (int i = 0; i < hist_vec_in.size(); i++){
0312
0313         std::pair<std::string, TH1D*> pair = hist_vec_in[i];
0314
0315         // todo: can only do the ratio w.r.t the first data hist
0316         if (isData_more){
0317             ratio_hist_vec_in.push_back( // note : hist_vec_in: data
0318                 (TH1D*) MC_hist_vec_in[0].second->Clone(Form("ratio_hist_%s", pair.first.c_str()))
0319             );
0320
0321             ratio_hist_vec_in.back() -> Sumw2(true);
0322             if (bins_are_same) {ratio_hist_vec_in.back() -> Divide(pair.second, ratio_hist_vec_in.back());}
0323
0324         }
0325         else { // note : hist_vec_in: MC
0326             ratio_hist_vec_in.push_back(
0327                 (TH1D*) data_hist_vec_in[0].second->Clone(Form("ratio_hist_%s", pair.first.c_str()))
0328             );
0329
0330             ratio_hist_vec_in.back() -> Sumw2(true);
0331             if (bins_are_same) {ratio_hist_vec_in.back() -> Divide(pair.second);}
0332         }
0333
0334
0335         // if (ratio_hist_vec_in.back()->GetXaxis()->GetLabelOffset() > 10)
0336         // {
0337         //     ratio_hist_vec_in.back()->GetXaxis()->SetLabelOffset(0.01);
0338         // }
0339
0340         ratio_hist_vec_in.back()->SetMarkerStyle(20);
0341         ratio_hist_vec_in.back()->SetMarkerSize(0.8);
0342         ratio_hist_vec_in.back()->SetMarkerColor(pair.second->GetLineColor());
0343         ratio_hist_vec_in.back()->SetLineColor(pair.second->GetLineColor());
0344
0345         if (i == 0){
0346             h1D_pad2->GetXaxis()->SetLabelSize(0.11);
0347             h1D_pad2->GetXaxis()->SetTitleOffset(0.8);
0348             h1D_pad2->GetXaxis()->SetTitleSize(0.14);
0349             h1D_pad2->GetXaxis()->SetNdivisions();
0350
0351             h1D_pad2->GetYaxis()->SetLabelSize(0.11);
0352             h1D_pad2->GetYaxis()->SetTitleOffset(0.4);
0353             h1D_pad2->GetYaxis()->SetTitleSize(0.14);
0354             h1D_pad2->GetYaxis()->SetNdivisions(505);
0355
0356             h1D_pad2 -> SetMinimum(ratio_y_low);
0357             h1D_pad2 -> SetMaximum(ratio_y_top);
0358
0359             h1D_pad2->GetYaxis()->SetTitle("Data/MC");
0360
0361
0362             h1D_pad2->Draw("ep");
0363             line -> DrawLine(h1D_pad2->GetXaxis()->GetXmin(), 1, h1D_pad2->GetXaxis()->GetXmax(), 1);
0364
0365             ratio_hist_vec_in.back()->Draw("ep same");
0366         }
0367         else {
0368             ratio_hist_vec_in.back()->Draw("ep same");
0369         }
0370     }
0371
0372
0373
0374     // std::string data_hist_name = data_hist->GetName();
0375
0376
0377     // 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)
0378     // {
0379     //     std::cout<<"what? : "<<data_hist_name<<std::endl;
0380
0381     //     // c1 -> Print(Form("%s/Comp_%s.pdf", (output_directory + "/" + sensor_folder).c_str(), output_filename.c_str()));
0382     // }
0383     if (true)
0384     {
0385         c1 -> Print(Form("%s/Comp_%s.pdf", output_directory_in.c_str(), output_filename_in.c_str()));
0386     }
0387
0388     return c1;
0389
0390     // 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 {};}
0391
0392     // return GetHistYInfo(ratio_hist);
0393
0394
0395     // pad1 -> Clear();
0396     // pad2 -> Clear();
0397     // c1 -> Clear();
0398
0399     // delete pad1;
0400     // delete pad2;
0401     // delete c1;
0402     // delete ratio_hist;
0403
0404     // return;
0405 }
0406
0407 TH1D * GetHist(std::string directory, std::string file_name, std::string hist_name)
0408 {
0409     TFile * file = new TFile(Form("%s/%s", directory.c_str(), file_name.c_str()), "READ");
0410     TH1D * hist = (TH1D*)file->Get(hist_name.c_str());
0411
0412     return hist;
0413 }
0414
0415 int TwoPlot_ratio()
0416 {
0417     std::string black_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run6_EvtZFitWidthChange/EvtVtxZ/old_folder/FinalResult_10cm_Pol2BkgFit/completed/vtxZ_-10_10cm_MBin14/Folder_BaseLine/Run_0/completed";
0418     std::string black_file_name = "Data_PreparedNdEtaEach_AlphaCorr_AllSensor_VtxZ10_Mbin14_00054280_00000_dNdEta.root";
0419     std::string black_hist_name = "h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC";
0420     std::string black_str = "Old";
0421
0422     std::string red_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run6_EvtZFitWidthChange/EvtVtxZ/FinalResult_10cm_Pol2BkgFit/completed/vtxZ_-10_10cm_MBin14/Folder_DeltaPhiCut/Run_0/completed";
0423     std::string red_file_name = "Data_PreparedNdEtaEach_AlphaCorr_AllSensor_VtxZ10_Mbin14_00054280_00000_dNdEta.root";
0424     std::string red_hist_name = "h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC";
0425     std::string red_str = "New";
0426
0427     std::string output_directory = black_directory;
0428     // std::string output_filename = "comp_h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC.root";
0429     std::string output_filename_pdf = "comp_data_h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC";
0430
0431     bool DoNormalize = false;
0432     bool set_log_y = false;
0433
0434     // TFile * file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0435
0436     TH1D * black_hist = GetHist(black_directory, black_file_name, black_hist_name);
0437     TH1D * red_hist = GetHist(red_directory, red_file_name, red_hist_name);
0438
0439     TCanvas * c1 = make_comparison(
0440         {{black_str, black_hist}},
0441         {{red_str, red_hist}},
0442
0443         output_directory,
0444         output_filename_pdf,
0445         "", // note : "Data", "MC", "Comp"
0446         DoNormalize, // note : DoNormalize
0447         set_log_y,  // note : set_log_y
0448         false, // note : isData_more
0449         0,  // note : MainPlot_y_low
0450         -99,  // note : MainPlot_y_top
0451         0.5, // note : ratio_y_low
0452         1.5 // note : ratio_y_top
0453     );
0454
0455     // file_out -> cd();
0456
0457     // c1 -> Write();
0458
0459     // file_out -> Close();
0460
0461     return 0;
0462 }