-
Notifications
You must be signed in to change notification settings - Fork 0
/
Scint.cc
96 lines (86 loc) · 2.43 KB
/
Scint.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include "Scint.hh"
extern bool bounce_info;
#include <TRandom3.h>
extern TRandom3 rndm;
#include <cassert>
Line Scint::Reflection(Line photon)
{
if( rndm.Uniform() > reflectivity ) return null;
TVector3 intersect;
int s=-1;
bool yplane;
for(int i=0;i<4;++i)
{
// bounce
intersect = surfaces[i].Intersection( photon );
// avoids numerical errors
if( (intersect - photon.point).Mag() < 1.e-6 ) continue;
// photon and surface are parallel
if( intersect == photon(photon.kParall) ) continue;
// ensure propagation (forward)
if( photon.slope.Z() > 0. && intersect.Z() < photon.point.Z() ) continue;
if( photon.slope.Z() < 0. && intersect.Z() > photon.point.Z() ) continue;
// indetify bouncing surface
yplane = ( (i%2) == 0 );
if( Inside( yplane, intersect) )
{
s=i;
break;
}
}
if( s < 0 )
{
std::cerr<<" Reflection fail"<<std::endl;
return null;
}
//assert(!(s<0));
TVector3 uprime;
if( yplane )
uprime.SetXYZ(photon.slope.X(),-1.*photon.slope.Y(),photon.slope.Z());
else
uprime.SetXYZ(-1.*photon.slope.X(),photon.slope.Y(),photon.slope.Z());
//double thetai = (surfaces[s].normal.Angle(uprime) - TMath::PiOver2())*TMath::RadToDeg();
double thetai = surfaces[s].normal.Angle(uprime) * TMath::RadToDeg();
h_internal_reflection->Fill( thetai );
if( bounce_info ) std::cout<<"Incidence angle "<<thetai<<" deg"<<std::endl;
Line reflray(uprime,intersect,photon.name);
return reflray;
}
std::vector<TVector3> Scint::Propagate(Line ray)
{
std::vector<TVector3> refelections;
int N=0;
while(1)
{
if( bounce_info )
{
std::cout<<"***** Bounce # "<<N<<" for "<<ray.name<<" *****"<<std::endl;
SmartPrint(ray.slope,"direction");
SmartPrint(ray.point,"position");
}
// if( N > 100000 )
// {
// std::cout<<"KILL ME"<<std::endl;
// refelections.clear();
// break;
// }
refelections.push_back( ray.point );
if( det.Hit( ray ) )
{
refelections.push_back( det.hit );
if( bounce_info ) SmartPrint(det.hit,"HIT");
break;
}
else
ray = Reflection( ray );
if( bounce_info ) std::cout<<"******************"<<std::endl;
if( ray.slope.Z() == null.slope.Z() )
{
if( bounce_info ) std::cout<<"Absorbed after "<<N<<" reflections\n"<<std::endl;
refelections.clear();
break;
}
++N;
}
return refelections;
}