Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 25 → Rev 54

/lightguide/trunk/src/RTUtil.cpp
13,10 → 13,10
style->SetStatColor(0);
style->SetPalette(1, 0);
style->SetPadTopMargin(0.13);
style->SetPadBottomMargin(0.13);
style->SetPadLeftMargin(0.13);
style->SetPadRightMargin(0.2);
style->SetPadTopMargin(0.10);
style->SetPadBottomMargin(0.10);
style->SetPadLeftMargin(0.10);
style->SetPadRightMargin(0.12);
style->SetTitleOffset(1.3, "y");
}
/lightguide/trunk/src/userFunctions.cpp
171,7 → 171,7
Init();
double izkoristek = Single(detector, parameters, TVector3(xh,yh,zh),theta, phi);
PrintGuideHead();
PrintGuideStat(detector, 1);
PrintGuideStat(izkoristek);
DrawData(detector, parameters, theta, izkoristek);
}
//------------------------------------------------------------------------------------------
183,7 → 183,7
double izkoristek = Grid(detector, parameters, NN, theta);
//printf("izkoristek = %.3lf\n", izkoristek);
PrintGuideHead();
PrintGuideStat(detector, (NN+1)*(NN+1));
PrintGuideStat(izkoristek);
DrawData(detector, parameters, theta, izkoristek);
}
//------------------------------------------------------------------------------------------
194,7 → 194,7
Init();
double izkoristek = RandYZ(detector, parameters, NN, theta, phi);
//printf("izkoristek = %.3lf\n", izkoristek);
PrintGuideHead(); PrintGuideStat(detector, NN);
PrintGuideHead(); PrintGuideStat(izkoristek);
DrawData(detector, parameters, theta, izkoristek);
}
//------------------------------------------------------------------------------------------
206,12 → 206,11
//CDetector detector = new CDetector();
double izkoristek = RandIso(detector, parameters, NN, theta, nnrays, showr);
//printf("izkoristek = %.3lf\n", izkoristek);
PrintGuideHead(); PrintGuideStat(detector, NN);
PrintGuideHead(); PrintGuideStat(izkoristek);
DrawData(detector, parameters, theta, izkoristek);
//TCanvas *canvasDetector = new TCanvas("canvasDetector","canvasDetector",500,500);
//canvasDetector->cd();
//detector->Draw();
//delete detector;
}
 
233,7 → 232,7
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
//printf("%6.2lf |%6.1lf\n", step[i], acc[i]*100.0);
PrintGuideStat(detector, NN);
PrintGuideStat(acc[i]);
}
//char sbuff[256];
263,8 → 262,8
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
acc[i] = RandYZ(detector, parameters, NN, step[i], phi);
//printf("%6.1lf |%6.1lf\n", step[i], acc[i]*100.0);
PrintGuideStat(detector, NN);
delete detector;
PrintGuideStat(acc[i]);
//delete detector;
}
//char sbuff[256];
313,8 → 312,8
//double acceptance = -1.0;
hAcceptance->Fill(d, a, acceptance);
//printf("%.2lf | %.2lf | ", x,y);
//PrintGuideStat(detector, NN);
//delete detector; //works fine, 50x50 grid takes ~4MB of RAM
//PrintGuideStat(izkoristek);
delete detector; //works fine, 50x50 grid takes ~4MB of RAM
}
}
330,7 → 329,7
TFile *file = new TFile("acceptance.root","RECREATE");
hAcceptance->Write();
file->Close();
delete file;
//delete file;
//hAcceptance->SetContour(100);
//gStyle->SetPalette(1,0);
370,7 → 369,7
printf("%.1lf | %.2lf | ", step[i], sipm_d);
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
PrintGuideStat(detector, NN);
PrintGuideStat(acc[i]);
sprintf(sbuff, "d%2d.gif", i);
//c3dview->SaveAs(sbuff);
}
429,7 → 428,7
printf("%.1lf | ", step[i]);
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
PrintGuideStat(detector, NN);
PrintGuideStat(acc[i]);
delete detector;
}
461,7 → 460,7
CDetector *detector = new CDetector(CENTER, parameters);
Init();
acc[i] = RandIso(detector, parameters, NN, theta);
PrintGuideStat(detector, NN);
PrintGuideStat(acc[i]);
}
//char sbuff[256];
/lightguide/trunk/src/raySimulator.cpp
153,7 → 153,7
printf("Acceptance\n");
}
//-----------------------------------------------------------------------------
void PrintGuideStat(CDetector *detector, double n_in = 1)
void PrintGuideStat(double acceptance)
{
/*int n_active = (detector->GetHActive())->GetEntries();
int n_enter = (detector->GetLG())->GetEnteranceHits();
168,8 → 168,9
printf("%5.1lf\n", 100.0*izkoristek);*/
double n_active = (detector->GetHActive())->GetEntries();
double r_acc = 100.0 * n_active / n_in;
//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);
}
//-----------------------------------------------------------------------------
268,8 → 269,8
detector->Propagate(*ray0, ray1, show_3d);
else
detector->Propagate(*ray0, ray1, 0);
delete ray0;
delete ray1;
//delete ray0;
//delete ray1;
}
return (detector->GetHActive())->GetEntries() / (double)NN;
344,8 → 345,8
detector->Propagate(*ray0, ray1, 0);
}
delete ray0;
delete ray1;
//delete ray0;
//delete ray1;
}
if(show_rand) {
/lightguide/trunk/src/guide.cpp
2,7 → 2,6
 
#include <iostream>
 
 
// vector output shortcut
void printv(TVector3 v)
{
448,9 → 447,13
cosTtotal = 0.0;
if(dbg) printf(" cosTtotal = %lf (Ttotal = %lf)\n", cosTtotal, TMath::ACos(cosTtotal)*DEGREE);
// reflection dependance on polarization missing
// reflection hardcoded to 0.96
p_ref = rand.Uniform(0.0, 1.0);
if (dbg) printf(" reflection probability = %f\n", p_ref);
p_ref = rand.Uniform(0.0, 1.0);
// 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;
463,8 → 466,9
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return 2;
} else { // ni tot. odboja
return REFLECTION;
} 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));
cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));
478,7 → 482,7
printf(" cosTN = %lf (TN = %lf) (Abs(TN) = %lf)\n", cosTN, TMath::ACos(cosTN)*DEGREE, TMath::ACos(TMath::Abs(cosTN))*DEGREE);
}
//if(cosTi<=cosTtotal) cosTt = TMath::Sqrt(1 - TMath::Power(N1_N2(sign_n), 2)*(1 - TMath::Power(cosTi, 2)));
if(fresnel) {
//if(fresnel) {
r_te = (n1*cosTi - n2*cosTt)/(n1*cosTi + n2*cosTt); // transverse
r_tm = (n2*cosTi - n1*cosTt)/(n1*cosTt + n2*cosTi); // paralel
497,23 → 501,26
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);
}
//p_ref = rand.Uniform(0.0, 1.0);
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");
out->Set(intersect, transmit);
out->SetPolarization(pol_t);
return 1;
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;
out->Set(intersect, transmit);
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return 2;
}
}
return REFLECTION;
}
//}
break;
// ----------------------------------------------------------------------------
// --------------- reflection at "reflection" probability ---------------------
// ----------------------------------------------------------------------------
523,11 → 530,11
cosTi = in.GetN() * n;
transmit = in.GetN() - 2*cosTi*n;
out->Set(intersect, transmit);
return 2; //sdhfvjhsdbfjhsdbcvjhsb
return REFLECTION; //sdhfvjhsdbfjhsdbcvjhsb
} else { // se ne odbije
transmit = in.GetN();
out->Set(intersect, transmit);
return -2;
return ABSORBED;
}
break;
542,12 → 549,12
} else { // ni tot. odboja
transmit = in.GetN();
out->Set(intersect, transmit);
return -2;
return ABSORBED;
}
} else { // se ne odbije
transmit = in.GetN();
out->Set(intersect, transmit);
return -2;
return ABSORBED;
}
break;
 
556,7 → 563,7
break;
}
 
return 1;
return REFRACTION;
}
//=================================================================================
 
605,8 → 612,8
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");
(hfate->GetXaxis())->SetBinLabel(2, "Ref Loss");
(hfate->GetXaxis())->SetBinLabel(3, "No Tot Refl");
(hfate->GetXaxis())->SetBinLabel(2, "No Ref");
(hfate->GetXaxis())->SetBinLabel(3, "Refrac");
(hfate->GetXaxis())->SetBinLabel(4, "LG Miss");
(hfate->GetXaxis())->SetBinLabel(5, "Exit");
(hfate->GetXaxis())->SetBinLabel(6, "Enter");
614,10 → 621,10
(hfate->GetXaxis())->SetBinLabel(8, "Absorb");
hnodb_all = (TH1F*)gROOT->FindObject("hnodb_all"); if(hnodb_all) delete hnodb_all;
hnodb_all = new TH1F("hnodb_all", "N reflected", VODNIK_MAX_ODB, -0.5, VODNIK_MAX_ODB-0.5);
hnodb_all = new TH1F("hnodb_all", "N reflected", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
hnodb_exit = (TH1F*)gROOT->FindObject("hnodb_exit"); if(hnodb_exit) delete hnodb_exit;
hnodb_exit = new TH1F("hnodb_exit", "N reflected and exit", VODNIK_MAX_ODB, -0.5, VODNIK_MAX_ODB-0.5);
hnodb_exit = new TH1F("hnodb_exit", "N reflected and exit", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
int nBins = nch + 1;
hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin;
644,7 → 651,6
CRay ray0;
CRay ray1;
TVector3 vec0, vec1;
//int fate = 11;
int inters_i = 0;
ray0 = in;
651,11 → 657,16
int n_odb = 0;
int last_hit = 0;
int propagation = 0;
if( !(s_side[0]->PropagateRay(ray0, &ray1, &vec1)) ) {
int result = s_side[0]->PropagateRay(ray0, &ray1, &vec1);
if( !(result) ) {
// ce -NI- presecisca z vstopno
if (dbg) printf(" GUIDE: missed the light guide\n");
fate = missed;
hfate->Fill(0);
//hfate->Fill(0);
} else if(result == REFLECTION) {
if (dbg) printf(" REFLECTED on the entry surface!\n");
fate = backreflected;
//hfate->Fill(-3);
} else {
if (dbg) printf(" GUIDE: ray entered\n");
points[0] = ray1.GetR();
662,7 → 673,8
hfate->Fill(2); // enter
hin->Fill(vec1.y(), vec1.z());
if (dbg) printf(" GUIDE: n_odb = %d\n", n_odb);
while (n_odb++ < VODNIK_MAX_ODB) {
while (n_odb++ < MAX_REFLECTIONS) {
if (dbg) printf(" GUIDE: Boundary test: %d\n",n_odb);
ray0 = ray1;
vec0 = vec1;
680,12 → 692,27
points[n_odb] = vec1;
if(inters_i == 0) {
fate = backreflected;
hfate->Fill(backreflected);
//hfate->Fill(backreflected);
break;
} // backreflection
// the passage is possible, test propagation
propagation = s_side[inters_i]->PropagateRay(ray0, &ray1, &vec1);
if (dbg) printf(" GUIDE: surface = %d, propagation = %d\n", inters_i, propagation);
if(propagation == REFRACTION) {
fate = refracted;
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
} // no total reflection when should be
if(propagation == ABSORBED) {
fate = noreflection;
break;
} //refraction due to finite reflectivity
if(inters_i == 5) { // successfull exit
// check on which side the vector is?
TVector3 ray = ray1.GetN();
701,7 → 728,7
ray0 = ray1;
break;
}
fate = hit;
fate = hitExit;
hout->Fill(vec1.y(), vec1.z());
hnodb_exit->Fill(n_odb-1);
n_odb++;
709,18 → 736,6
ray0 = ray1;
break;
}
if(propagation == REFRACTION) {
fate = refracted;
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
} // no total reflection when should be
if(propagation == -2) {
fate = noreflection;
break;
} //refraction due to finite reflectivity
last_hit = inters_i;
}
}
743,7 → 758,7
//--- material absorption ---
hfate->Fill(fate);
hfate->Fill(3);
hfate->Fill(rays);
hnodb_all->Fill(n_odb-2);
*n_points = n_odb+1;
*out = ray0;
858,7 → 873,9
guide(new Guide(center0, parameters)),
plate(new Plate(parameters)),
_plateWidth(parameters.getPlateWidth()),
_plateOn(parameters.getPlateOn())
_plateOn(parameters.getPlateOn()),
offsetY(parameters.getOffsetY()),
offsetZ(parameters.getOffsetZ())
{
// };
908,17 → 925,17
active = new CPlane4(plane_v);
hactive = (TH2F*)gROOT->FindObject("hactive"); if(hactive) delete hactive;
hactive = new TH2F("hactive", "Active area hits", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size);
//hactive = new TH2F("hactive", "Active area hits", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size);
hactive = new TH2F("hactive", "Active area hits", nBins, y_gap - p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ);
p_size = b/2.0;
//p_size = 2.5;
//p_size = M*0.6;
hlaser = (TH2F*)gROOT->FindObject("hlaser"); if(hlaser) delete hlaser;
hlaser = new TH2F("hlaser", "Hits at LG entrance", nBins, -p_size, p_size, nBins, -p_size, p_size);
hlaser = new TH2F("hlaser", ";x [mm]; y [mm]", nBins, -p_size+offsetY, p_size+offsetY, nBins, -p_size+offsetZ, p_size+offsetZ);
p_size = 2*b/2.0;
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);
plane_v[2].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap + p_size);
926,9 → 943,9
detector = new CPlane4(plane_v);
hdetector = (TH2F*)gROOT->FindObject("hdetector"); if(hdetector) delete hdetector;
hdetector = new TH2F("hdetector", "Hits detector plane", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size);
//hdetector = new TH2F("hdetector", "Hits detector plane", nBins, y_gap - p_size, y_gap + p_size, nBins, z_gap - p_size, z_gap + p_size);
hdetector = new TH2F("hdetector", ";x [mm]; y [mm]", nBins, y_gap-p_size + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ);
/*
window_circle = new CPlaneR(TVector3(x_offset+d+window_d, y_gap, z_gap), TVector3(-1.0, 0.0, 0.0), window_R);
946,10 → 963,12
histoPlate = (TH2F*)gROOT->FindObject("histoPlate"); if(histoPlate) delete histoPlate;
histoPlate = new TH2F("histoPlate", "Hits on glass plate", nBins, -p_size, +p_size, nBins, -p_size, +p_size);
}
 
//-----------------------------------------------------------------------------
// vrne 1 ce je zadel aktvino povrsino
// vrne <1 ce jo zgresi
int CDetector::Propagate(CRay in, CRay *out, int draw)
 
// Sledi zarku skozi vodnik. Vrne:
// 0, ce zgresi vstopno ploskev MISSED
// 1, ce zadane izstopno ploskev HIT
968,7 → 987,7
const int max_n_points = guide->GetMAXODB() + 2;
TVector3 pointsPlate[max_n_points];
//TVector3 intersection;
int fatePlate;
Fate fatePlate;
int nPointsPlate;
TPolyLine3D *line3d = new TPolyLine3D(2);
line3d->SetLineWidth(1);
983,7 → 1002,7
if(draw) rayin->DrawS(center.x()- _plateWidth, -10.0);
if(draw) {
if(fatePlate == 0) {
if(fatePlate == missed) {
rayout->SetColor(col_in);
rayout->DrawS(center.x() - _plateWidth, -10.0);
}
994,7 → 1013,7
line3d->SetPoint(1, pointsPlate[p_i+1].x(), pointsPlate[p_i+1].y(), pointsPlate[p_i+1].z());
line3d->DrawClone();
}
if(fatePlate != -2) {
if(fatePlate != noreflection) {
//rayout->SetColor(7);
rayout->DrawS(pointsPlate[p_i].x(), -0.1);
}
1001,11 → 1020,13
}
}
if(! (fatePlate == 1 or fatePlate == -1) ) {
if(! (fatePlate == hitExit or fatePlate == refracted) ) {
cout<<"CDetector::propagate Simulated ray missed the entry surface!"<<endl;
if (dbg) printf("CDetector::propagate Simulated ray missed the entry surface!\n");
return fatePlate;
}
// missing: if refracted at plate sides
//if (fatePlate == refracted) return fatePlate;
histoPlate->Fill(pointsPlate[0].y(), pointsPlate[0].z()); // entry point
}
else {
1019,15 → 1040,18
//const int max_n_points = guide->GetMAXODB() + 2;
TVector3 points[max_n_points];
TVector3 presecisce;
//int fate;
 
int n_points;
int fate_glass;
CRay *ray0 = new CRay(*rayout);
// delete rayout; -> creates dangling reference when tries to delete ray0!
//delete rayin;
delete rayin;
CRay *ray1 = new CRay;
fate = guide->PropagateRay(*ray0, ray1, &n_points, points);
if (dbg) {
if (fate == backreflected) printf("DETECTOR::backreflected\n");
}
line3d->SetLineColor(col_lg);
int p_i;
1057,8 → 1081,8
}
if(! (fate == hit or fate == refracted) ) {
if (dbg) printf("Detector: fate != hit, refracted\n");
if(! (fate == hitExit or fate == refracted) ) {
if (dbg) printf("Detector: fate != hit, refracted\n");
*out = *ray1;
return fate;
}
1091,8 → 1115,8
hglass->Fill(presecisce.y(), presecisce.z());
if(draw) ray1->DrawS(presecisce.x(), 10.0);
if(active->TestIntersection(&presecisce, *ray1)) {
fate = hit;
hactive->Fill(presecisce.y(), presecisce.z());
fate = hitExit;
hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
}
if(detector->TestIntersection(&presecisce, *ray1))
1101,13 → 1125,15
if(draw) ray1->DrawS(presecisce.x(), 10.0);
}
} else {
// Main test: ray and SiPM surface
if(active->TestIntersection(&presecisce, *ray1)) {
fate = hit;
hactive->Fill(presecisce.y(), presecisce.z());
hlaser->Fill((in.GetR()).y(), (in.GetR()).z());
fate = hitExit;
hactive->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
hlaser->Fill((in.GetR()).y() + offsetY, (in.GetR()).z() + offsetZ);
}
// If it is on the same plane as SiPM
if(detector->TestIntersection(&presecisce, *ray1))
hdetector->Fill(presecisce.y(), presecisce.z());
hdetector->Fill(offsetY + presecisce.y(), offsetZ + presecisce.z());
}
//} else {
//if(draw) ray1->Draw(presecisce.x(), center.x()+d+window_d);
1114,8 → 1140,9
//}
*out = *ray1;
//delete ray0;
//delete ray1;
delete ray0;
delete ray1;
delete rayout;
return fate;
}
//-----------------------------------------------------------------------------
1188,12 → 1215,12
}
}
 
int Plate::propagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
Fate Plate::propagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
{
CRay ray0;
CRay ray1;
TVector3 vec0, vec1;
int fate = 11;
Fate fate = enter;
int inters_i = 0;
ray0 = in;
1201,14 → 1228,18
int last_hit = 0;
int propagation = 0;
if( !(sides[0]->PropagateRay(ray0, &ray1, &vec1)) ) {
int result = sides[0]->PropagateRay(ray0, &ray1, &vec1);
if( !result ) {
// ce -NI- presecisca z vstopno
fate = 0;
fate = missed;
} else if(result == REFLECTION) {
if (dbg) printf("PLATE: reflected\n");
fate = backreflected;
} else {
points[0] = ray1.GetR();
//hfate->Fill(2);
//hin->Fill(vec1.y(), vec1.z());
while (n_odb++ < VODNIK_MAX_ODB) {
while (n_odb++ < MAX_REFLECTIONS) {
ray0 = ray1;
vec0 = vec1;
propagation = 11;
1219,12 → 1250,12
}
points[n_odb] = vec1;
if(inters_i == 0) {
fate = -3;
fate = backreflected;
break;} // backreflection
propagation = sides[inters_i]->PropagateRay(ray0, &ray1, &vec1);
if(inters_i == 5) { // successfull exit
fate = 1;
fate = hitExit;
//hout->Fill(vec1.y(), vec1.z());
//hnodb_exit->Fill(n_odb-1);
n_odb++;
1233,13 → 1264,17
break;
}
if(propagation == 1) {
fate = -1;
fate = refracted; //at side?
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;} // no total reflection when should be
if(propagation == -2) {fate = -2; break;} // absorption due to finite reflectivity
if(propagation == -2) {
fate = noreflection;
break;
} // absorption due to finite reflectivity
last_hit = inters_i;
}
}