Rev 54 | Go to most recent revision | Details | 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]; |
||
90 | sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf", |
||
91 | parameters.getA(), parameters.getM(), parameters.getD(), |
||
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"); |
||
100 | cdata->cd(cpc++); ((detector->GetLG())->GetHOut())->Draw("COLZ"); |
||
101 | |||
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 | //----------------------------------------------------------------------------- |
||
156 | void PrintGuideStat(CDetector *detector, double n_in = 1) |
||
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 | |||
171 | double n_active = (detector->GetHActive())->GetEntries(); |
||
172 | double r_acc = 100.0 * n_active / n_in; |
||
173 | printf("%7.3lf\n", r_acc); |
||
174 | } |
||
175 | //----------------------------------------------------------------------------- |
||
176 | |||
177 | //----------------------------------------------------------------------------- |
||
178 | // en zarek |
||
179 | double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0) |
||
180 | { |
||
181 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
182 | theta = Pi()*theta/180.0; |
||
183 | if(theta < 1e-6) theta = 1e-6; |
||
184 | phi = phi*Pi()/180.0; |
||
185 | |||
186 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
187 | |||
188 | //detector->Init(); |
||
189 | if(show_3d) detector->Draw(draw_width); |
||
190 | |||
191 | CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(), |
||
192 | TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi)); |
||
193 | CRay *ray1 = new CRay; |
||
194 | |||
195 | detector->Propagate(*ray0, ray1, show_3d); |
||
196 | |||
197 | return (detector->GetHActive())->GetEntries() / (double)(1); |
||
198 | } |
||
199 | //----------------------------------------------------------------------------- |
||
200 | // zarki, razporejeni v mrezi |
||
201 | double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0) |
||
202 | { |
||
203 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
204 | theta = Pi()*theta/180.0; |
||
205 | if(theta < 1e-6) theta = 1e-6; |
||
206 | |||
207 | //detector->Init(); |
||
208 | if(show_3d) detector->Draw(draw_width); |
||
209 | |||
210 | const double b = parameters.getB(); |
||
211 | |||
212 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
213 | |||
214 | for(int i = 0; i < NN; i++) { |
||
215 | for(int j = 0; j < NN; j++) { |
||
216 | CRay *ray0 = new CRay(CENTER.x() - offset, |
||
217 | 0.99*(i*b/NN + b/(2.0*NN) - b/2.0), |
||
218 | 0.99*(j*b/NN + b/(2.0*NN) - b/2.0), |
||
219 | TMath::Cos(theta), |
||
220 | 0.0, |
||
221 | TMath::Sin(theta)); |
||
222 | CRay *ray1 = new CRay; |
||
223 | if(i == (NN/2)) |
||
224 | detector->Propagate(*ray0, ray1, show_3d); |
||
225 | else |
||
226 | detector->Propagate(*ray0, ray1, 0); |
||
227 | } |
||
228 | } |
||
229 | |||
230 | return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN)); |
||
231 | } |
||
232 | //----------------------------------------------------------------------------- |
||
233 | // zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika) |
||
234 | // vsi pod kotom (theta, phi) |
||
235 | double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30) |
||
236 | { |
||
237 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
238 | theta = theta*3.14159265358979312/180.0; |
||
239 | if(theta < MARGIN) theta = MARGIN; |
||
240 | phi = phi*3.14159265358979312/180.0; |
||
241 | if(phi < MARGIN) phi = MARGIN; |
||
242 | |||
243 | TDatime now; |
||
244 | TRandom rand; |
||
245 | rand.SetSeed(now.Get()); |
||
246 | |||
247 | //detector->Init(); |
||
248 | if(show_3d) detector->Draw(draw_width); |
||
249 | |||
250 | double SiPM = parameters.getA(); |
||
251 | double M = parameters.getM(); |
||
252 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
253 | |||
254 | for(int i = 0; i < NN; i++) { |
||
255 | |||
256 | double rx = CENTER.x() - offset; |
||
257 | double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
258 | double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
259 | |||
260 | double rl = TMath::Cos(theta); |
||
261 | double rm = TMath::Sin(theta)*TMath::Sin(phi); |
||
262 | double rn = TMath::Sin(theta)*TMath::Cos(phi); |
||
263 | |||
264 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
265 | CRay *ray1 = new CRay; |
||
266 | |||
267 | if(i < show_rays) |
||
268 | detector->Propagate(*ray0, ray1, show_3d); |
||
269 | else |
||
270 | detector->Propagate(*ray0, ray1, 0); |
||
271 | delete ray0; |
||
272 | delete ray1; |
||
273 | } |
||
274 | |||
275 | return (detector->GetHActive())->GetEntries() / (double)NN; |
||
276 | } |
||
277 | //----------------------------------------------------------------------------- |
||
278 | // zarki, izotropno porazdeljeni znotraj kota theta |
||
279 | // = nakljucni vstopni polozaj in kot |
||
280 | double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0) |
||
281 | { |
||
282 | //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
||
283 | double pi = 3.14159265358979312; |
||
284 | theta = theta*3.14159265358979312/180.0; |
||
285 | if(theta < MARGIN) theta = MARGIN; |
||
286 | |||
287 | TDatime now; |
||
288 | TRandom rand; |
||
289 | rand.SetSeed(now.Get()); |
||
290 | double rx, ry, rz, rl, rm, rn; |
||
291 | double rphi, rtheta; |
||
292 | //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
||
293 | double theta_min_rad = TMath::Cos(theta); |
||
294 | |||
295 | TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
||
296 | hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
||
297 | hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi); |
||
298 | htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta; |
||
299 | htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0); |
||
300 | hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta; |
||
301 | hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0); |
||
302 | hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl; |
||
303 | hl = new TH1F("hl", "hl", 100, 0.0, 1.0); |
||
304 | hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm; |
||
305 | hm = new TH1F("hm", "hm", 100, 0.0, 1.0); |
||
306 | hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn; |
||
307 | hn = new TH1F("hn", "hn", 100, 0.0, 1.0); |
||
308 | |||
309 | |||
310 | //detector->Init(); |
||
311 | if (show_3d) detector->Draw(draw_width); |
||
312 | |||
313 | double SiPM = parameters.getA(); |
||
314 | double M = parameters.getM(); |
||
315 | double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
||
316 | |||
317 | for(int i = 0; i < NN; i++) { |
||
318 | |||
319 | rx = CENTER.x() - offset; |
||
320 | ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
321 | rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
||
322 | |||
323 | rphi = rand.Uniform(0.0, 2.0*pi); |
||
324 | rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0))); |
||
325 | |||
326 | rl = TMath::Cos(rtheta); |
||
327 | rm = TMath::Sin(rtheta)*TMath::Sin(rphi); |
||
328 | rn = TMath::Sin(rtheta)*TMath::Cos(rphi); |
||
329 | |||
330 | if(show_rand) { |
||
331 | htheta->Fill(rtheta); |
||
332 | hcostheta->Fill( TMath::Cos(rtheta) ); |
||
333 | hphi->Fill(rphi); |
||
334 | hl->Fill(rl); hm->Fill(rm); hn->Fill(rn); |
||
335 | } |
||
336 | |||
337 | CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
||
338 | CRay *ray1 = new CRay; |
||
339 | |||
340 | if(i < show_rays) { |
||
341 | detector->Propagate(*ray0, ray1, show_3d); |
||
342 | } |
||
343 | else { |
||
344 | detector->Propagate(*ray0, ray1, 0); |
||
345 | } |
||
346 | |||
347 | delete ray0; |
||
348 | delete ray1; |
||
349 | } |
||
350 | |||
351 | if(show_rand) { |
||
352 | TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand"); |
||
353 | if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700); |
||
354 | else c2rand->Clear(); |
||
355 | c2rand->Divide(2,3); |
||
356 | |||
357 | c2rand->cd(1); hphi->Draw(); |
||
358 | c2rand->cd(3); htheta->Draw(); |
||
359 | c2rand->cd(5); hcostheta->Draw(); |
||
360 | c2rand->cd(2); hl->Draw(); |
||
361 | c2rand->cd(4); hm->Draw(); |
||
362 | c2rand->cd(6); hn->Draw(); |
||
363 | } |
||
364 | |||
365 | |||
366 | double acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
||
367 | return acceptance; |
||
368 | } |
||
369 |