Subversion Repositories f9daq

Rev

Rev 70 | Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
25 f9daq 1
#include "include/guide.h"
2
 
3
#include <iostream>
4
 
5
 
6
// vector output shortcut
7
void printv(TVector3 v)
8
{
9
        printf("(x,y,z) = (%.4lf, %.4lf, %.4lf)\n", v.x(), v.y(), v.z());
10
}
11
// TVector3::Rotate does not seem accurate enough
12
TVector3 rotatey(TVector3 v, double theta)
13
{
14
        return TVector3(v.x() * TMath::Cos(theta) + v.z() * TMath::Sin(theta),
15
                                        v.y(),
16
                                        -v.x() * TMath::Sin(theta) + v.z() * TMath::Cos(theta));
17
}
18
// another shortcut not found in TMath
19
int sign(double in)
20
{
21
        if(in >= 0.0) return 1;
22
        else return -1;
23
}
24
//=================================================================================
25
 
26
//-----------------------------------------------------------------------------
27
void CRay::Set(TVector3 r0, TVector3 n0)
28
{
29
        r = r0; n = n0.Unit();
30
}
31
//-----------------------------------------------------------------------------
32
//void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0)
33
//{
34
        //r.SetXYZ(x0, y0, z0);
35
        //n.SetXYZ(l0, m0, n0); n = n.Unit();
36
//}
37
//-----------------------------------------------------------------------------
38
/*
39
CRay& CRay::operator = (const CRay& p)
40
{
41
        r.SetXYZ(p.GetR().x(), p.GetR().y(), p.GetR().z());
42
        //this->r.SetXYZ(p.x(), p.y(), p.z());
43
        n.SetXYZ(p.GetN().x(), p.GetN().y(), p.GetN().z());
44
        return *this;
45
} */
46
//-----------------------------------------------------------------------------
47
void CRay::Print()
48
{
49
        printf("---> CRay::Print() <---\n");
50
        printf("(x,y,z)=(%.2lf, %.2lf, %.2lf); (l,m,n)=(%.2lf, %.2lf, %.2lf)\n",
51
                r.x(), r.y(), r.z(), n.x(), n.y(), n.z());
52
}
53
//-----------------------------------------------------------------------------
54
void CRay::Draw()
55
{
56
double t = 50.0;
57
TPolyLine3D *line3d = new TPolyLine3D(2);
58
        //line3d->SetPoint(0, r.x() - t*n.x(), r.y() - t*n.y(), r.z() - t*n.z());
59
        line3d->SetPoint(0, r.x(), r.y(), r.z());
60
        line3d->SetPoint(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.z());
61
        line3d->SetLineWidth(1);
62
        line3d->SetLineColor(color);
63
 
64
        line3d->Draw();
65
}
66
//-----------------------------------------------------------------------------
67
void CRay::Draw(double x_from, double x_to)
68
{
69
double A1, A2;
70
  TPolyLine3D *line3d = new TPolyLine3D(2);
71
 
72
        if(n.x() < MARGIN) {
73
                A1 = A2 = 0.0;
74
        } else {
75
                A1 = (x_from - r.x())/n.x();   
76
                A2 = (x_to - r.x())/n.x();
77
        }
78
 
79
        line3d->SetPoint(0, x_from, A1*n.y()+r.y(), A1*n.z()+r.z());
80
        line3d->SetPoint(1, x_to, A2*n.y()+r.y(), A2*n.z()+r.z());
81
        line3d->SetLineWidth(1);
82
        line3d->SetLineColor(color);
83
 
84
        line3d->Draw();
85
}
86
//-----------------------------------------------------------------------------
87
void CRay::DrawS(double x_from, double t)
88
{
89
        double A1;
90
        TPolyLine3D *line3d = new TPolyLine3D(2);
91
 
92
        if(n.x() < MARGIN)
93
                A1 = 0.0;
94
        else
95
                A1 = (x_from - r.x())/n.x();
96
 
97
        line3d->SetPoint(0, x_from, A1*n.y()+r.y(), A1*n.z()+r.z());
98
        line3d->SetPoint(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.z());
99
        line3d->SetLineWidth(1);
100
        line3d->SetLineColor(color);
101
 
102
        line3d->Draw();
103
}
104
//=================================================================================
105
 
106
 
107
//=================================================================================
108
CPlane4::CPlane4() :
109
        n(TVector3(1.0, 0.0, 0.0)),
110
  A(0),
111
  B(0),
112
  C(0),
113
  D(0)
114
{ r[0] = TVector3(0.0,-1.0,-1.0);
115
        r[1] = TVector3(0.0,-1.0, 1.0);
116
        r[2] = TVector3(0.0, 1.0, 1.0);
117
        r[3] = TVector3(0.0, 1.0,-1.0);
118
        for(int i=0;i<4;i++) edge[i] = TVector3(0,0,0);
119
        for(int i=0;i<4;i++) angle_r[i] = 0;
120
};
121
//-----------------------------------------------------------------------------
122
CPlane4::CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
123
{
124
        //Set(r1, r2, r3, r4);
125
//}
126
//-----------------------------------------------------------------------------
127
// za izracun parametrov ravnine je en vektor prevec, vendar tega 
128
// rabim kot zadnji vogal poligona
129
//void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
130
//{
131
double x1,y1,z1, x2,y2,z2, x3,y3,z3;
132
 
133
        x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
134
        x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
135
        x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
136
 
137
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
138
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
139
        C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
140
        D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
141
 
142
        r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
143
        n.SetXYZ(A, B, C);
144
        n = n.Unit();
145
 
146
        for(int i=0;i<4;i++)
147
                edge[i] = r[i-3 ? i+1 : 0] - r[i];
148
 
149
        for(int i=0;i<4;i++)
150
                angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
151
};
152
 
153
void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
154
{
155
  double x1,y1,z1, x2,y2,z2, x3,y3,z3;
156
 
157
        x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
158
        x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
159
        x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
160
 
161
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
162
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
163
        C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
164
        D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
165
 
166
        r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
167
        n.SetXYZ(A, B, C);
168
        n = n.Unit();
169
 
170
        for(int i=0;i<4;i++)
171
                edge[i] = r[i-3 ? i+1 : 0] - r[i];
172
 
173
        for(int i=0;i<4;i++)
174
                angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
175
};
176
 
177
CPlane4::CPlane4(TVector3 *vr)
178
{
179
  double x1,y1,z1, x2,y2,z2, x3,y3,z3;
180
 
181
        x1 = vr[0].x(); y1 = vr[0].y(); z1 = vr[0].z();
182
        x2 = vr[1].x(); y2 = vr[1].y(); z2 = vr[1].z();
183
        x3 = vr[2].x(); y3 = vr[2].y(); z3 = vr[2].z();
184
 
185
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
186
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
187
        C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
188
        D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
189
 
190
        r[0] = vr[0]; r[1] = vr[1]; r[2] = vr[2]; r[3] = vr[3];
191
        n.SetXYZ(A, B, C);
192
        n = n.Unit();
193
 
194
        for(int i=0;i<4;i++)
195
                edge[i] = r[i-3 ? i+1 : 0] - r[i];
196
 
197
        for(int i=0;i<4;i++)
198
                angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
199
};
200
//-----------------------------------------------------------------------------
201
// posce presecisce !neskoncne! ravnine s premico (class CRay)
202
// ce najde presecisce vrne 1
203
int CPlane4::GetIntersection(TVector3 *vec, CRay ray)
204
{
205
TVector3 N; //nenormirani vektor (A,B,C)
206
double num, den; //stevec, imenovalec
207
double t;
208
TVector3 tmp;
209
 
210
        N.SetXYZ(A,B,C);
211
 
212
        num = N*ray.GetR() + D;
213
        den = N*ray.GetN();
214
 
215
        if (dbg) printf("t = %6.3lf / %6.3lf =  %6.3lf\n", num, den, num/den);
216
 
217
        //if(den == 0)
218
        if(TMath::Abs(den) < MARGIN) {
219
                //if(num == 0)
220
                if(TMath::Abs(num) < MARGIN) {
221
                  if (dbg) printf("The ray is on the surface!\n");
222
                        return 0; //return 2; // premica lezi na ravnini
223
                }
224
                else {
225
                        if (dbg) printf("The ray is parallel to the surface!\n");
226
                        return 0; // ni presecisca
227
                }
228
        }
229
 
230
        t = num / den;
231
 
232
        tmp = ray.GetR();
233
        tmp -= t*ray.GetN();
234
        *vec = tmp;
235
        return 1;
236
}
237
//-----------------------------------------------------------------------------
238
// ali je vektor vec, ki lezi na ravnini skupaj z e1 in e2, med njima
239
// angle_r je kot med e1 in e2, vsi vektorji imajo skupno izhodisce
240
int CPlane4::IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double angle)
241
{
242
  double angle_ve1, angle_ve2;
243
 
244
        if(dbg) printf("--- CPlane4::IsInTri ---\n");
245
 
246
        angle_ve1 = TMath::ACos(/*TMath::Abs*/( (e1.Unit()) * (vec.Unit()) ));
247
        angle_ve2 = TMath::ACos(/*TMath::Abs*/( (e2.Unit()) * (vec.Unit()) ));
248
 
249
        if(dbg)
250
        {
251
                //printf("angle_ve1 = %lf\n", angle_ve1*DEGREE);
252
                //printf("angle_ve2 = %lf\n", angle_ve2*DEGREE);
253
                //printf("angle_sum = %lf\n", (angle_ve1 + angle_ve2)*DEGREE);
254
                printf("  angle_r   = %lf\n", angle*DEGREE);
255
        }
256
 
257
  bool difference = (MARGIN < TMath::Abs(angle - (angle_ve1 + angle_ve2)));
258
  if (dbg) printf("  MARGIN < Difference = %d\n", difference);
259
  return (int) !difference;
260
}
261
//-----------------------------------------------------------------------------
262
// ali je vektor vec, ki lezi na ravnini!, znotraj meja, ki jih definirajo
263
// strije vogali te ravnine r[i]
264
int CPlane4::IsVectorIn(TVector3 vec)
265
{
266
  int status;
267
 
268
        if(dbg) printf("--- CPlane4::IsVectorIn ---\n");
269
 
270
        for(int i=0;i<3;i++)
271
        {
272
                status = IsInTri(vec - r[i], edge[i], -edge[i ? i-1 : 3], angle_r[i]);
273
                if(dbg) printf("  [%d] vec is %s\n", i, status ? "inside" : "outside");
274
                if(!status) return 0;
275
        }
276
 
277
        return 1;
278
}
279
//-----------------------------------------------------------------------------
280
int CPlane4::TestIntersection(CRay in)
281
{
282
        TVector3 tmp;
283
 
284
        if( GetIntersection(&tmp, in) )
285
                if( IsVectorIn(tmp) )
286
                        return 1;
287
 
288
        return 0;
289
}
290
//-----------------------------------------------------------------------------
291
int CPlane4::TestIntersection(TVector3 *vec, CRay in)
292
{
293
        TVector3 tmp;
294
 
295
        if( GetIntersection(&tmp, in) )
296
                if( IsVectorIn(tmp) ) {
297
                        *vec = tmp;
298
                        return 1;
299
                }
300
 
301
        return 0;
302
}
303
//-----------------------------------------------------------------------------
304
void CPlane4::Print()
305
{
306
        printf("--- CPlane4::Print() ---\n");
307
        printf("  r=(%.2lf, %.2lf, %.2lf); n=(%.2lf, %.2lf, %.2lf); ",
308
                r[0].x(), r[0].y(), r[0].z(), n.x(), n.y(), n.z());
309
        printf(  "(A,B,C,D)=(%.2lf, %.2lf, %.2lf, %.2lf) \n", A, B, C, D);
310
        for(int i=0;i<4;i++) printf("  edge[%d] = (%lf, %lf, %lf)\n", i, edge[i].x(), edge[i].y(), edge[i].z());
311
        for(int i=0;i<4;i++) printf("  angle[%d] = %lf\n", i, angle_r[i]*DEGREE);
312
}
313
//-----------------------------------------------------------------------------
314
void CPlane4::Draw(int color, int width)
315
{
316
TPolyLine3D *line3d = new TPolyLine3D(5);
317
 
318
        for(int i=0;i<4;i++) line3d->SetPoint(i, r[i].x(), r[i].y(), r[i].z());
319
        line3d->SetPoint(4, r[0].x(), r[0].y(), r[0].z());
320
        line3d->SetLineWidth(width); line3d->SetLineColor(color);
321
 
322
        line3d->Draw();
323
}
324
//=================================================================================
325
 
326
 
327
//=================================================================================
328
CSurface::CSurface(int type0):
329
  type(type0)
330
{
331
TVector3 vr[4];
332
TDatime now;
333
 
334
        vr[0].SetXYZ(0.0,-1.0,-1.0);
335
        vr[1].SetXYZ(0.0,-1.0, 1.0);
336
        vr[2].SetXYZ(0.0, 1.0, 1.0);
337
        vr[3].SetXYZ(0.0, 1.0,-1.0);
338
        //CPlane4::Set(vr);
339
        SetIndex(1.0, 1.5);
340
 
341
        reflection = c_reflectivity;   
342
        rand.SetSeed(now.Get());
343
 
344
        SetFresnel();
345
}
346
//-----------------------------------------------------------------------------
347
CSurface::CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
348
{
349
TDatime now;
350
 
351
        type = type0; CPlane4::Set(r1, r2, r3, r4);
352
        SetIndex(n10, n20);
353
 
354
        reflection = reflectivity;
355
        rand.SetSeed(now.Get());
356
 
357
        SetFresnel();
358
}
359
//-----------------------------------------------------------------------------
360
CSurface::CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
361
{
362
TDatime now;
363
 
364
        type = type0; CPlane4::Set(vr);
365
        SetIndex(n10, n20);
366
 
367
        reflection = reflectivity;
368
        rand.SetSeed(now.Get());
369
 
370
        SetFresnel();
371
}
372
//-----------------------------------------------------------------------------
373
void CSurface::SetIndex(double n10, double n20)
374
{
375
        n1 = n10; n2 = n20; n1_n2 = n1/n2;
376
 
377
        if(n1 > n2)
378
                cosTtotal = TMath::Sqrt( 1 - TMath::Power(n2/n1, 2) );
379
        else
380
                cosTtotal = 0.0;
381
}
382
//-----------------------------------------------------------------------------
383
// sprejme zarek, vrne uklonjen/odbit zarek in presecisce
384
// vrne 0 ce ni presecisca; 1 ce se je lomil
385
// 2 ce se je odbil; -2 ce se je absorbiral
386
int CSurface::PropagateRay(CRay in, CRay *out, TVector3 *intersection)
387
{
388
  if (dbg) printf("--- CSurface::PropagateRay ---\n");
389
  double cosTi, cosTt, p_ref;
390
  TVector3 intersect, transmit;
391
 
392
        if( !(GetIntersection(&intersect, in) == 1) )
393
                return 0;
394
 
395
        *intersection = intersect;
396
        if( !IsVectorIn(intersect) )
397
                return 0;
398
 
399
        // --------------- Fresnel ----------------------------------------------------
400
        // R_f = a_te * R_te  +  a_tm * R_tm
401
        double r_te=0;
402
        double r_tm=0;
403
        double R_te=0;
404
        double R_tm=0;
405
        double R_f = 0.0;
406
        double a_te = 0.0;
407
        double a_tm = 0.0;
408
        TVector3 v_te; // polarization perpendicular to the plane of incidence
409
        TVector3 v_tm; // inbound polarization parallel with the plane of incidence
410
        TVector3 v_tm_t;// transmited polarization parallel with the plane of incidence
411
        TVector3 pol_t = in.GetP(); // transmited polarization
412
        int sign_n; // sign of normal direction vs. inbound ray
413
        double cosTN; // debug
414
 
415
        if(fresnel) {  
416
                v_te = n.Cross(in.GetN()); v_te = v_te.Unit();
417
                v_tm = -v_te.Cross(in.GetN()); v_tm = v_tm.Unit();
418
                if(dbg) {
419
                        printf("  v_te = "); printv(v_te);
420
                        printf("  v_tm = "); printv(v_tm);
421
                }
422
 
423
                double cosAf = v_te * in.GetP();
424
                if(dbg) printf("  cosAf = %lf (Af = %lf)\n", cosAf, TMath::ACos(cosAf)*DEGREE);
425
 
426
                a_te = cosAf;
427
                a_tm = TMath::Sqrt(1 - cosAf*cosAf);
428
                if(dbg) printf("  a_te = %lf, a_tm = %lf\n", a_te, a_tm);
429
        }
430
        // ----------------------------------------------------------------------------
431
 
432
        if(type == 3) type = SURF_REFRA; //SURF_TOTAL -> SURF_REFRA
433
        switch(type){
434
                // ----------------------------------------------------------------------------
435
                // --------------- refraction from n1 to n2 -----------------------------------
436
                // ----------------------------------------------------------------------------
437
                case SURF_REFRA:
438
                        cosTi = in.GetN() * n;
439
                        if(dbg) printf("  cosTi = %lf (Ti = %lf)\n", cosTi, TMath::ACos(cosTi)*DEGREE);
440
                        sign_n = -sign(cosTi);
441
                        if(dbg) printf("  sign_n = %d\n", sign_n);
442
                        cosTi = TMath::Abs(cosTi);
443
 
444
                        // Check if there can be total reflection: n1 > n2
445
                        if(N1_N2(-sign_n) < 1.0)
446
                                cosTtotal = TMath::Sqrt( 1 - TMath::Power(N1_N2(-sign_n), 2) );
447
                        else
448
                                cosTtotal = 0.0;
449
 
450
                        if(dbg) printf("  cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
451
 
452
                        p_ref = rand.Uniform(0.0, 1.0);
453
 
454
                        if( (cosTi <= cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
455
                                if(dbg) printf("  TOTAL\n");
456
                                transmit = in.GetN() + sign_n*2*cosTi*n;
457
 
458
                                if(dbg) {
459
                                        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);
461
                                }      
462
                                out->Set(intersect, transmit);
463
 
464
                                pol_t = -in.GetP() + sign_n*2*cosTi*n;
465
                                out->SetPolarization(pol_t);
466
                                return 2;
467
                        } else { // ni tot. odboja
468
                                if(dbg) printf("  REFRACTION\n");
469
                                if(dbg) printf("  N1_N2(sign_n) = %lf\n", N1_N2(sign_n));      
470
                                cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));                          
471
                                if(dbg) printf("  cosTt = %lf (Tt = %lf) \n", cosTt, TMath::ACos(cosTt)*DEGREE);
472
 
473
                                transmit = N1_N2(sign_n)*in.GetN() + sign_n*(N1_N2(sign_n)*cosTi - cosTt)*n;
474
                                if(dbg) {printf("  transmit.Unit() = "); printv(transmit.Unit());}
475
 
476
                                if(dbg) {
477
                                        cosTN = transmit.Unit() * n;
478
                                        printf("  cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE); 
479
                                }
480
                                //if(cosTi<=cosTtotal) cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));
481
                                if(fresnel) {
482
                                        r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse
483
                                        r_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel
484
 
485
                                        if(dbg) printf("  r_te = %lf, r_tm = %lf\n", r_te, r_tm);
486
 
487
                                        v_tm_t = -v_te.Cross(transmit);
488
                                        v_tm_t = v_tm_t.Unit();
489
                                        pol_t = a_te * (1.0 -  TMath::Abs(r_te)) * v_te  +  a_tm * (1.0 -  TMath::Abs(r_tm)) * v_tm_t;
490
 
491
                                        if(dbg) {
492
                                                printf("  v_tm_t = "); printv(v_tm_t);
493
                                                printf("  pol_t = "); printv(pol_t);
494
                                        }
495
 
496
                                        R_te = TMath::Power(r_te, 2);
497
                                        R_tm = TMath::Power(r_tm, 2);
498
                                        R_f = a_te*a_te*R_te + a_tm*a_tm*R_tm;
499
 
500
                                        if(dbg) printf("  R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
501
                                }
502
 
503
                                //p_ref = rand.Uniform(0.0, 1.0);       
504
                                if(p_ref >= R_f) { // se lomi
505
                                        out->Set(intersect, transmit);
506
                                        out->SetPolarization(pol_t);
507
                                        return 1;
508
                                } else { // se odbije
509
                                        transmit = in.GetN() + sign_n*2*cosTi*n;
510
                                        out->Set(intersect, transmit);
511
                                        pol_t = -in.GetP() + sign_n*2*cosTi*n;
512
                                        out->SetPolarization(pol_t);
513
                                        return 2;
514
                                }                              
515
                        }
516
                        break;
517
                // ----------------------------------------------------------------------------
518
                // --------------- reflection at "reflection" probability ---------------------
519
                // ----------------------------------------------------------------------------
520
                case SURF_REFLE:
521
                        p_ref = rand.Uniform(0.0, 1.0);
522
                        if(p_ref < reflection) { // se odbije
523
                                cosTi = in.GetN() * n;
524
                                transmit = in.GetN() - 2*cosTi*n;
525
                                out->Set(intersect, transmit);
526
                                return 2; //sdhfvjhsdbfjhsdbcvjhsb
527
                        } else { // se ne odbije
528
                                transmit = in.GetN();
529
                                out->Set(intersect, transmit);         
530
                                return -2;
531
                        }                                              
532
                        break;
533
 
534
                // total reflection from n1 to n2 with R probbability
535
                case SURF_IMPER:
536
                        p_ref = rand.Uniform(0.0, 1.0);        
537
                        if(p_ref < reflection) { // se odbije
538
                                cosTi = in.GetN() * n;                 
539
                                if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj
540
                                        transmit = in.GetN() - 2*cosTi*n;
541
                                        out->Set(intersect, transmit);
542
                                } else { // ni tot. odboja
543
                                        transmit = in.GetN();
544
                                        out->Set(intersect, transmit);
545
                                        return -2;
546
                                }
547
                        } else { // se ne odbije
548
                                transmit = in.GetN();
549
                                out->Set(intersect, transmit);                 
550
                                return -2;
551
                        }                                              
552
                        break;
553
 
554
                default:
555
                        *out = in;
556
                        break;
557
        }
558
 
559
        return 1;
560
}
561
//=================================================================================
562
 
563
 
564
//=================================================================================
565
Guide::Guide(TVector3 center0, DetectorParameters &parameters)
566
{
567
double t;
568
 
569
        TDatime now; rand.SetSeed(now.Get());
570
 
571
        center = center0;
572
        double b = parameters.getB();
573
        double a = parameters.getA();
574
        _d = parameters.getD();
575
        n1 = parameters.getN1();
576
        n2 = parameters.getN2();
577
        // if PlateOn, then n0 = n3 (optical grease), else = n1 (air)
578
        double n0 = (parameters.getPlateOn() ? parameters.getN3(): n1);
579
        n3 = parameters.getN3();
580
        _r = c_reflectivity;
581
        int fresnel = parameters.getFresnel();
582
 
583
        t = b/2.0;
584
        vodnik_edge[0].SetXYZ(0.0, t,-t); vodnik_edge[1].SetXYZ(0.0, t, t);
585
        vodnik_edge[2].SetXYZ(0.0,-t, t); vodnik_edge[3].SetXYZ(0.0,-t,-t);
586
        t = a/2.0;
587
        vodnik_edge[4].SetXYZ(_d, t,-t); vodnik_edge[5].SetXYZ(_d, t, t);
588
        vodnik_edge[6].SetXYZ(_d,-t, t); vodnik_edge[7].SetXYZ(_d,-t,-t);
589
 
590
        for(int i = 0; i<8; i++) vodnik_edge[i] += center;
591
 
592
        s_side[0] = new CSurface(SURF_REFRA, vodnik_edge, n0, n2, _r);
593
        s_side[0]->FlipN();
594
 
595
        s_side[1] = new CSurface(SURF_REFRA, vodnik_edge[3], vodnik_edge[2], vodnik_edge[6], vodnik_edge[7], n2, n1, _r);
596
        s_side[2] = new CSurface(SURF_REFRA, vodnik_edge[2], vodnik_edge[1], vodnik_edge[5], vodnik_edge[6], n2, n1, _r);
597
        s_side[3] = new CSurface(SURF_REFRA, vodnik_edge[1], vodnik_edge[0], vodnik_edge[4], vodnik_edge[5], n2, n1, _r);
598
        s_side[4] = new CSurface(SURF_REFRA, vodnik_edge[0], vodnik_edge[3], vodnik_edge[7], vodnik_edge[4], n2, n1, _r);
599
 
600
        s_side[5] = new CSurface(SURF_REFRA, &vodnik_edge[4], n2, n3, _r);
601
        s_side[5]->FlipN();
602
 
603
        if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1);
604
 
605
        hfate = (TH1F*)gROOT->FindObject("hfate"); if(hfate) delete hfate;
606
        hfate = new TH1F("hfate", "Ray fate", 8, -3.5, 4.5);
607
        (hfate->GetXaxis())->SetBinLabel(1, "Back Ref");
608
        (hfate->GetXaxis())->SetBinLabel(2, "Ref Loss");
609
        (hfate->GetXaxis())->SetBinLabel(3, "No Tot Refl");
610
        (hfate->GetXaxis())->SetBinLabel(4, "LG Miss");
611
        (hfate->GetXaxis())->SetBinLabel(5, "Exit");
612
        (hfate->GetXaxis())->SetBinLabel(6, "Enter");
613
        (hfate->GetXaxis())->SetBinLabel(7, "Rays");
614
        (hfate->GetXaxis())->SetBinLabel(8, "Absorb");
615
 
616
        hnodb_all = (TH1F*)gROOT->FindObject("hnodb_all"); if(hnodb_all) delete hnodb_all;
617
        hnodb_all = new TH1F("hnodb_all", "N reflected", VODNIK_MAX_ODB, -0.5, VODNIK_MAX_ODB-0.5);
618
 
619
        hnodb_exit = (TH1F*)gROOT->FindObject("hnodb_exit"); if(hnodb_exit) delete hnodb_exit;
620
        hnodb_exit = new TH1F("hnodb_exit", "N reflected and exit", VODNIK_MAX_ODB, -0.5, VODNIK_MAX_ODB-0.5);
621
 
622
        int nBins = nch + 1;
623
        hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin;
624
        hin = new TH2F("hin", "Guide entrance window", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0);
625
 
626
        hout = (TH2F*)gROOT->FindObject("hout"); if(hout) delete hout;
627
        hout = new TH2F("hout", "Guide exit window", nBins, -a/2.0, +a/2.0, nBins, -a/2.0, +a/2.0);
628
 
629
        absorption = 0;
630
        A = 0;
631
}
632
//-----------------------------------------------------------------------------
633
// Sledi zarku skozi vodnik. Vrne:                                             
634
//  0, ce zgresi vstopno ploskev                                               
635
//  1, ce zadane izstopno ploskev                                              
636
// -1, ce se v vodniku ne odbije totalno 
637
//  2, enter the light guide, bin 2 of hfate = refraction                                     
638
// -2, ce se ne odbije zaradi koncnega R stranic                               
639
// -3, ce se odbije nazaj in gre nazaj ven skozi sprednjo ploskev              
640
// +4, ce se absorbira v materialu                                             
641
Fate Guide::PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
642
{
643
  if (dbg) printf("--- GUIDE::PropagateRay ---\n");
644
  CRay ray0;
645
  CRay ray1;
646
  TVector3 vec0, vec1;
647
  //int fate = 11;
648
  int inters_i = 0;
649
 
650
        ray0 = in;
651
        int n_odb = 0;
652
        int last_hit = 0;
653
        int propagation = 0;
654
        if( !(s_side[0]->PropagateRay(ray0, &ray1, &vec1)) ) {
655
                // ce -NI- presecisca z vstopno
656
                if (dbg) printf("  GUIDE: missed the light guide\n");
657
                fate = missed;
658
                hfate->Fill(0);
659
        } else {
660
          if (dbg) printf("  GUIDE: ray entered\n");
661
                points[0] = ray1.GetR();
662
                hfate->Fill(2); // enter
663
                hin->Fill(vec1.y(), vec1.z());
664
                if (dbg) printf("  GUIDE: n_odb = %d\n", n_odb);
665
                while (n_odb++ < VODNIK_MAX_ODB) {
666
                  if (dbg) printf("  GUIDE: Boundary test: %d\n",n_odb);
667
                        ray0 = ray1;
668
                        vec0 = vec1;
669
                        propagation = 11;
670
                        for(inters_i=0; inters_i<6; inters_i++) {
671
                          if (dbg) printf("  GUIDE: Test intersection with surface %d \n", inters_i);
672
                                if( inters_i != last_hit) {
673
                                  int testBoundary = s_side[inters_i]->TestIntersection(&vec1, ray1);
674
                                        if( testBoundary ) {
675
                                          if (dbg) printf("  GUIDE: ray intersects with LG surface %d\n",inters_i);
676
                                                break;
677
                                  }
678
                                }  
679
                        }                      
680
                        points[n_odb] = vec1;
681
                        if(inters_i == 0) {
682
                          fate = backreflected;
683
                          hfate->Fill(backreflected);
684
                          break;
685
                        } // backreflection
686
 
687
                        propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
688
                        if (dbg) printf("  GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
689
                        if(inters_i == 5) { // successfull exit
690
                        // check on which side the vector is?
691
                          TVector3 ray = ray1.GetN();
692
                          TVector3 exitNormal = s_side[5]->GetN();
693
                          //printf("theta(ray) = %lf, theta(normal5) = %lf ", ray.Theta()*DEGREE, exitNormal.Theta()*DEGREE);
694
                          //printf("phi(ray) = %lf, phi(normal5) = %lf\n", ray.Phi()*DEGREE, exitNormal.Phi()*DEGREE);
695
                          if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal);
696
                          if (ray*exitNormal > 0) {
697
                            if (dbg) printf("  GUIDE: ray is backreflected from exit window.\n");
698
                            fate = backreflected;
699
                            n_odb++;
700
                            points[n_odb] = vec1;
701
                            ray0 = ray1;
702
                            break;
703
                          }
704
                                fate =  hit;
705
                                hout->Fill(vec1.y(), vec1.z());
706
                                hnodb_exit->Fill(n_odb-1);
707
                                n_odb++;
708
                                points[n_odb] = vec1;
709
                                ray0 = ray1;
710
                                break;
711
                        }
712
                        if(propagation == REFRACTION) {
713
                          fate = refracted;
714
                          n_odb++;
715
                          points[n_odb] = vec1;
716
                          ray0 = ray1;
717
                          break;
718
                        } // no total reflection when should be
719
                        if(propagation == -2) {
720
                          fate = noreflection;
721
                          break;
722
                        } //refraction due to finite reflectivity
723
 
724
                        last_hit = inters_i;
725
                }      
726
        }
727
 
728
        //--- material absorption ---
729
        if(absorption) {
730
                double travel = 0.0;
731
                printf("n_odb = %d\n", n_odb); //dbg   
732
                for(int point = 0; point < n_odb-1; point++) {
733
                        travel += (points[point] - points[point+1]).Mag();
734
                        printf("travel = %lf\n", travel); //dbg   
735
                }
736
                double T_abs = TMath::Exp(-travel/A);
737
                printf("T_abs = %lf\n", T_abs); //dbg   
738
                double p_abs = rand.Uniform(0.0, 1.0); 
739
                printf("p_abs = %lf\n", p_abs); //dbg   
740
 
741
                if(p_abs > T_abs) fate = absorbed; // absorption
742
        }
743
        //--- material absorption ---
744
 
745
        hfate->Fill(fate);
746
        hfate->Fill(3);
747
        hnodb_all->Fill(n_odb-2);
748
        *n_points = n_odb+1;
749
        *out = ray0;
750
        return fate;
751
}
752
//-----------------------------------------------------------------------------
753
void Guide::GetVFate(int *out)
754
{
755
        for(int i=0;i<7;i++) out[i] = (int)hfate->GetBinContent(i+1);
756
}
757
//-----------------------------------------------------------------------------
758
void Guide::Draw(int color, int width)
759
{
760
        for(int i = 0; i<6; i++) s_side[i]->Draw(color, width);
761
}
762
//-----------------------------------------------------------------------------
763
void Guide::DrawSkel(int color, int width)
764
{
765
        TPolyLine3D *line3d = new TPolyLine3D(2);
766
        line3d->SetLineWidth(width); line3d->SetLineColor(color);
767
 
768
        for(int i=0; i<4; i++) {
769
                line3d->SetPoint(0, vodnik_edge[i+0].x(), vodnik_edge[i+0].y(), vodnik_edge[i+0].z());
770
                line3d->SetPoint(1, vodnik_edge[i+4].x(), vodnik_edge[i+4].y(), vodnik_edge[i+4].z());
771
                line3d->DrawClone();
772
        }
773
}
774
//=================================================================================
775
 
776
//=================================================================================
777
int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
778
{
779
        double num, den; //stevec, imenovalec
780
        double t;
781
        TVector3 tmp;
782
 
783
        if(dbg) printf("---> CPlaneR::TestIntersection <---\n");
784
        if(dbg) {printf("c = "); printv(center); printf(" | n = "); printv(n); printf("\n");}
785
 
786
        double D = - n*center; 
787
        num = n*ray.GetR() + D;
788
        den = n*ray.GetN();
789
 
790
        if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den);
791
 
792
        if(TMath::Abs(den) < MARGIN) {
793
                if(TMath::Abs(num) < MARGIN)
794
                        return 0;
795
                else
796
                        return 0;
797
        }
798
 
799
        t = num / den;
800
 
801
        if(dbg) printf("t = %.4lf | ", t);
802
 
803
        tmp = ray.GetR();
804
        tmp -= t*ray.GetN();
805
        *vec = tmp;
806
 
807
        if(dbg) {printv(tmp); printf(" | Rv = %.4lf <> R = %.4lf\n", ((tmp - center).Mag()), _r);}
808
 
809
 
810
        if( ((tmp - center).Mag()) < _r )
811
                return 1;
812
        else
813
                return 0;
814
}
815
//-----------------------------------------------------------------------------
816
void CPlaneR::Draw(int color, int width)
817
{
818
        const int NN = 32;     
819
        double phi, x, y;
820
 
821
        TPolyLine3D *arc;
822
        arc = new TPolyLine3D(NN+1);
823
        arc->SetLineWidth(width);
824
        arc->SetLineColor(color);
825
 
826
        for(int i=0; i<=NN; i++) {
827
                phi = i*2.0*TMath::Pi()/NN;
828
                x = _r*TMath::Cos(phi);
829
                y = _r*TMath::Sin(phi);
830
                arc->SetPoint(i, center.x(),  x,  y);
831
        }
832
        arc->Draw();
833
}
834
//=================================================================================
835
 
836
 
837
//=================================================================================
838
CDetector::CDetector(TVector3 center0, DetectorParameters& parameters) :
839
  center(center0),
840
  glass_on(parameters.getGlassOn()),
841
  glass_d(parameters.getGlassD()),
842
  //x_gap(parameters.getGap().X()),
843
  //y_gap(parameters.getGap().Y()),
844
  //z_gap(parameters.getGap().Z()),
845
        //glass(new CSurface),
846
        //glass_circle(new CPlaneR),
847
        //active(new CPlane4),
848
        col_in(2),
849
        col_lg(8),
850
        col_out(4),
851
        col_rgla(6),
852
        col_LG(1),
853
        col_glass(4),
854
        col_active(7),
855
        guide_on(parameters.getGuideOn()),
856
        //window_R( parameters.getB() ),
857
        //window_d(0),
858
  guide(new Guide(center0, parameters)),
859
  plate(new Plate(parameters)),
860
  _plateWidth(parameters.getPlateWidth()),
861
  _plateOn(parameters.getPlateOn())
862
  {
863
//  };
864
 
865
//-----------------------------------------------------------------------------
866
//void CDetector::Init()
867
//{
868
  double d = parameters.getD();
869
        double x_offset;
870
        if(guide_on) x_offset = center.x();
871
        else x_offset = center.x() - d;
872
 
873
        //guide = new CVodnik(center, SiPM, M, d, type_in, type_side, type_out, n1, n2, n3, reflectivity, fresnel, absorption, A);
874
 
875
        double b = parameters.getB();
876
        double n1 = parameters.getN1();
877
        double n2 = parameters.getN2();
878
        //double n3 = parameters.getN3();
879
        double reflectivity = c_reflectivity; // for faster simulation, not using Fresnel eqs.
880
        double x_gap = parameters.getGap().X();
881
        double y_gap = parameters.getGap().Y();
882
        double z_gap = parameters.getGap().Z();
883
 
884
        TVector3 plane_v[4];
885
        int nBins = nch + 1;
886
        double p_size = b/2.0;
887
        plane_v[0].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap - p_size);
888
        plane_v[1].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap + p_size);
889
        plane_v[2].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap + p_size);
890
        plane_v[3].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap - p_size);
891
        glass = new CSurface(SURF_REFRA, plane_v, n2, n1, reflectivity); glass->FlipN();
892
 
893
        glass_circle = new CPlaneR(TVector3(x_offset+d+glass_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), b);
894
 
895
        hglass = (TH2F*)gROOT->FindObject("hglass"); if(hglass) delete hglass;
896
        hglass = new TH2F("hglass", "Hits glass",
897
                                          nBins, y_gap - p_size, y_gap + p_size,
898
                  nBins, z_gap - p_size, z_gap + p_size);
899
 
900
 
901
        p_size = parameters.getActive()/2.0;
902
        //cout<<"SiPM active length "<<detectorActive<<endl;
903
        //p_size = 1.0/2.0;
904
        plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size);
905
        plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size);
906
        plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size);
907
        plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size);
908
        active = new CPlane4(plane_v);
909
 
910
        hactive = (TH2F*)gROOT->FindObject("hactive"); if(hactive) delete hactive;
911
        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);
912
 
913
 
914
        p_size = b/2.0;
915
        //p_size = 2.5;
916
        //p_size = M*0.6;
917
        hlaser = (TH2F*)gROOT->FindObject("hlaser"); if(hlaser) delete hlaser;
918
        hlaser = new TH2F("hlaser", "Hits at LG entrance", nBins, -p_size, p_size, nBins, -p_size, p_size);
919
 
920
 
921
        p_size = 2*b/2.0;
922
        plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size);
923
        plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size);
924
        plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size);
925
        plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size);
926
        detector = new CPlane4(plane_v);
927
 
928
        hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector;
929
        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);
930
 
931
 
932
        /*
933
        window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R);  
934
 
935
        p_size = M*a;
936
        plane_v[0].SetXYZ(x_offset+d+window_d, y_gap + p_size, z_gap - p_size);
937
        plane_v[1].SetXYZ(x_offset+d+window_d, y_gap + p_size, z_gap + p_size);
938
        plane_v[2].SetXYZ(x_offset+d+window_d, y_gap - p_size, z_gap + p_size);
939
        plane_v[3].SetXYZ(x_offset+d+window_d, y_gap - p_size, z_gap - p_size);
940
        window = new CSurface(SURF_REFRA, plane_v, n1, n2, reflectivity); window->FlipN();
941
 
942
        hwindow = (TH2F*)gROOT->FindObject("hwindow"); if(hwindow) delete hwindow;
943
        hwindow = new TH2F("hwindow", "Hits Window", nch, y_gap - window_R, y_gap + window_R, nch, z_gap - window_R, z_gap + window_R);
944
        */
945
        p_size = b/2.0;
946
        histoPlate = (TH2F*)gROOT->FindObject("histoPlate"); if(histoPlate) delete histoPlate;
947
        histoPlate = new TH2F("histoPlate", "Hits on glass plate", nBins, -p_size, +p_size, nBins, -p_size, +p_size);
948
}
949
//-----------------------------------------------------------------------------
950
// vrne 1 ce je zadel aktvino povrsino
951
// vrne <1 ce jo zgresi
952
int CDetector::Propagate(CRay in, CRay *out, int draw)
953
// Sledi zarku skozi vodnik. Vrne:                                             
954
//  0, ce zgresi vstopno ploskev MISSED                                              
955
//  1, ce zadane izstopno ploskev HIT                                             
956
// -1, ce se v vodniku ne odbije totalno REFRACTED
957
//  2, enter the light guide, bin 2 of hfate EXIT                                     
958
// -2, ce se ne odbije zaradi koncnega R stranic - no total reflection REFRACTED                             
959
// -3, ce se odbije nazaj in gre nazaj ven skozi sprednjo ploskev BACK_REFLECTED             
960
// +4, ce se absorbira v materialu ABSORBED
961
{
962
  if (dbg) printf("--- Detector::Propagate ---\n");
963
        //CRay *ray0 = new CRay; ray0->Set(in.GetR(), in.GetN()); ray0->SetColor(col_in);
964
        CRay *rayin = new CRay(in);
965
        rayin->SetColor(col_in);
966
        CRay *rayout = new CRay;
967
 
968
        const int max_n_points = guide->GetMAXODB() + 2;
969
        TVector3 pointsPlate[max_n_points];
970
        //TVector3 intersection;
971
        int fatePlate;
972
        int nPointsPlate;
973
        TPolyLine3D *line3d = new TPolyLine3D(2);
974
        line3d->SetLineWidth(1);
975
        line3d->SetLineColor(4);
976
 
977
        // Draw the plate and propagate the ray through
978
        // check if the ray should be reflected??
979
 
980
        if(_plateOn) {
981
 
982
          fatePlate = plate->propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
983
          if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
984
 
985
          if(draw) {
986
            if(fatePlate == 0) {
987
              rayout->SetColor(col_in);
988
              rayout->DrawS(center.x() - _plateWidth, -10.0);
989
              }
990
              else {
991
                int p_i;
992
                for(p_i = 0; p_i < nPointsPlate-1; p_i++) {
993
                                          line3d->SetPoint(0, pointsPlate[p_i].x(), pointsPlate[p_i].y(), pointsPlate[p_i].z());
994
                                          line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z());
995
                                          line3d->DrawClone();
996
                                  }
997
                                  if(fatePlate != -2) {
998
                                          //rayout->SetColor(7);
999
                                          rayout->DrawS(pointsPlate[p_i].x(), -0.1);
1000
                                  }
1001
                          }
1002
            }
1003
 
1004
          if(! (fatePlate == 1 or fatePlate == -1) ) {
1005
 
1006
                        cout<<"CDetector::propagate Simulated ray missed the entry surface!"<<endl;
1007
                        return fatePlate;
1008
                }
1009
                histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
1010
         }
1011
         else {
1012
          rayout = rayin;
1013
          if(draw) rayout->DrawS(center.x(), -10.0);
1014
          }
1015
 
1016
        // If the ray is not reflected in the plate
1017
        // Draw the light guide and propagate the ray through
1018
 
1019
        //const int max_n_points = guide->GetMAXODB() + 2;
1020
        TVector3 points[max_n_points];
1021
        TVector3 presecisce;
1022
        //int fate;
1023
        int n_points;
1024
        int fate_glass;
1025
        CRay *ray0 = new CRay(*rayout);
1026
        // delete rayout; -> creates dangling reference when tries to delete ray0!
1027
        //delete rayin;
1028
        CRay *ray1 = new CRay;
1029
 
1030
        fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
1031
 
1032
        line3d->SetLineColor(col_lg);  
1033
        int p_i;
1034
        if(guide_on) {
1035
                if(draw) {
1036
                        if(fate == missed) {
1037
                          if (dbg) printf("Detector: fate=missed\n");
1038
                          TVector3 r = ray1->GetR();
1039
                          TVector3 n = ray1->GetN();
1040
                          ray1->Set(r,n);
1041
                                ray1->DrawS(center.x(), 10.0);
1042
                        } else {
1043
                                for(p_i = 0; p_i < n_points-1; p_i++) {
1044
                                        line3d->SetPoint(0, points[p_i].x(), points[p_i].y(), points[p_i].z());
1045
                                        line3d->SetPoint(1, points[p_i+1].x(), points[p_i+1].y(), points[p_i+1].z());
1046
                                        line3d->DrawClone();
1047
                                }
1048
                                if(fate != noreflection) {
1049
                                  if (dbg) printf("Detector: fate != noreflection, fate = %d\n", (int)fate);
1050
                                        if(glass_on) {/*if(fate == 1)*/ ray1->Draw(points[p_i].x(), center.x() + guide->getD() + glass_d);}
1051
                                        else {
1052
                                        ray1->SetColor(col_out);
1053
                                        ray1->DrawS(points[p_i].x(), 10.0);
1054
                                        }
1055
                                }
1056
                        }
1057
                }
1058
 
1059
 
1060
                if(! (fate == hit or fate == refracted) ) {
1061
                if (dbg) printf("Detector: fate != hit, refracted\n");
1062
                        *out = *ray1;
1063
                        return fate;
1064
                }
1065
        } else {
1066
          if (dbg) printf("Detector: fate = hit or refracted");
1067
                ray1 = ray0;
1068
                if(draw) {
1069
                        if(glass_on) ray1->Draw(center.x(), center.x() /*+ window_d*/);
1070
                        else ray1->DrawS(center.x(), 10.0);
1071
                }
1072
        }
1073
 
1074
/*     
1075
        TVector3 pres_wind;
1076
        fate = window_circle->TestIntersection(&pres_wind, *ray1);
1077
        if(fate == 1) {
1078
                hwindow->Fill(pres_wind.y(), pres_wind.z());
1079
 
1080
                if(!guide_on) {
1081
                        window->PropagateRay(*ray0, ray1, &presecisce);
1082
                        if(draw) ray1->Draw(center.x() + window_d, center.x() + glass_d);
1083
                        *ray0 = *ray1;
1084
                }
1085
        */     
1086
                fate = missed; // zgresil aktivno povrsino
1087
                if(glass_on) {                 
1088
                        *ray0 = *ray1; ray1->SetColor(col_rgla);
1089
                        fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
1090
                        if(fate_glass == 1) {
1091
                                hglass->Fill(presecisce.y(), presecisce.z());
1092
                                if(draw) ray1->DrawS(presecisce.x(), 10.0);
1093
                                if(active->TestIntersection(&presecisce, *ray1)) {
1094
                                        fate = hit;
1095
                                        hactive->Fill(presecisce.y(), presecisce.z());
1096
                                        hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
1097
                                }
1098
                                if(detector->TestIntersection(&presecisce, *ray1))
1099
                                        hdetector->Fill(presecisce.y(), presecisce.z());
1100
                        } else if(fate_glass == 2) {
1101
                                if(draw) ray1->DrawS(presecisce.x(), 10.0);
1102
                        }
1103
                } else {
1104
                        if(active->TestIntersection(&presecisce, *ray1)) {
1105
                                fate = hit;
1106
                                hactive->Fill(presecisce.y(), presecisce.z());
1107
                                hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
1108
                        }
1109
                        if(detector->TestIntersection(&presecisce, *ray1))
1110
                                hdetector->Fill(presecisce.y(), presecisce.z());
1111
                }
1112
        //} else {
1113
                //if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
1114
        //}
1115
 
1116
        *out = *ray1;
1117
        //delete ray0;
1118
        //delete ray1;
1119
        return fate;
1120
}
1121
//-----------------------------------------------------------------------------
1122
void CDetector::Draw(int width)
1123
{
1124
        if(guide_on) {
1125
                if( TMath::Abs(guide->getN1()-guide->getN2()) < MARGIN ) {
1126
                  if(_plateOn) plate->drawSkel(col_LG, width);
1127
                        guide->DrawSkel(col_LG, width);
1128
                        }
1129
                else {
1130
                  if(_plateOn) plate->draw(4, width);
1131
                        guide->Draw(col_LG, width);
1132
                        }
1133
        }
1134
 
1135
        if(glass_on) glass_circle->Draw(col_glass, width);
1136
        //window_circle->Draw(col_glass, width);
1137
        active->Draw(col_active, width);
1138
}
1139
//=================================================================================
1140
 
1141
Plate::Plate(DetectorParameters& parameters)
1142
{
1143
  TVector3 center = CENTER;
1144
  const double b = parameters.getB();
1145
  const double n1 = parameters.getN1();
1146
        const double n2 = parameters.getN2();
1147
        const double n3 = parameters.getN3();
1148
        const double reflectivity = c_reflectivity;
1149
        const double t = b/2.;
1150
        const double plateWidth = parameters.getPlateWidth();
1151
        center.SetX( CENTER.X() - plateWidth );
1152
        plate_edge[0].SetXYZ(0.0, t,-t); plate_edge[1].SetXYZ(0.0, t, t);
1153
        plate_edge[2].SetXYZ(0.0,-t, t); plate_edge[3].SetXYZ(0.0,-t,-t);
1154
        plate_edge[4].SetXYZ(plateWidth, t,-t); plate_edge[5].SetXYZ(plateWidth, t, t);
1155
        plate_edge[6].SetXYZ(plateWidth,-t, t); plate_edge[7].SetXYZ(plateWidth,-t,-t);
1156
 
1157
        for(int i = 0; i<8; i++) plate_edge[i] += center;
1158
 
1159
        sides[0] = new CSurface(SURF_REFRA, plate_edge, n1, n2, reflectivity);
1160
        sides[0]->FlipN();
1161
 
1162
        sides[1] = new CSurface(SURF_REFRA, plate_edge[3], plate_edge[2], plate_edge[6], plate_edge[7], n2, n2, reflectivity);
1163
        sides[2] = new CSurface(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, reflectivity);
1164
        sides[3] = new CSurface(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, reflectivity);
1165
        sides[4] = new CSurface(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, reflectivity);
1166
 
1167
        sides[5] = new CSurface(SURF_REFRA, &plate_edge[4], n2, n3, reflectivity);
1168
        sides[5]->FlipN();
1169
 
1170
        for(int i=0; i<6; i++) sides[i]->SetFresnel(1);
1171
}
1172
 
1173
void Plate::draw(int color, int width)
1174
{
1175
        for(int i = 0; i<6; i++) sides[i]->Draw(color, width);
1176
}
1177
 
1178
void Plate::drawSkel(int color, int width)
1179
{
1180
        TPolyLine3D line3d(2);
1181
        line3d.SetLineWidth(width);
1182
        line3d.SetLineColor(color);
1183
 
1184
        for(int i=0; i<4; i++) {
1185
                line3d.SetPoint(0, plate_edge[i+0].x(), plate_edge[i+0].y(), plate_edge[i+0].z());
1186
                line3d.SetPoint(1, plate_edge[i+4].x(), plate_edge[i+4].y(), plate_edge[i+4].z());
1187
                line3d.DrawClone();
1188
        }
1189
}
1190
 
1191
int Plate::propagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
1192
{
1193
  CRay ray0;
1194
  CRay ray1;
1195
  TVector3 vec0, vec1;
1196
  int fate = 11;
1197
  int inters_i = 0;
1198
 
1199
        ray0 = in;
1200
        int n_odb = 0;
1201
        int last_hit = 0;
1202
        int propagation = 0;
1203
 
1204
        if( !(sides[0]->PropagateRay(ray0, &ray1, &vec1)) ) {
1205
                // ce -NI- presecisca z vstopno
1206
                fate = 0;
1207
        } else {
1208
                points[0] = ray1.GetR();
1209
                //hfate->Fill(2);
1210
                //hin->Fill(vec1.y(), vec1.z());
1211
                while (n_odb++ < VODNIK_MAX_ODB) {
1212
                        ray0 = ray1;
1213
                        vec0 = vec1;
1214
                        propagation = 11;
1215
                        for(inters_i=0; inters_i<6; inters_i++) {
1216
                                if( inters_i != last_hit) {
1217
                                        if( sides[inters_i]->TestIntersection(&vec1, ray1) ) break;
1218
                                        }
1219
                                }                                                      
1220
                        points[n_odb] = vec1;
1221
                        if(inters_i == 0) {
1222
                          fate = -3;
1223
                          break;} // backreflection
1224
 
1225
                        propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
1226
                        if(inters_i == 5) { // successfull exit
1227
                                fate =  1;
1228
                                //hout->Fill(vec1.y(), vec1.z()); 
1229
                                //hnodb_exit->Fill(n_odb-1); 
1230
                                n_odb++;
1231
                                points[n_odb] = vec1;
1232
                                ray0 = ray1;
1233
                                break;
1234
                        }
1235
                        if(propagation == 1) {
1236
                          fate = -1;
1237
                          n_odb++;
1238
                          points[n_odb] = vec1;
1239
                          ray0 = ray1;
1240
                          break;} // no total reflection when should be
1241
                        if(propagation == -2) {fate = -2; break;} // absorption due to finite reflectivity
1242
 
1243
                        last_hit = inters_i;
1244
                }      
1245
        }
1246
 
1247
        *n_points = n_odb+1;
1248
        *out = ray0;
1249
        return fate;
1250
};
1251
//=============================================================================================================================== <<<<<<<<
1252
 
1253