Subversion Repositories f9daq

Rev

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

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