/lightguide/trunk/include/guide.h |
---|
39,14 → 39,14 |
public: |
CRay() : |
r(TVector3(0,0,0)), |
n(TVector3(1,0,0)), |
k(TVector3(1,0,0)), |
//p(TVector3(0,1,0)), |
p(TVector3(0,0,1)), |
color(kBlack) |
{}; |
CRay(TVector3 r0, TVector3 n0) : |
CRay(TVector3 r0, TVector3 k0) : |
r(r0), |
n(n0.Unit()), |
k(k0.Unit()), |
//p(TVector3(0,1,0)), |
p(TVector3(0,0,1)), |
color(kBlack) |
53,21 → 53,24 |
{}; |
CRay(double x0, double y0, double z0, double l0, double m0, double n0) : |
r(TVector3(x0,y0,z0)), |
n(TVector3(l0,m0,n0).Unit()), |
k(TVector3(l0,m0,n0).Unit()), |
//p(TVector3(0,1,0)), |
p(TVector3(0,0,1)), |
color(kBlack) |
{}; |
void Set(TVector3 r0, TVector3 n0); |
void Set(TVector3 r0, TVector3 k0); |
//void Set(double x0, double y0, double z0, double l0, double m0, double n0); |
void SetColor(int c){color = c;}; |
void SetPolarization(TVector3 p0) {p = p0.Unit();}; |
void setPolarization(TVector3 p0) { |
p = p0.Unit(); |
if (p.Dot(k)>1e-3) printf("*** ERROR in CRay: E has component || with k\n"); |
}; |
//inline CRay & operator = (const CRay &); |
TVector3 GetR() const {return r;}; |
TVector3 GetN() const {return n;}; |
TVector3 GetK() const {return k;}; |
TVector3 GetP() const {return p;}; |
void Print(); |
75,13 → 78,13 |
void Draw(double x_from, double x_to); |
void DrawS(double x_from, double t); |
//r = position, k = unit wave vector, p = polarization |
private: |
TVector3 r; |
TVector3 n; |
TVector3 p; //r = point on line, n = normal, p = polarization |
TVector3 k; |
TVector3 p; |
int color; |
}; |
//================================================================================= |
//================================================================================= |
141,7 → 144,6 |
TVector3 n, center; |
double _r; |
}; |
//================================================================================= |
//================================================================================= |
// ravna opticna povrsina: refractor, zrcalo ali povrsina s totalnim odbojem |
188,7 → 190,6 |
int fresnel; // ali naj uposteva Fresnelove enacbe |
}; |
//================================================================================= |
//================================================================================= |
363,9 → 364,7 |
TH1F *hfate, *hnodb_all, *hnodb_exit; |
TH2F *hin, *hout; |
}; |
//================================================================================= |
//=============================================================================================================================== <<<<<<<< |
class Plate |
385,7 → 384,6 |
CSurface *sides[6]; |
}; |
// ================================================================================ |
class CDetector |
{ |
511,5 → 509,4 |
//================================================================================= |
#endif |
/lightguide/trunk/src/userFunctions.cpp |
---|
18,21 → 18,21 |
void getParameters(); |
//void SetLGType(int in = 1, int side = 1, int out = 0) |
//{detector->SetLGType(in, side, out);} |
//{detector->SetLGType(in, side, out);} |
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); } |
{ 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) |
{ parameters.setGap(x_gap0, y_gap0, z_gap0); } |
{ parameters.setGap(x_gap0, y_gap0, z_gap0); } |
void setGlass(double glassOn, double glassD) |
{ parameters.setGlass(glassOn, glassD); }; |
{ parameters.setGlass(glassOn, glassD); }; |
void setPlate(int plateOn = 1, double plateWidth = 1) |
{ parameters.setPlate(plateOn, plateWidth); } |
{ parameters.setPlate(plateOn, plateWidth); } |
void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); } |
// Refractive Indices |
45,77 → 45,72 |
//------------------------------------------------------------------------------------------ |
void Help() |
{ |
printf("void SetCenter(x, y, z)\n"); |
printf("void SetLGType(in, side, out)\n"); |
printf("SURF_DUMMY 0, "); |
printf("SURF_REFRA 1, "); |
printf("SURF_REFLE 2, "); |
printf("SURF_TOTAL 3, "); |
printf("SURF_IMPER 4\n"); |
printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n"); |
printf("void SetR(R0)\n"); |
printf("void SetGlass(glass_on0, glass_d0)\n"); |
printf("void SetGap(x_gap0, y_gap0, z_gap0)\n"); |
printf("LGG(NN, theta)\n"); |
printf("LGR(NN, theta, phi)\n"); |
printf("LGI(NN, theta)\n"); |
printf("LGI_gap(NN, min, max, steps, theta)\n"); |
printf("LGR_gap(NN, min, max, steps, theta, phi)\n"); |
printf("LGR_th(NN, min, max, steps, phi)\n"); |
printf("LGG_th(NN, min, max, steps)\n"); |
printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)"); |
printf("void SetCenter(x, y, z)\n"); |
printf("void SetLGType(in, side, out)\n"); |
printf("SURF_DUMMY 0, "); |
printf("SURF_REFRA 1, "); |
printf("SURF_REFLE 2, "); |
printf("SURF_TOTAL 3, "); |
printf("SURF_IMPER 4\n"); |
printf("void SetLG(SiPM0, M0, d0, n10, n20, n30, R0)\n"); |
printf("void SetR(R0)\n"); |
printf("void SetGlass(glass_on0, glass_d0)\n"); |
printf("void SetGap(x_gap0, y_gap0, z_gap0)\n"); |
printf("LGG(NN, theta)\n"); |
printf("LGR(NN, theta, phi)\n"); |
printf("LGI(NN, theta)\n"); |
printf("LGI_gap(NN, min, max, steps, theta)\n"); |
printf("LGR_gap(NN, min, max, steps, theta, phi)\n"); |
printf("LGR_th(NN, min, max, steps, phi)\n"); |
printf("LGG_th(NN, min, max, steps)\n"); |
printf("LGI_ad(NN, mina, maxa, mind, maxd, steps(memory dependent), theta)"); |
} |
//------------------------------------------------------------------------------------------ |
int show_in_steps = 0; |
void SetShowInSteps(int in = 1) |
{ |
show_in_steps = in; |
show_in_steps = in; |
} |
//------------------------------------------------------------------------------------------ |
int gs_set = 0; |
void SetGS() |
{ |
const UInt_t Number = 2; |
Double_t Red[Number] = { 1.00, 0.00}; |
Double_t Green[Number] = { 1.00, 0.00}; |
Double_t Blue[Number] = { 1.00, 0.00}; |
Double_t Stops[Number] = { 0.00, 1.00}; |
Int_t nb=50; |
if(!gs_set) { |
TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb); |
gs_set = 1; |
} |
const UInt_t Number = 2; |
Double_t Red[Number] = { 1.00, 0.00}; |
Double_t Green[Number] = { 1.00, 0.00}; |
Double_t Blue[Number] = { 1.00, 0.00}; |
Double_t Stops[Number] = { 0.00, 1.00}; |
Int_t nb=50; |
if(!gs_set) { |
TColor::CreateGradientColorTable(Number,Stops,Red,Green,Blue,nb); |
gs_set = 1; |
} |
} |
//------------------------------------------------------------------------------------------ |
void Draw1(CDetector *detector, int gs = 0, double range = 0.0) |
{ |
if(gs) SetGS(); |
TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500); |
gStyle->SetOptTitle(0); |
gStyle->SetOptStat(0); |
cdata->cd(0); |
if(range > 0.1) { |
((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range); |
((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range); |
} |
(detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]"); |
(detector->GetHLaser())->Draw("COLZ"); |
cdata->SaveAs("2dscan.gif"); |
//gStyle->SetOptTitle(1); |
//gStyle->SetOptStat(1); |
if(gs) SetGS(); |
TCanvas *cdata = new TCanvas("cdata", "2dscan", 350, 0, 500, 500); |
gStyle->SetOptTitle(0); |
gStyle->SetOptStat(0); |
cdata->cd(0); |
if(range > 0.1) { |
((detector->GetHLaser())->GetXaxis())->SetRangeUser(-range,range); |
((detector->GetHLaser())->GetYaxis())->SetRangeUser(-range,range); |
} |
(detector->GetHLaser())->SetTitle("Hits at LG entrance; x [mm]; y [mm]"); |
(detector->GetHLaser())->Draw("COLZ"); |
cdata->SaveAs("2dscan.gif"); |
//gStyle->SetOptTitle(1); |
//gStyle->SetOptStat(1); |
} |
//------------------------------------------------------------------------------------------ |
TVector3 p_pol0(0.0, 0.0, 1.0); |
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 |
123,62 → 118,66 |
// and surface normal (blue) |
void PolTest(double theta = 0.0) |
{ |
int p_type = 1; |
int p_type = SURF_REFRA; |
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; |
Init(); |
double cx = 0.0; |
TVector3 vodnik_edge[4]; |
double t = 3.0; |
vodnik_edge[0].SetXYZ(cx, t,-t); vodnik_edge[1].SetXYZ(cx, t, t); |
vodnik_edge[2].SetXYZ(cx,-t, t); vodnik_edge[3].SetXYZ(cx,-t,-t); |
/* |
#define SURF_DUMMY 0 |
#define SURF_REFRA 1 |
#define SURF_REFLE 2 |
#define SURF_TOTAL 3 |
#define SURF_IMPER 4 |
*/ |
CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN(); |
surf->SetFresnel(1); |
surf->Draw(); |
CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta)); |
ray->SetColor(kBlack); |
//p_pol = rotatey(p_pol0, -theta); |
p_pol = p_pol0; p_pol.RotateY(-theta); |
//printf("p_pol = "); printv(p_pol); |
ray->SetPolarization(p_pol); |
ray->DrawS(cx, -5.0); |
CRay *out = new CRay; out->SetColor(kRed); |
TVector3 *inters = new TVector3; |
surf->PropagateRay(*ray, out, inters); |
//if(fate == 1) out->DrawS(cx, 5.0); |
theta = 3.141593*theta/180.0; |
if(theta < 1e-6) theta = 1e-6; |
Init(); |
double cx = 0.0; |
TVector3 vodnik_edge[4]; |
double t = 3.0; |
vodnik_edge[0].SetXYZ(cx, t,-t); |
vodnik_edge[1].SetXYZ(cx, t, t); |
vodnik_edge[2].SetXYZ(cx,-t, t); |
vodnik_edge[3].SetXYZ(cx,-t,-t); |
CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); |
surf->FlipN(); |
surf->SetFresnel(1); |
surf->Draw(); |
CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta)); |
ray->SetColor(kRed); |
//p_pol = rotatey(p_pol0, -theta); |
TVector3 sE(0,1,0); |
TVector3 pE(0,0,1); |
TVector3 p_pol = pE; |
p_pol.RotateY(-theta); |
//printf("p_pol = "); printv(p_pol); |
ray->setPolarization(p_pol); |
ray->DrawS(cx, -5.0); |
CRay *out = new CRay; |
out->SetColor(kBlack); |
TVector3 *inters = new TVector3; |
surf->PropagateRay(*ray, out, inters); |
printf(" n1 = %f, n2 = %f\n", p_n1, p_n2); |
//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 *surfaceNormal = new CRay; |
surfaceNormal->SetColor(kBlue); |
surfaceNormal->Set(ray->GetR(), surf->GetN()); |
surfaceNormal->DrawS(cx, 1.0); |
CRay *incidentPolarization = new CRay; |
incidentPolarization->SetColor(kGreen); |
incidentPolarization->Set(ray->GetR(), p_pol); |
incidentPolarization->DrawS(cx, 1.0); |
printf(" GREEN: polarization vector\n"); |
CRay *surfaceNormal = new CRay; |
surfaceNormal->SetColor(kBlue); |
surfaceNormal->Set(ray->GetR(), surf->GetN()); |
surfaceNormal->DrawS(cx, 1.0); |
printf(" BLUE: surface normal vector\n"); |
} |
void ptt() |
{ |
for(double th = 0.0; th < 91.0; th += 5.0) { |
printf("%lf ", th); |
PolTest(th); |
} |
for(double th = 0.0; th < 91.0; th += 5.0) { |
printf("%lf ", th); |
PolTest(th); |
} |
} |
//------------------------------------------------------------------------------------------ |
// Propagate single ray and show the statistics |
185,11 → 184,11 |
void LGS(double xh = 0.0, double yh = 0.0, double zh = 0.0, double theta = 0.0, double phi = 0.0) |
{ |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
Init(); |
double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
} |
//------------------------------------------------------------------------------------------ |
// Propagate NN rays generated as grid and show the statistics |
196,13 → 195,13 |
void LGG(int NN=10, double theta=0.0, bool coupling=false) |
{ |
parameters.setCoupling(coupling); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
double izkoristek = Grid(detector, parameters, NN, theta); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
double izkoristek = Grid(detector, parameters, NN, theta); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
} |
//------------------------------------------------------------------------------------------ |
// Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics |
209,11 → 208,11 |
void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0) |
{ |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
Init(); |
double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
} |
//------------------------------------------------------------------------------------------ |
// Propagate NN rays isotropically generated in solid angle theta and show the statistics |
222,14 → 221,14 |
Init(); |
CDetector *detector = new CDetector(CENTER, parameters); |
//CDetector detector = new CDetector(); |
double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
//canvasDetector->cd(); |
//detector->Draw(); |
double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr); |
//printf("izkoristek = %.3lf\n", izkoristek); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, theta, izkoristek); |
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
//canvasDetector->cd(); |
//detector->Draw(); |
} |
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0) |
237,48 → 236,48 |
Init(); |
parameters.setCoupling(coupling); |
CDetector *detector = new CDetector(CENTER, parameters); |
const double theta = 18.5; |
double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, 18.5, izkoristek); |
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
//canvasDetector->cd(); |
//detector->Draw(); |
const double theta = 18.5; |
double izkoristek = beamtest(detector, parameters, NN, theta, phiMin, phiMax, nnrays, showr); |
PrintGuideHead(); |
PrintGuideStat(izkoristek); |
DrawData(detector, parameters, 18.5, izkoristek); |
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500); |
//canvasDetector->cd(); |
//detector->Draw(); |
} |
//------------------------------------------------------------------------------------------ |
void LGI_gap(int NN = 1e4, double min = 0.0, double max = 2.0, const int steps = 10, double theta = 30.0) |
{ |
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; |
//printf(" x_gap | acceptance \n"); |
PrintGuideHead(); |
for(int i=0; i<=steps; i++) { |
step[i] = min + i*(max - min)/(double)steps; |
parameters.setGap(step[i], 0.0, 0.0); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
//printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0); |
PrintGuideStat(acc[i]); |
} |
//char sbuff[256]; |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
show_3d = 0; show_data = 0; |
//printf(" x_gap | acceptance \n"); |
PrintGuideHead(); |
for(int i=0; i<=steps; i++) { |
step[i] = min + i*(max - min)/(double)steps; |
parameters.setGap(step[i], 0.0, 0.0); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
//printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0); |
PrintGuideStat(acc[i]); |
} |
//char sbuff[256]; |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
DrawAcc(steps+1, step, acc, (char*)"; x_gap;Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
288,28 → 287,28 |
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; |
//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; |
} |
320,29 → 319,29 |
//int tmp3d = show_3d, tmpdata = show_data; |
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, 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; |
} |
//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; |
//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, 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; |
} |
//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; |
} |
//------------------------------------------------------------------------------------------ |
354,19 → 353,19 |
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; |
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("."); |
} |
printf("\n"); |
//DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max); |
379,73 → 378,73 |
// 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) |
{ |
//char sbuff[256]; |
//char sbuff[256]; |
//show_3d = 0; // don't show simulations |
//show_data = 1; |
//double d = detector->GetD(); |
const double b = parameters.getB(); // upper side of LG |
//const double SiPM = 3.0; // the length of the detector itself |
//double reflectivity = detector->guide->getR(); |
TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max); |
// Use the Fresnel eq. instead of fixed reflectivity 96% |
//detector->SetFresnel(1); |
//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, b, d); |
CDetector *detector = new CDetector(CENTER, parameters); |
//detector->guide->setLG(y, M, x); |
//Init(); exclude simulation |
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 | ", d, a); |
//PrintGuideStat(acceptance); |
delete detector; //works fine, 50x50 grid takes ~4MB of RAM |
} |
printf("."); |
} |
printf("\n"); |
//show_3d = 0; // don't show simulations |
//show_data = 1; |
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(); |
file->Close(); |
//delete file; |
//hAcceptance->SetContour(100); |
//gStyle->SetPalette(1,0); |
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); |
//double d = detector->GetD(); |
const double b = parameters.getB(); // upper side of LG |
//const double SiPM = 3.0; // the length of the detector itself |
//double reflectivity = detector->guide->getR(); |
TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max); |
// Use the Fresnel eq. instead of fixed reflectivity 96% |
//detector->SetFresnel(1); |
//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, b, d); |
CDetector *detector = new CDetector(CENTER, parameters); |
//detector->guide->setLG(y, M, x); |
//Init(); exclude simulation |
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 | ", 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, 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(); |
file->Close(); |
//delete file; |
//hAcceptance->SetContour(100); |
//gStyle->SetPalette(1,0); |
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); |
} |
//------------------------------------------------------------------------------------------ |
459,55 → 458,55 |
{ |
TVector3 gapGuideSiPM(0.3, 0, 0); |
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters); |
int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
double sipm_d, M0; |
char sbuff[256]; |
int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
double sipm_d, M0; |
char sbuff[256]; |
show_3d = 1; show_data = 1; |
M0 = parameters.getM(); |
PrintGuideHead(); |
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]); |
//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, 0, theta); |
PrintGuideStat(acc[i]); |
sprintf(sbuff, "d%2d.gif", i); |
//c3dview->SaveAs(sbuff); |
} |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
//DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max); |
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480); |
TVirtualPad *pacc = cp->cd(0); |
pacc->SetRightMargin(0.05); |
pacc->SetLeftMargin(0.13); |
pacc->SetTopMargin(0.08); |
TGraph *gacceptance = new TGraph(steps+1, step, acc); |
gacceptance->SetTitle("; d [mm];Collection efficiency"); |
gacceptance->SetMarkerStyle(8); |
gacceptance->SetLineColor(12); |
gacceptance->SetLineWidth(1); |
(gacceptance->GetXaxis())->SetRangeUser(min, max); |
//(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
gacceptance->Draw("ACP"); |
show_3d = tmp3d; show_data = tmpdata; |
show_3d = 1; show_data = 1; |
M0 = parameters.getM(); |
PrintGuideHead(); |
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]); |
//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, 0, theta); |
PrintGuideStat(acc[i]); |
sprintf(sbuff, "d%2d.gif", i); |
//c3dview->SaveAs(sbuff); |
} |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
//DrawAcc(steps+1, step, acc, sbuff, (char*)"; d;Acceptance", min, max); |
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480); |
TVirtualPad *pacc = cp->cd(0); |
pacc->SetRightMargin(0.05); |
pacc->SetLeftMargin(0.13); |
pacc->SetTopMargin(0.08); |
TGraph *gacceptance = new TGraph(steps+1, step, acc); |
gacceptance->SetTitle("; d [mm];Collection efficiency"); |
gacceptance->SetMarkerStyle(8); |
gacceptance->SetLineColor(12); |
gacceptance->SetLineWidth(1); |
(gacceptance->GetXaxis())->SetRangeUser(min, max); |
//(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
gacceptance->Draw("ACP"); |
show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
// Acceptance vs. the length |
514,73 → 513,73 |
// The magnification ratio M rests the same all the time |
void LGI_d(int NN = 1e4, double min = 1, double max = 6.0, const int steps = 25, double theta = 30.0) |
{ |
int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
const double a = parameters.getA(); |
const double magnif = parameters.getM(); |
int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
const double a = parameters.getA(); |
const double magnif = parameters.getM(); |
//show_3d = show_in_steps; show_data = 1; |
// Use Fresnel equations |
//detector->SetFresnel(1); |
// Set glass (n=1.5) at the exit |
//detector->SetGlass(1,0); |
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, 0, theta); |
PrintGuideStat(acc[i]); |
delete detector; |
} |
//char sbuff[256]; |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
//show_3d = show_in_steps; show_data = 1; |
// Use Fresnel equations |
//detector->SetFresnel(1); |
// Set glass (n=1.5) at the exit |
//detector->SetGlass(1,0); |
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, 0, theta); |
PrintGuideStat(acc[i]); |
delete detector; |
} |
//char sbuff[256]; |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
DrawAcc(steps+1, step, acc, (char*)"; d [mm];Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
// Magnification optimization. a and d are fixed, M and b change in steps |
void LGI_M(int NN = 1e4, double min = 1.0, double max = 3.0, const int steps = 30, double theta = 30.0) |
{ |
int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
const double a = parameters.getA(); |
const double d = parameters.getD(); |
int tmp3d = show_3d, tmpdata = show_data; |
double step[steps], acc[steps]; |
show_3d = 0; show_data = 0; |
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, a*step[i], d); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
PrintGuideStat(acc[i]); |
} |
//char sbuff[256]; |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
const double a = parameters.getA(); |
const double d = parameters.getD(); |
show_3d = 0; show_data = 0; |
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, a*step[i], d); |
CDetector *detector = new CDetector(CENTER, parameters); |
Init(); |
acc[i] = RandIso(detector, parameters, NN, 0, theta); |
PrintGuideStat(acc[i]); |
} |
//char sbuff[256]; |
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap(y,z) = (%.1lf, %.1lf) | #theta = %.1lf", |
// detector->GetSiPM(), detector->GetM(), detector->GetD(), |
// (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta); |
DrawAcc(steps+1, step, acc, (char*)"; M;Acceptance", min, max); |
show_3d = tmp3d; show_data = tmpdata; |
} |
//------------------------------------------------------------------------------------------ |
588,11 → 587,11 |
void getParameters() |
{ |
printf("LIGHT GUIDE\n" |
" b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA()); |
" 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()); |
" 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()); |
" ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel()); |
return; |
} |
/lightguide/trunk/src/raySimulator.cpp |
---|
23,32 → 23,32 |
int show_3d = 0, show_data = 0; |
int draw_width = 2; |
// calls default constructor CVodnik.cpp |
// calls default constructor CVodnik.cpp |
//void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0) |
//{center.SetXYZ(x,y,z); detector = new CDetector(center);} |
//{center.SetXYZ(x,y,z); detector = new CDetector(center);} |
//void SetLGType(int in = 1, int side = 1, int out = 0) |
//{detector->SetLGType(in, side, out);} |
//{detector->SetLGType(in, side, out);} |
//void SetLG(double SiPM0 = 3.0, double M0 = 1.666, double d0 = 3.0, double n10 = 1.0, double n20 = 1.48, double n30 = 1.48, double R0 = 0.96) |
//{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);} |
//{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);} |
//void SetR(double R0) {detector->SetR(R0);} |
/* |
/* |
void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3) |
{detector->SetGlass(glass_on0, glass_d0);} |
void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0) |
{detector->SetGap(x_gap0 , y_gap0, z_gap0);} |
void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6) |
{detector->SetRCol(in, lg, out, gla);}; |
void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3) |
{detector->SetDCol(LG0, glass0, active0);}; |
void SetDWidth(int w = 1) {draw_width = w;} |
void SetFresnel(int b = 1) {detector->SetFresnel(b);} |
58,244 → 58,275 |
void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);} |
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);} |
*/ |
*/ |
//----------------------------------------------------------------------------- |
void Init(double rsys_scale = 7.0) |
void Init(double rsys_scale = 10.0) |
{ |
RTSetStyle(gStyle); |
gStyle->SetOptStat(10); |
if(show_3d) { |
double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale}; |
double rsys1[]={ rsys_scale, rsys_scale, rsys_scale}; |
c3dview = (TCanvas*)gROOT->FindObject("c3dview"); |
if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723); c3dview->SetFillColor(0);} |
else c3dview->Clear(); |
TView3D *view = new TView3D(1, rsys0, rsys1); |
view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]); |
//view->SetPerspective(); |
view->SetParallel(); |
view->FrontView(); |
view->Zoom();view->Zoom();view->Zoom();//view->Zoom(); |
//view->ToggleRulers(); |
} |
RTSetStyle(gStyle); |
gStyle->SetOptStat(10); |
if(show_3d) { |
double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale}; |
double rsys1[]={ rsys_scale, rsys_scale, rsys_scale}; |
c3dview = (TCanvas*)gROOT->FindObject("c3dview"); |
if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723); |
c3dview->SetFillColor(0);} |
else c3dview->Clear(); |
TView3D *view = new TView3D(1, rsys0, rsys1); |
view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]); |
//view->SetPerspective(); |
view->SetParallel(); |
view->FrontView(); |
view->Zoom();view->Zoom();view->Zoom();//view->Zoom(); |
//view->ToggleRulers(); |
} |
} |
//----------------------------------------------------------------------------- |
void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0) |
{ |
if(show_data) { |
char sbuff[256]; |
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) { |
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950); |
cdata->Divide(3,3); |
int cpc = 1; |
cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ"); |
cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ"); |
cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ"); |
(detector->GetHGlass())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ"); |
cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw(); |
}else{ |
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500); |
cdata->cd(0); (detector->GetHLaser())->Draw("COLZ"); |
cdata->SaveAs("cdata.gif"); |
} |
} |
if(show_data) { |
char sbuff[256]; |
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) { |
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950); |
cdata->Divide(3,3); |
int cpc = 1; |
cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ"); |
cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ"); |
cdata->cd(cpc++); //((detector->GetLG())->GetHOut())->Draw("COLZ"); |
(detector->GetHGlass())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ"); |
cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ"); |
cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw(); |
cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw(); |
}else{ |
RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500); |
cdata->cd(0); (detector->GetHLaser())->Draw("COLZ"); |
cdata->SaveAs("cdata.gif"); |
} |
} |
} |
//----------------------------------------------------------------------------- |
void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance", |
double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0) |
double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0) |
{ |
cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480); |
(cacc->cd(0))->SetRightMargin(0.05); |
(cacc->cd(0))->SetLeftMargin(0.13); |
(cacc->cd(0))->SetTopMargin(0.08); |
TGraph *gacceptance = new TGraph(n, x, y); |
gacceptance->SetTitle(htitle); |
gacceptance->SetMarkerStyle(8); |
gacceptance->SetLineColor(12); |
gacceptance->SetLineWidth(1); |
if(xmin < xmax) |
(gacceptance->GetXaxis())->SetRangeUser(xmin, xmax); |
(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
gacceptance->Draw("ACP"); |
//cacc->Print("acceptance.pdf"); |
//delete gacceptance; |
//delete cacc; |
cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480); |
(cacc->cd(0))->SetRightMargin(0.05); |
(cacc->cd(0))->SetLeftMargin(0.13); |
(cacc->cd(0))->SetTopMargin(0.08); |
TGraph *gacceptance = new TGraph(n, x, y); |
gacceptance->SetTitle(htitle); |
gacceptance->SetMarkerStyle(8); |
gacceptance->SetLineColor(12); |
gacceptance->SetLineWidth(1); |
if(xmin < xmax) |
(gacceptance->GetXaxis())->SetRangeUser(xmin, xmax); |
(gacceptance->GetYaxis())->SetRangeUser(ymin, ymax); |
gacceptance->Draw("ACP"); |
//cacc->Print("acceptance.pdf"); |
//delete gacceptance; |
//delete cacc; |
} |
//----------------------------------------------------------------------------- |
void PrintGlassStat(CDetector *detector) |
{ |
int glass_in = ((detector->GetLG())->GetHOut())->GetEntries(); |
int glass_out = (detector->GetHGlass())->GetEntries(); |
printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in); |
int glass_in = ((detector->GetLG())->GetHOut())->GetEntries(); |
int glass_out = (detector->GetHGlass())->GetEntries(); |
printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in); |
} |
//----------------------------------------------------------------------------- |
void PrintGuideHead() |
{ |
//printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n"); |
printf("Acceptance\n"); |
//printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n"); |
printf("Acceptance\n"); |
} |
//----------------------------------------------------------------------------- |
void PrintGuideStat(double acceptance) |
{ |
/*int n_active = (detector->GetHActive())->GetEntries(); |
/*int n_active = (detector->GetHActive())->GetEntries(); |
int n_enter = (detector->GetLG())->GetEnteranceHits(); |
double izkoristek = n_active/(double)n_enter; |
int fates[7]; (detector->GetLG())->GetVFate(fates); |
// printf(" gap (x,y,z) | theta phi | Back | Loss | Total| Miss | Exit | Acceptance\n"); |
printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |", |
(detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi); |
for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter); |
printf("%5.1lf\n", 100.0*izkoristek);*/ |
//double n_active = (detector->GetHActive())->GetEntries(); |
//double r_acc = 100.0 * n_active / n_in; |
double r_acc = 100.0 * acceptance; |
printf("%7.3lf\n", r_acc); |
//double n_active = (detector->GetHActive())->GetEntries(); |
//double r_acc = 100.0 * n_active / n_in; |
double r_acc = 100.0 * acceptance; |
printf("%7.3lf\n", r_acc); |
} |
//----------------------------------------------------------------------------- |
//----------------------------------------------------------------------------- |
// en zarek |
double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0) |
double Single(CDetector *detector, DetectorParameters& parameters, |
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0) |
{ |
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
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); |
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(), |
TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi)); |
CRay *ray1 = new CRay; |
detector->Propagate(*ray0, ray1, show_3d); |
delete ray0; |
delete ray1; |
return (detector->GetHActive())->GetEntries() / (double)(1); |
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); |
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(), |
TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), |
TMath::Sin(theta)*TMath::Cos(phi)); |
// Set z-polarization == vertical |
TVector3 polarization(0,0,1); |
polarization.RotateY(-theta); |
polarization.RotateZ(phi); |
if (polarization.Dot(ray0->GetK()) > 1e-5) printf("ERROR: pol not perep\n"); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
detector->Propagate(*ray0, ray1, show_3d); |
CRay *incidentPolarization = new CRay; |
incidentPolarization->SetColor(kGreen); |
TVector3 drawPosition = ray0->GetR(); |
drawPosition.SetX(drawPosition.X() - 5); |
incidentPolarization->Set(drawPosition, polarization); |
incidentPolarization->DrawS(drawPosition.X(), 1); |
TVector3 outPolarization = ray1->GetP(); |
drawPosition = ray1->GetR(); |
drawPosition.SetX(drawPosition.X() + 5); |
CRay* rayPol = new CRay(drawPosition, outPolarization); |
rayPol->SetColor(kBlack); |
rayPol->DrawS(drawPosition.X(), 1); |
delete ray0; |
delete ray1; |
return (detector->GetHActive())->GetEntries() / (double)(1); |
} |
//----------------------------------------------------------------------------- |
// zarki, razporejeni v mrezi |
double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0) |
double Grid(CDetector *detector, DetectorParameters& parameters, |
int NN = 10, double theta = 0.0) |
{ |
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
theta = Pi()*theta/180.0; |
if(theta < 1e-6) theta = 1e-6; |
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
theta = Pi()*theta/180.0; |
if(theta < 1e-6) theta = 1e-6; |
const double b = parameters.getB(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
for(int i = 0; i < NN; i++) { |
for(int j = 0; j < NN; j++) { |
CRay *ray0 = new CRay(CENTER.x() - offset, |
0.99*(i*b/NN + b/(2.0*NN) - b/2.0), |
0.99*(j*b/NN + b/(2.0*NN) - b/2.0), |
TMath::Cos(theta), |
0.0, |
TMath::Sin(theta)); |
CRay *ray1 = new CRay; |
if(i == (NN/2)) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
delete ray0; |
delete ray1; |
} |
} |
double acceptance = 0.0; |
/* |
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
const double b = parameters.getB(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
for(int i = 0; i < NN; i++) { |
for(int j = 0; j < NN; j++) { |
CRay *ray0 = new CRay(CENTER.x() - offset, |
0.99*(i*b/NN + b/(2.0*NN) - b/2.0), |
0.99*(j*b/NN + b/(2.0*NN) - b/2.0), |
TMath::Cos(theta), |
0.0, |
TMath::Sin(theta)); |
TVector3 polarization(0, 1, 0); |
polarization.RotateY(-theta); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
if(i == (NN/2)) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
delete ray0; |
delete ray1; |
} |
} |
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; |
*/ |
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, double theta, double phi, int show_rays) |
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; |
if(theta < MARGIN) theta = MARGIN; |
phi = phi*3.14159265358979312/180.0; |
if(phi < MARGIN) phi = MARGIN; |
TDatime now; |
TRandom rand; |
rand.SetSeed(now.Get()); |
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
for(int i = 0; i < NN; i++) { |
double rx = CENTER.x() - offset; |
double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
double rl = TMath::Cos(theta); |
double rm = TMath::Sin(theta)*TMath::Sin(phi); |
double rn = TMath::Sin(theta)*TMath::Cos(phi); |
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; |
} |
double acceptance = 0.0; |
/* |
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
theta = theta*3.14159265358979312/180.0; |
if(theta < MARGIN) theta = MARGIN; |
phi = phi*3.14159265358979312/180.0; |
if(phi < MARGIN) phi = MARGIN; |
TDatime now; |
TRandom rand; |
rand.SetSeed(now.Get()); |
//detector->Init(); |
if(show_3d) detector->Draw(draw_width); |
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
for(int i = 0; i < NN; i++) { |
double rx = CENTER.x() - offset; |
double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
double rl = TMath::Cos(theta); |
double rm = TMath::Sin(theta)*TMath::Sin(phi); |
double rn = TMath::Sin(theta)*TMath::Cos(phi); |
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn); |
TVector3 polarization(0, 0, 1); |
polarization.RotateY(-theta); |
polarization.RotateZ(phi); |
ray0->setPolarization(polarization); |
CRay *ray1 = new CRay; |
if(i < show_rays) |
detector->Propagate(*ray0, ray1, show_3d); |
else |
detector->Propagate(*ray0, ray1, 0); |
delete ray1; |
delete ray0; |
} |
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; |
acceptance = (detector->GetHActive())->GetEntries() / (double) NN; |
return acceptance; |
} |
//----------------------------------------------------------------------------- |
// zarki, izotropno porazdeljeni znotraj kota theta |
304,201 → 335,207 |
// 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) |
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; |
theta = theta*3.14159265358979312/180.0; |
if(theta < MARGIN) theta = MARGIN; |
TDatime now; |
TRandom rand; |
rand.SetSeed(now.Get()); |
double rx, ry, rz, rl, rm, rn; |
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); |
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48); |
double pi = 3.14159265358979312; |
theta = theta*3.14159265358979312/180.0; |
if(theta < MARGIN) theta = MARGIN; |
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*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); |
TDatime now; |
TRandom rand; |
rand.SetSeed(now.Get()); |
double rx, ry, rz, rl, rm, rn; |
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); |
//detector->Init(); |
if (show_3d) detector->Draw(draw_width); |
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn; |
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi; |
hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*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); |
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
for(int i = 0; i < NN; i++) { |
rx = CENTER.x() - offset; |
ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
rphi = rand.Uniform(0.0, 2.0*pi); |
//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); |
rn = TMath::Sin(rtheta)*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; |
/* |
//detector->Init(); |
if (show_3d) detector->Draw(draw_width); |
double SiPM = parameters.getA(); |
double M = parameters.getM(); |
double offset = (parameters._plateOn ? parameters._plateWidth : 0); |
for(int i = 0; i < NN; i++) { |
rx = CENTER.x() - offset; |
ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0); |
rphi = rand.Uniform(0.0, 2.0*pi); |
//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); |
rn = TMath::Sin(rtheta)*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; |
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 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; |
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); |
double pi = 3.14159265358979312; |
theta *= pi/180.0; |
if(theta < MARGIN) theta = MARGIN; |
phiMin *= pi/180.0; |
phiMax *= pi/180.0; |
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); |
TDatime now; |
TRandom rand; |
rand.SetSeed(now.Get()); |
double rx, ry, rz, rl, rm, rn; |
double rphi; |
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; |
/* |
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); |
// inital polarizaton |
TVector3 polarization(0, 0, 1); |
polarization.RotateY(-theta); |
polarization.RotateX(rphi); |
ray0->setPolarization(polarization); |
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; |
return acceptance; |
} |
/lightguide/trunk/src/guide.cpp |
---|
20,12 → 20,10 |
if(in >= 0.0) return 1; |
else return -1; |
} |
//================================================================================= |
//----------------------------------------------------------------------------- |
void CRay::Set(TVector3 r0, TVector3 n0) |
void CRay::Set(TVector3 r0, TVector3 k0) |
{ |
r = r0; n = n0.Unit(); |
r = r0; k = k0.Unit(); |
} |
//----------------------------------------------------------------------------- |
//void CRay::Set(double x0, double y0, double z0, double l0, double m0, double n0) |
42,14 → 40,12 |
n.SetXYZ(p.GetN().x(), p.GetN().y(), p.GetN().z()); |
return *this; |
} */ |
//----------------------------------------------------------------------------- |
void CRay::Print() |
{ |
printf("---> CRay::Print() <---\n"); |
printf("(x,y,z)=(%.2lf, %.2lf, %.2lf); (l,m,n)=(%.2lf, %.2lf, %.2lf)\n", |
r.x(), r.y(), r.z(), n.x(), n.y(), n.z()); |
r.x(), r.y(), r.z(), k.x(), k.y(), k.z()); |
} |
//----------------------------------------------------------------------------- |
void CRay::Draw() |
{ |
double t = 50.0; |
56,54 → 52,50 |
TPolyLine3D *line3d = new TPolyLine3D(2); |
//line3d->SetPoint(0, r.x() - t*n.x(), r.y() - t*n.y(), r.z() - t*n.z()); |
line3d->SetPoint(0, r.x(), r.y(), r.z()); |
line3d->SetPoint(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.z()); |
line3d->SetPoint(1, r.x() + t*k.x(), k.y() + t*k.y(), r.z() + t*k.z()); |
line3d->SetLineWidth(1); |
line3d->SetLineColor(color); |
line3d->Draw(); |
} |
//----------------------------------------------------------------------------- |
void CRay::Draw(double x_from, double x_to) |
{ |
double A1, A2; |
TPolyLine3D *line3d = new TPolyLine3D(2); |
if(n.x() < MARGIN) { |
if(k.x() < MARGIN) { |
A1 = A2 = 0.0; |
} else { |
A1 = (x_from - r.x())/n.x(); |
A2 = (x_to - r.x())/n.x(); |
A1 = (x_from - r.x())/k.x(); |
A2 = (x_to - r.x())/k.x(); |
} |
line3d->SetPoint(0, x_from, A1*n.y()+r.y(), A1*n.z()+r.z()); |
line3d->SetPoint(1, x_to, A2*n.y()+r.y(), A2*n.z()+r.z()); |
line3d->SetPoint(0, x_from, A1*k.y()+r.y(), A1*k.z()+r.z()); |
line3d->SetPoint(1, x_to, A2*k.y()+r.y(), A2*k.z()+r.z()); |
line3d->SetLineWidth(1); |
line3d->SetLineColor(color); |
line3d->Draw(); |
} |
//----------------------------------------------------------------------------- |
void CRay::DrawS(double x_from, double t) |
{ |
double A1; |
TPolyLine3D *line3d = new TPolyLine3D(2); |
if(n.x() < MARGIN) |
if(k.x() < MARGIN) |
A1 = 0.0; |
else |
A1 = (x_from - r.x())/n.x(); |
A1 = (x_from - r.x())/k.x(); |
line3d->SetPoint(0, x_from, A1*n.y()+r.y(), A1*n.z()+r.z()); |
line3d->SetPoint(1, r.x() + t*n.x(), r.y() + t*n.y(), r.z() + t*n.z()); |
line3d->SetPoint(0, x_from, A1*k.y()+r.y(), A1*k.z()+r.z()); |
line3d->SetPoint(1, r.x() + t*k.x(), r.y() + t*k.y(), r.z() + t*k.z()); |
line3d->SetLineWidth(1); |
line3d->SetLineColor(color); |
line3d->Draw(); |
} |
//================================================================================= |
//================================================================================= |
CPlane4::CPlane4() : |
n(TVector3(1.0, 0.0, 0.0)), |
A(0), |
117,7 → 109,6 |
for(int i=0;i<4;i++) edge[i] = TVector3(0,0,0); |
for(int i=0;i<4;i++) angle_r[i] = 0; |
}; |
//----------------------------------------------------------------------------- |
CPlane4::CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4) |
{ |
//Set(r1, r2, r3, r4); |
209,7 → 200,7 |
N.SetXYZ(A,B,C); |
num = N*ray.GetR() + D; |
den = N*ray.GetN(); |
den = N*ray.GetK(); |
if (dbg) printf("t = %6.3lf / %6.3lf = %6.3lf\n", num, den, num/den); |
229,7 → 220,7 |
t = num / den; |
tmp = ray.GetR(); |
tmp -= t*ray.GetN(); |
tmp -= t*ray.GetK(); |
*vec = tmp; |
return 1; |
} |
275,7 → 266,6 |
return 1; |
} |
//----------------------------------------------------------------------------- |
int CPlane4::TestIntersection(CRay in) |
{ |
TVector3 tmp; |
286,7 → 276,6 |
return 0; |
} |
//----------------------------------------------------------------------------- |
int CPlane4::TestIntersection(TVector3 *vec, CRay in) |
{ |
TVector3 tmp; |
299,7 → 288,6 |
return 0; |
} |
//----------------------------------------------------------------------------- |
void CPlane4::Print() |
{ |
printf("--- CPlane4::Print() ---\n"); |
309,7 → 297,6 |
for(int i=0;i<4;i++) printf(" edge[%d] = (%lf, %lf, %lf)\n", i, edge[i].x(), edge[i].y(), edge[i].z()); |
for(int i=0;i<4;i++) printf(" angle[%d] = %lf\n", i, angle_r[i]*DEGREE); |
} |
//----------------------------------------------------------------------------- |
void CPlane4::Draw(int color, int width) |
{ |
TPolyLine3D *line3d = new TPolyLine3D(5); |
320,10 → 307,8 |
line3d->Draw(); |
} |
//================================================================================= |
//================================================================================= |
CSurface::CSurface(int type0): |
type(type0) |
{ |
342,7 → 327,6 |
SetFresnel(); |
} |
//----------------------------------------------------------------------------- |
CSurface::CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity) |
{ |
TDatime now; |
355,7 → 339,6 |
SetFresnel(); |
} |
//----------------------------------------------------------------------------- |
CSurface::CSurface(int type0, TVector3 *vr, double n10, double n20, double reflectivity) |
{ |
TDatime now; |
368,7 → 351,6 |
SetFresnel(); |
} |
//----------------------------------------------------------------------------- |
void CSurface::SetIndex(double n10, double n20) |
{ |
n1 = n10; n2 = n20; n1_n2 = n1/n2; |
410,19 → 392,22 |
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 |
TVector3 pol_t = in.GetP(); // incident polarization |
int sign_n; // sign of normal direction vs. inbound ray |
double cosTN; // debug |
if(fresnel) { |
// p-polarization unit vector v_te |
// Decomposition of incident polarization vector |
// using unit vectors v_tm & v_te |
// in a_tm and a_te components |
//if(fresnel) { |
// s-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()); |
// k in this notation is in.GetK() |
v_te = n.Cross(in.GetK()); |
v_te = v_te.Unit(); |
v_tm = -v_te.Cross(in.GetN()); |
v_tm = -v_te.Cross(in.GetK()); |
v_tm = v_tm.Unit(); |
if(dbg) { |
printf(" v_te = "); printv(v_te); |
430,12 → 415,13 |
} |
double cosAf = v_te * in.GetP(); |
if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, TMath::ACos(cosAf)*DEGREE); |
double alpha = acos(cosAf); |
if(dbg) printf(" cosAf = %lf (Af = %lf)\n", cosAf, alpha*DEGREE); |
a_te = cosAf; |
a_tm = TMath::Sqrt(1 - cosAf*cosAf); |
if(dbg) printf(" a_te = %lf, a_tm = %lf\n", a_te, a_tm); |
} |
//} |
// ---------------------------------------------------------------------------- |
// reflection probability |
447,7 → 433,7 |
// --------------- refraction from n1 to n2 ----------------------------------- |
// ---------------------------------------------------------------------------- |
case SURF_REFRA: |
cosTi = in.GetN() * n; |
cosTi = in.GetK() * n; |
if(dbg) printf(" cosTi = %lf (Ti = %lf)\n", cosTi, TMath::ACos(cosTi)*DEGREE); |
sign_n = -sign(cosTi); |
if(dbg) printf(" sign_n = %d\n", sign_n); |
467,7 → 453,7 |
// If n1>n2 and theta>thetaCritical, total reflection |
if(cosTi < cosTtotal) { |
if(dbg) printf(" TOTAL\n"); |
transmit = in.GetN() + sign_n*2*cosTi*n; |
transmit = in.GetK() + sign_n*2*cosTi*n; |
if(dbg) { |
cosTN = TMath::Abs(transmit.Unit() * n); |
475,9 → 461,24 |
} |
out->Set(intersect, transmit); |
// Shift? |
pol_t = -in.GetP() + sign_n*2*cosTi*n; |
out->SetPolarization(pol_t); |
// Shift implemented, but only linear polarization is implemented |
if (dbg) printf("CSurface: Propagate TOTAL\n"); |
v_tm_t = -v_te.Cross(transmit); |
v_tm_t = v_tm_t.Unit(); |
// shift the p and s components |
double n12 = N1_N2(-sign_n); |
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; |
a_tm = sin(alpha); |
a_te = cos(alpha); |
if (dbg) printf(" deltaP = %f deltaS = %f; new a_tm = %f, a_te = %f", |
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); |
return REFLECTION; |
} else { |
// reflection or refraction according to Fresnel equations |
486,7 → 487,7 |
cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2))); |
if(dbg) printf(" cosTt = %lf (Tt = %lf) \n", cosTt, TMath::ACos(cosTt)*DEGREE); |
transmit = N1_N2(sign_n)*in.GetN() + sign_n*(N1_N2(sign_n)*cosTi - cosTt)*n; |
transmit = N1_N2(sign_n)*in.GetK() + sign_n*(N1_N2(sign_n)*cosTi - cosTt)*n; |
if(dbg) {printf(" transmit.Unit() = "); printv(transmit.Unit());} |
if(dbg) { |
cosTN = transmit.Unit() * n; |
503,7 → 504,7 |
// 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; |
pol_t = a_te*(1.0 - TMath::Abs(r_te))*v_te + a_tm*(1.0 - TMath::Abs(r_tm))*v_tm_t; |
if(dbg) { |
printf(" v_tm_t = "); printv(v_tm_t); |
521,14 → 522,17 |
if(p_ref >= R_f) { // se lomi |
if (dbg) printf(" SURFACE REFRACTED. Return.\n"); |
out->Set(intersect, transmit); |
out->SetPolarization(pol_t); |
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); |
transmit = in.GetN() + sign_n*2*cosTi*n; |
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); |
pol_t = -in.GetP() + sign_n*2*cosTi*n; |
out->SetPolarization(pol_t); |
//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); |
return REFLECTION; |
} |
541,12 → 545,12 |
case SURF_REFLE: |
p_ref = rand.Uniform(0.0, 1.0); |
if(p_ref < reflection) { // se odbije |
cosTi = in.GetN() * n; |
transmit = in.GetN() - 2*cosTi*n; |
cosTi = in.GetK() * n; |
transmit = in.GetK() - 2*cosTi*n; |
out->Set(intersect, transmit); |
return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb |
} else { // se ne odbije |
transmit = in.GetN(); |
transmit = in.GetK(); |
out->Set(intersect, transmit); |
return ABSORBED; |
} |
556,17 → 560,17 |
case SURF_IMPER: |
p_ref = rand.Uniform(0.0, 1.0); |
if(p_ref < reflection) { // se odbije |
cosTi = in.GetN() * n; |
cosTi = in.GetK() * n; |
if(TMath::Abs(cosTi) < cosTtotal) { // totalni odboj |
transmit = in.GetN() - 2*cosTi*n; |
transmit = in.GetK() - 2*cosTi*n; |
out->Set(intersect, transmit); |
} else { // ni tot. odboja |
transmit = in.GetN(); |
transmit = in.GetK(); |
out->Set(intersect, transmit); |
return ABSORBED; |
} |
} else { // se ne odbije |
transmit = in.GetN(); |
transmit = in.GetK(); |
out->Set(intersect, transmit); |
return ABSORBED; |
} |
579,10 → 583,7 |
return REFRACTION; |
} |
//================================================================================= |
//================================================================================= |
Guide::Guide(TVector3 center0, DetectorParameters ¶meters) : |
_d(parameters.getD()), |
_n1(parameters.getN1()), |
640,7 → 641,7 |
TVector3 activePosition(center); |
activePosition += TVector3(_d, 0, 0); |
TVector3 normal(1,0,0); |
grease = new CPlaneR(activePosition, normal, a/2.0); |
grease = new CPlaneR(activePosition, normal, 0.95*a/2.0); |
if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1); |
750,7 → 751,7 |
if (!hitActive) propagation = noCoupling->PropagateRay(ray0, &ray1, &vec1); |
} |
// check on which side the vector is? |
TVector3 ray = ray1.GetN(); |
TVector3 ray = ray1.GetK(); |
TVector3 exitNormal = s_side[5]->GetN(); |
if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal); |
if (ray*exitNormal > 0) { |
806,17 → 807,14 |
*out = ray0; |
return fate; |
} |
//----------------------------------------------------------------------------- |
void Guide::GetVFate(int *out) |
{ |
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); |
} |
//----------------------------------------------------------------------------- |
void Guide::DrawSkel(int color, int width) |
{ |
TPolyLine3D *line3d = new TPolyLine3D(2); |
828,9 → 826,8 |
line3d->DrawClone(); |
} |
} |
//================================================================================= |
//================================================================================= |
int CPlaneR::TestIntersection(TVector3 *vec, CRay ray) |
{ |
double num, den; //stevec, imenovalec |
842,7 → 839,7 |
double D = - n*center; |
num = n*ray.GetR() + D; |
den = n*ray.GetN(); |
den = n*ray.GetK(); |
if(dbg) printf("D = %.4lf | num = %.4lf | den = %.4lf\n", D, num, den); |
858,7 → 855,7 |
if(dbg) printf("t = %.4lf | ", t); |
tmp = ray.GetR(); |
tmp -= t*ray.GetN(); |
tmp -= t*ray.GetK(); |
*vec = tmp; |
if(dbg) {printv(tmp); printf(" | Rv = %.4lf <> R = %.4lf\n", ((tmp - center).Mag()), _r);} |
869,7 → 866,7 |
else |
return 0; |
} |
//----------------------------------------------------------------------------- |
void CPlaneR::Draw(int color, int width) |
{ |
const int NN = 32; |
888,10 → 885,8 |
} |
arc->Draw(); |
} |
//================================================================================= |
//================================================================================= |
CDetector::CDetector(TVector3 center0, DetectorParameters& parameters) : |
center(center0), |
glass_on(parameters.getGlassOn()), |
1048,6 → 1043,8 |
rayout->DrawS(center.x() - _plateWidth, -10.0); |
} |
else if(fatePlate == backreflected){ |
rayout->SetColor(kBlack); |
rayout->DrawS(center.x() - _plateWidth, 7.0); |
if (dbg) printf("Backreflected at plate!\n"); |
} |
else { |
1110,8 → 1107,8 |
if(fate == missed) { |
if (dbg) printf("Detector: fate=missed\n"); |
TVector3 r = ray1->GetR(); |
TVector3 n = ray1->GetN(); |
ray1->Set(r,n); |
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++) { |
1192,7 → 1189,7 |
delete rayin; |
return fate; |
} |
//----------------------------------------------------------------------------- |
void CDetector::Draw(int width) |
{ |
if(guide_on) { |
1210,8 → 1207,8 |
//window_circle->Draw(col_glass, width); |
active->Draw(col_active, width); |
} |
//================================================================================= |
Plate::Plate(DetectorParameters& parameters) |
{ |
TVector3 center = CENTER; |
1284,6 → 1281,7 |
fate = missed; |
} else if(result == REFLECTION) { |
if (dbg) printf("PLATE: reflected\n"); |
ray0 = ray1; |
fate = backreflected; |
} else { |
points[0] = ray1.GetR(); |
1300,6 → 1298,7 |
} |
points[n_odb] = vec1; |
if(inters_i == 0) { |
ray0 = ray1; |
fate = backreflected; |
break;} // backreflection |
1333,6 → 1332,5 |
*out = ray0; |
return fate; |
}; |
//=============================================================================================================================== <<<<<<<< |