Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot_ntp_g4hit.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot_ntp_g4hit.C
1 #include "TCanvas.h"
2 #include "TH1D.h"
3 #include "TH2D.h"
4 #include "TCut.h"
5 #include "TChain.h"
6 #include "TLegend.h"
7 #include "TFile.h"
8 
9 //#include "sPhenixStyle.C"
10 #include "/phenix/u/yuhw/RootMacros/sPHENIXStyle/sPhenixStyle.C"
11 
12 void plot_ntp_g4hit(
13  //const char* input = "g4svtx_eval.root",
14  const char* input = "/sphenix/user/frawley/pileup/macros/macros/g4simulations/eval_output/g4svtx_eval_1.root_g4svtx_eval.root",
15  const int min_nfromtruth = 30
16  ){
17 
18  //SetsPhenixStyle();
19  gROOT->SetStyle("Plain");
20  gStyle->SetOptStat(0);
21  gStyle->SetOptFit(0);
22  gStyle->SetOptTitle(1);
23 
24  bool acts = false; // is this output from acts tracking?
25 
26  // OK for embedded particles only
27  //TCut good_gtrack_cut = Form("gtrackID>=0 && gntpc > 20");
28  //TCut good_track_cut = Form("gtrackID>=0 && ntpc > 20 && quality < 10");
29 
30  // For Hijing events with embedded particles
31  TCut good_gtrack_cut = Form("gtrackID>=0 && gembed >= 2 && gntpc > 20 && gnmaps > 1");
32  TCut good_track_cut = Form("gtrackID>=0 && gembed >= 2 && ntpc > 20 && quality < 10 && nmaps > 1");
33 
34  //TCut good_gtrack_cut = Form("gtrackID>=0 && gembed == 2 && gntpc > 20");
35  //TCut good_track_cut = Form("gtrackID>=0 && gembed == 2 && ntpc > 20 && quality < 10");
36 
37 
38  // rough 4 sigma cut
39  //TCut pt_cut = Form("fabs((pt - gpt)/pt) < 4*(0.02+0.0012*gpt)"); // 4 times sigma (where sigma = 2.2% at 2 and 4.4% at 20
40 
41  TFile *fout = new TFile("root_files/ntp_track_out.root","recreate");
42 
43  TChain* ntp_track = new TChain("ntp_track","reco tracks");
44  TChain* ntp_gtrack = new TChain("ntp_gtrack","g4 tracks");
45 
46  bool use_list = false;
47  int n_list = 1000;
48 
49  int ifile = 0;
50  int nbadvtx = 0;
51  for(int i=0;i<n_list;i++)
52  {
53  char name[500];
54  ifile = i;
55 
56  //sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
57  sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output_2/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
58  //sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output_3/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
59 
60  // Skip any files where the event vertex was not reconstructed properly
61  TChain* ntp_vertex = new TChain("ntp_vertex","events");
62  ntp_vertex->Add(name);
63  TH1D *hzvtx = new TH1D("hzvtx","hzvtx",500,-10.0,10.0);
64  ntp_vertex->Draw("vz - gvz >>hzvtx");
65  double meanzdiff = hzvtx->GetMean() ;
66  cout << " file " << i << " has vertex mean Z difference = " << meanzdiff << endl;
67  if( fabs(meanzdiff) > 0.1 )
68  {
69  cout << " --- Bad event vertex, skip file number " << i << " with name " << name << endl;
70  nbadvtx++;
71  delete ntp_vertex;
72  delete hzvtx;
73  continue;
74  }
75  delete ntp_vertex;
76  delete hzvtx;
77 
78  cout << "Adding file number " << i << " with name " << name << endl;
79 
80  ntp_gtrack->Add(name);
81  ntp_track->Add(name);
82  }
83  cout << " Found " << nbadvtx << " events with badly reconstructed vertices, and skipped them" << endl << endl;
84 
85  cout << " ntp_gtrack chain entries = " << ntp_gtrack->GetEntries() << endl;
86  cout << " ntp_track chain entries = " << ntp_track->GetEntries() << endl;
87 
88  /*
89  TH1D *hch = new TH1D("hch","hch",100,-2,+2);
90  TH1D *hgch = new TH1D("hgch","hgch",100,-2,+2);
91  ntp_track->Draw("charge>>hch",good_track_cut);
92  ntp_gtrack->Draw("charge>>hgch",good_gtrack_cut);
93  */
94 
95  TCanvas *cmvtx = new TCanvas("cmvtx","cmvtx");
96  TH1D *hnmvtx = new TH1D("hnmvtx","hnmvtx", 100, -1, 4);
97  ntp_track->Draw("nmaps>>hnmvtx", good_track_cut);
98  hnmvtx->Draw();
99  hnmvtx->Write();
100 
101  int nbinptres = 80; // number of bins in pT
102  int nbinpteff = 120; // number of bins in pT
103  int nbinptdca = 40; // number of bins in pT
104  int nbin = 400; // number of bins in Y slice
105  double range = 0.4;
106  double ptmax = 40.0;
107 
108  TCanvas *cpt = new TCanvas("cpt","cpt");
109  TH2D *h1 = new TH2D("h1","h1",nbinptres, 0.0 ,ptmax,nbin,-range,range);
110  ntp_track->Draw("(pt-gpt)/gpt:gpt>>h1",good_track_cut);
111  //h1->FitSlicesY(0,0,-1,0,"qnrl");
112  h1->FitSlicesY();
113  TH1D*h1_1 = (TH1D*)gDirectory->Get("h1_1");
114  TH1D*h1_2 = (TH1D*)gDirectory->Get("h1_2");
115  h1_2->Draw("e");
116  //h1_1->Draw("same");
117 
118  h1_1->SetMarkerStyle(4);
119  h1_1->SetMarkerColor(kBlack);
120  h1_1->SetLineColor(kBlack);
121 
122  h1_2->SetMarkerStyle(20);
123  h1_2->SetMarkerColor(kBlack);
124  h1_2->SetLineColor(kBlack);
125  h1_2->GetYaxis()->SetRangeUser(0, 0.10);
126  h1_2->GetYaxis()->SetTitleOffset(1.5);
127  h1_2->SetStats(0);
128  h1_2->SetTitle(";p_{T} [GeV/c];#frac{#Delta p_{T}}{p_{T}}");
129 
130  h1->Write();
131  h1_2->Write();
132 
133  TCanvas * c3 = new TCanvas("c3","c3");
134  TH1D* h3_den = new TH1D("h3_den","; p_{T}; Efficiency",nbinpteff, 0, ptmax);
135  TH1D* h3_num = (TH1D*)h3_den->Clone("h3_num");;
136  TH1D* h3_eff = (TH1D*)h3_den->Clone("h3_eff");;
137 
138  cout<<__LINE__<<": "<< good_track_cut <<endl;
139  ntp_gtrack->Draw("gpt>>h3_den",good_gtrack_cut);
140  //ntp_track->Draw("gpt>>h3_num",good_track_cut && pt_cut);
141  ntp_track->Draw("gpt>>h3_num",good_track_cut);
142 
143  for(int i=1;i<=h3_den->GetNbinsX();++i){
144  double pass = h3_num->GetBinContent(i);
145  double all = h3_den->GetBinContent(i);
146  double eff = 0;
147  if(all > pass)
148  eff = pass/all;
149  else if(all > 0)
150  eff = 1.;
151 
152  //double err = BinomialError(pass, all);
153  h3_eff->SetBinContent(i, eff);
154  //h3_eff->SetBinError(i, err);
155  }
156 
157  //h3_eff->Draw("e,text");
158  h3_eff->Draw("p");
159  h3_eff->SetStats(0);
160  h3_eff->SetTitle("; p_{T} [GeV/c]; eff.");
161  h3_eff->SetMarkerStyle(20);
162  h3_eff->SetMarkerColor(kBlack);
163  h3_eff->SetLineColor(kBlack);
164  h3_eff->GetYaxis()->SetRangeUser(0.0, 1.1);
165 
166  h3_den->Write();
167  h3_num->Write();
168  h3_eff->Write();
169 
170  TCanvas *c4 = new TCanvas("c4","c4",5,5,800,800);
171 
172 
173  TH2D *h2 = new TH2D("h2","h2",nbinptdca, 0,ptmax,nbin,-0.06,0.06);
174  ntp_track->Draw("(dca3dxy):gpt>>h2",good_track_cut);
175  h2->FitSlicesY();
176  TH1D*h2_1 = (TH1D*)gDirectory->Get("h2_1");
177  TH1D*h2_2 = (TH1D*)gDirectory->Get("h2_2");
178  h2_2->Draw("e");
179 
180  h2_1->SetMarkerStyle(4);
181  h2_1->SetMarkerColor(kBlack);
182  h2_1->SetLineColor(kBlack);
183 
184  h2_2->SetMarkerStyle(20);
185  h2_2->SetMarkerColor(kBlack);
186  h2_2->SetLineColor(kBlack);
187  //h2_2->GetYaxis()->SetRangeUser(0, 0.1);
188  h2_2->GetYaxis()->SetRangeUser(0.,0.01);
189  h2_2->GetYaxis()->SetTitleOffset(1.5);
190  h2_2->SetStats(0);
191  h2_2->SetTitle(";p_{T} [GeV/c];dca2d (cm)}");
192 
193  h2->Write();
194  h2_2->Write();
195 
196 
197  TCanvas *c5 = new TCanvas("c5","c5",5,5,800,800);
198 
199  TH2D *h3 = new TH2D("h3","h3",nbinptdca, 0, ptmax, nbin, -0.1, 0.1);
200  if(acts)
201  ntp_track->Draw("(dca3dz-gvz):gpt>>h3",good_track_cut);
202  else
203  ntp_track->Draw("(dca3dz):gpt>>h3",good_track_cut);
204  h3->FitSlicesY();
205  TH1D*h3_1 = (TH1D*)gDirectory->Get("h3_1");
206  TH1D*h3_2 = (TH1D*)gDirectory->Get("h3_2");
207  h3_2->Draw("e");
208  //h3_1->Draw("same");
209 
210  h3_1->SetMarkerStyle(4);
211  h3_1->SetMarkerColor(kBlack);
212  h3_1->SetLineColor(kBlack);
213 
214  h3_2->SetMarkerStyle(20);
215  h3_2->SetMarkerColor(kBlack);
216  h3_2->SetLineColor(kBlack);
217  //h3_2->GetYaxis()->SetRangeUser(0, 0.1);
218  h3_2->GetYaxis()->SetRangeUser(0.,0.01);
219  h3_2->GetYaxis()->SetTitleOffset(1.5);
220  h3_2->SetStats(0);
221  h3_2->SetTitle(";p_{T} [GeV/c];dca2d (cm)}");
222 
223  h3->Write();
224  h3_2->Write();
225 
226  TCanvas *c6 = new TCanvas("c6","c6", 5,5,1200, 800);
227 
228  TH2D *h6 = new TH2D("h6","h6",nbinptdca, 0, ptmax, nbin, -1.0, 4.0);
229  ntp_track->Draw("ntrumaps:gpt>>h6",good_track_cut);
230  h6->Draw();
231 
232  TH1D *h6_1 = new TH1D("h6_1","h6_1",nbin, 0.0, 1.2);
233  int binlo = h6->GetYaxis()->FindBin(2.5);
234  int binhi = h6->GetYaxis()->FindBin(3.5);
235  h6->ProjectionX("h6_1",binlo,binhi);
236  h6_1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
237  h6_1->SetLineWidth(2.0);
238  h6_1->SetLineColor(kBlue);
239 
240 
241  TH1D *h6_2 = new TH1D("h6_2","h6_2",nbin, 0.0, 1.2);
242  binlo = h6->GetYaxis()->FindBin(-0.5);
243  binhi = h6->GetYaxis()->FindBin(0.5);
244  h6->ProjectionX("h6_2",binlo,binhi);
245  h6_2->GetXaxis()->SetTitle("p_{T} (GeV/c)");
246  h6_2->SetLineColor(kRed);
247  h6_2->SetLineWidth(2.0);
248 
249 
250  h6_1->Draw();
251  h6_2->Draw("same");
252 
253  TLegend *l = new TLegend(0.15,0.75,0.45,0.85,"","NDC");
254  l->SetBorderSize(1);
255  l->AddEntry(h6_2,"ntrumaps = 0","l");
256  l->AddEntry(h6_1,"ntrumaps = 3","l");
257  l->Draw();
258 
259  h6->Write();
260  h6_1->Write();
261  h6_2->Write();
262 
263  /*
264  TCanvas *c7 = new TCanvas("c7","c7", 5,5,1200, 800);
265 
266  TH2D *h7 = new TH2D("h7","h7",nbin, 0, 1.1, nbin, -1.0, 50.0);
267  ntp_track->Draw("ntrutpc:gpt>>h7","ntrumaps == 0");
268  h7->Draw();
269  */
270  /*
271  TH1D *h7_1 = new TH1D("h7_1","h7_1",100, 0.0, 1.2);
272  h7->ProjectionX("h7_1");
273  h7_1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
274  h7_1->SetLineWidth(2.0);
275  h7_1->SetLineColor(kBlue);
276  h7_1->Draw();
277  */
278  /*
279  TH2D *h8 = new TH2D("h8","h8",100, 0, 1.1, 100, -1.0, 50.0);
280  ntp_track->Draw("ntrutpc:gpt>>h7","ntrumaps == 3");
281  h8->Draw();
282 
283  TH1D *h8_1 = new TH1D("h8_1","h8_1",100, 0.0, 1.2);
284  h8->ProjectionX("h8_1");
285  h8_1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
286  h8_1->SetLineColor(kRed);
287  h8_1->SetLineWidth(2.0);
288  h8_1->Draw("same");
289 
290  TLegend *l = new TLegend(0.15,0.75,0.45,0.85,"","NDC");
291  l->SetBorderSize(1);
292  l->AddEntry(h7_1,"ntrumaps = 0","l");
293  l->AddEntry(h8_1,"ntrumaps = 3","l");
294  l->Draw();
295  */
296 }