Rev 70 | Rev 73 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
25 | f9daq | 1 | // own include |
2 | //#include "include/raySimulator.h" |
||
3 | #include "include/RTUtil.h" |
||
4 | #include "include/guide.h" |
||
5 | |||
6 | // general include |
||
7 | #include "TROOT.h" |
||
8 | #include "TSystem.h" |
||
9 | #include "TStyle.h" |
||
10 | #include "TCanvas.h" |
||
11 | #include "TMath.h" |
||
12 | #include "TView3D.h" |
||
13 | #include "TH1F.h" |
||
14 | #include "TH2F.h" |
||
15 | #include "TBenchmark.h" |
||
16 | #include "TPolyMarker.h" |
||
17 | #include "TGraph.h" |
||
18 | #include "TF1.h" |
||
19 | |||
20 | //================================================================================= |
||
21 | TCanvas *c3dview; |
||
22 | TCanvas *cacc; |
||
23 | int show_3d = 0, show_data = 0; |
||
24 | int draw_width = 2; |
||
25 | |||
26 | // calls default constructor CVodnik.cpp |
||
27 | |||
28 | //void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0) |
||
29 | //{center.SetXYZ(x,y,z); detector = new CDetector(center);} |
||
30 | |||
31 | |||
32 | //void SetLGType(int in = 1, int side = 1, int out = 0) |
||
33 | //{detector->SetLGType(in, side, out);} |
||
34 | |||
35 | //void SetLG(double SiPM0 = 3.0, double M0 = 1.666, double d0 = 3.0, double n10 = 1.0, double n20 = 1.48, double n30 = 1.48, double R0 = 0.96) |
||
36 | //{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);} |
||
37 | |||
38 | //void SetR(double R0) {detector->SetR(R0);} |
||
39 | /* |
||
40 | void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3) |
||
41 | {detector->SetGlass(glass_on0, glass_d0);} |
||
42 | |||
43 | void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0) |
||
44 | {detector->SetGap(x_gap0 , y_gap0, z_gap0);} |
||
45 | |||
46 | void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6) |
||
47 | {detector->SetRCol(in, lg, out, gla);}; |
||
48 | void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3) |
||
49 | {detector->SetDCol(LG0, glass0, active0);}; |
||
50 | |||
51 | |||
52 | void SetDWidth(int w = 1) {draw_width = w;} |
||
53 | |||
54 | void SetFresnel(int b = 1) {detector->SetFresnel(b);} |
||
55 | void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);} |
||
56 | |||
57 | void SetGuideOn(int b = 0) {detector->SetGuideOn(b);} |
||
58 | |||
59 | void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);} |
||
60 | void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);} |
||
61 | */ |
||
62 | //----------------------------------------------------------------------------- |
||
63 | void Init(double rsys_scale = 7.0) |
||
64 | { |
||
65 | RTSetStyle(gStyle); |
||
66 | gStyle->SetOptStat(10); |
||
67 | |||
68 | if(show_3d) { |
||
69 | double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale}; |
||
70 | double rsys1[]={ rsys_scale, rsys_scale, rsys_scale}; |
||
71 | |||
72 | c3dview = (TCanvas*)gROOT->FindObject("c3dview"); |
||
73 | if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723); c3dview->SetFillColor(0);} |
||
74 | else c3dview->Clear(); |
||
75 | |||
76 | TView3D *view = new TView3D(1, rsys0, rsys1); |
||
77 | view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]); |
||
78 | //view->SetPerspective(); |
||
79 | view->SetParallel(); |
||
80 | view->FrontView(); |
||
81 | view->Zoom();view->Zoom();view->Zoom();//view->Zoom(); |
||
82 | //view->ToggleRulers(); |
||
83 | } |
||
84 | } |
||
85 | //----------------------------------------------------------------------------- |
||
86 | void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0) |
||
87 | { |
||
88 | if(show_data) { |
||
89 | char sbuff[256]; |
||
70 | f9daq | 90 | sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf", |
91 | parameters.getA(), parameters.getLightYield()*acc, parameters.getD(), |
||
25 | f9daq | 92 | parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc); |
93 | |||
94 | if(!only2d) { |
||
95 | RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950); |
||
96 | cdata->Divide(3,3); |
||
97 | int cpc = 1; |
||
98 | cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ"); |
||
99 | cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ"); |
||
71 | f9daq | 100 | cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ"); |
101 | (detector->GetHGlass())->Draw("COLZ"); |
||
25 | f9daq | 102 | |
103 | cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ"); |
||
104 | cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ"); |
||
105 | cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ"); |
||
106 | |||
107 | cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw(); |
||
108 | cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw(); |
||
109 | cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw(); |
||
110 | }else{ |
||
111 | RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500); |
||
112 | cdata->cd(0); (detector->GetHLaser())->Draw("COLZ"); |
||
113 | cdata->SaveAs("cdata.gif"); |
||
114 | } |
||
115 | } |
||
116 | } |
||
117 | //----------------------------------------------------------------------------- |
||
118 | void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance", |
||
119 | double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0) |
||
120 | { |
||
121 | |||
122 | cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480); |
||
123 | |||
124 | (cacc->cd(0))->SetRightMargin(0.05); |
||
125 | (cacc->cd(0))->SetLeftMargin(0.13); |
||
126 | (cacc->cd(0))->SetTopMargin(0.08); |
||
127 | |||
128 | TGraph *gacceptance = new TGraph(n, x, y); |
||
129 | gacceptance->SetTitle(htitle); |
||
130 | gacceptance->SetMarkerStyle(8); |
||
131 | gacceptance->SetLineColor(12); |
||
132 | gacceptance->SetLineWidth(1); |
||
133 | if(xmin < xmax) |
||
134 | (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax); |
||
135 | (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
||
136 | |||
137 | gacceptance->Draw("ACP"); |
||
138 | //cacc->Print("acceptance.pdf"); |
||
139 | //delete gacceptance; |
||
140 | //delete cacc; |
||
141 | } |
||
142 | //----------------------------------------------------------------------------- |
||
143 | void PrintGlassStat(CDetector *detector) |
||
144 | { |
||
145 | int glass_in = ((detector->GetLG())->GetHOut())->GetEntries(); |
||
146 | int glass_out = (detector->GetHGlass())->GetEntries(); |
||
147 | printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in); |
||
148 | } |
||
149 | //----------------------------------------------------------------------------- |
||
150 | void PrintGuideHead() |
||
151 | { |
||
152 | //printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n"); |
||
153 | printf("Acceptance\n"); |
||
154 | } |
||
155 | //----------------------------------------------------------------------------- |
||
54 | f9daq | 156 | void PrintGuideStat(double acceptance) |
25 | f9daq | 157 | { |
158 | /*int n_active = (detector->GetHActive())->GetEntries(); |
||
159 | int n_enter = (detector->GetLG())->GetEnteranceHits(); |
||
160 | double izkoristek = n_active/(double)n_enter; |
||
161 | int fates[7]; (detector->GetLG())->GetVFate(fates); |
||
162 | |||
163 | // printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n"); |
||
164 | |||
165 | printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |", |
||
166 | (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi); |
||
167 | for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter); |
||
168 | |||
169 | printf("%5.1lf\n", 100.0*izkoristek);*/ |
||
170 | |||
54 | f9daq | 171 | //double n_active = (detector->GetHActive())->GetEntries(); |
172 | //double r_acc = 100.0 * n_active / n_in; |
||
173 | double r_acc = 100.0 * acceptance; |
||
25 | f9daq | 174 | printf("%7.3lf\n", r_acc); |
175 | } |
||
176 | //----------------------------------------------------------------------------- |
||
177 | |||
178 | //----------------------------------------------------------------------------- |
||
179 | // en zarek |
||
180 | double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0) |
||
181 | { |
||
182 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
183 | theta = Pi()*theta/180.0; |
||
184 | if(theta < 1e-6) theta = 1e-6; |
||
185 | phi = phi*Pi()/180.0; |
||
70 | f9daq | 186 | if(phi < 1e-6) phi = 1e-6; |
25 | f9daq | 187 | |
188 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
189 | |||
190 | //detector->Init(); |
||
191 | if(show_3d) detector->Draw(draw_width); |
||
192 | |||
193 | CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(), |
||
194 | TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi)); |
||
195 | CRay *ray1 = new CRay; |
||
196 | |||
197 | detector->Propagate(*ray0, ray1, show_3d); |
||
198 | |||
70 | f9daq | 199 | delete ray0; |
200 | delete ray1; |
||
201 | |||
25 | f9daq | 202 | return (detector->GetHActive())->GetEntries() / (double)(1); |
203 | } |
||
204 | //----------------------------------------------------------------------------- |
||
205 | // zarki, razporejeni v mrezi |
||
206 | double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0) |
||
207 | { |
||
208 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
209 | theta = Pi()*theta/180.0; |
||
210 | if(theta < 1e-6) theta = 1e-6; |
||
211 | |||
212 | //detector->Init(); |
||
213 | if(show_3d) detector->Draw(draw_width); |
||
214 | |||
215 | const double b = parameters.getB(); |
||
216 | |||
217 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
218 | |||
219 | for(int i = 0; i < NN; i++) { |
||
220 | for(int j = 0; j < NN; j++) { |
||
221 | CRay *ray0 = new CRay(CENTER.x() - offset, |
||
222 | 0.99*(i*b/NN + b/(2.0*NN) - b/2.0), |
||
223 | 0.99*(j*b/NN + b/(2.0*NN) - b/2.0), |
||
224 | TMath::Cos(theta), |
||
225 | 0.0, |
||
226 | TMath::Sin(theta)); |
||
227 | CRay *ray1 = new CRay; |
||
228 | if(i == (NN/2)) |
||
229 | detector->Propagate(*ray0, ray1, show_3d); |
||
230 | else |
||
70 | f9daq | 231 | detector->Propagate(*ray0, ray1, 0); |
232 | delete ray0; |
||
233 | delete ray1; |
||
25 | f9daq | 234 | } |
70 | f9daq | 235 | |
25 | f9daq | 236 | } |
70 | f9daq | 237 | double acceptance = 0.0; |
238 | /* |
||
239 | if( !(parameters.getPlateOn()) ) |
||
240 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN; |
||
241 | else |
||
242 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
243 | */ |
||
244 | acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN)); |
||
245 | return acceptance; |
||
25 | f9daq | 246 | } |
247 | //----------------------------------------------------------------------------- |
||
248 | // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika) |
||
249 | // vsi pod kotom (theta, phi) |
||
70 | f9daq | 250 | double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phi, int show_rays) |
25 | f9daq | 251 | { |
252 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
253 | theta = theta*3.14159265358979312/180.0; |
||
254 | if(theta < MARGIN) theta = MARGIN; |
||
255 | phi = phi*3.14159265358979312/180.0; |
||
256 | if(phi < MARGIN) phi = MARGIN; |
||
257 | |||
258 | TDatime now; |
||
259 | TRandom rand; |
||
260 | rand.SetSeed(now.Get()); |
||
261 | |||
262 | //detector->Init(); |
||
263 | if(show_3d) detector->Draw(draw_width); |
||
264 | |||
265 | double SiPM = parameters.getA(); |
||
266 | double M = parameters.getM(); |
||
267 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
268 | |||
269 | for(int i = 0; i < NN; i++) { |
||
270 | |||
271 | double rx = CENTER.x() - offset; |
||
272 | double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
273 | double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
274 | |||
275 | double rl = TMath::Cos(theta); |
||
276 | double rm = TMath::Sin(theta)*TMath::Sin(phi); |
||
277 | double rn = TMath::Sin(theta)*TMath::Cos(phi); |
||
278 | |||
279 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
280 | CRay *ray1 = new CRay; |
||
281 | |||
282 | if(i < show_rays) |
||
283 | detector->Propagate(*ray0, ray1, show_3d); |
||
284 | else |
||
285 | detector->Propagate(*ray0, ray1, 0); |
||
54 | f9daq | 286 | //delete ray0; |
287 | //delete ray1; |
||
25 | f9daq | 288 | } |
70 | f9daq | 289 | |
290 | double acceptance = 0.0; |
||
291 | /* |
||
292 | if( !(parameters.getPlateOn()) ) |
||
293 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
294 | else |
||
295 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
296 | //return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/ |
||
297 | acceptance = (detector->GetHActive())->GetEntries() / (double) NN; |
||
298 | return acceptance; |
||
25 | f9daq | 299 | } |
300 | //----------------------------------------------------------------------------- |
||
301 | // zarki, izotropno porazdeljeni znotraj kota theta |
||
302 | // = nakljucni vstopni polozaj in kot |
||
70 | f9daq | 303 | // NN - number of rays to be simulated with angles theta distributed uniformly |
304 | // in cos theta and phi uniformly: |
||
305 | // theta [0, theta] |
||
306 | // phi [0,360] |
||
307 | double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0) |
||
25 | f9daq | 308 | { |
309 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
310 | double pi = 3.14159265358979312; |
||
311 | theta = theta*3.14159265358979312/180.0; |
||
312 | if(theta < MARGIN) theta = MARGIN; |
||
313 | |||
314 | TDatime now; |
||
315 | TRandom rand; |
||
316 | rand.SetSeed(now.Get()); |
||
317 | double rx, ry, rz, rl, rm, rn; |
||
318 | double rphi, rtheta; |
||
319 | //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
||
320 | double theta_min_rad = TMath::Cos(theta); |
||
70 | f9daq | 321 | double theta_max_rad = TMath::Cos(thetaMin); |
25 | f9daq | 322 | |
323 | TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
||
324 | hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
||
325 | hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi); |
||
326 | htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta; |
||
327 | htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0); |
||
328 | hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta; |
||
329 | hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0); |
||
330 | hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl; |
||
331 | hl = new TH1F("hl", "hl", 100, 0.0, 1.0); |
||
332 | hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm; |
||
333 | hm = new TH1F("hm", "hm", 100, 0.0, 1.0); |
||
334 | hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn; |
||
335 | hn = new TH1F("hn", "hn", 100, 0.0, 1.0); |
||
336 | |||
337 | |||
338 | //detector->Init(); |
||
339 | if (show_3d) detector->Draw(draw_width); |
||
340 | |||
341 | double SiPM = parameters.getA(); |
||
342 | double M = parameters.getM(); |
||
343 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
344 | |||
345 | for(int i = 0; i < NN; i++) { |
||
346 | |||
347 | rx = CENTER.x() - offset; |
||
348 | ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
349 | rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
350 | |||
351 | rphi = rand.Uniform(0.0, 2.0*pi); |
||
70 | f9daq | 352 | //double theta_max_rad = TMath::Cos(17.0*pi/180.0); |
353 | //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0))); |
||
354 | rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad))); |
||
25 | f9daq | 355 | |
356 | rl = TMath::Cos(rtheta); |
||
357 | rm = TMath::Sin(rtheta)*TMath::Sin(rphi); |
||
358 | rn = TMath::Sin(rtheta)*TMath::Cos(rphi); |
||
359 | |||
360 | if(show_rand) { |
||
361 | htheta->Fill(rtheta); |
||
362 | hcostheta->Fill( TMath::Cos(rtheta) ); |
||
363 | hphi->Fill(rphi); |
||
364 | hl->Fill(rl); hm->Fill(rm); hn->Fill(rn); |
||
365 | } |
||
366 | |||
367 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
368 | CRay *ray1 = new CRay; |
||
369 | |||
370 | if(i < show_rays) { |
||
371 | detector->Propagate(*ray0, ray1, show_3d); |
||
372 | } |
||
373 | else { |
||
374 | detector->Propagate(*ray0, ray1, 0); |
||
375 | } |
||
376 | |||
70 | f9daq | 377 | delete ray0; |
378 | delete ray1; |
||
25 | f9daq | 379 | } |
380 | |||
381 | if(show_rand) { |
||
382 | TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand"); |
||
383 | if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700); |
||
384 | else c2rand->Clear(); |
||
385 | c2rand->Divide(2,3); |
||
386 | |||
387 | c2rand->cd(1); hphi->Draw(); |
||
388 | c2rand->cd(3); htheta->Draw(); |
||
389 | c2rand->cd(5); hcostheta->Draw(); |
||
390 | c2rand->cd(2); hl->Draw(); |
||
391 | c2rand->cd(4); hm->Draw(); |
||
392 | c2rand->cd(6); hn->Draw(); |
||
393 | } |
||
394 | |||
70 | f9daq | 395 | double acceptance = 0.0; |
396 | /* |
||
397 | if( !(parameters.getPlateOn()) ) |
||
398 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
399 | else |
||
400 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
401 | */ |
||
402 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
403 | return acceptance; |
||
404 | } |
||
405 | |||
406 | // Beamtest distribution |
||
407 | // with fixed theta and phi in interval phiMin, phiMax |
||
408 | double beamtest(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand) |
||
409 | { |
||
410 | double pi = 3.14159265358979312; |
||
411 | theta *= pi/180.0; |
||
412 | if(theta < MARGIN) theta = MARGIN; |
||
413 | phiMin *= pi/180.0; |
||
414 | phiMax *= pi/180.0; |
||
25 | f9daq | 415 | |
70 | f9daq | 416 | TDatime now; |
417 | TRandom rand; |
||
418 | rand.SetSeed(now.Get()); |
||
419 | double rx, ry, rz, rl, rm, rn; |
||
420 | double rphi; |
||
421 | //double rtheta; |
||
422 | //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
||
423 | //double theta_min_rad = TMath::Cos(theta); |
||
424 | |||
425 | TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
||
426 | hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
||
427 | hphi = new TH1F("hphi", "hphi", 100, -pi, pi); |
||
428 | htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta; |
||
429 | htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0); |
||
430 | hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta; |
||
431 | hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0); |
||
432 | hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl; |
||
433 | hl = new TH1F("hl", "hl", 100, 0.0, 1.0); |
||
434 | hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm; |
||
435 | hm = new TH1F("hm", "hm", 100, 0.0, 1.0); |
||
436 | hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn; |
||
437 | hn = new TH1F("hn", "hn", 100, 0.0, 1.0); |
||
438 | |||
439 | //detector->Init(); |
||
440 | if (show_3d) detector->Draw(draw_width); |
||
441 | |||
442 | double b = parameters.getB(); |
||
443 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
444 | |||
445 | for(int i = 0; i < NN; i++) { |
||
446 | |||
447 | rx = CENTER.x() - offset; |
||
448 | ry = rand.Uniform(-b/2.0, b/2.0); |
||
449 | rz = rand.Uniform(-b/2.0, b/2.0); |
||
450 | |||
451 | rphi = rand.Uniform(phiMin, phiMax); |
||
452 | //double theta_max_rad = TMath::Cos(17.0*pi/180.0); |
||
453 | //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0))); |
||
454 | //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad))); |
||
455 | |||
456 | rl = TMath::Cos(theta); |
||
457 | rm = TMath::Sin(theta)*TMath::Sin(rphi); |
||
458 | rn = TMath::Sin(theta)*TMath::Cos(rphi); |
||
459 | |||
460 | if(show_rand) { |
||
461 | //htheta->Fill(rtheta); |
||
462 | //hcostheta->Fill( TMath::Cos(rtheta) ); |
||
463 | hphi->Fill(rphi); |
||
464 | hl->Fill(rl); hm->Fill(rm); hn->Fill(rn); |
||
465 | } |
||
466 | |||
467 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
468 | CRay *ray1 = new CRay; |
||
469 | |||
470 | if(i < show_rays) { |
||
471 | detector->Propagate(*ray0, ray1, show_3d); |
||
472 | } |
||
473 | else { |
||
474 | detector->Propagate(*ray0, ray1, 0); |
||
475 | } |
||
476 | |||
477 | delete ray0; |
||
478 | delete ray1; |
||
479 | } |
||
480 | |||
481 | if(show_rand) { |
||
482 | TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand"); |
||
483 | if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700); |
||
484 | else c2rand->Clear(); |
||
485 | c2rand->Divide(2,3); |
||
486 | |||
487 | c2rand->cd(1); hphi->Draw(); |
||
488 | c2rand->cd(3); htheta->Draw(); |
||
489 | c2rand->cd(5); hcostheta->Draw(); |
||
490 | c2rand->cd(2); hl->Draw(); |
||
491 | c2rand->cd(4); hm->Draw(); |
||
492 | c2rand->cd(6); hn->Draw(); |
||
493 | } |
||
494 | |||
495 | double acceptance = 0.0; |
||
496 | /* |
||
497 | if( !(parameters.getPlateOn()) ) |
||
498 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
499 | else |
||
500 | acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
||
501 | */ |
||
502 | acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
25 | f9daq | 503 | return acceptance; |
504 | } |