Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 69 → Rev 70

/lightguide/trunk/src/RTUtil.cpp
1,6 → 1,9
//##########################################################################################
#include "TGaxis.h"
#include "TColor.h"
 
#include "include/RTUtil.h"
 
//##########################################################################################
void RTSetStyle(TStyle *style)
{
style->SetStatBorderSize(1);
13,12 → 16,23
style->SetStatColor(0);
style->SetPalette(1, 0);
style->SetMarkerStyle(kFullDotLarge);
//style->SetMarkerSize(7);
style->SetOptStat("ne");
style->SetOptFit(1);
style->SetPadTopMargin(0.10);
style->SetPadBottomMargin(0.10);
style->SetPadLeftMargin(0.10);
style->SetPadRightMargin(0.12);
style->SetPadBottomMargin(0.12);
style->SetPadLeftMargin(0.12);
style->SetPadRightMargin(0.15);
style->SetTitleOffset(1.3, "y");
style->SetTitleOffset(1.5, "y");
style->SetPalette(1, 0);
style->SetPaperSize(TStyle::kA4);
TGaxis::SetMaxDigits(4);
}
//##########################################################################################
RTCanvas::RTCanvas()
49,6 → 63,19
pad->Divide(nx, ny, 0.003, 0.005);
}
//------------------------------------------------------------------------------------------
void RTCanvas::Divide(int np)
{
if( np==2 ) pad->Divide(1, 2, 0.003, 0.005);
else if( 2<np && np<=4 ) pad->Divide(2, 2, 0.003, 0.005);
else if( 4<np && np<=6 ) pad->Divide(2, 3, 0.003, 0.005);
else if( 6<np && np<=8 ) pad->Divide(2, 4, 0.003, 0.005);
else if( np==9 ) pad->Divide(3, 3, 0.003, 0.005);
else if( 9<np && np<=12) pad->Divide(3, 4, 0.003, 0.005);
else if(12<np && np<=16) pad->Divide(4, 4, 0.003, 0.005);
else if(16<np && np<=25) pad->Divide(5, 5, 0.003, 0.005);
else if(25<np && np<=32) pad->Divide(4, 8, 0.003, 0.005);
}
//------------------------------------------------------------------------------------------
TPad* RTCanvas::cd(int i)
{
return (TPad*)(pad->cd(i));
58,6 → 85,22
{
can->SaveAs(filename);
}
//------------------------------------------------------------------------------------------
void RTCanvas::Update()
{
can->Update();
}
 
void SetGS()
{
const Int_t Number = 2;
Double_t Red[Number] = {1.0, 0.0};
Double_t Green[Number] = {1.0, 0.0};
Double_t Blue[Number] = {1.0, 0.0};
Double_t Stops[Number] = {0.0, 1.0};
Int_t nb = 50;
TColor::CreateGradientColorTable(Number, Stops, Red, Green, Blue, nb);
}
//##########################################################################################
 
 
/lightguide/trunk/src/userFunctions.cpp
9,29 → 9,38
//extern int show_data;
 
// Set the simulation parameters
void Show3D(int b) {show_3d = b;}
void ShowData(int b) {show_data = b;}
void showVisual(int b) {show_3d = b;}
void showData(int b) {show_data = b;}
 
// Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap)
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.52, TVector3(0.3, 0, 0));
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0));
// Print the detector parameters
void getParameters();
 
//void SetLGType(int in = 1, int side = 1, int out = 0)
//{detector->SetLGType(in, side, out);}
 
void SetLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n10 = 1.0, double n20 = 1.53, double n30 = 1.52)
{ parameters.setGuide(SiPM0, b0, d0, n10, n20, n30); }
void setLG(double SiPM0 = 3.0, double b0 = 5.0, double d0 = 3.0, double n1 = 1.0, double n2 = 1.53, double n3 = 1.46)
{ parameters.setGuide(SiPM0, b0, d0);
parameters.setIndices(n1, n2, n3); }
void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
void setGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
{ parameters.setGap(x_gap0, y_gap0, z_gap0); }
void SetGlass(double glassOn, double glassD)
void setGlass(double glassOn, double glassD)
{ parameters.setGlass(glassOn, glassD); };
 
void SetPlate(int plateOn = 1, double plateWidth = 1)
void setPlate(int plateOn = 1, double plateWidth = 1)
{ parameters.setPlate(plateOn, plateWidth); }
void SetFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
void setFresnel(int fresnel = 1) { parameters.setFresnel(fresnel); }
 
// Refractive Indices
// n1 - around light guide - air
// n2 - light guide (plate) material - k9 glass
// n3 - the material at the exit - optical grease, epoxy, air, etc.
void setIndices(double n1, double n2, double n3) {parameters.setIndices(n1, n2, n3);}
 
int save_ary = 0;
//------------------------------------------------------------------------------------------
void Help()
107,11 → 116,16
void SetPol(double x, double y, double z)
{p_pol0.SetXYZ(x,y,z);}
 
// Test function
// creates an optical boundary surface
// shows the propagation of light ray
// and polarization state of the incident ray (green)
// and surface normal (blue)
void PolTest(double theta = 0.0)
{
int p_type = 1;
double p_n1 = 1.5;
double p_n2 = 1.0;
double p_n1 = parameters.getN1();
double p_n2 = parameters.getN2();
theta = 3.141593*theta/180.0; if(theta < 1e-6) theta = 1e-6;
TVector3 p_pol;
130,11 → 144,11
#define SURF_IMPER 4
*/
CSurface *surf = new CSurface(p_type, vodnik_edge, p_n1, p_n2, 0.96); //surf->FlipN();
surf->SetFresnel(0);
surf->SetFresnel(1);
surf->Draw();
CRay *ray = new CRay(cx, 0.0, 0.0, TMath::Cos(theta), 0.0, TMath::Sin(theta));
ray->SetColor(1);
ray->SetColor(kBlack);
//p_pol = rotatey(p_pol0, -theta);
p_pol = p_pol0; p_pol.RotateY(-theta);
//printf("p_pol = "); printv(p_pol);
142,18 → 156,21
ray->DrawS(cx, -5.0);
CRay *out = new CRay; out->SetColor(2);
CRay *out = new CRay; out->SetColor(kRed);
TVector3 *inters = new TVector3;
int fate = surf->PropagateRay(*ray, out, inters);
if(fate == 1) out->DrawS(cx, 5.0);
surf->PropagateRay(*ray, out, inters);
//if(fate == 1) out->DrawS(cx, 5.0);
out->DrawS(cx, 5.0);
CRay *incidentPolarization = new CRay;
incidentPolarization->SetColor(kGreen);
incidentPolarization->Set(ray->GetR(), p_pol);
incidentPolarization->DrawS(cx, 1.0);
CRay *pp = new CRay; pp->SetColor(3);
pp->Set(ray->GetR(), p_pol);
pp->DrawS(cx, 2.0);
CRay *pn = new CRay; pn->SetColor(4);
pn->Set(ray->GetR(), surf->GetN());
pn->DrawS(cx, 3.0);
CRay *surfaceNormal = new CRay;
surfaceNormal->SetColor(kBlue);
surfaceNormal->Set(ray->GetR(), surf->GetN());
surfaceNormal->DrawS(cx, 1.0);
}
 
void ptt()
192,7 → 209,7
{
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
Init();
double izkoristek = RandYZ(detector, parameters, NN, theta, phi);
double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
//printf("izkoristek = %.3lf\n", izkoristek);
PrintGuideHead(); PrintGuideStat(izkoristek);
DrawData(detector, parameters, theta, izkoristek);
204,9 → 221,10
Init();
CDetector *detector = new CDetector(CENTER, parameters);
//CDetector detector = new CDetector();
double izkoristek = RandIso(detector, parameters, NN, theta, nnrays, showr);
double izkoristek = RandIso(detector, parameters, NN, 0.0, theta, nnrays, showr);
//printf("izkoristek = %.3lf\n", izkoristek);
PrintGuideHead(); PrintGuideStat(izkoristek);
PrintGuideHead();
PrintGuideStat(izkoristek);
DrawData(detector, parameters, theta, izkoristek);
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
//canvasDetector->cd();
213,6 → 231,20
//detector->Draw();
}
 
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, int nnrays = 30, int showr = 0)
{
Init();
CDetector *detector = new CDetector(CENTER, parameters);
//CDetector detector = new CDetector();
double izkoristek = beamtest(detector, parameters, NN, 18.5, phiMin, phiMax, nnrays, showr);
//printf("izkoristek = %.3lf\n", izkoristek);
PrintGuideHead();
PrintGuideStat(izkoristek);
DrawData(detector, parameters, 18.5, izkoristek);
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
//canvasDetector->cd();
//detector->Draw();
}
//------------------------------------------------------------------------------------------
230,7 → 262,7
parameters.setGap(step[i], 0.0, 0.0);
CDetector *detector = new CDetector(CENTER, parameters);
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
acc[i] = RandIso(detector, parameters, NN, 0, theta);
//printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
PrintGuideStat(acc[i]);
}
246,13 → 278,45
}
//------------------------------------------------------------------------------------------
 
void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
{
//int tmp3d = show_3d, tmpdata = show_data;
double show_rays = 10;
double step[steps], acc[steps];
 
//show_3d = 0; show_data = 0;
//printf(" theta | acceptance \n");
PrintGuideHead();
for(int i=0; i<=steps; i++) {
Init();
step[i] = min + i*(max - min)/(double)steps;
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
acc[i] = RandYZ(detector, parameters, NN, step[i], phi, show_rays);
//printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
PrintGuideStat(acc[i]);
delete detector;
}
//char sbuff[256];
//sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf)",
// detector->GetSiPM(), detector->GetM(), detector->GetD(),
// (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z());
DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
//show_3d = tmp3d; show_data = tmpdata;
}
 
 
//------------------------------------------------------------------------------------------
void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 5.0)
//void LGR_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15, double phi = 0)
void LGI_th(int NN = 1e4, double min = 0.0, double max = 30.0, const int steps = 15)
{
int tmp3d = show_3d, tmpdata = show_data;
double step[steps], acc[steps];
//int tmp3d = show_3d, tmpdata = show_data;
double step[steps], acc[steps];
 
show_3d = 0; show_data = 0;
//show_3d = 0; show_data = 0;
//printf(" theta | acceptance \n");
PrintGuideHead();
260,10 → 324,11
Init();
step[i] = min + i*(max - min)/(double)steps;
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
acc[i] = RandYZ(detector, parameters, NN, step[i], phi);
//acc[i] = RandYZ(detector, parameters, NN, step[i], phi, 30);
acc[i] = RandIso(detector, parameters, NN, 0, step[i], 30, show_3d);
//printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
PrintGuideStat(acc[i]);
//delete detector;
delete detector;
}
//char sbuff[256];
273,19 → 338,47
DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
show_3d = tmp3d; show_data = tmpdata;
//show_3d = tmp3d; show_data = tmpdata;
}
//------------------------------------------------------------------------------------------
 
void LGR_phi(int NN = 1e4, double maxTheta = 10.0, double maxPhi = 10.0, int steps = 10)
{
int tmp3d = show_3d, tmpdata = show_data;
show_3d = 1; show_data = 0;
//PrintGuideHead();
TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance;#theta;#phi", steps, 0, maxTheta, steps, 0, maxPhi);
double min = 0.0;
printf("\nWait, this takes a while ");
for(int i=0; i<=steps; i++) {
double theta = min + i*(maxTheta - min) / (double)steps;
for (int j=0; j<=steps; j++) {
Init();
double phi = min + j*(maxPhi - min) / (double)steps;
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
double acc = RandYZ(detector, parameters, NN, theta, phi, 30);
//PrintGuideStat(acc);
hAcceptance->Fill(i, j, acc);
//printf("Acc: %f ", acc);
delete detector;
}
printf(".");
}
printf("\n");
//DrawAcc(steps+1, step, acc, (char*)"; #theta [deg];Acceptance", min, max);
show_3d = tmp3d; show_data = tmpdata;
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1000, 1000);
cp->cd();
hAcceptance->Draw("COLZ");
}
 
// a vs. d
void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0)
{
int tmp3d = show_3d;
int tmpdata = show_data;
//char sbuff[256];
 
show_3d = 0; // don't show simulations
show_data = 1;
//show_3d = 0; // don't show simulations
//show_data = 1;
//double d = detector->GetD();
const double b = parameters.getB(); // upper side of LG
297,34 → 390,43
// Use the Fresnel eq. instead of fixed reflectivity 96%
//detector->SetFresnel(1);
printf(" d | a | Acceptance\n");
//printf(" d | a | Acceptance\n");
getParameters();
printf("Wait, this takes a while ");
for(int i=0; i<steps; i++) {
const double d = hAcceptance->GetXaxis()->GetBinCenter(i);
for(int j=0; j<steps; j++) {
Init();
const double a = hAcceptance->GetYaxis()->GetBinCenter(j);
const double M = b/a;
parameters.setGuide(a, M, d);
//const double M = b/a;
parameters.setGuide(a, b, d);
CDetector *detector = new CDetector(CENTER, parameters);
//detector->guide->setLG(y, M, x);
//Init(); exclude simulation
double acceptance = RandIso(detector, parameters, NN, theta, 0, 0);
double acceptance = RandIso(detector, parameters, NN, 0, theta, 30, show_3d);
//double acceptance = Grid(NN, theta);
//double acceptance = -1.0;
hAcceptance->Fill(d, a, acceptance);
//printf("%.2lf | %.2lf | ", x,y);
//PrintGuideStat(izkoristek);
//printf("%.2lf | %.2lf | ", d, a);
//PrintGuideStat(acceptance);
delete detector; //works fine, 50x50 grid takes ~4MB of RAM
}
printf(".");
}
printf("\n");
 
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);
TVirtualPad *pacc = cp->cd(0);
TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 1200, 1200);
cp->cd();
//TVirtualPad *pacc = cp->cd(0);
/*
pacc->SetRightMargin(0.10);
pacc->SetLeftMargin(0.10);
pacc->SetTopMargin(0.10);
pacc->SetBottomMargin(0.10);
*/
TFile *file = new TFile("acceptance.root","RECREATE");
hAcceptance->Write();
336,8 → 438,10
hAcceptance->SetTitle(";d [mm];a [mm]");
hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);
hAcceptance->Draw("COLZ");
char filename[128];
sprintf(filename,"LGI_ad%d.C", steps);
cp->SaveAs(filename);
show_3d = tmp3d; show_data = tmpdata;
}
 
//------------------------------------------------------------------------------------------
364,11 → 468,11
for(int i=0; i<=steps; i++) {
step[i] = min + i*(max - min)/(double)steps;
sipm_d = M0 - 2.0*step[i]*0.199819;//tan 11.3 deg
parameters.setGuide(sipm_d, M0/sipm_d, step[i]);
parameters.setGuide(sipm_d, M0*sipm_d, step[i]);
//printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
printf("%.1lf | %.2lf | ", step[i], sipm_d);
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
acc[i] = RandIso(detector, parameters, NN, 0, theta);
PrintGuideStat(acc[i]);
sprintf(sbuff, "d%2d.gif", i);
//c3dview->SaveAs(sbuff);
422,12 → 526,13
PrintGuideHead();
for(int i=0; i<=steps; i++) {
step[i] = min + i*(max - min)/(double)steps;
//parameters.setGuide(a, magnif, step[i]);
parameters.setGuide(a, magnif, step[i]);
CDetector *detector = new CDetector(CENTER, parameters);
//printf("sipm = %lf, M = %lf, d = %lf\n", detector->GetSiPM(), detector->GetM(), step[i]);
printf("%.1lf | ", step[i]);
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
acc[i] = RandIso(detector, parameters, NN, 0, theta);
PrintGuideStat(acc[i]);
delete detector;
}
456,10 → 561,11
PrintGuideHead();
for(int i=0; i<=steps; i++) {
step[i] = min + i*(max - min)/(double)steps;
parameters.setGuide(a, step[i], d);
//parameters.setGuide(a, step[i], d);
parameters.setGuide(a, a*step[i], d);
CDetector *detector = new CDetector(CENTER, parameters);
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
acc[i] = RandIso(detector, parameters, NN, 0, theta);
PrintGuideStat(acc[i]);
}
472,5 → 578,17
show_3d = tmp3d; show_data = tmpdata;
}
 
//------------------------------------------------------------------------------------------
 
void getParameters()
{
printf("LIGHT GUIDE\n"
" b=%f mm, d=%f mm, a=%f mm\n", parameters.getB(), parameters.getD(), parameters.getA());
printf("MATERIAL REFRACITVE INDICES\n"
" n1=%f, n2=%f, n3=%f\n", parameters.getN1(), parameters.getN2(), parameters.getN3());
printf("PLATE\n"
" ON: %d, width=%f mm\nUSES FRESNEL EQUATIONS: %d\n", parameters.getPlateOn(), parameters.getPlateWidth(), parameters.getFresnel());
return;
}
/lightguide/trunk/src/raySimulator.cpp
87,8 → 87,8
{
if(show_data) {
char sbuff[256];
sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
parameters.getA(), parameters.getM(), parameters.getD(),
sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
parameters.getA(), parameters.getLightYield()*acc, parameters.getD(),
parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
if(!only2d) {
183,6 → 183,7
theta = Pi()*theta/180.0;
if(theta < 1e-6) theta = 1e-6;
phi = phi*Pi()/180.0;
if(phi < 1e-6) phi = 1e-6;
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
195,6 → 196,9
detector->Propagate(*ray0, ray1, show_3d);
delete ray0;
delete ray1;
return (detector->GetHActive())->GetEntries() / (double)(1);
}
//-----------------------------------------------------------------------------
224,16 → 228,26
if(i == (NN/2))
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
detector->Propagate(*ray0, ray1, 0);
delete ray0;
delete ray1;
}
}
return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
return acceptance;
}
//-----------------------------------------------------------------------------
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
// vsi pod kotom (theta, phi)
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phi, int show_rays)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
theta = theta*3.14159265358979312/180.0;
272,13 → 286,25
//delete ray0;
//delete ray1;
}
return (detector->GetHActive())->GetEntries() / (double)NN;
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
//return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/
acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
return acceptance;
}
//-----------------------------------------------------------------------------
// zarki, izotropno porazdeljeni znotraj kota theta
// = nakljucni vstopni polozaj in kot
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
// NN - number of rays to be simulated with angles theta distributed uniformly
// in cos theta and phi uniformly:
// theta [0, theta]
// phi [0,360]
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
{
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
double pi = 3.14159265358979312;
292,6 → 318,7
double rphi, rtheta;
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
double theta_min_rad = TMath::Cos(theta);
double theta_max_rad = TMath::Cos(thetaMin);
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
322,7 → 349,9
rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
rphi = rand.Uniform(0.0, 2.0*pi);
rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
//double theta_max_rad = TMath::Cos(17.0*pi/180.0);
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
rl = TMath::Cos(rtheta);
rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
345,8 → 374,8
detector->Propagate(*ray0, ray1, 0);
}
//delete ray0;
//delete ray1;
delete ray0;
delete ray1;
}
if(show_rand) {
363,8 → 392,113
c2rand->cd(6); hn->Draw();
}
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
return acceptance;
}
 
// Beamtest distribution
// with fixed theta and phi in interval phiMin, phiMax
double beamtest(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
{
double pi = 3.14159265358979312;
theta *= pi/180.0;
if(theta < MARGIN) theta = MARGIN;
phiMin *= pi/180.0;
phiMax *= pi/180.0;
double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
TDatime now;
TRandom rand;
rand.SetSeed(now.Get());
double rx, ry, rz, rl, rm, rn;
double rphi;
//double rtheta;
//double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
//double theta_min_rad = TMath::Cos(theta);
 
TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
hphi = new TH1F("hphi", "hphi", 100, -pi, pi);
htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
hl = new TH1F("hl", "hl", 100, 0.0, 1.0);
hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
hm = new TH1F("hm", "hm", 100, 0.0, 1.0);
hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
//detector->Init();
if (show_3d) detector->Draw(draw_width);
 
double b = parameters.getB();
double offset = (parameters._plateOn ? parameters._plateWidth : 0);
for(int i = 0; i < NN; i++) {
rx = CENTER.x() - offset;
ry = rand.Uniform(-b/2.0, b/2.0);
rz = rand.Uniform(-b/2.0, b/2.0);
rphi = rand.Uniform(phiMin, phiMax);
//double theta_max_rad = TMath::Cos(17.0*pi/180.0);
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
//rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
rl = TMath::Cos(theta);
rm = TMath::Sin(theta)*TMath::Sin(rphi);
rn = TMath::Sin(theta)*TMath::Cos(rphi);
if(show_rand) {
//htheta->Fill(rtheta);
//hcostheta->Fill( TMath::Cos(rtheta) );
hphi->Fill(rphi);
hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
}
CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
CRay *ray1 = new CRay;
if(i < show_rays) {
detector->Propagate(*ray0, ray1, show_3d);
}
else {
detector->Propagate(*ray0, ray1, 0);
}
delete ray0;
delete ray1;
}
if(show_rand) {
TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
else c2rand->Clear();
c2rand->Divide(2,3);
c2rand->cd(1); hphi->Draw();
c2rand->cd(3); htheta->Draw();
c2rand->cd(5); hcostheta->Draw();
c2rand->cd(2); hl->Draw();
c2rand->cd(4); hm->Draw();
c2rand->cd(6); hn->Draw();
}
double acceptance = 0.0;
/*
if( !(parameters.getPlateOn()) )
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
else
acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
*/
acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
return acceptance;
}
 
/lightguide/trunk/src/guide.cpp
247,9 → 247,9
 
if(dbg)
{
//printf("angle_ve1 = %lf\n", angle_ve1*DEGREE);
//printf("angle_ve2 = %lf\n", angle_ve2*DEGREE);
//printf("angle_sum = %lf\n", (angle_ve1 + angle_ve2)*DEGREE);
printf("angle_ve1 = %lf\n", angle_ve1*DEGREE);
printf("angle_ve2 = %lf\n", angle_ve2*DEGREE);
printf("angle_sum = %lf\n", (angle_ve1 + angle_ve2)*DEGREE);
printf(" angle_r = %lf\n", angle*DEGREE);
}
 
397,23 → 397,32
// --------------- Fresnel ----------------------------------------------------
// R_f = a_te * R_te + a_tm * R_tm
// e - electrical/perependicular
// m - magnetic polarization/parallel
double r_te=0;
double r_tm=0;
double R_te=0;
double R_tm=0;
double R_te=0; // s reflection coefficient
double R_tm=0; // p reflection coefficient
double R_f = 0.0;
double a_te = 0.0;
double a_tm = 0.0;
TVector3 v_te; // polarization perpendicular to the plane of incidence
TVector3 v_tm; // inbound polarization parallel with the plane of incidence
double a_te = 0.0; // s-wave amplitude, cos Alpha
double a_tm = 0.0; // p-wave amplitude, sin Alpha
TVector3 v_te; // unit s-polarization vector
TVector3 v_tm; // unit p-polarization vector
TVector3 v_tm_t;// transmited polarization parallel with the plane of incidence
TVector3 pol_t = in.GetP(); // transmited polarization
int sign_n; // sign of normal direction vs. inbound ray
double cosTN; // debug
if(fresnel) {
v_te = n.Cross(in.GetN()); v_te = v_te.Unit();
v_tm = -v_te.Cross(in.GetN()); v_tm = v_tm.Unit();
if(fresnel) {
// p-polarization unit vector v_te
// is in the plane orthogonal to the plane of incidence
// defined as the plane spanned by
// incident surface vector n and wave vector k
// k in this notation is in.GetN()
v_te = n.Cross(in.GetN());
v_te = v_te.Unit();
v_tm = -v_te.Cross(in.GetN());
v_tm = v_tm.Unit();
if(dbg) {
printf(" v_te = "); printv(v_te);
printf(" v_tm = "); printv(v_tm);
428,7 → 437,7
}
// ----------------------------------------------------------------------------
if(type == 3) type = SURF_REFRA; //SURF_TOTAL -> SURF_REFRA
if(type == SURF_TOTAL) type = SURF_REFRA;
switch(type){
// ----------------------------------------------------------------------------
// --------------- refraction from n1 to n2 -----------------------------------
452,9 → 461,8
p_ref = rand.Uniform(0.0, 1.0);
if (dbg) printf(" reflection probability = %f\n", p_ref);
// Total reflection
/*
if( (cosTi <= cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
// If n1>n2 and theta>thetaCritical, total reflection
if( (cosTi < cosTtotal) && (p_ref < reflection) ) { // totalni odboj z verjetnostjo "reflection"
if(dbg) printf(" TOTAL\n");
transmit = in.GetN() + sign_n*2*cosTi*n;
467,7 → 475,7
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
} else { */
} else {
// reflection or refraction according to Fresnel equations
if(dbg) printf(" REFRACTION\n");
if(dbg) printf(" N1_N2(sign_n) = %lf\n", N1_N2(sign_n));
488,6 → 496,7
if(dbg) printf(" r_te = %lf, r_tm = %lf\n", r_te, r_tm);
// transmited polarization
v_tm_t = -v_te.Cross(transmit);
v_tm_t = v_tm_t.Unit();
pol_t = a_te * (1.0 - TMath::Abs(r_te)) * v_te + a_tm * (1.0 - TMath::Abs(r_tm)) * v_tm_t;
496,13 → 505,14
printf(" v_tm_t = "); printv(v_tm_t);
printf(" pol_t = "); printv(pol_t);
}
 
// Fresnel coefficients
R_te = TMath::Power(r_te, 2);
R_tm = TMath::Power(r_tm, 2);
R_f = a_te*a_te*R_te + a_tm*a_tm*R_tm;
if (dbg) printf(" R_te = %lf, R_tm = %lf, R_f = %lf\n", R_te, R_tm, R_f);
//}
}
if(p_ref >= R_f) { // se lomi
if (dbg) printf(" SURFACE REFRACTED. Return.\n");
573,7 → 583,8
{
double t;
 
TDatime now; rand.SetSeed(now.Get());
TDatime now;
rand.SetSeed(now.Get());
 
center = center0;
double b = parameters.getB();
582,33 → 593,41
n1 = parameters.getN1();
n2 = parameters.getN2();
// if PlateOn, then n0 = n3 (optical grease), else = n1 (air)
double n0 = (parameters.getPlateOn() ? parameters.getN3(): n1);
//double n0 = (parameters.getPlateOn() ? parameters.getN3(): n1);
double n0 = (parameters.getPlateOn() ? n2 : n1);
n3 = parameters.getN3();
_r = c_reflectivity;
int fresnel = parameters.getFresnel();
 
// light guide edges
t = b/2.0;
vodnik_edge[0].SetXYZ(0.0, t,-t); vodnik_edge[1].SetXYZ(0.0, t, t);
vodnik_edge[2].SetXYZ(0.0,-t, t); vodnik_edge[3].SetXYZ(0.0,-t,-t);
vodnik_edge[0].SetXYZ(0.0, t,-t);
vodnik_edge[1].SetXYZ(0.0, t, t);
vodnik_edge[2].SetXYZ(0.0,-t, t);
vodnik_edge[3].SetXYZ(0.0,-t,-t);
t = a/2.0;
vodnik_edge[4].SetXYZ(_d, t,-t); vodnik_edge[5].SetXYZ(_d, t, t);
vodnik_edge[6].SetXYZ(_d,-t, t); vodnik_edge[7].SetXYZ(_d,-t,-t);
vodnik_edge[4].SetXYZ(_d, t,-t);
vodnik_edge[5].SetXYZ(_d, t, t);
vodnik_edge[6].SetXYZ(_d,-t, t);
vodnik_edge[7].SetXYZ(_d,-t,-t);
for(int i = 0; i<8; i++) vodnik_edge[i] += center;
// light guide surfaces
s_side[0] = new CSurface(SURF_REFRA, vodnik_edge, n0, n2, _r);
s_side[0]->FlipN();
s_side[1] = new CSurface(SURF_REFRA, vodnik_edge[3], vodnik_edge[2], vodnik_edge[6], vodnik_edge[7], n2, n1, _r);
s_side[2] = new CSurface(SURF_REFRA, vodnik_edge[2], vodnik_edge[1], vodnik_edge[5], vodnik_edge[6], n2, n1, _r);
s_side[3] = new CSurface(SURF_REFRA, vodnik_edge[1], vodnik_edge[0], vodnik_edge[4], vodnik_edge[5], n2, n1, _r);
s_side[4] = new CSurface(SURF_REFRA, vodnik_edge[0], vodnik_edge[3], vodnik_edge[7], vodnik_edge[4], n2, n1, _r);
s_side[5] = new CSurface(SURF_REFRA, &vodnik_edge[4], n2, n3, _r);
s_side[5] = new CSurface(SURF_REFRA, &vodnik_edge[4], n2, n3, _r); // n3 - ref ind at the exit, grease, air, epoxy
s_side[5]->FlipN();
if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1);
// statistics histograms
hfate = (TH1F*)gROOT->FindObject("hfate"); if(hfate) delete hfate;
hfate = new TH1F("hfate", "Ray fate", 8, -3.5, 4.5);
(hfate->GetXaxis())->SetBinLabel(1, "Back Ref");
670,7 → 689,7
} else {
if (dbg) printf(" GUIDE: ray entered\n");
points[0] = ray1.GetR();
hfate->Fill(2); // enter
hfate->Fill(enter); // enter
hin->Fill(vec1.y(), vec1.z());
if (dbg) printf(" GUIDE: n_odb = %d\n", n_odb);
890,14 → 909,17
//guide = new CVodnik(center, SiPM, M, d, type_in, type_side, type_out, n1, n2, n3, reflectivity, fresnel, absorption, A);
double b = parameters.getB();
double n1 = parameters.getN1();
double n2 = parameters.getN2();
//double n3 = parameters.getN3();
double reflectivity = c_reflectivity; // for faster simulation, not using Fresnel eqs.
//double n1 = parameters.getN1();
//double n2 = parameters.getN2();
double n3 = parameters.getN3();
double reflectivity = c_reflectivity;
double x_gap = parameters.getGap().X();
double y_gap = parameters.getGap().Y();
double z_gap = parameters.getGap().Z();
// additional glass between at top of SiPM
// example: epoxy n=1.60
double n4 = 1.60;
TVector3 plane_v[4];
int nBins = nch + 1;
double p_size = b/2.0;
905,8 → 927,10
plane_v[1].SetXYZ(x_offset+d+glass_d, y_gap + p_size, z_gap + p_size);
plane_v[2].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap + p_size);
plane_v[3].SetXYZ(x_offset+d+glass_d, y_gap - p_size, z_gap - p_size);
glass = new CSurface(SURF_REFRA, plane_v, n2, n1, reflectivity); glass->FlipN();
glass = new CSurface(SURF_REFRA, plane_v, n3, n4, reflectivity);
glass->FlipN();
// additional circular glass between LG and SiPM
glass_circle = new CPlaneR(TVector3(x_offset+d+glass_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), b);
hglass = (TH2F*)gROOT->FindObject("hglass"); if(hglass) delete hglass;
914,7 → 938,7
nBins, y_gap - p_size, y_gap + p_size,
nBins, z_gap - p_size, z_gap + p_size);
// SiPM active surface
p_size = parameters.getActive()/2.0;
//cout<<"SiPM active length "<<detectorActive<<endl;
//p_size = 1.0/2.0;
934,7 → 958,7
hlaser = (TH2F*)gROOT->FindObject("hlaser"); if(hlaser) delete hlaser;
hlaser = new TH2F("hlaser", ";x [mm]; y [mm]", nBins, -p_size+offsetY, p_size+offsetY, nBins, -p_size+offsetZ, p_size+offsetZ);
// collection surface in SiPM plane
p_size = 1.4*b/2.0;
plane_v[0].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap - p_size);
plane_v[1].SetXYZ(x_offset+d+x_gap, y_gap + p_size, z_gap + p_size);
982,7 → 1006,8
//CRay *ray0 = new CRay; ray0->Set(in.GetR(), in.GetN()); ray0->SetColor(col_in);
CRay *rayin = new CRay(in);
rayin->SetColor(col_in);
CRay *rayout = new CRay;
CRay *rayout = new CRay(in);
rayout->SetColor(col_in);
 
const int max_n_points = guide->GetMAXODB() + 2;
TVector3 pointsPlate[max_n_points];
999,14 → 1024,16
if(_plateOn) {
fatePlate = plate->propagateRay(*rayin, rayout, &nPointsPlate, pointsPlate);
if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
if(draw) {
if(fatePlate == missed) {
rayout->SetColor(col_in);
rayout->DrawS(center.x() - _plateWidth, -10.0);
}
else {
else if(fatePlate == backreflected){
if (dbg) printf("Backreflected at plate!\n");
}
else {
int p_i;
for(p_i = 0; p_i < nPointsPlate-1; p_i++) {
line3d->SetPoint(0, pointsPlate[p_i].x(), pointsPlate[p_i].y(), pointsPlate[p_i].z());
1013,24 → 1040,30
line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z());
line3d->DrawClone();
}
if(fatePlate != noreflection) {
//rayout->SetColor(7);
rayout->DrawS(pointsPlate[p_i].x(), -0.1);
rayout->DrawS(pointsPlate[p_i].x(), -0.1);
if(fatePlate == noreflection) { // lost on plate side
rayout->SetColor(col_out);
rayout->DrawS(pointsPlate[p_i].x(), 10.0);
}
}
}
}
}
if(! (fatePlate == hitExit or fatePlate == refracted) ) {
if (dbg) printf("CDetector::propagate Simulated ray missed the entry surface!\n");
guide->GetHFate()->Fill(rays);
if (dbg)printf("CDetector::propagate Simulated ray missed the entry surface!\n");
if (fatePlate == backreflected)
guide->GetHFate()->Fill(fatePlate); // reflected back
else
guide->GetHFate()->Fill(noreflection); //lost on plate side
return fatePlate;
}
// missing: if refracted at plate sides
//if (fatePlate == refracted) return fatePlate;
//Ray hits light guide
histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
 
}
else {
rayout = rayin;
//rayout = rayin;
if(draw) rayout->DrawS(center.x(), -10.0);
}
1045,7 → 1078,7
int fate_glass;
CRay *ray0 = new CRay(*rayout);
// delete rayout; -> creates dangling reference when tries to delete ray0!
delete rayin;
//delete rayin; -> delete rayout!
CRay *ray1 = new CRay;
fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
1084,6 → 1117,10
if(! (fate == hitExit or fate == refracted) ) {
if (dbg) printf("Detector: fate != hit, refracted\n");
*out = *ray1;
delete ray0;
delete ray1;
delete rayout;
delete rayin;
return fate;
}
} else {
1095,18 → 1132,6
}
}
/*
TVector3 pres_wind;
fate = window_circle->TestIntersection(&pres_wind, *ray1);
if(fate == 1) {
hwindow->Fill(pres_wind.y(), pres_wind.z());
if(!guide_on) {
window->PropagateRay(*ray0, ray1, &presecisce);
if(draw) ray1->Draw(center.x() + window_d, center.x() + glass_d);
*ray0 = *ray1;
}
*/
fate = missed; // zgresil aktivno povrsino
if(glass_on) {
*ray0 = *ray1; ray1->SetColor(col_rgla);
1143,6 → 1168,7
delete ray0;
delete ray1;
delete rayout;
delete rayin;
return fate;
}
//-----------------------------------------------------------------------------
1171,27 → 1197,30
const double b = parameters.getB();
const double n1 = parameters.getN1();
const double n2 = parameters.getN2();
const double n3 = parameters.getN3();
const double reflectivity = c_reflectivity;
const double t = b/2.;
const double plateWidth = parameters.getPlateWidth();
center.SetX( CENTER.X() - plateWidth );
plate_edge[0].SetXYZ(0.0, t,-t); plate_edge[1].SetXYZ(0.0, t, t);
plate_edge[2].SetXYZ(0.0,-t, t); plate_edge[3].SetXYZ(0.0,-t,-t);
plate_edge[4].SetXYZ(plateWidth, t,-t); plate_edge[5].SetXYZ(plateWidth, t, t);
plate_edge[6].SetXYZ(plateWidth,-t, t); plate_edge[7].SetXYZ(plateWidth,-t,-t);
plate_edge[0].SetXYZ(0.0, t,-t);
plate_edge[1].SetXYZ(0.0, t, t);
plate_edge[2].SetXYZ(0.0,-t, t);
plate_edge[3].SetXYZ(0.0,-t,-t);
plate_edge[4].SetXYZ(plateWidth, t,-t);
plate_edge[5].SetXYZ(plateWidth, t, t);
plate_edge[6].SetXYZ(plateWidth,-t, t);
plate_edge[7].SetXYZ(plateWidth,-t,-t);
for(int i = 0; i<8; i++) plate_edge[i] += center;
sides[0] = new CSurface(SURF_REFRA, plate_edge, n1, n2, reflectivity);
sides[0] = new CSurface(SURF_REFRA, plate_edge, n1, n2, c_reflectivity);
sides[0]->FlipN();
sides[1] = new CSurface(SURF_REFRA, plate_edge[3], plate_edge[2], plate_edge[6], plate_edge[7], n2, n2, reflectivity);
sides[2] = new CSurface(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, reflectivity);
sides[3] = new CSurface(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, reflectivity);
sides[4] = new CSurface(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, reflectivity);
sides[1] = new CSurface(SURF_REFRA, plate_edge[3], plate_edge[2], plate_edge[6], plate_edge[7], n2, n2, c_reflectivity);
sides[2] = new CSurface(SURF_REFRA, plate_edge[2], plate_edge[1], plate_edge[5], plate_edge[6], n2, n2, c_reflectivity);
sides[3] = new CSurface(SURF_REFRA, plate_edge[1], plate_edge[0], plate_edge[4], plate_edge[5], n2, n2, c_reflectivity);
sides[4] = new CSurface(SURF_REFRA, plate_edge[0], plate_edge[3], plate_edge[7], plate_edge[4], n2, n2, c_reflectivity);
sides[5] = new CSurface(SURF_REFRA, &plate_edge[4], n2, n3, reflectivity);
sides[5] = new CSurface(SURF_REFRA, &plate_edge[4], n2, n2, c_reflectivity);
sides[5]->FlipN();
for(int i=0; i<6; i++) sides[i]->SetFresnel(1);
1237,7 → 1266,7
fate = backreflected;
} else {
points[0] = ray1.GetR();
//hfate->Fill(2);
//hfate->Fill(enter);
//hin->Fill(vec1.y(), vec1.z());
while (n_odb++ < MAX_REFLECTIONS) {
ray0 = ray1;
1264,7 → 1293,7
break;
}
if(propagation == 1) {
fate = refracted; //at side?
fate = noreflection; //at side
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;