Rev 73 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 73 | Rev 84 | ||
---|---|---|---|
Line 30... | Line 30... | ||
30 | //{ |
30 | //{ |
31 | //r.SetXYZ(x0, y0, z0); |
31 | //r.SetXYZ(x0, y0, z0); |
32 | //n.SetXYZ(l0, m0, n0); n = n.Unit(); |
32 | //n.SetXYZ(l0, m0, n0); n = n.Unit(); |
33 | //} |
33 | //} |
34 | //----------------------------------------------------------------------------- |
34 | //----------------------------------------------------------------------------- |
35 | /* |
35 | |
36 |
|
36 | CRay& CRay::operator = (const CRay& rhs) |
37 | { |
37 | { |
38 |
|
38 | r.SetXYZ(rhs.GetR().x(), rhs.GetR().y(), rhs.GetR().z()); |
39 | //this->r.SetXYZ(p.x(), p.y(), p.z()); |
39 | //this->r.SetXYZ(p.x(), p.y(), p.z()); |
- | 40 | k.SetXYZ(rhs.GetK().x(), rhs.GetK().y(), rhs.GetK().z()); |
|
40 |
|
41 | p.SetXYZ(rhs.GetP().x(), rhs.GetP().y(), rhs.GetP().z()); |
41 | return *this; |
42 | return *this; |
42 |
|
43 | } |
43 | void CRay::Print() |
44 | void CRay::Print() |
44 | { |
45 | { |
45 | printf("---> CRay::Print() <---\n"); |
46 | printf("---> CRay::Print() <---\n"); |
46 | printf("(x,y,z)=(%.2lf, %.2lf, %.2lf); (l,m,n)=(%.2lf, %.2lf, %.2lf)\n", |
47 | printf("(x,y,z)=(%.2lf, %.2lf, %.2lf); (l,m,n)=(%.2lf, %.2lf, %.2lf)\n", |
47 | r.x(), r.y(), r.z(), k.x(), k.y(), k.z()); |
48 | r.x(), r.y(), r.z(), k.x(), k.y(), k.z()); |
48 | } |
49 | } |
49 | void CRay::Draw() |
50 | void CRay::Draw() |
50 | { |
51 | { |
51 | double t = 50.0; |
52 | double t = 50.0; |
52 | TPolyLine3D *line3d = new TPolyLine3D(2); |
53 | TPolyLine3D *line3d = new TPolyLine3D(2); |
53 | //line3d->SetPoint(0, r.x() - t*n.x(), r.y() - t*n.y(), r.z() - t*n.z()); |
54 | //line3d->SetPoint(0, r.x() - t*n.x(), r.y() - t*n.y(), r.z() - t*n.z()); |
54 | line3d->SetPoint(0, r.x(), r.y(), r.z()); |
55 | line3d->SetPoint(0, r.x(), r.y(), r.z()); |
55 | line3d->SetPoint(1, r.x() + t*k.x(), k.y() + t*k.y(), r.z() + t*k.z()); |
56 | line3d->SetPoint(1, r.x() + t*k.x(), k.y() + t*k.y(), r.z() + t*k.z()); |
56 | line3d->SetLineWidth(1); |
57 | line3d->SetLineWidth(1); |
57 | line3d->SetLineColor(color); |
58 | line3d->SetLineColor(color); |
Line 92... | Line 93... | ||
92 | line3d->SetLineWidth(1); |
93 | line3d->SetLineWidth(1); |
93 | line3d->SetLineColor(color); |
94 | line3d->SetLineColor(color); |
94 | 95 | ||
95 | line3d->Draw(); |
96 | line3d->Draw(); |
96 | } |
97 | } |
97 | 98 | ||
98 | 99 | ||
99 | CPlane4::CPlane4() : |
100 | CPlane4::CPlane4() : |
100 | n(TVector3(1.0, 0.0, 0.0)), |
101 | n(TVector3(1.0, 0.0, 0.0)), |
101 | A(0), |
102 | A(0), |
102 | B(0), |
103 | B(0), |
Line 116... | Line 117... | ||
116 | //----------------------------------------------------------------------------- |
117 | //----------------------------------------------------------------------------- |
117 | // za izracun parametrov ravnine je en vektor prevec, vendar tega |
118 | // za izracun parametrov ravnine je en vektor prevec, vendar tega |
118 | // rabim kot zadnji vogal poligona |
119 | // rabim kot zadnji vogal poligona |
119 | //void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4) |
120 | //void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4) |
120 | //{ |
121 | //{ |
121 | double x1,y1,z1, x2,y2,z2, x3,y3,z3; |
122 | double x1,y1,z1, x2,y2,z2, x3,y3,z3; |
122 | 123 | ||
123 | x1 = r1.x(); y1 = r1.y(); z1 = r1.z(); |
124 | x1 = r1.x(); y1 = r1.y(); z1 = r1.z(); |
124 | x2 = r2.x(); y2 = r2.y(); z2 = r2.z(); |
125 | x2 = r2.x(); y2 = r2.y(); z2 = r2.z(); |
125 | x3 = r3.x(); y3 = r3.y(); z3 = r3.z(); |
126 | x3 = r3.x(); y3 = r3.y(); z3 = r3.z(); |
126 | 127 | ||
127 | A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1); |
128 | A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1); |
128 | B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3); |
129 | B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3); |
129 | C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1); |
130 | C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1); |
130 | D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2); |
131 | D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2); |
131 | 132 | ||
132 | r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4; |
133 | r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4; |
133 | n.SetXYZ(A, B, C); |
134 | n.SetXYZ(A, B, C); |
134 | n = n.Unit(); |
135 | n = n.Unit(); |
135 | 136 | ||
136 | for(int i=0;i<4;i++) |
137 | for(int i=0;i<4;i++) |
137 | edge[i] = r[i-3 ? i+1 : 0] - r[i]; |
138 | edge[i] = r[i-3 ? i+1 : 0] - r[i]; |
138 | 139 | ||
139 | for(int i=0;i<4;i++) |
140 | for(int i=0;i<4;i++) |
140 | angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) )); |
141 | angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) )); |
141 | }; |
142 | }; |
142 | 143 | ||
143 | void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4) |
144 | void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4) |
144 | { |
145 | { |
145 | double x1,y1,z1, x2,y2,z2, x3,y3,z3; |
146 | double x1,y1,z1, x2,y2,z2, x3,y3,z3; |
146 | 147 | ||
147 | x1 = r1.x(); y1 = r1.y(); z1 = r1.z(); |
148 | x1 = r1.x(); y1 = r1.y(); z1 = r1.z(); |
148 | x2 = r2.x(); y2 = r2.y(); z2 = r2.z(); |
149 | x2 = r2.x(); y2 = r2.y(); z2 = r2.z(); |
149 | x3 = r3.x(); y3 = r3.y(); z3 = r3.z(); |
150 | x3 = r3.x(); y3 = r3.y(); z3 = r3.z(); |
150 | 151 | ||
151 | A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1); |
152 | A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1); |
152 | B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3); |
153 | B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3); |
153 | C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1); |
154 | C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1); |
154 | D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2); |
155 | D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2); |
155 | 156 | ||
156 | r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4; |
157 | r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4; |
157 | n.SetXYZ(A, B, C); |
158 | n.SetXYZ(A, B, C); |
158 | n = n.Unit(); |
159 | n = n.Unit(); |
159 | 160 | ||
160 | for(int i=0;i<4;i++) |
161 | for(int i=0;i<4;i++) |
161 | edge[i] = r[i-3 ? i+1 : 0] - r[i]; |
162 | edge[i] = r[i-3 ? i+1 : 0] - r[i]; |
162 | 163 | ||
163 | for(int i=0;i<4;i++) |
164 | for(int i=0;i<4;i++) |
164 | angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) )); |
165 | angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) )); |
Line 168... | Line 169... | ||
168 | { |
169 | { |
169 | double x1,y1,z1, x2,y2,z2, x3,y3,z3; |
170 | double x1,y1,z1, x2,y2,z2, x3,y3,z3; |
170 | 171 | ||
171 | x1 = vr[0].x(); y1 = vr[0].y(); z1 = vr[0].z(); |
172 | x1 = vr[0].x(); y1 = vr[0].y(); z1 = vr[0].z(); |
172 | x2 = vr[1].x(); y2 = vr[1].y(); z2 = vr[1].z(); |
173 | x2 = vr[1].x(); y2 = vr[1].y(); z2 = vr[1].z(); |
173 | x3 = vr[2].x(); y3 = vr[2].y(); z3 = vr[2].z(); |
174 | x3 = vr[2].x(); y3 = vr[2].y(); z3 = vr[2].z(); |
174 | 175 | ||
175 | A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1); |
176 | A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1); |
176 | B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3); |
177 | B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3); |
177 | C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1); |
178 | C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1); |
178 | D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2); |
179 | D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2); |
179 | 180 | ||
Line 199... | Line 200... | ||
199 | 200 | ||
200 | N.SetXYZ(A,B,C); |
201 | N.SetXYZ(A,B,C); |
201 | 202 | ||
202 | num = N*ray.GetR() + D; |
203 | num = N*ray.GetR() + D; |
203 | den = N*ray.GetK(); |
204 | den = N*ray.GetK(); |
204 | 205 | ||
205 | if (dbg) printf("t = %6.3lf / %6.3lf = %6.3lf\n", num, den, num/den); |
206 | if (dbg) printf("t = %6.3lf / %6.3lf = %6.3lf\n", num, den, num/den); |
206 | 207 | ||
207 | //if(den == 0) |
208 | //if(den == 0) |
208 | if(TMath::Abs(den) < MARGIN) { |
209 | if(TMath::Abs(den) < MARGIN) { |
209 | //if(num == 0) |
210 | //if(num == 0) |
Line 214... | Line 215... | ||
214 | else { |
215 | else { |
215 | if (dbg) printf("The ray is parallel to the surface!\n"); |
216 | if (dbg) printf("The ray is parallel to the surface!\n"); |
216 | return 0; // ni presecisca |
217 | return 0; // ni presecisca |
217 | } |
218 | } |
218 | } |
219 | } |
219 | 220 | ||
220 | t = num / den; |
221 | t = num / den; |
221 | 222 | ||
222 | tmp = ray.GetR(); |
223 | tmp = ray.GetR(); |
223 | tmp -= t*ray.GetK(); |
224 | tmp -= t*ray.GetK(); |
224 | *vec = tmp; |
225 | *vec = tmp; |
Line 271... | Line 272... | ||
271 | TVector3 tmp; |
272 | TVector3 tmp; |
272 | 273 | ||
273 | if( GetIntersection(&tmp, in) ) |
274 | if( GetIntersection(&tmp, in) ) |
274 | if( IsVectorIn(tmp) ) |
275 | if( IsVectorIn(tmp) ) |
275 | return 1; |
276 | return 1; |
276 | 277 | ||
277 | return 0; |
278 | return 0; |
278 | } |
279 | } |
279 | int CPlane4::TestIntersection(TVector3 *vec, CRay in) |
280 | int CPlane4::TestIntersection(TVector3 *vec, CRay in) |
280 | { |
281 | { |
281 | TVector3 tmp; |
282 | TVector3 tmp; |
282 | 283 | ||
283 | if( GetIntersection(&tmp, in) ) |
284 | if( GetIntersection(&tmp, in) ) |
Line 362... | Line 363... | ||
362 | } |
363 | } |
363 | //----------------------------------------------------------------------------- |
364 | //----------------------------------------------------------------------------- |
364 | // sprejme zarek, vrne uklonjen/odbit zarek in presecisce |
365 | // sprejme zarek, vrne uklonjen/odbit zarek in presecisce |
365 | // vrne 0 ce ni presecisca; 1 ce se je lomil |
366 | // vrne 0 ce ni presecisca; 1 ce se je lomil |
366 | // 2 ce se je odbil; -2 ce se je absorbiral |
367 | // 2 ce se je odbil; -2 ce se je absorbiral |
367 | int CSurface::PropagateRay(CRay in, CRay |
368 | int CSurface::PropagateRay(CRay in, CRay& out, TVector3 *intersection) |
368 | { |
369 | { |
369 | if (dbg) printf("--- CSurface::PropagateRay ---\n"); |
370 | if (dbg) printf("--- CSurface::PropagateRay ---\n"); |
370 | double cosTi; // incident ray angle |
371 | double cosTi; // incident ray angle |
371 | double cosTt; // transmited ray angle |
372 | double cosTt; // transmited ray angle |
372 | TVector3 intersect, transmit; |
373 | TVector3 intersect, transmit; |
Line 413... | Line 414... | ||
413 | printf(" v_te = "); printv(v_te); |
414 | printf(" v_te = "); printv(v_te); |
414 | printf(" v_tm = "); printv(v_tm); |
415 | printf(" v_tm = "); printv(v_tm); |
415 | } |
416 | } |
416 | 417 | ||
417 | double cosAf = v_te * in.GetP(); |
418 | double cosAf = v_te * in.GetP(); |
418 | double alpha = |
419 | double alpha = TMath::ACos(cosAf); |
- | 420 | if(alpha < MARGIN) alpha = MARGIN; |
|
419 | if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE); |
421 | if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE); |
420 | 422 | ||
421 | a_te = cosAf; |
423 | a_te = cosAf; |
422 | a_tm = TMath::Sqrt(1 - cosAf*cosAf); |
424 | a_tm = TMath::Sqrt(1 - cosAf*cosAf); |
423 | if(dbg) printf(" a_te = %lf, a_tm = %lf\n", a_te, a_tm); |
425 | if(dbg) printf(" a_te = %lf, a_tm = %lf\n", a_te, a_tm); |
Line 457... | Line 459... | ||
457 | 459 | ||
458 | if(dbg) { |
460 | if(dbg) { |
459 | cosTN = TMath::Abs(transmit.Unit() * n); |
461 | cosTN = TMath::Abs(transmit.Unit() * n); |
460 | printf(" cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE); |
462 | printf(" cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE); |
461 | } |
463 | } |
462 | out |
464 | out.Set(intersect, transmit); |
463 | 465 | ||
464 | // Shift implemented, but only linear polarization is implemented |
466 | // Shift implemented, but only linear polarization is implemented |
465 | if (dbg) printf("CSurface: Propagate TOTAL\n"); |
467 | if (dbg) printf("CSurface: Propagate TOTAL\n"); |
466 | v_tm_t = -v_te.Cross(transmit); |
468 | v_tm_t = -v_te.Cross(transmit); |
467 | v_tm_t = v_tm_t.Unit(); |
469 | v_tm_t = v_tm_t.Unit(); |
468 | // shift the p and s components |
470 | // shift the p and s components |
469 | double n12 = N1_N2(-sign_n); |
471 | double n12 = N1_N2(-sign_n); |
470 | double deltaP = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/(pow(n12,2)*cosTi)); |
472 | double deltaP = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/(pow(n12,2)*cosTi)); |
471 | double deltaS = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/cosTi); |
473 | double deltaS = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/cosTi); |
472 | double delta = deltaP - deltaS; |
474 | double delta = deltaP - deltaS; |
473 | alpha |
475 | //alpha += delta; |
474 | a_tm = sin(alpha); |
476 | a_tm = sin(alpha); |
475 | a_te = cos(alpha); |
477 | a_te = cos(alpha); |
476 | if (dbg) printf(" deltaP = %f deltaS = %f; new a_tm = %f, a_te = %f", |
478 | if (dbg) printf(" deltaP = %f deltaS = %f; new a_tm = %f, a_te = %f", |
477 | deltaP, deltaS, a_tm, a_te); |
479 | deltaP, deltaS, a_tm, a_te); |
478 | pol_t = a_tm*v_tm_t + a_te*v_te; |
480 | pol_t = a_tm*v_tm_t + a_te*v_te; |
479 | if (dbg) printv(pol_t); |
481 | if (dbg) printv(pol_t); |
480 | out |
482 | out.setPolarization(pol_t); |
481 | 483 | ||
482 | return REFLECTION; |
484 | return REFLECTION; |
483 | } else { |
485 | } else { |
484 | // reflection or refraction according to Fresnel equations |
486 | // reflection or refraction according to Fresnel equations |
485 | if(dbg) printf(" REFRACTION\n"); |
487 | if(dbg) printf(" REFRACTION\n"); |
Line 519... | Line 521... | ||
519 | if (dbg) printf(" R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f); |
521 | if (dbg) printf(" R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f); |
520 | } |
522 | } |
521 | 523 | ||
522 | if(p_ref >= R_f) { // se lomi |
524 | if(p_ref >= R_f) { // se lomi |
523 | if (dbg) printf(" SURFACE REFRACTED. Return.\n"); |
525 | if (dbg) printf(" SURFACE REFRACTED. Return.\n"); |
524 | out |
526 | out.Set(intersect, transmit); |
525 | out |
527 | out.setPolarization(pol_t); |
526 | return REFRACTION; |
528 | return REFRACTION; |
527 | } else { // se odbije |
529 | } else { // se odbije |
528 | if (dbg) printf(" SURFACE REFLECTED. p_ref=%f, R_f=%f\n", p_ref, R_f); |
530 | if (dbg) printf(" SURFACE REFLECTED. p_ref=%f, R_f=%f\n", p_ref, R_f); |
529 | transmit = in.GetK() + sign_n*2*cosTi*n; |
531 | transmit = in.GetK() + sign_n*2*cosTi*n; |
530 | TVector3 v_tm_r = -v_te.Cross(transmit); |
532 | TVector3 v_tm_r = -v_te.Cross(transmit); |
531 | v_tm_r = v_tm_r.Unit(); |
533 | v_tm_r = v_tm_r.Unit(); |
532 | out |
534 | out.Set(intersect, transmit); |
533 | //pol_t = -in.GetP() + sign_n*2*cosTi*n; |
535 | //pol_t = -in.GetP() + sign_n*2*cosTi*n; |
534 | pol_t = a_te*(1.0 - TMath::Abs(r_te))*v_te + a_tm*(1.0 - TMath::Abs(r_tm))*v_tm_r; |
536 | pol_t = a_te*(1.0 - TMath::Abs(r_te))*v_te + a_tm*(1.0 - TMath::Abs(r_tm))*v_tm_r; |
535 | out |
537 | out.setPolarization(pol_t); |
536 | return REFLECTION; |
538 | return REFLECTION; |
537 | } |
539 | } |
538 | 540 | ||
539 | //} |
541 | //} |
540 | break; |
542 | break; |
Line 545... | Line 547... | ||
545 | case SURF_REFLE: |
547 | case SURF_REFLE: |
546 | p_ref = rand.Uniform(0.0, 1.0); |
548 | p_ref = rand.Uniform(0.0, 1.0); |
547 | if(p_ref < reflection) { // se odbije |
549 | if(p_ref < reflection) { // se odbije |
548 | cosTi = in.GetK() * n; |
550 | cosTi = in.GetK() * n; |
549 | transmit = in.GetK() - 2*cosTi*n; |
551 | transmit = in.GetK() - 2*cosTi*n; |
550 | out |
552 | out.Set(intersect, transmit); |
551 | return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb |
553 | return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb |
552 | } else { // se ne odbije |
554 | } else { // se ne odbije |
553 | transmit = in.GetK(); |
555 | transmit = in.GetK(); |
554 | out |
556 | out.Set(intersect, transmit); |
555 | return ABSORBED; |
557 | return ABSORBED; |
556 | } |
558 | } |
557 | break; |
559 | break; |
558 | 560 | ||
559 | // total reflection from n1 to n2 with R probbability |
561 | // total reflection from n1 to n2 with R probbability |
Line 561... | Line 563... | ||
561 | p_ref = rand.Uniform(0.0, 1.0); |
563 | p_ref = rand.Uniform(0.0, 1.0); |
562 | if(p_ref < reflection) { // se odbije |
564 | if(p_ref < reflection) { // se odbije |
563 | cosTi = in.GetK() * n; |
565 | cosTi = in.GetK() * n; |
564 | if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj |
566 | if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj |
565 | transmit = in.GetK() - 2*cosTi*n; |
567 | transmit = in.GetK() - 2*cosTi*n; |
566 | out |
568 | out.Set(intersect, transmit); |
567 | } else { // ni tot. odboja |
569 | } else { // ni tot. odboja |
568 | transmit = in.GetK(); |
570 | transmit = in.GetK(); |
569 | out |
571 | out.Set(intersect, transmit); |
570 | return ABSORBED; |
572 | return ABSORBED; |
571 | } |
573 | } |
572 | } else { // se ne odbije |
574 | } else { // se ne odbije |
573 | transmit = in.GetK(); |
575 | transmit = in.GetK(); |
574 | out |
576 | out.Set(intersect, transmit); |
575 | return ABSORBED; |
577 | return ABSORBED; |
576 | } |
578 | } |
577 | break; |
579 | break; |
578 | 580 | ||
579 | default: |
581 | default: |
580 |
|
582 | out = in; |
581 | break; |
583 | break; |
582 | } |
584 | } |
583 | 585 | ||
584 | return REFRACTION; |
586 | return REFRACTION; |
585 | } |
587 | } |
Line 618... | Line 620... | ||
618 | vodnik_edge[7].SetXYZ(_d,-t,-t); |
620 | vodnik_edge[7].SetXYZ(_d,-t,-t); |
619 | 621 | ||
620 | for(int i = 0; i<8; i++) vodnik_edge[i] += center; |
622 | for(int i = 0; i<8; i++) vodnik_edge[i] += center; |
621 | 623 | ||
622 | // light guide surfaces |
624 | // light guide surfaces |
623 | s_side[0] |
625 | s_side[0].Set(SURF_REFRA, vodnik_edge, n0, _n2, _r); |
624 | s_side[0] |
626 | s_side[0].FlipN(); |
625 | 627 | ||
626 | s_side[1] |
628 | s_side[1].Set(SURF_REFRA, vodnik_edge[3], vodnik_edge[2], |
627 | vodnik_edge[6], vodnik_edge[7], _n2, _n1, _r); |
629 | vodnik_edge[6], vodnik_edge[7], _n2, _n1, _r); |
628 | s_side[2] |
630 | s_side[2].Set(SURF_REFRA, vodnik_edge[2], vodnik_edge[1], |
629 | vodnik_edge[5], vodnik_edge[6], _n2, _n1, _r); |
631 | vodnik_edge[5], vodnik_edge[6], _n2, _n1, _r); |
630 | s_side[3] |
632 | s_side[3].Set(SURF_REFRA, vodnik_edge[1], vodnik_edge[0], |
631 | vodnik_edge[4], vodnik_edge[5], _n2, _n1, _r); |
633 | vodnik_edge[4], vodnik_edge[5], _n2, _n1, _r); |
632 | s_side[4] |
634 | s_side[4].Set(SURF_REFRA, vodnik_edge[0], vodnik_edge[3], |
633 | vodnik_edge[7], vodnik_edge[4], _n2, _n1, _r); |
635 | vodnik_edge[7], vodnik_edge[4], _n2, _n1, _r); |
634 | // n3 - ref ind at the exit, grease, air |
636 | // n3 - ref ind at the exit, grease, air |
635 | s_side[5] |
637 | s_side[5].Set(SURF_REFRA, &vodnik_edge[4], _n2, _n3, _r); |
636 | s_side[5] |
638 | s_side[5].FlipN(); |
637 | // exit surface in the case of bad coupling |
639 | // exit surface in the case of bad coupling |
638 | noCoupling |
640 | noCoupling.Set(SURF_REFRA, &vodnik_edge[4], _n2, 1.0, _r); |
639 | noCoupling |
641 | noCoupling.FlipN(); |
640 | // grease = specific pattern area of coupling |
642 | // grease = specific pattern area of coupling |
641 | TVector3 activePosition(center); |
643 | TVector3 activePosition(center); |
642 | activePosition += TVector3(_d, 0, 0); |
644 | activePosition += TVector3(_d, 0, 0); |
643 | TVector3 normal(1,0,0); |
645 | TVector3 normal(1,0,0); |
644 | grease |
646 | grease.Set(activePosition, normal, 0.95*a/2.0); |
645 | 647 | ||
646 | if(fresnel) for(int i=0; i<6; i++) s_side[i] |
648 | if(fresnel) for(int i=0; i<6; i++) s_side[i].SetFresnel(1); |
647 | 649 | ||
648 | // statistics histograms |
650 | // statistics histograms |
649 | hfate |
651 | //hfate = (TH1F)gROOT->FindObject("hfate"); //if(hfate) delete hfate; |
650 | hfate = |
652 | hfate = TH1F("hfate", "Ray fate", 8, -3.5, 4.5); |
651 | (hfate |
653 | (hfate.GetXaxis())->SetBinLabel(1, "Back Refl"); |
652 | (hfate |
654 | (hfate.GetXaxis())->SetBinLabel(2, "Side"); |
653 | (hfate |
655 | (hfate.GetXaxis())->SetBinLabel(3, "Refracted"); |
654 | (hfate |
656 | (hfate.GetXaxis())->SetBinLabel(4, "LG Miss"); |
655 | (hfate |
657 | (hfate.GetXaxis())->SetBinLabel(5, "Exit"); |
656 | (hfate |
658 | (hfate.GetXaxis())->SetBinLabel(6, "Enter"); |
657 | (hfate |
659 | (hfate.GetXaxis())->SetBinLabel(7, "Simulated"); |
658 | (hfate |
660 | (hfate.GetXaxis())->SetBinLabel(8, "Absorb"); |
659 | 661 | ||
660 | hnodb_all |
662 | //hnodb_all = (TH1F*)gROOT->FindObject("hnodb_all"); if(hnodb_all) delete hnodb_all; |
661 | hnodb_all = |
663 | hnodb_all = TH1F("hnodb_all", "", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5); |
662 | 664 | ||
663 | hnodb_exit |
665 | //hnodb_exit = (TH1F*)gROOT->FindObject("hnodb_exit"); if(hnodb_exit) delete hnodb_exit; |
664 | hnodb_exit = |
666 | hnodb_exit = TH1F("hnodb_exit", "", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5); |
665 | 667 | ||
666 | int nBins = nch |
668 | int nBins = nch; |
667 | hin |
669 | //hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin; |
668 | hin = |
670 | hin = TH2F("hin", ";x [mm]; y[mm]", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0); |
669 | 671 | ||
670 | hout |
672 | //hout = (TH2F*)gROOT->FindObject("hout"); if(hout) delete hout; |
671 | hout = |
673 | hout = TH2F("hout", ";x [mm];y [mm]", nBins, -a/2.0, +a/2.0, nBins, -a/2.0, +a/2.0); |
672 | } |
674 | } |
673 | //----------------------------------------------------------------------------- |
675 | //----------------------------------------------------------------------------- |
674 | // Sledi zarku skozi vodnik. Vrne: |
676 | // Sledi zarku skozi vodnik. Vrne: |
675 | // 0, ce zgresi vstopno ploskev |
677 | // 0, ce zgresi vstopno ploskev |
676 | // 1, ce zadane izstopno ploskev |
678 | // 1, ce zadane izstopno ploskev |
677 | // -1, ce se v vodniku ne odbije totalno |
679 | // -1, ce se v vodniku ne odbije totalno |
678 | // 2, enter the light guide, bin 2 of hfate = refraction |
680 | // 2, enter the light guide, bin 2 of hfate = refraction |
679 | // -2, ce se ne odbije zaradi koncnega R stranic |
681 | // -2, ce se ne odbije zaradi koncnega R stranic |
680 | // -3, ce se odbije nazaj in gre nazaj ven skozi sprednjo ploskev |
682 | // -3, ce se odbije nazaj in gre nazaj ven skozi sprednjo ploskev |
681 | // +4, ce se absorbira v materialu |
683 | // +4, ce se absorbira v materialu |
682 | Fate Guide::PropagateRay(CRay in, CRay |
684 | Fate Guide::PropagateRay(CRay in, CRay& out, int *n_points, TVector3 *points) |
683 | { |
685 | { |
684 | if (dbg) printf("--- GUIDE::PropagateRay ---\n"); |
686 | if (dbg) printf("--- GUIDE::PropagateRay ---\n"); |
685 | // ray0 - incident ray |
687 | // ray0 - incident ray |
686 | // ray1 - trans/refl ray |
688 | // ray1 - trans/refl ray |
687 | CRay ray0; |
689 | CRay ray0; |
Line 691... | Line 693... | ||
691 | 693 | ||
692 | ray0 = in; |
694 | ray0 = in; |
693 | int n_odb = 0; |
695 | int n_odb = 0; |
694 | int last_hit = 0; |
696 | int last_hit = 0; |
695 | int propagation = 0; |
697 | int propagation = 0; |
696 | int result = s_side[0] |
698 | int result = s_side[0].PropagateRay(ray0, ray1, &vec1); |
697 | if( !(result) ) { |
699 | if( !(result) ) { |
698 | // ce -NI- presecisca z vstopno |
700 | // ce -NI- presecisca z vstopno |
699 | if (dbg) printf(" GUIDE: missed the light guide\n"); |
701 | if (dbg) printf(" GUIDE: missed the light guide\n"); |
700 | fate = missed; |
702 | fate = missed; |
701 | //hfate->Fill(0); |
703 | //hfate->Fill(0); |
Line 704... | Line 706... | ||
704 | fate = backreflected; |
706 | fate = backreflected; |
705 | //hfate->Fill(-3); |
707 | //hfate->Fill(-3); |
706 | } else { |
708 | } else { |
707 | if (dbg) printf(" GUIDE: ray entered\n"); |
709 | if (dbg) printf(" GUIDE: ray entered\n"); |
708 | points[0] = ray1.GetR(); |
710 | points[0] = ray1.GetR(); |
709 | hfate |
711 | hfate.Fill(enter); // enter |
710 | hin |
712 | hin.Fill(vec1.y(), vec1.z()); |
711 | if (dbg) printf(" GUIDE: n_odb = %d\n", n_odb); |
713 | if (dbg) printf(" GUIDE: n_odb = %d\n", n_odb); |
712 | 714 | ||
713 | while (n_odb++ < MAX_REFLECTIONS) { |
715 | while (n_odb++ < MAX_REFLECTIONS) { |
714 | if (dbg) printf(" GUIDE: Boundary test: %d\n",n_odb); |
716 | if (dbg) printf(" GUIDE: Boundary test: %d\n",n_odb); |
715 | ray0 = ray1; |
717 | ray0 = ray1; |
716 | vec0 = vec1; |
718 | vec0 = vec1; |
717 | propagation = 11; |
719 | propagation = 11; |
718 | for(inters_i=0; inters_i<6; inters_i++) { |
720 | for(inters_i=0; inters_i<6; inters_i++) { |
719 | if (dbg) printf(" GUIDE: Test intersection with surface %d \n", inters_i); |
721 | if (dbg) printf(" GUIDE: Test intersection with surface %d \n", inters_i); |
720 | if( inters_i != last_hit) { |
722 | if( inters_i != last_hit) { |
721 | int testBoundary = s_side[inters_i] |
723 | int testBoundary = s_side[inters_i].TestIntersection(&vec1, ray1); |
722 | if( testBoundary ) { |
724 | if( testBoundary ) { |
723 | if (dbg) printf(" GUIDE: ray intersects with LG surface %d\n",inters_i); |
725 | if (dbg) printf(" GUIDE: ray intersects with LG surface %d\n",inters_i); |
724 | break; |
726 | break; |
725 | } |
727 | } |
726 | } |
728 | } |
Line 731... | Line 733... | ||
731 | //hfate->Fill(backreflected); |
733 | //hfate->Fill(backreflected); |
732 | break; |
734 | break; |
733 | } // backreflection |
735 | } // backreflection |
734 | 736 | ||
735 | // the passage is possible, test propagation |
737 | // the passage is possible, test propagation |
736 | propagation = s_side[inters_i] |
738 | propagation = s_side[inters_i].PropagateRay(ray0, ray1, &vec1); |
737 | 739 | ||
738 | if (dbg) printf(" GUIDE: surface = %d, propagation = %d\n", inters_i, propagation); |
740 | if (dbg) printf(" GUIDE: surface = %d, propagation = %d\n", inters_i, propagation); |
739 | 741 | ||
740 | 742 | ||
741 | if(propagation == ABSORBED) { |
743 | if(propagation == ABSORBED) { |
Line 744... | Line 746... | ||
744 | } //refraction due to finite reflectivity |
746 | } //refraction due to finite reflectivity |
745 | 747 | ||
746 | if(inters_i == 5) { |
748 | if(inters_i == 5) { |
747 | if (_badCoupling) { |
749 | if (_badCoupling) { |
748 | TVector3 hitVector(0,0,0); |
750 | TVector3 hitVector(0,0,0); |
749 | bool hitActive = grease |
751 | bool hitActive = grease.TestIntersection(&hitVector, ray0); |
750 | if (hitActive and dbg) printf(" GUIDE: hit grease\n"); |
752 | if (hitActive and dbg) printf(" GUIDE: hit grease\n"); |
751 | if (!hitActive) propagation = noCoupling |
753 | if (!hitActive) propagation = noCoupling.PropagateRay(ray0, ray1, &vec1); |
752 | } |
754 | } |
753 | // check on which side the vector is? |
755 | // check on which side the vector is? |
754 | TVector3 ray = ray1.GetK(); |
756 | TVector3 ray = ray1.GetK(); |
755 | TVector3 exitNormal = s_side[5] |
757 | TVector3 exitNormal = s_side[5].GetN(); |
756 | if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal); |
758 | if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal); |
757 | if (ray*exitNormal > 0) { |
759 | if (ray*exitNormal > 0) { |
758 | if (dbg) printf(" GUIDE: ray is backreflected from exit window.\n"); |
760 | if (dbg) printf(" GUIDE: ray is backreflected from exit window.\n"); |
759 | fate = backreflected; |
761 | fate = backreflected; |
760 | n_odb++; |
762 | n_odb++; |
761 | points[n_odb] = vec1; |
763 | points[n_odb] = vec1; |
762 | ray0 = ray1; |
764 | ray0 = ray1; |
763 | break; |
765 | break; |
764 | } |
766 | } |
765 | fate = hitExit; |
767 | fate = hitExit; |
766 | hout |
768 | hout.Fill(vec1.y(), vec1.z()); |
767 | hnodb_exit |
769 | hnodb_exit.Fill(n_odb-1); |
768 | n_odb++; |
770 | n_odb++; |
769 | points[n_odb] = vec1; |
771 | points[n_odb] = vec1; |
770 | ray0 = ray1; |
772 | ray0 = ray1; |
771 | break; |
773 | break; |
772 | } |
774 | } |
Line 798... | Line 800... | ||
798 | 800 | ||
799 | if(p_abs > T_abs) fate = absorbed; // absorption |
801 | if(p_abs > T_abs) fate = absorbed; // absorption |
800 | } |
802 | } |
801 | //--- material absorption --- |
803 | //--- material absorption --- |
802 | 804 | ||
803 | hfate |
805 | hfate.Fill(fate); |
804 | hfate |
806 | hfate.Fill(rays); |
805 | hnodb_all |
807 | hnodb_all.Fill(n_odb-2); |
806 | *n_points = n_odb+1; |
808 | *n_points = n_odb+1; |
807 |
|
809 | out = ray0; |
808 | return fate; |
810 | return fate; |
809 | } |
811 | } |
810 | void Guide::GetVFate(int *out) |
812 | void Guide::GetVFate(int *out) |
811 | { |
813 | { |
812 | for(int i=0;i<7;i++) out[i] = (int)hfate |
814 | for(int i=0;i<7;i++) out[i] = (int)hfate.GetBinContent(i+1); |
813 | } |
815 | } |
814 | void Guide::Draw(int color, int width) |
816 | void Guide::Draw(int color, int width) |
815 | { |
817 | { |
816 | for(int i = 0; i<6; i++) s_side |
818 | for(int i = 0; i<6; i++) (&s_side[i])->Draw(color, width); |
817 | } |
819 | } |
818 | void Guide::DrawSkel(int color, int width) |
820 | void Guide::DrawSkel(int color, int width) |
819 | { |
821 | { |
820 | TPolyLine3D *line3d = new TPolyLine3D(2); |
822 | TPolyLine3D *line3d = new TPolyLine3D(2); |
- | 823 | line3d->SetLineWidth(width); |
|
821 |
|
824 | line3d->SetLineColor(color); |
822 | 825 | ||
823 | for(int i=0; i<4; i++) { |
826 | for(int i=0; i<4; i++) { |
824 | line3d->SetPoint(0, vodnik_edge[i+0].x(), vodnik_edge[i+0].y(), vodnik_edge[i+0].z()); |
827 | line3d->SetPoint(0, vodnik_edge[i+0].x(), vodnik_edge[i+0].y(), vodnik_edge[i+0].z()); |
825 | line3d->SetPoint(1, vodnik_edge[i+4].x(), vodnik_edge[i+4].y(), vodnik_edge[i+4].z()); |
828 | line3d->SetPoint(1, vodnik_edge[i+4].x(), vodnik_edge[i+4].y(), vodnik_edge[i+4].z()); |
826 | line3d->DrawClone(); |
829 | line3d->DrawClone(); |
Line 897... | Line 900... | ||
897 | col_rgla(6), |
900 | col_rgla(6), |
898 | col_LG(1), |
901 | col_LG(1), |
899 | col_glass(4), |
902 | col_glass(4), |
900 | col_active(7), |
903 | col_active(7), |
901 | guide_on(parameters.getGuideOn()), |
904 | guide_on(parameters.getGuideOn()), |
902 | guide( |
905 | guide(Guide(center0, parameters)), |
903 | plate( |
906 | plate(Plate(parameters)), |
904 | _plateWidth(parameters.getPlateWidth()), |
907 | _plateWidth(parameters.getPlateWidth()), |
905 | _plateOn(parameters.getPlateOn()), |
908 | _plateOn(parameters.getPlateOn()), |
906 | offsetY(parameters.getOffsetY()), |
909 | offsetY(parameters.getOffsetY()), |
907 | offsetZ(parameters.getOffsetZ()) |
910 | offsetZ(parameters.getOffsetZ()) |
908 | { |
911 | { |
Line 927... | Line 930... | ||
927 | 930 | ||
928 | // additional glass between at top of SiPM |
931 | // additional glass between at top of SiPM |
929 | // example: epoxy n=1.60 |
932 | // example: epoxy n=1.60 |
930 | double n4 = 1.57; |
933 | double n4 = 1.57; |
931 | TVector3 plane_v[4]; |
934 | TVector3 plane_v[4]; |
932 | int nBins = nch |
935 | int nBins = nch; |
933 | double p_size = b/2.0; |
936 | double p_size = b/2.0; |
934 | plane_v[0].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap - p_size); |
937 | plane_v[0].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap - p_size); |
935 | plane_v[1].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap + p_size); |
938 | plane_v[1].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap + p_size); |
936 | plane_v[2].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap + p_size); |
939 | plane_v[2].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap + p_size); |
937 | plane_v[3].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap - p_size); |
940 | plane_v[3].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap - p_size); |
938 | glass |
941 | glass.Set(SURF_REFRA, plane_v, n3, n4, reflectivity); |
939 | glass |
942 | glass.FlipN(); |
940 | 943 | ||
941 | // additional circular glass between LG and SiPM |
944 | // additional circular glass between LG and SiPM |
942 | glass_circle |
945 | glass_circle.Set(TVector3(x_offset+d+glass_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), b); |
943 | 946 | ||
944 | hglass |
947 | //hglass = (TH2F*)gROOT->FindObject("hglass"); if(hglass) delete hglass; |
945 | hglass = |
948 | hglass = TH2F("hglass", "", |
946 | nBins, y_gap - p_size, y_gap + p_size, |
949 | nBins, y_gap - p_size, y_gap + p_size, |
947 | nBins, z_gap - p_size, z_gap + p_size); |
950 | nBins, z_gap - p_size, z_gap + p_size); |
948 | 951 | ||
949 | // SiPM active surface |
952 | // SiPM active surface |
950 | p_size = parameters.getActive()/2.0; |
953 | p_size = parameters.getActive()/2.0; |
Line 952... | Line 955... | ||
952 | 955 | ||
953 | plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size); |
956 | plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size); |
954 | plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size); |
957 | plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size); |
955 | plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size); |
958 | plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size); |
956 | plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size); |
959 | plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size); |
957 | active |
960 | active.Set(plane_v); |
958 | //active surface in case of bad coupling is circle d=a |
961 | //active surface in case of bad coupling is circle d=a |
959 | TVector3 activePosition(center); |
962 | TVector3 activePosition(center); |
960 | activePosition += TVector3(d + x_gap, 0, 0); |
963 | activePosition += TVector3(d + x_gap, 0, 0); |
961 | TVector3 normal(1,0,0); |
964 | TVector3 normal(1,0,0); |
962 | grease |
965 | grease.Set(activePosition, normal, 1.0*p_size); |
963 | 966 | ||
964 | hactive |
967 | //hactive = (TH2F*)gROOT->FindObject("hactive"); if(hactive) delete hactive; |
965 | //hactive = new TH2F("hactive", "Active area hits", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size); |
968 | //hactive = new TH2F("hactive", "Active area hits", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size); |
966 | hactive = |
969 | hactive = TH2F("hactive", ";x [mm];y [mm]", nBins, y_gap - p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ); |
967 | 970 | ||
968 | p_size = b/2.0; |
971 | p_size = b/2.0; |
969 | //p_size = 2.5; |
972 | //p_size = 2.5; |
970 | //p_size = M*0.6; |
973 | //p_size = M*0.6; |
971 | hlaser |
974 | //hlaser = (TH2F*)gROOT->FindObject("hlaser"); if(hlaser) delete hlaser; |
972 | hlaser = |
975 | hlaser = TH2F("hlaser", ";x [mm]; y [mm]", nBins, -p_size+offsetY, p_size+offsetY, nBins, -p_size+offsetZ, p_size+offsetZ); |
973 | 976 | ||
974 | // collection surface in SiPM plane |
977 | // collection surface in SiPM plane |
975 | p_size = 1.4*b/2.0; |
978 | p_size = 1.4*b/2.0; |
976 | plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size); |
979 | plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size); |
977 | plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size); |
980 | plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size); |
978 | plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size); |
981 | plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size); |
979 | plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size); |
982 | plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size); |
980 | detector |
983 | detector.Set(plane_v); |
981 | 984 | ||
982 | hdetector |
985 | //hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector; |
983 | //hdetector = new TH2F("hdetector", "Hits detector plane", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size); |
986 | //hdetector = new TH2F("hdetector", "Hits detector plane", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size); |
984 | hdetector = |
987 | hdetector = TH2F("hdetector", ";x [mm]; y [mm]", nBins, y_gap-p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ); |
985 | 988 | ||
986 | /* |
989 | /* |
987 | window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R); |
990 | window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R); |
988 | 991 | ||
989 | p_size = M*a; |
992 | p_size = M*a; |
Line 995... | Line 998... | ||
995 | 998 | ||
996 | hwindow = (TH2F*)gROOT->FindObject("hwindow"); if(hwindow) delete hwindow; |
999 | hwindow = (TH2F*)gROOT->FindObject("hwindow"); if(hwindow) delete hwindow; |
997 | hwindow = new TH2F("hwindow", "Hits Window", nch, y_gap - window_R, y_gap + window_R, nch, z_gap - window_R, z_gap + window_R); |
1000 | hwindow = new TH2F("hwindow", "Hits Window", nch, y_gap - window_R, y_gap + window_R, nch, z_gap - window_R, z_gap + window_R); |
998 | */ |
1001 | */ |
999 | p_size = b/2.0; |
1002 | p_size = b/2.0; |
1000 | histoPlate |
1003 | //histoPlate = (TH2F*)gROOT->FindObject("histoPlate"); if(histoPlate) delete histoPlate; |
1001 | histoPlate = |
1004 | histoPlate = TH2F("histoPlate", ";x [mm];y [mm]", nBins, -p_size, +p_size, nBins, -p_size, +p_size); |
- | 1005 | ||
- | 1006 | generated = TH2F("generated", ";x [mm];y [mm]", nBins, -p_size, +p_size, nBins, -p_size, +p_size); |
|
1002 | } |
1007 | } |
1003 | 1008 | ||
1004 | //----------------------------------------------------------------------------- |
1009 | //----------------------------------------------------------------------------- |
1005 | // vrne 1 ce je zadel aktvino povrsino |
1010 | // vrne 1 ce je zadel aktvino povrsino |
1006 | // vrne <1 ce jo zgresi |
1011 | // vrne <1 ce jo zgresi |
1007 | int CDetector::Propagate(CRay in, CRay |
1012 | int CDetector::Propagate(CRay in, CRay& out, int draw) |
1008 | // Sledi zarku skozi vodnik. Vrne: |
1013 | // Sledi zarku skozi vodnik. Vrne: |
1009 | // 0, ce zgresi vstopno ploskev MISSED |
1014 | // 0, ce zgresi vstopno ploskev MISSED |
1010 | // 1, ce zadane izstopno ploskev HIT |
1015 | // 1, ce zadane izstopno ploskev HIT |
1011 | // -1, ce se v vodniku ne odbije totalno REFRACTED |
1016 | // -1, ce se v vodniku ne odbije totalno REFRACTED |
1012 | // 2, enter the light guide, bin 2 of hfate EXIT |
1017 | // 2, enter the light guide, bin 2 of hfate EXIT |
1013 | // -2, ce se ne odbije zaradi koncnega R stranic - no total reflection REFRACTED |
1018 | // -2, ce se ne odbije zaradi koncnega R stranic - no total reflection REFRACTED |
1014 | // -3, ce se odbije nazaj in gre nazaj ven skozi sprednjo ploskev BACK_REFLECTED |
1019 | // -3, ce se odbije nazaj in gre nazaj ven skozi sprednjo ploskev BACK_REFLECTED |
1015 | // +4, ce se absorbira v materialu ABSORBED |
1020 | // +4, ce se absorbira v materialu ABSORBED |
1016 | { |
1021 | { |
1017 | if (dbg) printf("--- Detector::Propagate ---\n"); |
1022 | if (dbg) printf("--- Detector::Propagate ---\n"); |
- | 1023 | generated.Fill(in.GetR().y(), in.GetR().z()); |
|
1018 | //CRay *ray0 = new CRay; ray0->Set(in.GetR(), in.GetN()); ray0->SetColor(col_in); |
1024 | //CRay *ray0 = new CRay; ray0->Set(in.GetR(), in.GetN()); ray0->SetColor(col_in); |
1019 | CRay *rayin = new CRay(in); |
1025 | CRay *rayin = new CRay(in); |
1020 | rayin->SetColor(col_in); |
1026 | rayin->SetColor(col_in); |
1021 | CRay |
1027 | CRay rayout(in); |
1022 | rayout |
1028 | rayout.SetColor(col_in); |
1023 | 1029 | ||
1024 | const int max_n_points = guide |
1030 | const int max_n_points = guide.GetMAXODB() + 2; |
1025 | TVector3 pointsPlate[max_n_points]; |
1031 | TVector3 pointsPlate[max_n_points]; |
1026 | //TVector3 intersection; |
1032 | //TVector3 intersection; |
1027 | Fate fatePlate; |
1033 | Fate fatePlate; |
1028 | int nPointsPlate; |
1034 | int nPointsPlate; |
1029 | TPolyLine3D *line3d = new TPolyLine3D(2); |
1035 | TPolyLine3D *line3d = new TPolyLine3D(2); |
Line 1033... | Line 1039... | ||
1033 | // Draw the plate and propagate the ray through |
1039 | // Draw the plate and propagate the ray through |
1034 | // check if the ray should be reflected?? |
1040 | // check if the ray should be reflected?? |
1035 | 1041 | ||
1036 | if(_plateOn) { |
1042 | if(_plateOn) { |
1037 | 1043 | ||
1038 | fatePlate = plate |
1044 | fatePlate = plate.propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate); |
1039 | if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0); |
1045 | if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0); |
1040 | if(draw) { |
1046 | if(draw) { |
1041 | if(fatePlate == missed) { |
1047 | if(fatePlate == missed) { |
1042 | rayout |
1048 | rayout.SetColor(col_in); |
1043 | rayout |
1049 | rayout.DrawS(center.x() - _plateWidth, -10.0); |
1044 | } |
1050 | } |
1045 | else if(fatePlate == backreflected){ |
1051 | else if(fatePlate == backreflected){ |
1046 | rayout |
1052 | rayout.SetColor(kBlack); |
1047 | rayout |
1053 | rayout.DrawS(center.x() - _plateWidth, 7.0); |
1048 | if (dbg) printf("Backreflected at plate!\n"); |
1054 | if (dbg) printf("Backreflected at plate!\n"); |
1049 | } |
1055 | } |
1050 | else { |
1056 | else { |
1051 | int p_i; |
1057 | int p_i; |
1052 | for(p_i = 0; p_i < nPointsPlate-1; p_i++) { |
1058 | for(p_i = 0; p_i < nPointsPlate-1; p_i++) { |
1053 | line3d->SetPoint(0, pointsPlate[p_i].x(), pointsPlate[p_i].y(), pointsPlate[p_i].z()); |
1059 | line3d->SetPoint(0, pointsPlate[p_i].x(), pointsPlate[p_i].y(), pointsPlate[p_i].z()); |
1054 | line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z()); |
1060 | line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z()); |
1055 | line3d->DrawClone(); |
1061 | line3d->DrawClone(); |
1056 | } |
1062 | } |
1057 | rayout |
1063 | rayout.DrawS(pointsPlate[p_i].x(), -0.1); |
1058 | if(fatePlate == noreflection) { // lost on plate side |
1064 | if(fatePlate == noreflection) { // lost on plate side |
1059 | rayout |
1065 | rayout.SetColor(col_out); |
1060 | rayout |
1066 | rayout.DrawS(pointsPlate[p_i].x(), 10.0); |
1061 | } |
1067 | } |
1062 | } |
1068 | } |
1063 | } |
1069 | } |
1064 | 1070 | ||
1065 | if(! (fatePlate == hitExit or fatePlate == refracted) ) { |
1071 | if(! (fatePlate == hitExit or fatePlate == refracted) ) { |
1066 | guide |
1072 | guide.GetHFate()->Fill(rays); |
1067 | if (dbg)printf("CDetector::propagate Simulated ray missed the entry surface!\n"); |
1073 | if (dbg)printf("CDetector::propagate Simulated ray missed the entry surface!\n"); |
1068 | if (fatePlate == backreflected) |
1074 | if (fatePlate == backreflected) |
1069 | guide |
1075 | guide.GetHFate()->Fill(fatePlate); // reflected back |
1070 | else |
1076 | else |
1071 | guide |
1077 | guide.GetHFate()->Fill(noreflection); //lost on plate side |
1072 | return fatePlate; |
1078 | return fatePlate; |
1073 | } |
1079 | } |
1074 | 1080 | ||
1075 | //Ray hits light guide |
1081 | //Ray hits light guide |
1076 | histoPlate |
1082 | histoPlate.Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point |
1077 | 1083 | ||
1078 | } |
1084 | } |
1079 | else { |
1085 | else { |
1080 | //rayout = rayin; |
1086 | //rayout = rayin; |
1081 | if(draw) rayout |
1087 | if(draw) rayout.DrawS(center.x(), -10.0); |
1082 | } |
1088 | } |
1083 | 1089 | ||
1084 | // If the ray is not reflected in the plate |
1090 | // If the ray is not reflected in the plate |
1085 | // Draw the light guide and propagate the ray through |
1091 | // Draw the light guide and propagate the ray through |
1086 | 1092 | ||
Line 1088... | Line 1094... | ||
1088 | TVector3 points[max_n_points]; |
1094 | TVector3 points[max_n_points]; |
1089 | TVector3 presecisce; |
1095 | TVector3 presecisce; |
1090 | 1096 | ||
1091 | int n_points; |
1097 | int n_points; |
1092 | int fate_glass; |
1098 | int fate_glass; |
1093 | CRay *ray0 = new CRay( |
1099 | CRay *ray0 = new CRay(rayout); |
1094 | // delete rayout; -> creates dangling reference when tries to delete ray0! |
1100 | // delete rayout; -> creates dangling reference when tries to delete ray0! |
1095 | //delete rayin; -> delete rayout! |
1101 | //delete rayin; -> delete rayout! |
1096 | CRay |
1102 | CRay ray1; |
1097 | 1103 | ||
1098 | fate = guide |
1104 | fate = guide.PropagateRay(*ray0, ray1, &n_points, points); |
1099 | if (dbg) { |
1105 | if (dbg) { |
1100 | if (fate == backreflected) printf("DETECTOR::backreflected\n"); |
1106 | if (fate == backreflected) printf("DETECTOR::backreflected\n"); |
1101 | } |
1107 | } |
1102 | 1108 | ||
1103 | line3d->SetLineColor(col_lg); |
1109 | line3d->SetLineColor(col_lg); |
1104 | int p_i; |
1110 | int p_i; |
1105 | if(guide_on) { |
1111 | if(guide_on) { |
1106 | if(draw) { |
1112 | if(draw) { |
1107 | if(fate == missed) { |
1113 | if(fate == missed) { |
1108 | if (dbg) printf("Detector: fate=missed\n"); |
1114 | if (dbg) printf("Detector: fate=missed\n"); |
1109 | TVector3 r = ray1 |
1115 | TVector3 r = ray1.GetR(); |
1110 | TVector3 k = ray1 |
1116 | TVector3 k = ray1.GetK(); |
1111 | ray1 |
1117 | ray1.Set(r,k); |
1112 | ray1 |
1118 | ray1.DrawS(center.x(), 10.0); |
1113 | } else { |
1119 | } else { |
1114 | for(p_i = 0; p_i < n_points-1; p_i++) { |
1120 | for(p_i = 0; p_i < n_points-1; p_i++) { |
1115 | line3d->SetPoint(0, points[p_i].x(), points[p_i].y(), points[p_i].z()); |
1121 | line3d->SetPoint(0, points[p_i].x(), points[p_i].y(), points[p_i].z()); |
1116 | line3d->SetPoint(1, points[p_i+1].x(), points[p_i+1].y(), points[p_i+1].z()); |
1122 | line3d->SetPoint(1, points[p_i+1].x(), points[p_i+1].y(), points[p_i+1].z()); |
1117 | line3d->DrawClone(); |
1123 | line3d->DrawClone(); |
1118 | } |
1124 | } |
1119 | if(fate != noreflection) { |
1125 | if(fate != noreflection) { |
1120 | if (dbg) printf("Detector: fate != noreflection, fate = %d\n", (int)fate); |
1126 | if (dbg) printf("Detector: fate != noreflection, fate = %d\n", (int)fate); |
1121 | if(glass_on) {/*if(fate == 1)*/ ray1 |
1127 | if(glass_on) {/*if(fate == 1)*/ ray1.Draw(points[p_i].x(), center.x() + guide.getD() + glass_d);} |
1122 | else { |
1128 | else { |
1123 | ray1 |
1129 | ray1.SetColor(col_out); |
1124 | ray1 |
1130 | ray1.DrawS(points[p_i].x(), 10.0); |
1125 | } |
1131 | } |
1126 | } |
1132 | } |
1127 | } |
1133 | } |
1128 | } |
1134 | } |
1129 | 1135 | ||
1130 | 1136 | ||
1131 | if(! (fate == hitExit or fate == refracted) ) { |
1137 | if(! (fate == hitExit or fate == refracted) ) { |
1132 | if (dbg) printf("Detector: fate != hit, refracted\n"); |
1138 | if (dbg) printf("Detector: fate != hit, refracted\n"); |
1133 |
|
1139 | out = ray1; |
1134 | delete ray0; |
1140 | delete ray0; |
1135 |
|
1141 | //delete ray1; |
1136 |
|
1142 | //delete rayout; |
1137 | delete rayin; |
1143 | delete rayin; |
1138 | return fate; |
1144 | return fate; |
1139 | } |
1145 | } |
1140 | } else { |
1146 | } else { |
1141 | if (dbg) printf("Detector: fate = hit or refracted"); |
1147 | if (dbg) printf("Detector: fate = hit or refracted"); |
1142 | ray1 = ray0; |
1148 | ray1 = *ray0; |
1143 | if(draw) { |
1149 | if(draw) { |
1144 | //double epoxy = parameters->getGlassD(); |
1150 | //double epoxy = parameters->getGlassD(); |
1145 | if(glass_on) ray1 |
1151 | if(glass_on) ray1.Draw(center.x(), center.x() + glass_d); |
1146 | else ray1 |
1152 | else ray1.DrawS(center.x(), 10.0); |
1147 | } |
1153 | } |
1148 | } |
1154 | } |
1149 | 1155 | ||
1150 | fate = missed; // zgresil aktivno povrsino |
1156 | fate = missed; // zgresil aktivno povrsino |
1151 | if(glass_on) { |
1157 | if(glass_on) { |
1152 | *ray0 = |
1158 | *ray0 = ray1; |
1153 | ray1 |
1159 | ray1.SetColor(col_rgla); |
1154 | fate_glass = glass |
1160 | fate_glass = glass.PropagateRay(*ray0, ray1, &presecisce); |
1155 | if(fate_glass == REFRACTION) { |
1161 | if(fate_glass == REFRACTION) { |
1156 | hglass |
1162 | hglass.Fill(presecisce.y(), presecisce.z()); |
1157 | if(draw) ray1 |
1163 | if(draw) ray1.DrawS(presecisce.x(), 10.0); |
1158 | //if(active->TestIntersection(&presecisce, *ray1)) { |
1164 | //if(active->TestIntersection(&presecisce, *ray1)) { |
1159 | //fate = hitExit; |
1165 | //fate = hitExit; |
1160 | // |
1166 | //h-Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
1161 | //hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ); |
1167 | //hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ); |
1162 | //} |
1168 | //} |
1163 | //if(detector->TestIntersection(&presecisce, *ray1)) |
1169 | //if(detector->TestIntersection(&presecisce, *ray1)) |
1164 | //hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
1170 | //hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
1165 | //} else if(fate_glass == REFLECTION) { |
1171 | //} else if(fate_glass == REFLECTION) { |
1166 | else |
1172 | else |
1167 | if(draw) ray1 |
1173 | if(draw) ray1.DrawS(presecisce.x(), 10.0); |
1168 | } |
1174 | } |
1169 | } |
1175 | } |
1170 | 1176 | ||
1171 | // Main test: ray and SiPM surface |
1177 | // Main test: ray and SiPM surface |
1172 | if(active |
1178 | if(active.TestIntersection(&presecisce, ray1)) { |
1173 | fate = hitExit; |
1179 | fate = hitExit; |
1174 | hactive |
1180 | hactive.Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
1175 | hlaser |
1181 | hlaser.Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ); |
1176 | } |
1182 | } |
1177 | // If it is on the same plane as SiPM |
1183 | // If it is on the same plane as SiPM |
1178 | if(detector |
1184 | if(detector.TestIntersection(&presecisce, ray1)) |
1179 | hdetector |
1185 | hdetector.Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
1180 | //} |
1186 | //} |
1181 | //} else { |
1187 | //} else { |
1182 | //if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d); |
1188 | //if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d); |
1183 | //} |
1189 | //} |
1184 | 1190 | ||
1185 |
|
1191 | out = ray1; |
1186 | delete ray0; |
1192 | delete ray0; |
1187 |
|
1193 | //delete ray1; |
1188 |
|
1194 | //delete rayout; |
1189 | delete rayin; |
1195 | delete rayin; |
1190 | return fate; |
1196 | return fate; |
1191 | } |
1197 | } |
1192 | 1198 | ||
1193 | void CDetector::Draw(int width) |
1199 | void CDetector::Draw(int width) |
1194 | { |
1200 | { |
1195 | if(guide_on) { |
1201 | if(guide_on) { |
1196 | if( TMath::Abs(guide |
1202 | if( TMath::Abs(guide.getN1()-guide.getN2()) < MARGIN ) { |
1197 | if(_plateOn) plate |
1203 | if(_plateOn) plate.drawSkel(col_LG, width); |
1198 | guide |
1204 | guide.DrawSkel(col_LG, width); |
1199 | } |
1205 | } |
1200 | else { |
1206 | else { |
1201 | if(_plateOn) plate |
1207 | if(_plateOn) plate.draw(4, width); |
1202 | guide |
1208 | guide.Draw(col_LG, width); |
1203 | } |
1209 | } |
1204 | } |
1210 | } |
1205 | 1211 | ||
1206 | if(glass_on) glass_circle |
1212 | if(glass_on) glass_circle.Draw(col_glass, width); |
1207 | //window_circle->Draw(col_glass, width); |
1213 | //window_circle->Draw(col_glass, width); |
1208 | active |
1214 | active.Draw(col_active, width); |
1209 | } |
1215 | } |
1210 | 1216 | ||
1211 | 1217 | ||
1212 | Plate::Plate(DetectorParameters& parameters) |
1218 | Plate::Plate(DetectorParameters& parameters) |
1213 | { |
1219 | { |
Line 1228... | Line 1234... | ||
1228 | plate_edge[6].SetXYZ(plateWidth,-t, t); |
1234 | plate_edge[6].SetXYZ(plateWidth,-t, t); |
1229 | plate_edge[7].SetXYZ(plateWidth,-t,-t); |
1235 | plate_edge[7].SetXYZ(plateWidth,-t,-t); |
1230 | 1236 | ||
1231 | for(int i = 0; i<8; i++) plate_edge[i] += center; |
1237 | for(int i = 0; i<8; i++) plate_edge[i] += center; |
1232 | 1238 | ||
1233 | sides[0] |
1239 | sides[0].Set(SURF_REFRA, plate_edge, n1, n2, c_reflectivity); |
1234 | sides[0] |
1240 | sides[0].FlipN(); |
1235 | 1241 | ||
1236 | sides[1] |
1242 | sides[1].Set(SURF_REFRA, plate_edge[3], plate_edge[2], plate_edge[6], plate_edge[7], n2, n2, c_reflectivity); |
1237 | sides[2] |
1243 | sides[2].Set(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, c_reflectivity); |
1238 | sides[3] |
1244 | sides[3].Set(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, c_reflectivity); |
1239 | sides[4] |
1245 | sides[4].Set(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, c_reflectivity); |
1240 | 1246 | ||
1241 | sides[5] |
1247 | sides[5].Set(SURF_REFRA, &plate_edge[4], n2, n2, c_reflectivity); |
1242 | sides[5] |
1248 | sides[5].FlipN(); |
1243 | 1249 | ||
1244 | for(int i=0; i<6; i++) sides[i] |
1250 | for(int i=0; i<6; i++) sides[i].SetFresnel(1); |
1245 | } |
1251 | } |
1246 | 1252 | ||
1247 | void Plate::draw(int color, int width) |
1253 | void Plate::draw(int color, int width) |
1248 | { |
1254 | { |
1249 | for(int i = 0; i<6; i++) sides[i] |
1255 | for(int i = 0; i<6; i++) sides[i].Draw(color, width); |
1250 | } |
1256 | } |
1251 | 1257 | ||
1252 | void Plate::drawSkel(int color, int width) |
1258 | void Plate::drawSkel(int color, int width) |
1253 | { |
1259 | { |
1254 | TPolyLine3D line3d(2); |
1260 | TPolyLine3D line3d(2); |
Line 1260... | Line 1266... | ||
1260 | line3d.SetPoint(1, plate_edge[i+4].x(), plate_edge[i+4].y(), plate_edge[i+4].z()); |
1266 | line3d.SetPoint(1, plate_edge[i+4].x(), plate_edge[i+4].y(), plate_edge[i+4].z()); |
1261 | line3d.DrawClone(); |
1267 | line3d.DrawClone(); |
1262 | } |
1268 | } |
1263 | } |
1269 | } |
1264 | 1270 | ||
1265 | Fate Plate::propagateRay(CRay in, CRay |
1271 | Fate Plate::propagateRay(CRay in, CRay& out, int *n_points, TVector3 *points) |
1266 | { |
1272 | { |
1267 | CRay ray0; |
1273 | CRay ray0; |
1268 | CRay ray1; |
1274 | CRay ray1; |
1269 | TVector3 vec0, vec1; |
1275 | TVector3 vec0, vec1; |
1270 | Fate fate = enter; |
1276 | Fate fate = enter; |
Line 1273... | Line 1279... | ||
1273 | ray0 = in; |
1279 | ray0 = in; |
1274 | int n_odb = 0; |
1280 | int n_odb = 0; |
1275 | int last_hit = 0; |
1281 | int last_hit = 0; |
1276 | int propagation = 0; |
1282 | int propagation = 0; |
1277 | 1283 | ||
1278 | int result = sides[0] |
1284 | int result = sides[0].PropagateRay(ray0, ray1, &vec1); |
1279 | if( !result ) { |
1285 | if( !result ) { |
1280 | // ce -NI- presecisca z vstopno |
1286 | // ce -NI- presecisca z vstopno |
1281 | fate = missed; |
1287 | fate = missed; |
1282 | } else if(result == REFLECTION) { |
1288 | } else if(result == REFLECTION) { |
1283 | if (dbg) printf("PLATE: reflected\n"); |
1289 | if (dbg) printf("PLATE: reflected\n"); |
Line 1291... | Line 1297... | ||
1291 | ray0 = ray1; |
1297 | ray0 = ray1; |
1292 | vec0 = vec1; |
1298 | vec0 = vec1; |
1293 | propagation = 11; |
1299 | propagation = 11; |
1294 | for(inters_i=0; inters_i<6; inters_i++) { |
1300 | for(inters_i=0; inters_i<6; inters_i++) { |
1295 | if( inters_i != last_hit) { |
1301 | if( inters_i != last_hit) { |
1296 | if( sides[inters_i] |
1302 | if( sides[inters_i].TestIntersection(&vec1, ray1) ) break; |
1297 | } |
1303 | } |
1298 | } |
1304 | } |
1299 | points[n_odb] = vec1; |
1305 | points[n_odb] = vec1; |
1300 | if(inters_i == 0) { |
1306 | if(inters_i == 0) { |
1301 | ray0 = ray1; |
1307 | ray0 = ray1; |
1302 | fate = backreflected; |
1308 | fate = backreflected; |
1303 | break;} // backreflection |
1309 | break;} // backreflection |
1304 | 1310 | ||
1305 | propagation = sides[inters_i] |
1311 | propagation = sides[inters_i].PropagateRay(ray0, ray1, &vec1); |
1306 | if(inters_i == 5) { // successfull exit |
1312 | if(inters_i == 5) { // successfull exit |
1307 | fate = hitExit; |
1313 | fate = hitExit; |
1308 | //hout->Fill(vec1.y(), vec1.z()); |
1314 | //hout->Fill(vec1.y(), vec1.z()); |
1309 | //hnodb_exit->Fill(n_odb-1); |
1315 | //hnodb_exit->Fill(n_odb-1); |
1310 | n_odb++; |
1316 | n_odb++; |
Line 1327... | Line 1333... | ||
1327 | last_hit = inters_i; |
1333 | last_hit = inters_i; |
1328 | } |
1334 | } |
1329 | } |
1335 | } |
1330 | 1336 | ||
1331 | *n_points = n_odb+1; |
1337 | *n_points = n_odb+1; |
1332 |
|
1338 | out = ray0; |
1333 | return fate; |
1339 | return fate; |
1334 | }; |
1340 | }; |
1335 | 1341 | ||
1336 | 1342 |