CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
houghmap2D.cxx
Go to the documentation of this file.
1int houghmap2D()
2{
3 const double M_PI = 3.1415926;
4 //////////////////////////////////////////////////////////
5 // This file has been automatically generated
6 // (Wed Jan 10 14:20:10 2018 by ROOT version5.34/09)
7 // from TTree hit/hit
8 // found on file: hough_hist_test.root
9 //////////////////////////////////////////////////////////
10
11
12 //Reset ROOT and connect tree file
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 //Declaration of leaves types
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 // Set branch addresses.
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 // This is the loop skeleton
61 // To read only selected branches, Insert statements like:
62 // hit->SetBranchStatus("*",0); // disable all branches
63 // TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
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;
75 double x_max = M_PI;
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 //for (Long64_t i=0; i<nentries;i++)
89 //for (Long64_t i=0; i<1000;i++)
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 //cout<<i<<" "<<hit_nhit<<endl;
97 for(int j=0;j<hit_nhit;j++){
98 if(hit_flag[j] !=0)continue;
99 nhit++;
100 //if(hit_layer[j] >35)continue;
101 //cout<<hit_layer[j]<<endl;
102 //fill_histogram(hit_x[j],hit_y[j],hit_driftdist[j],x_bin*2,hough);
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 //if(hit_layer[j]==11){
111 //c_hit = new TCanvas(name,name,700,500);
112 //houghhit[ihit]->Draw();
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 //c_map->SaveAs(name);
128 ssname.str("");
129 ssname <<"hough_map.pdf";
130 sname = ssname.str();
131 name = sname.c_str();
132 //c_map->SaveAs(name);
133 ssname.str("");
134 ssname <<"hough_map.eps";
135 sname = ssname.str();
136 name = sname.c_str();
137 //c_map->SaveAs(name);
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 //cout<<theta*180/3.1415<<" "<<rho<<endl;
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);
158 if(alpha<0.)alpha+=2.*M_PI;
159 if(alpha>2.*M_PI)alpha-=2.*M_PI;
160 if(alpha>alpha_max)alpha_max=alpha;
161 //cout<<alpha<<" "<<alpha_max<<endl;
162 }
163 //TEllipse *e = new TEllipse(Xc,Yc,Rc,Rc,180-alpha_max*180/M_PI,180);
164 TEllipse *e = new TEllipse(Xc,Yc,Rc,Rc);
165 //e->Draw("Asame");
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 //TGraph *gr_xy_DriftCircle= new TGraph();
174 //TGraph *gr_uv_DriftCircle= new TGraph();
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 //cout<<sqrt(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j])<<endl;
181
182 //double u= 2*hit_x[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
183 //double v= 2*hit_y[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
184 //double w= 2*hit_driftdist[j]/(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]-hit_driftdist[j]*hit_driftdist[j]);
185 //gr_uv->SetPoint(point,u,v);
186 //TEllipse *euv = new TEllipse(u,v,w,w);
187 //euv->Draw("Csame");
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 double min = x1<y1?x1:y1;
205 double max = x2>y2?x2:y2;
206 x_min = min-0.1*(max-min);
207 y_min = min-0.1*(max-min);
208 x_max = max+0.1*(max-min);
209 y_max = max+0.1*(max-min);
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 //gr_uv->Draw("Psame");
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 //cout<<layersOnPeak[ihit]<<" "<<houghhit[ihit]->GetBinContent(binx,biny)<<endl;
241 //delete houghhit[ihit];
242 }
243 ssname.str("");
244 ssname <<"2d_track.eps";
245 sname = ssname.str();
246 name = sname.c_str();
247 //c_xy->SaveAs(name);
248 ssname.str("");
249 ssname <<"2d_track.pdf";
250 sname = ssname.str();
251 name = sname.c_str();
252 //c_xy->SaveAs(name);
253 ssname.str("");
254 ssname <<"2d_track.png";
255 sname = ssname.str();
256 name = sname.c_str();
257 //c_xy->SaveAs(name);
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 //TEllipse *euv = new TEllipse(u,v,w,w);
267 //euv->Draw("Csame");
268 //euv->SetFillStyle(0);
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 //gr_uv->SetPoint(point,u,v);
283 TEllipse *euv = new TEllipse(u,v,w,w);
284 euv->Draw("Csame");
285 euv->SetFillStyle(0);
286 //point++;
287 }
288 //delete hough;
289 }
290/*
291 ssname.str("");
292 ssname <<nbin<<"bin each layer deltaD";
293 sname = ssname.str();
294 name = sname.c_str();
295 TCanvas *c2 = new TCanvas(name,name,2500,2500 );
296 c2->Divide(5,5);
297 ssname.str("");
298 ssname <<nbin<<"residual_each_layer.eps";
299 sname = ssname.str();
300 name = sname.c_str();
301 c2->SaveAs(name);
302 ssname.str("");
303 ssname <<nbin<<"residual_each_layer.pdf";
304 sname = ssname.str();
305 name = sname.c_str();
306 c2->SaveAs(name);
307 ssname.str("");
308 ssname <<nbin<<"residual_each_layer.png";
309 sname = ssname.str();
310 name = sname.c_str();
311 c2->SaveAs(name);
312
313/*
314 ssname.str("");
315 ssname <<nbin<<"bin all layers deltaD";
316 sname = ssname.str();
317 name = sname.c_str();
318 TCanvas *c3 = new TCanvas(name,name,700,500 );
319 int yMax = h_deltaD_all->GetMaximum() > h_deltaD_MC_all->GetMaximum()? h_deltaD_all->GetMaximum():h_deltaD_MC_all->GetMaximum();
320 h_deltaD_all->Draw("same");
321 h_deltaD_all->SetLineColor(2);
322 h_deltaD_MC_all->Draw("same");
323 h_deltaD_MC_all->SetLineColor(4);
324 x_axis = h_deltaD_all->GetXaxis();
325 y_axis = h_deltaD_all->GetYaxis();
326 x_axis->SetTitle("residual (cm)");
327 x_axis->CenterTitle();
328 x_axis->SetTitleSize(0.05);
329 x_axis->SetLabelSize(0.04);
330 x_axis->SetTitleOffset(1.05);
331 y_axis->SetTitle("Entries");
332 y_axis->CenterTitle();
333 y_axis->SetTitleSize(0.05);
334 y_axis->SetLabelSize(0.04);
335 y_axis->SetTitleOffset(1.25);
336 y_axis->SetRangeUser(0,yMax*1.1);
337 TLegend* leg3= new TLegend(0.12 ,0.80,0.33,0.88);
338 leg3->AddEntry(h_deltaD_all,"measurement","l");
339 leg3->AddEntry(h_deltaD_MC_all,"MC truth","l");
340 leg3->SetFillColor(kWhite);
341 leg3->SetLineColor(kWhite);
342 leg3->Draw();
343
344 ssname.str("");
345 ssname <<nbin<<"residual_all_layer.eps";
346 sname = ssname.str();
347 name = sname.c_str();
348 c3->SaveAs(name);
349 ssname.str("");
350 ssname <<nbin<<"residual_all_layer.pdf";
351 sname = ssname.str();
352 name = sname.c_str();
353 c3->SaveAs(name);
354 ssname.str("");
355 ssname <<nbin<<"residual_all_layer.png";
356 sname = ssname.str();
357 name = sname.c_str();
358 c3->SaveAs(name);
359
360
361 ssname.str("");
362 ssname <<nbin<<"bin deltaD";
363 sname = ssname.str();
364 name = sname.c_str();
365 TCanvas *c4 = new TCanvas(name,name,1500,500 );
366 c4->Divide(3,1);
367 c4->cd(1);
368 h_deltaD[22]->Draw("same");
369 h_deltaD[22]->Fit("gaus","Q");
370 c4->cd(2);
371 h_deltaD[23]->Draw("same");
372 h_deltaD[23]->Fit("gaus","Q");
373 c4->cd(3);
374 h_deltaD[24]->Draw("same");
375 h_deltaD[24]->Fit("gaus","Q");
376
377 ssname.str("");
378 ssname <<nbin<<"residual.eps";
379 sname = ssname.str();
380 name = sname.c_str();
381 c4->SaveAs(name);
382 ssname.str("");
383 ssname <<nbin<<"residual.pdf";
384 sname = ssname.str();
385 name = sname.c_str();
386 c4->SaveAs(name);
387 ssname.str("");
388 ssname <<nbin<<"residual.png";
389 sname = ssname.str();
390 name = sname.c_str();
391 c4->SaveAs(name);
392
393
394 double sigma_1_3 = h_deltaD[22]->GetFunction("gaus")->GetParameter(2);
395 double sigma_err_1_3 = h_deltaD[22]->GetFunction("gaus")->GetParError(2);
396 double sigma_9_20 = h_deltaD[23]->GetFunction("gaus")->GetParameter(2);
397 double sigma_err_1_9_20 = h_deltaD[23]->GetFunction("gaus")->GetParError(2);
398 double sigma_37_43 = h_deltaD[24]->GetFunction("gaus")->GetParameter(2);
399 double sigma_err_1_37_43 = h_deltaD[24]->GetFunction("gaus")->GetParError(2);
400 double RMS_1_3 = h_deltaD[22]->GetRMS();
401 double RMS_err_1_3 = h_deltaD[22]->GetRMSError();
402 double RMS_9_20 = h_deltaD[23]->GetRMS();
403 double RMS_err_1_9_20 = h_deltaD[23]->GetRMSError();
404 double RMS_37_43 = h_deltaD[24]->GetRMS();
405 double RMS_err_1_37_43 = h_deltaD[24]->GetRMSError();
406 cout<<nbin<<" "<<sigma_1_3<<" "<<sigma_err_1_3<<" "<<sigma_9_20<<" "<<sigma_err_1_9_20<<" "<<sigma_37_43<<" "<<sigma_err_1_37_43<<endl;
407 gr_sigma_1_3->SetPoint(ibin,nbin,sigma_1_3);
408 gr_sigma_1_3->SetPointError(ibin,0,sigma_err_1_3);
409 gr_sigma_9_20->SetPoint(ibin,nbin,sigma_9_20);
410 gr_sigma_9_20->SetPointError(ibin,0,sigma_err_1_9_20);
411 gr_sigma_37_43->SetPoint(ibin,nbin,sigma_37_43);
412 gr_sigma_37_43->SetPointError(ibin,0,sigma_err_1_37_43);
413
414 gr_RMS_1_3->SetPoint(ibin,nbin,RMS_1_3);
415 gr_RMS_1_3->SetPointError(ibin,0,RMS_err_1_3);
416 gr_RMS_9_20->SetPoint(ibin,nbin,RMS_9_20);
417 gr_RMS_9_20->SetPointError(ibin,0,RMS_err_1_9_20);
418 gr_RMS_37_43->SetPoint(ibin,nbin,RMS_37_43);
419 gr_RMS_37_43->SetPointError(ibin,0,RMS_err_1_37_43);
420//*/
421
422 ibin++;
423 int bin(0);
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;
428 else bin = 1000;
429 nbin += bin;
430 }
431/*
432 TCanvas *c5 = new TCanvas("sigma_deltaD","sigma_deltaD",700,500 );
433 gr_sigma_1_3->SetMinimum(0.0);
434 gr_sigma_1_3->SetMaximum(1.20);
435 gr_sigma_1_3->Draw("APsame");
436 gr_sigma_1_3->SetMarkerStyle(20);
437 gr_sigma_1_3->SetMarkerSize(1.0);
438 gr_sigma_1_3->SetMarkerColor(2);
439 gr_sigma_9_20->Draw(" Psame");
440 gr_sigma_9_20->SetMarkerStyle(20);
441 gr_sigma_9_20->SetMarkerSize(1.0);
442 gr_sigma_9_20->SetMarkerColor(3);
443 gr_sigma_37_43->Draw(" Psame");
444 gr_sigma_37_43->SetMarkerStyle(20);
445 gr_sigma_37_43->SetMarkerSize(1.0);
446 gr_sigma_37_43->SetMarkerColor(4);
447
448 gr_RMS_1_3->Draw("Psame");
449 gr_RMS_1_3->SetMarkerStyle(25);
450 gr_RMS_1_3->SetMarkerSize(1.0);
451 gr_RMS_1_3->SetMarkerColor(2);
452 gr_RMS_9_20->Draw(" Psame");
453 gr_RMS_9_20->SetMarkerStyle(25);
454 gr_RMS_9_20->SetMarkerSize(1.0);
455 gr_RMS_9_20->SetMarkerColor(3);
456 gr_RMS_37_43->Draw(" Psame");
457 gr_RMS_37_43->SetMarkerStyle(25);
458 gr_RMS_37_43->SetMarkerSize(1.0);
459 gr_RMS_37_43->SetMarkerColor(4);
460
461 TLegend* leg5= new TLegend(0.55 ,0.55,0.88,0.88);
462 leg5->AddEntry(gr_sigma_1_3," 1~3 layer sigma","p");
463 leg5->AddEntry(gr_RMS_1_3," 1~3 layer RMS","p");
464 leg5->AddEntry(gr_sigma_9_20,"21~36 layer sigma","p");
465 leg5->AddEntry(gr_RMS_9_20,"21~36 layer RMS","p");
466 leg5->AddEntry(gr_sigma_37_43,"37~43 layer sigma","p");
467 leg5->AddEntry(gr_RMS_37_43,"37~43 layer RMS","p");
468 leg5->SetFillColor(kWhite);
469 leg5->SetLineColor(kWhite);
470 leg5->Draw();
471 c5->SaveAs("sigma_deltaD.eps");
472 c5->SaveAs("sigma_deltaD.png");
473 c5->SaveAs("sigma_deltaD.pdf");
474//*/
475}
476
477double intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
478{
479 const double M_PI = 3.1415926;
480 double phi;
481 double phi_center = atan2(y_center,x_center);
482 double r_center = sqrt(x_center*x_center+y_center*y_center);
483 if(charge==0)return 9999;
484 while(phi_center<0) phi_center += 2*M_PI;
485 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
486 if(r_center<r_cylinder/2) return -9999;
487 double dphi = acos( r_cylinder/(2*r_center) );
488 if(charge<0) phi = phi_center + dphi;
489 else phi = phi_center - dphi;
490 while(phi<0) phi += 2*M_PI;
491 while(phi>2*M_PI) phi -= 2*M_PI;
492 return phi;
493}
494
495void fill_histogram(double x, double y, double r, int nbin, TH2D* hough)
496{
497 const double M_PI = 3.1415926;
498 double U = 2*x/(x*x+y*y-r*r);
499 double V = 2*y/(x*x+y*y-r*r);
500 double R = 2*r/(x*x+y*y-r*r);
501 double delta_phi = 2.*M_PI/nbin;
502 double phi = -delta_phi/2;
503 for(int i =0; i<nbin; i++){
504 phi += delta_phi;
505 double u = U + R*cos(phi);
506 double v = V + R*sin(phi);
507 double normal = (v-V)/(u-U);
508 double k = -1./normal;
509 double b = v - k*u;
510 double x_cross = -b/(k+1/k);
511 double y_cross = b/(1+k*k);
512
513 double rho=sqrt(x_cross*x_cross+y_cross*y_cross);
514 double theta=atan2(y_cross,x_cross);
515 double slantOfLine = V*cos(theta)-U*sin(theta);
516 int chargeOfHitOnCir = (slantOfLine/fabs(slantOfLine)) * (rho/fabs(rho));
517 if( chargeOfHitOnCir == -1) continue;
518 //if( fabs(slantOfLine)<0.03 ) continue;
519 if(theta<0) {
520 theta=theta+M_PI;
521 rho=-rho;
522 }
523 if( normal ==0 && x>0) {
524 rho = fabs(x);
525 theta = 0;
526 }
527 if( normal ==0 && x<0) {
528 rho = -fabs(x);
529 theta = M_PI;
530 }
531 if(fabs(rho)>0.1)continue;
532 hough->Fill(theta,rho);
533 }
534}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
nentries
Int_t nentries
Definition: DataBase/tau_mode.c:8
alpha
const double alpha
Definition: FastVertexFit.cxx:4
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
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
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
M_PI
#define M_PI
Definition: TConstant.h:4
fill_histogram
void fill_histogram(double x, double y, double r, int nbin, TH2D *hough)
Definition: houghmap2D.cxx:495
houghmap2D
int houghmap2D()
Definition: houghmap2D.cxx:1
intersect_cylinder
double intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
Definition: houghmap2D.cxx:477