/lightguide/trunk/include/guide.h |
---|
199,14 → 199,14 |
_n2(n2), |
_n3(n3), |
_gap(gap), |
_fresnel(1), |
_guideOn(1), |
offsetY(b/2.0), |
offsetZ(b/2.0), |
_plateOn(1), |
_plateWidth(1), |
_fresnel(1), |
_glassOn(0), |
_glassD(0), |
_guideOn(1), |
_plateOn(1), |
offsetY(2.5), |
offsetZ(2.5) |
_glassD(0) |
{}; |
DetectorParameters(double a, double b, double d): |
_a(a), |
214,40 → 214,35 |
_d(d), |
_active(a), |
_n1(1.0), |
_n2(1.50), |
_n3(1.50), |
_n2(1.53), |
_n3(1.46), |
_gap(TVector3(0,0,0)), |
_plateWidth(1), |
_fresnel(1), |
_glassOn(0), |
_glassD(0), |
_guideOn(1), |
offsetY(b/2.0), |
offsetZ(b/2.0), |
_plateOn(1), |
offsetY(2.5), |
offsetZ(2.5) |
_plateWidth(1), |
_glassOn(0), |
_glassD(0) |
{}; |
~DetectorParameters() {}; |
void setGuide(double a, double b, double d, double n1, double n2, double n3) { |
void setGuide(double a, double b, double d) { |
_a = a; |
_b = b; |
//_M = b/a; |
_d = d; |
_n1 = n1; |
_n2 = n2; |
_n3 = n3; |
}; |
void setGuide(double a, double M, double d) { |
_a = a; |
_b = M*a; |
_d = d; |
}; |
void setGap(double x, double y, double z) { _gap = TVector3(x,y,z); }; |
void setFresnel(int fresnel) { _fresnel = fresnel; }; |
void setGlass(int glassOn, double glassD) { _glassOn = glassOn; _glassD = glassD; }; |
void setGuideOn(int guideOn) { _guideOn = guideOn; }; |
void setPlate(int plateOn, double plateWidth) { _plateOn = plateOn; _plateWidth = plateWidth; }; |
void setIndices(double n1, double n2, double n3) { _n1 = n1; _n2 = n2; _n3 = n3; }; |
double getM() { return _b/_a; }; |
double getLightYield() { return ( Power(_b,2)/Power(_active,2)); }; |
double getM() { return (_b/_a); }; |
double getA() {return _a;}; |
double getB() {return _b;}; |
double getD() {return _d;}; |
274,15 → 269,18 |
double _n2; |
double _n3; |
TVector3 _gap; |
double _plateWidth; |
int _fresnel; |
int _guideOn; |
double offsetY; |
double offsetZ; |
int _fresnel; |
int _plateOn; |
double _plateWidth; |
int _glassOn; |
double _glassD; |
int _guideOn; |
int _plateOn; |
}; |
435,6 → 433,7 |
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;}; |
/lightguide/trunk/include/RTUtil.h |
---|
8,35 → 8,32 |
#include "TStyle.h" |
#include "TDatime.h" |
//########################################################################################## |
// nastavi nastavitve izrisa, klicem z RTSetStyle(gStyle) |
void RTSetStyle(TStyle *style); |
//########################################################################################## |
// poenostavi delo s TCanvas |
class RTCanvas |
{ |
protected: |
TCanvas *can; |
TPad *pad; |
TPaveLabel *title, *date; |
public: |
RTCanvas(); |
RTCanvas(char *w_title, char *c_title, int x_min, int y_min, int x_size, int y_size); |
~RTCanvas() { |
delete can; |
delete pad; |
delete title; |
delete date; |
} |
void Divide(int nx, int ny); |
void Divide(int np); |
TPad* cd(int i); |
void SaveAs(const char *filename); |
void Update(); |
protected: |
TCanvas *can; |
TPad *pad; |
TPaveLabel *title, *date; |
}; |
//########################################################################################## |
#endif |
/lightguide/trunk/src/RTUtil.cpp |
---|
1,6 → 1,9 |
//########################################################################################## |
#include "TGaxis.h" |
#include "TColor.h" |
#include "include/RTUtil.h" |
//########################################################################################## |
void RTSetStyle(TStyle *style) |
{ |
style->SetStatBorderSize(1); |
13,12 → 16,23 |
style->SetStatColor(0); |
style->SetPalette(1, 0); |
style->SetMarkerStyle(kFullDotLarge); |
//style->SetMarkerSize(7); |
style->SetOptStat("ne"); |
style->SetOptFit(1); |
style->SetPadTopMargin(0.10); |
style->SetPadBottomMargin(0.10); |
style->SetPadLeftMargin(0.10); |
style->SetPadRightMargin(0.12); |
style->SetPadBottomMargin(0.12); |
style->SetPadLeftMargin(0.12); |
style->SetPadRightMargin(0.15); |
style->SetTitleOffset(1.3, "y"); |
style->SetTitleOffset(1.5, "y"); |
style->SetPalette(1, 0); |
style->SetPaperSize(TStyle::kA4); |
TGaxis::SetMaxDigits(4); |
} |
//########################################################################################## |
RTCanvas::RTCanvas() |
49,6 → 63,19 |
pad->Divide(nx, ny, 0.003, 0.005); |
} |
//------------------------------------------------------------------------------------------ |
void RTCanvas::Divide(int np) |
{ |
if( np==2 ) pad->Divide(1, 2, 0.003, 0.005); |
else if( 2<np && np<=4 ) pad->Divide(2, 2, 0.003, 0.005); |
else if( 4<np && np<=6 ) pad->Divide(2, 3, 0.003, 0.005); |
else if( 6<np && np<=8 ) pad->Divide(2, 4, 0.003, 0.005); |
else if( np==9 ) pad->Divide(3, 3, 0.003, 0.005); |
else if( 9<np && np<=12) pad->Divide(3, 4, 0.003, 0.005); |
else if(12<np && np<=16) pad->Divide(4, 4, 0.003, 0.005); |
else if(16<np && np<=25) pad->Divide(5, 5, 0.003, 0.005); |
else if(25<np && np<=32) pad->Divide(4, 8, 0.003, 0.005); |
} |
//------------------------------------------------------------------------------------------ |
TPad* RTCanvas::cd(int i) |
{ |
return (TPad*)(pad->cd(i)); |
58,6 → 85,22 |
{ |
can->SaveAs(filename); |
} |
//------------------------------------------------------------------------------------------ |
void RTCanvas::Update() |
{ |
can->Update(); |
} |
void SetGS() |
{ |
const Int_t Number = 2; |
Double_t Red[Number] = {1.0, 0.0}; |
Double_t Green[Number] = {1.0, 0.0}; |
Double_t Blue[Number] = {1.0, 0.0}; |
Double_t Stops[Number] = {0.0, 1.0}; |
Int_t nb = 50; |
TColor::CreateGradientColorTable(Number, Stops, Red, Green, Blue, nb); |
} |
//########################################################################################## |
/lightguide/trunk/src/userFunctions.cpp |
---|
9,29 → 9,38 |
//extern int show_data; |
// Set the simulation parameters |
void Show3D(int b) {show_3d = b;} |
void ShowData(int b) {show_data = b;} |
void showVisual(int b) {show_3d = b;} |
void showData(int b) {show_data = b;} |
// Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap) |
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.52, TVector3(0.3, 0, 0)); |
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0)); |
// Print the detector parameters |
void getParameters(); |
//void SetLGType(int in = 1, int side = 1, int out = 0) |
//{detector->SetLGType(in, side, out);} |
void SetLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n10 = 1.0, double n20 = 1.53, double n30 = 1.52) |
{ parameters.setGuide(SiPM0, b0, d0, n10, n20, n30); } |
void setLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n1 = 1.0, double n2 = 1.53, double n3 = 1.46) |
{ parameters.setGuide(SiPM0, b0, d0); |
parameters.setIndices(n1, n2, n3); } |
void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0) |
void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0) |
{ parameters.setGap(x_gap0, y_gap0, z_gap0); } |
void SetGlass(double glassOn, double glassD) |
void setGlass(double glassOn, double glassD) |
{ parameters.setGlass(glassOn, glassD); }; |
void SetPlate(int plateOn = 1, double plateWidth = 1) |
void setPlate(int plateOn = 1, double plateWidth = 1) |
{ parameters.setPlate(plateOn, plateWidth); } |
void SetFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); } |
void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); } |
// Refractive Indices |
// n1 - around light guide - air |
// n2 - light guide (plate) material - k9 glass |
// n3 - the material at the exit - optical grease, epoxy, air, etc. |
void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);} |
int save_ary = 0; |
//------------------------------------------------------------------------------------------ |
void Help() |
107,11 → 116,16 |
void SetPol(double x, double y, double z) |
{p_pol0.SetXYZ(x,y,z);} |
// Test function |
// creates an optical boundary surface |
// shows the propagation of light ray |
// and polarization state of the incident ray (green) |
// and surface normal (blue) |
void PolTest(double theta = 0.0) |
{ |
int p_type = 1; |
double p_n1 = 1.5; |
double p_n2 = 1.0; |
double p_n1 = parameters.getN1(); |
double p_n2 = parameters.getN2(); |
theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6; |
TVector3 p_pol; |
130,11 → 144,11 |
#define SURF_IMPER 4 |
*/ |
CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN(); |
surf->SetFresnel(0); |
surf->SetFresnel(1); |
surf->Draw(); |
CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta)); |
ray->SetColor(1); |
ray->SetColor(kBlack); |
//p_pol = rotatey(p_pol0, -theta); |
p_pol = p_pol0; p_pol.RotateY(-theta); |
//printf("p_pol = "); printv(p_pol); |
142,18 → 156,21 |
ray->DrawS(cx, -5.0); |
CRay *out = new CRay; out->SetColor(2); |
CRay *out = new CRay; out->SetColor(kRed); |
TVector3 *inters = new TVector3; |
int fate = surf->PropagateRay(*ray, out, inters); |
if(fate == 1) out->DrawS(cx, 5.0); |
surf->PropagateRay(*ray, out, inters); |
//if(fate == 1) out->DrawS(cx, 5.0); |
out->DrawS(cx, 5.0); |
CRay *incidentPolarization = new CRay; |
incidentPolarization->SetColor(kGreen); |
incidentPolarization->Set(ray->GetR(), p_pol); |
incidentPolarization->DrawS(cx, 1.0); |
CRay *pp = new CRay; pp->SetColor(3); |
pp->Set(ray->GetR(), p_pol); |
pp->DrawS(cx, 2.0); |
CRay *pn = new CRay; pn->SetColor(4); |
pn->Set(ray->GetR(), surf->GetN()); |
pn->DrawS(cx, 3.0); |
CRay *surfaceNormal = new CRay; |
surfaceNormal->SetColor(kBlue); |
surfaceNormal->Set(ray->GetR(), surf->GetN()); |
surfaceNormal->DrawS(cx, 1.0); |
} |
void ptt() |
192,7 → 209,7 |
{ |
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
Init(); |
double izkoristek = RandYZ(detector, parameters, NN, theta, phi); |
double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
204,9 → 221,10 |
Init(); |
CDetector *detector = new CDetector(CENTER, parameters); |
//CDetector detector = new CDetector(); |
double izkoristek = RandIso(detector, parameters, NN, theta, nnrays, showr); |
double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); PrintGuideStat(izkoristek); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
//canvasDetector->cd(); |
213,6 → 231,20 |
//detector->Draw(); |
} |
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, int nnrays = 30, int showr = 0) |
{ |
Init(); |
CDetector *detector = new CDetector(CENTER, parameters); |
//CDetector detector = new CDetector(); |
double izkoristek = beamtest(detector, parameters, NN, 18.5, phiMin, phiMax, nnrays, showr); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, 18.5, izkoristek); |
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
//canvasDetector->cd(); |
//detector->Draw(); |
} |
//------------------------------------------------------------------------------------------ |
230,7 → 262,7 |
parameters.setGap(step[i], 0.0, 0.0); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, theta); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
//printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0); |
PrintGuideStat(acc[i]); |
} |
246,13 → 278,45 |
} |
//------------------------------------------------------------------------------------------ |
void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0) |
{ |
//int tmp3d = show_3d, tmpdata = show_data; |
double show_rays = 10; |
double step[steps], acc[steps]; |
//show_3d = 0; show_data = 0; |
//printf(" theta | acceptance \n"); |
PrintGuideHead(); |
for(int i=0; i<=steps; i++) { |
Init(); |
step[i] = min + i*(max - min)/(double)steps; |
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays); |
//printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0); |
PrintGuideStat(acc[i]); |
delete detector; |
} |
//char sbuff[256]; |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z()); |
DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
//show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 5.0) |
//void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0) |
void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15) |
{ |
int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
//int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
show_3d = 0; show_data = 0; |
//show_3d = 0; show_data = 0; |
//printf(" theta | acceptance \n"); |
PrintGuideHead(); |
260,10 → 324,11 |
Init(); |
step[i] = min + i*(max - min)/(double)steps; |
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
acc[i] = RandYZ(detector, parameters, NN, step[i], phi); |
//acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30); |
acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d); |
//printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0); |
PrintGuideStat(acc[i]); |
//delete detector; |
delete detector; |
} |
//char sbuff[256]; |
273,19 → 338,47 |
DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
//show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10) |
{ |
int tmp3d = show_3d, tmpdata = show_data; |
show_3d = 1; show_data = 0; |
//PrintGuideHead(); |
TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi); |
double min = 0.0; |
printf("\nWait, this takes a while "); |
for(int i=0; i<=steps; i++) { |
double theta = min + i*(maxTheta - min) / (double)steps; |
for (int j=0; j<=steps; j++) { |
Init(); |
double phi = min + j*(maxPhi - min) / (double)steps; |
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
double acc = RandYZ(detector, parameters, NN, theta, phi, 30); |
//PrintGuideStat(acc); |
hAcceptance->Fill(i, j, acc); |
//printf("Acc: %f ", acc); |
delete detector; |
} |
printf("."); |
} |
printf("\n"); |
//DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000); |
cp->cd(); |
hAcceptance->Draw("COLZ"); |
} |
// a vs. d |
void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0) |
{ |
int tmp3d = show_3d; |
int tmpdata = show_data; |
//char sbuff[256]; |
show_3d = 0; // don't show simulations |
show_data = 1; |
//show_3d = 0; // don't show simulations |
//show_data = 1; |
//double d = detector->GetD(); |
const double b = parameters.getB(); // upper side of LG |
297,34 → 390,43 |
// Use the Fresnel eq. instead of fixed reflectivity 96% |
//detector->SetFresnel(1); |
printf(" d | a | Acceptance\n"); |
//printf(" d | a | Acceptance\n"); |
getParameters(); |
printf("Wait, this takes a while "); |
for(int i=0; i<steps; i++) { |
const double d = hAcceptance->GetXaxis()->GetBinCenter(i); |
for(int j=0; j<steps; j++) { |
Init(); |
const double a = hAcceptance->GetYaxis()->GetBinCenter(j); |
const double M = b/a; |
parameters.setGuide(a, M, d); |
//const double M = b/a; |
parameters.setGuide(a, b, d); |
CDetector *detector = new CDetector(CENTER, parameters); |
//detector->guide->setLG(y, M, x); |
//Init(); exclude simulation |
double acceptance = RandIso(detector, parameters, NN, theta, 0, 0); |
double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d); |
//double acceptance = Grid(NN, theta); |
//double acceptance = -1.0; |
hAcceptance->Fill(d, a, acceptance); |
//printf("%.2lf | %.2lf | ", x,y); |
//PrintGuideStat(izkoristek); |
//printf("%.2lf | %.2lf | ", d, a); |
//PrintGuideStat(acceptance); |
delete detector; //works fine, 50x50 grid takes ~4MB of RAM |
} |
printf("."); |
} |
printf("\n"); |
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480); |
TVirtualPad *pacc = cp->cd(0); |
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200); |
cp->cd(); |
//TVirtualPad *pacc = cp->cd(0); |
/* |
pacc->SetRightMargin(0.10); |
pacc->SetLeftMargin(0.10); |
pacc->SetTopMargin(0.10); |
pacc->SetBottomMargin(0.10); |
*/ |
TFile *file = new TFile("acceptance.root","RECREATE"); |
hAcceptance->Write(); |
336,8 → 438,10 |
hAcceptance->SetTitle(";d [mm];a [mm]"); |
hAcceptance->GetXaxis()->SetRangeUser(minD,maxD); |
hAcceptance->Draw("COLZ"); |
char filename[128]; |
sprintf(filename,"LGI_ad%d.C", steps); |
cp->SaveAs(filename); |
show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
364,11 → 468,11 |
for(int i=0; i<=steps; i++) { |
step[i] = min + i*(max - min)/(double)steps; |
sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg |
parameters.setGuide(sipm_d, M0/sipm_d, step[i]); |
parameters.setGuide(sipm_d, M0*sipm_d, step[i]); |
//printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]); |
printf("%.1lf | %.2lf | ", step[i], sipm_d); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, theta); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
PrintGuideStat(acc[i]); |
sprintf(sbuff, "d%2d.gif", i); |
//c3dview->SaveAs(sbuff); |
422,12 → 526,13 |
PrintGuideHead(); |
for(int i=0; i<=steps; i++) { |
step[i] = min + i*(max - min)/(double)steps; |
//parameters.setGuide(a, magnif, step[i]); |
parameters.setGuide(a, magnif, step[i]); |
CDetector *detector = new CDetector(CENTER, parameters); |
//printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]); |
printf("%.1lf | ", step[i]); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, theta); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
PrintGuideStat(acc[i]); |
delete detector; |
} |
456,10 → 561,11 |
PrintGuideHead(); |
for(int i=0; i<=steps; i++) { |
step[i] = min + i*(max - min)/(double)steps; |
parameters.setGuide(a, step[i], d); |
//parameters.setGuide(a, step[i], d); |
parameters.setGuide(a, a*step[i], d); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, theta); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
PrintGuideStat(acc[i]); |
} |
472,5 → 578,17 |
show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
void getParameters() |
{ |
printf("LIGHT GUIDE\n" |
" b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA()); |
printf("MATERIAL REFRACITVE INDICES\n" |
" n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3()); |
printf("PLATE\n" |
" ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel()); |
return; |
} |
/lightguide/trunk/src/raySimulator.cpp |
---|
87,8 → 87,8 |
{ |
if(show_data) { |
char sbuff[256]; |
sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf", |
parameters.getA(), parameters.getM(), parameters.getD(), |
sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf", |
parameters.getA(), parameters.getLightYield()*acc, parameters.getD(), |
parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc); |
if(!only2d) { |
183,6 → 183,7 |
theta = Pi()*theta/180.0; |
if(theta < 1e-6) theta = 1e-6; |
phi = phi*Pi()/180.0; |
if(phi < 1e-6) phi = 1e-6; |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
195,6 → 196,9 |
detector->Propagate(*ray0, ray1, show_3d); |
delete ray0; |
delete ray1; |
return (detector->GetHActive())->GetEntries() / (double)(1); |
} |
//----------------------------------------------------------------------------- |
224,16 → 228,26 |
if(i == (NN/2)) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
detector->Propagate(*ray0, ray1, 0); |
delete ray0; |
delete ray1; |
} |
} |
return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN)); |
double acceptance = 0.0; |
/* |
if( !(parameters.getPlateOn()) ) |
acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN; |
else |
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
*/ |
acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN)); |
return acceptance; |
} |
//----------------------------------------------------------------------------- |
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika) |
// vsi pod kotom (theta, phi) |
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30) |
double RandYZ(CDetector *detector, DetectorParameters& parameters, 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; |
272,13 → 286,25 |
//delete ray0; |
//delete ray1; |
} |
return (detector->GetHActive())->GetEntries() / (double)NN; |
double acceptance = 0.0; |
/* |
if( !(parameters.getPlateOn()) ) |
acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
else |
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
//return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/ |
acceptance = (detector->GetHActive())->GetEntries() / (double) NN; |
return acceptance; |
} |
//----------------------------------------------------------------------------- |
// zarki, izotropno porazdeljeni znotraj kota theta |
// = nakljucni vstopni polozaj in kot |
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0) |
// NN - number of rays to be simulated with angles theta distributed uniformly |
// in cos theta and phi uniformly: |
// theta [0, theta] |
// phi [0,360] |
double RandIso(CDetector *detector, DetectorParameters& parameters, 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; |
292,6 → 318,7 |
double rphi, rtheta; |
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
double theta_min_rad = TMath::Cos(theta); |
double theta_max_rad = TMath::Cos(thetaMin); |
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
322,7 → 349,9 |
rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
rphi = rand.Uniform(0.0, 2.0*pi); |
rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0))); |
//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(rtheta); |
rm = TMath::Sin(rtheta)*TMath::Sin(rphi); |
345,8 → 374,8 |
detector->Propagate(*ray0, ray1, 0); |
} |
//delete ray0; |
//delete ray1; |
delete ray0; |
delete ray1; |
} |
if(show_rand) { |
363,8 → 392,113 |
c2rand->cd(6); hn->Draw(); |
} |
double acceptance = 0.0; |
/* |
if( !(parameters.getPlateOn()) ) |
acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
else |
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
*/ |
acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
return acceptance; |
} |
// 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) |
{ |
double pi = 3.14159265358979312; |
theta *= pi/180.0; |
if(theta < MARGIN) theta = MARGIN; |
phiMin *= pi/180.0; |
phiMax *= pi/180.0; |
double acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
TDatime now; |
TRandom rand; |
rand.SetSeed(now.Get()); |
double rx, ry, rz, rl, rm, rn; |
double rphi; |
//double rtheta; |
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2); |
//double theta_min_rad = TMath::Cos(theta); |
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
hphi = new TH1F("hphi", "hphi", 100, -pi, pi); |
htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta; |
htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0); |
hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta; |
hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0); |
hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl; |
hl = new TH1F("hl", "hl", 100, 0.0, 1.0); |
hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm; |
hm = new TH1F("hm", "hm", 100, 0.0, 1.0); |
hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn; |
hn = new TH1F("hn", "hn", 100, 0.0, 1.0); |
//detector->Init(); |
if (show_3d) detector->Draw(draw_width); |
double b = parameters.getB(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
for(int i = 0; i < NN; i++) { |
rx = CENTER.x() - offset; |
ry = rand.Uniform(-b/2.0, b/2.0); |
rz = rand.Uniform(-b/2.0, b/2.0); |
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(rphi); |
rn = TMath::Sin(theta)*TMath::Cos(rphi); |
if(show_rand) { |
//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); |
CRay *ray1 = new CRay; |
if(i < show_rays) { |
detector->Propagate(*ray0, ray1, show_3d); |
} |
else { |
detector->Propagate(*ray0, ray1, 0); |
} |
delete ray0; |
delete ray1; |
} |
if(show_rand) { |
TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand"); |
if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700); |
else c2rand->Clear(); |
c2rand->Divide(2,3); |
c2rand->cd(1); hphi->Draw(); |
c2rand->cd(3); htheta->Draw(); |
c2rand->cd(5); hcostheta->Draw(); |
c2rand->cd(2); hl->Draw(); |
c2rand->cd(4); hm->Draw(); |
c2rand->cd(6); hn->Draw(); |
} |
double acceptance = 0.0; |
/* |
if( !(parameters.getPlateOn()) ) |
acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
else |
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries(); |
*/ |
acceptance = (detector->GetHActive())->GetEntries() / (double)NN; |
return acceptance; |
} |
/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; |
/lightguide/trunk/examples/compile.cxx |
---|
0,0 → 1,54 |
{ |
TVector3 dum; |
gROOT->ProcessLine(".L src/RTUtil.cpp+"); |
gROOT->ProcessLine(".L src/guide.cpp+"); |
gROOT->ProcessLine(".L src/raySimulator.cpp+"); |
gROOT->ProcessLine(".L src/userFunctions.cpp+"); |
// Draw the detector and propagate certain number of rays |
gROOT->ProcessLine("showVisual(1)"); |
// Draw the statistics |
gROOT->ProcessLine("showData(1)"); |
/* |
// Set the parameters of the light guide: |
// SiPM0: light guide exit surface linear dimension (surface = SiPM0^2) |
// M0: Magnification ratio = b / SiPM |
// b: Light guide entry surface linear dimension = M0*SiPM (surface = (M0*SiPM0)^2) |
// d0: length of light array |
// n10: refractive index outside of the light guide - entry surface and lateral sides |
// n20: refractive index of light guide material and the glass plate |
// n30: refractive index outside of the light guide at the exit surface (ie optical cement, detector) - it is used only by Fresenl eq, not for the refraction |
//gROOT->ProcessLine("SetLG(3.0, 1.666, 3.0, 1.0, 1.48, 1.48)"); |
// Set light guide surface types |
// |
// void SetLGType(int in = 1, int side = 1, int out = 1) |
// 0 = dummy; 1 = refractor; 2 = reflector; 3 = total reflection (oboslete, total reflection now included with refractor type) |
// |
//gROOT->ProcessLine("SetLGType(1, 1, 1)"); |
// Gap between light guide exit surface and sensitive area (SiPM) |
// supposed it is epoxy of the same refractive index as light guide |
gROOT->ProcessLine("SetGap(0.3, 0, 0)"); |
// Glass surface of the detector on/off; width |
gROOT->ProcessLine("SetGlass(0, 0)"); |
// Glass plate above the light guide on/off |
gROOT->ProcessLine("SetPlate(1, 1)"); |
// Draw the detector and propagate certain number of rays |
gROOT->ProcessLine("Show3D(1)"); |
// Draw the statistics |
gROOT->ProcessLine("ShowData(1)"); |
// Propagate 1e4 rays, isotropicaly distributed in 0-30 degree cone (relative to entry surface normal) |
// draws first 20 rays (slow) |
gROOT->ProcessLine("LGI(10000,30,20)"); |
//gROOT->ProcessLine("LGI_ad()"); |
//gROOT->ProcessLine("LGS()"); |
*/ |
} |