// $Id: OpticalAdaptor.cc,v 1.24 2005/04/27 14:22:12 thea Exp $
// Author: D.Demarco, M.Pallavicini
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: OpticalAdaptor *
* Package: Optics *
* Coordinator: Alessandro.Thea *
* *
*****************************************************************************/
//______________________________________________________________________________
//
// Optical Adapter abstract interface
// ==================================
//
// Abstarct base class providing the interface of an Optical Adaptor.
// It provides some general methods that can be used in the actual
// implementation.
//
#include <math.h>
#include "TVector.h"
#include "OpticalAdaptor.hh"
#include "euso.hh"
using namespace TMath;
ClassImp(OpticalAdaptor)
Interpolate* OpticalAdaptor::fgBG3 = 0;
Interpolate* OpticalAdaptor::fgMultiLayer = 0;
//______________________________________________________________________________
OpticalAdaptor::OpticalAdaptor() : pmt(0), fSide(0), fHeight(0), fReflectivity(0), fRadiomEff(1.) {
//
// Constructor
//
fFilterType = kSquare; //default
if (!fgBG3) {
fgBG3 = new Interpolate("config/Optics/OpticalAdaptor/BG3.dat",2);
fgBG3->SetXUnit(nm);
}
if (!fgMultiLayer) {
fgMultiLayer = new Interpolate("config/Optics/OpticalAdaptor/Multilayer.dat",2);
fgMultiLayer->SetXUnit(nm);
}
}
//______________________________________________________________________________
void OpticalAdaptor::SetGeometry( const PmtGeometry *g ) {
//
// Associate this with PmtGeometry g
//
pmt = g;
fXaxis=pmt->GetX();
fYaxis=pmt->GetY();
fZaxis=pmt->GetZ();
fPos=pmt->Position();
// fRot is the rotation matrix to go from global to local
// fInvRot is the rotation matrix to go from global to local
fRot=fRot.RotateAxes(fXaxis, fYaxis, fZaxis);
fInvRot=fRot.Inverse();
}
//______________________________________________________________________________
Bool_t OpticalAdaptor::IsAbsorbed(Photon *p) const {
//
// Applies the filter to the photon
//
Double_t value = 0;
switch ( fFilterType ) {
case kSquare:
// really stupid filter
return ((p->wl < 300*nm || p->wl > 400*nm) &&
EsafRandom::Get()->Rndm() > fRadiomEff);
case kBG3:
if ( p->wl < fgBG3->GetXmin() || p->wl > fgBG3->GetXmax() )
return kTRUE;
// value = Exp(fFilterThickness*Log(fgBG3->GetValue(p->wl))/(1.*mm));
value = Power(fgBG3->GetValue(p->wl),
fFilterThickness/(1.*mm));
return ( EsafRandom::Get()->Rndm() > value );
case kMultilayer:
if ( p->wl < fgMultiLayer->GetXmin() || p->wl > fgMultiLayer->GetXmax() )
return kTRUE;
// value = Exp(fFilterThickness*Log(fgMultiLayer->GetValue(p->wl))/(1.*mm));
value = Power(fgMultiLayer->GetValue(p->wl),
fFilterThickness/(1.*mm));
return ( EsafRandom::Get()->Rndm() > value );
default:
FatalError("Unknown filter type.");
}
return kTRUE;
}
//______________________________________________________________________________
Double_t OpticalAdaptor::IsHit( const Photon& ph) const {
//
// checks if ph is going to hit the top face of the OA and, in that case,
// it returns the distance it has to travel. Otherwise returns -1
//
if (fZaxis.Dot(ph.dir) > -kTolerance) {
cerr << "OpticalAdaptor::IsHit: photon is not going toward this oa" << endl;
return -1;
}
EVector topface_pos = fPos+fZaxis.Unit()*GetThickness();
EVector intPoint = ph.pos-ph.dir*((ph.pos-topface_pos).Dot(fZaxis)/ph.dir.Dot(fZaxis));
if ( isInside(goLocal(intPoint-fPos), TOP) )
return (intPoint-ph.pos).Mag();
return -1;
}
//______________________________________________________________________________
int OpticalAdaptor::whichFace(Photon *p) const {
// lpos is the position of photon in the local reference
EVector lpos=p->pos - fPos;
lpos=goLocal(lpos);
if(fabs(lpos[Z]-fHeight)<kTolerance) {
// photon on the top surface
#ifdef DEBUG
cout<<"TOP"<<endl;
#endif /* DEBUG */
return TOP;
} else if(fabs(lpos[X]-fSide)<kTolerance) {
// photon on the right surface
#ifdef DEBUG
cout<<"RIGHT"<<endl;
#endif /* DEBUG */
return RIGHT;
} else if(fabs(lpos[X])<kTolerance) {
// photon on the left surface
#ifdef DEBUG
cout<<"LEFT"<<endl;
#endif /* DEBUG */
return LEFT;
} else if(fabs(lpos[Y]+fSide)<kTolerance) {
// lpos[Y] < 0
// photon on the front surface
#ifdef DEBUG
cout<<"FRONT"<<endl;
#endif /* DEBUG */
return FRONT;
} else if(fabs(lpos[Y])<kTolerance) {
// photon on the back surface
#ifdef DEBUG
cout<<"BACK"<<endl;
#endif /* DEBUG */
return BACK;
} else if(fabs(lpos[Z])<kTolerance) {
// photon on the bottom surface
#ifdef DEBUG
cout<<"BOTTOM"<<endl;
#endif /* DEBUG */
return BOTTOM;
} else throw runtime_error("whichFace: not of surface");
}
//______________________________________________________________________________
vector<EVector> OpticalAdaptor::intPoints(Photon *p) const {
EVector lpos, ldir;
lpos=p->pos - fPos;
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);
// TOP
if(face==TOP || ldir[Z] < 0) ips[TOP]=out;
else {
dist=fHeight-lpos[Z];
ips[TOP]=lpos + ldir*(dist/ldir[Z]);
#ifdef DEBUG
cout<<"dist, ips[TOP]: "<<dist<<", "<<ips[TOP]<<endl;
#endif /* DEBUG */
}
// BOTTOM
if(face==BOTTOM || ldir[Z] > 0) ips[BOTTOM]=out;
else {
dist=-lpos[Z];
ips[BOTTOM]=lpos + ldir*(dist/ldir[Z]);
#ifdef DEBUG
cout<<"dist, ips[BOTTOM]: "<<dist<<", "<<ips[BOTTOM]<<endl;
#endif /* DEBUG */
}
// FRONT
if(face==FRONT || ldir[Y] > 0) ips[FRONT]=out;
else {
dist=-lpos[Y]-fSide;
ips[FRONT]=lpos + ldir*(dist/ldir[Y]);
#ifdef DEBUG
cout<<"dist, ips[FRONT]: "<<dist<<", "<<ips[FRONT]<<endl;
#endif /* DEBUG */
}
// BACK
if(face==BACK || ldir[Y] < 0) ips[BACK]=out;
else {
dist=-lpos[Y];
ips[BACK]=lpos + ldir*(dist/ldir[Y]);
#ifdef DEBUG
cout<<"dist, ips[BACK]: "<<dist<<", "<<ips[BACK]<<endl;
#endif /* DEBUG */
}
// RIGHT
if(face==RIGHT || ldir[X] < 0) ips[RIGHT]=out;
else {
dist=fSide-lpos[X];
ips[RIGHT]=lpos + ldir*(dist/ldir[X]);
#ifdef DEBUG
cout<<"dist, ips[RIGHT]: "<<dist<<", "<<ips[RIGHT]<<endl;
#endif /* DEBUG */
}
// LEFT
if(face==LEFT || ldir[X] > 0) ips[LEFT]=out;
else {
dist=0-lpos[X];
ips[LEFT]=lpos + ldir*(dist/ldir[X]);
#ifdef DEBUG
cout<<"dist, ips[LEFT]: "<<dist<<", "<<ips[LEFT]<<endl;
#endif /* DEBUG */
}
return ips;
}
//______________________________________________________________________________
bool OpticalAdaptor::isInside(const EVector &v, int face) const {
switch(face) {
case TOP:
if( v[X] > kTolerance && v[X] < fSide-kTolerance &&
-v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
fabs(v[Z]-fHeight) < kTolerance ) return true;
break;
case BOTTOM:
if( v[X] > kTolerance && v[X] < fSide-kTolerance &&
-v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
fabs(v[Z]) < kTolerance ) return true;
break;
case RIGHT:
if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
-v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
fabs(v[X]-fSide) < kTolerance ) return true;
break;
case LEFT:
if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
-v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
fabs(v[X]) < kTolerance ) return true;
break;
case FRONT:
if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
v[X] > kTolerance && v[X] < fSide-kTolerance &&
fabs(v[Y]+fSide) < kTolerance ) return true;
break;
case BACK:
if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
v[X] > kTolerance && v[X] < fSide-kTolerance &&
fabs(-v[Y]) < kTolerance ) return true;
break;
}
return false;
}