CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
bin_deltaD.cxx
Go to the documentation of this file.
1int bin_deltaD()
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[1000];
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 TGraphErrors *gr_sigma_1_3 = new TGraphErrors();
71 TGraphErrors *gr_sigma_9_20 = new TGraphErrors();
72 TGraphErrors *gr_sigma_37_43 = new TGraphErrors();
73 TGraphErrors *gr_RMS_1_3 = new TGraphErrors();
74 TGraphErrors *gr_RMS_9_20 = new TGraphErrors();
75 TGraphErrors *gr_RMS_37_43 = new TGraphErrors();
76 TCanvas *c = new TCanvas("sigma_layer","sigma_layer",1000,1000 );
77 TLegend* leg = new TLegend(0.12 ,0.80,0.33,0.88);
78 TGraphErrors *gr_sigma[100];
79 int ibin = 0;
80 int nbin = 100 ;
81 while(nbin <= 3000){
82 int x_bin = nbin;
83 int y_bin = nbin;
84 double x_max = M_PI;
85 double y_max = 0.1;
86
87 ssname.str("");
88 ssname <<nbin<<"bin, all layers";
89 sname = ssname.str();
90 name = sname.c_str();
91 TH1D *h_deltaD_all = new TH1D(name,name,100,-1,1);
92 ssname.str("");
93 ssname <<nbin<<"bin, MC all layers";
94 sname = ssname.str();
95 name = sname.c_str();
96 TH1D *h_deltaD_MC_all = new TH1D(name,name,100,-1,1);
97 TH1D *h_deltaD[25];
98 TH1D *h_deltaD_MC[25];
99 for(int l=1;l<=25;l++){
100 double cut;
101 ssname.str("");
102 ssname <<nbin<<"bin, ";
103 if(l<=3){ssname <<"layer "<<l;cut=1+50/nbin;}
104 else if(l<=15){ssname <<"layer "<<l+5;cut=0.5+(1000/nbin)*0.1;}
105 else if(l<=22){ssname <<"layer "<<l+21;cut=1+500/nbin;}
106 else if(l==23){ssname <<"1~3 layer";cut=1+50/nbin;}
107 else if(l==24){ssname <<"9~20 layer";cut=0.5+(1000/nbin)*0.1;}
108 else if(l==25){ssname <<"37~43 layer";cut=1+500/nbin;}
109 sname = ssname.str();
110 name = sname.c_str();
111 h_deltaD[l-1]= new TH1D(name,name,100,-cut,cut);
112
113 ssname.str("");
114 ssname <<nbin<<"bin, ";
115 if(l<=3){ssname <<"MC layer "<<l;cut=1+50/nbin;}
116 //else if(l<=15){ssname <<"MC layer "<<l+5;cut=1+100/nbin;}
117 else if(l<=15){ssname <<"MC layer "<<l+5;cut=0.5+(1000/nbin)*0.1;}
118 else if(l<=22){ssname <<"MC layer "<<l+21;cut=1+500/nbin;}
119 else if(l==23){ssname <<"MC 1~3 layer ";cut=1+50/nbin;}
120 else if(l==24){ssname <<"MC 9~20 layer ";cut=0.5+(1000/nbin)*0.1;}
121 else if(l==25){ssname <<"MC 37~43 layer ";cut=1+500/nbin;}
122 sname = ssname.str();
123 name = sname.c_str();
124 h_deltaD_MC[l-1]= new TH1D(name,name,100,-cut,cut);
125 }
126
127 int event =11 ;
128 Long64_t nentries = hit->GetEntries();
129 Long64_t nbytes = 0;
130 for (Long64_t i=0; i<nentries;i++)
131 //for (Long64_t i=0; i<1000;i++)
132 //for (Long64_t i=event; i<event+1;i++)
133 {
134 nbytes += hit->GetEntry(i);
135 TH2D* hough = new TH2D("hough","hough",x_bin,0,x_max,y_bin,-y_max,y_max);
136 TH2D *houghhit[100];
137 int layersOnPeak[100];
138 int nhit(0),ihit(0);
139 for(int j=0;j<hit_nhit;j++){
140 if(hit_flag[j] !=0)continue;
141 nhit++;
142 //if(hit_layer[j] >35)continue;
143 //cout<<hit_layer[j]<<endl;
144 //fill_histogram(hit_x[j],hit_y[j],hit_driftdist[j],x_bin*2,hough);
145 ssname.str("");
146 ssname <<"layer:"<<hit_layer[j]<<" ,wire"<<hit_wire[j];
147 sname = ssname.str();
148 name = sname.c_str();
149 layersOnPeak[ihit] = hit_layer[j];
150 houghhit[ihit] = new TH2D(name,name,x_bin,0,x_max,y_bin,-y_max,y_max);
151 fill_histogram(hit_x[j],hit_y[j],hit_driftdist[j],x_bin*2,houghhit[ihit]);
152 hough->Add(houghhit[ihit]);
153 ihit++;
154 }
155 /*
156 TCanvas *c1 = new TCanvas("houghspace","houghspace",1000,500 );
157 c1->Divide(2,1);
158 c1->cd(1);
159 hough->Draw();
160 c1->cd(2);
161 hough->Draw("LEGO2Z");
162 //*/
163 int binx,biny,binz;
164 hough->GetMaximumBin(binx,biny,binz);
165 double theta = hough->GetXaxis()->GetBinCenter(binx);
166 double rho = hough->GetYaxis()->GetBinCenter(biny);
167 for(int ihit=0;ihit<nhit;ihit++){
168 //cout<<layersOnPeak[ihit]<<" "<<houghhit[ihit]->GetBinContent(binx,biny)<<endl;
169 delete houghhit[ihit];
170 }
171 double Rc = 1./fabs(rho);
172 double Xc = cos(theta)/rho;
173 double Yc = sin(theta)/rho;
174 double deltaD(0);
175 for(int j=0;j<hit_nhit;j++){
176 if(hit_flag[j] !=0)continue;
177 if(hit_layer[j] <3){
178 double r_cgem = sqrt(hit_x[j]*hit_x[j]+hit_y[j]*hit_y[j]);
179 double phi_cross = intersect_cylinder(-1,Xc,Yc,r_cgem);
180 double phih_cluster = atan2(hit_y[j],hit_x[j]);
181 double dphi = phih_cluster - phi_cross;
182 if(dphi>M_PI)dphi-=2*M_PI;
183 if(dphi<-1*M_PI)dphi+=2*M_PI;
184 deltaD = r_cgem*dphi;
185 }else{
186 double l1l2 = sqrt((hit_x[j]-Xc)*(hit_x[j]-Xc)+(hit_y[j]-Yc)*(hit_y[j]-Yc));
187 if(l1l2>Rc)deltaD = l1l2-Rc-hit_driftdist[j];
188 else deltaD = l1l2-Rc+hit_driftdist[j];
189 }
190 //cout<<hit_layer[j]<<" "<<deltaD<<endl;
191 h_deltaD_all->Fill(deltaD);
192 if(hit_layer[j] <3) {h_deltaD[hit_layer[j]]->Fill(deltaD);h_deltaD[22]->Fill(deltaD);}
193 else if(hit_layer[j] <20) {h_deltaD[hit_layer[j]-5]->Fill(deltaD);h_deltaD[23]->Fill(deltaD);}
194 else {h_deltaD[hit_layer[j]-21]->Fill(deltaD);h_deltaD[24]->Fill(deltaD);}
195 }
196
197 for(int j=0;j<hit_nhit;j++){
198 if(hit_flag[j] !=0)continue;
199 if(hit_layer[j] <3){
200 double r_cgem = sqrt(hit_truth_x[j]*hit_truth_x[j]+hit_truth_y[j]*hit_truth_y[j]);
201 double phi_cross = intersect_cylinder(-1,Xc,Yc,r_cgem);
202 double phih_cluster = atan2(hit_truth_y[j],hit_truth_x[j]);
203 double dphi = phih_cluster - phi_cross;
204 if(dphi>M_PI)dphi-=2*M_PI;
205 if(dphi<-1*M_PI)dphi+=2*M_PI;
206 deltaD = r_cgem*dphi;
207 }else{
208 double l1l2 = sqrt((hit_truth_x[j]-Xc)*(hit_truth_x[j]-Xc)+(hit_truth_y[j]-Yc)*(hit_truth_y[j]-Yc));
209 deltaD = l1l2-Rc;
210 }
211 //cout<<hit_layer[j]<<" "<<deltaD<<endl;
212 h_deltaD_MC_all->Fill(deltaD);
213 if(hit_layer[j] <3) {h_deltaD_MC[hit_layer[j]]->Fill(deltaD);h_deltaD_MC[22]->Fill(deltaD);}
214 else if(hit_layer[j] <20) {h_deltaD_MC[hit_layer[j]-5]->Fill(deltaD);h_deltaD_MC[23]->Fill(deltaD);}
215 else {h_deltaD_MC[hit_layer[j]-21]->Fill(deltaD);h_deltaD_MC[24]->Fill(deltaD);}
216 }
217 delete hough;
218 }
219 //TCanvas *c2 = new TCanvas("deltaD","deltaD",700,500 );
220 //h_deltaD[0]->Draw();
221 //h_deltaD[0]->SetLineColor(2);
222 //h_deltaD_MC[0]->Draw("same");
223 //h_deltaD_MC[0]->SetLineColor(4);
224 ssname.str("");
225 ssname <<nbin<<"bin each layer deltaD";
226 sname = ssname.str();
227 name = sname.c_str();
228 TCanvas *c2 = new TCanvas(name,name,2500,2500 );
229 c2->Divide(5,5);
230 gr_sigma[ibin] = new TGraphErrors();
231 double sigma[25];
232 double error[25];
233 int layer[25]={1,2,3,8,9,10,11,12,13,14,15,16,17,18,19,20,37,38,39,40,41,42,43,0,0,0};
234 cout<<"-------------------------------------------------- "<<nbin<<" --------------------------------------------------"<<endl;
235 for(int l=0;l<25;l++){
236 c2->cd(l+1);
237 int yMax = h_deltaD[l]->GetMaximum() > h_deltaD_MC[l]->GetMaximum()? h_deltaD[l]->GetMaximum():h_deltaD_MC[l]->GetMaximum();
238 h_deltaD[l]->Draw("same");
239 h_deltaD[l]->SetLineColor(4);
240 h_deltaD[l]->Fit("gaus","Q");
241 h_deltaD[l]->GetFunction("gaus")->GetParameter(2);
242 sigma[l] = h_deltaD[l]->GetFunction("gaus")->GetParameter(2);
243 error[l] = h_deltaD[l]->GetFunction("gaus")->GetParError(2);
244 cout<<sigma[l]<<", ";
245 x_axis = h_deltaD[l]->GetXaxis();
246 y_axis = h_deltaD[l]->GetYaxis();
247 x_axis->SetTitle("residual (cm)");
248 x_axis->CenterTitle();
249 //x_axis->SetTitleSize(0.05);
250 //x_axis->SetLabelSize(0.04);
251 //x_axis->SetTitleOffset(1.05);
252 //x_axis->SetLimits(0,80 );
253 y_axis->SetTitle("Entries");
254 y_axis->CenterTitle();
255 //y_axis->SetTitleSize(0.05);
256 //y_axis->SetLabelSize(0.04);
257 //y_axis->SetTitleOffset(1.25);
258 y_axis->SetRangeUser(0,yMax*1.1);
259 h_deltaD_MC[l]->Draw("same");
260 h_deltaD_MC[l]->SetLineColor(2);
261 TLegend* leg2= new TLegend(0.12 ,0.80,0.33,0.88);
262 leg2->AddEntry(h_deltaD[l],"measurement","l");
263 leg2->AddEntry(h_deltaD_MC[l],"MC truth","l");
264 leg2->SetFillColor(kWhite);
265 leg2->SetLineColor(kWhite);
266 leg2->Draw();
267 if(l<22){
268 gr_sigma[ibin]->SetPoint(l,layer[l],sigma[l]);
269 gr_sigma[ibin]->SetPointError(l,0,error[l]);
270 }
271
272
273 }
274 c->cd();
275 if(ibin==0)gr_sigma[ibin]->Draw("APLsame");
276 else gr_sigma[ibin]->Draw("PLsame");
277 gr_sigma[ibin]->SetMarkerStyle(20);
278 gr_sigma[ibin]->SetMarkerSize(1.0);
279 gr_sigma[ibin]->SetMarkerColor(ibin+1);
280 gr_sigma[ibin]->SetLineColor(ibin+1);
281 ssname.str("");
282 ssname <<nbin<<"bin ";
283 sname = ssname.str();
284 name = sname.c_str();
285 leg->AddEntry(gr_sigma[ibin],name,"p");
286
287 ssname.str("");
288 ssname <<nbin<<"residual_each_layer.eps";
289 sname = ssname.str();
290 name = sname.c_str();
291 c2->SaveAs(name);
292 ssname.str("");
293 ssname <<nbin<<"residual_each_layer.pdf";
294 sname = ssname.str();
295 name = sname.c_str();
296 c2->SaveAs(name);
297 ssname.str("");
298 ssname <<nbin<<"residual_each_layer.png";
299 sname = ssname.str();
300 name = sname.c_str();
301 c2->SaveAs(name);
302
303 ssname.str("");
304 ssname <<nbin<<"bin all layers deltaD";
305 sname = ssname.str();
306 name = sname.c_str();
307 TCanvas *c3 = new TCanvas(name,name,700,500 );
308 int yMax = h_deltaD_all->GetMaximum() > h_deltaD_MC_all->GetMaximum()? h_deltaD_all->GetMaximum():h_deltaD_MC_all->GetMaximum();
309 h_deltaD_all->Draw("same");
310 h_deltaD_all->SetLineColor(2);
311 h_deltaD_MC_all->Draw("same");
312 h_deltaD_MC_all->SetLineColor(4);
313 x_axis = h_deltaD_all->GetXaxis();
314 y_axis = h_deltaD_all->GetYaxis();
315 x_axis->SetTitle("residual (cm)");
316 x_axis->CenterTitle();
317 x_axis->SetTitleSize(0.05);
318 x_axis->SetLabelSize(0.04);
319 x_axis->SetTitleOffset(1.05);
320 y_axis->SetTitle("Entries");
321 y_axis->CenterTitle();
322 y_axis->SetTitleSize(0.05);
323 y_axis->SetLabelSize(0.04);
324 y_axis->SetTitleOffset(1.25);
325 y_axis->SetRangeUser(0,yMax*1.1);
326 TLegend* leg3= new TLegend(0.12 ,0.80,0.33,0.88);
327 leg3->AddEntry(h_deltaD_all,"measurement","l");
328 leg3->AddEntry(h_deltaD_MC_all,"MC truth","l");
329 leg3->SetFillColor(kWhite);
330 leg3->SetLineColor(kWhite);
331 leg3->Draw();
332
333 ssname.str("");
334 ssname <<nbin<<"residual_all_layer.eps";
335 sname = ssname.str();
336 name = sname.c_str();
337 c3->SaveAs(name);
338 ssname.str("");
339 ssname <<nbin<<"residual_all_layer.pdf";
340 sname = ssname.str();
341 name = sname.c_str();
342 c3->SaveAs(name);
343 ssname.str("");
344 ssname <<nbin<<"residual_all_layer.png";
345 sname = ssname.str();
346 name = sname.c_str();
347 c3->SaveAs(name);
348
349 ssname.str("");
350 ssname <<nbin<<"bin deltaD";
351 sname = ssname.str();
352 name = sname.c_str();
353 TCanvas *c4 = new TCanvas(name,name,1500,500 );
354 c4->Divide(3,1);
355 c4->cd(1);
356 h_deltaD[22]->Draw("same");
357 h_deltaD[22]->Fit("gaus","Q");
358 c4->cd(2);
359 h_deltaD[23]->Draw("same");
360 h_deltaD[23]->Fit("gaus","Q");
361 c4->cd(3);
362 h_deltaD[24]->Draw("same");
363 h_deltaD[24]->Fit("gaus","Q");
364
365 ssname.str("");
366 ssname <<nbin<<"residual.eps";
367 sname = ssname.str();
368 name = sname.c_str();
369 c4->SaveAs(name);
370 ssname.str("");
371 ssname <<nbin<<"residual.pdf";
372 sname = ssname.str();
373 name = sname.c_str();
374 c4->SaveAs(name);
375 ssname.str("");
376 ssname <<nbin<<"residual.png";
377 sname = ssname.str();
378 name = sname.c_str();
379 c4->SaveAs(name);
380
381 double sigma_1_3 = h_deltaD[22]->GetFunction("gaus")->GetParameter(2);
382 double sigma_err_1_3 = h_deltaD[22]->GetFunction("gaus")->GetParError(2);
383 double sigma_9_20 = h_deltaD[23]->GetFunction("gaus")->GetParameter(2);
384 double sigma_err_1_9_20 = h_deltaD[23]->GetFunction("gaus")->GetParError(2);
385 double sigma_37_43 = h_deltaD[24]->GetFunction("gaus")->GetParameter(2);
386 double sigma_err_1_37_43 = h_deltaD[24]->GetFunction("gaus")->GetParError(2);
387 double RMS_1_3 = h_deltaD[22]->GetRMS();
388 double RMS_err_1_3 = h_deltaD[22]->GetRMSError();
389 double RMS_9_20 = h_deltaD[23]->GetRMS();
390 double RMS_err_1_9_20 = h_deltaD[23]->GetRMSError();
391 double RMS_37_43 = h_deltaD[24]->GetRMS();
392 double RMS_err_1_37_43 = h_deltaD[24]->GetRMSError();
393 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;
394 gr_sigma_1_3->SetPoint(ibin,nbin,sigma_1_3);
395 gr_sigma_1_3->SetPointError(ibin,0,sigma_err_1_3);
396 gr_sigma_9_20->SetPoint(ibin,nbin,sigma_9_20);
397 gr_sigma_9_20->SetPointError(ibin,0,sigma_err_1_9_20);
398 gr_sigma_37_43->SetPoint(ibin,nbin,sigma_37_43);
399 gr_sigma_37_43->SetPointError(ibin,0,sigma_err_1_37_43);
400
401 gr_RMS_1_3->SetPoint(ibin,nbin,RMS_1_3);
402 gr_RMS_1_3->SetPointError(ibin,0,RMS_err_1_3);
403 gr_RMS_9_20->SetPoint(ibin,nbin,RMS_9_20);
404 gr_RMS_9_20->SetPointError(ibin,0,RMS_err_1_9_20);
405 gr_RMS_37_43->SetPoint(ibin,nbin,RMS_37_43);
406 gr_RMS_37_43->SetPointError(ibin,0,RMS_err_1_37_43);
407
408 for(int l=0;l<25;l++){
409 cout<<sigma[l]<<endl;
410 delete h_deltaD[l];
411 delete h_deltaD_MC[l];
412 }
413 delete h_deltaD_all;
414 delete h_deltaD_MC_all;
415 cout<<endl;
416
417 int bin;
418 if(nbin<500) bin = 50;
419 else if(nbin<1000) bin =100;
420 else if(nbin<2000)bin =200;
421 else if(nbin<5000)bin =500;
422 else bin = 1000;
423 nbin += bin;
424 ibin++;
425 }
426 TCanvas *c5 = new TCanvas("sigma_deltaD","sigma_deltaD",700,500 );
427 gr_sigma_1_3->SetMinimum(0.0);
428 gr_sigma_1_3->SetMaximum(1.20);
429 gr_sigma_1_3->Draw("APsame");
430 gr_sigma_1_3->SetMarkerStyle(20);
431 gr_sigma_1_3->SetMarkerSize(1.0);
432 gr_sigma_1_3->SetMarkerColor(2);
433 gr_sigma_9_20->Draw(" Psame");
434 gr_sigma_9_20->SetMarkerStyle(20);
435 gr_sigma_9_20->SetMarkerSize(1.0);
436 gr_sigma_9_20->SetMarkerColor(3);
437 gr_sigma_37_43->Draw(" Psame");
438 gr_sigma_37_43->SetMarkerStyle(20);
439 gr_sigma_37_43->SetMarkerSize(1.0);
440 gr_sigma_37_43->SetMarkerColor(4);
441
442 gr_RMS_1_3->Draw("Psame");
443 gr_RMS_1_3->SetMarkerStyle(25);
444 gr_RMS_1_3->SetMarkerSize(1.0);
445 gr_RMS_1_3->SetMarkerColor(2);
446 gr_RMS_9_20->Draw(" Psame");
447 gr_RMS_9_20->SetMarkerStyle(25);
448 gr_RMS_9_20->SetMarkerSize(1.0);
449 gr_RMS_9_20->SetMarkerColor(3);
450 gr_RMS_37_43->Draw(" Psame");
451 gr_RMS_37_43->SetMarkerStyle(25);
452 gr_RMS_37_43->SetMarkerSize(1.0);
453 gr_RMS_37_43->SetMarkerColor(4);
454
455 TLegend* leg5= new TLegend(0.55 ,0.55,0.88,0.88);
456 leg5->AddEntry(gr_sigma_1_3," 1~3 layer sigma","p");
457 leg5->AddEntry(gr_RMS_1_3," 1~3 layer RMS","p");
458 leg5->AddEntry(gr_sigma_9_20," 9~20 layer sigma","p");
459 leg5->AddEntry(gr_RMS_9_20," 9~20 layer RMS","p");
460 leg5->AddEntry(gr_sigma_37_43,"37~43 layer sigma","p");
461 leg5->AddEntry(gr_RMS_37_43,"37~43 layer RMS","p");
462 leg5->SetFillColor(kWhite);
463 leg5->SetLineColor(kWhite);
464 leg5->Draw();
465 c5->SaveAs("sigma_deltaD.eps");
466 c5->SaveAs("sigma_deltaD.png");
467 c5->SaveAs("sigma_deltaD.pdf");
468
469 c->cd();
470 leg->Draw();
471 ssname.str("");
472 ssname <<"sigma_layer.png";
473 sname = ssname.str();
474 name = sname.c_str();
475 c->SaveAs(name);
476 ssname.str("");
477 ssname <<"sigma_layer.pdf";
478 sname = ssname.str();
479 name = sname.c_str();
480 c->SaveAs(name);
481 ssname.str("");
482 ssname <<"sigma_layer.eps";
483 sname = ssname.str();
484 name = sname.c_str();
485 c->SaveAs(name);
486}
487
488double intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
489{
490 const double M_PI = 3.1415926;
491 double phi;
492 double phi_center = atan2(y_center,x_center);
493 double r_center = sqrt(x_center*x_center+y_center*y_center);
494 if(charge==0)return 9999;
495 while(phi_center<0) phi_center += 2*M_PI;
496 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
497 if(r_center<r_cylinder/2) return -9999;
498 double dphi = acos( r_cylinder/(2*r_center) );
499 if(charge<0) phi = phi_center + dphi;
500 else phi = phi_center - dphi;
501 while(phi<0) phi += 2*M_PI;
502 while(phi>2*M_PI) phi -= 2*M_PI;
503 return phi;
504}
505
506void fill_histogram(double x, double y, double r, int nbin, TH2D* hough)
507{
508 const double M_PI = 3.1415926;
509 double U = 2*x/(x*x+y*y-r*r);
510 double V = 2*y/(x*x+y*y-r*r);
511 double R = 2*r/(x*x+y*y-r*r);
512 double delta_phi = 2.*M_PI/nbin;
513 double phi = -delta_phi/2;
514 for(int i =0; i<nbin; i++){
515 phi += delta_phi;
516 double u = U + R*cos(phi);
517 double v = V + R*sin(phi);
518 double normal = (v-V)/(u-U);
519 double k = -1./normal;
520 double b = v - k*u;
521 double x_cross = -b/(k+1/k);
522 double y_cross = b/(1+k*k);
523
524 double rho=sqrt(x_cross*x_cross+y_cross*y_cross);
525 double theta=atan2(y_cross,x_cross);
526 if(theta<0) {
527 theta=theta+M_PI;
528 rho=-rho;
529 }
530 if( normal ==0 && x>0) {
531 rho = fabs(x);
532 theta = 0;
533 }
534 if( normal ==0 && x<0) {
535 rho = -fabs(x);
536 theta = M_PI;
537 }
538 hough->Fill(theta,rho);
539 }
540}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
nentries
Int_t nentries
Definition: DataBase/tau_mode.c:8
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
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
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: bin_deltaD.cxx:506
bin_deltaD
int bin_deltaD()
Definition: bin_deltaD.cxx:1
intersect_cylinder
double intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
Definition: bin_deltaD.cxx:488