// ESAF : Euso Simulation and Analysis Framework
// class TestOpticalAdaptor
// $Id: OpticalSystem.cc,v 1.20 2005/04/15 18:54:42 thea Exp $
//
#include <math.h>
#include "TRotation.h"
#include "OpticalSystem.hh"
#include "EsafRandom.hh"
#include "EEvent.hh"
ClassImp(OpticalSystem)
ClassImp(TestOpticalSystem)
//______________________________________________________________________________
TestOpticalSystem::TestOpticalSystem() : OpticalSystem() {
fPos=EVector(0,0,Conf()->GetNum("TestOpticalSystem.pos.Z"))*mm;
fDZup=Conf()->GetNum("TestOpticalSystem.DZup")*mm;
fDZdown=Conf()->GetNum("TestOpticalSystem.DZdown")*mm;
fFocalDistance=Conf()->GetNum("TestOpticalSystem.FocalDistance")*mm;
fR=Conf()->GetNum("TestOpticalSystem.Radius")*mm;
psf=new Interpolate("config/Optics/"+Conf()->GetStr("TestOpticalSystem.psf.filename"));
tr=new Interpolate("config/Optics/"+Conf()->GetStr("TestOpticalSystem.tr.filename"));
}
//______________________________________________________________________________
Photon *TestOpticalSystem::Transport(Photon *p) const {
// theta is the angle between the photon direction and the lens axis
double theta=p->dir.Angle(EVector(0.,0.,1.));
if(theta > TMath::PiOver2()) theta=TMath::Pi() - theta;
#ifdef DEBUG
cout<<"position = "<<p->pos<<endl;
cout<<"direction = "<<p->dir<<endl;
cout<<"theta: "<<theta/TMath::Pi()*180.<<endl;
#endif /* DEBUG */
// check if photon is in FOV (FOV=30 deg)
if( theta > 30*deg) {
// photon outside FOV
// destroy it
#ifdef DEBUG
cout<<"photon outside FOV: destroyed"<<endl;
#endif /* DEBUG */
p=0;
return p;
}
// photon inside FOV
// check if it is absorbed or refracted
TRandom* rndm = EsafRandom::Get();
if(rndm->Rndm() > tr->GetValue(theta/deg)) {
#ifdef DEBUG
cout<<"photon absorbed by opticalSystem: destroyed"<<endl;
#endif /* DEBUG */
p=0;
return p;
}
EVector dir=p->dir.Unit();
// propagate photon to the center of the lens
double DZ=(dir[Z]>0 ? fDZdown : -fDZup);
EVector dummy=dir*(DZ/dir[Z]);
p->time+=dummy.Mag()/c_light;
p->pos+=dummy;
if(p->pos.Perp() > fR) {
// hit wall
#ifdef DEBUG
cout<<"Hit Wall: "<<p->pos<<endl;
#endif /* DEBUG */
p=0;
return p;
}
// calculate new direction
// This is for a plane focalplane
// EVector newDir=dir*(fFocalDistance/fabs(dir[Z])) - (p->pos - fPos);
// This is for a spherical focalplane with radius fFocalDistance
// newDir is the point on the focal surface
double thetain=p->dir.Theta(); // save original theta
// used later for PSF
EVector newDir=dir*fFocalDistance - (p->pos - fPos);
p->dir=newDir.Unit();
// Point Spread Function
EVector PSF=getPSF(p, thetain);
newDir+=PSF;
p->dir=newDir.Unit();
// propagate photon to the end of the lens
dir=p->dir;
DZ=(dir[Z]>0 ? fDZup : -fDZdown);
dummy=dir*(DZ/dir[Z]);
p->time+=dummy.Mag()/c_light;
p->pos+=dummy;
p->pos[Z]=(dir[Z]>0 ? Top() : Bottom());
if(p->pos.Perp() > fR) {
// hit wall
#ifdef DEBUG
cout<<"Hit Wall: "<<p->pos<<endl;
#endif /* DEBUG */
p=0;
return p;
}
#ifdef DEBUG
cout<<"new position = "<<p->pos<<endl;
cout<<"new direction = "<<p->dir<<endl;
#endif /* DEBUG */
return p;
}
//______________________________________________________________________________
EVector TestOpticalSystem::getPSF(Photon *p, double thetain) const {
// theta is the angle between the photon direction and the lens axis
double theta=p->dir.Angle(EVector(0.,0.,1.));
if(theta > TMath::PiOver2()) theta=TMath::Pi() - theta;
double psf_t=psf->GetValue(thetain/deg);
#ifdef DEBUG
cout<<"psf("<<thetain/deg<<"): "<<psf_t<<endl;
#endif /* DEBUG */
TRandom* rndm=EsafRandom::Get();
double r=rndm->Gaus(0., psf_t);
double phi=rndm->Rndm()*360.*deg;
EVector psf_v(1,0,0);
psf_v.SetPhi(phi);
psf_v.SetMag(r);
#ifdef DEBUG
cout<<"psf_v: "<<psf_v<<endl;
cout<<"p->dir: "<<p->dir<<endl;
#endif /* DEBUG */
TRotation rot;
EVector axis=p->dir.Cross(EVector(0.,0.,1.));
rot=rot.Rotate(-asin(axis.Mag()), axis.Unit());
psf_v.Transform(rot);
#ifdef DEBUG
cout<<"psf_v = "<<psf_v<<endl;
cout<<"dot = "<<psf_v.Dot(p->dir)<<endl;
#endif /* DEBUG */
return psf_v;
}
//______________________________________________________________________________
Bool_t OpticalSystem::HitSide(Photon *p) const{
// checks if p is going to hit the dead border of the optics when it's goint toward
// the optics i.e. if it's going to hit the walls between lenses from outside
Double_t dt = CylinderIntersection( p );
if (dt < 0 || (p->pos+p->dir*dt)[Z]-fDZup < kTolerance
|| (p->pos+p->dir*dt)[Z]-fDZdown < kTolerance )
return kFALSE;
return kTRUE;
}