source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/optics/src/LightPipe.cc @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 7.7 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: LightPipe.cc 2250 2005-10-21 13:33:39Z moreggia $
3// A.Thea created Oct, 26 2002
4
5#include <TMath.h> 
6#include "LightPipe.hh"
7#include "EVector.hh"
8#include "EConst.hh"
9
10using namespace TMath;
11using namespace sou;
12using namespace EConst;
13
14ClassImp(LightPipe)
15
16// constructor
17LightPipe::LightPipe(const EVector c1, const EVector c2, double t_x, double t_y, double b_x, double b_y, const PmtGeometry *p):
18DetectorPhotonTransporter(), pmt(p) {
19
20    par_X=pmt->GetX();
21    par_Y=pmt->GetY();
22    par_Z=pmt->GetZ();
23    pos=pmt->Position();
24
25    // rot is the rotation matrix to go from global to local
26    // inv_rot is the rotation matrix to go from global to local
27    rot=rot.RotateAxes(par_X, par_Y, par_Z);
28    inv_rot=rot.Inverse();
29
30    EVector loc_par_X(1.,0.,0.);
31    EVector loc_par_Y(0.,1.,0.);
32    EVector loc_par_Z(0.,0.,1.);
33
34    // set the correct size of the vectors;
35    corner.resize(2);
36    normal.resize(6);
37    edge.resize(6);
38   
39    // corners is the array of EVectros the EVectors filled
40    // with the position of the two main corners of the pipe
41    corner[0]=c1-pos;
42    corner[1]=c2-pos;
43    corner[0]=goLocal(corner[0]);
44    corner[1]=goLocal(corner[1]);
45   
46    // main edges of the pipes
47   
48    edge[0]=-t_x*loc_par_X;
49    edge[1]=-t_y*loc_par_Y;
50    edge[3]=b_x*loc_par_X;
51    edge[4]=b_y*loc_par_Y;
52    edge[2]=corner[1]-corner[0]+edge[3]+edge[4];
53    edge[5]=corner[0]-corner[1]+edge[0]+edge[1];
54   
55    // normvec contains the normal vectors to the faces of the pipe
56    // calculate the normal vectors directed inward
57
58    // TOP face 
59    normal[TOP]=edge[1].Cross(edge[0]); 
60    normal[TOP].SetMag(1.); 
61    // BACK
62    normal[BACK]=edge[0].Cross(edge[2]);
63    normal[BACK].SetMag(1.);
64    // RIGHT
65    normal[RIGHT]=edge[2].Cross(edge[1]);
66    normal[RIGHT].SetMag(1.);
67    // BOTTOM
68    normal[BOTTOM]=edge[3].Cross(edge[4]);
69    normal[BOTTOM].SetMag(1.);
70    // FRONT
71    normal[FRONT]=edge[5].Cross(edge[3]);
72    normal[FRONT].SetMag(1.);
73    // LEFT
74    normal[LEFT]=edge[4].Cross(edge[5]);
75    normal[LEFT].SetMag(1.);
76   
77#ifdef DEBUG
78    cout << "normal[TOP]: " << normal[TOP] << "\n" <<
79            "normal[BOTTOM]: " << normal[BOTTOM] << "\n" <<
80            "normal[LEFT]: " << normal[LEFT] << "\n" <<
81            "normal[RIGHT]: " << normal[RIGHT] << "\n" <<
82            "normal[FRONT]: " << normal[FRONT] << "\n" <<
83            "normal[BACK]: " << normal[BACK] << "\n" << endl;
84#endif /* DEBUG */
85#ifdef DEBUG
86    cout << "edge[0]=" <<edge[0] << "\n"
87         << "edge[1]=" <<edge[1] << "\n"
88         << "edge[2]=" <<edge[2] << "\n"
89         << "edge[3]=" <<edge[3] << "\n"
90         << "edge[4]=" <<edge[4] << "\n"
91         << "edge[5]=" <<edge[5] << "\n"
92         << endl;
93#endif /* DEBUG */
94   
95}
96
97Photon *LightPipe::Transport(Photon *p) const {
98        vector<EVector> ips=intPoints(p);
99
100        // find the hit position
101        int hitFace=0;
102        EVector oldpos=p->pos-pos;
103        oldpos=goLocal(oldpos);
104        double minlen=(ips[0]-oldpos).Mag();
105        for(int k=1;k<6;k++)
106           if ((ips[k]-oldpos).Mag()<minlen) {
107                   minlen=(ips[k]-oldpos).Mag();
108                   hitFace=k;
109           }
110       
111#ifdef DEBUG
112        cout<<"hit face: ";
113        switch (hitFace){
114                case TOP:
115                     cout << "TOP";
116                break;
117                case BOTTOM:
118                     cout << "BOTTOM";
119                break;
120                case LEFT:
121                     cout << "LEFT";
122                break;
123                case RIGHT:
124                     cout << "RIGHT";
125                break;
126                case FRONT:
127                     cout << "FRONT";
128                break;
129                case BACK:
130                     cout << "BACK";
131                break;
132                       
133        }
134        cout<<endl;
135        cout<<"hit position (local): "<<ips[hitFace]<<endl;
136#endif /* DEBUG */
137
138        EVector dummy=p->pos;
139        if (hitFace==TOP || hitFace==BOTTOM) {
140                p->pos=goGlobal(ips[hitFace]);
141                p->pos+=pos;
142                p->time+=(p->pos-dummy).Mag()/Clight();
143                return p;
144        }
145       
146        EVector ldir=goLocal(p->dir);
147        ldir=ldir-2*(normal[hitFace]*ldir)*normal[hitFace];
148        p->pos=goGlobal(ips[hitFace])+pos;
149        p->dir=goGlobal(ldir);
150        p->time+=(p->pos-dummy).Mag()/Clight();
151        return Transport(p);
152 
153}
154
155
156int LightPipe::whichFace(Photon *p) const{
157
158       
159
160    EVector lpos=p->pos-pos;
161    lpos=goLocal(lpos);
162
163    if((lpos-corner[0]).Dot(normal[TOP])<=kTolerance){
164        // photon on the top surface
165  #ifdef DEBUG
166        cout<<"TOP"<<endl;
167  #endif /* DEBUG */
168    return TOP;
169    }
170
171    else if((lpos-corner[0]).Dot(normal[RIGHT])<=kTolerance){
172        // photon on the top surface
173  #ifdef DEBUG
174        cout<<"RIGHT"<<endl;
175  #endif /* DEBUG */
176    return RIGHT;
177    }
178   
179    else if((lpos-corner[0]).Dot(normal[BACK])<=kTolerance){
180        // photon on the top surface
181  #ifdef DEBUG
182        cout<<"BACK"<<endl;
183  #endif /* DEBUG */
184    return BACK;
185    }
186   
187    else if((lpos-corner[1]).Dot(normal[BOTTOM])<=kTolerance){
188        // photon on the top surface
189  #ifdef DEBUG
190        cout<<"BOTTOM"<<endl;
191  #endif /* DEBUG */
192    return BOTTOM;
193    }
194   
195    else if((lpos-corner[1]).Dot(normal[LEFT])<=kTolerance){
196        // photon on the top surface
197  #ifdef DEBUG
198        cout<<"LEFT"<<endl;
199  #endif /* DEBUG */
200    return LEFT;
201    }
202   
203    else if((lpos-corner[1]).Dot(normal[FRONT])<=kTolerance){
204        // photon on the top surface
205  #ifdef DEBUG
206        cout<<"FRONT"<<endl;
207  #endif /* DEBUG */
208    return FRONT;
209    }
210    else throw runtime_error("whichFace: not of surface");
211
212}
213
214vector<EVector> LightPipe::intPoints(Photon *p) const {
215
216    EVector lpos, ldir;
217    lpos=p->pos-pos;
218    ldir=p->dir;
219    lpos=goLocal(lpos);
220    ldir=goLocal(ldir);
221   
222    int face=whichFace(p);
223   
224#ifdef DEBUG
225        cout<<"lpos: "<<lpos<<endl;
226        cout<<"ldir: "<<ldir<<endl;
227#endif /* DEBUG */
228   
229    vector<EVector> ips(6);
230    double dist;
231    EVector out(1e6*mm,1e6*mm,1e6*mm);
232   
233
234#ifdef DEBUG
235    cout << "ldir*normal[TOP]: " << ldir*normal[TOP] << "\n" <<
236            "ldir*normal[BOTTOM]: " << ldir*normal[BOTTOM] << "\n" <<
237            "ldir*normal[LEFT]: " << ldir*normal[LEFT] << "\n" <<
238            "ldir*normal[RIGHT]: " << ldir*normal[RIGHT] << "\n" <<
239            "ldir*normal[FRONT]: " << ldir*normal[FRONT] << "\n" <<
240            "ldir*normal[BACK]: " << ldir*normal[BACK] << "\n" << endl;
241#endif /* DEBUG */
242
243    if (face==TOP || ldir.Dot(normal[TOP]) > 0) ips[TOP]=out;
244    else{
245            dist=(lpos-corner[0]).Dot(normal[TOP]);
246            ips[TOP]=lpos-(dist/ldir.Dot(normal[TOP]))*ldir;
247#ifdef DEBUG
248                cout<<"dist, ips[TOP]: "<<dist<<", "<<ips[TOP]<<endl;
249#endif /* DEBUG */
250    } 
251
252    if (face==BOTTOM || ldir.Dot(normal[BOTTOM]) > 0) ips[BOTTOM]=out;
253    else{
254            dist=(lpos-corner[1]).Dot(normal[BOTTOM]);
255            ips[BOTTOM]=lpos-ldir*(dist/ldir.Dot(normal[BOTTOM]));
256#ifdef DEBUG
257                cout<<"dist, ips[BOTTOM]: "<<dist<<", "<<ips[BOTTOM]<<endl;
258#endif /* DEBUG */
259    } 
260
261    if (face==LEFT || ldir.Dot(normal[LEFT]) > 0) ips[LEFT]=out;
262    else{
263            dist=(lpos-corner[1]).Dot(normal[LEFT]);
264            ips[LEFT]=lpos-ldir*(dist/ldir.Dot(normal[LEFT]));
265#ifdef DEBUG
266                cout<<"dist, ips[LEFT]: "<<dist<<", "<<ips[LEFT]<<endl;
267#endif /* DEBUG */
268    } 
269   
270    if (face==RIGHT || ldir.Dot(normal[RIGHT]) > 0) ips[RIGHT]=out;
271    else{
272            dist=(lpos-corner[0]).Dot(normal[RIGHT]);
273            ips[RIGHT]=lpos-ldir*(dist/ldir.Dot(normal[RIGHT]));
274#ifdef DEBUG
275                cout<<"dist, ips[RIGHT]: "<<dist<<", "<<ips[RIGHT]<<endl;
276#endif /* DEBUG */
277    } 
278   
279    if (face==FRONT || ldir.Dot(normal[FRONT]) > 0) ips[FRONT]=out;
280    else{
281            dist=(lpos-corner[1]).Dot(normal[FRONT]);
282            ips[FRONT]=lpos-ldir*(dist/ldir.Dot(normal[FRONT]));
283#ifdef DEBUG
284                cout<<"dist, ips[FRONT]: "<<dist<<", "<<ips[FRONT]<<endl;
285#endif /* DEBUG */
286    }
287   
288    if (face==BACK || ldir.Dot(normal[BACK]) > 0) ips[BACK]=out;
289    else{
290            dist=(lpos-corner[0]).Dot(normal[BACK]);
291            ips[BACK]=lpos-ldir*(dist/ldir.Dot(normal[BACK]));
292#ifdef DEBUG
293                cout<<"dist, ips[BACK]: "<<dist<<", "<<ips[BACK]<<endl;
294#endif /* DEBUG */
295    }
296   
297   
298    return ips;
299}
300
301// destructor
302LightPipe::~LightPipe() {
303}
Note: See TracBrowser for help on using the repository browser.