1 | // $Id: SphericalIdealFocalSurface.cc 2789 2008-02-11 07:26:05Z naumov $ |
---|
2 | // Author: A.Thea 2005/01/16 |
---|
3 | |
---|
4 | /***************************************************************************** |
---|
5 | * ESAF: Euso Simulation and Analysis Framework * |
---|
6 | * * |
---|
7 | * Id: SphericalIdealFocalSurface * |
---|
8 | * Package: Optics * |
---|
9 | * Coordinator: Alessandro.Thea * |
---|
10 | * * |
---|
11 | *****************************************************************************/ |
---|
12 | |
---|
13 | //_____________________________________________________________________________ |
---|
14 | // |
---|
15 | // SphericalIdealFocalSurface |
---|
16 | // |
---|
17 | // Describes a spherical focal surface. |
---|
18 | // |
---|
19 | // Config file parameters |
---|
20 | // ====================== |
---|
21 | // |
---|
22 | // fR [mm] : Maximum extension of the ideal focal surface from the optical |
---|
23 | // axis. |
---|
24 | // |
---|
25 | // fPos [mm] : z coodinate of the tip of the surface in detector coordinate |
---|
26 | // system. |
---|
27 | // |
---|
28 | // fSphRadius [mm] : radius of the sphere the focal surface is part of. |
---|
29 | // |
---|
30 | // fPrec [mm] : precision required to claim that the photon has hit the FS. |
---|
31 | // |
---|
32 | |
---|
33 | #include "SphericalIdealFocalSurface.hh" |
---|
34 | #include <TMath.h> |
---|
35 | |
---|
36 | using namespace sou; |
---|
37 | |
---|
38 | ClassImp(SphericalIdealFocalSurface) |
---|
39 | |
---|
40 | const Int_t kMaxIter=50; |
---|
41 | |
---|
42 | //_____________________________________________________________________________ |
---|
43 | SphericalIdealFocalSurface::SphericalIdealFocalSurface() { |
---|
44 | // |
---|
45 | // Constructor |
---|
46 | // |
---|
47 | |
---|
48 | fR = Conf()->GetNum("SphericalIdealFocalSurface.fR")*mm; |
---|
49 | fPos = EVector(0,0,Conf()->GetNum("SphericalIdealFocalSurface.fPos.Z"))*mm; |
---|
50 | fSphRadius = Conf()->GetNum("SphericalIdealFocalSurface.fSphRadius")*mm; |
---|
51 | fPrec = Conf()->GetNum("KIdealFocalSurface.fPrec")*mm; |
---|
52 | fDZdown = TMath::Abs(Zfs(fR)); //TOCHECK |
---|
53 | |
---|
54 | } |
---|
55 | |
---|
56 | //_____________________________________________________________________________ |
---|
57 | SphericalIdealFocalSurface::~SphericalIdealFocalSurface() { |
---|
58 | // |
---|
59 | // Destructor |
---|
60 | // |
---|
61 | } |
---|
62 | |
---|
63 | //______________________________________________________________________________ |
---|
64 | Bool_t SphericalIdealFocalSurface::HitPosition(Photon *p) |
---|
65 | { |
---|
66 | // |
---|
67 | // Find whether the photons hits the ideal surface, and in such a case |
---|
68 | // calculated the final position |
---|
69 | // |
---|
70 | |
---|
71 | |
---|
72 | // reject photons going downward |
---|
73 | if ( p->dir.Z() < 0 ) return kFALSE; |
---|
74 | |
---|
75 | TVector3 in_pos = p->pos-fPos; |
---|
76 | EVector ax_pos(0,0,-fDZdown), ax_dir(0,0,1); |
---|
77 | |
---|
78 | EVector bottom_hit = in_pos+p->dir*((-fDZdown-in_pos[Z])/p->dir[Z]); |
---|
79 | |
---|
80 | if (bottom_hit.Perp() > fR) return kFALSE; |
---|
81 | |
---|
82 | //the minimum distance between the photon and the Zfs axis |
---|
83 | EVector dummy=ax_dir.Cross(p->dir); |
---|
84 | double min_dist=TMath::Abs((bottom_hit-ax_pos).Dot(dummy))/dummy.Mag(); |
---|
85 | |
---|
86 | //interaction point on top of the cylinder |
---|
87 | //limit inclination of the photon |
---|
88 | |
---|
89 | EVector top_hit; |
---|
90 | if (p->dir[Z]/p->dir.Perp() < (fDZup+fDZdown)/(2*fR)) { |
---|
91 | // photon hit the lateral surface of the cylinder |
---|
92 | top_hit=bottom_hit+p->dir* |
---|
93 | (fR*sqrt(1-min_dist*min_dist/(fR*fR))/p->dir.Perp()); |
---|
94 | } else { |
---|
95 | // check if the photon hit the top of the cylinder |
---|
96 | top_hit=bottom_hit+p->dir*((fDZup+fDZdown)/p->dir[Z]); |
---|
97 | if (top_hit.Perp() > fR) { |
---|
98 | top_hit=bottom_hit+p->dir* |
---|
99 | (fR*sqrt(1-min_dist*min_dist/(fR*fR))/p->dir.Perp()); |
---|
100 | } |
---|
101 | } |
---|
102 | #ifdef DEBUG |
---|
103 | cout << "p->dir " << p->dir << endl; |
---|
104 | cout << "top_hit= " << top_hit << " bottom_hit=" << bottom_hit << endl; |
---|
105 | #endif |
---|
106 | EVector step=(top_hit-bottom_hit)*.5; |
---|
107 | dummy=bottom_hit+step; |
---|
108 | for(int i=0; i < kMaxIter ; i++) { |
---|
109 | #ifdef DEBUG |
---|
110 | cout << "i=" << i << " dummy[Z]=" << dummy[Z] << " Zfs(dummy.Perp())=" << Zfs(dummy.Perp()) << endl; |
---|
111 | #endif |
---|
112 | if (TMath::Abs(Zfs(dummy.Perp())-dummy[Z]) < fPrec) { |
---|
113 | break; |
---|
114 | } |
---|
115 | step*=.5; |
---|
116 | if (Zfs(dummy.Perp()) > dummy[Z]) dummy+=step; |
---|
117 | else dummy-=step; |
---|
118 | } |
---|
119 | |
---|
120 | if ( dummy.Perp() < fR ) { |
---|
121 | p->posOnIfs=dummy+fPos; |
---|
122 | p->hitIfs=kTRUE; |
---|
123 | return kTRUE; |
---|
124 | } |
---|
125 | |
---|
126 | return kFALSE; |
---|
127 | |
---|
128 | } |
---|
129 | |
---|
130 | |
---|
131 | //______________________________________________________________________________ |
---|
132 | Double_t SphericalIdealFocalSurface::Zfs(Double_t r) { |
---|
133 | // |
---|
134 | // Zfs returns the z value of the focal surface as a function of radius in |
---|
135 | // local coordinates |
---|
136 | // |
---|
137 | Double_t x = r/fSphRadius; |
---|
138 | return -fSphRadius*(1-TMath::Sqrt(1.-x*x)); |
---|
139 | } |
---|
140 | |
---|
141 | //______________________________________________________________________________ |
---|
142 | Double_t SphericalIdealFocalSurface::Profile(Double_t r) { |
---|
143 | // |
---|
144 | // Profile returns the z value of the focal surface as a function of radius |
---|
145 | // in detector coordinates |
---|
146 | // |
---|
147 | |
---|
148 | return Zfs(r)+fPos.Z(); |
---|
149 | } |
---|