// ESAF : Euso Simulation and Analysis Framework
// $Id: LightPipe.cc,v 1.9 2005/04/15 18:54:41 thea Exp $
// A.Thea created Oct, 26 2002
#include <math.h>
#include "LightPipe.hh"
#include "EVector.hh"
#include "SystemOfUnits.hh"
ClassImp(LightPipe)
// constructor
LightPipe::LightPipe(const EVector c1, const EVector c2, double t_x, double t_y, double b_x, double b_y, const PmtGeometry *p):
DetectorPhotonTransporter(), pmt(p) {
par_X=pmt->GetX();
par_Y=pmt->GetY();
par_Z=pmt->GetZ();
pos=pmt->Position();
// rot is the rotation matrix to go from global to local
// inv_rot is the rotation matrix to go from global to local
rot=rot.RotateAxes(par_X, par_Y, par_Z);
inv_rot=rot.Inverse();
EVector loc_par_X(1.,0.,0.);
EVector loc_par_Y(0.,1.,0.);
EVector loc_par_Z(0.,0.,1.);
// set the correct size of the vectors;
corner.resize(2);
normal.resize(6);
edge.resize(6);
// corners is the array of EVectros the EVectors filled
// with the position of the two main corners of the pipe
corner[0]=c1-pos;
corner[1]=c2-pos;
corner[0]=goLocal(corner[0]);
corner[1]=goLocal(corner[1]);
// main edges of the pipes
edge[0]=-t_x*loc_par_X;
edge[1]=-t_y*loc_par_Y;
edge[3]=b_x*loc_par_X;
edge[4]=b_y*loc_par_Y;
edge[2]=corner[1]-corner[0]+edge[3]+edge[4];
edge[5]=corner[0]-corner[1]+edge[0]+edge[1];
// normvec contains the normal vectors to the faces of the pipe
// calculate the normal vectors directed inward
// TOP face
normal[TOP]=edge[1].Cross(edge[0]);
normal[TOP].SetMag(1.);
// BACK
normal[BACK]=edge[0].Cross(edge[2]);
normal[BACK].SetMag(1.);
// RIGHT
normal[RIGHT]=edge[2].Cross(edge[1]);
normal[RIGHT].SetMag(1.);
// BOTTOM
normal[BOTTOM]=edge[3].Cross(edge[4]);
normal[BOTTOM].SetMag(1.);
// FRONT
normal[FRONT]=edge[5].Cross(edge[3]);
normal[FRONT].SetMag(1.);
// LEFT
normal[LEFT]=edge[4].Cross(edge[5]);
normal[LEFT].SetMag(1.);
#ifdef DEBUG
cout << "normal[TOP]: " << normal[TOP] << "\n" <<
"normal[BOTTOM]: " << normal[BOTTOM] << "\n" <<
"normal[LEFT]: " << normal[LEFT] << "\n" <<
"normal[RIGHT]: " << normal[RIGHT] << "\n" <<
"normal[FRONT]: " << normal[FRONT] << "\n" <<
"normal[BACK]: " << normal[BACK] << "\n" << endl;
#endif /* DEBUG */
#ifdef DEBUG
cout << "edge[0]=" <<edge[0] << "\n"
<< "edge[1]=" <<edge[1] << "\n"
<< "edge[2]=" <<edge[2] << "\n"
<< "edge[3]=" <<edge[3] << "\n"
<< "edge[4]=" <<edge[4] << "\n"
<< "edge[5]=" <<edge[5] << "\n"
<< endl;
#endif /* DEBUG */
}
Photon *LightPipe::Transport(Photon *p) const {
vector<EVector> ips=intPoints(p);
// find the hit position
int hitFace=0;
EVector oldpos=p->pos-pos;
oldpos=goLocal(oldpos);
double minlen=(ips[0]-oldpos).Mag();
for(int k=1;k<6;k++)
if ((ips[k]-oldpos).Mag()<minlen) {
minlen=(ips[k]-oldpos).Mag();
hitFace=k;
}
#ifdef DEBUG
cout<<"hit face: ";
switch (hitFace){
case TOP:
cout << "TOP";
break;
case BOTTOM:
cout << "BOTTOM";
break;
case LEFT:
cout << "LEFT";
break;
case RIGHT:
cout << "RIGHT";
break;
case FRONT:
cout << "FRONT";
break;
case BACK:
cout << "BACK";
break;
}
cout<<endl;
cout<<"hit position (local): "<<ips[hitFace]<<endl;
#endif /* DEBUG */
EVector dummy=p->pos;
if (hitFace==TOP || hitFace==BOTTOM) {
p->pos=goGlobal(ips[hitFace]);
p->pos+=pos;
p->time+=(p->pos-dummy).Mag()/c_light;
return p;
}
EVector ldir=goLocal(p->dir);
ldir=ldir-2*(normal[hitFace]*ldir)*normal[hitFace];
p->pos=goGlobal(ips[hitFace])+pos;
p->dir=goGlobal(ldir);
p->time+=(p->pos-dummy).Mag()/c_light;
return Transport(p);
}
int LightPipe::whichFace(Photon *p) const{
EVector lpos=p->pos-pos;
lpos=goLocal(lpos);
if((lpos-corner[0]).Dot(normal[TOP])<=kTolerance){
// photon on the top surface
#ifdef DEBUG
cout<<"TOP"<<endl;
#endif /* DEBUG */
return TOP;
}
else if((lpos-corner[0]).Dot(normal[RIGHT])<=kTolerance){
// photon on the top surface
#ifdef DEBUG
cout<<"RIGHT"<<endl;
#endif /* DEBUG */
return RIGHT;
}
else if((lpos-corner[0]).Dot(normal[BACK])<=kTolerance){
// photon on the top surface
#ifdef DEBUG
cout<<"BACK"<<endl;
#endif /* DEBUG */
return BACK;
}
else if((lpos-corner[1]).Dot(normal[BOTTOM])<=kTolerance){
// photon on the top surface
#ifdef DEBUG
cout<<"BOTTOM"<<endl;
#endif /* DEBUG */
return BOTTOM;
}
else if((lpos-corner[1]).Dot(normal[LEFT])<=kTolerance){
// photon on the top surface
#ifdef DEBUG
cout<<"LEFT"<<endl;
#endif /* DEBUG */
return LEFT;
}
else if((lpos-corner[1]).Dot(normal[FRONT])<=kTolerance){
// photon on the top surface
#ifdef DEBUG
cout<<"FRONT"<<endl;
#endif /* DEBUG */
return FRONT;
}
else throw runtime_error("whichFace: not of surface");
}
vector<EVector> LightPipe::intPoints(Photon *p) const {
EVector lpos, ldir;
lpos=p->pos-pos;
ldir=p->dir;
lpos=goLocal(lpos);
ldir=goLocal(ldir);
int face=whichFace(p);
#ifdef DEBUG
cout<<"lpos: "<<lpos<<endl;
cout<<"ldir: "<<ldir<<endl;
#endif /* DEBUG */
vector<EVector> ips(6);
double dist;
EVector out(1e6*mm,1e6*mm,1e6*mm);
#ifdef DEBUG
cout << "ldir*normal[TOP]: " << ldir*normal[TOP] << "\n" <<
"ldir*normal[BOTTOM]: " << ldir*normal[BOTTOM] << "\n" <<
"ldir*normal[LEFT]: " << ldir*normal[LEFT] << "\n" <<
"ldir*normal[RIGHT]: " << ldir*normal[RIGHT] << "\n" <<
"ldir*normal[FRONT]: " << ldir*normal[FRONT] << "\n" <<
"ldir*normal[BACK]: " << ldir*normal[BACK] << "\n" << endl;
#endif /* DEBUG */
if (face==TOP || ldir.Dot(normal[TOP]) > 0) ips[TOP]=out;
else{
dist=(lpos-corner[0]).Dot(normal[TOP]);
ips[TOP]=lpos-(dist/ldir.Dot(normal[TOP]))*ldir;
#ifdef DEBUG
cout<<"dist, ips[TOP]: "<<dist<<", "<<ips[TOP]<<endl;
#endif /* DEBUG */
}
if (face==BOTTOM || ldir.Dot(normal[BOTTOM]) > 0) ips[BOTTOM]=out;
else{
dist=(lpos-corner[1]).Dot(normal[BOTTOM]);
ips[BOTTOM]=lpos-ldir*(dist/ldir.Dot(normal[BOTTOM]));
#ifdef DEBUG
cout<<"dist, ips[BOTTOM]: "<<dist<<", "<<ips[BOTTOM]<<endl;
#endif /* DEBUG */
}
if (face==LEFT || ldir.Dot(normal[LEFT]) > 0) ips[LEFT]=out;
else{
dist=(lpos-corner[1]).Dot(normal[LEFT]);
ips[LEFT]=lpos-ldir*(dist/ldir.Dot(normal[LEFT]));
#ifdef DEBUG
cout<<"dist, ips[LEFT]: "<<dist<<", "<<ips[LEFT]<<endl;
#endif /* DEBUG */
}
if (face==RIGHT || ldir.Dot(normal[RIGHT]) > 0) ips[RIGHT]=out;
else{
dist=(lpos-corner[0]).Dot(normal[RIGHT]);
ips[RIGHT]=lpos-ldir*(dist/ldir.Dot(normal[RIGHT]));
#ifdef DEBUG
cout<<"dist, ips[RIGHT]: "<<dist<<", "<<ips[RIGHT]<<endl;
#endif /* DEBUG */
}
if (face==FRONT || ldir.Dot(normal[FRONT]) > 0) ips[FRONT]=out;
else{
dist=(lpos-corner[1]).Dot(normal[FRONT]);
ips[FRONT]=lpos-ldir*(dist/ldir.Dot(normal[FRONT]));
#ifdef DEBUG
cout<<"dist, ips[FRONT]: "<<dist<<", "<<ips[FRONT]<<endl;
#endif /* DEBUG */
}
if (face==BACK || ldir.Dot(normal[BACK]) > 0) ips[BACK]=out;
else{
dist=(lpos-corner[0]).Dot(normal[BACK]);
ips[BACK]=lpos-ldir*(dist/ldir.Dot(normal[BACK]));
#ifdef DEBUG
cout<<"dist, ips[BACK]: "<<dist<<", "<<ips[BACK]<<endl;
#endif /* DEBUG */
}
return ips;
}
// destructor
LightPipe::~LightPipe() {
}