Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cluster_resolution.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file cluster_resolution.C
1 
2 #include <TChain.h>
3 #include <TCanvas.h>
4 #include <TH1.h>
5 #include <TH1F.h>
6 #include <TH2D.h>
7 #include <TF1.h>
8 #include <TFile.h>
9 #include <TGraph.h>
10 #include <TStyle.h>
11 #include <TROOT.h>
12 #include <TLegend.h>
13 #include <TLine.h>
14 #include <TLatex.h>
15 #include <TRandom1.h>
16 #include <TPolyLine.h>
17 #include <TColor.h>
18 #include <iostream>
19 #include <map>
20 #include <fstream>
21 #include <TMath.h>
22 #include <TLorentzVector.h>
23 #include <TEfficiency.h>
24 #include <TGraphAsymmErrors.h>
25 #include <string.h>
26 #include <list>
27 
28 
29 using namespace std;
30 
31 //#include "sPhenixStyle.C"
32 
33 void cluster_resolution()
34 {
35 
36  gROOT->SetStyle("Plain");
37  gStyle->SetOptStat(0);//(0)
38  gStyle->SetOptFit(0);
39  gStyle->SetOptTitle(1); //(0)
40  gStyle->SetStatW(0.3);
41  gStyle->SetStatH(0.3);
42 
43  bool verbose = false;
44  int embed_flag = 2; // embedded pions
45 
46  // Setup parameters
47  //=================
48 
49  // These should match the setup that was simulated!
50  int n_maps_layers = 3;
51  int n_intt_layers = 6;
52 
53  int n_tpc_layers_inner = 16;
54  int n_tpc_layers_mid = 16;
55  int n_tpc_layers_outer = 16;
56  int n_tpc_layers = n_tpc_layers_inner + n_tpc_layers_mid + n_tpc_layers_outer;;
57  // This should be inner maps layers (3) + intt(4) + TPC (48)
58  int nlayers = n_maps_layers+n_intt_layers+n_tpc_layers;
59 
60  double ptmax = 35.0;
61 
62  int a = 0;
63  int b = 0;
64  int c = 0;
65 
66  //===============================
67 
68  // Open the evaluator output file
69  cout << "Reading ntuple " << endl;
70 
71  TH2D *rphi = new TH2D("rphi","Cluster map",2000, -79.0, 79.0, 2000, -79.0, 79.0);
72  rphi->GetYaxis()->SetTitle("Y (cm)");
73  rphi->GetXaxis()->SetTitle("X (cm)");
74 
75  TH2D *zphi = new TH2D("zphi","Cluster z-phi map",120, -120.0, 120.0, 500, -2.0, 5.0);
76  zphi->GetYaxis()->SetTitle("#phi");
77  zphi->GetXaxis()->SetTitle("Z (cm)");
78 
79  TH2D *delta_rphi = new TH2D("delta_rphi","cluster r#phi errors by layer",60.0, 0.0, 60.0, 2000, -0.1, 0.1);
80  delta_rphi->GetYaxis()->SetTitle("Cluster r#phi Error (cm)");
81  delta_rphi->GetXaxis()->SetTitle("Tracking Layer");
82 
83  TH2D *delta_phi = new TH2D("delta_phi","cluster #phi errors by layer",60.0, 0.0, 60.0, 2000, -0.002, 0.002);
84  delta_phi->GetYaxis()->SetTitle("Cluster #phi Error (rad)");
85  delta_phi->GetXaxis()->SetTitle("Tracking Layer");
86 
87  TH2D *delta_phi_gphi = new TH2D("delta_phi_gphi","cluster #phi errors by gphi",2000.0, -0.1, 0.1, 2000, -0.002, 0.002);
88  delta_phi_gphi->GetYaxis()->SetTitle("Cluster #phi Error (rad)");
89  delta_phi_gphi->GetXaxis()->SetTitle("Truth #phi");
90 
91  TH2D *delta_z = new TH2D("delta_z","cluster Z errors by layer",60.0, 0.0, 60.0, 20000, -1.0, 1.0);
92  delta_z->GetYaxis()->SetTitle("Cluster Z Error (cm)");
93  delta_z->GetXaxis()->SetTitle("Tracking Layer");
94 
95  TH2D *cluster_size = new TH2D("cluster_size","cluster size by layer",60.0, 0.0, 60.0, 200, 0.0, 15.0);
96  cluster_size->GetYaxis()->SetTitle("Cluster Size (hits)");
97  cluster_size->GetXaxis()->SetTitle("Tracking Layer");
98 
99  TH2D *clusters_layer = new TH2D("clusters_layer","cluster eta by layer",50.0, 0.0, 50.0, 200, -3.0, 3.0);
100  clusters_layer->GetYaxis()->SetTitle("Cluster eta");
101  clusters_layer->GetXaxis()->SetTitle("Tracking Layer");
102 
103  TH2D *hzsize = new TH2D("hzsize","cluster Z size vs nparticles",100, 0.0, 6.0, 200, 0.0, 7.0);
104  hzsize->GetXaxis()->SetTitle("Cluster Z Size (cm)");
105  hzsize->GetYaxis()->SetTitle("nparticles");
106 
107  int pick_layer = 7;
108 
109  double nclusters = 0;
110 
111  // The condor job output files
112  for(int i=0;i <100; i++)
113  {
114  // Declare all these so that the include file does not cause complaints
115  TChain *ntp_vertex = new TChain("ntp_vertex","clusters");
116  TChain *ntp_cluster = new TChain("ntp_cluster","clusters");
117  TChain *ntp_gtrack = new TChain("ntp_gtrack","clusters");
118  TChain *ntp_track = new TChain("ntp_track","clusters");
119  TChain *ntp_g4hit = new TChain("ntp_g4hit","clusters");
120 
121  char name[500];
122 
123  // input files
124 
125  sprintf(name,"/sphenix/user/frawley/cluster_efficiency/macros/macros/g4simulations/screwed_up_tracks_out/g4svtx_eval_%i.root_g4svtx_eval.root",i);
126  //sprintf(name,"/sphenix/user/frawley/cluster_efficiency/macros/macros/g4simulations/test_fudge_06_075_2kevts_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
127  //sprintf(name,"/sphenix/user/frawley/cluster_efficiency/macros/macros/g4simulations/test_fudge_12_15_2kevts_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
128  //sprintf(name,"/sphenix/user/frawley/cluster_efficiency/macros/macros/g4simulations/test_fudge_085_105_2kevts_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
129  cout << "Enter file loop with name = " << name << endl;
130 
131  // This include file contains the definitions of the ntuple variables
132 #include "ntuple_variables.C"
133 
134 
135  // eliminate events with bad vertices here if needed
136  ntp_vertex->Add(name);
137  int toss_event = 0;
138  for(int cl=0;cl<ntp_vertex->GetEntries();cl++)
139  {
140  ntp_vertex->GetEntry(cl);
141  if( fabs(evz-egvz) > 0.1)
142  toss_event++;
143  }
144  if(toss_event > 0)
145  {
146  cout << " -- bad reco event vertex in one of the events, skip this file. vz = " << evz << " gvz = " << egvz << endl;
147  delete ntp_vertex;
148  continue;
149  }
150 
151 
152  ntp_cluster->Add(name);
153 
154 
155  cout << " ntp_cluster entries = " << ntp_cluster->GetEntries() << endl;
156 
157 
158  for(int p=0;p < ntp_cluster->GetEntries(); p++)
159  {
160  ntp_cluster->GetEntry(p);
161 
162  if(cgembed != embed_flag)
163  continue;
164 
165  if(gprimary != 1)
166  continue;
167 
168  tphi = atan(gy/gx);
169  r = sqrt(gx*gx+gy*gy);
170 
171  double dphi = trphi - tphi;
172  double drphi = r * dphi;
173 
174  // Extract the cluster Z resolution
175  double dz = z - gz;
176 
177  //if(drphi > -0.1 && drphi < 0.1)
178  {
179  delta_rphi->Fill( (double) layer, drphi);
180  delta_phi->Fill( (double) layer, dphi);
181  delta_z->Fill( (double) layer, dz);
182  if(layer > 38) delta_phi_gphi->Fill(tphi,dphi);
183  }
184  }
185 
186  delete ntp_vertex;
187  delete ntp_cluster;
188 
189  }// for i
190 
191  char label[500];
192 
193  // for fitting the r-phi cluster distributions
194  TF1 *fg = new TF1("fg","gaus(0)",-0.05,0.05); // restrict range of fit to this
195  fg->SetLineColor(kRed);
196  fg->SetParameter(0, 100.0);
197  fg->SetParameter(1, 0.0);
198  fg->SetParameter(2, 2e-02);
199  //fg->SetParLimits(2, 0.0, 5e-02);
200 
201  TCanvas *c7 = new TCanvas("c7","c7",50,50,1200,800);
202  c7->Divide(3,1);
203  int rebin = 1;
204  c7->cd(1);
205  gPad->SetLeftMargin(0.12);
206  gPad->SetRightMargin(0.01);
207  TH1D *hpy3 = new TH1D("hpy3","TPC inner clusters",2000,-0.10, 0.10);
208  delta_rphi->ProjectionY("hpy3",n_maps_layers+n_intt_layers+1,n_maps_layers+n_intt_layers+n_tpc_layers_inner);
209  hpy3->GetXaxis()->SetRangeUser(-0.04, 0.04);
210  hpy3->GetXaxis()->SetNdivisions(506);
211  hpy3->GetXaxis()->SetTitle("r#phi cluster error (cm)");
212  hpy3->GetXaxis()->SetTitleOffset(1.1);
213  hpy3->GetXaxis()->SetTitleSize(0.05);
214  hpy3->GetXaxis()->SetLabelSize(0.06);
215  hpy3->GetYaxis()->SetLabelSize(0.06);
216  hpy3->Rebin(rebin);
217  hpy3->Draw();
218  hpy3->Fit(fg,"R");
219  //double rms3 = 10000 * hpy3->GetRMS();
220  double rms3 = fg->GetParameter(2) * 10000.0;
221  double rms3_err = fg->GetParError(2) * 10000.0;
222  sprintf(label,"RMS %.1f #pm %.1f #mu m",rms3, rms3_err);
223  TLatex *l3 = new TLatex(0.45,0.92,label);
224  l3->SetNDC(1);
225  l3->Draw();
226 
227  c7->cd(2);
228  gPad->SetLeftMargin(0.12);
229  gPad->SetRightMargin(0.01);
230  TH1D *hpy4 = new TH1D("hpy4","TPC mid clusters",2000,-0.10, 0.10);
231  delta_rphi->ProjectionY("hpy4",n_maps_layers+n_intt_layers+n_tpc_layers_inner+1,n_maps_layers+n_intt_layers+n_tpc_layers_inner+n_tpc_layers_mid);
232  hpy4->GetXaxis()->SetRangeUser(-0.04, 04);
233  hpy4->GetXaxis()->SetNdivisions(506);
234  hpy4->GetXaxis()->SetTitle("r#phi cluster error (cm)");
235  hpy4->GetXaxis()->SetTitleOffset(1.1);
236  hpy4->GetXaxis()->SetTitleSize(0.05);
237  hpy4->GetXaxis()->SetLabelSize(0.06);
238  hpy4->GetYaxis()->SetLabelSize(0.06);
239  hpy4->Rebin(rebin);
240  hpy4->Draw();
241  hpy4->Fit(fg,"R");
242  //double rms4 = 10000 * hpy4->GetRMS();
243  double rms4 = fg->GetParameter(2) * 10000.0;
244  double rms4_err = fg->GetParError(2) * 10000.0;
245  sprintf(label,"RMS %.1f #pm %.1f #mu m",rms4, rms4_err);
246  TLatex *l4 = new TLatex(0.45,0.92,label);
247  l4->SetNDC(1);
248  l4->Draw();
249 
250  c7->cd(3);
251  gPad->SetLeftMargin(0.12);
252  gPad->SetRightMargin(0.01);
253  TH1D *hpy5 = new TH1D("hpy5","TPC outer clusters",2000,-0.10, 0.10);
254  delta_rphi->ProjectionY("hpy5",n_maps_layers+n_intt_layers+n_tpc_layers_inner+n_tpc_layers_mid+1,n_maps_layers+n_intt_layers+n_tpc_layers_inner+n_tpc_layers_mid+n_tpc_layers_outer);
255  //hpy5->GetXaxis()->SetRangeUser(-0.04, 04);
256  hpy5->GetXaxis()->SetNdivisions(506);
257  hpy5->GetXaxis()->SetTitle("r#phi cluster error (cm)");
258  hpy5->GetXaxis()->SetTitleOffset(1.1);
259  hpy5->GetXaxis()->SetTitleSize(0.05);
260  hpy5->GetXaxis()->SetLabelSize(0.06);
261  hpy5->GetYaxis()->SetLabelSize(0.06);
262  hpy5->Rebin(rebin);
263  hpy5->Draw();
264  hpy5->Fit(fg,"R");
265  //double rms5 = 10000 * hpy5->GetRMS();
266  double rms5 = fg->GetParameter(2) * 10000.0;
267  double rms5_err = fg->GetParError(2) * 10000.0;
268  sprintf(label,"RMS %.1f #pm %.1f #mu m",rms5, rms5_err);
269  TLatex *l5 = new TLatex(0.45,0.92,label);
270  l5->SetNDC(1);
271  l5->Draw();
272 
273 
274  TCanvas *c27 = new TCanvas("c27","c27",50,50,1200,800);
275  c27->Divide(3,2);
276 
277  c27->cd(1);
278  gPad->SetLeftMargin(0.12);
279  gPad->SetRightMargin(0.01);
280  TH1D *hpz1 = new TH1D("hpz1","MVTX clusters Z",500, -1.0, 1.0);
281  delta_z->ProjectionY("hpz1",1, n_maps_layers);
282  hpz1->GetXaxis()->SetTitle("Z cluster error (cm)");
283  hpz1->SetTitleOffset(0.1,"X");
284  hpz1->GetXaxis()->SetTitleSize(0.05);
285  hpz1->GetXaxis()->SetLabelSize(0.06);
286  hpz1->GetYaxis()->SetLabelSize(0.06);
287  hpz1->GetXaxis()->SetNdivisions(506);
288  hpz1->GetXaxis()->SetTitleOffset(1.1);
289  hpz1->GetXaxis()->SetRangeUser(-0.0016, 0.0016);
290  hpz1->Draw();
291  double zrms1 = 10000.0 * hpz1->GetRMS();
292  //char label[500];
293  sprintf(label,"RMS %.1f #mu m",zrms1);
294  TLatex *lz1 = new TLatex(0.55,0.92,label);
295  lz1->SetNDC(1);
296  lz1->Draw();
297 
298  c27->cd(2);
299  gPad->SetLeftMargin(0.12);
300  gPad->SetRightMargin(0.01);
301  //TH1D *hpz2 = new TH1D("hpz2","INTT clusters Z",500, -1.0, 1.0);
302  TH1D *hpz2 = new TH1D("hpz2","INTT clusters Z",500, -0.020, 0.020);
303  delta_z->ProjectionY("hpz2",n_maps_layers+1,n_maps_layers+1);
304  hpz2->GetXaxis()->SetRangeUser(-0.02, 0.02);
305  hpz2->GetXaxis()->SetTitle("Z cluster error (cm)");
306  hpz2->GetXaxis()->SetTitleOffset(0.6);
307  hpz2->GetXaxis()->SetTitleSize(0.05);
308  hpz2->GetXaxis()->SetLabelSize(0.06);
309  hpz2->GetYaxis()->SetLabelSize(0.06);
310  hpz2->GetXaxis()->SetNdivisions(506);
311  hpz2->GetXaxis()->SetTitleOffset(1.1);
312  hpz2->Draw();
313  double zrms2 = 10000 * hpz2->GetRMS();
314  sprintf(label,"RMS %.1f #mu m",zrms2);
315  TLatex *lz2 = new TLatex(0.55,0.92,label);
316  lz2->SetNDC(1);
317  lz2->Draw();
318 
319  c27->cd(3);
320  gPad->SetLeftMargin(0.12);
321  gPad->SetRightMargin(0.01);
322  TH1D *hpz3 = new TH1D("hpz3","TPC inner clusters Z",500,-0.2, 0.2);
323  delta_z->ProjectionY("hpz3",n_maps_layers+n_intt_layers+1, n_maps_layers+n_intt_layers+n_tpc_layers_inner);
324  hpz3->GetXaxis()->SetNdivisions(506);
325  hpz3->GetXaxis()->SetTitle("Z cluster error (cm)");
326  hpz3->GetXaxis()->SetTitleOffset(1.1);
327  hpz3->GetXaxis()->SetTitleSize(0.05);
328  hpz3->GetXaxis()->SetLabelSize(0.06);
329  hpz3->GetYaxis()->SetLabelSize(0.06);
330  hpz3->GetXaxis()->SetRangeUser(-0.2, 0.2);
331  hpz3->Draw();
332  double zrms3 = 10000 * hpz3->GetRMS();
333  sprintf(label,"RMS %.1f #mu m",zrms3);
334  TLatex *lz3 = new TLatex(0.55,0.92,label);
335  lz3->SetNDC(1);
336  lz3->Draw();
337 
338  c27->cd(4);
339  gPad->SetLeftMargin(0.12);
340  gPad->SetRightMargin(0.01);
341  TH1D *hpz4 = new TH1D("hpz4","TPC mid clusters Z",500,-0.2, 0.2);
342  delta_z->ProjectionY("hpz4",n_maps_layers+n_intt_layers+n_tpc_layers_inner+1, n_maps_layers+n_intt_layers+n_tpc_layers_inner+n_tpc_layers_mid);
343  hpz4->GetXaxis()->SetNdivisions(506);
344  hpz4->GetXaxis()->SetTitle("Z cluster error (cm)");
345  hpz4->GetXaxis()->SetTitleOffset(1.1);
346  hpz4->GetXaxis()->SetTitleSize(0.05);
347  hpz4->GetXaxis()->SetLabelSize(0.06);
348  hpz4->GetYaxis()->SetLabelSize(0.06);
349  hpz4->GetXaxis()->SetRangeUser(-0.2, 0.2);
350  hpz4->Draw();
351  double zrms4 = 10000 * hpz4->GetRMS();
352  sprintf(label,"RMS %.1f #mu m",zrms4);
353  TLatex *lz4 = new TLatex(0.55,0.92,label);
354  lz4->SetNDC(1);
355  lz4->Draw();
356 
357  c27->cd(5);
358  gPad->SetLeftMargin(0.12);
359  gPad->SetRightMargin(0.01);
360  TH1D *hpz5 = new TH1D("hpz5","TPC outer clusters Z",500,-0.2, 0.2);
361  delta_z->ProjectionY("hpz5",n_maps_layers+n_intt_layers+n_tpc_layers_inner+n_tpc_layers_mid+1, n_maps_layers+n_intt_layers+n_tpc_layers_inner+n_tpc_layers_mid+n_tpc_layers_outer);
362  hpz5->GetXaxis()->SetNdivisions(506);
363  hpz5->GetXaxis()->SetTitle("Z cluster error (cm)");
364  hpz5->GetXaxis()->SetTitleOffset(1.1);
365  hpz5->GetXaxis()->SetTitleSize(0.05);
366  hpz5->GetXaxis()->SetLabelSize(0.06);
367  hpz5->GetYaxis()->SetLabelSize(0.06);
368  hpz5->GetXaxis()->SetRangeUser(-0.2, 0.2);
369  hpz5->Draw();
370  double zrms5 = 10000 * hpz5->GetRMS();
371  sprintf(label,"RMS %.1f #mu m",zrms5);
372  TLatex *lz5 = new TLatex(0.55,0.92,label);
373  lz5->SetNDC(1);
374  lz5->Draw();
375 
376 
377 
378  TCanvas *c6 = new TCanvas("c6","c6",50,50,1200,800);
379  c6->Divide(2,1);
380  c6->cd(1);
381  gPad->SetLeftMargin(0.12);
382  gPad->SetRightMargin(0.01);
383  TH1D *hpy1 = new TH1D("hpy1","MVTX clusters",2000, -0.05, 0.05);
384  delta_rphi->ProjectionY("hpy1",1,n_maps_layers);
385  //hpy1->GetXaxis()->SetRangeUser(-0.0016, 0.0016);
386  //hpy1->GetXaxis()->SetRangeUser(-0.011, 0.011);
387  hpy1->GetXaxis()->SetRangeUser(-0.002, 0.002);
388  hpy1->GetXaxis()->SetTitle("r#phi cluster error (cm)");
389  hpy1->SetTitleOffset(0.1,"X");
390  hpy1->GetXaxis()->SetTitleSize(0.05);
391  hpy1->GetXaxis()->SetLabelSize(0.06);
392  hpy1->GetYaxis()->SetLabelSize(0.06);
393  hpy1->GetXaxis()->SetNdivisions(506);
394  hpy1->GetXaxis()->SetTitleOffset(1.1);
395  hpy1->Draw();
396  double rms1 = 10000.0 * hpy1->GetRMS();
397  sprintf(label,"RMS %.1f #mu m",rms1);
398  TLatex *l1 = new TLatex(0.55,0.92,label);
399  l1->SetNDC(1);
400  l1->Draw();
401 
402  //delta_rphi->Draw();
403 
404  c6->cd(2);
405  gPad->SetLeftMargin(0.12);
406  gPad->SetRightMargin(0.01);
407  TH1D *hpy2 = new TH1D("hpy2","INTT clusters (type 1)",2000, -0.05, 0.05);
408  delta_rphi->ProjectionY("hpy2",n_maps_layers+2,n_maps_layers + n_intt_layers);
409  hpy2->GetXaxis()->SetRangeUser(-0.011, 0.011);
410  hpy2->GetXaxis()->SetTitle("r#phi cluster error (cm)");
411  hpy2->GetXaxis()->SetTitleOffset(0.6);
412  hpy2->GetXaxis()->SetTitleSize(0.05);
413  hpy2->GetXaxis()->SetLabelSize(0.06);
414  hpy2->GetYaxis()->SetLabelSize(0.06);
415  hpy2->GetXaxis()->SetNdivisions(506);
416  hpy2->GetXaxis()->SetTitleOffset(1.1);
417  hpy2->Draw();
418  double rms2 = 10000 * hpy2->GetRMS();
419  sprintf(label,"RMS %.1f #mu m",rms2);
420  TLatex *l2 = new TLatex(0.55,0.92,label);
421  l2->SetNDC(1);
422  l2->Draw();
423 
424  /*
425  c6->cd(3);
426  gPad->SetLeftMargin(0.12);
427  gPad->SetRightMargin(0.01);
428  TH1D *hpy9 = new TH1D("hpy9","INTT clusters",200, -1.0, 1.0);
429  delta_rphi->ProjectionY("hpy9",n_maps_layers+1,n_maps_layers + n_intt_layers-3);
430  hpy9->GetXaxis()->SetRangeUser(-0.011, 0.011);
431  hpy9->GetXaxis()->SetTitle("r#phi cluster error (cm)");
432  hpy9->GetXaxis()->SetTitleOffset(0.6);
433  hpy9->GetXaxis()->SetTitleSize(0.05);
434  hpy9->GetXaxis()->SetLabelSize(0.06);
435  hpy9->GetYaxis()->SetLabelSize(0.06);
436  hpy9->GetXaxis()->SetNdivisions(506);
437  hpy9->GetXaxis()->SetTitleOffset(1.1);
438  hpy9->Draw();
439  double rms9 = 10000 * hpy9->GetRMS();
440  sprintf(label,"RMS %.1f #mu m",rms9);
441  TLatex *l9 = new TLatex(0.55,0.92,label);
442  l9->SetNDC(1);
443  l9->Draw();
444  */
445  /*
446  TCanvas *c8 = new TCanvas("C8","C8",50,50,1200,700);
447  c8->Divide(2,1);
448  c8->cd(1);
449  gPad->SetLeftMargin(0.15);
450  gPad->SetGrid();
451  delta_rphi->GetYaxis()->SetTitleOffset(1.7);
452  delta_rphi->Draw("p");
453  c8->cd(2);
454  gPad->SetLeftMargin(0.18);
455  gPad->SetGrid();
456  delta_phi->GetYaxis()->SetTitleOffset(2.3);
457  delta_phi->Draw("p");
458 
459 
460  TCanvas *c8z = new TCanvas("C8Z","C8Z",50,50,800,800);
461  c8z->SetLeftMargin(0.15);
462  delta_z->GetYaxis()->SetTitleOffset(1.7);
463  delta_z->Draw("p");
464 
465  TCanvas *c9 = new TCanvas("c9","c9",5,5,800,600);
466  delta_phi_gphi->Draw();
467  */
468 
469 }