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

         angle_ve1 = TMath::ACos(/*TMath::Abs*/( (e1.Unit()) * (vec.Unit()) ));
         angle_ve2 = TMath::ACos(/*TMath::Abs*/( (e2.Unit()) * (vec.Unit()) ));
         if(dbg)
+        {
-                //printf("angle_ve1 = %lf\n", angle_ve1*DEGREE);
+                printf("angle_ve1 = %lf\n", angle_ve1*DEGREE);
-                //printf("angle_ve2 = %lf\n", angle_ve2*DEGREE);
+                printf("angle_ve2 = %lf\n", angle_ve2*DEGREE);
-                //printf("angle_sum = %lf\n", (angle_ve1 + angle_ve2)*DEGREE);
+                printf("angle_sum = %lf\n", (angle_ve1 + angle_ve2)*DEGREE);
                 printf("  angle_r   = %lf\n", angle*DEGREE);
+        }
   bool difference = (MARGIN < TMath::Abs(angle - (angle_ve1 + angle_ve2)));
   if (dbg) printf("  MARGIN < Difference = %d\n", difference);
         if( !IsVectorIn(intersect) )
                 return 0;
         // --------------- Fresnel ----------------------------------------------------
         // R_f = a_te * R_te  +  a_tm * R_tm
+        // e - electrical/perependicular
+        // m - magnetic polarization/parallel
         double r_te=0;
         double r_tm=0;
-        double R_te=0;
+        double R_te=0; // s reflection coefficient
-        double R_tm=0;
+        double R_tm=0; // p reflection coefficient
         double R_f = 0.0;
-        double a_te = 0.0;
+        double a_te = 0.0; // s-wave amplitude, cos Alpha
-        double a_tm = 0.0;
+        double a_tm = 0.0; // p-wave amplitude, sin Alpha
-        TVector3 v_te; // polarization perpendicular to the plane of incidence
+        TVector3 v_te; // unit s-polarization vector
-        TVector3 v_tm; // inbound polarization parallel with the plane of incidence
+        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
         int sign_n; // sign of normal direction vs. inbound ray
         double cosTN; // debug
         if(fresnel) {
+          // p-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()
-                v_te = n.Cross(in.GetN()); v_te = v_te.Unit();
+                v_te = n.Cross(in.GetN());
+                v_te = v_te.Unit();
-                v_tm = -v_te.Cross(in.GetN()); v_tm = v_tm.Unit();
+                v_tm = -v_te.Cross(in.GetN());
+                v_tm = v_tm.Unit();
                 if(dbg) {
                         printf("  v_te = "); printv(v_te);
                         printf("  v_tm = "); printv(v_tm);
+                }
                 a_tm = TMath::Sqrt(1 - cosAf*cosAf);
                 if(dbg) printf("  a_te = %lf, a_tm = %lf\n", a_te, a_tm);
+        }
         // ----------------------------------------------------------------------------
-        if(type == 3) type = SURF_REFRA; //SURF_TOTAL -> SURF_REFRA
+        if(type == SURF_TOTAL) type = SURF_REFRA;
         switch(type){
                 // ----------------------------------------------------------------------------
                 // --------------- refraction from n1 to n2 -----------------------------------
                 // ----------------------------------------------------------------------------
                 case SURF_REFRA:
                         // 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);
-                        // Total reflection
+                        // If n1>n2 and theta>thetaCritical, total reflection
-                        /*
-                        if( (cosTi <= cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
+                        if( (cosTi < cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
                                 if(dbg) printf("  TOTAL\n");
                                 transmit = in.GetN() + sign_n*2*cosTi*n;
                                 if(dbg) {
                                         cosTN = TMath::Abs(transmit.Unit() * n);
                                 out->Set(intersect, transmit);
                                 pol_t = -in.GetP() + sign_n*2*cosTi*n;
                                 out->SetPolarization(pol_t);
                                 return REFLECTION;
-                        } else { */
+                        } 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);
                                         r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse
                                         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);
+                                        // 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;
                                         if(dbg) {
                                                 printf("  v_tm_t = "); printv(v_tm_t);
                                                 printf("  pol_t = "); printv(pol_t);
+                                        }
+          // Fresnel coefficients
                                         R_te = TMath::Power(r_te, 2);
                                         R_tm = TMath::Power(r_tm, 2);
                                         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(p_ref >= R_f) { // se lomi
                                   if (dbg) printf("   SURFACE REFRACTED. Return.\n");
                                         out->Set(intersect, transmit);
                                         out->SetPolarization(pol_t);
 //=================================================================================
 Guide::Guide(TVector3 center0, DetectorParameters &parameters)
+{
 double t;
+        TDatime now;
-        TDatime now; rand.SetSeed(now.Get());
+        rand.SetSeed(now.Get());
         center = center0;
         double b = parameters.getB();
         double a = parameters.getA();
         _d = parameters.getD();
         n1 = parameters.getN1();
         n2 = parameters.getN2();
         // if PlateOn, then n0 = n3 (optical grease), else = n1 (air)
-        double n0 = (parameters.getPlateOn() ? parameters.getN3(): n1);
+        //double n0 = (parameters.getPlateOn() ? parameters.getN3(): n1);
+        double n0 = (parameters.getPlateOn() ? n2 : n1);
         n3 = parameters.getN3();
         _r = c_reflectivity;
         int fresnel = parameters.getFresnel();
+  // light guide edges
         t = b/2.0;
-        vodnik_edge[0].SetXYZ(0.0, t,-t); vodnik_edge[1].SetXYZ(0.0, t, t);
+        vodnik_edge[0].SetXYZ(0.0, t,-t);
+        vodnik_edge[1].SetXYZ(0.0, t, t);
-        vodnik_edge[2].SetXYZ(0.0,-t, t); vodnik_edge[3].SetXYZ(0.0,-t,-t);
+        vodnik_edge[2].SetXYZ(0.0,-t, t);
+        vodnik_edge[3].SetXYZ(0.0,-t,-t);
         t = a/2.0;
-        vodnik_edge[4].SetXYZ(_d, t,-t); vodnik_edge[5].SetXYZ(_d, t, t);
+        vodnik_edge[4].SetXYZ(_d, t,-t);
+        vodnik_edge[5].SetXYZ(_d, t, t);
-        vodnik_edge[6].SetXYZ(_d,-t, t); vodnik_edge[7].SetXYZ(_d,-t,-t);
+        vodnik_edge[6].SetXYZ(_d,-t, t);
+        vodnik_edge[7].SetXYZ(_d,-t,-t);
         for(int i = 0; i<8; i++) vodnik_edge[i] += center;
+        // light guide surfaces
         s_side[0] = new CSurface(SURF_REFRA, vodnik_edge, n0, n2, _r);
         s_side[0]->FlipN();
         s_side[1] = new CSurface(SURF_REFRA, vodnik_edge[3], vodnik_edge[2], vodnik_edge[6], vodnik_edge[7], n2, n1, _r);
         s_side[2] = new CSurface(SURF_REFRA, vodnik_edge[2], vodnik_edge[1], vodnik_edge[5], vodnik_edge[6], n2, n1, _r);
         s_side[3] = new CSurface(SURF_REFRA, vodnik_edge[1], vodnik_edge[0], vodnik_edge[4], vodnik_edge[5], n2, n1, _r);
         s_side[4] = new CSurface(SURF_REFRA, vodnik_edge[0], vodnik_edge[3], vodnik_edge[7], vodnik_edge[4], n2, n1, _r);
-        s_side[5] = new CSurface(SURF_REFRA, &vodnik_edge[4], n2, n3, _r);
+        s_side[5] = new CSurface(SURF_REFRA, &vodnik_edge[4], n2, n3, _r); // n3 - ref ind at the exit, grease, air, epoxy
         s_side[5]->FlipN();
         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;
         hfate = new TH1F("hfate", "Ray fate", 8, -3.5, 4.5);
         (hfate->GetXaxis())->SetBinLabel(1, "Back Ref");
         (hfate->GetXaxis())->SetBinLabel(2, "No Ref");
         (hfate->GetXaxis())->SetBinLabel(3, "Refrac");
           fate = backreflected;
           //hfate->Fill(-3);
         } else {
           if (dbg) printf("  GUIDE: ray entered\n");
                 points[0] = ray1.GetR();
-                hfate->Fill(2); // enter
+                hfate->Fill(enter); // enter
                 hin->Fill(vec1.y(), vec1.z());
                 if (dbg) printf("  GUIDE: n_odb = %d\n", n_odb);
                 while (n_odb++ < MAX_REFLECTIONS) {
                   if (dbg) printf("  GUIDE: Boundary test: %d\n",n_odb);
         if(TMath::Abs(den) < MARGIN) {
                 if(TMath::Abs(num) < MARGIN)
                         return 0;
                 else
                         return 0;
         }
         t = num / den;
         if(dbg) printf("t = %.4lf | ", t);
         tmp = ray.GetR();
         tmp -= t*ray.GetN();
         *vec = tmp;
         else x_offset = center.x() - d;
         //guide = new CVodnik(center, SiPM, M, d, type_in, type_side, type_out, n1, n2, n3, reflectivity, fresnel, absorption, A);
         double b = parameters.getB();
-        double n1 = parameters.getN1();
+        //double n1 = parameters.getN1();
-        double n2 = parameters.getN2();
+        //double n2 = parameters.getN2();
-        //double n3 = parameters.getN3();
+        double n3 = parameters.getN3();
-        double reflectivity = c_reflectivity; // for faster simulation, not using Fresnel eqs.
+        double reflectivity = c_reflectivity;
         double x_gap = parameters.getGap().X();
         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;
         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);
         plane_v[2].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap + p_size);
         plane_v[3].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap - p_size);
-        glass = new CSurface(SURF_REFRA, plane_v, n2, n1, reflectivity); glass->FlipN();
+        glass = new CSurface(SURF_REFRA, plane_v, n3, n4, reflectivity);
+        glass->FlipN();
+        // additional circular glass between LG and SiPM
         glass_circle = new CPlaneR(TVector3(x_offset+d+glass_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), b);
         hglass = (TH2F*)gROOT->FindObject("hglass"); if(hglass) delete hglass;
         hglass = new TH2F("hglass", "Hits glass",
                                           nBins, y_gap - p_size, y_gap + p_size,
                   nBins, z_gap - p_size, z_gap + p_size);
+        // SiPM active surface
         p_size = parameters.getActive()/2.0;
         //cout<<"SiPM active length "<<detectorActive<<endl;
         //p_size = 1.0/2.0;
         plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size);
         plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size);
         plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size);
         plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size);
         active = new CPlane4(plane_v);
         hactive = (TH2F*)gROOT->FindObject("hactive"); if(hactive) delete hactive;
         //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);
         hactive = new TH2F("hactive", "Active area hits", nBins, y_gap - p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ);
         p_size = b/2.0;
         //p_size = 2.5;
         //p_size = M*0.6;
         hlaser = (TH2F*)gROOT->FindObject("hlaser"); if(hlaser) delete hlaser;
         hlaser = new TH2F("hlaser", ";x [mm]; y [mm]", nBins, -p_size+offsetY, p_size+offsetY, nBins, -p_size+offsetZ, p_size+offsetZ);
+        // collection surface in SiPM plane
         p_size = 1.4*b/2.0;
         plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size);
         plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size);
         plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size);
         plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size);
+{
   if (dbg) printf("--- Detector::Propagate ---\n");
         //CRay *ray0 = new CRay; ray0->Set(in.GetR(), in.GetN()); ray0->SetColor(col_in);
         CRay *rayin = new CRay(in);
         rayin->SetColor(col_in);
-        CRay *rayout = new CRay;
+        CRay *rayout = new CRay(in);
+        rayout->SetColor(col_in);
         const int max_n_points = guide->GetMAXODB() + 2;
         TVector3 pointsPlate[max_n_points];
         //TVector3 intersection;
         Fate fatePlate;
         // check if the ray should be reflected??
         if(_plateOn) {
           fatePlate = plate->propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
           if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
           if(draw) {
             if(fatePlate == missed) {
               rayout->SetColor(col_in);
               rayout->DrawS(center.x() - _plateWidth, -10.0);
+              }
+            else if(fatePlate == backreflected){
+              if (dbg) printf("Backreflected at plate!\n");
+              }
-              else {
+            else {
                 int p_i;
                 for(p_i = 0; p_i < nPointsPlate-1; p_i++) {
                                           line3d->SetPoint(0, pointsPlate[p_i].x(), pointsPlate[p_i].y(), pointsPlate[p_i].z());
                                           line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z());
                                           line3d->DrawClone();
+                                  }
+                                  rayout->DrawS(pointsPlate[p_i].x(), -0.1);
-                                  if(fatePlate != noreflection) {
+                                  if(fatePlate == noreflection) { // lost on plate side
-                                          //rayout->SetColor(7);
+                                    rayout->SetColor(col_out);
-                                          rayout->DrawS(pointsPlate[p_i].x(), -0.1);
+                                          rayout->DrawS(pointsPlate[p_i].x(), 10.0);
                                   }
-                          }
+                        }
-            }
+        }
           if(! (fatePlate == hitExit or fatePlate == refracted) ) {
+                        guide->GetHFate()->Fill(rays);
-                  if (dbg) printf("CDetector::propagate Simulated ray missed the entry surface!\n");
+                  if (dbg)printf("CDetector::propagate Simulated ray missed the entry surface!\n");
+                  if (fatePlate == backreflected)
+                    guide->GetHFate()->Fill(fatePlate); // reflected back
+                  else
+                    guide->GetHFate()->Fill(noreflection); //lost on plate side
                         return fatePlate;
+                }
-                // missing: if refracted at plate sides
+    //Ray hits light guide
-                //if (fatePlate == refracted) return fatePlate;
                 histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
+         }
          else {
-          rayout = rayin;
+          //rayout = rayin;
           if(draw) rayout->DrawS(center.x(), -10.0);
+          }
         // If the ray is not reflected in the plate
         // Draw the light guide and propagate the ray through
         int n_points;
         int fate_glass;
         CRay *ray0 = new CRay(*rayout);
         // delete rayout; -> creates dangling reference when tries to delete ray0!
-        delete rayin;
+        //delete rayin; -> delete rayout!
         CRay *ray1 = new CRay;
         fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
         if (dbg) {
           if (fate == backreflected) printf("DETECTOR::backreflected\n");
                 if(! (fate == hitExit or fate == refracted) ) {
                   if (dbg) printf("Detector: fate != hit, refracted\n");
                         *out = *ray1;
+                        delete ray0;
+                        delete ray1;
+                        delete rayout;
+                        delete rayin;
                         return fate;
+                }
         } else {
           if (dbg) printf("Detector: fate = hit or refracted");
                 ray1 = ray0;
                         if(glass_on) ray1->Draw(center.x(), center.x() /*+ window_d*/);
                         else ray1->DrawS(center.x(), 10.0);
+                }
+        }
-/*
-        TVector3 pres_wind;
-        fate = window_circle->TestIntersection(&pres_wind, *ray1);
-        if(fate == 1) {
-                hwindow->Fill(pres_wind.y(), pres_wind.z());
-                if(!guide_on) {
-                        window->PropagateRay(*ray0, ray1, &presecisce);
-                        if(draw) ray1->Draw(center.x() + window_d, center.x() + glass_d);
-                        *ray0 = *ray1;
-                }
-        */
                 fate = missed; // zgresil aktivno povrsino
                 if(glass_on) {
                         *ray0 = *ray1; ray1->SetColor(col_rgla);
                         fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
                         if(fate_glass == 1) {
                                 hglass->Fill(presecisce.y(), presecisce.z());
                                 if(draw) ray1->DrawS(presecisce.x(), 10.0);
                                 if(active->TestIntersection(&presecisce, *ray1)) {
                                         fate = hitExit;
                                         hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
                                         hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
+                                }
                   // 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 {
         *out = *ray1;
         delete ray0;
         delete ray1;
         delete rayout;
+        delete rayin;
         return fate;
+}
 //-----------------------------------------------------------------------------
 void CDetector::Draw(int width)
+{
+{
   TVector3 center = CENTER;
   const double b = parameters.getB();
   const double n1 = parameters.getN1();
         const double n2 = parameters.getN2();
-        const double n3 = parameters.getN3();
-        const double reflectivity = c_reflectivity;
         const double t = b/2.;
         const double plateWidth = parameters.getPlateWidth();
         center.SetX( CENTER.X() - plateWidth );
-        plate_edge[0].SetXYZ(0.0, t,-t); plate_edge[1].SetXYZ(0.0, t, t);
+        plate_edge[0].SetXYZ(0.0, t,-t);
+        plate_edge[1].SetXYZ(0.0, t, t);
+        plate_edge[2].SetXYZ(0.0,-t, t);
+        plate_edge[3].SetXYZ(0.0,-t,-t);
-        plate_edge[2].SetXYZ(0.0,-t, t); plate_edge[3].SetXYZ(0.0,-t,-t);
+        plate_edge[4].SetXYZ(plateWidth, t,-t);
-        plate_edge[4].SetXYZ(plateWidth, t,-t); plate_edge[5].SetXYZ(plateWidth, t, t);
+        plate_edge[5].SetXYZ(plateWidth, t, t);
-        plate_edge[6].SetXYZ(plateWidth,-t, t); plate_edge[7].SetXYZ(plateWidth,-t,-t);
+        plate_edge[6].SetXYZ(plateWidth,-t, t);
+        plate_edge[7].SetXYZ(plateWidth,-t,-t);
         for(int i = 0; i<8; i++) plate_edge[i] += center;
-        sides[0] = new CSurface(SURF_REFRA, plate_edge, n1, n2, reflectivity);
+        sides[0] = new CSurface(SURF_REFRA, plate_edge, n1, n2, c_reflectivity);
         sides[0]->FlipN();
-        sides[1] = new CSurface(SURF_REFRA, plate_edge[3], plate_edge[2], plate_edge[6], plate_edge[7], n2, n2, reflectivity);
+        sides[1] = new CSurface(SURF_REFRA, plate_edge[3], plate_edge[2], plate_edge[6], plate_edge[7], n2, n2, c_reflectivity);
-        sides[2] = new CSurface(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, reflectivity);
+        sides[2] = new CSurface(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, c_reflectivity);
-        sides[3] = new CSurface(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, reflectivity);
+        sides[3] = new CSurface(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, c_reflectivity);
-        sides[4] = new CSurface(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, reflectivity);
+        sides[4] = new CSurface(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, c_reflectivity);
-        sides[5] = new CSurface(SURF_REFRA, &plate_edge[4], n2, n3, reflectivity);
+        sides[5] = new CSurface(SURF_REFRA, &plate_edge[4], n2, n2, c_reflectivity);
         sides[5]->FlipN();
         for(int i=0; i<6; i++) sides[i]->SetFresnel(1);
+}
         } else if(result == REFLECTION) {
           if (dbg) printf("PLATE: reflected\n");
           fate = backreflected;
         } else {
                 points[0] = ray1.GetR();
-                //hfate->Fill(2);
+                //hfate->Fill(enter);
                 //hin->Fill(vec1.y(), vec1.z());
                 while (n_odb++ < MAX_REFLECTIONS) {
                         ray0 = ray1;
                         vec0 = vec1;
                         propagation = 11;
                                 points[n_odb] = vec1;
                                 ray0 = ray1;
                                 break;
+                        }
                         if(propagation == 1) {
-                          fate = refracted; //at side?
+                          fate = noreflection; //at side
                           n_odb++;
                           points[n_odb] = vec1;
                           ray0 = ray1;
                           break;} // no total reflection when should be

Subversion Repositories f9daq

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