385,7 → 385,8 |
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) ) |
437,6 → 438,9 |
} |
// ---------------------------------------------------------------------------- |
|
// reflection probability |
double p_ref = rand.Uniform(0.0, 1.0); |
|
if(type == SURF_TOTAL) type = SURF_REFRA; |
switch(type){ |
// ---------------------------------------------------------------------------- |
458,11 → 462,10 |
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; |
|
484,38 → 487,38 |
|
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_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel |
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); |
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; |
// 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) { |
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; |
// 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 (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"); |
if(p_ref >= R_f) { // se lomi |
if (dbg) printf(" SURFACE REFRACTED. Return.\n"); |
out->Set(intersect, transmit); |
out->SetPolarization(pol_t); |
return REFRACTION; |
919,7 → 922,7 |
|
// 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; |
1127,7 → 1130,8 |
if (dbg) printf("Detector: fate = hit or refracted"); |
ray1 = ray0; |
if(draw) { |
if(glass_on) ray1->Draw(center.x(), center.x() /*+ window_d*/); |
//double epoxy = parameters->getGlassD(); |
if(glass_on) ray1->Draw(center.x(), center.x() + glass_d); |
else ray1->DrawS(center.x(), 10.0); |
} |
} |
1134,22 → 1138,25 |
|
fate = missed; // zgresil aktivno povrsino |
if(glass_on) { |
*ray0 = *ray1; ray1->SetColor(col_rgla); |
*ray0 = *ray1; |
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)) { |
fate = hitExit; |
hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
hlaser->Fill((in.GetR()).y(), (in.GetR()).z()); |
} |
if(detector->TestIntersection(&presecisce, *ray1)) |
hdetector->Fill(presecisce.y(), presecisce.z()); |
} else if(fate_glass == 2) { |
//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(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); |
} |
} else { |
} |
|
// Main test: ray and SiPM surface |
if(active->TestIntersection(&presecisce, *ray1)) { |
fate = hitExit; |
1159,7 → 1166,7 |
// 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); |
//} |