Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 72 → Rev 73

/lightguide/trunk/src/guide.cpp
20,12 → 20,10
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)
42,14 → 40,12
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;
56,54 → 52,50
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();
A2 = (x_to - r.x())/n.x();
A1 = (x_from - r.x())/k.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(1, x_to, A2*n.y()+r.y(), A2*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*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(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.z());
line3d->SetPoint(0, x_from, A1*k.y()+r.y(), A1*k.z()+r.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),
117,7 → 109,6
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);
209,7 → 200,7
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);
 
229,7 → 220,7
t = num / den;
 
tmp = ray.GetR();
tmp -= t*ray.GetN();
tmp -= t*ray.GetK();
*vec = tmp;
return 1;
}
275,7 → 266,6
 
return 1;
}
//-----------------------------------------------------------------------------
int CPlane4::TestIntersection(CRay in)
{
TVector3 tmp;
286,7 → 276,6
 
return 0;
}
//-----------------------------------------------------------------------------
int CPlane4::TestIntersection(TVector3 *vec, CRay in)
{
TVector3 tmp;
299,7 → 288,6
 
return 0;
}
//-----------------------------------------------------------------------------
void CPlane4::Print()
{
printf("--- CPlane4::Print() ---\n");
309,7 → 297,6
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);
320,10 → 307,8
 
line3d->Draw();
}
//=================================================================================
 
 
//=================================================================================
CSurface::CSurface(int type0):
type(type0)
{
342,7 → 327,6
 
SetFresnel();
}
//-----------------------------------------------------------------------------
CSurface::CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
{
TDatime now;
355,7 → 339,6
 
SetFresnel();
}
//-----------------------------------------------------------------------------
CSurface::CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
{
TDatime now;
368,7 → 351,6
 
SetFresnel();
}
//-----------------------------------------------------------------------------
void CSurface::SetIndex(double n10, double n20)
{
n1 = n10; n2 = n20; n1_n2 = n1/n2;
410,19 → 392,22
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
 
if(fresnel) {
// p-polarization unit vector v_te
// Decomposition of incident polarization vector
// using unit vectors v_tm & v_te
// in a_tm and a_te components
//if(fresnel) {
// 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()
v_te = n.Cross(in.GetN());
// k in this notation is in.GetK()
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);
430,12 → 415,13
}
 
double cosAf = v_te * in.GetP();
if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, TMath::ACos(cosAf)*DEGREE);
double alpha = acos(cosAf);
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
447,7 → 433,7
// --------------- 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);
467,7 → 453,7
// 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);
475,9 → 461,24
}
out->Set(intersect, transmit);
 
// Shift?
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
// 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;
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 = a_tm*v_tm_t + a_te*v_te;
if (dbg) printv(pol_t);
out->setPolarization(pol_t);
 
return REFLECTION;
} else {
// reflection or refraction according to Fresnel equations
486,7 → 487,7
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;
503,7 → 504,7
// 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);
521,14 → 522,17
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;
out->SetPolarization(pol_t);
//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);
return REFLECTION;
}
 
541,12 → 545,12
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;
cosTi = in.GetK() * 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;
}
556,17 → 560,17
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;
}
579,10 → 583,7
 
return REFRACTION;
}
//=================================================================================
 
 
//=================================================================================
Guide::Guide(TVector3 center0, DetectorParameters &parameters) :
_d(parameters.getD()),
_n1(parameters.getN1()),
640,7 → 641,7
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);
 
750,7 → 751,7
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) {
806,17 → 807,14
*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);
828,9 → 826,8
line3d->DrawClone();
}
}
//=================================================================================
 
//=================================================================================
 
int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
{
double num, den; //stevec, imenovalec
842,7 → 839,7
 
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);
 
858,7 → 855,7
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);}
869,7 → 866,7
else
return 0;
}
//-----------------------------------------------------------------------------
 
void CPlaneR::Draw(int color, int width)
{
const int NN = 32;
888,10 → 885,8
}
arc->Draw();
}
//=================================================================================
 
 
//=================================================================================
CDetector::CDetector(TVector3 center0, DetectorParameters& parameters) :
center(center0),
glass_on(parameters.getGlassOn()),
1048,6 → 1043,8
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 {
1110,8 → 1107,8
if(fate == missed) {
if (dbg) printf("Detector: fate=missed\n");
TVector3 r = ray1->GetR();
TVector3 n = ray1->GetN();
ray1->Set(r,n);
TVector3 k = ray1->GetK();
ray1->Set(r,k);
ray1->DrawS(center.x(), 10.0);
} else {
for(p_i = 0; p_i < n_points-1; p_i++) {
1192,7 → 1189,7
delete rayin;
return fate;
}
//-----------------------------------------------------------------------------
 
void CDetector::Draw(int width)
{
if(guide_on) {
1210,8 → 1207,8
//window_circle->Draw(col_glass, width);
active->Draw(col_active, width);
}
//=================================================================================
 
 
Plate::Plate(DetectorParameters& parameters)
{
TVector3 center = CENTER;
1284,6 → 1281,7
fate = missed;
} else if(result == REFLECTION) {
if (dbg) printf("PLATE: reflected\n");
ray0 = ray1;
fate = backreflected;
} else {
points[0] = ray1.GetR();
1300,6 → 1298,7
}
points[n_odb] = vec1;
if(inters_i == 0) {
ray0 = ray1;
fate = backreflected;
break;} // backreflection
 
1333,6 → 1332,5
*out = ray0;
return fate;
};
//=============================================================================================================================== <<<<<<<<