Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 84 → Rev 83

/lightguide/trunk/include/guide.h
5,7 → 5,7
#include "TVector3.h"
#include "TMath.h"
#include "TPolyLine3D.h"
#include "TRandom2.h"
#include "TRandom.h"
#include "TDatime.h"
#include "TH1F.h"
#include "TH2F.h"
40,21 → 40,22
CRay() :
r(TVector3(0,0,0)),
k(TVector3(1,0,0)),
p(TVector3(0,1,0)),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
CRay(TVector3 r0, TVector3 k0) :
r(r0),
k(k0.Unit()),
p(TVector3(0,1,0)),
//p(TVector3(0,0,1)),
//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)),
k(TVector3(l0,m0,n0).Unit()),
p(TVector3(0,1,0)),
//p(TVector3(0,0,1)),
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
 
66,7 → 67,7
if (p.Dot(k)>1e-3) printf("*** ERROR in CRay: E has component || with k\n");
};
 
inline CRay & operator = (const CRay &);
//inline CRay & operator = (const CRay &);
 
TVector3 GetR() const {return r;};
TVector3 GetK() const {return k;};
102,7 → 103,6
CPlane4();
CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);
CPlane4(TVector3 *vr);
~CPlane4() { };
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;};
119,8 → 119,7
void Draw(int color = 1, int width = 1);
 
private:
TVector3 r[4];
TVector3 n;
TVector3 r[4], n;
double A, B, C, D;
TVector3 edge[4]; // vektorji stranic
double angle_r[4]; // koti ob posameznem vogalu
133,8 → 132,6
public:
CPlaneR(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
CPlaneR(){ n=TVector3(0,0,0); center=TVector3(0,0,0); _r=0; };
~CPlaneR() {};
 
void Set(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
166,17 → 163,6
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() {};
CSurface& operator=(const CSurface& rhs) {
type=rhs.type;
n1=rhs.n1;
n2=rhs.n2;
n1_n2=rhs.n1_n2;
reflection=rhs.reflection;
cosTtotal=rhs.cosTtotal;
fresnel=rhs.fresnel;
return *this;
}
 
void SetV(TVector3 *vr){CPlane4::Set(vr);};
void SetType(int type0){type = type0;};
189,7 → 175,7
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);
int PropagateRay(CRay in, CRay *out, TVector3 *intersection);
 
double N1_N2(int sign) { return ((sign > 0) ? n1/n2 : n2/n1);};
 
199,7 → 185,7
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
TRandom2 rand; // za racunanje verjetnosti odboja od zrcala
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
332,28 → 318,29
public:
Guide(TVector3 center0, DetectorParameters& parameters);
~Guide() {
//for (int jk=0; jk<6; jk++) delete s_side[jk];
//delete [] s_side;
//delete hfate;
//delete hnodb_all;
//delete hnodb_exit;
//delete hin;
//delete hout;
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);
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() {return &hfate;};
TH1F* GetHNOdb_all() {return &hnodb_all;};
TH1F* GetHNOdb_exit() {return &hnodb_exit;};
TH2F* GetHIn() {return &hin;};
TH2F* GetHOut() {return &hout;};
int GetExitHits() {return (int)hfate.GetBinContent(5);};
int GetEnteranceHits() {return (int)hfate.GetBinContent(6);};
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;};
 
368,18 → 355,14
double _absorption;
double _A;
bool _badCoupling;
TRandom2 rand; // for material absorption
CSurface s_side[6];
CPlaneR grease;
CSurface noCoupling;
TVector3 center;
TVector3 vodnik_edge[8];
TRandom rand; // for material absorption
CSurface *s_side[6];
CPlaneR* grease;
CSurface* noCoupling;
TVector3 center, vodnik_edge[8];
 
TH1F hfate;
TH1F hnodb_all;
TH1F hnodb_exit;
TH2F hin;
TH2F hout;
TH1F *hfate, *hnodb_all, *hnodb_exit;
TH2F *hin, *hout;
};
 
 
388,15 → 371,17
{
public:
Plate(DetectorParameters& parameters);
~Plate() { };
~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*);
Fate propagateRay(CRay, CRay*, int*, TVector3*);
 
private:
TVector3 plate_edge[8];
CSurface sides[6];
CSurface *sides[6];
};
 
 
405,17 → 390,17
public:
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 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;
452,16 → 437,15
 
//void Init();
 
int Propagate(CRay, CRay&, int);
int Propagate(CRay, CRay*, int);
 
Guide* GetLG() {return &guide;};
Guide* GetLG() const {return guide;};
//TH2F* GetHWindow() const {return hwindow;};
TH2F* GetHGlass() {return &hglass;};
TH2F* GetHActive() {return &hactive;};
TH2F* GetHLaser() {return &hlaser;};
TH2F* GetHDetector() {return &hdetector;};
TH2F* GetHPlate() {return &histoPlate;};
TH2F* GetGenerated() {return &generated;};
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;};
486,18 → 470,17
 
int glass_on;
double glass_d;
CSurface glass;
CPlaneR glass_circle;
TH2F hglass;
CSurface *glass;
CPlaneR *glass_circle;
TH2F *hglass;
 
//double x_gap, y_gap, z_gap;
CPlane4 active;
CPlaneR grease;
TH2F hactive;
TH2F hlaser;
CPlane4* active;
CPlaneR* grease;
TH2F *hactive, *hlaser;
 
CPlane4 detector;
TH2F hdetector;
CPlane4 *detector;
TH2F *hdetector;
 
int col_in, col_lg, col_out, col_rgla;
int col_LG, col_glass, col_active;
512,18 → 495,16
//CPlaneR *window_circle;
//TH2F *hwindow;
 
Guide guide;
Plate plate;
Guide *guide;
Plate *plate;
 
double _plateWidth;
int _plateOn;
 
TH2F histoPlate;
TH2F *histoPlate;
 
double offsetY;
double offsetZ;
 
TH2F generated;
};
 
 
/lightguide/trunk/src/userFunctions.cpp
128,18 → 128,16
 
double cx = 0.0;
TVector3 vodnik_edge[4];
double t = 4.0;
double t = 3.0;
vodnik_edge[0].SetXYZ(cx, t,-t);
vodnik_edge[1].SetXYZ(cx, t, t);
vodnik_edge[2].SetXYZ(cx,-t, t);
vodnik_edge[3].SetXYZ(cx,-t,-t);
 
CSurface surf;
//surf = CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96);
surf.Set(p_type, vodnik_edge, p_n1, p_n2, 0.96);
surf.FlipN();
surf.SetFresnel(1);
surf.Draw();
CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96);
surf->FlipN();
surf->SetFresnel(1);
surf->Draw();
 
CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
ray->SetColor(kRed);
153,13 → 151,13
 
ray->DrawS(cx, -5.0);
 
CRay out;
out.SetColor(kBlack);
CRay *out = new CRay;
out->SetColor(kBlack);
TVector3 *inters = new TVector3;
surf.PropagateRay(*ray, out, inters);
surf->PropagateRay(*ray, out, inters);
printf(" n1 = %f, n2 = %f\n", p_n1, p_n2);
//if(fate == 1) out->DrawS(cx, 5.0);
out.DrawS(cx, 5.0);
out->DrawS(cx, 5.0);
 
CRay *incidentPolarization = new CRay;
incidentPolarization->SetColor(kGreen);
169,7 → 167,7
 
CRay *surfaceNormal = new CRay;
surfaceNormal->SetColor(kBlue);
surfaceNormal->Set(ray->GetR(), surf.GetN());
surfaceNormal->Set(ray->GetR(), surf->GetN());
surfaceNormal->DrawS(cx, 1.0);
printf(" BLUE: surface normal vector\n");
}
/lightguide/trunk/src/raySimulator.cpp
96,42 → 96,18
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
cdata->Divide(3,3);
int cpc = 1;
cdata->cd(cpc++);
TH2F* generated = detector->GetGenerated();
int nGenerated = generated->GetEntries();
int minimum = generated->GetBinContent(generated->GetMinimumBin());
int maximum = generated->GetBinContent(generated->GetMaximumBin());
generated->GetZaxis()->SetRangeUser(0, 1.05*maximum);
double variation = (maximum-minimum)/(double)nGenerated;
printf("Statistical variation (max-min)/all = %f perc. \n", variation*100);
generated->SetTitle("Generated");
generated->Draw("colz");
cdata->cd(cpc++); //((detector->GetLG())->GetHIn())->Draw("COLZ");
TH2F* histoPlate = detector->GetHPlate();
histoPlate->Draw("COLZ");
cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ");
TH2F* histoGlass = (detector->GetHGlass());
histoGlass->Draw("COLZ");
(detector->GetHGlass())->Draw("COLZ");
 
cdata->cd(cpc++);
TH2F* histoActive = (detector->GetHActive());
histoActive->Draw("COLZ");
cdata->cd(cpc++);
TH2F* histoLaser = (detector->GetHLaser());
histoLaser->Draw("COLZ");
cdata->cd(cpc++);
TH2F* histoDetector = (detector->GetHDetector());
histoDetector->Draw("COLZ");
cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
 
cdata->cd(cpc++);
TH1F* fate=((detector->GetLG())->GetHFate());
fate->Draw();
cdata->cd(cpc++);
TH1F* reflections = ((detector->GetLG())->GetHNOdb_all());
reflections->Draw();
cdata->cd(cpc++);
TH1F *reflectedExit = ((detector->GetLG())->GetHNOdb_exit());
reflectedExit->Draw();
cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
}else{
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
203,7 → 179,7
//-----------------------------------------------------------------------------
// en zarek
double Single(CDetector *detector, DetectorParameters& parameters,
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
{
theta = Pi()*theta/180.0;
if(theta < 1e-6) theta = 1e-6;
218,14 → 194,13
CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi),
TMath::Sin(theta)*TMath::Cos(phi));
// Set y-polarization == vertical
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
if (polarization.Dot(ray0->GetK()) > 1e-2) printf("ERROR: pol not perep\n");
// Set z-polarization == vertical
TVector3 polarization(0,0,1);
polarization.RotateY(-theta);
polarization.RotateZ(phi);
if (polarization.Dot(ray0->GetK()) > 1e-5) printf("ERROR: pol not perep\n");
ray0->setPolarization(polarization);
CRay ray1;
CRay *ray1 = new CRay;
 
detector->Propagate(*ray0, ray1, show_3d);
 
236,8 → 211,8
incidentPolarization->Set(drawPosition, polarization);
incidentPolarization->DrawS(drawPosition.X(), 1);
 
TVector3 outPolarization = ray1.GetP();
drawPosition = ray1.GetR();
TVector3 outPolarization = ray1->GetP();
drawPosition = ray1->GetR();
drawPosition.SetX(drawPosition.X() + 5);
CRay* rayPol = new CRay(drawPosition, outPolarization);
rayPol->SetColor(kBlack);
244,7 → 219,7
rayPol->DrawS(drawPosition.X(), 1);
 
delete ray0;
//delete ray1;
delete ray1;
 
return (detector->GetHActive())->GetEntries() / (double)(1);
}
251,7 → 226,7
//-----------------------------------------------------------------------------
// zarki, razporejeni v mrezi
double Grid(CDetector *detector, DetectorParameters& parameters,
int NN = 10, double theta = 0.0)
int NN = 10, double theta = 0.0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = Pi()*theta/180.0;
261,32 → 236,27
if(show_3d) detector->Draw(draw_width);
 
const double b = parameters.getB();
double simulated = 0;
 
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
const int GRID = 100;
for(int i = 0; i < GRID; i++) {
for(int j = 0; j < GRID; j++) {
for(int jk = 0; jk<NN; jk++){
CRay *ray0 = new CRay(CENTER.x() - offset,
(i*b/GRID + b/(2.0*GRID) - b/2.0),
(j*b/GRID + b/(2.0*GRID) - b/2.0),
TMath::Cos(theta),
0.0,
TMath::Sin(theta));
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
ray0->setPolarization(polarization);
CRay ray1;
if(i == (GRID/2) and jk == 0)
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
delete ray0;
simulated++;
//delete ray1;
}
 
for(int i = 0; i < NN; i++) {
for(int j = 0; j < NN; j++) {
CRay *ray0 = new CRay(CENTER.x() - offset,
0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
TMath::Cos(theta),
0.0,
TMath::Sin(theta));
TVector3 polarization(0, 1, 0);
polarization.RotateY(-theta);
ray0->setPolarization(polarization);
CRay *ray1 = new CRay;
if(i == (NN/2))
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
delete ray0;
delete ray1;
}
 
}
297,7 → 267,7
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / simulated;
acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
return acceptance;
}
//-----------------------------------------------------------------------------
304,7 → 274,7
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
// vsi pod kotom (theta, phi)
double RandYZ(CDetector *detector, DetectorParameters& parameters,
int NN, double theta, double phi, int show_rays)
int NN, double theta, double phi, int show_rays)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = theta*3.14159265358979312/180.0;
313,20 → 283,21
if(phi < MARGIN) phi = MARGIN;
 
TDatime now;
TRandom2 rand;
TRandom rand;
rand.SetSeed(now.Get());
 
//detector->Init();
if(show_3d) detector->Draw(draw_width);
 
double b = parameters.getB();
double SiPM = parameters.getA();
double M = parameters.getM();
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 
for(int i = 0; i < NN; i++) {
 
double rx = CENTER.x() - offset;
double ry = rand.Uniform(-b/2.0, b/2.0);
double rz = rand.Uniform(-b/2.0, b/2.0);
double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
 
double rl = TMath::Cos(theta);
double rm = TMath::Sin(theta)*TMath::Sin(phi);
333,18 → 304,17
double rn = TMath::Sin(theta)*TMath::Cos(phi);
 
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
TVector3 polarization(0, 0, 1);
polarization.RotateY(-theta);
polarization.RotateZ(phi);
ray0->setPolarization(polarization);
CRay ray1;
CRay *ray1 = new CRay;
 
if(i < show_rays)
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
//delete ray1;
delete ray1;
delete ray0;
}
 
366,20 → 336,20
// theta [0, theta]
// phi [0,360]
double RandIso(CDetector *detector, DetectorParameters& parameters,
int NN = 1e3, double thetaMin=0.0, double thetaMax = 30.0, int show_rays = 30, int show_rand = 0)
int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
double pi = 3.14159265358979312;
thetaMax = thetaMax*pi/180.0;
thetaMin = thetaMin*pi/180.0;
if(thetaMin < MARGIN) thetaMin = MARGIN;
if(thetaMax < MARGIN) thetaMax = MARGIN;
theta = theta*3.14159265358979312/180.0;
if(theta < MARGIN) theta = MARGIN;
 
TDatime now;
TRandom2 rand;
TRandom rand;
rand.SetSeed(now.Get());
double rx, ry, rz, rl, rm, rn;
double rphi, rtheta;
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
double theta_min_rad = TMath::Cos(thetaMax);
double theta_min_rad = TMath::Cos(theta);
double theta_max_rad = TMath::Cos(thetaMin);
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
400,40 → 370,34
//detector->Init();
if (show_3d) detector->Draw(draw_width);
 
double b = parameters.getB();
double SiPM = parameters.getA();
double M = parameters.getM();
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 
for(int i = 0; i < NN; i++) {
 
double rx = CENTER.x() - offset;
double ry = rand.Uniform(-b/2.0, b/2.0);
double rz = rand.Uniform(-b/2.0, b/2.0);
rx = CENTER.x() - offset;
ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
 
double phi = rand.Uniform(0.0, 2.0*pi);
rphi = rand.Uniform(0.0, 2.0*pi);
//double theta_max_rad = TMath::Cos(17.0*pi/180.0);
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
double theta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
 
double rl = TMath::Cos(theta);
double rm = TMath::Sin(theta)*TMath::Sin(phi);
double rn = TMath::Sin(theta)*TMath::Cos(phi);
rl = TMath::Cos(rtheta);
rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
 
if(show_rand) {
htheta->Fill(theta);
hcostheta->Fill( TMath::Cos(theta) );
hphi->Fill(phi);
hl->Fill(rl);
hm->Fill(rm);
hn->Fill(rn);
htheta->Fill(rtheta);
hcostheta->Fill( TMath::Cos(rtheta) );
hphi->Fill(rphi);
hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
}
 
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
ray0->setPolarization(polarization);
CRay ray1;
CRay *ray1 = new CRay;
 
if(i < show_rays) {
detector->Propagate(*ray0, ray1, show_3d);
442,8 → 406,8
detector->Propagate(*ray0, ray1, 0);
}
 
//delete ray1;
delete ray0;
delete ray1;
}
 
if(show_rand) {
474,7 → 438,7
// Beamtest distribution
// with fixed theta and phi in interval phiMin, phiMax
double beamtest(CDetector *detector, DetectorParameters& parameters,
int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
{
double pi = 3.14159265358979312;
theta *= pi/180.0;
483,10 → 447,10
phiMax *= pi/180.0;
 
TDatime now;
TRandom2 rand;
TRandom rand;
rand.SetSeed(now.Get());
double rx, ry, rz, rl, rm, rn;
double phi;
double rphi;
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
514,19 → 478,19
ry = rand.Uniform(-b/2.0, b/2.0);
rz = rand.Uniform(-b/2.0, b/2.0);
 
phi = rand.Uniform(phiMin, phiMax);
rphi = rand.Uniform(phiMin, phiMax);
//double theta_max_rad = TMath::Cos(17.0*pi/180.0);
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
 
rl = TMath::Cos(theta);
rm = TMath::Sin(theta)*TMath::Sin(phi);
rn = TMath::Sin(theta)*TMath::Cos(phi);
rm = TMath::Sin(theta)*TMath::Sin(rphi);
rn = TMath::Sin(theta)*TMath::Cos(rphi);
 
if(show_rand) {
//htheta->Fill(rtheta);
//hcostheta->Fill( TMath::Cos(rtheta) );
hphi->Fill(phi);
hphi->Fill(rphi);
hl->Fill(rl);
hm->Fill(rm);
hn->Fill(rn);
534,12 → 498,11
 
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
// inital polarizaton
TVector3 k(ray0->GetK());
TVector3 polarization(k.Orthogonal());
polarization.Unit();
polarization.Rotate(phi, k);
TVector3 polarization(0, 0, 1);
polarization.RotateY(-theta);
polarization.RotateX(rphi);
ray0->setPolarization(polarization);
CRay ray1;
CRay *ray1 = new CRay;
 
if(i < show_rays) {
detector->Propagate(*ray0, ray1, show_3d);
548,8 → 511,8
detector->Propagate(*ray0, ray1, 0);
}
 
//delete ray1;
delete ray0;
delete ray1;
}
 
if(show_rand) {
/lightguide/trunk/src/guide.cpp
32,15 → 32,14
//n.SetXYZ(l0, m0, n0); n = n.Unit();
//}
//-----------------------------------------------------------------------------
 
CRay& CRay::operator = (const CRay& rhs)
/*
CRay& CRay::operator = (const CRay& p)
{
r.SetXYZ(rhs.GetR().x(), rhs.GetR().y(), rhs.GetR().z());
r.SetXYZ(p.GetR().x(), p.GetR().y(), p.GetR().z());
//this->r.SetXYZ(p.x(), p.y(), p.z());
k.SetXYZ(rhs.GetK().x(), rhs.GetK().y(), rhs.GetK().z());
p.SetXYZ(rhs.GetP().x(), rhs.GetP().y(), rhs.GetP().z());
n.SetXYZ(p.GetN().x(), p.GetN().y(), p.GetN().z());
return *this;
}
} */
void CRay::Print()
{
printf("---> CRay::Print() <---\n");
365,7 → 364,7
// sprejme zarek, vrne uklonjen/odbit zarek in presecisce
// vrne 0 ce ni presecisca; 1 ce se je lomil
// 2 ce se je odbil; -2 ce se je absorbiral
int CSurface::PropagateRay(CRay in, CRay& out, TVector3 *intersection)
int CSurface::PropagateRay(CRay in, CRay *out, TVector3 *intersection)
{
if (dbg) printf("--- CSurface::PropagateRay ---\n");
double cosTi; // incident ray angle
416,8 → 415,7
}
 
double cosAf = v_te * in.GetP();
double alpha = TMath::ACos(cosAf);
if(alpha < MARGIN) alpha = MARGIN;
double alpha = acos(cosAf);
if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE);
 
a_te = cosAf;
461,7 → 459,7
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);
out->Set(intersect, transmit);
 
// Shift implemented, but only linear polarization is implemented
if (dbg) printf("CSurface: Propagate TOTAL\n");
472,7 → 470,7
double deltaP = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/(pow(n12,2)*cosTi));
double deltaS = 2 * atan(sqrt(1 - cosTi*cosTi - pow(n12,2))/cosTi);
double delta = deltaP - deltaS;
//alpha += delta;
alpha += delta;
a_tm = sin(alpha);
a_te = cos(alpha);
if (dbg) printf(" deltaP = %f deltaS = %f; new a_tm = %f, a_te = %f",
479,7 → 477,7
deltaP, deltaS, a_tm, a_te);
pol_t = a_tm*v_tm_t + a_te*v_te;
if (dbg) printv(pol_t);
out.setPolarization(pol_t);
out->setPolarization(pol_t);
 
return REFLECTION;
} else {
523,8 → 521,8
 
if(p_ref >= R_f) { // se lomi
if (dbg) printf(" SURFACE REFRACTED. Return.\n");
out.Set(intersect, transmit);
out.setPolarization(pol_t);
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);
531,10 → 529,10
transmit = in.GetK() + sign_n*2*cosTi*n;
TVector3 v_tm_r = -v_te.Cross(transmit);
v_tm_r = v_tm_r.Unit();
out.Set(intersect, transmit);
out->Set(intersect, transmit);
//pol_t = -in.GetP() + sign_n*2*cosTi*n;
pol_t = a_te*(1.0 - TMath::Abs(r_te))*v_te + a_tm*(1.0 - TMath::Abs(r_tm))*v_tm_r;
out.setPolarization(pol_t);
out->setPolarization(pol_t);
return REFLECTION;
}
 
549,11 → 547,11
if(p_ref < reflection) { // se odbije
cosTi = in.GetK() * n;
transmit = in.GetK() - 2*cosTi*n;
out.Set(intersect, transmit);
out->Set(intersect, transmit);
return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
} else { // se ne odbije
transmit = in.GetK();
out.Set(intersect, transmit);
out->Set(intersect, transmit);
return ABSORBED;
}
break;
565,21 → 563,21
cosTi = in.GetK() * n;
if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj
transmit = in.GetK() - 2*cosTi*n;
out.Set(intersect, transmit);
out->Set(intersect, transmit);
} else { // ni tot. odboja
transmit = in.GetK();
out.Set(intersect, transmit);
out->Set(intersect, transmit);
return ABSORBED;
}
} else { // se ne odbije
transmit = in.GetK();
out.Set(intersect, transmit);
out->Set(intersect, transmit);
return ABSORBED;
}
break;
 
default:
out = in;
*out = in;
break;
}
 
622,55 → 620,55
for(int i = 0; i<8; i++) vodnik_edge[i] += center;
 
// light guide surfaces
s_side[0].Set(SURF_REFRA, vodnik_edge, n0, _n2, _r);
s_side[0].FlipN();
s_side[0] = new CSurface(SURF_REFRA, vodnik_edge, n0, _n2, _r);
s_side[0]->FlipN();
 
s_side[1].Set(SURF_REFRA, vodnik_edge[3], vodnik_edge[2],
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].Set(SURF_REFRA, vodnik_edge[2], vodnik_edge[1],
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].Set(SURF_REFRA, vodnik_edge[1], vodnik_edge[0],
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].Set(SURF_REFRA, vodnik_edge[0], vodnik_edge[3],
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].Set(SURF_REFRA, &vodnik_edge[4], _n2, _n3, _r);
s_side[5].FlipN();
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.Set(SURF_REFRA, &vodnik_edge[4], _n2, 1.0, _r);
noCoupling.FlipN();
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.Set(activePosition, normal, 0.95*a/2.0);
grease = new CPlaneR(activePosition, normal, 0.95*a/2.0);
 
if(fresnel) for(int i=0; i<6; i++) s_side[i].SetFresnel(1);
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 = TH1F("hfate", "Ray fate", 8, -3.5, 4.5);
(hfate.GetXaxis())->SetBinLabel(1, "Back Refl");
(hfate.GetXaxis())->SetBinLabel(2, "Side");
(hfate.GetXaxis())->SetBinLabel(3, "Refracted");
(hfate.GetXaxis())->SetBinLabel(4, "LG Miss");
(hfate.GetXaxis())->SetBinLabel(5, "Exit");
(hfate.GetXaxis())->SetBinLabel(6, "Enter");
(hfate.GetXaxis())->SetBinLabel(7, "Simulated");
(hfate.GetXaxis())->SetBinLabel(8, "Absorb");
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 = TH1F("hnodb_all", "", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
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 = TH1F("hnodb_exit", "", 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;
//hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin;
hin = TH2F("hin", ";x [mm]; y[mm]", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0);
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 = TH2F("hout", ";x [mm];y [mm]", nBins, -a/2.0, +a/2.0, nBins, -a/2.0, +a/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:
681,7 → 679,7
// -2, ce se ne odbije zaradi koncnega R stranic
// -3, ce se odbije nazaj in gre nazaj ven skozi sprednjo ploskev
// +4, ce se absorbira v materialu
Fate Guide::PropagateRay(CRay in, CRay& out, int *n_points, TVector3 *points)
Fate Guide::PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
{
if (dbg) printf("--- GUIDE::PropagateRay ---\n");
// ray0 - incident ray
695,7 → 693,7
int n_odb = 0;
int last_hit = 0;
int propagation = 0;
int result = s_side[0].PropagateRay(ray0, ray1, &vec1);
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");
708,8 → 706,8
} else {
if (dbg) printf(" GUIDE: ray entered\n");
points[0] = ray1.GetR();
hfate.Fill(enter); // enter
hin.Fill(vec1.y(), vec1.z());
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) {
720,7 → 718,7
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);
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;
735,7 → 733,7
} // backreflection
 
// the passage is possible, test propagation
propagation = s_side[inters_i].PropagateRay(ray0, ray1, &vec1);
propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
 
if (dbg) printf(" GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
 
748,13 → 746,13
if(inters_i == 5) {
if (_badCoupling) {
TVector3 hitVector(0,0,0);
bool hitActive = grease.TestIntersection(&hitVector, ray0);
bool hitActive = grease->TestIntersection(&hitVector, ray0);
if (hitActive and dbg) printf(" GUIDE: hit grease\n");
if (!hitActive) propagation = noCoupling.PropagateRay(ray0, ray1, &vec1);
if (!hitActive) propagation = noCoupling->PropagateRay(ray0, &ray1, &vec1);
}
// check on which side the vector is?
TVector3 ray = ray1.GetK();
TVector3 exitNormal = s_side[5].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");
765,8 → 763,8
break;
}
fate = hitExit;
hout.Fill(vec1.y(), vec1.z());
hnodb_exit.Fill(n_odb-1);
hout->Fill(vec1.y(), vec1.z());
hnodb_exit->Fill(n_odb-1);
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
802,26 → 800,25
}
//--- material absorption ---
 
hfate.Fill(fate);
hfate.Fill(rays);
hnodb_all.Fill(n_odb-2);
hfate->Fill(fate);
hfate->Fill(rays);
hnodb_all->Fill(n_odb-2);
*n_points = n_odb+1;
out = ray0;
*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);
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());
902,8 → 899,8
col_glass(4),
col_active(7),
guide_on(parameters.getGuideOn()),
guide(Guide(center0, parameters)),
plate(Plate(parameters)),
guide(new Guide(center0, parameters)),
plate(new Plate(parameters)),
_plateWidth(parameters.getPlateWidth()),
_plateOn(parameters.getPlateOn()),
offsetY(parameters.getOffsetY()),
932,20 → 929,20
// example: epoxy n=1.60
double n4 = 1.57;
TVector3 plane_v[4];
int nBins = nch;
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.Set(SURF_REFRA, plane_v, n3, n4, reflectivity);
glass.FlipN();
glass = new CSurface(SURF_REFRA, plane_v, n3, n4, reflectivity);
glass->FlipN();
 
// additional circular glass between LG and SiPM
glass_circle.Set(TVector3(x_offset+d+glass_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), b);
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 = TH2F("hglass", "",
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);
 
957,22 → 954,22
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.Set(plane_v);
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.Set(activePosition, normal, 1.0*p_size);
grease = new CPlaneR(activePosition, normal, 1.0*p_size);
 
//hactive = (TH2F*)gROOT->FindObject("hactive"); if(hactive) delete hactive;
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 = 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);
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 = TH2F("hlaser", ";x [mm]; y [mm]", nBins, -p_size+offsetY, p_size+offsetY, nBins, -p_size+offsetZ, p_size+offsetZ);
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;
980,11 → 977,11
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.Set(plane_v);
detector = new CPlane4(plane_v);
 
//hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector;
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 = 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);
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);
1000,16 → 997,14
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 = TH2F("histoPlate", ";x [mm];y [mm]", nBins, -p_size, +p_size, nBins, -p_size, +p_size);
 
generated = TH2F("generated", ";x [mm];y [mm]", nBins, -p_size, +p_size, nBins, -p_size, +p_size);
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);
}
 
//-----------------------------------------------------------------------------
// vrne 1 ce je zadel aktvino povrsino
// vrne <1 ce jo zgresi
int CDetector::Propagate(CRay in, CRay& out, int draw)
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
1020,14 → 1015,13
// +4, ce se absorbira v materialu ABSORBED
{
if (dbg) printf("--- Detector::Propagate ---\n");
generated.Fill(in.GetR().y(), in.GetR().z());
//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(in);
rayout.SetColor(col_in);
CRay *rayout = new CRay(in);
rayout->SetColor(col_in);
 
const int max_n_points = guide.GetMAXODB() + 2;
const int max_n_points = guide->GetMAXODB() + 2;
TVector3 pointsPlate[max_n_points];
//TVector3 intersection;
Fate fatePlate;
1041,16 → 1035,16
 
if(_plateOn) {
 
fatePlate = plate.propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
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);
rayout->SetColor(col_in);
rayout->DrawS(center.x() - _plateWidth, -10.0);
}
else if(fatePlate == backreflected){
rayout.SetColor(kBlack);
rayout.DrawS(center.x() - _plateWidth, 7.0);
rayout->SetColor(kBlack);
rayout->DrawS(center.x() - _plateWidth, 7.0);
if (dbg) printf("Backreflected at plate!\n");
}
else {
1060,31 → 1054,31
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);
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);
rayout->SetColor(col_out);
rayout->DrawS(pointsPlate[p_i].x(), 10.0);
}
}
}
 
if(! (fatePlate == hitExit or fatePlate == refracted) ) {
guide.GetHFate()->Fill(rays);
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
guide->GetHFate()->Fill(fatePlate); // reflected back
else
guide.GetHFate()->Fill(noreflection); //lost on plate side
guide->GetHFate()->Fill(noreflection); //lost on plate side
return fatePlate;
}
 
//Ray hits light guide
histoPlate.Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
 
}
else {
//rayout = rayin;
if(draw) rayout.DrawS(center.x(), -10.0);
if(draw) rayout->DrawS(center.x(), -10.0);
}
 
// If the ray is not reflected in the plate
1096,12 → 1090,12
 
int n_points;
int fate_glass;
CRay *ray0 = new CRay(rayout);
CRay *ray0 = new CRay(*rayout);
// delete rayout; -> creates dangling reference when tries to delete ray0!
//delete rayin; -> delete rayout!
CRay ray1;
CRay *ray1 = new CRay;
 
fate = guide.PropagateRay(*ray0, ray1, &n_points, points);
fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
if (dbg) {
if (fate == backreflected) printf("DETECTOR::backreflected\n");
}
1112,10 → 1106,10
if(draw) {
if(fate == missed) {
if (dbg) printf("Detector: fate=missed\n");
TVector3 r = ray1.GetR();
TVector3 k = ray1.GetK();
ray1.Set(r,k);
ray1.DrawS(center.x(), 10.0);
TVector3 r = ray1->GetR();
TVector3 k = ray1->GetK();
ray1->Set(r,k);
ray1->DrawS(center.x(), 10.0);
} else {
for(p_i = 0; p_i < n_points-1; p_i++) {
line3d->SetPoint(0, points[p_i].x(), points[p_i].y(), points[p_i].z());
1124,10 → 1118,10
}
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);}
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);
ray1->SetColor(col_out);
ray1->DrawS(points[p_i].x(), 10.0);
}
}
}
1136,34 → 1130,34
 
if(! (fate == hitExit or fate == refracted) ) {
if (dbg) printf("Detector: fate != hit, refracted\n");
out = ray1;
*out = *ray1;
delete ray0;
//delete ray1;
//delete rayout;
delete ray1;
delete rayout;
delete rayin;
return fate;
}
} else {
if (dbg) printf("Detector: fate = hit or refracted");
ray1 = *ray0;
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);
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);
*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);
hglass->Fill(presecisce.y(), presecisce.z());
if(draw) ray1->DrawS(presecisce.x(), 10.0);
//if(active->TestIntersection(&presecisce, *ray1)) {
//fate = hitExit;
//h-Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
//hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
//hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ);
//}
//if(detector->TestIntersection(&presecisce, *ray1))
1170,28 → 1164,28
//hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
//} else if(fate_glass == REFLECTION) {
else
if(draw) ray1.DrawS(presecisce.x(), 10.0);
if(draw) ray1->DrawS(presecisce.x(), 10.0);
}
}
 
// Main test: ray and SiPM surface
if(active.TestIntersection(&presecisce, ray1)) {
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);
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());
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;
*out = *ray1;
delete ray0;
//delete ray1;
//delete rayout;
delete ray1;
delete rayout;
delete rayin;
return fate;
}
1199,19 → 1193,19
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);
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(_plateOn) plate->draw(4, width);
guide->Draw(col_LG, width);
}
}
 
if(glass_on) glass_circle.Draw(col_glass, width);
if(glass_on) glass_circle->Draw(col_glass, width);
//window_circle->Draw(col_glass, width);
active.Draw(col_active, width);
active->Draw(col_active, width);
}
 
 
1236,23 → 1230,23
 
for(int i = 0; i<8; i++) plate_edge[i] += center;
 
sides[0].Set(SURF_REFRA, plate_edge, n1, n2, c_reflectivity);
sides[0].FlipN();
sides[0] = new CSurface(SURF_REFRA, plate_edge, n1, n2, c_reflectivity);
sides[0]->FlipN();
 
sides[1].Set(SURF_REFRA, plate_edge[3], plate_edge[2], plate_edge[6], plate_edge[7], n2, n2, c_reflectivity);
sides[2].Set(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, c_reflectivity);
sides[3].Set(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, c_reflectivity);
sides[4].Set(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, c_reflectivity);
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].Set(SURF_REFRA, &plate_edge[4], n2, n2, c_reflectivity);
sides[5].FlipN();
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);
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)
1268,7 → 1262,7
}
}
 
Fate Plate::propagateRay(CRay in, CRay& out, int *n_points, TVector3 *points)
Fate Plate::propagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
{
CRay ray0;
CRay ray1;
1281,7 → 1275,7
int last_hit = 0;
int propagation = 0;
 
int result = sides[0].PropagateRay(ray0, ray1, &vec1);
int result = sides[0]->PropagateRay(ray0, &ray1, &vec1);
if( !result ) {
// ce -NI- presecisca z vstopno
fate = missed;
1299,7 → 1293,7
propagation = 11;
for(inters_i=0; inters_i<6; inters_i++) {
if( inters_i != last_hit) {
if( sides[inters_i].TestIntersection(&vec1, ray1) ) break;
if( sides[inters_i]->TestIntersection(&vec1, ray1) ) break;
}
}
points[n_odb] = vec1;
1308,7 → 1302,7
fate = backreflected;
break;} // backreflection
 
propagation = sides[inters_i].PropagateRay(ray0, ray1, &vec1);
propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
if(inters_i == 5) { // successfull exit
fate = hitExit;
//hout->Fill(vec1.y(), vec1.z());
1335,7 → 1329,7
}
 
*n_points = n_odb+1;
out = ray0;
*out = ray0;
return fate;
};