Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

LightPipe - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: LightPipe.cc,v 1.12 2005/10/21 13:33:39 moreggia Exp $
// A.Thea created Oct, 26 2002

#include <TMath.h> 
#include "LightPipe.hh"
#include "EVector.hh"
#include "EConst.hh"

using namespace TMath;
using namespace sou;
using namespace EConst;

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()/Clight();
		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()/Clight();
	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() {
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005