Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cluster_efficiency.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file cluster_efficiency.C
1 #include "TChain.h"
2 #include "TString.h"
3 #include "TH1D.h"
4 #include "TH2D.h"
5 #include "TCanvas.h"
6 #include "TROOT.h"
7 #include "TStyle.h"
8 #include "TLine.h"
9 #include "TProfile.h"
10 #include "TCut.h"
11 #include "TLegend.h"
12 #include "TLatex.h"
13 
14 #include <iostream>
15 
16 void cluster_efficiency()
17 {
18  gROOT->SetStyle("Plain");
19  gStyle->SetOptStat(0);
20  gStyle->SetOptTitle(1);
21 
22  TCut good_g4cluster_cut = Form("layer >= 0 && gr > 0 && gembed > 0 && gtrackID > 0");
23  TCut good_cluster_cut = Form("layer >= 0 && r > 0 && gembed > 0 && gtrackID > 0");
24  TCut layer_cut = Form("layer > 6");
25  TCut layer_cut_tpc[3];
26  layer_cut_tpc[0] = Form("layer > 6 && layer <= 22");
27  layer_cut_tpc[1] = Form("layer > 22 && layer <= 38");
28  layer_cut_tpc[2] = Form("layer > 38 && layer <= 54");
29  TCut layer_cut_intt = Form("layer > 2 && layer < 7");
30  TCut layer_cut_mvtx = Form("layer < 3");
31  TCut etacut = Form("geta > -1.0");
32  //TCut etacut = Form("geta < 0.2 && geta > 0.1");
33  //TCut etacut = Form("geta < 1.0 && geta > -1.0");
34  //TCut ng4hitscut = Form("ng4hits > 0");
35  TCut ng4hitscut = Form("");
36 
37  TChain* ntp_g4cluster = new TChain("ntp_g4cluster","truth clusters");
38  TChain* ntp_cluster = new TChain("ntp_cluster","reco clusters");
39 
40  for(int ifile = 0; ifile < 100; ++ifile)
41  {
42  char name[500];
43 
44  //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/repo_100_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
45  //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/2k_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
46  //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/4k_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
47  //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/100muon_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
48 
49  sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/100_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
50  //sprintf(name,"/sphenix/user/frawley/current_repo/macros/macros/g4simulations/genfit_pion_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
51  //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
52 
53  //sprintf(name,"/sphenix/user/frawley/cluster_efficiency/macros/macros/g4simulations/low_occupancy_fiducial0.5_search2_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
54 
55  std::cout << "Open file " << name << std::endl;
56 
57  ntp_g4cluster->Add(name);
58  ntp_cluster->Add(name);
59  }
60 
61  // cluster efficiency vs layer
62 
63  TCanvas * ceff = new TCanvas("ceff","ceff",5,5,1200,800);
64  TH1D* h3_den = new TH1D("h3_den","; layer; Efficiency",54, 0, 54);
65  TH1D* h3_num = (TH1D*)h3_den->Clone("h3_num");;
66  TH1D* h3_eff = (TH1D*)h3_den->Clone("h3_eff");;
67 
68  ntp_g4cluster->Draw("layer>>h3_den",good_g4cluster_cut);
69  ntp_g4cluster->Draw("layer>>h3_num",good_cluster_cut);
70 
71  for(int i=1;i<=h3_den->GetNbinsX();++i){
72  double pass = h3_num->GetBinContent(i);
73  double all = h3_den->GetBinContent(i);
74  double eff = 0;
75  if(all > pass)
76  eff = pass/all;
77  else if(all > 0)
78  eff = 1.;
79  //std::cout << " i " << i << " pass " << pass << " all " << all << " eff " << eff << std::endl;
80  //double err = BinomialError(pass, all);
81  h3_eff->SetBinContent(i, eff);
82  //h3_eff->SetBinError(i, err);
83  }
84 
85  //h3_eff->Draw("e,text");
86  h3_eff->Draw("p");
87  h3_eff->SetStats(0);
88  h3_eff->SetTitle("truth clusters with matched reco cluster; layer; cluster efficiency");
89  h3_eff->SetMarkerStyle(20);
90  h3_eff->SetMarkerColor(kBlack);
91  h3_eff->SetLineColor(kBlack);
92  //h3_eff->GetYaxis()->SetRangeUser(0.0, 1.1);
93  h3_eff->GetYaxis()->SetRangeUser(0.0, 1.04);
94 
95  TLine *unit = new TLine(0,1.0,54,1.0);
96  unit->SetLineStyle(2);
97  unit->Draw();
98 
99  TCanvas *crestpc = new TCanvas("crestpc","crestpc",5,5,1600,800);
100  crestpc->Divide(3,1);
101  TH1D *hreslayer[3];
102  for(int i=0;i<3;++i)
103  {
104  crestpc->cd(i+1);
105  gPad->SetLeftMargin(0.15);
106  char name[500];
107  sprintf(name,"hreslayer%i",i);
108  if(i==0)
109  hreslayer[i] = new TH1D(name,"cluster r-phi resolution TPC inner",300,-0.095, 0.095);
110  else if(i==1)
111  hreslayer[i] = new TH1D(name,"cluster r-phi resolution TPC mid",300,-0.095, 0.095);
112  else
113  hreslayer[i] = new TH1D(name,"Cluster r-phi resolution TPC outer",300,-0.095, 0.095);
114  char drawit[500];
115  sprintf(drawit,"gr*(phi-gphi)>>hreslayer%i",i);
116  ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_tpc[i]);
117  hreslayer[i]->GetXaxis()->SetTitle("gr*(phi-gphi) (cm)");
118  hreslayer[i]->GetXaxis()->SetNdivisions(504);
119  hreslayer[i]->Draw("l");
120 
121  cout << " RMS for TPC r-phi for i = " << i << " is " << hreslayer[i]->GetRMS() << endl;;
122 
123  char latstr[500];
124  sprintf(latstr,"RMS = %.0f (#mum)", 1.0e04 * hreslayer[i]->GetRMS());
125  TLatex *lat = new TLatex(0.6, 0.915, latstr);
126  lat->SetNDC();
127  lat->Draw();
128  }
129 
130  TCanvas *crestpcz = new TCanvas("crestpcz","crestpcz",5,5,1600,800);
131  crestpcz->Divide(3,1);
132  TH1D *hreslayerz[3];
133  for(int i=0;i<3;++i)
134  {
135  crestpcz->cd(i+1);
136  gPad->SetLeftMargin(0.15);
137  char name[500];
138  sprintf(name,"hreslayerz%i",i);
139  if(i==0)
140  hreslayerz[i] = new TH1D(name,"Cluster z resolution TPC inner",300,-0.4, 0.4);
141  else if(i==1)
142  hreslayerz[i] = new TH1D(name,"Cluster z resolution TPC mid",300,-0.4, 0.4);
143  else
144  hreslayerz[i] = new TH1D(name,"Cluster z resolution TPC outer",300,-0.4, 0.4);
145  char drawit[500];
146  sprintf(drawit,"z-gz>>hreslayerz%i",i);
147  ntp_g4cluster->Draw(drawit, good_cluster_cut && layer_cut_tpc[i] && etacut && ng4hitscut);
148  hreslayerz[i]->GetXaxis()->SetTitle("z-gz (cm)");
149  hreslayerz[i]->GetXaxis()->SetNdivisions(504);
150  hreslayerz[i]->Draw("l");
151 
152  char latstr[500];
153  sprintf(latstr,"RMS = %.0f (#mum)", 1.0e04 * hreslayerz[i]->GetRMS());
154  TLatex *lat = new TLatex(0.6, 0.915, latstr);
155  lat->SetNDC();
156  lat->Draw();
157  }
158 
159  TCanvas *cressil = new TCanvas("cressil","cressil",5,5,1600,800);
160  cressil->Divide(2,1);
161  TH1D *hreslayer_sil[2];
162  for(int i=0;i<2;++i)
163  {
164  cressil->cd(i+1);
165  gPad->SetLeftMargin(0.15);
166  char name[500];
167  sprintf(name,"hreslayer_sil%i",i);
168  if(i==0)
169  hreslayer_sil[i] = new TH1D(name,"Cluster resolution INTT",300,-0.0049, 0.0049);
170  else
171  hreslayer_sil[i] = new TH1D(name,"Cluster resolution MVTX",300,-0.0025, 0.0025);
172  char drawit[500];
173  sprintf(drawit,"gr*(phi-gphi)>>hreslayer_sil%i",i);
174  if(i==0)
175  ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_intt);
176  else
177  ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_mvtx);
178  hreslayer_sil[i]->GetXaxis()->SetTitle("gr*(phi-gphi) (cm)");
179  hreslayer_sil[i]->GetXaxis()->SetNdivisions(504);
180  hreslayer_sil[i]->Draw("l");
181 
182  char latstr[500];
183  sprintf(latstr,"RMS = %.1f (#mum)", 1.0e04 * hreslayer_sil[i]->GetRMS());
184  TLatex *lat = new TLatex(0.6, 0.915, latstr);
185  lat->SetNDC();
186  lat->Draw();
187  }
188 
189  TCanvas *cressilz = new TCanvas("cressilz","cressilz",5,5,1600,800);
190  cressilz->Divide(2,1);
191  TH1D *hreslayer_silz[2];
192  for(int i=0;i<2;++i)
193  {
194  cressilz->cd(i+1);
195  gPad->SetLeftMargin(0.15);
196  char name[500];
197  sprintf(name,"hreslayer_silz%i",i);
198  if(i==0)
199  hreslayer_silz[i] = new TH1D(name,"Cluster resolution INTT z",300,-1.0, 1.0);
200  else
201  hreslayer_silz[i] = new TH1D(name,"Cluster resolution MVTX z",300,-0.0025, 0.0025);
202  char drawit[500];
203  sprintf(drawit,"z-gz>>hreslayer_silz%i",i);
204  if(i==0)
205  ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_intt);
206  else
207  ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_mvtx);
208  hreslayer_silz[i]->GetXaxis()->SetTitle("z-gz (cm)");
209  hreslayer_silz[i]->GetXaxis()->SetNdivisions(504);
210  hreslayer_silz[i]->Draw("l");
211 
212  char latstr[500];
213  if(i==0)
214  sprintf(latstr,"RMS = %.0f (#mum)", 1.0e04 * hreslayer_silz[i]->GetRMS());
215  else
216  sprintf(latstr,"RMS = %.1f (#mum)", 1.0e04 * hreslayer_silz[i]->GetRMS());
217  TLatex *lat = new TLatex(0.6, 0.915, latstr);
218  lat->SetNDC();
219  lat->Draw();
220  }
221 
222  TCanvas *ctruth = new TCanvas("ctruth","ctruth",5,5,1200,800);
223  ctruth->SetLeftMargin(0.15);
224  TH1D *htruth = new TH1D("htruth","Truth clusters per layer",540,0,60);
225  ntp_g4cluster->Draw("layer>>htruth", good_g4cluster_cut);
226  htruth->SetMarkerStyle(20);
227  htruth->SetMarkerSize(2);
228  htruth->GetXaxis()->SetTitle("layer");
229  htruth->GetYaxis()->SetTitle("truth clusters / layer");
230  htruth->GetYaxis()->SetTitleOffset(1.7);
231  htruth->Draw("p");
232 
233  TCanvas *cphi = new TCanvas("cphi","cphi", 5,5,800,800);
234  TH1D *hphi = new TH1D("hphi","hphi",10000, -4.0, 4.0);
235  ntp_g4cluster->Draw("phi>>hphi", good_cluster_cut && layer_cut);
236  hphi->Draw();
237 
238  TCanvas *ceta = new TCanvas("ceta","ceta", 5,5,800,800);
239  TH1D *heta = new TH1D("heta","heta",10000, -1.2, 1.2);
240  ntp_g4cluster->Draw("eta>>heta", good_cluster_cut && layer_cut);
241  heta->Draw();
242 
243  TCanvas * cetaz = new TCanvas("cetaz","cetaz",5,5,800,800);
244  TH2D *hetaz = new TH2D("hetaz","hetaz",200, -1.2, 1.2, 200, -0.5, 0.5);
245  ntp_g4cluster->Draw("(z-gz):eta>>hetaz", good_cluster_cut && layer_cut);
246  hetaz->Draw("colz");
247 
248  /*
249  TCanvas * cetazsize = new TCanvas("cetazsize","cetazsize",5,5,800,800);
250  TH2D *hetazsize = new TH2D("hetazsize","hetazsize",200, -1.2, 1.2, 200, 0, 5.0);
251  ntp_g4cluster->Draw("zsize:eta>>hetazsize", good_cluster_cut && layer_cut);
252  hetazsize->Draw("colz");
253  */
254 
255 
256  // zsize
257 
258  TCanvas *c2 = new TCanvas("c2","c2", 5,5,800,800);
259 
260  //TH2D *h21 = new TH2D("h21","h21",800, -110, 110, 5000, 0, 8);
261  TH2D *h21 = new TH2D("h21","h21",800, -1.2, 1.2, 5000, 0, 8);
262  //ntp_g4cluster->Draw("zsize*.98:z>>h21",good_cluster_cut);
263  ntp_g4cluster->Draw("zsize*.98:geta>>h21",good_cluster_cut);
264  h21->SetMarkerStyle(20);
265  h21->SetMarkerSize(0.5);
266  //h21->GetXaxis()->SetTitle("Z (cm)");
267  h21->GetXaxis()->SetTitle("#eta");
268  h21->GetYaxis()->SetTitle("zsize (cm)");
269 
270  //TH2D *h22 = new TH2D("h22","h22",800, -110, 110, 5000, 0, 8);
271  TH2D *h22 = new TH2D("h22","h22",800, -1.2, 1.2, 5000, 0, 8);
272  //ntp_g4cluster->Draw("gzsize:gz>>h22",good_g4cluster_cut);
273  ntp_g4cluster->Draw("gzsize:geta>>h22",good_g4cluster_cut);
274  h22->SetMarkerStyle(20);
275  h22->SetMarkerSize(0.5);
276  h22->SetMarkerColor(kRed);
277  //h21->GetXaxis()->SetTitle("Z (cm)");
278  h21->GetXaxis()->SetTitle("#eta");
279  h21->GetYaxis()->SetTitle("zsize (cm)");
280 
281  h21->Draw();
282  h22->Draw("same");
283 
284  // phisize
285 
286  TCanvas *c4 = new TCanvas("c4","c4", 5,5,800,800);
287 
288  TH2D *h41 = new TH2D("h41","h41",500, 0, 80.0, 5000, 0, 2.);
289  ntp_g4cluster->Draw("phisize*0.98:r>>h41",good_cluster_cut);
290 
291  h41->SetMarkerStyle(20);
292  h41->SetMarkerSize(0.5);
293  h41->GetXaxis()->SetTitle("radius (cm)");
294  h41->GetYaxis()->SetTitle("radius*phisize (cm)");
295 
296  TH2D *h42 = new TH2D("h42","h42",500, 0, 80.0, 5000, 0, 2.);
297  ntp_g4cluster->Draw("gphisize:gr>>h42",good_g4cluster_cut);
298  h42->SetMarkerStyle(20);
299  h42->SetMarkerSize(0.5);
300  h42->SetMarkerColor(kRed);
301  h42->GetXaxis()->SetTitle("radius (cm)");
302  h42->GetYaxis()->SetTitle("radius*phisize (cm)");
303 
304 
305  h41->Draw();
306  h42->Draw("same");
307 
308 
309  //h42->Draw();
310  //h41->Draw("same");
311 
312 }