2,7 → 2,6 |
|
#include <iostream> |
|
|
// vector output shortcut |
void printv(TVector3 v) |
{ |
448,9 → 447,13 |
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 |
p_ref = rand.Uniform(0.0, 1.0); |
if (dbg) printf(" reflection probability = %f\n", p_ref); |
|
p_ref = rand.Uniform(0.0, 1.0); |
|
// 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; |
463,8 → 466,9 |
|
pol_t = -in.GetP() + sign_n*2*cosTi*n; |
out->SetPolarization(pol_t); |
return 2; |
} else { // ni tot. odboja |
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))); |
478,7 → 482,7 |
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) { |
//if(fresnel) { |
r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse |
r_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel |
|
497,23 → 501,26 |
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); |
} |
|
//p_ref = rand.Uniform(0.0, 1.0); |
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 1; |
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 2; |
} |
} |
return REFLECTION; |
} |
|
//} |
break; |
|
// ---------------------------------------------------------------------------- |
// --------------- reflection at "reflection" probability --------------------- |
// ---------------------------------------------------------------------------- |
523,11 → 530,11 |
cosTi = in.GetN() * n; |
transmit = in.GetN() - 2*cosTi*n; |
out->Set(intersect, transmit); |
return 2; //sdhfvjhsdbfjhsdbcvjhsb |
return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb |
} else { // se ne odbije |
transmit = in.GetN(); |
out->Set(intersect, transmit); |
return -2; |
return ABSORBED; |
} |
break; |
|
542,12 → 549,12 |
} else { // ni tot. odboja |
transmit = in.GetN(); |
out->Set(intersect, transmit); |
return -2; |
return ABSORBED; |
} |
} else { // se ne odbije |
transmit = in.GetN(); |
out->Set(intersect, transmit); |
return -2; |
return ABSORBED; |
} |
break; |
|
556,7 → 563,7 |
break; |
} |
|
return 1; |
return REFRACTION; |
} |
//================================================================================= |
|
605,8 → 612,8 |
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, "Ref Loss"); |
(hfate->GetXaxis())->SetBinLabel(3, "No Tot Refl"); |
(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"); |
614,10 → 621,10 |
(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", VODNIK_MAX_ODB, -0.5, VODNIK_MAX_ODB-0.5); |
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", VODNIK_MAX_ODB, -0.5, VODNIK_MAX_ODB-0.5); |
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; |
644,7 → 651,6 |
CRay ray0; |
CRay ray1; |
TVector3 vec0, vec1; |
//int fate = 11; |
int inters_i = 0; |
|
ray0 = in; |
651,11 → 657,16 |
int n_odb = 0; |
int last_hit = 0; |
int propagation = 0; |
if( !(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"); |
fate = missed; |
hfate->Fill(0); |
//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(); |
662,7 → 673,8 |
hfate->Fill(2); // enter |
hin->Fill(vec1.y(), vec1.z()); |
if (dbg) printf(" GUIDE: n_odb = %d\n", n_odb); |
while (n_odb++ < VODNIK_MAX_ODB) { |
|
while (n_odb++ < MAX_REFLECTIONS) { |
if (dbg) printf(" GUIDE: Boundary test: %d\n",n_odb); |
ray0 = ray1; |
vec0 = vec1; |
680,12 → 692,27 |
points[n_odb] = vec1; |
if(inters_i == 0) { |
fate = backreflected; |
hfate->Fill(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(); |
701,7 → 728,7 |
ray0 = ray1; |
break; |
} |
fate = hit; |
fate = hitExit; |
hout->Fill(vec1.y(), vec1.z()); |
hnodb_exit->Fill(n_odb-1); |
n_odb++; |
709,18 → 736,6 |
ray0 = ray1; |
break; |
} |
if(propagation == REFRACTION) { |
fate = refracted; |
n_odb++; |
points[n_odb] = vec1; |
ray0 = ray1; |
break; |
} // no total reflection when should be |
if(propagation == -2) { |
fate = noreflection; |
break; |
} //refraction due to finite reflectivity |
|
last_hit = inters_i; |
} |
} |
743,7 → 758,7 |
//--- material absorption --- |
|
hfate->Fill(fate); |
hfate->Fill(3); |
hfate->Fill(rays); |
hnodb_all->Fill(n_odb-2); |
*n_points = n_odb+1; |
*out = ray0; |
858,7 → 873,9 |
guide(new Guide(center0, parameters)), |
plate(new Plate(parameters)), |
_plateWidth(parameters.getPlateWidth()), |
_plateOn(parameters.getPlateOn()) |
_plateOn(parameters.getPlateOn()), |
offsetY(parameters.getOffsetY()), |
offsetZ(parameters.getOffsetZ()) |
{ |
// }; |
|
908,17 → 925,17 |
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, 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", "Hits at LG entrance", nBins, -p_size, p_size, nBins, -p_size, p_size); |
hlaser = new TH2F("hlaser", ";x [mm]; y [mm]", nBins, -p_size+offsetY, p_size+offsetY, nBins, -p_size+offsetZ, p_size+offsetZ); |
|
|
p_size = 2*b/2.0; |
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); |
926,9 → 943,9 |
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", "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); |
|
946,10 → 963,12 |
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) |
|
// Sledi zarku skozi vodnik. Vrne: |
// 0, ce zgresi vstopno ploskev MISSED |
// 1, ce zadane izstopno ploskev HIT |
968,7 → 987,7 |
const int max_n_points = guide->GetMAXODB() + 2; |
TVector3 pointsPlate[max_n_points]; |
//TVector3 intersection; |
int fatePlate; |
Fate fatePlate; |
int nPointsPlate; |
TPolyLine3D *line3d = new TPolyLine3D(2); |
line3d->SetLineWidth(1); |
983,7 → 1002,7 |
if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0); |
|
if(draw) { |
if(fatePlate == 0) { |
if(fatePlate == missed) { |
rayout->SetColor(col_in); |
rayout->DrawS(center.x() - _plateWidth, -10.0); |
} |
994,7 → 1013,7 |
line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z()); |
line3d->DrawClone(); |
} |
if(fatePlate != -2) { |
if(fatePlate != noreflection) { |
//rayout->SetColor(7); |
rayout->DrawS(pointsPlate[p_i].x(), -0.1); |
} |
1001,11 → 1020,13 |
} |
} |
|
if(! (fatePlate == 1 or fatePlate == -1) ) { |
if(! (fatePlate == hitExit or fatePlate == refracted) ) { |
|
cout<<"CDetector::propagate Simulated ray missed the entry surface!"<<endl; |
if (dbg) printf("CDetector::propagate Simulated ray missed the entry surface!\n"); |
return fatePlate; |
} |
// missing: if refracted at plate sides |
//if (fatePlate == refracted) return fatePlate; |
histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point |
} |
else { |
1019,15 → 1040,18 |
//const int max_n_points = guide->GetMAXODB() + 2; |
TVector3 points[max_n_points]; |
TVector3 presecisce; |
//int fate; |
|
int n_points; |
int fate_glass; |
CRay *ray0 = new CRay(*rayout); |
// delete rayout; -> creates dangling reference when tries to delete ray0! |
//delete rayin; |
delete rayin; |
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; |
1057,8 → 1081,8 |
} |
|
|
if(! (fate == hit or fate == refracted) ) { |
if (dbg) printf("Detector: fate != hit, refracted\n"); |
if(! (fate == hitExit or fate == refracted) ) { |
if (dbg) printf("Detector: fate != hit, refracted\n"); |
*out = *ray1; |
return fate; |
} |
1091,8 → 1115,8 |
hglass->Fill(presecisce.y(), presecisce.z()); |
if(draw) ray1->DrawS(presecisce.x(), 10.0); |
if(active->TestIntersection(&presecisce, *ray1)) { |
fate = hit; |
hactive->Fill(presecisce.y(), presecisce.z()); |
fate = hitExit; |
hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
hlaser->Fill((in.GetR()).y(), (in.GetR()).z()); |
} |
if(detector->TestIntersection(&presecisce, *ray1)) |
1101,13 → 1125,15 |
if(draw) ray1->DrawS(presecisce.x(), 10.0); |
} |
} else { |
// Main test: ray and SiPM surface |
if(active->TestIntersection(&presecisce, *ray1)) { |
fate = hit; |
hactive->Fill(presecisce.y(), presecisce.z()); |
hlaser->Fill((in.GetR()).y(), (in.GetR()).z()); |
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(presecisce.y(), presecisce.z()); |
hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z()); |
} |
//} else { |
//if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d); |
1114,8 → 1140,9 |
//} |
|
*out = *ray1; |
//delete ray0; |
//delete ray1; |
delete ray0; |
delete ray1; |
delete rayout; |
return fate; |
} |
//----------------------------------------------------------------------------- |
1188,12 → 1215,12 |
} |
} |
|
int 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; |
TVector3 vec0, vec1; |
int fate = 11; |
Fate fate = enter; |
int inters_i = 0; |
|
ray0 = in; |
1201,14 → 1228,18 |
int last_hit = 0; |
int propagation = 0; |
|
if( !(sides[0]->PropagateRay(ray0, &ray1, &vec1)) ) { |
int result = sides[0]->PropagateRay(ray0, &ray1, &vec1); |
if( !result ) { |
// ce -NI- presecisca z vstopno |
fate = 0; |
fate = missed; |
} else if(result == REFLECTION) { |
if (dbg) printf("PLATE: reflected\n"); |
fate = backreflected; |
} else { |
points[0] = ray1.GetR(); |
//hfate->Fill(2); |
//hin->Fill(vec1.y(), vec1.z()); |
while (n_odb++ < VODNIK_MAX_ODB) { |
while (n_odb++ < MAX_REFLECTIONS) { |
ray0 = ray1; |
vec0 = vec1; |
propagation = 11; |
1219,12 → 1250,12 |
} |
points[n_odb] = vec1; |
if(inters_i == 0) { |
fate = -3; |
fate = backreflected; |
break;} // backreflection |
|
propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1); |
if(inters_i == 5) { // successfull exit |
fate = 1; |
fate = hitExit; |
//hout->Fill(vec1.y(), vec1.z()); |
//hnodb_exit->Fill(n_odb-1); |
n_odb++; |
1233,13 → 1264,17 |
break; |
} |
if(propagation == 1) { |
fate = -1; |
fate = refracted; //at side? |
n_odb++; |
points[n_odb] = vec1; |
ray0 = ray1; |
break;} // no total reflection when should be |
if(propagation == -2) {fate = -2; break;} // absorption due to finite reflectivity |
|
if(propagation == -2) { |
fate = noreflection; |
break; |
} // absorption due to finite reflectivity |
|
last_hit = inters_i; |
} |
} |