Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:17 

0001 #include <TCanvas.h>
0002 #include <TFile.h>
0003 #include <TH2F.h>
0004 #include <TStyle.h>
0005 #include <TSystem.h>
0006 #include <TTree.h>
0007 #include <algorithm>
0008 #include <cstdlib>
0009 #include <iostream>
0010 #include <string>
0011 #include <vector>
0012 
0013 float TickSize = 0.03;
0014 float AxisTitleSize = 0.05;
0015 float AxisLabelSize = 0.04;
0016 float LeftMargin = 0.15;
0017 float RightMargin = 0.08;
0018 float TopMargin = 0.08;
0019 float BottomMargin = 0.13;
0020 float textscale = 2.7;
0021 
0022 std::string plotdir = "ClusPhiSize";
0023 
0024 template <typename T> T *gethist(bool dohadd, std::string histdir, std::string histname)
0025 {
0026     std::string mergedfile = Form("%s/hists_merged.root", histdir.c_str());
0027     if (dohadd)
0028     {
0029         std::string haddcmd = Form("hadd -f -j 20 %s %s/hists_00*.root", mergedfile.c_str(), histdir.c_str());
0030         std::cout << haddcmd << std::endl;
0031         system(haddcmd.c_str());
0032     }
0033 
0034     TFile *f = new TFile(mergedfile.c_str(), "READ");
0035     T *h = (T *)f->Get(histname.c_str());
0036     return h;
0037 }
0038 
0039 void draw_comparison(vector<TH1F *> vhist, bool logy, std::string xAxisName, vector<std::string> vleg, std::string plotname)
0040 {
0041     // colors to cycle through
0042     // std::vector<int> colors = {kTBriBlue, kTBriCyan, kTBriGreen, kTBriYellow, kTVibOrange, kTBriRed, kTBriPurple};
0043     gStyle->SetPalette(kLightTemperature);
0044 
0045     float r = vhist[1]->Integral(-1, -1) / vhist[0]->Integral(-1, -1);
0046     vhist[0]->Sumw2();
0047     vhist[0]->Scale(1. / vhist[0]->Integral(-1, -1));
0048     vhist[1]->Sumw2();
0049     vhist[1]->Scale(r / vhist[1]->Integral(-1, -1)); // normalize relative to default 
0050     // normalize histograms of simulation to 1
0051     for (int i = 2; i < vhist.size(); i++)
0052     {
0053         vhist[i]->Sumw2();
0054         vhist[i]->Scale(1. / vhist[i]->Integral(-1, -1));
0055     }
0056 
0057     float maxbinc = -1, minbinc = 1e9;
0058     for (int i = 0; i < vhist.size(); i++)
0059     {
0060         if (vhist[i]->GetMaximum() > maxbinc)
0061         {
0062             maxbinc = vhist[i]->GetMaximum();
0063         }
0064         if (vhist[i]->GetMinimum(0) < minbinc)
0065         {
0066             minbinc = vhist[i]->GetMinimum(0);
0067         }
0068     }
0069 
0070     // two pads: top pad for the histogram, bottom pad for the ratio plot
0071     TCanvas *c = new TCanvas("c", "c", 800, 700);
0072     TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1);
0073     pad1->SetBottomMargin(0);
0074     pad1->SetRightMargin(RightMargin);
0075     pad1->Draw();
0076     TPad *pad2 = new TPad("pad2", "pad2", 0, 0, 1, 0.3);
0077     pad2->SetRightMargin(RightMargin);
0078     pad2->SetTopMargin(0);
0079     pad2->SetBottomMargin(0.31);
0080     pad2->Draw();
0081     pad1->cd();
0082     pad1->SetLogy(logy);
0083     // vhist[0]->GetXaxis()->SetTitle(xAxisName.c_str());
0084     vhist[0]->GetXaxis()->SetTitleSize(0);
0085     vhist[0]->GetXaxis()->SetLabelSize(0);
0086     // draw overflow bin
0087     vhist[0]->GetXaxis()->SetRange(1, vhist[0]->GetNbinsX());
0088     vhist[0]->GetYaxis()->SetTitle("Normalized counts");
0089     vhist[0]->GetYaxis()->SetTitleSize(AxisTitleSize);
0090     vhist[0]->GetYaxis()->SetLabelSize(AxisTitleSize);
0091     float yscale = (logy) ? 750 : 2.0;
0092     vhist[0]->GetYaxis()->SetRangeUser(minbinc*0.5, maxbinc * yscale);
0093     vhist[0]->SetLineWidth(2);
0094     vhist[0]->SetLineColor(kBlack);
0095     vhist[0]->SetMarkerStyle(20);
0096     vhist[0]->SetMarkerSize(0.8);
0097     vhist[0]->SetMarkerColor(kBlack);
0098     vhist[0]->Draw("PE");
0099     vhist[1]->SetLineWidth(2);
0100     vhist[1]->SetLineColor(kGray+2);
0101     vhist[1]->SetMarkerStyle(21);
0102     vhist[1]->SetMarkerSize(0.8);
0103     vhist[1]->SetMarkerColor(kGray+2);
0104     vhist[1]->Draw("PE same");
0105     for (int i = 2; i < vhist.size(); i++)
0106     {
0107         vhist[i]->SetLineWidth(2);
0108         // vhist[i]->SetLineColor(colors[i - 1]);
0109         vhist[i]->Draw("hist PLC same");
0110     }
0111     vhist[0]->Draw("PE same");
0112 
0113     TLegend *leg = new TLegend(0.55, 0.81, 0.8, 0.91);
0114     leg->AddEntry((TObject *)0, "#it{#bf{sPHENIX}} Internal", "");
0115     leg->AddEntry((TObject *)0, "Au+Au #sqrt{s_{NN}}=200 GeV", "");
0116     leg->SetBorderSize(0);
0117     leg->SetFillStyle(0);
0118     leg->SetTextSize(0.04);
0119     leg->Draw();
0120     float leg_ysize = 0.04 * vleg.size();
0121     TLegend *leg2 = new TLegend(0.55, 0.8 - leg_ysize, 0.8, 0.80);
0122     for (int i = 0; i < vleg.size(); i++)
0123     {
0124         if (i<2)
0125         {
0126             leg2->AddEntry(vhist[i], vleg[i].c_str(), "PE");
0127         }
0128         else
0129         {
0130             leg2->AddEntry(vhist[i], vleg[i].c_str(), "l");
0131         }
0132     }
0133     leg2->SetBorderSize(0);
0134     leg2->SetFillStyle(0);
0135     leg2->SetTextSize(0.035);
0136     leg2->Draw();
0137 
0138     pad2->cd();
0139     for (int i = 1; i < vhist.size(); i++)
0140     {
0141         TH1F *hratio = (TH1F *)vhist[i]->Clone(Form("hratio%d", i));
0142         hratio->Divide(vhist[0]);
0143         hratio->SetLineWidth(2);
0144         hratio->GetYaxis()->SetRangeUser(0.001, 1.999);
0145         hratio->GetYaxis()->SetTitle("Data/Sim.");
0146         hratio->GetYaxis()->SetTitleSize(AxisTitleSize*textscale);
0147         hratio->GetYaxis()->SetTitleOffset(0.5);
0148         hratio->GetYaxis()->SetLabelSize(AxisTitleSize*textscale);
0149         hratio->GetXaxis()->SetTitle(xAxisName.c_str());
0150         hratio->GetXaxis()->SetTitleSize(AxisTitleSize*textscale);
0151         hratio->GetXaxis()->SetTitleOffset(1.1);
0152         hratio->GetXaxis()->SetLabelSize(AxisTitleSize*textscale);
0153         if (i == 1)
0154         {
0155             hratio->GetYaxis()->SetNdivisions(505);
0156             hratio->GetXaxis()->SetRange(1, hratio->GetNbinsX());
0157             hratio->Draw("hist L");
0158         }
0159         else
0160         {
0161             hratio->Draw("hist L PLC same");
0162         }
0163     }
0164     // redraw axis
0165     pad2->RedrawAxis();
0166     // draw a line at y=1
0167     TLine *line = new TLine(vhist[0]->GetXaxis()->GetXmin(), 1, vhist[0]->GetXaxis()->GetXmax(), 1);
0168     line->SetLineStyle(2);
0169     line->Draw();
0170 
0171     c->SaveAs(Form("%s.pdf", plotname.c_str()));
0172     c->SaveAs(Form("%s.png", plotname.c_str()));
0173 }
0174 
0175 void comp_diffusionradius()
0176 {
0177     std::vector<std::string> infiledirs = {
0178         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_Run54280_20241113/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0179         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_Run54280_20241114_Zclustering/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0180         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241102/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0181         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241105_diffusionradius10um/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0182         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241107_diffusionradius20um/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0183         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241107_diffusionradius30um/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0184         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241108_diffusionradius50um/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0185         "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241109_diffusionradius100um/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0/",
0186     };
0187 
0188     // std::vector<std::string> infiledirs = {
0189     //     // "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_CombinedNtuple_Run54280_HotChannel_BCO/Cluster/",
0190     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_Run54280_20241113/Cluster/",
0191     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_Run54280_20241114_Zclustering/Cluster/",
0192     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241102/Cluster/",
0193     //     // "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241105_diffusionradius10um/Cluster/",
0194     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241107_diffusionradius20um/Cluster/",
0195     //     // "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_Ntuple_HIJING_ana443_20241107_diffusionradius30um/Cluster/",
0196     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241108_diffusionradius50um/Cluster/",
0197     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241109_diffusionradius100um/Cluster/",
0198     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241109_diffusionradius200um/Cluster/",
0199     //     "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_HIJING_ana443_20241114_diffusionradius5um_Zclustering/Cluster/",
0200     // };
0201 
0202     std::vector<std::string> leg = {
0203         "Data, Run54280", 
0204         "Data, Run54280, Z-clustering",
0205         "HIJING, diffusion radius 5#mum", 
0206         "HIJING, diffusion radius 10#mum",
0207         "HIJING, diffusion radius 20#mum",
0208         "HIJING, diffusion radius 30#mum",
0209         "HIJING, diffusion radius 50#mum",
0210         "HIJING, diffusion radius 100#mum",
0211         // "HIJING, diffusion radius 200#mum",
0212         // "HIJING, diffusion radius 5#mum, Z-clustering",
0213     };
0214 
0215     std::vector<std::pair<std::string, std::string>> histvarnames = {
0216         // {"hM_ClusPhiSize_all", "Cluster #phi size"}, 
0217         // {"hM_ClusADC_all", "Cluster ADC"},
0218         // {"hM_ClusPhiSize_all_etale0p1", "Cluster #phi size, |#eta|<0.1"},
0219         // {"hM_ClusADC_all_etale0p1", "Cluster ADC, |#eta|<0.1"},
0220         // {"hM_ClusPhiSize_all_etale0p1_ClusZSize1", "Cluster #phi size, |#eta|<0.1, cluster Z size=1"},
0221         // {"hM_ClusADC_all_etale0p1_ClusZSize1", "Cluster ADC, |#eta|<0.1, cluster Z size=1"}
0222         // {"hM_ClusZ_all", "Cluster Z [cm]"}
0223         {"hM_Eta_reco", "Reco-tracklet #eta"},
0224         {"hM_dPhi_reco_altrange", "Reco-tracklet #Delta#phi"},
0225         {"hM_dEta_reco_altrange", "Reco-tracklet #Delta#eta"},
0226         {"hM_cluseta", "Cluster #eta"}
0227     };
0228 
0229     std::map<std::pair<std::string, std::string>, std::vector<TH1F *>> histmap;
0230     for (int i = 0; i < infiledirs.size(); i++)
0231     {
0232         for (int j = 0; j < histvarnames.size(); j++)
0233         {
0234             std::string histname = histvarnames[j].first;
0235             bool dohadd = (j == 0);
0236             TH1F *h = gethist<TH1F>(dohadd, infiledirs[i], histname);
0237 
0238             if (histvarnames[j].first == "hM_ClusADC_all" || histvarnames[j].first == "hM_ClusADC_all_etale0p1" || histvarnames[j].first == "hM_ClusADC_all_etale0p1_ClusZSize1")
0239             {
0240                 h->Rebin(20);
0241                 h->GetXaxis()->SetMaxDigits(3);
0242                 h->GetXaxis()->SetNdivisions(-9);
0243             }
0244             histmap[histvarnames[j]].push_back(h);
0245         }
0246     }
0247 
0248     for (auto varname : histvarnames)
0249     {
0250         if (varname.first == "hM_ClusZ_all" || varname.first == "hM_Eta_reco" || varname.first == "hM_cluseta_zvtxwei")
0251             draw_comparison(histmap[varname], false, varname.second, leg, Form("%s/%s_compdiffusionradius", plotdir.c_str(), varname.first.c_str()));
0252         else
0253             draw_comparison(histmap[varname], true, varname.second, leg, Form("%s/%s_compdiffusionradius", plotdir.c_str(), varname.first.c_str()));
0254     }
0255 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!