| 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; |
| } |
| } |