Subversion Repositories f9daq

Rev

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