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

 // vrne 0 ce ni presecisca; 1 ce se je lomil
 // 2 ce se je odbil; -2 ce se je absorbiral
 int CSurface::PropagateRay(CRay in, CRay *out, TVector3 *intersection)
+{
   if (dbg) printf("--- CSurface::PropagateRay ---\n");
-  double cosTi, cosTt, p_ref;
+  double cosTi; // incident ray angle
+  double cosTt; // transmited ray angle
   TVector3 intersect, transmit;
         if( !(GetIntersection(&intersect, in) == 1) )
                 return 0;
                 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);
         if(type == SURF_TOTAL) type = SURF_REFRA;
         switch(type){
                 // ----------------------------------------------------------------------------
                 // --------------- refraction from n1 to n2 -----------------------------------
                                 cosTtotal = 0.0;
                         if(dbg) printf("  cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
                         // reflection dependance on polarization missing
                         // reflection hardcoded to 0.96
-                        p_ref = rand.Uniform(0.0, 1.0);
                         if (dbg) printf("   reflection probability = %f\n", p_ref);
                         // If n1>n2 and theta>thetaCritical, total reflection
-                        if( (cosTi < cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
+                        if(cosTi < cosTtotal) {
                                 if(dbg) printf("  TOTAL\n");
                                 transmit = in.GetN() + sign_n*2*cosTi*n;
                                 if(dbg) {
                                         cosTN = TMath::Abs(transmit.Unit() * 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;
                                 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(cosTi<=cosTtotal) cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));
                                 //if(fresnel) {
-                                        r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse
+                                r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse
-                                        r_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel
+                                r_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel
-                                        if(dbg) printf("  r_te = %lf, r_tm = %lf\n", r_te, r_tm);
+                                if(dbg) printf("  r_te = %lf, r_tm = %lf\n", r_te, r_tm);
-                                        // transmited polarization
+                                // transmited polarization
-                                        v_tm_t = -v_te.Cross(transmit);
+                                v_tm_t = -v_te.Cross(transmit);
-                                        v_tm_t = v_tm_t.Unit();
+                                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) {
+                                if(dbg) {
                                                 printf("  v_tm_t = "); printv(v_tm_t);
                                                 printf("  pol_t = "); printv(pol_t);
-                                        }
+                                }
-          // Fresnel coefficients
+        // Fresnel coefficients
-                                        R_te = TMath::Power(r_te, 2);
+                                R_te = TMath::Power(r_te, 2);
-                                        R_tm = TMath::Power(r_tm, 2);
+                                R_tm = TMath::Power(r_tm, 2);
-                                        R_f = a_te*a_te*R_te + a_tm*a_tm*R_tm;
+                                R_f = a_te*a_te*R_te + a_tm*a_tm*R_tm;
-                                        if (dbg) printf("  R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
+                                if (dbg) printf("  R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
-                                }
+                        }
-                                if(p_ref >= R_f) { // se lomi
+                        if(p_ref >= R_f) { // se lomi
-                                  if (dbg) printf("   SURFACE REFRACTED. Return.\n");
+                                if (dbg) printf("   SURFACE REFRACTED. Return.\n");
                                         out->Set(intersect, transmit);
                                         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);
+                                }
                         //}
                         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;
                                 transmit = in.GetN() - 2*cosTi*n;
                                 out->Set(intersect, transmit);
                                 return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
                         } else { // se ne odbije
                                 transmit = in.GetN();
                                 out->Set(intersect, transmit);
                                 return ABSORBED;
+                        }
                         break;
                         if(inters_i == 0) {
                           fate = backreflected;
                           //hfate->Fill(backreflected);
                           break;
                         } // backreflection
                         // the passage is possible, test propagation
                         propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
                         if (dbg) printf("  GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
         double y_gap = parameters.getGap().Y();
         double z_gap = parameters.getGap().Z();
         // additional glass between at top of SiPM
         // example: epoxy n=1.60
-        double n4 = 1.60;
+        double n4 = 1.57;
         TVector3 plane_v[4];
         int nBins = nch + 1;
         double p_size = b/2.0;
         plane_v[0].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap - p_size);
         plane_v[1].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap + p_size);
+                }
         } else {
           if (dbg) printf("Detector: fate = hit or refracted");
                 ray1 = ray0;
                 if(draw) {
+                  //double epoxy = parameters->getGlassD();
-                        if(glass_on) ray1->Draw(center.x(), center.x() /*+ window_d*/);
+                        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;
-                        *ray0 = *ray1; ray1->SetColor(col_rgla);
+                        ray1->SetColor(col_rgla);
                         fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
-                        if(fate_glass == 1) {
+                        if(fate_glass == REFRACTION) {
                                 hglass->Fill(presecisce.y(), presecisce.z());
                                 if(draw) ray1->DrawS(presecisce.x(), 10.0);
-                                if(active->TestIntersection(&presecisce, *ray1)) {
+                                //if(active->TestIntersection(&presecisce, *ray1)) {
-                                        fate = hitExit;
+                                        //fate = hitExit;
-                                        hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
+                                        //hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
-                                        hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
+                                        //hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ);
-                                }
+                                //}
-                                if(detector->TestIntersection(&presecisce, *ray1))
+                                //if(detector->TestIntersection(&presecisce, *ray1))
-                                        hdetector->Fill(presecisce.y(), presecisce.z());
+                                        //hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
-                        } else if(fate_glass == 2) {
+                        //} else if(fate_glass == REFLECTION) {
+                                else
                                 if(draw) ray1->DrawS(presecisce.x(), 10.0);
+                        }
-                } else {
+                        }
                   // Main test: ray and SiPM surface
                         if(active->TestIntersection(&presecisce, *ray1)) {
                                 fate = hitExit;
                                 hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
                                 hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ);
+                        }
                         // If it is on the same plane as SiPM
                         if(detector->TestIntersection(&presecisce, *ray1))
                                 hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
-                }
+                //}
         //} else {
                 //if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
         //}
         *out = *ray1;

Subversion Repositories f9daq

(root)/lightguide/trunk/src/guide.cpp – Rev 70 → 71