Subversion Repositories f9daq

Compare Revisions

Regard whitespace Rev 73 → Rev 72

/lightguide/trunk/include/guide.h
39,14 → 39,14
public:
CRay() :
r(TVector3(0,0,0)),
k(TVector3(1,0,0)),
n(TVector3(1,0,0)),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
CRay(TVector3 r0, TVector3 k0) :
CRay(TVector3 r0, TVector3 n0) :
r(r0),
k(k0.Unit()),
n(n0.Unit()),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
53,24 → 53,21
{};
CRay(double x0, double y0, double z0, double l0, double m0, double n0) :
r(TVector3(x0,y0,z0)),
k(TVector3(l0,m0,n0).Unit()),
n(TVector3(l0,m0,n0).Unit()),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
 
void Set(TVector3 r0, TVector3 k0);
void Set(TVector3 r0, TVector3 n0);
//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();
if (p.Dot(k)>1e-3) printf("*** ERROR in CRay: E has component || with k\n");
};
void SetPolarization(TVector3 p0) {p = p0.Unit();};
 
//inline CRay & operator = (const CRay &);
 
TVector3 GetR() const {return r;};
TVector3 GetK() const {return k;};
TVector3 GetN() const {return n;};
TVector3 GetP() const {return p;};
 
void Print();
78,13 → 75,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 k;
TVector3 p;
TVector3 n;
TVector3 p; //r = point on line, n = normal, p = polarization
int color;
};
//=================================================================================
 
 
//=================================================================================
144,6 → 141,7
TVector3 n, center;
double _r;
};
//=================================================================================
 
//=================================================================================
// ravna opticna povrsina: refractor, zrcalo ali povrsina s totalnim odbojem
190,6 → 188,7
 
int fresnel; // ali naj uposteva Fresnelove enacbe
};
//=================================================================================
 
 
//=================================================================================
364,7 → 363,9
TH1F *hfate, *hnodb_all, *hnodb_exit;
TH2F *hin, *hout;
};
//=================================================================================
 
//=============================================================================================================================== <<<<<<<<
 
 
class Plate
384,6 → 385,7
CSurface *sides[6];
};
 
// ================================================================================
 
class CDetector
{
509,4 → 511,5
 
 
 
//=================================================================================
#endif
/lightguide/trunk/src/userFunctions.cpp
111,6 → 111,11
 
//------------------------------------------------------------------------------------------
 
 
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
118,11 → 123,11
// and surface normal (blue)
void PolTest(double theta = 0.0)
{
int p_type = SURF_REFRA;
int p_type = 1;
double p_n1 = parameters.getN1();
double p_n2 = parameters.getN2();
theta = 3.141593*theta/180.0;
if(theta < 1e-6) theta = 1e-6;
theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6;
TVector3 p_pol;
 
Init();
 
129,33 → 134,31
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);
 
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);
/*
#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();
surf->SetFresnel(1);
surf->Draw();
 
CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
ray->SetColor(kRed);
ray->SetColor(kBlack);
//p_pol = rotatey(p_pol0, -theta);
TVector3 sE(0,1,0);
TVector3 pE(0,0,1);
TVector3 p_pol = pE;
p_pol.RotateY(-theta);
p_pol = p_pol0; 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(kBlack);
CRay *out = new CRay; out->SetColor(kRed);
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);
 
163,13 → 166,11
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/guide.cpp
20,10 → 20,12
if(in >= 0.0) return 1;
else return -1;
}
//=================================================================================
 
void CRay::Set(TVector3 r0, TVector3 k0)
//-----------------------------------------------------------------------------
void CRay::Set(TVector3 r0, TVector3 n0)
{
r = r0; k = k0.Unit();
r = r0; n = n0.Unit();
}
//-----------------------------------------------------------------------------
//void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0)
40,12 → 42,14
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(), k.x(), k.y(), k.z());
r.x(), r.y(), r.z(), n.x(), n.y(), n.z());
}
//-----------------------------------------------------------------------------
void CRay::Draw()
{
double t = 50.0;
52,50 → 56,54
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*k.x(), k.y() + t*k.y(), r.z() + t*k.z());
line3d->SetPoint(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.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(k.x() < MARGIN) {
if(n.x() < MARGIN) {
A1 = A2 = 0.0;
} else {
A1 = (x_from - r.x())/k.x();
A2 = (x_to - r.x())/k.x();
A1 = (x_from - r.x())/n.x();
A2 = (x_to - r.x())/n.x();
}
 
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->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->SetLineWidth(1);
line3d->SetLineColor(color);
 
line3d->Draw();
}
//-----------------------------------------------------------------------------
void CRay::DrawS(double x_from, double t)
{
double A1;
TPolyLine3D *line3d = new TPolyLine3D(2);
 
if(k.x() < MARGIN)
if(n.x() < MARGIN)
A1 = 0.0;
else
A1 = (x_from - r.x())/k.x();
A1 = (x_from - r.x())/n.x();
 
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->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->SetLineWidth(1);
line3d->SetLineColor(color);
 
line3d->Draw();
}
//=================================================================================
 
 
//=================================================================================
CPlane4::CPlane4() :
n(TVector3(1.0, 0.0, 0.0)),
A(0),
109,6 → 117,7
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);
200,7 → 209,7
N.SetXYZ(A,B,C);
 
num = N*ray.GetR() + D;
den = N*ray.GetK();
den = N*ray.GetN();
 
if (dbg) printf("t = %6.3lf / %6.3lf = %6.3lf\n", num, den, num/den);
 
220,7 → 229,7
t = num / den;
 
tmp = ray.GetR();
tmp -= t*ray.GetK();
tmp -= t*ray.GetN();
*vec = tmp;
return 1;
}
266,6 → 275,7
 
return 1;
}
//-----------------------------------------------------------------------------
int CPlane4::TestIntersection(CRay in)
{
TVector3 tmp;
276,6 → 286,7
 
return 0;
}
//-----------------------------------------------------------------------------
int CPlane4::TestIntersection(TVector3 *vec, CRay in)
{
TVector3 tmp;
288,6 → 299,7
 
return 0;
}
//-----------------------------------------------------------------------------
void CPlane4::Print()
{
printf("--- CPlane4::Print() ---\n");
297,6 → 309,7
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);
307,8 → 320,10
 
line3d->Draw();
}
//=================================================================================
 
 
//=================================================================================
CSurface::CSurface(int type0):
type(type0)
{
327,6 → 342,7
 
SetFresnel();
}
//-----------------------------------------------------------------------------
CSurface::CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
{
TDatime now;
339,6 → 355,7
 
SetFresnel();
}
//-----------------------------------------------------------------------------
CSurface::CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
{
TDatime now;
351,6 → 368,7
 
SetFresnel();
}
//-----------------------------------------------------------------------------
void CSurface::SetIndex(double n10, double n20)
{
n1 = n10; n2 = n20; n1_n2 = n1/n2;
392,22 → 410,19
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(); // incident polarization
TVector3 pol_t = in.GetP(); // transmited polarization
int sign_n; // sign of normal direction vs. inbound ray
double cosTN; // debug
 
// 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
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.GetK()
v_te = n.Cross(in.GetK());
// k in this notation is in.GetN()
v_te = n.Cross(in.GetN());
v_te = v_te.Unit();
v_tm = -v_te.Cross(in.GetK());
v_tm = -v_te.Cross(in.GetN());
v_tm = v_tm.Unit();
if(dbg) {
printf(" v_te = "); printv(v_te);
415,13 → 430,12
}
 
double cosAf = v_te * in.GetP();
double alpha = acos(cosAf);
if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE);
if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, TMath::ACos(cosAf)*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
433,7 → 447,7
// --------------- refraction from n1 to n2 -----------------------------------
// ----------------------------------------------------------------------------
case SURF_REFRA:
cosTi = in.GetK() * n;
cosTi = in.GetN() * 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);
453,7 → 467,7
// If n1>n2 and theta>thetaCritical, total reflection
if(cosTi < cosTtotal) {
if(dbg) printf(" TOTAL\n");
transmit = in.GetK() + sign_n*2*cosTi*n;
transmit = in.GetN() + sign_n*2*cosTi*n;
 
if(dbg) {
cosTN = TMath::Abs(transmit.Unit() * n);
461,24 → 475,9
}
out->Set(intersect, transmit);
 
// 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);
 
// Shift?
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
} else {
// reflection or refraction according to Fresnel equations
487,7 → 486,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.GetK() + sign_n*(N1_N2(sign_n)*cosTi - cosTt)*n;
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;
522,17 → 521,14
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.GetK() + sign_n*2*cosTi*n;
TVector3 v_tm_r = -v_te.Cross(transmit);
v_tm_r = v_tm_r.Unit();
transmit = in.GetN() + sign_n*2*cosTi*n;
out->Set(intersect, transmit);
//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);
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
}
 
545,12 → 541,12
case SURF_REFLE:
p_ref = rand.Uniform(0.0, 1.0);
if(p_ref < reflection) { // se odbije
cosTi = in.GetK() * n;
transmit = in.GetK() - 2*cosTi*n;
cosTi = in.GetN() * n;
transmit = in.GetN() - 2*cosTi*n;
out->Set(intersect, transmit);
return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
} else { // se ne odbije
transmit = in.GetK();
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
560,17 → 556,17
case SURF_IMPER:
p_ref = rand.Uniform(0.0, 1.0);
if(p_ref < reflection) { // se odbije
cosTi = in.GetK() * n;
cosTi = in.GetN() * n;
if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj
transmit = in.GetK() - 2*cosTi*n;
transmit = in.GetN() - 2*cosTi*n;
out->Set(intersect, transmit);
} else { // ni tot. odboja
transmit = in.GetK();
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
} else { // se ne odbije
transmit = in.GetK();
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
583,7 → 579,10
 
return REFRACTION;
}
//=================================================================================
 
 
//=================================================================================
Guide::Guide(TVector3 center0, DetectorParameters &parameters) :
_d(parameters.getD()),
_n1(parameters.getN1()),
641,7 → 640,7
TVector3 activePosition(center);
activePosition += TVector3(_d, 0, 0);
TVector3 normal(1,0,0);
grease = new CPlaneR(activePosition, normal, 0.95*a/2.0);
grease = new CPlaneR(activePosition, normal, a/2.0);
 
if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1);
 
751,7 → 750,7
if (!hitActive) propagation = noCoupling->PropagateRay(ray0, &ray1, &vec1);
}
// check on which side the vector is?
TVector3 ray = ray1.GetK();
TVector3 ray = ray1.GetN();
TVector3 exitNormal = s_side[5]->GetN();
if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal);
if (ray*exitNormal > 0) {
807,14 → 806,17
*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);
826,8 → 828,9
line3d->DrawClone();
}
}
//=================================================================================
 
 
//=================================================================================
int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
{
double num, den; //stevec, imenovalec
839,7 → 842,7
 
double D = - n*center;
num = n*ray.GetR() + D;
den = n*ray.GetK();
den = n*ray.GetN();
 
if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den);
 
855,7 → 858,7
if(dbg) printf("t = %.4lf | ", t);
 
tmp = ray.GetR();
tmp -= t*ray.GetK();
tmp -= t*ray.GetN();
*vec = tmp;
 
if(dbg) {printv(tmp); printf(" | Rv = %.4lf <> R = %.4lf\n", ((tmp - center).Mag()), _r);}
866,7 → 869,7
else
return 0;
}
 
//-----------------------------------------------------------------------------
void CPlaneR::Draw(int color, int width)
{
const int NN = 32;
885,8 → 888,10
}
arc->Draw();
}
//=================================================================================
 
 
//=================================================================================
CDetector::CDetector(TVector3 center0, DetectorParameters& parameters) :
center(center0),
glass_on(parameters.getGlassOn()),
1043,8 → 1048,6
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 {
1107,8 → 1110,8
if(fate == missed) {
if (dbg) printf("Detector: fate=missed\n");
TVector3 r = ray1->GetR();
TVector3 k = ray1->GetK();
ray1->Set(r,k);
TVector3 n = ray1->GetN();
ray1->Set(r,n);
ray1->DrawS(center.x(), 10.0);
} else {
for(p_i = 0; p_i < n_points-1; p_i++) {
1189,7 → 1192,7
delete rayin;
return fate;
}
 
//-----------------------------------------------------------------------------
void CDetector::Draw(int width)
{
if(guide_on) {
1207,8 → 1210,8
//window_circle->Draw(col_glass, width);
active->Draw(col_active, width);
}
//=================================================================================
 
 
Plate::Plate(DetectorParameters& parameters)
{
TVector3 center = CENTER;
1281,7 → 1284,6
fate = missed;
} else if(result == REFLECTION) {
if (dbg) printf("PLATE: reflected\n");
ray0 = ray1;
fate = backreflected;
} else {
points[0] = ray1.GetR();
1298,7 → 1300,6
}
points[n_odb] = vec1;
if(inters_i == 0) {
ray0 = ray1;
fate = backreflected;
break;} // backreflection
 
1332,5 → 1333,6
*out = ray0;
return fate;
};
//=============================================================================================================================== <<<<<<<<
 
 
/lightguide/trunk/src/raySimulator.cpp
60,7 → 60,7
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
*/
//-----------------------------------------------------------------------------
void Init(double rsys_scale = 10.0)
void Init(double rsys_scale = 7.0)
{
RTSetStyle(gStyle);
gStyle->SetOptStat(10);
70,8 → 70,7
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);
178,9 → 177,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;
192,32 → 191,11
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));
// 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);
TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
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;
 
225,8 → 203,7
}
//-----------------------------------------------------------------------------
// 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;
247,9 → 224,6
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);
273,8 → 247,7
//-----------------------------------------------------------------------------
// 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;
304,10 → 277,6
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)
314,8 → 283,8
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
delete ray1;
delete ray0;
//delete ray0;
//delete ray1;
}
 
double acceptance = 0.0;
335,8 → 304,7
// 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;
437,8 → 405,7
 
// 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;
451,6 → 418,9
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;
491,17 → 461,10
//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) {