Subversion Repositories f9daq

Rev

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