Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 69 → Rev 70

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