CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
bin_deltaZ.cxx
Go to the documentation of this file.
1int bin_deltaZ()
2{
3 double M_PI = 3.1415926;
4///*
5//////////////////////////////////////////////////////////
6// This file has been automatically generated
7// (Thu Feb 1 20:29:58 2018 by ROOT version5.34/09)
8// from TTree hit/hit
9// found on file: hough_hist_test.root
10//////////////////////////////////////////////////////////
11
12
13//Reset ROOT and connect tree file
14 //gROOT->Reset();
15 //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
16 //if (!f) {
17 //f = new TFile("hough_hist_test.root");
18 //}
19 //f->GetObject("hit",tree);
20 TChain * hit = new TChain("hit");
21 hit->Add("hough_hist_test.root");
22
23//Declaration of leaves types
24 Int_t hit_run;
25 Int_t hit_evt;
26 Int_t hit_nhit;
27 Int_t hit_hitid[10000];
28 Int_t hit_layer[10000];
29 Int_t hit_wire[10000];
30 Double_t hit_x[10000];
31 Double_t hit_y[10000];
32 Double_t hit_z[10000];
33 Double_t hit_driftdist[10000];
34 Double_t hit_drifttime[10000];
35 Int_t hit_flag[10000];
36 Double_t hit_truth_x[10000];
37 Double_t hit_truth_y[10000];
38 Double_t hit_truth_z[10000];
39 Double_t hit_truth_drift[10000];
40 Int_t hit_truth_ambig[10000];
41
42 // Set branch addresses.
43 hit->SetBranchAddress("hit_run",&hit_run);
44 hit->SetBranchAddress("hit_evt",&hit_evt);
45 hit->SetBranchAddress("hit_nhit",&hit_nhit);
46 hit->SetBranchAddress("hit_hitid",hit_hitid);
47 hit->SetBranchAddress("hit_layer",hit_layer);
48 hit->SetBranchAddress("hit_wire",hit_wire);
49 hit->SetBranchAddress("hit_x",hit_x);
50 hit->SetBranchAddress("hit_y",hit_y);
51 hit->SetBranchAddress("hit_z",hit_z);
52 hit->SetBranchAddress("hit_driftdist",hit_driftdist);
53 hit->SetBranchAddress("hit_drifttime",hit_drifttime);
54 hit->SetBranchAddress("hit_flag",hit_flag);
55 hit->SetBranchAddress("hit_truth_x",hit_truth_x);
56 hit->SetBranchAddress("hit_truth_y",hit_truth_y);
57 hit->SetBranchAddress("hit_truth_z",hit_truth_z);
58 hit->SetBranchAddress("hit_truth_drift",hit_truth_drift);
59 hit->SetBranchAddress("hit_truth_ambig",hit_truth_ambig);
60
61// This is the loop skeleton
62// To read only selected branches, Insert statements like:
63// hit->SetBranchStatus("*",0); // disable all branches
64// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
65//*/
66
67//////////////////////////////////////////////////////////
68// This file has been automatically generated
69// (Thu Feb 1 20:29:46 2018 by ROOT version5.34/09)
70// from TTree hot/hot
71// found on file: hough_hist_test.root
72//////////////////////////////////////////////////////////
73
74
75//Reset ROOT and connect tree file
76 //gROOT->Reset();
77 //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
78 //if (!f) {
79 //f = new TFile("hough_hist_test.root");
80 //}
81 //f->GetObject("hot",tree);
82 TChain * hot = new TChain("hot");
83 hot->Add("hough_hist_test.root");
84
85//Declaration of leaves types
86 Int_t hot_run;
87 Int_t hot_evt;
88 Int_t hot_trk;
89 Int_t hot_nhot;
90 Int_t hot_hitid[1000];
91 Int_t hot_layer[1000];
92 Int_t hot_wire[1000];
93 Double_t hot_x[1000];
94 Double_t hot_y[1000];
95 Double_t hot_z[1000];
96 Double_t hot_x0[1000];
97 Double_t hot_y0[1000];
98 Double_t hot_z0[1000];
99 Double_t hot_s0[1000];
100 Double_t hot_x1[1000];
101 Double_t hot_y1[1000];
102 Double_t hot_z1[1000];
103 Double_t hot_s1[1000];
104 Double_t hot_drift[1000];
105 Int_t hot_flag[1000];
106 Double_t hot_deltaD[1000];
107 Double_t hot_truth_x[1000];
108 Double_t hot_truth_y[1000];
109 Double_t hot_truth_z[1000];
110 Double_t hot_truth_drift[1000];
111 Int_t hot_truth_ambig[1000];
112
113 // Set branch addresses.
114 hot->SetBranchAddress("hot_run",&hot_run);
115 hot->SetBranchAddress("hot_evt",&hot_evt);
116 hot->SetBranchAddress("hot_trk",&hot_trk);
117 hot->SetBranchAddress("hot_nhot",&hot_nhot);
118 hot->SetBranchAddress("hot_hitid",hot_hitid);
119 hot->SetBranchAddress("hot_layer",hot_layer);
120 hot->SetBranchAddress("hot_wire",hot_wire);
121 hot->SetBranchAddress("hot_x",hot_x);
122 hot->SetBranchAddress("hot_y",hot_y);
123 hot->SetBranchAddress("hot_z",hot_z);
124 hot->SetBranchAddress("hot_x0",hot_x0);
125 hot->SetBranchAddress("hot_y0",hot_y0);
126 hot->SetBranchAddress("hot_z0",hot_z0);
127 hot->SetBranchAddress("hot_s0",hot_s0);
128 hot->SetBranchAddress("hot_x1",hot_x1);
129 hot->SetBranchAddress("hot_y1",hot_y1);
130 hot->SetBranchAddress("hot_z1",hot_z1);
131 hot->SetBranchAddress("hot_s1",hot_s1);
132 hot->SetBranchAddress("hot_drift",hot_drift);
133 hot->SetBranchAddress("hot_flag",hot_flag);
134 hot->SetBranchAddress("hot_deltaD",hot_deltaD);
135 hot->SetBranchAddress("hot_truth_x",hot_truth_x);
136 hot->SetBranchAddress("hot_truth_y",hot_truth_y);
137 hot->SetBranchAddress("hot_truth_z",hot_truth_z);
138 hot->SetBranchAddress("hot_truth_drift",hot_truth_drift);
139 hot->SetBranchAddress("hot_truth_ambig",hot_truth_ambig);
140
141// This is the loop skeleton
142// To read only selected branches, Insert statements like:
143// hot->SetBranchStatus("*",0); // disable all branches
144// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
145
146
147//////////////////////////////////////////////////////////
148// This file has been automatically generated
149// (Thu Feb 1 19:24:40 2018 by ROOT version5.34/09)
150// from TTree trk/trk
151// found on file: hough_hist_test.root
152//////////////////////////////////////////////////////////
153
154///*
155//Reset ROOT and connect tree file
156 //gROOT->Reset();
157 //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
158 //if (!f) {
159 //f = new TFile("hough_hist_test.root");
160 //}
161 //TTree* tree;
162 //f->GetObject("trk",tree);
163 TChain * trk = new TChain("trk");
164 trk->Add("hough_hist_test.root");
165
166//Declaration of leaves types
167 Int_t trk_run;
168 Int_t trk_evt;
169 Int_t trk_ntrk;
170 Int_t trk_trackId[100];
171 Int_t trk_charge[100];
172 Double_t trk_dr[100];
173 Double_t trk_phi0[100];
174 Double_t trk_kappa[100];
175 Double_t trk_dz[100];
176 Double_t trk_tanl[100];
177 Double_t trk_pxy[100];
178 Double_t trk_px[100];
179 Double_t trk_py[100];
180 Double_t trk_pz[100];
181 Double_t trk_p[100];
182 Double_t trk_theta[100];
183 Double_t trk_phi[100];
184 Double_t trk_x[100];
185 Double_t trk_y[100];
186 Double_t trk_z[100];
187 Double_t trk_r[100];
188 Double_t trk_chi2[100];
189 Double_t trk_fiTerm[100];
190 Double_t trk_matchChi2[100];
191 Int_t trk_nhit[100];
192 Int_t trk_ncluster[100];
193 Int_t trk_stat[100];
194 Int_t trk_ndof[100];
195 Int_t trk_nster[100];
196 Int_t trk_nlayer[100];
197 Int_t trk_firstLayer[100];
198 Int_t trk_lastLayer[100];
199 Int_t trk_nCgemXClusters[100];
200 Int_t trk_nCgemVClusters[100];
201 Int_t trk_nhop[100];
202 Int_t trk_nhot[100];
203 Double_t trk_Xc[100];
204 Double_t trk_Yc[100];
205 Double_t trk_R[100];
206 Int_t trk_truth_charge[100];
207 Double_t trk_truth_dr[100];
208 Double_t trk_truth_phi0[100];
209 Double_t trk_truth_kappa[100];
210 Double_t trk_truth_dz[100];
211 Double_t trk_truth_tanl[100];
212 Double_t trk_truth_pxy[100];
213 Double_t trk_truth_px[100];
214 Double_t trk_truth_py[100];
215 Double_t trk_truth_pz[100];
216 Double_t trk_truth_p[100];
217 Double_t trk_truth_theta[100];
218 Double_t trk_truth_phi[100];
219 Double_t trk_truth_x[100];
220 Double_t trk_truth_y[100];
221 Double_t trk_truth_z[100];
222 Double_t trk_truth_r[100];
223 Double_t trk_truth_cosTheta[100];
224 Double_t trk_truth_Xc[100];
225 Double_t trk_truth_Yc[100];
226 Double_t trk_truth_R[100];
227
228 // Set branch addresses.
229 trk->SetBranchAddress("trk_run",&trk_run);
230 trk->SetBranchAddress("trk_evt",&trk_evt);
231 trk->SetBranchAddress("trk_ntrk",&trk_ntrk);
232 trk->SetBranchAddress("trk_trackId",trk_trackId);
233 trk->SetBranchAddress("trk_charge",trk_charge);
234 trk->SetBranchAddress("trk_dr",trk_dr);
235 trk->SetBranchAddress("trk_phi0",trk_phi0);
236 trk->SetBranchAddress("trk_kappa",trk_kappa);
237 trk->SetBranchAddress("trk_dz",trk_dz);
238 trk->SetBranchAddress("trk_tanl",trk_tanl);
239 trk->SetBranchAddress("trk_pxy",trk_pxy);
240 trk->SetBranchAddress("trk_px",trk_px);
241 trk->SetBranchAddress("trk_py",trk_py);
242 trk->SetBranchAddress("trk_pz",trk_pz);
243 trk->SetBranchAddress("trk_p",trk_p);
244 trk->SetBranchAddress("trk_theta",trk_theta);
245 trk->SetBranchAddress("trk_phi",trk_phi);
246 trk->SetBranchAddress("trk_x",trk_x);
247 trk->SetBranchAddress("trk_y",trk_y);
248 trk->SetBranchAddress("trk_z",trk_z);
249 trk->SetBranchAddress("trk_r",trk_r);
250 trk->SetBranchAddress("trk_chi2",trk_chi2);
251 trk->SetBranchAddress("trk_fiTerm",trk_fiTerm);
252 trk->SetBranchAddress("trk_matchChi2",trk_matchChi2);
253 trk->SetBranchAddress("trk_nhit",trk_nhit);
254 trk->SetBranchAddress("trk_ncluster",trk_ncluster);
255 trk->SetBranchAddress("trk_stat",trk_stat);
256 trk->SetBranchAddress("trk_ndof",trk_ndof);
257 trk->SetBranchAddress("trk_nster",trk_nster);
258 trk->SetBranchAddress("trk_nlayer",trk_nlayer);
259 trk->SetBranchAddress("trk_firstLayer",trk_firstLayer);
260 trk->SetBranchAddress("trk_lastLayer",trk_lastLayer);
261 trk->SetBranchAddress("trk_nCgemXClusters",trk_nCgemXClusters);
262 trk->SetBranchAddress("trk_nCgemVClusters",trk_nCgemVClusters);
263 trk->SetBranchAddress("trk_nhop",trk_nhop);
264 trk->SetBranchAddress("trk_nhot",trk_nhot);
265 trk->SetBranchAddress("trk_Xc",trk_Xc);
266 trk->SetBranchAddress("trk_Yc",trk_Yc);
267 trk->SetBranchAddress("trk_R",trk_R);
268 trk->SetBranchAddress("trk_truth_charge",trk_truth_charge);
269 trk->SetBranchAddress("trk_truth_dr",trk_truth_dr);
270 trk->SetBranchAddress("trk_truth_phi0",trk_truth_phi0);
271 trk->SetBranchAddress("trk_truth_kappa",trk_truth_kappa);
272 trk->SetBranchAddress("trk_truth_dz",trk_truth_dz);
273 trk->SetBranchAddress("trk_truth_tanl",trk_truth_tanl);
274 trk->SetBranchAddress("trk_truth_pxy",trk_truth_pxy);
275 trk->SetBranchAddress("trk_truth_px",trk_truth_px);
276 trk->SetBranchAddress("trk_truth_py",trk_truth_py);
277 trk->SetBranchAddress("trk_truth_pz",trk_truth_pz);
278 trk->SetBranchAddress("trk_truth_p",trk_truth_p);
279 trk->SetBranchAddress("trk_truth_theta",trk_truth_theta);
280 trk->SetBranchAddress("trk_truth_phi",trk_truth_phi);
281 trk->SetBranchAddress("trk_truth_x",trk_truth_x);
282 trk->SetBranchAddress("trk_truth_y",trk_truth_y);
283 trk->SetBranchAddress("trk_truth_z",trk_truth_z);
284 trk->SetBranchAddress("trk_truth_r",trk_truth_r);
285 trk->SetBranchAddress("trk_truth_cosTheta",trk_truth_cosTheta);
286 trk->SetBranchAddress("trk_truth_Xc",trk_truth_Xc);
287 trk->SetBranchAddress("trk_truth_Yc",trk_truth_Yc);
288 trk->SetBranchAddress("trk_truth_R",trk_truth_R);
289
290// This is the loop skeleton
291// To read only selected branches, Insert statements like:
292// trk->SetBranchStatus("*",0); // disable all branches
293// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
294
295 TGraphErrors *gr_sigma_1_3_kb = new TGraphErrors();
296 TGraphErrors *gr_sigma_21_36_kb = new TGraphErrors();
297 TGraphErrors *gr_RMS_1_3_kb = new TGraphErrors();
298 TGraphErrors *gr_RMS_21_36_kb = new TGraphErrors();
299
300 TGraphErrors *gr_sigma_1_3_tr = new TGraphErrors();
301 TGraphErrors *gr_sigma_21_36_tr = new TGraphErrors();
302 TGraphErrors *gr_RMS_1_3_tr = new TGraphErrors();
303 TGraphErrors *gr_RMS_21_36_tr = new TGraphErrors();
304
305 TGraphErrors *gr_tanl_kb = new TGraphErrors();
306 TGraphErrors *gr_tanl_tr = new TGraphErrors();
307 TGraphErrors *gr_dz_kb = new TGraphErrors();
308 TGraphErrors *gr_dz_tr = new TGraphErrors();
309
310 TGraph* gr_sigma_kb = new TGraph();
311 TGraph* gr_sigma1_kb= new TGraph();
312 TGraph* gr_sigma2_kb= new TGraph();
313 TGraph* gr_sigma_tr = new TGraph();
314 TGraph* gr_sigma1_tr= new TGraph();
315 TGraph* gr_sigma2_tr= new TGraph();
316
317 TAxis *x_axis,*y_axis;
318 stringstream ssname;
319 string sname;
320 const char * name;
321 int point = 10;
322 int yMax;
323 int binx,biny,binz;
324 int layer[20] = {1,2,3,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,0};
325 int ibin = 0;
326 int nbin = 100;
327 while(nbin <= 2000){
328 int x_bin = nbin;
329 int y_bin = nbin;
330 double x_max = 0.93/sqrt(1-0.93*0.93);
331 double y_max = 50;
332
333 ssname.str("");
334 ssname <<nbin<<"bin,tanl";
335 sname = ssname.str();
336 name = sname.c_str();
337 TH1D *h_tanl_kb = new TH1D(name,name,100,-0.2,0.2);
338 ssname.str("");
339 ssname <<nbin<<"bin,dz";
340 sname = ssname.str();
341 name = sname.c_str();
342 TH1D *h_dz_kb = new TH1D(name,name,100,-2,2);
343
344 ssname.str("");
345 ssname <<nbin<<"bin,tanl,tr";
346 sname = ssname.str();
347 name = sname.c_str();
348 TH1D *h_tanl_tr = new TH1D(name,name,100,-0.2,0.2);
349 ssname.str("");
350 ssname <<nbin<<"bin,dz,tr";
351 sname = ssname.str();
352 name = sname.c_str();
353 TH1D *h_dz_tr = new TH1D(name,name,100,-2,2);
354
355 ssname.str("");
356 ssname <<nbin<<"bin,1~3 layer";
357 sname = ssname.str();
358 name = sname.c_str();
359 TH1D *h_1_3_kb = new TH1D(name,name,100,-1,1);
360 ssname.str("");
361 ssname <<nbin<<"bin,21~36 layer";
362 sname = ssname.str();
363 name = sname.c_str();
364 TH1D *h_21_36_kb = new TH1D(name,name,100,-10,10);
365
366 ssname.str("");
367 ssname <<nbin<<"bin,1~3 layer,tr";
368 sname = ssname.str();
369 name = sname.c_str();
370 TH1D *h_1_3_tr = new TH1D(name,name,100,-1,1);
371 ssname.str("");
372 ssname <<nbin<<"bin,21~36 layer,tr";
373 sname = ssname.str();
374 name = sname.c_str();
375 TH1D *h_21_36_tr = new TH1D(name,name,100,-10,10);
376
377 TH1D *h_deltaZ_kb[20];
378 TH1D *h_deltaZ_tr[20];
379 for(int l=1;l<=20;l++){
380 double cut;
381 ssname.str("");
382 ssname <<nbin<<"bin, ";
383 if(l<=3){ssname <<"layer "<<l;cut=1;}
384 else if(l<=19){ssname <<"layer "<<l+17;cut=10;}
385 else if(l==20){ssname <<"all layer";cut=10;}
386 sname = ssname.str();
387 name = sname.c_str();
388 h_deltaZ_kb[l-1]= new TH1D(name,name,100,-cut,cut);
389 //h_deltaZ_kb[l-1]= new TH1D(name,name,100,-10,10);
390
391 ssname.str("");
392 ssname <<nbin<<"bin, ";
393 if(l<=3){ssname <<"tr,layer "<<l;cut=1;}
394 else if(l<=19){ssname <<"tr,layer "<<l+17;cut=10;}
395 else if(l==20){ssname <<"tr,all layer ";cut=10;}
396 sname = ssname.str();
397 name = sname.c_str();
398 h_deltaZ_tr[l-1]= new TH1D(name,name,100,-cut,cut);
399 //h_deltaZ_tr[l-1]= new TH1D(name,name,100,-10,10);
400 }
401
402 int event =11 ;
403 Long64_t nentries = hot->GetEntries();
404 Long64_t nbytes = 0;
405 for (Long64_t i=0; i<nentries;i++)
406 //for (Long64_t i=0; i<1000;i++)
407 //for (Long64_t i=event; i<event+1;i++)
408 {
409 nbytes += hot->GetEntry(i);
410 TH2D* hough_kb = new TH2D("hough_kb","hough_kb",x_bin,-x_max,x_max,y_bin,-y_max,y_max);
411 TH2D* hough_tr = new TH2D("hough_tr","hough_tr",x_bin,0,M_PI,y_bin,-y_max,y_max);
412 TH2D *houghhit_kb[1000];
413 TH2D *houghhit_tr[1000];
414 int nhit(0),ihit(0);
415 for(int j=0;j<hot_nhot;j++){
416 if(hot_flag[j] ==0)continue;
417 if(hot_layer[j] <3 && hot_flag[j] ==1)continue;
418 nhit++;
419 //cout<<hot_layer[j]<<endl;
420 ssname.str("");
421 ssname <<"layer:"<<hot_layer[j]<<" ,wire,kb"<<hot_wire[j];
422 sname = ssname.str();
423 name = sname.c_str();
424 houghhit_kb[ihit] = new TH2D(name,name,x_bin,-x_max,x_max,y_bin,-y_max,y_max);
425 fill_histogram(hot_s0[j],hot_z0[j],x_max,nbin*point,houghhit_kb[ihit]);
426 if(hot_layer[j]>3)fill_histogram(hot_s1[j],hot_z1[j],x_max,nbin*point,houghhit_kb[ihit]);
427 hough_kb->Add(houghhit_kb[ihit]);
428
429 ssname.str("");
430 ssname <<"layer:"<<hot_layer[j]<<" ,wire,tr"<<hot_wire[j];
431 sname = ssname.str();
432 name = sname.c_str();
433 houghhit_tr[ihit] = new TH2D(name,name,x_bin,0,M_PI,y_bin,-y_max,y_max);
434 fill_histogram2(hot_s0[j],hot_z0[j],M_PI,nbin*point,houghhit_tr[ihit]);
435 if(hot_layer[j] >3)fill_histogram2(hot_s1[j],hot_z1[j],M_PI,nbin*point,houghhit_tr[ihit]);
436 hough_tr->Add(houghhit_tr[ihit]);
437 ihit++;
438 }
439
440 //TCanvas *c1 = new TCanvas("houghspace","houghspace",1000,1000 );
441 //c1->Divide(2,2);
442 //c1->cd(1);
443 //hough_kb->Draw();
444 //c1->cd(2);
445 //hough_kb->Draw("LEGO2Z");
446 //c1->cd(3);
447 //hough_tr->Draw();
448 //c1->cd(4);
449 //hough_tr->Draw("LEGO2Z");
450 //c1->SaveAs("houghspace.pdf");
451
452///*
453 hough_kb->GetMaximumBin(binx,biny,binz);
454 //int height = hough_kb->GetBinContent(binx,biny);
455 double tanl_kb = hough_kb->GetXaxis()->GetBinCenter(binx);
456 double dz_kb = hough_kb->GetYaxis()->GetBinCenter(biny);
457
458 hough_tr->GetMaximumBin(binx,biny,binz);
459 //int height = hough_tr->GetBinContent(binx,biny);
460 double tanl_tr = -1./tan(hough_tr->GetXaxis()->GetBinCenter(binx));
461 double dz_tr = hough_tr->GetYaxis()->GetBinCenter(biny)/sin(hough_tr->GetXaxis()->GetBinCenter(binx));
462//*/
463
464///*
465 //TCanvas *c;
466 int nhot(0);
467 for(int ihit=0;ihit<nhit;ihit++){
468 //ssname.str("");
469 //ssname <<ihit<<".pdf";
470 //sname = ssname.str();
471 //name = sname.c_str();
472 //c = new TCanvas(name,name,700,500);
473 //houghhit_kb[ihit]->Draw();
474 //cout<<houghhit_kb[ihit]->GetBinContent(binx,biny)<<endl;
475 //cout<<layersOnPeak[ihit]<<endl;
476 //if(houghhit_kb[ihit]->GetBinContent(binx,biny)>0){
477 //nhot++;
478 //}
479 //cout<<endl;
480 //cout<<layersOnPeak[ihit]<<" "<<houghhit_kb[ihit]->GetBinContent(binx,biny)<<endl;
481 //c->SaveAs(name);
482 delete houghhit_kb[ihit];
483 delete houghhit_tr[ihit];
484 }
485
486//*/
487
488///*
489 double deltaZ(0);
490 for(int j=0;j<hot_nhot;j++){
491 if(hot_flag[j] ==0)continue;
492 if(hot_layer[j] <3 && hot_flag[j] ==1)continue;
493
494 double dz0(0),dz1(0);
495 dz0 = tanl_kb*hot_s0[j] + dz_kb - hot_z0[j];
496 dz1 = tanl_kb*hot_s1[j] + dz_kb - hot_z1[j];
497 deltaZ = fabs(dz0) < fabs(dz1) ? dz0 : dz1;
498 if(hot_layer[j] <3) {h_deltaZ_kb[hot_layer[j]]->Fill(deltaZ);h_1_3_kb->Fill(deltaZ);}
499 else {h_deltaZ_kb[hot_layer[j]-17]->Fill(deltaZ);h_21_36_kb->Fill(deltaZ);}
500 h_deltaZ_kb[19]->Fill(deltaZ);
501
502 dz0 = tanl_tr*hot_s0[j] + dz_tr - hot_z0[j];
503 dz1 = tanl_tr*hot_s1[j] + dz_tr - hot_z1[j];
504 deltaZ = fabs(dz0) < fabs(dz1) ? dz0 : dz1;
505 if(hot_layer[j] <3) {h_deltaZ_tr[hot_layer[j]]->Fill(deltaZ);h_1_3_tr->Fill(deltaZ);}
506 else {h_deltaZ_tr[hot_layer[j]-17]->Fill(deltaZ);h_21_36_tr->Fill(deltaZ);}
507 h_deltaZ_tr[19]->Fill(deltaZ);
508
509 }
510
511///*
512 for (Long64_t itrk=0; itrk<trk->GetEntries();itrk++)
513 {
514 trk->GetEntry(itrk);
515 if(trk_evt != hot_evt)continue;
516 for(int jtrk=0;jtrk<trk_ntrk;jtrk++){
517 if(jtrk != hot_trk)continue;
518 h_tanl_kb->Fill(trk_truth_tanl[jtrk] - tanl_kb);
519 h_dz_kb->Fill(trk_truth_dz[jtrk] - dz_kb);
520 h_tanl_tr->Fill(trk_truth_tanl[jtrk] - tanl_tr);
521 h_dz_tr->Fill(trk_truth_dz[jtrk] - dz_tr);
522 }
523 break;
524 }
525 delete hough_kb;
526 delete hough_tr;
527//*/
528 }
529
530///*
531 //TCanvas *c2 = new TCanvas("deltaZ","deltaZ",700,500 );
532 //h_deltaZ_kb[0]->Draw();
533 //h_deltaZ_kb[0]->SetLineColor(2);
534 //h_deltaZ_tr[0]->Draw("same");
535 //h_deltaZ_tr[0]->SetLineColor(4);
536 ssname.str("");
537 ssname <<nbin<<"bin each layer deltaZ";
538 sname = ssname.str();
539 name = sname.c_str();
540 TCanvas *c2 = new TCanvas(name,name,2500,2000 );
541 c2->Divide(5,4);
542 double sigma1_kb[20],sigma2_kb[20];
543 double sigma1_tr[20],sigma2_tr[20];
544 for(int l=0;l<20;l++){
545 c2->cd(l+1);
546 if(l<3){
547 h_deltaZ_kb[l]->Fit("gaus","Q");
548 sigma1_kb[l] = h_deltaZ_kb[l]->GetFunction("gaus")->GetParameter(2);
549 sigma2_kb[l] = h_deltaZ_kb[l]->GetFunction("gaus")->GetParameter(2);
550 h_deltaZ_tr[l]->Fit("gaus","Q");
551 sigma1_tr[l] = h_deltaZ_tr[l]->GetFunction("gaus")->GetParameter(2);
552 sigma2_tr[l] = h_deltaZ_tr[l]->GetFunction("gaus")->GetParameter(2);
553 }else{
554 TF1 *f_kb = new TF1("f_kb", fitterFunction2, -10.0, 10.0, 5);
555 TF1 *f_tb = new TF1("f_tr", fitterFunction2, -10.0, 10.0, 5);
556 int xbin,ybin,zbin;
557 h_deltaZ_kb[l]->GetMaximumBin(xbin,ybin,zbin);
558 double skew = h_deltaZ_kb[l]->GetSkewness();
559 double p0 = h_deltaZ_kb[l]->GetMaximum();
560 double p1 = h_deltaZ_kb[l]->GetMean();
561 double p2 = skew > 0 ? 2.0*h_deltaZ_kb[l]->GetRMS() : 0.5*h_deltaZ_kb[l]->GetRMS();
562 double p3 = skew < 0 ? 2.0*h_deltaZ_kb[l]->GetRMS() : 0.5*h_deltaZ_kb[l]->GetRMS();
563 double p4 = h_deltaZ_kb[l]->GetBinCenter(xbin);
564 f_kb->SetParameters(p0,p1,p2,p3,p4);
565 f_kb->SetLineColor(4);
566 f_kb->FixParameter(4, p4);
567 h_deltaZ_kb[l]->Fit("f1");
568 sigma1_kb[l] = f_kb->GetParameter(2);
569 sigma2_kb[l] = f_kb->GetParameter(3);
570 //TF1 *g1 = new TF1("g1","gaus",-10,0);
571 //TF1 *g2 = new TF1("g2","gaus",0,10);
572 //h_deltaZ_kb[l]->Fit(g1,"R");
573 //h_deltaZ_kb[l]->Fit(g2,"R+");
574
575 h_deltaZ_tr[l]->GetMaximumBin(xbin,ybin,zbin);
576 skew = h_deltaZ_tr[l]->GetSkewness();
577 p0 = h_deltaZ_tr[l]->GetMaximum();
578 p1 = h_deltaZ_tr[l]->GetMean();
579 p2 = skew > 0 ? 2.0*h_deltaZ_tr[l]->GetRMS() : 0.5*h_deltaZ_tr[l]->GetRMS();
580 p3 = skew < 0 ? 2.0*h_deltaZ_tr[l]->GetRMS() : 0.5*h_deltaZ_tr[l]->GetRMS();
581 p4 = h_deltaZ_tr[l]->GetBinCenter(xbin);
582 f_tr->SetParameters(p0,p1,p2,p3,p4);
583 f_tr->SetLineColor(2);
584 f_tr->FixParameter(4, p4);
585 h_deltaZ_tr[l]->Fit("f1");
586 sigma1_tr[l] = f_tr->GetParameter(2);
587 sigma2_tr[l] = f_tr->GetParameter(3);
588 //TF1 *g3 = new TF1("g3","gaus",-10,0);
589 //TF1 *g4 = new TF1("g4","gaus",0,10);
590 //h_deltaZ_kb[l]->Fit(g3,"R");
591 //h_deltaZ_kb[l]->Fit(g4,"R+");
592 }
593 int yMax = h_deltaZ_kb[l]->GetMaximum() > h_deltaZ_tr[l]->GetMaximum()? h_deltaZ_kb[l]->GetMaximum():h_deltaZ_tr[l]->GetMaximum();
594 h_deltaZ_kb[l]->Draw("same");
595 h_deltaZ_kb[l]->SetLineColor(4);
596 h_deltaZ_tr[l]->Draw("same");
597 h_deltaZ_tr[l]->SetLineColor(2);
598 x_axis = h_deltaZ_kb[l]->GetXaxis();
599 y_axis = h_deltaZ_kb[l]->GetYaxis();
600 x_axis->SetTitle("residual (cm)");
601 x_axis->CenterTitle();
602 //x_axis->SetTitleSize(0.05);
603 //x_axis->SetLabelSize(0.04);
604 //x_axis->SetTitleOffset(1.05);
605 //x_axis->SetLimits(0,80 );
606 y_axis->SetTitle("Entries");
607 y_axis->CenterTitle();
608 //y_axis->SetTitleSize(0.05);
609 //y_axis->SetLabelSize(0.04);
610 //y_axis->SetTitleOffset(1.25);
611 y_axis->SetRangeUser(0,yMax*1.1);
612 TLegend* leg2= new TLegend(0.12 ,0.80,0.33,0.88);
613 leg2->AddEntry(h_deltaZ_kb[l],"k - b","l");
614 leg2->AddEntry(h_deltaZ_tr[l],"#theta - #rho","l");
615 leg2->SetFillColor(kWhite);
616 leg2->SetLineColor(kWhite);
617 leg2->Draw();
618
619 }
620
621 for(int l=0;l<20;l++){
622 cout<<setw(10)<<sigma1_kb[l]<<",";
623 }
624 cout<<endl;
625 for(int l=0;l<20;l++){
626 cout<<setw(10)<<sigma2_kb[l]<<",";
627 }
628
629 for(int l=0;l<20;l++){
630 cout<<setw(10)<<sigma1_tr[l]<<",";
631 }
632 cout<<endl;
633 for(int l=0;l<20;l++){
634 cout<<setw(10)<<sigma2_tr[l]<<",";
635 }
636
637
638 ssname.str("");
639 ssname <<nbin<<"bin_residualZ_layer.eps";
640 sname = ssname.str();
641 name = sname.c_str();
642 c2->SaveAs(name);
643 ssname.str("");
644 ssname <<nbin<<"bin_residualZ_layer.pdf";
645 sname = ssname.str();
646 name = sname.c_str();
647 c2->SaveAs(name);
648 ssname.str("");
649 ssname <<nbin<<"bin_residualZ_layer.png";
650 sname = ssname.str();
651 name = sname.c_str();
652 c2->SaveAs(name);
653
654 ssname.str("");
655 ssname <<nbin<<"bin,residual";
656 sname = ssname.str();
657 name = sname.c_str();
658 TCanvas *c3 = new TCanvas(name,name,1000, 500);
659 c3->Divide(2,1);
660
661 c3->cd(1);
662 h_1_3_kb->Fit("gaus","Q");
663 h_1_3_tr->Fit("gaus","Q");
664 h_1_3_kb->Draw();
665 h_1_3_kb->SetLineColor(4);
666 //h_1_3_kb->SetFillStyle(3345);
667 //h_1_3_kb->SetFillColor(4);
668 h_1_3_kb->GetXaxis()->SetTitle("residual (cm)");
669 h_1_3_kb->GetYaxis()->SetTitle("Entries");
670 h_1_3_tr->Draw("same");
671 h_1_3_tr->SetLineColor(2);
672 //h_1_3_tr->SetFillStyle(3354);
673 //h_1_3_tr->SetFillColor(2);
674 yMax = h_1_3_kb->GetMaximum() > h_1_3_tr->GetMaximum()? h_1_3_kb->GetMaximum() : h_1_3_tr->GetMaximum();
675 h_1_3_kb->GetYaxis()->SetRangeUser(0,yMax*1.1);
676 TLegend* leg31= new TLegend(0.12 ,0.80,0.33,0.88);
677 leg31->AddEntry(h_1_3_kb,"k - b","l");
678 leg31->AddEntry(h_1_3_tr,"#theta - #rho","l");
679 leg31->SetFillColor(kWhite);
680 leg31->SetLineColor(kWhite);
681 leg31->Draw();
682
683 c3->cd(2);
684 TF1 *f1 = new TF1("f1", fitterFunction, -10.0, 10.0, 5);
685 f1->SetParameters(h_21_36_kb->GetMaximum(),h_21_36_kb->GetMaximum()/2,h_21_36_kb->GetMean(),h_21_36_kb->GetRMS()/2.,h_21_36_kb->GetRMS()*2.);
686 f1->SetLineColor(4);
687
688 TF1 *f2 = new TF1("f2", fitterFunction, -10.0, 10.0, 5);
689 f2->SetParameters(h_21_36_tr->GetMaximum(),h_21_36_tr->GetMaximum()/2,h_21_36_tr->GetMean(),h_21_36_tr->GetRMS()/2.,h_21_36_tr->GetRMS()*2.);
690 f2->SetLineColor(2);
691 //double scale = 1.0/h_21_36_kb->Integral();
692 //h_21_36_kb->Scale(scale);
693 h_21_36_kb->Fit("f1");
694 h_21_36_tr->Fit("f2");
695 h_21_36_kb->Draw();
696 h_21_36_kb->SetLineColor(4);
697 h_21_36_kb->GetXaxis()->SetTitle("residual (cm)");
698 h_21_36_kb->GetYaxis()->SetTitle("Entries");
699 h_21_36_tr->Draw("same");
700 h_21_36_tr->SetLineColor(2);
701 yMax = h_21_36_kb->GetMaximum() > h_21_36_tr->GetMaximum()? h_21_36_kb->GetMaximum() : h_21_36_tr->GetMaximum();
702 h_21_36_kb->GetYaxis()->SetRangeUser(0,yMax*1.1);
703 TLegend* leg32= new TLegend(0.12 ,0.80,0.33,0.88);
704 leg32->AddEntry(h_21_36_kb,"k - b","l");
705 leg32->AddEntry(h_21_36_tr,"#theta - #rho","l");
706 leg32->SetFillColor(kWhite);
707 leg32->SetLineColor(kWhite);
708 leg32->Draw();
709
710 ssname.str("");
711 ssname <<nbin<<"bin_residualZ.eps";
712 sname = ssname.str();
713 name = sname.c_str();
714 c3->SaveAs(name);
715 ssname.str("");
716 ssname <<nbin<<"bin_residualZ.pdf";
717 sname = ssname.str();
718 name = sname.c_str();
719 c3->SaveAs(name);
720 ssname.str("");
721 ssname <<nbin<<"bin_residualZ.png";
722 sname = ssname.str();
723 name = sname.c_str();
724 c3->SaveAs(name);
725
726
727 ssname.str("");
728 ssname <<nbin<<"bin,tanl,dz";
729 sname = ssname.str();
730 name = sname.c_str();
731 TCanvas *c4 = new TCanvas(name,name,1000, 500);
732 c4->Divide(2,1);
733 c4->cd(1);
734 h_tanl_kb->Fit("gaus","Q");
735 h_tanl_tr->Fit("gaus","Q");
736 h_tanl_kb->Draw();
737 h_tanl_kb->SetLineColor(4);
738 h_tanl_kb->GetXaxis()->SetTitle("residual (cm)");
739 h_tanl_kb->GetYaxis()->SetTitle("Entries");
740 h_tanl_tr->Draw("same");
741 h_tanl_tr->SetLineColor(2);
742 yMax = h_tanl_kb->GetMaximum() > h_tanl_tr->GetMaximum()? h_tanl_kb->GetMaximum() : h_tanl_tr->GetMaximum();
743 h_tanl_kb->GetYaxis()->SetRangeUser(0,yMax*1.1);
744 TLegend* leg41= new TLegend(0.12 ,0.80,0.33,0.88);
745 leg41->AddEntry(h_tanl_kb,"k - b","l");
746 leg41->AddEntry(h_tanl_tr,"#theta - #rho","l");
747 leg41->SetFillColor(kWhite);
748 leg41->SetLineColor(kWhite);
749 leg41->Draw();
750 c4->cd(2);
751 h_dz_kb->Fit("gaus","Q");
752 h_dz_tr->Fit("gaus","Q");
753 h_dz_kb->Draw();
754 h_dz_kb->SetLineColor(4);
755 h_dz_kb->GetXaxis()->SetTitle("residual (cm)");
756 h_dz_kb->GetYaxis()->SetTitle("Entries");
757 h_dz_tr->Draw("same");
758 h_dz_tr->SetLineColor(2);
759 yMax = h_dz_kb->GetMaximum() > h_dz_tr->GetMaximum()? h_dz_kb->GetMaximum() : h_dz_tr->GetMaximum();
760 h_dz_kb->GetYaxis()->SetRangeUser(0,yMax*1.1);
761 TLegend* leg42= new TLegend(0.12 ,0.80,0.33,0.88);
762 leg42->AddEntry(h_dz_kb,"k - b","l");
763 leg42->AddEntry(h_dz_tr,"#theta - #rho","l");
764 leg42->SetFillColor(kWhite);
765 leg42->SetLineColor(kWhite);
766 leg42->Draw();
767 ssname.str("");
768 ssname <<nbin<<"bin_residual_tanl_dz.pdf";
769 sname = ssname.str();
770 name = sname.c_str();
771 c4->SaveAs(name);
772 ssname.str("");
773 ssname <<nbin<<"bin_residual_tanl_dz.png";
774 sname = ssname.str();
775 name = sname.c_str();
776 c4->SaveAs(name);
777 ssname.str("");
778 ssname <<nbin<<"bin_residual_tanl_dz.eps";
779 sname = ssname.str();
780 name = sname.c_str();
781 c4->SaveAs(name);
782
783
784 double sigma_tanl = h_tanl_kb->GetFunction("gaus")->GetParameter(2);
785 double sigma_tanl_err = h_tanl_kb->GetFunction("gaus")->GetParError(2);
786 double sigma_dz = h_dz_kb->GetFunction("gaus")->GetParameter(2);
787 double sigma_dz_err = h_dz_kb->GetFunction("gaus")->GetParError(2);
788 gr_tanl_kb->SetPoint(ibin,nbin,sigma_tanl);
789 gr_tanl_kb->SetPointError(ibin,0,sigma_tanl_err);
790 gr_dz_kb->SetPoint(ibin,nbin,sigma_dz);
791 gr_dz_kb->SetPointError(ibin,0,sigma_dz_err);
792
793 sigma_tanl = h_tanl_tr->GetFunction("gaus")->GetParameter(2);
794 sigma_tanl_err = h_tanl_tr->GetFunction("gaus")->GetParError(2);
795 sigma_dz = h_dz_tr->GetFunction("gaus")->GetParameter(2);
796 sigma_dz_err = h_dz_tr->GetFunction("gaus")->GetParError(2);
797 gr_tanl_tr->SetPoint(ibin,nbin,sigma_tanl);
798 gr_tanl_tr->SetPointError(ibin,0,sigma_tanl_err);
799 gr_dz_tr->SetPoint(ibin,nbin,sigma_dz);
800 gr_dz_tr->SetPointError(ibin,0,sigma_dz_err);
801//*/
802
803///*
804 double C1 =f1->GetParameter(0);
805 double C2 =f1->GetParameter(1);
806 double sigma1 = f1->GetParameter(3);
807 double sigma2 = f1->GetParameter(4);
808 double sigma_err1 = f1->GetParError(3);
809 double sigma_err2 = f1->GetParError(4);
810 double sigma_1_3 = h_1_3_kb->GetFunction("gaus")->GetParameter(2);
811 double sigma_err_1_3 = h_1_3_kb->GetFunction("gaus")->GetParError(2);
812 double sigma_21_36 = sqrt(C1/(C1+C2)*sigma1*sigma1 + C2/(C1+C2)*sigma2*sigma2);
813 double sigma_err_21_36 = sqrt(C1/(C1+C2)*sigma_err1*sigma_err1 + C2/(C1+C2)*sigma_err2*sigma_err2);
814 double RMS_1_3 = h_1_3_kb->GetRMS();
815 double RMS_err_1_3 = h_1_3_kb->GetRMSError();
816 double RMS_21_36 = h_21_36_kb->GetRMS();
817 double RMS_err_21_36 = h_21_36_kb->GetRMSError();
818 //cout<<nbin<<" "<<f1->GetParameter(2)<<" "<<f1->GetParameter(5)<<" "<<f1->GetParameter(2)*f1->GetParameter(6)+f1->GetParameter(5)*(1-f1->GetParameter(6))<<endl;
819 cout<<"kb sigma "<<setw(6)<<nbin<<setw(12)<<sigma_1_3<<setw(12)<<sigma_err_1_3<<setw(12)<<sigma_21_36<<setw(12)<<sigma_err_21_36<<endl;
820 cout<<"kb RMS "<<setw(6)<<nbin<<setw(12)<<RMS_1_3<<setw(12)<<RMS_err_1_3<<setw(12)<<RMS_21_36<<setw(12)<<RMS_err_21_36<<endl;
821 gr_sigma_1_3_kb->SetPoint(ibin,nbin,sigma_1_3);
822 gr_sigma_1_3_kb->SetPointError(ibin,0,sigma_err_1_3);
823 gr_sigma_21_36_kb->SetPoint(ibin,nbin,sigma_21_36);
824 gr_sigma_21_36_kb->SetPointError(ibin,0,sigma_err_21_36);
825
826 gr_RMS_1_3_kb->SetPoint(ibin,nbin,RMS_1_3);
827 gr_RMS_1_3_kb->SetPointError(ibin,0,RMS_err_1_3);
828 gr_RMS_21_36_kb->SetPoint(ibin,nbin,RMS_21_36);
829 gr_RMS_21_36_kb->SetPointError(ibin,0,RMS_err_21_36);
830
831 gr_sigma_tr->SetPoint(ibin,nbin,sigma_21_36);
832 gr_sigma1_tr->SetPoint(ibin,nbin,sigma1);
833 gr_sigma2_tr->SetPoint(ibin,nbin,sigma2);
834
835 C1 =f2->GetParameter(0);
836 C2 =f2->GetParameter(1);
837 sigma1 = f2->GetParameter(3);
838 sigma2 = f2->GetParameter(4);
839 sigma_err1 = f2->GetParError(3);
840 sigma_err2 = f2->GetParError(4);
841 sigma_1_3 = h_1_3_tr->GetFunction("gaus")->GetParameter(2);
842 sigma_err_1_3 = h_1_3_tr->GetFunction("gaus")->GetParError(2);
843 sigma_21_36 = sqrt(C1/(C1+C2)*sigma1*sigma1 + C2/(C1+C2)*sigma2*sigma2);
844 sigma_err_21_36 = sqrt(C1/(C1+C2)*sigma_err1*sigma_err1 + C2/(C1+C2)*sigma_err2*sigma_err2);
845 RMS_1_3 = h_1_3_tr->GetRMS();
846 RMS_err_1_3 = h_1_3_tr->GetRMSError();
847 RMS_21_36 = h_21_36_tr->GetRMS();
848 RMS_err_21_36 = h_21_36_tr->GetRMSError();
849 //cout<<nbin<<" "<<f2->GetParameter(2)<<" "<<f2->GetParameter(5)<<" "<<f2->GetParameter(2)*f2->GetParameter(6)+f2->GetParameter(5)*(1-f2->GetParameter(6))<<endl;
850 cout<<"tr sigma "<<setw(6)<<nbin<<setw(12)<<sigma_1_3<<setw(12)<<sigma_err_1_3<<setw(12)<<sigma_21_36<<setw(12)<<sigma_err_21_36<<endl;
851 cout<<"tr RMS "<<setw(6)<<nbin<<setw(12)<<RMS_1_3<<setw(12)<<RMS_err_1_3<<setw(12)<<RMS_21_36<<setw(12)<<RMS_err_21_36<<endl;
852 gr_sigma_1_3_tr->SetPoint(ibin,nbin,sigma_1_3);
853 gr_sigma_1_3_tr->SetPointError(ibin,0,sigma_err_1_3);
854 gr_sigma_21_36_tr->SetPoint(ibin,nbin,sigma_21_36);
855 gr_sigma_21_36_tr->SetPointError(ibin,0,sigma_err_21_36);
856
857 gr_RMS_1_3_tr->SetPoint(ibin,nbin,RMS_1_3);
858 gr_RMS_1_3_tr->SetPointError(ibin,0,RMS_err_1_3);
859 gr_RMS_21_36_tr->SetPoint(ibin,nbin,RMS_21_36);
860 gr_RMS_21_36_tr->SetPointError(ibin,0,RMS_err_21_36);
861//*/
862 gr_sigma_tr->SetPoint(ibin,nbin,sigma_21_36);
863 gr_sigma1_tr->SetPoint(ibin,nbin,sigma1);
864 gr_sigma2_tr->SetPoint(ibin,nbin,sigma2);
865 //cout<<setw(6)<<nbin<<setw(12)<<c1<<setw(12)<<c2<<setw(12)<<f1->GetParameter(3)<<setw(12)<<sigma1<<setw(12)<<sigma2<<setw(12)<<sig ma<<endl;
866 ibin++;
867 int bin;
868 if(nbin<500) bin = 50;
869 else if(nbin<1000) bin =100;
870 else if(nbin<2000)bin =200;
871 else if(nbin<5000)bin =500;
872 else bin = 1000;
873 nbin += bin;
874 }
875///*
876 TCanvas *c5 = new TCanvas("sigma_deltaZ","sigma_deltaZ",1000,1000);
877 //c5->Divide(2,1);
878 //c5->cd(1);
879 gr_sigma_1_3_kb->GetXaxis()->SetTitle("number of bins");
880 gr_sigma_1_3_kb->GetYaxis()->SetTitle("#sigma (cm)");
881 gr_sigma_1_3_kb->SetMinimum(0.0);
882 gr_sigma_1_3_kb->SetMaximum( 5.0);
883 gr_sigma_1_3_kb->Draw("APsame");
884 gr_sigma_1_3_kb->SetMarkerStyle(20);
885 gr_sigma_1_3_kb->SetMarkerSize(1.0);
886 gr_sigma_1_3_kb->SetMarkerColor(4);
887 gr_sigma_21_36_kb->Draw(" Psame");
888 gr_sigma_21_36_kb->SetMarkerStyle(21);
889 gr_sigma_21_36_kb->SetMarkerSize(1.0);
890 gr_sigma_21_36_kb->SetMarkerColor(4);
891
892 gr_RMS_1_3_kb->Draw("Psame");
893 gr_RMS_1_3_kb->SetMarkerStyle(24);
894 gr_RMS_1_3_kb->SetMarkerSize(1.0);
895 gr_RMS_1_3_kb->SetMarkerColor(4);
896 gr_RMS_21_36_kb->Draw(" Psame");
897 gr_RMS_21_36_kb->SetMarkerStyle(25);
898 gr_RMS_21_36_kb->SetMarkerSize(1.0);
899 gr_RMS_21_36_kb->SetMarkerColor(4);
900
901 //gr_sigma_1_3_tr->SetMinimum(0.0);
902 //gr_sigma_1_3_tr->SetMaximum(10.0);
903 //gr_sigma_1_3_tr->GetXaxis()->SetTitle("number of bins");
904 //gr_sigma_1_3_tr->GetYaxis()->SetTitle("#sigma (cm)");
905 gr_sigma_1_3_tr->Draw("Psame");
906 gr_sigma_1_3_tr->SetMarkerStyle(20);
907 gr_sigma_1_3_tr->SetMarkerSize(1.0);
908 gr_sigma_1_3_tr->SetMarkerColor(2);
909 gr_sigma_21_36_tr->Draw(" Psame");
910 gr_sigma_21_36_tr->SetMarkerStyle(21);
911 gr_sigma_21_36_tr->SetMarkerSize(1.0);
912 gr_sigma_21_36_tr->SetMarkerColor(2);
913
914 gr_RMS_1_3_tr->Draw("Psame");
915 gr_RMS_1_3_tr->SetMarkerStyle(24);
916 gr_RMS_1_3_tr->SetMarkerSize(1.0);
917 gr_RMS_1_3_tr->SetMarkerColor(2);
918 gr_RMS_21_36_tr->Draw(" Psame");
919 gr_RMS_21_36_tr->SetMarkerStyle(25);
920 gr_RMS_21_36_tr->SetMarkerSize(1.0);
921 gr_RMS_21_36_tr->SetMarkerColor(2);
922
923 //TLegend* leg52= new TLegend(0.55 ,0.60,0.88,0.88);
924 TLegend* leg5= new TLegend(0.55 ,0.60,0.88,0.88);
925 leg5->AddEntry(gr_sigma_1_3_kb," 1~3 layer sigma, k-b","p");
926 leg5->AddEntry(gr_sigma_1_3_tr," 1~3 layer sigma, #rho-#theta","p");
927 leg5->AddEntry(gr_sigma_21_36_kb,"21~36 layer sigma, k-b","p");
928 leg5->AddEntry(gr_sigma_21_36_tr,"21~36 layer sigma, #rho-#theta","p");
929 leg5->AddEntry(gr_RMS_1_3_kb," 1~3 layer RMS, k-b","p");
930 leg5->AddEntry(gr_RMS_1_3_tr," 1~3 layer RMS, #rho-#theta","p");
931 leg5->AddEntry(gr_RMS_21_36_kb,"21~36 layer RMS, k-b","p");
932 leg5->AddEntry(gr_RMS_21_36_tr,"21~36 layer RMS, #rho-#theta","p");
933 leg5->SetFillColor(kWhite);
934 leg5->SetLineColor(kWhite);
935 leg5->Draw();
936 c5->SaveAs("sigma_deltaZ.eps");
937 c5->SaveAs("sigma_deltaZ.png");
938 c5->SaveAs("sigma_deltaZ.pdf");
939
940
941 TCanvas *c6 = new TCanvas("sigma_tanl,dz","sigma_tanl,dz",1000, 500);
942 c6->Divide(2,1);
943 c6->cd(1);
944 //gr_tanl_kb->SetMinimum(0.0);
945 //gr_tanl_kb->SetMaximum(1.0);
946 gr_tanl_kb->GetXaxis()->SetTitle("number of bins");
947 gr_tanl_kb->GetYaxis()->SetTitle("#sigma (cm)");
948 gr_tanl_kb->SetTitle("tan #lambda");
949 gr_tanl_kb->Draw("APsame");
950 gr_tanl_kb->SetMarkerStyle(24);
951 gr_tanl_kb->SetMarkerSize(1.0);
952 gr_tanl_kb->SetMarkerColor(4);
953 gr_tanl_kb->SetLineColor(4);
954 gr_tanl_tr->Draw(" Psame");
955 gr_tanl_tr->SetMarkerStyle(24);
956 gr_tanl_tr->SetMarkerSize(1.0);
957 gr_tanl_tr->SetMarkerColor(2);
958 gr_tanl_tr->SetLineColor(2);
959 TLegend* leg61= new TLegend(0.55 ,0.75,0.88,0.88);
960 leg61->AddEntry(gr_tanl_kb,"k - b","l");
961 leg61->AddEntry(gr_tanl_tr,"#theta - #rho","l");
962 leg61->SetFillColor(kWhite);
963 leg61->SetLineColor(kWhite);
964 leg61->Draw();
965 c6->cd(2);
966 //gr_tanl_kb->SetMinimum(0.0);
967 //gr_tanl_kb->SetMaximum(1.0);
968 gr_dz_kb->GetXaxis()->SetTitle("number of bins");
969 gr_dz_kb->GetYaxis()->SetTitle("#sigma (cm)");
970 gr_dz_kb->SetTitle("dz");
971 gr_dz_kb->Draw("APsame");
972 gr_dz_kb->SetMarkerStyle(24);
973 gr_dz_kb->SetMarkerSize(1.0);
974 gr_dz_kb->SetMarkerColor(4);
975 gr_dz_kb->SetLineColor(4);
976 gr_dz_tr->Draw(" Psame");
977 gr_dz_tr->SetMarkerStyle(24);
978 gr_dz_tr->SetMarkerSize(1.0);
979 gr_dz_tr->SetMarkerColor(2);
980 gr_dz_tr->SetLineColor(2);
981 TLegend* leg62= new TLegend(0.55 ,0.75,0.88,0.88);
982 leg62->AddEntry(gr_dz_kb,"k - b","l");
983 leg62->AddEntry(gr_dz_tr,"#theta - #rho","l");
984 leg62->SetFillColor(kWhite);
985 leg62->SetLineColor(kWhite);
986 leg62->Draw();
987 c6->SaveAs("sigma_tanl.eps");
988 c6->SaveAs("sigma_tanl.png");
989 c6->SaveAs("sigma_tanl.pdf");
990
991
992 TCanvas *c7 = new TCanvas("c7","c7",1000, 500);
993 c7->Divide(2,1);
994 c7->cd(1);
995 gr_sigma_kb->SetMinimum(0.0);
996 gr_sigma_kb->SetMaximum(5.0);
997 gr_sigma_kb->GetXaxis()->SetTitle("number of bins");
998 gr_sigma_kb->GetYaxis()->SetTitle("#sigma (cm)");
999 gr_sigma_kb->Draw("APsame");
1000 gr_sigma_kb->SetMarkerStyle(24);
1001 gr_sigma_kb->SetMarkerSize(1.0);
1002 gr_sigma_kb->SetMarkerColor(2);
1003 gr_sigma1_kb->Draw("Psame");
1004 gr_sigma1_kb->SetMarkerStyle(24);
1005 gr_sigma1_kb->SetMarkerSize(1.0);
1006 gr_sigma1_kb->SetMarkerColor(3);
1007 gr_sigma2_kb->Draw("Psame");
1008 gr_sigma2_kb->SetMarkerStyle(24);
1009 gr_sigma2_kb->SetMarkerSize(1.0);
1010 gr_sigma2_kb->SetMarkerColor(4);
1011 TLegend* leg41= new TLegend(0.70 ,0.75,0.88,0.88);
1012 leg41->AddEntry(gr_sigma1_kb,"#sigma_{1}","p");
1013 leg41->AddEntry(gr_sigma2_kb,"#sigma_{2}","p");
1014 leg41->AddEntry(gr_sigma_kb ,"#sigma_{mean}","p");
1015 leg41->SetFillColor(kWhite);
1016 leg41->SetLineColor(kWhite);
1017 leg41->Draw();
1018 c7->cd(1);
1019 gr_sigma_tr->SetMinimum(0.0);
1020 gr_sigma_tr->SetMaximum(5.0);
1021 gr_sigma_tr->GetXaxis()->SetTitle("number of bins");
1022 gr_sigma_tr->GetYaxis()->SetTitle("#sigma (cm)");
1023 gr_sigma_tr->Draw("APsame");
1024 gr_sigma_tr->SetMarkerStyle(24);
1025 gr_sigma_tr->SetMarkerSize(1.0);
1026 gr_sigma_tr->SetMarkerColor(2);
1027 gr_sigma1_tr->Draw("Psame");
1028 gr_sigma1_tr->SetMarkerStyle(24);
1029 gr_sigma1_tr->SetMarkerSize(1.0);
1030 gr_sigma1_tr->SetMarkerColor(3);
1031 gr_sigma2_tr->Draw("Psame");
1032 gr_sigma2_tr->SetMarkerStyle(24);
1033 gr_sigma2_tr->SetMarkerSize(1.0);
1034 gr_sigma2_tr->SetMarkerColor(4);
1035 TLegend* leg41= new TLegend(0.70 ,0.75,0.88,0.88);
1036 leg41->AddEntry(gr_sigma1_tr,"#sigma_{1}","p");
1037 leg41->AddEntry(gr_sigma2_tr,"#sigma_{2}","p");
1038 leg41->AddEntry(gr_sigma_tr ,"#sigma_{mean}","p");
1039 leg41->SetFillColor(kWhite);
1040 leg41->SetLineColor(kWhite);
1041 leg41->Draw();
1042 ssname.str("");
1043 ssname <<"bin_sigma.pdf";
1044 sname = ssname.str();
1045 name = sname.c_str();
1046 c7->SaveAs(name);
1047 ssname.str("");
1048 ssname <<"bin_sigma.png";
1049 sname = ssname.str();
1050 name = sname.c_str();
1051 c7->SaveAs(name);
1052 ssname.str("");
1053 ssname <<"bin_sigma.eps";
1054 sname = ssname.str();
1055 name = sname.c_str();
1056 c7->SaveAs(name);
1057//*/
1058}
1059///*
1060void fill_histogram(double s, double z, double x_max, int nbin, TH2D* hough)
1061{
1062
1063 double dx = 2*x_max/nbin;
1064 double x = -x_max - 0.5*dx;
1065 double y = -1*s*x+z;
1066 double dy = -1*s*dx;
1067 for(int i =0; i<nbin; i++){
1068 x += dx;
1069 y += dy;
1070 hough->Fill(x,y);
1071 }
1072}
1073//*/
1074///*
1075void fill_histogram2(double s, double z, double x_max, int nbin, TH2D* hough)
1076{
1077
1078 double step = 2*x_max/nbin;
1079 double theta = -x_max - 0.5*step;
1080 double rho = 0;
1081 for(int i =0; i<nbin; i++){
1082 theta += step;
1083 rho = s*cos(theta) + z*sin(theta);
1084 hough->Fill(theta,rho);
1085 }
1086}
1087//*/
1088
1089double fitterFunction(double *x, double *par)
1090{
1091 Double_t arg1 = 0;
1092 if (par[3] != 0) arg1 = (x[0] - par[2])/par[3];
1093 Double_t fitval1 = par[0]*TMath::Exp(-0.5*arg1*arg1);
1094
1095 Double_t arg2 = 0;
1096 if (par[4] != 0) arg2 = (x[0] - par[2])/par[4];
1097 Double_t fitval2 = par[1]*TMath::Exp(-0.5*arg2*arg2);
1098
1099 Double_t pdfValue = fitval1+fitval2;
1100 return pdfValue;
1101}
1102double fitterFunction2(double *x, double *par)
1103{
1104 Double_t arg = 0;
1105 Double_t fitval = 0;
1106 if(x[0]<par[4]){
1107 if (par[2] != 0) arg = (x[0] - par[1])/par[2];
1108 fitval = par[0]*TMath::Exp(-0.5*arg*arg);
1109 return fitval;
1110 }else{
1111 if (par[3] != 0) arg = (x[0] - par[1])/par[3];
1112 fitval = par[0]*TMath::Exp(-0.5*arg*arg);
1113 return fitval;
1114 }
1115
1116}
1117/*
1118double find_cut(double percent, TH1D* histogram)
1119{
1120 int nbin = histogram->GetNbinsX();
1121 double sum = 0;
1122 double total = histogram->Integral(1,nbin);
1123 for(int ibin=1;ibin<nbin/2;ibin++){
1124 sum = histogram->Integral(nbin/2-ibin+1,nbin/2+ibin);
1125 if(sum/total>percent){
1126 cut= histogram->GetBinCenter(nbin/2+ibin);
1127 break;
1128 }
1129 }
1130 return cut;
1131}
1132Double_t fitf(Double_t *x, Double_t *par)
1133{
1134 Double_t arg = 0;
1135 if (par[2] != 0) arg = (x[0] - par[1])/par[2];
1136 Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg);
1137 return fitval;
1138}
1139double intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
1140{
1141 const double M_PI = 3.1415926;
1142 double phi;
1143 double phi_center = atan2(y_center,x_center);
1144 double r_center = sqrt(x_center*x_center+y_center*y_center);
1145 if(charge==0)return 9999;
1146 while(phi_center<0) phi_center += 2*M_PI;
1147 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
1148 if(r_center<r_cylinder/2) return -9999;
1149 double dphi = acos( r_cylinder/(2*r_center) );
1150 if(charge<0) phi = phi_center + dphi;
1151 else phi = phi_center - dphi;
1152 while(phi<0) phi += 2*M_PI;
1153 while(phi>2*M_PI) phi -= 2*M_PI;
1154 return phi;
1155}
1156
1157void fill_histogram(double x, double y, double r, int nbin, TH2D* hough)
1158{
1159 const double M_PI = 3.1415926;
1160 double U = 2*x/(x*x+y*y-r*r);
1161 double V = 2*y/(x*x+y*y-r*r);
1162 double R = 2*r/(x*x+y*y-r*r);
1163 double delta_phi = 2.*M_PI/nbin;
1164 double phi = -delta_phi/2;
1165 for(int i =0; i<nbin; i++){
1166 phi += delta_phi;
1167 double u = U + R*cos(phi);
1168 double v = V + R*sin(phi);
1169 double normal = (v-V)/(u-U);
1170 double k = -1./normal;
1171 double b = v - k*u;
1172 double x_cross = -b/(k+1/k);
1173 double y_cross = b/(1+k*k);
1174
1175 double rho=sqrt(x_cross*x_cross+y_cross*y_cross);
1176 double theta=atan2(y_cross,x_cross);
1177 if(theta<0) {
1178 theta=theta+M_PI;
1179 rho=-rho;
1180 }
1181 if( normal ==0 && x>0) {
1182 rho = fabs(x);
1183 theta = 0;
1184 }
1185 if( normal ==0 && x<0) {
1186 rho = -fabs(x);
1187 theta = M_PI;
1188 }
1189 hough->Fill(theta,rho);
1190 }
1191}
1192//*/
f1
TFile * f1
Definition: DataBase/tau_mode.c:4
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
nentries
Int_t nentries
Definition: DataBase/tau_mode.c:8
arg
double arg(const EvtComplex &c)
Definition: EvtComplex.hh:227
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85
s
XmlRpcServer s
Definition: HelloServer.cpp:11
tan
double tan(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:216
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
cut
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared cut
Definition: KarFin.h:27
M_PI
#define M_PI
Definition: TConstant.h:4
fitterFunction2
double fitterFunction2(double *x, double *par)
Definition: bin_deltaZ.cxx:1102
fill_histogram2
void fill_histogram2(double s, double z, double x_max, int nbin, TH2D *hough)
Definition: bin_deltaZ.cxx:1075
fill_histogram
void fill_histogram(double s, double z, double x_max, int nbin, TH2D *hough)
Definition: bin_deltaZ.cxx:1060
fitterFunction
double fitterFunction(double *x, double *par)
Definition: bin_deltaZ.cxx:1089
bin_deltaZ
int bin_deltaZ()
Definition: bin_deltaZ.cxx:1