Subversion Repositories f9daq

Rev

Rev 73 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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