2{
3 const double M_PI = 3.1415926;
4
5
6
7
8
9
10
11
12
13 gROOT->Reset();
14 gStyle->SetOptFit(0111);
15 TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hough_hist_test.root");
16 if (!f) {
17 f = new TFile("hough_hist_test.root");
18 }
19 TTree* tree;
20 f->GetObject("hit",tree);
21
22
23 Int_t hit_run;
24 Int_t hit_evt;
25 Int_t hit_nhit;
26 Int_t hit_hitid[1000];
27 Int_t hit_layer[1000];
28 Int_t hit_wire[1000];
29 Double_t hit_x[1000];
30 Double_t hit_y[1000];
31 Double_t hit_z[1000];
32 Double_t hit_driftdist[10000];
33 Double_t hit_drifttime[10000];
34 Int_t hit_flag[1000];
35 Double_t hit_truth_x[1000];
36 Double_t hit_truth_y[1000];
37 Double_t hit_truth_z[1000];
38 Double_t hit_truth_drift[1000];
39 Int_t hit_truth_ambig[1000];
40
41
42 hit->SetBranchAddress("hit_run",&hit_run);
43 hit->SetBranchAddress("hit_evt",&hit_evt);
44 hit->SetBranchAddress("hit_nhit",&hit_nhit);
45 hit->SetBranchAddress("hit_hitid",hit_hitid);
46 hit->SetBranchAddress("hit_layer",hit_layer);
47 hit->SetBranchAddress("hit_wire",hit_wire);
48 hit->SetBranchAddress("hit_x",hit_x);
49 hit->SetBranchAddress("hit_y",hit_y);
50 hit->SetBranchAddress("hit_z",hit_z);
51 hit->SetBranchAddress("hit_driftdist",hit_driftdist);
52 hit->SetBranchAddress("hit_drifttime",hit_drifttime);
53 hit->SetBranchAddress("hit_flag",hit_flag);
54 hit->SetBranchAddress("hit_truth_x",hit_truth_x);
55 hit->SetBranchAddress("hit_truth_y",hit_truth_y);
56 hit->SetBranchAddress("hit_truth_z",hit_truth_z);
57 hit->SetBranchAddress("hit_truth_drift",hit_truth_drift);
58 hit->SetBranchAddress("hit_truth_ambig",hit_truth_ambig);
59
60
61
62
63
64
65 TAxis *x_axis,*y_axis;
66 stringstream ssname;
67 string sname;
68 const char * name;
69
70 int ibin = 0;
71 int nbin = 100;
72 while(nbin <= 100 ){
73 int x_bin = nbin;
74 int y_bin = nbin;
76 double y_max = 0.1;
77
78 ssname.str("");
79 ssname <<nbin<<"bin, all layers";
80 sname = ssname.str();
81 name = sname.c_str();
82 TCanvas *c_hit ;
83
84
85 int event =1 ;
86 Long64_t
nentries = hit->GetEntries();
87 Long64_t nbytes = 0;
88
89
90 for (Long64_t i=event; i<event+1;i++)
91 {
92 nbytes += hit->GetEntry(i);
93 TH2D* hough = new TH2D("hough","hough",x_bin,0,x_max,y_bin,-y_max,y_max);
94 TH2D *houghhit[1000];
95 int nhit(0),ihit(0);
96
97 for(int j=0;j<hit_nhit;j++){
98 if(hit_flag[j] !=0)continue;
99 nhit++;
100
101
102
103 ssname.str("");
104 ssname <<"layer:"<<hit_layer[j]<<" ,wire"<<hit_wire[j];
105 sname = ssname.str();
106 name = sname.c_str();
107 houghhit[ihit] = new TH2D(name,name,x_bin,0,x_max,y_bin,-y_max,y_max);
108 fill_histogram(hit_x[j],hit_y[j],hit_driftdist[j],x_bin*2,houghhit[ihit]);
109 hough->Add(houghhit[ihit]);
110
111
112
113
114 ihit++;
115 }
116
117 TCanvas *c_map = new TCanvas("houghspace","houghspace",1000,500 );
118 c_map->Divide(2,1);
119 c_map->cd(1);
120 hough->Draw();
121 c_map->cd(2);
122 hough->Draw("LEGO2Z");
123 ssname.str("");
124 ssname <<"hough_map.png";
125 sname = ssname.str();
126 name = sname.c_str();
127
128 ssname.str("");
129 ssname <<"hough_map.pdf";
130 sname = ssname.str();
131 name = sname.c_str();
132
133 ssname.str("");
134 ssname <<"hough_map.eps";
135 sname = ssname.str();
136 name = sname.c_str();
137
138
139 int binx,biny,binz;
140 hough->GetMaximumBin(binx,biny,binz);
141 double theta = hough->GetXaxis()->GetBinCenter(binx);
142 double rho = hough->GetYaxis()->GetBinCenter(biny);
143
144
145 double Rc = 1./fabs(rho);
146 double Xc =
cos(theta)/rho;
147 double Yc =
sin(theta)/rho;
148 cout<<Xc<<" "<<Yc<<endl;
149 double x1(999),y1(999),x2(-999),y2(-999);
150 double alpha_max(0);
151 for(int j=0;j<hit_nhit;j++){
152 if(hit_flag[j] !=0)continue;
153 x1 = hit_x[j] < x1 ? hit_x[j] : x1;
154 y1 = hit_y[j] < y1 ? hit_y[j] : y1;
155 x2 = hit_x[j] > x2 ? hit_x[j] : x2;
156 y2 = hit_y[j] > y2 ? hit_y[j] : y2;
157 double alpha =
M_PI-atan2(hit_y[j]-Yc,hit_x[j]-Xc)+atan2(Yc,Xc);
161
162 }
163
164 TEllipse *e = new TEllipse(Xc,Yc,Rc,Rc);
165
166
167 TGraph *gr_xy = new TGraph();
168 TGraph *gr_MC = new TGraph();
169 TGraph *gr_O = new TGraph();
170 TGraph *gr_C = new TGraph();
171 gr_O->SetPoint(0,0,0);
172 gr_C->SetPoint(0,Xc,Yc);
173
174
175 int point(0);
176 for(int j=0;j<hit_nhit;j++){
177 if(hit_flag[j] !=0)continue;
178 gr_xy->SetPoint(point,hit_x[j],hit_y[j]);
179 gr_MC->SetPoint(point,hit_truth_x[j],hit_truth_y[j]);
180
181
182
183
184
185
186
187
188 point++;
189 }
190 ssname.str("");
191 ssname <<"event "<<i;
192 sname = ssname.str();
193 name = sname.c_str();
194 TCanvas *c_xy = new TCanvas(name,name,1000,1000 );
195 gr_xy->Draw("APsame");
196 double x_min(-81),x_max(81),y_min(-81),y_max(81);
197
198 double scale = (x2-x1)>(y2-y1)?(x2-x1):(y2-y1);
199 x_min = (x2+x1)/2.-scale*0.6;
200 x_max = (x2+x1)/2.+scale*0.6;
201 y_min = (y2+y1)/2.-scale*0.6;
202 y_max = (y2+y1)/2.+scale*0.6;
203
204
205
206
207
208
209
210
211 gr_xy->GetXaxis()->SetLimits(x_min,x_max);
212 gr_xy->SetMinimum(y_min);
213 gr_xy->SetMaximum(y_max);
214 gr_xy->SetMarkerStyle(5 );
215 gr_xy->SetMarkerSize(0.7);
216 gr_xy->SetMarkerColor(1);
217 gr_MC->Draw("Psame");
218 gr_MC->SetMarkerStyle(20);
219 gr_MC->SetMarkerSize(0.7);
220 gr_MC->SetMarkerColor(2);
221 gr_O->Draw("Psame");
222 gr_O->SetMarkerStyle(21);
223 gr_O->SetMarkerSize(1.5);
224 gr_O->SetMarkerColor(1);
225 gr_C->Draw("Psame");
226 gr_C->SetMarkerStyle(20);
227 gr_C->SetMarkerSize(1.5);
228 gr_C->SetMarkerColor(1);
229 e->Draw("same");
230 e->SetFillStyle(0);
231
232 for(int j=0;j<hit_nhit;j++){
233 if(hit_flag[j] !=0)continue;
234 TEllipse *exy = new TEllipse(hit_x[j],hit_y[j],hit_driftdist[j],hit_driftdist[j]);
235 exy->Draw("Csame");
236 exy->SetFillStyle(0);
237 }
238
239 for(int ihit=0;ihit<nhit;ihit++){
240
241
242 }
243 ssname.str("");
244 ssname <<"2d_track.eps";
245 sname = ssname.str();
246 name = sname.c_str();
247
248 ssname.str("");
249 ssname <<"2d_track.pdf";
250 sname = ssname.str();
251 name = sname.c_str();
252
253 ssname.str("");
254 ssname <<"2d_track.png";
255 sname = ssname.str();
256 name = sname.c_str();
257
258 TGraph *gr_uv= new TGraph();
259 int point(0);
260 for(int j=0;j<hit_nhit;j++){
261 if(hit_flag[j] !=0)continue;
262 double u= 2*hit_x[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
263 double v= 2*hit_y[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
264 double w= 2*hit_driftdist[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
265 gr_uv->SetPoint(point,u,
v);
266
267
268
269 point++;
270 }
271 ssname.str("");
272 ssname <<"uv "<<i;
273 sname = ssname.str();
274 name = sname.c_str();
275 TCanvas *c_uv = new TCanvas(name,name,1000,1000 );
276 gr_uv->Draw("APsame");
277 for(int j=0;j<hit_nhit;j++){
278 if(hit_flag[j] !=0)continue;
279 double u= 2*hit_x[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
280 double v= 2*hit_y[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
281 double w= 2*hit_driftdist[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
282
283 TEllipse *euv =
new TEllipse(u,
v,w,w);
284 euv->Draw("Csame");
285 euv->SetFillStyle(0);
286
287 }
288
289 }
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422 ibin++;
424 if(nbin<500)
bin = 50;
425 else if(nbin<1000)
bin =100;
426 else if(nbin<2000)
bin =100;
427 else if(nbin<5000)
bin =500;
430 }
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475}
*******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
void fill_histogram(double x, double y, double r, int nbin, TH2D *hough)