Subversion Repositories f9daq

Rev

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

  1. #include "include/guide.h"
  2.  
  3. #include <iostream>
  4.  
  5. // vector output shortcut
  6. void printv(TVector3 v)
  7. {
  8.   printf("(x,y,z) = (%.4lf, %.4lf, %.4lf)\n", v.x(), v.y(), v.z());
  9. }
  10. // TVector3::Rotate does not seem accurate enough
  11. TVector3 rotatey(TVector3 v, double theta)
  12. {
  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));
  16. }
  17. // another shortcut not found in TMath
  18. int sign(double in)
  19. {
  20.   if(in >= 0.0) return 1;
  21.   else return -1;
  22. }
  23.  
  24. void CRay::Set(TVector3 r0, TVector3 k0)
  25. {
  26.   r = r0; k = k0.Unit();
  27. }
  28. //-----------------------------------------------------------------------------
  29. //void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0)
  30. //{
  31. //r.SetXYZ(x0, y0, z0);
  32. //n.SetXYZ(l0, m0, n0); n = n.Unit();
  33. //}
  34. //-----------------------------------------------------------------------------
  35.  
  36. CRay& CRay::operator = (const CRay& rhs)
  37. {
  38.         r.SetXYZ(rhs.GetR().x(), rhs.GetR().y(), rhs.GetR().z());
  39.         //this->r.SetXYZ(p.x(), p.y(), p.z());
  40.         k.SetXYZ(rhs.GetK().x(), rhs.GetK().y(), rhs.GetK().z());
  41.         p.SetXYZ(rhs.GetP().x(), rhs.GetP().y(), rhs.GetP().z());
  42.         return *this;
  43. }
  44. void CRay::Print()
  45. {
  46.   printf("---> CRay::Print() <---\n");
  47.   printf("(x,y,z)=(%.2lf, %.2lf, %.2lf); (l,m,n)=(%.2lf, %.2lf, %.2lf)\n",
  48.       r.x(), r.y(), r.z(), k.x(), k.y(), k.z());
  49. }
  50. void CRay::Draw()
  51. {
  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());
  56.   line3d->SetPoint(1, r.x() + t*k.x(), k.y() + t*k.y(), r.z() + t*k.z());
  57.   line3d->SetLineWidth(1);
  58.   line3d->SetLineColor(color);
  59.  
  60.   line3d->Draw();
  61. }
  62. void CRay::Draw(double x_from, double x_to)
  63. {
  64.   double A1, A2;
  65.   TPolyLine3D *line3d = new TPolyLine3D(2);
  66.  
  67.   if(k.x() < MARGIN) {
  68.       A1 = A2 = 0.0;
  69.   } else {
  70.       A1 = (x_from - r.x())/k.x();
  71.       A2 = (x_to - r.x())/k.x();
  72.   }
  73.  
  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());
  76.   line3d->SetLineWidth(1);
  77.   line3d->SetLineColor(color);
  78.  
  79.   line3d->Draw();
  80. }
  81. void CRay::DrawS(double x_from, double t)
  82. {
  83.   double A1;
  84.   TPolyLine3D *line3d = new TPolyLine3D(2);
  85.  
  86.   if(k.x() < MARGIN)
  87.     A1 = 0.0;
  88.   else
  89.     A1 = (x_from - r.x())/k.x();
  90.  
  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());
  93.   line3d->SetLineWidth(1);
  94.   line3d->SetLineColor(color);
  95.  
  96.   line3d->Draw();
  97. }
  98.  
  99.  
  100. CPlane4::CPlane4() :
  101.             n(TVector3(1.0, 0.0, 0.0)),
  102.             A(0),
  103.             B(0),
  104.             C(0),
  105.             D(0)
  106. { r[0] = TVector3(0.0,-1.0,-1.0);
  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;
  112. };
  113. CPlane4::CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
  114. {
  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;
  123.  
  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();
  127.  
  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);
  132.  
  133.   r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
  134.   n.SetXYZ(A, B, C);
  135.   n = n.Unit();
  136.  
  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()) ));
  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.  
  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();
  151.  
  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);
  156.  
  157.   r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
  158.   n.SetXYZ(A, B, C);
  159.   n = n.Unit();
  160.  
  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()) ));
  166. };
  167.  
  168. CPlane4::CPlane4(TVector3 *vr)
  169. {
  170.   double x1,y1,z1, x2,y2,z2, x3,y3,z3;
  171.  
  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();
  175.  
  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);
  180.  
  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();
  184.  
  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()) ));
  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. {
  196.   TVector3 N; //nenormirani vektor (A,B,C)
  197.   double num, den; //stevec, imenovalec
  198.   double t;
  199.   TVector3 tmp;
  200.  
  201.   N.SetXYZ(A,B,C);
  202.  
  203.   num = N*ray.GetR() + D;
  204.   den = N*ray.GetK();
  205.  
  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();
  224.   tmp -= t*ray.GetK();
  225.   *vec = tmp;
  226.   return 1;
  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.  
  235.   if(dbg) printf("--- CPlane4::IsInTri ---\n");
  236.  
  237.   angle_ve1 = TMath::ACos(/*TMath::Abs*/( (e1.Unit()) * (vec.Unit()) ));
  238.   angle_ve2 = TMath::ACos(/*TMath::Abs*/( (e2.Unit()) * (vec.Unit()) ));
  239.  
  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.  
  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.  
  259.   if(dbg) printf("--- CPlane4::IsVectorIn ---\n");
  260.  
  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;
  269. }
  270. int CPlane4::TestIntersection(CRay in)
  271. {
  272.   TVector3 tmp;
  273.  
  274.   if( GetIntersection(&tmp, in) )
  275.     if( IsVectorIn(tmp) )
  276.       return 1;
  277.  
  278.   return 0;
  279. }
  280. int CPlane4::TestIntersection(TVector3 *vec, CRay in)
  281. {
  282.   TVector3 tmp;
  283.  
  284.   if( GetIntersection(&tmp, in) )
  285.     if( IsVectorIn(tmp) ) {
  286.         *vec = tmp;
  287.         return 1;
  288.     }
  289.  
  290.   return 0;
  291. }
  292. void CPlane4::Print()
  293. {
  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);
  300. }
  301. void CPlane4::Draw(int color, int width)
  302. {
  303.   TPolyLine3D *line3d = new TPolyLine3D(5);
  304.  
  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);
  308.  
  309.   line3d->Draw();
  310. }
  311.  
  312.  
  313. CSurface::CSurface(int type0):
  314.       type(type0)
  315. {
  316.   TVector3 vr[4];
  317.   TDatime now;
  318.  
  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();
  330. }
  331. CSurface::CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
  332. {
  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();
  342. }
  343. CSurface::CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
  344. {
  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();
  354. }
  355. void CSurface::SetIndex(double n10, double n20)
  356. {
  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;
  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
  368. int CSurface::PropagateRay(CRay in, CRay& out, TVector3 *intersection)
  369. {
  370.   if (dbg) printf("--- CSurface::PropagateRay ---\n");
  371.   double cosTi; // incident ray angle
  372.   double cosTt; // transmited ray angle
  373.   TVector3 intersect, transmit;
  374.  
  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
  396.   TVector3 pol_t = in.GetP(); // incident polarization
  397.   int sign_n; // sign of normal direction vs. inbound ray
  398.   double cosTN; // debug
  399.  
  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
  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
  408.       // k in this notation is in.GetK()
  409.       v_te = n.Cross(in.GetK());
  410.       v_te = v_te.Unit();
  411.       v_tm = -v_te.Cross(in.GetK());
  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();
  419.       double alpha = TMath::ACos(cosAf);
  420.       if(alpha < MARGIN) alpha = MARGIN;
  421.       if(dbg) printf("  cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE);
  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);
  426.   //}
  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:
  438.     cosTi = in.GetK() * n;
  439.     if(dbg) printf("  cosTi = %lf (Ti = %lf)\n", cosTi, TMath::ACos(cosTi)*DEGREE);
  440.     sign_n = -sign(cosTi);
  441.     if(dbg) printf("  sign_n = %d\n", sign_n);
  442.     cosTi = TMath::Abs(cosTi);
  443.  
  444.     // Check if there can be total reflection: n1 > n2
  445.     if(N1_N2(-sign_n) < 1.0)
  446.       cosTtotal = TMath::Sqrt( 1 - TMath::Power(N1_N2(-sign_n), 2) );
  447.     else
  448.       cosTtotal = 0.0;
  449.  
  450.     if(dbg) printf("  cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
  451.     // 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");
  458.         transmit = in.GetK() + sign_n*2*cosTi*n;
  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.         }
  464.         out.Set(intersect, transmit);
  465.  
  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;
  475.         //alpha += delta;
  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);
  482.         out.setPolarization(pol_t);
  483.  
  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.  
  492.         transmit = N1_N2(sign_n)*in.GetK() + sign_n*(N1_N2(sign_n)*cosTi - cosTt)*n;
  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();
  509.         pol_t = a_te*(1.0 - TMath::Abs(r_te))*v_te + a_tm*(1.0 - TMath::Abs(r_tm))*v_tm_t;
  510.  
  511.         if(dbg) {
  512.             printf("  v_tm_t = "); printv(v_tm_t);
  513.             printf("  pol_t = "); printv(pol_t);
  514.         }
  515.  
  516.         // Fresnel coefficients
  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;
  520.  
  521.         if (dbg) printf("  R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
  522.     }
  523.  
  524.     if(p_ref >= R_f) { // se lomi
  525.         if (dbg) printf("   SURFACE REFRACTED. Return.\n");
  526.         out.Set(intersect, transmit);
  527.         out.setPolarization(pol_t);
  528.         return REFRACTION;
  529.     } else { // se odbije
  530.         if (dbg) printf("   SURFACE REFLECTED. p_ref=%f, R_f=%f\n", p_ref, R_f);
  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();
  534.         out.Set(intersect, transmit);
  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;
  537.         out.setPolarization(pol_t);
  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
  550.         cosTi = in.GetK() * n;
  551.         transmit = in.GetK() - 2*cosTi*n;
  552.         out.Set(intersect, transmit);
  553.         return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
  554.     } else { // se ne odbije
  555.         transmit = in.GetK();
  556.         out.Set(intersect, transmit);
  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
  565.         cosTi = in.GetK() * n;
  566.         if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj
  567.             transmit = in.GetK() - 2*cosTi*n;
  568.             out.Set(intersect, transmit);
  569.         } else { // ni tot. odboja
  570.             transmit = in.GetK();
  571.             out.Set(intersect, transmit);
  572.             return ABSORBED;
  573.         }
  574.     } else { // se ne odbije
  575.         transmit = in.GetK();
  576.         out.Set(intersect, transmit);
  577.         return ABSORBED;
  578.     }
  579.     break;
  580.  
  581.   default:
  582.     out = in;
  583.     break;
  584.   }
  585.  
  586.   return REFRACTION;
  587. }
  588.  
  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())
  598. {
  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();
  609.  
  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);
  621.  
  622.   for(int i = 0; i<8; i++) vodnik_edge[i] += center;
  623.  
  624.   // light guide surfaces
  625.   s_side[0].Set(SURF_REFRA, vodnik_edge, n0, _n2, _r);
  626.   s_side[0].FlipN();
  627.  
  628.   s_side[1].Set(SURF_REFRA, vodnik_edge[3], vodnik_edge[2],
  629.       vodnik_edge[6], vodnik_edge[7], _n2, _n1, _r);
  630.   s_side[2].Set(SURF_REFRA, vodnik_edge[2], vodnik_edge[1],
  631.       vodnik_edge[5], vodnik_edge[6], _n2, _n1, _r);
  632.   s_side[3].Set(SURF_REFRA, vodnik_edge[1], vodnik_edge[0],
  633.       vodnik_edge[4], vodnik_edge[5], _n2, _n1, _r);
  634.   s_side[4].Set(SURF_REFRA, vodnik_edge[0], vodnik_edge[3],
  635.       vodnik_edge[7], vodnik_edge[4], _n2, _n1, _r);
  636.   // n3 - ref ind at the exit, grease, air
  637.   s_side[5].Set(SURF_REFRA, &vodnik_edge[4], _n2, _n3, _r);
  638.   s_side[5].FlipN();
  639.   // exit surface in the case of bad coupling
  640.   noCoupling.Set(SURF_REFRA, &vodnik_edge[4], _n2, 1.0, _r);
  641.   noCoupling.FlipN();
  642.   // grease = specific pattern area of coupling
  643.   TVector3 activePosition(center);
  644.   activePosition += TVector3(_d, 0, 0);
  645.   TVector3 normal(1,0,0);
  646.   grease.Set(activePosition, normal, 0.95*a/2.0);
  647.  
  648.   if(fresnel) for(int i=0; i<6; i++) s_side[i].SetFresnel(1);
  649.  
  650.   // statistics histograms
  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");
  661.  
  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);
  664.  
  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);
  667.  
  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);
  671.  
  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);
  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                                            
  684. Fate Guide::PropagateRay(CRay in, CRay& out, int *n_points, TVector3 *points)
  685. {
  686.   if (dbg) printf("--- GUIDE::PropagateRay ---\n");
  687.   // ray0 - incident ray
  688.   // ray1 - trans/refl ray
  689.   CRay ray0;
  690.   CRay ray1;
  691.   TVector3 vec0, vec1;
  692.   int inters_i = 0;
  693.  
  694.   ray0 = in;
  695.   int n_odb = 0;
  696.   int last_hit = 0;
  697.   int propagation = 0;
  698.   int result = s_side[0].PropagateRay(ray0, ray1, &vec1);
  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();
  711.       hfate.Fill(enter); // enter
  712.       hin.Fill(vec1.y(), vec1.z());
  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) {
  723.                   int testBoundary = s_side[inters_i].TestIntersection(&vec1, ray1);
  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
  738.           propagation = s_side[inters_i].PropagateRay(ray0, ray1, &vec1);
  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);
  751.                   bool hitActive = grease.TestIntersection(&hitVector, ray0);
  752.                   if (hitActive and dbg) printf("   GUIDE: hit grease\n");
  753.                   if (!hitActive) propagation = noCoupling.PropagateRay(ray0, ray1, &vec1);
  754.               }
  755.               // check on which side the vector is?
  756.               TVector3 ray = ray1.GetK();
  757.               TVector3 exitNormal = s_side[5].GetN();
  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;
  768.               hout.Fill(vec1.y(), vec1.z());
  769.               hnodb_exit.Fill(n_odb-1);
  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.  
  805.   hfate.Fill(fate);
  806.   hfate.Fill(rays);
  807.   hnodb_all.Fill(n_odb-2);
  808.   *n_points = n_odb+1;
  809.   out = ray0;
  810.   return fate;
  811. }
  812. void Guide::GetVFate(int *out)
  813. {
  814.   for(int i=0;i<7;i++) out[i] = (int)hfate.GetBinContent(i+1);
  815. }
  816. void Guide::Draw(int color, int width)
  817. {
  818.   for(int i = 0; i<6; i++) (&s_side[i])->Draw(color, width);
  819. }
  820. void Guide::DrawSkel(int color, int width)
  821. {
  822.   TPolyLine3D *line3d = new TPolyLine3D(2);
  823.   line3d->SetLineWidth(width);
  824.   line3d->SetLineColor(color);
  825.  
  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.   }
  831. }
  832.  
  833.  
  834. int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
  835. {
  836.   double num, den; //stevec, imenovalec
  837.   double t;
  838.   TVector3 tmp;
  839.  
  840.   if(dbg) printf("---> CPlaneR::TestIntersection <---\n");
  841.   if(dbg) {printf("c = "); printv(center); printf(" | n = "); printv(n); printf("\n");}
  842.  
  843.   double D = - n*center;
  844.   num = n*ray.GetR() + D;
  845.   den = n*ray.GetK();
  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();
  861.   tmp -= t*ray.GetK();
  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;
  871. }
  872.  
  873. void CPlaneR::Draw(int color, int width)
  874. {
  875.   const int NN = 32;
  876.   double phi, x, y;
  877.  
  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();
  890. }
  891.  
  892.  
  893. CDetector::CDetector(TVector3 center0, DetectorParameters& parameters) :
  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()),
  905.       guide(Guide(center0, parameters)),
  906.       plate(Plate(parameters)),
  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.   //{
  917.   double d = parameters.getD();
  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];
  935.   int nBins = nch;
  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);
  941.   glass.Set(SURF_REFRA, plane_v, n3, n4, reflectivity);
  942.   glass.FlipN();
  943.  
  944.   // additional circular glass between LG and SiPM
  945.   glass_circle.Set(TVector3(x_offset+d+glass_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), b);
  946.  
  947.   //hglass = (TH2F*)gROOT->FindObject("hglass"); if(hglass) delete hglass;
  948.   hglass = TH2F("hglass", "",
  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);
  960.   active.Set(plane_v);
  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);
  965.   grease.Set(activePosition, normal, 1.0*p_size);
  966.  
  967.   //hactive = (TH2F*)gROOT->FindObject("hactive"); if(hactive) delete hactive;
  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);
  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);
  970.  
  971.   p_size = b/2.0;
  972.   //p_size = 2.5;
  973.   //p_size = M*0.6;
  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);
  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);
  983.   detector.Set(plane_v);
  984.  
  985.   //hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector;
  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);
  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);
  988.  
  989.   /*
  990.         window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R);  
  991.  
  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();
  998.  
  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);
  1001.    */
  1002.   p_size = b/2.0;
  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);
  1007. }
  1008.  
  1009. //-----------------------------------------------------------------------------
  1010. // vrne 1 ce je zadel aktvino povrsino
  1011. // vrne <1 ce jo zgresi
  1012. int CDetector::Propagate(CRay in, CRay& out, int draw)
  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");
  1023.   generated.Fill(in.GetR().y(), in.GetR().z());
  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);
  1027.   CRay rayout(in);
  1028.   rayout.SetColor(col_in);
  1029.  
  1030.   const int max_n_points = guide.GetMAXODB() + 2;
  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);
  1038.  
  1039.   // Draw the plate and propagate the ray through
  1040.   // check if the ray should be reflected??
  1041.  
  1042.   if(_plateOn) {
  1043.  
  1044.       fatePlate = plate.propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
  1045.       if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
  1046.       if(draw) {
  1047.           if(fatePlate == missed) {
  1048.               rayout.SetColor(col_in);
  1049.               rayout.DrawS(center.x() - _plateWidth, -10.0);
  1050.           }
  1051.           else if(fatePlate == backreflected){
  1052.               rayout.SetColor(kBlack);
  1053.               rayout.DrawS(center.x() - _plateWidth, 7.0);
  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.               }
  1063.               rayout.DrawS(pointsPlate[p_i].x(), -0.1);
  1064.               if(fatePlate == noreflection) { // lost on plate side
  1065.                   rayout.SetColor(col_out);
  1066.                   rayout.DrawS(pointsPlate[p_i].x(), 10.0);
  1067.               }
  1068.           }
  1069.       }
  1070.  
  1071.       if(! (fatePlate == hitExit or fatePlate == refracted) ) {
  1072.           guide.GetHFate()->Fill(rays);
  1073.           if (dbg)printf("CDetector::propagate Simulated ray missed the entry surface!\n");
  1074.           if (fatePlate == backreflected)
  1075.             guide.GetHFate()->Fill(fatePlate); // reflected back
  1076.           else
  1077.             guide.GetHFate()->Fill(noreflection); //lost on plate side
  1078.           return fatePlate;
  1079.       }
  1080.  
  1081.       //Ray hits light guide
  1082.       histoPlate.Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
  1083.  
  1084.   }
  1085.   else {
  1086.       //rayout = rayin;
  1087.       if(draw) rayout.DrawS(center.x(), -10.0);
  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;
  1099.   CRay *ray0 = new CRay(rayout);
  1100.   // delete rayout; -> creates dangling reference when tries to delete ray0!
  1101.   //delete rayin; -> delete rayout!
  1102.   CRay ray1;
  1103.  
  1104.   fate = guide.PropagateRay(*ray0, ray1, &n_points, points);
  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");
  1115.               TVector3 r = ray1.GetR();
  1116.               TVector3 k = ray1.GetK();
  1117.               ray1.Set(r,k);
  1118.               ray1.DrawS(center.x(), 10.0);
  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);
  1127.                   if(glass_on) {/*if(fate == 1)*/ ray1.Draw(points[p_i].x(), center.x() + guide.getD() + glass_d);}
  1128.                   else {
  1129.                       ray1.SetColor(col_out);
  1130.                       ray1.DrawS(points[p_i].x(), 10.0);
  1131.                   }
  1132.               }
  1133.           }
  1134.       }
  1135.  
  1136.  
  1137.       if(! (fate == hitExit or fate == refracted) ) {
  1138.           if (dbg) printf("Detector: fate != hit, refracted\n");
  1139.           out = ray1;
  1140.           delete ray0;
  1141.           //delete ray1;
  1142.           //delete rayout;
  1143.           delete rayin;
  1144.           return fate;
  1145.       }
  1146.   } else {
  1147.       if (dbg) printf("Detector: fate = hit or refracted");
  1148.       ray1 = *ray0;
  1149.       if(draw) {
  1150.           //double epoxy = parameters->getGlassD();
  1151.           if(glass_on) ray1.Draw(center.x(), center.x() + glass_d);
  1152.           else ray1.DrawS(center.x(), 10.0);
  1153.       }
  1154.   }
  1155.  
  1156.   fate = missed; // zgresil aktivno povrsino
  1157.   if(glass_on) {
  1158.       *ray0 = ray1;
  1159.       ray1.SetColor(col_rgla);
  1160.       fate_glass = glass.PropagateRay(*ray0, ray1, &presecisce);
  1161.       if(fate_glass == REFRACTION) {
  1162.           hglass.Fill(presecisce.y(), presecisce.z());
  1163.           if(draw) ray1.DrawS(presecisce.x(), 10.0);
  1164.           //if(active->TestIntersection(&presecisce, *ray1)) {
  1165.           //fate = hitExit;
  1166.           //h-Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
  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
  1173.             if(draw) ray1.DrawS(presecisce.x(), 10.0);
  1174.       }
  1175.   }
  1176.  
  1177.   // Main test: ray and SiPM surface
  1178.   if(active.TestIntersection(&presecisce, ray1)) {
  1179.       fate = hitExit;
  1180.       hactive.Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
  1181.       hlaser.Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ);
  1182.   }
  1183.   // If it is on the same plane as SiPM
  1184.   if(detector.TestIntersection(&presecisce, ray1))
  1185.     hdetector.Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
  1186.   //}
  1187.   //} else {
  1188.   //if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
  1189.   //}
  1190.  
  1191.   out = ray1;
  1192.   delete ray0;
  1193.   //delete ray1;
  1194.   //delete rayout;
  1195.   delete rayin;
  1196.   return fate;
  1197. }
  1198.  
  1199. void CDetector::Draw(int width)
  1200. {
  1201.   if(guide_on) {
  1202.       if( TMath::Abs(guide.getN1()-guide.getN2()) < MARGIN ) {
  1203.           if(_plateOn) plate.drawSkel(col_LG, width);
  1204.           guide.DrawSkel(col_LG, width);
  1205.       }
  1206.       else {
  1207.           if(_plateOn) plate.draw(4, width);
  1208.           guide.Draw(col_LG, width);
  1209.       }
  1210.   }
  1211.  
  1212.   if(glass_on) glass_circle.Draw(col_glass, width);
  1213.   //window_circle->Draw(col_glass, width);
  1214.   active.Draw(col_active, width);
  1215. }
  1216.  
  1217.  
  1218. Plate::Plate(DetectorParameters& parameters)
  1219. {
  1220.   TVector3 center = CENTER;
  1221.   const double b = parameters.getB();
  1222.   const double n1 = parameters.getN1();
  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.  
  1239.   sides[0].Set(SURF_REFRA, plate_edge, n1, n2, c_reflectivity);
  1240.   sides[0].FlipN();
  1241.  
  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);
  1246.  
  1247.   sides[5].Set(SURF_REFRA, &plate_edge[4], n2, n2, c_reflectivity);
  1248.   sides[5].FlipN();
  1249.  
  1250.   for(int i=0; i<6; i++) sides[i].SetFresnel(1);
  1251. }
  1252.  
  1253. void Plate::draw(int color, int width)
  1254. {
  1255.   for(int i = 0; i<6; i++) sides[i].Draw(color, width);
  1256. }
  1257.  
  1258. void Plate::drawSkel(int color, int width)
  1259. {
  1260.   TPolyLine3D line3d(2);
  1261.   line3d.SetLineWidth(width);
  1262.   line3d.SetLineColor(color);
  1263.  
  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.   }
  1269. }
  1270.  
  1271. Fate Plate::propagateRay(CRay in, CRay& out, int *n_points, TVector3 *points)
  1272. {
  1273.   CRay ray0;
  1274.   CRay ray1;
  1275.   TVector3 vec0, vec1;
  1276.   Fate fate = enter;
  1277.   int inters_i = 0;
  1278.  
  1279.   ray0 = in;
  1280.   int n_odb = 0;
  1281.   int last_hit = 0;
  1282.   int propagation = 0;
  1283.  
  1284.   int result = sides[0].PropagateRay(ray0, ray1, &vec1);
  1285.   if( !result ) {
  1286.       // ce -NI- presecisca z vstopno
  1287.       fate = missed;
  1288.   } else if(result == REFLECTION) {
  1289.       if (dbg) printf("PLATE: reflected\n");
  1290.       ray0 = ray1;
  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) {
  1302.                   if( sides[inters_i].TestIntersection(&vec1, ray1) ) break;
  1303.               }
  1304.           }
  1305.           points[n_odb] = vec1;
  1306.           if(inters_i == 0) {
  1307.               ray0 = ray1;
  1308.               fate = backreflected;
  1309.               break;} // backreflection
  1310.  
  1311.           propagation = sides[inters_i].PropagateRay(ray0, ray1, &vec1);
  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;
  1338.   out = ray0;
  1339.   return fate;
  1340. };
  1341.  
  1342.  
  1343.