Subversion Repositories f9daq

Rev

Rev 25 | Rev 70 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 25 Rev 54
Line 2... Line 2...
2 
 
 2 
 
3 
#include <iostream>
 3 
#include <iostream>
4 
 
 - 
 
5 
 
 4 
 
6 
// vector output shortcut
 5 
// vector output shortcut
7 
void printv(TVector3 v)
 6 
void printv(TVector3 v)
8 
{
 7 
{
9 
        printf("(x,y,z) = (%.4lf, %.4lf, %.4lf)\n", v.x(), v.y(), v.z());
 8 
        printf("(x,y,z) = (%.4lf, %.4lf, %.4lf)\n", v.x(), v.y(), v.z());
Line 31... Line 30...
31 
//-----------------------------------------------------------------------------
 30 
//-----------------------------------------------------------------------------
32 
//void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0)
 31 
//void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0)
33 
//{
 32 
//{
34 
        //r.SetXYZ(x0, y0, z0);
 33 
        //r.SetXYZ(x0, y0, z0);
35 
        //n.SetXYZ(l0, m0, n0); n = n.Unit();
 34 
        //n.SetXYZ(l0, m0, n0); n = n.Unit();
36 
//}
 35 
//}
37 
//-----------------------------------------------------------------------------
 36 
//-----------------------------------------------------------------------------
38 
/*
 37 
/*
39 
CRay& CRay::operator = (const CRay& p)
 38 
CRay& CRay::operator = (const CRay& p)
40 
{
 39 
{
41 
        r.SetXYZ(p.GetR().x(), p.GetR().y(), p.GetR().z());
 40 
        r.SetXYZ(p.GetR().x(), p.GetR().y(), p.GetR().z());
42 
        //this->r.SetXYZ(p.x(), p.y(), p.z());
 41 
        //this->r.SetXYZ(p.x(), p.y(), p.z());
Line 47... Line 46...
47 
void CRay::Print()
 46 
void CRay::Print()
48 
{
 47 
{
49 
        printf("---> CRay::Print() <---\n");
 48 
        printf("---> CRay::Print() <---\n");
50 
        printf("(x,y,z)=(%.2lf, %.2lf, %.2lf); (l,m,n)=(%.2lf, %.2lf, %.2lf)\n",
 49 
        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());
 50 
                r.x(), r.y(), r.z(), n.x(), n.y(), n.z());
52 
}
 51 
}
53 
//-----------------------------------------------------------------------------
 52 
//-----------------------------------------------------------------------------
54 
void CRay::Draw()
 53 
void CRay::Draw()
55 
{
 54 
{
56 
double t = 50.0;
 55 
double t = 50.0;
57 
TPolyLine3D *line3d = new TPolyLine3D(2);
 56 
TPolyLine3D *line3d = new TPolyLine3D(2);
58 
        //line3d->SetPoint(0, r.x() - t*n.x(), r.y() - t*n.y(), r.z() - t*n.z());
 57 
        //line3d->SetPoint(0, r.x() - t*n.x(), r.y() - t*n.y(), r.z() - t*n.z());
Line 63... Line 62...
63 
 
 62 
 
64 
        line3d->Draw();
 63 
        line3d->Draw();
65 
}
 64 
}
66 
//-----------------------------------------------------------------------------
 65 
//-----------------------------------------------------------------------------
67 
void CRay::Draw(double x_from, double x_to)
 66 
void CRay::Draw(double x_from, double x_to)
68 
{
 67 
{
69 
double A1, A2;
 68 
double A1, A2;
70 
  TPolyLine3D *line3d = new TPolyLine3D(2);
 69 
  TPolyLine3D *line3d = new TPolyLine3D(2);
71 
 
 70 
 
72 
        if(n.x() < MARGIN) {
 71 
        if(n.x() < MARGIN) {
73 
                A1 = A2 = 0.0;
 72 
                A1 = A2 = 0.0;
74 
        } else {
 73 
        } else {
75 
                A1 = (x_from - r.x())/n.x();   
 74 
                A1 = (x_from - r.x())/n.x();   
76 
                A2 = (x_to - r.x())/n.x();
 75 
                A2 = (x_to - r.x())/n.x();
77 
        }
 76 
        }
78 
       
 77 
       
79 
        line3d->SetPoint(0, x_from, A1*n.y()+r.y(), A1*n.z()+r.z());
 78 
        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());
 79 
        line3d->SetPoint(1, x_to, A2*n.y()+r.y(), A2*n.z()+r.z());
81 
        line3d->SetLineWidth(1);
 80 
        line3d->SetLineWidth(1);
82 
        line3d->SetLineColor(color);
 81 
        line3d->SetLineColor(color);
83 
 
 82 
 
84 
        line3d->Draw();
 83 
        line3d->Draw();
85 
}
 84 
}
86 
//-----------------------------------------------------------------------------
 85 
//-----------------------------------------------------------------------------
87 
void CRay::DrawS(double x_from, double t)
 86 
void CRay::DrawS(double x_from, double t)
88 
{
 87 
{
89 
        double A1;
 88 
        double A1;
Line 127... Line 126...
127 
// za izracun parametrov ravnine je en vektor prevec, vendar tega
126 
// za izracun parametrov ravnine je en vektor prevec, vendar tega
128 
// rabim kot zadnji vogal poligona
 127 
// rabim kot zadnji vogal poligona
129 
//void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
 128 
//void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
130 
//{
 129 
//{
131 
double x1,y1,z1, x2,y2,z2, x3,y3,z3;
 130 
double x1,y1,z1, x2,y2,z2, x3,y3,z3;
132 
       
 131 
       
133 
        x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
 132 
        x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
134 
        x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
 133 
        x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
135 
        x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
 134 
        x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
136 
 
 135 
 
137 
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
 136 
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
138 
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
 137 
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
139 
        C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
 138 
        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);
 139 
        D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
141 
 
 140 
 
142 
        r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
 141 
        r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
143 
        n.SetXYZ(A, B, C);
 142 
        n.SetXYZ(A, B, C);
144 
        n = n.Unit();
 143 
        n = n.Unit();
145 
 
 144 
 
146 
        for(int i=0;i<4;i++)
 145 
        for(int i=0;i<4;i++)
147 
                edge[i] = r[i-3 ? i+1 : 0] - r[i];
 146 
                edge[i] = r[i-3 ? i+1 : 0] - r[i];
148 
 
 147 
 
149 
        for(int i=0;i<4;i++)
 148 
        for(int i=0;i<4;i++)
150 
                angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
 149 
                angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
151 
};
 150 
};
152 
 
 151 
 
153 
void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
 152 
void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
154 
{
 153 
{
155 
  double x1,y1,z1, x2,y2,z2, x3,y3,z3;
 154 
  double x1,y1,z1, x2,y2,z2, x3,y3,z3;
156 
       
 155 
       
157 
        x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
 156 
        x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
158 
        x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
 157 
        x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
159 
        x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
 158 
        x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
160 
 
 159 
 
161 
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
 160 
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
162 
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
 161 
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
163 
        C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
 162 
        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);
 163 
        D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
165 
 
 164 
 
166 
        r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
 165 
        r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
167 
        n.SetXYZ(A, B, C);
 166 
        n.SetXYZ(A, B, C);
168 
        n = n.Unit();
 167 
        n = n.Unit();
169 
 
 168 
 
170 
        for(int i=0;i<4;i++)
 169 
        for(int i=0;i<4;i++)
171 
                edge[i] = r[i-3 ? i+1 : 0] - r[i];
 170 
                edge[i] = r[i-3 ? i+1 : 0] - r[i];
172 
 
 171 
 
173 
        for(int i=0;i<4;i++)
 172 
        for(int i=0;i<4;i++)
174 
                angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
 173 
                angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
175 
};
 174 
};
176 
 
 175 
 
177 
CPlane4::CPlane4(TVector3 *vr)
 176 
CPlane4::CPlane4(TVector3 *vr)
178 
{
 177 
{
179 
  double x1,y1,z1, x2,y2,z2, x3,y3,z3;
 178 
  double x1,y1,z1, x2,y2,z2, x3,y3,z3;
180 
       
 179 
       
181 
        x1 = vr[0].x(); y1 = vr[0].y(); z1 = vr[0].z();
 180 
        x1 = vr[0].x(); y1 = vr[0].y(); z1 = vr[0].z();
Line 183... Line 182...
183 
        x3 = vr[2].x(); y3 = vr[2].y(); z3 = vr[2].z();
 182 
        x3 = vr[2].x(); y3 = vr[2].y(); z3 = vr[2].z();
184 
 
 183 
 
185 
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
 184 
        A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
186 
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
 185 
        B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
187 
        C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
 186 
        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);
 187 
        D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
189 
 
 188 
 
190 
        r[0] = vr[0]; r[1] = vr[1]; r[2] = vr[2]; r[3] = vr[3];
 189 
        r[0] = vr[0]; r[1] = vr[1]; r[2] = vr[2]; r[3] = vr[3];
191 
        n.SetXYZ(A, B, C);
 190 
        n.SetXYZ(A, B, C);
192 
        n = n.Unit();
 191 
        n = n.Unit();
193 
 
 192 
 
194 
        for(int i=0;i<4;i++)
 193 
        for(int i=0;i<4;i++)
Line 240... Line 239...
240 
int CPlane4::IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double angle)
 239 
int CPlane4::IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double angle)
241 
{
 240 
{
242 
  double angle_ve1, angle_ve2;
 241 
  double angle_ve1, angle_ve2;
243 
       
 242 
       
244 
        if(dbg) printf("--- CPlane4::IsInTri ---\n");
 243 
        if(dbg) printf("--- CPlane4::IsInTri ---\n");
245 
 
 244 
 
246 
        angle_ve1 = TMath::ACos(/*TMath::Abs*/( (e1.Unit()) * (vec.Unit()) ));
 245 
        angle_ve1 = TMath::ACos(/*TMath::Abs*/( (e1.Unit()) * (vec.Unit()) ));
247 
        angle_ve2 = TMath::ACos(/*TMath::Abs*/( (e2.Unit()) * (vec.Unit()) ));
 246 
        angle_ve2 = TMath::ACos(/*TMath::Abs*/( (e2.Unit()) * (vec.Unit()) ));
248 
 
 247 
 
249 
        if(dbg)
 248 
        if(dbg)
250 
        {
 249 
        {
Line 276... Line 275...
276 
 
 275 
 
277 
        return 1;
 276 
        return 1;
278 
}
 277 
}
279 
//-----------------------------------------------------------------------------
 278 
//-----------------------------------------------------------------------------
280 
int CPlane4::TestIntersection(CRay in)
 279 
int CPlane4::TestIntersection(CRay in)
281 
{
 280 
{
282 
        TVector3 tmp;
 281 
        TVector3 tmp;
283 
 
 282 
 
284 
        if( GetIntersection(&tmp, in) )
 283 
        if( GetIntersection(&tmp, in) )
285 
                if( IsVectorIn(tmp) )
 284 
                if( IsVectorIn(tmp) )
286 
                        return 1;
 285 
                        return 1;
287 
       
 286 
       
Line 294... Line 293...
294 
 
 293 
 
295 
        if( GetIntersection(&tmp, in) )
 294 
        if( GetIntersection(&tmp, in) )
296 
                if( IsVectorIn(tmp) ) {
 295 
                if( IsVectorIn(tmp) ) {
297 
                        *vec = tmp;
 296 
                        *vec = tmp;
298 
                        return 1;
 297 
                        return 1;
299 
                }
 298 
                }
300 
       
 299 
       
301 
        return 0;
 300 
        return 0;
302 
}
 301 
}
303 
//-----------------------------------------------------------------------------
 302 
//-----------------------------------------------------------------------------
304 
void CPlane4::Print()
 303 
void CPlane4::Print()
305 
{
 304 
{
Line 446... Line 445...
446 
                                cosTtotal = TMath::Sqrt( 1 - TMath::Power(N1_N2(-sign_n), 2) );
 445 
                                cosTtotal = TMath::Sqrt( 1 - TMath::Power(N1_N2(-sign_n), 2) );
447 
                        else
 446 
                        else
448 
                                cosTtotal = 0.0;
 447 
                                cosTtotal = 0.0;
449 
                       
 448 
                       
450 
                        if(dbg) printf("  cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
 449 
                        if(dbg) printf("  cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
- 
 
 450 
                        // reflection dependance on polarization missing
- 
 
 451 
                        // reflection hardcoded to 0.96
- 
 
 452 
                        p_ref = rand.Uniform(0.0, 1.0);
- 
 
 453 
                        if (dbg) printf("   reflection probability = %f\n", p_ref);    
451 
                       
 454 
                       
452 
                        p_ref = rand.Uniform(0.0, 1.0);
 455 
                        // Total reflection
453 
                       
 456 
                        /*
454 
                        if( (cosTi <= cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
 457 
                        if( (cosTi <= cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
455 
                                if(dbg) printf("  TOTAL\n");
 458 
                                if(dbg) printf("  TOTAL\n");
456 
                                transmit = in.GetN() + sign_n*2*cosTi*n;
 459 
                                transmit = in.GetN() + sign_n*2*cosTi*n;
457 
                               
 460 
                               
458 
                                if(dbg) {
 461 
                                if(dbg) {
Line 461... Line 464...
461 
                                }      
 464 
                                }      
462 
                                out->Set(intersect, transmit);
 465 
                                out->Set(intersect, transmit);
463 
                               
 466 
                               
464 
                                pol_t = -in.GetP() + sign_n*2*cosTi*n;
 467 
                                pol_t = -in.GetP() + sign_n*2*cosTi*n;
465 
                                out->SetPolarization(pol_t);
 468 
                                out->SetPolarization(pol_t);
466 
                                return 2;
 469 
                                return REFLECTION;
467 
                        } else { // ni tot. odboja
 470 
                        } else { */
- 
 
 471 
                          // reflection or refraction according to Fresnel equations
468 
                                if(dbg) printf("  REFRACTION\n");
 472 
                                if(dbg) printf("  REFRACTION\n");
469 
                                if(dbg) printf("  N1_N2(sign_n) = %lf\n", N1_N2(sign_n));      
 473 
                                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)));                          
 474 
                                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);
 475 
                                if(dbg) printf("  cosTt = %lf (Tt = %lf) \n", cosTt, TMath::ACos(cosTt)*DEGREE);
472 
                                                                               
 476 
                                                                               
Line 476... Line 480...
476 
                                if(dbg) {
 480 
                                if(dbg) {
477 
                                        cosTN = transmit.Unit() * n;
 481 
                                        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); 
 482 
                                        printf("  cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE); 
479 
                                }
 483 
                                }
480 
                                //if(cosTi<=cosTtotal) cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));
 484 
                                //if(cosTi<=cosTtotal) cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));
481 
                                if(fresnel) {
 485 
                                //if(fresnel) {
482 
                                        r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse
 486 
                                        r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse
483 
                                        r_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel
 487 
                                        r_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel
484 
                                       
 488 
                                       
485 
                                        if(dbg) printf("  r_te = %lf, r_tm = %lf\n", r_te, r_tm);
 489 
                                        if(dbg) printf("  r_te = %lf, r_tm = %lf\n", r_te, r_tm);
486 
                                       
 490 
                                       
Line 495... Line 499...
495 
 
 499 
 
496 
                                        R_te = TMath::Power(r_te, 2);
 500 
                                        R_te = TMath::Power(r_te, 2);
497 
                                        R_tm = TMath::Power(r_tm, 2);
 501 
                                        R_tm = TMath::Power(r_tm, 2);
498 
                                        R_f = a_te*a_te*R_te + a_tm*a_tm*R_tm;
 502 
                                        R_f = a_te*a_te*R_te + a_tm*a_tm*R_tm;
499 
                                       
 503 
                                       
500 
                                        if(dbg) printf("  R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
 504 
                                        if (dbg) printf("  R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
501 
                                }
 505 
                                //}
502 
                               
 506 
                                       
503 
                                //p_ref = rand.Uniform(0.0, 1.0);      
- 
 
504 
                                if(p_ref >= R_f) { // se lomi
 507 
                                if(p_ref >= R_f) { // se lomi
- 
 
 508 
                                  if (dbg) printf("   SURFACE REFRACTED. Return.\n");
505 
                                        out->Set(intersect, transmit);
 509 
                                        out->Set(intersect, transmit);
506 
                                        out->SetPolarization(pol_t);
 510 
                                        out->SetPolarization(pol_t);
507 
                                        return 1;
 511 
                                        return REFRACTION;
508 
                                } else { // se odbije
 512 
                                } else { // se odbije
- 
 
 513 
                                  if (dbg) printf("   SURFACE REFLECTED. p_ref=%f, R_f=%f\n", p_ref, R_f);
509 
                                        transmit = in.GetN() + sign_n*2*cosTi*n;
 514 
                                        transmit = in.GetN() + sign_n*2*cosTi*n;
510 
                                        out->Set(intersect, transmit);
 515 
                                        out->Set(intersect, transmit);
511 
                                        pol_t = -in.GetP() + sign_n*2*cosTi*n;
 516 
                                        pol_t = -in.GetP() + sign_n*2*cosTi*n;
512 
                                        out->SetPolarization(pol_t);
 517 
                                        out->SetPolarization(pol_t);
513 
                                        return 2;
 518 
                                        return REFLECTION;
- 
 
 519 
                                }      
514 
                                }                              
 520 
                                                       
515 
                        }
 521 
                        //}
516 
                        break;
 522 
                        break;
- 
 
 523 
                       
517 
                // ----------------------------------------------------------------------------
 524 
                // ----------------------------------------------------------------------------
518 
                // --------------- reflection at "reflection" probability ---------------------
 525 
                // --------------- reflection at "reflection" probability ---------------------
519 
                // ----------------------------------------------------------------------------
 526 
                // ----------------------------------------------------------------------------
520 
                case SURF_REFLE:
 527 
                case SURF_REFLE:
521 
                        p_ref = rand.Uniform(0.0, 1.0);
 528 
                        p_ref = rand.Uniform(0.0, 1.0);
522 
                        if(p_ref < reflection) { // se odbije
 529 
                        if(p_ref < reflection) { // se odbije
523 
                                cosTi = in.GetN() * n;
 530 
                                cosTi = in.GetN() * n;
524 
                                transmit = in.GetN() - 2*cosTi*n;
 531 
                                transmit = in.GetN() - 2*cosTi*n;
525 
                                out->Set(intersect, transmit);
 532 
                                out->Set(intersect, transmit);
526 
                                return 2; //sdhfvjhsdbfjhsdbcvjhsb
 533 
                                return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
527 
                        } else { // se ne odbije
 534 
                        } else { // se ne odbije
528 
                                transmit = in.GetN();
 535 
                                transmit = in.GetN();
529 
                                out->Set(intersect, transmit);         
 536 
                                out->Set(intersect, transmit);         
530 
                                return -2;
 537 
                                return ABSORBED;
531 
                        }                                              
 538 
                        }                                              
532 
                        break;
 539 
                        break;
533 
                       
 540 
                       
534 
                // total reflection from n1 to n2 with R probbability
 541 
                // total reflection from n1 to n2 with R probbability
535 
                case SURF_IMPER:
 542 
                case SURF_IMPER:
Line 540... Line 547...
540 
                                        transmit = in.GetN() - 2*cosTi*n;
 547 
                                        transmit = in.GetN() - 2*cosTi*n;
541 
                                        out->Set(intersect, transmit);
 548 
                                        out->Set(intersect, transmit);
542 
                                } else { // ni tot. odboja
 549 
                                } else { // ni tot. odboja
543 
                                        transmit = in.GetN();
 550 
                                        transmit = in.GetN();
544 
                                        out->Set(intersect, transmit);
 551 
                                        out->Set(intersect, transmit);
545 
                                        return -2;
 552 
                                        return ABSORBED;
546 
                                }
 553 
                                }
547 
                        } else { // se ne odbije
 554 
                        } else { // se ne odbije
548 
                                transmit = in.GetN();
 555 
                                transmit = in.GetN();
549 
                                out->Set(intersect, transmit);                 
 556 
                                out->Set(intersect, transmit);                 
550 
                                return -2;
 557 
                                return ABSORBED;
551 
                        }                                              
 558 
                        }                                              
552 
                        break;
 559 
                        break;
553 
 
 560 
 
554 
                default:
 561 
                default:
555 
                        *out = in;
 562 
                        *out = in;
556 
                        break;
 563 
                        break;
557 
        }
 564 
        }
558 
 
 565 
 
559 
        return 1;
 566 
        return REFRACTION;
560 
}
 567 
}
561 
//=================================================================================
 568 
//=================================================================================
562 
 
 569 
 
563 
 
 570 
 
564 
//=================================================================================
 571 
//=================================================================================
Line 603... Line 610...
603 
        if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1);
 610 
        if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1);
604 
       
 611 
       
605 
        hfate = (TH1F*)gROOT->FindObject("hfate"); if(hfate) delete hfate;
 612 
        hfate = (TH1F*)gROOT->FindObject("hfate"); if(hfate) delete hfate;
606 
        hfate = new TH1F("hfate", "Ray fate", 8, -3.5, 4.5);
 613 
        hfate = new TH1F("hfate", "Ray fate", 8, -3.5, 4.5);
607 
        (hfate->GetXaxis())->SetBinLabel(1, "Back Ref");
 614 
        (hfate->GetXaxis())->SetBinLabel(1, "Back Ref");
608 
        (hfate->GetXaxis())->SetBinLabel(2, "Ref Loss");
 615 
        (hfate->GetXaxis())->SetBinLabel(2, "No Ref");
609 
        (hfate->GetXaxis())->SetBinLabel(3, "No Tot Refl");
 616 
        (hfate->GetXaxis())->SetBinLabel(3, "Refrac");
610 
        (hfate->GetXaxis())->SetBinLabel(4, "LG Miss");
 617 
        (hfate->GetXaxis())->SetBinLabel(4, "LG Miss");
611 
        (hfate->GetXaxis())->SetBinLabel(5, "Exit");
 618 
        (hfate->GetXaxis())->SetBinLabel(5, "Exit");
612 
        (hfate->GetXaxis())->SetBinLabel(6, "Enter");
 619 
        (hfate->GetXaxis())->SetBinLabel(6, "Enter");
613 
        (hfate->GetXaxis())->SetBinLabel(7, "Rays");
 620 
        (hfate->GetXaxis())->SetBinLabel(7, "Rays");
614 
        (hfate->GetXaxis())->SetBinLabel(8, "Absorb");
 621 
        (hfate->GetXaxis())->SetBinLabel(8, "Absorb");
615 
       
 622 
       
616 
        hnodb_all = (TH1F*)gROOT->FindObject("hnodb_all"); if(hnodb_all) delete hnodb_all;
 623 
        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);
 624 
        hnodb_all = new TH1F("hnodb_all", "N reflected", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
618 
       
 625 
       
619 
        hnodb_exit = (TH1F*)gROOT->FindObject("hnodb_exit"); if(hnodb_exit) delete hnodb_exit;
 626 
        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);
 627 
        hnodb_exit = new TH1F("hnodb_exit", "N reflected and exit", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
621 
       
 628 
       
622 
        int nBins = nch + 1;
 629 
        int nBins = nch + 1;
623 
        hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin;
 630 
        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);
 631 
        hin = new TH2F("hin", "Guide entrance window", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0);
625 
       
 632 
       
Line 642... Line 649...
642 
{
 649 
{
643 
  if (dbg) printf("--- GUIDE::PropagateRay ---\n");
 650 
  if (dbg) printf("--- GUIDE::PropagateRay ---\n");
644 
  CRay ray0;
 651 
  CRay ray0;
645 
  CRay ray1;
 652 
  CRay ray1;
646 
  TVector3 vec0, vec1;
 653 
  TVector3 vec0, vec1;
647 
  //int fate = 11;
 - 
 
648 
  int inters_i = 0;
 654 
  int inters_i = 0;
649 
               
 655 
               
650 
        ray0 = in;
 656 
        ray0 = in;
651 
        int n_odb = 0;
 657 
        int n_odb = 0;
652 
        int last_hit = 0;
 658 
        int last_hit = 0;
653 
        int propagation = 0;
 659 
        int propagation = 0;
654 
        if( !(s_side[0]->PropagateRay(ray0, &ray1, &vec1)) ) {
 660 
        int result = s_side[0]->PropagateRay(ray0, &ray1, &vec1);
- 
 
 661 
        if( !(result) ) {
655 
                // ce -NI- presecisca z vstopno
 662 
                // ce -NI- presecisca z vstopno
656 
                if (dbg) printf("  GUIDE: missed the light guide\n");
 663 
                if (dbg) printf("  GUIDE: missed the light guide\n");
657 
                fate = missed;
 664 
                fate = missed;
658 
                hfate->Fill(0);
 665 
                //hfate->Fill(0);
- 
 
 666 
        } else if(result == REFLECTION) {
- 
 
 667 
          if (dbg) printf(" REFLECTED on the entry surface!\n");
- 
 
 668 
          fate = backreflected;
- 
 
 669 
          //hfate->Fill(-3);
659 
        } else {
 670 
        } else {
660 
          if (dbg) printf("  GUIDE: ray entered\n");
 671 
          if (dbg) printf("  GUIDE: ray entered\n");
661 
                points[0] = ray1.GetR();
 672 
                points[0] = ray1.GetR();
662 
                hfate->Fill(2); // enter
 673 
                hfate->Fill(2); // enter
663 
                hin->Fill(vec1.y(), vec1.z());
 674 
                hin->Fill(vec1.y(), vec1.z());
664 
                if (dbg) printf("  GUIDE: n_odb = %d\n", n_odb);
 675 
                if (dbg) printf("  GUIDE: n_odb = %d\n", n_odb);
- 
 
 676 
               
665 
                while (n_odb++ < VODNIK_MAX_ODB) {
 677 
                while (n_odb++ < MAX_REFLECTIONS) {
666 
                  if (dbg) printf("  GUIDE: Boundary test: %d\n",n_odb);
 678 
                  if (dbg) printf("  GUIDE: Boundary test: %d\n",n_odb);
667 
                        ray0 = ray1;
 679 
                        ray0 = ray1;
668 
                        vec0 = vec1;
 680 
                        vec0 = vec1;
669 
                        propagation = 11;
 681 
                        propagation = 11;
670 
                        for(inters_i=0; inters_i<6; inters_i++) {
 682 
                        for(inters_i=0; inters_i<6; inters_i++) {
Line 678... Line 690...
678 
                                }  
 690 
                                }  
679 
                        }                      
 691 
                        }                      
680 
                        points[n_odb] = vec1;
 692 
                        points[n_odb] = vec1;
681 
                        if(inters_i == 0) {
 693 
                        if(inters_i == 0) {
682 
                          fate = backreflected;
 694 
                          fate = backreflected;
683 
                          hfate->Fill(backreflected);
 695 
                          //hfate->Fill(backreflected);
684 
                          break;
 696 
                          break;
685 
                        } // backreflection
 697 
                        } // backreflection
686 
                                               
 698 
                       
- 
 
 699 
                        // the passage is possible, test propagation                   
687 
                        propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
 700 
                        propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
- 
 
 701 
                       
688 
                        if (dbg) printf("  GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
 702 
                        if (dbg) printf("  GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
- 
 
 703 
                       
- 
 
 704 
                        if(propagation == REFRACTION) {
- 
 
 705 
                          fate = refracted;
- 
 
 706 
                          n_odb++;
- 
 
 707 
                          points[n_odb] = vec1;
- 
 
 708 
                          ray0 = ray1;
- 
 
 709 
                          break;
- 
 
 710 
                        } // no total reflection when should be
- 
 
 711 
                        if(propagation == ABSORBED) {
- 
 
 712 
                          fate = noreflection;
- 
 
 713 
                          break;
- 
 
 714 
                        } //refraction due to finite reflectivity
- 
 
 715 
                       
689 
                        if(inters_i == 5) { // successfull exit
 716 
                        if(inters_i == 5) { // successfull exit
690 
                        // check on which side the vector is?
 717 
                        // check on which side the vector is?
691 
                          TVector3 ray = ray1.GetN();
 718 
                          TVector3 ray = ray1.GetN();
692 
                          TVector3 exitNormal = s_side[5]->GetN();
 719 
                          TVector3 exitNormal = s_side[5]->GetN();
693 
                          //printf("theta(ray) = %lf, theta(normal5) = %lf ", ray.Theta()*DEGREE, exitNormal.Theta()*DEGREE);
 720 
                          //printf("theta(ray) = %lf, theta(normal5) = %lf ", ray.Theta()*DEGREE, exitNormal.Theta()*DEGREE);
Line 699... Line 726...
699 
                            n_odb++;
 726 
                            n_odb++;
700 
                            points[n_odb] = vec1;
 727 
                            points[n_odb] = vec1;
701 
                            ray0 = ray1;
 728 
                            ray0 = ray1;
702 
                            break;
 729 
                            break;
703 
                          }
 730 
                          }
704 
                                fate =  hit;
 731 
                                fate =  hitExit;
705 
                                hout->Fill(vec1.y(), vec1.z());
 732 
                                hout->Fill(vec1.y(), vec1.z());
706 
                                hnodb_exit->Fill(n_odb-1);
 733 
                                hnodb_exit->Fill(n_odb-1);
707 
                                n_odb++;
 734 
                                n_odb++;
708 
                                points[n_odb] = vec1;
 735 
                                points[n_odb] = vec1;
709 
                                ray0 = ray1;
 736 
                                ray0 = ray1;
710 
                                break;
 737 
                                break;
711 
                        }
 738 
                        }
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;
 739 
                        last_hit = inters_i;
725 
                }      
 740 
                }      
726 
        }
 741 
        }
727 
       
 742 
       
728 
        //--- material absorption ---
 743 
        //--- material absorption ---
Line 741... Line 756...
741 
                if(p_abs > T_abs) fate = absorbed; // absorption
 756 
                if(p_abs > T_abs) fate = absorbed; // absorption
742 
        }
 757 
        }
743 
        //--- material absorption ---
 758 
        //--- material absorption ---
744 
       
 759 
       
745 
        hfate->Fill(fate);
 760 
        hfate->Fill(fate);
746 
        hfate->Fill(3);
 761 
        hfate->Fill(rays);
747 
        hnodb_all->Fill(n_odb-2);
 762 
        hnodb_all->Fill(n_odb-2);
748 
        *n_points = n_odb+1;
 763 
        *n_points = n_odb+1;
749 
        *out = ray0;
 764 
        *out = ray0;
750 
        return fate;
 765 
        return fate;
751 
}
 766 
}
Line 762... Line 777...
762 
//-----------------------------------------------------------------------------
 777 
//-----------------------------------------------------------------------------
763 
void Guide::DrawSkel(int color, int width)
 778 
void Guide::DrawSkel(int color, int width)
764 
{
 779 
{
765 
        TPolyLine3D *line3d = new TPolyLine3D(2);
 780 
        TPolyLine3D *line3d = new TPolyLine3D(2);
766 
        line3d->SetLineWidth(width); line3d->SetLineColor(color);
 781 
        line3d->SetLineWidth(width); line3d->SetLineColor(color);
767 
 
 782 
 
768 
        for(int i=0; i<4; i++) {
 783 
        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());
 784 
                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());
 785 
                line3d->SetPoint(1, vodnik_edge[i+4].x(), vodnik_edge[i+4].y(), vodnik_edge[i+4].z());
771 
                line3d->DrawClone();
 786 
                line3d->DrawClone();
772 
        }
 787 
        }
773 
}
 788 
}
774 
//=================================================================================
 789 
//=================================================================================
775 
 
 790 
 
776 
//=================================================================================
 791 
//=================================================================================
777 
int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
 792 
int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
778 
{
 793 
{
779 
        double num, den; //stevec, imenovalec
 794 
        double num, den; //stevec, imenovalec
780 
        double t;
 795 
        double t;
781 
        TVector3 tmp;
 796 
        TVector3 tmp;
782 
 
 797 
 
783 
        if(dbg) printf("---> CPlaneR::TestIntersection <---\n");
 798 
        if(dbg) printf("---> CPlaneR::TestIntersection <---\n");
784 
        if(dbg) {printf("c = "); printv(center); printf(" | n = "); printv(n); printf("\n");}
 799 
        if(dbg) {printf("c = "); printv(center); printf(" | n = "); printv(n); printf("\n");}
785 
       
 800 
       
786 
        double D = - n*center; 
 801 
        double D = - n*center; 
787 
        num = n*ray.GetR() + D;
 802 
        num = n*ray.GetR() + D;
788 
        den = n*ray.GetN();
 803 
        den = n*ray.GetN();
789 
       
 804 
       
790 
        if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den);
 805 
        if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den);
Line 856... Line 871...
856 
        //window_R( parameters.getB() ),
 871 
        //window_R( parameters.getB() ),
857 
        //window_d(0),
 872 
        //window_d(0),
858 
  guide(new Guide(center0, parameters)),
 873 
  guide(new Guide(center0, parameters)),
859 
  plate(new Plate(parameters)),
 874 
  plate(new Plate(parameters)),
860 
  _plateWidth(parameters.getPlateWidth()),
 875 
  _plateWidth(parameters.getPlateWidth()),
861 
  _plateOn(parameters.getPlateOn())
 876 
  _plateOn(parameters.getPlateOn()),
- 
 
 877 
  offsetY(parameters.getOffsetY()),
- 
 
 878 
  offsetZ(parameters.getOffsetZ())
862 
  {
 879 
  {
863 
//  };
 880 
//  };
864 
 
 881 
 
865 
//-----------------------------------------------------------------------------
 882 
//-----------------------------------------------------------------------------
866 
//void CDetector::Init()
 883 
//void CDetector::Init()
Line 906... Line 923...
906 
        plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size);
 923 
        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);
 924 
        plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size);
908 
        active = new CPlane4(plane_v);
 925 
        active = new CPlane4(plane_v);
909 
       
 926 
       
910 
        hactive = (TH2F*)gROOT->FindObject("hactive"); if(hactive) delete hactive;
 927 
        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);
 928 
        //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 
       
 - 
 
- 
 
 929 
        hactive = new TH2F("hactive", "Active area hits", nBins, y_gap - p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ);
913 
       
 930 
       
914 
        p_size = b/2.0;
 931 
        p_size = b/2.0;
915 
        //p_size = 2.5;
 932 
        //p_size = 2.5;
916 
        //p_size = M*0.6;
 933 
        //p_size = M*0.6;
917 
        hlaser = (TH2F*)gROOT->FindObject("hlaser"); if(hlaser) delete hlaser;
 934 
        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);
 935 
        hlaser = new TH2F("hlaser", ";x [mm]; y [mm]", nBins, -p_size+offsetY, p_size+offsetY, nBins, -p_size+offsetZ, p_size+offsetZ);
919 
       
 936 
       
920 
       
 937 
       
921 
        p_size = 2*b/2.0;
 938 
        p_size = 1.4*b/2.0;
922 
        plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size);
 939 
        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);
 940 
        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);
 941 
        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);
 942 
        plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size);
926 
        detector = new CPlane4(plane_v);
 943 
        detector = new CPlane4(plane_v);
927 
       
 944 
       
928 
        hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector;
 945 
        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);
 946 
        //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 
       
 - 
 
- 
 
 947 
        hdetector = new 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);
931 
       
 948 
       
932 
        /*
 949 
        /*
933 
        window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R);  
 950 
        window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R);  
934 
       
 951 
       
935 
        p_size = M*a;
 952 
        p_size = M*a;
Line 944... Line 961...
944 
        */
 961 
        */
945 
        p_size = b/2.0;
 962 
        p_size = b/2.0;
946 
        histoPlate = (TH2F*)gROOT->FindObject("histoPlate"); if(histoPlate) delete histoPlate;
 963 
        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);
 964 
        histoPlate = new TH2F("histoPlate", "Hits on glass plate", nBins, -p_size, +p_size, nBins, -p_size, +p_size);
948 
}
 965 
}
- 
 
 966 
 
949 
//-----------------------------------------------------------------------------
 967 
//-----------------------------------------------------------------------------
950 
// vrne 1 ce je zadel aktvino povrsino
 968 
// vrne 1 ce je zadel aktvino povrsino
951 
// vrne <1 ce jo zgresi
 969 
// vrne <1 ce jo zgresi
952 
int CDetector::Propagate(CRay in, CRay *out, int draw)
 970 
int CDetector::Propagate(CRay in, CRay *out, int draw)
- 
 
 971 
 
953 
// Sledi zarku skozi vodnik. Vrne:                                            
972 
// Sledi zarku skozi vodnik. Vrne:                                            
954 
//  0, ce zgresi vstopno ploskev MISSED                                              
 973 
//  0, ce zgresi vstopno ploskev MISSED                                              
955 
//  1, ce zadane izstopno ploskev HIT                                            
974 
//  1, ce zadane izstopno ploskev HIT                                            
956 
// -1, ce se v vodniku ne odbije totalno REFRACTED
 975 
// -1, ce se v vodniku ne odbije totalno REFRACTED
957 
//  2, enter the light guide, bin 2 of hfate EXIT                                    
976 
//  2, enter the light guide, bin 2 of hfate EXIT                                    
Line 966... Line 985...
966 
        CRay *rayout = new CRay;
 985 
        CRay *rayout = new CRay;
967 
 
 986 
 
968 
        const int max_n_points = guide->GetMAXODB() + 2;
 987 
        const int max_n_points = guide->GetMAXODB() + 2;
969 
        TVector3 pointsPlate[max_n_points];
 988 
        TVector3 pointsPlate[max_n_points];
970 
        //TVector3 intersection;
 989 
        //TVector3 intersection;
971 
        int fatePlate;
 990 
        Fate fatePlate;
972 
        int nPointsPlate;
 991 
        int nPointsPlate;
973 
        TPolyLine3D *line3d = new TPolyLine3D(2);
 992 
        TPolyLine3D *line3d = new TPolyLine3D(2);
974 
        line3d->SetLineWidth(1);
 993 
        line3d->SetLineWidth(1);
975 
        line3d->SetLineColor(4);
 994 
        line3d->SetLineColor(4);
976 
       
 995 
       
Line 981... Line 1000...
981 
           
 1000 
           
982 
          fatePlate = plate->propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
 1001 
          fatePlate = plate->propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
983 
          if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
 1002 
          if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
984 
         
 1003 
         
985 
          if(draw) {
 1004 
          if(draw) {
986 
            if(fatePlate == 0) {
 1005 
            if(fatePlate == missed) {
987 
              rayout->SetColor(col_in);
 1006 
              rayout->SetColor(col_in);
988 
              rayout->DrawS(center.x() - _plateWidth, -10.0);
 1007 
              rayout->DrawS(center.x() - _plateWidth, -10.0);
989 
              }
 1008 
              }
990 
              else {
 1009 
              else {
991 
                int p_i;
 1010 
                int p_i;
992 
                for(p_i = 0; p_i < nPointsPlate-1; p_i++) {
 1011 
                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());
 1012 
                                          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());
 1013 
                                          line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z());
995 
                                          line3d->DrawClone();
 1014 
                                          line3d->DrawClone();
996 
                                  }
 1015 
                                  }
997 
                                  if(fatePlate != -2) {
 1016 
                                  if(fatePlate != noreflection) {
998 
                                          //rayout->SetColor(7);
 1017 
                                          //rayout->SetColor(7);
999 
                                          rayout->DrawS(pointsPlate[p_i].x(), -0.1);
 1018 
                                          rayout->DrawS(pointsPlate[p_i].x(), -0.1);
1000 
                                  }
 1019 
                                  }
1001 
                          }
 1020 
                          }
1002 
            }
 1021 
            }
1003 
 
 1022 
 
1004 
          if(! (fatePlate == 1 or fatePlate == -1) ) {
 1023 
          if(! (fatePlate == hitExit or fatePlate == refracted) ) {
1005 
                       
 1024 
                       
1006 
                        cout<<"CDetector::propagate Simulated ray missed the entry surface!"<<endl;
 1025 
                  if (dbg) printf("CDetector::propagate Simulated ray missed the entry surface!\n");
1007 
                        return fatePlate;
 1026 
                        return fatePlate;
1008 
                }
 1027 
                }
- 
 
 1028 
                // missing: if refracted at plate sides
- 
 
 1029 
                //if (fatePlate == refracted) return fatePlate;
1009 
                histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
 1030 
                histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
1010 
         }
 1031 
         }
1011 
         else {
 1032 
         else {
1012 
          rayout = rayin;
 1033 
          rayout = rayin;
1013 
          if(draw) rayout->DrawS(center.x(), -10.0);
 1034 
          if(draw) rayout->DrawS(center.x(), -10.0);
Line 1017... Line 1038...
1017 
        // Draw the light guide and propagate the ray through
 1038 
        // Draw the light guide and propagate the ray through
1018 
       
 1039 
       
1019 
        //const int max_n_points = guide->GetMAXODB() + 2;
 1040 
        //const int max_n_points = guide->GetMAXODB() + 2;
1020 
        TVector3 points[max_n_points];
 1041 
        TVector3 points[max_n_points];
1021 
        TVector3 presecisce;
 1042 
        TVector3 presecisce;
1022 
        //int fate;
 1043 
 
1023 
        int n_points;
 1044 
        int n_points;
1024 
        int fate_glass;
 1045 
        int fate_glass;
1025 
        CRay *ray0 = new CRay(*rayout);
 1046 
        CRay *ray0 = new CRay(*rayout);
1026 
        // delete rayout; -> creates dangling reference when tries to delete ray0!
 1047 
        // delete rayout; -> creates dangling reference when tries to delete ray0!
1027 
        //delete rayin;
 1048 
        delete rayin;
1028 
        CRay *ray1 = new CRay;
 1049 
        CRay *ray1 = new CRay;
1029 
         
 1050 
         
1030 
        fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
 1051 
        fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
- 
 
 1052 
        if (dbg) {
- 
 
 1053 
          if (fate == backreflected) printf("DETECTOR::backreflected\n");
- 
 
 1054 
          }
1031 
       
 1055 
       
1032 
        line3d->SetLineColor(col_lg);  
 1056 
        line3d->SetLineColor(col_lg);  
1033 
        int p_i;
 1057 
        int p_i;
1034 
        if(guide_on) {
 1058 
        if(guide_on) {
1035 
                if(draw) {
 1059 
                if(draw) {
Line 1055... Line 1079...
1055 
                                }
 1079 
                                }
1056 
                        }
 1080 
                        }
1057 
                }
 1081 
                }
1058 
               
 1082 
               
1059 
               
 1083 
               
1060 
                if(! (fate == hit or fate == refracted) ) {
 1084 
                if(! (fate == hitExit or fate == refracted) ) {
1061 
                if (dbg) printf("Detector: fate != hit, refracted\n");
 1085 
                  if (dbg) printf("Detector: fate != hit, refracted\n");
1062 
                        *out = *ray1;
 1086 
                        *out = *ray1;
1063 
                        return fate;
 1087 
                        return fate;
1064 
                }
 1088 
                }
1065 
        } else {
 1089 
        } else {
1066 
          if (dbg) printf("Detector: fate = hit or refracted");
 1090 
          if (dbg) printf("Detector: fate = hit or refracted");
Line 1089... Line 1113...
1089 
                        fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
 1113 
                        fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
1090 
                        if(fate_glass == 1) {
 1114 
                        if(fate_glass == 1) {
1091 
                                hglass->Fill(presecisce.y(), presecisce.z());
 1115 
                                hglass->Fill(presecisce.y(), presecisce.z());
1092 
                                if(draw) ray1->DrawS(presecisce.x(), 10.0);
 1116 
                                if(draw) ray1->DrawS(presecisce.x(), 10.0);
1093 
                                if(active->TestIntersection(&presecisce, *ray1)) {
 1117 
                                if(active->TestIntersection(&presecisce, *ray1)) {
1094 
                                        fate = hit;
 1118 
                                        fate = hitExit;
1095 
                                        hactive->Fill(presecisce.y(), presecisce.z());
 1119 
                                        hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
1096 
                                        hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
 1120 
                                        hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
1097 
                                }
 1121 
                                }
1098 
                                if(detector->TestIntersection(&presecisce, *ray1))
 1122 
                                if(detector->TestIntersection(&presecisce, *ray1))
1099 
                                        hdetector->Fill(presecisce.y(), presecisce.z());
 1123 
                                        hdetector->Fill(presecisce.y(), presecisce.z());
1100 
                        } else if(fate_glass == 2) {
 1124 
                        } else if(fate_glass == 2) {
1101 
                                if(draw) ray1->DrawS(presecisce.x(), 10.0);
 1125 
                                if(draw) ray1->DrawS(presecisce.x(), 10.0);
1102 
                        }
 1126 
                        }
1103 
                } else {
 1127 
                } else {
- 
 
 1128 
                  // Main test: ray and SiPM surface
1104 
                        if(active->TestIntersection(&presecisce, *ray1)) {
 1129 
                        if(active->TestIntersection(&presecisce, *ray1)) {
1105 
                                fate = hit;
 1130 
                                fate = hitExit;
1106 
                                hactive->Fill(presecisce.y(), presecisce.z());
 1131 
                                hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
1107 
                                hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
 1132 
                                hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ);
1108 
                        }
 1133 
                        }
- 
 
 1134 
                        // If it is on the same plane as SiPM
1109 
                        if(detector->TestIntersection(&presecisce, *ray1))
 1135 
                        if(detector->TestIntersection(&presecisce, *ray1))
1110 
                                hdetector->Fill(presecisce.y(), presecisce.z());
 1136 
                                hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
1111 
                }
 1137 
                }
1112 
        //} else {
 1138 
        //} else {
1113 
                //if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
 1139 
                //if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
1114 
        //}
 1140 
        //}
1115 
                               
 1141 
                               
1116 
        *out = *ray1;
 1142 
        *out = *ray1;
1117 
        //delete ray0;
 1143 
        delete ray0;
1118 
        //delete ray1;
 1144 
        delete ray1;
- 
 
 1145 
        delete rayout;
1119 
        return fate;
 1146 
        return fate;
1120 
}
 1147 
}
1121 
//-----------------------------------------------------------------------------
 1148 
//-----------------------------------------------------------------------------
1122 
void CDetector::Draw(int width)
 1149 
void CDetector::Draw(int width)
1123 
{
 1150 
{
Line 1186... Line 1213...
1186 
                line3d.SetPoint(1, plate_edge[i+4].x(), plate_edge[i+4].y(), plate_edge[i+4].z());
 1213 
                line3d.SetPoint(1, plate_edge[i+4].x(), plate_edge[i+4].y(), plate_edge[i+4].z());
1187 
                line3d.DrawClone();
 1214 
                line3d.DrawClone();
1188 
        }
 1215 
        }
1189 
}
 1216 
}
1190 
 
 1217 
 
1191 
int Plate::propagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
 1218 
Fate Plate::propagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
1192 
{
 1219 
{
1193 
  CRay ray0;
 1220 
  CRay ray0;
1194 
  CRay ray1;
 1221 
  CRay ray1;
1195 
  TVector3 vec0, vec1;
 1222 
  TVector3 vec0, vec1;
1196 
  int fate = 11;
 1223 
  Fate fate = enter;
1197 
  int inters_i = 0;
 1224 
  int inters_i = 0;
1198 
               
 1225 
               
1199 
        ray0 = in;
 1226 
        ray0 = in;
1200 
        int n_odb = 0;
 1227 
        int n_odb = 0;
1201 
        int last_hit = 0;
 1228 
        int last_hit = 0;
1202 
        int propagation = 0;
 1229 
        int propagation = 0;
1203 
       
 1230 
       
1204 
        if( !(sides[0]->PropagateRay(ray0, &ray1, &vec1)) ) {
 1231 
        int result = sides[0]->PropagateRay(ray0, &ray1, &vec1);
- 
 
 1232 
        if( !result ) {
1205 
                // ce -NI- presecisca z vstopno
 1233 
                // ce -NI- presecisca z vstopno
1206 
                fate = 0;
 1234 
                fate = missed;
- 
 
 1235 
        } else if(result == REFLECTION) {
- 
 
 1236 
          if (dbg) printf("PLATE: reflected\n");
- 
 
 1237 
          fate = backreflected;
1207 
        } else {
 1238 
        } else {
1208 
                points[0] = ray1.GetR();
 1239 
                points[0] = ray1.GetR();
1209 
                //hfate->Fill(2);
 1240 
                //hfate->Fill(2);
1210 
                //hin->Fill(vec1.y(), vec1.z());
 1241 
                //hin->Fill(vec1.y(), vec1.z());
1211 
                while (n_odb++ < VODNIK_MAX_ODB) {
 1242 
                while (n_odb++ < MAX_REFLECTIONS) {
1212 
                        ray0 = ray1;
 1243 
                        ray0 = ray1;
1213 
                        vec0 = vec1;
 1244 
                        vec0 = vec1;
1214 
                        propagation = 11;
 1245 
                        propagation = 11;
1215 
                        for(inters_i=0; inters_i<6; inters_i++) {
 1246 
                        for(inters_i=0; inters_i<6; inters_i++) {
1216 
                                if( inters_i != last_hit) {
 1247 
                                if( inters_i != last_hit) {
1217 
                                        if( sides[inters_i]->TestIntersection(&vec1, ray1) ) break;
 1248 
                                        if( sides[inters_i]->TestIntersection(&vec1, ray1) ) break;
1218 
                                        }
 1249 
                                        }
1219 
                                }                                                      
 1250 
                                }                                                      
1220 
                        points[n_odb] = vec1;
 1251 
                        points[n_odb] = vec1;
1221 
                        if(inters_i == 0) {
 1252 
                        if(inters_i == 0) {
1222 
                          fate = -3;
 1253 
                          fate = backreflected;
1223 
                          break;} // backreflection
 1254 
                          break;} // backreflection
1224 
                                               
 1255 
                                               
1225 
                        propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
 1256 
                        propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
1226 
                        if(inters_i == 5) { // successfull exit
 1257 
                        if(inters_i == 5) { // successfull exit
1227 
                                fate =  1;
 1258 
                                fate = hitExit;
1228 
                                //hout->Fill(vec1.y(), vec1.z());
1259 
                                //hout->Fill(vec1.y(), vec1.z());
1229 
                                //hnodb_exit->Fill(n_odb-1);
1260 
                                //hnodb_exit->Fill(n_odb-1);
1230 
                                n_odb++;
 1261 
                                n_odb++;
1231 
                                points[n_odb] = vec1;
 1262 
                                points[n_odb] = vec1;
1232 
                                ray0 = ray1;
 1263 
                                ray0 = ray1;
1233 
                                break;
 1264 
                                break;
1234 
                        }
 1265 
                        }
1235 
                        if(propagation == 1) {
 1266 
                        if(propagation == 1) {
1236 
                          fate = -1;
 1267 
                          fate = refracted; //at side?
1237 
                          n_odb++;
 1268 
                          n_odb++;
1238 
                          points[n_odb] = vec1;
 1269 
                          points[n_odb] = vec1;
1239 
                          ray0 = ray1;
 1270 
                          ray0 = ray1;
1240 
                          break;} // no total reflection when should be
 1271 
                          break;} // no total reflection when should be
- 
 
 1272 
                       
- 
 
 1273 
                        if(propagation == -2) {
- 
 
 1274 
                          fate = noreflection;
- 
 
 1275 
                          break;
1241 
                        if(propagation == -2) {fate = -2; break;} // absorption due to finite reflectivity
 1276 
                          } // absorption due to finite reflectivity
1242 
                       
 1277 
                       
1243 
                        last_hit = inters_i;
 1278 
                        last_hit = inters_i;
1244 
                }      
 1279 
                }      
1245 
        }
 1280 
        }
1246 
       
 1281