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