WebSVN – f9daq – Diff – /lightguide/trunk/src/guide.cpp

 int sign(double in)
+{
   if(in >= 0.0) return 1;
   else return -1;
+}
-//=================================================================================
-//-----------------------------------------------------------------------------
-void CRay::Set(TVector3 r0, TVector3 n0)
+void CRay::Set(TVector3 r0, TVector3 k0)
+{
-  r = r0; n = n0.Unit();
+  r = r0; k = k0.Unit();
+}
 //-----------------------------------------------------------------------------
 //void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0)
 //{
 //r.SetXYZ(x0, y0, z0);
         r.SetXYZ(p.GetR().x(), p.GetR().y(), p.GetR().z());
         //this->r.SetXYZ(p.x(), p.y(), p.z());
         n.SetXYZ(p.GetN().x(), p.GetN().y(), p.GetN().z());
         return *this;
 } */
-//-----------------------------------------------------------------------------
 void CRay::Print()
 {
   printf("---> CRay::Print() <---\n");
   printf("(x,y,z)=(%.2lf, %.2lf, %.2lf); (l,m,n)=(%.2lf, %.2lf, %.2lf)\n",
-      r.x(), r.y(), r.z(), n.x(), n.y(), n.z());
+      r.x(), r.y(), r.z(), k.x(), k.y(), k.z());
+}
-//-----------------------------------------------------------------------------
 void CRay::Draw()
+{
   double t = 50.0;
   TPolyLine3D *line3d = new TPolyLine3D(2);
   //line3d->SetPoint(0, r.x() - t*n.x(), r.y() - t*n.y(), r.z() - t*n.z());
   line3d->SetPoint(0, r.x(), r.y(), r.z());
-  line3d->SetPoint(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.z());
+  line3d->SetPoint(1, r.x() + t*k.x(), k.y() + t*k.y(), r.z() + t*k.z());
   line3d->SetLineWidth(1);
   line3d->SetLineColor(color);
   line3d->Draw();
+}
-//-----------------------------------------------------------------------------
 void CRay::Draw(double x_from, double x_to)
+{
   double A1, A2;
   TPolyLine3D *line3d = new TPolyLine3D(2);
-  if(n.x() < MARGIN) {
+  if(k.x() < MARGIN) {
       A1 = A2 = 0.0;
   } else {
-      A1 = (x_from - r.x())/n.x();
+      A1 = (x_from - r.x())/k.x();
-      A2 = (x_to - r.x())/n.x();
+      A2 = (x_to - r.x())/k.x();
+  }
-  line3d->SetPoint(0, x_from, A1*n.y()+r.y(), A1*n.z()+r.z());
+  line3d->SetPoint(0, x_from, A1*k.y()+r.y(), A1*k.z()+r.z());
-  line3d->SetPoint(1, x_to, A2*n.y()+r.y(), A2*n.z()+r.z());
+  line3d->SetPoint(1, x_to, A2*k.y()+r.y(), A2*k.z()+r.z());
   line3d->SetLineWidth(1);
   line3d->SetLineColor(color);
   line3d->Draw();
 }
-//-----------------------------------------------------------------------------
 void CRay::DrawS(double x_from, double t)
+{
   double A1;
   TPolyLine3D *line3d = new TPolyLine3D(2);
-  if(n.x() < MARGIN)
+  if(k.x() < MARGIN)
     A1 = 0.0;
   else
-    A1 = (x_from - r.x())/n.x();
+    A1 = (x_from - r.x())/k.x();
-  line3d->SetPoint(0, x_from, A1*n.y()+r.y(), A1*n.z()+r.z());
+  line3d->SetPoint(0, x_from, A1*k.y()+r.y(), A1*k.z()+r.z());
-  line3d->SetPoint(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.z());
+  line3d->SetPoint(1, r.x() + t*k.x(), r.y() + t*k.y(), r.z() + t*k.z());
   line3d->SetLineWidth(1);
   line3d->SetLineColor(color);
   line3d->Draw();
+}
-//=================================================================================
-//=================================================================================
 CPlane4::CPlane4() :
             n(TVector3(1.0, 0.0, 0.0)),
             A(0),
             B(0),
             C(0),
 r[2] = TVector3(0.0, 1.0, 1.0);
 r[3] = TVector3(0.0, 1.0,-1.0);
 for(int i=0;i<4;i++) edge[i] = TVector3(0,0,0);
 for(int i=0;i<4;i++) angle_r[i] = 0;
 };
-//-----------------------------------------------------------------------------
 CPlane4::CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
+{
   //Set(r1, r2, r3, r4);
   //}
   //-----------------------------------------------------------------------------
   TVector3 tmp;
   N.SetXYZ(A,B,C);
   num = N*ray.GetR() + D;
-  den = N*ray.GetN();
+  den = N*ray.GetK();
   if (dbg) printf("t = %6.3lf / %6.3lf =  %6.3lf\n", num, den, num/den);
   //if(den == 0)
   if(TMath::Abs(den) < MARGIN) {
+  }
   t = num / den;
   tmp = ray.GetR();
-  tmp -= t*ray.GetN();
+  tmp -= t*ray.GetK();
   *vec = tmp;
   return 1;
+}
 //-----------------------------------------------------------------------------
 // ali je vektor vec, ki lezi na ravnini skupaj z e1 in e2, med njima
       if(!status) return 0;
+    }
   return 1;
+}
-//-----------------------------------------------------------------------------
 int CPlane4::TestIntersection(CRay in)
+{
   TVector3 tmp;
   if( GetIntersection(&tmp, in) )
     if( IsVectorIn(tmp) )
       return 1;
   return 0;
+}
-//-----------------------------------------------------------------------------
 int CPlane4::TestIntersection(TVector3 *vec, CRay in)
+{
   TVector3 tmp;
   if( GetIntersection(&tmp, in) )
         return 1;
+    }
   return 0;
+}
-//-----------------------------------------------------------------------------
 void CPlane4::Print()
 {
   printf("--- CPlane4::Print() ---\n");
   printf("  r=(%.2lf, %.2lf, %.2lf); n=(%.2lf, %.2lf, %.2lf); ",
       r[0].x(), r[0].y(), r[0].z(), n.x(), n.y(), n.z());
   printf(  "(A,B,C,D)=(%.2lf, %.2lf, %.2lf, %.2lf) \n", A, B, C, D);
   for(int i=0;i<4;i++) printf("  edge[%d] = (%lf, %lf, %lf)\n", i, edge[i].x(), edge[i].y(), edge[i].z());
   for(int i=0;i<4;i++) printf("  angle[%d] = %lf\n", i, angle_r[i]*DEGREE);
+}
-//-----------------------------------------------------------------------------
 void CPlane4::Draw(int color, int width)
+{
   TPolyLine3D *line3d = new TPolyLine3D(5);
   for(int i=0;i<4;i++) line3d->SetPoint(i, r[i].x(), r[i].y(), r[i].z());
   line3d->SetPoint(4, r[0].x(), r[0].y(), r[0].z());
   line3d->SetLineWidth(width); line3d->SetLineColor(color);
   line3d->Draw();
 }
-//=================================================================================
-//=================================================================================
 CSurface::CSurface(int type0):
       type(type0)
+{
   TVector3 vr[4];
   TDatime now;
   vr[3].SetXYZ(0.0, 1.0,-1.0);
   //CPlane4::Set(vr);
   SetIndex(1.0, 1.5);
   reflection = c_reflectivity;
   rand.SetSeed(now.Get());
   SetFresnel();
 }
-//-----------------------------------------------------------------------------
 CSurface::CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
+{
   TDatime now;
   type = type0; CPlane4::Set(r1, r2, r3, r4);
   SetIndex(n10, n20);
   reflection = reflectivity;
   rand.SetSeed(now.Get());
   SetFresnel();
 }
-//-----------------------------------------------------------------------------
 CSurface::CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
+{
   TDatime now;
   type = type0; CPlane4::Set(vr);
   SetIndex(n10, n20);
   reflection = reflectivity;
   rand.SetSeed(now.Get());
   SetFresnel();
+}
-//-----------------------------------------------------------------------------
 void CSurface::SetIndex(double n10, double n20)
 {
   n1 = n10; n2 = n20; n1_n2 = n1/n2;
   if(n1 > n2)
     cosTtotal = TMath::Sqrt( 1 - TMath::Power(n2/n1, 2) );
   else
   double a_te = 0.0; // s-wave amplitude, cos Alpha
   double a_tm = 0.0; // p-wave amplitude, sin Alpha
   TVector3 v_te; // unit s-polarization vector
   TVector3 v_tm; // unit p-polarization vector
   TVector3 v_tm_t;// transmited polarization parallel with the plane of incidence
-  TVector3 pol_t = in.GetP(); // transmited polarization
+  TVector3 pol_t = in.GetP(); // incident polarization
   int sign_n; // sign of normal direction vs. inbound ray
   double cosTN; // debug
+  // Decomposition of incident polarization vector
+  // using unit vectors v_tm & v_te
+  // in a_tm and a_te components
-  if(fresnel) {
+  //if(fresnel) {
-      // p-polarization unit vector v_te
+      // s-polarization unit vector v_te
       // is in the plane orthogonal to the plane of incidence
       // defined as the plane spanned by
       // incident surface vector n and wave vector k
-      // k in this notation is in.GetN()
+      // k in this notation is in.GetK()
-      v_te = n.Cross(in.GetN());
+      v_te = n.Cross(in.GetK());
       v_te = v_te.Unit();
-      v_tm = -v_te.Cross(in.GetN());
+      v_tm = -v_te.Cross(in.GetK());
       v_tm = v_tm.Unit();
       if(dbg) {
           printf("  v_te = "); printv(v_te);
           printf("  v_tm = "); printv(v_tm);
+      }
       double cosAf = v_te * in.GetP();
+      double alpha = acos(cosAf);
-      if(dbg) printf("  cosAf = %lf (Af = %lf)\n", cosAf, TMath::ACos(cosAf)*DEGREE);
+      if(dbg) printf("  cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE);
       a_te = cosAf;
       a_tm = TMath::Sqrt(1 - cosAf*cosAf);
       if(dbg) printf("  a_te = %lf, a_tm = %lf\n", a_te, a_tm);
-  }
+  //}
   // ----------------------------------------------------------------------------
   // reflection probability
   double p_ref = rand.Uniform(0.0, 1.0);
   switch(type){
   // ----------------------------------------------------------------------------
   // --------------- refraction from n1 to n2 -----------------------------------
   // ----------------------------------------------------------------------------
   case SURF_REFRA:
-    cosTi = in.GetN() * n;
+    cosTi = in.GetK() * n;
     if(dbg) printf("  cosTi = %lf (Ti = %lf)\n", cosTi, TMath::ACos(cosTi)*DEGREE);
     sign_n = -sign(cosTi);
     if(dbg) printf("  sign_n = %d\n", sign_n);
     cosTi = TMath::Abs(cosTi);
     if (dbg) printf("   reflection probability = %f\n", p_ref);
     // If n1>n2 and theta>thetaCritical, total reflection
     if(cosTi < cosTtotal) {
         if(dbg) printf("  TOTAL\n");
-        transmit = in.GetN() + sign_n*2*cosTi*n;
+        transmit = in.GetK() + sign_n*2*cosTi*n;
         if(dbg) {
             cosTN = TMath::Abs(transmit.Unit() * n);
             printf("  cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE);
+        }
         out->Set(intersect, transmit);
+        // Shift implemented, but only linear polarization is implemented
+        if (dbg) printf("CSurface: Propagate TOTAL\n");
+        v_tm_t = -v_te.Cross(transmit);
+        v_tm_t = v_tm_t.Unit();
+        // shift the p and s components
+        double n12 = N1_N2(-sign_n);
+        double deltaP = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/(pow(n12,2)*cosTi));
+        double deltaS = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/cosTi);
+        double delta = deltaP - deltaS;
-        // Shift?
+        alpha += delta;
+        a_tm = sin(alpha);
+        a_te = cos(alpha);
+        if (dbg) printf("  deltaP = %f deltaS = %f; new a_tm = %f, a_te = %f",
+            deltaP, deltaS, a_tm, a_te);
-        pol_t = -in.GetP() + sign_n*2*cosTi*n;
+        pol_t = a_tm*v_tm_t + a_te*v_te;
+        if (dbg) printv(pol_t);
-        out->SetPolarization(pol_t);
+        out->setPolarization(pol_t);
         return REFLECTION;
     } else {
         // reflection or refraction according to Fresnel equations
         if(dbg) printf("  REFRACTION\n");
         if(dbg) printf("  N1_N2(sign_n) = %lf\n", N1_N2(sign_n));
         cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));
         if(dbg) printf("  cosTt = %lf (Tt = %lf) \n", cosTt, TMath::ACos(cosTt)*DEGREE);
-        transmit = N1_N2(sign_n)*in.GetN() + sign_n*(N1_N2(sign_n)*cosTi - cosTt)*n;
+        transmit = N1_N2(sign_n)*in.GetK() + sign_n*(N1_N2(sign_n)*cosTi - cosTt)*n;
         if(dbg) {printf("  transmit.Unit() = "); printv(transmit.Unit());}
         if(dbg) {
             cosTN = transmit.Unit() * n;
             printf("  cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE);
+        }
         if(dbg) printf("  r_te = %lf, r_tm = %lf\n", r_te, r_tm);
         // transmited polarization
         v_tm_t = -v_te.Cross(transmit);
         v_tm_t = v_tm_t.Unit();
-        pol_t = a_te * (1.0 -  TMath::Abs(r_te)) * v_te  +  a_tm * (1.0 -  TMath::Abs(r_tm)) * v_tm_t;
+        pol_t = a_te*(1.0 - TMath::Abs(r_te))*v_te + a_tm*(1.0 - TMath::Abs(r_tm))*v_tm_t;
         if(dbg) {
             printf("  v_tm_t = "); printv(v_tm_t);
             printf("  pol_t = "); printv(pol_t);
+        }
+    }
     if(p_ref >= R_f) { // se lomi
         if (dbg) printf("   SURFACE REFRACTED. Return.\n");
         out->Set(intersect, transmit);
-        out->SetPolarization(pol_t);
+        out->setPolarization(pol_t);
         return REFRACTION;
     } else { // se odbije
         if (dbg) printf("   SURFACE REFLECTED. p_ref=%f, R_f=%f\n", p_ref, R_f);
-        transmit = in.GetN() + sign_n*2*cosTi*n;
+        transmit = in.GetK() + sign_n*2*cosTi*n;
+        TVector3 v_tm_r = -v_te.Cross(transmit);
+        v_tm_r = v_tm_r.Unit();
         out->Set(intersect, transmit);
-        pol_t = -in.GetP() + sign_n*2*cosTi*n;
+        //pol_t = -in.GetP() + sign_n*2*cosTi*n;
+        pol_t = a_te*(1.0 - TMath::Abs(r_te))*v_te + a_tm*(1.0 - TMath::Abs(r_tm))*v_tm_r;
-        out->SetPolarization(pol_t);
+        out->setPolarization(pol_t);
         return REFLECTION;
+    }
     //}
     break;
     // --------------- reflection at "reflection" probability ---------------------
     // ----------------------------------------------------------------------------
   case SURF_REFLE:
     p_ref = rand.Uniform(0.0, 1.0);
     if(p_ref < reflection) { // se odbije
-        cosTi = in.GetN() * n;
+        cosTi = in.GetK() * n;
-        transmit = in.GetN() - 2*cosTi*n;
+        transmit = in.GetK() - 2*cosTi*n;
         out->Set(intersect, transmit);
         return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
     } else { // se ne odbije
-        transmit = in.GetN();
+        transmit = in.GetK();
         out->Set(intersect, transmit);
         return ABSORBED;
+    }
     break;
     // total reflection from n1 to n2 with R probbability
   case SURF_IMPER:
     p_ref = rand.Uniform(0.0, 1.0);
     if(p_ref < reflection) { // se odbije
-        cosTi = in.GetN() * n;
+        cosTi = in.GetK() * n;
         if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj
-            transmit = in.GetN() - 2*cosTi*n;
+            transmit = in.GetK() - 2*cosTi*n;
             out->Set(intersect, transmit);
         } else { // ni tot. odboja
-            transmit = in.GetN();
+            transmit = in.GetK();
             out->Set(intersect, transmit);
             return ABSORBED;
+        }
     } else { // se ne odbije
-        transmit = in.GetN();
+        transmit = in.GetK();
         out->Set(intersect, transmit);
         return ABSORBED;
+    }
     break;
     break;
+  }
   return REFRACTION;
+}
-//=================================================================================
-//=================================================================================
 Guide::Guide(TVector3 center0, DetectorParameters &parameters) :
     _d(parameters.getD()),
     _n1(parameters.getN1()),
     _n2(parameters.getN2()),
     _n3(parameters.getN3()),
   noCoupling->FlipN();
   // grease = specific pattern area of coupling
   TVector3 activePosition(center);
   activePosition += TVector3(_d, 0, 0);
   TVector3 normal(1,0,0);
-  grease = new CPlaneR(activePosition, normal, a/2.0);
+  grease = new CPlaneR(activePosition, normal, 0.95*a/2.0);
   if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1);
   // statistics histograms
   hfate = (TH1F*)gROOT->FindObject("hfate"); if(hfate) delete hfate;
 Fate Guide::PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
+{
   if (dbg) printf("--- GUIDE::PropagateRay ---\n");
   // ray0 - incident ray
   // ray1 - trans/refl ray
   CRay ray0;
   CRay ray1;
   TVector3 vec0, vec1;
   int inters_i = 0;
   ray0 = in;
   int n_odb = 0;
                   bool hitActive = grease->TestIntersection(&hitVector, ray0);
                   if (hitActive and dbg) printf("   GUIDE: hit grease\n");
                   if (!hitActive) propagation = noCoupling->PropagateRay(ray0, &ray1, &vec1);
+              }
               // check on which side the vector is?
-              TVector3 ray = ray1.GetN();
+              TVector3 ray = ray1.GetK();
               TVector3 exitNormal = s_side[5]->GetN();
               if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal);
               if (ray*exitNormal > 0) {
                   if (dbg) printf("  GUIDE: ray is backreflected from exit window.\n");
                   fate = backreflected;
   hnodb_all->Fill(n_odb-2);
   *n_points = n_odb+1;
   *out = ray0;
   return fate;
+}
-//-----------------------------------------------------------------------------
 void Guide::GetVFate(int *out)
+{
   for(int i=0;i<7;i++) out[i] = (int)hfate->GetBinContent(i+1);
+}
-//-----------------------------------------------------------------------------
 void Guide::Draw(int color, int width)
 {
   for(int i = 0; i<6; i++) s_side[i]->Draw(color, width);
 }
-//-----------------------------------------------------------------------------
 void Guide::DrawSkel(int color, int width)
+{
   TPolyLine3D *line3d = new TPolyLine3D(2);
   line3d->SetLineWidth(width); line3d->SetLineColor(color);
       line3d->SetPoint(0, vodnik_edge[i+0].x(), vodnik_edge[i+0].y(), vodnik_edge[i+0].z());
       line3d->SetPoint(1, vodnik_edge[i+4].x(), vodnik_edge[i+4].y(), vodnik_edge[i+4].z());
       line3d->DrawClone();
+  }
+}
-//=================================================================================
-//=================================================================================
 int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
+{
   double num, den; //stevec, imenovalec
   double t;
   TVector3 tmp;
   if(dbg) printf("---> CPlaneR::TestIntersection <---\n");
   if(dbg) {printf("c = "); printv(center); printf(" | n = "); printv(n); printf("\n");}
   double D = - n*center;
   num = n*ray.GetR() + D;
-  den = n*ray.GetN();
+  den = n*ray.GetK();
   if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den);
   if(TMath::Abs(den) < MARGIN) {
       if(TMath::Abs(num) < MARGIN)
   t = num / den;
   if(dbg) printf("t = %.4lf | ", t);
   tmp = ray.GetR();
-  tmp -= t*ray.GetN();
+  tmp -= t*ray.GetK();
   *vec = tmp;
   if(dbg) {printv(tmp); printf(" | Rv = %.4lf <> R = %.4lf\n", ((tmp - center).Mag()), _r);}
   if( ((tmp - center).Mag()) < _r )
     return 1;
   else
     return 0;
+}
-//-----------------------------------------------------------------------------
 void CPlaneR::Draw(int color, int width)
+{
   const int NN = 32;
   double phi, x, y;
       y = _r*TMath::Sin(phi);
       arc->SetPoint(i, center.x(),  x,  y);
+  }
   arc->Draw();
+}
-//=================================================================================
-//=================================================================================
 CDetector::CDetector(TVector3 center0, DetectorParameters& parameters) :
       center(center0),
       glass_on(parameters.getGlassOn()),
       glass_d(parameters.getGlassD()),
       col_in(2),
           if(fatePlate == missed) {
               rayout->SetColor(col_in);
               rayout->DrawS(center.x() - _plateWidth, -10.0);
+          }
           else if(fatePlate == backreflected){
+              rayout->SetColor(kBlack);
+              rayout->DrawS(center.x() - _plateWidth, 7.0);
               if (dbg) printf("Backreflected at plate!\n");
+          }
           else {
               int p_i;
               for(p_i = 0; p_i < nPointsPlate-1; p_i++) {
   if(guide_on) {
       if(draw) {
           if(fate == missed) {
               if (dbg) printf("Detector: fate=missed\n");
               TVector3 r = ray1->GetR();
-              TVector3 n = ray1->GetN();
+              TVector3 k = ray1->GetK();
-              ray1->Set(r,n);
+              ray1->Set(r,k);
               ray1->DrawS(center.x(), 10.0);
           } else {
               for(p_i = 0; p_i < n_points-1; p_i++) {
                   line3d->SetPoint(0, points[p_i].x(), points[p_i].y(), points[p_i].z());
                   line3d->SetPoint(1, points[p_i+1].x(), points[p_i+1].y(), points[p_i+1].z());
       ray1 = ray0;
       if(draw) {
           //double epoxy = parameters->getGlassD();
           if(glass_on) ray1->Draw(center.x(), center.x() + glass_d);
           else ray1->DrawS(center.x(), 10.0);
       }
   }
   fate = missed; // zgresil aktivno povrsino
   if(glass_on) {
       *ray0 = *ray1;
       ray1->SetColor(col_rgla);
       fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
           //if(detector->TestIntersection(&presecisce, *ray1))
           //hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
           //} else if(fate_glass == REFLECTION) {
           else
             if(draw) ray1->DrawS(presecisce.x(), 10.0);
       }
+  }
   // Main test: ray and SiPM surface
   if(active->TestIntersection(&presecisce, *ray1)) {
       fate = hitExit;
   delete ray1;
   delete rayout;
   delete rayin;
   return fate;
+}
-//-----------------------------------------------------------------------------
 void CDetector::Draw(int width)
+{
   if(guide_on) {
       if( TMath::Abs(guide->getN1()-guide->getN2()) < MARGIN ) {
           if(_plateOn) plate->drawSkel(col_LG, width);
   if(glass_on) glass_circle->Draw(col_glass, width);
   //window_circle->Draw(col_glass, width);
   active->Draw(col_active, width);
+}
-//=================================================================================
 Plate::Plate(DetectorParameters& parameters)
+{
   TVector3 center = CENTER;
   const double b = parameters.getB();
   if( !result ) {
       // ce -NI- presecisca z vstopno
       fate = missed;
   } else if(result == REFLECTION) {
       if (dbg) printf("PLATE: reflected\n");
+      ray0 = ray1;
       fate = backreflected;
   } else {
       points[0] = ray1.GetR();
       //hfate->Fill(enter);
       //hin->Fill(vec1.y(), vec1.z());
                   if( sides[inters_i]->TestIntersection(&vec1, ray1) ) break;
+              }
+          }
           points[n_odb] = vec1;
           if(inters_i == 0) {
+              ray0 = ray1;
               fate = backreflected;
               break;} // backreflection
           propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
           if(inters_i == 5) { // successfull exit
   *n_points = n_odb+1;
   *out = ray0;
   return fate;
 };
-//=============================================================================================================================== <<<<<<<<

Subversion Repositories f9daq

(root)/lightguide/trunk/src/guide.cpp – Rev 72 → 73