Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 69 → Rev 70

/lightguide/trunk/src/guide.cpp
247,9 → 247,9
 
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_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);
}
 
397,23 → 397,32
// --------------- 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;
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;
double a_tm = 0.0;
TVector3 v_te; // polarization perpendicular to the plane of incidence
TVector3 v_tm; // inbound polarization parallel with the plane of incidence
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) {
v_te = n.Cross(in.GetN()); v_te = v_te.Unit();
v_tm = -v_te.Cross(in.GetN()); v_tm = v_tm.Unit();
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);
428,7 → 437,7
}
// ----------------------------------------------------------------------------
if(type == 3) type = SURF_REFRA; //SURF_TOTAL -> SURF_REFRA
if(type == SURF_TOTAL) type = SURF_REFRA;
switch(type){
// ----------------------------------------------------------------------------
// --------------- refraction from n1 to n2 -----------------------------------
452,9 → 461,8
p_ref = rand.Uniform(0.0, 1.0);
if (dbg) printf(" reflection probability = %f\n", p_ref);
// Total reflection
/*
if( (cosTi <= cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
// If n1>n2 and theta>thetaCritical, total reflection
if( (cosTi < cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
if(dbg) printf(" TOTAL\n");
transmit = in.GetN() + sign_n*2*cosTi*n;
467,7 → 475,7
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
} else { */
} 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));
488,6 → 496,7
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;
496,13 → 505,14
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");
573,7 → 583,8
{
double t;
 
TDatime now; rand.SetSeed(now.Get());
TDatime now;
rand.SetSeed(now.Get());
 
center = center0;
double b = parameters.getB();
582,33 → 593,41
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() ? parameters.getN3(): n1);
double n0 = (parameters.getPlateOn() ? n2 : n1);
n3 = parameters.getN3();
_r = c_reflectivity;
int fresnel = parameters.getFresnel();
 
// 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);
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);
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);
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");
670,7 → 689,7
} else {
if (dbg) printf(" GUIDE: ray entered\n");
points[0] = ray1.GetR();
hfate->Fill(2); // enter
hfate->Fill(enter); // enter
hin->Fill(vec1.y(), vec1.z());
if (dbg) printf(" GUIDE: n_odb = %d\n", n_odb);
890,14 → 909,17
//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; // for faster simulation, not using Fresnel eqs.
//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.60;
TVector3 plane_v[4];
int nBins = nch + 1;
double p_size = b/2.0;
905,8 → 927,10
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, n2, n1, reflectivity); glass->FlipN();
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;
914,7 → 938,7
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;
934,7 → 958,7
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);
982,7 → 1006,8
//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;
CRay *rayout = new CRay(in);
rayout->SetColor(col_in);
 
const int max_n_points = guide->GetMAXODB() + 2;
TVector3 pointsPlate[max_n_points];
999,14 → 1024,16
if(_plateOn) {
fatePlate = plate->propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
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 {
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());
1013,24 → 1040,30
line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z());
line3d->DrawClone();
}
if(fatePlate != noreflection) {
//rayout->SetColor(7);
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);
}
}
}
}
}
if(! (fatePlate == hitExit or fatePlate == refracted) ) {
if (dbg) printf("CDetector::propagate Simulated ray missed the entry surface!\n");
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;
}
// missing: if refracted at plate sides
//if (fatePlate == refracted) return fatePlate;
//Ray hits light guide
histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
 
}
else {
rayout = rayin;
//rayout = rayin;
if(draw) rayout->DrawS(center.x(), -10.0);
}
1045,7 → 1078,7
int fate_glass;
CRay *ray0 = new CRay(*rayout);
// delete rayout; -> creates dangling reference when tries to delete ray0!
delete rayin;
//delete rayin; -> delete rayout!
CRay *ray1 = new CRay;
fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
1084,6 → 1117,10
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 {
1095,18 → 1132,6
}
}
/*
TVector3 pres_wind;
fate = window_circle->TestIntersection(&pres_wind, *ray1);
if(fate == 1) {
hwindow->Fill(pres_wind.y(), pres_wind.z());
if(!guide_on) {
window->PropagateRay(*ray0, ray1, &presecisce);
if(draw) ray1->Draw(center.x() + window_d, center.x() + glass_d);
*ray0 = *ray1;
}
*/
fate = missed; // zgresil aktivno povrsino
if(glass_on) {
*ray0 = *ray1; ray1->SetColor(col_rgla);
1143,6 → 1168,7
delete ray0;
delete ray1;
delete rayout;
delete rayin;
return fate;
}
//-----------------------------------------------------------------------------
1171,27 → 1197,30
const double b = parameters.getB();
const double n1 = parameters.getN1();
const double n2 = parameters.getN2();
const double n3 = parameters.getN3();
const double reflectivity = c_reflectivity;
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);
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, reflectivity);
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, reflectivity);
sides[2] = new CSurface(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, reflectivity);
sides[3] = new CSurface(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, reflectivity);
sides[4] = new CSurface(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, 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] = new CSurface(SURF_REFRA, &plate_edge[4], n2, n3, 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);
1237,7 → 1266,7
fate = backreflected;
} else {
points[0] = ray1.GetR();
//hfate->Fill(2);
//hfate->Fill(enter);
//hin->Fill(vec1.y(), vec1.z());
while (n_odb++ < MAX_REFLECTIONS) {
ray0 = ray1;
1264,7 → 1293,7
break;
}
if(propagation == 1) {
fate = refracted; //at side?
fate = noreflection; //at side
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;