Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 72 → Rev 73

/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 &parameters) :
_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;
};
//=============================================================================================================================== <<<<<<<<