Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 83 → Rev 84

/lightguide/trunk/src/guide.cpp
32,14 → 32,15
//n.SetXYZ(l0, m0, n0); n = n.Unit();
//}
//-----------------------------------------------------------------------------
/*
CRay& CRay::operator = (const CRay& p)
 
CRay& CRay::operator = (const CRay& rhs)
{
r.SetXYZ(p.GetR().x(), p.GetR().y(), p.GetR().z());
r.SetXYZ(rhs.GetR().x(), rhs.GetR().y(), rhs.GetR().z());
//this->r.SetXYZ(p.x(), p.y(), p.z());
n.SetXYZ(p.GetN().x(), p.GetN().y(), p.GetN().z());
k.SetXYZ(rhs.GetK().x(), rhs.GetK().y(), rhs.GetK().z());
p.SetXYZ(rhs.GetP().x(), rhs.GetP().y(), rhs.GetP().z());
return *this;
} */
}
void CRay::Print()
{
printf("---> CRay::Print() <---\n");
364,7 → 365,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
415,7 → 416,8
}
 
double cosAf = v_te * in.GetP();
double alpha = acos(cosAf);
double alpha = TMath::ACos(cosAf);
if(alpha < MARGIN) alpha = MARGIN;
if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE);
 
a_te = cosAf;
459,7 → 461,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");
470,7 → 472,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",
477,7 → 479,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 {
521,8 → 523,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);
529,10 → 531,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;
}
 
547,11 → 549,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;
563,21 → 565,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;
}
 
620,55 → 622,55
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[0].Set(SURF_REFRA, vodnik_edge, n0, _n2, _r);
s_side[0].FlipN();
 
s_side[1] = new CSurface(SURF_REFRA, vodnik_edge[3], vodnik_edge[2],
s_side[1].Set(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],
s_side[2].Set(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],
s_side[3].Set(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],
s_side[4].Set(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();
s_side[5].Set(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();
noCoupling.Set(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, 0.95*a/2.0);
grease.Set(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 = 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");
//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");
 
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_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_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);
//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);
 
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);
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);
 
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);
//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);
}
//-----------------------------------------------------------------------------
// Sledi zarku skozi vodnik. Vrne:
679,7 → 681,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
693,7 → 695,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");
706,8 → 708,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) {
718,7 → 720,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;
733,7 → 735,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);
 
746,13 → 748,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");
763,8 → 765,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;
800,25 → 802,26
}
//--- 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());
899,8 → 902,8
col_glass(4),
col_active(7),
guide_on(parameters.getGuideOn()),
guide(new Guide(center0, parameters)),
plate(new Plate(parameters)),
guide(Guide(center0, parameters)),
plate(Plate(parameters)),
_plateWidth(parameters.getPlateWidth()),
_plateOn(parameters.getPlateOn()),
offsetY(parameters.getOffsetY()),
929,20 → 932,20
// example: epoxy n=1.60
double n4 = 1.57;
TVector3 plane_v[4];
int nBins = nch + 1;
int nBins = nch;
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();
glass.Set(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);
glass_circle.Set(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", "",
//hglass = (TH2F*)gROOT->FindObject("hglass"); if(hglass) delete hglass;
hglass = TH2F("hglass", "",
nBins, y_gap - p_size, y_gap + p_size,
nBins, z_gap - p_size, z_gap + p_size);
 
954,22 → 957,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 = new CPlane4(plane_v);
active.Set(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);
grease.Set(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 = 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);
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);
 
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);
//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);
 
// collection surface in SiPM plane
p_size = 1.4*b/2.0;
977,11 → 980,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 = new CPlane4(plane_v);
detector.Set(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 = 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);
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);
 
/*
window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R);
997,14 → 1000,16
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);
//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);
}
 
//-----------------------------------------------------------------------------
// 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
1015,13 → 1020,14
// +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 = new CRay(in);
rayout->SetColor(col_in);
CRay rayout(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;
1035,16 → 1041,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 {
1054,31 → 1060,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
1090,12 → 1096,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 = new CRay;
CRay ray1;
 
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");
}
1106,10 → 1112,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());
1118,10 → 1124,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);
}
}
}
1130,34 → 1136,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;
//hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
//h-Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
//hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ);
//}
//if(detector->TestIntersection(&presecisce, *ray1))
1164,28 → 1170,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;
}
1193,19 → 1199,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);
}
 
 
1230,23 → 1236,23
 
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[0].Set(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[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[5] = new CSurface(SURF_REFRA, &plate_edge[4], n2, n2, c_reflectivity);
sides[5]->FlipN();
sides[5].Set(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)
1262,7 → 1268,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;
1275,7 → 1281,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;
1293,7 → 1299,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;
1302,7 → 1308,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());
1329,7 → 1335,7
}
 
*n_points = n_odb+1;
*out = ray0;
out = ray0;
return fate;
};