Subversion Repositories f9daq

Compare Revisions

Regard whitespace Rev 72 → Rev 73

/lightguide/trunk/include/guide.h
39,14 → 39,14
public:
CRay() :
r(TVector3(0,0,0)),
n(TVector3(1,0,0)),
k(TVector3(1,0,0)),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
CRay(TVector3 r0, TVector3 n0) :
CRay(TVector3 r0, TVector3 k0) :
r(r0),
n(n0.Unit()),
k(k0.Unit()),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
53,21 → 53,24
{};
CRay(double x0, double y0, double z0, double l0, double m0, double n0) :
r(TVector3(x0,y0,z0)),
n(TVector3(l0,m0,n0).Unit()),
k(TVector3(l0,m0,n0).Unit()),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
 
void Set(TVector3 r0, TVector3 n0);
void Set(TVector3 r0, TVector3 k0);
//void Set(double x0, double y0, double z0, double l0, double m0, double n0);
void SetColor(int c){color = c;};
void SetPolarization(TVector3 p0) {p = p0.Unit();};
void setPolarization(TVector3 p0) {
p = p0.Unit();
if (p.Dot(k)>1e-3) printf("*** ERROR in CRay: E has component || with k\n");
};
 
//inline CRay & operator = (const CRay &);
 
TVector3 GetR() const {return r;};
TVector3 GetN() const {return n;};
TVector3 GetK() const {return k;};
TVector3 GetP() const {return p;};
 
void Print();
75,13 → 78,13
void Draw(double x_from, double x_to);
void DrawS(double x_from, double t);
 
//r = position, k = unit wave vector, p = polarization
private:
TVector3 r;
TVector3 n;
TVector3 p; //r = point on line, n = normal, p = polarization
TVector3 k;
TVector3 p;
int color;
};
//=================================================================================
 
 
//=================================================================================
141,7 → 144,6
TVector3 n, center;
double _r;
};
//=================================================================================
 
//=================================================================================
// ravna opticna povrsina: refractor, zrcalo ali povrsina s totalnim odbojem
188,7 → 190,6
 
int fresnel; // ali naj uposteva Fresnelove enacbe
};
//=================================================================================
 
 
//=================================================================================
363,9 → 364,7
TH1F *hfate, *hnodb_all, *hnodb_exit;
TH2F *hin, *hout;
};
//=================================================================================
 
//=============================================================================================================================== <<<<<<<<
 
 
class Plate
385,7 → 384,6
CSurface *sides[6];
};
 
// ================================================================================
 
class CDetector
{
511,5 → 509,4
 
 
 
//=================================================================================
#endif
/lightguide/trunk/src/userFunctions.cpp
111,11 → 111,6
 
//------------------------------------------------------------------------------------------
 
 
TVector3 p_pol0(0.0, 0.0, 1.0);
void SetPol(double x, double y, double z)
{p_pol0.SetXYZ(x,y,z);}
 
// Test function
// creates an optical boundary surface
// shows the propagation of light ray
123,11 → 118,11
// and surface normal (blue)
void PolTest(double theta = 0.0)
{
int p_type = 1;
int p_type = SURF_REFRA;
double p_n1 = parameters.getN1();
double p_n2 = parameters.getN2();
theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6;
TVector3 p_pol;
theta = 3.141593*theta/180.0;
if(theta < 1e-6) theta = 1e-6;
Init();
134,31 → 129,33
double cx = 0.0;
TVector3 vodnik_edge[4];
double t = 3.0;
vodnik_edge[0].SetXYZ(cx, t,-t); vodnik_edge[1].SetXYZ(cx, t, t);
vodnik_edge[2].SetXYZ(cx,-t, t); vodnik_edge[3].SetXYZ(cx,-t,-t);
/*
#define SURF_DUMMY 0
#define SURF_REFRA 1
#define SURF_REFLE 2
#define SURF_TOTAL 3
#define SURF_IMPER 4
*/
CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN();
vodnik_edge[0].SetXYZ(cx, t,-t);
vodnik_edge[1].SetXYZ(cx, t, t);
vodnik_edge[2].SetXYZ(cx,-t, t);
vodnik_edge[3].SetXYZ(cx,-t,-t);
 
CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96);
surf->FlipN();
surf->SetFresnel(1);
surf->Draw();
CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
ray->SetColor(kBlack);
ray->SetColor(kRed);
//p_pol = rotatey(p_pol0, -theta);
p_pol = p_pol0; p_pol.RotateY(-theta);
TVector3 sE(0,1,0);
TVector3 pE(0,0,1);
TVector3 p_pol = pE;
p_pol.RotateY(-theta);
//printf("p_pol = "); printv(p_pol);
ray->SetPolarization(p_pol);
ray->setPolarization(p_pol);
ray->DrawS(cx, -5.0);
CRay *out = new CRay; out->SetColor(kRed);
CRay *out = new CRay;
out->SetColor(kBlack);
TVector3 *inters = new TVector3;
surf->PropagateRay(*ray, out, inters);
printf(" n1 = %f, n2 = %f\n", p_n1, p_n2);
//if(fate == 1) out->DrawS(cx, 5.0);
out->DrawS(cx, 5.0);
166,11 → 163,13
incidentPolarization->SetColor(kGreen);
incidentPolarization->Set(ray->GetR(), p_pol);
incidentPolarization->DrawS(cx, 1.0);
printf(" GREEN: polarization vector\n");
CRay *surfaceNormal = new CRay;
surfaceNormal->SetColor(kBlue);
surfaceNormal->Set(ray->GetR(), surf->GetN());
surfaceNormal->DrawS(cx, 1.0);
printf(" BLUE: surface normal vector\n");
}
 
void ptt()
/lightguide/trunk/src/raySimulator.cpp
60,7 → 60,7
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
*/
//-----------------------------------------------------------------------------
void Init(double rsys_scale = 7.0)
void Init(double rsys_scale = 10.0)
{
RTSetStyle(gStyle);
gStyle->SetOptStat(10);
70,7 → 70,8
double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
c3dview = (TCanvas*)gROOT->FindObject("c3dview");
if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723); c3dview->SetFillColor(0);}
if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);
c3dview->SetFillColor(0);}
else c3dview->Clear();
TView3D *view = new TView3D(1, rsys0, rsys1);
177,9 → 178,9
//-----------------------------------------------------------------------------
// en zarek
double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
double Single(CDetector *detector, DetectorParameters& parameters,
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = Pi()*theta/180.0;
if(theta < 1e-6) theta = 1e-6;
phi = phi*Pi()/180.0;
191,11 → 192,32
if(show_3d) detector->Draw(draw_width);
CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi),
TMath::Sin(theta)*TMath::Cos(phi));
// Set z-polarization == vertical
TVector3 polarization(0,0,1);
polarization.RotateY(-theta);
polarization.RotateZ(phi);
if (polarization.Dot(ray0->GetK()) > 1e-5) printf("ERROR: pol not perep\n");
ray0->setPolarization(polarization);
CRay *ray1 = new CRay;
detector->Propagate(*ray0, ray1, show_3d);
CRay *incidentPolarization = new CRay;
incidentPolarization->SetColor(kGreen);
TVector3 drawPosition = ray0->GetR();
drawPosition.SetX(drawPosition.X() - 5);
incidentPolarization->Set(drawPosition, polarization);
incidentPolarization->DrawS(drawPosition.X(), 1);
 
TVector3 outPolarization = ray1->GetP();
drawPosition = ray1->GetR();
drawPosition.SetX(drawPosition.X() + 5);
CRay* rayPol = new CRay(drawPosition, outPolarization);
rayPol->SetColor(kBlack);
rayPol->DrawS(drawPosition.X(), 1);
 
delete ray0;
delete ray1;
203,7 → 225,8
}
//-----------------------------------------------------------------------------
// zarki, razporejeni v mrezi
double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
double Grid(CDetector *detector, DetectorParameters& parameters,
int NN = 10, double theta = 0.0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = Pi()*theta/180.0;
224,6 → 247,9
TMath::Cos(theta),
0.0,
TMath::Sin(theta));
TVector3 polarization(0, 1, 0);
polarization.RotateY(-theta);
ray0->setPolarization(polarization);
CRay *ray1 = new CRay;
if(i == (NN/2))
detector->Propagate(*ray0, ray1, show_3d);
247,7 → 273,8
//-----------------------------------------------------------------------------
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
// vsi pod kotom (theta, phi)
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phi, int show_rays)
double RandYZ(CDetector *detector, DetectorParameters& parameters,
int NN, double theta, double phi, int show_rays)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = theta*3.14159265358979312/180.0;
277,6 → 304,10
double rn = TMath::Sin(theta)*TMath::Cos(phi);
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
TVector3 polarization(0, 0, 1);
polarization.RotateY(-theta);
polarization.RotateZ(phi);
ray0->setPolarization(polarization);
CRay *ray1 = new CRay;
if(i < show_rays)
283,8 → 314,8
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
//delete ray0;
//delete ray1;
delete ray1;
delete ray0;
}
double acceptance = 0.0;
304,7 → 335,8
// in cos theta and phi uniformly:
// theta [0, theta]
// phi [0,360]
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
double RandIso(CDetector *detector, DetectorParameters& parameters,
int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
double pi = 3.14159265358979312;
405,7 → 437,8
 
// Beamtest distribution
// with fixed theta and phi in interval phiMin, phiMax
double beamtest(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
double beamtest(CDetector *detector, DetectorParameters& parameters,
int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
{
double pi = 3.14159265358979312;
theta *= pi/180.0;
418,9 → 451,6
rand.SetSeed(now.Get());
double rx, ry, rz, rl, rm, rn;
double rphi;
//double rtheta;
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
//double theta_min_rad = TMath::Cos(theta);
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
461,10 → 491,17
//htheta->Fill(rtheta);
//hcostheta->Fill( TMath::Cos(rtheta) );
hphi->Fill(rphi);
hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
hl->Fill(rl);
hm->Fill(rm);
hn->Fill(rn);
}
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
// inital polarizaton
TVector3 polarization(0, 0, 1);
polarization.RotateY(-theta);
polarization.RotateX(rphi);
ray0->setPolarization(polarization);
CRay *ray1 = new CRay;
if(i < show_rays) {
/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;
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;
};
//=============================================================================================================================== <<<<<<<<