Subversion Repositories f9daq

Compare Revisions

Regard whitespace Rev 71 → Rev 72

/lightguide/trunk/include/guide.h
40,20 → 40,23
CRay() :
r(TVector3(0,0,0)),
n(TVector3(1,0,0)),
p(TVector3(0,1,0)),
color(1)
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
CRay(TVector3 r0, TVector3 n0) :
r(r0),
n(n0.Unit()),
p(TVector3(0,1,0)),
color(1)
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
CRay(double x0, double y0, double z0, double l0, double m0, double n0) :
r(TVector3(x0,y0,z0)),
n(TVector3(l0,m0,n0).Unit()),
p(TVector3(0,1,0)),
color(1)
//p(TVector3(0,1,0)),
p(TVector3(0,0,1)),
color(kBlack)
{};
 
void Set(TVector3 r0, TVector3 n0);
88,7 → 91,7
// ^ |
// | |
// | |
// | �
// |
// r0<-----r3
class CPlane4
{
119,9 → 122,28
double angle_r[4]; // koti ob posameznem vogalu
};
//=================================================================================
// ravnina - krog
class CPlaneR
{
 
public:
CPlaneR(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
 
void Set(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
 
int TestIntersection(TVector3 *vec, CRay in);
 
void Draw(int color = 1, int NN = 32);
 
private:
TVector3 n, center;
double _r;
};
//=================================================================================
 
//=================================================================================
// ravna opticna povrsina: refractor, zrcalo ali povrsina s totalnim odbojem
#define SURF_DUMMY 0
#define SURF_REFRA 1
190,7 → 212,7
class DetectorParameters
{
public:
DetectorParameters(double a, double b, double d, double active, double n1, double n2, double n3, TVector3 gap):
DetectorParameters(double a, double b, double d, double active, double n1, double n2, double n3, TVector3 gap, bool coupling):
_a(a),
_b(b),
_d(d),
206,7 → 228,8
_plateOn(1),
_plateWidth(1),
_glassOn(0),
_glassD(0)
_glassD(0),
_coupling(coupling)
{};
DetectorParameters(double a, double b, double d):
_a(a),
224,7 → 247,8
_plateOn(1),
_plateWidth(1),
_glassOn(0),
_glassD(0)
_glassD(0),
_coupling(false)
{};
~DetectorParameters() {};
240,6 → 264,7
void setGuideOn(int guideOn) { _guideOn = guideOn; };
void setPlate(int plateOn, double plateWidth) { _plateOn = plateOn; _plateWidth = plateWidth; };
void setIndices(double n1, double n2, double n3) { _n1 = n1; _n2 = n2; _n3 = n3; };
void setCoupling(bool coupling) { _coupling = coupling; };
 
double getLightYield() { return ( Power(_b,2)/Power(_active,2)); };
double getM() { return (_b/_a); };
259,6 → 284,8
int getPlateOn() { return _plateOn; };
double getOffsetY() { return offsetY; };
double getOffsetZ() { return offsetZ; };
// bad coupling in the case the small amount of grease was applied
bool badCoupling() { return _coupling; };
private:
double _a;
281,6 → 308,7
int _glassOn;
double _glassD;
bool _coupling;
};
 
290,6 → 318,8
Guide(TVector3 center0, DetectorParameters& parameters);
~Guide() {
for (int jk=0; jk<6; jk++) delete s_side[jk];
delete grease;
delete noCoupling;
delete hfate;
delete hnodb_all;
delete hnodb_exit;
300,9 → 330,9
Fate PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points);
double getD() { return _d; }
double getN1() { return n1; };
double getN2() { return n2; };
double getN3() { return n3; };
double getN1() { return _n1; };
double getN2() { return _n2; };
double getN3() { return _n3; };
TH1F* GetHFate() const {return hfate;};
TH1F* GetHNOdb_all() const {return hnodb_all;};
TH1F* GetHNOdb_exit() const {return hnodb_exit;};
319,12 → 349,15
private:
Fate fate;
double _d; // parameters
double n1, n2, n3; // refractive index: n1 above entry surface, n2 inside, n3 after exit
double _n1, _n2, _n3; // refractive index: n1 above entry surface, n2 inside, n3 after exit
double _r;
double absorption;
double A;
double _absorption;
double _A;
bool _badCoupling;
TRandom rand; // for material absorption
CSurface *s_side[6];
CPlaneR* grease;
CSurface* noCoupling;
TVector3 center, vodnik_edge[8];
TH1F *hfate, *hnodb_all, *hnodb_exit;
333,27 → 366,8
//=================================================================================
 
//=============================================================================================================================== <<<<<<<<
// ravnina - krog
class CPlaneR
{
 
public:
CPlaneR(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
void Set(TVector3 c, TVector3 n0, double R0)
{center = c; n = n0; _r = R0;};
int TestIntersection(TVector3 *vec, CRay in);
void Draw(int color = 1, int NN = 32);
private:
TVector3 n, center;
double _r;
};
//=================================================================================
 
class Plate
{
public:
464,6 → 478,7
 
//double x_gap, y_gap, z_gap;
CPlane4 *active;
CPlaneR* grease;
TH2F *hactive, *hlaser;
 
CPlane4 *detector;
/lightguide/trunk/src/userFunctions.cpp
12,8 → 12,8
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.46, TVector3(0.3, 0, 0));
// Set the initial detector parameters (a, b, d, activeDetector, n1, n2, n3, gap, badCoupling)
DetectorParameters parameters(3.0, 5.0, 3.0, 3.0, 1, 1.53, 1.46, TVector3(0.3, 0, 0), false);
// Print the detector parameters
void getParameters();
 
193,9 → 193,10
}
//------------------------------------------------------------------------------------------
// Propagate NN rays generated as grid and show the statistics
void LGG(int NN = 10, double theta = 0.0)
void LGG(int NN=10, double theta=0.0, bool coupling=false)
{
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
parameters.setCoupling(coupling);
CDetector *detector = new CDetector(CENTER, parameters);
Init();
double izkoristek = Grid(detector, parameters, NN, theta);
//printf("izkoristek = %.3lf\n", izkoristek);
207,7 → 208,7
// Propagate NN rays genarated under the same angle theta, phi with random spacing and statistics
void LGR(int NN = 1e4, double theta = 0.0, double phi = 0.0)
{
CDetector *detector = new CDetector(TVector3(-2,0,0), parameters);
CDetector *detector = new CDetector(CENTER, parameters);
Init();
double izkoristek = RandYZ(detector, parameters, NN, theta, phi, 30);
//printf("izkoristek = %.3lf\n", izkoristek);
231,16 → 232,19
//detector->Draw();
}
 
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, int nnrays = 30, int showr = 0)
void LGB(int NN = 1e4, double phiMin=-19.4, double phiMax=19.4, bool coupling=true, int nnrays = 30, int showr = 0)
{
Init();
parameters.setCoupling(coupling);
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);
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();
/lightguide/trunk/src/guide.cpp
475,6 → 475,7
}
out->Set(intersect, transmit);
// Shift?
pol_t = -in.GetP() + sign_n*2*cosTi*n;
out->SetPolarization(pol_t);
return REFLECTION;
582,24 → 583,25
 
 
//=================================================================================
Guide::Guide(TVector3 center0, DetectorParameters &parameters)
Guide::Guide(TVector3 center0, DetectorParameters &parameters) :
_d(parameters.getD()),
_n1(parameters.getN1()),
_n2(parameters.getN2()),
_n3(parameters.getN3()),
_r(c_reflectivity),
_absorption(0),
_A(0),
_badCoupling(parameters.badCoupling())
{
double t;
 
TDatime now;
rand.SetSeed(now.Get());
 
center = center0;
double b = parameters.getB();
double a = parameters.getA();
_d = parameters.getD();
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() ? n2 : n1);
n3 = parameters.getN3();
_r = c_reflectivity;
double n0 = (parameters.getPlateOn() ? _n2 : _n1);
int fresnel = parameters.getFresnel();
 
// light guide edges
617,16 → 619,28
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] = 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); // n3 - ref ind at the exit, grease, air, epoxy
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);
// n3 - ref ind at the exit, grease, air
s_side[5] = new CSurface(SURF_REFRA, &vodnik_edge[4], _n2, _n3, _r);
s_side[5]->FlipN();
// exit surface in the case of bad coupling
noCoupling = new CSurface(SURF_REFRA, &vodnik_edge[4], _n2, 1.0, _r);
noCoupling->FlipN();
// grease = specific pattern area of coupling
TVector3 activePosition(center);
activePosition += TVector3(_d, 0, 0);
TVector3 normal(1,0,0);
grease = new CPlaneR(activePosition, normal, a/2.0);
if(fresnel) for(int i=0; i<6; i++) s_side[i]->SetFresnel(1);
643,20 → 657,17
(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", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
hnodb_all = new TH1F("hnodb_all", "", 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", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
hnodb_exit = new TH1F("hnodb_exit", "", MAX_REFLECTIONS, -0.5, MAX_REFLECTIONS-0.5);
int nBins = nch + 1;
hin = (TH2F*)gROOT->FindObject("hin"); if(hin) delete hin;
hin = new TH2F("hin", "Guide entrance window", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0);
hin = new TH2F("hin", ";x [mm]; y[mm]", nBins, -b/2.0, +b/2.0, nBins, -b/2.0, +b/2.0);
hout = (TH2F*)gROOT->FindObject("hout"); if(hout) delete hout;
hout = new TH2F("hout", "Guide exit window", nBins, -a/2.0, +a/2.0, nBins, -a/2.0, +a/2.0);
absorption = 0;
A = 0;
hout = new TH2F("hout", ";x [mm];y [mm]", nBins, -a/2.0, +a/2.0, nBins, -a/2.0, +a/2.0);
}
//-----------------------------------------------------------------------------
// Sledi zarku skozi vodnik. Vrne:
670,6 → 681,8
Fate Guide::PropagateRay(CRay in, CRay *out, int *n_points, TVector3 *points)
{
if (dbg) printf("--- GUIDE::PropagateRay ---\n");
// ray0 - incident ray
// ray1 - trans/refl ray
CRay ray0;
CRay ray1;
TVector3 vec0, vec1;
723,24 → 736,22
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
if(inters_i == 5) {
if (_badCoupling) {
TVector3 hitVector(0,0,0);
bool hitActive = grease->TestIntersection(&hitVector, ray0);
if (hitActive and dbg) printf(" GUIDE: hit grease\n");
if (!hitActive) propagation = noCoupling->PropagateRay(ray0, &ray1, &vec1);
}
// check on which side the vector is?
TVector3 ray = ray1.GetN();
TVector3 exitNormal = s_side[5]->GetN();
//printf("theta(ray) = %lf, theta(normal5) = %lf ", ray.Theta()*DEGREE, exitNormal.Theta()*DEGREE);
//printf("phi(ray) = %lf, phi(normal5) = %lf\n", ray.Phi()*DEGREE, exitNormal.Phi()*DEGREE);
if (dbg) printf("ray*n_5 = %lf\n", ray*exitNormal);
if (ray*exitNormal > 0) {
if (dbg) printf(" GUIDE: ray is backreflected from exit window.\n");
758,22 → 769,31
ray0 = ray1;
break;
}
 
if(propagation == REFRACTION) {
fate = refracted;
n_odb++;
points[n_odb] = vec1;
ray0 = ray1;
break;
} // no total reflection when should be
 
last_hit = inters_i;
}
}
//--- material absorption ---
if(absorption) {
if(_absorption) {
double travel = 0.0;
printf("n_odb = %d\n", n_odb); //dbg
if (dbg) printf("n_odb = %d\n", n_odb);
for(int point = 0; point < n_odb-1; point++) {
travel += (points[point] - points[point+1]).Mag();
printf("travel = %lf\n", travel); //dbg
if (dbg) printf("travel = %lf\n", travel);
}
double T_abs = TMath::Exp(-travel/A);
printf("T_abs = %lf\n", T_abs); //dbg
double T_abs = TMath::Exp(-travel/_A);
if(dbg)printf("T_abs = %lf\n", T_abs);
double p_abs = rand.Uniform(0.0, 1.0);
printf("p_abs = %lf\n", p_abs); //dbg
if(dbg)printf("p_abs = %lf\n", p_abs);
if(p_abs > T_abs) fate = absorbed; // absorption
}
876,12 → 896,6
center(center0),
glass_on(parameters.getGlassOn()),
glass_d(parameters.getGlassD()),
//x_gap(parameters.getGap().X()),
//y_gap(parameters.getGap().Y()),
//z_gap(parameters.getGap().Z()),
//glass(new CSurface),
//glass_circle(new CPlaneR),
//active(new CPlane4),
col_in(2),
col_lg(8),
col_out(4),
890,8 → 904,6
col_glass(4),
col_active(7),
guide_on(parameters.getGuideOn()),
//window_R( parameters.getB() ),
//window_d(0),
guide(new Guide(center0, parameters)),
plate(new Plate(parameters)),
_plateWidth(parameters.getPlateWidth()),
909,8 → 921,6
if(guide_on) x_offset = center.x();
else x_offset = center.x() - d;
//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();
937,23 → 947,28
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;
hglass = new TH2F("hglass", "Hits glass",
hglass = new TH2F("hglass", "",
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;
if (dbg) cout<<"SiPM active length "<<parameters.getActive()<<endl;
 
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);
plane_v[3].SetXYZ(x_offset+d+x_gap, y_gap - p_size, z_gap - p_size);
active = new CPlane4(plane_v);
//active surface in case of bad coupling is circle d=a
TVector3 activePosition(center);
activePosition += TVector3(d + x_gap, 0, 0);
TVector3 normal(1,0,0);
grease = new CPlaneR(activePosition, normal, 1.0*p_size);
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 + offsetY, y_gap + p_size + offsetY, nBins, z_gap - p_size + offsetZ, z_gap + p_size + offsetZ);
hactive = new TH2F("hactive", ";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);
p_size = b/2.0;
//p_size = 2.5;
995,7 → 1010,6
// 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