Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 71 → Rev 72

/lightguide/trunk/include/guide.h
37,46 → 37,49
class CRay
{
public:
CRay() :
r(TVector3(0,0,0)),
n(TVector3(1,0,0)),
p(TVector3(0,1,0)),
color(1)
{};
CRay(TVector3 r0, TVector3 n0) :
r(r0),
n(n0.Unit()),
p(TVector3(0,1,0)),
color(1)
{};
CRay(double x0, double y0, double z0, double l0, double m0, double n0) :
CRay() :
r(TVector3(0,0,0)),
n(TVector3(1,0,0)),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
CRay(TVector3 r0, TVector3 n0) :
r(r0),
n(n0.Unit()),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
CRay(double x0, double y0, double z0, double l0, double m0, double n0) :
r(TVector3(x0,y0,z0)),
n(TVector3(l0,m0,n0).Unit()),
p(TVector3(0,1,0)),
color(1)
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
 
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();};
//inline CRay & operator = (const CRay &);
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();};
 
TVector3 GetR() const {return r;};
TVector3 GetN() const {return n;};
TVector3 GetP() const {return p;};
void Print();
void Draw();
void Draw(double x_from, double x_to);
void DrawS(double x_from, double t);
//inline CRay & operator = (const CRay &);
 
TVector3 GetR() const {return r;};
TVector3 GetN() const {return n;};
TVector3 GetP() const {return p;};
 
void Print();
void Draw();
void Draw(double x_from, double x_to);
void DrawS(double x_from, double t);
 
private:
TVector3 r;
TVector3 n;
TVector3 p; //r = point on line, n = normal, p = polarization
int color;
TVector3 r;
TVector3 n;
TVector3 p; //r = point on line, n = normal, p = polarization
int color;
};
//=================================================================================
 
88,40 → 91,59
// ^ |
// | |
// | |
// | �
// |
// r0<-----r3
class CPlane4
{
 
public:
CPlane4();
CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);
CPlane4(TVector3 *vr);
void Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);
void Set(TVector3 *vr) {Set(vr[0], vr[1], vr[2], vr[3]);};
void FlipN(){n = -n;};
CPlane4();
CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);
CPlane4(TVector3 *vr);
void Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);
void Set(TVector3 *vr) {Set(vr[0], vr[1], vr[2], vr[3]);};
void FlipN(){n = -n;};
 
int GetIntersection(TVector3 *vec, CRay ray);
int IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double);
int IsVectorIn(TVector3 vec);
int TestIntersection(CRay in); // ray in & ??? plane
int TestIntersection(TVector3 *vec, CRay in); // plane defined with normal vec and ray in
TVector3 GetN() {return n;};
int GetIntersection(TVector3 *vec, CRay ray);
int IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double);
int IsVectorIn(TVector3 vec);
int TestIntersection(CRay in); // ray in & ??? plane
int TestIntersection(TVector3 *vec, CRay in); // plane defined with normal vec and ray in
 
void Print();
void Draw(int color = 1, int width = 1);
private:
TVector3 r[4], n;
double A, B, C, D;
TVector3 edge[4]; // vektorji stranic
double angle_r[4]; // koti ob posameznem vogalu
TVector3 GetN() {return n;};
 
void Print();
void Draw(int color = 1, int width = 1);
 
private:
TVector3 r[4], n;
double A, B, C, D;
TVector3 edge[4]; // vektorji stranic
double angle_r[4]; // koti ob posameznem vogalu
};
//=================================================================================
// ravnina - krog
class CPlaneR
{
 
public:
CPlaneR(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
 
void Set(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
 
int TestIntersection(TVector3 *vec, CRay in);
 
void Draw(int color = 1, int NN = 32);
 
private:
TVector3 n, center;
double _r;
};
//=================================================================================
 
//=================================================================================
// ravna opticna povrsina: refractor, zrcalo ali povrsina s totalnim odbojem
#define SURF_DUMMY 0
#define SURF_REFRA 1
135,40 → 157,40
{
 
public:
CSurface(int type0 = 0);
CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4,
double n10, double n20, double reflectivity);
CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity);
CSurface(int type0 = 0);
CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4,
double n10, double n20, double reflectivity);
CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity);
 
void SetV(TVector3 *vr){CPlane4::Set(vr);};
void SetType(int type0){type = type0;};
void SetIndex(double n10, double n20);
void SetReflection(double reflectivity){reflection = reflectivity;};
void SetFresnel(int f = 1) {fresnel = f;};
void SetV(TVector3 *vr){CPlane4::Set(vr);};
void SetType(int type0){type = type0;};
void SetIndex(double n10, double n20);
void SetReflection(double reflectivity){reflection = reflectivity;};
void SetFresnel(int f = 1) {fresnel = f;};
 
void Set(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
{type = type0; CPlane4::Set(vr); SetIndex(n10, n20); reflection = reflectivity;};
void Set(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
{type = type0; CPlane4::Set(r1, r2, r3, r4); SetIndex(n10, n20); reflection = reflectivity;};
void Set(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
{type = type0; CPlane4::Set(vr); SetIndex(n10, n20); reflection = reflectivity;};
void Set(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
{type = type0; CPlane4::Set(r1, r2, r3, r4); SetIndex(n10, n20); reflection = reflectivity;};
 
int PropagateRay(CRay in, CRay *out, TVector3 *intersection);
double N1_N2(int sign) { return ((sign > 0) ? n1/n2 : n2/n1);};
int PropagateRay(CRay in, CRay *out, TVector3 *intersection);
 
void Print(){printf("Type = %d\n", type); CPlane4::Print();};
private:
int type; //0 = dummy; 1 = refractor; 2 = reflector; 3 = total reflection
double n1, n2, n1_n2; //index of refraction, n1 @ +normal, n2 @ -normal
double reflection; // odbojnost stranic
TRandom rand; // za racunanje verjetnosti odboja od zrcala
double cosTtotal; //cosinus mejnega kota totalnega odboja za dana n1 in n2
int fresnel; // ali naj uposteva Fresnelove enacbe
double N1_N2(int sign) { return ((sign > 0) ? n1/n2 : n2/n1);};
 
void Print(){printf("Type = %d\n", type); CPlane4::Print();};
 
private:
int type; //0 = dummy; 1 = refractor; 2 = reflector; 3 = total reflection
double n1, n2, n1_n2; //index of refraction, n1 @ +normal, n2 @ -normal
double reflection; // odbojnost stranic
TRandom rand; // za racunanje verjetnosti odboja od zrcala
double cosTtotal; //cosinus mejnega kota totalnega odboja za dana n1 in n2
 
int fresnel; // ali naj uposteva Fresnelove enacbe
};
//=================================================================================
 
 
//=================================================================================
// +------- n1
// | -------
190,50 → 212,52
class DetectorParameters
{
public:
DetectorParameters(double a, double b, double d, double active, double n1, double n2, double n3, TVector3 gap):
_a(a),
_b(b),
_d(d),
_active(active),
_n1(n1),
_n2(n2),
_n3(n3),
_gap(gap),
_fresnel(1),
_guideOn(1),
offsetY(b/2.0),
offsetZ(b/2.0),
DetectorParameters(double a, double b, double d, double active, double n1, double n2, double n3, TVector3 gap, bool coupling):
_a(a),
_b(b),
_d(d),
_active(active),
_n1(n1),
_n2(n2),
_n3(n3),
_gap(gap),
_fresnel(1),
_guideOn(1),
offsetY(b/2.0),
offsetZ(b/2.0),
_plateOn(1),
_plateWidth(1),
_glassOn(0),
_glassD(0)
_plateWidth(1),
_glassOn(0),
_glassD(0),
_coupling(coupling)
{};
DetectorParameters(double a, double b, double d):
_a(a),
_b(b),
_d(d),
_active(a),
_n1(1.0),
_n2(1.53),
_n3(1.46),
_gap(TVector3(0,0,0)),
_fresnel(1),
_guideOn(1),
offsetY(b/2.0),
offsetZ(b/2.0),
_plateOn(1),
_plateWidth(1),
_glassOn(0),
_glassD(0)
_a(a),
_b(b),
_d(d),
_active(a),
_n1(1.0),
_n2(1.53),
_n3(1.46),
_gap(TVector3(0,0,0)),
_fresnel(1),
_guideOn(1),
offsetY(b/2.0),
offsetZ(b/2.0),
_plateOn(1),
_plateWidth(1),
_glassOn(0),
_glassD(0),
_coupling(false)
{};
~DetectorParameters() {};
 
void setGuide(double a, double b, double d) {
_a = a;
_b = b;
//_M = b/a;
_d = d;
};
};
void setGap(double x, double y, double z) { _gap = TVector3(x,y,z); };
void setFresnel(int fresnel) { _fresnel = fresnel; };
void setGlass(int glassOn, double glassD) { _glassOn = glassOn; _glassD = glassD; };
240,6 → 264,7
void setGuideOn(int guideOn) { _guideOn = guideOn; };
void setPlate(int plateOn, double plateWidth) { _plateOn = plateOn; _plateWidth = plateWidth; };
void setIndices(double n1, double n2, double n3) { _n1 = n1; _n2 = n2; _n3 = n3; };
void setCoupling(bool coupling) { _coupling = coupling; };
 
double getLightYield() { return ( Power(_b,2)/Power(_active,2)); };
double getM() { return (_b/_a); };
259,7 → 284,9
int getPlateOn() { return _plateOn; };
double getOffsetY() { return offsetY; };
double getOffsetZ() { return offsetZ; };
// bad coupling in the case the small amount of grease was applied
bool badCoupling() { return _coupling; };
 
private:
double _a;
double _b;
269,90 → 296,77
double _n2;
double _n3;
TVector3 _gap;
 
int _fresnel;
 
int _guideOn;
double offsetY;
double offsetZ;
 
int _plateOn;
double _plateWidth;
 
int _glassOn;
double _glassD;
bool _coupling;
 
};
 
class Guide
{
public:
Guide(TVector3 center0, DetectorParameters& parameters);
~Guide() {
for (int jk=0; jk<6; jk++) delete s_side[jk];
delete hfate;
delete hnodb_all;
delete hnodb_exit;
delete hin;
delete hout;
}
Fate PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points);
double getD() { return _d; }
double getN1() { return n1; };
double getN2() { return n2; };
double getN3() { return n3; };
TH1F* GetHFate() const {return hfate;};
TH1F* GetHNOdb_all() const {return hnodb_all;};
TH1F* GetHNOdb_exit() const {return hnodb_exit;};
TH2F* GetHIn() const {return hin;};
TH2F* GetHOut() const {return hout;};
int GetExitHits() {return (int)hfate->GetBinContent(5);};
int GetEnteranceHits() {return (int)hfate->GetBinContent(6);};
void GetVFate(int *out);
int GetMAXODB() const {return MAX_REFLECTIONS;};
void Draw(int color = 1, int width = 1);
void DrawSkel(int color = 1, int width = 1);
Guide(TVector3 center0, DetectorParameters& parameters);
~Guide() {
for (int jk=0; jk<6; jk++) delete s_side[jk];
delete grease;
delete noCoupling;
delete hfate;
delete hnodb_all;
delete hnodb_exit;
delete hin;
delete hout;
}
 
Fate PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points);
 
double getD() { return _d; }
double getN1() { return _n1; };
double getN2() { return _n2; };
double getN3() { return _n3; };
TH1F* GetHFate() const {return hfate;};
TH1F* GetHNOdb_all() const {return hnodb_all;};
TH1F* GetHNOdb_exit() const {return hnodb_exit;};
TH2F* GetHIn() const {return hin;};
TH2F* GetHOut() const {return hout;};
int GetExitHits() {return (int)hfate->GetBinContent(5);};
int GetEnteranceHits() {return (int)hfate->GetBinContent(6);};
void GetVFate(int *out);
int GetMAXODB() const {return MAX_REFLECTIONS;};
 
void Draw(int color = 1, int width = 1);
void DrawSkel(int color = 1, int width = 1);
 
private:
Fate fate;
double _d; // parameters
double n1, n2, n3; // refractive index: n1 above entry surface, n2 inside, n3 after exit
double _r;
double absorption;
double A;
TRandom rand; // for material absorption
CSurface *s_side[6];
TVector3 center, vodnik_edge[8];
TH1F *hfate, *hnodb_all, *hnodb_exit;
TH2F *hin, *hout;
double _d; // parameters
double _n1, _n2, _n3; // refractive index: n1 above entry surface, n2 inside, n3 after exit
double _r;
double _absorption;
double _A;
bool _badCoupling;
TRandom rand; // for material absorption
CSurface *s_side[6];
CPlaneR* grease;
CSurface* noCoupling;
TVector3 center, vodnik_edge[8];
 
TH1F *hfate, *hnodb_all, *hnodb_exit;
TH2F *hin, *hout;
};
//=================================================================================
 
//=============================================================================================================================== <<<<<<<<
// ravnina - krog
class CPlaneR
{
 
public:
CPlaneR(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
void Set(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
int TestIntersection(TVector3 *vec, CRay in);
void Draw(int color = 1, int NN = 32);
private:
TVector3 n, center;
double _r;
};
//=================================================================================
 
class Plate
{
361,14 → 375,14
~Plate() {
for (int jk=0; jk<6; jk++) delete sides[jk]; // the same, needs solution
};
 
void draw(int color, int width);
void drawSkel(int color, int width);
Fate propagateRay(CRay, CRay*, int*, TVector3*);
 
private:
TVector3 plate_edge[8];
CSurface *sides[6];
TVector3 plate_edge[8];
CSurface *sides[6];
};
 
// ================================================================================
376,28 → 390,28
class CDetector
{
public:
CDetector(TVector3 center0, DetectorParameters& parameters);
~CDetector() {
delete glass;
delete glass_circle;
delete hglass;
CDetector(TVector3 center0, DetectorParameters& parameters);
~CDetector() {
delete glass;
delete glass_circle;
delete hglass;
delete active;
delete hactive;
delete hlaser;
delete detector;
delete hdetector;
delete guide;
delete plate;
delete histoPlate;
//delete window;
//delete window_circle;
//delete hwindow;
}
delete hactive;
delete hlaser;
delete detector;
delete hdetector;
delete guide;
delete plate;
delete histoPlate;
//delete window;
//delete window_circle;
//delete hwindow;
}
 
//void SetLGType(int in = SURF_REFRA, int side = SURF_REFRA, int out = SURF_REFRA)
// {type_in = in; type_side = side; type_out = out;};
//void SetLGType(int in = SURF_REFRA, int side = SURF_REFRA, int out = SURF_REFRA)
// {type_in = in; type_side = side; type_out = out;};
 
/* void SetLG(double SiPM0, double M0, double d0, double n10, double n20, double n30, double R0)
/* void SetLG(double SiPM0, double M0, double d0, double n10, double n20, double n30, double R0)
{ SiPM=SiPM0;
M=M0;
d=d0;
406,92 → 420,93
n3=n30;
R=R0;
}; */
//void SetR(double R0) {R = R0;};
//void SetGlass(int glass_on0, double glass_d0)
//{glass_on = glass_on0; glass_d = glass_d0;};
//void SetGap(double x_gap0, double y_gap0, double z_gap0)
//{x_gap = x_gap0; y_gap = y_gap0; z_gap = z_gap0;};
void SetRCol(int in, int lg, int out, int gla)
{col_in = in; col_lg = lg; col_out = out; col_rgla = gla;};
void SetDCol(int LG0, int glass0, int active0)
{col_LG = LG0; col_glass = glass0; col_active = active0;};
//void SetFresnel(int b)
//{fresnel = b;};
//void SetAbsorption(int b, double A0)
//{absorption = b; A = A0;};
//void SetWindow(double wR, double d0) {window_R = wR; window_d = d0;};
//void SetGuideOn(int b)
//{guide_on = b;};
//void SetR(double R0) {R = R0;};
//void SetGlass(int glass_on0, double glass_d0)
//{glass_on = glass_on0; glass_d = glass_d0;};
//void SetGap(double x_gap0, double y_gap0, double z_gap0)
//{x_gap = x_gap0; y_gap = y_gap0; z_gap = z_gap0;};
void SetRCol(int in, int lg, int out, int gla)
{col_in = in; col_lg = lg; col_out = out; col_rgla = gla;};
void SetDCol(int LG0, int glass0, int active0)
{col_LG = LG0; col_glass = glass0; col_active = active0;};
//void SetFresnel(int b)
//{fresnel = b;};
//void SetAbsorption(int b, double A0)
//{absorption = b; A = A0;};
//void SetWindow(double wR, double d0) {window_R = wR; window_d = d0;};
//void SetGuideOn(int b)
//{guide_on = b;};
 
//void Init();
//void Init();
 
int Propagate(CRay, CRay*, int);
int Propagate(CRay, CRay*, int);
 
Guide* GetLG() const {return guide;};
//TH2F* GetHWindow() const {return hwindow;};
TH2F* GetHGlass() const {return hglass;};
TH2F* GetHActive() const {return hactive;};
TH2F* GetHLaser() const {return hlaser;};
TH2F* GetHDetector() const {return hdetector;};
TH2F* GetHPlate() const {return histoPlate;};
//double GetSiPM() const {return SiPM;};
//double GetM() const {return M;};
//double GetD() const {return d;};
//double GetR() const {return R;};
//TVector3 GetVGap() const {return (TVector3(x_gap, y_gap, z_gap));};
//int IsGlass() {return glass_on;};
//double getGlassWidth() {return glass_d;};
//double GetWindowR() {return window_R;};
//double GetWindowD() {return window_d;};
Guide* GetLG() const {return guide;};
//TH2F* GetHWindow() const {return hwindow;};
TH2F* GetHGlass() const {return hglass;};
TH2F* GetHActive() const {return hactive;};
TH2F* GetHLaser() const {return hlaser;};
TH2F* GetHDetector() const {return hdetector;};
TH2F* GetHPlate() const {return histoPlate;};
 
void Draw(int width = 2);
private:
Fate fate;
TVector3 center;
//double GetSiPM() const {return SiPM;};
//double GetM() const {return M;};
//double GetD() const {return d;};
//double GetR() const {return R;};
//TVector3 GetVGap() const {return (TVector3(x_gap, y_gap, z_gap));};
//int IsGlass() {return glass_on;};
//double getGlassWidth() {return glass_d;};
//double GetWindowR() {return window_R;};
//double GetWindowD() {return window_d;};
 
//int type_in, type_side, type_out;
//double SiPM, M, d, n1, n2, n3, R;
//double detectorActive;
void Draw(int width = 2);
 
int glass_on;
double glass_d;
CSurface *glass;
CPlaneR *glass_circle;
TH2F *hglass;
private:
Fate fate;
TVector3 center;
 
//double x_gap, y_gap, z_gap;
CPlane4 *active;
TH2F *hactive, *hlaser;
//int type_in, type_side, type_out;
//double SiPM, M, d, n1, n2, n3, R;
//double detectorActive;
 
CPlane4 *detector;
TH2F *hdetector;
int col_in, col_lg, col_out, col_rgla;
int col_LG, col_glass, col_active;
//int fresnel, absorption;
//double A;
int guide_on;
//double window_R, window_d;
//CSurface *window;
//CPlaneR *window_circle;
//TH2F *hwindow;
Guide *guide;
Plate *plate;
double _plateWidth;
int _plateOn;
TH2F *histoPlate;
double offsetY;
double offsetZ;
 
int glass_on;
double glass_d;
CSurface *glass;
CPlaneR *glass_circle;
TH2F *hglass;
 
//double x_gap, y_gap, z_gap;
CPlane4* active;
CPlaneR* grease;
TH2F *hactive, *hlaser;
 
CPlane4 *detector;
TH2F *hdetector;
 
int col_in, col_lg, col_out, col_rgla;
int col_LG, col_glass, col_active;
 
//int fresnel, absorption;
//double A;
 
int guide_on;
 
//double window_R, window_d;
//CSurface *window;
//CPlaneR *window_circle;
//TH2F *hwindow;
 
Guide *guide;
Plate *plate;
 
double _plateWidth;
int _plateOn;
 
TH2F *histoPlate;
 
double offsetY;
double offsetZ;
};
 
 
/lightguide/trunk/src/userFunctions.cpp
12,8 → 12,8
void showVisual(int b) {show_3d = b;}
void showData(int b) {show_data = b;}
 
// Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap)
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0));
// Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap, badCoupling)
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0), false);
// Print the detector parameters
void getParameters();
 
193,9 → 193,10
}
//------------------------------------------------------------------------------------------
// Propagate NN rays generated as grid and show the statistics
void LGG(int NN = 10, double theta = 0.0)
void LGG(int NN=10, double theta=0.0, bool coupling=false)
{
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
parameters.setCoupling(coupling);
CDetector *detector = new CDetector(CENTER, parameters);
Init();
double izkoristek = Grid(detector, parameters, NN, theta);
//printf("izkoristek = %.3lf\n", izkoristek);
207,7 → 208,7
// Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
{
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
CDetector *detector = new CDetector(CENTER, parameters);
Init();
double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
//printf("izkoristek = %.3lf\n", izkoristek);
231,19 → 232,22
//detector->Draw();
}
 
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, int nnrays = 30, int showr = 0)
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0)
{
Init();
parameters.setCoupling(coupling);
CDetector *detector = new CDetector(CENTER, parameters);
//CDetector detector = new CDetector();
double izkoristek = beamtest(detector, parameters, NN, 18.5, phiMin, phiMax, nnrays, showr);
//printf("izkoristek = %.3lf\n", izkoristek);
PrintGuideHead();
PrintGuideStat(izkoristek);
DrawData(detector, parameters, 18.5, izkoristek);
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
//canvasDetector->cd();
//detector->Draw();
const double theta = 18.5;
double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr);
PrintGuideHead();
PrintGuideStat(izkoristek);
DrawData(detector, parameters, 18.5, izkoristek);
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
//canvasDetector->cd();
//detector->Draw();
}
/lightguide/trunk/src/guide.cpp
5,20 → 5,20
// vector output shortcut
void printv(TVector3 v)
{
printf("(x,y,z) = (%.4lf, %.4lf, %.4lf)\n", v.x(), v.y(), v.z());
printf("(x,y,z) = (%.4lf, %.4lf, %.4lf)\n", v.x(), v.y(), v.z());
}
// TVector3::Rotate does not seem accurate enough
TVector3 rotatey(TVector3 v, double theta)
{
return TVector3(v.x() * TMath::Cos(theta) + v.z() * TMath::Sin(theta),
v.y(),
-v.x() * TMath::Sin(theta) + v.z() * TMath::Cos(theta));
return TVector3(v.x() * TMath::Cos(theta) + v.z() * TMath::Sin(theta),
v.y(),
-v.x() * TMath::Sin(theta) + v.z() * TMath::Cos(theta));
}
// another shortcut not found in TMath
int sign(double in)
{
if(in >= 0.0) return 1;
else return -1;
if(in >= 0.0) return 1;
else return -1;
}
//=================================================================================
 
25,13 → 25,13
//-----------------------------------------------------------------------------
void CRay::Set(TVector3 r0, TVector3 n0)
{
r = r0; n = n0.Unit();
r = r0; n = n0.Unit();
}
//-----------------------------------------------------------------------------
//void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0)
//{
//r.SetXYZ(x0, y0, z0);
//n.SetXYZ(l0, m0, n0); n = n.Unit();
//r.SetXYZ(x0, y0, z0);
//n.SetXYZ(l0, m0, n0); n = n.Unit();
//}
//-----------------------------------------------------------------------------
/*
45,60 → 45,60
//-----------------------------------------------------------------------------
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());
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());
}
//-----------------------------------------------------------------------------
void CRay::Draw()
{
double t = 50.0;
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->SetLineWidth(1);
line3d->SetLineColor(color);
double t = 50.0;
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->SetLineWidth(1);
line3d->SetLineColor(color);
 
line3d->Draw();
line3d->Draw();
}
//-----------------------------------------------------------------------------
void CRay::Draw(double x_from, double x_to)
{
double A1, A2;
double A1, A2;
TPolyLine3D *line3d = new TPolyLine3D(2);
 
if(n.x() < MARGIN) {
A1 = A2 = 0.0;
} else {
A1 = (x_from - r.x())/n.x();
A2 = (x_to - r.x())/n.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->SetLineWidth(1);
line3d->SetLineColor(color);
if(n.x() < MARGIN) {
A1 = A2 = 0.0;
} else {
A1 = (x_from - r.x())/n.x();
A2 = (x_to - r.x())/n.x();
}
 
line3d->Draw();
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);
double A1;
TPolyLine3D *line3d = new TPolyLine3D(2);
 
if(n.x() < MARGIN)
A1 = 0.0;
else
A1 = (x_from - r.x())/n.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->SetLineWidth(1);
line3d->SetLineColor(color);
if(n.x() < MARGIN)
A1 = 0.0;
else
A1 = (x_from - r.x())/n.x();
 
line3d->Draw();
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();
}
//=================================================================================
 
105,96 → 105,96
 
//=================================================================================
CPlane4::CPlane4() :
n(TVector3(1.0, 0.0, 0.0)),
A(0),
B(0),
C(0),
D(0)
n(TVector3(1.0, 0.0, 0.0)),
A(0),
B(0),
C(0),
D(0)
{ r[0] = TVector3(0.0,-1.0,-1.0);
r[1] = TVector3(0.0,-1.0, 1.0);
r[2] = TVector3(0.0, 1.0, 1.0);
r[3] = TVector3(0.0, 1.0,-1.0);
for(int i=0;i<4;i++) edge[i] = TVector3(0,0,0);
for(int i=0;i<4;i++) angle_r[i] = 0;
r[1] = TVector3(0.0,-1.0, 1.0);
r[2] = TVector3(0.0, 1.0, 1.0);
r[3] = TVector3(0.0, 1.0,-1.0);
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);
//}
//-----------------------------------------------------------------------------
// za izracun parametrov ravnine je en vektor prevec, vendar tega
// rabim kot zadnji vogal poligona
//void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
//{
double x1,y1,z1, x2,y2,z2, x3,y3,z3;
x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
//Set(r1, r2, r3, r4);
//}
//-----------------------------------------------------------------------------
// za izracun parametrov ravnine je en vektor prevec, vendar tega
// rabim kot zadnji vogal poligona
//void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
//{
double x1,y1,z1, x2,y2,z2, x3,y3,z3;
 
A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
 
r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
n.SetXYZ(A, B, C);
n = n.Unit();
A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
 
for(int i=0;i<4;i++)
edge[i] = r[i-3 ? i+1 : 0] - r[i];
r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
n.SetXYZ(A, B, C);
n = n.Unit();
 
for(int i=0;i<4;i++)
angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
for(int i=0;i<4;i++)
edge[i] = r[i-3 ? i+1 : 0] - r[i];
 
for(int i=0;i<4;i++)
angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
};
 
void CPlane4::Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4)
{
double x1,y1,z1, x2,y2,z2, x3,y3,z3;
x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
 
A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
x1 = r1.x(); y1 = r1.y(); z1 = r1.z();
x2 = r2.x(); y2 = r2.y(); z2 = r2.z();
x3 = r3.x(); y3 = r3.y(); z3 = r3.z();
 
r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
n.SetXYZ(A, B, C);
n = n.Unit();
A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
 
for(int i=0;i<4;i++)
edge[i] = r[i-3 ? i+1 : 0] - r[i];
r[0] = r1; r[1] = r2; r[2] = r3; r[3] = r4;
n.SetXYZ(A, B, C);
n = n.Unit();
 
for(int i=0;i<4;i++)
angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
for(int i=0;i<4;i++)
edge[i] = r[i-3 ? i+1 : 0] - r[i];
 
for(int i=0;i<4;i++)
angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
};
 
CPlane4::CPlane4(TVector3 *vr)
{
double x1,y1,z1, x2,y2,z2, x3,y3,z3;
x1 = vr[0].x(); y1 = vr[0].y(); z1 = vr[0].z();
x2 = vr[1].x(); y2 = vr[1].y(); z2 = vr[1].z();
x3 = vr[2].x(); y3 = vr[2].y(); z3 = vr[2].z();
 
A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
x1 = vr[0].x(); y1 = vr[0].y(); z1 = vr[0].z();
x2 = vr[1].x(); y2 = vr[1].y(); z2 = vr[1].z();
x3 = vr[2].x(); y3 = vr[2].y(); z3 = vr[2].z();
 
r[0] = vr[0]; r[1] = vr[1]; r[2] = vr[2]; r[3] = vr[3];
n.SetXYZ(A, B, C);
n = n.Unit();
A = y3*(z1 - z2) + y1*(z2 - z3) + y2*(z3 - z1);
B = x3*(z2 - z1) + x1*(z3 - z2) + x2*(z1 - z3);
C = x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1);
D = y3*(x1*z2 - x2*z1) + x3*(y2*z1 - y1*z2) + z3*(x2*y1 - x1*y2);
 
for(int i=0;i<4;i++)
edge[i] = r[i-3 ? i+1 : 0] - r[i];
r[0] = vr[0]; r[1] = vr[1]; r[2] = vr[2]; r[3] = vr[3];
n.SetXYZ(A, B, C);
n = n.Unit();
 
for(int i=0;i<4;i++)
angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
for(int i=0;i<4;i++)
edge[i] = r[i-3 ? i+1 : 0] - r[i];
 
for(int i=0;i<4;i++)
angle_r[i] = TMath::ACos(/*TMath::Abs*/( ((-edge[i ? i-1 : 3]).Unit()) * (edge[i].Unit()) ));
};
//-----------------------------------------------------------------------------
// posce presecisce !neskoncne! ravnine s premico (class CRay)
201,37 → 201,37
// ce najde presecisce vrne 1
int CPlane4::GetIntersection(TVector3 *vec, CRay ray)
{
TVector3 N; //nenormirani vektor (A,B,C)
double num, den; //stevec, imenovalec
double t;
TVector3 tmp;
TVector3 N; //nenormirani vektor (A,B,C)
double num, den; //stevec, imenovalec
double t;
TVector3 tmp;
 
N.SetXYZ(A,B,C);
num = N*ray.GetR() + D;
den = N*ray.GetN();
N.SetXYZ(A,B,C);
 
if (dbg) printf("t = %6.3lf / %6.3lf = %6.3lf\n", num, den, num/den);
num = N*ray.GetR() + D;
den = N*ray.GetN();
 
//if(den == 0)
if(TMath::Abs(den) < MARGIN) {
//if(num == 0)
if(TMath::Abs(num) < MARGIN) {
if (dbg) printf("The ray is on the surface!\n");
return 0; //return 2; // premica lezi na ravnini
}
else {
if (dbg) printf("The ray is parallel to the surface!\n");
return 0; // ni presecisca
}
}
t = num / den;
tmp = ray.GetR();
tmp -= t*ray.GetN();
*vec = tmp;
return 1;
if (dbg) printf("t = %6.3lf / %6.3lf = %6.3lf\n", num, den, num/den);
 
//if(den == 0)
if(TMath::Abs(den) < MARGIN) {
//if(num == 0)
if(TMath::Abs(num) < MARGIN) {
if (dbg) printf("The ray is on the surface!\n");
return 0; //return 2; // premica lezi na ravnini
}
else {
if (dbg) printf("The ray is parallel to the surface!\n");
return 0; // ni presecisca
}
}
 
t = num / den;
 
tmp = ray.GetR();
tmp -= t*ray.GetN();
*vec = tmp;
return 1;
}
//-----------------------------------------------------------------------------
// ali je vektor vec, ki lezi na ravnini skupaj z e1 in e2, med njima
239,20 → 239,20
int CPlane4::IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double angle)
{
double angle_ve1, angle_ve2;
if(dbg) printf("--- CPlane4::IsInTri ---\n");
 
angle_ve1 = TMath::ACos(/*TMath::Abs*/( (e1.Unit()) * (vec.Unit()) ));
angle_ve2 = TMath::ACos(/*TMath::Abs*/( (e2.Unit()) * (vec.Unit()) ));
if(dbg) printf("--- CPlane4::IsInTri ---\n");
 
if(dbg)
{
printf("angle_ve1 = %lf\n", angle_ve1*DEGREE);
printf("angle_ve2 = %lf\n", angle_ve2*DEGREE);
printf("angle_sum = %lf\n", (angle_ve1 + angle_ve2)*DEGREE);
printf(" angle_r = %lf\n", angle*DEGREE);
}
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_ve2 = %lf\n", 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);
return (int) !difference;
264,61 → 264,61
{
int status;
 
if(dbg) printf("--- CPlane4::IsVectorIn ---\n");
for(int i=0;i<3;i++)
{
status = IsInTri(vec - r[i], edge[i], -edge[i ? i-1 : 3], angle_r[i]);
if(dbg) printf(" [%d] vec is %s\n", i, status ? "inside" : "outside");
if(!status) return 0;
}
if(dbg) printf("--- CPlane4::IsVectorIn ---\n");
 
return 1;
for(int i=0;i<3;i++)
{
status = IsInTri(vec - r[i], edge[i], -edge[i ? i-1 : 3], angle_r[i]);
if(dbg) printf(" [%d] vec is %s\n", i, status ? "inside" : "outside");
if(!status) return 0;
}
 
return 1;
}
//-----------------------------------------------------------------------------
int CPlane4::TestIntersection(CRay in)
{
TVector3 tmp;
TVector3 tmp;
 
if( GetIntersection(&tmp, in) )
if( IsVectorIn(tmp) )
return 1;
return 0;
if( GetIntersection(&tmp, in) )
if( IsVectorIn(tmp) )
return 1;
 
return 0;
}
//-----------------------------------------------------------------------------
int CPlane4::TestIntersection(TVector3 *vec, CRay in)
{
TVector3 tmp;
TVector3 tmp;
 
if( GetIntersection(&tmp, in) )
if( IsVectorIn(tmp) ) {
*vec = tmp;
return 1;
}
return 0;
if( GetIntersection(&tmp, in) )
if( IsVectorIn(tmp) ) {
*vec = tmp;
return 1;
}
 
return 0;
}
//-----------------------------------------------------------------------------
void CPlane4::Print()
{
printf("--- CPlane4::Print() ---\n");
printf(" r=(%.2lf, %.2lf, %.2lf); n=(%.2lf, %.2lf, %.2lf); ",
r[0].x(), r[0].y(), r[0].z(), n.x(), n.y(), n.z());
printf( "(A,B,C,D)=(%.2lf, %.2lf, %.2lf, %.2lf) \n", A, B, C, D);
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);
printf("--- CPlane4::Print() ---\n");
printf(" r=(%.2lf, %.2lf, %.2lf); n=(%.2lf, %.2lf, %.2lf); ",
r[0].x(), r[0].y(), r[0].z(), n.x(), n.y(), n.z());
printf( "(A,B,C,D)=(%.2lf, %.2lf, %.2lf, %.2lf) \n", A, B, C, D);
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);
TPolyLine3D *line3d = new TPolyLine3D(5);
 
for(int i=0;i<4;i++) line3d->SetPoint(i, r[i].x(), r[i].y(), r[i].z());
line3d->SetPoint(4, r[0].x(), r[0].y(), r[0].z());
line3d->SetLineWidth(width); line3d->SetLineColor(color);
for(int i=0;i<4;i++) line3d->SetPoint(i, r[i].x(), r[i].y(), r[i].z());
line3d->SetPoint(4, r[0].x(), r[0].y(), r[0].z());
line3d->SetLineWidth(width); line3d->SetLineColor(color);
 
line3d->Draw();
line3d->Draw();
}
//=================================================================================
 
325,58 → 325,58
 
//=================================================================================
CSurface::CSurface(int type0):
type(type0)
type(type0)
{
TVector3 vr[4];
TDatime now;
TVector3 vr[4];
TDatime now;
 
vr[0].SetXYZ(0.0,-1.0,-1.0);
vr[1].SetXYZ(0.0,-1.0, 1.0);
vr[2].SetXYZ(0.0, 1.0, 1.0);
vr[3].SetXYZ(0.0, 1.0,-1.0);
//CPlane4::Set(vr);
SetIndex(1.0, 1.5);
reflection = c_reflectivity;
rand.SetSeed(now.Get());
SetFresnel();
vr[0].SetXYZ(0.0,-1.0,-1.0);
vr[1].SetXYZ(0.0,-1.0, 1.0);
vr[2].SetXYZ(0.0, 1.0, 1.0);
vr[3].SetXYZ(0.0, 1.0,-1.0);
//CPlane4::Set(vr);
SetIndex(1.0, 1.5);
 
reflection = c_reflectivity;
rand.SetSeed(now.Get());
 
SetFresnel();
}
//-----------------------------------------------------------------------------
CSurface::CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
{
TDatime now;
type = type0; CPlane4::Set(r1, r2, r3, r4);
SetIndex(n10, n20);
reflection = reflectivity;
rand.SetSeed(now.Get());
SetFresnel();
TDatime now;
 
type = type0; CPlane4::Set(r1, r2, r3, r4);
SetIndex(n10, n20);
 
reflection = reflectivity;
rand.SetSeed(now.Get());
 
SetFresnel();
}
//-----------------------------------------------------------------------------
CSurface::CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
{
TDatime now;
type = type0; CPlane4::Set(vr);
SetIndex(n10, n20);
reflection = reflectivity;
rand.SetSeed(now.Get());
SetFresnel();
TDatime now;
 
type = type0; CPlane4::Set(vr);
SetIndex(n10, n20);
 
reflection = reflectivity;
rand.SetSeed(now.Get());
 
SetFresnel();
}
//-----------------------------------------------------------------------------
void CSurface::SetIndex(double n10, double n20)
{
n1 = n10; n2 = n20; n1_n2 = n1/n2;
if(n1 > n2)
cosTtotal = TMath::Sqrt( 1 - TMath::Power(n2/n1, 2) );
else
cosTtotal = 0.0;
n1 = n10; n2 = n20; n1_n2 = n1/n2;
 
if(n1 > n2)
cosTtotal = TMath::Sqrt( 1 - TMath::Power(n2/n1, 2) );
else
cosTtotal = 0.0;
}
//-----------------------------------------------------------------------------
// sprejme zarek, vrne uklonjen/odbit zarek in presecisce
389,274 → 389,285
double cosTt; // transmited ray angle
TVector3 intersect, transmit;
 
if( !(GetIntersection(&intersect, in) == 1) )
return 0;
*intersection = intersect;
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; // s reflection coefficient
double R_tm=0; // p reflection coefficient
double R_f = 0.0;
double a_te = 0.0; // s-wave amplitude, cos Alpha
double a_tm = 0.0; // p-wave amplitude, sin Alpha
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
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_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);
}
double cosAf = v_te * in.GetP();
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
double p_ref = rand.Uniform(0.0, 1.0);
if(type == SURF_TOTAL) type = SURF_REFRA;
switch(type){
// ----------------------------------------------------------------------------
// --------------- refraction from n1 to n2 -----------------------------------
// ----------------------------------------------------------------------------
case SURF_REFRA:
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);
cosTi = TMath::Abs(cosTi);
// Check if there can be total reflection: n1 > n2
if(N1_N2(-sign_n) < 1.0)
cosTtotal = TMath::Sqrt( 1 - TMath::Power(N1_N2(-sign_n), 2) );
else
cosTtotal = 0.0;
if(dbg) printf(" cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
// reflection dependance on polarization missing
// reflection hardcoded to 0.96
if (dbg) printf(" reflection probability = %f\n", p_ref);
// If n1>n2 and theta>thetaCritical, total reflection
if(cosTi < cosTtotal) {
if(dbg) printf(" TOTAL\n");
transmit = in.GetN() + sign_n*2*cosTi*n;
if(dbg) {
cosTN = TMath::Abs(transmit.Unit() * n);
printf(" cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE);
}
out->Set(intersect, transmit);
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
} 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);
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
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);
}
if( !(GetIntersection(&intersect, in) == 1) )
return 0;
 
*intersection = intersect;
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; // s reflection coefficient
double R_tm=0; // p reflection coefficient
double R_f = 0.0;
double a_te = 0.0; // s-wave amplitude, cos Alpha
double a_tm = 0.0; // p-wave amplitude, sin Alpha
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
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_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);
}
 
double cosAf = v_te * in.GetP();
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
double p_ref = rand.Uniform(0.0, 1.0);
 
if(type == SURF_TOTAL) type = SURF_REFRA;
switch(type){
// ----------------------------------------------------------------------------
// --------------- refraction from n1 to n2 -----------------------------------
// ----------------------------------------------------------------------------
case SURF_REFRA:
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);
cosTi = TMath::Abs(cosTi);
 
// Check if there can be total reflection: n1 > n2
if(N1_N2(-sign_n) < 1.0)
cosTtotal = TMath::Sqrt( 1 - TMath::Power(N1_N2(-sign_n), 2) );
else
cosTtotal = 0.0;
 
if(dbg) printf(" cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
// reflection dependance on polarization missing
// reflection hardcoded to 0.96
if (dbg) printf(" reflection probability = %f\n", p_ref);
 
// If n1>n2 and theta>thetaCritical, total reflection
if(cosTi < cosTtotal) {
if(dbg) printf(" TOTAL\n");
transmit = in.GetN() + sign_n*2*cosTi*n;
 
if(dbg) {
cosTN = TMath::Abs(transmit.Unit() * n);
printf(" cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE);
}
out->Set(intersect, transmit);
 
// Shift?
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
} 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);
 
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
 
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);
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;
out->Set(intersect, transmit);
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
}
//}
break;
// ----------------------------------------------------------------------------
// --------------- reflection at "reflection" probability ---------------------
// ----------------------------------------------------------------------------
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;
out->Set(intersect, transmit);
return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
} else { // se ne odbije
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
break;
// total reflection from n1 to n2 with R probbability
case SURF_IMPER:
p_ref = rand.Uniform(0.0, 1.0);
if(p_ref < reflection) { // se odbije
cosTi = in.GetN() * n;
if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj
transmit = in.GetN() - 2*cosTi*n;
out->Set(intersect, transmit);
} else { // ni tot. odboja
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
} else { // se ne odbije
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
break;
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;
 
default:
*out = in;
break;
}
if (dbg) printf(" R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
}
 
return REFRACTION;
if(p_ref >= R_f) { // se lomi
if (dbg) printf(" SURFACE REFRACTED. Return.\n");
out->Set(intersect, transmit);
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;
out->Set(intersect, transmit);
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
}
 
//}
break;
 
// ----------------------------------------------------------------------------
// --------------- reflection at "reflection" probability ---------------------
// ----------------------------------------------------------------------------
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;
out->Set(intersect, transmit);
return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
} else { // se ne odbije
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
break;
 
// total reflection from n1 to n2 with R probbability
case SURF_IMPER:
p_ref = rand.Uniform(0.0, 1.0);
if(p_ref < reflection) { // se odbije
cosTi = in.GetN() * n;
if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj
transmit = in.GetN() - 2*cosTi*n;
out->Set(intersect, transmit);
} else { // ni tot. odboja
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
} else { // se ne odbije
transmit = in.GetN();
out->Set(intersect, transmit);
return ABSORBED;
}
break;
 
default:
*out = in;
break;
}
 
return REFRACTION;
}
//=================================================================================
 
 
//=================================================================================
Guide::Guide(TVector3 center0, DetectorParameters &parameters)
Guide::Guide(TVector3 center0, DetectorParameters &parameters) :
_d(parameters.getD()),
_n1(parameters.getN1()),
_n2(parameters.getN2()),
_n3(parameters.getN3()),
_r(c_reflectivity),
_absorption(0),
_A(0),
_badCoupling(parameters.badCoupling())
{
double t;
double t;
TDatime now;
rand.SetSeed(now.Get());
center = center0;
double b = parameters.getB();
double a = parameters.getA();
// if PlateOn, then n0 = n3 (optical grease), else = n1 (air)
//double n0 = (parameters.getPlateOn() ? parameters.getN3(): n1);
double n0 = (parameters.getPlateOn() ? _n2 : _n1);
int fresnel = parameters.getFresnel();
 
TDatime now;
rand.SetSeed(now.Get());
// 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[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[6].SetXYZ(_d,-t, t);
vodnik_edge[7].SetXYZ(_d,-t,-t);
 
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() ? n2 : n1);
n3 = parameters.getN3();
_r = c_reflectivity;
int fresnel = parameters.getFresnel();
for(int i = 0; i<8; i++) vodnik_edge[i] += center;
 
// 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[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[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); // 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");
(hfate->GetXaxis())->SetBinLabel(4, "LG Miss");
(hfate->GetXaxis())->SetBinLabel(5, "Exit");
(hfate->GetXaxis())->SetBinLabel(6, "Enter");
(hfate->GetXaxis())->SetBinLabel(7, "Rays");
(hfate->GetXaxis())->SetBinLabel(8, "Absorb");
hnodb_all = (TH1F*)gROOT->FindObject("hnodb_all"); if(hnodb_all) delete hnodb_all;
hnodb_all = new TH1F("hnodb_all", "N reflected", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
hnodb_exit = (TH1F*)gROOT->FindObject("hnodb_exit"); if(hnodb_exit) delete hnodb_exit;
hnodb_exit = new TH1F("hnodb_exit", "N reflected and exit", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
int nBins = nch + 1;
hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin;
hin = new TH2F("hin", "Guide entrance window", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0);
hout = (TH2F*)gROOT->FindObject("hout"); if(hout) delete hout;
hout = new TH2F("hout", "Guide exit window", nBins, -a/2.0, +a/2.0, nBins, -a/2.0, +a/2.0);
absorption = 0;
A = 0;
// 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);
// n3 - ref ind at the exit, grease, air
s_side[5] = new CSurface(SURF_REFRA, &vodnik_edge[4], _n2, _n3, _r);
s_side[5]->FlipN();
// exit surface in the case of bad coupling
noCoupling = new CSurface(SURF_REFRA, &vodnik_edge[4], _n2, 1.0, _r);
noCoupling->FlipN();
// grease = specific pattern area of coupling
TVector3 activePosition(center);
activePosition += TVector3(_d, 0, 0);
TVector3 normal(1,0,0);
grease = new CPlaneR(activePosition, normal, a/2.0);
 
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");
(hfate->GetXaxis())->SetBinLabel(4, "LG Miss");
(hfate->GetXaxis())->SetBinLabel(5, "Exit");
(hfate->GetXaxis())->SetBinLabel(6, "Enter");
(hfate->GetXaxis())->SetBinLabel(7, "Rays");
(hfate->GetXaxis())->SetBinLabel(8, "Absorb");
 
hnodb_all = (TH1F*)gROOT->FindObject("hnodb_all"); if(hnodb_all) delete hnodb_all;
hnodb_all = new TH1F("hnodb_all", "", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
 
hnodb_exit = (TH1F*)gROOT->FindObject("hnodb_exit"); if(hnodb_exit) delete hnodb_exit;
hnodb_exit = new TH1F("hnodb_exit", "", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
 
int nBins = nch + 1;
hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin;
hin = new TH2F("hin", ";x [mm]; y[mm]", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0);
 
hout = (TH2F*)gROOT->FindObject("hout"); if(hout) delete hout;
hout = new TH2F("hout", ";x [mm];y [mm]", nBins, -a/2.0, +a/2.0, nBins, -a/2.0, +a/2.0);
}
//-----------------------------------------------------------------------------
// Sledi zarku skozi vodnik. Vrne:
670,143 → 681,152
Fate Guide::PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
{
if (dbg) printf("--- GUIDE::PropagateRay ---\n");
// ray0 - incident ray
// ray1 - trans/refl ray
CRay ray0;
CRay ray1;
TVector3 vec0, vec1;
int inters_i = 0;
ray0 = in;
int n_odb = 0;
int last_hit = 0;
int propagation = 0;
int result = s_side[0]->PropagateRay(ray0, &ray1, &vec1);
if( !(result) ) {
// ce -NI- presecisca z vstopno
if (dbg) printf(" GUIDE: missed the light guide\n");
fate = missed;
//hfate->Fill(0);
} else if(result == REFLECTION) {
if (dbg) printf(" REFLECTED on the entry surface!\n");
fate = backreflected;
//hfate->Fill(-3);
} else {
if (dbg) printf(" GUIDE: ray entered\n");
points[0] = ray1.GetR();
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);
ray0 = ray1;
vec0 = vec1;
propagation = 11;
for(inters_i=0; inters_i<6; inters_i++) {
if (dbg) printf(" GUIDE: Test intersection with surface %d \n", inters_i);
if( inters_i != last_hit) {
int testBoundary = s_side[inters_i]->TestIntersection(&vec1, ray1);
if( testBoundary ) {
if (dbg) printf(" GUIDE: ray intersects with LG surface %d\n",inters_i);
break;
}
}
}
points[n_odb] = vec1;
if(inters_i == 0) {
fate = backreflected;
//hfate->Fill(backreflected);
break;
} // backreflection
// the passage is possible, test propagation
propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
if (dbg) printf(" GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
if(propagation == REFRACTION) {
fate = refracted;
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
} // no total reflection when should be
if(propagation == ABSORBED) {
fate = noreflection;
break;
} //refraction due to finite reflectivity
if(inters_i == 5) { // successfull exit
// check on which side the vector is?
TVector3 ray = ray1.GetN();
TVector3 exitNormal = s_side[5]->GetN();
//printf("theta(ray) = %lf, theta(normal5) = %lf ", ray.Theta()*DEGREE, exitNormal.Theta()*DEGREE);
//printf("phi(ray) = %lf, phi(normal5) = %lf\n", ray.Phi()*DEGREE, exitNormal.Phi()*DEGREE);
if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal);
if (ray*exitNormal > 0) {
if (dbg) printf(" GUIDE: ray is backreflected from exit window.\n");
fate = backreflected;
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
}
fate = hitExit;
hout->Fill(vec1.y(), vec1.z());
hnodb_exit->Fill(n_odb-1);
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
}
last_hit = inters_i;
}
}
//--- material absorption ---
if(absorption) {
double travel = 0.0;
printf("n_odb = %d\n", n_odb); //dbg
for(int point = 0; point < n_odb-1; point++) {
travel += (points[point] - points[point+1]).Mag();
printf("travel = %lf\n", travel); //dbg
}
double T_abs = TMath::Exp(-travel/A);
printf("T_abs = %lf\n", T_abs); //dbg
double p_abs = rand.Uniform(0.0, 1.0);
printf("p_abs = %lf\n", p_abs); //dbg
if(p_abs > T_abs) fate = absorbed; // absorption
}
//--- material absorption ---
hfate->Fill(fate);
hfate->Fill(rays);
hnodb_all->Fill(n_odb-2);
*n_points = n_odb+1;
*out = ray0;
return fate;
 
ray0 = in;
int n_odb = 0;
int last_hit = 0;
int propagation = 0;
int result = s_side[0]->PropagateRay(ray0, &ray1, &vec1);
if( !(result) ) {
// ce -NI- presecisca z vstopno
if (dbg) printf(" GUIDE: missed the light guide\n");
fate = missed;
//hfate->Fill(0);
} else if(result == REFLECTION) {
if (dbg) printf(" REFLECTED on the entry surface!\n");
fate = backreflected;
//hfate->Fill(-3);
} else {
if (dbg) printf(" GUIDE: ray entered\n");
points[0] = ray1.GetR();
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);
ray0 = ray1;
vec0 = vec1;
propagation = 11;
for(inters_i=0; inters_i<6; inters_i++) {
if (dbg) printf(" GUIDE: Test intersection with surface %d \n", inters_i);
if( inters_i != last_hit) {
int testBoundary = s_side[inters_i]->TestIntersection(&vec1, ray1);
if( testBoundary ) {
if (dbg) printf(" GUIDE: ray intersects with LG surface %d\n",inters_i);
break;
}
}
}
points[n_odb] = vec1;
if(inters_i == 0) {
fate = backreflected;
//hfate->Fill(backreflected);
break;
} // backreflection
 
// the passage is possible, test propagation
propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
 
if (dbg) printf(" GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
 
 
if(propagation == ABSORBED) {
fate = noreflection;
break;
} //refraction due to finite reflectivity
 
if(inters_i == 5) {
if (_badCoupling) {
TVector3 hitVector(0,0,0);
bool hitActive = grease->TestIntersection(&hitVector, ray0);
if (hitActive and dbg) printf(" GUIDE: hit grease\n");
if (!hitActive) propagation = noCoupling->PropagateRay(ray0, &ray1, &vec1);
}
// check on which side the vector is?
TVector3 ray = ray1.GetN();
TVector3 exitNormal = s_side[5]->GetN();
if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal);
if (ray*exitNormal > 0) {
if (dbg) printf(" GUIDE: ray is backreflected from exit window.\n");
fate = backreflected;
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
}
fate = hitExit;
hout->Fill(vec1.y(), vec1.z());
hnodb_exit->Fill(n_odb-1);
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
}
 
if(propagation == REFRACTION) {
fate = refracted;
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
} // no total reflection when should be
 
last_hit = inters_i;
}
}
 
//--- material absorption ---
if(_absorption) {
double travel = 0.0;
if (dbg) printf("n_odb = %d\n", n_odb);
for(int point = 0; point < n_odb-1; point++) {
travel += (points[point] - points[point+1]).Mag();
if (dbg) printf("travel = %lf\n", travel);
}
double T_abs = TMath::Exp(-travel/_A);
if(dbg)printf("T_abs = %lf\n", T_abs);
double p_abs = rand.Uniform(0.0, 1.0);
if(dbg)printf("p_abs = %lf\n", p_abs);
 
if(p_abs > T_abs) fate = absorbed; // absorption
}
//--- material absorption ---
 
hfate->Fill(fate);
hfate->Fill(rays);
hnodb_all->Fill(n_odb-2);
*n_points = n_odb+1;
*out = ray0;
return fate;
}
//-----------------------------------------------------------------------------
void Guide::GetVFate(int *out)
{
for(int i=0;i<7;i++) out[i] = (int)hfate->GetBinContent(i+1);
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);
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);
line3d->SetLineWidth(width); line3d->SetLineColor(color);
TPolyLine3D *line3d = new TPolyLine3D(2);
line3d->SetLineWidth(width); line3d->SetLineColor(color);
 
for(int i=0; i<4; i++) {
line3d->SetPoint(0, vodnik_edge[i+0].x(), vodnik_edge[i+0].y(), vodnik_edge[i+0].z());
line3d->SetPoint(1, vodnik_edge[i+4].x(), vodnik_edge[i+4].y(), vodnik_edge[i+4].z());
line3d->DrawClone();
}
for(int i=0; i<4; i++) {
line3d->SetPoint(0, vodnik_edge[i+0].x(), vodnik_edge[i+0].y(), vodnik_edge[i+0].z());
line3d->SetPoint(1, vodnik_edge[i+4].x(), vodnik_edge[i+4].y(), vodnik_edge[i+4].z());
line3d->DrawClone();
}
}
//=================================================================================
 
813,60 → 833,60
//=================================================================================
int CPlaneR::TestIntersection(TVector3 *vec, CRay ray)
{
double num, den; //stevec, imenovalec
double t;
TVector3 tmp;
double num, den; //stevec, imenovalec
double t;
TVector3 tmp;
 
if(dbg) printf("---> CPlaneR::TestIntersection <---\n");
if(dbg) {printf("c = "); printv(center); printf(" | n = "); printv(n); printf("\n");}
double D = - n*center;
num = n*ray.GetR() + D;
den = n*ray.GetN();
if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den);
if(dbg) printf("---> CPlaneR::TestIntersection <---\n");
if(dbg) {printf("c = "); printv(center); printf(" | n = "); printv(n); printf("\n");}
 
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;
if(dbg) {printv(tmp); printf(" | Rv = %.4lf <> R = %.4lf\n", ((tmp - center).Mag()), _r);}
if( ((tmp - center).Mag()) < _r )
return 1;
else
return 0;
double D = - n*center;
num = n*ray.GetR() + D;
den = n*ray.GetN();
 
if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den);
 
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;
 
if(dbg) {printv(tmp); printf(" | Rv = %.4lf <> R = %.4lf\n", ((tmp - center).Mag()), _r);}
 
 
if( ((tmp - center).Mag()) < _r )
return 1;
else
return 0;
}
//-----------------------------------------------------------------------------
void CPlaneR::Draw(int color, int width)
{
const int NN = 32;
double phi, x, y;
TPolyLine3D *arc;
arc = new TPolyLine3D(NN+1);
arc->SetLineWidth(width);
arc->SetLineColor(color);
const int NN = 32;
double phi, x, y;
 
for(int i=0; i<=NN; i++) {
phi = i*2.0*TMath::Pi()/NN;
x = _r*TMath::Cos(phi);
y = _r*TMath::Sin(phi);
arc->SetPoint(i, center.x(), x, y);
}
arc->Draw();
TPolyLine3D *arc;
arc = new TPolyLine3D(NN+1);
arc->SetLineWidth(width);
arc->SetLineColor(color);
 
for(int i=0; i<=NN; i++) {
phi = i*2.0*TMath::Pi()/NN;
x = _r*TMath::Cos(phi);
y = _r*TMath::Sin(phi);
arc->SetPoint(i, center.x(), x, y);
}
arc->Draw();
}
//=================================================================================
 
873,109 → 893,104
 
//=================================================================================
CDetector::CDetector(TVector3 center0, DetectorParameters& parameters) :
center(center0),
glass_on(parameters.getGlassOn()),
glass_d(parameters.getGlassD()),
//x_gap(parameters.getGap().X()),
//y_gap(parameters.getGap().Y()),
//z_gap(parameters.getGap().Z()),
//glass(new CSurface),
//glass_circle(new CPlaneR),
//active(new CPlane4),
col_in(2),
col_lg(8),
col_out(4),
col_rgla(6),
col_LG(1),
col_glass(4),
col_active(7),
guide_on(parameters.getGuideOn()),
//window_R( parameters.getB() ),
//window_d(0),
guide(new Guide(center0, parameters)),
plate(new Plate(parameters)),
_plateWidth(parameters.getPlateWidth()),
_plateOn(parameters.getPlateOn()),
offsetY(parameters.getOffsetY()),
offsetZ(parameters.getOffsetZ())
{
// };
//-----------------------------------------------------------------------------
//void CDetector::Init()
//{
center(center0),
glass_on(parameters.getGlassOn()),
glass_d(parameters.getGlassD()),
col_in(2),
col_lg(8),
col_out(4),
col_rgla(6),
col_LG(1),
col_glass(4),
col_active(7),
guide_on(parameters.getGuideOn()),
guide(new Guide(center0, parameters)),
plate(new Plate(parameters)),
_plateWidth(parameters.getPlateWidth()),
_plateOn(parameters.getPlateOn()),
offsetY(parameters.getOffsetY()),
offsetZ(parameters.getOffsetZ())
{
// };
 
//-----------------------------------------------------------------------------
//void CDetector::Init()
//{
double d = parameters.getD();
double x_offset;
if(guide_on) x_offset = center.x();
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 n2 = parameters.getN2();
double n3 = parameters.getN3();
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.57;
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, 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);
detector = new CPlane4(plane_v);
hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector;
//hdetector = new TH2F("hdetector", "Hits detector plane", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size);
hdetector = new TH2F("hdetector", ";x [mm]; y [mm]", nBins, y_gap-p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ);
/*
double x_offset;
if(guide_on) x_offset = center.x();
else x_offset = center.x() - d;
 
double b = parameters.getB();
//double n1 = parameters.getN1();
//double n2 = parameters.getN2();
double n3 = parameters.getN3();
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.57;
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, 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", "",
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;
if (dbg) cout<<"SiPM active length "<<parameters.getActive()<<endl;
 
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);
//active surface in case of bad coupling is circle d=a
TVector3 activePosition(center);
activePosition += TVector3(d + x_gap, 0, 0);
TVector3 normal(1,0,0);
grease = new CPlaneR(activePosition, normal, 1.0*p_size);
 
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", ";x [mm];y [mm]", 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);
detector = new CPlane4(plane_v);
 
hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector;
//hdetector = new TH2F("hdetector", "Hits detector plane", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size);
hdetector = new TH2F("hdetector", ";x [mm]; y [mm]", nBins, y_gap-p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ);
 
/*
window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R);
 
p_size = M*a;
plane_v[0].SetXYZ(x_offset+d+window_d, y_gap + p_size, z_gap - p_size);
plane_v[1].SetXYZ(x_offset+d+window_d, y_gap + p_size, z_gap + p_size);
982,13 → 997,13
plane_v[2].SetXYZ(x_offset+d+window_d, y_gap - p_size, z_gap + p_size);
plane_v[3].SetXYZ(x_offset+d+window_d, y_gap - p_size, z_gap - p_size);
window = new CSurface(SURF_REFRA, plane_v, n1, n2, reflectivity); window->FlipN();
 
hwindow = (TH2F*)gROOT->FindObject("hwindow"); if(hwindow) delete hwindow;
hwindow = new TH2F("hwindow", "Hits Window", nch, y_gap - window_R, y_gap + window_R, nch, z_gap - window_R, z_gap + window_R);
*/
p_size = b/2.0;
histoPlate = (TH2F*)gROOT->FindObject("histoPlate"); if(histoPlate) delete histoPlate;
histoPlate = new TH2F("histoPlate", "Hits on glass plate", nBins, -p_size, +p_size, nBins, -p_size, +p_size);
*/
p_size = b/2.0;
histoPlate = (TH2F*)gROOT->FindObject("histoPlate"); if(histoPlate) delete histoPlate;
histoPlate = new TH2F("histoPlate", "Hits on glass plate", nBins, -p_size, +p_size, nBins, -p_size, +p_size);
}
 
//-----------------------------------------------------------------------------
995,7 → 1010,6
// vrne 1 ce je zadel aktvino povrsino
// vrne <1 ce jo zgresi
int CDetector::Propagate(CRay in, CRay *out, int draw)
 
// Sledi zarku skozi vodnik. Vrne:
// 0, ce zgresi vstopno ploskev MISSED
// 1, ce zadane izstopno ploskev HIT
1006,195 → 1020,195
// +4, ce se absorbira v materialu ABSORBED
{
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(in);
rayout->SetColor(col_in);
//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(in);
rayout->SetColor(col_in);
 
const int max_n_points = guide->GetMAXODB() + 2;
TVector3 pointsPlate[max_n_points];
//TVector3 intersection;
Fate fatePlate;
int nPointsPlate;
TPolyLine3D *line3d = new TPolyLine3D(2);
line3d->SetLineWidth(1);
line3d->SetLineColor(4);
// Draw the plate and propagate the ray through
// 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 {
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) { // lost on plate side
rayout->SetColor(col_out);
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 (fatePlate == backreflected)
guide->GetHFate()->Fill(fatePlate); // reflected back
else
guide->GetHFate()->Fill(noreflection); //lost on plate side
return fatePlate;
}
//Ray hits light guide
histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
const int max_n_points = guide->GetMAXODB() + 2;
TVector3 pointsPlate[max_n_points];
//TVector3 intersection;
Fate fatePlate;
int nPointsPlate;
TPolyLine3D *line3d = new TPolyLine3D(2);
line3d->SetLineWidth(1);
line3d->SetLineColor(4);
 
}
else {
//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
//const int max_n_points = guide->GetMAXODB() + 2;
TVector3 points[max_n_points];
TVector3 presecisce;
// Draw the plate and propagate the ray through
// check if the ray should be reflected??
 
int n_points;
int fate_glass;
CRay *ray0 = new CRay(*rayout);
// delete rayout; -> creates dangling reference when tries to delete ray0!
//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");
}
line3d->SetLineColor(col_lg);
int p_i;
if(guide_on) {
if(draw) {
if(fate == missed) {
if (dbg) printf("Detector: fate=missed\n");
TVector3 r = ray1->GetR();
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++) {
line3d->SetPoint(0, points[p_i].x(), points[p_i].y(), points[p_i].z());
line3d->SetPoint(1, points[p_i+1].x(), points[p_i+1].y(), points[p_i+1].z());
line3d->DrawClone();
}
if(fate != noreflection) {
if (dbg) printf("Detector: fate != noreflection, fate = %d\n", (int)fate);
if(glass_on) {/*if(fate == 1)*/ ray1->Draw(points[p_i].x(), center.x() + guide->getD() + glass_d);}
else {
ray1->SetColor(col_out);
ray1->DrawS(points[p_i].x(), 10.0);
}
}
}
}
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(draw) {
//double epoxy = parameters->getGlassD();
if(glass_on) ray1->Draw(center.x(), center.x() + glass_d);
else ray1->DrawS(center.x(), 10.0);
}
}
fate = missed; // zgresil aktivno povrsino
if(glass_on) {
*ray0 = *ray1;
ray1->SetColor(col_rgla);
fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
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() + 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);
}
}
if(_plateOn) {
 
// 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 {
//if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
//}
*out = *ray1;
delete ray0;
delete ray1;
delete rayout;
delete rayin;
return fate;
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 {
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) { // lost on plate side
rayout->SetColor(col_out);
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 (fatePlate == backreflected)
guide->GetHFate()->Fill(fatePlate); // reflected back
else
guide->GetHFate()->Fill(noreflection); //lost on plate side
return fatePlate;
}
 
//Ray hits light guide
histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
 
}
else {
//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
 
//const int max_n_points = guide->GetMAXODB() + 2;
TVector3 points[max_n_points];
TVector3 presecisce;
 
int n_points;
int fate_glass;
CRay *ray0 = new CRay(*rayout);
// delete rayout; -> creates dangling reference when tries to delete ray0!
//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");
}
 
line3d->SetLineColor(col_lg);
int p_i;
if(guide_on) {
if(draw) {
if(fate == missed) {
if (dbg) printf("Detector: fate=missed\n");
TVector3 r = ray1->GetR();
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++) {
line3d->SetPoint(0, points[p_i].x(), points[p_i].y(), points[p_i].z());
line3d->SetPoint(1, points[p_i+1].x(), points[p_i+1].y(), points[p_i+1].z());
line3d->DrawClone();
}
if(fate != noreflection) {
if (dbg) printf("Detector: fate != noreflection, fate = %d\n", (int)fate);
if(glass_on) {/*if(fate == 1)*/ ray1->Draw(points[p_i].x(), center.x() + guide->getD() + glass_d);}
else {
ray1->SetColor(col_out);
ray1->DrawS(points[p_i].x(), 10.0);
}
}
}
}
 
 
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(draw) {
//double epoxy = parameters->getGlassD();
if(glass_on) ray1->Draw(center.x(), center.x() + glass_d);
else ray1->DrawS(center.x(), 10.0);
}
}
 
fate = missed; // zgresil aktivno povrsino
if(glass_on) {
*ray0 = *ray1;
ray1->SetColor(col_rgla);
fate_glass = glass->PropagateRay(*ray0, ray1, &presecisce);
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() + 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);
}
}
 
// 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 {
//if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
//}
 
*out = *ray1;
delete ray0;
delete ray1;
delete rayout;
delete rayin;
return fate;
}
//-----------------------------------------------------------------------------
void CDetector::Draw(int width)
{
if(guide_on) {
if( TMath::Abs(guide->getN1()-guide->getN2()) < MARGIN ) {
if(_plateOn) plate->drawSkel(col_LG, width);
guide->DrawSkel(col_LG, width);
}
else {
if(_plateOn) plate->draw(4, width);
guide->Draw(col_LG, width);
}
}
if(glass_on) glass_circle->Draw(col_glass, width);
//window_circle->Draw(col_glass, width);
active->Draw(col_active, width);
if(guide_on) {
if( TMath::Abs(guide->getN1()-guide->getN2()) < MARGIN ) {
if(_plateOn) plate->drawSkel(col_LG, width);
guide->DrawSkel(col_LG, width);
}
else {
if(_plateOn) plate->draw(4, width);
guide->Draw(col_LG, width);
}
}
 
if(glass_on) glass_circle->Draw(col_glass, width);
//window_circle->Draw(col_glass, width);
active->Draw(col_active, width);
}
//=================================================================================
 
1203,52 → 1217,52
TVector3 center = CENTER;
const double b = parameters.getB();
const double n1 = parameters.getN1();
const double n2 = parameters.getN2();
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[2].SetXYZ(0.0,-t, t);
plate_edge[3].SetXYZ(0.0,-t,-t);
plate_edge[4].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);
for(int i = 0; i<8; i++) plate_edge[i] += center;
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, c_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, c_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, n2, c_reflectivity);
sides[5]->FlipN();
for(int i=0; i<6; i++) sides[i]->SetFresnel(1);
const double n2 = parameters.getN2();
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[2].SetXYZ(0.0,-t, t);
plate_edge[3].SetXYZ(0.0,-t,-t);
plate_edge[4].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);
 
for(int i = 0; i<8; i++) plate_edge[i] += center;
 
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, c_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, c_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, n2, c_reflectivity);
sides[5]->FlipN();
 
for(int i=0; i<6; i++) sides[i]->SetFresnel(1);
}
 
void Plate::draw(int color, int width)
{
for(int i = 0; i<6; i++) sides[i]->Draw(color, width);
for(int i = 0; i<6; i++) sides[i]->Draw(color, width);
}
 
void Plate::drawSkel(int color, int width)
{
TPolyLine3D line3d(2);
line3d.SetLineWidth(width);
line3d.SetLineColor(color);
TPolyLine3D line3d(2);
line3d.SetLineWidth(width);
line3d.SetLineColor(color);
 
for(int i=0; i<4; i++) {
line3d.SetPoint(0, plate_edge[i+0].x(), plate_edge[i+0].y(), plate_edge[i+0].z());
line3d.SetPoint(1, plate_edge[i+4].x(), plate_edge[i+4].y(), plate_edge[i+4].z());
line3d.DrawClone();
}
for(int i=0; i<4; i++) {
line3d.SetPoint(0, plate_edge[i+0].x(), plate_edge[i+0].y(), plate_edge[i+0].z());
line3d.SetPoint(1, plate_edge[i+4].x(), plate_edge[i+4].y(), plate_edge[i+4].z());
line3d.DrawClone();
}
}
 
Fate Plate::propagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
1258,66 → 1272,66
TVector3 vec0, vec1;
Fate fate = enter;
int inters_i = 0;
ray0 = in;
int n_odb = 0;
int last_hit = 0;
int propagation = 0;
int result = sides[0]->PropagateRay(ray0, &ray1, &vec1);
if( !result ) {
// ce -NI- presecisca z vstopno
fate = missed;
} else if(result == REFLECTION) {
if (dbg) printf("PLATE: reflected\n");
fate = backreflected;
} else {
points[0] = ray1.GetR();
//hfate->Fill(enter);
//hin->Fill(vec1.y(), vec1.z());
while (n_odb++ < MAX_REFLECTIONS) {
ray0 = ray1;
vec0 = vec1;
propagation = 11;
for(inters_i=0; inters_i<6; inters_i++) {
if( inters_i != last_hit) {
if( sides[inters_i]->TestIntersection(&vec1, ray1) ) break;
}
}
points[n_odb] = vec1;
if(inters_i == 0) {
fate = backreflected;
break;} // backreflection
propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
if(inters_i == 5) { // successfull exit
fate = hitExit;
//hout->Fill(vec1.y(), vec1.z());
//hnodb_exit->Fill(n_odb-1);
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
}
if(propagation == 1) {
fate = noreflection; //at side
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;} // no total reflection when should be
if(propagation == -2) {
fate = noreflection;
break;
} // absorption due to finite reflectivity
last_hit = inters_i;
}
}
*n_points = n_odb+1;
*out = ray0;
return fate;
 
ray0 = in;
int n_odb = 0;
int last_hit = 0;
int propagation = 0;
 
int result = sides[0]->PropagateRay(ray0, &ray1, &vec1);
if( !result ) {
// ce -NI- presecisca z vstopno
fate = missed;
} else if(result == REFLECTION) {
if (dbg) printf("PLATE: reflected\n");
fate = backreflected;
} else {
points[0] = ray1.GetR();
//hfate->Fill(enter);
//hin->Fill(vec1.y(), vec1.z());
while (n_odb++ < MAX_REFLECTIONS) {
ray0 = ray1;
vec0 = vec1;
propagation = 11;
for(inters_i=0; inters_i<6; inters_i++) {
if( inters_i != last_hit) {
if( sides[inters_i]->TestIntersection(&vec1, ray1) ) break;
}
}
points[n_odb] = vec1;
if(inters_i == 0) {
fate = backreflected;
break;} // backreflection
 
propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
if(inters_i == 5) { // successfull exit
fate = hitExit;
//hout->Fill(vec1.y(), vec1.z());
//hnodb_exit->Fill(n_odb-1);
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
}
if(propagation == 1) {
fate = noreflection; //at side
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;} // no total reflection when should be
 
if(propagation == -2) {
fate = noreflection;
break;
} // absorption due to finite reflectivity
 
last_hit = inters_i;
}
}
 
*n_points = n_odb+1;
*out = ray0;
return fate;
};
//=============================================================================================================================== <<<<<<<<