Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot.C
1 void NormByBinWidth(TH1F *inh) {
2  for(int i=0; i!=inh->GetXaxis()->GetNbins(); ++i) {
3  float bef = inh->GetBinContent( i+1 );
4  float err = inh->GetBinError( i+1 );
5  float bw = inh->GetXaxis()->GetBinWidth( i+1 );
6  inh->SetBinContent( i+1, bef/bw );
7  inh->SetBinError( i+1, err/bw );
8  }
9 }
10 
11 TH1F* getHist(TString nom, TFile *fi, const char *na, int co, int mk=20, float mks=1) {
12  TH1F *ori = (TH1F*) fi->Get(na);
13  TH1F *ret = (TH1F*) ori->Clone( nom.Data() );
14  ret->SetLineColor( co );
15  ret->SetMarkerColor( co );
16  ret->SetMarkerStyle( mk );
17  ret->SetMarkerSize( mks );
18  ret->Sumw2();
19  return ret;
20 }
21 
22 TH2F* getHist2(TString nom, TFile *fi, const char *na) {
23  TH2F *ori = (TH2F*) fi->Get(na);
24  TH2F *ret = (TH2F*) ori->Clone( nom.Data() );
25  ret->Sumw2();
26  return ret;
27 }
28 
29 TCanvas* MakeEfficiency(TFile *fa, TString ti, float c0=0.01250, float c1=0.001858, bool plotOld=false) {
30  TCanvas *cTN = new TCanvas();
31  cTN->SetGridx(1);
32  cTN->SetGridy(1);
33  TString geem = "Truths2_Pt";
34  TString gere = "Truths3_Pt";
35  TString re = "Tracks2_ResPt";
36  TH1F *fa_GEEM = getHist("EffGenEM", fa, geem.Data(), kBlue-3);
37  TH1F *fa_GERE = getHist("EffGenRE", fa, gere.Data(), kGreen-3);
38  TH1F *fa_RE = (TH1F*) fa_GEEM->Clone("EffRec");
39  fa_RE->Reset();
40  TH2F *fa_RE2D = getHist2("EffRec2D", fa, re.Data());
41  for(int b=1; b!=fa_RE->GetXaxis()->GetNbins()+1; ++b) {
42  float pt = fa_RE->GetXaxis()->GetBinCenter(b);
43  float sg = TMath::Sqrt(c0*c0+c1*pt*c1*pt);
44  int x = fa_RE2D->GetXaxis()->FindBin( pt );
45  int ya = fa_RE2D->GetYaxis()->FindBin( -5*sg );
46  int yb = fa_RE2D->GetYaxis()->FindBin( +5*sg );
47  float inte = fa_RE2D->Integral(x,x,ya,yb);
48  cout << x << " " << ya << " " << yb << " " << inte << endl;
49  fa_RE->SetBinContent(b,inte);
50  }
51  TGraphAsymmErrors *eff1 = new TGraphAsymmErrors();
52  TGraphAsymmErrors *eff2 = new TGraphAsymmErrors();
53  eff1->Divide( fa_RE, fa_GEEM );
54  eff2->Divide( fa_RE, fa_GERE );
55  eff1->SetLineColor( kBlue-3 );
56  eff1->SetMarkerColor( kBlue-3 );
57  eff1->SetMarkerStyle( 20 );
58  eff2->SetLineColor( kGreen-3 );
59  eff2->SetMarkerColor( kGreen-3 );
60  eff2->SetMarkerStyle( 20 );
61 
62  TH2F *axis = new TH2F("axis",";p_{T} [GeV];Efficiency",10,0.15,36.1,10,0,1.05);
63  axis->Draw();
64  axis->GetXaxis()->SetMoreLogLabels();
65  if(plotOld) eff1->Draw( "PSAME" );
66  eff2->Draw( "PSAME" );
67 
68  if(0) { //debug
69  fa_GE->Draw();
70  fa_RE->Draw("SAME");
71  TLegend *leg = new TLegend(0.6,0.6,0.9,0.9);
72  leg->AddEntry(fa_GE, Form("GEN %.0f",fa_GE->GetEntries()) );
73  leg->AddEntry(fa_RE, Form("REC %.0f",fa_RE->GetEntries()) );
74  leg->Draw();
75  }
76  return cTN;
77 }
78 
79 TCanvas* MakeEfficiencyFast(TFile *fa, TString ti, TString hw, bool plotOld=false) {
80  TCanvas *cTN = new TCanvas();
81  cTN->SetGridx(1);
82  cTN->SetGridy(1);
83  TString geem = Form("Truths2%s",hw.Data());
84  TString gere = Form("Truths3%s",hw.Data());
85  TString re = Form("Truths4%s",hw.Data());
86  TH1F *fa_RE = getHist("EffRec", fa, re.Data(), kOrange);
87  TH1F *fa_GEEM = getHist("EffGenEM", fa, geem.Data(), kRed-3);
88  TH1F *fa_GERE = getHist("EffGenRE", fa, gere.Data(), kBlue-3);
89  TGraphAsymmErrors *eff1 = new TGraphAsymmErrors();
90  TGraphAsymmErrors *eff2 = new TGraphAsymmErrors();
91  eff1->Divide( fa_RE, fa_GEEM );
92  eff2->Divide( fa_RE, fa_GERE );
93  eff1->SetLineColor( kBlue-3 );
94  eff1->SetMarkerColor( kBlue-3 );
95  eff1->SetMarkerStyle( 20 );
96  eff2->SetLineColor( kOrange-3 );
97  eff2->SetMarkerColor( kOrange-3 );
98  eff2->SetMarkerStyle( 20 );
99 
100  TH2F *axis = new TH2F("axis",";p_{T} [GeV];Efficiency",10,0.15,36.1,10,0,1.05);
101  axis->Draw();
102  axis->GetXaxis()->SetMoreLogLabels();
103  if(plotOld) eff1->Draw( "PSAME" );
104  eff2->Draw( "PSAME" );
105  return cTN;
106 }
107 
108 TCanvas* MakePtResolution(TFile *fa, TString ti, bool showBoth=false) {
109  TString re1 = "Tracks2_ResPt"; //dpt/pt
110  TString re2 = "Tracks2_Res2Pt"; //dpt/pt
111  TH2F *fa_RE1 = getHist2("EffRec1", fa, re1.Data());
112  TH2F *fa_RE2 = getHist2("EffRec2", fa, re2.Data());
113 
114  TCanvas *qa = new TCanvas();
115  qa->Divide(5,5);
116  TCanvas *qb = new TCanvas();
117  qb->Divide(5,5);
118  const int nbx1 = fa_RE1->GetXaxis()->GetNbins();
119  float x1[60], ex1[60];
120  float mu1[60], emu1[60], sg1[60], esg1[60];
121  float mu2[60], emu2[60], sg2[60], esg2[60];
122  for( int b=0; b!=nbx1; ++b) {
123  x1[b] = fa_RE1->GetXaxis()->GetBinCenter(b+1);
124  ex1[b] = fa_RE1->GetXaxis()->GetBinWidth(b+1)/2.0;
125  TH1D *tmp;
126  TF1 *fun;
127 
128  tmp = (TH1D*) fa_RE1->ProjectionY( Form("bin1_%d",b),b+1,b+1 );
129  if( tmp->GetEntries() < 30 ) continue;
130  if(b<25) qa->cd(b+1);
131  tmp->GetXaxis()->SetRangeUser(-0.2,+0.2);
132  tmp->Fit("gaus","RWQI","",-0.5,+0.5);
133  fun = (TF1*) tmp->GetListOfFunctions()->At(0);
134  mu1[b] = fun->GetParameter(1);
135  emu1[b] = fun->GetParError(1);
136  sg1[b] = fun->GetParameter(2);
137  esg1[b] = fun->GetParError(2);
138 
139  tmp = (TH1D*) fa_RE2->ProjectionY( Form("bin2_%d",b),b+1,b+1 );
140  if( tmp->GetEntries() < 30 ) continue;
141  if(b<25) qb->cd(b+1);
142  tmp->GetXaxis()->SetRangeUser(-0.2,+0.2);
143  tmp->Fit("gaus","RWQI","",-0.5,+0.5);
144  fun = (TF1*) tmp->GetListOfFunctions()->At(0);
145  mu2[b] = fun->GetParameter(1)/x1[b];
146  emu2[b] = fun->GetParError(1)/x1[b];
147  sg2[b] = fun->GetParameter(2)/x1[b];
148  esg2[b] = fun->GetParError(2)/x1[b];
149  //cout << mu1[b] << " " << emu1[b] << " || " << sg1[b] << " " << esg1[b] << endl;
150  }
151 
152  TCanvas *cTN = new TCanvas();
153  cTN->SetGridx(1);
154  cTN->SetGridy(1);
155 
156  TH2F *axis = new TH2F("axisPtRes",";p_{T} [GeV]",10,0.15,36.1,10,-0.01,0.09);
157  axis->Draw();
158 
159  TGraphErrors *gRes1 = new TGraphErrors(nbx1,x1,sg1,ex1,esg1);
160  gRes1->Draw("PSAME");
161  gRes1->SetLineColor(kRed-3);
162  gRes1->SetMarkerColor(kRed-3);
163  gRes1->SetMarkerStyle(20);
164 
165  TGraphErrors *gRes1M = new TGraphErrors(nbx1,x1,mu1,ex1,emu1);
166  gRes1M->Draw("PSAME");
167  gRes1M->SetLineColor(kRed-3);
168  gRes1M->SetMarkerColor(kRed-3);
169  gRes1M->SetMarkerStyle(24);
170 
171  TGraphErrors *gRes2 = new TGraphErrors(20,x1,sg2,ex1,esg2);
172  if(showBoth) gRes2->Draw("PSAME");
173  gRes2->SetLineColor(kBlue-3);
174  gRes2->SetMarkerColor(kBlue-3);
175  gRes2->SetMarkerStyle(20);
176 
177  TGraphErrors *gRes2M = new TGraphErrors(20,x1,mu2,ex1,emu2);
178  if(showBoth) gRes2M->Draw("PSAME");
179  gRes2M->SetLineColor(kBlue-3);
180  gRes2M->SetMarkerColor(kBlue-3);
181  gRes2M->SetMarkerStyle(24);
182 
183  TF1 *fun = new TF1("fun","TMath::Sqrt([0]*[0]+[1]*x*[1]*x)",0.2,36);
184  //fun->SetParameter(0,0.01); fun->SetParLimits(0,0.001,0.02);
185  fun->SetParameter(0,1.23e-2); fun->SetParLimits(0,+1,-1);
186  fun->SetParameter(1,0.01); fun->SetParLimits(1,0.001,0.02);
187  fun->SetLineColor(kOrange-3);
188  gRes1->Fit(fun,"R","",0.2,36);
189 
190  gRes1->SetFillColor(kWhite);
191  gRes2->SetFillColor(kWhite);
192  gRes1M->SetFillColor(kWhite);
193  gRes2M->SetFillColor(kWhite);
194 
195  TLegend *leg = new TLegend(0.1,0.7,0.4,0.9);
196  leg->AddEntry(gRes1,"p_{T} Res #sigma_{1}");
197  if(showBoth) leg->AddEntry(gRes2,"p_{T} Res #sigma_{0}");
198  leg->AddEntry(gRes1M,"p_{T} Res #mu_{1}");
199  if(showBoth) leg->AddEntry(gRes2M,"p_{T} Res #mu_{0}");
200  leg->Draw();
201  return cTN;
202 }
203 
204 TCanvas* MakeDCA2DResolution(TFile *fa, TString ti) {
205  TString re1 = "Tracks2_Dca2D";
206  TH2F *fa_RE1 = getHist2("Dca2DRec1", fa, re1.Data());
207 
208  TCanvas *qa = new TCanvas();
209  qa->Divide(5,5);
210  const int nbx1 = fa_RE1->GetXaxis()->GetNbins();
211  float x1[60], ex1[60];
212  float mu1[60], emu1[60], sg1[60], esg1[60];
213  for( int b=0; b!=nbx1; ++b) {
214  x1[b] = fa_RE1->GetXaxis()->GetBinCenter(b+1);
215  ex1[b] = fa_RE1->GetXaxis()->GetBinWidth(b+1)/2.0;
216  TH1D *tmp;
217  TF1 *fun;
218 
219  tmp = (TH1D*) fa_RE1->ProjectionY( Form("dca2d_bin1_%d",b),b+1,b+1 );
220  if( tmp->GetEntries() < 30 ) continue;
221  if(b<25) qa->cd(b+1);
222  tmp->Fit("gaus","RWQI","",-0.5,+0.5);
223  fun = (TF1*) tmp->GetListOfFunctions()->At(0);
224  mu1[b] = fun->GetParameter(1)*1e4;
225  emu1[b] = fun->GetParError(1)*1e4;
226  sg1[b] = fun->GetParameter(2)*1e4;
227  esg1[b] = fun->GetParError(2)*1e4;
228  }
229 
230  TCanvas *cTN = new TCanvas();
231  cTN->SetGridx(1);
232  cTN->SetGridy(1);
233 
234  TH2F *axis = new TH2F("axisDCA2DRes",";p_{T}^{reco} [GeV];[#mu m]",10,0.15,36.1,10,-10,+80);
235  axis->Draw();
236 
237  TGraphErrors *gRes1 = new TGraphErrors(nbx1,x1,sg1,ex1,esg1);
238  gRes1->Draw("PSAME");
239  gRes1->SetLineColor(kRed-3);
240  gRes1->SetMarkerColor(kRed-3);
241  gRes1->SetMarkerStyle(20);
242 
243  TGraphErrors *gRes1M = new TGraphErrors(nbx1,x1,mu1,ex1,emu1);
244  gRes1M->Draw("PSAME");
245  gRes1M->SetLineColor(kRed-3);
246  gRes1M->SetMarkerColor(kRed-3);
247  gRes1M->SetMarkerStyle(24);
248 
249  gRes1->SetFillColor(kWhite);
250  gRes1M->SetFillColor(kWhite);
251 
252  TLegend *leg = new TLegend(0.55,0.7,0.9,0.9);
253  leg->AddEntry(gRes1,"DCA2D Gaus #sigma");
254  leg->AddEntry(gRes1M,"DCA2D Gaus #mu");
255  leg->Draw();
256  return cTN;
257 }
258 
259 TCanvas* MakeTH1F(TFile *fa,TString ti,TString wh,TString hw, bool norm=true, bool bw=false, bool log=false) {
260  TCanvas *cTN = new TCanvas();
261  TString c0 = Form("%s0%s",wh.Data(),hw.Data());
262  TString c1 = Form("%s1%s",wh.Data(),hw.Data());
263  TString c2 = Form("%s2%s",wh.Data(),hw.Data());
264  TString c3 = Form("%s3%s",wh.Data(),hw.Data());
265  TH1F *fa_T0N = getHist(Form("%s0",ti.Data()), fa, c0.Data(), kBlack, 20, 1.3);
266  TH1F *fa_T1N = getHist(Form("%s1",ti.Data()), fa, c1.Data(), kRed-3, 20, 1.3);
267  TH1F *fa_T2N = getHist(Form("%s2",ti.Data()), fa, c2.Data(), kBlue-3, 20, 1.3);
268  TH1F *fa_T3N = getHist(Form("%s3",ti.Data()), fa, c3.Data(), kGreen-3, 20, 1.3);
269  if(norm) {
270  int evA = ((TH1F*) fa->Get("Events"))->GetBinContent(1);
271  fa_T0N->Scale(1.0/evA);
272  fa_T1N->Scale(1.0/evA);
273  fa_T2N->Scale(1.0/evA);
274  fa_T3N->Scale(1.0/evA);
275  }
276  if(bw) {
277  NormByBinWidth(fa_T0N);
278  NormByBinWidth(fa_T1N);
279  NormByBinWidth(fa_T2N);
280  NormByBinWidth(fa_T3N);
281  }
282  if(log) cTN->SetLogy(1);
283  fa_T0N->Draw();
284  fa_T1N->Draw("SAME");
285  fa_T2N->Draw("SAME");
286  fa_T3N->Draw("SAME");
287  TLegend *leg = new TLegend(0.6,0.6,0.9,0.9);
288  leg->AddEntry(fa_T0N,fa_T0N->GetTitle());
289  leg->AddEntry(fa_T1N,fa_T1N->GetTitle());
290  leg->AddEntry(fa_T2N,fa_T2N->GetTitle());
291  leg->AddEntry(fa_T3N,fa_T3N->GetTitle());
292  leg->Draw();
293  fa_T0N->SetTitle( ti.Data() );
294  return cTN;
295 }
296 
297 void plot() {
298  gStyle->SetOptStat(0);
299  //TString file = "newTPCgeo";
300  TString file = "oldUpdate";
301 
302  TFile *fa = new TFile( Form("%s.root",file.Data()) );
303 
304  fa->ls();
305  MakeTH1F(fa,"Eta","Truths","_Eta",true,false,true)->Draw();
306  MakeTH1F(fa,"Pt [GeV]","Truths","_Pt",true,true,true)->Draw();
307 
308  MakeTH1F(fa,"Eta","Tracks","_Eta",true,false,true)->Draw();
309  MakeTH1F(fa,"Pt [GeV]","Tracks","_Pt",true,true,true)->Draw();
310 
311  TCanvas *cres = MakePtResolution(fa,"PtResolution");
312  TCanvas *cdca = MakeDCA2DResolution(fa,"DCA2DResolution");
313  //TCanvas *cef1 = MakeEfficiency(fa,"Efficiency",1.23e-2,2.17e-03);
314  TCanvas *cef1 = MakeEfficiency(fa,"Efficiency",1.23e-2,1.85e-3);
315  TCanvas *cef2 = MakeEfficiencyFast(fa,"Efficiency","_Pt");
316 
317  cres->SaveAs( Form("%s_cres.png",file.Data()), "png"); cres->Draw();
318  cdca->SaveAs( Form("%s_cdca.png",file.Data()), "png"); cdca->Draw();
319  cef1->SaveAs( Form("%s_cef1.png",file.Data()), "png"); cef1->Draw();
320  cef2->SaveAs( Form("%s_cef2.png",file.Data()), "png"); cef2->Draw();
321 }